第五节 应用MATLAB控制系统仿真

提纲 一、弹簧-重物-阻尼器系统 二、传递函数 三、结构图模型

引言 MATLAB是一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算和图形显示于一体，构成了一个方便的界面友好的用户环境。

一、弹簧-重物-阻尼器系统 弹簧—重物—阻尼器动力学系统如图2-1所示。重物M的位 移由y(t)表示，用微分方程描述如下： 该系统在初始位移作用下的瞬态响应为： 其中q =cos-1z ，初始位移是y(0)。 系统的瞬态响应当z＜1时为欠阻尼，当z＞１时为过阻尼， 当z＝1时为临界阻尼。

过阻尼情况：y(0)=0.15 m wn= (弧度/秒) （ ） 欠阻尼情况：y(0)=0.15 m wn= (弧度/秒) （ ） 利用MATLAB程序—unforced.m，可以显示初始位移为y(0) 的物体自由运动曲线，如图2-63所示。 在unforced.m程序中，变量y(0)，wn，t，z 1和z 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)); * 计算系统在给定初始条件下的自由运动 t1=acos(zeta1)*ones(1,length(t)); t2=acos(zeta2)*ones(1,length(t)); c1=(y0/sqrt(1-zeta1^2)); c2=(y0/sqrt(1-zeta2^2)); y1=c1*exp(-zeta1*wn*t)sin(wn*sqrt(1-zeta1^2)*t+t1); y2=c2*exp(-zeta2*wn*t)sin(wn*sqrt(1-zeta2^2)*t+t2); * 计算运动曲线的包络线 bu=c2*exp(-zeta2*wn*t);bl=-bu; * 画图 plot(t,y1, ‘-’,t,y2,‘-’,t,bu, ‘--’,bl, ‘--’),grid xlabel(‘Time[sec]’), ylabel(‘y(t) Displacement[m]’) text(0.2,0.85,[‘oeverdamped zeta1=’,num2str(zeta1),] ) text(0.2,0.80,[‘underdamped zeta2=’,num2str(zeta2),] ) （b）分析弹簧—重物—阻尼器的MATLAB程序unforced.m 图2-63 分析弹簧—重物—阻尼器的MATLAB指令

在欠阻尼和过阻尼情况下的响应曲线如图2-64所示 ： 图2-64 弹簧—重物—阻尼器的自由运动曲线   MATLAB可分析以传递函数形式描述的系统。分子多项式和 分母多项式都必须在MATLAB指令中指定。

在MATLAB中多项式由行向量组成，这些行向量包含了降次 排列的多项式系数。例如多项式p(s)=1s3+3s2+0s1+4s0，按图2-65 >>r=roots(p) r= -3.3553e+00 1.7765e-01+1.0773e+00j 1.7765e-01-1.0773e+00j >>p=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 >>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为零点位 置列向量。该指令执行后自动生成零极点分布图。 考虑传递函数： 和

传递函数G(s)/H(s)的零极点图如图2-67所示，相应的MATLAB 指令如图2-68所示。 图2-67 零极点图

