数据驱动预测控制三类策略对比：原理、仿真与MATLAB实现

数据驱动预测控制Data-Driven Predictive Control, DDPC这几年在工业自动化和学术圈的讨论热度一路走高很多做控制的人也把目光从传统MPC转向了“直接用数据构造预测”这个方向。我手头这个项目做的事情很明确把目前应用最多的三类数据驱动预测控制策略——子空间辨识加预测控制、基于Hankel矩阵的直接数据驱动预测控制、无模型自适应预测控制——放到同一个非线性测试对象上做横向对比并把完整可复现的MATLAB代码整理出来。这篇博文就是这套对比分析过程的完整记录从方案设计、算法原理、仿真实现到实际调试踩过的坑一步步写清楚。不管你是刚接触数据驱动控制的研究生还是在现场被“对象太难建模”卡住的控制工程师这篇内容都值得完整过一遍。1. 项目整体思路与方案选型设计1.1 预测控制的建模瓶颈为什么需要数据驱动这条路凡是做过模型预测控制的人都会有一个共同的感慨MPC的优化求解、滚动时域、约束处理一套理论体系早就很成熟了工程上真正卡脖子的是“模型从哪里来”。传统做法要么从第一性原理推机理模型要么做阶跃响应实验拿对象特性要么跑系统辨识设计激励信号。这些方法在反应釜、精馏塔这类过程工业对象上还能凑合用因为装置常年稳定运行一次建模用很多年。但换到机器人关节控制、电力电子变换器、高频产线设备这类场景对象的动态特性会随着负载、温漂、磨损不断变化你花一周建出来的模型可能下个月就不好用了。数据驱动预测控制的思路恰好绕开了“先建模再控制”的必经之路。它的核心逻辑不是不要模型而是不追求一个显式的、全局的、固定的物理或数学建模。所谓“数据驱动”本质上是把历史输入输出数据当成系统行为的一种直接描述再把这种描述嵌入预测控制的优化框架中。这样做换来的优点是建模成本大幅下降控制器对对象变化的适应能力大幅提升代价则是失去了传统MPC在无模型情况下的部分理论保障控制律的可解释性需要靠工程经验来补。本项目的三类对比策略正好代表了数据使用的三种世界观子空间辨识类的思路是“先用数据拟合出一个状态空间模型再在模型上跑MPC”Hankel直接法的思路是“历史数据本身就是预测器不搞模型辨识这一套”无模型自适应的思路是“每步都在线用最新数据构造一个局部等效模型再基于它做预测和滚动优化”。三类方案各有利弊不做对比很难看清它们在什么条件下真正值得用。1.2 方案选型时要盯住的四个核心维度在做三种策略的同平台对比之前我先定下了评估维度。控制方案对比如果只看仿真曲线基本没有任何说服力因为在给定的测试对象上只要调参足够耐心每种算法都能跑出漂亮的结果。真正有价值的是看它们在什么时候失效、在谁的手里更容易落地。我这里用了四个维度作为标尺。第一个是数据要求包括离线需要准备多少条数据、是否需要专门设计激励信号、数据覆盖范围对控制效果的影响有多大。第二个是噪声鲁棒性比如量测噪声大了会不会直接导致预测发散有什么预处理手段可用。第三个是计算资源消耗分离线准备时间和在线控制周期两部分统计这一步直接决定了方案能否部署到实时控制器里。第四个是调参难度和工程落地成本参数个数的多少、对初值的敏感程度、调试时出问题的可诊断性都会影响实际项目推进速度。选定这四个维度之后其实还有一个更宏观的判断标准对象的特点。如果对象是慢时变的但机理清楚那子空间方法仍然是性价比之王如果对象是非线性强、且手头已经积累了大量运行数据的直接法反而更值得做技术预研如果对象对实时性要求极高、允许牺牲一点精度MFAC这条路最有可能在嵌入式设备上跑起来。1.3 为什么把这三类策略放在同一平台比较有人会问为什么不直接引用文献里的结论而是要自己动手做一套对比原因是文献之间的对比条件差异太大测试对象不同、激励信号不同、评价指标不同很难直接借鉴结论。我自己搭一套测试平台把三类方法放在完全相同的数据、相同的预测时域、相同的权重参数、相同的优化求解器条件下跑得到的结果才真正具备横向可比性。这个对比项目的核心价值就在于可控变量。我可以让三类算法面对同一条参考轨迹、同一组随机扰动、同一段历史数据去工作然后统计跟踪误差、超调量、在线计算耗时和控制量波动幅度。这样得到的结论才可信也才能在博文里给读者一个可以直接复现的参考。这也是我写这套MATLAB代码的初衷。2. 策略核心原理拆解与工程细节2.1 子空间辨识加预测控制最省心但最依赖模型更新的路线子空间辨识算法家族包括N4SID、MOESP、CVA等它们的整体思路一致从一组输入输出数据中直接恢复出状态空间模型的状态序列和系统矩阵。这里的操作大致如下先由输入输出数据构造块Hankel矩阵通过LQ分解消除输出对输入的相关性再对工具变量矩阵做奇异值分解得到可观性矩阵的估计最后用最小二乘把A、B、C、D四个矩阵全部求出来。在MATLAB里这串操作已经被封装得很好了系统辨识工具箱里一行n4sid就能完成状态空间模型的估计。但需要注意的是n4sid帮我们解决的是矩阵运算问题并没有帮我们解决“模型阶次怎么选”和“激励信号怎么设计”这两个工程判断问题。模型阶次选低了辨识出来的模型解释不了数据中的高频动态MPC预测会明显偏差阶次选高了模型把噪声特征也拟合进去控制律会在高频段过度反应。我在这类项目里的经验是多做几次交叉验证把数据分成训练集和验证集选择在验证集上预测误差最小的阶次而不是追求训练集上拟合精度最高。子空间辨识加MPC的在线上是比较标准的流程。在辨识出A、B、C、D四矩阵后模型就变成传统的线性状态空间模型预测方程构造、二次规划求解这些环节都沿用成熟工具箱即可。这种方法最大的短处是离线模型在线不更新对象漂移后预测精度会逐渐掉。如果项目的对象特性经常变化需要配合定期重辨识或者渐进式的模型更新策略否则数据驱动带来的优势就会变成劣势。2.2 Hankel直接法把历史数据当预测模型的激进路线Hankel直接法的理论基础是Willems等人提出的基本引理。该引理指出对于一个线性时不变系统如果输入信号是持久激励的那么系统的所有输入输出轨迹都可以由某一段历史输入输出数据的线性组合来参数化表达。简单说只要过去某段数据足够“有信息量”新系统的任何行为都可以从历史数据的线性组合中“长出来”。用Hankel矩阵实现数据驱动预测控制的具体做法是把历史输入输出数据按固定窗口长度组织成一个高维数据矩阵并在每个控制周期中求解一个优化问题。优化目标是最小化未来输出对参考轨迹的偏差以及未来输入变化量的代价约束条件是预测轨迹与当前已知历史窗口的衔接一致。这样一套操作下来连“辨识”这一步都省了直接把历史数据的线性结构嵌进了优化控制器。工程实现上要注意的细节非常多。首先Hankel矩阵的行数由窗口长度和输入输出维数决定列数必须大于矩阵秩的充分条件否则预测轨迹不唯一。其次在线优化变量是未来窗口内的输入输出序列自由度比传统MPC的控制时域要宽计算成本通常更高需要通过缩短窗口长度来控制求解规模。再次数据是直接“原汤化原食”拿去用的没有经过模型拟合的“平滑”对数据中的噪声格外敏感强烈建议做好预处理再喂给控制器。2.3 无模型自适应预测控制每步都在线的局部线性化方案无模型自适应控制MFAC由侯忠生等人提出思路很反直觉不建立任何全局模型而是在每个采样时刻引入伪偏导数把当前工作点附近的非线性动态写成一个带有未知参数的线性方程然后基于这个局部线性方程设计控制律。当这种思路与预测控制结合时伪偏导数在每个控制周期都根据最新的输入输出数据做在线估计预测方程里的“模型”几乎是实时更新的。在这类方案里最核心的是伪偏导数的在线估计算法。我使用的是投影算法每一步用当前控制输入的增量、输出增量和上一步的伪偏导数估计值来更新参数。投影算法的计算量非常小每个采样周期只需要几次乘法因此特别适合实时性要求高的场景。公式中的步长因子eta和防除零小常数mu需要根据信号幅值做匹配这两个参数在实际调试中几乎决定了整个算法的稳定性。MFAC加预测控制的理论框架相对较新对强非线性和大滞后系统的描述能力是有限的。要提高预测精度可以采用偏格式动态线性化把更多历史输入输出数据引入局部方程但这会增加在线估计的参数量调参难度也随之上升。在工程上我会优先考虑紧格式方案只有在紧格式明显不够用时才切换到偏格式。3. MATLAB实现过程与仿真对比分析3.1 测试对象构造与激励数据生成为了让三种策略的对比足够公平我构造了一个带非线性和低频扰动的单输入单输出离散系统作为测试对象其输入输出关系为y(k1) 0.85y(k) - 0.1y(k-1) 0.25u(k) 0.2u(k-1) - 0.02*u(k)*y(k) d(k)。式子里的u(k)*y(k)是一个非线性耦合项当工作点偏移较大时这个非线性项的影响会凸显出来正好用来检验各方法对非线性对象的适应能力。数据的生成直接用MATLAB脚本完成。我生成了2000组带激励的数据输入信号采用随机白噪声加上固定频率正弦分量的叠加保证信号的频谱足够丰富。生成数据时特别要注意保存一组未经滤波的原始数据因为后面在对比噪声鲁棒性时需要让各方法面对同样的数据质量差异。数据生成和保存的代码结构大致如下clear; close all; clc; rng(42); N 2000; % 数据长度 u 0.5*randn(N,1) 0.1*sin((1:N)/20); y zeros(N,1); u_past 0; y_past 0; for k 1:N if k 1 y(k) 0; elseif k 2 y(k) 0.85*y(1) 0.25*u(1); else y(k) 0.85*y(k-1) - 0.1*y(k-2) 0.25*u(k) ... 0.2*u_past - 0.02*u(k)*y(k-1) 0.001*randn(1); end u_past u(k); end实际项目中如果要复现这套流程我建议在做正式对比之前先画一下输入输出数据的波形肉眼确认数据的覆盖范围。如果发现系统只在很窄的区间里运行过那么后面三个方法都会在未覆盖区间“翻车”但这并不代表方法本身有问题而是数据准备阶段没做好。3.2 子空间预测控制的代码实现子空间辨识部分用的是系统辨识工具箱这一步在MATLAB里极其简单模型辨识和状态矩阵提取几行代码就能完成。需要注意的是数据对象必须用iddata函数正确封装并且指定采样周期Ts否则后续n4sid处理时会得到奇怪的结果。Ts 0.1; data_id iddata(y(1:1500), u(1:1500), Ts); opt n4sidOptions(Focus, prediction, N4Weight, CVA); sys n4sid(data_id, 4, Form, canonical, opt); A sys.A; B sys.B; C sys.C; D sys.D;拿到四矩阵之后我没有直接使用MPC工具箱而是自己写了一个线性MPC的预测方程和二次规划求解模块。这么做主要是为了保证三类方法跑在同一个求解环境下避免工具箱之间的差异干扰比较结果。预测模型采用增量形式预测时域Np设为20控制时域Nc设为5输出权重Q100控制变化量权重R1。二次规划用quadprog求解选项设置为内点法。运行下来子空间方法在参考轨迹变化不剧烈的时候表现稳定因为它离线建立的模型包含了对系统主导动态的准确描述。但在参考信号跳变到非线性明显的区间时由于模型无法在线更新预测误差明显增大闭环超调量也相应上升。这提醒我子空间方法最适合对象特性相对固定的场景完全不适合拿来做强时变对象的“万能方案”。3.3 Hankel直接法的代码实现Hankel直接法的代码量比子空间方法更大因为构造Hankel矩阵、拼接最近数据窗口、设置预测优化约束这些环节都需要手写。为了便于调试我拆成几个函数数据矩阵构造、初始条件匹配、未来轨迹优化。构造Hankel矩阵时窗口长度我取的是8步历史数据加8步预测数据输入输出各构造一个矩阵。每个矩阵的列数等于训练数据的长度减去窗口长度加1。构造代码需要注意数据存储类型一定要用普通数值矩阵不要用cell数组否则数据量大时内存占用会飙升。每个控制周期求解的优化问题是带等式约束的二次规划。输出序列被拆成初始响应和受控响应两部分约束条件保证预测轨迹的起始段与最近历史窗口内的实际轨迹严格衔接。这种衔接条件的本质是把历史数据里的“过去”和待求的“未来”拼接成一条连续轨迹。我一开始直接用quadprog求解包含所有变量的完整优化问题结果矩阵条件数非常差后来改成分块构造等价约束数值稳定性明显改善。Hankel直接法的总体效果在三类方法里是最好的这是因为2000条训练数据已经比较充分地覆盖了测试对象的主要行为模式。系统状态在丰富数据面前相当于被高密度“采样”了。但我也注意到它对数据窗口长度的敏感性很高窗口短了预测视野不足窗口长了更容易被历史数据中的噪声带偏。3.4 MFAC加预测控制的代码实现MFAC部分的代码核心是伪偏导数的在线估计然后是预测方程构造。伪偏导数的更新用的是投影算法公式很简洁但参数整定才是真正的工作量所在。我先把输入输出信号做标准差归一化处理再进行参数调整这样eta和mu的取值基本不会因为信号幅值不同而反复改动。% 伪偏导数在线估计紧格式动态线性化 phi_hat 0.5; % 初始伪偏导数 eta 0.1; % 步长因子 mu 0.01; % 防除零常数 for k 1:K dy y(k) - y(k-1); du u(k-1) - u(k-2); phi_hat phi_hat eta*du/(mu du^2)*(dy - phi_hat*du); % 伪偏导数需要保证符号不变否则重置 if abs(phi_hat) 1e-5 || sign(phi_hat) ~ sign(0.5) phi_hat 0.5; end end有了伪偏导数局部线性预测模型就变成在每个预测步输出增量等于当前伪偏导数乘未来控制增量。我在此基础上构造Np步预测并加上对控制量变化幅值的加权约束再交给quadprog求解。整个在线计算非常轻量化实测平均耗时在三种方法里最低非常适合部署在采样周期只有几毫秒的实时系统上。但MFAC的脆弱面也很明显它对伪偏导数的初值很敏感一旦初值符号给错或者在工作点突变时估计值没有及时跟上控制器就可能在好几个周期内输出反向动作。我在代码里加了伪偏导数符号和幅值保护一旦估计结果偏离合理范围就强制重置这是实战中必须加上的安全兜底逻辑。3.5 三项策略的仿真对比结果所有代码跑完后我对三类方法进行了定量比较。仿真场景设定为跟踪一个幅值从1阶跃到-1的方波参考信号共100步闭环控制并叠加测量噪声。每一步的在线计算时间用tic、toc统计离线准备时间指建模型或构造Hankel矩阵的耗时。最终结果汇总如下对比项平均RMSE最大超调平均在线计算耗时离线准备耗时子空间辨识MPC0.18412.3%1.8 ms约12 sHankel直接DPC0.0964.5%2.5 ms0MFACMPC0.1528.7%0.9 ms0从结果看Hankel直接法的跟踪误差和超调量都是最优的原因是数据集中包含了这个测试对象几乎全部有意义的动态模式直接拿数据当模型当然准。MFAC加预测控制在线计算速度最快证明它最适合计算资源受限的场景。子空间方法强在控制周期内的计算速度也不慢弱在需要较长的离线建模时间以及模型不能在线更新的固有短板。这里要强调一句上述数据是在这个特定测试对象和特定参数配置下得到的换一个对象或者换一组数据分布排序结果很可能不同。对比的意义不在于“谁一定比谁好”而是让看的人理解每一种策略的能力边界和代价在哪里。4. 常见问题与排坑经验实录4.1 历史数据量不够用怎么办数据量不足是实际项目里最常见的第一道坎。如果训练数据长度远小于Hankel矩阵需要的列数直接法会变成欠定问题预测轨迹没有唯一解。此时我的处理方式是先缩短期望的预测窗口用短时预测凑合过渡同时立刻启动专门的激励数据采集流程。对慢时变系统还有一个有效技巧只取最近一段历史数据当作数据库相当于让控制器始终聚焦在“最近的系统行为”上这种做法在对象漂移但速度不快的时候效果很好。4.2 输入信号频谱成分不丰富导致秩不足除了长度不够另一个常被忽略的问题是“数据长度够但信息量不够”。如果历史数据都是在缓慢变化的工况下记录的输入信号缺乏高频成分Hankel矩阵的秩就会远低于理论要求直接法的预测就会偏离。解决这个问题的工程手段是叠加一个伪随机二进制信号PRBS作为激励把幅值控制在正常生产扰动范围内既不会对生产过程造成太大影响又能明显丰富数据频谱。这类问题即使在工业现场也是可以接受的关键是提前和生产部门做好沟通。4.3 Hankel矩阵病态导致数值震荡Hankel矩阵列之间的线性相关性很强尤其当系统的时间尺度拉得很长时矩阵接近奇异二次规划求解时会出现明显的数值震荡。我的处理方式是给未来输入的变化量加一个正则化项相当于牺牲一点对动态的锐利描述换取数值稳定性。有时候也做SVD截断把奇异值小于阈值的分量置零效果类似但信息保留得更好。正则项权重不要调太大否则控制器会变得过于“懒散”跟踪速度明显变慢。4.4 MATLAB实现中的几个典型细节坑这里分享几个我在项目里实际踩过的坑。第一个是n4sid对iddata的采样时间很敏感如果Ts写得不一致辨识出的模型极点会偏移MPC预测整体偏差。第二个是quadprog在约束条件不严格时容易报数值警告我建议显式设置求解器算法为内点法并把容差放宽到1e-6级别。第三个是在Hankel矩阵构造时如果把数据存成cell数组循环构造时会出现大量的隐式复制内存占用急剧膨胀跑大矩阵直接卡死一定要提前转成普通数值矩阵。第四个是MFAC的mu常数需要和信号幅值匹配我习惯先用数据标准差做归一化再整定这样mu对不同的对象都能保持在一个可预测的量级内。4.5 如何公平地比较不同数据驱动策略最后再说一个方法论上的建议。对比不同数据驱动控制策略时最容易犯的错误是“用同一个参数配置去跑所有算法”。事实上不同算法对预测时域、控制时域、权重参数的敏感度完全不同直接统一配置看起来是公平实际上反而对某些算法不公平。我的做法是把每类算法分别调参到“在验证集上表现最优”的状态然后在测试场景里进行最终对比。这样得到的比较结论反映的是方法本身的能力上限而不是某一组参数的运气。写在最后一点实际项目体会这三类数据驱动预测控制策略跑完一轮横向对比之后我最大的感受是这套技术并不是来“取代”传统模型预测控制的它更像是把建模环节的灵活性重新还给了一线工程师。对于机理清、模型成熟的对象完全没必要绕路去做数据驱动。但面对非线性强、机理不清晰、工况经常变化的系统数据驱动预测控制确实能打开新的局面。如果让我说一个最实用的组合方式我可能会建议用Hankel直接法做离线阶段的性能评估和边界探索用MFAC加预测控制做在线实时控制的主力把子空间辨识结果当作定期校准的参照基准。这种组合把三者的优势都利用上了同时把各自的短板用另外的方法补位。最后再分享一个我自己坚持了很久的经验数据驱动控制里最容易低估的是数据质量而非算法复杂度。多花时间分析历史数据覆盖范围、频谱成分和噪声水平比反复调参更有价值。方法的选择固然重要但你喂给方法的数据最终决定了控制性能的天花板。