Repository: PyPSA/linopy

Issue:
Add automatic merging of variables
At the moment linopy will create two coefficients if the same variable is used twice in the same expression. Gurobi can resolve this automatically but it prints a warning, but HiGHS seems to struggle with such equations (see example below). A solution to this would be to merge identical variables and calculate the resulting coefficient in linopy before creating the LP file or passing the equations to the solver.

I am not entirely sure whether this is the right approach nor how other solvers are dealing with this. I don't think that HiGHS is handling this correctly, so considering this as an upstream issue in HiGHS would be also valid. But given that Gurobi prints warnings, I guess solvers prefer to have this handled by the optimization language.

## Example

Consider the following example:

```Python
N = 5
m = linopy.Model()
var = m.add_variables(
    name="var", lower=xr.DataArray(np.zeros(N), dims="x", coords={"x": np.arange(N)})
)

m.add_constraints((var - 0.5 * var) == 21.0)
m.add_objective(var.sum())
m.solve(solver_name="gurobi", keep_files=True)
```

Gurobi prints:

```
...
Warning: row 0 (name "c0") contains 1 repeated variable(s) in constraints, grouped them, 0 cancellation
Warning: row 1 (name "c1") contains 1 repeated variable(s) in constraints, grouped them, 0 cancellation
Warning: row 2 (name "c2") contains 1 repeated variable(s) in constraints, grouped them, 0 cancellation
Warning: row 3 (name "c3") contains 1 repeated variable(s) in constraints, grouped them, 0 cancellation
Warning: row 4 (name "c4") contains 1 repeated variable(s) in constraints, grouped them, 0 cancellation
Warning: lp file contains 5 repeated variable(s) in constraints, grouped them, 0 cancellation
...
```

Highs cannot solve this model at all. I aborted the solving after 10 minutes (Gurobi can solve it in less than a second).

This equivalent problem can be solved with Gurobi and HiGHS in a 100-200ms without warnings printed:

```
N = 5
m = linopy.Model()
var = m.add_variables(
    name="var", lower=xr.DataArray(np.zeros(N), dims="x", coords={"x": np.arange(N)})
)

m.add_constraints((0.5 * var) == 21.0)
m.add_objective(var.sum())
m.solve(solver_name="highs", keep_files=True)
```

So instead having coefficients `1` and `-0.5`, a preprocessing step would need to compute `0.5` as the coefficient for `var`.

At the moment these are the coefficients:

