How to improve robustness regarding initial guesses?

[willigott]

Let's say I want to maximize

cos(x1) + cos(x2) + cos(x3)

where

x1 + x2 + x3 = n

with 1 <= xi <= n and xi should be an integer.

Then I can do:

from gekko import GEKKO

m = GEKKO()
m.options.SOLVER = 1

# optional solver settings with APOPT
m.solver_options = ['minlp_maximum_iterations 500', \
                    # minlp iterations with integer solution
                    'minlp_max_iter_with_int_sol 10', \
                    # treat minlp as nlp
                    'minlp_as_nlp 0', \
                    # nlp sub-problem max iterations
                    'nlp_maximum_iterations 50', \
                    # 1 = depth first, 2 = breadth first
                    'minlp_branch_method 2', \
                    # maximum deviation from whole number
                    'minlp_integer_tol 0.05', \
                    # covergence tolerance
                    'minlp_gap_tol 0.01']

n = 50

# Initialize variables
x1 = m.Var(value=1, lb=1, ub=n, integer=True)
x2 = m.Var(value=1, lb=1, ub=n, integer=True)
x3 = m.Var(value=1, lb=1, ub=n, integer=True)

# the constraint
m.Equation(x1 + x2 + x3 == n)

# the objective function; we multiply by -1 as we want to maximize
m.Obj(-1*(m.cos(x1) + m.cos(x2) + m.cos(x3)))

m.solve(disp=True)
print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('Objective: ' + str(m.options.objfcnval))

This will print:

--Integer Solution:  -1.96E+00 Lowest Leaf:  -1.96E+00 Gap:   0.00E+00
Iter:     6 I:  0 Tm:      0.00 NLPi:    2 Dpth:    2 Lvs:    1 Obj: -1.96E+00 Gap:  0.00E+00
 Successful solution

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.420000000507571E-002 sec
 Objective      :   -1.95568412506861     
 Successful solution
 ---------------------------------------------------

Results
x1: [1.0]
x2: [24.0]
x3: [25.0]
Objective: -1.95568413

That seems odd as I would expect a value close to 3.

When I then try

x1 = m.Var(value=2, lb=1, ub=n, integer=True)
x2 = m.Var(value=3, lb=1, ub=n, integer=True)
x3 = m.Var(value=4, lb=1, ub=n, integer=True)

I get

 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:    3 Dpth:    0 Lvs:    0 Obj: -4.40E-01 Gap:  0.00E+00
 Successful solution

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.350000000093132E-002 sec
 Objective      :  -0.440460272267080     
 Successful solution
 ---------------------------------------------------

Results
x1: [1.0]
x2: [48.0]
x3: [1.0]
Objective: -0.440460272

And finally, when using

x1 = m.Var(value=3, lb=1, ub=n, integer=True)
x2 = m.Var(value=4, lb=1, ub=n, integer=True)
x3 = m.Var(value=5, lb=1, ub=n, integer=True)

I get

Iter:     7 I:  0 Tm:      0.00 NLPi:    2 Dpth:    2 Lvs:    1 Obj: -2.94E+00 Gap:  0.00E+00
 Successful solution

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.550000000861473E-002 sec
 Objective      :   -2.94007771670051     
 Successful solution
 ---------------------------------------------------

Results
x1: [6.0]
x2: [19.0]
x3: [25.0]
Objective: -2.94007772

So, for three different initial guesses I receive three different objective values (which is, of course, not unusual for numeric optimization).

How do you usually deal with this? Is there a multi-start option or something like this? Thanks!

[APMonitor]

All of the solvers in GEKKO produce local optima for convex nonlinear or mixed integer nonlinear problems. As you referenced, multi-start methods are a popular method for finding a global optimum for non-convex problems. There are also dedicated solvers such as COUENNE and BARON that attempt to find a global solution. Another option in GEKKO is to use a multithreaded application to run the multi-start in parallel. Here is an example of running GEKKO in parallel: http://apmonitor.com/me575/index.php/Main/ParallelComputing

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

class ThreadClass(threading.Thread):
    def __init__(self, id, server, ai, bi):
        s = self
        s.id = id
        s.server = server
        s.m = GEKKO()
        s.a = ai
        s.b = bi
        s.objective = float('NaN')

        # initialize variables
        s.m.x1 = s.m.Var(1,lb=1,ub=5)
        s.m.x2 = s.m.Var(5,lb=1,ub=5)
        s.m.x3 = s.m.Var(5,lb=1,ub=5)
        s.m.x4 = s.m.Var(1,lb=1,ub=5)

        # Equations
        s.m.Equation(s.m.x1*s.m.x2*s.m.x3*s.m.x4>=s.a)
        s.m.Equation(s.m.x1**2+s.m.x2**2+s.m.x3**2+s.m.x4**2==s.b)

        # Objective
        s.m.Obj(s.m.x1*s.m.x4*(s.m.x1+s.m.x2+s.m.x3)+s.m.x3)

        # Set global options
        s.m.options.IMODE = 3 # steady state optimization
        s.m.options.SOLVER = 1 # 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)

        # Results
        #print('')
        #print('Results')
        #print('x1: ' + str(self.m.x1.value))
        #print('x2: ' + str(self.m.x2.value))
        #print('x3: ' + str(self.m.x3.value))
        #print('x4: ' + str(self.m.x4.value))

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

# Select server
server = 'http://byu.apmonitor.com'

# Optimize at mesh points
x = np.arange(20.0, 30.0, 2.0)
y = np.arange(30.0, 50.0, 2.0)
a, b = np.meshgrid(x, y)

# Array of threads
threads = []

# Calculate objective at all meshgrid points

# Load applications
id = 0
for i in range(a.shape[0]):
    for j in range(b.shape[1]):
        # Create new thread
        threads.append(ThreadClass(id, server, a[i,j], b[i,j]))
        # 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_like(a)

# Retrieve objective results
id = 0
for i in range(a.shape[0]):
    for j in range(b.shape[1]):
        obj[i,j] = threads[id].objective
        id += 1

# plot 3D figure of results
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(a, b, obj, \
                       rstride=1, cstride=1, cmap=cm.coolwarm, \
                       vmin = 12, vmax = 22, linewidth=0, antialiased=False)
ax.set_xlabel('a')
ax.set_ylabel('b')
ax.set_zlabel('obj')
ax.set_title('Multi-Threaded GEKKO')
plt.show()