仿真方法的介绍和对比

现代工程控制理论实验报告学生姓名： 任课老师：学 号： 班 级：实验六：仿真方法的介绍与对比、实验内容及目的实验性：选用五种仿真方法对以下系统进行仿真，比较仿真精 度和仿真时间的不同其中各项系数如下kp1kiikp2ki2kiTinik2T2n20.320.0015200.00080.9373.322.08696.14前四种仿真方法为：欧拉法、局部离散法、梯形法、龙格一库塔法，之后利用整体离散法作为衡量标准实验目的：1.分析比较各种仿真方法的特点以便将来能够选择最符合要求的仿真手段2.练习使用多种仿真方法，提高仿真程序的书写能力二、实验内容及仿真结果方法一：欧拉法（1）仿真原理及程序：°1首先将系统方框图化成积分图，求出系统的空间状态方程积分图如下：i=T4xr-4-T从积分图中可以得到系统的状态空间方程x = Ax + BR，其中A=：0 000000 -kil;kl2 0 0 -ki2 000 -kpl*kl2;kp2*kl/Tl kl/11 ■ 1/T] -kp2*kl/Tl 000 -kp2*kpl*kl/ll;0 0 1/TL ■ L/TL 0000;DOO L2/T2 -1/T2 0 0 Q;0000 1/T2 -1；T2 0 0.0 0 0 0 0 L/T2 -L/r2 0;000000 1/T2 -L/T2J;B=：kil; kpl*ki2; kpl+kp2*kl/TL; 0; QP 0; Qr O'.之后利用欧拉公式：x 一 xX — —k+1 k. dt，得 x = x + X * dt。

x = 土二 dt气 1 = (I + dt * A)气 + BdtR将状态空间方程代入可得：、广(I + dt * A)xk + BdtR通过循环、迭代可以得到每次的Xk+1，对应的程序如下：R=L, dt-0. 5; st=3000; Lp=flxCst/dt): x=~erDS ^S, 1); f「： i

局部离散相似法仿真程序的主要语句如下：ta=clock;fal3=esp(-dt/Tl);falw3=L-fai3:fai5=Eip(-dt/T2);raijn5=L-rai-5;for ' -1:Ipe=R-x (S);工⑴=i Cl)+dt4kil*a;ul = i(l) + kpl*^;E ⑵=E ⑵ +dt*kl2*(Lll-E(4) >;u2=x(2)+(ul ~i (4))fkp2；e ⑶-fai3+i (3) +fai[n3+kL+Li2;k⑷=fa!3*K14J +faLffl3+K(3);z (5) =fai (5) +f 3im5+lc2tx (i);i (6:2)=fai5+i(6:B]+fain5+x(5:7);t2=[t2 dt*l];y_J.ubu=[y_jLtou xOJ];Endtinuel-c time (clock i to);(2)系统的仿真曲线:方法三：梯形法(1)仿真原理及部分语句:首先观察欧拉法的主要公式x(k +1) = x(k) + x(k)dt，其中的误差主 要来源在于此时的x(k)，误差在于将上一时刻的变化率当作了这一 时刻的增长率，称此时的误差为截断误差。