![Coefficients](https://user-images.githubusercontent.com/691772/197225481-1723547a-6ea6-4409-93f1-dc10d476f86a.png)


Discussion:

Patch:
diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index 9158982..21e1cbf 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -5,6 +5,7 @@ Upcoming Release
 ----------------
 
 * Linopy now defines and uses a fixed set of solver status and termination codes. This allows for a more consistent and reliable handling of solver results. The new codes are defined in the `linopy.consants` module. The implementation is inspired by https://github.com/0b11001111 and the implementation in this `PyPSA fork <https://github.com/0b11001111/PyPSA/blob/innoptem-lopf/pypsa/linear_program/solver.py>`_
+* The automated summation of repeated variables in one constraint is now supported. Before the implementation for constraints like `x + x + x <= 5` was only working for solvers with a corresponding fallback computation. This is now fixed.
 * Integer variables are now fully supported.
 * Support exporting problems to MPS file via fast highspy MPS-writer.
 * The internal data structure of linopy classes were updated to a safer design. Instead of being defined as inherited xarray classes, the class `Variable`, `LinearExpression` and `Constraint` are now dataclasses with containing the xarray objects in the data field. This allows the package to have more flexible function design and a reduced set of wrapped functions that are sensible to use in the optimization context.
diff --git a/linopy/constraints.py b/linopy/constraints.py
index 55ea651..f1d5972 100644
--- a/linopy/constraints.py
+++ b/linopy/constraints.py
@@ -541,7 +541,11 @@ class Constraints:
 
             shape = (ncons, nvars)  # same as model.nvars/ncons but already there
 
-        return coo_matrix((data, (rows, cols)), shape=shape)
+        # Group repeated variables in the same constraint
+        df = pd.DataFrame({"data": data, "rows": rows, "cols": cols})
+        df = df.groupby(["rows", "cols"], as_index=False).sum()
+
+        return coo_matrix((df.data, (df.rows, df.cols)), shape=shape)
 
 
 class AnonymousConstraint:
diff --git a/linopy/expressions.py b/linopy/expressions.py
index 7539290..3e9e9f0 100644
--- a/linopy/expressions.py
+++ b/linopy/expressions.py
@@ -700,6 +700,7 @@ class LinearExpression:
 
         if not np.issubdtype(self.vars.dtype, np.integer):
             return self.assign(vars=self.vars.fillna(-1).astype(int))
+
         return self
 
     def equals(self, other: "LinearExpression"):
diff --git a/linopy/io.py b/linopy/io.py
index f9bc36f..a38a0b5 100644
--- a/linopy/io.py
+++ b/linopy/io.py
@@ -11,6 +11,7 @@ from pathlib import Path
 from tempfile import TemporaryDirectory
 
 import numpy as np
+import pandas as pd
 import xarray as xr
 from numpy import asarray, concatenate, ones_like, zeros_like
 from tqdm import tqdm
@@ -102,6 +103,10 @@ def constraints_to_file(m, f, log=False):
 
         if not c.size:
             continue
+        # Group repeated variables in the same constraint
+        df = pd.DataFrame({"coeffs": c, "labels": l, "vars": v})
+        df = df.groupby(["labels", "vars"], as_index=False).sum()
+        c, l, v = df.coeffs.values, df.labels.values, df.vars.values
 
         diff_con = l[:-1] != l[1:]
         new_con_b = concatenate([asarray([True]), diff_con])
@@ -323,6 +328,9 @@ def to_highspy(m):
 def to_block_files(m, fn):
     """
     Write out the linopy model to a block structured output.
+
+    Experimental: This function does not support grouping duplicated variables
+    in one constraint row yet!
     """
     if fn is None:
         fn = TemporaryDirectory(prefix="linopy-problem-", dir=m.solver_dir).name
@@ -389,7 +397,7 @@ def to_block_files(m, fn):
 
     def filtered(arr, mask, key):
         """
-        Set coefficients to zero for mask, keep others unchanged.
+        Set coefficients to zero where mask is False, keep others unchanged.
 
         PIPS requires this information to set the shape of sub-matrices.
         """


Test patch:
diff --git a/test/test_matrices.py b/test/test_matrices.py
index 61939da..403692a 100644
--- a/test/test_matrices.py
+++ b/test/test_matrices.py
@@ -50,3 +50,16 @@ def test_basic_matrices_masked():
     assert m.matrices.clabels.shape == m.matrices.sense.shape
     assert m.matrices.vlabels.shape == m.matrices.ub.shape
     assert m.matrices.vlabels.shape == m.matrices.lb.shape
+
+
+def test_matrices_duplicated_variables():
+    m = Model()
+
+    x = m.add_variables(pd.Series([0, 0]), 1, name="x")
+    y = m.add_variables(4, pd.Series([8, 10]), name="y")
+    z = m.add_variables(0, pd.DataFrame([[1, 2], [3, 4], [5, 6]]).T, name="z")
+    m.add_constraints(x + x + y + y + z + z == 0)
+
+    A = m.matrices.A.todense()
+    assert A[0, 0] == 2
+    assert np.isin(np.unique(A), [0, 2]).all()
diff --git a/test/test_optimization.py b/test/test_optimization.py
index b3998b0..5a52432 100644
--- a/test/test_optimization.py
+++ b/test/test_optimization.py
@@ -79,6 +79,19 @@ def model_with_inf():
     return m
 
 
+@pytest.fixture
+def model_with_duplicated_variables():
+    m = Model()
+
+    lower = pd.Series(0, range(10))
+    x = m.add_variables(coords=[lower.index], name="x")
+
+    m.add_constraints(x + x, GREATER_EQUAL, 10)
+    m.objective = 1 * x
+
+    return m
+
+
 @pytest.fixture
 def milp_binary_model():
     m = Model()
@@ -220,6 +233,13 @@ def test_default_settings_chunked(model_chunked, solver, io_api):
     assert np.isclose(model_chunked.objective_value, 3.3)
 
 
+@pytest.mark.parametrize("solver,io_api", params)
+def test_duplicated_variables(model_with_duplicated_variables, solver, io_api):
+    status, condition = model_with_duplicated_variables.solve(solver, io_api=io_api)
+    assert status == "ok"
+    assert all(model_with_duplicated_variables.solution["x"] == 5)
+
+
 @pytest.mark.parametrize("solver,io_api", params)
 def test_set_files(tmp_path, model, solver, io_api):
     status, condition = model.solve(