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

此外，研究采用 KNN 算法对采集的数据集进行了缺失值处理与离群值清洗，剔除了占总数据集不足 1% 的包含离群值的记录。

模型开发与评估

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

模型开发在 Anaconda 平台的 JupyterLab 环境中基于 Python 3.8 完成，并采用提前停止（Early Stopping）策略防止过拟合。通过网格搜索（Grid Search）进行超参数优化，并结合 5 折交叉验证评估模型性能以降低结果方差。评价指标选用决定系数和均方根误差。研究对比了多元线性回归（MLR）、支持向量机（SVM）、轻量级梯度提升机（LightGBM）、随机森林（RF）、极端梯度提升（XGBoost）和类别特征提升（CatBoost）六种机器学习算法。

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

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

利用 Python 的 shap 和 pdpbox 库进行数据分析与绘图，并将每个数据集对应的最优模型指定为解释器。研究流程在利用 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）的特征重要性最低。考虑到进水 pH 值长期稳定在 6.5–7.5 之间，将其判定为冗余特征并予以剔除；而进水氨氮因与出水总氮存在潜在关联，故予以保留。

通过 Pearson 相关系数分析各特征与出水总氮的相关性发现，出水 COD 和废水温度与出水总氮的相关性较弱。因此，为了降低模型复杂度并增强泛化能力，研究最终剔除了废水温度、进水 pH 和出水 COD 这三个冗余特征。

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

图 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 则是具有可控潜力的工艺参数。

为了进一步探究特征对模型输出的具体影响，研究利用 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 组合有助于降低出水总氮。综上所述，强化脱氮效率的最优运行区间为：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 浓度升高，即便添加外加碳源也难以发生反硝化作用。相比之下，7 号区（第二缺氧区）通过控制 DO 浓度并添加外加碳源，促进了显著的反硝化过程，因此 TN_7 成为影响出水总氮浓度最关键的参数。

此外，5 号区和 7 号区不受控的 DO 浓度导致了过度曝气，从而增加了电力消耗；同时，在进水水质波动的情况下，外加碳源的不规范投加也造成了化学品的过量使用，这两者共同推高了整体能耗，使 EC 成为影响脱氮效率的第二大关键参数。

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

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

溶解氧被认为是影响 ECP 工艺污水处理厂能耗最关键的因素。图 5c 证明，在 DO_5 处于 6.3–8.8 mg/L 时，出水总氮不仅能保持在较低水平，且能在较低能耗下实现高效脱氮。图 5d 显示，当 DO_7 在 0.5–3.0 mg/L 范围内时，无论能耗高低均可实现有效脱氮，这预示了该厂具有巨大的节能潜力。

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

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

该 GUI 基于 PDP 分析推导出的最优控制参数区间，并集成了简化版的 NSGA-II 和 TOPSIS 算法。通过该优化控制策略，用户可以灵活调整关键运行参数；在输入进水水质数据及核心工艺参数后，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 计）。

图 6b 展示了 365 天内优化前后的出水总氮浓度对比，优化后的总氮值始终低于优化前，表明 GUI 提供的参数组合有效提升了脱氮效果。能耗是污水处理厂脱氮优化中的关键因素。尽管已有研究建立了出水总氮和能耗的预测模型，但未能实现两者的同步优化。本研究采用 NSGA-II 模型对出水化学需氧量和总能耗进行多目标优化，实现了水质改善与节能的协同增效。

补充材料

附录内容清单

• 文本 S1：机器学习模型详述
• 文本 S2：多目标优化框架的建立
• 文本 S3：碳减排量的计算方法
• 图 S1：生态组合塘（ECPs）的平面布局图
• 图 S2：以出水总氮浓度为目标特征（基于 16 个初始参数）的六种机器学习模型训练图
• 图 S3：(a) 用于预测出水总氮（TN_eff）的 XGBoost 模型 SHAP 特征重要性分析（16 个初始参数）；(b) 该模型的 SHAP 显著性分析
• 图 S4：特征与标签（目标变量）之间的 Pearson 相关性分析
• 图 S5：特征'外加碳源'的部分依赖图（PDP）
• 图 S6：其他特征的部分依赖图
• 图 S7：以能耗（EC）为目标特征的六种机器学习模型训练图
• 图 S8：用于预测与优化的图形用户界面（GUI）应用示例
• 表 S1：采用生态组合塘工艺的大型城市污水处理厂进水水质及设计出水水质指标
• 表 S2：训练及评估数据集的统计特性
• 表 S3：不同模型的性能指标对比
• 表 S4：基于 XGBoost 模型确定的强化脱氮且平衡能耗与 COD 减量的优化参数
• 表 S5：碳排放计算所采用的相关参数

