Revise network simplification to account for DC lines

doc/release_notes.rst

@@ -14,6 +14,8 @@ E.g. if a new rule becomes available describe how to use it `snakemake -j1 run_t
 
 **New Features and major Changes**
 
+* Improve network simplification routine to account for representation HVDC as Line component `PR #743 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/743>`__
+
 PyPSA-Earth 0.2.2
 =================
 

scripts/add_electricity.py

@@ -253,40 +253,52 @@ def attach_load(n, demand_profiles):
     n.madd("Load", demand_df.columns, bus=demand_df.columns, p_set=demand_df)
 
 
-def update_transmission_costs(n, costs, length_factor=1.0, simple_hvdc_costs=False):
-    n.lines["capital_cost"] = (
-        n.lines["length"] * length_factor * costs.at["HVAC overhead", "capital_cost"]
-    )
-
-    if n.links.empty:
+def attach_dc_costs(lines_or_links, costs, length_factor=1.0, simple_hvdc_costs=False):
+    if lines_or_links.empty:
         return
 
-    dc_b = n.links.carrier == "DC"
-    # If there are no "DC" links, then the 'underwater_fraction' column
-    # may be missing. Therefore we have to return here.
-    # TODO: Require fix
-    if n.links.loc[n.links.carrier == "DC"].empty:
+    if lines_or_links.loc[lines_or_links.carrier == "DC"].empty:
         return
 
+    dc_b = lines_or_links.carrier == "DC"
     if simple_hvdc_costs:
         costs = (
-            n.links.loc[dc_b, "length"]
+            lines_or_links.loc[dc_b, "length"]
             * length_factor
             * costs.at["HVDC overhead", "capital_cost"]
         )
     else:
         costs = (
-            n.links.loc[dc_b, "length"]
+            lines_or_links.loc[dc_b, "length"]
             * length_factor
             * (
-                (1.0 - n.links.loc[dc_b, "underwater_fraction"])
+                (1.0 - lines_or_links.loc[dc_b, "underwater_fraction"])
                 * costs.at["HVDC overhead", "capital_cost"]
-                + n.links.loc[dc_b, "underwater_fraction"]
+                + lines_or_links.loc[dc_b, "underwater_fraction"]
                 * costs.at["HVDC submarine", "capital_cost"]
             )
             + costs.at["HVDC inverter pair", "capital_cost"]
         )
-    n.links.loc[dc_b, "capital_cost"] = costs
+    lines_or_links.loc[dc_b, "capital_cost"] = costs
+
+
+def update_transmission_costs(n, costs, length_factor=1.0, simple_hvdc_costs=False):
+    n.lines["capital_cost"] = (
+        n.lines["length"] * length_factor * costs.at["HVAC overhead", "capital_cost"]
+    )
+
+    attach_dc_costs(
+        lines_or_links=n.links,
+        costs=costs,
+        length_factor=length_factor,
+        simple_hvdc_costs=simple_hvdc_costs,
+    )
+    attach_dc_costs(
+        lines_or_links=n.lines,
+        costs=costs,
+        length_factor=length_factor,
+        simple_hvdc_costs=simple_hvdc_costs,
+    )
 
 
 def attach_wind_and_solar(

scripts/augmented_line_connections.py

@@ -37,7 +37,7 @@
 import pypsa
 from _helpers import configure_logging
 from add_electricity import load_costs
-from base_network import _set_links_underwater_fraction
+from base_network import _set_dc_underwater_fraction
 from networkx.algorithms import complement
 from networkx.algorithms.connectivity.edge_augmentation import k_edge_augmentation
 from pypsa.geo import haversine_pts
@@ -86,6 +86,7 @@ def haversine(p):
 
     # TODO: Currently only AC lines are read in and meshed. One need to combine
     # AC & DC lines and then move on.
+    # rather there is a need to filter-out DC lines
     network_lines = n.lines
     sel = network_lines.s_nom > 100  # TODO: Check, should be all selected or filtered?
     attrs = ["bus0", "bus1", "length"]
@@ -168,7 +169,8 @@ def haversine(p):
             lifetime=costs.at["HVAC overhead", "lifetime"],
         )
 
