Linopy transition

[ekatef]

Closes #494

Integrate linopy following https://github.com/PyPSA/pypsa-eur/pull/625

Checklist

• [x] 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.
• [x] 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]

Currently works locally with cbc. Points which need attention:

• [ ] there is a mismatch between batteries efficiency and links which has temporary resolved by using a single efficiency value, but needs additional investigation
• [x] what does mean replacement p with p_set in add_operational_reserve_margin_constraint()?
• [ ] the solver log has disappeared from the logger files, while output of solver logging into a console presents if running solving from VSCode but is skipped is using a plain terminal
• [ ] there is a fuzzy feeling that double-check is needed on transforming to xarray; in particular, the current calculation of ext_carrier_i differ from PyPSA-Eur implementation https://github.com/PyPSA/pypsa-eur/blob/6505d78f4e9a1362e8309bcd5b4b8b99ac60337c/scripts/solve_network.py#L350

[ekatef]

I can reproduce CI error locally when using lgpk, while the workflow works with cbc

[pz-max]

Cool! Just a comment. I would not follow my Linopy PR only but also check if the PyPSA-EUR main changed some parts.

[davide-f]

Hello @ekatef, that's great news! :D To answer to a question, the use of p_set is the more appropriate way to get access to the demand time series of load.

Regarding solver log, it may make sense to check the latest pypsa-eur implementation as I guess they have accounted for that.

As a remark, this morning we talked with Leon and it turns out that to merge this PR requires changes to the -sec version, which may be delicate also for some of their projects. So, it is a good idea to investigate and prepare this PR, but it can be merged only when the alignment can be performed also with the Sec version. What is your feeling?

[ekatef]

Hey @davide-f!

Thanks a lot for the detailed analysis :)

Ah, I see: the load is an input, not a result. So, n.loads_t.p_set is in fact more appropriate than n.loads_t.p. Thanks for the explanation!

Regarding a potential merge of this PR, absolutely agree that it's crucial to ensure that the changes introduced don't break anything. Actually, I perceive this work as a kind of experiment to estimate how much effort would be needed for transition to linopy interface. If you see a potential to align this with requirements of PyPSA-Sec at some stage, that is great and happy to adjust this PR accordingly.

My general impression is that linopy transition is not as challenging as it could be expected from the first sight (thanks to @pz-max PR in PyPSA-Eur and clear instructions)

[ekatef]

Some adjustments were done to get closer to PyPSA-Eur implementation. However, glpk solution still breaks. The reason seems to be the following command is given to a solver: 'glpsol --lp /var/folders/qn/vpndfm21795ckkq89np1ckp40000gn/T/linopy-problem-p4tiu2uu.lp --output /var/folders/qn/vpndfm21795ckkq89np1ckp40000gn/T/linopy-solve-7tkku64p.sol --solver_logfile /Users/ek/_github_/pypsa-earth/logs/NG_BJ/solve_network/elec_s_6_ec_lcopt_Co2L-4H_solver.log'

According to a warning, there is no --solver_logfile option available for glpsol, while it is transferred intorun_glpk as solver_options (not sure yet about the reason).

[ekatef]

Results of additional testing:

1. now highs is captured and used;
2. some introduced changes make a testing problem infeasible, and it looks like the reason for it is that the problem formulation somehow gets lost:

INFO:linopy.solvers:Log file at /tmp/highs.log.
ERROR:   getOptionIndex: Option "solver_logfile" is unknown
Presolving model
Problem status detected on presolve: Infeasible
Model   status      : Infeasible
Objective value     :  0.0000000000e+00
HiGHS run time      :          0.00
WARNING:linopy.constants:Optimization failed: 
Status: warning
Termination condition: infeasible
Solution: 0 primals, 0 duals
Objective: nan
Solver model: available
Solver message: infeasible

3. there are some troubles with logfile options for cbc:

No match for solver_logfile - ? for list of commands
No match for logs/MY/solve_network/elec_s_10_ec_lcopt_Co2L-2H_solver.log - ? for list of commands

[davide-f]

IT may be good to revive this PR and finalize it when possible, to be able to use the latest pypsa version. How far do you think we are to finalize?

I am aware that you have a lot on your desk right now, no pressure at all :) we can give priorities and tackle them one by one

[ekatef]

IT may be good to revive this PR and finalize it when possible, to be able to use the latest pypsa version. How far do you think we are to finalize?

