☰
一阶倒立摆Simulink建模与PID/LQR控制实战
2026/10/26 22:52:44 网站建设 项目流程

简介:一阶倒立摆系统是力学与控制领域研究动态稳定性的经典模型,在机械工程、机器人学和自动化中具有广泛的应用价值。这份资源面向需要借助MATLAB/Simulink完成建模、仿真与控制设计的工科学生与自学者,内容围绕动力学方程构建、Simulink模型搭建、控制器设计与仿真验证形成完整流程。压缩包大小约14KB,文件总数显示为0,文件类型明细暂未提供,内容集中在建模与仿真分析,适合在已有MATLAB环境的条件下参考学习。目前已有3844人浏览学习,具有不错的资料热度。通过该资源,读者可掌握用符号计算工具箱将一阶非线性微分方程转换为仿真模型的方法,理解PID控制中比例、积分、微分参数对稳定性和响应速度的影响,并借助仿真阶段的可视化输出评估控制效果,进一步了解LQR、滑模控制、自适应控制等高级策略,为机器人平衡控制等项目实践提供从理论到仿真的支持。

1. 一阶倒立摆建模仿真与控制:从动力学方程到 Simulink 里稳住那根杆

一阶倒立摆是控制理论里最典型的“看着简单、调起来要命”的实验对象:一个支点、一根杆,目标是让它保持直立。很多做机器人、自动化方向的人第一次接触 PID 和状态反馈,都是从这个系统入手的。这份资源的核心价值在于把一阶倒立摆的完整链路串了起来——从手推动力学方程,到 Simulink 里搭出仿真模型,再到用 PID 和 LQR 把杆子稳住,全程有可复现的步骤和参数。它适合正在做课程设计、考研复试项目或者刚接触控制仿真的从业者,直接照着搭模型、调参数,而不是对着教材干看公式。

2. 先把数学模型写对:拉格朗日方程与线性化的四个关键点

2.1 为什么必须从拉格朗日方程入手

一阶倒立摆在 Simulink 里搭模型之前,数学表达必须是对的。常见做法是用拉格朗日方程推导,而不是直接在 Simulink 里凭感觉连积分块。原因是倒立摆是一个典型的非保守系统,重力、惯性力、摩擦力耦合在一起,直接用牛顿第二定律拆受力分析容易遗漏力矩项。

一阶倒立摆的理想化模型由一个沿水平方向移动的小车和一根固定在小车上的摆杆组成。设小车质量为 M,摆杆质量为 m,摆杆长度为 l,摆杆与竖直方向的夹角为 θ,小车受到水平外力 F。系统的广义坐标取 x 和 θ,拉格朗日量 L = T - V,展开后得到两个耦合的动力学方程。

常见推导结果是:

  • 小车平动方程:关于 x 的二阶导数与 θ 的二阶导数、θ 的一阶导数平方项耦合
  • 摆杆转动方程:关于 θ 的二阶导数、θ 的一阶导数、x 的二阶导数耦合,含重力项

这组方程里存在 sin 和 cos 项,属于非线性微分方程,直接用于控制设计会很麻烦。工程上的标准处理方式是:在 θ = 0 附近做小角度近似,即 sin θ ≈ θ、cos θ ≈ 1,并把高次项和耦合的非线性项忽略掉,得到一个线性化的状态空间表达式。资源中给出了完整的推导结果和化简过程。

2.2 状态空间表达式与参数选取

线性化处理之后,系统可以写成标准的 x_dot = A x + B u 形式。取状态变量为 x、x_dot、θ、θ_dot,其中 x_dot 是小车速度、θ_dot 是摆杆角速度,控制输入 u 为小车加速度。A 矩阵中包含了小车质量、摆杆质量、摆杆长度、重力加速度等参数。

在 MATLAB 里敲这段代码可以把 A、B 矩阵直接算出来:

% 物理参数 M = 0.5; % 小车质量,单位 kg m = 0.2; % 摆杆质量,单位 kg b = 0.1; % 小车摩擦系数,单位 N/(m/s) l = 0.3; % 摆杆长度,单位 m g = 9.8; % 重力加速度,单位 m/s^2 I = 0.006; % 摆杆转动惯量,单位 kg*m^2 % 线性化后的状态矩阵 % 状态: [x; x_dot; theta; theta_dot] p = I*(M+m) + M*m*l^2; A = [0 1 0 0; 0 -(I+m*l^2)*b/p (m^2*g*l^2)/p 0; 0 0 0 1; 0 -(m*l*b)/p m*g*l*(M+m)/p 0]; B = [0; (I+m*l^2)/p; 0; m*l/p]; C = eye(4); % 观测全部状态 D = zeros(4,1); % 检查开环极点 eig(A)

