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

第二章 经典的连续系统仿真建模方法学本章讨论经典的连续系统数字仿真的原理与方法，内容包括连续系统数字仿真的基本概念、经典的数值积分法、经典的线性多步法等。

在数字计算机上进行连续系统仿真，首先要将连续模型离散化，因此，2.1节首先讨论离散化原理及要求，这是连续系统仿真的基础。

然后，2.2节对经典的数值积分法----龙格-库塔法及其它典型的数值积分法仿真建模原理进行详细分析，并通过实例说明其应用要点；而2.3节对经典的线性多步法进行了介绍．2.1 离散化原理及要求在数字计算机上对连续系统进行仿真时，首先遇到的问题是如何解决数字计算机在数值及时间上的离散性与被仿真系统数值及时间上的连续性这一基本问题。

从根本意义上讲，数字计算机所进行的数值计算仅仅是“数字”计算，它表示数值的精度受限于字长，这将引入舍入误差；另一方面，这种计算是按指令一步一步进行的，因而，还必须将时间离散化，这样就只能得到离散时间点上系统性能。

用数字仿真的方法对微分方程的数值积分是通过某种数值计算方法来实现的。

任何一种计算方法都只能是原积分的一种近似。

因此，连续系统仿真，从本质上是对原连续系统从时间、数值两个方面对原系统进行离散化，并选择合适的数值计算方法来近似积分运算，由此得到的离散模型来近似原连续模型。

如何保证离散模型的计算结果从原理上确能代表原系统的行为，这是连续系统数字仿真首先必须解决的问题。