I am aware that you have a lot on your desk right now, no pressure at all :) we can give priorities and tackle them one by one

@davide-f agree :) I'd expect that something about a week of focused work should be enough to resolve the remained questions. As discussed, #919 and #903 seem to have higher priorities now. Also, not sure what is the current status of PyPSA-Earth-Sec in terms of pyomo/linopy dependencies.

Anyway, after #919 and #903 will be finished, happy to get back to this PR. Ping me please in case this would get urgent.

[ekatef]

Testing with the updated PyPSA version (0.25.2 is used locally).

Clustering is failing currently:

1. the default SCIP solver (copy-pasted from PyPSA-Eur) is not included into PyPSA-Earth environment;
2. ipopt is installed but not seen by linopy for some reason.

CI is failing due to conda being unhappy: Operation too slow. Less than 30 bytes/sec transferred the last 60 seconds.

[davide-f]

Great work :) I see it is in work in progress; the CI is complying about the number of clusters; if/once you need a review/check, ping me :)

[ekatef]

Great work :) I see it is in work in progress; the CI is complying about the number of clusters; if/once you need a review/check, ping me :)

Thanks :D Yeah, there have been linopy-related changes in clustering, and some clustering parameters have been shifted. However, as we discussed, I'd focus on making the latest PyPSA version work first (hopefully, we are almost there :D), and adjusted the clustering details after.

Great to have you support and I'll definitely ask your assessment, once there will be some updates here! 🙂

[ekatef]

Update: The PR allows to run the workflow with the latest version of PyPSA (0.27.1), but gplk error persists:INFO:linopy.solvers:Invalid option '--solver_logfile'; try glpsol --help.

I suspect, that leads to troubles when linopy tries to read solver settings: the line info = pd.read_csv(info, sep=":", index_col=0, header=None)[1]) triggers pandas issues pandas.errors.EmptyDataError: No columns to parse from file.

It looks like some glpk setup should be adjusted.

[ekatef]

Testing with cbc also fails in some enigmatic way: optimisation itself is successful, while the workflow crashes after with TypeError: cannot unpack non-iterable NoneType object in status, condition = n.optimize.optimize_transmission_expansion_iteratively(**kwargs). Interestingly, the solver log is empty. Could it possibly mean some troubles when writing the solving results?

The full listing in terminal looks as follows:

rule solve_network:
    input: networks/NG_BJ/elec_s_8_ec_lcopt_Co2L-4H.nc
    output: results/NG_BJ/networks/elec_s_8_ec_lcopt_Co2L-4H.nc
    log: logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_solver.log, logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_python.log
    jobid: 1
    benchmark: benchmarks/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H
    reason: Missing output files: results/NG_BJ/networks/elec_s_8_ec_lcopt_Co2L-4H.nc; Input files updated by another job: networks/NG_BJ/elec_s_8_ec_lcopt_Co2L-4H.nc
    wildcards: simpl=, clusters=8, ll=copt, opts=Co2L-4H
    resources: tmpdir=/var/folders/qn/vpndfm21795ckkq89np1ckp40000gn/T, mem=8670
INFO:snakemake.logging:rule solve_network:
    input: networks/NG_BJ/elec_s_8_ec_lcopt_Co2L-4H.nc
    output: results/NG_BJ/networks/elec_s_8_ec_lcopt_Co2L-4H.nc
    log: logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_solver.log, logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_python.log
    jobid: 1
    benchmark: benchmarks/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H
    reason: Missing output files: results/NG_BJ/networks/elec_s_8_ec_lcopt_Co2L-4H.nc; Input files updated by another job: networks/NG_BJ/elec_s_8_ec_lcopt_Co2L-4H.nc
    wildcards: simpl=, clusters=8, ll=copt, opts=Co2L-4H
    resources: tmpdir=/var/folders/qn/vpndfm21795ckkq89np1ckp40000gn/T, mem=8670

INFO:snakemake.logging:
Restricted license - for non-production use only - expires 2025-11-24
~/miniconda3/envs/pypsa-earth-upg4/lib/python3.10/site-packages/pulp/apis/gurobi_api.py:170: UserWarning: GUROBI error: Unknown parameter '_test'.
  warnings.warn("GUROBI error: {}.".format(e))
