机电系统的计算机辅助分析

3.机电系统的计算机辅助分析与设计3.1机电系统的数学模型及其转换方法机电系统计算机仿真与辅助设计是建立在机电系统数学模型基础之上的。对于各类机电系统，利用仿真手段对其进行分析与设计，首先就需要建立相应的系统数学模型，此后，就需要研究如何将系统的数学模型转变为适合于计算机进行分析计算的仿真模型，即数值算法模型。在此基础上，即可通过对数学模型的求解分析，实现对系统动静态特性的分析与设计。显然，进行上述工作的重要基础就是系统的数学模型。因此本章首先介绍系统的几种典型数学描述，然后介绍各种数学模型之间的相互转换，以及系统环节不同形式的相互连接的MATLAB实现。3.1.1连续系统的数学描述连续系统的数学模型通常可以用微分方程、传递函数、状态空间表达式三种形式对系统加以描述。下面将简单对这几类数学模型加以回顾，同时给出MATLAB对它们的表示方法。.系统的微分方程形式模型一个系统的动态特性通常可用高阶微分方程加以描述，因此描述一个系统最常用的数学模型就是微分方程的形式。假设连续系统为单入单出(简称SISO)系统，其输入与输出分别用u(t)、y(t)加以表示，则描述系统的高阶微分方程为：dnydn-lydn-2y dy(3-1) 4 t+% v+…+凡1 h%ydtn1dtn-X2dtn-27dt(3-1)dn~xudn~2u—C. - +C2 7+…+CUdtn-x2dt''-2其初始条件为：y(t0)=y0,夕(。)=夕0，…，“(％)=“()，如果引入微分算子p= 贝式可以写作:p”y+%p”"yH Fan_xpy+any=cxpn~xu+qp〃，4 FcnuTOC\o"1-5"\h\z〃 ”一i即Xa>'-jPJy=YCn-iP，Uj=0 i=0对上式稍加整理并令〃o=1，可以得到Z*p'）= （3-2）7=0.系统传递函数形式模型1）传递函数模型对（3-1）式等号两边取拉氏变换，并假设y与u的各阶导数的初值均为零，则存在s"y（s）+6s"Ty（s）+…+*sY（s）+a“y（s） 口3）=clsn-lU（s）+c2sn~2U（s）+---+cnU（s）式中：Y（s）——输出y（f）的拉氏变换；U（s）——输入的拉氏变换。从而（3-1）式所描述的系统的传递函数为G（s）=，（s）=js"1+QS”2+...+c“_]S+c“ （34）U（s）sn+a15n-1+a2sn~2H \-an_ts+an对照（3-2）式与（3-4）式可以清楚地看出，当描述系统的微分方程的初始值为零时,用算子p所表示的式子与传递函数G（S）在形式上完全相同。传递函数是经典控制论描述系统的数学模型之一，它表达了系统输入量和输出量之间的关系。它只与系统本身的结构、特性和参数有关，而与输入量的变化无关,传递函数是研究线性系统动态响应和性能的重要手段与方法。在MATLAB语言中，可以利用分别定义的传递函数分子、分母多项式系数向量方便地对其加以描述。例如对于（3-4）式，系统可以别定义传递函数的分子、分母多项式系数向量为num=[cic2-Cn-icn]den=[1a\上…%-i%]这里分子、分母多项式系数向量中的系数均按s的降幕排列，由于传递函数G⑸的最高次项系数为I.所以分母多项式系数向量de”中第一个元素为1。在MATLAB5.x中，可以用tf来建立传递函数的系统模型，其基本格式为sys=tf(num,den)例3-1已知系统传递函数为G(s)=—— s4+3.J+a2s3+4S+6解：可以利用MATLAB将上述系统模型表示出来，并将其建立在工作空间(Workspace)中，写成num=[29]；den=[l2346]；printsys(num,den,*s*)执行上述语句后得到num/den=2s+9sA4+3sA3+2sA2+4s+6对于MATLAB5.x可以利用tf直接建立系统模型，即num=[29]；den=[13246J；model=tf(num,den)执行上述语句后得Transferfunction：2s+9sA4+3sA3+2sA2+4s+6由上例可以看出，对于MATLAB的不同版本建立的传递函数有所不同，但结果是一致的。下面考虑一种较复杂的传递函数形式。例3・2已知系统传递函数为--(3s+iXs+2)2(5s3+3s+8)解：可以利用MATLAB建立其相应的传递函数系统模型。num=7*[23]；den=conv(conv(conv([1,0,0],[3,1]),conv([12],[1,2])),[5,0,3,8])；model=tf(num,den)

运行结果Transferfunction：14s+2115sA8+65sA7+89sA6+83sA5+152sA4+140sA3+32sA2在这里使用了conv函数，它是MATLAB中的标准函数，用来求取两个向量的卷积分。因此当两个多项式相乘，需两个多项式系数向量相乘时，就可利用conv函数来进行。conv函数允许多重嵌套，由上例已清楚地看到这一点。对于连续时间系统可以用传递函数对其加以表示，而对于离散时间系统则应采用脉冲传递函数对其进行描述。脉冲传递函数一般可表示为关于z的降基多项式分式形式，即G(z)= =JZ1"+C吁遥mT+…+C]Z+Co~u(z)anzn+an_lzn-l+---+a,z+a0在MATLAB5.0中，对于离散系统同样可以建立相应的系统模型，其基本格式为num=[cm,cm.i，…，Ci，c0]：den=[an,a„.i，…，a”a。]；sys=tf(num,den.T)其中T为系统采样周期。2)系统的零极点形式模型系统的传递函数还可表示成另一种形式，即零极点形式。这种形式的系统传递函数比标准形式传递函数更加直观，可清楚地看到系统零极点分布情况。系统的零极点模型一般可表示为G(s)=KG(s)=K(s-ZiXs-Z2>Ys-Zm)(s-Pl\s-p2)-\s-pn)(3-5)其中，z,G=L2…,加)和Pj。=1,2,)分别为系统的零点和极点，Zi其中，z,G=L2可以是实数也可是复数，K为系统增益。MATLAB可以使用zpk函数建立零极点形式的系统模型.其基本格式为sys=zpk([z].[p],[k])其中，[z]、[p]、[k]分分别为系统的零极点和增益。3.系统的状态空间表达式状态方程是研究系统的最为有效的系统数学描述，不论是单入单出系统还是多入多出(简称MIMO)系统，若可用一组一阶微分方程对其加以表示，在引进相应的状态变量后，则可将这一组一阶微分方程写成紧凑形式，即状态空间表达式。通常一个线性定常系统可以表示为心AX+8U 36)Y=CX+DU其中上述第一式由n个一阶微分方程构成，称为系统的状态方程表达式，下式由/个线性代数程组构成，称为系统的输出方程。X为n维状态向量；U为机维输入向量；y为/维输出向量；A为"X"维的系统状态阵，由系统的参数所决定；B为nXm维的系统输入阵;C为IX”维输出阵；D为IXm维直接传输阵。应用MATLAB可以方便地表示系统的状态方程模型，只要按照矩阵输入方式建立相应的系统系数阵即可，考虑SISO系统，容易在MATLAB工作空间中建立系统的系数阵，形式为A=[〃l1。12…。；”21。22…。2〃；•••；an\an2***ann]\b2；•••；bn]；C=[clc2…由；D=d；当然，也完全可以在MATLAB工作空间中建立MIMO系统的系数阵。根据系统状态方程的系数阵，可在MATLAB中建立相应的系统模型，其基本格式为sys=ss(A,B.C.D)上述是MATLAB5.x中的格式，在MATLAB4.2中其格式为printsys(A,B,C,D)对于离散系统，其状态空间表达式可表示成X(%十1)=AX(k)+BU(k)Y(k)=CX(it)+DU(k)在MATLAB5.x中同样也可建立相应的系统模型，其格式为sys=ss(A,B,C.D,T)其中，T为系统采样周期。在MATLAB4.2中也可建立类似的系统模型。1.2系统模型的相互转换在前一节中已经介绍了描述一个系统的数学模型主要有三种形式；微分方程模型、传递函数模型(包括零极点形式模型)、状态空间模型等。显然在这些不同形式的系统数学模型中存在着内在的联系，虽然它们外在形式不同，但它们的实质内容是等价的。人们在进行系统分析研究时，往往根据不同的要求选择不同形式的系统数学模型，因此研究不同形式的数学模型之间的转换具有重要意义。.机电系统模型向传递函数或零极点增益形式的转换1).状态空间方程向传递函数形式的转换系统的状态空间方程可表示为X=AX+BUY=CX+DU据此，可以得到等效的系统传递函数模型。G(s)=^=C(s/-A)'B+Ddenys)显然在进行这种变换过程中，求取(sI-A)阵的逆比较困难，然而MATLAB有一系列的函数可以完成各种变换，其中就包括进行这种变换的ss2tf函数，其基本格式为[num,den]=ss2tf(A.B,C,D,iu)利用该函数即可实现将状态空间方程转换为传递函数的形式，iu用于指定变换所使用的输入量。为了获得传递函数的系统形式，还可以采用下述方式进行，即Gl=ss(A,B,C,D)；G2=tf(Gl)可以证明，由给定的状态方程模型转换为传递函数形式其结果是惟一的。例3-3 已知连续系统工(A,B,C,£»的系数矩阵是'2001 rfA=041,B=0,C=[11o],£)=0004j [1求取该系统相应的传递函数模型。解：应用MATLAB的ss2tf函数可以方便实现这种转换。A=[200；041；004]；B=[l01])：C=[l10]；D=0；[numden]=ss2tf(A,B,C>D.1)G=tf(num.den)运行上述语句即可实现要求的转换，利用G=tf(num,den)将转换后的传递函数分子，分母多项式系数向量构造成传递函数系数模型，下面给出了运算结果；num=01.0000-7.0000 14.0000den=1-10 32 -32Transferfunction：sA2-7s+14sA3-10sA2+32s-32对于多输入系统，应用ss2tf函数可以进行指定要求的模型转换。例3-4已知连续系统软£(A,B,C,D)的系数矩阵是-2.25 -5 -1.25-0.5'46A=2.25-4.25-1.25-0.250.25-0.5-1.25 -1，B=2242,C=[0202],D=01.25-1.75-0.25-0.7502求取系统相应的传递函数模型。解：该系统为一双输入系统，利用ss2tf函数可以进行指定输入的状态方程向传递函数的转换。ss2tf(A,B,C,D,T)函数中的T,指定了要转换的系统所对应的输入信号，下面给出了相应于第二输入信号的转换方法。A=[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]；B=[46；24；22；02]；C=l0202j；D=[00]；T=2；[num2,den2]=ss2tf(A,B,C,D,T)G=tf(num2,den2)运行结果num2=0 12.0000 32.000037.000 17.0000den2=1.00004.00006.2500 5.2500 2.2500Transferfunction：12sA3+32sA2+37s+17sA4+4sA3+6.25sA2+5.25s+2.25当T=1时，将指定第输入信号所对应的系统转换，下面给出了相应的转换结里y|\»Onuml=0 4.000014.0000 22.0000 15.0000denl=1.0000 4.0000 6.2500 5.2500 2.2500Transferfunction：4sA3+14sA2+22s+15sA4+4sA3+6.25sA2+5.25s+2.25.模型向零极点形式的转换系统的零极点模型实际上是传递函数模型的另一种形式，也是对系统进行分析的一类常用模型，因此在MATLAB中，也提供了实现将各类系统模型转换为零极点形式模型的函数。其基本格式为[z,p.k]=ss2zp(A,B,D,IU)[z,p,k]=tf2zp(num,den)Gzp=zpk(sys)上述第一式是将以状态方程形式给出的模型根据指定研究的输入，转换为零极点模型形式；第二式是将传递函数形式给出的模型转换为零极点形式；第三式十分简洁,用该函数可将非零极点形式的模型转换为零极点系统模型。例3-5已知系统状态方程的系数阵为A=] 0,B=0>C=[1.52],D=0试将其转换成零极点模型形式。解：利用转换函数[z,P,k]=ss2zp(A,B,C.D)即可实现所希望的转换，运行结果为L-1.3333P=kJ1.5000也可利用Gzp=zpk(sys)函数进行转换，即sys=ss(A,B»C,D)Gzp=zpk(sys)运行结果Continuous-timesystem.Zero/pole/gain；2(s+0.5)(s十2)(s+1)显然，两种方法所得到的结果是相关的。.系统模型向状态方程形式的转换前面介绍了系统模型转换为传递函数的方法，同样也可以利用MATLAB函数实现所需要的系统模型向状态方程的转换。其基本格式为[A,B,C,D]=tf2ss(num,den)[A,B,C,D]=zp2ss(Z,p,k)syss=ss(sys)例3・6已知系统传递函数为G(.)=-+40.4,V2+391.V+150应用MATLAB的模型转换函数将其转换为状态方程形式的模型。解：应用tf2ss转换函数很容易实现所要求的转换。但应注意，山于其状态变量选取的不同，转换结果是不惟一的。下面给出了应用tf2ss函数的转换程序与结果。num=[18,36]；den=[140.4391150]；[A,B,C,D]=tf2ss(num,den)运行结果为A=-40.4000 -391.0000 -150.00001.0000 0 00 1.0000 0B=01836D=利用［A,B,C,D］=zp2ss(z,p,k)可以将零极点形式给出的模型转换为状态方程。对于sys=ss(sys),可以将任意ED系统模型转换为状态方程。例3・7考虑例3・6,将程序进行修改。解：修改程序如下。num=L18,36］；den=L140.4391150］；G=tf(num»den)Sys=ss(G)运行结果A=-40.4000-391.0000-150.0000TOC\o"1-5"\h\z1.00000 0 00 1.0000 0B=C=0 18 36D=.1.3系统状态方程的变换与实现.状态方程的相似变换在前一节已经述及将系统的状态方程转换为传递函数时，其结果是惟一的；然而由于状态变量的选取不同，将系统的传递函数转换为状态方程则不能保证结果的惟一。换言之，对于一个系统的传递函数，可以存在着众多的状态方程实现。这样,对于同一个系统有着不同的状态方程的描述，因此也就存在着它们之间的相互变换,MATLAB中的ss2ss函数可以实现对系统状态方程的相似变换。其基本格式为GT=ss2ss(G»T)其中，G为原系统的状态方程模型，T为非奇异变换阵。该变换实现下列相似变换Z=7x,Z=AtU+BtUJ=Gz+DfU其中At=TAT',Bt=TBT',Ct=TCT',Dt=TDT'下面通过一个例子来说明这种变换。例3-8已知系统状态方程为试应用ss2ss函数进行状态方程的相似变换。解：由于变换阵7可以任意选择，只要保证其非奇异即可。在此选择单位反对角阵作为变换阵，下面给出了变换方法与结果。A=[10-42；800；020]；B=[200]'；C=[0.5-0.4375 0.4375]；D=0；G1=ss(A,B,C,D)T=fliplr(eye(3))；GT=ss2ss(Gl,T)得到结果为a=xlx2x3xl02.000000x2008.00000x32.00000-4.0000010.00000TOC\o"1-5"\h\zb=ulxl 0x2 0x32.00000c=xl x2 x3yl0.43750 0.43750 0.50000d=ulyl 0.规范型状态方程的实现在MATLAB中提供给用户一个状态方程的规范实现函数cannon,以进行ETI系统模型sys的规范状态空间表达式的实现。其基本格式为G1=cannon(sys,type)同时，状态方程的规范实现函数cannon还具有可以返回状态变换阵的形式，即[Gl,T]=cannon(sys.type)式中，sys表示原系统状态方程模型，字串type确定规范形式的类型，它可以是模态(modal)规范型(约当标准型)，也可以是伴随矩阵(companion)形式。T是状态变换阵返回变量，满足z=7x关系，其中要求sys为状态空间模型。例3-9已知系统X(A,B,C,0的系数阵为对其进行规范型变换(约当变换)，并给出变换阵。解：应用状态方程的规范实现函数cannon容易进行该变换。A=[5210：0460；0-3-50：0-3-6-1]：B=[1234]'；C=[l252]：D=0；sys=ss(A.B,C,D)；[G,T]=canon(sys,'modal')运行结果a=xlx2x3x4xl5.0000000x20-1.0000000x300-2.000000x40001.00000b=ulxl3.60714x2-20.00000x326.55760x411.79248c=xl x2x3x4yl1.00000 2.00000 -2.75412 0.74200d=ulyl 0Continuous-timesystem.T=1.00000.60710.464300-3.0000-6.00001.00000-3.3197-6.639400-2.3585-2.35850对于系统的其他规范型变换，如可控规范型和可观规范型变换，将放到后面章节中阐述。.系统的均衡实现考查系统£(A,B,C,0),其系统的系数阵分别为A=10 ° , B=10 ,C=[10s 10s], D=0_0 -25J[10s JLJ在该系统中，系数阵的各个元素的值相差极为悬殊，显然对这类问题直接进行求解必然会在数值运算过程中由于舍入的处理而带来严重的误差，一般将这类系统称为不均衡系统。在MATLAB中，为用户提供了可以进行系统均衡变换的函数，从而使不均衡系统经过变换，使其系数变为相对较为均衡，以消除计算中舍入对其造成的严重:误差影响。其均衡变换函数的具体格式为[Ab,Bb,Cb,G,T]=balreal(A,B,C)其中T为均衡变换阵，G为均衡系统的Gram阵，并且满足下述变换关系Ab=T'AT,Bb=T'B,Cb=CT,Db=D同时系统的状态变量也满足Xb=TX关系。以前面给出的不均衡系统作为示例，利用均衡变换函数实现系统的均衡变换，即[Ab,Bb,Cb,G,T]=balreal(A,B,C)运行结果Ab=

-14.013745745757560-6.64048427128458-6.64048427128458-14.013745745757560-6.64048427128458-6.64048427128458-20.98625424242441Bb=1.373097921795070.33852931507043Cb=0.33852931507043-0.000000000517288.558136184327520.33852931507043-0.000000000517288.558136184327520.067269591423080.00273040857692T=1.0e+004*O.OOOOOOOO8555.17284303362319由上述运行结果可以清楚地看出，经过均衡变换后，系统的系数阵实现了预期的系数均衡处理，从而大大降低了数值计算时的舍入误差。应该注意的是，只有稳定的系统才能进行均衡变换。.系统的降阶实现在控制系统的研究中，模型的降阶技术是简化系统分析的重要手段，其降阶实质就是山相对低阶的模型近似成一个高阶原系统，从而使高阶模型可以按照低阶的仿真与设计方法加以进行。在MATLAB中，为用户提供了实现系统降阶处理的专用函数,如modred函数。其基本格式为RSYS=modred(sys,ELIM)RSYS=modred(sys,ELIM,，mdc')RSYS=modred(sys,ELIM,del')其中，ELIM为待消去的状态；hdd表示在降阶中保证增益的匹配；deF表示在降阶中不能保证增益的匹配。例3/0已知系统的传递函数为G(.)=--J ' 54+20s3+136,V2+380s+343应用modred函数进行降阶处理，保留前两个状态，降为二阶系统。解：先构造modred所需要的函数，再进行降阶处理。num=180；den=[120136380343][a,b,c,d]=tf2ss(num,den)sys=ss(a,b»c,d)sysm=modred(sys,3：4,'del')执行上述语句得到系统降阶后的结果为a=xlx2xl-20.00000-136.00000x21.000000b=ulxl1.00000x20c=xlx2yi00d=ulyl 0显然在直接利用modred函数进行系统降阶处理时具有一定的盲目性，为此往往将balreal函数与modred函数相结合加以使用。由balreal函数先进行均衡变换，依据Gram阵确定对系统影响较小的状态，再应用modrea函数求出降阶后的系统。例3-11已知系统Z(4,B,C,O)的系数阵为--3 1 0-11 「1-1.52 1-4J [0在尽可能保持系统基本特性的情况下进行降阶处理。解：对系统进行均衡变换。A=[-310-1；-0.5-11-1；-1.51-20；-1.521-4]B=[l000]'C=[l0-10][Ab,Bb,Cb,G,T]=balreal(A,B,C)得到均衡变换后的系统模型，以及变换阵T、Gram阵。Ab=1.28290.4033-0.19000.03590.4033-1.77481.5444-0.3144-0.19001.5444-5.92012.8156-0.03590.3144-2.8156-1.0223Bb=0.9849-0.15770.07300.0138Cb=0.9849-0.1577 0.0730 -0.0138G=0.38710.00700.00050.0001T=0.6403-1.69381.8373-2.3068-0.1591-1.0219-0.0109-0.2640-0.3446-1.53611.7643-2.2930-0.3092-1.31701.5967-1.4304由Gram阵可以看出，变换后的模型中，第3、第4状态变量对系统的作用较小，因此可以利用modred函数进行降阶处理，保留第1、第2状态变量。键入下述命令sys=ss(Ab,Bb,Cb,Db)sysr=modred(sys,3：4,'mdc')得到降阶模型为a=xl x2xl -1.27804 0.36341x2 0.36341 -1.44664b=ulxl 0.98303x2 -0.14235c=xl x2yl0.98303 -0.14235d=ulyl0.00071在MATLAB中还给出最小实现函数minreal,它的基本格式为[Am,Bm»Cm,Dm]=minreal(A,B,C,D)[numm,denmj=minreal(num,den)该函数消去了不必要的状态，从而得到系统的最小实现。有关它的具体应用可参见相关帮助文件，在此不再详述。3.1.4控制系统模型的建立与典型连接1.基本系统模型的建立MATLAB为用户提供了一些基本系统或模型建立的函数，如二阶系统、随机n阶系统的模型建立等。下面分别介绍。1)二阶系统的生成在控制系统中，二阶系统占有相当的比例，即或是高阶系统往往也需要对其进行简化降阶处理，然后再对其进行分析研究。因此，研究与讨论二阶系统具有重要的意义。MATLAB提供了二阶系统生成函数。rd2,其基本格式为[A,B,C,D]=ord2(Wn,z)式中Wn为自然角频率，z为阻尼因子，返回变量A,B,C,D描述了连续二阶系统。应用该函数就可以生成所期望的以状态方程形式描述的二阶系统，上述自然角频率与阻尼因子的意义表示为s?+2勒+Wn表示0“，z表示g。同时MATLAB也提供了生成以传递函数形式描述的二阶系统的ord2函数，其格式为[num.den]=ord2(Wn,z)返回变量为传递函数的分子、分母多项式系数向量。例3-12分别生成以状态方程形式和传递函数形式描述的二阶系统，要求J=0.8,(Dn=2.2rad/s。解：直接应用。rd2函数进行。Wn=2.2；z=0.8；[A,B,C,D]=ord2(Wn,z)得A=0 1.0000-4.8400 -3.5200B=1 0D=0求取传递函数形式如下：Wn=2.2；z=0.8；[num,den]=ord2(Wn,z)；G=tf(num»den)结果为Transferfunction：sA2+3.52s+4.842)随机n阶系统的模型建立MATLAB为用户提供了建立随机n阶系统模型的函数，其基本格式为[num,den]=rmodel(n)[num,den]=rmodel(n,p)[A,B,C,D]=rmodei(n)[A,B,C,D]=rmodel(n»p,m)[num,den]=drmodel(n)[num,den]=drmodel(n,p)[A,B,C,D]=drmodel(n)[A,B,C,D]=drmodel(n,p,m)其中[num,den]=rmodel(n)可以随机牛.成n阶稳定传递函数模型。[num,den]=rmodel(n,p)可以随机生成单入p出的n阶稳定传递函数模型。

[A,B,C,D]=rmodel(n)可以随机生成n阶稳定SISO状态方程模型。[A,B,C,D]=rmodel(n,p,m)可以随机生成n阶稳定p出m入状态空间模型。drmodel(n)函数将生成相应的离散模型。如果欲生成一随机三阶二输入单输出状态空间表达式，则可应用[A,B,C,D]=rmodel(n,p,m)o运行IA,B,C,D]=rmodel(3,1,2)得到A=-0.9688-4.09320.69414.0638-1.0010-0.4993-0.84930.0983-2.0227B=000.71430.85801.62361.2540C=-1.593700.5711D=00.69003)系统模型的重构(1).子系统选取与删除MATLAB提供了一个从大系统中选择子系统的函数ssselect,其基本格式为[Ae,Be,Ce,De]=ssselect(A,B,C；D,Inputs,outputs)利用指定的输入、输出向量建立状态方程子系统。向量inputs指定系统输入，向量outputs指定系统输出。inputs,outputs,states)[Ae,Be,Ce,De]=ssselect(A,B.inputs,outputs,states)利用指定的输入、输出及状态向量构建子系统。例3-13已知系统2(4,B,C,D)的系数阵为A-7342例3-13已知系统2(4,B,C,D)的系数阵为A-73425160253722748496123412584734262>C=22461223542显然该系统为四阶三输入、二输出系统，现利用ssselect函数在原系统的基础上构造新系统，保留1、3输入信号，1、2输出以及1、2、4状态变量。解：根据题意要求可以编写相应的程序为A=[l534；71678；30224；422579]；B=[6413；12254；3826；6472]；C=[1593；2352]；D=[224612；23542]；；inputs=[13]；outputs=[l2]；states=[124]；[Al,Bl,Cl,DI]=ssselect(A,B,C,D,inputs*outputs,states)运行上述语句得到新构建的系统模型为Al=TOC\o"1-5"\h\z1 5 47 16 842259Bl=6 312 46 2Cl=5 33 2Dl=221223 2ssselect函数除了可应用于连续函数，也可应用于离散系统。MATLAB还提供了与ssselect函数进行相反操作的一类函数，它们是[Ar,Br,Cr»Dr]=ssdelete(A»B,C»D»inputs*outputs)[Ar,Br,Cr,Dr]=ssdelete(A,B,C,D,inputs,outputs,states)利用该函数可从状态空间系统(A,B,C,0)中删除分别由inputsxOutputs和states指定的输入、输出与状态。(2).状态的增广在对系统进行分析研究时，往往需要对状态在系统(输出)中加以增广。如对系统进行全状态反馈研究时，考虑到U=Kx,因此一般就需要在输出方程中增广状态。所以状态增广有时具有十分重要的实际用途，在MATLAB中特别为用户提供了一种

