第2章 连续系统数字仿真的基本算法

第2章连续系统数字仿真的基本算法2.1数值积分算法2.2数值积分算法的基本分析2.3连续系统仿真的离散相似算法2.4常用快速数字仿真算法2.5实时数字仿真算法小结2.1数值积分算法2.1.1

数值积分算法的基本原理2.1.2

欧拉法2.1.3

龙格-库塔法2.1.4

微分方程数值积分的矩阵分析上一页下一页返回2.1.1数值积分算法的基本原理连续系统仿真的数值积分算法就是利用数值积分法将常微分方程（组）描述的连续系统变换成离散形式的仿真模型——差分方程（组）。为了能在计算机上进行求解，首先要把被仿真系统的数学模型表示为一阶微分方程组或状态空间模型。上一页下一页返回假设一阶向量微分方程及初值问题为可以写成一阶微分方程组和初值问题形式上一页下一页返回(2.1)(2.2)在上的解析解为希望用公式来代替解析解，其中，，分别是，，的近似解。

上一页下一页返回

(2.3)(2.4)所谓数值解法就是寻求初值问题的解在一系列离散时间点上的近似解。相邻两个离散时间点的间距称为步长或步距。数值积分法的主要问题归结为如何对向量函数进行数值积分，求出在区间上定积分的近似解

，并且数值方法的共同特点是步进式的（即所谓递推式的）。上一页下一页返回常用的基本算法可以分为单步法、多步法和预估-校正法三大类，而每一类算法又可以分为显式算法和隐式算法两种类型。先从标量形式开始讨论。考虑一阶微分方程及初值问题解析解为希望用公式来代替解析解，其中分别是,,的近似值。

上一页下一页返回

(2.5)

(2.6)

(2.7)2.1.2欧拉法欧拉（Euler）法是一种最简单的数值积分算法，对于

在区间上求积分，有若区间足够小，则上的可以近似地看成是常数上一页下一页返回(2.8)即有欧拉法为式中上一页下一页返回

(2.9)

(2.10)欧拉法的图形表示：为图2.1中的非曲面面积actk+1tk，欧拉法用矩形面积abtk+1tk代替。显然，欧拉法仅适用于步长h很小的场合。上一页下一页返回图2.12.1.3龙格-库塔法龙格-库塔（Runge-Kutta）法是求解常微分方程初值问题(2.5)式的各种数值积分算法中应用得最广泛的一种，包括许多不同的公式。思路：用若干个时间点上f

的函数值的线性组合来代替f的各阶导数项，然后按泰勒公式展开确定其中的系数。上一页下一页返回

RK法包含有显式、隐式和半隐式等算法。仅介绍显式RK法。一般形式：式中，

Wi为待定的权因子；r为使用的k值的个数；ki为不同时间点上导数f的值；ci，aij为待定系数。上一页下一页返回

(2.13)当r=1时，只有一个k1，就得到了欧拉法当r=2时，将代入中，并将在附近用泰勒级数展开,可得上一页下一页返回(2.15)(2.14)于是，有将其与泰勒级数式逐项进行比较，令h，h2项的对应系数相等，得到以下关系取c2=1，有上一页下一页返回(2.16)(2.17)从而，有RK2法当r=3时，可得RK3法上一页下一页返回(2.18)(2.19)当r=4时，可得到如下著名的四阶RK法（亦称为经典RK法，简记为RK4法）。该算法是数字仿真中最常用的一种算法。计算量较大，但计算精度较高。在实际仿真中得到了广泛的应用。上一页下一页返回

(2.20)几种数值积分法都可以看成是在附近用泰勒级数展开而产生的，只不过是泰勒级数所取项数的多少不同而已。欧拉法只取前两项，

RK2法取了前三项，RK4法取了前五项。从理论上讲，可以构造任意高阶的计算方法。但是，精度的阶数与计算函数f值的次数之间的关系并非等量增加的。上一页下一页返回由此可见，RK4法有其优越性。每步计算f的次数234567r≥8算法精度阶数234456r-2表2.1f的计算次数与算法精度阶数的关系

上一页下一页返回2.1.4微分方程数值积分的矩阵分析对于一阶向量微分方程及初值问题RK4法的向量形式为上一页下一页返回

(2.21)

(2.22)具体写出来就是上一页下一页返回(2.23)在控制系统的仿真中，最常见的向量微分方程是线性定常系统的状态方程即有RK4法的4个可表示为上一页下一页返回

(2.24)(2.25)(2.26)2.2数值积分算法的基本分析2.2.1

单步法和多步法2.2.2

显式算法和隐式算法2.2.3

截断误差和舍入误差2.2.4

数值积分算法的计算稳定性2.2.5

数值积分算法的计算精度、速度、稳定性与步长的关系2.2.6

数值积分算法的选择原则2.2.7

误差估计与步长控制2.2.8

数值积分算法仿真实例上一页下一页返回2.2.1单步法和多步法上一页下一页返回上节中介绍的几种数值积分算法都有一个共同特点：在计算yk+1时只用到了yk

,而不直接用yk-1,yk-2

，…

等项的值，即在本次计算中仅仅用到了前一次的计算结果。这类算法称为单步法。

单步法运算有如下优点：需要存储的数据量少，占用的存储空间少；只需要知道初值，就可以启动递推公式进行运算，即可以自启动；容易实现变步长运算。

与单步法相对应的还有一类数值积分算法，称为多步法。在它的计算公式中，本次计算不仅要用到前一次的计算结果，还要用到更前面的若干次结果。例如AB4法:式中多步法与单步法相比，欲达到相同的精度，计算工作量要少得多。因此，在相同条件下多步法比单步法要快。多步法的主要缺点是，它不是自启动方法，必须用其它方法求初始几步的值。另外，多步法不容易实现变步长运算。上一页下一页返回

