Python问题：ODE系统的参数估计，使得系数是函数

Question

我的模型是 ODE 系统，系数是具有参数的函数。

我的目标是系数中的参数估计。

dC1dt = c_11(t)*C1 + c_12(t)*C2

dC2dt = c_22(t)*C2 + c_22(t)*C2

我知道如何用 R 解决它，但我必须用 Python 解决它。

首先，我尝试在固定条件下使用 odeint。

完整代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# fixed values
V_room = 40 #room volume[m^3]
v = 0.01 #velocity of material [m/s]
h = 2.5 # height of room [m/s]
Q = 1 # average flux [m^3/min]
V_0 =1 # initial value of nf [m^3]
Q_0 = 0.1 #initial value of Q_ff

#initial values 
C_0 = [0.0002, 2e-05]

#parameter 
k= 0.06       
beta = 0.001
alpha = 0.005

# measured time
ts = np.linspace(0,60,60)
# fixed functions

def V_nf(t):
    return V_0 + (V_room - V_0)*(1-np.exp(-alpha*t))

def V_ff(t):
    return V_room - V_nf(t)

def dVnf(t):
    return alpha*(V_nf(t)-V_0)

def dVff(t):
    return -dVnf(t)

def S_nf(t):
    return V_nf(t)**(2/3)

def Q_ff(t):
    return Q_0*np.exp(-beta*t)

=====

coefficients = {'coef_11' : coef11, 'coef_12' : coef12, 'coef_21' : coef21, 'coef_22' : coef22}

def coef11(t):
    return (-k*S_nf(t)-dVnf(t))/V_nf(t) - Q/V_room - v/h

def coef12(t):
    return (k*S_nf(t)+Q_ff(t))/V_nf(t) - Q*V_ff(t)/(V_room * V_nf(t))

def coef21(t):
    return k*S_nf(t)/V_ff(t)

def coef22(t):
    return (-dVff(t)-Q_ff(t)-k*S_nf(t))/V_ff(t) - v/h

def dCdt(C,t):
    
    C_nf = C[0]
    C_ff = C[1]
    c_11 = coefficients['coef_11'](t)
    c_12 = coefficients['coef_12'](t)
    c_21 = coefficients['coef_21'](t)
    c_22 = coefficients['coef_22'](t)
    
    dCnf = c_11* C_nf + c_12 * C_ff
    dCff = c_21* C_nf + c_22 * C_ff
    
    return [dCnf, dCff]

C = odeint(dCdt, C_0, ts, args= (coefficients,))

但是，我得到了

---

TypeError                                 Traceback (most recent call last)
<ipython-input-114-7cece3e6a166> in <module>
     13     return [dCnf, dCff]
     14 
---> 15 C = odeint(dCdt, C_0, ts, args= (coefficients,))

~\Anaconda3\lib\site-packages\scipy\integrate\odepack.py in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst)
    242                              full_output, rtol, atol, tcrit, h0, hmax, hmin,
    243                              ixpr, mxstep, mxhnil, mxordn, mxords,
--> 244                              int(bool(tfirst)))
    245     if output[-1] < 0:
    246         warning_msg = _msgs[output[-1]] + " Run with full_output = 1 to get quantitative information."

TypeError: dCdt() takes 2 positional arguments but 3 were given

 

不知道这个方法是否正确，为什么会出现这个问题。

正如我所说，我的目标是估计参数（k，alpha，beta）。

因此，如果您知道如何做到这一点，请使用数据。

数据：

C1 = array([2.00000e-04, 3.52006e+00, 1.01378e+00, 8.47760e-01, 6.19000e-01,
   6.09940e-01, 4.35010e-01, 3.49150e-01, 3.87830e-01, 3.24830e-01,
   1.97040e-01, 1.84630e-01, 1.72520e-01, 1.35980e-01, 1.38430e-01,
   1.21520e-01, 1.46680e-01, 9.08900e-02, 1.09650e-01, 8.71000e-02,
   9.24400e-02, 1.11200e-01, 7.96600e-02, 9.26300e-02, 8.08700e-02,
   8.04000e-02, 7.69200e-02, 9.06400e-02, 7.30600e-02, 7.22200e-02,
   6.57900e-02, 7.87200e-02, 7.67700e-02, 7.28100e-02, 6.65600e-02,
   6.87300e-02, 7.65600e-02, 6.87200e-02, 7.52700e-02, 8.66000e-02,
   6.90700e-02, 6.51900e-02, 5.39200e-02, 6.46600e-02, 5.96400e-02,
   6.90100e-02, 6.23700e-02, 1.00230e-01, 9.05900e-02, 4.96300e-02,
   6.75400e-02, 5.50200e-02, 4.75700e-02, 5.82800e-02, 5.34400e-02,
   5.30200e-02, 4.10700e-02, 4.10800e-02, 4.91300e-02, 4.16300e-02])

C2 = array([5.000e-05, 2.000e-05, 1.000e-04, 1.250e-03, 1.500e-03, 2.630e-03,
   2.600e-04, 5.540e-03, 2.160e-03, 9.420e-03, 8.030e-03, 1.369e-02,
   9.620e-03, 1.527e-02, 1.209e-02, 9.600e-03, 1.081e-02, 1.414e-02,
   1.522e-02, 1.244e-02, 2.223e-02, 2.312e-02, 2.683e-02, 2.398e-02,
   2.658e-02, 2.841e-02, 2.123e-02, 3.052e-02, 3.349e-02, 4.027e-02,
   3.110e-02, 3.134e-02, 3.185e-02, 2.964e-02, 4.078e-02, 2.530e-02,
   4.443e-02, 2.856e-02, 2.459e-02, 2.568e-02, 2.104e-02, 2.477e-02,
   2.297e-02, 2.512e-02, 2.133e-02, 2.002e-02, 1.427e-02, 3.033e-02,
   1.946e-02, 2.173e-02, 1.800e-02, 1.346e-02, 2.039e-02, 2.132e-02,
   1.416e-02, 1.376e-02, 1.079e-02, 6.640e-03, 1.324e-02, 1.056e-02])

Answer 1

您收到 TypeError 是因为您忘记了 dCdt 函数的定义中的参数。应该是def dCdt(C,t,coefficients):。

至于参数优化，有很多不同的方法和途径。你在R中使用了什么方法？在 python 中找到良好参数的最简单方法之一是使用scipy.optimize.minimize 最小化误差平方和。

def SumSquaredErr(x):
    k,alpha,beta = x
    
    # Insert all your function definitions here #
    
    C = odeint(dCdt, C_0, ts, args= (coefficients,))

    return ((C[1:,0] - C1)**2).sum() + ((C[1:,1] - C2)**2).sum()

x0 = [k,alpha,beta] # your best guess for the parameter values
bounds = [(None,None),(0,None),(None,None)] # alpha>0 so V_nf!=0
res = minimize(SumSquaredErr, x0, bounds=bounds)
print('best [k,alpha,beta]: ', res.x)

我在这里假设 ts = np.linspace(0,60,61) 和 C1 和 C2 是时间 ts[1:] 的数据，即 1,2,3,...,59,60。