问题描述
我的模型是 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,))
~\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,243 ixpr,--> 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])
解决方法
您收到 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,beta] # your best guess for the parameter values
bounds = [(None,None),(0,(None,None)] # alpha>0 so V_nf!=0
res = minimize(SumSquaredErr,x0,bounds=bounds)
print('best [k,beta]: ',res.x)
我在这里假设 ts = np.linspace(0,60,61) 以及 C1 和 C2 是时间 ts[1:] 的数据，即 1,2,3,...,59,60。