(2.28)2.2.2显式算法和隐式算法如果数值积分算法在计算yk+1时所用到数据均已求出，则称其为显式算法。

如果数值积分算法的计算公式右端中隐含有未知量yk+1，则称其为隐式算法。例如AM4法就是隐式算法:上一页下一页返回(2.29)隐式算法也不是自启动的算法，需要用另一个显式算法估计一个初值，然后再用隐式算法进行迭代运算，这就构成了预估-校正算法。

四阶Adams预估-校正算法：上一页下一页返回

(2.30)2.2.3截断误差和舍入误差在分析数值积分算法的精度时，通常用泰勒级数作为工具。假设前一步求得结果是准确的，即有则用泰勒级数求得在处的解析解为

不同的数值积分算法相当于在上式中取了不同项数之和而得到的近似解。上一页下一页返回(2.31)(2.32)这种由数值积分算法单独一步引进的附加误差称为局部截断误差。它是数值积分算法给出的解与微分方程的解析解之间的差，故又称为局部离散误差。步长h越小，局部截断误差就越小。若数值积分算法的局部截断误差为，则方法为r阶的。算法的阶数可以作为衡量算法精确度的一个重要标志。上一页下一页返回以上的分析是在假设前一步所得结果是准确的前提下得出的，即成立时，有但在求解过程中，实际上只有当k=0时，才成立。而当k=1,2,3,…时，并不成立。因而r阶算法求得的解的实际误差要大一些的。设yk是在无舍入误差情况下由r阶方法算出的微分方程式的近似解，y(t)为微分方程的解析解，为算法的整体截断误差，则可以证明，。即整体截断误差比局部截断误差低一阶。

上一页下一页返回(2.33)2.2.4数值积分算法的计算稳定性1.

计算稳定性的概念2.

欧拉法的计算稳定性3.

龙格-库塔法的计算稳定性4.

多步法的计算稳定性上一页下一页返回1.计算稳定性的概念连续系统的数字仿真，实质上就是将给定的微分方程（组）变换成差分方程（组），然后从初值开始，逐步进行递推计算。

【例2.1】已知微分方程及初值试比较在取不同步长时，其精确解与数值积分算法解之间的差异。上一页下一页返回

【解】该初值问题的解析解为上一页下一页返回对应的结果曲线如图2.2（a）所示。图2.2

(a)解析解取，分别用欧拉法和RK4法计算处的，所得结果曲线如图2.2

（b）—（g）所示。

图2.2

(b)欧拉法（h=0.1）图2.2

(c)RK4法（h=0.1）上一页下一页返回图2.2

(d)欧拉法（h=0.075）

图2.2

(e)RK4法（h=0.075）

图2.2

(f)欧拉法（h=0.05）

图2.2

(g)RK4法（h=0.05）

上一页下一页返回

从图中可以看出，解析解单调下降并迅速收敛到0。当时，欧拉法和RK4法的解曲线均发散，数值积分算法的解是错误的。当时，欧拉法的解曲线仍然发散，对应的解是错误的；RK4法的解曲线单调下降并收敛到0，对应的解是正确的。当时，欧拉法和RK4法的解均收敛到0（虽然欧拉法的解曲线是振荡收敛的）。如果只要求得到处的的解，则两种数值积分算法的解都可以认为是正确的。事实上，此时有

上一页下一页返回欧拉法RK4法解析解为什么会出现上述情况呢？这是因为数值积分算法只是一种近似积分方法，在反复的递推计算中会引进误差（这种误差通常是由初始数据的误差及计算过程中的舍入误差产生）。如果误差的累积越来越大（见【例2.1】）,将使计算出现不稳定，从而得到错误的结果。系统的稳定性与计算稳定性是两个不同的概念。由于选用的数值积分算法不同，即使对于同一系统，差分方程也各不相同，计算稳定性也就各不一样。上一页下一页返回通常用一个简单的一阶微分方程来考查数值积分算法的计算稳定性。微分方程及初值问题称为测试方程（TestEquation），其中，为方程的特征根。根据稳定性理论，当特征根位于左半s平面，即时，原系统稳定。上一页下一页返回

(2.34)2.欧拉法的计算稳定性用欧拉法对测试方程进行计算，有假定仅仅在初始值y0处引进了初始误差，而在递推计算过程没有引进任何其它误差。显然，此时的误差

仅由的误差引起的，所以有用（2.36）式减去（2.35）式，得到误差方程（2.35）（2.36）上一页下一页返回以此类推，有当时，。表明若在递推计算过程中的某步引入了误差，则随着递推步数的增加，这个误差将逐渐扩大，最终导致差分方程的解完全失真。反之，当时，。表明随着递推步数的增加，引入的误差会被逐渐减小或保持有界。在这种情况下，称差分方程是计算稳定的。上一页下一页返回

(2.37)

(2.38)显然，对于欧拉法，合理地选择步长h使其满足是保持其计算稳定性的充要条件。令，由，可知欧拉法在实轴上的稳定区域为（-2，0），即有因此，越大，步长h就应该取得越小。为了保证计算稳定性而对步长h有所限制的数值积分算法称为条件稳定算法。在任何步长h>0的情况下计算都稳定的数值积分算法称为绝对稳定算法（亦称为无条件稳定算法或恒稳算法）。欧拉法是一种条件稳定的计算方法。上一页下一页返回3.龙格-库塔法的计算稳定性将RK法应用于测试方程，容易得到下列误差方程令，代入上式，可得该式稳定的条件为可得如表2.2所示的各阶RK法在实轴上的稳定区域。上一页下一页返回（2.40）(2.39)显式RK法都是条件稳定算法，步长h的大小除了与所选用的算法有关外，还与方程本身的性质有关。对于实际系统，由于其特征值不一定为实数，因此满足（2.40）式的也是复数。一般而言，步长h必须满足不等式

