Limit H2 feedin to a share of CH4 load

etrago/appl.py

@@ -75,6 +75,7 @@
     # Scenario variations:
     "scn_extension": None,  # None or array of extension scenarios
     "scn_decommissioning": None,  # None or decommissioning scenario
+    "H2_vol_share": 15,  # in % [50/20/15/10/5/2/1/0] allowed H2 volumetric share for feedin
     # Export options:
     "lpfile": False,  # save pyomo's lp file: False or /path/to/lpfile.lp
     "csv_export": "results",  # save results as csv: False or /path/tofolder
@@ -243,6 +244,13 @@ def run_etrago(args, json_path):
             'nep2035_b2' includes all lines that will be replaced in
             NEP-scenario 2035 B2
 
+    H2_vol_share : int
+        15,
+        Allowed H2 volumetric share of the CH4 loads that could be fed
+        into the CH4 grid if H2_feedin links are present in the network
+        Possible values are: [50/20/15/10/5/2/1/0]
+        If 0 is set, the H2_feedin links are deleted of the network.
+
     lpfile : obj
         False,
         State if and where you want to save pyomo's lp file. Options:

etrago/cluster/gas.py

@@ -7,21 +7,21 @@
 
 if "READTHEDOCS" not in os.environ:
 
-    import numpy as np
-    import pandas as pd
-    import pypsa.io as io
     from pypsa import Network
     from pypsa.networkclustering import (
         aggregatebuses,
         aggregateoneport,
         busmap_by_kmeans,
     )
     from six import iteritems
+    import numpy as np
+    import pandas as pd
+    import pypsa.io as io
 
     from etrago.cluster.spatial import (
-        sum_with_inf,
         group_links,
         kmedoids_dijkstra_clustering,
+        sum_with_inf,
     )
     from etrago.tools.utilities import *
 
@@ -126,10 +126,13 @@ def weighting_for_scenario(ch4_buses, save=None):
             ].index
         # get all generators and loads related to ch4_buses
         generators_ = pd.Series(
-            etrago.network.generators.index, index=etrago.network.generators.bus
+            etrago.network.generators.index,
+            index=etrago.network.generators.bus,
         )
         buses_CH4_gen = generators_.index.intersection(rel_links.keys())
-        loads_ = pd.Series(etrago.network.loads.index, index=etrago.network.loads.bus)
+        loads_ = pd.Series(
+            etrago.network.loads.index, index=etrago.network.loads.bus
+        )
         buses_CH4_load = loads_.index.intersection(rel_links.keys())
 
         # sum up all relevant entities and cast to integer
@@ -142,7 +145,9 @@ def weighting_for_scenario(ch4_buses, save=None):
                 ].p_nom.sum()
             if i in buses_CH4_load:
                 rel_links[i] += (
-                    etrago.network.loads_t.p_set.loc[:, loads_.loc[i]].mean().sum()
+                    etrago.network.loads_t.p_set.loc[:, loads_.loc[i]]
+                    .mean()
+                    .sum()
                 )
             rel_links[i] = min(int(rel_links[i]), MAX_WEIGHT)
         weightings = pd.DataFrame.from_dict(rel_links, orient="index")
@@ -161,7 +166,9 @@ def weighting_for_scenario(ch4_buses, save=None):
     elif settings["gas_weight_fromcsv"] is not None:
         # create DataFrame with uniform weightings for all ch4_buses
         weight_ch4 = pd.DataFrame([1] * len(buses_ch4), index=buses_ch4.index)
-        loaded_weights = pd.read_csv(settings["gas_weight_fromcsv"], index_col=0)
+        loaded_weights = pd.read_csv(
+            settings["gas_weight_fromcsv"], index_col=0
+        )
         # load weights into previously created DataFrame
         loaded_weights.index = loaded_weights.index.astype(str)
         weight_ch4.loc[loaded_weights.index] = loaded_weights
@@ -190,7 +197,9 @@ def kmean_clustering_gas(etrago, network_ch4, weight, n_clusters):
         busmap_ch4.index.name = 'bus_id'
         busmap_ch4.name = 'cluster'
         busmap_ch4.to_csv(
-            "kmeans_ch4_busmap_" + str(settings["n_clusters_gas"]) + "_result.csv"
+            "kmeans_ch4_busmap_"
+            + str(settings["n_clusters_gas"])
+            + "_result.csv"
         )
     else:
 
@@ -203,9 +212,12 @@ def kmean_clustering_gas(etrago, network_ch4, weight, n_clusters):
 
 
 def get_h2_clusters(etrago, busmap_ch4):