文本 S1：机器学习模型详述

本研究对比了六种不同的机器学习模型，以评估其在污水处理厂总氮预测和能耗分析中的表现：

1. 多元线性回归 (MLR)：一种通过独立变量集捕捉因变量之间线性关系的统计方法。适用性：在变量间相互作用主要为线性时特别有效。
2. 支持向量机 (SVM)：其核心原则是识别能够最大化不同类别数据点之间间隔的最佳超平面。优势：在处理原始特征空间中难以解决的高维数据任务时表现出色。
3. 随机森林 (RF)：一种基于分类和决策树回归的集成学习算法。特性：相比传统装袋法，在处理小数据集和建模非线性关系方面更具优势。
4. 轻量级梯度提升机 (LightGBM)：由微软开发的基于直方图决策树的梯度提升框架。表现：训练速度快、准确度高，并能有效缓解模型过拟合问题。
5. 类别特征提升 (CatBoost)：基于梯度提升决策树（GBDT），专为有效处理异构数据和类别型数据而设计。稳定性：采用随机排列策略对数据排序，并利用权重系数减少低频和噪声数据的影响。
6. 极端梯度提升 (XGBoost)：集成多个梯度提升决策树，通过二阶梯度优化和正则化技术减少过拟合。适用范围：在小样本数据集上表现强劲，性能优于单体机器学习方法。

文本 S2：多目标优化框架的建立

本研究通过集成进化算法、评价模型与预测模型，构建了一个协同优化框架，旨在平衡出水水质与运行成本。

1. 多目标优化（NSGA-II）：采用简化的带策略非支配排序遗传算法，解决出水总氮、能耗与碳源消耗之间的多目标优化问题。通过优化溶解氧、各区域 COD 以及外加碳源等运行参数，协调各目标间的复杂关系。
2. TOPSIS 决策方法：利用逼近理想解排序法从 NSGA-II 生成的 Pareto 最优解集中筛选出最佳运行参数组合。基于目标函数值计算每个方案的 TOPSIS 评分，分值最高者即被确定为最优方案。
3. 模型预测与优化流程：采用 XGBoost 回归模型预测出水总氮和能耗。通过加载预训练模型，并对输入数据进行标准化处理，以确保预测的准确性。初始化种群并设置交叉、变异及选择等遗传操作。当最优解的适应度在连续若干代内提升比例低于设定阈值时，进化停止。从最终的 Pareto 前沿中抽取样本，通过 TOPSIS 计算得出最优解。
4. 图形用户界面 (GUI)：构建了图形化界面，方便用户输入进水水质和工艺参数。系统首先预测原始的出水总氮和能耗数据，随后运行多次优化算法以锁定最优解。界面最终呈现优化后的总氮浓度、能耗、节能效果、碳源节省量以及最优参数组合。

代码复现

污水处理数据集特征说明表

项目结构

