nonlinear-optimization - 为什么这个 GEKKO 脚本没有产生更好的解决方案？

Question

这是我的代码。我在给定约束 abs(expr1)=abs(expr2) 的情况下最大化表达式 abs(expr1)。

import numpy as np
from gekko import GEKKO

#init
m = GEKKO(remote=False)

x2,x3,x4,x5,x6,x7,x8 = [m.Var(lb=-2*np.pi, ub=2*np.pi) for i in range(7)]

#constraint
m.Equation((1/8)*m.abs(m.cos((1/2)*x6)*(m.sin(x2)*m.sin(x4)*m.sin((1/2)*x7)*((-4)*m.cos(x3)*m.cos((1/2)*x5)*m.cos(x8)+(3+m.cos(x5))*m.sin(x3)*m.sin(x8))+m.cos(x2)*m.cos((1/2)*x7)*m.sin((1/2)*x4)*(8*m.cos(x3)*m.cos(x5)*m.cos(x8)+(-1)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*m.sin(x3)*m.sin(x8)))) == (1/8)*m.abs(m.sin((1/2)*x6)*(4*m.cos(x3)*m.cos(x8)*(m.cos((1/2)*x5)*m.cos((1/2)*x7)*m.sin(x2)*m.sin(x4)+2*m.cos(x2)*m.cos(x5)*m.sin((1/2)*x4)*m.sin((1/2)*x7))+(-1)*m.sin(x3)*((3+m.cos(x5))*m.cos((1/2)*x7)*m.sin(x2)*m.sin(x4)+m.cos(x2)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*m.sin((1/2)*x4)*m.sin((1/2)*x7))*m.sin(x8))))

#objective
m.Obj(-((1/8)*m.abs(m.cos((1/2)*x6)*(m.sin(x2)*m.sin(x4)*m.sin((1/2)*x7)*((-4)*m.cos(x3)*m.cos((1/2)*x5)*m.cos(x8)+(3+m.cos(x5))*m.sin(x3)*m.sin(x8))+m.cos(x2)*m.cos((1/2)*x7)*m.sin((1/2)*x4)*(8*m.cos(x3)*m.cos(x5)*m.cos(x8)+(-1)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*m.sin(x3)*m.sin(x8))))))

#Set global options
m.options.IMODE = 3 

#execute
m.solve()

#output
print('')
print('Results')
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('x4: ' + str(x4.value))
print('x5: ' + str(x5.value))
print('x6: ' + str(x6.value))
print('x7: ' + str(x7.value))
print('x8: ' + str(x8.value))

我得到的解决方案都是 x_i=0，这是一个有效的解决方案，但不是最好的解决方案。例如

x2,x3,x4,x5,x6,x7,x8 = 0.9046, 1.9540, 1.8090, 0, 1.8090, 6.2832, 4.3291

满足约束条件(小数点后第 5 位)并且目标达到 -0.3003(这仍然不是最好的，但它是一个例子)。我尝试使用公差选项但无济于事。请注意，如果我删除等式约束，则目标会正确地达到最大值 -1。

为什么求解器卡在零解上？

Answer 1

Gekko 中的求解器是局部最小化器，而不是全局最小化器。你的问题有很多局部极小值与 sin 和 cos 函数。您可以使用多起点方法或全局方法(例如 simulated annealing)获得全局最小值.我使用多启动方法对您的脚本进行了修改。 -2*np.pi 和2*np.pi 之间的随机值用于初始化变量。

for xi in x:
    xi.value = np.random.rand(20)*4*np.pi - 2*np.pi

IMODE=2 同时解决所有这些情况。

m.options.IMODE = 2

如果您需要执行很多情况，那么您可以 parallelize this calculation with multiple threads .您还应该切换到 m.abs3 而不是 m.abs 以避免零处的非连续导数出现问题。另一种策略是对等式两边进行平方以避免绝对值。这是一个完整的版本:

import numpy as np
from gekko import GEKKO

#init
m = GEKKO(remote=False)

x = [m.Var(lb=-2*np.pi, ub=2*np.pi) for i in range(7)]
for xi in x:
    xi.value = np.random.rand(20)*4*np.pi - 2*np.pi
x2,x3,x4,x5,x6,x7,x8 = x

#constraint
m.Equation((1/8)*m.abs3(m.cos((1/2)*x6)*(m.sin(x2)*m.sin(x4)*m.sin((1/2)*x7)*\
    ((-4)*m.cos(x3)*m.cos((1/2)*x5)*m.cos(x8)+(3+m.cos(x5))*m.sin(x3)*m.sin(x8))+\
    m.cos(x2)*m.cos((1/2)*x7)*m.sin((1/2)*x4)*(8*m.cos(x3)*m.cos(x5)*m.cos(x8)+\
    (-1)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*m.sin(x3)*m.sin(x8)))) == \
    (1/8)*m.abs3(m.sin((1/2)*x6)*(4*m.cos(x3)*m.cos(x8)*(m.cos((1/2)*x5)*\
    m.cos((1/2)*x7)*m.sin(x2)*m.sin(x4)+2*m.cos(x2)*m.cos(x5)*m.sin((1/2)*\
    x4)*m.sin((1/2)*x7))+(-1)*m.sin(x3)*((3+m.cos(x5))*m.cos((1/2)*x7)*\
    m.sin(x2)*m.sin(x4)+m.cos(x2)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*\
    m.sin((1/2)*x4)*m.sin((1/2)*x7))*m.sin(x8))))

#objective
obj = m.Intermediate(-((1/8)*m.abs3(m.cos((1/2)*x6)*(m.sin(x2)*m.sin(x4)*\
    m.sin((1/2)*x7)*((-4)*m.cos(x3)*m.cos((1/2)*x5)*m.cos(x8)+(3+m.cos(x5))*\
    m.sin(x3)*m.sin(x8))+m.cos(x2)*m.cos((1/2)*x7)*m.sin((1/2)*x4)*(8*m.cos(x3)*\
    m.cos(x5)*m.cos(x8)+(-1)*(5*m.cos((1/2)*x5)+3*m.cos((3/2)*x5))*\
    m.sin(x3)*m.sin(x8))))))
m.Obj(obj)

#Set global options
m.options.IMODE = 2

#execute
m.solve()

#output
print('')
print('Best Result')
i = np.argmin(obj.value)
print('x2: ' + str(x2.value[i]))
print('x3: ' + str(x3.value[i]))
print('x4: ' + str(x4.value[i]))
print('x5: ' + str(x5.value[i]))
print('x6: ' + str(x6.value[i]))
print('x7: ' + str(x7.value[i]))
print('x8: ' + str(x8.value[i]))
print('obj: ' + str(obj.value[i]))

有多个目标为 -0.5 的最佳解决方案。

Best Result
x2: -3.1415936876
x3: 6.2545093655
x4: 3.1415896007
x5: -2.0848973806e-05
x6: -4.7122128433
x7: -4.712565114
x8: 0.029076147797
obj: -0.50000008426

Best Result
x2: -3.1416640191
x3: 3.1415941185
x4: 3.1415948958
x5: -3.1416088732
x6: 1.5708701192
x7: -4.7124627728
x8: -3.1415893349
obj: -0.5000000992