基于机器学习的生态组合塘强化城市污水处理厂脱氮优化

模型开发与评估

在模型开发阶段，预处理后的数据集按 8:2 的比例划分为训练集（80%）和测试集（20%）。为提升模型的泛化能力，研究首先通过 Pearson 相关性分析剔除预测增益低且存在共线性的特征，在减少特征维度的同时保留核心信息；同时，引入 SHAP（SHapley Additive exPlanations）值系统评估包括出水总氮在内的 16 个初始特征的重要性，剔除冗余项。基于对特征显著性、相关性及模型复杂度的综合平衡，最终构建了包含出水总氮在内的 13 个关键特征的精简数据集。

模型开发在 Anaconda 平台的 JupyterLab 环境中基于 Python 3.8 完成，并采用提前停止（Early Stopping）策略防止过拟合——若验证集均方误差连续 50 轮未改善则停止训练。通过网格搜索（Grid Search）进行超参数优化，并结合 5 折交叉验证评估模型性能以降低结果方差。评价指标选用决定系数和均方根误差，其中 RMSE 衡量预测值与观测值的差异，R² 反映拟合优度。研究对比了多元线性回归（MLR）、支持向量机（SVM）、轻量级梯度提升机（LightGBM）、随机森林（RF）、极端梯度提升（XGBoost）和类别特征提升（CatBoost）六种机器学习算法，旨在筛选出处理总氮去除任务的最佳模型。

总氮（TN）去除的优化策略

基于预测模型的输出结果，本研究通过调整工艺参数制定了针对性的脱氮优化策略。研究采用 SHAP 方法和部分依赖图（PDP）来评估最优模型中各特征的相对重要性。SHAP 是一种基于博弈论的机器学习解释方法，通过计算每个特征对模型输出的贡献度（即 SHAP 值），不仅能量化单个变量及其相互作用对预测值的正负影响，还能揭示这些影响在不同观测值间的变化规律。

为进一步衡量特征对模型整体预测的重要性，研究计算了平均绝对 Shapley（MAS）值。通过 SHAP 摘要图可以清晰阐明影响总氮去除的因素排序，而 PDP 则用于深入理解各特征的影响力及其相关性，并利用二维部分依赖图（2D-PDP）分析两个特征间的交互作用对预测结果的影响。

在研究实施中，利用 Python 的 shap 和 pdpbox 库进行数据分析与绘图，并将每个数据集对应的最优模型指定为解释器（Explainer）。研究流程在利用 PDP 分析出同时降低出水总氮与能耗的最优参数范围后，构建了用于出水总氮预测与优化的机器学习框架。通过 SHAP 分析确定的最优参数组合，旨在平衡脱氮效率、能耗及 COD 投加量。为便于实际应用，研究基于 tkinter 库开发了用户友好的图形用户界面（GUI）。此外，采用简化的带策略非支配排序遗传算法（NSGA-II）解决出水总氮、能耗与外加碳源消耗之间的多目标优化问题，并结合逼近理想解排序法（TOPSIS）获取最佳运行参数组合。

结果与讨论

基于出水总氮（TN）浓度的模型评估与比较

结果显示，MLR 的预测能力有限（测试集 R² = 0.459），而其余五种非线性模型的性能表现优异，测试集 R² 均在 0.775 至 0.829 之间。其中，XGBoost 模型表现最为突出，训练集准确率接近完美（R² = 0.999，RMSE = 0.073），但在测试集上的泛化能力略弱（R² = 0.829）。

随后，研究通过 SHAP 分析对特征贡献度进行了量化。结果发现，进水 pH 值和进水氨氮（NH₄⁺-N）的特征重要性最低（0.020），对模型性能提升的贡献微乎其微。考虑到进水 pH 值长期稳定在 6.5–7.5 之间，将其判定为冗余特征并予以剔除；而进水氨氮因与出水总氮存在潜在关联，故予以保留。

通过 Pearson 相关系数分析各特征与出水总氮的相关性发现，出水 COD（0.24）和废水温度（-0.26）与出水总氮的相关性较弱。此外，出水 COD 浓度远低于国家排放标准，对出水总氮浓度的影响有限，且在实际运行中，废水温度的调节空间较小。因此，为了降低模型复杂度并增强泛化能力，研究最终剔除了废水温度、进水 pH 和出水 COD 这三个冗余特征。

经过特征修剪后，对各模型进行了重新训练。优化后的 XGBoost 模型性能显著提升，训练集和测试集的 R²分别达到 0.997 和 0.911，RMSE 分别为 0.196 和 1.283。结果表明，由 SHAP 特征重要性评估引导的降维策略极大地改善了模型表现，最优模型（XGBoost）的测试集 R²从 0.829 提高到了 0.911。

图 2f 显示了 XGBoost 模型的拟合曲线及其 95% 预测区间。基于 13 个关键特征的组合，XGBoost 模型在所有模型中脱颖而出，成为预测总氮去除效果的最佳选择。由于绝大多数测试数据点均落在 95% 置信区间内，充分证明了该模型的可靠性与鲁棒性。