-    
+
     # Mapping of H2 buses to new CH4 cluster IDs
-    busmap_h2 = pd.Series(busmap_ch4.loc[etrago.ch4_h2_mapping.index].values, index = etrago.ch4_h2_mapping.values)
+    busmap_h2 = pd.Series(
+        busmap_ch4.loc[etrago.ch4_h2_mapping.index].values,
+        index=etrago.ch4_h2_mapping.values,
+    )
 
     # Create unique H2 cluster IDs
     n_gas = etrago.args["network_clustering"]["n_clusters_gas"]
@@ -226,7 +238,9 @@ def gas_postprocessing(etrago, busmap, medoid_idx):
     # Add all other buses to busmap
 
     missing_idx = list(
-        etrago.network.buses[(~etrago.network.buses.index.isin(busmap.index))].index
+        etrago.network.buses[
+            (~etrago.network.buses.index.isin(busmap.index))
+        ].index
     )
     next_bus_id = highestInteger(etrago.network.buses.index) + 1
     new_gas_buses = [str(int(x) + next_bus_id) for x in busmap]
@@ -309,33 +323,71 @@ def gas_postprocessing(etrago, busmap, medoid_idx):
     )
 
     # aggregation of the links and links time series
-    network_gasgrid_c.links, network_gasgrid_c.links_t = group_links(network_gasgrid_c)
+    network_gasgrid_c.links, network_gasgrid_c.links_t = group_links(
+        network_gasgrid_c
+    )
 
-    # Overwrite p_nom of links with carrier "H2_feedin" (eGon2035 only)
-    if etrago.args["scn_name"] == "eGon2035":
-        H2_energy_share = 0.05053  # H2 energy share via volumetric share outsourced in a mixture of H2 and CH4 with 15 %vol share
-        feed_in = network_gasgrid_c.links.loc[
-            network_gasgrid_c.links.carrier == "H2_feedin"
-        ]
-        pipeline_capacities = network_gasgrid_c.links.loc[
-            network_gasgrid_c.links.carrier == "CH4"
+    # Overwrite p_nom of links with carrier "H2_feedin"
+    if "H2_feedin" in network_gasgrid_c.links.carrier.to_list():
+
+        H2_vol_share = etrago.args["H2_vol_share"]
+
+        def att_H2_energy_share(H2_vol_share):
+            """
+            Return the fraction of H2 with respect to energy in a H2 CH4 mixture
+
+            The calculation of the values in the dictionary has been
+            made using the function H2_CH4_mix_energy_fractions of
+            https://github.com/openego/eGon-data/blob/dev/src/egon/data/datasets/hydrogen_etrago/h2_to_ch4.py
+            with T=25 (°C) and p=50 (bar).
+
+            Parameters
+            ----------
+            H2_vol_share : float
+                Volumetric fraction of H2 in the mixture
+
+            Returns
+            -------
+            H2_vol2en[H2_vol_share] : float
+                Fraction of H2 in mixture with respect to energy (LHV)
+
+            """
+            H2_vol2en = {
+                1: 0.00304,
+                2: 0.00612,
+                5: 0.01562,
+                10: 0.03242,
+                15: 0.05053,
+                20: 0.07011,
+                50: 0.23170,
+            }
+            return H2_vol2en[H2_vol_share]
+
+        rel_ch4_loads_carriers = ["rural_gas_boiler", "CH4_for_industry"]
+        ch4_loads = network_gasgrid_c.loads.loc[
+            network_gasgrid_c.loads.carrier.isin(rel_ch4_loads_carriers)
         ]
 
-        for bus in feed_in["bus1"].values:
-            # calculate the total pipeline capacity connected to a specific bus
-            nodal_capacity = pipeline_capacities.loc[
-                (pipeline_capacities["bus0"] == bus)
-                | (pipeline_capacities["bus1"] == bus),
-                "p_nom",
-            ].sum()
-            # multiply total pipeline capacity with H2 energy share corresponding to volumetric share
+        for bus in ch4_loads["bus"].unique():
+            # Calculate the sum of the CH4 loads at each CH4 bus (this sum is a time series)
+            load_names = ch4_loads.loc[ch4_loads["bus"].values == bus].index
+            ch4_loads_set = network_gasgrid_c.loads_t.p_set.loc[
+                :, load_names
+            ].sum(axis=1)
+            # Overwrite the "p_nom" values of the H2_link with a share of the CH4 total load
+            feedin_link = network_gasgrid_c.links.loc[
+                (network_gasgrid_c.links["carrier"].values == "H2_feedin")
+                & (network_gasgrid_c.links["bus1"].values == bus)
+            ]
+            # p_nom is the max of the share of the CH4 total load
             network_gasgrid_c.links.loc[
-                (network_gasgrid_c.links["bus1"].values == bus)
-                & (network_gasgrid_c.links["carrier"].values == "H2_feedin"),
-                "p_nom",
-            ] = (
-                nodal_capacity * H2_energy_share
-            )
+                feedin_link.index, "p_nom"
+            ] = ch4_loads_set.max() * att_H2_energy_share(H2_vol_share)
+            # and p_max_pu the time serie of the total CH4 load, normalized by its max value
+            network_gasgrid_c.links_t.p_max_pu.loc[
+                ch4_loads_set.index, feedin_link.index
+            ] = (ch4_loads_set / ch4_loads_set.max())
+
     # Insert components not related to the gas clustering
     other_components = ["Line", "StorageUnit", "ShuntImpedance", "Transformer"]
 
