第五章 快速数字仿真法

1、第五章第五章 快速数字仿真法快速数字仿真法n前两章从数值积分和面向结构图的仿真两方面讨论了控制系统数字仿真的基本原理方法、方法和程序。这些方法用于对控制系统进行非实时的仿真研究有很大方便，尤其是在具备一些通用的仿真程序时，这些方法就显得更为方便。但在一般情况下，为了达到一定的计算精度，这些方法的计算量比较大，因此计算的速度受到一定的限制，往往在实际应用中，还不能满足实时仿真的要求。有必要寻找一些能加快仿真速度的方法。能解决这个问题的方法很多，但各有各的特点和局限性。第五章第五章 快速数字仿真法快速数字仿真法 本章第一节到第三节介绍了三种常用的方法，以便根据情况选用。这些方法不仅在连续系统仿真时

2、可以使用，也可以用于连续控制器在计算机上的离散实现。本章还在第四节介绍了计算机控制系统的仿真。 5.1 增广矩阵法 5.2 替换法 5.3 零极点匹配法 5.4 计算机控制系统仿真5.1增广矩阵法 一、基本思想n假定一个连续系统的状态方程为 (5.1.1)n这是一个齐次方程，它的解是 (5.1.2)n已知n可以证明：如果取前五项，则计算精度与四阶龙格-库塔法相同。这就是说，如果被仿真的系统是一个齐次方程，在选定计算步距为n以后，若只取 的前五项，则有 axtx)( )0()(xetxatnatatntataatie)(!1! 3! 23322ate5.1增广矩阵法(5.1.3)n由于 都可以在

3、仿真前计算出，所以(5.1.3)式所示的递推计算公式中右端的系数项可以事先求出，而仿真计算就变成每次只做一个十分简单的递推推算。n但是，实际的物理系统模型大多是一个非齐次方程，即 (5.1.4)n其中为系统的控制量，假定它是一个单输入系统。根据控制理论可知，(5.1.4)式的解为n t0 (5.1.5) ) 1()2462()(443322tnxtatataatintx432,aaaabuaxxdbuexetxttaat)()0()(0)(5.1增广矩阵法n显见，求解(5.1.5)式的非齐次方程时，它的解，除了一个自由项之外，还有一个强制项： 。由于 的任意性，(5.1.5)式一般不容易求解。

4、但是，对于某些特殊的输入函数，如果能将控制量 增广到状态变量中去，使(5.1.5)式这样的非齐次方程，变成一个齐次方程(5.1.1)式，就可以避免计算复杂的强制项，而利用类似(5.1.3)式的计算方法。二、典型输入函数时的增广矩阵n假定被仿真的系统为(5.1.6)dbuetta)(0)(u)(tubuaxtx)( 0)0(xx)()(tcxty5.1增广矩阵法n其中a为维矩阵，即表示有n个状态变量。(1) 阶跃输入时n设n则定义第n+1个状态变量为n故n可得增广后的状态方程即输出方程为)( 1)(0tutu)( 1)()(01tututxn0)(1txn01)0(uxn5.1增广矩阵法(2)

5、斜坡输入n设n 则定义)()(00)()(11txtxbatxtxnn001)0()(uxxtxn)()(0)(1txtxctyntutu0)(tututxn01)()(012)()(utxtxnn0)0(1nx02)0(uxn5.1增广矩阵法n因此，系统增广后的状态方程为n初始条件为)()()(0001000)()()(2121txtxtxbatxtxtxnnnntnntxtxtxcty)()()(00)(2100210)0()0()0(uxxxxnn5.1增广矩阵法(3) 指数输入n设 n 定义 n 则有 n故系统增广后的状态方程为 teutu0)(tneutx01)()()(101txe

6、utxntn01)0(uxn)()(10)()(11txtxbatxtxnn001)0()(uxxtxn)()(0)(1txtxctyn5.2替换法 n一个连续物理系统最常见的数学表现形式就是s域的传递函数 。替换法的基本思想就是，设法找到s域（连续域）的某种对应关系，然后将 中的变量s转化成变量z，由此得到与系统传递函数 相对应的离散系统脉冲传递函数 （脉冲传递函数即采样系统输出脉冲序列的z变换与输入脉冲序列z变换之比），进而获得进行数字仿真用的递推算式，以便在计算机上求解计算。)(sg)(sg)(sg)(sg5.2替换法 n5.2.1 简单替换法ns域与z域的基本关系是 ，式中t为采样周期

7、，即 仿真计算时的计算步长；或者 ，这是一个超越 函数，不可直接用来进行替换。实际上必须要寻找其他近似的表示关系。一种最简单的替换关系可以从一阶差分方程中得到。n用传递函数 表示系统，在时域内可以用一个微分方程来表示，例如，系统(5.2.1) stez ztsln1)(sgssgsxsy1)()()(5.2替换法 的时域表示为 (5.2.2)n假若导数计算用下述差分来近似：n或间写成 (5.2.3)n则微分方程(5.2.2)式即等价为下述差分方程 ：)()(txdttdyttkyktydttdyktt) 1()()(tkykydttdyktt) 1()()(5.2替换法 或 (5.2.4)n当

8、微分方程中导数的计算采用(5.2.3)式差分表达式时，称为向后差分法。当然，导数的近似计算还可以采用向前差分法，此时，导数按下述差分来近似计算：n方程(5.2.2)式可等价为下述差分方程：)() 1()(kxtkyky)() 1()(ktxkykytkykydttdyktt)() 1()()()() 1(kxtkyky5.2替换法 或 (5.2.5)n现对差分方程(5.2.4)进行z变换，则是 (5.2.6)n比较(5.2.6)与(5.2.1)式，有 (5.2.7)或 (5.2.8) )()() 1(ktxkyky1)()()(ztzzgzxzy11ztzstzzs15.2替换法 n关系式(5

9、.2.7)或(5.2.8)就是一种最简单的替换方法。它表明，如果将传递函数 中的用(5.2.8)式替换，就可以将s域中的传递函数 变换为z域中脉冲传递函数 ，从而可以得到递推算式(5.2.4)式。这种替换相当于向后差分法。从(5.2.4)式可见，这种向后差分法也就是数值积分法中的超前欧拉法。n 若对(5.2.5)式做z变换，则可得要下述脉冲传递函数： (5.2.9)将(5.2.9)式与(5.2.1)式比较，有 )(sg)(zg1)()()(ztzgzxzy)(sg5.2替换法 (5.2.10)或 (5.2.11)n 关系式(5.2.10)或(5.2.11)也是一种最简单的替换方式，这种替换式相

10、当于向前差分法，即数值积分中的欧拉法。n 上述两种替换法比较简单，但局限性很大，实际工程中很少采用。下边对两种变换作简单的讨论。实际上替换式(5.2.8)及(5.2.11)均可以看做是s平面与z平面之间的相互映射。对(5.2.11)式来说，由于 ，设 ，所以11ztstzs11 tszjs5.2替换法 n对于z平面上的单位圆(及稳定域)， ，故n即等于 (5.2.12)n因此将z域中的单位圆映射到s域平面上，正好是以 为圆心，以 为半径的一个圆，如图5.1.1(a)所示。这就 是说z平面上的单位圆，按(5.2.11)式变换，它将是s平面上以为 半径的一个圆。 2222)1 (ttz12z1)1

