1

我正在尝试在 python 中解决模型并将模型的未知参数拟合到实验数据。该模型由 2 个 ODE 组成，我使用 scipy.integrate.solve_ivp 对其进行求解。模型的参数是未知的，因此，我想使用多种方法来拟合它们。我的第一选择是差分进化，因为它可以为更复杂的模型提供非常好的结果（之前将它用作另一个包的一部分）。然而，问题是，当我给出我的模型（我希望它找到计算点和实验点之间的最小二乘的全局最小值）时，它会找到模型变得不稳定的参数（LSODA 无法整合它，因为内部步长为 0)。我试图捕捉 LSODA 在这些情况下抛出的运行时警告，但这并没有帮助。

我的模型和代码如下：

def aggregation_model(time, z, k1, k2, k3, k_1, k_2):
    GPVI_0 = 55.5
    GPVI, GPVI_Clust = z

    dGPVIdt = - k1 * GPVI_Clust * GPVI + k_1 * (GPVI_0 - GPVI) - 2 * k2 * GPVI * GPVI
    dGPVI_Clustdt = - (k_2 * GPVI_Clust + k3) * GPVI_Clust + k2 * GPVI * GPVI

    return [dGPVIdt, dGPVI_Clustdt]


def res_squares(parameters):

    time = np.linspace(0, 300, 1000)

    timepoints = [0, 25, 50, 75, 100, 125, 150, 300]
    val = [0, 2, 5, 10, 15, 20, 20, 18]

    model_calc = solve_ivp(aggregation_model, [0, 300], [55.5, 0], args=parameters, max_step=100000,
                           dense_output=True, method='LSODA', rtol=1e-12, atol=1e-6)

    solution = model_calc.sol(time)
    transposed = list(map(list, zip(*solution.T)))

    size = [0, ]
    for i in range(1, len(transposed[0])):
        size.append((55.5 - transposed[0][i]) / transposed[1][i])

    diff = []
    for i in range(len(timepoints)):
        indexval = min(range(len(time)), key=lambda j: abs(time[j] - timepoints[i]))
        diff.append((val[i] - size[indexval])**2)

    # print(np.sqrt(np.sum(diff)))

    output = np.sum(diff)

    return output

if __name__ == '__main__':

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


    parameterBounds = []
    parameterBounds.append([-1000000, 1000000])  # parameter bounds for a
    parameterBounds.append([-1000000, 1000000])  # parameter bounds for a
    parameterBounds.append([-1000000, 1000000])  # parameter bounds for a
    parameterBounds.append([-1000000, 1000000])  # parameter bounds for a
    parameterBounds.append([-1000000, 1000000])  # parameter bounds for a

    result = differential_evolution(res_squares, parameterBounds, seed=3, disp=True, init='random')
    print(result.x)

    sol = solve_ivp(aggregation_model, [0, 300], [55.5, 0], args=parvalues,
                    dense_output=True, method='LSODA', rtol=1e-6, atol=1e-12)
4

1 回答 1

1

有几种方法可以解决这个问题，可能需要将它们组合起来以获得稳定的解决方案：

  1. 避免导致不稳定数值行为的参数组合。您为参数设置的电流范围非常广泛。我相信您至少对参数的数量级有所了解，以便您可以更多地限制参数空间。如果您对设置硬边界感到不舒服，我建议您采用贝叶斯方法并设置先验分布，这些分布不会完全排除非常大或非常小的值，但确实会降低它们的可能性。（如果你选择后者，我可以推荐PyMC3包来构建贝叶斯模型，并结合Sunode进行 ODE 集成。）
  2. 尝试使解在数值上更稳定。以前在类似情况下对我有很大帮助的是 ODE 的对数转换。这意味着，而不是解决dy/dt = f(y)你解决dx/dt = exp(-x) f(exp(x))where x = log(y)。（这可以通过将微分的链式规则应用于dlog(y)/dt。）
  3. 在您的情况下，一个问题似乎是求解器在无法进一步集成时不会停止并抛出错误。相反，它似乎会永远继续尝试下一个集成步骤（？）。使用minstepLSODA 求解器的参数可能有助于解决这个问题。设置min_step=1e-14至少使优化器通过差分进化算法的前 11 个阶段。之后，求解器和优化器都不会打印任何消息，所以我不确定发生了什么。

顺便说一句，您可能想要传递t_eval=timepoints给，solve_ivp以便求解器已经在您需要的时间点返回解决方案。这将节省您当前所做的插值。

于 2020-12-15T14:12:33.787 回答