@@ -366,11 +418,19 @@ def gas_postprocessing(etrago, busmap, medoid_idx):
         ].index:
             cluster = str(i)
             if cluster in busmap[medoid_idx].values:
-                medoid = busmap[medoid_idx][busmap[medoid_idx] == cluster].index
+                medoid = busmap[medoid_idx][
+                    busmap[medoid_idx] == cluster
+                ].index
                 h2_idx = network_gasgrid_c.buses.loc[
                     (network_gasgrid_c.buses.carrier == "H2_grid")
-                    & (network_gasgrid_c.buses.y == network_gasgrid_c.buses.at[i, "y"])
-                    & (network_gasgrid_c.buses.x == network_gasgrid_c.buses.at[i, "x"])
+                    & (
+                        network_gasgrid_c.buses.y
+                        == network_gasgrid_c.buses.at[i, "y"]
+                    )
+                    & (
+                        network_gasgrid_c.buses.x
+                        == network_gasgrid_c.buses.at[i, "x"]
+                    )
                 ]
                 if len(h2_idx) > 0:
                     h2_idx = h2_idx.index.tolist()[0]
@@ -413,7 +473,9 @@ def highestInteger(potentially_numbers):
     return highest
 
 
-def simultaneous_sector_coupling(network, busmap, carrier_based, carrier_to_cluster):
+def simultaneous_sector_coupling(
+    network, busmap, carrier_based, carrier_to_cluster
+):
     """Cluster sector coupling technology based on multiple connected carriers.
 
     The topology of the sector coupling technology must be in a way, that the
@@ -438,8 +500,12 @@ def simultaneous_sector_coupling(network, busmap, carrier_based, carrier_to_clus
         Busmap for the sector coupling cluster.
     """
     next_bus_id = highestInteger(busmap.values) + 1
-    buses_clustered = network.buses[network.buses["carrier"].isin(carrier_based)]
-    buses_to_cluster = network.buses[network.buses["carrier"] == carrier_to_cluster]
+    buses_clustered = network.buses[
+        network.buses["carrier"].isin(carrier_based)
+    ]
+    buses_to_cluster = network.buses[
+        network.buses["carrier"] == carrier_to_cluster
+    ]
     buses_to_skip = network.buses[
         network.buses["carrier"] == carrier_to_cluster + "_store"
     ]