TNeff-GUI-master/
├── Data/                          # 数据文件
│    ├── Data after deleting features.csv
│    └── Data.csv
├── Graphical User Interace (GUI)/ # 用户界面相关文件
│    ├── GUI.py                   # 基础 GUI 实现
│    ├── GUI_layout.py            # 优化版 GUI 布局
│    ├── NSGA-II.py               # 多目标优化算法
│    ├── TOPSIS.py                # 多准则决策方法
│    ├── model_effluent_TN.json     # 总氮预测模型
│    ├── model_energy_consumption.json # 能耗预测模型
│    ├── step 1.png                 # 使用步骤截图
│    └── step 2.png
├── Model Development - Visualization/ # 模型开发与可视化
│    ├── Catboost.py              # CatBoost 模型
│    ├── LGBM.py                  # LightGBM 模型
│    ├── MLR.py                   # 多元线性回归模型
│    ├── RF.py                    # 随机森林模型
│    ├── SVM.py                   # 支持向量机模型
│    └── XGBoost.py               # XGBoost 模型
├── Model Interpretability/        # 模型可解释性分析
│    ├── One-dimensional PDP.py   # 一维部分依赖图
│    ├── SHAP.py                  # SHAP 值分析
│    └── Two-dimensional PDP.py   # 二维部分依赖图
└── README.md                      # 项目说明文档

模型开发与可视化

#!/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. 数据加载与初步切分 ---
# 从 CSV 文件读取脱氮工艺数据集
file_path = r'data.csv'
df = pd.read_csv(file_path, encoding='GBK')
# 提取特征变量 (X) 和目标变量 (y，即出水 TN 浓度)
X = df.drop(columns='Effluent TN')
y = df['Effluent TN']
# 将数据划分为训练集 (80%) 和测试集 (20%)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.20, random_state=42)

# --- 2. 特征标准化 ---
# 使用 StandardScaler 进行标准化处理，确保不同量纲的参数对模型贡献公平
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# 将标准化模型保存为 pickle 文件，以便后续在 GUI 界面中调用
with open('scaler.pkl', 'wb') as f:
    pickle.dump(scaler, f)

# --- 3. 超参数优化 (Grid Search) ---
# 定义 XGBoost 的超参数搜索空间，包括树的数量、深度、学习率及正则化参数
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]
}
# 初始化 XGBoost 回归模型
base_model = XGBRegressor(objective='reg:squarederror', random_state=42)
# 使用网格搜索和 5 折交叉验证寻找最优超参数组合
grid_search = GridSearchCV(
    estimator=base_model,
    param_grid=param_grid,
    scoring='neg_mean_squared_error',
    cv=5,
    verbose=2,
    n_jobs=-1 # 使用所有 CPU 核心并行加速
)
# 执行网格搜索
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)
# 设定 5 折交叉验证，用于评估模型泛化能力并降低方差
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]
    # 将当前折数据再次划分，用于应用提前停止 (Early Stopping) 策略，防止过拟合
    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)
    # 训练模型，若验证集误差连续 50 轮不下降则提前停止
    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)
    # 计算本折的决定系数 R² 和均方根误差 RMSE
    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)))

# --- 5. 结果性能评估 ---
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}")

# --- 6. 模型拟合可视化 (观测值 vs 预测值) ---
# 绘制散点图，展示模型预测值与实际观测值的匹配度
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='测试集')
# 拟合并绘制回归线及 95% 置信区间 (CI)
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()

# --- 7. 模型解释性分析 (SHAP) ---
# 计算 SHAP 值，量化每个特征对模型输出的具体贡献度
explainer = shap.Explainer(xgboost_model, X_train_scaled)
shap_values = explainer(X_train_scaled)
# 绘制 SHAP 摘要图，识别影响 TN 去除的关键参数 (如 TN7, EC, DO5 等)
shap.summary_plot(shap_values, X_train_scaled, feature_names=X_train.columns)

模型可解释性

SHAP 可解释性

#!/usr/bin/env python # coding: utf-8 
# In[ ]:
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

# Load dataset
file_path = r'Data.csv'
df = pd.read_csv(file_path, encoding='GBK')
# Split features and target
X = df.drop(columns='Effluent TN')
y = df['Effluent TN']
# Train-test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Initialize XGBoost model
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
)
# Model training
xgboost_model.fit(X_train, y_train)
# Plot feature importance
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()
# SHAP value analysis
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 
# In[ ]:
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

