做控制的都知道一个道理模型预测控制MPC这东西成也模型败也模型。模型给得准滚动优化就是降维打击模型给得不准预测窗口越长系统越容易被带偏。这两年在做非线性系统控制项目时我把重心放到了数据驱动的MPC研究上——不依赖精确机理模型直接用历史数据驱动预测和优化仿真平台选的是Matlab。这套思路对非线性对象尤其有用因为机理模型往往难建、易漂移而数据是现成的。这篇文章把我在“非线性与数据驱动的模型预测控制MPC研究Matlab代码实现”这个项目里的建模思路、算法拆解和可运行代码一起整理出来重点讲清楚非线性MPC为什么难做、数据驱动MPC有哪些技术路线、Hankel矩阵和优化问题怎么构建、以及Matlab代码怎么落地。适合正在做控制算法仿真、毕业论文或者先进控制Demo的朋友参考。就算你之前没碰过MPC照着代码跑一遍也能把从数据到预测再到闭环控制的全链路摸清楚。1. MPC最大的坑模型不准后面全白费1.1 滚动优化再漂亮模型误差一票否决先说说MPC的基本机制。标准的模型预测控制在每个采样时刻都会基于当前状态和历史信息利用一个预测模型去推算未来N步的系统行为然后求解一个带约束的有限时域优化问题得到一组最优控制序列。这组序列只执行第一步下一个时刻再滚动优化一遍。这个“滚动”机制看着很聪明但它有个致命前提预测模型必须和真实系统足够接近。举个简单的数字例子。假设模型认为系统增益是2实际增益是3那么模型预测的未来输出就会和真实输出差出50%。如果控制器要跟踪一个目标值它会按照“预测输出偏低”来加大控制量结果实际输出就会超调甚至发散。预测时域N拉得越长模型误差带来的误差累积越严重。工业现场的经验是控制周期可以快一点、预测时域可以短一点但模型误差这道坎过不去什么都白搭。我在早期调试MPC的时候就遇到过机理模型在实验室验证很好、一到现场就跑偏的情况。最后排查来排查去发现是工艺参数漂移了模型里的传热系数和实际差了一大截。那段时间我最大的体会就是MPC项目最耗时的不是求解器调参而是模型的获取和维护。1.2 非线性系统建模的三重困境如果对象是线性系统情况还好办状态空间模型用辨识工具箱一把梭就行。但碰见非线性系统建模难度直接翻倍。第一重困境是工作点依赖。非线性系统在A工作点辨识出来的线性模型换个幅度或者换个工况增益和动态特性可能完全变样。最典型的是带有执行器饱和、摩擦死区或者化学反应强放热的过程对象。如果你用固定模型做MPC操作范围一旦超出设计区间控制器性能就是断崖式下跌。第二重困境是机理建模成本太高。工艺机理清楚的话还能建立一阶或二阶微分方程组但参数辨识需要大量实验设计而且不少参数比如反应活化能、传热系数很难直接测量。机理模型建出来之后还需要大量的修正表来补误差维护成本极高。第三重困境是模型老化。设备磨损、催化剂失活、环境温度变化都会让模型参数随运行时间发生变化。即使模型初始精度很好三个星期之后也可能不再可靠。也是这重重困境逼着大家开始认真看待数据驱动的MPC思路。既然精确模型难建、难维护干脆直接让数据说话——用系统的历史输入输出轨迹来表达未来行为把“建模”和“控制”放在同一个优化框架里解决。2. 数据驱动MPC的三条技术路线怎么选2.1 系统辨识加MPC工程上最稳妥的第一步数据驱动MPC最成熟的技术路线其实是“系统辨识MPC”的两步法。先用历史数据拟合一个显式模型常见的有ARX模型、状态空间模型、NARX非线性模型然后把模型塞进MPC框架里做预测控制。这条路线最大的优点是工具链完整Matlab里面System Identification Toolbox和Model Predictive Control Toolbox都可以直接拿来用对工程人员来说几乎没有门槛。我在项目里先做的就是这个方案用PRBS激励信号跑了一段开环实验拿数据拟合了一个二阶ARX模型然后搭线性MPC。系统在小信号范围内工作得很不错但把目标值拉高到非线性区间之后稳态误差就出来了。原因很简单ARX模型本质上是全局平均意义上的线性近似非线性系统的局部增益变化它表达不出来。不过作为第一步验证它的意义在于把MPC的整个链路先跑通数据采集、模型辨识、QP求解、闭环仿真每一步都有成熟工具支撑适合新手建立整体认知。2.2 进阶路线DeePC用数据直接做预测控制DeePCData-enabled Predictive Control数据赋能预测控制这两年热度很高它跟“辨识MPC”最大的区别在于不显式建模直接用历史输入输出轨迹构成Hankel矩阵把未来的预测表达成历史轨迹的线性组合。大白话讲就是把收集到的历史输入输出数据当成一本“字典”DeePC在优化的时候要从字典里找出一段和当前状态最匹配的轨迹通过对历史轨迹进行线性组合拼出一段符合控制目标的未来预测。整个过程把数据匹配、噪声补偿、未来输入求解、未来输出预测放在同一个优化问题里同时解决不需要先拟合出A、B、C、D矩阵。原始DeePC算法面向的是线性时不变系统理论依据是行为理论里的一个结论只要输入信号满足持续激励条件系统的所有可能轨迹都可以由历史轨迹的线性组合线性表示。但把它直接拿到非线性系统上效果会打折扣。所以我在代码里做了两个处理一是用滑动窗口只取最近一小段数据参与构建Hankel矩阵相当于局部线性化二是在闭环过程中不断把最新数据追加进历史库让“字典”始终覆盖当前工作点附近。2.3 非线性扩展局部窗口、在线更新和残差补偿对于非线性系统的数据驱动MPC可以按“成本从低到高”分成三个层次去扩展。最低成本的做法就是局部数据窗口这也是我主推的方法。既然全局线性近似覆盖不了非线性那就只用“当前工作点附近”的数据来构建Hankel矩阵。非线性系统在局部范围内总是近似线性的所以只要窗口够短数据表达能力就跟得上。窗口长度的选择有讲究太长会引入过多非线性失真太短则数据不足、优化问题欠定。中等成本的做法是在线更新历史库。每执行一步控制就把实测的输入输出追加到数据集中。这样即使系统参数缓慢漂移数据驱动控制器也能自适应跟上。历史库也不是无限增长的一般维持最近几百条数据就够用不然求解规模越来越大实时性会崩。更高成本的做法是用黑箱模型做残差补偿。比如先用线性数据驱动MPC做预测再用神经网络或者高斯过程回归去拟合线性模型预测残差把残差补偿项加进目标函数或约束里。这样对强非线性效果好但训练和调试成本都明显上升可解释性也会差一些。一般来说先试局部窗口和在线更新大多数场景已经够用。3. 核心算法拆解Hankel矩阵与QP优化3.1 Hankel矩阵到底怎么构造理解DeePC的关键是先理解Hankel矩阵的构造逻辑。假设过去窗口长度T_ini2预测时域N3数据序列长度为L7u [1; 2; 3; 4; 5; 6; 7]那么控制输入的Hankel矩阵可以拆成两块。上面T_ini行是过去窗口部分Up下面N行是未来窗口部分Uf。具体来说第一列Up取u(1:2)Uf取u(3:5)也就是[1;2]和[3;4;5] 第二列Up取u(2:3)Uf取u(4:6)也就是[2;3]和[4;5;6] 第三列Up取u(3:4)Uf取u(5:7)也就是[3;4]和[5;6;7]列数的计算公式是 K L - T_ini - N 1。这个公式推导起来很简单每列覆盖T_iniN个数据点相邻两列滑动一个数据点所以能滑出来的列数就是L减去窗口长度再加1。输出的Hankel矩阵构造方式和输入完全一样。在Matlab里可以用循环来构造因为内置的hankel函数对多通道数据处理并不方便尤其是要把同一份数据拆成Up和Uf两块时循环更直观也不容易出错。3.2 DeePC优化问题的标准形式DeePC在每个采样时刻要求解一个优化问题。决策变量分成四块数据组合系数g、初始匹配误差sigma、未来控制输入u_f、未来预测输出y_f。目标函数写成min ||g||^2 lambda_sigma ||sigma||^2 R ||u_f||^2 Qy ||y_f - y_ref||^2第一项是数据组合系数的正则项用来保证数值稳定性第二项是初始匹配误差的惩罚项用来吸收噪声和非线性偏差第三项和第四项是标准的MPC目标分别限制控制能量和跟踪误差。约束条件分成两组。第一组是初始条件匹配约束[Up; Yp] * g [u_ini; y_ini] sigma意思是投影到过去窗口上数据组合出来的轨迹必须和当前测量到的历史输入输出一致允许偏差sigma。第二组是未来轨迹一致性约束[Uf; Yf] * g [u_f; y_f]意思是未来预测完全由数据组合表达不依赖任何显式模型。再加上输入输出上下限约束整个优化问题就是一个标准的二次规划可以直接交给quadprog求解。这里有个参数需要重点说明lambda_sigma的取值。太小的话优化器会把非线性偏差全部归给sigmag就会失去预测能力太大的话初始匹配变成硬约束有测量噪声时容易无解。我在不同系统上试下来一般取1e3到1e6之间具体看信号幅值大小。3.3 T_ini、N和局部窗口Lw怎么定这几个参数是DeePC能否跑通的命门。T_ini是过去窗口长度理论上要大于系统的状态维度工程经验是取系统“明显记忆长度”的1.5到2倍。比如我后面代码里的系统是二阶差分方程记忆长度是2步那么T_ini取4就足够覆盖动态特性。如果T_ini取得太小初始匹配约束无法唯一确定当前状态预测会飘。N是预测时域取值逻辑和普通MPC一样要覆盖被控对象上升时间的60%到80%。取值太短会牺牲预见性太长则求解规模变大而且DeePC的Hankel矩阵列数会变少。Lw是构建Hankel矩阵的局部数据窗口长度。Lw决定列数KK必须远大于N否则优化问题缺乏自由度。我一般要求K至少是10倍的N也就是说Lw至少要在T_iniN的基础上再加10倍N。以N15为例Lw取200左右比较稳妥。4. Matlab代码实现非线性系统数据驱动MPC完整流程4.1 仿真对象选取整个项目里我用的仿真对象是这样一个二阶非线性差分方程x(k1) 0.9x(k) - 0.2x(k-1) 0.05sin(3x(k)) 0.4*u(k)输出测量为y(k) x(k) 测量噪声。这个系统有三个优点第一开环稳定控制器调试过程中不容易跑飞第二正弦非线性项在状态偏离零点后会让系统增益明显变化能真实反映出线性模型MPC的局限性第三二阶记忆长度考验数据驱动方法对动态的表达能力。控制任务是让系统输出跟踪一系列方波参考信号目标值在0.5、1.0、0.2、0.8之间切换。在参考值较大的区间正弦非线性造成的增益压缩效应会变得明显这时线性MPC模型的预测误差就会被暴露出来。4.2 离线数据生成代码离线数据是DeePC的“字典”必须包含足够丰富的信息。这里用PRBS叠加低频正弦作为激励信号PRBS保证幅值域覆盖低频正弦保证频域激励充分。%% 离线数据生成 rng(10); L 800; % 离线数据长度 u_hist 0.8*(2*(rand(L,1)0.5)-1) 0.3*sin(0.05*(1:L)); y_hist zeros(L,1); x 0; xp 0; for k 1:L xn 0.9*x - 0.2*xp 0.05*sin(3*x) 0.4*u_hist(k); xp x; x xn; y_hist(k) x 0.005*randn; % 测量噪声 end data_all struct(u, u_hist, y, y_hist);这段代码的逻辑很简单先生成激励信号然后跑一遍系统模型把输入输出记录保存到data_all结构体里。运行完之后data_all.u和data_all.y就是后面构建Hankel矩阵的数据源。4.3 DeePC求解函数求解函数是整套代码的核心。它接收当前历史的u_buf和y_buf长度T_ini、参考轨迹ref_seq、预测时域N以及目标权重输出最优控制序列。function u_opt dee_mpc_solve(data_all, u_buf, y_buf, ref_seq, N, Qy, R) T_ini length(u_buf); Lw 200; % 局部数据窗口长度 if length(data_all.u) Lw error(历史数据长度不足请增加离线数据或减小Lw); end u_d data_all.u(end - Lw 1 : end); y_d data_all.y(end - Lw 1 : end); K Lw - T_ini - N 1; % 构造Hankel矩阵 Up zeros(T_ini, K); Uf zeros(N, K); Yp zeros(T_ini, K); Yf zeros(N, K); for i 1:K Up(:,i) u_d(i : iT_ini-1); Uf(:,i) u_d(iT_ini : iT_iniN-1); Yp(:,i) y_d(i : iT_ini-1); Yf(:,i) y_d(iT_ini : iT_iniN-1); end % 决策变量 z [g; sigma; u_f; y_f] lambda_sigma 1e4; nz K 2*T_ini 2*N; H blkdiag(2*eye(K), 2*lambda_sigma*eye(2*T_ini), ... 2*R*eye(N), 2*Qy*eye(N)); f [zeros(K 2*T_ini N, 1); -2*Qy*ref_seq]; Aeq [Up, -eye(T_ini), zeros(T_ini,T_ini), zeros(T_ini,N), zeros(T_ini,N); Yp, zeros(T_ini,T_ini), -eye(T_ini), zeros(T_ini,N), zeros(T_ini,N); Uf, zeros(N,T_ini), zeros(N,T_ini), -eye(N), zeros(N,N); Yf, zeros(N,T_ini), zeros(N,T_ini), zeros(N,N), -eye(N)]; beq [u_buf; y_buf; zeros(N,1); zeros(N,1)]; % 输入输出上下限 umin -2; umax 2; ymin -3; ymax 3; Aineq [zeros(N,K2*T_ini), eye(N), zeros(N,N); zeros(N,K2*T_ini), -eye(N), zeros(N,N); zeros(N,K2*T_ini), zeros(N,N), eye(N); zeros(N,K2*T_ini), zeros(N,N), -eye(N)]; bineq [umax*ones(N,1); -umin*ones(N,1); ... ymax*ones(N,1); -ymin*ones(N,1)]; options optimoptions(quadprog, Display, off); z quadprog(H, f, Aineq, bineq, Aeq, beq, [], [], [], options); u_opt z(K 2*T_ini 1 : K 2*T_ini N); end这里要注意一个细节目标函数里的||y_f - ref_seq||^2展开后线性项是-2Qyref_seqy_f对应f向量的最后N个元素。quadprog默认求解的是0.5zHz f*z所以H矩阵里都乘了2。4.4 闭环主循环主循环负责把DeePC求解、系统模型更新、历史数据追加这三件事串联起来。重点在于每次求解完要把实际执行的控制量和测量输出追加进data_all让滑动窗口始终包含最新工况数据。%% 闭环仿真数据驱动DeePC Ts 0.1; Nsim 300; T_ini 4; N 15; Qy 10; R 0.1; % 参考轨迹 yref 0.2*ones(Nsim,1); yref(20:70) 0.5; yref(71:140) 1.0; yref(141:220) 0.2; yref(221:end) 0.8; % 预热系统到工作点附近 x 0; xp 0; u_buf zeros(T_ini,1); y_buf zeros(T_ini,1); for k 1:T_ini u_apply 0; xn 0.9*x - 0.2*xp 0.05*sin(3*x) 0.4*u_apply; xp x; x xn; u_buf(k) u_apply; y_buf(k) x 0.005*randn; data_all.u(end1) u_apply; data_all.y(end1) y_buf(k); end u_history zeros(Nsim,1); y_history zeros(Nsim,1); for k 1:Nsim ref_seq yref(k)*ones(N,1); u_opt dee_mpc_solve(data_all, u_buf, y_buf, ref_seq, N, Qy, R); u_apply u_opt(1); % 更新真实系统 xn 0.9*x - 0.2*xp 0.05*sin(3*x) 0.4*u_apply; xp x; x xn; y_measure x 0.005*randn; % 追加历史数据在线更新 data_all.u(end1) u_apply; data_all.y(end1) y_measure; % 更新缓冲区 u_buf [u_buf(2:end); u_apply]; y_buf [y_buf(2:end); y_measure]; u_history(k) u_apply; y_history(k) y_measure; end % 画图 t (1:Nsim)*Ts; figure; subplot(2,1,1); stairs(t, yref, --, LineWidth, 1.2); hold on; plot(t, y_history, LineWidth, 1.5); xlabel(时间/s); ylabel(输出); legend(参考,DeePC); grid on; subplot(2,1,2); stairs(t, u_history, LineWidth, 1.5); xlabel(时间/s); ylabel(控制量); grid on;把前面两个代码块和这个主循环合在一起就是一个完整的可运行DeePC仿真脚本。我实际跑下来的结果是DeePC在参考值切换后大概3到4个控制周期就能跟上稳态误差普遍在0.02以内。作为对照组用同样的离线数据拟合一个二阶ARX模型再做线性MPC在小目标值区间问题不大但在目标值为1.0或者0.8的区间稳态误差会明显变大这就是非线性增益变化给固定模型带来的伤害。4.5 对照实现线性模型MPC线性模型MPC的代码逻辑相对传统核心是状态空间模型加二次规划。先用离线数据辨识ARX模型data_id iddata(y_hist, u_hist, Ts); m_arx arx(data_id, [2 2 1]); [A_m, B_m, C_m, D_m] idssdata(m_arx);然后把预测方程写成矩阵形式。假设状态为x_mpc预测时域内状态序列可以表达为X Fx0 GU目标函数展开成关于U的二次型再用quadprog求解。这部分是标准操作网上资料很多我这里就不贴完整代码了。我想强调的是当你把两种控制器放在同一个非线性系统上对比时才能直观感受到“模型错误传递到预测”这件事有多要命。5. 调试实录DeePC常见问题与避坑指南5.1 典型症状和排查方向调试DeePC主循环的过程中我记录了几个高频问题整理成速查表放在下面症状可能原因建议排查方向控制量剧烈抖动R权重太小局部窗口Lw太小参考轨迹直接给大阶跃调大R增大Lw给参考轨迹加缓变斜坡预测输出和实际输出有明显偏差T_ini不够长覆盖不了系统记忆sigma惩罚太大导致过度吸收偏差增大T_ini调低lambda_sigma检查数据激励覆盖范围quadprog报不可行输入/输出约束太紧初始匹配约束和系统真实轨迹不兼容放宽约束边界增加历史数据检查T_ini和N是否匹配数据长度稳态跟踪误差大参考值超出了离线数据的覆盖范围补充目标值附近的数据在线更新历史库增大激励幅值求解越来越慢历史库无限增长Hankel矩阵列数过多限制data_all长度只保留最近Lw条数据定期清理5.2 关于激励、正则和局部窗口的几条经验lambda_sigma这个参数常规文档里只会告诉你“用来吸收噪声”但实际调试你会发现它对结果的影响非常敏感。我踩过的坑是第一次跑的时候把lambda_sigma调得特别大结果优化器为了保证初始匹配约束严格成立硬把g系数拉到一个很奇怪的值预测出来的轨迹倒是跟初始条件对上了但未来的预测完全失真。后来我把lambda_sigma从1e6降到1e4效果反而好了很多。经验是只要系统噪声不是特别大lambda_sigma往小里调一些给匹配误差留一点余地预测会更稳定。激励信号的数据覆盖范围是另一个容易被忽略的坑。如果离线数据全是在原点附近的小信号激励而参考轨迹突然要求系统跑到x1.0的地方DeePC的字典里根本找不到接近这个工作点的轨迹片段线性组合再聪明也拼不出现实动态。所以在生成离线数据的时候激励幅值要覆盖控制任务关心的整个操作区间。局部窗口Lw的长度选择与前两个参数都有关联。Lw太长非线性失真大Lw太短K太小优化问题欠定控制量会抖。我试过从100到500之间的不同取值综合来看200是一个比较舒服的中间值列数足够多时间跨度又不会大到让非线性效应平均化。5.3 走向更大规模非线性系统的三个扩展方向如果你的对象比这个仿真系统更复杂比如强耦合、高维、甚至含未建模动态这套基础DeePC代码就需要升级。第一个方向是Koopman算子。它的核心思想是把非线性系统提升到一个高维线性空间然后在那个空间里做数据驱动MPC。相当于先用坐标变换把系统“拉直”再复用经典的线性控制理论。代价是要额外维护一个提升函数库维度会变高。第二个方向是残差补偿。主控制器仍用线性DeePC同时训练一个神经网络或高斯过程模型去拟合DeePC预测和实际输出之间的残差把残差预测加到目标函数里。这个方案对强非线性效果不错但机器学习的调参经验又回来了。第三个方向是多模型调度。把操作空间划分成若干子区域在每个子区域分别建立Hankel矩阵控制器根据当前状态在工作点附近选矩阵或者做加权融合。这相当于把“局部窗口”的思想从时间域搬到了工作点空间域理论上更贴合强非线性场景。我个人做下来的体会是先从时间域局部窗口和在线更新开始绝大多数非线性过程都能压得住。真要到了压不住的那天再考虑Koopman或者学习型残差补偿也不迟。这个项目做完我最大的感受是数据驱动MPC真正的工程价值不是省掉了建模而是把“模型维护”变成了“数据维护”。你不需要花三个星期去校一套机理参数但需要认真设计激励实验、维护数据覆盖范围。Matlab里这套代码跑通之后换对象只需要改系统函数和激励信号参数性价比确实高。建议新手从这套代码入手先跑通闭环再慢慢改T_ini、N、Lw这几个参数观察它们对跟踪曲线和控制量的影响比看十篇论文都有用。