PyPSA: Python for Power System Analysis

Introduce optimization without usage of pyomo, based on nomopyomo which is way faster and memory efficient for especially large networks.

What happened so far:

New Modules

• linopf.py - includes all code for setting up the problem
• linopt.py - includes all solver (io to solver) relevant code
• stats.py - include code for post-optimization analysis

Ordinary features of lopf without pyomo:

• [x] include ramp limits
• [x] include Storage Units
• [x] calculate nodal imbalances
• [x] constant term in objective function
• [x] calculate voltage angles
• [x] implement glpk solver
• [x] logfile for cbc

Additional features

• [x] global constraints for transmission expansion volume/cost limit
• [x] ability to fix values via '_set' for arbitrary variables, also static variables (so f.e. n.generators.p_nom_set can fix values for extendable p_nom variables)
• [x] warm start
• [x] add iterative lopf (ilopf), updating impedances after each solving process
• [x] docstrings for linop{t,f}.py
• [x] add docstrings to documentation (+ incorporate changes)
• [x] ensure windows compatibility
• [x] add windows and osx in travis CI

Other changes

• determine_network_topology now includes find_slackbus and find_cycles

• add new lopf to tests in

  1. acdc test
  2. storage unit test

• skip python2.7 test in travis file. I had to much python3 code syntax, and didn't see the nessecity to rewrite it as so many package give up on python2 support. But perhaps you have a different view.

The network.lopf function got a new argument as proposed in #96, so it becomes possible now to call

network.lopf(pyomo=False)

In this pull request the licence for PyPSA is changed from the copyleft GNU General Public License Version 3 (GPLv3) to the more liberal MIT License. In addition, all copyright notices refer to the new list of PyPSA Developers in the documentation.

This change is suggested by the main PyPSA developers because the benefits of a more liberal license seem to outweigh the benefits of a copyleft approach.

Benefits of copyleft for PyPSA: The copyleft GPLv3 requires users who release software binaries based on PyPSA, even when used as a library, to also release their source code. This would benefit PyPSA if for example somebody developed a big extension, like a GUI, then they would have to release it under the GPLv3. However, we know no cases where anyone is planning such a large extension to sell commercially. Most user developments are either contributed back or kept private, rather than being sold in software bundles. If someone were to develop a commercial product on top of PyPSA, they would mostly likely rather choose liberal-licensed calliope or oemof. So the copyleft offers few advantages.

Benefits of liberal license: We have heard of several users avoiding PyPSA because of its copyleft GPLv3 license, which forces them to release their code that depends on PyPSA under the GPLv3. We would like to maximise the usage and impact of PyPSA to accelerate the energy transition.

The original justification for choosing GPLv3 was that it allowed us the freedom to change the license later. It makes more sense to step from a copyleft to a liberal license than the other way around.

We also refer to the very informative discussion on the switch from GPLv3 to MIT license in the oemof developer community.

This change requires the consent of all people who have contributed code to PyPSA. If you consent, please add a :heavy_check_mark: after your name or respond below. If you do not consent, please add a :x: after your name or respond below. Please respond by 1st August 2021. Please feel free to also discuss below. Since the major developers have agreed, we will either remove or rewrite contributions from those who do not consent. Thanks for your help!

• @nworbmot :heavy_check_mark:
• @coroa :heavy_check_mark:
• @FabianHofmann :heavy_check_mark:
• @fneum :heavy_check_mark:
• @lisazeyen :heavy_check_mark:
• @martacki :heavy_check_mark:
• @bourbl :heavy_check_mark:
• @schlachtberger :heavy_check_mark:
• @jdedecca :heavy_check_mark:
• @stianchris :heavy_check_mark:
• @fonsp
• @dcradu :heavy_check_mark:
• @hailiangliu89 :heavy_check_mark:
• @koen-vg :heavy_check_mark:
• @BurningTrapezoid :heavy_check_mark:
• @kaiweiy
• @nmartensen :heavy_check_mark:
• @rockstaedt :heavy_check_mark:
• @ebbekyhl :heavy_check_mark:
• @MinhHaDuong :heavy_check_mark:
• @mariiabilousova :heavy_check_mark:
• @chiaroo :heavy_check_mark:
• @jankaeh :heavy_check_mark:
• @miek770 :heavy_check_mark:
• @merraksh :heavy_check_mark:
• @euronion :heavy_check_mark:

