IPOPT solution by Gekko in comparison to GRG algorithm used within the Excel-Solver

The aim is to compute the thermodynamically equilibrated composition of the mixture at 1000K based on Gibbs' energy of formation of the reaction products and educts (steam+C2H6 in a molar ratio of 4:1) for the ethane steam gasification reaction as following:

from gekko import GEKKO     
m = GEKKO(remote=True)
x = m.Array(m.Var,13,value=1,lb=1e-27,ub=20.0)
H2,H2O,CO,O2,CO2,CH4,C2H6,C2H4,C2H2,lamda1,lamda2,lamda3,summe= x
H2.value = 3.0
H2O.value = 1.0
CO.value = 0.5
O2.value = 0.001
CO2.value = 0.5
CH4.value = 0.1
C2H6.value = 0.000215
C2H4.value = 0.00440125
C2H2.value = 0.0041294
summe.value = 8.0
lamda1.value = 1.0
lamda2.value = 1.0
lamda3.value = 1.0
eq1 = m.Param(value=14)
eq2 = m.Param(value=4)
eq3 = m.Param(value=2)
summe = m.Var(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2)
lamda2 = m.Var((-1)*m.log(H2 / summe) / 2)
lamda1 = m.Var(46.03 / 1.9872 - m.log(H2O / summe) + 2 * lamda2)
lamda3 = m.Var(47.942 / 1.9872 - m.log(CO / summe) + lamda1)
m.Equation(m.exp(-4.61 / 1.9872 - 4 * lamda2 - lamda3) * summe == CH4)
m.Equation(m.exp(-28.249 / 1.9872 - 4 * lamda2 - 2 * lamda3) * summe == C2H4)
m.Equation(m.exp(-40.604 / 1.9872 - 2 * lamda2 - 2 * lamda3) * summe == C2H2)
m.Equation(m.exp(-26.13 / 1.9872 - 6 * lamda2 - 2 * lamda3) * summe == C2H6)
m.Equation(m.exp(94.61 / 1.9872 - 2 * lamda1 - lamda3) * summe == CO2)
m.Equation(m.exp(-2 * lamda1) * summe == O2)
m.Equation(2 * CO2 + CO + 2 * O2 + H2O  == eq2)
m.Equation(4 * CH4 + 4 * C2H4 + 2 * C2H2 + 2 * H2 + 2 * H2O + 6 * C2H6  == eq1)
m.Equation(CH4 + 2 * C2H4 + 2 * C2H2 + CO2 + CO + 2 * C2H6  == eq3)
m.Minimize((summe-(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2))**2)
m.options.IMODE = 3 #IPOPT
m.options.MAX_ITER = 1000
m.options.OTOL = 1e-10
m.options.RTOL = 1e-10
m.solve()
print('x: ', x)
print('Objective: ',m.options.OBJFCNVAL)

---

EXIT: Optimal Solution Found. 
The solution was found. 
The final value of the objective function is    81.0000000000000     

---------------------------------------------------
Solver         :  IPOPT (v3.12)
Solution time  :   1.169999998819549E-002 sec
Objective      :    81.0000000000000     
Successful solution
---------------------------------------------------

x:  [[5.797458326] [1.202541674] [1.202541674] [1.6791020526e-21]
[0.79745832603] [1.6503662452e-22] [3.2694282596e-27] [1.1255712085e-27]
[1e-27] [2.185089969] [2.185089969] [2.185089969] [9.5605307091]]

Then using GRG as optimiser:

Both solution still differs. And interestingly, the output from a (constraint) root finder is still different, but converges at the same time:

Total number of equations:  13  
Number of implicit equations:  4  
Number of explicit equations:  9  
Solution method  CONSTRAINED  
Convergence tolerance:  1e-07  
# of iterations used:  8  

