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