Checklist

• [x] Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
• [x] A note for the release notes doc/release_notes.rst of the upcoming release is included.

As of Version 13.1, PyPSA only considers investment in one horizon, e.g. for a single year. We'd like to implement investment over multiple horizons simultaneously, e.g. for 2020, 2025, 2030, 2035, etc., to optimise investments dynamically over several years e.g. under tightening CO2 constraints. The objective function then becomes the discounted net present value (NPV) of energy system costs. Other open models such as OSeMOSYS, Genesys, DIETER (I think) and SWITCH do this already; we can copy their implementations and avoid reinventing the wheel.

On the data side:

i) We need a second index for snapshots that gives the investment year, e.g. snapshots = [(2020, "windy day"), (2020, "sunny day"), (2030, "windy day"), (2030, "sunny day")] The representative weather for each investment year may change with the year (e.g. due to climate change); the number of representative snapshots for each investment year may also change.

ii) We need different CO2/renewable constraints for different investment years.

iii) For assets, we need to take account of build years and lifetimes for generation, storage and transmission assets, e.g. we need additional fields

name, build year, life time, cost wind-in-2020, 2020, 20, 1200 wind-in-2025, 2025, 20, 1000 wind-in-2030, 2030, 25, 900

Then in the code we need to activate assets in the constraints depending on whether they are around in the investment year or not. We also need to take account of the discount rate as a pre-factor to get the NPV for each investment year. Perhaps this should be linked to issue #37 of separating the snapshot_weighting from the duration of each snapshot, which was discussed on the PyPSA mailing list before. BEWARE 1): The final investment horizon has to be treated carefully, to avoid end-of-horizon effects whereby assets are only built for the final horizon and have no value for the rest of their lifetimes. Either model 20 years further (with lower intra-year resolution) or have a "salvage value" for investments made in the final investment horizon (I think OSeMOSYS had this). BEWARE 2): p_nom_max may become linked between different generators as potentials are used up.

Changes proposed in this Pull Request

A new sub-module "optimization" for the linopy-based optimization is introduced. The optimization module is sub-structured in

• variables.py - for all variable definitions
• constraints.py - for all local constraint definitions
• global_constraints.py - for all global constraints
• common.py - for common utility functions
• optimize.py - for functions running the sanitiy checks, optimization and solution assignment
• abstract.py - for optimization functions abstracted from the normal lopf, (like iterative lopf, possibly also sclopf and others)

Remaining Todo's

• [x] sanitize check
  • [X] nans in nominal bounds
  • [x] nans in operational bounds
  • [x] ~link efficiency < 1 #319~
  • [x] intersection of committable and extendable assets
• [x] support multi-port links
• [x] support min_up_time min_down_time constraints
• [x] make p_set argument available for all static and time-dependent variables #202
• [x] documentation

Run it

Run the optimization by

n.optimize()

you can also create the model and optimize afterwards

n.optimize.create_model()
# do whatever stuff
n.optimize.solve_model()

You can access variables in

n.model.variables

and constraints in

n.model.constraints

From there you can look up and modify the LP, like for example

n.model.constraints['my_constraint'].lhs = n.model.constraints['my_constraint'].lhs + another_variable

Happy to have some feedback @nworbmot @fneum @lisazeyen @pz-max @p-glaum

Preview Documentation

https://pypsa--336.org.readthedocs.build/en/336/

Checklist

• [x] Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
• [x] Unit tests for new features were added.
• [x] Newly introduced dependencies are added to environment.yaml, environment_docs.yaml and setup.py (if applicable).
• [x] A note for the release notes doc/release_notes.rst of the upcoming release is included.
• [x] I consent to the release of this PR's code under the MIT license.

Closes # (if applicable).

Changes proposed in this Pull Request

Add open source HiGHS solver:

Further performance check:

• 1year, 3h resolution, 3 nodes: HiGHS 8min, Gurobi ~1min
• 1year, 1h resolution, 3 nodes: HiGHS 50min, Gurobi ~1min

Future work can:

