简介：这份开源代码对应论文《模型预测控制在实时水质调节中有多有效？状态空间建模和可扩展控制》，面向给排水、市政水务与环境工程方向的研究生、科研人员及工程师，用于解决供水网络中消毒剂水平实时调节这一难题。资源针对复杂组件（交汇处、水库、水箱、管道、泵和阀门）缺乏面向控制建模与可扩展控制算法的问题，给出水质控制问题的状态空间表示，并设计高度可扩展的模型预测控制算法，展示快速响应与对不确定性的适应能力。压缩包共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 代码：从数据到控制策略的完整落地路径
做环境工程或过程控制的朋友大概率遇到过这种场景：手里有一组监测站的历史水质数据，领导要你三天内给出未来一周的氨氮、溶解氧变化趋势，还要附带一套可执行的曝气控制方案。用 Excel 拉趋势线显然交不了差，Python 写 LSTM 又来不及调参。这时候一套结构清晰、开箱能跑的 MATLAB 水质建模与控制代码就成了救命稻草。Water-Quality-Modeling-and-Control 这个资源包正是冲着这个需求去的——它把水质预测模型和后续控制策略串成了一条线，不是单纯给你一个预测脚本，而是从数据预处理、模型训练、参数辨识一路铺到控制器设计。适合环境工程专业的研究生、做污水处理过程控制的工程师，以及需要快速搭建水质仿真平台验证算法的人。系统开源的好处在于，你能看到每一行状态方程怎么写的，而不是面对一个黑匣子干瞪眼。
2. 拆开资源包：水质模型选型与 MATLAB 实现骨架
2.1 为什么是机理模型加数据驱动，而不是纯神经网络
拿到这份代码之前，先想清楚一个选型问题：水质预测到底该用机理模型还是数据驱动模型。纯数据驱动的方法，比如 LSTM 或 GRU，在数据量充足时预测精度确实好看，但一旦进水水质突变或者传感器漂移，模型输出可能完全偏离物理常识——溶解氧预测出负值这种事，血泪经验告诉我们并不罕见。而纯机理模型，比如经典的 ASM1 或 BSM1，参数标定工作量巨大，且对水厂实际工况的适配性差。
这份资源包走的是中间路线：用机理模型搭建水质变化的动力学骨架，再用实测数据去校正关键参数。代码里能看到状态方程以微分方程组形式写在单独的 .m 文件中，而参数辨识部分则用最小二乘或粒子群算法去拟合。这种混合策略的好处是，预测结果既有物理约束兜底，又能贴合实际数据。对于做控制的人来说，机理模型还能直接拿来做控制器设计，不用再费劲去线性化一个神经网络。
2.2 代码目录结构与核心模块拆解
资源包解压后，目录结构大致如下（不同版本可能略有差异，以实际为准）：
Water-Quality-Modeling-and-Control/ ├── data/ # 存放实测水质数据 │ ├── raw/ # 原始监测数据 │ └── processed/ # 预处理后的数据 ├── models/ # 模型定义 │ ├── water_quality_model.m # 水质机理模型 │ └── parameter_id.m # 参数辨识脚本 ├── control/ # 控制策略 │ ├── pid_controller.m # PID 控制器 │ └── mpc_controller.m # 模型预测控制 ├── utils/ # 工具函数 │ ├── data_preprocess.m # 数据清洗与归一化 │ └── plot_results.m # 结果可视化 ├── main_simulation.m # 主仿真入口 └── README.md # 简要说明核心模块有三个：水质机理模型、参数辨识、控制器。机理模型部分通常包含几个关键状态变量——溶解氧、氨氮、硝态氮、异养菌和自养菌生物量。状态方程的形式类似：
function dXdt = water_quality_model(t, X, u, p) % 水质机理模型状态方程 % X: 状态变量 [SO, SNH, SNO, XB_H, XB_A] % u: 输入变量 [Q_in, SO_in, SNH_in] % p: 模型参数结构体 % 溶解氧平衡 dSO = (u(1)*(u(2) - X(1))/V) - ... (p.mu_H * X(4) * X(1))/(p.K_O + X(1)) / p.Y_H - ... (p.mu_A * X(5) * X(1))/(p.K_O_A + X(1)) / p.Y_A + ... p.KLa * (p.SO_sat - X(1)); % 氨氮平衡 dSNH = (u(1)*(u(3) - X(2))/V) - ... p.mu_A * X(5) * X(2)/(p.K_NH + X(2)) / p.Y_A - ... p.i_NXB * p.mu_H * X(4) * X(1)/(p.K_O + X(1)); % 硝态氮平衡 dSNO = -u(1)*X(3)/V + ... p.mu_A * X(5) * X(2)/(p.K_NH + X(2)) / p.Y_A; % 异养菌和自养菌生物量 dXB_H = p.mu_H * X(4) * X(1)/(p.K_O + X(1)) - p.b_H * X(4); dXB_A = p.mu_A * X(5) * X(2)/(p.K_NH + X(2)) - p.b_A * X(5); dXdt = [dSO; dSNH; dSNO; dXB_H; dXB_A]; end这段代码的逻辑是：每个状态变量的变化率由进水稀释、生物反应消耗/生成、以及曝气或沉淀等物理过程共同决定。参数结构体 p 里包含最大比增长速率 mu_H、mu_A，半饱和常数 K_O、K_NH，衰减系数 b_H、b_A 等。这些参数不是拍脑袋定的，需要用实测数据去辨识。代码里参数辨识部分通常用lsqnonlin或fmincon做优化，目标函数是模型输出与实测值的残差平方和。
提示：如果你拿到的数据时间步长和模型积分步长不一致，先用
interp1做重采样，否则 ode45 会报错或者结果失真。
2.3 数据预处理与参数辨识的实操步骤
数据预处理这一步，很多人会跳过或者随便填个缺失值，结果后面模型怎么调都不对。我一般会强制走一遍流程：先剔除传感器明显漂移的异常段，再用线性插值补缺失，最后做归一化。代码里data_preprocess.m大致长这样：
function [data_clean, data_norm] = data_preprocess(raw_data, method) % 数据预处理：异常值剔除、缺失值插补、归一化 % raw_data: 原始数据矩阵，每列一个变量 % method: 'zscore' 或 'minmax' % 1. 用 3-sigma 准则剔除异常值 mu = mean(raw_data, 'omitnan'); sigma = std(raw_data, 'omitnan'); outlier_idx = abs(raw_data - mu) > 3 * sigma; raw_data(outlier_idx) = NaN; % 2. 线性插值补缺失 for i = 1:size(raw_data, 2) nan_idx = isnan(raw_data(:, i)); if any(nan_idx) valid_idx = find(~nan_idx); raw_data(nan_idx, i) = interp1(valid_idx, raw_data(valid_idx, i), ... find(nan_idx), 'linear', 'extrap'); end end data_clean = raw_data; % 3. 归一化 if strcmp(method, 'zscore') data_norm = (data_clean - mean(data_clean)) ./ std(data_clean); elseif strcmp(method, 'minmax') data_norm = (data_clean - min(data_clean)) ./ (max(data_clean) - min(data_clean)); end end参数说明：raw_data的列顺序要和模型状态变量对应，不然辨识出来的参数张冠李戴。method选 zscore 适合数据分布接近正态的情况，minmax 适合有明确上下界的变量比如溶解氧浓度。参数辨识阶段，把预处理后的数据喂给parameter_id.m，它会调用优化器迭代调整参数，直到模型输出和实测值的误差降到可接受范围。常见做法是先用文献值做初值，再限定参数上下界，避免优化器跑到物理上不合理的区域。
3. 从预测到控制：PID 与 MPC 控制器的 MATLAB 落地
3.1 控制器设计前必须搞清的被控变量与执行器
预测模型跑通之后，下一步就是拿它做控制。但很多人在这里翻车——直接把预测模型的输出当成被控量，却忘了执行器是什么。在水质控制场景里，最常见的执行器是曝气风机和回流泵。被控变量通常是溶解氧浓度和出水氨氮浓度。曝气量增大，溶解氧上升，硝化反应加快，氨氮下降，但能耗也上去了。这是一个典型的多人多出耦合系统。
代码里提供了两种控制器：PID 和 MPC。PID 结构简单，适合单回路控制，比如只控溶解氧。MPC 则能处理多变量约束，比如在保证出水氨氮达标的前提下最小化曝气能耗。选哪个取决于你的执行器数量和算力预算。如果只是做课程设计或者快速验证，PID 足够；如果要做能耗优化或者应对进水波动，MPC 更合适。
3.2 PID 控制器的参数整定与抗积分饱和
PID 控制器在pid_controller.m里实现，核心代码不复杂，但有几个坑要注意：
function [u, e_int] = pid_controller(setpoint, y, e_int, Kp, Ki, Kd, dt, u_min, u_max) % 带抗积分饱和的 PID 控制器 % setpoint: 设定值 % y: 当前测量值 % e_int: 积分项累积值 % Kp, Ki, Kd: PID 参数 % dt: 采样周期 % u_min, u_max: 执行器上下限 e = setpoint - y; e_int = e_int + e * dt; % 计算未限幅输出 u_raw = Kp * e + Ki * e_int + Kd * (e - (setpoint - y)) / dt; % 抗积分饱和：如果输出饱和，停止积分累积 if u_raw > u_max u = u_max; e_int = e_int - e * dt; % 回退积分 elseif u_raw < u_min u = u_min; e_int = e_int - e * dt; else u = u_raw; end end逻辑说明：抗积分饱和是必须加的，否则设定值突变时积分项会累积到爆炸，等误差反向时控制器半天退不出来，溶解氧会像过山车一样震荡。参数整定方面，我一般先用 Ziegler-Nichols 法粗调，再根据仿真曲线微调。Kp 太大会导致曝气量频繁切换，风机寿命受影响；Ki 太大会让溶解氧响应变慢，氨氮可能超标。
注意：PID 的采样周期要和传感器更新周期匹配。如果传感器每 5 分钟才给一个数，你设 1 秒的 dt 就是自欺欺人。
3.3 MPC 控制器的预测时域与约束设置
MPC 控制器在mpc_controller.m里，代码量比 PID 大不少，但核心思想就三步：预测、优化、滚动。预测时域 Np 和控制时域 Nc 是两个关键参数。Np 太短，控制器看不到未来约束，容易短视；Np 太长，计算量上去，实时性受影响。对于水质控制这种慢过程，Np 取 10 到 20 个采样周期比较常见，Nc 取 3 到 5。
约束设置是 MPC 的精髓。比如曝气量有物理上限，出水氨氮有排放标准上限，这些都要写成约束条件。代码里用quadprog求解二次规划问题，目标函数通常是跟踪误差加权和加控制增量加权和。权重矩阵 Q 和 R 需要调，Q 大则跟踪精度优先，R 大则控制动作平滑优先。没有万能参数，只能根据你的工艺特点去试。
4. 避坑与排查：水质建模控制中那些让人头疼的问题
4.1 模型输出发散或出现负浓度
现象：仿真跑着跑着溶解氧变成负的，或者氨氮浓度指数上升。原因通常是参数辨识时没有加物理约束，导致反应速率常数取到了不合理的大值，或者积分步长太大，ode45 的误差控制没起作用。解决方法是给参数加上下界，比如 mu_H 限制在 0.1 到 6 之间，同时在 ode45 的 options 里把 RelTol 和 AbsTol 调小，比如 1e-6。如果还不行，换用 ode15s 处理刚性方程。
4.2 参数辨识不收敛或拟合结果无物理意义
现象：lsqnonlin跑完提示达到最大迭代次数，或者拟合出的半饱和常数比文献值大两个数量级。原因可能是初值选得太离谱，或者目标函数存在多个局部极小值。解决方法是先用文献典型值做初值，再用多起点优化，比如跑 10 次不同初值的lsqnonlin，取残差最小的那组。另外检查数据预处理有没有把异常值清干净，一个坏点能把整个辨识带偏。
4.3 控制器震荡或执行器频繁动作
现象：溶解氧设定值 2 mg/L，实际值在 1.5 到 2.5 之间来回摆，曝气阀门一直咔咔响。原因通常是 PID 参数太激进，或者 MPC 的权重矩阵 R 设得太小。解决方法是先降 Kp 和 Ki，观察响应变慢但震荡减小，再逐步加回来。MPC 的话把 R 调大，让控制增量代价更高。另外检查传感器有没有滤波，原始信号噪声大也会导致控制器误动作。
4.4 仿真速度慢到无法接受
现象：跑一次 7 天的仿真要十几分钟，调参调到崩溃。原因可能是模型里用了变步长求解器但事件检测太频繁，或者 MPC 每步都在线求解 QP 问题。解决方法是把机理模型的积分步长固定，用定步长求解器；MPC 可以改成显式 MPC 或者查表法，牺牲一点精度换速度。如果只是做参数扫描，用parfor并行跑。
4.5 数据时间戳错位导致预测滞后
现象：模型预测的峰值总是比实测晚几个小时。原因通常是数据对齐没做好，进水流量和出水水质的时间戳差了一个水力停留时间。解决方法是根据反应池容积和流量算出理论停留时间，对输入数据做相应的时间平移。这个坑很隐蔽，因为模型本身没问题，就是数据对错了。
5. 进阶技巧：用状态观测器补全不可测变量与仿真验证
实际水厂里，异养菌和自养菌的生物量浓度是没法在线测量的，但模型里又需要这些状态变量。这时候就得用状态观测器，比如扩展卡尔曼滤波或者龙伯格观测器，用可测的溶解氧和氨氮去估计不可测的生物量。代码包里虽然没有直接给观测器实现，但你可以基于现有的机理模型自己扩展。我一般会先线性化模型，在几个典型工况点求雅可比矩阵，然后设计观测器增益。验证方法是把估计值和仿真真值对比，看收敛速度和稳态误差。
另一个进阶用法是把这套代码当成测试平台，验证你自己的想法。比如你想试试模糊控制或者强化学习，不用从头搭模型，直接替换control/目录下的控制器文件就行。主仿真入口main_simulation.m负责调度，你只要保证输入输出接口一致。我习惯在改控制器之前，先把机理模型在开环状态下跑一遍，确认状态变量都在合理范围，再闭环。从那以后我每次换控制器都强制走一遍开环验证，省得浪费时间在错误的基础上调参。希望帮到你。
本文还有配套的精品资源，点击获取
