控制系统计算机仿真作业.docx
- 文档编号:7305206
- 上传时间:2023-01-22
- 格式:DOCX
- 页数:31
- 大小:104.25KB
控制系统计算机仿真作业.docx
《控制系统计算机仿真作业.docx》由会员分享,可在线阅读,更多相关《控制系统计算机仿真作业.docx(31页珍藏版)》请在冰豆网上搜索。
控制系统计算机仿真作业
兰州理工大学
《控制系统计算机仿真》
上机报告Ⅰ
院系:
电气工程与信息工程学院
班级:
14级自动化3班
姓名:
孙悦
学号:
1405220323
时间:
2017年6月15日
电气工程与信息工程学院
《控制系统计算机仿真》上机实验任务书Ⅰ(2017)
一、上机实验内容及要求
1.matlab软件
要求利用课余时间熟悉掌握matlab软件的基本数值运算、基本符号运算、基本程序设计方法及常用的图形命令操作。
2.各章节仿真实验内容及要求
具体实验内容及要求请详见上机实验报告。
二、上机实验时间安排及相关事宜
1.依据课程教学大纲要求,上机实验学时共16学时,学生须在每次上机之前做好相应的准备工作,以确保在有限的机时内完成仿真实验要求的内容;
2.实验完成后按规定完成相关的仿真实验报告;
3.仿真实验报告请按有关样本制作并A4打印,侧面装订,作为成绩评定的一部分。
自动化系《控制系统计算机仿真》课程组
2017年3月
一、Matlab基础操作
1-1用MATLAB语言求下列系统的状态方程、传递函数、零极点增益和部分分式形式的模型参数,并分别写出其相应的数学模型表达式:
(1)
程序如下:
num=[7,24,24]
den=[10,35,50,24]
[A,B,C,D]=tf2ss(num,den)
系统的状态方程:
A=
-3.5000-5.0000-2.4000
1.000000
01.00000
B=
1
0
0
C=
0.70002.40002.4000
D=
0
零极点增益形式:
[Z,P,K]=tf2zp(num,den)
Z=
-1.7143+0.6999i
-1.7143-0.6999i
P=
-1.2973+0.9838i
-1.2973-0.9838i
-0.9053
K=
0.7000
部分分式:
[R,P,H]=residue(num,den)
R=
-0.0071-0.2939i
-0.0071+0.2939i
0.7141
P=
-1.2973+0.9838i
-1.2973-0.9838i
-0.9053
H=
[]
(2)
A=[2.25-5-1.25-0.5;2.25-4.25-1.25-0.25;0.25-0.5-1.25-1;1.25-1.75-0.25-0.75]
A=
2.2500-5.0000-1.2500-0.5000
2.2500-4.2500-1.2500-0.2500
0.2500-0.5000-1.2500-1.0000
1.2500-1.7500-0.2500-0.7500
>>B=[4;2;2;0]
B=
4
2
2
0
>>C=[0202]
C=
0202
>>D=0
D=
0
零极点增益形式:
>>[Z,P,K]=ss2zp(A,B,C,D)
Z=
-1.0000+1.2247i
-1.0000-1.2247i
-1.5000
P=
-0.5000+0.8660i
-0.5000-0.8660i
-1.5000+0.0000i
-1.5000-0.0000i
K=
4.0000
传递函数形式:
>>num=[04142215]
num=
04142215
>>den=[146.255.252.25]
den=
1.00004.00006.25005.25002.2500
部分分式:
>>[R,P,H]=residue(num,den)
R=
4.0000
-0.0000
0.0000-2.3094i
0.0000+2.3094i
P=
-1.5000
-1.5000
-0.5000+0.8660i
-0.5000-0.8660i
H=
[]
1-2用殴拉法matlab编程实现下列系统的输出响应
在
上,
时的数值解。
,
要求保留4位小数,并将结果以图形的方式与真解
比较。
t=0:
0.1:
1
t=
Columns1through9
00.10000.20000.30000.40000.50000.60000.70000.8000
Columns10through11
0.90001.0000
h=0.1;
y
(1)=1;
t=0:
0.1:
1;
h=0.1;
y
(1)=1;
fori=1:
10
y(i+1)=y(i)+h*(-1*y(i));
end
plot(t,y,'r')
holdon
m=exp(-1*t)
m=
Columns1through9
1.00000.90480.81870.74080.67030.60650.54880.49660.4493
Columns10through11
0.40660.3679
plot(t,m,'bo')
1-3用四阶龙格-库塔梯形法matlab编程实现1-2题的数值解,要求以图形的方式通过与真值及殴拉法的比较,分析其精度。
h=0.1;
y
(1)=1;
t=0:
0.1:
1
t=
Columns1through9
00.10000.20000.30000.40000.50000.60000.70000.8000
Columns10through11
0.90001.0000
fori=1:
10
k1=-1*y(i)
k2=-1*(y(i)+k1*h/2)
k3=-1*(y(i)+k2*h/2)
k4=-1*(y(i)+h*k3)
y(i+1)=y(i)+(k1+2*k2+2*k3+k4)*h/6
end
k1=
-1
k2=
-0.9500
k3=
-0.9525
k4=
-0.9047
y=
1.00000.9048
k1=
-0.9048
k2=
-0.8596
k3=
-0.8619
k4=
-0.8187
y=
1.00000.90480.8187
k1=
-0.8187
k2=
-0.7778
k3=
-0.7798
k4=
-0.7407
y=
1.00000.90480.81870.7408
k1=
-0.7408
k2=
-0.7038
k3=
-0.7056
k4=
-0.6703
y=
1.00000.90480.81870.74080.6703
k1=
-0.6703
k2=
-0.6368
k3=
-0.6385
k4=
-0.6065
y=
1.00000.90480.81870.74080.67030.6065
k1=
-0.6065
k2=
-0.5762
k3=
-0.5777
k4=
-0.5488
y=
1.00000.90480.81870.74080.67030.60650.5488
k1=
-0.5488
k2=
-0.5214
k3=
-0.5227
k4=
-0.4965
y=
1.00000.90480.81870.74080.67030.60650.54880.4966
k1=
-0.4966
k2=
-0.4718
k3=
-0.4730
k4=
-0.4493
y=
1.00000.90480.81870.74080.67030.60650.54880.49660.4493
k1=
-0.4493
k2=
-0.4269
k3=
-0.4280
k4=
-0.4065
y=
Columns1through9
1.00000.90480.81870.74080.67030.60650.54880.49660.4493
Column10
0.4066
k1=
-0.4066
k2=
-0.3862
k3=
-0.3873
k4=
-0.3678
y=
Columns1through9
1.00000.90480.81870.74080.67030.60650.54880.49660.4493
Columns10through11
0.40660.3679
plot(t,y,'o')
holdon
m=exp(-1*t)
plot(t,m,'r*')
lea=y-m
plot(t,lea,'g')
holdoff
m=
Columns1through8
1.00000.90480.81870.74080.67030.60650.54880.4966
Columns9through11
0.44930.40660.3679
lea=
1.0e-006*
Columns1through8
00.08200.14830.20130.24290.27470.29830.3149
Columns9through11
0.32560.33150.3332
1-4采用matlab语言编程实现
。
程序:
disp('y=')
h=1;
y=0;
fori=0:
1:
63
y=y+2^i;
end
disp(y);
运行结果:
y=
1.8447e+019
y=
1.8447e+019
1-5编写matlab的M-函数,以实现
。
要求在函数中给出必要的解释和说明,同时检测输入和返回变量的个数。
functionsum=zuoye5(m)
sum=0;
formatlong
fori=0:
m
t=1;
forj=1:
i
t=t*2;
end
sum=sum+t;
end
二、控制系统分析
2-1设典型闭环结构控制系统如下图所示,当阶跃输入幅值
时,用sp3_1.m求取输出y(t)的响应。
进一步考虑:
当反馈通道为
时,如何通过对sp3_1.m程序的修改,以实现此时所对应的y(t)。
解:
程序1:
clearall;closeall;
a=[0.0160.8643.273.421];
b=[3025];V=2;R=20;
X0=[0000];n=4;T0=0;Tf=10;h=0.01;
b=b/a
(1);a=a/a
(1);A=a(2:
n+1);
A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];
B=[zeros(1,n-1),1]';
m1=length(b);
C=[fliplr(b),zeros(1,n-m1)];
Ab=A-B*C*V;
X=X0';y=0;t=T0;
N=round(Tf-T0)/h;
fori=1:
N
K1=Ab*X+B*R;
K2=Ab*(X+h*K1/2)+B*R;
K3=Ab*(X+h*K2/2)+B*R;
K4=Ab*(X+h*K3)+B*R;
X=X+h*(K1+K2*2+K3*2+K4)/6;
y=[y,C*X];
t=[t,t(i)+h];
end
[t',y']
plot(t,y,'b-')
title('反馈系数为2时的阶跃响应曲线');
holdon;gridon
程序2:
clearall;closeall;
a=[0.0160.8643.273.421];b=[3025];
a1=conv(a,[0.11]);b1=b;
b1=b1/a1
(1);a1=a1/a1
(1);m1=length(b1);
b=b/a
(1);a=a/a
(1);m=length(b);
V=1;R=20;
T0=0;Tf=15;h=0.01;
X10=[00000];n1=5;
X0=[0000];n=4;
A1=a1(2:
n1+1);
A=a(2:
n+1);
C1=[fliplr(b1),zeros(1,n1-m1)];
C=[fliplr(b),zeros(1,n-m)];
A1=[rot90(rot90(eye(n1-1,n1)));-fliplr(A1)];
A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];
B1=[zeros(1,n1-1),1]';
B=[zeros(1,n-1),1]';
Ab=A1-B1*C1*V;
X1=X10';y1=0;t=T0;
N=round(Tf-T0)/h;
fori=1:
N;
K1=Ab*X1+B1*R;
K2=Ab*(X1+h*K1/2)+B1*R;
K3=Ab*(X1+h*K2/2)+B1*R;
K4=Ab*(X1+h*K3)+B1*R;
X1=X1+h*(K1+K2*2+K3*2+K4)/6;
y1=[y1,C1*X1];
t=[t,t(i)+h];
end
X=X0';y=0;t=T0;
fori=1:
N
K1=A*X+B*(R-y1(i));
K2=A*(X+h*K1/2)+B*(R-y1(i));
K3=A*(X+h*K2/2)+B*(R-y1(i));
K4=A*(X+h*K3)+B*(R-y1(i));
X=X+h*(K1+K2*2+K3*2+K4)/6;
y=[y,C*X];
t=[t,t(i)+h];
end
[t',y']
plot(t,y,'r-');
title('反馈通道为惯性环节时的阶跃响应曲线');
holdon;gridon
2-2下图中,若各环节传递函数已知为:
,
,
,
,
,
,
,
,
,
;试列写链接矩阵W、W0和非零元素阵WIJ,将程序sp4_2完善后,应用此程序求输出
的响应曲线。
解:
程序
clearall;closeall;
a=[101011011]';
b=[0.010.0850.010.0510.00670.1510.010.01]';
c=[1111700.211300.10.0044]';
d=[00.1700.1500000]';
P=[a,b,c,d];
WIJ=[101;211;29-1;321;431;
48-1;541;651;67-0.212;761;
861;971];
n=9;Y0=1;
Yt0=[000000000];
h=0.001;L1=10;T0=0;Tf=3;
nout=7;
A=diag(P(:
1));B=diag(P(:
2));
C=diag(P(:
3));D=diag(P(:
4));
m=length(WIJ(:
1));
W0=zeros(n,1);W=zeros(n,n);
fork=1:
m
if(WIJ(k,2)==0)
W0(WIJ(k,1))=WIJ(k,3);
else
W(WIJ(k,1),WIJ(k,2))=WIJ(k,3);
end
end
Q=B-D*W;
R=C*W-A;V1=C*W0;
Ab=Q\R;b1=Q\V1;
Y=Yt0';y=Y(nout);t=T0;
N=round((Tf-T0)/(h*L1));
fori=1:
N
forj=1:
L1
K1=Ab*Y+b1*Y0;
K2=Ab*(Y+h*K1/2)+b1*Y0;
K3=Ab*(Y+h*K2/2)+b1*Y0;
K4=Ab*(Y+h*K3)+b1*Y0;
Y=Y+h*(K1+K2*2+K3*2+K4)/6;
end
y=[y,Y(nout)];
t=[t,t(i)+h*L1];
end
[t',y']
plot(t,y,'r');
title('单位阶跃响应曲线');
gridon
2-3用离散相似法仿真程序sp4_3.m重求上题输出
的数据与曲线,并与四阶龙格-库塔法比较精度,同时分析两种方法对仿真步长选取的区别。
解:
程序
clearall;closeall;
a=[101011011]';
b=[0.010.0850.010.0510.00670.1510.010.01]';
c=[1111700.211300.10.0044]';
d=[00.1700.1500000]';
P=[a,b,c,d];
WIJ=[101;211;29-1;321;431;
48-1;541;651;67-0.212;761;
861;971];
n=9;Y0=1;
h=0.001;L1=10;T0=0;Tf=3;
nout=7;
A=P(:
1);B=P(:
2);
C=P(:
3);D=P(:
4);
m=length(WIJ(:
1));
W0=zeros(n,1);W=zeros(n,n);
fork=1:
m
if(WIJ(k,2)==0)
W0(WIJ(k,1))=WIJ(k,3);
else
W(WIJ(k,1),WIJ(k,2))=WIJ(k,3);
end
end
fori=1:
n
if(A(i)==0);
FI(i)=1;
FIM(i)=h*C(i)/B(i);
FIJ(i)=h*h*C(i)/B(i)/2;
FIC(i)=1;FID(i)=0;
if(D(i)~=0)
FID(i)=D(i)/B(i);
else
end
else
FI(i)=exp(-h*A(i)/B(i));
FIM(i)=(1-FI(i))*C(i)/A(i);
FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i);
FIC(i)=1;FID(i)=0;
if(D(i)~=0)
FIM(i)=(1-FI(i))*D(i)/A(i)
FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)/A(i)
FIC(i)=C(i)/D(i)-A(i)/B(i);
FID(i)=D(i)/B(i);
else
end
end
end
Y=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;
t=T0:
h*L1:
Tf;N=length(t);
fork=1:
N-1
forl=1:
L1
Ub=Uk;
Uk=W*Y+W0*Y0;
Udot=(Uk-Ub)/h;
Uf=2*Uk-Ub;
X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;
Y=FIC'.*X+FID'.*Uf;
end
y=[y,Y(nout)];
end
[t',y']
plot(t,y)
2-4求下图非线性系统的输出响应y(t),并与无非线性环节情况进行比较。
P=[0.110.51;01200;2110;10110];
WIJ=[101;14-1;211;321;431];
Z=[0000];
S=[0000];
h=0.01;
L1=25;
n=4;
T0=0;
Tf=20;
nout=4;
Y0=10;
sp4_4;
plot(t,y,'r')
holdon
Z=[4000];
S=[5000];
sp4_4;
plot(t,y,'g')
A=diag(P(:
1));B=diag(P(:
2));
C=diag(P(:
3));D=diag(P(:
4));
m=length(WIJ(:
1));
W0=zeros(n,1);W=zeros(n.n);
fork=1:
m
if(WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);
elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);
end;
end;
fori=1:
n
if(A(i,i)==0);
FI(i,i)=1;
FIM(i,i)=h*C(i,i)/B(i,i);
FIJ(i,i)=h*h*C(i,i)/B(i,i)/2;
FIC(i,i)=1;FID(i,i)=0;
if(D(i,i)~=0);
FID(i,i)=D(i,i)/B(i,i);
else
end
else
FI(i,i)=exp(-h*A(i,i)/B(i,i));
FIM(i,i)=(1-FI(i,i))*C(i,i)/A(i,i);
FIJ(i,i)=h*C(i,i)/A(i,i)-FIM(i,i)*B(i,i)/A(i,i);
FIC(i,i)=1;FID(i,i)=0;
if(D(i,i)~=0);
FIC(i,i)=C(i,i)/D(i,i)-A(i,i)/B(i,i);
FID(i,i)=D(i,i)/B(i,i);
else
end
end
end
Y=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;
t=T0:
h*L1:
Tf;N=length(t);
fork=1:
N-1
fori=1:
L1
Ub=Uk;
Uk=W*Y+W0*Y0;
fori=1:
n
if(Z(i)~=0)
if(Z(i)==1)
Uk(i,i)=satu(Uk(i,i),S(i));
end
if(Z(i)==2)
Uk(i,i)=dead(Uk(i,i),S(i));
end
if(Z(i)==3)
[Uk(i,i),Ub
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 控制系统 计算机 仿真 作业