第6章系统时间响应仿真研究报告

本章将给出动力学系统仿真算法的设计思想和分析方法，并介绍由这些思想得到的一些常用仿真算法。根据实际问题的需要，灵活应用本章给出的常用算法的构造思想，将它们适当的组合，可构造出合适的算法，实现对复杂动力学系统有效的数字仿真。

在工程领域中，连续系统是最常见的系统，其仿真方法是系统仿真技术中最基本、最常用和最成熟的。进行数字仿真首先要建立被仿真系统的数学模型，并将此模型转换成计算机可接受的、与原模型等价的仿真模型，然后编制仿真程序，使模型在计算机上运转。计算机仿真1/9/20251第6章系统时间响应及其仿真仿真算法系统仿真MATLAB的函数采样控制系统仿真1/9/20252引言：对象与工具的矛盾连续系统数字计算机被仿真系统的数值及时间均具有连续性数字计算机的数值及时间均具有离散性对象与工具的矛盾前者如何用后者来实现？如何保证离散模型的计算结果从原理上的确能代表原系统的行为，这是连续系统数字仿真首先必须解决的问题。相似原理

如何将连续系统的数字模型转换成计算机可接受的等价仿真模型，采用何种方法在计算机上求解此模型，这是连续系统数字仿真算法要解决的问题。1/9/20253数字仿真的本质和基本要求用数字仿真的方法对微分方程的数值积分是通过某种数值计算方法来实现的。任何一种计算方法都只能是原积分的一种近似。因此，连续系统仿真，从本质上来说，是从时间、数值两个方面对原系统进行离散化，并选择合适的数值计算方法来近似积分运算，由此得到离散模型来近似原连续模型。相似原理用于仿真时，对仿真建模方法有三个基本要求：(1)稳定性(2)准确性(3)快速性对系统进行时域仿真分析，实际上就是要求解微分方程的“初值问题”。求的解精度问题效率问题存在问题1/9/20255连续系统数字仿真的两种方法⑶实时半实物仿真原理基于离散相似原理建立的欧拉法、梯形法、Adams法基于Taylor级数匹配原理建立的Runge—Kutta法、线性多步法⑴数值积分法：⑵离散相似法：采样控制系统的仿真方法⑸数值积分法的选择与计算步长的确定

⑷数值积分法稳定性分析

离散相似法数值积分法基于离散相似原理基于泰勒级数匹配原理连续系统近似离散化1/9/202566.1仿真算法积分的几何意义：曲线下面的面积积分的含义：离散和连续和1/9/202576.1.1数值积分法的基本原理⑴已知描述某系统的一阶微分方程及其初值为

在微分方程理论中称为初值问题

方程的解为

时的连续解为

在差分方程

问题的关键：如何计算此积分？数值积分法的说法是从数学观点提出的，离散相似原理的说法揭示了本质，与工程实际更接近，两者其实是统一的。1/9/20258数值积分法的基本原理⑴数值积分法（数值解法），就是对一阶常微分方程（组）建立离散形式的数学模型——差分方程，并求出其数值解。

⑵关键是如何计算Qk！围绕Qk，产生了各种各样的数值积分法！不同的积分方法，对系统求解的精度、速度和数值稳定性等，都有不同的影响。⑶根据已知的初值y0，可逐步递推算出以后各时刻的数值yi。采用不同的递推算法，就出现了各种各样的数值积分法。⑷数值积分是解决初值已知，对f(t,y)进行近似积分，对y(t)进行数值求解的方法。所谓数值解法，就是寻求初值问题的解在一系列离散点的近似解（数值解）结论计算步长：相邻两个离散点的间距称为计算步长或步距：h=tk+1－tk

1/9/20259欧拉法很少实用，但能说明构造数值解法一般计算公式的基本思想。对微分方程积分，写作图示曲线下的面积就是y(t)。在一个步距内，有在t＞t0时，f(t,y)是未知的，上式右端的积分是求不出来的。为了求此积分，把积分间隔取得足够小，使得在tk与tk+1之间的f(t,y)可以近似看成常数，这样便得到用矩形公式计算积分得近似公式：当t=t2时，对于任意时刻，注意：f(tk,yk)也就是y(tk)的导数。一般递推差分方程形式

1/9/202510梯形法

为了提高精度，可考虑用梯形代替矩形来近似小区间的曲线积分表示的曲面面积。梯形法近似积分形式式中隐含有未知量fk+1，梯形法是隐函数形式。一般用欧拉法估计初值，用梯形法校正：欧拉法估计梯形法校正通过反复迭代，直到满足误差ε要求梯形法实质上是采用了连续两点斜率平均值，以提高计算精度。估计—校正方法1/9/202511数值积分统一公式