状态增广函数augstate。其基本格式为asys=augstate(sys)augstate函数将状态增广到状态空间系统中，产生一个新的状态空间系统，其输入与状态和原系统相同，但其系统输出增加了全部的状态，构造成新系统为X=AX+BUUU例3-14已知系统状态方程为10024=80010024=800020B=0,C=[0.5-0.43750.4375],D=10将所有的状态增广到系统输出中。解：应用augstate函数进行输出状态增广。A=[10-42；800；020]；B=[200]'；C=[0.5-0.43750.4375]；D=l；SYS=ss(A,B,C,D)ASYS=augstate(SYS)执行上述语句得到新构建的系统为Continuoustimesystem.xlx2x3xl10.00000-4.000002.00000x28.0000000x302.000000b=ulxl2.00000x20x30c=xlx2x3yi0.50000-0.437500.43750y21.0000000y301.000000y4d=0 0 1.00000ulyiy2y3y41.00000000可见在新系统中，状态方程的系数阵没有变化，而输出方程中的C与D阵都出现了增广，满足了将全部状态添加至系统输出的要求。2.系统组合与连接所谓系统组合，就是将两个或多个子系统按一定方式加以连接形成新的系统。一般这种连接组合方式主要有串联、并联、反馈等形式。MATLAB提供了进行这类组合连接的相关函数。.模型串联在MATLAB中提供了模型串联连接函数series,该函数用于两个线性模型的串联，其基本格式为sys=series(sys1,sys2)上述串联函数实现了sysl和sys2相串联形成新系统sys。运行后形成sys=sysl*sys2如果相串联的两个系统(环节)sysl(s),sys2(s)的状态方程系数阵分别为(BPC”D()和［A2，B2,C2,D2J,则串联后整个系统的系数阵将变为B=,C=[DjCjC2],D=DjD2B=对于多入多出系统，上述串联函数的形式为sys=series(sysl>sys2,ouputsl,inputs?)应用上述函数将实现由outputsl指定的sysl输出端连接到由inputs2指定的sys2输入端，其中inputs2与outputsl分别为sys和sysl的输入、输出向量。.模型并联在MATLAB中提供了并联连接函数parallel,其基本格式为