图 2. 以出水总氮浓度为目标特征的六种机器学习模型训练图：(a) 多元线性回归 (MLR)；(b) 支持向量机 (SVM)；(c) 随机森林 (RF)；(d) 轻量级梯度提升机 (LightGBM)；(e) 类别特征提升 (CatBoost)；(f) 极端梯度提升 (XGBoost)。图中红线表示测试集的回归拟合线，粉色区域则代表这些回归线周围的 95% 置信区间。

基于 SHAP 和 PDP 的模型灵敏度分析

研究采用 SHAP 和 PDP（部分依赖图）方法分析特征重要性，阐明 XGBoost 模型的预测机制，并识别影响出水总氮浓度的关键参数。结果显示，TN_7、EC（能耗）、DO_5、TN_6、进水流量以及 DO_7 是决定出水总氮浓度的最重要特征。其中，TN_7 和 TN_6 为水质指标，而 EC、DO_5 和 DO_7 则是具有可控潜力的工艺参数。值得注意的是，DO_5 和 DO_7 与能耗（EC）显著相关，这为后续制定智能节能策略提供了逻辑基础。

为了进一步探究特征对模型输出的具体影响，研究利用 PDP 考察了单特征及多特征交互作用对出水总氮的影响。出水总氮受硝化和反硝化过程的共同调节。DO_5 部分依赖图显示，维持较高的 DO_5 浓度有利于提高总氮去除效率；5 号区作为好氧区，较低的溶解氧会抑制微生物的硝化作用。然而，由于该工艺中 5 号区的曝气兼具搅拌功能，导致其 DO_5 浓度时常超过 10 mg/L，这凸显了利用机器学习优化 DO_5 以降低能耗的必要性。

DO_7 与出水总氮浓度呈正相关。这是因为 7 号区为缺氧区，较高的溶解氧会干扰反硝化效率。同时，COD_6 和 COD_7 的 PDP 结果表明，适量的 COD 投加能增强脱氮性能，但过量投加则会增加废水负荷与化学品消耗，反而对反硝化产生不利影响。

在多特征交互分析方面，表明较低的 DO_7 与较高的 DO_5 结合可实现更高的脱氮效率。外加碳源（ECS）与 COD_6、COD_7 之间的转换关系，结果暗示较低浓度的 ECS 和 COD 组合有助于降低出水总氮。综上所述，强化脱氮效率的最优运行区间为：DO_5 控制在 5.8–8.8 mg/L，DO_7 在 0.3–2.8 mg/L，COD_6 在 28.0–35.5 mg/L，COD_7 在 0–26.4 mg/L，且 ECS 投加量控制在 1.9–2.3 m³/d。

识别并排序影响脱氮性能的关键参数

污水处理厂的生态组合塘（ECPs）旨在通过好氧硝化和缺氧反硝化实现脱氮。如图 3 所示，TN_7、EC（能耗）和 DO_5 是影响脱氮效率的三大核心参数。由于 5 号区作为最后的好氧区且溶解氧（DO）未受控制，导致后续的 6 号区（第一缺氧区）DO 浓度升高，即便添加外加碳源也难以发生反硝化作用；同理，8 号区（第三缺氧区）因缺乏碳源也无法进行反硝化，而 9-12 号好氧区主要负责将残余氨氮转化为硝态氮，对总氮去除的贡献极小。相比之下，7 号区（第二缺氧区）通过控制 DO 浓度并添加外加碳源，促进了显著的反硝化过程，因此 TN_7 成为影响出水总氮浓度最关键的参数。

此外，5 号区和 7 号区不受控的 DO 浓度导致了过度曝气，从而增加了电力消耗；同时，在进水水质波动的情况下，外加碳源的不规范投加也造成了化学品的过量使用，这两者共同推高了整体能耗，使 EC 成为影响脱氮效率的第二大关键参数。值得注意的是，5 号区的最高 DO 浓度超过了 9.0 mg/L，这种缺乏控制的过度曝气使得本应作为第一缺氧区的 6 号区呈现出好氧特性，这也解释了为何 DO_5 会成为影响脱氮性能的第三大重要参数。

脱氮效果、能耗与 COD 投加量之间的关系

在基于生态组合塘（ECP）的污水处理厂中，曝气和碳源通过直接影响微生物的硝化与反硝化过程，显著决定了总氮（TN）的去除效率。利用二维部分依赖图（2D-PDP）揭示多特征间潜在关系的优势，展示了 DO_5 与 COD_6 对出水总氮（TN_eff）的共同影响：当 DO_5 维持在 6.3–8.8 mg/L 且 COD_6 处于 29–35 mg/L 范围时，预测的 TN_eff 达到最小值；虽然更高的 DO_5 浓度（>9.5 mg/L）也能增强脱氮，但会导致过度曝气并显著增加电耗。

