简介：面向供水管网水质建模与控制的Matlab代码包，源自论文《模型预测控制在实时水质调节中有多有效？》，聚焦输配水网络中消毒剂浓度的实时优化问题。代码给出水质控制问题的新型状态空间表示，以及高度可扩展的模型预测控制（MPC）算法，通过增压站氯剂量与管网氯浓度间的显式关系实现仿真与调节。包内共542个文件，包括m脚本（核心算法）、inp管网模型、dat数据文件、pdf论文及msx多组分水质文件等，压缩包约12.14MB，结构紧凑。代码基于EPANET-Matlab-Toolkit运行，适合水质安全、智能供水控制方向的高年级本科生、研究生及工程师，可复现论文实验并在此基础上进行二次开发。目前已有851人学习，是理解水质建模与MPC控制结合的高质量开源参考。 水质预测这件事，我印象最深刻的是新手和非本专业的朋友几乎都会犯同一个错误：拿到监测数据就直接开跑神经网络。代码跑通的那一刻很兴奋，但放到实际场景里，预测值和真实值偏差大到没法用。做水质建模与控制，关键从来不是堆一个多复杂的模型，而是你有没有把“为什么选这个模型”“控制策略怎么跟预测结果联动”这套逻辑理顺。这个“Water-Quality-Modeling-and-Control”项目，用MATLAB把水质预测建模和反馈控制放到一个框架里做了闭环打通——既不是光做预测曲线好看，也不是只搭一套PID仿真糊弄事，而是把传感器数据、模型预测和控制器输出串成了一条完整的链路。这篇博文我会从整体设计思路、模型选择、控制策略到MATLAB代码实现逐一拆解，把我实际踩过的坑和调试经验也一并写出来，给你一条能直接落地的路径。
1. 项目整体设计与方案选型
1.1 为什么用MATLAB做水质建模和控制
水质问题本质上是一个多变量、强耦合、大迟延的复杂系统。pH、溶解氧(DO)、浊度、氨氮、总磷这些指标互相影响，而且水体本身有惯性，你投进去的药剂、曝气量要过一段时间才能看到效果。这种系统最适合用MATLAB来做，原因有三：
- 矩阵运算和向量化处理效率高，水质监测数据天生是时间序列和空间断面的矩阵结构，MATLAB处理起来比通用语言顺手得多；
- 工具箱齐全，神经网络工具箱、系统辨识工具箱、优化工具箱、Simulink控制仿真全套都有，不用在多个软件之间来回切换；
- 调试和可视化一体，
plot、scatter、heatmap几分钟就能出一张能直接拿给甲方或导师看的图。
另外一个很实在的考虑是：水质控制领域很多传统研究代码都是MATLAB写的，你拿MATLAB做，跟历史方法对标、验证模型的准确性和控制效果时，对比起来要方便得多。Python当然也可以做，但核心在控制系统仿真和实时性验证上，Simulink的物理建模能力优势就明显了。
1.2 预测与控制闭环架构设计
这个项目最关键的设计理念不是“预测完就结束”，而是把预测结果作为控制器的输入，形成闭环反馈。整个架构分成三层：
第一是数据采集层。实际的传感器监测站（pH探头、DO电极、浊度仪、氨氮分析仪）按固定的采样周期获取水质指标，这些数据在MATLAB中统一读入并转成训练集和实时输入流。
第二是预测建模层。基于历史数据训练的数据驱动模型（比如BP神经网络），或者机理模型与数据驱动结合的混合模型，对关键水质指标进行短期预测。这里要理解预测的目的：不是为了看趋势图，而是为了“提前知道未来一段时间水质会怎么变化”，给控制动作争取时间。
第三是控制执行层。控制器根据预测值与目标值之间的偏差，计算出合适的控制量。常见的比如：溶解氧偏低，就增大曝气量；pH偏离设定区间，就调整加酸或加碱泵的开度；浊度异常就加大絮凝剂投加量。控制器的输出反馈到执行机构，传感器再次采集数据，形成循环。
我实际用的方案更偏保守：先单独验证预测模型的精度，达标后再把预测模块接入Simulink仿真环境做闭环测试，最后才考虑移植到PLC或嵌入式设备。这样每一步出问题都能快速定位到是建模的锅还是控制的锅。
2. 水质预测模型的核心技术拆解
2.1 模型选型：机理模型与数据驱动模型的取舍
水质预测模型大致分两类，各有各的适用场景：
一类是基于物理化学机理的模型，比如QUAL2E、WASP这类河流水质模型，基于质量守恒、反应动力学方程来建模。这类模型的好处是物理解释性强，可以外推，但缺点是参数太多了——生化需氧量衰减系数、复氧系数、沉降速率，每一个都需要通过实验标定，标定过程又耗时间又烧钱，对于小规模项目来说根本不划算。
另一类是数据驱动模型，比如BP神经网络、LSTM、支持向量回归、随机森林等。这类模型不需要详细了解水体内部复杂的生化反应机理，只要你历史数据质量够高、样本量够大，模型就能自动学习输入输出之间的映射关系。这是目前水质预测领域绝对的主流。我在这个项目中时序性数据用得比较多，BP神经网络和LSTM都有涉及，但考虑到工程落地难度，主干模型走了BP神经网络。
选择BP神经网络的原因很务实：第一，在样本量不大的情况下，传统的BP网络比LSTM更稳，不容易出现训练不收敛的问题；第二，MATLAB的神经网络工具箱对BP的支持非常成熟，几行代码就能完成训练和仿真；第三，BP网络的推理速度极快，做实时控制时不会有计算延迟压力。LSTM适合数据很长、时序依赖特别强的场景，但调参成本和算力要求都会上一个台阶。
2.2 输入变量确定与关键参数设置
模型输入变量的选择是个关键中的关键。这里不能盲目地把所有监测指标全塞进去，一定要基于水质指标之间的相关性分析来做粗筛。操作上，我在MATLAB里用corrcoef函数计算各指标之间的相关系数矩阵，然后用heatmap可视化展示。比如溶解氧和水温通常呈负相关，pH和藻类活动强度相关，氨氮和溶解氧在一定条件下呈负相关关系。把相关性极低的指标剔除，减少输入维数，既能降低模型复杂度，又能避免引入无关噪声导致过拟合。
以溶解氧预测为例，我最终确定的输入变量包括：水温、pH、电导率、浊度、氨氮和上一时刻的溶解氧值。这里一定要包含被预测变量的历史值，因为水质过程有强惯性，当前状态很大程度上由之前的状态决定，这是时间序列预测的基本逻辑。
参数设置上我踩过的坑不少，关键参数说几个：
- 隐含层节点数：很多人直接拍脑袋选，我的经验做法是先用经验公式[n_1 = \sqrt{n + m} + a]（n为输入节点数，m为输出节点数，a取1~10之间的整数）算一个初始范围，然后从最小值开始逐次增加，观察训练集和验证集的误差变化来确定最优值。盲目设置过多节点会导致严重的过拟合，过少又欠拟合，同样一个数据集，隐含层节点数从10改成18，验证集误差可能翻倍。
- 学习率：0.01到0.1之间相对安全。学习率太大，损失函数来回震荡；太小，收敛极慢甚至卡在局部极小值。我在项目中用的是自适应学习率策略，训练初期0.05，后期衰减到0.005。
- 训练函数：我用过
trainlm（Levenberg-Marquardt），收敛速度确实快，但在数据量较大的时候内存占用较高；后来数据量大了以后改用trainscg（SCG共轭梯度法），稳定性更好。 - 归一化：必须做。MATLAB里
mapminmax函数一键搞定，把数据映射到[-1,1]区间。不做归一化的话，pH是7~8的数值，氨氮是0.5~2的数值，量纲差异会让训练过程变得极其不稳定。
2.3 模型评估：不能只看训练集误差
模型训练完，要用留出法或交叉验证来评估泛化能力。我习惯按时间顺序划分数据集：前70%做训练，中间15%做验证（用于调整超参数），最后15%做测试（评估最终模型效果）。注意水质监测数据是典型的时间序列数据，必须按时间顺序切分，不能用随机打乱的方式去划分，否则会造成数据泄漏，测试集里“偷看了”未来信息，评估结果虚高得离谱。
评价指标我用三个：
- 均方根误差（RMSE）：反映预测值与真实值的整体偏差水平，单位同原始数据，直观；
- 平均绝对百分比误差（MAPE）：反映相对误差，方便跟不同量纲的指标做对比；
- 决定系数（R²）：反映模型对真实数据变异的解释程度，越接近1越好。
实际结果中，我的溶解氧预测模型测试集R²达到了0.92左右，RMSE约0.35 mg/L，这个精度已经足以支持控制决策了。这里还想强调一下：水质预测精度不是越高越好，你要看在控制回路里对执行机构的动作影响。如果模型过于灵敏，把瞬时噪声也当成趋势预测进去，控制器就会频繁误动作，执行机构磨损加剧，投药量忽高忽低，反而把系统的稳定性破坏了。
3. 控制策略设计与Simulink仿真实现
3.1 控制目标与执行逻辑
水质控制系统跟典型的工业过程控制在逻辑上没有本质区别，核心就是感知、决策、执行三件事。控制目标需要根据实际应用场景来定，比如污水处理厂主要控制出水COD和氨氮达标，养殖水体主要控制溶解氧和pH在适宜区间。
我在项目中以溶解氧控制为一个演示对象：目标是把水体溶解氧浓度维持在6 mg/L左右。当预测模型计算出未来15分钟的溶解氧有跌破阈值的风险时，控制器就提前增加曝气设备的输出功率；反之预测偏高时，就适当降低。这个“提前量”就是预测控制的核心价值所在。
控制算法上，如果是入门级场景，PID就够用了。增量式PID的公式是：
[ \Delta u(k) = K_p [e(k) - e(k-1)] + K_i e(k) + K_d [e(k) - 2e(k-1) + e(k-2)] ]
其中e(k)是当前时刻溶解氧设定值与预测值之间的偏差，Kp、Ki、Kd分别对应比例、积分、微分系数。最终输出的是控制量的增量Δu，叠加到上一时刻的控制量上得到当前控制量，这种增量式PID的优势在于不会产生大幅的阶跃扰动，而且即使控制器出错，对执行机构的冲击也有限。
3.2 Simulink模型搭建步骤
Simulink仿真模型的搭建是项目里最能体现工程能力的地方。我搭建的模型分为四个环节：
第一步，输入信号模块。用“From Workspace”模块把实测的水质扰动信号读入模型，模拟外部环境变化，比如暴雨导致的入流负荷变化、温度骤降导致的水体复氧能力变化等。
第二步，预测模型模块。这一步有个关键技巧：直接用MATLAB Function模块嵌入训练好的BP神经网络权重矩阵，而不需要把整个神经网络工具箱都塞到Simulink里。训练好模型后，用net.IW和net.LW导出权重和偏置，然后在MATLAB Function里用纯矩阵运算实现前向传播。这样做的优势是仿真速度快，而且后续部署到嵌入式环境时几乎零移植成本。
第三步，PID控制器模块。直接用Simulink自带的PID Controller模块，调参时绑定波特图工具看相位裕度和增益裕度，比纯靠经验试凑靠谱得多。当然也可以自己搭增量式PID逻辑，两种方式我都试过，自带的模块更省事，自搭的更容易深入理解原理。
第四步，被执行对象模型。用传递函数或状态空间方程描述曝气系统对溶解氧浓度的影响。这里我用的简化一阶惯性加纯迟延模型：
[ G(s) = \frac{K}{Ts + 1} e^{-\tau s} ]
K是过程增益，T是时间常数，τ是纯迟延时间。水质过程大迟延的特性在这里体现得特别明显，τ值设置不合理会导致PID控制严重震荡。参数来源可以通过MATLAB系统辨识工具箱，对阶跃响应实验数据做拟合得到。没有实测数据的情况下，至少也要根据经验值设置合理范围，不能随意拍脑袋。
3.3 控制器参数整定经验
PID参数整定我建议分两步走。先在Simulink里用自动整定工具跑一轮，得到一个初始参数；然后在初始参数附近手动微调，观察系统的动态响应曲线。我常用的经验参数范围（针对溶解氧这类大惯性被控对象）：Kp取0.5~1.5，Ki取0.02~0.08，Kd取5~15。比例系数太大容易震荡，积分系数太大会产生明显超调，微分系数是对迟延的补偿，但太大会放大噪声。
这个环节我想重点提醒：扩散系数要跟采样周期匹配。如果你的预测模型输出频率是每15分钟一次，那PID控制周期也应该是15分钟级别，绝对不能1秒控制一次。控制周期太短，执行机构会疯狂动作，还可能导致系统不稳定。我一开始在这个问题上栽了跟头，直到看控制输出曲线才发现曝气阀门像抽筋一样来回开合，把控制周期调整到和预测频率一致后，系统立刻稳下来了。
4. 水质预测与控制的MATLAB代码实现
4.1 数据预处理与训练集构建
数据预处理的正确姿势，先用代码说明白。我读入的原始监测数据是Excel表格，每一列是一个水质指标，每一行是一个采样时刻的数据。第一步是用readtable读入数据，随后检查缺失值和异常值。缺失值我用线性插值处理，异常值则先通过3σ准则判断，剔除掉明显偏离正常范围的毛刺数据。举个例子，水温数据正常情况下在5~30℃区间，突然冒出一个50，就要警惕传感器故障或信号受干扰了。
% 读取原始监测数据 data = readtable('water_quality_data.xlsx'); % 提取关键变量 temp = data.WaterTemp; ph = data.pH; do = data.DO; % 溶解氧 turb = data.Turbidity; nh3n = data.NH3N; % 处理缺失值 - 线性插值 temp = fillmissing(temp, 'linear'); do = fillmissing(do, 'linear'); % 异常值处理 - 3σ原则 mu = mean(do); sigma = std(do); do(abs(do - mu) > 3 * sigma) = NaN; do = fillmissing(do, 'linear'); % 构建输入特征矩阵和输出向量 X = [temp, ph, turb, nh3n]; % 当前时刻及历史时刻状态 Y = do; % 预测目标 % 构造时序样本: 用前5个时刻预测下一时刻 lag = 5; X_seq = []; Y_seq = []; for i = (lag+1):length(Y) X_seq(end+1, :) = [X(i-lag:i-1, :), do(i-lag:i-1)]; Y_seq(end+1, :) = Y(i); end这里我把前5个时刻的数据组合成一个输入向量，等于把历史时序当作特征拼进了一维向量里，这是一种简单有效的处理时序问题的方式。如果你想兼顾短时和长时依赖，可以同时包含前1个时刻和前5个时刻的采样值。还有一点，如果样本量充足，可以把更多历史时刻加入输入，但是输入维度会迅速膨胀，训练时间变长。在样本量有限的情况下，不要一味增大滞后数。
4.2 BP神经网络训练核心代码
用MATLAB神经网络工具箱训练BP网络，代码极其简洁，但简洁的背后需要注意几个容易出问题的点：
% 划分训练集、验证集和测试集(按时间顺序) numSamples = size(X_seq, 1); trainNum = floor(numSamples * 0.7); valNum = floor(numSamples * 0.15); testNum = numSamples - trainNum - valNum; % 数据归一化到[-1, 1] [X_norm, ps_input] = mapminmax(X_seq', -1, 1); [Y_norm, ps_output] = mapminmax(Y_seq', -1, 1); % 创建BP神经网络 net = feedforwardnet([15, 10]); % 两层隐含层: 15个节点 + 10个节点 % 设置训练参数 net.trainFcn = 'trainlm'; net.trainParam.epochs = 1000; net.trainParam.goal = 1e-5; net.trainParam.min_grad = 1e-7; net.trainParam.max_fail = 20; % 验证集连续20次不下降就停止 % 配置数据划分 net.divideFcn = 'divideind'; net.divideParam.trainInd = 1:trainNum; net.divideParam.valInd = trainNum+1:trainNum+valNum; net.divideParam.testInd = trainNum+valNum+1:numSamples; % 训练 [net, tr] = train(net, X_norm, Y_norm); % 测试 Y_pred_norm = sim(net, X_norm(:, trainNum+valNum+1:end)); Y_pred = mapminmax('reverse', Y_pred_norm, ps_output); % 计算指标 Y_test = Y_seq(trainNum+valNum+1:end); rmse = sqrt(mean((Y_pred - Y_test).^2)); mape = mean(abs((Y_test - Y_pred) ./ Y_test)) * 100; r2 = 1 - sum((Y_test - Y_pred).^2) / sum((Y_test - mean(Y_test)).^2); fprintf('RMSE: %.4f, MAPE: %.2f%%, R2: %.4f\n', rmse, mape, r2);隐藏层节点数15和10是我经过网格搜索确定的结果，从5开始逐次加到20，对比验证集误差后选定的。另外要特别注意divideind这个划分方式：默认的dividerand是随机划分，对时间序列数据来说是灾难，必须改成按索引顺序划分，否则测试集里会出现和训练集时间重叠的样本，评估结果虚高。
还有一个小细节：max_fail的默认值通常是6，如果训练集误差持续下降但验证集误差连续6轮都不下降，训练就会提前终止。在水质数据这个场景下6往往太紧，训练早停导致模型欠拟合。我把它调整到20后，模型的稳定性明显提升了。
4.3 前向传播代码与Simulink嵌入
在Simulink中嵌入模型时，要导出神经网络的权重和偏置，用纯矩阵运算重新实现前向传播：
% 导出网络权重和偏置 W1 = net.IW{1}; % 输入层到隐含层1的权重 b1 = net.b{1}; W2 = net.LW{2,1}; % 隐含层1到隐含层2的权重 b2 = net.b{2}; W3 = net.LW{3,2}; % 隐含层2到输出层的权重 b3 = net.b{3};然后在Simulink的MATLAB Function模块里写：
function y_pred = nn_predict(x_input) % 输入x_input为归一化后的1xN特征向量 % 隐含层1: tansig激活 h1 = tansig(W1 * x_input' + b1); % 隐含层2: tansig激活 h2 = tansig(W2 * h1 + b2); % 输出层: purelin线性激活 y_pred_norm = W3 * h2 + b3; % 反归一化 y_pred = y_pred_norm * out_std + out_mean; end这里面的ps_output反归一化，我建议手动记录out_mean和out_std，直接写成代码里的常量，就不用每次都在Workspace和Simulink之间传递结构体变量了，省去不少麻烦。有个细节要注意：mapminmax的默认映射范围是[-1,1]，对应的转换公式是归一化后的值 = (原始值 - 均值) / 标准差，但mapminmax的实际计算用的是2*(x - xmin)/(xmax - xmin) - 1，所以如果你手动写反归一化，必须用mapminmax('reverse')来验证，防止公式出错。
4.4 PID控制与仿真结果解读
Simulink里PID控制器的参数设定好后，我跑了一组对比仿真。一组是纯PID控制（无预测模块），另一组是模型预测+PID控制。扰动信号我用的是一段持续上升的污染负荷，模拟污水处理厂进水的突发变化。结果显示：纯PID控制的溶解氧最低值掉到了4.2 mg/L，超过了允许下限；加入了预测模块的，溶解氧最低值保持在5.3 mg/L附近。从超调量和恢复时间看，预测+PID的方案超调量降低了约30%，恢复时间缩短了约40%。数据说明了一个道理——预测的意义不是替代反馈控制，而是给反馈控制争取时间，让控制器不至于总是“事后补救”。
5. 常见问题排查与避坑实录
5.1 模型层面的典型问题
训练不收敛或者收敛极慢，是最常见的问题。排查顺序是这样的：先看数据是否归一化，再看学习率是否过大或过小，然后检查隐含层节点数是否设置合理，最后看激活函数有没有选错。如果训练集误差能降但验证集误差飙升，那就是过拟合症状，解决方式是增加训练数据量、添加正则化项（MATLAB里performParam.regularization设置为0.01~0.1）或者减少隐含层节点数。
数据量不足在水质领域尤其常见。很多小型监测站一天一条数据甚至一周一条数据，几百条样本想训练深层网络纯属异想天开。我的对策是用数据增强手段：对原始时间序列做滑窗增采样，或者在同一量测点收集不同季节的数据合并使用。还有一个稳妥办法：如果数据量实在不够，就别碰神经网络了，老老实实改用多元线性回归或随机森林，稳定性反而更高。
5.2 控制层面的典型问题
控制效果震荡发散，优先确认被控对象的模型参数是否准确，尤其是时间常数T和纯迟延时间τ。我在调试中的一个真实经历是：把τ从5分钟改成15分钟后，原本震荡的系统一下就稳定了。原因在于控制器必须“适应”对象的迟延特性，迟延估算偏差会导致PID整定完全失真。
执行机构频繁动作的问题，前面已经提过跟控制周期有关。还有一个有效手段是给PID输出加一个“死区”：当偏差绝对值小于某个阈值（比如溶解氧偏差小于0.2 mg/L）时不改变控制输出，这样能显著减少执行机构动作频率，延长设备寿命。
5.3 工程落地时的几个提醒
从仿真走向实际应用的时候，要处理传感器噪声带来的毛刺。我的做法是加入滑动平均滤波或者一阶低通滤波，对传感器原始信号做平滑后再送入预测模型和控制模块。但注意滤波不能太重，否则信号相位滞后加剧，反而会让控制效果恶化。滤波器的时间常数建议控制在采样周期的1~2倍。
工具箱方面的兼容性问题，不同MATLAB版本之间，神经网络工具箱的函数名和参数设置有时会有差异。我用的版本里feedforwardnet是正常的，但换到旧版本可能是newff。代码迁移时先跑一遍文档自带的示例程序，能避开很多低级坑。另外，神经网络训练结果受初始权重影响极大，每次运行结果都有差异。为了解决这个问题，我在训练前用rng(42)固定了随机种子，保证可复现，这个细节在学术研究和工程交付中都很重要。
6. 一些实际的体会
这个项目做完，我最大的感受是：水质建模与控制，最终的瓶颈不在算法多先进，而在数据质量和控制对象的理解深度上。模型不需要一步到位上LSTM或强化学习那种花活，BP神经网络配合PID控制，在大多数实际场景里已经能解决90%的问题。我的建议是先将整个闭环跑通，再逐步引入更复杂的模型和算法。另一个心得是：MATLAB做这个方向确实趁手，从数据处理、模型训练到控制仿真全链路打通，基本不需要切换工具。最后分享一个我每次都会用的小技巧：在Simulink里加一个To Workspace模块把控制输出和预测值都记录下来，仿真结束后用plot对比曲线，查找问题的效率会高很多——有时候一帧图就能告诉你模型到底在哪个环节出了问题。
本文还有配套的精品资源，点击获取
