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

目录TOC\o"1-5"\h\z概述 1一、 实验目的 1二、 简述MATLAB语言的特点及其主要功能 1三、 控制系统仿真时常用的方法和指令 21、 控制系统仿真时常用的方法 2a、 数学仿真 •b、 半物理仿真 •c、 全物理仿真 •2、 控制系统仿真时常用的指令 •1） 、Bode图 •、绘制Bode图 •、系统的增益裕度和相角裕度 •2） ,Nichols图 •3） 、Nyquist图 •4） 、一般频率响应图 •5） 、频率响应的奇异值图 •6） 、绘制根轨迹 •四、 实验内容 •五、 心得体会 22六、 参考文献 22MATLAB是一种直观、高效的计算机语言，同时也是一个科学计算平台。它的伴随工具Simulink是用来对真实世界的动力学系统建模、模拟仿真和分析的软件。我们可将综合性和设计性实验项目通过MATLAB在计算机上仿真，使系统的观察实验的动态过程。目前，MATLAB已经成为我们当代大学生必须掌握的基本技能，在设计研究单位和工业部门，MATLAB已经成为研究和解决各种具体工程问题的一种标准软件。在完成了验证性、综合性和设计性实验后，课程设计必不可少。课程设计是工科实践教学的一个重要的环节，目的是培养我们综合运用理论知识分析和解决实际问题的方法和能力,实现由知识向技能的初步化。所以课程设计是培养我们思维创造能力最有效的途径。一、 实验目的1、 培养理论联系实际的科学态度，训练综合运用经典控制理论和相关课程知识的能力。2、 掌握自动控制原理的时域分析法，根轨迹法，频域分析法，以及各种（矫正）装置的作用及用法，能够利用不同的分析方法对给定系统进行性能分析，能根据不同的系统性能指标要求进行合理的系统设计，并调试满足系统的指标。3、 学会使用MATLAB语言及Simulink动态仿真工具进行系统的仿真与调试。4、 锻炼独立思考和动手解决控制系统实际问题的能力。二、 简述MATLAB语言的特点及其主要功能MATLAB是矩阵实验室（MatrixLaboratory）的简称，是美国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',)、系统的增益裕度和相角裕度增益裕度定义：正好使系统频率响应穿过临界点的增益与实际系统增益的比值。相角裕度：为了使系统的频率响应通过临界点必须引入的纯相位滞后量。[调用格式][gmpmwcgwcp]=margin(sys)[gmpmwcgwcp]=margin(mag,phase,w)Gm—系统增益裕度 pm—系统相角裕度wcgwcp---表示交叉频率、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)、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)、一般频率响应图[调用格式]h=freqs(num,den,w)[h,w]=freqs(num,den)[h,w]=freqs(num,den,f)freqs(num,den)Num一传递函数分子矢量 den一传递函数分母矢量W—频率区间矢量、频率响应的奇异值图[调用格式]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—绘制奇异值曲线的特性、绘制根轨迹、求系统的零点、极点、增益[调用格式] p=pole(sys)计算控制系统的极点Z=zero(sys)计算控制系统的零点[z,gain]=zero(sys)计算控制系统的零点、增益、绘制系统极点、零点图[调用格式] pzmap(sys)pzmap(sys1,sys2,„sysn)[pz]=pzmap(sys)、绘制根轨迹图[调用格式] rlocus(sys)rlocus(sys,k)Rlocus(sys1,sys2,„)[rk]=rlocus(sys)r=rlocus(sys,k)k—表示增益 r—系统的极点四、实验内容用matlab语言编制程序，实现以下系统：【实验原理】：线性系统的传递函数模型可一般地表示为：

