控制系统数字仿真与cad》第7章控制系统的计算机辅助分析.ppt

《控制系统数字仿真与cad》第7章控制系统的计算机辅助分析.ppt》由会员分享，可在线阅读，更多相关《控制系统数字仿真与cad》第7章控制系统的计算机辅助分析.ppt（127页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第5章 控制系统的计算机辅助分析l系统仿真实质上就是对描述系统的数学模型进行求解。对控制系统来说，系统的数学模型实际上就是某种微分方程或差分方程，因而在仿真过程中需要根据某种数值算法从系统给定的初始值出发，逐步地计算出每一个时刻系统的响应，最后绘制出系统的响应曲线，由此来分析系统的性能。在前面曾经介绍过一般常微分方程的数值解法，该方法是系统仿真的基础。其实，对于各种线性系统模型在典型输入信号作用下来说，当然没有必要采用那些通用的算法来完成这种任务，而是应该充分地利用线性系统的特点，采取更简单的方法来得到问题的解。这样做不但会大大提高运算的效率，而且可以提高仿真的精度和可靠性。本章主要介绍利用M

2、ATLAB的控制系统工具箱所提供的函数对线性系统进行计算机分析和处理。7.1 控制系统的稳定性分析在分析控制系统时，首先遇到的问题就是系统的稳定性。对线性系统来说，如果一个连续系统的所有极点都位于左半s平面，则该系统是稳定的。对离散系统来说，如果一个系统的全部极点都位于单位圆内，则此系统可以被认为是稳定的。由此可见，线性系统的稳定性完全取决于系统的极点在根平面上的位置。本节主要介绍几种利用MATLAB来判断系统稳定性的方法。l1.利用极点判断系统的稳定性l 判断一个线性系统稳定性的一种最有效的方法是直接求出系统所有的极点，然后根据极点的分布情况来确定系统的稳定性，对于极点的求取我们在上节中已作