11、 (222tt2222)1()1(tt)0 ,1(tt1t15.2替换法 n反过来说，s平面上只有部分面积通过(5.2.11)式才能映射到z平面的单位圆内。显然，若一个系统 是稳定的，其极点分布如图5.1.1(a)所示，通过(5.2.11)变化后有 两个极点能映射到在z平面的单位圆内；而极点 则映射到z平面的单位圆外。这样，原来为稳定的系统，通过(5.2.11)式的替换，仿真模型变得不稳定了，因此这种仿真模型失真太大。由(5.2.12)式还可以看出，若要使 稳定，就要求增大半径 ，即减 小计算步长t，从而增加了计算工作量，故采用(5.2.11)式不适合快速数字仿真。 )(sg21,ss3s)(

12、zgt15.2替换法 图5.2.1 简单替换法的映射关系 5.2替换法 由(5.2.8)式的替换关系，可以推得，s平面的左半部将映射到z平面单位圆的局部范围内，如图5.1.1(b)所示。为说明这一点，(5.2.8)可以改写为 n所以 (5.2.13)n显然，若 ， ，则 ；若 ，则(5.2.13)右端分母大于分子，即 ； tstststsz1121212111211122222)()1 ()()1 (4121ttttz0js 21212z021212z5.2替换法 n若 ，则(5.2.13)右端分母小于分子，即 。这表明，s平面右半部全部映射在该圆之内。由于这个小圆位于z平面的单位圆内，所以原

13、连续系统是稳定的，利用(5.2.8)式的替换关系得到的数学模型也一定是稳定的，并且与计算步长t无关。但是，由于整个s平面的左半部不是映射到z平面的整个单位圆内，而只是映射到其中一个小圆内，显然与准确z变换的结果相差较大，从而使仿真模型失真较为严重。 n例5.2.1 现给定一个二阶系统的传递函数，使用两种替换，求仿真数学模型。计算步长为t=1s。021212z5.2替换法 n解：解：(1)利用前差替换公式(5.2.11)式，得n对该式做z反比变换 n 或 12 . 01)()()(2sssxsysg)2 . 01 ()2 . 02(1) 1( 2 . 0) 1(1)()()()(222221tt

14、ztzttztzsgzxzyzgtzs)()()2 . 01 () 1()2 . 02()2(22kxtkyttkytky)()(8 . 1) 1(8 . 1kxkyky) 1() 1(8 . 1)(8 . 1) 1(kxkykyky5.2替换法 (2)利用后差替换公式(5.2.8)式，得nz反变换后得 1) 1(2 . 0) 1(1)()()(222tzzztzzxzyzg)2 . 01 (1)2 . 01 ()22 . 0()2 . 01 (222222ttzttzzzttt)()2 . 01 (1) 1()2 . 01 ()2 . 02()2(22kyttkytttky)2()2 . 0

