简介面向机械工程、机器人学及自动化领域学习者的MATLAB/Simulink建模仿真实战资料聚焦一阶倒立摆这一经典力学模型的动态稳定性与控制问题。内容围绕动力学方程构建、Simulink模型搭建及PID反馈控制策略展开涵盖考虑重力、惯性、摩擦力等因素的一阶非线性微分方程建模思路以及设定初始条件与外部扰动后利用Scope观察摆杆角度变化的方法适合需要快速上手倒立摆仿真练习的高校学生与初阶研究者也可作为课程设计与毕业设计的参考资料。压缩包整体约14KB小巧精简便于直接下载研习。目前已有3844人学习下载反映出较高的实用参考价值。从这份资料可掌握基础PID调参到LQR、滑模控制等进阶方案的优化方向为机器人稳定行走等实际应用打下基础。1. 一阶倒立摆的建模仿真与控制为什么它是稳控算法最佳的“试验台”把一根摆杆竖直放在小车上不给外力它一秒钟内就会倒下去。一阶倒立摆要做的就是让小车在轨道上来回移动把摆杆稳稳扶在竖直位置。这个系统结构简单、数学模型不过二阶却因为本质不稳定、强耦合、纯滞后极小成了控制理论课、仿真竞赛和嵌入式开发里最常见的验证对象。无论你是在做课程设计还是在准备研究生数学建模竞赛或者想把手里的 PID、状态反馈算法从仿真搬到真实电机上一阶倒立摆都是成本最低、最能看出控制器“真功夫”的试验台。本文按照建模、控制器设计、闭环仿真、避坑、实物移植的顺序把一阶倒立摆从方程到可运行代码完整走一遍每个环节都会给出参数设定依据和失败时的排查思路。2. 一阶倒立摆的数学建模从运动方程到 MATLAB/Simulink 状态空间模型2.1 物理模型选型与运动方程推导一阶倒立摆的物理模型有两种常见形态一种是经典的小车-摆杆系统小车在导轨上平移电机通过皮带驱动小车控制量是作用在小车上的水平力另一种是旋转臂式电机直接驱动摆杆转轴摆杆在竖直平面内绕轴转动控制量是电机输出的力矩。两者在数学上都叫一阶倒立摆区别在于状态变量和控制量的定义。做建模仿真分析时我一般先用旋转臂式把问题简化因为状态变量少、参数好标定后续移植到实物时也更容易用 STM32 或 Arduino 驱动直流电机和舵机来复现。旋转臂式的建模从摆杆的转动开始。设摆杆质量为 m质心到转轴的距离为 l摆杆绕转轴的转动惯量为 J重力加速度为 g控制力矩为 τ。以摆角 θ 为广义坐标θ0 表示摆杆竖直向上正方向取摆杆向右倾倒。对转轴取力矩平衡得到J·θ m·g·l·sinθ τ这个方程里没有小车位移项因为摆杆只做单轴转动。如果换成小车式方程会变成一组耦合的二阶微分方程还要额外处理小车位移。对于入门级的建模与仿真分析旋转臂式已经足够暴露一阶倒立摆的核心矛盾重力项 m·g·l·sinθ 始终想把摆杆拉倒而控制力矩 τ 必须连续对抗这个趋势。注意 θ0 的定义必须写清楚否则后面线性化和控制器符号很容易整体反号。2.2 线性化与状态空间矩阵符号和单位最容易出错的地方一阶倒立摆的原始方程含有 sinθ是强非线性系统。在摆角较小的范围内即 |θ| 不超过 0.3 rad约 17°时可以用 sinθ≈θ 做线性化。这个近似决定了仿真的适用范围控制器设计完成后如果初始摆角超过 0.3 rad线性化模型算出的增益往往无法让真实非线性系统稳定这也是后面仿真发散的一个重要来源。线性化方程改写为J·θ m·g·l·θ τ定义状态变量 x1θx2θ控制输入 uτ输出 yθ则状态空间方程为x1 x2 x2 (m·g·l / J)·x1 (1 / J)·u写成矩阵形式就是A [0 1; m·g·l/J 0]B [0; 1/J]C [1 0]D 0这里最容易翻车的是转动惯量 J 的计算。很多教材直接给 J m·l²/3这是把摆杆当成质量均匀的细长杆、且转轴在杆端点时的结果。如果你的摆杆是金属杆加末端配重质心位置和转动惯量都要重新算。m 是摆杆总质量l 是质心到转轴的距离不是摆杆全长。单位上全部用国际单位制长度用米质量用千克力矩用牛·米否则算出来的 A 矩阵特征值数量级会离谱到看不出系统是否稳定。2.3 在 MATLAB 和 Simulink 中搭出第一个可运行模型用 MATLAB 实现上述状态空间模型非常直接核心代码如下% 一阶倒立摆参数定义旋转臂式 m 0.05; % 摆杆质量kg l 0.12; % 质心到转轴距离m注意不是摆杆全长 g 9.81; % 重力加速度m/s^2 J m * l^2 / 3; % 转动惯量细长杆绕端点近似kg*m^2 % 线性化状态空间矩阵 A [0 1; m*g*l/J 0]; B [0; 1/J]; C [1 0]; D 0; sys ss(A, B, C, D); % 开环零输入响应初始角度 0.1 rad x0 [0.1; 0]; initial(sys, x0, 2); % 观察 2 秒摆角会快速发散 grid on;代码中 initial() 的作用是让系统在没有任何控制输入的情况下从初始角度 0.1 rad 开始自由演化。因为 A 矩阵存在右半平面极点曲线会呈指数发散这是验证模型符号是否正确的最快方法。如果发散方向反了比如角度先变小再反向增大说明重力项的正负号取错了。ss() 是 MATLAB 构建状态空间模型的函数后面所有控制器设计都基于这个 sys 对象。在 Simulink 里搭建同款模型同样简单。从 Continuous 库拖入 State-Space 模块把 A、B、C、D 四个矩阵填到模块参数里Initial conditions 填 [0.1; 0]输入接 Constant 0输出接 Scope仿真 2 秒就能看到同样的发散曲线。需要注意的是Simulink 模块里的矩阵名如果直接写 A、B、C、D会优先读取基础工作区的变量。为了避免公式和模块参数两处维护我习惯在模型的 InitFcn 回调里写参数计算代码让模型每次启动前自动刷新工作区变量确保改动参数后仿真结果同步更新。3. 控制器设计PID、极点配置与 LQR 的参数计算逻辑3.1 本质不稳定系统的控制需求为什么纯比例控制稳不住拿到线性化模型后最直观的思路是试 PID 控制。先写一个纯比例控制器 u Kp·θ代入闭环系统后特征方程变成 J·s² - m·g·l Kp 0。因为方程里没有一次项闭环极点落在虚轴上系统会持续等幅振荡无论 Kp 取多大都稳定不了。这个结果经常让第一次接触倒立摆的人困惑以为是参数没调好其实问题出在系统本质上缺少阻尼项。一阶倒立摆的 A 矩阵有一个正实部极点和负实部极点属于鞍点型不稳定系统。控制器的核心任务不仅是把正极点拉回左半平面还要给系统注入足够的阻尼。纯比例控制只能改变系统“刚度”相当于把摆杆拉回竖直位置的力但它无法消耗能量所以角度会一直振荡。加入微分项 u Kp·θ Kd·θ 后闭环特征方程才有一次项系统才可能出现实数极点或稳定共轭极点。这就是为什么网上大部分倒立摆 PID 例程强调“PD 控制够用PID 里的积分项反而容易引入振荡”。如果用 PID 控制器我一般先令 Ki0只调 PD 参数。调法跟普通二阶系统类似Kd 初始给 0Kp 从 0 慢慢增大直到角度曲线出现等幅振荡记下临界增益 Kp_crit然后取 Kp 0.6·Kp_critKd 从 0 开始增大直到曲线收敛且超调量小于 20%。这个过程在 Simulink 里可以用 Signal Editor 或手动改 Gain 模块实现但说实话靠肉眼观察 Scope 调参效率不高后面会讲到用极点配置直接算出增益省去大量试凑。3.2 极点配置法用 place() 直接算出反馈增益工程上更可靠的做法是状态反馈。取控制律 u -K·x -k1·θ - k2·θ闭环系统变成 x (A - B·K)·x。这里 K 是一个 1×2 的行向量k1 相当于比例增益k2 相当于微分增益。因为状态反馈利用了完整的系统状态信息它能同时改变闭环极点的实部和虚部控制效果远比单独调 PID 直观。极点配置的核心思路是先决定期望的闭环极点再反解 K。对于二阶系统期望极点通常写成一对共轭复数% 继续使用第 2 节定义的 A、B 矩阵 p [-3 2i, -3 - 2i]; % 期望闭环极点实部为负保证稳定 K place(A, B, p); % 检查闭环极点是否落在期望位置 eig(A - B*K) % 闭环零输入响应 sys_cl ss(A - B*K, B, C, D); x0 [0.1; 0]; initial(sys_cl, x0, 3); grid on;place() 函数会根据期望极点反算出状态反馈增益。代码里极点取 -3 ± 2i实部 -3 决定收敛速度虚部 ±2 决定振荡频率。实部绝对值越大摆杆回正越快但控制力矩也越大执行器更容易饱和。一般取系统开环极点实部的 2 到 5 倍作为参考旋转臂式倒立摆的开环极点约为 ±√(m·g·l/J)代入参数后大约是 ±6.6所以期望极点实部取 -10 到 -20 都不算激进。虚部如果太大控制量会像 PID 微分项一样高频抖动反而触发摩擦和噪声问题。place() 和 acker() 都是 MATLAB 提供的极点配置函数。acker 基于 Ackermann 公式适合二阶这种低阶系统place 使用数值算法在高阶或者期望极点有重根时数值更稳。如果 place 报错说系统不可控优先检查 B 矩阵是否写错尤其是控制输入对应哪个状态。一阶倒立摆的 B 矩阵是 [0; 1/J]只有第二个分量非零表示控制力矩直接作用于角加速度这一点不能记错。3.3 LQR 最优控制用加权矩阵替代人肉试凑极点配置需要人为指定极点而极点和系统性能之间不是线性关系调起来仍然靠经验。LQR 线性二次型调节器把问题转换成优化目标设计 K 使代价函数 J ∫(x·Q·x u·R·u)dt 最小。Q 是状态权重矩阵R 是控制权重它们比极点直观得多也是我在实物倒立摆上首选的控制设计方法。% LQR 控制器设计 Q diag([100, 1]); % 角度权重 100角速度权重 1 R 0.01; % 控制力矩权重 K_lqr lqr(A, B, Q, R); % 查看闭环极点 eig(A - B*K_lqr) % 仿真对比 sys_lqr ss(A - B*K_lqr, B, C, D); x0 [0.1; 0]; initial(sys_lqr, x0, 3); grid on;Q 和 R 的物理含义非常清楚。Q 的第一个对角元素对应角度误差的惩罚设得越大系统越 aggressively 把角度拉回零第二个对角元素对应角速度惩罚设得越大摆动速度就越平缓。R 是控制能量惩罚R 越大控制力矩越小系统响应越慢。常见调法是从 Qdiag([100,1])、R0.01 开始如果控制量超出执行器能力就把 R 调大如果收敛太慢就把 Q 的第一个元素调大。这种调参方式比极点配置少一层抽象也是业界对倒立摆、四旋翼这类欠驱动系统做控制时的常见做法。对比 PD、极点配置、LQR 三种方案PD 适合快速验证和资源受限的嵌入式环境极点配置适合课程设计里要求明确给出极点的情况LQR 最接近工程实际。在仿真阶段我习惯用 LQR 设计出 K 矩阵再回到 Simulink 里验证非线性模型下的表现而不是直接在 Scope 上反复改 PID 参数。4. 闭环仿真验证Simulink 参数设定与结果判定4.1 在 Simulink 中搭建含饱和与扰动的闭环回路模型和控制器都准备好后闭环仿真不能只把两个模块一接就完事必须把执行器饱和、测量噪声和外部扰动放进去否则仿真结果完全没有迁移价值。我在 Simulink 里搭闭环的常用结构是State-Space 模块作为被控对象前面接一个 Gain 模块实现 -K 反馈在反馈回路里加一个 Saturation 模块限制控制力矩上限再在对象输入端用 Pulse Generator 叠加周期性扰动用来模拟外力冲击。搭建步骤可以按四步走第一步拖入 State-Space 模块把线性化模型的 A、B、C、D 填好初始条件设 [0.1; 0]第二步用两个 Gain 模块分别实现 -k1·θ 和 -k2·θ求和后作为控制输入第三步在控制输入通道上加 Saturation上下限都设成 0.5 N·m第四步加入 Pulse Generator振幅设 0.02 N·m周期设 2 秒占空比设 10%这样能模拟一个短暂的冲击力矩。Scope 同时观察角度和控制量曲线就能看到扰动出现后系统如何恢复。如果不想手拉一堆 Gain 模块可以用 MATLAB Function 块写控制律代码量更少function u fcn(theta, theta_dot, K) % 状态反馈控制律 u -K * [theta; theta_dot]; end这里 K 直接作为外部参数传入 MATLAB Function 块修改方便。逻辑上它等价于两个 Gain 模块的求和但代码可读性更好尤其是 K 矩阵含有两个以上状态变量时。参数说明u 是电机力矩K 是 LQR 或 place 计算出的反馈增益theta 是摆角theta_dot 是角速度。注意函数输出 u 会被 Saturation 模块限幅所以即使 K 算得很大实际作用到对象上的力矩也不会超过设定值。4.2 求解器、步长和初始角度的设定原则Simulink 仿真参数设置对倒立摆这种本质不稳定系统影响极大。默认的求解器是 ode45变步长模式对大多数仿真够用但一阶倒立摆在初始阶段发散速度很快如果最大步长不限制误差会指数放大导致曲线看起来像随机震荡。表里是我验证过的保守参数组合参数推荐值原因求解器ode45线性模型无刚性变步长够用最大步长0.001 s限制数值误差防止发散仿真时长5 s足够观察收敛和扰动恢复初始角度0.1 rad约 5.7°在线性化有效范围内控制限幅±0.5 N·m模拟真实电机最大输出扰动幅度0.02 N·m近似随机外力冲击如果出现高频振荡或者“Simulink cannot solve”的报错第一步不是调控制器而是把求解器换成 ode15s最大步长降到 1e-4 再试。ode15s 是刚性求解器适合系统动力学在局部区域变化剧烈的情况。另一个要点是如果后面要移植到 STM32 或 Arduino务必把仿真改成固定步长离散求解器步长选 0.001 s和实物控制周期一致。连续仿真的控制律在离散化后性能会有损失早早在仿真里暴露这个问题能省很多调试时间。4.3 从仿真曲线判定控制器是否合格仿真跑完后不能只看“曲线最终平了”就认为控制器合格。我一般按三个指标判定稳定时间、超调量、稳态误差。稳定时间指角度从初始值进入 ±0.01 rad 误差带后不再超出的时间一阶倒立摆 LQR 控制器通常能在 1.5 秒内完成超调量指角度曲线越过零点后的最大反向幅度超过 20% 说明阻尼偏小稳态误差看曲线最终是否严格趋近于零如果存在偏移说明模型里存在未补偿的常量干扰比如摩擦力或重力偏置。控制量曲线同样需要检查。理想情况下控制力矩在初始阶段会出现一个较大的恢复脉冲随后迅速回落到接近零的小幅波动区间。如果控制量全程贴着限幅最大值走说明期望极点或 LQR 权重设计得过于激进执行器已经饱和控制器实际上失去了调整余地。这种情况下即使仿真曲线看似稳定实物上也极易翻车因为真实电机的响应延迟会让饱和问题更严重。把仿真曲线和实物曲线叠在一起对比这也是判断模型可信度的最终方法。5. 避坑指南一阶倒立摆建模与仿真中最常见的五个翻车点这是全文最值得反复看的部分。以下每条都是我用 MATLAB/Simulink 做倒立摆仿真时实实在在踩过的坑按“现象、原因、解决”写清楚遇到类似问题可以直接对照排查。5.1 现象初始发散曲线方向不对系统越控越不收敛现象initial() 画出的开环响应不是指数发散而是先往回收再反向倒下去或者在仿真一开始就报矩阵维度错误。原因角度零点定义和重力力矩符号不一致。θ0 若定义在竖直向上重力力矩是 m·g·l·sinθ线性化后为正项若在推导时把正方向取反A 矩阵变成 [0 1; -mgl/J 0]闭环极点跑到虚轴附近系统看起来像“稳定”实则没有对抗重力的能力。另一个常见原因是 l 取了摆杆全长而不是质心到转轴距离导致开环极点数值偏大。解决把 θ0 定义为竖直向上正方向为摆杆向右倾倒重力项写成 m·g·l·θ。检查后如果发现符号一直反就在最终控制律 u -K·x 前面加负号相当于把控制器方向反转。最稳妥的办法是先用 initial(sys, x0, 2) 验证开环响应确认角度在 2 秒内增大到 0.3 rad 以上再进入控制器设计。这一步花一分钟能省掉后面所有因为符号错误导致的无效调试。5.2 现象加控制器后曲线高频振荡Scope 上全是毛刺现象角度曲线大致在零附近但波形带有很多高频小锯齿控制量曲线抖动明显感觉像系统在抖。原因反馈回路里没有滤波角速度信号被测量噪声污染。连续仿真虽然没有真实传感器噪声但变步长求解器在快速变化区域的数值噪声会被微分项放大等效于一个高频噪声源。极点配置或 PID 的微分增益越大高频放大越严重。解决在角速度反馈支路加一阶低通滤波器截止频率设在控制器带宽的 5 到 10 倍比如控制器闭环极点实部是 10滤波器截止频率选 100 rad/s 左右。Simulink 里直接用 Transfer Fcn 模块分子填 [1]分母填 [0.01 1]时间常数 0.01 秒。改用 LQR 也能从根源上规避这个问题因为 LQR 的状态反馈不依赖显式微分角速度通过观测器估计时天然带有滤波效果。5.3 现象Simulink 报错提示“无法求解”或仿真步长无限缩小现象仿真运行到某时刻突然报错提示无法自动选择步长或者步长不断缩小到 1e-8 以下仿真卡死。原因系统进入数值刚性区域一般是状态变量超出线性化有效范围后sinθ 的非线性作用被线性模型放大导致反馈控制失效角度迅速撞到物理极限。也可能是控制器输出过大饱和模块设置不当让系统在饱和边界来回切换数值求解器无法收敛。解决先检查初始角度是否保持在 0.1 rad 以内如果偏大就减小再检查 Saturation 上下限是否合理0.5 N·m 对应的小电机能否输出这么大的力矩。把求解器换成 ode15s最大步长手动设 1e-4 秒通常能绕过数值死区。还有一个隐蔽原因是 C 矩阵填成 [0 1] 而不是 [1 0]导致输出是角速度反馈回路正反馈系统发散后自然无法求解。5.4 现象PID 怎么调都稳不住角度总是来回摆现象比例增益从很小调到很大角度始终等幅振荡或者发散Kd 加得再多也没用。原因这是本质不稳定系统的典型表现纯比例控制只能把正极点拉到虚轴无法提供阻尼而 PID 的积分项对倒立摆是灾难性的积分作用会不断累积误差让系统出现低频极限环。很多人的埋点是 PID 里的 I 项默认不为零导致问题雪上加霜。解决先用 PD 控制Ki 设 0Kp 从小到大调找出临界振荡点再把 Kp 降到临界值的 60%Kd 从 0 开始增大。如果调了很久仍不稳定直接放弃手调 PID用 place() 或 lqr() 计算 K 矩阵再返回 PID 等效参数。事实上状态反馈里的 k1 等价于 Kpk2 等价于 Kd用 MATLAB 算一次就够比手工试凑快得多。5.5 现象仿真成功但实物移植后电机乱撞小车反复向两端跑现象仿真曲线很完美角度稳定在零点但同一套增益烧进实物后电机嗡嗡响、小车快速来回冲甚至撞上限位开关。原因仿真模型没有执行器饱和、电机响应延迟、死区摩擦和测量量化误差。一阶倒立摆对控制量极其敏感仿真中的无限带宽力矩在实物里根本不存在真实电机有电流环带宽限制、减速箱背隙粘滞这些非线性因素叠加后仿真里稳定的增益在实物上就变成了极限环振荡。解决回到仿真在控制输入通道加 Saturation 和 Rate Limiter模拟电机最大力矩和力矩变化率在传感器支路加量化模块模拟编码器分辨率把控制周期改成实际中断周期 1 ms用离散求解器重新仿真。仿真里各项指标仍然合格后再把增益和代码一起烧进 STM32 或 Arduino。养成这个习惯后实物调试时间能从几天压缩到半天以内。6. 从仿真到实物状态观测器、摩擦模型与离散化验证6.1 用状态观测器补上角速度信号实物系统通常只有编码器或电位计测量摆角角速度要靠差分或者观测器估计。直接对角度差分会严重放大噪声尤其在低速小角度区域所以工程上普遍用降维状态观测器以角度为输入、估计角速度和可能存在的重力偏置。观测器极点应选在控制器极点的 3 到 5 倍左边保证观测器响应速度远快于控制器但又不会因为太快而放大测量噪声。在 Simulink 里可以先用一个 Observer 模块验证效果再通过 MATLAB Coder 生成 C 代码嵌入控制循环。6.2 摩擦、死区与执行器延时让模型和实物曲线对得上仿真和实物的最大偏差来源就是摩擦。摆杆转轴的库仑摩擦可以建模成 τ_f μ·sign(θ)μ 用空载摆动法实测给摆杆一个初始角度看每次摆动衰减的幅度反推出阻尼系数。电机减速箱的死区等效成控制死区低于某个阈值力矩的输出无效。这些非线性项不需要精确建模只要加入模型里实物追踪曲线的误差就会显著减小。执行器延时通常用一阶惯性环节 1/(T·s1) 近似时间常数 T 从电机铭牌或电流环带宽估算比如带宽 200 Hz 的电机T 大约取 5 ms。6.3 控制周期离散化与 C 代码生成从仿真到实物的最后一步是离散化。连续设计得到的 K 矩阵是针对连续系统的实物控制周期如果是 1 ms必须把控制器离散化后重新验证稳定性。用 MATLAB 的 c2d(sys, 0.001) 将开环状态空间模型离散化再用 dlqr() 设计离散 LQR或者用 place() 在 z 域配置极点。生成 C 代码时MATLAB Coder 会把状态更新方程和控制律直接转成定时器中断里的计算不用手写矩阵乘法减少移植错误。如果你用的是 Arduino也可以直接写一个简单的状态反馈函数把 K 矩阵和上一次状态存成全局变量定时器中断里调一次。我做倒立摆最大的教训就是“仿真越完美实物越容易翻车”。后来养成的习惯是每次仿真的输入通道都加饱和反馈通路上加测量噪声把最坏情况跑一遍才敢写控制器的 C 代码。控制周期固定、扰动有界、输出限幅这三个条件是仿真结果能迁移到实物的前提。还有一个个人习惯每次调完参数都会把 K 矩阵和对应的仿真曲线截图存下来这样修改后续参数时能快速回滚到上一组稳定参数毕竟是多变量调参后悔药不是每次都有。希望这篇文章能帮你在建模仿真和实物控制之间少踩几个坑。本文还有配套的精品资源点击获取
