第四章 快速数字仿真

1、 第四章第四章 快速数字仿真快速数字仿真 前面介绍了应用数值积分法和离散相似法进行控制系统数字仿真的基本原理、方法和程序，对于控制系统的数字仿真是比较方便的。但是，这些数字仿真程序的一个共同特点之一是：为使仿真达到一定的精度，计算步长不能偏大，因此计算量加大，计算速度较低，这对要求进行实时仿真(在线仿真)的场合无法满足要求，就是离线仿真也嫌太慢。如在对一些复杂系统或高阶系统进行数字仿真时，往往仿真实际运行时间需要数小时甚至几十小时，速度太慢而限制了数字仿真的推广和应用。 如：在液压系统的数字仿真中，由于参数的限制，计算步长要取的极小(如 h5E6)，这将使数字仿真时间很长，尤其是如果系统中的两

2、个环节的时间常数相差达到几十倍或上百倍时，如按小的时间常数来确定计算步长，将会引起计算不稳定以致整个系统不稳定,仿真无法进行。 q为了解决数字仿真速度问题，一般采取一下三种措施： 硬件措施：采用高速计算机 软件措施：采用执行速度快的语言，如汇编语言. 算法措施：研究并采用适合于快速仿真的仿真算法. q本节只讲述快速仿真算法。一般情况下，在解决计算速度与精度这一矛盾时，以计算速度作为矛盾的主要方面，在满足计算稳定性及工程要求的精度条件下，尽可能提高仿真计算速度，降低不必要的精度以满足速度要求。一般对快速仿真算法有两点基本要求： 每一步的计算量要少； 算法要具有良好的计算稳定性，在保证一定的仿真精

3、 度情况下，允许采用较大的计算步长。 前面我们介绍的各种仿真方法基本上都是基于时域的范围内对系统进行研究的，也就是说系统的输入、输出、以及状态量等主要是在一种连续状态下的时间函数（在研究过程中主要以数值积分或虚拟的离散状态分析为主），但在研究过程中，这些方法在各方面多少都有局限性。频域仿真建模方法：面向频域仿真建模方法：面向 S 域的传递函数域的传递函数G (s)根据相似原根据相似原 理，得到与它理，得到与它相匹配相匹配的的 Z 域传递函数域传递函数G(z)。若若G (s)是稳定是稳定的，而脉冲传递函数的，而脉冲传递函数G(z)也应是稳定的，对于同一输入信也应是稳定的，对于同一输入信号，由号，

4、由G(z)求得的求得的y(kT)输出与由输出与由G(s)求出的求出的y(t)应具有相一应具有相一致的时域特性或具有相一致的频率特性。致的时域特性或具有相一致的频率特性。 第一节第一节 替换法替换法快速仿真解决的主要矛盾是快速性，而对精度可以降低一些要求。因此，对于一个高阶连续系统若能将其传递函数G (s)比较方便的转换为与其相匹配的，并允许采样周期 T 较大的脉冲传递函数G (z) ，然后由G (z)获得差分方程进行数字仿真，不但能够减少计算步数，也能使每一步的计算量大为减少，从而获得快速仿真的效果。 v本节介绍的替换法是根据相匹配原理由直接导出的变换方法。 一、替换法基本原理 替换法的基本思

