Multi-start implementation with GEKKO

user_11968625 0
Viewed 111

Could not find an implementation of a multi-start approach for the solution of nlp optimization problems with GEKKO. Here there is an example using the six-hump function as a case study. The six-hump function is difficult to optimize due to the presence of multiple local optima. The multi-start approach works well and increases the chances to solve optimisation problems globally when combined with robust derivative based solvers as the ones included in GEKKO.

import numpy as np
from gekko import GEKKO
import sobol_seq

# General definition of the problem
lb = np.array([-3.0, -2.0]) # lower bounds
ub = np.array([3.0, 2.0]) # upper bounds
n_dim = lb.shape[0] # number of dimensions
 
# matrix of initial values
multi_start = 10 # number of times the optimisation problem is to be solved
# Sobol  
sobol_space = sobol_seq.i4_sobol_generate(n_dim, multi_start)
x_initial = lb + (ub-lb)*sobol_space # array containing the initial points

# Multi-start optimisation
localsol = [0]*multi_start # storage of local solutions                         
localval = np.zeros((multi_start))

for i in range(multi_start):
    print('multi-start optimization, iteration =', i+1)
    # definition of the problem class with GEKKO
    m = GEKKO(remote=False)
    m.options.SOLVER = 3 
    x = m.Array(m.Var, n_dim)
    # definition of the initial values and bounds for the optimizer
    for j in range(n_dim):
        x[j].value = x_initial[i,j]
        x[j].lower = lb[j]
        x[j].upper = ub[j]
    # Definition of the objective function
    f = (4 - 2.1*x[0]**2 + (x[0]**4)/3)*x[0]**2 + x[0]*x[1] \
          + (-4 + 4*x[1]**2)*x[1]**2 # six-hump function
    # Solving the problem
    m.Obj(f)
    m.solve(disp=False)
    localval[i] = m.options.OBJFCNVAL
    x_local = [x[j].value for j in range(n_dim)]
    localsol[i] = np.array(x_local)

# selecting the best solution
minindex = np.argmin(localval)
x_opt    = localsol[minindex]
f_opt    = localval[minindex]
gekko mathematical-optimization python
1 Answers

Thanks for posting the multi-start code. That is an interesting problem with a simple objective function. A contour plot shows the six local minima of the objective function.

contour plot

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
# Variables at mesh points
x = np.arange(-3, 3, 0.05)
y = np.arange(-2, 2, 0.05)
X,Y = np.meshgrid(x, y)
obj=(4-2.1*X**2+(X**4)/3)*X**2+X*Y \
     +(-4+4*Y**2)*Y**2 # six-hump function
# Create a contour plot
plt.figure()
# Weight contours
CS = plt.contour(X, Y, obj,np.linspace(-1,100,102))
plt.clabel(CS, inline=1, fontsize=8)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

Gekko can run in parallel with multiple threads when searching for a global solution. The optimal solution is the orange dot.

Optimal solution

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

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

        # initialize variables
        s.m.x = s.m.Var(xi,lb=-3,ub=3)
        s.m.y = s.m.Var(yi,lb=-2,ub=2)

        # Objective
        s.m.Minimize((4-2.1*s.m.x**2+(s.m.x**4)/3)*s.m.x**2+s.m.x*s.m.y \
                     + (-4+4*s.m.y**2)*s.m.y**2)

        # Set global options
        s.m.options.SOLVER = 3 # APOPT solver

        threading.Thread.__init__(s)

    def run(self):

        # Don't overload server by executing all scripts at once
        sleep_time = random.random()
        time.sleep(sleep_time)

        print('Running application ' + str(self.id) + '\n')

        # Solve
        self.m.solve(disp=False)

        # Retrieve objective if successful
        if (self.m.options.APPSTATUS==1):
            self.objective = self.m.options.objfcnval
        else:
            self.objective = float('NaN')
        self.m.cleanup()

# General definition of the problem
lb = np.array([-3.0, -2.0]) # lower bounds
ub = np.array([3.0, 2.0]) # upper bounds
n_dim = lb.shape[0] # number of dimensions
 
# matrix of initial values
multi_start = 10 # number of times the optimisation problem is to be solved
# Sobol  
sobol_space = sobol_seq.i4_sobol_generate(n_dim, multi_start)
x_initial = lb + (ub-lb)*sobol_space # array containing the initial points

# Array of threads
threads = []

# Calculate objective for each initial value
id = 0
for i in range(multi_start):
    # Create new thread
    threads.append(ThreadClass(id, x_initial[i,0], x_initial[i,1]))
    # 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 = 3.0 # max time
it = 0.0 # incrementing time
st = 1.0 # sleep time
while (threading.activeCount()>=1):
    time.sleep(st)
    it = it + st
    print('Active Threads: ' + str(threading.activeCount()))
    # Terminate after max time
    if (it>=mt):
        break

# Wait for all threads to complete
#for t in threads:
#    t.join()
#print('Threads complete')

# Initialize array for objective
obj = np.empty(multi_start)
xs = np.empty_like(obj)
ys = np.empty_like(obj)

# Retrieve objective results
id = 0
for i in range(multi_start):
    xs[i] = threads[id].m.x.value[0]
    ys[i] = threads[id].m.y.value[0]
    obj[i] = threads[id].objective
    id += 1

print('Objective',obj)
print('x',xs)
print('y',ys)

best = np.argmin(obj)
print('Lowest Objective',best)
print('Best Objective',obj[best])

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
# Variables at mesh points
x = np.arange(-3, 3, 0.05)
y = np.arange(-2, 2, 0.05)
X,Y = np.meshgrid(x, y)
obj=(4-2.1*X**2+(X**4)/3)*X**2+X*Y \
     +(-4+4*Y**2)*Y**2 # six-hump function
# Create a contour plot
plt.figure()
# Weight contours
CS = plt.contour(X, Y, obj,np.linspace(-1,100,102))
plt.plot(xs,ys,'rx')
plt.plot(xs[best],ys[best],color='orange',marker='o',alpha=0.7)
plt.clabel(CS, inline=1, fontsize=8)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

In additional to finding all six local optimal values, one of the solutions is at a local maximum at (0,0). This can happen when local solvers are tasked with solving a global optimum. Here is the solution:

Objective [-1.80886276e-35 -2.15463824e-01 -2.15463824e-01 -2.15463824e-01
  2.10425031e+00 -1.03162845e+00 -2.15463824e-01  2.10425031e+00
 -1.03162845e+00 -2.15463824e-01]
x [-1.32794585e-19  1.70360672e+00 -1.70360672e+00  1.70360672e+00
  1.60710475e+00  8.98420119e-02 -1.70360672e+00 -1.60710477e+00
 -8.98420136e-02  1.70360671e+00]
y [ 2.11414394e-18 -7.96083565e-01  7.96083565e-01 -7.96083569e-01
  5.68651455e-01 -7.12656403e-01  7.96083569e-01 -5.68651452e-01
  7.12656407e-01 -7.96083568e-01]
Lowest Objective 5
Best Objective -1.0316284535
user_23669410
Related