自动控制原理课程设计MATLAB仿真.doc

目录概述1一、实验目的1二、简述MATLAB语言的特点及其主要功能1三、控制系统仿真时常用的方法和指令21、控制系统仿真时常用的方法2a、数学仿真2b、半物理仿真2c、全物理仿真22、控制系统仿真时常用的指令21)、Bode图2、绘制Bode图2、系统的增益裕度和相角裕度22）、Nichols图33）、Nyquist图34）、一般频率响应图35）、频率响应的奇异值图36）、绘制根轨迹4四、实验内容4五、心得体会22六、参考文献22概述MATLAB 是一种直观、高效的计算机语言，同时也是一个科学计算平台。它的伴随工具Simulink 是用来对真实世界的动力学系统建模、模拟仿真和分析的软件。我们可将综合性和设计性实验项目通过MATLAB 在计算机上仿真，使系统的观察实验的动态过程。目前，MATLAB 已经成为我们当代大学生必须掌握的基本技能，在设计研究单位和工业部门，MATLAB 已经成为研究和解决各种具体工程问题的一种标准软件。在完成了验证性、综合性和设计性实验后，课程设计必不可少。课程设计是工科实践教学的一个重要的环节，目的是培养我们综合运用理论知识分析和解决实际问题的方法和能力,实现由知识向技能的初步化。所以课程设计是培养我们思维创造能力最有效的途径。一、实验目的1、培养理论联系实际的科学态度，训练综合运用经典控制理论和相关课程知识的能力。2、掌握自动控制原理的时域分析法，根轨迹法，频域分析法，以及各种（矫正）装置的作用及用法，能够利用不同的分析方法对给定系统进行性能分析，能根据不同的系统性能指标要求进行合理的系统设计，并调试满足系统的指标。3、学会使用MATLAB语言及Simulink动态仿真工具进行系统的仿真与调试。4、锻炼独立思考和动手解决控制系统实际问题的能力。二、简述MATLAB语言的特点及其主要功能MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连 matlab开发工作界面接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。三、控制系统仿真时常用的方法和指令1、控制系统仿真时常用的方法a、数学仿真 也称计算机仿真，就是在计算机上实现描写系统物理过程的数学模型，并在这个模型上对系统进行定量的研究和实验。这种仿真方法常用于系统的方案设计阶段和某些不适合做实物仿真的场合（包括某些故障模式）。它的特点是重复性好、精度高、灵活性大、使用方便、成本较低、可以是实时的、也可以是非实时的。数学仿真的逼真度和精度取决于仿真计算机的精度和数学模型的正确性与精确性。数学仿真可采用模拟计算机、数字计算机和数字-模拟混合计算机。b、半物理仿真即采用部分物理模型和部分数学模型的仿真。其中物理模型采用控制系统中的实物，系统本身的动态过程则采用数学模型。半物理仿真系统通常由满足实时性要求的仿真计算机、运动模拟器（一般采用三轴机械转台）、目标模拟器、控制台和部分实物组成。c、全物理仿真 即全部采用物理模型的仿真，又称实物模拟。全物理仿真技术复杂，一般只在必要时才采用。2、控制系统仿真时常用的指令1)、Bode图、绘制Bode图bode函数绘制控制系统幅频和相频图调用格式 bode(sys) bode(sys, w) bode(sys1,sys2.sysn)bode(sys1,sys2.sysn,w) bode(sys1,plotstyle1,sys2,plotstyle2,.sysn,plotstylen,)mag, phase, w=bode(sys)w频率区间矢量，控制频率起止范围mag, phase, w输出幅值矢量 输出相角 输出频率区间绘制系统幅频图调用格式 bodemag(sys) bodemag(sys, w) bodemag(sys1,sys2.sysn) bodemag(sys1,sys2.sysn,w) bodemag(sys1,plotstyle1,sys2,plotstyle2,.sysn,plotstylen,)、系统的增益裕度和相角裕度增益裕度定义：正好使系统频率响应穿过临界点的增益与实际系统增益的比值。相角裕度：为了使系统的频率响应通过临界点必须引入的纯相位滞后量。调用格式 gm pm wcg wcp=margin(sys) gm pm wcg wcp =margin(mag, phase, w)Gm系统增益裕度 pm系统相角裕度wcg wcp-表示交叉频率2）、Nichols图对数幅相特性图(Nichols图)是描述系统频率特性的图示方法。该图纵坐标表示频率特性的对数幅值，以分贝为单位；横坐标表示频率特性的相位角。调用格式 nichols(sys) nichols(sys, w) nichols(sys1,sys2.sysn) nichols(sys1,sys2.sysn,w) nichols(sys1,plotstyle1,sys2,plotstyle2,.sysn,plotstylen,)mag, phase, w=nichols(sys)mag, phase=nichols(sys, w)3）、Nyquist图系统频率特性的极坐标图Nyquist图调用格式 nyquist(sys) nyquist(sys, w) nyquist(sys1,sys2.sysn) nyquist(sys1,sys2.sysn,w) nyquist(sys1,plotstyle1,sys2,plotstyle2,.sysn,plotstylen,)mag, phase, w=nyquist(sys)mag, phase=nyquist(sys, w)4）、一般频率响应图调用格式 h=freqs(num , den, w)h ,w=freqs(num, den)h ,w=freqs(num ,den, f) freqs(num ,den)Num传递函数分子矢量 den传递函数分母矢量W频率区间矢量5）、频率响应的奇异值图调用格式 sigma(sys) sigma(sys, w) sigma(sys ,w, type)sigma(sys1,sys2.sysn) sigma(sys1,sys2.sysn,w)sigma(sys1,sys2.sysn,w, type)sigma(sys1,plotstyle1,sys2,plotstyle2,.sysn,plotstylen,)sv, w=sigma(sys) sv=sigma(sys, w)Type绘制奇异值图的类型Type=1 控制系统的频率响应H的奇异值图Type=2 控制系统的频率响应I+H的奇异值图Type=3 控制系统的频率响应I+H-1的奇异值图Sv频率的奇异值 plotstyle绘制奇异值曲线的特性6）、绘制根轨迹、求系统的零点、极点、增益调用格式 p=pole(sys) 计算控制系统的极点 Z=zero(sys) 计算控制系统的零点 z, gain=zero(sys) 计算控制系统的零点、增益、绘制系统极点、零点图调用格式 pzmap(sys) pzmap(sys1,sys2,sysn) p z=pzmap(sys)、绘制根轨迹图调用格式 rlocus(sys) rlocus(sys ,k)Rlocus(sys1,sys2,)r k=rlocus(sys) r=rlocus(sys,k)k表示增益 r系统的极点23四、实验内容1.用matlab语言编制程序，实现以下系统：【实验原理】：线性系统的传递函数模型可一般地表示为： 将系统的分子和分母多项式的系数按降幂的方式以向量的形式输入给两个变量和，就可以轻易地将传递函数模型输入到MATLAB环境中。命令格式为：; ;在MATLAB控制系统工具箱中，定义了tf() 函数，它可由传递函数分子分母给出的变量构造出单个的传递函数对象。从而使得系统模型的输入和处理更加方便。 该函数的调用格式为： Gtf(num，den);将系统增益、零点和极点以向量的形式输入给三个变量、Z和P，就可以将系统的零极点模型输入到MATLAB工作空间中，命令格式为： 在MATLAB控制工具箱中，定义了zpk()函数，由它可通过以上三个MATLAB变量构造出零极点对象，用于简单地表述零极点模型。该函数的调用格式为：G=zpk(Z,P,KGain) 1） 【程序】 num=5,24,0,18;den=1,4,6,2,2;sys=tf(num,den)Transfer function: 5 s3 + 24 s2 + 18-s4 + 4 s3 + 6 s2 + 2 s + 22）【程序】 A=1,6,6;roots(A) 注释：求解的根ans = -4.7321 -1.2679 B=1,3,2,5;roots(B) 注释：求解的根ans = -2.9042 -0.0479 + 1.3112i -0.0479 - 1.3112iz=-2,-4.7321,-4.7321,-1.2679,-1.2679;p=0,-1,-1,-1,-2.9042,-0.0479+ 1.3112i,-0.0479 - 1.3112i;k=4;sys=zpk(z,p,k) Zero/pole/gain: 4 (s+2) (s+1.268)2 (s+4.732)2-s (s+1)3 (s+2.904) (s2 + 0.0958s + 1.722)2.两环节G1、G2串联，求等效的整体传递函数G(s)【实验原理】：若假定两环节均为单输入单输出的系统SA和SB。两个环节级联：sysseries(SA，SB) 方法一：【程序】 G1=tf(2,1,3);G2=tf(7,1,2,1);sys=G1*G2 Transfer function: 14-s3 + 5 s2 + 7 s + 3方法二：【程序】 G1=tf(2,1,3);G2=tf(7,1,2,1);sys=series(G1,G2) Transfer function: 14-s3 + 5 s2 + 7 s + 33.两环节G1、G2并联，求等效的整体传递函数G(s)【实验原理】：若假定两环节均为单输入单输出的系统SA和SB两个环节并联：sys=parallel(SA，SB) 方法一：【程序】 G1=tf(2,1,3);G2=tf(7,1,2,1);sys=G1+G2 nsfer function: 2 s2 + 11 s + 23-s3 + 5 s2 + 7 s + 3方法二：【程序】 G1=tf(2,1,3);G2=tf(7,1,2,1);sys=parallel(G1,G2) Transfer function: 2 s2 + 11 s + 23-s3 + 5 s2 + 7 s + 34.已知系统结构如图，求闭环传递函数。其中的两环节G1、G2分别为【实验原理】：若假定两环节均为单输入单输出的系统SA和SB。A环节前向，B环节反馈：S=feedback(SA,SB)负反馈 【程序】G1=tf(3,100,1,2,81);G2=tf(2,2,5);sys=feedback(G1,G2,-1) Transfer function: 6 s2 + 215 s + 500-2 s3 + 9 s2 + 178 s + 605正反馈 【程序】 G1=tf(3,100,1,2,81);G2=tf(2,2,5);sys=feedback(G1,G2,+1) Transfer function: 6 s2 + 215 s + 500-2 s3 + 9 s2 + 166 s + 2055.已知某闭环系统的传递函数为，求其单位阶跃响应曲线，单位脉冲响应曲线。【实验原理】：LTI模型的阶跃响应函数step( )格式：step(sys)功能：绘制系统sys(sys由函数tf、zpk或ss产生)的阶跃响应，结果不返回数据，只返回图形。对多输入多输出模型，将自动求每一输入的阶跃响应。LTI模型的单位冲激响应函数impulse( )格式：impulse(sys)功能：绘制系统sys（sys由函数tf、zpk或ss产生）的单位冲激响应，结果不返回数据，只返回图形。（1）、单位阶跃响应【程序】G=tf(10,25,0.16,1.96,10,25)Transfer function: 10 s + 25-0.16 s3 + 1.96 s2 + 10 s + 25step(G); %计算并绘制系统的单位阶跃响应title(单位阶跃响应)（2）、单位脉冲响应 【程序】impulse(G); %计算并绘制系统的单位脉冲响应title(单位脉冲响应)6.典型二阶系统的传递函数为 ， 为自然频率， 为阻尼比，试绘出当=0.5，分别取0、2、4、6、8、10时该系统的单位阶跃响应曲线；分析阻尼比分别为0.5、1时系统的稳定性。【实验原理】：=0.5时，利用阶跃响应函数step( )绘制分别取0、2、4、6、8、10时的阶跃响应曲线。(1)、=0.5，分别取0、2、4、6、8、10时该系统的单位阶跃响应曲线【程序】 G1=tf(0,4,1,2,4); %=0.5 Wn=2 G2=tf(0,16,1,4,16); %=0.5 Wn=4 G3=tf(0,36,1,6,36); %=0.5 Wn=6 G4=tf(0,64,1,8,64); %=0.5 Wn=8 G5=tf(0,100,1,10,100); %=0.5 Wn=10 step(G1,G2,G3,G4,G5); %计算并绘制系统的单位阶跃响应title(单位阶跃响应)（2）、=-0.5,=-1时的阶跃响应曲线 【程序】 G6=tf(0,100,1,-10,100);%=-0.5,Wn=10时的单位阶跃响应 G7=tf(0,100,1,-20,100);%=-1,Wn=10时的单位阶跃响应 step(G6,G7);title(=-0.5，=-1时的阶跃响应)分析稳定性： 从=-0.5，=-1时的阶跃响应曲线图可以看出，曲线不收敛，而是发散的，说明系统是不稳定的。7. 设有一高阶系统开环传递函数为，试绘制该系统的零极点图和闭环根轨迹图。【实验原理】：、求系统的零点、极点、增益调用格式 p=pole(sys) 计算控制系统的极点 Z=zero(sys) 计算控制系统的零点 z, gain=zero(sys) 计算控制系统的零点、增益、绘制系统极点、零点图调用格式 pzmap(sys) pzmap(sys1,sys2,sysn) p z=pzmap(sys)、绘制根轨迹图调用格式 rlocus(sys) rlocus(sys ,k)Rlocus(sys1,sys2,)r k=rlocus(sys) r=rlocus(sys,k)k表示增益 r系统的极点【程序】 G=tf(0.016,0.218,1.436,9.359,0.06,0.268,0.635,6.271) Transfer function:0.016 s3 + 0.218 s2 + 1.436 s + 9.359-0.06 s3 + 0.268 s2 + 0.635 s + 6.271 pzmap(G);title(零极点图) 【程序】 rlocus(G); grid ;title(闭环根轨迹图)8.单位反馈系统前向通道的传递函数为： ，试绘制该系统的Bode图和Nyquist曲线，说明软件绘制曲线与手动绘制曲线的异同。【实验原理】：1)、Bode图调用格式 bode(sys) bode(sys, w) bode(sys1,sys2.sysn)bode(sys1,sys2.sysn,w) bode(sys1,plotstyle1,sys2,plotstyle2,.sysn,plotstylen,)mag, phase, w=bode(sys)w频率区间矢量，控制频率起止范围mag, phase, w输出幅值矢量 输出相角 输出频率区间增益裕度定义：正好使系统频率响应穿过临界点的增益与实际系统增益的比值。相角裕度：为了使系统的频率响应通过临界点必须引入的纯相位滞后量。调用格式 gm pm wcg wcp=margin(sys) gm pm wcg wcp =margin(mag, phase, w)Gm系统增益裕度 pm系统相角裕度wcg wcp-表示交叉频率2）、Nyquist图 调用格式 nyquist(sys) nyquist(sys, w) nyquist(sys1,sys2.sysn) nyquist(sys1,sys2.sysn,w) nyquist(sys1,plotstyle1,sys2,plotstyle2,.sysn,plotstylen,)mag, phase, w=nyquist(sys)mag, phase=nyquist(sys, w)【程序】 G=tf(2,8,12,8,2,1,5,10,10,5,1,0) Transfer function: 2 s4 + 8 s3 + 12 s2 + 8 s + 2-s6 + 5 s5 + 10 s4 + 10 s3 + 5 s2 + s （1）、Bode图【程序】 bode(G);title(Bode图)（2）、Nyquist图【程序】 nyquist(G);title(nyquist图)软件绘制曲线与手动绘制曲线的异同: 软件绘制曲线能准确的显示函数的图像，有利于分析计算。手动绘制曲线是采用对数幅频渐进线的方法得到。这种方法省去了逐点连线的繁琐，有利于分析系统稳定性等方面的问题。9已知某控制系统的开环传递函数，试绘制系统的开环频率特性曲线，并求出系统的幅值与相位裕量。【实验原理】：Nichols图对数幅相特性图(Nichols图)是描述系统频率特性的图示方法。该图纵坐标表示频率特性的对数幅值，以分贝为单位；横坐标表示频率特性的相位角。调用格式 nichols(sys) nichols(sys, w) nichols(sys1,sys2.sysn) nichols(sys1,sys2.sysn,w) nichols(sys1,plotstyle1,sys2,plotstyle2,.sysn,plotstylen,)mag, phase, w=nichols(sys)mag, phase=nichols(sys, w)【程序】 z=;p=0,-1,-2;k=1.5;sys=zpk(z,p,k) Zero/pole/gain: 1.5-s (s+1) (s+2)(1)、开环频率特性曲线【程序】 nichols(sys);title(开环频率特性曲线)（2）、系统的幅值与相位裕量【程序】 Gm,Pm,Wcg,Wcp=margin(sys)Gm = 注释：Gm表示系统的增益裕度 4.0000Pm = 注释：Pm表示系统的相角裕度 41.5114Wcg = 注释：Wcg表示Nyquist曲线与负实轴交点处频率 1.4142 Wcp = 注释：Wcp表示截止频率 0.611010.在SIMULINK中建立系统，该系统阶跃输入时的连接示意图如下。k为学生学号后三位。绘制其单位阶跃响应曲线，分析其峰值时间tp、延迟时间td、上升时间tr、调节时间ts及超调量。 【实验原理】：y,t=step(G) 该函数还同时返回了自动生成的时间变量t，对返回的这一对变量y和t的值进行计算，可以得到时域性能指标。 峰值时间(timetopeak)可由以下命令获得： Y,k=max(y); timetopeak=t(k) 应用取最大值函数max()求出y的峰值及相应的时间，并存于变量Y和k中。然后在变量t中取出峰值时间，并将它赋给变量timetopeak。 最大(百分比)超调量(percentovershoot)可由以下命令得到：C=dcgain(G);Y,k=max(y); percentovershoot=100*(Y-C)/C dcgain( )函数用于求取系统的终值，将终值赋给变量C，然后依据超调量的定义，由Y和C计算出百分比超调量。 上升时间(risetime)可利用MATLAB中控制语句编制M文件来获得。首先简单介绍一下循环语句while的使用。C=dcgain(G);n=1;while y(n)0.98*C)&(y(i) G1=tf(1,1,0);G2=tf(16,1,9);G3=tf(1,1);G4=series(G1,G2);G=feedback(G4,G3,-1)Transfer function: 16-s2 + 9 s + 16 、绘制单位阶跃响应曲线【程序】 step(G);title(单位阶跃响应)、峰值时间 tp【程序】 y,t=step(G); Y,k=max(y);tp=t(k) %计算峰值时间tptp =2.4908 、上升时间tr 由于K=16很小，上升时间无法计算、调节时间ts【程序】 c=dcgain(G); i=length(t); while(y(i)0.98*c)&(y(i) ts=t(i)ts = 1.7888、超调量percentovershoot【程序】 C=dcgain(G);Y,k=max(y); percentovershoot=100*(Y-C)/Cpercentovershoot = -0.3665 超调量为负值说明该系统无超调量因此可近似认为是一阶系统*11. 给定系统如下图所示，试设计一个串联校正装置，使幅值裕度大于h10分贝、相位裕度45（1）、原系统的幅值裕度和相位裕度【程序】 G=tf(100, 0.04, 1, 0); Gw, Pw, Wcg, Wcp=mar