3 MATLAB控制系统仿真讲的很多.docx

3 MATLAB控制系统仿真讲的很多.docx

《3 MATLAB控制系统仿真讲的很多.docx》由会员分享，可在线阅读，更多相关《3 MATLAB控制系统仿真讲的很多.docx（42页珍藏版）》请在冰豆网上搜索。

3MATLAB控制系统仿真讲的很多

2.5应用MATLAB控制系统仿真

MATLAB是一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算和图形显示于一体，构成了一个方便的界面友好的用户环境。

由控制领域专家推出的MATLAB工具箱之一的控制系统（ControlSystem），在控制系统计算机辅助分析与设计方面获得了广泛的应用，并且MATLAB工具箱的内容还在不断增加，应用范围也越来越宽。

控制系统的分析与设计方法，不论是古典的还是现代的，都是以数学模型为基础进行的。

MATLAB可以用于以传递函数形式描述的控制系统。

在本节中，我们首先以一个典型的动力学系统“弹簧—重物—阻尼器”的数学模型为例，说明如何使用MATLAB进行辅助分析。

之后，我们讨论传递函数和结构图。

特别的，主要介绍以下内容：

如何使用MATLAB求解多项式，计算传递函数的零点和极点，计算闭环传递函数，计算结构图的等效变换以及闭环系统对单位阶跃输入的响应等。

这部分所包含的MATLAB函数有：

roots，roots1，series，parallel，feedback，cloop，poly，conv，polyval，printsys，pzmap和step等。

2.5.1.举例

一个弹簧—重物—阻尼器动力学系统如图2-1所示。

重物M的位移由y（t）表示，用微分方程描述如下：

该系统在初始位移作用下的瞬态响应为

其中θ=cos-1ζ，初始位移是y（0）。

系统的瞬态响应当ζ＜1时为欠阻尼，当ζ＞１时为过阻尼，当ζ＝1时为临界阻尼。

考虑

过阻尼情况：

y（0）=0.15m ωn=

（弧度/秒）

 （

）

欠阻尼情况：

y（0）=0.15m  ωn=

（弧度/秒）

 （

）

利用MATLAB程序—unforced.m，可以显示初始位移为y（0）的物体自由运动曲线，如图2-63所示。

在unforced.m程序中，变量y（0），ωn，t，ζ1和ζ2的值由指令直接输入工作区，然后运行unforced.m程序就可以产生响应曲线。

MATLAB提供了一种交互式的应用环境，我们可以在指令窗口随时修改ωn，ζ1和ζ2值并再次运行unforced.m程序。

在上述欠阻尼和过阻尼情况下的响应曲线如图2-64所示。

可以看到程序运行后阻尼系数的值被标注在不同的曲线上，以防止混淆。

对于弹簧—重物—阻尼器系统，利用MATLAB求微分方程的解是容易而有效的。

一般说来，当分析一个具有多样性的输入和初始条件以及元件参数的闭环反馈控制系统时，欲直接得到这些因素与系统响应的关系是比较困难的。

在这种情况下，可以利用MATLAB非常有效的数值计算能力，反复运行MATLAB程序并绘出曲线图，就可以得出结论。

 >>y0=0.15;wn=sqrt

（2）;

  >>zeta1=3/（2*sqrt

（2）;zeta2=1/（2*sqrt

（2）;

  >>t=[0:

0.1:

10];

  >>unforced

（a）MATLAB指令窗口

  %unforced.m

  %计算系统在给定初始条件下的自由运动

  t1=acos（zeta1）*ones（1,length（t）;

  t2=acos（zeta2）*ones（1,length（t）;

  c1=（y0/sqrt（1-zeta1^2）;c2=（y0/sqrt（1-zeta2^2）;

  y1=c1*exp（-zeta1*wn*t）sin（wn*sqrt（1-zeta1^2）*t+t1）;

  y2=c2*exp（-zeta2*wn*t）sin（wn*sqrt（1-zeta2^2）*t+t2）;

  %计算运动曲线的包络线

  bu=c2*exp（-zeta2*wn*t）;bl=-bu;

  %画图

  plot（t,y1,‘-’,t,y2,‘-’,t,bu,‘--’,bl,‘--’）,grid

  xlabel（‘Time[sec]’）,ylabel（‘y（t）Displacement[m]’）

  text（0.2,0.85,[‘oeverdampedzeta1=’,num2str（zeta1）,]）

  text（0.2,0.80,[‘underdampedzeta2=’,num2str（zeta2）,]）

  （b）分析弹簧—重物—阻尼器的MATLAB程序unforced.m

