- 1 -MATLAB 实现控制系统稳定性分析稳定是控制系统的重要性能,也是系统能够工作的首要条件,因此,如何分析系统的稳定 性并找出保证系统稳定的措施,便成为自动控制理论的一个基本任务.线性系统的稳定性取 决于系统本身的结构和参数,而与输入无关.线性系统稳定的条件是其特征根均具有负实部. 在实际工程系统中,为避开对特征方程的直接求解,就只好讨论特征根的分布,即看其是 否全部具有负实部,并以此来判别系统的稳定性,由此形成了一系列稳定性判据,其中最重要 的一个判据就是 Routh 判据.Routh 判据给出线性系统稳定的充要条件是:系统特征方程式 不缺项,且所有系数均为正,劳斯阵列中第一列所有元素均为正号,构造 Routh 表比用求根判 断稳定性的方法简单许多,而且这些方法都已经过了数学上的证明,是完全有理论根据的,是 实用性非常好的方法. 但是,随着计算机功能的进一步完善和 Matlab 语言的出现,一般在工程实际当中已经不 再采用这些方法了.本文就采用 Matlab 对控制系统进行稳定性分析作一探讨. 1 系统稳定性分析的 Matlab 实现 1.1 直接判定法 根据稳定的充分必要条件判别线性系统的稳定性,最简单的方法是求出系统所有极点,并 观察是否含有实部大于 0 的极点,如果有,系统则不稳定.然而实际的控制系统大部分都是高 阶系统,这样就面临求解高次方程,求根工作量很大,但在 Matlab 中只需分别调用函数 roots(den)或 eig(A)即可,这样就可以由得出的极点位置直接判定系统的稳定性. 已知控制系统的传递函数为（1） 245035102424723423ssssssssG若判定该系统的稳定性,输入如下程序: G=tf(1,7,24,24,1,10,35,50,24); roots(G.den1) 运行结果: ans =-4.0000-3.0000-2.0000-1.0000 由此可以判定该系统是稳定系统. 1.2 用根轨迹法判断系统的稳定性 根轨迹法是一种求解闭环特征方程根的简便图解法,它是根据系统的开环传递函数极点、 零点的分布和一些简单的规则,研究开环系统某一参数从零到无穷大时闭环系统极点在 s 平 面的轨迹.控制工具箱中提供了 rlocus 函数,来绘制系统的根轨迹,利用 rlocfind 函数,在图 形窗口显示十字光标,可以求得特殊点对应的 K 值. 已知一控制系统,H(s)=1,其开环传递函数为:(2) 21sssKsG绘制系统的轨迹图. 程序为: G=tf(1,1 3 2 0);rlocus(G); k,p=rlocfind(G) 根轨迹图如图 1 所示,光标选定虚轴临界点,程序 结果为: - 2 -图 1 系统的根轨迹图 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.0158 上述数据显示了增益及对应的闭环极点位置.由此可得出如下结论:- 3 -(1)06 时,系统的一对复根的实部为正,系统处于不稳定状态. 1.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) 奈氏稳定判据的内容是:若开环传递函数在 s 平半平面上有 P 个极点,则当系统角频率 X 由-变到+时,如果开环频率特性的轨迹在复平面上时针围绕(-1,j0)点转 P 圈,则闭环系统 稳定,否则,是不稳定的.图 2 阶跃响应曲线 当 k=4 时,从图 3 中 k=4 可以看出,Nyquist 曲不包围(-1,j0)点,同时开环系统所有极点 都位于平面左半平面,因此,根据奈氏判据判定以此构成闭环系统是稳定的,这一点也可以从 图 2 中 k=4 系统单位阶跃响应得到证实,从图 2 中 k=4 可以看出系统约 23 s 后就渐渐趋于 稳定.当 k=10 时,从图 3 中 k=10 可以看- 4 -图 3 Nyquist 曲线 出,Nyquist 曲线按逆时针包围(-1,j0)点 2 圈,但此时 P=0,所以据奈氏判据判定以此构成的闭 环系统是不稳定的,图 2 中 k=10 的系统阶跃响应曲线也证实了这一点,系统振荡不定。 1.4 bode 图法判断系统的稳定性 bode 判据,实质上是 Nyquist 判据的引伸.本开环系统是最小相位系统,即 P=0,用 Xc 表 示对数幅频特性曲线与横轴(0dB)交点的频率,Xg 表示对数相频特性曲线与横轴(-180o)交点 的频率,则对数判据可表述如下： 在 P=0 时,若开环对数幅频特性比其对数相频特性先交于横轴,即 XcXg,则闭环系统不稳定;若 Xc=Xg,则闭环系统临界稳定. 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)- 5 -图 4 k=4 时开环系统的 bode 图 由图 4 开环系统的 bode 图可知,Xc p=1 1 2 2 3 5; roots(p)二）用根轨迹法判别系统稳定性：对给定的系统的开环传递函数，进行仿真。1某系统的开环传递函数为 ，在 command window 窗口输入程序，记录系统闭环零极点图及零极点数据，判断该闭环系统是否稳定。 clear n1=0.25 1; d1=0.5 1 0; s1=tf(n1,d1); sys=feedback(s1,1); P=sys.den1;p=roots(P) pzmap(sys) p,z=pzmap(sys)- 6 -2某系统的开环传递函数为 ，在 command window 窗口输入程序，记录系统开环根轨迹图、系统开环增益及极点，确定系统稳定时 K 的取值范围。 clear n=1;d=conv(1 1 0,0.5 1); sys=tf(n,d); rlocus(sys) k,poles=rlocfind(sys)三）频率法判别系统稳定性：对给定的系统的开环传递函数，进行仿真。1已知系统开环传递函数 ，在 command window 窗口输入程序，用 Bode 图法判稳，记录运行结果，并用阶跃相应曲线验证（记录相应曲线）1）绘制开环系统 Bode 图，记录数据。 num=75*0 0 0.2 1; den=conv(1 0,1 16 100); sys=tf(num,den); Gm,Pm,Wcg,Wcp=margin(sys) margin(sys)2）绘制系统阶跃响应曲线，证明系统的稳定性。 num=75*0 0 0.2 1; den=conv(1 0,1 16 100); s=tf(num,den); sys=feedback(s,1); t=0:0.01:30; step(sys,t)2已知系统开环传递函数 ，在 command window 窗口输入程序，用 Nyquist 图法判稳，记录运行结果，并用阶跃相应曲线验证（记录相应曲线）。1）绘制 Nyquist 图，判断系统稳定性。 clear num=10000; den=1 5 100 0; GH=tf(num,den); nyquist(GH)2）用阶跃响应曲线验证系统的稳定性 num=10000;den=1 5 100 0; s=tf(num,den);- 7 - sys=feedback(s,1); t=0:0.01:0.6; step(sys,t)tftf 用法：用法： G=tf(2G=tf(2 1,11,1 2 2 2)2) 或或 s=tf(s)s=tf(s) ; ; % %定义定义 s s 为传递函数拉普拉斯算子；为传递函数拉普拉斯算子； G=G=（2s+12s+1）/(s2+2s+2);/(s2+2s+2); % %定义传递函数；定义传递函数； 其实生成的传递函数可以任意计算。其实生成的传递函数可以任意计算。 set(G)set(G)可以得到传递函数对象的属性，可以修改或预设其属性，如可以得到传递函数对象的属性，可以修改或预设其属性，如 G=tf(2G=tf(2 1,11,1 2 2 2,variable,p);%2,variable,p);%修改使用的变量修改使用的变量 G=tf(2G=tf(2 1,11,1 2 2 2,inputdelay,0.25);%2,inputdelay,0.25);%设置输入延迟，即设置输入延迟，即 G=exp(-0.25s)G=exp(-0.25s)（2s+12s+1） /(s2+2s+2)/(s2+2s+2) G=tf(1G=tf(1 3 3 2,12,1 5 5 7 7 3,0.1);%3,0.1);%设置离散情况的采样周期设置离散情况的采样周期tfdatatfdata： 功能描述：获得功能描述：获得 tftf 模型传递函数的参数模型传递函数的参数 对于对于 SISOSISO 系统：系统： G=tf(2G=tf(2 1,11,1 2 2 2);2); num,den=tfdata(G,v);num,den=tfdata(G,v); 对于离散系统：对于离散系统： num,den,Ts=tfdata(G)num,den,Ts=tfdata(G)； 其实要得到系统的参数，可以直接引用传递函数的属性，如其实要得到系统的参数，可以直接引用传递函数的属性，如 G.denG.den 等等函数名 zpk 功能描述：生成零极点增益传递函数模型或转换成零极点模型 G=zpk(-1,-3,0,-2,-5,10); 可以用于转化 G=tf(-10 20 0,1 7 20 28 19 5) sys=zpk(G);zpkdata 功能描述：获取零极点增益模型的参数 z ,p ,k=zpkdata(G,v);filt()filt() 功能描述：生成 DSP 形式的离散传递函数： 例子：生成采样时间为 0.5 的 DSP 形式传递函数：H=filt(2 1,1 0.4 2,0.5)%求闭环传函- 8 -num1,den1=series(1,1 1,1 0,1 0 2);num2,den2=feedback(1,1 0