sys=parallel（sysl»sys2»inLin2,outl»out2）该函数将两个线性系统以并联方式加以连接，ini与in2指定了相连接的输入端；outl与out2指定了进行信号相加的输出端。最终得到以vl、u>v2、作为输入zl、y、z2作为输出的系统sysoinl>in2分别为sysRsys2的输入向量，outl>out2分别为sysksys2的输出向量。当sysl、sys2为SISO系统时，该函数简化为sys=parallel（sysLsys2）这是标准并联连接的形式。sysl与sys2这两个系统在共同的输入信号作用下，将产生两个输出信号，而并联系统的输出就是这两个系统输出之和。若用传递函数对系统加以描述，系统总的传递函数G（5）=G|（5）+G2（5）o若考虑(\num\s)/\GO(\num\s)/\GO可，G«)=num2[s}de〃2(s)则系统总传递函数为G()nutnx(s\len2(s)+num1(5卜町(5)denx(s)den2[s}例3・15已知两个线性系统G-)=共应用parallel函数进行系统并联连接。解：程序如下。numl=[124]；denl=[152]；num2=[16]；den2=[171]；sysl=tf(numl,deni)；sys2=tf(num2»den2)；sys=parallel(sysl,sys2)运行上述语句得到连接后的总系统传递函数为Transferfunction：13sA3+99s八2+72s+16s八4+12s人3+38s八2+19s+2如果系统用状态方程的形式给出，则并联后的系统模型为0A2x2

Yy=[gC2]''+(D1+D2)w若己知两个线性系统,应用parallel函数即可方便地实现系统的连接Al=[12；34]；A2=[24；35]；Bl=[26J'；B2=[25]'；Cl=[l5]；C2=[23]；DI=4；D2=l；sysl=ss(Al,Bl,Cl,DI)；sys2=ss(A2,B2,C2,D2)；sys=parallel(sysLsys2)运行结果a=xlx2x3x4xl1.000002.0000000x23.000004.0000000x3002.000004.00000x4003.000005.00000b=ulxl2.00000x26.00000x32.00000x45.00000c=xlx2x3x4yi1.000005.000002.00000 3.00000d=ulyi5.00000.反馈连接反馈系统是系统中最为堂要与常见的一类系统，在MATLAB中为使用者提供了一种方便的构造反馈系统的函数feedback,其一般格式为sys=feedback(sysl,sys2»sign)执行该语句将实现两个系统的反馈连接，sign缺省时即为负反馈，sign=l时为正反馈。如果由sysl与sys2表示的前向系统和反馈系统用传递函数描述，则反馈系统的传递函数为1士s*i(s)sysi(s)若由sysl与sys2表示的前向系统和反馈系统分别以状态方程的形式给出(41,B\,IC1,DI),(A2,52,C2,02),则反馈系统的状态方程系数阵可表示为「4一8|。2。1卬-BtC2W]FBtW一一[B1ClQ Az-%。]。2Hd -[名打卬/C=[<2C,-D,c2w],D=DtW式中w=(/+£>]£>2)t，Q=I-DlWD1对于MIMO系统，可以建立更加复杂的反馈系统，其格式为sys=feedback(sys1,sys2,feedin,feedout,sign)式中的feedin与feedout均为向量，feedin为sysl的输入向量，指定哪些sysl的输入与反馈环相连接；与此类似的是，feedout指定sysl的哪些输出端用于反馈。最终实现的反馈系统与sysl具有相同的输入、输出端。例3-16已知前向环节与反馈环节的状态空间表达式的系数阵分别为-20 1 「「…1人A?=1 0,斗=0,。2=[。1],=。试将前向环节的输入1和输出2与反馈环节构成负反馈系统。解：显然由系统可知，环节1为一个多入多出系统，根据要求可以应用反馈函数feedback加以实现。其相应程序如下。Al=[10；01]；Bl=[l1；01]；Cl=[l3；20]；Dl=[10；25]；A2=[-20：10]；B2=[l0]；C2=[01]；D2=0；sysl=ss(Al,Bl,Cl,DI)；sys2=(A2»B2,C2,D2)；feedin=1：feedout=2；

sign=-l；sys=feedback(sysLsys2,feedin,feedout,sign)运行结果a=xlx2x3x4xl1.0000000-1.00000x201.0000000x32.000000-2.00000-2.00000x4001.000000b=ulu2xl1.000001.00000x201.00000x32.000005.00000x400c=xlx2x3x4yi1.000003.000000-1.00000y22.0000000-2.00000d=ulu2yi1.000000y22.000005.00000由此可以看出，经过构建的新反馈系统，其输入输出端数量没有改变,而状态发生了变化。.系统的复杂连接在系统分析与研究的实际问题中，往往会遇到一些十分复杂的控制系统，采用上面介绍的方法有时会感到难于处理，需要更加有力的手段和方法。MATLAB为进行这种工作提供了必要的选择，如可利用connect函数以及nblocks—blkbuild所构造的程序模块，实现山众多环节组合连接成复杂系统。这部分内容将在下面章节中介绍。3.2连续线性系统的计算机辅助分析与设计在经典控制理论中，频率特性及根轨迹分析方法是分析和设计连续系统的基本方法。在分析自动控制系统的工作性能时，最直观的方法是求出它的时域响应特性。但是，对于高阶控制系统的时域特性，很难用分析方法确定，而频域中的•些图解法则可以比较方便地用于控制系统的分析和设计。用频率特性方法对控制系统作出了分析和设计之后，再根据时域和频域之间的关系就能确定系统的时域响应特性。因此，用频率特性方法分析和设计控制系统有以下优点：①不必求解闭环系统的特征方程式，只要用一些较为简单的图解法就可以研究它的绝对稳定性和相对稳定性；②系统的频率特性可以用实验方法测出，这对于某些很难确定动态方程的系统或环节是非常有用的；③用频率法设计系统，对不希望的噪声可忽略不计，在规定的程度之内，频率法还可以用于某些非线性系统。在控制系统分析中，为了避开直接求解高阶特征方程式的根时遇到的困难，在实践中提出了一种图解求根法，即根轨迹法。所谓根轨迹是指当系统的某一个（或几个）参数从-8变化到+00时，闭环特征方程的根在根平面上描绘的一些曲线，应用这些曲线，可以根据某个参数确定相应的特征根。在根轨迹法中，一般取系统的开环传递系数K（又称开环系统的放大倍数）作为可变参数。由于根轨迹是以K为可变参数，根据开环系统的零极点画出来的，因而它能反映出开环系统零极点与闭环系统的极点（特征根）之间的关系。利用根轨迹可以分析系统参数和结构一定的闭环系统的时域响应特性，以及参数变化对时域响应特性的影响，而且还可以根据对时域响应特性的要求确定可变参数及调整开环系统的零极点位置，并改变它们的个数。也就是说，根轨迹法可用于解决线性系统的分析和综合问题。根轨迹法和频率特性法都具有直观的特点。由于根轨迹和闭环系统的动态响应有着直接联系，所以只要对根轨迹进行观察，用不着进行复杂计算就可以看出动态响应的主要特征。下面首先介绍用根轨迹法及频率特性法分析系统的简单原理，然后介绍利用MATLAB绘制很轨迹、奈奎斯特图和伯德图，并求系统稳定裕量的方法，最后利用这些工具对系统进行校正设计。3.3线性系统根轨迹的绘制极轨迹法是分析和设计线性定常控制系统的图解方法，使用十分简便。特别是在多回路系统研究中，应用根轨迹比用其他方法更方便，因此在工程实践中获得了

广泛应用。本节主要介绍根轨迹的基本概念，并从闭环零极点与开环零极点之间的关系推导出根轨迹方程，然后将向量形式的根轨迹方程转化为常用的相角条件和模值条件形式，利用MATLAB提供的根轨迹函数绘制跟轨迹，最后介绍设计工具。3.1根轨迹方程 G&) H" G&) H") 图3-1一般控制系统(3-8)①($)二——软)(、 (3-7)''1+G(s)”(s)为了用图解法确定所有闭环极点，令闭环传递函数表达式(3-7)的分母为零，得闭环系统特征方程l+G(s)〃(s)=O当系统有m个开环零点和n个开环极点时，(3-8)式等价为IUm)K旦 =-1 (3-9)在(S-P,)i=l式中，芍为已知的开环零点，p,为已知的开环极点，K可从零变为无穷。把(3-9)式称为根轨迹方程。根据(3-9)式，可以画出当K从零变到无穷时系统的连续根轨迹。应当指出，只要闭环特征方程能化成(3-9)式的形式，都可以绘制根轨迹，其中处于变动地位的实参数不限定是根轨迹增益K,也可以是系统其他变化参数。但是，用(3-9)式表达的开环零点和开环极点，在S平面上的位置必须是确定的，否则无法绘制根轨迹。此外，如果需要绘制一个以上参数变化时的根轨迹图，那么画出的不再是简单的根轨迹，而是根轨迹簇。根轨迹方程实质上是一个向量方程，直接使用很不方便。考虑到-1=1/供+1)" (*=0,±1,±2,…)因此，根轨迹方程(43)可用如下两个方程描述s-Zj-»s-Pi= +(m=0,±1,±2,…) (3-10)j=l /=1和flilK=T (ADn\s-zj\j=l(3—10)式和(3-11)式是根轨迹上点应该同时满足的两个条件，前者称为相角条件，后者叫做模值条件。根据这两个条件，可以完全确定s平面上的极轨迹和根轨迹上对应的K值。应当指出，相角条件是确定s平面上根轨迹的充分必要条件。这就是说，绘制根轨迹时，只需要使用相角条件；而当需要确定根轨迹上各点的K值时，才使用模值条件。3.2用MATLAB绘制极轨迹MATLAB5.3的控制系统工具箱给出了一系列关于根轨迹的函数，如表3-1所示。使用这些函数可以方便地绘制根轨迹图，以及对系统进行分析和设计。表3—1根轨迹函数函数功能函数功能pzamp绘制零极点图rlocus计算并绘制Evans根轨迹rtool根轨迹设计GUI工具sgrid在根轨迹或零极点图中添加S平面网格rlocfind从根轨迹图中选择反馈增益zgrid在根轨迹或零极点图中添加Z平面网格.计算并绘制极轨迹函数rlocus计算SISO开环模型的Evans根轨迹，根轨迹以反馈增益的函数形式给出了闭环极点的轨道(假定为负反馈)。该函数的一般用法为rlocus(G)其功能是计算并绘制开环SISO模型G的根轨迹,G可以是如图4-2所示的三种模型中的任何一种(其中K为常数)。对于图(a)模型,；对于图(b)模型，G=FGp；对于图(c)模型，G=CGp。由于系统闭环传递函数分别为"=Gp,TF=―羯一，1+KGp 1+KFGpTF=―&—.因此若G能表达为1+KCGp喇嘴则 O(s)+KN(s)=0方程的根就是系统闭环极点。rlocus函数可以自适应选择一组正增益K,产生平滑的曲线，用户也可以自己指定向量K,格式为rlocus(G,K)函数rlocus还可以带输出参数，格式为[r,K]=rlocus(G)或r=rlocus(G,K)这时，返回自适应选择的增益K和对于这些增益的复数根的位置八矩阵r的列数和增益K的长度相同，它的第m列的元素是对于增益K(m)的闭环系统的根。.选择反馈增益函数rlocfind返回在根轨迹中和一些极点相关的反馈增益。用法如下。[K,poles]=rlocfind(G)从函数rlocus生成的模型G的根轨迹图中选择反馈增益。运行rlocfind函数后，在根轨迹图上出现十字形光标，用来选择极点位置。与选中的点相关联的根轨迹增益返回到K中，列向量poles包含了该增益的闭环极点。使用该函数时，模型G的根轨迹留必须在当前图形窗口。[K,Poles]=rlocfind(Cip)用向量产指定期望的极点位置，并计算这些位置上

