diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index c2be3909..c0fb745d 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -83,6 +83,7 @@ jobs: snakemake -call solve_elec_networks --configfile config/test/config.electricity.yaml --rerun-triggers=mtime snakemake -call all --configfile config/test/config.overnight.yaml --rerun-triggers=mtime snakemake -call all --configfile config/test/config.myopic.yaml --rerun-triggers=mtime + snakemake -call all --configfile config/test/config.perfect.yaml --rerun-triggers=mtime - name: Upload artifacts uses: actions/upload-artifact@v3 diff --git a/.gitignore b/.gitignore index 459ef860..b79cd7ab 100644 --- a/.gitignore +++ b/.gitignore @@ -8,6 +8,7 @@ __pycache__ *dconf gurobi.log .vscode +*.orig /bak /resources diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 803eb61f..94a332b5 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -33,7 +33,7 @@ repos: rev: v2.2.6 hooks: - id: codespell - args: ['--ignore-regex="(\b[A-Z]+\b)"', '--ignore-words-list=fom,appartment,bage,ore,setis,tabacco,berfore'] # Ignore capital case words, e.g. country codes + args: ['--ignore-regex="(\b[A-Z]+\b)"', '--ignore-words-list=fom,appartment,bage,ore,setis,tabacco,berfore,vor'] # Ignore capital case words, e.g. country codes types_or: [python, rst, markdown] files: ^(scripts|doc)/ diff --git a/.readthedocs.yml b/.readthedocs.yml index 900dba1e..30684052 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -14,4 +14,3 @@ build: python: install: - requirements: doc/requirements.txt - system_packages: false diff --git a/Snakefile b/Snakefile index 32a1b125..83530df7 100644 --- a/Snakefile +++ b/Snakefile @@ -66,6 +66,11 @@ if config["foresight"] == "myopic": include: "rules/solve_myopic.smk" +if config["foresight"] == "perfect": + + include: "rules/solve_perfect.smk" + + rule all: input: RESULTS + "graphs/costs.pdf", diff --git a/config/config.default.yaml b/config/config.default.yaml index 81a26a0b..63f8f9c7 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -461,6 +461,7 @@ sector: years_of_storage: 25 co2_sequestration_potential: 200 co2_sequestration_cost: 10 + co2_sequestration_lifetime: 50 co2_spatial: false co2network: false cc_fraction: 0.9 @@ -501,6 +502,20 @@ sector: OCGT: gas biomass_to_liquid: false biosng: false + limit_max_growth: + enable: false + # allowing 30% larger than max historic growth + factor: 1.3 + max_growth: # unit GW + onwind: 16 # onshore max grow so far 16 GW in Europe https://www.iea.org/reports/renewables-2020/wind + solar: 28 # solar max grow so far 28 GW in Europe https://www.iea.org/reports/renewables-2020/solar-pv + offwind-ac: 35 # offshore max grow so far 3.5 GW in Europe https://windeurope.org/about-wind/statistics/offshore/european-offshore-wind-industry-key-trends-statistics-2019/ + offwind-dc: 35 + max_relative_growth: + onwind: 3 + solar: 3 + offwind-ac: 3 + offwind-dc: 3 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#industry industry: @@ -553,11 +568,13 @@ industry: hotmaps_locate_missing: false reference_year: 2015 + # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#costs costs: year: 2030 version: v0.6.0 rooftop_share: 0.14 # based on the potentials, assuming (0.1 kW/m2 and 10 m2/person) + social_discountrate: 0.02 fill_values: FOM: 0 VOM: 0 @@ -771,6 +788,7 @@ plotting: gas pipeline new: '#a87c62' # oil oil: '#c9c9c9' + imported oil: '#a3a3a3' oil boiler: '#adadad' residential rural oil boiler: '#a9a9a9' services rural oil boiler: '#a5a5a5' @@ -902,6 +920,7 @@ plotting: H2 for shipping: "#ebaee0" H2: '#bf13a0' hydrogen: '#bf13a0' + retrofitted H2 boiler: '#e5a0d9' SMR: '#870c71' SMR CC: '#4f1745' H2 liquefaction: '#d647bd' diff --git a/config/config.perfect.yaml b/config/config.perfect.yaml new file mode 100644 index 00000000..f355763c --- /dev/null +++ b/config/config.perfect.yaml @@ -0,0 +1,43 @@ +# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: CC0-1.0 +run: + name: "perfect" + +# docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#foresight +foresight: perfect + +# docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#scenario +# Wildcard docs in https://pypsa-eur.readthedocs.io/en/latest/wildcards.html +scenario: + simpl: + - '' + ll: + - v1.0 + clusters: + - 37 + opts: + - '' + sector_opts: + - 1p5-4380H-T-H-B-I-A-solar+p3-dist1 + - 1p7-4380H-T-H-B-I-A-solar+p3-dist1 + - 2p0-4380H-T-H-B-I-A-solar+p3-dist1 + planning_horizons: + - 2020 + - 2030 + - 2040 + - 2050 + + +# docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#co2-budget +co2_budget: + # update of IPCC 6th AR compared to the 1.5SR. (discussed here: https://twitter.com/JoeriRogelj/status/1424743828339167233) + 1p5: 34.2 # 25.7 # Budget in Gt CO2 for 1.5 for Europe, global 420 Gt, assuming per capita share + 1p6: 43.259666 # 35 # Budget in Gt CO2 for 1.6 for Europe, global 580 Gt + 1p7: 51.4 # 45 # Budget in Gt CO2 for 1.7 for Europe, global 800 Gt + 2p0: 69.778 # 73.9 # Budget in Gt CO2 for 2 for Europe, global 1170 Gt + + +sector: + min_part_load_fischer_tropsch: 0 + min_part_load_methanolisation: 0 diff --git a/config/test/config.perfect.yaml b/config/test/config.perfect.yaml new file mode 100644 index 00000000..c99f4122 --- /dev/null +++ b/config/test/config.perfect.yaml @@ -0,0 +1,91 @@ +# SPDX-FileCopyrightText: : 2017-2023 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: CC0-1.0 + +tutorial: true + +run: + name: "test-sector-perfect" + disable_progressbar: true + shared_resources: true + shared_cutouts: true + +foresight: perfect + +scenario: + ll: + - v1.0 + clusters: + - 5 + sector_opts: + - 8760H-T-H-B-I-A-solar+p3-dist1 + planning_horizons: + - 2030 + - 2040 + - 2050 + +countries: ['BE'] + +snapshots: + start: "2013-03-01" + end: "2013-03-08" + +electricity: + co2limit: 100.e+6 + + extendable_carriers: + Generator: [OCGT] + StorageUnit: [battery] + Store: [H2] + Link: [H2 pipeline] + + renewable_carriers: [solar, onwind, offwind-ac, offwind-dc] + +sector: + min_part_load_fischer_tropsch: 0 + min_part_load_methanolisation: 0 +atlite: + default_cutout: be-03-2013-era5 + cutouts: + be-03-2013-era5: + module: era5 + x: [4., 15.] + y: [46., 56.] + time: ["2013-03-01", "2013-03-08"] + +renewable: + onwind: + cutout: be-03-2013-era5 + offwind-ac: + cutout: be-03-2013-era5 + max_depth: false + offwind-dc: + cutout: be-03-2013-era5 + max_depth: false + solar: + cutout: be-03-2013-era5 + +industry: + St_primary_fraction: + 2020: 0.8 + 2030: 0.6 + 2040: 0.5 + 2050: 0.4 + +solving: + solver: + name: glpk + options: glpk-default + mem: 4000 + +plotting: + map: + boundaries: + eu_node_location: + x: -5.5 + y: 46. + costs_max: 1000 + costs_threshold: 0.0000001 + energy_max: + energy_min: + energy_threshold: 0.000001 diff --git a/doc/foresight.rst b/doc/foresight.rst index c1be3443..dd1e0ecc 100644 --- a/doc/foresight.rst +++ b/doc/foresight.rst @@ -41,10 +41,10 @@ Perfect foresight scenarios .. warning:: - Perfect foresight is currently under development and not yet implemented. + Perfect foresight is currently implemented as a first test version. -For running perfect foresight scenarios, in future versions you will be able to -set in the ``config/config.yaml``: +For running perfect foresight scenarios, you can adjust the + ``config/config.perfect.yaml``: .. code:: yaml diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 211cde63..505c747e 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -20,6 +20,8 @@ Upcoming Release * Files extracted from sector-coupled data bundle have been moved from ``data/`` to ``data/sector-bundle``. +* New feature multi-decade optimisation with perfect foresight. + * It is now possible to specify years for biomass potentials which do not exist in the JRC-ENSPRESO database, e.g. 2037. These are linearly interpolated. @@ -32,6 +34,7 @@ Upcoming Release * A bug preventing custom powerplants specified in ``data/custom_powerplants.csv`` was fixed. (https://github.com/PyPSA/pypsa-eur/pull/732) + PyPSA-Eur 0.8.1 (27th July 2023) ================================ diff --git a/doc/tutorial.rst b/doc/tutorial.rst index f0ded3fb..40c1907c 100644 --- a/doc/tutorial.rst +++ b/doc/tutorial.rst @@ -133,89 +133,82 @@ This triggers a workflow of multiple preceding jobs that depend on each rule's i graph[bgcolor=white, margin=0]; node[shape=box, style=rounded, fontname=sans, fontsize=10, penwidth=2]; edge[penwidth=2, color=grey]; - 0[label = "solve_network", color = "0.21 0.6 0.85", style="rounded"]; - 1[label = "prepare_network\nll: copt\nopts: Co2L-24H", color = "0.02 0.6 0.85", style="rounded"]; - 2[label = "add_extra_components", color = "0.37 0.6 0.85", style="rounded"]; - 3[label = "cluster_network\nclusters: 6", color = "0.39 0.6 0.85", style="rounded"]; - 4[label = "simplify_network\nsimpl: ", color = "0.11 0.6 0.85", style="rounded"]; - 5[label = "add_electricity", color = "0.23 0.6 0.85", style="rounded"]; - 6[label = "build_renewable_profiles\ntechnology: onwind", color = "0.57 0.6 0.85", style="rounded"]; - 7[label = "base_network", color = "0.09 0.6 0.85", style="rounded"]; - 8[label = "build_shapes", color = "0.41 0.6 0.85", style="rounded"]; - 9[label = "retrieve_databundle", color = "0.28 0.6 0.85", style="rounded"]; - 10[label = "retrieve_natura_raster", color = "0.62 0.6 0.85", style="rounded"]; - 11[label = "build_bus_regions", color = "0.53 0.6 0.85", style="rounded"]; - 12[label = "retrieve_cutout\ncutout: europe-2013-era5", color = "0.05 0.6 0.85", style="rounded,dashed"]; - 13[label = "build_renewable_profiles\ntechnology: offwind-ac", color = "0.57 0.6 0.85", style="rounded"]; - 14[label = "build_ship_raster", color = "0.64 0.6 0.85", style="rounded"]; - 15[label = "retrieve_ship_raster", color = "0.07 0.6 0.85", style="rounded,dashed"]; - 16[label = "retrieve_cutout\ncutout: europe-2013-sarah", color = "0.05 0.6 0.85", style="rounded,dashed"]; - 17[label = "build_renewable_profiles\ntechnology: offwind-dc", color = "0.57 0.6 0.85", style="rounded"]; - 18[label = "build_renewable_profiles\ntechnology: solar", color = "0.57 0.6 0.85", style="rounded"]; - 19[label = "build_hydro_profile", color = "0.44 0.6 0.85", style="rounded"]; - 20[label = "retrieve_cost_data", color = "0.30 0.6 0.85", style="rounded"]; - 21[label = "build_powerplants", color = "0.16 0.6 0.85", style="rounded"]; - 22[label = "build_electricity_demand", color = "0.00 0.6 0.85", style="rounded"]; - 23[label = "retrieve_electricity_demand", color = "0.34 0.6 0.85", style="rounded,dashed"]; - 1 -> 0 - 2 -> 1 - 20 -> 1 - 3 -> 2 - 20 -> 2 - 4 -> 3 - 20 -> 3 - 5 -> 4 - 20 -> 4 - 11 -> 4 - 6 -> 5 - 13 -> 5 - 17 -> 5 - 18 -> 5 - 19 -> 5 - 7 -> 5 - 20 -> 5 - 11 -> 5 - 21 -> 5 - 9 -> 5 - 22 -> 5 - 8 -> 5 - 7 -> 6 - 9 -> 6 - 10 -> 6 - 8 -> 6 - 11 -> 6 - 12 -> 6 - 8 -> 7 - 9 -> 8 - 8 -> 11 - 7 -> 11 - 7 -> 13 - 9 -> 13 - 10 -> 13 - 14 -> 13 - 8 -> 13 - 11 -> 13 - 12 -> 13 - 15 -> 14 - 12 -> 14 - 16 -> 14 - 7 -> 17 - 9 -> 17 - 10 -> 17 - 14 -> 17 - 8 -> 17 - 11 -> 17 - 12 -> 17 - 7 -> 18 - 9 -> 18 - 10 -> 18 - 8 -> 18 - 11 -> 18 - 16 -> 18 - 8 -> 19 - 12 -> 19 - 7 -> 21 - 23 -> 22 + 0[label = "solve_network", color = "0.33 0.6 0.85", style="rounded"]; + 1[label = "prepare_network\nll: copt\nopts: Co2L-24H", color = "0.03 0.6 0.85", style="rounded"]; + 2[label = "add_extra_components", color = "0.45 0.6 0.85", style="rounded"]; + 3[label = "cluster_network\nclusters: 6", color = "0.46 0.6 0.85", style="rounded"]; + 4[label = "simplify_network\nsimpl: ", color = "0.52 0.6 0.85", style="rounded"]; + 5[label = "add_electricity", color = "0.55 0.6 0.85", style="rounded"]; + 6[label = "build_renewable_profiles\ntechnology: solar", color = "0.15 0.6 0.85", style="rounded"]; + 7[label = "base_network", color = "0.37 0.6 0.85", style="rounded,dashed"]; + 8[label = "build_shapes", color = "0.07 0.6 0.85", style="rounded,dashed"]; + 9[label = "retrieve_databundle", color = "0.60 0.6 0.85", style="rounded"]; + 10[label = "retrieve_natura_raster", color = "0.42 0.6 0.85", style="rounded"]; + 11[label = "build_bus_regions", color = "0.09 0.6 0.85", style="rounded,dashed"]; + 12[label = "build_renewable_profiles\ntechnology: onwind", color = "0.15 0.6 0.85", style="rounded"]; + 13[label = "build_renewable_profiles\ntechnology: offwind-ac", color = "0.15 0.6 0.85", style="rounded"]; + 14[label = "build_ship_raster", color = "0.02 0.6 0.85", style="rounded"]; + 15[label = "retrieve_ship_raster", color = "0.40 0.6 0.85", style="rounded"]; + 16[label = "build_renewable_profiles\ntechnology: offwind-dc", color = "0.15 0.6 0.85", style="rounded"]; + 17[label = "build_line_rating", color = "0.32 0.6 0.85", style="rounded"]; + 18[label = "retrieve_cost_data\nyear: 2030", color = "0.50 0.6 0.85", style="rounded"]; + 19[label = "build_powerplants", color = "0.64 0.6 0.85", style="rounded,dashed"]; + 20[label = "build_electricity_demand", color = "0.13 0.6 0.85", style="rounded,dashed"]; + 21[label = "retrieve_electricity_demand", color = "0.31 0.6 0.85", style="rounded"]; + 22[label = "copy_config", color = "0.23 0.6 0.85", style="rounded"]; + 1 -> 0 + 22 -> 0 + 2 -> 1 + 18 -> 1 + 3 -> 2 + 18 -> 2 + 4 -> 3 + 18 -> 3 + 5 -> 4 + 18 -> 4 + 11 -> 4 + 6 -> 5 + 12 -> 5 + 13 -> 5 + 16 -> 5 + 7 -> 5 + 17 -> 5 + 18 -> 5 + 11 -> 5 + 19 -> 5 + 9 -> 5 + 20 -> 5 + 8 -> 5 + 7 -> 6 + 9 -> 6 + 10 -> 6 + 8 -> 6 + 11 -> 6 + 8 -> 7 + 9 -> 8 + 8 -> 11 + 7 -> 11 + 7 -> 12 + 9 -> 12 + 10 -> 12 + 8 -> 12 + 11 -> 12 + 7 -> 13 + 9 -> 13 + 10 -> 13 + 14 -> 13 + 8 -> 13 + 11 -> 13 + 15 -> 14 + 7 -> 16 + 9 -> 16 + 10 -> 16 + 14 -> 16 + 8 -> 16 + 11 -> 16 + 7 -> 17 + 7 -> 19 + 21 -> 20 } | diff --git a/envs/environment.yaml b/envs/environment.yaml index c3af36bb..09c209a6 100644 --- a/envs/environment.yaml +++ b/envs/environment.yaml @@ -55,5 +55,5 @@ dependencies: - pip: - - tsam>=1.1.0 - - pypsa>=0.25.1 + - git+https://github.com/fneum/tsam.git@performance + - pypsa>=0.25.2 diff --git a/rules/collect.smk b/rules/collect.smk index 0c713e5a..c9bb10ea 100644 --- a/rules/collect.smk +++ b/rules/collect.smk @@ -60,6 +60,15 @@ rule solve_sector_networks: ), +rule solve_sector_networks_perfect: + input: + expand( + RESULTS + + "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years.nc", + **config["scenario"] + ), + + rule plot_networks: input: expand( diff --git a/rules/postprocess.smk b/rules/postprocess.smk index cf0038a3..c34fe27e 100644 --- a/rules/postprocess.smk +++ b/rules/postprocess.smk @@ -8,31 +8,62 @@ localrules: copy_conda_env, -rule plot_network: - params: - foresight=config["foresight"], - plotting=config["plotting"], - input: - network=RESULTS - + "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", - regions=RESOURCES + "regions_onshore_elec_s{simpl}_{clusters}.geojson", - output: - map=RESULTS - + "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}-costs-all_{planning_horizons}.pdf", - today=RESULTS - + "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}-today.pdf", - threads: 2 - resources: - mem_mb=10000, - benchmark: - ( +if config["foresight"] != "perfect": + + rule plot_network: + params: + foresight=config["foresight"], + plotting=config["plotting"], + input: + network=RESULTS + + "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", + regions=RESOURCES + "regions_onshore_elec_s{simpl}_{clusters}.geojson", + output: + map=RESULTS + + "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}-costs-all_{planning_horizons}.pdf", + today=RESULTS + + "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}-today.pdf", + threads: 2 + resources: + mem_mb=10000, + benchmark: + ( + BENCHMARKS + + "plot_network/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}" + ) + conda: + "../envs/environment.yaml" + script: + "../scripts/plot_network.py" + + +if config["foresight"] == "perfect": + + rule plot_network: + params: + foresight=config["foresight"], + plotting=config["plotting"], + input: + network=RESULTS + + "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years.nc", + regions=RESOURCES + "regions_onshore_elec_s{simpl}_{clusters}.geojson", + output: + **{ + f"map_{year}": RESULTS + + "maps/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}-costs-all_" + + f"{year}.pdf" + for year in config["scenario"]["planning_horizons"] + }, + threads: 2 + resources: + mem_mb=10000, + benchmark: BENCHMARKS - + "plot_network/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}" - ) - conda: - "../envs/environment.yaml" - script: - "../scripts/plot_network.py" + +"postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years_benchmark" + conda: + "../envs/environment.yaml" + script: + "../scripts/plot_network.py" rule copy_config: diff --git a/rules/solve_perfect.smk b/rules/solve_perfect.smk new file mode 100644 index 00000000..ef4e367d --- /dev/null +++ b/rules/solve_perfect.smk @@ -0,0 +1,194 @@ +# SPDX-FileCopyrightText: : 2023 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: MIT +rule add_existing_baseyear: + params: + baseyear=config["scenario"]["planning_horizons"][0], + sector=config["sector"], + existing_capacities=config["existing_capacities"], + costs=config["costs"], + input: + network=RESULTS + + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", + powerplants=RESOURCES + "powerplants.csv", + busmap_s=RESOURCES + "busmap_elec_s{simpl}.csv", + busmap=RESOURCES + "busmap_elec_s{simpl}_{clusters}.csv", + clustered_pop_layout=RESOURCES + "pop_layout_elec_s{simpl}_{clusters}.csv", + costs="data/costs_{}.csv".format(config["scenario"]["planning_horizons"][0]), + cop_soil_total=RESOURCES + "cop_soil_total_elec_s{simpl}_{clusters}.nc", + cop_air_total=RESOURCES + "cop_air_total_elec_s{simpl}_{clusters}.nc", + existing_heating="data/existing_infrastructure/existing_heating_raw.csv", + existing_solar="data/existing_infrastructure/solar_capacity_IRENA.csv", + existing_onwind="data/existing_infrastructure/onwind_capacity_IRENA.csv", + existing_offwind="data/existing_infrastructure/offwind_capacity_IRENA.csv", + output: + RESULTS + + "prenetworks-brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", + wildcard_constraints: + planning_horizons=config["scenario"]["planning_horizons"][0], #only applies to baseyear + threads: 1 + resources: + mem_mb=2000, + log: + LOGS + + "add_existing_baseyear_elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.log", + benchmark: + ( + BENCHMARKS + + "add_existing_baseyear/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}" + ) + conda: + "../envs/environment.yaml" + script: + "../scripts/add_existing_baseyear.py" + + +rule add_brownfield: + params: + H2_retrofit=config["sector"]["H2_retrofit"], + H2_retrofit_capacity_per_CH4=config["sector"]["H2_retrofit_capacity_per_CH4"], + threshold_capacity=config["existing_capacities"]["threshold_capacity"], + input: + network=RESULTS + + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", + network_p=solved_previous_horizon, #solved network at previous time step + costs="data/costs_{planning_horizons}.csv", + cop_soil_total=RESOURCES + "cop_soil_total_elec_s{simpl}_{clusters}.nc", + cop_air_total=RESOURCES + "cop_air_total_elec_s{simpl}_{clusters}.nc", + output: + RESULTS + + "prenetworks-brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", + threads: 4 + resources: + mem_mb=10000, + log: + LOGS + + "add_brownfield_elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.log", + benchmark: + ( + BENCHMARKS + + "add_brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}" + ) + conda: + "../envs/environment.yaml" + script: + "../scripts/add_brownfield.py" + + +rule prepare_perfect_foresight: + input: + **{ + f"network_{year}": RESULTS + + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_" + + f"{year}.nc" + for year in config["scenario"]["planning_horizons"][1:] + }, + brownfield_network=lambda w: ( + RESULTS + + "prenetworks-brownfield/" + + "elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_" + + "{}.nc".format(str(config["scenario"]["planning_horizons"][0])) + ), + output: + RESULTS + + "prenetworks-brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years.nc", + threads: 2 + resources: + mem_mb=10000, + log: + LOGS + + "prepare_perfect_foresight{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}.log", + benchmark: + ( + BENCHMARKS + + "prepare_perfect_foresight{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}" + ) + conda: + "../envs/environment.yaml" + script: + "../scripts/prepare_perfect_foresight.py" + + +rule solve_sector_network_perfect: + params: + solving=config["solving"], + foresight=config["foresight"], + sector=config["sector"], + planning_horizons=config["scenario"]["planning_horizons"], + co2_sequestration_potential=config["sector"].get( + "co2_sequestration_potential", 200 + ), + input: + network=RESULTS + + "prenetworks-brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years.nc", + costs="data/costs_2030.csv", + config=RESULTS + "config.yaml", + output: + RESULTS + + "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years.nc", + threads: 4 + resources: + mem_mb=config["solving"]["mem"], + shadow: + "shallow" + log: + solver=RESULTS + + "logs/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years_solver.log", + python=RESULTS + + "logs/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years_python.log", + memory=RESULTS + + "logs/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years_memory.log", + benchmark: + ( + BENCHMARKS + + "solve_sector_network/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years}" + ) + conda: + "../envs/environment.yaml" + script: + "../scripts/solve_network.py" + + +rule make_summary_perfect: + input: + **{ + f"networks_{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}": RESULTS + + f"postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_brownfield_all_years.nc" + for simpl in config["scenario"]["simpl"] + for clusters in config["scenario"]["clusters"] + for opts in config["scenario"]["opts"] + for sector_opts in config["scenario"]["sector_opts"] + for ll in config["scenario"]["ll"] + }, + costs="data/costs_2020.csv", + output: + nodal_costs=RESULTS + "csvs/nodal_costs.csv", + nodal_capacities=RESULTS + "csvs/nodal_capacities.csv", + nodal_cfs=RESULTS + "csvs/nodal_cfs.csv", + cfs=RESULTS + "csvs/cfs.csv", + costs=RESULTS + "csvs/costs.csv", + capacities=RESULTS + "csvs/capacities.csv", + curtailment=RESULTS + "csvs/curtailment.csv", + energy=RESULTS + "csvs/energy.csv", + supply=RESULTS + "csvs/supply.csv", + supply_energy=RESULTS + "csvs/supply_energy.csv", + prices=RESULTS + "csvs/prices.csv", + weighted_prices=RESULTS + "csvs/weighted_prices.csv", + market_values=RESULTS + "csvs/market_values.csv", + price_statistics=RESULTS + "csvs/price_statistics.csv", + metrics=RESULTS + "csvs/metrics.csv", + co2_emissions=RESULTS + "csvs/co2_emissions.csv", + threads: 2 + resources: + mem_mb=10000, + log: + LOGS + "make_summary_perfect.log", + benchmark: + (BENCHMARKS + "make_summary_perfect") + conda: + "../envs/environment.yaml" + script: + "../scripts/make_summary_perfect.py" + + +ruleorder: add_existing_baseyear > add_brownfield diff --git a/scripts/_benchmark.py b/scripts/_benchmark.py new file mode 100644 index 00000000..4e3413e9 --- /dev/null +++ b/scripts/_benchmark.py @@ -0,0 +1,256 @@ +# -*- coding: utf-8 -*- +# SPDX-FileCopyrightText: : 2020-2023 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: MIT +""" + +""" + +from __future__ import absolute_import, print_function + +import logging +import os +import sys +import time + +logger = logging.getLogger(__name__) + +# TODO: provide alternative when multiprocessing is not available +try: + from multiprocessing import Pipe, Process +except ImportError: + from multiprocessing.dummy import Process, Pipe + +from memory_profiler import _get_memory, choose_backend + + +# The memory logging facilities have been adapted from memory_profiler +class MemTimer(Process): + """ + Write memory consumption over a time interval to file until signaled to + stop on the pipe. + """ + + def __init__( + self, monitor_pid, interval, pipe, filename, max_usage, backend, *args, **kw + ): + self.monitor_pid = monitor_pid + self.interval = interval + self.pipe = pipe + self.filename = filename + self.max_usage = max_usage + self.backend = backend + + self.timestamps = kw.pop("timestamps", True) + self.include_children = kw.pop("include_children", True) + + super(MemTimer, self).__init__(*args, **kw) + + def run(self): + # get baseline memory usage + cur_mem = _get_memory( + self.monitor_pid, + self.backend, + timestamps=self.timestamps, + include_children=self.include_children, + ) + + n_measurements = 1 + mem_usage = cur_mem if self.max_usage else [cur_mem] + + if self.filename is not None: + stream = open(self.filename, "w") + stream.write("MEM {0:.6f} {1:.4f}\n".format(*cur_mem)) + stream.flush() + else: + stream = None + + self.pipe.send(0) # we're ready + stop = False + while True: + cur_mem = _get_memory( + self.monitor_pid, + self.backend, + timestamps=self.timestamps, + include_children=self.include_children, + ) + + if stream is not None: + stream.write("MEM {0:.6f} {1:.4f}\n".format(*cur_mem)) + stream.flush() + + n_measurements += 1 + if not self.max_usage: + mem_usage.append(cur_mem) + else: + mem_usage = max(cur_mem, mem_usage) + + if stop: + break + stop = self.pipe.poll(self.interval) + # do one more iteration + + if stream is not None: + stream.close() + + self.pipe.send(mem_usage) + self.pipe.send(n_measurements) + + +class memory_logger(object): + """ + Context manager for taking and reporting memory measurements at fixed + intervals from a separate process, for the duration of a context. + + Parameters + ---------- + filename : None|str + Name of the text file to log memory measurements, if None no log is + created (defaults to None) + interval : float + Interval between measurements (defaults to 1.) + max_usage : bool + If True, only store and report the maximum value (defaults to True) + timestamps : bool + Whether to record tuples of memory usage and timestamps; if logging to + a file timestamps are always kept (defaults to True) + include_children : bool + Whether the memory of subprocesses is to be included (default: True) + + Arguments + --------- + n_measurements : int + Number of measurements that have been taken + mem_usage : (float, float)|[(float, float)] + All memory measurements and timestamps (if timestamps was True) or only + the maximum memory usage and its timestamp + + Note + ---- + The arguments are only set after all the measurements, i.e. outside of the + with statement. + + Example + ------- + with memory_logger(filename="memory.log", max_usage=True) as mem: + # Do a lot of long running memory intensive stuff + hard_memory_bound_stuff() + + max_mem, timestamp = mem.mem_usage + """ + + def __init__( + self, + filename=None, + interval=1.0, + max_usage=True, + timestamps=True, + include_children=True, + ): + if filename is not None: + timestamps = True + + self.filename = filename + self.interval = interval + self.max_usage = max_usage + self.timestamps = timestamps + self.include_children = include_children + + def __enter__(self): + backend = choose_backend() + + self.child_conn, self.parent_conn = Pipe() # this will store MemTimer's results + self.p = MemTimer( + os.getpid(), + self.interval, + self.child_conn, + self.filename, + backend=backend, + timestamps=self.timestamps, + max_usage=self.max_usage, + include_children=self.include_children, + ) + self.p.start() + self.parent_conn.recv() # wait until memory logging in subprocess is ready + + return self + + def __exit__(self, exc_type, exc_val, exc_tb): + if exc_type is None: + self.parent_conn.send(0) # finish timing + + self.mem_usage = self.parent_conn.recv() + self.n_measurements = self.parent_conn.recv() + else: + self.p.terminate() + + return False + + +class timer(object): + level = 0 + opened = False + + def __init__(self, name="", verbose=True): + self.name = name + self.verbose = verbose + + def __enter__(self): + if self.verbose: + if self.opened: + sys.stdout.write("\n") + + if len(self.name) > 0: + sys.stdout.write((".. " * self.level) + self.name + ": ") + sys.stdout.flush() + + self.__class__.opened = True + + self.__class__.level += 1 + + self.start = time.time() + return self + + def print_usec(self, usec): + if usec < 1000: + print("%.1f usec" % usec) + else: + msec = usec / 1000 + if msec < 1000: + print("%.1f msec" % msec) + else: + sec = msec / 1000 + print("%.1f sec" % sec) + + def __exit__(self, exc_type, exc_val, exc_tb): + if not self.opened and self.verbose: + sys.stdout.write(".. " * self.level) + + if exc_type is None: + stop = time.time() + self.usec = usec = (stop - self.start) * 1e6 + if self.verbose: + self.print_usec(usec) + elif self.verbose: + print("failed") + sys.stdout.flush() + + self.__class__.level -= 1 + if self.verbose: + self.__class__.opened = False + return False + + +class optional(object): + def __init__(self, variable, contextman): + self.variable = variable + self.contextman = contextman + + def __enter__(self): + if self.variable: + return self.contextman.__enter__() + + def __exit__(self, exc_type, exc_val, exc_tb): + if self.variable: + return self.contextman.__exit__(exc_type, exc_val, exc_tb) + return False diff --git a/scripts/add_electricity.py b/scripts/add_electricity.py index 27793396..782b873c 100755 --- a/scripts/add_electricity.py +++ b/scripts/add_electricity.py @@ -406,6 +406,7 @@ def attach_wind_and_solar( capital_cost=capital_cost, efficiency=costs.at[supcar, "efficiency"], p_max_pu=ds["profile"].transpose("time", "bus").to_pandas(), + lifetime=costs.at[supcar, "lifetime"], ) diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py index 08810470..01c69997 100644 --- a/scripts/add_existing_baseyear.py +++ b/scripts/add_existing_baseyear.py @@ -305,6 +305,18 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas if "EU" not in vars(spatial)[carrier[generator]].locations: bus0 = bus0.intersection(capacity.index + " gas") + # check for missing bus + missing_bus = pd.Index(bus0).difference(n.buses.index) + if not missing_bus.empty: + logger.info(f"add buses {bus0}") + n.madd( + "Bus", + bus0, + carrier=generator, + location=vars(spatial)[carrier[generator]].locations, + unit="MWh_el", + ) + already_build = n.links.index.intersection(asset_i) new_build = asset_i.difference(n.links.index) lifetime_assets = lifetime.loc[grouping_year, generator].dropna() @@ -605,6 +617,10 @@ def add_heating_capacities_installed_before_baseyear( ], ) + # drop assets which are at the end of their lifetime + links_i = n.links[(n.links.build_year + n.links.lifetime <= baseyear)].index + n.mremove("Link", links_i) + # %% if __name__ == "__main__": @@ -613,13 +629,13 @@ if __name__ == "__main__": snakemake = mock_snakemake( "add_existing_baseyear", - configfiles="config/test/config.myopic.yaml", + # configfiles="config/test/config.myopic.yaml", simpl="", - clusters="5", - ll="v1.5", + clusters="37", + ll="v1.0", opts="", - sector_opts="24H-T-H-B-I-A-solar+p3-dist1", - planning_horizons=2030, + sector_opts="1p7-4380H-T-H-B-I-A-solar+p3-dist1", + planning_horizons=2020, ) logging.basicConfig(level=snakemake.config["logging"]["level"]) diff --git a/scripts/build_industrial_distribution_key.py b/scripts/build_industrial_distribution_key.py index bc4a26bc..888a1a96 100644 --- a/scripts/build_industrial_distribution_key.py +++ b/scripts/build_industrial_distribution_key.py @@ -96,6 +96,23 @@ def prepare_hotmaps_database(regions): gdf.rename(columns={"index_right": "bus"}, inplace=True) gdf["country"] = gdf.bus.str[:2] + # the .sjoin can lead to duplicates if a geom is in two regions + if gdf.index.duplicated().any(): + import pycountry + + # get all duplicated entries + duplicated_i = gdf.index[gdf.index.duplicated()] + # convert from raw data country name to iso-2-code + s = df.loc[duplicated_i, "Country"].apply( + lambda x: pycountry.countries.lookup(x).alpha_2 + ) + # Get a boolean mask where gdf's country column matches s's values for the same index + mask = gdf["country"] == gdf.index.map(s) + # Filter gdf using the mask + gdf_filtered = gdf[mask] + # concat not duplicated and filtered gdf + gdf = pd.concat([gdf.drop(duplicated_i), gdf_filtered]).sort_index() + # the .sjoin can lead to duplicates if a geom is in two overlapping regions if gdf.index.duplicated().any(): # get all duplicated entries @@ -147,6 +164,7 @@ def build_nodal_distribution_key(hotmaps, regions, countries): return keys +# %% if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake @@ -154,7 +172,7 @@ if __name__ == "__main__": snakemake = mock_snakemake( "build_industrial_distribution_key", simpl="", - clusters=48, + clusters=128, ) logging.basicConfig(level=snakemake.config["logging"]["level"]) diff --git a/scripts/build_line_rating.py b/scripts/build_line_rating.py index 7f842d43..d3a970bd 100755 --- a/scripts/build_line_rating.py +++ b/scripts/build_line_rating.py @@ -41,7 +41,7 @@ The following heat gains and losses are considered: - heat gain through resistive losses - heat gain through solar radiation -- heat loss through radiation of the trasnmission line +- heat loss through radiation of the transmission line - heat loss through forced convection with wind - heat loss through natural convection diff --git a/scripts/make_summary_perfect.py b/scripts/make_summary_perfect.py new file mode 100644 index 00000000..1a3391a9 --- /dev/null +++ b/scripts/make_summary_perfect.py @@ -0,0 +1,745 @@ +# -*- coding: utf-8 -*- +# SPDX-FileCopyrightText: : 2020-2023 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: MIT +""" +Create summary CSV files for all scenario runs with perfect foresight including +costs, capacities, capacity factors, curtailment, energy balances, prices and +other metrics. +""" + + +import numpy as np +import pandas as pd +import pypsa +from make_summary import ( + assign_carriers, + assign_locations, + calculate_cfs, + calculate_nodal_cfs, + calculate_nodal_costs, +) +from prepare_sector_network import prepare_costs +from pypsa.descriptors import get_active_assets, nominal_attrs +from six import iteritems + +idx = pd.IndexSlice + +opt_name = {"Store": "e", "Line": "s", "Transformer": "s"} + + +def calculate_costs(n, label, costs): + investments = n.investment_periods + cols = pd.MultiIndex.from_product( + [ + costs.columns.levels[0], + costs.columns.levels[1], + costs.columns.levels[2], + investments, + ], + names=costs.columns.names[:3] + ["year"], + ) + costs = costs.reindex(cols, axis=1) + + for c in n.iterate_components( + n.branch_components | n.controllable_one_port_components ^ {"Load"} + ): + capital_costs = c.df.capital_cost * c.df[opt_name.get(c.name, "p") + "_nom_opt"] + active = pd.concat( + [ + get_active_assets(n, c.name, inv_p).rename(inv_p) + for inv_p in investments + ], + axis=1, + ).astype(int) + capital_costs = active.mul(capital_costs, axis=0) + discount = ( + n.investment_period_weightings["objective"] + / n.investment_period_weightings["years"] + ) + capital_costs_grouped = capital_costs.groupby(c.df.carrier).sum().mul(discount) + + capital_costs_grouped = pd.concat([capital_costs_grouped], keys=["capital"]) + capital_costs_grouped = pd.concat([capital_costs_grouped], keys=[c.list_name]) + + costs = costs.reindex(capital_costs_grouped.index.union(costs.index)) + + costs.loc[capital_costs_grouped.index, label] = capital_costs_grouped.values + + if c.name == "Link": + p = ( + c.pnl.p0.multiply(n.snapshot_weightings.generators, axis=0) + .groupby(level=0) + .sum() + ) + elif c.name == "Line": + continue + elif c.name == "StorageUnit": + p_all = c.pnl.p.multiply(n.snapshot_weightings.stores, axis=0) + p_all[p_all < 0.0] = 0.0 + p = p_all.groupby(level=0).sum() + else: + p = ( + round(c.pnl.p, ndigits=2) + .multiply(n.snapshot_weightings.generators, axis=0) + .groupby(level=0) + .sum() + ) + + # correct sequestration cost + if c.name == "Store": + items = c.df.index[ + (c.df.carrier == "co2 stored") & (c.df.marginal_cost <= -100.0) + ] + c.df.loc[items, "marginal_cost"] = -20.0 + + marginal_costs = p.mul(c.df.marginal_cost).T + # marginal_costs = active.mul(marginal_costs, axis=0) + marginal_costs_grouped = ( + marginal_costs.groupby(c.df.carrier).sum().mul(discount) + ) + + marginal_costs_grouped = pd.concat([marginal_costs_grouped], keys=["marginal"]) + marginal_costs_grouped = pd.concat([marginal_costs_grouped], keys=[c.list_name]) + + costs = costs.reindex(marginal_costs_grouped.index.union(costs.index)) + + costs.loc[marginal_costs_grouped.index, label] = marginal_costs_grouped.values + + # add back in all hydro + # costs.loc[("storage_units","capital","hydro"),label] = (0.01)*2e6*n.storage_units.loc[n.storage_units.group=="hydro","p_nom"].sum() + # costs.loc[("storage_units","capital","PHS"),label] = (0.01)*2e6*n.storage_units.loc[n.storage_units.group=="PHS","p_nom"].sum() + # costs.loc[("generators","capital","ror"),label] = (0.02)*3e6*n.generators.loc[n.generators.group=="ror","p_nom"].sum() + + return costs + + +def calculate_cumulative_cost(): + planning_horizons = snakemake.config["scenario"]["planning_horizons"] + + cumulative_cost = pd.DataFrame( + index=df["costs"].sum().index, + columns=pd.Series(data=np.arange(0, 0.1, 0.01), name="social discount rate"), + ) + + # discount cost and express them in money value of planning_horizons[0] + for r in cumulative_cost.columns: + cumulative_cost[r] = [ + df["costs"].sum()[index] / ((1 + r) ** (index[-1] - planning_horizons[0])) + for index in cumulative_cost.index + ] + + # integrate cost throughout the transition path + for r in cumulative_cost.columns: + for cluster in cumulative_cost.index.get_level_values(level=0).unique(): + for lv in cumulative_cost.index.get_level_values(level=1).unique(): + for sector_opts in cumulative_cost.index.get_level_values( + level=2 + ).unique(): + cumulative_cost.loc[ + (cluster, lv, sector_opts, "cumulative cost"), r + ] = np.trapz( + cumulative_cost.loc[ + idx[cluster, lv, sector_opts, planning_horizons], r + ].values, + x=planning_horizons, + ) + + return cumulative_cost + + +def calculate_nodal_capacities(n, label, nodal_capacities): + # Beware this also has extraneous locations for country (e.g. biomass) or continent-wide (e.g. fossil gas/oil) stuff + for c in n.iterate_components( + n.branch_components | n.controllable_one_port_components ^ {"Load"} + ): + nodal_capacities_c = c.df.groupby(["location", "carrier"])[ + opt_name.get(c.name, "p") + "_nom_opt" + ].sum() + index = pd.MultiIndex.from_tuples( + [(c.list_name,) + t for t in nodal_capacities_c.index.to_list()] + ) + nodal_capacities = nodal_capacities.reindex(index.union(nodal_capacities.index)) + nodal_capacities.loc[index, label] = nodal_capacities_c.values + + return nodal_capacities + + +def calculate_capacities(n, label, capacities): + investments = n.investment_periods + cols = pd.MultiIndex.from_product( + [ + capacities.columns.levels[0], + capacities.columns.levels[1], + capacities.columns.levels[2], + investments, + ], + names=capacities.columns.names[:3] + ["year"], + ) + capacities = capacities.reindex(cols, axis=1) + + for c in n.iterate_components( + n.branch_components | n.controllable_one_port_components ^ {"Load"} + ): + active = pd.concat( + [ + get_active_assets(n, c.name, inv_p).rename(inv_p) + for inv_p in investments + ], + axis=1, + ).astype(int) + caps = c.df[opt_name.get(c.name, "p") + "_nom_opt"] + caps = active.mul(caps, axis=0) + capacities_grouped = ( + caps.groupby(c.df.carrier).sum().drop("load", errors="ignore") + ) + capacities_grouped = pd.concat([capacities_grouped], keys=[c.list_name]) + + capacities = capacities.reindex( + capacities_grouped.index.union(capacities.index) + ) + + capacities.loc[capacities_grouped.index, label] = capacities_grouped.values + + return capacities + + +def calculate_curtailment(n, label, curtailment): + avail = ( + n.generators_t.p_max_pu.multiply(n.generators.p_nom_opt) + .sum() + .groupby(n.generators.carrier) + .sum() + ) + used = n.generators_t.p.sum().groupby(n.generators.carrier).sum() + + curtailment[label] = (((avail - used) / avail) * 100).round(3) + + return curtailment + + +def calculate_energy(n, label, energy): + investments = n.investment_periods + cols = pd.MultiIndex.from_product( + [ + energy.columns.levels[0], + energy.columns.levels[1], + energy.columns.levels[2], + investments, + ], + names=energy.columns.names[:3] + ["year"], + ) + energy = energy.reindex(cols, axis=1) + + for c in n.iterate_components(n.one_port_components | n.branch_components): + if c.name in n.one_port_components: + c_energies = ( + c.pnl.p.multiply(n.snapshot_weightings.generators, axis=0) + .groupby(level=0) + .sum() + .multiply(c.df.sign) + .groupby(c.df.carrier, axis=1) + .sum() + ) + else: + c_energies = pd.DataFrame( + 0.0, columns=c.df.carrier.unique(), index=n.investment_periods + ) + for port in [col[3:] for col in c.df.columns if col[:3] == "bus"]: + totals = ( + c.pnl["p" + port] + .multiply(n.snapshot_weightings.generators, axis=0) + .groupby(level=0) + .sum() + ) + # remove values where bus is missing (bug in nomopyomo) + no_bus = c.df.index[c.df["bus" + port] == ""] + totals[no_bus] = float( + n.component_attrs[c.name].loc["p" + port, "default"] + ) + c_energies -= totals.groupby(c.df.carrier, axis=1).sum() + + c_energies = pd.concat([c_energies.T], keys=[c.list_name]) + + energy = energy.reindex(c_energies.index.union(energy.index)) + + energy.loc[c_energies.index, label] = c_energies.values + + return energy + + +def calculate_supply(n, label, supply): + """ + Calculate the max dispatch of each component at the buses aggregated by + carrier. + """ + + bus_carriers = n.buses.carrier.unique() + + for i in bus_carriers: + bus_map = n.buses.carrier == i + bus_map.at[""] = False + + for c in n.iterate_components(n.one_port_components): + items = c.df.index[c.df.bus.map(bus_map).fillna(False)] + + if len(items) == 0: + continue + + s = ( + c.pnl.p[items] + .max() + .multiply(c.df.loc[items, "sign"]) + .groupby(c.df.loc[items, "carrier"]) + .sum() + ) + s = pd.concat([s], keys=[c.list_name]) + s = pd.concat([s], keys=[i]) + + supply = supply.reindex(s.index.union(supply.index)) + supply.loc[s.index, label] = s + + for c in n.iterate_components(n.branch_components): + for end in [col[3:] for col in c.df.columns if col[:3] == "bus"]: + items = c.df.index[c.df["bus" + end].map(bus_map).fillna(False)] + + if len(items) == 0: + continue + + # lots of sign compensation for direction and to do maximums + s = (-1) ** (1 - int(end)) * ( + (-1) ** int(end) * c.pnl["p" + end][items] + ).max().groupby(c.df.loc[items, "carrier"]).sum() + s.index = s.index + end + s = pd.concat([s], keys=[c.list_name]) + s = pd.concat([s], keys=[i]) + + supply = supply.reindex(s.index.union(supply.index)) + supply.loc[s.index, label] = s + + return supply + + +def calculate_supply_energy(n, label, supply_energy): + """ + Calculate the total energy supply/consuption of each component at the buses + aggregated by carrier. + """ + + investments = n.investment_periods + cols = pd.MultiIndex.from_product( + [ + supply_energy.columns.levels[0], + supply_energy.columns.levels[1], + supply_energy.columns.levels[2], + investments, + ], + names=supply_energy.columns.names[:3] + ["year"], + ) + supply_energy = supply_energy.reindex(cols, axis=1) + + bus_carriers = n.buses.carrier.unique() + + for i in bus_carriers: + bus_map = n.buses.carrier == i + bus_map.at[""] = False + + for c in n.iterate_components(n.one_port_components): + items = c.df.index[c.df.bus.map(bus_map).fillna(False)] + + if len(items) == 0: + continue + + if c.name == "Generator": + weightings = n.snapshot_weightings.generators + else: + weightings = n.snapshot_weightings.stores + + if i in ["oil", "co2", "H2"]: + if c.name == "Load": + c.df.loc[items, "carrier"] = [ + load.split("-202")[0] for load in items + ] + if i == "oil" and c.name == "Generator": + c.df.loc[items, "carrier"] = "imported oil" + s = ( + c.pnl.p[items] + .multiply(weightings, axis=0) + .groupby(level=0) + .sum() + .multiply(c.df.loc[items, "sign"]) + .groupby(c.df.loc[items, "carrier"], axis=1) + .sum() + .T + ) + s = pd.concat([s], keys=[c.list_name]) + s = pd.concat([s], keys=[i]) + + supply_energy = supply_energy.reindex( + s.index.union(supply_energy.index, sort=False) + ) + supply_energy.loc[s.index, label] = s.values + + for c in n.iterate_components(n.branch_components): + for end in [col[3:] for col in c.df.columns if col[:3] == "bus"]: + items = c.df.index[c.df["bus" + str(end)].map(bus_map).fillna(False)] + + if len(items) == 0: + continue + + s = ( + (-1) + * c.pnl["p" + end] + .reindex(items, axis=1) + .multiply(n.snapshot_weightings.objective, axis=0) + .groupby(level=0) + .sum() + .groupby(c.df.loc[items, "carrier"], axis=1) + .sum() + ).T + s.index = s.index + end + s = pd.concat([s], keys=[c.list_name]) + s = pd.concat([s], keys=[i]) + + supply_energy = supply_energy.reindex( + s.index.union(supply_energy.index, sort=False) + ) + + supply_energy.loc[s.index, label] = s.values + + return supply_energy + + +def calculate_metrics(n, label, metrics): + metrics = metrics.reindex( + pd.Index( + [ + "line_volume", + "line_volume_limit", + "line_volume_AC", + "line_volume_DC", + "line_volume_shadow", + "co2_shadow", + ] + ).union(metrics.index) + ) + + metrics.at["line_volume_DC", label] = (n.links.length * n.links.p_nom_opt)[ + n.links.carrier == "DC" + ].sum() + metrics.at["line_volume_AC", label] = (n.lines.length * n.lines.s_nom_opt).sum() + metrics.at["line_volume", label] = metrics.loc[ + ["line_volume_AC", "line_volume_DC"], label + ].sum() + + if hasattr(n, "line_volume_limit"): + metrics.at["line_volume_limit", label] = n.line_volume_limit + metrics.at["line_volume_shadow", label] = n.line_volume_limit_dual + + if "CO2Limit" in n.global_constraints.index: + metrics.at["co2_shadow", label] = n.global_constraints.at["CO2Limit", "mu"] + + return metrics + + +def calculate_prices(n, label, prices): + prices = prices.reindex(prices.index.union(n.buses.carrier.unique())) + + # WARNING: this is time-averaged, see weighted_prices for load-weighted average + prices[label] = n.buses_t.marginal_price.mean().groupby(n.buses.carrier).mean() + + return prices + + +def calculate_weighted_prices(n, label, weighted_prices): + # Warning: doesn't include storage units as loads + + weighted_prices = weighted_prices.reindex( + pd.Index( + [ + "electricity", + "heat", + "space heat", + "urban heat", + "space urban heat", + "gas", + "H2", + ] + ) + ) + + link_loads = { + "electricity": [ + "heat pump", + "resistive heater", + "battery charger", + "H2 Electrolysis", + ], + "heat": ["water tanks charger"], + "urban heat": ["water tanks charger"], + "space heat": [], + "space urban heat": [], + "gas": ["OCGT", "gas boiler", "CHP electric", "CHP heat"], + "H2": ["Sabatier", "H2 Fuel Cell"], + } + + for carrier in link_loads: + if carrier == "electricity": + suffix = "" + elif carrier[:5] == "space": + suffix = carrier[5:] + else: + suffix = " " + carrier + + buses = n.buses.index[n.buses.index.str[2:] == suffix] + + if buses.empty: + continue + + if carrier in ["H2", "gas"]: + load = pd.DataFrame(index=n.snapshots, columns=buses, data=0.0) + else: + load = n.loads_t.p_set.reindex(buses, axis=1) + + for tech in link_loads[carrier]: + names = n.links.index[n.links.index.to_series().str[-len(tech) :] == tech] + + if names.empty: + continue + + load += ( + n.links_t.p0[names].groupby(n.links.loc[names, "bus0"], axis=1).sum() + ) + + # Add H2 Store when charging + # if carrier == "H2": + # stores = n.stores_t.p[buses+ " Store"].groupby(n.stores.loc[buses+ " Store","bus"],axis=1).sum(axis=1) + # stores[stores > 0.] = 0. + # load += -stores + + weighted_prices.loc[carrier, label] = ( + load * n.buses_t.marginal_price[buses] + ).sum().sum() / load.sum().sum() + + if carrier[:5] == "space": + print(load * n.buses_t.marginal_price[buses]) + + return weighted_prices + + +def calculate_market_values(n, label, market_values): + # Warning: doesn't include storage units + + carrier = "AC" + + buses = n.buses.index[n.buses.carrier == carrier] + + ## First do market value of generators ## + + generators = n.generators.index[n.buses.loc[n.generators.bus, "carrier"] == carrier] + + techs = n.generators.loc[generators, "carrier"].value_counts().index + + market_values = market_values.reindex(market_values.index.union(techs)) + + for tech in techs: + gens = generators[n.generators.loc[generators, "carrier"] == tech] + + dispatch = ( + n.generators_t.p[gens] + .groupby(n.generators.loc[gens, "bus"], axis=1) + .sum() + .reindex(columns=buses, fill_value=0.0) + ) + + revenue = dispatch * n.buses_t.marginal_price[buses] + + market_values.at[tech, label] = revenue.sum().sum() / dispatch.sum().sum() + + ## Now do market value of links ## + + for i in ["0", "1"]: + all_links = n.links.index[n.buses.loc[n.links["bus" + i], "carrier"] == carrier] + + techs = n.links.loc[all_links, "carrier"].value_counts().index + + market_values = market_values.reindex(market_values.index.union(techs)) + + for tech in techs: + links = all_links[n.links.loc[all_links, "carrier"] == tech] + + dispatch = ( + n.links_t["p" + i][links] + .groupby(n.links.loc[links, "bus" + i], axis=1) + .sum() + .reindex(columns=buses, fill_value=0.0) + ) + + revenue = dispatch * n.buses_t.marginal_price[buses] + + market_values.at[tech, label] = revenue.sum().sum() / dispatch.sum().sum() + + return market_values + + +def calculate_price_statistics(n, label, price_statistics): + price_statistics = price_statistics.reindex( + price_statistics.index.union( + pd.Index(["zero_hours", "mean", "standard_deviation"]) + ) + ) + + buses = n.buses.index[n.buses.carrier == "AC"] + + threshold = 0.1 # higher than phoney marginal_cost of wind/solar + + df = pd.DataFrame(data=0.0, columns=buses, index=n.snapshots) + + df[n.buses_t.marginal_price[buses] < threshold] = 1.0 + + price_statistics.at["zero_hours", label] = df.sum().sum() / ( + df.shape[0] * df.shape[1] + ) + + price_statistics.at["mean", label] = n.buses_t.marginal_price[buses].mean().mean() + + price_statistics.at["standard_deviation", label] = ( + n.buses_t.marginal_price[buses].droplevel(0).unstack().std() + ) + + return price_statistics + + +def calculate_co2_emissions(n, label, df): + carattr = "co2_emissions" + emissions = n.carriers.query(f"{carattr} != 0")[carattr] + + if emissions.empty: + return + + weightings = n.snapshot_weightings.generators.mul( + n.investment_period_weightings["years"] + .reindex(n.snapshots) + .fillna(method="bfill") + .fillna(1.0), + axis=0, + ) + + # generators + gens = n.generators.query("carrier in @emissions.index") + if not gens.empty: + em_pu = gens.carrier.map(emissions) / gens.efficiency + em_pu = ( + weightings["generators"].to_frame("weightings") + @ em_pu.to_frame("weightings").T + ) + emitted = n.generators_t.p[gens.index].mul(em_pu) + + emitted_grouped = ( + emitted.groupby(level=0).sum().groupby(n.generators.carrier, axis=1).sum().T + ) + + df = df.reindex(emitted_grouped.index.union(df.index)) + + df.loc[emitted_grouped.index, label] = emitted_grouped.values + + if any(n.stores.carrier == "co2"): + co2_i = n.stores[n.stores.carrier == "co2"].index + df[label] = n.stores_t.e.groupby(level=0).last()[co2_i].iloc[:, 0] + + return df + + +outputs = [ + "nodal_costs", + "nodal_capacities", + "nodal_cfs", + "cfs", + "costs", + "capacities", + "curtailment", + "energy", + "supply", + "supply_energy", + "prices", + "weighted_prices", + "price_statistics", + "market_values", + "metrics", + "co2_emissions", +] + + +def make_summaries(networks_dict): + columns = pd.MultiIndex.from_tuples( + networks_dict.keys(), names=["cluster", "lv", "opt"] + ) + df = {} + + for output in outputs: + df[output] = pd.DataFrame(columns=columns, dtype=float) + + for label, filename in iteritems(networks_dict): + print(label, filename) + try: + n = pypsa.Network(filename) + except OSError: + print(label, " not solved yet.") + continue + # del networks_dict[label] + + if not hasattr(n, "objective"): + print(label, " not solved correctly. Check log if infeasible or unbounded.") + continue + assign_carriers(n) + assign_locations(n) + + for output in outputs: + df[output] = globals()["calculate_" + output](n, label, df[output]) + + return df + + +def to_csv(df): + for key in df: + df[key] = df[key].apply(lambda x: pd.to_numeric(x)) + df[key].to_csv(snakemake.output[key]) + + +# %% +if __name__ == "__main__": + # Detect running outside of snakemake and mock snakemake for testing + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + snakemake = mock_snakemake("make_summary_perfect") + + run = snakemake.config["run"]["name"] + if run != "": + run += "/" + + networks_dict = { + (clusters, lv, opts + sector_opts): "results/" + + run + + f"postnetworks/elec_s{simpl}_{clusters}_l{lv}_{opts}_{sector_opts}_brownfield_all_years.nc" + for simpl in snakemake.config["scenario"]["simpl"] + for clusters in snakemake.config["scenario"]["clusters"] + for opts in snakemake.config["scenario"]["opts"] + for sector_opts in snakemake.config["scenario"]["sector_opts"] + for lv in snakemake.config["scenario"]["ll"] + } + + print(networks_dict) + + nyears = 1 + costs_db = prepare_costs( + snakemake.input.costs, + snakemake.config["costs"], + nyears, + ) + + df = make_summaries(networks_dict) + + df["metrics"].loc["total costs"] = df["costs"].sum().groupby(level=[0, 1, 2]).sum() + + to_csv(df) diff --git a/scripts/plot_network.py b/scripts/plot_network.py index ae1d0e0a..c0b35411 100644 --- a/scripts/plot_network.py +++ b/scripts/plot_network.py @@ -24,7 +24,7 @@ from make_summary import assign_carriers from plot_summary import preferred_order, rename_techs from pypsa.plot import add_legend_circles, add_legend_lines, add_legend_patches -plt.style.use(["ggplot", "matplotlibrc"]) +plt.style.use(["ggplot"]) def rename_techs_tyndp(tech): @@ -913,6 +913,159 @@ def plot_series(network, carrier="AC", name="test"): ) +def plot_map_perfect( + network, + components=["Link", "Store", "StorageUnit", "Generator"], + bus_size_factor=1.7e10, +): + n = network.copy() + assign_location(n) + # Drop non-electric buses so they don't clutter the plot + n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True) + # investment periods + investments = n.snapshots.levels[0] + + costs = {} + for comp in components: + df_c = n.df(comp) + if df_c.empty: + continue + df_c["nice_group"] = df_c.carrier.map(rename_techs_tyndp) + + attr = "e_nom_opt" if comp == "Store" else "p_nom_opt" + + active = pd.concat( + [n.get_active_assets(comp, inv_p).rename(inv_p) for inv_p in investments], + axis=1, + ).astype(int) + capital_cost = n.df(comp)[attr] * n.df(comp).capital_cost + capital_cost_t = ( + (active.mul(capital_cost, axis=0)) + .groupby([n.df(comp).location, n.df(comp).nice_group]) + .sum() + ) + + capital_cost_t.drop("load", level=1, inplace=True, errors="ignore") + + costs[comp] = capital_cost_t + + costs = pd.concat(costs).groupby(level=[1, 2]).sum() + costs.drop(costs[costs.sum(axis=1) == 0].index, inplace=True) + + new_columns = preferred_order.intersection(costs.index.levels[1]).append( + costs.index.levels[1].difference(preferred_order) + ) + costs = costs.reindex(new_columns, level=1) + + for item in new_columns: + if item not in snakemake.config["plotting"]["tech_colors"]: + print( + "Warning!", + item, + "not in config/plotting/tech_colors, assign random color", + ) + snakemake.config["plotting"]["tech_colors"] = "pink" + + n.links.drop( + n.links.index[(n.links.carrier != "DC") & (n.links.carrier != "B2B")], + inplace=True, + ) + + # drop non-bus + to_drop = costs.index.levels[0].symmetric_difference(n.buses.index) + if len(to_drop) != 0: + print("dropping non-buses", to_drop) + costs.drop(to_drop, level=0, inplace=True, axis=0, errors="ignore") + + # make sure they are removed from index + costs.index = pd.MultiIndex.from_tuples(costs.index.values) + + # PDF has minimum width, so set these to zero + line_lower_threshold = 500.0 + line_upper_threshold = 1e4 + linewidth_factor = 2e3 + ac_color = "gray" + dc_color = "m" + + line_widths = n.lines.s_nom_opt + link_widths = n.links.p_nom_opt + linewidth_factor = 2e3 + line_lower_threshold = 0.0 + title = "Today's transmission" + + line_widths[line_widths < line_lower_threshold] = 0.0 + link_widths[link_widths < line_lower_threshold] = 0.0 + + line_widths[line_widths > line_upper_threshold] = line_upper_threshold + link_widths[link_widths > line_upper_threshold] = line_upper_threshold + + for year in costs.columns: + fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()}) + fig.set_size_inches(7, 6) + fig.suptitle(year) + + n.plot( + bus_sizes=costs[year] / bus_size_factor, + bus_colors=snakemake.config["plotting"]["tech_colors"], + line_colors=ac_color, + link_colors=dc_color, + line_widths=line_widths / linewidth_factor, + link_widths=link_widths / linewidth_factor, + ax=ax, + **map_opts, + ) + + sizes = [20, 10, 5] + labels = [f"{s} bEUR/a" for s in sizes] + sizes = [s / bus_size_factor * 1e9 for s in sizes] + + legend_kw = dict( + loc="upper left", + bbox_to_anchor=(0.01, 1.06), + labelspacing=0.8, + frameon=False, + handletextpad=0, + title="system cost", + ) + + add_legend_circles( + ax, + sizes, + labels, + srid=n.srid, + patch_kw=dict(facecolor="lightgrey"), + legend_kw=legend_kw, + ) + + sizes = [10, 5] + labels = [f"{s} GW" for s in sizes] + scale = 1e3 / linewidth_factor + sizes = [s * scale for s in sizes] + + legend_kw = dict( + loc="upper left", + bbox_to_anchor=(0.27, 1.06), + frameon=False, + labelspacing=0.8, + handletextpad=1, + title=title, + ) + + add_legend_lines( + ax, sizes, labels, patch_kw=dict(color="lightgrey"), legend_kw=legend_kw + ) + + legend_kw = dict( + bbox_to_anchor=(1.52, 1.04), + frameon=False, + ) + + fig.savefig( + snakemake.output[f"map_{year}"], transparent=True, bbox_inches="tight" + ) + + +# %% if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake @@ -921,10 +1074,9 @@ if __name__ == "__main__": "plot_network", simpl="", opts="", - clusters="5", - ll="v1.5", - sector_opts="CO2L0-1H-T-H-B-I-A-solar+p3-dist1", - planning_horizons="2030", + clusters="37", + ll="v1.0", + sector_opts="4380H-T-H-B-I-A-solar+p3-dist1", ) logging.basicConfig(level=snakemake.config["logging"]["level"]) @@ -938,16 +1090,23 @@ if __name__ == "__main__": if map_opts["boundaries"] is None: map_opts["boundaries"] = regions.total_bounds[[0, 2, 1, 3]] + [-1, 1, -1, 1] - plot_map( - n, - components=["generators", "links", "stores", "storage_units"], - bus_size_factor=2e10, - transmission=False, - ) + if snakemake.params["foresight"] == "perfect": + plot_map_perfect( + n, + components=["Link", "Store", "StorageUnit", "Generator"], + bus_size_factor=2e10, + ) + else: + plot_map( + n, + components=["generators", "links", "stores", "storage_units"], + bus_size_factor=2e10, + transmission=False, + ) - plot_h2_map(n, regions) - plot_ch4_map(n) - plot_map_without(n) + plot_h2_map(n, regions) + plot_ch4_map(n) + plot_map_without(n) # plot_series(n, carrier="AC", name=suffix) # plot_series(n, carrier="heat", name=suffix) diff --git a/scripts/plot_summary.py b/scripts/plot_summary.py index a501dcf7..86f40225 100644 --- a/scripts/plot_summary.py +++ b/scripts/plot_summary.py @@ -49,6 +49,10 @@ def rename_techs(label): # "H2 Fuel Cell": "hydrogen storage", # "H2 pipeline": "hydrogen storage", "battery": "battery storage", + "H2 for industry": "H2 for industry", + "land transport fuel cell": "land transport fuel cell", + "land transport oil": "land transport oil", + "oil shipping": "shipping oil", # "CC": "CC" } @@ -157,11 +161,11 @@ def plot_costs(): df.index.difference(preferred_order) ) - new_columns = df.sum().sort_values().index + # new_columns = df.sum().sort_values().index fig, ax = plt.subplots(figsize=(12, 8)) - df.loc[new_index, new_columns].T.plot( + df.loc[new_index].T.plot( kind="bar", ax=ax, stacked=True, @@ -213,17 +217,22 @@ def plot_energy(): logger.info(f"Total energy of {round(df.sum()[0])} TWh/a") + if df.empty: + fig, ax = plt.subplots(figsize=(12, 8)) + fig.savefig(snakemake.output.energy, bbox_inches="tight") + return + new_index = preferred_order.intersection(df.index).append( df.index.difference(preferred_order) ) - new_columns = df.columns.sort_values() + # new_columns = df.columns.sort_values() fig, ax = plt.subplots(figsize=(12, 8)) - logger.debug(df.loc[new_index, new_columns]) + logger.debug(df.loc[new_index]) - df.loc[new_index, new_columns].T.plot( + df.loc[new_index].T.plot( kind="bar", ax=ax, stacked=True, @@ -267,8 +276,6 @@ def plot_balances(): i for i in balances_df.index.levels[0] if i not in co2_carriers ] - fig, ax = plt.subplots(figsize=(12, 8)) - for k, v in balances.items(): df = balances_df.loc[v] df = df.groupby(df.index.get_level_values(2)).sum() @@ -279,7 +286,7 @@ def plot_balances(): # remove trailing link ports df.index = [ i[:-1] - if ((i not in ["co2", "NH3"]) and (i[-1:] in ["0", "1", "2", "3"])) + if ((i not in ["co2", "NH3", "H2"]) and (i[-1:] in ["0", "1", "2", "3"])) else i for i in df.index ] @@ -313,6 +320,8 @@ def plot_balances(): new_columns = df.columns.sort_values() + fig, ax = plt.subplots(figsize=(12, 8)) + df.loc[new_index, new_columns].T.plot( kind="bar", ax=ax, @@ -345,8 +354,6 @@ def plot_balances(): fig.savefig(snakemake.output.balances[:-10] + k + ".pdf", bbox_inches="tight") - plt.cla() - def historical_emissions(countries): """ @@ -354,8 +361,7 @@ def historical_emissions(countries): """ # https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-16 # downloaded 201228 (modified by EEA last on 201221) - fn = "data/bundle-sector/eea/UNFCCC_v23.csv" - df = pd.read_csv(fn, encoding="latin-1") + df = pd.read_csv(snakemake.input.co2, encoding="latin-1", low_memory=False) df.loc[df["Year"] == "1985-1987", "Year"] = 1986 df["Year"] = df["Year"].astype(int) df = df.set_index( @@ -379,18 +385,21 @@ def historical_emissions(countries): e["waste management"] = "5 - Waste management" e["other"] = "6 - Other Sector" e["indirect"] = "ind_CO2 - Indirect CO2" - e["total wL"] = "Total (with LULUCF)" - e["total woL"] = "Total (without LULUCF)" + e["other LULUCF"] = "4.H - Other LULUCF" pol = ["CO2"] # ["All greenhouse gases - (CO2 equivalent)"] if "GB" in countries: countries.remove("GB") countries.append("UK") - # remove countries which are not included in eea historical emission dataset - countries_to_remove = {"AL", "BA", "ME", "MK", "RS"} - countries = list(set(countries) - countries_to_remove) - year = np.arange(1990, 2018).tolist() + year = df.index.levels[0][df.index.levels[0] >= 1990] + + missing = pd.Index(countries).difference(df.index.levels[2]) + if not missing.empty: + logger.warning( + f"The following countries are missing and not considered when plotting historic CO2 emissions: {missing}" + ) + countries = pd.Index(df.index.levels[2]).intersection(countries) idx = pd.IndexSlice co2_totals = ( @@ -453,18 +462,12 @@ def plot_carbon_budget_distribution(input_eurostat): plt.rcParams["xtick.labelsize"] = 20 plt.rcParams["ytick.labelsize"] = 20 - plt.figure(figsize=(10, 7)) - gs1 = gridspec.GridSpec(1, 1) - ax1 = plt.subplot(gs1[0, 0]) - ax1.set_ylabel("CO$_2$ emissions (Gt per year)", fontsize=22) - ax1.set_ylim([0, 5]) - ax1.set_xlim([1990, snakemake.params.planning_horizons[-1] + 1]) - - path_cb = "results/" + snakemake.params.RDIR + "csvs/" - countries = snakemake.params.countries emissions_scope = snakemake.params.emissions_scope report_year = snakemake.params.eurostat_report_year input_co2 = snakemake.input.co2 + + # historic emissions + countries = snakemake.params.countries e_1990 = co2_emissions_year( countries, input_eurostat, @@ -474,15 +477,37 @@ def plot_carbon_budget_distribution(input_eurostat): input_co2, year=1990, ) - CO2_CAP = pd.read_csv(path_cb + "carbon_budget_distribution.csv", index_col=0) - - ax1.plot(e_1990 * CO2_CAP[o], linewidth=3, color="dodgerblue", label=None) - emissions = historical_emissions(countries) + # add other years https://sdi.eea.europa.eu/data/0569441f-2853-4664-a7cd-db969ef54de0 + emissions.loc[2019] = 2.971372 + emissions.loc[2020] = 2.691958 + emissions.loc[2021] = 2.869355 + + if snakemake.config["foresight"] == "myopic": + path_cb = "results/" + snakemake.params.RDIR + "/csvs/" + co2_cap = pd.read_csv(path_cb + "carbon_budget_distribution.csv", index_col=0)[ + ["cb"] + ] + co2_cap *= e_1990 + else: + supply_energy = pd.read_csv( + snakemake.input.balances, index_col=[0, 1, 2], header=[0, 1, 2, 3] + ) + co2_cap = ( + supply_energy.loc["co2"].droplevel(0).drop("co2").sum().unstack().T / 1e9 + ) + co2_cap.rename(index=lambda x: int(x), inplace=True) + + plt.figure(figsize=(10, 7)) + gs1 = gridspec.GridSpec(1, 1) + ax1 = plt.subplot(gs1[0, 0]) + ax1.set_ylabel("CO$_2$ emissions \n [Gt per year]", fontsize=22) + # ax1.set_ylim([0, 5]) + ax1.set_xlim([1990, snakemake.params.planning_horizons[-1] + 1]) ax1.plot(emissions, color="black", linewidth=3, label=None) - # plot committed and uder-discussion targets + # plot committed and under-discussion targets # (notice that historical emissions include all countries in the # network, but targets refer to EU) ax1.plot( @@ -499,7 +524,7 @@ def plot_carbon_budget_distribution(input_eurostat): [0.45 * emissions[1990]], marker="*", markersize=12, - markerfacecolor="white", + markerfacecolor="black", markeredgecolor="black", ) @@ -523,21 +548,7 @@ def plot_carbon_budget_distribution(input_eurostat): ax1.plot( [2050], - [0.01 * emissions[1990]], - marker="*", - markersize=12, - markerfacecolor="white", - linewidth=0, - markeredgecolor="black", - label="EU under-discussion target", - zorder=10, - clip_on=False, - ) - - ax1.plot( - [2050], - [0.125 * emissions[1990]], - "ro", + [0.0 * emissions[1990]], marker="*", markersize=12, markerfacecolor="black", @@ -545,14 +556,19 @@ def plot_carbon_budget_distribution(input_eurostat): label="EU committed target", ) + for col in co2_cap.columns: + ax1.plot(co2_cap[col], linewidth=3, label=col) + ax1.legend( fancybox=True, fontsize=18, loc=(0.01, 0.01), facecolor="white", frameon=True ) - path_cb_plot = "results/" + snakemake.params.RDIR + "graphs/" - plt.savefig(path_cb_plot + "carbon_budget_plot.pdf", dpi=300) + plt.grid(axis="y") + path = snakemake.output.balances.split("balances")[0] + "carbon_budget.pdf" + plt.savefig(path, bbox_inches="tight") +# %% if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake @@ -571,6 +587,7 @@ if __name__ == "__main__": for sector_opts in snakemake.params.sector_opts: opts = sector_opts.split("-") - for o in opts: - if "cb" in o: - plot_carbon_budget_distribution(snakemake.input.eurostat) + if any(["cb" in o for o in opts]) or ( + snakemake.config["foresight"] == "perfect" + ): + plot_carbon_budget_distribution(snakemake.input.eurostat) diff --git a/scripts/prepare_perfect_foresight.py b/scripts/prepare_perfect_foresight.py new file mode 100644 index 00000000..ca94dedc --- /dev/null +++ b/scripts/prepare_perfect_foresight.py @@ -0,0 +1,553 @@ +# -*- coding: utf-8 -*- +# SPDX-FileCopyrightText: : 2020-2023 The PyPSA-Eur Authors +# +# SPDX-License-Identifier: MIT +""" +Concats pypsa networks of single investment periods to one network. +""" + +import logging +import re + +import numpy as np +import pandas as pd +import pypsa +from _helpers import update_config_with_sector_opts +from add_existing_baseyear import add_build_year_to_new_assets +from pypsa.descriptors import expand_series +from pypsa.io import import_components_from_dataframe +from six import iterkeys + +logger = logging.getLogger(__name__) + + +# helper functions --------------------------------------------------- +def get_missing(df, n, c): + """ + Get in network n missing assets of df for component c. + + Input: + df: pandas DataFrame, static values of pypsa components + n : pypsa Network to which new assets should be added + c : string, pypsa component.list_name (e.g. "generators") + Return: + pd.DataFrame with static values of missing assets + """ + df_final = getattr(n, c) + missing_i = df.index.difference(df_final.index) + return df.loc[missing_i] + + +def get_social_discount(t, r=0.01): + """ + Calculate for a given time t and social discount rate r [per unit] the + social discount. + """ + return 1 / (1 + r) ** t + + +def get_investment_weighting(time_weighting, r=0.01): + """ + Define cost weighting. + + Returns cost weightings depending on the the time_weighting + (pd.Series) and the social discountrate r + """ + end = time_weighting.cumsum() + start = time_weighting.cumsum().shift().fillna(0) + return pd.concat([start, end], axis=1).apply( + lambda x: sum([get_social_discount(t, r) for t in range(int(x[0]), int(x[1]))]), + axis=1, + ) + + +def add_year_to_constraints(n, baseyear): + """ + Add investment period to global constraints and rename index. + + Parameters + ---------- + n : pypsa.Network + baseyear : int + year in which optimized assets are built + """ + + for c in n.iterate_components(["GlobalConstraint"]): + c.df["investment_period"] = baseyear + c.df.rename(index=lambda x: x + "-" + str(baseyear), inplace=True) + + +def hvdc_transport_model(n): + """ + Convert AC lines to DC links for multi-decade optimisation with line + expansion. + + Losses of DC links are assumed to be 3% per 1000km + """ + + logger.info("Convert AC lines to DC links to perform multi-decade optimisation.") + + n.madd( + "Link", + n.lines.index, + bus0=n.lines.bus0, + bus1=n.lines.bus1, + p_nom_extendable=True, + p_nom=n.lines.s_nom, + p_nom_min=n.lines.s_nom, + p_min_pu=-1, + efficiency=1 - 0.03 * n.lines.length / 1000, + marginal_cost=0, + carrier="DC", + length=n.lines.length, + capital_cost=n.lines.capital_cost, + ) + + # Remove AC lines + logger.info("Removing AC lines") + lines_rm = n.lines.index + n.mremove("Line", lines_rm) + + # Set efficiency of all DC links to include losses depending on length + n.links.loc[n.links.carrier == "DC", "efficiency"] = ( + 1 - 0.03 * n.links.loc[n.links.carrier == "DC", "length"] / 1000 + ) + + +def adjust_electricity_grid(n, year, years): + """ + Add carrier to lines. Replace AC lines with DC links in case of line + expansion. Add lifetime to DC links in case of line expansion. + + Parameters + ---------- + n : pypsa.Network + year : int + year in which optimized assets are built + years: list + investment periods + """ + n.lines["carrier"] = "AC" + links_i = n.links[n.links.carrier == "DC"].index + if n.lines.s_nom_extendable.any() or n.links.loc[links_i, "p_nom_extendable"].any(): + hvdc_transport_model(n) + links_i = n.links[n.links.carrier == "DC"].index + n.links.loc[links_i, "lifetime"] = 100 + if year != years[0]: + n.links.loc[links_i, "p_nom_min"] = 0 + n.links.loc[links_i, "p_nom"] = 0 + + +# -------------------------------------------------------------------- +def concat_networks(years): + """ + Concat given pypsa networks and adds build_year. + + Return: + n : pypsa.Network for the whole planning horizon + """ + + # input paths of sector coupling networks + network_paths = [snakemake.input.brownfield_network] + [ + snakemake.input[f"network_{year}"] for year in years[1:] + ] + # final concatenated network + n = pypsa.Network() + + # iterate over single year networks and concat to perfect foresight network + for i, network_path in enumerate(network_paths): + year = years[i] + network = pypsa.Network(network_path) + adjust_electricity_grid(network, year, years) + add_build_year_to_new_assets(network, year) + + # static ---------------------------------- + # (1) add buses and carriers + for component in network.iterate_components(["Bus", "Carrier"]): + df_year = component.df + # get missing assets + missing = get_missing(df_year, n, component.list_name) + import_components_from_dataframe(n, missing, component.name) + # (2) add generators, links, stores and loads + for component in network.iterate_components( + ["Generator", "Link", "Store", "Load", "Line", "StorageUnit"] + ): + df_year = component.df.copy() + missing = get_missing(df_year, n, component.list_name) + + import_components_from_dataframe(n, missing, component.name) + + # time variant -------------------------------------------------- + network_sns = pd.MultiIndex.from_product([[year], network.snapshots]) + snapshots = n.snapshots.drop("now", errors="ignore").union(network_sns) + n.set_snapshots(snapshots) + + for component in network.iterate_components(): + pnl = getattr(n, component.list_name + "_t") + for k in iterkeys(component.pnl): + pnl_year = component.pnl[k].copy().reindex(snapshots, level=1) + if pnl_year.empty and ~(component.name == "Load" and k == "p_set"): + continue + if component.name == "Load": + static_load = network.loads.loc[network.loads.p_set != 0] + static_load_t = expand_series(static_load.p_set, network_sns).T + pnl_year = pd.concat( + [pnl_year.reindex(network_sns), static_load_t], axis=1 + ) + columns = (pnl[k].columns.union(pnl_year.columns)).unique() + pnl[k] = pnl[k].reindex(columns=columns) + pnl[k].loc[pnl_year.index, pnl_year.columns] = pnl_year + + else: + # this is to avoid adding multiple times assets with + # infinite lifetime as ror + cols = pnl_year.columns.difference(pnl[k].columns) + pnl[k] = pd.concat([pnl[k], pnl_year[cols]], axis=1) + + n.snapshot_weightings.loc[year, :] = network.snapshot_weightings.values + + # (3) global constraints + for component in network.iterate_components(["GlobalConstraint"]): + add_year_to_constraints(network, year) + import_components_from_dataframe(n, component.df, component.name) + + # set investment periods + n.investment_periods = n.snapshots.levels[0] + # weighting of the investment period -> assuming last period same weighting as the period before + time_w = n.investment_periods.to_series().diff().shift(-1).fillna(method="ffill") + n.investment_period_weightings["years"] = time_w + # set objective weightings + objective_w = get_investment_weighting( + n.investment_period_weightings["years"], social_discountrate + ) + n.investment_period_weightings["objective"] = objective_w + # all former static loads are now time-dependent -> set static = 0 + n.loads["p_set"] = 0 + n.loads_t.p_set.fillna(0, inplace=True) + + return n + + +def adjust_stores(n): + """ + Make sure that stores still behave cyclic over one year and not whole + modelling horizon. + """ + # cyclic constraint + cyclic_i = n.stores[n.stores.e_cyclic].index + n.stores.loc[cyclic_i, "e_cyclic_per_period"] = True + n.stores.loc[cyclic_i, "e_cyclic"] = False + # non cyclic store assumptions + non_cyclic_store = ["co2", "co2 stored", "solid biomass", "biogas", "Li ion"] + co2_i = n.stores[n.stores.carrier.isin(non_cyclic_store)].index + n.stores.loc[co2_i, "e_cyclic_per_period"] = False + n.stores.loc[co2_i, "e_cyclic"] = False + # e_initial at beginning of each investment period + e_initial_store = ["solid biomass", "biogas"] + co2_i = n.stores[n.stores.carrier.isin(e_initial_store)].index + n.stores.loc[co2_i, "e_initial_per_period"] = True + # n.stores.loc[co2_i, "e_initial"] *= 10 + # n.stores.loc[co2_i, "e_nom"] *= 10 + e_initial_store = ["co2 stored"] + co2_i = n.stores[n.stores.carrier.isin(e_initial_store)].index + n.stores.loc[co2_i, "e_initial_per_period"] = True + + return n + + +def set_phase_out(n, carrier, ct, phase_out_year): + """ + Set planned phase outs for given carrier,country (ct) and planned year of + phase out (phase_out_year). + """ + df = n.links[(n.links.carrier.isin(carrier)) & (n.links.bus1.str[:2] == ct)] + # assets which are going to be phased out before end of their lifetime + assets_i = df[df[["build_year", "lifetime"]].sum(axis=1) > phase_out_year].index + build_year = n.links.loc[assets_i, "build_year"] + # adjust lifetime + n.links.loc[assets_i, "lifetime"] = (phase_out_year - build_year).astype(float) + + +def set_all_phase_outs(n): + # TODO move this to a csv or to the config + planned = [ + (["nuclear"], "DE", 2022), + (["nuclear"], "BE", 2025), + (["nuclear"], "ES", 2027), + (["coal", "lignite"], "DE", 2030), + (["coal", "lignite"], "ES", 2027), + (["coal", "lignite"], "FR", 2022), + (["coal", "lignite"], "GB", 2024), + (["coal", "lignite"], "IT", 2025), + (["coal", "lignite"], "DK", 2030), + (["coal", "lignite"], "FI", 2030), + (["coal", "lignite"], "HU", 2030), + (["coal", "lignite"], "SK", 2030), + (["coal", "lignite"], "GR", 2030), + (["coal", "lignite"], "IE", 2030), + (["coal", "lignite"], "NL", 2030), + (["coal", "lignite"], "RS", 2030), + ] + for carrier, ct, phase_out_year in planned: + set_phase_out(n, carrier, ct, phase_out_year) + # remove assets which are already phased out + remove_i = n.links[n.links[["build_year", "lifetime"]].sum(axis=1) < years[0]].index + n.mremove("Link", remove_i) + + +def set_carbon_constraints(n, opts): + """ + Add global constraints for carbon emissions. + """ + budget = None + for o in opts: + # other budgets + m = re.match(r"^\d+p\d$", o, re.IGNORECASE) + if m is not None: + budget = snakemake.config["co2_budget"][m.group(0)] * 1e9 + if budget != None: + logger.info("add carbon budget of {}".format(budget)) + n.add( + "GlobalConstraint", + "Budget", + type="Co2Budget", + carrier_attribute="co2_emissions", + sense="<=", + constant=budget, + investment_period=n.investment_periods[-1], + ) + + # drop other CO2 limits + drop_i = n.global_constraints[n.global_constraints.type == "co2_limit"].index + n.mremove("GlobalConstraint", drop_i) + + n.add( + "GlobalConstraint", + "carbon_neutral", + type="co2_limit", + carrier_attribute="co2_emissions", + sense="<=", + constant=0, + investment_period=n.investment_periods[-1], + ) + + # set minimum CO2 emission constraint to avoid too fast reduction + if "co2min" in opts: + emissions_1990 = 4.53693 + emissions_2019 = 3.344096 + target_2030 = 0.45 * emissions_1990 + annual_reduction = (emissions_2019 - target_2030) / 11 + first_year = n.snapshots.levels[0][0] + time_weightings = n.investment_period_weightings.loc[first_year, "years"] + co2min = emissions_2019 - ((first_year - 2019) * annual_reduction) + logger.info( + "add minimum emissions for {} of {} t CO2/a".format(first_year, co2min) + ) + n.add( + "GlobalConstraint", + f"Co2Min-{first_year}", + type="Co2min", + carrier_attribute="co2_emissions", + sense=">=", + investment_period=first_year, + constant=co2min * 1e9 * time_weightings, + ) + + return n + + +def adjust_lvlimit(n): + """ + Convert global constraints for single investment period to one uniform if + all attributes stay the same. + """ + c = "GlobalConstraint" + cols = ["carrier_attribute", "sense", "constant", "type"] + glc_type = "transmission_volume_expansion_limit" + if (n.df(c)[n.df(c).type == glc_type][cols].nunique() == 1).all(): + glc = n.df(c)[n.df(c).type == glc_type][cols].iloc[[0]] + glc.index = pd.Index(["lv_limit"]) + remove_i = n.df(c)[n.df(c).type == glc_type].index + n.mremove(c, remove_i) + import_components_from_dataframe(n, glc, c) + + return n + + +def adjust_CO2_glc(n): + c = "GlobalConstraint" + glc_name = "CO2Limit" + glc_type = "primary_energy" + mask = (n.df(c).index.str.contains(glc_name)) & (n.df(c).type == glc_type) + n.df(c).loc[mask, "type"] = "co2_limit" + + return n + + +def add_H2_boilers(n): + """ + Gas boilers can be retrofitted to run with H2. + + Add H2 boilers for heating for all existing gas boilers. + """ + c = "Link" + logger.info("Add H2 boilers.") + # existing gas boilers + mask = n.links.carrier.str.contains("gas boiler") & ~n.links.p_nom_extendable + gas_i = n.links[mask].index + df = n.links.loc[gas_i] + # adjust bus 0 + df["bus0"] = df.bus1.map(n.buses.location) + " H2" + # rename carrier and index + df["carrier"] = df.carrier.apply( + lambda x: x.replace("gas boiler", "retrofitted H2 boiler") + ) + df.rename( + index=lambda x: x.replace("gas boiler", "retrofitted H2 boiler"), inplace=True + ) + # todo, costs for retrofitting + df["capital_costs"] = 100 + # set existing capacity to zero + df["p_nom"] = 0 + df["p_nom_extendable"] = True + # add H2 boilers to network + import_components_from_dataframe(n, df, c) + + +def apply_time_segmentation_perfect( + n, segments, solver_name="cbc", overwrite_time_dependent=True +): + """ + Aggregating time series to segments with different lengths. + + Input: + n: pypsa Network + segments: (int) number of segments in which the typical period should be + subdivided + solver_name: (str) name of solver + overwrite_time_dependent: (bool) overwrite time dependent data of pypsa network + with typical time series created by tsam + """ + try: + import tsam.timeseriesaggregation as tsam + except: + raise ModuleNotFoundError( + "Optional dependency 'tsam' not found." "Install via 'pip install tsam'" + ) + + # get all time-dependent data + columns = pd.MultiIndex.from_tuples([], names=["component", "key", "asset"]) + raw = pd.DataFrame(index=n.snapshots, columns=columns) + for c in n.iterate_components(): + for attr, pnl in c.pnl.items(): + # exclude e_min_pu which is used for SOC of EVs in the morning + if not pnl.empty and attr != "e_min_pu": + df = pnl.copy() + df.columns = pd.MultiIndex.from_product([[c.name], [attr], df.columns]) + raw = pd.concat([raw, df], axis=1) + raw = raw.dropna(axis=1) + sn_weightings = {} + + for year in raw.index.levels[0]: + logger.info(f"Find representative snapshots for {year}.") + raw_t = raw.loc[year] + # normalise all time-dependent data + annual_max = raw_t.max().replace(0, 1) + raw_t = raw_t.div(annual_max, level=0) + # get representative segments + agg = tsam.TimeSeriesAggregation( + raw_t, + hoursPerPeriod=len(raw_t), + noTypicalPeriods=1, + noSegments=int(segments), + segmentation=True, + solver=solver_name, + ) + segmented = agg.createTypicalPeriods() + + weightings = segmented.index.get_level_values("Segment Duration") + offsets = np.insert(np.cumsum(weightings[:-1]), 0, 0) + timesteps = [raw_t.index[0] + pd.Timedelta(f"{offset}h") for offset in offsets] + snapshots = pd.DatetimeIndex(timesteps) + sn_weightings[year] = pd.Series( + weightings, index=snapshots, name="weightings", dtype="float64" + ) + + sn_weightings = pd.concat(sn_weightings) + n.set_snapshots(sn_weightings.index) + n.snapshot_weightings = n.snapshot_weightings.mul(sn_weightings, axis=0) + + return n + + +def set_temporal_aggregation_SEG(n, opts, solver_name): + """ + Aggregate network temporally with tsam. + """ + for o in opts: + # segments with package tsam + m = re.match(r"^(\d+)seg$", o, re.IGNORECASE) + if m is not None: + segments = int(m[1]) + logger.info(f"Use temporal segmentation with {segments} segments") + n = apply_time_segmentation_perfect(n, segments, solver_name=solver_name) + break + return n + + +# %% +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + snakemake = mock_snakemake( + "prepare_perfect_foresight", + simpl="", + opts="", + clusters="37", + ll="v1.5", + sector_opts="1p7-4380H-T-H-B-I-A-solar+p3-dist1", + ) + + update_config_with_sector_opts(snakemake.config, snakemake.wildcards.sector_opts) + # parameters ----------------------------------------------------------- + years = snakemake.config["scenario"]["planning_horizons"] + opts = snakemake.wildcards.sector_opts.split("-") + social_discountrate = snakemake.config["costs"]["social_discountrate"] + for o in opts: + if "sdr" in o: + social_discountrate = float(o.replace("sdr", "")) / 100 + + logger.info( + "Concat networks of investment period {} with social discount rate of {}%".format( + years, social_discountrate * 100 + ) + ) + + # concat prenetworks of planning horizon to single network ------------ + n = concat_networks(years) + + # temporal aggregate + opts = snakemake.wildcards.sector_opts.split("-") + solver_name = snakemake.config["solving"]["solver"]["name"] + n = set_temporal_aggregation_SEG(n, opts, solver_name) + + # adjust global constraints lv limit if the same for all years + n = adjust_lvlimit(n) + # adjust global constraints CO2 limit + n = adjust_CO2_glc(n) + # adjust stores to multi period investment + n = adjust_stores(n) + + # set phase outs + set_all_phase_outs(n) + + # add H2 boiler + add_H2_boilers(n) + + # set carbon constraints + opts = snakemake.wildcards.sector_opts.split("-") + n = set_carbon_constraints(n, opts) + + # export network + n.export_to_netcdf(snakemake.output[0]) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index cd5d9570..998f954e 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -220,7 +220,7 @@ def co2_emissions_year( # TODO: move to own rule with sector-opts wildcard? -def build_carbon_budget(o, input_eurostat, fn, emissions_scope, report_year, input_co2): +def build_carbon_budget(o, input_eurostat, fn, emissions_scope, report_year): """ Distribute carbon budget following beta or exponential transition path. """ @@ -577,6 +577,7 @@ def add_co2_tracking(n, options): capital_cost=options["co2_sequestration_cost"], carrier="co2 stored", bus=spatial.co2.nodes, + lifetime=options["co2_sequestration_lifetime"], ) n.add("Carrier", "co2 stored") @@ -3374,7 +3375,7 @@ if __name__ == "__main__": spatial = define_spatial(pop_layout.index, options) - if snakemake.params.foresight == "myopic": + if snakemake.params.foresight in ["myopic", "perfect"]: add_lifetime_wind_solar(n, costs) conventional = snakemake.params.conventional_carriers @@ -3451,7 +3452,7 @@ if __name__ == "__main__": if "cb" not in o: continue limit_type = "carbon budget" - fn = "results/" + snakemake.params.RDIR + "csvs/carbon_budget_distribution.csv" + fn = "results/" + snakemake.params.RDIR + "/csvs/carbon_budget_distribution.csv" if not os.path.exists(fn): emissions_scope = snakemake.params.emissions_scope report_year = snakemake.params.eurostat_report_year @@ -3507,7 +3508,7 @@ if __name__ == "__main__": ) n.mremove("Line", idx) - first_year_myopic = (snakemake.params.foresight == "myopic") and ( + first_year_myopic = (snakemake.params.foresight in ["myopic", "perfect"]) and ( snakemake.params.planning_horizons[0] == investment_year ) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 5e8c0356..83281284 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -33,7 +33,9 @@ import numpy as np import pandas as pd import pypsa import xarray as xr +from _benchmark import memory_logger from _helpers import configure_logging, update_config_with_sector_opts +from pypsa.descriptors import get_activity_mask logger = logging.getLogger(__name__) pypsa.pf.logger.setLevel(logging.WARNING) @@ -47,6 +49,69 @@ def add_land_use_constraint(n, planning_horizons, config): _add_land_use_constraint(n) +def add_land_use_constraint_perfect(n): + """ + Add global constraints for tech capacity limit. + """ + logger.info("Add land-use constraint for perfect foresight") + + def compress_series(s): + def process_group(group): + if group.nunique() == 1: + return pd.Series(group.iloc[0], index=[None]) + else: + return group + + return s.groupby(level=[0, 1]).apply(process_group) + + def new_index_name(t): + # Convert all elements to string and filter out None values + parts = [str(x) for x in t if x is not None] + # Join with space, but use a dash for the last item if not None + return " ".join(parts[:2]) + (f"-{parts[-1]}" if len(parts) > 2 else "") + + def check_p_min_p_max(p_nom_max): + p_nom_min = n.generators[ext_i].groupby(grouper).sum().p_nom_min + p_nom_min = p_nom_min.reindex(p_nom_max.index) + check = ( + p_nom_min.groupby(level=[0, 1]).sum() + > p_nom_max.groupby(level=[0, 1]).min() + ) + if check.sum(): + logger.warning( + f"summed p_min_pu values at node larger than technical potential {check[check].index}" + ) + + grouper = [n.generators.carrier, n.generators.bus, n.generators.build_year] + ext_i = n.generators.p_nom_extendable + # get technical limit per node and investment period + p_nom_max = n.generators[ext_i].groupby(grouper).min().p_nom_max + # drop carriers without tech limit + p_nom_max = p_nom_max[~p_nom_max.isin([np.inf, np.nan])] + # carrier + carriers = p_nom_max.index.get_level_values(0).unique() + gen_i = n.generators[(n.generators.carrier.isin(carriers)) & (ext_i)].index + n.generators.loc[gen_i, "p_nom_min"] = 0 + # check minimum capacities + check_p_min_p_max(p_nom_max) + # drop multi entries in case p_nom_max stays constant in different periods + # p_nom_max = compress_series(p_nom_max) + # adjust name to fit syntax of nominal constraint per bus + df = p_nom_max.reset_index() + df["name"] = df.apply( + lambda row: f"nom_max_{row['carrier']}" + + (f"_{row['build_year']}" if row["build_year"] is not None else ""), + axis=1, + ) + + for name in df.name.unique(): + df_carrier = df[df.name == name] + bus = df_carrier.bus + n.buses.loc[bus, name] = df_carrier.p_nom_max.values + + return n + + def _add_land_use_constraint(n): # warning: this will miss existing offwind which is not classed AC-DC and has carrier 'offwind' @@ -82,7 +147,6 @@ def _add_land_use_constraint(n): def _add_land_use_constraint_m(n, planning_horizons, config): # if generators clustering is lower than network clustering, land_use accounting is at generators clusters - planning_horizons = param["planning_horizons"] grouping_years = config["existing_capacities"]["grouping_years"] current_horizon = snakemake.wildcards.planning_horizons @@ -116,7 +180,7 @@ def _add_land_use_constraint_m(n, planning_horizons, config): n.generators.p_nom_max.clip(lower=0, inplace=True) -def add_co2_sequestration_limit(n, limit=200): +def add_co2_sequestration_limit(n, config, limit=200): """ Add a global constraint on the amount of Mt CO2 that can be sequestered. """ @@ -130,16 +194,146 @@ def add_co2_sequestration_limit(n, limit=200): limit = float(o[o.find("seq") + 3 :]) * 1e6 break - n.add( + if not n.investment_periods.empty: + periods = n.investment_periods + names = pd.Index([f"co2_sequestration_limit-{period}" for period in periods]) + else: + periods = [np.nan] + names = pd.Index(["co2_sequestration_limit"]) + + n.madd( "GlobalConstraint", - "co2_sequestration_limit", + names, sense="<=", constant=limit, type="primary_energy", carrier_attribute="co2_absorptions", + investment_period=periods, ) +def add_carbon_constraint(n, snapshots): + glcs = n.global_constraints.query('type == "co2_limit"') + if glcs.empty: + return + for name, glc in glcs.iterrows(): + rhs = glc.constant + carattr = glc.carrier_attribute + emissions = n.carriers.query(f"{carattr} != 0")[carattr] + + if emissions.empty: + continue + + # stores + n.stores["carrier"] = n.stores.bus.map(n.buses.carrier) + stores = n.stores.query("carrier in @emissions.index and not e_cyclic") + time_valid = int(glc.loc["investment_period"]) + if not stores.empty: + last = n.snapshot_weightings.reset_index().groupby("period").last() + last_i = last.set_index([last.index, last.timestep]).index + final_e = n.model["Store-e"].loc[last_i, stores.index] + time_i = pd.IndexSlice[time_valid, :] + lhs = final_e.loc[time_i, :] - final_e.shift(snapshot=1).loc[time_i, :] + + n.model.add_constraints(lhs <= rhs, name=f"GlobalConstraint-{name}") + + +def add_carbon_budget_constraint(n, snapshots): + glcs = n.global_constraints.query('type == "Co2Budget"') + if glcs.empty: + return + for name, glc in glcs.iterrows(): + rhs = glc.constant + carattr = glc.carrier_attribute + emissions = n.carriers.query(f"{carattr} != 0")[carattr] + + if emissions.empty: + continue + + # stores + n.stores["carrier"] = n.stores.bus.map(n.buses.carrier) + stores = n.stores.query("carrier in @emissions.index and not e_cyclic") + time_valid = int(glc.loc["investment_period"]) + weighting = n.investment_period_weightings.loc[time_valid, "years"] + if not stores.empty: + last = n.snapshot_weightings.reset_index().groupby("period").last() + last_i = last.set_index([last.index, last.timestep]).index + final_e = n.model["Store-e"].loc[last_i, stores.index] + time_i = pd.IndexSlice[time_valid, :] + lhs = final_e.loc[time_i, :] * weighting + + n.model.add_constraints(lhs <= rhs, name=f"GlobalConstraint-{name}") + + +def add_max_growth(n, config): + """ + Add maximum growth rates for different carriers. + """ + + opts = snakemake.params["sector"]["limit_max_growth"] + # take maximum yearly difference between investment periods since historic growth is per year + factor = n.investment_period_weightings.years.max() * opts["factor"] + for carrier in opts["max_growth"].keys(): + max_per_period = opts["max_growth"][carrier] * factor + logger.info( + f"set maximum growth rate per investment period of {carrier} to {max_per_period} GW." + ) + n.carriers.loc[carrier, "max_growth"] = max_per_period * 1e3 + + for carrier in opts["max_relative_growth"].keys(): + max_r_per_period = opts["max_relative_growth"][carrier] + logger.info( + f"set maximum relative growth per investment period of {carrier} to {max_r_per_period}." + ) + n.carriers.loc[carrier, "max_relative_growth"] = max_r_per_period + + return n + + +def add_retrofit_gas_boiler_constraint(n, snapshots): + """ + Allow retrofitting of existing gas boilers to H2 boilers. + """ + c = "Link" + logger.info("Add constraint for retrofitting gas boilers to H2 boilers.") + # existing gas boilers + mask = n.links.carrier.str.contains("gas boiler") & ~n.links.p_nom_extendable + gas_i = n.links[mask].index + mask = n.links.carrier.str.contains("retrofitted H2 boiler") + h2_i = n.links[mask].index + + n.links.loc[gas_i, "p_nom_extendable"] = True + p_nom = n.links.loc[gas_i, "p_nom"] + n.links.loc[gas_i, "p_nom"] = 0 + + # heat profile + cols = n.loads_t.p_set.columns[ + n.loads_t.p_set.columns.str.contains("heat") + & ~n.loads_t.p_set.columns.str.contains("industry") + & ~n.loads_t.p_set.columns.str.contains("agriculture") + ] + profile = n.loads_t.p_set[cols].div( + n.loads_t.p_set[cols].groupby(level=0).max(), level=0 + ) + # to deal if max value is zero + profile.fillna(0, inplace=True) + profile.rename(columns=n.loads.bus.to_dict(), inplace=True) + profile = profile.reindex(columns=n.links.loc[gas_i, "bus1"]) + profile.columns = gas_i + + rhs = profile.mul(p_nom) + + dispatch = n.model["Link-p"] + active = get_activity_mask(n, c, snapshots, gas_i) + rhs = rhs[active] + p_gas = dispatch.sel(Link=gas_i) + p_h2 = dispatch.sel(Link=h2_i) + + lhs = p_gas + p_h2 + + n.model.add_constraints(lhs == rhs, name="gas_retrofit") + + def prepare_network( n, solve_opts=None, @@ -200,9 +394,14 @@ def prepare_network( if foresight == "myopic": add_land_use_constraint(n, planning_horizons, config) + if foresight == "perfect": + n = add_land_use_constraint_perfect(n) + if snakemake.params["sector"]["limit_max_growth"]["enable"]: + n = add_max_growth(n, config) + if n.stores.carrier.eq("co2 stored").any(): limit = co2_sequestration_potential - add_co2_sequestration_limit(n, limit=limit) + add_co2_sequestration_limit(n, config, limit=limit) return n @@ -612,12 +811,19 @@ def extra_functionality(n, snapshots): add_battery_constraints(n) add_lossy_bidirectional_link_constraints(n) add_pipe_retrofit_constraint(n) + if n._multi_invest: + add_carbon_constraint(n, snapshots) + add_carbon_budget_constraint(n, snapshots) + add_retrofit_gas_boiler_constraint(n, snapshots) def solve_network(n, config, solving, opts="", **kwargs): set_of_options = solving["solver"]["options"] cf_solving = solving["options"] + kwargs["multi_investment_periods"] = ( + True if config["foresight"] == "perfect" else False + ) kwargs["solver_options"] = ( solving["solver_options"][set_of_options] if set_of_options else {} ) @@ -664,18 +870,20 @@ def solve_network(n, config, solving, opts="", **kwargs): return n +# %% if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake snakemake = mock_snakemake( - "solve_network", + "solve_sector_network_perfect", + configfiles="../config/test/config.perfect.yaml", simpl="", - opts="Ept", - clusters="37", - ll="v1.0", - sector_opts="", - planning_horizons="2020", + opts="", + clusters="5", + ll="v1.5", + sector_opts="8760H-T-H-B-I-A-solar+p3-dist1", + planning_horizons="2030", ) configure_logging(snakemake) if "sector_opts" in snakemake.wildcards.keys(): @@ -702,13 +910,18 @@ if __name__ == "__main__": co2_sequestration_potential=snakemake.params["co2_sequestration_potential"], ) - n = solve_network( - n, - config=snakemake.config, - solving=snakemake.params.solving, - opts=opts, - log_fn=snakemake.log.solver, - ) + with memory_logger( + filename=getattr(snakemake.log, "memory", None), interval=30.0 + ) as mem: + n = solve_network( + n, + config=snakemake.config, + solving=snakemake.params.solving, + opts=opts, + log_fn=snakemake.log.solver, + ) + + logger.info("Maximum memory usage: {}".format(mem.mem_usage)) n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) n.export_to_netcdf(snakemake.output[0])