这段代码的逻辑是先把转动惯量 I 纳入计算,避免把摆杆当成质点导致模型偏简单;p 是分母项,由系统全部质量和几何参数组合而成。代码运行后 eig(A) 会得到一个正实部极点,这正是倒立摆开环不稳定的数学体现——任何微小的角度偏移都会被重力放大。参数里的 b 是摩擦系数,很多资料直接设为零,但实际仿真中保留摩擦会让模型更接近真实装置,调出来的控制器参数也有参考价值。

2.3 模型验证:开环响应必须发散

拿到 A、B 矩阵之后,不要急着搭 Simulink,先用 MATLAB 验证一下开环响应是否符合物理直觉。常见的做法是给一个微小的初始角度,比如 θ = 0.1 rad,然后观察系统状态随时间的变化。

% 开环仿真验证 sys = ss(A, B, C, D); x0 = [0; 0; 0.1; 0]; % 初始角度 0.1 rad t = 0:0.01:2; % 无控制输入 [u, ~] = size(B); u_zero = zeros(length(t), u); [y, t_out, x] = lsim(sys, u_zero, t, x0); figure; plot(t_out, x(:,3)); xlabel('时间 (s)'); ylabel('摆杆角度 (rad)'); title('开环响应:角度随时间发散');

这里用 lsim 做线性系统的仿真,u_zero 是一个与时间等长的零输入序列,代表没有施加任何控制力。运行后角度曲线应当迅速偏离零点并单调发散,如果出现振荡收敛或者曲线贴合零点的情况,要么是 A 矩阵里摩擦项设得过大,要么是状态矩阵拼错了。这一步很重要,它确认了数学模型没有原则性错误,后续 Simulink 仿真才有意义。

3. Simulink 模型搭建:从积分块到子系统封装的完整流程

3.1 Simulink 里复现动力学方程的方式

拿到线性化模型后,Simulink 里的实现有两种主流方式:一种是用积分块直接搭微分方程,另一种是直接用 State-Space 块填 A、B、C、D 矩阵。前者能看清信号流,适合教学和确认物理结构;后者简洁、不容易连错线,适合做控制设计。资源里两种都有,平时调试我一般先用 State-Space 块快速打通闭环,再搭一个积分块版本做交叉验证。

先看 State-Space 块的用法,双击模块填入 A、B、C、D 四个变量名,Simulink 会自动识别工作区里的矩阵。输入端接控制信号 u,输出端接 Demux 后分别引出 x、x_dot、θ、θ_dot 四路信号。这里的 C 矩阵取单位阵,是为了把全部状态都暴露出来,方便后续 PID 或 LQR 使用。

如果选择积分块搭法,结构上分成三个回路:小车速度回路、摆杆角速度回路、小车位移与摆杆角度积分回路。θ_dot 经过积分变成 θ,θ 反馈回去参与计算;x_dot 积分变成 x。这种搭法对初学者更直观,但信号线多,容易把反馈符号搞反。

3.2 参数加载脚本与仿真配置

Simulink 模型一般不写死参数,而是从工作区读取。资源里提供了一个参数加载脚本,打开模型前先运行它:

% 参数初始化脚本 init_par.m clear; clc; % 物理参数 M = 0.5; % 小车质量 kg m = 0.2; % 摆杆质量 kg b = 0.1; % 小车摩擦系数 l = 0.3; % 摆杆长度 m I = 0.006; % 摆杆转动惯量 kg*m^2 g = 9.8; % 重力加速度 m/s^2 % 线性化状态空间矩阵 p = I*(M+m) + M*m*l^2; A = [0 1 0 0; 0 -(I+m*l^2)*b/p (m^2*g*l^2)/p 0; 0 0 0 1; 0 -(m*l*b)/p m*g*l*(M+m)/p 0]; B = [0; (I+m*l^2)/p; 0; m*l/p]; C = eye(4); D = zeros(4,1); % PID 初值 Kp = 80; Ki = 30; Kd = 10; % LQR 权重矩阵 Q = diag([1 1 100 10]); R = 0.1; % 仿真时间 stop_time = 10;

把所有参数集中到这个脚本里,好处是改参数不用去 Simulink 里一个个翻对话框,改完脚本重新运行一下就行。Q 矩阵里的权重分配值得多说一句:第三个元素是 100,对应 θ 状态,代表对摆杆角度的惩罚最重;第四个元素是 10,对应角速度,用来抑制抖动。R 取 0.1 表示控制量成本较低,允许控制器输出较大加速度。这套权重在初调时比较合理,后续按需调整。

