如何用Scipy minimize优化ODE参数？解决ValueError报错

问题背景 #

定义了Langmuir吸附的ODE函数：

def langmuir_ode(t, AL, Rmax, kon, koff):
    A = np.interp(t, time, concs)
    L = Rmax - AL
    return kon * A * L - koff * AL

其中time和concs为已知的时间、浓度数组，Rmax、kon、koff是需要优化的参数。

目标函数用于计算模拟数据ysim与ODE求解结果的误差：

def obj(t, *args):
    sol = solve_ivp(fun=langmuir_ode, 
                    t_span=[0,t[-1]], 
                    y0=[0.0], 
                    t_eval=t, 
                    args=args
                    )
    sol = sol.y.ravel()
    error = ysim - sol
    res_sd = np.std(error)
    print('Standard deviation of fit residuals:', format(res_sd,'.2f'))
    SSQ = np.sum(error**2)
    return SSQ

用真实参数调用obj(time, 250, 1e5, 5e-3)能得到合理结果，但调用minimize时触发报错：

initial_values = 300, 5e4, 4e-5
fit = minimize(obj, initial_values)

报错信息：

ValueError: Values in `t_eval` are not within `t_span`.

错误原因 #

scipy.optimize.minimize要求目标函数第一个参数必须是待优化的参数数组，但你当前的obj函数第一个参数是固定的时间数组t。调用minimize时，程序会把initial_values作为第一个参数传给obj，也就是把[300,5e4,4e-5]当成t传入，导致t_span=[0,4e-5]，而t_eval是[300,5e4,4e-5]，其中300、5e4远大于4e-5，触发t_eval不在t_span范围内的报错。

修正方案 #

1. 调整目标函数参数顺序 #

将待优化参数放在obj的第一个位置，固定参数（t、ysim、time、concs）通过minimize的args参数传入：

def obj(params, t, ysim, time, concs):
    Rmax, kon, koff = params
    # 封装符合solve_ivp要求的ODE函数
    def ode_fun(t, AL):
        A = np.interp(t, time, concs)
        L = Rmax - AL
        return kon * A * L - koff * AL
    
    sol = solve_ivp(fun=ode_fun, 
                    t_span=[0, t[-1]], 
                    y0=[0.0], 
                    t_eval=t)
    sol = sol.y.ravel()
    error = ysim - sol
    res_sd = np.std(error)
    print('Standard deviation of fit residuals:', format(res_sd,'.2f'))
    return np.sum(error**2)

2. 正确调用minimize #

传入初始值，同时通过args传递所有固定参数：

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import minimize

# 假设time、concs、ysim为已定义的变量
initial_values = [300, 5e4, 4e-5]
fit = minimize(obj, initial_values, args=(time, ysim, time, concs))

# 查看优化结果
print('最优参数:', fit.x)
print('最优残差平方和:', format(fit.fun, '.2f'))

额外优化建议 #

• 若需要复用原langmuir_ode函数，可直接将参数通过solve_ivp的args传递，无需重新封装：

  def obj(params, t, ysim, time, concs):
      Rmax, kon, koff = params
      sol = solve_ivp(fun=langmuir_ode, 
                      t_span=[0, t[-1]], 
                      y0=[0.0], 
                      t_eval=t,
                      args=(Rmax, kon, koff))
      sol = sol.y.ravel()
      error = ysim - sol
      res_sd = np.std(error)
      print('Standard deviation of fit residuals:', format(res_sd,'.2f'))
      return np.sum(error**2)
  

• 添加参数边界约束，避免优化出无意义的负值：

  # 设置Rmax、kon、koff的下界均为正数
  bounds = [(1e-3, None), (1e-6, None), (1e-6, None)]
  fit = minimize(obj, initial_values, args=(time, ysim, time, concs), bounds=bounds)
  

内容的提问来源于stack exchange，提问作者Alex van der Kooi