连续系统的数字仿真.doc

第4章 连续系统的数字仿真数值积分法上一章从仿真原理方面讨论了连续系统的仿真方法。本章将从构造积分器的角度再对仿真方法做进一步的讨论。4.1 欧拉法数值积分法是把微分方程化成积分运算，再进一步化成代数运算的过程，主要解决如何构造一个积分器，然后求出积分器的差分方程的问题。有了积分器就能很容易地对系统进行仿真了。数值积分法最初是从数值计算56的角度得到的。但是为了和第3章所讨论的方法统一起来，我们用插入离散再现环节的方法，以状态空间描述为基础，推出线性系统数值积分法的仿真模型9。仍假设线性系统的状态空间描述为 （4-1） （4-2）式中 维状态向量；维输入向量；维系统矩阵；维输入矩阵；维输出向量；维输出矩阵；维传递矩阵。此系统的方框图如图4.1所示。在系统积分器入口e处加入离散再现环节，则可得到连续系统的另一离散相似系统，如图4.1所示图4.1 线性定常系统的另一种离散相似系统框图从图4.1中，可得到 （4-3）当 时， （4-4）当 时， （4-5）用式（4-5）减式（4-4）可得到(4-6)当取 时，即 （4-7）把式（4-7）代入式（4-6）得 简记为 （4-8）式（4-8）被称为欧拉公式。欧拉公式可以从图4.2的几何图形中得到解释。有了式（4-8），就很容易求出系统式（4-1）、式（4-2）的差分方程 （4-9） （4-10）从图4.2可以看出，这种方法精度低（截断误差与成正比56），但是它的公式非常简单，不用求，因此在实时仿真中它的应用非常广泛。图4.2 欧拉公式的几何解释4 2 梯形法为了提高仿真精度，离散再现环节采取图4.3的形式。即 （4-11）图4.3 梯形公式的离散-再现环节框图把式（4-11）代入式（4-6）可得 （4-12）式（4-12）称为梯形公式，其几何解释如图4.4所示。 图4.4 梯形公式的几何解释由图4.1可知 （4-13） （4-14）把式（4-13）、（4-14）代入式（4-12），得到系统的差分方程 （4-15） （4-16） 式（4-15）是一个隐式，因为求 时，等式右边还有未知数。为了得到该式的显式形式，可把含的项移到方程左边，再整理而得到，即系统解的显式公式为 （4-17） 显然，式（4-17）要比式（4-9）的精度高些。但是，当系统阶次较高时，求是比较困难的。为此，在计算式（4-14）时，需要先估算的值。此时，可以设定离散再现过程只有第1路信号（见图4.3），根据式（4-8）则有 把上式代入式（4-14），用代替，则有 （4-18）把式（4-13）、式（4-18）代入式（4-12）可得到系统解的简化显式公式为 （4-19）此式没有矩阵求逆的运算，所以比式（4-17）容易计算。在仿真时，也可直接使用式（4-13）、式（4-18）、式（4-12）进行计算。 【例 4-1】 已知一多变量系统的结构框图如图4.5所示，请用梯形公式对此系统进行仿真，并输出的仿真结果。 图4.5 多变量系统结构框图【解】 根据图4.5，可得到该系统的状态方程及输出方程为依据式（4-13）有根据式（4-18）有根据式（4-12）有根据式（4-16）有计算步距及仿真时间的估算仿真程序清单如下：#include graphics.h /*头文件声明*/#include math.h#include string.hint 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; 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=x10+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所示的离散再现环节，即 （4-20）把式（4-20）代入式（4-6）可得图4.7 四阶龙格库塔公式的离散-再现环节框图 (4-21)由图4.1及图4.7 可有 （4-22） （4-23） （4-24） （4-25）由于计算式（4-21）中的时，需要已知、中的、。因此，在计算式（4-21）之前应先估算、中和的近似值和。在估算中的时，设定离散再现过程只有第1 路信号（见图4.7），则 把上式代入式（4-23）得 （4-26）估算中的时，设定离散再现过程只有第2 路信号（见图4.7），则 把上式代入式（4-24）得 （4-27）估算中的时，设定离散再现过程只有第3 路信号（见图4.7），则 把上式代入式（4-25）得 （4-28）式（4-21）、式（4-22）、式（4-26）式（4-28）称为四阶龙格库塔公式。此公式的总体截断误差与56成正比。如果系统的输入U为阶跃函数或斜坡函数，则 （4-29）把式（4-22）、式（4-26）式（4-29）代入式（4-21）可得（4-30）式（4-30）是四阶龙格库塔公式在输入为阶跃或斜坡函数时的总体表示，在仿真前，先求出方程的系数，仿真时就只有简单的代数运算了，这样仿真的速度要比用分离公式的速度快。非常有趣的是式（4-30）与式（3-51）完全相同。由此说明，当系统的输入为阶跃或斜坡函数时，四阶龙格库塔公式即是把离散-再现环节加在系统的入口处，使用三角保持器，取的定义式到的4次项所得到的差分方程。由此可见，使用三角保持器比使用四阶龙格库塔公式要精确。同理也可证明，取的定义式到的2次项，所得到的差分方程即是梯形公式，取的定义式到的1次项，所得到的差分方程即是欧拉公式9。4.4 阿达姆斯（Adams）法阿达姆斯方法是线性多步法，在计算时，需要知道该时刻以前几步的结果。相比之下，欧拉法、梯形法、龙格库塔法都是单步法。对于图4.1所示的线性定常系统，取图4.8所示的离散再现环节，即 （4-31） 图4.8 阿达姆斯法的离散再现环节框图把式（4-31）代入式（4-6）可得 （4-32）此式称为阿达姆斯四阶精度显式。由图4.1可知 （4-33） （4-34） （4-35） （4-36）四阶阿达姆斯法和四阶龙格库塔法均有四阶精度。但是用阿达姆斯公式时，第一步计算需要知道、各个时刻的值。在初始条件不为零的情况下，、应该用其他方法求出，所以使用起来不太方便。但是，对于热工系统的仿真，一般是零初始条件，所以这种方法还是实用的。【例 4-2】 已知一主汽压力系统，其框图如图4.9所示，试用四阶龙格库塔公式，四阶阿达姆斯公式对此系统进行仿真，输出的仿真结果（对系统做定压扰动，扰动量为1，系统的初始条件为零）。 图4.9 主汽压力调节系统框图【解】 把图4.9转换成图4.10所示的形式，并按图4.10所示设状态变量为、。 图4-10 主汽压力调节系统的变化框图则有根据式（4-21）、式（4-22）、式（4-26）式（4-28），即可得到使用四阶龙格库塔公式时的系统差分方程由于为单位阶跃函数，所以 （时）。根据式（4-32）、式（4-33）、式（4-34）、式（4-35）、式（4-36），即可得到使用四阶阿达姆斯方法时的系统差分方程 由于为单位阶跃函数，所以当时， ，而第一步运算时， ，。由于零初始条件，所以。由于求龙格库塔公式及阿达姆斯公式系数时所用到的公式形式相同，故在仿真程序中使用子程序计算。在四阶龙格库塔法及阿达姆斯法的运行程序块中、分别代表、；、分别代表、在四阶龙格库塔法中，数组代表，在四阶阿达姆斯法的运行程序块中数组代表。、用下式估算：仿真程序清单如下：#include graphics.h /*头文件声明*/#include math.h#include string.hconst VN=2; /*输出变量的个数*/float Outputy2501; /*输出变量的最大维数*/float e45; /*定义变量的类型*/float kk;int LP,i,j,L;float ST,DT;float TAT,Ti;float x10,x20,x30, x11,x21,x31;float p0;float z1,z2,z3;void Dispcurve(); /*画图子程序*/ void Firoutput (); /*第一个输出函数子程序块*/void Secoutput(); /*第二个输出函数子程序块*/main() TAT=1;Ti=60; p0=1; /*参数赋初值*/ ST=600;DT=2; LP=ST/DT; Firoutput(); /*调用第一个输出函数子程序*/ Secoutput(); /*调用第二个输出函数子程序*/ Dispcurve(); /*调用画图子程序*/void Firoutput() x10=0;x20=0;x30=0; Outputy00=x30; for(i=1;iLP;i+) z1=x10;z2=x20;z3=x30; /*变量赋初值*/ /*计算e11,e21,e31*/ e11= -1.0*z3/TAT/Ti+1.0*p0/TAT/Ti; e21=0.6369*z1/54.0-1.0*z2/54.0-0.6369*z3/TAT/54.0; e21=e21+0.6369*p0/54.0/TAT; e31=1.0*z2/53.6-1.0*z3/53.6; z1=x10+DT*e11/2.0; z2=x20+DT*e21/2.0; z3=x30+DT*e31/2.0; /*计算e12,e22,e32*/ e12= -1.0*z3/TAT/Ti+1.0*p0/TAT/Ti; e22=0.6369*z1/54.0-1.0*z2/54.0-0.6369*z3/TAT/54.0; e22=e22+0.6369*p0/54.0/TAT; e32=1.0*z2/53.6-1.0*z3/53.6; z1=x10+DT*e12/2.0; z2=x20+DT*e22/2.0; z3=x30+DT*e32/2.0; /*计算e13,e23,e33*/ e13= -1.0*z3/TAT/Ti+1.0*p0/TAT/Ti; e23=0.6369*z1/54.0-1.0*z2/54.0-0.6369*z3/TAT/54.0; e23=e23+0.6369*p0/54.0/TAT; e33=1.0*z2/53.6-1.0*z3/53.6; z1=x10+DT*e13; z2=x20+DT*e23; z3=x30+DT*e33; /*计算e14,e24,e34*/ e14= -1.0*z3/TAT/Ti+1.0*p0/TAT/Ti; e24=0.6369*z1/54.0-1.0*z2/54.0-0.6369*z3/TAT/54.0; e24=e24+0.6369*p0/54.0/TAT; e34=1.0*z2/53.6-1.0*z3/53.6; x11=x10+DT*(e11+2*e12+2*e13+e14)/6.0; x21=x20+DT*(e21+2*e22+2*e23+e24)/6.0; x31=x30+DT*(e31+2*e32+2*e33+e34)/6.0; Outputy0i=x31; x10=x11;x20=x21;x30=x31; void Secoutput() Outputy10=0; /*变量赋初值*/ for (i=1;i=3;i+) for(j=1;j=4;j+) eij=0; x10=0;x20=0;x30=0; for(L=1;L=LP;L+) /*计算e*/ e11= -1.0*x30/TAT/Ti+1.0*p0/TAT/Ti; e21=0.6369*x10/54.0-1.0*x20/54.0-0.6369*x30/TAT/54.0; e21=e21+0.6369*p0/54.0/TAT; e31=1.0*x20/53.6-1.0*x30/53.6; kk=55.0*e11-59.0*e12+37.0*e13-9.0*e14; x11=x10+DT*kk/24.0; kk=55.0*e21-59.0*e22+37.0*e23-9.0*e24; x21=x20+DT*kk/24.0; kk=55.0*e31-59.0*e32+37.0*e33-9.0*e34; x31=x30+DT*kk/24.0; Outputy1L=x31; x10=x11;x20=x21;x30=x31;for (i=1;i=3;i+) for(j=1;j4;j+) ei5-j=ei4-j; 主汽压力调节系统仿真结果如图4.11所示。从图中可以看出，两种方法的仿真曲线完全重合。图4.11主气压力系统的仿真结果四阶阿达姆斯法与四阶龙格库塔法一样，是通过数值计算的方法推导而来，这里选择如此复杂的离散再现环节的目的是为了说明，怎样用离散相似法得到四阶阿达姆斯公式。通过上面介绍的各种方法可以看出，用离散相似法很容易构成各种积分公式，也由此说明，微分方程的数值计算，实质上就是怎样构造积分器的问题。也可以看出，在构造积分器时，补偿环节综合的信号越多，精度也越高，但无论怎样，这些信号的权系数之和都等于1。因此可以说，数值积分法只是对再现信号的相位进行了补偿。根据上述结论，可以改变补偿环节的相位及幅值（权系数），设计出更多的数值积分公式。45 非线性系统数值积分公式假设非线性系统的状态空间描述为 （4-37）式（4-37）的离散相似系统框图如图4.12所示。图4.12 非线性系统的离散相似系统图4.12与图4.1的形式相同，因此按着线性系统同样的方法加入不同的离散再现环节，即可得到各种积分公式，这些公式列于表4.1中表4.1 的几种数值积分公式算法名称积分公式（差分方程）总体截断误差欧拉公式（可用于实时仿真）预报校正梯形公式三阶龙格库塔公式（可用于实时仿真）四阶龙格库塔公式实时四阶龙格库塔公式（可用于实时仿真）四阶阿达姆斯公式（可用于实时仿真）注：总体截断误差阶次越高，系统的解越精确；中的r称为各种算法的阶次56实时仿真所用的算法与非实时仿真所用的算法稍有不同，前者求时不能用到X和U的某些先前值。在实时仿真计算时，如果当前时刻还没计算就要求知道该时刻或更前时刻的值，这显然是不可能的。因此，并非所有的算法都能直接用于实时仿真。为了便于使用，在表4.1中注明了哪些算法可用于实时仿真。【例 4-3】 某随动系统简化整理后的框图如图4.13所示。系统的摩擦非线性特性如图4.14所示。系统输入为一斜坡函数。试用欧拉法对此系统进行仿真，并输出r及c 的仿真结果。已知条件为 当时 ；当时 ；当时 且 （）；图4.13 随动系统简化框图图4.14 系统的摩擦非线性特性【解】由图4.13可得到该系统的状态方程根据表4.1中的欧拉公式可得到系统的差分方程仿真程序清单如下：#include graphics.h /*头文件声明*/ #include math.h #include string.h const VN=2; /*输出变量的最大个数*/ float Outputy2501; /*输出变量的个数和维数*/ int LP,i,j; /*定义仿真步数类型*/ float ST,DT; /*仿真时间和仿真步距*/ float x10,x11,x20,x21,c;void Dispcurve(); /*画图子程序*/ main() /*主运行程序块*/ x10=x20=0; c=x10; /*参数初始化*/ ST=10;DT=0.1; LP=ST/DT; Outputy00=0; /*输出参数的初始值*/ Outputy10=c; for(i=1;i0) x21=x20+DT*(10*i*DT-x10-0.5*x20-1); if (x200) x21=x20+DT*(10*i*DT-x10-0.5*x20+1); if (x20=0)&(abs(10*i*DT-x10)2) x21=x20+DT*(10*i*DT-x10-0.5*x20-2); if (x20=0)&(10*i*DT-x10)-2) x21=x20+DT*(10*i*DT-x10-0.5*x20+2); c=x11; Outputy0i=10*i*DT; Outputy1i=c; x10=x11;x20=x21; /*把X(k+1)放入X(k)*/ Dispcurve(); /*调用画图子程序*/ 随动系统的仿真结果如图4.13所示。图4.15 随动系统的仿真结果通过本节的叙述可以看出，数值积分法和离散相似法在本质上没有区别，仅仅是采样开关的位置不同，保持器、补偿器的结构和参数不一样。这两种方法原本是从不同的角度推导出来的，它们却又有内在的联系，充分体现了仿真方法的统一性。表4.1所列举的是传统的数值积分公式，它不仅适合于线性系统，也可用于非线性系统。但本节开始推导出来的公式是建立在线性状态方程的基础上的，当然仅适用于线性系统。图4.2和图4.3是从传统数值积分法的几何解释借鉴过来的，显然仅对一阶系统才是正确的，因为高阶系统无法用平面图来描述，但道理是一样的。 46 关于数字仿真计算的稳定性分析连续系统的数字仿真，实质上就是将给定的微分方程变换为差分方程，然后将该差分方程在数字计算机上求解。尽管原来系统的差分方程是稳定的，但是由于变换的方法及所取的计算步长不同，最后求得的差分方程解可能是不稳定的。这样，在数字计算机上求得的数值解也将是发散的。下面对这个问题做一些专门讨论。差分方程的解与微分方程的解类似，均可分为特解及一般解。因此，在讨论差分方程稳定性时，可以只研究与自由运动有关的这一部分。现考察一阶常系数微分方程式 （4-38）如果用数值积分法中的欧拉公式，可得到微分方程的数值解 （4-39） 为了研究这个解的稳定性，考查齐次解 （4-40）当时，式（4-40）的差分方程是不稳定的。表4-2列出了时的解序列。表4.2 差分方程式 在 、 时的不稳定响应Ky(k+1)0-2142-8316但是当，差分方程（4-40）式是稳定的。表4-3列出了时的解序列。表4.3 差分方程式 在 、 时的稳定响应Ky(k+1)0-0.910.812-0.72930.6561由此可见，采用古典数值积分法中的欧拉法仿真时，如果步长超过某一个限度，数字解将是不稳定的。同理，也可以证明显型古典数值积分公式都存在数值解不稳定问题，它们对步长都是有限制的。但是差分方程的阶次越高，允许的计算步长越大，请参阅文献56。对于方程式（4-38），如果采用离散相似法，在输入处加入采样开关及零阶再现环节，可得到系统的差分方程 （4-41）为了研究系统解的稳定性，考查齐次解 （4-42）不难看出，无论取何值，都不会大于1。即对任意的计算步长，解序列都是稳定的。所以，从仿真系统稳定性角度看，状态方程离散化方法要比数值积分方法好，其特点是，如果原来的微分方程是稳定的，则通过这种方法所得到的差分方程也是稳定的。因此，许多快速仿真算法都是由离散相似法推导出来的。4.7连续系统数字仿真小结对于数字计算机来说，非线性环节的仿真是既简单又准确的，这在3.5节中已作了详细的说明。线性系统则不一样。由于数字计算机不能对微分方程直接求解，必须将微分方程化成与其近似的差分方程才能求得数值解。即差分方程数值解在其采样点上的值和原微分方程在同一时刻的解值近似相等。因此，线性系统的仿真问题，实际上就是将描述该系统的微分方程（或状态方程）化成相似的差分方程（或离散状态方程）的问题。我们把后者称作仿真模型。本书介绍了两种线性系统仿真方法，即离散相似法和数值积分法