Building DAG of jobs...
WARNING:snakemake.logging:Building DAG of jobs...
Using shell: /bin/bash
WARNING:snakemake.logging:Using shell: /bin/bash
Provided cores: 1 (use --cores to define parallelism)
WARNING:snakemake.logging:Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
WARNING:snakemake.logging:Rules claiming more threads will be scaled down.
Provided resources: mem=8670
WARNING:snakemake.logging:Provided resources: mem=8670
Select jobs to execute...
WARNING:snakemake.logging:Select jobs to execute...
Changing to shadow directory: _github_/pypsa-earth/.snakemake/shadow/tmp1q_7j4bm
WARNING:snakemake.logging:Changing to shadow directory: _github_/pypsa-earth/.snakemake/shadow/tmp1q_7j4bm
INFO:pypsa.io:Imported network elec_s_8_ec_lcopt_Co2L-4H.nc has buses, carriers, generators, global_constraints, links, loads, storage_units, stores
No expandable lines found. Skipping iterative solving.
INFO:__main__:No expandable lines found. Skipping iterative solving.
INFO:linopy.model: Solve problem using Cbc solver
INFO:linopy.model:Solver options:
 - solver_logfile: logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_solver.log
 - store_basis: True
INFO:linopy.io: Writing time: 0.13s
INFO:linopy.solvers:No match for solver_logfile - ? for list of commands
No match for logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_solver.log - ? for list of commands
No match for store_basis - ? for list of commands
No match for True - ? for list of commands
Total time (CPU seconds):       0.07   (Wallclock seconds):       0.08


INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 3779 primals, 7761 duals
Objective: 3.15e+09
Solver model: not available
Solver message: Optimal - objective value 3154337070.51124573


~/miniconda3/envs/pypsa-earth-upg4/lib/python3.10/site-packages/pypsa/optimization/optimize.py:357: FutureWarning: A value is trying to be set on a copy of a DataFrame or Series through chained assignment using an inplace method.
The behavior will change in pandas 3.0. This inplace method will never work because the intermediate object on which we are setting values always behaves as a copy.

For example, when doing 'df[col].method(value, inplace=True)', try using 'df.method({col: value}, inplace=True)' or df[col] = df[col].method(value) instead, to perform the operation inplace on the original object.


  n.df(c)[attr + "_opt"].update(df)
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.
INFO:pypsa.optimization.abstract:Mean square difference after iteration 1 is nan
INFO:pypsa.optimization.abstract:Running last lopf with fixed branches (HVDC links and HVAC lines)
INFO:linopy.model: Solve problem using Cbc solver
INFO:linopy.model:Solver options:
 - solver_logfile: logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_solver.log
 - store_basis: True
INFO:linopy.io: Writing time: 0.13s
INFO:linopy.solvers:No match for solver_logfile - ? for list of commands
No match for logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_solver.log - ? for list of commands
No match for store_basis - ? for list of commands
No match for True - ? for list of commands
Total time (CPU seconds):       0.07   (Wallclock seconds):       0.08


INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 3779 primals, 7761 duals
Objective: 3.15e+09
Solver model: not available
Solver message: Optimal - objective value 3154337070.51124573


~/miniconda3/envs/pypsa-earth-upg4/lib/python3.10/site-packages/pypsa/optimization/optimize.py:357: FutureWarning: A value is trying to be set on a copy of a DataFrame or Series through chained assignment using an inplace method.
The behavior will change in pandas 3.0. This inplace method will never work because the intermediate object on which we are setting values always behaves as a copy.

