处理Gekko的非最佳解
我遇到了一些情况，感觉Gekko在局部最优解上卡住了。我想知道有什么方法可以绕过这个问题，或者更深入地了解原因（包括下面的默认设置）。
比如，运行下面的场景会得到一个目标值为“-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
现在如果我简单地把“m.Equation(sum(x8)<=6)”改成“m.Equation(sum(x8)==6)”，它会返回一个更好的解（-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
考虑到“6”在“<=6”的范围内，为什么Gekko不尝试把值提升到6呢？由于问题的规模和复杂性，发布完整的代码和数值也很困难，所以希望能根据这些信息得到一些反馈。
1 个回答
1
Gekko求解器是一种基于梯度的非线性规划（NLP）求解器，主要用于寻找局部最小值。为了帮助Gekko找到全局最优解，有几种策略可以使用。
下面是一个例子，可以帮助理解局部最小值和全局最小值之间的重要区别。以下脚本会产生一个局部解（而不是全局解），结果是（7,0,0），目标值为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}')
在一些求解器中，比如BARON、遗传算法、模拟退火等，都有基于梯度的全局优化方法。一种简单的方法是使用多次启动，尝试不同的初始条件（猜测），可以通过网格搜索的方式来进行，或者用贝叶斯方法更智能地搜索，特别是在初始猜测数量较少时。
多次启动与并行处理
网格搜索很容易进行并行处理，可以同时从多个位置开始。这里是同一个优化问题，通过并行化的Gekko优化找到了全局解。
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}')
贝叶斯优化
另一种方法是通过智能搜索，将初始条件与优化解的表现进行映射。它会在预期表现最佳的区域进行搜索，或者在尚未测试且不确定性较高的地方进行探索。
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']}")
这两种多次启动的方法都能找到全局解：
Objective: 936.00
Solution: [[0.0] [0.0] [8.0]]
如果你发现等式约束总是能产生全局最优解，那么你也可以考虑将不等式约束改为等式约束，以确保约束的边界。
关于这些多次启动方法的更多信息，可以在工程优化课程的全局优化和求解器调优页面找到。