控制系统仿真课件2

1、l连续系统仿真的数值积分算法就是利用数值积分法将常微分方程（组）描述的连续系统变换成离散形式的仿真模型差分方程（组）。l为了能在计算机上进行求解，首先要把被仿真系统的数学模型表示为一阶微分方程组或状态空间模型。 假设一阶向量微分方程及初值问题为00)(),(xtxxtxf可以写成一阶微分方程组和初值问题形式00200210012121222111)()()(),(),(),(nnnnnnnxtxxtxxtxxxxtfxxxxtfxxxxtfx (2.1)(2.2)110,kkttttt110d),()(d),()()(01kkkttkttktxttxtxttxtxff在上的解析解为 (2.3)

2、希望用公式 kkkxxq11kxkxkq)(1ktx)(ktx1d),(kktttxtf来代替解析解，其中 ， ， 分别是 ， ， 的近似解。 (2.4) 所谓数值解法就是寻求初值问题的解在一系列离散时间点 121,kktttt,21xx1,kkxx上的近似解。相邻两个离散时间点的间距 kktth1称为步长或步距 。数值积分法的主要问题归结为如何对向量函数 ),(yf t进行数值积分，求出 ),(yf t在区间 ),(1kktt上定积分的近似解 kq ，并且数值方法的共同特点是步进式的（即所谓递推式的）。 常用的基本算法可以分为单步法、多步法和预估-校正法三大类，而每一类算法又可以分为显式算法

3、和隐式算法两种类型。.)(!.)(! 2)()()()(0)(0)2(20)1(001txkhtxhhtxtxhtxtxkk00)(),(xtxxtfx 已知0 xhtthtthttnn11201,.,nxxxx,.,321从初值开始，在一系列时刻点求解的近似解。(h是计算步长）若x对t的各阶导数都存在，)(txhtt01)(0htx在时的解用泰勒级数表示为:)(1tx)(2tx)(txnxxxx,.,321 在h很小的情况下，可以只取上式前两项、三项或四项之和作近似计算，从而得到)(1tx的计算值，再把当做初值计算，依此类推逐步计算出的数值解00)(),(xtxxtfx 欧拉（Euler）法

4、是一种最简单的数值积分算法，对于 在区间 上求积分，有 ),(1kktt1d),()()(1kkttkktxtftxtx (2.8)其取泰勒公式的前两项来近似计算11),()(kkkkkxxthfxtx即有),(kkkxthfq 欧拉法为 kkkhfxx1式中),(kkkxtff (2.9) (2.10) 欧拉法的图形表示： 为图2.1中的非曲面面积actk+1tk，欧拉法用矩形面积abtk+1tk代替。1d),(kktttxtf显然，欧拉法仅适用于步长h很小的场合。 图2.1txt0 x0t1x1t2x22h欧拉法是将泰勒公式截断 以上的高阶项而得到的，)(1phO)(2hORn),(!.)

5、,(! 3),(! 200)(00)3(300)2(2xtxkhxtxhxtxhRkhn通常把截断项称截断误差，用其表示数值计算方法的精度，截断误差为时称为p阶精度。对于欧拉法说明p=1具有1阶精度。截断误差l由数学家卡尔龙格和马丁威尔海姆库塔于1900年左右发明。l龙格-库塔（Runge-Kutta）法是求解常微分方程初值问题(2.5)式的各种数值积分算法中应用得最广泛的一种，包括许多不同的公式。l思路：用若干个时间点上 f 的函数值的线性组合来代替 f 的各阶导数项，然后按泰勒公式展开确定其中的系数。 以二阶龙格库塔法为例 考虑泰勒展开式的前三项之和作近似计算 为避免上式右边微分计算，为此

6、考虑用函数 ),( xtf在点 ),(00 xtf和点 ),(12010hkbxhbtf的线性组合 来代替（2-11）的右端后两项，即假设 ),(001xtfk ),(120102hkbxhbtfk0)2(2001! 2xhxhxx (2.11)(221101kakahxx则 (2.12)1.求式（2-11）分解)*(! 2),(00,20001xtxfftfhxthfxx2.求式（2-12）的分解0012100120102),(),(xxttxfhkbtfhbxtfhkbxhbtfk),(),(0012100200101xxttxfhkbtfhbxtfhaxtfhaxx3.比较两式得2/12

7、/11221221babaaa为此设 21bb ，得到 12/12121bbaa至此二阶龙格-库塔法的计算公式为),(),()(210020012101hkxhtfkxtfkkkhxx若写成一般递推形式，即为),(),()(2121211hkxhtfkxtfkkkhxxnnnnnn (2.19) 根据上述原理，在泰勒展开时保留h项，称之为一阶龙格-库塔公式（即欧拉法） ; 若保留h、h2、h3项，称之为三阶龙格-库塔公式，若保留h、h2、h3、h4项，称之为四阶龙格-库塔公式，余类推。)2,()2,2(),()4(62131213211hkhkxhtfkhkxhtfkxtfkkkkhxxnnn

8、nnnnn三阶龙格三阶龙格-库塔公式库塔公式四阶龙格-库塔公式（古典形式） 该算法是数字仿真中最常用的一种算法。计算量较大，但计算精度较高。在实际仿真中得到了广泛的应用。 ),()2/, 2/()2,2(),()22(6342312143211hkxhtfkhkxhtfkhkxhtfkxtfkkkkkhxxnnnnnnnnnn四阶龙格-库塔公式（吉尔公式）)221 (22,()221 () 12(2, 2/()2,2(),()22()22(632421312143211kkxhtfkkkhxhtfkhkxhtfkxtfkkkkkhxxnnnnnnnnnnl几种数值积分法都可以看成是 在 附近用

9、泰勒级数展开而产生的，只不过是泰勒级数所取项数的多少不同而已。欧拉法只取前两项， RK2法取了前三项，RK4法取了前五项。从理论上讲，可以构造任意高阶的计算方法。但是，精度的阶数与计算函数f值的次数之间的关系并非等量增加的。 )(tyktt 龙格-库塔法的几点说明1.公式采用了平均斜率的概念，提高了计算精度。2.计算公式由两部分组成；3.步长在整个计算中并不要求固定， 但在每一步的计算必须用同一步长。4.单步法。优点：自启动、占用存贮空间小5.龙格-库塔法的精度取决于步长及计算方法。 以四阶为例，达到相同的精度，四阶的步长可以比二阶大10倍，而四阶的每步计算量仅比二阶的大一倍，所以常用四阶。高

10、于四阶是否更好？h是否取得越小精度就越高？ h越小计算次数越多，引入较大的舍人误差。更高阶时精度提高不明显 ;【举例】：已知方程5 . 0, 1 . 01)0(1thxtxx计算到取本题的精确解为tettx)(方法一：采用欧拉法hthxhtxhxtxtfhxxkkkkkkkkk*)1 () 1( *)(,(*109049. 1.5 . 00561. 1.4 . 0029. 1.3 . 001. 11 . 0*1 . 0) 1 . 01 (.2 . 011 . 0*1 . 0) 1 . 01 (.1 . 01.04544331122001100 xtxtxttxxttxxtxt方法二：经典四阶R

11、-K方法求解。当t=0时，09525. 0)*1 ,(0475. 0)*2/11 ,2/1 (05. 02/11 ,2/111 , 0,342312001khhfkkhhfkhkhfkfxtfk于是00483750. 1)22(6432101kkkkhxx10653094. 107032029. 104081842. 101873090. 15432xxxx依次递推得tEuler法梯形法四阶R-K法精确解011110.11.0000001.0047621.004837501.004837410.21.0100001.0185941.01873.901.018730750.31.0290001.

12、0406331.040818421.040818220.41.0561001.0700971.070320291.070320050.51.0904901.1062781.106530941.10653066欧拉法计算程序清单：x=1;h=0.1；y=0;for t=0:0.1:0.5; f=-x+t+1 x=x+h*f y=y,xend四阶R-K方法程序清单：x=1;h=0.1;y=0;for t=0:0.1:0.5; k1=-x+t+1 k2=-(x+k1*h/2)+(t+h/2)+1 k3=-(x+k2*h/2)+(t+h/2)+1 k4=-(x+k3*h)+(t+h)+1 x=x+h/

13、6*(k1+2*k2+2*k3+k4) y=y,xend若单变量系统的状态空间表达式为yCxBuAxx其中，A：n*n;B:n*1;C:1*n.假设在仿真中，数值积分法采用四阶龙格-库塔法。因对于高阶系统，状态方程式可写成以下个一阶微分方程1 122123( ,)iiiinniinxa xa xa xbuf t x x xx求解一阶微分方程组式的四阶龙格-库塔公式如下:11234(22)6kkiiiiiihxxkkkk(1,2, )in11 122211113112241133( )()()()222()()()222()()()kkkiiiinnikkkiiiinniikkkiiiinnii

14、kkkiiiinniikka xa xa xbu thhhkaxkaxkbu thhhkaxkaxkbu tkaxhkaxhkbu th 式中, 为 时刻的 值； 为 时刻 值 kixkttix1kixktthix111112Tkkkknxxxx12Tkkkknxxxx111211TnKkkk212222TnKkkk313233TnKkkk414244TnKkkk 令则式（2-4）可写成如下矩阵的形式11234(22)6kkhxxKKKK1213243( )()()22()()22()()kkkkkkkkKAxbu thhKA xKbu thhKA xKbu tKA xhKbu th 例2 假

15、设单变量系统框图如图所示。试根据四阶龙格-库塔法，求输出量y的阶跃响应。解 仿真程序Example2_2.m如下 %Example2_2.m r=1; numo=1,2; deno=1 3 6 2; numh=3; denh=1; num,den=feedback(numo,deno,numh,denh); A,b,C,D=tf2ss(num,den); Tf=input(仿真时间Tf=); h=input(计算步长h=); x=zeros(length(A),1); y=0; t=0; for i=1:Tf/h K1=A*x+b*r; K2=A*(x+h*K1/2)+b*r; K3=A*(x

16、+h*K2/2)+b*r; K4=A*(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 plot(t,y)5fT 0.02h 取仿真时间：计算步长：在MATLAB环境下执行以上程序可得如图所示仿真曲线。l如果在欧拉法中积分值改用梯形来代替，可提如果在欧拉法中积分值改用梯形来代替，可提高精度，则高精度，则),(),(2111kkkkkkxtfxtfhxx此式称为二阶隐式亚当姆斯法公式或梯形积分公式。 这种隐式公式在计算机上通常采用迭代法来求解，即先猜一个初值 ，求出 ，然后用 再去求 )0(1kx)1(1kx)1(1kx)2(1kx，直到 与 相差甚小为止 ，此时即认为 )(1nkx)1(1nkx1)(1knkxx 这种方法计算复杂，为了简化计算，通常用显式多步法来计算初值，然后用隐式多步法只做一次校正计算。比如采用一阶显式计算初值，然后用二阶