Skip to content

feat: add sos reformulations into linopy to simplify adoption of new sos features #549

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
FBumann wants to merge 14 commits into
base: master
Choose a base branch
from
Open

feat: add sos reformulations into linopy to simplify adoption of new sos features #549

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
d62762f
The SOS constraint reformulation feature has been implemented success…
FBumann Jan 20, 2026
80d5752
Documentation Summary
FBumann Jan 20, 2026
3f6db06
Added Tests for Multi-dimensional SOS
FBumann Jan 20, 2026
a906e71
Add custom big_m parameter for SOS reformulation
FBumann Jan 20, 2026
2258692
Add custom big_m parameter for SOS reformulation
FBumann Jan 20, 2026
0a3dde5
Simplification summary:
FBumann Jan 20, 2026
34b00b0
Revert some docs changes to be more surgical
FBumann Jan 20, 2026
da1dc61
Add math to docs
FBumann Jan 20, 2026
9b0744c
Improve docs
FBumann Jan 20, 2026
7bac84d
Code simplifications:
FBumann Jan 20, 2026
33f7148
Fix mypy
FBumann Jan 20, 2026
cfd7c06
Fix mypy
FBumann Jan 20, 2026
a1f5c5d
Add constants for sos attr keys
FBumann Jan 20, 2026
69610eb
Add release notes
FBumann Jan 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions doc/release_notes.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ Release Notes
.. Upcoming Version

* Fix LP file writing for negative zero (-0.0) values that produced invalid syntax like "+-0.0" rejected by Gurobi
* Add SOS1 and SOS2 reformulations for solvers not supporting them.

Version 0.6.0
--------------
Expand Down
81 changes: 80 additions & 1 deletion doc/sos-constraints.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ Method Signature

.. code-block:: python

Model.add_sos_constraints(variable, sos_type, sos_dim)
Model.add_sos_constraints(variable, sos_type, sos_dim, big_m=None)

**Parameters:**

Expand All @@ -85,6 +85,8 @@ Method Signature
Type of SOS constraint (1 or 2)
- ``sos_dim`` : str
Name of the dimension along which the SOS constraint applies
- ``big_m`` : float | None
Custom Big-M value for reformulation (see :ref:`sos-reformulation`)

**Requirements:**

Expand Down Expand Up @@ -254,6 +256,83 @@ SOS constraints are supported by most modern mixed-integer programming solvers t
- MOSEK
- MindOpt

For unsupported solvers, use automatic reformulation (see below).

.. _sos-reformulation:

SOS Reformulation
-----------------

For solvers without native SOS support, linopy can reformulate SOS constraints
as binary + linear constraints using the Big-M method.

.. code-block:: python

# Automatic reformulation during solve
m.solve(solver_name="highs", reformulate_sos=True)

# Or reformulate manually
m.reformulate_sos_constraints()
m.solve(solver_name="highs")

**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_m`` parameter

.. code-block:: python

# 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: :math:`\text{SOS1}(\{x_1, x_2, \ldots, x_n\})` means at most one
:math:`x_i` can be non-zero.

Given :math:`x = (x_1, \ldots, x_n) \in \mathbb{R}^n_+`, introduce binary
:math:`y = (y_1, \ldots, y_n) \in \{0,1\}^n`:

.. math::

x_i &\leq M_i \cdot y_i \quad \forall i \in \{1, \ldots, n\} \\
x_i &\geq 0 \quad \forall i \in \{1, \ldots, n\} \\
\sum_{i=1}^{n} y_i &\leq 1 \\
y_i &\in \{0, 1\} \quad \forall i \in \{1, \ldots, n\}

where :math:`M_i \geq \sup\{x_i\}` (upper bound on :math:`x_i`).

**SOS2 Reformulation:**

Original constraint: :math:`\text{SOS2}(\{x_1, x_2, \ldots, x_n\})` means at most two
:math:`x_i` can be non-zero, and if two are non-zero, they must have consecutive indices.

Given :math:`x = (x_1, \ldots, x_n) \in \mathbb{R}^n_+`, introduce binary
:math:`y = (y_1, \ldots, y_{n-1}) \in \{0,1\}^{n-1}`:

.. math::

x_1 &\leq M_1 \cdot y_1 \\
x_i &\leq M_i \cdot (y_{i-1} + y_i) \quad \forall i \in \{2, \ldots, n-1\} \\
x_n &\leq M_n \cdot y_{n-1} \\
x_i &\geq 0 \quad \forall i \in \{1, \ldots, n\} \\
\sum_{i=1}^{n-1} y_i &\leq 1 \\
y_i &\in \{0, 1\} \quad \forall i \in \{1, \ldots, n-1\}

