连续系统仿真的基础算法

1、2.1 2.1 数值积分算法数值积分算法 2.2 2.2 数值积分算法的基本分析数值积分算法的基本分析 2.3 2.3 连续系统仿真的离散相似算法连续系统仿真的离散相似算法2.4 2.4 常用快速数字仿真算法常用快速数字仿真算法2.5 2.5 实时数字仿真算法实时数字仿真算法 小结小结3.1.1 数值积分算法的基本原理3.1.2 欧拉法3.1.3 龙格-库塔法3.1.4 微分方程数值积分的矩阵分析 上一页下一页返回连续系统仿真的数值积分算法就是利用数值积分法将常微分方程（组）描述的连续系统变换成离散形式的仿真模型差分方程（组）。为了能在计算机上进行求解，首先要把被仿真系统的数学模型表示为一阶微

2、分方程组或状态空间模型。上一页下一页返回 假设一阶向量微分方程及初值问题为00)(),(yyyfytt可以写成一阶微分方程组和初值问题形式00200210012121222111)()()(),(),(),(nnnnnnnytyytyytyyyytfyyyytfyyyytfy上一页下一页返回 (2.1)(2.2)110,kkttttt110d),()(d),()()(01kkkttkttktttttttyfyyfyy在上的解析解为希望用公式 kkkqyy11kykykq)(1kty)(kty1d),(kktttt yf来代替解析解，其中 ， ， 分别是 ， ， 的近似解。 上一页下一页返回 (

3、2.3)(2.4) 所谓数值解法就是寻求初值问题的解在一系列离散时间点 121,kktttt,21yy1,kkyy上的近似解。相邻两个离散时间点的间距 kktth1称为步长或步距 。数值积分法的主要问题归结为如何对向量函数 ),(yf t进行数值积分，求出 ),(yf t在区间 ),(1kktt上定积分的近似解 kq ，并且数值方法的共同特点是步进式的（即所谓递推式的）。 上一页下一页返回 常用的基本算法可以分为单步法、多步法和预估-校正法三大类，而每一类算法又可以分为显式算法和隐式算法两种类型。 先从标量形式开始讨论。考虑一阶微分方程及初值问题 00)(),(ytyytfy 解析解为 11d

4、),()(d),()()(001kkkttkttktytftytytftyty希望用公式 kkkqyy1来代替解析解，其中 分别是 , , 的近似值。kkkqyy,1)(1kty)(kty1d),(kktttytf 上一页下一页返回 (2.5) (2.6) (2.7)00)(),(ytyytfy ),(1kktt 欧拉（Euler）法是一种最简单的数值积分算法，对于 在区间 上求积分，有 1d),()()(1kkttkktytftyty若区间 ),(1kktt足够小，则 ),(1kktt上的 ),( ytf可以近似地看成是常数),(kkytf上一页下一页返回 (2.8)11),()(kkkkk

5、yythfyty即有),(kkkythfq 欧拉法为 kkkhfyy1式中),(kkkytff 上一页下一页返回 (2.9) (2.10) 欧拉法的图形表示： 为图2.1中的非曲面面积actk+1tk，欧拉法用矩形面积abtk+1tk代替。1d),(kktttytf显然，欧拉法仅适用于步长h很小的场合。 上一页下一页返回图2.1龙格-库塔（Runge-Kutta）法是求解常微分方程初值问题(2.5)式的各种数值积分算法中应用得最广泛的一种，包括许多不同的公式。思路：用若干个时间点上 f 的函数值的线性组合来代替 f 的各阶导数项，然后按泰勒公式展开确定其中的系数。 上一页下一页返回 RK法包含

6、有显式、隐式和半隐式等算法。仅介绍显式RK法。 一般形式： rickahyhctfkkWhyyijjijkikiriiikk, 2 , 10),(11111 式中， Wi为待定的权因子；r为使用的k值的个数；ki为不同时间点上导数 f 的值；ci，aij为待定系数。上一页下一页返回 (2.13)当r =1时，只有一个k1，就得到了欧拉法 kkkhfyy1当r =2时，),(),()(12122122111khayhctfkytfkkWkWhyykkkkkk将 kfk 1代入 2k中，并将 2k在 ),(kkyt 附近用泰勒级数展开,可得)()(),(),(),(221222122122hOff

7、 haf hcfhOyfytfatfchytffhayhctfkkykkkttyykkkkkkkkk上一页下一页返回 (2.15)(2.14)于是，有22111khWkhWyykk)()(3212222221hOffaWhfcWhfWWhykykkkk 将其与泰勒级数式逐项进行比较，令h，h2项的对应系数相等，得到以下关系 212112122221aWcWWW取c2=1，有 2121WW121a上一页下一页返回 (2.16)(2.17)从而，有RK2法 ),(),()(2121211hkyhtfkytfkkkhyykkkkkk当r=3时，可得RK3 法)32,32()3,3(),()3(423

8、121311hkyhtfkkhyhtfkytfkkkhyykkkkkkkk上一页下一页返回(2.18) (2.19) 当r=4时，可得到如下著名的四阶RK法（亦称为经典RK法，简记为RK4法）。 ),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyykkkkkkkkkk 该算法是数字仿真中最常用的一种算法。计算量较大，但计算精度较高。在实际仿真中得到了广泛的应用。 上一页下一页返回 (2.20)几种数值积分法都可以看成是 在 附近用泰勒级数展开而产生的，只不过是泰勒级数所取项数的多少不同而已。欧拉法只取前两项， RK2

