第4章-控制系统的数学模型MATLAB实现

4.4系统的根轨迹分析4.5系统的频域分析4.6系统的性质分析4.7离散系统的分析第四章控制系统数学模型的MATLAB实现4.2系统的组合和连接4.1反馈控制系统的数学模型及其转换4.3系统的时域分析14.1

反馈控制系统的数学模型及其转换

控制系统的分析是系统设计的重要步骤之一

(1)在设计控制器前要分析系统的不可变部分，确定原系统在哪些方面的性能指标还不满足设计要求，有针对性的设计控制器；

(2)控制器设计完成后要验证整个闭环系统的性能指标是否满足设计要求。在控制系统基本理论和控制系统工具箱函数的基础上，利用MATLAB语言及其工具箱来解决控制系统的分析问题，包括系统模型的建立、模型的转换以及线性系统的时域分析、频域分析、根轨迹分析和系统的稳定性分析，为系统的仿真和设计做准备.下页上页2

为了对系统的性能进行分析首先要建立其数学模型，在MATLAB中提供了3种数学模型描述的形式：（1）传递函数模型tf(）（2）零极点形式的数学模型zpk()（3）状态空间模型ss()

本节首先介绍利用MATLAB提供的3个函数来建立系统的数学模型，下一节在此基础上介绍各种数学模型之间的相互转换。下页上页34.1.1

系统的数学模型(3+1)种格式：sys＝tf（num，den）功能：建立系统的传递函数模型说明：假设系统是单输入单输出系统（简称SISO），其输入输出分别用u（t），y（t）来表示，则得到线性系统的传递函数模型：

在MATLAB语言中，可以利用传递函数分子、分母多项式的系数向量进行描述。分子num、分母den多项式的系数向量分别为:这里分子、分母多项式系数按s的降幂排列。1．tf

传递函数模型下页上页4例4-1：已知系统的传递函数为:试建立系统的传递函数模型。解:程序如下num=[29];den=[13246];sys=tf(num,den)下页上页程序运行结果:Transferfunction:2s+9-----------------------------s^4+3s^3+2s^2+4s+65解:程序如下num=7*[23];den=conv(conv(conv([100],[31]),conv([12],[12])),

[538]);sys=tf(num,den)下页上页例4-2：已知系统传递函数如下应用Matlab语言建立系统的传递函数模型。程序运行结果:Transferfunction:14s+21------------------------------------------------------15s^7+74s^6+143s^5+172s^4+140s^3+32s^262．zpk

零极点形式的数学模型格式：sys＝zpk(z,p,k),

功能：建立零极点形式的数学模型说明：系统的传递函数零极点模型一般表示为：下页上页其中Zi（i＝1,2…,m）和Pi（i＝1,2…,n）分别为系统的零点和极点，K为系统的增益。[z]、[p]、[k]分别为系统的零点向量、极点向量和增益向量。7例4-3：已知系统传递函数如下应用Matlab语言建立系统的零极点形式模型。解:程序如下z=-4;p=[-1-2-3];k=5sys=zpk(z,p,k)上面的程序也可用下面的一行语句替换,sys=zpk([-4],[-1-2-3],[5])下页上页程序运行结果:Zero/pole/gain:5(s+4)-----------------(s+1)(s+2)(s+3)83．SS状态空间模型格式：sys＝ss（A,B,C,D），sys＝ss（A,B,C,D,T）功能：建立系统的状态空间模型，其中，T为取样时间说明：状态方程是研究系统的最为有效的系统数学描述，在引进相应的状态变量后，可将一组一阶微分方程表示成状态方程的形式。下页上页其中:

B为n×m维系统输入阵；C为l

×n维输出阵；D为l×m维直接传输阵。X为n×1维状态向量,U为m×1维输入矩阵;Y为l×1维输出向量;其中:A为n×n的系统状态阵，由系统参数决定。9例4-3a已知状态方程如下：

其中，

A=[12;34]，B=[56;78]，C=[14]，D=[69]。建立系统的状态空间模型。解：程序如下：A=[12;34]B=[56;78]C=[14]D=[69]sys=ss(A,B,C,D)下页上页结果如下：A=

x1x2x112x234B=u1u2x156x278C=

x1x2y114D=u1u2y169Continuous-timemodel.10或用语句

sys=ss([12;34],[56;78],[14],[69])来替代上述程序。显示结果相同。若用语句

sys=ss([12;34],[56;78],[14],[69]，5)显示结果和上述相同外，并且显示：

Samplingtime:5下页上页114．Simulink模型

不一定非要用数学表达式描述的数学模型对象才能对系统进行仿真，恰恰相反，就有不用数学表达式描述的数学模型对象，也能对系统进行仿真，这就是MATLAB特有的一种数学模型，即就是Simulink模型窗口里的动态结构图。

只要在Simulink工作窗里，按其规则画出动态结构图，再按规则将结构图的参量用实际系统的数据进行设置，就是对系统建立了数学模型，而且用Simulink模型窗口里的动态结构图来求系统的传递函数，还用来直接方便地对系统进行各种仿真，这一方法使用的极其广泛。124.1.2模型的转换

在进行系统分析时，往往根据不同的要求选择不同形式的数学模型，因此经常要在不同形式数学模型之间相互转换，下面介绍三种模型之间的相互转换函数。1、ss2tf函数将状态空间形式转换为传递函数形式格式：

[num,den]=ss2tf(A,B,C,D,iu)

说明：

ss2tf函数可以将状态空间表示通过

其中：iu用于指定变换所使用的输入量（个数）；num和den分别为传递函数的分子、分母多项式系数向量。ss2tf还可以应用离散时间系统，这时得到的是Z变换表示。转换为传递函数形式。下页上页13

例4-4：已知系统Σ(A,B,C,D)的系数矩阵是求取该系统相应的传递函数模型。C＝[110],D=0程序如下：

A=[200;041;004],B=[101]‘,C=[110],D=0,[num,den]=ss2tf(A,B,C,D,1)显示结果：num=01.0000-7.000014.0000den=1-1032-32→传递函数的分子→传递函数的分母其中，“1”指1个输入,可缺省不写下页上页142.ss2zp函数将系统的状态空间模型转换为零极点增益模型格式：[z,p,k]=ss2zp(A,B,C,D,iu)

式中,iu用于指定变换所使用的输入量（个数）.3.tf2ss函数将系统的传递函数模型转换为状态空间模型。格式：[A,B,C,D]=tf2ss（num，den）下页上页4.tf2zp

将系统的传递函数模型转换为零极点增益模型格式：[z,p,k]=tf2zp(num,den)5.zp2ss

将系统的零极点增益模型转换为状态空间模型。格式：[A,B,C,D]=zp2ss(z,p,k)6.zp2tf

将系统零极点增益模型转换为传递函数模型。

格式：[num,den]=zp2tf(z,p,k)15例4-5：已知系统的传递函数为应用MATLAB的模型转换函数将其转换为状态方程形式的模型。

解:num=[1836];den=[140.4391];[A,B,C,D]=tf2ss(num,den）运行程序显示结果:A=-40.4000-391.00001.00000B=10C=1836D=016应用MATLAB的模型转换函数将其转换为零极点形式的模型。例4-6：已知系统的传递函数为解：应用tf2ZP转换函数将传递函数转换为零极点模型.程序如下mun=[1836]den=[140.439150][z,p,k]=tf2zp(num,den)下页上页z=3.5000+1.3229i3.5000-1.3229ip=-39.5090-0.4455+1.8969i-0.4455-1.8969ik=1.0000177、由方框图建模控制系统工具箱中提供了方框图建模函数connect()

格式：

[A,B,C,D]=connect(a,b,c,d,q,iu,iy)

blkbuild

%从方框图中构造状态空间系统

功能：将方框图模型描述转换成状态空间模型。式中：a、b、c三个参数不需用户输入。

q矩阵规定了各框之间的连接关系。其每一行的第一个元素是框号，其余的元素依次是与该框连接的框号；

iu、iy分别表示输入、输出施加的框号。下页上页18例4-6：如下图所示，将由框图表示的系统转换成状态空间描述和传递函数程序如下：n1=1;d1=1;n2=0.5;d2=1;n3=4;d3=[14];n4=1;d4=[12];n5=1;d5=[13];n6=2;d6=1;n7=5;d7=1;n8=1;d8=1;nblocks=8;%指定块的总数blkbuild

%从方框图(blockbuild)

中构造状态空间系统

下页上页19q=[10000;21-6-7-8;32000;43000;54000;63000;74000;85000]%q矩阵表示框图的结构。由于前向通道有5个框，所以q矩阵有5列元素。共有8个框，每行第1个元素为框号，如第2个框与第1个框按1的关系连接，与第6、7、8个框按-1关系连接，…….依次类推。重画框图iu=[1];%输入施加于第1个框上iy=[8];%输出施加第8个框上[A,B,C,D]=connect(a,b,c,d,q,iu,iy)[num,den]=ss2tf(A,B,C,D,1)%转换成传递函数，式中”1”是指正反馈下页上页20显示结果如下：Statemodel[a,b,c,d]oftheblockdiagramhas8inputsand8outputs.A=-8.0000-2.5000-0.50004.0000-2.0000001.0000-3.0000B=0.500000C=001D=0num=00.00000.00002.0000den=1.000013.000056.000080.0000214.2

系统的组合和连接

所谓系统组合，就是将两个或多个子系统按一定方式加以连接形成新的系统。这种连接组合方式主要有串联、并联、反馈等形式。MATLAB提供了进行这类组合连接的相关函数。1.series系统的串联格式1：sys＝series（sys1,sys2），格式2：sys＝series（sys1,sys2,outputs1,inputs2）功能：用于将两个线性模型串联形成新的系统即

sys＝sys1*sys2说明：格式1：对应于SISO系统的串联连接。格式2：对应于MIMO系统的串联连接；其中sys1的输出向量为outputs1sys2的输入向量为inputs2下页上页222．Parallel----系统的并联格式1：sys=parallel(sys1,sys2)格式2：sys=parallel(sys1,sys2,in1,in2,out1,out2)功能：将两个系统以并联方式连接成新的系统，即sys=sys1+sys2。说明：并联连接时，输入信号相同，并联后其输出为sys1和sys2

这两个系统的输出之和。若用传递函数来描述，则系统输出

Y(S)=Y1(S)+Y2(S)=G1(S)U(S)+G2(S)U(S)=[G1(S)+G2(S)]U(S)

所以总的传递函数为G(s)=G1(s)+G2(s)。格式1：对应于SISO系统的并联连接。并联后其输出为sys1

和

sys2这两个系统的输出之和。格式2：对应于MIMO系统的并联连接。in1与in2指定了相连接的输入端，out1和out2指定了进行信号相加的输出端下页上页23例4-7已知两个线性系统，分别应用series和parallel函数进行系统的串并联连接。num1=[124];den1=[152];sys1=tf(num1,den1)%创建传函1num2=[16];den2=[171];sys2=tf(num2,den2)%创建传函2Sys_series=series(sys1,sys2)%两系统的串联，或用sys1*sys2Sys_parallel=parallel(sys1,sys2)%两系统的并联,或用sys1+sys2程序如下：下页上页24显示结果：

12s+4Transferfunction:-------------s^2+5s+2

s+6Transferfunction:-------------s^2+7s+1

Transferfunction:12s^2+76s+24--------------------------------s^4+12s^3+38s^2+19s+2

13s^3+99s^2+72s+16Transferfunction:------------------------------------s^4+12s^3+38s^2+19s+2下页上页传函sys1传函sys2两系统串联后的传函两系统并联后的传函253．feedback系统的反馈连接格式1：sys=feedback(sys1,sys2,sign)

格式2：sys=feedback(sys1,sys2,feedin,feedout,sign)

功能：

实现两个系统的反馈连接。说明：格式1：对于SISO系统，sys1表示前向通道传函，

sys2表示反馈通道传函,sign=1,正反馈.sign=-1,负反馈(默认值,可省略)

格式2：在已确立的MIMO系统sys1中,sys2表示反馈通道传函。其中feedin和feedout

分别指定了sys1的输入、输出端口号。最终实现的反馈系统与sys1具有相同的输入、输出端.sign含义同格式1

下页上页26例4-8已知前向环节和反馈环节的状态空间表达式的系数阵分别为试将前向环节的输入1和输出2与反馈环节构成负反馈系统。y1y2SYS2u2u1_SYS1u3y3程序如下:

A1=[10;01];B1=[11;01];C1=[13;20];D1=[10;25];A2=[-20;10];B2=[10]';C2=[01];D2=0;sys1=ss(A1,B1,C1,D1);sys2=ss(A2,B2,C2,D2);%建立两系统的状态空间模型feedin=1;feedout=2;sign=-1;sys=feedback(sys1,sys2,feedin,feedout,sign)下页上页由sys=feedback(sys1,sys2,feedin,feedout,sign)得：27显示结果如下：即为整个系统状态方程的系数矩阵下页上页a=x1x2x3x4x1100-1x20100x320-2-2x40010b=u1u2x111x201x325x400c=x1x2x3x4y1130-1y2200-2d=u1u2y110y225Continuous-timemodel.284.3线性系统的时域分析系统的时域分析是指输入信号采用单位阶跃或单位冲激函数，其响应是时间t的函数，称为时域响应。从时域响应可以获得系统的各个方面的性能。

1、impulse求连续系统的单位冲激响应。

格式1：impulse(sys)[Y,X,T]=impulse(sys)

格式2：impulse(sys,t)[Y,X]=impulse(sys,t)

格式3：impulse(sys,iu)[Y,X,T]=impulse(sys,iu)

格式4：impulse(sys,iu,t)[Y,X]=impulse(sys,iu,t)

说明：sys为tf(),zpk(),ss()中任一种模型。其中T是时间向量对于不带返回参数的该函数在当前窗口中绘制出响应曲线。对于带有返回参数的将不绘制曲线,其中Y是输出向量X是状态向量。t为用户设定的时间向量,即仿真时间。对于MIMO系统,iu表示第iu个输入到所有输出的冲激响应曲线.注意函数内的时间参数t写在sys之后.下页上页29例4-9：考虑如图所示的典型反馈控制系统结构，已知H(S)G(S)Gc(S)其中，

_解：G=tf(4,[1234]);Gc=tf([1-3],[13]);H=tf(1,[0.011]);G1=G*Gc;G2=feedback(G1,H);impulse(G2);%作整个反馈控制系统的冲激响应曲线,figure,impulse(G1)%作前向通道传函所代表的系统的冲激响应曲线,求系统的开环和闭环单位脉冲响应。下页上页30开环系统的单位冲激响应闭环系统的单位冲激响应运行结果如下：下页上页31注意：上述程序也可由plot()作图函数实现.G=tf(4,[1234]);Gc=tf([1-3],[13]);H=tf(1,[0.011]);G1=G*Gc;G2=feedback(G1,H);[y,x]=impulse(G2);plot(x,y)[y1,x1]=impulse(G1);figure,plot(x1,y1)下页上页32注意：若要求两个波形显示在同一图形窗口中,并加有图例说明.G=tf(4,[1234]);Gc=tf([1-3],[13]);H=tf(1,[0.011]);G1=G*Gc;G2=feedback(G1,H);impulse(G1)%作前向通道传函所代表的系统的冲激响应曲线

holdon%保持作图impulse(G2);%作整个反馈控制系统的冲激响应曲线Legend(‘y1’,‘y2’)%图例说明grid%生成网格下页上页332、step

求连续系统的单位阶跃响应。格式1：step(sys)[Y,X,T]=step(sys)格式2：step(sys,t)[Y,X]=step(sys,t)格式3：step(sys,iu)[Y,X,T]=step(sys,iu)格式4：step(sys,iu,t)[Y,X]=step(sys,iu,t)说明：step()中的参数意义和implse()函数相同。sys为tf(),zpk(),ss()中任一种模型。其中T是时间向量.对于不带返回参数的该函数在当前窗口中绘制出响应曲线。对于带有返回参数的将不绘制曲线,其中X是状态向量,Y是输出向量。t为用户设定的时间向量,即仿真时间。对于MIMO系统,iu表示第iu个输入到所有输出的冲激响应曲线.下页上页34例4-10考虑下面传递函数模型：试绘制其单位阶跃响应曲线。G=tf([172424],[110355024]);y=step(G)y=size(y)%可得y是101X1阶列向量t=[0:0.1:10]’;%t与y有相同的长度plot(t,y)%绘制传递函数模型的单位阶跃响应曲线grid解:程序如下:下页上页35用plot(t,y)所绘制的阶跃响应曲线注意:可以不定义t,直接用语句step(G)

绘制传递函数G的单位阶跃响应曲线例4-10图(a)下页上页36直接用step(G)所绘制的阶跃响应曲线G=tf([172424],[110355024]);step(G)grid注意:

横坐标最大值是0～

4.5,纵坐标仍然0～1

例4-10图(b)下页上页37G=tf([172424],[110355024]);t=[0:0.1:10]';

step(G,t)%注意t的位置在后grid直接用step(G)给图时,用户也可设定仿真时间t例4-10图(c)下页上页38例4-11：求下面的零极点模型的单位阶跃响应曲线。

解:程序如下:z=[-12];p=[-0.5-1.5-3-4-5];k=6;G=zpk(z,p,k);%建立零极点模型Step(G)%绘制零极点模型的单位阶跃响应曲线下页上页39系统的特征方程为

过阻尼系统的特征根为

（即有

>1）称

为阻尼比

二阶系统的阶跃响应下页上页40（2）欠阻尼系统的特征根为

（即0<

<1

）,称

为阻尼比

下页上页413）临界阻尼系统的特征根为（即阻尼比

=1）下页上页42例4-12：典型二阶系统传递函数为：试分析不同参数下的系统单位阶跃响应解：1、假设将自然频率固定为＝1，ζ＝0,0.1,…1,2,3,5。程序如下:wn=1;zetas=[0:0.1:1,2,3,5];t=0:0.1:12;fori=1:length(zetas)

Gc=tf(wn^2,[1,2*zetas(i)*wn,wn^2]);holdon

step(Gc,t);endholdoffζ---称阻尼比下页上页43显示结果如下:下页上页ζ=0ζ=0.1ζ=1…ζ=5由图可知,阻尼系数ζ从0值增加到5的过程,系统从等幅振荡到减幅振荡到非振荡变化44当ζ=0时,即R=0为等幅振荡当ζ=0.1时,R≠0为欠阻尼振荡放电过程step(Gc,‘g’,t);作图为绿色step(Gc,t)默认作图颜色为蓝色下页上页45当ζ=1时,为临界阻尼非振荡放电过程下页上页462、将阻尼比ζ的值固定在ζ＝0.55，自然频率变化范围为0.1～1程序如下:wn=[0.1:0.1:1];zetas=0.55;t=0:0.1:12;fori=1:length(wn)

Gc=tf(wn(i)^2,[1,2*zetas*wn(i),wn(i)^2]);holdon

step(Gc,t);endholdoff注意ζ＝0.55<1,此时电路为欠阻尼振荡放电过程下页上页47显示结果如下:从图中可看出，自然频率W越小，起振越慢；自然频率W越大，起振越快。下页上页483.initial

求连续系统的零输入响应。格式1：initial(sys,x0)[Y,X,T]=initial(sys,x0)格式2：initial(sys,x0,t)[Y,X,T]=initial(sys,x0,t)说明：

initial函数可计算出连续时间线性系统由于初始状态所引起的响应（故而称为零输入响应)。同前述一样，

t为用户设定的时间向量,即仿真时间，

T是返回的时间参数向量。

X0是初始状态值。注意：sys只能为ss()模型下页上页49例4-13：二阶系统当初始状态为x0＝[1;0]时，求系统的零输入响应。解:程序如下:a=[-0.5572-0.7814;0.78140];b=[1;0];c=[1.96916.4493];d=0;x0=[1;0];t=0:0.1:20;initial(a,b,c,d,x0,t)%绘制系统的零输入响应下页上页注意：如果将上述状态空间模型转换为传递函数模型：并设为sys,若用initial(sys,x0)来求零状态响应时,MATLAB就会显示如下错误信息:???Errorusing==>rfinputsOnlyavailableforstate-spacemodels.即指出:只能是状态空间模型504lsim

求任意输入信号时系统的响应

格式1：lsim(sys1,u,t)[Y,X]=lsim(sys1,u,t)

格式2：lsim(sys2,u,t,x0)[Y,X]=lsim(sys2,u,t,x0)

说明：

u为输入信号.

t

为用户设定的时间向量,即仿真时间.sys1为tf()或zpk()模型。

sys2为ss()模型。其中

x0为初始条件下页上页51程序如下：

t=0:0.1:5;u=4*exp(-3*t);num=[1,7,24,24];den=[1,10,35,50,24];

lsim(num,den,u,t),legend(‘系统曲线’,‘输入信号曲线’)grid例4-13考虑下面传递函数模型：试绘制当输入信号为u=4e-3t时，该系统的响应曲线。注意：图释语句和网格线语句均必须置于作图语句后，才能生效。5253下页上页4.4线性系统的根轨迹一、定义

根轨迹指当系统开环放大倍数K从0到＋∞时,闭环特征方程的根在复平面上的轨迹。是美国学者W.R.埃文斯在1948年提出的根轨迹方法。可见，根轨迹法是一种直接由开环传递函数求取闭环特征根的方法.二、意义利用根轨迹来分析系统的暂态和稳态性能。三、根轨迹的分类或根轨迹法是由系统的开环传递函数的零极点分布情况画出系统闭环特征根轨迹。54下页上页四、根轨迹方程如图所示，设开环传递函数为kG(S)=k×num/den,其闭环系统的传递函数为闭环系统的根轨迹方程（即闭环系统的特征方程）如下：五、根轨迹的画法规则如下：n为开环传递函数分母的阶数，m为分子阶数H(S)kG(S)551、pzmap

绘制系统的极、零点图。格式1：pzmap(A,B,C,D)[p,z]=pzmap(A,B,C,D)格式2：pzmap(num,den)[p,z]=pzmap(num,den)格式3：pzmap(p,z)说明：极点用“×”表示，零点用“o”表示。对于不带返回参数的将绘制零极点图。对于带有返回参数的将不作图,其中返回参数P为极点的列向量，z为零点的列向量。格式3

是将已知的零点z极点p绘制在复平面上。

MATLAB专门提供了绘制根轨迹的函数：

pzmap()绘制极零点图

rlocus()绘制根轨迹,

rlocfind()计算根轨迹的增益,下页上页56例4-15有连续系统要求绘制出零极点图。解:程序如下:num=[0.050.045];den=conv([1-1.80.9];[156]);[p,z]=pzmap(num,den)pzmap(num,den)下页上页得极点：p1=-2.0000，p2=-3.0000

p3=0.9000+0.3000i,p4=0.9000-0.3000i，零点z=-0.9000p1p2p3p4z根据图可知：系统是不稳定的572、

rlocus

求系统根轨迹。格式1：rlocus(num,den)rlocus(num,den,k)[R,K]=rlocus(num,den)[R,K]=rlocus(num,den，k)格式2：rlocus(A,B,C,D)rlocus(A,B,C,D,k)[R,K]=rlocus(A,B,C,D)[R,K]=rlocus(A,B,C,D,k)说明:k为用户设定的开环增益值，若省略机器自动生成.

对于不带返回参数的将绘制根轨迹，对于带有返回参数的将不作图。其中返回参数R是对应开环增益K的闭环极点的位置。

R的列数和增益K的长度相同，R的第m列元素是对于第m个开环增益K的闭环系统的根注意：格式2可绘出连续系统和离散时间SISO状态空间表示的根轨迹。下页上页583、rlocfind

计算根轨迹上指定点的开环根轨迹增益值k,并将该增益下所有的闭环极点显示出来。当这个函数启动起来之后，在图形窗口上出现要求用户使用鼠标定位的提示，这时用户用鼠标点击根轨迹上所要求的点后，将返回一个k值（此即为开环根轨迹增益值）,同时返回该k值下的所有闭环极点p的值，并将此闭环极点直接在根轨迹曲线上显示出来。

格式1：[K,poles]=rlocfind(A,B,C,D)格式2：[K,poles]=rlocfind(num,den)下页上页说明：式中，

A、B、C、D为开环系统的状态方程系数矩阵；

num、den为开环系统传递函数的分子、分母；

k为相对应的开环增益向量；

Poles为与开环增益k相对应的闭环极点。59下页上页格式：[K,poles]=rlocfind(A,B,C,D,P)

[K,poles]=rlocfind(num,den,P)也可通过指定开环极点p得到开环增益的向量K

。

式中，p为指定的开环极点；

A、B、C、D为开环系统的状态方程系数矩阵；

num、den为开环系统传递函数的分子、分母；

k为相对应的开环增益向量；

Poles为与开环增益k相对应的闭环极点。说明：求出根轨迹上离p点很近的一个根及所对应的增益K60下页上页解:程序如下:num=[0.050.045];den=conv([1-1.80.9],[156]);G=tf(num,den)rlocus(G)%绘制由G构成的闭环系统根轨迹图[K,Poles]=rlocfind(G)%在根轨迹上选择极点位置，得到开环增益K的值。注意:上述程序最后两行不能互换位置.例4-16：设开环传递函数G(s)如下所示，绘制闭环系统的根轨迹图61运行程序显示图(a)后，同时在命令窗口中显示”Selectapointinthegraphicswindow“selected_point=-1.9882+1.9938i

%用户在根轨迹上任意选择的点K=484.7197%为上述选择的点对应的增益Poles=-4.3122，1.1428+2.0179i，1.1428-2.0179i，-1.1734Poles为根轨迹上增益K所对应的所有极点。得到图（b)图（a)根轨迹图图（b)根轨迹上增益K所对应的极点下页上页p1p2p3p4z62例4-17图4-17某闭环系统的开环传递函数为利用MATLAB命令画出该系统的根轨迹。程序如下：num=[11];den1=conv([1,0],[1-1]);den=conv(den1,[1416]);rlocus(num,den)%绘制根轨迹，如图所示。下页上页63p=1.5i;%p为离根轨迹较近的虚轴上的一个点。

