MATLAB与控制系统仿真实验

自动化学院（院、系）

姓名

工业自动化黄国盛

专业

144班

组

课

教师评定

学号201421714406

实验日期

2017年4月6号

实验一MATLAB绘图基础

一、实验目的

1.2.3.了解MATLAB常用命令和常见的内建函数使用。熟悉矩阵基本运算以及点运算。掌握MATLAB绘图的基本操作：向量初始化、向量基本运算、绘图命令plot,plot3,mesh,surf使用、绘制多个图形的方法。二、实验设备及软件

PC机一台、MATLABR2015b三、实验内容

建立并执行M文件multi_plot.m，使之画出如图的曲线。四、实现程序如下

x=linspace(0,4*pi);%向量初始化y1=sin(2*x);%定义y1函数y2=2*cos(2*x);%定义y2函数plot(x,y1,x,y2);%在一个plot坐标图中绘制两条曲线axis([04*pi-22]);%定义坐标范围title('Plotofy=sin(2x)anditsderivative');%插入标题xlabel('x');%定义x轴名称ylabel('y');%定义y轴名称legend('sin(2x)','2cos(2x)',-1);%曲线标签gridon;%添加网格第1页共27页五、实现思考

1、实现同时画出多图还有其它方法，请思考怎样实现，并给出一种实现方法。(参考程序如下)%%subplotx=linspace(0,4*pi);%向量初始化y1=sin(2*x);%定义y1函数subplot(2,1,1)%两行一列个图，第一个plot(x,y1);%绘制y1曲线axis([04*pi-22]);%定义坐标范围title('Plotofy=sin(2x)');%插入标题xlabel('x');%定义x轴名称ylabel('y');%定义y轴名称gridon;%添加网格subplot(2,1,2)%两行一列个图，第二个y2=2*cos(2*x);%定义y2函数plot(x,y2);%绘制y2曲线axis([04*pi-22]);%定义坐标范围title('Plotderivativeofy=sin(2x):y=2cos(2x)');%插入标题xlabel('x');%定义x轴名称ylabel('y');%定义y轴名称gridon;%添加网格%holdon;holdoffx=linspace(0,4*pi);%向量初始化y1=sin(2*x);%定义y1函数plot(x,y1,'r--');%绘制y1曲线holdon;%保持图形y2=2*cos(2*x);%定义y2函数plot(x,y2);%绘制y2曲线axis([04*pi-22]);%定义坐标范围title('Plotofy=sin(2x)anditsderivative');%插入标题xlabel('x');%定义x轴名称第2页共27页ylabel('y');%定义y轴名称legend('sin(2x)','2cos(2x)',-1);%曲线标签gridon;%添加网格holdoff;%取消图形保持%%figure1,figure2...x=linspace(0,4*pi);%向量初始化%%figure1,figure2...x=linspace(0,4*pi);%向量初始化figure(1);%定义第一个figure窗体y1=sin(2*x);%定义y1函数plot(x,y1,'r--');%绘制y1曲线axis([04*pi-11]);%定义坐标范围title('Plotofy=sin(2x)');%插入标题xlabel('x');%定义x轴名称ylabel('y');%定义y轴名称gridon;%添加网格figure(2);%定义第二个figure窗体y2=2*cos(2*x);%定义y2函数plot(x,y2);%绘制y2曲线axis([04*pi-22]);%定义坐标范围title('Plotderivativeofy=sin(2x):y=2cos(2x)');%插入标题xlabel('x');%定义x轴名称ylabel('y');%定义y轴名称gridon;%添加网格第3页共27页2、思考三维曲线（plot3）与曲面(mesh,surf)的用法，（1）绘制参数方程x(t)t3sin(3t)et,yt3cos(3t)et,zt2的三维曲线；（2）绘制二元函数zf(x,y)(x22x)ex2y2xy，在XOY平面内选择一个区域(-3:0.1:3,-2:0.1:2)，然后绘制出其三维表面图形。(以下给出PLOT3和SURF的示例)3、%eg.plot3：绘制参数方程x=sin(t),y=cos(t),z=8*t的三维曲线t=0:pi/50:10*pi;%向量初始化plot3(sin(t),cos(t),8*t);%绘制参数方程x=sin(t),y=cos(t),z=8*t的三维曲线%绘制参数方程x(t)tsin(3t)e,ytcos(3t)e,zt的三维曲线t=0:pi/50:10*pi;%向量初始化plot3((t.^3).*sin(3*t).*exp(-t),(t.^3).*cos(3*t).*exp(-t),t.^2);%绘制参数方程的三维曲线3t3t2第4页共27页%eg.mesh:绘制二元方程z=sin(x2y2)

