Add automatic merging of variables

[lumbric]

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:

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:

[FabianHofmann]

@lumbric thanks, I saw that issue coming :)

It would be good to solve that via a pandas.groupby operation that groups coefficients of variable labels together. I'll put it on my todo list.