[k,poles]=rlocfind(num,den,p)％求出根轨迹上离p点很近的一个根及所对应的增益k和其它三个根。

显示结果:k=22.5031，

poles=-1.5229+2.7454i-1.5229-2.7454i

0.0229+1.5108i

0.0229-1.5108ip1=1.5108i;%再以上述极点的虚部作为新的一点[k,poles]=rlocfind(num,den,p1)

％求出根轨迹上离p1点很近的一个根及所对应的增益k和其它三个根显示结果:k=22.6464，

poles=-1.5189+2.7382i-1.5189-2.7382i

0.0189+1.5197i

0.0189-1.5197i

下页上页64用同样的方法可得到根轨迹由左半平面穿过虚轴时的增益和四个根如下:k=35.6853

poles=0.0000+2.5616i0.0000-2.5616i-1.5000+1.7856i1.5000-1.7856增益的另一个临界值为35.6853,

由此可得增益的稳定范围为23.3160<k<35.6853再以此根的虚部为新的根，重复上述步骤，逐步地使极点向虚轴接近，几步后可得到下面的结果：k=23.316，poles=-1.5000+2.7040i-1.5000-2.7040i0.0000+1.5616i

0.0000-1.5616i这就是根轨迹由右半平面穿过虚轴时的增益及四个根。这时增益的临界值为23.3160，654.5系统的频域分析