的根轨迹增益。向量K的第，"项是根据极点位置p(⑼计算的增益，矩阵poles的第m列是相应的闭环极点。例3-17至绘制如下系统的根轨迹图(0.055+0.045"飞2一1.85+0.9卜+5s+6)解：下面的程序可以计算并绘制根轨迹图，图3-3所示，然后用rlocfind函数在图中选择极点位置(十字光标如图所示)，得到反馈增益的值。G=tf([0,05,0,045],conv([l,-1.8,0.9][1,5,6]))rlocus(G) 图3-3根轨迹图和交互选择反馈增益K=rlocfind(G) %执行后出现下面一条信息，等待选择权点位置selectapointInthegraphicswindowselected_point= %选择极点位置后得到位置和增益值-0.8894-0.0292iK=4.9655e+003.3.3根轨迹设计工具控制系统工具箱中提供了控制系统设计的GUI工具——根轨迹设计GUKRootLocusDesignGUI)«该工具用根轨迹法交互式设计SISO补偿器。调用根轨迹设计工具的方法是在MATLAB命令窗口输入命令rltool,其窗口如图3-4所示。.窗口功能注意在窗口右上角的系统框图，图中的模块“p”为被控对象的模型，模块“K”就是设计补偿器模型，窗口左上角显示着设计过程中补偿器模型(ZPK模型).设计的最终目的就是为被控对象“p”设计出一个适当的补偿器“K”，使得整个系统(右上角的系统框图)的性能满足要求。根轨迹设计GUI中间的窗口是进行设计的主窗口，它显示整个系统的根轨迹留(由于还没有设计补偿器，而且“F”和“H”都等于1,因而图4-4中的根轨迹图就是被控对象的根轨迹)。主窗口左上角有四个按钮，分别对应不同的编辑功能。