3.3 从 Simulink 导出数据做定量分析

仿真跑完之后不能只盯着 Scope 看曲线,要导出数据做定量评估。Simulink 里可以用 To Workspace 块导出时间序列,也可以直接在 MATLAB 中用 sim 命令运行模型:

% 运行模型并导出结果 out = sim('cart_pole_sim.slx', 10); % 提取输出信号 t = out.tout; % 时间序列 theta = out.theta.Data(:); % 摆杆角度 x = out.x.Data(:); % 小车位移 xdot = out.xdot.Data(:); % 小车速度 % 稳态性能指标 theta_settle = theta(abs(theta) < 0.01); % 收敛到 0.01 rad 以内的数据 settle_time = t(length(theta_settle)); % 近似调节时间 fprintf('调节时间: %.2f s\n', settle_time); % 超调量 overshoot = (max(abs(theta)) / 0.1) * 100; fprintf('超调量: %.2f %%\n', overshoot);

这里假设模型设置了三个 To Workspace 块,输出变量名分别为 theta、x、xdot。sim 函数的第二个参数是仿真时长,返回的 out 对象里包含全部导出信号。settle_time 的计算逻辑是找到角度首次持续落在 ±0.01 rad 区间内的时刻,这个方法虽然粗糙但够用于比较不同参数下的控制效果。如果发现 settle_time 一直是 0,多半是角度阈值设得太宽,初始偏移 0.1 rad 还没被控制住就被算进去了。

4. PID 控制与参数整定:先把 Kp、Ki、Kd 调到不晃,再谈优化

4.1 PID 在倒立摆里的控制结构

一阶倒立摆的 PID 控制结构通常是双环的:内环控制摆杆角度,外环控制小车位移。如果只控制角度不管位移,小车会一路漂移直到撞上轨道限位。资源里的单级 PID 结构只反馈 θ 这一个状态,适合验证基本原理;实际仿真中更推荐用角度环加位置环的串级结构。

角度环的输入是 θ_ref - θ,其中 θ_ref 为 0,即竖直向上。由于倒立摆角度在平衡点附近可以近似线性,PID 的 P 项起主导作用:角度偏差大了,输出一个反向加速度把它推回去;D 项提供阻尼,防止来回振荡;I 项用来消除稳态误差,但调大了容易引起低频振荡,要谨慎。串级结构里,外环位置 PID 的输出作为内环角度环的目标值,位置偏了就让摆杆轻微倾斜产生加速度,位置回到中心后角度目标恢复为 0。

4.2 用 MATLAB 脚本完成常规整定

在 Simulink 里手动试参数效率低,调试时通常写一个循环脚本,批量跑不同参数组合。下面给一个简化的参数扫描脚本:

% PID 参数扫描 kp_list = [50 80 120]; ki_list = [10 20 30]; kd_list = [5 8 12]; best_ae = inf; best_param = []; for kp = kp_list for ki = ki_list for kd = kd_list % 写入工作区 assignin('base', 'Kp', kp); assignin('base', 'Ki', ki); assignin('base', 'Kd', kd); % 运行仿真 out = sim('cart_pole_sim_pid.slx', 5); theta = out.theta.Data(:); % 计算绝对误差积分 ae = sum(abs(theta)) * 0.01; % 0.01 是步长 if ae < best_ae best_ae = ae; best_param = [kp ki kd]; end end end end fprintf('最优参数: Kp=%g, Ki=%g, Kd=%g\n', best_param(1), best_param(2), best_param(3));

循环里用 assignin 往基础工作区写入 Kp、Ki、Kd,这样 Simulink 模型里的 PID 块引用工作区同名变量时就会自动更新。ae 指标是角度绝对误差乘以步长再累加,等价于误差曲线面积。这套脚本的逻辑是用粗粒度网格搜索找到好的初值,之后再用人工微调。要注意的是 sim 每次调用都会重新初始化模型,所以状态清零没问题,但如果模型里有随机扰动源,多次仿真之间需要固定随机种子,否则比较参数优劣就不公平了。

4.3 参数整定后的边界观察

参数调好后,要做的是边界测试,否则换个初始条件就翻车。我一般会测试三组场景:小扰动(θ0 = 0.05),中等扰动(θ0 = 0.2),大扰动(θ0 = 0.4)。PID 参数是根据小角度线性化模型设计的,在大角度下非线性效应明显,控制器可能无法将摆杆拉回。理想的 PID 参数应该保证 0.2 rad 内能稳定回收,0.4 rad 允许轻微超调但最终收敛。

