简介一份面向控制理论与自动化专业学生、研究者的MATLAB/Simulink仿真源码用于一级倒立摆系统的建模与控制器设计可直接应用于相关课程设计、毕业设计或算法验证。资源共3个.m脚本压缩包仅3KB包括主仿真程序、LMI控制器求解脚本与参数定义文件可覆盖倒立摆动力学建模、状态反馈镇定、H∞鲁棒控制及线性矩阵不等式求解等典型环节代码结构紧凑便于逐行研读与修改。已有5697人学习下载适合正在学习现代控制理论、需要快速上手Simulink仿真的读者。通过运行代码可直观观察摆杆角度响应与控制效果并依据LMI条件调整控制器参数理解从系统模型建立、控制器设计到闭环仿真的完整流程同时为深入探索更复杂倒立摆控制问题提供可扩展的基础模板。1. 一级倒立摆为什么是 LMI 控制器的试金石一级倒立摆的开环模型里有一个右半平面极点线性化之后这一项直接进入 A 矩阵任何一点初始偏差都会被重力放大。H∞ 控制要处理的就是在这种开环不稳定前提下把扰动到输出的能量增益压到指定阈值以下而 LMI 是解决这类带性能约束的凸优化问题的标准工具。相比直接极点配置LMI 能同时把稳定性、H∞ 扰动抑制水平和控制能量上限写进同一组矩阵不等式求解一次就得到控制器。用 Simulink 做闭环仿真能直观看到积分器链路、执行器限幅和信号记录对理论控制器的影响。这套资源里的 Parameters.m、LMID.m、LMIC.m 正好对应建模、LMI 求解、控制器生成三步适合用来做控制理论验证也适合入门鲁棒控制的人反复修改参数找感觉。2. 从动力学方程到 Simulink 可用的状态空间模型2.1 摆杆运动方程与线性化一级倒立摆通常被写成小车与摆杆耦合的二阶方程组。资源描述里给出的方程是摆杆子系统的简化写法直接用于 Simulink 积分器搭建时会漏掉车体与摆杆之间的惯性耦合项仿真出来的响应和实际物理过程对不上。实际搭建仿真时更常用下面这组线性化方程(M m) * xddot m * L * thetaddot u m * L * xddot (J m * L^2) * thetaddot m * g * L * theta其中 theta 是摆杆与竖直向上方向的夹角u 是作用在小车上的水平控制力M 是底座质量m 是摆杆质量L 是摆杆质心到转轴的距离J 是摆杆绕质心轴的转动惯量。这组方程假设 theta 工作在小角度范围sin(theta) 约等于 thetacos(theta) 约等于 1同时暂时忽略摩擦和非线性阻尼。第二个方程右侧的 mgL*theta 来自重力在摆杆偏离竖直位置后产生的切向力矩符号为正意味着它促使角度继续增大这就是倒立摆开环不稳定的物理来源。把两个方程写成矩阵形式加速度项是耦合在一起的[Mm m*L ] [xddot] [u ] [m*L Jm*L^2] [thetaddot] [m*g*L*theta]这个 2x2 系数矩阵通常被命名为 Dmat在搭 Simulink 模型之前先求一次逆把 [xddot; thetaddot] 显式解出来。后续只需要四个积分器做两次积分分别得到速度和位置模型链路会非常清晰也方便后面插入扰动信号。2.2 状态变量选择与矩阵形式状态向量取 x [x; xdot; theta; thetadot]控制输入是 u。写成状态空间形式就是 xdot Ax Bu。注意 A 矩阵不是随意填数它由 Dmat 的逆乘以耦合矩阵得到。下面的代码可以直接放在 Parameters.m 末尾也可以独立存成一个 setup 脚本在运行 LMID.m 之前执行。% Parameters.m 运行后得到 M, m, g, L, J % 构造状态空间矩阵供 LMID.m 和 Simulink 使用 Dmat [M m, m * L; m * L, J m * L^2]; invD inv(Dmat); A zeros(4, 4); A(1, 2) 1; % xdot 到 x 的积分链 A(2, 3) invD(1, 2) * m * g * L; % theta 通过耦合项驱动 xddot A(3, 4) 1; % thetadot 到 theta 的积分链 A(4, 3) invD(2, 2) * m * g * L; % theta 直接进入 thetaddot B [0; invD(1, 1); % 控制力 u 对 xddot 的通道 0; invD(2, 1)]; % 控制力 u 对 thetaddot 的通道 C [1 0 0 0; 0 0 1 0]; % 可测输出小车位置和摆角 D zeros(2, 1);这里 A(2,3) 和 A(4,3) 都是正的说明角度项不是阻尼而是驱动项系统的特征值中会有一个右半平面极点这正是倒立摆开环不稳定的体现。B 矩阵第二行来自 invD(1,1)大致相当于等效质量的倒数第四行来自 invD(2,1)表示同一个控制力对摆角加速度的交叉影响。C 矩阵设置为位置和摆角两路输出是为了仿真时直接观察物理量但状态反馈控制器还需要 xdot 和 thetadot所以 Simulink 模型里要把状态从积分器中引出来而不是只使用 C 的输出这一点在第 4 章会再次强调。2.3 Parameters.m 中的参数约定与典型值Parameters.m 是整个流程的入口必须先运行后面的 LMID.m 和 LMIC.m 才能拿到物理参数。文件本身只是赋值语句但参数单位不一致会直接导致控制器设计失败。最容易出错的是 L 到底用摆杆全长还是质心距离以及 J 是绕质心还是绕转轴。下面这组参数是我搭课程模型时常用的一组可以直接替代默认值。参数物理含义典型值说明M小车质量1.0 kg含底座、导轨折算质量m摆杆质量0.2 kg含连接件折算g重力加速度9.81 m/s^2固定环境L摆杆质心到转轴距离0.4 m不是摆杆全长J摆杆绕质心转动惯量0.006 kg*m^2均匀摆杆可近似 m*L^2/3这些数值不是资源里写死的实际实验台上通常需要标定。如果 LMI 求解出的 K 在 Simulink 里表现为发散优先检查 L 是否误用了全长J 是否错用了绕转轴的值。转动惯量差两倍A(4,3) 就差两倍控制器的工作频率会明显偏移。Parameters.m 除了物理参数一般还会预留一组加权变量给后续 H∞ 性能指标使用。常见的是 rho 和 Wzrho 表示扰动通道的权重Wz 表示被调输出的权重。如果原始文件里没有这些变量在 LMID.m 里自己补几行即可不影响参数主流程。3. LMI 控制器设计从 Lyapunov 不等式到 LMID.m / LMIC.m3.1 为什么状态反馈会写成双线性不等式设计目标是找状态反馈 u Kx使得闭环系统 xdot (A B2K)x 渐近稳定同时从扰动 w 到被调输出 z 的 H∞ 范数小于一个预设标量 gamma。根据有界实引理闭环系统满足 H∞ 性能 gamma 的充分必要条件是存在对称正定矩阵 P使得一个分块矩阵负定。这个矩阵里同时出现 PK 和 K*P形成双线性项而 LMI 求解器只能处理线性的矩阵不等式的可行性问题。解决的思路是引入变量替换。令 X P^(-1)对原来的不等式做合同变换再定义新变量 W KX原来带 PK 的双线性项就变成了 AX XA B2*W W*B2所有项对 X 和 W 都是线性的。解出 X 和 W 之后再用 K W * X^(-1) 恢复实际反馈增益。LMID.m 里的核心变量就是 X 和 W理解了这一步整个文件的可读性会提升很多后续在 LMI 里增加极点区域约束或控制能量约束也只是在相同框架上补块。3.2 LMID.m 中的核心求解段拿到资源时可以先看 LMID.m 用的是什么求解框架。如果看到setlmis([])、lmivar、lmiterm三件套那是 LMI Control Toolbox 的写法。它的优点是无需安装第三方求解器缺点是lmiterm里矩阵块的行列位置很容易写错尤其是带s的对称项位置不对会在运行时报 The LMI is not feasible 或找不到决策变量。现在更常见的做法是用 YALMIP 重写同样的 LMI逻辑完全等价排查问题也更方便。下面这段可以直接替换 LMID.m 里原来的求解段。% LMID.m 核心段状态反馈 H∞ 控制器 % 前置条件Parameters.m 已运行A, B 已构造 B1 [0; 0; 0; 1]; % 扰动通道作用在摆角加速度上 B2 B; % 控制通道 C1 [1 0 0 0; % 被调输出 1小车位置 0 0 1 0; % 被调输出 2摆角 0 0 0 0]; % 被调输出 3占位配合控制惩罚 D12 [0; 0; 0.1]; % 控制信号惩罚权重 X sdpvar(4, 4, symmetric); % Lyapunov 矩阵的逆 W sdpvar(1, 4); % 变量替换后的增益 gamma sdpvar(1); % H∞ 性能水平作为优化目标 M [A*X X*A B2*W W*B2, B1, (C1*X D12*W); B1, -gamma, zeros(1, 3); C1*X D12*W, zeros(3, 1), -gamma*eye(3)]; LMIs [X 0, M 0]; optimize(LMIs, gamma, sdpsettings(solver, sedumi, verbose, 0)); K value(W) / value(X); gamma_opt value(gamma);代码里 C1 和 D12 的维度需要严格匹配。C1 是 3x4前三行对应位置、摆角和一个占位输出D12 是 3x1第三行的 0.1 会把控制力的能量计入 H∞ 指标。D12 的数值越大求解器倾向于给出增益更小的 K控制力更平缓但抗扰能力会变差。gamma 是优化变量optimize 的第二个参数把它作为最小化目标最终 gamma_opt 就是当前加权下能达到的最小扰动抑制比。B1 表示扰动进入状态方程的位置这里放在摆角加速度通道模拟外力直接推动摆杆如果扰动来自小车底座可以改成 B1 [0;1;0;0]设计结果会不一样。3.3 LMIC.m 的职责生成控制器并校验闭环LMIC.m 放在 LMID.m 之后执行主要工作是把上一步得到的 K 写回工作区并校验闭环系统是否真的稳定。调试时这一步非常重要因为 LMI 求解器偶尔会在数值边界上返回一个看似可行、实际上闭环有极点在虚轴附近的解。% LMIC.m 生成控制器并校验闭环极点 Acl A B2 * K; eig_cl eig(Acl); if any(real(eig_cl) 0) error(闭环系统仍有右半平面极点请检查参数或放宽 gamma); end fprintf(design gamma: %.4f\n, gamma_opt); fprintf(closed-loop poles: ); fprintf(%.3f , sort(real(eig_cl))); fprintf(\n); save(controller.mat, K, gamma_opt, Acl);这里的 Acl 是状态反馈后的闭环状态矩阵Acl 的特征值实部必须全部为负。打印出的实部越负响应越快但控制力峰值也会更大这是因为反馈增益 K 的模随着期望极点左移而增大。保存到 controller.mat 之后Simulink 可以直接从基工作区读取 K避免每次打开模型重新执行 LMI 求解。如果资源里还涉及离散化LMIC.m 里会先用 c2d 把连续控制器转成离散增益再生成 Simulink 可用的离散状态空间一级倒立摆课程设计通常默认状态连续可测K 是常数增益就够用了。4. Simulink 仿真模型搭建与联合调用4.1 用积分器链路搭被控对象而不是直接用 State-Space 模块很多教程直接拖一个 State-Space 模块把 A、B、C、D 填进去就运行。缺点很明显State-Space 模块默认只输出 y Cx Du如果 C 是 2x4 矩阵那 xdot 和 thetadot 这两个状态就没有外部端口状态反馈控制器拿不到完整状态。虽然模块参数里可以勾选 state 输出端口但用四个积分器组成的积分链路对调试更友好每一路信号都能拉到 Scope 上观察故障定位快很多。链路结构是加速度计算模块输出 xddot 和 thetaddotxddot 经过两个积分器得到 xdot 和 xthetaddot 经过两个积分器得到 thetadot 和 theta。四个状态按 [x; xdot; theta; thetadot] 的次序用 Mux 合成一个四维向量送到 Gain 模块乘以 K得到控制力 u。u 同时反馈回加速度计算模块。加速度计算用 MATLAB Function 模块实现代码如下。function [xdd, thdd] accel(xd, th, thd, u, invD, m, g, L) % 输入xd 小车速度th 摆角thd 摆角速度u 控制力 % 输出小车加速度 xdd摆角角加速度 thdd acc invD * [u; m * g * L * th]; xdd acc(1); thdd acc(2); % 说明thd 在线性化模型中没有阻尼项 % 保留输入端口是为了后续加摩擦或风阻时不用改接口这里把 invD 作为外部参数传入而不是每次计算时重新 inv(Dmat)能减少每个仿真步的矩阵求逆开销。thd 输入虽然没参与运算但保留端口可以在后期加入阻尼项时不用再改模块接口。积分器的初始条件设置为x 初始 0xdot 初始 0theta 初始 0.1thetadot 初始 0。theta 初始 0.1 rad 可以测试控制器在非平衡点能不能把摆杆拉回竖直位置这是判断控制器有效性的最低要求。4.2 工作区变量与 Simulink 参数的绑定Parameters.m、LMID.m、LMIC.m 依次运行后基础工作区里已经有 M、m、g、L、J、invD、A、B、K 这些变量。Simulink 模型里 Gain 模块的参数填 KMATLAB Function 的参数表填 invD、m、g、L积分器初始值可以直接写数字。仿真开始时会自动从基础工作区读取这些变量但有一个常见问题如果模型设置了 Model Workspace 覆盖模型里的变量来自 .mat 文件而 LMIC.m 生成的 K 在基础工作区两边会不一致。最常见做法是把 Parameters.m、LMID.m、LMIC.m 三个调用放到模型的 InitFcn 回调中。在模型属性对话框的 Callbacks 标签页找到 InitFcn填入Parameters; LMID; LMIC;。这样每次启动仿真都会自动按顺序执行不需要手动运行脚本。如果想在命令行里批量修改模块参数可以用 set_param这里给出一个示例。% 从命令行设置 Simulink 模块参数 mdl invpendulum_lmi; open_system(mdl); set_param([mdl /Gain], Gain, K); set_param([mdl /Integrator1], InitialCondition, 0); set_param([mdl /Integrator3], InitialCondition, 0.1);这段代码主要解决两个问题一是在校准 K 时不用每次点开 Gain 模块手动粘贴矩阵二是可以写进循环里做批量调参实验。模块名必须与模型里的实际显示名称完全一致包括大小写。Integrator3 对应 theta 状态的积分器0.1 表示仿真开始时摆角已经偏了 0.1 rad这个值也是后续验证控制器鲁棒性的基准输入。4.3 运行仿真与结果记录模型搭好并启用信号记录后用 sim 命令回到 MATLAB 工作区比手动点 Run 更容易复现结果。假设在模型里给 theta、x、u 三条信号线分别启用了 Logging仿真结束后数据会挂在 out.logsout 下。% 运行 10 秒仿真并取回记录信号 out sim(mdl, StopTime, 10); t out.tout; theta out.logsout.get(theta).Values.Data; x out.logsout.get(x).Values.Data; u out.logsout.get(u).Values.Data;out.tout 是仿真时间序列。get 方法的字符串参数是模型里信号的名字必须一致。theta 是摆角轨迹x 是位置轨迹u 是控制力轨迹。取回数据后可以先画在一张图里确认摆角在 2 秒内回到 0 附近控制力峰值没有长时间顶在饱和限幅上。如果仿真直接发散大概率不是 LMI 求解问题而是反馈增益方向反了把 K 改成 -K 再试通常能快速定位。5. 调参与抗扰验证把 LMI 控制器推到边界5.1 从 gamma 看鲁棒性与控制力之间的取舍LMID.m 解出的 gamma_opt 是当前加权矩阵下能达到的最小扰动抑制水平但这个数值不是越小越好。gamma 越小反馈增益 K 的模越大仿真中的控制力越容易碰到饱和模块产生极限环。典型参数下gamma_opt 压到 1 以下时K 的角度通道会超过 500.1 rad 的初始偏差就能产生 5N 以上的力遇到执行器上限后摆角会出现持续小幅振荡。遇到这种情况把 D12 从 0.1 改到 0.3重新求解gamma_opt 会变大但控制力会明显平缓。调 D12 是 LMI 参数里最直接影响大闭环行为的旋钮。5.2 扰动注入与实测 H∞ 增益为了验证设计指标在 Simulink 的 u 信号汇合点叠加一个扰动 w。用 Signal Builder 或 Step 模块设置在 2 秒时阶跃到 2N持续 0.5 秒后恢复。扰动信号和被调输出 z 都记录到 logsout然后按下式计算实际扰动到输出的能量增益。% 计算实际扰动到输出的能量增益 Ew trapz(t, w.^2); Ez trapz(t, z.^2); gamma_actual sqrt(Ez / Ew); fprintf(design gamma: %.4f\n, gamma_opt); fprintf(actual gamma: %.4f\n, gamma_actual);trapz 是梯形积分用来算离散时间序列的能量。如果 gamma_actual 明显大于 gamma_opt说明仿真中引入了模型里没有的饱和、摩擦或积分误差如果小很多说明加权矩阵留了余量可以继续减小 D12 或增大 C1 中位置和角度的权重。实际增益小于等于设计值的含义是扰动到输出的能量放大倍数被 LMI 真实约束住了。5.3 三个容易被忽视的坑第一J 必须绕质心不要用绕转轴的惯量。用错之后 LMI 求解照样返回 K但闭环极点在虚轴附近抗扰曲线会有很长的拖尾。第二不要在反馈通路里直接形成代数环。加速度计算模块使用当前时刻的 u而 u 来自状态反馈中间若没有 Memory 或 Unit DelaySimulink 会警告代数环gamma_actual 会比设计值偏大。在反馈路径上加 Unit Delay 是最快的解决办法代价是控制器变成一步滞后。第三记录信号的名字不要带空格和中文某些 MATLAB 版本下 logsout 解析带空格的信号名会返回空数组。调试顺序建议从固定 D12 开始用二分法在 0.5 到 0.05 之间搜索每次重跑 LMID.m 并记录 gamma_opt 和实际响应这条曲线能直接告诉你权重与鲁棒性的折中关系。本文还有配套的精品资源点击获取