如何在不设置/假设初始条件的情况下将微分方程组拟合到数据？

Question

我是一名生物学学生，正在研究一些人口动态问题。我所掌握的数据只是一段时间内环境中捕食者和猎物数量的总和。我有两个耦合微分方程来描述每个群体的动态。我想做的是将这些方程与数据进行拟合，以使预测和测量（总）人口之间的误差最小化。

方程式如下：

dx/dt = Ax - Bx - Cx dy/dt = Ay + B*x

我想最小化如下定义的成本函数：

((x+y) - 数据)²

请注意，我没有单独的 x 和 y 数据，只有 T=x+y。因此我想在初始时间 t0 找到 x 和 y 的最佳猜测。

我目前使用的代码是：

    days = np.array([0,50,90,130,190,240,281,310,350,450])
    data = np.array([51,48,46,44,43,52,52,54,56,97])

    def system_ode(populations, time, params):
        x = populations[0]
        y = populations[1]

        A = params[0]
        B = params[1]
        C = params[2]

        dx_dt = x*(A - B - C)
        dy_dt = x*A + y*B

        return [dx_dt, dy_dt]

    def cost_func(params, time, data):

        init_populations = [data[0],0]

        time = np.linspace(days[0], days[-1], num = len(days))

        predicted_data = odeint(system_ode, init_populations, time, args=(params,))

        sum_pops = np.sum(predicted_data, axis=1)
        cost = np.sum(sum_pops - data)**2

        return cost


    init_params = [0.3, 0.3, 0.3] #initial parameter guesses

    opt_result = minimize(cost_func, init_params, args=(days, data))
    opt_params = opt_result.x.tolist()

    # now integrating with the optimized parameters
    init_populations = [data[0],0]
    time = np.linspace(days[0], days[-1], 100)

    pred_data_opt = odeint(system_ode, init_populations, time, args=(opt_params,))

    # Plotting 
    plt.plot(days, data[:], 'bo', label='Measured Total')
    plt.plot(time, pred_data_opt[:, 0], 'r-', label='Predicted X')
    plt.plot(time, pred_data_opt[:, 1], 'b-', label='Predicted Y')
    plt.xlabel('Time')
    plt.ylabel('Population')
    plt.legend()
    plt.show()

如图所示，这会导致非常次优的拟合，但我不确定为什么。此外，我目前正在强制曲线通过第一个点（因为我提供了初始条件并假设 Y 的总体在 t0 时为 0）。有没有办法在没有这个假设的情况下拟合这些方程？

Answer 1

我做了几件事：

• 我手动尝试了一些参数以进入正确的区域。参数的缩放很重要。
• dy/dt 的 ODE 实施错误。 A 和 B 交换了。
• 正如 Michael Cao 在对该问题的评论中所指出的，我将总和的平方修改为平方和。
• 同样根据 Michael 的输入，我添加了初始条件作为参数。
• 然后我将最小化方法从默认更改为Powell。如果我从头开始设置它，我建议使用

  least_squares

  而不是

  minimize

  。它更针对这个问题，因此更稳健。
• 我更改了目标中的时间评估点以匹配测量的时间实例。
• 我将总 X + Y 的图表添加到图中。

也许更多......这是代码：

import numpy as np from scipy.integrate import odeint from scipy.optimize import minimize import matplotlib.pylab as plt

days = np.array([0, 50, 90, 130, 190, 240, 281, 310, 350, 450]) data = np.array([51, 48, 46, 44, 43, 52, 52, 54, 56, 97], float)


def system_ode(populations, time, params):
    x, y = populations
    A, B, C = params

    dx_dt = x * (A - B - C)
    dy_dt = x * B + y * A
    return [dx_dt, dy_dt]


def cost_func(params, time, data):
    init_populations = params[3:]

    predicted_data = odeint(system_ode, init_populations, days,
                            args=(params[:3],))
    sum_pops = np.sum(predicted_data, axis=1)
    return np.sum((sum_pops - data) ** 2)


# A, B, x[0], y[0] init_params = [2e-3, 2e-3, 1e-2, data[0] / 2, data[0] / 2]

opt_result = minimize(cost_func, init_params, args=(days, data),
                      method="Powell") opt_params = opt_result.x.tolist() print("Optimal parameters and initial state: ", opt_params)

# now integrating with the optimized parameters init_populations = opt_params[3:] time = np.linspace(days[0], days[-1], 100)

pred_data_opt = odeint(system_ode, init_populations, time,
                       args=(opt_params[:3], )) sum_pops = np.sum(pred_data_opt, axis=1)

# Plotting plt.plot(days, data[:], 'ko', label='Measured Total') plt.plot(time, pred_data_opt[:, 0], 'r-', label='Predicted X') plt.plot(time, pred_data_opt[:, 1], 'b-', label='Predicted Y') plt.plot(time, sum_pops, 'k-', label='Predicted Total') plt.xlabel('Time') plt.ylabel('Population') plt.legend() plt.show()