边界观察的价值在于暴露模型的线性适用范围。比如某个参数组合在小扰动下表现很好,一加大初始角度就发散,那说明 D 项提供的阻尼不够,或者 P 项饱和了。遇到这种情况优先加 Kd 而不是 Kp,因为加 P 会让系统更“硬”,对大角度有强对抗,但也容易激发高频振荡。

5. 一阶倒立摆的五个高频坑:从参数异常到仿真发散的现象与解决

5.1 结构奇异:Simulink 里报「Simulink cannot solve the algebraic loop」

现象:运行仿真直接报代数环错误,模型卡住跑不动。

原因:把角度信号直接用于计算加速度,而加速度又积分回去产生新的角度,形成瞬时反馈环。Simulink 无法在同一个仿真步内解出满足全部代数关系的解,就报错了。常见触发场景是把 PID 输出直接加在加速度输入端,同时又从同一个角度信号取反馈。

解决:在反馈路径上串一个 Memory 块或者 Unit Delay 块,打破瞬时依赖。Memory 块在离散求解器里相当于一个零阶保持器,输出上一时刻的输入值。加上之后模型变成有记忆的离散系统,代数环消除。

另外,检查求解器类型。如果用的是 ode45 这类可变步长连续求解器,代数环更容易出现;换个固定步长离散求解器,比如 ode4,也能绕开部分代数环问题。

5.2 发散但数值不大:角度曲线整体缓慢漂移

现象:仿真能跑完,但是 θ 曲线不振荡,而是缓慢地往一个方向漂,10 秒后偏了 0.5 rad。

原因:控制器积分项累积误差。倒立摆的稳态误差理论上接近零,但如果 I 项增益过大,积分器会持续累积小的角度偏差,输出一个额外的偏置加速度,把摆杆推离平衡点。这个问题在初始角度不是恰好为 0 时特别明显。

解决:把 Ki 下调 50% 再跑一次,如果漂移速度明显下降,说明确实是积分饱和。也可以给 PID 块的积分项加一个限幅,比如限在 ±5 以内,防止积分输出无限增长。

5.3 输出剧烈抖动但角度不发散

现象:角度曲线看上去在零附近,但是放大后全是高频毛刺,控制器输出的加速度曲线更是锯齿状。

原因:D 项对噪声敏感。仿真里的噪声来源一般是量化误差和求解器容差,D 项直接把高频分量放大,产生了过强的修正作用。另一个可能是固定步长求解器的步长太大,导致角度导数计算不准。

解决:先检查求解器步长,固定步长模式下步长不要超过 0.001,否则动力学方程的分辨率不够。如果步长已经很小还抖动,给 D 项的信号路径上加一个低通滤波器,把高频噪声先滤掉再微分。MATLAB 的 PID 控制块里有滤波器系数 N,默认 N = 100,把它调到 20 左右能明显改善抖动。

5.4 小车位移持续增大但摆杆角度稳定

现象:摆杆角度控制得很好,一直稳定在 ±0.01 rad 以内,但小车位置不断往一个方向跑,仿真到后面小车位置已经超出模型定义域。

原因:控制器只用了角度反馈去控制加速度,完全忽略了小车的位移。系统是一个四阶模型,只控制一个状态必然导致其他状态失控。这个现象在单级 PID 里几乎是必然发生的,因为一个 PID 控制器只能控制一个输出。

解决:要么改用串级 PID——外环位置环的输出作为内环角度环的目标;要么直接上状态反馈控制,用 LQR 一次处理四个状态。推荐后者,LQR 天然考虑所有状态,不会出现这种单状态失控的问题。

5.5 Scope 曲线正常但导出数据为空

现象:Scope 里波形显示正常,点击导出数据后得到的变量是空的或者维度不对。

原因:Scope 块自带的日志功能默认可能关闭,或者 To Workspace 变量的输出格式没有设置成数组。资源里如果用的是 Scope 加信号日志的组合,导出时容易出这种问题。

解决:统一用 To Workspace 块导出,双击设置变量名和输出格式为 Array。如果必须用 Scope 保存,需要在 Simulink 配置里把日志数据保存选项打开,并且勾选单仿真输出。最简单的办法还是在模型里加几个 To Workspace 块,写清楚变量名,后续处理用这些变量而不是 Scope 数据。

6. LQR 实战:用状态反馈一次稳住四个状态,参数矩阵怎么设不翻车