溶解氧被认为是影响 ECP 工艺污水处理厂能耗最关键的因素。证明，在 DO_5 处于 6.3–8.8 mg/L 时，出水总氮不仅能保持在较低水平，且能在较低能耗下实现高效脱氮。当 DO_7 在 0.5–3.0 mg/L 范围内时，无论能耗高低均可实现有效脱氮，这预示了该厂具有巨大的节能潜力。由于最优参数范围的研究表明能耗受曝气及其他参数变化的共同影响，研究进一步开发了六种机器学习模型来预测能耗。对比结果显示，SVM 表现较差（测试集 R²=0.523），而 XGBoost 模型凭借最高的 R²（0.915）和最低的测试集 RMSE（883.446）再次表现最优。因此，XGBoost 预测模型被认为能够精准捕捉运行参数变化驱动的能耗波动。

强化脱氮与平衡能耗及 COD 投药量的参数优化

机器学习模型的研究结果强调了优化工艺参数以同步降低出水总氮（TN_eff）浓度、能耗（EC）和 COD 投药量的重要性。为了便于实际工程应用，本研究开发了一个集成 TN_eff、EC 和 COD 投药量模型的图形用户界面（GUI）。

该 GUI 基于 PDP 分析推导出的最优控制参数区间，并集成了简化版的 NSGA-II（带策略非支配排序遗传算法）和 TOPSIS（逼近理想解排序法）算法。通过该优化控制策略，用户可以灵活调整关键运行参数；在输入进水水质数据（如进水氨氮、进水总氮、5/6/7 号区总氮浓度）及核心工艺参数（如进水流量、5/7 号区溶解氧、6/7 号区 COD 以及外加碳源量）后，GUI 能够自动预测并优化出水总氮、能耗和外加碳源（ECS）的使用。

水质参数与可控参数范围

Water Quality Parameters（水质参数）

Controllable Parameters（可控参数）

Suggested Parameters（建议参数）

操作按钮与结果输出

操作按钮

• Predict and Optimize（预测与优化）：绿色按钮
• Reset（重置）：橙色按钮

结果输出区域

• Original Effluent TN (mg/L)：原始出水总氮
• Optimized Effluent TN (mg/L)：优化后出水总氮
• Original Energy Consumption (kWh)：原始能耗
• Optimized Energy Consumption (kWh)：优化后能耗
• Energy Savings (kWh)：节能量
• Carbon Source Savings (m³/d)：碳源节省量
• Optimal Parameters Combination:（最优参数组合）：下方预留显示区域

图 S8 展示了该 GUI 系统在实际应用场景下的分析与优化策略。研究收集了未参与模型训练和测试的独立数据，通过 GUI 进行了进一步的模型验证。为了在强化脱氮的同时兼顾能耗与 COD 投药量的减少，建议根据 PDP 识别出的区间进行工艺优化。例如，通过优化关键参数，出水总氮浓度从 14.29 mg/L 降至 11.59 mg/L，污水处理厂能耗每天减少了 128.66 kWh，COD 投加浓度降低了 34.29 mg/L，对应外加碳源投加量从 3.72 m³/d 降至 2.25 m³/d（以 1200 kg COD/m³ ECS 计）。随后，GUI 基于 365 天内不同特性的进水数据自动计算了脱氮提升、能耗及 COD 减量效果，实现了出水总氮、能耗和外加碳源使用的同步降低。值得注意的是，全年的出水总氮浓度平均下降了 17.50%，而 COD 投药量减少了 33.29%。

图 6b 展示了 365 天内优化前后的出水总氮浓度对比，优化后的总氮值始终低于优化前，表明 GUI 提供的参数组合有效提升了脱氮效果。能耗是污水处理厂脱氮优化中的关键因素。尽管已有研究建立了出水总氮和能耗的预测模型，但未能实现两者的同步优化。本研究采用 NSGA-II 模型对出水化学需氧量和总能耗进行多目标优化，实现了水质改善与节能的协同增效。结果表明，通过调整工艺运行参数，可以有效平衡能耗与 COD 投药量的削减，实现协调优化。这些研究成果证明，该优化控制策略能有效提升采用生态组合塘（ECPs）工艺的污水处理厂的脱氮性能，同时显著降低运行能耗和外加碳源的成本。

代码实现

模型开发与可视化

#!/usr/bin/env python # coding: utf-8 
# 导入机器学习、数据处理及可视化所需的库
from sklearn.model_selection import KFold, GridSearchCV, train_test_split
from sklearn.metrics import mean_squared_error, r2_score
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.preprocessing import StandardScaler
import pickle
from xgboost import XGBRegressor
import shap

# --- 1. 数据加载与初步切分 ---
file_path = r'data.csv'
df = pd.read_csv(file_path, encoding='GBK')
X = df.drop(columns='Effluent TN')
y = df['Effluent TN']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.20, random_state=42)

# --- 2. 特征标准化 ---
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
with open('scaler.pkl', 'wb') as f:
    pickle.dump(scaler, f)

# --- 3. 超参数优化 (Grid Search) ---
param_grid = {
    'n_estimators': [100, 300, 500, 1000],
    'max_depth': [3, 5, 7, 9, 12],
    'learning_rate': [0.01, 0.05, 0.1, 0.2],
    'subsample': [0.6, 0.8, 1.0],
    'colsample_bytree': [0.6, 0.8, 1.0],
    'gamma': [0, 1, 5],
    'reg_alpha': [0, 0.1, 1],
    'reg_lambda': [1, 1.5, 2]
}
base_model = XGBRegressor(objective='reg:squarederror', random_state=42)
grid_search = GridSearchCV(
    estimator=base_model,
    param_grid=param_grid,
    scoring='neg_mean_squared_error',
    cv=5,
    verbose=2,
    n_jobs=-1
)
grid_search.fit(X_train_scaled, y_train)
best_params = grid_search.best_params_
print("\n找到的最佳参数组合:")
print(best_params)

