Auto-convert hard to soft constraints

API — operate on an already-defined constraint

The phrase "generalize softening of a defined constraint" points at a method on the constraint object, with a convenience kwarg at creation:

class Constraint:
    def soften(
        self,
        penalty,                 # ConstantLike: scalar, or DataArray aligned to (a subset of) the constraint's coords
        *,
        max_violation=None,      # ConstantLike | None — caps the slack (upper bound); None = unbounded
        name=None,               # slack variable name; default f"{self.name}_slack"
    ) -> Variable | tuple[Variable, Variable]:
        ...

# soften after the fact
c = m.add_constraints(emissions <= budget, name="co2_limit")
slack = c.soften(penalty=1e4)

# or inline at creation (thin wrapper that calls .soften)
m.add_constraints(emissions <= budget, name="co2_limit", penalty=1e4)

# post-solve: how much each element was violated
slack.solution

Returns the slack Variable (a (s_pos, s_neg) tuple for equalities) so the user can read violations and reuse the handle. penalty >= 0 is asserted.

Implementation sketch (uses only existing primitives)

def soften(self, penalty, *, max_violation=None, name=None):
    m = self.model
    assert (penalty >= 0).all() if hasattr(penalty, "all") else penalty >= 0
    name = name or f"{self.name}_slack"
    upper = inf if max_violation is None else max_violation

    s = m.add_variables(lower=0, upper=upper, coords=self.lhs.coords,
                        mask=self.mask, name=name)

    if self.sign == "<=":
        self.lhs = self.lhs - s            # .lhs setter already exists
    elif self.sign == ">=":
        self.lhs = self.lhs + s
    else:                                   # "="
        s_neg = m.add_variables(lower=0, upper=upper, coords=self.lhs.coords,
                                mask=self.mask, name=f"{name}_neg")
        self.lhs = self.lhs - s + s_neg
        s = (s, s_neg)

    viol = s[0] + s[1] if isinstance(s, tuple) else s
    direction = 1 if m.sense == "min" else -1
    m.objective += direction * (penalty * viol).sum()
    return s

No new solver machinery — it's variables + an lhs rewrite + an objective term, all of which linopy already supports. Rewriting .lhs in place (rather than remove_constraints + re-add) keeps the constraint's name and dual intact.

Design points worth flagging in the issue

1. max_violation = the "penalty quantity" knob — a hard cap on the cheap slack. Calling soften twice with increasing penalty and caps yields a piecewise penalty curve; a future penalty=[(qty, price), …] form could express that directly. Keep v1 to single-tier + optional cap.
2. Masks and per-element penalties are first-class because slack inherits c.coords/c.mask — masked-out entries get no slack, and penalty broadcasts over any coord subset.
3. Parameter, not variable — penalty (and max_violation) must be data; allowing a Variable penalty makes the objective bilinear. Reject for now.
4. Round-trip — record softened constraints (slack name + penalty) in model metadata so to_netcdf/from_netcdf preserves them, and expose c.is_soft / m.constraints.softened.
5. QP-safe — the penalty term is linear and adds cleanly to a quadratic objective.