CO  1.685236
H2  5.118105
H2O  1.387357
SUM  8.899115  
C2H2  0.0041294  
C2H4  0.0044224  
C2H6  0.0092403  
CH4  0.2269213  
CO2  0.0522585  
lamda1  1.537016  
lamda2  0.2765838  
lamda3  2.53957  
O2  0.4114448  

Here is a comparison of the solutions:

The variables lamda1 ,lamda2, lamda3, summe are defined twice and the equations associated with those variables are only used to initialize the second definition of the variables.

H2,H2O,CO,O2,CO2,CH4,C2H6,C2H4,C2H2,lamda1,lamda2,lamda3,summe= x
summe = m.Var(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2)
lamda2 = m.Var((-1)*m.log(H2 / summe) / 2)
lamda1 = m.Var(46.03 / 1.9872 - m.log(H2O / summe) + 2 * lamda2)
lamda3 = m.Var(47.942 / 1.9872 - m.log(CO / summe) + lamda1)

Switching them to Intermediate() definitions is an easy way to fix the problem.

summe = m.Intermediate(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2)
lamda2 = m.Intermediate((-1)*m.log(H2 / summe) / 2)
lamda1 = m.Intermediate(46.03 / 1.9872 - m.log(H2O / summe) + 2 * lamda2)
lamda3 = m.Intermediate(47.942 / 1.9872 - m.log(CO / summe) + lamda1)

Here is the complete script:

from gekko import GEKKO   
import numpy as np  
m = GEKKO(remote=True)
x = m.Array(m.Var,9,value=1,lb=1e-27,ub=20.0)
H2,H2O,CO,O2,CO2,CH4,C2H6,C2H4,C2H2=x
H2.value = 3.0
H2O.value = 1.0
CO.value = 0.5
O2.value = 0.001
CO2.value = 0.5
CH4.value = 0.1
C2H6.value = 0.000215
C2H4.value = 0.00440125
C2H2.value = 0.0041294
eq1 = m.Param(value=14)
eq2 = m.Param(value=4)
eq3 = m.Param(value=2)
summe = m.Intermediate(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2)
lamda2 = m.Intermediate((-1)*m.log(H2 / summe) / 2)
lamda1 = m.Intermediate(46.03 / 1.9872 - m.log(H2O / summe) + 2 * lamda2)
lamda3 = m.Intermediate(47.942 / 1.9872 - m.log(CO / summe) + lamda1)
m.Equation(m.exp(-4.61 / 1.9872 - 4 * lamda2 - lamda3) * summe == CH4)
m.Equation(m.exp(-28.249 / 1.9872 - 4 * lamda2 - 2 * lamda3) * summe == C2H4)
m.Equation(m.exp(-40.604 / 1.9872 - 2 * lamda2 - 2 * lamda3) * summe == C2H2)
m.Equation(m.exp(-26.13 / 1.9872 - 6 * lamda2 - 2 * lamda3) * summe == C2H6)
m.Equation(m.exp(94.61 / 1.9872 - 2 * lamda1 - lamda3) * summe == CO2)
m.Equation(m.exp(-2 * lamda1) * summe == O2)
m.Equation(2 * CO2 + CO + 2 * O2 + H2O  == eq2)
m.Equation(4 * CH4 + 4 * C2H4 + 2 * C2H2 + 2 * H2 + 2 * H2O + 6 * C2H6  == eq1)
m.Equation(CH4 + 2 * C2H4 + 2 * C2H2 + CO2 + CO + 2 * C2H6  == eq3)
m.Minimize((summe-(H2 + O2 + H2O + CO + CO2 + CH4 + C2H6 + C2H4 + C2H2))**2)
m.options.IMODE = 3 #IPOPT
m.options.MAX_ITER = 1000
m.options.OTOL = 1e-10
m.options.RTOL = 1e-10
m.solve()
print('x: ', x)
print('Objective: ',m.options.OBJFCNVAL)
print(f'H2: {H2.value[0]}')
print(f'H2O: {H2O.value[0]}')
print(f'CO: {CO.value[0]}')
print(f'O2: {O2.value[0]}')
print(f'CO2: {CO2.value[0]}')
print(f'CH4: {CH4.value[0]}')
print(f'C2H6: {C2H6.value[0]}')
print(f'C2H4: {C2H4.value[0]}')
print(f'C2H2: {C2H2.value[0]}')