# --- 4. 交叉验证训练与提前停止策略 ---
xgboost_model = XGBRegressor(**best_params, objective='reg:squarederror', random_state=42)
kf = KFold(n_splits=5, shuffle=True, random_state=42)
train_r2_scores, train_rmse_scores = [], []
test_r2_scores, test_rmse_scores = [], []
all_y_train_pred, all_y_test_pred = [], []
all_y_train_true, all_y_test_true = [], []

for train_index, test_index in kf.split(X_train_scaled):
    X_train_fold, X_test_fold = X_train_scaled[train_index], X_train_scaled[test_index]
    y_train_fold, y_test_fold = y_train.iloc[train_index], y_train.iloc[test_index]
    X_train_part, X_val_part, y_train_part, y_val_part = train_test_split(
        X_train_fold, y_train_fold, test_size=0.2, random_state=42)
    xgboost_model.fit(X_train_part, y_train_part, eval_set=[(X_val_part, y_val_part)], early_stopping_rounds=50, verbose=False)
    best_iteration = xgboost_model.best_iteration
    xgboost_model.n_estimators = best_iteration
    xgboost_model.fit(X_train_fold, y_train_fold)
    y_train_pred = xgboost_model.predict(X_train_fold)
    y_test_pred = xgboost_model.predict(X_test_fold)
    all_y_train_pred.extend(y_train_pred)
    all_y_test_pred.extend(y_test_pred)
    all_y_train_true.extend(y_train_fold.values)
    all_y_test_true.extend(y_test_fold.values)
    train_r2_scores.append(r2_score(y_train_fold, y_train_pred))
    train_rmse_scores.append(np.sqrt(mean_squared_error(y_train_fold, y_train_pred)))
    test_r2_scores.append(r2_score(y_test_fold, y_test_pred))
    test_rmse_scores.append(np.sqrt(mean_squared_error(y_test_fold, y_test_pred)))

print(f"\n训练集平均结果：R²: {np.mean(train_r2_scores):.3f}, RMSE: {np.mean(train_rmse_scores):.3f}")
print(f"测试集平均结果：R²: {np.mean(test_r2_scores):.3f}, RMSE: {np.mean(test_rmse_scores):.3f}")

# --- 5. 结果性能评估 ---
plt.figure(figsize=(10, 6), dpi=600)
plt.scatter(all_y_train_true, all_y_train_pred, edgecolors='black', c='darkgreen', marker='^', s=100, alpha=0.6, label='训练集')
plt.scatter(all_y_test_true, all_y_test_pred, edgecolors='black', c='lightyellow', marker='o', s=100, alpha=0.6, label='测试集')
X_plot = np.array(all_y_test_true).reshape(-1, 1)
y_plot_pred = np.array(all_y_test_pred)
sorted_idx = np.argsort(X_plot.flatten())
X_sorted = X_plot[sorted_idx]
y_sorted = y_plot_pred[sorted_idx]
linear_model = sm.OLS(y_sorted, sm.add_constant(X_sorted)).fit()
y_pred_line = linear_model.predict(sm.add_constant(X_sorted))
predictions_summary = linear_model.get_prediction(sm.add_constant(X_sorted)).summary_frame(alpha=0.05)
plt.plot(X_sorted, y_pred_line, color='red', linewidth=2, label='拟合回归线')
plt.fill_between(X_sorted.flatten(), predictions_summary['obs_ci_lower'], predictions_summary['obs_ci_upper'], color='lightpink', alpha=0.3, label='95% 置信区间')
plt.plot([0, 17.5], [0, 17.5], c='blue', linestyle='--', label='理想线 (y=x)')
plt.legend(); plt.show()

# --- 6. 模型解释性分析 (SHAP) ---
explainer = shap.Explainer(xgboost_model, X_train_scaled)
shap_values = explainer(X_train_scaled)
shap.summary_plot(shap_values, X_train_scaled, feature_names=X_train.columns)

模型可解释性

SHAP 可解释性

#!/usr/bin/env python # coding: utf-8
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from xgboost import XGBRegressor
import numpy as np
import matplotlib.pyplot as plt
import shap

file_path = r'D:\vscode-water\TNeff-GUI-master\Modelset\Data.csv'
df = pd.read_csv(file_path, encoding='GBK')
X = df.drop(columns='Effluent TN')
y = df['Effluent TN']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

xgboost_model = XGBRegressor(
    objective='reg:squarederror',
    n_estimators=1000,
    max_depth=6,
    learning_rate=0.01,
    subsample=0.8,
    colsample_bytree=0.8,
    random_state=42
)
xgboost_model.fit(X_train, y_train)