—、bsm+bsm-ih fbs+bG(s)=— 2 m m+1 n>msn+asn-if fas+a1 n-1 n将系统的分子和分母多项式的系数按降幂的方式以向量的形式输入给两个变量num和den，就可以轻易地将传递函数模型输入到MATLAB环境中。命令格式为：num二[b,b，…,b,b];12mmf1den=[1,a,a，…,a,a];12 n-1n在MATLAB控制系统工具箱中，定义了tf()函数，它可由传递函数分子分母给出的变量构造出单个的传递函数对象。从而使得系统模型的输入和处理更加方便。该函数的调用格式为：G=tf(num,den);将系统增益、零点和极点以向量的形式输入给三个变量KGain、Z和P,就可以将系统的零极点模型输入到MATLAB工作空间中，命令格式为：KGain二K;Z二[—z；—z；…；—z]；1 2 mP二[—P；—P；•••；—P]；1 2 n在MATLAB控制工具箱中，定义了zpk()函数，由它可通过以上三个MATLAB变量构造出零极点对象，用于简单地表述零极点模型。该函数的调用格式为：G=zpk(Z,P,KGain)1) 5s3f24s2f18G(s)=—s4f4s3f6s2f2sf2【程序】>>num=[5,24,0,18];den=[1,4,6,2,2];sys=tf(num,den)Transferfunction:5s3+24s2+18s4+4s3+6s2+2s+24(sf2)(s2f6sf6)2丿G(s)=注释：求解s2+6s+6的根注释：求解s3+3s2+2s+5的根s(s注释：求解s2+6s+6的根注释：求解s3+3s2+2s+5的根-4.7321-1.2679>>B=[1,3,2,5];roots(B)ans=-2.9042-0.0479+1.3112i-0.0479-1.3112i〉〉z二[-2,-4.7321,-4.7321,-1.2679,T.2679];p二[0,-l,-l,T,-2.904-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厂2s(s+1厂3(s+2.904)(sJ+0.0958s+1.722)两环节G1、G2串联，求等效的整体传递函数G(s)【实验原理】：若假定两环节均为单输入单输出的系统SA和SB。两个环节级联：sys=series(SA,SB)27G(s)二 G(s)二-1 s+3 2 s2+2s+1方法一：【程序】〉〉G1=tf(2,[1,3]);G2=tf(7,[1,2,1]);sys=G1*G2Transferfunction:14s3+5s2+7s+3方法二：【程序】〉〉G1=tf(2,[1,3]);G2=tf(7,[1,2,1]);sys=series(G1,G2)Transferfunction:14s3+5s2+7s+3两环节G1、G2并联，求等效的整体传递函数G(s)【实验原理】：若假定两环节均为单输入单输出的系统SA和SB两个环节并联：sys=parallel(SA，SB)27G(s)二 G(s)二1s+3 2 s2+2s+1方法一：【程序】〉〉G1=tf(2,[1,3]);G2=tf(7,[1,2,1]);sys=G1+G2nsferfunction:

2s2+11s+23s3+5s2+7s+3方法二：【程序】>>G1=tf(2,[1,3]);G2=tf(7,[1,2,1]);sys=parallel(G1,G2)Transferfunction:2s2+11s+23两G2(S)二s3+5s2+两G2(S)二已知系统结构如图，求闭环传递函数。其中的环节G1、G2分别为q(s)二s'：；Hi【实验原理】：若假定两环节均为单输入单输出的系统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)Transferfunction:6s八2+215s+5002s八3+9s八2+178s+605正反馈【程序】>>G1=tf([3,100],[1,2,81]);G2=tf(2,[2,5]);sys=feedback(G1,G2,+1)Transferfunction:6s八2+215s+5002s八3+9s八2+166s+205已知某闭环系统的传递函数为g(s)二 100+25 ，求其单位阶跃响应0.16s3+1.96s2+100+25曲线，单位脉冲响应曲线。【实验原理】：LTI模型的阶跃响应函数step()格式：step(sys)功能：绘制系统sys(sys由函数tf、zpk或ss产生)的阶跃响应，结果不返回数据，只返回图形。对多输入多输出模型，将自动求每一输入的阶跃响应。LTI模型的单位冲激响应函数impulse()格式：impulse(sys)功能：绘制系统sys(sys由函数tf、zpk或ss产生)的单位冲激响应，结果不返回数据，只返回图形。、单位阶跃响应【程序】>>G=tf([10,25],[0.16,1.96,10,25])Transferfunction:10s+250.16s八3+1.96s八2+10s+25>>step(G); %计算并绘制系统的单位阶跃响应title('单位阶跃响应')应曲线；分析阻尼比分别为-0.5、-1时系统的稳定性。实验原理】：£=0.5时，利用阶跃响应函数step()绘制o分别取0、2、4、6、8、10时的阶跃响应曲n线。、£=0.5，o分别取0、2、4、6、8、10时该系统的单位阶跃响应曲线n程序】>>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('单位阶跃响应')%计算并绘制系统的单位阶跃响应单位阶跃响应0.2System:G3System:G2Time(sec):1.Time(sec):1.85Amplitude:0.£Amplitude:0.9740System:G4Time(sec):0.456Amplitude:1.16apnt$Systeti:G5Time(sec):1294Amplitude>>step(G1,G2,G3,G4,G5);title('单位阶跃响应')%计算并绘制系统的单位阶跃响应单位阶跃响应0.2System:G3System:G2Time(sec):1.Time(sec):1.85Amplitude:0.£Amplitude:0.9740System:G4Time(sec):0.456Amplitude:1.16apnt$Systeti:G5Time(sec):1294Amplitude1.12System:G1Time(sec):3.69Amplitude:0.9743Time(sec)程序】>>G6=tf([0,100],[1,-10,100]);%>>G7=tf([0,100],[1,-20,100]);%§=-0.5,Wn=10时的单位阶跃响应§=-1,Wn=10时的单位阶跃响应§=-0.5,Wn=10时的单位阶跃响应§=-1,Wn=10时的单位阶跃响应§=-1时的阶跃响应')分析稳定性：从§=-0.5，§=-1时的阶跃响应曲线图可以看出，曲线不收敛，而是发散的，说明系统是不稳定的。设有一高阶系统开环传递函数为G(s)=°.°16s3+。勺8s2+1'436s+9出9，试绘制0.06s3+0.268s2+0.635s+6.271该系统的零极点图和闭环根轨迹图。【实验原理】：、求系统的零点、极点、增益［调用格式］p=pole(sys)计算控制系统的极点Z=zero(sys)计算控制系统的零点［z,gain］=zero(sys)计算控制系统的零点、增益、绘制系统极点、零点图［调用格式］ pzmap(sys)pzmap(sys1,sys2,„sysn)[pz]=pzmap(sys)、绘制根轨迹图[调用格式]rlocus(sys)rlocus(sys,k)Rlocus(sys1,sys2,„)[rk]=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])Transferfunction:0.016s八3+0.218sJ+1.436s+9.3590.06s八3+0.268s八2+0.635s+6.271>>pzmap(G);title('零极点图')零极点图wsA」皿u'oltnE--8-6wsA」皿u'oltnE--8-6【程序】>>rlocus(G);grid;title('闭环根轨迹图’)0：92System:GGan:2.63Poe:-0.265+6.251Damping:0.0423Frequency(rad/sOvershoot(0：92System:GGan:2.63Poe:-0.265+6.251Damping:0.0423Frequency(rad/sOvershoot(%):80.98「.：System:GGain:2.59Pole:-7.7Damping:1Overshoot(%):0Frequency(rad/sec):7.74…艺..3 \………■……■..…;…•…1…“二::::壬®?..••.-•••••-System:GGain:1.83Pole:-0.0646・5.97iDamping:0.0108Overshoot(%):96.7Frequency(rad/sec):5.97单位反馈系统前向通道的传递函数为：G(s)= 2s4+8s3+12s2+8s+2 ,试绘s6+5s5+10s4+10s3+5s2+S制该系统的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]—输出幅值矢量输出相角 输出频率区间增益裕度定义：正好使系统频率响应穿过临界点的增益与实际系统增益的比值。相角裕度：为了使系统的频率响应通过临界点必须引入的纯相位滞后量。[调用格式][gmpmwcgwcp]=margin(sys)

[gmpmwcgwcp]=margin(mag,phase,w)Gm—系统增益裕度pm—系统相角裕度wcgwcp---表示交叉频率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])Transferfunction:2s4+8s3+12s2+8s+2

s6+5s5+10s4+10s3+5s2+s(1)、Bode图程序】>>bode(G);title('Bode图')BodeffioIJ i i iiiiiiI i i iiiiiiI6020n--2LI40-6LISystem:GFrequency(rad/sec):0.0763Magnitude(dB):28.3-HL■10_110106020n--2LI40-6LISystem:GFrequency(rad/sec):0.0763Magnitude(dB):28.3-HL■10_11010Frequency(rad/sec)2)、Nyquist图程序】>>nyquist(G);title('nyquist图')切一遵切一遵A-一河匚一511311二软件绘制曲线与手动绘制曲线的异同:软件绘制曲线能准确的显示函数的图像，有利于分析计算。手动绘制曲线是采用对数幅频渐进线的方法得到。这种方法省去了逐点连线的繁琐，有利于分析系统稳定性等方面的问题。9.已知某控制系统的开环传递函数G(s)= K ,k=1.5，试绘制系统的开环频s(s+1)(s+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)【程序】>>z=[];p=[0,-1,-2];k=[1.5];sys=zpk(z,p,k)Zero/pole/gain:1.5s(s+1)(s+2)(1)、开环频率特性曲线【程序】>>nichols(sys);title('开环频率特性曲线')-225-135-90-180Open-LoopPhase(deg)-120-225-135-90-180Open-LoopPhase(deg)-120丄-270开环频率特性曲线2)、系统的幅值与相位裕量程序】>>[Gm,Pm,Wcg,Wcp]=margin(sys)Gm= 注释：Gm表示系统的增益裕度4.0000Pm= 注释：Pm表示系统的相角裕度41.5114Wcg= 注释：Wcg表示Nyquist曲线与负实轴交点处频率1.4142Wcp= 注释：Wcp表示截止频率0.611010.在SIMULINK中建立系统，该系统阶跃输入时的连接示意图如下。k为学生学号后三位。绘制其单位阶跃响应曲线，分析其峰值时间t、延迟时间t、上升时pd间t、调节时间t及超调量。rs【实验原理】：[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)/Cdcgain()函数用于求取系统的终值，将终值赋给变量C，然后依据超调量的定义，由Y和C计算出百分比超调量。上升时间(risetime)可利用MATLAB中控制语句编制M文件来获得。首先简单介绍一下循环语句while的使用。C=dcgain(G);

