简介：这份资源是论文《模型预测控制在实时水质调节中有多有效？状态空间建模和可扩展控制》的配套MATLAB源代码，面向给排水、市政供水与自动控制方向的研究生、科研人员及工程师，用于解决供水网络中消毒剂水平实时调节这一难题。代码提出水质控制问题的状态空间表示，明确氯投加量与管网氯浓度之间的关系，并设计可扩展的模型预测控制算法，兼顾响应速度与不确定性适应能力。压缩包共542个文件，约12.14MB，包含130个m脚本、110个inp管网拓扑、75个dat数据、21个msx水质模型及若干shp、xlsx、pdf等，覆盖建模、仿真与结果分析全流程。资源依赖EPANET-Matlab-Toolkit 2.1.8.1，加载main.m后即可运行，配置集中于Constants4Concentrat相关文件。目前已有853人学习，适合希望复现论文实验、掌握管网水质建模与MPC控制实现细节的读者参考。
1. 水质预测 MATLAB 代码：从数据到控制策略的完整落地路径
做环境工程或过程控制的朋友大概率遇到过这种场景：手里有一堆监测站的历史数据，领导要你下周给出未来三天的水质预判，还要顺带说清楚该怎么调曝气量。用 Excel 硬拟合吧，精度不够；上深度学习吧，数据量又撑不起来。这套 Water-Quality-Modeling-and-Control 的 MATLAB 代码包，恰好卡在中间那个位置——它用机理模型加系统辨识的思路，把水质预测和后续控制策略串成了一条线。代码结构清晰，主脚本加函数模块的经典布局，适合有 MATLAB 基础、想快速搭出可运行原型的人。你不需要从头推公式，但得知道每个模块在算什么，不然调参时就是盲人摸象。
2. 水质模型选型与数据预处理：为什么不能直接上神经网络
2.1 机理模型和黑箱模型的边界在哪
水质系统有个特点：变量之间不是纯统计关系，而是有物质守恒约束的。比如溶解氧的消耗和有机物降解速率挂钩，总磷的变化受底泥释放影响。如果你直接拿 LSTM 去拟合，数据量少的时候它会把噪声当规律学，换一个季节就崩。这套代码走的是灰箱路线——用一阶动力学方程搭骨架，再用最小二乘辨识关键系数。常见做法是先把 BOD、COD、氨氮、溶解氧这几个核心指标的变化写成微分方程，然后离散化求解。
代码里model_definition.m这个文件定义了状态变量和速率方程。我一般会先打开它看三件事：状态变量有几个、每个方程里哪些参数是待辨识的、有没有饱和项。饱和项很关键，比如溶解氧低于某个阈值时，硝化反应速率会骤降，这个非线性特征如果被忽略，预测曲线在低氧区间会偏得离谱。
% model_definition.m 片段：状态变量与速率方程 % 状态向量 x = [BOD; NH4; DO; TP] % 参数向量 p = [k_d; k_n; k_r; theta] function dxdt = water_quality_ode(t, x, p, u) % x(1): BOD浓度, x(2): 氨氮浓度, x(3): 溶解氧, x(4): 总磷 % p(1): BOD降解系数, p(2): 硝化速率常数 % p(3): 复氧系数, p(4): 温度修正系数 % u: 外部输入，如曝气量、进水流量 T = 20; % 参考温度，实际使用时从数据里读 theta = p(4); k_d = p(1) * theta^(T - 20); k_n = p(2) * theta^(T - 20); k_r = p(3); % 溶解氧饱和浓度，用温度经验公式 DO_sat = 14.6 - 0.4 * T + 0.008 * T^2; % 速率方程 r_d = k_d * x(1); % BOD降解耗氧 r_n = k_n * x(2) * x(3) / (x(3) + 2); % 硝化反应，带半饱和项 r_r = k_r * (DO_sat - x(3)); % 复氧 dxdt = zeros(4,1); dxdt(1) = -r_d + u(1); % BOD: 降解 + 进水 dxdt(2) = -r_n + u(2); % 氨氮: 硝化消耗 + 进水 dxdt(3) = -r_d - 4.57*r_n + r_r + u(3); % 溶解氧: 消耗 + 复氧 + 曝气 dxdt(4) = -0.1*r_d + u(4); % 总磷: 与BOD耦合释放 + 进水 end这段代码的逻辑是：先算温度修正后的速率常数，再算三个核心反应速率，最后按质量平衡写微分方程。参数theta一般取 1.02 到 1.08 之间，如果你拿到的数据跨越了冬夏两季，这个值不调准，预测曲线在温度突变时会明显滞后。u向量是外部输入，实际用的时候需要把监测站的进水流量和曝气量按时间对齐后插值到求解步长上。
2.2 数据清洗：缺失值和异常值的处理顺序不能反
原始监测数据里，缺失值和异常值往往同时存在。我见过有人先填缺失再剔异常，结果把插值出来的假数据当真实值保留下来。正确顺序是：先识别异常并置为 NaN，再统一插值。异常识别用滑动窗口的 3σ 准则，窗口宽度取 24 小时比较稳，太短会误杀真实波动，太长会漏掉突发排放。
% data_preprocess.m 片段：异常检测与插值 function data_clean = preprocess_raw(data_raw, window_hours) % data_raw: 时间戳 + 各指标浓度 % window_hours: 滑动窗口小时数 data_clean = data_raw; vars = {'BOD','NH4','DO','TP'}; for i = 1:length(vars) col = data_raw.(vars{i}); % 滑动窗口3σ检测 win = window_hours; mu = movmean(col, win, 'omitnan'); sigma = movstd(col, win, 'omitnan'); outlier = abs(col - mu) > 3 * sigma; col(outlier) = NaN; % 线性插值，限制最大连续缺失长度 col = fillmissing(col, 'linear', 'MaxGap', 6); data_clean.(vars{i}) = col; end end参数MaxGap设成 6 的意思是：连续缺失超过 6 个采样点就不插值了，直接留 NaN，后面建模时跳过这段。这个阈值根据你的采样频率来定，如果是小时数据，6 小时以上的缺失用线性插值就是自欺欺人。滑动窗口的movstd计算量不小，数据超过十万条时建议先降采样再检测，或者改用中位数绝对偏差，对异常值更鲁棒。
3. 参数辨识与预测验证：最小二乘的坑比你想的多
3.1 用 lsqnonlin 做参数辨识的实操配置
参数辨识的本质是找一组 p 让模型输出和实测值的残差平方和最小。MATLAB 里lsqnonlin是常用选择，但默认配置在 Wasserstein 距离意义下可能收敛到局部极小。我一般会做三件事：设合理的上下界、给初值、开并行计算。
% parameter_estimation.m 片段 function p_est = estimate_params(data_clean, p0, lb, ub) % 提取输入输出 t_span = data_clean.time; y_meas = [data_clean.BOD, data_clean.NH4, data_clean.DO, data_clean.TP]; u_input = [data_clean.inflow, data_clean.ammonia_load, ... data_clean.aeration, data_clean.phosphorus_load]; % 目标函数：仿真输出与实测的残差 obj_fun = @(p) simulate_and_residual(p, t_span, y_meas, u_input); % 配置优化选项 options = optimoptions('lsqnonlin', ... 'Display', 'iter', ... 'Algorithm', 'trust-region-reflective', ... 'MaxFunctionEvaluations', 5000, ... 'UseParallel', true); % 多初值尝试，避免局部极小 p_est = []; best_resnorm = inf; for k = 1:5 p_try = p0 .* (1 + 0.2 * randn(size(p0))); % 初值扰动 [p_temp, resnorm] = lsqnonlin(obj_fun, p_try, lb, ub, options); if resnorm < best_resnorm best_resnorm = resnorm; p_est = p_temp; end end endAlgorithm选trust-region-reflective是因为它能处理边界约束，比levenberg-marquardt更适合有物理意义的参数范围。MaxFunctionEvaluations设 5000 是经验值，水质模型一般 2000 到 3000 次就能收敛，设太大浪费时间。多初值扰动那步很关键，我遇到过 BOD 降解系数在 0.1 和 0.5 两个初值下收敛到完全不同的结果，残差差了三倍。扰动幅度 0.2 是试出来的，太大容易跑到不合理区域，太小起不到探索作用。
3.2 预测验证：别只看 RMSE
很多人验证模型就报一个 RMSE，这在水质预测里不够。因为水质指标的量纲和波动范围差异大，BOD 可能几百，溶解氧只有个位数，同一个 RMSE 对两者的意义完全不同。我习惯用三个指标：NSE、PBIAS 和 RSR。NSE 看整体拟合，PBIAS 看系统性偏差，RSR 是 RMSE 和标准差的比值，能反映模型是否比直接用均值预测更好。
% validate_model.m 片段 function metrics = calc_metrics(y_true, y_pred) % Nash-Sutcliffe效率系数 NSE = 1 - sum((y_true - y_pred).^2) / sum((y_true - mean(y_true)).^2); % 百分比偏差 PBIAS = 100 * sum(y_true - y_pred) / sum(y_true); % RMSE-标准差比值 RSR = sqrt(mean((y_true - y_pred).^2)) / std(y_true); metrics = struct('NSE', NSE, 'PBIAS', PBIAS, 'RSR', RSR); end判断标准：NSE 大于 0.65 算可接受，大于 0.8 算好；PBIAS 绝对值小于 15% 说明没有明显系统偏差；RSR 小于 0.7 表示模型优于均值预测。如果 NSE 高但 PBIAS 也大，说明模型在某个区间系统性偏高或偏低，这时候要回去检查速率方程里有没有漏掉某个源汇项。
4. 控制策略实现：从预测到执行的最后一公里
4.1 模型预测控制框架的简化实现
预测做完了，接下来是怎么用。这套代码里带了一个简化版的模型预测控制（MPC）框架，核心思路是：用辨识好的模型往前预测 N 步，在控制步长内优化曝气量等操作变量，让溶解氧维持在目标区间，同时 minimize 能耗。
% mpc_controller.m 片段 function u_opt = mpc_step(x_current, p_est, target_DO, horizon) % x_current: 当前状态 % target_DO: 溶解氧目标值 % horizon: 预测时域步数 n_u = 1; % 控制变量：曝气量 u_min = 0; u_max = 100; % 曝气量范围 % 目标函数：跟踪误差 + 控制增量惩罚 cost_fun = @(u_seq) mpc_cost(u_seq, x_current, p_est, target_DO, horizon); % 用fmincon求解带约束优化 options = optimoptions('fmincon', 'Display', 'off', ... 'Algorithm', 'sqp', 'MaxIterations', 100); u0 = 50 * ones(horizon, 1); % 初值 u_opt = fmincon(cost_fun, u0, [], [], [], [], ... u_min*ones(horizon,1), u_max*ones(horizon,1), [], options); u_opt = u_opt(1); % 只取第一步执行 end function J = mpc_cost(u_seq, x0, p, target, N) J = 0; x = x0; for k = 1:N % 用欧拉法前进一步 dxdt = water_quality_ode(0, x, p, [0;0;u_seq(k);0]); x = x + dxdt * 0.1; % 步长0.1小时 % 跟踪误差 J = J + (x(3) - target)^2; % 控制平滑惩罚 if k > 1 J = J + 0.1 * (u_seq(k) - u_seq(k-1))^2; end end endhorizon一般取 10 到 20 步，对应未来 1 到 2 小时。步长 0.1 小时是折中，太小计算量大，太大精度不够。控制增量惩罚系数 0.1 是调出来的，设大了曝气量变化太慢跟不上水质波动，设小了阀门频繁动作影响设备寿命。fmincon用sqp算法是因为目标函数有非线性约束，虽然这里只写了边界约束，但实际工程中可能还要加溶解氧变化率限制。
4.2 控制效果评估：别忽略执行器延迟
仿真里控制效果很好，一上实际系统就振荡，八成是执行器延迟没建模。鼓风机从收到信号到实际出气量变化，通常有几十秒到几分钟的滞后。代码里actuator_delay.m提供了一个一阶惯性加纯延迟的模型，我一般会把延迟时间设成实际设备说明书上的响应时间，惯性时间常数取延迟的三分之一到二分之一。
% actuator_delay.m 片段 function u_actual = apply_delay(u_cmd, delay_time, tau, dt, state) % 一阶惯性 + 纯延迟 % state: 保存延迟队列和当前输出 % 延迟队列 state.buffer = [state.buffer(2:end); u_cmd]; u_delayed = state.buffer(1); % 一阶惯性 u_actual = state.y + (dt/tau) * (u_delayed - state.y); state.y = u_actual; enddelay_time和tau这两个参数如果拍脑袋填，控制效果会差很多。我的做法是：先做一次阶跃响应实验，记录从给信号到溶解氧开始变化的时间，那就是纯延迟；然后看溶解氧变化到新稳态的 63.2% 用了多久，减去纯延迟，就是惯性时间常数。没有实验条件的话，保守一点，延迟取 2 分钟，惯性取 1 分钟，至少不会振荡。
5. 避坑与排查：那些让我熬夜的报错和异常
5.1 求解器报错“无法满足积分容差”
现象：运行ode45时提示Unable to meet integration tolerances without reducing the step size。原因通常是速率方程在某个参数组合下出现了刚性，比如硝化速率常数辨识得太大，导致氨氮浓度在一步内变成负数，下一时刻又反弹。解决：换ode15s求解器，或者在速率方程里加非负约束，浓度低于零时强制截断。
5.2 参数辨识结果物理意义不合理
现象：辨识出的复氧系数是负值，或者温度修正系数大于 2。原因：初值离真实值太远，优化跑到了无物理意义的区域。解决：在lsqnonlin里设紧的上下界，复氧系数下界设 0.01，温度修正系数限制在 1.0 到 1.1 之间。如果还是跑偏，检查输入输出数据的时间对齐，错位一小时就能让辨识结果完全失真。
5.3 预测曲线在雨天突变
现象：晴天预测很准，一到降雨事件就偏得离谱。原因：模型没有考虑降雨带来的面源污染冲击和稀释效应。解决：在输入向量里加入降雨量作为前馈变量，或者在检测到降雨时切换到另一组辨识参数。代码里预留了rainfall_flag接口，但默认没启用，需要自己把降雨数据接进去。
5.4 MPC 控制量频繁抖动
现象：曝气量指令在相邻时刻反复跳变。原因：目标函数里控制增量惩罚系数太小，或者预测时域太短。解决：把增量惩罚系数从 0.1 调到 0.5 以上，同时把horizon从 10 增加到 15。如果还抖，检查状态反馈有没有噪声，加一个一阶低通滤波再送进控制器。
5.5 仿真速度慢到无法接受
现象：一次参数辨识跑半小时以上。原因：ode45在每次目标函数调用时都重新求解整个时间序列，而lsqnonlin默认要调用几百次。解决：把UseParallel打开，同时把仿真时间分段，用上一次的末状态作为下一段的初值，减少重复计算。数据量特别大时，先降采样到小时均值再辨识，精度损失很小但速度快一个数量级。
6. 进阶技巧：用滚动窗口辨识应对水质季节性变化
水质系统的参数不是恒定的。温度、微生物群落结构、底泥特性都会随季节变，一套参数用全年，春秋季的预测误差会明显大于冬夏。我后来养成的习惯是：每 30 天用最近 90 天的数据重新辨识一次参数，然后对比新旧参数的 NSE，如果新参数在验证集上 NSE 下降超过 0.05，就说明系统发生了结构性变化，需要检查是不是有新的排放源接入。
具体操作上，把parameter_estimation.m包一层循环，按时间窗口切片调用。注意每次辨识的初值用上一窗口的结果，这样收敛快，而且参数变化轨迹平滑，不会出现某个月突然跳变的情况。
% rolling_identification.m 片段 function p_history = rolling_identify(data_clean, window_days, step_days) % window_days: 辨识窗口长度 % step_days: 滚动步长 t_start = data_clean.time(1); t_end = data_clean.time(end); p_history = []; p_init = [0.3; 0.1; 0.5; 1.05]; % 初始猜测 current = t_start; while current + days(window_days) <= t_end idx = data_clean.time >= current & ... data_clean.time < current + days(window_days); data_win = data_clean(idx, :); p_est = estimate_params(data_win, p_init, ... [0.01; 0.01; 0.01; 1.0], [2.0; 1.0; 5.0; 1.1]); p_history = [p_history; current, p_est']; p_init = p_est; % 用当前结果作为下一窗口初值 current = current + days(step_days); end endwindow_days取 90 是因为水质系统的记忆长度大概在一个季度，太短辨识不稳定，太长跟不上变化。step_days取 30 是折中，每周滚动一次计算量太大，每季度滚动一次又太滞后。跑完这个循环，把p_history画出来，如果某个参数在特定月份持续漂移，那就是系统在告诉你：该去现场看看是不是有新的污染源了。
从那以后我每次拿到新的监测数据，都强制走一遍滚动辨识，哪怕只是看看参数有没有异常漂移。这套代码给的是一个骨架，真正让它在你手里跑出价值，靠的是把现场工况和模型参数对应起来的那点直觉。希望帮到你。
本文还有配套的精品资源，点击获取
