倒立摆建模与Matlab仿真：从运动方程到LQR控制

简介面向控制理论与机器人技术学习者的一份倒立摆系统建模与Matlab仿真PDF文档系统讲解如何通过控制小车水平运动使摆杆保持垂直这一经典非线性问题。资料完整呈现从物理模型、数学推导到控制器设计、仿真验证的全过程给定摆杆质量、长度、小车质量与重力加速度等参数基于牛顿第二定律建立运动方程在小角度假设下进行线性化并得到含位移、速度、摆角、角速度的状态空间表达式。此后对系统展开能控性、能观性与稳定性分析结合超调量≤10%、调节时间≤4s的指标要求利用状态反馈u-Kx进行极点配置并给出全维观测器与降维观测器的设计方法以解决部分状态不可直接测量的问题。Matlab仿真用于验证系统在不同初始条件下的动态响应帮助读者理解控制器与观测器增益对稳定性的影响。资源共1个PDF文件大小649KB公式推导和矩阵求解过程完整清晰。已有369人学习适合自动化、机器人等相关专业本科生用于课程设计、毕业设计或工程师自学参考。1. 倒立摆建模与 Matlab 仿真为什么先把模型写对倒立摆系统建模及 Matlab 仿真是控制领域里被反复验证过的一条学习路径。很多人以为这是数学建模竞赛或“自动控制原理”课程里的一个例题真正做过一遍才会意识到倒立摆的失稳特性、非线性耦合和实时控制要求远不是一两个公式能覆盖的。你要解决的第一道坎是把直线一级倒立摆的运动方程从物理坐标系里“抠”出来再到 Matlab/Simulink 里用数值积分跑出符合直观的响应曲线最后才谈得上控制器设计。本文围绕“建模 - 仿真 - 控制”这条主线展开适合做机器人、电机驱动、飞控相关仿真验证的工程师也适合准备数学建模竞赛时需要一个可复现被控对象模型的同学。全文不依赖任何现成工具箱脚本可以逐行照抄。2. 倒立摆系统建模的数学模型从牛顿方程到状态空间2.1 参数定义与坐标系先把 M、m、l、b 摆到桌面上倒立摆系统的建模第一步不是列方程而是明确坐标系。这里采用最常见的小车-摆杆模型小车在水平轨道上运动摆杆通过转轴安装在小车上摆杆初始状态是竖直向上我们把这个位置定义为 θ 0。水平向右为 x 正方向逆时针为 θ 正方向。模型中涉及的关键参数如下表所示。符号物理含义典型仿真取值单位M小车质量0.5kgm摆杆质量0.2kgl摆杆质心到转轴的距离0.3mb小车与导轨的粘性摩擦系数0.1N/(m/s)I摆杆绕质心的转动惯量m·l²/3 或 m·l²kg·m²g重力加速度9.8m/s²F施加在小车上的水平外力一般是控制量 uN提示I 在不同文献里有两种取法。如果摆杆被简化为末端集中质量则 I m·l²如果视为均匀细杆则 I m·l²/3。仿真时两者都能用只是自然频率和控制增益会有差异。转轴处没有额外弹性阻尼摆杆与小车之间的铰接视为理想。这样做的目的是先建立一个可控的基准模型后面如果做实物再把摩擦力矩和间隙造成的误差当作扰动加回到模型里。2.2 非线性运动方程小车与摆杆的耦合关系从牛顿第二定律出发可以对小车和摆杆分别做受力分析但更省力的方式是使用拉格朗日方程。拉格朗日函数 L T - V其中 T 是系统总动能V 是总势能。小车只有水平位移速度为 摆杆质心的速度和角度有关写出T 1/2·M·ẋ² 1/2·m·(ẋ² 2·l·θ̇·ẋ·cosθ l²·θ̇²)V m·g·l·cosθ对广义坐标 x 和 θ 分别代入欧拉-拉格朗日方程后得到两个耦合的二阶微分方程(M m)·ẍ b·ẋ m·l·θ̈·cosθ - m·l·θ̇²·sinθ F(I m·l²)·θ̈ m·l·cosθ·ẍ - m·g·l·sinθ 0注意第二个式子里-m·g·l·sinθ 这一项决定了系统在 θ 0 附近的失稳行为。当摆杆稍微向右偏转sinθ 为正角度加速度会把摆杆继续往同一方向推。这个“正反馈”是倒立摆系统建模后第一个要和实物经验对齐的地方不是所有仿真发散都是算法问题模型本身就包含不稳定极点。2.3 在平衡点附近线性化化简为 4 维状态空间模型非线性模型能用于时域仿真但做控制器设计和稳定性分析时必须线性化。取工作点 θ 0、ẋ 0、θ̇ 0并假设小角度范围 sinθ ≈ θcosθ ≈ 1。代入非线性方程后得到(M m)·ẍ b·ẋ m·l·θ̈ Fm·l·ẍ (I m·l²)·θ̈ - m·g·l·θ 0把两个方程联立解出 ẍ 和 θ̈写成状态空间形式。定义状态向量 x_state [x, ẋ, θ, θ̇]ᵀ控制输入为 F则线性定常系统为ẋ_state A·x_state B·u其中系数矩阵按参数表数值代入后是数值矩阵。书写代码时不必手工展开太长的符号公式直接调用 solve 或写矩阵运算即可。3. 用 Matlab 脚本实现倒立摆仿真ode45 跑通开环模型3.1 将二阶方程写成一阶状态方程拆出可执行代码Matlab 里做倒立摆仿真最常见的做法是用 ode45 求解一阶 ODE 组。我们需要把上面的两个二阶方程降阶。设 x(1) xx(2) ẋx(3) θx(4) θ̇则 ẋ(1) x(2)ẋ(3) x(4)剩下 ẋ(2) 和 ẋ(4) 通过联立方程求出。function dx pendCart(t, x, M, m, l, b, g, isLinear) % 状态定义: x(1)小车位置, x(2)小车速度 % x(3)摆杆角度, x(4)摆杆角速度 dx zeros(4, 1); u 0; % 开环仿真, 控制力为零 if isLinear % 线性化模型: A*x B*u I m * l^2; denom (M m) * (I m*l^2) - (m*l)^2; dx(1) x(2); dx(2) ((I m*l^2) * (u - b*x(2)) - m^2*g*l^2*x(3)) / denom; dx(3) x(4); dx(4) ((M m) * m*g*l*x(3) - m*l*(u - b*x(2))) / denom; else % 非线性模型: 求解 2x2 线性方程组 I m * l^2; A [(Mm), m*l*cos(x(3)); m*l*cos(x(3)), I m*l^2]; rhs [u - b*x(2) m*l*x(4)^2*sin(x(3)) m*g*l*sin(x(3))]; acc A \ rhs; dx(1) x(2); dx(2) acc(1); dx(3) x(4); dx(4) acc(2); end end这段代码的核心是质量矩阵 A。非线性情况下的质量矩阵里含有 cosθ所以不再是常数矩阵每次积分步进都需要重新计算。rhs 第一行里多出的 m·l·θ̇²·sinθ 是摆杆摆动带来的向心加速度项在小角度假设下通常被忽略但若要仿真大角度摆动必须保留。线性化与非线性两种模型都写在同一函数里用最后布尔参数 isLinear 切换方便后续对比。运行仿真时使用以下脚本M 0.5; m 0.2; l 0.3; b 0.1; g 9.8; x0 [0; 0; 0.2; 0]; % 初始摆角 0.2 rad tspan [0 5]; opts odeset(RelTol, 1e-6, AbsTol, 1e-8); [t_nonlin, y_nonlin] ode45((t, x) pendCart(t, x, M, m, l, b, g, false), tspan, x0, opts); [t_lin, y_lin] ode45((t, x) pendCart(t, x, M, m, l, b, g, true), tspan, x0, opts); figure plot(t_lin, y_lin(:, 3), b-, LineWidth, 1.5) hold on plot(t_nonlin, y_nonlin(:, 3), r--, LineWidth, 1.5) grid on xlabel(时间/s), ylabel(摆角 θ/rad) legend(线性模型, 非线性模型)ode45 是变步长四阶五阶 Runge-Kutta 算法的实现对于连续非刚性问题足够。RelTol 设置到 1e-6 可以保证稳态轨迹不因积分误差产生虚假振荡代价是计算时间增加仿真 5 秒的物理时间在普通笔记本上通常只有几十毫秒。3.2 开环仿真结果怎么看初始误差后模型会不会发散运行上面的脚本把初始角度设到 0.2 rad大约是 11 度。观察第三列输出也就是摆杆角度的轨迹曲线。开环情况下线性模型会在前 0.5 秒内迅速使角度增大到超出小角度范围随后数值结果趋于无穷非线性模型则在角度增大到一定程度后由于 sinθ 不再近似等于 θ曲线可能会出现回落或周期性摆动。两种模型出现差异不是 bug而是线性化边界本身的体现。时间范围线性模型行为非线性模型行为0~0.3s角度指数增长角度指数增长趋势相近0.3~1s角度持续发散角度增速放缓或翻转1s 以后超出数值范围可能进入摆动或震荡这里有个概念容易混淆仿真发散和系统不稳定不是一回事。开环倒立摆在物理上确实不稳定所以角度增大是模型正确的表现。而数值发散则表现为看似无界但毫无物理规律的跳跃值比如角度超过 10⁶ rad这是积分步长过大或模型切换错误导致的。3.3 把同一模型搬进 Simulink避免代数环的 3 个操作除了脚本倒立摆仿真在 Simulink 里也很常见。基本思路是用两个积分器分别构造小车和摆杆的角度路径再用 Sum 和 Gain 模块搭建耦合方程。搭建时不建议直接使用 Transfer Fcn 模块因为它只支持线性模型无法体现非线性。推荐按以下方式操作。第一使用积分器模块起于初始值两个积分器串联给角度、角速度另两个积分器串联给小车位置、速度。第二把 cosθ 和 sinθ 计算部分拆成独立模块用 Trigonometric Function 引用角度信号再把它们送进矩阵乘法和 Gain 组。第三求解 ẍ 和 θ̈ 时不要人工展开公式而是用 2×2 Mass Matrix 封装积矩阵求逆部分减少符号混乱。仿真配置里需要注意求解器的选择。如果模型是纯连续非刚性的用默认的 ode45 即可如果加入了 PWM 信号、量化区间或者离散 PID 控制器就应该切换到离散求解器并设置固定步长否则可能出现“仿真结果时好时坏”的现象。4. 倒立摆仿真的控制闭环LQR 状态反馈的参数设计4.1 为什么要优先用 LQR而不是单纯调 PID对倒立摆进行控制最直观的想法是调 PID。实际调过的工程师大概率会碰见一组令人头疼的现象比例增益给的太小摆杆刚爬起来就开始倒给的太大会引发剧烈抖动电机电流一路顶到饱和。PID 的问题在于它是单输入控制当摆杆角度、小车位置和速度同时参与约束时缺乏直接处理多输出耦合的手段。LQR 状态反馈控制器天然适合这个场景。设计目标是找到反馈增益矩阵 K使控制量 u -K·x_state 在最小化二次性能指标时同时稳定摆角和小车位置。性能指标是∫(xᵀ·Q·x uᵀ·R·u)dtQ 里每一项对应状态变量的权重R 对应控制能量的惩罚。权重选得合理就能在“站得稳”和“动作不猛”之间找到平衡。4.2 LQR 权重 Q/R 的物理含义与对照表权重设置是倒立摆仿真里最值得花时间的环节。Q 矩阵往往取对角形式对角线上的四个元素分别对应位置误差、速度、角度和角速度的惩罚。它们的物理含义如下表所示。Q 对角线元素对应状态设大后观察到的现象常见初值Q(1,1)小车位置 x小车快速回中但控制力变大10~50Q(2,2)小车速度 ẋ抑制快速移动减少冲击1~5Q(3,3)摆杆角度 θ角度偏差被快速纠正100~500Q(4,4)摆杆角速度 θ̇抑制角速度振荡10~40R控制力 u设大会让控制变软响应变慢0.1~1从控制理论角度讲Q(3,3) 通常要比 Q(1,1) 大一个数量级因为摆杆角度是决定系统稳定性的主导状态。角度权重不足时控制器会为了把小车位置拉回原点而大尺度摆动静止误差仿真曲线呈现明显拖尾。4.3 闭环仿真代码从 K 矩阵到抗扰动测试在 Matlab 中设计 LQR 并完成闭环仿真的代码如下。注意权重矩阵的输入顺序要和状态向量定义严格一致。% 线性化状态空间矩阵 M 0.5; m 0.2; l 0.3; b 0.1; g 9.8; I m * l^2; Delta (M m) * (I m*l^2) - (m*l)^2; A [0 1 0 0; 0 -(I m*l^2)*b/Delta -m^2*g*l^2/Delta 0; 0 0 0 1; 0 m*l*b/Delta m*g*l*(M m)/Delta 0]; B [0; (I m*l^2)/Delta; 0; -m*l/Delta]; C eye(4); D zeros(4, 1); % LQR 权重 Q diag([10, 1, 100, 10]); R 0.1; K lqr(A, B, Q, R); % 闭环仿真, 使用非线性模型, 控制强度限制在 /-15N sat_force 15; u_fun (x) max(-sat_force, min(sat_force, -K * x)); [t, y] ode45((t, x) pendCart_closed(t, x, M, m, l, b, g, u_fun), ... [0 5], [0; 0; 0.1; 0], odeset(RelTol, 1e-6));pusher 自定义函数 pendCart_closed 的主体和开环版一致只是把 rhs 里的 u 替换为 u_fun。代码里 LQR 求解得到 K 后没有直接把 -K·x 施加到系统里而是经过一个饱和限幅。这样做更有工程意义真实电机输出是有上限的仿真里如果不设置饱和权重很大时控制器会指挥电机输出几十牛顿的力看起来曲线很完美实际却不可能实现。提示饱和限幅值设置得太小LQR 在初始误差较大时可能无法恢复直立。可以先跑一个不加限幅的仿真观察 -K·x 的最大瞬时值再取它的 80% 作为限幅值。5. Matlab 仿真的最后一公里用离散控制量替换连续量前的验证LQR 仿真跑通以后不要急着把它当做最终模型。从连续时间仿真转向半实物或嵌入式实现还有三个步骤可以在 Matlab 里直接完成。针对倒立摆系统建模及 Matlab 仿真这条线我认为这三步能省去不少实物调试时间。第一用阶跃响应和初值扰动交叉验证。给位置目标一个 0.1 m 的阶跃观察摆角是否在过渡过程中越过 ±0.3 rad。如果没有越限说明 LQR 增益处于合适范围。这个数据还可以折算成系统带宽为后续数字控制器采样频率的选择提供依据。第二对执行器延迟和噪声做蒙特卡洛仿真。倒立摆仿真发散再收敛的情况经常能在加入 5 ms 执行器延迟后出现。把延迟简化成一个 memory 模块或 state 变量在 ode45 里用 dde23 这类延迟微分方程求解器也可以处理。第三量化仿真结果与实物差异来自哪一层。仿真阶段可以把摩擦系数 b 从 0.1 调到 0.01观察摆杆是否会在平衡点附近出现极限环振荡。此时可以检查模型里缺少的库伦摩擦项或直接加入一个简单的死区补偿函数输出曲线会立刻贴近真实电机响应。把这三项补完由模型驱动的控制设计就可以放心交付到下一步工程验证中。本文还有配套的精品资源点击获取