微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

为什么此GEKKO脚本无法提供更好的解决方案?

如何解决为什么此GEKKO脚本无法提供更好的解决方案?

这是我的代码。给定约束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,x8 = 0.9046,1.9540,1.8090,6.2832,4.3291

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

为什么求解器卡在零解中?

解决方法

Gekko中的求解器是局部最小化器,而不是全局最小化器。您遇到的问题是sincos函数有很多局部最小值。您可以使用多重启动方法或全局方法(例如simulated annealing)来获得全局最小值。我使用多启动方法对脚本进行了修改。 -2*np.pi2*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

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。