最近完整复现了一套水下机器人AUV的路径规划与MPC模型预测控制跟踪工作,这也是我在IEEE顶刊复现系列里做得最扎实的一个。整套流程用Matlab代码实现,覆盖了从AUV运动学建模、障碍环境下的路径搜索,到MPC跟踪控制器设计的所有环节,最后看着仿真曲线跑通、AUV沿着规划路径绕过所有障碍物稳稳到达目标点时,那种“终于把论文变成代码”的踏实感,真不是随便跑个demo能比的。如果你是正在做AUV、无人船、移动机器人路径跟踪相关课题的学生或工程师,或者想系统学习MPC在强非线性、强约束系统中怎么落地,这篇文章值得花半小时看完。我会先把整个复现方案的逻辑拆清楚,再讲建模、路径规划、MPC设计的具体实现,最后列出我调试过程中真正踩过的坑和排查记录。
1. 项目复现的整体思路拆解
1.1 顶刊方法的核心逻辑:为什么是“规划+MPC跟踪”
先讲清楚AUV任务的大闭环:任务下发后,AUV需要在已知或部分已知的障碍环境里找一条从起点到目标点的无碰路径,然后沿着这条路径运动。路径规划是慢决策,跟踪控制是快动作,两者缺一不可。只做规划不做跟踪,开环控制遇到初始误差和海流干扰就会偏航;只做跟踪不做规划,控制器没有参考轨迹可用,等于让一个没有地图的人闭着眼赶路。
顶刊文章通常把“规划+跟踪”放在一起呈现:规划层在几何空间生成路径,控制层在时间轴上完成轨迹跟踪。有的文章还会在规划层考虑能耗、海流影响,在控制层显式处理推进器饱和与障碍约束。复现时不要急着抄代码,先把文章拆成五个模块:路径生成器、轨迹参数化、AUV模型、MPC优化器、仿真环境。每个模块单独验证,再合起来联调,这是所有复现工作里最值得先做的事。
为什么MPC特别适合AUV跟踪?因为AUV是大惯性系统,推进器有输出饱和限制,水下还有流干扰和测量噪声。MPC基于预测模型做有限时域滚动优化,能提前“看到”前方轨迹的弯曲程度,提前减速或转向,这是纯PID做不好的。但MPC必须要有一条参考轨迹作为输入,所以前面必须配路径规划。这个逻辑理顺了,整个复现项目的架构就清晰了。
1.2 复现方案选型:Matlab技术栈的选择依据
这个复现项目用Matlab而不是C++或Python,不是拍脑袋选的,而是对比过的。Matlab在矩阵运算、优化求解、可视化上的优势太明显了,尤其是quadprog这类QP求解器,接口稳定,调试方便,适合学术复现和参数迭代。C++更适合工程部署,Python生态在机器人领域也不错,但做论文复现时,Matlab是最省事的路径。
我见过不少同学一上来就搭Simulink模型,结果被模型配置、采样时间一致性、S-Function调试折腾到放弃。其实顶刊复现阶段用纯脚本足够,Simulink更多是快速控制原型阶段才需要的工具。我最终选择了纯脚本方案,代码可读性好、可移植性强,后面想转成C++也容易。
补充一点环境相关的经验:Matlab建议R2020b以上,因为quadprog的内点法算法和新版优化工具箱对稀疏问题支持更好;如果要用fmincon做非线性MPC,还需要Global Optimization Toolbox。工具版本不一定要最新,但别用太老的版本,否则有些求解器选项会报错。
| 方案 | 学习成本 | 开发效率 | 实时性 | 适用阶段 |
|---|---|---|---|---|
| Matlab纯脚本 | 低 | 高 | 中 | 学术复现、参数调试 |
| Simulink | 中 | 中 | 中 | 快速控制原型 |
| C++/Python | 高 | 中 | 高 | 工程部署、硬件在环 |
2. 核心环节实现:AUV数学模型与路径规划
2.1 AUV运动学与动力学建模
MPC的预测靠模型,路径规划也要用状态空间,所以建模是整个复现的地基。水下机器人真实运动有六个自由度,但路径规划和水平面跟踪最关心的是x、y和偏航角ψ这三个量,深度方向通常由单独的控制器负责。我做的水平面三自由度简化,在顶刊场景里是非常常见的处理方式。
运动学方程很简单,但离散化时容易出错,写在这里供你对比:
% 状态量: x, y, psi % 输入量: u(前进速度), v(横漂速度), r(偏航角速度) x_dot = u * cos(psi) - v * sin(psi); y_dot = u * sin(psi) + v * cos(psi); psi_dot = r;动力学模型用水平面三自由度形式:M·ν̇ + C(ν)·ν + D(ν)·ν = τ + τ_dist。展开之后可以写成如下方程组:
% 简化水平面动力学(示例参数) m11 * udot = m22 * v * r - d11 * u + tau_u + w_u; m22 * vdot = -m11 * u * r - d22 * v + w_v; m33 * rdot = (m11 - m22) * u * v - d33 * r + tau_r + w_r;这里m11、m22、m33是惯性质量项(包含附加质量),d11、d22、d33是阻尼系数,w_u、w_v、w_r是海流等扰动。示例参数我取的是m11=25、m22=30、m33=12、d11=3、d22=5、d33=4,量级上接近中轻型AUV的水动力系数。真实项目里这些参数需要靠CFD计算或实验辨识,仿真复现阶段用合理量级的值就行。
建模里最需要注意的一个工程折中:MPC控制器内部的预测模型不必和仿真被控对象完全一致。顶刊复现中常见做法是控制器用线性化或简化模型,仿真用完整非线性模型,这样能顺带测试控制器的鲁棒性。我把非线性动力学放在auv_dynamics.m里作为被控对象,MPC内部只用线性化状态方程,这个模型失配检验非常有效。
2.2 路径规划算法拆解与路径平滑
路径规划算法我选择RRT作为基础,再针对AUV场景做扩展。为什么不用A?因为A需要栅格化地图,在连续坐标水下环境中会引入栅格尺寸导致的分辨率问题,且AUV在开阔水域移动时,RRT这种基于采样的方法更容易产生平滑路径,也更容易在代价函数里加入能耗、海流等连续代价项。
RRT*的核心流程是采样、扩展、邻域搜索和重连。我给的代码是单次迭代的核心结构:
% RRT* 单次迭代核心流程 q_rand = sample_free_space(map); % 工作空间随机采样 q_near = nearest_node(tree, q_rand); % 搜索最近节点 q_new = steer(q_near, q_rand, step_size); % 向随机点扩展一步 if is_collision_free(map, q_near, q_new) near_nodes = near_neighbors(tree, q_new, r_neighbor); parent = choose_min_cost_parent(near_nodes, q_new); % 选择代价最小的父节点 tree = rewire(tree, near_nodes, q_new, parent); % 重连优化 end实操参数方面,我建议步长取0.5到1.0米,最大迭代次数3000次,邻域半径取步长的1.5倍左右。碰撞检测不要用栅格法,直接把障碍物简化为圆形,判断新路径段到障碍圆心的距离是否小于障碍半径加AUV安全半径:
function free = is_collision_free(p1, p2, obstacles) free = true; for i = 1:length(obstacles) % 障碍物: cx, cy, r d = distance_point_to_segment(obstacles(i).center, p1, p2); if d < obstacles(i).r + 0.3 free = false; return; end end endRRT*直接输出的路径是折线,如果直接把折线丢给MPC当参考轨迹,控制器会在每个折点处产生剧烈控制量,前几条路径全部因为这样跑飞了。解决办法是路径平滑,我用三次样条对路径点做插值,Matlab自带的cscvn可以快速做到。平滑之后还要做时间参数化,把几何路径变成时间轨迹:计算路径累计弧长,按期望速度1.0米每秒赋值时间戳,再在统一时间轴上插值得到x_ref(t)、y_ref(t),航向参考psi_ref(t)由相邻路径点的切线方向计算。这一步是“路径规划”和“轨迹跟踪”之间的桥梁,很多复现代码跑不通,问题都出在这里。
3. MPC模型预测控制设计与实现
3.1 MPC的核心理念与AUV中的适用场景
MPC的原理用一句话就能讲明白:每走一步,看未来几步,算出一个局部最优控制序列,但只执行第一步,然后重新计算。这就像做饭时先想好接下来半小时的流程,每完成一步就重新看一眼时间、重新调整后面安排,而不是一开始就把三小时的计划全部锁死。
具体到AUV上,MPC的控制算法核心是三个东西:预测模型、滚动优化、反馈校正。预测模型用线性化状态方程推算未来状态,滚动优化每个控制周期求解一个带约束的优化问题,反馈校正用当前实际状态覆盖预测初值。这套机制天然适合处理AUV的推进器饱和限制、速度边界、障碍物距离约束,这是PID完全做不到的。
设计代价函数时,我采用典型的形式:位置偏差和航向偏差的加权平方和,加上控制增量惩罚。
J = sum_{i=1}^{Np} (y_ref(k+i) - y_pred(k+i))' * Q * (y_ref(k+i) - y_pred(k+i)) ... + sum_{j=0}^{Nc-1} delta_u(k+j)' * R * delta_u(k+j);这个形式看起来简单,实际落地时有不少门道,后面会详细讲权重怎么设、约束怎么加。
3.2 状态空间预测模型离散化与QP构造
MPC要落地实现,第一步是把连续状态方程离散化。我用的状态量x=[x, y, psi, u, v, r]共6维,控制输入是推进器推力tau_u和转艏力矩tau_r二维。在工作点处做线性化得到连续状态矩阵A_c和B_c,然后按仿真步长h=0.1秒做欧拉离散:
h = 0.1; A = eye(6) + h * A_c; B = h * B_c;预测时域Np=20,控制时域Nc=5,这意味着一次要预测未来2秒的轨迹、优化未来0.5秒的控制序列。这个时间尺度对AUV这种大惯性系统是比较合适的,太短预测不到轨迹弯道,太长求解压力大。
然后构造标准QP问题。决策变量是控制序列U=[u(k), ..., u(k+Nc-1)],通过预测模型可以推出状态序列与决策变量的线性关系,把代价函数写成QP标准型。
function [u_opt, exitflag] = solve_mpc(x0, ref_seq, U_prev) [H, f, Aineq, bineq, lb, ub] = build_qp(x0, ref_seq, U_prev); opts = optimoptions('quadprog', 'Display', 'off', 'Algorithm', 'interior-point-convex'); [U_star, ~, exitflag] = quadprog(H, f, Aineq, bineq, [], [], lb, ub, U_prev, opts); u_opt = U_star(1:2); % 只取第一个控制量执行 end约束部分,我分别加了输入饱和约束和变化率约束:
% 控制量上下限 lb = repmat([-tau_u_max; -tau_r_max], Nc, 1); ub = repmat([tau_u_max; tau_r_max], Nc, 1); % 控制增量约束通过 Aineq 表达 Aineq = build_delta_u_matrix(Nc, 2, delta_max); bineq = zeros(2 * (Nc - 1), 1);关于避障约束,我建议在MPC层暂时用软约束处理,而不是把非线性障碍约束直接塞进QP。位置避障要写成||p - p_obs|| >= d_safe,这个非凸约束会让求解器频繁报不可行。我的做法是把障碍物当作惩罚项加到代价函数里,或者先保证规划层已经绕开了障碍物,MPC层专注于跟踪,工程上更稳妥。
3.3 求解器配置与控制器参数调试经验
参数调试顺序很重要。我建议先固定R矩阵,调Q矩阵;然后固定Q,微调R;最后再动预测时域和控制时域。千万不要一上来同时改好几个参数,否则你根本不知道是哪个改动让曲线变好的。
我最终使用的权重矩阵如下:
Q = diag([10, 10, 5, 1, 1, 0.5]); % 位置误差权重最大 R = diag([0.1, 0.1]); % 控制增量惩罚注意Q矩阵里位置误差是“米”的量级,航向误差是“弧度”的量级,直接放一起会导致弧度误差被淹没。解决办法有两个:要么按物理量级归一化,要么把Q的航向权重加大。我实测下来,Q中x、y取10,psi取5就能有不错的平衡。
预测时域的调试经验是:Np太小会退化成近似静态反馈,Np太大计算时间暴涨。我试过把Np从20调到30,跟踪误差从0.3米降到0.08米,但单步求解时间从12毫秒涨到28毫秒。仿真阶段这个时间还能接受,如果后面要做实时部署,就要在误差和时间之间重新权衡。
求解器配置方面,quadprog建议选interior-point-convex算法,对大矩阵问题更稳。如果遇到Hessian矩阵不正定的报错,常见原因是模型离散化后矩阵条件数过大,解决办法是对状态量做归一化,或者把R矩阵稍微增大一点,给优化问题加点正则项。
4. 整体联调:仿真环境搭建与跑通流程
4.1 仿真环境初始化:障碍物地图、AUV参数与代码结构
联调之前先把仿真环境搭好。我用的是50米乘50米的水域,起点放在(2,2),目标点放在(48,48),在中间随机布置3到5个半径2到4米的圆形障碍物。环境初始化脚本只需要几十行:
% init_map.m map.bounds = [0, 50, 0, 50]; map.obstacles = [25, 25, 3.5; 20, 15, 2.5; 35, 35, 3.0; 15, 40, 2.8]; map.start = [2; 2; 0]; map.goal = [48; 48; 0];代码文件组织上,我按功能拆分成清晰的模块,方便逐模块验证和后续扩展:
main.m:主仿真循环,负责数据流串联init_map.m/init_auv.m:地图和AUV参数初始化plan_path_RRTstar.m:全局路径规划smooth_path.m:路径平滑ref_traj_param.m:轨迹时间参数化,生成参考轨迹mpc_controller.m:MPC控制器主函数build_qp.m:QP问题构造auv_dynamics.m:非线性被控对象模型plot_results.m:结果可视化
这种模块划分的好处是,每个函数都可以单独测试。例如先给mpc_controller.m一个人为构造的正弦参考轨迹,确认跟踪效果没问题,再接上路径规划的实时输出,最后再把障碍物加进去,排查问题的范围能缩小很多。
4.2 联调三步走:先跟踪、后避障、再完整闭环
整个联调我分了三个阶段,每个阶段都有明确的通过标准。
第一阶段是固定轨迹跟踪。随便给一条光滑参考轨迹,比如直线加圆弧的组合,让MPC去跟踪。这个阶段不需要路径规划,只要验证控制器本身没问题。如果这个阶段都在发散,那不是路径规划的问题,是MPC的参数或离散化的问题。
第二阶段是规划器的独立测试。调用plan_path_RRTstar.m生成路径,再用smooth_path.m和ref_traj_param.m转成时间轨迹,但先不要让MPC参与闭环,直接在图上对比规划路径、平滑后的路径和时间轨迹重投影路径是否一致。这一步能快速暴露平滑插值的时间参数化错误。
第三阶段才是完整闭环。主循环里每个控制周期执行:当前状态 → 从参考轨迹切片取未来Np步参考序列 → 求解MPC → 把第一个控制量输入给AUV非线性动力学 → 积分一步得到下一时刻状态。
for k = 1:N_sim x_cur = x_history(:, k); ref_seq = get_ref_slice(ref_traj, t_k, dt * Np); [u_k, flag] = solve_mpc(x_cur, ref_seq, U_pre); [~, x_next] = ode45(@(t, x) auv_dynamics(x, u_k), [t_k, t_k + dt], x_cur); x_history(:, k+1) = x_next(end, :)'; U_pre = [u_k; u_k]; % 为下一次迭代准备初值 end完整闭环跑通后,我还加了一个测试:在仿真中途突然把AUV位置沿y方向推偏0.5米,观察MPC能否把轨迹拉回参考线。实测在2秒左右就恢复了,这个测试对判断控制器鲁棒性很有说服力,录像截图放到汇报材料里也好看。
5. 复现过程中的坑与排查技巧
5.1 高频问题速查表
复现过程中遇到的问题,我整理成了一张速查表,基本覆盖了新手能遇到的大部分情况:
| 问题表现 | 可能原因 | 解决办法 |
|---|---|---|
| 轨迹发散 | 模型离散化错误、状态初值给错、权重过大 | 检查A/B矩阵、初值状态重置、降低Q值 |
| QP求解失败 | Np过大、约束过紧、Hessian不正定 | 降低Np/Nc、放宽控制约束、增大R正则项 |
| 规划路径穿障碍 | 采样步长过大、碰撞检测精度不足、平滑后偏移 | 缩小步长、增加安全半径、平滑后重新碰撞检测 |
| 跟踪误差恒定大 | 参考轨迹未做时间参数化、权重不匹配 | 均匀弧长采样再插值、调整Q矩阵 |
| 单步计算超时 | Np/Nc过大、用了fmincon | 减小Np/Nc、换quadprog或OSQP |
| 每次运行结果不一样 | RRT*随机采样未固定种子 | 设置rng(42)固定随机种子 |
其中“参考轨迹未做时间参数化”是我见过最多的问题。很多代码里直接把规划出的路径点当参考轨迹,但路径点之间的时间间隔是随机的,MPC预测时参考序列对应的时刻是乱的,跟踪误差自然降不下来。正确做法是先把路径按弧长均匀采样,再按期望速度分配时间戳。
5.2 三个影响复现成败的细节教训
第一个教训:MPC里位置状态约束要谨慎。我一开始想在MPC里直接加避障约束,结果求解器频繁报“不可行”。后来把避障从硬约束改成惩罚项,或者干脆交给规划层处理,问题才解决。MPC的硬约束越多,可行域越窄,一旦预测时域内有约束相冲突,整个QP直接无解,比你想象的脆弱得多。
第二个教训:RRT*是一种随机算法,每次跑的结果都不一样。如果没有固定随机种子,你会发现同一个场景每次路径都不一样,复现结果也讲不清楚。设置rng(42)让采样序列可复现,这只是个细节,但直接影响复现工作的可信度。我后来在main.m开头固定种子后,路径规划结果完全一致,后续调参也能严格对比了。
第三个教训:仿真被控对象和MPC内部模型要用不同的函数文件。我一开始图省事,MPC预测和实际仿真共用同一个模型函数,结果控制器“惊讶地”获得了完美的模型先验,跟踪效果看起来好得离谱,但一加扰动就崩。真正合理的做法是仿真用非线性模型,MPC内部用线性化模型,保留模型失配,这样复现出来的结果才真实可信,也更接近论文里面对比实验的设置。
5.3 怎么判断复现是否成功
复现不只是“能跑”,还要给出可信的量化结果。我建议从三个维度评价:稳态精度、避障能力、计算效率。
稳态精度我用的指标是到达目标点后AUV最终位置与目标点的距离,我这次复现结果稳定在0.05米以内。避障能力是随机生成50个不同障碍场景,统计无碰到达目标的成功率,我测下来是48/50,两个失败场景都是因为障碍物把可行通道堵得过于狭窄,这是规划问题的极限场景,不是控制器的锅。计算效率方面,MPC单步求解平均时间在15毫秒以内,低于0.1秒的仿真步长,满足实时性要求。
另一个判断复现是否成功的标准是看原文的实验设置。很多顶刊文章会给出特定的障碍布局、初始条件和参考轨迹形式,复现时不要自己随便改场景来“凑效果”。场景改动会导致整个对比失去意义,这一点务必守住。
6. 从仿真到应用:后续扩展空间
6.1 算法层面的扩展方向
这套复现代码的扩展空间很大。算法层面,我现在用的是线性MPC,可以升级到非线性MPC,直接对非线性动力学做预测,用fmincon或SQP求解,能处理大角度机动和强非线性工况。代价是计算量成倍增加,通常需要配合实时求解器才能落地。
另一个扩展方向是自适应MPC,把模型参数在线辨识嵌入到MPC更新中。AUV的水动力系数在不同航速和深度下会变化,固定参数的MPC到真实场景肯定吃亏。我在仿真里加了一个简单的递推最小二乘估计模块,在线更新阻尼系数,跟踪误差比固定参数时又降了一截。
路径规划侧也能升级。基础RRT只优化路径长度,实际AUV任务还要考虑能耗、海流、安全裕度,可以把多目标代价函数嵌入到RRT的扩展和重连阶段。比如顶刊里常见的做法是给每条路径计算一个“流场中的实际抵达时间”,用时间代价替代距离代价,效果比单纯最短路径好得多。
6.2 工程落地层面的建议
如果要把这套算法搬出Matlab,往实际AUV或无人船硬件平台走,有几个问题要提前规划。
代码移植到C++或Python时,算法结构可以保持不变,但求解器要换成更适合实时计算的方案,比如OSQP或写一个定制的内点法求解器。我在Matlab里用quadprog,移植到C++后如果继续用通用QP库,性能可能不够。一个工程优化思路是:预计算与迭代无关的矩阵(比如Hessian矩阵中与控制量无关的部分),在控制循环里只更新与当前状态、参考轨迹相关的向量。
硬件在环测试时必须把状态估计误差加进来。仿真里所有状态都是直接可得的,但真实水下环境中定位靠惯导、声呐或视觉,位置和航向估计都有延迟和噪声。建议先在仿真环境里给状态加高斯噪声和一拍延迟,测试控制器的容错能力,再上硬件,否则很容易出现“仿真效果好、下水就没法用”的尴尬情况。
最后再多说一句关于“真实感”的问题:水下环境里海流是时变的、非均匀的,MPC的预测模型里如果能把流场估计值作为前馈项补偿进去,效果会好很多。我后面计划在预测模型里增加一个海流速度状态,把洋流当干扰估计器处理,这样MPC的预测轨迹会明显更贴近实际情况。这条路值得做水下方向的朋友尝试。
复现完这套流程之后,我自己最大的感受是:顶刊文章的仿真视图看起来很漂亮,但任何一个参数改动都可能让曲线完全变样。所以拿到代码后,先从最简单场景跑通,再一步步加约束、加扰动,每次只改一个变量。另外,多花时间在模型离散化和参考轨迹时间对齐上,比盲目调Q、R矩阵的性价比高得多。这是我踩过几次坑之后最想提醒你的经验。