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: 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.