①缺省时按钮固被按下，这时在主窗口中可以用鼠标交互式地移动闭环极点（根轨迹上的红日方证）和设计出的补偿器的零极点。②按钮圆或匐被按下时，可以在主窗口中添加补偿器的极点或零点。③按钮囚被按下时，可以在主窗口中橡除补偿器的极点或零点。主窗口的左下方是坐标轴设置按钮，它们储存当前的坐标范围;②载入存储坐标范围;把坐标轴设置为正方形;设置坐标轴的纵横比相等。沿X方向放大;沿丫方向放大;观察全部图形。图3-4根轨迹设计GUI主窗口右下方是放大工具按钮，它的作用分别为：③按钮囚被按下时，可以在主窗口中橡除补偿器的极点或零点。主窗口的左下方是坐标轴设置按钮，它们储存当前的坐标范围;②载入存储坐标范围;把坐标轴设置为正方形;设置坐标轴的纵横比相等。沿X方向放大;沿丫方向放大;观察全部图形。图3-4根轨迹设计GUI主窗口右下方是放大工具按钮，它的作用分别为：沿x方向和丫方向都放大的作用分别是:「Stepfleto「8。”厂型㈤厂N卧点I，用广选杼需要显示的系统时赢扁和频率晌应分析。选中J0P的•个或几个靛后,将调用LTIViewer来髭小①镐④响应的曲线，随着设计过程中系统模型的变化，响应也随之变化，这样便于在设计过程中随时观察系统的特性。.输入输出模型根轨迹设计GUI窗口菜单File下的ImportModel用来输入LTI设计模型，其窗口如图3-5所示。窗口左上角显示的是系统的结构框图，它允许两种反馈结构，另种如图3d右上角所示。在这里输入的模型包括被控对象P、前置滤波器F、传感器动态特性H。缺省情况下，P、F、H的值是1。图3-5载入LTI设计模型在窗口左下角指定输入模型的位置，包括工作空间、MAT文件和Simulink框图。窗口中间显示了指定位置（如工作空间）中的所有模型名称。如果要使被控对象P的模型为GL则首先选定模型G1,然后在窗口右下方P前面的箭头上点下，被控对象P的模型就成为G1。同样可以为H和F指定相应的模型。系统设置完成后，按下“0K”按钮返回主窗口，这时主窗口中间会显示出输入的LTI模型系统的根轨迹图（图3—4）,其中红色的小方块表示闭环极点。类似地，菜单File下的ImportCompensator选项用来输入要进行设计的补偿器K的模型。模型输入后，主窗口中显示出整个系统的根轨迹图，在图中可以对补偿器K的零极点进行增删和编辑，但不能编辑原系统的零极点，只能沿轨迹移动闭环极点。模型的输出用菜单File下的Export选项，其窗口如图3—6所示。模型可以输出到磁盘MAT文件，也可以输出到工作空间。窗口左侧显示的几项内容的含义分别为：DM指设计的模型；0L指开环模型；CL指闭环模型；K指补偿器模型。

