Fix variables via bound collapse instead of equality constraint

doc/release_notes.rst

@@ -7,6 +7,8 @@ Upcoming Version
 * Add documentation about `LinearExpression.where` with `drop=True`. Add `BaseExpression.variable_names` property.
 * Add ``BaseExpression.has_terms`` property: boolean array, true at slots with at least one live term (`#741 <https://github.com/PyPSA/linopy/issues/741>`_).
 * ``Variable.fix(value)`` now places ``value`` correctly on variables with named dimensions; previously array values could be misaligned.
+* ``Variable.fix()`` now fixes a variable by collapsing its bounds (``lower = upper = value``) instead of adding a ``__fix__`` equality constraint; ``unfix()`` restores the original bounds (`#769 <https://github.com/PyPSA/linopy/issues/769>`_). A fix outside the current bounds now warns and overrides instead of raising, and its shadow price appears as the variable's reduced cost rather than a constraint dual.
+* Binary variable bounds are now respected by the solver, so fixing a binary works (they were previously forced to ``[0, 1]``).
 
 **Features**
 

examples/manipulating-models.ipynb

@@ -357,7 +357,7 @@
    "id": "30",
    "metadata": {},
    "source": [
-    "Calling `unfix()` on all variables removes the fix constraints and `unrelax()` restores the integrality of `z`."
+    "Calling `unfix()` on all variables restores their original bounds and `unrelax()` restores the integrality of `z`."
    ]
   },
   {

linopy/constants.py

@@ -37,7 +37,9 @@ class PerformanceWarning(UserWarning):
     short_LESS_EQUAL: LESS_EQUAL,
 }
 
-FIX_CONSTRAINT_PREFIX = "__fix__"
+STASHED_LOWER = "_stashed_lower"
+STASHED_UPPER = "_stashed_upper"
+STASHED_ATTRS: list[str] = [STASHED_LOWER, STASHED_UPPER]
 
 TERM_DIM = "_term"
 STACKED_TERM_DIM = "_stacked_term"

linopy/io.py

@@ -249,6 +249,9 @@ def bounds_to_file(
         list(m.variables.continuous)
         + list(m.variables.integers)
         + list(m.variables.semi_continuous)
+        + [
+            n for n in m.variables.binaries if m.variables[n].fixed
+        ]  # fixed binaries need bounds
     )
     if not len(list(names)):
         return

linopy/model.py

@@ -1149,16 +1149,8 @@ def remove_variables(self, name: str) -> None:
         -------
         None.
         """
-        from linopy.constants import FIX_CONSTRAINT_PREFIX
-
         variable = self.variables[name]
 
-        # Clean up fix constraint if present
-        fix_name = f"{FIX_CONSTRAINT_PREFIX}{name}"
-        if fix_name in self.constraints:
-            self.constraints.remove(fix_name)
-
-        # Clean up relaxed registry if present
         self._relaxed_registry.pop(name, None)
 
         to_remove = [

linopy/solvers.py

@@ -1268,12 +1268,6 @@ def _build_solver_model(
                 [integrality_map[v] for v in vtypes[int_mask]], dtype=np.int32
             )
             h.changeColsIntegrality(len(labels), labels, integrality)
-            if len(model.binaries):
-                labels = np.arange(len(vtypes))[vtypes == "B"]
-                n = len(labels)
-                h.changeColsBounds(
-                    n, labels, np.zeros_like(labels), np.ones_like(labels)
-                )
 
         c = M.c
         h.changeColsCost(len(c), np.arange(len(c), dtype=np.int32), c)
@@ -2827,9 +2821,6 @@ def _build_solver_model(
         if len(model.binaries.labels) + len(model.integers.labels) > 0:
             idx = [i for (i, v) in enumerate(M.vtypes) if v in ["B", "I"]]
             task.putvartypelist(idx, [mosek.variabletype.type_int] * len(idx))
-            if len(model.binaries.labels) > 0:
-                bidx = [i for (i, v) in enumerate(M.vtypes) if v == "B"]
-                task.putvarboundlistconst(bidx, mosek.boundkey.ra, 0.0, 1.0)
 
         if len(model.constraints) > 0:
             if set_names:

linopy/variables.py

@@ -56,10 +56,12 @@
 )
 from linopy.config import options
 from linopy.constants import (
-    FIX_CONSTRAINT_PREFIX,
     HELPER_DIMS,
     SOS_DIM_ATTR,
     SOS_TYPE_ATTR,
+    STASHED_ATTRS,
+    STASHED_LOWER,
+    STASHED_UPPER,
     TERM_DIM,
 )
 from linopy.types import (
@@ -1013,7 +1015,7 @@ def flat(self) -> DataFrame:
         -------
         df : pandas.DataFrame
         """
-        ds = self.data
+        ds = self.data.drop_vars(STASHED_ATTRS, errors="ignore")
 
         def mask_func(data: pd.DataFrame) -> pd.Series:
             return data["labels"] != -1