15、1 (22kxttt)2(455. 0)(455. 0) 1(kxkyky5.2替换法 n从上述计算结果可见，同一结果采用不同的计算公式，所得到的数学模型差别是很大的，不仅差分方程系数不同，而且形式也常常是不同的。5.2.2 双线性变换双线性变换 n 上述简单替换法使用起来问题较多，比较适用的替换法是双线性变换法（又称tustin法）。n这种替换关系可以从s到z变量准确的映射关系推得。n根据定义，z变量与s变量的关系是或者 (5.2.14)stez ztsln15.2替换法 n(5.2.14)式可以展成无穷级数 n如果取该级数的第一项作为它的近似，则 (5.2.15) 或 (5.2.16)n数

16、学上称这种变换为双线性变换。利用(5.2.15)就可以把传递函 转换为脉冲传递函数 。33) 1() 1(31112lnzzzzz112zzts2/12/1tstsz) sg) zg5.2替换法 n应当指出，这种变换关系相当于数值积分中的梯形法。事实上，积分环节(5.2.1)式对应的微分方程形式，若采用梯形积分法可得 (5.2.17)n对(5.2.17)式的差分方程进行z变换得n 可见，积分环节从梯形积分法所推得的递推公式与从双线性变换所得的递推公式是相同的。)1()(2) 1()(kxkxtkykyszztzgzxzy1) 1(2) 1()()()(5.2替换法 n采用双线性变换的一个和重要

17、的优点是，在s域里稳定的传递函数，通过这种变换得到的脉冲传递函数 也一定是稳定的。因为这种变换将s域左半平面准确地映射到z域的单位圆内。实际上， 将代入关系式(5.2.16)式中，得22222221221ttttz)(zgjs5.2替换法 n由该式可知，若 ，即 ，则 ；若 ，则 ；若 ，则 。这就是说，采用(5.2.16)式映射，s域左半平面映射到z域的单位圆内，如图5.2.2所示，所以若原系统 稳定，那么利用这种变换得到的脉冲传递函数 也一定是稳定的。 图5.2.2 双线性变换的映射关系0js 1z01z01z)(sg)(zg5.2替换法 n例5.2.2 对例5-1传递函数 n利用双线性变

18、换求其仿真数学模型，计算步长仍取 。n解：将替换关系(5.2.15)式代入 式中，则得12 . 01)()()(2sssxsysgst1)(sg5.2替换法n反变换后可得ttztzttzt4 . 04824 . 04122222211122 . 01121)()()(2zztzztzgzxzy852. 011. 1) 1(185. 06 . 464 . 512222zzzzzz)1()(2) 1(185. 0) 1(852. 0)(11. 1) 1(kxkxkxkykyky5.2替换法n对任意传递函数采用这种变换，它的脉冲传递函数即为 (5.2.18)n由(5.2.18)式可见， 与 是同阶的

19、。这种方法在概念上是很清楚的。实际运算时可能很麻烦，为此对常用的典型环节，与先求出 与 之间的各系数关系式，列成表格，在使用时查找是很方便的。此外目前许多程序都具有这种变换功能，在实际应用中可供选用。另外， 将这部分工作编成计算机程序，交付计算机来完成，也不是一件难事。112)()(zztssgzg)(zg)(sg)(zg)(sg5.2替换法5.2.3 状态方程的双线性变换状态方程的双线性变换n双线性变换公式不仅可以方便地用于传递函数，同时也可用于系统的状态方程。若系统的状态方程为 (5.2.19)n式中 分别为 维列向量；a为 维矩阵；b为 维矩阵；c为 维矩阵；为 维矩阵。n对(5.2.1

20、9)式左边进行拉氏变换，得)()()()()()(tdutctytbutaxtx yux,rmn,nnmnnr mr )()()(1sbuasisx5.2替换法n该式又可写成(5.2.20)n式中 t为计算步长。n现将 用(5.2.15)式代替，则得n进一步整理可得 )(222)(1subtatitssxits2)(2211)(1zubtatzzzx5.2替换法(5.2.21) )(2) 1() 1(2) 1()(1zubtzzatizzx)(2) 1(2)2(1zubtzatiatiz)() 1(1121zugzifzf)()(111121zugfiziffzi5.2替换法 式中 (5.2.

21、22) (5.2.23) n从(5.2.22)式及(5.2.23)式中消去，可得即atif21atif2221btg iff221211112fffi5.2替换法n将该式代换(5.2.21)式中的 的第二项i，则得 (5.2.24)n将(5.2.16)式中第二式做z变换并将(5.2.24)式代入，则得 (5.2.25) )(izi )()2()()(111211111211zugffffziffzizx)(2)()()(12112111112111211zugfffzigfffziffzi)(2)(1211211111zugfffzigf)()(2)()(1111211211zudgcfgff

22、fziczy5.2替换法n由(5.2.25)式可见，该式相当于一个离散状态方程的输出方程的z变换，若令 (5.2.26) n则以(5.2.26)式系数组成的离散状态方程如下： (5.2.27)211fffch 1212gfg111gcfde)()()()()() 1(keukhxkykgukfxkxdddd5.2替换法n 进一步研究证明，方程(5.2.27)式的脉冲传递函数就是(5.2.25)式。而(5.2.27)式即为(5.2.19)式的等价差分方程组，并且各 矩阵均为系数状态方程(5.2.19)式的系数矩阵 和计算步长 的函数，它可在计算机上迭代求解。n下面我们对双线性变换的一些性质作简单

23、的讨论。n前面已经说明了，若原连续系统是稳定的，采用双线性变换所得离散数学模型也是稳定的。现在再进一步研究其稳定特性。 对单位阶跃输入的稳态响应可以由(5.2.24)式求得ehgf,dcba,t)(zx5.2替换法n利用(5.2.22)式及(5.2.23)式，可得 112112111112)(zsgfffzigfx1121111)( 2gffff111121)(2gfffi1111111122)(2gfiatatgfatixdbabtat11225.2替换法n原连续系统的稳态值可以从(5.2.20)式令 求得，即 n显然，连续系统的稳态值与离散数学模型的稳态值是相同的。这说明采用双向性变换并不

24、会改变系统的稳态值。n采用双向性变换不仅可以保证系统的稳定性及稳态值不变，同时也具有一定的精度。这精度主要是指离散数学模型的频率特性或输出序列与连续系统的频率特性或输出序列是相近的。为了说明这一点，利用例5.2.2双线性变换后的系统来进行计算。对该题易得到0sbaxd1852. 01 . 1) 1(185. 0)(22zzzzg5.2替换法n将 代入后即可求得脉冲传递函数的频率特性为 n 将 代入上式得n若令 ，则可以得到 的频率特性。连续系统的频率特性可将 代入例5.2.2的 中即可求得。比较这两种频率特性可见，在低频段相差较小，当 超过0.8rad/s时，两个频率特性相差较大，但频率在高频

25、段，又重复对称出现。 tjez852. 01 . 1)() 1(185. 0)(22tjtjtjeeejgttetjsincos)sin11. 12(sin)852. 0cos11. 12(cos)sin22(cos) 1sin22(cos185. 0)(ttjttttjttjgst1)(zgjs )(sg5.2替换法这表明在一定的频率范围内双线性变换是有一定精度的，但在某频段内，双线性变换的频率特性会发生一定的畸变。产生这种畸变的主要原因是，通过(5.2.15)式变换后，s平面上无限长的 虚轴被压缩成平面上的单位圆，等效的说被压缩在 之内（ ， 为计算步长），即产生频率翘曲。n 应当指出，根

26、据(5.2.15)式进行变换，虽也将s平面左半部映射到z平面的单位圆内内，但它与 的映射关系并不相同，即按照(5.2.15)式将s平面上的一点 映射到z平面上的位置与按 映射位置并不相同。 j2sts2tistjeztjez5.3 零极点匹配法n零极点匹配法是获得连续传递函数等价离散数学模型的一种简单有效的方法。时间域函数的采样运算就是直接将s域的函数通过 (5.3.1)n关系映射到z平面。零极点匹配法（又称零极点映射法）不仅通过(5.3.1)将极点映射到z平面，而且还将 的零点映射到z平面。通常，传递函数的极点多于零点数，此时可以认为该传递函数在无穷远处有一定数目的零点，这些零点通过(5.3

27、.1)式映射关系被映射在z=1处。 tjez)(sg5.3 零极点匹配法n若给定的连续系统传递函数为 (5.3.2)n式中: 分别为重复零极点的实部与虚部。若 ，即极点数多于零点数 个。按上述意思，通过(5.3.1)式映射可得下述脉冲传递函数： (5.3.3) rjqjjjiminiiiidcssbasasksg11221122)()()()()(rjqjtcjtctnititmitpjjjjjjeztdezezeztbezezzkzg11221221)cos( 2)()cos2()() 1()(pmnrq)2()2(iiiidcba,p5.3 零极点匹配法n式中 是离散传递函数的增益， 可以

28、根据z域频率响应在某些关键频率上能较好地等价连续频率响应的原则来选择。通常对低通滤波器性质的传递函数或一般控制系统来说，它可以按直流( )增益相等的原则来选取。n 归纳起来，零极点匹配法的步骤如下： (1) 将 的所有极点按 映射到z平面，如 的实极点数在 处，则 有极点；倘若 有复数极点 ，则 有极点 。(2) 用同样的方法将 全部有限零点按 映射到z 平面。(3) 将 所有在 处的零点映射为 。 k0)(sgstez jjs)(zgtjjezjijjdcstjdtcjjjeezk)(sg)(sg)(zg)(zgstez s1z)(zg5.3 零极点匹配法(4) 选择离散数学模型的增益 ，选

29、择的方法是在给定的关键频率处使连续系统的增益与离散系统的增益相等。若按直流增益相等选择，则有 n例例5.3.1 现将零极点法用于前述例题，求差分方程的表达式。 n 解：解：该传递函数没有限零点，他的两个极点为 k10)()(zszgsg12 . 01)(2sssg995. 01 . 02, 1js5.3 零极点匹配法n对计算步长t=1s，经过映射后，在z与内的两个极点为 n由于没有限零点，故需补充两个在z=-1处的零点，所以 n 可以通过直流增益相等的原则来求得。因为 n所以 9 758. 07 492. 0995. 01 . 02, 1jeeeezjjbtt)9 j0.758-7 492.

