数据驱动和模型预测控制这两个词放在一起很多人第一反应是“是不是以后不用建模了给一堆数据就能做MPC了”说实话我第一次接触这个概念时也这么想过但真正动手在Matlab里跑通非线性系统的数据驱动MPC之后才意识到这里有个关键误区数据驱动并不是抛弃模型而是换一种方式获取模型——从输入输出数据中直接提取一个可用的预测模型再在这个模型上做滚动优化。这篇文章我打算完整讲一遍我在模拟项目X里实现非线性与数据驱动MPC的过程包括为什么在非线性系统上要用数据驱动方法、Hankel矩阵和子空间辨识是怎么回事、Matlab代码怎么组织、以及我实际跑实验时踩过的那些坑。如果你正在做基于Matlab的MPC仿真或者需要处理一个机理模型难写、非线性又明显的被控对象这篇文章应该能给你一个从原理到代码的完整参考。1. 为什么非线性MPC需要数据驱动模型获取的痛点1.1 机理建模在非线性系统上有多难传统MPC的核心是一个预测模型这个模型通常来自机理建模比如机械臂的拉格朗日方程、化工反应釜的能量守恒和质量守恒、飞行器的刚体动力学方程。机理模型的优点是物理意义清晰外推能力好但缺点也非常显著——在非线性系统上机理建模的难度会指数级上升。拿我在模拟项目X里用到的一个被控对象举例一个带有显著摩擦非线性和饱和特性的机电执行机构。理论推导可以得到它的标称模型但实际运行中摩擦力矩随温度漂移、负载变化带来的参数摄动、还有执行器死区这些因素会让机理模型和实际对象之间出现明显的失配。这种情况下一旦工况偏离标称点MPC预测的轨迹就会失真控制性能下降还是小事严重的会产生振荡甚至发散。还有一个更现实的场景很多系统的机理模型根本写不出来或者写出来之后复杂到无法用于实时优化。比如一个存在复杂传热过程的温控系统内部的热交换系数、热容参数很难精确标定你花了一个月推公式最后得到的模型精度可能还不如黑箱辨识的结果。这时候数据驱动方法的价值就体现出来了。1.2 线性化MPC的局限和线性参数时变方法的困境很多做工程的朋友会说非线性系统没关系我在工作点附近做线性化不就行了确实在工作点附近做线性MPC是工程上非常成熟的做法但它有两个硬伤首先是工作点切换的问题。一个化工过程可能需要在多个工况之间切换每个工况做一次线性化就意味着需要多套线性模型而且切换边界上如何平滑过渡是个很麻烦的问题。其次是强非线性问题比如磁滞、饱和、死区这类不连续或不可微特性线性化本身就会丢失关键的系统行为。线性参数时变LPV方法可以部分解决这个问题它把非线性系统看成一系列局部线性模型的插值但LPV模型的调度变量怎么选、子模型怎么划分、稳定性怎么保证在实际操作中每一步都是挑战。我在模拟项目X里最开始也试过LPV方案调度变量的选择就对最终控制效果影响很大稍不留神插值出来的模型就会失去对实际对象的代表性。1.3 数据驱动MPC到底改变了什么数据驱动MPC的思路是绕开“机理建模”这个大坑直接从输入输出数据中学习系统的预测模型。这个思路并不是新的子空间辨识方法从上世纪九十年代就开始发展了但近年来数据驱动MPC重新受到关注核心原因是两个一是被控对象的数据获取变得越来越容易传感器和执行器的成本大幅下降二是算力提升之后在线更新数据模型变得现实了。需要特别澄清的是数据驱动MPC不是“无模型控制”。它依然有模型只不过这个模型是数据构建出来的。我自己的理解是数据驱动是为了解决“模型从哪来”的问题而不是解决“要不要模型”的问题。很多文献里提到的数据驱动MPC本质上是用Hankel矩阵存储历史轨迹通过子空间辨识或持久激励数据直接构造预测器然后在预测器之上做滚动优化。2. 从高维轨迹数据到可用的预测器核心方法拆解2.1 一条时间序列里能挖出什么假设我们有一个单输入单输出的非线性系统采集到了一批输入输出数据。数据驱动MPC要做的第一件事是从这批数据里提取一个能预测未来输出的模型。最朴素的想法是直接建立一个回归模型用过去的输入输出去预测未来的输出。这个思想其实和ARX模型很像但在非线性系统上单纯增加过去数据的长度并不能显著提升预测精度——非线性系统的状态信息往往无法仅仅由有限长的输入输出窗口来近似表征。这也是为什么数据驱动MPC在实现上通常要借助子空间方法或动态模式分解这类工具而不是简单地做回归。我在模拟项目X里采用的是子空间辨识的做法利用过去的输入输出窗口构造回归量再通过最小二乘或SVD分解直接构建预测器。子空间方法的优势在于它不要求你对系统结构做太多假设只需要数据具备充分的激励性就能把系统的可观子空间和可控子空间信息提取出来。2.2 Hankel矩阵的构造看起来简单但细节不少Hankel矩阵是整个数据驱动MPC的载体。它的构造逻辑是把采集到的输入序列和输出序列按照指定窗口长度切成一块块“历史片段”每一行代表一个时间片段的轨迹信息行数足够多时这些行构成的子空间就能大致覆盖系统在激励信号下的动态行为范围。在Matlab里构造Hankel矩阵的代码非常短但要注意边界处理function H buildHankel(data, p, f) % data: 输入或输出序列列向量 % p: 过去窗口长度 % f: 未来窗口长度 % 返回: 块Hankel矩阵每块是[pf]长的窗口按时间滑动 N length(data); numCols N - p - f 1; H zeros(p f, numCols); for k 1:numCols H(:, k) data(k:kpf-1); end end窗口长度p和f的选择有讲究。p至少要比系统的有效状态维度大否则过去信息不够充分f则决定了预测时域的长度。在模拟项目X里我尝试过把p从5增加到15系统的预测精度在p10之后提升不明显但矩阵条件数开始恶化。这涉及一个本质矛盾窗口越长历史信息越丰富但数据矩阵的维度越高对激励信号的要求也越高。激励不充分时Hankel矩阵会秩亏提取出来的模型就会出现较大偏差。2.3 SVD降维从数据海洋里捞系统动态Hankel矩阵构造完成后接下来就是把它拆开找到系统的主要动态模式。这一步的核心数学工具是奇异值分解SVD。将过去的输入输出数据组合成回归矩阵将未来的输出数据作为目标矩阵通过最小二乘求取预测器的系数矩阵。在做最小二乘之前通常需要先对回归矩阵做SVD截断较小的奇异值目的是提升数值稳定性同时排除数据中噪声主导的方向。这一段我给出模拟项目X中实际使用的预测器构造代码function [Lp, Lf] buildPredictor(u_past, y_past, u_future, y_future, r) % 输入: 过去输入输出窗口矩阵、未来输入输出窗口矩阵 % r: 截断阶次根据奇异值能量保留 % 组合回归量过去输出和输入 Xi [y_past; u_past]; % 目标未来输出 Yf y_future; % SVD分解 [U, S, V] svd(Xi, econ); % 按能量保留前r个奇异值 singularValues diag(S); energy cumsum(singularValues.^2) / sum(singularValues.^2); r_actual find(energy 0.99, 1); if nargin 4 ~isempty(r) r_actual min(r, size(Xi, 1)); end % 截断后的伪逆 Xi_pinv V(:, 1:r_actual) * diag(1 ./ singularValues(1:r_actual)) * U(:, 1:r_actual); % 最小二乘得到预测器 L Yf * Xi_pinv; % 拆分为Lp作用于过去数据Lf作用于未来输入 Lp L(:, 1:size(y_past, 1)); Lf L(:, size(y_past, 1)1:end); end这里有个被我反复验证过的细节SVD截断阶次的选取直接影响预测质量。截断阶次太低会损失系统动态信息预测器表达力不足截断阶次太高又会把噪声子空间放进来导致预测结果在噪声影响下抖动。我通常的做法是看奇异值能量占比曲线找到一个“肘部位置”一般保留能量达到95%到99%的奇异值。但在处理强非线性系统时要小心非线性系统的动态往往分布在一个较宽的奇异值谱上强行截断到很低阶次会影响对剧烈动态的预测能力。3. Matlab实现链条离线辨识、在线滚动优化、闭环更新3.1 总体结构设计一个完整的数据驱动MPC仿真在Matlab里的实现大致分为四个环节数据采集、离线辨识预测器构造、在线滚动优化、闭环性能评估。这四个环节的代码结构如下% 环节1: 数据采集 —— 在设计好的激励信号下运行被控对象或高精度仿真模型 % 得到 u_data (输入), y_data (输出) % 环节2: 离线辨识打印 —— 构建Hankel矩阵SVD提取预测器 % 环节3: 在线滚动 —— 每个采样周期求解一个有限时域优化问题 % 环节4: 评估 —— 绘制跟踪曲线、控制量曲线统计误差指标在模拟项目X里我把环节2和环节3封装成了两个独立的脚本中间用.mat文件传递预测器参数。这样做的原因是调试起来更灵活如果你在做实时仿真的话建议把环节3单独做成一个函数方便在Simulink里以MATLAB Function模块的方式调用。3.2 数据采集阶段的激励信号设计数据驱动方法有个基本前提数据必须充分激励系统动态。我最初犯的错误就是为了图省事直接用了阶跃信号作为采集激励结果采回来的数据里系统的动态信息占很少SVD得到的模型在后续MPC预测时偏差很大。后来自查原因才明白问题所在阶跃信号只在一个方向上激励系统非线性系统的多个工作点根本没有充分覆盖。正确的做法是使用持续激励信号最常见的是伪随机二进制信号PRBS或者幅值逐渐变化的扫频信号。PRBS在高层次上相当于把一个幅值在正负之间不断切换的二进制信号输入系统切换间隔需要根据系统的时间常数来选择。切换太频繁系统响应跟不上采集到的输出几乎是滤波后的噪声有效信息不多切换太慢则数据大部分是稳态响应动态信息不足。以模拟项目X的机电系统为例系统时间常数约为0.3秒我最终采用的PRBS采样周期是0.05秒保持时间设为6个采样周期这样能在30秒的采集时间内覆盖足够的频率区间。% PRBS激励信号生成示例 N 600; % 采样点数 Ts 0.05; % 采样周期 clockPeriod 6; % PRBS保持周期采样周期整数倍 amplitude 1.0; u_prbs idinput(N, prbs, [0, 1/clockPeriod], [-amplitude amplitude]); u_data u_prbs; % 将被控对象在u_data激励下运行记录y_data注意在采集数据时还要考虑初始瞬态的影响。我的做法是预留一段“预热时间”先把系统运行到稳态之后再开始记录数据否则初始状态带来的瞬态响应会被当成系统动态的一部分后续预测会引入偏差。3.3 在线优化问题的数学形式与Matlab表达预测器构建好之后MPC的核心就变成了一个滚动优化问题。在每个采样时刻控制器根据当前已知的过去输入输出信息预测未来一段时域内的系统输出然后求解一个优化问题得到最优控制序列并且只执行第一步下一时刻重新求解。优化问题的数学形式如下在这个问题中( \hat{y}(ki) ) 是由数据驱动预测器生成的预测输出r是参考轨迹Δu是控制增量。第一项衡量输出跟踪性能第二项惩罚控制动作的剧烈变化最后一项是终端代价项用来近似表征无限时域的性能在非线性系统上它对稳定性有实际作用。在Matlab里我使用fmincon来求解这个优化问题。虽然fmincon是个通用约束优化求解器在线性MPC场景里通常会抱怨它速度慢但在非线性数据驱动MPC里由于预测器本身可能是非线性的比如使用了核方法或者局部加权模型fmincon这种通用求解器反而成了最稳妥的选择。% 滚动优化核心代码单步示意 function u_opt solveMPC(Lp, Lf, X_past, U_future_ref, params) % X_past: 过去窗口的输入输出组合 % params: 包含权重、约束、预测时域等 Np params.Np; Nu params.Nu; % 决策变量未来控制序列 u0 zeros(Nu, 1); % 约束控制量限幅 lb repmat(params.u_min, Nu, 1); ub repmat(params.u_max, Nu, 1); % 定义代价函数 costFun (u) mpcCostFunction(u, Lp, Lf, X_past, params); % 求解 options optimoptions(fmincon, Display, off, ... Algorithm, sqp, MaxIterations, 100); u_opt fmincon(costFun, u0, [], [], [], [], lb, ub, [], options); end代价函数内部做的事情是用决策变量u下一时刻的控制序列和当前X_past一起通过预测器Lp、Lf计算未来输出预测再和参考轨迹比较累加跟踪误差和控制增量惩罚。这里面有个重要实现细节预测器Lp作用于过去窗口输出Lf作用于未来输入但两者在预测器内部是纠缠在一起的所以每次求解时都要根据当前测量值重新构造X_past然后调用预测器。3.4 为什么要在代价函数里加控制增量惩罚加控制增量惩罚有什么实际意义我一开始做数据驱动MPC时没有加这一项结果控制量出现了明显的高频抖动。原因在于数据驱动预测器本质上是从有限数据里提取的动态特性对高频成分的预测并不准而MPC优化器发现“用高频控制量可以更好地匹配参考轨迹”时它就会倾向于输出高频控制信号。控制增量惩罚的作用是给这些高频分量一个代价让优化器在跟踪精度和控制平滑度之间做权衡。实际的调参经验是控制增量权重( R_{\Delta u} )太小控制量抖动明显太大系统响应变慢。我在模拟项目X里的做法是先设一个中等值比如输出权重设为1控制增量权重设为0.1然后看仿真结果如果控制量有抖动就把权重以0.05步长逐步增大直到抖动消失且跟踪性能还能接受。4. 实测中的调参与踩坑记录从失效到稳定4.1 变量缩放数据驱动MPC的隐形陷阱这个问题在机理模型MPC里也存在但在数据驱动MPC里会被放大因为预测器是直接从数据中计算的。如果输入和输出的数量级差很大比如输出是温度几百摄氏度控制量是电压0到10V那么SVD分解时大数量级的变量会在矩阵中占主导小数量级的变量的动态信息容易被淹没。解决办法是在构造Hankel矩阵之前对数据做归一化处理。我在模拟项目X里的做法是对输入输出分别减去均值、除以标准差然后在预测器构造完成后再把缩放系数记录下来在线运行时对预测输出做反向缩放。% 数据归一化 u_mean mean(u_data); u_std std(u_data); y_mean mean(y_data); y_std std(y_data); u_norm (u_data - u_mean) / u_std; y_norm (y_data - y_mean) / y_std; % 预测器在归一化数据上构造 % 在线使用时 % y_pred_real y_pred_norm * y_std y_mean;这个改动看起来不起眼但实际效果非常明显。在未做归一化时我的数据驱动MPC在跟踪阶跃参考时经常出现稳态偏差排查了很久发现是预测器模型的数值条件太差导致的。归一化之后同样的数据和参数系统变得稳定多了。4.2 预测时域、窗口长度和正则化之间的耦合关系在调MPC参数时预测时域( N_p )、过去窗口长度p、未来输入窗口长度Nu这三者之间有很强的耦合关系。一个典型问题是如果预测时域太长预测误差会累积导致优化问题为了“弥补远期误差”而输出激进的控制动作如果预测时域太短系统来不及响应跟踪效果差。我在模拟项目X里最终选择的参数组合是采样周期0.05秒过去窗口长度p10预测时域Np15控制时域Nu3。这个组合的逻辑是预测时域覆盖约0.75秒大约是系统时间常数的2.5倍既能覆盖主要的动态响应过程又不会因为预测过远而让误差主导优化。还有一个容易被忽略的细节——数据矩阵的正则化。实际采集的数据总会有噪声Hankel矩阵的条件数可能很大。在构建预测器时我会在最小二乘求解中加入一个小的正则化项防止噪声引起的系数矩阵爆炸。这个做法的本质类似于岭回归系数很小但能显著提升预测器的稳定性。4.3 约束处理非线性系统上硬约束会咬人数据驱动MPC的优化问题里控制量约束通常比较好处理——直接作为决策变量的上下界。但输出约束比如系统温度不能超过某个上限处理起来要特别小心。原因是数据驱动预测器对输出的预测存在偏差如果在预测中还强制输出约束严格满足优化器可能会为了满足一个本身就不准确的约束而给出奇怪的控制序列。我在模拟项目X里遇到过这种情况给输出加了一个上限约束结果优化器在接近上限时控制动作出现剧烈的不连续变化系统实际输出反而发生了超调。后来我把约束改成了软约束形式——在代价函数里加一个对输出越界的惩罚项而不是硬性约束。这样优化器在权衡之后如果越界的代价小于强行约束带来的控制恶化代价它允许适当越界同时整体性能更好。软约束的实现是在代价函数里加入一个越界惩罚function J mpcCostFunction(u, Lp, Lf, X_past, params) % 计算预测输出 y_pred predictOutput(Lp, Lf, X_past, u); % 跟踪误差 trackError y_pred - params.reference; % 控制增量 u_prev X_past(end-params.Nu1:end, 1); % 简化示意 deltaU [u(1) - u_prev(end); diff(u)]; % 输出软约束惩罚 softPenalty params.rho * sum(max(0, y_pred - params.y_max).^2) ... params.rho * sum(max(0, params.y_min - y_pred).^2); % 总代价 J sum(trackError.^2) params.Rdu * sum(deltaU.^2) softPenalty; end这个rho系数的选取也很有意思。rho太小约束起不到作用rho太大等价于硬约束又回到了前面说的问题。我的经验是从一个中间值开始比如输出权重相同的量级然后看系统在接近约束边界时的行为来调整。4.4 终端代价让数据驱动MPC更有稳定性保证前面提过终端代价项很多做线性MPC的人会忽略它因为线性MPC可以通过LQR求解一个解析的终端代价矩阵。但在非线性数据驱动MPC里稳定性的理论分析更加困难终端代价的作用就变得更加实际——它能告诉优化器“眼光放远一点不要只顾眼前几步”。数据驱动场景下做终端代价的常用办法是在预测末端给一个基于局部线性化模型或者经验模型的终端惩罚。虽然不如理论上那么漂亮但在实际仿真中能明显改善闭环响应。我在模拟项目X里采用的方案是用数据驱动预测器的末端状态乘以一个常数矩阵作为终端惩罚。这个常数矩阵我是在离线时通过对平衡点附近的线性近似模型求解黎卡提方程得到的。5. 进阶方向与扩展思考从离线到在线5.1 在线模型更新的价值与代价离线数据驱动MPC有个天然的局限如果系统特性随运行时间发生变化比如设备老化、环境温度变化离线的预测器就会逐渐失配。解决思路是让模型能够在线更新——每隔一段时间利用新采样的输入输出数据重新训练或更新预测器。在Matlab里实现时更新方式有两种一是每次更新都完全重新构造Hankel矩阵和预测器计算量较大适合采样周期比较长的系统二是递归更新SVD的低秩分量计算量小但实现复杂度高。对于大多数仿真和原型验证来说第一种方式已经足够。我在模拟项目X里设置了一个简单的模型更新策略每30秒用滚动窗口内最近的数据重新构造一次预测器结果显示当模拟的参数摄动发生之后控制性能在1到2个更新周期内就能恢复。5.2 从单变量系统到多变量系统的扩展数据驱动MPC的优势在多变量系统上体现得更明显。多变量系统的机理建模非常繁琐耦合关系难描述但数据驱动方法只需要把输入输出矩阵的维度扩展就行。Hankel矩阵的构造逻辑不变SVD分解同样适用代价函数里把输出跟踪误差改为向量形式即可。当然多变量系统的数据需求会显著增加特别是各输入通道之间要设计相互独立的激励信号否则辨识出来的模型无法区分各通道的独立影响。我建议在采集数据时对各输入通道采用不同的PRBS序列并且错开时间偏移这样数据信息量会大很多。5.3 数据驱动MPC和扰动观测器的配合数据驱动MPC本身没有显式地处理外部扰动但实际系统里扰动是常态。一个实用的扩展是在MPC外层加一个扰动观测器估计系统受到的等效扰动然后把扰动估计值叠加到预测器的输出上。这种做法的本质是“数据驱动模型扰动补偿”的组合在工程上比单纯依赖数据驱动模型去拟合所有扰动要稳健得多。我在模拟项目X里给系统加了一个正弦扰动纯数据驱动MPC的跟踪误差明显增大加了扰动观测器之后误差恢复到无扰动时的水平。这个扩展只需要在MPC代价函数里把参考轨迹做一个前馈修正实现成本很低但收益非常明显。5.4 最后给新手的一些建议如果你正准备在Matlab里开始搭建自己的数据驱动MPC我的建议是不要一上来就追求复杂的算法。先把最简单的流程跑通用高精度仿真模型替代真实对象采一批充分激励的数据构建Hankel矩阵SVD得到预测器然后在一个简单的跟踪任务上验证预测精度。预测精度达标之后再把它和fmincon联起来做滚动优化。我走过的弯路包括数据激励不充分就急着做预测器、忘记归一化导致数值问题、硬约束导致优化器行为异常、没有加控制增量惩罚导致控制抖动。这几个问题每一个都会让结果看起来很糟糕但只要前置检查和参数调整做到位数据驱动MPC在非线性系统上的表现其实可以非常稳定。总结一下数据驱动MPC的价值在于给非线性系统的控制问题提供了一条绕过机理建模的务实路径但它对数据质量和数值处理的要求比想象中更高。Matlab作为一个快速原型验证环境非常适合用来理清这套方法的每一个细节。希望我踩过的这些坑能帮你少走一些弯路。