Python：SciPy Curve_Fit用于将ODE系统拟合到化学反应动力学数据

如何解决Python：SciPy Curve_Fit用于将ODE系统拟合到化学反应动力学数据

我正在尝试建立一种将ODE系统适合我的模拟化学反应数据的方法。我正在尝试找到一种构建工具来从实验数据中查找某些参数（速率常数，例如k1_f，k1_r，k2_f等）的方法。

据我了解，这些速率常数/参数是针对所追踪的每种物质/化学物质的化学反应速率/ ODE定义的。因此，如果我们将原始数据适合ODE，则可以提取这些参数。而且，我必须能够在没有每种物质/化学物质数据的情况下进行此操作（因为有时无法通过实验确定这些数据）。

我在下面尽了最大的努力，从Kitchin组网站上构建了一个简单的示例，但无法弄清楚如何使用curve_fit将ODE系统拟合到数据中。

我正在尝试解决此问题，部分是为了找到可行的方法，另一方面是为了更好地了解合适的ODE的工作原理。

下面是我正在处理的代码：

#Defining initial Data
tspan = [0.,0.6,1.2,1.8,2.4,3.,3.6,4.2,4.8]
Prod_data = [0.,0.60420736,0.73895125,0.7591921,0.76204478,0.76244296,0.76249842,0.76250616,0.76250722]

#Defining Function to get fit
def fitfunc1(t,k):
    k1f = k[0]
    k1r = k[1]
    k2f = k[2]
    k2r = k[3]
    
    def oderxn(C,t):
        #Defining differential form for each species
        Ce = C[0] 
        Cs = C[1]
        Ces = C[2]
        Cp = C[3]
        
        #Defining rate laws for forward and reverse reaction
        frxn1 = k1f * Ce * Cs
        rrxn1 = k1r * Ces
        frxn2 = k2f * Ces
        rrxn2 = k2r * Ce * Cp

        dEdt = -frxn1+rrxn1+frxn2-rrxn2
        dSdt = -frxn1+rrxn1
        dESdt = frxn1-rrxn1-frxn2+rrxn2
        dPdt = frxn2-rrxn2
        return [dEdt,dSdt,dESdt,dPdt]

    C0 = [0.4,1,0]
    Csol = odeint(oderxn,C0,t) 
    return Csol[:,0] 

C_int = (0.4,Prod_data)

print(C_int)
    
k_fit_a = curve_fit(fitfunc1,tspan,C_int)
print (k_fit_a)


tfit = np.linspace(0,5);
fit = fitfunc1(tfit,k_fit_a)

import matplotlib.pyplot as plt
plt.plot(tspan,A_data,'ro',label='data')
plt.plot(tfit,fit,'b-',label='fit')
plt.legend(loc='best')

以下是使用Symfit的示例（我必须开始工作），以帮助解决问题：

from symfit import variables,parameters,ODEModel,D,Fit
from symfit.core.support import key2str
import numpy as np
import matplotlib.pyplot as plt

#Defining initial Data
tspan = [0.,0.76250722]

# Define our ODE model
t,E,S,ES,P = variables('t,P')
k1f,k1r,k2f,k2r = parameters('k1f,k2r')

S_0 = 1 #Inital [S]0 from above
E_0 = 0.4 #Initial [E]0 from above

frxn1 = k1f * E * S
rrxn1 = k1r * ES
frxn2 = k2f * ES
rrxn2 = k2r * E * P

model_dict = {
    
    D(E,t): -frxn1+rrxn1+frxn2-rrxn2,D(S,t): -frxn1+rrxn1,D(ES,t): frxn1-rrxn1-frxn2+rrxn2,D(P,t): frxn2-rrxn2,}
model = ODEModel(
    model_dict,initial={t: 0.0,E: E_0,S: S_0,ES: 0.0,P: 0.0}
)
print(model)

#Fitting the Data
#k1_f.min,k1_f.max = 1,20
#k1_r.min,k1_r.max = 1,20
#k2_f.min,k2_f.max = 0,20
#k2_r.min,k2_r.max = 0,20

fit = Fit(model,t=tspan,E= None,S=None,ES=None,P=Prod_data) 
fit_result = fit.execute()
print(fit_result)

#Graphing Results
taxis = np.linspace(0,5,1000)
model_fit = model(t=taxis,**fit_result.params)._asdict()
plt.scatter(tspan,Prod_data,label='[P] Data',color='blue')
colors = ('red','blue','green','yellow')
for var in model:
    plt.plot(taxis,model_fit[var],label='[{}]'.format(var.name),linewidth=2.0)
plt.legend()
plt.show()