Repository: PyPSA/linopy

Issue:
Shuffled constraints and variables causing infeasibilities
Dear linopy team,

I am currently using linopy 0.1.4 with PyPSA 0.22.1

While trying to solve a small optimisation model via PyPSA, I encountered some infeasibilities on a model which shouldn't have any. The issue happens in the specific case where I combine unit commitment, ramp constraints and rolling horizon - the first iteration always goes well, and the second is found infeasible. 

I tried to dig into the PyPSA and linopy code and it seems like the constraints are well-written all along the way, except at the very end during the export to the lp file. In the LP file, right and left hand terms are mixed, grouping part of constraints that are not related, hence causing the infeasibilities.
When rolling horizon is disabled, it does not lead to infeasibilities but constraints are still not written properly...
However, the to_file linopy function is a bit beyond my capabilities, and it is likely I missed something upstream, too. 

Here is the code used for the tests, the logs, and the dataset. Thanks a lot for your help!

```
n = pypsa.Network()
    n.import_from_csv_folder("main/test")

    n_periods = 1000
    for i in range(n_periods):
        print("Iteration " + str(i+1) + " over " + str(n_periods))
        begin = i * int(len(n.snapshots) / n_periods)
        end = min((i + 1) * int(len(n.snapshots) / n_periods + 1) + 2, len(n.snapshots))

        m = n.optimize.create_model(snapshots=n.snapshots[begin: end], linearized_unit_commitment=False)

        m.to_file("pypsa.lp")
        n.optimize.solve_model(solver_name='gurobi')

        try:
            se = m.compute_set_of_infeasible_constraints()
            print(se)
        except:
            print("worked")
```
```
IIS runtime: 0.00 seconds (0.00 work units)
 c112: x8 <= 0

 c113: x1 <= 10000

 c118: x4 <= 10000

 c294: x0 <= 10000

 c295: x2 <= 5000

 c296: x3 <= 5000

 c297: x5 <= 7500
[test (1).zip](https://github.com/PyPSA/linopy/files/11145974/test.1.zip)


 c299: x6 <= 7500

 c298: x7 <= 7500

 c574: x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 = 200000
```



Discussion:

Patch:
diff --git a/linopy/expressions.py b/linopy/expressions.py
index ae0be11..21ebe5f 100644
--- a/linopy/expressions.py
+++ b/linopy/expressions.py
@@ -886,6 +886,8 @@ class LinearExpression:
 
     reindex = exprwrap(Dataset.reindex, fill_value=_fill_value)
 
+    reindex_like = exprwrap(Dataset.reindex_like, fill_value=_fill_value)
+
     rename = exprwrap(Dataset.rename)
 
     rename_dims = exprwrap(Dataset.rename_dims)
diff --git a/linopy/model.py b/linopy/model.py
index 79e8427..cb4b60f 100644
--- a/linopy/model.py
+++ b/linopy/model.py
@@ -472,7 +472,7 @@ class Model:
         labels_reindexed = labels.reindex_like(self.variables.labels, fill_value=-1)
         if not labels.equals(labels_reindexed):
             warnings.warn(
-                f"Reindexing variable {name} to match existing coordinates.",
+                f"Reindexing variable `{name}` to match existing coordinates.",
                 UserWarning,
             )
             labels = labels_reindexed
@@ -593,6 +593,29 @@ class Model:
         if labels is None:
             labels = (lhs.vars.chunk() + rhs).sum("_term")
 
+        if mask is not None:
+            mask = DataArray(mask).astype(bool)
+            mask, _ = xr.align(mask, labels, join="right")
+            assert set(mask.dims).issubset(
+                labels.dims
+            ), "Dimensions of mask not a subset of resulting labels dimensions."
+
+        # It is important to end up with monotonically increasing labels in the
+        # model's constraints container as we use it for indirect indexing.
+        labels_reindexed = labels.reindex_like(self.constraints.labels, fill_value=-1)
+        if not labels.equals(labels_reindexed):
+            warnings.warn(
+                f"Reindexing constraint `{name}` to match existing coordinates.",
+                UserWarning,
+            )
+            labels = labels_reindexed
+            lhs = lhs.reindex_like(labels)
+            rhs = rhs.reindex_like(labels)
+            if mask is None:
+                mask = labels != -1
+            else:
+                mask = mask.reindex_like(labels_reindexed, fill_value=False)
+
         self.check_force_dim_names(labels)
 
         start = self._cCounter
