《控制系统数字仿真》PPT课件.ppt

返回总目录,第4章 控制系统数字仿真,数字仿真是在数字机上建立系统模型并利用模型做实验，所以，进行数字仿真首先要建立描述被仿真系统的数学模型，并将此模型转换成计算机可接受的、与原模型等价的仿真模型，然后编制程序，使模型在计算机上运行。本章主要讲述数字仿真的基本理论与方法。,连续系统数值积分方法,连续系统的数学模型，一般都能以微分方程的形式给出，所以，连续系统数学仿真算法问题通常可归结为如何用计算机来求解微分方程的问题。数值积分法是解决该问题的重要方法之一，例如已知,(4.1),求,连续系统数值积分方法,解：对式(4.1)两边积分，则,当,时,令,则,(4.2),数值积分法是解决在已知初值的情况下，对,进行近似积分，对,进行数值求解的方法，在数学上称为微分方程初值问题的数值方法。,连续系统数值积分方法,一.欧拉法,表示函数,在,和,相邻两次采样时刻之间的积分。,若将此定积分中的,近似看成常数，即,则,由式(4.2)得,即,(4.3),连续系统数值积分方法,式中，,为计算步距。,，1，2，,式为(4.3)著名的欧拉积分公式。欧拉公式计算简单，但精度较低，其原因是将,看成常数,，从而用矩形面积代替准确的曲面面积,，形成了较大的误差面积。为了提高计算精度产生了梯形法。,连续系统数值积分方法,二. 梯形法,梯形法是用梯形的面积近似代替定积分 ，即 式中， ； ，所 以 即,(4.4),连续系统数值积分方法,从式(4.4)可以看出，用梯形法进行数值积分，会出现一个问题：在计算,时，先要用,去计算式(4.4)右端的,，而此时 还未求出，,显然，这是不可能实现的。所以，一般采用欧拉公式先预报一个,然后将预报 的代入式(4.4)进行校正，求出 ，即,(4.5),式中，,一般称式(4.5)为预报校正公式。显然，梯形法比欧拉法精度要高，同时，每前进一个步距，计算工作量也比欧拉法约多一倍。,连续系统数值积分方法,三. 龙格-库塔法,对于微分方程式(4.1)，若在其初值附近展开成泰勒级数，并只取前三项，则有,(4.6),设式(4.1)的解可以写成如下形式,(4.7),连续系统数值积分方法,将用 二元函数泰勒级数展开式展开，并只取前三项，则,将 、 代入式(4.7),(4.8),比较式(4.6)与式(4.8)，得,连续系统数值积分方法,三个方程中有四个未知数，因而解不是唯一的。所以，限定 ，则求出,。代入式(4.7)得,写成一般形式,(4.9),连续系统数值积分方法,若在 的假设下， 时的准确解为 ，,为用式(4.9)求得的近似解，则它们之差 ，称为此 时计算的截断误差。,式(4.6)只取到泰勒级数展开式中 的二阶导数项，略去了三阶以上高阶导数项。为纪念提出该方法的德国数学家C. Runge和M. W .Kutta，称这种计算方法为二阶龙格-库塔法。其截断误差正比于步长 的三次方，同理若在式(4.6)的计算中，取到的 三阶或四阶导数项，则有相应的三阶或四阶龙格-库塔法，相应的截断误差也应正比于 或 。,连续系统数值积分方法,一般在计算精度要求较高的情况下，多使用四阶龙格-库塔法。其计算公式为,(4.10),连续系统数值积分方法,将欧拉公式与 的泰勒级数展开式比较可知，欧拉公式即泰勒级数展开式只取前两项。而将预报校正公式与二阶龙格库塔公式比较，发现它们完全是一样的。这样，通过对龙格库塔法的介绍，可将前面讲的几种数值积分方法统一起来，都可看成是 在初值附近展开成泰勒级数所产生的，只是其泰勒级数所取项数的多少不同。欧拉法只取前两项，梯形法与二阶龙格库塔法取了前三项，四阶龙格库塔法取了前五项。从理论上说，取的项数越多，计算精度越高，但相应的计算公式更复杂，计算工作量也更大。,连续系统数值积分方法,以上介绍的几种数值积分公式，有一个共同的特点，在计算 时只用到 ，而不直接用 ， ， 各项的值，即在本次计算中，仅仅用到前一次的计算结果，而不需要利用更前面各步的结果。这类计算方法称为单步法。单步法运算有下列优点：,(1) 需要存储的数据量少，占用的存储空间少。 (2) 只需要知道初值，就可启动递推公式进行运算，具有这种能力的计算方法称为自启动的计算方法。 (3) 容易实现变步长运算。 与单步法相对应的还有一类数值积分方法，在它的数值积分公式中，本次计算不仅利用前一次的计算结果，还必须利用更前面各次结果，此类方法称为多步法。如四阶阿达姆斯(Adams)积分公式,连续系统数值积分方法,即是多步法计算公式。 多步法与单步法相比，欲达到相同精度，计算工作量较小，从式(4.11)和四阶龙格库塔式(4.10)的比较中可明显看出。在计算 时，式(4.11)只需计算 ，而 ， ， 已由前三次计算求出，而四阶龙格库塔公式每前进一步都要计算 ， ， ， ，相当于四次计算右端函数 ，因此在相同的条件下多步法比单步法要快。,(4.11),连续系统数值积分方法,四. 微分方程数值积分的矩阵分析方法,前面介绍的数值积分公式，都是以求解单个典型微分方程 进行介绍的。而工程实际系统中大量的仿真对象是以一阶微分方程组或矩阵微分方程的形式给出，如,或,连续系统数值积分方法,式中，,在这种情况下，各数值积分公式显然应采用相应的矩阵形式,连续系统数值积分方法,(1) 欧拉公式,(4.12),(2) 梯形公式,(4.13),(3) 二阶龙格-库塔公式,(4.14),连续系统数值积分方法,(4) 四阶龙格-库塔公式,(4.15),连续系统数值积分方法,连续系统数值积分方法,对一个 维向量 ，每前进一个步距 ，至少要求 个 之值。,对常见的线性定常系统 ，四阶龙格库塔法的四个可表示为,连续系统数值积分方法,连续系统数值积分方法,对 维向量 ，取 是有四个分量的一维矢量，其中 再取 为一零向量，其中 ， 。,则四个求,的公式可合并为一个公式,(4.16),在以后介绍的程序中，四阶龙格库塔法的计算将建立在上式的基础上。,连续系统数值积分方法,五. 数值积分方法的计算稳定性,这里所说的数值积分方法的计算稳定性，在数学上是指微分方程初值问题算法的数值稳定性或计算稳定性。由于在系统仿真运算阶段大量求解微分方程问题时，实际上是用不定数值积分方法求解微分方程的问题。如欧拉法就是在已知 时，用 完成对微分方程 的求解。习惯上将这些方法称为数值方法。,从稳定性理论中，可以知道如何从系统的微分方程或传递函数去判断该系统的稳定性。那么，对于一个稳定的连续系统，当用某数值积分方法进行仿真计算时，是否仍然稳定呢？先看下面的例子。,连续系统数值积分方法,【例4.1】 已知微分方程,其精确解为,取 ，用欧拉法和四阶龙格-库塔法计算 时的 ：,连续系统数值积分方法,欧拉法,四阶龙格-库塔法,精确解,连续系统数值积分方法,显然，此时数值积分法计算的结果是错误的。为什么会出现这种情况呢？这是因为数值积分方法只是一种近似方法，它在反复的递推运算中将引入误差。若误差的积累越来越大，将使计算出现不稳定，从而得出错误的结果。所以，原系统的稳定性与数值积分法计算的稳定性是不同的两个概念。前者用原系统的微分方程、传递函数来讨论，后者用逼近微分方程的差分方程来讨论。由于选用的数值积分法不同，即使对同一系统，差分方程也各不相同，稳定性也各不一样。如何分析数值积分法的数值计算稳定性呢？一般来讲，对高阶微分方程的数值计算稳定性作全面分析是比较困难的，通常用简单的一阶微分方程来考查其相应差分方程的计算稳定性。,连续系统数值积分方法,微分方程,(4.17),称测试方程。据稳定性理论，当其特征方程的根在平面的左半平面，即根的实部 时，则原系统稳定。此时相应的数值积分法的计算稳定性如何呢？,连续系统数值积分方法,1. 欧拉法的计算稳定性,对式(4.17)按欧拉公式计算有,(4.18),用式(4.18)计算的 ，并不是 时的真值，而只是真值包含了各种误差的近似值。随着递推次数的增加，此误差是否会不断扩大，使 完全不能表示此时的真值呢？这决定了此差分方程计算的稳定性。为简化问题的讨论，假设只在 时刻有误差引入，这样 的误差 仅由 的误差引起，所以,(4.19),连续系统数值积分方法,式(4.19)减式(4.18)得误差方程,同理,当 ，则 ，表明此算法若在计算中的某步引入了误差，随着计算步数的增加，这个误差将逐渐扩大，以致差分方程的解完全失真。,连续系统数值积分方法,反之，当 ，则 随着计算步数的增加，误差 逐渐减小并趋于零或有界，在这种情况下，称此差分方程的算法是计算稳定的。显然，合理的选择步长使其满足 ，是保证欧拉法计算稳定的重要条件。,连续系统数值积分方法,2. 梯形法的计算稳定性,对测试方程按梯形法计算公式，则,(4.20),所以,(4.21),连续系统数值积分方法,若原系统稳定，则根据稳定性理论， 应是小于零的实数。对式(4.21)， 取任何正数都可以保证下式成立：,所以,这就表明，某步计算引入的误差，将随着计算步数的增加而减小。也就是说，梯形法的计算在任何步长下都是稳定的，是一种绝对稳定的计算方法。,连续系统数值积分方法,六. 数值方法的计算精度、速度、稳定性与步距的关系,从上面的分析可以看出，一般数值积分方法的稳定性与所选用的计算步距有关。各种数值积分方法的差分方程是对原微分方程的近似逼近，存在明显的截断误差。且计算机的字长有限，计算只能限制在有限位，将引入舍入误差。这些误差都和计算步距密切相关，所以计算步距是影响计算精度、速度和稳定性的重要因素。若步距取得较大，截断误差就会相应地增大；反之，若步距取得较小，截断误差就会减小，但在给定时间范围内，计算次数必然增加，使舍入误差积累，相应地增加。,连续系统数值积分方法,这两种误差与步距的关系可用图4.1表示。从图中可知，两种误差对步距的要求是矛盾的，但两者之和有一个最小值，步距最好能选在最小值。然而，实际要做到这一点是很困难的。一般只能根据经验确定一个 附近的合理步长区，如可将 限制在系统的最小时间常数数量级上。,图4.1 误差与步距关系,控制系统的结构及其描述,一. 控制系统的典型结构形式,1. 单输入-单输出开环控制结构 控制系统结构如图4.2所示。其中控制器、控制对象等环节可用第3章所述任何一种数学模型描述, 为参考输入量, 为控制量, 为输出量。,图4.2 单输入-单输出开环控制结构图,控制系统的结构及其描述,2. 单输入-单输出前馈控制结构 控制系统结构如图4.3所示。这种结构形式本身仍为开环控制形式，但为补偿输入引起的误差，在已知误差变化规律的情况下，加入补偿环节，对误差作提前修正。补偿环节同样可用任一种数学模型描述。,图4.3 单输入-单输出前馈控制结构图,控制系统的结构及其描述,3. 单输入单输出闭环控制结构 这是控制系统中应用最广泛的控制结构形式，大多数控制系统为保证其控制精度、性能指标要求，都采用这种闭环负反馈形式，如图4.4所示。 该系统控制原理主要是根据偏差(或误差)确定控制量，使输出量按希望的精度变化，减小或消除偏差，即所谓的偏差控制系统。图中 为偏差量(或误差量)， 为反馈量。 反馈环节可用任一种数学模型描述。当反馈为单位反馈时，有,则偏差,控制系统的结构及其描述,图4.4 单输入-单输出闭环控制结构图,称作误差，而此时的闭环系统称作单位反馈系统。,控制系统的结构及其描述,凡是单输入-单输出控制系统结构均能方便地表示为图论中的拓扑结构形式，即使系统环节再增多，反馈和前馈联系再复杂，表述也很清晰，如图4.5所示。图中各环节用序号表示，每个环节都有自己的输入和输出变量 和 (1,2,.，6)，对应其相互关系的数学描述，可以是第3章所述形式的任一种。,图4.5 控制系统拓扑结构图,控制系统的结构及其描述,4. 多输入-多输出控制系统结构 若干个单输入-单输出控制系统结构形式通过一定方式组合在一起，构成了多输入-多输出控制系统结构形式，用来描述较复杂的多变量控制系统。 根据不同的环节组合形式，也有不同的控制结构。如变量和环节有单向耦合作用关系的结构形式见图4.6(a)、交叉反馈耦合作用关系的结构形式见图4.6(b)等。,控制系统的结构及其描述,(a) 单向耦合结构 (b) 交叉反馈耦合结构 图4.6 多输入-多输出控制结构图,控制系统的结构及其描述,二.控制系统的典型环节描述,任何一个复杂的控制系统，都是由一些简单的不同类型的典型环节采用串联、并联、反馈等形式组合而成的。若对常见的一些典型环节能准确地加以定量描述，则复杂系统的描述也只是复杂在各部分的相互连接关系上。采用前述系统关系描述，把这种连接关系用相应的数学关系表达出来，就可以得到在计算机上能方便运行的“仿真模型”。经典控制理论中常见的典型环节如下所述：,控制系统的结构及其描述,(1) 比例环节,(2) 惯性环节,控制系统的结构及其描述,(3) 惯性比例环节,(4) 积分环节,控制系统的结构及其描述,(5) 比例积分环节,(6) 二阶振荡环节,式中，为阻尼比，01时，系统为欠阻尼二阶振荡情况； 为无阻尼自然振荡频率。,控制系统的结构及其描述,(7) 纯滞后环节,式中，,(8) 高阶线性环节,为纯滞后时间。,控制系统的结构及其描述,式中， , 为分母多项式系数，亦称特征多项式系数； , ,为分子多项式系数；它们均为反映高阶线性环节动态响应性能的重要参数。,以上典型环节，无论其阶次高低，均可应用MATLAB语言方便地求出输入量作用下输出量的响应性能。 和 之间的数学表达关系可以由传递函数 方便地转化为其他各种形式，如状态方程、零极点增益、部分分式等。,控制系统的结构及其描述,三. 控制系统的连接矩阵,对图4.7所表示的线性系统，各环节均为线性的，在各自的输入量 作用下，给出各自的输出量 ，这种作用关系是通过各环节的数学描述体现出来的，其数学关系可以是(4.2.2)节所述的任何一种。但各环节之间存在相互作用， 不是孤立的，只要与其他环节有连接关系，就要受到相应 变化的影响。因此，要完整地将系统描述出来，还应该分析各环节输出 对其他环节有无输入作用，才能完整地进行仿真分析。,控制系统的结构及其描述,图4.7 控制系统拓扑结构,控制系统的结构及其描述,根据图4.7中 、 的连接关系，可逐个写出每个环节输入 受哪些环节输出 的影响，现列写如下：,(4.22),控制系统的结构及其描述,引入向量,则式(4.22)可表示为,控制系统的结构及其描述,即,式中，,控制系统的结构及其描述,式中，W为连接矩阵，nn型，矩阵中元素清楚地表示出各环节之间的连接关系，W0为输入连接矩阵，nm型(当参考输入为m维时)，阵中元素表示环节与参考输入之间的连接关系。仔细研究连接矩阵W，可从其元素值直接看出各环节之间的连接情况，即,wij0，环节j不与环节i相连； wij 0，环节j与环节i有连接关系； wij 0，环节j与环节i直接相连(wij=1)或通过比例系数相连(wij为任意正实数)； wij 0，环节j与环节i直接负反馈相连(wij=-1或通过比例系数负反馈相连(wij为任意负实数)； wii 0，环节i单位自反馈(wii=l或wii=-1)或通过比例系数自反馈(wii为任意实数)。,面向结构图的数字仿真,一. 典型闭环系统的数字仿真,1. 典型闭环系统结构形式 控制系统最常见的典型闭环系统结构如图4.8所示。其中,面向结构图的数字仿真,图4.8 典型闭环系统结构图,面向结构图的数字仿真,2. 系统仿真模型与求解思路 所谓仿真模型是指经一定方式把数学模型转化为便于在计算机上运行的表达形式。这种表达形式往往是一些适合于具体编程实现的数学关系式。 对图4.8所示系统的开环传递函数，可按照能控标准型写出其开环状态方程,(4.23),式中,面向结构图的数字仿真,其中， ， 矩阵中， ， ，( 1,2,， 0,1,，, )为式(3.4)首 一化后分母、分子各系数，即 ， , 且 。,由图又知，控制量 ，代入式(4.23)得,再由，,，则,(4.24),即得系统闭环状态方程，其中,为系统闭环系数矩阵，而输入矩阵B和输出矩阵C不变。这就是图4.8系统的仿真模,面向结构图的数字仿真,仿真模型一旦确立，就可以考虑求解与编程实现。观察式(4.24)可知，该式其实为一个一阶微分方程组的矩阵表达形式，而数值积分法最适宜解一阶微分方程，当采用四阶龙格库塔法求解此闭环状态方程时，其步骤如下：,对 知,此为对应 个状态变量的一阶导数 的n维向量表达式。其中，r为随时间变化的已知输入函数。于是，当求 解过程进行到 时刻，欲求 时刻各量，需先求,面向结构图的数字仿真,以上各系数 、 、 、 可认为是对应维状态变量在时刻的四组斜率，每组为维斜率向量，即每组个龙格库塔系数。,面向结构图的数字仿真,最后，再由,求得 时刻状态 ，立即可得输出相应时刻值,按以上算式，取k0，1，2，N不断递推，即求得所需时刻各点的状态变量 和输出量 。,面向结构图的数字仿真,3仿真程序框图与实现,基于以上仿真模型和求解思路，我们可以着手考虑面向这类典型闭环控制系统结构图的仿真程序的编制与实现。 作为系统仿真程序，使用时应尽可能方便，使用者只要将开环传递函数 的分母、分子各系数 ， 和反馈系数输入计算机，计算机就掌握了有关该系统的基本信息模型参数。然后求取首一化表达形式，以及形成开、闭环状态方程各矩阵等步骤均由仿真程序自动完成，无需人工干预。因此，程序应有输入数据模块和初始化程序模块。,面向结构图的数字仿真,在参考输入函数 的作用下，系统输出 开始随时间变化，仿真程序应能按照给定的计算步长，采用已确定的数值算法，对系统中各状态变量和输出逐点变化情况进行求解运算。这部分模块是整个仿真程序的核心部分，计算速度、精度误差均取决于它。通常称为运行程序模块。 在使用者规定的时间范围内，将计算数据按照一定要求存储，并在仿真结束时，按使用者指定格式输出仿真结果，以便对系统进行分析研究，这就是所谓的输出程序模块。 综上所述，构成一个完整的仿真程序必须至少建立：输入数据块、初始化块、运行计算块、输出结果块四个子程序模块，才能正常地完成仿真任务。,面向结构图的数字仿真,1) 程序框图 程序框图如图4.9所示。,图4.9 典型闭环系统的仿真程序框图,面向结构图的数字仿真,2) 程序语句 用MATLAB语言编制的程序语句如下： EXP4-1.M % 输入数据 aao,a1，an； %n1维分母系数向量 bb0,b1,.,bm； %ml维分子系数向量 XOx10,x20,.,xn0； %状态向量初值 VVo ； %反馈系数 nno ； %系统阶次 T0=t0； %起始时间 Tf=tf； %终止时间 h=h0； %计算步长 R=r %阶跃输入函数的幅值 % 形成开、闭环系数矩阵,面向结构图的数字仿真,B=b/a(1)；a=a/a(1)；A=a(2:n+1) %首一化处理 A=rot90(rot90(eye(n-1,n)；-fliplr(A)；%形成能控标准形A阵 B=zeros(1,n-1)，1； %形成输入矩阵B(n维列向量) m1=length(b)； %分子系数向量维数m+1 C=fliplr(b),zeros(1,n-m1)； %形成输出矩阵C(n维行向量) =A-B*C*V； %形成闭换系数矩阵X=X0；y=0；t=T0； %设初值，准备开始递推运算 % 运算求解 N=round(Tf-T0)/h)； %确定输出点数 For i=1:N %四阶龙格-库塔法 K1=*X+B*R K2=*(X+h*K1/2)+B*R； K3=*(X+h*K2/2)+B*R； K4=*(X+h*K3)+B*R %求各系数K X=X+h*(K1+2*K2+2*K3+K4)/6 %求状态 y=y,C*X； %求输出并以向量形式保存 t=t,t(i)+h； %输出对应时刻以列向量形式保存 end % 输出结果 t,y %输出数据形式结果 plot(t,y) %输出曲线形式结果,面向结构图的数字仿真,以上即为采用MATLAB语言实现面向典型闭环系统数字仿真的程序。 上面程序中用到了几个MATLAB语言的特殊功能，使得程序较为简练。简要说明如下： eye(m,n) 产生mn型 阵，即主对角元素均为1，其余元素为零的长方阵 length() 求取括号内向量的维数 zeros(m,n) 产生mn型零阵 rot90() 将括号内矩阵左转90 fliplr() 将括号内矩阵或向量左右翻转 round() 对括号内运算取整 plot(t,y) 将打印为横轴，为纵轴的曲线,面向结构图的数字仿真,【例4.2】 求图4.10所示系统的阶跃响应 的数值解。,图4.10 控制系统结构图,面向结构图的数字仿真,解：该系统结构形式为典型闭环控制系统，用EXP4-1.M求解过程如下： (1) 取开环放大系数k=1，反馈系数v=1，阶跃输入幅值r=1。 (2) 利用MATLAB语言中conv()卷积函数功能，先将系统开环传递函数 化为要求的 传递函数形式的分母、分子多项式，系数向量为 和 。,(3) 设系统状态向量初值 均为零。 (4) 系统运行参数 ， ， ， 。,面向结构图的数字仿真,(5) 按以上步骤和参数，在MATLAB语言环境下，输入以下命令语句： k=1； a=conv(1 0 0,conv(0.25 1,0.25 1)； a= 0.0625 0.5000 1.0000 0.0000 0.0000 b=2*k k b= 2 1 X0=0 0 0 0 v=1；n=4；T0=0；Tf=10；h=0.25；R=1； EXP4-1.M % 运行EXP4-1.M仿真程序 Plot(t,y) % 输出曲线形式结果(见图4.11),面向结构图的数字仿真,图4.11 输出响应 曲线,面向结构图的数字仿真,二. 复杂连接的闭环系统数字仿真,实际工程中常常遇到的是复杂结构形式的控制系统，这些控制系统是由若干典型环节按照一定规律连接而成的，可用前述的形式将其描述为复杂连接的闭环系统结构图。要采用上节所述方法仿真，必须先将复杂结构图简化成为图4.8的典型结构式，求出开环传递函数，然后运用上节程序进行仿真分析。但存在如下问题：系统结构复杂，内部存在交叉耦合、局部回环等情况时，化简并非易事，尤其通过手工化简使工作量陡增；有时在分析中，还需要得知结构图中某些环节的输出变量情况，若经过化简消去这些环节，则不便进行观察分析；在分析中常常需要改变某个参数，观察其对输出的影响，但改动一个参数值，就有可能需将所有开环传递函数分子、分母系数统统重新计算再次输入计算机，很不方便；对实际系统中存在的非线性情况，更是无法加以考虑。鉴于以上原因，编制面向复杂连接的闭环系统结构图的系统仿真程序，必须使其能克服上节仿真程序的不足之处且应该具有以下特点： (1) 可按照系统结构图输入各环节参数，对应关系明确，改变参数方便。 (2) 可方便地观察各环节输出动态响应。 (3) 各环节存在非线性特性时易于处理。,面向结构图的数字仿真,1典型环节的二次模型化 复杂连接闭环系统数字仿真程序的基本思路是：与实际系统的结构图相对应，在计算机程序中也应能方便地给出表示各实际环节的典型环节，并将环节之间的连接关系输入计算机，由计算机程序自动形成闭环状态方程，运用数值积分方法求解响应。因此，选定典型环节是很重要的，要使其既具有代表性，又不致于造成输入数据复杂烦琐。考虑控制系统典型环节有以下几种情况。,比例环节：,积分环节：,面向结构图的数字仿真,比例积分：,惯性环节：,一阶超前滞后环节：,二阶振荡环节：,面向结构图的数字仿真,可见，除二阶振荡环节外，都是一阶环节，完全可用一个通用一阶环节(见图4.12)表示，即,图4.12 典型一阶环节,(4.25),式中， 为第i个环节的输出； 为第i个环节的输入；n为系统中的环节数(也就是系统的阶次)。而二阶振荡环节可以化为如图4.13方式连接而成的等效结构图。,面向结构图的数字仿真,由图可见，二阶振荡环节完全可以用一阶环节等效连接得到。因此，选定式(4.25)作为组成复杂系统仿真模型的典型环节是适合的。,图4.13 二阶振荡环节的等效结构图,面向结构图的数字仿真,设：输入向量 ，其中各分量表示各环节输入量；,输出向量 ，其中各分量表示各环节输出量；,模型参数矩阵为,于是系统中所有环节输出、输入关系用矩阵表示如下:,面向结构图的数字仿真,2系统的连接矩阵与仿真求解 1) 连接矩阵 仅把系统中各环节描述出来还不够，要进行数值积分求解，还必须把各环节之间的相互作用关系清楚地表达出来，这种表达是通过建立本章第4.2节所述的连接矩阵实现的。为说明问题，现举例说明如下：,图4.14 某闭环系统结构图,面向结构图的数字仿真,设系统结构如图4.14所示，图中编号1、2、3表示三个典型环节， 、 为环节间直接作用的比例系数。各环节输入 与输出 的关系为,整理为矩阵形式得,(4.27),面向结构图的数字仿真,式中,(4.28),面向结构图的数字仿真,2) 求系统的解 将式(4.27)代入式(4.26)，则,整理，得,为避免输入函数为阶跃函数时,初始点导数 , 为阶跃函数)为无穷大，要求使 ，则,即,面向结构图的数字仿真,式中，,若Q矩阵的逆矩阵 存在，则式(4.29)两边同时左乘 ，得,即,式中，,(4.30),面向结构图的数字仿真,式(4.30)是典型的一阶微分方程组矩阵形式，利用4.3.1节介绍的求解方法，可方便地求出各环节的输出响应。然而，建立系统仿真模型还应注意以下两点： (1) 存在的条件。 为了分析 存在的条件，分析图4.15所示例子。,图4.15 系统结构图,面向结构图的数字仿真,面向结构图的数字仿真,由于Q阵中有两列全为零，所以 不存在。其原因在于它的第一个环节为纯比例环节，第三个环节为纯微分环节，按典型环节给出的方程分别为,显然不是微分方程而是代数方程，因而使 不存在。为避免系统中出现代数方程，必须保证系统中所有环节的 ，从而使环节能用微分方程给出。从数学上说，因为 ， 保证了B非奇异，即使此时 ，也可保证Q非奇异。所以当系统中出现纯比例、纯微分这类使 的环节时，则应设法与其他环节合并处理，或设法化为系统可接受的环节。,面向结构图的数字仿真,(2),的条件,当被仿真系统的输入函数有阶跃函数时，阶跃函数在t=0点的 ,此时若 ,将使微分方程中的 为无穷大。使,无解。为避免出现这种情况，要求 。由于D为NN对角阵， 为NR矩阵，而 中只有承受输入作用的那些环节的 ，其余都为零。,面向结构图的数字仿真,要使 ，只要被输入函数作用的那些环节的 即可。,面向结构图的数字仿真,3仿真程序框图与实现 根据前面分析，在已建立的仿真模型基础上，用MATLAB语言编程实现系统仿真。具体实现时应注意采取合理的方法。 l) 系统参数输入方法 为利用MATLAB语言的矩阵处理功能，系统参数输入按各环节 输入参数矩阵P。,如,然后由程序自动形成式(4.26)中A、B、C、D各矩阵和式(4.30)中闭环状态方程各系数矩阵。,面向结构图的数字仿真,2) 连接矩阵输入方法 连接矩阵W和 中大量的元素是零元素，而表示连接关系的 非零元素并不多。因此，程序中采用只输入非零元素的方式，可大大加快输入速度，输入后由程序内部自动形成完整的连接矩阵W和 ，具体方法是建立非零元素矩阵 (m3型)，将非零元素按照i，j， 次序逐行输入。其中，i为被作用环节号；j为作用环节号； 为作用关系值(连接系数)；m为非零元素个数(包括 阵中非零元素)。,例如，按照图4.14和式(4.28)，相应 的应按如下格式写出：,面向结构图的数字仿真,3) 程序框图 程序框图如图4.16所示。框图中四阶龙格-库塔法求解部分完全可调用4.2 节程序EXP4-1.M中的运行程序模块。,面向结构图的数字仿真,面向结构图的数字仿真,4) 程序实现 采用MATLAB语言编程实现。 % 输入数据 Pa1，b1，cl，dl； a2，b2，c2，d2； ； an，bn，cn，dn； %各环节参数输入,； %连接阵非零元素输入 nno； %环节个数系统阶次 Y0y0 %阶跃输入幅值 Yt0y10，y20，yn0 %各环节初值 hh0 %计算步长 Lll1 %打印间隔点数(每隔点l1输出一次) TO=t0 Tftf noutnput %输出环节编号 % 形成闭环各系数阵,面向结构图的数字仿真,Adiag(P(:,1)；Bdiag(P(:,2)； Cdiag(P(:,3)；Ddiag(P(:,4)； mlength(Wij(:,1)； %求非零元素个数 W0zeros(n,1)；Wzeros(n,n)； %建立初始W(nn型方阵)、W0阵(n维列向量) For k=l : m if (Wij(k,2)O): WO(Wij(k,1)Wij(k,3)； else W(Wij(k,l)，(Wij(k,2)Wij(k,3)； end； end； %求Wo阵和W阵 ； ； %求和阵 ； ； %求和阵 ； ； %形成闭环系数阵,面向结构图的数字仿真,% 数值积分求解 YYtO；yY(nout)；t=T0； %置初值，做好求解准备 Nround(Tf-TO)(h*Ll)； %总输出点数 for il: N； %每循环一次，输出一点数据 for j1: L1； %每输出点之间计算Ll次 (调用四阶龙格-库塔法程序块，求解 end； yy,Y(nout)； %保存输出环节动态响应值 tt,t(i)+h*L1； %保存时间向量 end ； % 输出结果 t,y %输出响应数据 plot(t,y) %输出响应图形 同样，在以上程序语句中用到了以下MATLAB语言提供的特殊功能函数： diag( )； %产生以括号内向量元素作为主对角元素的对角阵 P(:,I)； %取P阵中第I列的所有元素作为单独列向量,若需要观察其他环节动态响应，可重新输入 值，再次运行程序，即可得到另一环节输出变化曲线、数据。,面向结构图的数字仿真,4计算步长的选择 程序中采用固定步长计算方法，即计算步长 是固定不变的，这样计算过程简便，误差也在工程设计允许范围之内。 通常可按以下经验数据选择四阶龙格-库塔法的步长值。,为系统开环频率特性的剪切频率；,为系统阶跃响应的上升时间；,为系统阶跃响应的调节时间(过渡过程时间)。,或,面向结构图的数字仿真,若系统中有局部闭环，则以上各值应按响应速度最快的局部闭环考虑。当 、 大致能估计范围时，以上经验数据意味着，为充分反映系统响应开始阶段变化较快的情况，在 内至少应计算10个点；或为全面反映系统响应整个过渡过程变化情况，在 内至少应计算40个点。 此外，还应注意到， 的选取应小于系统中最小时间常数 的两倍，即 以保证数值计 算的稳定性，得到较可靠的结果。,连续系统的离散相似法,差分方程是描述离散(采样)控制系统的数学模型，其主要特点是方程中各变量由各相邻时刻的变化量(差分及高阶差分关系)制约，方程一旦列出，就相当于得到递推方程。从初始时刻开始，可以递推求出各离散时刻的状态变量值。因此，在一定条件下，将连续系统的微分方程形式的状态方程转化为差分方程形式的状态方程组，就可以达到不必采用数值积分求解的目的。 连续系统状态方程的离散化方法在控制理论课程中已有详细叙述，这里从数字仿真的角度简要重述如下。,连续系统的离散相似法,1连续系统状态方程的离散化 设连续系统状态方程为,为状态初始值,(4.31),则由现代控制理论基础可知，状态变量,的解为,即,(4.32),连续系统的离散相似法,这就是控制理论中介绍的线性时不变系统的运动方程， 称为状态转移矩阵。若将该系统进行离散化处理，并取相邻两采样时刻为KT和(k+1)T，则其状态变量的值为,(4.33),(4.34),为使此关系式能用递推形式给出，用式(4.34)- *式(4.33)得,连续系统的离散相似法,取 进行变量代换，则,若采用零阶保持器，在相邻两个采样时刻输入信号保持不变， 即 ， ，如图4.17所示。则此时,(4.35),连续系统的离散相似法,图4.17 两采样时刻之间 的变化,若希望递推公式精度更高些，应该考虑到在两次采样时刻kT和(k+1)T之间 一直在变化，用一阶保持器近似更为合理，如图4.17所示。,连续系统的离散相似法,而,所以,(4.36),于是将式(4.36)代入式(4.35)中得,连续系统的离散相似法,已知,令,(4.37),则采用零阶保持器时有,(4.38),采用一阶保持器时有,(4.39),连续系统的离散相似法,利用MATLAB语言控制系统工具箱中的功能函数，可使离散化系数矩阵 、 、 的求取过程大为简化。若已知连续系统状态方程各阵模型参数(A、B、C、D)以及采样周期T，则语句,返回的矩阵F、G就是所要求的 、 。,如果考虑精度高一些的、输入加了一阶保持器的算法，则在求得F、G后，再用一条组合语句,得到的矩阵H就是所要求的 。语句中所用的求取公式为,=,是对式(4.37)通过分步积分得到的。,连续系统的离散相似法,MATLAB中还提供了功能更强的求取连续系统离散化矩阵的函数c2dm( )。它与c2d ( )的区别在于：允许使用者在调用时自行选择确定离散化变换方式(见第3章)，并且所得到的是标准的离散化状态方程,式中，各系数阵(Ad,Bd,Cd,Dd)也可由语句来求得。,语句调用格式如下：,选项)；,连续系统的离散相似法,与其他转换方式类似，语句,选项)；,是离散化过程的逆过程，它是用以完成从离散化系统转换为连续系统各系数矩阵求取过程的功能函数。有关这些语句使用的详细情况请参阅有关参考文献。,连续系统的离散相似法,2典型环节的离散化 对阶次较高的系统，系统离散化后求取 、 、 并非易事，遇到复杂连接的系统结构图时，欲采用上述方法，事先需要设法化简、合并，还要做大量的工作，这就失去了计算机仿真运算的意义。因此，下面考虑如何把典型环节连续模型化为离散模型，使离散化仿真模型也能面向复杂连续系统。,连续系统的离散相似法,仍以式(4.25)定义的典型环节为前提,即,为方便起见，略去下标i，并改写为,(4.40),令,(4.41),连续系统的离散相似法,则,(4.42),显然，由式(4.41)和式(4.42)易得对应典型环节的连续微分方程,(4.43),连续系统的离散相似法,按离散化步骤，可推得,(4.44),其中,(4.45),连续系统的离散相似法,对典型一阶环节来说，它们均为标量系数，且求法固定。只要将已知A、B、C、D各值代入式(4.45)即可求得离散化系数 、 、 以及 、 。,(1) 积分环节。,，则,A=0,B=1,C=K,D=0,连续系统的离散相似法,所以，状态与输出递推公式为,连续系统的离散相似法,(2) 比例积分环节。,A=0,B=1,C=K,D=bK,由于A，B，C均与积分环节相同，故 、 、 和 与积分环节完全相同，相应状态方程也完全相同。,但因 ，只有 不同，所以输出方程为,连续系统的离散相似法,(3) 惯性环节。,A=a,B=1,C=K,D=0 则,状态与输出递推公式为,连续系统的离散相似法,(4) 比例惯性环节：,A=a,B=1,C=K,很明显，因为A、B、C均与惯性环节相同，状态方程离散化系数 、 、也相同，故状态方程递推公式与惯性环节的状态方程完全相同。又因为 ， 、 与惯性