@@ -1033,7 +1035,8 @@ def to_polars(self) -> pl.DataFrame:
         -------
         pl.DataFrame
         """
-        df = to_polars(self.data)
+        ds = self.data.drop_vars(STASHED_ATTRS, errors="ignore")
+        df = to_polars(ds)
         df = filter_nulls_polars(df)
         check_has_nulls_polars(df, name=f"{self.type} {self.name}")
         return df
@@ -1341,10 +1344,16 @@ def fix(
         overwrite: bool = True,
     ) -> None:
         """
-        Fix the variable to a given value by adding an equality constraint.
+        Fix the variable to a given value by collapsing its bounds.
+
+        Sets ``lower = upper = value``.
 
         If no value is given, the current solution value is used.
 
+        A fix value outside the variable's current bounds emits a warning, but
+        does not cause infeasibilities (the bounds are overridden). Fixing a
+        binary variable to anything other than 0 or 1 raises.
+
         Parameters
         ----------
         value : float/array_like, optional
@@ -1354,8 +1363,9 @@ def fix(
             Integer and binary variables are always rounded to 0 decimal places.
             Default is 8.
         overwrite : bool, optional
-            If True (default), overwrite an existing fix constraint for this
-            variable. If False, raise an error if the variable is already fixed.
+            If True (default), re-fix a variable that is already fixed to the
+            new value (the originally stashed bounds are kept). If False, raise
+            an error if the variable is already fixed.
         """
         if value is None:
             try:
@@ -1368,50 +1378,72 @@ def fix(
                 )
                 raise ValueError(msg) from None
 
+        is_fixed = self.fixed
+        is_binary = self.attrs["binary"]
+        is_integer = self.attrs["integer"]
+
+        if is_fixed and not overwrite:
+            msg = (
+                f"Variable '{self.name}' is already fixed. Use "
+                "overwrite=True to replace the existing fix value."
+            )
+            raise ValueError(msg)
+
         value = broadcast_to_coords(
             value, self.coords, label=f"fix() for variable '{self.name}'"
         )
 
-        if self.attrs.get("integer") or self.attrs.get("binary"):
+        if is_binary and not (np.isclose(value, 0) | np.isclose(value, 1)).all():
+            msg = (
+                f"Cannot fix binary variable '{self.name}' to a value "
+                "other than 0 or 1."
+            )
+            raise ValueError(msg)
+
+        if is_integer or is_binary:
             value = value.round(0)
         else:
             value = value.round(decimals)
 
-        if (value < self.lower).any() or (value > self.upper).any():
-            msg = (
-                f"Fix values for variable '{self.name}' are outside the "
-                "variable bounds."
-            )
-            raise ValueError(msg)
+        if is_fixed:
+            lower, upper = self.data[STASHED_LOWER], self.data[STASHED_UPPER]
+        else:
+            lower, upper = self.data.lower, self.data.upper
 
-        constraint_name = f"{FIX_CONSTRAINT_PREFIX}{self.name}"
+        if not is_binary and ((value < lower).any() or (value > upper).any()):
+            warn(
+                f"Fix values for variable '{self.name}' lie outside its current "
+                "bounds; the bounds are overridden by the fix value.",
+                UserWarning,
+                stacklevel=2,
+            )
 
-        if constraint_name in self.model.constraints:
-            if not overwrite:
-                msg = (
-                    f"Variable '{self.name}' is already fixed. Use "
-                    "overwrite=True to replace the existing fix constraint."
-                )
-                raise ValueError(msg)
-            self.model.remove_constraints(constraint_name)
+        if not is_fixed:
+            self._data = assign_multiindex_safe(
+                self.data,
+                **{STASHED_LOWER: lower, STASHED_UPPER: upper},
+            )
 
-        self.model.add_constraints(self, "=", value, name=constraint_name)
+        self.lower = value
+        self.upper = value
 
     def unfix(self) -> None:
         """
-        Remove the fix constraint for this variable.
+        Unfix the variable, restoring the bounds it had before :meth:`fix`.
         """
-        constraint_name = f"{FIX_CONSTRAINT_PREFIX}{self.name}"
-        if constraint_name in self.model.constraints:
-            self.model.remove_constraints(constraint_name)
+        if not self.fixed:
+            return
+
+        self.lower = self.data[STASHED_LOWER]
+        self.upper = self.data[STASHED_UPPER]
+        self._data = self.data.drop_vars(STASHED_ATTRS)
 
     @property
     def fixed(self) -> bool:
         """
         Return whether the variable is currently fixed.
         """
-        constraint_name = f"{FIX_CONSTRAINT_PREFIX}{self.name}"
-        return constraint_name in self.model.constraints
+        return all(attr in self.data for attr in STASHED_ATTRS)
 
 
 class AtIndexer:
@@ -1727,7 +1759,7 @@ def fix(
         decimals : int, optional
             Number of decimal places to round continuous variables to.
         overwrite : bool, optional
-            If True, overwrite existing fix constraints.
+            If True, re-fix variables that are already fixed.
         """
         for var in self.data.values():
             var.fix(value=value, decimals=decimals, overwrite=overwrite)