3、过介绍，下面举例说明其判断方法。l例5-1 已知闭环系统的传递函数为l判断系统的稳定性，并给出不稳定极点。解：可以利用下面的MATLAB程序%ex5_1.mnum=3 2 1 4 2;den=3 5 1 2 2 1;z,p=tf2zp(num,den);ii=find(real(p)0);n1=length(ii);if(n10)disp(The Unstable Poles are:);disp(p(ii);else disp(System id Stable);end执行结果显示：The Unstable Poles are:0.4103+0.6801i 0.4103-0.6801i当然，

4、如果增加以下两条语句，则可画出例5-1系统的零极点图，如图5-1所示。系统的零极点图，如图7-1所示.图5-1 零极点图【例5-2】已知离散系统的开环脉冲传递函数为：判断单位负反馈系统的稳定性。解：则可利用下面的MATLAB程序：%ex7_2.mnum0=5 4 1 0.6-3 0.5;den0=1 0 0 0 0 0;numc,denc=cloop(num0,den0);r=roots(denc);ii=find(abs(r)1);n1=length(ii);if(n10)disp(system is Unstable,with,int2str(n1),unstable pole);else

5、 disp(System is Stable);End执行结果显示：system is Unstable,with 1 unstable pole2.利用特征值判断系统的稳定性 对于线性定常系统称多项式为系统的特征多项式。其中，称为系统的特征多项式系数。令特征多项式等于零，即得系统的特征方程|sI-A|=sn+a1sn-1+an-1s+an0的根称为系统的特征值，即系统的闭环极点。当然判断系统的稳定性同样可利用特征值来判断。【例5-3】已知系统的状态方程为判断系统的稳定性。解：可利用以下的MATLAB程序。%ex5_3.mA=2.25-5-1.25-0.5;2.25-4.25-1.25-0.2

6、5;0.25-0.5-1.25-1;1.25-1.75-0.25-0.75;P=poly(A);r=roots(P);ii=find(real(r)0);n=length(ii);if(n0)disp(system is Unstable);else disp(System is Stable);end执行结果显示：System is Stable对于例5-3，利用下列命令可得同样的结果。r=eig(A);ii=find(real(r)0);n=length(ii);if(n0)disp(System is Unstable);else disp(System is Stable);end3.

7、用李雅普诺夫第二法来判断系统的稳定性l在高阶系统或者特征多项式中，当某些系数不是数值时，利用求闭环极点或特征值的方法来判断系统的稳定性是比较困难的。在这种情况下，利用李雅普诺夫第二法比较有效，尤其在系统含有非线性环节时更是如此。l线性定常连续系统l (5-2)l在平衡状态xe=0处，渐近稳定的充要条件是：对任给的一个正定对称矩阵Q，存在一个正定的对称矩阵P，且满足矩阵方程l ATPPAQ (5-3)l而标量函数V(x)=xTPx是这个系统的一个二次型形式的李雅普诺夫函数。l MATLAB提供了李雅普诺夫方程的求解函数lyap()，其调用格式为l P=lyap(A,Q)l式中，A,Q和P矩阵与式

8、（5-3）中各矩阵相对应。l更一般的，利用函数P=lyap(A,Q)可以求解下面给出的李雅普诺夫方程。l AP+PB=-Q (5-4)l对于离散系统的李雅普诺夫方程的求解函数为ldlyap().l【例5-4】设系统的状态方程为l其平衡状态在坐标原点处，试判断该系统的稳定性。l解：MATLAB程序为：l%ex5_4.mlA=0 1;-1-1;Q=eye(size(A);P=lyap(A,Q);li1=find(P(1,1)0);n1=length(i1);li2=find(det(P)0);n2=length(i2);lif(n10&n20)l disp(P0,正定，系统在原点处的平衡状态是渐进

9、稳定的);lelsel disp(系统不稳定);lend执行结果显示：P0,正定，系统在原点处的平衡状态是渐进稳定的课后作业lP216，715.2 控制系统的时域分析利用时域分析方法能够了解控制系统的动态性能，如系统的上升时间，调节时间，超调量和稳态误差都可以通过系统在给定输入信号作用下的过渡过程来评价。Matlab控制系统工具箱中提供了多种求取多种线性系统在特定输入下的时间响应曲线的函数，如表5所示。函数名功能函数名功能gensig()输入信号产生initial()求连续系统的零输入响应step()求连续系统的单位阶跃响应dinitial求离散系统的零输入响应dstep()求离散系统的单位阶

10、跃响应lsim()求连续系统对任意输入响应impulse()求连续系统的单位脉冲响应dlism()求离散系统对任意输入响应dimpulse()求离散系统的单位脉冲响应表5时域响应函数l1.任意信号函数l生成任意信号函数gensig()的调用格式为l u,t=gensig(type,Ta)l或 u,t=gensig(type,Ta,Tf,T)l其中,第一式产生一个类型为type的信号序列u(t)，周期为Ta，type为以下标识字符串之一：sin正 弦 波；square方 波；pulse脉冲序列；第二式同时定义信号序列u(t)的持续时间Tf和采样时间T。l例例5-5 生成一个周期为5秒，持续时间为

11、30秒，采样时间为0.1秒的方波。l解解 Matlab窗口中执行以下命令可得图5-2所示结果。lu,t=gensig(square,5,30,0.1);lplot(t,u)laxis(0,30,0.5,1.5)图5-2 方波信号l2.连续系统的单位阶跃响应l单位阶跃响应函数step()的调用格式为l y,x,t=step(num,den,t)l或 y,x,t=step(A,B,C,D,iu,t)l式中,t为选定的仿真时间向量;函数返回值y为系统在各个仿真时刻的输出所组成的矩阵；而x为自动选择的状态变量的时间响应数据。l 如只想绘制出系统的阶跃响应曲线，则可以由如下的格式调用此函数l step(

12、num,den,t)l step(A,B,C,D,t)l当 然，时 间 向 量 t也 可 以 省 略，此 时，MATLAB自动选择一个比较合适的仿真时间。l【例56】假设系统的开环传递函数为l试求该系统在单位负反馈下的阶跃响应曲线和最大超调量。l解：MATLAB程序为：l%ex5_6.mlnum0=20;den0=1 8 36 40 0;num,den=cloop(num0,den0);lt=0:0.1:10;y,x,t=step(num,den,t);plot(t,y)lM=(max(y)-1)/1)*100;disp(最大超调量M=num2str(M)%)l执行结果为：最大超调量M2.55

13、46%，单位阶跃响应曲线如图53中曲线所示。l另外，对于例56结果，也可利用step(num,den,t)命令得到。在该命令自动绘制的系统单位阶跃曲线的窗口中，首先在空白处单击鼠标右键，利用弹出菜单中的Characteristics Peak Response选项，获得系统阶跃响应曲线上的峰值标记点；然后单击此标记点即可获得该系统的响应峰值、最大超调量(%)和峰值时间分别为1.03、2.55和5.8，如图53所示。图5-3 例56的单位阶跃响应曲线l例例5-7 对于典型二阶系统l l试绘制出无阻尼自然振荡频率n=6，阻尼比分别为0.2,0.4,1.0,2.0时系统的单位阶跃响应曲线。解解 MA

14、TLAB程序为%Example5_7.mwn=6;zeta=0.2:0.2:1.0,2.0;figure(1);hold onfor k=zeta;num=wn.2;den=1,2*k*wn,wn.2;step(num,den);endtitle(Step Response);hold off 执行后可得如图5-5所示的单位阶跃响应曲线。从图中可以看出，在过阻尼()和临界阻尼()响应曲线中，临界阻尼响应应具有最短的上升时间，响应速度最快；在欠阻尼（）响应曲线中，阻尼系数越小，超调量越大，上升时间越短，通常取l为宜，这时超调量适度，调节时间较短。图5-4 例57的单位阶跃响应曲线l【例58】对例

15、57中典型二阶系统，绘制出l ，取2，4，6，8，12时的单位阶跃响应。l解：MATLAB程序为：l%Example5_8.mlw=2:2:12;zeta=0.7;lfigure(1);hold onlfor wn=wlnum=wn.2;lden=1,2*zeta*wn,wn.2;lstep(num,den);endltitle(Step Response);hold off图5-5 例58的单位阶跃响应曲线l3.离散系统的单位阶跃响应l 离散系统的单位阶跃响应函数dstep()的调用格式为l y,x=dstep(num,den,n)l或 y,x=dstep(G,H,C,D,iu,n)l式中，

16、n为选定的取样点个数，当n省略时，取样点数由函数自动选取，其余参数定义同前。例例5-9 已知二阶离散系统试求其单位阶跃响应。解：解：MATLAB程序为%Example5_9.mnum=2-3.4 1.5;den=1-1.6 0.8;dstep(num,den);title(Piscrete Step Response)执行后得如图5-7所示的单位阶跃响应曲线。图5-7 例59的单位阶跃响应曲线l例例5-10 对于多入多出系统对于多入多出系统l求单位阶跃响应。求单位阶跃响应。l解：这是双入双出系统，因此其阶跃响解：这是双入双出系统，因此其阶跃响应应有应应有4个。个。MATLAB程序如下：程序如下

17、：l%Example5_10.mlA=2.25-5-1.25-0.5;2.25-4.25-1.25-0.25;0.25-0.5-1.25-1;1.25-1.75-0.25-0.75;lB=4 6;2 4;2 2;0 2;C=0 0 0 1;0 2 0 2;D=zeros(2,2);lfigure(1);step(A,B,C,D)图5-7 例510的结果曲线l4.单位脉冲响应 l 单位脉冲响应函数impulse()和dimpulse()与单位阶跃函数step()和dstep()的调用格式完全一致，这里就不一一列写了。l5.系统的零输入响应l对于连续系统由初始状态所引起的响应，即零输入响应，可由函

18、数initial()来求得，其调用格式为l y,x,t=initial(A,B,C,D,x0)l或 y,x,t=initial(A,B,C,D,x0,t)l其中 x0为初始状态，其余参数定义同前。l同样对于离散系统的零输入响应函数dinitial()的调用格式为：ly,x=dinitial(A,B,C,D,x0)l或 y,x=dinitial(A,B,C,D,x0,n)l其中,n为取样点数，省略时由函数自动选取。l【例511】已知系统的状态空间表达式为：l以T0.5为采样周期，采用双线性变换算法转换成离散系统，然后求出离散系统的单位阶跃响应、单位脉冲响应及零输入响应（设初始状态x0=1 1 1

19、-1T）l解：MATLAB的程序为：l%Example5_11.mlA1=-1.6-0.9 0 0;0.9 0 0 0;0.4 0.5-5.0-2.45;0 0 2.45 0;lB1=1;0;1;0;C1=1 1 1 1;D1=0;lT=0.5;A,B,C,D=c2dm(A1,B1,C1,D1,T,tustin);lfigure(1);subplot(2,2,1);dstep(A,B,C,D),title(Discrete step response)lsubplot(2,2,2);dimpulse(A,B,C,D),title(Discrete impulse response)lsubpl

20、ot(2,2,3);x0=1;1;1;-1;dinitial(A,B,C,D,x0)laxis(0 6-0.5 2.5);title(Discrete Initial Response)lsubplot(2,2,4);z,p,k=ss2zp(A,B,C,D);pzmap(z,p);ltitle(Discrete Pole-Zero Map)图5-8 例511的曲线l6.任意输入函数的响应l连续系统对任意输入函数的响应可利用MATLAB的函数lsim()求取，其调用格式为l y,x=lsim(num,den,u,t)l或 y,x=lsim(A,B,C,D,iu,u,t)l其中 u为给定输入序列构

21、成的矩阵，它的每列对应一个输入，每行对应一个新的时间点，其行数与时间t的长度相等。其他用法同step()函数。l【例512】已知系统l试求周期为4s的方波输出响应。l解：利用函数rem()编写的MATLAB程序为：l%ex5_12.mlnum=2 5 1;den=1 2 3;lt=(0:0.1:10);period=4;u=1*(rem(t,period)=period./2);llsim(num,den,u,t);title(Square Ware Response)l执行后可得如图59所示的输出响应曲线。图5-9 例512的输出响应曲线l同样，离散系统对任意输入函数的响应，可利用dlsim

22、()函数求得，其基本调用格式为：ly,x=dlsim(num,den,u,n)l或y,x=dlsim(A,B,C,D,iu,u,t)l其中，n为取样点数。l【例513】对离散系统l求当输入为幅值1的方波信号时，系统的输出响应。l解：MATLAB程序为：l%Example5_13.mlnum=0.632;den=1-1.368 0.568;lu1=ones(1,50),-1*ones(1,50);lu=u1 u1 u1;lfigure(1),dlsim(num,den,u),title(Discrete System Simulation)l执行后可得如图510所示的仿真结果。图5-10 例51

23、3的仿真结果5.3 根轨迹法l根轨迹是分析和设计线性定常系统的一种非常简便的图解方法，特别适用于多回路系统的研究。l所谓根轨迹是指，当开环系统的某一参数从零变化到无穷大时，闭环系统的特征方程根在s平面上所形成的轨迹。一般地，将这一参数选作开环系统的增益K，而在无零极点对消时，闭环系统特征方程的根就是闭环传递函数的极点。l关于控制系统的根轨迹分析，在MATLAB控制系统工具箱中提供了几个函数，如表52所示。函数名说明pzmap()绘制系统的零极点图rlocus()绘制根轨迹rlocfind()求给定根的根轨迹增益sgrid()绘制连续系统的 网格根轨迹。zgrid()绘制离散系统的 网格根轨迹。

24、表52根轨迹函数l1.绘制系统的零极点图l pzmap()函数可绘制系统的零极点图，其调用格式为l p,z=pzmap(A,B,C,D)l或 p,z=pzmap(p,z)，l p,z=pzmap(num,den)l其中，列向量p为系统极点位置，列向量z为系统的零点位置。l 对于单变量系统，用pzmap()函数在复平面内可求出系统的的零极点；对单多变量系统，用pzmap()可求出系统的特征向量和传递零点。当带有输出变量时，可得到零极点位置，如需要可通过pzmap(p,z)绘制出零极点图，图中的极点用“”表示，零点用“o”表示。2.绘制系统的根轨迹l对于图5-11所示的负反馈系统，其特征方程可表示

25、为l l或 l l利用rlocus()函数可绘制出当开环增益K由0至变化时，闭环系统的特征根在s平面变化的轨迹，该函数的调用格式为l r,K=rlocus(num,den)l r,K=rlocus(num,den,K)l或 r,K=rlocus(A,B,C,D)l r,K=rlocus(A,B,C,D,K)l其中 返回值r为系统的闭环极点，lK为相应的增益。图511lrlocus函数既适用于连续系统，也适用于离散系统。rlocus（num,den）或rlocus（A,B,C,D）绘制系统根轨迹时，增益K是自动选取的，rlocus（num,den,K）或rlocus（A,B,C,D,K）可利用指

26、定的增益K来绘制系统的根轨迹。在不带有输出变量的引用函数时，rlocus()可在当前图形窗口中绘制出系统的根轨迹图。当带有输出变量的引用函数时，可得到根轨迹的位置列向量r及相应的增益K列向量，再利用plot(r,x)可绘制出根轨迹。l 在系统分析过程中，常常希望确定根轨迹上某一点处的增益值K，这时可利用MATLAB中的rlocfind()函数，在使用此函数前要首先得到系统的根轨迹，然后再执行如下命令l K,poles=rlocfind(num,den)l或 K,poles=rlocfind(A,B,C,D)l执行上述命令后，将在屏幕上的图形中生成一个十字光标，使用鼠标移动它至所希望的位置，然后

27、敲击鼠标左键即可得到该极点的位置坐标值poles以及它所对应的增益K值。l例例5-14 已知某负反馈系统的开环传递函数为l试绘制系统根轨迹，并分析系统稳定的K值范围。l解：解：MATLAB的程序为l%Example5_20.mlnum=1;den=conv(1,0,conv(1,1,1,2);lrlocus(num,den),lK,poles=rlocfind(num,den);执行以上程序，并移动鼠标到根轨迹与虚轴的交点处单击鼠标左键后可得如图5-12所示的根轨迹和如下结果：selected a point in the graphics windowlselected_point=l 0.

28、0000-1.4142ilK=l 6.0000lpoles=l -3.0000 l 0.0000+1.4142i l 0.0000-1.4142i 图5-12根轨迹图l 由此可见根轨迹与虚轴交点处的增益K=6，这说明当K 6时，系统不稳定；利用rlocfind()函数也可找出根轨迹从实轴上的分离点处的增益K=0.38,这说明当0 K 0.38时，系统为单调衰减稳定，当0.38 K subplot(2,1,1);semilogx(w,20*log10(mag)l subplot(2,1,2);semilogx(w,phase)l 在同一个窗口上同时绘制出系统的Bode图了，其中前一条命令中对幅值

29、向量mag求取分贝(dB)值。l 如果只想绘制出系统的Bode图，而对获得幅值和相位的具体数值并不感兴趣，则可以采用如下简单的调用格式l bode(num,den,)l bode(A,B,C,D,iu,）l或更简单地l bode(num,den)l bode(A,B,C,D,iu)例例5-16 已知二阶系统的开环传递函数为绘制出当a分别取0.2,0.4,0.6,0.8,1.0时系统的Bode图。解：解：当a取0.2,0.4,0.6,0.8,1.0时二阶系统的Bode图可直接采用bode()函数得到。MATLAB程序为：%ex5_16.m w=logspace(0,1);a=0.2:0.2:1.

30、0;figure(1);num=9;for k=aden=1,6*k,9;bode(num,den,w);hold on end grid;title(Bode plot);hold off 执行后得如图5-14所示Bode图。从图中可以看出，当0时，相角趋于0，当时，相角趋于-1800，当=n时，相角等于-900，此时的幅值也最大。图514 Bode图l3.幅值裕量和相位裕量l在判断系统稳定性时，常常需要求出系统的幅值裕量和相位裕量。利用MATLAB控制系统工具箱提供的margin()函数可以求出系统的幅值裕量与相位裕量，该函数的调用格式为l Gm,Pm,Wcg,Wcp=margin(num

31、,den)l或 Gm,Pm,Wcg,Wcp=margin(A,B,C,D)l式中 Gm和Pm分别为系统的幅值裕量和相位裕量，而Wcg和Wcp分别为幅值裕量和相位裕量处相应的频率值，即相位穿越频率 和幅值穿越 。l例例5-17 给定系统的开环状态空间表达式为 l求系统的幅值裕量和相位裕量，并画出Bode图。l解：MATLAB程序为：l%ex5_17.mlw=logspace(-1,1);lA=0 1 0 0;0 0 1 0;0 0 0 1;-62.5-213.8-204.2-54;lB=0;0;0;1;lC=1562 1875 0 0;D=0;lGm,Pm,Wcg,Wcp=margin(A,B,

32、C,D),bode(A,B,C,D,1,w);ldisp(幅值裕量=,num2str(20*log10(Gm),dB,相位裕量=,num2str(Pm),度。);ltitle(Gain margin=,num2str(Gm),Phase margin=num2str(Pm)执行后得如下数据及图执行后得如下数据及图5-155-15所示所示BodeBode图。图。lGm=l 4.4922lPm=l 23.0705lWcg=l 12.6466lWcp=l 5.8275 l幅值裕量=13.0492dB，相位裕量=23.0705度，即l12.6466时，相位曲线穿越180 线，幅值裕量为Gm=4.492

33、2=13.0492dB。l5.8275时，系统的幅值为1，相位裕量为Pm为23.0705度。图5-15 Bode 图l除了根据系统模型直接求取幅值和相位裕量之外，MATLAB的控制系统工具箱中还提供了由幅值和相位相应数据来求取裕量的方法，这时函数的调用格式为lGm,Pm,Wcg,Wcp=margin(mag,phase,w)l式中,频率响应可以是由bode()函数获得的幅值和相位向量，也可以是系统的实测幅值与相位向量，w为相应的频率点向量。l4.频率响应值lMATLAB也提供了直接的求取频率响应数据的函数freqresp(),其调用格式为l F=freqresp(num,den,sqrt(-1

34、)*w)l或 F=freqresp(A,B,C,D,iu,sqrt(-1)*w)l式中,F为频率响应，w为给定的频率范围向量。例【518】已知系统的传递函数为求在频率为0.110之间的频率响应。解：MATLAB程序为：l%ex5_18.mlw=logspace(-1,1);num=1 7 24 24;den=1 10 35 50 24;lF=freqresp(num,den,sqrt(-1)*w)l结果显示：lF=l 0.9892-0.1073il 0.9870-0.1176il 0.9843-0.1289il l 0.0266-0.0983il 0.0239-0.0919il 0.0212-

35、0.0857ilMATLAB也提供了直接求取LTI系统（LTI系统的定义参见第九章）在单个复数频率点处的频率响应数据的函数evalfr()，其调用格式为：lF=evalfr(sys,w)l式中，F为频率响应，w为给定的频率向量。l例【519】求例518所示系统在频率等于0.1的频率向量。l解：MATLAB程序为：l%ex5_19.mlw=0+0.1j;num=1 7 24 24;den=1 10 35 50 24;sys=tf(num,den);F=evalfr(sys,w);l结果显示：lFl0.9892-0.1073i5.系统的奈奎斯特图（Nyquist图）lNyquist图是根据开环频率

36、特性G(j)H(j)在复平面上绘制其极坐标图，利用开环系统的Nyquist曲线，可判断变换系统的稳定性。lNyquist稳定性判据可表示为：当从+变化时，Nyquist曲线G(j)H(j)逆时针包围(-1,j0)点的次数N，等于系统开环传递函数G(s)H(s)位于右半s平面的极点数P，即NP，则闭环系统稳定，否则闭环系统不稳定。nyquist()函数的调用格式为 Re,Im,w=nyquist(num,den)Re,Im,w=nyquist(num,den,w)Re,Im,w=nyquist(A,B,C,D)Re,Im,w=nyquist(A,B,C,D,iu)Re,Im,w=nyquist(

37、A,B,C,D,iu,w)其中 返回值Re,Im和w分别为频率特性的实部向量、虚部向量和对应的频率向量，有了这些值就可利用命令plot(Re,Im)来直接绘出系统的奈奎斯特图。l 当然也可使用下面的简单命令来直接绘出系统的奈奎斯特图。l nyquist(num,den,w)l或 nyquist(A,B,C,D)l更简单地 nyquist(num,den)l或 nyquist(A,B,C,D,iu)l它的使用方法基本同bode()函数的用法。l例例5-26 已知系统的开环传递函数为l绘制Nyquist图，并判断系统的稳定性。l解解 MATLAB命令为lnum=0.5;den=1 2 1 0.5;

38、lnyquist(num,den)l执行后可得如图5-16所示的曲线，由于Nyquist曲线没有包围(-1,j0)点，且P0，所以由G(s)H(s)构成的单位负反馈闭环系统稳定。图5-16 例520的Nyquist曲线l例【521】已知系统的开环传递函数为：l绘制系统的Nyquist曲线，并判别闭环系统的稳定性。l解：MATLAB命令为：lk=50;z=;p=-5 2;lnum,den=zp2tf(z,p,k);nyquist(num,den)l执行后可得如图521所示的Nyquist曲线。由图可知，lNyquist曲线按逆时针方向包围(-1,j0)点1次，而开环系统包含右半s平面上的一个极点

39、，所以以此构成的闭环系统稳定。图5-17 例521的Nyquist曲线l例【522】已知多环系统l其系统方框图如图518所示，试用Nyquist曲线判断系统的稳定性。图518 例522的图l解：先算出内环传递函数，l然后以G0(s)为开环传递函数绘制出Nyquist曲线，但这里不能直接采用奈氏判据，因为在前向通道上有一放大系数k=10,因此，奈氏判据中的临界点应改成(-1/k,j0)点。MATLAB程序为：l%ex5_22.mlk1=16.7/(0.8*0.25*0.0625);z1=0;p1=-1/0.8-1/0.25-1/0.0625;lnum1,den1=zp2tf(z1,p1,k1);

40、num,den=cloop(num1,den1);nyquist(num,den);lz,p,k=tf2zp(num,den);num1,den1=zp2tf(z,p,10*k);num,den=cloop(num1,den1);lz,p,k=tf2zp(num,den);pl执行后可得系统的闭环极点为：lPl 1.0e+002*l -0.1062+1.2113il -0.1062-1.2113il -0.0001l并得到如图519所示的Nyquist曲线。由图可知，Nyquist曲线不包围(-0.1,j0)点，而开环系统的三个极点均位于左半s平面，因此这个系统是稳定的。图5-19 例522的

41、Nyquist曲线l6.尼柯尔斯图(Nichols图)l 函数调用格式为l mag,phase,w=nichols(num,den,w)l或 mag,phase,w=nichols(A,B,C,D,iu,w)l可见该函数的调用格式以及返回的值与bode()函数完全一致，事实上虽然它们使用的算法不同,但这两个函数得出的结果还是基本一致的。但Nichols图的绘制方式和Bode图是不同的，它可由以下命令绘制lplot(phase,20*log10(mag)l 当然，Nichols图也可采用与Bode图类似的简单命令来直接绘制。例例5-23 已知单位负反馈的开环传递函数为试绘制Nichols图。解解

42、 MATLAB程序为%Example5_23.mnum=1;den=conv(1,0,conv(1,1,0.2,1);w=logspace(-1,1,400);mag,phase=nichols(num,den,w);plot(phase,20*log10(mag)ngrid 执行后可得如图5-20所示的Nichols图。图5-21 例523的Nichols图l7.离散系统的频率分析 l离散时间系统的频率分析也可以调用相应的MATLAB控制系统工具箱函数来完成，这些函数是以连续系统的函数名前加一字母d来命名的，例如离散时间系统的Bode图可以由dbode()函数求出，离散时间系统的Nyquis

43、t图可以由dnyquist()函数来求出，以及离散时间系统Nichols图可以由dnichols()函数来求出，其实在MATLAB的控制系统工具箱中这样的函数命名方式是相当普遍的，它们的调用格式与连续系统类似。例如，dbode()函数的调用格式为：mag,phase=dbode(num,den,Ts,w)l或 mag,phase=dbode(G,H,C,D,Ts,iu,w)l式中，num,den和G,H,C,D分别为离散系统传递函数和状态方程模型的参数。Ts为采样周期；iu为输入序号；w仍为频率向量。当w省略时，频率点在0 弧度之间自动选取。返回值mag,phase则分别为该系统的幅值和相位向

44、量，通过它们可以绘制出系统的频率响应曲线。当然这一函数的调用格式也可以简化成下面格式:l dbode(num,den,Ts,w)l dbode(G,H,C,D,Ts,iu,w)l更简单地 dbode(num,den,Ts)l或 dbode(G,H,C,D,Ts,iu)l在这种情况下，将直接绘制出离散时间系统的Bode图。l【例524】已知开环系统的离散时间状态空间表达式为：l试绘制出Bode图（设采样周期为Ts=0.1）l解：MATLAB命令为l%ex5_24.mlG=-1-2-2;0-1 1;1 0-1;H=2;0;1;C=1 2 0;D=0;ldbode(G,H,C,D,0.1,1);gr

45、idl如图521所示。图521 例524的Bode图l8.时间延迟系统的频率分析l带有时间延迟的连续控制系统的传递函数模型可以写成：l式中，为延迟时间常数；为不带有时间延迟的传递函数模型。由此可见，带有时间延迟的系统从某种意义上说，相当于在一个不带时间延迟的传递函数模型后面串接一个纯时间延迟环节 。l带时间延迟的状态方程模型可以写成：l l 函数的调用格式为l num,den=pade(,n)（55）l8.时间延迟系统的频率分析l带有时间延迟的连续控制系统的传递函数模型可以写成：l式中，为延迟时间常数；为不带有时间延迟的传递函数模型。由此可见，带有时间延迟的系统从某种意义上说，相当于在一个不带

46、时间延迟的传递函数模型后面串接一个纯时间延迟环节 。l带时间延迟的状态方程模型可以写成：（56）其中，矩阵和不含有时间延迟是类似的。带时间延迟的连续控制系统的频率响应可以由以下两种方法直接求得，即精确法和近似法。（1）首先求出不含时间延迟的传递函数模型 的幅频响应和相频响应,这样整个系统的幅频响应和 的幅频响应是一致的，而相频响应等于 的相位再减去 ，亦即（2）先将纯时间延迟环节 近似为有理函数（即传递函数类的分子分母多项式）的形式，再根据串联后总的传递函数来算出频率响应。l纯时间延迟环节 常采用 近似方法，它是1892年由法国数学家 提出的一种著名的有理近似方法，其表达式为：l式中，称为 近

47、似系数。它利用MATLAB控制系统工具箱中提供的函数pade()求得，该函数的调用格式为lnum,den=pade(tau,n)l式中，tau为延迟时间常数 ；n为要求拟合的阶数；调用该函数后将返回pade近似的传递函数模型num,den。一般情况下，取pade近似的拟合阶次为3或4就可以获得相当满意的精度。l【例525】系统的开环传递函数为：l采用精确法和近似法分别计算出系统的频率响应。l解：MATLAB程序为：l%ex5_25.mlnum1=1 1;lden1=conv(1,2,conv(1,2,1,2);lw=logspace(-1,2);tau=0.5;lm1,p1=bode(num1

48、,den1,w);lp1=p1-tau*w*180/pi;ln2,d2=pade(tau,4);lnumT=conv(n2,num1);denT=conv(den1,d2);lm2,p2=bode(numT,denT,w);lsubplot(2,1,1);lsemilogx(w,20*log10(m1),w,20*log10(m2),-)lsubplot(2,1,2);semilogx(w,p1,w,p2,-)l执行后可得如图522所示的Bode图，图中实线为系统采用精确法后的Bode图，虚线为系统采用近似法的Bode图。图522 Bode图9.利用频率响应数据辩识系统模型 对于某些比较复杂的

49、控制系统的数学模型往往没有办法从物理上得出，但可以通过适当的实验手段测试出系统的某种响应信息，如可以通过频率响应测试出系统的频率响应数据，或通过数据采集系统来测试出系统时间响应的输入与输出数据。有了系统的某种响应数据，就可以根据它来获得系统的数学模型。这种获得数学模型的过程称为系统辨识。在matlab信号处理工具箱中，给出了一个根据系统的频率响应数据来辨识系统传递函数模型的函数invfreqs()，该函数的调用格式为：num,den=invfreqs(F,w,m,n)式中，w为由频率点构成的向量；F为复数向量，其实部和虚部为辨识时用到的频率响应数据的实部和虚部数据；m和n分别为待辨识系统的分子

50、和分母阶次。返回的num和den分别为辨识出传递函数分子和分母的系数向量，即系统的传递函数模型。如果给出系统的幅频和相频phase响应数据，则可以由下面的方式来调用invfreqs()函数。lnum,den=invfreqs(mag.*exp(sqrt(-1)*phase),w,m,n)l除了由频率响应数据来辨识系统模型外，还可以根据阶跃响应及脉冲响应数据对系统的传递函数进行辨识，其具体的做法是，首先将阶跃响应数据或脉冲响应数据转化为响应的频率响应数据，然后再根据上面的方法来辨识原系统的模型。由脉冲响应数据求取频率响应数据的方法是很显然的，因为脉冲响应函数g(t)和频率响应函数G(j )满足下