The objective value is now 0 with 0 degrees of freedom as required for root finding. The solution is:

 The solution was found.
 
 The final value of the objective function is   0.000000000000000E+000
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   1.099999999860302E-002 sec
 Objective      :   0.000000000000000E+000
 Successful solution
 ---------------------------------------------------
 
x:  [[5.0] [2.0] [2.0] [1.0421441742e-21] [3.9704669403e-23]
 [2.1920286233e-23] [1e-27] [1e-27] [1e-27]]
Objective:  0.0
H2: 5.0
H2O: 2.0
CO: 2.0
O2: 1.0421441742e-21
CO2: 3.9704669403e-23
CH4: 2.1920286233e-23
C2H6: 1e-27
C2H4: 1e-27
C2H2: 1e-27

This solution doesn't look quite right for the water-gas shift reaction. I had to go back and refresh my memory for the equilibrium composition of a reacting mixture containing ethane (C₂H₆) and steam (H₂O) at 1000 K.

---

1. Steam Reforming Reaction

The primary reaction for hydrogen production:

C2H6 + 2H2O → 2CO + 5H2

• Purpose: Converts hydrocarbons into carbon monoxide (CO) and hydrogen (H₂) in the presence of steam.
• Thermodynamic Effect: Highly endothermic, favored at high temperatures.

---

2. Water-Gas Shift Reaction

A secondary reaction between carbon monoxide and steam:

CO + H2O ↔ CO2 + H2

• Purpose: Converts CO to CO₂ and produces additional H₂.
• Thermodynamic Effect: Exothermic, favored at lower temperatures but occurs at equilibrium even at high temperatures.

---

3. Other Hydrocarbon Reactions

The system also includes the decomposition and partial oxidation of hydrocarbons:

• Ethylene Formation (C₂H₄):

  C2H6 → C2H4 + H2
  

• Acetylene Formation (C₂H₂):

  C2H4 → C2H2 + H2
  

---

4. Carbon Deposition

Carbon deposition can occur under certain conditions:

CO → C (solid) + CO2

or

CH4 → C (solid) + 2H2

• Purpose: Represents undesirable side reactions where carbon deposits as a solid.
• Thermodynamic Effect: Can occur at lower hydrogen-to-carbon or oxygen-to-carbon ratios, depending on system conditions.

---

Thermodynamic Modeling

The equilibrium composition is calculated by minimizing the Gibbs free energy of the system:

G = Σ (ni * (gi° + RT * ln(ni / ntotal)))

Where:

• ni: Molar amount of species i.
• gi°: Standard Gibbs free energy of formation for species i.
• ntotal: Total moles of all species.

Constraints:

• Mass balances for carbon, hydrogen, and oxygen.
• Bounds on mole amounts to avoid non-physical solutions.

Here is a complete script based on these equations:

from gekko import GEKKO

# Create GEKKO model
m = GEKKO(remote=True)

# Define variables for species molar amounts
n_C2H6 = m.Var(value=1.0, lb=0, ub=1)  # Ethane
n_H2O = m.Var(value=4.0, lb=0, ub=1)   # Water
n_CO = m.Var(value=0.0, lb=0, ub=1)    # Carbon monoxide
n_CO2 = m.Var(value=0.0, lb=0, ub=1)   # Carbon dioxide
n_H2 = m.Var(value=0.0, lb=0, ub=1)    # Hydrogen
n_CH4 = m.Var(value=0.0, lb=0, ub=1)   # Methane
n_C = m.Var(value=0.0, lb=0, ub=1)     # Solid carbon
n_C2H4 = m.Var(value=0.0, lb=0, ub=1)  # Ethylene
n_C2H2 = m.Var(value=0.0, lb=0, ub=1)  # Acetylene

