使用 ODE 系统进行参数优化

Question

我有一对 ODE，我目前正试图将它们拟合到我拥有的一个小数据集，但是我在优化两个参数（a 和 c）时遇到了一些问题。 ODE 采用稍微改变的 Lotka-Volterra 形式，由以下公式给出：

dT/dt = aT - bTL/(T+L+G)

dL/dt = cTL/(T+L+G) - dL

其中 b、G 和 d 是已知的，并且有一个小数据集可用于 T(t)（但不是 L(t)）。

目前我已经尝试使用 odeint 来求解方程组，并定义了要在 lmfit.Minimize 中使用的残差函数，但是无论在 lmfit.Minimize 中选择哪种方法，我的解决方案都不太适合数据改变边界。

目前的代码如下：

def eqns(y, t, paras):
    T,L=y
    
    try:
        c = paras['c'].value
        a = paras['a'].value

    except KeyError:
        c, a = paras

    b = 60*24    
    G = 1.7E9    
    d = 0.068
    
    return [a*T-b*L*(T/(G+T+L)),c*T*(L/(G+T+L))-d*L]

# Solution to differential equations T'(t) = model(t,x,paras), given initial condition T0
def sol(t, x0, paras):
    return odeint(eqns, x0, t, args=(paras,),rtol=1e-8,hmin=0.001,hmax=0.1)


# Define function to compute residuals
def residual(paras, t, data):
    arg0 = paras['T0'].value, paras['L0'].value
    model = sol(t, arg0, paras)
    x2_model = model[:, 0]
    return ((x2_model - data)**2).ravel()


# Set initial conditions
T0 = 50000
L0 = 1
y0 = [T0, L0]

# Measured data
t_measured = np.array([18,21,26,28,33])
T_measured = np.array([12667366.43,24917043.97,74910183.58,122959334.2,157298406])

# Set parameters including bounds
params = Parameters()
params.add('T0', value=T0, vary=False)
params.add('L0', value=L0, vary=False)
params.add('c', value=3, min=0.1, max=100)
params.add('a', value=0.2, min=0.01, max=10)

# Fit model
result = minimize(residual, params, args=(t_measured, T_measured), method='leastsq')

# Check effectiveness of fit
data_fitted = sol(np.linspace(0., 100, 1000), y0, result.params)

# Statistics of Fit
report_fit(result)

我相当肯定 a 和 c 的限制是合理的给定系统，但是拟合仍然不起作用。残差函数有问题吗？

Answer 1

有几个问题：

• odeint 积分器和任何其他数值积分器将初始条件 x0 应用于时间跨度或时间数组中给定的第一个数字。这里是t=18，从上下文来看是错误的，因为显然模拟应该从t=0 开始。有两种解决方案，

  • 将0 添加到时间数组t=np.insert(t,0,0.0) 并从结果return odeint(...)[1:] 中删除此值。这可能会导致绘图时间数组实际从 0 开始的时间出现问题。

  • 切换到solve_ivp，其中时间跨度和评估时间是不同的参数。

    def sol(t, x0, paras):
        res = solve_ivp(eqns, (0,max(t)+10), x0, t_eval=t, args=(paras,),\
                        atol=[1e-2,1e-4],rtol=1e-9)
        return res.y
    

    注意，结果是转置为odeint一个，所以你需要下一个x2_model = model[0]

• 残差不必要地平方。这可能会导致lmfit.minimize 内部的算法失效，因此请仅使用return (x2_model - data).ravel()

• 大家可能会注意到，我设置了atol 参数。偏离当前设置的范围，无论是作为单个值还是作为每个组件的值元组，都会导致误差估计有相当大的增加。

然后参数被估计为

    c:   19.7725209 +/- 0.62555507 (3.16%) (init = 3)
    a:   0.27960923 +/- 0.00146551 (0.52%) (init = 0.2)

看起来还不错，剧情也很合适