# Read data from CSV file
file_path = r'Data.csv'
df = pd.read_csv(file_path, encoding='GBK')
# Split features and target
X = df.drop(columns='Effluent TN')
y = df['Effluent TN']
# Split dataset into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Initialize XGBoost model
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
)
# Train the model
xgboost_model.fit(X_train, y_train)
# List of features to plot
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'
]
# Set Seaborn theme
sns.set_theme(style="ticks", palette="deep", font_scale=1.1)
def plot_pdp(feature):
    # Calculate Partial Dependence and Individual Conditional Expectation
    pdp = partial_dependence(xgboost_model, X_train, [feature], kind="both", grid_resolution=50)
    # Prepare data for plotting
    plot_x = pd.Series(pdp.grid_values[0]).rename('x')
    plot_y = pd.Series(pdp.average[0]).rename('y')
    plot_i = pdp.individual[0]
    # Smooth interpolation
    tck = splrep(plot_x, plot_y, s=30)
    xnew = np.linspace(plot_x.min(), plot_x.max(), 300)
    ynew = splev(xnew, tck, der=0)
    # Create plot
    fig, ax = plt.subplots(figsize=(8, 6))
    # Plot ICE curves
    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)
    # Plot smoothed PDP
    ax.plot(xnew, ynew, color='peru', linewidth=2, label='Smoothed PDP')
    # Add confidence interval
    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')
    # Add rug plot
    sns.rugplot(data=X_train.sample(100), x=feature, height=0.05, color='k', alpha=0.3, ax=ax)
    # Set labels and limits
    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)
    # Add legend
    ax.legend()
    # Save and show plot
    plt.savefig(f'./pdpplot_{feature}.png', dpi=900, bbox_inches='tight')
    plt.show()
# Plot PDP for each feature
for feature in features_to_plot:
    plot_pdp(feature)

#!/usr/bin/env python # coding: utf-8 
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.patches as patches
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'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 模型
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)
# 需要绘制 PDP 图的特征列表
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'
]
# 设置 Seaborn 主题，使图表更具学术感
sns.set_theme(style="ticks", palette="deep", font_scale=1.1)
def plot_pdp_journal_style(feature):
    """ 绘制期刊级别的局部依赖图 (PDP) """
    # 1. 计算局部依赖 (PDP) 和个体条件期望 (ICE)
    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]
    # 2. 计算动态标准差 (沿着样本轴计算每个网格点的不确定性)
    std_dev = np.std(plot_i, axis=0)
    # 平滑插值 (针对平均 PDP，用于绘制平滑曲线)
    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))
    # 3. 绘制双层动态阴影区间
    # 外层灰色阴影 (±1 倍标准差，表示较大波动范围)
    ax.fill_between(plot_x, plot_y - std_dev, plot_y + std_dev, color='grey', alpha=0.3, linewidth=0)
    # 内层浅黄色阴影 (±0.5 倍标准差，表示较集中波动范围)
    ax.fill_between(plot_x, plot_y - 0.5 * std_dev, plot_y + 0.5 * std_dev, color='khaki', alpha=0.5, linewidth=0)
    # 4. 绘制原始 PDP 均值折线 (深灰色虚线)
    ax.plot(plot_x, plot_y, color='dimgrey', linestyle='--', linewidth=1.5, zorder=4)
    # 5. 绘制平滑后的 PDP 实线 (橙棕色实线)
    ax.plot(xnew, ynew, color='peru', linewidth=2, label='Smoothed PDP', zorder=5)
    # 6. 添加底部地毯图 (Rug Plot)，展示数据分布密度
    sns.rugplot(data=X_train.sample(min(150, len(X_train))), x=feature, height=0.03, color='grey', alpha=0.7, ax=ax)
    # 7. 添加红色虚线强调框 (作为示例，仅针对 DO of Zone 5 添加)
    if feature == 'DO of Zone 5':
        box_x_start = 6.4
        box_width = 2.2
        box_y_start = plot_y.min() - 0.15
        box_height = (plot_y.max() - plot_y.min()) + 0.3
        rect = patches.Rectangle((box_x_start, box_y_start), box_width, box_height, linewidth=2.5, edgecolor='red', facecolor='none', linestyle='--', zorder=6)
        ax.add_patch(rect)
    # 设置标签和坐标轴限制
    ax.set_ylabel('Partial Dependence')
    # 如果特征名包含特定物质，可以自动添加单位，提升图表专业度
    if 'DO' in feature or 'COD' in feature or 'TN' in feature or 'NH4' in feature:
        ax.set_xlabel(f'{feature} (mg/L)')
    else:
        ax.set_xlabel(feature)
    # 让 X 轴两端稍微留白，避免图形顶格
    x_range = plot_x.max() - plot_x.min()
    ax.set_xlim(plot_x.min() - 0.05 * x_range, plot_x.max() + 0.05 * x_range)
    # 添加图例（可选，期刊图中有时不需要图例，如果不需要可以注释掉下面这行）
    # ax.legend(loc='upper right')
    # 保存并展示图表
    plt.savefig(f'./journal_pdp_{feature}.png', dpi=600, bbox_inches='tight')
    plt.show()
    plt.close()
