Python GEKKO: моделирование химической реакции - PullRequest
Новая
       11

Python GEKKO: моделирование химической реакции

4 голосов
/ 20 февраля 2020

Я использую Python GEKKO для моделирования химической реакции, которую можно описать так:

1 -> 2 -> 3 -> 4

с побочными реакциями как следует:

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 обозначает прогнозируемые значения, пунктирные линии представляют собой кривые, которые рассчитываются с использованием единой рассчитанной константы скорости.

Concentration plot

Вопрос в том, откуда возникает это отклонение? Я не могу найти источник.

Спасибо! 

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()

1 Ответ

2 голосов
/ 22 февраля 2020

Проблема решается путем увеличения количества узлов для каждого временного шага. 

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

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

Improved accuracy

Gekko требует от пользователей указывать временные точки где они хотели бы вычислить решение. В вашем случае они соответствуют измерениям. Вам просто нужен хотя бы один внутренний узел (NODES=3), чтобы повысить точность интеграции ODE. Точность, как правило, увеличивается с увеличением числа узлов (2-6), но при этом также увеличиваются вычислительные затраты, поскольку модель рассчитывается на всех узлах. Более подробную информацию об ортогональном расположении конечных элементов можно найти в курсе Машинное обучение и динамическая оптимизация c .

...