# Universal gas constant
R = 1.9872 # cal/(K mol)

# Gibbs free energies of formation at 1000K
gibbs_energies = {
    "C2H6": -26.13,    # Ethane
    "H2O": 46.03,      # Water (vapor)
    "CO": 47.942,      # Carbon monoxide
    "CO2": 94.61,      # Carbon dioxide
    "H2": 0.0,         # Hydrogen gas
    "CH4": -4.61,      # Methane
    "C": 0.0,          # Solid carbon
    "C2H4": -28.249,   # Ethylene
    "C2H2": -40.604    # Acetylene
}


# Define total moles
total_moles = (
    n_C2H6 + n_H2O + n_CO + n_CO2 + n_H2 + n_CH4 + n_C + n_C2H4 + n_C2H2
)

# Define the objective function (Gibbs free energy minimization)
m.Obj(
    n_C2H6 * (gibbs_energies["C2H6"]/R + m.log(n_C2H6 / total_moles + 1e-10)) +
    n_H2O * (gibbs_energies["H2O"]/R + m.log(n_H2O / total_moles + 1e-10)) +
    n_CO * (gibbs_energies["CO"]/R + m.log(n_CO / total_moles + 1e-10)) +
    n_CO2 * (gibbs_energies["CO2"]/R + m.log(n_CO2 / total_moles + 1e-10)) +
    n_H2 * (gibbs_energies["H2"]/R + m.log(n_H2 / total_moles + 1e-10)) +
    n_CH4 * (gibbs_energies["CH4"]/R + m.log(n_CH4 / total_moles + 1e-10)) +
    n_C * (gibbs_energies["C"]/R + m.log(n_C / total_moles + 1e-10)) +
    n_C2H4 * (gibbs_energies["C2H4"]/R + m.log(n_C2H4 / total_moles + 1e-10)) +
    n_C2H2 * (gibbs_energies["C2H2"]/R + m.log(n_C2H2 / total_moles + 1e-10))
)

# Mass balance constraints
# Carbon balance: 4(C2H6) = CO + CO2 + CH4 + C + C2H4 + C2H2
m.Equation(4*n_C2H6 == n_CO + n_CO2 + n_CH4 + n_C + n_C2H4 + n_C2H2)

# Oxygen balance: 3*H2O = CO + CO2
m.Equation(3*n_H2O == n_CO + n_CO2)

# Hydrogen balance with 4:1 ratio: C2H6 + 4*H2O = 2*H2 + CH4 + C2H4 + C2H2
m.Equation(4*n_H2O + n_C2H6 == 2*n_H2 + n_CH4 + n_C2H4 + n_C2H2)

# Total mole constraint for mole fraction result
m.Equation(total_moles == 1)

# Solve the model
m.solve(disp=True)

# Extract and display results
results = {
    "C2H6 (mol)": n_C2H6.value[0],
    "H2O (mol)": n_H2O.value[0],
    "CO (mol)": n_CO.value[0],
    "CO2 (mol)": n_CO2.value[0],
    "H2 (mol)": n_H2.value[0],
    "CH4 (mol)": n_CH4.value[0],
    "C (mol)": n_C.value[0],
    "C2H4 (mol)": n_C2H4.value[0],
    "C2H2 (mol)": n_C2H2.value[0]
}
[print(f'{r}: {results[r]:.4}') for r in results]

with results:

C2H6 (mol): 0.1993
H2O (mol): 0.003599
CO (mol): 0.0108
CO2 (mol): 0.0
H2 (mol): 0.0
CH4 (mol): 2.915e-09
C (mol): 0.5726
C2H4 (mol): 0.0004254
C2H2 (mol): 0.2133