# 遍历所有特征并绘图
for feature in features_to_plot:
    print(f"正在绘制 {feature} 的期刊风格 PDP 图...")
    plot_pdp_journal_style(feature)
print("所有图表绘制完成！")

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

# 1. 读取数据 (请确保路径正确)
file_path = r'Data.csv'
df = pd.read_csv(file_path, encoding='GBK')
# 2. 分离特征和目标变量
X = df.drop(columns='Effluent TN')
y = df['Effluent TN']
# 3. 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# 4. 初始化并训练 XGBoost 模型
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)
# 5. 定义需要绘制 PDP 图的特征列表
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'
]
# 6. 设置 Seaborn 主题，提升学术感
sns.set_theme(style="ticks", palette="deep", font_scale=1.1)
# 7. 定义核心绘图函数
def plot_pdp_journal_style(feature):
    # 计算局部依赖 (PDP) 和个体条件期望 (ICE)
    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]
    # 计算动态标准差 (核心：针对每个网格点计算所有样本的不确定性)
    std_dev = np.std(plot_i, axis=0)
    # 平滑插值 (针对平均 PDP)
    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))
    # 绘制双层动态阴影区间 (外层灰色 ±1 Std, 内层浅黄 ±0.5 Std)
    ax.fill_between(plot_x, plot_y - std_dev, plot_y + std_dev, color='grey', alpha=0.3, linewidth=0)
    ax.fill_between(plot_x, plot_y - 0.5 * std_dev, plot_y + 0.5 * std_dev, color='khaki', alpha=0.5, linewidth=0)
    # 绘制原始均值折线 (灰色虚线) 与平滑曲线 (橙棕色实线)
    ax.plot(plot_x, plot_y, color='dimgrey', linestyle='--', linewidth=1.5, zorder=4)
    ax.plot(xnew, ynew, color='peru', linewidth=2, label='Smoothed PDP', zorder=5)
    # 添加底部地毯图 (限制最多展示 150 个点以防重叠严重)
    sns.rugplot(data=X_train.sample(min(150, len(X_train))), x=feature, height=0.03, color='grey', alpha=0.7, ax=ax)
    # 添加红色虚线强调框 (以 DO of Zone 5 为例)
    if feature == 'DO of Zone 5':
        box_x_start = 6.4
        box_width = 2.2
        box_y_start = plot_y.min() - 0.15
        box_height = (plot_y.max() - plot_y.min()) + 0.3
        rect = patches.Rectangle((box_x_start, box_y_start), box_width, box_height, linewidth=2.5, edgecolor='red', facecolor='none', linestyle='--', zorder=6)
        ax.add_patch(rect)
    # 设置标签与 X 轴单位智能匹配
    ax.set_ylabel('Partial Dependence')
    if any(keyword in feature for keyword in ['DO', 'COD', 'TN', 'NH4', 'Carbon']):
        ax.set_xlabel(f'{feature} (mg/L)')
    elif 'Flowrate' in feature:
        ax.set_xlabel(f'{feature} (m³/d)')
    elif 'Energy' in feature:
        ax.set_xlabel(f'{feature} (kWh)')
    else:
        ax.set_xlabel(feature)
    # 设置 X 轴留白范围
    x_range = plot_x.max() - plot_x.min()
    ax.set_xlim(plot_x.min() - 0.05 * x_range, plot_x.max() + 0.05 * x_range)
    # 保存、展示并关闭画布
    plt.savefig(f'./journal_pdp_{feature.replace(" ", "_").replace("/", "_")}.png', dpi=600, bbox_inches='tight')
    plt.show()
    plt.close()