(2.41)

r所在的实负稳定区域1（-2，0）2（-2，0）3（-2.51，0）4（-2.785，0）

表2.2RK法在实轴上的稳定区域上一页下一页返回4.多步法的计算稳定性 将多步法应用于测试方程，可以得到类似的误差方程。

结论（见表2.3）：在同阶的多步法中，隐式算法的稳定区域比显式算法的大得多；随着阶数的增加，多步法的稳定区域逐渐减小。上一页下一页返回方法的阶数1234显式（-2，0）（-1，0）（-6/11，0）（-3/10，0）隐式（-∞，0）（-∞，0）（-6，0）（-3，0）表2.3Adams法在实轴上的稳定区域

除了AM1法和AM2法（隐式Adams法）为绝对稳定算法外，其它算法都是条件稳定算法。这就是说，步长h必须满足不等式式中，Ｍ为由积分算法确定的常数。上一页下一页返回|λh|<M(2.42)2.2.5数值积分算法的计算精度、速度、稳定性与步长的关系各种数值积分算法对应的差分方程是对原微分方程的近似，存在着明显的截断误差。并且，由于计算机的字长有限，计算只能限制在有限的位数，将引入舍入误差。步长取得过大了会带来较大的截断误差，甚至会出现计算不稳定现象；步长取得过小了又会增加计算步数，而使舍入误差累积，总的误差反而变大。所以，仿真的总误差与步长的关系不是单调函数关系，而是一个具有极值的函数关系，如图2.3所示。

上一页下一页返回返回上一页下一页返回图2.3

【例2.2】考虑如下二阶系统在s上的数字仿真，并将不同步长下仿真结果与解析解进行精度比较。

【解】采用RK4法进行计算，步长选取为0.001s到1s之间变化，并且将计算结果与解析解进行比较，求出最大误差数据列于表2.4中。步长h0.0010.010.020.050.10.20.51.0R=0最大误差0.007970.000780.000380.000210.001040.015010.594709.48049R=0.5最大误差0.001040.000100.00050.000010.000030.001070.016450.79498上一页下一页返回表2.4从表2.4中可以看出，当步长h=0.05时，总误差最小；当步长h小于0.05时，由于舍入误差变大而使总误差增加；当步长h大于0.05时，则由于截断误差的增加也使得总误差加大。另外，当系统的解变化激烈时（如R=0），误差对步长的变化较为敏感；当系统的解变化平稳时，步长的变化对误差的影响就要缓和得多。数值积分算法确定以后，在选择步长时，需要综合考虑。对于RK4法，步长的选择应满足上一页下一页返回

(2.43)对于一般工程系统的分析和设计而言，仿真结果的误差不超过0.5%就可以满足精度要求了。经验表明，对于RK4法，通常可以根据系统方程中的最小时间常数来选择步长，一般取或者根据系统中反应最快的小闭环的开环剪切频率来选择，取上述两种选择步长的经验公式都是针对RK4法而言的。如果想要用RK2法进行仿真，则为了达到相同的精度，步长一般应取为RK4法的1/10。上一页下一页返回(2.44)(2.45)如果采用欧拉法，则为了达到相同精度，步长还要取得更小。显然RK4法的计算效率较高，这就是控制系统数字仿真中通常选用RK4法的重要原因之一。上一页下一页返回2.2.6数值积分算法的选择原则正确选择数值积分算法是一个十分重要而复杂的问题。到目前为止还没有一个普遍适用的规律，这里仅给出一些选择时应考虑的原则。精度要求 当需要进行高精度仿真时，可以采用高阶的多步隐式算法和较小的步长。经验表明，低精度问题最好用低阶算法来处理。如果选用高阶算法，则应在保证计算稳定的前提下，把步长取得大一些。上一页下一页返回计算速度 计算速度取决于在给定积分时间内的计算步数和每步计算所需的时间。在右端函数复杂、计算量大而精度要求又高的问题中，可以采用Adams预估-校正法。为了加快计算速度，在积分算法选定后，应在保证精度的前提下，尽量选择较大的步长，以减少计算步数。

计算稳定性 保证数值解的计算稳定性，是进行数字仿真的先决条件。从计算稳定性的角度来看，同阶的RK法优于显式Adams法，但又不如隐式Adams法。最好避免使用显式Adams法。上一页下一页返回自启动能力 单步法具有自启动能力。多步法没有自启动能力。一般简单的仿真程序多采用单步法。

步长变化能力 如果要求仿真时步长可变，最好选用单步法。

数值积分算法的选择与多种因素有关，要根据具体情况而定。经验表明，当微分方程（组）右端函数的计算量不大而精度要求又不很高时，RK4法是很好的选择。在控制系统的数字仿真中，最常见的是线性系统的状态方程，右端函数中线性部分一般占了绝大部分。因此，对于一般控制系统的仿真，通常都采用RK4法。

上一页下一页返回2.2.7误差估计与步长控制一个成熟的仿真程序或仿真语言通常将步长的自动控制作为必要的手段。实现步长自动控制的前提是要有一个好的局部截断误差估计公式。为了估计一种RK法的误差，通常采用的方法是设法找到另一个低阶（一般是低一阶）的RK公式，要求这两个公式中的ki相同，则两个公式计算结果之差就可以作为误差估计。上一页下一页返回第一个四阶变步长RK法是Merson在1957年给出的，称之为龙格-库塔-默森法（简称为RKM34法），四阶RKM公式为（2.46）上一页下一页返回三阶RKM公式为式中，由（2.46）式定义。用作为误差估计，即有该算法可得到四阶精度的计算结果，三阶精度的误差估计，是目前被广泛应用的一种数值积分算法，它在实轴上的稳定域为（-3.548，0）。其缺点是计算量较大，每步需要计算5次f。上一页下一页返回(2.48)

