3 MATLAB控制系统仿真(讲的很多)

2.5应用MATLAB控制系统仿真MATLAB是一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算和图形显示于一体，构成了一个方便的界面友好的用户环境。由控制领域专家推出的MATLAB工具箱之一的控制系统(ControlSystem)，在控制系统计算机辅助分析与设计方面获得了广泛的应用，并且MATLAB工具箱的内容还在不断增加，应用范围也越来越宽。基础进行的。MATLAB可以用于以传递函数形式描述的控制系统。在本节中，我们首先以一个典型的动力学系统“弹簧—重物—阻尼器”的数学模型为例，说明如何使用MATLAB进行辅助分析。之后，我们讨论传递函数和结构图。特别的，主要介绍以下内容：如何使用MATLAB求解多项式，计算传递函数的零点和极点，计算闭环传递函数，计算结构图的等效变换以及闭环系统对单位阶跃输入的响应等。MATLAB函数有：roots，roots1，series，parallel，feedback，cloop，poly，conv，polyval，printsys，pzmap和step等。2.5.1.举例一个弹簧—重物—阻尼器动力学系统如图2-1所示。重物M的位移由y(t)表示，用微分方程描述如下：该系统在初始位移作用下的瞬态响应为其中cos-1，初始位移是y(0)。系统的瞬态响应当＜1时为欠阻尼，当＞１时为过阻尼，当＝1时为临界阻尼。考虑过阻尼情况：y(0)=0.15mn=(弧度/秒)()欠阻尼情况：y(0)=0.15mn=(弧度/秒)()利用MATLAB程序—unforced.m，可以显示初始位移为y(0)的物体自由运动曲线，如图2-63所示。在unforced.m程序中，变量y(0)，n，t，1和2的值由指令直接输入工作区，然后运行unforced.m程序就可以产生响应曲线。MATLAB提供了一种交互式的应用环境，我们可以在指令窗口随时修改n，1和2值并再次运行unforced.m程序。在上述欠阻尼和过阻尼情况下的响应曲线如图2-64所示。可以看到程序运行后阻尼系数的值被标注在不同的曲线上，以防止混淆。对于弹簧—重物—阻尼器系统，利用MATLAB求微分方程的解是容易而有效的。一般说来，当分析一个具有多样性的输入和初始条件以及元件参数的闭环反馈控制系统时，欲直接得到这些因素与系统响应的关系是比较困难的。在这种情况下，可以利用MATLAB非常有效的数值计算能力，反复运行MATLAB程序并绘出曲线图，就可以得出结论。y0=0.15;wn=sqrt(2);zeta1=3/(2*sqrt(2）;zeta2=1/(2*sqrt(2）;t=[0:0.1:10];unforceda）MATLAB指令窗口%unforced.m%计算系统在给定初始条件下的自由运动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,‘--’),gridxlabel(‘Time[sec]’),ylabel(‘y(t)Displacement[m]’)text(0.2,0.85,[‘oeverdampedzeta1=’,num2str(zeta1),])text(0.2,0.80,[‘underdampedzeta2=’,num2str(zeta2),])（b）分析弹簧—重物—阻尼器的MATLAB程序unforced.m图2-64弹簧—重物—阻尼器的自由运动曲线MATLAB可以用来分析以传递函数形式描述的系统。由于传递函数是多项式的比值，所以如何处理多项式是MATLAB首先需要解决的问题。记住，分子多项式和分母多项式都必须在MATLAB指令中指定。MATLAB中多项式由行向量组成，而这些行向量包含了降次排列的多项式系数。例如多项式p(s)=1s3+3s2+0s1+4s0，按图2-65的格式输入p=[1304]，请注意，尽管s的系数为0，它也一定要包含在确定p(s)的行向量中。p是一个包含降幂排列的p(s)系数的行向量，那么roots(p)是一个包含多项式根的列向量。相反的，如果r是一个包含多项式根的列向量，那么，poly(r)是一个包含降幂排列多项式系数的行向量，见图2-65。我们可以用上述的roots()函数计算多项式p(s)的根，但是当多项式有重根时，函数roots1()能给出更精确的结果。p=[1304];r=roots(p)r=-3.3553e+001.7765e-01+1.0773e+00j1.7765e-01-1.0773e+00jp=poly(r)p=1.0003.0000.000-0.000j4.000+0.000j图2-65输入多项式并求根MATLAB的conv()函数完成。假设我们要把两个多项式相乘合并成一个多项式n(s)，即n(s)=(3s2+2s+1)(s+4)=3s3+14s2+9s+4与此运算相关的MATLAB函数就是conv()，如图2-66所示。函数polyval()用来计算多项式的值。多项式n(s)在s=-5处的值为n(-5)=-66，见图2-66。p=[321];q=[14];n=conv(p,q)n=31494value=polyval(n,-5)value=-66图2-66MATLAB的conv()函数和polyval()函数2.5.2.传递函数MATLAB函数获得传递函数在复平面内的零极点分布图。设传递函数为G(s)=num/den，其中num和den均为多项式。利用函数[P,Z]=pzmap(num,den)可以获得G(s)的零极点位置，即P为极点位置列向量，Z为零点位置列向量。该指令执行后自动生成零极点分布图。考虑传递函数和利用一系列MATLAB指令和函数，我们可以计算传递函数的零极点、特征方程和两个传递函数相除，还可以在复平面上获得G(s)/H(s)的零极点分布图。传递函数G(s)/H(s)的零极点图如图2-67所示，相应的MATLAB指令如图2-68所示。图2-67中清楚的表示出了5个零点的分布，但是看上去只有两个极点，实际上这是不可能的，因为我们知道极点的个数一定大于或等于零点的个数。应用函数roots1()可以看到，实际上有四个极点位于s=-1。因此，在同一位置上重复的极点或零点在零极点分布图上是区分不出的。2-67零极点图2.5.3.结构图模型MATLAB函数将这些部分联接起来构成一个闭环控制系统，实现结构图的转换，计算从输入R(s)到输出C(s)的传递函数。为了方便介绍MATLAB函数，我们从结构图的基本变换开始。一个简单的开环控制系统可以通过G1(s)与G2(s)两个环节的串联而得到，利用series()函数可以求串联连接的传递函数，函数的具体形式为[num,den]=series(num1,den1,num2,den2)例如G1(s)和G2(s)的传递函数分别为，则串联函数的用法示于图2-69。numg=[601];deng=[1331];z=roots(numg)z=0+0.4082j0-0.4082jp=roots1(deng)p=-1-1-1n1=[11];n2=[12];d1=[12*j];d2=[1–2*j];d3=[13];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+75s^2+4s+12ａs^5+6s^4+14s^3+16s^2+9s+2pzmap(num,den)title(‘Pole-ZeroMap’)图2-68绘制零极点图指令num1=[1];den1=[50000];num2=[11];den2=[12];[num,den]=series(num1,den1,num2,den2);printsys(num,den)num/den=s+1500s^3+1000s^2图2-69series函数的用法，利用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为负反馈（默认值）。2-70闭环反馈系统的结构图feedback()函数的具体用法为[num,den]=feedback(numg,deng,numh,denh,sign)其中numh为H(s)的分子多项式，denh为分母多项式。假设闭环反馈系统的结构图如图2-70所示，被控对象G(s)和控制部分Gc(s)以及测量环节H(s)的传递函数分别为，，应用series()函数和feedback()函数求解闭环传递函数的MATLAB指令如图2-71所示。numg=[1];deng=[500];numc=[11];denc=[12];numh=[1];denh=[110];[num1,den1]=series(numc,denc,numg,deng);[num,den]=feedback(num1,den1,numh,denh,-1);printsys(num,den)num/den=s^2+11s+105s^4+60s^3+100s^2+s+1图2-71feedback()函数的应用例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=[110];ng2=[1];dg2=[11];ng3=[101];dg3=[144];ng4=[11];dg4=[16];nh1=[1];dh1=[1];nh2=[2];dh2=[1];nh3=[11];dh3=[12];[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+22s^6+38s^5+261s^4+1001s^3+1730s^2+1546s+732图2-72多环结构图简化pzmap()或roots()函