diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 7b9009c30..78e70b579 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -50,7 +50,7 @@ repos: - id: blackdoc # Formatting with "black" coding style -- repo: https://github.com/psf/black +- repo: https://github.com/psf/black-pre-commit-mirror rev: 23.12.1 hooks: # Format Python files diff --git a/Snakefile b/Snakefile index 83530df7f..14ce0e405 100644 --- a/Snakefile +++ b/Snakefile @@ -14,7 +14,7 @@ from snakemake.utils import min_version min_version("7.7") -if not exists("config/config.yaml"): +if not exists("config/config.yaml") and exists("config/config.default.yaml"): copyfile("config/config.default.yaml", "config/config.yaml") diff --git a/config/config.default.yaml b/config/config.default.yaml index b8945f754..1598ca1cd 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -158,6 +158,7 @@ renewable: resource: method: wind turbine: Vestas_V112_3MW + add_cutout_windspeed: true capacity_per_sqkm: 3 # correction_factor: 0.93 corine: @@ -166,13 +167,13 @@ renewable: distance_grid_codes: [1, 2, 3, 4, 5, 6] natura: true excluder_resolution: 100 - potential: simple # or conservative clip_p_max_pu: 1.e-2 offwind-ac: cutout: europe-2013-era5 resource: method: wind - turbine: NREL_ReferenceTurbine_5MW_offshore + turbine: NREL_ReferenceTurbine_2020ATB_5.5MW + add_cutout_windspeed: true capacity_per_sqkm: 2 correction_factor: 0.8855 corine: [44, 255] @@ -181,13 +182,13 @@ renewable: max_depth: 50 max_shore_distance: 30000 excluder_resolution: 200 - potential: simple # or conservative clip_p_max_pu: 1.e-2 offwind-dc: cutout: europe-2013-era5 resource: method: wind - turbine: NREL_ReferenceTurbine_5MW_offshore + turbine: NREL_ReferenceTurbine_2020ATB_5.5MW + add_cutout_windspeed: true capacity_per_sqkm: 2 correction_factor: 0.8855 corine: [44, 255] @@ -196,7 +197,6 @@ renewable: max_depth: 50 min_shore_distance: 30000 excluder_resolution: 200 - potential: simple # or conservative clip_p_max_pu: 1.e-2 solar: cutout: europe-2013-sarah @@ -211,7 +211,6 @@ renewable: corine: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 26, 31, 32] natura: true excluder_resolution: 100 - potential: simple # or conservative clip_p_max_pu: 1.e-2 hydro: cutout: europe-2013-era5 @@ -448,7 +447,6 @@ sector: solar_cf_correction: 0.788457 # = >>> 1/1.2683 marginal_cost_storage: 0. #1e-4 methanation: true - helmeth: false coal_cc: false dac: true co2_vent: false @@ -477,14 +475,28 @@ sector: - nearshore # within 50 km of sea # - offshore ammonia: false - min_part_load_fischer_tropsch: 0.9 - min_part_load_methanolisation: 0.5 + min_part_load_fischer_tropsch: 0.7 + min_part_load_methanolisation: 0.3 + min_part_load_methanation: 0.3 use_fischer_tropsch_waste_heat: true + use_haber_bosch_waste_heat: true + use_methanolisation_waste_heat: true + use_methanation_waste_heat: true use_fuel_cell_waste_heat: true - use_electrolysis_waste_heat: false + use_electrolysis_waste_heat: true electricity_distribution_grid: true electricity_distribution_grid_cost_factor: 1.0 electricity_grid_connection: true + transmission_efficiency: + DC: + efficiency_static: 0.98 + efficiency_per_1000km: 0.977 + H2 pipeline: + efficiency_per_1000km: 1 # 0.979 + compression_per_1000km: 0.019 + gas pipeline: + efficiency_per_1000km: 1 #0.977 + compression_per_1000km: 0.01 H2_network: true gas_network: false H2_retrofit: false @@ -494,6 +506,7 @@ sector: gas_distribution_grid_cost_factor: 1.0 biomass_spatial: false biomass_transport: false + biogas_upgrading_cc: false conventional_generation: OCGT: gas biomass_to_liquid: false @@ -544,8 +557,8 @@ industry: MWh_NH3_per_tNH3: 5.166 MWh_CH4_per_tNH3_SMR: 10.8 MWh_elec_per_tNH3_SMR: 0.7 - MWh_H2_per_tNH3_electrolysis: 6.5 - MWh_elec_per_tNH3_electrolysis: 1.17 + MWh_H2_per_tNH3_electrolysis: 5.93 + MWh_elec_per_tNH3_electrolysis: 0.2473 MWh_NH3_per_MWh_H2_cracker: 1.46 # https://github.com/euronion/trace/blob/44a5ff8401762edbef80eff9cfe5a47c8d3c8be4/data/efficiencies.csv NH3_process_emissions: 24.5 petrochemical_process_emissions: 25.5 @@ -776,6 +789,7 @@ plotting: fossil gas: '#e05b09' natural gas: '#e05b09' biogas to gas: '#e36311' + biogas to gas CC: '#e51245' CCGT: '#a85522' CCGT marginal: '#a85522' allam: '#B98F76' @@ -877,6 +891,7 @@ plotting: # heat demand Heat load: '#cc1f1f' heat: '#cc1f1f' + heat vent: '#aa3344' heat demand: '#cc1f1f' rural heat: '#ff5c5c' residential rural heat: '#ff7c7c' @@ -946,7 +961,6 @@ plotting: Sabatier: '#9850ad' methanation: '#c44ce6' methane: '#c44ce6' - helmeth: '#e899ff' # synfuels Fischer-Tropsch: '#25c49a' liquid: '#25c49a' diff --git a/doc/configtables/offwind-ac.csv b/doc/configtables/offwind-ac.csv index 6b756799c..c3512a9e7 100644 --- a/doc/configtables/offwind-ac.csv +++ b/doc/configtables/offwind-ac.csv @@ -12,5 +12,4 @@ ship_threshold,--,float,"Ship density threshold from which areas are excluded." max_depth,m,float,"Maximum sea water depth at which wind turbines can be build. Maritime areas with deeper waters are excluded in the process of calculating the AC-connected offshore wind potential." min_shore_distance,m,float,"Minimum distance to the shore below which wind turbines cannot be build. Such areas close to the shore are excluded in the process of calculating the AC-connected offshore wind potential." max_shore_distance,m,float,"Maximum distance to the shore above which wind turbines cannot be build. Such areas close to the shore are excluded in the process of calculating the AC-connected offshore wind potential." -potential,--,"One of {'simple', 'conservative'}","Method to compute the maximal installable potential for a node; confer :ref:`renewableprofiles`" clip_p_max_pu,p.u.,float,"To avoid too small values in the renewables` per-unit availability time series values below this threshold are set to zero." diff --git a/doc/configtables/offwind-dc.csv b/doc/configtables/offwind-dc.csv index 1f72228aa..350955973 100644 --- a/doc/configtables/offwind-dc.csv +++ b/doc/configtables/offwind-dc.csv @@ -12,5 +12,4 @@ ship_threshold,--,float,"Ship density threshold from which areas are excluded." max_depth,m,float,"Maximum sea water depth at which wind turbines can be build. Maritime areas with deeper waters are excluded in the process of calculating the AC-connected offshore wind potential." min_shore_distance,m,float,"Minimum distance to the shore below which wind turbines cannot be build." max_shore_distance,m,float,"Maximum distance to the shore above which wind turbines cannot be build." -potential,--,"One of {'simple', 'conservative'}","Method to compute the maximal installable potential for a node; confer :ref:`renewableprofiles`" clip_p_max_pu,p.u.,float,"To avoid too small values in the renewables` per-unit availability time series values below this threshold are set to zero." diff --git a/doc/configtables/onwind.csv b/doc/configtables/onwind.csv index ba9482e56..b7e823b36 100644 --- a/doc/configtables/onwind.csv +++ b/doc/configtables/onwind.csv @@ -9,7 +9,6 @@ corine,,, -- distance,m,float,"Distance to keep from areas specified in ``distance_grid_codes``" -- distance_grid_codes,--,"Any subset of the `CORINE Land Cover code list `_","Specifies areas according to CORINE Land Cover codes to which wind turbines must maintain a distance specified in the setting ``distance``." natura,bool,"{true, false}","Switch to exclude `Natura 2000 `_ natural protection areas. Area is excluded if ``true``." -potential,--,"One of {'simple', 'conservative'}","Method to compute the maximal installable potential for a node; confer :ref:`renewableprofiles`" clip_p_max_pu,p.u.,float,"To avoid too small values in the renewables` per-unit availability time series values below this threshold are set to zero." correction_factor,--,float,"Correction factor for capacity factor time series." excluder_resolution,m,float,"Resolution on which to perform geographical elibility analysis." diff --git a/doc/configtables/sector.csv b/doc/configtables/sector.csv index 856ea0749..938c120ad 100644 --- a/doc/configtables/sector.csv +++ b/doc/configtables/sector.csv @@ -71,7 +71,6 @@ solar_thermal,--,"{true, false}",Add option for using solar thermal to generate solar_cf_correction,--,float,The correction factor for the value provided by the solar thermal profile calculations marginal_cost_storage,currency/MWh ,float,The marginal cost of discharging batteries in distributed grids methanation,--,"{true, false}",Add option for transforming hydrogen and CO2 into methane using methanation. -helmeth,--,"{true, false}",Add option for transforming power into gas using HELMETH (Integrated High-Temperature ELectrolysis and METHanation for Effective Power to Gas Conversion) coal_cc,--,"{true, false}",Add option for coal CHPs with carbon capture dac,--,"{true, false}",Add option for Direct Air Capture (DAC) co2_vent,--,"{true, false}",Add option for vent out CO2 from storages to the atmosphere. @@ -108,6 +107,11 @@ electricity_distribution _grid,--,"{true, false}",Add a simplified representatio electricity_distribution _grid_cost_factor,,,Multiplies the investment cost of the electricity distribution grid ,,, electricity_grid _connection,--,"{true, false}",Add the cost of electricity grid connection for onshore wind and solar +transmission_efficiency,,,Section to specify transmission losses or compression energy demands of bidirectional links. Splits them into two capacity-linked unidirectional links. +-- {carrier},--,str,The carrier of the link. +-- -- efficiency_static,p.u.,float,Length-independent transmission efficiency. +-- -- efficiency_per_1000km,p.u. per 1000 km,float,Length-dependent transmission efficiency ($\eta^{\text{length}}$) +-- -- compression_per_1000km,p.u. per 1000 km,float,Length-dependent electricity demand for compression ($\eta \cdot \text{length}$) implemented as multi-link to local electricity bus. H2_network,--,"{true, false}",Add option for new hydrogen pipelines gas_network,--,"{true, false}","Add existing natural gas infrastructure, incl. LNG terminals, production and entry-points. The existing gas network is added with a lossless transport model. A length-weighted `k-edge augmentation algorithm `_ can be run to add new candidate gas pipelines such that all regions of the model can be connected to the gas network. When activated, all the gas demands are regionally disaggregated as well." H2_retrofit,--,"{true, false}",Add option for retrofiting existing pipelines to transport hydrogen. @@ -118,6 +122,7 @@ gas_distribution_grid _cost_factor,,,Multiplier for the investment cost of the g ,,, biomass_spatial,--,"{true, false}",Add option for resolving biomass demand regionally biomass_transport,--,"{true, false}",Add option for transporting solid biomass between nodes +biogas_upgrading_cc,--,"{true, false}",Add option to capture CO2 from biomass upgrading conventional_generation,,,Add a more detailed description of conventional carriers. Any power generation requires the consumption of fuel from nodes representing that fuel. biomass_to_liquid,--,"{true, false}",Add option for transforming solid biomass into liquid fuel with the same properties as oil biosng,--,"{true, false}",Add option for transforming solid biomass into synthesis gas with the same properties as natural gas diff --git a/doc/configtables/solar.csv b/doc/configtables/solar.csv index 803445d5d..7da1281bc 100644 --- a/doc/configtables/solar.csv +++ b/doc/configtables/solar.csv @@ -10,6 +10,5 @@ capacity_per_sqkm,:math:`MW/km^2`,float,"Allowable density of solar panel placem correction_factor,--,float,"A correction factor for the capacity factor (availability) time series." corine,--,"Any subset of the `CORINE Land Cover code list `_","Specifies areas according to CORINE Land Cover codes which are generally eligible for solar panel placement." natura,bool,"{true, false}","Switch to exclude `Natura 2000 `_ natural protection areas. Area is excluded if ``true``." -potential,--,"One of {'simple', 'conservative'}","Method to compute the maximal installable potential for a node; confer :ref:`renewableprofiles`" clip_p_max_pu,p.u.,float,"To avoid too small values in the renewables` per-unit availability time series values below this threshold are set to zero." excluder_resolution,m,float,"Resolution on which to perform geographical elibility analysis." diff --git a/doc/index.rst b/doc/index.rst index d30dd8b97..909a96a2d 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -116,7 +116,7 @@ of the individual parts. topics we are working on. Please feel free to help or make suggestions. This project is currently maintained by the `Department of Digital -Transformation in Energy Systems `_ at the +Transformation in Energy Systems `_ at the `Technische Universität Berlin `_. Previous versions were developed within the `IAI `_ at the `Karlsruhe Institute of Technology (KIT) `_ which was funded by diff --git a/doc/release_notes.rst b/doc/release_notes.rst index a7644682c..a4f20b86c 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -10,6 +10,13 @@ Release Notes Upcoming Release ================ +* Add option to specify losses for bidirectional links, e.g. pipelines or HVDC + links, in configuration file under ``sector: transmission_efficiency:``. Users + can specify static or length-dependent values as well as a length-dependent + electricity demand for compression, which is implemented as a multi-link to + the local electricity buses. The bidirectional links will then be split into + two unidirectional links with linked capacities. + * Pin ``snakemake`` version to below 8.0.0, as the new version is not yet supported by ``pypsa-eur``. @@ -38,12 +45,30 @@ Upcoming Release * Split configuration to enable SMR and SMR CC. +* Bugfix: The unit of the capital cost of Haber-Bosch plants was corrected. + * The configuration setting for country focus weights when clustering the network has been moved from ``focus_weights:`` to ``clustering: focus_weights:``. Backwards compatibility to old config files is maintained. +* Extend options for waste usage from Haber-Bosch, methanolisation and methanation. + +* Use electrolysis waste heat by default. + +* Add new ``sector_opts`` wildcard option "nowasteheat" to disable all waste heat usage. + +* Set minimum part loads for PtX processes to 30% for methanolisation and methanation, and to 70% for Fischer-Tropsch synthesis. + +* Add VOM as marginal cost to PtX processes. + +* Add pelletizing costs for biomass boilers. + * The ``mock_snakemake`` function can now be used with a Snakefile from a different directory using the new ``root_dir`` argument. +* Switch to using hydrogen and electricity inputs for Haber-Bosch from https://github.com/PyPSA/technology-data. + +* Add option to capture CO2 contained in biogas when upgrading (``sector: biogas_to_gas_cc``). + * Merged option to extend geographical scope to Ukraine and Moldova. These countries are excluded by default and is currently constrained to power-sector only parts of the workflow. A special config file @@ -59,6 +84,17 @@ Upcoming Release default setting points to an empty hull at ``data/custom_extra_functionality.py``. +* Validate downloads from Zenodo using MD5 checksums. This identifies corrupted + or incomplete downloads. + +* Add locations, capacities and costs of existing gas storage using Global + Energy Monitor's `Europe Gas Tracker + `_. + +* Remove HELMETH option. + +* Print Irreducible Infeasible Subset (IIS) if model is infeasible. Only for + solvers with IIS support. **Bugs and Compatibility** @@ -190,6 +226,8 @@ PyPSA-Eur 0.8.1 (27th July 2023) (https://github.com/PyPSA/pypsa-eur/pull/672) +* Addressed deprecation warnings for ``pandas=2.0``. ``pandas=2.0`` is now minimum requirement. + PyPSA-Eur 0.8.0 (18th March 2023) ================================= diff --git a/rules/build_sector.smk b/rules/build_sector.smk index dd49fc6f0..ab8ff4edf 100644 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -85,12 +85,12 @@ if config["sector"]["gas_network"] or config["sector"]["H2_retrofit"]: rule build_gas_input_locations: input: - lng=HTTP.remote( + gem=HTTP.remote( "https://globalenergymonitor.org/wp-content/uploads/2023/07/Europe-Gas-Tracker-2023-03-v3.xlsx", keep_local=True, ), entry="data/gas_network/scigrid-gas/data/IGGIELGN_BorderPoints.geojson", - production="data/gas_network/scigrid-gas/data/IGGIELGN_Productions.geojson", + storage="data/gas_network/scigrid-gas/data/IGGIELGN_Storages.geojson", regions_onshore=RESOURCES + "regions_onshore_elec_s{simpl}_{clusters}.geojson", regions_offshore=RESOURCES diff --git a/rules/common.smk b/rules/common.smk index d34160507..2c8cf69c1 100644 --- a/rules/common.smk +++ b/rules/common.smk @@ -2,6 +2,11 @@ # # SPDX-License-Identifier: MIT +import os, sys + +sys.path.insert(0, os.path.abspath("scripts")) +from _helpers import validate_checksum + def memory(w): factor = 3.0 diff --git a/rules/retrieve.smk b/rules/retrieve.smk index 0483a8335..4ded2a46b 100644 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -77,6 +77,7 @@ if config["enable"]["retrieve"] and config["enable"].get("retrieve_cutout", True retries: 2 run: move(input[0], output[0]) + validate_checksum(output[0], input[0]) if config["enable"]["retrieve"] and config["enable"].get("retrieve_cost_data", True): @@ -121,6 +122,7 @@ if config["enable"]["retrieve"] and config["enable"].get( retries: 2 run: move(input[0], output[0]) + validate_checksum(output[0], input[0]) if config["enable"]["retrieve"] and config["enable"].get( @@ -167,6 +169,7 @@ if config["enable"]["retrieve"] and ( "IGGIELGN_LNGs.geojson", "IGGIELGN_BorderPoints.geojson", "IGGIELGN_Productions.geojson", + "IGGIELGN_Storages.geojson", "IGGIELGN_PipeSegments.geojson", ] @@ -226,6 +229,7 @@ if config["enable"]["retrieve"]: retries: 2 run: move(input[0], output[0]) + validate_checksum(output[0], input[0]) if config["enable"]["retrieve"]: @@ -242,6 +246,7 @@ if config["enable"]["retrieve"]: "data/Copernicus_LC100_global_v3.0.1_2019-nrt_Discrete-Classification-map_EPSG-4326.tif", run: move(input[0], output[0]) + validate_checksum(output[0], input[0]) if config["enable"]["retrieve"]: diff --git a/scripts/_helpers.py b/scripts/_helpers.py index 398f3a30f..9945f70f0 100644 --- a/scripts/_helpers.py +++ b/scripts/_helpers.py @@ -4,6 +4,7 @@ # SPDX-License-Identifier: MIT import contextlib +import hashlib import logging import os import urllib @@ -11,6 +12,7 @@ import pandas as pd import pytz +import requests import yaml from pypsa.components import component_attrs, components from pypsa.descriptors import Dict @@ -318,3 +320,63 @@ def update_config_with_sector_opts(config, sector_opts): if o.startswith("CF+"): l = o.split("+")[1:] update_config(config, parse(l)) + + +def get_checksum_from_zenodo(file_url): + parts = file_url.split("/") + record_id = parts[parts.index("record") + 1] + filename = parts[-1] + + response = requests.get(f"https://zenodo.org/api/records/{record_id}", timeout=30) + response.raise_for_status() + data = response.json() + + for file in data["files"]: + if file["key"] == filename: + return file["checksum"] + return None + + +def validate_checksum(file_path, zenodo_url=None, checksum=None): + """ + Validate file checksum against provided or Zenodo-retrieved checksum. + Calculates the hash of a file using 64KB chunks. Compares it against a + given checksum or one from a Zenodo URL. + + Parameters + ---------- + file_path : str + Path to the file for checksum validation. + zenodo_url : str, optional + URL of the file on Zenodo to fetch the checksum. + checksum : str, optional + Checksum (format 'hash_type:checksum_value') for validation. + + Raises + ------ + AssertionError + If the checksum does not match, or if neither `checksum` nor `zenodo_url` is provided. + + + Examples + -------- + >>> validate_checksum("/path/to/file", checksum="md5:abc123...") + >>> validate_checksum( + ... "/path/to/file", + ... zenodo_url="https://zenodo.org/record/12345/files/example.txt", + ... ) + + If the checksum is invalid, an AssertionError will be raised. + """ + assert checksum or zenodo_url, "Either checksum or zenodo_url must be provided" + if zenodo_url: + checksum = get_checksum_from_zenodo(zenodo_url) + hash_type, checksum = checksum.split(":") + hasher = hashlib.new(hash_type) + with open(file_path, "rb") as f: + for chunk in iter(lambda: f.read(65536), b""): # 64kb chunks + hasher.update(chunk) + calculated_checksum = hasher.hexdigest() + assert ( + calculated_checksum == checksum + ), "Checksum is invalid. This may be due to an incomplete download. Delete the file and re-execute the rule." diff --git a/scripts/add_brownfield.py b/scripts/add_brownfield.py index 741025801..9ddd3d999 100644 --- a/scripts/add_brownfield.py +++ b/scripts/add_brownfield.py @@ -120,6 +120,33 @@ def add_brownfield(n, n_p, year): n.links.loc[new_pipes, "p_nom_min"] = 0.0 +def disable_grid_expansion_if_LV_limit_hit(n): + if not "lv_limit" in n.global_constraints.index: + return + + total_expansion = ( + n.lines.eval("s_nom_min * length").sum() + + n.links.query("carrier == 'DC'").eval("p_nom_min * length").sum() + ).sum() + + lv_limit = n.global_constraints.at["lv_limit", "constant"] + + # allow small numerical differences + if lv_limit - total_expansion < 1: + logger.info( + f"LV is already reached (gap {diff} MWkm), disabling expansion and LV limit" + ) + extendable_acs = n.lines.query("s_nom_extendable").index + n.lines.loc[extendable_acs, "s_nom_extendable"] = False + n.lines.loc[extendable_acs, "s_nom"] = n.lines.loc[extendable_acs, "s_nom_min"] + + extendable_dcs = n.links.query("carrier == 'DC' and p_nom_extendable").index + n.links.loc[extendable_dcs, "p_nom_extendable"] = False + n.links.loc[extendable_dcs, "p_nom"] = n.links.loc[extendable_dcs, "p_nom_min"] + + n.global_constraints.drop("lv_limit", inplace=True) + + if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake @@ -150,5 +177,7 @@ def add_brownfield(n, n_p, year): add_brownfield(n, n_p, year) + disable_grid_expansion_if_LV_limit_hit(n) + n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) n.export_to_netcdf(snakemake.output[0]) diff --git a/scripts/build_biomass_potentials.py b/scripts/build_biomass_potentials.py index d7c467cf6..aae1fb989 100644 --- a/scripts/build_biomass_potentials.py +++ b/scripts/build_biomass_potentials.py @@ -134,7 +134,7 @@ def disaggregate_nuts0(bio): # get population in nuts2 pop_nuts2 = pop.loc[pop.index.str.len() == 4] by_country = pop_nuts2.total.groupby(pop_nuts2.ct).sum() - pop_nuts2["fraction"] = pop_nuts2.total / pop_nuts2.ct.map(by_country) + pop_nuts2.loc[:, "fraction"] = pop_nuts2.total / pop_nuts2.ct.map(by_country) # distribute nuts0 data to nuts2 by population bio_nodal = bio.loc[pop_nuts2.ct] diff --git a/scripts/build_clustered_population_layouts.py b/scripts/build_clustered_population_layouts.py index 083f3de46..2f2376569 100644 --- a/scripts/build_clustered_population_layouts.py +++ b/scripts/build_clustered_population_layouts.py @@ -25,10 +25,7 @@ cutout = atlite.Cutout(snakemake.input.cutout) clustered_regions = ( - gpd.read_file(snakemake.input.regions_onshore) - .set_index("name") - .buffer(0) - .squeeze() + gpd.read_file(snakemake.input.regions_onshore).set_index("name").buffer(0) ) I = cutout.indicatormatrix(clustered_regions) diff --git a/scripts/build_electricity_demand.py b/scripts/build_electricity_demand.py index f7b6cddd0..d7d9927d6 100755 --- a/scripts/build_electricity_demand.py +++ b/scripts/build_electricity_demand.py @@ -81,7 +81,7 @@ def rename(s): return s[: -len(pattern)] return ( - pd.read_csv(fn, index_col=0, parse_dates=[0]) + pd.read_csv(fn, index_col=0, parse_dates=[0], date_format="%Y-%m-%dT%H:%M:%SZ") .tz_localize(None) .filter(like=pattern) .rename(columns=rename) diff --git a/scripts/build_energy_totals.py b/scripts/build_energy_totals.py index 6f9585c1d..67b864663 100644 --- a/scripts/build_energy_totals.py +++ b/scripts/build_energy_totals.py @@ -189,12 +189,12 @@ def idees_per_country(ct, year, base_dir): ct_totals["total residential water"] = df.at["Water heating"] assert df.index[23] == "Electricity" - ct_totals["electricity residential water"] = df[23] + ct_totals["electricity residential water"] = df.iloc[23] ct_totals["total residential cooking"] = df["Cooking"] assert df.index[30] == "Electricity" - ct_totals["electricity residential cooking"] = df[30] + ct_totals["electricity residential cooking"] = df.iloc[30] df = pd.read_excel(fn_residential, "RES_summary", index_col=0)[year] @@ -202,13 +202,13 @@ def idees_per_country(ct, year, base_dir): ct_totals["total residential"] = df[row] assert df.index[47] == "Electricity" - ct_totals["electricity residential"] = df[47] + ct_totals["electricity residential"] = df.iloc[47] assert df.index[46] == "Derived heat" - ct_totals["derived heat residential"] = df[46] + ct_totals["derived heat residential"] = df.iloc[46] assert df.index[50] == "Thermal uses" - ct_totals["thermal uses residential"] = df[50] + ct_totals["thermal uses residential"] = df.iloc[50] # services @@ -222,12 +222,12 @@ def idees_per_country(ct, year, base_dir): ct_totals["total services water"] = df["Hot water"] assert df.index[24] == "Electricity" - ct_totals["electricity services water"] = df[24] + ct_totals["electricity services water"] = df.iloc[24] ct_totals["total services cooking"] = df["Catering"] assert df.index[31] == "Electricity" - ct_totals["electricity services cooking"] = df[31] + ct_totals["electricity services cooking"] = df.iloc[31] df = pd.read_excel(fn_tertiary, "SER_summary", index_col=0)[year] @@ -235,13 +235,13 @@ def idees_per_country(ct, year, base_dir): ct_totals["total services"] = df[row] assert df.index[50] == "Electricity" - ct_totals["electricity services"] = df[50] + ct_totals["electricity services"] = df.iloc[50] assert df.index[49] == "Derived heat" - ct_totals["derived heat services"] = df[49] + ct_totals["derived heat services"] = df.iloc[49] assert df.index[53] == "Thermal uses" - ct_totals["thermal uses services"] = df[53] + ct_totals["thermal uses services"] = df.iloc[53] # agriculture, forestry and fishing @@ -282,28 +282,28 @@ def idees_per_country(ct, year, base_dir): ct_totals["total two-wheel"] = df["Powered 2-wheelers (Gasoline)"] assert df.index[19] == "Passenger cars" - ct_totals["total passenger cars"] = df[19] + ct_totals["total passenger cars"] = df.iloc[19] assert df.index[30] == "Battery electric vehicles" - ct_totals["electricity passenger cars"] = df[30] + ct_totals["electricity passenger cars"] = df.iloc[30] assert df.index[31] == "Motor coaches, buses and trolley buses" - ct_totals["total other road passenger"] = df[31] + ct_totals["total other road passenger"] = df.iloc[31] assert df.index[39] == "Battery electric vehicles" - ct_totals["electricity other road passenger"] = df[39] + ct_totals["electricity other road passenger"] = df.iloc[39] assert df.index[41] == "Light duty vehicles" - ct_totals["total light duty road freight"] = df[41] + ct_totals["total light duty road freight"] = df.iloc[41] assert df.index[49] == "Battery electric vehicles" - ct_totals["electricity light duty road freight"] = df[49] + ct_totals["electricity light duty road freight"] = df.iloc[49] row = "Heavy duty vehicles (Diesel oil incl. biofuels)" ct_totals["total heavy duty road freight"] = df[row] assert df.index[61] == "Passenger cars" - ct_totals["passenger car efficiency"] = df[61] + ct_totals["passenger car efficiency"] = df.iloc[61] df = pd.read_excel(fn_transport, "TrRail_ene", index_col=0)[year] @@ -312,39 +312,39 @@ def idees_per_country(ct, year, base_dir): ct_totals["electricity rail"] = df["Electricity"] assert df.index[15] == "Passenger transport" - ct_totals["total rail passenger"] = df[15] + ct_totals["total rail passenger"] = df.iloc[15] assert df.index[16] == "Metro and tram, urban light rail" assert df.index[19] == "Electric" assert df.index[20] == "High speed passenger trains" - ct_totals["electricity rail passenger"] = df[[16, 19, 20]].sum() + ct_totals["electricity rail passenger"] = df.iloc[[16, 19, 20]].sum() assert df.index[21] == "Freight transport" - ct_totals["total rail freight"] = df[21] + ct_totals["total rail freight"] = df.iloc[21] assert df.index[23] == "Electric" - ct_totals["electricity rail freight"] = df[23] + ct_totals["electricity rail freight"] = df.iloc[23] df = pd.read_excel(fn_transport, "TrAvia_ene", index_col=0)[year] assert df.index[6] == "Passenger transport" - ct_totals["total aviation passenger"] = df[6] + ct_totals["total aviation passenger"] = df.iloc[6] assert df.index[10] == "Freight transport" - ct_totals["total aviation freight"] = df[10] + ct_totals["total aviation freight"] = df.iloc[10] assert df.index[7] == "Domestic" - ct_totals["total domestic aviation passenger"] = df[7] + ct_totals["total domestic aviation passenger"] = df.iloc[7] assert df.index[8] == "International - Intra-EU" assert df.index[9] == "International - Extra-EU" - ct_totals["total international aviation passenger"] = df[[8, 9]].sum() + ct_totals["total international aviation passenger"] = df.iloc[[8, 9]].sum() assert df.index[11] == "Domestic and International - Intra-EU" - ct_totals["total domestic aviation freight"] = df[11] + ct_totals["total domestic aviation freight"] = df.iloc[11] assert df.index[12] == "International - Extra-EU" - ct_totals["total international aviation freight"] = df[12] + ct_totals["total international aviation freight"] = df.iloc[12] ct_totals["total domestic aviation"] = ( ct_totals["total domestic aviation freight"] @@ -364,7 +364,7 @@ def idees_per_country(ct, year, base_dir): df = pd.read_excel(fn_transport, "TrRoad_act", index_col=0)[year] assert df.index[85] == "Passenger cars" - ct_totals["passenger cars"] = df[85] + ct_totals["passenger cars"] = df.iloc[85] return pd.Series(ct_totals, name=ct) diff --git a/scripts/build_gas_input_locations.py b/scripts/build_gas_input_locations.py index a3b945abc..9ad3760d9 100644 --- a/scripts/build_gas_input_locations.py +++ b/scripts/build_gas_input_locations.py @@ -23,11 +23,10 @@ def read_scigrid_gas(fn): return df -def build_gem_lng_data(lng_fn): - df = pd.read_excel(lng_fn[0], sheet_name="LNG terminals - data") +def build_gem_lng_data(fn): + df = pd.read_excel(fn[0], sheet_name="LNG terminals - data") df = df.set_index("ComboID") - remove_status = ["Cancelled"] remove_country = ["Cyprus", "Turkey"] remove_terminal = ["Puerto de la Luz LNG Terminal", "Gran Canaria LNG Terminal"] @@ -42,9 +41,50 @@ def build_gem_lng_data(lng_fn): return gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326") -def build_gas_input_locations(lng_fn, entry_fn, prod_fn, countries): +def build_gem_prod_data(fn): + df = pd.read_excel(fn[0], sheet_name="Gas extraction - main") + df = df.set_index("GEM Unit ID") + + remove_country = ["Cyprus", "Türkiye"] + remove_fuel_type = ["oil"] + + df = df.query( + "Status != 'shut in' \ + & 'Fuel type' != 'oil' \ + & Country != @remove_country \ + & ~Latitude.isna() \ + & ~Longitude.isna()" + ).copy() + + p = pd.read_excel(fn[0], sheet_name="Gas extraction - production") + p = p.set_index("GEM Unit ID") + p = p[p["Fuel description"] == "gas"] + + capacities = pd.DataFrame(index=df.index) + for key in ["production", "production design capacity", "reserves"]: + cap = ( + p.loc[p["Production/reserves"] == key, "Quantity (converted)"] + .groupby("GEM Unit ID") + .sum() + .reindex(df.index) + ) + # assume capacity such that 3% of reserves can be extracted per year (25% quantile) + annualization_factor = 0.03 if key == "reserves" else 1.0 + capacities[key] = cap * annualization_factor + + df["mcm_per_year"] = ( + capacities["production"] + .combine_first(capacities["production design capacity"]) + .combine_first(capacities["reserves"]) + ) + + geometry = gpd.points_from_xy(df["Longitude"], df["Latitude"]) + return gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326") + + +def build_gas_input_locations(gem_fn, entry_fn, sto_fn, countries): # LNG terminals - lng = build_gem_lng_data(lng_fn) + lng = build_gem_lng_data(gem_fn) # Entry points from outside the model scope entry = read_scigrid_gas(entry_fn) @@ -55,25 +95,30 @@ def build_gas_input_locations(lng_fn, entry_fn, prod_fn, countries): | (entry.from_country == "NO") # malformed datapoint # entries from NO to GB ] + sto = read_scigrid_gas(sto_fn) + remove_country = ["RU", "UA", "TR", "BY"] + sto = sto.query("country_code != @remove_country") + # production sites inside the model scope - prod = read_scigrid_gas(prod_fn) - prod = prod.loc[ - (prod.geometry.y > 35) & (prod.geometry.x < 30) & (prod.country_code != "DE") - ] + prod = build_gem_prod_data(gem_fn) mcm_per_day_to_mw = 437.5 # MCM/day to MWh/h + mcm_per_year_to_mw = 1.199 # MCM/year to MWh/h mtpa_to_mw = 1649.224 # mtpa to MWh/h - lng["p_nom"] = lng["CapacityInMtpa"] * mtpa_to_mw - entry["p_nom"] = entry["max_cap_from_to_M_m3_per_d"] * mcm_per_day_to_mw - prod["p_nom"] = prod["max_supply_M_m3_per_d"] * mcm_per_day_to_mw + mcm_to_gwh = 11.36 # MCM to GWh + lng["capacity"] = lng["CapacityInMtpa"] * mtpa_to_mw + entry["capacity"] = entry["max_cap_from_to_M_m3_per_d"] * mcm_per_day_to_mw + prod["capacity"] = prod["mcm_per_year"] * mcm_per_year_to_mw + sto["capacity"] = sto["max_cushionGas_M_m3"] * mcm_to_gwh lng["type"] = "lng" entry["type"] = "pipeline" prod["type"] = "production" + sto["type"] = "storage" - sel = ["geometry", "p_nom", "type"] + sel = ["geometry", "capacity", "type"] - return pd.concat([prod[sel], entry[sel], lng[sel]], ignore_index=True) + return pd.concat([prod[sel], entry[sel], lng[sel], sto[sel]], ignore_index=True) if __name__ == "__main__": @@ -83,7 +128,7 @@ def build_gas_input_locations(lng_fn, entry_fn, prod_fn, countries): snakemake = mock_snakemake( "build_gas_input_locations", simpl="", - clusters="37", + clusters="128", ) logging.basicConfig(level=snakemake.config["logging"]["level"]) @@ -104,9 +149,9 @@ def build_gas_input_locations(lng_fn, entry_fn, prod_fn, countries): countries = regions.index.str[:2].unique().str.replace("GB", "UK") gas_input_locations = build_gas_input_locations( - snakemake.input.lng, + snakemake.input.gem, snakemake.input.entry, - snakemake.input.production, + snakemake.input.storage, countries, ) @@ -116,9 +161,13 @@ def build_gas_input_locations(lng_fn, entry_fn, prod_fn, countries): gas_input_nodes.to_file(snakemake.output.gas_input_nodes, driver="GeoJSON") + ensure_columns = ["lng", "pipeline", "production", "storage"] gas_input_nodes_s = ( - gas_input_nodes.groupby(["bus", "type"])["p_nom"].sum().unstack() + gas_input_nodes.groupby(["bus", "type"])["capacity"] + .sum() + .unstack() + .reindex(columns=ensure_columns) ) - gas_input_nodes_s.columns.name = "p_nom" + gas_input_nodes_s.columns.name = "capacity" gas_input_nodes_s.to_csv(snakemake.output.gas_input_nodes_simplified) diff --git a/scripts/build_heat_demand.py b/scripts/build_heat_demand.py index 734942600..777684043 100644 --- a/scripts/build_heat_demand.py +++ b/scripts/build_heat_demand.py @@ -31,10 +31,7 @@ cutout = atlite.Cutout(snakemake.input.cutout).sel(time=time) clustered_regions = ( - gpd.read_file(snakemake.input.regions_onshore) - .set_index("name") - .buffer(0) - .squeeze() + gpd.read_file(snakemake.input.regions_onshore).set_index("name").buffer(0) ) I = cutout.indicatormatrix(clustered_regions) diff --git a/scripts/build_line_rating.py b/scripts/build_line_rating.py index 032ba39cf..c53d2899b 100755 --- a/scripts/build_line_rating.py +++ b/scripts/build_line_rating.py @@ -119,7 +119,7 @@ def calculate_line_rating(n, cutout): .apply(lambda x: int(re.findall(r"(\d+)-bundle", x)[0])) ) # Set default number of bundles per line - relevant_lines["n_bundle"].fillna(1, inplace=True) + relevant_lines["n_bundle"] = relevant_lines["n_bundle"].fillna(1) R *= relevant_lines["n_bundle"] R = calculate_resistance(T=353, R_ref=R) Imax = cutout.line_rating(shapes, R, D=0.0218, Ts=353, epsilon=0.8, alpha=0.8) diff --git a/scripts/build_renewable_profiles.py b/scripts/build_renewable_profiles.py index 3a1c525e0..c579f588d 100644 --- a/scripts/build_renewable_profiles.py +++ b/scripts/build_renewable_profiles.py @@ -26,20 +26,9 @@ renewable: {technology}: - cutout: - corine: - grid_codes: - distance: - natura: - max_depth: - max_shore_distance: - min_shore_distance: - capacity_per_sqkm: - correction_factor: - potential: - min_p_max_pu: - clip_p_max_pu: - resource: + cutout: corine: grid_codes: distance: natura: max_depth: + max_shore_distance: min_shore_distance: capacity_per_sqkm: + correction_factor: min_p_max_pu: clip_p_max_pu: resource: .. seealso:: Documentation of the configuration file ``config/config.yaml`` at @@ -48,21 +37,30 @@ Inputs ------ -- ``data/bundle/corine/g250_clc06_V18_5.tif``: `CORINE Land Cover (CLC) `_ inventory on `44 classes `_ of land use (e.g. forests, arable land, industrial, urban areas). +- ``data/bundle/corine/g250_clc06_V18_5.tif``: `CORINE Land Cover (CLC) + `_ inventory on `44 + classes `_ of + land use (e.g. forests, arable land, industrial, urban areas). .. image:: img/corine.png :scale: 33 % -- ``data/bundle/GEBCO_2014_2D.nc``: A `bathymetric `_ data set with a global terrain model for ocean and land at 15 arc-second intervals by the `General Bathymetric Chart of the Oceans (GEBCO) `_. +- ``data/bundle/GEBCO_2014_2D.nc``: A `bathymetric + `_ data set with a global terrain + model for ocean and land at 15 arc-second intervals by the `General + Bathymetric Chart of the Oceans (GEBCO) + `_. .. image:: img/gebco_2019_grid_image.jpg :scale: 50 % - **Source:** `GEBCO `_ + **Source:** `GEBCO + `_ - ``resources/natura.tiff``: confer :ref:`natura` - ``resources/offshore_shapes.geojson``: confer :ref:`shapes` -- ``resources/regions_onshore.geojson``: (if not offshore wind), confer :ref:`busregions` +- ``resources/regions_onshore.geojson``: (if not offshore wind), confer + :ref:`busregions` - ``resources/regions_offshore.geojson``: (if offshore wind), :ref:`busregions` - ``"cutouts/" + params["renewable"][{technology}]['cutout']``: :ref:`cutout` - ``networks/base.nc``: :ref:`base` @@ -128,25 +126,25 @@ This script functions at two main spatial resolutions: the resolution of the network nodes and their `Voronoi cells `_, and the resolution of the -cutout grid cells for the weather data. Typically the weather data grid is -finer than the network nodes, so we have to work out the distribution of -generators across the grid cells within each Voronoi cell. This is done by -taking account of a combination of the available land at each grid cell and the -capacity factor there. +cutout grid cells for the weather data. Typically the weather data grid is finer +than the network nodes, so we have to work out the distribution of generators +across the grid cells within each Voronoi cell. This is done by taking account +of a combination of the available land at each grid cell and the capacity factor +there. First the script computes how much of the technology can be installed at each cutout grid cell and each node using the `GLAES -`_ library. This uses the CORINE land use data, -Natura2000 nature reserves and GEBCO bathymetry data. +`_ library. This uses the CORINE land use +data, Natura2000 nature reserves and GEBCO bathymetry data. .. image:: img/eligibility.png :scale: 50 % :align: center -To compute the layout of generators in each node's Voronoi cell, the -installable potential in each grid cell is multiplied with the capacity factor -at each grid cell. This is done since we assume more generators are installed -at cells with a higher capacity factor. +To compute the layout of generators in each node's Voronoi cell, the installable +potential in each grid cell is multiplied with the capacity factor at each grid +cell. This is done since we assume more generators are installed at cells with a +higher capacity factor. .. image:: img/offwinddc-gridcell.png :scale: 50 % @@ -164,20 +162,14 @@ :scale: 50 % :align: center -This layout is then used to compute the generation availability time series -from the weather data cutout from ``atlite``. +This layout is then used to compute the generation availability time series from +the weather data cutout from ``atlite``. -Two methods are available to compute the maximal installable potential for the -node (`p_nom_max`): ``simple`` and ``conservative``: - -- ``simple`` adds up the installable potentials of the individual grid cells. - If the model comes close to this limit, then the time series may slightly - overestimate production since it is assumed the geographical distribution is - proportional to capacity factor. - -- ``conservative`` assertains the nodal limit by increasing capacities - proportional to the layout until the limit of an individual grid cell is - reached. +The maximal installable potential for the node (`p_nom_max`) is computed by +adding up the installable potentials of the individual grid cells. +If the model comes close to this limit, then the time series may slightly +overestimate production since it is assumed the geographical distribution is +proportional to capacity factor. """ import functools import logging @@ -210,7 +202,6 @@ resource = params["resource"] # pv panel params / wind turbine params correction_factor = params.get("correction_factor", 1.0) capacity_per_sqkm = params["capacity_per_sqkm"] - p_nom_max_meth = params.get("potential", "conservative") if isinstance(params.get("corine", {}), list): params["corine"] = {"grid_codes": params["corine"]} @@ -277,15 +268,14 @@ snakemake.input.country_shapes, buffer=buffer, invert=True ) + logger.info("Calculate landuse availability...") + start = time.time() + kwargs = dict(nprocesses=nprocesses, disable_progressbar=noprogress) - if noprogress: - logger.info("Calculate landuse availabilities...") - start = time.time() - availability = cutout.availabilitymatrix(regions, excluder, **kwargs) - duration = time.time() - start - logger.info(f"Completed availability calculation ({duration:2.2f}s)") - else: - availability = cutout.availabilitymatrix(regions, excluder, **kwargs) + availability = cutout.availabilitymatrix(regions, excluder, **kwargs) + + duration = time.time() - start + logger.info(f"Completed landuse availability calculation ({duration:2.2f}s)") # For Moldova and Ukraine: Overwrite parts not covered by Corine with # externally determined available areas @@ -304,8 +294,19 @@ func = getattr(cutout, resource.pop("method")) if client is not None: resource["dask_kwargs"] = {"scheduler": client} + + logger.info("Calculate average capacity factor...") + start = time.time() + capacity_factor = correction_factor * func(capacity_factor=True, **resource) layout = capacity_factor * area * capacity_per_sqkm + + duration = time.time() - start + logger.info(f"Completed average capacity factor calculation ({duration:2.2f}s)") + + logger.info("Calculate weighted capacity factor time series...") + start = time.time() + profile, capacities = func( matrix=availability.stack(spatial=["y", "x"]), layout=layout, @@ -315,17 +316,13 @@ **resource, ) - logger.info(f"Calculating maximal capacity per bus (method '{p_nom_max_meth}')") - if p_nom_max_meth == "simple": - p_nom_max = capacity_per_sqkm * availability @ area - elif p_nom_max_meth == "conservative": - max_cap_factor = capacity_factor.where(availability != 0).max(["x", "y"]) - p_nom_max = capacities / max_cap_factor - else: - raise AssertionError( - 'Config key `potential` should be one of "simple" ' - f'(default) or "conservative", not "{p_nom_max_meth}"' - ) + duration = time.time() - start + logger.info( + f"Completed weighted capacity factor time series calculation ({duration:2.2f}s)" + ) + + logger.info(f"Calculating maximal capacity per bus") + p_nom_max = capacity_per_sqkm * availability @ area logger.info("Calculate average distances.") layoutmatrix = (layout * availability).stack(spatial=["y", "x"]) diff --git a/scripts/build_retro_cost.py b/scripts/build_retro_cost.py index f5313c21f..3ca2b1741 100644 --- a/scripts/build_retro_cost.py +++ b/scripts/build_retro_cost.py @@ -836,9 +836,9 @@ def calculate_heat_losses(u_values, data_tabula, l_strength, temperature_factor) F_red_temp = map_to_lstrength(l_strength, F_red_temp) Q_ht = ( - heat_transfer_perm2.groupby(level=1, axis=1) + heat_transfer_perm2.T.groupby(level=1) .sum() - .mul(F_red_temp.droplevel(0, axis=1)) + .T.mul(F_red_temp.droplevel(0, axis=1)) .mul(temperature_factor.reindex(heat_transfer_perm2.index, level=0), axis=0) ) @@ -878,7 +878,7 @@ def calculate_gain_utilisation_factor(heat_transfer_perm2, Q_ht, Q_gain): Calculates gain utilisation factor nu. """ # time constant of the building tau [h] = c_m [Wh/(m^2K)] * 1 /(H_tr_e+H_tb*H_ve) [m^2 K /W] - tau = c_m / heat_transfer_perm2.groupby(level=1, axis=1).sum() + tau = c_m / heat_transfer_perm2.T.groupby(axis=1).sum().T alpha = alpha_H_0 + (tau / tau_H_0) # heat balance ratio gamma = (1 / Q_ht).mul(Q_gain.sum(axis=1), axis=0) diff --git a/scripts/build_ship_raster.py b/scripts/build_ship_raster.py index 90e006b0b..02f4d5d5d 100644 --- a/scripts/build_ship_raster.py +++ b/scripts/build_ship_raster.py @@ -64,7 +64,7 @@ with zipfile.ZipFile(snakemake.input.ship_density) as zip_f: zip_f.extract("shipdensity_global.tif") with rioxarray.open_rasterio("shipdensity_global.tif") as ship_density: - ship_density = ship_density.drop(["band"]).sel( + ship_density = ship_density.drop_vars(["band"]).sel( x=slice(min(xs), max(Xs)), y=slice(max(Ys), min(ys)) ) ship_density.rio.to_raster(snakemake.output[0]) diff --git a/scripts/build_solar_thermal_profiles.py b/scripts/build_solar_thermal_profiles.py index d285691a4..ee6ed881d 100644 --- a/scripts/build_solar_thermal_profiles.py +++ b/scripts/build_solar_thermal_profiles.py @@ -33,10 +33,7 @@ cutout = atlite.Cutout(snakemake.input.cutout).sel(time=time) clustered_regions = ( - gpd.read_file(snakemake.input.regions_onshore) - .set_index("name") - .buffer(0) - .squeeze() + gpd.read_file(snakemake.input.regions_onshore).set_index("name").buffer(0) ) I = cutout.indicatormatrix(clustered_regions) diff --git a/scripts/build_temperature_profiles.py b/scripts/build_temperature_profiles.py index 9db37c257..a13ec3c21 100644 --- a/scripts/build_temperature_profiles.py +++ b/scripts/build_temperature_profiles.py @@ -31,10 +31,7 @@ cutout = atlite.Cutout(snakemake.input.cutout).sel(time=time) clustered_regions = ( - gpd.read_file(snakemake.input.regions_onshore) - .set_index("name") - .buffer(0) - .squeeze() + gpd.read_file(snakemake.input.regions_onshore).set_index("name").buffer(0) ) I = cutout.indicatormatrix(clustered_regions) diff --git a/scripts/plot_network.py b/scripts/plot_network.py index f44bb6de3..674811208 100644 --- a/scripts/plot_network.py +++ b/scripts/plot_network.py @@ -31,7 +31,7 @@ def rename_techs_tyndp(tech): tech = rename_techs(tech) if "heat pump" in tech or "resistive heater" in tech: return "power-to-heat" - elif tech in ["H2 Electrolysis", "methanation", "helmeth", "H2 liquefaction"]: + elif tech in ["H2 Electrolysis", "methanation", "H2 liquefaction"]: return "power-to-gas" elif tech == "H2": return "H2 storage" @@ -495,7 +495,7 @@ def plot_ch4_map(network): # make a fake MultiIndex so that area is correct for legend fossil_gas.index = pd.MultiIndex.from_product([fossil_gas.index, ["fossil gas"]]) - methanation_i = n.links[n.links.carrier.isin(["helmeth", "Sabatier"])].index + methanation_i = n.links.query("carrier == 'Sabatier'").index methanation = ( abs( n.links_t.p1.loc[:, methanation_i].mul( diff --git a/scripts/plot_summary.py b/scripts/plot_summary.py index 5804e7854..67ac9b553 100644 --- a/scripts/plot_summary.py +++ b/scripts/plot_summary.py @@ -121,7 +121,6 @@ def rename_techs(label): "gas boiler", "gas", "natural gas", - "helmeth", "methanation", "ammonia", "hydrogen storage", diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index d5ca27a7b..b77733062 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -95,12 +95,14 @@ def define_spatial(nodes, options): spatial.gas.industry = nodes + " gas for industry" spatial.gas.industry_cc = nodes + " gas for industry CC" spatial.gas.biogas_to_gas = nodes + " biogas to gas" + spatial.gas.biogas_to_gas_cc = nodes + "biogas to gas CC" else: spatial.gas.nodes = ["EU gas"] spatial.gas.locations = ["EU"] spatial.gas.biogas = ["EU biogas"] spatial.gas.industry = ["gas for industry"] spatial.gas.biogas_to_gas = ["EU biogas to gas"] + spatial.gas.biogas_to_gas_cc = ["EU biogas to gas CC"] if options.get("co2_spatial", options["co2network"]): spatial.gas.industry_cc = nodes + " gas for industry CC" else: @@ -452,10 +454,11 @@ def add_carrier_buses(n, carrier, nodes=None): n.add("Carrier", carrier) unit = "MWh_LHV" if carrier == "gas" else "MWh_th" + # preliminary value for non-gas carriers to avoid zeros + capital_cost = costs.at["gas storage", "fixed"] if carrier == "gas" else 0.02 n.madd("Bus", nodes, location=location, carrier=carrier, unit=unit) - # capital cost could be corrected to e.g. 0.2 EUR/kWh * annuity and O&M n.madd( "Store", nodes + " Store", @@ -463,8 +466,7 @@ def add_carrier_buses(n, carrier, nodes=None): e_nom_extendable=True, e_cyclic=True, carrier=carrier, - capital_cost=0.2 - * costs.at[carrier, "discount rate"], # preliminary value to avoid zeros + capital_cost=capital_cost, ) n.madd( @@ -805,14 +807,13 @@ def add_ammonia(n, costs): bus2=nodes + " H2", p_nom_extendable=True, carrier="Haber-Bosch", - efficiency=1 - / ( - cf_industry["MWh_elec_per_tNH3_electrolysis"] - / cf_industry["MWh_NH3_per_tNH3"] - ), # output: MW_NH3 per MW_elec - efficiency2=-cf_industry["MWh_H2_per_tNH3_electrolysis"] - / cf_industry["MWh_elec_per_tNH3_electrolysis"], # input: MW_H2 per MW_elec - capital_cost=costs.at["Haber-Bosch", "fixed"], + efficiency=1 / costs.at["Haber-Bosch", "electricity-input"], + efficiency2=-costs.at["Haber-Bosch", "hydrogen-input"] + / costs.at["Haber-Bosch", "electricity-input"], + capital_cost=costs.at["Haber-Bosch", "fixed"] + / costs.at["Haber-Bosch", "electricity-input"], + marginal_cost=costs.at["Haber-Bosch", "VOM"] + / costs.at["Haber-Bosch", "electricity-input"], lifetime=costs.at["Haber-Bosch", "lifetime"], ) @@ -1023,7 +1024,7 @@ def insert_gas_distribution_costs(n, costs): f"Inserting gas distribution grid with investment cost factor of {f_costs}" ) - capital_cost = costs.loc["electricity distribution grid"]["fixed"] * f_costs + capital_cost = costs.at["electricity distribution grid", "fixed"] * f_costs # gas boilers gas_b = n.links.index[ @@ -1100,6 +1101,7 @@ def add_storage_and_grids(n, costs): efficiency=costs.at["OCGT", "efficiency"], capital_cost=costs.at["OCGT", "fixed"] * costs.at["OCGT", "efficiency"], # NB: fixed cost is per MWel + marginal_cost=costs.at["OCGT", "VOM"], lifetime=costs.at["OCGT", "lifetime"], ) @@ -1160,7 +1162,7 @@ def add_storage_and_grids(n, costs): if options["gas_network"]: logger.info( - "Add natural gas infrastructure, incl. LNG terminals, production and entry-points." + "Add natural gas infrastructure, incl. LNG terminals, production, storage and entry-points." ) if options["H2_retrofit"]: @@ -1205,10 +1207,25 @@ def add_storage_and_grids(n, costs): remove_i = n.generators[gas_i & internal_i].index n.generators.drop(remove_i, inplace=True) - p_nom = gas_input_nodes.sum(axis=1).rename(lambda x: x + " gas") + input_types = ["lng", "pipeline", "production"] + p_nom = gas_input_nodes[input_types].sum(axis=1).rename(lambda x: x + " gas") n.generators.loc[gas_i, "p_nom_extendable"] = False n.generators.loc[gas_i, "p_nom"] = p_nom + # add existing gas storage capacity + gas_i = n.stores.carrier == "gas" + e_nom = ( + gas_input_nodes["storage"] + .rename(lambda x: x + " gas Store") + .reindex(n.stores.index) + .fillna(0.0) + * 1e3 + ) # MWh_LHV + e_nom.clip( + upper=e_nom.quantile(0.98), inplace=True + ) # limit extremely large storage + n.stores.loc[gas_i, "e_nom_min"] = e_nom + # add candidates for new gas pipelines to achieve full connectivity G = nx.Graph() @@ -1343,6 +1360,7 @@ def add_storage_and_grids(n, costs): bus2=spatial.co2.nodes, p_nom_extendable=True, carrier="Sabatier", + p_min_pu=options.get("min_part_load_methanation", 0), efficiency=costs.at["methanation", "efficiency"], efficiency2=-costs.at["methanation", "efficiency"] * costs.at["gas", "CO2 intensity"], @@ -1351,23 +1369,6 @@ def add_storage_and_grids(n, costs): lifetime=costs.at["methanation", "lifetime"], ) - if options["helmeth"]: - n.madd( - "Link", - spatial.nodes, - suffix=" helmeth", - bus0=nodes, - bus1=spatial.gas.nodes, - bus2=spatial.co2.nodes, - carrier="helmeth", - p_nom_extendable=True, - efficiency=costs.at["helmeth", "efficiency"], - efficiency2=-costs.at["helmeth", "efficiency"] - * costs.at["gas", "CO2 intensity"], - capital_cost=costs.at["helmeth", "fixed"], - lifetime=costs.at["helmeth", "lifetime"], - ) - if options.get("coal_cc"): n.madd( "Link", @@ -1630,7 +1631,7 @@ def build_heat_demand(n): electric_nodes = n.loads.index[n.loads.carrier == "electricity"] n.loads_t.p_set[electric_nodes] = ( n.loads_t.p_set[electric_nodes] - - electric_heat_supply.groupby(level=1, axis=1).sum()[electric_nodes] + - electric_heat_supply.T.groupby(level=1).sum().T[electric_nodes] ) return heat_demand @@ -1697,6 +1698,7 @@ def add_heat(n, costs): n.madd( "Generator", nodes[name] + f" {name} heat vent", + bus=nodes[name] + f" {name} heat", location=nodes[name], carrier=name + " heat vent", p_nom_extendable=True, @@ -1723,15 +1725,17 @@ def add_heat(n, costs): if sector in name: heat_load = ( heat_demand[[sector + " water", sector + " space"]] - .groupby(level=1, axis=1) - .sum()[nodes[name]] + .T.groupby(level=1) + .sum() + .T[nodes[name]] .multiply(factor) ) if name == "urban central": heat_load = ( - heat_demand.groupby(level=1, axis=1) - .sum()[nodes[name]] + heat_demand.T.groupby(level=1) + .sum() + .T[nodes[name]] .multiply( factor * (1 + options["district_heating"]["district_heating_loss"]) ) @@ -2165,12 +2169,42 @@ def add_biomass(n, costs): bus1=spatial.gas.nodes, bus2="co2 atmosphere", carrier="biogas to gas", - capital_cost=costs.loc["biogas upgrading", "fixed"], - marginal_cost=costs.loc["biogas upgrading", "VOM"], + capital_cost=costs.at["biogas", "fixed"] + + costs.at["biogas upgrading", "fixed"], + marginal_cost=costs.at["biogas upgrading", "VOM"], + efficiency=costs.at["biogas", "efficiency"], efficiency2=-costs.at["gas", "CO2 intensity"], p_nom_extendable=True, ) + if options.get("biogas_upgrading_cc"): + # Assuming for costs that the CO2 from upgrading is pure, such as in amine scrubbing. I.e., with and without CC is + # equivalent. Adding biomass CHP capture because biogas is often small-scale and decentral so further + # from e.g. CO2 grid or buyers. This is a proxy for the added cost for e.g. a raw biogas pipeline to a central upgrading facility + + n.madd( + "Link", + spatial.gas.biogas_to_gas_cc, + bus0=spatial.gas.biogas, + bus1=spatial.gas.nodes, + bus2="co2 stored", + bus3="co2 atmosphere", + carrier="biogas to gas CC", + capital_cost=costs.at["biogas CC", "fixed"] + + costs.at["biogas upgrading", "fixed"] + + costs.at["biomass CHP capture", "fixed"] + * costs.at["biogas CC", "CO2 stored"], + marginal_cost=costs.at["biogas CC", "VOM"] + + costs.at["biogas upgrading", "VOM"], + efficiency=costs.at["biogas CC", "efficiency"], + efficiency2=costs.at["biogas CC", "CO2 stored"] + * costs.at["biogas CC", "capture rate"], + efficiency3=-costs.at["gas", "CO2 intensity"] + - costs.at["biogas CC", "CO2 stored"] + * costs.at["biogas CC", "capture rate"], + p_nom_extendable=True, + ) + if options["biomass_transport"]: # add biomass transport transport_costs = pd.read_csv( @@ -2296,6 +2330,7 @@ def add_biomass(n, costs): efficiency=costs.at["biomass boiler", "efficiency"], capital_cost=costs.at["biomass boiler", "efficiency"] * costs.at["biomass boiler", "fixed"], + marginal_cost=costs.at["biomass boiler", "pelletizing cost"], lifetime=costs.at["biomass boiler", "lifetime"], ) @@ -2315,7 +2350,7 @@ def add_biomass(n, costs): + costs.at["BtL", "CO2 stored"], p_nom_extendable=True, capital_cost=costs.at["BtL", "fixed"], - marginal_cost=costs.at["BtL", "efficiency"] * costs.loc["BtL", "VOM"], + marginal_cost=costs.at["BtL", "efficiency"] * costs.at["BtL", "VOM"], ) # TODO: Update with energy penalty @@ -2336,7 +2371,7 @@ def add_biomass(n, costs): p_nom_extendable=True, capital_cost=costs.at["BtL", "fixed"] + costs.at["biomass CHP capture", "fixed"] * costs.at["BtL", "CO2 stored"], - marginal_cost=costs.at["BtL", "efficiency"] * costs.loc["BtL", "VOM"], + marginal_cost=costs.at["BtL", "efficiency"] * costs.at["BtL", "VOM"], ) # BioSNG from solid biomass @@ -2355,7 +2390,7 @@ def add_biomass(n, costs): + costs.at["BioSNG", "CO2 stored"], p_nom_extendable=True, capital_cost=costs.at["BioSNG", "fixed"], - marginal_cost=costs.at["BioSNG", "efficiency"] * costs.loc["BioSNG", "VOM"], + marginal_cost=costs.at["BioSNG", "efficiency"] * costs.at["BioSNG", "VOM"], ) # TODO: Update with energy penalty for CC @@ -2379,7 +2414,7 @@ def add_biomass(n, costs): capital_cost=costs.at["BioSNG", "fixed"] + costs.at["biomass CHP capture", "fixed"] * costs.at["BioSNG", "CO2 stored"], - marginal_cost=costs.at["BioSNG", "efficiency"] * costs.loc["BioSNG", "VOM"], + marginal_cost=costs.at["BioSNG", "efficiency"] * costs.at["BioSNG", "VOM"], ) @@ -2612,6 +2647,8 @@ def add_industry(n, costs): p_min_pu=options.get("min_part_load_methanolisation", 0), capital_cost=costs.at["methanolisation", "fixed"] * options["MWh_MeOH_per_MWh_H2"], # EUR/MW_H2/a + marginal_cost=options["MWh_MeOH_per_MWh_H2"] + * costs.at["methanolisation", "VOM"], lifetime=costs.at["methanolisation", "lifetime"], efficiency=options["MWh_MeOH_per_MWh_H2"], efficiency2=-options["MWh_MeOH_per_MWh_H2"] / options["MWh_MeOH_per_MWh_e"], @@ -2729,6 +2766,8 @@ def add_industry(n, costs): efficiency=costs.at["Fischer-Tropsch", "efficiency"], capital_cost=costs.at["Fischer-Tropsch", "fixed"] * costs.at["Fischer-Tropsch", "efficiency"], # EUR/MW_H2/a + marginal_cost=costs.at["Fischer-Tropsch", "efficiency"] + * costs.at["Fischer-Tropsch", "VOM"], efficiency2=-costs.at["oil", "CO2 intensity"] * costs.at["Fischer-Tropsch", "efficiency"], p_nom_extendable=True, @@ -2934,8 +2973,13 @@ def add_waste_heat(n): if not urban_central.empty: urban_central = urban_central.str[: -len(" urban central heat")] + link_carriers = n.links.carrier.unique() + # TODO what is the 0.95 and should it be a config option? - if options["use_fischer_tropsch_waste_heat"]: + if ( + options["use_fischer_tropsch_waste_heat"] + and "Fischer-Tropsch" in link_carriers + ): n.links.loc[urban_central + " Fischer-Tropsch", "bus3"] = ( urban_central + " urban central heat" ) @@ -2943,8 +2987,48 @@ def add_waste_heat(n): 0.95 - n.links.loc[urban_central + " Fischer-Tropsch", "efficiency"] ) + if options["use_methanation_waste_heat"] and "Sabatier" in link_carriers: + n.links.loc[urban_central + " Sabatier", "bus3"] = ( + urban_central + " urban central heat" + ) + n.links.loc[urban_central + " Sabatier", "efficiency3"] = ( + 0.95 - n.links.loc[urban_central + " Sabatier", "efficiency"] + ) + + # DEA quotes 15% of total input (11% of which are high-value heat) + if options["use_haber_bosch_waste_heat"] and "Haber-Bosch" in link_carriers: + n.links.loc[urban_central + " Haber-Bosch", "bus3"] = ( + urban_central + " urban central heat" + ) + total_energy_input = ( + cf_industry["MWh_H2_per_tNH3_electrolysis"] + + cf_industry["MWh_elec_per_tNH3_electrolysis"] + ) / cf_industry["MWh_NH3_per_tNH3"] + electricity_input = ( + cf_industry["MWh_elec_per_tNH3_electrolysis"] + / cf_industry["MWh_NH3_per_tNH3"] + ) + n.links.loc[urban_central + " Haber-Bosch", "efficiency3"] = ( + 0.15 * total_energy_input / electricity_input + ) + + if ( + options["use_methanolisation_waste_heat"] + and "methanolisation" in link_carriers + ): + n.links.loc[urban_central + " methanolisation", "bus4"] = ( + urban_central + " urban central heat" + ) + n.links.loc[urban_central + " methanolisation", "efficiency4"] = ( + costs.at["methanolisation", "heat-output"] + / costs.at["methanolisation", "hydrogen-input"] + ) + # TODO integrate usable waste heat efficiency into technology-data from DEA - if options.get("use_electrolysis_waste_heat", False): + if ( + options.get("use_electrolysis_waste_heat", False) + and "H2 Electrolysis" in link_carriers + ): n.links.loc[urban_central + " H2 Electrolysis", "bus2"] = ( urban_central + " urban central heat" ) @@ -2952,7 +3036,7 @@ def add_waste_heat(n): 0.84 - n.links.loc[urban_central + " H2 Electrolysis", "efficiency"] ) - if options["use_fuel_cell_waste_heat"]: + if options["use_fuel_cell_waste_heat"] and "H2 Fuel Cell" in link_carriers: n.links.loc[urban_central + " H2 Fuel Cell", "bus2"] = ( urban_central + " urban central heat" ) @@ -3310,6 +3394,57 @@ def set_temporal_aggregation(n, opts, solver_name): return n +def lossy_bidirectional_links(n, carrier, efficiencies={}): + "Split bidirectional links into two unidirectional links to include transmission losses." + + carrier_i = n.links.query("carrier == @carrier").index + + if ( + not any((v != 1.0) or (v >= 0) for v in efficiencies.values()) + or carrier_i.empty + ): + return + + efficiency_static = efficiencies.get("efficiency_static", 1) + efficiency_per_1000km = efficiencies.get("efficiency_per_1000km", 1) + compression_per_1000km = efficiencies.get("compression_per_1000km", 0) + + logger.info( + f"Specified losses for {carrier} transmission " + f"(static: {efficiency_static}, per 1000km: {efficiency_per_1000km}, compression per 1000km: {compression_per_1000km}). " + "Splitting bidirectional links." + ) + + n.links.loc[carrier_i, "p_min_pu"] = 0 + n.links.loc[ + carrier_i, "efficiency" + ] = efficiency_static * efficiency_per_1000km ** ( + n.links.loc[carrier_i, "length"] / 1e3 + ) + rev_links = ( + n.links.loc[carrier_i].copy().rename({"bus0": "bus1", "bus1": "bus0"}, axis=1) + ) + rev_links["length_original"] = rev_links["length"] + rev_links["capital_cost"] = 0 + rev_links["length"] = 0 + rev_links["reversed"] = True + rev_links.index = rev_links.index.map(lambda x: x + "-reversed") + + n.links = pd.concat([n.links, rev_links], sort=False) + n.links["reversed"] = n.links["reversed"].fillna(False) + n.links["length_original"] = n.links["length_original"].fillna(n.links.length) + + # do compression losses after concatenation to take electricity consumption at bus0 in either direction + carrier_i = n.links.query("carrier == @carrier").index + if compression_per_1000km > 0: + n.links.loc[carrier_i, "bus2"] = n.links.loc[carrier_i, "bus0"].map( + n.buses.location + ) # electricity + n.links.loc[carrier_i, "efficiency2"] = ( + -compression_per_1000km * n.links.loc[carrier_i, "length_original"] / 1e3 + ) + + if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake @@ -3387,6 +3522,15 @@ def set_temporal_aggregation(n, opts, solver_name): if "nodistrict" in opts: options["district_heating"]["progress"] = 0.0 + if "nowasteheat" in opts: + logger.info("Disabling waste heat.") + options["use_fischer_tropsch_waste_heat"] = False + options["use_methanolisation_waste_heat"] = False + options["use_haber_bosch_waste_heat"] = False + options["use_methanation_waste_heat"] = False + options["use_fuel_cell_waste_heat"] = False + options["use_electrolysis_waste_heat"] = False + if "T" in opts: add_land_transport(n, costs) @@ -3476,6 +3620,18 @@ def set_temporal_aggregation(n, opts, solver_name): if options["electricity_grid_connection"]: add_electricity_grid_connection(n, costs) + for k, v in options["transmission_efficiency"].items(): + lossy_bidirectional_links(n, k, v) + + # Workaround: Remove lines with conflicting (and unrealistic) properties + # cf. https://github.com/PyPSA/pypsa-eur/issues/444 + if snakemake.config["solving"]["options"]["transmission_losses"]: + idx = n.lines.query("num_parallel == 0").index + logger.info( + f"Removing {len(idx)} line(s) with properties conflicting with transmission losses functionality." + ) + n.mremove("Line", idx) + first_year_myopic = (snakemake.params.foresight in ["myopic", "perfect"]) and ( snakemake.params.planning_horizons[0] == investment_year ) diff --git a/scripts/retrieve_databundle.py b/scripts/retrieve_databundle.py index 75d8519e1..25894063c 100644 --- a/scripts/retrieve_databundle.py +++ b/scripts/retrieve_databundle.py @@ -36,7 +36,7 @@ import tarfile from pathlib import Path -from _helpers import configure_logging, progress_retrieve +from _helpers import configure_logging, progress_retrieve, validate_checksum logger = logging.getLogger(__name__) @@ -65,6 +65,8 @@ disable_progress = snakemake.config["run"].get("disable_progressbar", False) progress_retrieve(url, tarball_fn, disable=disable_progress) + validate_checksum(tarball_fn, url) + logger.info("Extracting databundle.") tarfile.open(tarball_fn).extractall(to_fn) diff --git a/scripts/retrieve_gas_infrastructure_data.py b/scripts/retrieve_gas_infrastructure_data.py index 42b726dbd..d984b9feb 100644 --- a/scripts/retrieve_gas_infrastructure_data.py +++ b/scripts/retrieve_gas_infrastructure_data.py @@ -11,7 +11,7 @@ import zipfile from pathlib import Path -from _helpers import progress_retrieve +from _helpers import progress_retrieve, validate_checksum logger = logging.getLogger(__name__) @@ -35,6 +35,8 @@ disable_progress = snakemake.config["run"].get("disable_progressbar", False) progress_retrieve(url, zip_fn, disable=disable_progress) + validate_checksum(zip_fn, url) + logger.info("Extracting databundle.") zipfile.ZipFile(zip_fn).extractall(to_fn) diff --git a/scripts/retrieve_sector_databundle.py b/scripts/retrieve_sector_databundle.py index 0d172c8d1..cb6cc969f 100644 --- a/scripts/retrieve_sector_databundle.py +++ b/scripts/retrieve_sector_databundle.py @@ -13,7 +13,7 @@ import tarfile from pathlib import Path -from _helpers import configure_logging, progress_retrieve +from _helpers import configure_logging, progress_retrieve, validate_checksum if __name__ == "__main__": if "snakemake" not in globals(): @@ -34,6 +34,8 @@ disable_progress = snakemake.config["run"].get("disable_progressbar", False) progress_retrieve(url, tarball_fn, disable=disable_progress) + validate_checksum(tarball_fn, url) + logger.info("Extracting databundle.") tarfile.open(tarball_fn).extractall(to_fn) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index e09bda9bd..5a0455771 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -689,6 +689,35 @@ def add_battery_constraints(n): n.model.add_constraints(lhs == 0, name="Link-charger_ratio") +def add_lossy_bidirectional_link_constraints(n): + if not n.links.p_nom_extendable.any() or not "reversed" in n.links.columns: + return + + n.links["reversed"] = n.links.reversed.fillna(0).astype(bool) + carriers = n.links.loc[n.links.reversed, "carrier"].unique() + + forward_i = n.links.query( + "carrier in @carriers and ~reversed and p_nom_extendable" + ).index + + def get_backward_i(forward_i): + return pd.Index( + [ + re.sub(r"-(\d{4})$", r"-reversed-\1", s) + if re.search(r"-\d{4}$", s) + else s + "-reversed" + for s in forward_i + ] + ) + + backward_i = get_backward_i(forward_i) + + lhs = n.model["Link-p_nom"].loc[backward_i] + rhs = n.model["Link-p_nom"].loc[forward_i] + + n.model.add_constraints(lhs == rhs, name="Link-bidirectional_sync") + + def add_chp_constraints(n): electric = ( n.links.index.str.contains("urban central") @@ -747,9 +776,13 @@ def add_pipe_retrofit_constraint(n): """ Add constraint for retrofitting existing CH4 pipelines to H2 pipelines. """ - gas_pipes_i = n.links.query("carrier == 'gas pipeline' and p_nom_extendable").index + if "reversed" not in n.links.columns: + n.links["reversed"] = False + gas_pipes_i = n.links.query( + "carrier == 'gas pipeline' and p_nom_extendable and ~reversed" + ).index h2_retrofitted_i = n.links.query( - "carrier == 'H2 pipeline retrofitted' and p_nom_extendable" + "carrier == 'H2 pipeline retrofitted' and p_nom_extendable and ~reversed" ).index if h2_retrofitted_i.empty or gas_pipes_i.empty: @@ -788,6 +821,7 @@ def extra_functionality(n, snapshots): if "EQ" in o: add_EQ_constraints(n, o) add_battery_constraints(n) + add_lossy_bidirectional_link_constraints(n) add_pipe_retrofit_constraint(n) if n._multi_invest: add_carbon_constraint(n, snapshots) @@ -848,6 +882,9 @@ def solve_network(n, config, solving, opts="", **kwargs): f"Solving status '{status}' with termination condition '{condition}'" ) if "infeasible" in condition: + labels = n.model.compute_infeasibilities() + logger.info("Labels:\n" + labels) + n.model.print_infeasibilities() raise RuntimeError("Solving status 'infeasible'") return n