# 8. 循环遍历所有特征并执行绘图
print("开始生成特征依赖图...")
for feature in features_to_plot:
    print(f"正在绘制：{feature}")
    plot_pdp_journal_style(feature)
print("所有图表已成功生成并保存至当前目录！")

双变量的 PDP

#!/usr/bin/env python # coding: utf-8 
# In[ ]:
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

# Load dataset
file_path = r'Data.csv'
df = pd.read_csv(file_path, encoding='GBK')
# Split features and target
X = df.drop(columns='Effluent TN')
y = df['Effluent TN']
# Split dataset into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Initialize XGBoost model
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
)
# Train the model
xgboost_model.fit(X_train, y_train)
# Feature selection and plot 2D PDP，Complementing other target features and related features
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
# ==========================================
# 下面是调用示例
# ==========================================
# 1. 准备备选方案数据
# 每行代表一个方案（如 NSGA-II 的不同解）
# 列分别是：[出水 TN, 能耗，外加碳源]
# 注意：基于你的代码逻辑，这里的数值都是越小越好
data = np.array([
    [12.5, 300, 50], # 策略 1：TN 偏高，但能耗和碳源极省
    [10.0, 450, 80], # 策略 2：各项指标都处于中等水平
    [8.5, 600, 120], # 策略 3：TN 极低（水质极好），但极其费电费药
    [11.0, 350, 60]  # 策略 4：相对均衡，略偏向节能
])
# 2. 赋予权重
# 假设我们认为出水达标最重要 (0.5)，其次是控制能耗 (0.3)，最后是碳源用量 (0.2)
weights = np.array([0.5, 0.3, 0.2])
# 3. 调用代码计算得分
scores = topsis(data, weights)
# 4. 打印结果
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.py

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

# ==========================================
# 1. 默认数据与变量范围设定
# ==========================================
# 默认工况数据：代表污水厂当前的实际运行状态（包含不可控的进水参数，和当前的操作参数）
default_data = {
    'Influent NH4+-N': 29.6, # 进水氨氮 (固定参数)
    'Influent TN': 35.9, # 进水总氮 (固定参数)
    'Influent Flowrate': 38808.4, # 进水流量 (固定参数)
    'DO of Zone 5': 7.1, # 5 区溶解氧 (控制变量)
    'TN of Zone 5': 22.2, # 5 区总氮 (状态变量)
    'COD of Zone 6': 55.0, # 6 区 COD (控制变量)
    'TN of Zone 6': 13.8, # 6 区总氮 (状态变量)
    'DO of Zone 7': 4.06, # 7 区溶解氧 (控制变量)
    'COD of Zone 7': 31.9, # 7 区 COD (控制变量)
    'TN of Zone 7': 11.7, # 7 区总氮 (状态变量)
    'External Carbon Source': 3.3 # 外加碳源 (控制变量)
}
# 优化变量边界：定义遗传算法在搜索最佳操作参数时，不能超过的工艺极限
variable_ranges = [
    (5.8, 7.5), # DO of Zone 5 的寻优范围
    (28.0, 35.5), # COD of Zone 6 的寻优范围
    (0.3, 2.8), # DO of Zone 7 的寻优范围
    (1.0, 26.4), # COD of Zone 7 的寻优范围
    (1.9, 2.3) # External Carbon Source 的寻优范围
]
# 初始化标准化工具（为模型预测前的数据预处理做准备）
df = pd.DataFrame([default_data])
scaler = StandardScaler()
scaler.fit(df)
# 初始化加载两个 XGBoost 预测模型（一个测能耗，一个测出水 TN）
model_tn = xgb.Booster()
model_power = xgb.Booster()
# 预测缓存字典：记忆已经算过的参数组合，极大地加速遗传算法
prediction_cache = {}

