用GEKKO实现多个启动

Question

找不到用GEKKO解决nlp优化问题的多启动方法的实现方法。这里有一个以六峰函数为例的例子.六峰函数由于存在多个局部最优而难以优化.多起点方法工作良好，并增加了在全球范围内解决优化问题的机会，如果与基于导数的稳健求解器相结合，将其包括在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]   

Answer 1

谢谢你张贴的多重启动代码。这是一个简单的目标函数的有趣问题。等高线图显示了目标函数的六个局部极小。

​

​

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可以与多个线程并行运行。最优解是橙色点。

​

​

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()

除了找到所有六个局部最优值外，其中一个解在(0,0)处的局部最大值。当局部求解者被指派求解全局最优时，这种情况就会发生。以下是解决办法：

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