Solution of a nonlinear system highly dependent on initial guess

Question

I have a system of nonlinear equations

1. xF - xA * exp(-(wAA * xA + wBB * xB)) = 0
2. xF - xB * exp(-(wBB * xB + wAB * xA)) = 0
3. xA + xB + xF -1 = 0

where the variables are xA, xB and xF, while wAA, wBB and wAB are given parameters. The variables are non negative and in the interval [0,1] Here there is a minimal example code

import math
from scipy.optimize import least_squares

coupleInteractions = [10000., 10., 0.]

def equations(p):
    xA, xB, xF = p
    wAA, wBB, wAB = coupleInteractions

    eq1 = xF - xA * math.exp(-(wAA * xA + wBB * xB))
    eq2 = xF - xB * math.exp(-(wBB * xB + wAB * xA))
    eq3 = xA + xB + xF -1

    return (eq1, eq2, eq3)

guess = [0., 1., 0.]
bounds = ([0., 0., 0.], [1., 1., 1.])
root = least_squares(equations, x0=guess, bounds=bounds, xtol=1.e-15, gtol = 1.e-15, ftol = 1.e-15, loss="cauchy")

print(root)
print(root.x)

I want to test the code in some limit cases. For example if all the parameters are zero, wAA = wAB = wBB = 0, I have a simple linear system, with xA = xB = xF = 1./3. In this case the code works perfectly without dependence on the guess values.

I want to test, like in the code posted, the case wAA >> wBB and wAA >> wAB. From a mathematical point of view I obtain xA = 1, xB = 0 and xF = 0. But with the least_square function I obtain the wrong answer. I deliberately use a "wrong" guess to test the method, because in the intermediate cases I do not want that the results higly depends on the initial choice. Is there an error in my code?

I have performed a new test to see if a global optimization could work. I have defined the function:

def scalarEquations(p):
    xA, xB, xF = p
    wAA, wBB, wAB = coupleInteractions

    eq1 = xF - xA * math.exp(-(wAA * xA + wBB * xB))
    eq2 = xF - xB * math.exp(-(wBB * xB + wAB * xA))
    eq3 = xA + xB + xF -1

    return (eq1**2 + eq2**2 + eq3 **2)

It returns the sum of the squares for the 3 equations. Then I performed a brute force approach.

from scipy.optimize import brute
rranges = (slice(0,1, 0.05), slice(0,1, 0.005), slice(0,1, 0.005))
t = brute(scalarEquations, rranges)
print(t)

I obtain the same "wrong" solution. I know that a global optimizer should find the absolute minimum. I do not understand.

Answer 1

Your code looks ok. Looks like there is a local minimum near [0,1,0] for wAA large and wBB greater than e and the optimization is converging to it.