电子科大-工程控制基础上机(高端版)

基于Matlab的控制系统

的计算机辅助设计、分析及仿真1、Matlab简介2、控制系统的数学描述与建模3、连续系统模型的离散化4、控制系统分析方法（稳定性、时域和频域分析）5、控制系统的设计与仿真6、Simulink仿真1MATLAB简介MATLAB语言的简洁高效性MATLAB语言的科学运算功能MATLAB语言的绘图功能MATLAB庞大的工具箱与模块集MATLAB强大的动态系统仿真功能

MATLAB语言是当前国际上控制领域的首选计算机语言，也是很多理工科专业最适合的计算机数学语言。通过学习可更深入理解和掌握数学问题的求解思想，提高求解数学问题的能力，为今后其他专业课程的学习提供帮助。直接赋值语句赋值变量＝赋值表达式

例：>>a=pi^2a=9.8696

例：表示矩阵>>B=[1+9i,2+8i,3+7j;4+6j5+5i,6+4i;7+3i,8+2j1i]B=1.0000+9.0000i2.0000+8.0000i3.0000+7.0000i4.0000+6.0000i5.0000+5.0000i6.0000+4.0000i7.0000+3.0000i8.0000+2.0000i0+1.0000iMATLAB的基本语句函数调用语句[返回变量列表］＝函数名（输入变量列表）例：［a,b,c]=my_fun(d,e,f,c)冒号表达式

v=s1:s2:s3

该函数生成一个行向量v，其中s1是起始值，s2是步长（若省略步长为1），s3是最大值。

例：用不同的步距生成(0,p)间向量。>>v1=0:0.2:piv1=Columns1through900.20000.40000.60000.80001.00001.20001.40001.6000Columns10through161.80002.00002.20002.40002.60002.80003.0000循环结构for结构while结构MATLAB语言流程控制二维图形绘制二维图形绘制基本语句构造向量:2、控制系统的数学描述与建模

控制系统的数学建模函数命令tf

()，用法：sys=tf(num,den)orsysd=tf(num,den,T)

——建立控制系统的传递函数模型，或者将零极点模型或者状态空间模型转换为传递函数模型；函数命令zpk()，用法：sys=zpk(z,p,k)orsysd=zpk(z,p,k,T)

——建立控制系统的零极点模型;函数命令ss()，用法：sys=ss(A,B,C,D)orsysd=ss(A,B,C,D,T)

——建立控制系统的状态空间模型;函数命令tf2ss(),tf2zp(),zp2ss(),zp2tf(),ss2tf(),ss2zp()——数学模型的转化;基于Simulink的动态结构模型.传递函数

单输入单输出线性连续系统的传递函数为其中m≤n。G(s)的分子多项式的根称为系统的零点,分母多项式的根称为系统的极点。令分母多项式等于零,得系统的特征方程:D(s)=a0sn+a1sn－1+……+an－1s+an=0

因传递函数为多项式之比,所以我们先研究MATLAB是如何处理多项式的。MATLAB中多项式用行向量表示,行向量元素依次为降幂排列的多项式各项的系数,例如多项式P(s)=s3+2s+4,其输入为

>>P=［1024］

注意尽管s2项系数为0,但输入P(s)时不可缺省0。

MATLAB下多项式乘法处理函数调用格式为C=conv(A,B)

例如给定两个多项式A(s)=s+3和B(s)=10s2+20s+3,求C(s)=A(s)B(s),则应先构造多项式A(s)和B(s),然后再调用conv()函数来求C(s)>>A=［1,3］;B=［10,20,3］；>>C=conv(A,B)

C=1050639即得出的C(s)多项式为10s3

+50s2+63s+9

MATLAB提供的conv()函数的调用允许多级嵌套,例如

G(s)=4(s+2)(s+3)(s+4)可由下列的语句来输入

>>G=4*conv(［1,2］,conv(［1,3］,［1,4］))对于其它复杂的表达式,如可由下列语句来输入

>>num=conv(［1,1］,conv(［1,2,6］,［1,2,6］));>>den=conv(［1,0,0］,conv(［1,3］,［1,2,3,4］));>>G=tf(num,den)Transferfunction:

