君，已阅读到文档的结尾了呢~~
扫扫二维码，随身浏览文档
手机或平板扫扫即可继续访问
Matlab仿真在自动控制理论中的应用
举报该文档为侵权文档。
举报该文档含有违规或不良信息。
反馈该文档无法正常浏览。
举报该文档为重复文档。
推荐理由：
将文档分享至：
分享完整地址
文档地址：
粘贴到BBS或博客
flash地址：
支持嵌入FLASH地址的网站使用
html代码：
&embed src='/DocinViewer--144.swf' width='100%' height='600' type=application/x-shockwave-flash ALLOWFULLSCREEN='true' ALLOWSCRIPTACCESS='always'&&/embed&
450px*300px480px*400px650px*490px
支持嵌入HTML代码的网站使用
您的内容已经提交成功
您所提交的内容需要审核后才能发布，请您等待！
3秒自动关闭窗口基于matlab-simulink的系统建模预防针技术与应用-第六章_图文_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
评价文档：
基于matlab-simulink的系统建模预防针技术与应用-第六章
上传于||文档简介
&&基​于​m​a​t​l​a​b​-​s​i​m​u​l​i​n​k​的​系​统​建​模​预​防​针​技​术​与​应​用​-​第​六​章
大小：8.60MB
登录百度文库，专享文档复制特权，财富值每天免费拿！
你可能喜欢《自动控制理论》课程设计指导书霍爱清 薛朝妹自动化教研室目录1控制系统的数学描述????????????????????????????1 1.1 1.2 1.3 1.4 微分方程??????????????????????????????1 传递函数??????????????????????????????6 状态空间描述????????????????????????????9 模型转换??????????????????????????
????11 单变量系统的两种主要校正方式????????????????????15 串联校正举例????????????????????????????18 概述????????????????????????????????24 控制系统函数全集??????????????????????????30 应用实例??????????????????????????????31 引导????????????????????????????????412控制系统的校正??????????????????????????????15 2.1 2.2 PI、PD、PID 校正??????????????????????????15 2.3 MATLAB 在自动控制系统中的应用?????????????????????24 3.1 3.2 3.334SIMULINK 简介??????????????????????????????41 4.1 4.2 SIMULINK 模型的构建????????????????????????455自动控制系统设计任务???????????????????????????50 5.1 5.2 5.3 5.4 任务一???????????????????????????????50 任务二???????????????????????????????52 任务三???????????????????????????????53 任务四???????????????????????????????54 时域分析例题及程序?????????????????????????56 根轨迹分析例题及程序????????????????????????64 频域分析例题及程序?????????????????????????666附录??????????????????????????????????56 6.1 6.2 6.3第一章 控制系统的数学描述1.1 微分方程1.1.1 物理系统的微分方程 利用机械学 、电学、流体力学和热力学等的物理规律，我们可以得到物理系统的动态 方程。它们通常用常系数线性微分方程来描述。 1.1.2 数值解 通过拉普拉斯变换和反变换，可得到线性时不变方程的解析解，也可用状态转移矩阵 φ (t)求解。这些分析方法通常只限于常系数的线性微分方程。解析解是精确的，然而通常寻 找解析解是困难的，甚至是不可能的。而数值分析方法直接在时域里求解微分方程，不仅适 用于线性时不变方程，也适用于非线性以及时变微分方程。 MATLAB 提供了两个求微分方程数值解的函数，它们采用龙格-库塔(Runge-kutta)法。 Ode23 和 ode45 分别表示采用 2 阶和 4 阶龙格-库塔公式，后者具有更高的精度。 n 阶微分方程必须化为 n 个首 1 的一阶微分方程组，且放入 M-文件中，以便返回方程 状态变量的导数，下面的例子介绍这些函数的用法。 例 1.1 对图 1-1 的机械系统，已知三个量――拉力、摩擦力、以及弹簧力都影响质量 M 的加速度。 利用牛顿运动定理，建立系统的力平衡方程式M令d2x dx ?B ? Kx ? f(t) d2t dtx 1 ? x, x 2 ?dx ，有 dtdx1 ? x2 dtdx 2 1 ? [f(t) ? Bx 2 ? Kx1 ] dt M设质量 M=1kg，摩擦系数 B=5N/m/sec，弹簧常数 K=25N/m。在 t=0 时刻，施加 25N 的 拉力。上述方程及已知量在 M-文件 mechsys.m 中定义如下： function xdot=mechsys(t, x); F=25;1M=1;B=5;K=25; xdot=[x(2);1/M*(F-B*x(2)-K*x(1))]; 下面的 M-文件使用 ode23 对系统在零初始条件下进行仿真： t0=0; x0=[0,0]; tol=0.001; trace=0; tfinal=3; ％时间间隔 0～3 秒 ％零初始条件 ％精度 ％如果非零，则打印出每一步的计算值[t, x]=ode23(’mechsys’,t0,tfinal,x0,tol,trace) subplot(211),plot(t, x); title (’Time response of mechanical translational system’) xlabel (’Time-sec’) text (2,1.2,’displacement’) text (2,.2,’veloclty’) d=x(:,1);v=x(:,2); subplot(212),plot(d,v); title (’velocity versus displacement’) xlabel (’displacement’) ylabel (’velocity’) subplot(111) 仿真结果如图 1-2。Time response of mechanical translational system 3 2 1 0 -1 displacement veloclty00.51.5 2 Time-sec velocity versus displacement12.533 2velocity1 0 -100.20.40.6 0.8 displacement11.21.4图 1-22例 1.2 电路图见图 1-3。R=1.4Ω ，L=2H，C=0.32F。初始状态：电感电流为零，电容 电压为 0.5V。ｔ＝０时刻接入１伏的电压。求 0&t&15 秒时 i(t)和 v(t)的值，并且画出电流与 电容电压的关系曲线。 根据基尔霍夫定律[RLC]有：Ri ? L和 令 得di ? v0 ? vs dtdv0 dti ?Cx1 ? v 0 , x 2 ? i? x1 ? 1 x2 C图 1-3? x2 ?1 (v s ? x 1 ? Rx 2 ) L以上方程式在 M-文件 electsys.m 中定义如下： function xdot=electsys(t, x); %状态导数 V=1; R=1.4; L=2; C=0.32; xdot=[x(2)/C;1/L*(V-x(1)-R*x(2))]; 下面的 M-文件使用函数 ode23 对系统进行仿真。 t0=0;tfinal=15; x0=[0.5,0]; tol=0.001; trace=0; ~时间间隔 0～15 秒 ％初值条件 x1=0.5,x2=0 ％精度 ％若非零值则打印出每一步的计算值 %阶跃输入[t, x]=ode23('electsys',t0,tfinal,x0,tol,trace); subplot(211), plot(t, x); title('Time response of a RLC series circuit') xlabel('Time-sec') text(8,1.05,'Capacitor voltage'),text(8,0.05,'current') vc=x(:,1), i=x(:,2); subplot(212),plot(vc, i); title('current versus capacitor voltage') xlabel('Capacitor Voltage') ylabel('current') subplot(111) 仿真结果见图 1-4。3Time response of a RLC series circuit 1.5 1 0.5 0 -0.5 current Capacitor voltage010 Time-sec current versus capacitor voltage5150.15 0.1current0.05 0 -0.05 -0.1 0.5 0.6 0.7 0.8 0.9 1 Capacitor Voltage 1.1 1.2 1.3图 1-4 1.1.3 非线性系统 在变量的一定范围内大多数物理系统是线性的。然而，随着变量范围的延伸，所有的系 统最终都是非线性的。对于非线性系统，叠加原理不再成立。Ode23 和 ode45 可以求解非线 性微分方程，见例 1-3。 例1.3 图 1-5 的单摆，一长Ｌ米的无重量绳悬挂重为 W=mg 千克 的物体,阻尼系数为 B 千克/米/秒。通常，该系统用线性微分方程近似描 述，实际上它是非线性的。 如果θ 为绳的倾角，物体的速度将为 Lθ ，那么离心力为? F ? ? Wsin? - BL?根据牛顿定理有F ? mL ?? ?联结上两式，得mL ?? ? BL ?? ? Wsin? ? 0 ? ?? 令 x 1 ? ?, x 2 ? ? 得4? x1 ? x 2 ,? x2 ? ?B W x2 ? sinx1 m mL％状态导数上面方程式在 M-文件 pendulum.m 中定义如下： function xdot=pendulum(t, x); W=2; L=.6; B=0.02; g=9.81; m=W/g; xdot=[x(2); -B/m*x(2)-W/(m*L)*sin(x(1))]; 下面的 M-文件对系统进行仿真。 t0=0; tfinal=5; ％时间间隔 0～5 秒 ％初始条件 x1=1,x2=0 trace=0; ％精度，不打印每一步的计算值x0=[1,0]; tol=0.0001;[t, x]=ode23('pendulum',t0,tfinal,x0,tol,trace); subplot(211),plot(t, x) title('Time response of a rigid pendulum'),xlabel('Time-sec') text(3.2,3.5,'Velocity'),text(3.2,1.2,'Angle-Rad') th=x(:,1);w=x(:,2); subplot(212),plot(th, w) title('Phase plane plot of pendulum') xlabel('Position-Rad'),ylabel('Angular velocity') subplot(111) 仿真结果见图 1-6。Time response of a rigid pendulum 4 2 Angle-Rad 0 -2 -4 Velocity00.511.52.5 3 3.5 Time-sec Phase plane plot of pendulum244.554Angular velocity2 0 -2 -4 -1-0.8-0.6-0.4-0.2 0 0.2 Position-Rad0.40.60.81图 1-651.1.4 线性化 在小信号条件下，非线性系统可以线性化。例１.３中单摆的运动方程在初始角较小的 情况下，可被线性化。设θ =θ +Δ θ ，单摆运动方程如下：? ? mL( ?? 0 ? ???) ? BL( ? 0 ? ??) ? Wsin(? 0 ? ??) ? 0 ? ?性方程(1.1)当Δ θ 〈 〈∞时，有 sinΔ θ ≈0,cosΔ θ ≈1,又因为θ 0 较小，有 sinθ 0≈Δ θ ，那么有线? mL ??? ? BL?? ? W ?? ? 0 ?当Δ θ 很小时，线性方程（1.2）的解与（1.1）一致。(1.2)1.2 传递函数线性时不变系统的传递函数定义为零初始条件下输出量的拉普拉斯变换与输入量的拉 普拉斯变换像函数之比。尽管传递函数只能用于线性系统，但它比微分方程提供更为直观的 信息。令传递函数的分母多项式等于零，便得到特征方程。特征方程的根是系统的极点，分 子多项式的根是系统的零点。那么传递函数便可由常数项与系统的零、极点确定。常数项， 通常记作 k，是系统的增益。利用传递函数，我们可以方便的研究系统参数的改变对响应的 影响。 通过拉普拉斯反变换可得到系统在时域的响应。 通常需要用有理函数的部分分式展开。 在这部分举几个例子介绍 MATLLAB 中求特征多项式的根，求传递函数的零、极点， 部分分式展开以及已知零、极点求传递函数等函数的功能。 1.2.1 多项式的根和特征多项式 如果 P 是包含多项式系数的行向量，roots(P)得到一个列向量，其元素为多项式的根。 如果 r 是包含多项式根的一个行/列向量，poly(r)得到一个行向量，其元素为多项式的系数。 例1.4 求多项式 s6+9s5+31.25s4+61.25s3+67.75s2+14.75s+15 的根。 多项式系数以降幂次序排列在一行向量中。用 roots 求根。 &&P=[1 9 31.25 61.25 67.75 14.75 15]; &&r=roots(P) 多项式的根从列向量 r 中得到 r= -4.0 -1.0000 + 2.0000i -1.0000 - 2.0000i -0.0000 + 5.0000i -0.0000 - 5.0000i 例 1.5 多项式的根为-1,-2,-3±j4。求多项式方程。6为了输入复数，必须首先建立虚数单位。然后在行/列向量中输入根。使用 poly 得到多 项式方程。 && i=sqrt(-1); && r=[-1;-2;-3+4*i;-3-4*i]; && p=poly(r) 多项式的系数从行向量中得到 p= 1 9 45 87 50 因此，多项式方程为 s4+9s3+45s2+87s+50=0 例 1.6 求下列矩阵的特征方程的根1 - 1? ?0 ?- 6 - 11 6 ? A?? ? ?- 6 - 11 5 ? ? ?用 ploy 求矩阵的特征方程的根。用 roots 求方程的根。 && A=[0 1 -1;-6 -11 6;-6 -11 5]; && P=poly(A) && r=roots(P) 结果如下： P= 1.0 11.0 r= -3.0 -1..2 传递函数的零点和极点 ? 函数 tf2zp 求传递函数的零点，极点和增益。 例 1.7 求下列传递函数的零点，极点和增益。H(s) ?&& num=[1 11 30 0];s 3 ? 11s 2 ? 30s s 4 ? 9s 3 ? 45s 2 ? 87s ? 50&& den=[1 9 45 87 50]； && [z, p, k]=tf2zp(num, den)7z= 0 -6.0 p= -3.0000 + 4.0000i -3.0000 - 4.0000i -2.0 因而有H(s) ??s(s? 5)(s ? 6) (s ? 1)(s ? 2)(s ? 3 ? j 4)(s ? 3 ? j 4)函数 zp2tf 根据给定零点，极点和增益求传递函数。 系统的零点为-6,-5,0,极点为-3±j4,-2,-1,增益为 1。求其传递函数。例1.8&& z=[-6;-5;0];k=1; && i=sqrt(-1); && p=[-3+4*i;-3-4*i;-2;-1]; && [num, den]=zp2tf(z, p, k) 上面程序的结果为 num = 0 den = 1 9 45 87 50 因此，传递函数为 1 11 30 0H(s) ?1.2.3 部分分式展开s 3 ? 11 s 2 ? 30 s s 4 ? 9 s 3 ? 45 s 2 ? 87 s? 50函数[r, p, k]=residue(b, a)，对两个多项式的比进行部分分式展开，如P(s) b m s m ? b m ?1s m ?1 ? ? ? ? b1s ? b 0 ? Q(s) a n s n ? a n ?1s n ?1 ? ? ? ? ? a 1s ? a(1.3)向量 b, a 以 s 的降幂顺序排列多项式的系数。部分分式展开后余数送入列向量 r，极点 送入列向量 p，常数项送入 k。 例1.9 对 F(s)进行部分分式展开82 s 3 ? 9 s? 1 F(s) ? 3 2 s ? s ? 4 s? 4&& b=[2 0 9 1]; && a=[1 1 4 4]; && [r, p, k]=residue(b, a) 结果如下： r= 0.0000 - 0.0 + 0.2500i -2.0000 p= -0.0000 + 2.0000i -0.0000 - 2.0000i -1.0000 k= 2 因而，部分分式展开为2?? 2 j 0.25 ? j 0.25 ?2 1 ? ? ? 2? ? 2 s? 1 s? j 2 s? j 2 s? 1 s ? 4函数[b, a]=residue(r, p, k)将部分分式转化为多项式比 P(s)/Q(s)。1.3状态空间描述集总参数的线性网络可用微分方程表示为? x(t) ? Ax(t) ? Bu(t)(1.4)该系统的一阶微分方程即为状态方程，X 是状态向量。状态空间方法易采用数字或模拟 计算机求解。 另外， 状态空间方法容易拓展到非线性系统。 状态方程可从 n 阶微分方程得到， 或者在系统模型中选用合适的状态变量直接写出。 1.3.1 将微分方程化成状态方程 设一个 n 阶线性系统由微分方程描述，我们讨论如何选择状态变量，使该系统由状态方 程描述。dn y d n ?1 y dy ? a n ?1 n ?1 ? ? ? ? ? a 1 ? a 0 y ? u(t) n dt dt dty(t),u(t)分别为系统的输出、输入。9(1.5)状态模型不是唯一的，它依赖于状态变量的选择。一个有用的选择方法是定义? x 1 ? y, x 2 ? y, x 3 ? ??,? ? ?, x n ? y (n ?1) y那么，y 的 n 阶导数可用各状态变量的线性组合来表示。 于是有? x1 ? x 2 ? x2 ? x3┇? x n ?1 ? x n ? x n ? ? a 0 x 1 ? a 1 x 2 ? ? ? ? ? a n ?1 x n ? u(t)写成矩阵形式为? ? x1 ? ? 0 ?x ? ? 0 ? ? 2 ? ? ? ? ??? ? ? ? ? ? ? x n ?1 ? ? 0 ? x n ? ?? a 0 ? ? ? ?1 0 0 ? a10 1 0 ? a2? 0]x0 ? ? x 1 ? ?0 ? ? 0 ? ? x 2 ? ?0 ? ?? ? ? ? ? ? ? ? ? ? ? ? u ?t ? ?? ? ? ? ? 1 ? ? x n ?1 ? ?0? ? ? a n ?1 ? ? x n ? ?1? ? ? ? ?? ?(1.7)输出方程为 y=[1 0 0 即(1.8)? 0 ? 0 ? A?? ? ? ? 0 ?? a 0 ?1 0 0 ? a10 1 0 ? a20 ? ? 0 ? ? ? ? ? 1 ? ? ? a n ?1 ? ? ??0 ? ?0 ? ? ? B ? ??? ? ? ?0 ? ?1 ? ? ?D=0(1.9)C ? ?1 0 0 ? 0?例1.10 给定系统的微分方程描述为 y(3)+2y(2)+3y(1)+4y=5u(t) 利用(1.9)即可写出其相应的一个状态空间描述为：1 0 ? ?0 ?0 A?? 0 1 ? ? ?? 4 ? 3 ? 2? ? ?1.3.2 矩阵的对角化?0 ? B ? ?0 ? ? ? ?5 ? ? ?C ? ?1 0 0?上述状态方程，称为能控标准型。 为什么要研究 A 矩阵的对角化问题呢？原因之一是将矩阵对角化后， 对角元即为矩阵 A10的特征值λ 1,λ 2,?,λ n。从而，状态转移矩阵亦可化为对角元为 e （假定 A 有 n 个互不相同的特征值） 。λ1t,eλ2t,?,eλnt的对角矩阵。给定线性系统。当 A 有 n 个互不相同的特征值时，我们可以找到非奇异变换矩阵 P，令 x(t)=py(t) 将上面状态方程化为对角线规范型 (1.10)? y ? ay? bu(t)其中 a=p-1AP, b=p-1B 有很多方法可以求得 P，如利用 A 的特征向量可构造 P。 例1.11 定系统的状态空间描述为(1.11)? 1 - 1? ? x 1 ? ?0? ? x 1 ? ?0 ? x ? ? ?- 6 - 11 6 ? ? x ? ? ?0? u ? ? 2? ? ?? 2 ? ? ? ? x 3 ? ?- 6 - 11 5 ? ? x 3 ? ?1 ? ?? ? ? ? ?? ? ?求它的对角规范型和变换矩阵 P。 && A=[0 1 -1;-6 -11 6;-6 -11 5]; && [P,L]=eig(A); B=[0;0;1] ％L 是一个对角元为特征值的对角矩阵 ％P 是一个变换矩阵，其列是相应于特征值的特征向量 && a=inv(P)*A*P && b=inv(P)*B 结果为： p= -0.0 -0.7071 a= -1.0 -0.0000 b= 2.7 10.0 -2.0 -0.0 -3.2 0.9 -0.3 -0.8285 ％A 矩阵对角化1.4模型转换1.4.1 传递函数向状态空间描述的转换11控制系统工具箱包含一组模型转换的函数。 [A,B,C,D]=tf2ss(num, den)将传递函数转换为 状态空间描述。 例 1.12 求下面传递函数的状态空间描述&& num=[1 7 2]; den=[1 状态方程各矩阵如下：9 2624];&& [A, B, C, D]=tf2ss(num, den)?? 9 ? 26 ? 24 ? A??1 0 0 ? ? ? ?0 1 0 ? ? ??1 ? B ? ?0 ? ? ? ?0 ? ? ?已知状态方程和输出方程C ? ?1 7 2?D=01.4.2 状态空间描述向传递函数的转换? x ? Ax? Buy=Cx+Du 采用拉普拉斯变换有 Y(s)=C(sI-A)-1Bu(s)+Du(s) 则(1.12) (1.13)G(s) ?Y(s) ? C(sI? A)?1 B? D U(s)(1.14)函数 ss2tf(A,B,C,D,i)是将状态空间描述转换为对第一个输入的传递函数。 [num,den]=ss2tf(A,B,C,D,i)是将状态空间描述化为分子、分母多项式形式的传递函数。 [z,p,k]=ss2zp(A,B,C,D,i)将状态空间描述化为零极点形式表示的传递函数。 例 1.13 一个系统的状态空间描述如下? ? x 1 ? ?0 1 0 ? ? x 1 ? ?10 ? ? x ? ? ?0 0 1 ? ? x ? ? ?0 ? u ? ? 2? ? ?? 2 ? ? ? ? x 3 ? ?- 1 - 2 - 3? ? x 3 ? ?0 ? ?? ? ? ? ?? ? ?y=[1 0 0]x 求传递函数 G(s)=Y(s)/U(s)12&&A=[0 1 0; 0 00 1; -1 -2 &&C=[1 0 0];D=[0]; &&[num,den]=ss2tf(A,B,C,D,1) &&[z,p,k]=ss2zp(A,B,C,D,1)-3];B=[10; 0; 0];其中，ss2tf(A,B,C,D,1)中“1”表示对第一个输入。 传递函数的分子、分母多项式系数如下： num= 0 den= 1.0 2.0 传递函数的零、极点如下： z= -1 -2 p= -0.3i -0.3i -2.3247 k=10 因而传递函数为 10.0 20.000010 (s 2 ? 3 s ? 2) G(s) ? 3 s ? 3 s 2 ? 2 s? 1 10 (s ? 1)(s ? 2) ? (s ? 0.3376 ? 0.5623 i)(s ? 0.3376 ? 0.5623 i)(s ? 2.3247 )1.4.3 由方框图求状态空间描述和传递函数 控制系统工具箱中提供了函数[A,B,C,D]=connect(a, b, c, q, iu, iy)。将方框图描述转换成 状态空间描述和传递函数。其中 q 矩阵规定了各框之间的连接关系。其每一行的第一个元素 是框号，其余的元素依次是于该框连接的框号，iu,iy 分别表示输入，输出施加的框号。 例 1.14 将图 1-7 由框图表示的系统转换成状态空间描述和传递函数。 &&n1=1; d1=1; n2=0.5; d2=1; n3=4; d3=[1 4]; &&n4=1; d4=[1 2]; n5=1; d5=[1 3]; n6=2; d6=1; &&n7=5; d7=1; n8=1; d8=1; &&nblocks=8; blkbuild &&q=[1 2 0 0 0 -7 0 -8 ％q 矩阵表示框图的结构。 如第 2 个框于第 1 个框按131 -63 4 5 6 7 82 3 4 3 4 50 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0];1 的关系连接，于第 6.7.8 个框按-1 关系连接，依次类推。&&iu=[1]; &&iy=[8]; &&[A, B, C, D]=connect(a, b, c, d, q, iu, iy) &&[num, den]=ss2tf(A,B,C,D,1) 结果为 A= -8.0 0.4 0 B= 0.5 0 0 C= 0 0 1 D= 0 num= 0 den= 1.0 即 13.0 56.0 80.0 0 0 2 -2.5 -2.0 1.0 -0.5 0 -3.0％输入施加于第 1 个框上 ％由第 8 个框输出 ％转换成传递函数C(s) 2 ? 3 2 R(s) s ? 13 s ? 56 s? 8014第二章 控制系统的校正2.1 单变量系统的两种主要校正方式单变量系统常用的校正主要有两种方式。一种是校正装置与被控对象串联，如图 2-1 所 示。这种校正方式称为串联校正。另一种校正方式是从被控对象中引出反馈信号，与被控对 象或其一部分构成反馈回路，并在局部反馈回路设置校正装置。这种校正方式称为局部反馈 校正，如图 2-2 所示。图 2-1图 2-2 串联校正和局部反馈校正应用都相当普遍， 究竟选择哪一种， 取决于系统中信号的性质， 可供采用的元件以及其他条件。两种校正方式结合起来可以收到更好的效果。2.2PI、PD、PID 校正2.2.1 PD 校正（超前校正） 超前校正（亦称 PD 校正）的传递函数为G(s) ?? Ts? 1Ts ? 1, (? ? 1)（2.1）其对数频率特性如图 2-3 所示。超前校正能够产生相位超前角，它的强度可由参数α 表 征。 超前校正的相频特性函数是 θ （ω ）=arctgα ω T-arctgω T 最大相移点位于对数频率的中心点，即15（2.2）?m ?最大相移量为11 ? T ?(2.3)? m ? ? (? m ) ? arctg ? ? arctg或者1?? arcsin? ?1 ? ?1（2.4）sin ? m ?? ?1 ? ?1（2.5）??容易求出，在ω ｍ点有1 ? sin ? m 1 ? sin ? mL(ω ｍ)=10lgα(2.6)图 2-3 基于频率法综合超前校正的步骤是： 1. 首先根据静态指标要求，确定开环比例系数 K，并按已确定的 K 画出系统固有部分 根据静态指标要求预选ω c，从 Bode 图上求出系统固有部分在ω c 点的相角。 根据性能指标要求的相角裕量，确定在ω c 点是否需要提供相角超前量。如需要， 如果所需相角超前量不大于 60?,按（2.5）求出超前校正强度α 。 令 ? m ? ?c ? 1 的 Bode 图。 2. 3.算出需要提供的相角超前量θ m。 4. 5. 6. ，从而求出超前校正的两个转折频率 1/α T 和 1/T。( ? T)计算系统固有部分在ω c 点的增益 Lg(dB)及超前校正装置在ω c 点的增益 Lc(dB)。 如16果 Lg+Lc&0，则校正后系统的截止角频率ω c′比预选的值要高。如果高出较多，应采用滞 后超前校正，如果只是略高一些，则只需核算ω c′点的相角裕量，若满足要求，综合完毕， 否则转第 3 步。 如果 Lg+Lc&0， 则实际的ω c′低于预选的ω c。 可将系统的开环增益提高到 Lg+Lc=0 （即 将系统的开环比例系数提高 lg-1[-(Lg+Lc)]/20 倍） 。 超前校正的主要作用是产生超前相角，可用于补偿系统固有部分在截止角频率ω c 附近 的相角滞后，以提高系统的相角稳定裕量，改善系统的动态特性。 2.2.2 PI 校正（滞后校正） 滞后校正（亦称 PI 校正）的传递函数的G 0 (s) ?Ts? 1 , ( ? ? 1) ? Ts? 1(2.7)其对数频率特性如图 2-4 所示。参数β 表征滞后校正的强度。图 2-4 基于频率法的滞后校正指标的综合步骤是： 1. 首先根据静态指标要求确定开环比例系数 K， 按照所确定的 K 画出系统固有部分的 根据动态指标要求试选ω c，从图上求出在试选的ω c 点的相角，判断是否满足相位 Bode 图。 2. 裕量的要求（注意计入滞后校正将会带来的 50?~120?的滞后量） ，如果满足，转向下一步。 否则，如果允许降低ω c，就适当重选较低的ω c。 3.c)=20lgβ从图上求出系统固有部分在ω c 点的开环增益 Lg(ω c)。如果 Lg(ω c)&0，令 Lg(ω ，求出β ，就是滞后校正的强度，如果 Lg(ω c)&0，则无须校正，且可将开环比例 选择ω 2=1/T=（1/5~1/10）ω c，进而确定ω 1=1/(β T)。 画出校正后系统的 Bode 图，校核相位裕量。17系数提高。 4. 5.滞后校正的主要作用是降低中频段和高频段的开环增益， 但同时使低频段的开环增益不 受影响，从而达到坚固静态性能和稳定性。它的副作用是会在ω c 点产生一定的相角滞后。 2.2.3 PID 校正（滞后超前校正） 超前校正的主要作用是增加相角裕量，改善系统的动态响应特性。 滞后校正的主要作用 是改善系统的静态特性，两种校正结合起来就能同时改善系统的动态和静态特性特性。滞后 超前校正（亦称 PID 校正）综合了前面两种校正的功能。 滞后超前校正的传递函数为G 0 (s) ?(T2 s ? 1)(? T1s ? 1) , ( ? ? ? ? 1), (T2 ? ? T1 ) ( ? T2 s ? 1)(T1s ? 1)(2.8)它相当于一个滞后校正与一个超前校正相串联，其对数频率特性如图 2-5 所示。图 2-5 基于频率法的滞后超前校正的综合步骤是： 1. 首先根据静态指标要求确定开环比例系数，按照所确定的 K 画出系统固有部分的 按指标要求确定ω c，检查系统固有部分在ω c 的对数幅频的斜率是否为-2，如果是， 按综合超前校正的步骤 3~6 综合超前部分 Gc1(s)(注意在确定θ ｍ时要计入滞后校正 Bode 图。 2. 求出ω c 点的相角。 3. 将造成的 50~120 的相角滞后量)。在第 6 步时注意，通常 Lg(ω c)+Lc(ω c)比 0 高出很多，所 以要引进滞后校正。 4. 5. 6. 令 20lgβ = Lg(ω c)+Lc(ω c),求出β 。 按综合校正的步骤 4~5 综合滞后部分 Gc2(s)。 将滞后校正与超前校正串联在一起，构成滞后超前校正：Gc(s)=Gc1(s)*Gc2(s)。2.3串联校正举例对于一结构如图 2-6 所示的系统，给定固有部分的传递函数 Gg(s)和性能指标要求，试设计串联校正装置 K(s)。设18G g (s) ?1. 2. 3.100 s(0.1s? 1)(0.01 s? 1)若要求开环比例系数 K≥100 1/秒，相角裕量γ ≥30??，ω c≥45 1/秒。 若要求开环比例系数 K≥100 1/秒，相角裕量γ ≥40?，ω c=5 1/秒。 若要求开环比例系数 K≥100 1/秒，相角裕量γ ≥40?，ω c≥20 1/秒，分别设计校正装置，并比较它们的作用效果。图 2-6 对于 1，需进行超前校正，这里给出一个参考的 K(s)=(0.05s+1)/(0.005s+1)，于是系统的 开环传递函数为G 0 (s) ?闭环传递函数为0.05 s? 1 100
(s ? 20) ? ? 0.005 s? 1 s(0.1s? 1)(0.01 s? 1) s(s? 10)(s ? 100 )(s ? 200 )G 0 (s) ? s ? 20,000 ,000 s(s? 10)(s ? 100 )(s ? 200 ) ?
s ? 20.000 .000 s ? 20,000,000 ? 4 3 s ? 310s ? 23,000s 2 ? s ? 20,000,000其单位阶跃响应的实现为： num=[ den=[1 310 20,000,000]; 23,000 ,000,000];t=xx: xx:（起始时间，步长，终止时间） c=step(num, den, t); plot(t, c); xlabel (’t-sec.’), ylabel(’c(t)’),grid, pause 对于 2 需进行滞后校正，这里给出一参考的K(s) ?0.5 s? 1 10 s? 1仿照上述，可求出系统的阶跃响应。 对于 3，需进行滞后超前校正，这里给出一参考的K(s) ?(0.25 s? 1)(0.1s? 1) (1.25 s ? 1)(0.02 s ? 1)19同理可得到系统的阶跃响应。 比较上述三种不同的校正装置 K(s)所产生的不同效果，可以看出三种校正方式的作用。 这里需要特别指出阶跃响应的起始时间，步长及终止时间的给定问题。 起始时间通常给 0，终止时间要根据系统的过渡过程时间的长短来给定。我们可按估计 过渡过程时间的经验公式 ts≈(4~9)/ω c 给定此值：如 2，ω c=5 1/秒，则过渡过程时间大约在 0.8~1.8 秒之间，于是终止时间 可给 2 秒，我们希望在 0 至 2 秒内显示出 20 个点，则步长为 0.1 秒。若已知系统的ω c，建 议步长按 1/(2ω c)给定。对ω c 未知，也不便估计过渡过程时间的系统只好用试探法，待曲线 显示出以后再调整步长和终止时间。 图 2-7 为系统固有部分的频率特性图（幅频特性和相频特性）可以看出该系统已接近临 界稳定状态了。 图 2-8 至图 2-11 分别表示不加校正，加超前校正，加滞后校正，和加滞后超前校正的 系统阶跃响应。Bode Diagrams 100 50Phase (deg); Magnitude (dB)0 -50 -100 0-100-200-300 -1 10100101102103Frequency (rad/sec)图 2-7num=100; den=conv([1,0], conv([0.1,1], [0.01,1])); bode(num, den)202 1.8 1.6 1.4 1.2c(t)1 0.8 0.6 0.4 0.2 000.20.40.60.81 t-sec1.21.41.61.82图 2-8numo=100; deno=conv([1,0],conv([0.1,1],[0.01,1])); t=0:0.01:2; numc= denc=[zeros(1,length(deno)-length(numo)),numo]+ c=step(numc,denc,t); plot(t,c) xlabel('t-sec'),ylabel('c(t)'),grid,pause1.41.210.8c(t)0.6 0.4 0.2 0 00.20.40.60.81 t-sec1.21.41.61.82图 2-921numo=conv(100,[0.05,1]); deno=conv([0.005,1],conv([1,0],conv([0.1,1],[0.01,1]))); t=0:0.01:2; numc= denc=[zeros(1,length(deno)-length(numo)),numo]+ c=step(numc,denc,t); plot(t,c) xlabel('t-sec'),ylabel('c(t)'),grid,pause1.41.210.8c(t)0.60.40.2000.20.40.60.81 t-sec1.21.41.61.82图 2-10numo=conv(100,[0.5,1]); deno=conv([10,1],conv([1,0],conv([0.1,1],[0.01,1]))); t=0:0.01:2; numc= denc=[zeros(1,length(deno)-length(numo)),numo]+ c=step(numc,denc,t); plot(t,c) xlabel('t-sec'),ylabel('c(t)'),grid,pause221.41.210.8c(t)0.60.40.2000.20.40.60.81 t-sec1.21.41.61.82图 2-11numo=conv(conv([0.25,1],[0.1,1]),[100]); deno=conv(conv([1.25,1],[0.02,1]),conv([1,0],conv([0.1,1],[0.01,1]))); t=0:0.01:2; numc= denc=[zeros(1,length(deno)-length(numo)),numo]+ c=step(numc,denc,t); plot(t,c) xlabel('t-sec'),ylabel('c(t)'),grid,pause23第三章 MATLAB 在控制理论中的应用3.1 概述MATLAB 提供了大量的控制工程计算、设计库函数。其中，控制系统软件包包括复数 运算、特征值计算、方程求解、矩阵变换以及 FFT 等重要计算工具及举例。MATLAB 的线 性代数处理， 矩阵运算和数值分析的能力为控制系统工程设计及其它学科研究提供了可靠的 基础和强有力的研究工具。 控制系统软件包利用 MATLAB 矩阵功能提供了适用于控制工程的专用函数，这些函数 大部分用 M 文件表示。控制系统软件包可以方便地用于控制系统设计、分析和建模。 在控制系统软件包中，控制系统通常采用传递函数与状态空间两种形式建模，允许“经 典”和“现代”技术并用，既可处理连续时间系统也可处理离散时间系统，并且可以进行不 同模型表示形式之间的相互转换，也可以计算和绘制时间响应、频率响应及根轨迹图。此外 M 文件还能够进行极点配置和最优控制器的参数计算。 即使在软件包中没有提供的功能， 也 可以通过编写新的 M 文件方式来构造。 3.1.1 系统模型 控制系统软件包可用于线性时不变（简称 LTI）系统模型。时不变系统模型包括：状态 空间模型；传递函数模型；零一极点增益模型；部分分式模型；离散时间模型。 下面介绍如何用矩阵表示不同类型的线性时不变系统模型。 〔1〕状态空间模型 一个 LTI 微分方程系统总可以用一组一阶微分方程组来表示。 按矩阵或状态空间表示形 式，LTI 系统模型的一般形式如下? x =Ax＋Buy=Cx＋Du 其中：u 是一个 nu 维控制输入向量，x 是一个 ns 维状态向量，y 是 ny 维输出向量。 采用 MATLAB 表示状态空间系统十分容易。A，B，C，D 都是矩阵，均作为独立变量 处理。 例 设有一个具有两个极点的二阶系统，自然角频率ω n=1.5，阻尼系数ξ =0.2。按下面 Wn=1.5； z= 0.2； a=［0 124状态空间模型表示形式输入该系统-Wn^2 b=［0 Wn^2]； c=[1 d=0； 0]；-2*z*Wn]；在 MATLAB 中，采用状态空间表示系统模型是最通用的 LTI 系统模型的表示形式。对 于多输入多输出（简称 MIMO）系统，状态空间表示是唯一容易处理的模型表示形式之一。 （2）传递函数模型 状态空间系统模型的等效表示是拉氏变换的传递函数形式如下 Y（S）=H（S）U（S） 其中：H（s）=C（SI－A）-1B＋D 在通常情况下，H（S）需要用三维矩阵表示，即 H（S）的维数是 ny 行，nu 列，ns＋1 阶， 其中 ns＋1 是多项式的系数。由于 MATLAB 变量都是两维的，因此，只限于用单输入（输 入是一维的 u）多输出（简称为 SIMO）系统描述方式表示该系统H(s) ?N(s) N(1)s nn ?1 ? N(2)s nn ? 2 ? ? ? N(nn) ? q(s) q(1)s nq?1 ? q(2)s nq? 2 ? ? ? q(nq)行向量 q 是按 s 降幂排列的分母系数。 矩阵 N 包含了分子系数， 输出向量的个数与向量 行数相等。 例 考虑单输入多输出（简称 SIMO）系统H(S) ??33s ? 5s ? 2s ? 13S? 2 S3 ? 2S? 5 2?按下面形式输入此系统到 MATLAB 中 num=［ 0 0 3 2 1 0 2 5]； den=［ 3 5 2 1］ ； 系统输入表明，多项式中的缺省项系数值用零替补。 （3）零-极点增益模型 传递函数可以表示为系数或零一极点增益形式，SIMO 系统表示为H(s) ?Z(s) (s ? Z(1))(s ? Z(2)) ? (s ? Z(n )) ?K P(s) (s ? P(1))(s ? P(2)) ? (s ? P(n ))MATLAB 中规定，把多项式的根存放到列向量中，而行向量包含多项式系数。所以在 因子形式下， 列向量 P 包含有传递函数分母的极点分布， 分子零点存放到矩阵 Z 的列向量中。 矩阵的列数与向量 y 的输出数相等。分子传递函数的增益存放到列向量 k 中、对于 SIMO 系25统，k 是标量。 在 MATLAB 中，把多项式的根存储于列向量中，而多项式的系数存储于行向量中，两 者可以互相转换。 例 if P= ［l 3 5 2] then r=rOOts（P） r= -1.1i -1.1i -0.5436 and PP=poly（r） PP= l.0 5.0 使其又返回到原多项式。以后将会看到，对于 SIMO 的 LTI 系统，tf2zp 将多项式 传递函数形式转换为因子零-极点增益形式，而 zp2tf 则进行相反形式的转换。 例 考虑一个 SIMO 系统Z(s) 4 ( S?1)( S? 2 ) H(S) ? ? P(s) (s ? 3)(s ? 4)(s ? 5)输入形式如下 k=[3；4]； z=[-12 -1 inf -2]； P=［-3 -4 -5]; 在这个多输出系统中，第 1 个分子比第 2 个分子阶数低，因此，采用 inf 来补充第 1 个 分子的那个无限零点，使第 1 个分子与第 2 个分子有相同的阶数。 4）部分分式模型 传递函数也可以采用部分分式展开式或留数的形式表示。对于 SIMO 系统表示为?3( S?12)?H(s) ?r (1) r (2) r (n ) ? ??? ? k (s) s ? p(1) s ? p(2) s ? p( n )行向量 P 表示极点，行向量 r 表示对应于 P 中极点的留数，而且行向量 k 包含原传递函数中26分子阶次大于分母阶次的部分。 传递函数形式与部分分式形式可以互相转换。 5）离散时间系统模型 在 MATLAB 中，离散时间线性时不变系统与连续时间线性时不变系统有相同的表示方 法，即状态空间法，多项式传递函数法，零-极点增益法与部分分式法。 离散时间 LTI 系统模型可以表示为一组一阶差分方程组形式。 在状态空间或矩阵表示形 式中可以写成如下形式 x（n＋l）=Ax(n)+Bu(n) y（n）=Cx(n）＋Du(n) 其中：u 是控制输入向量，x 是状态向量, y 是输出向量。 另一种等效表示是 Z 变换传递函数描述法，其形式如下 Y（z）=H（z）U（z） H（z）=C（zI 一 A） ?1 B＋D 在单输入为 u 的 SIMO 系统中，Z 变换传递函数描述形式为H(z) ?N(z) N(1) ? N(2)z ?1 ? ? ? N(nn)z ? ( nn ?1) ? q(z) q (1) ? q (2)z ?1 ? ? ? q (nq)z ?( nq?1)其中：向量 q 中包含按单位延迟 1/z 升幕排列的分母系数，矩阵 N 包含分子的系数，其行数 与 y 输出的数目相等。 零一极点增益表示形式为H(z) ?3.1.2 模型转换Z(z) (z ?1 ? Z(1))( z ?1 ? Z(2)) ? (z ?1 ? Z(n )) ? K ?1 P( z ) (z ? P(1))( z ?1 ? P(2)) ? (z ?1 ? P(n ))下面一组函数允许 LTI 系统模型的不同表示形式之间可以互相转换。 （1）状态空间模型到传递函数模型的转换。 ［num，den]=ss2tf（a，b，c，d，u） （2）状态空间模型到零-极点增益模型的转换。 [z，p，k]= sstzp（a，b，c，d，u） (3)传递函数模型到状态空间模型的转换。 ［a，b，c，d］= tf2ss（num，den） （4）传递函数模型到零-极点增益模型的转换。 ［z，p，k］= tf2zp（num，den） （5）零一极点增益模型到状态空间模型的转换。 [a，b，c，d]= zp2ss（z，p，k） (6) 零-极点增益模型到传递函数模型的转换。 [num，den]= zp2tf（z，p，k）27（7）传递函数模型到部分分式模型的转换。 ［z，p，k］= residue（num，den） （8）部分分式模型到传递函数模型的转换。 [num，den]= residue（z，p，k） （9）连续系统模型到离散系统模型的转换。 [ad，bd]=c2d（a，b，Ts） （10）离散系统模型到连续系统模型的转换。 [a，b]=d2c（ab，bd，Ts） 3.1.3 分析函数 控制系统软件包提供了控制系统工程需要的基本的时域与频域分析工具函数。 连续时间系统分析函数 impulse 脉冲响应 step 阶跃响应 lsim 任意输入的仿真 bode 波特图 nyquist 奈奎斯特图 lyap 李雅普诺夫方程 gram 可控性与可观性 离散时间系统分析函数 dimpulse 单位采样响应 dstep 阶跃响应 filter SISO 系统 z 变换仿真 dbode 离散波特图 freqz SISO 系统 Z 变换频域响应 dlyap 李雅普诺夫方程 dgram 离散可控性与可观性 3.1.4 闭环系统建模 上节给出的函数为连续系统和离散系统提供了频域和时域分析工具，适合于 4 参数 (A,B,C,D)系统、3 参数（z，p，k）系统和 2 参数（N，q） 。但是还没有用于分析闭环系统的 专门命令工具。闭环系统必须满足完全的闭环系统动力学。 例 设状态空间描述的开环系统模型如下:? x =Ax＋Buy=Cx＋Du 并具有参考输入 r 的全状态反馈控制准则 u= -Kx＋Nr 为了给该闭环系统建模，求解该闭环系统矩阵? x = Ax＋Bu = Ax-BKx＋BNr =（A-BK）x＋BNr28y= Cx+Du = CX-DKx＋DNr =（C-DK）x＋DNr 组成闭环系统矩阵简单的 M 函数为 aa=a-b*k； bb= b*n； cc＝c-d*k； dd=d*n； 这个新函数允许用标准分析工具研究其闭环系统特性。 建立描述闭环系统（A，B，C，D）矩阵的这一基本方法可以扩展到建立更加复杂系统 的模型。对复杂系统来说， （A，B，C，D）矩阵伴随着子系统增加而急剧增大。由于采用 这种方法，建模和分析工具完全通用，可适用于任意 LTI 系统。 建立模型的其它函数 append connect parallel series minreal 3.1.5 设计函数 为了参考闭环控制系统选择反馈增益的过程，采用了控制系统设计（design）术语。设 计也包括控制器结构选择和可能性估计器结构。大部分设计方法是反复的，带有分析的组合 参数选择、仿真及物理观察。 控制系统软件包有一套帮助实现增益选择工具的函数。 对于这些函数更多的信息可以通 过在线帮助进一步了解。 增益选择函数 margin place rlocus lqe lqr dlqe dlqr 增益与相位裕量 极点配置 根轨迹 线性平方估计器设计 线性平方调节器设计 离散线性平方估计器设计 离散线性平方调节器设计 两个子系统构合成 方框图建模 系统并联连接后的等效系统生成 系统串联连后的等效系统生成 最小实现和零-极点相消特别是 connect 函数，是一种寻找状态空间模型的综合性函数。293．2 控制系统函数全集本节包含了控制系统软件包中函数的全集， 分类列出， 以便快速查阅。 对全部函数功能、 格式的详细描述可以利用在线帮助功能得到。 3.2.1 模型建立 append connect parallel series ord2 3.2.2 模型转换 ss2tf ss2zp tf2ss zp2tf zp2ss residue c2d d2c tf2zp 3.2.3 模型实现 ctrbf obsvf mineral balreal modred dbalreal dmodreal 3.2.4 模型特性 damp gram gramians dgram ctrb obsv tzero 3.2.5 时间响应 impulse step 两系统合成函数 方框图建模函数 系统并联后的等效系统生成函数 系统串联后的等效系统生成函数 形成二阶系统的 A, B，C，D 函数 状态空间模型到传递函数模型的转换函数 状态空间模型到零一极点模型的转换函数 传递函数模型到状态空间模型的转换函数 零一极点模型到传递函数模型的转换函数 零一极点模型到状态空间模型的转换函数 部分分式展开函数 连续时间模型到离散时间模型的转换函数 离散时间模型到连续时间模型的转换函数 传递函数模型到零一极点模型的转换函数 可控性阶梯形式函数 可观性阶梯形式函数 最小实现及零一极点相消函数 平衡实现函数 模型降价函数 离散平衡实现函数 离散模型降阶函数 阻尼系数及自然频率函数 可控性与可观性函数 用于时变系统的可控性与可观性函数 离散系统可控性与可观性 可控性矩阵函数 可观性矩阵函数 传输零点函数 冲击响应 阶跃响应30lsim dimpulse dstep dlsim filter 3.2.6 频率响应 bode nyquist dbode freqz freqs 3.2.7 增益选择 lqr lqe dlqr dlqe margin place rlocus 3.2.8 应用 lyap dlyap fixphase abcdcheck nargcheck任意输入的连续系统仿真 离散时间单位脉冲响应 离散时间阶跃响应 任意输入的离散时间系统仿真 SIMO Z 变换仿真 波特图 奈奎斯特图 离散波特图 Z 变换频率响应 拉氏变换频率响应 线性二次调节器设计 线性二次估测器设计 离散线性二次调节器设计 离散线性二次估测器设计 幅值和相角裕量 极点配置 根轨迹 李雅鲁诺夫方程 离散李雅鲁诺夫方程 波特图展开相角 检查（A，B,C，D）的一致性 检查 M 文件幅角数3.3应用实例3.3.1 简单系统示例 控制系统如图 3-1 所示，试分析其闭环特性。下面是该系统的演示过程，带符号％语句为注释语句。 %清文字屏 %clc %清图形屏31%clf %系统 1 描述 a1=[1 3 0]; b1=2; %系统 2 描述 a2=[1 2 2]; b2=[1 2]; %两系统串联后的传递函数： [b12,a12]=series(b1,a1,b2,a2); %串联后频率响应（波特图） figure(1) bode(b12,a12);pause %clf %求振幅及相位裕度 figure(2) [m,p,w]=bode(b12,a12) %求开环振幅 Gm 和相位裕度 Pm 及其对应频率 margin(m,p,w); [Gm,Pm,Wg,Wp]=margin(m,p,w);pause %画 nichols 曲线，先用 ngrid 函数画坐标网 figure(3) ngrid('new'); nichols(b12,a12);pause %clf %画根轨迹，先消除画 nichols 曲线时冻结的坐标系， %用 axis 命令恢复自动定标。 figure(4) rlocus(b12,a12); %用鼠标器找适当的极点 [k1,p1]= rlocfind(b12,a12);pause %将环路闭合 b1=k1*b1; b121=k1*b12; %再看闭环振幅及相位裕度32figure(5) [m,p,w]=bode(b121,a12); margin(m,p,w); [Gm,Pm,Wg,Wp]=margin(m,p,w);pause %求闭环传递函数，[ bb，ab]指单位反馈情况， %[ bb1，ab1]指将系统 2 置于反馈通道上： [bb,ab]=cloop(b121,a12);pause [bb1,ab1]=feedback(b1,a1,b2,a2,-1);pause %求闭环频率响应 figure(6) bode(bb1,ab1);pause % %闭环过渡过程 figure(7) step(bb1,ab1);pause %clg figure(8) impulse(bb1,ab1);pause % %再求闭环零一极点，作为校核： [zb,pb,kb]=tf2zp(bb,ab);pause printsys(zb,pb,kb); format echo off %演示结束 运行结果： 开环系统幅值裕量及相位裕量：Gm=3.5435,Pm=55.022,Wg=1.6159,Wp=0.66898Bode DiagramB o d e D ia g r a m 50Gm = 10.989 dB (at 1.6159 rad/sec), Pm = 55.022 deg (at 0.66898 rad/sec) 5000Magnitude (dB) Phase (deg)-1 0 1 2M a g n itu d e ( d B )-50-50-100-100-150 -90-150 -90-135-135Ph a s e ( d e g )-180-180-225-225-270 10 10 10 F r e q u e n c y ( r a d /s e c ) 10-270 -1 1010010 Frequency (rad/sec)110233Nichols Chart 40 0 dB 0.25 dB 0.5 dB 1 dB 3 dB 6 dB 0Root Locus 320-1 dB -3 dB -6 dB -12 dB21-20Open-Loop Gain (dB)-20 dBImag Axis-40-40 dB0-60-60 dB-1-80-80 dB-2-100 -100 dB-3-120 -360 -120 dB -315 -270 -225 -180 Open-Loop Phase (deg) -135 -90 -45 0-5-4-3-2 Real Axis-101在根轨迹上选择一点：selected_point = -4.0158 + 0.0309i 此时闭环系统幅值裕量及相位裕量：Gm= 0.34677,Pm= -29.162,Wg= 1.6159,Wp= 2.5192Bode Diagram Gm = -9.1992 dB (at 1.6159 rad/sec), Pm = -29.162 deg (at 2.5192 rad/sec) 50 20 10 0Magnitude (dB) Magnitude (dB)Bode Diagram0-10 -20 -30 -40 -50-50-100 -90-60 225-135Phase (deg) Phase (deg)0 1 2180135-18090-22545-270 -1 100 10 10 Frequency (rad/sec) 10 10-110010 Frequency (rad/sec)1102Step Response 10 50Imp u l s e Re s p o n s e40 5 3020 0AmplitudeA mp l i t u d e10-5 0-10 -10 -20-150123 Time (sec)456-30 0123 Ti me ( s e c )456根据选定根轨迹上一点的位置可计算出相应的闭环零点，闭环极点及闭环增益值。 3.3.2 综合设计示例 该设计举例是控制系统工具箱中的一个演示程序 diskdemo.m 的译文。它是通过计算机 硬盘读写磁头位置控制器的设计来说明设计经典数字控制系统的方法。 MATLAB 环境下， 在 键入 diskdemo 可以看到它的运行结果。对照本程序可以加深理解如何综合利用各种 MATLAB 的函数完成特定的系统设计任务。 设计程序如下：34echo off %Diskdemo %echo on %pause ％接任意键继续 ％运行开始%根据牛顿定律，读写磁头最简单的数学模型可表示为如下的微分方程形式： % I*theta_ddot＋ C*theta_dot＋ K*theta= Ki*i %其中： %I 为读／写磁头的惯量 %C 为支架的粘滞阻尼系数 %K 为回动弹簧倔强系数 %Ki 为马达转矩常数 %theta_ddot 读／写磁头的角加速度 %theta_dot 读／写磁头的角速度 %theta 读／写磁头的位置 %i 输入电流 %按任意键继续 %进行拉普拉斯变换，得传递函数为： %令：I=0.0lKgm^ 2，C＝0.004Nm／（rad／sec） ，K＝10Nm／rad %Ki＝0.05Nm／rad，则该系统的传递函数描述为： %I=0.01； C ＝0.004；K ＝ 10； Ki ＝0.05； I=0.01; C=0.004; K=10; Ki=0.05; num=[Ki]; den=[I,C,K]; %printsys(NUM,DEN,'s'); %pause %按任意键继续 %给定的任务是设计一数字控制器，用来精确控制读写磁头的位置。我们将在数字域中 %进行设计。由于被控的模型是连续的，我们首先必须将之离散化。MATLAB 中的函 %数 c2d%有数种方法进行离散化（例如' zoh'是零阶保持器，' foh'是一阶保持器等，可 %用 help c2dm 查询） ，并画出其 BODE 图。 %设采样时间 Ts = 0. 005（5ms） 。 Ts=0.005; figure(1)35[m,p,w]=bode(num,den); %margin(m,p,w); %[Gm,Pm,Wg,Wp]=margin(m,p,w); %[mag,phase,w]=bode(num,den);pause [num,den]=c2dm(num,den,Ts,'zoh'); %figure(2) %dbode(num,den,Ts); [mzoh,pzoh,wzoh]=dbode(num,den,Ts); %margin(mzoh,pzoh,wzoh); subplot(211) semilogx(w,20*log10(m)), semilogx(wzoh, 20*log10(mzoh)), xlabel('Frequency(rad/sec)'), ylabel('Gain db') title('c2d comparision plot') subplot(212) semilogx(w,p), hold on semilogx(wzoh,pzoh), hold off xlabel('Frequency(rad/sec)'), ylabel('Phase deg') %subplot(111), pause %按任意键继续? %现在分析离散时间系统 figure(2) dstep(num,den);pause %任意键结束绘图... %系统振荡颇大，可能是由于阻尼太小，可以通过计算开环特征值来检验这一点。 disp('Open loop discrete eigenvalues'); figure(3)36ddamp(den,Ts); zgrid('new'), pzmap(1,den); pause %按任意键结束绘图... %注意极点靠近单位圆，阻尼很小，需要设计一个补偿器以增大该系统的阻尼，我们 %尝试来设计一个补偿器。最简单的补偿器是加一个放大器构成负反馈，现在来看它 %的根轨迹： figure(4) rlocus(num,den);pause %hold off； %按任意键结束绘图... %如根轨迹所示，极点很快离开单位圆而使系统不稳定，可见需要引入一些超前网络 %或一些有零点的补偿器。引入下面的补偿器：D(Z)=k(z+a)/(z+b) ,其中：a&b %按任意键继续... %将设计补偿器与模型串联，并设：a=-0.85，b=0.0 % j＝sqrt（－ 1） ;则： [numc,denc]=zp2tf([0.85]',[0]',1); [num2,den2]=series(num,den,numc,denc); pause figure(5) %按任意键继续... %现在看频率响应的变化 [mag,phase,w]=dbode(num,den,1); %使用归一化频率 [mag2,phase2]=dbode(num2,den2,1,w); %按任意键继续 %现在绘制一幅比较图? echo off subplot(211) semilogx(w,20*log10(mag),w,20*log10(mag2)); xlabel('Frequency（rad／sec）'); ylabel('Gain db'); subplot(212);37semilogx(w,phase,w,phase2); xlabel('Frequency（rad／sec）'); ylabel('Phase deg'); %hold on %[s,limits]=inquire('axis'); %plot(limits(1:2),[-18O,-180],'w--'); %hold off； %绘-180°线 %按任意键继续? echo on %可以看出比较器确实引入了超前 pause %按任意键继续? %现在看加入补偿器之后的根轨迹。 figure(6) zgrid('new'), rlocus(num2,den2); hold off pause %按任意键结束绘图... %极点在单位圆中经过了一些区域。现在应使用 RLOCFIND 来选择具有最大阻尼系数的 %极点。 （ZGRID 可以 0.l 的步长画出 Z 从 0 到 1 之间的阻尼系数曲线） 。 pause %按任意键继续... %在曲线上选择一点。 [k,poles]=rlocfind(num2,den2); disp(num2str(k)); ddamp(poles,Ts); %生成闭环系统以便对设计进行分析。 [numc,denc]=feedback(num2,den2,k,1); %sys=feedback(num2,den2,k,1); %[numc,denc]=tfdata(sys,'v'); %这些特征值应与前面选择的特征值相一致。 disp('Closed loop eigenvalues'),38ddamp(denc,Ts); pause %按任意键继续... %绘制闭环时间响应曲线 dstep(numc,denc); pause %按任意键结束绘图... %可以看出时间响应相当理想，经过约 15 个采样点，即： %15 X TS 秒就稳定下来。 %disp(('Our disc drive will have a seek time'),num2str(15*Ts),'seconds') %pause %按任意健继续... %现在看一看系统的鲁棒性。鲁棒性最常见的经典指标是增益和相位裕量。该指标确 %定如下，变成单位反馈系统，计算 Bode 响应并寻找相位和增益的越界频率 Wc， %MATLAB 中的函数 MARGIN 可确定 Bode 响应的相位裕量和增益裕量。把设计的系统 %与前面选择的增益联结起来，以便组成一个单位反馈系统，首先计算开环 Bode 响应。 %[num2,den2]=series(num2,den2,1,1,Ts); %计算 Bode 响应和裕量。 [mag,phase,w]=dbode(num2,den2,Ts); [Gm,Pm,Wcg,Wcp]=margin(mag,phase,w); %绘制 Bode 曲线并显示裕量 figure(7) margin(mag,phase,w) pause %按任意键结束绘图.... %显示增益裕量（db） 、相位裕量和角频率。 %margins=[20*log10(Gm),Wcg,Pm,Wcp] %设计的系统有足够的鲁棒性，能够承受 10db 的增益增加和 40 度的相位滞后而保持稳 %定。继续这一设计过程，就能得到一个既能保持开环系统稳定又能节约寻找时间的 %补偿器（选用更大的阻尼可减少阶跃响应的稳定时间） 。 %echo off %程序结束。 运行结果：3940第四章 SIMULINK 动态仿真集成环境4.l 引导SIMULINK 是实现动态系统建模、仿真的一个集成环境。它的存在使 MATLAB 的功能 得到进一步的扩展。这种扩展的意义表现在： （1）实现了可视化建模。在 Windows 视窗里， 用户通过简单的鼠标操作就可建立起直观的系统模型，并进行仿真； （2）实现了多工作环境 间文件互用和数据交换，如 SIMULINK 与 MATLAB，SIMULINK 与 C、FORTRAN， SIMULINK 与 DSP，SIMULINK 与实时硬件工作环境等的信息交换都可以方便地实现； （3） 把理论研究和工程实现有机地结合在一起。 本章从 SIMULINK 工具包安装开始，先通过一个简单模型的创建示例向读者展现框图 模型创建、仿真的大致过程，进而对主要的操作做更详细的介绍。 4.1.1 系统要求 SIMULINK 的广泛应用是从 MATLAB 4．0 版开始的。在 4．0 版中，SIMULINK 被包 含在 MATLAB 的核心执行文件中；在 MATLAB 4．2 版中，SIMULINK 则以 MATLAB 的 工具包形式单独出现。 因此， SIMULINK 工具包需专门购买和安装。 本文介绍 MATLAB5． 1 版的 SIMULINK。 SIMULINK 运行对系统的要求，除要求增加硬盘空间外，其余和 MATLAB 完全一样。 凡是能够运行 MATLAB 核心的配置均能很好地运行 SIMULINK。 4.1.2 SIMULINK 的安装 SIMULINK 的安装与 MATLAB 的安装类似。关于它在 MSWindows 下的安装方法，可 以参考关于 MATLAB 安装的其他有关的书籍。安装完后，会自动在 MATLAB 的目录结构 中加入 SIMULINK 子目录，并自动修改 MATLAB 的启动文件 matlabrc．m。 4.1.3 SIMULINK 入门 为了让同学们较快认识 SIMULINK，以一简单系统为例，从建模到仿真分步叙述如下： （1）在 MATLAB 指令窗运行指令 simulink 便打开图 4-1 所示的 SIMULINK 模块库。图 4-1 SIMULINK 标准模块库41（2）选择［File］菜单中的[New]命令创建一个“Untitled”空白窗。 （3）用鼠标双击信号源（Source）子库方框打开信号源模块库，并调整该信号源窗口 和“Untitled”空白窗互不重叠，如图 4-2 所示。图 4-2 信号发生器模块的复制 （4）将信号发生器（Signal Generator）模块从子库窗口拷贝到“Untitled”窗口。方法 是： 先把鼠标的光标移到信号发生器的方框图中， 按住鼠标的左按钮， 将它拖到新的窗口中， 然后松开按钮，复制过程就完成了。该操作完成后的情况如图 4-2 所示。 （5）模块内部参数的设置。被拷贝到新窗口中的模块包含着与原始模块一样的内部参 数设置。如果想改变信号发生器的参数设置，用鼠标双击它的图标，就会弹出一个如图 4-3 所示的对话框，显示出对信号波形、幅值和频率的控制。 在图 4-3 中，亮的方波是当前选取的波形，频率（Frequency）是 1，幅值（Amplitude） 为 1，其[Frequency][ Amplitude]和[Units]输入框中的数值表示当前的设置值，可以在其中输 入数值重新设置它。另外还能用鼠标通过控制滑标来设置频率和幅值。在完成了参数的设置 后，单击［Apply］按钮保存设置。图 4-3 信号发生器的控制面板42用同样的方法打开信号显示模块库（Sinks） ，从中选取示波器（Scope） ，然后把它拷贝 到新的系统窗口中，如图 4-4 所示。图 4-4 复制示波器 Signal Generator 和 Scope 外面的大于号（＞）分别表示信号的输出和输入。为了连接 这两个模块，使用鼠标的任一个按钮，点击输入或输出端口;看到光标变为（＋）以后，拖 动十字图符到另外一个端口，然后释放鼠标按钮，则带箭头的连线会表示信号的流向，如图 4-5 所示。至此，一个最简单的模型创建完成。图 4-5 绘制连线和模型的建成 （6）仿真操作一：信号发生器的参数设置。图 4-6 示波器的控制面板双击示波器图标后，出现一个虚拟示波器（图 4-6） 。通过鼠标操作，使虚拟示波器与图 4-4 模型窗互不覆盖。选择[Simulation]菜单中的［Start］命令启动仿真过程。在虚拟示波器 上将看到波形。 如果看不到希望的波形，应该先仔细选择示波器的时间轴范围（Horizontal Range）和波 形幅值范围（ Vertical Range） 。 在本例中，因为信号是频率为 1 的单位幅值正弦波。所以，若想在示波器上看到一个周43期以上的波形，就必须把示波器的时间范围设得大于 2，而幅值范围应不小于 1；否则所显 示的波形将被削截。 （7）仿真操作二：算法参数的设置。 在 SIMULINK 中，仿真中动态数据的计算都是由数值积分实现的。尽管本例从信号源 到示波器没通过其他环节，但动态数据仍是经数值积分计算而得的。因此，适当选取算法参 数是必需的。 算法参数的设置方法是：在图 4-5 所示的模型窗口中，在［Simulation］untitled 下拉菜 单中选取[Parameters]选项。于是，出现一个名为“simulink parameters”的对话框。把该框 中的［Max Step Size］的默认值 10 改为 0.1（确保最大步长小于周期的 1/10） 。 经过以上设置后，就可以在示波器上看到如图 4-6 所示的逼真波形。关于算法参数设置 将在下一节作更详细介绍。 4.1.4 界面与菜单 SIMULINK 工 作 窗 是 MATLAB 窗 口 系 列 中 的 一 种 新 类 型 。 下 面 分 别 简 要 说 明 SIMULINK 窗口中各菜单选项的含义， 其中有些菜单的操作方法在后面的有关部分中将陆续 详细介绍。 File 文件操作菜单 New Ctrl＋N 创建新的 SIMULINK 窗口 Open Ctrl＋O 打开已经存在的 SIMULINK 模型文件 Close Ctrl＋W 关闭当前的 SIMUILINK 窗口 Save Ctrl＋S 保存当前的模型文件 Save As 将文件另外保存 Print Ctrl＋P 打印 Print Setup 打印设置 Exit MATLAB Ctrl＋Q 退出 MATLAB Clipboard 剪贴板菜单 Copy 进行拷贝 Copy Options 拷贝选项 在拷贝选项运作时，会弹出对话框提供两种选项： （1）将方框图以 Windows 图元文件 （文件的扩展名是 wmf）的形式剪贴； （2）将方框图以位图（文件的扩展名是 bmp）格式剪 贴。 Edit 编辑菜单 Cut Ctrl＋X 剪切选定的内容 Copy Ctrl＋C 将当前选定的内容拷贝到剪贴板上 Paste Ctrl＋V 将剪贴板上的内容粘贴到当前光标所在位置 Clear 清除选定的内容 Select All 选择整个窗口 Options 建图操作菜单 Group Ctrl＋G 将选定的内容集合成组44Ungroup Ctrl＋U 把集合组还原为分立状态 Mask Ctrl＋M 封装一个新的 SIMULINK 模块 Unmask 打开一个封装好的模块 Flip Horizontal Ctrl＋F 将模块图水平旋转 180 度 Rotate Ctrl＋R 将模块图顺时针旋转 90 度 Reroute Line Ctrl＋L 重布连线，将模块间连线规范化 Simulation 仿真操作菜单 Start／Stop Ctrl＋T 仿真启动或停止 Restart 仿真重启动 Continue／Pause 仿真继续或暂停 Parameters 仿真参数的设置 Style 方框图样式菜单 Drop Shadows 加阴影效果 Orientation 模块输入输出口方向设置 Title 模块名称的设置 Font 字体的设置 Foreground Color 前景颜色 Background Color 背景颜色 Screen Color 屏幕颜色 Sample Time Color 给不同采样时间序列填加颜色 Wide Vector Lines 设置向量的线宽 Update Diagram Ctrl＋D 更新方框图4.2模型的构建为了掌握 SIMULINK，应学习如何操作模块以及如何建立模型。对于 SIMULINK 的模块库 （如图 4-1 所示） 的内容这里不做介绍。 本节着重对构造模型的过程作更为详细的介绍。 SIMULINK 完全采用方框图的“抓取”功能来构造动态系统。系统的创建过程就是绘制 方框图的过程。在 SIMULINK 环境中方框图的绘制完全依赖于鼠标操作。鼠标的不同操作 体现在光标的形状上。在缺省状态下鼠标是一个箭头，但在图形的操作中鼠标呈现不同的形 状。 4.2.1 创建模型文件 建立 SIMULINK 模型通常有三种方式： （1）直接从 MATLAB 指令窗中选取[File]中[New]子菜单的[Model]命令，MATLAB 会 打开一个新方框图窗。当然也可像第 4．1.3 节介绍的那样，在 MATLAB 指令窗下输入 simulink 命令，打开 SIMULINK 模块库窗口，然后再从它的［File]菜单中选取［New]命令， 创建新方框图窗。 （2）如果方框图模型已经存在，那么在 MATLAB 指令窗下直接键入模型文件名字，便 会打开该模型的方框图窗口。用户可以对它进行编辑、修改和仿真。45（3）由于 SIMULINK 模型是以 ASCll 码形式保存的 S 一函数文件，所以还可以使用字 处理软件对它进行创建、修改和编辑，这同建立一个普通 M 文件的过程一样，但要符合创 建 S 一函数的语法规则。 4.2.2 标准模块的选取 由图 4-1 可见，SIMULINK 模块库按模块类型分为七个子库。用鼠标双击子库方框便可 见到库中存放的各种标准模块。这些模块都可被复制到用户的模型窗中（见第 4．1．3 节中 的第四步） 。 打开标准模块的另一个办法是：在 MATLAB 指令窗里键人 blocklib 命令。该指令执行 后，可以见到一组常用标准模块。 4.2.3 模块的移动、删除和拷贝 模块的移动 该操作非常简单：将光标置于待移动模块图标上，然后按住鼠标将模块拖到合适的地方 即可。模块移动时，它与其他模块的连线也随之移动。 模块的选定 模块选定操作是许多其他操作（如删除、剪切、拷贝）的“前导性”操作。选定模块的 方法有两种： (1）用鼠标单击待选模块，模块四个角处便出现小黑块，表示已经选定。 （2）如果选择一组模块，可以按住鼠标按钮拉出一个矩形虚线框，将所有待选模块包在 其中，然后松开按钮，则矩形里所有的模块同时被选中。 模块的删除、剪切和拷贝 对选定模块的删除、剪切和拷贝操作分述如下： (1)按[Delete]键，把选定模块删除。 （2）选择［Edit］菜单中［Cut］命令后，便将选定模块移到 Windows 的剪贴板上，可 供 Pastel 命令重新粘贴。 （3）运行[Edit]菜单中[Copy]命令；然后将光标移到将粘贴的地方，按一下鼠标；看到选 定的模块恢复原状，再运行［Edit］菜单中的 [Paste］命令，就会在选定的位置上复制出相 应的模块。新复制的模块和原模块的名称也会自动编号，以资区别。 另一种更简单的复制操作是：先按下［Ctrl］键不放，然后用鼠标把待拷贝模块拖到希 望位置后，松开鼠标左键，便完成拷贝工作。 4.2.4 模块的连接 关于模块连接，在本章一开始就已经涉及。值得提醒的是：连接模块时，要注意模块的 输入输出口和各模块间的信号流向。在 SIMULINK 中，模块总是由输入口接受信号，从输 出口发送信号。 [例 1]建立如图 4-7 所示的模型46图 4-7 一个简单系统的构造 （1）利用 SIMULINK 窗中［File］菜单的[New]选项开辟一个新的“Untitled”窗口。 （2）分别从信号源库（Sourse） 、线性模块库(Linear)、数据接受器库（Sinks）中，用鼠 标把信号发生器（Signal Generator） 、求和模块（Sum）和传递函数模块（Transfer Fcn） 、示 波器（Scope）拖取到“Untitled”窗口。各模块的位置如图 4-7 所示。 （3）图 4-7 中模块间的连线有两类：一类是“直”线，它从某模块的输出口出发直指另 一模块的输入口。这种连线的生成方法在第 4．1．3 节中已经讲过。另一类连线是“折”线。 这类“折”线的生成方法是：把鼠标光标移到传递函数模块和示波器间连线的中点附近，按 下鼠标右键，光标由箭头变为十字，往下拉动鼠标到适当位置后放开右键，屏幕上就出现一 条由中点引出的箭头线，再从此箭头开始用鼠标水平向左画线到适当位置，再松开鼠标键， 照此操作，直到整个“折”线绘完。 4.2.5 模块属性的改变 模块的属性可以分为两种：模块的标题和模块的内部参数。为适合自己的需要，那些被 用户所复制的标准模块的标题和参数常需作必要的修改。下面就介绍修改方法。 标题的修改 模块标题是指标识模块图标的字符串。 对用户所建模型窗中模块标题进行修改的具体方 法如下： （1）用鼠标单击标题，使之增亮反显。 （2）输入新的标题（中西文均可） 。 （3）然后用鼠标点去窗口中的任一地方，修改工作结束。 模块内部参数的修改 双击待修改参数模块的图标，打开参数设置对话框，然后通过改变对话框适当栏目中的 数据便可。下面通过具体算例来说明具体操作。 ［例 1］把图 4-7 方框图模型修改成图 4-8 所示的模型。注意到本系统是以文件 sim1.m 保存的。47图 4-8系统 sim1．m（1）对标题的修改 用鼠标点亮传递函数块标题“Transfer Fen” ，输人汉字字符“传递函数” 。再点亮示波器 模块标题“Scope” ，输人汉字字符“示波器” ，再在窗口空白处按点鼠标左键，标题修改结 束。 （2）对求和块输入极性的修改 用鼠标双击求和块（Sum） ，就会弹出对话框。 把［List of signs］栏中的缺省极性改成如图 4-9 所示形式，再点[Apply]键，原求和模 块图标使自动改成图 4-8 所示形式。图 4-9 求和模块对话框 （3）对传递函数块参数的修改 用鼠标双击传递函数模块图标，弹出如图 4-10 所示对话框。图 4-10 传递函数模块对话框48［Denominator］栏表达的是传递函数的分母多项式系数向量。把该栏中的原缺省值改 成图 4-10 所示的行向量，再点[Apply]键，原传递函数模块图标中的函数表达式就自动变成 图 4-8 中的形式。 值得指出的两点是： （1）在参数设置时，任何 MATLAB 工作内存中已有的变量、合法 表达式、MATLAB 语句等都可以填写在设置栏中； （2）模块图标的大小是可以用鼠标操作 调整的。 因此， 假如传递函数表达式太长， 原方框容纳不下， 可以用鼠标把它拉到适当大小， 使整个方框图标美观易读。 4.2.6 模型文件的保存 构造好一个模型后，选择［File］菜单中的［Save］命令将模型存盘。文件是以 ASCll 码形式存储的 M 文件。文件包含了该模型的所有信息（数学模型内涵和外部框图表现）. 例如，当上小节例 2 的全部步骤执行完后，再在“Untitled”模型视窗中的［File］菜单 中选[Save]子项，并在弹出对话框的［文件名］栏填写 sim1．m，再点［确定］便完成了文 件的保存。以后在 MATLAB 指令窗下运行 sim1 时，MATLAB 会重新打开此系统窗口。 由 SIMULINK 创建的 M 文件也可以用 type 指令查看其内容， 当然也可以用字处理器进 行修改和编辑。此文件是 MATLAB 自动创建的标准 S 一函数，前面以字符“％”开头的物 理行是函数的注释行，紧接着程序定义模型的方框图形信息，程序最后部分是实现数学模型 仿真运算的指令。49第五章自动控制系统设计任务 任务一1．*某单位负反馈系统，前向通道的传递函数为 G(s) ?9 ，试设计校正装置使 S(S ? 2)系统的相角裕量 ? ≥45° ，截止频率 ? C 不小于 4rad/s。 单位负反馈系统的开环传递函数为 角裕量 ? ≥65° ，幅值裕量为 34db。 3.** 对于结构如图 5-1 所示系统，给定固有部分的传递函数 Gg（S）和性能指标要求，试 分别设计串联校正装置 K（s） ，并比较它们的作用效果。 设 Gg(s)2.*G(s) ?20 ，设计系统使相 s(0.001s ? 1)(0.5s ? 1)?100 s 2 (0.1s ? 1)(0.01s ? 1)1）若要求开环比例系数 K≥100 S-1，相角 裕量 ? ≥30° ? C ≥45 S-1； ， 2）若要求开环比例系数 K≥100 S-1，相角裕量 ? ≥40° ? C =5 S-1； ， 3）若要求开环比例系数 K≥100 S-1，相角裕量 ? ≥40° ? C ≥20 S-1； ，4．* 被控对象传递函数为： G ( s ) ?10 ( s ? 1)( s ? 2)( s ? 3)设计反馈控制器 u=-kx，使闭环系统的极点为 S1, 2 ? ?2 ? j 2 3 ，s3=-10。 5．** 含积分环节的伺服系统设计。设对象为 G ( s) ?1 ， s( s ? 2)( s ? 1)设计反馈控制器 u=-kx+k1r，使闭环系统的极点为 ? 2 ? j 2 3,?10 。506．***不含积分环节的伺服系统设计。对车载倒立摆系统实现控制，使小车位置作为输出的 闭环系统具有极点： ? 1 ? j 3,?5,?5,?5 。51任务二1． *画出开环传递函数为 G(s) ? 时的 BODE 图。K(s ? 2) 的控制系统的根轨迹和 K=16 s(s ? 4)(s ? 8)(s 2 ? 2s ? 5)2． *已知单位负反馈系统的开环传递函数为 G(s) ? BODE 图，并求出系统的相对稳定裕量。200 ，试画出系统的 S(S ? S ? 100 )23． **试用 MATLAB 命令连接图示系统， （编写 m 文件） ，求出系统的开环和闭环传函，并 画出系统的 BODE 图。52任务三1. *已知单位反馈系统的前向通道的传递函数为 G(s) ?35 ，试建立该系统的仿真连接 S ?S2图，通过编写 M 文件作出单位阶跃输入下的系统响应。 2．**建立如下控制系统模型，并研究单位阶跃输入下的系统响应。画出仿真连接图，通过 编写 M 文件作出响应曲线。3．***已知系统如图下所示，对象的传函为 G P (s) ?1 ，设采样周期为 T=1s， S(S ? 1)G c (z) ?5.4(z - 0.5)(z - 0.37) ，试研究输入为单位速度信号时的系统响应。 (z - 1)(z ? 0.72)53任务四**1. 已知单位负反馈系统被控对象的传递函数为: G0 ( s ) ? K 01 s(0.1s ? 1)(0.001s ? 1)试用 Bode 图设计法对系统进行超前串联校正设计，使之满足： （1）在斜坡信号 r(t)=R t 作用下，系统的稳态误差 ess≤0.001R; （2）系统校正后,相角稳定裕度 43?&γ&48? 。 **2．已知单位负反馈系统被控对象传递函数为： G0 ( s) ? K 01 s(0.1s ? 1)(0.2s ? 1)试用 Bode 图设计法对系统进行滞后串联校正设计，使之满足： （1）在单位斜坡信号 r(t)=t 作用下，系统的速度误差系数 Kv≥30s-1; （2）系统校正后剪切频率 ωc≥2.3s-1; （3）系统校正后相角稳定裕度 γ&40? 。 ****3.已知某系统的开环传递函数为: G0 ( s ) ? K 01 , s( s ? 1)( s ? 2)试用 Bode 图设计法对系统进行超前-滞后串联校正设计，使之满足： （1）在单位斜坡信号 r(t)=t 作用下，系统的速度误差系数 Kv=10s-1; （2）系统校正后剪切频率 ωc≥1.5s-1; （3）系统校正后相角稳定裕度 γ≥45? 。 （4）计算校正后系统的时域指标:σ%≤25%, ***4. 已知某系统的开环传递函数为: G0 ( s ) ? K 0 tp≤2s, ts≤6s。1 , s ( s ? 2)试用根轨迹几何设计法对系统进行超前串联校正设计，使之满足： （1）阶跃响应的超调量：σ%≤20%； （2）阶跃响应的调节时间：ts≤1.2s。 ***5. 已知单位负反馈系统被控对象传递函数为： G0 ( s ) ? K 02500 s( s ? 25)试用根轨迹几何设计法对系统进行滞后串联校正设计，使之满足： （1）阶跃响应的超调量：σ%≤15%； （2）阶跃响应的调节时间：ts≤0.3s； （3）单位斜坡响应稳态误差 ess≤0.01。54**6. 已 知 过 程 控 制 系 统 的 被 控 广 义 对 象 为 一 带 延 迟 的 惯 性 环 节 ， 其 传 递 函 数 为 ：G0 ( s) ?8 e ?180s 360 s ? 1试用 Ziegler-Nichols 整定公式计算系统 P、PI、PID 校正器的参数，并进行单位阶跃响 应的仿真。 **7. 已知一系统的传递函数为：Y ( s) 10 ? U ( s) s( s ? 1)( s ? 2)试判别系统的可控性并设计反馈控制器，使得闭环系统的极点为-2，-1±i。? **8. 已知一系统的状态方程为： x ? Ax ? Bu 其中 A=[0 1 0 0；0 0 C1 0；0 0 0 1；0 0 11 0],B=[0；1；0；-1]。 试判别系统的可控性并设计状态反馈控制器，使得闭环系统的极点为-1、-2、-1±i。? x ? Ax ? Bu y ? cx**9. 已知一系统的状态方程为：，方程中 A=[0 1；-2 -3]; B=[0；1]; c=[2 0]。 试判别系统的可观性并设计状态观测器，使得闭环系统的极点为-10，-10。***10. 已知一系统的状态方程为：? x ? Ax ? Bu y ? cx， B=[0；1；0；-1], c=[1 0 0 0].其中 A=[0 1 0 0；0 0 C1 0；0 0 0 1；0 0 11 0], （1）试判别系统的可观性；（2）若系统可观，试设计全状态观测器，使得闭环系统的极点为-2，-3，-2±i。 （3）设计三阶状态观测器，使得闭环系统的极点为 -3，-2±i。55附 录附录一 时域分析例1?2 n 已知二阶系统的传递函数为： ? (s) ? 2 s ? 2?? n s ? ? 2 n当 ωn=1 时，试计算 ξ 从 0.1 变至 1 时二阶系统的响应，并绘制一簇阶跃响应曲线。程序及结果如下： %zeta 与时域响应关系 % num=1;y=zeros(200,1);i=0; for bc=0.1:0.1:1 den=[1,2*bc,1]; t=[0:0.1:19.9]'; sys=tf(num,den); i=i+1; y(:,i)=step(sys,t); end mesh(flipud(y),[-100 20])561.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 1050 200180160140120100806040200例2已知二阶系统为： G (s) ?k ， s ? cs ? k2c=｛1,2,4｝,k=｛1.25,2,29｝ 试绘制该系统所对应的三组不同参数下的阶跃响应曲线（在同一坐标下） 。 程序及结果如下： c=[1 2 4]; k=[1.25 2 29]; t=linspace(0,10,100)'; for j=1:3 num=k(j); den=[1 c(j) k(j)]; sys=tf(num,den); y(:,j)=step(sys,t); end plot(t,y(:,1:3)),grid gtext('a=1 b=1.25'), gtext('a=2 b=2'), gtext('a=4 b=29')57例 22 已知单位反馈系统开环传函为： G (s) ?k s（0.5s ? 1)(4s ? 1)试绘制 k=1.4,2.3,3.5 时单位阶跃响应曲线（在同一坐标下） 。 %作出不同开环增益时的阶跃响应曲线，并求某一 K 值时的性能指标 num=1;den=conv(conv([1 0],[0.5 1]),[4,1]); rangek=[1.4 2.3 3.5]; t=linspace(0,20,200)'; for i=1:3 s1=tf(num*rangek(i),den); sys=feedback(s1,1); y(:,i)=step(sys,t); end plot(t,y(:,1:3)),grid gtext('k=1.4'), gtext('k=2.3'), gtext('k=3.5')583.53k=3.52.5 k=2.3 2 k=1.41.510.50-0.5-1-1.502468101214161820例3设控制系统开环传递函数为： G(s) ?1.25 ， s2 ? s试绘制该闭环系统的单位阶跃响应曲线，并计算系统的性能指标。 程序及结果如下： % %计算时域性能指标(tp,ess,b1,b2,sigma,n,阻尼比，T,f) global y t sys=tf(1.25,[1 1 0]); Gc=feedback(sys,1); step(Gc); [y,t]=step(Gc); [mp,tf]=max(y); tp=t(tf); ct=length(t); tm=max(t); yss=y(ct); q=1; m=q-1; while m&3, for a=(tm/100):.01:tm59j=[0:a:tm]; for i=1:length(j); if (y(i+1)-y(i))&0&(y(i)-y(i-1))&0, m=m+1; pm(m)=y(i); tp(m)=t(i); end end end end yss=y(ct);ess=1- b1=pm(1)- b2=pm(2)- sigma=100*b1/ n=b1/b2; pusi=(b1-b2)/b1; T=(tp(2)-tp(1)); f=1/T; 结 果 ： yss=0.99871 ， ess=0.0012948, n=20.336， pusi=0.95083, b1=0.20885, T=6.2945, b2=0.01027, f=0.15887 sigma=20.91260Step Response 1.41.210.8Amplitude0.60.40.200246 Time (sec)81012例4已知系统开环传递函数为： (s) ? 100 G(s ? 2) ， 试对系统的稳定性进行判断。 s(s ? 1)(s ? 20 )程序及结果如下： k=100; z=[-2]; p=[0,-1,-20]; [n1,d1]=zp2tf(z,p,k); P=n1+d1; roots(P) 根为 ans = -12.0 -3.1010 特征根实部均为负值，所以闭环系统稳定。61例5已知单位反馈开环零极点增益模型为： G (s) ?3(2s ? 1) ， s(s ? 2)(s ? 1)试绘制该系统单位斜坡响应曲线，并计算单位斜坡响应的稳态误差。 解：1）判稳： （用 roots()函数） %P214L4-16 J1 k=6;z=-0.5;p=[-2 1 0]; [n1,d1]=ZP2TF(z,p,k); s=tf(n1,d1); sys=feedback(s,1); roots(sys.den{1}) ans =-0.1084 + 1.9541i -0.1084 - 1.9541i -0.7832 所有根均有负实部，故闭环系统稳定。 2）求单位阶跃响应曲线及稳态误差： % P214L4-16 J2 k=6; z=-0.5; p=[-2 1 0]; [n1,d1]=zp2tf(z,p,k); s=tf(n1,d1); sys=feedback(s,1); t=[0:0.1:30]'; y=step(sys,t); subplot(121),plot(t,y),grid es=1-y; subplot(122),plot(t,es),grid622.51.5211.50.5100.5-0.50-1-0.50102030-1.50102030?× ?? ??? ? ?°?× ?? ? ???? ? ì ?ú ? ó ?ì ?ú ?单位阶跃响应曲线及其误差响应曲线 3）求单位斜坡响应曲线及稳态误差： % % 求单位斜坡响应曲线及稳态误差 k=6; z=-0.5; p=[-2 1 0]; [n1,d1]=zp2tf(z,p,k); s=tf(n1,d1); sys=feedback(s,1); t=[0:0.1:50]'; num=sys.num{1}; den=[sys.den{1},0]; sys=tf(num,den); y=step(sys,t); subplot(121),plot(t,[t y]),grid es=t-y; subplot(122),plot(t,es),grid ess=es(length(es))63ess = -0.667860 0.60.4 50 0.2400-0.2 30 -0.420-0.6-0.8 10 -1001020304050-1.201020304050? ? ?± ? ? ??? ? ?°? ? ? ? ??? ? ?? ì ?ú ? ó ?ì ?ú ?单位斜坡响应曲线及其误差响应曲线64附录二根轨迹分析(s ? 3) s(s ? 5)(s ? 6)(s 2 ? 2s ? 2)例1 已知一单位反馈系统开环传函为 G (s) ? k试在根轨迹上选择一点，求出该点的增益 k 及其闭环极点的位置，并判断在该点系统的稳定 性。 程序： num=[1,3]; den=conv(conv(conv([1 0],[1 5]),[1 6]),[1 2 2]); rlocus(num,den); [k,poles]=rlocfind(sys); range=[33:1:37]'; cpole=rlocus(num,den,range); [range,cpole] 结果： selected_point = -5.3780 - 0.0476i ans = Columns 1 through 5 33.5 + 0.6697i -5.5745 - 0.6697i -1.0 + 1.0 -5.5768 + 0.6850i -5.5768 - 0.6850i -1.5 + 1.0 -5.5791 + 0.7001i -5.5791 - 0.7001i -1.2 + 1.0 -5.5815 + 0.7147i -5.5815 - 0.7147i -1.8 + 1.0 -5.5838 + 0.7291i -5.5838 - 0.7291i -1.6 + 1.3712i Column 6 -0.0260 - 1.3210i -0.0155 - 1.3340i -0.0052 - 1.8 - 1.6 - 1.3712i &&例2已知带有延迟因子的系统开环传递函数为： G (s) ? 1） 试绘制根轨迹图； 2） 求系统临界稳定时根轨迹增益；651 e ?s s(s ? 1)( 0.5s ? 1)3） 求系统 k=0.5 时单位阶跃响应曲线。 程序： n1=[1]; d1=conv(conv([1 0],[1 1]),[0.5 1]); %s1=tf(n1,d1); [np,dp]=pade(1,3); G=tf(n1,d1)*tf(np,dp);figure(1); rlocus(G);hold on [k,p]=rlocfind(G); %sys=feedback(G,1);figure(2) %step(sys); selected_point =-2.7605 - 0.1612iRoot Locus 15105Imag Axis0-5-10 -12 -10 -8 -6 -4 -2 Real Axis 0 2 4 6 866附录三频域分析例1 已知单位负反馈系统前向通道的传递函数为：G (s) ?2s 4 ? 8s 3 ? 12s 2 ? 8s ? 2 s 6 ? 5s 5 ? 10s 4 ? 10s 3 ? 5s 2 ? s试绘制出 Bode 图并计算系统的频域性能指标。 程序及结果如下： num=[0 0 2 8 12 8 2];den=[1 5 10 10 5 1 0]; sys=tf(num,den); [mag,phase,w]=bode(sys); [gm,pm,wcp,wcg]=margin(mag,phase,w); margin(mag,phase,w) kg=20*log(gm) 结果：gm=332.17, pm=38.692, wg=25.847, kg=50.4272wc=1.249Bode Diagram Gm = 50.427 dB (at 25.847 rad/sec), Pm = 38.692 deg (at 1.249 rad/sec) 40 30 20Magnitude (dB) Phase (deg)10 0 -10 -20 -30 -40 -90-135-180 -1 10100101Frequency (rad/sec)67例2 已知一带延迟因子的系统开环传函为： G (s) ? 10s ? 5 -0.8s e ， (s ? 3) 2试求其有理传递函数的频率响应，同时在同一坐标中绘制以 Pade 近似延迟因子式系 统的 Bode 图，并求此时系统的频域性能指标。 程序及结果如下： n1=[10 50]; d1=conv([1 3],[1 3]); s1=tf(n1,d1); w=logspace(-1,2); [mag1,phase1]=bode(n1,d1,w); [np,dp]=pade(0.8,4); %s1=tf(n1,d1)*tf(np,dp); numt=conv(n1,np); dent=conv(d1,dp); s2=tf(numt,dent); [mag2,phase2]=bode(numt,dent,w); subplot(211); semilogx(w,20*log10(mag1),w,20*log10(mag2),'r--');grid on title('bode plot'); xlabel('Frequency（rad／sec）'); ylabel('Gain db'); subplot(212);grid semilogx(w,phase1,w,phase2,'r--');grid on xlabel('Frequency（rad／sec）'); ylabel('Phase deg'); [gm1,pm1,wcp1,wcg1]=margin(s1); [gm2,pm2,wcp2,wcg2]=margin(s2);68例3 已知两个单位负反馈系统开环传递函数分别为：G 1 (s) ?2.7 s ? 5s 2 ? 4s3， G 2 (s) ?2.7 s ? 5s 2 ? 4s3试用 Bode 图法判断闭环系统的稳定性。 程序如下： num1=[0 0 0 2.7]; den1=[1 5 4 0]; sys1=tf(num1,den1); figure(1);hold on [gm,pm,wcp,wcg]=margin(sys1); margin(sys1); title('对数频率特性图 1'); xlabel('频率 rad/sec'); ylabel('Gain dB');num2=[0 0 0 2.7]; den2=[1 5 -4 0]; sys2=tf(num2,den2); figure(2); [gm,pm,wcp,wcg]=margin(sys2); margin(sys2); title('对数频率特性图 2'); xlabel('频率 rad/sec'); ylabel('Gain dB'); 结论：系统 1 稳定，系统 2 不稳定（Warning: The closed-loop system is unstable.）69? ? ? ? ?? ? ?????? ? Gm = 17.393 dB (at 2 rad/sec), Pm = 51.732 deg (at 0.57831 rad/sec) 500-50-100Gain dB-150 -90-135-180-225-270 -1 1010010 ? ? rad/sec ??1102? ? ? ? ?? ? ?????? ? Gm = Inf, Pm = -58.05 deg (at 0.53456 rad/sec) 500-50-100Gain dB-150 18013590 -1 1010010 ? ? rad/sec ??110270
更多相关文档