控制系统计算机仿真上机报告.docx
《控制系统计算机仿真上机报告.docx》由会员分享，可在线阅读，更多相关《控制系统计算机仿真上机报告.docx（20页珍藏版）》请在冰豆网上搜索。
控制系统计算机仿真上机报告
******理工大学
《控制系统计算机仿真》上机报告
院系：
班级：
姓名：
学号：
时间：
年月日
2-2用MATLAB语言求下列系统的状态方程、传递函数、零极点增益和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：
（1）
（2）
解：
（1）源程序：
num=[172424];
den=[110355024];
G=tf（num,den）
[Z,P,K]=tf2zp（num,den）
[A,B,C,D]=tf2ss（num,den）
[R,P,H]=residue（num,den）
结果：
Transferfunction:
s^3+7s^2+24s+24
---------------------------------
s^4+10s^3+35s^2+50s+24
Z=
-2.7306+2.8531i
-2.7306-2.8531i
-1.5388
P=
-4.0000
-3.0000
-2.0000
-1.0000
K=
1
A=
-10-35-50-24
1000
0100
0010
B=
1
0
0
0
C=
172424
D=
0
R=
4.0000
-6.0000
2.0000
1.0000
P=
-4.0000
-3.0000
-2.0000
-1.0000
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]
B=[4;2;2;0]
C=[0202]
D=[0]
[num,den]=ss2tf（A,B,C,D）
[Z,P,K]=tf2zp（num,den）
[R,P,H]=residue（num,den）
结果：
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
C=
0202
D=
0
num=
04.000014.000022.000015.0000
den=
1.00004.00006.25005.25002.2500
Z=
-1.0000+1.2247i
-1.0000-1.2247i
-1.5000
P=
-1.5000
-1.5000
-0.5000+0.8660i
-0.5000-0.8660i
K=
4.0000
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=
[]
2-3用殴拉法求下列系统的输出响应
在
上，
时的数值解。
，
要求保留4位小数，并将结果以图形的方式与真解
比较。
解：
源程序：
t=0:
0.1:
1;
h=0.1;
y
（1）=1;
fori=1:
1:
10;
y（i+1）=y（i）+h*（-2*y（i））;
Y=[y,y（i）]
end
plot（t,y,'r*'）;
holdon
m=exp（-2*t）;
plot（t,m）
结果：
Y=
Columns1through8
1.00000.80000.64000.51200.40960.32770.26210.2097
Columns9through12
0.16780.13420.10740.1342
分析：
经由图像比较，用欧拉法求得的解与真值
大体一致。
2-5用四阶龙格－库塔梯形法求解2-3的数值解，并通过与真值及殴拉法的比较，分析其精度。
解：
源程序：
t=0:
0.1:
1;
h=0.1;
y
（1）=1;
fori=1:
1:
10;
k1=-2*（y（i））;
k2=-2*（y（i）+h/2*k1）;
k3=-2*（y（i）+h/2*k2）;
k4=-2*（y（i）+h*k3）;
y（i+1）=y（i）+h/6*（k1+2*k2+2*k3+k4）;
Y=[y,y（i）]
end
plot（t,y,'r*'）;
holdon
m=exp（-2*t）;
plot（t,m）
结果：
Y=
Columns1through8
1.00000.81870.67030.54880.44930.36790.30120.2466
Columns9through12
0.20190.16530.13530.1653
分析：
经由图像比较，用四阶龙格－库塔梯形法求解题2-3的数值解与真值
一致，比用欧拉法求得的解更接近于真值。
4-2设典型闭环结构控制系统如下图所示，当阶跃输入幅值
时，用sp4_1.m求取输出y（t）的响应。
解：
源程序：
b=[1.875e+6,1.562e+6];
a=[1,54,204.2,213.8,63.5];
X0=[0,0,0,0];
V=0.002;
n=4;
T0=0;
Tf=10;h=0.01;R=20;
b=b/a
（1）;a=a/a
（1）;A=a（2:
5）;
A=[rot90（rot90（eye（3,4）））;-fliplr（A）];
B=[zeros（1,3）,1]';
m1=length（b）;
C=[fliplr（b）,zeros（1,4-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+2*K2+2*K3+K4）/6;
y=[y,C*X];
t=[t,t（i）+h];
end
[t',y']
plot（t,y）
b=[1.875e+6,1.562e+6];
a=[1,54,204.2,213.8,63.5];
X0=[0,0,0,0];
V=0.002;
n=4;
T0=0;
Tf=10;h=0.01;R=20;
b=b/a
（1）;a=a/a
（1）;A=a（2:
5）;
A=[rot90（rot90（eye（3,4）））;-fliplr（A）];
B=[zeros（1,3）,1]';
m1=length（b）;
C=[fliplr（b）,zeros（1,4-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+2*K2+2*K3+K4）/6;
y=[y,C*X];
t=[t,t（i）+h];
end
[t',y']
plot（t,y）
结果：
分析：
由结论及图形反映了开环放大洗漱K与反馈系数V对结果的影响，K值增大会使系统震荡加剧，甚至可能造成不稳定。
4-5下图中，若各环节传递函数已知为：
，，
，
，
，
，但
；试列写链接矩阵W、W0和非零元素阵WIJ，将程序sp4_2完善后，应用此程序求输出
的响应曲线。
解：
源程序：
P=[1,0.01,1,0;
0,0.085,1,0.017;
0,0.051,1,0.15;
1,0.15,0.21,0;
1,0.01,0.1,0;
1,0.01,0.0044,0;
];
WIJ=[101;
211;
26-1;
321;
35-1;
431;
44-0.212;
541;
641
];
n=6;
Y0=1;
Yt0=[000000];
h=0.01;
L1=10;
T0=0;
Tf=40;
nout=5;
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;
Q=B-D*W;Qn=inv（Q）;
R=C*W-A;V1=C*W0;
Ab=Qn*R;b1=Qn*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+2*K2+2*K3+K4）/6;
end
y=[y,Y（nout）];
t=[t,t（i）+h*L1];
end
[t',y']
plot（t,y）
结果：
4-7用离散相似法仿真程序sp4_3.m重求上题输出
的数据与曲线，并与四阶龙格－库塔法比较精度。
解：
源程序：
P=[1,0.01,1,0;0,0.085,1,0.17;0,0.051,1,0.15;1,0.15,0.21,0;1,0.01,0.1,0;1,0.01,0.0044,0];
W=[000000;10000-1;0100-10;001-0.21200;000100;000100];
W0=[1;0;0;0;0;0];
n=6;
Y0=1;
Yt0=[000000];
h=0.01;
L1=10;
T0=0;
Tf=40;
nout=4;
A=diag（P（:
1））;B=diag（P（:
2））;
C=diag（P（:
3））;D=diag（P（:
4））;
fori=1:
6
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（6,1）;X=Y;y=0;Uk=zeros（6,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
plot（t,y,'r*'）
结果：
分析：
用离散相似法与四阶龙格—库塔法比较来看，两者精度很接近，并且与四阶龙格-库塔法相比，离散相似法更容易推广到解决非线性系统的仿真问题。
4-8求下图非线性系统的输出响应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'）
子程序sp4_4.m文件：
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）,Ubb（i,i）]=backlash（Ubb（i,i）,Uk（i,i）,Ub（i,i）,S（i））;
end
end
end
Udot=（Uk-Ub）/h;
Uf=2*Uk-Ub;
X=FI'*X+FIM'*Uk+FIJ'*Udot;
Yb=Y;
Y=FIC'*X+FID'*Uf;
fori=1:
n
if（Z（i）~=0）
if（Z（i）==4）
Y（i,i）=satu（Y（i,i）,S（i））;
end
if（Z（i）==5）
Y（i,i）=dead（Y（i,i）,S（i））;
end
if（Z（i）==6）
[Y（i,i）,Ubb（i,i）]=backlash（Ubb（i,i）,Y（i,i）,Yb（i,i）,S（i））;
end
end
end
end
y=[y,Y（nout）];
end
backlash函数：
function[Uc,Ubb]=backlash（Urb,Ur,Ucb,S1）
if（Ur>Urb）
if（（Ur-S1）>=Ucb）
Uc=Ur-S1;
elseUc=Ucb;
end
elseif（Urif（（Ur+S1）<=Ucb）
Uc=Ur+S1;
elseUc=Ucb;
end
elseUc=Ucb;
end
end
Ubb=Ur;
Satu函数：
functionUc=satu（Ur,S1）
if（abs（Ur）>=S1）
if（Ur>0）
Uc=S1;
elseUc=-S1;
end
elseUc=Ur;
end
dead函数：
functionUc=dead（Ur,S1）
if（abs（Ur）>=S1）
if（Ur>0）
Uc=Ur-S1;
elseUc=Ur+S1;
end
elseUc=0;
end
结果：
分析：
仿真结果显示滞环非线性缓解N对线性系统输出响应的影响，N值越大，滞后越明显。
升级会员 