MATLAB控制系统仿真

MATLAB控制系统仿真提纲一、弹簧-重物-阻尼器系统二、传递函数三、结构图模型引言MATLAB是一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算和图形显示于一体，构成了一个方便的界面友好的用户环境。控制系统的分析与设计方法，都是以数学模型为基础进行的。MATLAB可以用于以传递函数形式描述的控制系统。在本节中，首先举例说明如何使用MATLAB进行辅助分析。然后讨论传递函数和结构图。一、弹簧-重物-阻尼器系统弹簧—重物—阻尼器动力学系统如图2-1所示。重物M的位移由y(t)表示，用微分方程描述如下：22d()d()()()ddytytMBKytfttt该系统在初始位移作用下的瞬态响应为：22(0)()sin(1)1ntnyytet其中cos-1，初始位移是y(0)。系统的瞬态响应当＜1时为欠阻尼，当＞１时为过阻尼，当＝1时为临界阻尼。过阻尼情况：y(0)=0.15mn=(弧度/秒)（）欠阻尼情况：y(0)=0.15mn=(弧度/秒)（）213222,3KBMM221222,1KBMM利用MATLAB程序—unforced.m，可以显示初始位移为y(0)的物体自由运动曲线，如图2-63所示。在unforced.m程序中，变量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-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-63分析弹簧—重物—阻尼器的MATLAB指令图2-64弹簧—重物—阻尼器的自由运动曲线在欠阻尼和过阻尼情况下的响应曲线如图2-64所示：MATLAB可分析以传递函数形式描述的系统。分子多项式和分母多项式都必须在MATLAB指令中指定。在MATLAB中多项式由行向量组成，这些行向量包含了降次排列的多项式系数。例如多项式p(s)=1s3+3s2+0s1+4s0，按图2-65的格式输入p=[1304]，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()。函数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()函数设传递函数为G(s)=num/den，其中num和den均为多项式。利用函数：二、传递函数[P,Z]=pzmap(num,den)可得G(s)的零极点位置，即P为极点位置列向量，Z为零点位置列向量。该指令执行后自动生成零极点分布图。考虑传递函数：23261()331sGssss和(1)(2)()(2)(2)(3)ssHssisis图2-67零极点图传递函数G(s)/H(s)的零极点图如图2-67所示，相应的MATLAB指令如图2-68所示。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+图2-68绘制零极点图指令三、结构图模型一个开环控制系统可以通过G1(s)与G2(s)两个环节的串联而得到，利用series()函数可以求串联连接的传递函数，函数的具体形式为：[num,den]=series(num1,den1,num2,den2)例如G1(s)和G2(s)的传递函数分别为：11()2sGss221()500Gss12232()111()()()25005001000CsssGsGsRsssss则串联函数的用法示于图2-69：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()函数可得到系统的传递函数。指令的具体形式为：系统以反馈方式构成闭环，则系统的闭环传递函数为：()()1()()BGsGsGsHs[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为负反馈（默认值）。feedback()函数的用法为：[num,den]=feedback(numg,deng,numh,denh,sign)其中numh为H(s)的分子多项式，denh为分母多项式。闭环反馈系统的结构图如图2-70所示，被控对象G(s)和控制部分Gc(s)以及测量环节H(s)的传递函数分别为：，，1()2csGssnumcdenc21()5Gssnumgdeng1()10Hssnumhdenh图2-70闭环反馈系统的结构图应用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所示，给定各环节的传递函数为：11()10Gss21()1Gss2321()44sGsss41()6sGss1()1Hs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=[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()函数可查看传递函数是否有相同的零极点，还可使用minreal()函数除去传递函数共同的零极点因子。如图2-73所示。numg=[16116];deng=[1712115];printsys(numg,deng)numg/deng=s^3+6s^2+11s+6s^4+7s^3+12s^2+11s+5[num,den]=minreal(numg,deng);printsys(num,den)1pole-zeroscancellednum/den=s^2+4s+3s^3+6s^2+6s+5图2-73minreal()函数的应用例2.2所示的位置随动系统，在给定各元件参数并忽略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=[20.50];num3=[0.20];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=54002s^2+2.5s+5400t=[0:0.005:3];[y,t]=step(num,den,t);plot(t,y),grid图2-75位置随动