n=n+1;endn=n+1;risetime=t(n)在阶跃输入条件下，y的值由零逐渐增大，当以上循环满足y=C时，退出循环，此时对应的时刻，即为上升时间。调节时间(setllingtime)可由while语句编程得到：C=dcgain(G);i=length(t);while(y(i)>0.98*C)&(y(i)<1.02*C)i=i-1;endsetllingtime=t(i)用向量长度函数length()可求得t序列的长度，将其设定为变量i的上限值。(1)K=16时的仿真模型■5o■5o'2102)、示波器显示单位阶跃响应曲线fil^iQ更忆盹詞囲0、时域分析.一.一从仿真示波器中的波形可以看出该二阶系统可近似看成是一阶系统①、闭1);程序】>>G1=tf(1,[1,0]);G2=tf(16,[1,9]);G3=tf(11);>>G4=series(G1,G2);>>G=feedback(G4,G3,-1)Transferfunction:16s2+9s+16②、绘制单位阶跃响应曲线【程序】>>step(G);title('单位阶跃响应')0 0.5 1 1.5 2 2.5Time(sec)单位阶舐响应ooO0 0.5 1 1.5 2 2.5Time(sec)单位阶舐响应ooOOoOoOO<LIP「WQ.E-=l、峰值时间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)<1.02*c)i=i-1;end>>ts=t(i)ts=1.7888、超调量percentovershoot【程序】>>C=dcgain(G);[Y,k]=max(y);>>percentovershoot=100*(Y-C)/Cpercentovershoot=-0.3665超调量为负值说明该系统无超调量因此可近似认为是一阶系统*11.给定系统如下图所示，试设计一个串联校正装置，使幅值裕度大于h〉10分贝、相位裕度Y$45°(1