1、MATLABMATLAB控制系统仿真控制系统仿真提纲v一、弹簧一、弹簧-重物重物-阻尼器系统阻尼器系统v二、传递函数二、传递函数v三、结构图模型三、结构图模型引言 MATLAB是一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算和图形显示于一体，构成了一个方便的界面友好的用户环境。 控制系统的分析与设计方法，都是以数学模型为基础进行的。MATLAB可以用于以传递函数形式描述的控制系统。 在本节中，首先举例说明如何使用MATLAB进行辅助分析。然后讨论传递函数和结构图。 一、弹簧一、弹簧-重物重物-阻尼器系统阻尼器系统 弹簧重物阻尼器动力学系统如图2-1所示。重物M的位 移由y(t)表示，
2、用微分方程描述如下：22d( )d ( )( )( )ddy ty tMBKy tf ttt该系统在初始位移作用下的瞬态响应为：22(0)( )sin(1)1ntnyy tet 其中 cos-1 ，初始位移是y(0)。 系统的瞬态响应当1时为欠阻尼，当时为过阻尼，当1时为临界阻尼。过阻尼情况：y(0)=0.15 m n= (弧度/秒) （ ）欠阻尼情况：y(0)=0.15 m n= (弧度/秒) （ ）2132 22 ,3KBMM2212 22 ,1KBMM 利用MATLAB程序unforced.m，可以显示初始位移为y(0)的物体自由运动曲线，如图2-63所示。 在unforced.m程序中
3、，变量y(0)，n，t， 1和 2的值由指令直接输入工作区，然后运行unforced.m程序就可以产生响应曲线。 y0=0.15;wn=sqrt(2);zeta1=3/(2*sqrt(2); zeta2=1/(2*sqrt(2);t=0:0.1:10;unforced （a）MATLAB指令窗口 * 计算系统在给定初始条件下的自由运动t1=acos(zeta1)*ones(1,length(t);t2=acos(zeta2)*ones(1,length(t);c1=(y0/sqrt(1-zeta12); c2=(y0/sqrt(1-zeta22);y1=c1*exp(-zeta1*wn*t)s
4、in(wn*sqrt(1-zeta12)*t+t1);y2=c2*exp(-zeta2*wn*t)sin(wn*sqrt(1-zeta22)*t+t2); * 计算运动曲线的包络线bu=c2*exp(-zeta2*wn*t);bl=-bu; * 画图plot(t,y1, -,t,y2,-,t,bu, -,bl, -),gridxlabel(Timesec), ylabel(y(t) Displacementm)text(0.2,0.85,oeverdamped zeta1=,num2str(zeta1), )text(0.2,0.80,underdamped zeta2=,num2str(ze
5、ta2), ) （b）分析弹簧重物阻尼器的MATLAB程序unforced.m图2-63 分析弹簧重物阻尼器的MATLAB指令图2-64 弹簧重物阻尼器的自由运动曲线 在欠阻尼和过阻尼情况下的响应曲线如图2-64所示 ： MATLAB可分析以传递函数形式描述的系统。分子多项式和分母多项式都必须在MATLAB指令中指定。 在MATLAB中多项式由行向量组成，这些行向量包含了降次排列的多项式系数。例如多项式p(s)=1s3+3s2+0s1+4s0，按图2-65的格式输入p=1 3 0 4， p=1 3 0 4;r=roots(p)r= -3.3553e+00 1.7765e-01+1.0773e+
6、00j 1.7765e-01-1.0773e+00jp=poly(r)p= 1.000 3.000 0.000-0.000j 4.000+0.000j图2-65 输入多项式并求根 矩阵乘法由MATLAB的conv()函数完成。把两个多项式相乘合并成一个多项式n(s)，即： n(s)= (3s2 +2s +1) (s +4) = 3s3 +14s2 +9s +4 与此运算相关的MATLAB函数就是conv()。函数polyval()用来计算多项式的值。多项式n(s)在s = -5处值为n (-5) = -66，见图2-66。p=3 2 1;q=1 4;n=conv(p,q)n= 3 14 9 4
7、value=polyval(n,-5)value= -66图2-66 MATLAB的conv()函数和polyval()函数 设传递函数为G(s)=num/den，其中num和den均为多项式。利用函数：二、传递函数二、传递函数P , Z=pzmap(num , den) 可得G(s)的零极点位置，即P为极点位置列向量，Z为零点位置列向量。该指令执行后自动生成零极点分布图。 考虑传递函数：23261( )331sG ssss 和(1)(2)( )(2 )(2 )(3)ssH ssi si s图2-67 零极点图 传递函数G(s)/H(s)的零极点图如图2-67所示，相应的MATLAB指令如图2
8、-68所示。 numg=6 0 1;deng=1 3 3 1;z=roots(numg)z= 0+0.4082j0-0.4082j p=roots1(deng)p= -1 -1 -1n1=1 1; n2=1 2; d1=1 2*j; d2=1 2*j; d3=1 3;numh=conv(n1,n2); denh=conv(d1,conv(d2,d3);num=conv(numg,denh); den=conv(deng,numh);printsys(num,den)num/den= 6s5+18s4+25s3+ 图2-68 绘制零极点图指令三、三、结构图模型结构图模型 一个开环控制系统可以通过
9、G1 (s)与G2 (s)两个环节的串联而得到，利用series()函数可以求串联连接的传递函数，函数的具体形式为： num,den= series(num1,den1, num2,den2) 例如G1 (s)和G2 (s)的传递函数分别为：11( )2sG ss221( )500G ss12232( )111( )( )( )2 5005001000C sssG s G sR sssss则串联函数的用法示于图2-69： num1=1;den1=500 0 0;num2=1 1;den2=1 2;num,den=series(num1,den1,num2,den2);printsys(num,
10、den)num/den= s+1 500s3+1000s2图2-69 series函数的用法 当系统是以并联的形式连接时，利用parallel()函数可得到系统的传递函数。指令的具体形式为： 系统以反馈方式构成闭环，则系统的闭环传递函数为：( )( )1( )( )BG sGsG s H snum,den= parallel (num1,den1, num2,den2) 求闭环传递函数的MATLAB函数有两个：cloop()和feedback()其中cloop()函数只能用于H (s)=1（即单位反馈）的情况。 cloop()函数的具体用法为：num,den= cloop (numg,deng
11、, sign)其中numg和deng分别为G (s)的分子和分母多项式，sign=1为正反馈，sign= -1为负反馈（默认值）。 feedback()函数的用法为：num,den= feedback (numg,deng,numh,denh, sign) 其中numh为H (s)的分子多项式，denh为分母多项式。 闭环反馈系统的结构图如图2-70所示，被控对象G(s)和控制部分Gc (s)以及测量环节H (s)的传递函数分别为：，1( )2csG ssnumcdenc21( )5G ssnumgdeng1( )10H ssnumhdenh图2-70 闭环反馈系统的结构图 应用series(
12、)函数和feedback()函数求闭环传递函数的MATLAB指令如图2-71 所示： numg=1;deng=5 0 0;numc=1 1;denc=1 2;numh=1;denh=1 10;num1,den1=series(numc,denc,numg,deng);num,den=feedback(num1,den1,numh,denh,-1);printsys(num,den)num/den= s2+11s+10 5s4+60s3+100s2+s+1图2-71 feedback()函数的应用 例例2.12 一个多环的反馈系统如图2-49所示，给定各环节的传递函数为：11( )10G ss2
13、1( )1G ss2321( )44sG sss41( )6sG ss1( )1H s 2( )2Hs 31( )2sHss试求闭环传递函数GB(s)=C(s)/R(s)。解解 求解步骤如下：步骤1：输入系统各环节的传递函数；步骤2：将H2的综合点移至G2后；步骤3：消去G3，G2，H2环；步骤4：消去包含H3的环；步骤5：消去其余的环，计算GB (s)。根据上述步骤的MATLAB指令以及计算结果在图2-72中。ng1=1;dg1=1 10;ng2=1;dg2=1 1;ng3=1 0 1;dg3=1 4 4;ng4=1 1;dg4=1 6;nh1=1;dh1=1;nh2=2;dh2=1;nh3
14、=1 1;dh3=1 2;n1,d1=series(ng2,dg2,nh2,dh2);n2,d2=feedback(ng3,dg3,n1,d1,-1);n3,d3=series(n2,d2,ng4,dg4);n4,d4=feedback(n3,d3,nh3,dh3,-1);n5,d5=series(ng1,dg1,ng2,dg2);n6,d6=series(n5,d5,n4,d4);n7,d7=cloop(n6,d6,-1);printsys(n7,d7)num/den= s4+ 3s3+ 3s2+3s+2 2s6+38s5+261s4+1001s3+1730s2+1546s+732图2-72
15、 多环结构图简化 通过pzmap()或roots()函数可查看传递函数是否有相同的零极点，还可使用minreal()函数除去传递函数共同的零极点因子。 如图2-73所示。numg=1 6 11 6;deng=1 7 12 11 5;printsys(numg,deng)numg/deng= s3+6s2+11s+6 s4+7s3+12s2+11s+5num,den=minreal(numg,deng);printsys(num,den)1 pole-zeros cancellednum/den= s2+4s+3 s3+6s2+6s+5图2-73 minreal()函数的应用 例2.2所示的位置
16、随动系统，在给定各元件参数并忽略La和令ML = 0的情况下，其结构图如图2-74所示： 图2-74 位置随动系统的结构图 第一步求闭环传递函数GB (s)= c(s) / r(s)，求解过程及结果如图2-75所示。第二步利用step()函数计算参考输入 r (t)为单位阶跃信号时输出 c (t)的响应。 num1=200; den1=20;num2=1; den2=2 0.5 0;num3=0.2 0; den3=1; num4=540; den4=1;na,da=series(num1,den1,num2,den2);nb,db=feedback(na,da,num3,den3,-1);n
17、c,dc=series(nb,db,num4,den4);num,den=cloop(nc,dc,-1);printsys(num,den)num/den= 5400 2s2+2.5s+5400 t=0:0.005:3;y,t=step(num,den,t);plot(t,y),grid图2-75 位置随动系统的结构图简化及阶跃响应指令图2-76 位置随动系统的阶跃响应曲线 图2-76给出了位置随动系统的阶跃响应曲线。用plot()函数用于画出y(t)曲线，grid函数用于给图形加上网格。2.7 循序渐进设计示例：磁盘驱动读取系统循序渐进设计示例：磁盘驱动读取系统我们指出了磁盘驱动系统的基本设
18、计目标：尽可能将磁头准确定位在指定的磁道上，并且磁头从1个磁道转移到另1个磁道所花的时间不超过10ms。在这里，我们将完成设计流程的第4、5步。首先应选定执行机构、传感器和控制器，然后建立控制对象和传感器等元部件的模型。根据表2.1磁盘驱动读取系统的典型参数，我们有)1000)(20(5000)()(sssRLsbJssKsGmG(s)还可以改写成) 1)(1(/)(sssbRKsGLm图2-8 磁盘驱动器读取系统框图模型表2.1磁盘驱动器读取系统典型参数图2-78 磁盘驱动器读取系统框图模型其中 。由于 ，因此常被略去，有：sRLsbJL1/,50/LL或)20(5)(sssG该闭环系统的框
19、图模型见图2-9。利用框图变换化简规则，有：)(1)()()(sGKsGKsRsYaa利用G(s)的2阶近似表示，可以有：aaKssKsRsY5205)()(2) 105. 0(25. 0) 1(/)(ssssbRKsGLm当取Ka=40时，最后可得：)(20020200)(2sRsssY2-79 闭环系统的框图模型 使用MATLAB的函数step，可以得到 时，如图2-80所示的系统 阶跃响应。radssR1 . 0)(小结 本章讨论了如何建立控制系统以及元部件的数学模型问题。本章所涉及的数学模型共有三种，即微分方程、传递函数、结构图或信号流图。 利用传递函数研究线性系统，可根据传递函数的极点和零点分布判定系统对不同输入信号的响应特性。 利用结构图或信号流图可以了解系统中的每个变量，还可以通过梅逊（Mason）公式，方便地求得系统输入输出间的传递函数。 利用MATLAB软件可求解系统在不同参数和输入情况下的响应。