图2-64弹簧—重物—阻尼器的自由运动曲线

MATLAB可以用来分析以传递函数形式描述的系统。

由于传递函数是多项式的比值，所以如何处理多项式是MATLAB首先需要解决的问题。

记住，分子多项式和分母多项式都必须在MATLAB指令中指定。

在MATLAB中多项式由行向量组成，而这些行向量包含了降次排列的多项式系数。

例如多项式p（s）=1s3+3s2+0s1+4s0，按图2-65的格式输入p=[1304]，请注意，尽管s的系数为0，它也一定要包含在确定p（s）的行向量中。

如果p是一个包含降幂排列的p（s）系数的行向量，那么roots（p）是一个包含多项式根的列向量。

相反的，如果r是一个包含多项式根的列向量，那么，poly（r）是一个包含降幂排列多项式系数的行向量，见图2-65。

我们可以用上述的roots（）函数计算多项式p（s）的根，但是当多项式有重根时，函数roots1（）能给出更精确的结果。

  >>p=[1304];

   >>r=roots（p）

   r=-3.3553e+00

   1.7765e-01+1.0773e+00j

   1.7765e-01-1.0773e+00j

   >>p=poly（r）

   p=1.0003.0000.000-0.000j4.000+0.000j

   图2-65输入多项式并求根

矩阵乘法由MATLAB的conv（）函数完成。

假设我们要把两个多项式相乘合并成一个多项式n（s），即

n（s）=（3s2+2s+1）（s+4）

 =3s3+14s2+9s+4

    与此运算相关的MATLAB函数就是conv（），如图2-66所示。

函数polyval（）用来计算多项式的值。

多项式n（s）在s=-5处的值为n（-5）=-66，见图2-66。

 >>p=[321];q=[14];

 >>n=conv（p,q）

   n=31494

 >>value=polyval（n,-5）

   value=-66

      图2-66MATLAB的conv（）函数和polyval（）函数

2.5.2.传递函数

下面我们将介绍如何利用MATLAB函数获得传递函数在复平面内的零极点分布图。

设传递函数为G（s）=num/den，其中num和den均为多项式。

利用函数

[P,Z]=pzmap（num,den）

可以获得G（s）的零极点位置，即P为极点位置列向量，Z为零点位置列向量。

该指令执行后自动生成零极点分布图。

考虑传递函数

  和  

利用一系列MATLAB指令和函数，我们可以计算传递函数的零极点、特征方程和两个传递函数相除

，还可以在复平面上获得G（s）/H（s）的零极点分布图。

　　传递函数G（s）/H（s）的零极点图如图2-67所示，相应的MATLAB指令如图2-68所示。

图2-67中清楚的表示出了5个零点的分布，但是看上去只有两个极点，实际上这是不可能的，因为我们知道极点的个数一定大于或等于零点的个数。

应用函数roots1（）可以看到，实际上有四个极点位于s=-1。

因此，在同一位置上重复的极点或零点在零极点分布图上是区分不出的。

图2-67零极点图

2.5.3.结构图模型

假设我们已为某系统的各部分建立了传递函数，下一步的任务就是利用MATLAB函数将这些部分联接起来构成一个闭环控制系统，实现结构图的转换，计算从输入R（s）到输出C（s）的传递函数。

为了方便介绍MATLAB函数，我们从结构图的基本变换开始。

一个简单的开环控制系统可以通过G1（s）与G2（s）两个环节的串联而得到，利用series（）函数可以求串联连接的传递函数，函数的具体形式为

[num,den]=series（num1,den1,num2,den2）

例如G1（s）和G2（s）的传递函数分别为

，   

则

串联函数的用法示于图2-69。

>>numg=[601];deng=[1331];

>>z=roots（numg）

z=

 0+0.4082j

