简介面向嘈杂环境下的数据驱动模型预测控制MPC研究这份Matlab程序包提供了完整的可运行实验环境。程序包以Matlab 2014、2019a、2024a为运行平台内置附赠案例数据采用参数化编程方式用户可快速修改参数以适配不同控制场景适合控制工程领域的学生、工程师及研究人员用来验证算法或开展二次开发。整个压缩包共包含2002个文件总大小约13.85MB其中以1697个Matlab脚本.m为主体另有269个HTML辅助文档、C语言源文件、PDF说明等覆盖模型构建、控制器求解与结果分析等环节代码注释详细、结构清晰降低了学习门槛。目前该资源已有57人学习感兴趣者可借此直观感受数据驱动MPC在噪声干扰下的鲁棒性与自适应能力为复杂系统控制方案的设计与优化提供参考。1. 嘈杂环境中的数据驱动 MPC数据越“脏”越不能上来就调包现场工程师拿到一个标着“嘈杂环境中的数据驱动 MPC.rar”的压缩包时心情通常是矛盾的传统模型预测控制最折磨人的环节就是机理建模而数据驱动 MPC 号称能从历史数据里直接把预测模型学出来省掉一次次阶跃测试和辨识实验但“嘈杂环境”四个字又提醒你这批代码要面对的是传感器噪声、工艺扰动和执行器死区叠加在一起的真实工业数据。结论先说数据驱动 MPC 在嘈杂环境下完全能用但直接拿开环数据训练再闭环保真大概率前几步就翻车。这篇笔记把从数据生成、Hankel 矩阵构造、DeePC 求解到参数整定和闭环验证的完整路径拆开讲适合锅炉燃烧优化、压缩机控制、机器人轨迹跟踪这些被机理模型逼疯、又怕噪声毁掉闭环的从业者。2. 数据驱动 MPC 的三条技术路线子空间、DeePC 与伪数据驱动陷阱先讲原理。数据驱动 MPC 并不是丢掉优化框架而是把“预测模型”这个环节从机理白盒换成数据黑匣子。标准 MPC 在每个控制周期要做的事是“基于当前状态预测未来 N 步输出在输入输出约束下求使代价函数最小的控制序列”传统做法需要被控对象的传递函数或状态空间模型数据驱动 MPC 的差异只在于预测这一动作不再依赖物理方程而是由历史输入输出数据决定的线性映射或参数化表示。2.1 路线一子空间辨识 标准 MPC控制器内核不动只换模型来源最保守、也最适合工程师上手的数据驱动 MPC 路线是子空间辨识。它的思想是直接从一批足够激励的输入输出数据里估计出一个低阶状态空间模型A, B, C, D然后用这个模型继续跑你已经写好的标准 MPC 求解器。常见做法是先构造由过去和未来的输入输出构成的块 Hankel 矩阵做一次 LQ 分解或 SVD把可观测矩阵分离出来进而得到状态序列估计最后用最小二乘解出系统矩阵。这条路的优点在于你已有的 MPC 代码、约束和权重全都不用改等于只把“模型来源”换掉了。在噪声环境下子空间辨识最容易踩的坑是将来噪声当作状态通道导致控制量频繁抖动。常规解决办法是使用工具变量法用过去的输入作为工具变量来消除噪声项与状态的关联本质是在 SVD 之前先做一个投影去偏。做完这步后续 MPC 的鲁棒性基本等同你原有控制器比较可控。2.2 路线二DeePC用数据本身当预测模型不是先学模型再用模型近五年学术界和工程界讨论“数据驱动 MPC”时很大一部分指的是 DeePCData-enabled Predictive Control数据驱动预测控制。它依据 Willems 基本引理对一个可控的线性系统只要输入序列持续激励那么任意一条足够长的历史输入输出轨迹都可以表示为同一个 Hankel 矩阵行空间的线性组合。也就是说数据块本身就可以当预测模型用不需要先提炼出 A, B, C, D。DeePC 每一拍求解的是一个带约束的二次规划设 T_ini 步历史窗口、N 步预测窗口把输入和输出拼成一个块 Hankel 矩阵再引入一个系数向量 g让“Hankel 矩阵乘 g”同时复现过去 T_ini 步和预测未来 N 步然后对预测输出施加入力约束和代价函数。噪声环境下这个等式约束不能精确满足因此要加松弛变量并由正则化项惩罚。比起子空间辨识DeePC 的优势是省去中间模型辨识环节劣势是矩阵规模随数据量和窗口长度膨胀实时性压力明显更大。2.3 路线三伪数据驱动的端到端学习控制第三种常见做法是端到端策略学习即不构造 Hankel 矩阵也不设计预测时域直接把当前状态送到神经网络或高斯过程回归里输出控制量。它在演示视频里效果惊人但在嘈杂的工业环境里我一般不会优先推荐原因有三个一是它没有预测时域的概念也没有约束的可行性保证本质上已经不是 MPC二是噪声会被当作状态特征直接映射到控制量噪声方向的小扰动都可能带来执行机构的大动作三是这类方法的数据需求量和闭环验证成本远高于前两者。与其说它是数据驱动 MPC不如说是“伪数据驱动”的强化学习应用。2.4 选型对照数据量、噪声强度与实时性路线需要的数据量噪声鲁棒性约束处理实时性子空间辨识 标准 MPC中等半小时开环激励足够依赖辨识去偏手段中高靠原有 MPC 层保障高模型阶次低DeePC大Hankel 矩阵需要足够列依赖正则化与松弛变量中等直接在优化中施加中等矩阵规模大端到端学习很大且需要闭环试错低噪声易被放大难以严格保证推理虽快训练与验证成本高对绝大多数带噪声的过程控制场景我会优先建议子空间辨识打底、DeePC 做对比验证如果追求论文级效果或系统强烈非线性再考虑更复杂的变体。选型时还有一个容易被忽略的现实因素你的团队后续要维护的是哪个代码栈。子空间辨识加标准 MPC 的维护成本最低DeePC 一旦需要调窗长和正则化排错难度会陡增。3. 最小复现跑通一个带测量噪声的 mpc 模型预测控制闭环标题里的“rar”通常意味着拿到的是一批脚本、数据和说明文档而不是一个可以 pip install 的库。这类工程包常见的组织方式是一个生成开环数据的脚本、一个构造 Hankel 矩阵并求解优化问题的脚本、一个闭环仿真主循环以及一个集中存放参数的配置文件。建议先把这个骨架建起来再往里填你的真实数据比直接改网上陌生人贴的单个文件要稳得多。3.1 先梳理包结构把运行入口和参数入口分开一个我常用的结构大致是这样data_driven_mpc/ ├── data_gen.py # 开环激励 噪声注入生成 u_history.npy / y_history.npy ├── deepc_core.py # 构造块 Hankel 矩阵定义并求解 DeePC 优化问题 ├── run_closed_loop.py # 闭环主循环调用 deepc_core 里的求解器 └── params.yaml # 所有窗口长度、正则化权重、约束边界集中管理这样组织有两个好处换数据时不用改求解器代码调参数时不破坏代码结构。无论你拿到的压缩包内部怎么命名建议先按这个思路把文件和参数重新归位否则后续更换采样周期或数据文件时到处找超参数的体验会非常痛苦。这看起来像流程问题实际是很多工程师最初的“玄学”烦恼来源参数明明改了代码里却被另一处硬编码覆盖掉了。3.2 第一步生成带噪声的开环数据集假设被控对象是一个二阶系统状态方程为x[k1] A x[k] b u[k] w[k]测量输出为y[k] c x[k] v[k]其中w是过程噪声、v是测量噪声。用随机二进制序列做输入激励模拟工厂里常见的 PRBS 测试代码保留随机种子以保证结果可复现import numpy as np # 被控对象二阶离散状态空间采样周期 0.1s A np.array([[0.95, 0.10], [-0.05, 0.85]]) b np.array([[0.10], [0.20]]) c np.array([[1.0, 0.5]]) T 2000 # 开环数据长度 rng np.random.default_rng(42) u rng.choice(np.array([-1.0, 1.0]), sizeT) # PRBS 式激励 u u 0.01 * rng.standard_normal(T) # 输入端小幅扰动 x np.zeros((2, T)) y np.zeros((T, 1)) w_std 0.05 # 过程噪声标准差 v_std 0.10 # 测量噪声标准差 for k in range(T - 1): x[:, [k 1]] A x[:, [k]] b * u[k] w_std * rng.standard_normal((2, 1)) y[k] c x[:, [k]] v_std * rng.standard_normal() u u.reshape(-1, 1) np.save(u_history.npy, u) np.save(y_history.npy, y) print(SNR(dB):, 20 * np.log10(np.std(y) / v_std))这段代码的要点是激励信号要足够“有效”。随机二进制序列相比正弦扫频的优势是能同时激发多个频率段让 Hankel 矩阵行满秩的概率更高。测量噪声标准差 v_std 决定了信噪比工程上一般先调到能看到输出趋势、又不至于淹没信号的程度建议设成输出标准差的 5% 到 15%。过程噪声 w_std 不要设成零否则闭环仿真的“嘈杂环境”就名不副实了。3.3 第二步实现 DeePC 核心求解器DeePC 的求解逻辑是把“历史轨迹必须被复现”和“未来轨迹要尽量接近参考”同时塞进一个带约束的二次规划。假设 T_ini 为过去窗口长度N 为预测时域块 Hankel 矩阵由输入轨迹和输出轨迹上下堆叠而成总行数为 2倍的T_ini N。先构建矩阵def build_hankel(u, y, T_ini, N, T_data): L T_ini N # 每条轨迹的长度 n_cols T_data - L 1 H_u np.zeros((L * u.shape[1], n_cols)) H_y np.zeros((L * y.shape[1], n_cols)) for j in range(n_cols): H_u[:, j] u[j:j L].flatten() H_y[:, j] y[j:j L].flatten() return np.vstack([H_u, H_y]) # 上半块是输入轨迹下半块是输出轨迹接下来是每一拍的优化求解。把 Hankel 矩阵乘上系数向量 g就得到一整条长度为 T_ini N 的输入轨迹和输出轨迹。上半块前 T_ini 行应等于过去输入下半块前 T_ini 行应等于过去输出允许松弛上半块后 N 行是未来的控制序列下半块后 N 行是预测输出import cvxpy as cp import numpy as np def solve_deepc(H, T_ini, N, u_past, y_past, U_ub, U_lb, lambda_g, lambda_y): L T_ini N n_g H.shape[1] g cp.Variable(n_g) sigma_y cp.Variable(N) # 输出等式松弛变量吸收噪声影响 y_ref np.zeros(N) # 参考轨迹实际场景由设定值平移后得到 u_traj H[0:L, :] g y_traj H[L:2 * L, :] g u_past_pred u_traj[0:T_ini] y_past_pred y_traj[0:T_ini] u_fut_pred u_traj[T_ini:L] y_fut_pred y_traj[T_ini:L] cost cp.sum_squares(y_fut_pred - y_ref) \ lambda_g * cp.sum_squares(g) \ lambda_y * cp.sum_squares(sigma_y) constraints [ u_past_pred u_past.flatten(), y_past_pred y_past.flatten() sigma_y, U_lb u_fut_pred, u_fut_pred U_ub, ] prob cp.Problem(cp.Minimize(cost), constraints) prob.solve(solvercp.CLARABELLE) # 也可换成 OSQP / ECOS return g.value, sigma_y.value这里 cost 的三项分别对应跟踪误差、系数向量范数、输出松弛变量范数。lambda_g 是抵抗噪声过拟合的关键lambda_y 是允许“输出历史轨迹不必被严格复现”的退路。代码中行索引的规则是H 的上半块是输入轨迹下半块是输出轨迹每一块内部前 T_ini 行是过去、后 N 行是未来。真实验证时建议先在开环数据上做一步预测测试确认索引对齐了再接入闭环否则会出现看似正常、其实完全失效的“静默翻车”。3.4 第三步闭环主循环把求解器嵌进被控对象闭环仿真的每一拍要做三件事取当前窗口的输入输出作为初始条件、调用求解器得到系数 g、把 g 对应的第一个未来输入施加给被控对象然后更新状态T_sim 800 T_ini 50 N 30 H build_hankel(u, y, T_ini, N, T) u_sim np.zeros((T_sim, 1)) y_sim np.zeros((T_sim, 1)) # 前 T_ini 步用零输入建立初始窗口实际工程中来自上一段运行历史 for k in range(T_sim): if k T_ini: u_sim[k] 0.0 # 更新被控对象状态y_sim[k] 为带噪声测量 continue u_past u_sim[k - T_ini:k] y_past y_sim[k - T_ini:k] g_vec, _ solve_deepc(H, T_ini, N, u_past, y_past, 2.0, -2.0, 1e-3, 1e2) u_sim[k] H[T_ini, :] g_vec # 只取未来控制序列的第一个值 # 在这里把 u_sim[k] 施加给被控对象更新状态并记录 y_sim[k]关键不在循环本身而在三个细节。第一每次“只取第一个控制量”就是滚动优化的意义如果为了省事把全部序列都执行一遍就不是 MPC 了。第二T_ini 长度要大于系统记忆尤其对含积分过程的系统太短会导致初始条件信息不足。第三如果输入端有执行器饱和约束不能只加在优化变量上仿真环境里的被控对象也要加同样的限幅否则控制器会以为它给了一个 2.0 的输入实际执行机构只动到 1.5闭环表现自然对不上。4. 嘈杂环境的建模链路数据清洗、正则化与参数整定这一篇集中回答“数据脏了怎么办”。数据驱动 MPC 的上限由数据决定不是由求解器决定的。哪怕优化器写得再漂亮喂进去的是一段带趋势漂移、异常脉冲和时变增益的数据学出来的预测也一样飘。所以进入参数整定之前先把数据处理做扎实。4.1 数据清洗异常值剔除、趋势项消除与均值去偏工业数据的第一次清洗是异常值检测。传感器卡死、通信抖动会产生单点跳变这种脉冲在 Hankel 矩阵里会污染一整列。常见做法是先做中位数滤波再做阈值判别如果某一时刻的测量偏离局部中位数超过 k 倍标准差就把该点标记为异常并按线性插值替换。第二步是去趋势因为开环激励期间工艺基线可能缓慢漂移这会让数据不满足时不变假设。工程上我会先做一阶差分并检查差分序列是否平稳再决定是减趋势项还是高通滤波。最后一定要减均值因为 DeePC 的初始条件里包含直流分量如果不把均值单独建模控制器会为一个实际上不存在的直流工作点浪费执行机构行程。import numpy as np from scipy.ndimage import median_filter def clean_data(y_raw, k_thresh8.0): y_med median_filter(y_raw, size7) resid y_raw - y_med std_r np.std(resid) mask np.abs(resid) k_thresh * std_r y_clean y_raw.copy() y_clean[mask] y_med[mask] y_clean y_clean - np.mean(y_clean) # 去均值直流工作点另行建模 return y_clean, mask参数 k_thresh 一般取 6 到 10。取 6 会更激进适合传感器经常毛刺的现场取 10 更保守适合本身较干净但偶尔有周期扰动的数据。中位数滤波器窗口 size 取 5 到 9窗口太小去不掉成串的尖峰窗口太大会把真实阶跃变化也抹平。清洗完成后建议打印被标记的样本占总样本的比例如果超过 2%就要怀疑是传感器硬件问题而不是单纯的偶发噪声。4.2 正则化与秩约束对抗噪声的后悔药噪声环境下DeePC 的 Hankel 矩阵几乎必然满秩因为噪声给矩阵注入了无关的奇异值。此时不加正则优化器会利用这些虚假自由度把历史轨迹“完美复现”但预测出的未来轨迹对噪声极度敏感。两个后悔药一个是给系数 g 加 2 范数正则逼着解不要依赖过大的系数组合另一个更本质是直接对 Hankel 矩阵做截断 SVD预估系统阶次 n只保留前 n 个奇异值把剩下的当作噪声子空间丢掉。实际操作中我一般两个都用先截断把有效秩降下来再保留一个小量 lambda_g 防止截断边界处数值抖动。截断 SVD 的 Python 写法很简短def truncate_hankel(H, rank): U, S, Vt np.linalg.svd(H, full_matricesFalse) S[rank:] 0.0 return (U * S) Vtrank 怎么选一个实用办法是画出奇异值曲线找到明显“膝盖”的位置。以第 3 章的二阶系统为例前两个奇异值通常比后面的高一个数量级取 rank2 是合理的如果系统有积分环节或明显延迟rank 适当增加到 4 到 6。注意截断之后 H 不再严格等于原始数据矩阵但求解器并不关心这一点它只关心行空间是否被合理压缩。4.3 参数整定第一轮可直接抄的准则初学者拿到窗长 T_ini、预测时域 N、数据长度 T_data、lambda_g 和 lambda_y通常无从下手。下面是一组在多个仿真工程上验证得比较顺的初始值参数第一轮取值调整方向T_data2000 点以上噪声大就加到 5000保证 Hankel 列数远大于有效秩T_ini对象阶次估值的 3 到 5 倍系统记忆长就加大过大会让可用列数变少N覆盖上升时间的 1.5 倍太短会短视激进太长松弛变量权力过大lambda_g1e-3 到 1e-2噪声大则调大先抑制振荡再回调lambda_y1e1 到 1e3噪声大应允许更大松弛但过大会放弃历史信息一个实用的整定顺序是先用截断 SVD 确定有效阶次接着固定 T_ini 和 N然后单独观察给定值阶跃下的闭环响应。如果输出在设定值附近高频抖动优先加大 lambda_g如果输出迟缓甚至不跟踪再检查 lambda_y 是否太大。噪声大的系统里lambda_y 偏大比偏小的后果更具隐蔽性——每一步求解都合理但长时间看就是“躺平不干活”控制器把一切偏差都归咎于松弛变量而 g 几乎没有作用。5. 避坑记录数据驱动 MPC 最容易翻车的 5 个场景下面这些坑不是理论里的假想是实际操作中最常被反复问到的五类典型故障。每条按现象、原因、解决的顺序写对应到代码和参数上可以直接排查。5.1 现象训练误差很小闭环仿真前几步就发散这是最让人沮丧的一种场景离线重构历史轨迹时误差只有 0.01一接闭环第一个时刻输出就飞了。原因通常是 Hankel 矩阵的行空间包含了噪声对应的高频方向离线预测用的是同一批数据的内插所以“预测”天然容易过拟合闭环时初始条件是从在线测量取的噪声方向与训练数据不一致问题立刻暴露。解决对 H 做截断 SVD、降低有效秩并把 lambda_g 从 1e-3 往上调到 1e-1 左右观察同时确认初始条件里没有把噪声放大测量滤波如果已经显著滞后也不要用过大的 T_ini。5.2 现象开环激励数据明明很丰富求解时却报“不可行”DeePC 里最容易被忽略的约束是输入输出的初始轨迹等式约束。如果 T_ini 步的过去数据在 Hankel 矩阵里没有对应的列能匹配等式约束就会冲突。原因通常是激励不充分激励信号在某些幅值区间停留时间太短导致 Hankel 矩阵中某些行的线性组合覆盖不到当前工作点。解决检查持续激励条件把输入改为幅值在 [-2, 2] 之间随机切换、每段保持 3 到 5 个采样周期的 PRBS如果仍然不可行可以先把初始轨迹等式约束加松弛变量但这种方式只能用来排查不能作为长期方案留在控制器里。5.3 现象输出跟踪看着不错控制量却像锯齿一样疯狂切换控制量高频抖动的直接原因是代价函数里缺少对控制量变化量的惩罚或者有效秩没有被截断。很多工程师第一反应是调权重但真正的问题是模型把噪声当成了可控动态被控对象本质上对高频控制不那么敏感而噪声让预测输出对高频输入非常敏感。解决顺序是先做截断 SVD 降低有效秩再在代价函数里增加对u_fut_pred[1:] - u_fut_pred[:-1]的惩罚项。如果你现在的代码没有构造对控制增量的惩罚即便把 lambda_g 调大也只能缓解系数爆炸不能直接抑制控制量切换频率。5.4 现象把训练数据换到另一个工况闭环性能明显下降数据驱动 MPC 的前提是工作点附近的线性化。嘈杂环境还叠加了另一个问题同一个过程在不同负荷下增益和动态特性变化很大一组数据训练出来的 Hankel 矩阵不可能覆盖所有工况。解决用分段建模每个工况段单独构造 Hankel 矩阵与训练数据在控制器里根据当前工作点切换数据块或者在数据生成阶段加入双脉冲测试以覆盖不同幅值区间的增益再用离线闭环仿真验证。如果你的现场工况切换很频繁那么纯线性的数据驱动 MPC 也许不该是首选子空间辨识加增益调度会更务实。5.5 现象按论文默认参数跑通实测性能比不过 PID很多初次尝试者忽略了一个事实论文里的 lambda_g 和 lambda_y 是针对某个特定系统调出来的换系统就必须重新整定PID 的整定经验也不能直接照搬。这里面有个深层原因DeePC 的代价函数各项量纲不一致输出误差的平方和与 g 的平方和可能相差好几个数量级论文里设成 1e-3 的 lambda_g 放到实际系统里相当于没有约束。解决先把所有信号做归一化再整定然后从一次开环预测实验的误差尺度出发让 lambda_g 与跟踪误差项的数值量级处于同一范围。不要盲从压缩包里的默认配置那些配置通常只对生成它的那一台设备有效。6. 验证闭环性能蒙特卡罗扰动测试与三个可复现指标数据驱动控制的汇报容易让人“看体验”但工程决策需要数字。单次仿真具有偶然性噪声环境下我一般至少跑 30 种不同的随机种子每次使用不同的噪声序列然后统计三个指标RMS 跟踪误差、控制量变化总量、约束违反率。计算代码很简短err_rms np.sqrt(np.mean(err_seq ** 2)) act_var np.mean(np.abs(np.diff(u_seq, axis0)), axis0) viol_rate np.sum(viol_count) / len(viol_count)判断这一技术方向值不值得投入我的经验是先花两天把开环数据生成与清洗跑通再用现有 PID 或常规 MPC 做基线与数据驱动 MPC 在同样的扰动序列下对比。如果数据驱动版本在约束条件相同的情况下RMS 跟踪误差能下降 20% 以上并且控制量动作量没有翻倍就值得继续投入精力和预算否则再天花乱坠的“无模型”故事也不值得为它付出额外的维护成本。黑匣子轻便但必须建立在可验证的性能增益上否则就会变成后期排错的地狱入口。我自己的教训是曾经花了两周做一版端到端学习控制最后还是在闭环里加了 PID 垫底才敢交给现场白白多付出一倍时间。希望这篇笔记里的整定顺序和参数准则能帮你省掉那一周多的返工也希望你下次从压缩包里解出来的数据驱动 MPC不再是玄学而是你手里握得住的一套测量、整定与验证流程。希望帮到你。本文还有配套的精品资源点击获取