feature_importances = xgboost_model.feature_importances_
sorted_idx = np.argsort(feature_importances)[::-1]
plt.figure(figsize=(12, 8))
plt.bar(range(len(feature_importances)), feature_importances[sorted_idx], color='deepskyblue', alpha=0.3)
plt.xticks(range(len(feature_importances)), [X.columns[i] for i in sorted_idx], rotation=90, fontsize=18)
plt.title('Feature Importance', fontsize=20)
plt.xlabel('Features', fontsize=20)
plt.ylabel('Importance', fontsize=20)
plt.show()

explainer = shap.TreeExplainer(xgboost_model)
shap_values = explainer.shap_values(X_train)
shap.summary_plot(shap_values, X_train, feature_names=X.columns)

单变量的 PDP

#!/usr/bin/env python # coding: utf-8
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.inspection import partial_dependence
from scipy.interpolate import splev, splrep

file_path = r'D:\vscode-water\TNeff-GUI-master\Modelset\Data.csv'
df = pd.read_csv(file_path, encoding='GBK')
X = df.drop(columns='Effluent TN')
y = df['Effluent TN']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

xgboost_model = XGBRegressor(
    objective='reg:squarederror',
    n_estimators=1000,
    max_depth=6,
    learning_rate=0.01,
    subsample=0.8,
    colsample_bytree=0.8,
    random_state=42
)
xgboost_model.fit(X_train, y_train)

features_to_plot = [
    'External Carbon Source', 'COD of Zone 6', 'COD of Zone 7',
    'DO of Zone 7', 'DO of Zone 5', 'Influent NH4+-N',
    'Influent TN', 'Influent Flowrate', 'TN of Zone 5',
    'TN of Zone 6', 'TN of Zone 7', 'Energy Consumption'
]

sns.set_theme(style="ticks", palette="deep", font_scale=1.1)
def plot_pdp(feature):
    pdp = partial_dependence(xgboost_model, X_train, [feature], kind="both", grid_resolution=50)
    plot_x = pd.Series(pdp.grid_values[0]).rename('x')
    plot_y = pd.Series(pdp.average[0]).rename('y')
    plot_i = pdp.individual[0]
    tck = splrep(plot_x, plot_y, s=30)
    xnew = np.linspace(plot_x.min(), plot_x.max(), 300)
    ynew = splev(xnew, tck, der=0)
    fig, ax = plt.subplots(figsize=(8, 6))
    for a in plot_i:
        a_series = pd.Series(a)
        df_i = pd.concat([plot_x, a_series.rename('y')], axis=1)
        sns.lineplot(data=df_i, x="x", y="y", color='k', linewidth=1.5, linestyle='--', alpha=0.6, ax=ax)
    ax.plot(xnew, ynew, color='peru', linewidth=2, label='Smoothed PDP')
    std_error = np.std(plot_y) / np.sqrt(len(plot_y))
    lower_bound = plot_y - 1.96 * std_error
    upper_bound = plot_y + 1.96 * std_error
    ax.fill_between(plot_x, lower_bound, upper_bound, color='khaki', alpha=0.3, label='95% CI')
    sns.rugplot(data=X_train.sample(100), x=feature, height=0.05, color='k', alpha=0.3, ax=ax)
    ax.set_ylabel('Partial Dependence')
    ax.set_xlabel(feature)
    x_min = plot_x.min() - 0.1*(plot_x.max() - plot_x.min())
    x_max = plot_x.max() + 0.1*(plot_x.max() - plot_x.min())
    ax.set_xlim(x_min, x_max)
    ax.legend()
    plt.savefig(f'./pdpplot_{feature}.png', dpi=900, bbox_inches='tight')
    plt.show()

for feature in features_to_plot:
    plot_pdp(feature)

双变量的 PDP

#!/usr/bin/env python # coding: utf-8
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.inspection import PartialDependenceDisplay

file_path = r'D:\vscode-water\TNeff-GUI-master\Modelset\Data.csv'
df = pd.read_csv(file_path, encoding='GBK')
X = df.drop(columns='Effluent TN')
y = df['Effluent TN']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

xgboost_model = XGBRegressor(
    objective='reg:squarederror',
    n_estimators=1000,
    max_depth=6,
    learning_rate=0.01,
    subsample=0.8,
    colsample_bytree=0.8,
    random_state=42
)
xgboost_model.fit(X_train, y_train)

target_feature = 'External Carbon Source'
related_features = ['COD of Zone 6', 'COD of Zone 7']

for feature in related_features:
    if feature == target_feature:
        continue
    feature_idx = X.columns.get_loc(feature)
    target_idx = X.columns.get_loc(target_feature)
    plt.figure(figsize=(8, 6))
    PartialDependenceDisplay.from_estimator(
        xgboost_model, X_train, [(target_idx, feature_idx)],
        feature_names=X.columns, kind='average'
    )
    plt.title(f'2D Partial Dependence: {target_feature}_feature vs {feature}')
    plt.show()

TOPSIS 算法

import numpy as np

