考试用Word格式文档下载.docx

1、 y1（i+1）=y1（i）+h*（2*h+y1（i）; t1=t1,t1（i）+h;end 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; y3（i+1）=2*exp（t3（i）-t3（i）-1; t3=t3,t3（i）+h;plot（t1,y1,r,t2,y2,g,t3,y3,b系统响应曲线如下：图1 h=0.1时系统响应结果图2 h=0.01时系统响

2、应结果分析：红线-欧拉法，绿线-RK4法，蓝线-解析解。通过图中的结果我们可以看出RK4法的解与解析解更接近，其原因是欧拉法是用一阶微分方程计算得到的。再通过两个图1和图2的比较可以知道步长越短结果越靠近。2、a=0 1 0; 0 0 1; -22.06 -27 -10;b=0;0;1;c=40.6 0 0;X1=0;0;Y1=0;u=1;Y2=0;Y3=0;X2=0;t0=0;t3=0;N=round（tf-t0）/h;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

3、+2*k3+k4）/6; Y1=Y1,c*X1; x=X2（:,i）+h*（a*X2（:,i）+b*u）; y=c*x; X2=X2,x; Y2=Y2,y; 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）plot（t1,Y1,t2,Y2,t3,Y3,） 图3 h=0.1时系统的响应曲线由图可得出用RK4法得到的结果与解析方程得到的结果更接近。图4 h=0.05时系统的响应曲线步长越小，结果越精确。图6 h=0.25时系统的响应曲线：大步长越小，RK4法得到的结果与解析解更接近，但是

4、使用欧拉法得到的结果与解析方程得到的结果仍然差距很大；加大步长时，得到的结果不稳定，不能够很好的对系统进行仿真，另外，由于系统步长选择偏大，根据解析解得到的结果也与实际值有了一定的差距。3、k=1;a=conv（1 0 0,conv（0.25 1,0.25 1）b=2*k kX0=0 0 0 0 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

5、=length（b）; C=fliplr（b）,zeros（1,n-m1）; Ab=A-B*C*V; X=X0y=0;t=T0; N=round（Tf-T0）/h; 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;t,yplot（t,y）当K=1，V=1时，系统响应曲线如下：当K=2，V=1时，系统响应曲线如下：当K=1，V=2时，系统响应曲线如下：当k取值增大，v值不变时，或者当v值增大，k值不

6、变时，系统输出的波头增多，而且也变陡，稳态精度降低，当k（v）增加到一定程度时系统便发散了（即不稳定了）。4、A=-21 19 -20;19 -21 20;40 -40 -40;x=1;-1;X=x;t=0;h=0.01; 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; y1=X（1,:）;y2=X（2,:y3=X（3,: plot（t,y1,t,y2,t,y3,当h=0.01时，系统响应曲线如下：当h=0.04时，系统响应曲线如下：当h=0.06时，系统响应曲

7、线如下：如图，当h=0.01、0.04时，在t=0.2s以后系统输出便趋于平稳，当取h=0.06时，系统输出呈发散振荡形式，原先稳定的系统变得不稳定了，这便是病态系统。四思考题：1.不对，阶次越高计算机计算速度越慢，而且每个阶次都会有一个误差，那么系统的稳定性将受到影响。2.（1）精度 1）截断误差：由算法本身的精度阶次所决定 2）舍入误差：由步长决定 3）累计误差：由以上两项误差随时间积累决定（2）计算速度（3）稳定性实验二第一部分面向结构图的数值积分法仿真1、用面向方框图的数字仿真方法对下列系统进行仿真。2、求解下图所示系统在f=-1（t）阶跃扰动作用下第、第环节的动态过程。分别用面向框图

8、的数值积分法（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,1n=5;Y0=1;Yt0=0 0 0 0 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（:W0=zeros（n,1）;W=

9、zeros（n,n）;for k=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;Qn=inv（Q）;R=C*W-A;V1=C*W0;Ab=Qn*R;b1=Qn*V1;Y=Yt0y=Y（nout）;N=round（Tf-T0）/（h*L1）;N; for j=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

10、+2*k3+k4）/6; y=y,Y（nout）; t=t,t（i）+L1*h; plot（t,y）2.P=1 0.0149 1 0; 1 0.00254 26.6667 0; 1 0.391 0 0.4199; 0 0.248 1 0; 1 2.88 1 0;WIJ=1 5 1; 2 1 1; 3 4 1; 4 2 1; 4 3 -1; 5 0 1; 5 4 -1;Y0=-1;h=0.005;L1=10;Tf=10;nout=5; if （WIJ（k,2）=0）;W0（WIJ（k,1）=WIJ（k,3）; K1=Ab*Y+b1*Y0; K2=Ab*（Y+h*K1/2）+b1*Y0; K3=A