6.1 从 PID 到 LQR 的切换理由

PID 控制一阶倒立摆时,默认只反馈角度,忽略位移,结果就是杆子虽然立住了但小车到处跑。LQR 不一样,它天生就是多变量控制器,一次输入四个状态,不存在顾此失彼的问题。LQR 设计的关键是两个权重矩阵 Q 和 R:Q 里不同对角线元素代表对对应状态惩罚的权重,R 代表控制量的惩罚。权重设得不同,控制行为天差地别。

不夸张地说,Q 和 R 的调法有一定的玄学成分,但基础规律是清楚的:

  • Q(3,3) 对应摆杆角度惩罚,权重越高,杆子越“硬”,对扰动回得快,但控制量也大
  • Q(4,4) 对应角速度惩罚,权重高会压制摆动频率,让系统更平稳
  • Q(1,1) 对应位移惩罚,权重高会让小车更快回中心
  • R 是控制量惩罚,R 越小控制越猛,R 越大输出越柔和

资源里给了一组基础 Q、R 值,直接跑能稳定,但未必适合你的物理参数。调整方法:先固定 R = 0.1,Q 取对角矩阵从 diag([1 1 100 10]) 开始,看角度超调量超过 20% 就加大 Q(3,3),控制量太大就加大 R。

6.2 用 MATLAB 计算 LQR 增益并接入 Simulink

代码如下:

% 计算 LQR 增益矩阵 K Q = diag([1 1 100 10]); R = 0.1; [K, S, e] = lqr(A, B, Q, R); % 打印结果 K e % 仿真闭环系统 sys_cl = ss(A - B*K, B, C, D); % 初始状态: 位移 0.1m, 角度 0.2rad x0 = [0.1; 0; 0.2; 0]; t = 0:0.01:5; [y, t_out, x] = initial(sys_cl, x0, t); figure; subplot(2,1,1); plot(t_out, x(:,1)); hold on; plot(t_out, x(:,3)); xlabel('时间 (s)'); ylabel('状态'); legend('位移 x', '角度 \theta'); subplot(2,1,2); plot(t_out, x(:,2)); hold on; plot(t_out, x(:,4)); xlabel('时间 (s)'); ylabel('速度'); legend('dx/dt', 'd\theta/dt');

lqr 返回的 K 是 1×4 的增益矩阵,A - B*K 是闭环状态矩阵,其极点决定了控制性能。e 是闭环极点的向量,所有元素都应该有负实部,否则系统仍然不稳定。initial 函数模拟的是零输入下从初始状态出发的自由响应,用于检查控制器能否把系统拉回平衡点。运行后观察 x 和 θ 两条曲线,它们应当同步收敛到零,并且位移曲线不应该单调偏离——如果 x 先变大再变小,说明位置权重偏低,需要加大 Q(1,1)。

6.3 从仿真曲线量化 LQR 控制效果

LQR 仿真结束后,除了看曲线形状,直接输出三个核心指标用来对比不同参数下的表现:

% 计算性能指标 theta_all = x(:,3); % 角度超调量,初始角度 0.2 rad theta_peak = max(abs(theta_all)); overshoot_pct = (theta_peak - 0.2) / 0.2 * 100; % 调节时间:角度进入 ±0.01 rad 并保持 idx_settle = find(abs(theta_all) > 0.01, 1, 'last'); settle_time = t_out(idx_settle); % 控制能量(近似):用状态反馈 u = -K*x u_all = -K * x'; u_energy = sum(u_all.^2) * 0.01; % 积分近似 fprintf('角度超调量: %.1f%%\n', overshoot_pct); fprintf('调节时间: %.3f s\n', settle_time); fprintf('控制能量: %.3f\n', u_energy);

这段代码的三个指标恰好对应控制器设计的三个维度:超调量衡量动态品质,调节时间衡量快速性,控制能量衡量代价。实际调参时就是一个三角权衡——想让角度回得快,调节时间短,控制能量必然变大;想省控制能量,超调就大,调节时间就长。LQR 的好处是 Q 和 R 的权重调整对这三种指标的影响方向基本可预期,不会出现 PID 那种调了一个参数另两个一起恶化的局面。

从那以后,我做倒立摆仿真前都会先跑一遍开环响应确认模型没问题,再进控制设计。这个习惯帮我避免了很多次在错误模型上调参数的无效劳动。每次换物理参数,我也会强制自己把所有 Q、R、PID 初值写在脚本顶部统一管理,而不是散落在模型各个模块里。希望帮到你。

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

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

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

立即咨询