From ff4e6ab0116f04b701ee45f901d739d45baa2daa Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Sun, 22 Feb 2026 11:10:41 -0600 Subject: [PATCH] Add indicator constraints --- doc/release_notes.rst | 2 +- linopy/constants.py | 4 + linopy/io.py | 111 ++++++++++++++++ linopy/model.py | 139 +++++++++++++++++++ linopy/solver_capabilities.py | 3 + test/test_indicator_constraints.py | 205 +++++++++++++++++++++++++++++ 6 files changed, 463 insertions(+), 1 deletion(-) create mode 100644 test/test_indicator_constraints.py diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 979b2263..72e5a637 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -5,7 +5,7 @@ Upcoming Version ---------------- * Add SOS1 and SOS2 reformulations for solvers not supporting them. - +* Add indicator constraints for solvers that support them. Version 0.6.4 -------------- diff --git a/linopy/constants.py b/linopy/constants.py index 2e1ef47a..fa6b07d8 100644 --- a/linopy/constants.py +++ b/linopy/constants.py @@ -54,6 +54,10 @@ SOS_DIM_ATTR = "sos_dim" SOS_BIG_M_ATTR = "big_m_upper" +# Indicator constraint attribute keys +INDICATOR_BINARY_VAR_ATTR = "indicator_binary_var" +INDICATOR_BINARY_VAL_ATTR = "indicator_binary_val" + class ModelStatus(Enum): """ diff --git a/linopy/io.py b/linopy/io.py index 54090e87..5863736b 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -399,6 +399,66 @@ def sos_to_file( _format_and_write(df, columns, f) +def indicator_constraints_to_file( + m: Model, + f: BufferedWriter, + explicit_coordinate_names: bool = False, +) -> None: + """ + Write indicator constraints to the s.t. section of an LP file. + + Indicator constraints appear in the Subject To section with the format: + ``ic0: b = 1 -> +1.0 x <= 5.0`` + """ + if not m.indicator_constraints: + return + + # If no regular constraints were written, we need the s.t. header + if not len(m.constraints): + f.write(b"\n\ns.t.\n\n") + + print_variable_scalar, _ = get_printers_scalar( + m, explicit_coordinate_names=explicit_coordinate_names + ) + + for ic_name, ic_data in m.indicator_constraints.items(): + labels_flat = ic_data.labels.values.flatten() + binary_var_flat = ic_data.binary_var.values.flatten() + binary_val_flat = np.broadcast_to( + ic_data.binary_val.values, labels_flat.shape + ).flatten() + coeffs_flat = ic_data.coeffs.values.reshape(len(labels_flat), -1) + vars_flat = ic_data.vars.values.reshape(len(labels_flat), -1) + sign_flat = np.broadcast_to(ic_data.sign.values, labels_flat.shape).flatten() + rhs_flat = np.broadcast_to(ic_data.rhs.values, labels_flat.shape).flatten() + + for i in range(len(labels_flat)): + if labels_flat[i] == -1: + continue + + bvar_name = print_variable_scalar(int(binary_var_flat[i])) + bval = int(binary_val_flat[i]) + + # Build LHS string from coeffs and vars + terms = [] + for c, v in zip(coeffs_flat[i], vars_flat[i]): + if v == -1: + continue + var_name = print_variable_scalar(int(v)) + coeff = float(c) + if coeff >= 0: + terms.append(f"+{coeff} {var_name}") + else: + terms.append(f"{coeff} {var_name}") + + lhs_str = " ".join(terms) + sign_str = str(sign_flat[i]) + rhs_val = float(rhs_flat[i]) + + line = f"ic{labels_flat[i]}: {bvar_name} = {bval} -> {lhs_str} {sign_str} {rhs_val}\n" + f.write(line.encode()) + + def constraints_to_file( m: Model, f: BufferedWriter, @@ -487,6 +547,11 @@ def to_lp_file( slice_size=slice_size, explicit_coordinate_names=explicit_coordinate_names, ) + indicator_constraints_to_file( + m, + f=f, + explicit_coordinate_names=explicit_coordinate_names, + ) bounds_to_file( m, f=f, @@ -757,6 +822,46 @@ def add_sos(s: xr.DataArray, sos_type: int, sos_dim: str) -> None: for _, s in stacked.groupby("_sos_group"): add_sos(s.unstack("_sos_group"), sos_type, sos_dim) + if m.indicator_constraints: + import gurobipy + + sense_map = { + "<=": gurobipy.GRB.LESS_EQUAL, + ">=": gurobipy.GRB.GREATER_EQUAL, + "=": gurobipy.GRB.EQUAL, + } + + for ic_name, ic_data in m.indicator_constraints.items(): + labels_flat = ic_data.labels.values.flatten() + binary_var_flat = ic_data.binary_var.values.flatten() + binary_val_flat = np.broadcast_to( + ic_data.binary_val.values, labels_flat.shape + ).flatten() + coeffs_flat = ic_data.coeffs.values.reshape(len(labels_flat), -1) + vars_flat = ic_data.vars.values.reshape(len(labels_flat), -1) + sign_flat = np.broadcast_to( + ic_data.sign.values, labels_flat.shape + ).flatten() + rhs_flat = np.broadcast_to(ic_data.rhs.values, labels_flat.shape).flatten() + + x_list = x.tolist() + for i in range(len(labels_flat)): + if labels_flat[i] == -1: + continue + + lhs = gurobipy.LinExpr() + for c, v in zip(coeffs_flat[i], vars_flat[i]): + if v != -1: + lhs.add(x_list[v], float(c)) + + model.addGenConstrIndicator( + x_list[int(binary_var_flat[i])], + bool(binary_val_flat[i]), + lhs, + sense_map[str(sign_flat[i])], + float(rhs_flat[i]), + ) + model.update() return model @@ -784,6 +889,12 @@ def to_highspy(m: Model, explicit_coordinate_names: bool = False) -> Highs: "Use io_api='lp' instead." ) + if m.indicator_constraints: + raise NotImplementedError( + "Indicator constraints are not supported by the HiGHS direct API. " + "Use a solver that supports them (gurobi, cplex)." + ) + import highspy print_variable, print_constraint = get_printers_scalar( diff --git a/linopy/model.py b/linopy/model.py index e72b3efa..8f6b6803 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -126,6 +126,7 @@ class Model: # containers "_variables", "_constraints", + "_indicator_constraints", "_objective", "_parameters", "_solution", @@ -136,8 +137,10 @@ class Model: # TODO: move counters to Variables and Constraints class "_xCounter", "_cCounter", + "_icCounter", "_varnameCounter", "_connameCounter", + "_icnameCounter", "_blocks", # TODO: check if these should not be mutable "_chunk", @@ -185,6 +188,7 @@ def __init__( """ self._variables: Variables = Variables({}, model=self) self._constraints: Constraints = Constraints({}, model=self) + self._indicator_constraints: dict[str, Dataset] = {} self._objective: Objective = Objective(LinearExpression(None, self), self) self._parameters: Dataset = Dataset() @@ -192,8 +196,10 @@ def __init__( self._termination_condition: str = "" self._xCounter: int = 0 self._cCounter: int = 0 + self._icCounter: int = 0 self._varnameCounter: int = 0 self._connameCounter: int = 0 + self._icnameCounter: int = 0 self._blocks: DataArray | None = None self._chunk: T_Chunks = chunk @@ -219,6 +225,16 @@ def constraints(self) -> Constraints: """ return self._constraints + @property + def indicator_constraints(self) -> dict[str, Dataset]: + """ + Indicator constraints assigned to the model. + + Each entry is a Dataset with fields: coeffs, vars, sign, rhs, + labels, binary_var, binary_val. + """ + return self._indicator_constraints + @property def objective(self) -> Objective: """ @@ -814,6 +830,122 @@ def add_constraints( self.constraints.add(constraint) return constraint + def add_indicator_constraints( + self, + binary_var: Variable, + binary_val: int, + lhs: ConstraintLike | ExpressionLike | VariableLike, + sign: SignLike | None = None, + rhs: ConstantLike | None = None, + name: str | None = None, + ) -> Dataset: + """ + Add indicator constraints to the model. + + An indicator constraint has the form: + (binary_var == binary_val) => (linear_constraint) + + The linear constraint is only enforced when binary_var equals + binary_val. These constraints are handled natively by solvers + like Gurobi and CPLEX via general constraints. + + Parameters + ---------- + binary_var : linopy.Variable + Binary variable serving as the indicator. Must have binary=True. + binary_val : int + Triggering value, must be 0 or 1. + lhs : linopy.Constraint, linopy.LinearExpression, or linopy.Variable + The conditionally enforced constraint. If a LinearExpression or + Variable is passed, ``sign`` and ``rhs`` must also be provided. + sign : str, optional + Constraint sign ('<=', '>=', '='). Required when ``lhs`` is an + expression. + rhs : numeric, optional + Right-hand side. Required when ``lhs`` is an expression. + name : str, optional + Name for the indicator constraint group. + + Returns + ------- + xarray.Dataset + The stored indicator constraint data. + """ + if not binary_var.attrs.get("binary", False): + raise ValueError( + "Indicator variable must be binary. " + f"Variable '{binary_var.name}' is not binary." + ) + + if binary_val not in (0, 1): + raise ValueError(f"binary_val must be 0 or 1, got {binary_val}.") + + if name in self._indicator_constraints: + raise ValueError( + f"Indicator constraint '{name}' already assigned to model." + ) + if name is None: + name = f"indcon{self._icnameCounter}" + self._icnameCounter += 1 + + # Build constraint data from lhs + if isinstance(lhs, Constraint): + if sign is not None or rhs is not None: + raise ValueError("sign and rhs must be None when lhs is a Constraint.") + data = lhs.data + elif isinstance(lhs, LinearExpression): + if sign is None or rhs is None: + raise ValueError( + "sign and rhs are required when lhs is a LinearExpression." + ) + sign = maybe_replace_signs(as_dataarray(sign)) + data = lhs.to_constraint(sign, rhs).data + elif isinstance(lhs, Variable | ScalarVariable | ScalarLinearExpression): + if sign is None or rhs is None: + raise ValueError("sign and rhs are required when lhs is a Variable.") + sign = maybe_replace_signs(as_dataarray(sign)) + data = lhs.to_linexpr().to_constraint(sign, rhs).data + elif isinstance(lhs, AnonymousScalarConstraint): + if sign is not None or rhs is not None: + raise ValueError("sign and rhs must be None when lhs is a Constraint.") + data = lhs.to_constraint().data + else: + raise TypeError( + f"lhs must be a Constraint, LinearExpression, or Variable, " + f"got {type(lhs)}." + ) + + # Add binary variable labels and value + binary_var_labels = binary_var.labels + data["binary_var"] = binary_var_labels + data["binary_val"] = binary_val + + # Broadcast all fields together + (data,) = xr.broadcast(data, exclude=[TERM_DIM]) + + # Assign unique labels with the same shape as rhs (non-term dims) + labels = DataArray( + np.full(data.rhs.shape, -1, dtype=int), + coords=data.rhs.coords, + dims=data.rhs.dims, + ) + data["labels"] = labels + start = self._icCounter + end = start + data.labels.size + data.labels.values = np.arange(start, end).reshape(data.labels.shape) + self._icCounter += data.labels.size + + data = data.assign_attrs(label_range=(start, end), name=name) + + self._indicator_constraints[name] = data + return data + + def remove_indicator_constraints(self, name: str) -> None: + """ + Remove indicator constraint by name. + """ + del self._indicator_constraints[name] + def add_objective( self, expr: Variable @@ -1408,6 +1540,13 @@ def solve( "Use reformulate_sos=True or a solver that supports SOS (gurobi, cplex)." ) + if self.indicator_constraints: + if not solver_supports(solver_name, SolverFeature.INDICATOR_CONSTRAINTS): + raise ValueError( + f"Solver {solver_name} does not support indicator constraints. " + "Use a solver that supports them (gurobi, cplex)." + ) + try: solver_class = getattr(solvers, f"{solvers.SolverName(solver_name).name}") # initialize the solver as object of solver subclass diff --git a/linopy/solver_capabilities.py b/linopy/solver_capabilities.py index f0507317..9e543861 100644 --- a/linopy/solver_capabilities.py +++ b/linopy/solver_capabilities.py @@ -49,6 +49,7 @@ class SolverFeature(Enum): # Special constraint types SOS_CONSTRAINTS = auto() # Special Ordered Sets (SOS1/SOS2) constraints + INDICATOR_CONSTRAINTS = auto() # Indicator (conditional) constraints # Solver-specific SOLVER_ATTRIBUTE_ACCESS = auto() # Direct access to solver variable attributes @@ -86,6 +87,7 @@ def supports(self, feature: SolverFeature) -> bool: SolverFeature.SOLUTION_FILE_NOT_NEEDED, SolverFeature.IIS_COMPUTATION, SolverFeature.SOS_CONSTRAINTS, + SolverFeature.INDICATOR_CONSTRAINTS, SolverFeature.SOLVER_ATTRIBUTE_ACCESS, } ), @@ -134,6 +136,7 @@ def supports(self, feature: SolverFeature) -> bool: SolverFeature.LP_FILE_NAMES, SolverFeature.READ_MODEL_FROM_FILE, SolverFeature.SOS_CONSTRAINTS, + SolverFeature.INDICATOR_CONSTRAINTS, } ), ), diff --git a/test/test_indicator_constraints.py b/test/test_indicator_constraints.py new file mode 100644 index 00000000..2b404168 --- /dev/null +++ b/test/test_indicator_constraints.py @@ -0,0 +1,205 @@ +"""Tests for indicator constraint support.""" + +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from linopy import Model, available_solvers + + +def test_add_indicator_constraints_basic() -> None: + """Indicator constraint is created with correct fields.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + assert "ic0" in m.indicator_constraints + ic = m.indicator_constraints["ic0"] + assert "coeffs" in ic + assert "vars" in ic + assert "sign" in ic + assert "rhs" in ic + assert "binary_var" in ic + assert "binary_val" in ic + assert "labels" in ic + + +def test_add_indicator_constraints_from_constraint() -> None: + """Indicator constraint accepts a Constraint object as lhs.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + con = x <= 5 + m.add_indicator_constraints(b, 1, con, name="ic0") + assert "ic0" in m.indicator_constraints + + +def test_indicator_constraints_validation_non_binary() -> None: + """Non-binary variable raises ValueError.""" + m = Model() + y = m.add_variables(lower=0, upper=1, name="y") + x = m.add_variables(lower=0, upper=10, name="x") + with pytest.raises(ValueError, match="must be binary"): + m.add_indicator_constraints(y, 1, x, "<=", 5) + + +def test_indicator_constraints_validation_bad_value() -> None: + """Invalid binary_val raises ValueError.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + with pytest.raises(ValueError, match="must be 0 or 1"): + m.add_indicator_constraints(b, 2, x, "<=", 5) + + +def test_indicator_constraints_duplicate_name() -> None: + """Duplicate name raises ValueError.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + with pytest.raises(ValueError, match="already assigned"): + m.add_indicator_constraints(b, 0, x, ">=", 0, name="ic0") + + +def test_indicator_constraints_not_in_regular_constraints() -> None: + """Indicator constraints are separate from regular constraints.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_constraints(x >= 0, name="regular") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + assert "regular" in m.constraints + assert "ic0" not in m.constraints + assert "ic0" in m.indicator_constraints + + +def test_remove_indicator_constraints() -> None: + """Indicator constraints can be removed.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + assert "ic0" in m.indicator_constraints + m.remove_indicator_constraints("ic0") + assert "ic0" not in m.indicator_constraints + + +def test_indicator_constraints_lp_file(tmp_path: Path) -> None: + """LP file contains general constraints section.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_constraints(x >= 0, name="dummy") + m.add_objective(x) + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + fn = tmp_path / "test.lp" + m.to_file(fn) + content = fn.read_text() + assert "= 1 ->" in content + assert "ic0:" in content + + +@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") +def test_indicator_constraints_solve_gurobi_active() -> None: + """ + Indicator constraint enforced when binary is at trigger value. + + Maximize x subject to: + b = 1 (fixed) + (b == 1) => (x <= 5) + x in [0, 10] + Since b=1 and the indicator enforces x<=5, optimal x=5. + """ + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + # Force b = 1 + m.add_constraints(b >= 1, name="fix_b") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + m.add_objective(x, sense="max") + m.solve(solver_name="gurobi") + assert m.objective.value is not None + assert np.isclose(m.objective.value, 5, atol=1e-6) + + +@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") +def test_indicator_constraints_solve_gurobi_inactive() -> None: + """ + Indicator constraint NOT enforced when binary differs from trigger. + + Maximize x subject to: + b = 1 (fixed) + (b == 0) => (x <= 5) + x in [0, 10] + Since b=1 but trigger is 0, the constraint is inactive. Optimal x=10. + """ + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + # Force b = 1 + m.add_constraints(b >= 1, name="fix_b") + m.add_indicator_constraints(b, 0, x, "<=", 5, name="ic0") + m.add_objective(x, sense="max") + m.solve(solver_name="gurobi") + assert m.objective.value is not None + assert np.isclose(m.objective.value, 10, atol=1e-6) + + +@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") +def test_indicator_constraints_solve_gurobi_multiple() -> None: + """ + Multiple indicator constraints work together. + + b1=1, b2=1 (forced), x in [0, 20]. + (b1 == 1) => (x <= 10) + (b2 == 1) => (x <= 5) + Maximize x. + Both constraints active, so x <= min(10, 5) = 5. Optimal x=5. + """ + m = Model() + b1 = m.add_variables(name="b1", binary=True) + b2 = m.add_variables(name="b2", binary=True) + x = m.add_variables(lower=0, upper=20, name="x") + m.add_constraints(b1 >= 1, name="fix_b1") + m.add_constraints(b2 >= 1, name="fix_b2") + m.add_indicator_constraints(b1, 1, x, "<=", 10, name="ic1") + m.add_indicator_constraints(b2, 1, x, "<=", 5, name="ic2") + m.add_objective(x, sense="max") + m.solve(solver_name="gurobi") + assert m.objective.value is not None + assert np.isclose(m.objective.value, 5, atol=1e-6) + + +def test_unsupported_solver_raises() -> None: + """Solvers without indicator support raise ValueError.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_constraints(x >= 0, name="dummy") + m.add_objective(x) + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + for solver in ["glpk", "highs", "mosek", "mindopt"]: + if solver in available_solvers: + with pytest.raises( + ValueError, match="does not support indicator constraints" + ): + m.solve(solver_name=solver) + + +def test_indicator_constraints_with_coords() -> None: + """Indicator constraints work with multi-dimensional coords.""" + m = Model() + idx = pd.RangeIndex(3, name="i") + b = m.add_variables(coords=[idx], name="b", binary=True) + x = m.add_variables(coords=[idx], lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + ic = m.indicator_constraints["ic0"] + # Should have 3 indicator constraints (one per index) + assert ic.labels.size == 3