控制系统的典型连接串联环节，用法：

sys=series(sys1,sys2)orsys=sys1*sys2*sys3*…*sysn

并联环节，用法：

sys=parallel(sys1,sys2)orsys=sys1+sys2+sys3+…+sysn

反馈环节，用法：等效传递函数——G=feedback(G1,G2,sign)解编程如下：

num1=[0.1,1];den1=[0.4,0];sys1=tf(num1,den1);

num2=[15];den2=[0.054,1];sys2=tf(num2,den2);

num3=[1.5];den3=[0.12,1];sys3=tf(num3,den3);

sys123=sys1*sys2*sys3

运行结果所以系统等效传递函数为

Transferfunction:2.25s+22.5---------------------------------0.002592s^3+0.0696s^2+0.4s已知某系统前向通道三个模块的传递函数分别为。试求串联连接的等效传递函数解编程如下

num1=[3];den1=[1,1];sys1=tf(num1,den1);

num2=[6,10];den2=[1,2,1];sys2=tf(num2,den2);

sys=sys1+sys2;num=sys.num{1}

den=sys.den{1}

运行结果所以系统的等效传递函数为

num=092213den=1331已知两子系统传递函数分别为：试求两系统并联连接的等效传递函数的num与den向量。3、连续系统模型的离散化

函数命令c2d()，用法:sysd=c2d(sys,T)orsysd=c2d(sys,T,method)

输入参量sys为连续时间模型对象,T,采样周期;

Sysd,带采样时间t的离散事件模型;Method,用来指定离散化采用的方法:

′zoh′——采用零阶保持器法;′foh′——采用一阶保持器法;′tustin′——采用双线形变换法;′prewarp′——采用频率预曲折的双线性变化

法;′matched′——采用零极点匹配法.解编程如下：

k=10;z=[-5];p=[-1,-3,-8];sys=zpk(z,p,k);t=0.1;sys1=tf(sys);

disp('Discretesystem--usingc2dwithzoh')s1=c2d(sys1,t,'zoh')

disp('Discretesystem--usingc2dwithmatched')s2=c2d(sys1,t,'matched')运行结果Discretesystem--usingc2dwithzoh

Transferfunction:0.0401z^2+0.007485z-0.01929---------------------------------z^3-2.095z^2+1.41z-0.3012

Samplingtime:0.1Discretesystem--usingc2dwithmatched

Transferfunction:0.03596z^2+0.01415z-0.02181---------------------------------z^3-2.095z^2+1.41z-0.3012

Samplingtime:0.1已知连续系统零极点增益模型为试采用零阶保持器与零极点匹配法求其离散传递函数.设采样周期T=0.1s.4控制系统的分析方法

控制系统的稳定性分析特征方程的根，即闭环极点应具有负实部或分布在左半s平面上。——直接判断困难。劳斯（Routh）稳定性判据：由特征方程的系数来判断。根轨迹法频率响应特性 根据S平面和Z平面之间的关系，离散系统的稳定性可以由特征方程稳定区域的根在Z平面中的位置来确定——必须位于Z平面中单位园的内部。如果有一个根恰好位于单位园上则系统处于临界稳定，临界稳定在实践中属于不稳定。稳定区域不稳定区域临界稳定Z平面的稳定性条件运行结果：已知系统的闭环传递函数为：

试判断该闭环系统的稳定性。解：根据题意，运行下列MATLAB程序：num=[5410.630.5];den=[100000];[z,p]=tf2zp(num,den)ii=find(real(p)>0);n1=length(ii);if(n1>0),disp('SystemisUnstable');else,disp('SystemisStable');endz=-0.7822+0.5660i-0.7822-0.5660i0.4681+0.6367i0.4681-0.6367i-0.1718p=00000SystemisStable解：根据题意，运行下列MATLAB程序：num=[1];den=[31-1-21];[z,p]=tf2zp(num,den)ii=find(abs(p)>1);n1=length(ii);if(n1>0),disp('SystemisUnstable');else,disp('SystemisStable');endz=

