ch2.1 2.2 连续系统仿真.ppt

第3章连续系统的离散化 数值积分法 用数字计算机来仿真或模拟一个连续控制系统的目的就是求解系统的数学模型 连续系统的数学模型通常是用微分方程 组 或状态空间来描述的 所以求取系统的数值解 需要把连续系统的微分方程 组 等化为近似的差分方程 这个过程称为离散化 微分方程离散化的方法有 化导数为差商的方法泰勒展开式法数值积分法 1 3 1连续系统的离散化 3 1 1化导数为差商的方法 例如 对于一阶微分方程的初值问题 在区间 a b 上取若干个离散的等距离的点 第i个点处的导数可近似表示成差商 从而把微分方程化为差分方程 h为步长 2 在点附近的另一个点处 可以用泰勒展开式逼近为 因此可以把一阶微分方程利用泰勒展开式转化为差分方程 用泰勒展开式建立差分方程时 如果p较大 则差分方程会非常复杂 p 1呢 3 其中的高阶导数项由一阶导数逐次求导得出 3 1 2泰勒展开式法 3 1 3数值积分法 把微分方程 对于方程右边的定积分 可以采用矩形公式 或梯形公式 4 求定积分 3 2常用的数值积分方法 1 欧拉法 对微分方程 对比泰勒级数 5 如果步长h足够小 则f x y 的变化也很小 可以用矩形面积 近似代替曲线积分 例 6 欧拉法的特点欧拉法是一种最简单的数值积分算法 但是其计算精度低 因此实际中很少使用 欧拉法产生较大误差的原因 欧拉法只适用于步长h很小 或者f t y 变化很小的场合 7 2 龙格 库塔法龙格 库塔法是求解常微分方程的各种数值积分方法中应用的最广泛的一种 泰勒展开式法中 取的项数越多 计算精度越高 但是计算f的高阶导数很不方便 德国数学家C Runge和M W Kutta提出间接利用泰勒级数法 用若干个时间点处f的函数值的线性组合来代替f的各阶导数项 然后按泰勒展开式确定其中的系数 8 对于泰勒展开式 由于 所以 为了避免计算f的导数项 采用待定系数法 令 其中 W和a c是待定的系数 k为不同时间点处导数f的值 r为使用k值的个数 即计算导数f的次数 9 将上式和泰勒展开式 3 1 逐项对比确定待定系数 当r 1时 只有一项 即 将及以上项都略去 得 当r 2时 10 11 与式 3 1 比较 可得系数 上式中3个方程 4个未知数 因此解不是唯一的 如果取 则 因此 得到二阶龙格 库塔法的一般迭代公式 若取到y的三阶或四阶导数项 则有相应的三阶或四阶龙格 库塔法 四阶龙格 库塔迭代公式 12 每一步需要计算4次f的值 截断误差 龙格 库塔法的特点 从理论上说 取的项数越多 计算精度越高 但相应的计算公式更复杂 计算工作量也更大 在计算时只用到 而不直接用 各项的值 即在本次计算中 仅仅用到前一次的计算结果 而不需要利用更前面各步的结果 这类计算方法称为单步法 精度较高 对于常用的4阶龙格 库塔法 其截断误差为 13 3 3关于数值积分法的讨论 单步法和多步法单步法如4阶龙格 库塔法 它有以下特点 1 需要存储的数据量少 占用的存储空间少 2 只需要知道初值 就可启动递推公式进行运算 具有这种能力的计算方法称为自启动的计算方法 3 容易实现变步长运算 4 缺点是随着精度的提高 计算工作量会增加 与单步法相对应的还有一类数值积分方法 在它的数值积分公式中 本次计算不仅利用前一次的计算结果 还必须利用更前面时刻的值 此类方法称为多步法 如四阶亚当斯 Adams 积分公式 在相同的条件下多步法比单步法要快 但是 多步法不是自启动的 而且多步法不容易实现变步长运算 14 显式法和隐式法如果数值积分算法在计算时所用到数据均已求出 则称其为显式算法 例如 欧拉法和龙格 库塔法 如果方程右端含有未知量 则称为隐式算法 例如 Adams隐式法 尽管隐式算法计算过程复杂 计算速度慢 但有时基于对精度 稳定性等考虑 仍然经常使用 15 1 计算稳定性的概念 连续系统的数字仿真 实质上就是将给定的微分方程 组 变换成差分方程 组 然后从初值开始 逐步进行递推计算 例 已知微分方程及初值 试比较在取不同步长时 其精确解与数值积分算法解之间的差异 16 解 该初值问题的解析解为 对应的结果曲线如图所示 图1解析解 17 取 分别用欧拉法和RK4法计算 处的 所得结果曲线如图2 b g 所示 图2 b 欧拉法 h 0 1 图2 c RK4法 h 0 1 18 图2 d 欧拉法 h 0 075 图2 e RK4法 h 0 075 图2 f 欧拉法 h 0 05 图2 g RK4法 h 0 05 19 从图中可以看出 解析解单调下降并迅速收敛到0 当时 欧拉法和RK4法的解曲线均发散 数值积分算法的解是错误的 当时 欧拉法的解曲线仍然发散 对应的解是错误的 RK4法的解曲线单调下降并收敛到0 对应的解是正确的 当时 欧拉法和RK4法的解均收敛到0 虽然欧拉法的解曲线是振荡收敛的 如果只要求得到处的的解 则两种数值积分算法的解都可以认为是正确的 事实上 此时有 20 欧拉法 RK4法 解析解 为什么会出现上述情况呢 这是因为数值积分算法只是一种近似积分方法 在反复的递推计算中会引进误差 这种误差通常是由初始数据的误差及计算过程中的舍入误差产生 如果误差的累积越来越大 将使计算出现不稳定 从而得到错误的结果 系统的稳定性与计算稳定性是两个不同的概念 前者用原系统的微分方程 传递函数来讨论 后者用逼近微分方程的差分方程来讨论 由于选用的数值积分算法不同 即使对于同一系统 差分方程也各不相同 计算稳定性也就各不一样 21 通常用一个简单的一阶微分方程来考查数值积分算法的计算稳定性 微分方程及初值问题 称为测试方程 TestEquation 其中 为方程的特征根 根据稳定性理论 当特征根位于左半s平面 即时 原系统稳定 22 欧拉法的计算稳定性 考虑微分方程 使用欧拉法进行计算 有 为了简化讨论 假定仅仅引入初始误差 而在递推过程中没有引进任何其他误差 所以 a b b a 得误差方程 因此有 讨论 表明若在递推计算中引入了误差 则随着计算步数的增加 误差将逐渐扩大 表明若在递推计算中引入了误差 则随着计算步数的增加 误差会被抑制 结论 对于欧拉法 合理选择步长使 1 h 1 是保持其计算稳定性的充要条件 当 1 h 1时 当 1 h 1时 23 梯形法的计算稳定性对梯形计算公式 24 25 若原系统稳定 根据稳定性理论 应是小于零的实数 h取任何正数都可以保证下式成立 所以梯形法的计算公式在任何步长下都是稳定的 这是一种绝对稳定的计算方法 龙格 库塔法的计算稳定性1二阶龙格 库塔法由二阶龙格 库塔法计算公式得 26 2四阶龙格 库塔法按四阶龙格 库塔法计算公式得 27 yn 1 yn h K1 2K2 2K3 K4 6 28 对于其它阶数的龙格 库塔法 可以得到误差方程 令 因此稳定的条件是 29 在分析数值积分算法的精度时 通常用泰勒级数作为工具 假设前一步求得结果是准确的 即有 则用泰勒级数求得在 处的解析解为 不同的数值积分算法相当于在上式中取了不同项数之和而得到的近似解 截断误差与舍入误差 30 这种由数值积分算法单独一步引进的附加误差称为局部截断误差 它是数值积分算法给出的解与微分方程的解析解之间的差 故又称为局部离散误差 步长h越小 局部截断误差就越小 若数值积分算法的局部截断误差为 则方法为r阶的 算法的阶数可以作为衡量算法精确度的一个重要标志 31 以上的分析是在假设前一步所得结果是准确的前提下得出的 即 成立时 有 但在求解过程中 实际上只有当k 0时 才成立 而当k 1 2 3 时 并不成立 因而r阶算法求得的解的实际误差要大一些的 设yk是在无舍入误差情况下由r阶方法算出的微分方程式的近似解 y t 为微分方程的解析解 为算法的整体截断误差 则可以证明 即整体截断误差比局部截断误差低一阶 32 在计算机逐次计算时 初始数据的误差及计算过程的误差都会使误差不断积累 如果这种仿真误差能够被抑制 不会随计算时间而无限增大 则可以认为相应的计算方法是稳定的 否则是数值不稳定的 仿真误差与数值计算方法 计算机的精度以及计算步长的选择有关 仿真误差一般有如下两种 1 各种数值积分方法的差分方程都是对原微分方程的近似逼近 存在明显的截断误差 2 计算机的字长有限 计算只能限制在有限位 将引入舍入误差 若步距取得较大 截断误差就会相应地增大 反之 若步距取得较小 截断误差就会减小 但在给定时间范围内 计算次数必然增加 使舍入误差积累 相应地增加 截断误差与舍入误差 33 从图中可知 两种误差对步距的要求是矛盾的 但两者之和有一个最小值 步距最好能选在最小值 然而 实际要做到这一点是很困难的 一般只能根据经验确定一个h0附近的合理步长区 如可将h限制在系统的最小时间常数数量级上 截断误差与舍入误差 34 仿真的总误差与步长的关系不是单调函数关系 而是一个具有极值的函数关系 34 例 考虑如下二阶系统 在上的数字仿真 并将不同步长下仿真结果与解析解进行精度比较 解 采用RK4法进行计算 步长选取为0 001s到1s之间变化 并且将计算结果与解析解进行比较 求出最大误差数据列于表1中 表1 35 从表1中可以看出 当步长h 0 05时 总误差最小 当步长h小于0 05时 由于舍入误差变大而使总误差增加 当步长h大于0 05时 则由于截断误差的增加也使得总误差加大 另外 当系统的解变化激烈时 如R 0 误差对步长的变化较为敏感 当系统的解变化平稳时 步长的变化对误差的影响就要缓和得多 数值积分算法确定以后 在选择步长时 需要综合考虑 对于RK4法 步长的选择应满足 36 对于一般工程系统的分析和设计而言 仿真结果的误差不超过0 5 就可以满足精度要求了 经验表明 对于RK4法 通常可以根据系统方程中的最小时间常数来选择步长 一般取 或者根据系统中反应最快的小闭环的开环剪切频率来选择 取 上述两种选择步长的经验公式都是针对RK4法而言的 如果想要用RK2法进行仿真 则为了达到相同的精度 步长一般应取为RK4法的1 10 37 如果采用欧拉法 则为了达到相同精度 步长还要取得更小 显然RK4法的计算效率较高 这就是控制系统数字仿真中通常选用RK4法的重要原因之一 38 数值积分算法的选择原则 数值积分算法与仿真计算的精度 速度 误差积累 计算稳定性 自启动能力等因素都有密切关系 因此选择合适的数值积分算法非常重要 目前还没有普遍适用的规律 下面仅给出一些参考的原则 1 精度要求数值积分算法的仿真精度主要受截断误差 舍入误差和误差累积的影响 它们与积分算法 步长 计算时间及所用的计算机精度等有关 降低计算精度 经验表明 低精度问题用低阶算法来处理 如果选用高阶算法 则应在保证计算稳定性的前提下 把步长取得大一些 39 2 计算速度计算速度取决于在给定积分时间内的计算步数和每步计算所需的时间 为了加快计算速度 在积分算法选定后 应在保证精度的前提下 尽量选择较大的步长 以减少计算步数 40 3 计算稳定性保证数值解的计算稳定性 是进行数字仿真的先决条件 从计算稳定性的角度来看 同阶的RK法优于显式Adams法 但又不如隐式Adams法 3 自启动能力单步法具有自启动能力 多步法没有自启动能力 必须借助于单步法启动后才能开始工作 因此实现起来要困难一些 一般简单的仿真程序多采用单步法 4 步长变化能力单步法在整个仿真过程中 步长可以在一定的范围内变化 而多步法对步长的变化有严格的要求 如果要求仿真时步长可变 最好采用单步法 对于一般控制系统的仿真 通常都采用4阶龙格 库塔法 41 41 变步长仿真策略 加倍 减半 法 主要思想 设定最小允许误差和最大允许误差 当估计的局部截断误差大于最大允许误差时 将步长减半 并重新计算本步 当误差在最大允许误差和最小允许误差之间时 本步计算结果有效 下步步长不变 当误差小于最小允许误差时 本步计算结果有效 下步步长加倍 变步长控制策略虽然增加了每一步的计算量 但从总体上考虑往往是合算的 因为它可以很好地解决仿真中计算精度与计算量之间的矛盾 42 病态方程及处理 对于状态方程 状态矩阵A的特征值的实部对应于振动振幅的增减 虚部对应于振动的角频率 令 为时间常数 系统最大的时间常数决定了系统总的过渡过程时间 即决定了积分求解总的时间 当最大时间常数与最小时间常数相差悬殊时 称这样的系统或微分方程为 病态的 在仿真算法中 积分步长受限于最小时间常数 于是当病态很严重时 积分步数往往很大 造成计算时间很长 这是病态方程带来的计算困难问题 对病态微分方程往往采用变步长 变阶 隐式的仿真算法 43 3 4数值积分算法仿真实例 在MATLAB Simulink环境下 利用数值积分算法对系统进行仿真的途径有两种 对于以微分方程给出的数学模型 通常用MATLAB语言编程实现 而对于控制系统分析和设计中最常见的以系统结构图描述的数学模型 则采用Simulink实现更为方便 44 1 数值积分算法的编程实现 用MATLAB语言编程实现仿真的主要步骤是调用MATLAB中的ODE OrdinaryDifferentialEquation 解函数 MATLAB提供的常用ODE解函数如下 ode45此算法被推荐为首选算法 ode23这是一个比ode45低阶的算法 ode113用于更高阶或大的标量计算 ode23t用于解决难度适中的问题 ode23s用于解决难度较大的问题 ode15s与ode23相同 但要求的精度更高 ode23tb用于解决难度较大的问题 45 这些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 指令绘制出来 46 例 某地区某病菌传染的系统动力学模型为 式中 x1表示可能受到传染的人数 x2表示已经被传染的病人数 x3表示已治愈的人数 试用ode45编程 对其进行仿真研究 并绘制出对应的时间响应曲线 47 解 采用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 48 其中 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所示结果曲线 49 图2 4 50 2 Simulink模型的构建与仿真 对于以系统结构图描述的数学模型 采用Simulink构模并仿真是十分方便的 在大多数情况下 用户甚至无需编制一条仿真程序就可以进行仿真 Simulink模型的构建 仿真参数的设置和ODE算法的选择 Simulink仿真结果输出 与M函数的组合仿真 51 Simulink模型的构建 为了能用Simulink对系统进行仿真 首先要在Simulink环境下打开一个空白模型窗口 然后依据系统结构图给定的环节和信号流程 从Simulink模块库的各个子库中选择相应的模块 并用鼠标左键将它们拖入模型窗口 双击选择的模块 设置需要的参数 对各模块进行连接 这就构成了需要的Simulink模型 即仿真结构图 52 53 仿真参数的设置和ODE算法的选择 在模型窗口中 见下页图 的Simulation菜单下选择其中的仿真参数子菜单 SimulationParameters 就会弹出一个仿真参数对话框 如图2 5所示 图2 5 54 在仿真参数对话框中有5个标签 默认的标签为微分方程求解程序Solver的设置 在该标签下的对话框主要接受微分方程求解算法及仿真参数的设置 算法的选择通过该对话框可以由Solveroptions栏目选择不同的求解算法 定步长 Fixed step 下支持的算法 变步长 Variable step 下支持的算法如图2 6所示 一般情况下 连续系统仿真应该选择ode45变步长算法或者定步长的ode4 即RK4 算法 而离散系统一般默认地选择定步长的discrete nocontinuousstates 算法 55 图2 6 56 计算步长的选择定步长算法的计算步长应该由Fixedstepsize框填入参数进行选择 一般可以选择auto 而变步长算法建议步长范围使用auto选项 仿真时间的设置在仿真参数对话框的相关对话框中还可以修改仿真的初始时间和终止时间 57 Simulink仿真结果输出 在完成了仿真参数的设置和ODE算法的选择后 就可以启动仿真 Simulink会自动将系统结构图转换成状态空间模型并调用所选择的算法进行计算 为了得到所需要的仿真结果 除了可以直接采用Scope模块显示仿真结果曲线外 还可以将仿真结果数据传送到MATLAB工作空间中 利用plot指令绘制相应的图线 58 与M函数的组合仿真 如果Simulink模型中存在复杂的非线性环节或复杂的逻辑运算 而在MATLAB提供的所有工具箱中都找不到相应的模块 可以自己编制一个M函数 嵌入Simulink模型中 59 例 某非线性系统如图所示 试求r t 2 1 t 时系统的动态响应 解 构建系统的Simulink模型 如下图所示 为了便于研究问题的方便起见 不采用Discontinuities子库中的Saturation模块 而选择User Defined Functions子库中的MATLABFcn模块 并将参数MATLABFunction设置为 saturation zone 60 图2 18 exam2 7 mdl 61 然后编制M函数内容如下 saturation zonefunctionfunction uo saturation zone ui ifui 1uo 1 elseifui 1uo 1 elseuo ui end启动仿真后 得到如右图所示仿真结果 62 3 5面对微分方程的数字仿真 3 5 1仿真程序的主要功能 连系统数字仿真的主要过程是 将描述原系统的各种数学模型转换成仿真模型 用选定的仿真算法进行计算来获取所要求的信息 为此 应将某一类仿真问题编制成仿真程序 一般通用程序具有以下基本功能 1 具有模型转换功能 能将输入的数学模型转换成仿真模型 2 具有良好 提示显示清晰的输入界面 能方便地输入系统结构参数和仿真参数 63 3 能提供一种或数种数值积分法进行求解 4 能提供多种可选择的输出功能 如显示 打印仿真结果 绘制图形及输出文件等 3 5 2仿真程序的结构仿真程序采用自顶向下和模块化的设计方法将程序系统分成三层 主控模块 功能模块 基本子模块 64 65 主控模块 是仿真程序的主体 由主控模块进行仿真逻辑控制 调用各工程模块完成仿真的过程控制 功能模块 由若干个独立的小功能模块组成 如输入块 运行块等 小功能模块是一个或几个基本独立单一 定义明确功能的模块 它们处于程序的最下层 完成仿真过程某一阶段的任务 输入块负责完成对系统模型的参数 初值及仿真用参数 如计算步长 打印间隔 仿真时间 输出参数个数等 的设置及修改 66 运行块负责进行运动状态的计算 反复计算微分方程右端函数 调用积分子模块后 完成各步的计算 运行块是整个仿真程序的核心 它决定了仿真的精度及速度 输出块负责输出仿真结果 按使用者要求提供各种输出格式 如输出列表文本 绘制图形曲线或存放在数据文中 无论原系统的数学模型以那种形式给出 最后都需要转化为一阶微分方程组的形式 才能用数值积分法求解 这是一个由外部模型到内部 67 模型的转化问题 数字仿真程序要在计算机上完成模型转化 建立相应的内部状态模型 并运行求解 下面给出一个数字仿真的例子 例 已知一微分方程组为 68 初始条件 x1 0 0 x2 0 0用四阶龙格 库塔法的仿真结果如下 69 70 3 6面向结构框图的数字仿真程序 面向微分方程的数字仿真过程 存在以下不足之处 完整的控制系统 常常是以框图的形式来描述的 这样 将各环节分开 再改为状态方程或传递函数形式是相当麻烦的 如果要研究参数变化的影响 就显得非常的不方便在实际系统中 常常有子闭环 若系统用状态空间或微分方程描述 则子闭环的影响就难以描述系统中往往存在非线性环节 上述方法很难处理 71 在控制系统中 我们看到的往往是系统的框图 而框图是由各个典型环节连接而成的 于是可利用这些典型环节及其连接方式 进行仿真 它有以下几个特点 便于研究各环节参数对系统的影响可以得到各个环节的动态响应可对多变量系统进行仿真 3 6 1典型仿真模型1 典型环节的确定一个控制系统的结构图通常是由一些不同的环节组成 一般常见的环节 有积分环节惯性环节振荡环节 比例惯性 超前 滞后 环节 比例积分环节 72 为了编程方便 各常用环节可用一个典型环节来表示 选定的典型环节 首先要求它有典型性 即能用它来描述较多环节类型 保证所编制的环节有较大的通用性 其次要求环节的结构简单 用它编制的程序是简短易于计算的 常见的环节中积分环节的结构最为简单 但在构成系统时的组件多 增加了系统结构的复杂性 选比例 惯性环节作为典型环节 能够做到不增加结构的复杂性 所编制的程序具有简单通用性强的优点 73 比例惯性环节的传递函数形式为利用这个典型环节 只要改变系数便可组成其它常用环节 积分环节对应于比例积分环节对应于惯性环节对应于比例惯性环节对应于 i 1 2 n 74 对于振荡环节 如果其特征方程可以因式分解 实根时 用两个串联的惯性环节替代 共轭复根时 用两个串联的惯性环节加负反馈来实现等效环节如下图所示 一般复根时则化成积分环节和惯性环节加负反馈的结构图见下页图 如 其特征根为和 75 图中 2 传递函数的矩阵表达式 用n个环节描述n阶控制系统的传达函数的表达式为 76 将其转换成一阶微分方程组形式为对每个环节列写出对应的传递函数为 或写成一阶微分方程组形式 77 写成方程形式式中 78 A B C D都是维的系数矩阵 u和y分别为输入输出的单列向量 3 连接矩阵 表示了各环节输入与输出之间的数量关系 但是一个环节的输入可能由几个环节的输出构成 因此要根据具体的情况找出系统各环节的连接关系即环节的输入量是由哪些环节的输出量组合而成的 79 下图表示一个控制系统的结构图 该系统由五个环节组成 系统中的比例系数已知 各环节的输入用表示 输