3 MATLAB控制系统仿真(讲的很多)

2.5 应用MATLAB控制系统仿真

MATLAB是一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算和图形显示于一体，构成了一个方便的界面友好的用户环境。由控制领域专家推出的MATLAB工具箱之一的控制系统(Control System)，在控制系统计算机辅助分析与设计方面获得了广泛的应用，并且MATLAB工具箱的内容还在不断增加，应用范围也越来越宽。

控制系统的分析与设计方法，不论是古典的还是现代的，都是以数学模型为基础进行的。MATLAB可以用于以传递函数形式描述的控制系统。在本节中，我们首先以一个典型的动力学系统“弹簧—重物—阻尼器”的数学模型为例，说明如何使用MATLAB进行辅助分析。之后，我们讨论传递函数和结构图。特别的，主要介绍以下内容：如何使用MATLAB求解多项式，计算传递函数的零点和极点，计算闭环传递函数，计算结构图的等效变换以及闭环系统对单位阶跃输入的响应等。

这部分所包含的MATLAB函数有：roots，roots1，series，parallel，feedback，cloop，poly，conv，polyval，printsys，pzmap和step等。

2.5.1. 举例

一个弹簧—重物—阻尼器动力学系统如图2-1所示。重物M的位移由y(t)表示，用微分方程描述如下：

该系统在初始位移作用下的瞬态响应为

其中cos-1 ，初始位移是y(0)。系统的瞬态响应当＜1时为欠阻尼，当＞１时为过阻尼，当＝1时为临界阻尼。考虑

过阻尼情况：y(0)=0.15 m n= (弧度/秒) ( )

欠阻尼情况：y(0)=0.15 m n= (弧度/秒) ( )

利用MATLAB程序—unforced.m，可以显示初始位移为y(0)的物体自由运动曲线，如图2-63所示。在unforced.m程序中，变量y(0)，n，t，1和2的值由指令直接输入工作区，然后运行unforced.m程序就可以产生响应曲线。MATLAB提供了一种交互式的应用环境，我们可以在指令窗口随时修改n，1和2值并再次运行unforced.m程序。在上述欠阻尼和过阻尼情况下的响应曲线如图2-所示。可以看到程序运行后阻尼系数的值被标注在不同的曲线上，以防止混淆。 对于弹簧—重物—阻尼器系统，利用MATLAB求微分方程的解是容易而有效的。一般说来，当分析一个具有多样性的输入和初始条件以及元件参数的闭环反馈控制系统时，欲直接得到这些因素与系统响应的关系是比较困难的。在这种情况下，可以利用MATLAB非常有效的数值计算能力，反复运行MATLAB程序并绘出曲线图，就可以得出结论。