频域分析法是利用系统开环的奈氏图、波特图、尼氏图分析系统的性能，如系统的稳态性能、动态性能、稳定性。一、幅相频率特性(奈氏图)66

奈氏图判定系统稳定的充要条件：如果开环系统有P个极点在右半平面，相应于频率ω从-∞→+∞变化时，开环频率特性G(jω)H(jω)曲线逆时针方向环绕（－1，j0）点的次数N等于右半根平面内的开环系统的极点数P，那么闭环系统就是稳定的，否则是不稳定的。下页上页同时奈氏判据指出，如果系统在开环状态下是稳定的，闭环系统稳定的充要条件是：它的开环幅相频率特性曲线不包围(-1，j0)点。反之，若曲线包围(-1，j0)点，则闭环系统将是不稳定的。若曲线通过(-1，j0)点，则闭环系统处于稳定边界。671、nyquist

求连续系统的Nyquist曲线格式1：nyquist(sys)[re,im,w]=nyquist(sys)格式2：nyquist(sys,w)[re,im,w]=nyquist(sys,w)格式3：nyquist(sys,iu,w)[re,im,w]=nyquist(sys,iuw)说明：对于不带返回参数的将绘制Nyquist曲线。对于带有返回参数的将不绘制曲线,其中，返回参数re、im分别为开环G(jw)在各频率点的实部和虚部，即：re=Re(G(jw)),im=Im(G(jw

))。

