第五节
应用 MATLAB控制系统仿真
提纲
?一、弹簧 -重物 -阻尼器系统
?二、传递函数
?三、结构图模型
引言
MATLAB是一套高性能的数值计算和可视化软件,
它集数值分析、矩阵运算和图形显示于一体,构成了
一个方便的界面友好的用户环境。
控制系统的分析与设计方法,都是以数学模型为基
础进行的。 MATLAB可以用于以传递函数形式描述的
控制系统。
在本节中, 首先举例说明如何使用 MATLAB进行辅
助分析 。 然后讨论传递函数和结构图 。
一、弹簧 -重物 -阻尼器系统
弹簧 —重物 —阻尼器动力学系统如图 2-1所示 。 重物 M的位
移由 y(t)表示, 用微分方程描述如下:
2
2
d ( ) d ( ) ( ) ( )
dd
y t y tM B Ky t f t
tt? ? ?
该系统在初始位移作用下的瞬态响应为:
2
2
(0 )( ) s i n ( 1 )
1
n t
n
yy t e t?? ? ? ?
?
?? ? ?
?
其中 ??cos-1?, 初始位移是 y(0)。
系统的瞬态响应当 ?< 1时为欠阻尼, 当 ?>1时为过阻尼,
当 ?= 1时为临界阻尼 。
过阻尼情况,y(0)=0.15 m ?n= (弧度 /秒 )
( )
欠阻尼情况,y(0)=0.15 m ?n= (弧度 /秒 )
( )
2
1
3
22? ?
2,3KBMM??
2
2
1
22? ?
2,1KBMM??
利用 MATLAB程序 —unforced.m,可以显示初始位移为 y(0)
的物体自由运动曲线,如图 2-63所示。
在 unforced.m程序中,变量 y(0),?n,t,? 1和 ? 2的值由指令
直接输入工作区,然后运行 unforced.m程序就可以产生响应曲
线。
>>y0=0.15;wn=sqrt(2);
>>zeta1=3/(2*sqrt(2)); zeta2=1/(2*sqrt(2));
>>t=[0:0.1:10];
>>unforced
( a) MATLAB指令窗口
* 计算系统在给定初始条件下的自由运动
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,[‘oeverdamped zeta1=’,num2str(zeta1),] )
text(0.2,0.80,[‘underdamped zeta2=’,num2str(zeta2),] )
( b)分析弹簧 —重物 —阻尼器的 MATLAB程序 unforced.m
图 2-63 分析弹簧 —重物 —阻尼器的 MATLAB指令
图 2-64 弹簧 —重物 —阻尼器的自由运动曲线
在欠阻尼和过阻尼情况下的响应曲线如图 2-64所示,
MATLAB可分析以传递函数形式描述的系统 。 分子多项式和
分母多项式都必须在 MATLAB指令中指定 。
在 MATLAB中多项式由行向量组成,这些行向量包含了降次
排列的多项式系数。例如多项式 p(s)=1s3+3s2+0s1+4s0,按图 2-65
的格式输入 p=[1 3 0 4],
>>p=[1 3 0 4];
>>r=roots(p)
r=
-3.3553e+00
1.7765e-01+1.0773e+00j
1.7765e-01-1.0773e+00j
>>p=poly(r)
p=
1.000 3.000 0.000-0.000j 4.000+0.000j
图 2-65 输入多项式并求根
矩阵乘法由 MATLAB的 conv()函数完成 。 把两个多项式相乘
合并成一个多项式 n(s),即:
n(s)= (3s2 +2s +1) (s +4) = 3s3 +14s2 +9s +4
与此运算相关的 MATLAB函数就是 conv()。 函数 polyval()用来计
算多项式的值 。 多项式 n(s)在 s = -5处值为 n (-5) = -66,见图 2-66。
>>p=[3 2 1];q=[1 4];
>>n=conv(p,q)
n=
3 14 9 4
>>value=polyval(n,-5)
value=
-66
图 2-66 MATLAB的 conv()函数和 polyval()函数
设传递函数为 G(s)=num/den,其中 num和 den均为多项式 。 利
用函数:
二、传递函数
[P,Z]=pzmap(num,den)
可得 G(s)的零极点位置, 即 P为极点位置列向量, Z为零点位
置列向量 。 该指令执行后自动生成零极点分布图 。
考虑传递函数:
2
32
61()
3 3 1
sGs
s s s
??
? ? ?
和 ( 1 ) ( 2)()
( 2 ) ( 2 ) ( 3 )
ssHs
s i s i s
???
? ? ?
图 2-67零极点图
传递函数 G(s)/H(s)的零极点图如图 2-67所示,相应的 MATLAB
指令如图 2-68所示。
>>numg=[6 0 1];deng=[1 3 3 1];
>>z=roots(numg)
z=
0+0.4082j
0-0.4082j
>> p=roots1(deng)
p=
-1
-1
-1
>>n1=[1 1]; n2=[1 2]; d1=[1 2*j]; d2=[1 –2*j]; d3=[1 3];
>>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+
图 2-68 绘制零极点图指令
三,结构图模型
一个开环控制系统可以通过 G1 (s)与 G2 (s)两个环节的
串联而得到, 利用 series()函数可以求串联连接的传递
函数, 函数的具体形式为:
[num,den]= series(num1,den1,num2,den2)
例如 G1 (s)和 G2 (s)的传递函数分别为:
1
1()
2
sGs
s
??
? 2 2
1()
500Gs s?
12 2 3 2
( ) 1 1 1( ) ( )
( ) 2 5 0 0 5 0 0 1 0 0 0
C s s sG s G s
R s s s s s
??? ? ? ?