def topsis(data, weights):
    if len(data) == 0:
        return np.array([])
    if data.ndim == 1:
        data = data.reshape(1, -1)
    norm_data = (data - data.min(axis=0)) / (data.max(axis=0) - data.min(axis=0) + 1e-10)
    weighted_data = norm_data * weights
    ideal_best = np.min(weighted_data, axis=0)
    ideal_worst = np.max(weighted_data, axis=0)
    dist_best = np.sqrt(np.sum((weighted_data - ideal_best) ** 2, axis=1))
    dist_worst = np.sqrt(np.sum((weighted_data - ideal_worst) ** 2, axis=1))
    scores = dist_worst / (dist_best + dist_worst + 1e-10)
    return scores

# 示例调用
data = np.array([
    [12.5, 300, 50],
    [10.0, 450, 80],
    [8.5, 600, 120],
    [11.0, 350, 60]
])
weights = np.array([0.5, 0.3, 0.2])
scores = topsis(data, weights)
print("各策略的 TOPSIS 综合得分 (得分越高越优):")
for i, score in enumerate(scores):
    print(f"策略 {i+1}: {score:.4f}")
best_index = np.argmax(scores)
print(f"\n最佳推荐方案是：策略 {best_index + 1}")

NSGA-II 优化

import numpy as np
import random
from deap import base, creator, tools, algorithms
import xgboost as xgb
import pandas as pd
from sklearn.preprocessing import StandardScaler

# 默认工况数据
default_data = {
    'Influent NH4+-N': 29.6,
    'Influent TN': 35.9,
    'Influent Flowrate': 38808.4,
    'DO of Zone 5': 7.1,
    'TN of Zone 5': 22.2,
    'COD of Zone 6': 55.0,
    'TN of Zone 6': 13.8,
    'DO of Zone 7': 4.06,
    'COD of Zone 7': 31.9,
    'TN of Zone 7': 11.7,
    'External Carbon Source': 3.3
}

variable_ranges = [
    (5.8, 7.5),
    (28.0, 35.5),
    (0.3, 2.8),
    (1.0, 26.4),
    (1.9, 2.3)
]

df = pd.DataFrame([default_data])
scaler = StandardScaler()
scaler.fit(df)

model_tn = xgb.Booster()
model_power = xgb.Booster()
prediction_cache = {}

def predict_outputs(data):
    if data.ndim == 1:
        key = tuple(data)
        data = data.reshape(1, -1)
    else:
        return model_power.predict(xgb.DMatrix(data)), model_tn.predict(xgb.DMatrix(data))
    if key in prediction_cache:
        return np.array([prediction_cache[key][0]]), np.array([prediction_cache[key][1]])
    dmatrix = xgb.DMatrix(data)
    pred_energy = model_power.predict(dmatrix)
    pred_tn = model_tn.predict(dmatrix)
    prediction_cache[key] = (pred_energy[0], pred_tn[0])
    return pred_energy, pred_tn

def evaluate(individual, fixed_water_params, true_energy, true_tn):
    full_input = np.concatenate([fixed_water_params, individual])
    energy, tn = predict_outputs(np.array(full_input))
    delta_energy = energy[0] - true_energy
    delta_tn = tn[0] - true_tn
    return delta_energy, delta_tn

def optimize_parameters(fixed_water_params, initial_operating_params):
    if not hasattr(creator, 'FitnessMin'):
        creator.create("FitnessMin", base.Fitness, weights=(-1.0, -1.0))
    if not hasattr(creator, 'Individual'):
        creator.create("Individual", list, fitness=creator.FitnessMin)
    toolbox = base.Toolbox()
    toolbox.register("attr_float", lambda: [random.uniform(r[0], r[1]) for r in variable_ranges])
    toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.attr_float)
    toolbox.register("population", tools.initRepeat, list, toolbox.individual)
    toolbox.register("mate", tools.cxSimulatedBinaryBounded, low=[r[0] for r in variable_ranges], up=[r[1] for r in variable_ranges], eta=15.0)
    def custom_mutate(individual):
        individual = tools.mutPolynomialBounded(individual, low=[r[0] for r in variable_ranges], up=[r[1] for r in variable_ranges], eta=20.0, indpb=0.3)[0]
        for i in range(len(individual)):
            low, up = variable_ranges[i]
            individual[i] = max(low, min(up, individual[i]))
            if isinstance(individual[i], complex) or np.isnan(individual[i]) or np.isinf(individual[i]):
                individual[i] = (low + up) / 2
        return individual
    toolbox.register("mutate", custom_mutate)
    toolbox.register("select", tools.selNSGA2)
    full_input = np.concatenate([fixed_water_params, initial_operating_params])
    true_energy, true_tn = predict_outputs(full_input)
    def fitness_func(ind):
        return evaluate(ind, fixed_water_params, true_energy[0], true_tn[0])
    toolbox.register("evaluate", fitness_func)
    pop = toolbox.population(n=100)
    pop[0][:] = initial_operating_params.copy()
    for ind in pop:
        ind.fitness.values = toolbox.evaluate(ind)
    hof = tools.ParetoFront()
    hof.update(pop)
    best_fitness = min([ind.fitness.values[0] for ind in pop])
    stall_count = 0
    for gen in range(20):
        offspring = algorithms.varOr(pop, toolbox, lambda_=100, cxpb=0.5, mutpb=0.5)
        for ind in offspring:
            if not ind.fitness.valid:
                ind.fitness.values = toolbox.evaluate(ind)
        pop = toolbox.select(pop + offspring, k=20)
        hof.update(pop)
        current_best = min([ind.fitness.values[0] for ind in pop])
        if current_best < best_fitness * 0.99:
            best_fitness = current_best
            stall_count = 0
        else:
            stall_count += 1
            if stall_count >= 5:
                break
    n_samples = min(60, len(hof))
    pareto_samples = [hof[i] for i in np.linspace(0, len(hof) - 1, n_samples, dtype=int)]
    return pareto_samples, true_energy[0], true_tn[0]

