Lossy Bidirectional Links

Snakefile

@@ -732,7 +732,10 @@ if config["augmented_line_connection"].get("add_to_snakefile", False) == False:
 
 
 rule add_extra_components:
+    params:
+        transmission_efficiency=config["sector"]["transmission_efficiency"],
     input:
+        overrides="data/override_component_attrs",
         network="networks/" + RDIR + "elec_s{simpl}_{clusters}.nc",
         tech_costs=COSTS,
     output:
@@ -806,6 +809,12 @@ if config["monte_carlo"]["options"].get("add_to_snakefile", False) == False:
             solving=config["solving"],
             augmented_line_connection=config["augmented_line_connection"],
         input:
+            overrides=(
+                os.getcwd() + "/data/override_component_attrs"
+                if os.name == "nt"
+                else "data/override_component_attrs"
+            ),
+            # on Windows os.getcwd() is required because of the "copy-minimal" shadow directory
             network="networks/" + RDIR + "elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
         output:
             "results/" + RDIR + "networks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}.nc",
@@ -872,6 +881,12 @@ if config["monte_carlo"]["options"].get("add_to_snakefile", False) == True:
             solving=config["solving"],
             augmented_line_connection=config["augmented_line_connection"],
         input:
+            overrides=(
+                os.getcwd() + "/data/override_component_attrs"
+                if os.name == "nt"
+                else "data/override_component_attrs"
+            ),
+            # on Windows os.getcwd() is required because of the "copy-minimal" shadow directory
             network="networks/"
             + RDIR
             + "elec_s{simpl}_{clusters}_ec_l{ll}_{opts}_{unc}.nc",
@@ -1615,7 +1630,12 @@ if config["foresight"] == "overnight":
             solving=config["solving"],
             augmented_line_connection=config["augmented_line_connection"],
         input:
-            overrides="data/override_component_attrs",
+            overrides=(
+                os.getcwd() + "/data/override_component_attrs"
+                if os.name == "nt"
+                else "data/override_component_attrs"
+            ),
+            # on Windows os.getcwd() is required because of the "copy-minimal" shadow directory
             # network=RESDIR
             # + "prenetworks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}_{sopts}_{planning_horizons}_{discountrate}.nc",
             network=RESDIR
@@ -2081,7 +2101,12 @@ if config["foresight"] == "myopic":
                 "co2_sequestration_potential", 200
             ),
         input:
-            overrides="data/override_component_attrs",
+            overrides=(
+                os.getcwd() + "/data/override_component_attrs"
+                if os.name == "nt"
+                else "data/override_component_attrs"
+            ),
+            # on Windows os.getcwd() is required because of the "copy-minimal" shadow directory
             network=RESDIR
             + "prenetworks-brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sopts}_{planning_horizons}_{discountrate}_{demand}_{h2export}export.nc",
             costs=CDIR + "costs_{planning_horizons}.csv",

config.default.yaml

@@ -602,6 +602,9 @@ sector:
   transmission_efficiency:
     electricity distribution grid:
       efficiency_static: 0.97 # efficiency of distribution grid (i.e. 3% loses)
+    H2 pipeline:
+      efficiency_per_1000km: 1
+      compression_per_1000km: 0.017 # DEA technology data. Mean of  Energy losses, lines 5000-20000 MW and lines >20000 MW for 2020, 2030 and 2050, [%/1000 km]
 
   dynamic_transport:
     enable: false # If "True", then the BEV and FCEV shares are obtained depending on the "Co2L"-wildcard (e.g. "Co2L0.70: 0.10"). If "False", then the shares are obtained depending on the "demand" wildcard and "planning_horizons" wildcard as listed below (e.g. "DF_2050: 0.08")
@@ -695,7 +698,6 @@ sector:
     biomass: biomass
   keep_existing_capacities: true
 
-
 solving:
   options:
     formulation: kirchhoff

doc/release_notes.rst

@@ -17,6 +17,8 @@ This part of documentation collects descriptive release notes to capture the mai
 
 * In alternative clustering, generate hydro inflows by shape and avoid hydro inflows duplication for plants installed in the same node `PR #1120 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1120>`
 
+* Add a function to calculate length-based efficiencies and apply it to the H2 pipelines. `PR #1192 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1192>`__
+
 **Minor Changes and bug-fixing**
 
 * Prevent computation of powerplantmatching if replace option is selected for custom_powerplants `PR #1281 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1281>`__
@@ -25,6 +27,8 @@ This part of documentation collects descriptive release notes to capture the mai
 
 * Fix readthedocs by explicitly specifying the location of the Sphinx config `PR #1292 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1292>`__
 
+* Fix lossy bidirectional links, especially H2 pipelines, which would sometimes gain H2 instead of losing it.  `PR #1192 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1192>`__
+
 PyPSA-Earth 0.6.0
 =================
 

scripts/_helpers.py

