控制系统的设计与仿真PPT课件

第4章控制系统的设计与仿真 4 1系统建模与仿真框图的创建4 2控制系统设计4 3控制系统的时域仿真4 4实例 倒摆系统的建模与仿真 4 1系统建模与仿真框图的创建 4 1 1问题描述 众所周知 现实世界中存在着各种不同的控制系统 对于线性时不变 LTI 系统 一般可以分为连续和离散系统 MATLAB中为用户提供了丰富的针对各种系统的建模手段 图4 1显示了MATLAB中各种线性时不变 LTI 系统之间的转换关系 图4 1连续与离散系统的关系示意图 图4 1中显示了MATLAB可以完成离散和连续系统的建模 并且同一系统可以表示成连续系统 也可以表示成离散系统 它们之间可以以状态方程形式进行转化 这一节将结合一个具体实例来演示MATLAB中各种模型创建和相互之间进行转化的方法 以及如何用Simulink进行连续系统的仿真 首先给出实例的源程序MODLDEMO M 然后根据不同的主题结合实例进行讲述 例4 1对于MassSpringDashpot机械系统 4 1 试建立该系统的连续和离散模型 并进行时域和频域仿真 解 程序源代码如下 MODLDEMO M演示各种建模与仿真 时域和频域 技巧 clearall closeall 程序开始 清空工作空间 关闭所有窗口 deletemodldemo out diarymodldemo out 打开二进制文件 disp MODLDEMO OUT DiaryFileforMODLDEMO M disp m 1 各种系统常数定义 k 1 单位kg s 2 c 2 02 51 20 0 单位kg s 创建系统的状态空间模型 disp StateSpaceModels km k m A1 01 km c 1 m A2 01 km c 2 m A3 01 km c 3 m A4 01 km c 4 m B 01 m C 10 D 0 sys1s ss A1 B C D sys2s ss A2 B C D sys3s ss A3 B C D sys4s ss A4 B C D 仿真系统的脉冲和阶跃响应 时域 t 0 2 15 y1 impulse sys1s t y2 impulse sys2s t y3 impulse sys3s t y4 impulse sys4s t figure 1 subplot 221 plot t y1 r title CriticalDamping Impulse grid xlabel Time ylabel SystemResponse subplot 222 plot t y2 r title OverDamping Impulse grid xlabel Time ylabel SystemResponse subplot 223 plot t y3 r title UnderDamping Impulse grid xlabel Time ylabel SystemResponse subplot 224 plot t y4 r title NoDamping Impulse grid xlabel Time ylabel SystemResponse y1 step sys1s t y2 step sys2s t y3 step sys3s t y4 step sys4s t figure 2 subplot 221 plot t y1 r title CriticalDamping Step grid xlabel Time ylabel SystemResponse subplot 222 plot t y2 r title OverDamping Step grid xlabel Time ylabel SystemResponse subplot 223 plot t y3 r title UnderDamping Step grid xlabel Time ylabel SystemResponse subplot 224 plot t y4 r title NoDamping Step grid xlabel Time ylabel SystemResponse disp hitanykeytocontinue pause 在Matlab中进行模型转化 对于m k 1 系统简化的传递函数为 G s 1 s 2 cs 1 z s p s disp TransferFunctionForm sys1t tf sys1s sys2t tf sys2s sys3t tf sys3s sys4t tf sys4s disp hitanykeytocontinue pause 将传递函数形式转化成零极点形式 disp Zero Pole GainForm sys1z zpk sys1t sys2z zpk sys2t sys3z zpk sys3t sys4z zpk sys4t disp hitanykeytocontinue pause 演示留数定理的使用 disp ResidueForm n1 d1 tfdata sys1t n2 d2 tfdata sys2t n3 d3 tfdata sys3t n4 d4 tfdata sys4t 注意tfdata命令定义元胞数组 将元胞数组转化成向量形式 n1 n1 1 d1 d1 1 n2 n2 1 d2 d2 1 n3 n3 1 d3 d3 1 n4 n4 1 d4 d4 1 使用留数定理 r1 pr1 kr residue n1 d1 r2 pr2 kr residue n2 d2 r3 pr3 kr residue n3 d3 r4 pr4 kr residue n4 d4 disp hitanykeytocontinue pause 进行频域仿真 这里使用nyquist函数 参考bode freqs等函数 w logspace 2 2 100 re im nyquist sys1s w SS形式 re1 1 re 1 1 im1 1 im 1 1 g1 re1 i im1 mag1 20 log10 abs g1 phase1 angle g1 180 pi re im nyquist sys2s w SS形式 re2 1 re 1 1 im2 1 im 1 1 g2 re2 i im2 mag2 20 log10 abs g2 phase2 angle g2 180 pi re im nyquist sys3t w TF形式 re3 1 re 1 1 im3 1 im 1 1 g3 re3 i im3 mag3 20 log10 abs g3 phase3 angle g3 180 pi re im nyquist sys4t w TF形式 re4 1 re 1 1 im4 1 im 1 1 g4 re4 i im4 mag4 20 log10 abs g4 phase4 angle g4 180 pi 幅值曲线 figure 3 subplot 221 semilogx w mag1 r title CriticalDamping xlabel frequency ylabel G jw indb subplot 222 semilogx w mag2 r title OverDamping xlabel frequency ylabel G jw indb subplot 223 semilogx w mag3 r title UnderDamping xlabel frequency ylabel G jw indb subplot 224 semilogx w mag4 r title NoDamping xlabel frequency ylabel G jw indb disp hitanykeytocontinue pause 频率曲线 figure 4 subplot 221 semilogx w phase1 r title CriticalDamping xlabel frequency ylabel angle subplot 222 semilogx w phase2 r title OverDamping xlabel frequency ylabel angle subplot 223 semilogx w phase3 r title UnderDamping xlabel frequency ylabel angle subplot 224 semilogx w phase4 r title NoDamping xlabel frequency ylabel angle disp hitanykeytocontinue pause 绘制系统Nichols图 对数坐标形式 figure 5 subplot 221 plot phase1 mag1 r title CriticalDamping xlabel angle ylabel G jw indb gtext w 0 subplot 222 plot phase2 mag2 r title OverDamping xlabel angle ylabel G jw indb gtext w 0 subplot 223 plot phase3 mag3 r title UnderDamping xlabel angle ylabel G jw indb gtext w 0 subplot 224 plot phase4 mag4 r title NoDamping xlabel angle ylabel G jw indb gtext w 0 disp hitanykeytocontinue pause 绘制Nyquist图 实部 虚部形式 figure 6 subplot 221 plot re1 im1 r title CriticalDamping xlabel Real ylabel Imag gtext w 0 subplot 222 plot re2 im2 r title OverDamping xlabel Real ylabel Imag gtext w 0 subplot 223 plot re3 im3 r title UnderDamping xlabel Real ylabel Imag gtext w 0 subplot 224 plot re4 im4 r title NoDamping xlabel Real ylabel Imag gtext w 0 disp hitanykeytocontinue pause 直接使用Nyquist函数绘制Nyquist图 figure 7 subplot 2 2 1 nyquist sys1s w title CriticalDamping subplot 2 2 2 nyquist sys2s w title OverDamping subplot 2 2 3 nyquist sys3s w title UnderDamping subplot 2 2 4 nyquist sys4s w title NoDamping disp hitanykeytocontinue disp hitcntrl ctogetoutofthefileforinteractiveanalysis pause 使用Simulink框图 文件名为mdemosl mdl 如图4 2所示 进行仿真 clearall 清空所有内存 m 1 k 1 c 1 2 定义常数 disp DatafromgraphicalSimulinkmodel 计算线性化模型 A B C D linmod mdemosl syss ss A B C D sysz1 zpk syss syssm minreal syss sysz2 zpk syssm diaryoff 关闭二进制文件 图4 2MassSpringDashpot系统仿真框图 4 1 2物理模型 例4 1中研究的对象是一个简单的物质交换机械系统 可以用微分方程表示成 4 2 其中 y t 是系统的瞬时交换的质量 k和c分别为比例常数 如果取状态变量x1 y和x2 dy dt 则可以得到系统的状态方程形式 可以看出 矩阵A B C和D可以看成LTI系统的标准状态方程形式 于是式 4 3 可以写成 4 3 4 1 3时域仿真 一般对控制系统进行时域仿真可以采用impulse step和lsim函数 这些函数用来处理系统的状态空间描述 例如 对于输入信号u t 系统的动态仿真可表示成 sys ss A B C D Y T X lsim sys U t xo 当然 这些函数也可以处理系统的频域表达形式 这时 LTI系统可以通过下式表示 sys tf num den 作为一个典型的例子 例4 1使用了impulse和step函数来仿真系统的时域特性 图4 3和图4 4分别为例4 1执行的仿真结果 图4 3典型二阶系统的脉冲响应曲线 图4 4典型二阶系统的阶跃响应曲线 4 1 4模型转换 LTI系统的频域描述可以用传递函数表示为 Y s G s U s 4 5 系统的传递函数矩阵为 4 6 对于例4 1研究的简单SISO机械系统 其传递函数可以写成简单的标量形式 4 7 但是对于多输入多输出 MIMO 系统而言 其传递函数描述就有些复杂了 这时MATLAB中的模型转换函数可以发挥作用 它可以完成系统在状态方程形式与传递函数形式之间的互换 同时也可以将传递函数形式转换成零极点 增益形式 相关的函数包括sys1 ss A B C D sys2 tf sys1 sys3 zpk sys2 我们也可以采用ssdata tfdata和zpkdata等命令将存储在与一个指定LTI对象相联的数据结构中的信息抽取出来 例如 num den tfdata sys2 返回LTI对象sys2的分子和分母多项式系数 num与den为相应的元胞数组 其行数为输出的维数 列数等于输入的维数 其中第i行第j列元素表示从第j个输入到第i个输出的传递函数 另一种从数据结构中得到元胞数组的方法是使用MATLAB的celldisp命令 简单的显示数据的信息 可以使用如下的命令 fieldnames sys2 num1 sys2 num den1 sys2 dencelldisp num1 celldisp den1 也可以对零极点 增益形式完成显示的操作 例如 Z P K zpkdata sys3 将获取LTI系统sys3每一个IO通道的零极点和增益大小 元胞数组Z P和矩阵K的行与列分别与输出和输入的维数相同 其中第i行第j列元素表示从第j个输入到第i个输出传递函数的零极点和增益 对于单输入单输出 SISO 系统 其传递函数与零极点 增益形式可以简化成普通的分数形式 即 4 8 4 1 5留数定理 我们也可以将系统写成几个分数相加的形式 例如对于SISO的机械系统 G s 可以写成 4 9 4 10 r1 r2可以通过各种不同的方法计算得到 对于上述问题 MATLAB可用residue函数来完成这一运算 即 r p k residue B A 其中B和A为包含多项式系数的行向量 而r和p是包含留数和极点的列向量 如果B s 比A s 的维数大 则k s 不为零 4 1 6频域仿真 系统的频域仿真在概念上是非常直观的 但是计算起来常常比较复杂 作为练习 读者可以针对不同的 值 计算下面的表达式 G j C j I A 1B D 4 11 然后通常采用下列三种方法来绘制频域曲线 Bode图 Nichols图和Nyquist曲线 MATLAB中的一些函数用来获取系统的频域信号 首先必须产生一个频率的向量 采样点的坐标通常采用对数形式 即从10d1到10d2共n个点 可以通过下面的命令来完成 w logspace d1 d2 n bode和nyquist函数可以用来计算每一个频率 所对应的G j 即 MAG PHASE bode sys w RE IM nyquist sys w 图4 5典型二阶系统的Bode幅值曲线 图4 6典型二阶系统的Bode频率曲线 图4 7 图4 8 4 1 7仿真框图及系统建立 许多设计系统都可以由一些基本的组件和框图中的反馈回路组成 在有些情况下 寻找系统的等价描述和相应的状态空间矩阵是比较困难的 幸运的是 我们可以借助MATLAB从系统的Simulink仿真框图直接建立它的状态空间描述 这些工作可以通过MATLAB中的控制工具箱或Simulink的图形仿真界面来完成 为了演示这一过程 同样考虑一下例4 1所述的简单机械系统 首先建立该系统的Simulink仿真框图 然后自动创建原系统的状态空间和整个系统的传递函数形式 下面将状态方程展开 并且进行Laplace变换 4 12 上述拉普拉斯变换可以对应于如图4 10所示的基本模块 现在将这些模块连接起来 定义输出为Y s X1 s 最后得到如图4 2所示的仿真框图 下面我们可以使用linmod函数来计算LTI系统的状态矩阵 A B C D linmod mdemosl 计算的结果与前面得到的结果一致 图4 9MATLAB计算产生的典型二阶系统的Nyquist曲线 图4 10Laplace变换下的基本模块 4 2控制系统设计 在经典控制系统设计中通常以线性系统模型为研究对象 对于一个线性时不变 LTI 系统 其状态方程可以描述为 4 13 4 14 这里已经假定系统的输出没有显式地包含输入变量u 即D 0 系统 4 14 也可以表示为传递函数形式 Y s G s U s G s C sI A 1B 4 15 一个LTI系统的控制系统方框图如图4 11所示 图4 11系统的线性状态方程模型 4 2 1经典比例控制器 在经典控制系统的例子中 首先来看一下图4 12所示的一个简单的闭环系统 对于一个SISO系统而言 系统传递函数G s 仅仅是式 4 15 所示的标量函数 该传递函数嵌入在图4 12所示的方框图中 反馈回路包含传感器传递函数H s 而控制器部分只有简单的增益环节Kc组成 rd是闭环系统期望的响应或参考点 图4 12SISO系统的经典比例控制器框图 图4 12所示闭环系统的传递函数可以写成 4 16 其中 Gc s 为闭环传递函数 Kc是经典比例增益 对于单位反馈情况有H s 1 Gc s 可以简化为 4 17 下面是标量输入函数的时域表示 u t Kc rd t y t Kc rd t CTX t 4 18 从而式 4 14 可以写成 4 19 这里的参考点rd成为系统的一个独立输入变量 既然控制器只有唯一的参数Kc需要确定 因此该系统的控制器设计比较简单 闭环系统的暂态响应由状态方程系数矩阵的特征值或者整个系统的根极点确定 我们可以在时域中通过选择合适的控制参数Kc 使得 A KcBCT 的特征值产生期望的暂态响应 上升时间 最大超调量等 与此类似 也可以在传递函数中通过选择合适的控制参数Kc来设计式 4 17 系统的根极点位置 这两种设计方法是等价的 我们知道Gc s 的极点是1 KcG s 的根 因此可以将极点配置方程看作控制增益Kc的根 运用根轨迹方法可以确定满足设计要求的控制参数 4 2 2状态反馈控制器 上述经典控制器的主要不足是系统仅有唯一的控制参数Kc可供调整 而对于N维控制系统 系统开环矩阵具有N个特征值或者开环传递函数具有N个极点 即 det A I 0或det sI A 0 4 20 要想将所有这些系统根极点调整到需要的位置 控制器至少需要N个独立变量 因此仅仅将系统输出信号进行反馈将不能满足控制器设计的要求 一个自然的想法就是将系统的所有状态变量X都进行反馈 这就产生了状态反馈控制器 对于SISO系统 状态反馈后的系统输入变成 u t rd t KTsX t 4 21 Ks称为系统的反馈系数 这样 闭环系统的状态方程可以写成 4 22 闭环系统的框图如图4 13所示 同时 图4 11所示的状态反馈系统变成图4 14所示的仿真框图 图4 13SISO系统的状态反馈控制器 图4 14 4 2 3完全可控性 为了设计具有状态观测器的状态反馈控制器 让我们首先熟悉有关系统可控性的定义 假设一个SISOLTI系统由式 4 23 描述 4 23 如果该系统能够构造一个无约束的输入信号u t 使得系统能够在有限的时间间隔内 t0 t tf 由初始状态运动到任何其它的状态 则可以说系统在t0时刻是可控的 如果系统的每个状态都是可控的 则称该系统是完全可控的 不失一般性 假设X tf 0 t0 0 则 4 24 根据完全可控性的定义 有 4 25 或者 根据Sylvester积分公式 有 4 26 4 27 4 28 当如下矩阵非奇异时 系统满足完全可控的条件 M BABA2B AN 1B 4 30 或者 4 29 4 2 4极点配置 设计状态反馈控制器的最简单方法是采用极点配置 其基本思想是首先确定闭环系统N个根极点的期望位置 然后设计适当的反馈增益 从而将系统的极点调整到期望的位置 如果系统是完全可控的 则这一过程完全可以表示成包含N个未知参数的N个方程组的求解 所需要设计的反馈控制增益就是该方程组的解 如果系统比较简单 则完全可以通过手工计算完成系统的极点配置 但无论是手工计算 还是通过MATLAB函数自动计算 其基本步骤都是相同的 如下所示 1 检查系统的可控矩阵是否满秩 2 确定闭环系统的期望极点 1 2 N 3 确定希望配置的极点位置后 可以建立期望的特征方程 4 最后建立闭环系统的特征方程 即 sI A BKTs 0 将 3 4 步建立的方程联立 由于其多项式的系数相等 由此可以建立N个位置参数的N个方程组 从而可以唯一地确定系统的反馈增益矩阵KTs 例4 2假定SISOLTI系统的状态方程为 闭环系统的期望极点为 1 2 1 8 2 4j 试设计确定系统状态反馈的增益矩阵 解 首先观察开环系统的极点 可以看出 系统开环极点为s1 2 4 539 系统是不稳定的 闭环系统的期望极点是由期望的系统暂态响应特性 上升时间 读者可以验证 1 2 1 8 2 4j的闭环极点将产生较好的动态特性 大约10 的最大超调量和大约0 6s的上升时间 因此 期望的闭环系统极点是不唯一的 下面在已经确定期望闭环系统极点的情况下来设计系统的反馈增益矩阵 Step1 验证系统的可控性 矩阵M的秩等于N 因此系统满足完全可控性条件 Step2 3 计算期望的特征方程 s 1 s 2 s 1 8 j2 4 s 1 8 j2 4 s2 1s 2 0 Step4 计算闭环系统的特征方程 因此 4 2 5带全观测器的状态反馈控制 设计状态反馈控制器的主要问题是要求系统的所有状态变量都是可测的 然而对于一个实际系统而言 有些状态的信号值很难测量甚至不可能直接通过传感器进行测量 或者虽然可以进行直接测量 但在经济上却要增加相应的成本 这样 如果不能得到系统的全状态向量 前面讲述的状态反馈控制就不可能实现 解决以上问题的方法是利用系统某种数学形式的仿真来估计不能测量的状态值 这种方法称之为系统的状态观测器设计 下面假定以SISOLTI系统为研究对象 这意味着系统有唯一的可控变量和唯一的可测量 同时 假定系统输出y t 是唯一能够测量的量 它将被引入到状态观测器中来提高状态值的估计过程 这里采用来表示状态向量X t 的在t时刻的估计值 状态观测器的框图如图4 15所示 注意变量xc表示 该观测器使用u t 和y t 作为输入量 并且输出系统状态关于时间的估计值 从框图中可以看出 4 31 这里的L为未知的增益 它是根据该子系统期望的暂态响应特性确定的 称为状态观测器的增益矩阵 对于SISO系统 L是长度为N的列向量 图4 15SISO系统的状态观测器模型 观测器的设计过程与前面讲述的标准状态反馈控制器类似 这里的观测器的增益选择应使状态观测器的特征值是稳定的 同时使得观测器的动态变化快于整个闭环系统的动态属性 观测器的特征值由下式确定 det sI A LCT 0 4 32 在状态反馈控制系统中加入状态观测器 可得到图4 16所示的系统框图 对于该系统 系统输入为 4 33 如果系统模型与状态观测器模型都采用相同的状态空间矩阵A B CT进行描述 则对于被研究对象有 将式 4 33 的输入代入上面的方程 则得到系统的完整模型 4 34 4 35 图4 16具有全状态观测器的状态反馈控制系统框图 SISO系统 对于状态观测器 将式 4 33 代入式 4 31 可以得到观测器的完整模型 或者 4 36 4 37 定义误差向量 将式 4 37 代入得到误差向量的动态模型 4 38 4 39 4 2 6完全可观性 如果系统的每个状态X t0 都可以通过y t 一段时间的观测值来确定 则该系统被称为是完全可观的 考虑式 4 23 定义的SISOLTI系统 其时域解为 4 40 假设u 0 为方便计算 上式可以写成 4 41 其中 CTeAt已知 y t 可测 因此状态向量X 0 可以通过y t 的观测值间接计算得到 对于SISO系统 方程 4 41 仅有一个方程 但包含N个未知参数 然而 由于该方程独立于时间变量 因此 在多个时刻对y t 进行测量 可以获得多个类似的方程 将它们联立 就可以唯一确定系统的初始状态X 0 下面给方程 4 41 两边同时乘以已知的系数矩阵 得到 4 42 4 43 4 44 变换方程有 将整个观测时间内的所有方程进行联立 得到 4 45 其中 4 46 最后求解方程 4 45 得到 如果W是非奇异的 X 0 可以由y t 的观测值唯一确定 从而系统是完全可观的 再次使用Sylvester的积分公式 得到 定义 4 48 4 49 4 50 4 2 7观测器增益的确定 确定观测增益矩阵同样采用极点配置的方法 然而在这里 我们确定的是状态观测器误差方程的极点位置 误差极点位置的选择比较随意 但误差动态变化应该比被控系统的动态变化快一些 如果系统完全可观 则 A LCT 的N个特征值的位置应该唯一确定观测器增益矩阵的N个元素 设计的过程如下 1 检查系统可观矩阵是否奇异 2 为误差方程指定期望的极点位置 1 2 N 这些极点位置与系统的主导极点相比较应该更靠近复平面的左手边 3 根据期望的极点位置创建期望的特征方程 4 最后创建误差方程的特征方程 从而得到含有N个位置参数的方程组 同样 以例4 2为例来说明SISO系统状态观测器的设计过程 例4 4用MATLAB的place函数重新设计例4 2的状态反馈控制器与例4 3中的全状态观测器 解 程序代码如下 SFSOTEST MSISOLTI系统的状态反馈控制器与全状态观测器的设计 clearall closeall nfig 0 打开二进制文件保存结果 deletesfsotest out diarysfsotest out disp disp SFSOTEST OUT DiaryFileforSFSOTEST M disp PartI 创建线性系统模型 显示其开环系统是不稳定的Setupbasedataforthelinear A 01 20 60 B 01 建立系统的状态空间模型 C 10 D 0 disp StateSpaceMatricesforthePlant A B C D computeeigenvaluesofstatematrixforopenloopplant disp Eigenvaluesofthe OpenLoopPlant 计算开环系统状态矩阵特征值 ev eig A PartII 加入状态反馈控制器以稳定系统 对状态1的输出量进行仿真 检查系统的可控性 disp ControllabilityMatrixforthissystem M ctrb A B disp RankofControllabilityMatrix rank M clp 1 8 2 4j 1 8 2 4j 计算状态反馈控制增益 Ks place A B clp disp Desiredclosedlooppolesforstatefeedbackcontroller clpdisp Statefeedbackgainsneededtogivedesiredpoles Ks disp Calculatedeigenvaluesofsystemwithstatefeedback eig A B Ks Nv 1 0 C inv A B Ks B 计算Nv disp SetpointgainforzeroSSerror Nv to 0 tf 5 nfig 0 对被控系统 控制器进行仿真 t linspace to tf 101 syscl1 ss A B Ks B Nv C D y1 t x1 step syscl1 t nfig nfig 1 figure nfig 绘制相关结果曲线 subplot 2 1 1 plot t x1 1 r t x1 2 g grid title StatesforStateFeedbackTestCase xlabel Time ylabel StateVariables legend x1 t x2 t PartIII 加入状态反馈控制器与全状态观测器 仿真状态1的阶跃响应特性 检查系统的完全可观性 disp ObservabilityMatrixforthissystem H obsv A C disp RankofObservabilityMatrix rank H 针对指定的观测器极点计算观测器增益 op 3 clp 观测器的速度是闭环系统的3倍 L place A C op L L disp Desiredobserverpolesforstatefeedbackcontroller op disp Estimatorgainsneededtogivedesiredpoles L disp Calculatedeigenvaluesofestimatorsystem eig A L C A11 A A12 B Ks B1 B Nv A21 L C A22 A L C B Ks B2 B Nv zz 0 AB A11A12 A21A22 BB B1 B2 CB Czz C 对控制系统进行仿真 syscl2 ss AB BB CB D y2 t x2 step syscl2 t 分离系统状态与观测器的估计状态 nn max size A xp2 x2 1 nn xe2 x2 nn 1 2 nn 绘制全观测器情况下的计算结果 subplot 2 1 2 plot t xp2 1 r t xp2 2 g grid title StatesforStateFeedbackwithFullObserverTestCase xlabel Time ylabel StateVariables legend x1 t x2 t 绘制状态估计的误差 nfig nfig 1 figure nfig plot t xp2 1 xe2 1 r t xp2 2 xe2 2 g grid title DifferenceBetweenPlantandObserverStates xlabel Time ylabel ErrorinStateVariables legend e1 t e2 t Diaryoff 关闭二进制文件 图4 17系统模型与观测器动态性能演示 图4 17系统模型与观测器动态性能演示 图4 18状态误差的动态曲线 4 2 8对偶系统 回忆前面论述的系统可控与可观性的定义 可知系统的可控性要求其可控矩阵满秩 反映的是状态矩阵A与输入矩阵B之间的关系 而系统的可观性要求其可观矩阵满秩 反映的是状态矩阵A与输出矩阵CT之间的关系 这两个概念实际上反映了控制系统的对偶原则 对于一般的MIMO系统 4 52 4 54 4 53 对于系统1 其完全状态可控的充要条件是M矩阵满秩 即 M1 BABA2B AN 1B 4 55 系统完全状态可观性的充要条件是H 矩阵满秩 即 H 1 C A C A 2C A N 1C 4 56 对于系统2 其完全状态可控的充要条件是M矩阵满秩 即 M2 C A C A 2C A N 1C 4 57 系统完全状态可观性的充要条件是H 矩阵满秩 即 H 2 BABA2B AN 1B 4 58 综上所述 给定系统的可观性可以通过其对偶系统的可控性来检验 而研究系统的可控性则可以通过其对偶系统的可观性来研究 这些性质称为系统的对偶原则 4 3控制系统的时域仿真 前一节以经典控制器和状态反馈控制器为例讲述了控制系统控制参数的设计过程 设计过程主要依据的是系统的时域特性 即闭环系统状态矩阵的特征值决定了闭环系统的暂态响应特性 控制器设计的目标是选择控制器的增益 使得闭环系统状态矩阵的特征值位于期望的极点位置 控制器参数确定以后 下一步需要对闭环系统进行仿真 在设计过程中研究的对象一般是系统的线性化模型 而在仿真过程中 应该尽可能准确地再现实际的系统模型 这就常常要求以实际的时变或非线性系统为仿真的对象 这一节主要论述控制器参数确定以后的系统仿真过程 为简单起见 仍然以SISO系统为研究对象 其线性与非线性模型分别为 4 59 4 60 4 3 1经典比例控制器设计 单位反馈回路的简单比例控制系统框图如图4 19所示 系统控制输入为 4 61 使用线性化模型的闭环仿真方程为 4 62 写成标准的状态方程形式 4 63 如果采用非线性模型进行仿真 必须使用MATLAB中标准的ODE求解器来完成 为此 用户必须定义一个函数文件 MATLAB中的ODE求解器调用该函数来完成非线性系统的仿真 该函数文件包含下面的内容 1 指定t时刻的参考输入rd t 2 计算系统t时刻的输出y t CTX t 3 确定t时刻的输入u t Kc rd t y t 4 计算t时刻状态变量的导数 图4 19SISO系统的经典比例单位反馈控制框图 4 3 2状态反馈控制器 带全状态观测器的状态反馈控制框图如图4 20所示 该框图类似于图4 16 不同之处在于这里的模型包括一个附加的稳态状态增益模块 其中包含一个归一化的增益变量Nr 我们再来看看该系统线性化模型与非线性模型的仿真方程 并且重新调整控制器的增益 KTs L和Nr等 控制规律写成 4 64 该系统的闭环模型包括实际的系统状态X t 以及估计的系统状态 这样 该系统具有2N个未知参数 对于线性系统模型 该系统的完整模型为被控系统模型 4 65 观测器模型 4 66 写成标准的状态方程形式为 图4 20具有全状态观测器的状态反馈控制系统 4 3 3混合控制器 下面将研究一种将状态反馈控制运用到经典控制器中的例子 在这个例子中 状态反馈制器为经典控制器提供修正的参考信号 例如 对于图4 18的经典控制系统 如果选择的增益Kc不满足系统的设计要求 我们可能会手动调节参考点rd 以使系统具有更好的暂态响应特性 可以设想一种控制方案来自动为参考信号rd t 提供实时的修正信号rdm t 控制框图如图4 20所示 图中左边的sum模块的输出为修正后的参考点信号 在此之后 该系统看上去就像是具有单位反馈回路和比例增益Kc的经典控制器 在这个控制方案中 状态观测器和状态反馈回路的目标是提供一个修正的参考信号 来改善经典控制系统的控制效果 这种方法称为状态反馈辅助的经典控制方法 SFACC 也称为混合的控制器 图4 21显示的控制系统实际上可以看作是前面讲述的状态反馈控制系统 其控制框图可以认为是从图4 20的标准形式通过某些模块变换发展而来 图4 22 4 24显示了其中的变换 主要步骤包括 1 在系统输入前加上一个比例增益模块 并且在相应的引入回路中抵消它的影响 如图4 22所示 图4 21具有嵌入经典控制器的状态反馈控制框图 SFACC 2 对系统输出y构造一个负反馈 它可以对估计输出的正反馈加以抵消 如图4 23所示 3 分离出经典控制器部分 图4 24中的阴影部分 4 最后定义修正的状态反馈增益Ksm和修正的稳态状态归一因子Nnn 4 68 4 69 图4 22在SFC中加入经典控制增益 图4 23加入经典负反馈回路 图4 24分离出经典比例控制器部分 图4 21所示的SFACC控制系统与传统的状态反馈控制系统SFC 如图4 20所示 在功能上是一致的 仅有的近似存在于第 2 步 这一步假设估计的输出与实际的系统输出完全抵消 如果状态观测器设计成相对控制系统有足够快的反应时间 这种近似是可以满足要求的 对于线性模型 从以上的讨论中可以看出 SFC和SFACC具有相同的暂态响应 然而 SFACC算法比前者具有更好的控制效果 并且SFACC可以方便地在原先经典控制系统的基础上进行扩展 而不需要破坏原来已经存在的控制系统 这让用户可以从经典控制器直接过渡到现代控制系统 为了对SFACC控制系统进行仿真 我们假定控制规律为 修正后的参考信号为 将式 4 71 代入式 4 70 有 4 70 4 71 4 72 与SAC控制系统类似 SFACC闭环系统既包含实际状态也包含估计状态 因此 可以获得2N个未知参数 对于线性系统模型 有系统模型为 4 73 4 74 4 75 或者 4 76 写成标准的状态方程形式为 其中 4 78 4 3 4参考输入增益的调整 在前面的框图中 为了减小比例控制器带来的稳态误差 引入了一个输入归一化模块 该模块的增益值可以根据下面的算式进行确定 图4 25显示了系统的简化模型 线性模型 4 79 4 80 控制律 系统稳态时 我们希望yss与参考点rdss相同 简化系统稳态时的状态方程为 4 81 求解该方程 得稳态时的状态向量为 4 82 图4 25具有SS增益模块的状态反馈控制系统 稳态输出为 4 83 最后将yss rdss代入 得到Nr的计算公式为 4 84 4 4实例 倒摆系统的建模与仿真 4 4 1引言 这一节我们先看一个实际的例子 图4 26是某个倒摆系统的示意图 倒摆通过转动关节安装在驱动小车上 通过对小车施加一定的驱动 使倒摆保持一定的姿态 这是姿态控制问题的典型代表 图4 26驱动小车上的倒摆示意图 4 4 2基本方程 假定倒摆由无质量的轻杆和质量为m的小球组成 小车的质量为M 系统所受的外力包括小球受到的重力和对小车水平方向的驱动力u t x t 和 t 分别表示小车的水平坐标和倒摆偏离垂直方向的角度 如图4 27所示 根据牛顿运动学第二定律 在水平x轴方向上系统满足下面的方程 4 85 图4 27小球受到的力矢量图 其中 小球受到的重力为常值 其中心的坐标表示为 xG yG 满足xG x lsin yG lcos 4 86 l为倒摆的杆长 将式 4 86 代入式 4 85 得到 4 87 由下面的基本关系 4 88 4 89 得到 同样 通过小球的力矩平衡关系 可以得到小球的平衡方程 4 90 4 91 由式 4 86 4 89 可以计算得到Fx和Fy的表达式 4 92 4 93 将式 4 92 4 93 代入式 4 91 中 两边消去l 最后得到 4 94 方程 4 90 和 4 94 构成了倒摆系统的动力学模型 显然 从数学的角度上看 该系统为明显的非线性系统 但是对小车施加驱动力的目的是保持倒摆在垂直方向上的姿态 因此 我们关注的是倒摆在垂直方向附近的动态变化 为此 可以将系统在该参考位置 0 进行线性化 这样可以对简化得到的线性模型进行控制器的设计 在后面的内容中 将分别对线性化前后的模型进行比较和研究 4 4 3非线性状态方程描述 为了对倒摆系统的非线性模型进行数值仿真 需要首先将模型表示成标准的状态方程形式 由式 4 94 可以得到 4 95 4 96 将其代入式 4 90 有 4 97 4 98 将其代入式 4 90 中 同样 可以由式 4 94 得到 4 99 上式满足式 4 95 的标准状态方程的形式 假定系统输出为倒摆的角度和小车的x轴坐标 则系统的输出方程为 4 102 4 4 4线性状态方程描述 为了得到倒摆在垂直位置附近的线性化模型 我们只要对式 4 97 描述的非线性状态方程进行线性化即可 根据前几节论述的方法 该系统的线性化模型为 4 103 这里的参考状态为倒摆静止在垂直方向 并且小车的驱动力为零 即Z0 0 u0 0 下面来逐项计算Jacobian矩阵 Jz Z0 u0 的第一列由决定 其中第一和第三个函数与z1不直接相关 其偏导数等于零 对于第二项 有 同样 对第四项有 Jz Z0 u0 的第二列由fi z2 z0 u0决定 式 4 101 显示f1 z2 1 f3 z2 0 对于第二和第四项有 Jz Z0 u0 的第三和第四列中仅有的非零项为 4 4 5倒摆系统的MATLAB仿真 1 倒摆系统的开环与闭环仿真这一部分的仿真主要研究下面三方面的内容 1 对倒摆非线性模型与线性化模型不稳定动态特性的比较 2 线性化模型可以由MATLAB函数linmod直接计算得到 以及如何用S函数表示非线性模型 3 列举一些不够理想的经典反馈控制方案 仿真计算的M文件包括INVPN1 M INVPNNL1 M和INVPNNL2 M等等 其中 INVPN1 M为主文件 主要功能是控制仿真计算的流向同时完成以上的分析过程 第 1 步比较将检验倒摆线性化系统与非线性系统的阶跃响应 倒摆的非线性模型在INVPNNL1 M文件中进行定义 仿真采用MATLAB的ode23函数计算模型中状态变量的导数 由于模型是不稳定的 因此只仿真系统1s内的动态行为 仿真结果如图4 28所示 从图中可以发现 随着小车水平指标的增长 倒摆将沿逆时针方向偏离垂直状态 这与实际经验相吻合 同时可以观察到 在仿真的最后 非线性模型与线性化后的模型动态行为的差异越来越大 这是因为对系统的线性化是在倒摆垂直方向 Z0 0 进行的 因此这种近似只是在该点附近有意义 图4 28倒摆开环系统的阶跃响应 INVPN1 M文件的第 2 步将演示如何在MATLAB Simulink环境中进行动态系统的仿真 在MATLAB仿真中 一个Simulink仿真方框图可以被转换成相应的S函数 该函数具备所有的仿真方程 可以用来计算状态变量的导数 输出函数以及系统的初始状态 随后 转换后的S函数可以用来进行仿真 使用sim命令 或者计算线性化模型 使用linmod 在这个例子中 S函数直接在INVPNNL2 M文件中构造 它可以通过Simulink中的Sfunction模块连接到仿真方框图中 尽管某些具体语法不同 但S函数文件INVPNNL2 M所包含的基本方程与ode23函数文件INVPNNL1 M中所使用的方程基本相同 有关S函数的使用说明请参考MATLAB Simulink的有关文档和帮助 可以参考sfuntmpl模板文件 在这个例子中 sim和linmod函数分别被用来仿真非线性系统的阶跃响应和计算参考点附近的线性模型 将其结果与前面MATLAB的标准分析结果相比较 可以发现二者是