5、想替换法的基本思想是，设法找到s域（连续域）与z域（离散域）的某种对应关系，然后将G(s)中的变量s转换成变量z，由此得到与连续系统传递函数G(s)相对应的离散系统的脉冲传递函数G(z) ，进而获得进行数字仿真的递推算式，以便在计算机上求解计算。 根据拉普拉斯变换与z变换的关系，可得： （4.1-1） 或 式中，T采样周期 式（4.1-1）是一个超越函数，不能直接用它来替换 sTez zTsln11.欧拉替换 )(tfdtdx已知微分方程：kkkkkxTxTfxx1根据欧拉法公式：xTxz1所以可得：sxx1因为1 zTxx即：11zTs故有：欧拉替换：简单，但是欧拉替换：简单，但是稳定性差稳

6、定性差，并不实用。，并不实用。下面分析其稳定性：下面分析其稳定性： )1()1( 11 1 z 1 则 1Tsz 设22222222222TTTTzTTzjs也就是：故有平面上的单位圆对于因为，Z平面上的单位圆按该替换式反映射到平面上的单位圆按该替换式反映射到S平面上，将是一个平面上，将是一个 以（以（1/T，0）为圆心，以为圆心，以1/T为半径的圆为半径的圆。 图图4.1 4.1 S S域到域到Z Z域的映射关系域的映射关系l 一个一个原来稳定的系统原来稳定的系统G(s)，通过替换得到的仿真通过替换得到的仿真模型模型G(z)却却可能是不稳定的可能是不稳定的。 2. 双线性变换双线性变换112

7、 kkkkxxTxx观察梯形积分公式：112 121 zzTsxzTxz即：可得：）（一种近似映射关系：平面的平面到，则可得由若只取式中第一项近似为：的一种幂级数展开形式或由：2-4.1 112 )11(121)11(3111 2lnlnln1 123zzTszszznzzzzzzzTsn）（也可以写作：3-4.1 )1z-1(2 1-1zTs 式（4.1-2）和式（4.1-3）都称为算子的屠斯汀公式。用这一关系替换G(s)G(s)中的s s的方法称为双线性变换法或屠斯汀法。可以证明，双线性变换法符合相匹配原理。 ）（）得：代入式（将4-4.1 )2()2(1)2()2(1 3-4.12222

8、2TTTTzjs。，则；若则，；若，则，即）可知，若由式（1010104-4.1zzzjs 这就是说，这就是说，Z平面上的单位圆，映射到平面上的单位圆，映射到S平面上将是整个左平面上将是整个左 半平面，其逆也真。半平面，其逆也真。 即即如果原来如果原来G(s)稳定，那么稳定，那么G(z)也是稳定的也是稳定的。二、替换法步骤采用双线性变换法求仿真模型的具体步骤如下： 设已知一线性系统的传递函数为： ）（5-4.1 )( )(1111110nmasasasbsbsbsbsGnnnnmmmm）（的形式：，并整理成有理分式）替代得到按式（中的将6-4.1 1)( )(3-4.1)() 1 ()1(12

9、211)1(1110nnnnnnnnzzzzzzzzGzGssG)0() 1 () 1()( )0() 1 ()2() 1()( )()()()2(110121uununuyynynynyzUzGzYznnnn方程形式：的差分转换成便于计算机计算变换，将利用逆。为输出序列，为输入序列，式中，nkkyku,2, 1 ,0)()( 。真模型。计算步长取利用双线性变换求其仿）（函数为：已知一线性系统的传递例sssssG1T 7-4.1 12 . 01)( 1-4.124 . 566 . 432 )4 . 04()82()4 . 04()2(2)2( 1)112(2 . 0)112(1112)( 7-

10、4.13-4.11 121221222212221121111zzzzTTzTzTTTTzTzTTzzTzzTzzTzGs可得：）中的）替换式（）由式（解：)2() 1(2)(3185. 0)2(852. 0) 1(111. 1)( )()()()2(nunununynynyzUzGzY转换成差分方程：将根匹配法也是一种从系统传递函数G (s)导出脉冲传递函数G(z) 方法，是根据相匹配原理提出的又一种常见的快速仿真方法。一、根匹配法原理根匹配的基本思想基本思想是使离散化所得的脉冲传递函数与连续系统的传递函数的零点和极点相匹配，从而使离散化模型的动态、静态特性与连续系统相同。所以，这种方法又称

11、零极点零极点匹配法匹配法。连续系统的动态、静态特性由传递函数G(s)G(s)的零、极点分布及增益大小所决定，采样系统的动态、静态特性由脉冲传递函数G(z)G(z)的零、极点分布及增益大小所决定。如果能找到一种简单的映射关系z=f(s)z=f(s)，将连续系统G(s)G(s)的零、极点映射为脉冲传递函数G(z)G(z)的零、极点，再根据一定的原则确定G(z)G(z)的增益，则可方便地由G(s)G(s)求得与之等价的G(z)G(z) 。 第二节第二节 根匹配法根匹配法匹配原理。的原理，显然它符合相。这就是根匹配法等价的，从而得到与的增益确定的其他特征然后根据的零、极点确定相对应的点的零、极这一映射

12、关系直接由用根据根匹配的思想，利)()()()(), 2 , 1, 2 , 1,()(), 2 , 1, 2 , 1,()(zGsGKzGsGnimjezezzGnimjrsGezzTiTrjijsTij为系统增益。极点，个为系统的，个零点，为系统的，式中，）（KnmrrrmnsssrsrsrsKsUsYsGnmnm212121211-4.2 )()()()()()()( 二、根匹配法步骤二、根匹配法步骤 由根匹配法求仿真模型的具体步骤如下： 设原系统的传递函数为： ）（）的零、极点，得到：替换式（利用2-4.2 )()()()()( 1-4.2,) 1 (2121mnezezezezezez

13、KzGezezTTTTrTrTrzTiTrjnmij）（冲传递函数为：应配在原点处，于是脉知，平面上相应的零点由，那么无穷远处，即平面负实轴的为零点位于个相应的零点。如果认须再配上平面必个零点。因此，在存在着平面无穷远处实际上还在次高于分子时，即连续系统的分母阶）中，当在式（3-4.2 )()()()()( 01-4.22121mnezezezzezezezKzGezzssmnzmnsmnTTTmnTrTrTrzTnm使中频增益相同。而对于带通滤波器则应应使其低频增益相同，。一般对于低通滤波器从而可确定某一临界频率处相同，分别得到的频率特性在与使入函数。，则应另选输。若，从而确定相等分别得到的

14、输出变量与令由选择适当的输入函数种方法：，一般有两的原则确定的终值，根据终值相等散系统的终值及离值定理求出连续系统在典型输入下，根据终 )()(B. 0)()0)()( )()(,A. )( )()2(zzzKzGsGyKyyzGsGKzGsG差分方程形式。的转换成便于计算机计算变换，将利用逆 )()()()3(zUzGzYz差分方程。试用根匹配法求系统的）（数为：设一连续系统的传递函例 4-4.2 )( 1-4.2asKsG为采样周期。式中，）（可得脉冲传递函数：），故由式（分子阶次，分母阶次，一个无穷远处零点，有一个极点零、极点确定）根据（解：TezzKzGmnasGzGsGaTz 5-4

15、.2 )( 3-4.20 1 )( )()(1 1）（所以，应该相等，故得：上述两个终值）（变换的终值定理可得：根据）（理知：根据拉氏变换的终值定。的终值，进而确定和）用终值定理分别求（8-4.2 )1 ( 1 )( 111lim 7-4.2 )()(1lim)( 6-4.2 1lim )()(lim)( )()(21100aTzaTzaTzaTzzzsseaKKeKaKyeKzzezzKzzzUzGzzyzaKsasKssUssGyKzGsG）（）得：）代入式（将式（9-4.2 )()()1 ()( 5-4.28-4.2 zUzYezzeaKzGaTaT）（变换，可得差分方程）进行逆将式（求

16、差分方程。）根据（10-4.2 )()1 () 1()( 9-4.2 )(3nueaKnyenyzzGaTaT差分方程。试用根匹配法求系统的）（数为：一个连续系统的传递函例 11-4.2 1)()( 2-4.22sssG为采样周期。式中，）（）有根据式（零、极点确定）根据（解：TezzzKzGzGsGTz 12-2 . 4 )() 1()( 3-4.2 )()(1 2。来确定斜坡输入所以选择入时，输出终值为零，有零值零点，加阶跃输因。的终值，进而确定和）用终值定理分别求（zKttusGKzGsG)( )( )()(2 13-4.2 )1 ( )( )1 ( ) 1()() 1(1lim )()

