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

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

模型开发与评估

在模型开发阶段，预处理后的数据集按 8:2 的比例划分为训练集（80%）和测试集（20%）。为了提升模型的泛化能力，我们首先通过 Pearson 相关性分析剔除预测增益低且存在共线性的特征，同时引入 SHAP 值系统评估特征的重要性。

模型对比与选择

研究对比了多元线性回归（MLR）、支持向量机（SVM）、轻量级梯度提升机（LightGBM）、随机森林（RF）、极端梯度提升（XGBoost）和类别特征提升（CatBoost）六种机器学习算法。

结果显示，MLR 的预测能力有限（测试集 R² = 0.459），而其余五种非线性模型的性能表现优异。其中，XGBoost 模型表现最为突出，训练集准确率接近完美（R² = 0.999），但在测试集上的泛化能力略弱（R² = 0.829）。经过特征修剪（剔除冗余特征如进水 pH、废水温度、出水 COD）后，优化后的 XGBoost 模型性能显著提升，训练集和测试集的 R²分别达到 0.997 和 0.911，RMSE 分别为 0.196 和 1.283。

图 2f 显示了 XGBoost 模型的拟合曲线及其 95% 预测区间。由于绝大多数测试数据点均落在 95% 置信区间内，充分证明了该模型的可靠性与鲁棒性。因此，后续研究将采用 XGBoost 模型进行出水总氮浓度的计算与分析。

图 2. 以出水总氮浓度为目标特征的六种机器学习模型训练图。(f) 极端梯度提升 (XGBoost)。

代码实现：模型训练与评估

下面展示了核心的模型训练流程，包括数据加载、标准化、网格搜索超参数优化以及交叉验证策略。这里特别加入了提前停止（Early Stopping）机制，防止过拟合——若验证集均方误差连续 50 轮未改善则停止训练。

#!/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 与 PDP 分析

利用 SHAP 摘要图和部分依赖图（PDP），我们可以深入理解各特征的影响力及其相关性。

1. 关键特征识别：TN_7、EC（能耗）、DO_5、TN_6、进水流量以及 DO_7 是决定出水总氮浓度的最重要特征。其中，DO_5 和 DO_7 与能耗（EC）显著相关，这为后续制定智能节能策略提供了逻辑基础。
2. 单特征影响：维持较高的 DO_5 浓度有利于提高总氮去除效率，但该工艺中 5 号区的曝气兼具搅拌功能，导致其 DO_5 浓度时常超过 10 mg/L，这凸显了利用机器学习优化 DO_5 以降低能耗的必要性。
3. 交互作用：较低的 DO_7 与较高的 DO_5 结合可实现更高的脱氮效率。适量的 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。

多目标优化与 GUI 开发

为了便于实际应用，我们开发了图形用户界面（GUI），集成了 NSGA-II（带策略非支配排序遗传算法）和 TOPSIS（逼近理想解排序法）算法。用户可以通过输入进水水质数据和核心工艺参数，自动预测并优化出水总氮、能耗和外加碳源的使用。

核心优化代码片段

以下是 NSGA-II 优化主流程和 TOPSIS 决策方法的实现逻辑。

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,
    '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),   # 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)

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

# ==========================================
# 2. 核心评估函数定义
# ==========================================
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

# ==========================================
# 3. 遗传算法优化主流程
# ==========================================
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 交互界面

为了方便现场工程师使用，我们基于 Tkinter 开发了图形用户界面。它允许用户输入当前工况参数，系统会自动预测原始出水总氮和能耗，并根据优化算法推荐新的参数组合，从而同步降低出水总氮、能耗和外加碳源用量。

#!/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. 默认数据字典
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')

run_count = 0
first_external_carbon_source = None
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()
feature_entry_map = {}
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"
}

def validate_input(value):
    try:
        return float(value)
    except ValueError:
        return 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}")

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")]

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")

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()

实施效果

通过优化关键参数，出水总氮浓度从 14.29 mg/L 降至 11.59 mg/L，污水处理厂能耗每天减少了 128.66 kWh，COD 投加浓度降低了 34.29 mg/L，对应外加碳源投加量从 3.72 m³/d 降至 2.25 m³/d。

全年数据显示，出水总氮浓度平均下降了 17.50%，而 COD 投药量减少了 33.29%。值得注意的是，仅通过强化总氮去除、降低能耗与碳源投加量，即可实现年碳减排量 788.40 吨二氧化碳。

本研究成果为进水水质波动压力下、采用生态组合塘工艺的城镇污水处理厂提供了总氮强化去除的优化模型，助力其既满足严格的总氮排放标准，又助力节能与碳减排目标的实现。