.设置网格和边界线利用Tools菜单下的AddGrid/Boundary选项可以在根轨迹图中添加网格并设置边界线。其设置窗口如图3—7所示。例如，需要显示阻尼比和固有频率的网格，则只要选中该项前的单选按钮，然后选择下面的“GridOn?”前的复选框即可。同时，设计的系统需要满足条件：稳定时间小于0.04s,阻尼比小于0.4,那么只要在该窗口的Setting 图3—6输出模型Time后输入0.04,在DampingRatio后输入0.4即可。度置的内容见图3—7,设友CaronjFsJNwWFfwwcyPtMcBvwthociPOJPSdUlry：r：e-代闲四手［I飞f：..「■”友CaronjFsJNwWFfwwcyPtMcBvwthociPOJPSdUlry：r：e-代闲四手［I飞f：..「■” -：图3—7网格和边界线设置窗图3—8添加了网格和边界的根轨迹图.转换为Simulink框目主窗口File菜单下的DrawSimulinkDiagram选项能够把LT1系统模型转换为相应的Simulink仿真框图，这样便于在Simulink中对,系统作进〜步的仿真分析。

3.4线性系统频率特性与奈奎斯特图频域分析法是利用频率特性研究线性系统的一种经典方法。绘制奈奎斯特(Nyquist)图是频率法中利用控制系统的开环幅相频率特性，判断其闭环系统稳定性的重要方法之一，它是由著名的奈奎斯特判据而来的。开环系统的幅相频率特性容易计算，且可通过实验求得，因而奈氏判据使用方便，同时物理意义明确。这个判据确定了开环系统的频率持性与闭环系统动态响应之间的联系，它不仅能判断闭环系统的稳定性，而且可利用它找到改善闭环系统动态响应的方法，因而它是构成经典设计方法的重要基础。3.4.1奈奎斯特稳定判据奈奎斯特判指定义如下：①开环系统稳定时，即开环系统没有极点在右半根平面，如果相应于频率3从-8f+8变化时，开环系统频率特性G(/0)”(j0)曲线不环绕(-1,jO)点，那么闭环系统就是稳定的，否则是不稳定的；②开环系统不稳定时.即开环系统有P个极点在右半根平面时，如果相应于频率。从-8f+8变化时.开环系统频率特性G(jco)H(jco)曲线反时针方向环绕(-1,jO)点的次数N等于位于右半根平面内的开环系统的极点数P,那么闭环系统就是稳定的，否则是不稳定的。奈奎斯特图就是根据奈氏判据,在开环系统频率特性的平面上绘制开环系统幅相频率特性的极坐标图。例3-18设有-控制系统的开环系统传递函数G(s)〃(G(s)〃(s)=(1+7>)(1+7^)试判别其闭环系统的稳定性。解：开环系统的频率特性为(1+(期(1+72%)

当3=0时,G(jo)/7(，＜y)=k/o°当3f+8时，G(J(o)H(Ja))=K/180。上面这两个特殊点确定了奈奎斯特曲线的变化趋势，再计算对应几个。的G(的值，便能绘出如图3-9所示的开环系统幅相频率特性曲线图，即奈氏图。因为例中给出的开环传递函数没有极点在右半根平面，而且不论K怎样变化，G(j0)”(Oy)曲线都不会包围(-1J0)点，因此闭环系统总是稳定的。下面通过例题准确绘制奈奎斯特曲线。G(s)”(s)=K(l+7>)G(s)”(s)=K(l+7>)s(\+T2s)(\+T3s)(\+T4s)(3-12)求其模值和相角的表达式。解：(3-6)式可以写成 K(l+.」o) e-Ttj(0=A(a))e讷⑻+工力。)。+73，。)(1+T4j(o)这里A3)=K|l+7Jo|

1M中这里A3)=。(0)=--一arctan(T269)-arctan(7^69)-arctan(7^69)+arctan(Tg)-T5co利用上面的表达式编制计算程序，指定具体数据，就可以绘制出奈奎斯特图了。3.4.2利用奈奎斯特图来分析系统的稳定裕量绘制控制系统的奈奎斯特图，不仅可以确定系统

的稳定性，还可以确定系统的稳定裕量。奈氏曲线离G(ja))H(jco)平面上的(J,jO)点越远，说明稳定裕星越大。稳定裕量是以G(j0)”(jo)轨迹上两个特殊点的位置来度量的。所谓特殊点是指如图4-10所示奈奎。叫,拜公1—S斯特图的A点和B点。A点规定为的轨迹和单位圆的交点(交点频率设为3C)。显然，在A点，G("y)"("y)的幅值为A(0c)=l，而相角为 图3-10由奈氏图示稳定裕量-以0°)(相角规定为顺时针为负，逆时针为正)。相位裕量定义为向量OA与负实轴间的夹角，并用，表示。根据图3-10,可由下式确定：/=180—(p〈a)c)(3-13)对最小相位系统来说；(3-13)①r>o,表示相位裕量有正值，说明系统趋向稳定；②r<0,表示相位裕量为负值，说明系统趋向不稳定；③r=o,表示A点在负实轴上，系统为临界稳定，这时奈氏曲线穿越(-1J0)点。通常y的值在30°-180°之间。B点规定为叼轨迹和负实轴的交点，交点频率设为3g。增益裕量(又称为幅值裕量)定义为G(j3)H(j3)轨迹在B点的幅值43m的倒数，并Kg表示，即(3-14)其中，Kg>l,表示系统趋向于稳定；Kg<A，表不系统趋向于不稳定；3=201gA(%)Kg=l,表示奈氏曲线穿越(-l,j0)3=201gA(%)3.4.3利用MATLAB绘制奈奎斯特(Nyquist)曲线在MATLAB中给出如下绘制频率特性曲线的相关函数，见表3—2。利用这些函数可以绘制出nyquist,bode,nichols等曲线，并进行增益裕量和相位格量的分析。表3-2频域分析函数函数功能bode计算并绘制bode图evalfr计算模型在某个指定复频率点上的频率响应(不包括FRD模型)freqresp计算在一组实频率点上的频率响应margin计算增益裕度和相位裕度ngrid在nichols图上加nichols图nichols计算并绘制nichols图nyquist计算并绘制nyquist图sigma计算并绘制奇异值图控制系统工具箱中提供了一个MATLAB函数nyquist,该函数可以用来直接求解来奎斯传阵列，或绘制出来奎斯特国。该函数可以由下面的格式调用[rx,ry]=nyquist(C,w)其中G为系统的LTI对象模型。变量给w出要分析的频率点构成的频率向量，返回的向量分别为系统奈奎斯传阵列的实部和虚部。如果只给出一个返回变量，则返回的变量为复数矩阵，其实部和虚部可以用来给制系统的奈奎斯特图。应该注意，这两个返回的变量均是三维矩阵。所谓三维矩阵，实际上是针对多变量系统而设计的，第i个输出变量针对第j输入变量的频域响应实部和虚部分别由rx(l,j,:)和ry(l,j,:)给出。对于单变量系统，实部和虚部分别由rx(U,:)和给出,更简单地，可以山rx=rx(:)和ry=ry(:)来提取实部和虚部向量。该函数还可以由下面各种其他形式来调用[rx,ry,w]=nyquist(G),或nyquist(G)与其他控制系统工具箱中的函数类似，前一种格式会自动选择频率向量，而后一种格式能自动绘制出系统的条奎斯特图。例3-20系统开环传递函数G(s)=1000G(s)=1000(52+35+2)(5+5)绘制系统的条奎斯特图，并讨论其稳定性。解：用下面的MATLAB命令直接绘制出系统的条奎斯特图，如图3—11(a)所示。G=tf(1000,conv([l,3,2],[1,5]))；nyqulst(G)；axis(,square")该图中(-1,jO)点附近奈奎斯特图的情况并不很清楚，可以用CtrlLAB程序对得

出的奈奎斯特图进行局部放大(选择图形窗口中的Zoomin菜单项，或用快捷按钮侬),则可以得出如图3—11(b)所示的局部放大结果。从局部放大的图形中可以看出.来奎斯特图逆时件包围(-1,jO)点2次，而原开环系统中没有不稳定极点，从而可以得出结论，闭环系统有2个不稳定极点。这也可以由下面的MATLAB语句来进一步验证。G_close=feedback(G.1)；roots(G_clos.den{1})-12.81962.4098+8.5427i2.4098-8.5427i-30-20NyquistDiaqrams40402010010-NyquistDiagrams-30-20NyquistDiaqrams40402010010-一400 20 40 60 80100 -20-10 0 1020 30Rea]Axis Rea!Axis(a) (b)图3—11例3-20的奈奎斯特囹3.4.4幅值裕量与相位裕量控制系统工具箱中提供了margin函数来求取给定线性系统的幅值裕量和相位裕量，该函数可以由[Gm,Pm,wcg，Wcp]=margin(G)格式来调用。可以看出，幅值裕量与相位裕量可以由ETI对象G求卜H。返回的变量(Gm,Wcg)为幅值裕量的值与相应的频率，而(Pm,Wcg)则为系统的相位裕量与相应的频率。若得出的裕量为无穷大，则得出的值为inf,这时相应的频率值为NaN(表示非数值，这两个量均是MATLAB保留的常数)。如果已知系统的频域响应数据，还可以由下面的格式