xy

22的三维曲面，x,y的取值范围为[-8，8]。[x,y]=meshgrid(-8:0.5:8,-8:0.5:8);c=sqrt(x.^2+y.^2)+eps;z=sin(c)./c;surf(x,y,z);shadingflat;%定义x,y轴范围和网格%定义方程分子%定义z方程%绘制曲面图%添加彩色效果，在faceted的基础上去掉图上的网格线2xzf(x,y)(x2x)e%绘制绘制二元函数2y2xy[x,y]=meshgrid(-3:0.1:3,-2:0.1:2);c=x.^2-2*x;z=(x.^2-2*x).*exp(-(x.^2)-(y.^2)-(x.*y));surf(x,y,z);shadingflat;%定义x,y轴范围和网格%定义方程分子%定义z方程%绘制曲面图%添加彩色效果，在faceted的基础上去掉图上的网格线第5页共27页第6页共27页仲恺农业工程学院实验报告纸

自动化学院（院、系）

姓名

工业自动化黄国盛

专业

144班

组

教师评定

课

学号201421714406

实验日期

2017年4月13号

实验二基于Simulink的控制系统仿真

一、实验目的

1．掌握MATLAB软件的Simulink平台的基本操作；2．能够利用Simulink平台研究PID控制器对系统的影响；三、实验设备及软件

PC机一台、MATLABR2015b三、实验原理

