1、系统仿真技术 第6章 病态系统仿真,剡昌锋 刘军 兰州理工大学机电工程学院,6.1病态系统的定义,系统中各环节的时间常数差异巨大。为保证仿真计算的稳定性，由于仿真步长必须限制在最小时间常数的数量级而选得很小，然而仿真结束的时间则决定于系统中的最大时间常数，若按满足稳定性要求所选择的步长进行仿真，则不仅整个仿真所花费的时间非常长，甚至由于计算的舍入误差而导致整个仿真的失败。这就是所谓“病态系统仿真”问题。,病态系统的定义（续）,病态系统定义: （1） 令 （2）J称为系统的雅可比矩阵。 若J的特征值全部具有负实部，且有： ，则该系统称为病态系统，在某些文献中也叫做刚性系统(Stiff)，而: 称
2、为病态比，一般在50以上。,6.2 线性病态系统仿真,对线性定常系统，我们可用如下状态方程进行一般性描述：(3) 1. 增广矩阵法 将u(t)作为系统的增广状态,线性病态系统仿真(续）,其中: 由于病态系统特征值相差倍数很大，必须用加速收敛的方法计算该状态转移矩阵的值。 2.蛙跳算法 基本思想： （1）考虑作用函数为阶跃函数，则增广状态十分简单。 选择q，使得：,线性病态系统仿真(续）,(3) 仿真计算时采用如下“蛙跳”方式：,线性病态系统仿真(续）,即在qh以前采用加倍跳跃式计算，而在qh以后每隔qh计算一次。 优点是： h可以取得很小（可按最小时间常数考虑），从而保证初始阶段的精度而计算量
3、却不大，而到qh以后，小时间常数的作用完成，则加大步长计算，从而加快仿真计算速度。 从x(h)到x(qh)都是以x(0)为基础进行计算，所以误差传播比较均匀（仅仅是状态转移矩阵的误差）。,6.3 非线性病态系统仿真,一般非线性系统的仿真大多采用数值积分法。而数值积分法一般又只具有有限的稳定域，典型的如龙格库塔法，仿真步长限定在系统最小时间常数的数量级，才能保证计算的稳定性，而系统的过渡过程时间却决定于最大时间常数，因而对病态系统来说计算量极大，加上存在误差传播，仿真的精度甚至稳定性也会受到影响。,6.3.1 吉尔(Gear)法,6.3.1.1.Stiff稳定域 Gear研究后发现，并不要求一定
4、采用恒稳方法，而只要具有所谓Stiff稳定域就可以了。 Stiff稳定域定义：对实际的物理系统，时间常数一般小于零。选择仿真步长h若满足：与可保证仿真的稳定性，称该算法具有stiff稳定域。,Stiff稳定域（续）,实际上具有Stiff稳定域的方法与恒稳方法只在近虚轴处有一点差别，即如果系统中的极点全部为实极点，那么无论选择多大的步长，计算是恒稳的。如果系统中有复极点（实部仍为负数），只要步长的选择满足上述条件，也能保证算法稳定。,Stiff稳定域（续）,Stiff域中与的确定： 按病态系统的大特征值来选择步长：该特征值所对应的模态大约要经过4倍左右时间常数的时间才能有效地衰减掉，即，也就是
5、这样，此时即使加大步长h，也能保持计算的稳定性。基于这一考虑，可设- 4。,Stiff稳定域（续）,另一方面，考虑到系统特征值为复数，它所对应的瞬态响应呈振荡型。一个振荡周期内至少计算N个点。最小振荡周期为： 其中h为计算步长，若选择N8，则有： ，因此可选/4。综上所述，如果选择某一种方法，其稳定域/4，且|4，则从使用的角度来看，图6.1所示的稳定域与恒稳域没有差别，从而完全可以用于病态系统的仿真。,6.3.1.2吉尔（Gear）法的基本原理,设系统： 满足Stiff稳定域的多步法：Gear提出的用于病态系统仿真的计算公式是：(1),用于病态系统仿真的Gear公式的系数表,吉尔（Gear）
6、法的基本原理(续）,稳定域如图6.2所示。 从图上可以看出，该方法在5阶以下（包括5阶）的稳定域满足Stiff稳定域的条件（/4,|4，而且还可能穿过负实轴。,吉尔（Gear）法的基本原理(续）,在用Gear法仿真非线性病态系统时，有以下三个基本问题需要解决： 1）启动问题 上述Gear法本质上是隐式多步法。 对于初值问题，困难：隐式方法一般用显式方法启动，即先进行预报，然后通过迭代进行校正。如果迭代方法的收敛性不好，可能引起计算发散或计算量加大。 即使选择的迭代方法收敛性满足要求，显式多步法预报，仍然难以启动，必须采用单步法启动，由于单步法不具有Stiff稳定域，因而很难保证计算的稳定性。
7、2）变步长策略 非线性病态系统仿真往往采用变步长策略，如何适时地将步长调整到合适长度，以同时满足仿真精度和速度的要求。 3）加速迭代 为了提高计算效能，加速迭代也是非线性病态系统仿真中重要问题。,6.3.1.3 单步多值法,以三阶为例，采用显式多步法进行预报，然后用隐式法校正。 其显式预报的公式是： （2） 三阶隐式Gear公式校正：（3）,单步多值法（续）,其中等式右边第4项为导函数项，它是通过将第i次迭代所得到的y的预报值代入导函数后计算得到的。为便于程序实现，由（2）式，并令：,单步多值法（续）,迭代的校正公式可表示成：,单步多值法（续）,对一般情形，令：,单步多值法（续）,其中p表示G
8、ear法的阶次。对三阶预报迭代校正算法：用矩阵的方法加以表示：（4）（5） 其中B，C分别是相应的系数矩阵。,单步多值法（续）,然而，对初值问题，上式是不能自启动的。 解决方法-单步多值法：用高阶导数值来取代前几步的y及f的值。先定义一个向量，称之为Nordsieck向量：（6）（7） 需要确定Z向量与Y向量之间的关系。,单步多值法（续）,采用多项式逼近，以三阶为例：（8）在 处，有 （9）,单步多值法（续）,同样，也可以得到：由（8）式及（9）式消去 ，整理后可得到：,单步多值法（续）,写成矩阵形式，就是：简记为： （10）Q阵就是Y向量与Z向量之间的变换阵。对p阶Gear法，Q阵为（p+1
9、）阶非奇异方阵。,单步多值法（续）,将（10）式代入（4）及（5）式，可得到单步多值法的计算公式：(11) 说明：（11）式中计算 时以 为自变量，但由于 是一个标量，所以实际上只是用到它们的第一个分量 。,单步多值法（续）,若记 （12）则（11）式可简写为：（13）在三阶的情况下，P及L的值如下：,6.3.1.4 误差估计与控制,当用单值多步的Gear法对非线性病态系统进行仿真时，它要求从初值开始，必须依靠显式法来启动，即先从 及 开始，按一阶公式计算，然后逐次升阶，计算出 在这种升阶过程中必须满足误差要求，对k阶多步法，其截断误差为: (14) 其中 是t=点上y(t)的（p+1）阶导数
10、值，而是所讨论区间中的某一个点。,误差估计（续）,例如，对三阶Gear法，就是 区间上的某一个点。若设 ，则第k步的截断误差为： (15) 为了估计 ，首先要估计 ，已知：现在要用 的差分来近似估计 ，即：,误差估计（续）,两边同乘以 则可得：（16） 将（16）式代入(15)式可得： （17） 仿真中一般采用的是相对误差，若要求每一步的相对误差不大于 ，即 （18）其中 为到目前为止已出现过的y的最大绝对值（注：若y的初值为0，则应取 1为宜。,误差估计（续）,若系统为微分方程组，状态变量y为N个（y(1)，y(2)，y(N)），相对误差可定义为： （19）若 ，则本步计算结果有效，进入下一
11、步；如果不满足，则需要减小步长或者采取其它措施。,误差控制技术,误差控制：一是改变步长，其二是改变阶次。 无论是 或 ，均需要考虑变阶或（与）变步长，即通过改变p及h以使 。 变阶或（与）变步长的原则： 首先考虑仅改变步长h（阶次不变），设新的步长为 ，且 。那么， 应取多大为宜呢？,误差控制技术（续）,为简便起见，我们根据单变量表达式（15）来分析，即: 即根据（16）式， 故有：,误差控制技术（续）,可得到：考虑到误差仅仅是估计值，考虑经验系数1/1.2，这样， 表达式如（20）式：（20）若改变步长的效果不理想，则考虑要提升仿真方法的阶次。,误差控制技术（续）,设当前为p阶，考虑用p1阶
12、计算。假设此时步长 能使相对误差接近规定的要求，即：采用二阶差分来近似 ，即：,误差控制技术（续）,考虑到： 则不难得到：也考虑经验系数（这里取其为1/1.4，），则可得到升阶时 的表达式如（21）式：（21） 如果仿真步长缩短到相当小而其误差仍然达不到要求，则可能是因阶次过高而稳定域达不到要求的缘故。,误差控制技术（续）,若将当前仿真的阶次p降低一阶，仿真步长h变为 ，此时要使仿真误差接近规定要求，则：考虑经验系数为1/1.3，则可得到：(22) 式中 ，它是向量 的最后一个分量。,误差控制技术（续）,步长发生变化时，如何由当前的 产生在 下的 已知 只变步长：,误差控制技术（续）,降阶变步
13、长：从p阶降为p-1阶，步长由h变为 ，这时， 、 用 表示：,误差控制技术（续）,升阶变步长：当由p阶升到p+1阶，步长由h变为 ，则这时， 、 用 表示：然而， 是未知的，为此，必须由已知的数据来估计。若采用差分法，即：,误差控制技术（续）,则得到 的表达式如下：,6.3.1.5 加速收敛问题,Gear法用于病态系统进行仿真时，必须先用显式公式启动，然后用隐式法进行校正，经过多次迭代，达到适当精度后，再往前推进。显然，迭代的收敛速度极大地影响着仿真效能。 三阶Gear法为例(1) 其中 就是 的简写，它一般是 的非线性函数，（1）式是一非线性方程。,加速收敛问题(续）,迭代计算时，我们采用
14、的迭代公式为：（2） 其中 是向量L的第一个分量，这称为Picard迭代。对该迭代过程：（3）,加速收敛问题(续）,为保证收敛，则要求： 即： 一般 大约为13之间，可见，h不能太大，否则迭代将不收敛或收敛速度十分缓慢。这大大限制了Gear法的有效性。 牛顿迭代法： （1）式可改写为： （4）,加速收敛问题(续）,简记为： （5） 令 （6） 显然，（5）式与g(y)=0同解。由牛顿迭代法，不难得到：（7）,加速收敛问题(续）,从(6)式求得 的表达式，并将 代入 及（7）式，可得：(8)将 代入(5)式，所得结果在代入(8)式，可得：,加速收敛问题(续）,(9),加速收敛问题(续）,多了一个
15、“加速”因子： 对微分方程组来说，y是N维向量，(9)式成为如下形式：(10) 其中Y，F，G为N维向量，I为N阶单位阵，而 为NN的方阵，亦称为雅可比矩阵：,加速收敛问题(续）,(11) 采用牛顿法：(12)其中Z为（K1）1的向量，而G为标量。,加速收敛问题(续）,对于y为N维的情况，则有：(13)式中的Z为（K1）N的矩阵，而G为N1的向量。 牛顿法能加快迭代过程的收敛，这是该方法的优点，然而，因每次迭代时都必须计算雅可比阵，还要计算矩阵求逆，计算量将增加较多。,6.3.1.6 病态性探测,系统并不是每一步均处于病态。从数值计算的角度进行定义“病态性”和“非病态性”：如果用非Stiff法
16、仿真，其步长限制是由于计算精度原因引起的而不是由于稳定性引起的，则称系统呈现“非病态性”，反之，若步长是因稳定性而受到限制时，则系统呈现“病态性”。当系统呈现病态性时应采用Stiff法，而当系统呈现非病态性时应采用非Stiff法。,病态性探测 (续）,病态性探测的三种方法： 1）稳定半径法 非Stiff法的稳定域的形状接近半圆，其稳定域与实轴的交点称为稳定半径（ ）。若记为系统的最大特征值，为保证计算稳定，则要求： (1) 各种方法的稳定半径已知，而的值可由雅可比阵的范数 来估计，即： 因此，在用stiff方法对病态系统仿真过程中，雅可比矩阵的值可以得到，则在下一步计算时，先用（1）式进行判断
17、病态性，再决定是否仍采用stiff方法仿真。,病态性探测 (续）,2）嵌入低阶大稳定域法 如果当前采用的是非Stiff法，而计算误差不能满足要求，为判断其稳定性，可降低阶次，如果此时误差减少，则说明当前已处于非稳定状态（非Stiff法的低阶的稳定域大于高阶的稳定域）。因而下一步应采用Stiff法计算。,病态性探测 (续）,3）小稳定域试探法 嵌入低阶大稳定域法对于高阶的非Stiff法是有效的，然而，如果当前使用的已经是低阶方法，则很难再嵌入一个更低阶的算法。小稳定域试探法采用与嵌入低阶大稳定域法相反的思路来进行病态性探测。若当前采用的是非Stiff法，为探测下一步的病态性，步长不变，用比其高一阶的方法进行探测，若误差加大，说明当前已处于病态或病态边缘，因此下一步要么减少步长，要么采用Stiff方法计算。,