Emptymatrix:0-by-1

p=

-0.7357+0.6859i-0.7357-0.6859i0.5690+0.0753i0.5690-0.0753iSystemisUnstable运行结果：已知系统的特征方程为如下所示,试判断该系统稳定性:控制系统的时域分析方法直接在时域内对系统进行分析；从输出量的时域表达式分析系统的稳定性、动态性能和稳态误差；常用命令函数，step(),impulse(),initial(),dstep,dimpulse,dinitial,nec.任意函数作用下系统的响应

用线性仿真函数lsim来实现，其调用格式为

［y，x］=lsim(num，den，u，t)

式中；y(t)为系统输出响应；x(t)为系统状态响应；u为系统输入信号；t为仿真时间。

求的单位阶跃响应。解：编制Matlab程序如下：num=[25];den=[1,4,25];g=tf(num,den);step(g)或num=[25];den=[1,4,25];step(num,den)求当时的单位冲击响应。解：

编制Matlab程序如下：wn=5;zeta=0.6;num=wn.^2;den=[1,2*zeta*wn,wn.^2];impulse(num,den)设线性离散系统如图所示，且a=1，K=1，T=1s，输入为单位阶跃序列。试分析系统的过渡过程。解：被控对象Z传递函数编制Matlab程序如下：num=[0.3680.264];den=[1-10.632];t=0:1:25;dstep(num,den,t)或[y,x]=dstep(num,den,t)程序运行结果：y=00.36801.00001.39941.39941.14700.89460.80170.86830.99371.07691.08091.03230.98120.96080.97270.99751.01471.01631.00700.99670.99230.99430.99921.00281.0033x=[]调整时间ts=12s（12个采样周期），超调量p=40％，峰值时间tp=3s

，振荡次数N=1.5次，衰减比=5:1

，稳态误差ess=0。例：已知一单位负反馈系统，其开环零极点增益模型为试分析在单位阶跃及斜坡信号作用下系统的稳态误差。利用Matlab进行稳态误差分析单位阶跃信号t=0.1;k=6;z=-0.5;p=[-210];[n1,d1]=zp2tf(z,p,k);s=tf(n1,d1);sys=feedback(s,1);roots(sys.den{1})ans=-0.1084+1.9541i-0.1084-1.9541i-0.7832系统稳定（a）系统稳定性判断sysd=c2d(sys,t)step(sys);t1=[0:t:300]';y=step(sysd,t1);subplot(121),plot(t1,y),grid;subplot(122),ess=1-y;plot(t1,ess),gridess(length(ess))阶跃响应及稳态误差ans=-8.5931e-014系统为一型系统，其Kp=，即ess=0。t=0.1;k=6;z=-0.5;p=[-210];[n1,d1]=zp2tf(z,p,k);s=tf(n1,d1);sys=feedback(s,1);num=sys.num{1};den=[sys.den{1},0];sys=tf(num,den);sysd=c2d(sys,t)step(sys);t1=[0:t:50]';y=step(sysd,t1);subplot(121),plot(t1,[t1,y]),grid;subplot(122),ess=t1-y;plot(t1,ess),gridess(length(ess))斜坡响应及稳态误差系统为一型系统，其Kv=C，即ess=C。ans=-0.6678控制系统的频域分析方法以系统的频率特性来研究系统的性能；对数频率特性图（Bode图），

[mag,phase,w]=dbode(sys,w)极坐标图（Nyquist图）

