y (0) k 1
yk
hf
(tk ,
yk )
yk
1
yk
1 2
h[
f
(t
k
,
y
k
)
f (tk 1 , y (0) k 1 )]
龙格-库塔法的基本原理
在连续系统仿真中，主要的数值计算
工作是求解一阶微分方程：dy
dt
f
(t,
y)
若已知y的初值y0,再按离散化原理，对上
式我们可以写成：y(t ) k1
y(tk)
tk1
tk
f
(t,
y)dt
再对上式的右端函数f(t,y)(为任意非线性函 数)在tk附近展开成泰勒级数，依照展开的 阶次不同我们就构成了不同的龙格-库塔公 式。
二阶龙格—库塔公式,记在tk时刻的状态变量为yk,并 定义两个附加向量型变量 :
yk1 k1 f
yk
h 2
(k1
(tk , yk )
2. 仿真程序流程框图
开始
输入系统阶次、计算步长、 阶跃函数幅值
输入传递函数分子、分母系数
求状态方程系数矩阵 A,B,C
求四阶龙格库塔法各系数 Ki，j
计算状态变量
计算输出值
N
时间到？
Y 输出仿真结果
结束
MATLAB提供了两个常微分方程求解的函数ode23()、ode45(),分别采用了二阶三级的 RKF方法和四阶五级的RKF法，并采用自适应变步长的求解方法，即当解的变化较慢 时采用较大的计算步长，从而使得计算速度很快；当解的变化较快时，步长会自动变 小，从而使计算精度提高。
允许
误差
E k
步长控制
改变下一步 计算步长
误差估计
积分算法
本步 计算
2、误差估计
• 通常采用的方法是设法找到一个比目前使用的 龙格-库塔公式低一阶的R-K公式，将两式计算
结果之差视为误差估计值。
• 例如：现以RKM4-5为计算公式
y k 1
y k
h 6
(k1
4k 4
k
)
5
k1
f
(tk
,
y) k
6.1.2 离散化原理
在数字计算机上对连续系统进行仿真时，首先遇到的问题是，数字 计算机的数值及时间都是离散的（计算精度，指令执行时间）， 而被仿真系统的数值和时间是连续的，后者如何用前者来实现？
设系统模型为：y f ( y,u,t) ,其中u(t)为输入变量，y(t)为系统状态变 量 t似=k。原h。令理如仿。果真时u’(间tk)间≈隔u(为tk)h, y,离’(t散k) 化≈ y后(t的k)，输则入认变为量两为模u’型(tk等), 其价中，t称k表为示相
6.2.1 仿真过程的三类误差 通常，系统仿真的最终精度与现场原始数据的采集、使用的 计算机设备档次、仿真计算时的数值积分公式等均有相应的 关系，可以分为以下3种情况。 1. 初始误差
在对系统仿真时，要采集现场的原始数据，而计算时要提供 初始条件，这样由于数据的采集不一定很准，会造成仿真过 程中产生一定的误差，此类误差称为初始误差。 要消除或减 小初始误差，就应对现场数据进行准确的检测，也可以多次 采集，以其平均值作为参考初始数据。
控制系统的数学模型经过合理的近似及简化，大多数建立为常微分方程 的表达形式。由于数学计算的难度和实际系统的复杂程度，在实际中遇到 的大部分微分方程难以得到其解析解，通常都是通过数字计算机采用数值 计算的方法来求取其数值解。在高级仿真软件（例如MATLAB）环境下， 已提供了功能十分强大、且能保证相应精度的数值求解的功能函数或程序 段，使用者仅需要按规定的语言规格调用即可，而无需从数值算法的底层 考虑其编程实现过程。
第6章 数值积分法仿真
本章主要教学内容
本章主要介绍控制系统数学模型的相关知识，通 过本章的学习，读者应掌握以下内容：
➢求解常微分方程数值解的一般方法 ➢数值积分法的基本概念及其常用方法 ➢以系统微分方程或传递函数作为数学模型的仿真过程及程序 设计方法 ➢以系统动态结构图作为数学模型的仿真过程及程序设计方法 ➢仿真步长的选择与系统仿真精度和稳定性的对应关系 ➢快速仿真算法的概念、特点及其应用
6.2.2 稳定性分析
由于系统仿真时存在误差，对仿真结果产生会影响。若计 算结果对系统仿真的计算误差反应不敏感，那么称之为算法稳 定，否则称算法不稳定。对于不稳定的算法，误差会不断积累， 最终可能导致仿真计算达不到系统要求而失败。
1. 系统的稳定性与仿真步长的关系
一个数值解是否稳定，取决于该系统微分方程的特征根是 否满足稳定性要求，而不同的数值积分公式具有不同的稳定区 域，在仿真时要保证稳定就要合理选择仿真步长，使微分方程 的解处于稳定区域之中。
k2)
k2
f
(tk
h,
yk
hk1 )
四阶龙格—库塔公式 ：
y
k
1
yk
h 6
(k1
2k2
2k3
k4 )
k1 f (tk , yk )
k 2
f (tk
h 2
,
yk
h 2
k1
)
k
3
f (tk
h 2
,
yk
h 2
k2 )
k4 f (tk h, yk hk3 )
不论几阶RK法，它们的计算公式都是由两部分组成，即上一步的结果yk和 步长h乖以tk至tk+1时间间隔间各外推点的导数ki的加权平均和
2. 舍入误差
目前，系统仿真大都采用计算机程序处理和数 值计算，由于计算机的字长有限，不同档次的计 算机其计算结果的有效值不一致，导致仿真过程 出现舍入误差。 一般情况下，要降低舍入误差应 选择挡次高些的计算机，其字长越长，仿真数值 结果尾数的舍入误差就越小。
3. 截断误差
当仿真步距确定后，采用的数值积分公式的阶 次将导致系统仿真时产生截断误差，阶次越高， 截断误差越小。通常仿真时多采用四阶龙格—库 塔法，其原因就是这种计算公式的截断误差较小。
6.1 数值积分法
6.1.1 概述
数字仿真模型、算法及仿真工具
控制系统的数字仿真是利用数字计算机作为仿真工具，采用数学上的各 种数值算法求解控制系统运动的微分方程，得到被控物理量的运动规律。 通常，计算机模拟被控对象是用一定的仿真算法来实现被控对象的运动规 律，这是基于被控对象的数学模型来完成的。
微分方程数值解的方法主要是数值积分法。 对于形如 y f (y,u,t) 的系统，已知系统状态变量y的初值y0,现
要计算y随时间变化的过程y(t)，对微分方程的积分可以写作： y t
y(t) f (t, y)dt t 0
0
右图所示曲线下的面积就 是y(t),由于难以得到积分 的数值表达式，所以采用 近似的方法，常用有三种形式：
3、最优步长法
• 基本方法是根据本步的误差估计来近似地确定下一步可能的最大
步1)给长定，允步许骤的如相下对: 误差ε0，设本步步长为hk,本步相对估计误差为ek，
e E y ek由下式求得:
/(| | 1)
k
k
k
E h 若采用的RK公式是m阶，则上式中的Ek可以表示为
( ) m
k
k
e t h y 通常取ζ=tk,则有:
在仿真计算中，h值不宜选的太小，因为这样会加大计算量； 也不能选的过大，否则会导致仿真系统不稳定或误差积累过大。
通常掌握的原则是： 在保证计算稳定性及计算精度的要求下，尽可能选较大的 仿真步长。对于一般工程系统的仿真处理，采用四阶龙格—库 塔法居多 。
龙格-库塔法的步长控制策略
• 控制策略由误差估计和步长控制两方面 组成，其基本思想如下图所示:
Syntax
[T,Y] = solver(odefun,tspan,y0)
[T,Y] = solver(odefun,tspan,y0,options)
where solver is one of ode45, ode23, ode113, ode15s, ode23s, ode23t, or ode23tb. ArgumentsodefunA function that evaluates the right-hand side of the differential equations. All solvers solve systems of equations in the form or problems that involve a mass matrix, . The ode23s solver can solve only equations with constant mass matrices. ode15s and ode23t can solve problems with a mass matrix that is singular, i.e., differential-algebraic equations (DAEs). tspanA vector specifying the interval of integration, [t0,tf]. To obtain solutions at specific times (all increasing or all decreasing), use tspan = [t0,t1,...,tf].y0A vector of initial conditions.optionsOptional integration argument created using the odeset function. See odeset for details.p1,p2...Optional parameters that the solver passes to odefun and all the functions specified in options