Issue with integer optimization with Gekko

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.

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.

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()