Python 实现活性污泥模型 ASM2d 参数校准与优化

Python 实现活性污泥模型 ASM2d 参数校准与优化

活性污泥模型（ASM2d）是污水处理工艺模拟的核心工具。在实际工程中，模型参数的准确性直接决定了仿真结果的可信度。本文基于 QSDSAN 库与 Optuna 框架，详细解析如何从零开始实现 ASM2d 的参数校准流程。我们将通过 Python 代码复现论文中的参数调优方法，涵盖进水浓度修正、灵敏度分析筛选、单目标优化及动态仿真等关键环节。

1. 进水浓度调整

在运行仿真前，必须确保输入的水质数据与模型内部计算一致。由于 ASM2d 组分复杂，直接输入的 COD/TN/TP 往往无法直接映射到具体的组分浓度。我们需要通过迭代调整 WasteStream 的浓度输入，使得计算得到的实际值与目标值在小数位上保持一致。

核心逻辑如下：

1. 设定初始调整值为目标值。
2. 构建 WasteStream 并计算实际 COD/TN/TP。
3. 比较计算值与目标值的差异。
4. 若未收敛，则根据误差更新调整值，继续循环。

import qsdsan as qs
from qsdsan import WasteStream

def Adjusting_influent_concentration(first_inflow_value, input_COD_value, input_TN_value, input_TP_value):
    """
    目的：根据用户输入的目标进水 COD/TN/TP 值，迭代调整 WasteStream 的浓度输入，
    使得通过 WasteStream 计算得到的实际 COD/TN/TP 与目标值一致。
    """
    adjusting_COD_value = input_COD_value
    adjusting_TN_value = input_TN_value
    adjusting_TP_value = input_TP_value

    def simulate_influent_internal(first_inflow_value, adjusting_COD_value, adjusting_TN_value, adjusting_TP_value):
        Q_inf = first_inflow_value
        Temp = 273.15 + 20
        influent = WasteStream('influent', T=Temp)
        default_inf_kwargs = {
            'concentrations': {
                'S_I': 0, 'X_I': 0, 'S_F': 0, 'S_A': 0,
                'X_S': adjusting_COD_value,
                'S_NH4': 0, 'S_N2': 0, 'S_NO3': adjusting_TN_value,
                'S_PO4': adjusting_TP_value,
                'X_PP': 0, 'X_PHA': 0, 'X_H': 0.15,
                'X_AUT': 0, 'X_PAO': 0, 'S_ALK': 7 * 12,
            },
            'units': ('m3/d', 'mg/L'),
        }
        influent.set_flow_by_concentration(Q_inf, **default_inf_kwargs)
        return round(influent.COD, 2), round(influent.TN, 2), round(influent.TP, 2)

    iteration_count = 0
    while True:
        iteration_count += 1
        simulate_influent_COD, simulate_influent_TN, simulate_influent_TP = simulate_influent_internal(
            first_inflow_value, adjusting_COD_value, adjusting_TN_value, adjusting_TP_value)

        if (round(simulate_influent_COD, 3) == round(input_COD_value, 3) and
            round(simulate_influent_TN, 3) == round(input_TN_value, 3) and
            round(simulate_influent_TP, 3) == round(input_TP_value, 3)):
            break

        difference_COD = input_COD_value - simulate_influent_COD
        adjusting_COD_value += difference_COD
        difference_TN = input_TN_value - simulate_influent_TN
        adjusting_TN_value += difference_TN
        difference_TP = input_TP_value - simulate_influent_TP
        adjusting_TP_value += difference_TP

        if iteration_count > 100:
            print(f"警告：迭代次数超过 100 次，停止调整")
            break

    return adjusting_COD_value, adjusting_TN_value, adjusting_TP_value

2. 参数选择（OSA 方法）

为了减少优化维度，我们采用 OSA（基于统计/灵敏度文件）方法进行参数筛选。该步骤读取预先计算的敏感性分析文件，选取对目标指标（如 COD）影响最大的前 7 个参数作为优化变量。

这不仅能加快收敛速度，还能避免无关参数干扰优化结果。代码主要涉及读取 Excel 文件，提取参数名及平均值，并与原始默认参数字典进行匹配。

import pandas as pd

def OSA(name_1):
    """
    OSA: 读取预先计算的参数灵敏度文件，选取灵敏度最高（平均值最大）的前 7 个参数，
    并返回这些参数及其原始默认值作为优化或调参的输入。
    """
    file_path = fr'/root/autodl-tmp/data/raw/{name_1}_OSA_FP.xlsx'
    df = pd.read_excel(file_path)
    top_parameters = df[['Parameter', 'average_value']].sort_values(by='average_value', ascending=False).head(7)

    original_params = {
        'f_SI': 0, 'Y_H': 0.625, 'f_XI_H': 0.1, 'Y_PAO': 0.625, 'Y_PO4': 0.4,
        # ... 更多默认参数 ...
        'mu_AUT': 1, 'b_AUT': 0.15, 'K_O2_AUT': 0.5
    }

    params = []
    for param in top_parameters['Parameter']:
        if param in original_params:
            params.append({'name': param, 'initial_value': original_params[param]})

    param_names = [p['name'] for p in params]
    param_initial_values = [p['initial_value'] for p in params]
    return param_names, param_initial_values

3. 单目标优化

选定参数后，利用 Optuna 进行单目标优化。我们以最小化出水 COD 或 TN 的百分比误差为目标函数。