图形用户界面 (GUI)

#!/usr/bin/env python # coding: utf-8
import tkinter as tk
from tkinter import ttk
import xgboost as xgb
import pandas as pd
from sklearn.preprocessing import StandardScaler

default_data = {
    'Influent NH4+-N': 29.6,
    'Influent TN': 35.9,
    'Influent Flowrate': 38808.4,
    'DO of Zone 5': 7.1,
    'TN of Zone 5': 22.2,
    'COD of Zone 6': 55.0,
    'TN of Zone 6': 13.8,
    'DO of Zone 7': 4.06,
    'COD of Zone 7': 31.9,
    'TN of Zone 7': 11.7,
    'External Carbon Source': 3.3
}

df = pd.DataFrame([default_data])
scaler = StandardScaler()
scaler.fit(df)

model_tn = xgb.Booster()
model_tn.load_model('model_effluent_TN.json')
model_power = xgb.Booster()
model_power.load_model('model_energy_consumption.json')

def validate_input(value):
    try:
        return float(value)
    except ValueError:
        return None

run_count = 0
first_external_carbon_source = None

def clear_inputs():
    global run_count, first_external_carbon_source
    run_count = 0
    first_external_carbon_source = None
    for entry in feature_entry_map.values():
        entry.delete(0, tk.END)
    original_tn.set("")
    original_power.set("")
    optimized_tn.set("")
    optimized_power.set("")
    energy_savings_var.set("")
    carbon_savings_var.set("")
    output_text.set("")

def predict():
    global run_count, first_external_carbon_source
    if run_count >= 2:
        output_text.set("The prediction limit has been reached. Please click the Reset button to reset and try again.")
        return
    inputs = []
    for feature, entry in feature_entry_map.items():
        value = entry.get().strip()
        if not value:
            inputs.append(default_data[feature])
        else:
            validated = validate_input(value)
            if validated is None:
                output_text.set(f"Invalid input for {feature}. Please enter a valid number.")
                return
            inputs.append(validated)
    try:
        inputs_df = pd.DataFrame([inputs], columns=default_data.keys())
        inputs_scaled = scaler.transform(inputs_df)
        dmatrix = xgb.DMatrix(inputs_scaled)
        effluent_tn_pred = model_tn.predict(dmatrix)[0]
        power_consumption_pred = model_power.predict(dmatrix)[0]
        if run_count == 0:
            original_tn.set(f"{effluent_tn_pred:.4f}")
            original_power.set(f"{power_consumption_pred:.4f}")
            first_external_carbon_source = inputs_df["External Carbon Source"].iloc[0]
        elif run_count == 1:
            optimized_tn.set(f"{effluent_tn_pred:.4f}")
            optimized_power.set(f"{power_consumption_pred:.4f}")
            original_power_value = float(original_power.get()) if original_power.get() else 0.0
            energy_savings = original_power_value - power_consumption_pred
            energy_savings_var.set(f"{energy_savings:.4f}")
            second_external_carbon_source = inputs_df["External Carbon Source"].iloc[0]
            carbon_savings = first_external_carbon_source - second_external_carbon_source
            carbon_savings_var.set(f"{carbon_savings:.4f}")
            run_count += 1
    except Exception as e:
        output_text.set(f"An error occurred during prediction:\n{e}")

units = {
    "Influent NH4+-N": "5.4 - 60.7 mg/L",
    "Influent TN": "9.8 - 108.0 mg/L",
    "TN of Zone 5": "4.8 - 44.6 mg/L",
    "TN of Zone 6": "3.0 - 42.3 mg/L",
    "TN of Zone 7": "2.8 - 38.9 mg/L",
    "Influent Flowrate": "0 - 69970.0 m³/d",
    "DO of Zone 5": "5.8 - 8.8 mg/L",
    "COD of Zone 6": "28.0 - 35.5 mg/L",
    "DO of Zone 7": "0.3 - 2.8 mg/L",
    "COD of Zone 7": "0 - 26.4 mg/L",
    "External Carbon Source": "1.9 - 2.3 m³/d"
}

root = tk.Tk()
root.title("Effluent TN Prediction and Optimization")
root.geometry("1200x800")
root.configure(bg="#e8f0f2")
default_font = ("Arial", 12)
main_frame = tk.Frame(root, bg="#e8f0f2", padx=20, pady=20)
main_frame.pack(expand=True)
title_label = tk.Label(main_frame, text="Effluent TN Prediction and Optimization", font=("Arial", 18, "bold"), bg="#e8f0f2")
title_label.grid(row=0, column=0, columnspan=8, pady=15)

