I am trying to solve nine nonlinear equations with nine variables using fsolve in python. But I don't know how to set the constraints to the variables. Besides, the iteration of fsolve is not making good progress with the current code.
from scipy.optimize import fsolve
import numpy as np
KH = 3.39E-02
Kw = 1.00E-14
KC1 = 4.47E-07
KC2 = 4.69E-11
KP1 = 6.92E-03
KP2 = 6.17E-08
KP3 = 4.79E-13
K = 1
Pt = 0.1
pCO2 = 0.0004
def equations(p):
q1, q2, q3, q4, q5, q6, q7, q8, q9 = p
e1 = q1/pCO2-KH
e2 = q2*q3-Kw
e3 = q2*q4/q1-KC1
e4 = q2*q5/q4-KC2
e5 = q2*q7/q6-KP1
e6 = q2*q8/q7-KP2
e7 = q2*q9/q8-KP3
e8 = q3+q4+2*q5+q7+2*q8+3*q9-q2-K
e9 = q6+q7+q8+q9-Pt
return [e1, e2, e3, e4, e5, e6, e7, e8, e9]
root = fsolve(equations, (1, 1, 1, 1, 1, 1, 1, 1, 1))
print(root)
All the equations are actually equal to 0, but apparently fsolve could not take 0 as the first guess. The constraint is that all the variables should be positive. How could I find the solutions?
Thanks for the suggestion of using scipy.optimize.minimize.
from scipy.optimize import minimize
import numpy as np
KH = 3.39E-02
Kw = 1.00E-14
KC1 = 4.47E-07
KC2 = 4.69E-11
KP1 = 6.92E-03
KP2 = 6.17E-08
KP3 = 4.79E-13
K = 1
Pt = 0.1
pCO2 = 0.0004
def equations(p):
q1, q2, q3, q4, q5, q6, q7, q8, q9 = p
E = np.empty((9))
E[0] = q1/pCO2-KH
E[1] = q2*q3-Kw
E[2] = q2*q4/q1-KC1
E[3] = q2*q5/q4-KC2
E[4] = q2*q7/q6-KP1
E[5] = q2*q8/q7-KP2
E[6] = q2*q9/q8-KP3
E[7] = q3+q4+2*q5+q7+2*q8+3*q9-q2-K
E[8] = q6+q7+q8+q9-Pt
return E
#initial point
pGuess = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
#bounds q1, q2, q3, q4, q5, q6, q7, q8, q9 >= 0
bounds = [(0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None)]
# Solve the constrained optimization problem
p = minimize(lambda p: np.linalg.norm(equations(p)), x0=pGuess, bounds=bounds)
print(p.x)
But I encountered dividing by 0 during the calculation.