[Gm,Pm,Wcg,wcp]=margin(mag,phse,w)；调用此函数，其中(mag,phase,w)分别为频域响应的幅值、相位与频率向量。例3-21考虑如下系统模型,求它的幅值裕量和相位裕量，并求其闭环阶跃响应。G(s)=G(s)=3.5§3+2sl+3s+2解：由下面的MATLAB语句直接求出系统的幅值裕量与相位裕量及闭环阶跃响应。G=tf(3.5,[l,2,3,2])[Gm,Pm,Wcg,Wcp]=margin(G)G_close=feedback(G,1);step(G_close)得出闭环系统的阶跃响应曲线如图3—12(a)所示，系统的幅值裕量很接近稳定的边界点1,且相位裕量只有7.578。所以尽管闭环系统稳定，但其性能将不会太好。可以看出在闭环系统的阶跃响应中有较强的振荡，这就验证了前面得出的结论。0.2 0.20.2 0.2 o'_1~1~- 0 50 100 150 00.05OJ(M5Q.20250303504Q460.5Time(s) Time(s)(a) (b)图3-12闭环系统的阶跃响应如果系统的相位裕量y>45°,则称该系统有较好的相位格量，然而这样的数值并不很绝对。考虑一个新的系统模型G(s)=G(s)=(5+1)(52+S+9)可以由下面的MATLAB语句直接求出系统的幅值裕量与相位裕量为G=tf(100*conv([l,5],[l,5]),conv([l,l],[l,l,9]);[Gm,Pm,Wcg,WcP|=margin(G)；[Gm,Pm,Wcg,WcpJ