# ==========================================
# 2. 核心评估函数定义
# ==========================================
def predict_outputs(data):
    """使用 XGBoost 预测给定参数下的能耗和 TN，并利用缓存机制加速"""
    # 将一维输入转为二维，并生成用于缓存的特征键值 (tuple)
    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]])
    # 如果没预测过，调用 XGBoost 进行预测
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):
    """
    适应度函数：计算当前虚拟个体 (individual) 与原始实际工况的差值
    差值越小（甚至为负），说明该方案比当前实际运行策略更好
    """
    # 拼接进水参数 (固定) + 当前尝试的操作参数 (可变) -> 完整的模型输入
    full_input = np.concatenate([fixed_water_params, individual])
    energy, tn = predict_outputs(np.array(full_input))
    # 计算当前方案预测值与基准值的差（Delta）
    delta_energy = energy[0] - true_energy
    delta_tn = tn[0] - true_tn
    return delta_energy, delta_tn

# ==========================================
# 3. 遗传算法优化主流程
# ==========================================
def optimize_parameters(fixed_water_params, initial_operating_params):
    # --- A. DEAP 算法框架初始化 ---
    # 定义优化目标：权重 (-1.0, -1.0) 代表这是一个双目标【最小化】问题
    if not hasattr(creator, 'FitnessMin'):
        creator.create("FitnessMin", base.Fitness, weights=(-1.0, -1.0))
    # 定义个体：继承自 list，包含上述的适应度属性
    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)
    # 注册交叉算子：模拟二进制交叉 (SBX)，带边界约束
    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]))
            # 截断越界值
            # 修复 NaN/Inf 等会导致模型崩溃的异常值
            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)
    # 注册选择算子：采用经典的 NSGA-II 拥挤度距离排序法
    toolbox.register("select", tools.selNSGA2)
    # --- B. 计算基准状态 ---
    # 评估当前实际工况的表现，作为优化的参考标尺
    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)
    # --- C. 初始化种群与帕累托前沿 ---
    pop = toolbox.population(n=100) # 生成 100 个初始个体
    pop[0][:] = initial_operating_params.copy() # 【启发式策略】将当前实际工况混入初始种群，确保底线
    # 计算初始种群适应度
    for ind in pop:
        ind.fitness.values = toolbox.evaluate(ind)
    # 名人堂（ParetoFront）：用于一直保存进化过程中找到的最优不被支配解
    hof = tools.ParetoFront()
    hof.update(pop)
    # 早停机制变量初始化
    best_fitness = min([ind.fitness.values[0] for ind in pop])
    stall_count = 0
    # --- D. 进化循环 (核心过程) ---
    for gen in range(20): # 最多繁衍 20 代
        # 产生子代（交叉概率 50%，变异概率 50%）
        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) + 子代 (offspring) 中挑选最优秀的 20 个个体进入下一代
        pop = toolbox.select(pop + offspring, k=20)
        hof.update(pop) # 更新全局最优解集
        # --- 早停判定逻辑 ---
        # 如果模型已经连续 5 代没有找到更好的解，提前结束以节省时间
        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
    # --- E. 结果输出 ---
    # 为了避免输出上百个结果让用户无法挑选，均匀采样提取最多 60 个代表性方案
    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]

模型可视化

#!/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

# 1. 默认数据字典：定义了进水参数、各生化池状态以及外部碳源的基准值
# 当用户在 GUI 中没有输入任何值时，系统会默认使用这些数值进行预测
default_data = {
    'Influent NH4+-N': 29.6, # 进水氨氮
    'Influent TN': 35.9, # 进水总氮
    'Influent Flowrate': 38808.4, # 进水流量
    'DO of Zone 5': 7.1, # 5 区溶解氧
    'TN of Zone 5': 22.2, # 5 区总氮
    'COD of Zone 6': 55.0, # 6 区化学需氧量
    'TN of Zone 6': 13.8, # 6 区总氮
    'DO of Zone 7': 4.06, # 7 区溶解氧
    'COD of Zone 7': 31.9, # 7 区化学需氧量
    'TN of Zone 7': 11.7, # 7 区总氮
    'External Carbon Source': 3.3 # 外加碳源量
}
# 将默认数据转换为 Pandas DataFrame 格式，方便后续处理
df = pd.DataFrame([default_data])
# 2. 数据标准化处理
scaler = StandardScaler()
scaler.fit(df)
# 3. 加载预训练的 XGBoost 模型
# 加载预测出水总氮 (Effluent TN) 的模型
model_tn = xgb.Booster()
model_tn.load_model('model_effluent_TN.json')
# 加载预测能耗 (Energy Consumption) 的模型
model_power = xgb.Booster()
model_power.load_model('model_energy_consumption.json')
# 4. 辅助函数：验证用户输入是否为合法的浮点数
def validate_input(value):
    try:
        return float(value)
    except ValueError:
        return None