• While the above results look incredible, generally all open source solvers struggle when 1 year is simulated in hourly resolution. I.e. Gurobi solves 5 nodes for 1 year in hourly resolution in 1min, while HiGHS and others may need 2+ hours. This performance drop occurs probably because open-source solvers don't address properly the existence of "dense columns".
• Julian Hall, corrects this statement by saying the HiGHS IPM deals with dense columns. Bottlenecks for the performance drops need to be investigated. Further, Cbc is the MIP solver, one might be better referring to Clp when only solving LPs.
• Other work could implement problem specific tailored scaling techniques into PyPSA (to be useful for all solvers). See some ongoing discussion here: https://github.com/PyPSA/PyPSA/issues/309

More discussions: https://forum.openmod.org/t/open-source-highs-solver-performance-boost-for-energy-system-models/2922/2

Checklist

• [x] Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
• [x] A note for the release notes doc/release_notes.rst of the upcoming release is included.
• [x] I consent to the release of this PR's code under the MIT license.

@FabianHofmann and I discussed a few ideas for improving the PyPSA documentation:

• [x] add autodocs for documented functions
• [x] change theme to sphinx_rtd_theme from alabaster
• [x] readthedocs automatic deployment
• [x] remove some duplicate content
• [x] add shields.io stickers
• [x] API reference

What one could do another time:

