第6章 线性系统matlab特性分析与参数计算

控制相关专业研究生选修课程系统建模方法马宏军东北大学信息学院控制理论与导航技术研究所2013年3月逸夫楼203MATLAB实现控制系统稳定性分析•稳定是控制系统的重要性能,也是系统能够工作的首要条件,因此,如何分析系统的稳定性并找出保证系统稳定的措施,便成为自动控制理论的一个基本任务.•线性系统的稳定性取决于系统本身的结构和参数,而与输入无关.1直接判定法•根据稳定的充分必要条件判别线性系统的稳定性,最简单的方法是求出系统所有极点,并观察是否含有实部大于0的极点,如果有,系统则不稳定.•然而实际的控制系统大部分都是高阶系统,这样就面临求解高次方程,求根工作量很大•在Matlab中只需分别调用函数roots(den)或eig(A)即可,这样就可以由得出的极点位置直接判定系统的稳定性.1直接判定法已知控制系统的传递函数为若判定该系统的稳定性,输入如下程序:G=tf([1,7,24,24],[1,10,35,50,24]);roots(G.den{1})运行结果:ans=-4.0000-3.0000-2.0000-1.0000由此可以判定该系统是稳定系统.245035102424723423ssssssssG2根轨迹法判断系统的稳定性•根轨迹法是一种求解闭环特征方程根的简便图解法,它是根据系统的开环传递函数极点、零点的分布和一些简单的规则,研究开环系统某一参数从零到无穷大时闭环系统极点在s平面的轨迹.•控制工具箱中提供了rlocus函数,来绘制系统的根轨迹,利用rlocfind函数,在图形窗口显示十字光标,可以求得特殊点对应的K值.2根轨迹法判断系统的稳定性已知一控制系统,H(s)=1,其开环传递函数为:绘制系统的轨迹图.程序为:G=tf(1,[1320]);rlocus(G);[k,p]=rlocfind(G)21sssKsG2根轨迹法判断系统的稳定性根轨迹图如图1所示,2根轨迹法判断系统的稳定性光标选定虚轴临界点,程序结果为:selected_point=0-0.0124ik=0.0248p=-2.0122-0.9751-0.0127光标选定分离点,程序结果为:selected_point=-1.9905-0.0124ik=0.0308p=-2.0151-0.9692-0.01582根轨迹法判断系统的稳定性上述数据显示了增益及对应的闭环极点位置.由此可得出如下结论:(1)0k0.4时,闭环系统具有不同的实数极点,表明系统处于过阻尼状态;(2)k=0.4时,对应为分离点,系统处于临界阻尼状态;(3)0.4k6时,系统主导极点为共轭复数极,系统为欠阻尼状态;(4)k=6时,系统有一对虚根,系统处于临界稳定状态;(5)k6时,系统的一对复根的实部为正,系统处于不稳定状态.3用Nyquist曲线判断系统的稳定性Matlab提供了函数Nyquist来绘制系统的Nyquist曲线,若方法2系统分别取k=4和k=10(图2为阶跃响应曲线),通过Nyquist曲线判断系统的稳定性,程序如下:num1=[4];num2=[10];den1=[1,3,2,0];gs1=tf(num1,den1);gs2=tf(num2,den1);hs=1;gsys1=feedback(gs1,hs);gsys2=feedback(gs2,hs);t=[0:0.1:25];figure(1);subplot(2,2,1);step(gsys1,t)subplot(2,2,3);step(gsys2,t)subplot(2,2,2);nyquist(gs1)subplot(2,2,4);nyquist(gs2)3用Nyquist曲线判断系统的稳定性图2阶跃响应曲线3用Nyquist曲线判断系统的稳定性•奈氏稳定判据的内容是:若开环传递函数在s平半平面上有P个极点,则当系统角频率X由-∞变到+∞时,如果开环频率特性的轨迹在复平面上时针围绕(-1,j0)点转P圈,则闭环系统稳定,否则,是不稳定的.3用Nyquist曲线判断系统的稳定性图3Nyquist曲线3用Nyquist曲线判断系统的稳定性当k=4时,从图3中k=4可以看出,Nyquist曲不包围(-1,j0)点,同时开环系统所有极点都位于平面左半平面,因此,根据奈氏判据判定以此构成闭环系统是稳定的,这一点也可以从图2中k=4系统单位阶跃响应得到证实,从图2中k=4可以看出系统约23s后就渐渐趋于稳定.当k=10时,从图3中k=10可以看出,Nyquist曲线按逆时针包围(-1,j0)点2圈,但此时P=0,所以据奈氏判据判定以此构成的闭环系统是不稳定的,图2中k=10的系统阶跃响应曲线也证实了这一点,系统振荡不定。4bode图法判断系统的稳定性•bode判据,实质上是Nyquist判据的引伸.本开环系统是最小相位系统,即P=0,用Xc表示对数幅频特性曲线与横轴(0dB)交点的频率,Xg表示对数相频特性曲线与横轴(-180o)交点的频率,则对数判据可表述如下：•在P=0时,若开环对数幅频特性比其对数相频特性先交于横轴,即XcXg,则闭环系统稳定;若开环对数幅频特性比其对数相频特性后交于横轴,即XcXg,则闭环系统不稳定;若Xc=Xg,则闭环系统临界稳定.4bode图法判断系统的稳定性num1=[4];num2=[10];den1=[1,3,2,0];gs1=tf(num1,den1);gs2=tf(num2,den1);hs=1;gsys1=feedback(gs1,hs);gsys2=feedback(gs2,hs);t=[0:0.1:25];figure(1);subplot(1,1,1);bode(gs1)4bode图法判断系统的稳定性图4k=4时开环系统的bode图由图4开环系统的bode图可知,XcXg,故当k=4时,此时的闭环系统必然稳定.实际上,系统的控制bode图还可用于系统相对稳定性的分析.时域指标的求取（单位）阶跃函数（Stepfunction）1()0()00ttrtt（单位）斜坡函数（Rampfunction）0()00ttrtt（单位）加速度函数（Accelerationfun）210()200ttrtt（单位）脉冲函数（Impulsefunction）0()00ttt正弦函数（Simusoidalfunction）()sinrtAt典型的试验信号时域响应性能指标动态性能指标0tMp超调量允许误差10.90.50.1trtpts表示性能指标td,tr,tp,Mp和ts的单位阶跃响应曲线tdc(t)0.02或0.05)(c)(c)(c)(c上升时间响应曲线从稳态值的10%上升到90%，所需的时间。上升时间越短，响应速度越快dt:rt峰值时间：响应曲线达到峰值所需要的时间。pt时域响应性能指标动态性能指标0tMp超调量允许误差10.90.50.1trtpts表示性能指标td,tr,tp,Mp和ts的单位阶跃响应曲线tdc(t)0.02或0.05)(c)(c)(c)(c调节时间响应曲线达到并永远保持在一个允许误差范围内，所需的最短时间。用稳态值的百分数（通常取5%或2%）表示。st%100)()()(%hhthp超调量%：指响应的最大偏离量h(tp)于终值之差的百分比，即()()()ctcc即：ttr时域响应性能指标求取1.峰值时间（timetopeak）峰值时间可由下述命令获得：[Y,k]=max(y)%求出y的峰值及相应的时间timetopeak=t(k)%获得峰值时间时域响应性能指标求取2.超调量（Percentovershoot）超调量可由以下命令获得：C=dcgain(G)%求取系统的终值[Y,k]=max(y)%求出y的峰值及相应的时间percentovershoot=100*(Y-C)/C%计算超调量时域响应性能指标求取C=dcgain(G)%求取系统的终值n=1whiley(n)0.1*C%通过循环，求取输出第一n=n+1次到达终值的10%的时间endm=1;whiley(n)0.9*C%通过循环，求取输出第一m=m+1次到达终值的90%的时间endrisetime=t(m)-t(n)3.上升时间时域响应性能指标求取4.调节时间（settingtime）C=dcgain(G)i=length(t)while(y(i)0.98*C)&(y(i)1.02*C)i=i-1endSettingtime=t(i)时域响应性能指标求取例已知二阶系统传递函数为3(),(13)(13)Gssisi编程求取系统的性能指标。时域响应性能指标求取G=zpk([],[-1+3*i,-1-3*i],3);%计算最大峰值时间和它对应的超调量C=dcgain(G)[y,t]=step(G);plot(t,y)grid[Y,k]=max(y);timeopeak=t(k)%取得最大峰值时间percentovershoot=100*(Y-C)/C%计算超调量%计算上升时间n=1whiley(n)Cn=n+1endrisetime=t(n)时域响应性能指标求取%计算稳态响应时间i=length(t);while(y(i)0.98*C)&(y(i)1.02*C)i=i-1;endsettingtime=t(i)时域响应性能指标求取C=0.3000timeopeak=1.0928percentovershoot=34.7385risetime=0.6954settingtime=3.4771运行程序后，结果为：时域响应性能指标求取频域指标：增益裕度和相角裕度•margin()[Gm,Pm,Wcg,Wcp]=margin(num,den)Gm=增益裕度（dB）Pm=相角裕度（deg）G(s)=num/den[Gm,Pm,Wcg,Wcp]=margin(sys)•Gm_dB=20*log10(Gm)[Gm,Pm,Wcg,Wcp]=margin(mag,phase,w)margin(num,den)直接在Bode图上标明margin(sys)Wcg=相角-180对应频率Wcp=增益0dB对应频率bode(num,den)margin(num,den)Gm=26.9405Pm=48.5831Wcg=47.5445Wcp=2.939bode(num,den)[Gm,Pm,Wcg,Wcp]=margin(num,den)title(['Gm=',num2str(Gm),'Pm=',num2str(Pm)])