Special Ordered Sets (SOS) Constraints#
Special Ordered Sets (SOS) are a constraint type used in mixed-integer programming to model situations where only one or two variables from an ordered set can be non-zero. Linopy supports both SOS Type 1 and SOS Type 2 constraints.
Overview#
SOS constraints are particularly useful for:
SOS1: Modeling mutually exclusive choices (e.g., selecting one facility from multiple locations)
SOS2: Piecewise linear approximations of nonlinear functions
Improving branch-and-bound efficiency in mixed-integer programming
Types of SOS Constraints#
SOS Type 1 (SOS1)#
In an SOS1 constraint, at most one variable in the ordered set can be non-zero.
Example use cases: - Facility location problems (choose one location among many) - Technology selection (choose one technology option) - Mutually exclusive investment decisions
SOS Type 2 (SOS2)#
In an SOS2 constraint, at most two adjacent variables in the ordered set can be non-zero. Adjacency follows the order the members are declared in (see Ordering).
Example use cases: - Piecewise linear approximation of nonlinear functions - Portfolio optimization with discrete risk levels - Production planning with discrete capacity levels
Basic Usage#
Adding SOS Constraints#
To add SOS constraints to variables in linopy:
import linopy
import pandas as pd
import xarray as xr
# Create model
m = linopy.Model()
# Create variables with numeric coordinates
coords = pd.Index([0, 1, 2], name="options")
x = m.add_variables(coords=[coords], name="x", lower=0, upper=1)
# Add SOS1 constraint
m.add_sos_constraints(x, sos_type=1, sos_dim="options")
# For SOS2 constraint
breakpoints = pd.Index([0.0, 1.0, 2.0], name="breakpoints")
lambdas = m.add_variables(coords=[breakpoints], name="lambdas", lower=0, upper=1)
m.add_sos_constraints(lambdas, sos_type=2, sos_dim="breakpoints")
Method Signature#
Model.add_sos_constraints(variable, sos_type, sos_dim, big_m=None)
Parameters:
variableVariableThe variable to which the SOS constraint should be applied
sos_type{1, 2}Type of SOS constraint (1 or 2)
sos_dimstrName of the dimension along which the SOS constraint applies
big_mfloat | NoneCustom Big-M value for reformulation (see SOS Reformulation)
Requirements:
The specified dimension must exist in the variable
Only one SOS constraint can be applied per variable
Ordering#
A set runs in the order its members are declared along sos_dim. For
sos_type=2 that decides which members are adjacent; an SOS1 set is
unordered, so it does not apply.
The coordinates need not be numeric. Where they are numeric but neither ascend nor descend, linopy warns: such a set was previously ordered by coordinate value, so its meaning changes. Sort the index to keep that ordering. Descending coordinates are not warned about, because reversing a set does not change which of its members are adjacent.
The weights written to an LP file, or passed to a solver that supports SOS sets natively, are the coordinates themselves where they ascend, and member positions otherwise.
Examples#
Example 1: Facility Location (SOS1)#
import linopy
import pandas as pd
import xarray as xr
# Problem data
locations = pd.Index([0, 1, 2, 3], name="locations")
costs = xr.DataArray([100, 150, 120, 80], coords=[locations])
benefits = xr.DataArray([200, 300, 250, 180], coords=[locations])
# Create model
m = linopy.Model()
# Decision variables: build facility at location i
build = m.add_variables(coords=[locations], name="build", lower=0, upper=1)
# SOS1 constraint: at most one facility can be built
m.add_sos_constraints(build, sos_type=1, sos_dim="locations")
# Objective: maximize net benefit
net_benefit = benefits - costs
m.add_objective(-((net_benefit * build).sum()))
# Solve
m.solve(solver_name="gurobi")
if m.status == "ok":
solution = build.solution.to_pandas()
selected_location = solution[solution > 0.5].index[0]
print(f"Build facility at location {selected_location}")
Example 2: Piecewise Linear Approximation (SOS2)#
import numpy as np
# Approximate f(x) = x² over [0, 3] with breakpoints
breakpoints = pd.Index([0, 1, 2, 3], name="breakpoints")
x_vals = xr.DataArray(breakpoints.to_series())
y_vals = x_vals**2
# Create model
m = linopy.Model()
# SOS2 variables (interpolation weights)
lambdas = m.add_variables(lower=0, upper=1, coords=[breakpoints], name="lambdas")
m.add_sos_constraints(lambdas, sos_type=2, sos_dim="breakpoints")
# Interpolated coordinates
x = m.add_variables(name="x", lower=0, upper=3)
y = m.add_variables(name="y", lower=0, upper=9)
# Constraints
m.add_constraints(lambdas.sum() == 1, name="convexity")
m.add_constraints(x == lambdas @ x_vals, name="x_interpolation")
m.add_constraints(y == lambdas @ y_vals, name="y_interpolation")
m.add_constraints(x >= 1.5, name="x_minimum")
# Objective: minimize approximated function value
m.add_objective(y)
# Solve
m.solve(solver_name="gurobi")
Working with Multi-dimensional Variables#
SOS constraints are created for each dimension that is not sos_dim.
# Multi-period production planning
periods = pd.Index(range(3), name="periods")
modes = pd.Index([0, 1, 2], name="modes")
# 2D variables: periods × modes
period_modes = m.add_variables(
lower=0, upper=1, coords=[periods, modes], name="use_mode"
)
# Adds SOS1 constraint for each period
m.add_sos_constraints(period_modes, sos_type=1, sos_dim="modes")
Accessing SOS Variables#
You can easily identify and access variables with SOS constraints:
# Get all variables with SOS constraints
sos_variables = m.variables.sos
print(f"SOS variables: {list(sos_variables.keys())}")
# Check SOS properties of a variable
for var_name in sos_variables:
var = m.variables[var_name]
sos_type = var.attrs["sos_type"]
sos_dim = var.attrs["sos_dim"]
print(f"{var_name}: SOS{sos_type} on dimension '{sos_dim}'")
Variable Representation#
Variables with SOS constraints show their SOS information in string representations:
print(build)
# Output: Variable (locations: 4) - sos1 on locations
# -----------------------------------------------
# [0]: build[0] ∈ [0, 1]
# [1]: build[1] ∈ [0, 1]
# [2]: build[2] ∈ [0, 1]
# [3]: build[3] ∈ [0, 1]
LP File Export#
The generated LP file will include a SOS section:
sos
s0: S1 :: x0:0 x1:1 x2:2
s3: S2 :: x3:0.0 x4:1.0 x5:2.0
Solver Compatibility#
SOS constraints are supported by most modern mixed-integer programming solvers through the LP file format:
Supported solvers (via LP file):
Gurobi
CPLEX
COIN-OR CBC
SCIP
Xpress
Direct API support:
Gurobi (via
gurobipy)
Unsupported solvers:
HiGHS (does not support SOS constraints)
GLPK
MOSEK
MindOpt
For unsupported solvers, use automatic reformulation (see below).
SOS Reformulation#
For solvers without native SOS support, linopy can reformulate SOS constraints as binary + linear constraints using the Big-M method.
# Automatic reformulation during solve (apply / undo bracketed by Model.solve)
m.solve(solver_name="highs", reformulate_sos=True)
# Or stage the reformulation manually — e.g. to inspect or export the MILP
m.apply_sos_reformulation()
m.solve(solver_name="highs")
m.undo_sos_reformulation()
Requirements:
Variables must have non-negative lower bounds (lower >= 0)
Big-M values are derived from variable upper bounds
For infinite upper bounds, specify custom values via the
big_mparameter
# Finite bounds (default)
x = m.add_variables(lower=0, upper=100, coords=[idx], name="x")
m.add_sos_constraints(x, sos_type=1, sos_dim="i")
# Infinite upper bounds: specify Big-M
x = m.add_variables(lower=0, upper=np.inf, coords=[idx], name="x")
m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=10)
The reformulation uses the tighter of big_m and variable upper bound.
Mathematical Formulation#
SOS1 Reformulation:
Original constraint: \(\text{SOS1}(\{x_1, x_2, \ldots, x_n\})\) means at most one \(x_i\) can be non-zero.
Given \(x = (x_1, \ldots, x_n) \in \mathbb{R}^n_+\), introduce binary \(y = (y_1, \ldots, y_n) \in \{0,1\}^n\):
where \(M_i \geq \sup\{x_i\}\) (upper bound on \(x_i\)).
SOS2 Reformulation:
Original constraint: \(\text{SOS2}(\{x_1, x_2, \ldots, x_n\})\) means at most two \(x_i\) can be non-zero, and if two are non-zero, they must have consecutive indices.
Given \(x = (x_1, \ldots, x_n) \in \mathbb{R}^n_+\), introduce binary \(y = (y_1, \ldots, y_{n-1}) \in \{0,1\}^{n-1}\):
where \(M_i \geq \sup\{x_i\}\). Interpretation: \(y_i = 1\) activates interval \([i, i+1]\), allowing \(x_i\) and \(x_{i+1}\) to be non-zero.
Common Patterns#
Piecewise Linear Cost Function#
Note
For a higher-level API that handles all the SOS2 bookkeeping automatically, see Piecewise Linear Constraints.
def add_piecewise_cost(model, variable, breakpoints, costs):
"""Add piecewise linear cost function using SOS2."""
n_segments = len(breakpoints)
lambda_coords = pd.Index(range(n_segments), name="segments")
lambdas = model.add_variables(
coords=[lambda_coords], name="cost_lambdas", lower=0, upper=1
)
model.add_sos_constraints(lambdas, sos_type=2, sos_dim="segments")
cost_var = model.add_variables(name="cost", lower=0)
x_vals = xr.DataArray(breakpoints, coords=[lambda_coords])
c_vals = xr.DataArray(costs, coords=[lambda_coords])
model.add_constraints(lambdas.sum() == 1, name="cost_convexity")
model.add_constraints(variable == (x_vals * lambdas).sum(), name="cost_x_def")
model.add_constraints(cost_var == (c_vals * lambdas).sum(), name="cost_def")
return cost_var
Mutually Exclusive Investments#
def add_exclusive_investments(model, projects, costs, returns):
"""Add mutually exclusive investment decisions using SOS1."""
project_coords = pd.Index(projects, name="projects")
invest = model.add_variables(
coords=[project_coords], name="invest", binary=True
)
model.add_sos_constraints(invest, sos_type=1, sos_dim="projects")
total_cost = (invest * costs).sum()
total_return = (invest * returns).sum()
return invest, total_cost, total_return
See Also#
Creating Variables: Creating variables with coordinates
Creating Constraints: Adding regular constraints
Overview: General linopy usage patterns