0-0.4082j

>>p=roots1（deng）

p=

 -1

 -1

 -1

>>n1=[11];n2=[12];d1=[12*j];d2=[1–2*j];d3=[13];

>>numh=conv（n1,n2）;denh=conv（d1,conv（d2,d3）;

>>num=conv（numg,denh）;den=conv（deng,numh）;

>>printsys（num,den）

num/den=

  6s^5+18s^4+25s^3+75s^2+4s+12   ａ

  s^5+6s^4+14s^3+16s^2+9s+2

>>pzmap（num,den）

>>title（‘Pole-ZeroMap’）

图2-68绘制零极点图指令 

>>num1=[1];den1=[50000];

>>num2=[11];den2=[12];

>>[num,den]=series（num1,den1,num2,den2）;

>>printsys（num,den）

num/den=

        s+1    

  500s^3+1000s^2

图2-69series函数的用法

当系统是以并联的形式连接时，利用parallel（）函数可得到系统的传递函数。

指令的具体形式为

[num,den]=parallel（num1,den1,num2,den2）

如果系统以反馈方式构成闭环，则系统的闭环传递函数为



求闭环传递函数的MATLAB函数有两个：

cloop（）和feedback（），其中cloop（）函数只能用于H（s）=1（即单位反馈）的情况。

cloop（）函数的具体用法为

[num,den]=cloop（numg,deng,sign）

其中numg和deng分别为G（s）的分子和分母多项式，sign=1为正反馈，sign=-1为负反馈（默认值）。

                                    

图2-70闭环反馈系统的结构图  

feedback（）函数的具体用法为

[num,den]=feedback（numg,deng,numh,denh,sign）

其中numh为H（s）的分子多项式，denh为分母多项式。

假设闭环反馈系统的结构图如图2-70所示，被控对象G（s）和控制部分Gc（s）以及测量环节H（s）的传递函数分别为

，

，

应用series（）函数和feedback（）函数求解闭环传递函数的MATLAB指令如图2-71所示。

>>numg=[1];deng=[500];

>>numc=[11];denc=[12];

>>numh=[1];denh=[110];

>>[num1,den1]=series（numc,denc,numg,deng）;

>>[num,den]=feedback（num1,den1,numh,denh,-1）;

>>printsys（num,den）

num/den=

        s^2+11s+10    

  5s^4+60s^3+100s^2+s+1

图2-71feedback（）函数的应用

例2.12 一个多环的反馈系统如图2-49所示，给定各环节的传递函数为

，

，

，

，

，

试求闭环传递函数GB（s）=C（s）/R（s）。

解 求解过程可按如下步骤进行：

步骤1：

输入系统各环节的传递函数

步骤2：

将H2的综合点移至G2后

步骤3：

消去G3，G2，H2环

步骤4：

消去包含H3的环

步骤5：

消去其余的环，计算GB（s）

根据上述步骤的MATLAB指令以及计算结果在图2-72中。

>>ng1=[1];dg1=[110];

>>ng2=[1];dg2=[11];

>>ng3=[101];dg3=[144];

>>ng4=[11];dg4=[16];

>>nh1=[1];dh1=[1];

>>nh2=[2];dh2=[1];

>>nh3=[11];dh3=[12];

>>[n1,d1]=series（ng2,dg2,nh2,dh2）;

>>[n2,d2]=feedback（ng3,dg3,n1,d1,-1）;

>>[n3,d3]=series（n2,d2,ng4,dg4）;

>>[n4,d4]=feedback（n3,d3,nh3,dh3,-1）;

>>[n5,d5]=series（ng1,dg1,ng2,dg2）;

>>[n6,d6]=series（n5,d5,n4,d4）;

>>[n7,d7]=cloop（n6,d6,-1）;

>>printsys（n7,d7）

num/den=

               s^4+3s^3+3s^2+3s+2          

  2s^6+38s^5+261s^4+1001s^3+1730s^2+1546s+732

图2-72多环结构图简化

有时我们需要关心闭环传递函数是否有零极点对消的情况出现。

当然通过pzmap（）或roots（）函数可以查看传递函数是否有相同的零极点，另外还可以使用minreal（）函数除去传递函数共同的零极点因子。