sys为tf(),zpk(),ss()中任一种模型。

w设定频率范围，省略时由机器自动产生。

对于MIMO系统，iu表示用第iu个输入变量来绘制系统的Nyquist曲线。返回参数为第i个输出变量针对第j个输入变量的频率响应的实部和虚部，即re(i,j,:)和im(i,j,:)。

68例4-17:系统开环传递函数绘制系统的Nyquist图，并讨论其稳定性下页上页程序如下：G=tf(1000,conv([132],[15]));nyquist(G);%绘制奈奎斯特图Axis(‘square’)%将两个坐标轴设定相等的长度（但刻度的增量可以不一定相等），即绘制正方坐标下的图形，通常在两个坐标轴具有相同比例时使用，如画园时使用。语句中单引号和园括号都可省略不写。运行程序显示如下的图形。69由图可知，Nyquist图曲线包围了点（－1，j0），而原开环系统中极点分别为p1=-1,p2=-2,p3=-5,没有不稳定极点，从而可以得出结论，闭环系统是不稳定的。下页上页70或者说：奈奎斯特稳定性判据（又称奈氏判据）

1、对于开环稳定的系统，闭环系统稳定的充分必要条件是开环系统的奈氏曲线G(jω)H(jω)不包围(-1,j0)点。反之，则闭环系统是不稳定的。

