简介：这份资源是论文《模型预测控制在实时水质调节中有多有效？状态空间建模和可扩展控制》的配套MATLAB源代码，面向给排水、市政供水与自动控制方向的研究生、科研人员及工程师，用于解决供水网络中消毒剂氯浓度的实时调节难题。代码基于EPANET-Matlab-Toolkit 2.1.8.1构建，提出水质控制的状态空间表示，并实现可扩展的模型预测控制算法，覆盖交汇处、水库、水箱、管道、泵与阀门等复杂组件。压缩包共542个文件，约12.14MB，包含130个m脚本、110个inp管网拓扑、75个dat数据、21个msx水质文件及若干shp、xlsx、pdf等，兼顾建模、仿真与结果分析。加载main.m后即可在Matlab中运行，配置集中于Constants4Concentrat相关文件。目前已有853人学习，适合希望复现论文实验、掌握状态空间建模与MPC水质控制流程的读者参考。
1. 水质预测 MATLAB 代码：从数据到控制策略的完整落地路径
做水质预测这件事，很多人第一反应是找一份能跑的 MATLAB 代码，把历史数据丢进去，看能不能拟合出未来几天的氨氮或溶解氧变化。但真正做过一轮的人会发现，预测只是中间产物，最终要回答的问题是：预测结果怎么反过来指导曝气、加药、调回流。水质建模与控制这条链路，核心难点不在模型本身，而在数据预处理、参数辨识、控制器整定这三个环节的衔接。我见过不少方案，模型 R² 能到 0.9 以上，但一接控制器就震荡，问题往往出在预测输出没有做平滑和约束。这篇内容面向的是手里有监测数据、想在 MATLAB 里跑通“预测 + 控制”闭环的从业者，不涉及任何特定项目源码，只讲可复现的通用做法。
2. 水质数据预处理与特征工程：为什么原始监测数据不能直接喂给模型
2.1 缺失值与异常值的处理逻辑
水质监测数据最常见的三个问题：传感器漂移导致的缓慢偏移、通讯中断导致的整段缺失、清洗维护导致的尖峰跳变。直接删缺失行会破坏时间连续性，直接插值会把尖峰抹平。我一般用“分段处理”策略：先标记缺失段长度，小于 3 个采样周期的用线性插值，大于 3 个周期的用同工况历史均值填充，同时生成一个缺失标记列作为模型输入。
% 假设 rawData 为 timetable，变量为 NH3N、DO、Temp、pH % 步骤1：标记缺失 missIdx = ismissing(rawData.NH3N); missLen = zeros(size(missIdx)); count = 0; for i = 1:length(missIdx) if missIdx(i) count = count + 1; else if count > 0 missLen(i-count:i-1) = count; end count = 0; end end % 步骤2：短缺失线性插值，长缺失用同月均值 filledNH3N = fillmissing(rawData.NH3N, 'linear', 'MaxGap', 3); monthIdx = month(rawData.Time); for m = 1:12 mMask = (monthIdx == m) & isnan(filledNH3N); if any(mMask) monthMean = mean(rawData.NH3N(monthIdx == m), 'omitnan'); filledNH3N(mMask) = monthMean; end end % 步骤3：异常值用 Hampel 滤波识别并替换 [cleanNH3N, outlierIdx] = hampel(filledNH3N, 5, 3);这段代码的关键参数是MaxGap和hampel的窗口与阈值。MaxGap=3意味着连续缺失不超过 3 个点才做线性插值，超过就交给月均值。hampel的窗口设为 5 是经验值，对应约 1 小时的采样频率；阈值 3 表示超过 3 倍标准差才判为异常。如果采样频率是 15 分钟一次，窗口可以缩到 3 或 4，否则会把真实突变当成异常抹掉。
2.2 特征构造：滞后项与滚动统计量
水质变量之间有明显的时滞效应，比如氨氮升高后溶解氧会在数小时后下降。只把当前时刻的 pH、温度、DO 作为输入，模型学不到这种动态。我通常构造三类特征：滞后项（前 1、2、4、8 个时刻的值）、滚动均值（窗口 4 和 12）、滚动标准差（窗口 4）。特征数量会膨胀到原来的 5 到 8 倍，但树模型和神经网络都能处理。
% 构造滞后和滚动特征 lagSteps = [1, 2, 4, 8]; rollWins = [4, 12]; features = table(); features.Time = rawData.Time; for v = {'NH3N', 'DO', 'Temp', 'pH'} vName = v{1}; for L = lagSteps features.([vName '_lag' num2str(L)]) = lag(rawData.(vName), L); end for W = rollWins features.([vName '_rollMean' num2str(W)]) = movmean(rawData.(vName), W); features.([vName '_rollStd' num2str(W)]) = movstd(rawData.(vName), W); end end % 删除因滞后产生的 NaN 行 features = rmmissing(features);这里有个容易翻车的地方：movmean和movstd默认是中心对齐，会引入未来信息。做预测任务时必须改成'Endpoints','discard'或者手动右对齐，否则离线评估指标会虚高，上线后直接崩。我一般用movmean(x, W, 'Endpoints', 'discard')然后手动补 NaN，确保每个时刻只用到过去的数据。
2.3 训练集与测试集的划分陷阱
时间序列不能随机打乱划分。常见做法是按时间顺序取前 70% 训练、后 30% 测试。但如果数据有明显的季节性，比如夏季藻类爆发导致水质突变，随机划分会让训练集和测试集分布不一致。更稳妥的方式是按完整周期划分：如果有两年数据，第一年训练，第二年测试；如果只有一年，按季度分层，每个季度取前 70% 做训练。
% 按时间顺序划分 n = height(features); splitIdx = floor(0.7 * n); trainData = features(1:splitIdx, :); testData = features(splitIdx+1:end, :); % 检查分布差异：用简单统计量对比 vars = {'NH3N_lag1', 'DO_lag1', 'Temp_lag1'}; for v = vars fprintf('%s 训练均值 %.3f 测试均值 %.3f\n', v{1}, ... mean(trainData.(v{1})), mean(testData.(v{1}))); end如果训练和测试的均值差异超过 20%，说明划分有问题，需要重新调整。这个检查步骤很多人跳过，等到模型上线后才发现预测偏差大，回头查数据已经来不及了。
3. 水质预测模型选型与 MATLAB 实现：从 ARIMA 到 LSTM 的取舍
3.1 为什么我优先试 ARIMA 而不是直接上深度学习
水质数据通常只有几千到几万个样本点，深度学习容易过拟合。ARIMA 及其变体（SARIMA、ARIMAX）在小样本上表现稳定，参数可解释，而且 MATLAB 的arima函数支持外生变量输入，能把温度、pH 作为回归项。我的一般流程是：先用 ARIMA 建立基线，如果残差还有明显自相关，再考虑 LSTM 或 GRU。
% 建立 ARIMAX 模型：NH3N 为响应，DO 和 Temp 为外生变量 y = trainData.NH3N_lag1; X = [trainData.DO_lag1, trainData.Temp_lag1]; % 自动定阶：尝试 p=1:3, d=0:1, q=1:3 bestAIC = Inf; bestModel = []; for p = 1:3 for d = 0:1 for q = 1:3 try Mdl = arima(p, d, q); EstMdl = estimate(Mdl, y, 'X', X, 'Display', 'off'); [~, ~, logL] = infer(EstMdl, y, 'X', X); [aic, ~] = aicbic(logL, p+q+1); if aic < bestAIC bestAIC = aic; bestModel = EstMdl; end catch continue; end end end end % 预测未来 24 步 futureX = [testData.DO_lag1(1:24), testData.Temp_lag1(1:24)]; [yForecast, yMSE] = forecast(bestModel, 24, 'X', futureX);这段代码的核心是arima(p,d,q)的定阶和forecast的外生变量输入。p是自回归阶数，d是差分阶数，q是移动平均阶数。水质数据一般d=0或1，p和q不超过 3。forecast的第二个参数是预测步长，24 步对应 6 小时（15 分钟采样）或 24 小时（1 小时采样）。yMSE是预测方差，可以用来构造置信区间，后续做控制时可以把置信区间作为约束。
3.2 LSTM 的 MATLAB 实现与超参数设置
如果 ARIMA 残差还有结构，或者数据量超过 2 万条，我会试 LSTM。MATLAB 的 Deep Learning Toolbox 提供了lstmLayer和trainNetwork，但直接套用默认参数容易翻车。关键超参数：隐藏单元数、初始学习率、dropout 概率、mini-batch 大小。
% 构造 LSTM 输入：序列长度 12，特征数 8 seqLen = 12; numFeatures = 8; numResponses = 1; numHiddenUnits = 64; layers = [ ... sequenceInputLayer(numFeatures) lstmLayer(numHiddenUnits, 'OutputMode', 'last') dropoutLayer(0.2) fullyConnectedLayer(numResponses) regressionLayer]; options = trainingOptions('adam', ... 'MaxEpochs', 100, ... 'MiniBatchSize', 32, ... 'InitialLearnRate', 0.005, ... 'LearnRateSchedule', 'piecewise', ... 'LearnRateDropPeriod', 30, ... 'LearnRateDropFactor', 0.5, ... 'GradientThreshold', 1, ... 'ValidationData', {XVal, YVal}, ... 'ValidationFrequency', 10, ... 'Verbose', false); % 构造训练数据：滑动窗口 XTrain = cell(size(trainData, 1) - seqLen, 1); YTrain = zeros(size(trainData, 1) - seqLen, 1); for i = 1:length(XTrain) XTrain{i} = trainData{i:i+seqLen-1, 1:numFeatures}'; YTrain(i) = trainData.NH3N_lag1(i+seqLen); end net = trainNetwork(XTrain, YTrain, layers, options);隐藏单元 64 是起点，如果验证损失震荡就降到 32，如果欠拟合就升到 128。初始学习率 0.005 配合每 30 轮减半，是我在多个水质数据集上试出来的稳定组合。GradientThreshold=1防止梯度爆炸，水质数据里偶尔有极端值，不加这个容易训练失败。dropoutLayer(0.2)是正则化，样本少于 1 万条时可以提到 0.3。
3.3 模型评估：别只看 RMSE
RMSE 对极端值敏感，水质预测里偶尔的峰值会被放大。我一般同时看三个指标：RMSE、MAE、以及峰值处的相对误差。峰值相对误差的计算方式是：取测试集里 NH3N 最高的 10% 样本，算预测值与真实值的平均绝对百分比误差。
% 评估指标 yTrue = testData.NH3N_lag1(seqLen+1:end); yPred = predict(net, XTest); rmse = sqrt(mean((yTrue - yPred).^2)); mae = mean(abs(yTrue - yPred)); % 峰值相对误差 threshold = prctile(yTrue, 90); peakIdx = yTrue >= threshold; peakMape = mean(abs((yTrue(peakIdx) - yPred(peakIdx)) ./ yTrue(peakIdx))) * 100; fprintf('RMSE: %.4f, MAE: %.4f, 峰值MAPE: %.2f%%\n', rmse, mae, peakMape);如果峰值 MAPE 超过 15%，说明模型对突变捕捉不够，需要检查特征里有没有包含足够的滞后项，或者考虑加入天气数据作为外生变量。
4. 从预测到控制：MPC 在水质调节中的参数整定与避坑
4.1 为什么用 MPC 而不是 PID
水质控制的对象有大滞后和大惯性，PID 在这种场景下要么响应慢，要么超调严重。MPC（模型预测控制）能显式处理约束，比如曝气量不能为负、加药泵有最大流量限制。MATLAB 的 Model Predictive Control Toolbox 提供了mpc对象，可以直接用预测模型作为内部模型。
% 用 ARIMAX 模型构造 MPC 内部模型 % 将 ARIMAX 转为状态空间形式 sys = ss(bestModel); % 注意：arima 对象需先转为 idpoly 再转 ss % 实际做法：用 estimate 得到的参数手动构造 idpoly % 这里假设已得到离散状态空间模型 sys % 创建 MPC 对象 Ts = 1; % 采样时间 1 小时 mpcObj = mpc(sys, Ts, 10, 3); % 预测时域 10，控制时域 3 % 设置约束 mpcObj.ManipulatedVariables(1).Min = 0; mpcObj.ManipulatedVariables(1).Max = 100; % 曝气量上限 mpcObj.ManipulatedVariables(1).RateMin = -10; mpcObj.ManipulatedVariables(1).RateMax = 10; % 设置权重 mpcObj.Weights.OutputVariables = 1; mpcObj.Weights.ManipulatedVariablesRate = 0.1; % 仿真 [ySim, uSim] = sim(mpcObj, 50, mpcObj.Model.Nominal.Y);预测时域 10 和控制时域 3 是起点。如果系统滞后大，预测时域要覆盖滞后时间，比如滞后 5 小时就设 15 以上。控制时域太小会导致控制量变化剧烈，太大会增加计算量。权重方面，输出权重 1 是基准，控制量变化率权重 0.1 是为了平滑控制动作，如果执行机构允许快速变化可以降到 0.05。
4.2 预测误差补偿：别让模型偏差直接传给执行器
MPC 依赖内部模型，但模型预测总有偏差。我一般加一个反馈校正项：用当前时刻的实际测量值与模型预测值的差，作为未来预测的补偿。
% 反馈校正 yMeas = 2.5; % 当前实测 NH3N yPredCurrent = predict(net, currentFeatures); bias = yMeas - yPredCurrent; % 将 bias 加到 MPC 的预测输出上 % 方法：修改 mpcObj.Model.Nominal.Y 或使用 output disturbance model mpcObj.Model.Nominal.Y = mpcObj.Model.Nominal.Y + bias;这个 bias 不能无限累积，我一般设一个上限，比如 ±0.5 mg/L，超过就说明模型需要重新训练，而不是继续补偿。
4.3 避坑与常见问题排查
现象一：MPC 输出震荡，执行器频繁动作。原因通常是控制时域太小或权重设置不合理。解决：把控制时域从 3 提到 5，同时把ManipulatedVariablesRate权重从 0.1 提到 0.3。
现象二：预测值一直偏高或偏低。原因可能是训练集和测试集分布不一致，或者传感器漂移。解决：检查训练集和测试集的统计量差异，如果差异大，重新划分；如果差异小，检查传感器是否需要校准。
现象三：LSTM 训练损失不下降。原因可能是学习率太大或序列长度不合适。解决：学习率降到 0.001，序列长度从 12 调到 24 或 6，观察验证损失。
现象四：ARIMA 定阶循环跑不完。原因是有组合导致estimate报错但没被捕获。解决：在try-catch里加continue，并限制最大迭代次数。
现象五：MPC 仿真报错“模型不可控”。原因通常是状态空间模型有不可控模态。解决：检查sys的可控性矩阵秩，如果秩不足，简化模型阶数或改用 ARX 模型。
5. 进阶技巧：用滚动预测和自适应权重提升长期稳定性
5.1 滚动预测：每步都用最新数据修正
离线训练好的模型上线后，输入分布会慢慢漂移。我一般用滚动预测：每来一个新样本，就把它加入输入窗口，重新预测下一步。MATLAB 里可以用predictAndUpdateState对 LSTM 做增量更新。
% LSTM 滚动预测 net = resetState(net); for t = 1:length(testData) x = testData{t, 1:numFeatures}'; [net, yPred] = predictAndUpdateState(net, x); % 存储 yPred 并与实际值对比 endresetState在每次滚动开始时调用，确保状态不跨序列累积。如果数据有长期趋势，可以每 100 步用最近的数据微调网络，但微调学习率要设得很小，比如 0.0001。
5.2 自适应权重：根据预测置信度调整控制力度
MPC 的输出权重可以动态调整。当预测方差大时，降低输出权重，让控制器更保守；预测方差小时，提高输出权重，让控制器更积极。
% 根据预测方差调整权重 forecastVar = yMSE(end); % 最后一步的预测方差 if forecastVar > 0.5 mpcObj.Weights.OutputVariables = 0.5; else mpcObj.Weights.OutputVariables = 1.5; end这个策略在雨季或水质突变期特别有用，能避免控制器因为模型不准而过度调节。
5.3 一个具体技巧：用历史相似工况做预测校正
水质系统有周期性，比如每天早晚的用水高峰。我一般维护一个历史工况库，当当前工况与某个历史片段相似时，用历史片段的实际后续值来校正模型预测。
% 相似工况检索 currentPattern = testData{end-seqLen+1:end, 1:numFeatures}; minDist = Inf; bestMatch = []; for k = 1:length(historyPatterns) dist = norm(currentPattern - historyPatterns{k}); if dist < minDist minDist = dist; bestMatch = k; end end % 用历史后续值校正 if minDist < 0.5 correction = historyOutcomes{bestMatch}(1:24) - yForecast; yForecast = yForecast + 0.3 * correction; % 0.3 是校正系数 end校正系数 0.3 是经验值，太大容易过拟合历史，太小没效果。这个技巧在数据量不足时特别管用，相当于用历史案例做数据增强。
我自己的习惯是：每次模型上线前，先用滚动预测跑一遍最近三个月的数据，看预测误差有没有随时间增大的趋势。如果有，说明模型需要重新训练或者加入新特征。这个检查花不了多少时间，但能避免很多线上翻车。希望帮到你。
本文还有配套的精品资源，点击获取