正如图2-73所示，如果传递函数有相同的零极点，应用minreal（）函数后，传递函数的分子和分母多项式各减少了一阶，消去了相同的零极点。

>>numg=[16116];deng=[1712115];

>>printsys（numg,deng）

numg/deng=

      s^3+6s^2+11s+6  

   s^4+7s^3+12s^2+11s+5

>>[num,den]=minreal（numg,deng）;

>>printsys（num,den）

1pole-zeroscancelled

num/den=

     s^2+4s+3 

   s^3+6s^2+6s+5

图2-73minreal（）函数的应用

最后重新考虑例2.2所示的位置随动系统，我们的目的是计算闭环系统在输入作用下的响应。

在给定各元件参数并忽略La和令ML=0的情况下，其结构图如图2-74所示。

计算的第一步是求闭环传递函数GB（s）=θc（s）/θr（s），求解过程及结果如图2-75所示。

第二步是利用step（）函数计算参考输入θr（t）为单位阶跃信号时输出θc（t）的响应。

可见，特征方程是2阶的，且ωn=52，ζ=0.012，由于阻尼比很小，可以预料响应会强烈震荡。

　　　　



　　　　图2-74位置随动系统的结构图

　　　　　



　　　　　图2-76位置随动系统的阶跃响应曲线

>>num1=[200];den1=[20];num2=[1];den2=[20.50];

>>num3=[0.20];den3=[1];num4=[540];den4=[1];

>>[na,da]=series（num1,den1,num2,den2）;

>>[nb,db]=feedback（na,da,num3,den3,-1）;

>>[nc,dc]=series（nb,db,num4,den4）;

>>[num,den]=cloop（nc,dc,-1）;

>>printsys（num,den）

num/den=

        5400    

2s^2+2.5s+5400

>>t=[0:

0.005:

3];

>>[y,t]=step（num,den,t）;

>>plot（t,y）,grid

图2-75位置随动系统的结构图简化及阶跃响应指令

图2-76给出了位置随动系统的阶跃响应曲线。

y（t），即θc（t）的离散时间点t将以行向量的形式给出，它从0秒开始，按0.005秒的步长增加，直到3秒为止。

使用plot（）函数用于画出y（t）曲线，grid函数用于给图形加上网格。

由于控制系统的性能指标通常以阶跃响应的形式给出，因此step（）函数是非常重要的。

有关step（）函数的内容将在后续章节中进一步介绍。

3.4应用MATLAB分析控制系统的性能

这一节将用两个例子描述反馈控制的优点，同时说明如何利用MATLAB来分析控制系统。

系统分析的主要内容包括如何抑制干扰、如何减小稳态误差、如何调节瞬态响应以及如何减少系统对参数变化的影响等。

    第一个例子是带有负载转矩干扰信号的电枢控制直流电动机。

开环系统结构图如图3-37（a）所示，为了改善系统性能，加入速度反馈如图3-37（b）所示。

系统的各元器件参数值在表3-6中给出。

表3-6 速度控制系统的参数

参数名

Ra

Km

J

B

Ke

Ka

Ks

参数值

1

10

2

0.5

0.1

54

1

图3-37速度控制系统结构图

从图中可以看出，系统有Ua（s）（或Vr（s）和ML（s）两个输入。

由于这是一个线性系统，按叠加定理可以分别考虑两个输入的独立作用结果。

为了研究干扰对系统的作用，可令Ua（s）=0（或Vr（s）=0），此时只有干扰ML（s）起作用。

相反地，为了研究参考输入对系统的响应，可令ML（s）=0。

如果系统具有很好的抗干扰能力，则干扰信号ML（s）对输出ω（s）的影响就应该很小，下面就来验证此结论。

　　首先，考虑图3-37（a）所示的开环系统，从ML（s）到ωo（s）（此处的下标“o”表示开环）的传递函数为

假设干扰信号为单位阶跃信号，即ML（s）=1/s。

利用MATLAB可以计算系统的单位阶跃响应如图3-38（a）所示，而用于分析此开环控制系统的MATLAB程序文本opentach.m示于图3-38（b）。

