diff --git a/doc/release_notes.rst b/doc/release_notes.rst index dae3a6dd4..626bf1cf9 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -29,6 +29,11 @@ Upcoming Version * Added support for the GPU-accelerated `NVIDIA 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) diff --git a/linopy/constraints.py b/linopy/constraints.py index f3b301cce..747377df3 100644 --- a/linopy/constraints.py +++ b/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. diff --git a/linopy/model.py b/linopy/model.py index 614adf387..bf6f4c009 100644 --- a/linopy/model.py +++ b/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, diff --git a/linopy/objective.py b/linopy/objective.py index e7bb6dc76..dc2e8ad41 100644 --- a/linopy/objective.py +++ b/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 diff --git a/test/test_constraint.py b/test/test_constraint.py index 24ae7f936..c0f0e41be 100644 --- a/test/test_constraint.py +++ b/test/test_constraint.py @@ -30,6 +30,7 @@ ConstraintBase, Constraints, ) +from linopy.testing import assert_linequal, assert_varequal @pytest.fixture @@ -95,6 +96,41 @@ def test_add_constraints_uses_model_freeze_default() -> None: ) +def test_add_constraints_penalty_softens_constraint( + m: Model, y: linopy.Variable +) -> None: + """`penalty=` on add_constraints is a shortcut for calling `.soften()`.""" + m.add_objective(y.sum(), sense="min") + penalty_coeff = 10 + constraint = m.add_constraints(y >= 10, name="constraint", penalty=penalty_coeff) + + assert isinstance(constraint, linopy.constraints.Constraint) + # The slack term was added to the lhs, same effect as calling .soften() directly: + assert constraint.lhs.nterm == 2 + expected_objective = ( + y.sum() + penalty_coeff * m.variables["constraint_slack_pos"].sum() + ) + assert_linequal(m.objective.expression, expected_objective) + + +def test_add_constraints_penalty_with_freeze_true_raises( + m: Model, x: linopy.Variable +) -> None: + with pytest.raises(ValueError, match="`penalty` cannot be combined"): + m.add_constraints(x >= 0, name="frozen_penalized", freeze=True, penalty=10) + + +def test_add_constraints_penalty_with_model_freeze_default_raises() -> None: + """ + `freeze=None` resolves to the model's `freeze_constraints` default, which must + also be checked against `penalty`, not just an explicit `freeze=True`. + """ + m = Model(freeze_constraints=True) + x = m.add_variables(coords=[pd.RangeIndex(10, name="first")], name="x") + with pytest.raises(ValueError, match="`penalty` cannot be combined"): + m.add_constraints(x >= 0, name="frozen_by_default_penalized", penalty=10) + + def test_constraint_name(c: linopy.constraints.CSRConstraint) -> None: assert c.name == "c" @@ -1119,3 +1155,260 @@ def bound(m: Model, i: int) -> AnonymousScalarConstraint: r = repr(con) assert "≥" in r assert "=" in r + + +# Constraint.soften method's tests + + +def test_constraint_soften_returns_slack_for_le_and_ge( + m: Model, y: linopy.Variable +) -> None: + """ + Checks that a constraint of type 'less or equal' or 'great and equal' returns + one slack each. + """ + z = m.variables["z"] + m.add_objective((z + y).sum(), sense="min") + + # create new constraint, and add slack for 'greater or equal' + z_constraint = m.add_constraints(z >= -10, name="constraint_over_z") + z_slack = z_constraint.soften(penalty=10) + assert isinstance(z_slack.positive, linopy.Variable) + assert z_slack.negative is None + + # create new constraint, and add slack for 'less or equal': + y_constraint = m.add_constraints(y <= 0, name="constraint_over_y") + y_slack = y_constraint.soften(penalty=10) + assert isinstance(y_slack.positive, linopy.Variable) + assert y_slack.negative is None + + +def test_constraint_soften_returns_slack_for_eq(m: Model, y: linopy.Variable) -> None: + """Checks that a constraint of type 'equal'' returns two slack variables""" + m.add_objective(y.sum(), sense="min") + y_contraint = m.add_constraints(y == 10, name="equality_contraint") + y_slack = y_contraint.soften(penalty=10) + + assert isinstance(y_slack.positive, linopy.Variable) + assert isinstance(y_slack.negative, linopy.Variable) + + +def test_constraint_soften_raises_on_detached_mutable_constraint( + m: Model, x: linopy.Variable, mc: linopy.Constraint +) -> None: + """ + `.mutable()` on a frozen constraint returns a detached copy that is not + registered in `model.constraints`; softening it would silently fail to + affect the actual model, so `soften` must raise instead. + """ + m.add_objective(x.sum(), sense="min") + + with pytest.raises(ValueError, match="not the constraint registered"): + # mc := m.constraints["c"].mutable(), detached from the model + mc.soften(penalty=10) + + +def test_constraint_soften_twice_raises(m: Model, y: linopy.Variable) -> None: + """ + Softening an already-softened constraint must raise instead of silently + stacking a second, redundant slack term onto the same lhs. + """ + m.add_objective(y.sum(), sense="min") + constraint = m.add_constraints(y >= 10, name="constraint") + constraint.soften(penalty=10, name="first_slack") + + with pytest.raises(ValueError, match="already softened"): + constraint.soften(penalty=5, name="second_slack") + + +def test_constraint_soften_twice_raises_via_add_constraints_penalty( + m: Model, y: linopy.Variable +) -> None: + """The same guard applies when the first soften came from `add_constraints(penalty=...)`.""" + m.add_objective(y.sum(), sense="min") + constraint = m.add_constraints(y >= 10, name="constraint", penalty=10) + + with pytest.raises(ValueError, match="already softened"): + constraint.soften(penalty=5, name="second") + + +def test_constraint_soften_updates_lhs(m: Model, y: linopy.Variable) -> None: + """ + Tests that the left hand side gains the slack term(s) with the expected sign + per constraint direction. + """ + m.add_objective(y.sum(), sense="min") + + # '<=' : lhs -> lhs - positive_slack + le_constraint = m.add_constraints(y <= 10, name="le_constraint") + le_slack = le_constraint.soften(penalty=10) + + # Assert that the new lhs has 1 extra term: + assert le_constraint.lhs.nterm == 2 + + # Assert tht the label for the constraint's new variable is the same as the slack: + assert_equal( + le_constraint.vars.isel({le_constraint.term_dim: -1}), + le_slack.positive.labels, + ) + + # Since the it's less or equal, we expect the coeff of the variable to be -1 + assert (le_constraint.coeffs.isel({le_constraint.term_dim: -1}) == -1).all() + + # '>=' : lhs -> lhs + positive_slack + ge_constraint = m.add_constraints(y >= -10, name="ge_constraint") + ge_slack = ge_constraint.soften(penalty=10) + + assert ge_constraint.lhs.nterm == 2 + + assert_equal( + ge_constraint.vars.isel({ge_constraint.term_dim: -1}), + ge_slack.positive.labels, + ) + assert (ge_constraint.coeffs.isel({ge_constraint.term_dim: -1}) == 1).all() + + # '==' : lhs -> lhs - positive_slack + negative_slack + eq_constraint = m.add_constraints(y == 0, name="eq_constraint") + eq_slack = eq_constraint.soften(penalty=10) + assert eq_slack.negative is not None + + # Assert there's two extra terms: + assert eq_constraint.lhs.nterm == 3 + + # The positive slack is added first to the lhs, therefore, it should be on position + # -2 of the lhs: + assert_equal( + eq_constraint.vars.isel({eq_constraint.term_dim: -2}), + eq_slack.positive.labels, + ) + + # The negative slack is added next to the lhs, therefore, it should be on position + # -1 of the lhs: + assert_equal( + eq_constraint.vars.isel({eq_constraint.term_dim: -1}), + eq_slack.negative.labels, + ) + + assert (eq_constraint.coeffs.isel({eq_constraint.term_dim: -2}) == -1).all() + assert (eq_constraint.coeffs.isel({eq_constraint.term_dim: -1}) == 1).all() + + +def test_constraint_soften_updates_objective_min_sense( + m: Model, y: linopy.Variable +) -> None: + penalty_coeff = 10 + original_objective = y.sum() + m.add_objective(original_objective, sense="min") + constraint = m.add_constraints(y >= 10, name="constraint") + slack = constraint.soften(penalty=penalty_coeff) + + # For sense='min', the penalty must added with positive sign: + expected_objective = original_objective + penalty_coeff * slack.positive.sum() + assert_linequal(m.objective.expression, expected_objective) + + +def test_constraint_soften_updates_objective_max_sense( + m: Model, y: linopy.Variable +) -> None: + penalty_coeff = 10 + original_objective = y.sum() + m.add_objective(original_objective, sense="max") + constraint = m.add_constraints(y <= 10, name="constraint") + slack = constraint.soften(penalty=penalty_coeff) + + # For sense='max', the penalty must added with negative sign: + expected_objective = original_objective - penalty_coeff * slack.positive.sum() + assert_linequal(m.objective.expression, expected_objective) + + +def test_constraint_soften_max_violation_bounds_slack( + m: Model, x: linopy.Variable +) -> None: + """max_violation sets the slack variable's upper bound; default is unbounded (inf).""" + m.add_objective(x.sum(), sense="min") + + bounded_constraint = m.add_constraints(x >= 0, name="bounded_constraint") + bounded_slack = bounded_constraint.soften(penalty=10, max_violation=5) + assert (bounded_slack.positive.upper == 5).all() + + unbounded_constraint = m.add_constraints(x >= 0, name="unbounded_constraint") + unbounded_slack = unbounded_constraint.soften(penalty=10) + assert np.isinf(unbounded_slack.positive.upper).all() + + +def test_constraint_soften_negative_penalty_raises( + m: Model, y: linopy.Variable +) -> None: + """Setting penalty < 0 raises ValueError.""" + m.add_objective(y.sum(), sense="min") + constraint = m.add_constraints(y >= 10, name="constraint") + with pytest.raises(ValueError, match="Penalty is not positive"): + constraint.soften(penalty=-10) + + +def test_constraint_soften_without_objective_raises( + m: Model, y: linopy.Variable +) -> None: + """Calling soften before model.add_objective raises ValueError.""" + constraint = m.add_constraints(y <= 10, name="constraint") + with pytest.raises(ValueError, match="Objective must be defined"): + constraint.soften(penalty=10) + + +def test_constraint_soften_respects_mask(m: Model, x: linopy.Variable) -> None: + """Masked constraint entries produce masked slack variables (mask propagated).""" + m.add_objective(x.sum(), sense="min") + mask = pd.Series([False] * 5 + [True] * 5) + constraint = m.add_constraints(x >= 0, name="masked_constraint", mask=mask) + slack = constraint.soften(penalty=10) + + assert_equal(slack.positive.mask, constraint.mask) + assert constraint.mask is not None + assert (constraint.mask.values == mask.values).all() + + +def test_constraint_slack_none_before_soften(m: Model, x: linopy.Variable) -> None: + """A constraint that was never softened exposes no slack variable.""" + constraint = m.add_constraints(x >= 0, name="constraint") + assert constraint.slack is None + + +def test_constraint_slack_matches_returned_slack_for_le_and_ge( + m: Model, x: linopy.Variable, y: linopy.Variable +) -> None: + """ + `.slack` resolves to the same variables `soften` returned, for both + inequality directions (no negative slack). + """ + m.add_objective((x + y).sum(), sense="min") + + le_constraint = m.add_constraints(y <= 0, name="le_constraint") + le_slack = le_constraint.soften(penalty=10) + assert le_constraint.slack is not None + assert_varequal(le_constraint.slack.positive, le_slack.positive) + assert le_constraint.slack.negative is None + + ge_constraint = m.add_constraints(x >= -10, name="ge_constraint") + ge_slack = ge_constraint.soften(penalty=10) + assert ge_constraint.slack is not None + assert_varequal(ge_constraint.slack.positive, ge_slack.positive) + assert ge_constraint.slack.negative is None + + +def test_constraint_slack_matches_returned_slack_for_eq( + m: Model, y: linopy.Variable +) -> None: + """ + For an equality constraint, `.slack` carries both the positive and the + negative slack variable. + """ + m.add_objective(y.sum(), sense="min") + constraint = m.add_constraints(y == 10, name="eq_constraint") + slack = constraint.soften(penalty=10) + + resolved = constraint.slack + assert resolved is not None + assert_varequal(resolved.positive, slack.positive) + assert slack.negative is not None + assert resolved.negative is not None + assert_varequal(resolved.negative, slack.negative)