Add functionality to merge each isolated bus to the backbone network

config.default.yaml

@@ -64,7 +64,8 @@ cluster_options:
     remove_stubs: true
     remove_stubs_across_borders: true
     p_threshold_drop_isolated: 20 # [MW] isolated buses are being discarded if bus mean power is below the specified threshold
-    p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if bus mean power is below the specified threshold
+    p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if a bus mean power is below the specified threshold
+    s_threshold_fetch_isolated: 0.05 # [-] a share of the national load for merging an isolated network into a backbone network
   cluster_network:
     algorithm: kmeans
     feature: solar+onwind-time

config.tutorial.yaml

@@ -77,7 +77,8 @@ cluster_options:
     remove_stubs: true
     remove_stubs_across_borders: true
     p_threshold_drop_isolated: 20 # [MW] isolated buses are being discarded if bus mean power is below the specified threshold
-    p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if bus mean power is below the specified threshold
+    p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if a bus mean power is below the specified threshold
+    s_threshold_fetch_isolated: 0.05 # [-] a share of the national load for merging an isolated network into a backbone network
   cluster_network:
     algorithm: kmeans
     feature: solar+onwind-time

doc/configtables/cluster_options.csv

@@ -8,6 +8,7 @@ simplify_network,,,
 -- remove_stubs_across_borders, bool, "{True, False}", "True: Stub lines and links can be removed across borders."
 -- p_threshold_drop_isolated, MW, positive number, "Isolated buses are discarded if bus mean power is below the `p_threshold_drop_isolated`."
 -- p_threshold_merge_isolated, MW, positive number, "Isolated buses are merged into a single isolated bus if bus mean power is below `p_threshold_merge_isolated`."
+-- s_threshold_fetch_isolated, [-], positive number, "Isolated networks are merged into a backbone network of a respective country if the network load comprises a share of the national load less than p_threshold_fetch_isolated."
 cluster_network,,,
 -- algorithm,,"{hac, kmeans}", "Clustering algorithm used in the cluster_network rule. Options available are Hierarchical Agglomerative Clustering (HAC) or k-means."
 -- feature,,"Str in the format ‘carrier1+carrier2+...+carrierN-X’, where CarrierI can be from {‘solar’, ‘onwind’, ‘offwind’, ‘ror’} and X is one of {‘cap’, ‘time’}. Examples: solar+offwind-cap, solar-time", "Only for Hierarchical Agglomerative Clustering (HAC). Feature(s) used to do the clustering."

doc/release_notes.rst

@@ -20,6 +20,8 @@ E.g. if a new rule becomes available describe how to use it `snakemake -j1 run_t
 
 * Generalize line types for AC and DC networks. `PR #999 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/999>`__
 
+* Add an option to merge isolated networks into respective backbone networks by countries. `PR #903 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/903>`__
+
 **Minor Changes and bug-fixing**
 
 * Minor bug-fixing to get the generalised line types work for DC lines. `PR #1008 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1008>`__

scripts/simplify_network.py

@@ -89,6 +89,7 @@
 import sys
 from functools import reduce
 
+import geopandas as gpd
 import numpy as np
 import pandas as pd
 import pypsa
@@ -743,6 +744,136 @@ def drop_isolated_nodes(n, threshold):
     return n
 
 
