考试用.docx
《考试用.docx》由会员分享，可在线阅读，更多相关《考试用.docx（29页珍藏版）》请在冰豆网上搜索。
考试用
实验一面向微分方程的数值积分法仿真
二、实验内容
1、已知系统微分方程及初值条件
取步长
，试分别用欧拉方程法和RK4法求
时的
值，并将求得的值与解析解
比较（将三个解绘于同一坐标中，且用数值进行比较），说明造成差异的原因。
（①编程完成；②选用MATLABode函数完成。
）
2、已知系统的传递函数为
在单位阶跃输入下，系统响应的解析解为
试分别用欧拉方程法和RK4法对系统进行仿真（编程完成）：
1）比较两种数值积分解与解析解得逼近程度；（绘图）
2）改变步长，分析步长对数值解精度的影响；
3）不断加大步长，分析计算稳定性的变化。
3、求下图所示系统的阶跃响应
的数值解。
（，
，
，
，
）分析
、
对系统响应的影响。
（①编程用RK4求解；②Simulink）
4、已知系统状态状态方程为
采用RK4法，步长分别取
，求系统的零输入响应，并绘图分析各状态变量的响应状态及产生的原因。
（提示：
病态系统）
三．实验程序及结果：
1.t0=0;
tf=2;
h=0.1;
y1=1;
y2=1;
y3=1;
t1=0;
t2=0;
t3=0
n=round（tf-t0）/h;
fori=1:
n
y1（i+1）=y1（i）+h*（2*h+y1（i））;
t1=[t1,t1（i）+h];
end
fori=1:
n
k1=y2（i）+t2（i）;
k2=y2（i）+h*k1/2+t2（i）+h/2;
k3=y2（i）+h*k2/2+t2（i）+h/2;
k4=y2（i）+h*k3+t2（i）+h;
y2（i+1）=y2（i）+h*（k1+2*k2+2*k3+k4）/6;
t2=[t2,t2（i）+h];
end
fori=1:
n
y3（i+1）=2*exp（t3（i））-t3（i）-1;
t3=[t3,t3（i）+h];
end
plot（t1,y1,'r',t2,y2,'g',t3,y3,'b'）
系统响应曲线如下：
图1h=0.1时系统响应结果
图2h=0.01时系统响应结果
分析：
红线-欧拉法，绿线-RK4法，蓝线-解析解。
通过图中的结果我们可以看出RK4法的解与解析解更接近，其原因是欧拉法是用一阶微分方程计算得到的。
再通过两个图1和图2的比较可以知道步长越短结果越靠近。
2、
a=[010;001;-22.06-27-10];
b=[0;0;1];
c=[40.600];
X1=[0;0;0];Y1=0;
u=1;
Y2=0;Y3=0;
X2=[0;0;0];
h=0.1;
t0=0;
tf=2;
t1=0;t2=0;t3=0;
N=round（tf-t0）/h;
fori=1:
N
k1=a*X1+b;
k2=b+a*（h*k1/2+X1）;
k3=b+a*（h*k2/2+X1）;
k4=b+a*（h*k3+X1）;
X1=X1+h*（k1+2*k2+2*k3+k4）/6;
Y1=[Y1,c*X1];
t1=[t1,t1（i）+h];
end
fori=1:
N
x=X2（:
i）+h*（a*X2（:
i）+b*u）;
y=c*x;
X2=[X2,x];
Y2=[Y2,y];
t2=[t2,t2（i）+h];
end
fori=1:
N
Y3（i+1）=1.84-4.95*t3（i）*exp（-1.88*t3（i））-1.5*exp（-1.88*t3（i））-0.34*exp（-6.24*t3（i））
t3=[t3,t3（i）+h];
end
plot（t1,Y1,'r',t2,Y2,'g',t3,Y3,'b'）
图3h=0.1时系统的响应曲线
分析：
由图可得出用RK4法得到的结果与解析方程得到的结果更接近。
图4h=0.05时系统的响应曲线
分析：
步长越小，结果越精确。
图6h=0.25时系统的响应曲线：
大步长越小，RK4法得到的结果与解析解更接近，但是使用欧拉法得到的结果与解析方程得到的结果仍然差距很大；加大步长时，得到的结果不稳定，不能够很好的对系统进行仿真，另外，由于系统步长选择偏大，根据解析解得到的结果也与实际值有了一定的差距。
3、
k=1;
a=conv（[100],conv（[0.251],[0.251]））
b=[2*kk]
X0=[0000]
v=1;n0=4;tf=10;h0=0.25;r=1;
V=v;
n=n0;
T0=0;
Tf=tf;
h=h0;
R=r;
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+2*k2+2*k3+k4）/6;
y=[y,C*X];
t=[t,t（i）+h];
end
[t',y']
plot（t,y）
当K=1，V=1时，系统响应曲线如下：
，
当K=2，V=1时，系统响应曲线如下：
当K=1，V=2时，系统响应曲线如下：
分析：
当k取值增大，v值不变时，或者当v值增大，k值不变时，系统输出的波头增多，而且也变陡，稳态精度降低，当k（v）增加到一定程度时系统便发散了（即不稳定了）。
4、
A=[-2119-20;19-2120;40-40-40];
x=[1;0;-1];X=x;t=0;
t0=0;tf=2;h=0.01;
n=round（tf-t0）/h;
fori=1:
n
k1=A*x;
k2=A*（x+h*k1/2）;
k3=A*（x+h*k2/2）;
k4=A*（x+h*k3）;
x=x+h*（k1+2*k2+2*k3+k4）/6;
X=[X,x];
t=[t,t（i）+h];
end
y1=X（1,:
）;y2=X（2,:
）;y3=X（3,:
）;
plot（t,y1,'r',t,y2,'g',t,y3,'b'）
当h=0.01时，系统响应曲线如下：
当h=0.04时，系统响应曲线如下：
当h=0.06时，系统响应曲线如下：
分析：
如图，当h=0.01、0.04时，在t=0.2s以后系统输出便趋于平稳，当取h=0.06时，系统输出呈发散振荡形式，原先稳定的系统变得不稳定了，这便是病态系统。
四．思考题：
1.不对，阶次越高计算机计算速度越慢，而且每个阶次都会有一个误差，那么系统的稳定性将受到影响。
2.
（1）精度
1）截断误差：
由算法本身的精度阶次所决定
2）舍入误差：
由步长决定
3）累计误差：
由以上两项误差随时间积累决定
（2）计算速度
（3）稳定性
实验二
第一部分 面向结构图的数值积分法仿真
二、实验内容
1、用面向方框图的数字仿真方法对下列系统进行仿真。
2、求解下图所示系统在f=-1（t）阶跃扰动作用下第④、第⑤环节的动态过程。
分别用面向框图的数值积分法（RK4法）、MATLAB中有关系统建模的命令和Simulink三种方法求解。
三．实验程序及结果
1.
P=[0,0.07,1,0.14;1,0.012,1,0;0,0.05,1,0.15;10,1,1,0;1,0.01,0,0.0008];
WIJ=[1,0,1;1,4,-1;2,1,1;3,2,1;3,5,-1;4,3,1;5,4,1]
n=5;Y0=1;
Yt0=[00000];
h=0.01;t=0;L1=5;
T0=0;Tf=2;nout=4;
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*1;
k2=Ab*（Y+h*k1/2）+b1*1;
k3=Ab*（Y+h*k2/2）+b1*1;
k4=Ab*（Y+h*k3）+b1*1;
Y=Y+h*（k1+2*k2+2*k3+k4）/6;
end;
y=[y,Y（nout）];
t=[t,t（i）+L1*h];
end;
plot（t,y）
2.
P=[10.014910;
10.0025426.66670;
10.39100.4199;
00.24810;
12.8810];
WIJ=[151;
211;
341;
421;
43-1;
501;
54-1;];
n=5;
Y0=-1;
Yt0=[00000];
h=0.005;
L1=10;
T0=0;
Tf=10;
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）
gridon;
xlabel（'t/s'）;
ylabel（'Amplitude'）;
axis（[010-0.080.02]）;
结果：
第二部分 面向结构图的离散相似法仿真
二、实验内容
1、已知控制系统结构图如图所示，设输入阶跃函数幅值Y0=10，滞环非线性参数s=1（滞环宽度），请用离散相似法编程和Simulink法对系统进行如下分析：
1）不考虑非线性环节影响时，求解y（t）的阶跃响应；
2）考虑非线性环节影响，其余参数不变，求解y（t）并与线性情况所得结果进行比较；
3）改变的滞环非线性参数s，分析该非线性对系统的影响。
2、系统结构图如图所示，先理论分析该系统是否会产生自激振荡，若会求出振荡的振幅和频率，并用离散相似法编程和Simulink法验证分析的结果。
三．实验程序及结果
1.
（1）Simulink法
结果：
（2）MATLAB法
P=[110510;10.510;10.110;0110];
WIJ=[101;211;321;431;14-1];
X0=[0000];
Z=[0006];S=[0001];
h=0.01;L1=25;
n=4;
T0=0;
Tf=10;
nout=4;
Y0=10;
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）;
elseW（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=（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）;Ubb=Uk;
t=T0:
h*L1:
Tf;N=length（t）;
fork=1:
N-1
forl=1:
L1
Ub=Uk;
Uk=W*Y+W0*Y0;
fori=1:
n
if（Z（i）~=0）
if（Z（i）==1）
Uk（i）=satu（Uk（i）,S（i））;
end
if（Z（i）==2）
Uk（i）=dead（Uk（i）,S（i））;
end
if（Z（i）==3）
[Uk（i）,Ubb（i）]=backlash（Ubb（i）,Uk（i）,Ub（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）=satu（Y（i）,S（i））;
end
if（Z（i）==5）
Y（i）=dead（Y（i）,S（i））;
end
if（Z（i）==6）
[Y（i）,Ubb（i）]=backlash（Ubb（i）,Y（i）,Yb（i）,S（i））;
end
end
end
end
y=[y,Y（nout）];
end
plot（t,y）
结果：
S4=1
S4=5
2.
（1）Simulink法：
结果：
（2）Matlab法
P=[0,1,1,0;1,1,1,0;2,1,1,0];
WIJ=[101;13-1;211;321];
n=3;Y0=10;
X0=[0000];
Yt0=[0000];
h=0.01;
L1=25;
t0=0;tf=10;
nout=3;Z=[006];S=[001];
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）;
elseW（WIJ（k,1）,WIJ（k,2））=WIJ（k,3）;
end;
end;
sp3_3
sp3_4
plot（t,y）
实验三采样控制系统的数字仿真
二、实验内容
1、已知采样系统结构如图所示，分别用程序设计方法和Simulink法求系统的输出响应。
（取仿真时间:
Tf=10；采样周期:
T=0.2；计算步长：
h=0.01）（注：
第一个框图中，分子分母中z的指数均为-1；第二个框图中e的指数为-Ts）
2、设某数字控制系统如图所示，采样周期为T=0.1s，初始状态
。
（注：
第二个框图中e的指数为-Ts）
（1）数字控制器为PI调节器 ，其中
（2）数字控制器为PID调节器
，其中
要求：
1）用程序设计法、MATLAB控制工具箱时域响应分析函数、Simulink法求系统在单位阶跃输入信号
作用下的输出响应；
2）比较两种控制器作用下系统的响应，由此得出微分调节的作用。
三，程序及结果
1.
（1）Simulink法
（2）matlab法：
clc;
clear;
G=[2.72-1];
F=[0.717];
P=[0110;
1110];
WIJ=[101;
211];
n=2;
Y0=1;
X0=[00];
h=0.01;
T0=0;
Tf=10;
Ts=0.2;
nout=2;
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）;
elseW（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;
U=0;uk=0;ek=0;E=zeros（n,1）;x2=0;
t=T0:
h:
Tf;
t0=T0:
h:
Ts;N=length（t0）;
t1=T0:
Ts:
Tf;M=length（t1）;
fork=1:
M-1;
U=uk;
ek=Y0-x2;
E=[ek;E
（1）];
uk=-F*U+G*E;
forl=1:
N-1
Ub=Uk;
Uk=W*Y+W0*uk;
Udot=（Uk-Ub）/h;
Uf=2*Uk-Ub;
X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;
Yb=Y;
Y=FIC'.*X+FID'.*Uf;
y（（（N-1）*（k-1）+l）,1）=Y（nout）;
end
x2=Y（nout）;
end
plot（1:
1000,y）;
结果：
2.
（1）
G=[15.2-15];
F=[-1];
P=[1110;
1410];
WIJ=[101;
211];
n=2;
Y0=1;
X0=[00];
h=0.01;
L1=25;
T0=0;
Tf=10;
Ts=0.1;
nout=2;
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）;
elseW（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;
U=0;uk=0;ek=0;E=zeros（2,1）;x2=0;
t=T0:
h:
Tf;
t0=T0:
h:
Ts;N=length（t0）;
t1=T0:
Ts:
Tf;M=length（t1）;
fork=1:
M-1;
U=uk;
ek=Y0-x2;
E=[ek;E
（1）];
uk=-F*U+G*E;
forl=1:
N-1
Ub=Uk;
Uk=W*Y+W0*uk;
Udot=（Uk-Ub）/h;
Uf=2*Uk-Ub;
X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;
Yb=Y;
Y=FIC'.*X+FID'.*Uf;
y（（（N-1）*（k-1）+l）,1）=Y（nout）;
end
x2=Y（nout）;
end
t=0:
0.01:
9.99;
plo
升级会员 