feat(constraint.py): Auto-convert hard to soft constraints through soften method

doc/release_notes.rst

@@ -29,6 +29,11 @@ Upcoming Version
 
 * Added support for the GPU-accelerated `NVIDIA cuOpt <https://docs.nvidia.com/cuopt/>`__ solver for linear, mixed-integer and convex quadratic problems, via ``model.solve("cuopt", io_api="direct")``. Install it with ``pip install "linopy[gpu]"`` — Linux only, and requires an NVIDIA GPU of compute capability 7.0 or higher with a CUDA 12 driver (525.60.13 or newer). See :doc:`gpu-acceleration` for the supported problem classes and the known limitations.
 
+*New feature: constraint softening*
+
+* A constraint can now be softened with ``Constraint.soften(penalty, max_violation=None, name=None)``, which adds a slack variable (a positive/negative pair for equality constraints) to the constraint's ``lhs`` and a penalty term to the objective, returning a ``Slack`` named tuple. It is not supported on frozen constraints or on detached copies from ``.mutable()``, ``.sel()``, or ``.isel()``. ``model.add_constraints(..., penalty=...)`` is a shortcut that softens the constraint right after creation and cannot be combined with ``freeze=True``.
+* The slack variable(s) created by ``soften()`` can be retrieved afterwards via the new ``Constraint.slack`` property.
+
 *Other*
 
 * ``add_piecewise_formulation`` gained a ``mask`` parameter declaring which breakpoint slots hold a real breakpoint. It is needed for **ragged** curves — entities with different numbers of breakpoints — which are stored densely with the surplus slots left absent. Under v1 that absence must be declared (``mask=x_pts.notnull()``) rather than read off the NaN padding. (https://github.com/PyPSA/linopy/issues/884)

linopy/constraints.py

@@ -16,6 +16,7 @@
 from typing import (
     TYPE_CHECKING,
     Any,
+    NamedTuple,
     overload,
 )
 from warnings import warn
@@ -1431,6 +1432,18 @@ def _rhs_grid_values(expr: CSRLinearExpression, rhs: DataArray) -> np.ndarray:
     return rhs.transpose(*expr.grid.dims).to_numpy().reshape(-1)
 
 
+class Slack(NamedTuple):
+    """
+    Slack variable(s) added by :meth:`Constraint.soften`.
+
+    `negative` is None for inequality constraints, since those only
+    need one slack variable to absorb a violation in a single direction.
+    """
+
+    positive: variables.Variable
+    negative: variables.Variable | None
+
+
 class Constraint(ConstraintBase):
     """
     Constraint backed by an xarray Dataset.
@@ -1615,6 +1628,21 @@ def lhs(self, value: ExpressionLike | VariableLike | ConstantLike) -> None:
         )
         self.update(lhs=value)
 
+    @property
+    def slack(self) -> Slack | None:
+        """
+        Slack variable(s) added via :meth:`soften`, or ``None`` if the
+        constraint has never been softened.
+        """
+        positive = self.data.attrs.get("slack_positive")
+        if positive is None:
+            return None
+        negative = self.data.attrs.get("slack_negative", "")
+        return Slack(
+            positive=self.model.variables[positive],
+            negative=self.model.variables[negative] if negative else None,
+        )
+
     def _assign_lhs(
         self, expr: expressions.LinearExpression, rhs: DataArray | None = None
     ) -> None:
@@ -2015,6 +2043,135 @@ def from_rule(cls, model: Model, rule: Callable, coords: CoordsLike) -> Constrai
         data = lhs.data.assign(sign=sign, rhs=rhs)
         return cls(data, model=model)
 
+    def soften(
+        self,
+        penalty: ConstantLike,
+        *,
+        max_violation: ConstantLike | None = None,
+        name: str | None = None,
+    ) -> Slack:
+        """
+        Soften a constraint, adding a slack variable and a penalty to the objective function.
+
+        Parameters
+        ----------
+        penalty : constant-like
+            The penalty that will match the slack variable inside the objective function. Must be bigger than 0.
+        max_violation: constant-like
+            The max violation possible that caps the slack (upper bound). If None, the slack will be unbounded.
+        name: string
+            The name for the slack variable. If None, it well reuse the constraint name and add a '_slack'.
+
+        Returns
+        -------
+        Slack
+            Named tuple with the `positive` slack variable, and the `negative` one for equality constraints (`None`
+            for inequality constraints).
+
+        Notes
+        -----
+        Not supported on frozen constraints (e.g. a CSRConstraint from add_constraints(..., freeze=True) or
+        Model(freeze_constraints=True)). This method is only defined on Constraint and calling it on a CSRConstraint
+        raises AttributeError. Calling .mutable() first does not help either, since the resulting Constraint is a
+        detached copy not registered in model.constraints, so soften raises ValueError on it instead.
+
+        Softening an already-softened constraint raises ValueError instead of stacking a second, redundant slack term
+        onto the same lhs.
+
+        Examples
+        --------
+        >>> from linopy import Model
+        >>> import pandas as pd
+
+        >>> m = Model()
+        >>> investments = pd.Index(["A", "B", "C"], name="investments")
+        >>> expected_return = pd.Series(
+        ...     [0.08, 0.03, 0.1], index=investments, name="expected_return"
+        ... )
+        >>> w = m.add_variables(lower=0, upper=1, coords=[investments], name="weights")
+        >>> m.add_objective((expected_return * w).sum(), sense="max")
+        >>> budget_penalty = 2
+
+        >>> budget_constraint = m.add_constraints(w.sum() == 1, name="budget")
+        >>> slack = budget_constraint.soften(penalty=budget_penalty)
+        """
+        # Verify valid penalty to continue:
+        if not bool(np.all(np.asarray(penalty) > 0)):
+            raise ValueError("Penalty is not positive.")
+
+        # Require the objective function to exist before using soften method (this is to avoid
+        # `add_objective` overwriting the penalty term added below, since it replaces rather than merges):
+        model = self.model
+        if model.objective.expression.empty:
+            raise ValueError(
+                "Objective must be defined via `model.add_objective` before calling `soften` on constraints."
+            )
+
+        # A detached copy of the constraint (e.g. from `.mutable()`, `.sel()`, `.isel()`) isn't in
+        # `model.constraints`, so .soften would silently do nothing on the real model. This check is to avoid that:
+        if model.constraints.data.get(self.name) is not self:
+            raise ValueError(
+                f"Constraint {self.name!r} is not the constraint registered in the model, so "
+                "`soften` would not affect it (it may be a detached copy from `.mutable()`, "
+                "`.sel()`, or `.isel()`). Call `soften` on `model.constraints[name]` directly."
+            )
+
+        if self.slack is not None:
+            raise ValueError(
+                f"Constraint {self.name!r} was already softened (existing slack "
+                f"variable {self.slack.positive.name!r})"
+            )
+
+        name = name or f"{self.name}_slack"
+        upper = np.inf if max_violation is None else max_violation
+
+        sign_values = pd.unique(self.sign.values.ravel())
+        if len(sign_values) > 1:
+            raise NotImplementedError(
+                "Constraint.soften does not support constraints with mixed signs."
+            )
+
+        positive_slack = model.add_variables(
+            lower=0,
+            upper=upper,
+            coords=self.lhs.coords,
+            mask=self.mask,
+            name=f"{name}_pos",
+        )
+        negative_slack = None
+
+        # Update left hand side depending on the sign of the constraint:
+        sign = sign_values.item()
+        if sign == "<=":
+            self.update(lhs=self.lhs - positive_slack)
+        elif sign == ">=":
+            self.update(lhs=self.lhs + positive_slack)
+        else:
+            negative_slack = model.add_variables(
+                lower=0,
+                upper=upper,
+                coords=self.lhs.coords,
+                mask=self.mask,
+                name=f"{name}_neg",
+            )
+            self.update(lhs=self.lhs - positive_slack + negative_slack)
+
+        # Update objective function:
+        constraint_violation = (
+            positive_slack + negative_slack
+            if negative_slack is not None
+            else positive_slack
+        )
+        direction = 1 if model.sense == "min" else -1
+        model.objective += direction * (penalty * constraint_violation).sum()
+
+        self._data = self._data.assign_attrs(
+            slack_positive=positive_slack.name,
+            slack_negative=negative_slack.name if negative_slack is not None else "",
+        )
+
+        return Slack(positive=positive_slack, negative=negative_slack)
+
     def to_polars(self) -> pl.DataFrame:
         """
         Convert the constraint to a polars DataFrame.

linopy/model.py

@@ -1190,6 +1190,7 @@ def add_constraints(
         mask: MaskLike | None = ...,
         freeze: Literal[False] = ...,
         scaling: ConstantLike = ...,
+        penalty: ConstantLike | None = ...,
     ) -> Constraint: ...
 
     @overload
@@ -1207,6 +1208,7 @@ def add_constraints(
         mask: MaskLike | None = ...,
         freeze: Literal[True] = ...,
         scaling: ConstantLike = ...,
+        penalty: None = ...,
     ) -> CSRConstraint: ...
 
     @overload
@@ -1224,6 +1226,7 @@ def add_constraints(
         mask: MaskLike | None = ...,
         freeze: bool | None = ...,
         scaling: ConstantLike = ...,
+        penalty: ConstantLike | None = ...,
     ) -> ConstraintBase: ...
 
     def add_constraints(
@@ -1240,6 +1243,7 @@ def add_constraints(
         mask: MaskLike | None = None,
         freeze: bool | None = None,
         scaling: ConstantLike = 1,
+        penalty: ConstantLike | None = None,
     ) -> ConstraintBase:
         """
         Assign a new, possibly multi-dimensional array of constraints to the
@@ -1282,6 +1286,16 @@ def add_constraints(
             Positive finite scaling factor(s) for constraint rows. Solver-side
             left-hand-side coefficients and right-hand-side values are multiplied
             by this factor. The default is 1.
+        penalty : constant-like, optional
+            If given, soften the constraint right away by calling
+            :meth:`Constraint.soften` with this penalty, adding a slack variable
+            and a penalty term to the objective. Not allowed together with
+            ``freeze=True`` (or a model default of ``freeze_constraints=True``),
+            since softening requires a mutable, registered ``Constraint``.
+            The resulting Slack is not returned by this shortcut; retrieve
+            the slack variable(s) from model.variables using the derived
+            name f"{name}_slack_pos" (and f"{name}_slack_neg" for
+            equality constraints).
 
         Returns
         -------
@@ -1294,6 +1308,12 @@ def add_constraints(
             freeze = self.freeze_constraints
         freeze = freeze and not self.chunk
 
+        if penalty is not None and freeze:
+            raise ValueError(
+                "`penalty` cannot be combined with `freeze=True` (or a model default of `freeze_constraints=True`),"
+                "since `soften` is not supported on frozen constraints."
+            )
+
         if isinstance(sign, str):
             sign = maybe_replace_sign(sign)
         elif sign is not None:
@@ -1383,7 +1403,10 @@ def add_constraints(
 
         enforce_no_multiindex(data, context=f"constraint {name!r}")
         constraint = Constraint(data, name=name, model=self, skip_broadcast=True)
-        return self.constraints.add(constraint, freeze=freeze)
+        added = self.constraints.add(constraint, freeze=freeze)
+        if penalty is not None:
+            constraint.soften(penalty=penalty)
+        return added
 
     def add_indicator_constraints(
         self,

linopy/objective.py

@@ -273,7 +273,7 @@ def to_matrix(self, *args: Any, **kwargs: Any) -> csc_matrix:
     sel = objwrap(expressions.LinearExpression.sel)
 
     def __add__(
-        self, expr: int | QuadraticExpression | LinearExpression | Objective
+        self, expr: ConstantLike | QuadraticExpression | LinearExpression | Objective
     ) -> Objective:
         if isinstance(expr, Objective):
             expr = expr.expression