倒立摆建模与Matlab仿真:从运动方程到LQR控制
2026/9/20 2:00:17 网站建设 项目流程

简介:面向控制理论与机器人技术学习者的一份倒立摆系统建模与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.5kg
m摆杆质量0.2kg
l摆杆质心到转轴的距离0.3m
b小车与导轨的粘性摩擦系数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·θ̈ = F

m·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) = x,x(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 = [(M+m), 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)dt,Q 里每一项对应状态变量的权重,R 对应控制能量的惩罚。权重选得合理,就能在“站得稳”和“动作不猛”之间找到平衡。

4.2 LQR 权重 Q/R 的物理含义与对照表

权重设置是倒立摆仿真里最值得花时间的环节。Q 矩阵往往取对角形式,对角线上的四个元素分别对应位置误差、速度、角度和角速度的惩罚。它们的物理含义如下表所示。

Q 对角线元素对应状态设大后观察到的现象常见初值
Q(1,1)小车位置 x小车快速回中,但控制力变大10~50
Q(2,2)小车速度 ẋ抑制快速移动,减少冲击1~5
Q(3,3)摆杆角度 θ角度偏差被快速纠正100~500
Q(4,4)摆杆角速度 θ̇抑制角速度振荡10~40
R控制力 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,观察摆杆是否会在平衡点附近出现极限环振荡。此时可以检查模型里缺少的库伦摩擦项,或直接加入一个简单的死区补偿函数,输出曲线会立刻贴近真实电机响应。把这三项补完,由模型驱动的控制设计就可以放心交付到下一步工程验证中。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询