连续系统的数字仿真

单击此处编辑母版标题样式,*,*,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,内蒙古科技大学,第三章 连续系统的数字仿真,本章内容,(1),熟悉在数字计算机仿真技术中常用的几种数值积分法，,特别是四阶龙格,-,库塔法；,(2),典型环节及其系数矩阵的确定；,(3),各连接矩阵的确定；,(4),利用,MATLAB,在四阶龙格,-,库塔法的基础上，对以状态,空间表达式和方框图描述的连续系统进行仿真；,(5),了解以增广矩阵法为基础的连续系统的快速仿真方法。,1,用数字计算机来仿真或模拟一个连续控制系统的目的就是求解系统的数学模型。由控制理论知，一个,n,阶连续系统可以被描述成由,n,个积分器组成的模拟结构图。因此利用数字计算机来进行连续系统的仿真，从本质上讲就是要在数字计算机上构造出,n,个数字积分器，也就是让数字计算机进行,n,次数值积分运算。可见，连续系统数字仿真中的最基本的算法是数值积分算法。,2,3.1,数值积分法,连续系统通常把数学模型化为状态空间表达式，为了对,n,阶连续系统在数字计算机上仿真及求解，就要采用数值积分法来求解系统数学模型中的,n,个一阶微分方程。,设,n,阶连续系统由以下,n,个一阶微分方程组成,（,3-1,）,所谓数值积分法，就是要逐个求出区间,a ,b,内若干个离散点,a,t,0,t,1,t,0,时，,x,(,t,),是未知的，因此式（,3-2,）右端的积分是求不出的。为了解决这个问题，把积分间隔取得足够小，使得在,t,k,与,t,k,+1,之间的,f,(,t,x,(,t,),可以近似看作常数,f,(,t,k,x,(,t,k,),，,这样便得到用矩形公式积分的近似公式,或简化为,这就是欧拉公式。,5,以,x,(,t,0,)=,x,0,作为初始值，应用欧拉公式，就可以一步步求出每一时刻,t,k,的,x,k,值，,即,k,=0,x,1,x,0,+,f,(,t,0,x,0,),h,k,=1,x,2,x,1,+,f(t,1,x,1,),h,k,=,n,-1,x,n,x,n-1,+,f,(,t,n-1,x,n-1,),h,这样式,(3-1),的解,x,(,t,),就求出来了。欧拉法的计算虽然比较简单，但精度较低。图,3-1,为欧拉法的几何解释。,6,3.1.2,梯形法,由上可知欧拉公式中的积分是用矩形面积,f,(,t,k,x,k,),h,来近似的。,由图,3-2,知，用矩形面积,t,k,abt,k,+1,代替积分，其误差就是图中阴影部分。为了提高精度现用梯形面积,t,k,act,k,+1,来代替积分，即,于是可得梯形法的计算公式为,7,由于上式右边包含未知量,x,k,+1,所以每一步都必须通过迭代求解，每一步迭代的初值,x,k,+1,(0),通常采用欧拉公式来计算，因此梯形法每一步迭代公式为,3-3,）,式中 迭代次数,R,=0,1,2,8,3.1.3,预估校正法,虽然梯形法比欧拉法精确，但是由于每一步都要进行多次叠代，计算量大，为了简化计算，有时只对式（,3-3,）进行一次叠代就可以了，因此可得,通常称这类方法为预估校正方法。它首先根据欧拉公式计算出,x,k+,1,的预估值,x,k+,1,(0),，,然后再对它进行校正，以得到更准确的近似值,x,k+,1,(1),。,9,3.1.4,龙格库塔法,根据泰勒级数将式（,3-1,）在,t,k+,1,=,t,k,+h,时刻的解,x,k+,1,=,x,(,t,k,+h,),在,t,k,附近展开，有,3-5,）,可以看出，提高截断误差的阶次，便可提高其精度，但是由于计算各阶导数相当麻烦，所以直接采用泰勒级数公式是不适用的，为了解决提高精度问题，龙格和库塔两人先后提出了间接使用泰勒级数公式的方法，即用函数值,f,(,t,x,),的线性组合来代替,f,(,t,x,),的导数，然后按泰勒公式确定其中的系数,这样既能避免计算,f,(,t,x,),的导数，又可以提高数值计算精度，其方法如下。,10,因,故式（,3-5,）可写成,（,3-6,）,为了避免计算式（,3-6,）中的各阶导数项，可令,x,k+,1,由以下多项式表示。,（,3-7,）,11,式中,a,m,为待定因子，,v,为使用,f,函数值的个数，,k,m,满足下列方程,（,3-8,）,即：,将式（,3-7,）展开成,h,的幂级数并与微分方程式（,3-1,）精确解式（,3-6,）逐项比较,便可求得式（,3-7,）和式（,3-8,）中的系数,a,m,b,mj,和,c,m,等。,12,现以,v,=2,为例，来说明这些参数的确定方法。设,v,=2,，,则有,（,3-9,）,将,k,1,和,k,2,在同一点（,t,k,x,k,）,上用二元函数展开为,13,将,k,1,和,k,2,代入式（,3-9,）整理后可得,（,3-10,）,将上式与式（,3-6,）逐项进行比较，可得以下关系式,若取,则,14,于是可得,(3-11),由于式（,3-11,）只取到泰勒级数展开式的,h,2,项，故称这种方法为两阶龙格库塔法，其截断误差为,0(,h,3,),。,15,同理当,v,=4,时，仿照上述方法可得如下四阶龙格,-,库塔公式,16,通过上述龙格,-,库塔法的介绍，可以把以上介绍的几种数值积分法统一起来，它们都是基于在初值附近展开成泰勒级数的原理，所不同的是取泰勒级数多少项。欧拉公式仅取到,h,项，梯形法与二阶龙格库塔法相同均取到,h,2,项，四阶龙格库塔法取到,h,4,项。从理论上讲，取得的项数愈多，计算精度愈高，但计算量愈大，愈复杂，计算误差也将增加，因此要适当的选择。目前在数字仿真中，最常用的是四阶龙格库塔法，其截断误差为,(,h,5,),已能满足仿真精度的要求。,17,3.2,连续系统的数字仿真程序,若系统的状态空间表达式为,（,3-14,）,（,3-15,）,其中,A,:,n,n; b,:,n,1,; C,:1,n,18,假设在仿真中，数值积分法采用四阶龙格库塔方法，因对于,n,阶系统，其状态方程式（,3-14,）可写成以下,n,个一阶微分方程,(3-16),19,故根据式（,3-12,）可得求解以上一阶微分方程组的四阶龙格库塔公式如下,式中,x,i,k,为,t=,t,k,时刻的,x,i,值，,x,i,k+,1,为,t=,t,k,+h,时刻的,x,i,值。,20,令,21,则式（,3-17,）可写成如下矩阵的形式,根据式（,3-15,）可得,t=t,k,+1,时刻的输出,22,23,例,3-1,假设单变量系统如下图所示。,试根据四阶龙格库塔法，求输出量,y,的动态响应。,24,解,仿真程序如下,Example3.1,取仿真时间：,T,f,=5,；,计算步长：,h,=0.02,在,MATLAB,环境下执行以上程序可得如图,3-6,所示仿真曲线。,25,图,3-6,26,3.3,面向系统结构图的仿真,这种方法与上节介绍的方法相比，有以下几个主要优点：,）,便于研究各环节参数对系统的影响；,）,可以得到每个环节的动态响应；,）,可对多变量系统进行仿真。,下面具体介绍面向结构图的仿真方法。,27,3.3.1,典型环节的确定,一个控制系统可能由各种各样的环节所组成，但比较常见环节有：,1,）,比例环节：,2,）,积分环节：,3,）,比例积分环节：,4,）惯性环节：,5,）,超前滞后环节：,6,）二阶振荡环节：,28,为了编制比较简单而且通用的仿真程序必须恰当选择仿真环节。在这里选用图,3-7,所示的典型环节作为仿真环节，即,式中,u,为典型环节的输入，,x,为典型环节的输出。,利用这个典型环节，只要改变,a ,b ,c,和,d,参数的值，便可分别表示以上所述的各一阶环节，至于二阶振荡环节，则可用,两个一阶环节等效连,接得到，如图,3-8,所示。,29,3.3.2,连接矩阵,一个控制系统用典型环节来描述时，必须用连接矩阵把各个典型环节连接起来。所谓连接矩阵，就是用矩阵的形式表示各个典型环节之间的关系。下面介绍连接矩阵的建立方法,假设多输入多输出系统的结构图如图,3-9,所示。,30,图中带数字的方框表示典型环节，,表示比例系数。,31,由图可得各环节的输入与各环节的输出间的关系以及系统的输出与环节的输出间关系分别为,和,32,写成矩阵形式,和,33,或写成,(3-20),定义式中的,W,，,W,0,和,W,c,阵为连接矩阵，,W,反映了各典型环节输入输出间的连接关系，,W,0,反映了系统的参考输入与各环节输入间的连接关系，,W,c,反映了系统的输出与各环节输出间的关系。,34,一般也将系统中各典型环节的系数写成如下矩阵的形式（假设系统由,n,个典型环节组成）,(3-21),35,3.3.3,确定系统的状态方程,典型环节和连接矩阵确定后，便可求得系统的状态空间表达式，推导过程如下。,假设系统由,n,个典型环节组成，则根据典型环节的传递函数有,（,i,=1,2,）,即,36,写成矩阵形,（,3-22,）,式中：,37,将式（,3-20,）中的上式进行拉氏变换后代入式（,3-22,）中可得,对上式两边取拉氏反变换得,（,3-23,）,若参考输入向量,r,=,r,1,r,2,r,m,T,中的,r,1,r,2,r,m,均为阶跃函数，则上式可简化为,（,3-24,）,38,令,则式（,3-24,）可写成,若,H,的逆存在，则有,再令,可得 （,3-25,）,上式即为闭环系统的状态方程，它是一个典型的状态方程，利用前面介绍的求解方法可方便地求出各典型环节的输出响应，最后根据式（,3-20,）中的第二式便可求出系统的输出响应。,39,3.3.4,面向结构图的数字仿真程序,面向结构图的数字仿真程序框图如图,-10,所示，其程序清单通过下例给出。,40,例,3-2,假设某一系统由四个典型环节组成，如图,3-11,所示。求输出量,y,的动态响应。,41,解,由图可得各环节的输入与输出以及系统的输出与环节的输出间关系为,42,根据以上两式和各典型环节的系数值，可得如下连接矩阵和系数矩阵,43,取仿真时间,:,T,f,=10,；,计算步长,:,h,=0.05,在,MATLAB,环境下执行以上程序可得如图,3-12,所示仿真曲线。,仿真程序如下,Example3.2,44,图,3-12,45,3.4,连续系统的快速仿真,前面介绍过的两种连续系统的数字仿真方法，当系统比较复杂并要求满足较高的计算精度时，计算工作量较大，计算速度较慢，有时不能满足实时仿真的要求，为了解决这个问题，下面介绍一种连续系统的快速数字仿真方法,-,增广矩阵法。,46,3.4.1,增广矩阵法的基本原理,设连续系统的状态方程为,（,3-26,）,则它的解为,将,e,At,展开成泰勒级数，即,则有,47,可以证明，如果取,e,At,的泰勒级数的前五项，则式（,3-27,）的计算精度与四阶龙格库塔法相同，设计算步长为,T,则式（,3-27,）可写成,上式括号中只是矩阵的相乘及加法，仿真计算十分简单，因此可大大加快数字仿真的计算速度，如果将矩阵的相乘在数字仿真前算好，则数字仿真的速度将更快。,48,如果要求解的状态方程为非齐次方程,为了能对其应用上述的快速计算方法，就要将控制量,Bu,(,t,),项增广到状态变量中去，将其转换为齐次方程，这就是增广矩阵法的基本原理。当,u,(,t,),是一些典型函数时，增广矩阵是很容易实现的。下面介绍三种典型输入的增广矩阵。,49,3.4.2,三种典型输入函数的增广矩阵,假设,n,阶连续系统的状态空间表达式为,（,3-29,）,1.,当输入信号为阶跃函数时,即,u,(,t,)=,r,0,定义第,n+1,个状态变量为,x,n+,1,(,t,),=u,(,t,),=r,0,则 （,3-30,）,50,将式（,3-30,）增广到式（,3-29,）中可得增广后的状态空间表达式,增广后的系统矩阵 称为增广矩阵,。,51,2.,当输入信号为斜坡输入时，即,定义,则,52,因此增广后的状态空间表达式为,53,3.,当输入信号为抛物线时，即,定义,则,54,因此增广后的状态空间表达式为,55,r=10;,P=0.1 1 0.5 1;0 1 1 0;2 1 2 0;10 1 10 0;,W=0 0 0 -1;1 0 0 0 ;0 1 0 0;0 0 1 0;,W0=1;0;0;0;,wc,=0 0 0 1;,tf,=input(,仿真时间,tf,=);,h=input(,计算步长,h=);,A1=diag(P(:,1);,B1=diag(P(:,2);,C1=diag(P(:,3);,D1=diag(P(:,4);,H=B1-D1*W;,H=,inv(H,);,Q=C1*W-A1;,A=H*Q;,B=H*C1*W0;,56,x=zeros(length(A),1);,y=0;,t=0;,for i=1:tf/h,k1=A*,x+B,*r;,k2=A*(,x+h,*k1/2)+B*r;,k3=A*(,x+h,*k2/2)+B*r;,k4=A*(,x+h,*k3)+B*r;,x=,x+h,*(k1+2*k2+2*k3+k4)/6;,y=,y;wc,*,x;t,=,t;t(i)+h,;,end,plot(t,y,),57,r=20.4;,num0=1.875e6 1.562e6;,den0=1 54 204.2 213.8 63.5;,numh,=0.002;denh=1;,num,den,=feedback(num0,den0,numh,denh);,A,b,C,D,=tf2ss(num,den);,Tf,=input(,仿真时间,Tf,=);,h=input(,计算步长,h=);,x=zeros(length(A),1);y=0;t=0;,58,for i=1:Tf/h,k1=A*,x+b,*r;,k2=A*(,x+h,*k1/2)+b*r;,k3=A*(,x+h,*k2/2)+b*r;,k4=A*(,x+h,*k3)+b*r;,x=,x+h,*(k1+2*k2+2*k3+k4)/6;,y=,y;C,*,x;t,=,t;t(i)+h,;,end,plot(t,y,),59,