-        _set_links_underwater_fraction(snakemake.input.regions_offshore, n)
+        _set_dc_underwater_fraction(n.links, snakemake.input.regions_offshore)
+        _set_dc_underwater_fraction(n.lines, snakemake.input.regions_offshore)
 
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
     n.export_to_netcdf(snakemake.output.network)

scripts/base_network.py

@@ -142,10 +142,7 @@ def _load_buses_from_osm(fp_buses, config):
     return buses
 
 
-def _set_links_underwater_fraction(fp_offshore_shapes, n):
-    if n.links.empty:
-        return
-
+def add_underwater_links(n, fp_offshore_shapes):
     if not hasattr(n.links, "geometry"):
         n.links["underwater_fraction"] = 0.0
     else:
@@ -159,6 +156,38 @@ def _set_links_underwater_fraction(fp_offshore_shapes, n):
             )
 
 
+def _set_dc_underwater_fraction(lines_or_links, fp_offshore_shapes):
+    # HVDC part always has some links as converters
+    # excluding probably purely DC networks which are currently somewhat exotic
+    if lines_or_links.empty:
+        return
+
+    if lines_or_links.loc[lines_or_links.carrier == "DC"].empty:
+        # Add "underwater_fraction" both to lines and links
+        lines_or_links["underwater_fraction"] = 0.0
+        return
+
+    if not hasattr(lines_or_links, "geometry"):
+        lines_or_links["underwater_fraction"] = 0.0
+    else:
+        offshore_shape = gpd.read_file(fp_offshore_shapes).unary_union
+        if offshore_shape is None or offshore_shape.is_empty:
+            lines_or_links["underwater_fraction"] = 0.0
+        else:
+            branches = gpd.GeoSeries(
+                lines_or_links.geometry.dropna().map(shapely.wkt.loads)
+            )
+            # fix to avoid NaN for links during augmentation
+            if branches.empty:
+                lines_or_links["underwater_fraction"] = 0
+            else:
+                lines_or_links["underwater_fraction"] = (
+                    # TODO Check assumption that all underwater lines are DC
+                    branches.intersection(offshore_shape).length
+                    / branches.length
+                )
+
+
 def _load_lines_from_osm(fp_osm_lines, config, buses):
     lines = (
         read_csv_nafix(
@@ -322,15 +351,12 @@ def _set_lines_s_nom_from_linetypes(n):
     n.lines["s_nom"] = (
         np.sqrt(3)
         * n.lines["type"].map(n.line_types.i_nom)
-        * n.lines["v_nom"]
-        * n.lines.num_parallel
+        * n.lines.eval("v_nom * num_parallel")
     )
     # Re-define s_nom for DC lines
-    n.lines.loc[n.lines["carrier"] == "DC", "s_nom"] = (
-        n.lines["type"].map(n.line_types.i_nom)
-        * n.lines["v_nom"]
-        * n.lines.num_parallel
-    )
+    n.lines.loc[n.lines["carrier"] == "DC", "s_nom"] = n.lines["type"].map(
+        n.line_types.i_nom
+    ) * n.lines.eval("v_nom * num_parallel")
 
 
 def _remove_dangling_branches(branches, buses):
@@ -499,7 +525,8 @@ def base_network(inputs, config):
 
     _set_countries_and_substations(inputs, config, n)
 
-    _set_links_underwater_fraction(inputs.offshore_shapes, n)
+    _set_dc_underwater_fraction(n.lines, inputs.offshore_shapes)
+    _set_dc_underwater_fraction(n.links, inputs.offshore_shapes)
 
     return n
 

scripts/cluster_network.py

@@ -614,6 +614,8 @@ def clustering_for_n_clusters(
             .dropna(),
             fill_value=0,
         )
+    if not n.lines.loc[n.lines.carrier == "DC"].empty:
+        clustering.network.lines["underwater_fraction"] = 0
 
     return clustering
 

scripts/simplify_network.py

@@ -129,6 +129,12 @@ def simplify_network_to_base_voltage(n, linetype, base_voltage):
     # Note: s_nom is set in base_network
     n.lines["num_parallel"] = n.lines.eval("s_nom / (sqrt(3) * v_nom * i_nom)")
 
+    # Re-define s_nom for DC lines
+    is_dc_carrier = n.lines["carrier"] == "DC"
+    n.lines.loc[is_dc_carrier, "num_parallel"] = n.lines.loc[is_dc_carrier].eval(
+        "s_nom / (v_nom * i_nom)"
+    )
+
     # Replace transformers by lines
     trafo_map = pd.Series(n.transformers.bus1.values, n.transformers.bus0.values)
     trafo_map = trafo_map[~trafo_map.index.duplicated(keep="first")]
@@ -155,15 +161,22 @@ def _prepare_connection_costs_per_link(n, costs, config):
 
     connection_costs_per_link = {}
 
+    if not n.links.loc[n.links.carrier == "DC"].empty:
+        dc_lengths = n.links.length
+        unterwater_fractions = n.links.underwater_fraction
+    elif not n.lines.loc[n.lines.carrier == "DC"].empty:
+        dc_lengths = n.lines.length
+        unterwater_fractions = n.lines.underwater_fraction
+
     for tech in config["renewable"]:
         if tech.startswith("offwind"):
             connection_costs_per_link[tech] = (
-                n.links.length
+                dc_lengths
                 * config["lines"]["length_factor"]
                 * (
-                    n.links.underwater_fraction
+                    unterwater_fractions
                     * costs.at[tech + "-connection-submarine", "capital_cost"]
-                    + (1.0 - n.links.underwater_fraction)
+                    + (1.0 - unterwater_fractions)
                     * costs.at[tech + "-connection-underground", "capital_cost"]
                 )
             )
@@ -287,10 +300,15 @@ def simplify_links(n, costs, config, output, aggregation_strategies=dict()):
             pass
         return n, n.buses.index.to_series()
 
+    dc_as_links = not (n.lines.carrier == "DC").any()
+
     # Determine connected link components, ignore all links but DC
     adjacency_matrix = n.adjacency_matrix(
-        branch_components=["Link"],
-        weights=dict(Link=(n.links.carrier == "DC").astype(float)),
+        branch_components=["Link", "Line"],
+        weights=dict(
+            Link=(n.links.carrier == "DC").astype(float),
+            Line=(n.lines.carrier == "DC").astype(float),
+        ),
     )
 
     _, labels = connected_components(adjacency_matrix, directed=False)
@@ -300,6 +318,15 @@ def simplify_links(n, costs, config, output, aggregation_strategies=dict()):
 
     # Split DC part by supernodes
     def split_links(nodes):
+        # only DC nodes are of interest for further supernodes treatment
+        nodes_links = n.links["bus0"].to_list() + n.links["bus1"].to_list()
+        nodes_dc_lines = [
+            d
+            for d in nodes
+            if n.lines.loc[(n.lines.bus0 == d) | (n.lines.bus1 == d)].dc.any()
+        ]
+        nodes = nodes_links + nodes_dc_lines
+
         nodes = frozenset(nodes)
 
         seen = set()
@@ -319,7 +346,6 @@ def split_links(nodes):
                     seen.add(m)
                     for m2, ls2 in G.adj[m].items():
                         # there may be AC lines which connect ends of DC chains
-                        # TODO remove after debug
                         if m2 in seen or m2 == u or contains_ac(ls2):
                             continue
                         buses.append(m2)
@@ -341,12 +367,12 @@ def split_links(nodes):
     )
 
     for lbl in labels.value_counts().loc[lambda s: s > 2].index:
-        for b, buses, links in split_links(labels.index[labels == lbl]):
+        for b, buses, dc_edges in split_links(labels.index[labels == lbl]):
             if len(buses) <= 2:
                 continue
 
             logger.debug("nodes = {}".format(labels.index[labels == lbl]))
-            logger.debug("b = {}\nbuses = {}\nlinks = {}".format(b, buses, links))
+            logger.debug("b = {}\nbuses = {}\nlinks = {}".format(b, buses, dc_edges))
 
             m = sp.spatial.distance_matrix(
                 n.buses.loc[b, ["x", "y"]], n.buses.loc[buses[1:-1], ["x", "y"]]
@@ -356,37 +382,62 @@ def split_links(nodes):
                 n, busmap, costs, config, connection_costs_per_link, buses
             )
 
-            all_links = [i for _, i in sum(links, [])]
+            # TODO revise a variable name for `dc_edges` variable
+            # `dc_edges` is a list containing dc-relevant graph elements like [('Line', '712308316-1_0')]
+            all_dc_branches = [i for _, i in sum(dc_edges, [])]
+            all_links = list(set(n.links.index).intersection(all_dc_branches))
+            all_dc_lines = list(set(n.lines.index).intersection(all_dc_branches))
+
+            # p_max_pu = config["links"].get("p_max_pu", 1.0)
+            all_dc_lengths = pd.concat(
+                [n.links.loc[all_links, "length"], n.lines.loc[all_dc_lines, "length"]]
+            )
+            name = all_dc_lengths.idxmax() + "+{}".format(len(all_dc_branches) - 1)
+
+            # HVDC part is represented as "Link" component
+            if dc_as_links:
+                p_max_pu = config["links"].get("p_max_pu", 1.0)
+                lengths = n.links.loc[all_links, "length"]
+                i_links = [i for _, i in dc_edges if _ == "Link"]
+                length = sum(n.links.loc[i_links, "length"].mean() for l in dc_edges)
+                p_nom = min(n.links.loc[i_links, "p_nom"].sum() for l in dc_edges)
+                underwater_fraction = (
+                    lengths * n.links.loc[all_links, "underwater_fraction"]
+                ).sum() / lengths.sum()
+            # HVDC part is represented as "Line" component
+            else:
+                p_max_pu = config["lines"].get("p_max_pu", 1.0)
+                lengths = n.lines.loc[all_dc_lines, "length"]
+                length = lengths.sum() / len(lengths) if len(lengths) > 0 else 0
+                p_nom = n.lines.loc[all_dc_lines, "s_nom"].min()
+                underwater_fraction = (
+                    (lengths * n.lines.loc[all_dc_lines, "underwater_fraction"]).sum()
+                    / lengths.sum()
+                    if len(lengths) > 0
+                    else 0
+                )
 
-            p_max_pu = config["links"].get("p_max_pu", 1.0)
-            lengths = n.links.loc[all_links, "length"]
-            name = lengths.idxmax() + "+{}".format(len(links) - 1)
             params = dict(
                 carrier="DC",
                 bus0=b[0],
                 bus1=b[1],
-                length=sum(
-                    n.links.loc[[i for _, i in l], "length"].mean() for l in links
-                ),
-                p_nom=min(n.links.loc[[i for _, i in l], "p_nom"].sum() for l in links),
-                underwater_fraction=sum(
-                    lengths
-                    / lengths.sum()
-                    * n.links.loc[all_links, "underwater_fraction"]
-                ),
+                length=length,
+                p_nom=p_nom,
+                underwater_fraction=underwater_fraction,
                 p_max_pu=p_max_pu,
                 p_min_pu=-p_max_pu,
                 underground=False,
                 under_construction=False,
             )
 
             logger.info(
-                "Joining the links {} connecting the buses {} to simple link {}".format(
-                    ", ".join(all_links), ", ".join(buses), name
+                "Joining the links and DC lines {} connecting the buses {} to simple link {}".format(
+                    ", ".join(all_dc_branches), ", ".join(buses), name
                 )
             )
 
             n.mremove("Link", all_links)
+            n.mremove("Line", all_dc_lines)
 
             static_attrs = n.components["Link"]["attrs"].loc[lambda df: df.static]
             for attr, default in static_attrs.default.items():