3．机电系统的计算机辅助分析与设计 3．1 机电系统的数学模型及其转换方法 机电系统计算机仿真与辅助设计是建立在机电系统数学模型基础之上的。对于各类机电系统，利用仿真手段对其进行分析与设计，首先就需要建立相应的系统数学模型，此后，就需要研究如何将系统的数学模型转变为适合于计算机进行分析计算的仿真模型，即数值算法模型。在此基础上，即可通过对数学模型的求解分析，实现对系统动静态特性的分析与设计。显然，进行上述工作的重要基础就是系统的数学模型。因此本章首先介绍系统的几种典型数学描述，然后介绍各种数学模型之间的相互转换，以及系统环节不同形式的相互连接的 MATLAB实现。 3．1．1 连续系统的数学描述 连续系统的数学模型通常可以用微分方程、传递函数、状态空间表达式三种形式对系统加以描述。下面将简单对这几类数学模型加以回顾，同时给出MATLAB对它们的表示方法。 1．系统的微分方程形式模型 一个系统的动态特性通常可用高阶微分方程加以描述，因此描述一个系统最常用的数学 模型就是微分方程的形式。假设连续系统为单入单出（简称SISO）系统，其输入与输出分别用u（t）、y（t）加以表示，则描述系统的高阶微分方程为：  （3-1） 其初始条件为：，，…，，，… 如果引入微分算子 ，则（3-1）式可以写作：  即  对上式稍加整理并令，可以得到  （3-2） 2．系统传递函数形式模型 1） 传递函数模型 对（3－1）式等号两边取拉氏变换，并假设y与u的各阶导数的初值均为零，则存在  （3-3） 式中：——输出的拉氏变换； ——输入的拉氏变换。 从而（3-1）式所描述的系统的传递函数为  （3-4） 对照（3-2）式与（3-4）式可以清楚地看出，当描述系统的微分方程的初始值为零时，用算子p所表示的式子与传递函数G（S）在形式上完全相同。 传递函数是经典控制论描述系统的数学模型之一，它表达了系统输入量和输出量之间的关系。它只与系统本身的结构、特性和参数有关，而与输入量的变化无关，传递函数是研究线性系统动态响应和性能的重要手段与方法。在MATLAB语言中，可以利用分别定义的传递函数分子、分母多项式系数向量方便地对其加以描述。例如对于（3-4）式，系统可以别定义传递函数的分子、分母多项式系数向量为 num=[ c1 c2…cn-1 cn] den=[1 a1 a2…an-1 an] 这里分子、分母多项式系数向量中的系数均按s的降幂排列，由于传递函数G(s)的最高次项系数为 1，所以分母多项式系数向量den中第一个元素为1。在 MATLAB 5.X中，可以用tf来建立传递函数的系统模型，其基本格式为 sys=tf(num,den) 例3-1已知系统传递函数为  解：可以利用MATLAB将上述系统模型表示出来，并将其建立在工作空间（Workspace）中，写成 num＝[2 9]； den＝[1 2 3 4 6]； printsys（num，den，'s'） 执行上述语句后得到 num/den= 2s+9 s ^ 4 + 3 s ^ 3 + 2 s ^ 2 + 4 s + 6 对于 MATLAB 5．x可以利用 tf直接建立系统模型，即 num=［2 9］； den＝［1 3 2 4 6]； model=tf(num,den) 执行上述语句后得 Transfer function： 2s+9 s ^ 4 + 3 s ^ 3 + 2 s ^ 2 + 4 s + 6 由上例可以看出，对于MATLAB的不同版本建立的传递函数有所不同，但结果是一致的。下面考虑一种较复杂的传递函数形式。 例3-2 已知系统传递函数为  解：可以利用MATLAB建立其相应的传递函数系统模型。 num= 7*［2 3］； den= conv（conv（conv（［1,0,0］,［3,1］）,conv（［1 2］,［1,2］））,［5,0,3,8］）； model=tf（num,den） 运行结果 Transfer function： 14s+21 15s ^ 8 + 65s ^ 7 + 89s ^ 6 + 83s ^ 5 + 152s ^ 4 + 140s ^ 3 + 32s ^ 2 在这里使用了conv函数，它是MATLAB中的标准函数，用来求取两个向量的卷积分。因此当两个多项式相乘，需两个多项式系数向量相乘时，就可利用conv函数来进行。conv函数允许多重嵌套，由上例已清楚地看到这一点。 对于连续时间系统可以用传递函数对其加以表示，而对于离散时间系统则应采用脉冲传递函数对其进行描述。脉冲传递函数一般可表示为关于z的降幂多项式分式形式，即  在 MATLAB 5．0中，对于离散系统同样可以建立相应的系统模型，其基本格式为 num=[ cm，cm-1，…，c1， c0]； den=[ an，an-1，…，a1， a0]； sys＝ tf（num，den，T） 其中T为系统采样周期。 2） 系统的零极点形式模型 系统的传递函数还可表示成另一种形式，即零极点形式。这种形式的系统传递函数比标准形式传递函数更加直观，可清楚地看到系统零极点分布情况。系统的零极点模型一般可表示为  （3-5） 其中，和分别为系统的零点和极点，、既可以是实数也可是复数，K为系统增益。MATLAB可以使用zpk函数建立零极点形式的系统模型．其基本格式为 sys= zpk（[z]，[p]，[k] ） 其中，[z]、[p]、[k]分分别为系统的零极点和增益。 3．系统的状态空间表达式 状态方程是研究系统的最为有效的系统数学描述，不论是单入单出系统还是多入多出(简称MIMO）系统，若可用一组一阶微分方程对其加以表示，在引进相应的状态变量后，则可将这一组一阶微分方程写成紧凑形式，即状态空间表达式。通常一个线性定常系统可以表示为  （3-6） 其中上述第一式由n个一阶微分方程构成，称为系统的状态方程表达式，下式由l个线性代数程组构成，称为系统的输出方程。X为n维状态向量；U为m维输入向量；Y为l维输出向量；A为n×n维的系统状态阵，由系统的参数所决定；B为n×m维的系统输入阵;C为l×n维输出阵；D为l×m维直接传输阵。 应用MATLAB可以方便地表示系统的状态方程模型，只要按照矩阵输入方式建立相应的系统系数阵即可，考虑SISO系统，容易在MATLAB工作空间中建立系统的系数阵，形式为 A=[a11 a12…a1n ；a21 a22…a2n ；…；an1 an2…ann]； B=[b1； b2；…；bn]； C=[c1 c2 …cn]； D=d； 当然，也完全可以在MATLAB工作空间中建立MIMO系统的系数阵。根据系统状态方程的系数阵，可在MATLAB中建立相应的系统模型，其基本格式为 sys=ss (A，B，C，D) 上述是 MATLAB 5.x中的格式，在 MATLAB 4.2中其格式为 printsys（A，B，C，D） 对于离散系统，其状态空间表达式可表示成 X（k十1）=AX（k）＋BU（k） Y（k）＝ CX（k）＋ DU（k） 在MATLAB 5.x中同样也可建立相应的系统模型，其格式为 sys=ss（A， B， C， D， T） 其中，T为系统采样周期。在 MATLAB 4．2中也可建立类似的系统模型。 3．1．2 系统模型的相互转换 在前一节中已经介绍了描述一个系统的数学模型主要有三种形式；微分方程模型、传递函数模型（包括零极点形式模型）、状态空间模型等。显然在这些不同形式的系统数学模型中存在着内在的联系，虽然它们外在形式不同，但它们的实质内容是等价的。人们在进行系统分析研究时，往往根据不同的要求选择不同形式的系统数学模型，因此研究不同形式的数学模型之间的转换具有重要意义。 1． 机电系统模型向传递函数或零极点增益形式的转换 1）．状态空间方程向传递函数形式的转换 系统的状态空间方程可表示为  据此，可以得到等效的系统传递函数模型。 =C(sI-A)-1B+D 显然在进行这种变换过程中，求取（sI－A）阵的逆比较困难，然而MATLAB有一系列的函数可以完成各种变换，其中就包括进行这种变换的ss2tf函数，其基本格式为 [num，den]＝ss2tf（A，B，C，D，iu） 利用该函数即可实现将状态空间方程转换为传递函数的形式，iu用于指定变换所使用的输入量。为了获得传递函数的系统形式，还可以采用下述方式进行，即 G1=ss（A，B，C，D）； G2=tf（G1） 可以证明，由给定的状态方程模型转换为传递函数形式其结果是惟一的。 例3-3 已知连续系统 Σ（A，B，C，D）的系数矩阵是 A= ，B= ， C=， D= 求取该系统相应的传递函数模型。 解： 应用MATLAB的ss2tf函数可以方便实现这种转换。 A=[2 0 0；0 4 1；0 0 4]； B=[1 0 1]’ ： C=[1 1 0]； D=0； [num den]=ss2tf（A，B，C，D，1） G=tf（num，den） 运行上述语句即可实现要求的转换，利用G=tf（num，den）将转换后的传递函数分子，分母多项式系数向量构造成传递函数系数模型，下面给出了运算结果； num= 0 1.0000 -7．0000 14．0000 den= l -10 32 -32 Transfer function： s ^ 2 – 7 s + 14 s ^ 3 – 10 s ^ 2 + 32 s - 32 对于多输入系统，应用ss2tf函数可以进行指定要求的模型转换。 例3-4 已知连续系统软Σ（A，B，C，D）的系数矩阵是 A=，B=，C=，D= 求取系统相应的传递函数模型。 解：该系统为一双输入系统，利用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=[4 6；2 4；2 2；0 2]； C=[0 2 0 2」； D=[0 0］； T=2； ［num2，den2］=ss2tf(A,B,C,D,T) G=tf（num2，den2） 运行结果 num2= 0 12.0000 32．0000 37．000 17．0000 den2= 1．0000 4．0000 6．2500 5.2500 2.2500 Transfer function： 12 s ^ 3 + 32 s ^ 2 + 37 s + 17 s ^ 4 + 4 s ^ 3 + 6.25 s ^ 2 + 5.25 s + 2.25 当T=1时，将指定第一输入信号所对应的系统转换，下面给出了相应的转换结果。 num1= 0 4.0000 14.0000 22.0000 15.0000 denl= 1.0000 4.0000 6.2500 5．2500 2．2500 Transfer function： 4 s ^ 3 + 14 s ^ 2 + 22 s + 15 s ^ 4 + 4 s ^ 3 + 6.25 s ^ 2 + 5.25 s + 2.25 2）．模型向零极点形式的转换 系统的零极点模型实际上是传递函数模型的另一种形式，也是对系统进行分析的一类常用模型，因此在MATLAB中，也提供了实现将各类系统模型转换为零极点形式模型的函数。其基本格式为 [z，p，k］=ss2zp（A， B，D，IU） [z，p，k]=tf2zp（num，den） Gzp=zpk（sys） 上述第一式是将以状态方程形式给出的模型根据指定研究的输入，转换为零极点模型形式；第二式是将传递函数形式给出的模型转换为零极点形式；第三式十分简洁，用该函数可将非零极点形式的模型转换为零极点系统模型。 例3-5 已知系统状态方程的系数阵为 A=，B=，C=，D= 试将其转换成零极点模型形式。 解：利用转换函数 [z ，P，k]=ss2zp（A，B，C，D） 即可实现所希望的转换，运行结果为 z= -1．3333 p= -2 -1 k＝ 1.5000 也可利用Gzp=zpk（sys)函数进行转换，即 sys＝ss（A， B， C， D） Gzp=zpk（sys） 运行结果 Continuous-time system． Zero/pole/gain； 2（s＋0.5） (s十2) (s+1) 显然，两种方法所得到的结果是相关的。 2． 系统模型向状态方程形式的转换 前面介绍了系统模型转换为传递函数的方法，同样也可以利用MATLAB函数实现所需要的系统模型向状态方程的转换。其基本格式为 [A，B，C，D]=tf2ss（num，den） [A，B，C，D]=zp2ss（Z，p，k） syss=ss（sys） 例 3-6已知系统传递函数为   应用MATLAB的模型转换函数将其转换为状态方程形式的模型。 解: 应用tf2ss转换函数很容易实现所要求的转换。但应注意，由于其状态变量选取的不同，转换结果是不惟一的。下面给出了应用tf2ss函数的转换程序与结果。 num＝[18,36]； den=［1 40.4 391 150］； ［A，B，C，D］=tf2ss（num，den） 运行结果为 A= -40.4000 -391.0000 -150.0000 1.0000 0 0 0 1.0000 0 B＝ 1 0 0 C＝ 0 18 36 D= 0 利用[A，B，C，D]=zp2ss(z,p,k)可以将零极点形式给出的模型转换为状态方程。对于sys=ss（sys），可以将任意LTI系统模型转换为状态方程。 例3-7 考虑例3-6，将程序进行修改。 解： 修改程序如下。 num=［18，36］； den=［1 40.4 391 150］； G=tf（num， den） Sys= ss（G） 运行结果 A= -40.4000 -391.0000 -150.0000 1.00000 0 0 0 1.0000 0 B= 1 0 0 C＝ 0 18 36 D＝ 0 3．1．3 系统状态方程的变换与实现 1．状态方程的相似变换 在前一节已经述及将系统的状态方程转换为传递函数时，其结果是惟一的；然而由于状态变量的选取不同，将系统的传递函数转换为状态方程则不能保证结果的惟一。换言之，对于一个系统的传递函数，可以存在着众多的状态方程实现。这样，对于同一个系统有着不同的状态方程的描述，因此也就存在着它们之间的相互变换，MATLAB中的ss2ss函数可以实现对系统状态方程的相似变换。其基本格式为 GT=ss2ss（G，T） 其中，G为原系统的状态方程模型，T为非奇异变换阵。 该变换实现下列相似变换 z = Tx， = AtU + BtU y = Ctz + Dtu 其中 At = TAT-1， Bt = TBT-1， Ct = TCT-1， Dt = TDT-1 下面通过一个例子来说明这种变换。 例3-8 已知系统状态方程为 A =  ，B =  ，C = ，D = 试应用ss2ss函数进行状态方程的相似变换。 解： 由于变换阵T可以任意选择，只要保证其非奇异即可。在此选择单位反对角阵作为变换阵，下面给出了变换方法与结果。 A =［10 -4 2；8 0 0；0 2 0］； B =［2 0 0］'； C =［0.5 -0.4375 0.4375］； D = 0； G1 = ss(A， B， C， D) T = fliplr(eye（3））； GT = ss2ss（G1，T） 得到结果为 a = xl x2 x3 xl 0 2.00000 0 x2 0 0 8. 00000 x3 2. 00000 -4.00000 10. 00000 b = u1 x1 0 x2 0 x3 2. 00000 c = x1 x2 x3 y1 0.43750 0.43750 0. 50000 d = ul y1 0 2．规范型状态方程的实现 在MATLAB中提供给用户一个状态方程的规范实现函数cannon，以进行LTI系统模型sys的规范状态空间表达式的实现。其基本格式为 G1 = cannon（sys,type） 同时，状态方程的规范实现函数cannon还具有可以返回状态变换阵的形式，即 [G1，T] = cannon（sys，type） 式中，sys表示原系统状态方程模型，字串type确定规范形式的类型，它可以是模态（modal）规范型（约当标准型），也可以是伴随矩阵（companion）形式。T是状态变换阵返回变量，满足z=Tx关系，其中要求sys为状态空间模型。 例3-9 已知系统Σ(A，B，C，D)的系数阵为 A=，B=，C=，D= 对其进行规范型变换（约当变换），并给出变换阵。 解：应用状态方程的规范实现函数cannon容易进行该变换。 A=[5 2 1 0；0 4 6 0；0 -3 -5 0；0 -3 -6 -1]； B=［1 2 3 4］'； C=[1 2 5 2］； D=0； sys=ss（A，B， C， D）； [G，T]＝ canon（sys，'modal'） 运行结果 a= xl x2 x3 x4 xl 5.0000 0 0 0 x2 0 -1.00000 0 0 x3 0 0 -2.00000 0 x4 0 0 0 1.00000 b= ul xl 3.60714 x2 -20.00000 x3 26.55760 x4 11.79248 c= xl x2 x3 x4 y1 l.00000 2.00000 -2.75412 0.74200 d= ul y1 0 Continuous-time system． T= 1.0000 0.6071 0.4643 0 0 -3.0000 -6.0000 1.0000 0 -3.3197 -6.6394 0 0 -2.3585 -2.3585 0 对于系统的其他规范型变换，如可控规范型和可观规范型变换，将放到后面章节中阐述。 3．系统的均衡实现 考查系统Σ(A，B，C，D)，其系统的系数阵分别为 A=，B=，C=，D= 在该系统中，系数阵的各个元素的值相差极为悬殊，显然对这类问题直接进行求解必然会在数值运算过程中由于舍入的处理而带来严重的误差，一般将这类系统称为不均衡系统。在MATLAB中, 为用户提供了可以进行系统均衡变换的函数，从而使不均衡系统经过变换，使其系数变为相对较为均衡，以消除计算中舍入对其造成的严重误差影响。其均衡变换函数的具体格式为 ［Ab，Bb，Cb，G，T]= balreal（A， B， C） 其中T为均衡变换阵，G为均衡系统的Gram阵，并且满足下述变换关系 Ab=T-1AT，Bb=T-1B， Cb=CT，Db=D 同时系统的状态变量也满足Xb= TX关系。 以前面给出的不均衡系统作为示例，利用均衡变换函数实现系统的均衡变换，即 [ Ab，Bb，Cb，G，T]=balreal（A，B，C） 运行结果 Ab= - 14.013745745757560 -6.64048427128458 -6.64048427128458 -20.98625424242441 Bb＝ 1.37309792179507 0.33852931507043 Cb＝ 1.37309792179507 0.33852931507043 G= 0.06726959142308 0.00273040857692 T= 1.0e+004 * 0.00000000855 -0.00000000051728 5.17284303362319 8.55813618432752 由上述运行结果可以清楚地看出，经过均衡变换后，系统的系数阵实现了预期的系数均衡处理，从而大大降低了数值计算时的舍入误差。应该注意的是，只有稳定的系统才能进行均衡变换。 4．系统的降阶实现 在控制系统的研究中，模型的降阶技术是简化系统分析的重要手段，其降阶实质就是由相对低阶的模型近似成一个高阶原系统，从而使高阶模型可以按照低阶的仿真与设计方法加以进行。在MATLAB中，为用户提供了实现系统降阶处理的专用函数，如modred函数。其基本格式为 RSYS=modred（sys，ELIM） RSYS=modred（sys，ELIM，'mdc'） RSYS=modred（sys， ELIM，'del'） 其中，ELIM为待消去的状态；'mdc'表示在降阶中保证增益的匹配；'del'表示在降阶中不能保证增益的匹配。 例3-10 已知系统的传递函数为  应用modred函数进行降阶处理，保留前两个状态，降为二阶系统。 解：先构造modred所需要的函数，再进行降阶处理。 num= 180； den= [1 20 136 380 343] ［a，b，c，d］= tf2ss（num，den） sys＝ss( a，b，c，d) sysm= modred（sys，3：4，'del'） 执行上述语句得到系统降阶后的结果为 a= xl x2 xl -20.00000 -136.00000 x2 1.00000 0 b= ul xl l.00000 x2 0 c= x1 x2 y1 0 0 d= u1 y1 0 显然在直接利用modred函数进行系统降阶处理时具有一定的盲目性，为此往往将balreal函数与modred函数相结合加以使用。由balreal函数先进行均衡变换，依据Grarn阵确定对系统影响较小的状态，再应用modrea函数求出降阶后的系统。 例3-11 已知系统Σ(A，B，C，D)的系数阵为 A=，B=，C= 在尽可能保持系统基本特性的情况下进行降阶处理。 解：对系统进行均衡变换。 A=[-3 1 0 -l；-0.5 -1 1 -1； -1.5 1 -2 0； -1.5 2 1 -4] B=[1 0 0 0］' C=[1 0 -1 0] [ Ab，Bb，Cb，G，T]＝balreal（A，B，C） 得到均衡变换后的系统模型，以及变换阵T、Gram阵。 Ab= -1.2829 0.4033 -0.1900 0.0359 0.4033 -1.7748 1.5444 -0.3144 -0.1900 1.5444 -5.9201 2.8156 -0.0359 0.3144 -2.8156 -1.0223 Bb= 0.9849 -0.1577 0.0730 0.0138 Cb= 0.9849 -0.1577 0.0730 -0.0138 G= 0.3871 0.0070 0.0005 0.0001 T= 0.6403 -1.6938 1.8373 -2.3068 -0.1591 -1.0219 -0.0109 -0.2640 -0.3446 -1.5361 1.7643 -2.2930 -0.3092 -1.3170 1.5967 -1.4304 由Gram阵可以看出，变换后的模型中，第3、第4状态变量对系统的作用较小，因此可以利用modred函数进行降阶处理，保留第1、第2状态变量。键入下述命令 sys=ss(Ab,Bb,Cb,Db) sysr=modred（sys，3：4，'mdc'） 得到降阶模型为 a= xl x2 xl -l.27804 0.36341 x2 0.36341 -1.44664 b= ul xl 0.98303 x2 -0.14235 c= xl x2 y1 0.98303 -0.14235 d= u1 y1 0.00071 在MATLAB中还给出最小实现函数minreal，它的基本格式为 [Am， Bm， Cm， Dm]= minreal（A，B， C， D） [numm，denm]=minreal（num，den） 该函数消去了不必要的状态，从而得到系统的最小实现。有关它的具体应用可参见相关帮助文件，在此不再详述。 3．1．4 控制系统模型的建立与典型连接 1．基本系统模型的建立 MATLAB为用户提供了一些基本系统或模型建立的函数，如二阶系统、随机n阶系统的模型建立等。下面分别介绍。 1） 二阶系统的生成 在控制系统中，二阶系统占有相当的比例，即或是高阶系统往往也需要对其进行简化降阶处理，然后再对其进行分析研究。因此，研究与讨论二阶系统具有重要的意义。 MATLAB提供了二阶系统生成函数ord2，其基本格式为 [A，B，C，D]=ord2（Wn，z） 式中Wn为自然角频率，z为阻尼因子，返回变量A，B，C，D描述了连续二阶系统。应用该函数就可以生成所期望的以状态方程形式描述的二阶系统，上述自然角频率与阻尼因子的意义表示为  Wn表示，z表示。 同时MATLAB也提供了生成以传递函数形式描述的二阶系统的ord2函数，其格式为 ［num，den］=ord2（Wn，z） 返回变量为传递函数的分子、分母多项式系数向量。 例3－12分别生成以状态方程形式和传递函数形式描述的二阶系统，要求=0 ．8，=2．2 rad／s。 解：直接应用ord2函数进行。 Wn=2．2； z＝0．8； ［A，B，C，D]＝ord2（Wn，z） 得 A＝ 0 1.0000 -4.8400 -3.5200 B= 0 1 C= 1 0 D= 0 求取传递函数形式如下： Wn=2.2； z=0.8； ［num，den］=ord2（Wn， z）； G=tf（num，den） 结果为 Transfer function： 1 s ^ 2 + 3.52 s + 4.84 2） 随机n阶系统的模型建立 MATLAB为用户提供了建立随机n阶系统模型的函数，其基本格式为 [ num， den]= rmodel（n） [num，den］= rmodel（n，p） [A，B，C，D]=rmodel（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）。 运行 [A，B，C，D]＝rmodel（3，l，2） 得到 A= - 0.9688 -4.0932 0.6941 4.0638 -1.0010 -0.4993 -0.8493 0.0983 -2.0227 B= 0 0 0.7143 0.8580 1.6236 1.2540 C= -1.5937 0 0.5711 D= 0 0.6900 3） 系统模型的重构 （1）．子系统选取与删除 MATLAB提供了一个从大系统中选择子系统的函数ssselect，其基本格式为 ［Ae，Be，Ce，De］＝ ssselect（A，B，C；D，Inputs，outputs） 利用指定的输入、输出向量建立状态方程子系统。向量inputs指定系统输入，向量outputs指定系统输出。 [Ae，Be，Ce，De]=ssselect（A，B，C，D，inputs，outputs，states） 利用指定的输入、输出及状态向量构建子系统。 例3-13 已知系统Σ(A，B，C，D)的系数阵为 A=，B=，C=，D= 显然该系统为四阶三输入、二输出系统，现利用ssselect函数在原系统的基础上构造新系统，保留1、3输入信号，1、2输出以及1、2、4状态变量。 解：根据题意要求可以编写相应的程序为 A=[1 5 3 4； 7 16 7 8； 3 0 22 4； 42 25 7 9］； B=[6 41 3；12 25 4；3 8 26；6 47 2]； C=［1 5 9 3；2 3 5 2］； D=[22 46 12；23 54 2］；； inputs＝［1 3］； outputs=[1 2]； states＝［1 2 4]； ［A1，B1，C1，D1］=ssselect（A，B，C，D，inputs，outputs，states） 运行上述语句得到新构建的系统模型为 A1= 1 5 4 7 16 8 42 25 9 B1= 6 3 12 4 6 2 C1= 1 5 3 2 3 2 D1= 22 12 23 2 ssselect函数除了可应用于连续函数，也可应用于离散系统。 MATLAB还提供了与ssse1ect函数进行相反操作的一类函数，它们是 ［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，D）中删除分别由inputs、Outputs和states指定的输入、输出与状态。 （2）．状态的增广 在对系统进行分析研究时，往往需要对状态在系统（输出）中加以增广。如对系统进行全状态反馈研究时，考虑到U=Kx，因此一般就需要在输出方程中增广状态。所以状态增广有时具有十分重要的实际用途，在MATLAB中特别为用户提供了一种状态增广函数augstate。其基本格式为 asys＝ augstate（sys） augstate函数将状态增广到状态空间系统中，产生一个新的状态空间系统，其输入与状态和原系统相同，但其系统输出增加了全部的状态，构造成新系统为 =AX+BU =X+U 例3-14 已知系统状态方程为 A=，B=，C=，D= 将所有的状态增广到系统输出中。 解：应用augstate函数进行输出状态增广。 A=［10 -4 2；8 0 0；0 2 0 ］； B=［2 0 0］'； C=[0.5 -0.4375 0.4375]； D=1； SYS=ss（A，B，C，D） ASYS=augstate（SYS） 执行上述语句得到新构建的系统为 Continuoustime system． a= x1 x2 x3 xl 10.00000 -4.00000 2.00000 x2 8.00000 0 0 x3 0 2.00000 0 b= ul x1 2. 00000 x2 0 x3 0 c= xl x2 x3 y1 0.50000 -0.43750 0.43750 y2 1.00000 0 0 y3 0 1.00000 0 y4 0 0 1.00000 d= ul y1 l.00000 y2 0 y3 0 y4 0 可见在新系统中，状态方程的系数阵没有变化，而输出方程中的C与D阵都出现了增广，满足了将全部状态添加至系统输出的要求。 2．系统组合与连接 所谓系统组合，就是将两个或多个子系统按一定方式加以连接形成新的系统。一般这种连接组合方式主要有串联、并联、反馈等形式。MATLAB提供了进行这类组合连接的相关函数。 1）．模型串联 在MATLAB中提供了模型串联连接函数series，该函数用于两个线性模型的串联，其基本格式为 sys=series（sys1，sys2） 上述串联函数实现了sysl和sys2相串联形成新系统sys。运行后形成 sys=sys1*sys2 如果相串联的两个系统（环节）sysl（s）、sys2（s）的状态方程系数阵分别为（A1，B1，C1，D1）和[A2，B2，C2，D2]，则串联后整个系统的系数阵将变为 A=，B=，C=，D= 对于多入多出系统，上述串联函数的形式为 sys=series（sysl， sys2， ouputsl， inputs2） 应用上述函数将实现由outputsl指定的sysl输出端连接到由inputs2指定的sys2输入端，其中inputs2与 outputs1分别为sys和 sys1的输入、输出向量。 2）．模型并联 在MATLAB中提供了并联连接函数parallel，其基本格式为 sys＝parallel（sysl，sys2，inl，in2，out1，out2） 该函数将两个线性系统以并联方式加以连接， in1与in2指定了相连接的输入端；outl与out2指定了进行信号相加的输出端。最终得到以v1、u、v2、作为输入z1、y、z2作为输出的系统sys。inl、in2分别为sysl、sys2的输入向量，outl、out2分别为sysl、sys2的输出向量。当sysl、sys2为SISO系统时，该函数简化为 sys=parallel（sysl，sys2） 这是标准并联连接的形式。sysl与sys2这两个系统在共同的输入信号作用下，将产生两个输出信号，而并联系统的输出就是这两个系统输出之和。若用传递函数对系统加以描述，系统总的传递函数 G(s)= G1(s)+ G2(s)。若考虑  ，  则系统总传递函数为  例3-15已知两个线性系统 ， 应用parallel函数进行系统并联连接。 解：程序如下。 num1=［12 4］； denl=［1 5 2］； num2=［1 6］； den2=［1 7 l］； sysl= tf（numl， denl）； sys2=tf（num2， den2） ； sys=parallel（sysl，sys2） 运行上述语句得到连接后的总系统传递函数为 Transfer function： 13 s ^ 3 + 99 s ^ 2 + 72 s + 16 s ^ 4 + 12 s ^ 3 + 38 s ^ 2 + 19 s + 2 如果系统用状态方程的形式给出，则并联后的系统模型为 =+ y=+ 若已知两个线性系统，应用parallel函数即可方便地实现系统的连接 A1=［1 2；3 4］；A2=[2 4；3 5]； B1=[2 6]'；B2=［2 5］'； C1=[1 5]；C2=［2 3］； D1=4；D2=1； sysl＝ss（A1，B1，C1，D1）； sys2= ss（A2，B2，C2，D2）； sys= parallel（ sysl， sys2） 运行结果 a= xl x2 x3 x4 xl l.00000 2.00000 0 0 x2 3.00000 4.00000 0 0 x3 0 0 2.00000 4.00000 x4 0 0 3.00000 5.00000 b= ul xl 2.00000 x2 6.00000 x3 2.00000 x4 5.00000 c= xl x2 x3 x4 y1 l.00000 5.00000 2.00000 3.00000 d= ul y1 5.00000 3）．反馈连接 反馈系统是系统中最为重要与常见的一类系统，在MATLAB中为使用者提供了一种方便的构造反馈系统的函数feedback，其一般格式为 sys=feedback（sysl，sys2，sign） 执行该语句将实现两个系统的反馈连接，sign缺省时即为负反馈，sign=1时为正反馈。 如果由sysl与sys2表示的前向系统和反馈系统用传递函数描述，则反馈系统的传递函数为  若由sysl与sys2表示的前向系统和反馈系统分别以状态方程的形式给出（A1，B1，C1，D1），（A2，B2，C2，D2），则反馈系统的状态方程系数阵可表示为 ，， ， 式中， 对于MIMO系统，可以建立更加复杂的反馈系统，其格式为 sys=feedback(sys1,sys2,feedin,feedout,sign) 式中的feedin与feedout均为向量，feedin为sysl的输入向量，指定哪些sysl的输入与反馈环相连接；与此类似的是，feedout指定sysl的哪些输出端用于反馈。最终实现的反馈系统与sys1具有相同的输入、输出端。 例3-16 已知前向环节与反馈环节的状态空间表达式的系数阵分别为 ，，，， ，，， 试将前向环节的输入1和输出2与反馈环节构成负反馈系统。 解：显然由系统可知，环节1为一个多入多出系统，根据要求可以应用反馈函数feedback加以实现。其相应程序如下。 A1=[ 1 0； 0 1]； B1=[1 1；0 1]； C1=[1 3；2 0]；D1=［1 0；2 5］； A2=[-2 0；1 0] ；B2=[1 0]； C2=[0 1]；D2=0； sys1=ss（A1，B1，C1，D1）； sys2=(A2，B2，C2，D2)； feedin=1； feedout=2 ； sign=-1； sys＝feedback（sysl，sys2，feedin，feedout，sign） 运行结果 a= xl x2 x3 x4 xl l.00000 0 0 -1.00000 x2 0 1.00000 0 0 x3 2. 00000 0 -2.00000 -2.00000 x4 0 0 1.00000 0 b= ul u2 x1 l.00000 1.00000 x2 0 1.00000 x3 2.00000 5.00000 x4 0 0 c= x1 x2 x3 x4 y1 l.00000 3.00000 0 -1.00000 y2 2.00000 0 0 -2.00000 d= ul u2 y1 l.00000 0 y2 2.00000 5.00000 由此可以看出，经过构建的新反馈系统，其输入输出端数量没有改变，而状态发生了变化。 4）．系统的复杂连接 在系统分析与研究的实际问题中，往往会遇到一些十分复杂的控制系统，采用上面介绍的方法有时会感到难于处理，需要更加有力的手段和方法。MATLAB为进行这种工作提供了必要的选择，如可利用connect函数以及nblocks—blkbuild所构造的程序模块，实现由众多环节组合连接成复杂系统。这部分内容将在下面章节中介绍。 3．2 连续线性系统的计算机辅助分析与设计 在经典控制理论中，频率特性及根轨迹分析方法是分析和设计连续系统的基本方法。在分析自动控制系统的工作性能时，最直观的方法是求出它的时域响应特性。但是，对于高阶控制系统的时域特性，很难用分析方法确定，而频域中的一些图解法则可以比较方便地用于控制系统的分析和设计。用频率特性方法对控制系统作出了分析和设计之后，再根据时域和频域之间的关系就能确定系统的时域响应特性。因此，用频率特性方法分析和设计控制系统有以下优点： ① 不必求解闭环系统的特征方程式，只要用一些较为简单的图解法就可以研究它的绝对稳定性和相对稳定性； ② 系统的频率特性可以用实验方法测出，这对于某些很难确定动态方程的系统或环节是非常有用的； ③ 用频率法设计系统，对不希望的噪声可忽略不计，在规定的程度之内，频率法还可以用于某些非线性系统。 在控制系统分析中，为了避开直接求解高阶特征方程式的根时遇到的困难，在实践中提出了一种图解求根法，即根轨迹法。所谓根轨迹是指当系统的某一个(或几个)参数从变化到时，闭环特征方程的根在根平面上描绘的一些曲线，应用这些曲线，可以根据某个参数确定相应的特征根。在根轨迹法中，一般取系统的开环传递系数K(又称开环系统的放大倍数)作为可变参数。 由于根轨迹是以K为可变参数，根据开环系统的零极点画出来的，因而它能反映出开环系统零极点与闭环系统的极点(特征根)之间的关系。利用根轨迹可以分析系统参数和结构已定的闭环系统的时域响应特性，以及参数变化对时域响应特性的影响，而且还可以根据对时域响应特性的要求确定可变参数及调整开环系统的零极点位置，并改变它们的个数。也就是说，根轨迹法可用于解决线性系统的分析和综合问题。 根轨迹法和频率特性法都具有直观的特点。由于根轨迹和闭环系统的动态响应有着直接联系，所以只要对根轨迹进行观察，用不着进行复杂计算就可以看出动态响应的主要特征。 下面首先介绍用根轨迹法及频率特性法分析系统的简单原理，然后介绍利用MATLAB绘制很轨迹、奈奎斯特图和伯德图，并求系统稳定裕量的方法，最后利用这些工具对系统进行校正设计。 3．3 线性系统根轨迹的绘制 极轨迹法是分析和设计线性定常控制系统的图解方法，使用十分简便。特别是在多回路系统研究中，应用根轨迹比用其他方法更方便，因此在工程实践中获得了广泛应用。本节主要介绍根轨迹的基本概念，并从闭环零极点与开环零极点之间的关系推导出根轨迹方程，然后将向量形式的根轨迹方程转化为常用的相角条件和模值条件形式，利用MATLAB提供的根轨迹函数绘制跟轨迹，最后介绍设计工具。 3．3．1根轨迹方程 根轨迹是系统所有闭环极点的集合。 设控制系统 如图3-1所示，其闭环传递函数为  (3-7) 为了用图解法确定所有闭环极点， 令闭环传递函数表达式(3-7)的分母为零， 得闭环系统特征方程  (3-8) 当系统有m个开环零点和n个开环极点时，(3-8)式等价为  (3-9) 式中，zj为已知的开环零点，pi为已知的开环极点，K可从零变为无穷。把(3-9)式称为根轨迹方程。根据(3-9)式，可以画出当K从零变到无穷时系统的连续根轨迹。应当指出，只要闭环特征方程能化成(3-9)式的形式，都可以绘制根轨迹，其中处于变动地位的实参数不限定是根轨迹增益K，也可以是系统其他变化参数。但是，用(3-9)式表达的开环零点和开环极点，在S平面上的位置必须是确定的，否则无法绘制根轨迹。此外，如果需要绘制一个以上参数变化时的根轨迹图，那么画出的不再是简单的根轨迹，而是根轨迹簇。 根轨迹方程实质上是一个向量方程，直接使用很不方便。考虑到   因此，根轨迹方程(43)可用如下两个方程描述   (3-10) 和  (3-11) (3－10)式和(3-11)式是根轨迹上点应该同时满足的两个条件，前者称为相角条件，后者叫做模值条件。根据这两个条件，可以完全确定s平面上的极轨迹和根轨迹上对应的K值。应当指出，相角条件是确定s平面上根轨迹的充分必要条件。这就是说，绘制根轨迹时，只需要使用相角条件；而当需要确定根轨迹上各点的K值时，才使用模值条件。 3．3．2用 MATLAB绘制极轨迹 MATLAB 5.3的控制系统工具箱给出了一系列关于根轨迹的函数，如表 3-1所示。使用这些函数可以方便地绘制根轨迹图，以及对系统进行分析和设计。 表3－1根轨迹函数 函数 功能 函数 功能  pzamp 绘制零极点图 rlocus 计算并绘制Evans根轨迹  rtool 根轨迹设计GUI工具 sgrid 在根轨迹或零极点图中添加s平面网格  rlocfind 从根轨迹图中选择反馈增益 zgrid 在根轨迹或零极点图中添加z平面网格  1．计算并绘制极轨迹 函数rlocus计算SISO开环模型的Evans根轨迹，根轨迹以反馈增益的函数形式给出了闭环极点的轨道(假定为负反馈)。该函数的一般用法为 rlocus(G) 其功能是计算并绘制开环SISO模型G的根轨迹，G可以是如图4-2所示的三种模型中的任何一种(其中 K为常数)。对于图(a)模型，；对于图(b)模型，；对于图(c)模型，。  由于系统闭环传递函数分别为，，。因此若G能表达为  则  方程的根就是系统闭环极点。rlocus函数可以自适应选择一组正增益K，产生平滑的曲线，用户也可以自己指定向量K，格式为 rlocus(G,K) 函数rlocus还可以带输出参数，格式为 [r，K]= rlocus(G)或r= rlocus(G，K) 这时，返回自适应选择的增益K和对于这些增益的复数根的位置r。矩阵r的列数和增益K的长度相同，它的第m列的元素是对于增益K(m)的闭环系统的根。 2．选择反馈增益 函数rlocfind返回在根轨迹中和一些极点相关的反馈增益。用法如下。 1) [K，poles]= rlocfind(G)从函数rlocus生成的模型G的根轨迹图中选择反馈增益。运行rlocfind函数后，在根轨迹图上出现十字形光标，用来选择极点位置。与选中的点相关联的根轨迹增益返回到K中，列向量poles包含了该增益的闭环极点。使用该函数时，模型G的根轨迹留必须在当前图形窗口。 2) [ K，Poles]=rlocfind(G, p)用向量P 指定期望的极点位置，并计算这些位置上的根轨迹增益。向量K的第m项是根据极点位置p(m)计算的增益，矩阵poles的第m列是相应的闭环极点。 例3-17 至绘制如下系统的根轨迹图  解：下面的程序可以计算并绘制根轨迹图， 图3－3所示，然后用rlocfind函数在图中选择极点 位置(十字光标如图所示)，得到反馈增益的值。 G= tf([ 0.05,0.045],conv([l ,-l.8,0.9 ] [l,5,6])) rlocus( G) K＝rlocfind(G) %执行后出现下面一条信息，等待选择权点位置 select a point In the graphics window selected_point＝ %选择极点位置后得到位置和增益值 -0.8894 -0.0292i K= 4.9655e+003 3．3．3根轨迹设计工具 控制系统工具箱中提供了控制系统设计的 GUI工具——根轨迹设计GUI(Root LocusDesign GUI)。该工具用根轨迹法交互式设计 SISO补偿器。 调用根轨迹设计工具的方法是在MATLAB命令窗口输入命令rltool，其窗口如图3-4所示。 1．窗口功能 注意在窗口右上角的系统框图，图中的模块“p”为被控对象的模型，模块“K”就是设计补偿器模型，窗口左上角显示着设计过程中补偿器模型(ZPK模型)。设计的最终目的就是为被控对象“p”设计出一个适当的补偿器“K”，使得整个系统(右上角的系统框图)的性能满足要求。 根轨迹设计GUI中间的窗口是进行设计的主窗口，它显示整个系统的根轨迹留(由于还没有设计补偿器，而且“F”和“H”都等于1，因而图4-4中的根轨迹图就是被控对象的根轨迹)。 主窗口左上角有四个按钮，分别对应不同的编辑功能。 ①缺省时按钮 被按下，这时在主窗口中可以用鼠标交互式地移动闭环极点(根轨迹上的红色方块)和设计出的补偿器的零极点。 ②按钮 或 被按下时，可以在主窗口 中添加补偿器的极点或零点。 ③按钮 被按下时，可以在主窗口中檫除补偿器的极点或零点。 主窗口的左下方是坐标轴设置按钮，它们 的作用分别是： ① 　 储存当前的坐标范围； ② 载入存储坐标范围； ③ 　 把坐标轴设置为正方形； 设置坐标轴的纵横比相等。 主窗口右下方是放大工具按钮， 它的作用分别为： ① 　沿X方向和Y方向都放大 ② 　 沿X方向放大； ③ 沿Y方向放大； 观察全部图形。 根轨迹设计 GUI下方的复选框 　　　　　　　　　 ，用于选择需要显示的系统时间响应和频率响应分析。选中其中的一个或几个响应后，将调用 LTI Viewer来显示响应的曲线，随着设计过程中系统模型的变化，响应也随之变化，这样便于在设计过程中随时观察系统的特性。 2．输入输出模型 根轨迹设计GUI窗口菜单File下的 Import Model用来输入 LTI设计模型，其窗口如图3-5所示。窗口左上角显示的是系统的结构框图，它允许两种反馈结构，另一种如图3-4右上角所示。在这里输入的模型包括被控对象P、前置滤波器F、传感器动态特性H。缺省情况下，P、F、H的值是1。  图3-5载入LTI设计模型 在窗口左下角指定输入模型的位置，包括工作空间、MAT文件和Simulink框图。窗口中间显示了指定位置(如工作空间)中的所有模型名称。 如果要使被控对象P的模型为G1，则首先选定模型G1，然后在窗口右下方P前面的箭头上点一下，被控对象P的模型就成为G1。同样可以为H和F指定相应的模型。系统设置完成后，按下“OK”按钮返回主窗口，这时主窗口中间会显示出输入的LTI模型系统的根轨迹图(图3－4)，其中红色的小方块表示闭环极点。 类似地，菜单File下的 Import Compensator选项用来输入要进行设计的补偿器K的模型。模型输入后，主窗口中显示出整个系统的根轨迹图，在图中可以对补偿器K的零极点进行增删和编辑，但不能编辑原系统的零极点，只能沿轨迹移动闭环极点。 模型的输出用菜单File下的Export选项，其窗口如图3－6所示。模型可以输出到磁盘MAT文件，也可以输出到工作空间。窗口左侧显示的几项内容的含义分别为 ：DM指设计的模型；OL指开环模型；CL指闭环模型；K指补偿器模型。 3．设置网格和边界线 利用 Tools菜单下的 Add Grid/Boundary选 项可以在根轨迹图中添加网格并设置边界线。其 设置窗口如图3－7所示。 例如，需要显示阻尼比和固有频率的网格，则只要选中该项前的单选按钮，然后选择下面的 “Grid On？”前的复选框即可。 同时，设计的系统需要满足条件：稳定时间小于0.04s，阻尼比小于0.4，那么只要在该窗口的 Setting Time后输入 0.04， 在 Damping Ratio后输入0.4即可。设置的内容见图 3－7，设置后的根轨迹图如图3－8所示。  图3－7网格和边界线设置窗 图3－8添加了网格和边界的根轨迹图 4．转换为Simulink框目 主窗口 File菜单下的 Draw Simulink Diagram选项能够把 LTI系统模型转换为相应的Simulink仿真框图，这样便于在Simulink中对系统作进一步的仿真分析。 3．4 线性系统频率特性与奈奎斯特图 频域分析法是利用频率特性研究线性系统的一种经典方法。绘制奈奎斯特(Nyquist)图是频率法中利用控制系统的开环幅相频率特性，判断其闭环系统稳定性的重要方法之一，它是由著名的奈奎斯特判据而来的。开环系统的幅相频率特性容易计算，且可通过实验求得，因而奈氏判据使用方便，同时物理意义明确。这个判据确定了开环系统的频率持性与闭环系统动态响应之间的联系，它不仅能判断闭环系统的稳定性，而且可利用它找到改善闭环系统动态响应的方法，因而它是构成经典设计方法的重要基础。 3．4．1奈奎斯特稳定判据 奈奎斯特判指定义如下： ①开环系统稳定时，即开环系统没有极点在右半根平面，如果相应于频率ω从-∞→+∞变化时，开环系统频率特性 曲线不环绕(-l，j0)点 ，那么闭环系统就是稳定的，否则是不稳定的； ②开环系统不稳定时．即开环系统有P个极点在右半根平面时，如果相应于频率ω从-∞→+∞变化时．开环系统频率特性 曲线反时针方向环绕(-1，j0)点的次数N等于位于右半根平面内的开环系统的极点数P，那么闭环系统就是稳定的，否则是不稳定的。 奈奎斯特图就是根据奈氏判据，在开环系统频率特性的平面上绘制开环系统幅相频率特性的极坐标图。 例3-18 设有一控制系统的开环系统传递函数  试判别其闭环系统的稳定性。 解：开环系统的频率特性为  当ω=0时, 0° 当ω→+∞时,  180° 上面这两个特殊点确定了奈奎斯特曲线的变化趋势，再计算对应几个ω的的值，便能绘出如图3－9所示的开环系统幅相频率特性曲线图，即奈氏图。因为例中给出的开环传递函数没有极点在右半根平面，而且不论K怎样变化，曲线都不会包围(-1,j0)点，因此闭环系统总是稳定的 。 下面通过例题准确绘制奈奎斯特曲线。 例3-19设一控制系统的开环传递函数为 　　　　　 　　(3-12) 求其模值和相角的表达式。 　　解：(3-6)式可以写成 = 这里 =  利用上面的表达式编制计算程序，指定具体数据，就可以绘制出奈奎斯特图了。 3．4．2利用奈奎斯特图来分析系统的稳定裕量 绘制控制系统的奈奎斯特图，不仅可以确定系统 的稳定性，还可以确定系统的稳定裕量。奈氏曲线离 平面上的(-1,j0)点越远，说明稳定裕量 越大。稳定裕量是以轨迹上两个特殊点 的位置来度量的。所谓特殊点是指如图4-10所示奈奎 斯特图的A点和B点。A点规定为的 轨迹和单位圆的交点(交点频率设为ωc)。显然，在A 点，的幅值为，而相角为 (相角规定为顺时针为负，逆时针为正)。 相位裕量定义为向量OA与负实轴间的夹角，并用γ表示。根据图3-10，可由下式确定： 　　　　　 　　　　　　　 　　　(3-13) 对最小相位系统来说； γ＞0，表示相位裕量有正值，说明系统趋向稳定； γ＜0，表示相位裕量为负值，说明系统趋向不稳定； γ＝0，表示A点在负实轴上，系统为临界稳定，这时奈氏曲线穿越(-1,j0)点。 通常γ的值在30°～180°之间。 B点规定为G(jω)H(jω)轨迹和负实轴的交点，交点频率设为ωg。 增益裕量(又称为幅值裕量)定义为G(jω)H(jω)轨迹在B点的幅值A(ωg)的倒数，并Kg表示，即  (3-14) 其中，Kg＞1，表示系统趋向于稳定； Kg＜1，表示系统趋向于不稳定； Kg＝1，表示奈氏曲线穿越(-1,j0)点，这时系统为临界稳定。 　　增益裕量Kg也可以用分贝表示，这时　　　　　　 　　　　　　　　　　　　　　　dB 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(1,j,:)和ry(1,j,:)给出。对于单变量系统，实部和虚部分别由rx(l,l,:)和ry(1.1,:)给出，更简单地，可以由rx＝rx(:)和ry=ry(:)来提取实部和虚部向量。 该函数还可以由下面各种其他形式来调用 [rx,ry,w]＝ nyquist(G)，或 nyquist(G) 与其他控制系统工具箱中的函数类似，前一种格式会自动选择频率向量，而后一种格式能自动绘制出系统的条奎斯特图。 例3-20系统开环传递函数  绘制系统的条奎斯特图，并讨论其稳定性。 解：用下面的MATLAB命令直接绘制出系统的条奎斯特图，如图3－11(a)所示。 G= tf(1000,conv([1,3,2] ,[ 1,5]))； nyqulst(G)；axis(‘square’) 该图中(-1，j0)点附近奈奎斯特图的情况并不很清楚，可以用CtrlLAB程序对得出的奈奎斯特图进行局部放大(选择图形窗口中的Zoom in菜单项，或用快捷按钮 )，则可以得出如图3－11(b)所示的局部放大结果。从局部放大的图形中可 以看出．来奎斯特图逆时件包围(-1，j0)点2次，而原开环系统中没有不稳定极点，从而可以得出结论，闭环系统有2个不稳定极点。这也可以由下面的MATLAB语句来进一步验证。 G_close= feedback(G，1)； roots(G_clos,den|1|) ans= -12.8196 2.4098 + 8.5427i 2.4098 - 8.5427i  图3－11例3-20的奈奎斯特囹 3．4．4幅值裕量与相位裕量 控制系统工具箱中提供了margin函数来求取给定线性系统的幅值裕量和相位裕量，该函数可以由 [Gm，Pm ，wcg ，Wcp]= margin(G) 格式来调用。可以看出，幅值裕量与相位裕量可以由LTI对象G求出。返回的变量(Gm，Wcg)为幅值裕量的值与相应的频率，而(Pm，Wcg)则为系统的相位裕量与相应的频率。若得出的裕量为无穷大，则得出的值为inf，这时相应的频率值为NaN(表示非数值，这两个量均是MATLAB保留的常数)。 如果已知系统的频域响应数据，还可以由下面的格式 [Gm，Pm，Wcg，wcp]= margin(mag，phse,w)； 调用此函数，其中(mag,phase，w)分别为频域响应的幅值、相位与频率向量。 例3-21考虑如下系统模型，求它的幅值裕量和相位裕量，并求其闭环阶跃响应。  解：由下面的MATLAB语句直接求出系统的幅值裕量与相位裕量及闭环阶跃响应。 G=tf(3.5,[1,2,3,2]) [Gm，Pm，Wcg, Wcp]＝margin(G)；]Gm，Pm，Wcg，Wcp] ans＝ 1.1429 7.1578 1.7321 1.6542 G_close＝feedback(G；1); step(G_close) 得出闭环系统的阶跃响应曲线如图3－12(a)所示，系统的幅值裕量很接近稳定的边界点1，且相位裕量只有7.578°所以尽管闭环系统稳定，但其性能将不会太好。可以看出在闭环系统的阶跃响应中有较强的振荡，这就验证了前面得出的结论。  如果系统的相位裕量y＞45°,则称该系统有较好的相位格量，然而这样的数值并不很绝对。 考虑一个新的系统模型  可以由下面的MATLAB语句直接求出系统的幅值裕量与相位裕量为 G＝tf(100*conv([1,5],[1,5]),conv([1,1],[1,1,9]); [Gm，Pm，Wcg，WcP]＝ margin(G)；[Gm，Pm，Wcg，Wcp] ans= Inf 85.4356 NaN 100.3281 G_close＝feedback(G,1)；step(G_close) 并可以得出系统的闭环阶跃响应如图4-12(b)所示。可以看出，该系统有无穷大的幅值裕量,且相位格量高达85.435 6°，所以系统的闭环响应是较理想的 3．5 线性系统的伯德图分析 3．5．1伯德图 利用奈奎斯特图可以从开环幅相频率特性分析闭环系统的稳定性及稳定裕量，但是奈氏曲线不能明显地表示系统增益及时间常数等参数变化时对系统性能的影响。为了克服奈奎斯特图分析设计系统的缺点，将系统的开环频率特性的幅频特性曲线和相频特性曲线分开画，并对幅频特性进行对数运算(相频特性不取对数)，横坐标表示频率ω，纵坐标分别表示对数幅值和相角，分别画出的对数幅频特性曲线图和相频特性曲线图称为伯德(Bode)图。 为了能把一个较宽频率范围的图形紧凑地表示在一张图上，横生标采用对数分度，而表示对数幅值和相角的纵坐标仍然用线性分度．这样就形成了一种半对数坐标。画在半对数坐标L的对数幅频特性和相频特性又称对数坐标图。因此，伯德图又称对数坐标图。 下面举例说明伯德图的画法。 例3-22 设一控制系统的开环传递函数为  (3－15) 试绘制相应的伯德图。 解：其幅频和相频特性分别为  (3－16)  (3－17) 绘制伯德图时．对幅值两边取常用对数．得  (3－18) 这样就把幅值的乘除运算变成了加减运算，使计算简化。 通常在对数前面还乘上系数20以使对数幅频的计算单位变成分贝dB表示、为了方便，以分贝为单位的对数幅频特性常用L(ω)来表示即(3－18)或可改写为  (3－19) 设(3－9)式中的、、， 则由(3－17)式和(3－19)式可画出系统的 伯德图，如图3-13所示。 3．5．2伯德图的渐近线 从前面的分析中可以看到，要想给制比较复 杂的伯德图，需要借助于计算机才能实现。事实 上，如果只想定性地研究系统的频域响应，则使 用伯德图的渐近线就足够了。系统伯德图的渐近 线可以由比较简单、规范的方法容易地绘制出来。 这里首先考虑幅频特性渐近线的切换点的规则。 当频率值达到G(s)的一个单独极点时，幅频传性的渐近线的斜率将减小20 dB/十倍频；若频率值到达了一对共轭的极点时，则渐近线的斜率值减小到40db/十倍频时，如果频率值到达一个p重极点时，福频特性的渐进线的斜率减小20p dB/十倍频；若频率值到达了一对共轭的零点时，则渐近线的低利率将增加40dB/十倍频；如果频率值到达一个p重的零点，则渐近线的斜率将增加20p dB/十倍频；如果系统含有m个微分和n个积分环节．则斯近线的起始频率应该是 20(m-n) dB／十倍频。 根据上述原则可以编写一个MATLAB函数bode_asym来求取系统的伯德幅频特性渐近线的数值，该函数的调用格式为 [wpos,ypos]＝bode_asym(G,w) 使用bode_sym函数的程序如下： function[wPos,ypos]＝bode_sym(G,w) G1＝zpk(G);wpos＝[];slope＝[]; If nargin=＝ 1, w＝ freqint2(G); end zero＝G1.z{1}; polar＝ G1.p{1};gain＝ G1.k; num d＝ 0; numl＝ 0; %找出非零极点及相应的斜率，记录零极点个数 for i＝1：length(zero); if isreal(zero( i)) if zero(i)==0 num_d＝ num_d＋ 1; else wpos＝[wpos,abs(zero(i))]; slope=[slope, 20]; end, else if imag(zero(i))＞ 0 wpos＝[wpos,abs(zero(i))]; slope＝[slope, 40]; end, end, end for i＝1：length(polar); if isreal(polar(i)) if polar(i)＝＝0 num_i＝ num_i＋ 1; else wpos＝[wpos,abs(polar(i))]; slope＝[slope,-20]; end,else if imag(polar( i))＞ 0 wpos＝[wpos,abs(polar(i)]; slope＝[slope,一40]; end,end,end %加入首末点，并对转折点排序 wpos＝[w(1),wpos w(length(w))]; slope＝[(num_d－num_i)*20,slope,0]; [Wpos,I_seq]=srot(wpos); slope=slope(i_seq); % 合并相同的零极点 Wposl(1)= wpos(1); slopel(1)＝slope(1); j＝2; for i=2:length(wpos) if wpos(i)==wpos(i-1) slopel(j－1)=slopel(j－1)＋ slope(i) i= i十1; else wpos1(j)=wpos(i)；Slope1(i)=slope(i) j＝j+1; end end wpos=wpos1; slope=slope1; [yslope,pp]=bode(G,w(1)); ypos(1)=20* log10(yslope); for i＝2:lengt(wpos) slopes=sum(slope(1:i－1)); ypos(i)＝ypos(i一l)＋slopes*log10(wpos(i)/ wpos(i-1)); end if nargout== 0 [x0,y0,wl]=bode(G); x=wpos;y=ypos; subplot(212),semilogx(wl,y0(:)) subplot(211),semilogx(wl,20*log10(x0(:)),x,y) end 其中G为系统的LTI对象模型，变量W为用户指定的频率点构成的向量，返回的变量为(wpos，ypos)渐近线的转折点。获得了这些数据之后，可以用semilogx(wpos,ypos)命令绘制出系统的伯德幅频特性渐近线。如果调用bode＿asym函数时不给出返回变量，则带有渐近线的幅频特性伯德图将自动地绘制出来。 3．5．3由伯德图判定系统稳定性 因为系统的伯德图只是频域响应的另外一种描述，它完全等效于一般的奈奎斯特图形表示，所以奈奎斯特定理的结论可以由伯德图的形式叙述。 ① 对开环稳定的模型，若假设为单位负反馈结构，则当且仅当相位曲线在幅值曲线进入负值之前，不穿越-180°线，则闭环系统是稳定的；若在相位穿越-180°线时。幅值曲线在变成负值之前以斜率为-20p dB/十倍频的速度下降，则闭环系统有p个不稳定假定。 ② 对含有p个不稳定极点的开环对象模型，若在相位穿越-180°线时，幅值曲线在变成负值之前以斜率为-20dB／十倍频的速度下降，则闭环系统稳定。 3．5．4应用MATLAB函数绘制伯德图 在表3－2中已经给出了MATLAB相关函数。MATLAB控制系统工具箱中提供了bode函数来求取、绘制给定线性系统的伯德图，该函数可以由下面的格式来调用 [mag，Pha]=bode(G,w) 其中G为系统的LTI对象模型，变量W为给定的频率点向量，该函数在这些频率点上对系统进行频域分析。返回的变量对[mag,pha]分别为系统的幅值和相位向量，其中相位的单位为角度。这里得出的[mag,pha]仍为二维矩阵。 在绘制伯德图时还经常需要使用MATLAB命令 magdB=20*log10(mag) 将得出的幅值向量变换成分贝的形式。 该函数还可以使用下面的格式进行调用 [ mag，pha，w] = bode(G) 其中频率点向量W可以自动地由bode函数根据系统模型产生。若在调用bode函数时不返回任何变量，则将自动绘制出系统的伯德图，而不再向MATLAB工作空间返回任何信息。 例3-23考虑二阶系统传递函数模型  试用MATLAB绘制不同ζ和ωn时的伯德图。 解：可以从模型直接得出系统频域响应的幅值和相位为  (3-20a)  (3-20b) 幅值用分贝表示成  (3-21) 将自然频率ωn的值固定为ωn=1，改变阻尼ζ的值，得出的系统伯德图如图3-14 (a)所示；再将阻尼比ζ的值固定为ζ=0.707，改变自然频率ωn的值，得到系统的伯德图如图3－14(b)所示。 程序如下 wn=1;zet=[0:0.1:1,2,3,5]; w=logspace(－1,l,100); mm=[]; pp=[] ; for i=l:length( zet) m=wn2./sqrt((wn^2-w.^2).^2 ＋ 4*zet(i)^2*wn^2*w.^2) ; p=atan2(-2*zet(i)*wn*w, wn^2- w.^2); pp=[pp,180* p’ /pi]; mm=[ mm,20*log10(m’)]; end subplot(2ll);semilogx(w,mm), Subplot(212);semilogx(w,pp) wn＝[0.1:0.1:1];zet＝0.707; w=logspace(－l,1,100) ; mm=[]; pp=[]; for i＝ 1:length(wu) nn= wn(i)^ 2;dd=[l,2* zet*wn(i),wn(i)^2]; [m,p]＝ bode(nn,dd, w); pp=[pp,p];mm=[mm,m」; end subplot(211);semilogx(w,20*log10(mm)), subplot(212);semilogx(w,pp) 从图(a)可以看出，当阻尼比ζ较小时，则系统频域响应在自然频率ωn附近将表现出较强的振荡，该现象称为谐振。在图(b)中，当自然频率ωn的值增加时，伯德图的带宽将增加，该现象将使得系统的时域响应速度变快。  还应注意，二阶系统的谐振频率ωr，总是等于系统的自然频率ωn，闭环系统阶跃响应的幅值Mr，也总是等于系统频域响应的最大幅值。 例3-24传递函数模型  绘制出系统的伯德图及其渐近线。 解：使用下面的MATLAB命令绘制出系统的伯德图及其渐近线，如图 4－15所示。 G=tf(1,[conv([1,1] [1,2],0 )]) [x0，y0，w]=bode(G);[x，y] ＝bode_asym (G．W)： subplot(211),semilogx(w,20*log10(x0(:)), x,y) subplot(212),semilogX(W,y0(：)) 可以看出，因为系统在S=0处有极点．故伯德图的渐近线的斜率初值为-20 dB/十倍频，在整个渐近线上还有两个转折点，分别在ω＝1 和ω=2处，这是因为系统在这两个位置含有极点。 例3-25系统开环模型  绘出伯德图和渐近线，并讨论其稳定性。 解：可以由下面的MATLAB语句绘制出系统的伯德图，如图4－16(a)所示。 G= tf(1000，conv([1,1],conv([1,2]),[1.5]))) ; [x0,y0,w]=bode(G) ; [x,y]=bode＿asym(G,w); subplot(211),semilogx(w,20(log10(x0),x,y) subplot(212),semilogx(w,y0) 由图可以看出。因为相频特性在4 rad/s处穿越180°线，该频率低于幅值穿越0 dB线时的值(其值为ω＝9 rad/s)，所以闭环系统是不稳定的。 通过进一步分析可以看出，由于ω=4 rad/s时幅频特性的斜率为-40 dB/十倍频，所以可以发现闭环系统将有两个不稳定极点。系统的稳定性分析结果可以由下面的MATLAB语句直接得出。 rolcus(G);[k,pp]=rlocfind(G) Select a point In the graphics window selected_point＝0.0231+4.13491 K＝ 135.5319 pp= -8.1151 0.0576 + 4.2345i 0.0576 - 4.2345i  这样绘制出的根轨迹如图4－16(b)所示。从该图可以看出，使得闭环系统稳定的最大增益为K＝135.5319，否则闭环系统将有两个不稳定复数共轭极点。 3．6 基于复域与顾域的线性系统设计与校正 前面各节的内容主要集中于解决反馈控制系统分析，本节将介绍反馈控制系统的设计问题。事实上，控制系统的设计可以认为是系统分析的逆问题，因为在系统分析中，常假设系统的控制器是已知的。在控制系统设计中，将研究如何对给定对象模型找出其控制器的策略。 首先介绍串联补偿器的频域设计方法，这样的频域补偿器中最有效的一种是所谓的超前/滞后补偿器，又称控制系统的频域综合问题。然后介绍最常用的工程校正方法——比例积分微分(PID)接制，并用 Ziegler－Nichols方法和解析法对 PID控制器进行参数整定。 3．6．1串联超前/滞后补偿器设计 相位超前和滞后补偿器通常可以等效地表示为由电阻和电容元件构成的无源网络形式，这样的网络又称为RC网络。串联补偿器主要有三种形式：相位超前补偿器、相位滞后补偿器和相位超前/滞后补偿器。下面详细分析这三种补偿器。 1．超前/滞后环节基本特性 1）．相位超前补偿器 该补偿器的等效超前(RC)网络表示如图 3－17(a)所示。简记阻抗力为，则这样网络的传递函数模型可以写为  (3-22) 式中  (3-23) 可以看出总有>1。  下面给出相位超前补偿器的重一般公式  　　　　　 　　　(3－24) 相位超前补偿器的零极点位置如图 3－17(b)所示。因为有α＞1，所以在s平面上极点的位置总在补偿器零点的左边。假设T=1，使用下面的MATLAB语句，可以绘制出不同值下补偿器的伯德图和奈奎斯特图，如图3-18(a)，(b)所示。 alpha＝1.5;T=1 ;alpha0=1.5：0.5：5 G=tf([alpha* T,1]/alpha,[T,1]); [x,y,w1]=nyquist(G,｛.01,1}); [m,p,w2]=bode(G,{.01,10}); x＝x(:); y＝y(:); m=m(:)’; p=p(:)’; for i=1:length( alpha0) G1=tf[alpha0(i)*T,1]/alpha0(i),[T,1]]; [xx, yy]＝nyquist(G1, w1) ; [mm,pp]=bode(G1,w2); x=[x,xx(:)];y=[y,yy(:)];  m＝[m;mm(:)]; p[p;pp(:)’]; end subplot(211); semilogx(w2,20*log10(m)), subplot(212); semilogx(w2,P) figure，plot([x，x]，[y,－y]), axis([0, 1,－0.4, 0.4]): axis(‘square’) 一般来说，超前补偿使系统的相角裕量增加，从而提高系统的相对稳定性。对于给定的系统增益K，超前补偿增加了系统的稳态误差。为减小稳态误差，必须使用大增益的补偿器。超前补偿器也使增益穿越频率ωc增加，这样会使阶跃响应过渡过程加快，增大了系统带宽。 2）．相让滞后补偿 相位滞后补偿器的等效RC网络如图3-19(a)所示，而该结构的零极点位置如图3-19(b)所示。简记阻抗Z1=R1与Z2=R2+1/(Cs),则这样网络的传递函数可以写成   式中: ＜1。相位滞后补偿器的更一般的形式可以写成  (3-26) 设T=1 ，可以由下面的程序绘制出各种α值下的伯德图与奈奎斯特图，如图3－20(a)、(b)所示。 alpha＝0.9;T=1 ;alpha0=0.8:-0.1:0.1 G=tf([alpha* T,1]/alpha,[T,1]); [x,y,w1]=nyquist(G,｛.01,100｝); [m,p,w2]=bode(G,{.01,100}); x＝x(:); y＝y(:); m=m(:)’; p=p(:)’; for i=1:length( alpha0) G1=tf[alpha0(i)*T,1]/alpha0(i),[T,1]]; [xx, yy]＝nyquist(G1, w1) ;[mm,pp]=bode(G1,w2); m＝[m;mm(:)’]; p=[p;pp(:)’]; end subplot(211); semilogx(w2，20*log10(m))， subplot(212); semilogx(w2，P) figure,plot([x,x],[y,－y]), axis([0,10,5,5]): axis(‘square’) 可以看出，和相位超前补偿器不同，在滞后作用下系统的相位减小了，同时低频段增益得到增加。  3）．相位超前/滞后补偿 这种补偿器的等效RC网络如图3－21(a)所示，零极点位置如图3－21(b)所示。简记阻抗，则这样网络的传递函数模型可以写成  (3-27) 式中： ,则有 可以看出，α＞ l且β＜ 1时，使得(3－21)式中的第一项具有超前性质，而第二项具有滞后性质。  若分别选择T1 =0.5,T2=0.005且 α=3,β=1/3,则补偿器的伯德图及其渐近线可以由下面的MATLAB命令得出，如图3－22(a)所示。 T1=0.5;T2=0.005;alpha=3;beta=1/3; nc=conv([alpha*alpha*T1,1],[beta*T2,1]); dc=conv([T1,1],[T2,1]); w=logspace(-2,4);G=tf(nc,dc); [m,p]=bode(G,w); [x,y]=bd_asymp(G,w);m=m(:);p=(:); subplot(211);semilogx(w2，20*log10(m),x,y)， subplot(212);semilogx(w2，p) 如果选择了T1 =0.5,T2=0.005，则可以由下面的MATLAB命令绘制出补偿器的伯德图及其渐近线，如图3－22(b)所示。 可以看出，滞后补偿器在低频时有较大的增益，这样就可以降低稳态误差。因此滞后补偿适用于需要经常限制闭环系统的带宽，以限制外部噪声信号进入控制环的情况。 例3－26 试为如下对象环节设计超前补偿器。  解：首先可以用下面的MATLAB语句得出原系统的幅值裕量与相位裕量。 G=tf(100,[0.04,1,0]); [Gm,Pm,Wcg,Wcp]=margin(G);[Gm,Pm,Wcg,Wcp] ans ＝ Inf 28.0202 NaN 46.9782 W=logspace(-1,3);bade(G,w);  可以看出，这个模型有无穷大的幅值裕量，而其相位裕量为γ= 28°其发生的频率为ω＝47 rad/s。绘制出原系统的伯德图，如图3-23(a)所示，图中标注了相位裕量。如果这时想引入一个超前补偿器来增大相位裕量，如假设控制器为Gc(s)=(0.0262s+1)/(0.0106s+1)，则通过下面的MATLAB语句得出校正后系统的幅值裕量和相位裕量。 Gc1=tf([0.0262,1],[0.0106,1]);bode(Gc1,w); G_op=Gc1*G;[Gm,Pm,Wcg,Wcp]=margin(G_op);[Gm,Pm,Wcg,Wcp] ans ＝ Inf 47.5917 NaN 60.3252 起前补偿器模型的伯德图如图3-23(b)所示。可以看出，在频率ω＝60 rad/ s处系统的幅值和相位均增加了，在这样的控制器下，校正后系统的相位裕量增加到47.6，而剪切频率增加到ω=60 rad/s。通过下面的 MATLAB语句绘制出补偿后系统的伯德图，如图3－24(a)所示，在图中还同时绘制出了原系统的伯德图。 [m,p]=bode(G,w);[m1,p1]=bode(G_op,w); subplot(211);semilogx(w,20*log10([m(:)’;m1(:)’])) subplot(212);semilogx(w,[p(:)’;p1(:)’]) G_cl1=feedback(G,1);G_cl2=feedback(G_op,1); [y,t]=step(G_cl1);y=[y,step(G_cl2,t)];figue,plot(x,y)  如果开环系统的响应不很理想，则可以引入超前补偿器，这种补偿器可以增加系统的增益，从而增加系统的剪切频率。同时还会提高该频率下的相位值，使得系统的相位裕量也增加。 补偿后闭环系统的阶跃响应曲线如图3－24(b)所示。可以看出，这样补偿下系统的阶跃向应特性显著地改进了。因为随着系统相位裕量的增加，超调量减小了，而随着系统剪切频率的增加，系统的响应速度也加快了。 2．串联超前滞后补偿器设计 利用MATLAB提供的强大功能，可以方便地选择一个串联补偿器，并立即得到校正后的系统性能指标，通过适当的调整就可以得到满足设计要求的补偿器。然而，如何通过解析的方法，直接得到符合设计要求的补偿呢？首先，设所求串联超前/滞后校正环节的传递函数为  (≤,≤) (3-28) 如果对象模型在剪切频率ω=ωc处的相位角为，则在该频率下补偿器的相位角应该为其中γ为期望的补偿后系统相位裕量。记ω=ωc频率下对象模型的幅值为，则可以采用下面的步骤对系统进行校正。 当＞0时，需要超前补偿器，该控制器可以设计成  (3-29) 且,  (3-30) 补偿后系统的稳态误差系数为  (3-31) 其中为对象模型G(s)在s= O极点的重数，而G0(S)为带有补偿器系统的开环传递函数模型。根据系统稳态误差的设计要求，可以确定补偿后系统满足这一指标的增益系数K(容许稳态误差系数)。如果K1≥Kv，则对指定的相位裕量用超前补偿就足够了；否则，还应该再设计相位超前滞后补偿器。 ② 超前滞后补偿器可以进一步设计成 , (3-32) ③ 如果φc(ωc)＜0，则需要按下面的方法设计相位滞后补偿器。 ,, (3-33) 式中。 根据上面的算法，可以编写一个 MATLAB的 leadlag函数进行超前/ 滞后补偿器设计，该函数的调用格式为： [Gc,Kc,Zc1,Pc1,Zc2,Pc2]=leadlag(G,Wce,gama,Kv,key) 或 [Gc]=leadlag(G,Wce,gama,Kv,key) 程序如下： function [Gc,Kc,Zc1,Pc1,Zc2,Pc2]=leadlag(G,Wce,gama,Kv,type) G=tf(G);[gain,phase]=bode(G,Wce); Phi_c=sin((gama-phase-180)*pi/180); den=G.den{1};a=den(length(den):-1:1); aa=find(abs(a)<=0) num=G.num{1}; G_n=num(length(num)); if length(aa)>0 if aa(1)>1,a=a(aa(1)+1); else,a=a(1):end alpha=sqrt((1-Phi_c)/(1+Phi_c)); Zc=alpha*Wce;Pc=Wce/alpha; Kc=sqrt((Wce*Wce+Pc*Pc)/(Wce*Wcd+Zc*Zc))/gain; K1=G_n*Kc*alpha/a; if nargin==4,type=1; if Phi_c<0,type=2; else,if K1<Kv,type=3;end,end end if nargout==1 switch type case 1, Gc=tf([1 Zc]*Kc,[1 Pc]); case 2,Kc=1/Gai;K1=G_n*Kc/a; Gc=tf([1 0.1*Wce],[1 K1*Gcn(2)/Kv]); case 3, Zc2=Wce*0.1;Pc2=K1*Zc2/Kv; Gcn=Kc*conv([1 Zc],[1 Zc2]); Gcd=conv([1 Pc],[1 Pc2]); end else switch type case 1, Zc1=Zc;Pc1=Pc; case 2,Kc=1/Gai;K1=G_n*Kc/a;Kc=K1; Zc2=0.1*Wce;Pc2=K1*Gcn(2)/Kv; case3, Zc1=Zc;Pc1=Pc; Zc2=Wce*0.1;Pc2=K1*Zc2/Kv; end,end 其中G为对象模型的LTI对象，Wce为期望的剪切频率，gama为期望的校正后系统相位裕量． Kv为容许稳态误差系数。如果给出变量type，则将按 type的值选择控制器结构，否则将自动选择控制器结构。返回的变量为补偿器的传递函数模型Gc及补偿器参数。 例3-27设被控对象的传递函数  试设计串联补偿器，使系统速度稳态误差小于10% ，相角裕量PM=45° 。 解：绘制原系统的伯德图如图3-25(a) 所示，由图可以知道，系统相角裕量为32°，小于期望值45°,需要进行超前补偿。 由速度稳态误差可知系统容许误差系数Kv=10，设系统期望穿越频率保传不变,即ωc=9 rad/s。用上面的补偿器设计函数leadlag可以得到系统的超前补偿器。程序如下。 G=tf(400,[[1,30,200],0]);gamma=45;we=9; f1=figure;f2=figure;m=[];p=[];y=[]; w=logspace(-1,2);t=0:0.001:5; Gc=leadlag(G,we,gamma,10); G_o=Gc*G;G_c=feedback(G_o,1); [m1,p1]=bode(G_o,w); [m2,p2]=bode(G,w); figure(f1) subplot(211),semilogx(w,20*log10([m1(:)’;m2(:)’])) subplot(212),semilogx(w,[p1(:)’;p2(:)’]) G_c1=feedback(G_o,1); G_c2=feedback(G,1); [y,t]=step(G_c2);y=[y,step(G_c1,t)]; figure(f2),plot(t,y)  补偿前后系统伯德国和阶跃响应如图3－25(a)、(b)所示。超前补偿器为  补偿后系统 GM、PM、ωc可由margin函数得到。GM=3.461 dB,PM=44.999°, ωg =19.291 rad/s, ωc=9.000 rad/s ，这些指标均满足设计要求。 例3-28设被控对象的传递函数  试设计串联补偿器，使系统速度稳态误差小于2%，相角裕量PM= 40°,截止频率不小于2.3 rad/s。. 解:首先绘制原系统的伯德图及阶跃响应，得到系统相角裕量15.6°,小于期望值40°;系统响应速度较快，但有较大超调量和稳态误差。在这种情况下超前补偿不能解决稳态误差与暂态响应问题，需要进行滞后补偿。 由于速度稳态误差为2%，因此系统容许误差系数为儿Kv=50；系统期望穿越频率为ωc =2.3 rad/s. 用上面的补偿器设计函数leadlag可以得到合适的滞后补偿器。程序如下。 g=tf(30,[conv([0.1,1],[0.02,1]),0]); gamma=40;wc=2.3;Kv=50;type=2; f1=figure;f2=figure;m=[];p=[];y=[]; w=logspace(-2,2);t=0:0.001:5; Gc=leadlag(G,wc,gamma,Kv,type); Go=Gc*G;G_c=feedback(G_o,1); G_c1=feedback(G,1); [m1,p1]=bode(G_o,w); [m2,p2]＝bode(G,w); y1=step(G_c,t); y2=step(G_c1,t); figure(f1); subplot(211);semilogx(w2,20*log10([m1(:),m2(:)])) subplot(212);semilogx(w,[p1(:),p2(:)]) figure(f2),plot(t,[y1,y2]) 补偿前后系统伯德图和阶跃响应如图4－26(a)、(b)所示。滞后补偿器为  补偿后系统GM,PM,ωt可由 margin函数得到。GM＝24.728 dB, PM= 68.927°, ωt=2.310 rad/s，校正后指标均满足设计要求。  图3-26 例3-28的伯德图和阶跃响应 例3-29设被控对象的传递函数  试设计串联补偿器，使系统速度稳态误差小于 5%，相角裕量PM＝45士3°，动态调节时间不超过2.0 s，幅值裕量不低于10 dB。 解：由原系统的伯德图及阶跃响应可知，原系统GM、PM均小于零，系统闭环不稳定。由于原系统截止频率处的相角滞后远小平-180°,且响应速度有一定要求，因此考虑对系统进行串联超前/滞后校正。 由速度稳态误差5%得容许误差系数为Kv=20；设截止频率为20 rad/s，用上面的补偿器设计函数leadlag可以得到合适的超前滞后补偿器。程序如下。 G=tf(180,[conv([0.167],[0.5,1]),0]); gamma=48;wc=20;Kv=20;type=3; f1=figure;f2=figure;m=[];p=[];y=[]; w=logspace(-2,2);t=0:0.001:5; Gc=leadlag(G,wc,gamma,Kv,type); G_o=Gc*G;G_c=feedback(G_o,1); G_c1=feedback(G,1); [m1,p1]=bode(G_o,w); [m2,p2]＝bode(G,w); y1=step(G_c,t); y2=step(G_c1,t); figure(f1); subplot(211);semilogx(w2,20*log10([m1(:),m2(:)])) subplot(212);semilogx(w,[p1(:),p2(:)]) figure(f2),plot(t,[y1,y2])  图3－27 例3-29的伯德图和阶跃响应 计算得滞后补偿器为  补偿后系统 GM、PM可以margin函数得到，GM=27.398 dB，PM=90.741°从图3－27阶跃响应曲线可以看出，调节时间约为1.5 s，小于设计要求的2.0 s，因此系统稳定，且性能指标满足设计要求。 3．6．2比例积分微分(PID)控制 与串联超前/滞后补偿相比，PID控制算法被更广泛地作为调节器或控制器用于工程实际中。PID控制器可以采用以下几种形式：比例(P)控制器KD=KI=0；比例微分(PD)控制器，KI=0；比例积分(PI)控制器KD=0；全PID控制器。比例控制是最简单的控制结构，但它也能使系统满足某一方面的特性要求，如GM、PM、稳态误差等。在闭环控制中，加入微分控制，增加阻尼作用；加人积分控制，则增加了系统的稳态误差猜度(系统增型)。PID控制器可以用于补偿系统达到大多数特性参数的要求，是目前在过程控制中用得最普遍的控制器，被广 泛地用于各类控制(模拟、数字和适应控制)中。 1．Ziegler－Nichols方法 由于在设计PID控制器中，要调整3个参数，根轨迹与伯德图设计方法通常不被直接采用。Ziegler与Nichols发展了PID调节器设计方法，该方法基于简单的稳定性分析方法。首先，置KD=KI=0，然后增加比例系数直至系统开始振荡(即闭环系统极点在jω 轴上)。有将该比例系数乘0.6，而其他参数按下式计算得出 ,, 式中：Km为系统开始振荡时的K值；ωm为振荡频率。尽管该设计方法在设计过程中没有考虑任何特性要求，但是实践证明，Ziegler_ Nichols方法得到的参数，使过程控制器具有良好的工作性能。 利用根轨迹法或伯德图方法可以确定Km和ωm。例如，对于给定的被控对象的传递函数，可以得到一根轨迹图。对应穿越jω 轴时的增益Km，而此点的ω值即为ωm。另外，对于给定的被控对象传递函数可绘制伯德图，在频率为ωg处确定增益裕量Km=10(GM/20)，而ωn=ωg。注意：伯德图只能给出近似的结果。 例3-30被控对象的传递函数  用 Ziegler－NichoIs方法设计 PID控制器。、 解：计算得出其主要参数为： 闭环极点｛－4.2士j0.93，－21.5｝；增益穿 越频率ωc=-1.95 rad/s；GM=23 dB；PM=73°。 该闭环阶跃响应没有超调，而其斜坡输入稳； 态误差为 0.5 。对应该系统的根轨迹(图 3－28)， 使用rlocus及rlocfind命令可以求得增益Km =14.7和穿越频率ωm=14.0 rad/s，使用Ziegler- Nichols方程可求出了下例增益参数约为： ,, 程序如下 G=tf(400,[[1,30,200],0]); rlocus(G); [km,pole]=rlocfind(G); wm=imag(pole(2)); kp=0.6*km;kd=kp*pi/(4*wm);ki=kp*wm/pi; nk=[kd kp ki],dk=[1 0]; 图3－29绘制了系统补偿前后的开环伯德图与闭环阶跃响应。被补偿的阶跃响应具有60%超调量，过渡过程时间为1.2 s，增益穿越频率为11 rad/s，相角裕量为25°。程序如下。 Gc=tf([0.5,9,40],[1,0]);G=tf(400,[[1,30,200],0]); G_op=Gc*G; [m,p,w]=bode(G,{.01,100});[m1,p1]=bode(G_op,w); subplot(211);semilogx(w2,20*log10([m(:)’;m1(:)’])) subplot(212);semilogx(w,[p(:)’;p1(:)’]) G_cl1=feedback(G,1);G_cl2=feedback(G_op,1); [y,t]=step(G_cl1,3,5);y=[y,step(G_cl2,t)]; figure,plot(t,y) 在传统的Ziegler－Nichols方法之后出现了各种各样的变形及类似方法，如Chien－Hrones- Reswick(CHR)方法，Cohen－Coon公式等，如果读者有兴趣可以参考相关文献。  图3－29开环伯德图与闭环阶跃响应 2．解析方法 如前所述，采用Ziegler－Nichols方法不能设计一个PID控制器去满足某一特定闭环系统的特性要求。解析法则可以根据给定的稳态误差和操作特性参数来确定PID的参数，PID控制系统的开环传递函数为 , 如果 G(s)是n 型系统，补偿后的系统则为对n+1型系统。误差常数Kn+1；等于稳态误差ess的倒数，即 ∣ 对于给定的稳态误差指标，由上式可以求得KI的值。由时域指标，如超调量和过渡过程时间，可以确定闭环阻尼系数和自然振荡频率。该系统的传递函数为  式中：,。对于K值较大的情况，增益穿越频率给出了闭环自然频率的精确测量值。这样就得到了闭环自然振荡频率与开环增益穿越频率的关系。 闭环阻尼系数与相角裕量PM的关系可由奈奎斯特图中相角裕量PM的定义，并利用简单的几何方法得到，即  因为，则或。 因此，自然振荡频率对应开环增益穿越频率ωc，而希望的相角裕量PM可以由闭环阻尼系数求出。所以在处，补偿的系统增益应为 1，相角。由上述分析结果(且K由知)可以写出  又可以导出  由此可以看出， 上述过程可以编制成MATLAB的函数，其调用方法为 [kp,kd,nk,dk]= pid(G,hi,dpm,wgc) 程序如下： function [kp,kd,nk,dk]=pid(G,ki,dpm,wgc) ngv=polyval(ng,j*wgc); dgv=polyval(dg,j*wgc); g=ngv/dgv; thetar=(dpm-180)*pi/180; ejtheta=cos(thetar)+j*sin(thetar); eqn=(ejtheta/g)+j*(ki/wgc); x=imag(eqn); r=real(eqn); kp=r kd=x/wgc if ki~=0, dk=[1 0];nk=[kd kp ki]; else dk=1;nk=[kd kp]; end 例3-31设计PID控制的解析方法，被控传递函数为  系统的技术指标要求：单位斜坡输入稳态误差为0.1；超调量为10%；过渡过程时间为2 s。 解：因为原系统为Ⅰ型系统，可以求得加入PID控制器后的稳态误差常数为 ∣  由超调量与过渡过程两项技术指标，可以确定希望的闭环阻尼系数与自然振荡频率分别为  rad/s。因此，给定rad/s, 调用函数 pid可以得到  由图3-30可以看出，期望的系统相角裕量是80°，而增益穿越频率是4 rad/s。闭环系统的阶跃响应见图3-31。由图可见，阶跃响应的过渡过程时间大约为2 s，然而超调量大于技术指标要求，为22%，这是由于PID控制器引入的零点所造成的。绘图程序与图3-29相同。  图3－30例3-31的系统伯德图 图3-31例3-31的系统阶跃响应 3．PD控制 积分控制多用于系统增型，以减小系统稳态误差。积分控制的弱点是系统增加了极点，使得系统向不稳定的方面变化。相反地．系统增加微分环节，相当于系统增加了零点，使得系统向稳定方向变化。因此，微分控制器常添加在反馈通道中，特别是在希望增加阻尼比的情况下。PD控制器设计方法同样可以使用前面的pid函数实现参数选择，因此就不再举例详细讨论。