@@ -461,7 +527,9 @@ def simultaneous_sector_coupling(network, busmap, carrier_based, carrier_to_clus
 
     # cluster sector coupling technologies
     busmap = sc_multi_carrier_based(buses_to_cluster, connected_links)
-    busmap = {bus_id: bus_num + next_bus_id for bus_id, bus_num in busmap.items()}
+    busmap = {
+        bus_id: bus_num + next_bus_id for bus_id, bus_num in busmap.items()
+    }
 
     # cluster appedices
     skipped_links = network.links.loc[
@@ -502,7 +570,9 @@ def simultaneous_sector_coupling(network, busmap, carrier_based, carrier_to_clus
     return busmap
 
 
-def consecutive_sector_coupling(network, busmap, carrier_based, carrier_to_cluster):
+def consecutive_sector_coupling(
+    network, busmap, carrier_based, carrier_to_cluster
+):
     """Cluster sector coupling technology based on single connected carriers.
 
     The topology of the sector coupling technology must be in a way, that the
@@ -530,8 +600,12 @@ def consecutive_sector_coupling(network, busmap, carrier_based, carrier_to_clust
     buses_to_skip = network.buses[
         network.buses["carrier"] == carrier_to_cluster + "_store"
     ]
-    buses_to_cluster = network.buses[network.buses["carrier"] == carrier_to_cluster]
-    buses_clustered = network.buses[network.buses["carrier"] == carrier_based[0]]
+    buses_to_cluster = network.buses[
+        network.buses["carrier"] == carrier_to_cluster
+    ]
+    buses_clustered = network.buses[
+        network.buses["carrier"] == carrier_based[0]
+    ]
     busmap_sc = {}
 
     for base in carrier_based:
@@ -632,7 +706,8 @@ def sc_multi_carrier_based(buses_to_cluster, connected_links):
         clusters.loc[bus_id] = tuple(
             sorted(
                 connected_links.loc[
-                    connected_links["bus1_clustered"] == bus_id, "bus0_clustered"
+                    connected_links["bus1_clustered"] == bus_id,
+                    "bus0_clustered",
                 ].unique()
             )
         )
@@ -705,13 +780,16 @@ def get_clustering_from_busmap(
             with_time=with_time,
             custom_strategies=one_port_strategies.get(one_port, {}),
         )
-        io.import_components_from_dataframe(network_gasgrid_c, new_df, one_port)
+        io.import_components_from_dataframe(
+            network_gasgrid_c, new_df, one_port
+        )
         for attr, df in iteritems(new_pnl):
             io.import_series_from_dataframe(network_gasgrid_c, df, one_port, attr)
     # Aggregate links
     new_links = (
         network.links.assign(
-            bus0=network.links.bus0.map(busmap), bus1=network.links.bus1.map(busmap)
+            bus0=network.links.bus0.map(busmap),
+            bus1=network.links.bus1.map(busmap),
         )
         .dropna(subset=["bus0", "bus1"])
         .loc[lambda df: df.bus0 != df.bus1]
@@ -725,7 +803,9 @@ def get_clustering_from_busmap(
     # bus1=12 and bus0=12, bus1=1) they are aggregated to a single pipeline.
     # therefore, the order of bus0/bus1 is adjusted
     pipeline_mask = new_links["carrier"] == "CH4"
-    sorted_buses = np.sort(new_links.loc[pipeline_mask, ["bus0", "bus1"]].values, 1)
+    sorted_buses = np.sort(
+        new_links.loc[pipeline_mask, ["bus0", "bus1"]].values, 1
+    )
     new_links.loc[pipeline_mask, ["bus0", "bus1"]] = sorted_buses
 
     # import the links and the respective time series with the bus0 and bus1
@@ -749,7 +829,9 @@ def run_spatial_clustering_gas(self):
         gas_network, weight, n_clusters = preprocessing(self)
 
         if method == "kmeans":
-            busmap = kmean_clustering_gas(self, gas_network, weight, n_clusters)
+            busmap = kmean_clustering_gas(
+                self, gas_network, weight, n_clusters
+            )
             medoid_idx = None
         elif method == "kmedoids-dijkstra":
 

etrago/tools/network.py

@@ -64,6 +64,7 @@
     convert_capital_costs,
     crossborder_capacity,
     delete_dispensable_ac_buses,
+    delete_h2_feedin,
     drop_sectors,
     export_to_csv,
     filter_links_by_carrier,
@@ -263,6 +264,8 @@ def __init__(
 
     delete_dispensable_ac_buses = delete_dispensable_ac_buses
 
+    delete_h2_feedin = delete_h2_feedin
+
     get_clustering_data = get_clustering_data
 
     adjust_CH4_gen_carriers = adjust_CH4_gen_carriers
@@ -345,5 +348,7 @@ def adjust_network(self):
 
         self.delete_dispensable_ac_buses()
 
+        self.delete_h2_feedin()
+
     def _ts_weighted(self, timeseries):
         return timeseries.mul(self.network.snapshot_weightings, axis=0)

etrago/tools/utilities.py

@@ -1162,6 +1162,15 @@ def count(bus):
     return
 
 
+def delete_h2_feedin(self):
+    """Delete H2_feedin links if H2_vol_share = 0"""
+
+    if self.args["H2_vol_share"] == 0:
+        self.network.links = self.network.links[
+            self.network.links.carrier != "H2_feedin"
+        ]
+
+
 def set_line_costs(self, cost110=230, cost220=290, cost380=85, costDC=375):
     """Set capital costs for extendable lines in respect to PyPSA [€/MVA]