Matlab求解阶跃响应性能指标
概述
阶跃响应性能指标定义
Matlab函数
示例
求一阶系统阶跃响应性能指标
求复杂系统阶跃响应性能指标
求解不同阻尼比时二阶系统的阶跃响应性能指标
联系作者

概述
工科的同学对应阶跃响应应该不会陌生，简单来说，阶跃响应是指输入量发生阶跃变化时动态系统的输出阶跃响应，通过分析输出阶跃响应的性能指标，可以分析和比较动态系统的动态性能和稳态性能。
阶跃响应性能指标定义
阶跃响应性能指标主要包括稳态值、上升时间、峰值时间和超调量，定义如下：

稳态值ys：当时间趋近于无穷大时，阶跃响应的输出值，ys=y(∞)。

上升时间tr：输出阶跃响应达到90%稳态值时所对应的时刻。

峰值时间tm：输出阶跃响应峰值ym所对应的时刻。

超调量σ：输出阶跃响应峰值ym与稳态值ys之差所占稳态值ys的百分比, σ%= (ym-ys)/ys。

调整时间ts：输出阶跃响应进入稳态值ys±Δ误差带范围内所对应的时刻, 一般取Δ=0.02或Δ=0.05。


Matlab函数
按照阶跃响应性能指标的定义，作者使用Matlab开发了函数Fun_Step_Performance.m，使用数值算法求出各类阶跃响应的性能指标值，函数简单、易用、通用性好。
function [ys,tr,ts,tm,ov] = Fun_Step_Performance(t,y,drawflag)
% [ys,tr,ts,ov] = Fun_Step_Performance(t,y) 标准阶跃响应的性能指标求解
% 本程序适用于标准阶跃响应曲线，末尾时间必须已经接近稳态值
% t-y 为阶跃响应的时间-输出配对序列，可由[y,t] = step(sys)求得
% drawflag为时候作图标志，不输入或输入非0值时，默认作图，输入0时不做图
% ys 稳态值
% tr 上升时间，默认为0-90%的上升时间
% ts 调整时间，默认为2%的调整时间
% tm 为峰值时间
% ov 超调量 %
% e.g.
%  sys = tf(1,[1 2*0.5*1 1]);
%  [y,t] = step(sys,15);
%  [ys,tr,ts,tm,ov] = Fun_Step_Performance(t,y,1);

示例
求一阶系统阶跃响应性能指标
% Eg 1 一阶系统
sys = tf(1,[3 1]);
[y,t] = step(sys,25);
[ys,tr,ts,tm,ov] = Fun_Step_Performance(t,y);


%% 阶跃响应指标结果：
上升时间：7s
调整时间：11.5s
峰值时间：25s，超调量：0%
稳态值：1
%% 阶跃响应指标结果显示结束

求复杂系统阶跃响应性能指标
sys = tf(1,[1 2*0.20*1 1]) * tf(1,[2 1]) * tf([1.5 1],[1 2*0.25*3 9]);  % 5阶系统
[y,t] = step(sys,35);
[ys,tr,ts,tm,ov] = Fun_Step_Performance(t,y);


%% 阶跃响应指标结果：
上升时间：2.0877s
调整时间：17.3158s
峰值时间：3.5614s，超调量：40.1285%
稳态值：0.111
%% 阶跃响应指标结果显示结束

求解不同阻尼比时二阶系统的阶跃响应性能指标
% Eg 3 求解不同阻尼比时二阶系统的阶跃响应特性
wn = 1;    % 固有频率
kes_vet = [0.2 0.4 0.6 0.8 1 1.2];   % 阻尼比序列
figure
hold on
for ii = 1:length(kes_vet)
    kes = kes_vet(ii);
    sys = tf(1,[1 2*kes*wn wn^2]);   % 二阶系统传递函数
    [y,t] = step(sys,50);            % 阶跃响应
    [ys(ii),tr(ii),ts(ii),tm(ii),ov(ii)] = Fun_Step_Performance(t,y,0);  % 求解阶跃响应,不绘图
    plot(t,y)
    Str{ii} = [ '\xi = '  num2str(kes)];
end
legend(Str)
xlabel('时间t/s')
ylabel('输出响应y')
 
