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

doc/release_notes.rst

@@ -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
 --------------

doc/sos-constraints.rst

@@ -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:**
 
@@ -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:**
 
@@ -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
 ---------------
 

linopy/constants.py

@@ -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):
     """

linopy/io.py

@@ -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:
@@ -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):
@@ -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()

linopy/model.py

@@ -38,6 +38,9 @@
     GREATER_EQUAL,
     HELPER_DIMS,
     LESS_EQUAL,
+    SOS_BIG_M_ATTR,
+    SOS_DIM_ATTR,
+    SOS_TYPE_ATTR,
     TERM_DIM,
     ModelStatus,
     TerminationCondition,
@@ -65,6 +68,7 @@
     IO_APIS,
     available_solvers,
 )
+from linopy.sos_reformulation import reformulate_all_sos
 from linopy.types import (
     ConstantLike,
     ConstraintLike,
@@ -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
@@ -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}"
             )
@@ -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,
@@ -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.
@@ -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]:
         """
@@ -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.
 
@@ -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}")