系统时间响应和动态仿真

1、第五章 系统时间响应和动态仿真5.1 概述系统的时间响应是指系统在输入信号或初始状态作用下，系统的输出随时间变化的状况。系统的时间响应反映了系统的特征和性能。充分了解和分析系统的响应，以及与系统结构、参数之间的关系，是我们设计和校正系统的基础。低阶系统：用解析法实现；高阶系统：计算机仿真。仿真主要过程：建立模型、仿真运行和分析研究仿真结果。仿真运行借助一定的仿真软件，而这些软件是某些算法的实现。本章介绍算法的基本原理。5.2 基于数值积分的连续系统仿真5.2.1 数值积分的基本原理一阶微分方程：（1）0( )( , ( ),( )y tf t y tatby ay 根据常微分方程理论可知,如果

2、f(t,y(t)在区域D:at b,-y0时，局部截断误差Rn与h2为同阶无穷小。（4）整体截断误差整体截断误差是个积累误差，计算繁杂。但有如下结论：如果一个方法的局部截断误差为1()pO h则其整体误差将为()pO h，称该方法为p阶方法，该方法具有p阶精度。二 梯度法欧拉法是用前一点的斜率值f(tm)确定下一点的y(tm+1)的值，精度较低。为了提高精度梯度法是两个点斜率的平均值来确定下一点的y值，即),(),(2111mmmmmmytfytfhyy(8)由上式可见，这是隐式算法，不能自启动。在实际计算中，首先用欧拉法预估ym+1的值，再进行校正这样可得一个预估一校正算法：hytfytfy

3、yhytfyymmmmmmmmmm),(),(21校正：),(预估：1111(9)为了提高计算精度，也可采用多次迭代校正，即预估：hytfyymmmm),()0(1校正：hytfytfyyhytfytfyyimmmmmimmmmmmm),(),(21.),(),(21)1(11)(1)0(11)1(1直到)1(1)(1imimyy为止三 龙格库达法欧拉法取泰勒公式的前两项进行近似计算求微分方程的解。为了提高计算精度，必须取泰勒公式(3)的前若干项，但公式直接利用高阶导数，计算不方便。为了避免求函数的高阶导数，可用计算区间内几个点的导数值的加权线性组合的数值积分计算方法，称为龙格一库达法(以下简

4、称RK法)。取泰勒公式(3)前三项，有221),(21),(21hyfytftfythfyhyhyyymmmmmmmmmmm （11）设式(11)的解具有下面形式：)(),(),(2211112121KaKahyyhKbyhbtfKytfKmmmmmm（12）式中，a1，a2，b1，b2为待定系数。将K2用二元函数泰勒级数展开式展开，并只取前三项，则有11212(,)(,)(,)mmmmmmmmmmffyya hf tya h f tyb hb hf tytymmmmyfhKbtfhbytfK1212),(21211221221babaaa(13)将K1和K2代入式(12)，则有(14)比较式

5、(11)和式(14)可得(15)3个方程包含4个未知数，解不唯一。若设a1=a2，则1212121bbaa(18)将上面结果代入式(12)得(17),(),()(2121211hKyhtfKytfKKKhyymmmmmm该式为二阶RK法公式。这种算法的局部截断误差为.!41!3143 hyhyRn或)(3hORn(16)显式RK法的一般计算公式为11111(,)(,)2rmmiiimmrjmjmjsssyyK hKf xyKf ta h yb K hjp式中，i为待定权系数，j为使用K的个数(即级数)，Kj为所取各点导数f的值，aj，bjs为待定系数。（19）当r1时，式(19)变为欧拉法数值

6、积分公式hytfyymmmm),(1当r2 时，式(19)可写为11122122211(,)(,)mmmmmmyyK hK hKf tyKf ta h yb K这是二阶RK法公式，其局部截断误差为O(h3)，故也称为二阶二级RK法公式。1123121312(4)6( ,)1(,)22(,2)mmmmmmmmhyyKKKKf tyhKf tyK hKf th yK hK h当r3 时，式(19)可写为这是三阶RK公式，其局部截断误差为O(h4)。当r4 时，式(19)可写为),()21,2()21,2(),()22(6342312143211hKyhtfKhKyhtfKhKyhtfKytfKKK

7、KKhyymmmmmmmmmm这是四阶RK公式，其局部截断误差为O(h3)。式(21)在系统仿真的数值积分计算中的应用最为广泛，因为它计算精度和速度都较好。由于RK法的约束方程少于待定系数补充的约束条件不同，上面公式的系数则不同，RK公式则有不同形式，但这对最后的计算结果影响不大。（21）由以上公式可知，RK法的阶数是以泰勒公式的截取的项数来定义的，或者说以局部截断误差来定义的。若局部截断误差为O(h(p+1)则称为p阶RK法。此外，RK法级数是以式(19)所取K值的个数定义的，即公式所取导数计算点的个数来定义的。若取K值计算点有r个，则称为r级RK法。因此式(19)也称为r级RK法计算公式。

8、同一阶的RK法也可以选用不同级数，则同一阶的RK公式原则上也有多种形式，但其中有一种计算量最少。根据理论推导，对显式RK法，当阶数p4时，最小级数rmin=p；当阶数p4时，最小级数rminp。如二阶二级RK公式除式(17)外，还有)2,2(),(12121KhyhtfKytfKhKyymmmmmm（22）四阶四级RK公式除式(21)外，还有112341213124123(33)8( ,)(23)(,)332(,(3)33(,)mmmmmmmmmmhyyKKKKKf t yhhKf tyKhKf th yKKKf th yhKhKhK对于阶数和级数不相同的显式RK法的局部截断误差估计可采用下面

9、方法：若一个K阶r级的RK公式，其输出是ym+1，而一个(K-1)阶r级的RK公式，其输出是 则这一步的截断误差Em可用下式估计：)24( ,11mmmyyE1my计算时，显然两个RK公式的阶数不相同，但级数相同，因此估计Em时并不增加计算工作量。如Runge-Kutta-Fehlberg公式是一个五阶六级公式：611iiimmKchyy)25( , ),(11ijiijmimiKbhyhatfK6,.,2, 1i它的系数ai，bij，ci值如表费尔格公式系数表所示系数ci*是四阶六级RK公式中的系数，因此误差估计为61*11)(iiiimmmKcchyyE它的精度等级为四阶，故上面的计算公式

10、称为RK-45。此外，还有RKM-34，RKS-34，RK23等公式，这些RK基本算法在MATLAB的求解微分方程的函数中大多被采用。RK法除显式外，还有隐式RK法。RK法可以用定步距也可用变步距。四 阿达姆斯(Adams)法Adams数值积分法是多步法，且有显式积分和隐式积分两种。1阿达姆斯显式积分(Adams-Bashforth公式)对于微分方程(1)，其数值积分精确解形式可写成式(2)。为了近似计算积分 ，取r+1个点(tm，ym)，(tm-1，ym-1)，(tm-r，ym-r)构成多项式Pr(t)逼近f(t，y)，则可得1),(mmttdtytf1)(1mmttrmmdttPyy或).

11、(22101mrrmmmmmfbfbfbfbhyy式中，系数bj计算结果列于表2-2，为向后差分。（26）32112232112132mmmmmmmmmmmmmmmmffffffffffffffff一阶向后差分：二阶向后差分：三阶向后差分：将式(27)代入式(26)，可得Adams显式一般形式（27）).(1101rmrmmmmfffhyy（28）式中， 为显式Adams公式系数，如表2-3所示。i由式(28)可知，计算ym+1所需的全部信息均为已知，故为显式。2阿达姆斯隐式积分公式(Adams-moulton公式)式中，系数cj可计算结果列于表2-4。和阿达姆斯显式积分解法类似，为了近似计算

12、积分 可用r+1个点(tm+1，ym+1)，(tm，ym)，(tm-r+1，ym-r+1)构成一个多项式Qr(t)来逼近函数f(t,y)，则可得1),(mmttdtytf).(11101mrrmmmmfcfcfchyy将式(27)代入(29)可得Adams隐式一般形式：（29）).(110111rmrmmmmfrfrfrhyy式中，r1为隐式公式Adams系数，如表2-5所示。（30）由式(30)可知，计算ym+1时所需的信息包含(k+1)计算点的导数值fm+1，故为隐式。阿达姆斯法是多步算法，不能自启动，开始几步要用单步法直至所需的前几点信息全部获得后才开始正式启动阿达姆斯法。阿达姆斯法的截

13、断误差可用下式来估计：)31( ,)(2)2(rrrmhyBE式中， 为所取逼近多项式时刻范围内函数的导数，h为计算步长，Br为截断误差函数。)()2(ry显式阿达姆斯法和隐式阿达姆斯法的截断误差系数列于表2-6。显然，隐式Adams法比显式的精度高。3预估校正公式(Adams-Bashforth-Moulton公式)(0)101(1)( )( )( )101111(.),(32)(.),(33)rmmmmrmiiirimmmmrmyyh b fbfbfyyh c fcfcf显式预估隐式预估隐式Adams法具有较高精度，但要提供首次估值，这可以由显式Adams公式来完成，称为“预估”。然后，用

14、隐式Adams公式进行迭代运算，直至达到一定精度要求为止，这称为“校正”。y(i)m+1为ym+1的第i次迭代计算结果，或者将上式差分展开和整理并项后得这里,),(),(1)(11)(1)(11)(1mjimjimjimmimmmmfffytffytff（34）（35）（36）式中，y(i)m+1为函数f(t，y)的第(m+1)时刻点的迭代值。例如，当r=3时，预估校正公式为：).(校正：).(预估：110)(11) 1(1110)0(1rmmimmimrmrmmmmfrfrfrhyyfffhyy(0)1123( 1)( )1112(5559379)24(9195)24mmmmmmiimmmm

15、mmhyyffffhyyffff预估 校正 ，（37），（38），（39），（40）5.2.2 数值积分方法的选择一 计算精度 数值积分法的离散数值解只是精确解的近似，必然存在误差。数值积分计算的误差来自两个方面：一是舍入误差；二是局部截断误差。 舍入误差是由于计算机的位数有限，计算时必然舍去精确值的某个小值而引起的误差。舍入误差每次计算时均会发生，因此计算次数增加，会使舍入误差的积累值增大，舍入误差和计算步长h成反比。因为过小的计算步长会引起计算次数增加，从而使舍入误差增大。 局部截断误差是由积分方法和阶次的限制而引起的误差。 数值积分计算的综合误差为舍入误差和局部截断误差之和，如图3可见，

16、存在一个最优的积分步长h，使计算误差e最小。图3 误差与积分步长2 积分步长的选择与控制积分步长h选择应在保证数值积分计算稳定性和精度前提下，尽可能选择较大的积分步长，以减少仿真计算次数，缩短仿真时间。积分步长选择还与被仿真系统的快速性有关，有一些推荐的经验公式，如tr为系统阶跃响应上升时间；wc为系统幅值穿越频率。 10rth 或ch51，（41），（42）数值积分计算时，积分步长有固定步长和变步长两种工作方式。固定步长就是在整个仿真计算过程中，积分步长h始终保持不变。变步长就是在仿真积分计算的每一步，根据计算误差的大小改变步长h。变步长的目的是在保证一定的计算精度的前提下，尽可能地选取较大

17、的步长。变步长方法是：首先估计计算误差；判断误差是否在允许的误差范围内；若在允许误差范围内，则该步计算有效，否则计算无效，改变步长，重新计算。因此，在变步长的积分计算中，必须解决误差估计方法和步长调整策略两个问题。关于误差估计，在不同的数值积分方法中采用的算法不同。在变步距的RK法中，误差估计用一个K阶RK法计算值ym+1和一个(K-1)阶RK计算值maxyEmm式中，Em为估计误差；ymax为方程解最大值 ； 为一小正数。在变步长前一般还需确定最大允许误差和最小允许误差。 ，（43）1my之差来获取，正如式(24)所示。Adams法的误差用式(31) 来估计。同时，这些误差估计方法也不是唯一

18、的。在变步长时，一般采用相对误差作为步长调整的控制参数。相对误差计算式为步长调整策略通常有两种：对分策略和最优控制策略。(1)对分策略 若 ，则本步成功；令hm+1=2hm，做下步计算； 若 ，则本步成功；令hm+1=hm ，做下步计算； 若 ，则本步失败；令hm=0.5hm，重新计算。minmmaxminmmaxm(2)最优步长控制策略 变步长计算公式为式中，a为小于1的系数，常取0.8或0.9；p为数值积分方法的阶次。 为了控制步长的变化范围，可以同时给出限制条件，如1max1()pmmmhh10.110,mmmhhh（44）数值计算时，若 ，本步计算失败，重新按式(44)修正步长，再做本

19、步计算；若 ，本步计算成功，按式(44)计算下一步hm+1，做下一步计算。 maxmmaxm（45）三 积分计算稳定性 常微分方程的数值积分法，实质上就是将微分方程差分化，然后从初值开始，逐步进行迭代运算。显然，要保持迭代运算正常进行下去，首先必须保持这一数值解法的稳定性。用差分方程来代替微分方程，用差分方程的解作为微分方程的近似解这里存在一个所选定的差分方程对存在的扰动(如初值误差、舍入误差等)是否敏感，即在这些扰动影响下，计算误差是逐渐减小还是无限增大。若计算误差随递推计算过程逐渐减小或有限的，称为稳定的，否则是不稳定的。同一微分方程，由于采用的数值积分方法不同，差分方程不同，稳定性也不同

20、。而一个数值解是否稳定，决定于该差分方程的特征根是否满足稳定要求。为了说明这个问题，讨论一简单的微分方程式中，为方程特征根且具有下面形式)( tydtdy0j（46）（47）下面采用几种数值积法求解并进行稳定性分析。1、欧拉前差公式差分方程 （48）代入式（46）得 （49）两边进行z变换，则得 )(1mmmthfyymmmhyyy1)()1 ()(zYhzzY，（50）差分方程的特征根为z-(1+ h)=0，（51）由差分方程稳定性条件得 (52)将式(47)代入式（52）可得11hz11)(hja（53）即式（55）表明，在 平面上，欧拉前差法稳定区域为一个圆，圆心为（- ，0），半径为

21、很明显，欧拉前差算法的稳定条件是 ，（56）或 ，（55）222222()()1111ahhahhh1h1ah2Th2，（54），（57）（66）（1）降阶处理 将方程特征值中，| Rei |特别大的去掉，其理论根据是系统主导极点决定了系统动态响应特征。这样做虽然给数值积分带来方便，但仿真结果不能完全真实地反映系统开始阶段的动态响应特征。 （2）变步长积分 积分步长的自动调整，使仿真阶段运行处于优化状态。Gear 数值积分法被公认为解决病态系统仿真行之有效的方法。 5.2.3 微分方程数值积分的矩阵分析方法),.,(.),.,(),.,(2121222111nnnnnxxxtfxxxxtfxx

22、xxtfx),(XtFX 以上讨论的是一阶微分方程的数值积分法。一般来说，一个连续时间系统模型要用一个高阶微分方程加以描述。 对于一个n阶微分方程描述的连续系统，可以转化为用n个一阶微分方程组米表示，形式如：其中其中其中或写成矩阵形式为其中TnxxxX,.,21),.,(.),.,(),.,(),(21212211nnnnxxxtfxxxtfxxxtfXtF(69)(68)这样，前述的数值积分公式可改写为矩阵形式，如二阶RK公式(17)改写为),(),()(2121211hKXhtFKXtFKKKhXXmmmmmm四阶RK公式(21)可改写为),()2,2()21,2(),()22(23423

23、12143211hKXhtFKKhXhtFKhKXhtFKXtFKKKKKhXXmmmmmmmmmm式中，K1，K2，K3，K4均为n维列向量。对于线性时不变系统，n阶微分方程转化为状态空间表达式：式中，状态变量由于状态方程实质上是一阶微分方程组，则四阶RK公式(21)可改写为系统输出：DUCXYBUAXXTnxxxX,.,21)()()2()2()2()2()()22(2342312143211htBUhKXAKhtBUKhXAKhtBUKhXAKtBUAXKKKKKhXXmmmmmmmmmm)(mmmtDUCXY（70）5.2.4 数值积分方法的MATLAB函数ODE45为一种显式R-K(

24、4，5)公式，它属于单步法，即计算y(tm)的值，只需要前一时刻y(tm-1)的值。变步距数值计算，误差估计为Dorman-Prince公式。通常为首先选用的最好函数，用于求解非刚性微分方程。对于大多数问题能获得满意的解。ODE23为一种显式R-K(2，3)公式，采用Bogacki-Shampine公式，它也属于单步法，变步距，适用于求解非刚性微分方程。在允许计算误差较大和解具有轻微刚性方程时效果比ODE45更好。MATLAB工具箱提供了各种数值积分方法的常用函数，这些函数的基本功能是用数值计算方法求解常系数微分方程(ODE)或微分方程组。MATLAB提供了七种解函数它们是：ODE45，ODE

25、23，ODE113，ODE15s，ODE23t，ODE23tb。下面分别介绍它们的特点：ODE113为变阶Adams-Bashforth-Moulton PECE算法，适用于求解非刚性微分方程，在允许误差较严格的场合，它比ODE45更有效：它属于多步法，需要前几步值计算当前值。ODE15s是基于数值微分公式(NDFs)变阶变步距算法，它常比反向差分公式(BDFs)即Gear算法更有效。ODE15s函数允许在两种方法之间进行选择。它属于多步法，适用于刚性微分方程。当所求微分方程是刚性方程或用ODE45函数求解失败情况下，可尝试使用ODE15s函数求解。ODE23s基于二阶修正的Rosebrock

26、公式的一种算法。它属于单步法，因为在计算精度不高场合下，比ODE15s更有效。若在使用函数ODE15s无效情况下，可尝试使用ODE23s来求解某些刚性微分方程。ODE23t采用“自由插值”实现梯度规则算法，适用于存在中等刚性微分方程并要求解无数值衰减的情况。ODE23tb TR-BDF方法的一种实现即隐式RK法。第一阶段采用梯度法(TR)，第二阶段采用二阶BDF公式。结构上，两个阶段采用相同的迭代矩阵。 以上七种解常微分方程MATLAB函数的调用格式完全相同：),0, (,)0, (,optionsytspanFsolverYTytspanfsolverYTsolver为解函数名，即ode23

27、，ode45，ode113，ode23s，ode23t，ode23tb，ode15s；F为常微分方程(组)(或系统模型)文件名，字符串，它应该是MATLAB m-函数来描述。所有这些解函数用以求解一阶微分方程(1)或一阶微分方程组(69) 。高阶微分方程必须首先化为一阶微分方程组或状态方程。tspan=t0，tfinal，t0为积分时间初值，tfinal为时间终值。或写成给定时刻向量，即tspan=t0，t1，tfinal。y0为初值；options是选择项，由函数odeset设定。选择项包括相对允许误差RelTol和绝对允许误差AbsTol。缺省时，RelTol=10-3，AbsTol=10

28、-6。MATLAB中，每一步的计算误差必须满足 ，该步计算结果才被接受。 T为计算点的时间向量。Y为微分方程向量或矩阵。),(*max(ReabToliyabslTolei 5.3基于离散相似法的连续系统仿真基本思想：将连续系统离散化。方法：离散化处理通过采样保持器实现。MATLAB提供了离散化函数。离散时间系统的数学模型有差分方程、z传递函数和离散状态空间表达式，它们可以相互转换。离散时间系统的状态空间模型形式如下：)()()()()() 1(kDUkCXkYkGUkFXkX这个方程本身就是递推公式，当已知初始条件X(0)及输入U(k)，则可由上式计算各离散点的值。5.3.1离散化方法duB

29、eXeTkXTkATkTkA)()()(01) 1() 1(0) 1( 稍加处理后，即可得 duBekTXeTkXTkATkkTAT)(1) 1() 1()()( 由于上式中积分项与 k 无关，故令 k，且采用零阶保持器，即积分式内)(u不变，)(uu(kT)，则上式变为 )()()(01kTuBdekTeTkXtATATX 由以上分析可见，从幅频特性来看，不论采用哪一种保持器，都会引起离散化模型值衰减。但对一般信号而言，采用三角形保持器精度较高。从相频特性而言，零阶保持器和一阶保持器会产生相位滞后， 且一阶保持器引起更大的相位滞后， 而三角形保持器无相位滞后，这意味着，它可无失真地恢复原信号

30、。三角形保持器是较理想的保持器。但实际计算时，它需要后一时刻（k+1）T 的信号值，而这个值尚未计算出来，使用较困难。 连续系统离散化必须在系统中加入保持器，但加入保持器后必使得离散化模型较之原连续产生幅值衰减和相位滞后， 产生失真， 且往往导致模型的稳定性变差甚至不稳定。 因此，进行系统仿真时，尽可能减少离散化的环节数。 二、离散系统连续化 控制工具箱还提供了离散系统连续化的函数，其主要功能是将离散时间系统模型转换为连续时间系统模型。函数调用格式为 sysc=d2c(sysd,method) 其中，sysd 为原离散时间系统的 MATLAB 模型；sysc 为等价的边连续时间系统模型；met

31、hod 为模型转换方法，其意义和函数 C2D 相同，缺省时，采用零阶保持器。 三、离散系统重新采样 MATLAB 控制工具箱还提供了离散时间系统重新采样函数 D2D，其功能是产生一个和原离散时间系统采样周期不同的离散时间系统模型。函数调用格式为 Sys2=d2d(sys1,Ts) 其中，sys1 为原离散时间系统模型；sys2 为重新采样后离散时间系统模型；Ts 为新的采样周期。 5.4 系统仿真的MATLAB函数MATLAB控制工具箱提供了直接求系统对给定输入下的时间响应函数用于系统动态仿真。这些函数既适用于连续时间系统，也适用于离散时间系统。对于连续时间系统，这些函数均采用离散相似法，连续

32、系统离散化函数c2d已在这些函数内部调用，用户不必再进行模型转换。一 阶跃响应仿真函数函数step用来计算LTI的单位阶跃响应。适用于连续时间系统，也适用于离散时间系统；适用于SISO系统，也可用于MIMO系统。(1)基本调用格式 用于绘制系统单位阶跃响应曲线，分为：),(),()(TsysstepTfinalsysstepsysstepsys为系统模型；Tfinal为仿真终止时间；T为用户指定的仿真时间向量，对于离散时间系统T=T0：Ts：Tfinal，Ts为采样周期，对于连续时间系统T=T0：dt：Tfinal，dt为连续系统离散化的采样周期，T箱为仿真开始时间。(2)多系统阶跃响应调用格式 用于在同一幅图中绘制多个系统的单位阶跃响应曲线。式中，sys1，sys2，为系统模型；T为仿真时间向量，意义和上面相同。 这种调用格式，用户还可以定义每