MATLAB 辅助教学.docx
《MATLAB 辅助教学.docx》由会员分享，可在线阅读，更多相关《MATLAB 辅助教学.docx（74页珍藏版）》请在冰豆网上搜索。
MATLAB辅助教学
MATLAB辅助教学
转到：
第一章第二章第三章第四章第五章第六章第七章第八章返回目录
第一章控制系统的数学模型
例1 阻尼系统如图所示。
c是阻尼系数，k是弹簧系数，m是运动
体的方程，f（t）是外作用力，x（t）是物体的机械位移，该系统原来处
于静止状态。
试列出系统的运动方程，并求其传递函数。
弹簧的弹性力为kx（t）,阻尼器阻力为
。
根据牛顿运动定律，可以得到该系统
的运动方程为：
在MATLAB中，建立此系统的传递函数的程序代码：
symsmkc
s=tf（'s'）;
G='1/（m*s^2+c*s+k）'
可得到该系统的闭环传递函数为：
例2 计算时域函数
的拉氏变换，并计算函数
的拉氏逆变换。
程序代码：
symstay;
y=laplace
（（t^2）*exp（a*t））,
symssF;
F=ilaplace（10/（s^2+4）-3*s/（s^2+4）+（s+4）/（（s+4）^2+16））
例3 已知一系统的传递函数
，提取其分子和
分母多项式，并绘制零极点图。
程序代码：
num=[328];
den=[13842];
G=tf（num,den）,
[tt,ff]=tfdata（G,'v'）%提取分子和分母多项式
pzmap（G）%绘制零极点图
gridon%打开绘图网格
例4 求一传递函数为
的零极点及其增益，
并绘制零极点图。
程序代码：
num=[328];
den=conv（[1612],[124]）;
G=tf（num,den）;
[z,p,k]=zpkdata（G,'v'）
pzmap（G）%绘制零极点图
axis（[-3.2,0,-2,2]）%设定零极点图坐标
gridon
例5已知控制系统的微分方程式为：
试
建立状态空间模型。
由本题可知
，
。
程序代码：
A=[010;001;-3-5-9];
B=[1;-5;41];
C=[100];
D=0;
G=ss（A,B,C,D）
例6
，
，求其TF模
型和零极点模型，并求其零极点和绘制零极点图。
程序代码：
A=[0100;0010;0001;-50-48-28.5-9];
B=[0;0;0;1];
C=[10200];
D=0;
Gss=ss（A,B,C,D）
Gtf=tf（Gss）
Gzpk=zpk（Gss）
[z,p,k]=zpkdata（Gss,'v'）
pzmap（Gss）
gridon
例7 已知一串联系统的三个传递函数
、
和
，求此系统的传递函数。
程序代码：
G1=tf（[265],[1452]）;
G2=tf（[141],[1980]）;
z=[-3-7];
p=[-1-4-6];
k=5;
G3=zpk（z,p,k）;
G=G3*G2*G1
例8 已知两状态方程描述的系统进行并联连接，其中状态系数矩阵为：
，
，
，
；
，
，
，
。
求此并联系统的
状态空间模型。
程序代码：
G1=ss（[-221;0-24;1-40],[0;0;-1],[1-11],1）;
G2=ss（[010;001;-5-3-2],[0;0;1],[1.510.5],2）;
G=G1+G2
例9 如图所示的连接框图，其中
，
，
，
，求
此系统总的传递函数。
程序代码：
G1=tf（[113],[1421]）;
G2=tf（[567],[1234]）;
G3=tf（[18171],[16116]）;
G4=tf（[56],[11080128]）;
G=（G2*G1+G3）*G4
例10 如图所示为控制系统的控制框图，求系统的闭环传递函数。
程序代码：
G1=tf（[10],[11]）;
G2=tf（[2],conv（[10],[11]））;
H2=tf（[12],[13]）;
H1=tf（[50],conv（[12],[14]））;
GH=feedback（G2,H2,+1）;
Gc=GH*G1;G=feedback（Gc,H1）
第二章时域分析法
例1 计算有初始条件的零输入响应。
设控制系统状空间表达是：
，
。
零输入时的状态值
。
程序代码：
A=[-2-2.5-0.5;100;010];
B=[100]';
C=[01.51];
D=0;
G=ss（A,B,C,D）;
Gs=tf（G）
y=initial（G,[102],0:
2:
100）;
h=plot（y）;
set（h,'linewidth',3*get（h,'linewidth'））;
grid;axis（[0100-0.22.5]）;
title（'零输入响应曲线'）;
xlabel（'t'）;
ylabel（'y'）;
可以看到，系统输出的初始值不为零；在本题中，输出为2。
在程序中增加零
状态输出就可得到实际的系统输出。
程序代码如下：
y1=step（G,100）;
yy=y1+y;
plot（[yy1yy]）;
grid;title（'响应输出曲线'）;
xlabel（'t'）;
ylabel（'y'）;
axis（[0120-0.23]）;
text（40,0.4,'零输入响应'）
text（40,0.7,'零状态输出响应'）
text（12,2.5,'总的输出响应'）
例2 某位置随动系统的方框图如图所示。
求此系统的单位阶跃响应曲线。
先求系统的传递函数，代码如下：
G=tf（[10],conv（[10],[0.51]））;
H=tf（[0.11],[1]）;
GH1=feedback（G,H）;
GH=feedback（GH1,1）
绘制阶跃响应曲线，代码如下：
GH=tf（[10],[0.5220]）;
step（GH）
title（'单位阶跃响应曲线'）;
gridon
例3 已知某控制系统的传递函数
，对于任意的
输入信号，求系统的输出响应曲线。
1 当输入信号是
时；
2 当输入信号是
时。
绘制正弦信号输出响应曲线的程序代码：
G=tf（conv（[5],[11]）,conv（[10],[1423]））;
t=[0:
.1:
150]';
u=sin（t+30*pi/180）;
y=lsim（G,u,t）;
plot（y（1:
150））;
set（gca,'ytick',-4:
1:
8）;
title（'正弦信号的输出响应曲线'）;
xlabel（'t'）;
ylabel（'y'）;
grid
绘制余弦信号输出响应曲线的程序代码：
G=tf（conv（[5],[11]）,conv（[10],[1423]））;
t=[0:
0.1:
150]';
u=2*cos（5*t+30*pi/180）;
y=lsim（G,u,t）;
plot（y（1:
150））;
title（'余弦信号的输出响应曲线'）;
xlabel（'t'）;
ylabel（'y'）
gridon
例4 已知单位负反馈控制系统的开环传递函数为
，试判断此闭环系统的稳定性。
首先求闭环系统的特征多项式，程序代码：
z=-2;
p=[0-0.5-0.8-3];
k=0.2;
Go=zpk（z,p,k）;
Gc=tf（feedback（Go,1））;
dc=Gc.den;dens=poly2str（dc{1},'s'）
根据特征多项式，求其特征根来判断系统的稳定性。
程序如下：
den=[14.34.31.40.4];
p=roots（den）
利用零极点图来判断系统的稳定性。
程序如下：
[z,p,k]=zpkdata（Gc,'v'）
pzmap（Gc）
例5 如图所示的控制系统方框图，计算其动态误差系数。
程序代码：
num=[112];
den=10;
G=tf（num,den）;
Es=feedback（1,G）;
symssnumsdens
nums=poly2sym（Es.num{1}）;
dens=poly2sym（Es.den{1}）;
y=sym2poly（taylor（nums/dens））;
y1=length（y）;
fori=1:
y1
yy（i）=y（y1-i+1）;
ifyy（i）==0
k（i）=inf;
else
k（i）=1/yy（i）;
end
end
k
例6已知二阶振荡环节的传递函数
，其中ωn＝
0.4，ζ从0变换到2，求此系统的单位阶跃响应曲线、脉冲响应曲线
和斜坡响应曲线。
系统单位阶跃响应曲线的程序代码：
symss
forzeta=[0:
0.2:
0.8,1:
0.5:
2]
wn=0.4;wn=sym（num2str（wn））;zet=sym（num2str（zeta））;
ifzeta==0
figure
（1）
ezplot（ilaplace（wn^2/s/（s^2+wn^2））,[080]）;
gridon
title（'\xi=0'）
elseifzeta==1
figure
（2）
ezplot（ilaplace（wn^2/s/（s+wn）^2）,[080]）;
holdon;
else
figure
（2）
ezplot（ilaplace（wn^2/s/（s^2+2*zet*wn*s+wn^2））,[080]）;
holdon;
end
end
gridon
title（'\xi:
0.2,0.4,0.6,0.8,1.0,1.5,2.0'）
axis（[08001.8]）
gtext（'0.4'）
gtext（'1.0'）
gtext（'2.0'）
系统脉冲响应曲线的程序代码：
symss
forzeta=[0:
0.2:
0.8,1:
0.5:
2]
wn=0.4;wn=sym（wn）;zet=sym（zeta）;
ifzeta==0
figure
（1）
ezplot（ilaplace（wn^2/（s^2+wn^2））,[060]）;
gridon
title（'\xi=0'）
elseifzeta==1
figure
（2）
ezplot（ilaplace（wn^2/（s+wn）^2）,[060]）;
holdon;
else
figure
（2）
ezplot（ilaplace（wn^2/（s^2+2*zet*wn*s+wn^2））,[060]）;
holdon;
end
end
gridon,title（'\xi:
0.2,0.4,0.6,0.8,1.0,1.5,2.0'）
axis（[060-0.20.4]）
gtext（'0.4'）
gtext（'1.0'）
gtext（'2.0'）
系统斜坡响应的程序代码：
当输入信号为斜坡信号时，R（s）=
。
symss
forzeta=[0:
0.2:
0.8,1:
0.5:
2]
wn=0.4;wn=sym（wn）;zet=sym（zeta）;
ifzeta==0
figure
（1）
ezplot（ilaplace（wn^2/s^2/（s^2+wn^2））,[060]）;
gridon
title（'\xi=0'）
elseifzeta==1
figure
（2）
ezplot（ilaplace（wn^2/s^2/（s+wn）^2）,[080]）;
holdon;
else
figure
（2）
ezplot（ilaplace（wn^2/s^2/（s^2+2*zet*wn*s+wn^2））,[080]）;
holdon;
end
end
gridon
title（'\xi:
0.2,0.4,0.6,0.8,1.0,1.5,2.0'）
axis（[080080]）
例7 已知控制系统的传递函数
，求此系统单位阶
跃响应曲线和脉冲响应曲线。
单位阶跃响应曲线程序代码：
num=[124];
den=[11254];
[r,p,k]=residue（num,den）;
symss
y=sum（ilaplace（r./（s*（s-p））））;
ezplot（y,[030]）;
axis（[03001.4]）;
ylabel（'y（t）'）;
title（'单位阶跃响应曲线'）;
grid
脉冲响应曲线程序代码：
num=[124];
den=[11254];
[r,p,k]=residue（num,den）;
symss
y=sum（ilaplace（r./（s-p）））;
ezplot（y,[025]）;
axis（[025-0.41.0]）;
ylabel（'y（t）'）;
title（'脉冲响应曲线 '）;
gridon
例8 求一阶惯性环节
的阶跃响应曲线和脉冲响应曲线。
绘制阶跃响应曲线程序代码：
forT=5:
5:
30
G=tf（[1],[T1]）;
step（G,160）;
holdon;
end
grid
axis（[016001.1]）;
title（'T:
5~30阶跃响应曲线'）;
set（gca,'ytick',0:
.1:
1.1）;
绘制脉冲响应曲线程序代码：
forT=5:
5:
30
G=tf（[1],[T1]）;
impulse（G,100）;
holdon;
end
grid
axis（[0100-0.010.22]）;
title（'T:
5~30脉冲响应曲线'）;
例9 已知三阶系统
，将此展开为部
分分式展开式，并求此控制系统的脉冲响应曲线。
程序代码：
num=[2115];
den=[173260];
[r,p,k]=residue（num,den）
symss
y=sum（ilaplace（r./（s-p）））;
ezplot（y,[04]）;
axis（[04-12.2]）;
set（gca,'ytick',-0.9:
.3:
2.1）;
title（'脉冲响应曲线'）;
ylabel（'y（t）'）;
gridon
第三章 根轨迹法
例1 根轨迹的概念
程序代码：
%G（s）H（s）=50*/s（s+1）
ng=1;
dg=poly（[0-1]）;
GH=tf（50*ng,dg）
GH=zpk（GH）
pzmap（ng,dg）
title（'Zero-PolePlot','FontWeight','bold'）
r=rlocus（ng,dg）;
rlocus（ng,dg）;
holdon
plot（r,'LineWidth',2）,
title（'Root-LocusPlotofG（s）H（s）=50/s（s+1）','FontWeight','bold'）
例2 分支数、起点与终点
程序代码：
%G（s）H（s）=k*（s+1）/s（s+2）
ng=[11];
dg=poly（[0-2]）;
GH=tf（ng,dg）
GH=zpk（GH）
[p,z]=pzmap（ng,dg）;
pmm=length（p）;zmm=length（z）;
fori=1:
pmm
pre（i）=real（p（i））;
pim（i）=imag（p（i））;
end
fori=1:
zmm
zre（i）=real（z（i））;
zim（i）=imag（z（i））;
end
plot（pre,pim,'x','LineWidth',2）
holdon
plot（zre,zim,'o','LineWidth',2）,
title（'Zero-PolePlot','FontWeight','bold'）
[r,k]=rlocus（ng,dg）;
k1=[0:
.02:
16];
rlocus（ng,dg,k1）;
%plot（r,'LineWidth',2）
title（'Root-LocusPlotofG（s）H（s）=k*（s+1）/s（s+2）','FontWeight',
'bold'）
例3渐近线、对称于实轴
程序代码：
%G（s）H（s）=k*/s（s+1）（s+2）
ng=1;
dg=poly（[0-1-2]）;
GH=tf（ng,dg）
GH=zpk（GH）
%pzmap（ng,dg）
title（'Zero-PolePlot'）
%r=rlocus（ng,dg）;
rlocus（ng,dg,':
'）;
axis（[-3.51.5-33]）
holdon
%plot（r,':
r'）
[xx,yy,pp,nex]=asymptote（ng,dg）;
plot（xx,yy,'b','LineWidth',2）
gridon
title（'Root-LocusPlotof
G（s）H（s）=k*/s（s+1）（s+2）','FontWeight','bold'）,
text（-1.44,-0.2,['\sigma=',num2str（pp）],'FontWeight','bold'）%渐近
线与实轴交点
text（-1.44,-0.5,['\0=\pm',num2str（180/nex）,'\circ',',',num2str
（180）,'\circ'],'FontWeight','bold'）%渐近线与实轴夹角
例4 起始角与终止角、与虚轴交点
程序代码：
%G（s）H（s）=k*/s（s^2+2s+2）
ng=1;
dg=[1220];
%GH=tf（ng,dg）
GH=zpk（GH）
r=rlocus（ng,dg）;
rlocus（ng,dg）;
sgrid
holdon
plot（r,'LineWidth',2）,[xx,yy]=asymptote（ng,dg）;
plot（xx,yy）,
title（'Root-LocusPlotofG（s）H（s）=k*/s（s^2+2s+2）
）','FontWeight','bold'）
%title（'Root-LocusPlotofG（s）H（s）=k*/s（s^2+2s+2）（用鼠标定点确定K
！
）','FontWeight','bold'）
%text（-1.1,1.1,['\Theta=',num2str（-45）,'\circ']）,
%text（-1.1,-1.1,['\Theta=',num2str（45）,'\circ']）,
%[ng,dg]=dispK（ng,dg）;
例5 分离点、实轴上的根轨迹
程序代码：
%G（s）H（s）=10*/s（s+1）（s+2）
%ng=[14];
dg=poly（[0-1-2-3]）;
ng=1;
dg=poly（[0-1-2]）;
GH=tf（ng,dg）
GH=zpk（GH）
r=rlocus（ng,dg）;
rlocus（ng,dg）;
axis（[-3.51.5-33]）
holdon
plot（r,'LineWidth',2）,
[xx,yy]=asymptote（ng,dg）;
plot（xx,yy）
gridon,
title（'G（s）H（s）=k*/s（s+1）（s+2）（用鼠标定点确定K！
每次选点的第一点
按在窗口外则结束操作）'）,%,'FontWeight','bold'
dispD1（GH）;
%xx=get（gca,'Xlim'）;
yy=get（gca,'Ylim'）;
%text（xx
（1）+0.1,0.93*yy
（2）,'点开环极点，'）
text（xx
（1）+0.1,0.84*yy
（2）,'则中断点取。
'）,
[ng,dg]=dispK（ng,dg）;
例6 与平行于虚轴的直线对称情况
程序代码：
%G（s）H（s）=k*/s（s+4）（s^2+4s+20）
%ng=1;
rr=roots（[146]）;
ng=1;
rr=roots（[1420]）;dg=poly（[0-4rr']）;
GH=tf（ng,dg）
GH=zpk（GH）,
[p,z]=pzmap（ng,dg）;
pmm=length（p）;
fori=1:
pmm
pre（i）=real（p（i））;
pim（i）=imag（p（i））;
end
plot（pre,pim,'x','LineWidth',2）
holdon
title（'Zero-PolePlot','FontWeight','bold'）
r=rlocus（ng,dg）;
rlocus（ng,dg）;
holdon
plot（r,'LineWidth',2）,
axis（[-62-55]）,[xx,yy]=asymptote（ng,dg）;
plot（xx,yy,-xx-4,-yy）,
title（'Root-LocusPlotofG（s）H（s）=k*/s（s+4）（s^2+4s+20）',
'FontWeight','bold'）,
%dispD2（GH）;
[ng,dg]=dispK（ng,dg）;
例7带开环零点的二阶系统
程序代码：
%G（s）H（s）=k*（s+1）/s^2
ng=[1
升级会员 