ans=Inf85.4356NaN100.3281G_close=feedback(G,1)；step(G_close)并可以得出系统的闭环阶跃响应如图4-12(b)所示。可以看出，该系统有无穷大的幅值裕量，且相位格量高达85.4356°,所以系统的闭环响应是较理想的3.5线性系统的伯德图分析3.5.1伯德图利用奈奎斯特图可以从开环幅相频率特性分析闭环系统的稳定性及稳定裕量，但是奈氏曲线不能明显地表示系统增益及时间常数等参数变化时对系统性能的影响。为了克服奈奎斯特图分析设计系统的缺点，将系统的开环频率特性的幅频特性曲线和相频特性曲线分开画，并对幅频特性进行对数运算(相频特性不取对数)，横坐标表示频率3,纵坐标分别表示对数幅值和相角，分别画出的对数幅频特性曲线图和相频特性曲线图称为伯德(Bode)图。为了能把一个较宽频率范围的图形紧凑地表示在一张图上，横生标采用对数分度，而表示对数幅值和相角的纵坐标仍然用线性分度.这样就形成了一种半对数坐标。画在半对数坐标L的对数幅频特性和相频特性又称对数坐标图。因此，伯德图又称对数坐标图。下面举例说明伯德图的画法。例3-22设一控制系统的开环传递函数为(3-15)(3-16)(3-17)G(s)H(s)=(3-15)(3-16)(3-17)试绘制相应的伯德图。解：其幅频和相频特性分别为,、 K4。)=/ ，—/ ,+•y/l+(T2G))0(0)=---arctan(Tg)-arctan(T26y)绘制伯德图时.对幅值4。)两边取常用对数.得lgA(0)=1gK-lg口一lgJl+(7»2-1gJl+W。)? (3-18)

这样就把幅值的乘除运算变成了加减运算，使计算简化。通常在对数前面还乘上系数20以使对数幅频的计算单位变成分贝dB表示、为了方便，以分贝为单位的对数幅频特性常用L(。，来表示即(3—18)或可改写为=201gA(⑼=201g/C-201g<y-101g[l+(Tl<y)2]-101g[l+(7'2ty)2](3-19)设(3—9)式中的K=20、(=1、T2=\,则由(3T7)式和(3T9)式可画出系统的伯德图，如图3-13所示。13.5.2伯德图的渐近线从前面的分析中可以看到，要想给制比较复杂的伯德图，需要借助于计算机才能实现。事实上，如果只想定性地研究系统的频域响应，则使用伯德图的渐近线就足够了。系统伯德图的渐近线可以由比较简单、规范的方法容易地绘制出来。这里首先考虑幅频特性渐近线的切换点的规则。当频率值达到G⑶的一个单独极点时，幅频传性的渐近线的斜率将减小20dB/十倍频；若频率值到达了一•对共版的极点时，则渐近线的斜率值减小到40db/十倍频时，如果频率值到达一个p重极点时，福频特性的渐进线的斜率减小20PdB/十倍频；若频率值到达了一对共瓢的零点时，则渐近线的低利率将增加40dB/十倍频；如果频率值到达一个p重的零点，则渐近线的斜率将增加20pdB/十倍频；如果系统含有m个微分和〃个积分环节.则斯近线的起始频率应该是20(m〃)dB/十倍频。根据上述原则可以编写一个MATLAB函数bode_asym来求取系统的伯德幅频特性渐近线的数值，该函数的调用格式为[wpos,ypos]=bode_asym(G,w)使用bode_sym函数的程序如下：function[wPos,ypos]=bode_sym(G,w)G1=zpk(G);wpos=[];slope=[];Ifnargin==1,w=freqint2(G);endzero=G1.z{1};polar