Potential bug in add_sos type 2

[nick-gorman]

Hi Team and thanks for all your hard work on mip, its awesome!

I think there may be a bug in add_sos type 2, where non consecutive weights are allowed to be non-zero. Or maybe I'm using it wrong.

I've provided a worked example below. (not important but the example is based off this blog post which I modified to reproduced the 'bug' I was getting in my own code base. https://medium.com/bcggamma/hands-on-modeling-non-linearity-in-linear-optimization-problems-f9da34c23c9a)

import numpy as np
from mip import Model, xsum, MAXIMIZE, maximize

"""In this example we use a special ordered set of type 2 to constraint y to be the piece linear aproximation of the square of x, and to the objective function. We also set the upper and lower bound of x to equal 346, one of the break points in the linearisation, this means y should perfectly approximate x ** 2 at this point, however it does not appear to do so."""

def square(x):
    return x ** 2

# Generate the data
beta = 1
x_low = 74
x_high = 448
x_samples = np.linspace(x_low, x_high, 12)
y_samples = square(x_samples)

# Instantiate a new model
m = Model("dummy_model", sense=MAXIMIZE)

# Declare the required variables
x = m.add_var(lb=346.0, ub=346.0)
y = m.add_var(lb=0, ub=448 ** 2)
weights = [(m.add_var(lb=0, ub=1), 0) for i in x_samples]

# Define the set of weights and link them to x and y
sos2 = m.add_sos(weights, 2)
weight_sum = m.add_constr(xsum([weights[i][0] for i in range(len(weights))]) == 1)
x_link = m.add_constr(xsum([weights[i][0]*x_samples[i] for i in range(len(x_samples))]) == x)
y_link = m.add_constr(xsum([weights[i][0]*y_samples[i] for i in range(len(y_samples))]) == y)

# Set the objective
m.objective = maximize(xsum([1 * y]))

# Optimize
m.optimize()

# Because x is constrained to be 346 y should equal 346 ** 2
print('y = {}'.format(y.x))
print('x ** 2 = {}'.format(346 ** 2))

[h-g-s]

Hi @nick-gorman , thanks for the detailed bug report. I'll investigate it.

[h-g-s]

Hi @nick-gorman , I added a BINARY var type to variables in weights

weights = [(m.add_var(lb=0, ub=1, var_type=BINARY), 0) for i in x_samples]

and it worked in CBC and Gurobi, could you check ?

[nick-gorman]

Hi @h-g-s, I can confirm this works for the example I gave, however if the weights are binary the piece-wise linear approximation won't work when the value of x isn't exactly equal to one of the linearisation's break points.

[h-g-s]

I never used type 2 SOS, thanks for clarifying. I'm adapting the code to C++ to make tests in Cbc easier.

[h-g-s]