设系统模型为：),,(t u y f y= ，其中u (t )为输入变量，y (t )为系统变量；令仿真时间间隔为h ，离散化后的输入变量为)(ˆk t u，系统变量为)(ˆk t y ，其中k t 表示t=kh 。

如果)()(ˆk k t u t u ≈，)()(ˆk k t y t y ≈，即0)()(ˆ)(≈-=k k k u t u t u t e ，0)()(ˆ)(≈-=k k k y t y t y t e （对所有k=0,1,2,…），则可认为两模型等价，这称为相似原理（参见图2.1）。

实际上，要完全保证0)(,0)(==k y k u t e t e 是很困难的。

进一步分析离散化引入的误差，随着计算机技术的发展，由计算机字长引入的舍入误差可以忽略，关键是数值积分算法，也称为仿真建模方法。

相似原理用于仿真时，对仿真建模方法有三个基本要求：（1）稳定性：若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。

关于稳定性的详细讨论将在2.4节中进行。

（2）准确性：有不同的准确性评价准则，最基本的准则是：图2.1 相 似 原 理绝对误差准则：δ≤-=)()(ˆ)(k k k y t y t yt e 相对误差准则：δ≤-=)(ˆ)()(ˆ)(k k k k y t yt y t yt e其中δ 规定精度的误差量。

（3）快速性：如前所述，数字仿真是一步一步推进的，即由某一初始值)(0t y 出发，逐步计算，得到)(,),(),(21k t y t y t y ，每一步计算所需时间决定了仿真速度。

若第k 步计算对应的系统时间间隔为,1k k k t t h -=+计算机由)(k t y 计算)(1+k t y 需要的时间为T k ，则，若T k =h k 称为实时仿真，T k <h k 称为超实时仿真，而大多数情况是T k >h k ，对应于离线仿真。

连续系统数字仿真中离散化最基本的算法是数值积分算法。

对于形如),,(t u y f y= 的系统，已知系统变量y 的初始条件y t y ()00=，现在要求y 随时间变化的过程y t ()。

计算过程可以这样考虑（参见图2.2）：首先求出初始点y t y ()00=的f t y ()00，，微分方程可以写作：y t y f t y dt t t()(,)=+⎰00（2.1）图2.2所示曲线下的面积就是y t ()，由于难以得到f(y,u,t)积分的数值表达式，人们对数值积分方法进行了长期探索，其中欧拉法是最经典的近似方法。

欧拉法用矩形面积近似表示积分结果，也就是当t=t 1时，y t ()1的近似值为y 1 ：y y t y t f t y 11000=≅+⋅()()∆， (2.2)重复上述作法，当t t =2时y y t y t t f t y 2212111=≅+-⋅()()()， 所以，对任意时刻t k+1，有：y y t y t t f t y κκκκκκκ+++=≅+-⋅111()()()， （2.3） 令t t h κκκ+-=1称为第κ步的计算步距。

若积分过程中步距不变h h κ=，可以证明，欧拉法的截断误差正比于h 2。

为进一步提高计算精度，人们提出了“梯形法”。

梯形法近似积分形式如式(2.4)所示，令：t t h h κκκ+-=1＝已知：t t =κ时y t ()κ的近似值y κ,那么：y y t y h f t y f t y κκκκκκκ++++=≅++111112()[(,)(,)] (2.4) 可见，梯形法是隐函数形式。

采用这种积分方法最简单的预报−校正方法是用欧拉法估计初值，用梯形法校正，即：y y h f t y f t y i i κκκκκκ++++≅++111112()()[(,)(,)] (2.5) y y h f t y i κκκκ+≅+⋅1()(,) (2.6)式(2.6)称作预报公式，采用欧拉法，式(2.5)为校正公式，采用梯形法。

用欧拉法估计一次y i κ+1()的值，代入校正公式得到y κ+1的校正值y i κ++11()。

设ε是规定的足够小正整数，称作允许误差，若i =0, i +1=1称作第一次校正；i i =+=112，称作第二次校正；通过反复迭代，直到满足y y i i κκε+++-≤111，这时y i κ++11()是满足误差要求的校正值。

上述方法是针对(2.3)式所示的微分方程在已知初值情况下进行求解，因此也称为微分方程初值问题数值计算法，为统一起见，本书中称为数值积分法。

连续系统数字仿真的离散化方法有两类，它们是数值积分方法和离散相似方法，本章讨论数值积分法。

数值积分方法采用递推方式进行运算，而采用不同的积分方法会引进不同的计算误差，为了提高计算精度，往往会增加运算量。

就同一种积分算法而言，为提高计算精度，减小积分步距h ,计算量增大，影响系统运算速度。

因此，计算精度和速度是连续系统仿真中常迂到的一对矛盾，也是数字仿真中要求解决的问题之一。

也就是说，选择合适的算法、合适的软、硬件环境，在保证计算精度的前提下，考虑怎样提高仿真的速度。

经典的数值积分法可分为两类：单步法与多步法，下面我们将分别来介绍这两类方法中的最常用的算法。

2.2 龙格库塔法2.2.1龙格-库塔法基本原理在上一节中我们已经讲过，在连续系统的仿真中，主要的数值计算工作是对dydtf t y =(,)的一阶微分方程进行求解。

因为 y t y t f t y dt t t ()()(,)κκκκ+=++⎰11若令： y y t κκ≅(),Q κκκ≅+⎰f t y dt t t k (,)1则有 y t y y ()κκκκ++≅=+11Q (2.7)因此主要问题就是如何对Q κ进行数值求解，即如何对f t y (,)进行积分。

通常称作“右端函数”计算问题。

已知y t y ()00=，假设我们从t 0跨出一步，t t h 10=+，t 1时刻为y y t h 10=+()，可以在t 0附近展开台劳级数，只保留h 2项，则有：200010)(21),(h tf dt dy y f h y t f y y t ∂∂∂∂+++= (2.8) 在式(2.8)中括号后的下标0表示括号中的函数将用t t y y ==00，均同。

我们假设这个解可以写成如下形式：y y a k a k h 101122=++() (2.9) 其中， k f t y 100=(,)k f t b h y b k h 201021=++()，对k 2式右端的函数在t t y y ==00，处展成台劳级数，保留h 项，可得：h k b tfb y t f k t 0)(),(121002++≅∂∂ 将k k 12，代入(2.9)式，则有：])(),([),(012100200101h yfk b t f b y t f h a y t hf a y y t ∂∂∂∂++++= 将上式与(2.8)式进行比较，可得：可见有四个未知数a a b b 1212，，，，但只有三个方程，因此有无穷多个解，若限定a a 12=，则可得其中一个解：a ab b 1212121====，将它们代入(2.9)式可得一组计算公式：y y h10122=++()K K (2.10) 其中，K K K 1002001==++f t y f t h y h (,)()，， 若写成一般递推形式，即为：y t y y hk k k ()()++≅=++11122K K (2.11) 其中， K K K 121==++f t y f t h y h k k k k (,)(,)，由于(2.8)式只取了h h ，2两项，而将h 3以上的高阶项略去了，所以这种递推公式的截断误差正比于h 的三次方，又由于计算时只取了h 及h 2项，故这种方法被称为二阶龙格−库塔法（简称RK-2）。

根据上述原理，若在展成台劳级数时保留h h h ，，23及h 4项，那么可得一套截断误差正比于h 5的四阶龙格--库塔法（简称RK-4）公式：y t y y k k k ()()++≅=++++1112341622K K K K (2.12)其中： K 1=f t y (,)κκK K 2122=++f t h y h()κκ， K K 3222=++f t h y h()κκ，K K 43=++f t h y h ()κκ，由于这组计算公式有较高的精度，所以在数字仿真中应用较为普遍。

为了帮助读者更好地掌握这种方法的使用，下面我们对龙格--库塔法的特点作一些介绍。

1. 由于在解a a b b 1212，，，时，可以得到许多种龙格--库塔公式，经常使用的除(2.11)及(2.12)式给出的两组公式外，还有：y y h κκ+=+12K (2.13) a a a b a b 12212211212+===,/,/其中 K 1=f t y (,)κκ K K 2122=++f t h y h()κκ， y y hκκ+=++++11234833()K K K K (2.14) 其中 K 1=f t y (,)κκK K 2133=++f t h y h()κκ， K K K 312233=+-+f t h y hh ()κκ，K K K K 4123=++-+f t h y h h h ()κκ，(2. 13)式为二阶龙格--库塔公式，(2.14)式是Shampine 提出的，称为RKS3-4公式。