>>y0=0.15;wn=sqrt(2);

 >>zeta1=3/(2*sqrt(2）; zeta2=1/(2*sqrt(2）;  >>t=[0:0.1:10];  >>unforced （a）MATLAB指令窗口  % unforced.m

 %计算系统在给定初始条件下的自由运动  t1=acos(zeta1)*ones(1,length(t）;  t2=acos(zeta2)*ones(1,length(t）;

 c1=(y0/sqrt(1-zeta1^2）; c2=(y0/sqrt(1-zeta2^2）;  y1=c1*exp(-zeta1*wn*t)sin(wn*sqrt(1-zeta1^2)*t+t1);  y2=c2*exp(-zeta2*wn*t)sin(wn*sqrt(1-zeta2^2)*t+t2);  %计算运动曲线的包络线

 bu=c2*exp(-zeta2*wn*t);bl=-bu;  %画图

 plot(t,y1, ‘-’,t,y2,‘-’,t,bu, ‘--’,bl, ‘--’),grid  xlabel(‘Time[sec]’), ylabel(‘y(t) Displacement[m]’)  text(0.2,0.85,[‘oeverdamped zeta1=’,num2str(zeta1),])  text(0.2,0.80,[‘underdamped zeta2=’,num2str(zeta2),])  （b）分析弹簧—重物—阻尼器的MATLAB程序unforced.m

图2- 弹簧—重物—阻尼器的自由运动曲线

MATLAB可以用来分析以传递函数形式描述的系统。由于传递函数是多项式的比值，所以如何处理多项式是MATLAB首先需要解决的问题。记住，分子多项式和分母多项式都必须在MATLAB指令中指定。

在MATLAB中多项式由行向量组成，而这些行向量包含了降次排列的多项式系数。例如多项式p(s)=1s3+3s2+0s1+4s0，按图2-65的格式输入p=[1 3 0 4]，请注意，尽管s的系数为0，它也一定要包含在确定p(s)的行向量中。 如果p是一个包含降幂排列的p(s)系数的行向量，那么roots(p)是一个包含多项式根的列向量。相反的，如果r是一个包含多项式根的列向量，那么，poly(r)是一个包含降幂排列多项式系数的行向量，见图2-65。我们可以用上述的roots()函数计算多项式p(s)的根，但是当多项式有重根时，函数roots1()能给出更精确的结果。

 >>p=[1 3 0 4];  >>r=roots(p)  r= -3.3553e+00

 1.7765e-01+1.0773e+00j  1.7765e-01-1.0773e+00j  >>p=poly(r)

 p=1.000 3.000 0.000-0.000j 4.000+0.000j  图2-65 输入多项式并求根

矩阵乘法由MATLAB的conv()函数完成。假设我们要把两个多项式相乘合并成一个多项式n(s)，即

n(s)= (3s2 +2s +1) (s +4)

= 3s+14s+9s +4

与此运算相关的MATLAB函数就是conv()，如图2-66所示。函数polyval()用来计算多项式的值。多项式n(s)在s=-5处的值为n(-5)= -66，见图2-66。  >>p=[3 2 1];q=[1 4];  >>n=conv(p,q)  n=3 14 9 4

 >>value=polyval(n,-5)  value=-66

图2-66 MATLAB的conv()函数和polyval()函数

3

2

2.5.2. 传递函数

下面我们将介绍如何利用MATLAB函数获得传递函数在复平面内的零极点分布图。 设传递函数为G(s)=num/den，其中num和den均为多项式。利用函数

[P , Z]=pzmap(num , den)

可以获得G(s)的零极点位置，即P为极点位置列向量，Z为零点位置列向量。该指令执行后自动生成零极点分布图。 考虑传递函数

和

利用一系列MATLAB指令和函数，我们可以计算传递函数的零极点、特征方程和两个传递函数相除

，还可以在复平面上获得G(s)/H(s)的零极点分布图。

传递函数G(s)/H(s)的零极点图如图2-67所示，相应的MATLAB指令如图2-68所示。图2-67中清楚的表示出了5个零点的分布，但是看上去只有两个极点，实际上这是不可能的，因为我们知道极点的个数一定大于或等于零点的个数。应用函数roots1()可以看到，实际上有四个极点位于s= -1。因此，在同一位置上重复的极点或零点在零极点分布图上是区分不出的。

图2-67 零极点图

2.5.3. 结构图模型

假设我们已为某系统的各部分建立了传递函数，下一步的任务就是利用MATLAB函数将这些部分联接起来构成一个闭环控制系统，实现结构图的转换，计算从输入R(s)到输出C(s)的传递函数。为了方便介绍MATLAB函数，我们从结构图的基本变换开始。一个简单的开环控制系统可以通过G1(s)与G2(s)两个环节的串联而得到，利用series()函数可以求串联连接的传递函数，函数的具体形式为

[num,den]= series(num1,den1, num2,den2)

例如G1 (s)和G2 (s)的传递函数分别为

，

则

串联函数的用法示于图2-69。

>>numg=[6 0 1];deng=[1 3 3 1];

>>z=roots(numg) z=

0+0.4082j 0-0.4082j

>> p=roots1(deng) p= -1 -1 -1

>>n1=[1 1]; n2=[1 2]; d1=[1 2*j]; d2=[1 –2*j]; d3=[1 3]; >>numh=conv(n1,n2); denh=conv(d1,conv(d2,d3）; >>num=conv(numg,denh); den=conv(deng,numh); >>printsys(num,den) num/den=

6s^5+18s^4+25s^3+75s^2+4s+12 ａ s^5+6s^4+14s^3+16s^2+9s+2

>>pzmap(num,den)

>>title(‘Pole-Zero Map’)

图2-68 绘制零极点图指令

>>num1=[1];den1=[500 0 0]; >>num2=[1 1];den2=[1 2];

>>[num,den]=series(num1,den1,num2,den2); >>printsys(num,den) num/den= s+1

500s^3+1000s^2

图2-69 series函数的用法

当系统是以并联的形式连接时，利用parallel()函数可得到系统的传递函数。指令的具体形式为

[num,den]= parallel (num1,den1, num2,den2)

如果系统以反馈方式构成闭环，则系统的闭环传递函数为



求闭环传递函数的MATLAB函数有两个：cloop()和feedback()，其中cloop()函数只能用于H (s)=1（即单位反馈）的情况。 cloop()函数的具体用法为

[num,den]= cloop (numg,deng, sign) 其中numg和deng分别为G (s)的分子和分母多项式，sign=1为正反馈，sign= -1为负反馈（默认值）。



图2-70 闭环反馈系统的结构图

feedback()函数的具体用法为

[num,den]= feedback (numg,deng,numh,denh, sign)

其中numh为H (s)的分子多项式，denh为分母多项式。 假设闭环反馈系统的结构图如图2-70所示，被控对象G(s)和控制部分Gc (s)以及测量环节H (s)的传递函数分别为

， ，

应用series()函数和feedback()函数求解闭环传递函数的MATLAB指令如图2-71所示。

>>numg=[1];deng=[5 0 0]; >>numc=[1 1];denc=[1 2]; >>numh=[1];denh=[1 10];

>>[num1,den1]=series(numc,denc,numg,deng); >>[num,den]=feedback(num1,den1,numh,denh,-1); >>printsys(num,den) num/den=

s^2+11s+10

5s^4+60s^3+100s^2+s+1

图2-71 feedback()函数的应用

例2.12 一个多环的反馈系统如图2-49所示，给定各环节的传递函数为

， ， ，

， ，

试求闭环传递函数GB(s)=C(s)/R(s)。 解 求解过程可按如下步骤进行： 步骤1：输入系统各环节的传递函数 步骤2：将H2的综合点移至G2后 步骤3：消去G3，G2，H2环 步骤4：消去包含H3的环

步骤5：消去其余的环，计算GB (s) 根据上述步骤的MATLAB指令以及计算结果在图2-72中。

>>ng1=[1];dg1=[1 10]; >>ng2=[1];dg2=[1 1]; >>ng3=[1 0 1];dg3=[1 4 4]; >>ng4=[1 1];dg4=[1 6];

>>nh1=[1];dh1=[1]; >>nh2=[2];dh2=[1]; >>nh3=[1 1];dh3=[1 2];

>>[n1,d1]=series(ng2,dg2,nh2,dh2); >>[n2,d2]=feedback(ng3,dg3,n1,d1,-1); >>[n3,d3]=series(n2,d2,ng4,dg4); >>[n4,d4]=feedback(n3,d3,nh3,dh3,-1); >>[n5,d5]=series(ng1,dg1,ng2,dg2); >>[n6,d6]=series(n5,d5,n4,d4); >>[n7,d7]=cloop(n6,d6,-1); >>printsys(n7,d7) num/den=

s^4+ 3s^3+ 3s^2+3s+2 2s^6+38s^5+261s^4+1001s^3+1730s^2+16s+732 图2-72 多环结构图简化

有时我们需要关心闭环传递函数是否有零极点对消的情况出现。当然通过pzmap()或roots()函数可以查看传递函数是否有相同的零极点，另外还可以使用minreal()函数除去传递函数共同的零极点因子。正如图2-73所示，如果传递函数有相同的零极点，应用minreal()函数后，传递函数的分子和分母多项式各减少了一阶，消去了相同的零极点。 

>>numg=[1 6 11 6];deng=[1 7 12 11 5]; >>printsys(numg,deng) numg/deng=

s^3+6s^2+11s+6 s^4+7s^3+12s^2+11s+5

>>[num,den]=minreal(numg,deng); >>printsys(num,den) 1 pole-zeros cancelled num/den=

s^2+4s+3

s^3+6s^2+6s+5

图2-73 minreal()函数的应用

最后重新考虑例2.2所示的位置随动系统，我们的目的是计算闭环系统在输入作用下的响应。在给定各元件参数并忽略La和令ML = 0的情况下，其结构图如图2-74所示。计算的第一步是求闭环传递函数GB (s)= c(s) / r(s)，求解过程及结果如图2-75所示。第二步是利用step()函数计算参考输入 r (t)为单位阶跃信号时输出 c (t)的响应。可见，特征方程是2阶的，且n =52，=0.012，由于阻尼比很小，可以预料响应会强烈震荡。







图2-74 位置随动系统的结构图

 





  图2-76 位置随动系统的阶跃响应曲线



>>num1=[200]; den1=[20];num2=[1]; den2=[2 0.5 0]; >>num3=[0.2 0]; den3=[1]; num4=[0]; den4=[1]; >>[na,da]=series(num1,den1,num2,den2); >>[nb,db]=feedback(na,da,num3,den3,-1); >>[nc,dc]=series(nb,db,num4,den4); >>[num,den]=cloop(nc,dc,-1); >>printsys(num,den) num/den= 00 2s^2+2.5s+00 >>t=[0:0.005:3];

>>[y,t]=step(num,den,t); >>plot(t,y),grid

图2-75 位置随动系统的结构图简化及阶跃响应指令

图2-76给出了位置随动系统的阶跃响应曲线。y(t)，即 c (t)的离散时间点t将以行向量的形式给出，它从0秒开始，按0.005秒的步长增加，直到3秒为止。使用plot()函数用于画出y(t)曲线，grid函数用于给图形加上网格。 由于控制系统的性能指标通常以阶跃响应的形式给出，因此step()函数是非常重要的。有关step()函数的内容将在后续章节中进一步介绍。

3.4 应用MATLAB分析控制系统的性能

这一节将用两个例子描述反馈控制的优点，同时说明如何利用MATLAB来分析控制系统。系统分析的主要内容包括如何抑制干扰、如何减小稳态误差、如何调节瞬态响应以及如何减少系统对参数变化的影响等。

第一个例子是带有负载转矩干扰信号的电枢控制直流电动机。开环系统结构图如图3-37(a)所示，为了改善系统性能，加入速度反馈如图3-37(b)所示。系统的各元器件参数值在表3-6中给出。

表3-6 速度控制系统的参数

参数名 参数值 Ra 1 Km 10 J 2 B 0.5 Ke 0.1 Ka Ks 1

图3-37 速度控制系统结构图

从图中可以看出，系统有Ua(s)(或Vr(s)和ML(s)两个输入。由于这是一个线性系统，按叠加定理可以分别考虑两个输入的作用结果。为了研究干扰对系统的作用，可令

Ua(s)=0(或Vr(s)=0)，此时只有干扰ML(s)起作用。相反地，为了研究参考输入对系统的响应，可令ML(s)=0。如果系统具有很好的抗干扰能力，则干扰信号ML(s)对输出(s)的影响就应该很小，下面就来验证此结论。

首先，考虑图3-37(a)所示的开环系统，从ML(s)到o(s)(此处的下标“o”表示开环)的传递函数为

假设干扰信号为单位阶跃信号，即ML(s)=1/s。利用MATLAB可以计算系统的单位阶跃响应如图3-38(a)所示，而用于分析此开环控制系统的MATLAB程序文本opentach.m示于图3-38(b)。

在输入信号Ua(s)=0的情况下，稳态误差就是干扰响应 o(t)的终值。在图3-38(a)的曲线中，干扰响应 o(t)在t = 7秒后已近似不变，所以近似稳态误差值为

 o(∞) ≈ -0.663(弧度/秒)

同样，通过计算从ML (s)到 c(s)(此处下标“c”表示闭环)的闭环传递函数可分析图3-37(b)所示闭环系统的抗干扰性能。对于干扰输入的闭环传递函数为

(a) 开环速度系统对阶跃干扰的响应曲线

(b) MATLAB程序文本:opentach.m

图3-38 开环速度控制系统分析

闭环系统对单位阶跃干扰输入的响应曲线 (t)和MATLAB程序文本

closedtach.m分别示于图3-39(a)(b)。同前，稳态误差就是 (t)的终值，稳态误差的近似值为

在本例中，闭环系统与开环系统对单位阶跃干扰信号的输出响应的稳态值之比为

可见通过引入负反馈已明显减小了干扰对输出的影响，这说明闭环反馈系统具有抑制噪声特性。

(a) 闭环系统对阶跃干扰的响应曲线

(b) MATLAB程序文本:opentach.m 图3-39 闭环速度控制系统分析

第二个例子是分析闭环控制系统的控制器增益K对瞬态响应的影响。图3-40是闭环控制系统的结构图。在参考输入R (s)和干扰输入N (s)同时作用下系统的输出为

图3-40 反馈控制系统的结构图

如果单纯考虑增益K对参考输入产生的瞬态响应的影响，可以预计增加K将导致超调量增加、调整时间减少和响应速度提高。在增益K=20和K=100时，系统对参考输入的单位阶跃响应曲线以及相应的MATLAB程序文本gain_kr.m示于图3-41。对比两条响应曲线，可以看出上述预计的正确性。尽管在图中不能明显看出增大K能减少调整时间，但是这一点可以通过观察MATLAB程序的运行数据得以验证。这个例子说明了控制器增益K是如何改变系统瞬态响应的。根据以上分析，选择K=20可能是一个比较好的方案。尽管如此，在做出最后决定之前还应该考虑其他因素。

(a) 阶跃响应曲线

% K=20和K=100时，参考输入的单位阶跃响应：gain_kr.m numg=[1];deng=[1 1 0]; K1=100;K2=20;

num1=[11 K1];num2=[11 K2];den=[0 1]; %简化结构图

[na,da]=series(num1,den,numg,deng); [nb,db]=series(num2,den,numg,deng); [numa,dena]=cloop(na,da); [numb,denb]=cloop(nb,db); %选择时间间隔 t=[0:0.01:2.0];

[c1,x,t]=step(numa,dena,t); [c2,x,t]=step(numb,denb,t); plot(t,c1,'--',t,c2)

xlabel('Time[sec]'),ylabel('Cr(t)'),grid

(b) MATLAB程序文本:gain_kr.m

图3-41 单位阶跃输入的响应分析

在对K做出最后选择之前，非常重要的是要研究系统对单位阶跃干扰的响应，有关结果和相应的MATLAB程序文本如图3-42所示。从中可以看到，增加K减少了单位干扰响应的幅值。对于K=20和K=100，响应的稳态值分别为0.05和0.01。对干扰输入的稳态值可按终值定理求得

如果仅从抗干扰的角度考虑，选择K=100更合适。

在本例中所求出的稳态误差、超调量和调整时间(2%误差)归纳于表3-7。

在控制系统设计中有很多成熟的经验，设计者常常要权衡利弊，综合考虑。在这个例子中，增加K导致了更好的抗干扰性，然而减少K可以使系统具有更好的瞬态响应性能。如何选择K的最终权力留给了设计者。尽管MATLAB软件对控制系统的分析和设计很有帮助，但是控制工程师的经验往往更重要。

(a) 阶跃响应曲线

% K=20和K=100时，干扰输入的单位阶跃响应：gain_kn.m numg=[1];deng=[1 1 0]; K1=100;K2=20;

num1=[11 K1];num2=[11 K2];den=[0 1]; %简化结构图

[numa,dena]=feedback(numg,deng,num1,den); [numb,denb]=feedback(numg,deng,num2,den); %选择时间间隔 t=[0:0.01:2.5];

[c1,x,t]=step(numa,dena,t); [c2,x,t]=step(numb,denb,t); plot(t,c1,'--',t,c2)

xlabel('Time[sec]'),ylabel('Cn(t)'),grid

(b) MATLAB程序文本:gain_kn.m

图3-42 单位阶跃干扰的响应分析

表3-7 K=20和K=100时，控制系统的响应特性 K值 超调量p% 调整时间ts 稳态误差ess K=20 4% 1.0秒 5%

K=100 22% 0.7秒 1%

(a) 被控对象变化时系统的灵敏度

%系统灵敏度分析

K=20;num=[1 1 0];den=[1 12 K]; %取s=j,的范围为10～10，共取200点

-1

3

w=logspace(-1,3,200);s=w*j; %S为灵敏度，S2为灵敏度的近似值 n=s.^2+s;d=s.^2+12*s+K;S=n./d; n2=s;d2=K;S2=n2./d2;

subplot(211),plot(real(S),imag(S) xlabel('Real(S)'),ylabel('imag(S)'),grid subplot(212),loglog(w,abs(S),w,abs(S2) xlabel('w[rad/sec]'),ylabel('abs(S)')

(b) MATLAB程序文本

图3-43 系统的灵敏度分析

最后来分析被控对象变化时系统的灵敏度。在本例中，被控对象的传递函数和闭环系统的传递函数分别为

，

系统的灵敏度可由式(3-78)得出

利用上式可计算不同s值所对应的灵敏度SG ，并绘制出频率—灵敏度曲线。图3-43(a)中给

-1 -4

出的是K=20，s=j，=10~10时，系统的灵敏度相对于频率的变化曲线，图3-43(b)是相应的MATLAB程序文本。在低频段，系统的灵敏度可近似为

可见，增大K值，可以减少系统的灵敏度。

另外，在MATLAB函数中还有impulse函数(脉冲响应函数)和lsim函数(任意输入下的响应函数)，这两个函数的用法与step函数相近，这里不再介绍

4．6 利用Matlab绘制系统的根轨迹

本章前面的内容介绍了控制系统根轨迹的绘制以及利用系统大致的根轨迹图分析系统性能的方法，若要由根轨迹获得系统在某一特定参数下准确的性能指标或者准确的闭环极点，需要依据幅值条件精确地作图。如果利用MATLAB工具箱中函数，则可方便、准确地作出根轨迹图，并利用图对系统进行分析。 MATLAB工具箱中，求系统根轨迹的几个常用函数有rlocus, rlocfind, sgrid，下面通过具体的例子来说明这些函数的应用。

例4-13 控制系统的开环传递函数为

绘制系统的根轨迹图。

解 利用函数rlocus函数可直接作出系统的根轨迹图，程序如下： % example4-13 %

num=[1,5];

dun=[1,6,11,6,0];

rlocus(num,dun)

执行该程序后，可得到如图4-20所示的根轨迹。

图4-20 例4-13题根轨迹图

利用函数rolcus可画出系统的根轨迹图后，可用rlocfind函数在根轨迹上选择任意极点，得到相应的开环增益 和其它闭环极点。 例4-14 控制系统的开环传递函数为

绘制系统的根轨迹图，并确定根轨迹的分离点及相应的开环增益

解 将开环传递函数写为

。

Matlab程序如下： % example4-14 %

num=[1];

den=[0.0002,0.03,1,0]; rlocus(num,den)

title(‘Root Locus’) [k,p]=rlocfind(num,den)

程序执行过程中，先绘出系统的根轨迹，并在图形窗口中出现十字光标，提示用户在根轨迹上选择一点，这时，将十字光标移到所选择的地方，可得到该处对应的系统开环增益及其它闭环极点。此例中，将十字光标移至根轨迹的分离点处，可得到 k = 9.6115 p =

-107.7277 -21.9341 -20.3383

若光标能准确定位在分离点处，则应有两个重极点，即 根轨迹图如图4-21所示。

相等。程序执行后，得到的

图4-21 例4-14 系统的根轨迹

例4-15 开环系统的传递函数为

绘制系统的根轨迹，并分析系统的稳定性。

解 Matlab程序如下 %example 4-15 %

num=[1,3]; den1=[1,6,5];

den=conv(den1,den1); figure(1)

rlocus(num,den)

[k,p]=rlocfind(num,den) % analyzing the stability figure(2) k=159;

num1=k*[1,3]; den=[1,6,5];

den1=conv(den,den);

[num,den]=cloop(num1,den1,-1); impulse(num,den)

title(‘Impulse Response (k=160)’) % analyzing the stability figure(3) k=161

num1=k*[1,3]; den=[1,6,5];

den1=conv(den,den);

[num,den]=cloop(num1,den1,-1); impulse(num,den)

由第1段程序得到根轨迹后，将十字线移到根轨迹与虚轴的交点上，可得到在交点处

，可知，使系统临界稳定的根轨迹增益为

当系统的根轨迹增益

100％，已接近临界稳定的状态。当

，根轨迹如图4-15(a)所示。

时，系统是稳定的，但系统的阻尼非常小，超调量近似为

时，系统具有正实部的复数极点，系统不稳

定。执行第2、3段程序后，得到图4-22(b)和(c)。由图4-22(b)和(c)可清楚看到，当 时，闭环系统的脉冲响应是收敛的，故系统稳定，而当 是发散的，故系统不稳定。

时，闭环系统的脉冲响应

图4-22(a) 例题4-15系统的根轨迹图

图4-22(b)

时系统的脉冲响应 图4-22(c)

时系统的脉冲

响应

例4-16 单位反馈系统的开环传递函数为

试绘制系统的根轨迹，确定当系统的阻尼比

解 Matlab程序如下： %example 4-16 %

num=[4 3 1];

时系统的闭环极点，并分析系统的性能。

den=[3 5 1 0]; sgrid

rlocus(num,den)

[k,p]=rlocfind(num,den)

执行以上程序后，可得到绘有由等阻尼比系数和自然频率构成的栅格线的根轨迹图，如图4-23所示。屏幕出现选择根轨迹上任意点的十字线，将十字线的交点移至根轨迹与 的等阻尼比线相交处，可得到

k =

0.2752

p =

-1.70

-0.1623 + 0.1653i -0.1623 - 0.1653i

此时系统有三个闭环极点，一个负实数极点，两个共轭复数极点，实数极点远离虚轴，其距虚轴的距离是复数极点的10倍，且复数极点附近无闭环零点，因此，这对共轭复数极点满足主导极点的条件，系统可简化为由主导极点决定的二阶系统，系统的性能可用二阶系统的分析方法得到。 系统的特征方程为

所以，系统的闭环传递函数为

图4-23 例4-16根轨迹图

5.6 MATLAB在系统频域分析中的应用

前面介绍了几种用曲线表示开环系统频率特性的方法，利用这些曲线可以分析系统的稳定性及其它的频域性能指标，还可由开环频率特性求取控制系统的闭环频率特性。利用MATLAB工具箱中函数，可以准确地作出系统的频率特性曲线，为控制系统的设计和分析提供了极大的方便。

MATLAB工具箱中，绘制系统频率特性曲线的几个常用函数有bode, nyquist，nichols, ngrid, margin等，下面通过具体的例子来说明这些函数的应用。

1． Bode函数

函数bode可以求出系统的Bode 图，其格式为 [mag,phase,w]=bode(num,den) [mag,phase,w]=bode(num,den,w)

bode(num,den)可以绘制传递函数为 时系统的Bode图。 bode(num,den,w)可利用指定的频率值 绘制系统的Bode图。

当带输出变量引用函数时，可以得到系统Bode图相应的幅值、相角及频率值。其中 mag=

,

phase=

2．Nyquis函数

函数nyquis的功能是求系统的奈氏曲线，格式为 [re,im,w]=nyquist(num,den) [re,im,w]=nyquist(num,den,w)

nyquist(num,den)可以得到开环传递函数为 时系统的nyquist曲线。 nyquist(num,den,w)可以根据指定的频率值 绘制系统的nyquist曲线。

当带输出变量引用函数时，可以得到系统nyquist曲线的数据，而不直接绘出nyquist曲线。

3．nichols函数

nichols函数可求出系统频率特性的nichols图线，格式为 [mag,phase,w]= nichols (num,den) [mag,phase,w]= nichols(num,den,w)

nichols (num,den)可以得到开环传递函数为 时系统的nichols 图线。 nichols (num,den,w)可以根据指定的频率值 绘制系统的nichols 图线。

当带输出变量引用函数时，可以得到系统nichols 图线的数据，而不直接绘出nichols 图线。

4．ngrid函数

ngrid 函数的功能是绘制nichols图线上的网格。格式为 ngrid

ngrid(‘new’)

ngrid函数可在nichols图线上加网格线，ngrid(‘new’) 可在绘制网格前清除原图，然后再设置程hold,on，这样，后续的nichols函数可与网格绘制在一起。

5． margin函数

margin可求出开环系统的幅值裕度和相角裕度，其格式为 margin(num,den)

[gm,pm,wcg,wcp]=margin(mag,phase,w)

margin(num,den)可计算系统的相角裕度和幅值裕度，并绘制出Bode图。

margin(mag,phase,w)可以由幅值裕度和相角裕度绘制出Bode图，其中，mag、 phase和 w是由bode得到的幅值裕度、相角裕度和频率。

当带输出变量引用函数时，仅计算幅值裕度、相角裕度及幅值穿越频率wcg和相角穿越频率wcp，不绘制Bode图。

例5－13 二阶振荡环节的传递函数为

绘制

取不同值时的Bode图。

解 取n=8，取[0.1;0.1;1.0]，由bode函数得到Bode图，MATLAB程序如下： %example 5-13 % wn=5;

kosi=[0.1:0.2:1.0]; w=logspace(-1,1,100); figure(1)

num=[wn.^2]; for kos=kosi

den=[1 2*kos*wn wn.^2];

[mag,pha,w1]=bode(num,den,w); subplot(2,1,1);hold on semilogx(w1,mag); subplot(2,1,2);hold on semilogx(w1,pha); end

subplot(2,1,1);grid on title('Bode Plot');

xlabel('Frequency(rad/sec)'); ylabel('Gain dB'); subplot(2,1,2);grid on

xlabel('Frequency(rad/sec)'); ylabel('Phase deg'); hold off

程序执行后，得到二阶振荡环节的Bode图如图5－63所示。

例5－14 开环系统的传递函数为

1． 绘制系统的奈氏曲线，并用奈氏判据判断系统的稳定性， 2． 求闭环系统的单位脉冲响应。

解 利用nyquist函数绘制奈氏曲线如图5－所示。由图可知，奈氏曲线不围绕(-1,j0)点，N=0,开环传递函数有一个右半s平面的极点，P=1，由奈氏判据

系统不稳定。又由系统的脉冲响应曲线图5－65，可知脉冲响应是发散的，也说明系统不稳定。 MATLAB程序如下： %example 5-14 % k=30; z=[1];

p=[-2 -3 2];

[num,den]=zp2tf(z,p,k); figure(1)

nyquist(num,den) title('Nyquist Plot'); figure(2)

[num1,den1]=cloop(num,den); impulse(num,den)

title('Impulse Response')

例5－15 系统的开环传递函数为

绘制Nichols图。 解 MATLAB程序为 %example 5-15 %

k=16.7/.0125; z=[0];

p=[-1.25 -4 -16];

[num,den]=zp2tf(z,p,k); figure(1) ngrid('new')

nichols(num,den) title('Nichols Plot')

程序执行后，可得到系统的Nichols图如图5－66所示。

6．5 MATLAB 在系统校正中的应用

MATLAB为系统校正提供了方便的工具，改变校正装置的参数，可清楚地看到校正对系统性能的影响。下面例题给出了MATLAB语句的应用，读者还可根据问题的需要，自己编写合适的函数，使程序更加简洁。

例6－9 单位反馈系统的开环传递函数为

试确定串联校正装置的特性，使系统满足在斜坡函数作用下系统的稳态误差小于0.1，相角裕度

。

解 根据系统稳态精度的要求，选择开环增益K=12，求原系统的相角裕度。 >> num=12; den=[2,1,0];

[gm,pm,wcg,wcp]=margin(num,den); [gm,pm,wcg,wcp] margin(num,den)

ans =

Inf 11.68 Inf 2.4240

可知，原系统相角裕度为 ，不满足指标要求，系统的Bode如图6－34所示。考虑采用串联超前校正装置，以增加系统的相角裕度。

选择超前校正装置的最大超前相角

，则有

为使超前装置的相角补偿作用最大，选择校正后系统的剪切频率在最大超前相角发生的频率上，由

图6－34，当幅值为时-10lglg4.77=-6.8dB时，相应的频率为3.6rad/s。选择此频率作为校正后系统的剪切频率

由式（6－12）确定参数

初选校正装置的传递函数为

校正后系统的开环传递函数为

校正后系统的Bode图如图6－35所示。 由图可知，校正后系统的相角裕度为 即为所求。

例6－10 设控制系统的开环传递函数为

，满足性能指标的要求，故该超前装置

试设计一串联校正装置，使校正后系统的相角裕度不小于40，幅值裕度不低于10dB,剪切频率大于(1rab/s)。

解 （1）作校正前系统的对数频率特性。 >> num=8;

den=[0.125,0.75,1,0]; margin(num,den) [gm,pm,wcg,wcp] ans =

0.7500 -7.5156 2.8284 3.2518

程序运行后，得到图6－36。

o

由图6－36可知，原系统具相角裕度和幅值裕度均为负值，故系统不稳定。考虑到系统的剪切频率为 ，大于系统性能指标要求的剪切频率，故采用滞后装置对系统进行校正。根据相角裕度 角的影响，选择校正后系统的相角裕度为

的要求和滞后装置对系统相

由图6－36知，对应相角为

时的频率为c'=1.1rad/s>1，幅值为15.7dB。

取20lgb=-15.7dB，得b=0.1。取滞后装置的第二个转折频率为0.1c'=0.11，有1/bT=0.11，

则T=55.43。初选校正装置的传递函数为

作出校正后系统的Bode图如图6－37中所示。由图，可得到校正后系统的相角裕度为

，幅值裕度为12.73dB，剪切频率为c=1.1rad/s，满足系统性能指标的要求，故

初选校正装置合适，校正后系统的开环传递函数为