简介这份资源是论文《模型预测控制在实时水质调节中有多有效状态空间建模和可扩展控制》的配套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:iseqLen-1, :); % 转置为 特征×时间 YTrain{i} y(iseqLen:iseqLenpredLen-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.00150 轮后降一半。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(k1) 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):(twindowSize-1); X_window X(idx, :); % 预测下一步 y_pred predict(net, {X_window}); predictions(t) y_pred{1}(1); % 实际值模拟中从数据取实际中从传感器读 actuals(t) y(twindowSize); % 可选在线更新每 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可以做集成但时序数据我更倾向手动加权。最后说一个我自己的习惯每次改完代码先跑一个最小可复现示例确认数据流和维度都对再上全量数据。水质数据维度多、缺失多维度对不上是最高频的翻车点。希望帮到你。本文还有配套的精品资源点击获取