(2.47)与RKM34法类似的还有E.Fehlberg在1969年提出的RKF45法（即龙格-库塔-费尔别格法），它每步计算6次f，获得五阶精度的计算结果，四阶精度的误差估计，被公认为是对非病态系统进行仿真最为有效的算法之一。1978年L.F.Shampine提出了RKS34法（即龙格-库塔-夏普勒法），它每步只计算4次f却能获得四阶精度的计算结果，三阶精度的误差估计。这两种算法的稳定性与RK4法相差不大。上一页下一页返回步长控制策略是“加倍-减半法。”设定最小允许误差和最大允许误差，当估计的局部误差大于最大允许误差时，将步长减半，并重新计算本步；当误差在最大允许误差和最小允许误差之间时，步长不变；当误差小于最小允许误差时，将步长加倍。每步的局部误差通常取以下形式上一页下一页返回

(2.49)对于向量微分方程，误差估计为并将ek修改成步长控制策略可以表示为上一页下一页返回

(2.50)(2.51)

(2.52)2.2.8数值积分算法仿真实例在MATLAB/Simulink环境下，利用数值积分算法对系统进行仿真的途径有两种。对于以微分方程给出的数学模型，通常用MATLAB语言编程实现；而对于控制系统分析和设计中最常见的以系统结构图描述的数学模型，则采用Simulink实现更为方便。上一页下一页返回1.

数值积分算法的编程实现2.

Simulink

模型的构建与仿真上一页下一页返回1.数值积分算法的编程实现 用MATLAB语言编程实现仿真的主要步骤是调用MATLAB中的ODE（OrdinaryDifferentialEquation）解函数。MATLAB提供的常用ODE解函数如下：ode45

此算法被推荐为首选算法；ode23

这是一个比ode45低阶的算法；ode113

用于更高阶或大的标量计算；ode23t

用于解决难度适中的问题；ode23s

用于解决难度较大的问题；ode15s

与ode23相同，但要求的精度更高；ode23tb

用于解决难度较大的问题。上一页下一页返回这些ODE解函数的调用格式基本相同。例如，ode45的基本调用格式为

[t,x]=ode45('方程函数名',tspan,x0,tol)其中，方程函数名为描述系统状态方程的M函数名称，tspan一般为仿真时间范围（例如，取tspan=[t0,tf]，t0为起始计算时间，tf为终止计算时间）；x0为系统状态变量初值，应使该向量元素个数等于系统状态变量的个数；tol用来指定的精度，其默认值为10-3（即0.1%的相对误差），一般应用中可以直接采用默认值。函数返回两个结果t向量和x阵。由于计算中采用了步长自动控制策略，因而t向量不一定是等间隔。但是，仿真结果可以用plot（t,x）指令绘制出来。上一页下一页返回

【例2.3】某地区某病菌传染的系统动力学模型为式中，x1表示可能受到传染的人数，x2

表示已经被传染的病人数，x3表示已治愈的人数。试用ode45编程，对其进行仿真研究，并绘制出对应的时间响应曲线。上一页下一页返回【解】采用MATLAB语言进行编程，文件名为exam2_3.m。程序如下：%这是例2.3的仿真程序clearx0=[620,10,70]; %置状态变量初值tspan=[0,30]; %置仿真时间区间[t,x]=ode45('fun2_3',tspan,x0); %调用ode45求仿真解plot(t,x(:,1),t,x(:,2),‘--’,t,x(:,3),‘:’); %用不同的线型绘制仿真结果曲线xlabel('time(天)t0=0,tf=30'); %对t-x轴进行标注ylabel('x(人):x1(0)=620,x2(0)=10;x3(0)=70');legend('x1','x2','x3');grid;上一页下一页返回其中，fun2_3是保存在名为fun2_3的M文件中的M函数：functionxdot=fun2_3(t,x)xdot1(1)=-0.001*x(1)*x(2); %第一个微分方程xdot1(2)=0.001*x(1)*x(2)-0.072*x(2); %第二个微分方程xdot1(3)=0.072*x(2); %第三个微分方程xdot=xdot1';运行exam2_3.m后，得到如图2.4所示结果曲线。上一页下一页返回返回上一页下一页返回图2.42.Simulink

模型的构建与仿真 对于以系统结构图描述的数学模型，采用Simulink构模并仿真是十分方便的。在大多数情况下，用户甚至无需编制一条仿真程序就可以进行仿真。上一页下一页返回Simulink

模型的构建仿真参数的设置和ODE算法的选择Simulink仿真结果输出仿真结构的参数化与M函数的组合仿真上一页下一页返回Simulink

模型的构建为了能用Simulink对系统进行仿真，首先要在Simulink环境下打开一个空白模型窗口；然后依据系统结构图给定的环节和信号流程，从Simulink模块库的各个子库中选择相应的模块，并用鼠标左键将它们拖入模型窗口；双击选择的模块，设置需要的参数；对各模块进行连接，这就构成了需要的Simulink模型（即仿真结构图，参见1.6.2节）。上一页下一页返回仿真参数的设置和ODE算法的选择在模型窗口中（见下页图）的Simulation菜单下选择其中的仿真参数子菜单（SimulationParameters），就会弹出一个仿真参数对话框，如图2.5所示。上一页下一页返回图2.5返回上一页下一页返回在仿真参数对话框中有5个标签，默认的标签为微分方程求解程序Solver的设置，在该标签下的对话框主要接受微分方程求解算法及仿真参数的设置。算法的选择

