系统仿真技术第6章--病态系统仿真

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