• update screenshots with new colors (https://github.com/PyPSA/pypsa-eur/pull/35) and cartopy
• sort examples

Current state is on https://pypsa-readthedocs.readthedocs.io/

Once merged, need to change settings on https://readthedocs.org and add correct webhook in github as well as forwarding https://pypsa.org/documentation to https://pypsa.readthedocs.io.

The last linopf iteration fixes the link and line capacities to ensure consistency of impedances and power flow. However, when using store components, it also fixes the capacities of their respective in- and outflow links. Filtering by carrier solves this problem. This is a quick fix which has worked well for me, I'm sure there's a more elegant solution than using a for loop - however, it is nice that no error messages are displayed, using the loc approach.

Closes # (if applicable).

Changes proposed in this Pull Request

Filtering by carrier before fixing p_nom links

Checklist

• [x] Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
• ~[ ] Unit tests for new features were added (if applicable).~
• ~[ ] Newly introduced dependencies are added to environment.yaml, environment_docs.yaml and setup.py (if applicable).~
• [x] A note for the release notes doc/release_notes.rst of the upcoming release is included.

I am solving a lopf using the Barrier method from Gurobi (with the crossover turned off, since it takes too long otherwise) and I get just a sub-optimal solution. Then PyPSA prints:

WARNING:pypsa.opf:WARNING! Optimization might be sub-optimal. Writing output anyway

but it does not write any output. It seems to be problem with the format that the sub-optimal solution is given (maybe having some differences with the way an optimal solution is given).

Side note: I am using "extendable" storage and I am mainly interested in the "p_nom_opt" that the solver would give me.

Changes proposed in this Pull Request

The pandapower network importer was not working correctly. @rbolgaryn and myself working currently on fixing the importer. In particular, we are fixing the transformer and the shunt components.

Checklist

• [x] Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
• [x] Unit tests for new features were added (if applicable).
• [x] A note for the release notes doc/release_notes.rst of the upcoming release is included.
• [x] I consent to the release of this PR's code under the MIT license.

Checklist

• [x] I am using the current master branch or the latest release. Please indicate:
• [x] test whether potential bug persists for more than the one snapshot @camiloromers.

Describe the Bug

When I try to run a OPF using the network folder attached below I got the following error:

ERROR: Rule failed when generating expression for constraint
    gen_up_time_force_0 with index 1: IndexError: index 1 is out of bounds for
    axis 0 with size 1
ERROR: Constructing component 'gen_up_time_force_0' from data=None failed:
    IndexError: index 1 is out of bounds for axis 0 with size 1

May I know what's missing in my network or any clue relating the error I am getting?

Thanks for your help.

[L2RPN_2020_ieee_118_pypsa.zip](https://github.com/PyPSA/PyPSA/files/4071575/L2RPN_2020_ieee_118_pypsa.zip)

We have been using PyPSA, and we noticed that for larger numbers of snapshots and branches, the time to create the kirchhoff cycles was quite high. This PR significantly reduces the time, primarily by moving the time-independent attribute calculations outside of the snapshots loop.

I have also added a new example build_model.py, which can show the speed increase. It simply creates a (randomized) larger network to induce cycles. On my computer, the new build_model script runs about 50% faster on this branch than when I run it versus master.

Let me know if you would like any changes. One easy additional change could be to use define_sub_network_cycle_constraints for the cycles formulation as well, but I haven't been using it, so I haven't been able to check for sure that that would be reasonable.

Currently the flow argument of n.plot() hijacks the line_widths argument such that it no longer controls the line width but the arrow size. I would like to control arrow size and line widths separately.

Checklist

• [x] I am using the current master branch or the latest release. Please indicate:

Describe the Bug

This problem arises in the function get_buses_linemap_and_lines in the following code when the lines dataframe is empty:

return (
        buses,
        linemap,
        linemap_p,
        linemap_n,
        lines.reset_index()
        .rename(columns={"bus0_s": "bus0", "bus1_s": "bus1"}, copy=False)
        .set_index("name"),
        lines_t,
    )

This seems to be due to the function aggregatelines called previously that does stuff on the lines dataframe including groupbys. That should normally drop such columns but that does not happen when the dataframe is empty. This is an error we experienced when running the pypsa-earth workflow on rwanda.

Error Message

If applicable, paste any terminal output to help illustrating your problem. In some cases it may also be useful to share your list of installed packages: conda list.

Exception has occurred: ValueError       (note: full exception trace is shown but execution is paused at: _run_module_as_main)
cannot insert bus1_s, already exists
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.10/site-packages/pandas/core/frame.py", line 4814, in insert
    raise ValueError(f"cannot insert {column}, already exists")
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.10/site-packages/pandas/core/frame.py", line 6358, in reset_index
    new_obj.insert(
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.10/site-packages/pandas/util/_decorators.py", line 331, in wrapper
    return func(*args, **kwargs)
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.10/site-packages/pypsa/networkclustering.py", line 351, in get_buses_linemap_and_lines
    lines.reset_index()
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.10/site-packages/pypsa/networkclustering.py", line 378, in get_clustering_from_busmap
    buses, linemap, linemap_p, linemap_n, lines, lines_t = get_buses_linemap_and_lines(
  File "/data/davidef/git_world/pypsa-earth/scripts/cluster_network.py", line 444, in clustering_for_n_clusters
    clustering = get_clustering_from_busmap(
  File "/data/davidef/git_world/pypsa-earth/scripts/cluster_network.py", line 577, in <module>
    clustering = clustering_for_n_clusters(
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/home/davidef/miniconda3/envs/pypsa-earth/lib/python3.10/runpy.py", line 196, in _run_module_as_main (Current frame)
    return _run_code(code, main_globals, None,

Provide an easy way to calculate the carbon intensity of components after optimisation.

Common units: g/kWh, kg/MWh, t/MWh

For generators:

n.generators_t.p / n.generators.efficiency * n.generators.carrier.map(n.carriers.co2_emissions)

Needs to be generalised for all components.

Hi,

Occasionally we had the situation that there's subnetworks that could be defined separately, maybe we are also interested in solving those subnetwork separatly. Then we want to combine them into one big network. Would it be an idea to implement the class method add such that network1+network2 makes sense and gives you the combined network? Could make some workflows quite tidy. I guess there's a ton of things to consider there, overlapping buses for instance to have them connected after the merger. So at this stage just an idea.

Best, Andreas

Hi,

We found a small bug, pypsa 0.21.1:

#%%
import pypsa

n = pypsa.Network()
n.set_snapshots(np.arange(10))
n.add('Generator', name='test', bus='nobus', p_nom=100)
n.add('Bus', name='bus')
n.add('Generator', name='test2', bus='bus', p_nom=100)
n.lopf(solver_name='glpk')

#%%
from pypsa.linopt import define_constraints, get_var, linexpr
def extra(network, snapshots):
    define_constraints(network, linexpr((1, get_var(network, 'Generator', 'p').loc[:, 'test'])), '>=', 40, 'test')

n.lopf(solver_name='glpk', extra_functionality=extra)

A generator is put at a bus that doesn't exist - I guess a frequent thing in models due to data issues, and certainly captured by consistency_check. First stage solves normally and you could assume it just creates a virtual bus, equalizing power to zero there. It might be good if that was a bit stricter there...

Second one is a bit more puzzling, you can force the generator to produce power into the void. It solves normally. Not the most frequent scenario, but maybe worth a look.

Best, Andreas