通过该对话框可以由Solveroptions栏目选择不同的求解算法，定步长（Fixed-step）下支持的算法如图2.6(a)所示，变步长（Variable-step）下支持的算法如图2.6(b)所示。一般情况下，连续系统仿真应该选择ode45变步长算法或者定步长的ode4（即RK4）算法，而离散系统一般默认地选择定步长的discrete（nocontinuousstates）算法。上一页下一页返回返回

图2.6(a)图2.6

(b)上一页下一页返回计算步长的选择定步长算法的计算步长应该由Fixedstepsize框填入参数进行选择，一般可以选择auto；而变步长算法建议步长范围使用auto选项。仿真时间的设置 在仿真参数对话框的相关对话框中还可以修改仿真的初始时间和终止时间。上一页下一页返回Simulink仿真结果输出在完成了仿真参数的设置和ODE算法的选择后，就可以启动仿真。Simulink会自动将系统结构图转换成状态空间模型并调用所选择的算法进行计算。为了得到所需要的仿真结果，除了可以直接采用Scope模块显示仿真结果曲线外，还可以将仿真结果数据传送到MATLAB工作空间中，利用plot指令绘制相应的图线。上一页下一页返回

【例2.4】考虑如图2.7所示的直流电机拖动系统，试研究外环PI控制器对系统阶跃响应的影响。上一页下一页返回图2.7【解】构建系统的Simulink模型，如图2.8所示。启动仿真后，可以立即得出仿真结果，并自动返回到MATLAB工作空间中，其中时间变量名默认为tout，输出信号的变量名默认为yout。在MATLAB命令窗中运行指令>>plot(tout,yout,’k’);grid;上一页下一页返回图2.8(exam2_4.mdl)得到如图2.9(a)所示的阶跃响应曲线。显然，该曲线不是很理想，超调量较大。为此，可以将外环的PI控制器参数调整为，并分别选择α=0.17，0.5，1，1.5，可以得到如图2.9

(b)所示的仿真结果。可以看出，如果选择PI控制器为，就能够得到较为满意的控制效果。图2.9

（a）图2.9

（b）上一页下一页返回

【例2.5】考虑著名的VandePol方程试绘制其相轨迹（）。【解】选择状态变量则有如下非线性状态方程上一页下一页返回第1个方程可以看成是将x2(t)信号作为一个积分器的输入端，积分器的输出将成为x1(t)信号。同样，x2(t)信号本身也可以看成是一个积分器的输出，而积分器的输入端信号应该为。于是，利用Simulink提供的各种模块可以得到如图2.10所示的仿真结构。上一页下一页返回图2.10(exam2_5.mdl)

Simulink模型中有些模块需要将输入端和输出端（通常用于反馈路径）掉换一下方向。为此，可以用鼠标单击需要掉换方向的模块选中它，则选中的模块的四个角出现黑点，表明它处于选中的状态。然后，打开Simulink的Format菜单(见图2.11)，选择其中的翻转子菜单（Flipblock）即可。上一页下一页返回图2.11启动仿真后，仿真结果将赋给MATLAB工作空间内的变量t，x1和x2。在MATLAB命令窗口中运行指令>>plot(x1,x2,’k’);grid;

得到相轨迹，如图2.12所示。上一页下一页返回图2.12仿真结构的参数化Simulink模型中的参数可以是实际数值，也可以是用字母表示的变量名。用字母表示的变量可以在MATLAB的工作空间中赋值，或用M文件赋值，然后进行仿真计算。上一页下一页返回

【例2.6】含有磁滞回环非线性环节的控制系统如图2.13

(a)所示，其中，磁滞回环非线性环节的特性如图2.13（b）所示。试研究该非线性环节对系统性能的影响。上一页下一页返回图2.13

【解】构建系统的Simulink模型，如图2.14所示。为了便于研究问题的方便起见，将Backlash模块的参数deadbandwidth设置为C1。在MATLAB命令窗口依次运行C1的不同值(C1=0.5,1,2)的指令后启动仿真，并在仿真结束后在MATLAB命令窗口中运行指令>>plot(tout,yout,’k’);

上一页下一页返回图2.14(exam2_6.mdl)得到如图2.15所示的仿真曲线。显然，不同磁滞宽度的阶跃响应曲线的形状不相同。需要注意的是，为了将不同C1值对应的阶跃响应曲线绘制在同一坐标轴下，在第一次绘制图形后，应该在MATLAB命令窗口中运行指令>>holdon上一页下一页返回图2.15(r(t)=2∙1(t))此外，如果输入的幅值增大或减小，则原来系统响应曲线的形状将可能出现不同。图2.16为C1=1时，输入幅值等于3和0.6的仿真结果。显然，这一点与线性系统不相同。上一页下一页返回图2.16与M函数的组合仿真

如果Simulink模型中存在复杂的非线性环节或复杂的逻辑运算，而在MATLAB提供的所有工具箱中都找不到相应的模块，可以自己编制一个M函数，嵌入Simulink模型中。上一页下一页返回

【例2.7】

某非线性系统如图2.17所示，试求r(t)=2.1(t)时系统的动态响应。

【解】构建系统的Simulink模型，如图2.18所示。为了便于研究问题的方便起见，不采用Discontinuities子库中的Saturation模块，而选择User-Defined-Functions子库中的MATLABFcn模块,并将参数MATLABFunction设置为：saturation_zone。上一页下一页返回图2.17返回上一页下一页返回图2.18(exam2_7.mdl)然后编制M函数内容如下：%saturation_zonefunctionfunction[uo]=saturation_zone(ui)ifui>=1

uo=1;elseif

ui<=-1

uo=-1;else

uo=ui;end

