feat: Add indicator constraints

doc/release_notes.rst

@@ -37,6 +37,11 @@ Most users should keep calling ``model.solve(...)``. If you want more control, y
 * New ``linopy.licensed_solvers``: the subset of installed solvers that currently pass a license check. Handy in tests and for picking a solver at runtime.
 * New helpers for explicit license checks: ``linopy.solvers.check_solver_licenses("gurobi", "mosek")``, ``Gurobi.license_status()``, ``Gurobi.is_available()``. They return a ``LicenseStatus`` dataclass (``name``, ``ok``, ``message``).
 
+*Constraints — indicator constraints*
+
+* Add indicator constraints for solvers that support them. They are part of the unified constraints container: ``model.add_indicator_constraints`` returns a ``Constraint`` and the constraint is stored in ``model.constraints`` (filterable via ``model.constraints.indicator`` / ``model.constraints.regular``), so it round-trips through netCDF and ``model.copy()``.
+
+
 *Constraints — CSR-backed storage*
 
 * Add ``CSRConstraint``: a memory-efficient immutable constraint representation backed by scipy CSR sparse matrices. Up to 90% memory savings for constraints with many terms and 30–120× faster matrix generation for direct solver APIs.

linopy/common.py

@@ -786,7 +786,11 @@ def __init__(self, constraints: Any) -> None:
     @cached_property
     def clabels(self) -> np.ndarray:
         """Active constraint labels in build order, shape (n_active_cons,)."""
-        label_lists = [c.active_labels() for c in self._constraints.data.values()]
+        label_lists = [
+            c.active_labels()
+            for c in self._constraints.data.values()
+            if not c.is_indicator
+        ]
         return (
             np.concatenate(label_lists) if label_lists else np.array([], dtype=np.intp)
         )

linopy/constants.py

@@ -92,6 +92,10 @@ class PerformanceWarning(UserWarning):
 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 EvolvingAPIWarning(FutureWarning):
     """

linopy/constraints.py

@@ -179,6 +179,21 @@ def dual(self) -> DataArray:
     def dual(self, value: DataArray) -> None:
         """Set the dual values DataArray."""
 
+    @property
+    @abstractmethod
+    def is_indicator(self) -> bool:
+        """Whether the constraint is an indicator constraint."""
+
+    @property
+    @abstractmethod
+    def binary_var(self) -> DataArray | None:
+        """Get the indicator binary variable labels, or None."""
+
+    @property
+    @abstractmethod
+    def binary_val(self) -> int | np.ndarray | None:
+        """Get the indicator triggering value(s), or None."""
+
     @abstractmethod
     def has_variable(self, variable: variables.Variable) -> bool:
         """Check if the constraint references any of the given variable labels."""
@@ -509,6 +524,8 @@ class CSRConstraint(ConstraintBase):
         "_name",
         "_cindex",
         "_dual",
+        "_binvar_labels",
+        "_binval",
     )
 
     def __init__(
@@ -522,6 +539,8 @@ def __init__(
         name: str = "",
         cindex: int | None = None,
         dual: np.ndarray | None = None,
+        binvar_labels: np.ndarray | None = None,
+        binval: int | np.ndarray | None = None,
     ) -> None:
         self._csr = csr
         self._con_labels = con_labels
@@ -532,6 +551,8 @@ def __init__(
         self._name = name
         self._cindex = cindex
         self._dual = dual
+        self._binvar_labels = binvar_labels
+        self._binval = binval
 
     @property
     def model(self) -> Model:
@@ -686,6 +707,20 @@ def dual(self, value: DataArray) -> None:
         vals = np.asarray(value).ravel()
         self._dual = vals[self.active_positions]
 
+    @property
+    def is_indicator(self) -> bool:
+        return self._binvar_labels is not None
+
+    @property
+    def binary_var(self) -> DataArray | None:
+        if self._binvar_labels is None:
+            return None
+        return self._active_to_dataarray(self._binvar_labels, fill=-1)
+
+    @property
+    def binary_val(self) -> int | np.ndarray | None:
+        return self._binval
+
     def _to_dataset(self, nterm: int) -> Dataset:
         """
         Reconstruct labels/coeffs/vars Dataset from the CSR matrix.
@@ -743,9 +778,18 @@ def _to_dataset(self, nterm: int) -> Dataset:
     def data(self) -> Dataset:
         """Reconstruct the xarray Dataset from the CSR representation."""
         ds = self._to_dataset(self.nterm)
-        extra = {"sign": self.sign, "rhs": self.rhs}
+        extra: dict[str, Any] = {"sign": self.sign, "rhs": self.rhs}
         if self._dual is not None:
             extra["dual"] = self._active_to_dataarray(self._dual, fill=np.nan)
