简介一阶倒立摆是控制理论中的经典动力学模型本资源为MATLAB/Simulink仿真实现基于牛顿第二定律的可简化微分方程模型面向控制工程、机器人学及自动化相关学习者用于掌握系统建模、仿真流程与稳定性控制方法。资源以rar压缩包发布共11个文件包括6张PNG结果/代码截图、3个M脚本、1个SLX模型及1个SLXC缓存压缩包大小仅159KB文件按模型、代码与结果的组合方式排列便于快速定位和对照验证。目前已有3419人学习下载压缩包内仿真覆盖了从动力学方程到Simulink模块的完整过程包含位移变化曲线、角度变化曲线和运算结果可以清晰观察摆杆在稳定与不稳定平衡点附近的动态行为并理解重力和支撑力对系统的作用。同时可通过调整LQR或PID参数分析不同反馈策略对系统稳定性的影响并对照截图理解参数设定与代码实现细节。整体内容轻量紧凑适合课程设计、毕业设计或控制算法入门其研究思路对无人机、机器人平衡等实际应用也具有参考价值是一份可直接运行和二次开发的实用资料。1. 一阶倒立摆的 MatlaB 仿真为什么这个简单模型值得反复折腾一阶倒立摆的 matlab 仿真几乎是每个控制方向学生都绕不过去的东西。它看起来只是一个滑块顶着一根杆子物理结构简单到不能再简单但它天然自带一个正实部极点——也就是说不施加控制摆杆必然倒下。把这个系统在 matlab 里建出来、用状态空间描述、再设计控制器把它“锁”在竖直位置这套流程覆盖了现代控制理论最核心的建模、分析和综合能力很多实验室招新人也会拿它当试金石。这份资源适合两类人一类是做课设或毕设、需要完整跑通建模到仿真全流程的人另一类是刚学完状态空间、想找一个具体对象把 LQR、极点配置这些概念落到代码上的从业者。它能帮你解决的问题很直接动力学方程怎么推、A 和 B 矩阵怎么写、LQR 权重怎么配、Simulink 里怎么搭出来。2. 动力学建模拉格朗日方程与状态空间矩阵的完整推导2.1 为什么用拉格朗日方程而不是牛顿法第一版我做倒立摆建模时习惯性地想用牛顿第二定律硬推。小车一个方程、摆杆一个方程中间还要处理铰链处的约束力写出来一大堆中间变量稍不注意符号就乱。更麻烦的是摆杆在转动同时质心又在平动平动和转动耦合在一起受力分析图画错一个方向后面全废。换成拉格朗日方程之后思路就清晰很多。拉格朗日法的核心是不需要拆解铰链约束力只需要写出系统的总动能 T 和总势能 V然后代入拉格朗日方程。系统的运动方程由能量唯一决定这对于带约束的多刚体系统尤其好用。一阶倒立摆的动能由三部分组成小车的平动动能、摆杆质心的平动动能、摆杆绕质心的转动动能。 势能只需要考虑摆杆质心高度变化。这里有一个常见的细节坑摆杆绕自身质心的转动惯量 J 是单独一项不能和摆杆平动混在一起否则后面矩阵里的耦合项会差一个量级。2.2 从运动方程到四阶状态空间推导与代码验证设小车质量 M摆杆质量 m摆杆质心到铰链距离 l摆杆绕质心转动惯量 J摆角 θ 定义为摆杆与竖直向上方向的夹角小车位移为 x外力为 F。系统的总动能和势能分别为T ½Mẋ² ½m(ẋ² 2lẋθ̇cosθ l²θ̇²) ½Jθ̇²V mgl·cosθ代入拉格朗日方程后得到两个耦合的非线性运动方程(Mm)ẍ mlθ̈cosθ - mlθ̇²sinθ Fmlẍcosθ (Jml²)θ̈ - mgl·sinθ 0这两个方程就是原始模型注意第二个方程的 -mgl·sinθ 项来自势能对 θ 的偏导。接下来做线性化在倒立平衡点 θ0 附近令 sinθ≈θ、cosθ≈1并忽略高阶小量 θ̇²方程组就变成(Mm)ẍ mlθ̈ Fmlẍ (Jml²)θ̈ - mglθ 0从这个线性方程组解出 ẍ 和 θ̈再定义状态向量 x [x, ẋ, θ, θ̇]ᵀ就得到标准的状态空间表达式 ẋ Ax Bu。这里引入中间量 Δ (Mm)J Mml²A 矩阵和 B 矩阵的推导结果如下A [0, 1, 0, 0; 0, 0, -m²gl²/Δ, 0; 0, 0, 0, 1; 0, 0, (Mm)mgl/Δ, 0]B [0; (Jml²)/Δ; 0; ml/Δ]注意 A 矩阵右上角的 2×2 块结构左上角是积分链 [0 1; 0 0]右下角是 [0 1; a₄₃ 0] 的形式其中 a₄₃ (Mm)mgl/Δ 是个正数这决定了系统必然有一个正实部极点。直接在 matlab 里建矩阵并验证开环极点%% 一阶倒立摆动力学参数 M 0.5; % 小车质量, kg m 0.2; % 摆杆质量, kg l 0.3; % 摆杆质心到铰链距离, m J 0.006; % 摆杆绕质心转动惯量, kg*m^2, 按细杆均质近似 J m*l^2/3 g 9.8; % 重力加速度, m/s^2 Delta (M m)*J M*m*l^2; % 推导中出现的公共分母项 A [0 1 0 0; 0 0 -m^2*g*l^2/Delta 0; 0 0 0 1; 0 0 (Mm)*m*g*l/Delta 0]; B [0; (J m*l^2)/Delta; 0; m*l/Delta]; C eye(4); % 四个状态全可观 D zeros(4,1); % 开环极点必然有一个正实部, 这是倒立摆不稳定的根源 openloop_poles eig(A)运行之后 openloop_poles 会显示一对共轭纯虚数零点外加一对实根 ±5.58 左右。正实部那个极点对应的就是摆杆倾倒模态时间常数只有 0.18 秒也就是说不控制的话一个小扰动在 0.2 秒内就会让摆角明显偏离。这个数据很有用后面判断控制器响应速度是否够快就靠它做基准。从线性方程组还可以推出摆角对作用力的传递函数Φ(s)/F(s) (ml/Δ) / (s² - (Mm)mgl/Δ)这个传递函数分母没有 s 的一次项说明系统结构上就是一个“重力负刚度 惯性”的二阶不稳定环节比典型二阶振荡环节的阻尼项缺了关键一项。分析和设计控制器时始终要记得这一条你面对的是一个结构不稳定对象不是调调增益就能稳定下来的。2.3 一组能直接用的物理参数资源里默认用的是下面这组参数尺寸接近实验室常见的直线倒立摆教学平台但不需要匹配特定硬件型号参数符号数值单位小车质量M0.5kg摆杆质量m0.2kg质心到铰链距离l0.3m转动惯量J0.006kg·m²重力加速度g9.8m/s²J 的选取按细杆均质模型 J m·l²/3 0.2×0.09/3 0.006 计算如果后续要用到摆杆是空心管或带配重的情况J 需要按实际质量分布重新算不能照抄。这组参数下 Δ 0.0132闭环反馈设计都在这个基准上做。3. 控制器设计LQR 权重矩阵怎么配才能一次稳柱3.1 为什么这里选 LQR 而不是一味调 PID倒立摆不是不能用 PID工程上也确实有人用。但倒立摆是四阶系统只有一个控制输入 FPID 的三个增益要同时稳住摆角和小车位移两个输出调参过程非常依赖经验和试凑。更本质的问题是 PID 是输出反馈而倒立摆的状态里有四个量都是控制需要的——摆角要拉回零、摆角速度要抑制、小车位移要回归零点、小车速度也要阻尼。PID 只看摆角误差相当于把一个四阶系统的信息压缩成一个标量再用丢掉了太多状态信息。LQR 用的是状态反馈 u -Kx设计过程是一次性把四个状态全部纳入考量通过权重矩阵 Q 来声明“哪个状态重要”再由黎卡提方程算出最优增益。这个过程有明确的数学依据不需要手工试凑而且得到的 K 矩阵天然满足一定的鲁棒性裕度。 对于课设和仿真验证来说LQR 是最稳妥、也最容易讲清楚设计逻辑的方案。3.2 Q、R 矩阵怎么设置四个权重的物理意义与整定顺序Q 矩阵是对角阵 diag(q₁, q₂, q₃, q₄)四个对角元分别对应小车位移、小车速度、摆角、摆角速度的权重。权重越大意味着闭环后该状态被压得越狠。我这里一个实际经验是权重不要按数值大小去“猜”要按物理量纲去理解q₃ 对应的是摆角单位是 rad权重 500 意味着控制器会优先保证摆角回到零——这对倒立摆是生死攸关的而 q₁ 对应小车位移权重 100 表示在摆杆稳定之后小车也要回到初始位置附近。R 是控制输入的权重R 越小控制力越不受约束系统响应越快但执行器的负担也越大。我一般习惯把 R 先固定为 1然后只调 Q。注意控制量 F 的量纲是牛顿如果执行器实际最大推力只有 ±10 NR 就不能设得太小否则算出来的 K 会让控制力频繁撞到饱和限幅。整定顺序上我的习惯是先固定 q₃、q₄ 这一组“摆杆组”权重初值取 q₃500、q₄50再调 q₁、q₂ 这组“小车组”权重让小车位移也能收敛最后回过头看控制量峰值。R 放在最后动只有在控制量长期处于饱和状态时才去调它。资源里的默认权重就是这样一组已经验证过能稳定工作的配置%% LQR 控制器设计 Q diag([100 10 500 50]); % 对应 [x, x_dot, theta, theta_dot] R 1; % 控制力权重 % 可控性校验状态不完全可控时 lqr 会直接报错或给出病态结果 if rank(ctrb(A, B)) 4 error(系统不完全可控, 请检查模型); end K lqr(A, B, Q, R) % 状态反馈增益, 1x4 向量 % 闭环状态矩阵 Acl A - B*K; closedloop_poles eig(Acl) % 验证所有极点都在左半平面这段代码的关键点是 lqr 函数直接返回黎卡提方程的解不需要手写求解过程。但有一个前提条件必须满足系统 (A, B) 可控。rank(ctrb(A, B)) 必须在 lqr 之前校验这是血泪经验——我做第一版的时候直接跳过可控性检查结果某一组参数下 lqr 返回了一个带 NaN 的 KSimulink 里一跑就发散半天没找到原因。K 矩阵的量纲值得留意K 的四个元素分别对应 [N/m, N/(m/s), N/rad, N/(rad/s)]因为 u -Kx 中 x 的物理量纲各不相同。如果你的状态向量把角度用了度而不是弧度那 K 的第三、四列数值会差 57.3 倍这是仿真“玄学”发散最常见的来源之一。3.3 完整闭环仿真代码从 lqr 到 ode45控制器设计完成后我习惯先用 ode45 做一次纯数值闭环验证而不是直接上 Simulink。原因很简单ode45 脚本形式跑起来快改一个参数重跑一次只需几秒钟能快速确认 K 矩阵到底能不能稳住系统。脚本写成一个函数文件加一个运行脚本代码组织如下%% 闭环动力学 - 单独存为 inverted_pendulum_dynamics.m function dxdt inverted_pendulum_dynamics(t, x, A, B, K) u -K * x; % LQR 状态反馈 dxdt (A - B*K) * x; % 闭环线性系统 end%% 主仿真脚本 t_span [0 10]; % 仿真时长 10 秒 x0 [0; 0; 0.1; 0]; % 初始摆角 0.1 rad, 约 5.7 度 [t, x] ode45((t, x) inverted_pendulum_dynamics(t, x, A, B, K), ... t_span, x0); % 控制量还原 u zeros(length(t), 1); for i 1:length(t) u(i) -K * x(i, :); end % 绘制摆角与小车位移 figure(Color, w); subplot(2,1,1); plot(t, x(:,3)*180/pi, LineWidth, 1.5); ylabel(摆角 (deg)); grid on; title(LQR 闭环响应); subplot(2,1,2); plot(t, x(:,1), LineWidth, 1.5); ylabel(小车位移 (m)); xlabel(时间 (s)); grid on;这段代码有几个点需要说明。第一初始摆角 0.1 rad 是有意为之的——它对应大约 5.7 度的小角度既在线性化近似范围内又足够让响应曲线有可读性能看出控制器把摆角拉回零的过程。第二控制量 u 是通过逐行计算还原的没有在 ode45 内部输出原因是ode45 的默认输出只在求解点处有值直接用 K 和状态矩阵乘会比这种循环方式简洁但循环方式更直观方便之后在这个位置加扰动或限幅。第三摆角绘制时乘了 180/pi 转成度方便观察工程意义。跑完之后看曲线摆角应该从 5.7 度快速拉回零超调量大约在 10% 到 15% 之间调节时间在 2 秒以内小车位移会先往反方向运动一段再回零位移峰值大约 0.02 到 0.05 米。如果看到摆角曲线以很高的频率震荡且不衰减优先怀疑 K 矩阵符号反了或者 Q、R 权重差距过大导致闭环极点实部过小。4. Simulink 搭建把状态空间装进方块图看动画输出4.1 两种搭建方式S-Function 还是直接用 State-Space 模块在 Simulink 里搭倒立摆有两种常见做法。第一种是使用 State-Space 模块把第 2 章算好的 A、B、C、D 矩阵直接填进去再用增益模块和反馈环组成闭环。第二种是写一个 S-Function把非线性动力学方程写进函数里。对一阶倒立摆来说如果没有特殊要求我推荐直接用 State-Space 模块理由很实际资源里的模型已经做了线性化LQR 也是基于线性模型设计的用线性模块搭建和理论完全一致调试时信号流向一目了然不会引入额外的数值问题。搭建步骤按顺序执行先从 Simulink Library Browser 拖入 State-Space 模块双击打开后把 A、B、C、D 四个矩阵填进去。这里有一个容易忽略的细节C 矩阵要设成 eye(4)因为状态反馈需要四个状态全部作为输出如果你只关心摆角C 取 [0 0 1 0] 也能跑但 Matrix Gain 模块得到的反馈量就不完整了。D 矩阵直接填 zeros(4,1)。State-Space 模块的每个端口都是向量信号四个状态输出用 Mux 模块合并成一个向量进入 Matrix Gain 模块乘以 K再取负号作为反馈。控制量 u 在反馈环里要经过一个 Saturation 限幅模块限幅值设为 ±10 N模拟实际执行器的推力上限模块关键参数说明State-SpaceA模型A, B模型B, Ceye(4), Dzeros(4,1)四个状态全部输出Initial conditions[0; 0; 0.1; 0]初始摆角 5.7 度Matrix GainK 矩阵 1×4LQR 增益Saturation上限 10, 下限 -10执行器限幅Sum两个输入, 符号 -计算 u -Kx 扰动反馈环最后从 Saturation 输出连回 State-Space 的输入端口就构成了一个标准的单位反馈闭环。这里建议把 State-Space 模块和反馈增益的名称改一下比如 Pendulum_Plant、LQR_Gain这样模型里多个信号时排查连线问题会快很多。4.2 动画与数据导出让结果直观可见Simulink 模型本身就带 Scope 模块但只看波形曲线对倒立摆这种运动对象来说不够直观。我更建议在模型里加一个 matlab Function 模块把每个时刻的 x 和 θ 传进去实时绘制小车和摆杆的运动动画。动画不参与控制计算只做可视化所以它对仿真的稳定性没有任何影响但它能让你一眼看出控制过程是否合理——比如摆杆是单侧收敛还是来回震荡小车是否在大范围漂移。%% 动画绘制函数 - 用于 Simulink Matlab Function 模块 function animate_pendulum(x, theta) % x: 小车位移, theta: 摆角(rad) cart_width 0.15; cart_height 0.08; % 小车矩形位置 x_left x - cart_width/2; x_right x cart_width/2; y_bottom -cart_height/2; y_top cart_height/2; % 摆杆终点坐标 rod_length 0.3; % 可视化摆杆长度 tip_x x rod_length * sin(theta); tip_y rod_length * cos(theta); clf; rectangle(Position, [x_left y_bottom cart_width cart_height], ... FaceColor, [0.5 0.5 0.5]); line([x tip_x], [0 tip_y], LineWidth, 3, Color, k); axis([-1 1 -0.3 0.7]); axis equal; drawnow; end动画函数里 rod_length 取的 0.3 是质心距离但可视化画整根摆杆时用这个长度视觉上不会和数学建模冲突因为动画只反映运动形态不反映质量分布。注意 clf 加 drawnow 的组合每帧都会清除重绘仿真速度为默认的十倍速仿真时画面会流畅很多。如果发现画图本身拖慢了仿真速度可以把动画回调的触发频率降低比如每 5 个仿真步长才刷新一次这个调整不影响结果数据。模型搭建好之后用 Simulink 的 To Workspace 模块把状态信号导出到 matlab 工作区变量名设为 simout。这样后续分析超调量、调节时间、控制量峰值全部可以用脚本处理不用在 Scope 上肉眼读数。导出数据还可以和纯 ode45 脚本的结果做交叉对比——两者的曲线应当完全一致这本身就是验证模型正确性的有效手段。5. 仿真避坑五个最容易翻车的细节5.1 现象模型跑起来直接发散幅值指数增长原因线性化模型符号错误。最常见的是 θ 的定义方向搞反。资源里默认 θ 是摆杆与竖直向上方向的夹角竖直位置 θ0。如果你把 θ 定义成从竖直向下方向起算那么势能项和重力项会出现一个负号A 矩阵右下角的 (Mm)mgl/Δ 会变成负值整个系统极点配置全部错位LQR 按这个模型算出来的 K 符号也反了。解决回到拉格朗日方程检查 V mgl·cosθ 这一项。θ 从竖直向上起算时平衡点处势能取极大值这是倒立摆“倒”的本质。建好 A 矩阵之后先跑一次 eig(A)确认存在正实部极点且数值约为 ±5.58这一步能拦住绝大部分符号错误。5.2 现象LQR 算出 K 后闭环仿真仍然发散原因执行器限幅导致大偏差区间内反馈失效。LQR 是在无约束条件下算的最优增益但 Simulink 模型里加了 ±10 N 的 Saturation 模块后初始摆角较大时控制量瞬间饱和等效增益大幅下降系统进入了线性控制设计时未考虑的工况。解决确认初始摆角在合理范围内。0.1 rad 以下 LQR 可以轻松拉住如果要验证大角度性能需要引入增益调度或改用非线性控制方法单纯调 Q 矩阵解决不了本质问题。仿真时先拿掉 Saturation 模块验证纯线性闭环稳定性再加限幅看控制量是否频繁撞限如果经常饱和适当增大 R 或减小 Q 中的摆角权重。5.3 现象波形高频振荡看起来像数值噪声原因ode45 默认误差容差在小刚度系统上表现良好但 LQR 闭环后系统带宽被推向较高频段默认的缺省容差可能不足以精确跟踪快速变化的状态分量。解决运行 ode45 时显式指定误差容差RelTol 设为 1e-6、AbsTol 设为 1e-8如果波形仍有毛刺就继续收紧一个量级。代价是求解步数显著增加、仿真时间变长但对 10 秒量级的倒立摆仿真来说耗时完全可接受。用 Simulink 仿真时对应位置在 Configuration Parameters 里设置 Solver 的相对误差和绝对误差。5.4 现象小角度下控制正常初始摆角超过 15 度后完全失控原因线性化模型只在 θ0 附近成立。sinθ≈θ 的近似在 15 度时误差已经达到 4%30 度时误差接近 14%摆杆的几何非线性效应已经不能忽略。LQR 基于线性模型设计天然不具备大角度下的稳定性保证。解决把仿真初始摆角限制在 0.2 rad约 11.5 度以内。做课设答辩演示时这个范围足够展示控制器性能。如果任务要求大角度起摆或大范围稳定就要先做摆起控制swing-up再切换 LQR 锁定这是另一个控制问题不在本次资源的线性模型范围内。5.5 现象代码跟资源里一致但结果就是和文档对不上原因单位混用。资源里所有角度运算都用弧度但很多人在设置 Simulink 初始条件或波形观察时用了度。LQR 的 K 矩阵第三、四列依赖角度量纲如果状态向量角度部分是度反馈增益直接就错了。解决全流程统一。初始条件写 0.1 而不是 5.7Scope 显示摆角时用 180/pi 做转换不要在任何中间环节使用角度制。资源里也有一个单位检查函数用来确认 A、B、K 矩阵的数值范围是否合理跑一遍能发现 80% 的量纲错误。6. 从仿真到实物前的最后一道验证鲁棒性测试仿真跑通只是第一步真正决定这套控制器能不能往实物上移植的是鲁棒性。电机输出力矩系数、摆杆长度、质量分布任何一项实测值和标称值有偏差闭环稳定性都可能被破坏。这里给一个实用的参数摄动测试方法把 M、m、l 三个参数各自在 ±20% 范围内随机扰动每组参数重新计算 LQR 增益然后跑闭环仿真统计摆角调节时间的变化范围。%% 参数摄动鲁棒性测试 rng(2024); N_test 200; settle_time zeros(N_test, 1); for i 1:N_test % 参数摄动 M_i M * (1 0.2*randn()); m_i m * (1 0.2*randn()); l_i l * (1 0.2*randn()); J_i m_i * l_i^2 / 3; % 转动惯量随质量长度同步变化 Delta_i (M_i m_i)*J_i M_i*m_i*l_i^2; A_i [0 1 0 0; 0 0 -m_i^2*g*l_i^2/Delta_i 0; 0 0 0 1; 0 0 (M_im_i)*m_i*g*l_i/Delta_i 0]; B_i [0; (J_i m_i*l_i^2)/Delta_i; 0; m_i*l_i/Delta_i]; % 用原始增益做固定反馈测试, 不重新整定 [t_i, x_i] ode45((t, x) inverted_pendulum_dynamics(t, x, A_i, B_i, K), ... [0 10], [0; 0; 0.1; 0]); % 摆角 2% 稳定时间估算 theta_final 0; % 稳定值为 0 idx find(abs(x_i(:,3)) 0.02*0.1 t_i 0.5, 1, first); if isempty(idx) settle_time(i) NaN; else settle_time(i) t_i(idx); end end sum_valid sum(~isnan(settle_time)); fprintf(稳定收敛试验占比: %.1f%%\n, sum_valid / N_test * 100); fprintf(调节时间中位数: %.2f s\n, median(settle_time(~isnan(settle_time))));这个测试的本质是验证固定增益 K 在模型失配下还能不能稳住系统——工程上叫参数鲁棒性。200 组随机试验里如果收敛比例低于 95%说明 Q 权重设计过于激进控制器对模型依赖太高这时候回去调大 R 或者减小小车的权重让控制量温和一些鲁棒性通常立竿见影地变好。我自己的判断标准是调节时间中位数不超过 2.5 秒且 200 组里失效样本不超过 5 个这套增益才值得往实物上移植。做完参数摄动再补一个脉冲扰动测试仿真进行到 5 秒时给小车施加一个 20 N、持续 0.1 秒的瞬时推力观察摆角最大偏移和恢复时间。这个测试模拟的是实际中常见的碰撞或地面不平导致的冲击对判断控制器的抗扰能力很有参考意义。资源里这套方法跑完之后我对每次仿真任务都会强制过一遍三道检查可控性矩阵满秩、所有闭环极点实部为负、参数摄动 200 组收敛率超过 95%。这三条都满足我才放心把结果写进报告或往下继续做实物。这套流程看着繁琐但比起在实物上翻车返工成本几乎可以忽略。希望你也能把这套检查变成自己的习惯——省下的时间都是自己的。希望帮到你。本文还有配套的精品资源点击获取