启动仿真后，得到如图2.19所示仿真结果。上一页下一页返回图2.192.3连续系统仿真的离散相似算法各种数值积分算法的公式都是一些差分方程，仅仅计算了在各个离散时间点处系统的近似解。这相当于对连续系统进行了离散化处理，把原来的连续系统模型近似等价为一个离散系统模型。从本质上讲，连续系统的数字仿真就是要找出一个与该系统等价的离散模型。数值积分算法也是离散化算法。上一页下一页返回本节将从连续系统离散化的角度出发，建立连续系统模型的等价离散化模型，并用采样系统的理论和方法介绍另一种常用的仿真算法。这种算法使得连续系统在进行（虚拟的）离散化处理后仍保持与原系统“相似”，故称之为离散相似算法。与数值积分算法相比，离散相似算法的每步计算量要小得多，稳定性也要好得多，因而允许采用较大的计算步长。然而，它通常只适合线性定常系统的仿真，具有一定的局限性。考虑到应用范围与实现的便利性，本节仅介绍时域离散相似算法。

上一页下一页返回2.3.1

时域离散相似算法的基本概念2.3.2

时域离散化模型的推导2.3.3

时域离散化模型的特性分析2.3.4

时域离散算法仿真实例上一页下一页返回2.3.1时域离散相似算法的基本概念连续系统的离散化是在原连续系统的输入、端出端分别加上虚拟的采样开关，使输入、输出信号离散化，此时系统模型即为离散化模型。为了使输入信号u(t)离散化后仍能保持原来的变化规律，在输入采样开关后设置信号保持器（亦称为信号重构器），复现原输入信号u(t)，其结构如图2.20所示。

上一页下一页返回图2.20当连续系统用状态空间模型表示时，在状态方程的输入、输出端分别加入虚拟的采样开关，并在输入采样开关后加上保持器，如图2.21所示。

上一页下一页返回图2.212.3.2

时域离散化模型的推导设系统的状态空间模型为

（2.53）

（2.54）对状态方程进行拉氏变换，得

（2.55）从而，有

（2.56）对（2.56）式进行拉氏反变换，并利用卷积积分，

（2.57）式中

上一页下一页返回（2.58）设采样周期为T，考察kT及(k+1)T处两个相邻采样时刻状态向量x(t)的值。将t=(k+1)T代入（2.57）式，得由（2.57）式，有于是，得

（2.59）对上式右端积分项中的输入向量u(t)作不同的处理，可以得到不同的时域等价离散化模型。

上一页下一页返回于是，（2.59）式变为上一页下一页返回假定对输入u(t)采用零阶保持器，如图2.22所示，则有图2.22取t=τ-kT进行变量代换，则

（2.60）式中

（2.61）这样，就得到了采用零阶保持器的等价离散化模型

（2.62）

（2.63）

上一页下一页返回如果对输入u(t)采用三角形保持器，即假定在两个相邻采样时间之间u(t)按u(kT)与u[(k+1)T]连线斜率变化，有由于故代入（2.59）式，整理可得

（2.64）令

（2.65）上一页下一页返回则（2.64）式可以写为（2.66）这是采用三角形保持器时的离散化状态方程。

由于在（2.62）式和（2.66）式中，当采样周期T选定之后，以及都是常数矩阵，可以一次性计算出来，不再变化。这样，利用（2.62）式或（2.66）式，在已知各状态变量初值的情况下，可以很容易地求出在各采样时刻处各状态变量的值，进而，求出相应的输出方程的值。

上一页下一页返回2.3.3时域离散化模型的特性分析

将离散化状态方程（2.62）式或（2.66）式代入测试方程（2.34）式，有

而测试方程的解析解有如下关系

由此可见，时域离散相似算法是绝对稳定算法。上一页下一页返回由图2.22可知，在输入端采样开关后加入零阶保持器相当于对输入u(t)作了零阶近似；加入三角形保持器相当于对输入u(t)作了一阶近似。因此，（2.62）式和（2.66）式的局部截断误差分别是和，分别相当于欧拉法和RK2法的精度。

对于单输入-单输出系统的状态方程，（2.62）式和（2.66）式每步的计算量主要是1次矩阵与向量的乘法（即n2次乘法），相当于欧拉法的每步计算量。

上一页下一页返回2.3.4时域离散算法仿真实例

MATLAB的控制系统工具箱中提供了实现连续系统和离散系统之间相互转换的函数，最常用的是将连续系统转换成离散系统的函数c2d和c2dm，其调用格式分别为[Ad,Bd]=c2d(A,B,T)[Ad,Bd,Cd,Dd]=c2dm(A,B,C,D,T,‘method’)c2d将用矩阵A、B描述的连续系统按照（2.62）式转换成离散系统，T为离散化步长。c2dm将用矩阵A、B、C、D描述的连续系统按照指定的离散化方法转换成离散系统，其中method为转换时所选用的离散化方法，其选项主要有以下几种：上一页下一页返回

zoh

假设对输入信号加一个零阶保持器；

tustin

采用双线性变换法（Tustin法）

prewarp

采用改进的Tustin法；

matched

采用根匹配法。

c2dm函数还可以实现连续系统传递函数到离散系统脉冲传递函数之间的转换，调用格式为

[numd,dend]=c2dm(num,den,T,‘method’)用于根据method所指定的转换方法，将由num和den描述的传递函数G(s)转换成由numd和dend描述的脉冲传递函数G(z)

。上一页下一页返回

【例2.8】考虑如图2.23所示的非线性控制系统，试用时域离散相似算法研究其单位阶跃响应。

【解】图中的传递函数对应的状态空间模型为在该模型前加上虚拟的采样开关和零阶保持器，如图2.24所示。

上一页下一页返回图2.23相应的仿真模型如图2.25所示。上一页下一页返回对应的离散化状态空间模型为式中图2.24图2.25为了避免手工计算和，采用如图2.26所示的Simulink模型，并将DiscreteState-Space模块的四个参数Parameters分别设置为：a1，b1，c1，d1，Sampletime为：0.05。在MATLAB命令窗口运行下列指令>>a=[00;1-1];b=[1;0];c=[01];d=0;T=0.05;>>[a1,b1,c1,d1]=c2dm(a,b,c,d,T,‘zoh’);

