连续系统的数字仿真.doc

1、第 4 章 连续系统的数字仿真数值积分法上一章从仿真原理方面讨论了连续系统的仿真方法。本章将从构造积分器的角度再对仿真方法做进一步的讨论。4.1 欧拉法数值积分法是把微分方程化成积分运算，再进一步化成代数运算的过程，主要解决如何构造一个积分器，然后求出积分器的差分方程的问题。有了积分器就能很容易地对系统进行仿真了。数值积分法最初是从数值计算 56的角度得到的。但是为了和第 3 章所讨论的方法统一起来，我们用插入离散再现环节的方法，以状态空间描述为基础，推出线性系统数值积分法的仿真模型 9。仍假设线性系统的状态空间描述为（4-1）BUAX（4-2）DCY式中 维状态向量；X1n 维输入向量；Ur

2、 维系统矩阵；A 维输入矩阵；B 维输出向量；Ym 维输出矩阵；Cn 维传递矩阵。Dr此系统的方框图如图 4.1 所示。在系统积分器入口 e 处加入离散再现环节，则可得到连续系统的另一离散相似系统，如图 4.1 所示BDAH(s) CC(s)td0X(0)ehTe(t)kT离 散 再 现 过 程图 4.1 线性定常系统的另一种离散相似系统框图从图 4.1 中，可得到（4-3）thdeXt0)()(当 时， （4-4 ）KTtkTtk)(当 时， （4-5）)1(Tkhte)1(01用式（4-5）减式（4-4）可得到(4-6)TkhdteXTk)1()1(当取 时，即)(,scesHT（4-7）

3、khTkt)1(把式（4-7）代入式（4-6）得)()()1(eXT简记为 （4-8）kk式（4-8）被称为欧拉公式。欧拉公式可以从图 4.2 的几何图形中得到解释。有了式（4-8） ，就很容易求出系统式（ 4-1） 、式（4-2）的差分方程（4-9）)()()1kTBUXAIkX（4-10）1(DkCY从图 4.2 可以看出，这种方法精度低（截断误差与 成正比 56 ） ，但是它的公式2T非常简单，不用求，因此在实时仿真)()(nm、中它的应用非常广泛。图 4.2 欧拉公式的几何解释4 2 梯形法为了提高仿真精度，离散再现环节采取图 4.3 的形式。即)1()21Tkeeh（4-11）tkT

4、kTX()e(kT)(k+1)T+1 误 差te零 阶保 持 1e(kT)e(t) )(teh2/TSe/1 12图 4.3 梯形公式的离散-再现环节框图把式（4-11）代入式（4-6）可得（4-12）)1()2)(1( keTkX式（4-12）称为梯形公式，其几何解释如图 4.4 所示。误 差 )1()21TkekTX()e(kT)(k+1)Tt+1)(te图 4.4 梯形公式的几何解释由图 4.1 可知（4-13）)()()(kBUAke（4-14）)11X把式（4-13） 、 （4-14）代入式（4-12） ，得到系统的差分方程（4-15）)1()2)(2)() kUBTkATIkX（4

5、-16）11( DUkCY式（4-15）是一个隐式，因为求 时，等式右边还有未知数 。为了)(kX)(kX得到该式的显式形式，可把含 的项移到方程左边，再整理而得到，即系统解的显式公式为（4-17）)(1()2()2()()1 11 kUBATIkATIIkX 显然，式（4-17）要比式（ 4-9）的精度高些。但是，当系统阶次较高时，求是比较困难的。为此，在计算式（4-14）时，需要先估算 的值1)2(ATI )1(kX。此时，可以设定离散再现过程只有第 1 路信号（见图 4.3） ，根据式（4-8）0kX则有)()1(0kTeXk把上式代入式（4-14） ，用 代替 ，则有0)1（4-18）

