CODEWITHSUNDEEP
pythonnumpyoptimizationintegergekko
youlyouz ali
youlyouz ali

Reputation: 91

Issue with integer optimization with Gekko

I tried with this code and the result was (2,2). Obviously it's not true, since (0,4) is a best solution. could you please tell me where is it the problem?

import numpy as np
from gekko import GEKKO
m = GEKKO()

x=[m.Var(lb=0,ub=4,integer=True) for i in range(2)]

m.options.SOLVER = 1

#Equations

m.Equation(m.sum(x)>=4)

#Objective

m.Obj(x[0]*x[1])
m.solve(disp=False)
print(x[0].value)
print(x[1].value)

Upvotes: 4

Views: 317

Answers (1)

John Hedengren
John Hedengren

Reputation: 14346

All of the solvers (m.options.SOLVER with 1=APOPT, 2=BPOPT, 3=IPOPT) find the local maximum point (2,2) with initial guess (0,0) and other initial guesses as well. There are two local optima at (0,4) and (4,0). This is a non-convex optimization problem. The solvers are nonlinear convex optimization solvers so they "should" find one of the local solutions. This problem reminds me of a saddle point problem if you rotate the graph to orient with the x0=x1 line. Optimizers can have problems with saddle points if there is no check for an indefinite Hessian.

3d plot

Better initial guesses help the solver find the optimal solution such as x0=1 and x1=3. Also, an upper bound <1.999 on one of the variables or a lower bound of >2.001 also finds the optimal value. The constraints limit the search region so that it is a convex optimization problem and avoids the saddle point.

from gekko import GEKKO
m = GEKKO()

x = m.Array(m.Var,2,value=1,lb=0,ub=4,integer=True)
x[0].value = 1; x[1].value = 3 # initial guesses
#x[0].UPPER=1
m.Equation(m.sum(x)>=4)
m.Minimize(x[0]*x[1])

m.options.SOLVER = 1
m.solve(disp=True)
print(x[0].value)
print(x[1].value)

I created a contour and 3D plot of the optimization problem. Because the problem is non-convex, you'd probably want to switch to a multi-start method or use a non-convex optimization solver.

contour plot

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Design variables at mesh points
x0 = np.arange(0.0, 4.01, 0.1)
x1 = np.arange(0.0, 4.01, 0.1)
x0,x1 = np.meshgrid(x0,x1)

# Equation sum(x)>4
xsum = x0+x1

# Objective Minimize(x0*x1)
xobj = x0*x1

# Create a contour plot
plt.figure()
CS = plt.contour(x0, x1, xobj)
plt.clabel(CS, inline=1, fontsize=10)
CS = plt.contour(x0, x1, xsum,[3.6,3.8,4.0],\
                 colors='k',linewidths=[0.5,1,4])
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel('x0')
plt.ylabel('x1')
plt.savefig('contour.png')

# Create 3D plot
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(x0, x1, xobj, \
            rstride=1, cstride=1, cmap=cm.coolwarm, \
            vmin = 0, vmax = 10, linewidth=0, \
            antialiased=False)
plt.show()

Upvotes: 2

Related Questions