figure
subplot(221)
plot( kes_vet,tr,'-o' )
xlabel('阻尼比\xi')
ylabel('上升时间/s')
subplot(222)
plot( kes_vet,ts,'-o' )
xlabel('阻尼比\xi')
ylabel('调整时间/s')
subplot(223)
plot( kes_vet,tm,'-o' )
xlabel('阻尼比\xi')
ylabel('峰值时间/s')
subplot(224)
plot( kes_vet,ov,'-o' )
xlabel('阻尼比\xi')
ylabel('超调量/%')




tr =
    1.8349    2.0737    2.4540    2.9954    4.0000    4.9407

ts =
   19.2661    8.2949    5.8282    3.6866    5.5000    7.9051

tm =
    3.2110    3.4562    3.9877    5.1843   40.0000   50.0000

ov =
   52.6622   25.3725    9.4610    1.5144         0         0

对于二阶系统，阻尼比的变化不影响输出稳态值，随着阻尼比增加，上升时间逐步增大、调整时间先减小再增大、峰值时间逐步变大、超调量逐步变小；当阻尼比在0.707左右时，上升时间和调整时间较快，且超调量很小，系统综合性能较好，工程上通常设计阻尼比在0.707左右，称之为最佳阻尼比。
联系作者
有Matlab/Simulink方面的技术问题，欢迎发送邮件至944077462@qq.com讨论。

更多Matlab/Simulink原创资料，欢迎关注微信公众号：Matlab Fans
源程序下载:

1. csdn资源: Matlab求解阶跃响应指标

2. 扫码关注微信公众号Matlab Fans，回复BK06获取百度网盘下载链接。

以下是第二个传递函数的做法，若想求第一个或其他的传递函数的参数，只需要改变分子矩阵num和分母矩阵即可
num1=[0.01];
den1=[1,0.002,0.01];
t=0:0.001:10;%为了确保精度，即各个参数的小数的位数，减小步长
y=step(num1,den1,t);
n=length(t);
[ymax,ind]=max(y);%y是原系统的阶跃响应，是一个二维矩阵，返回本身曲线的最大值和对应的时间矩阵变量的第几位
yss=y(n);%我们这里就设我们时间的最终值对应的就是稳态值
mp=(ymax-yss)/yss;%由超调量定义可以计算超调量
for k=1:n%该层循环是计算上升时间即第一次达到稳态值的时间
if y(k)<=yss&y(k+1)>=yss%第一次达到稳态值
tr=t(k+1);%用tr存储稳态值
break
end
end
for i=n:-1:1%该层循环用于计算调整时间，我们是从最终时间（无穷远的时间）往前开始遍历查找，否则从前开始查找的话会出错
    if y(i)>=1.02*yss|y(i)<=0.98*yss%当得到的值与稳态值差的绝对值与稳态值之比为2%时说明到达了调整时间
        ts=t(i);%ts存储调整时间
        break
    end
end
disp('峰值时间为:')  
disp(t(ind))
disp('上升时间为:')
disp(tr)
disp('调整时间为:')
disp(ts)
disp('超调量为:')
disp(mp)

1.Matlab作单位阶跃响应曲线的三维图



定义闭环系统传递函数如下：

H(s)=1/(s^2+2*zeta*s+1)
%Matlab作单位阶跃响应曲线的三维图
%标准二阶系统响应曲线 wn=1
t=0:0.2:10;
zeta=[0 0.2 0.4 0.8 1];
for n=1:6
num=1;
den=[1 2*zeta(n) 1];
[y(1:51,n),x,t]=step(num,den,t) ;
end


plot(t,y)