For example, when doing 'df[col].method(value, inplace=True)', try using 'df.method({col: value}, inplace=True)' or df[col] = df[col].method(value) instead, to perform the operation inplace on the original object.


  n.df(c)[attr + "_opt"].update(df)
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.
ERROR:_helpers:An error happened in module '_github_/pypsa-earth/scripts/solve_network.py', function 'solve_network': cannot unpack non-iterable NoneType object
Traceback (most recent call last):
  File "_github_/pypsa-earth/.snakemake/shadow/tmp1q_7j4bm/.snakemake/scripts/tmpidm9hrsp.solve_network.py", line 596, in <module>
    n = solve_network(
  File "_github_/pypsa-earth/.snakemake/shadow/tmp1q_7j4bm/.snakemake/scripts/tmpidm9hrsp.solve_network.py", line 544, in solve_network
    status, condition = n.optimize.optimize_transmission_expansion_iteratively(**kwargs)
TypeError: cannot unpack non-iterable NoneType object
Not cleaning up _github_/pypsa-earth/.snakemake/shadow/tmp1q_7j4bm/.snakemake/scripts/tmpidm9hrsp.solve_network.py
CRITICAL:snakemake.logging:Not cleaning up _github_/pypsa-earth/.snakemake/shadow/tmp1q_7j4bm/.snakemake/scripts/tmpidm9hrsp.solve_network.py
[Wed Apr 10 08:32:22 2024]
INFO:snakemake.logging:[Wed Apr 10 08:32:22 2024]
Error in rule solve_network:
    jobid: 0
    input: networks/NG_BJ/elec_s_8_ec_lcopt_Co2L-4H.nc
    output: results/NG_BJ/networks/elec_s_8_ec_lcopt_Co2L-4H.nc
    log: logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_solver.log, logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_python.log (check log file(s) for error details)

ERROR:snakemake.logging:Error in rule solve_network:
    jobid: 0
    input: networks/NG_BJ/elec_s_8_ec_lcopt_Co2L-4H.nc
    output: results/NG_BJ/networks/elec_s_8_ec_lcopt_Co2L-4H.nc
    log: logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_solver.log, logs/NG_BJ/solve_network/elec_s_8_ec_lcopt_Co2L-4H_python.log (check log file(s) for error details)

RuleException:
CalledProcessError in file _github_/pypsa-earth/Snakefile, line 837:
Command 'set -euo pipefail;  ~/miniconda3/envs/pypsa-earth-upg4/bin/python3.10 _github_/pypsa-earth/.snakemake/shadow/tmp1q_7j4bm/.snakemake/scripts/tmpidm9hrsp.solve_network.py' returned non-zero exit status 1.
  File "_github_/pypsa-earth/Snakefile", line 837, in __rule_solve_network
  File "~/miniconda3/envs/pypsa-earth-upg4/lib/python3.10/concurrent/futures/thread.py", line 58, in run
ERROR:snakemake.logging:RuleException:
CalledProcessError in file _github_/pypsa-earth/Snakefile, line 837:
Command 'set -euo pipefail;  ~/miniconda3/envs/pypsa-earth-upg4/bin/python3.10 _github_/pypsa-earth/.snakemake/shadow/tmp1q_7j4bm/.snakemake/scripts/tmpidm9hrsp.solve_network.py' returned non-zero exit status 1.
  File "_github_/pypsa-earth/Snakefile", line 837, in __rule_solve_network
  File "~/miniconda3/envs/pypsa-earth-upg4/lib/python3.10/concurrent/futures/thread.py", line 58, in run
Shutting down, this might take some time.
WARNING:snakemake.logging:Shutting down, this might take some time.
Exiting because a job execution failed. Look above for error message
ERROR:snakemake.logging:Exiting because a job execution failed. Look above for error message
Shutting down, this might take some time.

[ekatef]

Great work :) I see it is in work in progress; the CI is complying about the number of clusters; if/once you need a review/check, ping me :)

@davide-f When you have time, it would be great to have your opinion on what could go wrong with solver set-up for linopy. I have tested glpk and cbc with the latest PyPSA version of 0.27.1, but both fail due to some troubles which look like some misunderstanding between a solver or linopy and/or the system itself.

Any ideas on possible workarounds, would be very helpful!

[davide-f]

Ah @ekatef , we have merge conflicts; once they are solved, the CI is executed again

[pz-max]

Just for historical updates as I observed the PyPSA-Eur merge.

• we don't have to introduce Linopy here now and can do it in the merge (PyPSA-Earth-Sec would break if we introduce linopy here - but I am sure you tracking this)
• we can use the exact PyPSA-Eur equivalents, saving a lot of work

((so wouldn't spend too much time here))

[davide-f]

@ekatef It works now!!!! :D

[ekatef]

@ekatef It works now!!!! :D

Fantastic news!!! :)))) Still struggling with local testing (probably, an environment update is needed...), but trust to you and CI :D

[ekatef]

@ekatef It works now!!!! :D

@davide-f I confirm that it works!! 🎉 🎉 🎉 The trick has been to update the environment 😄

Thanks a lot for the final fixes!

Pre-commit is not happy as yaml-linter complains on the lines wider that 88. But I'm not sure it makes sense to fix it, as it is not linked with the current PR and would require to completelly re-structure the config.

[ekatef]