这一思想被广泛应用于许多算法之中。实际应用时，可采用加权平均，即在每一步中取若干点，分别求出其斜率，然后加不同的权。数值积分统一公式

f(t,y)——y(t)的导数权值权值计算步距前一步或多步计算的结果梯形法公式1/9/202512几个概念⑴⒈单步法与多步法解初值问题的各种数值方法的共同特点是：步进式，即从最初一点或几点出发，每一步根据yk一点或与前面几点yk-1，yk-2，…，来计算yk+1的值，这样逐步递推。单步法：从yk推进到yk+1只需用到时刻tk的数据时，称为单步法。多步法：从yk推进到yk+1需要用到时刻tk以及过去时刻tk-1，tk-2，…的数据时，称为多步法。⒉显式与隐式显式：计算yk+1所用到的数据均已解算出来；隐式：在算式中隐含有未知量。1/9/202513几个概念⑵⒊截断误差分析数值积分的精度常用泰勒级数作为工具。假设前一步的结果y(tk)是精确的，则在数学上，可用泰勒级数求得下一步的精确解为：表示高阶无穷小之意。若只取前两项之和来近似，则由这种方法单独一步引入的附加误差——局部截断误差（局部离散误差），是近似值与微分方程的解之间的误差。若某种方法的局部截断误差为，则称它有r阶精度，即该方法是r阶的。r是衡量精度的重要标志。欧拉法只是精确解的一次近似式，因此欧拉法的截断误差为

，欧拉法为一阶精度。1/9/202514几个概念⑶⒋舍入误差舍入误差与计算机的字长、所使用的数字系统、数的运算次序以及子程序的精确度等有关。舍入误差与计算步长成反比，步长小，计算次数就多，舍入误差就大。由于计算机技术的发展，现代计算机仿真可以忽略舍入误差。⒌数值稳定性采用数值方法解稳定的方程（稳定的系统）时，应保持系统稳定的特征，即要求用于计算的差分方程是稳定的。但是，由于计算机逐次计算时，初始数据的误差及计算过程的舍入误差对后面的计算结果将会产生影响（误差会传播）。所以带来计算数值是否稳定的问题。所谓稳定性问题就是指误差的积累是否受到控制的问题。一般地，如果计算结果对初始数据的误差以及计算过程的舍入误差不敏感的话，就说相应的计算方法是稳定的，否则称之为不稳定。1/9/202515泰勒级数匹配原理表示高阶无穷小直接应用Taylor展开式构造精度较高的方法是不现实的：计算y(t)的高阶导数很困难，求复合函数的各阶偏导数往往更复杂。即使能得到解析公式，计算量也是很可观的。如何计算y(tm+1)的近似值ym+1能否不计算高阶导数、各阶偏导数？计算精度与精确展开式的前若干项一致？应用这种原理构造的ym+1的计算过程将不需要计算y(t)的高阶导数，或者只需计算一些比较低阶的导数，但与精确值y(tm+1)之间的误差是h的高阶无穷小量。用点(tm，y(tm))邻域中的若干个点(t，y(t))上的值来构造近似值ym+1的计算公式。要求该公式在(tm，y(tm))处的Taylor展开式与精确展开式的前若干项一致。1/9/202516龙格—库塔法基本思想利用低阶导数构成的曲线去拟合含有高阶导数的曲线，避免计算高阶导数。

1/9/202517二阶龙格—库塔法设原微分的解具有如下的形式将K2展成泰勒一次近似式代回、比较，有多个解！a1=a2二阶RK法公式如何选取待求的参数a，b1和b2？1/9/202518四阶RK法公式RK—4公式是4点导数的加权平均和，由于第2点、第3点以h/2为步长，精度较高，所以加权平均时各取2份，第1、4点各取1份。龙格—库塔法的基本思想