在输入信号Ua（s）=0的情况下，稳态误差就是干扰响应ωo（t）的终值。

在图3-38（a）的曲线中，干扰响应ωo（t）在t=7秒后已近似不变，所以近似稳态误差值为

ωo（∞）≈-0.663（弧度/秒）

同样，通过计算从ML（s）到ωc（s）（此处下标“c”表示闭环）的闭环传递函数可分析图3-37（b）所示闭环系统的抗干扰性能。

对于干扰输入的闭环传递函数为

（a） 开环速度系统对阶跃干扰的响应曲线

（b）   MATLAB程序文本:

opentach.m

图3-38开环速度控制系统分析

闭环系统对单位阶跃干扰输入的响应曲线ω（t）和MATLAB程序文本closedtach.m分别示于图3-39（a）（b）。

同前，稳态误差就是ω（t）的终值，稳态误差的近似值为

在本例中，闭环系统与开环系统对单位阶跃干扰信号的输出响应的稳态值之比为

可见通过引入负反馈已明显减小了干扰对输出的影响，这说明闭环反馈系统具有抑制噪声特性。

（a） 闭环系统对阶跃干扰的响应曲线

（b）   MATLAB程序文本:

opentach.m

图3-39闭环速度控制系统分析

第二个例子是分析闭环控制系统的控制器增益K对瞬态响应的影响。

图3-40是闭环控制系统的结构图。

在参考输入R（s）和干扰输入N（s）同时作用下系统的输出为

图3-40反馈控制系统的结构图

如果单纯考虑增益K对参考输入产生的瞬态响应的影响，可以预计增加K将导致超调量增加、调整时间减少和响应速度提高。

在增益K=20和K=100时，系统对参考输入的单位阶跃响应曲线以及相应的MATLAB程序文本gain_kr.m示于图3-41。

对比两条响应曲线，可以看出上述预计的正确性。

尽管在图中不能明显看出增大K能减少调整时间，但是这一点可以通过观察MATLAB程序的运行数据得以验证。

这个例子说明了控制器增益K是如何改变系统瞬态响应的。

根据以上分析，选择K=20可能是一个比较好的方案。

尽管如此，在做出最后决定之前还应该考虑其他因素。

（a） 阶跃响应曲线

%K=20和K=100时，参考输入的单位阶跃响应：

gain_kr.m

numg=[1];deng=[110];

K1=100;K2=20;

num1=[11K1];num2=[11K2];den=[01];

%简化结构图

[na,da]=series（num1,den,numg,deng）;

[nb,db]=series（num2,den,numg,deng）;

[numa,dena]=cloop（na,da）;

[numb,denb]=cloop（nb,db）;

%选择时间间隔

t=[0:

0.01:

2.0];

[c1,x,t]=step（numa,dena,t）;

[c2,x,t]=step（numb,denb,t）;

plot（t,c1,'--',t,c2）

xlabel（'Time[sec]'）,ylabel（'Cr（t）'）,grid

（b）  MATLAB程序文本:

gain_kr.m

图3-41单位阶跃输入的响应分析

在对K做出最后选择之前，非常重要的是要研究系统对单位阶跃干扰的响应，有关结果和相应的MATLAB程序文本如图3-42所示。

从中可以看到，增加K减少了单位干扰响应的幅值。

对于K=20和K=100，响应的稳态值分别为0.05和0.01。

对干扰输入的稳态值可按终值定理求得

如果仅从抗干扰的角度考虑，选择K=100更合适。

在本例中所求出的稳态误差、超调量和调整时间（2%误差）归纳于表3-7。

在控制系统设计中有很多成熟的经验，设计者常常要权衡利弊，综合考虑。

在这个例子中，增加K导致了更好的抗干扰性，然而减少K可以使系统具有更好的瞬态响应性能。

如何选择K的最终权力留给了设计者。

尽管MATLAB软件对控制系统的分析和设计很有帮助，但是控制工程师的经验往往更重要。

（a） 阶跃响应曲线

%K=20和K=100时，干扰输入的单位阶跃响应：

gain_kn.m

numg=[1];deng=[110];

K1=100;K2=20;

num1=[11K1];num2=[11K2];den=[01];

%简化