11、b*（Y+h*K2/2）+b1*Y0; K4=Ab*（Y+h*K3）+b1*Y0; Y=Y+h*（K1+2*K2+2*K3+K4）/6;y=y,Y（nout）;t=t,t（i）+h*L1;grid on;xlabel（t/sylabel（Amplitudeaxis（0 10 -0.08 0.02）;结果：第二部分 面向结构图的离散相似法仿真1、已知控制系统结构图如图所示，设输入阶跃函数幅值Y0=10，滞环非线性参数s=1（滞环宽度），请用离散相似法编程和Simulink法对系统进行如下分析：1）不考虑非线性环节影响时，求解y（t）的阶跃响应；2）考虑非线性环节影响，其余参数不变，求解y（t）并

12、与线性情况所得结果进行比较；3）改变的滞环非线性参数s，分析该非线性对系统的影响。2、系统结构图如图所示，先理论分析该系统是否会产生自激振荡，若会求出振荡的振幅和频率，并用离散相似法编程和Simulink法验证分析的结果。（1）Simulink法（2）MATLAB法P=1 10 5 10;1 0.5 1 0;1 0.1 1 0;0 1 1 0;WIJ=1 0 1;2 1 1;3 2 1;4 3 1;1 4 -1;X0=0 0 0 0;Z=0 0 0 6;S=0 0 0 1; h=0.01;L1=25;n=4;Y0=10;A=（P（:B=（P（:C=（P（:D=（P（: end if（A（i）=

13、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 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）; 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）;Y=zeros（n,1）;

14、X=Y;Uk=zeros（n,1）;Ubb=Uk;t=T0:h*L1:Tf;N=length（t）;N-1 for l=1:L1 Ub=Uk; Uk=W*Y+W0*Y0; for i=1: if（Z（i）=0） if（Z（i）=1） Uk（i）=satu（Uk（i）,S（i）; if（Z（i）=2） Uk（i）=dead（Uk（i）,S（i）; if（Z（i）=3） Uk（i）,Ubb（i）=backlash（Ubb（i）,Uk（i）,Ub（i）,S（i）; end Udot=（Uk-Ub）/h; Uf=2*Uk-Ub; X=FI.*X+FIM.*Uk+FIJ.*Udot; Yb=Y; Y=F

15、IC.*X+FID.*Uf; if（Z（i）=4） Y（i）=satu（Y（i）,S（i）; if（Z（i）=5） Y（i）=dead（Y（i）,S（i）; if（Z（i）=6） Y（i）,Ubb（i）=backlash（Ubb（i）,Y（i）,Yb（i）,S（i）;S4=1S4=5 （1）Simulink法：（2）Matlab法P=0,1,1,0;1,1,1,0;2,1,1,0;1 3 -1;3 2 1;n=3;Yt0=0 0 0 0;nout=3;Z=0 0 6;S=0 0 1;sp3_3sp3_4实验三采样控制系统的数字仿真1、已知采样系统结构如图所示，分别用程序设计方法和Simulin

16、k法求系统的输出响应。（取仿真时间: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）比较两种控制器作用下系统的响应，由此得出微分调节的作用。三，程序及结果（2）matlab法：clc;clear;G=2.72 -1;F

17、=0.717;P=0 1 1 0; 1 1 1 0; 2 1 1;n=2;X0=0 0;Ts=0.2;nout=2;A=P（:,1）;B=P（:,2）;C=P（:,3）;D=P（:,4）; if（WIJ（k,2）=0）;end if（A（i）=0） FIM（i）=h*（C（i）/B（i）; if（D（i）=0） FI（i）=exp（-h）*A（i）/B（i）; FIM（i）=（1-FI（i）*（C（i）/A（i）; FIM（i）=（1-FI（i）*（D（i）/A（i）; FIJ（i）=h*D（i）/A（i）-FIM（i）*B（i）*A（i）;Ub=Uk;U=0;uk=0;ek=0;E=zeros（n,1）;x2=0;h:t0=T0:Ts;N=length（t0）;t1=T0:Ts:M=length（t1）;M-1; U=uk; ek=Y0-x2; E=ek;E（1）; uk=-F*U+G*E; Uk=W*Y+W0*uk; y（N-1）*（k-1）+l）,1）=Y（nout）; x2=Y（nout）;plot（1:1000,y）;（1）G=15.2 -15;F=-1;P=1 1 1 0; 1 4 1 0;Ts=0.1;E=zeros（2,1）;t=0:0.01:9.99;plo