@@ -1483,3 +1483,91 @@ def safe_divide(numerator, denominator, default_value=np.nan):
             f"Division by zero: {numerator} / {denominator}, returning NaN."
         )
         return np.nan
+
+
+def lossy_bidirectional_links(n, carrier):
+    """
+    Split bidirectional links of type carrier into two unidirectional links to include transmission losses.
+    """
+
+    # identify all links of type carrier
+    carrier_i = n.links.query("carrier == @carrier").index
+
+    if carrier_i.empty:
+        return
+
+    logger.info(f"Splitting bidirectional links with the carrier {carrier}")
+
+    # set original links to be unidirectional
+    n.links.loc[carrier_i, "p_min_pu"] = 0
+
+    # add a new links that mirror the original links, but represent the reversed flow direction
+    # the new links have a cost and length of 0 to not distort the overall cost and network length
+    rev_links = (
+        n.links.loc[carrier_i].copy().rename({"bus0": "bus1", "bus1": "bus0"}, axis=1)
+    )
+    rev_links["length_original"] = rev_links[
+        "length"
+    ]  # tracker for the length of the original links length
+    rev_links["capital_cost"] = 0
+    rev_links["length"] = 0
+    rev_links["reversed"] = True  # tracker for easy identification of reversed links
+    rev_links.index = rev_links.index.map(lambda x: x + "-reversed")
+
+    # add the new reversed links to the network and fill the newly created trackers with default values for the other links
+    n.links = pd.concat([n.links, rev_links], sort=False)
+    n.links["reversed"] = n.links["reversed"].fillna(False).infer_objects(copy=False)
+    n.links["length_original"] = n.links["length_original"].fillna(n.links.length)
+
+
+def set_length_based_efficiency(n, carrier, bus_suffix, transmission_efficiency):
+    """
+    Set the efficiency of all links of type carrier in network n based on their length and the values specified in the config.
+    Additionally add the length based electricity demand required for compression (if applicable).
+    The bus_suffix refers to the suffix that differentiates the links bus0 from the corresponding electricity bus, i.e. " H2".
+    Important:
+    Call this function AFTER lossy_bidirectional_links when creating links that are both bidirectional and lossy,
+    and have a length based electricity demand for compression. Otherwise the compression will not consistently take place at
+    the inflow bus and instead vary between the inflow and the outflow bus.
+    """
+
+    # get the links length based efficiency and required compression
+    if carrier not in transmission_efficiency:
+        raise KeyError(
+            f"An error occurred when setting the length based efficiency for the Links of type {carrier}."
+            f"The Link type {carrier} was not found in the config under config['sector']['transmission_efficiency']."
+        )
+    efficiencies = transmission_efficiency[carrier]
+    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)
+
+    # indetify all links of type carrier
+    carrier_i = n.links.loc[n.links.carrier == carrier].index
+
+    # identify the lengths of all links of type carrier
+    # use "length_original" for lossy bidirectional links and "length" for any other link
+    if ("reversed" in n.links.columns) and any(n.links.loc[carrier_i, "reversed"]):
+        lengths = n.links.loc[carrier_i, "length_original"]
+    else:
+        lengths = n.links.loc[carrier_i, "length"]
+
+    # set the links' length based efficiency
+    n.links.loc[carrier_i, "efficiency"] = (
+        efficiency_static * efficiency_per_1000km ** (lengths / 1e3)
+    )
+
+    # set the links's electricity demand for compression
+    if compression_per_1000km > 0:
+        # connect the links to their corresponding electricity buses
+        n.links.loc[carrier_i, "bus2"] = n.links.loc[
+            carrier_i, "bus0"
+        ].str.removesuffix(bus_suffix)
+        # TODO: use these lines to set bus 2 instead, once n.buses.location is functional and remove bus_suffix.
+        """
+        n.links.loc[carrier_i, "bus2"] = n.links.loc[carrier_i, "bus0"].map(
+            n.buses.location
+        )  # electricity
+        """
+        # set the required compression demand
+        n.links.loc[carrier_i, "efficiency2"] = -compression_per_1000km * lengths / 1e3

scripts/add_extra_components.py

@@ -57,7 +57,13 @@
 import numpy as np
 import pandas as pd
 import pypsa
-from _helpers import configure_logging, create_logger
+from _helpers import (
+    configure_logging,
+    create_logger,
+    lossy_bidirectional_links,
+    override_component_attrs,
+    set_length_based_efficiency,
+)
 from add_electricity import (
     _add_missing_carriers_from_costs,
     add_nice_carrier_names,
@@ -225,7 +231,7 @@ def attach_stores(n, costs, config):
         )
 
 
-def attach_hydrogen_pipelines(n, costs, config):
+def attach_hydrogen_pipelines(n, costs, config, transmission_efficiency):
     elec_opts = config["electricity"]
     ext_carriers = elec_opts["extendable_carriers"]
     as_stores = ext_carriers.get("Store", [])