在 objective 函数中，Optuna 会采样一组参数，调用仿真函数计算误差。如果仿真失败（如出现 NaN），则剪枝该试验。优化过程中，系统会自动记录参数重要性（灵敏度），并生成可视化图表（平行坐标图、优化历史图），帮助我们理解参数间的关联。

import optuna

def OT(param_names, param_initial_values, day, name_1):
    """
    使用 Optuna 对单一目标（COD 或 TN）进行参数优化。
    """
    # 定义目标函数
    def objective(trial):
        parameter_adjustments = []
        for i in range(len(param_names)):
            adjustment = trial.suggest_uniform(
                param_names[i], 
                param_initial_values[i] * 0.5, 
                param_initial_values[i] * 1.5
            )
            parameter_adjustments.append(adjustment)

        try:
            # 调用仿真函数获取误差
            cod_err, tn_err, tp_err = simulate_one(..., parameter_adjustments, ...)
            
            if np.isnan(cod_err) or np.isinf(cod_err):
                raise optuna.exceptions.TrialPruned()

            if name_1 == "COD":
                trial.set_user_attr('TN_percentage_difference', tn_err)
                return cod_err
            elif name_1 == "TN":
                trial.set_user_attr('COD_percentage_difference', cod_err)
                return tn_err
        except Exception as e:
            raise optuna.exceptions.TrialPruned()

    # 创建 Study 对象
    sqlite_path = fr"../data/process/{day}_{name}.db"
    storage_url = f"sqlite:///{sqlite_path.replace(os.sep, '/')}"
    study = optuna.create_study(direction="minimize", storage=storage_url, load_if_exists=True)

    # 执行优化
    study.optimize(objective, n_trials=70, n_jobs=-1)

    # 保存结果与可视化
    # ... 保存 HTML 图表与最优参数 ...
    return study

4. 仿真函数核心

仿真是整个流程的引擎。simulate_one 函数负责搭建反应器系统（厌氧池、好氧池、沉淀池），设置生物过程模型（ASM2d），并执行动态求解。

关键点在于：

• 正确配置反应器容积与回流比。
• 初始化各单元浓度（通常从 Excel 加载）。
• 设置求解器容差（BDF 方法）。
• 提取出水数据并计算误差。

def simulate_one(first_inflow_value, external_return_flow_value, input_COD_value, output_COD_value,
                 adjusting_COD_value, input_TN_value, output_TN_value, adjusting_TN_value,
                 input_TP_value, output_TP_value, adjusting_TP_value, param_names,
                 param_initial_values, parameter_adjustments, day, name, trial):
    # 1. 设置流量与温度
    Q_inf = first_inflow_value
    Q_was = 385
    Q_ext = external_return_flow_value
    Temp = 273.15 + 20

    # 2. 创建 WasteStream 对象
    influent = WasteStream('influent', T=Temp)
    effluent = WasteStream('effluent', T=Temp)
    # ... 其他流股 ...

    # 3. 设置进水组成
    default_inf_kwargs = { ... } # 同上文
    influent.set_flow_by_concentration(Q_inf, **default_inf_kwargs)

    # 4. 配置反应器系统
    V_an = 1000 * 33.33
    V_ae = 1333 * 33.33
    aer1 = pc.DiffusedAeration('aer1', DO_ID='S_O2', KLa=240, DOsat=8.0, V=V_ae)
    # ... 构建 CSTR 单元 A1, A2, O1, O2, O3, C1 ...

    # 5. 设置生物过程模型
    asm2d = pc.ASM2d()
    p2 = asm2d.aero_hydrolysis
    adjusted_params = {}
    for i in range(len(param_names)):
        adjusted_params[param_names[i]] = parameter_adjustments[i]
    p2.set_parameters(**adjusted_params)

    # 6. 构建系统并初始化
    sys = System('example_system', path=(...), recycle=(...))
    # ... 加载初始数据 ...

    # 7. 运行仿真
    sys.simulate(t_span=(0, 50), method='BDF')

    # 8. 计算误差
    simulate_COD = effluent.COD
    simulate_TN = effluent.TN
    simulate_TP = effluent.TP

    COD_pct_diff = abs((output_COD_value - simulate_COD) / output_COD_value * 100)
    TN_pct_diff = abs((output_TN_value - simulate_TN) / output_TN_value * 100)
    TP_pct_diff = abs((output_TP_value - simulate_TP) / output_TP_value * 100)

    return COD_pct_diff, TN_pct_diff, TP_pct_diff

5. 主流程整合

最后，将上述模块整合到主循环中。程序逐日读取进水与出水观测数据，自动执行参数选择与调优，并将最优结果保存至本地。通过这种方式，我们可以实现对污水处理厂长期运行的自动化参数校准。

def main():
    # 读取输入数据
    file_path = r'..\data\raw\input.xlsx'
    df = pd.read_excel(file_path)

    # 主循环：逐日执行
    for day in range(1, 50):
        # 提取当日数据
        first_inflow_value = df['inflow_flow'].iloc[day]
        # ... 其他参数 ...

        # 配置优化目标
        name_1 = 'TN'  # 优化目标
        name_2 = 'FP'  # 参数选择方法
        name_3 = 'TT'  # 调优方法
        name = f'{name_1}_{name_2}_{name_3}'

        # 动态加载函数
        param_function = globals().get(name_2)
        tuning_function = globals().get(name_3)

        # 执行流程
        param_names, param_initial_values = param_function(name_1)
        tuning_function(param_names, param_initial_values, day, name_1)

if __name__ == "__main__":
    main()

通过以上步骤，我们完成了一个完整的 ASM2d 参数校准框架。实际应用中，请根据具体项目路径修改文件引用，并根据水质特性调整优化策略。