Python GEKKO: Modeling large chemical reaction network in plug flow reactor

Question

I have a chemical reaction network including 65 reactions and 21 chemical species and I wanted to use GEKKO to model a simple plug flow reactor or to be more precise use GEKKO to solve the resulting DAE system.

The equations generally look like

dJ_i/dx = R_i or 0 = R_i

where J_i is the molar flux and R_i is the net reaction rate. R_i consist of multiple rates r_j calculated from the reaction network which generally looks like

R_i = sum( nu_ij * r_j ) with r = k * C_aⁿ * C_b^m

where C are my gas concentrations and also my unknowns I want to solve for. K, n and m are known from the reaction network and additional data files.

As the reaction network is too large and complicated to type every rate equation by hand, I wanted to use matrix and array operations like with numpy or write a "help function" to do intermediate calculations and formulate the final equations which I pass to the solver.

Normally with Scipy or comparable packages I would pass the following function to a solver:

P.L = ... # reactor length
P.C0 = ... # const. value for concentration
P.Jref= ... # reference molar flux
P.Rm =.. # given reaction matrix
k = ... # read rate constants from given additional data file 

P.a = P.RM.copy() # additional matrix to simplify rate calculation
P.a[P.a >0] = 0.
P.a = P.a * -1.

#J0 dimensionless molar flux of species i at x=0 (value between 0 and 1) 
J0 = np.zeros([21])
J0 [0] = ...
J0[1] = ...
...



def odefun(x,J):

    # dimless molar flux to concentration 
    c = ( (J[:]) / np.sum(J[:])) * P.C0
    C = np.tile(c,(np.shape(Para.Kin.RM)[0],1))

    r = P.k[:] * np.prod(np.power(X,P.a),1) # reaction rate

    R = np.transpose(P.RM) @ r # net rate of every chemical species

    dJ_dx =  P.L/P.J_ref * R[:] # ODE for every species

    return dJ_dx

Is there an equivalent why for GEKKO or a simple way to implement such a large reaction network ?

Answer 1

Many matrix operations in numpy are permitted with gekko arrays. The qualification is that numpy functions must return a symbolic result. Here are a few examples of matrix operations:

Matrix Operations with Numpy np.prod() and np.trace()

import numpy as np
from gekko import GEKKO

m = GEKKO(remote=False)

# Random 3x3
A = np.random.rand(3,3)
# Random 3x1
b = np.random.rand(3,1)
# Gekko array 3x3
p = m.Array(m.Param,(3,3))
# Gekko array 3x1
y = m.Array(m.Var,(3,1))

# Dot product of A p
x = np.dot(A,p)
# Dot product of x y
w = np.dot(x,y)
# Dot product of p y
z = np.dot(p,y)
# Trace (sum of diag) of p
t = np.trace(p)

# solve Ax = b
s = m.axb(A,b)
m.solve()

Matrix Operations to Define Constraints and Objective with m.sum()

from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
ni = 3; nj = 2; nk = 4
# solve AX=B
A = m.Array(m.Var,(ni,nj),lb=0)
X = m.Array(m.Var,(nj,nk),lb=0)
AX = np.dot(A,X)
B = m.Array(m.Var,(ni,nk),lb=0)
# equality constraints
m.Equations([AX[i,j]==B[i,j] for i in range(ni) \
                             for j in range(nk)])
m.Equation(5==m.sum([m.sum([A[i][j] for i in range(ni)]) \
                                    for j in range(nj)]))
m.Equation(2==m.sum([m.sum([X[i][j] for i in range(nj)]) \
                                    for j in range(nk)]))
# objective function
m.Minimize(m.sum([m.sum([B[i][j] for i in range(ni)]) \
                                 for j in range(nk)]))
m.solve()
print(A)
print(X)
print(B)

Here are related questions:

• GEKKO - optimization in matrix form
• Follow up question: GEKKO optimization in matrix form
• How to use gekko variables as componets of an array?
• Using Gekko to minimize norm of matrix

The Gekko underlying modeling language (APMonitor) has been used for other PDE systems such as a Solid Oxide Fuel Cells (SOFCs) with reaction kinetics, axial / radial spatial dependence, and time dependence. I recommend IMODE=7 to speed up the solution if you are only simulating and not optimizing. For the SOFC case, converting from Simulink to APMonitor (Gekko precursor) sped up the simulation significantly. We later applied the simulator in Model Predictive Control.

Jacobsen, L. T. and Hedengren, J. D., Model Predictive Control with a Rigorous Model of a Solid Oxide Fuel Cell, American Control Conference (ACC), Washington, DC, 2013.

One other strategy for large reaction networks is to focus on intialization strategies. Using m.options.COLDSTART=2 performs a decomposition of the problem to solve sub-sets of the model that can be successively solved independently. Additional details are provided here:

Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, Computers and Chemical Engineering, 2015, Vol. 78, pp. 39-50, DOI: 10.1016/j.compchemeng.2015.04.016.