left_labels = [
    ("Influent NH₄⁺-N", "Influent NH4+-N"),
    ("Influent TN", "Influent TN"),
    ("TN of Zone 5", "TN of Zone 5"),
    ("TN of Zone 6", "TN of Zone 6"),
    ("TN of Zone 7", "TN of Zone 7")
]
right_labels = [
    ("Influent Flowrate", "Influent Flowrate"),
    ("DO of Zone 5", "DO of Zone 5"),
    ("COD of Zone 6", "COD of Zone 6"),
    ("DO of Zone 7", "DO of Zone 7"),
    ("COD of Zone 7", "COD of Zone 7"),
    ("External Carbon Source", "External Carbon Source")
]

feature_entry_map = {}

tk.Label(main_frame, text="Water Quality Parameters", font=("Arial", 14, "bold"), bg="#e8f0f2").grid(row=1, column=1, columnspan=2, pady=10)
for i, (label, feature) in enumerate(left_labels, start=2):
    tk.Label(main_frame, text=label, font=default_font, bg="#e8f0f2").grid(row=i, column=0, padx=10, pady=5, sticky="e")
    entry = tk.Entry(main_frame, font=default_font, width=20)
    entry.grid(row=i, column=1, padx=10, pady=5, sticky="ew")
    feature_entry_map[feature] = entry
    unit_text = units.get(feature, "")
    tk.Label(main_frame, text=f"({unit_text})", font=default_font, bg="#e8f0f2").grid(row=i, column=2, padx=10, pady=5, sticky="w")

tk.Label(main_frame, text="Controllable Parameters", font=("Arial", 14, "bold"), bg="#e8f0f2").grid(row=1, column=3, columnspan=2, pady=10)
for i, (label, feature) in enumerate(right_labels, start=2):
    tk.Label(main_frame, text=label, font=default_font, bg="#e8f0f2").grid(row=i, column=3, padx=10, pady=5, sticky="e")
    entry = tk.Entry(main_frame, font=default_font, width=20)
    entry.grid(row=i, column=4, padx=10, pady=5, sticky="ew")
    feature_entry_map[feature] = entry
    unit_text = units.get(feature, "")
    tk.Label(main_frame, text=f"({unit_text})", font=default_font, bg="#e8f0f2").grid(row=i, column=5, padx=10, pady=5, sticky="w")

original_tn = tk.StringVar()
original_power = tk.StringVar()
optimized_tn = tk.StringVar()
optimized_power = tk.StringVar()
energy_savings_var = tk.StringVar()
carbon_savings_var = tk.StringVar()
output_text = tk.StringVar()

predict_button = tk.Button(main_frame, text="Predict", font=("Arial", 14), command=predict, bg="#4CAF50", fg="white", width=18)
predict_button.grid(row=8, column=0, columnspan=4, pady=20)
reset_button = tk.Button(main_frame, text="Reset", font=("Arial", 14), command=clear_inputs, bg="#FF5722", fg="white", width=18)
reset_button.grid(row=8, column=4, columnspan=4, pady=20)

output_frame = tk.Frame(main_frame, bg="#e8f0f2")
output_frame.grid(row=9, column=0, columnspan=8, padx=10, pady=20)

tk.Label(output_frame, text="Original Effluent TN (mg/L):", font=default_font, bg="#e8f0f2").grid(row=0, column=0, sticky="e")
tk.Entry(output_frame, textvariable=original_tn, font=default_font, state="readonly").grid(row=0, column=1)
tk.Label(output_frame, text="Optimized Effluent TN (mg/L):", font=default_font, bg="#e8f0f2").grid(row=0, column=2, sticky="e")
tk.Entry(output_frame, textvariable=optimized_tn, font=default_font, state="readonly").grid(row=0, column=3)

tk.Label(output_frame, text="Original Energy Consumption (kWh):", font=default_font, bg="#e8f0f2").grid(row=1, column=0, sticky="e")
tk.Entry(output_frame, textvariable=original_power, font=default_font).grid(row=1, column=1)
tk.Label(output_frame, text="Optimized Energy Consumption (kWh):", font=default_font, bg="#e8f0f2").grid(row=1, column=2, sticky="e")
tk.Entry(output_frame, textvariable=optimized_power, font=default_font, state="readonly").grid(row=1, column=3)

tk.Label(output_frame, text="Energy Savings (kWh):", font=default_font, bg="#e8f0f2").grid(row=2, column=0, sticky="e")
tk.Entry(output_frame, textvariable=energy_savings_var, font=default_font, state="readonly").grid(row=2, column=1)
tk.Label(output_frame, text="Carbon Source Savings (m³/d):", font=default_font, bg="#e8f0f2").grid(row=2, column=2, sticky="e")
tk.Entry(output_frame, textvariable=carbon_savings_var, font=default_font, state="readonly").grid(row=2, column=3)

output_label = tk.Label(output_frame, textvariable=output_text, font=default_font, bg="#e8f0f2", fg="#FF5722")
output_label.grid(row=3, column=0, columnspan=4)

root.mainloop()