简介一阶倒立摆系统是动态稳定控制中的经典研究对象兼顾理论深度与工程应用价值。这份资源聚焦其建模仿真与控制实现面向自动控制、机械电子及机器人方向的初学者与进阶者尤其适合需要借助MATLAB/Simulink快速搭建验证环境的人群。压缩包为RAR格式整体仅14KB虽未显示具体文件明细但核心内容围绕Simulink仿真模型展开便于直接融入个人工作区。已有3841人浏览学习在同类控制模型资源中具备一定参考热度资源所涉及的方法体系完整先借助动力学方程考虑重力、惯性、摩擦力等建立模型再在Simulink中搭建仿真结构通过PID控制器进行反馈整定还可拓展到滑模或LQR等高级控制策略。配套的理论梳理可帮助读者系统梳理从方程到仿真再到优化的每个环节在参数调整与扰动测试中理解直立控制的稳定性要求为后续研究或实际应用提供模板。1. 一阶倒立摆为什么说它是控制算法的试金石一阶倒立摆cart-pole在控制领域的位置类似快速排序在算法里的位置结构简单但把开环不稳定、非最小相位、欠驱动这三个让控制工程师最头疼的性质全占了。摆杆角加速度对水平力的响应方向和直觉相反——推小车向右杆顶反而先向左倒这种反直觉特性决定了它必须依赖闭环反馈开环控制无论怎么调都稳不住。这里说的一阶指的是单个摆杆对应串级倒立摆里的一级不是微分方程阶数它的动力学方程本身是二阶的。本文把这套系统从拉格朗日方程推导、Simulink 建模到 PD 与 LQR 控制器实现完整走一遍代码和参数可以直接抄走改。适合正在做倒立摆课程设计、准备机器人自平衡项目或者想找一个标准模型快速验证控制算法的工程师。2. 从拉格朗日方程到 Simulink 可执行模型先把动力学写对2.1 物理模型与坐标约定质量点简化不是偷懒一阶倒立摆的经典模型包含两部分可以在水平方向平动的小车质量 M以及铰接在小车上的摆杆。摆杆简化为质量 m 集中于杆顶的等效质点等效摆长 l 是铰链到质心的距离。这个简化假设在课程设计和仿真阶段完全够用它把转动惯量的计算折进等效摆长里避免在 Simulink 里多传一个 J 参数。坐标约定如下小车位移 x 向右为正摆杆与竖直向上方向夹角 θ 向右倒为正控制量 F 是小车受到的水平外力向右为正。系统的动能和势能代入欧拉-拉格朗日方程得到两组耦合的非线性微分方程(Mm)·x m·l·θ·cosθ - m·l·θ²·sinθ F m·l²·θ m·l·x·cosθ m·g·l·sinθ 0第一行是水平方向的力平衡第二行是摆杆绕铰链点的力矩平衡。两个方程通过 x 和 θ 互相耦合所以不能像处理线性系统那样把两个方程分开求解。θ² 项代表离心力它在摆杆高速摆动时影响显著当 θ 较小时这项是高阶小量线性化时会被忽略。这里要留意符号约定如果 θ 以竖直向下为 0方程形式完全不同控制器极性也要跟着反转。2.2 符号推导与线性化得到状态空间矩阵在竖直向上位置附近令 sinθ≈θ、cosθ≈1忽略 θ² 项上述方程组可以化简为线性形式θ -(Mm)·g·θ / (M·l) - F / (M·l) x m·g·θ / M F / M取状态向量[x, θ, x, θ]控制输入 uF输出为全部状态得到标准状态空间模型。这份线性模型是后续 LQR 设计的基础也是判断系统本质特性的依据。用 MATLAB 写出 A、B 矩阵% 一阶倒立摆线性化状态空间模型 M 1.0; % 小车质量 kg m 0.1; % 摆杆等效质量 kg l 0.5; % 铰链到质心距离 m g 9.81; % 重力加速度 m/s^2 A [0 0 1 0; 0 0 0 1; 0 m*g/M 0 0; 0 -(Mm)*g/(M*l) 0 0]; B [0; 0; 1/M; -1/(M*l)]; C eye(4); D zeros(4,1); sys ss(A,B,C,D); eig(A) % 查看开环极点参数说明A(4,2) 是摆角 θ 到角加速度 θ 的自激项数值约 -21.6说明 θ 偏差越大倒下的加速度越大B(4) 为负对应向右推车、杆顶先向左倒的非最小相位特性。开环极点分布在实轴两侧其中一个极点约 4.65对应时间常数约 0.2 秒意味着控制器必须在 200ms 内对偏差作出反应否则摆杆彻底倒下。2.3 Simulink 模型搭建用 MATLAB Function 封装非线性动力学Simulink 搭建倒立摆模型有两条路线一是直接用线性状态空间模块State-Space填 A、B 矩阵适合验证 LQR二是按非线性方程搭完整模型适合观察大角度行为和测试鲁棒控制。我一般先用非线性模型因为后续对比 PD 和 LQR 时线性模型的仿真结果会过于完美掩盖控制器的真实表现。非线性模型的搭建思路是从 θ 和 x 的解析式出发用积分器把角加速度积分成角速度和角度把加速度积分成速度和位移。首先要做一步代数化简把两个耦合方程解成显式表达式θ ( -(Mm)·g·sinθ - m·l·θ²·sinθ·cosθ - F·cosθ ) / ( l·(M m·sin²θ) ) x ( F - m·l·θ²·sinθ - m·l·θ·cosθ ) / (Mm)这里把 x 表达成 θ 的函数所以两个输出块之间有耦合关系。推荐用 MATLAB Function 块一次算出两个加速度避免 Simulink 因模块间即时反馈产生代数环function [theta_ddot, x_ddot] cartpole_dynamics(F, theta, theta_dot) % 一阶倒立摆非线性动力学 M 1.0; m 0.1; l 0.5; g 9.81; den l * (M m * sin(theta)^2); theta_ddot ( -(Mm)*g*sin(theta) ... - m*l*theta_dot^2*sin(theta)*cos(theta) ... - F*cos(theta) ) / den; x_ddot ( F - m*l*theta_dot^2*sin(theta) ... - m*l*theta_ddot*cos(theta) ) / (Mm); end逻辑说明den 是化简后的分母永远不会为 0M 和 m 都为正sin²θ 有界所以模型在全角度范围都不会出现奇异点这也是很多教材版本里没有说明的细节。x_ddot 的表达式里直接引用 theta_ddot这是两个方程解耦后的代数关系不需要额外状态。MATLAB Function 块的输入端口按 F、theta、theta_dot 顺序连接输出端口连到两组积分器即可。模型内部各模块配置可以参考以下清单模块库位置关键设置MATLAB FunctionUser-Defined Functions输入 F、theta、theta_dot输出 theta_ddot、x_ddotIntegrator两组Continuoustheta 初始值设为 0.15模拟初始倾斜ScopeSinks端口数设为 2分别观察 theta 和 xSaturationCommon Blocks限幅设为 ±10 N模拟执行器饱和PID ControllerContinuous预留后续接控制器输出2.4 线性参考模型给控制器设计用的另一份图纸非线性模型装进 Simulink 之后建议在旁边再放一个 State-Space 模块填写 2.2 的 A、B 矩阵。两套模型并行跑控制器设计阶段用线性模型计算增益验证阶段用非线性模型做最终确认。两者的差异主要来自被忽略的 cosθ 非线性项和 θ² 离心项当摆角超过 30° 时线性模型会明显偏离真实响应这也是判断控制器鲁棒性的一个指标。3. PD 与 LQR 控制策略Simulink 里的实现、整定与对比3.1 为什么先试 PD 而不是 PID积分项在直立控制里的副作用倒立摆直立控制的直觉方案是 PD用 θ 作为比例输入θ 作为微分输入控制律写成u -Kp·θ - Kd·θ比例项提供纠正当前角度偏差的力矩微分项提供阻止角速度继续增加的阻尼。为什么传统 PID 里的积分项在这里很少用因为倒立摆的重力矩在小角度下近似线性mg·l·sinθ ≈ mg·l·θPD 控制的静差可以通过比例项直接克服加入积分项反而引入相位滞后系统更容易振荡。如果你的被控对象换成伺服电机位置环I 项仍然有用但那是另一套整定逻辑。在 Simulink 里实现 PD 有两种方式直接用两个 Gain 块把 θ 和 θ 加权相加再取负或者用 PID Controller 块把 I 设为 0D 设为所需微分系数滤波系数 N 取 200 左右。第二种方式的好处是自带微分滤波可以抑制高频噪声放大推荐优先使用。3.2 PD 参数整定从临界增益到一组可用初值PD 参数如果从零开始试效率太低。我在仿真里常用的套路是先忽略小车动力学只把 θ 和 θ 接成闭环逐渐增大 Kp 直到系统出现等幅振荡记下临界增益 Ku 和振荡周期 Tu。然后按经验值取 Kp ≈ 0.6·KuKd ≈ Kp·Tu/8。以 2.1 节的参数为例一组能直接跑通的初值是 Kp80Kd12。% PD 控制器参数计算脚本 Kp 80; Kd 12; % 反馈回路 u -Kp*theta - Kd*theta_dot;这里的符号约定很关键u 是施加在小车上的水平力正向定义为向右。当 θ 为正向右倒时控制力必须为负向左推所以两个增益都带负号。如果你的模型里 θ 的方向定义反了控制器输出的极性反了仿真会立刻发散这不是参数问题是符号问题。3.3 LQR 状态反馈Q、R 矩阵的物理含义与代码实现LQR 的标准写法是设计一个状态反馈u -K·x使得代价函数最小。它和 PD 的本质区别在于PD 只用了 θ 和 θ 两个状态LQR 把 x、θ、x、θ 四个状态全部纳入反馈。代价函数里 Q 矩阵的对角线分别惩罚四个状态偏离零的程度R 惩罚控制能量。用 MATLAB 计算增益矩阵% LQR 控制器增益计算 Q diag([10, 100, 1, 10]); % 状态权重: x, theta, x_dot, theta_dot R 0.1; % 控制力权重 K_lqr lqr(A, B, Q, R); % 输出增益矩阵, 用于填入 Simulink 的 Gain 块 disp(mat2str(K_lqr, 4));逻辑说明lqr 函数自动求解代数黎卡提方程AP PA - PBR⁻¹BP Q 0不需要手写迭代过程。Q 四个元素分别对应 [x, θ, x, θ] 的权重。把 θ 对应权重调到 100是因为摆角是首要控制目标x 权重 10 意味着小车允许在有限范围内移动但不希望它无限漂移速度项权重小是为了避免对测量噪声敏感。R 取 0.1 表示我们允许控制器使用最大约 3N 的控制力如果执行器推力受限R 需要相应调大。把 K_lqr 复制到 Simulink 的 Gain 块时注意状态顺序必须与模型的信号排列一致Gain 块输入依次是 x、θ、x、θ输出直接连到 Saturation 前的求和点。LQR 反馈和 PD 反馈可以做成一个 Switch 切换方便对照组仿真。3.4 PD 与 LQR 的执行器需求对比两种控制策略仿真对比时除了看摆角是否收敛还要关注控制器输出幅值。PD 的 Kp 从 80 开始试时初始摆角 0.15 rad 对应的输出约 12N会触发 ±10N 的饱和限幅LQR 因为显式惩罚控制力输出通常在 5N 以内这对实际电机选型是有意义的参考。对比项PDLQR反馈状态θ、θx、θ、x、θ小车位置约束不控制会漂移Q 中 x 权重控制回零参数个数2Q、R 共 5 个权重整定方法临界比例度法Bryson 法则 试凑抗模型误差能力中依赖精确模型典型输出幅值易饱和可通过 R 约束3.5 滑模控制如何在这个模型上验证如果想把研究深度往上推一层滑模控制是经典选择。设计滑模面s θ λ·θλ0控制律中包含符号函数项η·sign(s)用不连续项抵消参数不确定性。在 Simulink 里用 MATLAB Function 实现时sign 函数会导致仿真步长变小求解器切换成 ode45 后持续推进速度变慢这个现象本身就是离散仿真的真实反馈。注意仿真里滑模的抖振不一定会造成实际硬件上的高频振动因为理想符号函数在离散求解器下等效于极限环所以仿真结果要谨慎外推到实验台。4. 扰动注入与仿真参数整定让模型暴露真实问题4.1 三类典型扰动的 Simulink 实现只给初始角度偏差验证的是基本稳定性真正检验控制器性能的是外部扰动。以下三类扰动覆盖了大部分课程设计和工程预研需求扰动类型Simulink 实现关注指标脉冲力Step 块从 0 跳到 10N0.1s 后跳回 0加到 F 信号前最大摆角、恢复时间白噪声Band-Limited White Noise噪声功率设 0.01稳态摆角方差参数失配将 MATLAB Function 里的 M 改为 1.2 倍是否仍能稳定脉冲扰动最直观。在 F 信号路径上并联一个 Step 块通过加法器混入控制力Simulink 的 Signal Editor 也可以实现更复杂的扰动序列但 Step 块更轻量。白噪声注入要用 Band-Limited White Noise 而不是 Uniform Random Number前者的功率谱密度参数可以直接对应实际传感器噪声水平。4.2 求解器与仿真步长选择倒立摆系统的自然频率约 4.6 rad/s加上控制器之后闭环带宽通常在 2~5 Hz仿真步长对这个量级非常宽容。我一般分两档处理快速验证控制器逻辑用 ode45 可变步长最大步长限制在 1e-2批量对比实验用 ode4 固定步长 1e-3保证不同控制参数下的数值积分误差一致。% 批量运行 Simulink 模型并读取摆角数据 model cartpole_sim; load_system(model); Q_cases {diag([10 100 1 10]), diag([50 200 2 20])}; R 0.1; for i 1:length(Q_cases) K lqr(A, B, Q_cases{i}, R); set_param([model /LQR_Gain], Gain, mat2str(K)); simOut sim(model, StopTime, 5); theta simOut.logsout.get(theta).Values.Data; fprintf(Q%d 最大摆角: %.4f rad\n, i, max(abs(theta))); end逻辑说明set_param 把计算好的 LQR 增益写回 Simulink 的 Gain 块mat2str 负责把行向量转成 Gain 块能识别的字符串格式。logsout 是模型里为 theta 信号勾选了日志记录才能读取如果模型里没开日志可以用simOut.get(theta)替代具体的变量名取决于 Scope 或 Outport 的命名。批量跑完只看最大摆角不够建议同时记录控制力积分评估能量消耗。4.3 从 Scope 曲线识别的假稳定现象仿真跑出来的曲线看着稳定不代表模型真正收敛。以下三种现象在倒立摆仿真里反复出现现象可能原因检查方式摆角高频小幅振颤D 项增益过大或固定步长太大Kd 减半观察步长改 1e-4小车单方向持续加速LQR 的 x 权重为 0位置不可观测Q(1,1) 设为 10 以上控制器输出长时间顶在限幅值初始偏移过大或 R 太小观察 Saturation 输出曲线还有一种隐蔽情况Scope 里摆角看起来在零点附近但小车已经漂移到了几百米外。因为 Scope 默认只显示最近一段仿真时间或者 y 轴范围没设置位移信号被压缩成一条直线。检查方法很简单把 Scope 的 y 轴范围固定到 [-1, 1]再跑一遍就暴露了。4.4 从仿真结果反推动力学模型的准确性仿真和实验台的偏差通常来自三个被简化掉的环节电机电流环延迟、编码器量化噪声、静摩擦力。Simulink 里用 Transport Delay 模拟传感器延迟用 Quantizer 模拟编码器分辨率用 Dead Zone Coulomb Friction 模拟摩擦加到模型里之后你会发现原来能稳的 LQR 增益开始出现极限环。这一步建议在所有课程设计里都做因为它把MATLAB 里能用推向实际硬件上能用的关键一步。5. 调参技巧与仿真到实验前的最后检查清单5.1 Bryson 法则快速初始化 Q、RLQR 的 Q 矩阵初始值不用拍脑袋。Bryson 法则的核心思想是每个状态或输入的权重取它的容许最大值的平方倒数。例如摆角容许 ±0.15 rad控制力容许 ±5N则 θ 的权重约 44R 约 0.04以此作为起点再微调。状态/输入容许上限权重计算初始值x0.10 m1/0.01100θ0.15 rad1/0.022544x0.5 m/s1/0.254θ1.0 rad/s1/1.01F5 N1/250.04代入后 Qdiag([100, 44, 4, 1])R0.04。跑完如果发现摆角收敛太慢把 Q 中 θ 对应元素倍增如果控制力超过执行器极限增大 R。每次只改一个参数这套流程比随手试数要高效得多。5.2 反馈信号的角度回绕处理积分器直接输出的 θ 会在多次穿越 ±π 之后累积到很大数值当 θ 从 179° 跳到 -179° 时反馈信号产生一个 358° 的假跳变微分项会瞬间输出一个巨大的脉冲。在反馈回路进入 PID Controller 或 Gain 之前加一个角度归一化函数function theta_wrap wrap_angle(theta) theta_wrap atan2(sin(theta), cos(theta)); end逻辑说明atan2(sin, cos) 把任意角度映射回 [-π, π] 区间跳变处变成连续过渡。注意这个函数放在反馈路径不要放在动力学计算里否则模型内部的 θ 和显示用的 θ 不一致Scope 曲线会让人困惑。5.3 能控性检查与延迟预算模型建好后先跑一行代码验证能控性LQR 的前提条件不满足时算出来的 K 毫无意义rank(ctrb(A, B)) % 满秩4, 才能继续设计 LQR传感器延迟对倒立摆的影响比想象中严重。用 Transport Delay 把 θ 反馈延迟 2msLQR 的相位裕度明显下降延迟加到 10ms 时即使调高 θd 对应权重也压不住振荡。这给实验台的硬件选型一个明确约束传感器采样周期和通信延迟必须控制在 10ms 以内否则再好的控制器设计都白搭。先把这层延迟在仿真里压住再回到实验台去调编码器滤波和电机驱动器参数能省下大量现场排查时间。把非线性模型、线性模型、PD 和 LQR 四套东西放在同一个 Simulink 工程里靠 Switch 切换批量跑对比脚本是我在倒立摆项目上一直沿用的工作流。这个流程的产物不只是一份仿真而是一套可以迁移到两轮自平衡小车、单级倒立摆实验台和欠驱动机械臂的参数整定框架。本文还有配套的精品资源点击获取