how to add a set of SOS1 constraint in pyomo

Question

I am solving an optimization problem where I need to assign a carton type to every product, and determine how many cartons of the corresponding type should be used. Every product can be assigned to only 1 carton type. How many units of a specific product a carton of a specific type can hold is known.

I have created decision variables as follows:

import pyomo.environ as pyo
model = pyo.ConcreteModel(name="Pack_size_optim")
model.dv_prod_carton = pyo.Var(prod, cartons, within=pyo.NonNegativeIntegers)

prod, cartons are the lists having unique products and carton types

Now, for my problem, what I want is decision variables : for each product across all carton types should be a SOS 1 set, i.e. for every product only one carton type should have > 0 value.

I was trying something like below, but it is not working :

for i in prod:    
    pyo.SOSConstraint(var = [model.dv_prod_carton[i, j] for j in cartons], sos = 1)

I think var argument in the pyo.SOSConstraint function above should be an 'IndexedVar', but instead I have a list of variables.

Can anybody please help me in what I am trying to achieve above. Thank you.

Answer 1

There are (at least) four ways to model this:

• SOS1. Pyomo has support for SOS1 variables: https://pyomo.readthedocs.io/en/stable/library_reference/kernel/sos.html.
• Binary variables with big-M constraints. This is the most portable (over different solvers) but requires good bounds. I use this when possible.
• Binary variables with indicator constraints.
• Disjunctive programming. This is probably a bit of overkill for this case. Pyomo has support for this: https://pyomo.readthedocs.io/en/stable/contributed_packages/gdpopt.html

I was able to get SOS1 variables to work by explicitly specifying .ref and .sosno suffices. I am sure there are better ways, but at least this seems to work:

# SOS1 test with concrete model

import pyomo.environ as pyo

model = pyo.ConcreteModel()

model.x = pyo.Var([1,2], domain=pyo.PositiveReals, bounds=(0,1))

model.OBJ = pyo.Objective(expr = model.x[1] + model.x[2],sense=pyo.maximize)

# don't exactly know how to get these approaches to work
#pyomo.core.kernel.sos.sos1([model.x[1],model.x[2]])
#pyo.SOSConstraint(var = [model.x[1], model.x[2]], sos = 1)
#model.sos_constraint = pyo.SOSConstraint(var = model.x, sos = 1)
#Message: Solver does not support SOS level 1 constraints 

model.sosno = pyo.Suffix(direction=pyo.Suffix.EXPORT)
model.ref = pyo.Suffix(direction=pyo.Suffix.EXPORT)

# Add entry for each index of model.x
model.sosno.set_value(model.x, 1)
model.ref[model.x[1]] = 1
model.ref[model.x[2]] = 2

# if SOS1 works correctly solution should have one of x[1],x[2] nonzero
# indeed we see with Cplex:
# 'Objective': {'OBJ': {'Value': 1.0}}, 'Variable': {'x[1]': {'Value': 1.0}}

This is a bit as you would do this in AMPL.