17、(1lim)( 11) 1(lim )()(lim)( 2222112200）（相等，就要求：为使终值根据终值定理知：TeKyeTKzTzezzzKzzzUzGzzysssssUssGyTzTzTzzzss）（变换，可得差分方程）进行逆将式（求差分方程。）根据（15-4.2 )1()()1 ()2() 1(2)( 14-4.2 )(322nunuTenyenyenyzzGTTT）（）得：）代入式（将式（14-4.2 )()()() 1()1 ()( 5-4.213-4.2 22zUzYezTzzezGTT的差分方程。配法求系统用双线性变换法及根匹 ) 115. 0)(1)(110() 15(1

18、00)( 4.1sssssG个点即可）。，手算输出量初值位阶跃下的仿真结果（秒，计算在单次，而采样频率为，若的差分方程。用根匹配法求出50 /10rad/s 154 . 0 2)( 4.20222ysssGnnnn习习 题题）（其解为：）（矩阵方程）中状态方程变为齐次时，式（当输入函数）（它的时域解为：）（间表达式为：出的连续系统的状态空假定一个单输入单输4.3.4 )0()( 4.3.3 4.4.10)(4.3.2 )()0()( 4.3.1 )0()()()()()( 0)(0 xetxAxxtudBuexetxxxtCxtytButAxtxAtttAAt一一.增广矩阵法的基本思想增广矩阵

19、法的基本思想第三节第三节 增广矩阵法增广矩阵法）（展开成级数形式：将）（时，则计算步长选为。当）就可以递推计算出可以求得，则根据式（如果4.3.6 ! 3! 2 4.3.5 )() 1( )(4.3.43322nTATATAATIeenTxeTnxTtxennAtAtAtAt）（）（）变为：则式（龙格库塔法相同。项，则计算精度与四阶可以证明：如果取前4.3.7 )(2462) 1(4.3.55443322nTxTATATAATITnx由于 都可以在仿真前算出，所以式(4.3.7)右端的系数项可以事先求出，而仿真计算就变成每次只作一个十分简单的递推计算，从而达到快速仿真。 然而一般系统的状态方程

20、都是非齐次的，按式(4.3.2)求解时，除了计算一个自由项外，还要计算一个强制项 使计算复杂化，所以应设法将非齐次项，即输入变量“增广”成为新的状态变量(即增广状态变量)，使原状态方程变成增广形式的齐次状态方程：432AAAA，dBuettA)(0)(0)0(4.3.9 )()(4.3.8 )(xxtxCtytxAx）（）（式(4.3.8)，称为增广状态方程，式(4.3.9)称为增广输出方程，其中 为增广状态向量，它包含了 及增广进去的增广状态变量。 为增广状态矩阵， 为增广输出矩阵。对于式(4.3.8)，就可以用式(4.3.7)进行仿真计算，这种方法就称为增广矩阵法。这时对于n阶系统就是把输

21、入 当作第n+1个状态变量，使原来解非齐次方程的问题变成求解齐次方程的问题。现在的关键问题是如何把 增广成新的状态变量构成新的增广状态向量，下面将分别就几种典型输入信号来说明增广矩阵的形成方法。 )(tx)(txAC)(tu)(tu）（初始条件）（）（程和输出方程为：故可得增广后的状态方则个状态变量为：为阶跃幅值，定义第，设阶跃输入信号4.3.12 )0()0( 4.3.11 )()( 0)( 4.3.10 )()(00)()( )0( 0)( 1)( 1 1)( ) 1 (0011110110100uxxxtxtxCtytxtxBAtxtxuxtxtutuxnututunnnnnnn二二.

22、典型输入信号的增广矩阵法典型输入信号的增广矩阵法）（程为：故可得增广后的状态方则个状态变量为：个，第，定义第设斜坡输入信号4.3.13 )()()(0001000)()()( )0( 0)0( 0)( )()( )()( 21)( )2(21210212012010txtxtxBAtxtxtxuxxtxutxtxtututxnntutunnnnnnnnnn）（初始条件）（输出方程为：4.3.15 0)0()0()0( 4.3.14 )()()(00)( 002121uxxxxtxtxtxCtynnnn 时域矩阵法是在时域中用无穷矩阵来进行系统分析的方法。这种方法不仅适用于分析采样控制系统，而且

23、可以通过引入适当的采样开关和保持器以分析连续系统。它的好处是，如果能获得系统的脉冲响应，那么采用该法可以允许采用较大的计算步长，每一步的计算量少且与系统的阶次无关，从而实现快速数字仿真。 一、时域矩阵法一、时域矩阵法1. 时域矩阵的概念时域矩阵的概念 设有如图4.4.1所示的开环采样系统，受控对象的传递函数为 ，保持器的传递函数为 ， 则系统的开环传递函数 。设初试条件为零时，对应 的脉冲传递函数为 。 )(0sG)(sGh)()()(0sGsGsGh)(0sG)()(1sGLtg第四节第四节 时域矩阵法时域矩阵法) 1 . 4 . 4( )()()(0knnunkgky的输出为：根据卷积和定

24、理，系统)(0sG)(sGh)(tu)(*tuT)(ty开环采样系统图 4.4.1)()0() 1 () 1()0()()() 1 ()0()0() 1 () 1 ()0()0()0(4.4.1kugukgukgkyugugyugy）展开得：将式（）（写成矩阵形式：4.4.2 )() 1 ()0()0()2() 1()(00)0() 1 (000)0()() 1 ()0(kuuugkgkgkggggkyyy)0()2() 1()(00)0() 1 (000)0( )()2() 1 ()0( )()2() 1 ()0( 4.4.3 gkgkgkggggGkuuuuUkyyyyYGUYTT其中：）

25、（或记作：q列矢量Y和U包含了输出和输入的各个离散值，即各采样时刻值。G是由(k+1)(k +1)个元素组成的方阵，其元素由g(t)在离散时刻的值确定，这个矩阵是下三角阵，只在对角线上、对角线下半部有值，其余元素为零。若响应 y(k)在各采样时刻都要被确定，则G的阶次是无穷大的，所以称为无穷矩阵，又称时域矩阵。实用时，一般 G的阶次有十几到二十几左右就能满足要求。 q从式(4.4.3)看出，若能在系统仿真前先测得u(k)和g(k)(g(k)也可通过拉普拉斯反变换求得)，以一定的数字形式送入计算机，在仿真时就只需计算一个矩阵相乘式，且这种算法与系统的阶次无关。因此这种算法允许采用较大的计算步长，从而提高了仿真计算速度。 2. 闭环连续系统的时域矩阵闭环连续系统的时域矩阵 闭环系统离散化后的方框图如图4.4.2所示。 框图闭环系统离散化后的方图 4.4.2)(0sG)(sGh)(tu)(*teT)(ty)(te) 5 . 4 . 4( )4 . 4 . 4( 4.4.2YUEGEY）可知：由图（）（）（写成矩阵形式：4.4.7 )() 1 ()0()() 1 ()0()() 1 ()0(4.4.6 )() 1 ()0()0()2(