# 如果输入包含字母或特殊符号导致转换失败，返回 None
# 5. 全局控制变量
run_count = 0 # 记录点击'Predict'按钮的次数（0 代表首次原始预测，1 代表二次优化预测）
first_external_carbon_source = None # 记录首次预测时的外加碳源量，用于后续计算差值
# 6. 重置功能：清空所有输入框、结果展示框，并重置全局变量
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)
    # 清空下方结果显示区的 StringVar 变量
    original_tn.set("")
    original_power.set("")
    optimized_tn.set("")
    optimized_power.set("")
    energy_savings_var.set("")
    carbon_savings_var.set("")
    output_text.set("")
# 7. 核心预测逻辑
def predict():
    global run_count, first_external_carbon_source
    # 限制预测次数最多为 2 次（原始状态 vs 优化状态）
    if run_count >= 2:
        output_text.set("The prediction limit has been reached. Please click the Reset button to reset and try again.\n(预测次数已达上限，请点击 Reset 按钮重置后再试。)")
        return
    # 收集用户输入
    inputs = []
    for feature, entry in feature_entry_map.items():
        value = entry.get().strip() # 获取输入框内容并去除首尾空格
        if not value: # 如果输入为空，则采用上方定义的 default_data 中的默认值
            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:
        # 将收集到的输入列表转换为 DataFrame，保持特征列名一致
        inputs_df = pd.DataFrame([inputs], columns=default_data.keys())
        # 使用 StandardScaler 对输入数据进行特征缩放（标准化）
        inputs_scaled = scaler.transform(inputs_df)
        # 转换为 XGBoost 专用的数据结构 DMatrix
        dmatrix = xgb.DMatrix(inputs_scaled)
        # 调用模型进行预测，[0] 取出预测数组中的第一个标量值
        effluent_tn_pred = model_tn.predict(dmatrix)[0]
        power_consumption_pred = model_power.predict(dmatrix)[0]
        # ---------------- 逻辑分支 ----------------
        # # 第一次点击 Predict：记录并显示'原始'运行状态的结果
        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]
        # 第二次点击 Predict：记录并显示'优化后'的运行结果，并计算节能和节碳量
        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}")
# 8. GUI 界面中各项参数的参考范围提示文本
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"
}
# ==========================================
# 9. 构建 Tkinter GUI 界面
# ==========================================
root = tk.Tk()
root.title("Effluent TN Prediction and Optimization")
root.geometry("1200x800")
root.configure(bg="#e8f0f2") # 设置主窗口背景色为浅灰蓝色
default_font = ("Arial", 12) # 创建主容器 Frame
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")
]
# 用于保存界面上所有 Entry(输入框) 的字典，键为特征名，值为对应的控件对象
feature_entry_map = {}
# 渲染左列 (Water Quality Parameters)
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")
# 渲染右列 (Controllable Parameters)
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")
# 10. 底部结果展示区变量绑定
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 函数
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)
# 重置按钮：绑定 clear_inputs 函数
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 = tk.Frame(main_frame, bg="#e8f0f2")
output_frame.grid(row=9, column=0, columnspan=8, padx=10, pady=20)
# 第一行：总氮 (TN) 结果对比
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)
# readonly 防止用户修改预测结果
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)
# 注意：这里没有设为 readonly，如果需要防止修改建议加上 state="readonly"
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()