串联函数的用法示于图 2-69:
>>num1=[1];den1=[500 0 0];
>>num2=[1 1];den2=[1 2];
>>[num,den]=series(num1,den1,num2,den2);
>>printsys(num,den)
num/den=
s+1
500s^3+1000s^2
图 2-69 series函数的用法
当系统是以并联的形式连接时, 利用 parallel()函数可得到系
统的传递函数 。 指令的具体形式为:
系统以反馈方式构成闭环, 则系统的闭环传递函数为:
()()
1 ( ) ( )B
GsGs
G s H s? ?
[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为负反馈(默认值)。
feedback()函数的用法为:
[num,den]= feedback (numg,deng,numh,denh,sign)
其中 numh为 H (s)的分子多项式, denh为分母多项式 。
闭环反馈系统的结构图如图 2-70所示, 被控对象 G(s)和控制
部分 Gc (s)以及测量环节 H (s)的传递函数分别为:
,,
1()
2c
sGs
s
???
?
numc
denc 2
1()
5Gs s??
numg
deng
1()
10Hs s?? ?
numh
denh
图 2-70 闭环反馈系统的结构图
应用 series()函数和 feedback()函数求闭环传递函数的 MATLAB
指令如图 2-71 所示:
>>numg=[1];deng=[5 0 0];
>>numc=[1 1];denc=[1 2];
>>numh=[1];denh=[1 10];
>>[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-71 feedback()函数的应用
例 2.12 一个多环的反馈系统如图 2-49所示, 给定各环节的传
递函数为:
1
1()
10Gs s? ? 2
1()
1Gs s? ?
2
3 2
1()
44
sGs
ss
??
?? 4 1() 6sGs s ?? ?
1 ( ) 1Hs ? 2 ( ) 2Hs ?
3
1()
2
sHs
s
??
?
试求闭环传递函数 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=[1 10];
>>ng2=[1];dg2=[1 1];
>>ng3=[1 0 1];dg3=[1 4 4];
>>ng4=[1 1];dg4=[1 6];
>>nh1=[1];dh1=[1];
>>nh2=[2];dh2=[1];
>>nh3=[1 1];dh3=[1 2];
>>[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所示 。
>>numg=[1 6 11 6];deng=[1 7 12 11 5];
>>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)
1 pole-zeros cancelled
num/den=
s^2+4s+3
s^3+6s^2+6s+5
图 2-73 minreal()函数的应用
例 2.2所示的位置随动系统,在给定各元件参数并忽略 La和令
ML = 0的情况下,其结构图如图 2-74所示:
图 2-74 位置随动系统的结构图
第一步求闭环传递函数 GB (s)=?c(s) /?r(s),求解过程及结果
如图 2-75所示。第二步利用 step()函数计算参考输入 ?r (t)为单位
阶跃信号时输出 ?c (t)的响应。
>>num1=[200]; den1=[20];num2=[1]; den2=[2 0.5 0];
>>num3=[0.2 0]; 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 位置随动系统的阶跃响应曲线
图 2-76给出了位置随动系统的阶跃响应曲线。用 plot()函数用于
画出 y(t)曲线,grid函数用于给图形加上网格。
小结
本章讨论了如何建立控制系统以及元部件的数学模
型问题。本章所涉及的数学模型共有三种,即微分方
程、传递函数、结构图或信号流图。
利用传递函数研究线性系统, 可根据传递函数的极点
和零点分布判定系统对不同输入信号的响应特性 。
利用结构图或信号流图可以了解系统中的每个变量,
还可以通过梅逊 ( Mason) 公式, 方便地求得系统输入
输出间的传递函数 。
利用 MATLAB软件可求解系统在不同参数和输入情况
下的响应 。