最近在复现一个很有意思的控制方向项目:把神经网络和模型预测控制揉在一起,分别在四旋翼无人机和非线性机器人汽车两套系统上做仿真验证,全程用Matlab实现。这个方向现在很热,因为传统MPC太依赖精确的机理模型,而神经网络刚好擅长从数据里学动态特性,两者一结合,既能保留MPC处理约束和滚动优化的能力,又能把模型建模的负担甩给数据。这篇文章我把自己复现过程中踩过的坑、核心原理的消化理解、以及代码组织的思路全部梳理出来,给正在做同类项目的同学提供一份可以直接上手的参考。
这个复现项目适合几类人看:一是搞控制理论、想验证NN-MPC联合算法效果的研究者;二是做无人机或移动机器人轨迹跟踪、但苦于动力学模型不精确的工程开发者;三是刚入门MPC、想找一个带完整代码和系统模型的练手案例的学生。内容会覆盖系统建模、神经网络训练、MPC求解器配置、以及联调过程中最常见的发散问题和参数调试技巧,你可以按顺序看,也可以直接跳到对应的坑。
1. 复现前先想清楚:这个项目到底在解决什么问题
1.1 为什么要把神经网络塞进MPC
模型预测控制的核心逻辑大家都知道:用被控对象的模型,往前预测未来一段时间的行为,然后求解一个带约束的最优控制序列,只把第一步下发出去,下一时刻再滚动重算。这个框架本身很优雅,但它最致命的短板就是预测模型要足够准。四旋翼飞行器有强耦合、强非线性的气动特性,参数还容易受负载、风速影响;轮式机器人汽车更不用说,轮胎摩擦、打滑这些因素很难用公式精确描述。机理模型一旦失配,MPC的预测就失真,控制品质直线下降。
神经网络在这时候就能顶上。它不需要你手动推导每一项物理关系,只需要你给足输入输出数据,它就能把系统的那套复杂映射关系学出来。所以NN-MPC的典型搭法就是把神经网络作为MPC的预测模型,或者作为机理模型的误差补偿项。我这次复现选用的是前者——完整用神经网络做一步预测器,在MPC的滚动优化里反复调用这个预测器来生成未来状态轨迹。
1.2 两个被控对象到底差在哪
这个项目同时包含四旋翼和机器人汽车,刚开始我也有点疑惑,后来想明白了:这两类系统代表了两种非常典型的控制难点,放在一起横向对比特别有说服力。
四旋翼属于典型的欠驱动系统——它只有四个输入(四个旋翼的拉力),但需要控制六个自由度的位置和姿态。姿态动力学和位置动力学之间又存在强耦合,绕任何一个轴的转动都会带动平移方向的加速度变化。这种系统对MPC的预测模型要求极高,因为模型里一点小的偏差,经由姿态环放大后,位置轨迹可能整个跑飞。
机器人汽车则完全不同。它是一门典型的非完整约束系统,也就是说,它不能像全向平台那样横着走,只能沿着车身方向前进、后退和转弯。这种系统的运动学方程是简单的[dx/dt; dy/dt; dtheta/dt] = [v cosθ; v sinθ; ω],非线性项主要来自三角函数耦合。它看起来比四旋翼简单,但实际做轨迹跟踪时,初始航向角、转弯半径、参考轨迹曲率都会影响到MPC的收敛性。两个系统放在一起复现,正好覆盖了欠驱动强耦合和非完整约束这两大类难啃的骨头。
1.3 复现路线图与文件组织
我在动手写代码之前先把复现路线拆成了四个阶段,每个阶段都有一个明确的检查点,不然一上来就混在一起很容易卡住不知道查谁。
| 阶段 | 核心任务 | 检查点 |
|---|---|---|
| 阶段一 | 建立两个被控对象的数学模型并验证开环行为 | 给定输入序列,状态不发散,符合物理直觉 |
| 阶段二 | 先用机理模型搭一套线性MPC或NMPC做基线 | 能在简单参考轨迹上跟踪成功,建立对照组 |
| 阶段三 | 采集数据训练神经网络,作为预测模型嵌入MPC | 纯NN预测的单步误差控制在1e-3量级 |
| 阶段四 | NN-MPC联调,对比两组系统下的性能差异 | 跟踪误差收敛、MPC求解时间无明显恶化 |
文件组织我用的是两个独立目录,一个放四旋翼相关的quad目录,一个放机器人汽车相关的vehicle目录,公共的神经网络训练、CasADi接口代码放在shared目录。这样两个系统互不干扰,跑仿真的时候也很清晰。接下来我按照这个路线,把每个阶段的技术细节展开讲。
2. NN-MPC核心原理:不把原理啃透,后面全是坑
2.1 MPC滚动优化的本质和公式记忆法
MPC的原理其实用一句话就能说透:每走一步,就重新规划一次剩下的路。就像你开车导航去一个陌生的地方,导航会给你一条规划路线,但你真开出去一段路之后,导航会根据当前实际位置重新计算新路线。这个"每隔几秒重新规划一次"的动作,就是滚动时域优化。
数学上,MPC每一时刻要解一个如下形式的最优化问题:
min J = Σ_{k=0}^{Np-1} (x(k) - x_ref(k))^T Q (x(k) - x_ref(k)) + u(k)^T R u(k) + Δu(k)^T S Δu(k)
约束:x(k+1) = f(x(k), u(k)),x(0) = x_current,u_min ≤ u(k) ≤ u_max
这里的Np就是预测时域,决定了你的控制器眼睛往未来看多远。Q和R是两个权重矩阵,分别决定状态误差和输入能量在代价函数里的地位。Q取得大,控制器就会更激进地修正状态偏差;R取得大,控制器就会更心疼输入能量,动作变得温柔。这个参数的性格调性,跟你实际调PID的P、I参数是一个道理,只不过MPC的做法更全局化——它会综合未来Np步的误差统一决定当前这一步怎么走。
我复现的时候用一个特别形象的类比帮助理解:MPC本质上是一个"每步都戴着镣铐的贪心算法",但又不完全贪心,因为它的成本函数看到了未来Np步,不算太短视。
2.2 神经网络在项目里扮演什么角色
在这个项目中,神经网络扮演的是"系统动态模拟器"的角色。MPC需要知道"当前状态x(k)、当前输入u(k)作用下,下一步状态x(k+1)是多少"。传统做法是用四旋翼的牛顿-欧拉方程去算,用机器人汽车的运动学方程去算。我这次复现的是用神经网络直接拟合这个映射关系。
具体来说,我训练的是一个多输入多输出的前馈网络。输入特征是当前状态加上当前控制输入,输出是下一步的状态增量或者直接是下一步状态:
x(k+1) = NN(x(k), u(k))
网络结构用的是单隐层+双隐层都试过,最后选择了双隐层,每层12个神经元,激活函数用tanh。之所以不用太深的网络,是因为后面要嵌入MPC优化循环里反复迭代,每多一层非线性计算,CasADi符号求解的负担就重一分。数据采集上,我在机理模型的仿真环境里给控制输入叠加了白噪声激励和扫频信号,确保机器人车和四旋翼的状态轨迹能够覆盖足够宽的工作范围。这里有个非常关键的经验:训练数据一定要覆盖控制器实际工作时可能到达的区域。如果训练数据只在小偏差范围内采集,MPC在误差稍大的情况下查询神经网络,就会得到外推结果,基本等于瞎猜。
2.3 NN与MPC的接口:梯度这关必须过
NN-MPC和传统NMPC最大的区别就在这里。传统NMPC里面,模型是明确的解析表达式,优化器可以轻松算出目标函数对控制输入的梯度;神经网络作为预测模型时,它是一个多层非线性复合函数,想要获得梯度,一个办法是手动去反传推导,另一个更省事的办法是借助符号计算工具。
我在Matlab里用的方案是CasADi工具箱。具体操作分三步走:
- 用Matlab的Deep Learning Toolbox训练神经网络,训练完把网络权重和偏置抽取出来,保存成结构体。
- 在CasADi符号环境里,用这些权重手工重建神经网络的正向传播过程——相当于把每一层的线性变换和tanh激活函数用符号变量重写一遍。
- 把这个符号化的神经网络封装成一个函数,用于MPC成本函数和约束表达式中的状态更新。
这样做最大的好处是,CasADi自动帮你求导。你不需要手算任何一条雅可比矩阵的链式法则,IPOPT优化器会自动通过CasADi的符号微分拿到梯度信息。我第一次踩坑就是天真地想把训练好的net对象直接塞进MPC优化里反复调用,结果Matlab报错或者慢到无法容忍。后来改成提取权重在CasADi里重建,求解速度和梯度精度都正常了。
3. Matlab代码实现:从动力学建模到仿真闭环
3.1 四旋翼模型怎么代码化
四旋翼的完整动力学是13维状态,但复现项目一般不需要全部拉满,可以用简化但保持核心耦合特性的模型。我的状态定义是12维:
x = [x, y, z, vx, vy, vz, φ, θ, ψ, p, q, r]
输入u定义成四个旋翼的拉力T和一个三轴力矩τ = [τφ, τθ, τψ]。在Matlab代码里,状态导数的计算是初始化仿真循环里最核心的部分:
function dz = quadrotor_ode(z, u, params) % 状态解包 pos = z(1:3); vel = z(4:6); euler = z(7:9); omega = z(10:12); phi = euler(1); theta = euler(2); psi = euler(3); % 旋转矩阵,用于将机体坐标系下的力转换到世界坐标系 R = rotation_matrix(phi, theta, psi); % 总拉力沿机体Z轴方向 thrust_inertial = R * [0; 0; u(1)]; % 位置加速度 = 推力加速度 + 重力 accel = thrust_inertial / params.m - [0; 0; params.g]; % 姿态角变化率,这里用简化的运动学关系,避免奇异点 T = [1, sin(phi)*tan(theta), cos(phi)*tan(theta); 0, cos(phi), -sin(phi); 0, sin(phi)/cos(theta), cos(phi)/cos(theta)]; euler_dot = T * omega; % 角加速度,简化忽略气动阻尼,保留惯性耦合项 omega_dot = params.I \ (u(2:4) - cross(omega, params.I * omega)); dz = [vel; accel; euler_dot; omega_dot]; end注意几个细节:姿态角用的欧拉角表示,在θ接近90度的时候会出现奇异性,仿真轨迹里尽量不要设计穿越这个姿态的动作;旋转矩阵转换推力方向这个环节最容易写错,我一开始坐标变换的符号反了,导致四旋翼一推油门就往地面砸,这个坑后面我会再详细说。
四旋翼的MPC里如果直接用全量12维状态做预测,计算负担会很大。我在复现中采用的是分层架构:位置环用NN-MPC计算期望的拉力向量,再反解出期望横滚角和俯仰角,姿态环用经典PID或者简单的几何控制器来跟踪。这样NN-MPC只需要预测三个位置和三个速度量,维度降下来之后,CasADi求解效率提升非常明显。
3.2 机器人汽车(非完整约束)模型怎么代码化
机器人汽车的模型要简单得多,它的核心是非完整约束:
function dz = unicycle_ode(z, u) % z = [x; y; theta] % u = [v; omega] theta = z(3); dx = u(1) * cos(theta); dy = u(1) * sin(theta); dtheta = u(2); dz = [dx; dy; dtheta]; end这看起来只是三个微分方程,但实际做MPC轨迹跟踪时,最大的坑在于:参考轨迹的航向角theta_ref和位置x_ref、y_ref之间必须保持一致性。很多初学MPC的同学参考轨迹只给x_ref和y_ref,不给theta_ref,然后优化器就不知道该把车头朝向哪边,结果就是车一边跟踪一边原地打转。解决办法是计算出参考轨迹的切线角度:theta_ref = atan2(dy_ref, dx_ref),并在MPC代价函数中把航向误差也罚进去。
机器人的MPC输入还有一层物理约束要处理:线速度v和角速度ω都有上下限,复现时我设定的范围是v ∈ [-0.5, 1.0] m/s,ω ∈ [-1.5, 1.5] rad/s。如果把这些约束直接在MPC里写成硬约束,碰到参考轨迹太激进的工况会直接导致优化问题无解,所以我在代价函数里对速度约束用了一种软约束处理方式,让约束越界时付出惩罚代价而不是直接宣布无解,这样控制器在现场总能给出一个可行的折中方案,后面我会展开讲。
3.3 神经网络训练与嵌入MPC的具体操作
这一步是整个复现过程最值得记录的一段,我把代码路径完整走一遍。
先说一下数据采集。我在Matlab里用开环扫频信号驱动两个机理模型,四旋翼的控制输入里混入了幅值很小的白噪声,机器人车的速度输入里做了变速正弦扫描。采集到的输入输出数据存成时序样本对:input = [x(k); u(k)],output = [x(k+1)]。数据量每套系统采了2万多样本,做了归一化和随机打乱划分,训练集和验证集大概6:4。
训练用的是一个双隐层前馈网络:
% 输入维度:状态维度 + 控制输入维度 input_dim = size(Xtrain, 1); output_dim = size(Ytrain, 1); net = feedforwardnet([12, 12]); net.layers{1}.transferFcn = 'tansig'; net.layers{2}.transferFcn = 'tansig'; net.layers{3}.transferFcn = 'purelin'; net.trainParam.epochs = 800; net.trainParam.goal = 1e-6; [net, tr] = train(net, Xtrain, Ytrain);训练完成后,提取权重结构体:
weight_struct = struct(); weight_struct.W1 = net.IW{1}; weight_struct.b1 = net.b{1}; weight_struct.W2 = net.LW{2,1}; weight_struct.b2 = net.b{2}; weight_struct.W3 = net.LW{3,2}; weight_struct.b3 = net.b{3}; save('nn_model_vehicle.mat', 'weight_struct');然后在CasADi里重建这个网络:
import casadi.* function pred_x_next = nn_predict_casadi(x_sym, u_sym, ws) x_u = [x_sym; u_sym]; % 归一化参数必须在训练时保存下来,这里同步使用 x_u_norm = (x_u - ws.x_mean) ./ ws.x_std; z1 = tanh(ws.W1 * x_u_norm + ws.b1); z2 = tanh(ws.W2 * z1 + ws.b2); y_norm = ws.W3 * z2 + ws.b3; % 反归一化 pred_x_next = y_norm .* ws.y_std + ws.y_mean; end这段代码看起来很平淡,但我实际调试时遇到过几个问题,值得单独说:
第一,归一化参数必须在训练之后立刻保存下来,并且嵌入MPC时要用完全相同的均值和标准差。第二次复现这个项目的时候我把归一化参数写成全局变量去手动复制,结果一个符号变量算出来的结果整体漂移了几个量级,位置轨迹直接上天。
第二,激活函数必须在CasADi里使用对应的符号函数。tanh在Matlab里直接写没问题,但如果你训练时用了logsig,这里也要换成对应的表达式,不能想当然。
第三,网络预测的是状态增量还是全量状态,这个选择影响很大。我的做法是预测状态增量Δx = x(k+1) - x(k),然后下一步状态用 x(k) + Δx 来算。这样做的好处是,神经网络只需要学习系统的变化趋势,而不是去硬拟合状态的绝对数值范围,训练误差和预测稳定性都会明显改善。这个细节我强烈建议你复现的时候保留。
3.4 MPC求解器配置与参数调试
NN-MPC的优化求解我用的工具链是CasADi + IPOPT,在Matlab环境下通过casadi脚本接口调用,不依赖Simulink,这样代码更轻量、调试也更直接。
MPC参数我按以下几组来配置,你先记这个基线值,后面再按系统特性微调:
| 参数 | 四旋翼位置环 | 机器人汽车 |
|---|---|---|
| 采样时间 dt | 0.05 s | 0.1 s |
| 预测时域 Np | 25 | 20 |
| 控制时域 Nc | 10 | 8 |
| Q(状态权重) | diag([50, 50, 80, 10, 10, 10]) | diag([40, 40, 20]) |
| R(输入权重) | diag([0.1, 0.1]) | diag([0.5, 0.5]) |
| 约束方式 | 输入硬约束+状态软约束 | 输入硬约束+状态软约束 |
核心的CasADi优化问题构造如下:
opti = casadi.Opti(); % 决策变量:控制序列 U = opti.variable(nu, Np); % 状态轨迹,由神经网络预测得到 X = opti.variable(nx, Np + 1); opti.subject_to(X(:, 1) == x0); for k = 1:Np opti.subject_to(X(:, k+1) == x0 + dt * nn_predict_casadi(X(:, k), U(:, k), weight_struct)); end % 输入约束 for k = 1:Np opti.subject_to(u_min <= U(:, k) <= u_max); end % 轨迹跟踪代价函数 opti.minimize(sum(traj_cost));注意我在实现时不是直接用状态更新公式x(k+1) = f(x(k), u(k)),而是用了上面提到的神经网络预测增量再积分的做法,也就是后一步状态等于前一步状态加上预测增量。这样可以显著降低网络输出噪声带来的数值抖动。
调参顺序上有个很实用的经验:先固定采样时间和预测时域,把Q、R调到轨迹基本能跟踪上,再回来调Np。Np越大,MPC越有远见,但是计算量呈线性增长,而且模型误差在长时间预测里不断累积——特别是神经网络本身就是从数据里学的,长时间展开时误差会被放大,所以Np不宜盲目求大。我试过把机器人车的Np从20调到40,结果单步求解时间从0.08秒猛增到0.35秒,而且跟踪性能提升非常有限,显然是不划算的。
4. 复现中我踩过的坑:常见问题与排查实录
4.1 问题:NN-MPC一上来就发散,控制量剧烈震荡
这是我第一次联调四旋翼位置环时遇到的最头疼的问题。刚把神经网络预测模型接进MPC里,仿真一开始,控制量就出现幅度巨大的高频率抖振,位置轨迹直接飞掉。
排查之后发现根子在于:我训练神经网络时用的数据是从开环扫频激励得到的,但扫频信号的频率成分和幅值范围跟闭环MPC实际工作时产生的轨迹模式差异很大。也就是说,MPC在工作点附近查询的输入输出组合,正好落在训练数据分布比较稀疏的区域,神经网络在这个区域里给出的预测几乎是乱猜的。
解决办法有两个方向,我都试过。第一个方向是加大数据覆盖,在闭环MPC底部跑一遍基线控制器,把这个控制器的实际输入输出数据也混进训练集;第二个方向是在MPC代价函数里增加Δu的惩罚项,也就是对控制变化率做约束,防止控制量在相邻采样周期内剧烈跳变。实际采用两个方向叠加之后,震荡大幅缓解。
4.2 问题:CasADi重建网络预测时结果全错,符号变量类型不匹配
这个问题的症状特别诡异:单独调用网络预测函数时,输出结果跟训练好的net预测完全一致,但一放进MPC优化的目标函数里,结果就乱掉。我排查了很久,最后发现是我把矩阵运算的维度搞错了——Matlab里训练好的网络权重矩阵的维度,不是我想当然的那个。net.IW{1}的维度是当前层神经元个数乘输入维度,net.LW{2,1}才是第二层到第一层之间的权重,我最初把LW和IW索引搞反了,CasADi符号表达式看着没错,但实际数值初始化时维度不匹配。
建议在重建完神经网络之后,先做一步数值校验:随机取几个样本点,让CasADi函数和Matlab原始net分别计算输出,对比误差是不是在1e-8量级。这个校验一定要做,否则后面MPC优化出来的控制量会让你根本不知道是哪一环出的错。
4.3 问题:机器人汽车轨迹跟踪出现"侧滑"和原地打转
机器人车模型是非完整约束的,按理说不会侧滑,但我的MPC仿真里出现了类似侧滑的运动轨迹。症结在于参考轨迹本身不满足车的运动学约束——我随手画了一条参考路径,x_ref和y_ref的时间导数跟theta_ref不一致,导致车在跟踪时被迫产生一个不可能的横向速度分量。
解决方式是先把参考轨迹生成过程标准化:明确指定机器人车前向速度v_ref(t)和角速度ω_ref(t),然后用运动学方程积分出x_ref、y_ref、theta_ref。这样参考轨迹天然满足车辆的运动学约束,MPC跟踪起来就顺了。如果你一定要用任意路径做参考,那么至少要让theta_ref等于路径切线角,并且路径的曲率不能超过机器人车的角速度上限。
4.4 常见问题速查表
最后整理一个速查表,把我在这个复现项目里遇到、以及能被预测到的常见问题都列进去,方便你复现时对照排查:
| 现象 | 可能原因 | 排查与解法 |
|---|---|---|
| 初始阶段控制量震荡 | NN训练数据覆盖不足;Δu惩罚权重太小 | 混入闭环数据重新训练;增大S矩阵权重到1~5 |
| 位置轨迹发散但姿态正常 | 旋转矩阵推力方向写错或符号反号 | 单独用恒定输入做开环验证,检查推力是否指向正确的方向 |
| 优化求解时间过长 | Np过大;网络层数过多;CasADi表达式未做预编译 | 降低Np到20左右;隐藏层不超过两层;用casadi.Function预编译 |
| 机器人车原地打转 | 参考轨迹缺少航向角约束;theta_ref与位置导数不一致 | 用运动学模型生成参考轨迹;对航向误差增加罚项 |
| MPC偶尔无解 | 硬约束太紧,控制量限幅又很小 | 输入限幅保留安全边界;状态约束改软约束,代价函数加松弛罚项 |
| 网络预测单步误差很小但轨迹误差大 | 多步展开时误差累积 | 训练时引入多步预测损失;MPC滚动时使用反馈校正 |
| CasADi函数赋值错误 | 权重矩阵维度、索引混淆 | 固定用数值校验函数对比net和Casadi重建输出;在Matlab命令行单步验证 |
这个项目复现完之后,我对NN-MPC的理解已经有了质的提升。中间走了不少弯路,尤其是数据采集和网络嵌入这两块,每一步都真是“看着简单做起来坑”。我个人的体会是,这类算法项目的核心难点其实不在控制理论本身,而在如何让机器学习模型和控制优化器这两个不同语法的系统真正顺畅地协作起来。如果你刚开始复现这个方向,建议一定先把基线MPC用机理模型跑通,再一步步换成神经网络预测,每一步都验证,不要试图一把梭。
最后再分享一个小技巧:用CasADi的时候,把训练好的网络保存成结构体之后再编写一个独立的nn_eval函数,这个函数能在普通数值模式下运行,也能在符号模式下运行。每次调整权重或者归一化参数之后,先跑数值模式验证若干样本,再进MPC联调,能帮你节省至少两天的排查时间。