MATLAB大作业.doc

2008 级MATAB语言与控制系统仿真课程大作业姓 名 学 号 所在院系 电气工程及其自动化 班 级 日 期 2011.1.16 同组人员 作业评分 评阅人 设计报告评分表项 目应包括的主要内容或考核要点常见问题扣分得分基本部分63分方案论证9分性能指标分析；控制方法及实现方案未分析指标没有说明方案理由未进行方案比较有分析指标，未用程序分析有方案说明，未用程序说明-5-1-1-1-1设计过程35分控制器设计与参数计算 有设计过程，未运用程序有简单程序应用，无算法有程序，没有程序说明未说明设计理由-20-5-5-2结果分析10分对设计结果的分析与核算，分析原因和改进没有误差分析没有对结果验证有量化指标，但没有分析有分析验证，没有运用程序没有总结或总结太虚-9-3-2-2-2格式规范9分重点考查完整性，图表，公式的规范性图号图名等问题没有对图表说明中英文参杂，不一致字体不一致图形截屏参考文献规范问题参考文献的引用问题-2-2-1-1-1-1-1提高部分27分第1项9分提出改进的性能指标，完成分析，设计并对结果予以验证未提出其他方案有，未用程序验证有程序，未对程序说明-9-4-2 第2项9分考虑参数变化，干扰影响等其他因素，完成分析，设计并对结果予以验证未对误差干扰进行分析有，未用程序验证有程序，未对程序说明-9-4-2第3项9分提出其他更完善的性能指标，完成分析，设计并对结果予以验证未提出其他指标要求有，未用程序验证有程序，未对程序说明-9-4-2报告得分90分合计特色加分0-30分报告的特色和难度系数，掌握程度予以评价总分报告得分+答辩/特色加分 目录1、 引言.42、 设计方法.53、 结果分析.184、 深入探讨.205、 总结.27六、致谢和参考文献.29七、附录（程序）.30铣床（the milling machine）的系统设计一、引言本学期开始学习MATLAB语言、控制系统分析与设计，其实最早接触matlab是在学习概率论时。在老师的讲解与自己平时的练习中，不断地使用和学习matlab软件。其实一开始学习的目标就不局限于控制系统领域，希望了解该软件的各种用途。当然了，最基本的用途是绘图。经过半个学期的汲取，觉得它用途广泛，功能强大。我最喜欢的就是它的仿真（simulink），也许是对于初学者来说最简单的操作。学习matlab的同时，也同时加深了对控制系统知识的掌握。自控中的绘图可以在matlab下的响应分析的图形用户界面（ltiview）中建立。此外还有sisotool等自控方面的超强工具箱。图10.1铣床模型图10.2铣床系统方框图铣床系统的开环传递函数 （10.1）设计要求：重新设计滞后校正器参数1）对于斜坡输入，稳态误差小于2）使校正后系统的阶跃响应超调量对于（1）参数a不影响系统的设计，在以下分析中作为1处理 本章的核心是滞后校正系统，但之所以对铣床系统采用滞后校正，理由见（四）深入探讨中频域下的分析。另外我们超前校正器、超前-滞后校正器、PI校正器、PD校正器、PID校正器。二、设计方法以运用matlab编程为主，其他各种工具箱为辅助根据极点、零点、增益输入原系统开环传递函数，并加入单位负反馈构成闭环系统程序1G=zpk(,0,-1,-5,2);%原系统开环传递函数sys=feedback(G,1);%单位负反馈，构成闭环得到 （10.2） （10.3）第一个要求稳态误差小于a/8程序2tf1=10;a=1;t=0:0.1:tf1;u=a*t;%斜坡输入y_t=lsim(sys,u,t);plot(t,y_t,b-,t,u,r:)yss=y_t(length(t);y_deta=t(length(t)-yssif y_deta8Kgang=2.02;%K/alphaKvu=2*Kgang/5;alpha=Kvc/Kvu%alpha即z/p得到零点与极点的倍数得到了校正后系统开环传递函数 (10.11)看系统是否符合要求（1）程序8G=zpk(,0,-1,-5,2);%原系统开环传递函数p=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=2.02*tf(nc,dc);G2=G*Gc;%引入K/alphasys2=feedback(G2,1);%单位负反馈，构成闭环系统建立后，输入斜坡信号程序9tf1=10000;a=1;t=0:0.1:tf1;u=a*t;%斜坡输入y_t=lsim(sys2,u,t);plot(t,y_t,b-,t,u,r:)yss=y_t(length(t);y_deta=u(length(t)-yssif y_deta8y=zeros(200,1);i=0;for K=0.1:0.1:2%K/alphaKvu=2*K/5;alpha=Kvc/Kvu;%alpha即z/pp=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=K*tf(nc,dc);%引入K/alphaG2=G*Gc;sys2=feedback(G2,1);%单位负反馈，构成闭环t=0:0.1:19.9;i=i+1;y(:,i)=step(sys2,t);endplot(y)legend(K=0.1,K=0.2,K=0.3,K=0.4,K=0.5,K=0.6,K=0.7,K=0.8,K=0.9,K=1.0,K=1.1,K=1.2,K=1.3,K=1.4,K=1.5,K=1.6,K=1.7,K=1.8,K=1.9,K=2.0)图10.12不同取值下的阶跃响应求各曲线超调量，先将时间范围扩大到200程序14G=zpk(,0,-1,-5,2);%原系统开环传递函数Kvc=10;%Kvcopm8y=zeros(200,1);i=0;for K=0.1:0.1:2%K/alphaKvu=2*K/5;alpha=Kvc/Kvu;%alpha即z/pp=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=K*tf(nc,dc);%引入K/alphaG2=G*Gc;sys2=feedback(G2,1);%单位负反馈，构成闭环t=0:1:199;%改变时间长度i=i+1;y(:,i)=step(sys2,t);endplot(y)legend(K=0.1,K=0.2,K=0.3,K=0.4,K=0.5,K=0.6,K=0.7,K=0.8,K=0.9,K=1.0,K=1.1,K=1.2,K=1.3,K=1.4,K=1.5,K=1.6,K=1.7,K=1.8,K=1.9,K=2.0)图10.13时间轴扩大100观察分布情况程序15sigma=zeros(20,20);j=0;Kgan=zeros(20,20);for i=1:1:20; mp,tf=max(y(:,i);%求最大幅值点 yss=y(length(t);%求稳态值 tp=t(tf);%峰值时间 j=j+1; Kgan(j,1)=j/10;%K/alpha sigma(j,1)=100*(mp-yss)/yss;%超调量endplot(Kgan,sigma,0 2,20 20)gridxlabel(K/alpha)ylabel(Overshoot)图10.14与的关系若找到的取值范围，问题就比较清晰程序16n=find(sigma(:,1)8y=zeros(200,1);i=0;y_deta=zeros(20,1);Kgan=zeros(20,1);for K=0.1:0.1:2%K/alphaKvu=2*K/5;alpha=Kvc/Kvu;%alpha即z/pp=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=K*tf(nc,dc);%引入K/alphaG2=G*Gc;sys2=feedback(G2,1);%单位负反馈，构成闭环t=0:1:199;i=i+1;Kgan(i,1)=i/10;u=t;y(:,i)=lsim(sys2,u,t);yss=y(length(t),i);y_deta(i,1)=u(length(t)-yssif y_deta(i,1)(1/8)%稳态误差小于a/8 disp(yes)%符合要求（1）else disp(no)%不符合endendplot(Kgan,y_deta,0 2,0.125 0.125)%找到误差小于a/8的直线gridxlabel(K/alpha)ylabel(E(s)图10.15与稳态误差的关系程序18m=find(y_deta8Kgang=1;%K/alphaKvu=2*Kgang/5;alpha=Kvc/Kvu;%alpha即z/p将参数带入校正后系统程序20G=zpk(,0,-1,-5,2);%原系统开环传递函数p=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=1*tf(nc,dc);G2=G*Gc;%引入K/alphasys2=feedback(G2,1);%单位负反馈，构成闭环求单位斜坡响应程序21tf1=1000;a=1;t=0:0.1:tf1;u=a*t;%斜坡输入y_t=lsim(sys2,u,t);plot(t,y_t,b-,t,u,r:)gridyss=y_t(length(t);y_deta=u(length(t)-yss;if y_deta8Kgang=0.5;%K/alphaKvu=2*Kgang/5;alpha=Kvc/Kvu;%alpha即z/pG=zpk(,0,-1,-5,2);%原系统开环传递函数p=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=Kgang*tf(nc,dc);G2=G*Gc;%引入K/alphasys2=feedback(G2,1);%单位负反馈，构成闭环tf1=1000;a=1;t=0:0.1:tf1;u=a*t;%斜坡输入y_t=lsim(sys2,u,t);plot(t,y_t,b-,t,u,r:)gridyss=y_t(length(t);y_deta=u(length(t)-yss;if y_deta=-135);w1=w(max(m)得到当，相位等于45deg频率穿越点系统已经满足要求二，而且有23.5dB的增益裕度若利用增益K来调节系统使其满足要求： (10.14)对于要求一： (10.15)得到K=20程序28G=zpk(,0,-1,-5,2);%原系统开环传递函数bode(G);%原系统伯德图hold onbode(20*G)%引入增益Kgridlegend(K=1,K=20)图10.19增益校正前后系统伯德图虽然K=20满足条件一可是相位裕度为负，闭环系统不稳定所以不能用使用单一的增益校正。由于要求一等价于而且 (10.16)于是将要求一等价于：相位裕度在引入的条件下，设之后校正器 (10.17)采用滞后校正的理由：因为原系统引入增益后，相位裕度不足，但低频下的相位均大于-135deg，所以有充足的裕度来实现滞后校正。相位滞后校正是通过衰减校正前系统伯德图的幅值部分来减小增益穿越频率，这样可以使系统产生必要的相位裕度求滞后校正系统的步骤及原理(10.18) (10.19)取得校正系统最小相位 图10.20相位滞后网络的伯德图步骤：(1) 根据稳态误差要求确定开环增益K。绘制未校正系统的伯德图，并求出其相位裕量和增益裕量。(2) 确定校正后系统的增益剪切频率，即要求的相位裕度下的频率，校正后系统此处应该幅值为0dB(3) (3) 求值。确定原系统频率特性在处幅值下降到0dB时所必需的衰减量M。由等式 (10.20)(4) 选取T值。为了使滞后校正装置产生的相位滞后对校正后系统的增益剪切频率处的影响足够小，应满足，取 (10.21)程序29Ko=20;G=zpk(,0,-1,-5,2*Ko);%原系统开环传递函数sys=tf(G);mag,phase,w=bode(sys)figure(1);margin(sys)gridhold ongama=50;%频域下要求二pha=gama+5-180;wgc=spline(phase,w,pha);%找校正后的频率穿越点na=polyval(sys.num1,j*wgc);da=polyval(sys.den1,j*wgc);g=na/da;g1=abs(g);h=20*log10(g1);%h的取值等于20lg（beta），才能使幅值等于0beta=10(h/20);%求betaT=10/wgc;betat=beta*T;Gc=tf(T 1,betat 1)%滞后校正系统函数bode(G*Gc)legend(G,G*Gc)图10.21滞后校正前后系统伯德图得到校正滞后系统传递函数 不含增益 (10.22)程序30sys2=feedback(Gc*sys,1);%校正后闭环传递函数step(sys2)%求阶跃响应图10.22滞后校正后系统阶跃响应没有达到要求二误差分析： 1）编程中很多赋值只是取得一些分散的点，因而求取参数时存在误差 2）预计相位裕度虽然理论足够，但相位滞后校正器的最大相位不为零，余量5可能不足改变相位裕度余量程序31Ko=20;G=zpk(,0,-1,-5,2*Ko);%原系统开环传递函数sys=tf(G);mag,phase,w=bode(sys)figure(1);margin(sys)gridhold ongama=50;%频域下要求二pha=gama+10-180;wgc=spline(phase,w,pha);%找校正后的频率穿越点na=polyval(sys.num1,j*wgc);da=polyval(sys.den1,j*wgc);g=na/da;g1=abs(g);h=20*log10(g1);%h的取值等于20lg（beta），才能使幅值等于0beta=10(h/20);%求betaT=10/wgc;betat=beta*T;Gc=tf(T 1,betat 1)%滞后校正系统函数bode(G*Gc)legend(G,G*Gc) (10.23)图10.23修改后系统的阶跃响应图10.24修正后系统伯德图加入增益后传递函数为 (10.24)与式（10.12）（10.13）相比较，其介于两者之间。又因为可以去0.5和1.0之间的值，则时域和频域下的结果可以得到统一5、 总结（暨心得）本次设计是结合了英文教材Binder10和课件Chapter10来设计的，两者的方向截然不同。英文教材虽然提及滞后校正，但其完全是基于时域下，用理论推导出传递函数，似乎只有在运用根轨迹处依赖了电脑操作，而课件里侧重于频域下编程。所以设计分析是时域，深入探讨是频域，结合着线性控制的理论。两种方法的优缺点。频域下的滞后校正也有很多不严格的地方，一些参数求取的近似。但其更加直观和易于操作。时域下省略了转化要求参数的步骤，但时域的下的高阶系统分析比较繁琐，若采用近似模型误差很大。铣床系统是一型系统，在斜坡输入下的稳态误差不为零。是否工程上一定要求非零的稳态误差。如果没有此类要求是否可以通过改善系统的型来消除稳态误差，更加符合要求一。当然了这一章的侧重点是滞后校正系统的设计（Lag Compensator Design）。这也与该系统的频率特性有关，充足的相位可以采用滞后校正。然而本次设计并没有考虑，滞后校正对于频带带宽的影响，如果转折频率1/T和1/T选取的过小，会造成带宽过小，是系统响应变得缓慢。那么由滞后校正引起的响应时间的增加，应如何再不影响斜坡输入稳态误差和阶跃响应超调量的前提下，得到改善。我分别从时域和频域下求取了滞后校正器，是否可以利用状态空间的方法。如果利用状态空间的话，要求又如何转变。学习matlab还是觉得只是纯粹的数学问题，虽然每一章都有实例，但最根本的还是围绕着传递函数来解决问题的。不理解系统原理的话，仍然可以掌握每一章节的精髓。个人比较喜欢sisotool，虽然用的很少，但对于根轨迹和零极点配置的理解和掌握有很大帮助。这次也用到了借助根轨迹来限制阻尼比，也使得自控上的这方面知识得到应用。六、致谢和参考文献线性控制系统工程（Linear Control Systems Engineering）七、附录（程序）程序1：G=zpk(,0,-1,-5,2);%原系统开环传递函数sys=feedback(G,1);%单位负反馈，构成闭环程序2tf1=10;a=1;t=0:0.1:tf1;u=a*t;%斜坡输入y_t=lsim(sys,u,t);plot(t,y_t,b-,t,u,r:)yss=y_t(length(t);y_deta=t(length(t)-yssif y_deta8Kgang=2.02;%K/alphaKvu=2*Kgang/5;alpha=Kvc/Kvu%alpha即z/p程序8G=zpk(,0,-1,-5,2);%原系统开环传递函数p=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=2.02*tf(nc,dc);G2=G*Gc;%引入K/alphasys2=feedback(G2,1);%单位负反馈，构成闭环程序9tf1=10000;a=1;t=0:0.1:tf1;u=a*t;%斜坡输入y_t=lsim(sys2,u,t);plot(t,y_t,b-,t,u,r:)yss=y_t(length(t);y_deta=u(length(t)-yssif y_deta8y=zeros(200,1);i=0;for K=0.1:0.1:2%K/alphaKvu=2*K/5;alpha=Kvc/Kvu;%alpha即z/pp=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=K*tf(nc,dc);%引入K/alphaG2=G*Gc;sys2=feedback(G2,1);%单位负反馈，构成闭环t=0:0.1:19.9;i=i+1;y(:,i)=step(sys2,t);endplot(y)legend(K=0.1,K=0.2,K=0.3,K=0.4,K=0.5,K=0.6,K=0.7,K=0.8,K=0.9,K=1.0,K=1.1,K=1.2,K=1.3,K=1.4,K=1.5,K=1.6,K=1.7,K=1.8,K=1.9,K=2.0)程序14G=zpk(,0,-1,-5,2);%原系统开环传递函数Kvc=10;%Kvcopm8y=zeros(200,1);i=0;for K=0.1:0.1:2%K/alphaKvu=2*K/5;alpha=Kvc/Kvu;%alpha即z/pp=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=K*tf(nc,dc);%引入K/alphaG2=G*Gc;sys2=feedback(G2,1);%单位负反馈，构成闭环t=0:1:199;%改变时间长度i=i+1;y(:,i)=step(sys2,t);endplot(y)legend(K=0.1,K=0.2,K=0.3,K=0.4,K=0.5,K=0.6,K=0.7,K=0.8,K=0.9,K=1.0,K=1.1,K=1.2,K=1.3,K=1.4,K=1.5,K=1.6,K=1.7,K=1.8,K=1.9,K=2.0)程序15sigma=zeros(20,20);j=0;Kgan=zeros(20,20);for i=1:1:20; mp,tf=max(y(:,i);%求最大幅值点 yss=y(length(t);%求稳态值 tp=t(tf);%峰值时间 j=j+1; Kgan(j,1)=j/10;%K/alpha sigma(j,1)=100*(mp-yss)/yss;%超调量endplot(Kgan,sigma,0 2,20 20)gridxlabel(K/alpha)ylabel(Overshoot)程序16n=find(sigma(:,1)8y=zeros(200,1);i=0;y_deta=zeros(20,1);Kgan=zeros(20,1);for K=0.1:0.1:2%K/alphaKvu=2*K/5;alpha=Kvc/Kvu;%alpha即z/pp=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=K*tf(nc,dc);%引入K/alphaG2=G*Gc;sys2=feedback(G2,1);%单位负反馈，构成闭环t=0:1:199;i=i+1;Kgan(i,1)=i/10;u=t;y(:,i)=lsim(sys2,u,t);yss=y(length(t),i);y_deta(i,1)=u(length(t)-yssif y_deta(i,1)(1/8)%稳态误差小于a/8 disp(yes)%符合要求（1）else disp(no)%不符合endendplot(Kgan,y_deta,0 2,0.125 0.125)%找到误差小于a/8的直线gridxlabel(K/alpha)ylabel(E(s)程序18m=find(y_deta8Kgang=1;%K/alphaKvu=2*Kgang/5;alpha=Kvc/Kvu;%alpha即z/p程序20G=zpk(,0,-1,-5,2);%原系统开环传递函数p=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=1*tf(nc,dc);G2=G*Gc;%引入K/alphasys2=feedback(G2,1);%单位负反馈，构成闭环程序21tf1=1000;a=1;t=0:0.1:tf1;u=a*t;%斜坡输入y_t=lsim(sys2,u,t);plot(t,y_t,b-,t,u,r:)gridyss=y_t(length(t);y_deta=u(length(t)-yss;if y_deta8Kgang=0.5;%K/alphaKvu=2*Kgang/5;alpha=Kvc/Kvu;%alpha即z/pG=zpk(,0,-1,-5,2);%原系统开环传递函数p=0.001;z=alpha*p;nc=1 z;dc=1,p;Gc=Kgang*tf(nc,dc);G2=G*Gc;%引入K/alphasys2=feedback(G2,1);%单位负反馈，构成闭环tf1=1000;a=1;t=0:0.1:tf1;u=a*t;%斜坡输入y_t=lsim(sys2,u,t);plot(