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

2. 参数选择 (OSA)

为了减少优化维度，我们采用 OSA (Optimal Sensitivity Analysis) 方法。通过读取预先计算的灵敏度文件，筛选出对特定目标（如 COD 或 TN）影响最大的前 7 个参数作为优化变量。

代码实现

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

    # 原始参数字典（备选）
    original_params = {
        'f_SI': 0, 'Y_H': 0.625, 'f_XI_H': 0.1, 'Y_PAO': 0.625, 'Y_PO4': 0.4,
        'Y_PHA': 0.2, 'f_XI_PAO': 0.1, 'Y_A': 0.24, 'f_XI_AUT': 0.1, 'K_h': 3,
        # ... 更多参数定义 ...
        'mu_AUT': 1, 'b_AUT': 0.15, 'K_O2_AUT': 0.5, 'K_NH4_AUT': 1
    }

    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)

选定参数后，使用 Optuna 进行自动化调优。这里以最小化 COD 或 TN 的模拟误差为目标。

目标函数与执行

def OT(param_names, param_initial_values, day, name_1):
    """
    使用 Optuna 对单一目标（COD 或 TN）进行参数优化。
    """
    # 调整进水浓度
    adjusting_COD_value, adjusting_TN_value, adjusting_TP_value = Adjusting_influent_concentration(
        first_inflow_value, input_COD_value, input_TN_value, input_TP_value)

    def objective(trial):
        parameter_adjustments = []
        # 在 [0.5*init, 1.5*init] 范围内随机采样
        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_diff, TN_diff, TP_diff = 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)

            # 异常剪枝
            if np.isnan(COD_diff) or np.isinf(COD_diff):
                raise optuna.exceptions.TrialPruned()

            # 根据目标返回误差
            if name_1 == "COD":
                trial.set_user_attr('TN_percentage_difference', TN_diff)
                return COD_diff
            elif name_1 == "TN":
                trial.set_user_attr('COD_percentage_difference', COD_diff)
                return TN_diff
        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", study_name="parameter_optimization",
                                storage=storage_url, load_if_exists=True)

    # 执行优化
    study.optimize(objective, n_trials=70, n_jobs=-1,
                   catch=(RuntimeError, LookupError, ValueError, AssertionError))

    # 记录结果与可视化
    importance_df = pd.DataFrame(list(optuna.importance.get_param_importances(study).items()),
                                 columns=['Parameter', 'Importance'])
    # ... 保存 Excel 及 HTML 图表逻辑 ...
    return study

4. 仿真函数 (simulate_one)

这是连接模型与优化器的核心桥梁。它负责构建完整的污水处理系统（厌氧池、好氧池、沉淀池等），运行动态仿真，并计算模拟出水与观测值的误差。

关键步骤

1. 系统构建：使用 qsdsan 库定义反应器单元（CSTR）和生物化学过程（ASM2d）。
2. 参数注入：将 Optuna 生成的参数调整值应用到具体的生化过程中。
3. 动态求解：设置时间步长和求解器（BDF 方法），运行仿真。
4. 误差计算：提取出水浓度，计算百分比差异。

def simulate_one(..., param_names, param_initial_values, parameter_adjustments, ...):
    # 1. 基础设置
    Q_inf = first_inflow_value
    Temp = 273.15 + 20
    influent = WasteStream('influent', T=Temp)
    # ... 设置其他流股 ...

    # 2. 配置曝气系统和反应器
    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)
    
    # 3. 设置生物过程模型
    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)

    # 4. 构建反应器系统
    A1 = su.CSTR('A1', ins=[influent, int_recycle, ext_recycle], V_max=V_an, suspended_growth_model=asm2d)
    # ... 构建 A2, O1, O2, O3, C1 ...
    sys = System('example_system', path=(A1, A2, O1, O2, O3, C1), recycle=(int_recycle, ext_recycle))

    # 5. 初始化与仿真
    sys.set_dynamic_tracker(influent, effluent, A1, A2, O1, O2, O3, C1)
    sys.simulate(t_span=(0, 50), method='BDF')

    # 6. 提取结果与误差
    simulate_COD = effluent.COD
    COD_percentage_difference = abs((output_COD_value - simulate_COD) / output_COD_value * 100)
    return COD_percentage_difference, TN_percentage_difference, TP_percentage_difference

5. 主流程整合

最后，将所有模块串联起来。通过修改配置项（name_1, name_2, name_3），可以灵活切换优化目标和算法策略。

# 主循环示例
for day in range(1, 50):
    # 读取当天数据
    first_inflow_value = inflow_data.iloc[day]
    # ... 获取其他观测值 ...

    # 配置策略
    name_1 = 'TN'       # 优化目标
    name_2 = 'FP'       # 参数选择方法 (OSA/FP/TSA)
    name_3 = 'TT'       # 参数调优方法 (OT/MO_OT/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)

通过上述流程，我们可以实现对 ASM2d 模型的自动化参数校准。在实际应用中，建议根据具体水质数据和工艺特点，适当调整灵敏度分析阈值和优化试验次数，以获得更稳健的模型参数。