如左图，s2的大小即为截断误差的大小设E(k)表示为从k到k+1时刻的平均增长率则提高仿真精度x(k) |的实质，即为提高E(k)的精度，使其 更接近实际值而梯形法的实质便是：令 E(k)= x(k) +；(k + D，则 x(k +1) = x(k) + T^^kl^^klH而x(k) = Ax(k) + Bu(k)，可以求得准确值但x (k +1)无法得到准确值，只能借助欧拉公式求其估计值: x (k +1) = x (k) + Tx(k)x (k +1) = Ax( k +1) + Bu(k +1)而整理以上几条公式，可以得到直接的递推关系式:…一, T 2 — T T 2 T 』x (k +1) = (I + TA + — A2) x (k) + (— + — A) Bu (k) + _ Bu (k +1)2 2 2 2与整体离散法相比较:x (k +1) = M (k) +u (k)8=eAT=I+AT+—T2A2 +... 2!8 =T eATBdt= (T+-T2A1 + - T3A2...+—TkAk)B m 2! 3! k!0可以发现梯形法实际为离散相似法中的8、8”取前三项和前两项。

对应的主要语句如下：R L dt= 1. st= 3000; lp=fli(5t7dt': 3erosW 1); ._T=[Q]. 4二 L^ye If): tc^clock: El=-：ercs 0 1 :E3=zeros (fl, 1 : for祯=r+dt*El：E2=A+sa4E+E;~=(El+E2j/2:s=x4-dt*E:心=“_T :endt imel =et 1 ne i.c 1 oak. to):Ng_ fal., m tp击,teit4j= value (y_rB dt);text4=r 梯£ 法 1&真时间time 1= ' ? nuui2^tr (timsl), char (13,10)!, text4](2)系统的仿真曲线：方法四：龙格一库塔法(1)仿真原理及部分语句:龙格一库塔法实际叫做四阶龙格一库塔法为了提高公式X (k +1) = x( k) + E( k )dt中E(k)的精度，这个方法通过四个不同的变化率(E)加权平均作为 最终的E(k)设k时刻系统的状态变量x的一阶导数为E1，则由状态空间方程 可知：E1 = Ax(k) + Bu(k )之后为了寻找k +1时刻系统的变化率，需要用到k +1时刻的状态2 2变量x(k + !)。

但x(k +1)在此时是未知的，需要预测这个方法的关 22键就在于x(k + ；)的求取既然k时刻x的变化率E1已知，可以通过类似欧拉法的方式,1 , Tx (k + —) = x(k) + — E1，1 2 2将x(k + i)求出之后借助x (k + i)可以求得k +1时刻系统的变化率2 1 2 2E211E1(k + 万)=Ax1(k + % + Bu而E2已知的同时又可以通过1Tx (k + -) = x(k) + — E22 2 2求得另一个x( k+!)，由此又可以得到k+1时刻大小不同的变化率 22E3最后借助E3可以求出k+1时刻系统的状态变量x(k+1)，x(k +1) = x(k) + T * E 3E 4 = Ax (k +1) + Bu(k +1)这样又可以得到k+1时刻系统变化率的估计值E4(k+1)这样集齐了 4个变化率E1、E2、E3、E4,为了求得从k到k(k +1) 时系统的平均变化率的近似值，需要将这四个变化率求平均而E2、 E3为中间时刻求到的值，更为逼近k到k(k + 2)时系统的平均变化率， 因此权重要大一点所以：厂 E1 + 2 E 2 + 2 E 3 + E 4E = 通过x(k +1) = x(k) + E(k)dt即可求出状态变量的递推关系式主要语句如下：R=：:dt 1;st=3000；(W, L);}_lg -：-：此=皿；I (N);tc-clock;f;r 1 1：/El=Ati+B+R;苫倾1=砂<1172*巳1;E2 A*sciL+B*R;ia2=i+dt/2+E2;E3=A^Ci2J-B+R；sci3=5+

eAT=I+AT+ 2t2A2 +...8 =T eATBdt= (T+- T2A1 + - T3A2...+-TkAk )B m 2! 3! k!0为了保证精度，可以取k=20对应的matlab中的程序如下nn=20;i.ni=7erci5 N, N):for n=I : nr.fal[U=raLm-i-dT'' /facnarlal (n)应(nT」；-enclfal=I+A* Cfal> i ;faLm=f接下来利用循环仿真，可以得到系统最终的输出far 1=1:lp(N.i 1 ；t3=Lt3 dtsiZ;endt imel -etime (clock, to);[ts,赎,貌1,切,如,跋，testBZ = value (y_whole, dt);tejit3=[,整体离散法仿真时间timel= 1 , numSstr (timel), char (13, 10)F, test3](2)仿真曲线输出三、不同仿真方法之间的比较1、五种仿真方法求得的系统的输出曲线作对比1)比较不同仿真方法的仿真精度(为了便于观察曲线之间的不同，在保证能稳定的同时将仿真步长dt尽量放大)从图像中还是难以看出曲线之间的差异，以整体离散法做标准，将不同种仿真方法得到的输出与整体离散法输出的 差值的平方作为误差值。

用不同仿真方法求得的仿真误差如下：仿真方法欧拉法局部离散法梯形法龙阶一库塔整体离散法误差0.01660.49802.6034e-069.6107e。