第2章连续系统数值积分法的时域数字仿真 - ke.ppt

第2章连续系统数值积分法数字仿真 要对系统进行仿真处理 首先应知道系统的数学模型 然后才可以对系统进行模拟 同时 如果知道系统的模型 就可以在此基础上设计一个合适的控制器 使得系统响应达到预期的效果 1 在线性系统理论中 一般常用的数学模型有 1 微分方程模型2 传递函数模型 系统外部模型 3 状态方程模型 系统内部模型 2 微分方程模型是控制系统模型的基础 一般来说 利用机械学 电学 力学等物理规律 得到控制系统的动态方程 这些方程对于线性定常连续系统而言是一种常系数的线性微分方程 3 静态和动态数学模型1 静态模型 表达式是不含时间变量t的代数方程 平衡状态下各变量间对应关系 变量不随时间变化 2 动态模型 表达式是含有时间变量t的微分方程 描述了系统的非平衡状态 变量随时间而变化 静态模型包含在动态模型之中 4 线性系统概念1 叠加原理叠加性 在u1 t 作用下 输出y1 t u2 t 作用下 输出y2 t 则在u1 t u2 t 作用下 输出为y1 t y2 t 2 线性系统既满足叠加性又满足齐次性的系统 5 非线性系统概念1 不满足叠加原理 2 系统中有一个环节为非线性环节 整个系统既是个非线性系统 6 建立数学模型的方法1 机理分析法通过对系统内部运动机理的分析 根据系统所遵循的物理或化学规律 在忽略一些次要因素或作出一些近似处理后进而得出系统特性方程 往往表现为微分方程或代数方程的形式 称为机理模型 2 系统辨识发 假定数学模型的结构 对实际系统加入某典型信号 得到实际系统的输出数据 按照一定的原则 由输入和输出数据来确定模型参数 例1RLC串联电路 列写图示RLC串联电路的微分方程 ui t 为输入量 uo t 为输出量 解 设回路电流为i t 则由基尔霍夫 Kirchihoff 定律得 消去中间变量i 将 2 代入 1 中 得 式 3 称为线性定常系数二阶微分方程 将刚才得到的RLC微分方程进行改写 令 T1 L R T2 RC则式 3 化简得 若利用u0 q C则式 3 可变为 时间常数 例2弹簧 质量 阻尼系统 求图示弹簧 质量 阻尼机械位移系统的微分方程 外作用力F为输入 位移x为输出量 解 设质量m相对于出示平衡状态的位移为x 根据牛顿第二定律 得 整理得 问 方程中为何不出现重力 惯性力 阻尼力 弹簧阻力 7 相似系统的概念1 电 力系统具有相同的微分方程形式 2 相似量 在微分方程中占据相同位置的物理量 如 L m R f 1 C k q x ui f3 意义 利用电路或其它简单系统研究复杂系统 在自动控制理论中 数学模型有多种形式 1 时域中常用的数学模型有 微分方程 差分方程和状态方程 2 复数域中有传递函数 结构图等 3 频域中有频域特性等 2 1 1控制系统的时域数学模型控制系统的微分方程是在时间域描述系统动态性能的数学模型 在给定外作用及初始条件下 求解微分方程可以得到系统的输出相应 这种方法比较直观 借助计算机可快速准确地求得结果 2 1连续系统的数学模型 连续系统一般是由微分方程来描述 而现行系统又是以线性常微分方程来描述 设系统的输入信号为u t 输出信号为y t 则线性定常系统由下述n阶线性常微分方程来描述 2 1 1 或 微分方程主要用于表征运动系统的动态特性 所以作为和时间联系最紧密的描述方式 微分方程及其时域解是最基本的描述方式 但应用古典方法解微分方程比较复杂 对于二阶以内系统研究较清楚 但对于高阶系统 没有一般封闭的解或解析解 所以一般采用拉普拉斯变换的方法求解 2 1 2控制系统的复数域数学模型用拉普拉斯 Laplace 变换法求解线性系统的微分方程时 可以得到控制系统在复数域用的数学模型 传递函数 1 基本思想在零初始条件 y及u的各阶导数的初值为零 下 通过拉普拉斯变换 将微分方程变为S域 复数域 内的代数方程 在S域内研究对象的运动 进行系统的综合 必要时 可通过反拉普拉斯变换变为时域形式 2 优点 将微分方程问题转换为代数方程问题 传递函数不仅可以表征系统的动态性能 而且可以用来研究系统的结构或参数变化时对系统性能的影响 3 定义传递函数 零初始条件下 系统或环节输出量的拉氏变换与输入量的拉氏变换之比 则对式 2 1 2 两边取拉普拉斯变换得 2 1 3 式中 Y s L y t U s L u t 于是系统的传递函数描述形式为 2 1 3 传递函数是经典控制理论研究的主要工具微分方程和传递函数是用系统的输入和输出之间的关系来描述系统的 表示了系统的外部特征 所以称其为 外部模型 微分方程是在时域内研究 传递函数是在S域 复域 内研究 二者可以相互转换 但存在很多局限性 当系统具有多个输入和输出 以及系统参数随时间变化 多变量 时变 时 研究起来就比较困难 2 1 3控制系统的状态空间表达式微分方程和传递函数都只描述了系统输入与输出之间的关系 而没有考虑系统内部状态的动态运行 仅仅实现系统输入与输出的关系是不够的 还必须复现模型的内部变量 即状态变量的动态变化规律 系统在时间域中的行为或运动信息的集合称为状态 确定系统状态的一组独立 数目最小 变量称为状态变量 用向量x表示 它将复杂系统化为一阶线性方程组 再运用矩阵理论 借助数字计算机获取其解 线性定常系统状态空间表达式的一般形式为 式 2 1 4 是由n个一阶微分方程组成 称为状态方程 式 2 1 5 是由m个线性代数方程组成 称为输出方程 2 1 4 2 1 5 用RLC网络说明式 2 1 4 2 1 5 1 状态方程由系统的状态变量构成的一阶微分方程组称为系统的状态方程 如图 通常以uo t 和i t 作为系统的状态变量 得 即 1 式 1 即为图示RLC系统的状态方程 若将状态变量用一般符号xi表示 即令x1 uo t x2 i t 并写成适量矩阵形式 则式 1 变为 则 一般式为 2 2 输出方程在指定系统输出的情况下 该输出与状态变量间的函数关系式 称为系统的输出方程 在图示RLC系统中 指定x1 uo t 作为输出 输出一般用y表示 则有 y uo t 或y x1写成矩阵形式 则一般形式 3 3 状态空间表达式状态方程 2 和输出方程 3 总和起来 构成一个系统完成的动态描述 称为系统的状态空间表达式 即 2 1 4 2 1 5 A为n n维的系统矩阵 表示系统内部各状态变量之间的关联情况 由控制对象的参数决定 B为n r维的输入矩阵 表示输入对每个状态变量的作用情况 C为m n维的输出矩阵 表示输出与每个状态变量之间的关系 D为m r维的前馈矩阵 表示输入对输出的直接传递关系 一般情况D 0 对于单输入 单输出系统 输入量u及输出量y均为标量 B为列向量 C为行向量 D为标量 状态空间表达式是在状态空间描述系统的 反映了系统的内部特性 所以称其为 内部模型 现代控制理论分析和设计系统的主要方法 适用于多变量系统 时变系统 非线性系统的研究 1 微分方程与传递函数之间的转换通过拉普拉斯及反拉普拉斯变换即能实现二者之间的转换 2 微分方程化为状态空间表达式状态空间表达式有两种建立方法 从系统的机理出发建立状态空间表达式和从已知系统的高阶微分方程或传递函数求相应的状态空间表达式 后者也称为 实现问题 从高阶微分方程或传递函数变换为状态方程 即分解为多个一阶微分方程 那么此时的状态方程可以有无穷多种形式 这是由于在同一系统中 状态变量选取的不同 状态方程也不同 2 2模型转化 实现问题 2 2 1模型转换 不含输入量u导数项的微分方程假定一个连续系统用微分方程描述为 选取y t 及其各阶导数为状态变量 即取 2 2 1 2 2 2 则可将时 2 2 1 化为一阶微分方程组 2 2 3 这个一阶微分方程组就是以状态变量x1 x2 xn表示的系统状态方程 写为矩阵方程形式 即 而系统的输出方程y x1 写成矩阵形式 即 2 2 4 2 2 5 由式 2 2 4 和 2 2 5 构成的系统状态空间表达式可简写为 包含输入量u的导数项的微分方程若系统微分方程中包含输入量u的导数项 如式 2 1 1 所示 则为了在状态方程中不出现导数项 可以选取如下的n个变量为状态变量 2 2 6 2 2 7 式中 0 1 2 n 1是n个待定常数 对xn求导数 并考虑 2 1 2 有 由 2 2 7 将y n 1 y y均以xi及u的各阶导数表示 整理得 令上式中u的各阶导数项的系数均为零 可确定各 的值 式中 1 2 n可根据式 2 1 1 的系数求得 即 2 2 8 根据所选状态变量得出一阶微分方程组 2 2 9 以及输出方程 2 2 10 写成矩阵方程形式 2 2 11 2 2 12 由式 2 2 11 和 2 2 12 构成的系统状态空间表达式可简写为 2 2 13 习题 将下列微分方程转换为状态空间表达式形式 例1 给定微分方程形式为 试将其转化为状态空间表达式形式 3 传递函数化为状态空间表达式对于式 2 1 3 描述的系统 引入一个中间变量x t 其拉普拉斯变换为X s 若x t 及其各阶导数的初值均为零 可将式 2 1 3 改写为 2 2 16 2 2 17 将式 2 2 14 和式 2 2 15 取反拉普拉斯变换 可得 2 2 14 2 2 15 令 比较式 2 2 1 和式 2 2 17 可以看出 两个方程是一样的 故可选状态变量 则可得可控标准型状态空间表达式为 式中 以上方法得出的是可控标准型 应用不同的方法还可以推导出可观标准型 约当标准型等形式 对系统进行分析 由状态方程 将状态变量代入上式 同时对输出方程两边取拉普拉斯变换得 在零初值条件下 两边取拉普拉斯变换 4 状态空间表达式化为传递函数 2 2 2MATLAB应用 1 模型在MATLAB中的表示方法1 传递函数的表示连续系统 tf num den 离散系统 tf num den T T 采样时间num 传递函数分子系数构成的向量 den 传递函数分母系数构成的向量 num b1 b2 bm den a1 a2 an 分子分母向量内容分别是原传递函数的分子分母系数的降幂排列 例 传递函数为 可以对传递函数进行首一化处理 形成首一式 即使a1 1 对于如有 s 2 s 3 乘积形式的传函 可以利用多项式相乘函数处理 格式 c conv a b 其中a和b分别表示一个多项式 c为a和b多项式相乘的结果多项式 例1 求c s A s B s s 2 2s 1 s 3 8s 2 6s 4 A 1 2 1 B 1 8 6 4 c conv A B 例2 求的传函表示形式 2 状态空间表达式的表示连续系统 ss A B C D 离散系统 ss A B C D T A B C D分别为系统状态空间表达式的状态矩阵 输入矩阵 输出矩阵 反馈矩阵 实际上 一组矩阵 A B C D 就能确定一个系统模型 例 已知系统状态空间表达式为 3 零极点模型的表示格式 zpk z p k 零极点模型对于研究系统稳定性 PID校正等辅助分析和综合方法十分有用 z表示系统零点向量 p表示系统极点向量 k表示系统增益矩阵 零极点模型的一般形式为 相应的增益为k 零极点向量分别为 2 模型在MATLAB中的转换1 化传递函数为状态空间表达式格式 A B C D tf2ss num den 输入参数num和den表示传递函数分子和分母系数构成的向量 输出参数A B C D表示状态矩阵 该函数用可控标准型方法实现转换过程 例 传递函数为 将其转换为状态空间表达式 2 化状态空间表达式为传递函数格式 num den ss2tf A B C D iu 输入参数A B C D表示状态矩阵 iu表示输入代号 对于SISO系统 iu 1 对于多输入系统 必须对各变量 输入信号逐个求取传函子矩阵 最后再获得整个的传递函数矩阵 求其传递函数形式 再反求状态方程 比较结果 试述原因 例 状态方程模型为 3 化零极点模型为传递函数格式 num den zp2tf z p k 4 化传递函数为零极点模型格式 z p k tf2zp num den 5 化零极点模型为状态空间表达式格式 A B C D zp2ss z p k 6 化状态空间表达式为零极点模型格式 z p k ss2zp A B C D 以上A B C D为系统状态矩阵 z p k为系统零极点和增益矩阵 num den为系统传函分子和分母向量 数值仿真算法的中心问题是如何将用微分方程描述的动力学系统模型转换为能在数字计算机上运算的类型一个连续系统的时域数学模型的基本形式为微分方程或状态空间表达式 高阶微分方程可转化为一阶方程组 即状态方程n阶连续系统的状态空间表达式的状态方程是由n个一阶微分方程组成的微分方程组 数值积分方法是数字仿真中对微分方程求解的基本算法 目前常用的数值积分方法有很多种 如 龙格 库塔法 C Runge M W Kutta 2 3数值积分法 任何微分方程的求解都是以积分计算为基础的 随着电子计算机的日益普遍 数值积分方法以成为近似求解的主要方法 假设一个系统的数学模型用状态方程表示为 2 3 1数值积分的基本概念 式中x为n维状态向量 其初始值 t0 0时 为x0 t为时间变量 2 3 1 对式 2 3 1 两边取积分 得 数值积分法的思想就是从初值x0开始 在一系列时刻求解x t1 x t2 x tn 的近似解x1 x2 xn 其中 t1 t0 t2 t1 tn 1 tn h h为计算步长 当t tn 1 t0 tn时 有 令 若f的各阶导数都存在 则Qn的近似解可由泰勒 Taylor 公式展开得到 如t t1时 2 3 3 即为Q0 数值积分法示意图 为精确计算x t1 除需计算x t 对t的各阶导数之外 还要计算无穷项之和 显然这十分困难 但在h很小的情况下 可以只式 2 3 3 的前两项 三项或四项之和作近似计算 从而得到x t1 的计算值 再把x t1 当做初值计算x t2 依次类推 就可逐步计算出x t 的数值解x1 x2 xn 实质上是用收敛的无穷级数的前几项代替了无穷级数 这等于抛弃了无穷级数的后段 所造成的误差叫截断误差 此外 由于受到机器字长的限制 用机器代码表示的数据必须舍掉一定的位数 所以对于计算位数很多 甚至是无穷小的数据 将引入舍入误差 少量运算的舍入误差一般是微不足道的 但在计算机上完成千百次运算之后 舍入误差的积累将是十分惊人的 因此在将来选择数值算法的时候 这两种误差是必须要考虑的 这两种误差中 截断误差是方法误差 舍入误差是工具误差 欧拉法是古典数值积分计算中最基本的方法 欧拉法只取式 2 3 3 的前两项之和来近似计算x t1 若初始条件为t t0 x t0 x0 则计算t t0 h时的x t1 x1 欧拉法的计算公式为 即 一般的 2 3 4 一般的 2 3 2欧拉法 Euler 式 2 3 4 是一个递推差分公式 只要给出初值x0 就可依据f t0 x0 和h计算出x1 再以x1计算出x2 依次递推计算出x3 x4 xn 完成数值积分过程 欧拉法的几何意义如上图所示 其实质是 在计算式 2 3 2 中的Qn 1式 将f t x 看做常数f tn 1 xn 1 从而用矩形的面积代替曲面面积Qn 1 产生面积误差 即 由式 2 3 2 得 即 2 3 4 欧拉法的计算公式是根据泰勒级数展开式截断以上的高阶项得到的 通常把截断项 称为截断误差 用截断误差可以表示数值计算方法的精度 截断误差为时 称方法为r阶的 对于欧拉法 r 1 故欧拉法是1阶精度的 特点 运算速度快 计算精度低 截断误差正比于 适用于精度要求低的场合 为提高精度就要减小h 这样不但运算量大 计算中的舍入误差也会加大 取式 2 3 3 泰勒级数中的前三项之和来近似计算 得 2 3 3梯形法 二阶龙格 库塔法 改写为 2 3 5 设 2 3 1 的解可以写成如下形式 2 3 6 将k2用二元函数泰勒级数展开式展开 取前三项 得 将k1 k2代入 2 3 6 得 2 3 7 比较式 2 3 5 与 2 3 7 得 显然 上式有非唯一解 假设a1 a2 则得到 a1 a2 1 2 b1 b2 1 代入 2 3 6 得 一般的 2 3 8 2 3 8 即为二阶龙格 库塔法 RK2 计算公式 预报公式 校正公式 一般在精度较高的情况下 多使用四阶龙格 库塔法 RK4 即 2 3 3 取前5项 得 2 3 4四阶龙格 库塔法 Runge Kutta 式中 1 一阶微分方程组 2 4一阶微分方程组的数字仿真程序 定义x1 x2 xn是n个状态变量 x10 x20 xn0是n状态变量的初值 2 MATLAB数值算法函数MATLAB提供了两个常微分求解函数 ode23 采用二阶 三阶龙格 库塔法 ode45 采用四阶 五阶龙格 库塔法 它们采用自动变步长的求解方法 即当求解变慢时 采用较大的计算步长 使计算速度变快 相反当方程的求解变化较快时 积分步长会自动变小 使计算精度提高 这种能自动变步长的方法称为自适应变步长方法 格式 T X ode23 系统函数名 Tspan X0 Options T X ode45 系统函数名 Tspan X0 Options 说明 1 输入参数中的 系统函数名 是一个字符串 表示微分方程的形式 函数的编写格式是固定的 不按格式写 会得到错误的结果 格式 functionxdot 函数名 t x 其中t为时间变量 x为状态变量 xdot为状态变量导数 2 Tspan T0 Tfinal 表示积分区间 T0为仿真起始时间 Tfinal为仿真终止时间 3 X0为系统状态变量的初值 函数表示在初始条件X0下 从T0到Tfinal对微分方程dx dt f t x 进行积分 函数f t x 必须返回一个列向量 两个输出参数是列向量t与矩阵x 其中向量t包含估计响应的积分点 而矩阵x的行数与向量t的长度相等 t向量中的积分点不是等间距的 因为为了保持所需的相对精度 积分算法改变了步长 4 Options为可选项 为积分参数 可由函数odeset来设置 Options最常用的是相对误差 RelTol 默认是1e 3 和绝对误差 AbsTol 默认是1e 6 格式 Options odeset RelTol 1e 3 AbsTol 1e 6 5 仿真结果可用函数plot t x 绘制出响应曲线来 6 用此方法要求已知系统的状态方程 要求解该方程组 首先要建立一个函数文件VDP M 函数的参数是时间t和一个二维向量 返回值是一个列向量 代表的是导数的值 最后用二维曲线表示输出响应曲线 解 首先对于上述微分方程 通过重新定义两个新的状态变量来实现从高阶微分方程到一阶微分方程组的转换 例1 描述振荡器的经典范德波方程 例2 已知系统模型为 出示状态为0 1 0 07 输入为单位阶跃 u 1 试绘制响应曲线 解 1 先求出状态方程 2 写系统函数 保存为 m文件 3 设置参数 并用ode45进行仿真 4 绘制输出曲线 plot t x 1 注意 由于龙格 库塔法的推导是在泰勒级数展开式的基础上 因而它要求所求的解具有较好的光滑性 假如解的光滑性较差 其精度可能反而不如改进的欧拉法公式 所以在实际计算中 应当针对问题的具体特点 选择合适的算法 2 6结构图法数字仿真 控制系统中常常是由许多环节构成的 应用前面所学内容 一般首先是用高阶微分方程或传递函数写出系统状态空间表达式 根据状态方程求取系统的动态响应 非常繁琐 尤其是当改变某一环节参数时 状态方程的系统矩阵A B C D都要重新计算 很不方便 这对于研究对象参数变化对整个控制系统的