Dealing with Non-Optimal Solutions from Gekko

Question

I'm running into some situations where it seems like Gekko is getting stuck in local maximums and was wondering what approaches could be used to get around this or dig deeper into the cause (including default settings below).

For example, running the scenario below yields an objective of "-5127.34945104756"

m = GEKKO(remote=False)
m.options.NODES = 3
m.options.IMODE = 3
m.options.MAX_ITER = 1000
m.options.SOLVER=1

#Limit max lnuc weeks
m.Equation(sum(x8)<=6)

m.Maximize(m.sum(simu_total_volume))

m.solve(disp = True)

#Objective      :   -5127.34945104756

Now If I simply change "m.Equation(sum(x8)<=6)" to "m.Equation(sum(x8)==6)", it returns a better solution (-5638.55528892101):

m = GEKKO(remote=False)
m.options.NODES = 3
m.options.IMODE = 3
m.options.MAX_ITER = 1000
m.options.SOLVER=1

#Limit max lnuc weeks
m.Equation(sum(x8)==6)

m.Maximize(m.sum(simu_total_volume))

m.solve(disp = True)
# Objective      :   -5638.55528892101

Given that "6" falls in the range of <=6, is there a reason why Gekko wouldn't try to go all the way up to 6 here? Posting the full code/values would also be difficult given size/scale of the problem, so appreciate any feedback based on this.

Answer 1

Gekko solvers are gradient-based Nonlinear Programming (NLP) solvers that find local minima. There are a few strategies to help Gekko find the global optimum.

Here is an example that can help with this important topic of local vs. global minima. The following script produces the local (not global) solution of (7,0,0) with objective 951.0.

from gekko import GEKKO
m = GEKKO(remote=False)
x = m.Array(m.Var,3,lb=0)
x1,x2,x3 = x
m.Minimize(1000-x1**2-2*x2**2-x3**2-x1*x2-x1*x3)
m.Equations([8*x1+14*x2+7*x3==56,
             x1**2+x2**2+x3**2>=25])
m.solve(disp=False)
res=[print(f'x{i+1}: {xi.value[0]}') for i,xi in enumerate(x)]
print(f'Objective: {m.options.objfcnval:.2f}')

There are gradient-based methods for global optimization found in solvers such as BARON, genetic algorithms, simulated annealing, etc. An easy approach is to perform a multi-start method with different initial conditions (guesses) over a grid search or intelligently with a Bayesian approach to more intelligently search if the number of initial guesses is small.

Multi-Start with Parallel Threading

A grid search is easy to parallelize to simultaneously start from multiple locations. Here is the same optimization problem where the global solution is found with a parallelized gekko optimizations.

import numpy as np
import threading
import time, random
from gekko import GEKKO

class ThreadClass(threading.Thread):
    def __init__(self, id, xg):
        s = self
        s.id = id
        s.m = GEKKO(remote=False)
        s.xg = xg
        s.objective = float('NaN')

        # initialize variables
        s.m.x = s.m.Array(s.m.Var,3,lb=0)
        for i in range(3):
            s.m.x[i].value = xg[i]
        s.m.x1,s.m.x2,s.m.x3 = s.m.x

        # Equations
        s.m.Equation(8*s.m.x1+14*s.m.x2+7*s.m.x3==56)
        s.m.Equation(s.m.x1**2+s.m.x2**2+s.m.x3**2>=25)

        # Objective
        s.m.Minimize(1000-s.m.x1**2-2*s.m.x2**2-s.m.x3**2
                     -s.m.x1*s.m.x2-s.m.x1*s.m.x3)

        # Set solver option
        s.m.options.SOLVER = 1

        threading.Thread.__init__(s)

    def run(self):
        print('Running application ' + str(self.id) + '\n')
        self.m.solve(disp=False,debug=0) # solve
        # Retrieve objective if successful
        if (self.m.options.APPSTATUS==1):
            self.objective = self.m.options.objfcnval
        else:
            self.objective = float('NaN')
        self.m.cleanup()

# Optimize at mesh points
x1_ = np.arange(0.0, 10.0, 3.0)
x2_ = np.arange(0.0, 10.0, 3.0)
x3_ = np.arange(0.0, 10.0, 3.0)
x1,x2,x3 = np.meshgrid(x1_,x2_,x3_)

threads = [] # Array of threads

# Load applications
id = 0
for i in range(x1.shape[0]):
    for j in range(x1.shape[1]):
        for k in range(x1.shape[2]):
            xg = (x1[i,j,k],x2[i,j,k],x3[i,j,k])
            # Create new thread
            threads.append(ThreadClass(id, xg))
            # Increment ID
            id += 1

# Run applications simultaneously as multiple threads
# Max number of threads to run at once
max_threads = 8
for t in threads:
    while (threading.activeCount()>max_threads):
        # check for additional threads every 0.01 sec
        time.sleep(0.01)
    # start the thread
    t.start()

# Check for completion
mt = 10.0 # max time (sec)
it = 0.0  # time counter
st = 1.0  # sleep time (sec)
while (threading.active_count()>=3):
    time.sleep(st)
    it = it + st
    print('Active Threads: ' + str(threading.active_count()))
    # Terminate after max time
    if (it>=mt):
        break

# Initialize array for objective
obj = np.empty_like(x1)

# Retrieve objective results
id = 0
id_best = 0; obj_best = 1e10
for i in range(x1.shape[0]):
    for j in range(x1.shape[1]):
        for k in range(x1.shape[2]):
            obj[i,j,k] = threads[id].objective
            if obj[i,j,k]<obj_best:
                id_best = id
                obj_best = obj[i,j,k]
            id += 1

print(obj)
print(f'Best objective {obj_best}')
print(f'Solution {threads[id_best].m.x}')

Bayesian Optimization

Another approach is to intelligently search by mapping the initial conditions to the performance of the optimized solution. It searches in areas where it expects the best performance or where it hasn't been tested and the uncertainty is high.

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

# Define the search space for the hyperparameters
space = {'x1': hp.quniform('x1', 0, 10, 3),
         'x2': hp.quniform('x2', 0, 10, 3),
         'x3': hp.quniform('x3', 0, 10, 3)}

def objective(params):
    m = GEKKO(remote=False)
    x = m.Array(m.Var,3,lb=0)
    x1,x2,x3 = x
    x1.value = params['x1']
    x2.value = params['x2']
    x3.value = params['x3']
    m.Minimize(1000-x1**2-2*x2**2-x3**2-x1*x2-x1*x3)
    m.Equations([8*x1+14*x2+7*x3==56,
                 x1**2+x2**2+x3**2>=25])
    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, 'x':x}

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"Solution: {sol['x']}")

Both multi-start methods find the global solution:

Objective: 936.00
Solution: [[0.0] [0.0] [8.0]]

If you determine that equality constraints always produce the global optimum, then you could also switch from inequality constraints to enforce the constraint boundary.

Additional information on these multi-start methods is in the Engineering Optimization course on the page for Global Optimization and Solver Tuning.