如何使用Python求解微分方程组优化间歇反应器目标参数

核心思路 #

• 首先将原有的微分方程求解逻辑封装为以[B0, HA, T]为输入的可调用函数，函数返回负的全反应周期内X+A的最大值（Python优化库默认求解最小值，要最大化目标只需对结果取负即可）
• 选择带边界约束的全局优化算法完成参数寻优，避免陷入局部最优，使用scipy自带的差分进化算法differential_evolution即可满足需求

完整可运行代码 #

import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import differential_evolution

# 固定公共参数（不需要参与优化的部分）
xc = 0.4068       # Cellulose
xh1 = 0.22136     # Xylans
xh2 = 0.03786     # Arabinans
xh3 = 0.0333      # Acetyls
R = 8.314/1000 # kJ/(mol K)
k10 = 2.37 # Xylans and arabinans preexponential factors
n1 = 1.51  # Xylans and arabinans acid concentration exponent
E1 = 83.3 # Xylans and arabinans activation energy (kJ/mol)
k20 = 2.17 # Xylose preexponential factor
n2 = 0.29 # Xylose acid concentration exponent
E2 = 143.5 # Xylose activation energy (kJ/mol)
k1ac0 = 2.37 # Acetyls preexponential factor
n1ac = 0.604 # Acetyls acid concentration exponent
E1ac = 83.3 # Acetyls activation energy (kJ/mol)
k30 = 2.37 # Cellulose preexponential factor
n3 = 1.359 # Cellulose acid concentration exponent
E3 = 94.962 # Cellulose activation energy (kJ/mol)
Bmax = 120 # Maximum Biomass Concentration (g/L)
Xylmax = xh1*Bmax  # Maximum Xylan Concentration (g/L)
Arabmax = xh2*Bmax # Maximum Arabinan Concentration (g/L)
Acetmax = xh3*Bmax # Maximum Acetyl Concentration (g/L)
Celmax = xc*Bmax   # Maximum Cellulose Concentration (g/L)
t = np.linspace(0, 240, 720) # Time (min)

# 重构动力学模型，支持传入HA和T参数
def model(z, t, HA, T):
    Xyl = z[0] # Xylan concentration (g/L)
    X = z[1] # Xylose concentration (g/L)
    Arab = z[2] # Arabinan concentration (g/L)
    A = z[3] # Arabinose concentration (g/L)
    Acet = z[4] # Acetyl concentration (g/L)
    w = z[5] # derivative of Acet
    Ac = z[6] # Acetic acid concentration (g/L)
    u = z[7] # derivative of acetic acid
    Cel = z[8] # Cellulose concentration (g/L)
    G = z[9] # Glucose concentration (g/L)
    F = z[10] # Furfural concentration (g/L)
    Xyl0 = xh1 * 25 # 取参数上限计算避免除零，不影响动力学结果
    K1 = k10*(10**10)*(HA**n1)*np.exp(-E1/(R*T)) # min-1
    K2 = k20*(10**15)*(HA**n2)*np.exp(-E2/(R*T)) # min-1
    K1ac = k1ac0*(10**10)*(HA**n1ac)*((Xyl/Xyl0)**2)*np.exp(-E1ac/(R*T)) # min-1
    K3 = k30*(10**10)*(HA**n3)*np.exp(-E3/(R*T)) # min-1
    ef1 = 1/(1+((Xyl/Xylmax)**20)) # Xylan effective coefficient
    ef2 = 1 / (1 + ((Arab / Arabmax)**20)) # Arabinan effective coefficient
    ef3 = 1 / (1 + ((Acet/Acetmax)**20)) # Acetyl effective coefficient
    ef4 = 1 / (1 + ((Cel/Celmax)**20)) # Cellulose effective coefficient
    Xylef = 0.95*ef1*Xyl+(1-ef1)*Xylmax # Xylan effective concentration
    Arabef = 0.95*ef2*Arab+(1-ef2)*Arabmax # Arabinan effective coefficient
    Acetef = 0.95*ef3*Acet+(1-ef3)*Acetmax # Acetyl effective coefficient
    Celef = 0.95*ef4*Cel+(1-ef4)*Celmax # Cellulose effective coefficient
    dXyldt = -K1*Xylef
    dXdt = K1*Xylef-K2*X
    dArabdt = -K1*Arabef
    dAdt = K1*Arabef-K2*A
    dAcetdt = w
    dwdt = (-K1ac*Acetef-10*w)/4
    dAcdt = u
    dudt = (K1ac*Acetef-10*u)/4
    dCeldt = -K3*Celef
    dGdt = K3*Celef
    dFdt = K2*(X+A)
    return [dXyldt, dXdt, dArabdt, dAdt, dAcetdt, dwdt, dAcdt, dudt, dCeldt, dGdt, dFdt]

# 目标函数：输入参数组，返回负的最大X+A
def objective(params):
    B0, HA, T = params
    # 计算初始条件
    Xyl0 = xh1*B0
    Arab0 = xh2*B0
    Acet0 = xh3*B0
    w0 = 0
    Cel0 = xc*B0
    X0 = 0
    A0 = 0
    Ac0 = 0
    u0 = 0
    G0 = 0
    F0 = 0
    z0 = [Xyl0, X0, Arab0, A0, Acet0, w0, Ac0, u0, Cel0, G0, F0]
    # 求解微分方程
    z = odeint(model, z0, t, args=(HA, T))
    X = z[:, 1]
    A = z[:, 3]
    max_total = np.max(X + A)
    return -max_total

# 定义参数边界
bounds = [(15, 25), (0.03, 0.22), (373, 453)]

# 调用差分进化算法寻优
result = differential_evolution(objective, bounds, seed=42, popsize=15, maxiter=100)

# 输出结果
print("===== 寻优结果 =====")
print(f"最优初始生物质浓度 B0 = {result.x[0]:.2f} g/L")
print(f"最优硫酸浓度 HA = {result.x[1]:.3f} mol/L")
print(f"最优反应温度 T = {result.x[2]:.1f} K（即 {result.x[2]-273.15:.1f} ℃）")
print(f"最大X+A浓度 = {-result.fun:.2f} g/L")

可选调优说明 #

• 若需要更快的求解速度，可以先降低popsize和maxiter得到初步最优解，再调用scipy.optimize.minimize做局部精修
• 若对X+A的取值时间有约束（比如要求反应120min内达到最大值），可以修改目标函数中取最大值的时间切片范围
• 优化完成后可将最优参数代入原代码，即可输出对应的浓度、收率曲线做验证

内容的提问来源于stack exchange，提问作者user17345875