where :math:`M_i \geq \sup\{x_i\}`. Interpretation: :math:`y_i = 1` activates interval
:math:`[i, i+1]`, allowing :math:`x_i` and :math:`x_{i+1}` to be non-zero.

Common Patterns
---------------

Expand Down
5 changes: 5 additions & 0 deletions linopy/constants.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,11 @@
CV_DIM,
]

# SOS constraint attribute keys
SOS_TYPE_ATTR = "sos_type"
SOS_DIM_ATTR = "sos_dim"
SOS_BIG_M_ATTR = "big_m_upper"


class ModelStatus(Enum):
"""
Expand Down
10 changes: 5 additions & 5 deletions linopy/io.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@

from linopy import solvers
from linopy.common import to_polars
from linopy.constants import CONCAT_DIM
from linopy.constants import CONCAT_DIM, SOS_DIM_ATTR, SOS_TYPE_ATTR
from linopy.objective import Objective

if TYPE_CHECKING:
Expand Down Expand Up @@ -374,8 +374,8 @@ def sos_to_file(

for name in names:
var = m.variables[name]
sos_type = var.attrs["sos_type"]
sos_dim = var.attrs["sos_dim"]
sos_type = var.attrs[SOS_TYPE_ATTR]
sos_dim = var.attrs[SOS_DIM_ATTR]

other_dims = [dim for dim in var.labels.dims if dim != sos_dim]
for var_slice in var.iterate_slices(slice_size, other_dims):
Expand Down Expand Up @@ -773,8 +773,8 @@ def to_gurobipy(
if m.variables.sos:
for var_name in m.variables.sos:
var = m.variables.sos[var_name]
sos_type: int = var.attrs["sos_type"] # type: ignore[assignment]
sos_dim: str = var.attrs["sos_dim"] # type: ignore[assignment]
sos_type: int = var.attrs[SOS_TYPE_ATTR] # type: ignore[assignment]
sos_dim: str = var.attrs[SOS_DIM_ATTR] # type: ignore[assignment]

def add_sos(s: xr.DataArray, sos_type: int, sos_dim: str) -> None:
s = s.squeeze()
Expand Down
105 changes: 93 additions & 12 deletions linopy/model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,9 @@
GREATER_EQUAL,
HELPER_DIMS,
LESS_EQUAL,
SOS_BIG_M_ATTR,
SOS_DIM_ATTR,
SOS_TYPE_ATTR,
TERM_DIM,
ModelStatus,
TerminationCondition,
Expand Down Expand Up @@ -65,6 +68,7 @@
IO_APIS,
available_solvers,
)
from linopy.sos_reformulation import reformulate_all_sos
from linopy.types import (
ConstantLike,
ConstraintLike,
Expand Down Expand Up @@ -556,6 +560,7 @@ def add_sos_constraints(
variable: Variable,
sos_type: Literal[1, 2],
sos_dim: str,
big_m: float | None = None,
) -> None:
"""
Add an sos1 or sos2 constraint for one dimension of a variable
Expand All @@ -569,15 +574,26 @@ def add_sos_constraints(
Type of SOS
sos_dim : str
Which dimension of variable to add SOS constraint to
big_m : float | None, optional
Big-M value for SOS reformulation. Only used when reformulating
SOS constraints for solvers that don't support them natively.

- None (default): Use variable upper bounds as Big-M
- float: Custom Big-M value

The reformulation uses the tighter of big_m and variable upper bound:
M = min(big_m, var.upper).

Tighter Big-M values improve LP relaxation quality and solve time.
"""
if sos_type not in (1, 2):
raise ValueError(f"sos_type must be 1 or 2, got {sos_type}")
if sos_dim not in variable.dims:
raise ValueError(f"sos_dim must name a variable dimension, got {sos_dim}")

if "sos_type" in variable.attrs or "sos_dim" in variable.attrs:
existing_sos_type = variable.attrs.get("sos_type")
existing_sos_dim = variable.attrs.get("sos_dim")
if SOS_TYPE_ATTR in variable.attrs or SOS_DIM_ATTR in variable.attrs:
existing_sos_type = variable.attrs.get(SOS_TYPE_ATTR)
existing_sos_dim = variable.attrs.get(SOS_DIM_ATTR)
raise ValueError(
f"variable already has an sos{existing_sos_type} constraint on {existing_sos_dim}"
)
Expand All @@ -589,7 +605,12 @@ def add_sos_constraints(
f"but got {variable.coords[sos_dim].dtype}"
)

variable.attrs.update(sos_type=sos_type, sos_dim=sos_dim)
# Process and store big_m value
attrs_update: dict[str, Any] = {SOS_TYPE_ATTR: sos_type, SOS_DIM_ATTR: sos_dim}
if big_m is not None:
attrs_update[SOS_BIG_M_ATTR] = float(big_m)

variable.attrs.update(attrs_update)

def add_constraints(
self,
Expand Down Expand Up @@ -830,18 +851,65 @@ def remove_sos_constraints(self, variable: Variable) -> None:
-------
None.
"""
if "sos_type" not in variable.attrs or "sos_dim" not in variable.attrs:
if SOS_TYPE_ATTR not in variable.attrs or SOS_DIM_ATTR not in variable.attrs:
raise ValueError(f"Variable '{variable.name}' has no SOS constraints")

sos_type = variable.attrs["sos_type"]
sos_dim = variable.attrs["sos_dim"]
sos_type = variable.attrs[SOS_TYPE_ATTR]
sos_dim = variable.attrs[SOS_DIM_ATTR]

del variable.attrs["sos_type"], variable.attrs["sos_dim"]
del variable.attrs[SOS_TYPE_ATTR], variable.attrs[SOS_DIM_ATTR]

# Also remove big_m attribute if present
variable.attrs.pop(SOS_BIG_M_ATTR, None)

logger.debug(
f"Removed sos{sos_type} constraint on {sos_dim} from {variable.name}"
)

def reformulate_sos_constraints(self, prefix: str = "_sos_reform_") -> list[str]:
"""
Reformulate SOS constraints as binary + linear constraints.

This converts SOS1 and SOS2 constraints into equivalent binary variable
formulations using the Big-M method. This allows solving models with SOS
constraints using solvers that don't support them natively (e.g., HiGHS, GLPK).

Big-M values are determined as follows:
1. If custom big_m was specified in add_sos_constraints(), use that
2. Otherwise, use the variable bounds (tightest valid Big-M)

Parameters
----------
prefix : str, optional
Prefix for naming auxiliary variables and constraints.
Default is "_sos_reform_".

Returns
-------
list[str]
List of variable names that were reformulated.

Raises
------
ValueError
If any SOS variable has infinite bounds and no custom big_m was specified.

Examples
--------
>>> m = Model()
>>> x = m.add_variables(
... lower=0, upper=1, coords=[pd.Index([0, 1, 2], name="i")], name="x"
... )
>>> m.add_sos_constraints(x, sos_type=1, sos_dim="i")
>>> m.reformulate_sos_constraints() # Now solvable with HiGHS
['x']

With custom big_m for tighter relaxation:

>>> m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=0.5)
"""
return reformulate_all_sos(self, prefix=prefix)

def remove_objective(self) -> None:
"""
Remove the objective's linear expression from the model.
Expand Down Expand Up @@ -1126,6 +1194,7 @@ def solve(
remote: RemoteHandler | OetcHandler = None, # type: ignore
progress: bool | None = None,
mock_solve: bool = False,
reformulate_sos: bool = False,
**solver_options: Any,
) -> tuple[str, str]:
"""
Expand Down Expand Up @@ -1195,6 +1264,11 @@ def solve(
than 10000 variables and constraints.
mock_solve : bool, optional
Whether to run a mock solve. This will skip the actual solving. Variables will be set to have dummy values
reformulate_sos : bool, optional
Whether to automatically reformulate SOS constraints as binary + linear
constraints for solvers that don't support them natively.
This uses the Big-M method and requires all SOS variables to have finite bounds.
Default is False.
**solver_options : kwargs
Options passed to the solver.

Expand Down Expand Up @@ -1293,10 +1367,17 @@ def solve(
)

# SOS constraints are not supported by all solvers
if self.variables.sos and not solver_supports(
solver_name, SolverFeature.SOS_CONSTRAINTS
):
raise ValueError(f"Solver {solver_name} does not support SOS constraints.")
if self.variables.sos:
if reformulate_sos and not solver_supports(
solver_name, SolverFeature.SOS_CONSTRAINTS
):
logger.info(f"Reformulating SOS constraints for solver {solver_name}")
self.reformulate_sos_constraints()
elif not solver_supports(solver_name, SolverFeature.SOS_CONSTRAINTS):
raise ValueError(
f"Solver {solver_name} does not support SOS constraints. "
"Use reformulate_sos=True or a solver that supports SOS (gurobi, cplex)."
)

try:
solver_class = getattr(solvers, f"{solvers.SolverName(solver_name).name}")
Expand Down
Loading