数据驱动和模型预测控制这两个词放在一起，很多人第一反应是“是不是以后不用建模了？给一堆数据就能做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矩阵，每块是[p+f]长的窗口，按时间滑动 N = length(data); numCols = N - p - f + 1; H = zeros(p + f, numCols); for k = 1:numCols H(:, k) = data(k:k+p+f-1); end end窗口长度p和f的选择有讲究。p至少要比系统的有效状态维度大，否则过去信息不够充分；f则决定了预测时域的长度。在模拟项目X里，我尝试过把p从5增加到15，系统的预测精度在p=10之后提升不明显，但矩阵条件数开始恶化。这涉及一个本质矛盾：窗口越长，历史信息越丰富，但数据矩阵的维度越高，对激励信号的要求也越高。激励不充分时，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}(k+i) ) 是由数据驱动预测器生成的预测输出，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秒，过去窗口长度p=10，预测时域Np=15，控制时域Nu=3。这个组合的逻辑是：预测时域覆盖约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.Nu+1: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作为一个快速原型验证环境，非常适合用来理清这套方法的每一个细节。希望我踩过的这些坑能帮你少走一些弯路。
