Solving non linear equations: add constraints to Gibbs free energy problem

Question

I'm trying to solve a non linear system that will minimize the Gibbs free energy using the Lagrange method, with the exponential formulation. The equations already have the Lagrangians in the exponential form Y1...Y6 that are later converted to the mole number of the chemical species n1...n9.

The problem is that fsolve() gives answers that vary a lot, even if rerunning the problem with the same guess, it gives different values. But the main problem is that all the solutions that I reach with different guesses don't make physical sense because after converting Ys to ns I get negative values for mass.

So by the physics involved I can determine that all [n1...n9] >= 0. Also can determine all the maximum values for the [n1...n9].

How can I add this to the code?

import numpy as np
import scipy
from scipy.optimize import fsolve
import time
#
# "B" is the energy potentials of the species [C_gr , CO , CO2 , H2 , CH4 , H2O , N2* , SiO2* , H2S]
B = [-11.0, -309.3632404425132, -613.3667287153355, -135.61840658777166, -269.52018727412405, -434.67499662354476, -193.0773646004259, -980.0, -230.02942769438977]
# "a_atoms" is the number of atoms in the reactants [C, H, O, N*, S, SiO2*]  
# * Elements that doesn't react. '
a_atoms = [4.27311296e-02, 8.10688756e-02, 6.17738749e-02, 1.32864225e-01, 3.18931655e-05, 3.74477901e-04]
P_zero = 100.0 # Standard energy pressure
P_eq = 95.0 # Reaction pressure
# Standard temperature 298.15K, reaction temperature 940K.
#
start_time = time.time()
def GibbsEq(z):
# Lambda's exponentials:
    Y1 = z[0]
    Y2 = z[1] 
    Y3 = z[2]
    Y4 = z[3] 
    Y5 = z[4] 
    Y6 = z[5]
# Number of moles in each phase:
    N1 = z[6]
    N2 = z[7]
    N3 = z[8]
# Equations of energy conservation and mass conservation:
    F = np.zeros(9) 
    F[0] = (P_zero/P_eq) * N1 * ((B[1] * (Y1 * Y3) + B[2] * (Y1 * Y3**2) + B[4] * (Y1 * Y2**2)) + N2 * (B[0] * Y1)) - a_atoms[0]
    F[1] = (P_zero/P_eq) * N1 * (2 * B[3] * Y2**2 + 4 * B[4] * (Y1 * Y2**4) + 2 * B[5] * ((Y2**2) * Y3) + 2 * B[8] * ((Y2**2) * Y5)) - a_atoms[1]
    F[2] = (P_zero/P_eq) * N1 * (B[1] * (Y1 * Y3) + 2 * B[2] * (Y1 * Y3**2) + B[5] * ((Y2**2) * Y3)) - a_atoms[2]
    F[3] = (P_zero/P_eq) * N1 * (2 * B[6]**2) - a_atoms[3]
    F[4] = (P_zero/P_eq) * N1 * (B[8] * ((Y2**2) * Y5)) - a_atoms[4]
    F[5] = N3 * (B[7] * Y5)  - a_atoms[5]
# 
    F[6] = (P_zero/P_eq) * (B[1] * (Y1 * Y3) + B[2] * (Y1 * Y3**2) + B[3] * Y2**2 + B[4] * (Y1 * Y2**4) + B[5] * ((Y2**2) * Y3) + B[6] * Y4 + B[8] * Y5) - 1 
    F[7] = B[0] * Y1 - 1 
    F[8] = B[7] * Y6 - 1
    return F
#
zGuess = np.ones(9)
z = scipy.optimize.fsolve(GibbsEq, zGuess)
end_time = time.time()
time_solution = (end_time - start_time)
print('Solving time: {} s'.format(time_solution))
#
n1 = z[7] * B[0] * z[0]
n2 = z[6] * B[1] * z[0] * z[2]
n3 = z[6] * B[2] * z[0] * z[2]**2
n4 = z[6] * B[3] * z[1]**2
n5 = z[6] * B[4] * z[0] * z[1]**4
n6 = z[6] * B[5] * z[1]**2 * z[4]
n7 = z[6] * B[6] * z[3]**2
n8 = z[8] * B[7] * z[5]
n9 = z[6] * B[8] * z[1]**2 * z[4]
N_T = [n1, n2, n3, n4, n5, n6, n7, n8, n9]
print(z)
print(z[6],z[7],z[8])
print(N_T)
for n in N_T:
    if n < 0:
        print('Error: there is negative values for mass in the solution!')
        break

