I am trying to do a dynamic simulation of ammonia reactor using GEKKO. Unfortunately, my code error says that it can't reach a solution. The equations I am trying to solve are steady-state in terms of species continuity equations but dynamic in terms of heat balance. I used the code shown below to first solve for the steady-state condition, which it does successfully by setting the IMODE = 1. However, when I import the steady-state results and used it to initialise my dynamic-state code, no solution is achieved. I do not understand how my code is not working when solving for the dynamic model since it worked fine for the steady-state model. I have tried IMODE = 4 and 7 but neither of them worked. Does anyone have any idea why my code is not working for the dynamic simulation?
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
from GEKKO_single import n_store, P_store, xH2_store, xN2_store, xNH3_store, xinert_store, T_store, nr_store, Pr_store, xr_H2_store, xr_N2_store, xr_NH3_store, xr_inert_store, Tr_store
m = GEKKO()
#REACTOR BED PARAMETERS
stoi = [-3, -1, 2, 0] # Reaction Stoichiometric coefficient [H2, N2, NH3, Ar]
Vbed_R1 = m.Const(0.7) # Bed Volume [m3]
Vbed_R2 = m.Const(1.3)
rhocat = m.Const(2200) # Density of catalyst [kg/m3]
Cp = m.Const(35500) # Heat capacity of gas [J/(kmol K)]
Cpcat = m.Const(1100) # Heat capacity of catalyst [J/(kg K)]
dHrx = m.Const(91.9836e6) # Heat of reaction [J/kmol]
R = m.Const(8.314) # Gas Constant [J/(mol K)]
#KINETIC PARAMETERS
Afor = m.Const(1.79e4) # Forward reaction activity coefficient [kmole/(m3.hr.atm1.5)]
Abac = m.Const(2.57e16) # Backward reaction activity coefficient [kmole.atm0.5/(m3.hr)]
Eafor = m.Const(87090) # Forward reaction activation energy [J/mol]
Eabac = m.Const(198464) # Backward reaction activation energy [J/mol]
# SPLITTER CONDITIONS
u1 = m.Const(0.2302)
u2 = m.Const(0.2659)
u3 = m.Const(0.5039)
#HEAT EXCHANGER
Cp_c = Cp
Cp_h = Cp
A_HX2 = m.Const(50)
U = m.Const(536)
# FEED CONDITIONS
nf = m.Const(0.4072)
Pf = m.Const(150.2780)
Tf = m.Const(403.08370) # Temperature [K]
xf_H2 = m.Const(0.7050) # Mole fraction H2 [-]
xf_N2 = m.Const(0.2398) # Mole fraction N2 [-]
xf_NH3 = m.Const(0.0351) # Mole fraction NH3 [-]
xf_inert = m.Const(0.0202) # Mole fraction inert [-]
#VOLUMETRIC COMPARTMENTS
ns_R1 = 20
ns_R2 = 20
ns = ns_R1 + ns_R2
#MASS OF THE CATALYST
mcat_1 = (Vbed_R1/ns_R1)*rhocat # Mass of the catalyst [kg]
mcat_2 = (Vbed_R2/ns_R2)*rhocat
#%% DEFINITIONS
#CONNECTION NUMBER
split_1 = 0
split_2 = 1
split_3 = 2
HX2_o = 3
R1_i = 4
R1_o = 5
R2_i = 6
R2_o = 7
total = R2_o + 1
nt = int(3600*2) #seconds
m.time = np.linspace(0,nt,nt) #seconds
Pf = m.MV()
Pf.value = np.ones(nt) * Pf
Pf.value[1800:] = 120 #bar
#%% DEFINITION AND INITIALIZATION OF PARAMETERS
# STEADY-STATE
n = [m.Var(value = n_store[i], lb = 0) for i in range(total)]
P = [m.Var(value = P_store[i], lb = 0) for i in range(total)]
xH2 = [m.Var(value = xH2_store[i], ub = 1, lb = 0) for i in range(total)]
xN2 = [m.Var(value = xN2_store[i], ub = 1, lb = 0) for i in range(total)]
xNH3 = [m.Var(value = xNH3_store[i], ub = 1, lb = 0) for i in range(total)]
xinert = [m.Var(value = xinert_store[i], ub = 1, lb = 0) for i in range(total)]
T = [m.Var(value = T_store[i], lb = 0) for i in range(total)]
# REACTOR BEDS
n_r = [m.Var(value = nr_store[i], lb = 0) for i in range(ns)]
P_r = [m.Var(value = Pr_store[i], lb = 0) for i in range(ns)]
xH2_r = [m.Var(value = xr_H2_store[i], ub = 1, lb = 0) for i in range(ns)]
xN2_r = [m.Var(value = xr_N2_store[i], ub = 1, lb = 0) for i in range(ns)]
xNH3_r = [m.Var(value = xr_NH3_store[i], ub = 1, lb = 0) for i in range(ns)]
xinert_r = [m.Var(value = xr_inert_store[i], ub = 1, lb = 0) for i in range(ns)]
T_r = [m.Var(value = Tr_store[i], lb = 0) for i in range(ns)]
#%% SPLITTER
#SPLITTER 1
m.Equation( n[split_1] == u1 * nf )
m.Equation( P[split_1] == Pf )
m.Equation( xH2[split_1] == xf_H2 )
m.Equation( xN2[split_1] == xf_N2 )
m.Equation( xNH3[split_1] == xf_NH3 )
m.Equation( xinert[split_1] == xf_inert )
m.Equation( T[split_1] == Tf )
#SPLITTER 2
m.Equation( n[split_2] == u2 * nf )
m.Equation( P[split_2] == Pf )
m.Equation( xH2[split_2] == xf_H2 )
m.Equation( xN2[split_2] == xf_N2 )
m.Equation( xNH3[split_2] == xf_NH3 )
m.Equation( xinert[split_2] == xf_inert )
m.Equation( T[split_2] == Tf )
#SPLITTER 3
m.Equation( n[split_3] == u3 * nf )
m.Equation( P[split_3] == Pf )
m.Equation( xH2[split_3] == xf_H2 )
m.Equation( xN2[split_3] == xf_N2 )
m.Equation( xNH3[split_3] == xf_NH3 )
m.Equation( xinert[split_3] == xf_inert )
m.Equation( T[split_3] == Tf )
#%% HEAT EXCHANGER 2
NTU = m.Intermediate( U*A_HX2 / (n[split_3] * Cp_c) )
Cr = m.Intermediate( n[split_3]*Cp_c / (n[R2_o]*Cp_h) )
Q = m.Intermediate( (n[split_3] * Cp_c * (T[R2_o] - T[split_3]) ) * (1 - m.exp(-NTU * (1-Cr)) ) / ( 1 - Cr * m.exp(-NTU * (1-Cr)) ) )
m.Equation(n[HX2_o] == n[split_3] )
m.Equation(P[HX2_o] == P[split_3] )
m.Equation(xH2[HX2_o] == xH2[split_3] )
m.Equation(xN2[HX2_o] == xN2[split_3] )
m.Equation(xNH3[HX2_o] == xNH3[split_3] )
m.Equation(xinert[HX2_o] == xinert[split_3] )
m.Equation(T[HX2_o] == T[split_3] + (Q / (n[split_3]*Cp_c)) )
#%% MIXER 1
m.Equation(n[R1_i] == n[HX2_o] + n[split_1] )
m.Equation(P[R1_i] == P[split_1] )
m.Equation(xH2[R1_i] == (n[HX2_o]*xH2[HX2_o] + n[split_1]*xH2[split_1]) / n[R1_i] )
m.Equation(xN2[R1_i] == (n[HX2_o]*xN2[HX2_o] + n[split_1]*xN2[split_1]) / n[R1_i] )
m.Equation(xNH3[R1_i] == (n[HX2_o]*xNH3[HX2_o] + n[split_1]*xNH3[split_1]) / n[R1_i] )
m.Equation(xinert[R1_i] == (n[HX2_o]*xinert[HX2_o] + n[split_1]*xinert[split_1]) / n[R1_i] )
m.Equation(T[R1_i] == (n[HX2_o]*T[HX2_o] + n[split_1]*T[split_1]) / n[R1_i])
#%% REACTOR BED 1
#FIRST COMPARTMENT
m.Equation(n_r[0] == n[R1_i] )
m.Equation(P_r[0] == P[R1_i] )
m.Equation(n_r[0]*xH2_r[0] == n[R1_i]*xH2[R1_i] )
m.Equation(n_r[0]*xN2_r[0] == n[R1_i]*xN2[R1_i] )
m.Equation(n_r[0]*xNH3_r[0] == n[R1_i]*xNH3[R1_i] )
m.Equation(n_r[0]*xinert_r[0] == n[R1_i]*xinert[R1_i] )
m.Equation(T_r[0] == T[R1_i] )
#MIDDLE COMPARTMENTS
m.Equations([n_r[i] == (n_r[i-1] + sum(stoi)*mcat_1* (4.75*( Afor*m.exp(-Eafor/(R*T_r[i])) * xN2_r[i]*P_r[i] * ((xH2_r[i]*P_r[i])**1.5) / (xNH3_r[i]*P_r[i]) - Abac*m.exp(-Eabac/(R*(T_r[i]))) * (xNH3_r[i]*P_r[i]) / ((xH2_r[i]*P_r[i])**1.5) )/3600/rhocat)) for i in range(1,ns_R1) ])
m.Equations([P_r[i] == P_r[i-1] for i in range(1,ns_R1)])
m.Equations([n_r[i]*xH2_r[i] == ((n_r[i-1]*xH2_r[i-1]+mcat_1* (4.75*( Afor*m.exp(-Eafor/(R*T_r[i])) * xN2_r[i]*P_r[i] * ((xH2_r[i]*P_r[i])**1.5) / (xNH3_r[i]*P_r[i]) - Abac*m.exp(-Eabac/(R*(T_r[i]))) * (xNH3_r[i]*P_r[i]) / ((xH2_r[i]*P_r[i])**1.5) )/3600/rhocat)*stoi[0])) for i in range(1,ns_R1)])
m.Equations([n_r[i]*xN2_r[i] == ((n_r[i-1]*xN2_r[i-1]+mcat_1* (4.75*( Afor*m.exp(-Eafor/(R*T_r[i])) * xN2_r[i]*P_r[i] * ((xH2_r[i]*P_r[i])**1.5) / (xNH3_r[i]*P_r[i]) - Abac*m.exp(-Eabac/(R*(T_r[i]))) * (xNH3_r[i]*P_r[i]) / ((xH2_r[i]*P_r[i])**1.5) )/3600/rhocat)*stoi[1])) for i in range(1,ns_R1)])
m.Equations([n_r[i]*xNH3_r[i] == ((n_r[i-1]*xNH3_r[i-1]+mcat_1* (4.75*( Afor*m.exp(-Eafor/(R*T_r[i])) * xN2_r[i]*P_r[i] * ((xH2_r[i]*P_r[i])**1.5) / (xNH3_r[i]*P_r[i]) - Abac*m.exp(-Eabac/(R*(T_r[i]))) * (xNH3_r[i]*P_r[i]) / ((xH2_r[i]*P_r[i])**1.5) )/3600/rhocat)*stoi[2])) for i in range(1,ns_R1)])
m.Equations([n_r[i]*xinert_r[i] == (n_r[i-1]*xinert_r[i-1]) for i in range(1,ns_R1)])
m.Equations([(mcat_1*Cpcat)* T_r[i].dt() == (Cp*(n_r[i-1]*T_r[i-1] - n_r[i]*T_r[i]) + (4.75*( Afor*m.exp(-Eafor/(R*T_r[i])) * xN2_r[i]*P_r[i] * ((xH2_r[i]*P_r[i])**1.5) / (xNH3_r[i]*P_r[i]) - Abac*m.exp(-Eabac/(R*(T_r[i]))) * (xNH3_r[i]*P_r[i]) / ((xH2_r[i]*P_r[i])**1.5) )/3600/rhocat)*mcat_1*dHrx) for i in range(1,ns_R1)])
#PASS THE LAST VALUE
m.Equation(n[R1_o] == n_r[ns_R1-1] )
m.Equation(P[R1_o] == P_r[ns_R1-1] )
m.Equation(xH2[R1_o] == xH2_r[ns_R1-1] )
m.Equation(xN2[R1_o] == xN2_r[ns_R1-1] )
m.Equation(xNH3[R1_o] == xNH3_r[ns_R1-1] )
m.Equation(xinert[R1_o] == xinert_r[ns_R1-1] )
m.Equation(T[R1_o] == T_r[ns_R1-1] )
#%% MIXER 2
m.Equation(n[R2_i] == n[R1_o] + n[split_2])
m.Equation(P[R2_i] == P[split_2] )
m.Equation(xH2[R2_i] == (n[R1_o]*xH2[R1_o] + n[split_2]*xH2[split_2]) / n[R2_i] )
m.Equation(xN2[R2_i] == (n[R1_o]*xN2[R1_o] + n[split_2]*xN2[split_2]) / n[R2_i] )
m.Equation(xNH3[R2_i] == (n[R1_o]*xNH3[R1_o] + n[split_2]*xNH3[split_2]) / n[R2_i] )
m.Equation(xinert[R2_i] == (n[R1_o]*xinert[R1_o] + n[split_2]*xinert[split_2]) / n[R2_i] )
m.Equation(T[R2_i] == (n[R1_o]*T[R1_o] + n[split_2]*T[split_2]) / n[R2_i])
#%% REACTOR BED 2
#FIRST COMPARTMENT
m.Equation(n_r[ns_R1] == n[R2_i] )
m.Equation(P_r[ns_R1] == P[R2_i] )
m.Equation(n_r[ns_R1]*xH2_r[ns_R1] == n[R2_i]*xH2[R2_i] )
m.Equation(n_r[ns_R1]*xN2_r[ns_R1] == n[R2_i]*xN2[R2_i] )
m.Equation(n_r[ns_R1]*xNH3_r[ns_R1] == n[R2_i]*xNH3[R2_i] )
m.Equation(n_r[ns_R1]*xinert_r[ns_R1] == n[R2_i]*xinert[R2_i] )
m.Equation(T_r[ns_R1] == T[R2_i])
#MIDDLE COMPARTMENTS
m.Equations([n_r[i] == (n_r[i-1] + sum(stoi)*mcat_2* (4.75*( Afor*m.exp(-Eafor/(R*T_r[i])) * xN2_r[i]*P_r[i] * ((xH2_r[i]*P_r[i])**1.5) / (xNH3_r[i]*P_r[i]) - Abac*m.exp(-Eabac/(R*(T_r[i]))) * (xNH3_r[i]*P_r[i]) / ((xH2_r[i]*P_r[i])**1.5) )/3600/rhocat)) for i in range(ns_R1+1, ns) ])
m.Equations([P_r[i] == P_r[i-1] for i in range(ns_R1+1, ns)])
m.Equations([n_r[i]*xH2_r[i] == ((n_r[i-1]*xH2_r[i-1]+mcat_2* (4.75*( Afor*m.exp(-Eafor/(R*T_r[i])) * xN2_r[i]*P_r[i] * ((xH2_r[i]*P_r[i])**1.5) / (xNH3_r[i]*P_r[i]) - Abac*m.exp(-Eabac/(R*(T_r[i]))) * (xNH3_r[i]*P_r[i]) / ((xH2_r[i]*P_r[i])**1.5) )/3600/rhocat)*stoi[0])) for i in range(ns_R1+1, ns)])
m.Equations([n_r[i]*xN2_r[i] == ((n_r[i-1]*xN2_r[i-1]+mcat_2* (4.75*( Afor*m.exp(-Eafor/(R*T_r[i])) * xN2_r[i]*P_r[i] * ((xH2_r[i]*P_r[i])**1.5) / (xNH3_r[i]*P_r[i]) - Abac*m.exp(-Eabac/(R*(T_r[i]))) * (xNH3_r[i]*P_r[i]) / ((xH2_r[i]*P_r[i])**1.5) )/3600/rhocat)*stoi[1])) for i in range(ns_R1+1, ns)])
m.Equations([n_r[i]*xNH3_r[i] == ((n_r[i-1]*xNH3_r[i-1]+mcat_2* (4.75*( Afor*m.exp(-Eafor/(R*T_r[i])) * xN2_r[i]*P_r[i] * ((xH2_r[i]*P_r[i])**1.5) / (xNH3_r[i]*P_r[i]) - Abac*m.exp(-Eabac/(R*(T_r[i]))) * (xNH3_r[i]*P_r[i]) / ((xH2_r[i]*P_r[i])**1.5) )/3600/rhocat)*stoi[2])) for i in range(ns_R1+1, ns)])
m.Equations([n_r[i]*xinert_r[i] == (n_r[i-1]*xinert_r[i-1]) for i in range(ns_R1+1, ns)])
m.Equations([(mcat_2*Cpcat)* T_r[i].dt() == (Cp*(n_r[i-1]*T_r[i-1] - n_r[i]*T_r[i]) + (4.75*( Afor*m.exp(-Eafor/(R*T_r[i])) * xN2_r[i]*P_r[i] * ((xH2_r[i]*P_r[i])**1.5) / (xNH3_r[i]*P_r[i]) - Abac*m.exp(-Eabac/(R*(T_r[i]))) * (xNH3_r[i]*P_r[i]) / ((xH2_r[i]*P_r[i])**1.5) )/3600/rhocat)*mcat_2*dHrx) for i in range(ns_R1+1, ns)])
#PASS THE LAST VALUE
m.Equation(n[R2_o] == n_r[ns-1])
m.Equation(P[R2_o] == P_r[ns-1] )
m.Equation(xH2[R2_o] == xH2_r[ns-1] )
m.Equation(xN2[R2_o] == xN2_r[ns-1] )
m.Equation(xNH3[R2_o] == xNH3_r[ns-1] )
m.Equation(xinert[R2_o] == xinert_r[ns-1] )
m.Equation(T[R2_o] == T_r[ns-1])
#%% SOLVE DYNAMIC MODEL
m.options.IMODE = 4
m.solve(debug=0)
#PLOT MY GEKKO RESULTS
plt.figure()
t = m.time
plt.plot(t/60, T_r[-1])
plt.xlabel("Time [minutes]")
plt.ylabel("Temperature [K]")
GEKKO_single is missing so the script is missing the configuration parameters. Here are some tips in absence of any verification:
nt = 60 #seconds
m.time = np.linspace(0,nt,nt) #seconds
COLDSTART=2:m.options.COLDSTART=2
m.solve()
This can take longer to solve but generally does a better job of finding an initial solution.
IMODE=1:m.options.IMODE=1
m.solve()
m.options.IMODE=4
m.solve()
This should also work but you may need to transfer values from the steady state solution if you previously set the value property for variables when you defined the problem.
If you can post the GEKKO_single configuration parameters then we can give more specific feedback by verifying the methods.