As a nice side-effect: PyPSA-Earth now works with HiGHS 🥳

[davide-f]

Amazing!

[ekatef]

The major points left:

• [x] do not revise implementation of add_operational_reserve_margin (in PyPSA-Earth decomposed on two functions which look a way more clear as compared with PyPSA-Eur);
• [x] check implementation of add_battery_constraints. Currently, it's different in PyPSA-Earth and PyPSA-Eur, which may it be the effect of the sec-merge (which imply the need to revise implementation after PyPSA-Earth merge);
• [x] add_RES_constraints must be refactored and reimplemented;
• [x] solve_network should be checked: that may be the point where PyPSA-Earth may need a slightly different approach as one used in PyPSA-Eur may be not the most clear one from the usability perspective.

[ekatef]

Update: testing successful in the base-case configuration, but there is a problem when fetching isolated networks (something with p key which is expected but not found in the provided data structures).

[ekatef]

A couple of notes on linopy grammar:

1. pypsa.Network() doesn't initialise model property -> an explicit call n.optimize.create_model() is needed
2. Multiplication is not commutative in the current linopy implementation: n.model["Generator-p"] * n.snapshot_weightings.generators works, while n.snapshot_weightings.generators * n.model["Generator-p"]` throws an error.
3. Not sure what is the way to make groupby(some_grouper).sum() work on linopy variables.

[ekatef]

Local testing on reproducibility of objective value

main version

INFO:main:Objective function: 5813422825.0 INFO:main:Objective constant: 2319283737.695124

linopy version

INFO:main:Objective function: 5473150322.268269 INFO:main:Objective constant: 2319359288.9025226

That means ~6% difference for the objective function and basically same value for the objective constant (1e-3%, but here is rather a question of numerical precision). So, the PR seems to work properly.

As a usability note, that is not obvious at all from the CI logs. The PR changes format of solver traceback due to different conventions of lopf and linopy, which makes finding correspondance quite tricky.

[ekatef]

@finozzifa would be very grateful if you would have time to check this PR.

Pre-commit is currently unhappy due to excessively long lines in the configuration file. But that won't be fixed in this PR.

Another concern is backward compatibility. I'm currently looking into that, but not sure if it'll be possible to find a code solution for that.

[ekatef]

As discussed with @davide-f, the backward compatibility issues can be tackled by maintenance measures and adding versioning of the configuration files. A possible implementation for versioning #1058.

This PR is ready for review.

[ekatef]

Added some first comments :) Great work. Indeed this PR is quite large too already, let's find the best way to finalize it :) Let's remember to add a line in the release_note too

@davide-f thanks a lot for the great review! You have made me recognise that PR mixes two different things: PyPSA update and actually linopy transition. Moreover, the updated PyPSA version leads to some odd behaviour in the case when custom columns are present in the components on the network, which should be tracked regardless linopy methods.

Have opened #1065 and moved there the changes from this PR which relate to PyPSA update. Would you mind to continue the discussion on the aggregation strategies there? Then, for this PR we can focus on the solving part only. What do you think?

[ekatef]

Thanks a lot for the in-depth review @davide-f!

Along the way, we have identified a number of the additional issues, both related to this PR and completely independent ones. The list looks like follows

1. Fetching isolated networks gets broken at the stage of applying get_switchable_as_dense( ) when aggregated storage units and loads. The issue looks quite weird and I have not been able to debug it properly. Currently, it's "fixed" by switching-off the functionality to merge isolated grids. That is not a solution which which I'm happy but it could make sense to tackle it separately.
2. add_RES_constraints( ) may need further revisions.
3. A revision seems to be needed for implementation of the hydro-related part in add_EQ_constraints( ) [both in PyPSA-Earth and PyPSA-Eur]
4. It may be a good idea to adjust a definition of the reserve margin constraint in add_operational_reserve_margin_constraint( ) to make it compatible with the current regulation status in the most countries of the world.

Regarding the overall implementation strategy, my feeling is that it makes sense to disentangle PyPSA upgrade from enabling linopy. Have opened #1065 to update PyPSA and a stacked PR to enable linopy. Both PRs are mainly cherry-picking from this PR which should hopefully facilitate the final checks.

The PyPSA update leads to ~1% change of the objective function, while enabling linopy leaves the objective function pretty much the same (~1e-5 is the relative difference). So, I think this work is close to be finished, keeping in mind the points above 🙂