From d62762fcab83395738ff7220ff6b65207cf82da4 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 09:07:37 +0100 Subject: [PATCH 01/14] The SOS constraint reformulation feature has been implemented successfully. Here's a summary: Implementation Summary New File: linopy/sos_reformulation.py Core reformulation functions: - validate_bounds_for_reformulation() - Validates that variables have finite bounds - compute_big_m_values() - Computes Big-M values from variable bounds - reformulate_sos1() - Reformulates SOS1 constraints using binary indicators and Big-M constraints - reformulate_sos2() - Reformulates SOS2 constraints using segment indicators and adjacency constraints - reformulate_all_sos() - Reformulates all SOS constraints in a model Modified: linopy/model.py - Added import for reformulate_all_sos - Added reformulate_sos_constraints() method to Model class - Added reformulate_sos: bool = False parameter to solve() method - Updated SOS constraint check to automatically reformulate when reformulate_sos=True and solver doesn't support SOS natively New Test File: test/test_sos_reformulation.py 36 comprehensive tests covering: - Bound validation (finite/infinite) - Big-M computation - SOS1 reformulation (basic, negative bounds, multi-dimensional) - SOS2 reformulation (basic, trivial cases, adjacency) - Integration with solve() and HiGHS - Equivalence with native Gurobi SOS support - Edge cases (zero bounds, multiple SOS, custom prefix) Usage Example m = linopy.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.add_objective(x.sum(), sense='max') # Works with HiGHS (which doesn't support SOS natively) m.solve(solver_name='highs', reformulate_sos=True) --- linopy/model.py | 61 +++- linopy/sos_reformulation.py | 342 ++++++++++++++++++++ test/test_sos_reformulation.py | 571 +++++++++++++++++++++++++++++++++ 3 files changed, 970 insertions(+), 4 deletions(-) create mode 100644 linopy/sos_reformulation.py create mode 100644 test/test_sos_reformulation.py diff --git a/linopy/model.py b/linopy/model.py index 657b2d45..1c343ab9 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -65,6 +65,7 @@ IO_APIS, available_solvers, ) +from linopy.sos_reformulation import reformulate_all_sos from linopy.types import ( ConstantLike, ConstraintLike, @@ -842,6 +843,45 @@ def remove_sos_constraints(self, variable: Variable) -> None: 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). + + Note: This requires all SOS variables to have finite bounds, as the Big-M + values are computed from the variable bounds. + + 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. + + 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'] + """ + return reformulate_all_sos(self, prefix=prefix) + def remove_objective(self) -> None: """ Remove the objective's linear expression from the model. @@ -1126,6 +1166,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 +1236,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 +1339,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}") diff --git a/linopy/sos_reformulation.py b/linopy/sos_reformulation.py new file mode 100644 index 00000000..8fedb0a5 --- /dev/null +++ b/linopy/sos_reformulation.py @@ -0,0 +1,342 @@ +""" +Linopy SOS constraint reformulation module. + +This module provides functions to reformulate SOS1 and SOS2 constraints +as binary + linear constraints for solvers that don't support them natively. +""" + +from __future__ import annotations + +import logging +from typing import TYPE_CHECKING + +import numpy as np +import pandas as pd +from xarray import DataArray + +if TYPE_CHECKING: + from linopy.model import Model + from linopy.variables import Variable + +logger = logging.getLogger(__name__) + + +def validate_bounds_for_reformulation(var: Variable) -> None: + """ + Validate that a variable has finite bounds required for SOS reformulation. + + Parameters + ---------- + var : Variable + Variable to validate. + + Raises + ------ + ValueError + If any bound is infinite (required for Big-M formulation). + """ + lower = var.lower + upper = var.upper + + if np.isinf(lower).any(): + raise ValueError( + f"Variable '{var.name}' has infinite lower bounds. " + "Finite bounds are required for SOS reformulation (Big-M method)." + ) + if np.isinf(upper).any(): + raise ValueError( + f"Variable '{var.name}' has infinite upper bounds. " + "Finite bounds are required for SOS reformulation (Big-M method)." + ) + + +def compute_big_m_values(var: Variable) -> tuple[DataArray, DataArray]: + """ + Compute Big-M values from variable bounds. + + Parameters + ---------- + var : Variable + Variable with finite bounds. + + Returns + ------- + tuple[DataArray, DataArray] + (M_upper, M_lower) - Big-M values computed from bounds. + M_upper = upper bound (for x <= U * y constraints) + M_lower = lower bound (for x >= L * y constraints) + """ + return var.upper, var.lower + + +def reformulate_sos1(model: Model, var: Variable, prefix: str) -> None: + """ + Reformulate an SOS1 constraint as binary + linear constraints. + + SOS1: At most one variable can be non-zero. + + Reformulation: + For each x[i] with bounds [L[i], U[i]]: + - Add binary y[i] + - x[i] <= U[i] * y[i] (if U[i] > 0) + - x[i] >= L[i] * y[i] (if L[i] < 0) + - sum(y) <= 1 + + Parameters + ---------- + model : Model + The model to add reformulation constraints to. + var : Variable + The variable with SOS1 constraint. + prefix : str + Prefix for naming auxiliary variables and constraints. + """ + sos_dim = var.attrs["sos_dim"] + var_name = var.name + + # Get bounds + M_upper, M_lower = compute_big_m_values(var) + + # Extract coords as list of DataArrays for add_variables + coords_list = [var.coords[d] for d in var.dims] + + # Create binary indicator variables with same dimensions as original variable + y_name = f"{prefix}{var_name}_y" + y = model.add_variables( + coords=coords_list, + name=y_name, + binary=True, + ) + + # Add upper bound constraints: x <= U * y (when U > 0) + # This ensures x can only be positive when y = 1 + has_positive_upper = (M_upper > 0).any() + if has_positive_upper: + model.add_constraints( + var <= M_upper * y, + name=f"{prefix}{var_name}_upper", + ) + + # Add lower bound constraints: x >= L * y (when L < 0) + # This ensures x can only be negative when y = 1 + has_negative_lower = (M_lower < 0).any() + if has_negative_lower: + model.add_constraints( + var >= M_lower * y, + name=f"{prefix}{var_name}_lower", + ) + + # Add cardinality constraint: sum(y) <= 1 over the SOS dimension + # This is summed over sos_dim, keeping other dimensions + model.add_constraints( + y.sum(dim=sos_dim) <= 1, + name=f"{prefix}{var_name}_card", + ) + + logger.debug(f"Reformulated SOS1 constraint for variable '{var_name}'") + + +def reformulate_sos2(model: Model, var: Variable, prefix: str) -> None: + """ + Reformulate an SOS2 constraint as binary + linear constraints. + + SOS2: At most two adjacent variables can be non-zero. + + Reformulation: + For ordered x[i], i = 0..n-1: + - Add binary z[i] for i = 0..n-2 (segment indicators) + - x[0] <= U[0] * z[0] + - x[i] <= U[i] * (z[i-1] + z[i]) for i = 1..n-2 + - x[n-1] <= U[n-1] * z[n-2] + - Similar for lower bounds if L[i] < 0 + - sum(z) <= 1 + + Parameters + ---------- + model : Model + The model to add reformulation constraints to. + var : Variable + The variable with SOS2 constraint. + prefix : str + Prefix for naming auxiliary variables and constraints. + """ + sos_dim = var.attrs["sos_dim"] + var_name = var.name + + # Get the size of the SOS dimension + n = var.sizes[sos_dim] + + # Trivial case: single element always satisfies SOS2 + if n <= 1: + logger.debug( + f"Skipping SOS2 reformulation for '{var_name}' (trivial case: n={n})" + ) + return + + # Get bounds + M_upper, M_lower = compute_big_m_values(var) + + # Create n-1 binary segment indicators + # z[i] indicates that the "active segment" is between positions i and i+1 + segment_coords_values = var.coords[sos_dim].values[:-1] # n-1 segment indicators + + # Build coords for z: same as var but with truncated sos_dim + z_coords_list = [] + for d in var.dims: + if d == sos_dim: + z_coords_list.append(pd.Index(segment_coords_values, name=sos_dim)) + else: + z_coords_list.append(var.coords[d]) + + z_name = f"{prefix}{var_name}_z" + z = model.add_variables( + coords=z_coords_list, + name=z_name, + binary=True, + ) + + # For SOS2, each x[i] can only be non-zero if an adjacent segment is active + # x[0] needs z[0] active + # x[i] (1 <= i <= n-2) needs z[i-1] or z[i] active + # x[n-1] needs z[n-2] active + + # Convert to LinearExpression to avoid SOS attribute validation issues during isel + # (Variable's isel preserves SOS attrs, which fail validation when dim is removed) + x_expr = 1 * var + z_expr = 1 * z + + # Process upper bound constraints (for positive bounds) + has_positive_upper = (M_upper > 0).any() + if has_positive_upper: + # First element: x[0] <= U[0] * z[0] + x_first = x_expr.isel({sos_dim: 0}) + z_first = z_expr.isel({sos_dim: 0}) + U_first = M_upper.isel({sos_dim: 0}) + model.add_constraints( + x_first <= U_first * z_first, + name=f"{prefix}{var_name}_upper_first", + ) + + # Middle elements: x[i] <= U[i] * (z[i-1] + z[i]) for i = 1..n-2 + if n > 2: + for i in range(1, n - 1): + x_i = x_expr.isel({sos_dim: i}) + z_prev = z_expr.isel({sos_dim: i - 1}) + z_curr = z_expr.isel({sos_dim: i}) + U_i = M_upper.isel({sos_dim: i}) + model.add_constraints( + x_i <= U_i * (z_prev + z_curr), + name=f"{prefix}{var_name}_upper_mid_{i}", + ) + + # Last element: x[n-1] <= U[n-1] * z[n-2] + x_last = x_expr.isel({sos_dim: n - 1}) + z_last = z_expr.isel({sos_dim: n - 2}) + U_last = M_upper.isel({sos_dim: n - 1}) + model.add_constraints( + x_last <= U_last * z_last, + name=f"{prefix}{var_name}_upper_last", + ) + + # Process lower bound constraints (for negative bounds) + has_negative_lower = (M_lower < 0).any() + if has_negative_lower: + # First element: x[0] >= L[0] * z[0] + x_first = x_expr.isel({sos_dim: 0}) + z_first = z_expr.isel({sos_dim: 0}) + L_first = M_lower.isel({sos_dim: 0}) + model.add_constraints( + x_first >= L_first * z_first, + name=f"{prefix}{var_name}_lower_first", + ) + + # Middle elements: x[i] >= L[i] * (z[i-1] + z[i]) for i = 1..n-2 + if n > 2: + for i in range(1, n - 1): + x_i = x_expr.isel({sos_dim: i}) + z_prev = z_expr.isel({sos_dim: i - 1}) + z_curr = z_expr.isel({sos_dim: i}) + L_i = M_lower.isel({sos_dim: i}) + model.add_constraints( + x_i >= L_i * (z_prev + z_curr), + name=f"{prefix}{var_name}_lower_mid_{i}", + ) + + # Last element: x[n-1] >= L[n-1] * z[n-2] + x_last = x_expr.isel({sos_dim: n - 1}) + z_last = z_expr.isel({sos_dim: n - 2}) + L_last = M_lower.isel({sos_dim: n - 1}) + model.add_constraints( + x_last >= L_last * z_last, + name=f"{prefix}{var_name}_lower_last", + ) + + # Add cardinality constraint: sum(z) <= 1 + model.add_constraints( + z.sum(dim=sos_dim) <= 1, + name=f"{prefix}{var_name}_card", + ) + + logger.debug(f"Reformulated SOS2 constraint for variable '{var_name}'") + + +def reformulate_all_sos(model: Model, prefix: str = "_sos_reform_") -> list[str]: + """ + Reformulate all SOS constraints in the model. + + Parameters + ---------- + model : Model + The model containing SOS constraints to reformulate. + prefix : str, optional + Prefix for naming auxiliary variables and constraints. + Default is "_sos_reform_". + + Returns + ------- + list[str] + List of variable names that were reformulated. + """ + reformulated_vars = [] + + # Get all SOS variables + sos_vars = list(model.variables.sos) + + for var_name in sos_vars: + var = model.variables[var_name] + sos_type = var.attrs.get("sos_type") + sos_dim = var.attrs.get("sos_dim") + + if sos_type is None or sos_dim is None: + continue + + # Skip single-element SOS (trivially satisfied) + if var.sizes[sos_dim] <= 1: + logger.debug( + f"Skipping SOS{sos_type} reformulation for '{var_name}' (single element)" + ) + continue + + # Validate bounds + validate_bounds_for_reformulation(var) + + # Check if all bounds are zero (variable already fixed to 0) + M_upper, M_lower = compute_big_m_values(var) + if (M_upper == 0).all() and (M_lower == 0).all(): + logger.debug( + f"Skipping SOS{sos_type} reformulation for '{var_name}' (fixed to zero)" + ) + continue + + # Perform reformulation based on SOS type + if sos_type == 1: + reformulate_sos1(model, var, prefix) + elif sos_type == 2: + reformulate_sos2(model, var, prefix) + + # Remove the SOS constraint from the variable + model.remove_sos_constraints(var) + reformulated_vars.append(var_name) + + logger.info(f"Reformulated {len(reformulated_vars)} SOS constraint(s)") + return reformulated_vars diff --git a/test/test_sos_reformulation.py b/test/test_sos_reformulation.py new file mode 100644 index 00000000..57458625 --- /dev/null +++ b/test/test_sos_reformulation.py @@ -0,0 +1,571 @@ +"""Tests for SOS constraint reformulation.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest + +from linopy import Model, available_solvers +from linopy.sos_reformulation import ( + compute_big_m_values, + reformulate_all_sos, + reformulate_sos1, + reformulate_sos2, + validate_bounds_for_reformulation, +) + + +class TestValidateBounds: + """Tests for validate_bounds_for_reformulation.""" + + def test_finite_bounds_pass(self) -> None: + """Finite bounds should pass validation.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=-1, upper=1, coords=[idx], name="x") + # Should not raise + validate_bounds_for_reformulation(x) + + def test_infinite_upper_bounds_raise(self) -> None: + """Infinite upper bounds should raise ValueError.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=np.inf, coords=[idx], name="x") + with pytest.raises(ValueError, match="infinite upper bounds"): + validate_bounds_for_reformulation(x) + + def test_infinite_lower_bounds_raise(self) -> None: + """Infinite lower bounds should raise ValueError.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=-np.inf, upper=1, coords=[idx], name="x") + with pytest.raises(ValueError, match="infinite lower bounds"): + validate_bounds_for_reformulation(x) + + def test_mixed_infinite_bounds_raise(self) -> None: + """Mixed finite/infinite bounds should raise ValueError.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables( + lower=np.array([0, -np.inf, 0]), + upper=np.array([1, 1, 1]), + coords=[idx], + name="x", + ) + with pytest.raises(ValueError, match="infinite lower bounds"): + validate_bounds_for_reformulation(x) + + +class TestComputeBigM: + """Tests for compute_big_m_values.""" + + def test_positive_bounds(self) -> None: + """Test Big-M computation with positive bounds.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=10, coords=[idx], name="x") + M_upper, M_lower = compute_big_m_values(x) + assert np.allclose(M_upper.values, [10, 10, 10]) + assert np.allclose(M_lower.values, [0, 0, 0]) + + def test_negative_bounds(self) -> None: + """Test Big-M computation with negative bounds.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=-10, upper=-1, coords=[idx], name="x") + M_upper, M_lower = compute_big_m_values(x) + assert np.allclose(M_upper.values, [-1, -1, -1]) + assert np.allclose(M_lower.values, [-10, -10, -10]) + + def test_mixed_bounds(self) -> None: + """Test Big-M computation with mixed bounds.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=-5, upper=10, coords=[idx], name="x") + M_upper, M_lower = compute_big_m_values(x) + assert np.allclose(M_upper.values, [10, 10, 10]) + assert np.allclose(M_lower.values, [-5, -5, -5]) + + def test_varying_bounds(self) -> None: + """Test Big-M computation with varying bounds.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables( + lower=np.array([-1, -2, -3]), + upper=np.array([1, 2, 3]), + coords=[idx], + name="x", + ) + M_upper, M_lower = compute_big_m_values(x) + assert np.allclose(M_upper.values, [1, 2, 3]) + assert np.allclose(M_lower.values, [-1, -2, -3]) + + +class TestSOS1Reformulation: + """Tests for SOS1 reformulation.""" + + def test_basic_sos1(self) -> None: + """Test basic SOS1 reformulation.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + # Reformulate + reformulate_sos1(m, x, "_test_") + m.remove_sos_constraints(x) + + # Check auxiliary variables and constraints were added + assert "_test_x_y" in m.variables + assert "_test_x_upper" in m.constraints + assert "_test_x_card" in m.constraints + + # Binary variable should have same dimensions + y = m.variables["_test_x_y"] + assert y.dims == x.dims + assert y.sizes == x.sizes + + def test_sos1_negative_bounds(self) -> None: + """Test SOS1 reformulation with negative bounds.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=-2, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + reformulate_sos1(m, x, "_test_") + m.remove_sos_constraints(x) + + # Should have both upper and lower constraints + assert "_test_x_upper" in m.constraints + assert "_test_x_lower" in m.constraints + assert "_test_x_card" in m.constraints + + def test_sos1_all_negative_bounds(self) -> None: + """Test SOS1 reformulation with all negative bounds (no upper constraint).""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=-2, upper=-1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + reformulate_sos1(m, x, "_test_") + m.remove_sos_constraints(x) + + # Should only have lower constraint (no positive upper bounds) + assert "_test_x_upper" not in m.constraints + assert "_test_x_lower" in m.constraints + assert "_test_x_card" in m.constraints + + def test_sos1_multidimensional(self) -> None: + """Test SOS1 reformulation with multi-dimensional variables.""" + m = Model() + idx_i = pd.Index([0, 1, 2], name="i") + idx_j = pd.Index([0, 1], name="j") + x = m.add_variables(lower=0, upper=1, coords=[idx_i, idx_j], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + reformulate_sos1(m, x, "_test_") + m.remove_sos_constraints(x) + + # Binary variable should have same dimensions + y = m.variables["_test_x_y"] + assert set(y.dims) == {"i", "j"} + + # Cardinality constraint should have reduced dimensions (summed over i) + card_con = m.constraints["_test_x_card"] + assert "j" in card_con.dims + # The constraint has j dimension preserved + + +class TestSOS2Reformulation: + """Tests for SOS2 reformulation.""" + + def test_basic_sos2(self) -> None: + """Test basic SOS2 reformulation.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=2, sos_dim="i") + + reformulate_sos2(m, x, "_test_") + m.remove_sos_constraints(x) + + # Check auxiliary variables and constraints were added + assert "_test_x_z" in m.variables + assert "_test_x_upper_first" in m.constraints + assert "_test_x_upper_last" in m.constraints + assert "_test_x_card" in m.constraints + + # Segment indicators should have n-1 elements + z = m.variables["_test_x_z"] + assert z.sizes["i"] == 2 # n-1 = 3-1 = 2 + + def test_sos2_trivial_single_element(self) -> None: + """Test SOS2 with single element (trivially satisfied, skipped).""" + m = Model() + idx = pd.Index([0], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=2, sos_dim="i") + + # Should not add anything for trivial case + reformulate_sos2(m, x, "_test_") + + assert "_test_x_z" not in m.variables + + def test_sos2_two_elements(self) -> None: + """Test SOS2 with two elements (one segment indicator).""" + m = Model() + idx = pd.Index([0, 1], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=2, sos_dim="i") + + reformulate_sos2(m, x, "_test_") + m.remove_sos_constraints(x) + + # Should have 1 segment indicator + z = m.variables["_test_x_z"] + assert z.sizes["i"] == 1 + + def test_sos2_with_middle_constraints(self) -> None: + """Test SOS2 with 4+ elements to ensure middle constraints are created.""" + m = Model() + idx = pd.Index([0, 1, 2, 3], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=2, sos_dim="i") + + reformulate_sos2(m, x, "_test_") + m.remove_sos_constraints(x) + + # Should have first, middle, and last constraints + assert "_test_x_upper_first" in m.constraints + assert "_test_x_upper_mid_1" in m.constraints # middle element at index 1 + assert "_test_x_upper_mid_2" in m.constraints # middle element at index 2 + assert "_test_x_upper_last" in m.constraints + + def test_sos2_negative_bounds(self) -> None: + """Test SOS2 reformulation with negative bounds.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=-2, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=2, sos_dim="i") + + reformulate_sos2(m, x, "_test_") + m.remove_sos_constraints(x) + + # Should have both upper and lower constraints + assert "_test_x_upper_first" in m.constraints + assert "_test_x_lower_first" in m.constraints + + +class TestReformulateAllSOS: + """Tests for reformulate_all_sos.""" + + def test_reformulate_single_sos1(self) -> None: + """Test reformulating a single SOS1 variable.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + reformulated = reformulate_all_sos(m) + + assert reformulated == ["x"] + assert len(list(m.variables.sos)) == 0 # No more SOS variables + + def test_reformulate_multiple_sos(self) -> None: + """Test reformulating multiple SOS variables.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + y = m.add_variables(lower=0, upper=2, coords=[idx], name="y") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + m.add_sos_constraints(y, sos_type=2, sos_dim="i") + + reformulated = reformulate_all_sos(m) + + assert set(reformulated) == {"x", "y"} + assert len(list(m.variables.sos)) == 0 + + def test_reformulate_skips_single_element(self) -> None: + """Test that single element SOS is skipped.""" + m = Model() + idx = pd.Index([0], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + reformulated = reformulate_all_sos(m) + + # Single element SOS is skipped (trivially satisfied) + assert reformulated == [] + + def test_reformulate_skips_zero_bounds(self) -> None: + """Test that variables with zero bounds are skipped.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=0, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + reformulated = reformulate_all_sos(m) + + # Zero bounds means variable is fixed to 0, no reformulation needed + assert reformulated == [] + + def test_reformulate_raises_on_infinite_bounds(self) -> None: + """Test that infinite bounds raise ValueError.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=np.inf, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + with pytest.raises(ValueError, match="infinite"): + reformulate_all_sos(m) + + +class TestModelReformulateSOS: + """Tests for Model.reformulate_sos_constraints method.""" + + def test_reformulate_inplace(self) -> None: + """Test inplace reformulation.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + result = m.reformulate_sos_constraints() + + assert result == ["x"] + assert len(list(m.variables.sos)) == 0 + assert "_sos_reform_x_y" in m.variables + + +@pytest.mark.skipif("highs" not in available_solvers, reason="HiGHS not installed") +class TestSolveWithReformulation: + """Tests for solving with SOS reformulation.""" + + def test_sos1_maximize_with_highs(self) -> None: + """Test SOS1 maximize problem with HiGHS using reformulation.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + m.add_objective(x * np.array([1, 2, 3]), sense="max") + + m.solve(solver_name="highs", reformulate_sos=True) + + # Should maximize by choosing x[2] = 1 + assert np.isclose(x.solution.values[2], 1, atol=1e-5) + assert np.isclose(x.solution.values[0], 0, atol=1e-5) + assert np.isclose(x.solution.values[1], 0, atol=1e-5) + assert m.objective.value is not None + assert np.isclose(m.objective.value, 3, atol=1e-5) + + def test_sos1_minimize_with_highs(self) -> None: + """Test SOS1 minimize problem with HiGHS using reformulation.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + m.add_objective(x * np.array([3, 2, 1]), sense="min") + + m.solve(solver_name="highs", reformulate_sos=True) + + # Should minimize to 0 by setting all x = 0 + assert np.isclose(m.objective.value, 0, atol=1e-5) + + def test_sos2_maximize_with_highs(self) -> None: + """Test SOS2 maximize problem with HiGHS using reformulation.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=2, sos_dim="i") + m.add_objective(x * np.array([1, 2, 3]), sense="max") + + m.solve(solver_name="highs", reformulate_sos=True) + + # SOS2 allows two adjacent non-zeros, so x[1] and x[2] can both be 1 + # Maximum is 2 + 3 = 5 + assert m.objective.value is not None + assert np.isclose(m.objective.value, 5, atol=1e-5) + # Check that at most two adjacent variables are non-zero + nonzero_count = (np.abs(x.solution.values) > 1e-5).sum() + assert nonzero_count <= 2 + + def test_sos2_different_coefficients(self) -> None: + """Test SOS2 with different coefficients.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=2, sos_dim="i") + m.add_objective(x * np.array([2, 1, 3]), sense="max") + + m.solve(solver_name="highs", reformulate_sos=True) + + # Best is x[1]=1 and x[2]=1 giving 1+3=4 + # or x[0]=1 and x[1]=1 giving 2+1=3 + assert m.objective.value is not None + assert np.isclose(m.objective.value, 4, atol=1e-5) + + def test_reformulate_sos_false_raises_error(self) -> None: + """Test that HiGHS without reformulate_sos raises error.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + m.add_objective(x.sum(), sense="max") + + with pytest.raises(ValueError, match="does not support SOS"): + m.solve(solver_name="highs", reformulate_sos=False) + + def test_multidimensional_sos1_with_highs(self) -> None: + """Test multi-dimensional SOS1 with HiGHS.""" + m = Model() + idx_i = pd.Index([0, 1, 2], name="i") + idx_j = pd.Index([0, 1], name="j") + x = m.add_variables(lower=0, upper=1, coords=[idx_i, idx_j], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + m.add_objective(x.sum(), sense="max") + + m.solve(solver_name="highs", reformulate_sos=True) + + # For each j, at most one x[i, j] can be non-zero + # Maximum is achieved by one non-zero per j column: 2 total + assert m.objective.value is not None + assert np.isclose(m.objective.value, 2, atol=1e-5) + + # Check SOS1 is satisfied for each j + for j in idx_j: + nonzero_count = (np.abs(x.solution.sel(j=j).values) > 1e-5).sum() + assert nonzero_count <= 1 + + +@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") +class TestEquivalenceWithGurobi: + """Tests comparing reformulated solutions with native Gurobi SOS.""" + + def test_sos1_equivalence(self) -> None: + """Test that reformulated SOS1 gives same result as native Gurobi.""" + gurobipy = pytest.importorskip("gurobipy") + + # Native Gurobi solution + m1 = Model() + idx = pd.Index([0, 1, 2], name="i") + x1 = m1.add_variables(lower=0, upper=1, coords=[idx], name="x") + m1.add_sos_constraints(x1, sos_type=1, sos_dim="i") + m1.add_objective(x1 * np.array([1, 2, 3]), sense="max") + + try: + m1.solve(solver_name="gurobi") + except gurobipy.GurobiError as exc: + pytest.skip(f"Gurobi environment unavailable: {exc}") + + # Reformulated solution with HiGHS + m2 = Model() + x2 = m2.add_variables(lower=0, upper=1, coords=[idx], name="x") + m2.add_sos_constraints(x2, sos_type=1, sos_dim="i") + m2.add_objective(x2 * np.array([1, 2, 3]), sense="max") + + if "highs" in available_solvers: + m2.solve(solver_name="highs", reformulate_sos=True) + assert np.isclose(m1.objective.value, m2.objective.value, atol=1e-5) + + def test_sos2_equivalence(self) -> None: + """Test that reformulated SOS2 gives same result as native Gurobi.""" + gurobipy = pytest.importorskip("gurobipy") + + # Native Gurobi solution + m1 = Model() + idx = pd.Index([0, 1, 2], name="i") + x1 = m1.add_variables(lower=0, upper=1, coords=[idx], name="x") + m1.add_sos_constraints(x1, sos_type=2, sos_dim="i") + m1.add_objective(x1 * np.array([1, 2, 3]), sense="max") + + try: + m1.solve(solver_name="gurobi") + except gurobipy.GurobiError as exc: + pytest.skip(f"Gurobi environment unavailable: {exc}") + + # Reformulated solution with HiGHS + m2 = Model() + x2 = m2.add_variables(lower=0, upper=1, coords=[idx], name="x") + m2.add_sos_constraints(x2, sos_type=2, sos_dim="i") + m2.add_objective(x2 * np.array([1, 2, 3]), sense="max") + + if "highs" in available_solvers: + m2.solve(solver_name="highs", reformulate_sos=True) + assert np.isclose(m1.objective.value, m2.objective.value, atol=1e-5) + + +class TestEdgeCases: + """Tests for edge cases.""" + + def test_zero_lower_bound(self) -> None: + """Test handling of zero lower bound.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + reformulate_all_sos(m) + + # Should only have upper constraints (no negative lower bounds) + assert "_sos_reform_x_upper" in m.constraints + assert "_sos_reform_x_lower" not in m.constraints + + def test_preserves_non_sos_variables(self) -> None: + """Test that non-SOS variables are preserved.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_variables(lower=0, upper=2, coords=[idx], name="y") # No SOS + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + reformulate_all_sos(m) + + # y should be unchanged + assert "y" in m.variables + assert "sos_type" not in m.variables["y"].attrs + + def test_custom_prefix(self) -> None: + """Test custom prefix for reformulation.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + reformulate_all_sos(m, prefix="_custom_") + + assert "_custom_x_y" in m.variables + assert "_custom_x_upper" in m.constraints + assert "_custom_x_card" in m.constraints + + def test_constraints_with_sos_variables(self) -> None: + """Test that existing constraints with SOS variables work after reformulation.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + y = m.add_variables(lower=0, upper=10, name="y") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + # Add constraint involving SOS variable + m.add_constraints(x.sum() <= y, name="linking") + + # Reformulate + reformulate_all_sos(m) + + # Original constraint should still exist + assert "linking" in m.constraints + + def test_float_coordinates(self) -> None: + """Test SOS with float coordinates (common for piecewise linear).""" + m = Model() + breakpoints = pd.Index([0.0, 0.5, 1.0], name="bp") + lambdas = m.add_variables(lower=0, upper=1, coords=[breakpoints], name="lambda") + m.add_sos_constraints(lambdas, sos_type=2, sos_dim="bp") + + reformulate_all_sos(m) + + # Should work with float coordinates + assert "_sos_reform_lambda_z" in m.variables + z = m.variables["_sos_reform_lambda_z"] + # Segment indicators have n-1 = 2 elements + assert z.sizes["bp"] == 2 From 80d57527867b5b121912fbf1ab5d36c353c155cb Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 09:15:40 +0100 Subject: [PATCH 02/14] Documentation Summary MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit New Section: "SOS Reformulation for Unsupported Solvers" Added a comprehensive section (~300 lines) covering: 1. Enabling Reformulation - Shows reformulate_sos=True parameter and manual reformulate_sos_constraints() method 2. Requirements - Explains finite bounds requirement for Big-M method 3. Mathematical Formulation - Clear LaTeX math for both: - SOS1: Binary indicators y_i, upper/lower linking constraints, cardinality constraint - SOS2: Segment indicators z_j, first/middle/last element constraints, cardinality constraint 4. Interpretation - Explains how the constraints work intuitively with examples 5. Auxiliary Variables and Constraints - Documents the naming convention (_sos_reform_ prefix) 6. Multi-dimensional Variables - Shows how broadcasting works 7. Edge Cases Table - Lists all handled edge cases (single-element, zero bounds, all-positive, etc.) 8. Performance Considerations - Trade-offs between native SOS and reformulation 9. Complete Example - Piecewise linear approximation of x² with HiGHS 10. API Reference - Added method signatures for: - Model.add_sos_constraints() - Model.remove_sos_constraints() - Model.reformulate_sos_constraints() - Variables.sos property --- doc/sos-constraints.rst | 275 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 275 insertions(+) diff --git a/doc/sos-constraints.rst b/doc/sos-constraints.rst index 37dd72d2..3fab97b2 100644 --- a/doc/sos-constraints.rst +++ b/doc/sos-constraints.rst @@ -254,6 +254,246 @@ SOS constraints are supported by most modern mixed-integer programming solvers t - MOSEK - MindOpt +For these solvers, linopy provides automatic reformulation (see :ref:`sos-reformulation` below). + +.. _sos-reformulation: + +SOS Reformulation for Unsupported Solvers +----------------------------------------- + +Linopy can automatically reformulate SOS constraints as binary + linear constraints +using the Big-M method. This allows you to use SOS constraints with solvers that +don't support them natively (HiGHS, GLPK, MOSEK, etc.). + +Enabling Reformulation +~~~~~~~~~~~~~~~~~~~~~~ + +Pass ``reformulate_sos=True`` to the ``solve()`` method: + +.. code-block:: python + + import linopy + import pandas as pd + + m = linopy.Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + m.add_objective(x.sum(), sense="max") + + # Now works with HiGHS! + m.solve(solver_name="highs", reformulate_sos=True) + +You can also reformulate manually before solving: + +.. code-block:: python + + # Reformulate in place + reformulated_vars = m.reformulate_sos_constraints() + print(f"Reformulated: {reformulated_vars}") + + # Then solve with any solver + m.solve(solver_name="highs") + +Requirements +~~~~~~~~~~~~ + +**Finite bounds are required.** The reformulation uses the Big-M method, which +derives M values from the variable bounds. If any SOS variable has infinite +bounds, a ``ValueError`` is raised: + +.. code-block:: python + + # This will raise ValueError + x = m.add_variables(lower=0, upper=np.inf, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + m.solve(solver_name="highs", reformulate_sos=True) + # ValueError: Variable 'x' has infinite upper bounds. + +Mathematical Formulation +~~~~~~~~~~~~~~~~~~~~~~~~ + +The reformulation converts SOS constraints into Mixed-Integer Linear Programming +(MILP) constraints using binary indicator variables. + +**SOS1 Reformulation** + +Given variables :math:`x_i` for :math:`i \in I` with bounds :math:`L_i \leq x_i \leq U_i`, +the SOS1 constraint "at most one :math:`x_i` is non-zero" is reformulated as: + +.. math:: + + \text{Binary indicators:} \quad & y_i \in \{0, 1\} \quad \forall i \in I \\[0.5em] + \text{Upper linking:} \quad & x_i \leq U_i \cdot y_i \quad \forall i \in I \text{ where } U_i > 0 \\[0.5em] + \text{Lower linking:} \quad & x_i \geq L_i \cdot y_i \quad \forall i \in I \text{ where } L_i < 0 \\[0.5em] + \text{Cardinality:} \quad & \sum_{i \in I} y_i \leq 1 + +**Interpretation:** + +- :math:`y_i = 1` means variable :math:`x_i` is "selected" (allowed to be non-zero) +- :math:`y_i = 0` forces :math:`x_i = 0` via the linking constraints +- The cardinality constraint ensures at most one :math:`y_i = 1` + +**Example:** For :math:`x \in [0, 10]`: + +- If :math:`y = 0`: constraint :math:`x \leq 10 \cdot 0 = 0` forces :math:`x = 0` +- If :math:`y = 1`: constraint :math:`x \leq 10 \cdot 1 = 10` allows :math:`x \in [0, 10]` + +**SOS2 Reformulation** + +Given ordered variables :math:`x_0, x_1, \ldots, x_{n-1}` with bounds :math:`L_i \leq x_i \leq U_i`, +the SOS2 constraint "at most two adjacent :math:`x_i` are non-zero" is reformulated using +segment indicators: + +.. math:: + + \text{Segment indicators:} \quad & z_j \in \{0, 1\} \quad \forall j \in \{0, \ldots, n-2\} \\[0.5em] + \text{First variable:} \quad & x_0 \leq U_0 \cdot z_0 \\[0.5em] + \text{Middle variables:} \quad & x_i \leq U_i \cdot (z_{i-1} + z_i) \quad \forall i \in \{1, \ldots, n-2\} \\[0.5em] + \text{Last variable:} \quad & x_{n-1} \leq U_{n-1} \cdot z_{n-2} \\[0.5em] + \text{Cardinality:} \quad & \sum_{j=0}^{n-2} z_j \leq 1 + +(Similar constraints for lower bounds when :math:`L_i < 0`) + +**Interpretation:** + +- :math:`z_j = 1` means "segment :math:`j`" is active (between positions :math:`j` and :math:`j+1`) +- Variable :math:`x_i` can only be non-zero if an adjacent segment is active +- The cardinality constraint ensures at most one segment is active +- This guarantees at most two adjacent variables :math:`(x_j, x_{j+1})` are non-zero + +**Example:** For :math:`n = 4` variables :math:`(x_0, x_1, x_2, x_3)`: + +- If :math:`z_1 = 1` (segment 1 active): :math:`x_1` and :math:`x_2` can be non-zero +- Constraints force :math:`x_0 = 0` (needs :math:`z_0`) and :math:`x_3 = 0` (needs :math:`z_2`) + +Auxiliary Variables and Constraints +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The reformulation creates auxiliary variables and constraints with a ``_sos_reform_`` prefix: + +.. code-block:: python + + m.reformulate_sos_constraints() + + # For SOS1 variable 'x': + # - Binary indicators: _sos_reform_x_y + # - Upper constraints: _sos_reform_x_upper + # - Lower constraints: _sos_reform_x_lower (if L < 0) + # - Cardinality: _sos_reform_x_card + + # For SOS2 variable 'x': + # - Segment indicators: _sos_reform_x_z + # - Upper constraints: _sos_reform_x_upper_first, _sos_reform_x_upper_mid_*, _sos_reform_x_upper_last + # - Lower constraints: _sos_reform_x_lower_first, _sos_reform_x_lower_mid_*, _sos_reform_x_lower_last + # - Cardinality: _sos_reform_x_card + +You can use a custom prefix: + +.. code-block:: python + + m.reformulate_sos_constraints(prefix="_my_sos_") + +Multi-dimensional Variables +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For multi-dimensional variables, the reformulation respects xarray broadcasting. +Constraints are created for each combination of non-SOS dimensions: + +.. code-block:: python + + # 2D variable: periods × options + periods = pd.Index([0, 1], name="periods") + options = pd.Index([0, 1, 2], name="options") + x = m.add_variables(lower=0, upper=1, coords=[periods, options], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="options") + + # After reformulation: + # - Binary y has shape (periods: 2, options: 3) + # - Cardinality constraint: sum over 'options' for each period + # - Result: at most one option selected PER period + +Edge Cases +~~~~~~~~~~ + +The reformulation handles several edge cases automatically: + +.. list-table:: + :header-rows: 1 + :widths: 30 70 + + * - Case + - Handling + * - Single-element SOS + - Skipped (trivially satisfied) + * - All-zero bounds (L=U=0) + - Skipped (variable already fixed to 0) + * - All-positive bounds + - Only upper linking constraints created + * - All-negative bounds + - Only lower linking constraints created + * - Mixed bounds + - Both upper and lower constraints created + +Performance Considerations +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +**Reformulation trade-offs:** + +- **Pros:** Works with any MIP solver; explicit constraints can sometimes help the solver +- **Cons:** Adds binary variables and constraints; may be slower than native SOS support + +**When to use native SOS:** + +- Use native SOS (Gurobi, CPLEX) when available—solvers have specialized branching strategies +- Native SOS often provides better performance for piecewise linear problems + +**When reformulation is useful:** + +- When you must use a solver without SOS support (HiGHS, GLPK) +- For model portability across different solvers +- When debugging—explicit constraints are easier to inspect + +Example: Piecewise Linear with HiGHS +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + import linopy + import pandas as pd + import numpy as np + + # Approximate f(x) = x² using piecewise linear with HiGHS + m = linopy.Model() + + breakpoints = pd.Index([0.0, 1.0, 2.0, 3.0], name="bp") + x_vals = breakpoints.to_numpy() + y_vals = x_vals**2 # [0, 1, 4, 9] + + # SOS2 interpolation weights + lambdas = m.add_variables(lower=0, upper=1, coords=[breakpoints], name="lambdas") + m.add_sos_constraints(lambdas, sos_type=2, sos_dim="bp") + + # Interpolated point + 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).sum(), name="x_interp") + m.add_constraints(y == (lambdas * y_vals).sum(), name="y_interp") + m.add_constraints(x >= 1.5, name="x_min") + + # Minimize y (approximated x²) + m.add_objective(y) + + # Solve with HiGHS using reformulation + m.solve(solver_name="highs", reformulate_sos=True) + + print(f"x = {x.solution.item():.2f}") + print(f"y ≈ x² = {y.solution.item():.2f}") + # Output: x = 1.50, y ≈ x² = 2.50 (exact: 2.25) + Common Patterns --------------- @@ -309,3 +549,38 @@ See Also - :doc:`creating-variables`: Creating variables with coordinates - :doc:`creating-constraints`: Adding regular constraints - :doc:`user-guide`: General linopy usage patterns + +API Reference +------------- + +.. py:method:: Model.add_sos_constraints(variable, sos_type, sos_dim) + + Add an SOS1 or SOS2 constraint for one dimension of a variable. + + :param variable: Variable to constrain + :type variable: Variable + :param sos_type: Type of SOS constraint (1 or 2) + :type sos_type: Literal[1, 2] + :param sos_dim: Dimension along which to apply the SOS constraint + :type sos_dim: str + +.. py:method:: Model.remove_sos_constraints(variable) + + Remove SOS constraints from a variable. + + :param variable: Variable from which to remove SOS constraints + :type variable: Variable + +.. py:method:: Model.reformulate_sos_constraints(prefix="_sos_reform_") + + Reformulate SOS constraints as binary + linear constraints using the Big-M method. + + :param prefix: Prefix for auxiliary variable and constraint names + :type prefix: str + :returns: List of variable names that were reformulated + :rtype: list[str] + :raises ValueError: If any SOS variable has infinite bounds + +.. py:attribute:: Variables.sos + + Property returning all variables with SOS constraints. From 3f6db066fa8413c0d0a6e892cb593f3d67731083 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 09:21:58 +0100 Subject: [PATCH 03/14] Added Tests for Multi-dimensional SOS Unit Tests - test_sos2_multidimensional: Tests that SOS2 reformulation with multi-dimensional variables (i, j) correctly creates: - Segment indicators z with shape (i: n-1, j: m) - Cardinality constraint preserves the j dimension Integration Tests - test_multidimensional_sos2_with_highs: Solves a multi-dimensional SOS2 problem with HiGHS and verifies: - Optimal objective value (4 total - two adjacent non-zeros per column) - SOS2 constraint satisfied for each j: at most 2 non-zeros, and if 2, they're adjacent Test Results test_sos1_multidimensional PASSED test_sos2_multidimensional PASSED test_multidimensional_sos1_with_highs PASSED test_multidimensional_sos2_with_highs PASSED The implementation correctly handles multi-dimensional variables by leveraging xarray's broadcasting - the SOS constraint is applied along the sos_dim for each combination of the other dimensions. --- test/test_sos_reformulation.py | 47 ++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/test/test_sos_reformulation.py b/test/test_sos_reformulation.py index 57458625..0103e44b 100644 --- a/test/test_sos_reformulation.py +++ b/test/test_sos_reformulation.py @@ -256,6 +256,27 @@ def test_sos2_negative_bounds(self) -> None: assert "_test_x_upper_first" in m.constraints assert "_test_x_lower_first" in m.constraints + def test_sos2_multidimensional(self) -> None: + """Test SOS2 reformulation with multi-dimensional variables.""" + m = Model() + idx_i = pd.Index([0, 1, 2], name="i") + idx_j = pd.Index([0, 1], name="j") + x = m.add_variables(lower=0, upper=1, coords=[idx_i, idx_j], name="x") + m.add_sos_constraints(x, sos_type=2, sos_dim="i") + + reformulate_sos2(m, x, "_test_") + m.remove_sos_constraints(x) + + # Segment indicator should have (n-1) elements in i dimension, same j dimension + z = m.variables["_test_x_z"] + assert set(z.dims) == {"i", "j"} + assert z.sizes["i"] == 2 # n-1 = 3-1 = 2 + assert z.sizes["j"] == 2 + + # Cardinality constraint should have j dimension preserved + card_con = m.constraints["_test_x_card"] + assert "j" in card_con.dims + class TestReformulateAllSOS: """Tests for reformulate_all_sos.""" @@ -437,6 +458,32 @@ def test_multidimensional_sos1_with_highs(self) -> None: nonzero_count = (np.abs(x.solution.sel(j=j).values) > 1e-5).sum() assert nonzero_count <= 1 + def test_multidimensional_sos2_with_highs(self) -> None: + """Test multi-dimensional SOS2 with HiGHS.""" + m = Model() + idx_i = pd.Index([0, 1, 2], name="i") + idx_j = pd.Index([0, 1], name="j") + x = m.add_variables(lower=0, upper=1, coords=[idx_i, idx_j], name="x") + m.add_sos_constraints(x, sos_type=2, sos_dim="i") + m.add_objective(x.sum(), sense="max") + + m.solve(solver_name="highs", reformulate_sos=True) + + # For each j, at most two adjacent x[i, j] can be non-zero + # Maximum is achieved by two adjacent non-zeros per j column: 4 total + assert m.objective.value is not None + assert np.isclose(m.objective.value, 4, atol=1e-5) + + # Check SOS2 is satisfied for each j + for j in idx_j: + sol_j = x.solution.sel(j=j).values + nonzero_indices = np.where(np.abs(sol_j) > 1e-5)[0] + # At most 2 non-zeros + assert len(nonzero_indices) <= 2 + # If 2 non-zeros, they must be adjacent + if len(nonzero_indices) == 2: + assert abs(nonzero_indices[1] - nonzero_indices[0]) == 1 + @pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") class TestEquivalenceWithGurobi: From a906e71c34b042c91f98ff75f746f59c27898abb Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 09:50:02 +0100 Subject: [PATCH 04/14] Add custom big_m parameter for SOS reformulation Allow users to specify custom Big-M values in add_sos_constraints() for tighter LP relaxations when variable bounds are conservative. - Add big_m parameter: scalar or tuple(upper, lower) - Store as variable attrs (big_m_upper, big_m_lower) - Skip bound validation when custom big_m provided - Scalar-only design ensures NetCDF persistence works correctly For per-element Big-M values, users should adjust variable bounds directly. --- doc/sos-constraints.rst | 69 +++++++++++++++++++--- linopy/model.py | 48 +++++++++++++-- linopy/sos_reformulation.py | 70 ++++++++++++++++------ test/test_sos_reformulation.py | 104 +++++++++++++++++++++++++++++++++ 4 files changed, 259 insertions(+), 32 deletions(-) diff --git a/doc/sos-constraints.rst b/doc/sos-constraints.rst index 3fab97b2..129ded3f 100644 --- a/doc/sos-constraints.rst +++ b/doc/sos-constraints.rst @@ -298,17 +298,61 @@ You can also reformulate manually before solving: Requirements ~~~~~~~~~~~~ -**Finite bounds are required.** The reformulation uses the Big-M method, which -derives M values from the variable bounds. If any SOS variable has infinite -bounds, a ``ValueError`` is raised: +**Finite bounds or custom Big-M required.** The reformulation uses the Big-M method. +By default, Big-M values are derived from variable bounds. If any SOS variable has +infinite bounds, you must either: + +1. Set finite bounds on the variable, or +2. Specify custom Big-M values via the ``big_m`` parameter .. code-block:: python - # This will raise ValueError - x = m.add_variables(lower=0, upper=np.inf, coords=[idx], name="x") + # Option 1: Finite bounds (default Big-M = bounds) + x = m.add_variables(lower=0, upper=100, coords=[idx], name="x") m.add_sos_constraints(x, sos_type=1, sos_dim="i") - m.solve(solver_name="highs", reformulate_sos=True) - # ValueError: Variable 'x' has infinite upper bounds. + + # Option 2: Custom Big-M (allows infinite bounds) + 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) + +Custom Big-M Values +~~~~~~~~~~~~~~~~~~~ + +The ``big_m`` parameter in ``add_sos_constraints()`` allows you to specify tighter +Big-M values than the variable bounds. This is useful when: + +- Variable bounds are conservatively large (e.g., ``upper=1e6``) +- Other constraints in your model imply tighter effective bounds +- You have domain knowledge about the actual feasible range + +**Syntax options:** + +.. code-block:: python + + # Scalar: symmetric Big-M (M_upper = 10, M_lower = -10) + m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=10) + + # Tuple of scalars: asymmetric (M_upper = 10, M_lower = -5) + m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=(10, -5)) + +**Why tighter Big-M matters:** + +Tighter Big-M values lead to: + +- Better LP relaxation bounds (closer to optimal integer solution) +- Fewer branch-and-bound nodes +- Faster solve times + +.. code-block:: python + + # Example: Variable has large bounds but we know effective range is smaller + x = m.add_variables(lower=0, upper=1000, coords=[idx], name="x") + + # Other constraints limit x to [0, 10] in practice + m.add_constraints(x <= 10) + + # Use tight Big-M for better performance + m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=10) Mathematical Formulation ~~~~~~~~~~~~~~~~~~~~~~~~ @@ -553,7 +597,7 @@ See Also API Reference ------------- -.. py:method:: Model.add_sos_constraints(variable, sos_type, sos_dim) +.. py:method:: Model.add_sos_constraints(variable, sos_type, sos_dim, big_m=None) Add an SOS1 or SOS2 constraint for one dimension of a variable. @@ -563,6 +607,13 @@ API Reference :type sos_type: Literal[1, 2] :param sos_dim: Dimension along which to apply the SOS constraint :type sos_dim: str + :param big_m: Custom Big-M value(s) for reformulation. Can be: + + - ``None`` (default): Use variable bounds + - ``float``: Symmetric Big-M (upper = big_m, lower = -big_m) + - ``tuple[float, float]``: Asymmetric (upper, lower) + + :type big_m: float | tuple[float, float] | None .. py:method:: Model.remove_sos_constraints(variable) @@ -579,7 +630,7 @@ API Reference :type prefix: str :returns: List of variable names that were reformulated :rtype: list[str] - :raises ValueError: If any SOS variable has infinite bounds + :raises ValueError: If any SOS variable has infinite bounds and no custom big_m .. py:attribute:: Variables.sos diff --git a/linopy/model.py b/linopy/model.py index 1c343ab9..d8bf80f3 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -557,6 +557,7 @@ def add_sos_constraints( variable: Variable, sos_type: Literal[1, 2], sos_dim: str, + big_m: float | tuple[float, float] | None = None, ) -> None: """ Add an sos1 or sos2 constraint for one dimension of a variable @@ -570,6 +571,21 @@ def add_sos_constraints( Type of SOS sos_dim : str Which dimension of variable to add SOS constraint to + big_m : float | tuple[float, float] | None, optional + Big-M value(s) for SOS reformulation. Only used when reformulating + SOS constraints for solvers that don't support them natively. + + - None (default): Use variable bounds as Big-M (tightest valid choice) + - float: Symmetric Big-M (M_upper = big_m, M_lower = -big_m) + - tuple (upper, lower): Asymmetric Big-M values + + The Big-M values are used in constraints: + - x <= M_upper * y (upper linking) + - x >= M_lower * y (lower linking, M_lower should be negative) + + Tighter Big-M values improve LP relaxation quality and solve time. + If you know tighter effective bounds from other constraints in your + model, specifying them here can significantly improve performance. """ if sos_type not in (1, 2): raise ValueError(f"sos_type must be 1 or 2, got {sos_type}") @@ -590,7 +606,22 @@ 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 values + attrs_update: dict[str, Any] = {"sos_type": sos_type, "sos_dim": sos_dim} + if big_m is not None: + if isinstance(big_m, tuple): + if len(big_m) != 2: + raise ValueError( + "big_m tuple must have exactly 2 elements (upper, lower)" + ) + attrs_update["big_m_upper"] = float(big_m[0]) + attrs_update["big_m_lower"] = float(big_m[1]) + else: + # Scalar: symmetric (upper = big_m, lower = -big_m) + attrs_update["big_m_upper"] = float(big_m) + attrs_update["big_m_lower"] = -float(big_m) + + variable.attrs.update(attrs_update) def add_constraints( self, @@ -839,6 +870,10 @@ def remove_sos_constraints(self, variable: Variable) -> None: del variable.attrs["sos_type"], variable.attrs["sos_dim"] + # Also remove big_m attributes if present + variable.attrs.pop("big_m_upper", None) + variable.attrs.pop("big_m_lower", None) + logger.debug( f"Removed sos{sos_type} constraint on {sos_dim} from {variable.name}" ) @@ -851,8 +886,9 @@ def reformulate_sos_constraints(self, prefix: str = "_sos_reform_") -> list[str] 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). - Note: This requires all SOS variables to have finite bounds, as the Big-M - values are computed from the variable bounds. + 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 ---------- @@ -868,7 +904,7 @@ def reformulate_sos_constraints(self, prefix: str = "_sos_reform_") -> list[str] Raises ------ ValueError - If any SOS variable has infinite bounds. + If any SOS variable has infinite bounds and no custom big_m was specified. Examples -------- @@ -879,6 +915,10 @@ def reformulate_sos_constraints(self, prefix: str = "_sos_reform_") -> list[str] >>> 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) diff --git a/linopy/sos_reformulation.py b/linopy/sos_reformulation.py index 8fedb0a5..06bfc7eb 100644 --- a/linopy/sos_reformulation.py +++ b/linopy/sos_reformulation.py @@ -25,6 +25,9 @@ def validate_bounds_for_reformulation(var: Variable) -> None: """ Validate that a variable has finite bounds required for SOS reformulation. + If custom big_m values are provided via variable attributes, bounds + validation is skipped for those dimensions. + Parameters ---------- var : Variable @@ -33,40 +36,69 @@ def validate_bounds_for_reformulation(var: Variable) -> None: Raises ------ ValueError - If any bound is infinite (required for Big-M formulation). + If any bound is infinite and no custom big_m is provided. """ - lower = var.lower - upper = var.upper + # If custom big_m is provided, skip bound validation + has_custom_upper = var.attrs.get("big_m_upper") is not None + has_custom_lower = var.attrs.get("big_m_lower") is not None + + if not has_custom_lower: + lower = var.lower + if np.isinf(lower).any(): + raise ValueError( + f"Variable '{var.name}' has infinite lower bounds. " + "Finite bounds are required for SOS reformulation (Big-M method). " + "Alternatively, specify big_m in add_sos_constraints()." + ) - if np.isinf(lower).any(): - raise ValueError( - f"Variable '{var.name}' has infinite lower bounds. " - "Finite bounds are required for SOS reformulation (Big-M method)." - ) - if np.isinf(upper).any(): - raise ValueError( - f"Variable '{var.name}' has infinite upper bounds. " - "Finite bounds are required for SOS reformulation (Big-M method)." - ) + if not has_custom_upper: + upper = var.upper + if np.isinf(upper).any(): + raise ValueError( + f"Variable '{var.name}' has infinite upper bounds. " + "Finite bounds are required for SOS reformulation (Big-M method). " + "Alternatively, specify big_m in add_sos_constraints()." + ) def compute_big_m_values(var: Variable) -> tuple[DataArray, DataArray]: """ - Compute Big-M values from variable bounds. + Compute Big-M values from variable bounds or custom big_m attributes. + + If custom big_m values were specified via add_sos_constraints(big_m=...), + those are used (broadcast to variable dimensions). Otherwise, the variable + bounds are used as the tightest valid Big-M values. Parameters ---------- var : Variable - Variable with finite bounds. + Variable with finite bounds (or custom big_m attributes). Returns ------- tuple[DataArray, DataArray] - (M_upper, M_lower) - Big-M values computed from bounds. - M_upper = upper bound (for x <= U * y constraints) - M_lower = lower bound (for x >= L * y constraints) + (M_upper, M_lower) - Big-M values for reformulation constraints. + M_upper is used for: x <= M_upper * y + M_lower is used for: x >= M_lower * y """ - return var.upper, var.lower + import xarray as xr + + # Check for custom big_m values in attributes (always scalars) + big_m_upper = var.attrs.get("big_m_upper") + big_m_lower = var.attrs.get("big_m_lower") + + if big_m_upper is not None: + # Broadcast scalar to variable dimensions using xr.full_like + M_upper = xr.full_like(var.labels, big_m_upper, dtype=float) + else: + M_upper = var.upper + + if big_m_lower is not None: + M_lower = xr.full_like(var.labels, big_m_lower, dtype=float) + else: + M_lower = var.lower + + return M_upper, M_lower def reformulate_sos1(model: Model, var: Variable, prefix: str) -> None: diff --git a/test/test_sos_reformulation.py b/test/test_sos_reformulation.py index 0103e44b..39f022e5 100644 --- a/test/test_sos_reformulation.py +++ b/test/test_sos_reformulation.py @@ -101,6 +101,37 @@ def test_varying_bounds(self) -> None: assert np.allclose(M_upper.values, [1, 2, 3]) assert np.allclose(M_lower.values, [-1, -2, -3]) + def test_custom_big_m_scalar(self) -> None: + """Test Big-M with custom scalar value.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=100, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=10) + M_upper, M_lower = compute_big_m_values(x) + assert np.allclose(M_upper.values, [10, 10, 10]) + assert np.allclose(M_lower.values, [-10, -10, -10]) + + def test_custom_big_m_tuple(self) -> None: + """Test Big-M with custom tuple (upper, lower).""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=-100, upper=100, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=(5, -3)) + M_upper, M_lower = compute_big_m_values(x) + assert np.allclose(M_upper.values, [5, 5, 5]) + assert np.allclose(M_lower.values, [-3, -3, -3]) + + def test_custom_big_m_allows_infinite_bounds(self) -> None: + """Test that custom big_m allows variables with infinite bounds.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + 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) + # Should not raise - custom big_m bypasses bound validation + validate_bounds_for_reformulation(x) + M_upper, M_lower = compute_big_m_values(x) + assert np.allclose(M_upper.values, [10, 10, 10]) + class TestSOS1Reformulation: """Tests for SOS1 reformulation.""" @@ -616,3 +647,76 @@ def test_float_coordinates(self) -> None: z = m.variables["_sos_reform_lambda_z"] # Segment indicators have n-1 = 2 elements assert z.sizes["bp"] == 2 + + def test_custom_big_m_removed_on_remove_sos(self) -> None: + """Test that custom big_m attributes are removed with SOS constraint.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=100, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=10) + + assert "big_m_upper" in x.attrs + assert "big_m_lower" in x.attrs + + m.remove_sos_constraints(x) + + assert "big_m_upper" not in x.attrs + assert "big_m_lower" not in x.attrs + + +@pytest.mark.skipif("highs" not in available_solvers, reason="HiGHS not installed") +class TestCustomBigM: + """Tests for custom Big-M functionality.""" + + def test_solve_with_custom_big_m(self) -> None: + """Test solving with custom big_m value.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + # Large bounds but tight effective constraint + x = m.add_variables(lower=0, upper=1000, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=1) + m.add_objective(x * np.array([1, 2, 3]), sense="max") + + m.solve(solver_name="highs", reformulate_sos=True) + + # With big_m=1, maximum should be 3 (x[2]=1) + assert m.objective.value is not None + assert np.isclose(m.objective.value, 3, atol=1e-5) + + def test_solve_with_infinite_bounds_and_custom_big_m(self) -> None: + """Test solving with infinite bounds but custom big_m.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + 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=5) + m.add_objective(x * np.array([1, 2, 3]), sense="max") + + # Should work because big_m is specified + m.solve(solver_name="highs", reformulate_sos=True) + + assert m.objective.value is not None + # Max is 15 (x[2]=5) + assert np.isclose(m.objective.value, 15, atol=1e-5) + + def test_custom_big_m_tuple(self) -> None: + """Test solving with asymmetric big_m tuple.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=-100, upper=100, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=(2, -1)) + m.add_objective(x * np.array([1, 2, 3]), sense="max") + + m.solve(solver_name="highs", reformulate_sos=True) + + # Max is 6 (x[2]=2, limited by big_m_upper=2) + assert m.objective.value is not None + assert np.isclose(m.objective.value, 6, atol=1e-5) + + def test_big_m_invalid_tuple_length(self) -> None: + """Test that invalid tuple length raises error.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") + + with pytest.raises(ValueError, match="exactly 2 elements"): + m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=(1, 2, 3)) From 2258692bf921db8de26a2fc856c863b7d6a2d952 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 09:54:08 +0100 Subject: [PATCH 05/14] Add custom big_m parameter for SOS reformulation Allow users to specify custom Big-M values in add_sos_constraints() for tighter LP relaxations when variable bounds are conservative. - Add big_m parameter: scalar or tuple(upper, lower) - Store as variable attrs (big_m_upper, big_m_lower) for NetCDF persistence - Use tighter of big_m and variable bounds: min() for upper, max() for lower - Skip bound validation when custom big_m provided (allows infinite bounds) Scalar-only design ensures NetCDF persistence works correctly. For per-element Big-M values, users should adjust variable bounds directly. --- doc/sos-constraints.rst | 9 +++++++++ linopy/model.py | 9 +++------ linopy/sos_reformulation.py | 24 +++++++++++++----------- test/test_sos_reformulation.py | 7 +++++-- 4 files changed, 30 insertions(+), 19 deletions(-) diff --git a/doc/sos-constraints.rst b/doc/sos-constraints.rst index 129ded3f..8ac70f42 100644 --- a/doc/sos-constraints.rst +++ b/doc/sos-constraints.rst @@ -335,6 +335,15 @@ Big-M values than the variable bounds. This is useful when: # Tuple of scalars: asymmetric (M_upper = 10, M_lower = -5) m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=(10, -5)) +.. note:: + + The reformulation uses the **tighter** of custom ``big_m`` and variable bounds: + + - ``M_upper = min(big_m_upper, var.upper)`` + - ``M_lower = max(big_m_lower, var.lower)`` + + This ensures a loose ``big_m`` won't make the relaxation worse than using bounds alone. + **Why tighter Big-M matters:** Tighter Big-M values lead to: diff --git a/linopy/model.py b/linopy/model.py index d8bf80f3..c7d0b876 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -575,17 +575,14 @@ def add_sos_constraints( Big-M value(s) for SOS reformulation. Only used when reformulating SOS constraints for solvers that don't support them natively. - - None (default): Use variable bounds as Big-M (tightest valid choice) + - None (default): Use variable bounds as Big-M - float: Symmetric Big-M (M_upper = big_m, M_lower = -big_m) - tuple (upper, lower): Asymmetric Big-M values - The Big-M values are used in constraints: - - x <= M_upper * y (upper linking) - - x >= M_lower * y (lower linking, M_lower should be negative) + The reformulation uses the tighter of big_m and variable bounds: + M_upper = min(big_m_upper, var.upper), M_lower = max(big_m_lower, var.lower). Tighter Big-M values improve LP relaxation quality and solve time. - If you know tighter effective bounds from other constraints in your - model, specifying them here can significantly improve performance. """ if sos_type not in (1, 2): raise ValueError(f"sos_type must be 1 or 2, got {sos_type}") diff --git a/linopy/sos_reformulation.py b/linopy/sos_reformulation.py index 06bfc7eb..fb7a1f39 100644 --- a/linopy/sos_reformulation.py +++ b/linopy/sos_reformulation.py @@ -63,11 +63,11 @@ def validate_bounds_for_reformulation(var: Variable) -> None: def compute_big_m_values(var: Variable) -> tuple[DataArray, DataArray]: """ - Compute Big-M values from variable bounds or custom big_m attributes. + Compute Big-M values from variable bounds and custom big_m attributes. - If custom big_m values were specified via add_sos_constraints(big_m=...), - those are used (broadcast to variable dimensions). Otherwise, the variable - bounds are used as the tightest valid Big-M values. + Uses the tighter of variable bounds and custom big_m values. This ensures + the best possible LP relaxation - a loose big_m won't make things worse + if variable bounds are already tight. Parameters ---------- @@ -87,16 +87,18 @@ def compute_big_m_values(var: Variable) -> tuple[DataArray, DataArray]: big_m_upper = var.attrs.get("big_m_upper") big_m_lower = var.attrs.get("big_m_lower") + # Start with variable bounds + M_upper = var.upper + M_lower = var.lower + + # Use tighter of custom big_m and variable bounds if big_m_upper is not None: - # Broadcast scalar to variable dimensions using xr.full_like - M_upper = xr.full_like(var.labels, big_m_upper, dtype=float) - else: - M_upper = var.upper + custom_upper = xr.full_like(var.labels, big_m_upper, dtype=float) + M_upper = xr.where(custom_upper < M_upper, custom_upper, M_upper) if big_m_lower is not None: - M_lower = xr.full_like(var.labels, big_m_lower, dtype=float) - else: - M_lower = var.lower + custom_lower = xr.full_like(var.labels, big_m_lower, dtype=float) + M_lower = xr.where(custom_lower > M_lower, custom_lower, M_lower) return M_upper, M_lower diff --git a/test/test_sos_reformulation.py b/test/test_sos_reformulation.py index 39f022e5..17455753 100644 --- a/test/test_sos_reformulation.py +++ b/test/test_sos_reformulation.py @@ -102,14 +102,17 @@ def test_varying_bounds(self) -> None: assert np.allclose(M_lower.values, [-1, -2, -3]) def test_custom_big_m_scalar(self) -> None: - """Test Big-M with custom scalar value.""" + """Test Big-M uses tighter of custom value and bounds.""" m = Model() idx = pd.Index([0, 1, 2], name="i") x = m.add_variables(lower=0, upper=100, coords=[idx], name="x") m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=10) M_upper, M_lower = compute_big_m_values(x) + # big_m=10 gives big_m_upper=10, big_m_lower=-10 + # M_upper = min(10, 100) = 10 (custom is tighter) + # M_lower = max(-10, 0) = 0 (bound is tighter) assert np.allclose(M_upper.values, [10, 10, 10]) - assert np.allclose(M_lower.values, [-10, -10, -10]) + assert np.allclose(M_lower.values, [0, 0, 0]) def test_custom_big_m_tuple(self) -> None: """Test Big-M with custom tuple (upper, lower).""" From 0a3dde514179f385073083cc7291768abe8e5d56 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 10:19:28 +0100 Subject: [PATCH 06/14] =?UTF-8?q?Simplification=20summary:=20=20=20?= =?UTF-8?q?=E2=94=8C=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=AC=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=AC=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=AC=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=90=20=20=20=E2=94=82=20=20=20?= =?UTF-8?q?=20=20=20=20=20=20File=20=20=20=20=20=20=20=20=20=E2=94=82=20?= =?UTF-8?q?=20Before=20=20=20=E2=94=82=20=20=20After=20=20=20=E2=94=82=20R?= =?UTF-8?q?eduction=20=E2=94=82=20=20=20=E2=94=9C=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=BC?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=BC=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=BC=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=A4=20=20=20=E2=94=82=20sos=5Freformulation.py=20?= =?UTF-8?q?=E2=94=82=20377=20lines=20=E2=94=82=20223=20lines=20=E2=94=82?= =?UTF-8?q?=2041%=20=20=20=20=20=20=20=E2=94=82=20=20=20=E2=94=9C=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=BC=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=BC=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=BC=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=A4=20=20=20=E2=94=82=20sos-constraints.rst=20?= =?UTF-8?q?=20=E2=94=82=20647=20lines=20=E2=94=82=20164=20lines=20?= =?UTF-8?q?=E2=94=82=2075%=20=20=20=20=20=20=20=E2=94=82=20=20=20=E2=94=94?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=B4=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=B4?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=B4=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80=E2=94=80?= =?UTF-8?q?=E2=94=80=E2=94=80=E2=94=98=20=20=20Code=20changes:=20=20=20-?= =?UTF-8?q?=20Merged=20validate=5Fbounds=5Ffor=5Freformulation=20into=20co?= =?UTF-8?q?mpute=5Fbig=5Fm=5Fvalues=20=20=20-=20Factored=20out=20add=5Flin?= =?UTF-8?q?king=5Fconstraints=20helper=20in=20SOS2=20=20=20-=20Used=20np.m?= =?UTF-8?q?inimum/np.maximum=20instead=20of=20xr.where=20=20=20-=20Kept=20?= =?UTF-8?q?proper=20docstrings=20with=20Parameters/Returns=20sections?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Doc changes: - Removed: Variable Representation, LP File Export, Common Patterns, Performance Considerations - Trimmed: Examples to one each, Mathematical formulation to equations only - Condensed: API reference, multi-dimensional explanation --- doc/sos-constraints.rst | 574 +++------------------------------ linopy/sos_reformulation.py | 354 ++++++-------------- test/test_sos_reformulation.py | 15 +- 3 files changed, 156 insertions(+), 787 deletions(-) diff --git a/doc/sos-constraints.rst b/doc/sos-constraints.rst index 8ac70f42..5a17c82b 100644 --- a/doc/sos-constraints.rst +++ b/doc/sos-constraints.rst @@ -3,7 +3,7 @@ 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. +SOS constraints model situations where only one (SOS1) or two adjacent (SOS2) variables from an ordered set can be non-zero. .. contents:: :local: @@ -12,500 +12,104 @@ Special Ordered Sets (SOS) are a constraint type used in mixed-integer programmi 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. The adjacency is determined by the ordering weights (coordinates) of the variables. - -**Example use cases:** -- Piecewise linear approximation of nonlinear functions -- Portfolio optimization with discrete risk levels -- Production planning with discrete capacity levels +- **SOS1**: At most one variable can be non-zero (mutually exclusive choices) +- **SOS2**: At most two adjacent variables can be non-zero (piecewise linear functions) Basic Usage ----------- -Adding SOS Constraints -~~~~~~~~~~~~~~~~~~~~~~~ - -To add SOS constraints to variables in linopy: - .. code-block:: python 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 + # SOS1: at most one option selected + options = pd.Index([0, 1, 2], name="options") + x = m.add_variables(coords=[options], name="x", lower=0, upper=1) 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") + # SOS2: at most two adjacent breakpoints active + breakpoints = pd.Index([0.0, 1.0, 2.0], name="bp") 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 -~~~~~~~~~~~~~~~~ - -.. code-block:: python - - Model.add_sos_constraints(variable, sos_type, sos_dim) - -**Parameters:** - -- ``variable`` : Variable - The variable to which the SOS constraint should be applied -- ``sos_type`` : {1, 2} - Type of SOS constraint (1 or 2) -- ``sos_dim`` : str - Name of the dimension along which the SOS constraint applies + m.add_sos_constraints(lambdas, sos_type=2, sos_dim="bp") **Requirements:** -- The specified dimension must exist in the variable -- The coordinates for the SOS dimension must be numeric (used as weights for ordering) -- Only one SOS constraint can be applied per variable - -Examples --------- - -Example 1: Facility Location (SOS1) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -.. code-block:: python - - 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) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -.. code-block:: python - - 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") +- The SOS dimension must have numeric coordinates (used as ordering weights) +- Only one SOS constraint per variable -Working with Multi-dimensional Variables ------------------------------------------ +Multi-dimensional Variables +~~~~~~~~~~~~~~~~~~~~~~~~~~~ -SOS constraints are created for each dimension that is not sos_dim. +For multi-dimensional variables, the SOS constraint applies independently for each combination of non-SOS dimensions: .. code-block:: python - # Multi-period production planning - periods = pd.Index(range(3), name="periods") + periods = pd.Index([0, 1], 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: - -.. code-block:: python - - # 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: - -.. code-block:: python - - 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: - -.. code-block:: text - - sos - - s0: S1 :: x0:0 x1:1 x2:2 - s3: S2 :: x3:0.0 x4:1.0 x5:2.0 + x = m.add_variables(lower=0, upper=1, coords=[periods, modes], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="modes") + # Result: at most one mode selected PER period 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:** +**Native SOS support:** Gurobi, CPLEX, CBC, SCIP, Xpress -- Gurobi (via ``gurobipy``) - -**Unsupported solvers:** - -- HiGHS (does not support SOS constraints) -- GLPK -- MOSEK -- MindOpt - -For these solvers, linopy provides automatic reformulation (see :ref:`sos-reformulation` below). +**No SOS support (use reformulation):** HiGHS, GLPK, MOSEK .. _sos-reformulation: -SOS Reformulation for Unsupported Solvers ------------------------------------------ +SOS Reformulation +----------------- -Linopy can automatically reformulate SOS constraints as binary + linear constraints -using the Big-M method. This allows you to use SOS constraints with solvers that -don't support them natively (HiGHS, GLPK, MOSEK, etc.). - -Enabling Reformulation -~~~~~~~~~~~~~~~~~~~~~~ - -Pass ``reformulate_sos=True`` to the ``solve()`` method: +For solvers without native SOS support, linopy can reformulate SOS constraints as binary + linear constraints using the Big-M method. .. code-block:: python - import linopy - import pandas as pd - - m = linopy.Model() - idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") - m.add_sos_constraints(x, sos_type=1, sos_dim="i") - m.add_objective(x.sum(), sense="max") - - # Now works with HiGHS! + # Automatic reformulation during solve m.solve(solver_name="highs", reformulate_sos=True) -You can also reformulate manually before solving: - -.. code-block:: python - - # Reformulate in place - reformulated_vars = m.reformulate_sos_constraints() - print(f"Reformulated: {reformulated_vars}") - - # Then solve with any solver + # Or reformulate manually + m.reformulate_sos_constraints() m.solve(solver_name="highs") -Requirements +Big-M Values ~~~~~~~~~~~~ -**Finite bounds or custom Big-M required.** The reformulation uses the Big-M method. -By default, Big-M values are derived from variable bounds. If any SOS variable has -infinite bounds, you must either: - -1. Set finite bounds on the variable, or -2. Specify custom Big-M values via the ``big_m`` parameter +Big-M values are derived from variable bounds by default. For infinite bounds, specify custom values: .. code-block:: python - # Option 1: Finite bounds (default Big-M = bounds) + # Finite bounds: Big-M = 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") - # Option 2: Custom Big-M (allows infinite bounds) + # Infinite bounds: specify Big-M explicitly 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) -Custom Big-M Values -~~~~~~~~~~~~~~~~~~~ - -The ``big_m`` parameter in ``add_sos_constraints()`` allows you to specify tighter -Big-M values than the variable bounds. This is useful when: - -- Variable bounds are conservatively large (e.g., ``upper=1e6``) -- Other constraints in your model imply tighter effective bounds -- You have domain knowledge about the actual feasible range - -**Syntax options:** - -.. code-block:: python - - # Scalar: symmetric Big-M (M_upper = 10, M_lower = -10) - m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=10) - - # Tuple of scalars: asymmetric (M_upper = 10, M_lower = -5) + # Asymmetric Big-M m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=(10, -5)) -.. note:: - - The reformulation uses the **tighter** of custom ``big_m`` and variable bounds: - - - ``M_upper = min(big_m_upper, var.upper)`` - - ``M_lower = max(big_m_lower, var.lower)`` - - This ensures a loose ``big_m`` won't make the relaxation worse than using bounds alone. - -**Why tighter Big-M matters:** - -Tighter Big-M values lead to: - -- Better LP relaxation bounds (closer to optimal integer solution) -- Fewer branch-and-bound nodes -- Faster solve times - -.. code-block:: python - - # Example: Variable has large bounds but we know effective range is smaller - x = m.add_variables(lower=0, upper=1000, coords=[idx], name="x") - - # Other constraints limit x to [0, 10] in practice - m.add_constraints(x <= 10) - - # Use tight Big-M for better performance - m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=10) +The reformulation uses the tighter of ``big_m`` and variable bounds. Mathematical Formulation ~~~~~~~~~~~~~~~~~~~~~~~~ -The reformulation converts SOS constraints into Mixed-Integer Linear Programming -(MILP) constraints using binary indicator variables. - -**SOS1 Reformulation** - -Given variables :math:`x_i` for :math:`i \in I` with bounds :math:`L_i \leq x_i \leq U_i`, -the SOS1 constraint "at most one :math:`x_i` is non-zero" is reformulated as: +**SOS1:** Binary indicators :math:`y_i`, linking constraints, cardinality :math:`\sum y_i \leq 1` .. math:: - \text{Binary indicators:} \quad & y_i \in \{0, 1\} \quad \forall i \in I \\[0.5em] - \text{Upper linking:} \quad & x_i \leq U_i \cdot y_i \quad \forall i \in I \text{ where } U_i > 0 \\[0.5em] - \text{Lower linking:} \quad & x_i \geq L_i \cdot y_i \quad \forall i \in I \text{ where } L_i < 0 \\[0.5em] - \text{Cardinality:} \quad & \sum_{i \in I} y_i \leq 1 - -**Interpretation:** - -- :math:`y_i = 1` means variable :math:`x_i` is "selected" (allowed to be non-zero) -- :math:`y_i = 0` forces :math:`x_i = 0` via the linking constraints -- The cardinality constraint ensures at most one :math:`y_i = 1` + y_i \in \{0, 1\}, \quad x_i \leq U_i \cdot y_i, \quad x_i \geq L_i \cdot y_i, \quad \sum y_i \leq 1 -**Example:** For :math:`x \in [0, 10]`: - -- If :math:`y = 0`: constraint :math:`x \leq 10 \cdot 0 = 0` forces :math:`x = 0` -- If :math:`y = 1`: constraint :math:`x \leq 10 \cdot 1 = 10` allows :math:`x \in [0, 10]` - -**SOS2 Reformulation** - -Given ordered variables :math:`x_0, x_1, \ldots, x_{n-1}` with bounds :math:`L_i \leq x_i \leq U_i`, -the SOS2 constraint "at most two adjacent :math:`x_i` are non-zero" is reformulated using -segment indicators: +**SOS2:** Segment indicators :math:`z_j` for :math:`j = 0, \ldots, n-2`, adjacency constraints .. math:: - \text{Segment indicators:} \quad & z_j \in \{0, 1\} \quad \forall j \in \{0, \ldots, n-2\} \\[0.5em] - \text{First variable:} \quad & x_0 \leq U_0 \cdot z_0 \\[0.5em] - \text{Middle variables:} \quad & x_i \leq U_i \cdot (z_{i-1} + z_i) \quad \forall i \in \{1, \ldots, n-2\} \\[0.5em] - \text{Last variable:} \quad & x_{n-1} \leq U_{n-1} \cdot z_{n-2} \\[0.5em] - \text{Cardinality:} \quad & \sum_{j=0}^{n-2} z_j \leq 1 - -(Similar constraints for lower bounds when :math:`L_i < 0`) - -**Interpretation:** - -- :math:`z_j = 1` means "segment :math:`j`" is active (between positions :math:`j` and :math:`j+1`) -- Variable :math:`x_i` can only be non-zero if an adjacent segment is active -- The cardinality constraint ensures at most one segment is active -- This guarantees at most two adjacent variables :math:`(x_j, x_{j+1})` are non-zero - -**Example:** For :math:`n = 4` variables :math:`(x_0, x_1, x_2, x_3)`: - -- If :math:`z_1 = 1` (segment 1 active): :math:`x_1` and :math:`x_2` can be non-zero -- Constraints force :math:`x_0 = 0` (needs :math:`z_0`) and :math:`x_3 = 0` (needs :math:`z_2`) - -Auxiliary Variables and Constraints -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -The reformulation creates auxiliary variables and constraints with a ``_sos_reform_`` prefix: - -.. code-block:: python - - m.reformulate_sos_constraints() - - # For SOS1 variable 'x': - # - Binary indicators: _sos_reform_x_y - # - Upper constraints: _sos_reform_x_upper - # - Lower constraints: _sos_reform_x_lower (if L < 0) - # - Cardinality: _sos_reform_x_card - - # For SOS2 variable 'x': - # - Segment indicators: _sos_reform_x_z - # - Upper constraints: _sos_reform_x_upper_first, _sos_reform_x_upper_mid_*, _sos_reform_x_upper_last - # - Lower constraints: _sos_reform_x_lower_first, _sos_reform_x_lower_mid_*, _sos_reform_x_lower_last - # - Cardinality: _sos_reform_x_card - -You can use a custom prefix: - -.. code-block:: python - - m.reformulate_sos_constraints(prefix="_my_sos_") - -Multi-dimensional Variables -~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -For multi-dimensional variables, the reformulation respects xarray broadcasting. -Constraints are created for each combination of non-SOS dimensions: - -.. code-block:: python - - # 2D variable: periods × options - periods = pd.Index([0, 1], name="periods") - options = pd.Index([0, 1, 2], name="options") - x = m.add_variables(lower=0, upper=1, coords=[periods, options], name="x") - m.add_sos_constraints(x, sos_type=1, sos_dim="options") - - # After reformulation: - # - Binary y has shape (periods: 2, options: 3) - # - Cardinality constraint: sum over 'options' for each period - # - Result: at most one option selected PER period - -Edge Cases -~~~~~~~~~~ - -The reformulation handles several edge cases automatically: - -.. list-table:: - :header-rows: 1 - :widths: 30 70 - - * - Case - - Handling - * - Single-element SOS - - Skipped (trivially satisfied) - * - All-zero bounds (L=U=0) - - Skipped (variable already fixed to 0) - * - All-positive bounds - - Only upper linking constraints created - * - All-negative bounds - - Only lower linking constraints created - * - Mixed bounds - - Both upper and lower constraints created - -Performance Considerations -~~~~~~~~~~~~~~~~~~~~~~~~~~ - -**Reformulation trade-offs:** - -- **Pros:** Works with any MIP solver; explicit constraints can sometimes help the solver -- **Cons:** Adds binary variables and constraints; may be slower than native SOS support - -**When to use native SOS:** - -- Use native SOS (Gurobi, CPLEX) when available—solvers have specialized branching strategies -- Native SOS often provides better performance for piecewise linear problems - -**When reformulation is useful:** - -- When you must use a solver without SOS support (HiGHS, GLPK) -- For model portability across different solvers -- When debugging—explicit constraints are easier to inspect + z_j \in \{0, 1\}, \quad x_0 \leq U_0 \cdot z_0, \quad x_i \leq U_i(z_{i-1} + z_i), \quad x_{n-1} \leq U_{n-1} \cdot z_{n-2}, \quad \sum z_j \leq 1 Example: Piecewise Linear with HiGHS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -514,133 +118,47 @@ Example: Piecewise Linear with HiGHS import linopy import pandas as pd - import numpy as np - # Approximate f(x) = x² using piecewise linear with HiGHS m = linopy.Model() - breakpoints = pd.Index([0.0, 1.0, 2.0, 3.0], name="bp") - x_vals = breakpoints.to_numpy() - y_vals = x_vals**2 # [0, 1, 4, 9] + # Approximate f(x) = x² over [0, 3] + bp = pd.Index([0.0, 1.0, 2.0, 3.0], name="bp") + x_vals, y_vals = bp.to_numpy(), bp.to_numpy() ** 2 - # SOS2 interpolation weights - lambdas = m.add_variables(lower=0, upper=1, coords=[breakpoints], name="lambdas") + lambdas = m.add_variables(lower=0, upper=1, coords=[bp], name="lambdas") m.add_sos_constraints(lambdas, sos_type=2, sos_dim="bp") - # Interpolated point 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).sum(), name="x_interp") m.add_constraints(y == (lambdas * y_vals).sum(), name="y_interp") m.add_constraints(x >= 1.5, name="x_min") - - # Minimize y (approximated x²) m.add_objective(y) - # Solve with HiGHS using reformulation m.solve(solver_name="highs", reformulate_sos=True) - print(f"x = {x.solution.item():.2f}") - print(f"y ≈ x² = {y.solution.item():.2f}") - # Output: x = 1.50, y ≈ x² = 2.50 (exact: 2.25) - -Common Patterns ---------------- - -Piecewise Linear Cost Function -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -.. code-block:: python - - 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 -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -.. code-block:: python - - 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 --------- - -- :doc:`creating-variables`: Creating variables with coordinates -- :doc:`creating-constraints`: Adding regular constraints -- :doc:`user-guide`: General linopy usage patterns - API Reference ------------- .. py:method:: Model.add_sos_constraints(variable, sos_type, sos_dim, big_m=None) - Add an SOS1 or SOS2 constraint for one dimension of a variable. + Add SOS constraint to a variable. :param variable: Variable to constrain - :type variable: Variable - :param sos_type: Type of SOS constraint (1 or 2) - :type sos_type: Literal[1, 2] - :param sos_dim: Dimension along which to apply the SOS constraint - :type sos_dim: str - :param big_m: Custom Big-M value(s) for reformulation. Can be: - - - ``None`` (default): Use variable bounds - - ``float``: Symmetric Big-M (upper = big_m, lower = -big_m) - - ``tuple[float, float]``: Asymmetric (upper, lower) - - :type big_m: float | tuple[float, float] | None + :param sos_type: 1 (at most one non-zero) or 2 (at most two adjacent non-zero) + :param sos_dim: Dimension for SOS ordering + :param big_m: Custom Big-M: ``float`` (symmetric) or ``tuple[float, float]`` (upper, lower) .. py:method:: Model.remove_sos_constraints(variable) Remove SOS constraints from a variable. - :param variable: Variable from which to remove SOS constraints - :type variable: Variable - .. py:method:: Model.reformulate_sos_constraints(prefix="_sos_reform_") - Reformulate SOS constraints as binary + linear constraints using the Big-M method. - - :param prefix: Prefix for auxiliary variable and constraint names - :type prefix: str - :returns: List of variable names that were reformulated - :rtype: list[str] - :raises ValueError: If any SOS variable has infinite bounds and no custom big_m + Convert SOS to binary + linear constraints. Returns list of reformulated variable names. .. py:attribute:: Variables.sos - Property returning all variables with SOS constraints. + All variables with SOS constraints. diff --git a/linopy/sos_reformulation.py b/linopy/sos_reformulation.py index fb7a1f39..37f0d963 100644 --- a/linopy/sos_reformulation.py +++ b/linopy/sos_reformulation.py @@ -1,8 +1,8 @@ """ -Linopy SOS constraint reformulation module. +SOS constraint reformulation using Big-M method. -This module provides functions to reformulate SOS1 and SOS2 constraints -as binary + linear constraints for solvers that don't support them natively. +Converts SOS1/SOS2 constraints to binary + linear constraints for solvers +that don't support them natively. """ from __future__ import annotations @@ -21,297 +21,168 @@ logger = logging.getLogger(__name__) -def validate_bounds_for_reformulation(var: Variable) -> None: - """ - Validate that a variable has finite bounds required for SOS reformulation. - - If custom big_m values are provided via variable attributes, bounds - validation is skipped for those dimensions. - - Parameters - ---------- - var : Variable - Variable to validate. - - Raises - ------ - ValueError - If any bound is infinite and no custom big_m is provided. - """ - # If custom big_m is provided, skip bound validation - has_custom_upper = var.attrs.get("big_m_upper") is not None - has_custom_lower = var.attrs.get("big_m_lower") is not None - - if not has_custom_lower: - lower = var.lower - if np.isinf(lower).any(): - raise ValueError( - f"Variable '{var.name}' has infinite lower bounds. " - "Finite bounds are required for SOS reformulation (Big-M method). " - "Alternatively, specify big_m in add_sos_constraints()." - ) - - if not has_custom_upper: - upper = var.upper - if np.isinf(upper).any(): - raise ValueError( - f"Variable '{var.name}' has infinite upper bounds. " - "Finite bounds are required for SOS reformulation (Big-M method). " - "Alternatively, specify big_m in add_sos_constraints()." - ) - - def compute_big_m_values(var: Variable) -> tuple[DataArray, DataArray]: """ Compute Big-M values from variable bounds and custom big_m attributes. - Uses the tighter of variable bounds and custom big_m values. This ensures - the best possible LP relaxation - a loose big_m won't make things worse - if variable bounds are already tight. + Uses the tighter of variable bounds and custom big_m values to ensure + the best possible LP relaxation. Parameters ---------- var : Variable - Variable with finite bounds (or custom big_m attributes). + Variable with bounds (and optionally big_m_upper/big_m_lower attrs). Returns ------- tuple[DataArray, DataArray] - (M_upper, M_lower) - Big-M values for reformulation constraints. - M_upper is used for: x <= M_upper * y - M_lower is used for: x >= M_lower * y - """ - import xarray as xr + (M_upper, M_lower) for reformulation constraints: + x <= M_upper * y and x >= M_lower * y - # Check for custom big_m values in attributes (always scalars) + Raises + ------ + ValueError + If resulting Big-M values are infinite. + """ big_m_upper = var.attrs.get("big_m_upper") big_m_lower = var.attrs.get("big_m_lower") - # Start with variable bounds M_upper = var.upper M_lower = var.lower - # Use tighter of custom big_m and variable bounds if big_m_upper is not None: - custom_upper = xr.full_like(var.labels, big_m_upper, dtype=float) - M_upper = xr.where(custom_upper < M_upper, custom_upper, M_upper) - + M_upper = np.minimum(M_upper, big_m_upper) if big_m_lower is not None: - custom_lower = xr.full_like(var.labels, big_m_lower, dtype=float) - M_lower = xr.where(custom_lower > M_lower, custom_lower, M_lower) + M_lower = np.maximum(M_lower, big_m_lower) + + # Validate finiteness + if np.isinf(M_upper).any(): + raise ValueError( + f"Variable '{var.name}' has infinite upper bounds. " + "Set finite bounds or specify big_m in add_sos_constraints()." + ) + if np.isinf(M_lower).any(): + raise ValueError( + f"Variable '{var.name}' has infinite lower bounds. " + "Set finite bounds or specify big_m in add_sos_constraints()." + ) return M_upper, M_lower def reformulate_sos1(model: Model, var: Variable, prefix: str) -> None: """ - Reformulate an SOS1 constraint as binary + linear constraints. + Reformulate SOS1 constraint as binary + linear constraints. - SOS1: At most one variable can be non-zero. - - Reformulation: - For each x[i] with bounds [L[i], U[i]]: - - Add binary y[i] - - x[i] <= U[i] * y[i] (if U[i] > 0) - - x[i] >= L[i] * y[i] (if L[i] < 0) - - sum(y) <= 1 + For each x[i] with bounds [L[i], U[i]]: + - Add binary indicator y[i] + - x[i] <= U[i] * y[i] (upper linking, if U > 0) + - x[i] >= L[i] * y[i] (lower linking, if L < 0) + - sum(y) <= 1 (cardinality) Parameters ---------- model : Model - The model to add reformulation constraints to. + Model to add reformulation constraints to. var : Variable - The variable with SOS1 constraint. + Variable with SOS1 constraint. prefix : str Prefix for naming auxiliary variables and constraints. """ sos_dim = var.attrs["sos_dim"] - var_name = var.name - - # Get bounds + name = var.name M_upper, M_lower = compute_big_m_values(var) - # Extract coords as list of DataArrays for add_variables - coords_list = [var.coords[d] for d in var.dims] - - # Create binary indicator variables with same dimensions as original variable - y_name = f"{prefix}{var_name}_y" - y = model.add_variables( - coords=coords_list, - name=y_name, - binary=True, - ) - - # Add upper bound constraints: x <= U * y (when U > 0) - # This ensures x can only be positive when y = 1 - has_positive_upper = (M_upper > 0).any() - if has_positive_upper: - model.add_constraints( - var <= M_upper * y, - name=f"{prefix}{var_name}_upper", - ) - - # Add lower bound constraints: x >= L * y (when L < 0) - # This ensures x can only be negative when y = 1 - has_negative_lower = (M_lower < 0).any() - if has_negative_lower: - model.add_constraints( - var >= M_lower * y, - name=f"{prefix}{var_name}_lower", - ) + coords = [var.coords[d] for d in var.dims] + y = model.add_variables(coords=coords, name=f"{prefix}{name}_y", binary=True) - # Add cardinality constraint: sum(y) <= 1 over the SOS dimension - # This is summed over sos_dim, keeping other dimensions - model.add_constraints( - y.sum(dim=sos_dim) <= 1, - name=f"{prefix}{var_name}_card", - ) + if (M_upper > 0).any(): + model.add_constraints(var <= M_upper * y, name=f"{prefix}{name}_upper") + if (M_lower < 0).any(): + model.add_constraints(var >= M_lower * y, name=f"{prefix}{name}_lower") - logger.debug(f"Reformulated SOS1 constraint for variable '{var_name}'") + model.add_constraints(y.sum(dim=sos_dim) <= 1, name=f"{prefix}{name}_card") def reformulate_sos2(model: Model, var: Variable, prefix: str) -> None: """ - Reformulate an SOS2 constraint as binary + linear constraints. - - SOS2: At most two adjacent variables can be non-zero. + Reformulate SOS2 constraint as binary + linear constraints. - Reformulation: - For ordered x[i], i = 0..n-1: - - Add binary z[i] for i = 0..n-2 (segment indicators) - - x[0] <= U[0] * z[0] - - x[i] <= U[i] * (z[i-1] + z[i]) for i = 1..n-2 - - x[n-1] <= U[n-1] * z[n-2] - - Similar for lower bounds if L[i] < 0 - - sum(z) <= 1 + For ordered x[0..n-1]: + - Add n-1 binary segment indicators z[i] + - x[0] <= U[0] * z[0] + - x[i] <= U[i] * (z[i-1] + z[i]) for middle elements + - x[n-1] <= U[n-1] * z[n-2] + - Similar for lower bounds if L < 0 + - sum(z) <= 1 Parameters ---------- model : Model - The model to add reformulation constraints to. + Model to add reformulation constraints to. var : Variable - The variable with SOS2 constraint. + Variable with SOS2 constraint. prefix : str Prefix for naming auxiliary variables and constraints. """ sos_dim = var.attrs["sos_dim"] - var_name = var.name - - # Get the size of the SOS dimension + name = var.name n = var.sizes[sos_dim] - # Trivial case: single element always satisfies SOS2 if n <= 1: - logger.debug( - f"Skipping SOS2 reformulation for '{var_name}' (trivial case: n={n})" - ) return - # Get bounds M_upper, M_lower = compute_big_m_values(var) - # Create n-1 binary segment indicators - # z[i] indicates that the "active segment" is between positions i and i+1 - segment_coords_values = var.coords[sos_dim].values[:-1] # n-1 segment indicators - - # Build coords for z: same as var but with truncated sos_dim - z_coords_list = [] - for d in var.dims: - if d == sos_dim: - z_coords_list.append(pd.Index(segment_coords_values, name=sos_dim)) - else: - z_coords_list.append(var.coords[d]) - - z_name = f"{prefix}{var_name}_z" - z = model.add_variables( - coords=z_coords_list, - name=z_name, - binary=True, - ) - - # For SOS2, each x[i] can only be non-zero if an adjacent segment is active - # x[0] needs z[0] active - # x[i] (1 <= i <= n-2) needs z[i-1] or z[i] active - # x[n-1] needs z[n-2] active - - # Convert to LinearExpression to avoid SOS attribute validation issues during isel - # (Variable's isel preserves SOS attrs, which fail validation when dim is removed) - x_expr = 1 * var - z_expr = 1 * z - - # Process upper bound constraints (for positive bounds) - has_positive_upper = (M_upper > 0).any() - if has_positive_upper: - # First element: x[0] <= U[0] * z[0] - x_first = x_expr.isel({sos_dim: 0}) - z_first = z_expr.isel({sos_dim: 0}) - U_first = M_upper.isel({sos_dim: 0}) - model.add_constraints( - x_first <= U_first * z_first, - name=f"{prefix}{var_name}_upper_first", - ) + # Create n-1 segment indicators + z_coords = [ + pd.Index(var.coords[sos_dim].values[:-1], name=sos_dim) + if d == sos_dim + else var.coords[d] + for d in var.dims + ] + z = model.add_variables(coords=z_coords, name=f"{prefix}{name}_z", binary=True) - # Middle elements: x[i] <= U[i] * (z[i-1] + z[i]) for i = 1..n-2 - if n > 2: - for i in range(1, n - 1): - x_i = x_expr.isel({sos_dim: i}) - z_prev = z_expr.isel({sos_dim: i - 1}) - z_curr = z_expr.isel({sos_dim: i}) - U_i = M_upper.isel({sos_dim: i}) - model.add_constraints( - x_i <= U_i * (z_prev + z_curr), - name=f"{prefix}{var_name}_upper_mid_{i}", - ) - - # Last element: x[n-1] <= U[n-1] * z[n-2] - x_last = x_expr.isel({sos_dim: n - 1}) - z_last = z_expr.isel({sos_dim: n - 2}) - U_last = M_upper.isel({sos_dim: n - 1}) - model.add_constraints( - x_last <= U_last * z_last, - name=f"{prefix}{var_name}_upper_last", - ) + # Convert to expressions to avoid SOS attr validation on isel + x, z_expr = 1 * var, 1 * z + + def add_linking_constraints(M: DataArray, sign: str, suffix: str) -> None: + """Add x <= M*z or x >= M*z constraints for first/middle/last elements.""" + if sign == "upper" and not (M > 0).any(): + return + if sign == "lower" and not (M < 0).any(): + return - # Process lower bound constraints (for negative bounds) - has_negative_lower = (M_lower < 0).any() - if has_negative_lower: - # First element: x[0] >= L[0] * z[0] - x_first = x_expr.isel({sos_dim: 0}) - z_first = z_expr.isel({sos_dim: 0}) - L_first = M_lower.isel({sos_dim: 0}) + op = (lambda a, b: a <= b) if sign == "upper" else (lambda a, b: a >= b) + + # First: x[0] op M[0] * z[0] model.add_constraints( - x_first >= L_first * z_first, - name=f"{prefix}{var_name}_lower_first", + op(x.isel({sos_dim: 0}), M.isel({sos_dim: 0}) * z_expr.isel({sos_dim: 0})), + name=f"{prefix}{name}_{suffix}_first", ) - - # Middle elements: x[i] >= L[i] * (z[i-1] + z[i]) for i = 1..n-2 - if n > 2: - for i in range(1, n - 1): - x_i = x_expr.isel({sos_dim: i}) - z_prev = z_expr.isel({sos_dim: i - 1}) - z_curr = z_expr.isel({sos_dim: i}) - L_i = M_lower.isel({sos_dim: i}) - model.add_constraints( - x_i >= L_i * (z_prev + z_curr), - name=f"{prefix}{var_name}_lower_mid_{i}", - ) - - # Last element: x[n-1] >= L[n-1] * z[n-2] - x_last = x_expr.isel({sos_dim: n - 1}) - z_last = z_expr.isel({sos_dim: n - 2}) - L_last = M_lower.isel({sos_dim: n - 1}) + # Middle: x[i] op M[i] * (z[i-1] + z[i]) + for i in range(1, n - 1): + model.add_constraints( + op( + x.isel({sos_dim: i}), + M.isel({sos_dim: i}) + * (z_expr.isel({sos_dim: i - 1}) + z_expr.isel({sos_dim: i})), + ), + name=f"{prefix}{name}_{suffix}_mid_{i}", + ) + # Last: x[n-1] op M[n-1] * z[n-2] model.add_constraints( - x_last >= L_last * z_last, - name=f"{prefix}{var_name}_lower_last", + op( + x.isel({sos_dim: n - 1}), + M.isel({sos_dim: n - 1}) * z_expr.isel({sos_dim: n - 2}), + ), + name=f"{prefix}{name}_{suffix}_last", ) - # Add cardinality constraint: sum(z) <= 1 - model.add_constraints( - z.sum(dim=sos_dim) <= 1, - name=f"{prefix}{var_name}_card", - ) + add_linking_constraints(M_upper, "upper", "upper") + add_linking_constraints(M_lower, "lower", "lower") - logger.debug(f"Reformulated SOS2 constraint for variable '{var_name}'") + model.add_constraints(z.sum(dim=sos_dim) <= 1, name=f"{prefix}{name}_card") def reformulate_all_sos(model: Model, prefix: str = "_sos_reform_") -> list[str]: @@ -321,56 +192,39 @@ def reformulate_all_sos(model: Model, prefix: str = "_sos_reform_") -> list[str] Parameters ---------- model : Model - The model containing SOS constraints to reformulate. + Model containing SOS constraints to reformulate. prefix : str, optional - Prefix for naming auxiliary variables and constraints. - Default is "_sos_reform_". + Prefix for auxiliary variables and constraints. Default: "_sos_reform_" Returns ------- list[str] - List of variable names that were reformulated. + Names of variables that were reformulated. """ - reformulated_vars = [] - - # Get all SOS variables - sos_vars = list(model.variables.sos) + reformulated = [] - for var_name in sos_vars: + for var_name in list(model.variables.sos): var = model.variables[var_name] sos_type = var.attrs.get("sos_type") sos_dim = var.attrs.get("sos_dim") if sos_type is None or sos_dim is None: continue - - # Skip single-element SOS (trivially satisfied) if var.sizes[sos_dim] <= 1: - logger.debug( - f"Skipping SOS{sos_type} reformulation for '{var_name}' (single element)" - ) continue - # Validate bounds - validate_bounds_for_reformulation(var) - - # Check if all bounds are zero (variable already fixed to 0) + # Check if fixed to zero M_upper, M_lower = compute_big_m_values(var) if (M_upper == 0).all() and (M_lower == 0).all(): - logger.debug( - f"Skipping SOS{sos_type} reformulation for '{var_name}' (fixed to zero)" - ) continue - # Perform reformulation based on SOS type if sos_type == 1: reformulate_sos1(model, var, prefix) elif sos_type == 2: reformulate_sos2(model, var, prefix) - # Remove the SOS constraint from the variable model.remove_sos_constraints(var) - reformulated_vars.append(var_name) + reformulated.append(var_name) - logger.info(f"Reformulated {len(reformulated_vars)} SOS constraint(s)") - return reformulated_vars + logger.info(f"Reformulated {len(reformulated)} SOS constraint(s)") + return reformulated diff --git a/test/test_sos_reformulation.py b/test/test_sos_reformulation.py index 17455753..0b1fa9f2 100644 --- a/test/test_sos_reformulation.py +++ b/test/test_sos_reformulation.py @@ -12,20 +12,18 @@ reformulate_all_sos, reformulate_sos1, reformulate_sos2, - validate_bounds_for_reformulation, ) class TestValidateBounds: - """Tests for validate_bounds_for_reformulation.""" + """Tests for bound validation in compute_big_m_values.""" def test_finite_bounds_pass(self) -> None: """Finite bounds should pass validation.""" m = Model() idx = pd.Index([0, 1, 2], name="i") x = m.add_variables(lower=-1, upper=1, coords=[idx], name="x") - # Should not raise - validate_bounds_for_reformulation(x) + compute_big_m_values(x) # Should not raise def test_infinite_upper_bounds_raise(self) -> None: """Infinite upper bounds should raise ValueError.""" @@ -33,7 +31,7 @@ def test_infinite_upper_bounds_raise(self) -> None: idx = pd.Index([0, 1, 2], name="i") x = m.add_variables(lower=0, upper=np.inf, coords=[idx], name="x") with pytest.raises(ValueError, match="infinite upper bounds"): - validate_bounds_for_reformulation(x) + compute_big_m_values(x) def test_infinite_lower_bounds_raise(self) -> None: """Infinite lower bounds should raise ValueError.""" @@ -41,7 +39,7 @@ def test_infinite_lower_bounds_raise(self) -> None: idx = pd.Index([0, 1, 2], name="i") x = m.add_variables(lower=-np.inf, upper=1, coords=[idx], name="x") with pytest.raises(ValueError, match="infinite lower bounds"): - validate_bounds_for_reformulation(x) + compute_big_m_values(x) def test_mixed_infinite_bounds_raise(self) -> None: """Mixed finite/infinite bounds should raise ValueError.""" @@ -54,7 +52,7 @@ def test_mixed_infinite_bounds_raise(self) -> None: name="x", ) with pytest.raises(ValueError, match="infinite lower bounds"): - validate_bounds_for_reformulation(x) + compute_big_m_values(x) class TestComputeBigM: @@ -130,8 +128,7 @@ def test_custom_big_m_allows_infinite_bounds(self) -> None: idx = pd.Index([0, 1, 2], name="i") 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) - # Should not raise - custom big_m bypasses bound validation - validate_bounds_for_reformulation(x) + # Should not raise - custom big_m makes result finite M_upper, M_lower = compute_big_m_values(x) assert np.allclose(M_upper.values, [10, 10, 10]) From 34b00b03015e2b0aac9e7b45c6a7b15ecb46f36b Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 13:05:08 +0100 Subject: [PATCH 07/14] Revert some docs changes to be more surgical --- doc/sos-constraints.rst | 336 +++++++++++++++++++++++++++++++--------- 1 file changed, 261 insertions(+), 75 deletions(-) diff --git a/doc/sos-constraints.rst b/doc/sos-constraints.rst index 5a17c82b..da4d4be0 100644 --- a/doc/sos-constraints.rst +++ b/doc/sos-constraints.rst @@ -3,7 +3,7 @@ Special Ordered Sets (SOS) Constraints ======================================= -SOS constraints model situations where only one (SOS1) or two adjacent (SOS2) variables from an ordered set can be non-zero. +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. .. contents:: :local: @@ -12,60 +12,259 @@ SOS constraints model situations where only one (SOS1) or two adjacent (SOS2) va Overview -------- -- **SOS1**: At most one variable can be non-zero (mutually exclusive choices) -- **SOS2**: At most two adjacent variables can be non-zero (piecewise linear functions) +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. The adjacency is determined by the ordering weights (coordinates) of the variables. + +**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: + .. code-block:: python import linopy import pandas as pd + import xarray as xr + # Create model m = linopy.Model() - # SOS1: at most one option selected - options = pd.Index([0, 1, 2], name="options") - x = m.add_variables(coords=[options], name="x", lower=0, upper=1) + # 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") - # SOS2: at most two adjacent breakpoints active - breakpoints = pd.Index([0.0, 1.0, 2.0], name="bp") + # 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="bp") + m.add_sos_constraints(lambdas, sos_type=2, sos_dim="breakpoints") + +Method Signature +~~~~~~~~~~~~~~~~ + +.. code-block:: python + + Model.add_sos_constraints(variable, sos_type, sos_dim, big_m=None) + +**Parameters:** + +- ``variable`` : Variable + The variable to which the SOS constraint should be applied +- ``sos_type`` : {1, 2} + Type of SOS constraint (1 or 2) +- ``sos_dim`` : str + Name of the dimension along which the SOS constraint applies +- ``big_m`` : float | tuple[float, float] | None + Custom Big-M values for reformulation (see :ref:`sos-reformulation`) **Requirements:** -- The SOS dimension must have numeric coordinates (used as ordering weights) -- Only one SOS constraint per variable +- The specified dimension must exist in the variable +- The coordinates for the SOS dimension must be numeric (used as weights for ordering) +- Only one SOS constraint can be applied per variable + +Examples +-------- + +Example 1: Facility Location (SOS1) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. code-block:: python + + 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) -Multi-dimensional Variables -~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # SOS1 constraint: at most one facility can be built + m.add_sos_constraints(build, sos_type=1, sos_dim="locations") -For multi-dimensional variables, the SOS constraint applies independently for each combination of non-SOS dimensions: + # 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) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code-block:: python - periods = pd.Index([0, 1], name="periods") + 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. + +.. code-block:: python + + # Multi-period production planning + periods = pd.Index(range(3), name="periods") modes = pd.Index([0, 1, 2], name="modes") - x = m.add_variables(lower=0, upper=1, coords=[periods, modes], name="x") - m.add_sos_constraints(x, sos_type=1, sos_dim="modes") - # Result: at most one mode selected PER period + + # 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: + +.. code-block:: python + + # 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: + +.. code-block:: python + + 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: + +.. code-block:: text + + sos + + s0: S1 :: x0:0 x1:1 x2:2 + s3: S2 :: x3:0.0 x4:1.0 x5:2.0 Solver Compatibility -------------------- -**Native SOS support:** Gurobi, CPLEX, CBC, SCIP, Xpress +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:** -**No SOS support (use reformulation):** HiGHS, GLPK, MOSEK +- Gurobi (via ``gurobipy``) + +**Unsupported solvers:** + +- HiGHS (does not support SOS constraints) +- GLPK +- 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. +For solvers without native SOS support, linopy can reformulate SOS constraints +as binary + linear constraints using the Big-M method. .. code-block:: python @@ -76,18 +275,16 @@ For solvers without native SOS support, linopy can reformulate SOS constraints a m.reformulate_sos_constraints() m.solve(solver_name="highs") -Big-M Values -~~~~~~~~~~~~ - -Big-M values are derived from variable bounds by default. For infinite bounds, specify custom values: +**Requirements:** Big-M values are derived from variable bounds. For infinite bounds, +specify custom values via the ``big_m`` parameter: .. code-block:: python - # Finite bounds: Big-M = bounds (default) + # 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 bounds: specify Big-M explicitly + # Infinite 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) @@ -96,69 +293,58 @@ Big-M values are derived from variable bounds by default. For infinite bounds, s The reformulation uses the tighter of ``big_m`` and variable bounds. -Mathematical Formulation -~~~~~~~~~~~~~~~~~~~~~~~~ - -**SOS1:** Binary indicators :math:`y_i`, linking constraints, cardinality :math:`\sum y_i \leq 1` - -.. math:: +Common Patterns +--------------- - y_i \in \{0, 1\}, \quad x_i \leq U_i \cdot y_i, \quad x_i \geq L_i \cdot y_i, \quad \sum y_i \leq 1 - -**SOS2:** Segment indicators :math:`z_j` for :math:`j = 0, \ldots, n-2`, adjacency constraints - -.. math:: - - z_j \in \{0, 1\}, \quad x_0 \leq U_0 \cdot z_0, \quad x_i \leq U_i(z_{i-1} + z_i), \quad x_{n-1} \leq U_{n-1} \cdot z_{n-2}, \quad \sum z_j \leq 1 - -Example: Piecewise Linear with HiGHS -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Piecewise Linear Cost Function +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code-block:: python - import linopy - import pandas as pd + 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") - m = linopy.Model() + 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") - # Approximate f(x) = x² over [0, 3] - bp = pd.Index([0.0, 1.0, 2.0, 3.0], name="bp") - x_vals, y_vals = bp.to_numpy(), bp.to_numpy() ** 2 + cost_var = model.add_variables(name="cost", lower=0) - lambdas = m.add_variables(lower=0, upper=1, coords=[bp], name="lambdas") - m.add_sos_constraints(lambdas, sos_type=2, sos_dim="bp") + x_vals = xr.DataArray(breakpoints, coords=[lambda_coords]) + c_vals = xr.DataArray(costs, coords=[lambda_coords]) - x = m.add_variables(name="x", lower=0, upper=3) - y = m.add_variables(name="y", lower=0, upper=9) + 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") - m.add_constraints(lambdas.sum() == 1, name="convexity") - m.add_constraints(x == (lambdas * x_vals).sum(), name="x_interp") - m.add_constraints(y == (lambdas * y_vals).sum(), name="y_interp") - m.add_constraints(x >= 1.5, name="x_min") - m.add_objective(y) + return cost_var - m.solve(solver_name="highs", reformulate_sos=True) - -API Reference -------------- - -.. py:method:: Model.add_sos_constraints(variable, sos_type, sos_dim, big_m=None) +Mutually Exclusive Investments +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - Add SOS constraint to a variable. +.. code-block:: python - :param variable: Variable to constrain - :param sos_type: 1 (at most one non-zero) or 2 (at most two adjacent non-zero) - :param sos_dim: Dimension for SOS ordering - :param big_m: Custom Big-M: ``float`` (symmetric) or ``tuple[float, float]`` (upper, lower) + def add_exclusive_investments(model, projects, costs, returns): + """Add mutually exclusive investment decisions using SOS1.""" + project_coords = pd.Index(projects, name="projects") -.. py:method:: Model.remove_sos_constraints(variable) + invest = model.add_variables( + coords=[project_coords], name="invest", binary=True + ) + model.add_sos_constraints(invest, sos_type=1, sos_dim="projects") - Remove SOS constraints from a variable. + total_cost = (invest * costs).sum() + total_return = (invest * returns).sum() -.. py:method:: Model.reformulate_sos_constraints(prefix="_sos_reform_") + return invest, total_cost, total_return - Convert SOS to binary + linear constraints. Returns list of reformulated variable names. -.. py:attribute:: Variables.sos +See Also +-------- - All variables with SOS constraints. +- :doc:`creating-variables`: Creating variables with coordinates +- :doc:`creating-constraints`: Adding regular constraints +- :doc:`user-guide`: General linopy usage patterns From da1dc61349adf13a2f958fabf3d71bcc68e4b9c2 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 13:07:56 +0100 Subject: [PATCH 08/14] Add math to docs --- doc/sos-constraints.rst | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/doc/sos-constraints.rst b/doc/sos-constraints.rst index da4d4be0..84eecb06 100644 --- a/doc/sos-constraints.rst +++ b/doc/sos-constraints.rst @@ -293,6 +293,34 @@ specify custom values via the ``big_m`` parameter: The reformulation uses the tighter of ``big_m`` and variable bounds. +Mathematical Formulation +~~~~~~~~~~~~~~~~~~~~~~~~ + +**SOS1 Reformulation:** +Given variables :math:`x_i` with bounds :math:`L_i \leq x_i \leq U_i`, the constraint +"at most one :math:`x_i` is non-zero" becomes: + +.. math:: + + y_i &\in \{0, 1\} \quad \forall i \\ + x_i &\leq U_i \cdot y_i \quad \forall i \text{ where } U_i > 0 \\ + x_i &\geq L_i \cdot y_i \quad \forall i \text{ where } L_i < 0 \\ + \sum_i y_i &\leq 1 + +**SOS2 Reformulation:** +Given ordered variables :math:`x_0, \ldots, x_{n-1}`, the constraint +"at most two adjacent variables are non-zero" uses segment indicators :math:`z_j`: + +.. math:: + + z_j &\in \{0, 1\} \quad j = 0, \ldots, n-2 \\ + x_0 &\leq U_0 \cdot z_0 \\ + x_i &\leq U_i \cdot (z_{i-1} + z_i) \quad i = 1, \ldots, n-2 \\ + x_{n-1} &\leq U_{n-1} \cdot z_{n-2} \\ + \sum_j z_j &\leq 1 + +(Analogous constraints for lower bounds when :math:`L_i < 0`.) + Common Patterns --------------- From 9b0744c6a9c7c41719b33e27c3543db5a53a7296 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 13:22:49 +0100 Subject: [PATCH 09/14] Improve docs --- doc/sos-constraints.rst | 40 ++++++++++++++++++++++++++-------------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/doc/sos-constraints.rst b/doc/sos-constraints.rst index 84eecb06..1e3bedbe 100644 --- a/doc/sos-constraints.rst +++ b/doc/sos-constraints.rst @@ -297,29 +297,41 @@ Mathematical Formulation ~~~~~~~~~~~~~~~~~~~~~~~~ **SOS1 Reformulation:** -Given variables :math:`x_i` with bounds :math:`L_i \leq x_i \leq U_i`, the constraint -"at most one :math:`x_i` is non-zero" becomes: + +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:: - y_i &\in \{0, 1\} \quad \forall i \\ - x_i &\leq U_i \cdot y_i \quad \forall i \text{ where } U_i > 0 \\ - x_i &\geq L_i \cdot y_i \quad \forall i \text{ where } L_i < 0 \\ - \sum_i y_i &\leq 1 + 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:** -Given ordered variables :math:`x_0, \ldots, x_{n-1}`, the constraint -"at most two adjacent variables are non-zero" uses segment indicators :math:`z_j`: + +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:: - z_j &\in \{0, 1\} \quad j = 0, \ldots, n-2 \\ - x_0 &\leq U_0 \cdot z_0 \\ - x_i &\leq U_i \cdot (z_{i-1} + z_i) \quad i = 1, \ldots, n-2 \\ - x_{n-1} &\leq U_{n-1} \cdot z_{n-2} \\ - \sum_j z_j &\leq 1 + 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\} -(Analogous constraints for lower bounds when :math:`L_i < 0`.) +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 --------------- From 7bac84d8edeb7282803274a564f0899432ab31c1 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 13:29:59 +0100 Subject: [PATCH 10/14] Code simplifications: MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 1. sos_reformulation.py (230 → 203 lines): - compute_big_m_values now returns single DataArray (not tuple) - Removed all lower bound handling - only supports non-negative variables - Removed add_linking_constraints helper function - Simplified SOS1/SOS2 to only add upper linking constraints 2. model.py: - Simplified big_m parameter from float | tuple[float, float] | None to float | None - Removed big_m_lower attribute handling 3. Documentation (sos-constraints.rst): - Updated big_m type signature - Removed asymmetric Big-M example - Added explicit requirement that variables must have non-negative lower bounds 4. Tests (46 → 38 tests): - Removed tests for negative bounds - Removed tests for tuple big_m - Added tests for negative lower bound validation error Rationale: The mathematical formulation in the docs assumes x ∈ ℝⁿ₊ (non-negative reals). This matches 99%+ of SOS use cases (selection indicators, piecewise linear weights). The simplified code is now consistent with the documented formulation. --- doc/sos-constraints.rst | 18 ++-- linopy/model.py | 32 ++----- linopy/sos_reformulation.py | 118 +++++++++-------------- test/test_sos_reformulation.py | 170 +++++++-------------------------- 4 files changed, 97 insertions(+), 241 deletions(-) diff --git a/doc/sos-constraints.rst b/doc/sos-constraints.rst index 1e3bedbe..a2731400 100644 --- a/doc/sos-constraints.rst +++ b/doc/sos-constraints.rst @@ -85,8 +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 | tuple[float, float] | None - Custom Big-M values for reformulation (see :ref:`sos-reformulation`) +- ``big_m`` : float | None + Custom Big-M value for reformulation (see :ref:`sos-reformulation`) **Requirements:** @@ -275,8 +275,11 @@ as binary + linear constraints using the Big-M method. m.reformulate_sos_constraints() m.solve(solver_name="highs") -**Requirements:** Big-M values are derived from variable bounds. For infinite bounds, -specify custom values via the ``big_m`` parameter: +**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 @@ -284,14 +287,11 @@ specify custom values via the ``big_m`` parameter: x = m.add_variables(lower=0, upper=100, coords=[idx], name="x") m.add_sos_constraints(x, sos_type=1, sos_dim="i") - # Infinite bounds: specify Big-M + # 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) - # Asymmetric Big-M - m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=(10, -5)) - -The reformulation uses the tighter of ``big_m`` and variable bounds. +The reformulation uses the tighter of ``big_m`` and variable upper bound. Mathematical Formulation ~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/linopy/model.py b/linopy/model.py index c7d0b876..9a8bebdd 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -557,7 +557,7 @@ def add_sos_constraints( variable: Variable, sos_type: Literal[1, 2], sos_dim: str, - big_m: float | tuple[float, float] | None = None, + big_m: float | None = None, ) -> None: """ Add an sos1 or sos2 constraint for one dimension of a variable @@ -571,16 +571,15 @@ def add_sos_constraints( Type of SOS sos_dim : str Which dimension of variable to add SOS constraint to - big_m : float | tuple[float, float] | None, optional - Big-M value(s) for SOS reformulation. Only used when reformulating + 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 bounds as Big-M - - float: Symmetric Big-M (M_upper = big_m, M_lower = -big_m) - - tuple (upper, lower): Asymmetric Big-M values + - None (default): Use variable upper bounds as Big-M + - float: Custom Big-M value - The reformulation uses the tighter of big_m and variable bounds: - M_upper = min(big_m_upper, var.upper), M_lower = max(big_m_lower, var.lower). + 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. """ @@ -603,20 +602,10 @@ def add_sos_constraints( f"but got {variable.coords[sos_dim].dtype}" ) - # Process and store big_m values + # Process and store big_m value attrs_update: dict[str, Any] = {"sos_type": sos_type, "sos_dim": sos_dim} if big_m is not None: - if isinstance(big_m, tuple): - if len(big_m) != 2: - raise ValueError( - "big_m tuple must have exactly 2 elements (upper, lower)" - ) - attrs_update["big_m_upper"] = float(big_m[0]) - attrs_update["big_m_lower"] = float(big_m[1]) - else: - # Scalar: symmetric (upper = big_m, lower = -big_m) - attrs_update["big_m_upper"] = float(big_m) - attrs_update["big_m_lower"] = -float(big_m) + attrs_update["big_m_upper"] = float(big_m) variable.attrs.update(attrs_update) @@ -867,9 +856,8 @@ def remove_sos_constraints(self, variable: Variable) -> None: del variable.attrs["sos_type"], variable.attrs["sos_dim"] - # Also remove big_m attributes if present + # Also remove big_m attribute if present variable.attrs.pop("big_m_upper", None) - variable.attrs.pop("big_m_lower", None) logger.debug( f"Removed sos{sos_type} constraint on {sos_dim} from {variable.name}" diff --git a/linopy/sos_reformulation.py b/linopy/sos_reformulation.py index 37f0d963..3ed09858 100644 --- a/linopy/sos_reformulation.py +++ b/linopy/sos_reformulation.py @@ -21,39 +21,40 @@ logger = logging.getLogger(__name__) -def compute_big_m_values(var: Variable) -> tuple[DataArray, DataArray]: +def compute_big_m_values(var: Variable) -> DataArray: """ - Compute Big-M values from variable bounds and custom big_m attributes. + Compute Big-M values from variable bounds and custom big_m attribute. - Uses the tighter of variable bounds and custom big_m values to ensure + Uses the tighter of variable upper bound and custom big_m to ensure the best possible LP relaxation. Parameters ---------- var : Variable - Variable with bounds (and optionally big_m_upper/big_m_lower attrs). + Variable with bounds (and optionally big_m_upper attr). Returns ------- - tuple[DataArray, DataArray] - (M_upper, M_lower) for reformulation constraints: - x <= M_upper * y and x >= M_lower * y + DataArray + M_upper for reformulation constraints: x <= M_upper * y Raises ------ ValueError - If resulting Big-M values are infinite. + If variable has negative lower bounds or infinite upper bounds. """ - big_m_upper = var.attrs.get("big_m_upper") - big_m_lower = var.attrs.get("big_m_lower") + # SOS reformulation requires non-negative variables + if (var.lower < 0).any(): + raise ValueError( + f"Variable '{var.name}' has negative lower bounds. " + "SOS reformulation requires non-negative variables (lower >= 0)." + ) + big_m_upper = var.attrs.get("big_m_upper") M_upper = var.upper - M_lower = var.lower if big_m_upper is not None: M_upper = np.minimum(M_upper, big_m_upper) - if big_m_lower is not None: - M_lower = np.maximum(M_lower, big_m_lower) # Validate finiteness if np.isinf(M_upper).any(): @@ -61,46 +62,36 @@ def compute_big_m_values(var: Variable) -> tuple[DataArray, DataArray]: f"Variable '{var.name}' has infinite upper bounds. " "Set finite bounds or specify big_m in add_sos_constraints()." ) - if np.isinf(M_lower).any(): - raise ValueError( - f"Variable '{var.name}' has infinite lower bounds. " - "Set finite bounds or specify big_m in add_sos_constraints()." - ) - return M_upper, M_lower + return M_upper def reformulate_sos1(model: Model, var: Variable, prefix: str) -> None: """ Reformulate SOS1 constraint as binary + linear constraints. - For each x[i] with bounds [L[i], U[i]]: + For each x[i] with upper bound M[i]: - Add binary indicator y[i] - - x[i] <= U[i] * y[i] (upper linking, if U > 0) - - x[i] >= L[i] * y[i] (lower linking, if L < 0) - - sum(y) <= 1 (cardinality) + - x[i] <= M[i] * y[i] + - sum(y) <= 1 Parameters ---------- model : Model Model to add reformulation constraints to. var : Variable - Variable with SOS1 constraint. + Variable with SOS1 constraint (must have non-negative lower bounds). prefix : str Prefix for naming auxiliary variables and constraints. """ sos_dim = var.attrs["sos_dim"] name = var.name - M_upper, M_lower = compute_big_m_values(var) + M = compute_big_m_values(var) coords = [var.coords[d] for d in var.dims] y = model.add_variables(coords=coords, name=f"{prefix}{name}_y", binary=True) - if (M_upper > 0).any(): - model.add_constraints(var <= M_upper * y, name=f"{prefix}{name}_upper") - if (M_lower < 0).any(): - model.add_constraints(var >= M_lower * y, name=f"{prefix}{name}_lower") - + model.add_constraints(var <= M * y, name=f"{prefix}{name}_upper") model.add_constraints(y.sum(dim=sos_dim) <= 1, name=f"{prefix}{name}_card") @@ -108,12 +99,11 @@ def reformulate_sos2(model: Model, var: Variable, prefix: str) -> None: """ Reformulate SOS2 constraint as binary + linear constraints. - For ordered x[0..n-1]: + For ordered x[0..n-1] with upper bounds M[i]: - Add n-1 binary segment indicators z[i] - - x[0] <= U[0] * z[0] - - x[i] <= U[i] * (z[i-1] + z[i]) for middle elements - - x[n-1] <= U[n-1] * z[n-2] - - Similar for lower bounds if L < 0 + - x[0] <= M[0] * z[0] + - x[i] <= M[i] * (z[i-1] + z[i]) for middle elements + - x[n-1] <= M[n-1] * z[n-2] - sum(z) <= 1 Parameters @@ -121,7 +111,7 @@ def reformulate_sos2(model: Model, var: Variable, prefix: str) -> None: model : Model Model to add reformulation constraints to. var : Variable - Variable with SOS2 constraint. + Variable with SOS2 constraint (must have non-negative lower bounds). prefix : str Prefix for naming auxiliary variables and constraints. """ @@ -132,7 +122,7 @@ def reformulate_sos2(model: Model, var: Variable, prefix: str) -> None: if n <= 1: return - M_upper, M_lower = compute_big_m_values(var) + M = compute_big_m_values(var) # Create n-1 segment indicators z_coords = [ @@ -146,41 +136,25 @@ def reformulate_sos2(model: Model, var: Variable, prefix: str) -> None: # Convert to expressions to avoid SOS attr validation on isel x, z_expr = 1 * var, 1 * z - def add_linking_constraints(M: DataArray, sign: str, suffix: str) -> None: - """Add x <= M*z or x >= M*z constraints for first/middle/last elements.""" - if sign == "upper" and not (M > 0).any(): - return - if sign == "lower" and not (M < 0).any(): - return - - op = (lambda a, b: a <= b) if sign == "upper" else (lambda a, b: a >= b) - - # First: x[0] op M[0] * z[0] + # First: x[0] <= M[0] * z[0] + model.add_constraints( + x.isel({sos_dim: 0}) <= M.isel({sos_dim: 0}) * z_expr.isel({sos_dim: 0}), + name=f"{prefix}{name}_upper_first", + ) + # Middle: x[i] <= M[i] * (z[i-1] + z[i]) + for i in range(1, n - 1): model.add_constraints( - op(x.isel({sos_dim: 0}), M.isel({sos_dim: 0}) * z_expr.isel({sos_dim: 0})), - name=f"{prefix}{name}_{suffix}_first", + x.isel({sos_dim: i}) + <= M.isel({sos_dim: i}) + * (z_expr.isel({sos_dim: i - 1}) + z_expr.isel({sos_dim: i})), + name=f"{prefix}{name}_upper_mid_{i}", ) - # Middle: x[i] op M[i] * (z[i-1] + z[i]) - for i in range(1, n - 1): - model.add_constraints( - op( - x.isel({sos_dim: i}), - M.isel({sos_dim: i}) - * (z_expr.isel({sos_dim: i - 1}) + z_expr.isel({sos_dim: i})), - ), - name=f"{prefix}{name}_{suffix}_mid_{i}", - ) - # Last: x[n-1] op M[n-1] * z[n-2] - model.add_constraints( - op( - x.isel({sos_dim: n - 1}), - M.isel({sos_dim: n - 1}) * z_expr.isel({sos_dim: n - 2}), - ), - name=f"{prefix}{name}_{suffix}_last", - ) - - add_linking_constraints(M_upper, "upper", "upper") - add_linking_constraints(M_lower, "lower", "lower") + # Last: x[n-1] <= M[n-1] * z[n-2] + model.add_constraints( + x.isel({sos_dim: n - 1}) + <= M.isel({sos_dim: n - 1}) * z_expr.isel({sos_dim: n - 2}), + name=f"{prefix}{name}_upper_last", + ) model.add_constraints(z.sum(dim=sos_dim) <= 1, name=f"{prefix}{name}_card") @@ -214,8 +188,8 @@ def reformulate_all_sos(model: Model, prefix: str = "_sos_reform_") -> list[str] continue # Check if fixed to zero - M_upper, M_lower = compute_big_m_values(var) - if (M_upper == 0).all() and (M_lower == 0).all(): + M = compute_big_m_values(var) + if (M == 0).all(): continue if sos_type == 1: diff --git a/test/test_sos_reformulation.py b/test/test_sos_reformulation.py index 0b1fa9f2..a5301179 100644 --- a/test/test_sos_reformulation.py +++ b/test/test_sos_reformulation.py @@ -19,10 +19,10 @@ class TestValidateBounds: """Tests for bound validation in compute_big_m_values.""" def test_finite_bounds_pass(self) -> None: - """Finite bounds should pass validation.""" + """Finite non-negative bounds should pass validation.""" m = Model() idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=-1, upper=1, coords=[idx], name="x") + x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") compute_big_m_values(x) # Should not raise def test_infinite_upper_bounds_raise(self) -> None: @@ -33,25 +33,25 @@ def test_infinite_upper_bounds_raise(self) -> None: with pytest.raises(ValueError, match="infinite upper bounds"): compute_big_m_values(x) - def test_infinite_lower_bounds_raise(self) -> None: - """Infinite lower bounds should raise ValueError.""" + def test_negative_lower_bounds_raise(self) -> None: + """Negative lower bounds should raise ValueError.""" m = Model() idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=-np.inf, upper=1, coords=[idx], name="x") - with pytest.raises(ValueError, match="infinite lower bounds"): + x = m.add_variables(lower=-1, upper=1, coords=[idx], name="x") + with pytest.raises(ValueError, match="negative lower bounds"): compute_big_m_values(x) - def test_mixed_infinite_bounds_raise(self) -> None: - """Mixed finite/infinite bounds should raise ValueError.""" + def test_mixed_negative_lower_bounds_raise(self) -> None: + """Mixed finite/negative lower bounds should raise ValueError.""" m = Model() idx = pd.Index([0, 1, 2], name="i") x = m.add_variables( - lower=np.array([0, -np.inf, 0]), + lower=np.array([0, -1, 0]), upper=np.array([1, 1, 1]), coords=[idx], name="x", ) - with pytest.raises(ValueError, match="infinite lower bounds"): + with pytest.raises(ValueError, match="negative lower bounds"): compute_big_m_values(x) @@ -63,41 +63,21 @@ def test_positive_bounds(self) -> None: m = Model() idx = pd.Index([0, 1, 2], name="i") x = m.add_variables(lower=0, upper=10, coords=[idx], name="x") - M_upper, M_lower = compute_big_m_values(x) - assert np.allclose(M_upper.values, [10, 10, 10]) - assert np.allclose(M_lower.values, [0, 0, 0]) - - def test_negative_bounds(self) -> None: - """Test Big-M computation with negative bounds.""" - m = Model() - idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=-10, upper=-1, coords=[idx], name="x") - M_upper, M_lower = compute_big_m_values(x) - assert np.allclose(M_upper.values, [-1, -1, -1]) - assert np.allclose(M_lower.values, [-10, -10, -10]) - - def test_mixed_bounds(self) -> None: - """Test Big-M computation with mixed bounds.""" - m = Model() - idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=-5, upper=10, coords=[idx], name="x") - M_upper, M_lower = compute_big_m_values(x) - assert np.allclose(M_upper.values, [10, 10, 10]) - assert np.allclose(M_lower.values, [-5, -5, -5]) + M = compute_big_m_values(x) + assert np.allclose(M.values, [10, 10, 10]) def test_varying_bounds(self) -> None: - """Test Big-M computation with varying bounds.""" + """Test Big-M computation with varying upper bounds.""" m = Model() idx = pd.Index([0, 1, 2], name="i") x = m.add_variables( - lower=np.array([-1, -2, -3]), + lower=np.array([0, 0, 0]), upper=np.array([1, 2, 3]), coords=[idx], name="x", ) - M_upper, M_lower = compute_big_m_values(x) - assert np.allclose(M_upper.values, [1, 2, 3]) - assert np.allclose(M_lower.values, [-1, -2, -3]) + M = compute_big_m_values(x) + assert np.allclose(M.values, [1, 2, 3]) def test_custom_big_m_scalar(self) -> None: """Test Big-M uses tighter of custom value and bounds.""" @@ -105,22 +85,9 @@ def test_custom_big_m_scalar(self) -> None: idx = pd.Index([0, 1, 2], name="i") x = m.add_variables(lower=0, upper=100, coords=[idx], name="x") m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=10) - M_upper, M_lower = compute_big_m_values(x) - # big_m=10 gives big_m_upper=10, big_m_lower=-10 - # M_upper = min(10, 100) = 10 (custom is tighter) - # M_lower = max(-10, 0) = 0 (bound is tighter) - assert np.allclose(M_upper.values, [10, 10, 10]) - assert np.allclose(M_lower.values, [0, 0, 0]) - - def test_custom_big_m_tuple(self) -> None: - """Test Big-M with custom tuple (upper, lower).""" - m = Model() - idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=-100, upper=100, coords=[idx], name="x") - m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=(5, -3)) - M_upper, M_lower = compute_big_m_values(x) - assert np.allclose(M_upper.values, [5, 5, 5]) - assert np.allclose(M_lower.values, [-3, -3, -3]) + M = compute_big_m_values(x) + # M = min(10, 100) = 10 (custom is tighter) + assert np.allclose(M.values, [10, 10, 10]) def test_custom_big_m_allows_infinite_bounds(self) -> None: """Test that custom big_m allows variables with infinite bounds.""" @@ -129,8 +96,8 @@ def test_custom_big_m_allows_infinite_bounds(self) -> None: 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) # Should not raise - custom big_m makes result finite - M_upper, M_lower = compute_big_m_values(x) - assert np.allclose(M_upper.values, [10, 10, 10]) + M = compute_big_m_values(x) + assert np.allclose(M.values, [10, 10, 10]) class TestSOS1Reformulation: @@ -157,36 +124,6 @@ def test_basic_sos1(self) -> None: assert y.dims == x.dims assert y.sizes == x.sizes - def test_sos1_negative_bounds(self) -> None: - """Test SOS1 reformulation with negative bounds.""" - m = Model() - idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=-2, upper=1, coords=[idx], name="x") - m.add_sos_constraints(x, sos_type=1, sos_dim="i") - - reformulate_sos1(m, x, "_test_") - m.remove_sos_constraints(x) - - # Should have both upper and lower constraints - assert "_test_x_upper" in m.constraints - assert "_test_x_lower" in m.constraints - assert "_test_x_card" in m.constraints - - def test_sos1_all_negative_bounds(self) -> None: - """Test SOS1 reformulation with all negative bounds (no upper constraint).""" - m = Model() - idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=-2, upper=-1, coords=[idx], name="x") - m.add_sos_constraints(x, sos_type=1, sos_dim="i") - - reformulate_sos1(m, x, "_test_") - m.remove_sos_constraints(x) - - # Should only have lower constraint (no positive upper bounds) - assert "_test_x_upper" not in m.constraints - assert "_test_x_lower" in m.constraints - assert "_test_x_card" in m.constraints - def test_sos1_multidimensional(self) -> None: """Test SOS1 reformulation with multi-dimensional variables.""" m = Model() @@ -205,7 +142,6 @@ def test_sos1_multidimensional(self) -> None: # Cardinality constraint should have reduced dimensions (summed over i) card_con = m.constraints["_test_x_card"] assert "j" in card_con.dims - # The constraint has j dimension preserved class TestSOS2Reformulation: @@ -273,20 +209,6 @@ def test_sos2_with_middle_constraints(self) -> None: assert "_test_x_upper_mid_2" in m.constraints # middle element at index 2 assert "_test_x_upper_last" in m.constraints - def test_sos2_negative_bounds(self) -> None: - """Test SOS2 reformulation with negative bounds.""" - m = Model() - idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=-2, upper=1, coords=[idx], name="x") - m.add_sos_constraints(x, sos_type=2, sos_dim="i") - - reformulate_sos2(m, x, "_test_") - m.remove_sos_constraints(x) - - # Should have both upper and lower constraints - assert "_test_x_upper_first" in m.constraints - assert "_test_x_lower_first" in m.constraints - def test_sos2_multidimensional(self) -> None: """Test SOS2 reformulation with multi-dimensional variables.""" m = Model() @@ -372,6 +294,16 @@ def test_reformulate_raises_on_infinite_bounds(self) -> None: with pytest.raises(ValueError, match="infinite"): reformulate_all_sos(m) + def test_reformulate_raises_on_negative_lower_bounds(self) -> None: + """Test that negative lower bounds raise ValueError.""" + m = Model() + idx = pd.Index([0, 1, 2], name="i") + x = m.add_variables(lower=-1, upper=1, coords=[idx], name="x") + m.add_sos_constraints(x, sos_type=1, sos_dim="i") + + with pytest.raises(ValueError, match="negative lower bounds"): + reformulate_all_sos(m) + class TestModelReformulateSOS: """Tests for Model.reformulate_sos_constraints method.""" @@ -576,19 +508,6 @@ def test_sos2_equivalence(self) -> None: class TestEdgeCases: """Tests for edge cases.""" - def test_zero_lower_bound(self) -> None: - """Test handling of zero lower bound.""" - m = Model() - idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") - m.add_sos_constraints(x, sos_type=1, sos_dim="i") - - reformulate_all_sos(m) - - # Should only have upper constraints (no negative lower bounds) - assert "_sos_reform_x_upper" in m.constraints - assert "_sos_reform_x_lower" not in m.constraints - def test_preserves_non_sos_variables(self) -> None: """Test that non-SOS variables are preserved.""" m = Model() @@ -649,19 +568,17 @@ def test_float_coordinates(self) -> None: assert z.sizes["bp"] == 2 def test_custom_big_m_removed_on_remove_sos(self) -> None: - """Test that custom big_m attributes are removed with SOS constraint.""" + """Test that custom big_m attribute is removed with SOS constraint.""" m = Model() idx = pd.Index([0, 1, 2], name="i") x = m.add_variables(lower=0, upper=100, coords=[idx], name="x") m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=10) assert "big_m_upper" in x.attrs - assert "big_m_lower" in x.attrs m.remove_sos_constraints(x) assert "big_m_upper" not in x.attrs - assert "big_m_lower" not in x.attrs @pytest.mark.skipif("highs" not in available_solvers, reason="HiGHS not installed") @@ -697,26 +614,3 @@ def test_solve_with_infinite_bounds_and_custom_big_m(self) -> None: assert m.objective.value is not None # Max is 15 (x[2]=5) assert np.isclose(m.objective.value, 15, atol=1e-5) - - def test_custom_big_m_tuple(self) -> None: - """Test solving with asymmetric big_m tuple.""" - m = Model() - idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=-100, upper=100, coords=[idx], name="x") - m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=(2, -1)) - m.add_objective(x * np.array([1, 2, 3]), sense="max") - - m.solve(solver_name="highs", reformulate_sos=True) - - # Max is 6 (x[2]=2, limited by big_m_upper=2) - assert m.objective.value is not None - assert np.isclose(m.objective.value, 6, atol=1e-5) - - def test_big_m_invalid_tuple_length(self) -> None: - """Test that invalid tuple length raises error.""" - m = Model() - idx = pd.Index([0, 1, 2], name="i") - x = m.add_variables(lower=0, upper=1, coords=[idx], name="x") - - with pytest.raises(ValueError, match="exactly 2 elements"): - m.add_sos_constraints(x, sos_type=1, sos_dim="i", big_m=(1, 2, 3)) From 33f714878535b105094c1fa46e56ebed8cc26a47 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 13:41:19 +0100 Subject: [PATCH 11/14] Fix mypy --- linopy/sos_reformulation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/linopy/sos_reformulation.py b/linopy/sos_reformulation.py index 3ed09858..de730708 100644 --- a/linopy/sos_reformulation.py +++ b/linopy/sos_reformulation.py @@ -54,7 +54,7 @@ def compute_big_m_values(var: Variable) -> DataArray: M_upper = var.upper if big_m_upper is not None: - M_upper = np.minimum(M_upper, big_m_upper) + M_upper = M_upper.clip(max=big_m_upper) # type: ignore[arg-type] # Validate finiteness if np.isinf(M_upper).any(): @@ -84,7 +84,7 @@ def reformulate_sos1(model: Model, var: Variable, prefix: str) -> None: prefix : str Prefix for naming auxiliary variables and constraints. """ - sos_dim = var.attrs["sos_dim"] + sos_dim = str(var.attrs["sos_dim"]) name = var.name M = compute_big_m_values(var) @@ -115,7 +115,7 @@ def reformulate_sos2(model: Model, var: Variable, prefix: str) -> None: prefix : str Prefix for naming auxiliary variables and constraints. """ - sos_dim = var.attrs["sos_dim"] + sos_dim = str(var.attrs["sos_dim"]) name = var.name n = var.sizes[sos_dim] From cfd7c062c4f97f93a82a129a398f4a1bfcc269be Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 13:42:54 +0100 Subject: [PATCH 12/14] Fix mypy --- test/test_sos_reformulation.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/test_sos_reformulation.py b/test/test_sos_reformulation.py index a5301179..c6445b03 100644 --- a/test/test_sos_reformulation.py +++ b/test/test_sos_reformulation.py @@ -354,6 +354,7 @@ def test_sos1_minimize_with_highs(self) -> None: m.solve(solver_name="highs", reformulate_sos=True) # Should minimize to 0 by setting all x = 0 + assert m.objective.value is not None assert np.isclose(m.objective.value, 0, atol=1e-5) def test_sos2_maximize_with_highs(self) -> None: @@ -476,6 +477,8 @@ def test_sos1_equivalence(self) -> None: if "highs" in available_solvers: m2.solve(solver_name="highs", reformulate_sos=True) + assert m1.objective.value is not None + assert m2.objective.value is not None assert np.isclose(m1.objective.value, m2.objective.value, atol=1e-5) def test_sos2_equivalence(self) -> None: @@ -502,6 +505,8 @@ def test_sos2_equivalence(self) -> None: if "highs" in available_solvers: m2.solve(solver_name="highs", reformulate_sos=True) + assert m1.objective.value is not None + assert m2.objective.value is not None assert np.isclose(m1.objective.value, m2.objective.value, atol=1e-5) From a1f5c5d6283c54c04acf35e422baebbf2510554e Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 19:07:33 +0100 Subject: [PATCH 13/14] Add constants for sos attr keys --- linopy/constants.py | 5 +++++ linopy/io.py | 10 +++++----- linopy/model.py | 23 +++++++++++++---------- linopy/sos_reformulation.py | 12 +++++++----- linopy/variables.py | 20 ++++++++++---------- 5 files changed, 40 insertions(+), 30 deletions(-) diff --git a/linopy/constants.py b/linopy/constants.py index 021a9a10..2e1ef47a 100644 --- a/linopy/constants.py +++ b/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): """ diff --git a/linopy/io.py b/linopy/io.py index 56fe033d..7dc4177e 100644 --- a/linopy/io.py +++ b/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() diff --git a/linopy/model.py b/linopy/model.py index 9a8bebdd..8828d819 100644 --- a/linopy/model.py +++ b/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, @@ -588,9 +591,9 @@ def add_sos_constraints( 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}" ) @@ -603,9 +606,9 @@ def add_sos_constraints( ) # Process and store big_m value - attrs_update: dict[str, Any] = {"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: - attrs_update["big_m_upper"] = float(big_m) + attrs_update[SOS_BIG_M_ATTR] = float(big_m) variable.attrs.update(attrs_update) @@ -848,16 +851,16 @@ 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("big_m_upper", None) + variable.attrs.pop(SOS_BIG_M_ATTR, None) logger.debug( f"Removed sos{sos_type} constraint on {sos_dim} from {variable.name}" diff --git a/linopy/sos_reformulation.py b/linopy/sos_reformulation.py index de730708..beb83ab5 100644 --- a/linopy/sos_reformulation.py +++ b/linopy/sos_reformulation.py @@ -14,6 +14,8 @@ import pandas as pd from xarray import DataArray +from linopy.constants import SOS_BIG_M_ATTR, SOS_DIM_ATTR, SOS_TYPE_ATTR + if TYPE_CHECKING: from linopy.model import Model from linopy.variables import Variable @@ -50,7 +52,7 @@ def compute_big_m_values(var: Variable) -> DataArray: "SOS reformulation requires non-negative variables (lower >= 0)." ) - big_m_upper = var.attrs.get("big_m_upper") + big_m_upper = var.attrs.get(SOS_BIG_M_ATTR) M_upper = var.upper if big_m_upper is not None: @@ -84,7 +86,7 @@ def reformulate_sos1(model: Model, var: Variable, prefix: str) -> None: prefix : str Prefix for naming auxiliary variables and constraints. """ - sos_dim = str(var.attrs["sos_dim"]) + sos_dim = str(var.attrs[SOS_DIM_ATTR]) name = var.name M = compute_big_m_values(var) @@ -115,7 +117,7 @@ def reformulate_sos2(model: Model, var: Variable, prefix: str) -> None: prefix : str Prefix for naming auxiliary variables and constraints. """ - sos_dim = str(var.attrs["sos_dim"]) + sos_dim = str(var.attrs[SOS_DIM_ATTR]) name = var.name n = var.sizes[sos_dim] @@ -179,8 +181,8 @@ def reformulate_all_sos(model: Model, prefix: str = "_sos_reform_") -> list[str] for var_name in list(model.variables.sos): var = model.variables[var_name] - sos_type = var.attrs.get("sos_type") - sos_dim = var.attrs.get("sos_dim") + sos_type = var.attrs.get(SOS_TYPE_ATTR) + sos_dim = var.attrs.get(SOS_DIM_ATTR) if sos_type is None or sos_dim is None: continue diff --git a/linopy/variables.py b/linopy/variables.py index e2570b5d..75349c4b 100644 --- a/linopy/variables.py +++ b/linopy/variables.py @@ -52,7 +52,7 @@ to_polars, ) from linopy.config import options -from linopy.constants import HELPER_DIMS, TERM_DIM +from linopy.constants import HELPER_DIMS, SOS_DIM_ATTR, SOS_TYPE_ATTR, TERM_DIM from linopy.solver_capabilities import SolverFeature, solver_supports from linopy.types import ( ConstantLike, @@ -195,10 +195,10 @@ def __init__( if "label_range" not in data.attrs: data.assign_attrs(label_range=(data.labels.min(), data.labels.max())) - if "sos_type" in data.attrs or "sos_dim" in data.attrs: - if (sos_type := data.attrs.get("sos_type")) not in (1, 2): + if SOS_TYPE_ATTR in data.attrs or SOS_DIM_ATTR in data.attrs: + if (sos_type := data.attrs.get(SOS_TYPE_ATTR)) not in (1, 2): raise ValueError(f"sos_type must be 1 or 2, got {sos_type}") - if (sos_dim := data.attrs.get("sos_dim")) not in data.dims: + if (sos_dim := data.attrs.get(SOS_DIM_ATTR)) not in data.dims: raise ValueError( f"sos_dim must name a variable dimension, got {sos_dim}" ) @@ -328,8 +328,8 @@ def __repr__(self) -> str: dim_names = self.coord_names dim_sizes = list(self.sizes.values()) masked_entries = (~self.mask).sum().values - sos_type = self.attrs.get("sos_type") - sos_dim = self.attrs.get("sos_dim") + sos_type = self.attrs.get(SOS_TYPE_ATTR) + sos_dim = self.attrs.get(SOS_DIM_ATTR) lines = [] if dims: @@ -1244,8 +1244,8 @@ def __repr__(self) -> str: if ds.coords else "" ) - if (sos_type := ds.attrs.get("sos_type")) in (1, 2) and ( - sos_dim := ds.attrs.get("sos_dim") + if (sos_type := ds.attrs.get(SOS_TYPE_ATTR)) in (1, 2) and ( + sos_dim := ds.attrs.get(SOS_DIM_ATTR) ): coords += f" - sos{sos_type} on {sos_dim}" r += f" * {name}{coords}\n" @@ -1401,8 +1401,8 @@ def sos(self) -> Variables: { name: self.data[name] for name in self - if self[name].attrs.get("sos_dim") - and self[name].attrs.get("sos_type") in (1, 2) + if self[name].attrs.get(SOS_DIM_ATTR) + and self[name].attrs.get(SOS_TYPE_ATTR) in (1, 2) }, self.model, ) From 69610eb495642dbb0bef7e81cbe64bd9b345d457 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Tue, 20 Jan 2026 21:27:38 +0100 Subject: [PATCH 14/14] Add release notes --- doc/release_notes.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/release_notes.rst b/doc/release_notes.rst index d9bc95aa..4df695c7 100644 --- a/doc/release_notes.rst +++ b/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 --------------