+        if self._binvar_labels is not None:
+            extra["binary_var"] = self._active_to_dataarray(
+                self._binvar_labels, fill=-1
+            )
+            binval = self._binval
+            if isinstance(binval, np.ndarray) and binval.ndim > 0:
+                extra["binary_val"] = self._active_to_dataarray(binval, fill=-1)
+            else:
+                extra["binary_val"] = binval
         return assign_multiindex_safe(ds, **extra).assign_attrs(self.attrs)
 
     def __repr__(self) -> str:
@@ -833,6 +877,14 @@ def to_netcdf_ds(self) -> Dataset:
         }
         if isinstance(self._sign, str):
             attrs["sign"] = self._sign
+        if self._binvar_labels is not None:
+            attrs["is_indicator"] = True
+            data_vars["_binvar_labels"] = DataArray(self._binvar_labels, dims=["_flat"])
+            binval = self._binval
+            if isinstance(binval, np.ndarray) and binval.ndim > 0:
+                data_vars["_binval"] = DataArray(binval, dims=["_flat"])
+            else:
+                attrs["binval"] = int(binval)  # type: ignore[arg-type]
         return Dataset(data_vars, attrs=attrs)
 
     @classmethod
@@ -859,8 +911,23 @@ def from_netcdf_ds(cls, ds: Dataset, model: Model, name: str) -> CSRConstraint:
             con_labels = np.arange(cindex, cindex + len(rhs), dtype=np.intp)
         else:
             con_labels = np.arange(len(rhs), dtype=np.intp)
+        binvar_labels: np.ndarray | None = None
+        binval: int | np.ndarray | None = None
+        if "_binvar_labels" in ds:
+            binvar_labels = ds["_binvar_labels"].values
+            binval = ds["_binval"].values if "_binval" in ds else attrs["binval"]
         return cls(
-            csr, con_labels, rhs, sign, coords, model, name, cindex=cindex, dual=dual
+            csr,
+            con_labels,
+            rhs,
+            sign,
+            coords,
+            model,
+            name,
+            cindex=cindex,
+            dual=dual,
+            binvar_labels=binvar_labels,
+            binval=binval,
         )
 
     def has_variable(self, variable: variables.Variable) -> bool:
@@ -1023,6 +1090,16 @@ def from_mutable(
         dual = (
             con.data["dual"].values.ravel()[active_mask] if "dual" in con.data else None
         )
+        binvar_labels: np.ndarray | None = None
+        binval: int | np.ndarray | None = None
+        binary_var = con.binary_var
+        if binary_var is not None:
+            binvar_labels = binary_var.values.ravel()[active_mask]
+            bv = con.binary_val
+            if isinstance(bv, np.ndarray) and bv.ndim > 0:
+                binval = bv.ravel()[active_mask]
+            else:
+                binval = bv
         return cls(
             csr,
             con_labels,
@@ -1033,6 +1110,8 @@ def from_mutable(
             con.name,
             cindex=cindex,
             dual=dual,
+            binvar_labels=binvar_labels,
+            binval=binval,
         )
 
 
@@ -1157,6 +1236,20 @@ def rhs(self, value: ExpressionLike) -> None:
         self.lhs = self.lhs - value.reset_const()
         self._data = assign_multiindex_safe(self.data, rhs=value.const)
 
+    @property
+    def is_indicator(self) -> bool:
+        return "binary_var" in self._data
+
+    @property
+    def binary_var(self) -> DataArray | None:
+        return self._data.get("binary_var")
+
+    @property
+    def binary_val(self) -> int | np.ndarray | None:
+        if "binary_val" in self._data:
+            return self._data["binary_val"].values
+        return None
+
     @property
     def lhs(self) -> expressions.LinearExpression:
         data = self.data[["coeffs", "vars"]].rename({self.term_dim: TERM_DIM})
@@ -1661,7 +1754,7 @@ def ncons(self) -> int:
         This excludes constraints with missing labels or where all variables
         are masked (vars == -1).
         """
-        return sum(con.ncons for con in self.data.values())
+        return sum(con.ncons for con in self.data.values() if not con.is_indicator)
 
     @property
     def inequalities(self) -> Constraints:
@@ -1677,6 +1770,20 @@ def equalities(self) -> Constraints:
         """
         return self[[n for n, s in self.items() if (s.sign == EQUAL).all()]]
 
+    @property
+    def indicator(self) -> Constraints:
+        """
+        Get the subset of constraints which are indicator constraints.
+        """
+        return self[[n for n, c in self.items() if c.is_indicator]]
+
+    @property
+    def regular(self) -> Constraints:
+        """
+        Get the subset of constraints which are not indicator constraints.
+        """
+        return self[[n for n, c in self.items() if not c.is_indicator]]
+
     def sanitize_zeros(self) -> None:
         """
         Filter out terms with zero and close-to-zero coefficient.
@@ -1865,9 +1972,13 @@ def to_matrix(self) -> tuple[scipy.sparse.csr_array, np.ndarray]:
         csrs = []
         con_labels_list = []
         for c in self.data.values():
+            if c.is_indicator:
+                continue
             csr, con_labels = c.to_matrix(label_index)
             csrs.append(csr)
             con_labels_list.append(con_labels)
+        if not csrs:
+            raise ValueError("No constraints available to convert to matrix.")
         return scipy.sparse.vstack(csrs, format="csr"), np.concatenate(con_labels_list)
 
     def reset_dual(self) -> None:

linopy/io.py

@@ -459,6 +459,61 @@ 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: x0 = 1 -> +1.0 x1 <= 5.0``
+    """
+    if not len(m.constraints.indicator):
+        return
+
+    if not len(m.constraints.regular):
+        f.write(b"\n\ns.t.\n\n")
+
+    print_variable_scalar, _ = get_printers_scalar(
+        m, explicit_coordinate_names=explicit_coordinate_names
+    )
+
+    for con in m.constraints.indicator.data.values():
+        ic_data = con.data
+        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(binary_var_flat[i : i + 1])[0]
+            valid = vars_flat[i] != -1
+            var_names = print_variable_scalar(vars_flat[i][valid])
+
+            terms = []
+            for coeff, var_name in zip(coeffs_flat[i][valid], var_names):
+                coeff = float(coeff)
+                prefix = "+" if coeff >= 0 else ""
+                terms.append(f"{prefix}{coeff} {var_name}")
+
+            lhs_str = " ".join(terms)
+            line = (
+                f"ic{labels_flat[i]}: {bvar_name} = {int(binary_val_flat[i])} -> "
+                f"{lhs_str} {sign_flat[i]} {float(rhs_flat[i])}\n"
+            )
+            f.write(line.encode())
+
+
 def constraints_to_file(
     m: Model,
     f: BufferedWriter,
@@ -467,26 +522,27 @@ def constraints_to_file(
     slice_size: int = 2_000_000,
     explicit_coordinate_names: bool = False,
 ) -> None:
-    if not len(m.constraints):
+    regular = m.constraints.regular
+    if not len(regular):
         return
 
     print_variable, print_constraint = get_printers(
         m, explicit_coordinate_names=explicit_coordinate_names
     )
 
     f.write(b"\n\ns.t.\n\n")
-    names = m.constraints
+    names = list(regular)
     if progress:
         names = tqdm(
-            list(names),
+            names,
             desc="Writing constraints.",
             colour=TQDM_COLOR,
         )
 
     # to make this even faster, we can use polars expression
     # https://docs.pola.rs/user-guide/expressions/plugins/#output-data-types
     for name in names:
-        con = m.constraints[name]
+        con = regular[name]
         for con_slice in con.iterate_slices(slice_size):
             df = con_slice.to_polars()
 
@@ -547,6 +603,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,
@@ -1092,7 +1153,7 @@ def copy(m: Model, include_solution: bool = False, deep: bool = True) -> Model:
     Model
         A deep or shallow copy of the model.
     """
-    from linopy.constraints import Constraint, Constraints
+    from linopy.constraints import Constraint, ConstraintBase, Constraints
     from linopy.expressions import LinearExpression
     from linopy.model import Model, Objective
     from linopy.variables import Variable, Variables
@@ -1122,15 +1183,18 @@ def copy(m: Model, include_solution: bool = False, deep: bool = True) -> Model:
         new_model,
     )
 
+    def _copy_con_data(con: ConstraintBase) -> xr.Dataset:
+        d = con.mutable().data
+        if include_solution:
+            return d.copy(deep=deep)
+        cols = m.constraints.dataset_attrs + (
+            ["binary_var", "binary_val"] if con.is_indicator else []
+        )
+        return d[cols].copy(deep=deep)
+
     new_model._constraints = Constraints(
         {
-            name: Constraint(
-                con.mutable().data.copy(deep=deep)
-                if include_solution
-                else con.mutable().data[m.constraints.dataset_attrs].copy(deep=deep),
-                new_model,
-                name,
-            )
+            name: Constraint(_copy_con_data(con), new_model, name)
             for name, con in m.constraints.items()
         },
         new_model,