I think I understood why the SOS2 is being ignored. Since there are no integral variables, the C Interface is (wrongly) assuming that there are no objects to branch and solving the LP relaxation only (I'll fix it) . Please try the following, just add an independent binary variable:

k = m.add_var(var_type=BINARY, obj = 1.0)

and see if it works

[nick-gorman]

Yeah that works now

[nick-gorman]

Hi @h-g-s

An interesting new result, the fix of adding the binary variable stopped working when I modified the test code, an example with more break points resulted in the bug coming back.

import numpy as np
from mip import Model, xsum, MAXIMIZE, maximize, BINARY, CONTINUOUS

"""In this example we use a special ordered set of type 2 to constraint y to be the piece linear 
    aproximation of the  square of x, and to the objective function. We also set the upper and lower 
    bound of x to equal 346, one of the break points in the linearisation, this means y should 
    perfectly approximate x ** 2 at this point, however it does not appear to do so."""

def square(x):
    return x ** 2

# Generate the data
x_low = 0.0
x_high = 448
x_samples = np.linspace(x_low, x_high, 50)
y_samples = square(x_samples)

# Instantiate a new model
m = Model("dummy_model", sense=MAXIMIZE)

# Declare the required variables
x = m.add_var(lb=380, ub=380)
y = m.add_var(lb=0, ub=448 ** 2)
weights = [(m.add_var(lb=0, ub=1, var_type=CONTINUOUS), 0) for i in x_samples]

# Define the set of weights and link them to x and y
sos2 = m.add_sos(weights, 2)
weight_sum = m.add_constr(xsum([weights[i][0] for i in range(len(weights))]) == 1)
x_link = m.add_constr(xsum([weights[i][0]*x_samples[i] for i in range(len(x_samples))]) == x)
y_link = m.add_constr(xsum([weights[i][0]*y_samples[i] for i in range(len(y_samples))]) == y)
k = m.add_var(var_type=BINARY, obj=1.0)

# Set the objective
m.objective = maximize(xsum([1 * y]))

# Optimize
m.optimize()

# Because x is constrained to be 346 y should equal 346 ** 2
print('y = {}'.format(y.x))
print('x ** 2 = {}'.format(380 ** 2))

[nick-gorman]

Hi @h-g-s ,

I've done some more thinking about this, and I read this documentation from XPress-MP, section 3.4.5 Special Ordered Sets of type 2 (SOS2), here they talk about the need for another set of constraints see just below Figure 3.5, that force the weight variables to only be non-zero in adjacent pairs.

I wonder should these constraints be created by mip/cbc automatically when an SOS of type 2 is create? Maybe that's not happening and causing the bug?

In the example below I created these constraints explicitly and then used SOS of type 1 to make sure only one binary variable is non-zero. This appears to make the linearisation work correctly, but I wonder if the solver behaves as efficiently as if an SOS of type 2 had been used.

Cheers, Nick

import numpy as np
from mip import Model, xsum, MAXIMIZE, maximize, BINARY, CONTINUOUS

def square(x):
    return x ** 2

# Generate the data
x_low = 0.0
x_high = 500
x_samples = np.linspace(x_low, x_high, 26)
y_samples = square(x_samples)

# Instantiate a new model
m = Model("dummy_model", sense=MAXIMIZE)

# Declare the required variables
x = m.add_var(lb=380, ub=380)
y = m.add_var(lb=0, ub=500 ** 2)
weights = [(m.add_var(lb=0, ub=1, var_type=CONTINUOUS), 0) for i in x_samples]

# Create a set of binary variables and use them to constraint the weights to only be non-zero in adjacent pairs.
binaries = [(m.add_var(lb=0, ub=1, var_type=BINARY), 1) for i in x_samples]
for i in range(0, len(weights)):
    if i == 0:
        m.add_constr(weights[i][0] <= binaries[i][0])
    elif i == len(weights) - 1:
        m.add_constr(weights[i][0] <= binaries[i - 1][0])
    else:
        m.add_constr(weights[i][0] <= binaries[i - 1][0] + binaries[i][0])

sos1 = m.add_sos(binaries, 1)

# Link weights to  x and y
weight_sum = m.add_constr(xsum([weights[i][0] for i in range(len(weights))]) == 1)
x_link = m.add_constr(xsum([weights[i][0]*x_samples[i] for i in range(len(x_samples))]) == x)
y_link = m.add_constr(xsum([weights[i][0]*y_samples[i] for i in range(len(y_samples))]) == y)
k = m.add_var(var_type=BINARY, obj=1.0)

# Set the objective
m.objective = maximize(xsum([1 * y]))

# Optimize
m.optimize()

# Because x is constrained to be 346 y should equal 346 ** 2
print('y = {}'.format(y.x))
print('x ** 2 = {}'.format(380 ** 2))

[h-g-s]

Hi @nick-gorman , I've talked to forrest and read the original papers on SOS. Some of the conclusions until now:

The weights on the SOS should all be different. In fact, they should be monotonically increasing, but it is not necessary to specify them in this format, since they will be sorted after. It seems that when used to linearize a function they must be the x values. 0, 1, 2... should also be suitable.

I'm improving the documentation on SOS and will include more tests/examples in the package.

[nick-gorman]

Hi @h-g-s , that makes a lot of sense, everything seems to work on my end when I use the weights input a you suggest. Thanks for your help and all your work on mip!