+def transform_to_gdf(n, network_crs):
+    buses_df = n.buses.copy()
+    buses_df["bus_id"] = buses_df.index
+
+    # load data are crucial to deal with sub-networks
+    buses_df = buses_df.join(n.loads_t.p_set.sum().T.rename("load"))
+    buses_df["load_in_subnetw"] = buses_df.groupby(["country", "sub_network"])[
+        "load"
+    ].transform("sum")
+    buses_df["load_in_country"] = buses_df.groupby(["country"])["load"].transform("sum")
+    buses_df["sbntw_share_of_country_load"] = buses_df.apply(
+        lambda row: (
+            row.load_in_subnetw / row.load_in_country if row.load_in_country > 0 else 0
+        ),
+        axis=1,
+    )
+    buses_df["is_backbone_sbntw"] = (
+        buses_df.groupby(["country"], as_index=True)[
+            "sbntw_share_of_country_load"
+        ].transform("max")
+        <= buses_df["sbntw_share_of_country_load"]
+    )
+
+    gdf_buses = gpd.GeoDataFrame(
+        buses_df,
+        geometry=gpd.points_from_xy(buses_df.x, buses_df.y),
+        crs=network_crs,
+    )
+    return gdf_buses
+
+
+def merge_into_network(n, threshold, aggregation_strategies=dict()):
+    """
+    Find isolated AC nodes and sub-networks in the network and merge those of
+    them which have load value and a number of buses below than the specified
+    thresholds into a backbone network.
+
+    Parameters
+    ----------
+    n : PyPSA.Network
+        Original network
+    threshold : float
+        Load power used as a threshold to merge isolated nodes
+    aggregation_strategies: dictionary
+        Functions to be applied to calculate parameters of the aggregated grid
+
+    Returns
+    -------
+    modified network
+    """
+    # keep original values of the overall load and generation in the network
+    # to track changes due to drop of buses
+    generators_mean_origin = n.generators.p_nom.mean()
+    load_mean_origin = n.loads_t.p_set.mean().mean()
+
+    network_crs = snakemake.params.geo_crs
+
+    n.determine_network_topology()
+
+    n_buses_gdf = transform_to_gdf(n, network_crs=network_crs)
+
+    # do not merge sub-networks spanned through a number of countries
+    n_buses_gdf["is_multicnt_subntw"] = n_buses_gdf.sub_network.map(
+        n_buses_gdf.groupby(["sub_network"]).country.nunique() > 1
+    )
+
+    gdf_islands = (
+        n_buses_gdf.query("~is_multicnt_subntw")
+        .query("carrier=='AC'")
+        .query("sbntw_share_of_country_load < @threshold")
+    )
+    # return the original network if no isolated nodes are detected
+    if len(gdf_islands) == 0:
+        return n, n.buses.index.to_series()
+
+    gdf_backbone_buses = n_buses_gdf.query("is_backbone_sbntw").query("carrier=='AC'")
+
+    # find the closest buses of the backbone networks for each isolated network and each country
+    islands_bcountry = {k: d for k, d in gdf_islands.groupby("country")}
+    gdf_map = (
+        gdf_backbone_buses.query("country in @islands_bcountry")
+        .groupby("country")
+        .apply(lambda d: gpd.sjoin_nearest(islands_bcountry[d["country"].values[0]], d))
+    )
+    nearest_bus_df = n.buses.loc[n.buses.index.isin(gdf_map.bus_id_right)]
+
+    i_lines_islands = n.lines.loc[n.lines.bus1.isin(gdf_islands.index)].index
+    n.mremove("Line", i_lines_islands)
+
+    isolated_buses_mapping = (
+        gdf_map[["bus_id_right"]].droplevel("country").to_dict()["bus_id_right"]
+    )
+
+    busmap = (
+        n.buses.index.to_series()
+        .replace(isolated_buses_mapping)
+        .astype(str)
+        .rename("busmap")
+    )
+
+    # return the original network if no changes are detected
+    if (busmap.index == busmap).all():
+        return n, n.buses.index.to_series()
+
+    bus_strategies, generator_strategies = get_aggregation_strategies(
+        aggregation_strategies
+    )
+
+    clustering = get_clustering_from_busmap(
+        n,
+        busmap,
+        bus_strategies=bus_strategies,
+        aggregate_generators_weighted=True,
+        aggregate_generators_carriers=None,
+        aggregate_one_ports=["Load", "StorageUnit"],
+        line_length_factor=1.0,
+        generator_strategies=generator_strategies,
+        scale_link_capital_costs=False,
+    )
+
+    load_mean_final = n.loads_t.p_set.mean().mean()
+    generators_mean_final = n.generators.p_nom.mean()
+
+    logger.info(
+        f"Fetched {len(gdf_islands)} isolated buses into the network. Load attached to a single bus with discrepancies of {(100 * ((load_mean_final - load_mean_origin)/load_mean_origin)):2.1E}% and {(100 * ((generators_mean_final - generators_mean_origin)/generators_mean_origin)):2.1E}% for load and generation capacity, respectively"
+    )
+
+    return clustering.network, busmap
+
+
 def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
     """
     Find isolated nodes in the network and merge those of them which have load
@@ -755,6 +886,8 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
         Original network
     threshold : float
         Load power used as a threshold to merge isolated nodes
+    aggregation_strategies: dictionary
+        Functions to be applied to calculate parameters of the aggregated grid
 
     Returns
     -------
@@ -775,19 +908,19 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
 
     # isolated buses with load below than a specified threshold should be merged
     i_load_islands = n.loads_t.p_set.columns.intersection(i_islands)
-    i_suffic_load = i_load_islands[
+    i_islands_merge = i_load_islands[
         n.loads_t.p_set[i_load_islands].mean(axis=0) <= threshold
     ]
 
     # all the nodes to be merged should be mapped into a single node
     map_isolated_node_by_country = (
         n.buses.assign(bus_id=n.buses.index)
-        .loc[i_suffic_load]
+        .loc[i_islands_merge]
         .groupby("country")["bus_id"]
         .first()
         .to_dict()
     )
-    isolated_buses_mapping = n.buses.loc[i_suffic_load, "country"].replace(
+    isolated_buses_mapping = n.buses.loc[i_islands_merge, "country"].replace(
         map_isolated_node_by_country
     )
     busmap = (
@@ -821,7 +954,7 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
     generators_mean_final = n.generators.p_nom.mean()
 
     logger.info(
-        f"Merged {len(i_suffic_load)} buses. Load attached to a single bus with discrepancies of {(100 * ((load_mean_final - load_mean_origin)/load_mean_origin)):2.1E}% and {(100 * ((generators_mean_final - generators_mean_origin)/generators_mean_origin)):2.1E}% for load and generation capacity, respectively"
+        f"Merged {len(i_islands_merge)} buses. Load attached to a single bus with discrepancies of {(100 * ((load_mean_final - load_mean_origin)/load_mean_origin)):2.1E}% and {(100 * ((generators_mean_final - generators_mean_origin)/generators_mean_origin)):2.1E}% for load and generation capacity, respectively"
     )
 
     return clustering.network, busmap
@@ -976,6 +1109,7 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
         0.0, cluster_config.get("p_threshold_drop_isolated", 0.0)
     )
     p_threshold_merge_isolated = cluster_config.get("p_threshold_merge_isolated", False)
+    s_threshold_fetch_isolated = cluster_config.get("s_threshold_fetch_isolated", False)
 
     n = drop_isolated_nodes(n, threshold=p_threshold_drop_isolated)
     if p_threshold_merge_isolated:
@@ -986,6 +1120,14 @@ def merge_isolated_nodes(n, threshold, aggregation_strategies=dict()):
         )
         busmaps.append(merged_nodes_map)
 
+    if s_threshold_fetch_isolated:
+        n, fetched_nodes_map = merge_into_network(
+            n,
+            threshold=s_threshold_fetch_isolated,
+            aggregation_strategies=aggregation_strategies,
+        )
+        busmaps.append(fetched_nodes_map)
+
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
     n.export_to_netcdf(snakemake.output.network)