用f(tm,ym)在几个不同点的数值加权平均代替φ(t，y(t)，h)的值，而使截断误差的阶数尽可能高。也就是说，取不同点的斜率加权平均作为平均斜率，从而提高方法的阶数。1/9/202519龙格—库塔法的一般形式不论几阶龙格—库塔法，它们的计算公式总是由两部分组成：上一步的计算结果ym及步长h乘以tm～tm+1中各点导数的加权平均和。1/9/2025206.1.2数值积分方法的选择从前面的分析可知，数值积分方法就是采用不同的差分方程来逼近原微分方程。数值积分方法的选择与仿真的精度、速度、计算稳定性、自启动能力等密切相关，涉及因素较多，至今尚无一种确定的方法来选择最好的积分形式。这里提出一些选择时应考虑的因素。数值积分的精度主要受三个因素的影响：截断误差、舍入误差、积累误差。其中截断误差取决于积分方法的阶次与步长。在同一种算法下(阶次相同)，步长越小，截断误差越小；同样的步长下，算法阶次越高，截断误差越小。舍入误差主要取决于计算机的字长，字长越长，舍入误差越小。而积累误差对计算精度的影响，又与积分算法和数值积分计算时间的长短有关。总之，应根据精度要求，合理地选择积分方法和阶数。(1)精度要求1/9/202521数值积分方法的选择⑴计算速度取决于积分步长及每一步积分所需的时间。每一步的计算时间与积分方法有关，它主要取决于计算右函数的次数。在进行数值积分中，最费时的部分是积分变量导数的计算。RK4在每一步中都要计算4次导数，费时较多；而显式Adams法中，每步只要求计算1次变量导数，计算速度相对较快。在积分方法已定的条件下，应在保证精度要求的前提下，尽量选用较大步长，以缩短积分时间。(2)计算速度进行数字仿真，首先应保证所选方法的数值解是稳定的，否则计算结果就失去了意义。数值方法的选择具有较大的灵活性，应根据具体情况而定。(3)数值计算的稳定性1/9/202522数值积分方法的选择⑵针对一些变量随时间变化过程中包含有变化较快及较慢的部分，采用根据给定计算误差，调整积分步长的策略，使快过程计算时步长会自动缩小，慢过程计算时步长又自动放大，可以获得较好的仿真速度和精度。(5)变步长能力(4)自启动能力

单步法具有自启动能力，多步法没有自启动能力，需用单步法启动后才能计算。1/9/202523步长选择原则⑴

(6)步长选择原则

选择步长除了应考虑上述诸因素的影响外，一般可根据采用的积分方法的稳定半径rk来确定积分步长。为了保证计算的稳定性，步长h一般应限制在最小时间常数Tmin（相当于最大特征值λmax的倒数）的数量级，因此h的选择首先应满足

或

还可根据经验公式选取步长。如工程中采用RK4时，h可选为

（见采样定理）或

ωc——系统开环频率特性的剪切频率（穿越频率）；

tr——系统在阶跃函数作用下的上升时间；

ts——系统在阶跃函数作用下的调整时间。1/9/202524步长选择原则⑵系统中最小时间常数的极点只影响瞬态过程起始段的形状，而瞬态过程则主要由那些靠近虚轴的主导极点所决定，由于固定步长是按起始段选取的，这就会造成后面阶段计算量的浪费，因而可采用变步长的方法。应当指出1/9/2025256.1.3时域离散相似法标量齐次微分方程齐次状态方程由转移而来，是一种向量变换关系。变换矩阵状态转移矩阵注意：状态转移矩阵与线性变换矩阵的区别！系统自由运动的性质完全由状态转移矩阵确定！1/9/202526状态方程的解状态方程的解为：如何解得状态转移矩阵！非齐次状态方程零输入响应自由响应零状态响应强迫响应关键适当地选择控制作用u(t)，可使系统的响应过程，即状态X(t)按预定性能指标的要求变化。结论1/9/202527时域离散相似法根据状态方程的解式，可以得到

1/9/202528时域离散相似法令τ＝kT＋t，并考虑B是常数阵，而在相邻的采样点之间保持不变，即当kT≤τ＜(k＋1)T

时，，可得

离散化模型实质上就是原来连续系统的解1/9/202529连续系统的离散化保持器采样器采样周期系统开环幅值穿越频率保持器传递函数sysd=c2d(sysc,Ts,method)零阶保持器‘zoh’1/9/2025306.2系统仿真的MATLAB函数[T,Y]=ode45(‘f’,tspan,y0,option)[T,Y]=ode23(‘f’,tspan,y0,option)[T,Y]=ode113(‘f’,tspan,y0,option)[T,Y]=ode15s(‘f’,tspan,y0,option)[T,Y]=ode23s(‘f’,tspan,y0,option)[T,Y]=ode23t(‘f’,tspan,y0,option)[T,Y]=ode23tb(‘f’,tspan,y0,option)1/9/202531例6—1求方程的解1/9/202532例6—2

设表示为一阶微分方程组求方程的解1/9/202533时间响应仿真的函数step(sys,T)T=[T0,Ts,Tfinal]

Ts:采样周期step(sys)系统默认例6—4基本调用格式阶跃响应仿真函数step(sys,Tfinal)1/9/202534step(sys1,sys2,...)step(sys1,'r',sys2,'y--',sys3,'gx',...)多系统调用格式例6—5绘制两个系统系统1系统21/9/202535时间响应仿真的函数impulse(sys)Initial(sys,X0)例6—6lsim(sys,U,T)例6—7gensig(Type,Tau)脉冲响应仿真函数impulse(sys,Tfinal)impulse