简介：这份资源是论文《模型预测控制在实时水质调节中有多有效？状态空间建模和可扩展控制》的配套MATLAB源代码，面向给排水、市政供水网络与自动控制方向的研究生、科研人员及工程师，用于解决供水网络中消毒剂氯浓度实时调节的建模与控制难题。压缩包共542个文件，约12.14MB，包含130个m脚本、110个inp管网拓扑、75个dat数据、30个pdf文献、21个msx与17个net模型文件，以及shp、xlsx、xml等地理与配置数据，覆盖建模、仿真与结果分析全流程。资源基于EPANET-Matlab-Toolkit 2.1.8.1开发，核心函数集中在WQMC文件夹，加载main.m后即可运行，配置参数位于Constants4Concentrat相关文件中。目前已有853人学习下载。读者可获取状态空间建模的完整实现、可扩展MPC算法代码及多场景管网算例，便于复现论文实验、对比控制策略并迁移到自身供水网络的水质调控研究中。
1. 水质预测 MATLAB 代码：从数据到控制策略的完整落地路径
做水质预测这件事，很多人第一反应是找现成的工具箱或者调包，但真正落到工程现场，你会发现核心难点根本不在模型本身，而在于数据怎么组织、参数怎么标定、控制策略怎么和预测结果闭环。水质建模与控制这个方向，本质上是用数学模型描述水体中污染物迁移转化规律，再基于预测结果设计控制方案——比如调节曝气量、投加药剂量、调整进水流量。它适合环境工程、自动化、市政水务方向的研究生和工程师，尤其是手头有监测数据但不知道怎么变成可运行 MATLAB 代码的人。
我见过太多人卡在第一步：数据格式不对、时间戳不统一、缺失值处理随意，后面模型跑出来全是玄学。这篇内容就按我实际做项目的顺序，把水质预测 MATLAB 代码从数据预处理、模型选型、参数辨识到控制策略实现的完整链路拆开讲，每一步都给可复现的代码和参数说明，最后说几个我踩过的坑。
2. 数据预处理与特征工程：把原始监测数据变成模型能吃的矩阵
2.1 水质数据的典型结构与读取方式
实际项目里拿到的水质数据通常来自在线监测站或实验室采样，格式无非几种：CSV、Excel、数据库导出。字段一般包括时间戳、pH、溶解氧（DO）、化学需氧量（COD）、氨氮、总磷、浊度、水温。常见做法是先用readtable读进来，统一时间列格式，再按时间排序。
% 读取水质监测数据，假设文件为 CSV 格式 opts = detectImportOptions('water_quality_data.csv'); opts.VariableNamingRule = 'preserve'; % 保留原始列名 data = readtable('water_quality_data.csv', opts); % 统一时间列格式，假设第一列是时间 data.Time = datetime(data.Time, 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); data = sortrows(data, 'Time'); % 按时间升序排列 % 查看前几行和变量类型 head(data); summary(data);这段代码的关键在于detectImportOptions自动识别列类型，但实际数据里经常有混合类型导致识别错误，所以我会手动检查summary输出。时间列必须转成datetime类型，后面做重采样和插值都依赖它。如果时间格式不统一，比如有的带秒有的不带，先用datetime的InputFormat参数强制解析，解析失败的用NaT标记，后续统一处理。
2.2 缺失值处理与异常值剔除的实操策略
水质数据缺失是常态，传感器故障、通讯中断、维护窗口都会导致缺值。我一般不用均值填充，因为水质参数有强时序相关性，均值填充会破坏趋势。常用做法是线性插值加前后向填充组合，异常值用滑动窗口的 MAD（中位数绝对偏差）检测。
% 假设 data 中需要处理的变量为 DO、COD、NH3N vars = {'DO', 'COD', 'NH3N'}; for i = 1:length(vars) col = data.(vars{i}); % 线性插值填补缺失值 col = fillmissing(col, 'linear', 'EndValues', 'nearest'); % 滑动窗口 MAD 异常检测，窗口大小 24（假设小时数据） win = 24; med = movmedian(col, win, 'omitnan'); mad = movmad(col, win, 'omitnan'); threshold = 3; % 3 倍 MAD 作为阈值 outliers = abs(col - med) > threshold * mad; % 将异常值替换为 NaN 后再插值 col(outliers) = NaN; col = fillmissing(col, 'linear', 'EndValues', 'nearest'); data.(vars{i}) = col; end % 检查处理后是否还有缺失 sum(ismissing(data(:, vars)))这里fillmissing的'linear'方法适合连续缺失不超过 3 个点的场景，如果缺失段很长，线性插值会失真，这时候我会考虑用同期历史数据回归填充或者直接标记为无效样本。MAD 阈值取 3 是经验值，水质数据波动大的时候可以放宽到 4，但太宽会漏掉真实异常。movmedian和movmad的窗口大小要根据采样频率定，小时数据用 24，分钟数据用 60 或 120。
2.3 特征构造：滞后项、滑动统计量与归一化
水质预测模型输入不能只有当前时刻值，污染物浓度有滞后效应，所以要把过去几个时刻的值作为特征。同时滑动均值和滑动标准差能捕捉局部波动趋势。归一化用 min-max 或 z-score 都行，但要注意训练集和测试集要用同一套参数。
% 构造滞后特征和滑动统计量 lagSteps = [1, 2, 3, 6, 12]; % 滞后步数 vars = {'DO', 'COD', 'NH3N', 'pH', 'Temp'}; for i = 1:length(vars) for lag = lagSteps data.([vars{i} '_lag' num2str(lag)]) = lagmatrix(data.(vars{i}), lag); end % 滑动均值与标准差，窗口 6 和 24 data.([vars{i} '_ma6']) = movmean(data.(vars{i}), 6, 'omitnan'); data.([vars{i} '_std6']) = movstd(data.(vars{i}), 6, 'omitnan'); data.([vars{i} '_ma24']) = movmean(data.(vars{i}), 24, 'omitnan'); end % 删除因滞后产生的 NaN 行 data = rmmissing(data); % 归一化，保存参数供测试集使用 featureNames = data.Properties.VariableNames(3:end); % 假设前两列是时间和其他非特征列 mu = mean(data{:, featureNames}); sigma = std(data{:, featureNames}); dataNorm = (data{:, featureNames} - mu) ./ sigma; % 转成矩阵供模型使用 X = dataNorm;lagmatrix生成的滞后特征会导致前几行 NaN，用rmmissing删掉。归一化参数必须保存，测试集用同样的mu和sigma，否则数据泄露。特征数量会膨胀得很快，5 个变量 × 5 个滞后 + 3 个滑动统计 = 40 个特征，后面做特征选择或者降维是必要的。
3. 水质预测模型选型与 MATLAB 实现：从机理模型到数据驱动
3.1 机理模型与数据驱动模型的适用边界
水质建模有两条路：机理模型和数据驱动模型。机理模型比如 ASM 系列（活性污泥模型）、QUAL 系列，优点是物理意义明确、可解释性强，缺点是参数多、标定困难、对数据要求高。数据驱动模型比如 LSTM、SVR、随机森林，优点是上手快、预测精度往往更高，缺点是黑匣子、外推能力差。
我的经验是：如果要做控制策略且对可解释性有要求，用简化机理模型加参数辨识；如果只是做短期预测且数据量充足，用 LSTM 或集成学习。MATLAB 里 Deep Learning Toolbox 和 Statistics and Machine Learning Toolbox 都能直接用。
3.2 用 LSTM 做多步水质预测的完整代码
LSTM 适合处理时序数据，MATLAB 里构建训练流程比较直接。下面是一个预测未来 6 步 COD 的示例。
% 准备 LSTM 训练数据，假设 X 是归一化后的特征矩阵，y 是目标 COD % 构造序列数据，输入序列长度 48，输出序列长度 6 seqLen = 48; predLen = 6; numSamples = size(X, 1) - seqLen - predLen + 1; XTrain = cell(numSamples, 1); YTrain = cell(numSamples, 1); for i = 1:numSamples XTrain{i} = X(i:i+seqLen-1, :)'; % 转置为 特征×时间 YTrain{i} = y(i+seqLen:i+seqLen+predLen-1)'; end % 划分训练集和验证集 idx = randperm(numSamples); trainRatio = 0.8; trainIdx = idx(1:round(trainRatio*numSamples)); valIdx = idx(round(trainRatio*numSamples)+1:end); XTrainSet = XTrain(trainIdx); YTrainSet = YTrain(trainIdx); XValSet = XTrain(valIdx); YValSet = YTrain(valIdx); % 定义 LSTM 网络结构 numFeatures = size(X, 2); numResponses = predLen; layers = [ sequenceInputLayer(numFeatures) lstmLayer(128, 'OutputMode', 'sequence') dropoutLayer(0.2) lstmLayer(64, 'OutputMode', 'last') fullyConnectedLayer(numResponses) regressionLayer ]; % 训练选项 options = trainingOptions('adam', ... 'MaxEpochs', 100, ... 'MiniBatchSize', 32, ... 'InitialLearnRate', 0.001, ... 'LearnRateSchedule', 'piecewise', ... 'LearnRateDropPeriod', 50, ... 'LearnRateDropFactor', 0.5, ... 'GradientThreshold', 1, ... 'ValidationData', {XValSet, YValSet}, ... 'ValidationFrequency', 10, ... 'Verbose', false, ... 'Plots', 'training-progress'); % 训练网络 net = trainNetwork(XTrainSet, YTrainSet, layers, options); % 预测 YPred = predict(net, XValSet);网络结构里第一层 LSTM 输出序列，第二层只输出最后时刻，这样能捕捉长期依赖后做多步预测。dropoutLayer防过拟合，水质数据噪声大时很有用。学习率用分段衰减，初始 0.001，50 轮后降一半。GradientThreshold设 1 防止梯度爆炸。验证集比例 20%，如果数据量少可以调到 30%。
3.3 模型评估与超参数调优的实操要点
评估指标用 RMSE、MAE、MAPE，但水质预测里 MAPE 在浓度接近零时会爆炸，所以我主要看 RMSE 和 NSE（纳什效率系数）。超参数调优用贝叶斯优化，MATLAB 的bayesopt可以直接接。
% 定义超参数优化变量 optimVars = [ optimizableVariable('numHidden1', [64, 256], 'Type', 'integer') optimizableVariable('numHidden2', [32, 128], 'Type', 'integer') optimizableVariable('dropout', [0.1, 0.5]) optimizableVariable('learnRate', [1e-4, 1e-2], 'Transform', 'log') ]; % 定义目标函数 objFcn = @(params) trainLSTMForBayes(XTrainSet, YTrainSet, XValSet, YValSet, params); % 贝叶斯优化 results = bayesopt(objFcn, optimVars, ... 'MaxObjectiveEvaluations', 30, ... 'IsObjectiveDeterministic', false, ... 'UseParallel', false); % 获取最优参数 bestParams = bestPoint(results);trainLSTMForBayes函数内部用当前参数构建网络、训练、返回验证集 RMSE。贝叶斯优化 30 次评估在单机上大概跑几小时，如果时间紧可以减到 15 次。UseParallel开并行需要 Parallel Computing Toolbox，没有就设 false。
4. 控制策略设计与 MATLAB 仿真：把预测结果变成执行动作
4.1 基于预测的曝气控制逻辑
水质控制里最典型的是曝气控制，目标是把溶解氧维持在设定范围。预测模型给出未来 6 步 DO 趋势，控制器根据趋势提前调整曝气量。我用的是预测控制（MPC）思路，MATLAB 里有 Model Predictive Control Toolbox，但也可以手写简化版。
% 简化 MPC 控制器，控制曝气量 % 状态：DO 浓度，控制量：曝气量 u % 预测模型：DO(k+1) = a*DO(k) + b*u(k) + c*COD(k) a = 0.85; b = 0.12; c = -0.03; % 模型参数，需辨识 DO_setpoint = 2.0; % 设定值 mg/L u_min = 0; u_max = 1; % 曝气量范围 N = 6; % 预测步长 % 当前时刻状态 DO_current = 1.5; COD_current = 30; % 滚动优化 u_opt = zeros(N, 1); DO_pred = zeros(N, 1); DO_pred(1) = a*DO_current + b*u_opt(1) + c*COD_current; for k = 2:N DO_pred(k) = a*DO_pred(k-1) + b*u_opt(k) + c*COD_current; end % 简单比例控制作为对比 Kp = 0.5; u_p = Kp * (DO_setpoint - DO_current); u_p = max(u_min, min(u_max, u_p)); fprintf('MPC 预测 DO 轨迹：%.2f %.2f %.2f %.2f %.2f %.2f\n', DO_pred); fprintf('比例控制输出：%.2f\n', u_p);这段代码只是演示逻辑，实际 MPC 需要在线求解优化问题，用quadprog或 MPC Toolbox。模型参数a、b、c要通过历史数据辨识，可以用最小二乘或nlarx。预测步长 6 对应未来 6 个采样周期，太长计算量大，太短控制效果差。
4.2 参数辨识：从数据里反推模型系数
机理模型参数辨识用lsqcurvefit或nlinfit，数据驱动模型参数就是训练出来的权重。这里说一个简化的一阶模型辨识。
% 假设有输入输出数据 u_data 和 y_data % 模型：y(k) = a*y(k-1) + b*u(k-1) % 构造回归矩阵 N = length(y_data); Y = y_data(2:end); Phi = [y_data(1:end-1), u_data(1:end-1)]; % 最小二乘辨识 theta = Phi \ Y; a_hat = theta(1); b_hat = theta(2); % 验证 y_pred = zeros(N, 1); y_pred(1) = y_data(1); for k = 2:N y_pred(k) = a_hat*y_pred(k-1) + b_hat*u_data(k-1); end rmse = sqrt(mean((y_data - y_pred).^2)); fprintf('辨识参数 a=%.4f, b=%.4f, RMSE=%.4f\n', a_hat, b_hat, rmse);Phi \ Y是 MATLAB 的左除，等价于最小二乘解。辨识前要对数据去均值或差分，否则会有直流偏置导致参数有偏。验证时用递推预测而不是一步预测，更能反映模型长期性能。
4.3 仿真验证：闭环控制效果怎么评估
闭环仿真要把预测模型和控制器接起来跑，看设定值跟踪和抗扰动能力。
% 闭环仿真 T_sim = 200; DO = zeros(T_sim, 1); u = zeros(T_sim, 1); DO(1) = 1.0; DO_setpoint = 2.0; for k = 2:T_sim % 控制器 error = DO_setpoint - DO(k-1); u(k) = max(0, min(1, 0.5*error + 0.1*sum(error))); % 被控对象 DO(k) = 0.85*DO(k-1) + 0.12*u(k) + 0.02*randn(); end % 计算性能指标 ISE = sum((DO_setpoint - DO).^2); IAE = sum(abs(DO_setpoint - DO)); fprintf('ISE=%.2f, IAE=%.2f\n', ISE, IAE);randn()模拟过程噪声，实际项目里噪声统计特性要从数据里估计。ISE 和 IAE 越小控制效果越好，但也要看控制量是否频繁波动，执行机构寿命也要考虑。
5. 避坑与排查：水质预测 MATLAB 代码里最容易翻车的五个地方
5.1 数据泄露：归一化参数用了全量数据
现象：模型在训练集和测试集上表现都很好，但上线后预测完全不准。原因：归一化时用了全部数据的均值和标准差，测试集信息泄露到训练过程。解决：划分训练集和测试集后再计算归一化参数，测试集用训练集的参数变换。
5.2 滞后特征构造错误：lagmatrix 方向搞反
现象：模型预测值总是滞后于真实值，或者超前。原因：lagmatrix默认正滞后是向过去取，但有人误用负值导致取到未来值。解决：明确lagmatrix(x, 1)表示 x(t-1)，构造完检查前几行 NaN 位置是否正确。
5.3 LSTM 训练不收敛：学习率太大或序列太长
现象：训练损失震荡不下降，或者直接 NaN。原因：初始学习率 0.01 以上容易震荡，序列长度超过 100 梯度容易爆炸。解决：学习率从 0.001 开始，加GradientThreshold，序列长度控制在 50 以内。
5.4 控制量饱和：执行机构物理限制没考虑
现象：仿真里控制效果很好，实际设备要么全开要么全关。原因：控制器输出没有限幅，或者限幅后没有抗积分饱和。解决：u = max(u_min, min(u_max, u))，积分项加限幅或条件积分。
5.5 模型外推失效：训练数据覆盖范围太窄
现象：遇到极端天气或突发污染，预测完全离谱。原因：数据驱动模型只在训练分布内有效。解决：监测输入特征是否超出训练范围，超出时切换到保守控制策略或机理模型兜底。
6. 进阶技巧：用滚动预测和在线更新让模型持续可用
水质过程是慢时变的，传感器会漂移，季节变化会改变生化反应速率。一个训练好就不管的模型，几个月后精度就会下降。我一般用滚动预测加在线更新来维持性能。
滚动预测的做法是每次只预测未来一步，然后把实际观测值加入输入窗口，再预测下一步。这样误差不会累积，但计算频率高。在线更新是每隔一段时间用新数据微调模型，LSTM 可以用trainNetwork继续训练，但要注意灾难性遗忘，学习率要设得很小。
% 滚动预测示例 windowSize = 48; predHorizon = 1; numSteps = 100; predictions = zeros(numSteps, 1); actuals = zeros(numSteps, 1); for t = 1:numSteps % 取当前窗口数据 idx = (t):(t+windowSize-1); X_window = X(idx, :)'; % 预测下一步 y_pred = predict(net, {X_window}); predictions(t) = y_pred{1}(1); % 实际值（模拟中从数据取，实际中从传感器读） actuals(t) = y(t+windowSize); % 可选：在线更新，每 20 步微调一次 if mod(t, 20) == 0 % 用最近数据微调，学习率设小 options = trainingOptions('adam', 'InitialLearnRate', 1e-5, ... 'MaxEpochs', 5, 'Verbose', false); net = trainNetwork({X_window}, {y_pred}, net.Layers, options); end end % 计算滚动预测 RMSE rmse_rolling = sqrt(mean((predictions - actuals).^2)); fprintf('滚动预测 RMSE=%.4f\n', rmse_rolling);在线更新时MaxEpochs设 5 以内，学习率 1e-5 量级，否则会把之前学到的模式覆盖掉。如果数据量足够，更好的做法是维护一个滑动窗口训练集，每次用最近几个月的数据重新训练，而不是在旧模型上微调。
另一个技巧是集成预测，用 LSTM、SVR、随机森林各跑一遍，取加权平均。权重按验证集误差倒数分配。这样单模型失效时整体不会崩。MATLAB 里fitrensemble可以做集成，但时序数据我更倾向手动加权。
最后说一个我自己的习惯：每次改完代码，先跑一个最小可复现示例，确认数据流和维度都对，再上全量数据。水质数据维度多、缺失多，维度对不上是最高频的翻车点。希望帮到你。
本文还有配套的精品资源，点击获取
