PyPSA
diff --git a/‎doc/release_notes.rst‎
Lines changed: 3 additions & 0 deletions b/‎doc/release_notes.rst‎
Lines changed: 3 additions & 0 deletions
diff --git a/‎doc/sos-constraints.rst‎
Lines changed: 80 additions & 1 deletion b/‎doc/sos-constraints.rst‎
Lines changed: 80 additions & 1 deletion
diff --git a/‎linopy/constants.py‎
Lines changed: 5 additions & 0 deletions b/‎linopy/constants.py‎
Lines changed: 5 additions & 0 deletions
diff --git a/‎linopy/io.py‎
Lines changed: 5 additions & 5 deletions b/‎linopy/io.py‎
Lines changed: 5 additions & 5 deletions
diff --git a/‎linopy/model.py‎
Lines changed: 102 additions & 46 deletions b/‎linopy/model.py‎
Lines changed: 102 additions & 46 deletions
@@ -4,6 +4,9 @@ Release Notes
 Upcoming Version
 ----------------
 
+* Add SOS1 and SOS2 reformulations for solvers not supporting them.
+
+
 Version 0.6.4
 --------------
 
 
@@ -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
 ---------------
 
 
@@ -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):
     """
 
@@ -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:
@@ -371,8 +371,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):
@@ -740,8 +740,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()
 
@@ -39,6 +39,9 @@
     GREATER_EQUAL,
     HELPER_DIMS,
     LESS_EQUAL,
+    SOS_BIG_M_ATTR,
+    SOS_DIM_ATTR,
+    SOS_TYPE_ATTR,
     TERM_DIM,
     ModelStatus,
     TerminationCondition,
@@ -66,6 +69,10 @@
     IO_APIS,
     available_solvers,
 )
+from linopy.sos_reformulation import (
+    reformulate_sos_constraints,
+    undo_sos_reformulation,
+)
 from linopy.types import (
     ConstantLike,
     ConstraintLike,
@@ -591,6 +598,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
@@ -604,15 +612,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}"
             )
@@ -624,7 +643,13 @@ def add_sos_constraints(
                 f"but got {variable.coords[sos_dim].dtype}"
             )
 
-        variable.attrs.update(sos_type=sos_type, sos_dim=sos_dim)
+        attrs_update: dict[str, Any] = {SOS_TYPE_ATTR: sos_type, SOS_DIM_ATTR: sos_dim}
+        if big_m is not None:
+            if big_m <= 0:
+                raise ValueError(f"big_m must be positive, got {big_m}")
+            attrs_update[SOS_BIG_M_ATTR] = float(big_m)
+
+        variable.attrs.update(attrs_update)
 
     def add_constraints(
         self,
@@ -891,18 +916,22 @@ 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_ATTR], variable.attrs[SOS_DIM_ATTR]
 
-        del variable.attrs["sos_type"], variable.attrs["sos_dim"]
+        variable.attrs.pop(SOS_BIG_M_ATTR, None)
 
         logger.debug(
             f"Removed sos{sos_type} constraint on {sos_dim} from {variable.name}"
         )
 
+    reformulate_sos_constraints = reformulate_sos_constraints
+
     def remove_objective(self) -> None:
         """
         Remove the objective's linear expression from the model.
@@ -1187,6 +1216,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]:
         """
@@ -1256,6 +1286,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.
 
@@ -1353,11 +1388,25 @@ def solve(
                 f"Solver {solver_name} does not support quadratic problems."
             )
 
-        # 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.")
+        sos_reform_result = None
+        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}")
+                sos_reform_result = reformulate_sos_constraints(self)
+            elif reformulate_sos and solver_supports(
+                solver_name, SolverFeature.SOS_CONSTRAINTS
+            ):
+                logger.warning(
+                    f"Solver {solver_name} supports SOS natively; "
+                    "reformulate_sos=True is ignored."
+                )
+            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}")
@@ -1406,44 +1455,51 @@ def solve(
                 if fn is not None and (os.path.exists(fn) and not keep_files):
                     os.remove(fn)
 
-        result.info()
-
-        self.objective._value = result.solution.objective
-        self.status = result.status.status.value
-        self.termination_condition = result.status.termination_condition.value
-        self.solver_model = result.solver_model
-        self.solver_name = solver_name
-
-        if not result.status.is_ok:
-            return result.status.status.value, result.status.termination_condition.value
+        try:
+            result.info()
+
+            self.objective._value = result.solution.objective
+            self.status = result.status.status.value
+            self.termination_condition = result.status.termination_condition.value
+            self.solver_model = result.solver_model
+            self.solver_name = solver_name
+
+            if not result.status.is_ok:
+                return (
+                    result.status.status.value,
+                    result.status.termination_condition.value,
+                )
 
-        # map solution and dual to original shape which includes missing values
-        sol = result.solution.primal.copy()
-        sol = set_int_index(sol)
-        sol.loc[-1] = nan
+            # map solution and dual to original shape which includes missing values
+            sol = result.solution.primal.copy()
+            sol = set_int_index(sol)
+            sol.loc[-1] = nan
 
-        for name, var in self.variables.items():
-            idx = np.ravel(var.labels)
-            try:
-                vals = sol[idx].values.reshape(var.labels.shape)
-            except KeyError:
-                vals = sol.reindex(idx).values.reshape(var.labels.shape)
-            var.solution = xr.DataArray(vals, var.coords)
-
-        if not result.solution.dual.empty:
-            dual = result.solution.dual.copy()
-            dual = set_int_index(dual)
-            dual.loc[-1] = nan
-
-            for name, con in self.constraints.items():
-                idx = np.ravel(con.labels)
+            for name, var in self.variables.items():
+                idx = np.ravel(var.labels)
                 try:
-                    vals = dual[idx].values.reshape(con.labels.shape)
+                    vals = sol[idx].values.reshape(var.labels.shape)
                 except KeyError:
-                    vals = dual.reindex(idx).values.reshape(con.labels.shape)
-                con.dual = xr.DataArray(vals, con.labels.coords)
+                    vals = sol.reindex(idx).values.reshape(var.labels.shape)
+                var.solution = xr.DataArray(vals, var.coords)
+
+            if not result.solution.dual.empty:
+                dual = result.solution.dual.copy()
+                dual = set_int_index(dual)
+                dual.loc[-1] = nan
+
+                for name, con in self.constraints.items():
+                    idx = np.ravel(con.labels)
+                    try:
+                        vals = dual[idx].values.reshape(con.labels.shape)
+                    except KeyError:
+                        vals = dual.reindex(idx).values.reshape(con.labels.shape)
+                    con.dual = xr.DataArray(vals, con.labels.coords)
 
-        return result.status.status.value, result.status.termination_condition.value
+            return result.status.status.value, result.status.termination_condition.value
+        finally:
+            if sos_reform_result is not None:
+                undo_sos_reformulation(self, sos_reform_result)
 
     def _mock_solve(
         self,