python - 曲线拟合常微分方程

Question

我想为如下所示的 ODE 拟合一条曲线-

dA/dt = k1*profit + k2

我有观察到的变量 A 和 profit 的时间序列，我想获得 k1 和 的最优值>k2 在 python 中使用曲线拟合技术。我能够为此编写下面的代码，但解决方案不太合适，或者我的方法可能是错误的。

import numpy as np
from scipy.optimize import curve_fit
from scipy.integrate import odeint

def fitfunc(t, k1, k2):
    'Function that returns A computed from an ODE for k1 and k2'
    def myode(area, t):
        profit = get_profit(t)
        return k1*profit + k2

    A0 = 200000    #initial value of A
    out = odeint(myode, A0, t)
    return out[:,0]

k_fit, kcov = curve_fit(fitfunc, time_span, obs_A)  #time span from 1999-2019 and obs_A is the observed values of A

modeled_A = fitfunc(time_span, k_fit[0], k_fit[1])

20年期间的利润和obs_A数据为:

profit = [ 7.65976374e+06, -6.13172279e+06,  1.03946093e+07,  2.59937877e+06,
       -7.88358386e+06, -1.38918115e+04, -3.13403157e+06, -4.74348806e+06,
        1.87296164e+07,  4.13680709e+07, -1.77191198e+07,  2.39249499e+06,
        1.38521564e+07,  6.52548348e+07, -5.78102494e+07, -5.72469988e+07,
       -5.99056006e+06, -1.72424523e+07,  1.78509987e+07,  9.27860105e+06,
       -9.96709853e+06]

obs_A = [200000., 165000., 150000., 180000., 190000., 195000., 200000.,
       165000., 280000., 235000., 250000., 250000., 250000., 295000.,
       295000., 285000., 245000., 315000., 235000., 245000., 305000.]

time_span = np.arange(1999,2020)

这里的 get_profit 是一个输出给定 t 利润值的函数，它是通过对观察到的 profit 数据进行插值而创建的，如下面-

profit_fun = interp1d(t, profit.values, 1, fill_value="extrapolate")

def get_profit(t):
    return profit_fin(t)

我不确定如何在此处使用 profit 变量，因为它会在每个时间步发生变化。我的方法正确吗？

Answer 1

(根据要求，这是代码)

首先，进行设置。仅添加了 fitfun2，它修改了 fitfunc，删除了对 get_profit 的调用(因此不插入数据)。

import numpy as np
from scipy.optimize import curve_fit
from scipy.integrate import odeint
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

def fitfunc(t, k1, k2):  # Original
    'Function that returns A computed from an ODE for k1 and k2'
    def myode(area, t):
        profit = get_profit(t)
        return k1*profit + k2

    A0 = 20000    #initial value of A
    out = odeint(myode, A0, t)
    return out[:,0]

def fitfunc2(t, k1, k2): # Modified
    'Modified fitfunc, removing the call to `profit_fun`'
    def myode(area, t):
        return k1*t+k2

    A0 = 20000    #initial value of A
    out = odeint(myode, A0, t)
    return out[:,0]

profit = np.array([ 7.65976374e+06, -6.13172279e+06,  1.03946093e+07,  2.59937877e+06,
       -7.88358386e+06, -1.38918115e+04, -3.13403157e+06, -4.74348806e+06,
        1.87296164e+07,  4.13680709e+07, -1.77191198e+07,  2.39249499e+06,
        1.38521564e+07,  6.52548348e+07, -5.78102494e+07, -5.72469988e+07,
       -5.99056006e+06, -1.72424523e+07,  1.78509987e+07,  9.27860105e+06,
       -9.96709853e+06])

obs_A = np.array([200000., 165000., 150000., 180000., 190000., 195000., 200000.,
       165000., 280000., 235000., 250000., 250000., 250000., 295000.,
       295000., 285000., 245000., 315000., 235000., 245000., 305000.])

time_span = np.arange(1999,2020)

profit_fun = interp1d(time_span, profit, 1, fill_value="extrapolate")

def get_profit(t):
    return profit_fun(t)

现在，拟合并绘制结果

p0 = (1E-2, 1E4)
k_fit, kcov = curve_fit(fitfunc, time_span, obs_A, p0=p0)
k_fit2, kcov2 = curve_fit(fitfunc2, time_span, obs_A, p0=p0)

modeled_A = fitfunc(time_span, *k_fit)
guess_A = fitfunc(time_span, *p0)

modeled_A2 = fitfunc2(time_span, *k_fit2)
guess_A2 = fitfunc2(time_span, *p0)

plt.plot(time_span, obs_A, marker='o', lw=0, label='data')
plt.plot(time_span, modeled_A, label='model A (original)')
plt.plot(time_span, modeled_A2, label='model A (modified)')
plt.plot(time_span, guess_A, label='initial guess (original)')
plt.plot(time_span, guess_A2, label='initial guess (modified)')

plt.legend()

这是图表:

如前所述，修改参数 k 不会影响原始模型的曲线形状。它看起来仍然有点“逐步”。删除对 get_profit 的调用，曲线变得更加平滑，但我不知道这是不是@Isr729 所期望的。