>>numg=[6 0 1];deng=[1 3 3 1]; >>z=roots(numg) z= 0+0.4082j 0-0.4082j >> p=roots1(deng) p= -1 >>n1=[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= 6s^5+18s^4+25s^3+ 图2-68 绘制零极点图指令

[num,den]= series(num1,den1, num2,den2) 三、结构图模型 一个开环控制系统可以通过G1 (s)与G2 (s)两个环节的串联而得到，利用series()函数可以求串联连接的传递函数，函数的具体形式为： [num,den]= series(num1,den1, num2,den2) 例如G1 (s)和G2 (s)的传递函数分别为： 则

>>num1=[1];den1=[500 0 0]; >>num2=[1 1];den2=[1 2]; 串联函数的用法示于图2-69：   >>num1=[1];den1=[500 0 0]; >>num2=[1 1];den2=[1 2]; >>[num,den]=series(num1,den1,num2,den2); >>printsys(num,den) num/den= s+1 500s^3+1000s^2 图2-69 series函数的用法

[num,den]= cloop (numg,deng, sign) 当系统是以并联的形式连接时，利用parallel()函数可得到系 统的传递函数。指令的具体形式为： 系统以反馈方式构成闭环，则系统的闭环传递函数为： [num,den]= parallel (num1,den1, num2,den2) 求闭环传递函数的MATLAB函数有两个：cloop()和feedback() 其中cloop()函数只能用于H (s)=1（即单位反馈）的情况。 cloop()函数的具体用法为： [num,den]= cloop (numg,deng, sign) 其中numg和deng分别为G (s)的分子和分母多项式，sign=1为正 反馈，sign= -1为负反馈（默认值）。

[num,den]= feedback (numg,deng,numh,denh, sign) 其中numh为H (s)的分子多项式，denh为分母多项式。 闭环反馈系统的结构图如图2-70所示，被控对象G(s)和控制 部分Gc (s)以及测量环节H (s)的传递函数分别为： ，， 图2-70 闭环反馈系统的结构图  

应用series()函数和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= s^2+11s+10 5s^4+60s^3+100s^2+s+1   图2-71 feedback()函数的应用

例2.12 一个多环的反馈系统如图2-49所示，给定各环节的传 递函数为： 试求闭环传递函数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=[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= s^4+ 3s^3+ 3s^2+3s+2 2s^6+38s^5+261s^4+1001s^3+1730s^2+1546s+732 图2-72 多环结构图简化

通过pzmap()或roots()函数可查看传递函数是否有相同的零极 点，还可使用minreal()函数除去传递函数共同的零极点因子。 如图2-73所示。 >>numg=[1 6 11 6];deng=[1 7 12 11 5]; >>printsys(numg,deng) numg/deng= s^3+6s^2+11s+6 s^4+7s^3+12s^2+11s+5 >>[num,den]=minreal(numg,deng); >>printsys(num,den) 1 pole-zeros cancelled num/den= s^2+4s+3 s^3+6s^2+6s+5 图2-73 minreal()函数的应用

例2.2所示的位置随动系统，在给定各元件参数并忽略La和令 ML = 0的情况下，其结构图如图2-74所示： 图2-74 位置随动系统的结构图 第一步求闭环传递函数GB (s)=q c(s) /q r(s)，求解过程及结果 如图2-75所示。第二步利用step()函数计算参考输入q r (t)为单位 阶跃信号时输出q c (t)的响应。

图2-75 位置随动系统的结构图简化及阶跃响应指令 >>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); >>[nc,dc]=series(nb,db,num4,den4); >>[num,den]=cloop(nc,dc,-1); >>printsys(num,den) num/den= 5400 2s^2+2.5s+5400 >>t=[0:0.005:3]; >>[y,t]=step(num,den,t); >>plot(t,y),grid 图2-75 位置随动系统的结构图简化及阶跃响应指令

画出y(t)曲线，grid函数用于给图形加上网格。 图2-76 位置随动系统的阶跃响应曲线 图2-76给出了位置随动系统的阶跃响应曲线。用plot()函数用于 画出y(t)曲线，grid函数用于给图形加上网格。

2.7 循序渐进设计示例：磁盘驱动读取系统 G(s)还可以改写成 根据表2.1磁盘驱动读取系统的典型参数，我们有 我们指出了磁盘驱动系统的基本设计目标：尽可能将磁头准确定位在指定的磁道上，并且磁头从1个磁道转移到另1个磁道所花的时间不超过10ms。在这里，我们将完成设计流程的第4、5步。首先应选定执行机构、传感器和控制器，然后建立控制对象和传感器等元部件的模型。 根据表2.1磁盘驱动读取系统的典型参数，我们有 G(s)还可以改写成

图2-8 磁盘驱动器读取系统框图模型 图2-78 磁盘驱动器读取系统框图模型 表2.1磁盘驱动器读取系统典型参数

其中 。由于 ，因此τ常被略去，有： 或 该闭环系统的框图模型见图2-9。利用框图变换化简规则，有： 利用G(s)的2阶近似表示，可以有：

使用MATLAB的函数step，可以得到 时，如图2-80所示的系统 阶跃响应。 当取Ka=40时，最后可得： 2-79 闭环系统的框图模型 使用MATLAB的函数step，可以得到 时，如图2-80所示的系统 阶跃响应。

小结 利用传递函数研究线性系统，可根据传递函数的极点和零点分布判定系统对不同输入信号的响应特性。 本章讨论了如何建立控制系统以及元部件的数学模型问题。本章所涉及的数学模型共有三种，即微分方程、传递函数、结构图或信号流图。 利用传递函数研究线性系统，可根据传递函数的极点和零点分布判定系统对不同输入信号的响应特性。 利用结构图或信号流图可以了解系统中的每个变量，还可以通过梅逊（Mason）公式，方便地求得系统输入输出间的传递函数。 利用MATLAB软件可求解系统在不同参数和输入情况下的响应。