mesh(t,zeta,y')




2.用Matlab求上升时间、峰值时间、最大过调量和调整时间


定义闭环系统传递函数如下：

H(s)=25/(s^2+6*s+25)


%求上升时间、峰值时间、最大过调量和调整时间
%可以应用到高阶系统
num=[0 0 25];
den=[1 6 25];
t=0:0.005:5;
[y,x,t]=step(num,den,t) ;
r=1;
while y(r)<1.0001;r=r+1;end;
rise_time=(r-1)*0.005


[ymax,tp]=max(y);
peak_time=(tp-1)*0.005


max_overshoot=ymax-1


s=1001;while y(s)>0.98&y(s)<1.02;s=s-1;end;
settle_time=(s-1)*0.005







运行后结果：

rise_time =    0.5550

peak_time =   0.7850

max_overshoot =    0.0948

settle_time =    1.1850



求单位脉冲响应的两种方法（求斜坡响应也可用此方法）



在初始条件为零时，H(s)的单位脉冲响应也就是sH(s)的单位阶跃响应

H(s)=1/（s^2+0.2*s+1）



num=[1];

den=[1 0.2 1]；

impulse(num,den)



或 



num=[0 1 0];

den=[1 0.2 1]；

step(num,den)

1. 准备
终值：c(∞)

上升时间 tr：响应从峰值的10%上升到峰值的90%所需要的时间；而阶跃响应则是从终值的10%上升到终值的90%所需要的时间；对有振荡的系统，也可以定义为从0到第一次到达终值所需的时间。

峰值时间 tp：响应到达第一个峰值所需的时间。

调节时间 ts：响应到达并保持在终值的正负5%误差带内所需的最短时间；有时也用终值的正负2%。

超调量 σ%：峰值c(tp )超出终值的百分比，即：


题目：
已知二阶系统传递函数为



试分别用游动鼠标法和编程法求取系统的峰值，超调量，上升时间以及调节时间。
2. 编程法：
s = tf('s');   % 用字母 s 构造传递函数
Gs = 3/(s^2+2*s+10); % 构造 Gs
t = 0:0.005:5; % 设置横轴范围和步长
[y,x,t] = step(Gs,t); % 根据步长逐步响应传递函数
[ymax,tp] = max(y); % 获取最大值的点的数据

r=1;
r1=0;
while y(r) < (ymax*0.9) % 过滤小于90%的
    if y(r) > (ymax*0.1) % 过滤大于10%的
        r1=r1+1;
    end;
    r=r+1;
end;

rise_time = r1*0.005 % 上升时间
peak_time = (tp-1)*0.005 % 峰值时间
ystable = dcgain(Gs); % 稳态值
max_overshoot = (ymax-ystable)/ystable % 超调量

r2=1001; % 由 (5-0)/0.005+1 求得
while y(r2) > ystable*0.98 && y(r2) < ystable*1.02 
    % 从稳态倒推回去，求得首次达到并维持在稳态值正负2%的范围的值
    r2 = r2 - 1;
end
settle_time = (r2-1)*0.005 % 调节时间

求得：
rise_time =   0.6100
peak_time = 1.0450
max_overshoot = 0.3509
settle_time = 3.5350

3. 游动鼠标法：
游动鼠标法不适合用于plot()命令画出的图形，只能在非plot()函数输出的曲线上求取。

1.获取必要数据
s = tf('s');
Gs = 3/(s^2+2*s+10); % 构造 Gs
t = 0:0.005:5;
step(Gs,t)  % 画出图像
[y,x,t] = step(Gs,t); 
[ymax,tp] = max(y); 
ymax	% 峰值
ym10 = ymax * 0.1  % 峰值 10%
ym90 = ymax * 0.9  % 峰值 90%
ystable = dcgain(Gs) % 稳态值
ys98 = ystable * 0.98  % 稳态 -2%
ys102 = ystable * 1.02   %  稳态 +2%

得到图像并求得：
ymax = 0.4053
ym10 = 0.0405
ym90 = 0.3647
ystable = 0.3000
ys98 = 0.2940
ys102 = 0.3060

2.点击工具栏中的数据游标，在图中分别标出纵坐标为 ymax = 0.4053，ym10 = 0.0405，ym90 = 0.3647，ys98 = 0.2940，ys102 = 0.3060 的五点：





3.收集图像中 Time 对应的 Amplitude 数据，并计算可得：

tr = 0.783 - 0.176 = 0.607

tp = 1.04

ts = 3.53    (因为0.295位于此波谷最低处，而0.295>0.294，故不取4.19)

σ% = (0.405 - 0.3)/0.3 × 100% = 35%
下面是另外两个博客，和这个关系比较密切，关于 PID 算法解释的：

PID算法模型分析：基于温度控制

临界比例度法整定P、PI、PID控制器的参数的matlab算法实现
>>> End <<<