Replace pyomo with linopy

[GbotemiB]

Closes #1320 .

Changes proposed in this Pull Request

This PR aim to remove pyomo from cluster_network.py and replace it with linopy. Scip was also introduced to be used in cases solver (glpk and cbc) does not support quadratic objectives.

The changes were adapted from pypsa-eur

Checklist

• I consent to the release of this PR's code under the AGPLv3 license and non-code contributions under CC0-1.0 and CC-BY-4.0.
• I tested my contribution locally and it seems to work fine.
• Code and workflow changes are sufficiently documented.
• Newly introduced dependencies are added to envs/environment.yaml and doc/requirements.txt.
• Changes in configuration options are added in all of config.default.yaml and config.tutorial.yaml.
• Add a test config or line additions to test/ (note tests are changing the config.tutorial.yaml)
• Changes in configuration options are also documented in doc/configtables/*.csv and line references are adjusted in doc/configuration.rst and doc/tutorial.rst.
• A note for the release notes doc/release_notes.rst is amended in the format of previous release notes, including reference to the requested PR.

[ekatef]

Hello @GbotemiB, and many thanks for finding a solution for #1320! A great catch and nice implementation.

I have only one major question which relates to a way how the objective is formulated. Could you please check it?

It could be also great to make sure that those parts in the code are well documented as the current solution looks like a temporary one on PyPSA/linopy level, and we should be ready to further revisions. Have added a couple of minor suggestions on it, but obviously feel free to add any insights you have found in course of your investigation

[ekatef]

A great idea to improve the comment focusing on the reason to leave out this term since it's constant!

Can the explanation contain a bit more details? My first thought has been that (L * n_clusters) ** 2 should be multiplied by constant 🙂

[GbotemiB]

I also think so, but for some reason it was commented out. Maybe @FabianHofmann can provide more insights here since he worked on it in pypsa-eur

[ekatef]

@GbotemiB ahh, I see! You have checked the PR 😄
It appears that Fabian has just used a quick formulation in the comment which has been supplemented with a bit more details after. Can we do the same? A possible revision:

# leave out constant in objective (L * n_clusters) ** 2 as it doesn't affect the clustering results

[ekatef]

Super @GbotemiB, almost there!
The only point left is to revert a change in the environment name and make sure that you are happy with the comment on the constant in the objective.

[GbotemiB]

Thank you @ekatef, I have revert the changes. The comments looks good to me.

#Edited
I still have one more action, which is to test highs solver.

[GbotemiB]

I tried with highs, I got a few errors. I got it to work by setting the option to "ipm", but the results werent the same.

Here is the error by using default highs

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options:
 - solver: simple
INFO:linopy.io: Writing time: 0.03s
Running HiGHS 1.9.0 (git hash: n/a): Copyright (c) 2024 HiGHS under MIT licence terms
INFO:linopy.solvers:Log file at /private/var/folders/c3/657p__892k72qxh6r8fdkqsh0000gn/T/highs.log
WARNING: Value "simple" for solver option is not one of "simplex", "choose", "ipm" or "pdlp"
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e-01, 7e+00]
  Bound  [1e+00, 8e+01]
  RHS    [4e+00, 4e+00]
ERROR:   Cannot solve MIQP problems with HiGHS
WARNING:linopy.solvers:Solution status unknown. Trying to parse solution.
WARNING:linopy.constants:Optimization failed: 
Status: unknown
Termination condition: unknown
Solution: 3 primals, 1 duals
Objective: 0.00e+00
Solver model: available
Solver message: not set

Here is the result using "ipm" option with highs

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options:
 - solver: ipm
INFO:linopy.io: Writing time: 0.03s
Running HiGHS 1.9.0 (git hash: n/a): Copyright (c) 2024 HiGHS under MIT licence terms
INFO:linopy.solvers:Log file at /private/var/folders/c3/657p__892k72qxh6r8fdkqsh0000gn/T/highs.log
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e-01, 7e+00]
  Bound  [1e+00, 8e+01]
  RHS    [4e+00, 4e+00]
Solving LP relaxation since solver = ipm
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-1); columns 0(-3); elements 0(-3) - Reduced to empty
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-agbiaazn
Model status        : Optimal
Objective value     : -1.5153803770e+01
HiGHS run time      :          0.00
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 3 primals, 1 duals
Objective: -1.52e+01
Solver model: available
Solver message: optimal

Here is the result from cplex

INFO:linopy.model: Solve problem using Cplex solver
INFO:linopy.io: Writing time: 0.03s
Version identifier: 22.1.1.0 | 2022-11-28 | 9160aff4d
CPXPARAM_Read_DataCheck                          1
Found incumbent of value -9.153804 after 0.00 sec. (0.00 ticks)
Tried aggregator 1 time.
MIQP Presolve eliminated 1 rows and 3 columns.
All rows and columns eliminated.
Presolve time = 0.00 sec. (0.00 ticks)

Root node processing (before b&c):
  Real time             =    0.00 sec. (0.00 ticks)
Parallel b&c, 10 threads:
  Real time             =    0.00 sec. (0.00 ticks)
  Sync time (average)   =    0.00 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    0.00 sec. (0.00 ticks)
WARNING:linopy.solvers:Dual values of MILP couldn't be parsed
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 3 primals, 0 duals
Objective: -9.15e+00
Solver model: available
Solver message: integer optimal solution

And here is the result with scip

INFO:linopy.io: Writing time: 0.02s
original problem has 4 variables (0 bin, 3 int, 0 impl, 1 cont) and 2 constraints
presolving:
   (0.0s) symmetry computation started: requiring (bin +, int +, cont +), (fixed: bin -, int -, cont -)
   (0.0s) no symmetry present (symcode time: 0.00)
presolving (0 rounds: 0 fast, 0 medium, 0 exhaustive):
 0 deleted vars, 0 deleted constraints, 0 added constraints, 0 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
 0 implications, 0 cliques
presolved problem has 4 variables (2 bin, 1 int, 0 impl, 1 cont) and 2 constraints
      1 constraints of type <linear>
      1 constraints of type <nonlinear>
Presolving Time: 0.00

 time | node  | left  |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr|  dualbound   | primalbound  |  gap   | compl. 
t 0.0s|     1 |     0 |     0 |     - | trivial|   0 |   4 |   1 |   0 |   0 |  0 |   0 |   0 |-9.153804e+00 |-9.153804e+00 |   0.00%| unknown

SCIP Status        : problem is solved [optimal solution found]
Solving Time (sec) : 0.00
Solving Nodes      : 1
Primal Bound       : -9.15380377099400e+00 (2 solutions)
Dual Bound         : -9.15380377099400e+00
Gap                : 0.00 %
[cons_linear.c:18687] ERROR: constraint is not linear
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 3 primals, 0 duals
Objective: -9.15e+00
Solver model: available
Solver message: optimal

[GbotemiB]

I think it is ok to leave the "highs" out of the solver for now.

All good from my end.

[ekatef]

Perfect, thank you so much @GbotemiB!

Regarding HiGHS, it appears that there may be some deeper issues: PyPSA-Eur also is not using it for clustering. But totally agree that it has been worth trying.

[ekatef]

Merged! Many thanks @GbotemiB for taking care about that. We are moving forward to make the model really stable! 😄