Why this GEKKO script does not yield a better solution?

Question

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?

score 1 · Accepted Answer · answered Oct 08 '20 at 10:14

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

Does squaring of both side have an impact on speed? What if the constraint is complex and I need to 'conjugate' square? Is it still faster than abs3 ? — optiman, Oct 08 '20 at 20:27
The squared terms make the model more nonlinear. The `abs3` function adds a binary variable that can be challenging for large-scale systems of equations. I recommend that you try both approaches. — John Hedengren, Oct 08 '20 at 22:17

Why this GEKKO script does not yield a better solution?

1 Answers1