并分别键入K=1，10，30。选择定步长的discrete（nocontinuousstates）算法，步长为0.05。启动仿真后，得到如图2.27所示的仿真结果。

上一页下一页返回图2.26(exam2_8.mdl)图2.27(a)K=1图2.27(b)K=10图2.27(c)

K=30上一页下一页返回2.4常用快速数字仿真算法在许多实际应用场合下，希望能尽可能快地获取仿真结果。这就构成了实际的快速需求和仿真速度之间的矛盾。当采用普通微型计算机对某些控制系统实施仿真时，这种矛盾往往变得非常突出，从而提出了快速仿真的需求。快速数字仿真算法的种类很多，它们各有优缺点，目前还没有一种算法能够适用于所有问题。本节将介绍常用快速数字仿真算法。上一页下一页返回2.4.1

仿真中对快速性的需求2.4.2

相匹配原理

2.4.3

替换法2.4.4

根匹配法

上一页下一页返回2.4.1

仿真中对快速性的需求通常需要采用快速仿真的问题可以分为下列几个方面：利用仿真技术进行控制系统的参数优化设计时，需要反复地对系统进行仿真计算，以获得满意的工程参数；在数学-物理混合仿真中，由于有实物系统介入仿真模型，要求仿真模型的时间比例尺与真实系统的时间比例尺完全相同（即所谓实时仿真），如果系统比较复杂或者方程个数很多，要在规定时间内完成一次系统中每一个方程的计算，就要求仿真的计算速度比较快；上一页下一页返回在复杂系统的控制中，常常需要在线用仿真方法对被控系统的状态进行预测，以确定系统的控制策略，这就要求仿真模型的时间比例尺小于系统的运行时间比例尺（即所谓超实时仿真），从而对仿真速度提出了更高的要求。一般而言，对于第一种情况，仿真速度慢意味着不能及时获得仿真结果，并且占用大量的机时，费用较高，效率低下；而对于后两种情况，仿真速度慢意味着仿真的目的达不到或者失去实际应用价值。上一页下一页返回算法的计算精度和速度是一对矛盾，在研究快速数字仿真算法时，计算速度成为矛盾的主要方面。通常在满足计算稳定性及工程精度要求的前提下，应当尽可能提高仿真计算速度。一般快速数字仿真算法有如下两点基本要求：每步计算量要小；算法具有良好的稳定性，允许采用较大的计算步长，同时又能保证必要的计算精度。应当根据实际情况，尽可能地减小仿真中附加信息的计算和处理。

上一页下一页返回2.4.2

相匹配原理相匹配的含义是，如果被仿真系统的数学模型是稳定的，则其仿真模型也应该是稳定的，并且二者的瞬态、稳态特性一致。如果对于同一输入信号，二者的输出具有相一致的时域特性，或者二者具有相一致的频率特性，则称仿真模型与原系统模型相匹配。将高阶系统的数学模型G(s)转换成与之相匹配、每步计算量较小、允许采用较大步长且具有合理精度的仿真模型G(z)，就可以利用对应的差分方程进行快速仿真。从G(s)直接推导出G(z)的方法有两种。一种称为替换法，另一种称为根匹配法。

上一页下一页返回2.4.3

替换法1．替换法的基本思想2．双线性变换法3．双线性变换法步骤4．双线性变换法特性分析5．双线性变换法仿真实例上一页下一页返回1．替换法的基本思想根据控制理论，s域和z域之间的准确映射关系为

（2.67）或

（2.68）式中，T为采样周期，亦为仿真计算时的步长。（2.68）式是超越方程，很难由得到一个关于变量u和y的线性差分方程。

上一页下一页返回2．双线性变换法ln(z)可以展开成下列的无穷级数

（2.70）取该级数的第一项，则可以得到s域与z域之间一种近似映射关系

（2.71）或

（2.72）利用（2.71）式可以把G(s)转换为G(z)。数学上称这种变换方法为双线性变换法（或Tustin法）。

上一页下一页返回双线性变换法是否满足相匹配原理呢？首先，根据相匹配原理，要求当G(s)稳定时，G(z)也应该是稳定的。将代入（2.72）式，得

即有

（2.73）由（2.73）式可知,双线性变换法将左半s平面映射到z平面的单位圆内。也就是说，如果原来系统的G(s)是稳定的，则通过双线性变换法得到的G(z)必然是稳定的。上一页下一页返回其次，双线性变换后的稳态增益不变。设连续系统的传递函数为

其稳态增益为bn/an

。若将G(s)进行双线性变换，有

显然，G(z)的稳态增益仍然是bn/an。

上一页下一页返回上一页下一页返回最后，双线性变换法具有一定的精度。设有方程（2.74）即（2.75）将（2.71）式代入上式，得（2.76）由此可得差分方程

（2.77）这相当于应用AM2法（也称为梯形法）于（2.74）式。综上所述，双线性变换法符合相匹配原理，是一种具有一定计算精度的绝对稳定算法。

双线性变换法适合于以分子分母多项式形式描述的G(s)。

上一页下一页返回【例2.9】已知连续系统的传递函数为

（2.78）试采用双线性变换法求出对应的脉冲传递函数和差分方程，并对所得结果进行分析。【解】将（2.71）式代入上式，得脉冲传递函数

（2.79）

上一页下一页返回于是,差分方程为

由（2.79）式可知，因为

所以，G(z)是稳定的。

