1、系统仿真技术 第1章 连续系统模型描述,陈无畏 合肥工业大学机械与汽车工程学院,1.1 连续系统模型描述,连续系统-系统状态变化在时间上是连续的，可以用方程式（常微分方程、偏微分方程、差分方程）描述系统模型。,一个系统可以定义成如下集合结构：T：时间基，描述系统变化的时间坐标。 T为整数则称为离散时间系统， T为实数则称为连续时间系统。 X：输入集， 代表外部环境对系统的作用。 X被定义为 ,其中 ，X即代表n个实值的输入变量。 ：输入段集，描述某个时间间隔内输入模式，是(X,T)的子集。 Q：内部状态集，是系统内部结构建模的核心。,：状态转移函数，定义系统内部状态是如何变化的。它是映射：其含
2、义：若系统在 时刻处于状态q，并施加 一个输入段 ，则 表示系统 处于 状态。 ：输出函数，它是映射： 输出函数给出了一个输出段集。 Y：输出段集，系统通过它作用于环境。,连续系统数学模型典型形式,常微分方程 传递函数 状态空间描述 权函数（脉冲过渡函数）,1.1.1常微分方程-输入/输出水平,（1） 其中n为系统的阶次， 为系统的 结构参数， 为输入函数的结构参 数，它们均为实常数 。,1.1.2 传递函数-输入/输出水平,若系统的初始条件为零，对（1）式两边取拉氏变换后稍加整理：(2)(2)式称为系统的传递函数。,1.1.3 状态空间描述-状态结构水平,系统内部模型状态空间模型。状态空 间
3、描述的一般形式为： 状态方程 : (3) 输出方程 : (4),1.2 模型结构变换,连续系统仿真要将这个系统的模型在计算机上实现出来，首先要把系统的各种描述形式转换成内部模型-状态空间模型，我们将其称为模型结构变换。,1.2.1 输入/输出水平模型到内部 模型的变换,假设一连续系统，它的数学模型如(5)式所示:（a0=1）(5) 今引进n个状态变量：， ， ，,输入/输出水平模型到内部模型的变换（续）,则有 将上述n个一阶微分方程写成矩阵形式可得（6）,输入/输出水平模型到内部模型的变换（续）,（7）外部模型变换到内部模型不唯一，所以仿 真模型也不唯一。一个系统有多种实现，最 小实现的充要条
4、件是(A、B、C)为完全能控 且完全能观测。,1.2.2 系统状态初始值变换,如果系统是非零初始条件，那么从外部模型变 换到内部模型还必须考虑如何将给定的初始条件转 变为相应的状态变量的初始值。若系统是由如下一般形式的n阶微分方程来描 述：系统初始条件为：,伴随方程法,一阶微分方程组的状态变量记为 ， 如果它们满足如下关系： （8）（9）（10）（11）该状态方程与原方程等价。,伴随方程法（续）,证明：将（8）两边分别进行微分n次，可得：（12） 其中p为微分算子符号。对（9）式两边分别进行n-j(j=1,2,n-1)次微分，可得：（13）对（10）式也引入微分算子：（14）将（12）、（13
5、）、（14）所包括的n+1个等式左 右两边分别相加，消去同类项，稍加整理后就得到原高阶微 分方程，表明了两者之间的等价关系。,伴随方程法（续）,伴随方程法明显地表示了状态变量与原输入/输出变量 及其高阶导数之间的关系，因而易于进行初始值的转换。这 样得到状态方程及输出方程：（15）其中,伴随方程法（续）,设a0=1,初值转换方程：伴随方程有多种形式，因而得到的状态方程也不 唯一。那么，实现这种初值转换的条件是什么呢？,伴随方程法（续）,考虑转换后得到的系统状态空间模型为: 即假定u的n阶导数项的系数c0=0，已知系统的初始条件为:则为了由上述初始值求出状态变量的初始值，可列出以下方 程:,伴随
6、方程法（续）,于是可得下列矩阵方程 (16)其中,伴随方程法（续）,由(16)式可得： (17) 即若 存在，则可由(17)式求出 x(t) 的初始 值。由控制理论可知，是(A、B、C)的能观判别 阵，若(A、B、C)是完全能观的，则非奇异。这 就是说，由高阶微分方程输入/输出变量初始值转变 为状态初始值的条件是：内部模型(A、B、C)是完 全能观的。,1.2.3 典型环节的传递函数,控制系统由许多元件组合而成，这些元件的物理结构和作用原理是多种多样的，但抛开具体结构和物理特点，从传递函数的数学模型来看，可以划分成几种典型环节，常用的典型环节有比例环节、惯性环节、积分环节、微分环节、振荡环节、
7、延迟环节等。,1. 比例环节,环节输出量与输入量成正比，不失真也无时间滞后的环节称为比例环节，也称无惯性环节。输入量与输出量之间的表达式为,c(t)=Kr(t),比例环节的传递函数为,式中K为常数，称为比例环节的放大系数或增益。,2. 惯性环节(非周期环节),惯性环节的动态方程是一个一阶微分方程,其传递函数为,式中 T 惯性环节的时间常数K 惯性环节的增益或放大系数,当输入为单位阶跃函数时，其单位阶跃响应为,单位阶跃响应曲线,惯性环节实例很多，如图所示的R-L网络，输入为电压u，输出为电感电流i，其传递函数,式中,3. 积分环节,输出量正比于输入量的积分的环节称为积分环节，其动态特性方程,其传
8、递函数,式中Ti为积分时间常数。,积分环节的单位阶跃响应为,它随时间直线增长，当输入突然消失，积分停止，输出维持不变，故积分环节具有记忆功能，如图所示。,上图为运算放大器构成的积分环节，输入ui(t)，输出u0(t)，其传递函数为,式中Ti = RC,4. 微分环节,理想微分环节的特征输出量正比于输入量的微分，其动态方程,其传递函数,式中 Td 称微分时间常数,它的单位阶跃响应曲线,如图所示，理想微分环节实际上难以实现，因此我们常采用带有惯性的微分环节，其传递函数,其单位阶跃响应为,曲线如下图所示，实际微分环节的阶跃响应是按指数规律下降，若K值很大而Td值很小时，实际微分环节就愈接近于理想微分
9、环节。,5. 二阶振荡环节(二阶惯性环节),二阶振荡环节的动态方程为,其传递函数,式中 为无阻尼自然振荡角频率，为阻尼比。,图中所示为RLC网络，输入为u0(t)、输出ui(t)，其动态特性方程,其传递函数,式中,6. 延迟环节(时滞环节),延迟环节是输入信号加入后，输出信号要延迟一段时间后才重现输入信号，其动态方程为,其传递函数是一个超越函数,式中称延迟时间,需要指出，在实际生产中，有很多场合是存在迟延的，比如皮带或管道输送过程、管道反应和管道混合过程，多个设备串联以及测量装置系统等。迟延过大往往会使控制效果恶化，甚至使系统失去稳定。,1.2.4 分解结构水平转换-面向结构图的模型变换,对于
10、任何一个组合的典型环节都可以用状态方程 表示：,（20）,面向结构图系统方程描述,面向结构图系统方程描述（续）,U=WY+W0y0 (21) W称为系统的连接矩阵，它描述了系统内部各环节连接情况，每个元 素Wij 表示第 j 个环节的输出到第 i 个环节的输入之间的联接系数。W0称为外部输入的连接矩阵，它描述了外部输入对系统的作用情 况。对单输入系统，W0是一个列矢量，Woj 表示外部输入信号 y0 作用在 第j个环节上的作用系数。在上图中，y0 只作用在第一个环节上，故W01 1。若为多输入系统则W0也是一个矩阵，它的列数等于输入量的个数。,系统方程转换,将(21)式代入(20)式，则可得：
11、（22）（23） 其中：Q=BDW，P=CWA ，V1=CW0 ， V2=DW0 如果Q 阵的逆存在，那么对(23)式两边左乘Q1，则得：（24） 这是一个标准的一阶常微分方程组。,系统方程转换（续）,说明：（1）矩阵方程的右端有两项与外加作用信号 有关，一 项是 ，另一项 。若外加作用函数是单位阶跃阵，此 时 ，为了便于计算，就要求V2是零向量。 如果外加作用信号是阶跃信号，那么必须限制外加作用信号 所用的那个环节Di=0。（2）只有当Q 阵能求逆时，才能获得(14)式。当系统中 各环节不存在纯微分环节或纯比例环节时就能保证Q阵可以求 逆。（3）关于Q 的逆阵不存在时的结构变换，以下例说明：
12、,结构变换例子,则 Q=BDW P=CW-A V1=CW0,因为Q 阵中出现1、3两列全零，所以Q -1不存在。其原因是1、3两个环节是比例和微分环节。对于上述系统，可以将其结构加以变换。,结构变换例子（续）,系统结构变换的程序法：设系统的状态方程为：Q 阵中有(NM)列元素为全零，这说明有(NM)个环 节的Y不出现在方程的左端。也就是说，系统中有(NM)个 代数方程。系统结构变换程序法就是通过矩阵的初等变换，把系统 中(NM)个代数方程分离出来。具体做法，就是对Q、 P、 V各矩阵线性变换，使上述方程中的各矩阵变为 ， 从而Y变为 ，方程变为：(25),结构变换例子（续）,其中： (26),
13、结构变换例子（续）,在变换的过程中，对Q、P、V作列变换，将使 Y原来的编号改变。所以程序中应将列变换的情况 记录下来，以便在最后求出 后，经过反变换而仍 能恢复Y。完成变换后，系统可以分为M个微分方程与(N- M)个代数方程。而该M个微分组的 的逆存在。 因此:(27) 由此解得 ，然后再解(N-M)个代数 方程：,结构变换例子（续）,(28)由当 时， ， 由此方程（28）可以写成：,1.3 微分方程的线性化,实际的物理系统往往有间隙、死区、饱和等各类非线性现象。严格地讲，几乎所有实际物理和化学系统都是非线性的。目前，线性系统的理论已经相当成熟，但非线性系统的理论还远不完善。因此，在工程允
14、许范围内，尽量对所研究的系统进行线性化处理，然后用线性理论进行分析不失为一种有效的方法。,当非线性因素对系统影响较小时，一般可直接将系统当作线性系统处理。另外，如果系统的变量只发生微小的偏移，则可通过切线法进行线性化，以求得其增量方程式。,非线性函数的线性化，是指将非线性函数在工作点附近展开成泰勒级数，忽略掉高阶无穷小量及余项，得到近似的线性化方程，来替代原来的非线性函数。,假如元件的输出与输入之间关系x2=f(x1)的曲线如图，元件的工作点为(x10，x20)。将非线性函数x2= f(x1)在工作点(x10，x20)附近展开成泰勒级数,当(x1x10)为微小增量时，可略去二阶以上各项，写成,
15、其中 为工作点(x10，x20)处的斜率，即此时以工作点处的切线代替曲线，得到变量在工作点的增量方程，经上述处理后，输出与输入之间就成为线性关系。,图2-8为一铁芯线圈，输出为ui(t)，输入为i(t)。 线圈的微分方程为,当工作过程中线圈的电压和电流只在工作点（u0,i0）附近变化时，即有,线圈中的磁通 对 也有增量变化，假如在i0附近连续可微，将在i0 附近展开成泰勒级数，即,因是微小增量，将高阶无穷小量略去，得近似式,这就是铁芯线圈的增量化方程，为简便起见，常略去增量符号而写成,1.4 方框图,在控制工程中，为了便于对系统进行分析和设计，常将各元件在系统中的功能及各部分之间的联系用图形来
16、表示，即方框图和信号流图。,1.4.1 方框图,方框图也称方块图或结构图，具有形象和直观的特点。系统方框图是系统中各元件功能和信号流向的图解，它清楚地表明了系统中各个环节间的相互关系。构成方框图的基本符号有四种，即信号线、比较点、传递环节的方框和引出点。,1.5 在MATLAB中数学模型的表示,线性系统理论中常用的数学模型有微分方程、传递函数、状态空间表达式等,而这些模型之间又有着某些内在的等效关系。MATLAB主要使用传递函数和状态空间表达式来描述线性时不变系统(Linear Time Invariant简记为LTI)。,1.5.1 传递函数,单输入单输出线性连续系统的传递函数为,其中mn。
17、G(s)的分子多项式的根称为系统的零点,分母多项式的根称为系统的极点。令分母多项式等于零,得系统的特征方程为,D(s)=a0sn+a1sn1+an1s+an=0,系统的传递函数在MATLAB下可由其分子和分母多项式唯一地确定出来,其格式为,sys=tf(num,den),其中num为分子多项式，den为分母多项式,num=b0,b1,b2,bm； den=a0,a1,a2,an；,对于其它复杂的表达式,如,可由下列语句来输入,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),Trans
18、fer function:,1.5.2 传递函数的特征根及零极点图,传递函数G(s)输入之后,分别对分子和分母多项式作因式分解,则可求出系统的零极点,MATLAB提供了多项式求根函数roots(),其调用格式为,roots(p),其中p为多项式。,例如,多项式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.00
19、00 3.0000 0.0000 4.0000,而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,p,z=pzmap(num,den) 其中, p传递函数G(s)= numden的极点z传递函数G(s)= numden的零点 例如,传递函数,传递函数在复平面上的零极点图,采用pzmap()函数来完成,零极点图上,零点用“。”表示,极点用“”表示。其调用格式为：,用
20、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, 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
21、(s)表达式,pzmap(num,den) %零极点图 title(pole-zero Map),零极点图如图所示 ：,1.5.3 控制系统的方框图模型,若已知控制系统的方框图,使用MATLAB函数可实现方框图转换。,1.串联 如图所示G1(s)和G2(s)相串联,在MATLAB中可用串联函数series( )来求G1(s)G2(s),其调用格式为 num,den=series(num1,den1,num2,den2) 其中：,2.并联如图所示G1(s)和G2(s)相并联,可由MATLAB的并联函数parallel( )来实现,其调用格式为,num,den=parallel(num1,den1
22、,num2,den2),其中：,3.反馈,反馈连接如图所示。使用MATLAB中的feedback( )函数来实现反馈连接,其调用格式为,num,den=feedback(numg,deng,numh,denh,sign),式中：,sign为反馈极性,若为正反馈其为1,若为负反馈其为1或缺省。,例如 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=,MATLAB中的函数series,parallel
23、和feedback可用来简化多回路方框图。另外,对于单位反馈系统,MATLAB可调用cloop( )函数求闭环传递函数,其调用格式为,num,den=cloop(num1,den1,sign),1.5.4 控制系统的零极点模型,传递函数可以是时间常数形式,也可以是零极点形式,零极点形式是分别对原系统传递函数的分子和分母进行因式分解得到的。MATLAB控制系统工具箱提供了零极点模型与时间常数模型之间的转换函数,其调用格式分别为,z,p,k= tf2zp(num,den) num,den= zp2tf(z,p,k),其中第一个函数可将传递函数模型转换成零极点表示形式,而第二个函数可将零极点表示方式
24、转换成传递函数模型。,例如 G(s)=,用MATLAB语句表示,num=12 24 12 20;den=2 4 6 2 2; z,p,k=tf2zp(num,den)z= 1.9294 0.03530.9287i0.03530.9287i,p=0.95671.2272i 0.95671.2272i 0.04330.6412i 0.04330.6412ik=6,即变换后的零极点模型为 G(s)=,可以验证MATLAB的转换函数,调用zp2tf()函数将得到原传递函数模型。 num,den=zp2tf(z,p,k)num = 0 6.0000 12.0000 6.0000 10.0000den =
25、 1.0000 2.0000 3.0000 1.0000 1.0000,即,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为输出量个数。,在一般情况下,控制系统的状态空间表达式项简记为(A,B,C,D)。,例如：设一个双输入双输出系统的状态空间表达式为,系统模型可由MATLAB命令直观地表示,A=1,
26、2,4;3,2,6;0,1,5 B=4,6;2,2;0,2 C=0,0,1;0,2,0 D= zeros(2,2),MATLAB的控制系统工具箱提供了由状态空间表达式转换成传递函数或由传递函数转换成状态空间表达式的转换函数ss2tf( )和tf2ss( )。其调用格式为,num,den=ss2tf(A,B,C,D,iu),反过来,若已知系统的传递函数,求取系统状态空间表达式的调用格式为 A,B,C,D=tf2ss(num,den),例如系统的传递函数为,系统的状态空间表达式为,num=1,2,3; den=1,3,6,1; A,B,C,D=tf2ss(num,den),A= -3 -6 -11
27、 0 0 0 1 0 B = 100 C = 1 2 3 D = 0,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的状态响应。,例 试求下列系统的单位脉冲响应MATLAB命令为： t=0：0.1：40； num=1； den=1，0.3，1； impulse(nu
28、m，den，t)； grid； title(Unit-impulse Response of G(s)=1/(s2+0.3s+1) 其响应结果如图所示。,例 系统传递函数为求取其单位脉冲响应 的MATLAB命令为 t=0：0.1：10；num=1； den=1，1，1； y，x，t=impulse(num，den，t) plot(t，y)；grid xlabel(t)； ylable(y)； 其响应结果如图所示。,1.6.2 单位阶跃响应,当输入为单位阶跃信号时，系统的输出为单位阶跃响应，在MATLAB中可用step( )函数实现，其调用格式为 y， x， t=step(num， den， t
29、) 或step(num， den),例 求系统传递函数为num=1； den=1，0.5，1； t=0：0.1：10； y，x，t=step(num，den，t)； plot(t，y)；grid； xlabel(Time sec t)； ylabel(y) 响应曲线如图所示,图 单位阶跃响应,1.6.3 斜坡响应,在MATLABA中没有斜坡响应命令，因此，需要 利用阶跃响应命令来求斜坡响应。根据单位斜坡响应 输入是单位阶跃输入的积分。当求传递函数为斜坡响 应时，可先用除得，再利用阶跃响应命令即可求得斜 坡响应。,例 已知闭环系统传递函数 对单位斜坡输入 则,系统单位斜坡响应的MATLAB命令：
30、 num=1； den=1，0.3，1，0； t=0：0.1：10； c=step(num，den，t)； plot(t，c)； grid； xlabel(t sec)； ylabel(Input and Output) 其响应结果如图所示。,1.6.4 任意函数作用下系统的响应,用线性仿真函数lsim来实现，其调用格式为 y，x=lsim(num，den，u，t) 式中 ；y(t)为系统输出响应；x(t)为系统状态响应；u为系统输入信号；t为 仿真时间。,例 反馈系统如图 (a)所示，系统输入信号为图 (b)所示的三角波，求取系统输出响应。,MATLAB实现指令 numg=10，20；deng=1，10，0； num，den=cloop(numg，deng，-1)； v1=0：0.1：2； v2=1.9：-0.1：-2； v3=-1.9：0.1：0； t=0：0.1：8； u=v1，v2，v3； y，x=lsim(num，den，u，t)； plot(t，y，t，u)； xlabel(Time sec)；,1.6.5 Simulink中的时域响应举例,例图的Simulink的仿真框图可演示系统对典型信号的时间响应曲线，图中给出阶跃响应曲线。,返回,