6、()()(kBUekAke把式（4-13） 、式（4-18）代入式（4-12）可得到系统解的简化显式公式为（4-)1(2)(2()2()1 kBUTATXTIX19）此式没有矩阵求逆的运算，所以比式（4-17）容易计算。在仿真时，也可直接使用式（4-13） 、式（4-18） 、式（4-12）进行计算。【例 4-1】 已知一多变量系统的结构框图如图 4.5 所示，请用梯形公式对此系统进行仿真，并输出 的仿真结果。21y、12.7-1.26 1/s1/s0.111U2U-0.42350.8534+-0.1270.1260.423-0.853tt1图 4.5 多变量系统结构框图【解】 根据图 4.5

7、，可得到该系统的状态方程及输出方程为 212121 8534.06.70853.62.0471 uxx2121xy依据式（4-13）有 )(8534.026.17)(0853.126.047)( 212121 kukxke根据式（4-18）有 )2(1.0357)()(.)( 221121 kkTexke根据式（4-12）有 )1()2)()1( 21212 keTkxkx根据式（4-16）有 )1(0.)(221kxky计算步距及仿真时间的估算 5.0)816.0(7.)05()( nTD2853.11S仿真程序清单如下：#include “graphics.h“ /*头文件声明*/#inc

8、lude “math.h“#include “string.h“int VN=2; /*定义输出变量的个数*/float Outputy2500,ST,DT; /*定义输出变量的最大维数*/int LP,i,j; /*定义仿真步数 */float x10,x11,x20,x21,e10,e20,e11,e21,y1,y2; /*定义参数类型 */float u1,u2;void Dispcurve(); /*画图子程序*/main() x10=0;x20=0; /*参数赋初值*/Outputy00=0;Outputy10=0;u1=1;u2=1;y1=0;y2=0; ST=100;DT=0.5

9、;LP=(int)(ST/DT);for(i=1;i=LP;i+) /*计算各个差分方程*/ e10=-0.12771*x10+0.04233*x20+12.771*u1-0.04235*u2;e20=0.01262*x10-0.08533*x20-1.262*u1+0.08534*u2;e11=-0.12771*(x10+DT*e10)+0.04233*(x20+DT*e20);e11=e11+12.771*u1-0.04235*u2;e21=0.01262*(x10+DT*e10)-0.08533*(x20+DT*e20);e21=e21-1.262*u1+0.08534*u2;x11=x

10、10+DT*e10/2.0+DT*e11/2.0;x21=x20+DT*e20/2.0+DT*e21/2.0;y1=0.01*x11;y2=x21;Outputy0i=y1;Outputy1i=y2; /*输出值*/x10=x11;x20=x21; /*把 X(k+1)放入 X(k)*/Dispcurve(); /*调用画图子程序 */多变量系统的仿真结果见图 4.6。图 4.6 多变量系统的仿真结果4 3 龙格库塔（Runge-Kutta）法在非实时仿真中，有时需要更高的精度。下面的两节中将再介绍两种更精确的方法及其离散再现环节。对于图 4.1 所示的系统，取如图 4.7 所示的离散再现环节

11、，即（4-20）)2(61)43eeth Tkt)1(把式（4-20）代入式（4-6）可得零 阶保 持 1e(kT) )(teh2/31 12TSe/1216/TSet) 2334图 4.7 四阶龙格库塔公式的离散-再现环节框图(4-21)2(6)1( 431eTkX由图 4.1 及图 4.7 可有（4-22）)()(1BUAe（4-23）212kk（4-24）)()(3Xe（4-25）14kBUkA由于计算式（4-21）中的 时，需要已知 、 、 中的 、 。)(2e34)21(kX)(k因此，在计算式（4-21）之前应先估算 、 、 中 和 的近似值234)(和 。在估算 中的 时，设定离

12、散再现过程只有第 1 路)21(0kX)(0ke)1(kX信号（见图 4.7） ，则102)(Txx把上式代入式（4-23）得（4-26）)()(12 kBUekAe估算 中的 时，设定离散再现过程只有第 2 路信号（见图 4.7） ，则3X20)(2(ekxx把上式代入式（4-24）得（4-27）)1()(23 kBUTAe估算 中的 时，设定离散再现过程只有第 3 路信号（见图 4.7） ，则41kX30)(1(Tekx把上式代入式（4-25）得（4-28）)1()(34BUAe式（4-21） 、式（4-22） 、式（4-26）式（4-28 ）称为四阶龙格库塔公式。此公式的总体截断误差与

13、56成正比。5T如果系统的输入 U 为阶跃函数或斜坡函数，则（4-29）2/)1()21(kk把式（4-22） 、式（4-26）式（4-29）代入式（4-21 ）可得)1()2462( )(2483()46)13 33232 kBUATI kBUATATIxATkx（4-30）式（4-30）是四阶龙格库塔公式在输入为阶跃或斜坡函数时的总体表示，在仿真前，先求出方程的系数，仿真时就只有简单的代数运算了，这样仿真的速度要比用分离公式的速度快。非常有趣的是式（4-30）与式（ 3-51）完全相同。由此说明，当系统的输入为阶跃或斜坡函数时，四阶龙格库塔公式即是把离散-再现环节加在系统的入口处，使用三角

14、保持器，取 的定义式到 的 4 次项所得到的差分方程。由此可见，使用三角保持器比使用四Atet阶龙格库塔公式要精确。同理也可证明，取 的定义式到 的 2 次项，所得到的差分方Atet程即是梯形公式，取 的定义式到 的 1 次项，所得到的差分方程即是欧拉公式 9。Atet4.4 阿达姆斯（Adams ）法阿达姆斯方法是线性多步法，在计算 时，需要知道该时刻以前几步的结果。)(kX相比之下，欧拉法、梯形法、龙格库塔法都是单步法。对于图 4.1 所示的线性定常系统，取图 4.8 所示的离散再现环节，即 )3(9)2(37)1(59)(241)( TkekeTkeTteh （4-31）k零 阶保 持

15、1e(kT) )(teh2459Te243729Te3t)图 4.8 阿达姆斯法的离散再现环节框图把式（4-31）代入式（4-6）可得（4-32 ）)3(9)2(37)1(59)(24)(1( kekeTkX此式称为阿达姆斯四阶精度显式。由图 4.1 可知（4-33）)()()(kBUAke（4-34 ）11X（4-35）)2()()2(kkke（4-36）33BUA四阶阿达姆斯法和四阶龙格库塔法均有四阶精度。但是用阿达姆斯公式时，第一步计算需要知道 、 、 、 各个时刻的值。在初始条件不为零的情况下，)0(x1)2(x、 、 应该用其他方法求出，所以使用起来不太方便。但是，对于热工)1(x2

16、3系统的仿真，一般是零初始条件，所以这种方法还是实用的。【例 4-2】 已知一主汽压力系统，其框图如图 4.9 所示，试用四阶龙格库塔公式，四阶阿达姆斯公式对此系统进行仿真，输出 的仿真结果（对系统做定压扰动，扰动量为tP1，系统的初始条件为零） 。 )1(sTi)6.531)(49.0s1 601 tp0p+-图 4.9 主汽压力调节系统框图【解】 把图 4.9 转换成图 4.10 所示的形式，并按图 4.10 所示设状态变量为 、1x、 。2x3 1x2x3x1sTi s54639.0s6.51tp0p+图 4-10 主汽压力调节系统的变化框图则有 0321321 5469.6.5.31049.05469.0 pTxTx ii 3pt根据式（4-21） 、式（4-22） 、式（4-26）式（4-28 ） ，即可得到使用四阶龙格库塔公式时的系统差分方程 )(26)1( 32143312132 kxeeTkx1 12 203 30().690.69.69()545454103.i iTTe xkpk 1112 2203 331()0.6910.69.691()5454542()3.i iTxkee pkxke