第六讲-连续系统仿真概述VIP

计算机仿真技术,0,Email:swpikyc,计算机仿真技术,主讲：况雨春,第六讲连续系统仿真概述,计算机仿真技术,1,作业讲解,x,y=meshgrid(-8:.5:8);R=sqrt(x.2+y.2)+eps;z=sin(R)./R;surf(x,y,z,FaceColor,red,EdgeColor,none);camlightleft,计算机仿真技术,2,连续系统的动态特性，一般可用微分方程、状态方程、传递函数以及系统的结构图等形式来描述。步骤：建立系统的数学模型-仿真本质：对微分方程求数值解中心问题：选择一种仿真算法将系统模型转化为能在计算机上运行的离散模型，即仿真模型。两大类算法：欧拉法、龙格库塔法离散相似法：离散状态法,一、连续系统的数字仿真,计算机仿真技术,3,1、建立数学模型2、建立仿真模型3、编写仿真程序,二、连续系统的数字仿真过程,计算机仿真技术,4,连续系统模型描述,系统状态变化在时间上是连续的，可以用方程式（常微分方程、偏微分方程、差分方程）描述系统模型。数学模型常用的三种形式：微分方程、传递函数、状态方程。系统的输入量u(t)，输出量y(t)，系统的内部状态变量x(t)都是时间的连续函数,计算机仿真技术,5,微分方程数学模型,一般形式,n为系统的阶次，ai为系统的结构参数，ci为输入函数的结构参数，它们均为实常数。y为输出量,u为输入量系统没有输入函数则上式变为齐次方程一阶系统：追击敌机问题二阶系统：汽车轮子悬置系统,计算机仿真技术,6,追击敌机问题,敌机在100KM高处，速度Vb=20km/min我机追击Vf=30km/min每分钟改变一次航向，使机头对准敌机，相距10公里可以攻击。12min内，敌机逃脱，则任务失败。,计算机仿真技术,7,XbVb*tYb=100Xf=Vf*cosYf=Vf*sintan=(100-Yf)/(Xb-Xf)描述方程龙格库塔法求解位置求解我机角度及位置判断距离时间end,计算机仿真技术,8,二阶系统：汽车轮子悬置系统,y(t)+p1y(t)+y(t)=1.0,主程序（系统参数及仿真参数）,运行子程序,微分方程子程序,积分子程序,判断,输出子程序,N,Y,计算机仿真技术,9,传递函数数学模型,表达系统输入量转化为输出量的传递函数关系若系统的初始条件为零，即系统在t=0时已处于一个稳定状态，那么对上式两边取拉氏变换后可得：,Y(s)=LY(t)，U(s)=LU(t)稍加整理，并记：,计算机仿真技术,10,状态方程数学模型,以上模型都只描述了系统输入与输出之间的关系，而没有描述系统内部的情况，所以这些模型称为外部模型。从仿真的角度来看，为在计算机上对系统的模型进行试验，就要在计算机上实现这个模型。有时，仅仅实现系统输入与输出之间的关系是不够的，还必须实现模型的内部变量状态变量，因此仿真要求采用系统内部模型，即状态方程模型。,计算机仿真技术,11,计算机仿真技术,12,计算机仿真技术,13,状态方程模型的一般形式为：,第一个方程称为状态方程，第二个方程称为输出方程。其中A是n*n维矩阵，B是n*1维矩阵，C是1*n维矩阵。对微分方程式的单输入/单输出的n阶系统，易于将其转换为上述形式的状态方程引进n个内部状态变量x1，x2，xn，作用函数为单输入u，输出变量为单输出y，仿真时，必须将系统的外部模型转换成内部模型，也就是建立与输入/输出特性等价的状态方程。,计算机仿真技术,14,y(t)+p1y(t)+y(t)=1.0,X1=yX2=yX3=y”X1=X2y”(t)=X2=-X1-p1X2+1,计算机仿真技术,15,数值积分法系统仿真,数值积分法,连续系统仿真首先遇到的问题：解决数字计算机在数值及时间上的离散性与被仿真系统数值及时间上的连续性这一基本问题。连续系统仿真，从本质上是对原连续系统从时间、数值两个方面对原系统进行离散化，并选择合适的数值计算方法来近似积分运算，由此得到的离散模型来近似原连续模型。如何保证离散模型的计算结果从原理上确能代表原系统的行为，这是连续系统数字仿真首先必须解决的问题。,计算机仿真技术,16,设系统模型为：，其中u(t)为输入变量，y(t)为系统变量；令仿真时间间隔为h，离散化后的输入变量为，系统变量为，其中tk表示t=kh。如，即，（对所有k=0,1,2,），则可认为两模型等价，这称为相似原理,计算机仿真技术,17,连续系统数字仿真中离散化最基本的算法是数值积分算法。对于形如的系统，已知系统变量的初始条件，现在要求y随时间变化的过程y(t)。计算过程可以这样考虑：首先求出初始点的，微分方程可以写作：,欧拉法用矩形面积近似表示积分结果，也就是当t=t1时，y(t1)的近似值为y1：,计算机仿真技术,18,所以，对任意时刻tk+1，有：,为进一步提高计算精度，人们提出了“梯形法”。梯形法近似积分形式如下所示，（没法自启动）,采用这种积分方法最简单的预报校正方法是用欧拉法估计初值，用梯形法校正，即：,计算机仿真技术,19,四阶龙格库塔法,计算机仿真技术,20,计算机仿真技术,21,龙格-库塔法的精度取决于步长h的大小及方法的阶次。许多计算实例表明：为达到相同的精度，四阶方法的h可以比二阶方法的h大10倍，而四阶方法的每步计算量仅比二阶方法大一倍，所以总的计算量仍比二阶方法小。正是由于上述原因，一般系统进行数字仿真常用四阶龙格-库塔公式。值得指出的是：高于四阶的方法由于每步计算量将增加较多，而精度提高不快，因此使用得也比较少。,计算机仿真技术,22,如果我们将步长取得十分小，那么欧拉公式也可以获得很高的计算精度吗？从理论上讲是这样，但是实际上由于计算机字长有限，在计算中存在着舍入误差，它与计算次数成正比。计算步长很小，则一个系统的过渡过程的计算就会增加许多，舍入误差就会十分明显地表现出来，这样一来，很难保证有较高的精度。,计算机仿真技术,23,在MATLAB里给出了许多数值积分仿真算法：欧拉法：euler只用来对仿真结果进行验证龙格库塔法：ode23,ode45系统为高度非线性或者不连续用法：t,x=solver(model,t0,tf,x0,tol,minstep,maxstep)初始时间t0；时间终点：tf；初始状态：x0后面依次：容许误差、最小步长、最大步长。,计算机仿真技术,24,面向方程的龙格库塔法仿真举例,为了加深对数值积分法的理解，以便能使用仿真程序来解决一两个实际问题，本小节将给出面向方程的龙格库塔法仿真实例。“面向方程”是指仿真对象是以一阶微分方程组形式描述的系统，如下式所示：,是n维状态向量是k维参数向量是m维控制向量是n维函数向量是给定的初始状态,计算机仿真技术,25,下面对仿真对象作一说明。设卫星在空中运行的运动方程为：,其中是k重力系数(k=401408km3/s)。卫星轨道采用极坐标表示，通过仿真，研究发射速度对卫星轨道的影响。,计算机仿真技术,26,这是一个二阶微分方程组，为此，首先要将其转换成一阶微分方程组。若设，则有：,四个一阶微分方程，有4个状态变量y(1)y(2).。,计算机仿真技术,27,今希望用直角坐标输出，故还引入两个定义变量，这样又得到两个代数方程：,根据卫星的发射速度，可以建立起方程组的初始条件：y(1)0=6400km(卫星到地心的距离，即为地球之半径)，仿真研究的目的是要计算出当发射速度为8km/s，10km/s及12km/s时的卫星轨道。,计算机仿真技术,28,已知y(4)=v/y(1)，由v及y(1)0可求出y(4)0，它们分别是：0.00125，0.0015625，0.001875(1/s)，为了仿真该系统，用户要输入系统的运动方程，即要将上面列的两个方程组写到程序中去。这部分写在系统模型输入(modsub)程序块中；比如例中要求输出y(21),y(22),y(1),y(2)。则A(1)，A(2)，A(3)，A(4)四个数组都要用到这部分写在输出转换(repsub)程序块中；状态变量的初始数据及要改变的参量数值，放在主程序中。此系统没有参数变量。,计算机仿真技术,29,计算机仿真技术,30,下面就是为仿真卫星轨道系统由用户写的一段仿真源程序。在系统模型输入程序块(modsub)中写入以下程序：Y(21)=Y(1)*cos(Y(2)Y(22)=Y(1)*sin(Y(2)G(1)=Y(3)G(2)=Y(4)G(3)=-401408/(Y(1)*Y(1)+Y(1)*Y(4)*Y(4)G(4)=-2*Y(3)*Y(4)/Y(1),计算机仿真技术,31,在输出转换(repsub)程序块中写入以下程序：A(1)=Y(21)A(2)=Y(22)A(3)=Y(1)A(4)=Y(2),如果需要，可以在输出打印(outsut)程序块中加入Y(4)0=Y(4)的打印语句。程序写好后则可按以下步骤进行操作：(1)将用户写好的仿真源程序输入到计算机中。(2)运行按计算机提问输入必要的数据。计算机执行并在终端上或打印机上输出运行结果。,计算机仿真技术,32,InputtotalsimulationtimeT1,stept2=10000,200InputtheorderofsystemN1=4InputthenumberofoutputpointsN3=51symbol191fSymbols10.5InputthenumberofrunsJ1=3Inputinitialvaluesofstastevariables=6400，0，0，00125Inputthevaluesofparameters=0,计算机仿真技术,33,面向结构图的系统仿真,工程上常常将系统描述为结构图的形式，因此，工程技术人员更习惯面向结构图的仿真方法。这里介绍一个面向结构图的线性系统模型变换方法，以便通过计算机自动将结构图模型转换成一阶微分方程组，,典型环节的选择结构图是由各种典型环节构成的通常一个系统中比较常见的动态环节有：积分环节：比例积分环节：,计算机仿真技术,34,惯性环节：一阶超前(或滞后)环节二阶振荡环节：,选择什麽样的环节做为模型描述的基础，这是模型结构变换十分重要的一步。选择典型环节的原则是：要有典型性，即由它可组成各种动态环节，另外由它组成系统简便，由计算机将它转换成系统的微分方程组容易实现。有人建议用的一阶超前(或滞后)环节来作为典型环节。这种典型环节可以很容易的表示上述常见的环节，,计算机仿真技术,35,假定有一个如图所示的系统，其中用方框框起来并标上1，2，符号的都是典型环节,计算机仿真技术,36,可以写出系统的动态方程为：(A+BS)Y=(C+DS)uu=WY+W0y0在方程中A、B、C、D阵中的元素反映各环节的参数，为对角矩阵。该方程描述了各环节的输入、输出及参数之间的关系，,W阵和W0阵均为系统的连接矩阵，描述了系统各环节之间连接关系W称为系统的连接矩阵，它描述了系统内部各环节连接情况，W0称为外部输入的连接矩阵，它描述了外部输入对系统的作用情况。,计算机仿真技术,37,系统方程转换系解决结构图如何转换成一组一阶线性微分方程。,计算机仿真技术,38,离散相似法将一个连续系统进行离散化处理，然后求得与它等价的离散模型。由于连续系统的模型可以用传递函数来表示，也可以用状态空间模型来表示，因此，与连续系统等价的离散模型可以通过两个途径获得，其一是对传递函数作离散化处理得离散传递函数(或脉冲传递函数)，称为频域离散相似模型。其二是基于状态方程离散化，得到时域离散相似模型。假设有一个连续系统，它由以下状态方程描述：,计算机仿真技术,39,其中，u(t)是系统输入；u(k)是加虚拟采样开关后，在kT时刻系统输入；x(k)是加虚拟采样开关后在kT时刻系统输出；等价的连续信号对该连续系统进行离散化处理后可以得到系统离散相似模型,计算机仿真技术,40,在MATLAB中数学模型的表示,线性系统理论中常用的数学模型有微分方程、传递函数、状态空间表达式等,而这些模型之间又有着某些内在的等效关系。MATLAB主要使用传递函数和状态空间表达式来描述线性时不变系统(LinearTimeInvariant简记为LTI)。,计算机仿真技术,41,1.5.1传递函数,单输入单输出线性连续系统的传递函数为,其中mn。G(s)的分子多项式的根称为系统的零点,分母多项式的根称为系统的极点。令分母多项式等于零,得系统的特征方程:,D(s)=a0sn+a1sn1+an1s+an=0,计算机仿真技术,42,系统的传递函数在MATLAB下可由其分子和分母多项式唯一地确定出来,其格式为,sys=tf(num,den),其中num为分子多项式，den为分母多项式,num=b0,b1,b2,bm；den=a0,a1,a2,an；,计算机仿真技术,43,对于其它复杂的表达式,如,可由下列语句来输入,num=conv(1,1,conv(1,2,6,1,2,6);den=conv(1,0,0,conv(1,3,1,2,3,4);G=tf(num,den),Transferfunction:,计算机仿真技术,44,1.5.2传递函数的特征根及零极点图,传递函数G(s)输入之后,分别对分子和分母多项式作因式分解,则可求出系统的零极点,MATLAB提供了多项式求根函数roots(),其调用格式为,roots(p),其中p为多项式。,计算机仿真技术,45,例如,多项式p(s)=s3+3s2+4,p=1,3,0,4;%p(s)=s3+3s2+4r=roots(p)%p(s)=0的根r=-3.35330.1777+1.0773i0.1777-1.0773i,反过来,若已知特征多项式的特征根,可调用MATLAB中的poly()函数,来求得多项式降幂排列时各项的系数,如上例,poly(r)p=1.00003.00000.00004.0000,计算机仿真技术,46,而polyval函数用来求取给定变量值时多项式的值,其调用格式为,polyval(p,a)其中p为多项式;a为给定变量值,例如,求n(s)=(3s2+2s+1)(s+4)在s=5时值：,n=conv(3,2,1,1,4);value=polyval(n,-5)value=66,计算机仿真技术,47,p,z=pzmap(num,den)其中,p传递函数G(s)=numden的极点z传递函数G(s)=numden的零点例如,传递函数,传递函数在复平面上的零极点图,采用pzmap()函数来完成,零极点图上,零点用“。”表示,极点用“”表示。其调用格式为,计算机仿真技术,48,用MATLAB求出G(s)的零极点,H(s)的多项式形式,及G(s)H(s)的零极点图,numg=6,0,1;deng=1,3,3,1;z=roots(numg)z=0+0.4082i00.4082i%G(s)的零点p=roots(deng)p=1.0000+0.0000i1.0000+0.0000i%G(s)的极点1.0000+0.0000i,计算机仿真技术,49,n1=1,1;n2=1,2;d1=1,2*i;d2=1,-2*i;d3=1,3;numh=conv(n1,n2);denh=conv(d1,conv(d2,d3);printsys(numh,denh),numh/denh=,%H(s)表达式,pzmap(num,den)%零极点图title(pole-zeroMap),计算机仿真技术,50,零极点图如图所示：,计算机仿真技术,51,1.5.3控制系统的方框图模型,若已知控制系统的方框图,使用MATLAB函数可实现方框图转换。,1.串联如图所示G1(s)和G2(s)相串联,在MATLAB中可用串联函数series()来求G1(s)G2(s),其调用格式为num,den=series(num1,den1,num2,den2)其中：,计算机仿真技术,52,2.并联如图所示G1(s)和G2(s)相并联,可由MATLAB的并联函数parallel()来实现,其调用格式为,num,den=parallel(num1,den1,num2,den2),其中：,计算机仿真技术,53,3.反馈,反馈连接如图所示。使用MATLAB中的feedback()函数来实现反馈连接,其调用格式为,num,den=feedback(numg,deng,numh,denh,sign),式中：,sign为反馈极性,若为正反馈其为1,若为负反馈其为1或缺省。,计算机仿真技术,54,例如G(s)=,H(s)=,负反馈连接。,numg=1,1;deng=1,2;numh=1;denh=1,0;num,den=feedback(numg,deng,numh,denh,1);printsys(num,den),num/den=,计算机仿真技术,55,MATLAB中的函数series,parallel和feedback可用来简化多回路方框图。另外,对于单位反馈系统,MATLAB可调用cloop()函数求闭环传递函数,其调用格式为,num,den=cloop(num1,den1,sign),计算机仿真技术,56,1.5.4控制系统的零极点模型,传递函数可以是时间常数形式,也可以是零极点形式,零极点形式是分别对原系统传递函数的分子和分母进行因式分解得到的。MATLAB控制系统工具箱提供了零极点模型与时间常数模型之间的转换函数,其调用格式分别为,z,p,k=tf2zp(num,den)num,den=zp2tf(z,p,k),其中第一个函数可将传递函数模型转换成零极点表示形式,而第二个函数可将零极点表示方式转换成传递函数模型。,计算机仿真技术,57,例如G(s)=,用MATLAB语句表示：,num=12241220;den=24622;z,p,k=tf2zp(num,den)z=1.92940.03530.9287i0.03530.9287i,计算机仿真技术,58,p=0.95671.2272i0.95671.2272i0.04330.6412i0.04330.6412ik=6,即变换后的零极点模型为G(s)=,计算机仿真技术,59,可以验证MATLAB的转换函数,调用zp2tf()函数将得到原传递函数模型。num,den=zp2tf(z,p,k)num=06.000012.00006.000010.0000den=1.00002.00003.00001.00001.0000,即,计算机仿真技术,60,1.5.5状态空间表达式,状态空间表达式是描述系统特性的又一种数学模型,它由状态方程和输出方程构成,即x(t)=Ax(t)+Bu(t)y(t)=Cx(t)+Du(t),式中x(t)Rn称为状态向量,n为系统阶次;,ARnn称为系统矩阵;BRnp称为控制矩阵,p为输入量个数;CRqn称为输出矩阵;DRqp称为连接矩阵,q为输出量个数。,计算机仿真技术,61,在一般情况下,控制系统的状态空间表达式项简记为(A,B,C,D)。,例如：设一个双输入双输出系统的状态空间表达式为,计算机仿真技术,62,系统模型可由MATLAB命令直观地表示：,A=1,2,4;3,2,6;0,1,5B=4,6;2,2;0,2C=0,0,1;0,2,0D=zeros(2,2),MATLAB的控制系统工具箱提供了由状态空间表达式转换成传递函数或由传递函数转换成状态空间表达式的转换函数ss2tf()和tf2ss()。其调用格式为,num,den=ss2tf(A,B,C,D,iu),计算机仿真技术,63,反过来,若已知系统的传递函数,求取系统状态空间表达式的调用格式为A,B,C,D=tf2ss(num,den),例如系统的传递函数为,系统的状态空间表达式为,num=1,2,3;den=1,3,6,1;A,B,C,D=tf2ss(num,den),计算机仿真技术,64,A=-3-6-1100010B=100C=123D=0,返回,计算机仿真技术,65,1.6用MATLAB和SIMULINK进行瞬态响应分析,1.6.1单位脉冲响应当输入信号为单位脉冲函数(t)时，系统输出为单位脉冲响应，MATLAB中求取脉冲响应的函数为impulse()，其调用格式为y，x，t=impulse(num，den，t)或impulse(num，den)式中G(s)=num/den；t为仿真时间；y为时间t的输出响应；x为时间t的状态响应。,计算机仿真技术,66,例试求下列系统的单位脉冲响应MATLAB命令为：t=0：0.1：40；num=1；den=1，0.3，1；impulse(num，den，t)；grid；title(Unit-impulseResponseofG(s)=1/(s2+0.3s+1)其响应结果如图所示。,计算机仿真技术,67,例系统传递函数为求取其单位脉冲响应的MATLAB命令为t=0：0.1：10；num=1；den=1，1，1；y，x，t=impulse(num，den，t)plot(t，y)；gridxlabel(t)；ylable(y)；其响应结果如图所示。,计算机仿真技术,68,1.6.2单位阶跃响应,当输入为单位阶跃信号时，系统的输出为单位阶跃响应，在MATLAB中可用step()函数实现，其调用格式为y，x，t=step(num，den，t)或step(num，den),计算机仿真技术,69,例求系统传递函数为num=1；den=1，0.5，1；t=0：0.1：10；y，x，t=step(num，den，t)；plot(t，y)；grid；xlabel(Timesect)；ylabel(y)响应曲线如图所示,图单位阶跃响应,计算机仿真技术,70,1.6.3斜坡响应,在MATLABA中没有斜坡响应命令，因此，需要利用阶跃响应命令来求斜坡响应。根据单位斜坡响应输入是单位阶跃输入的积分。当求传递函数为的斜坡响应时，可先用除得，再利用阶跃响应命令即可求得斜坡响应。,计算机仿真技术,71,例已知闭环系统传递函数对单位斜坡输入则,计算机仿真技术,72,系统单位斜坡响应的MATLAB命令：num=1；den=1，0.3，1，0；t=0：0.1：10；c=step(num，den，t)；plot(t，c)；grid；xlabel(tsec)；ylabel(InputandOutput)其响应结果如图所示。,计算机仿真技术,73,1.6.4任意函数作用下系统的响应,用线性仿真函数lsim来实现，其调用格式为y，x=lsim(num，den，u，t)式中；y(t)为系统输出响应；x(t)为系统状态响应；u为系统输入信号；t为仿真时间。,计算机仿真技术,74,例反馈系统如图(a)所示，系统输入信号为图(b)所示的三角波，求取系统输出响应。,计算机仿真技术,75,MATLAB实现指令numg=10，20；deng=1，10，0；num，den