[re,im,w]=dnyquist(sys,w)求系统增益裕度和相位裕度的函数，margin();[Gm,Pm,Wcp,wcg]=margin(sys)根轨迹图，[p,z]=pzmap(sys),[p,z]=rlocus(sys).对数频率特性Matlab提供的相关函数为：函数dbode()绘制线性离散系统对数频率特性图。函数dbode()的使用格式如下：[mag,phase,w]=dbode(num,den,Ts,w)当包含左方变量时，屏幕上不显示伯德图，返回角频率点w处的幅值和相位角（单位为度）。函数[Gm,Pm,ωcg,ωcp]=margin(sys)除绘制线性离散系统对数频率特性图外，还可以求出系统的幅值裕度Gm、相位裕度Pm以及所在点对应的频率值ωcg,ωcp试求所示离散系统，在T=1s，K=1时的开环对数频率特性（伯德图）。在T=1s，K=1时，编制Matlab程序如下：num=[0.368,0.264];den=[1,-1.368,0.386];sys=tf(num,den,1)sysw=d2c(sys,'tustin')syszpk=zpk(sysw)margin(syszpk)title('离散系统伯德图');grid解：可知，系统的开环传递函数为直接用dbode()

函数：num=[0.368,0.264];den=[1,-1.368,0.368];dbode(num,den,1)title('离散系统伯德图');grid绘制程序2系统的相位裕度幅值裕度系统性能预测：稳定，有一定超调。结论：极坐标频率特性Matlab提供的相关函数为：

Dnyquist(num,den,t)试求上例所示离散系统，在T=1s，K=1时的极坐标频率特性曲线。解：编制程序如下：num=[0.368,0.264];den=[1,-1.368,0.386];dnyquist(num,den,1)title('离散系统极坐标图');已知系统传递函数为绘制零极点图和根轨迹图.MATLAB命令为:num=[0,0,2,8,12,8,2];den=[1,5,10,10,5,1,0];sys=tf(num,den);subplot(1,2,1)pzmap(sys)[p,z]=pzmap(sys)subplot(1,2,2)rlocus(sys)运行结果：p=0-1.0010+0.0007i-1.0010-0.0007i-0.9996+0.0012i-0.9996-0.0012i-0.9988

z=-1.0002-1.0000+0.0002i-1.0000-0.0002i-0.9998

解：系统的开环Ｚ传递函数试求所示离散系统，在T=1s和2s时的根轨迹图，并求出相应的临界增益。例：编制Matlab程序如下：t1=1;t2=2;k=1;n=2;fori=1:nifi<=1t=t1;elseifi<=2t=t2;endb1=t-1+exp(-t);b2=1-exp(-t)-t*exp(-t);a0=1;a1=1+exp(-t);a2=exp(-t);num=[k*b1,k*b2];den=[a0,-a1,a2];figure(i)

rlocus(num,den)title('根轨迹图')holdp=0:0.07:2*pi;x=sin(p);y=cos(p);plot(x,y,'.')holdv=[-4,2,-2,2];axis(v)grid[k,poles]=rlocfind(num,den)endT=1秒T=2秒零界增益K=2.52零界增益K=0.58结论采样周期T与系统的临界增益K密切相关。T增长，K降低。系统的零极点改变时，根轨迹的形状也发生变化，系统的性能也跟着变化。为了满足性能指标的要求，可以采用零极点对消的办法，或加入新的零极点，来改善系统的性能。5控制系统的设计和仿真频域设计方法，根据系统的传递函数进行设计。时域设计方法，基于系统的状态空间模型进行设计。频域设计方法超前校正的设计步骤如下：

(1)根据稳态性能要求，确定系统开环增益K。

(2)根据外环增益K，画出校正前系统的Bode图，并计算未校正系统的相位裕度等性能指标．以检测是否满足要求.(3)确定需要增加的最大相位超前角：式中,γ0为期望相位裕度,γ1为原来开环系统相位裕度.(4)由确定α值及最大相位超前角所对应的频率ωm,并取新的幅值穿越频率ωcg=

ωm。

(5)确定超前校正器传递函数，并绘制校正后的Bode图和进行验证。经典控制理论设计已知单位正反馈系统被控对象的传递函数为试用Bode图设计方法对系统进行超前校正设计。性能指标要求如下：(1)在斜坡信号r(t)=v0t,ess≤0.0001v0;(2)在系统校正后,相位裕量48º<γ<53º.解:(1)求K0速度误差系数Kv=K=K0,(K是系统的开环增益)

Ess=v0/Kv=v0/K0≤0.0001v0所以,Kv=K0≥1000,取K0=1000(2)作