How can I add constraints in fsolve?
Are there other solvers in python that have more options of constrains to get stability and more independence of initial guess?

Thank you!

Siva-Sg · Accepted Answer · 2019-04-10T12:59:17.553

One answer for both the questions.

fsolve does not support constraints. You may provide initial estimate as positive values , but that does not guarantee positive roots . However, you can reformulate your problem as optimization problem and minimize the cost function imposing constraints using any optimization function such as scipy.optimize.minimize.

As a minimal example , if you want to find the positive root of the equation x*x -4 , you could do as below.

scipy.optimize.minimize(lambda x:(x*x-4)**2,x0= [5], bounds =((0,None),))

The bounds parameter which takes (min,max) pair can be used to impose the positive constraint on the root.

output :

 fun: array([1.66882981e-17])
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
      jac: array([1.27318954e-07])
  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 20
      nit: 9
   status: 0
  success: True
        x: array([2.])

Going by this, your code can be modified as below. Just add the bounds , change the function return statement , and change the fsolve to scipy.optimize.minimize with the bounds.

import numpy as np
import scipy
from scipy.optimize import fsolve
import time
#
# "B" is the energy potentials of the species [C_gr , CO , CO2 , H2 , CH4 , H2O , N2* , SiO2* , H2S]
B = [-11.0, -309.3632404425132, -613.3667287153355, -135.61840658777166, -269.52018727412405, -434.67499662354476, -193.0773646004259, -980.0, -230.02942769438977]
# "a_atoms" is the number of atoms in the reactants [C, H, O, N*, S, SiO2*]  
# * Elements that doesn't react. '
a_atoms = [4.27311296e-02, 8.10688756e-02, 6.17738749e-02, 1.32864225e-01, 3.18931655e-05, 3.74477901e-04]
P_zero = 100.0 # Standard energy pressure
P_eq = 95.0 # Reaction pressure
# Standard temperature 298.15K, reaction temperature 940K.
#
start_time = time.time()
def GibbsEq(z):
# Lambda's exponentials:
    Y1 = z[0]
    Y2 = z[1] 
    Y3 = z[2]
    Y4 = z[3] 
    Y5 = z[4] 
    Y6 = z[5]
# Number of moles in each phase:
    N1 = z[6]
    N2 = z[7]
    N3 = z[8]

    bounds =((0,None),)*9
# Equations of energy conservation and mass conservation:
    F = np.zeros(9) 
    F[0] = (P_zero/P_eq) * N1 * ((B[1] * (Y1 * Y3) + B[2] * (Y1 * Y3**2) + B[4] * (Y1 * Y2**2)) + N2 * (B[0] * Y1)) - a_atoms[0]
    F[1] = (P_zero/P_eq) * N1 * (2 * B[3] * Y2**2 + 4 * B[4] * (Y1 * Y2**4) + 2 * B[5] * ((Y2**2) * Y3) + 2 * B[8] * ((Y2**2) * Y5)) - a_atoms[1]
    F[2] = (P_zero/P_eq) * N1 * (B[1] * (Y1 * Y3) + 2 * B[2] * (Y1 * Y3**2) + B[5] * ((Y2**2) * Y3)) - a_atoms[2]
    F[3] = (P_zero/P_eq) * N1 * (2 * B[6]**2) - a_atoms[3]
    F[4] = (P_zero/P_eq) * N1 * (B[8] * ((Y2**2) * Y5)) - a_atoms[4]
    F[5] = N3 * (B[7] * Y5)  - a_atoms[5]
# 
    F[6] = (P_zero/P_eq) * (B[1] * (Y1 * Y3) + B[2] * (Y1 * Y3**2) + B[3] * Y2**2 + B[4] * (Y1 * Y2**4) + B[5] * ((Y2**2) * Y3) + B[6] * Y4 + B[8] * Y5) - 1 
    F[7] = B[0] * Y1 - 1 
    F[8] = B[7] * Y6 - 1
    return (np.sum(F)**2)
#
zGuess = np.ones(9)
z = scipy.optimize.minimize(GibbsEq, zGuess , bounds=bounds)
end_time = time.time()
time_solution = (end_time - start_time)
print('Solving time: {} s'.format(time_solution))
#

print(z.x)

print(N_T)
for n in N_T:
    if n < 0:
        print('Error: there is negative values for mass in the solution!')
        break

output:

Solving time: 0.012451648712158203 s
[1.47559173 2.09905553 1.71722403 1.01828262 1.17529548 1.08815712
 1.00294916 1.00104157 1.08815763]

Solving non linear equations: add constraints to Gibbs free energy problem

1 Answers1