Question: 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:    Solution Vector x GEKKO EXCEL GRG Root finder GEKKO-final     H2: 1 5.797458326 5.344360712 5.118105 5.219   H2O: 2 1.202541674 1.521995663 1.387357 1.681   CO: 3 1.202541674 1.388351672 1.685236 1.581   O2: 4 1.6791020526e-21 5.82E-21 0.4114448 1e-5   CO2: 5 0.79745832603 0.544826 0.0522585 0.3693   CH4: 6 1.6503662452e-22 0.0668213 0.2269213 0.05   C2H6: 7 3.2694282596e-27 1.70E-07 0.0092403 1e-5   C2H4: 8 1.1255712085e-27 9.74E-08 0.0044224 1e-5   C2H2: 9 1e-27 3.25E-10 0.0041294 1e-5   lamda1: 10 2.185089969 24.3878385 1.537016    lamda2: 11 2.185089969 0.253110984 0.2765838    lamda3: 12 2.185089969 1.558979624 2.53957    summe: 13 9.5605307091 8.866356107 8.899115 8.90   

Answer: 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