python GEKKO : Modelling a chemical reaction

Question

我正在使用 Python GEKKO 来模拟化学 react ，可以这样描述:

1 -> 2 -> 3 -> 4

副 react 如下:

2 -> 5

3 -> 5

产物(4)是稳定的。这导致了以下一组 ODE(速率方程)，具有速率常数 k 和组分 c(i) 的浓度。

m.Equation(c_1p.dt() == -k1*c_1p)
m.Equation(c_2p.dt() == k1*c_1p - (k21 + k22)*c_2p)
m.Equation(c_3p.dt() == k21*c_2p - (k31 + k32)*c_3p)
m.Equation(c_4p.dt() == k31*c_3p)

我在 GEKKO 中实现了这些方程来估计速率常数。我将测量值​​初始化为参数，将预测值初始化为变量，将速率常数初始化为固定值(有限制)。目标函数是最小二乘法。这工作正常，结果在预期范围内(R2 > 0.99)。

问题是，当我尝试通过使用计算出的速率常数来求解 ODE(使用 GEKKO 或 scipy odeint)来验证这些结果时，我得到了不同的结果(见图 1)。点是测量值，X 标记预测值，虚线表示使用 odeint 使用计算的速率常数计算的曲线。



问题是，这种偏差从何而来？我找不到来源。

谢谢!

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

m = GEKKO(remote=False)
p = GEKKO(remote=False)

##  Arrays: Measurement Data

c_1 = [29.99122982,  1.24267302,1,1,1,1]
c_2 = [92.163874642, 34.957203757, 15.868747086, 0.0, 0.0, 0.0]
c_3 = [1.5290222821, 7.3374483783, 3.1998472876, 0.66149919406, 0.069616016241, 0.0071280867007]
c_4 = [0.0, 38.487210404, 51.285418999, 57.66299199, 60.869003531, 64.89535785]

m.time = [0,15,30,60,120,180]
t = m.time

p.time = np.linspace(0,180,181)
t_p = p.time


##Variables

c_1m = m.Param(value=c_1)      
c_2m = m.Param(value=c_2)
c_3m = m.Param(value=c_3)
c_4m = m.Param(value=c_4)

c_1p = m.Var(value=c_1m)          
c_2p = m.Var(value=c_2m)
c_3p = m.Var(value=c_3m)
c_4p = m.Var(value=c_4m)

c_1pp = p.Var(value = c_1[0])    
c_2pp = p.Var(value = c_2[0])
c_3pp = p.Var(value = c_3[0])
c_4pp = p.Var(value = c_4[0])

k1 = m.FV(lb = 0, ub = 2)
k1.STATUS = 1
k21 = m.FV(lb = 0)
k21.status = 1
k22 = m.FV(lb = 0)
k22.status = 1
k31 = m.FV(lb = 0)
k31.status = 1
k32 = m.FV(lb = 0)
k32.status = 1

##m.Equations

m.Equation(c_1p.dt() == -k1*c_1p)
m.Equation(c_2p.dt() == k1*c_1p - k21*c_2p - k22*c_2p)
m.Equation(c_3p.dt() == k21*c_2p - k31*c_3p - k32*c_3p)
m.Equation(c_4p.dt() == k31*c_3p)


##Objective

m.Obj((c_4p - c_4m)**2 + (c_3p-c_3m)**2 + (c_2p-c_2m)**2 + (c_1p-c_1m)**2)

##Application options

m.options.IMODE = 5
m.options.solver = 3

p.options.IMODE = 4

##Solve

m.solve()

##p.Equations

p.Equation(c_1pp.dt() == -k1[0]*c_1pp)
p.Equation(c_2pp.dt() == k1[0]*c_1p - k21[0]*c_2pp - k22[0]*c_2pp)
p.Equation(c_3pp.dt() == k21[0]*c_2p - k31[0]*c_3pp - k32[0]*c_3pp)
p.Equation(c_4pp.dt() == k31[0]*c_3pp)

p.solve()

def rxn(C,t_p):

    c_10 = C[0]
    c_20 = C[1]
    c_30 = C[2]
    c_40 = C[3]

    d1dt = -k1[0] * c_10
    d2dt = k1[0] * c_10 - (k21[0] + k22[0]) * c_20
    d3dt = k21[0] * c_20 - (k31[0] + k32[0]) * c_30
    d4dt = k31[0] * c_30

    return [d1dt, d2dt, d3dt, d4dt]

C = [29.991229823,92.163874642,1.5290222821,0.0]
model = odeint(rxn,C,t_p)

##Plot/Print

print('c_4m = ' +  str(c_4m.VALUE))
print('c_4p = ' +  str(c_4p.VALUE))
print('c_3m = ' +  str(c_3p.VALUE))
print('c_3p = ' + str(c_3m.VALUE))
print('c_2m = ' + str(c_2m.VALUE))
print('c_2p = ' + str(c_2p.VALUE))
print('c_1p = ' +str(c_1p.VALUE))
print('c_1m = ' + str(c_1m.value))
print('k1 = ' + str(k1.VALUE))
print('k21 = ' + str(k21.VALUE))
print('k22 = ' + str(k22.VALUE))
print('k31 = ' + str(k31.VALUE))
print('k32 = ' + str(k32.VALUE))


plt.figure(1)
plt.plot(t,c_1m,'ro', label = "1 meas")
plt.plot(t,c_1p, 'rx', label = "1 pred")
plt.plot(p.time, model[:,0] ,'r--', label= "1 mod")
plt.plot(t,c_2m,'go', label = "2 meas")
plt.plot(t,c_2p,'gx', label = "2 pred")
plt.plot(p.time,model[:,1],'g--', label = "2 mod")
plt.plot(t,c_3m,'yo', label = "3 meas")
plt.plot(t, c_3p, 'yx', label= "3 pred")
plt.plot(p.time,model[:,2],'y--', label = "3 mod")
plt.plot(t,c_4m,'bo', label="4 meas")
plt.plot(t,c_4p,'bx', label = "4 pred")
plt.plot(p.time,model[:,3],'b--', label = "4 mod")
plt.xlabel('Time [min]')
plt.ylabel('Concentration [mmol/l]')
plt.legend(loc='best', ncol = 4)
plt.grid()
plt.show()

Answer 1

该问题通过增加每个时间步的节点数来解决。

##Application options
ND = 3
m.options.IMODE = 5
m.options.NODES = ND

p.options.IMODE = 4
p.options.NODES = ND

Gekko 要求用户指定他们想要计算解决方案的时间点。在您的情况下，它们对应于测量值。您只需要至少一个内部节点 ( NODES=3 ) 来提高 ODE 集成的准确性。准确度通常随着更多节点 (2-6) 而增加，但也会增加计算成本，因为模型是在所有节点上计算的。在 Machine Learning and Dynamic Optimization course 中有更多关于有限元正交搭配的信息。 .