2、对于开环不稳定的系统，有p个开环极点位于右半s平面，则闭环系统稳定的充分必要条件是,当ω从-∞→∞变化时，开环系统的奈氏曲线G(jω)H(jω)逆时针包围(-1,j0)点p次。下页上页712、nichols

求连续系统的Nichols（尼克尔斯）频率响应曲线。格式1：nichols(sys)[re,im,w]=nichols(sys)格式2：nichols(sys,w)[re,im,w]=nichols(sys,w)格式3：nichols(sys,iu,w)[re,im,w]=nichols(sys,iuw)说明：nicholsh函数的输入变量定义与nyquist相同

iu是指第iu个输入，w是频率点构成的向量

3bode求连续系统的Bode（伯德）频率响应。格式1：bode(sys)[mag,phase,w]=bode(sys)格式2：bode(sys,w)[mag,phase,w]=bode(sys,w)格式3：bode(sys,iu,w)[mag,phase,w]=bode(sys,iuw)说明：bode函数的输入变量定义与nyquist相同Bode图可用于分析系统的增益裕度、相位裕度、增益、带宽以及稳定性等特性。mag和phase分别是幅值和相位数组。iu表示从系统第iu个输入到所有输出的Bode图

下页上页72试用MATLAB绘制出不同和的伯德图。解：wn=1;zet=[0:0.1:1,2,3,5];holdonfori=1:length(zet)num=wn^2;den=[1,2*zet(i)*wn,wn^2]

