Search code examples
pythonnonlinear-optimizationgekko

Gekko not finding solution to non-linear program unless I predefine values of integer variables

I'm currently using Gekko to solve a non-linear program. The program involves sending some input volume through multiple exchanges, and every exchange along the way requires users to pay a one-time fee. Although these exchange networks can get quite complex, the network I have set up to test my solution simply goes 1 --> 2 --> 3, with the arrows representing the exchanges that are used.

To model the use of the exchange, I have added a binary variable Q which is 1 if the exchange is used and 0 if the exchange is not used. Given that the volume has to flow from 1 to 3 in the network above, both exchanges must be used, however, Gekko fails to discover this fact on its own and says that the solution is not found.

As more complex networks will require Gekko to decide for itself whether it should use some exchanges, I would like it to be able to do this on its own.

from gekko import GEKKO

m = GEKKO(remote=False)
m.options.SOLVER = 1

# Constants
vol = m.Const(10, 'Vol')
x = m.Const(10, 'x')
y = m.Const(10, 'y')
k = m.Const(100, 'k')
g = m.Const(0.5, 'gas_fee')



# Inputs and Outputs
i121, i231, o231 = [m.Var(lb=0) for i in range(3)]

# Activation
q121, q231 = [m.Var(lb=0, ub=1, integer=True) for i in range(2)]



# Constraints
m.Equation((y - (k/(x+i121))) == i231)
m.Equation((y - (k/(x+i231))) == o231)

m.Equation(i121 == i121*q121)
m.Equation(i231 == i231*q231)

m.Equation(i121 == vol)


# Objective Function
m.Minimize(g*q121 + g*q231 - o231)
m.solve(disp=True)

I have tried initiating the binary variables with values of 1, in which case the solution is found correctly - so for some reason it seems that Gekko does not try changing the binary variables to 1 on its own to find the solution.

q121, q231 = [m.Var(value = 1, lb=0, ub=1, integer=True) for i in range(2)]

FYI, the correct solution would be:

i121 = 10
q121 = 1
i231 = 5
o231 = 3.33333...
q231 = 1
objective function = -2.3333
Solution

  • Managed to get it working by editing the constraints for activating edges with the binary q variable.

    Instead of using:

    m.Equation(i121 == i121*q121)

    I'm now using:

    m.Equation(i121 <= H*q121)

    Where H is some very large value (larger than the max value for i121)

    This seems to be working even for more complex problems with more edges.