30、0( )9 j0.7587 492. 0() 1( )(2zzzkzgk01)()(szsgzg1819. 0985. 014zk1819. 0985. 014zk5.3 零极点匹配法n故可得 =0.209，最后脉冲传递函数 为 (5.3.4)n该式可以转化为下述差分方程： n 应当说明，将无穷远处的零点映射在z=-1处，这样的处理，理由并不充分，主要的根据是在双线性变换 中就在分子上引入了零点（z+1）。这样做的好处，是有助于消除标准z变换中产生的频率混叠现象。819. 0985. 0) 1(209. 0)(22zzzzgk)(zg)1()(2) 1(209. 0) 1(819. 0)(98

31、5. 0) 1(kxkxkxkykyky1121zzts5.3 零极点匹配法n当然，如若认为无穷远处的零点是负实数，则由于 ，故可以将其看作是映射在z=0处的零点。此时可以将(5.3.3)式中 用 来代替，但经验表明，这种处理方法精度较差。n此外，为了使离散模型与原连续系统更为接近，常常还可以将无穷远处的零点映射到（0，-1）之间某处，此时可以将(5.3.3)式中须用 来代替 ，式中 。这样在 中将有两个参数 及 ，进行实际仿真计算，则是很简单的。0 tezpzpz)(pz) 1( 10)(zgkpz) 1( 5.3 零极点匹配法n对于某些传递函数 ，其直流增益可能等于零或无穷大，无法利用直流

32、增益相等的原则来确定传递函数的 。此时可以通过加入适当类型的输入信号，利用使其稳态值相等的方法选择 值。例如n根据上述步骤可得到 n由于 在s=0处有一零点，稳态增益等于零，加阶跃输入信号时，终值将等于零。为使终值)(sgk2) 1()()()(sssxsysg2)() 1)(1( )(tezzzkzgk)(sg5.3 零极点匹配法n而由 可得 n若是稳态值相等，则应选择 为 11)1(lim)(220ssssys2221)1 (2 ) 1()() 1)(1( 1lim)(ttzetkztzezzzkzzytekt2)1(2)(zgk5.3 零极点匹配法n最后可得离散数学模型为n当然，对这种系

33、统还可以不依据稳态增益来选取，而依据某个特殊频率时的增益相等的原则来确定。22)() 1)(1(2)1 ()(ttezzztezg5.4 计算机控制系统仿真n计算机控制系统是由离散部分（数字计算机或数字控制器）和连续部分（保持器或数模转换器以及被控对象） 两部分组合而成。如图5.4.1所示。 图5.4.1 计算机控制系统组成图5.4 计算机控制系统仿真n这类系统的模型与第四章所述连续系统离散化模型形式相同，都是差分方程，所以从仿真方法来讲两者都是相同的。但是要注意的是连续系统离散化模型是一种近似模型，因为在建立模型过程中人为的加入了采样开关和保持器，所以它只能是近似描述原来的连续系统。而计算机

34、控制系统仿真的特点是，数字控制器是实际存在的离散量运算，采样开关也是实际存在的不是人为虚拟的。采样周期ts是确定的，它与连续部分离散化虚拟采样周期t是不同的。下面来说明如何来考虑这两者之间的关系。一、采样周期及计算步距一、采样周期及计算步距5.4 计算机控制系统仿真n将图5.4.1所示的计算机控制系统用传递函数形式表示如图5.4.2所示。 图5.4.2 计算机控制系统图图中 被控对象的传递函数； 保持器传递函数； 数字控制器的脉冲传递函数。)(sg)(sgh)(zd5.4 计算机控制系统仿真n 是该计算机控制系统实际存在的采样开关的采样周期，这是计算机控制系统的一个重要参数。通常它是根据被控对

35、象的反应快慢事先选定的。比如，当被控对象的阶跃响应曲线如图5.4.3所示时，若选择数字控制器为p，i，d形式，那么可选 为0.3l左右。n对于图像5.4.2所示的计算机控制系统仿真，一般有两种情况：(1) 如果 可以求出，而且只要求计算输出量 ，不要求计算中间变量，那么，可以用同样的采样间隔 ，分别写出数字控制器的差分方程，然后进行数字仿真。此时，采样周期 ，也就是计算步距 。 stst)()(sgsgzh)(tytstst5.4 计算机控制系统仿真(2) 如果 不太好求出，或者不仅要求计算输出量 ，还要求计算出中间的其他状态变量，或者被控对象中有非线性环节，因此在做数字仿真时必须将连续部分分

36、成几部分来分别计算。此时，就要在每个连续部分前分别加上虚拟的采样器和保持器。 图5.4.3 被控对象的阶跃响应曲线 )()(sgsgzh)(ty5.4 计算机控制系统仿真n为了使计算有一定的精度，这些虚拟采样器的采样间隔t（也就是计算步距）就不能选得过分大，一般来讲 。此时，在做仿真计算中，就应对连续部分与离散部分分别计算，而不能用同一步距。为了仿真计算的方便，通常选取t使 ，其中n为正整数。 此时对整个计算机系统做仿真计算时，仿真程序将分为两个循环，一个是对连续部分，一个是对离散部分。由于 ，要对连续部分做n次计算后，才对离散部分做一次计算，如图5.4.4所示。 图中，连续部分的计算步距了；

37、 离散部分的采样周期； 要求的仿真时间。ttsnttsnttststft5.4 计算机控制系统仿真图5.4.4 计算机控制系统数字仿真程序的流程图 5.4 计算机控制系统仿真二、计算机控制系统的仿真方法二、计算机控制系统的仿真方法n现以图5.4.2所示的计算机控制系统为例进行说明。要实现计算机控制系统仿真，一般是计算离散部分的输出 ，其输入为 ，用 作为连续部分的输入，计算其输入 ，这一部分的仿真方法和连续系统完全相同，只要在确定连接矩阵w的系数时，把 去掉，看成是一个阶跃信号 加在系统上即可。(1) 数字控制器模型 )(2ne)()(1nyune)(2ne)(ny)(zd)(2ne5.4 计

38、算机控制系统仿真n 数字控制可以有各种控制规律，既可以有不同的脉冲传递函数或差分方程，但建立模型的方法基本相同，这里先以常用的pid（比例、积分和微分） 控制规律为例来说明建立模型的方法。 (5.4.1)式中： 比例系数； 积分时间常数； 为分时间常数。n 现若用z变换方法把 变换 为，则可得ststksesesddip11)()()(12pkitdt)(sd)(zd5.4 计算机控制系统仿真 (5.4.2)n式中 n由(5.4.2)式进行反变换，可得差分方程如下： (5.4.3) (5.4.3)式及所求的数字控制器的一种模型形式。它与连续部分离散化所得的差分方程形式相同。因此计算机控制系统的

39、仿真就是要研究高阶差分方程的仿真方法。1211121)1 ()1 ()()()(zzkzkkzezezddpiipittkk ttkkdpd) 2() 1()2()()() 1()(21122neknekknekkkneneddpdpi5.4 计算机控制系统仿真(2) 差分方程仿真n设系统闭环脉冲传递函数的一般形式为 (5.4.4)n其差分方程为 n 解上式即可求出采样时刻的输出。nnnnmmmmzazazazabzbzbbzuzyzf)1(12211)1(11101)()()()() 1()2() 1()(121nkyankyakyakyakynn)() 1() 1()(110mkubmku

40、bkubkubmm5.4 计算机控制系统仿真n若需要求出控制器的输入和输出信号时，则可分别求出控制器部分 及连续部分 的差分方程。对于每一个采样点，这两部分各计算一次，然后对控制器的信号进行综合，以便得到下一次计算控制器的输入。如此循环迭代下去，便可逐步求出各采样时刻的响应。以上两种算法所得到的结果是相同的，但后一种算法的运算时间要长一些。n 以上两种算法，最后都归结为高阶差分方程。下面介绍用计算机解这种方程的方法及其程序实现。n假定高阶差分方程具有如下形式： )(zd)(zg5.4 计算机控制系统仿真 (5.4.5)n式中：输出变量；输入变量。n由(5.4.5)式可知，在解此差分方程时，要用

41、到计算时刻以前若干个采样时刻的输入值和输出值，这可在内存中设置若干个单元，见数据存储起来，以便在计算时使用。n下面首先讨论输出变量的处理方法。由(5.4.5)式可知，输出变量x的阶数n为，因此，需要在内存中设置n个单元以存放计算时刻以前n个采样时刻的输出量。这些单元的安排如图5.4.5(a)所示。)()2() 1()(21nkxakxakxakxn)() 1()(10mkebkebkebmxe5.4 计算机控制系统仿真 图5.4.5 变量存储单元的安排5.4 计算机控制系统仿真n在计算中，现在时刻前n个输出量 分别从 ，单元中取出。在每一次计算中其操作顺序都是“取出平移存入”。这使得“存入”，

42、是指计算所得的现时刻的输出量由存入到第n单元中去。n对于输入变量也可采用上述相似的方法进行处理。根据(5.4.5)式，e的阶数为m，在内存中设置m+1个单元存放e的数据。其安排如图5.4.5(b)所示。 ) 1(),2(,),1(),(kxkxnkxnkx) 1(),2(,),1(),(kxkxnkxnkxnn, 1, 2 , 15.4 计算机控制系统仿真n在计算中，现时刻前m个输入变量 分别从 单元中取出。在每一次计算中，其操作顺序都是“存入取出平移”。这里的“存入”，是指计算所得到的现时刻的输入 存入到第m+1单元中去。n差分方程仿真程序：n#includen#includen#inclu

43、denvoid main() ),1(,),1(),(kemkemke),1(,),1(),(kemkemkem, 2 , 1)(ke5.4 计算机控制系统仿真n n int i,j,n;n int n,m,k;n double *x,*e,*a,*b,e,x;n clrscr( ) ;n printf(“please input n,m,k,e:”);n scanf(“%d %d %d %lf”,&n,&m,*k,&e);n a= malloc(n*sizeof(double);b= malloc(m*sizeof(double); n x= malloc(n*siz

44、eof(double);e= malloc(m*sizeof(double);n for(i=0;in;i+) n n printf(“input a%d;”,i);n scanf(“%lf”%ai);n xi=0;n 5.4 计算机控制系统仿真n for(i=0;ik;i+)n n printf(“input b%d;”,i);n scanf(“%lf”%bi);n ei=0;n n em-1=e;nfor(i=0;ik;i+)n n x=0n for(j=0;jk;j+)n x=x+aj*em-j-1;5.4 计算机控制系统仿真nfor(j=0;jm;j+)n x=x+bj*em-j-1;

45、n for(j=0;jk-1;j+)n xj=xj+1;nfor(j=0;jm-1;j+)n ej=ej+1;nxk-1=x;n printf(“%d %lfk”,i,x);n nbioskey(1);n 5.4 计算机控制系统仿真三、纯时延环节的数字仿真三、纯时延环节的数字仿真n我们现在研究纯时延环节的数字仿真方法，我们将看到纯时延环节其实是一种特殊的差分方程描述的系统，因此可以利用上面讨论的方法解决纯时延环节的仿真。假设纯时延环节的输入量为 ，输出量为 ，延迟时间为 ，则时延环节的传递函数 (5.4.6)n再假设采样周期为 ，且 (5.4.7)(tu)(tysesusysg)()()(20

46、cctt5.4 计算机控制系统仿真n式中， 为正数部分， 为小数部分。则有 (5.4.8)n对(5.4.8)式作z变换 (5.4.9)n再作z反变换，得到差分方程 (5.4.10)0c2ctsccesusysg)(20)()()()(20)()()(cczzuzyzg)()(20ccnuny5.4 计算机控制系统仿真n若 为 的整数倍，则 n其仿真方法同差分方程一样，如图5.4.6所示 图5.4.6 纯时延环节变量存储单元的安排n利用basic语言的仿真程序为：nq(c0+1)=uny=q(1)nfor i=1 to c0n q(i)=q(i+1)nnext it)()(0cnuny5.4 计

47、算机控制系统仿真n当 时，则因 ，因此我们可以将 的值取在 和 两个时刻之间，故可以利用线性插补公式求得 (5.4.11) n增设一个内存单元存放 ，故可取 个内存单元，如如图5.4.7所示 图5.4.7 纯时延环节变量存储单元的安排02c10200cccc)(ny)(0cnu)1()()1 ()(0202cnuccnucny)1(0 cnu20c)1(0 cnu5.4 计算机控制系统仿真n利用basic语言的仿真程序为：nq(c0+2)=uny=(1-c2)*q(2)+c2*q(1)nfor i=1 to c0+1n q(i)=q(i+1)nnext i四、采样周期改变引起仿真模型的变化四、

48、采样周期改变引起仿真模型的变化n 为了减少计算量，加快仿真的速度，对计算机控制系统进行数字仿真时，通常会遇到这样的情况，即原来的 比较小，现在希望用一个比较大的 来进行仿真。这样要对离散部分的仿真模型进行修改，但当给出的离散部分仿真模型不是含 时，须将原z域零极点通过映射关系对应到s域中，再采用 映射到z域中找到 改变为 后新的零极点。 tssttssttsst5.4 计算机控制系统仿真n下面通过一个实际的例子来说明一下当采样周期改变时，如何改变仿真模型。n 假设在一个数字驾驶仪中有一个数字校正环节，它的采样周期 ，它的脉冲传递函数为 (5.4.12)n若用 来进行数字仿真，那么它的仿真模型可

49、以用差分方程来表示，即 (5.4.13)sts04. 064. 098. 062. 2)()()(zzzuzyzd)1(98. 0)(62. 2) 1(64. 0)(kukukykysts04. 05.4 计算机控制系统仿真n今要求仿真时 ，那么仿真模型应当怎样变化呢？n假设当 时， 可以表示成 ，它可以看做是 的仿真模型，为了使 时的仿真模型与 时的原模型具有相同的特性，首先要求 与 两个脉冲传递函数在s平面上映射具有相同的零极点，同时稳态值相同。n已知 在z平面上的一个极点 ，同时有一个极点 ，它们映射到s平面上，有一个极点 ，一个零点 。sst1 . 0 )(zd)( zdsts04.

50、0)( zd)(zd64. 0pz98. 0zzpszs16.1164. 0ln25ln1pspzts505. 098. 0ln251zszztssst1 . 0 sst1 . 0 sst1 . 0 )(zd5.4 计算机控制系统仿真n根据零极点相配原则，当 时，有 n 故 (5.4.14)n 现在还有一个增益系数 要加以确定。它可以根据原型与模型具有相同的稳态值这一原则来确定。n假设在 的输入端加一个单位阶跃函数，则根据终值定理，输出的终值为 sst1 . 0 9508. 0505. 01esezzstz327. 016.111 . 0esezpstp7 327. 08 950. 0)( z

51、zkzdzzk)(zd36. 002. 062. 2164. 098. 062. 21lim)(1zzzzzzyz5.4 计算机控制系统仿真n而仿真模型的终值为 n若要求有 ，则可得n因此可得改变采样周期后的仿真模型为 (5.4.15)3 672. 02 049. 0)( ky9 988. 136. 02 049. 002. 062. 23 672. 0zk)( )(yy7 327. 08 0.950-z9 988. 1)( zzd5.4 计算机控制系统仿真n相依的差分方程为 (5.4.16)n在matlab语言中有指令d2d()，该指令可以完成由于采样周期的改变而引起的模型变化。该指令的调用格式为： sys1 = d2d(sys,ts) sys1 = d2d(sys,nd)n 调用说明：该指令将原离散时间系统的采样时间改变为采样时间为ts秒的新模型sys1。其实该指令相当于指令： sys1 = c2d(d2c(sys),ts)n指令sys1 = d2d(sys,nd)是对系统sys在输入端加一个时延，该时延必须是采样时间的nd倍。如果系统sys是单输入，则nd