@@ -600,10 +623,6 @@ class Model:
         self._cCounter += labels.size
 
         if mask is not None:
-            mask = DataArray(mask)
-            assert set(mask.dims).issubset(
-                labels.dims
-            ), "Dimensions of mask not a subset of resulting labels dimensions."
             labels = labels.where(mask, -1)
 
         lhs = lhs.data.rename({"_term": f"{name}_term"})


Test patch:
diff --git a/test/test_constraints.py b/test/test_constraints.py
index 59938e8..43b0d33 100644
--- a/test/test_constraints.py
+++ b/test/test_constraints.py
@@ -22,8 +22,8 @@ def test_constraint_assignment():
 
     lower = xr.DataArray(np.zeros((10, 10)), coords=[range(10), range(10)])
     upper = xr.DataArray(np.ones((10, 10)), coords=[range(10), range(10)])
-    x = m.add_variables(lower, upper)
-    y = m.add_variables()
+    x = m.add_variables(lower, upper, name="x")
+    y = m.add_variables(name="y")
 
     m.add_constraints(1 * x + 10 * y, EQUAL, 0)
 
@@ -42,8 +42,8 @@ def test_anonymous_constraint_assignment():
 
     lower = xr.DataArray(np.zeros((10, 10)), coords=[range(10), range(10)])
     upper = xr.DataArray(np.ones((10, 10)), coords=[range(10), range(10)])
-    x = m.add_variables(lower, upper)
-    y = m.add_variables()
+    x = m.add_variables(lower, upper, name="x")
+    y = m.add_variables(name="y")
     con = 1 * x + 10 * y == 0
     m.add_constraints(con)
 
@@ -86,6 +86,26 @@ def test_constraint_assignment_chunked():
     assert isinstance(m.constraints.coeffs.c.data, dask.array.core.Array)
 
 
+def test_constraint_assignment_with_reindex():
+    m = Model()
+
+    lower = xr.DataArray(np.zeros((10, 10)), coords=[range(10), range(10)])
+    upper = xr.DataArray(np.ones((10, 10)), coords=[range(10), range(10)])
+    x = m.add_variables(lower, upper, name="x")
+    y = m.add_variables(name="y")
+
+    m.add_constraints(1 * x + 10 * y, EQUAL, 0)
+
+    shuffled_coords = pd.Index([2, 1, 3, 4, 6, 5, 7, 9, 8, 0], name="first")
+    con = m.variables["x"].loc[shuffled_coords] <= 10
+    with pytest.warns(UserWarning):
+        con_shuffled = m.add_constraints(con, name="con_shuffled")
+
+    assert con_shuffled.indexes["dim_0"].equals(m.constraints["con0"].indexes["dim_0"])
+    # check if labels are monotonically increasing
+    assert np.all(np.diff(np.ravel(con_shuffled.labels)) > 0)
+
+
 def test_wrong_constraint_assignment_repeated():
     # repeated variable assignment is forbidden
     m = Model()
diff --git a/test/test_variables.py b/test/test_variables.py
index a678fb6..079825e 100644
--- a/test/test_variables.py
+++ b/test/test_variables.py
@@ -45,6 +45,15 @@ def test_variables_assignment_with_merge():
     assert m.variables.labels.var0[-1].item() != -1
 
 
+def test_variables_assignment_with_reindex(m):
+    shuffled_coords = [pd.Index([2, 1, 3, 4, 6, 5, 7, 9, 8, 0], name="first")]
+    with pytest.warns(UserWarning):
+        a = m.add_variables(coords=shuffled_coords, name="a")
+    assert a.indexes["first"].equals(m["x"].indexes["first"])
+    # check if labels are monotonically increasing
+    assert np.all(np.diff(np.ravel(a.labels)) > 0)
+
+
 def test_scalar_variables_name_counter():
     m = Model()
     m.add_variables()