Reputation: 33
Here is my code. I'm maximizing an expression abs(expr1) given the constraint abs(expr1)=abs(expr2).
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))
The solution I get is all x_i=0 which is a valid solution but not the best one. For example
x2,x3,x4,x5,x6,x7,x8 = 0.9046, 1.9540, 1.8090, 0, 1.8090, 6.2832, 4.3291
satisfies the constraint (up to the 5th decimal place) and the objectives reaches -0.3003 (which is still not the best but it is an example). I tried to play with the tolerance options but to no avail. Note that if I remove the equality constraint the objective properly reaches the maximal value -1.
Why does the solver get stuck with the zero solution?
Upvotes: 3
Views: 119
Reputation: 14401
The solvers in Gekko are local minimizers, not global minimizers. You problem has many local minima with the sin
and cos
functions. You can get a global minimum by using a multi-start method or a global approach such as simulated annealing. I used a modification of your script with a multi-start method. Random values between -2*np.pi
and 2*np.pi
are used to initialize the variables.
for xi in x:
xi.value = np.random.rand(20)*4*np.pi - 2*np.pi
IMODE=2
solves all of these cases simultaneously.
m.options.IMODE = 2
If you need to perform many cases then you can parallelize this calculation with multiple threads. You should also switch to m.abs3
instead of m.abs
to avoid problems with a non-continuous derivative at zero. Another strategy is to square both sides of your equation to avoid the absolute value. Here is a complete version:
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]))
There are multiple best solutions that have an objective of -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
Upvotes: 1