控制系统的MATLAB仿真

通过本章学习 应该掌握以下内容 第6章MATLAB7控制系统分析 控制系统时域分析控制系统频域分析控制系统稳定性分析利用MATLAB语言进行编程仿真 1 控制系统时域响应仿真的主要问题 6 1 1概述 6 1控制系统时域分析的MATLAB实现 时域分析法是根据自动控制系统微分方程 用拉普拉斯变换求解动态响应的过程曲线 典型的动态过程响应有单位阶跃响应 单位斜坡响应 单位加速度响应与单位冲激响应等 时域分析的另外一个目的是求解响应的性能指标 2 时域分析MATLAB实现的方法 利用MATLAB所提供的求取连续系统的单位阶跃响应函数step 单位冲激响应函数impulse 零输入响应函数initial 以及其他函数都可以求出其相应的响应 MATLAB也提供了相应离散系统的单位阶跃响应函数dstep 单位冲激响应函数dimpulse 零输入响应函数dinitial 以及其他函数也可以求出其相应的响应 连续系统的单位阶跃响应函数为step 函数 该命令的调用格式如下 step sys step sys t step sys iu step sysl sys2 sysn y t x step sys 6 1 2求连续系统的单位阶跃响应 step sys 函数用来计算系统的单位阶跃响应 可用于SISO 单输入单输出 或者MIMO 多输入多输出 的连续时间系统或者离散时间系统 计算离散时间系统阶跃响应的函数为dstep 当函数命令为无等式左边输出变量的格式时 函数在当前图形窗口中直接绘制出系统的阶跃响应曲线 LTI对象sys可以是由函数tf zpk ss 中任何一个建立的闭环系统模型 step sys t 函数用于计算系统的阶跃响应 在函数中 可以指定为一个仿真终止时间 此时t为一标量 也可以设置为一个时间矢量 例如用t 0 dt Tfinal命令 若是离散系统 时间间隔出必须与采样周期匹配 step sys iu 函数可以绘制从第iu个输入到所有输出的单位阶跃响应曲线 step sysl sys2 sysn 函数可同时仿真多个LTI对象 注意 输出变量 y t x 三个元素的顺序不能错 y t x step sys 函数为带有输出变量引用的函数 可计算系统阶跃响应的输出数据 而不绘制出曲线 输出变量y是系统的输出响应矢量 输出变量t为取积分值的时间矢量 输出变量x是系统的状态轨迹数据 对于MIMO系统 y的维数为 1engthoft numberofoutputs numberofinputs x的维数为 1engthoft numberofstates numberofinputs 解 则可以将其用下列MATLAB语句表示 num 31 den 251 G tf num den step G 30 例 补充 现有如下SISO 单输入单输出 系统 其传递函数如下 试绘制其阶跃响应曲线 step G r 30 在绘制的系统阶跃响应曲线上 单击曲线上一点 则可以显示该点对应的时间信息和幅值信息 解则可以将其用下列MATLAB语句表示 num 1 2 21 den 31 251 h tf num den step h 例 补充 现有如下MIMO 多输入多输出 系统 其传递函数如下 试绘制其阶跃响应曲线 离散系统的单位阶跃响应函数dstep 的调用格式为 y x dstep num den n 或 y x dstep G H C D iu n 式中n为选定的取样点个数 当n省略时 取样点数由函数自动选取 其余参数定义同前 离散系统的单位阶跃响应 解1 可以将其用下列MATLAB语句表示 num 31 den 1251 dstep num den 40 例 补充 现有系统传递函数如下 采样周期为0 1 试绘制其阶跃响应曲线 求出了系统的阶跃响应曲线 要分析系统的稳定性能 还需要计算系统的性能指标 包括峰值时间tp 超调量 调节时间ts 及稳态误差ess等 这些参数可以直接从图上读取 也可通过MABLAB的M函数来求取 直接读取数据比较方便 但由于是人为读取的 存在一定的误差 这里介绍利用M函数来计算这些性能指标 function mp tp ess b1 b2 sigma n pusi T f targ y t targ函数中mp为曲线的峰值 tp为峰值时间 ess为阶跃响应的余差 b1为第一 正向 峰值的超调量 b2为第二 正向 峰值的超调量 sigma为阶跃响应的超调量 n为阶跃响应的衰减比 pusi为阶跃响应的衰减率 T为衰减振荡周期 f振荡频率 function mp tp ess b1 b2 sigma n pusi T f targ y t mp tf max y ct length t mp tp t tf tm max tf yss ct q 1 m q 1 whilem0 m m 1 pm m y i tp m t i endendendendyss y ct ess 1 yss b1 pm 1 yss b2 pm 2 yss sigma 100 b1 yss n b1 b2 pusi b1 b2 b1 T tp 2 tp 1 f 1 Ttp tp 1 例6 1已知典型二阶系统的传递函数为 试绘出当 0 5 分别为2 4 6 8 10 12时的单位阶跃响应曲线 解 依题意 编写MATLAB程序 M文件 如下 w 2 2 12 设定w的值kesai 0 5 figure 1 holdon 同时显示各条曲线forwn wnum wn 2 den 12 kesai wnwn 2 step num den endholdoff 重新设置坐标轴的属性gridon title 单位阶跃响应曲线 xlabel 时间 ylabel 振幅 6 1 3求连续系统的单位冲激响应 求连续系统单位冲激响应的函数为impulse 该函数命令的调用格式如下 impulse sys impulse sys t impulse sys iu impulse sysl sys2 sysN y t x impulse sys 例6 2已知单位正反馈系统前向通道的传递函数为 试绘出系统的单位阶跃响应和脉冲响应曲线 并计算系统的阶跃响应性能指标 解 1 求阶跃响应曲线 在命令窗口输入 num 4 den 110 sys tf num den closys feedback sys 1 求闭环传递函数 step closys 可得如图所示的阶跃响应曲线 2 求脉冲响应曲线在命令窗口继续输入 impulse closys 可得如图所示的脉冲响应曲线 3 计算系统的阶跃响应性能指标首先在MATLAB程序编辑窗口编辑前面给出的targ函数 并保存函数名为targ m 然后继续在MATLAB命令窗口输入以下语句 y t step closys mp tp ess b1 b2 sigma n pusi T f targ y t 运行完上述语句后 可得以下性能指标 mp 1 4441tp 1 6062ess 0 0017b1 0 4425b2 0 0857sigma 44 1728n 5 1651pusi 0 8064T 3 2125f 0 3113 6 1 4求连续系统的零输入响应 求连续系统零输入响应的函数为initial 该函数命令的调用格式如下 y t x initial sys initial sys x0 initial sys x0 t initial sysl sys2 sysN x0 例6 3已知一个二阶系统 试在同一个坐标系里绘制该系统所对应的三组不同参数配合下的阶跃响应曲线 解 由已知系统的数学模型 并根据题目要求 编制如下程序 c 124 k 1 25229 y zeros 100 3 构造一个100 3的零矩阵t linspace 0 10 100 在0到10之间构造100个线形分布的向量forj 1 3num k j den 1c j k j sys tf num den y j step sys t endplot t y 1 3 grid 绘制阶跃响应线 并添加网格线gtext a 1 25b 1 在鼠标指定的位置加文字gtext a 2b 2 gtext a 29b 4 c 124 k 1 25229 y zeros 100 3 构造一个100 3的零矩阵t linspace 0 10 100 在0到10之间构造100个线形分布的向量forj 1 3num k j den 1c j k j sys tf num den y j step sys t endplot t y 1 t y 2 o t y 3 x grid 绘制阶跃响应线 并添加网格线gtext a 1 25b 1 在鼠标指定的位置加文字gtext a 2b 2 gtext a 29b 4 程序执行后 即可得到系统的阶跃响应曲线 同时鼠标成大的十字形 将鼠标放到相应的位置 单击鼠标 则把说明文字加到相应的曲线上 如图所示 6 2控制系统频域分析的MATLAB实现 6 2 1概述 1 控制系统频域响应仿真的主要问题 对于频域分析 就是用先进的计算机仿真技术来精确绘制Bode图 Nyquist曲线图 Nichols曲线图等多种曲线 并计算系统的频域性能指标 剪切频率 也叫截止频率 穿越频率 相角稳定裕度 幅值稳定裕度 以便研究系统控制过程的稳定性 快速性及稳态精度等性能 2 频域分析MATLAB实现的方法 利用MATLAB软件 运行它提供频率分析的函数 能够方便 简单 快捷地绘制各种曲线 Bode图 Nyquist曲线 Nichols曲线图和根轨迹图等 并计算出频域性能指标 还可以借助于这些曲线对系统进行分析 6 2 2求连续系统Bode图 求连续系统Bode图的函数为bode 该函数命令的调用格式为 bode sys bode sys iu bode sys w bode sysl sys2 sysn bode sysl sys2 sysn w mag phase w bode sys 1 求连续系统Bode图的函数 bode sys 函数用来计算并显示绘制的系统Bode图 可用于SISO或者MIMO的连续时间系统或者离散时间系统 当函数命令为无等式左边输出变量的格式时 函数在当前图形窗口中直接绘制出系统的Bode图 LTI对象sys可以是由函数tf zpk ss 中任何一个函数建立的系统开环模型 bode sys iu 可得到从系统第iu个输入到所有输出的Bode图 bode sys w 函数用于显示绘制的系统Bode图 函数中输入参数w用来定义绘制Bode图时的频率范围或者频率点 若要定义频率范围 w必须为 wmin wmax 格式 如果定义频率点 则w必须为由需要频率点频率构成的向量 bode sysl sys2 sysn w 函数可同时在一个窗口绘制多个LTI对象的Bode图 这些系统必须有同样数量的输入与输出 但可以同时含有连续与离散时间系统 mag phase w bode sys 函数为带有输出变量引用的函数 可计算系统Bode图的输出数据 而不绘制出曲线 输出变量mag是系统Bode图的振幅值 输出变量phase为Bode图的相位值 输出变量w为Bode图的频率点 2 求离散系统Bode图的函数dbode 求离散系统Bode图的函数为dbode 该函数命令的调用格式为 mag phase w dbode a b c d Ts mag phase w dbode a b c d Ts iu mag phase w dbode a b c d Ts iu w mag phase w dbode num den Ts mag phase w dbode num den Ts w 6 2 3绘制系统Nyquist曲线的函数 计算 绘制 系统Nyquist曲线的函数为nyquist 该函数命令的调用格式为 re im w nyquist sys nyquist sys nyquist sys w nyquist sys iu nyquist sysl sys2 sysn nyquist sysl sys2 sysn w re im w nyquist sys 函数为带输出变量引用的函数 可计算系统在频率w处的频率响应输出数据 而不绘制出曲线 其中 输出变量re为频率响应的实部 im为频率响应的虚部 w是频率点 如果系统为离散系统 则相应的nyquist函数为dnyquist 函数命令调用为 re im w dnyquist a b c d Ts re im w dnyquist a b c d Ts iu re im w dnyquist a b c d Ts iu w re im w dnyquist num den Ts re im w dnyquist num den Ts w 6 2 4求连续系统Nichols曲线的函数 求连续系统Nichols曲线的函数为nichols 该函数命令的调用格式为 mag phase w nichols sys nichols sys nichols sys w nichols sysl sys2 sysn nichols sysl sys2 sysn w 6 2 5求系统幅值裕度与相位裕度 求系统幅值裕度与相位裕度的函数为margin 该函数命令的调用格式为 margin sys margin mag phase w Gm Pm Wcp Wcg margin sys Gm Pm Wcp Wcg margin mag phase w margin sys 函数可以从频率响应数据中计算出幅值稳定裕度 相位稳定裕度及其对应的角频率 输入参量sys一般是用系统的开环传递函数描述的系统模型 对于开环SISO系统 既可以是连续时间系统也可是离散时间系统 当不带输出变量引用函数时 margin 函数可在当前图形窗口中绘制出带有稳定裕度的Bode图 margin mag phase w 函数可以在当前窗口中绘制出带有系统的幅值裕度与相位裕度Bode图 其中 mag phase及w分别为由bode或dbode求出的幅值裕度 相位裕度及其对应的角频率 Gm Pm Wcp Wcg margin sys 函数命令可计算出系统的频率性能指标 返回的Gm是系统幅值裕度 不是分贝值 及其对应的角频率Wcp Pm是相位裕度及其对应的角频率Wcg 函数不绘制Bode图 Gm Pm Wcp Wcg margin mag phase w 函数是margin mag phase w 函数带输出变量的引用形式 该函数不绘制Bode图 输出变量返回的Gm是系统幅值裕度及其对应的角频率Wcp Pm是相位裕度及其对应的角频率Wcg 对于离散时间系统 可以先用函数dbode 计算频率响应 然后再调用margin 函数 即用margin mag phase w 函数命令 再求系统性能指标或绘制Bode图 例6 4已知一单位反馈系统前向通道的传递函数为 试绘制出Bode图和nyquist曲线 并计算系统的频域性能指标 解 根据题目要求 给出以下调用函数命令的程序段 num 281282 den 151010510 sys tf num den figure 1 mag phase w bode sys figure 2 nyquist sys figure 3 nichols sys gm pm wcp wcg margin mag phase w 求得的系统的性能指标如下 gm 332 1666 pm 38 6923 wcp 25 8471 wcg 1 2490 6 2 6求系统零极点和根轨迹的函数 绘制系统的零极点图的函数为pzmap 该函数命令的调用格式为 pzmap sys pzmap p z p z pzmap sys 1 求系统零极点函数 pzmap sys 函数可在复平面内绘制出系统的零极点图 sys可以是状态空间形式或传递函数形式 pzmap p z 函数可在复平面内制出的零极点图 其中列矢量p为极点位置 列矢量z为零点位置 这个命令用于直接绘制给定的零极点图 p z pzmap sys 函数可在复平面内绘制出系统的零极点图 并能求出系统所有零点和极点 2 求系统根轨迹函数 求连续系统根轨迹的函数为rlocus 该函数命令的调用格式为 rlocus sys rlocus sys k r k rlocus sys rlocus sys 函数可绘制系统的根轨迹 sys可以是状态空间形式或传递函数形式 rlocus sys k 函数可利用制定增益clocus来绘制系统的根轨迹 r k rlocus sys 函数可绘制系统的根轨迹 同时得到复根轨迹的位置矩阵r及相应的增益矢量k 例6 5设有一高阶系统开环传递函数为 试绘制其零极点图和闭环根轨迹图 解 1 绘制系统的零极点图在命令窗口输入 n1 0 010 2181 4369 359 d1 0 060 2680 6356 271 sys tf n1 d1 pzmap sys p z pzmap sys 计算出系统的三个极点和三个零点 p 5 77100 6522 4 2054i0 6522 4 2054iz 16 5389 2 6305 7 0476i 2 6305 7 0476i 2 绘制系统闭环根轨迹图继续在命令窗口输入 rlocus sys 小结 连续系统 控制系统时域分析MATLAB实现 单位阶跃响应函数dstep 单位冲激响应函数dimpulse 零输入响应函数dinitial 单位阶跃响应函数step 单位冲激响应函数impulse 零输入响应函数initial 离散系统 连续系统 控制系统频域分析MATLAB实现 求Bode图的函数为bode 求系统Nyquist曲线的函数为nyquist 求Nichols曲线的函数为nichols 离散系统 求Bode图的函数为dbode 求系统Nyquist曲线的函数为dnyquist 求Nichols曲线的函数为dnichols 求系统幅值裕度与相位裕度的函数为margin 求系统零极点和根轨迹 求零极点的函数为pzmap 求根轨迹的函数为rlocus 格式 不带输出变量引用的函数 可绘制出曲线 带有输出变量引用的函数 可计算系统Bode图的输出数据 而不绘制出曲线 6 3利用MATLAB语言进行编程仿真 6 3 1面向传递函数的线性系统仿真 1 仿真算法 控制系统最常见的典型闭环系统结构如图所示 由图可知 如果G s 为真分式 则系统的开环状态方程为 式中 经过推导 可得 2 仿真程序框图与实现 3 初始数据 4 程序语句 c c a 1 首一化处理a a a 1 A a 2 n 1 A rot90 rot90 eye n 1 n fliplr A 形成能控标准形A阵B zeros 1 n 1 1 形成输入阵B n维列向量 m1 length c 分子系数向量维数m 1C fliplr c zeros 1 n m1 形成输出阵C n维行向量 Ab A B C V 形成闭环系数阵AbX X0 y 0 t T0 设初值 准备递推运算 N round Tf T0 h 确定输出点数fori 1 N 四阶龙格 库塔法K1 Ab X B R K2 Ab X h K1 2 B R K3 Ab X h K2 2 B R K4 Ab X h K3 B R 求各次斜率KX X h K1 2 K2 2 K3 K4 6 求状态y y C X 求输出并以向量形式保存t t t i h 输出对应时刻以向量形式保存end t y 输出数据形式结果plot t y 输出曲线形式结果 例6 6图所示的单位负反馈系统 系统阶次为n 4 外部输入单位阶跃函数信号r t 1 t 取仿真步长h 0 01 各状态变量初值均为零 仿真系统的动态响应 解 1 准备系统参数及初始值系统特征式阶次 n 4仿真计算步长 h 0 01输入阶跃幅值 R 1传递函数中分子系数 c 2 3 2 8 传递函数中分母系数 a 1 5 4 8 2 反馈系数 V 1状态变量初值 X0 0 0 0 0 起始时间 T0 0终止时间 Tf 10 2 输入上述参数 并运行程序CSS1 在命令窗口键入以下命令语句 n 4 h 0 01 R 1 V 1 T0 0 Tf 10 X0 0 0 0 0 a 15482 c 2328 CSS1 6 3 2面向结构图的线性系统仿真 1 仿真算法 1 仿真系统的典型结构框图 由图可知 根据图中的拓扑连接关系 可逐个写出每个环节输入受哪些环节输出的制约和影响 现列写如下 把环节之间的关系和环节与参考输入的关系分别用矩阵表示出来 成为 即 连接关系矩阵 输入连接矩阵 各环节输入向量 各环节输出向量 从联接矩阵W的元素值直接看出各环节之间的联接情况 即 表示环节j不与环节i相连 表示环节j与环节i有连接关系 表示环节j不与环节i直接相连或通过比例系数相连 表示环节j不与环节i负反馈相连或通过比例系数负反馈相连 表示环节i单位自反馈或通过比例系数自反馈 2 典型环节的确定 这里选用下式作为典型环节 由于每个动态环节都可以用典型环节来表示 故系统中所有环节的输出 输入关系可以统一用矩阵来表示 即 式中 3 系统的求解 联立下面两式 整理得 其中 2 仿真程序框图与实现 3 仿真程序 1 参数输入 系统参数输入方法 由程序自动形成A B C D各矩阵 A diag P 1 B diag P 2 C diag P 3 D diag P 4 联接矩阵输入方法 建立非零元素矩阵WIJ m 3型 将非零元素按照i j wIJ次序逐行输入 其中 i为被作用环节号 第一列 j为作用环节号 第二列 wIJ为作用关系值 联接系数 第三列 而下标为0 表示系统外加参考输入作用 2 程序实现 A diag P 1 B diag P 2 求开环系统系数C diag P 3 D diag P 4 m length WIJ 1 求非零元素个数W0 zeros n 1 W zeros n n 建立初始W W0阵fork 1 mif WIJ k 2 0 W0 WIJ k 1 WIJ k 3 求W0阵elseW WIJ k 1 WIJ k 2 WIJ k 3 求W阵endendQ B D W Qn inv Q 求Q和Q逆阵R C W A V1 C W0 求R V1阵Ab Qn R b1 Qn V1 形成闭环系数阵Y Yt0 y Y nout t T0 置初值 N round Tf T0 h L1 总输出点数fori 1 N 每循环一次 输出一点数据forj 1 L1 两输出点之间计算L1次K1 Ab Y b1 Y0 K2 Ab Y h K1 2 b1 Y0 K3 Ab Y h K2 2 b1 Y0 K4 Ab Y h K3 b1 Y0 Y Y h K1 2 K2 2 K3 K4 6 龙格 库塔法求解微分方程endy y Y nout 保存输出环节动态响应值t t t i h L1 保存时间向量end t y 输出响应数据plot t y 输出响应图形 例6 7如图所示带有电流调节器和速度调节器的双闭环调速系统 外部输入信号为单位阶跃函数 状态变量初值为零 仿真该系统的动态响应 解 1 将系统按典型环节形式编号 共有9个环节 2 确定环节参数矩阵P及连接关系矩阵WIJ 3 确定其它仿真参数 环节个数n 9 阶跃输入幅值Y0 1 各环节初值Yt0 计算步长h 0 001 打印间隔点数L1 10 每隔10个点输出一次 仿真起始时间T0 0 仿真终止时间Tf 0 4 输出环节编号nout 7 4 在命令窗口输入以上数据 运行程序CSS2 得到所示的单位阶跃响应曲线 n 9 Y0 1P 10 0110 00 016410 0867 10 00210 00 0410 0128 10 00167300 10 01282 50 01700 10 0110 10 0210 WIJ 101 211 29 1 321 431 48 1 541 651 67 0 136 761 860 072 970 0067 Yt0 000000000 h 0 001 L1 10 T0 0 Tf 0 4 nout 7css2 习题6 4已知某随动系统的动态结构图如图所示 采用如下的方法仿真该系统在单位阶跃函数作用下的动态和稳态响应 状态变量初值为零 1 利用MATLAB的提供的环节方框图模型的化简函数及其它分析函数 2 调用CSS2程序仿真系统的阶跃响应 1 2 3 4 6 5 6 3 3线性系统的离散相似法仿真 1 仿真算法 1 连续系统状态方程的及其离散化 设连续系统状态方程为 若使用一阶保持器 则离散化状态方程的解 或 上面的式子中 2 典型环节状态方程的离散化 典型环节的传递函数为 略去下标i 并改写为 令 典型环节的连续微分方程为 则 将其离散化可推得 其中 3 常见动态环节的离散系数及离散化状态方程 积分环节的传递函数为 环节离散化系数 这里 比例积分环节的传递函数为 这里 惯性环节 环节离散化系数 比例惯性环节 4 的求取 2 仿真程序框图 3 仿真程序 输入数据 输入系统各环节参数阵P 连接阵WIJ 运行参数n h L1 T0 Tf 阶跃输入幅值Y0 输出环节编号nout 方法同CSS2 A diag P 1 B diag P 2 求开环系统系数C diag P 3 D diag P 4 m length WIJ 1 求非零元素个数W0 zeros n 1 W zeros n n 建立初始W W0阵fork 1 mif WIJ k 2 0 W0 WIJ k 1 WIJ k 3 求W0阵elseW WIJ k 1 WIJ k 2 WIJ k 3 求W阵endend fori 1 nif A i i 0 求积分或积分比例环节各系数FI i 1 FIM i h C i i B i i FIJ i h h C i i B i i 2 FIC i 1 FID i 0 求积分环节系数if D i i 0 若为积分比例环节 修正系数FID i D i i B i i elseendelse 求惯性或惯性比例环节各系数FI i exp h A i i B i i FIM i 1 FI i C i i A i i FIJ i h C i i A i i FIM i B i i A i i FIC i 1 FID i 0 求惯性环节系数if D i i 0 若为惯性比例环节 修正系数FIC i C i i D i i A i i B i i FID i D i i B i i elseendendendY zeros n 1 X Y y 0 Uk zeros n 1 t T0 置初值N round Tf T0 h L1 计算总输出点数 fork 1 Nforl 1 L1Ub Uk 保存前一次输入值Uk W Y W0 Y0 求当前UkUdot Uk Ub h 求当前Uk的导数Uf 2 Uk Ub 求下一步的输入UfX FI X FIM Uk FIJ Udot 求状态向量Y FIC X FID