1、系统仿真技术 第2章 经典的连续系统仿真建模方法学,陈无畏合肥工业大学机械与汽车工程学院,2.1 离散化原理及要求,问题：数字计算机在数值及时间上的离散性-被仿真系统数值及时间上的连续性? 连续系统的仿真，从本质上：对原连续系统从时间、数值两个方面进行离散化并选择合适的数值计算方法来近似积分运算 离散模型原连续模型？,相似原理,设系统模型为： ，其中u(t)为输 入变量，y(t)为系统变量；令仿真时间间隔为 h，离散化后的输入变量为 ，系统变量 为 ，其中 表示t=nh。如果 , 且 即 , （对所有n=0,1,2,） 则可认为两模型等价。,原连续模型,仿真模型,对仿真建模方法三个基本要求,（
2、1）稳定性：若原连续系统是稳定的，则离散化后 得到的仿真模型也应是稳定的。 （2）准确性：有不同的准确性评价准则，最基本的 准则是： 绝对误差准则： 相对误差准则： 其中 规定精度的误差量。,对仿真建模方法三个基本要求(续),（3）快速性：若第n步计算对应的系统时间间隔 为 计算机由计算需要的时间为Tn，若 Tn=hn 称为实时仿真，Tnhn称为超实时仿真 Tnhn 称为亚实时仿真，对应于离线仿真,数值积分算法,对 ，已知系统变量y的初始条件 要求y随时间变化的过程初值问题 计算过程：由初始点 的 欧拉法 对任意时刻tn+1 截断误差正比于,数值积分算法（续）,梯形法: 是隐函数形式。预报欧拉
3、法估计初值，校 正用梯形法校正：校正公式 预报公式 反复迭代，直到满足 经典的数值积分法分为两类：单步法与多步法,2.2 龙格库塔法,2.2.1龙格-库塔法基本原理 对若令：则有 的数值求解：称作“右端函数”计算问题。 在 附近展开泰勒级数，只保留 项，则有：,龙格-库塔法基本原理（续）,假设这个解可以写成如下形式： 其中 对 式右端的函数展成泰勒级数，保留h项， 可得：代入，则有：,龙格-库塔法基本原理（续）,将(2)式与(1)式进行比较，可得：四个未知数 但只有三个方程， 因此有无穷多个解。若限定 ，则 计算公式：其中,龙格-库塔法基本原理（续）,若写成一般递推形式，即为：其中（1）截断误
4、差正比于h3，称为二阶龙格-库塔法 （简称RK-2）。 （2）截断误差正比于h5的四阶龙格-库塔法（简称 RK-4）公式： 其中：,2.2.2龙格-库塔法的特点,1.形式多样性例： 非唯一解，可以得到许多 种龙格-库塔公式： (中点公式) 其中 各种龙格-库塔法可以写成如下一般形式：其中,龙格-库塔法的特点(续）,式中各系数满足以下关系 s称为级数，表示每步计算右端函数f的最少次数。 可以证明，1阶公式至少要计算一次，2阶公式 ; .；4阶公式 ；依此类推。有时为了某种特殊 需要，可以选择 的计算公式。,龙格-库塔法的特点(续）,2.单步法在计算 时只用到 ，而不直接用 等 项。优点：存储量减
5、小，可以自启动. 3.可变步长 步长h在整个计算中并不要求固定，可以根据精度 要求改变 ，但是在一步中，为计算若干个系数 ，则 必须用同一个步长h。,龙格-库塔法的特点(续）,4.速度与精度 四阶方法的h可以比二阶方法的h大10倍，每步计 算量仅比二阶方法大一倍，高于四阶的方法由于每步 计算量将增加较多，而精度提高不快。,2.2.3 实时龙格库塔法,实时仿真：要求仿真模型的运行速度往往与实际 系统运行的速度保持一致。一般的数值积分法难以满足实时仿真的要求，这 不仅仅是因为由这些方法所得到的模型的执行速度较 慢，而且这些方法的机理不符合实时仿真的特点。 考虑系统,实时龙格库塔法(续）,RK-2公
6、式如下：一个计算步内分两子步： tn时刻：利用当前的un，yn计算k1-计算一次右端函数f需 。 tn+h/2时刻：应计算k2，尽管此时yn +1/2已经得到，但un +1则无法得到。(若对un +1也进行预报加大仿真误差)。仿真执行延迟h/2输出要迟后半个计算步距。,实时龙格库塔法(续）,实时龙格库塔法(续）,实时2阶龙格库塔法：tn时刻，应计算k1，利用当前的un，yn，需要 ； tn+h/2时刻，应计算k2，此时yn +1/2已经得到，un +1/2也 可得到，k2的计算就不会引入新的误差。计算一次右 端函数 需要 ，可实时输出yn +1。,实时龙格库塔法(续）,2.3 线性多步法,2.
7、3.1 线性多步法基本原理 基本原理：利用一个多项式去匹配变量若干已知 值和它们的导数值。 设： 时刻的 和 已知;预报：由 和 来计算校正：若 也已知，由它们来计算,线性多步法(续),采用的多项式具有以下形式(m阶) 其中： 是待定系数， 在 时刻， ，可得到：（2-1）,线性多步法(续),由 和 确定 ， 需要m+1个独立方程。该m+1个方程可由以下 等式导出：（2-2）,线性多步法(续),1、预报公式令m=2k-1，从（2-2）式得到如下方程组：,（2-3）将其写成矩阵形式： （2-4） 其中上标p表示预报。其解为： （2-5）,由于 为常数阵，其逆存在，Z向量中的各元素为已知值，因而d
8、向量的各元素值可计算得到，从而由 ,得到下一时刻的预报值。,缺点：只有 是所需要的，其它元素的计算成为多余 ,得不到 与 和 显式表达式 。,线性多步法(续),定义： （m+1）1的列向量(2-6) 定义辅助变量 (2-7) 此式可改写为 (2-8) 向量 的元素可划分为两个组(2-9),例：k=3，则(2-8)式为：可计算得到：,只依赖于k，即先前 和 的个数，而 与它们的数值无关。这样，可以预先求解（2-8）式 得到从而得到 的显式表达式：,例：试推导用 预报 公式 条件：已知 由,2、校正公式,预报公式显式公式，未包括 。校正：对该预报值应进行校正，即先预报 得到 ，然后再用此值推出 。
9、由 和 以及 来预 报 ，可令m=2k-1，从（2-2）式得到如下 方程组：,将其写成矩阵形式： （2-10）,校正公式(续),其中上标c表示校正，可得(2-11),校正公式(续),定义： 为（m+1）1的列向量， 上标T表示转置。将 左乘（2-10）式可得：（2-12）定义 （2-13）可改写为 （2-14）（2-15）,例：k=3同样， 只依赖于k，即先前 和 的个数， 而与它们的数值无关。这样 （2-16）从而 （2-17）,例：已知 ，预估 ，然后用校正 。 预估 预报公式为 校正校正公式为,2.3.2 线性多步法误差分析,为了便于分析，对预报公式和校正公式，定义统一的表达 式：（2-
10、18）-显式预报 显式预报时称为后向差分公式（BDF）同时均不等于0时为隐式校正公式 k称为公式的阶次。假设变量各时间的精确值已经得到，将其代入（21）式，可得： （2-19）,线性多步法误差分析（续）,在 附近，将每个函数展开成泰勒级数：（2-20）对所有i(i=0,1,2,k)，将（2-20）式代入（2-19）式，合并 同类项，可得（2-21）,线性多步法误差分析（续）,其中（2-22）如果 均为0，则称为p阶公式（2-23）,线性多步法误差分析（续）,下面，如果我们能证明上一节推导出的公式若能满 足求 均为0的条件，则就得出了这些公式 的截断误差满足（2-23）式。以三阶公式为例，将（2
11、-22）式与 相关表达式表示成右端为零，可得：,先讨论预报公式，由于 和 ，这意味着要将上述矩阵的第1列移到等式的右边，并去掉第5列。为了使矩阵成为方阵，将其最后两行也去掉。,其结果与用于推导预报公式的矩阵方程完全一样。,对校正公式，可采用类似的办法，只是 ，这样要将第 第5列移到等式的右边；若假定 ，则可去掉第4列。同样 为得到方阵，去掉最后两行，结果就是3阶校正公式。,这就表明，上一节导出的预报与校正公式的截断误差系数可以用下式来计算,（2-24）,2.4 稳定性分析,仿真方法选择的基本要求：仿真计算不改 变原系统的绝对稳定性。原系统是稳定的。观察欧拉法仿真递推公式故有 （2-25） yn (n=0,1,2,)为它的一个仿真解。,稳定性分析（续）,设 为其准确解，即（2-26） 用(2-26)式减去(2-25)式，可得： 即 （2-27） 特征方程为 （2-28）显然，为了使扰动序列n不随n增加而增长，必须 要求： 我们称它所对应的域就是该算法的稳定域： h 21/，即h小于等于系统时间常数的两倍。,确定数值积分法稳定域的一般方法,测试方程：数值积分公式 (2-29) 其中 是一个关于 高阶多项式函数， 则只有当 时，算法才稳定。,
