Python ODE拟合疑问：最小二乘法无结果或计算耗时过长

问题描述 #

尝试复现论文《Modeling Plasma Virus Concentration during Primary HIV Infection》的结果，需将4个参数（d、k、delta、p）拟合到给定数据集对应的微分方程（ODE）中，但运行scipy.optimize.least_squares或其他优化函数时始终无结果，已尝试调整初始猜测值为[0.013, 0.46e-2, 0.40, 980]。

原始代码 #

c = 3 #global constant 
v0= 10e-6

#dados
timeDataPatient1=np.array([0,22,43,78,106,146,183,230,268,358,435,489,519,534,584,610,687,778])
VirDataPatient1=10e3*np.array([v0*10e-3,27.2,210,85.9,81.1,46.2,60.1,82.8,103,72.1,79.4,70.4,207,42.6,10.8,54.2,22.3,40.8])


xdata = timeDataPatient1
ydata = VirDataPatient1

#inicial condition and max time of data set
IC = [10,0,10e-6]
t_max = max(xdata)

def diffeqs(t,y,d,k,delta,p):
    '''Model of the differencial equation'''
    global lam, c
    
    T=y[0]
    I=y[1]
    V=y[2]
    
    dT = 10*d - d*T-k*T*V
    dI = k*T*V - delta*I
    dV = p*I-c*V
    
    return [dT,dI,dV]

def Vtime(theta):
    '''solves the differencial equation, and returns the values of the model on the time points corresponding to the data set'''

    global t_max, IC, xdata
    
    sol = solve_ivp(diffeqs, (0,t_max), IC, args = (theta),t_eval=xdata)
    
    return sol.y[2]

#now I define the objetive function to minimize. For a parameter theta, that corresponds to d,k,delta,p in the ODE model,it calculates the difference between the log of the data given and the predicts by the ODE model.
  
def calcErrorHIV(theta):
    '''objetive function given in the paper'''
    global ydata
    
    dif = np.log(ydata)-np.log(Vtime(theta))
    
    return sum(dif**2)

#should return theta parameters that best fits the data, but It doesn't compute.
sp.optimize.least_squares(calcErrorHIV,[0.013,0.46e-2,0.40,980])

问题排查与修正点 #

1. 未定义的全局变量 #

diffeqs函数中声明了global lam但从未定义或使用该变量，会导致运行报错，直接删除无用的global lam声明即可。

2. 参数传递格式错误 #

solve_ivp的args参数要求传入元组，原代码中args=(theta)等价于args=theta（列表），无法正确拆分为4个参数传递给diffeqs函数，需修改为args=tuple(theta)。

3. 目标函数不符合least_squares要求 #

scipy.optimize.least_squares要求目标函数返回残差数组，而非残差平方和的标量。原函数返回sum(dif**2)，需改为直接返回残差数组dif，优化器会自动处理平方和最小化逻辑。

4. 缺失必要导入 #

代码中使用了np、solve_ivp、sp，需在代码开头添加导入语句：

import numpy as np
from scipy.integrate import solve_ivp
import scipy as sp

5. 添加参数范围约束（可选但推荐） #

HIV模型的参数有生物学意义的合理范围，可通过bounds参数限制参数值，避免优化过程中出现无意义的数值导致ODE求解失败：

• d：T细胞死亡率，范围(1e-3, 0.1)
• k：感染率，范围(1e-9, 1e-4)
• delta：感染细胞死亡率，范围(0.01, 2)
• p：病毒产生率，范围(10, 2000)

修正后的完整代码 #

import numpy as np
from scipy.integrate import solve_ivp
import scipy as sp

c = 3  # global constant 
v0 = 10e-6

# dados
timeDataPatient1 = np.array([0,22,43,78,106,146,183,230,268,358,435,489,519,534,584,610,687,778])
VirDataPatient1 = 10e3 * np.array([v0*10e-3,27.2,210,85.9,81.1,46.2,60.1,82.8,103,72.1,79.4,70.4,207,42.6,10.8,54.2,22.3,40.8])

xdata = timeDataPatient1
ydata = VirDataPatient1

# inicial condition and max time of data set
IC = [10, 0, 10e-6]
t_max = max(xdata)

def diffeqs(t, y, d, k, delta, p):
    '''Model of the differential equation'''
    global c
    
    T = y[0]
    I = y[1]
    V = y[2]
    
    dT = 10*d - d*T - k*T*V
    dI = k*T*V - delta*I
    dV = p*I - c*V
    
    return [dT, dI, dV]

def Vtime(theta):
    '''solves the differential equation, returns model values at data time points'''
    global t_max, IC, xdata
    
    # 将参数列表转为元组传递给ODE函数
    sol = solve_ivp(diffeqs, (0, t_max), IC, args=tuple(theta), t_eval=xdata)
    
    return sol.y[2]

def calcErrorHIV(theta):
    '''Objective function: returns log-scale residuals'''
    global ydata
    
    model_v = Vtime(theta)
    # 避免取对数时出现非正数
    if np.any(model_v <= 0):
        return np.full_like(ydata, np.inf)
    
    dif = np.log(ydata) - np.log(model_v)
    return dif

# 定义参数范围（生物学合理区间）
bounds = (
    [1e-3, 1e-9, 0.01, 10],    # 下限
    [0.1, 1e-4, 2, 2000]       # 上限
)

# 运行优化
result = sp.optimize.least_squares(calcErrorHIV, [0.013, 0.46e-2, 0.40, 980], bounds=bounds)
print("拟合参数：", result.x)
print("优化成功状态：", result.success)
print("状态信息：", result.message)

额外说明 #

• 4参数拟合本身确实需要一定计算时间，但如果修正代码后仍长时间无结果，可检查初始猜测值是否在合理范围，或尝试使用scipy.optimize.minimize配合L-BFGS-B算法（支持边界约束）替代least_squares。
• 代码中添加了对模型输出非正数的判断，避免取对数时出现错误，这是ODE求解过程中可能出现的问题（比如参数不合理导致病毒浓度变为负数）。

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