Python GEKKO: Right use of 2D Intermediates

Question

I have recently asked a about modeling a plug flow reactor with a large reaction network using GEKKO where I had problems with the intermediate calculations. I now tried to use the m.Intermediate and also found an early post of how to create 2D arrays of Intermediates, but my code got very lengthy and complicated and still does not work, because the data arrays in my final equations don't change their length according to the time steps. I also think, that the exzessive and complicated use of Intermediates is not their intended use, but I can't find another way to do all my necessary calculations as the reaction network is very large.

I would be very thankful if anyone could help me or give me a hint. This would be my code:

dae = GEKKO(remote=False)
dae.time = np.linspace(0., 1.,100)

# parameters
k = dae.Param(Para.Kin.Rate_const) # rate constants numpy array with 21 elements
C0 = dae.Param(Para.C0) # initial concentration 
Sites = dae.Param(Para.Sites) # surface sites concentration
RMT = dae.Param(np.transpose(Para.Kin.RM)) # reaction matrix numpy 2D array
a = dae.Param(Para.a) # help matrix numpy 2D array
SV = dae.Param(Para.SV) # specific surface
Pflag_gas = dae.Param(Para.Kin.flag_gas)
Pflag_surface = dae.Param(Para.Kin.flag_surface)
Pflag_solid = dae.Param(Para.Kin.flag_solid)
Pflag_r_gas = dae.Param(Para.Kin.flag_r_gas)
dim = dae.Param(Para.L/Para.J_ref)


# intermediates
J = [None]*21
S = [None]*21
Si = [None]*21
c = [None]*21
x = [None]*21
X = [[None for x in range(21)] for y in range(65)]
XX = [None]*65
r = [None]*65
R_help=[None]*65
R = [None]*21

# variable
Y_gekko = dae.Array(dae.Var,21,lb=0.)

for i in range(21):
    Y_gekko[i].value = Para.Y0[i]

# flags are numpy arrays containing 0 or 1 and have the same length as Y_gekko
for i in range(21):
    J[i] = dae.Intermediate(Y_gekko[i] * Pflag_gas[i])
    S[i] = dae.Intermediate(Y_gekko[i] * Pflag_surface[i])
    Si[i] = dae.Intermediate(Y_gekko[i] * Pflag_solid[i])

for i in range(21):
    c[i] = dae.Intermediate( ( (J[i]) / dae.sum(J)) * C0)
                   
    x[i] = dae.Intermediate( c[i] + Si[i] + S[i] * Sites)

    for j in range(65):
        X[j][i] =  dae.Intermediate(x[i])

for i in range(65):
    for j in range(21):
        XX[i] = dae.Intermediate(X[i][j]**a[i][j])
        
for i in range(65):
    r[i] = dae.Intermediate(SV**(1-Pflag_r_gas) * k[i] * XX[i])
    
for i in range(21):   
    for j in range(65):
        R_help[j] = dae.Intermediate(RMT[i][j] * r[j])
        
    R[i] = dae.Intermediate(dae.sum(R_help))
        

# solve DAE system
dae.options.IMODE = 7
    
dae.Equation(Y_gekko[0].dt() == dim * R[0])
dae.Equation(1 == Y_gekko[1] + Y_gekko[2]+ Y_gekko[3] + Y_gekko[4])
dae.Equation(0 == R[2])
dae.Equation(0 == R[3])
dae.Equation(0 == R[4])
dae.Equation(Y_gekko[5].dt() == dim * R[5])
dae.Equation(Y_gekko[6].dt() == dim * R[6])
dae.Equation(Y_gekko[7].dt() == dim * R[7])
dae.Equation(Y_gekko[8].dt() == dim * R[8])
dae.Equation(Y_gekko[9].dt() == dim * R[9])
dae.Equation(Y_gekko[10].dt() == dim * R[10])
dae.Equation(Y_gekko[11].dt() == dim * R[11])
dae.Equation(Y_gekko[12].dt() == dim * R[12])
dae.Equation(Y_gekko[13].dt() == dim * R[13])
dae.Equation(Y_gekko[14].dt() == dim * R[14])
dae.Equation(Y_gekko[15].dt() == dim * R[15])
dae.Equation(Y_gekko[16].dt() == dim * R[16])
dae.Equation(Y_gekko[17].dt() == dim * R[17])
dae.Equation(Y_gekko[18].dt() == dim * R[18])
dae.Equation(Y_gekko[19].dt() == dim * R[19])
dae.Equation(Y_gekko[20].dt() == dim * R[20])

dae.solve()

Answer 1

Intermediates are only needed when an equation symbolic expression exceeds 15,000 characters (long equation), intermediate values are reused in multiple expressions, or when the intermediate values are of interest in the solution. There are a few issues with the script:

1. Use x[i] directly in next expression. Use XX[j][i] (no intermediate needed). Don't redefine intermediate variables once they are defined once.

Original

for i in range(21):
    c[i] = dae.Intermediate( ( (J[i]) / dae.sum(J)) * C0)
    x[i] = dae.Intermediate( c[i] + Si[i] + S[i] * Sites)
    for j in range(65):
        X[j][i] =  dae.Intermediate(x[i])
for i in range(65):
    for j in range(21):
        XX[i] = dae.Intermediate(X[i][j]**a[i][j])

Modified

X = [[None for x in range(21)] for y in range(65)]
for i in range(21):
    c[i] = dae.Intermediate( ( (J[i]) / dae.sum(J)) * C0)
    x[i] = dae.Intermediate( c[i] + Si[i] + S[i] * Sites)
for i in range(65):
    for j in range(21):
        XX[j][i] = x[i]**a[i][j])

2. Use a list comprehension instead of a loop.

Original

for i in range(21):   
    for j in range(65):
        R_help[j] = dae.Intermediate(RMT[i][j] * r[j])
    R[i] = dae.Intermediate(dae.sum(R_help))

Modified

for i in range(21):   
    R[i] = dae.sum([RMT[i][j]*r[j] for j in range(65)])

You can also try using sum() instead of m.sum() to improve the compile time. Use of intermediates is acceptable, but not needed in these cases.