控制系统计算机仿真 第05讲

第二章 连续系统数字仿真的基本算法授课人 李会军 内容回顾 连续系统离散化方法 化导数为差商的方法 2 欧拉法 改进欧拉法 梯形法 思考 h代表什么意思 和截断误差与舍入误差的关系是什么 内容回顾 连续系统离散化方法 3 RK1法 RK2法 RK3法 RK4法 2 2 数值积分法的基本分析 计算稳定性 什么是计算稳定性 在仿真计算时 由于截断误差和舍入误差的存在 会使仿真误差不断积累 如果仿真误差能够被抑制 不会随计算时间的增加而无限增大 则认为计算方法是计算稳定的 否则是计算不稳定的 示例1 已知微分方程及初值如下所示 试比较在取不同步长时 其精确解与数值解之间的差异 4 思考 什么是截断误差 什么是舍入误差 2 2 数值积分法的基本分析 计算稳定性 微分方程的解析解为 精确解的曲线形式 5 注意 仿真程序见Exercise1 m 2 2 数值积分法的基本分析 计算稳定性 取步长h 0 1 h 0 075 h 0 05 分别使用欧拉法和RK4法 计算t 1 5处的y t 值 得到的结果曲线如下所示 欧拉法h 0 1 RK4法h 0 1 6 思考 当h 0 1 h 0 075 h 0 05时仿真步数分别是多少 2 2 数值积分法的基本分析 计算稳定性 欧拉法h 0 075 RK4法h 0 075 欧拉法h 0 05 RK4法h 0 05 7 2 2 数值积分法的基本分析 计算稳定性 仿真结果分析 精确解单调下降并迅速收敛到0 当h 0 1时 欧拉法和RK4法的解曲线均发散 仿真结果是错误的 当h 0 075时 欧拉法的解曲线仍然发散 仿真结果是错误的 RK4法的解曲线单调下降并收敛到0 仿真结果是正确的 当h 0 05时 欧拉法和RK4法的解均收敛到0 虽然欧拉法的解曲线是振荡收敛的 如果只要求得到t 1 5处的y t 值 则两种数值积分算法的解都可以认为是正确的 8 2 2 数值积分法的基本分析 计算稳定性 为什么会出现这种情况呢 因为数值积分算法只是一种近似计算方法 在递推计算中会引进误差 如果误差累积越来越大 将使计算不稳定 从而得到错误结果 系统稳定性与计算稳定性思考 什么是系统稳定性 如何判断系统稳定性 系统的稳定性与计算稳定性是两个不同概念 前者根据系统的微分方程 传递函数来讨论 后者根据逼近微分方程的差分方程来讨论 对于同一系统 选用的数值积分算法不同 得到的差分方程也各不相同 计算稳定性也不同 9 2 2 数值积分法的基本分析 计算稳定性 如何分析数值积分算法的计算稳定性呢 通常用简单的一阶微分方程检验算法的计算稳定性 该一阶微分方程称为测试方程 TestEquation 如果为方程的特征根 根据稳定性理论 当特征根位于s平面的左半平面上时 即时 原系统稳定 10 思考 为什么用一阶微分方程作为测试方程 2 2 数值积分法的基本分析 欧拉法的计算稳定性 欧拉法的计算公式 按照欧拉法计算的并不是y t 在t k 1 h处的真值 而是包含了各种误差的近似值 如果随着递推次数的增加 近似值的误差越来越大 说明算法是计算不稳定的 如果误差始终可控 说明算法是计算稳定的 使用欧拉法计算测试函数 如下所示 为了简化讨论 假定仅仅引入初始误差 而在递推过程中没有引进其它误差 这样的误差仅由的误差决定 11 2 2 数值积分法的基本分析 欧拉法的计算稳定性 令 b a 计算误差方程 1 当 随着计算步数增加 误差越来越大 计算不稳定 2 当 随着计算步数增加 误差越来越小 计算稳定 12 2 2 数值积分法的基本分析 梯形法的计算稳定性 令 d c 计算误差方程 如果为负实数 而步长h大于零 梯形法是计算稳定的 13 如果一阶系统是稳定的 使用梯形法仿真时 一定是计算稳定的 2 2 数值积分法的基本分析 龙格 库塔法的计算稳定性 14 RK2的计算稳定性 如果RK2计算稳定 需满足 2 2 数值积分法的基本分析 龙格 库塔法的计算稳定性 15 RK4的计算稳定性 2 2 数值积分法的基本分析 龙格 库塔法的计算稳定性 16 令 h g 计算误差方程 如果RK4计算稳定 需满足 推论 如果RKn计算稳定 需满足 2 2 数值积分法的基本分析 数值积分算法的选择原则 17 数值积分算法的选择与仿真计算的精度 速度 误差积累 计算稳定性 自启动能力等因素都有密切关系 选择原则如下 精度要求 仿真精度主要受截断误差 舍入误差的影响 它们与积分算法 仿真步长 计算时间 计算机精度等有关 选择经验 低精度问题用低阶算法来处理 如果选用高阶算法 则应在保证计算稳定性的前提下 把步长取得大一些 思考 为什么采用高阶算法时 要把步长取得大一些 2 2 数值积分法的基本分析 数值积分算法的选择原则 18 计算速度要求 计算速度取决于在给定计算时间内的计算步数和每步计算所需的时间 为了加快计算速度 在数值积分算法选定后 应在保证精度的前提下 尽量选择较大的步长 以减少计算步数 计算稳定性要求 保证数值解的计算稳定性 是进行数字仿真的先决条件 从计算稳定性的角度来看 同阶的RK法优于显式Adams法 但又不如隐式Adams法 稳定 是基础 维稳 压倒一切 2 2 数值积分法的基本分析 数值积分算法的选择原则 19 自启动能力要求 单步法具有自启动能力 多步法没有自启动能力 一般简单系统仿真时 多采用单步法 步长变化能力要求 单步法在仿真过程中 步长可以在一定的范围内变化 而多步法对步长的变化有严格的要求 如果要求仿真时进行变步长运算 最好采用单步法 对于一般控制系统的仿真 采用RK4法基本可以满足要求 2 2 数值积分法的基本分析 变步长仿真策略 20 加倍 减半法 主要思想 设定最小允许误差和最大允许误差 当估计的局部截断误差大于最大允许误差时 将步长减半 并重新计算本步 当局部截断误差在最大允许误差和最小允许误差之间时 本步计算结果有效 下步步长不变 当局部截断误差小于最小允许误差时 本步计算结果有效 下步步长加倍 总结 变步长控制策略虽然增加了每一步的计算量 但它可以很好地解决仿真中计算精度与仿真速度之间的矛盾 2 2 数值积分法的基本分析 数值积分算法仿真实例 21 在MATLAB Simulink环境下 利用数值积分算法对系统进行仿真的途径有两种 对于以微分方程给出的数学模型 通常用MATLAB语言编程实现 而对于以结构框图给出的数学模型 可采用Simulink实现 用MATLAB语言编程实现仿真的主要步骤是调用ODE解函数 MATLAB提供的常用ODE解函数如下 ode45此算法被推荐为首选算法 ode23这是一个比ode45低阶的算法 ode113用于更高阶或大的标量计算 其它 ode45表示采用四阶 五阶Runge Kutta算法 它用4阶方法提供候选解 5阶方法控制误差 是一种自适应步长 变步长 的常微分方程数值解法 2 2 数值积分法的基本分析 数值积分算法仿真实例 22 ODE解函数的调用格式 t x ode45 方程函数名 tspan x0 tol 方程函数名 描述系统状态方程的M函数名称 tspan 为仿真时间范围 例如 取tspan t0 tf t0为起始计算时间 tf为终止计算时间 x0 为系统状态变量初值 应使该向量元素个数等于系统状态变量的个数 tol 指定精度 其默认值为10 3 注意 函数返回结果为t向量和x阵 由于计算中采用了步长自动控制策略 因而t向量不一定是等时间间隔 仿真结果可以用plot t x 指令绘制出来 2 2 数值积分法的基本分析 数值积分算法仿真实例 23 例1 某地区传染病的传播模型如下 式中 x1表示可能受到传染的人数 x2表示已经被传染的病人数 x3表示已治愈的人数 试用ode45编程 对其进行仿真研究 并绘制出对应的时间响应曲线 2 2 数值积分法的基本分析 数值积分算法仿真实例 24 采用M语言进行编程 程序如下 仿真程序clearx0 620 10 70 置状态变量初值tspan 0 30 置仿真时间区间 t x ode45 fun2 3 tspan x0 调用ode45求仿真解plot t x 1 k t x 2 k t x 3 k 用不同的线型绘制仿真结果曲线xlabel time 天 t0 0 tf 30 对t x轴进行标注ylabel x 人 x1 0 620 x2 0 10 x3 0 70 legend x1 x2 x3 grid 2 2 数值积分法的基本分析 数值积分算法仿真实例 25 其中 fun