bode(num,den)endholdoff当阻尼比比较小时，则系统的频域响应在自然频率附近将表现出比较强的振荡，该现象称为谐振。例4-18：考虑二阶系统传递函数模型下页上页1、为固定值，变化时732、为固定值，变化时

当自然频率的值增加时，伯德图的带宽将增加，该现象使得系统的时域响应速度变快。解：wn=[0.1:0.1:1];zet=0.707;holdonfori=1:length(wn)

num=wn(i)^2;den=[1,2*zet*wn(i),wn(i)^2];

bode(num,den)endholdoff下页上页744、margin求取给定线性定常系统的幅值裕量和相角的裕量格式1：margin（sys）格式2：margin（mag,phase,w）格式3：[Gm,Pm,Wcg,Wcp]=margin（sys）

[Gm,Pm,Wcg,Wcp]=margin（mag,phase,w）

margin函数可从频率响应数据中计算出幅值裕度(不是db)、相角裕度和剪切频率。说明：格式1画出bode图，并标注幅值裕度和对应频率，相角裕度和对应频率。格式2由给定的幅值mag、相位phase及频率w画出bode图格式3：不画图，返回幅值裕度Gm和对应频率Wcg,相角裕度Pm和对应频率Wcp

。下页上页75例4-19：考虑如下系统模型：求它的幅值裕度和相角裕度，并求其闭环阶跃响应。例3-20系统的数学模型为：求出系统的幅值裕度与相角裕度。下页上页764.6线性系统的性质分析