9、法取了前三项，RK4法取了前五项。从理论上讲，可以构造任意高阶的计算方法。但是，精度的阶数与计算函数f值的次数之间的关系并非等量增加的。 )(tyktt 上一页下一页返回每步计算f 的次数234567r8算法精度阶数234456r-2由此可见，RK4法有其优越性。 表2.1 f 的计算次数与算法精度阶数的关系 上一页下一页返回 对于一阶向量微分方程及初值问题 00)(),(yyyfyttRK4法的向量形式为),()2,2()2,2(),()22(6342312143211KyfKKyfKKyfKyfKKKKKyyhhthhthhtthkkkkkkkkkk上一页下一页返回 (2.21) (2.2

10、2)具体写出来就是上一页下一页返回),()2,2,2,2()2,2,2,2(),()22(63,23, 213, 142,22, 212, 131,21, 211, 12, 2, 114321,1,nknkkkiinknkkkiinknkkkiiknkkkiiiiiikikihkyhkyhkyhtfkkhykhykhyhtfkkhykhykhyhtfkyyytfkkkkkhyy(2.23) 在控制系统的仿真中，最常见的向量微分方程是线性定常系统的状态方程 uBAxx即有ubbbxxxaaaaaaaaaxxxnnnnnnnnn212121222211121121RK4法的4个iK可表示为)()2

