Trying to solve a system of non-linear equations using GEKKO. The system has multiple solutions, however GEKKO gives the wrong one. How do I fix this?

I am trying to solve systems of multiple non-linear equations similar to this one:

As you can see by the graph, there are two solutions to the system. Through much trial-and-error and different methods (fsolve, etc.) to accomplish this, I stumbled upon GEKKO, and thus the following:

from gekko import GEKKO
m = GEKKO()
x,y = [m.Var(1) for i in range(2)]
m.Equations([0.2618 == 0.2094 + x*(1-(1/y))**(1/y),\
             0.9200 == 1-math.e**((-((.5236-.2094)/x)**y))])
m.solve(disp=False)
print(x.value,y.value)

This code successfully solves the system, however it outputs the incorrect solution. I can't figure out quite how to get the other one (which is the one I need). Any GEKKO experts that can help? Thanks!

I've tried multiple different methods of solving non-linear systems (fsolve, Sympy nsolve, etc) and consistently ran into a division-by-zero error that I was unable to work around. Other methods of solving systems gave no solution. I cannot find the proper documentation for GEKKO to figure out how to get the other solution.

Adding bounds for the search region can help locate a particular solution.

from gekko import GEKKO
m = GEKKO(remote=False)
x,y = m.Array(m.Var,2,value=1)
m.Equations([0.2618 == 0.2094 + x*(1-(1/y))**(1/y),\
             0.9200 == 1-m.exp((-((.5236-.2094)/x)**y))])

y.LOWER = 0
m.solve(disp=False)
print(f'Solution: x={x.value[0]:0.2f} y={y.value[0]:0.2f}')

x.UPPER = 10; x.LOWER = 0
y.UPPER = 0; y.LOWER = -1
m.solve(disp=False)
print(f'Solution: x={x.value[0]:0.2f} y={y.value[0]:0.2f}')

This gives the solutions:

Solution: x=0.16 y=1.32
Solution: x=6.85 y=-0.30

If it isn't possible to specify constraints then use a multi-start method that tries different initial guess values. This is often required for problems that have many variables and where all solutions are not known. For the 2 variable problem, specify y.value initial guess values to find both solutions.

from gekko import GEKKO
m = GEKKO(remote=False)
x,y = m.Array(m.Var,2)
m.Equations([0.2618 == 0.2094 + x*(1-(1/y))**(1/y),\
             0.9200 == 1-m.exp((-((.5236-.2094)/x)**y))])

y.value = -1
m.solve(disp=False)
print(f'Solution: x={x.value[0]:0.2f} y={y.value[0]:0.2f}')

y.value = 1
m.solve(disp=False)
print(f'Solution: x={x.value[0]:0.2f} y={y.value[0]:0.2f}')

This also gives both solutions:

Solution: x=6.85 y=-0.30
Solution: x=0.16 y=1.32

Sympy can find all solutions for simple problems. See Module 10 (Jupyter Notebook) of the Data Science course.