4.6.1线性系统稳定性分析1.直接判定方法首先求出系统的所有极点，当其实部大于零，则系统为不稳定系统，否则称为稳定系统。若极点的实部等于0的，则系统称为临界稳定系统。对于传递函数模型tf（num，den）,利用求根函数roots（den）来求极点。对于状态方程模型SS（A,B,C,D）利用求特征值函数eig（A）来求特征值。这样根据极点或特征值直接判定系统的稳定性了。例4-21假设系统的传递函数模型为判断系统的稳定性。

下页上页77

1)routh

构造系统的Routh表格式:[rtab,msg]=routh(den)说明：其中den是系统的分母多项式向量，rtab是构造的Routh表矩阵，msg为字符串型变量，返回有关信息。系统不稳定极点的数目等于所产生的Routh表中第一列元素的符号变化次数。Routh（）函数

2间接判定方法

直接求取高阶代数方程的根是一件很困难的工作，所以出现了一些判定给定系统稳定性的间接方法，本节介绍其中两种方法，即Routh-Hurwitz（劳斯-霍尔维茨）判定方法和Lyapunov判定方法，它们分别适用于传递函数模型和状态空间模型。例4-22：同例4-20，应用的Routh表判断系统的稳定性。下页上页78

从运行结果可知Routh表第一列没有符号的变化，所以系统是稳定的。2）hurwitz(霍尔维茨)构造Hurwitz矩阵。格式:H＝hurwitz（den）说明：H为构造的Hurwitz矩阵，den为系统的分母多项式D(S)，D(S)=a0Sn+a1Sn-1+a2Sn-2+……+an-1S1+