11、(2)2(2)(3423121htuhhtuhhtuhtukkkkkkkkBKxAKBKxAKBKxAKBAxK上一页下一页返回 (2.24)(2.25) (2.26)3.2.1 单步法和多步法3.2.2 显式算法和隐式算法3.2.3 截断误差和舍入误差3.2.4 数值积分算法的计算稳定性3.2.5 数值积分算法的计算精度、速度、稳定性与步长的关系3.2.6 数值积分算法的选择原则3.2.7 误差估计与步长控制3.2.8 数值积分算法仿真实例上一页下一页返回上节中介绍的几种数值积分算法都有一个共同特点：在计算yk+1时只用到了yk ,而不直接用yk-1, yk-2 ， 等项的值，即在本次计算中

12、仅仅用到了前一次的计算结果。这类算法称为单步法。单步法运算有如下优点：需要存储的数据量少，占用的存储空间少；只需要知道初值，就可以启动递推公式进行运算，即可以自启动；容易实现变步长运算。 上一页下一页返回与单步法相对应的还有一类数值积分算法，称为多步法。在它的计算公式中，本次计算不仅要用到前一次的计算结果，还要用到更前面的若干次结果。例如AB4法:)9375955(243211kkkkkkffffhyy式中 3 , 2 , 1 , 0),(iyihtffikkikl 多步法与单步法相比，欲达到相同的精度，计算工作量要少得多。因此，在相同条件下多步法比单步法要快。多步法的主要缺点是，它不是自启动

13、方法，必须用其它方法求初始几步的值。另外，多步法不容易实现变步长运算。 上一页下一页返回 (2.28)如果数值积分算法在计算yk+1时所用到数据均已求出，则称其为显式算法。 如果数值积分算法的计算公式右端中隐含有未知量 yk+1，则称其为隐式算法。例如AM4法就是隐式算法 : )5199(242111kkkkkkffffhyy上一页下一页返回(2.29)隐式算法也不是自启动的算法，需要用另一个显式算法估计一个初值，然后再用隐式算法进行迭代运算，这就构成了预估-校正算法。 四阶Adams预估-校正算法： )5199(24)9375955(2421)0(11321)0(1kkkkkkkkkkkkf

14、fffhyyffffhyy上一页下一页返回 (2.30) 在分析数值积分算法的精度时，通常用泰勒级数作为工具。假设前一步求得结果是准确的，即有 )(kktyy 则用泰勒级数求得在1kt处的解析解为 )()(!1)(! 21)()()(1)(2rkrrkkkkhOtyhrtyhtyhtyhty 不同的数值积分算法相当于在上式中取了不同项数之和而得到的近似解。 上一页下一页返回 (2.31)(2.32)这种由数值积分算法单独一步引进的附加误差称为局部截断误差。它是数值积分算法给出的解与微分方程的解析解之间的差，故又称为局部离散误差。步长h越小，局部截断误差就越小。 若数值积分算法的局部截断误差为

15、，则方法为r阶的。)(1rhO 算法的阶数可以作为衡量算法精确度的一个重要标志。上一页下一页返回 以上的分析是在假设前一步所得结果是准确的前提下得出的，即 )(kktyy 成立时，有 )()(111rkkhOyty 但在求解过程中，实际上只有当k=0 时， 才成立。而当k=1,2,3,时， 并不成立。因而r阶算法求得的解的实际误差要大一些的。 )(kktyy 设yk是在无舍入误差情况下由r阶方法算出的微分方程式的近似解， y(t)为微分方程的解析解， 为算法的整体截断误差，则可以证明， 。即整体截断误差比局部截断误差低一阶。 kkkyty)()(rkhO上一页下一页返回(2.33)1. 计算稳

16、定性的概念 2. 欧拉法的计算稳定性 3. 龙格-库塔法的计算稳定性 4. 多步法的计算稳定性 上一页下一页返回连续系统的数字仿真，实质上就是将给定的微分方程（组）变换成差分方程（组），然后从初值开始，逐步进行递推计算。 【例2.1】 已知微分方程及初值 5 . 1031)0(30tyyy 试比较在取不同步长时，其精确解与数值积分算法解之间的差异。 上一页下一页返回 【解】该初值问题的解析解为 tety3031)(上一页下一页返回对应的结果曲线如图2.2（a）所示。图2.2 (a) 解析解5 . 0 ,075. 0 , 1 . 0h5 . 1t)(ty取，分别用欧拉法和RK4法计算处的，所得结

17、果曲线如图2.2 （b）（g）所示。 图2.2 (b) 欧拉法（h=0.1） 图2.2 (c) RK4法（h=0.1） 上一页下一页返回图2.2 (d) 欧拉法（h=0.075） 图2.2 (e) RK4法（h=0.075） 图2.2 (f) 欧拉法（h=0.05） 图2.2 (g) RK4法（h=0.05） 上一页下一页返回 从图中可以看出，解析解单调下降并迅速收敛到0。 当 时，欧拉法和RK4法的解曲线均发散，数值积分算法的解是错误的。 1 . 0h 当 时，欧拉法的解曲线仍然发散，对应的解是错误的；RK4法的解曲线单调下降并收敛到0，对应的解是正确的。 075. 0h 当 时，欧拉法和R

18、K4法的解均收敛到0（虽然欧拉法的解曲线是振荡收敛的）。如果只要求得到 处的 的解，则两种数值积分算法的解都可以认为是正确的。事实上，此时有 05. 0h5 . 1t)(ty 上一页下一页返回欧拉法 RK4法 10101044085. 3)5 . 1 ( 4)5 . 1 (y解析解 21105417286. 9)5 . 1 (y为什么会出现上述情况呢？这是因为数值积分算法只是一种近似积分方法，在反复的递推计算中会引进误差（这种误差通常是由初始数据的误差及计算过程中的舍入误差产生）。如果误差的累积越来越大（见【例2.1】）,将使计算出现不稳定，从而得到错误的结果。系统的

19、稳定性与计算稳定性是两个不同的概念。由于选用的数值积分算法不同，即使对于同一系统，差分方程也各不相同，计算稳定性也就各不一样。上一页下一页返回通常用一个简单的一阶微分方程来考查数值积分算法的计算稳定性。 微分方程及初值问题0)0(yyyy称为测试方程（Test Equation），其中， 为方程的特征根。根据稳定性理论，当特征根位于左半s平面，即 时，原系统稳定。 j0Re上一页下一页返回 (2.34) 用欧拉法对测试方程进行计算，有 , 1 , 0)1 (1kyhhyyykkkk 假定仅仅在初始值y0处引进了初始误差 ，而在递推计算过程没有引进任何其它误差。显然，此时 的误差01ky 1k仅

20、由 ky的误差 k引起的，所以有 )(1 ()(11kkkkkkkkyhyhyy用（2.36）式减去（2.35）式，得到误差方程 （2.35）（2.36）上一页下一页返回kkh)1 (1以此类推，有 0121)1 ()1 ()1 (kkkkhhh 当 时， 。表明若在递推计算过程中的某步引入了误差，则随着递推步数的增加，这个误差将逐渐扩大，最终导致差分方程的解完全失真。 11 hkk1 反之，当 时， 。表明随着递推步数的增加，引入的误差会被逐渐减小或保持有界。在这种情况下，称差分方程是计算稳定的。 11 hkk1上一页下一页返回 (2.37) (2.38) 显然，对于欧拉法，合理地选择步长h

21、使其满足 是保持其计算稳定性的充要条件。 11 h 令 ，由 ，可知欧拉法在实轴上的稳定区域为（-2，0），即有 hh11 h2hh因此， 越大，步长h就应该取得越小。 Re 为了保证计算稳定性而对步长h有所限制的数值积分算法称为条件稳定算法。在任何步长h0的情况下计算都稳定的数值积分算法称为绝对稳定算法（亦称为无条件稳定算法或恒稳算法）。欧拉法是一种条件稳定的计算方法。 上一页下一页返回 将RK法应用于测试方程，容易得到下列误差方程 )()(!1)(! 211121rkrkhOhrhh令 ，代入上式，可得该式稳定的条件为 hh1!1! 21121rhrhh可得如表2.2所示的各阶RK法在实轴

22、上的稳定区域。 上一页下一页返回 （2.40）(2.39)显式RK法都是条件稳定算法，步长h的大小除了与所选用的算法有关外，还与方程本身的性质有关。对于实际系统，由于其特征值不一定为实数，因此满足（2.40）式的也是复数。一般而言，步长h必须满足不等式 (2.41) 1hh12/12hh 6/2/132hhh24/6/2/1432hhhhr 所在的实负稳定区域1（-2，0）2（-2，0）3（-2.51，0）4（-2.785，0） 表2.2 RK法在实轴上的稳定区域 32h上一页下一页返回将多步法应用于测试方程，可以得到类似的误差方程。结论：在同阶的多步法中，隐式算法的稳定区域比显式算法的大得多

23、；随着阶数的增加，多步法的稳定区域逐渐减小。 上一页下一页返回方法的阶数1234显式（-2，0）（-1，0）（-6/11，0）（-3/10，0）隐式（-，0）（-，0）（-6，0）（-3，0） 表2.3 Adams法在实轴上的稳定区域 除了AM1法和AM2法（隐式Adams法）为绝对稳定算法外，其它算法都是条件稳定算法。这就是说，步长h必须满足不等式 式中，为由积分算法确定的常数。上一页下一页返回|h|plot(tout,yout,k);grid; 上一页下一页返回图2.8 (exam2_4.mdl) 得到如图2.9 (a)所示的阶跃响应曲线。显然，该曲线不是很理想，超调量较大。为此，可以将外

24、环的PI控制器参数调整为 ，并分别选择=0.17，0.5，1，1.5，可以得到如图2.9 (b)所示的仿真结果。可以看出，如果选择PI控制器为 ，就能够得到较为满意的控制效果。sas85. 0) 1(ss85. 015 . 1图2.9 （a）图2.9 （b）上一页下一页返回 【例2.5】考虑著名的Van de Pol方程 0) 1(2yyyy 25. 0)0(, 0)0(yy试绘制其相轨迹（ ）。【解】选择状态变量 yxyx21则有如下非线性状态方程 1221221) 1(xxxxxx上一页下一页返回200 t 第1个方程可以看成是将x2(t)信号作为一个积分器的输入端，积分器的输出将成为x1

25、(t)信号。同样，x2(t)信号本身也可以看成是一个积分器的输出，而积分器的输入端信号应该为 。于是，利用Simulink提供的各种模块可以得到如图2.10所示的仿真结构。 1212) 1(xxx上一页下一页返回 图2.10(exam2_5.mdl) Simulink模型中有些模块需要将输入端和输出端（通常用于反馈路径）掉换一下方向。为此，可以用鼠标单击需要掉换方向的模块选中它，则选中的模块的四个角出现黑点，表明它处于选中的状态。然后， 打开Simulink的Format菜单(见图2.11)，选择其中的翻转子菜单（Flip block）即可。 上一页下一页返回图2.11 启动仿真后，仿真结果将

26、赋给MATLAB工作空间内的变量t，x1和x2。在MATLAB命令窗口中运行指令 plot(x1,x2,k);grid; 得到相轨迹，如图2.12所示。 上一页下一页返回图2.12Simulink模型中的参数可以是实际数值，也可以是用字母表示的变量名。用字母表示的变量可以在MATLAB的工作空间中赋值，或用M文件赋值，然后进行仿真计算。 上一页下一页返回 【例2.6】 含有磁滞回环非线性环节的控制系统如图2.13 (a)所示，其中，磁滞回环非线性环节的特性如图2.13（b）所示。试研究该非线性环节对系统性能的影响。 上一页下一页返回图2.13 【解】 构建系统的Simulink模型，如图2.1

27、4所示。为了便于研究问题的方便起见，将Backlash模块的参数deadband width设置为C1。 在MATLAB命令窗口依次运行C1的不同值(C1=0.5,1,2)的指令后启动仿真，并在仿真结束后在MATLAB命令窗口中运行指令plot(tout,yout,k); 上一页下一页返回图2.14 (exam2_6.mdl)得到如图2.15所示的仿真曲线。显然，不同磁滞宽度的阶跃响应曲线的形状不相同。需要注意的是，为了将不同C1值对应的阶跃响应曲线绘制在同一坐标轴下，在第一次绘制图形后，应该在MATLAB命令窗口中运行指令hold on 上一页下一页返回图2.15 (r(t)=21(t) 此

28、外，如果输入的幅值增大或减小，则原来系统响应曲线的形状将可能出现不同。图2.16为C1=1时，输入幅值等于3和0.6的仿真结果。显然，这一点与线性系统不相同。 上一页下一页返回图2.16如果Simulink模型中存在复杂的非线性环节或复杂的逻辑运算，而在MATLAB提供的所有工具箱中都找不到相应的模块，可以自己编制一个M函数，嵌入Simulink模型中。上一页下一页返回 【例2.7】 某非线性系统如图2.17所示，试求r(t)=2.1(t) 时系统的动态响应。 【解】构建系统的Simulink模型，如图2.18所示。为了便于研究问题的方便起见，不采用Discontinuities 子库中的Sa

29、turation模块，而选择User-Defined-Functions子库中的MATLAB Fcn模块, 并将参数MATLAB Function设置为：saturation _zone。 上一页下一页返回图2.17返回上一页下一页返回图2.18 (exam2_7.mdl)然后编制M函数内容如下：% saturation_ zone functionfunction uo = saturation_ zone (ui)if ui=1 uo=1;elseif uia=0 0;1 -1; b=1;0; c=0 1; d=0; T=0.05; a1,b1,c1,d1= c2dm(a,b,c,d,T,

30、 zoh); 并分别键入K=1，10，30。选择定步长的discrete（no continuous states）算法，步长为0.05。启动仿真后，得到如图2.27所示的仿真结果。 )(T)(T上一页下一页返回 图2.26 (exam2_8.mdl)图2.27(a) K=1图2.27(b) K=10图2.27(c) K=30上一页下一页返回 在许多实际应用场合下，希望能尽可能快地获取仿真结果。这就构成了实际的快速需求和仿真速度之间的矛盾。当采用普通微型计算机对某些控制系统实施仿真时，这种矛盾往往变得非常突出，从而提出了快速仿真的需求。 快速数字仿真算法的种类很多，它们各有优缺点，目前还没有一

31、种算法能够适用于所有问题。本节将介绍常用快速数字仿真算法。上一页下一页返回3.4.1 仿真中对快速性的需求3.4.2 相匹配原理 3.4.3 替换法3.4.4 根匹配法 上一页下一页返回 通常需要采用快速仿真的问题可以分为下列几个方面： 利用仿真技术进行控制系统的参数优化设计时，需要反 复地对系统进行仿真计算，以获得满意的工程参数； 在数学-物理混合仿真中，由于有实物系统介入仿真模型，要求仿真模型的时间比例尺与真实系统的时间比例尺完全相同（即所谓实时仿真），如果系统比较复杂或者方程个数很多，要在规定时间内完成一次系统中每一个方程的计算，就要求仿真的计算速度比较快；上一页下一页返回 在复杂系统的

32、控制中，常常需要在线用仿真方法对被控系统的状态进行预测，以确定系统的控制策略，这就要求仿真模型的时间比例尺小于系统的运行时间比例尺（即所谓超实时仿真），从而对仿真速度提出了更高的要求。 一般而言，对于第一种情况，仿真速度慢意味着不能及时获得仿真结果，并且占用大量的机时，费用较高，效率低下；而对于后两种情况，仿真速度慢意味着仿真的目的达不到或者失去实际应用价值。上一页下一页返回 算法的计算精度和速度是一对矛盾，在研究快速数字仿真算法时，计算速度成为矛盾的主要方面。通常在满足计算稳定性及工程精度要求的前提下，应当尽可能提高仿真计算速度。 一般快速数字仿真算法有如下两点基本要求： 每步计算量要小；

33、算法具有良好的稳定性，允许采用较大的计算步长，同时又能保证必要的计算精度。 应当根据实际情况，尽可能地减小仿真中附加信息的计算和处理。 上一页下一页返回相匹配的含义是，如果被仿真系统的数学模型是稳定的，则其仿真模型也应该是稳定的，并且二者的瞬态、稳态特性一致。如果对于同一输入信号，二者的输出具有相一致的时域特性，或者二者具有相一致的频率特性，则称仿真模型与原系统模型相匹配。将高阶系统的数学模型G(s)转换成与之相匹配、每步计算量较小、允许采用较大步长且具有合理精度的仿真模型G(z) ，就可以利用对应的差分方程进行快速仿真。从G(s)直接推导出G(z)的方法有两种。一种称为替换法，另一种称为根匹

34、配法。 上一页下一页返回1替换法的基本思想2双线性变换法3双线性变换法步骤4双线性变换法特性分析5双线性变换法仿真实例上一页下一页返回根据控制理论，s域和z域之间的准确映射关系为 （2.67）或 （2.68）式中，T为采样周期，亦为仿真计算时的步长。（2.68）式是超越方程，很难由 得到一个关于变量u和y的线性差分方程。 Tsez zTsln1)()()(zUzGzY上一页下一页返回ln(z)可以展开成下列的无穷级数 （2.70） 取该级数的第一项，则可以得到s域与z域之间一种近似映射关系 （2.71）或 （2.72）利用（2.71）式可以把G(s)转换为G(z)。数学上称这种变换方法为双线性

35、变换法（或Tustin法）。 121233) 1() 1() 12(1) 1() 1(31112lnmmzzmzzzzz112zzTs2/12/1TsTsz上一页下一页返回 双线性变换法是否满足相匹配原理呢？ 首先，根据相匹配原理，要求当G(s)稳定时，G(z)也应该是稳定的。 将 代入（2.72）式，得 即有 （2.73） 由（2.73）式可知,双线性变换法将左半s平面映射到z平面的单位圆内。也就是说，如果原来系统的G(s)是稳定的，则通过双线性变换法得到的G(z)必然是稳定的。js)(21)(21jTjTz22222)2()21 ()2()21 (TTTTz上一页下一页返回 其次，双线性变

36、换后的稳态增益不变。 设连续系统的传递函数为 其稳态增益为bn/an 。若将G(s)进行双线性变换，有显然，G(z)的稳态增益仍然是bn/an。 nnnnmmnasasbsbsG11)(nnnnmmnazzTazzTbzzTbzG11)112()112()112()(上一页下一页返回uy ssUsYsG1)()()(11112)()()(zzTzUzYzG) 1()(2) 1()(kukuTkyky) 1()(2) 1(kykyTky上一页下一页返回 最后，双线性变换法具有一定的精度。 设有方程 （2.74）即 （2.75）将（2.71）式代入上式，得 （2.76）由此可得差分方程 （2.77

37、）这相当于应用AM2法（也称为梯形法）于（2.74）式 。 综上所述，双线性变换法符合相匹配原理，是一种具有一定计算精度的绝对稳定算法。 双线性变换法适合于以分子分母多项式形式描述的G(s)。 上一页下一页返回【例2.9】 已知连续系统的传递函数为 （2.78）试采用双线性变换法求出对应的脉冲传递函数和差分方程，并对所得结果进行分析。【解】 将（2.71）式代入上式，得脉冲传递函数 （2.79） 12)(2ssssG22)2()2() 1)(1(21)112(2)112(112)()()(TzTzzTzzTzzTzzTzUzYzG22122)22()22(211)2(2zTTzTTzTT上一页

38、下一页返回)2()()2(2)2()22() 1()22(2)(22kukuTTkyTTkyTTky122TT 于是,差分方程为 由（2.79）式可知，因为 所以，G(z)是稳定的。 G(s)的分子多项式为1阶，分母多项式为2阶，而G(z)的分子、分母多项式的阶次相同，均为2阶。 G(s)的稳态增益为0， G(z)的稳态增益也为0。 上一页下一页返回为了考虑双线性变换法所得仿真模型的精度，除了可以在时域中进行讨论外，也可以从频域的角度进行分析。事实上将，分别代入（2.78）式和（2.79）式，可得 （2.80） （2.81）取T=1s，分别求出（2.80）式和（2.81）式的幅频特性和相频特性

39、，如表2.5所示。 js Tjez1)(2)()(2jjjjG2232)1 ()(2j22)22() 1)(1()2(2)(TTeeeTTeGTjTjTjTj上一页下一页返回表2.5 G(s)与G(z)频率特性比较(rad/sec)0.81.01.21.5幅频特性连续系统模型0.099010.27520.44120.48780.500.49180.4615离散化模型0.09910.2770.4470.4930.4980.4760.417相频特性连续系统模型78.5856.6028.0712.680-10.39-22.62离散化模型78.5756.3626.519.57-5.0

40、6-17.67-33.56上一页下一页返回 设连续系统的传递函数为 （2.82）采用双线性变换法求取仿真模型的步骤如下： 将G(s)中的s按（2.71）式替换得到G(z)，并整理成有理分式形式 （2.83） 根据得到的G(z)，由Y(z)=G(z)U(z)，两边取Z反变换，得到便于计算机递推计算的差分方程。 （2.84） mnasasbsbsUsYsGnnnnmmn11)()()(nnnnzzzzzG111101)()() 1()()() 1()(101nkukukunkykykynn上一页下一页返回 双线性变换法是一种绝对稳定的算法，可采用大步长。 双线性变换法是由（2.70）式取一次项近似

41、得到的，所以对应的差分方程虽具有一定的精度，但当T取得太大时精度会降低。因此在选取T时，主要应考虑仿真精度和仿真速度的需求。经验表明，若取 ( 为被仿真系统的自然频率)，可以保证较好的仿真精度。 双线性变换法得到的差分方程的每项系数，都是预先一次计算确定的，在仿真递推过程中不需要重新计算，因此仿真递推的计算量较小。对于n阶连续系统，每步计算量约为(2n+1)次乘法。所以，该法适合于快速仿真。nT152n上一页下一页返回双线性变换法中由G(s)到G(z)的转换， 以及由G(z)到差分方程的转换，都可很容易地编程实现。双线性变换法具有串联性。如图2.29所示，G1(s)和G2(s)相串联。若令 G

42、(s)= G1(s). G2(s),则有 G(z)= G1(z). G2(z) 图2.29式中G(z) ,G1(z),G2(z) 分别是G(s),G1(s) ,G2(s)的双线性变换结果。这就是说，双线性变换没有破坏拉氏变换的串联性，可以用简单的低阶环节组成复杂的高阶系统。上一页下一页返回双线性变换公式不仅可以方便地应用于传递函数，也可以应用于状态空间模型。双线性变换后系统的阶次不变，且分子、分母多项式具有相同的阶次。将（2.71）式代入（2.82）式，有 （2.85）可见，G(z)的分子、分母多项式均为n次，在分子上增添了nm个z=-1的零点。由此得到的差分方程是多步关系式（当n1时）。由于

43、被仿真系统一般仅给出在初始时刻处nnnnmmnazzTazzTbzzTbzG11)112()112()112()(niimiimnnnpzzzzab11)()() 1(上一页下一页返回的系统的初始值，而在以前各采样点处的信息需要经过变换求得，这会给双线性变换法的应用带来一定的麻烦。可以证明，双线性变换后得到的脉冲传递函数频率特性与原连续系统的频率特性在低频段相差较小，在高频段差别较大。这一点也可以从表2.5看出。所以，双线性变换法主要用于有限带宽的系统。 上一页下一页返回 利用MATLAB的控制系统工具箱中的c2dm函数，很容易实现连续系统传递函数到离散系统脉冲传递函数之间的转换。【例2.10

44、】考虑如图2.30所示的单位反馈系统，试用双线性变换法对其进行仿真 (r(t)=1(t) 。上一页下一页返回图2.30 【解】采用如图2.31所示的参数化Simulink模型，其中Discrete Filter模块中参数Numerator设置为：nd1，Denominator为：dd1，Sample time为：0.05。其余两个Discrete Filter模块中的参数分别为nd2，dd2；nd3，dd3；Sample time均为0.05。 在MATLAB命令窗口运行下列指令T=0.05; n1=1 3.33; d1=1 24; n2=10; d2=1 0; n3=1; d3=1 10 1

45、6; nd1,dd1=c2dm(n1,d1,T,tustin);nd2,dd2=c2dm(n2,d2,T, tustin); nd3,dd3= c2dm(n3,d3,T,tustin);上一页下一页返回图2.31 (exam2_10.mdl)选择定步长的discrete（no continuous states）算法，步长为0.05，启动仿真，可以得到如图2.32所示的仿真结果。上一页下一页返回图2.321根匹配法原理 2根匹配法步骤 3根匹配法特性分析 上一页下一页返回根匹配法的基本思想是要使离散化模型的瞬态特性和稳态特性与原连续系统保持一致。更明确地说，就是要使离散化后所得脉冲传递函数的零

46、点和极点与原连续系统传递函数的零点和极点相匹配。所以，这种算法又称为零极点匹配法。先根据s域和z域之间的准确映射关系 ,直接由G(s)的零点和极点确定对应G(z)的零点和极点；再根据一定的原则确定G(z)的增益，则可以很方便地由G(s)求出与之对应的G(z)，使二者的瞬态特性和稳态特性保持一致。显然，根匹配法适合于以零极点形式描述的G(s) 。 Tsez 上一页下一页返回Tsez 求出连续系统的零点和极点，即将系统的传递函数写成下列的零极点形式 （2.87） 利用关系 在z平面上一一对应地确定零极点的位置。对于（2.87）式表示的系统，与z1,z2,zm相对应的零点为eTz1,eTz2,eTz

47、m;与p1,p2,pm相对应的极点为eTp1,eTp2,eTpm 。于是，有 （2.88） mnpspspszszszsKsUsYsGnm)()()()()()()(2121)()()()()(2121nmTpTpTpTzTzTzzezezezezezezKzG上一页下一页返回s0Tez)()()()()(2121nmTpTpTpmnTzTzTzzezezezzezezezKzG上一页下一页返回 在G(z)的分子上配上n-m个附加零点，使G(z)的分子多项式和分母多项式的阶次相同。因为当nm，即连续系统的分母多项式阶次高于分子多项式阶次时，在s平面的无穷远处实际上还存在n-m个零点。因此，在z

48、平面上必须再配上n-m个相应的零点。如果认为n-m个零点位于s平面上负实轴的无穷远处，即 ,那么z平面上相应的零点由可知，应配在原点处。于是，对应的脉冲传递函数为 （2.89） 根据G(s)的特征确定Kz，一般有两种方法：一种是选择适当的输入信号，应用终值定理，分别确定连续系统模型G(s)和离散化仿真模型G(z)的终值，然后根据终值（不等于零或无穷大）相等的原则确定G(z)的增益Kz。另一种是使由G(s)和G(z)分别得到的频率特性在某一临界频率处相同，从而确定Kz。一般对于低通滤波器应使其低频增益相同，而对于带通滤波器则应使中频增益相同。 根据得到的G(z)，由Y(z)=G(z)U(z)，两

49、边取Z反变换，得到便于计算机递推计算的差分方程。 上一页下一页返回)()()(1trtytyT)1(111)(111TsTsTsG1/)(TTzezzKzG上一页下一页返回 【例2.11】 设某连续系统的微分方程为 试用根匹配法确定其离散化模型。 【解】 首先写出系统的传递函数 它有一个极点p1=-1/T1和一个无穷远处的零点。 由（2.89）式可得脉冲传递函数 式中，T为采样周期。 根据拉氏变换的终值定理，有 而根据Z变换的终值定理，可得 令两终值相等，有 增益Kz为1111lim1)(lim)()(lim)(1000ssTssssGsRssGysss11/111111lim1)(1lim)

50、()(1lim)(TTzTTzzzzeKzzezzKzzzzzGzzzRzGzzy111/TTzeK1/1TTzeK上一页下一页返回 于是，求得的离散化模型为 根据可以进一步求出差分方程为 )()1 () 1()(11/krekyekyTTTT1/111111)1 ()(zeeezzezGTTTTTTTT上一页下一页返回102)()()(222nnnsssUsYsG)()(212pspssGn22, 11nnjp22211211TjTTpTjTTpnnnneeezeeez上一页下一页返回 【例2.12】 二阶系统的传递函数为试用根匹配法求其离散化模型。 【解】首先，将原系统的传递函数写成零极点

51、形式式中根据（2.86）式，得到z平面上的两个极点 原系统有两个无穷远处的零点，相应地应在z平面上配上两个原点处的零点，则对应的脉冲传递函数为根据终值定理分别求出单位阶跃输入作用下的终值为式中 TnTzTTjTjTzTjTTjTznnnnnnnnnnezTezzKezeeezzKeezeezzKzG222221122112)1cos(2)()()(2222baKzzzGzzysssGyzzs11)(1lim)(11)(lim)(10TnTnnebTea22)1cos(2上一页下一页返回11baKzbaKz1212211)1 ()(bzazbabazzzbazG)()1 ()2() 1()(ku

52、bakbykayky令两终值相等，得于是 最后，得到离散化模型为对应的差分方程为上一页下一页返回2) 1()()()(sssUsYsG2)() 1()(TzezzzKzG22) 1()(,1)(zTzzUssU上一页下一页返回Tsez 【例2.13】 连续系统的传递函数为 试用根匹配法求其离散化模型。 【解】根据 以及G(s)的零极点，可得 在利用终值定理确定Kz时，由于在s平面上存在一个原点处的零点，加上跃阶信号后，终值为零。为了使终值取得不为零的有限值，应当加上单位斜坡信号，即于是，有 从而得到最后，离散化模型为11) 1(lim1)(lim)(22020sssssssGyss222121

53、)1 () 1()() 1(1lim) 1()(1lim)(TzTzzzeTKzTzezzzKzzzTzzGzzyTeKTz2)1 (211222)1 (1)1 ()() 1()1 ()(zezTeezzzTezGTTTT上一页下一页返回由根匹配法所采用映射关系 可以看出，如果原系统模型G(s)是稳定的，则它的全部极点都位于左半s平面。根匹配法具有一定的精度。 以【例2.13】为例，比较G(s)和G(z)的频率特性。 取采样周期T=1s，由将G(s)和G(z)的幅频特性和相频特性列于表2.6。Tsez 上一页下一页返回2121)() 1()1 ()(ezzzezG2)3679. 0() 1(4

54、 . 0zzz （rad/sec）0.81.01.22.0幅频特性 连续系统模型0.09900.27520.44120.48780.50.49180.4 离散化模型0.09540.26590.43070.48120.50.50010.449相频特性 连续系统模型75.5856.628.0712.680-10.39-36.87 离散化模型80.5062.3739.627.919.1012.450.35表2.6 G(s)与G(z)频率特性比较上一页下一页返回 比较表2.6中G(s)和G(z)的频率特性可知，二者的幅频特性是比较接近的，而离散化模型G(z)在相位上超前较大。当=1时

55、，G(s)的相位为0，而G(z)的相位为19.1。一般而言，如果T 取得比较大，计算精度还会降低。 上一页下一页返回实时数字仿真通常是指把一个数字仿真过程嵌入到一个具有实物模型的实际系统或仿真系统的运行过程中，这就要求仿真模型的时间比例尺必须完全等于原始模型的时间比例尺。当仿真系统有实际装置或操作人员介入试验回路时，为了真实反映实际系统的工作情况，要求仿真系统必须按照实际系统运行的时序要求来完成数字仿真过程的各个步骤。 上一页下一页返回3.5.1 实时数字仿真的概念 3.5.2 实时数字仿真算法的特性 3.5.3 常用的实时数字仿真算法 上一页下一页返回【例3.14】 飞行器控制系统的半实物仿

56、真 一个飞行器控制系统可以分为被控对象和控制器两个子系统。在导弹、运载火箭的设计和飞行试验前的仿真中，控制对象的运动通常为弹体运动，包括质心运动，绕质心的旋转运动和弹体振动等。显然，在实验室里难以直接实现被控对象的运动。所以，通常采用一定的数值方法，将描述弹体运动的数学模型进行离散化，建立数字仿真模型，然后编程运行来代替弹体运动，并将其称为数字处理过程。而控制器及其执行机构，包括舵式控制发动机等则用实物表示，并将其称为实物系统过程。于是，由实物系统过程和数字处理过程等组成仿真回路处理过程等组成仿真上一页下一页返回回路，即半实物数字仿真模型。由于在这样的飞行器控制系统的仿真中介入了实物，因而在弹

57、体运动的数字处理过程中信息的采样输入、仿真模型的计算、数据的输出都必须在规定的时间内按照实际运行的时序完成，进行实时数字仿真。 可以将控制系统的实时数字仿真描述如下。 假设一个实际的控制系统是由实物系统过程A和实物系统过程B两部分组成，如图2.33所示。而在实时仿真过程中，常常将系统的一部分（例如，实物系统过程A）用数学模型代替。在实际实现时，用计算机的数字处理过程A来代替，如图2.34所示。 上一页下一页返回上一页下一页返回图2.33 图2.34图2.34中计算机数字处理过程A的输入是z(t)，输出是y(t)。对于任意一组输入值，经过采样控制和A/D转换，由计算机数字处理过程A的计算得到相应

58、的输出值ym。将ym称为数字处理过程A对于输入zm的响应。响应ym对于输入zm的时间延迟称为响应时间。过程B为实物系统，通过D/A转换和输出控制接受过程A的输出量ym，并输出z(t)。计算机数字处理过程A必须在与实物系统同步的条件下获取动态输入信号，并实时地产生动态输出响应。这时，仿真模型的输入和输出都是具有固定采样周期的数值序列。在数字处理过程A满足系统各项功能要求的情况下，对于任意特定的输入z(t)，响应时间都能满足系统所要求的时间限制，则这种数字仿真过程通常为实时数字仿真过程。 上一页下一页返回 实时仿真模型中包含有实物或人，它具有如下特性。 1实时性 实时性主要体现在如下几个方面。（1

59、）实物模型的运行时间 这是实物模型在输入作用下的运行时间。若用计算机来代替，则响应时间的要求就是实时要求，其具体数值必须满足随机尖峰负载时的处理要求。（2）数据转换（如数模D/A、模数A/D转换）的时间 从发出转换命令到完成的时间。D/A、A/D转换本身需要一定的时间。在仿真中，希望越快越好，但转换速度越快，系统成本越高。因此要根据实际需要确定。上一页下一页返回 （3）输入/输出响应时间 计算机接收动态输入到产生动态输出的响应时间，有固定的最小的要求。它是计算过程的延时，可以通过采用不同的算法和计算步长来达到要求。如果无法达到要求，则必须更换仿真系统所用的计算机。 2周期性 在实时仿真模型中，

60、模型和各个子模型都有规律性和周期要求。它们以设定的周期时间接受输入信息并一帧一帧地产生输出信息。每一个子系统都必须按照规定的顺序，在分配给它们的时间片内完成信息采样、计算、恢复和输出等任务。上一页下一页返回 3可靠性 由于有实物包含在内，与一般计算机仿真系统相比，实时仿真系统要复杂得多，并且一般均为专用的仿真系统。因此，仿真系统的可靠性是第一位的。在【例2.14】的飞行器半实物仿真中，系统运行的任何差错都可能导致仿真的失败。而引起差错的原因，一是数学模型或仿真模型或相应的计算机软件有错；二是外界干扰破坏了程序的正常运行；三是硬件系统故障。第一类错误可以通过修改模型、选择恰当的仿真算法或调整软件