PID（比例-积分-微分）控制器是目前在实际工程中应用最为广泛的一种控制策略。PID算法简单实用，不要求受控对象的精确数学模型。1．模拟PID控制器典型的PID控制结构如图1所示。PID控制器比例`r(t)e(t)积分微分u(t)y(t)对象模型图1连续系统PID控制器的表达式为典型PID控制结构tx(t)Kpe(t)KIe()dKD0de(t)dt（1）式中，KPKI

,和KD

分别为比例系数，积分系数和微分系数，分别是这些运算的加权系数。对式（7-21）进行拉普拉斯变换，整理后得到连续PID控制器的传递函数为GC(s)KP

KPKI

,KI1KDsKP(1TDs)sTIsTIKP/KI,TDKD/KP

第7页共27页（2）），PID控制器的性能也就确定下显然和KD

这3个参数一旦确定（注意来。为了避免微分运算，通常采用近似的PID控制器，气传递函数为GC(s)KP(1

TDs1)TIs0.1TDs1（3）四、实验过程

PID控制器的KP,KI和KD这3三个参数的大小决定了PID控制器的比例，积分和微分控制作用的强弱。下面请通过一个直流电动机调速系统，利用MATLAB软件中的Simulink平台，使用期望特性法来确定这3个参数的过程。并且分析这3个参数分别是如何影响控制系统性能的。【问题】某直流电动机速度控制系统如图2所示，采用PID控制方案，使用期望特性法来确定KP,KI和KD这3三个参数。期望系统对应的闭环特征根为：-300，-300，-30+j30和-30-j30。请建立该系统的Simulink模型，观察其单位阶跃响应曲线，并且分析这3个参数分别对控制性能的影响。转速给定GcPIDG126350s600晶闸管整流器G2+—G313.5s传动装置G45.3转速+—转速控制器60s60直流电动机电动势系数0.012测速反馈系数图2直流电动机PID控制系统（1）使用期望特性法来设计PID控制器。首先，假设PID控制器的传递函数为:所示的系统闭环的传递函数为GC(s)KP

KIKDs

KKKs,其中P,I和D这3个参数待定。图2113120550(KDs2KPsKI)

GB(s)4s660s3(36801357447KD)s2(4860001357447KP)s1357447KI

如果希望闭环极点为：-300，-300，-30+j30和-30-j30，则期望特征多项式为：s4660s3127800s26480000s162106。对应系数相等，可求得：KD0.067，KP4.4156，KI119.34。在命令窗口中输入这3个参数值，并且建立该系统的Simulink模型，如图3所示。第8页共27页图3直流电动机PID控制系统的Simulink仿真模型输入信号为单位阶跃信号，在t=0s时从0变化到1。系统响应曲线如图4所示。图4直流电动机PID控制系统响应曲线（2）分析比例系数在KP对控制性能的影响KI119.34和KD0.067保持不变的情况下，KP分别取值0.5，5和20，系统的响应曲线如图5所示。KP取值较小时系统的响应较慢，而当KP取值较大时系统的响应速度较快，但超调量增加。可见，当图5改变KP时的系统响应曲线第9页共27页（3）分析积分系数在KI对控制性能的影响KD0.067和KP4.4156保持不变的情况下，KI分别取值20，120，300，系统的响应曲线如图6KI取值较小时系统响应进入稳态的速度较慢。而当KI取较大值时系统的响应进入稳态的速所示。可见，当度较快，但超调量增加。图6改变（4）分析微分系数在KI时的系统响应曲线KD对控制性能的影响KP4.4156和KI119.34保持不变的情况下，KD分别取值0.01，0.07，0.2，系统的响应曲线KD取值较小时系统响应对变化趋势的调节较慢。超调量较大。而当KD取值较大时系KD取值过大时，对变化趋势的调节过强，阶跃响应的如图7所示。可见，当统的响应进入稳态的速度较快。但是超调量增加。当初期出现尖脉冲。图7改变KD时的系统响应曲线第10页共27页仲恺农业工程学院实验报告纸

自动化学院（院、系）

姓名

工业自动化黄国盛

专业

144班

组

课

教师评定

学号201421714406

实验日期

2017年4月20号

实验三控制系统数学模型转换及MATLAB实现

一、实验目的

1.熟悉MATLAB的实验环境。2.掌握MATLAB建立系统数学模型的方法。二、实验设备及软件

PC机一台、MATLABR2015b三、实验内容

（注：实验报告只提交第2题）2、用MATLAB做如下练习。（1）用2种方法建立系统G(s)（2）用2种方法建立系统G(s)

s2

的多项式模型。s25s1010(s1)

的零极点模型。(s1)(s5)(s10)dc(t)

10c(t)20e(t)dtdb(t)205b(t)10c(t)

dt（3）如图，已知G（s）和H（s）两方框对应的微分方程是：6

且初始条件为零。试求传递函数C(s)/R(s)及E(s)/R(s)。R10MEG（s）CBH（s）四、实现代码和结果

%2-1:getthetransferfunction%method1%SYS=TF(NUM,DEN)createsacontinuous-timetransferfunctionSYSwith第11页共27页%numerator(s)NUManddenominator(s)DEN.TheoutputSYSisaTFobject.num=[12];%定义分子多项式den=[1510];%定义分母多项式tfsys=tf(num,den)%生成传递函数%method2%S=TF('s')specifiesthetransferfunctionH(s)=s(Laplacevariable).%Z=TF('z',TS)specifiesH(z)=zwithsampletimeTS.%Youcanthenspecifytransferfunctionsdirectlyasrationalexpressions%inSorZ,e.g.,%s=tf('s');H=(s+1)/(s^2+3*s+1)s=tf('s');%s=tf('s')用拉普拉斯算符s的有理表达式来创建传递函数tfsys=(s+2)/(s^2+5*s+10)%生成传递函数结果如下：tfsys=s+2--------------s^2+5s+10Continuous-timetransferfunction.%2-2:Createzero-pole-gainmodels%method1%SYS=ZPK(Z,P,K)createsacontinuous-timezero-pole-gain(ZPK)%modelSYSwithzerosZ,polesP,andgainsK.TheoutputSYSis%aZPKobject.k=10;%定义比例系数z=[-1];%定义零点p=[-2-5-10];%定义极点zpksys=zpk(z,p,k)%零极点模型%method2%S=ZPK('s')specifiesH(s)=s(Laplacevariable).%Z=ZPK('z',TS)specifiesH(z)=zwithsampletimeTS.%YoucanthenspecifyZPKmodelsdirectlyasrationalexpressionsin%SorZ,e.g.,%z=zpk('z',0.1);H=(z+.1)*(z+.2)/(z^2+.6*z+.09)S=zpk('s');%用零极点函数定义zpksys=10*(S+1)/((S+2)*(S+5)*(S+10))%生成零极点模型结果如下：zpksys=10(s+1)------------------(s+2)(s+5)(s+10)Continuous-timezero/pole/gainmodel.第12页共27页%2-3S=tf('s');G=20/(6*s+10);%GETG(S)H=10/(20*s+5);%GETH(S)Gf=feedback(G,H,-1);%合并为闭环负反馈Gcr=series(10,Gf)%GETC/RGer=10/(1+G*H)%GETE/R结果如下：Gcr=4000s+1000---------------------120s^2+230s+250Ger=1200s^2+2300s+500-----------------------120s^2+230s+250Continuous-timetransferfunction.第13页共27页仲恺农业工程学院实验报告纸

自动化学院

（院、系）

姓名

工业自动化黄国盛

专业

144

班

组

教师评定

课

学号201421714406

实验日期2017年4月27日

实验四控制系统的MATLAB辅助分析

一、实验目的

1.2.3.4.熟悉基于MATLAB的系统稳定性分析。熟悉基于MATLAB的时域分析。熟悉基于MATLAB的根轨迹分析。熟悉基于MATLAB的频域分析。二、实验设备及软件

PC机一台、MATLABR2015b三、实验内容

（注：实验报告只提交第2题）2、基于MATLAB完成如下练习。（1）结合理论分析开环增益K0和时间常数T改变对系统稳定性及稳态误差的影响。系统传递函数为G(s)

10K0。（保持T＝0.1，改变K0（取1～5）分别绘制相应阶跃响应曲线，并分析K0变s(0.1s1)(Ts1)

化对系统稳定性的影响；在K0＝1（系统稳定）和K0＝2（系统临界稳定）2种情况下，分别绘制T＝0.1和T＝0.01时系统的阶跃响应，分析T变化对系统阶跃响应及稳定性的影响）（2）系统开环传递函数为G(s)

K

，要求：绘制并记录根轨迹；确定根轨迹的分离点与s(0.2s1)(0.5s1)相应的根轨迹增益；确定临界稳定时的根轨迹增益。s2（3）设系统如图所示。其中Gc(s)为改善系统性能而加入的校正装置。若Gc(s)可从s,s,三种传递s20

2函数中任选一种，分别做出相应的根轨迹图，分析说明应选哪一种校正装置。R－G1（s）－G2（s）CGc（s）第14页共27页注：图中G1（s）＝100/(s+2),G2（s）=10/(s(s+10))（4）已知系统开环传函G(s)单分析K值变化对系统的影响。（5）已知单位反馈系统G(s)稳定裕度。K0，画出当K0取10,500,1000时的Bode图和Nyquist图，并简2s10s50K

，设K分别为4和10，试确定系统是否稳定，如稳定时求出系统的(s1)3四、实验结果和分析

一、结合理论分析开环增益K0和时间常数T改变对系统稳定性及稳态误差的影响。

G(s)

设系统的开环传递函数为：10K0s(0.1s1)(Ts1)

K1，即令R1R2100K，建立系统数学模型，绘制其取T0.1，即令R100K,C1F；取0阶跃曲线。（1）首先理论上分析K0对稳定性的影响。保持T0.1不变，改变K0，令K0分别等于1，2，3，4，5，K0值范围，并对上述情况进行稳定性判断。用劳斯判据求出使系统稳定的键入程序：%定义元件参数R1=10^5;R=10^5;R2=[1,2,3,4,5]*10^5;C1=10^(-6);C2=10^(-7);T=[R*C1,R*C2];%建立系统传递函数，并绘制其阶跃响应曲线fori=1:5K0(i)=R2(i)/R1;第15页共27页num=10*K0(i);den=[0.1*T(1),0.1+T(1),1,0];Gopen=tf(num,den)Gclose=feedback(Gopen,1,-1)figure(i)step(Gclose)endK01时系统的阶跃响应曲线：

K02时系统的阶跃响应曲线：

第16页共27页K03时系统的阶跃响应曲线：

K04时系统的阶跃响应曲线：

K05时系统的阶跃响应曲线：运行结果可知，K02时系统临界稳定；随着K0的增加，系统将趋于不稳定。第17页共27页（2）在K01（系统稳定）和K02（系统临界稳定）两种情况下，分别绘制T0.1和T0.01（即保持R100K不变，C分别取1F和0.1F时）系统的阶跃响应，分析T值变化对系统阶跃响应及稳定性的影响。键入程序：%定义元件参数R1=10^5;R=10^5;R2=[1,2,3,4,5]*10^5;C1=10^(-6);C2=10^(-7);T=[R*C1,R*C2];%取K0=1，分别绘制T=0.1和T=0.01时的阶跃响应曲线K0=R2(1)/R1;fori=1:2num=10*K0;den=[0.1*T(i),0.1+T(i),1,0];Gopen(i)=tf(num,den)Gclose(i)=feedback(Gopen(i),1,-1)endfigure(1)step(Gclose(1),'r',Gclose(2),'g')由图可知，时间常数T减小时，系统动态性能得到改善。同理绘制K0=2，分别绘制T=0.1和T=0.01时的阶跃响应曲线。第18页共27页（3）取K0=1和T=0.01，改变系统的输入信号（分别取单位阶跃、单位斜坡、单位加速度），观察在不同输入下的响应曲线及相应的稳态误差。键入程序：K0=1;T=0.01;Num=10*K0;den=[0.1*T,0.1+T,1,0];Gopen=tf(num,den)Gclose=feedback(Gopen,1,-1)figure(1)step(Gclose)figure(2)t=0:0.01:5;u1=t;lsim(Gclose,u1,t)figure(3)t=0:0.01:5;l=length(t);fori=1:lu2(i)=t(i)^2/2;endlsim(Gclose,u2,t)第19页共27页系统单位阶跃响应曲线：由图可见，系统对于单位阶跃响应输入可以实现无差跟踪。系统单位斜坡响应曲线：

由图可见，系统对于单位斜坡输入可以跟踪，但存在一定稳态误差。系统单位加速度响应曲线：

由图可见，系统对于单位加速度输入随时间的推移，误差越来越大，即不能跟踪。第20页共27页（4）改变K0值，绘制系统在单位斜坡输入下的响应曲线，分析改变开环放大系数对系统稳态误差的影响。%定义元件参数K0=[12];T=0.01;t=0:0.01:3;u1=t;fori=1:2num(i)=10*K0(i);den=[0.1*T,0.1+T,1,0];Gopen(i)=tf(num(i),den);Gclose(i)=feedback(Gopen(i),1,-1);endfigure(1)lsim(Gclose(1),'r',Gclose(2),'g',u1,t)由图可见，开环增益K0越大，系统的稳态误差就越小，故可以通过增大开环增益K0来减小稳态误差。（5）改变系统型次，绘制系统在单位斜坡输入下的响应曲线，分析改变系统型次对系统稳态误差的影响。K0=[12];num=[10*K0(1)10*K0(1)];den1=[0.1*T,0.1+T,1];den2=[0.1*T,0.1+T,1,0];den3=[0.1*T,0.1+T,1,0,0];den=[den1,den2,den3];Gopen1=tf(num,den1)Gopen2=tf(num,den2)Gopen3=tf(num,den3)Gopen=[Gopen1,Gopen2,Gopen3]fori=1:3Gclose(i)=feedback(Gopen(i),1,-1);endfigure(2)holdonlsim(Gclose(1),'r',Gclose(2),'g',Gclose(3),'b',u1,t)第21页共27页由图可见，系统型次越高，对斜坡输入的稳态误差越小，故或通过提高系统型次来减小或消除稳态误差。内容小结：1.K0影响系统的稳定性。K0减小，稳定性增加。2.K0一定时，T决定系统的动态性能。T减小，系统调节时间减小。3.K0=1，T=0.01时，在各种输入下系统的响应。对于1型系统，在单位阶跃输入下无稳态误差；在单位斜坡输入下，有恒定的稳态误差；在单位加速度输入下，系统误差为无限大。4．消除稳态误差的方法。（1）K0增加，稳态误差减小；（2）系统型别增加，稳态误差减小。二、系统开环传递函数为G(s)

K

，要求：绘制并记录根轨迹；确定根轨迹的

s(0.2s1)(0.5s1)分离点与相应的根轨迹增益；确定临界稳定时的根轨迹增益。

实现代码和结果如下：>>clearall>>num=[1];den=conv(conv([10],[0.21]),[0.51]);sys=tf(num,den);figure(1);rlocus(sys);figure(2);第22页共27页margin(sys);[gm]=margin(sys)得：gm=7.0000s2三、设系统如图所示。其中Gc(s)为改善系统性能而加入的校正装置。若Gc(s)可从s,s,

s202三种传递函数中任选一种，分别做出相应的根轨迹图，分析说明应选哪一种校正装置。

R－G1（s）－G2（s）CGc（s）注：图中G1（s）＝100/(s+2),G2（s）=10/(s(s+10))实现代码和结果：G1=zpk([0-20],[-23.25-3.375-5.625i-3.375+5.625i],1);G2=zpk([00-20],[-23.25-3.375-5.625i-3.375+5.625i],1);第23页共27页G3=zpk([00],[-23.25-3.375-5.625i-3.375+5.625i],1);z=0.707;figure(1);rlocus(G1);sgrid(z,'new');K=3.02;Kt=K/10;holdon;rlocus(G1,K);figure(2);rlocus(G2);figure(3);rlocus(G3);num1=[100];den1=[120];num2=[10];den2=[1100];num3=[Kt0];den3=[001];[numf,denf]=feedback(num2,den2,num3,den3);[numc,denc]=series(num1,den1,numf,denf);[num,den]=cloop(numc,denc);sys=tf(num,den);t=0:0.001:5;figure(4);step(sys,t);gridon;第24页共27页结论：可见，加入校正装置后改变了系统闭环极点的分布。由根轨迹图可知，加第一种装置后，适当选择K，可使系统闭环阻尼为0.707，达到最佳工程阻尼比。K0，画出当K0取10,500,1000时的Bode图和Nyquist

s210s50

图，并简单分析K值变化对系统的影响。

四、已知系统开环传函G(s)

实现代码和结果如下：k=[105001000];第25页共27页fori=1:3G(i)=tf(k(i),[11050]);endfigure(1);bode(G(1),'r:',G(2),'b--',G(3))title('系统K/(s^2+10s+50)Bode图,K=10,500,1000','fontsize',16);gridonfigure(2);nyquist(G(1),'r',G(2),'b',G(3),'g')结论：随着K的增加，系统的幅频特性向上平移，相频特性不变，即K值只影响系统幅频特性的起点，不改变其形状，且对相频特性无影响。K值的增大不会影响nyquist曲线的形状，只会改变G(jw)向量幅值的大小，从而改变曲线的包围区域。五、已知单位反馈系统G(s)求出系统的稳定裕度。

实现代码和结果如下：K

，设K分别为4和10，试确定系统是否稳定，如稳定时3(s1)

第26页共27页%h幅值裕度%r相角裕度%wx与-180度线相交频率%wc剪切频率K1=[4]K2=[10]b1=conv([11],conv([11],[11]))[h1,r1,wx1,wc1]=margin(K1,b1)[h2,r2,wx2,wc2]=margin(K2,b1)nyquist(K1,b1);holdon;nyquist(K2,b1);得：K1=4h1=2.0003r1=27.1424wx1=1.7322wc1=1.2328K2=10h2=0.8001r2=-7.0310wx2=1.7322wc2=1.9083由Nyquist判据的K=4时，闭环系统稳定，此时h>1,γ>0;K=10时，闭环系统不稳定，此时h<1,γ<0;对于最小相位系统，当h>1（h>0dB）或γ>0时，闭环系统稳定，h、γ越大，系统相对稳定性越好。第27页共27页