使用 ODE 系统进行参数优化

Posted 2023-02-18

【中文标题】使用 ODE 系统进行参数优化

【英文标题】：Parameter Optimisation with system of ODEs

【发布时间】：2022-01-07 11:11:34

【问题描述】：

我有一对 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 的限制是合理的给定系统，但是拟合仍然不起作用。残差函数有问题吗？

【问题讨论】：

您说“我相当肯定 a 和 c 的限制在系统的情况下是合理的，但是适合仍然不起作用。我”您担心的依据是什么？除了帮我解决这个问题，还有什么问题需要解决？

抱歉，这个问题不太清楚。 a 和 c 的限制是合理的，因此不需要更改，但是在绘制时，解决方案与测量数据不符并产生巨大的误差线（来自 report_fit）。所以问题是，这可能是由于残差函数造成的吗？还是整合？

【参考方案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)

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

【讨论】：

谢谢！我设法获得了与您共享的值相似的值。一个问题是我正在使用的数据应该在初始斜率中包含所有点（即在峰值之前），这会是 atol 参数的问题吗？

不，我不这么认为。我怀疑您需要调整其他参数。或者更长的时间序列的测量数据。有时有助于对模型中的变量进行归一化处理，以便求解器看到大部分大小为 0.1 到 100 或类似以 1 为中心的分量。求解器，尤其是“元”求解器 minimize，可能是没有详尽地测试它们的启发式部分如何对状态分量和残差中截然不同的尺度做出反应。