@@ -261,10 +267,15 @@ def attach_hydrogen_pipelines(n, costs, config):
         p_nom_extendable=True,
         length=h2_links.length.values,
         capital_cost=costs.at["H2 pipeline", "capital_cost"] * h2_links.length,
-        efficiency=costs.at["H2 pipeline", "efficiency"],
         carrier="H2 pipeline",
     )
 
+    # split the pipeline into two unidirectional links to properly apply transmission losses in both directions.
+    lossy_bidirectional_links(n, "H2 pipeline")
+
+    # set the pipelines efficiency and the electricity required by the pipeline for compression
+    set_length_based_efficiency(n, "H2 pipeline", " H2", transmission_efficiency)
+
 
 if __name__ == "__main__":
     if "snakemake" not in globals():
@@ -274,8 +285,10 @@ def attach_hydrogen_pipelines(n, costs, config):
 
     configure_logging(snakemake)
 
-    n = pypsa.Network(snakemake.input.network)
+    overrides = override_component_attrs(snakemake.input.overrides)
+    n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides)
     Nyears = n.snapshot_weightings.objective.sum() / 8760.0
+    transmission_efficiency = snakemake.params.transmission_efficiency
     config = snakemake.config
 
     costs = load_costs(
@@ -287,7 +300,7 @@ def attach_hydrogen_pipelines(n, costs, config):
 
     attach_storageunits(n, costs, config)
     attach_stores(n, costs, config)
-    attach_hydrogen_pipelines(n, costs, config)
+    attach_hydrogen_pipelines(n, costs, config, transmission_efficiency)
 
     add_nice_carrier_names(n, config=snakemake.config)
 

scripts/solve_network.py

@@ -854,6 +854,62 @@ def add_existing(n):
             n.generators.loc[tech_index, tech] = existing_res
 
 
+def add_lossy_bidirectional_link_constraints(n: pypsa.components.Network) -> None:
+    """
+    Ensures that the two links simulating a bidirectional_link are extended the same amount.
+    """
+
+    if not n.links.p_nom_extendable.any() or "reversed" not in n.links.columns:
+        return
+
+    # ensure that the 'reversed' column is boolean and identify all link carriers that have 'reversed' links
+    n.links["reversed"] = n.links.reversed.fillna(0).astype(bool)
+    carriers = n.links.loc[n.links.reversed, "carrier"].unique()  # noqa: F841
+
+    # get the indices of all forward links (non-reversed), that have a reversed counterpart
+    forward_i = n.links.loc[
+        n.links.carrier.isin(carriers) & ~n.links.reversed & n.links.p_nom_extendable
+    ].index
+
+    # function to get backward (reversed) indices corresponding to forward links
+    # this function is required to properly interact with the myopic naming scheme
+    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
+            ]
+        )
+
+    # get the indices of all backward links (reversed)
+    backward_i = get_backward_i(forward_i)
+
+    # get the p_nom optimization variables for the links using the get_var function
+    links_p_nom = get_var(n, "Link", "p_nom")
+
+    # only consider forward and backward links that are present in the optimization variables
+    subset_forward = forward_i.intersection(links_p_nom.index)
+    subset_backward = backward_i.intersection(links_p_nom.index)
+
+    # ensure we have a matching number of forward and backward links
+    if len(subset_forward) != len(subset_backward):
+        raise ValueError("Mismatch between forward and backward links.")
+
+    # define the lefthand side of the constrain p_nom (forward) - p_nom (backward) = 0
+    # this ensures that the forward links always have the same maximum nominal power as their backward counterpart
+    lhs = linexpr(
+        (1, get_var(n, "Link", "p_nom")[backward_i].to_numpy()),
+        (-1, get_var(n, "Link", "p_nom")[forward_i].to_numpy()),
+    )
+
+    # add the constraint to the PySPA model
+    define_constraints(n, lhs, "=", 0, "Link-bidirectional_sync")
+
+
 def extra_functionality(n, snapshots):
     """
     Collects supplementary constraints which will be passed to
@@ -881,6 +937,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)
 
     if (
         snakemake.config["policy_config"]["hydrogen"]["temporal_matching"]
@@ -986,11 +1043,8 @@ def solve_network(n, config, solving={}, opts="", **kwargs):
 
     is_sector_coupled = "sopts" in snakemake.wildcards.keys()
 
-    if is_sector_coupled:
-        overrides = override_component_attrs(snakemake.input.overrides)
-        n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides)
-    else:
-        n = pypsa.Network(snakemake.input.network)
+    overrides = override_component_attrs(snakemake.input.overrides)
+    n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides)
 
     if snakemake.params.augmented_line_connection.get("add_to_snakefile"):
         n.lines.loc[n.lines.index.str.contains("new"), "s_nom_min"] = (

test/config.sector.yaml

@@ -20,8 +20,8 @@ countries: ["NG", "BJ"]
 
 electricity:
   extendable_carriers:
-    Store: []
-    Link: []
+    Store: [H2]
+    Link: [H2 pipeline]
 
   co2limit: 7.75e7