fix: ignore masked variable labels when matching constraints

doc/release_notes.rst

@@ -56,6 +56,7 @@ Version v0.9.0
 * ``linopy.testing.assert_linequal`` now aligns dimension order before comparing, so mathematically identical expressions built in different orders (e.g. ``x + y`` versus ``y + x``, which inherit different dimension orders from xarray broadcasting) are correctly treated as equal. Genuinely different expressions still fail. (`#801 <https://github.com/PyPSA/linopy/pull/801>`__)
 * Summing an expression over a dimension that carries an auxiliary (non-dimension) coordinate no longer leaks that coordinate onto the internal term dimension, where it broke later arithmetic with a ``CoordinateValidationError``. Auxiliary coordinates on the remaining dimensions still propagate. (`#295 <https://github.com/PyPSA/linopy/issues/295>`__)
 * ``Solver.close()`` (also triggered by ``model.solver = None`` and the next ``solve()`` call) now explicitly disposes the ``gurobipy`` model before the environment. Previously the model was only dereferenced, so a user-held ``model.solver_model`` reference silently kept the Gurobi license acquired after ``close()``. (`#459 <https://github.com/PyPSA/linopy/issues/459>`__)
+* ``Model.remove_variables`` no longer removes constraints that never reference the removed variable. A masked variable carries ``-1`` label entries, which matched the ``-1`` that marks an empty term slot in a constraint, so any masked variable looked like it was used by any constraint with padded terms. Models built with ``mask=`` could silently lose constraints and solve to a wrong optimum. (`#883 <https://github.com/PyPSA/linopy/issues/883>`__)
 
 Version 0.8.0
 -------------

linopy/common.py

@@ -471,6 +471,19 @@ def replace_by_map(ds: DataArray, mapping: np.ndarray) -> DataArray:
     )
 
 
+def assigned_labels(labels: np.ndarray | DataArray) -> np.ndarray:
+    """
+    Flatten labels and drop the -1 sentinels.
+
+    ``-1`` marks a masked entry in a variable's or constraint's labels, but an
+    empty term slot in the ``vars`` field of a constraint or expression.
+    Matching labels of one object against another must therefore ignore it,
+    otherwise every masked entry compares equal to every empty term slot.
+    """
+    flat = np.asarray(labels).ravel()
+    return flat[flat != -1]
+
+
 def to_path(path: str | Path | None) -> Path | None:
     """
     Convert a string to a Path object.

linopy/constraints.py

@@ -38,6 +38,7 @@
     VariableLabelIndex,
     align_lines_by_delimiter,
     assign_multiindex_safe,
+    assigned_labels,
     check_has_nulls,
     check_has_nulls_polars,
     coords_from_dataset,
@@ -202,7 +203,12 @@ def data_attrs(self) -> list[str]:
 
     @abstractmethod
     def has_variable(self, variable: variables.Variable) -> bool:
-        """Check if the constraint references any of the given variable labels."""
+        """
+        Check if the constraint references any of the given variable labels.
+
+        Masked variable entries (label -1) are ignored: they are not part of
+        the model and would otherwise match every empty term slot.
+        """
 
     @abstractmethod
     def sanitize_zeros(self) -> ConstraintBase:
@@ -943,7 +949,7 @@ def from_netcdf_ds(cls, ds: Dataset, model: Model, name: str) -> CSRConstraint:
     def has_variable(self, variable: variables.Variable) -> bool:
         vlabels = self._model.variables.label_index.vlabels
         return bool(
-            np.isin(vlabels[self._csr.indices], variable.labels.values.ravel()).any()
+            np.isin(vlabels[self._csr.indices], assigned_labels(variable.labels)).any()
         )
 
     def to_matrix_with_rhs(
@@ -1520,7 +1526,7 @@ def dual(self, value: ConstantLike) -> None:
         self._data = assign_multiindex_safe(self.data, dual=value)
 
     def has_variable(self, variable: variables.Variable) -> bool:
-        return bool(self.data["vars"].isin(variable.labels.values.ravel()).any())
+        return bool(self.data["vars"].isin(assigned_labels(variable.labels)).any())
 
     def _matrix_export_data(
         self, label_index: VariableLabelIndex

linopy/model.py

@@ -31,6 +31,7 @@
 from linopy.alignment import as_dataarray, broadcast_to_coords
 from linopy.common import (
     assign_multiindex_safe,
+    assigned_labels,
     best_int,
     maybe_replace_signs,
     replace_by_map,
@@ -1448,7 +1449,7 @@ def remove_variables(self, name: str) -> None:
         self.variables.remove(name)
 
         self.objective = self.objective.sel(
-            {TERM_DIM: ~self.objective.vars.isin(variable.labels)}
+            {TERM_DIM: ~self.objective.vars.isin(assigned_labels(variable.labels))}
         )
 
     def remove_constraints(self, name: str | list[str]) -> None:

test/test_model.py

@@ -11,6 +11,7 @@
 from tempfile import gettempdir
 
 import numpy as np
+import pandas as pd
 import pytest
 import xarray as xr
 
@@ -160,6 +161,36 @@ def test_remove_variable() -> None:
     assert not m.objective.vars.isin(x.labels).any()
 
 
+@pytest.mark.parametrize("freeze", [False, True])
+def test_remove_masked_variable_keeps_unrelated_constraints(freeze: bool) -> None:
+    # https://github.com/PyPSA/linopy/issues/883
+    m: Model = Model(freeze_constraints=freeze)
+
+    i = pd.Index(range(3), name="i")
+    mask = [True, False, True]
+    a = m.add_variables(coords=[i], name="a", mask=mask)
+    b = m.add_variables(coords=[i], name="b", mask=mask)
+    c = m.add_variables(coords=[i], name="c")
+
+    # `b` is masked, so the constraint carries empty term slots (-1), but it
+    # never references `a`
+    without_a = m.add_constraints(b.sum() + c, EQUAL, 0, name="without_a")
+    assert not without_a.has_variable(a)
+
+    with_a = m.add_constraints(a.sum() + c, EQUAL, 0, name="with_a")
+    assert with_a.has_variable(a)
+
+    m.add_objective((1 * a).sum() + (1 * c).sum())
+
+    with pytest.warns(UserWarning, match="with_a"):
+        m.remove_variables("a")
+
+    assert "without_a" in m.constraints
+    assert "with_a" not in m.constraints
+    assert not m.objective.vars.isin(a.labels[a.labels != -1]).any()
+    assert m.objective.vars.isin(c.labels).any()
+
+
 def test_remove_constraint() -> None:
     m: Model = Model()