ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

嘈杂环境下数据驱动MPC实战:从DeePC求解到参数整定与避坑

嘈杂环境下数据驱动MPC实战:从DeePC求解到参数整定与避坑

简介:面向嘈杂环境下的数据驱动模型预测控制(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”时,很大一部分指的是 DeePC(Data-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[k+1] = 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]), size=T) # 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(solver=cp.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_thresh=8.0): y_med = median_filter(y_raw, size=7) 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_matrices=False) S[rank:] = 0.0 return (U * S) @ Vt

rank 怎么选?一个实用办法是画出奇异值曲线,找到明显“膝盖”的位置。以第 3 章的二阶系统为例,前两个奇异值通常比后面的高一个数量级,取 rank=2 是合理的;如果系统有积分环节或明显延迟,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, axis=0)), axis=0) viol_rate = np.sum(viol_count) / len(viol_count)

判断这一技术方向值不值得投入,我的经验是:先花两天把开环数据生成与清洗跑通,再用现有 PID 或常规 MPC 做基线,与数据驱动 MPC 在同样的扰动序列下对比。如果数据驱动版本在约束条件相同的情况下,RMS 跟踪误差能下降 20% 以上,并且控制量动作量没有翻倍,就值得继续投入精力和预算;否则再天花乱坠的“无模型”故事也不值得为它付出额外的维护成本。黑匣子轻便,但必须建立在可验证的性能增益上,否则就会变成后期排错的地狱入口。我自己的教训是,曾经花了两周做一版端到端学习控制,最后还是在闭环里加了 PID 垫底才敢交给现场,白白多付出一倍时间。希望这篇笔记里的整定顺序和参数准则能帮你省掉那一周多的返工,也希望你下次从压缩包里解出来的数据驱动 MPC,不再是玄学,而是你手里握得住的一套测量、整定与验证流程。希望帮到你。

本文还有配套的精品资源,点击获取

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进