G(s)的分子多项式为1阶，分母多项式为2阶，而G(z)的分子、分母多项式的阶次相同，均为2阶。G(s)的稳态增益为0，G(z)的稳态增益也为0。上一页下一页返回为了考虑双线性变换法所得仿真模型的精度，除了可以在时域中进行讨论外，也可以从频域的角度进行分析。事实上将，分别代入（2.78）式和（2.79）式，可得

（2.80）

（2.81）取T=1s，分别求出（2.80）式和（2.81）式的幅频特性和相频特性，如表2.5所示。

上一页下一页返回表2.5G(s)与G(z)频率特性比较

ω(rad/sec)0.10.30.60.81.01.21.5幅频特性连续系统模型0.099010.27520.44120.48780.500.49180.4615离散化模型0.09910.2770.4470.4930.4980.4760.417相频特性连续系统模型78.58°56.60°28.07°12.68°0-10.39°-22.62°离散化模型78.57°56.36°26.51°9.57°-5.06°-17.67°-33.56°上一页下一页返回3．双线性变换法步骤设连续系统的传递函数为

（2.82）采用双线性变换法求取仿真模型的步骤如下：将G(s)中的s按（2.71）式替换得到G(z)，并整理成有理分式形式

（2.83）根据得到的G(z)，由Y(z)=G(z)U(z)，两边取Z反变换，得到便于计算机递推计算的差分方程。

（2.84）

上一页下一页返回4．双线性变换法特性分析双线性变换法是一种绝对稳定的算法，可采用大步长。双线性变换法是由（2.70）式取一次项近似得到的，所以对应的差分方程虽具有一定的精度，但当T取得太大时精度会降低。因此在选取T时，主要应考虑仿真精度和仿真速度的需求。经验表明，若取(为被仿真系统的自然频率)，可以保证较好的仿真精度。双线性变换法得到的差分方程的每项系数，都是预先一次计算确定的，在仿真递推过程中不需要重新计算，因此仿真递推的计算量较小。对于n阶连续系统，每步计算量约为(2n+1)次乘法。所以，该法适合于快速仿真。上一页下一页返回双线性变换法中由G(s)到G(z)的转换，以及由G(z)到差分方程的转换，都可很容易地编程实现。双线性变换法具有串联性。如图2.29所示，G1(s)和G2(s)相串联。若令

G(s)=G1(s).

G2(s),则有

G(z)=G1(z).

G2(z)图2.29式中G(z),G1(z),G2(z)分别是G(s),G1(s),G2(s)的双线性变换结果。这就是说，双线性变换没有破坏拉氏变换的串联性，可以用简单的低阶环节组成复杂的高阶系统。上一页下一页返回双线性变换公式不仅可以方便地应用于传递函数，也可以应用于状态空间模型。双线性变换后系统的阶次不变，且分子、分母多项式具有相同的阶次。将（2.71）式代入（2.82）式，有

（2.85）可见，G(z)的分子、分母多项式均为n次，在分子上增添了n–m个z=-1的零点。由此得到的差分方程是多步关系式（当n>1时）。由于被仿真系统一般仅给出在初始时刻处上一页下一页返回的系统的初始值，而在以前各采样点处的信息需要经过变换求得，这会给双线性变换法的应用带来一定的麻烦。可以证明，双线性变换后得到的脉冲传递函数频率特性与原连续系统的频率特性在低频段相差较小，在高频段差别较大。这一点也可以从表2.5看出。所以，双线性变换法主要用于有限带宽的系统。

上一页下一页返回5．双线性变换法仿真实例利用MATLAB的控制系统工具箱中的c2dm函数，很容易实现连续系统传递函数到离散系统脉冲传递函数之间的转换。【例2.10】考虑如图2.30所示的单位反馈系统，试用双线性变换法对其进行仿真(r(t)=1(t))。上一页下一页返回图2.30

【解】采用如图2.31所示的参数化Simulink模型，其中DiscreteFilter模块中参数Numerator设置为：nd1，Denominator为：dd1，Sampletime为：0.05。其余两个DiscreteFilter模块中的参数分别为nd2，dd2；nd3，dd3；Sampletime均为0.05。在MATLAB命令窗口运行下列指令>>T=0.05;>>n1=[13.33];d1=[124];n2=[10];d2=[10];n3=[1];d3=[11016];>>[nd1,dd1]=c2dm(n1,d1,T,'tustin');[nd2,dd2]=c2dm(n2,d2,T,'tustin');

>>[nd3,dd3]=c2dm(n3,d3,T,'tustin');上一页下一页返回图2.31(exam2_10.mdl)选择定步长的discrete（nocontinuousstates）算法，步长为0.05，启动仿真，可以得到如图2.32所示的仿真结果。上一页下一页返回图2.322.4.4根匹配法1．根匹配法原理

2．根匹配法步骤

3．根匹配法特性分析

上一页下一页返回1．根匹配法原理根匹配法的基本思想是要使离散化模型的瞬态特性和稳态特性与原连续系统保持一致。更明确地说，就是要使离散化后所得脉冲传递函数的零点和极点与原连续系统传递函数的零点和极点相匹配。所以，这种算法又称为零极点匹配法。先根据s域和z域之间的准确映射关系,直接由G(s)的零点和极点确定对应G(z)的零点和极点；再根据一定的原则确定G(z)的增益，则可以很方便地由G(s)求出与之对应的G(z)，使二者的瞬态特性和稳态特性保持一致。显然，根匹配法适合于以零极点形式描述的G(s)。上一页下一页返回2．根匹配法步骤求出连续系统的零点和极点，即将系统的传递函数写成下列的零极点形式

（2.87）利用关系在z平面上一一对应地确定零极点的位置。对于（2.87）式表示的系统，与z1,z2,…zm相对应的零点为eTz1,eTz2,…,eTzm;与p1,p2,…pm相对应的极点为eTp1,eTp2,…,eTpm

。于是，有

（2.88）

上一页下一页返回上