使用 GEKKO IPOPT 对氨反应器进行动态建模
Dynamic modelling of ammonia reactor using GEKKO IPOPT
我正在尝试使用 GEKKO 对氨反应器进行动态模拟。不幸的是，我的代码错误表明它无法找到解决方案。我试图求解的方程在物种连续性方程方面是稳态的，但在热平衡方面是动态的。我使用下面显示的代码首先解决了稳态条件，它通过设置 IMODE = 1 成功地解决了这个问题。但是，当我导入稳态结果并用它来初始化我的动态代码时，没有解决方案已完成。我不明白我的代码在求解动态模型时为什么不起作用，因为它对稳态模型运行良好。我已经尝试过 IMODE = 4 和 7，但它们都不起作用。有谁知道为什么我的代码不适用于动态模拟？
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 丢失，因此脚本缺少配置参数。以下是在没有任何验证的情况下的一些提示：
- 设置更小的时间范围，看看它是否会完成一个解决方案：
nt = 60 #seconds
m.time = np.linspace(0,nt,nt) #seconds
- 初始化
COLDSTART=2:
m.options.COLDSTART=2
m.solve()
这可能需要更长的时间来解决，但通常可以更好地找到初始解决方案。
- 初始化
IMODE=1:
m.options.IMODE=1
m.solve()
m.options.IMODE=4
m.solve()
这也应该有效，但如果您之前在定义问题时为变量设置了 value 属性，则可能需要从稳态解中转移值。
如果你可以post配置参数GEKKO_single那么我们可以通过验证方法给出更具体的反馈。
我正在尝试使用 GEKKO 对氨反应器进行动态模拟。不幸的是，我的代码错误表明它无法找到解决方案。我试图求解的方程在物种连续性方程方面是稳态的，但在热平衡方面是动态的。我使用下面显示的代码首先解决了稳态条件，它通过设置 IMODE = 1 成功地解决了这个问题。但是，当我导入稳态结果并用它来初始化我的动态代码时，没有解决方案已完成。我不明白我的代码在求解动态模型时为什么不起作用，因为它对稳态模型运行良好。我已经尝试过 IMODE = 4 和 7，但它们都不起作用。有谁知道为什么我的代码不适用于动态模拟？
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 丢失，因此脚本缺少配置参数。以下是在没有任何验证的情况下的一些提示：
- 设置更小的时间范围，看看它是否会完成一个解决方案：
nt = 60 #seconds
m.time = np.linspace(0,nt,nt) #seconds
- 初始化
COLDSTART=2:
m.options.COLDSTART=2
m.solve()
这可能需要更长的时间来解决，但通常可以更好地找到初始解决方案。
- 初始化
IMODE=1:
m.options.IMODE=1
m.solve()
m.options.IMODE=4
m.solve()
这也应该有效，但如果您之前在定义问题时为变量设置了 value 属性，则可能需要从稳态解中转移值。
如果你可以post配置参数GEKKO_single那么我们可以通过验证方法给出更具体的反馈。