经典的连续系统仿真建模方法学

32第二章经典的连续系统仿真建模方法学本章讨论经典的连续系统数字仿真的原理与方法，内容包括连续系统数字仿真的基本概念、经典的数值积分法、经典的线性多步法等。在数字计算机上进行连续系统仿真，首先要将连续模型离散化，因此，2.1节首先讨论离散化原理及要求，这是连续系统仿真的基础。然后，2.2节对经典的数值积分法----龙格-库塔法及其它典型的数值积分法仿真建模原理进行详细分析，并通过实例说明其应用要点；而2.3节对经典的线性多步法进行了介绍．2.1离散化原理及要求在数字计算机上对连续系统进行仿真时，首先遇到的问题是如何解决数字计算机在数值及时间上的离散性与被仿真系统数值及时间上的连续性这一基本问题。从根本意义上讲，数字计算机所进行的数值计算仅仅是“数字”计算，它表示数值的精度受限于字长，这将引入舍入误差；另一方面，这种计算是按指令一步一步进行的，因而，还必须将时间离散化，这样就只能得到离散时间点上系统性能。用数字仿真的方法对微分方程的数值积分是通过某种数值计算方法来实现的。任何一种计算方法都只能是原积分的一种近似。因此，连续系统仿真，从本质上是对原连续系统从时间、数值两个方面对原系统进行离散化，并选择合适的数值计算方法来近似积分运算，由此得到的离散模型来近似原连续模型。如何保证离散模型的计算结果从原理上确能代表原系统的行为，这是连续系统数字仿真首先必须解决的问题。设系统模型为：),,(tuyfy，其中u(t)为输入变量，y(t)为系统变量；令仿真时间间隔为h，离散化后的输入变量为)(ˆktu，系统变量为)(ˆkty，其中kt表示t=kh。如果)()(ˆkktutu，)()(ˆkktyty，即0)()(ˆ)(kkkututute，0)()(ˆ)(kkkytytyte（对所有k=0,1,2,…），则可认为两模型等价，这称为相似原理（参见图2.1）。实际上，要完全保证0)(,0)(kykutete是很困难的。进一步分析离散化引入的误差，随着计算机技术的发展，由计算机字长引入的舍入误差可以忽略，关键是数值积分算法，也称为仿真建模方法。相似原理用于仿真时，对仿真建模方法有三个基本要求：（1）稳定性：若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。关于稳定性的详细讨论将在2.4节中进行。（2）准确性：有不同的准确性评价准则，最基本的准则是：原连续模型),,(tuyfy仿真模型),ˆ,ˆ(ˆktuyfyu(t)hy(t))(ˆkty)(ˆktu-+0)(kyte图2.1相似原理33f(t,y)f(t0,yo)ttt0t1图2.2数值积分法原理绝对误差准则：)()(ˆ)(kkkytytyte相对误差准则：)(ˆ)()(ˆ)(kkkkytytytyte其中规定精度的误差量。（3）快速性：如前所述，数字仿真是一步一步推进的，即由某一初始值)(0ty出发，逐步计算，得到)(,),(),(21ktytyty，每一步计算所需时间决定了仿真速度。若第k步计算对应的系统时间间隔为,1kkktth计算机由)(kty计算)(1kty需要的时间为Tk，则，若Tk=hk称为实时仿真，Tkhk称为超实时仿真，而大多数情况是Tkhk，对应于离线仿真。连续系统数字仿真中离散化最基本的算法是数值积分算法。对于形如),,(tuyfy的系统，已知系统变量y的初始条件yty()00，现在要求y随时间变化的过程yt()。计算过程可以这样考虑（参见图2.2）：首先求出初始点yty()00的fty()00，，微分方程可以写作：ytyftydttt()(,)00（2.1）图2.2所示曲线下的面积就是yt()，由于难以得到f(y,u,t)积分的数值表达式，人们对数值积分方法进行了长期探索，其中欧拉法是最经典的近似方法。欧拉法用矩形面积近似表示积分结果，也就是当t=t1时，yt()1的近似值为y1：yytytfty11000()()，(2.2)重复上述作法，当tt2时yytyttfty2212111()()()，所以，对任意时刻tk+1，有：yytyttfty111()()()，（2.3）令tth1称为第步的计算步距。若积分过程中步距不变hh，可以证明，欧拉法的截断误差正比于h2。为进一步提高计算精度，人们提出了“梯形法”。梯形法近似积分形式如式(2.4)所示，令：tthh1＝已知：tt时yt()的近似值y,那么：yytyhftyfty111112()[(,)(,)](2.4)可见，梯形法是隐函数形式。采用这种积分方法最简单的预报校正方法是用欧拉法估计初值，用梯形法校正，即：yyhftyftyii111112()()[(,)(,)](2.5)yyhftyi1()(,)(2.6)34式(2.6)称作预报公式，采用欧拉法，式(2.5)为校正公式，采用梯形法。用欧拉法估计一次yi1()的值，代入校正公式得到y1的校正值yi11()。设是规定的足够小正整数，称作允许误差，若i=0,i+1=1称作第一次校正；ii112，称作第二次校正；通过反复迭代，直到满足yyii111，这时yi11()是满足误差要求的校正值。上述方法是针对(2.3)式所示的微分方程在已知初值情况下进行求解，因此也称为微分方程初值问题数值计算法，为统一起见，本书中称为数值积分法。连续系统数字仿真的离散化方法有两类，它们是数值积分方法和离散相似方法，本章讨论数值积分法。数值积分方法采用递推方式进行运算，而采用不同的积分方法会引进不同的计算误差，为了提高计算精度，往往会增加运算量。就同一种积分算法而言，为提高计算精度，减小积分步距h,计算量增大，影响系统运算速度。因此，计算精度和速度是连续系统仿真中常迂到的一对矛盾，也是数字仿真中要求解决的问题之一。也就是说，选择合适的算法、合适的软、硬件环境，在保证计算精度的前提下，考虑怎样提高仿真的速度。经典的数值积分法可分为两类：单步法与多步法，下面我们将分别来介绍这两类方法中的最常用的算法。2.2龙格库塔法2.2.1龙格-库塔法基本原理在上一节中我们已经讲过，在连续系统的仿真中，主要的数值计算工作是对dydtfty(,)的一阶微分方程进行求解。因为ytytftydttt()()(,)11若令：yyt(),Qftydtttk(,)1则有ytyy()11Q(2.7)因此主要问题就是如何对Q进行数值求解，即如何对fty(,)进行积分。通常称作“右端函数”计算问题。已知yty()00，假设我们从t0跨出一步，tth10，t1时刻为yyth10()，可以在t0附近展开台劳级数，只保留h2项，则有：200010)(21),(htfdtdyyfhytfyyt(2.8)在式(2.8)中括号后的下标0表示括号中的函数将用ttyy00，均同。我们假设这个解可以写成如下形式：yyakakh101122()(2.9)其中，kfty100(,)kftbhybkh201021()，对k2式右端的函数在ttyy00，处展成台劳级数，保留h项，可得：35hkbtfbytfkt0)(),(121002将kk12，代入(2.9)式，则有：])(),([),(012100200101hyfkbtfbytfhaythfayyt将上式与(2.8)式进行比较，可得：可见有四个未知数aabb1212，，，，但只有三个方程，因此有无穷多个解，若限定aa12，则可得其中一个解：aabb1212121，将它们代入(2.9)式可得一组计算公式：yyh10122()(2.10)其中，1002001ftyfthyh(,)()，，若写成一般递推形式，即为：ytyyhkkk()()11122(2.11)其中，121ftyfthyhkkkk(,)(,)，由于(2.8)式只取了hh，2两项，而将h3以上的高阶项略去了，所以这种递推公式的截断误差正比于h的三次方，又由于计算时只取了h及h2项，故这种方法被称为二阶龙格库塔法（简称RK-2）。根据上述原理，若在展成台劳级数时保留hhh，，23及h4项，那么可得一套截断误差正比于h5的四阶龙格--库塔法（简称RK-4）公式：ytyykkk()()1112341622(2.12)其中：1fty(,)2122fthyh()，3222fthyh()，43fthyh()，由于这组计算公式有较高的精度，所以在数字仿真中应用较为普遍。为了帮助读者更好地掌握这种方法的使用，下面我们对龙格--库塔法的特点作一些介绍。1.由于在解aabb1212，，，时，可以得到许多种龙格--库塔公式，经常使用的除(2.11)及(2.12)式给出的两组公式外，还有：yyh12(2.13)aaabab12212211212+===,/,/36其中1fty(,)2122fthyh()，yyh11234833()(2.14)其中1fty(,)2133fthyh()，312233fthyhh()，4123fthyhhh()，(2.13)式为二阶龙格--库塔公式，(2.14)式是Shampine提出的，称为RKS3-4公式。对上述公式可以这样解释：因为yt()满足dydtfty(,)这个微分方程，已知(,)ty处yt()的导数为fty(,)1。根据外推原理，假设在t--th这个区间中yt()的导数不变，则可得yt()在tth＋的估计值yhi0。令i为这个估计值处的导数，因此(2.11)式就是按()122这个导数线性外推而得的y1。而(2.12)式则是4点导数的加权平均和，由于第2及第3点是以h2为步长，精度较高，所以加权平均时各取2份，第1、4点各取1份。不论几阶龙格--库塔法，它们的计算公式总是由两部分组成的，即是上一步的结果y及步长h乘以t--t1中各点导数的加权平均和。因此各种龙格---库塔法可以写成如下一般形式：yyhCiiis11(2.15)其中iiijjjiftahyhb()，11is12，，，式中各系数满足以下关系aabisCiijjiiis11110231，，，(2.16)一般来讲，为了减少计算量总是希望减少每步计算右端函数f的次数（即减少计算的次数）。可以证明，1阶龙格--库塔公式至少要计算一次右端函数，即(2.11)式中的s(称为级数)至少等于1，而2阶公式smin2；3阶公式smin3；4阶公式smin4；依此类推。有时为了某种特殊需要，可以选择ssmin的计算公式。表2.1为常见的各阶龙格---库塔公式，其中第二项是4阶4级RK-4公式1，而第三项是374阶5级RK-4计算公式(简称RKM45)。其中，后三项有误差估计公式，所以可以用误差来控制积分步长，即每积分一步，利用误差估计公式估计出这一步的截断误差。若它大于允许的误差，则缩小步长，反之可以放大步长。表2.1常见的龙格--库塔公式名称公式特点二阶龙格---库塔法yyhftyfthyh1121212()(,)()，四阶龙格---库塔法yyhftyfthyhfthyhfthyh1123412132436222222()(,)()()()，，，四