Lookup values with sos1 in gekko

Question

The optimizer needs to select the appropriate discrete option to maximize an objective. It uses a lookup to correlate the option (0-5) to associated values that can then be used in the optimization definition. Here is a simplified version that demonstrates the issue:

Data

x_data = np.array([0, 1, 2, 3, 4, 5])
y_data = np.array([1000,2000,9000,4500,5000,900])
z_data = np.array([15,13,12,17,11,10])

I found this related post that demonstrates how to use a cspline function to create the relationship between the discrete choice x variable and the dependent y and z variables. I've created the x variable as an sos1 type that selects one from the list of values.

from gekko import GEKKO
import numpy as np

# Define Data Points
x_data = np.array([0, 1, 2, 3, 4, 5])
y_data = np.array([1000,2000,9000,4500,5000,900])
z_data = np.array([15,13,12,17,11,10])

# Create Variables
m = GEKKO()
x = m.sos1(x_data)
y,z = m.Array(m.Var,2)
m.cspline(x, y, x_data, y_data)
m.cspline(x, z, x_data, z_data)

# Define objective and solve
m.Maximize(y*z)
m.options.SOLVER = 1
m.solve(disp=True)

# Print the results
print(f'x: {x.value[0]}, y: {y.value[0]}, z: {z.value[0]}')
print('Objective:', y.value[0]*z.value[0])
print(y_data*z_data)

I get the wrong answer (maybe local solution) when solving this problem.

x: 0.0, y: 1000.0, z: 15.0
Objective: 15000.0

It may have something to do with a local minimum that is created by the cubic spline. Do you have any suggestions to find the maximum for this synthetic function?

Answer 1

Try initializing with a different guess value for x. A local minimum is found when the x<0.3 is used as the initial guess value. With x>=0.3, it finds the global maximum value. If the order of x is arbitrary, remove the local minima by first sorting y_data * z_data to create a new index for x_data. Because x_data are all integer options, you also don't need to use m.sos1() but can use an integer decision variable instead.

x = m.Var(value=1,lb=0,ub=5,integer=True)

The m.sos1() object creates a new binary variable for each option. It is best to use this type when the options are discrete, but not integer values such as select from a list of pipe diameters of 0.5, 1.5, 2.0, 5.0. Here is a complete script:

from gekko import GEKKO
import numpy as np

# Define Data Points
x_data = np.array([0, 1, 2, 3, 4, 5])
y_data = np.array([1000,2000,9000,4500,5000,900])
z_data = np.array([15,13,12,17,11,10])

# Create Variables
m = GEKKO(remote=False)
x = m.Var(value=1,lb=0,ub=5,integer=True)
y,z = m.Array(m.Var,2)
m.cspline(x, y, x_data, y_data)
m.cspline(x, z, x_data, z_data)

# Define objective and solve
m.Maximize(y*z)
m.options.SOLVER = 1
m.solve(disp=False)

# Print the results
print(f'x: {x.value[0]}, y: {y.value[0]}, z: {z.value[0]}')
print('Objective:', y.value[0]*z.value[0])
print(y_data*z_data)

When dealing with local minima, try creating a plot of the objective function with respect to one decision variable at a time. Here is the cspline objective plot of the objective function.

from gekko import GEKKO
import numpy as np

x_data = np.array([0, 1, 2, 3, 4, 5])
y_data = np.array([1000,2000,9000,4500,5000,900])
z_data = np.array([15,13,12,17,11,10])

m = GEKKO(remote=False)
x = m.Param()
y,z = m.Array(m.Var,2)
m.cspline(x, y, x_data, y_data)
m.cspline(x, z, x_data, z_data)

m.Maximize(y*z)
m.options.SOLVER = 1
m.solve(disp=False)

xc = np.linspace(0,5,50)
obj = []
for xi in xc:
    print('x = ', xi)
    x.value = xi
    m.solve(disp=False)
    obj.append(-m.options.OBJFCNVAL)
    
import matplotlib.pyplot as plt
plt.figure(figsize=(6,3))
plt.plot(x_data,y_data*z_data,'ro')
plt.plot(xc,obj,'k--')
plt.xlabel('x'); plt.ylabel('Objective')
plt.grid(); plt.tight_layout()
plt.savefig('cspline.png',dpi=300); plt.show()

Another option is to use a piece-wise linear function m.pwl() instead of cubic splines although the way the problem is posed with slack variables may cause a problem with more complex problems and multiple objective functions. For example, the "curved" interpolation points indicate a potential issue with this function.

from gekko import GEKKO
import numpy as np

# Define Data Points
x_data = np.array([0, 1, 2, 3, 4, 5])
y_data = np.array([1000,2000,9000,4500,5000,900])
z_data = np.array([15,13,12,17,11,10])

# Create Variables
m = GEKKO(remote=False)
x = m.Param()
y,z = m.Array(m.Var,2,lb=0,ub=10000)
m.pwl(x, y, x_data, y_data)
m.pwl(x, z, x_data, z_data)

m.options.SOLVER = 1

xc = np.linspace(0,5,50)
obj = []
for xi in xc:
    print('x = ', xi)
    x.value = xi
    m.solve(disp=False)
    obj.append(y.value[0]*z.value[0])
    
import matplotlib.pyplot as plt
plt.figure(figsize=(6,3))
plt.plot(x_data,y_data*z_data,'ro')
plt.plot(xc,obj,'k--')
plt.xlabel('x'); plt.ylabel('Objective')
plt.grid(); plt.tight_layout()
plt.savefig('pwl.png',dpi=300); plt.show()