GEKKO optimizer converging on wrong point

I am using the GEKKO optimization package (https://gekko.readthedocs.io/en/latest/index.html) in python. Since my problem sometimes includes integer variables I run with the following options to get optimization for a steady state mixed integer problem (APOPT):

m.options.IMODE = 3
m.options.SOLVER = 1

The optimizer finds a solution that is not optimal. Consider the following minimal working example:

from gekko import GEKKO

m = GEKKO(remote=False)
# Define objective
def obj(x1, x2):
    return x1*0.04824647331436083 + x2*0.1359023029124411 + x1*x2*(0.5659570743890336)

# Define variables
x1 = m.Var(lb=-0.7280376309420935, ub=0.2719623690579065)
x2 = m.Var(lb=-0.8912733608888134, ub=0.10872663911118663)

u = m.Intermediate(obj(x1,x2))
m.Maximize(u)

# Solve the problem
m.options.SOLVER = 'APOPT'
m.solve(disp=False)

# Print the results
print(f'x1: {x1.value[0]}')
print(f'x2: {x2.value[0]}')
print(f'u: {u.value[0]}')

This gives the output:

x1: 0.27196236906
x2: 0.10872663911
u: 0.044632524297

However, consider the following solution and its utility:

x1 = -0.7280376309420935
x2 = -0.8912733608888134
obj(x1, x2)

This results in a utility of 0.2109871851434535

A contour plot verifies that there are two local maxima. Gekko solvers report local minima or maxima when the Karush Kuhn Tucker (KKT) conditions are satisfied.

from gekko import GEKKO

m = GEKKO(remote=False)
# Define objective
def obj(x1, x2):
    return x1*0.04824647331436083 + x2*0.1359023029124411 + x1*x2*(0.5659570743890336)

# Define variables
x1 = m.Var(lb=-0.7280376309420935, ub=0.2719623690579065)
x2 = m.Var(lb=-0.8912733608888134, ub=0.10872663911118663)

u = m.Intermediate(obj(x1,x2))
m.Maximize(u)

# Solve the problem
m.options.SOLVER = 'APOPT'
m.solve(disp=False)

# Print the results
print(f'x1: {x1.value[0]}')
print(f'x2: {x2.value[0]}')
print(f'u: {u.value[0]}')

import numpy as np
import matplotlib.pyplot as plt

# Generate grid of values
x1v = np.linspace(-0.7280376309420935, 0.2719623690579065, 400)
x2v = np.linspace(-0.8912733608888134, 0.10872663911118663, 400)
x1g, x2g = np.meshgrid(x1v, x2v)
ug = obj(x1g, x2g)

# Create contour plot
plt.figure(figsize=(6, 3.5))
contour = plt.contourf(x1g, x2g, ug, levels=50, cmap='viridis')
plt.colorbar(contour); plt.xlabel('x1'); plt.ylabel('x2')
plt.tight_layout(); plt.show()

For problems with multiple local minima or maxima, there are multi-start methods to assist in finding the global optimum: Multi-start implementation with GEKKO Below is an implementation of hyperopt that is used to find the global optimum. hyperopt is a Python package for performing hyperparameter optimization with a variety of optimization algorithms including random search, Tree-structured Parzen Estimator (TPE), and adaptive TPE, as well as a simple and flexible way to define the search space for the hyperparameters.

from gekko import GEKKO
from hyperopt import fmin, tpe, hp
from hyperopt import STATUS_OK, STATUS_FAIL

# Define the search space for the multi-start parameters
space = {'x1': hp.quniform('x1', -0.72, 0.27, 0.1),
         'x2': hp.quniform('x2', -0.89, 0.10, 0.1)}

def objective(params):
    m = GEKKO(remote=False)
    x1 = m.Var(lb=-0.7280376309420935, ub=0.2719623690579065)
    x2 = m.Var(lb=-0.8912733608888134, ub=0.10872663911118663)
    x1.value = params['x1']
    x2.value = params['x2']
    u = m.Intermediate(x1*0.04824647331436083 + x2*0.1359023029124411 + x1*x2*(0.5659570743890336))
    m.Maximize(u)
    m.options.SOLVER = 1
    m.solve(disp=False,debug=False)
    obj = m.options.objfcnval
    if m.options.APPSTATUS==1:
        s=STATUS_OK
    else:
        s=STATUS_FAIL
    m.cleanup()
    return {'loss':obj, 'status': s, 'x1':x1.value[0],'x2':x2.value[0],'u':u.value[0]}

best = fmin(objective, space, algo=tpe.suggest, max_evals=50)
sol = objective(best)
print(f"Solution Status: {sol['status']}")
print(f"Objective: {-sol['loss']:.2f}")
print(f"x1: {sol['x1']}")
print(f"x2: {sol['x2']}")
print(f"u: {sol['u']}")

This produces the global maximum with the TPE algorithm to find the best optimal solution by changing the initial guess values using a Bayesian approach.

Solution Status: ok
Objective: 0.21
x1: -0.72803763094
x2: -0.89127336089
u: 0.21098718514