anH矩阵

hurwitz()函数下页上页793)posdef

判定矩阵的正定性格式：[key,sdet]=posdef(A)说明：其中key返回矩阵A正定性的标记，若key＝1则表示该

矩阵A为正定矩阵，否则矩阵A为非正定矩阵。sdet返回各个左上角子矩阵的行列式值。

posdef（）函数例4-23：考虑例4-20中的系统模型由运行结果可知Hurwitz阵是正定的，因此系统稳定。

用Hurwitz(霍尔维茨)判据判定系统的稳定性。下页上页804)Lyap

解Lyapunov方程格式：X＝lyap（A,B,C）X＝lyap（A,C）说明：lyap(A,B,C)求解矩阵方程AX+XB=－C的解X；lyap(A,C)求解矩阵方程AX+XA'=-C的解X。例4-24：考虑一个状态方程模型

利用Lyapunov判据判定出系统的稳定性。解：选定一个正定的对角矩阵w=diag（1,2,3,4）Lyapunov方程的解V不是正定矩阵，因为其左上角子矩阵的行列式值都为负，故而由Lyapunov判据可知，系统不稳定。下页上页814.6.2线性系统的能控性和能观性

1、能控性n阶系统的完全能控性只取决于状态方程中的（A,B）矩阵，系统完全能控的充要条件是能控矩阵Tc满秩即：rank(Tc)=n,Tc=[B，AB，…An-1B]能控矩阵Tc由ctrb（）函数自动产生出来，其调用格式为：Tc＝ctrb（A，B）2