Split simplify_network from cluster_network

Snakefile

@@ -4,6 +4,7 @@ localrules: all, prepare_links_p_nom, base_network, add_electricity, add_sectors
 
 wildcard_constraints:
     lv="[0-9\.]+",
+    simpl="[a-zA-Z0-9]*",
     clusters="[0-9]+",
     sectors="[+a-zA-Z0-9]+",
     opts="[-+a-zA-Z0-9]+"
@@ -72,16 +73,30 @@ rule add_electricity:
     resources: mem_mb=1000
     script: "scripts/add_electricity.py"
 
-rule cluster_network:
+rule simplify_network:
     input:
         network='networks/{network}.nc',
         regions_onshore="resources/regions_onshore.geojson",
         regions_offshore="resources/regions_offshore.geojson"
     output:
-        network='networks/{network}_{clusters}.nc',
+        network='networks/{network}_s{simpl}.nc',
-        regions_onshore="resources/regions_onshore_{network}_{clusters}.geojson",
+        regions_onshore="resources/regions_onshore_{network}_s{simpl}.geojson",
-        regions_offshore="resources/regions_offshore_{network}_{clusters}.geojson"
+        regions_offshore="resources/regions_offshore_{network}_s{simpl}.geojson"
-    benchmark: "benchmarks/cluster_network/{network}_{clusters}"
+    benchmark: "benchmarks/simplify_network/{network}_s{simpl}"
+    threads: 1
+    resources: mem_mb=1000
+    script: "scripts/simplify_network.py"
+
+rule cluster_network:
+    input:
+        network='networks/{network}_s{simpl}.nc',
+        regions_onshore="resources/regions_onshore_{network}_s{simpl}.geojson",
+        regions_offshore="resources/regions_offshore_{network}_s{simpl}.geojson"
+    output:
+        network='networks/{network}_s{simpl}_{clusters}.nc',
+        regions_onshore="resources/regions_onshore_{network}_s{simpl}_{clusters}.geojson",
+        regions_offshore="resources/regions_offshore_{network}_s{simpl}_{clusters}.geojson"
+    benchmark: "benchmarks/cluster_network/{network}_s{simpl}_{clusters}"
     threads: 1
     resources: mem_mb=1000
     script: "scripts/cluster_network.py"
@@ -97,20 +112,20 @@ rule add_sectors:
     script: "scripts/add_sectors.py"
 
 rule prepare_network:
-    input: 'networks/elec_{clusters}.nc'
+    input: 'networks/{network}_s{simpl}_{clusters}.nc'
-    output: 'networks/elec_{clusters}_lv{lv}_{opts}.nc'
+    output: 'networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc'
     threads: 1
     resources: mem_mb=1000
     script: "scripts/prepare_network.py"
 
 rule solve_network:
-    input: "networks/elec_{clusters}_lv{lv}_{opts}.nc"
+    input: "networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc"
-    output: "results/networks/{clusters}_lv{lv}_{opts}.nc"
+    output: "results/networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc"
     shadow: "shallow"
     log:
-        gurobi="logs/{clusters}_lv{lv}_{opts}_gurobi.log",
+        gurobi="logs/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_gurobi.log",
-        python="logs/{clusters}_lv{lv}_{opts}_python.log"
+        python="logs/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_python.log"
-    benchmark: "benchmarks/solve_network/{clusters}_lv{lv}_{opts}"
+    benchmark: "benchmarks/solve_network/{network}_s{simpl}_{clusters}_lv{lv}_{opts}"
     threads: 4
     resources: mem_mb=lambda w: 100000 * int(w.clusters) // 362
     script: "scripts/solve_network.py"

scripts/cluster_network.py

@@ -3,69 +3,30 @@
 import pandas as pd
 idx = pd.IndexSlice
 
+import logging
+logger = logging.getLogger(__name__)
+
 import os
 import numpy as np
 import scipy as sp
+from scipy.sparse.csgraph import connected_components
 import xarray as xr
 import geopandas as gpd
 import shapely
+import networkx as nx
 
+from six import iteritems
 from six.moves import reduce
 
 import pypsa
+from pypsa.io import import_components_from_dataframe, import_series_from_dataframe
 from pypsa.networkclustering import (busmap_by_stubs, busmap_by_kmeans,
-                                     _make_consense, get_clustering_from_busmap)
+                                     _make_consense, get_clustering_from_busmap,
+                                     aggregategenerators, aggregateoneport)
 def normed(x):
     return (x/x.sum()).fillna(0.)
 
-def simplify_network_to_380(n):
+def weighting_for_country(n, x):
-    ## All goes to v_nom == 380
-
-    n.buses['v_nom'] = 380.
-
-    linetype_380, = n.lines.loc[n.lines.v_nom == 380., 'type'].unique()
-    lines_v_nom_b = n.lines.v_nom != 380.
-    n.lines.loc[lines_v_nom_b, 'num_parallel'] *= (n.lines.loc[lines_v_nom_b, 'v_nom'] / 380.)**2
-    n.lines.loc[lines_v_nom_b, 'v_nom'] = 380.
-    n.lines.loc[lines_v_nom_b, 'type'] = linetype_380
-
-    # Replace transformers by lines
-    trafo_map = pd.Series(n.transformers.bus1.values, index=n.transformers.bus0.values)
-    trafo_map = trafo_map[~trafo_map.index.duplicated(keep='first')]
-    several_trafo_b = trafo_map.isin(trafo_map.index)
-    trafo_map.loc[several_trafo_b] = trafo_map.loc[several_trafo_b].map(trafo_map)
-    missing_buses_i = n.buses.index.difference(trafo_map.index)
-    trafo_map = trafo_map.append(pd.Series(missing_buses_i, missing_buses_i))
-
-    for c in n.one_port_components|n.branch_components:
-        df = n.df(c)
-        for col in df.columns:
-            if col.startswith('bus'):
-                df[col] = df[col].map(trafo_map)
-
-    n.mremove("Transformer", n.transformers.index)
-    n.mremove("Bus", n.buses.index.difference(trafo_map))
-
-    return n, trafo_map
-
-def remove_stubs(n):
-    n.determine_network_topology()
-
-    busmap = busmap_by_stubs(n, ['carrier', 'country'])
-
-    n.buses.loc[busmap.index, ['x','y']] = n.buses.loc[busmap, ['x','y']].values
-
-    clustering = get_clustering_from_busmap(
-        n, busmap,
-        bus_strategies=dict(country=_make_consense("Bus", "country")),
-        line_length_factor=snakemake.config['lines']['length_factor'],
-        aggregate_generators_weighted=True,
-        aggregate_one_ports=["Load", "StorageUnit"]
-    )
-
-    return clustering.network, busmap
-
-def weighting_for_country(x):
     conv_carriers = {'OCGT', 'PHS', 'hydro'}
     gen = (n
            .generators.loc[n.generators.carrier.isin(conv_carriers)]
@@ -84,8 +45,6 @@ def weighting_for_country(x):
     w= g + l
     return (w * (100. / w.max())).astype(int)
 
-    return weighting_for_country
-
 
 ## Plot weighting for Germany
 
@@ -98,7 +57,7 @@ def plot_weighting(n, country):
 
 # # Determining the number of clusters per country
 
-def distribute_clusters(n_clusters):
+def distribute_clusters(n, n_clusters):
     load = n.loads_t.p_set.mean().groupby(n.loads.bus).sum()
     loadc = load.groupby([n.buses.country, n.buses.sub_network]).sum()
     n_cluster_per_country = n_clusters * normed(loadc)
@@ -117,44 +76,41 @@ def distribute_clusters(n_clusters):
 
     return n_cluster_per_country.astype(int)
 
-def distribute_clusters_exactly(n_clusters):
+def distribute_clusters_exactly(n, n_clusters):
     for d in [0, 1, -1, 2, -2]:
-        n_cluster_per_country = distribute_clusters(n_clusters + d)
+        n_cluster_per_country = distribute_clusters(n, n_clusters + d)
         if n_cluster_per_country.sum() == n_clusters:
             return n_cluster_per_country
     else:
-        return distribute_clusters(n_clusters)
+        return distribute_clusters(n, n_clusters)
 
-def busmap_for_n_clusters(n_clusters):
+def busmap_for_n_clusters(n, n_clusters):
-    n_clusters = distribute_clusters_exactly(n_clusters)
+    n.determine_network_topology()
+
+    n_clusters = distribute_clusters_exactly(n, n_clusters)
     def busmap_for_country(x):
         prefix = x.name[0] + x.name[1] + ' '
         if len(x) == 1:
             return pd.Series(prefix + '0', index=x.index)
-        weight = weighting_for_country(x)
+        weight = weighting_for_country(n, x)
         return prefix + busmap_by_kmeans(n, weight, n_clusters[x.name], buses_i=x.index)
     return n.buses.groupby(['country', 'sub_network'], group_keys=False).apply(busmap_for_country)
 
-def plot_busmap_for_n_clusters(n_clusters=50):
+def plot_busmap_for_n_clusters(n, n_clusters=50):
-    busmap = busmap_for_n_clusters(n_clusters)
+    busmap = busmap_for_n_clusters(n, n_clusters)
     cs = busmap.unique()
     cr = sns.color_palette("hls", len(cs))
     n.plot(bus_colors=busmap.map(dict(zip(cs, cr))))
     del cs, cr
 
-def clustering_for_n_clusters(n_clusters):
+def clustering_for_n_clusters(n, n_clusters):
     clustering = get_clustering_from_busmap(
-        n, busmap_for_n_clusters(n_clusters),
+        n, busmap_for_n_clusters(n, n_clusters),
         bus_strategies=dict(country=_make_consense("Bus", "country")),
         aggregate_generators_weighted=True,
         aggregate_one_ports=["Load", "StorageUnit"]
     )
 
-    # set n-1 security margin to 0.5 for 37 clusters and to 0.7 from 200 clusters
-    # (there was already one of 0.7 in-place)
-    s_max_pu = np.clip(0.5 + 0.2 * (n_clusters - 37) / (200 - 37), 0.5, 0.7)
-    clustering.network.lines['s_max_pu'] = s_max_pu
-
     return clustering
 
 def save_to_geojson(s, fn):
@@ -162,40 +118,46 @@ def save_to_geojson(s, fn):
         os.unlink(fn)
     s.reset_index().to_file(fn, driver='GeoJSON')
 
-def cluster_regions(busmaps):
+def cluster_regions(busmaps, input=None, output=None):
+    if input is None: input = snakemake.input
+    if output is None: output = snakemake.output
+
     busmap = reduce(lambda x, y: x.map(y), busmaps[1:], busmaps[0])
 
     for which in ('regions_onshore', 'regions_offshore'):
-        regions = gpd.read_file(getattr(snakemake.input, which)).set_index('name')
+        regions = gpd.read_file(getattr(input, which)).set_index('name')
-        geom_c = regions.geometry.groupby(clustering.busmap).apply(shapely.ops.cascaded_union)
+        geom_c = regions.geometry.groupby(busmap).apply(shapely.ops.cascaded_union)
         regions_c = gpd.GeoDataFrame(dict(geometry=geom_c))
-        save_to_geojson(regions_c, getattr(snakemake.output, which))
+        regions_c.index.name = 'name'
+        save_to_geojson(regions_c, getattr(output, which))
 
 if __name__ == "__main__":
     # Detect running outside of snakemake and mock snakemake for testing
     if 'snakemake' not in globals():
-        from vresutils import Dict
+        from vresutils.snakemake import MockSnakemake, Dict
-        import yaml
+        snakemake = MockSnakemake(
-        snakemake = Dict()
+            path='..',
-        with open('../config.yaml') as f:
+            wildcards=Dict(network='elec', simpl='', clusters='45'),
-            snakemake.config = yaml.load(f)
+            input=Dict(
-        snakemake.wildcards = Dict(clusters='37')
+                network='networks/{network}_s{simpl}.nc',
-        snakemake.input = Dict(network='../networks/elec.nc',
+                regions_onshore='resources/regions_onshore_{network}_s{simpl}.geojson',
-                               regions_onshore='../resources/regions_onshore.geojson',
+                regions_offshore='resources/regions_offshore_{network}_s{simpl}.geojson'
-                               regions_offshore='../resources/regions_offshore.geojson')
+            ),
-        snakemake.output = Dict(network='../networks/elec_{clusters}.nc'.format(**snakemake.wildcards),
+            output=Dict(
-                                regions_onshore='../resources/regions_onshore_{clusters}.geojson'.format(**snakemake.wildcards),
+                network='networks/{network}_s{simpl}_{clusters}.nc',
-                                regions_offshore='../resources/regions_offshore_{clusters}.geojson'.format(**snakemake.wildcards))
+                regions_onshore='resources/regions_onshore_{network}_s{simpl}_{clusters}.geojson',
+                regions_offshore='resources/regions_offshore_{network}_s{simpl}_{clusters}.geojson'
+            )
+        )
+
+    logger = logging.getLogger()
+    logger.setLevel(snakemake.config['logging_level'])
 
     n = pypsa.Network(snakemake.input.network)
 
-    n, trafo_map = simplify_network_to_380(n)
-
-    n, stub_map = remove_stubs(n)
-
     n_clusters = int(snakemake.wildcards.clusters)
-    clustering = clustering_for_n_clusters(n_clusters)
+    clustering = clustering_for_n_clusters(n, n_clusters)
 
     clustering.network.export_to_netcdf(snakemake.output.network)
 
-    cluster_regions((trafo_map, stub_map, clustering.busmap))
+    cluster_regions((clustering.busmap,))

scripts/simplify_network.py

@@ -0,0 +1,218 @@
+# coding: utf-8
+
+import pandas as pd
+idx = pd.IndexSlice
+
+import logging
+logger = logging.getLogger(__name__)
+
+import os
+import re
+import numpy as np
+import scipy as sp
+from scipy.sparse.csgraph import connected_components
+import xarray as xr
+import geopandas as gpd
+import shapely
+import networkx as nx
+
+from six import iteritems
+from six.moves import reduce
+
+import pypsa
+from pypsa.io import import_components_from_dataframe, import_series_from_dataframe
+from pypsa.networkclustering import (busmap_by_stubs, busmap_by_kmeans,
+                                     _make_consense, get_clustering_from_busmap,
+                                     aggregategenerators, aggregateoneport)
+
+from cluster_network import clustering_for_n_clusters, cluster_regions
+
+def simplify_network_to_380(n):
+    ## All goes to v_nom == 380
+    logger.info("Mapping all network lines onto a single 380kV layer")
+
+    n.buses['v_nom'] = 380.
+
+    linetype_380, = n.lines.loc[n.lines.v_nom == 380., 'type'].unique()
+    lines_v_nom_b = n.lines.v_nom != 380.
+    n.lines.loc[lines_v_nom_b, 'num_parallel'] *= (n.lines.loc[lines_v_nom_b, 'v_nom'] / 380.)**2
+    n.lines.loc[lines_v_nom_b, 'v_nom'] = 380.
+    n.lines.loc[lines_v_nom_b, 'type'] = linetype_380
+
+    # Replace transformers by lines
+    trafo_map = pd.Series(n.transformers.bus1.values, index=n.transformers.bus0.values)
+    trafo_map = trafo_map[~trafo_map.index.duplicated(keep='first')]
+    several_trafo_b = trafo_map.isin(trafo_map.index)
+    trafo_map.loc[several_trafo_b] = trafo_map.loc[several_trafo_b].map(trafo_map)
+    missing_buses_i = n.buses.index.difference(trafo_map.index)
+    trafo_map = trafo_map.append(pd.Series(missing_buses_i, missing_buses_i))
+
+    for c in n.one_port_components|n.branch_components:
+        df = n.df(c)
+        for col in df.columns:
+            if col.startswith('bus'):
+                df[col] = df[col].map(trafo_map)
+
+    n.mremove("Transformer", n.transformers.index)
+    n.mremove("Bus", n.buses.index.difference(trafo_map))
+
+    return n, trafo_map
+
+def _aggregate_and_move_components(n, busmap, aggregate_one_ports={"Load", "StorageUnit"}):
+    def replace_components(n, c, df, pnl):
+        n.mremove(c, n.df(c).index)
+
+        import_components_from_dataframe(n, df, c)
+        for attr, df in iteritems(pnl):
+            if not df.empty:
+                import_series_from_dataframe(n, df, c, attr)
+
+    generators, generators_pnl = aggregategenerators(n, busmap)
+    replace_components(n, "Generator", generators, generators_pnl)
+
+    for one_port in aggregate_one_ports:
+        df, pnl = aggregateoneport(n, busmap, component=one_port)
+        replace_components(n, one_port, df, pnl)
+
+    buses_to_del = n.buses.index.difference(busmap)
+    n.mremove("Bus", buses_to_del)
+    for c in n.branch_components:
+        df = n.df(c)
+        n.mremove(c, df.index[df.bus0.isin(buses_to_del) | df.bus1.isin(buses_to_del)])
+
+def simplify_links(n):
+    ## Complex multi-node links are folded into end-points
+    logger.info("Simplifying connected link components")
+
+    # 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)))
+
+    _, labels = connected_components(adjacency_matrix, directed=False)
+    labels = pd.Series(labels, n.buses.index)
+
+    G = n.graph()
+
+    def split_links(nodes):
+        nodes = frozenset(nodes)
+
+        seen = set()
+        supernodes = {m for m in nodes
+                      if len(G.adj[m]) > 2 or (set(G.adj[m]) - nodes)}
+
+        for u in supernodes:
+            for m, ls in iteritems(G.adj[u]):
+                if m not in nodes or m in seen: continue
+
+                buses = [u, m]
+                links = [list(ls)] #[name for name in ls]]
+
+                while m not in (supernodes | seen):
+                    seen.add(m)
+                    for m2, ls in iteritems(G.adj[m]):
+                        if m2 in seen or m2 == u: continue
+                        buses.append(m2)
+                        links.append(list(ls)) # [name for name in ls])
+                        break
+                    else:
+                        # stub
+                        break
+                    m = m2
+                if m != u:
+                    yield pd.Index((u, m)), buses, links
+            seen.add(u)
+
+    busmap = n.buses.index.to_series()
+
+    for lbl in labels.value_counts().loc[lambda s: s > 2].index:
+
+        for b, buses, links 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))
+
+            m = sp.spatial.distance_matrix(n.buses.loc[b, ['x', 'y']],
+                                           n.buses.loc[buses[1:-1], ['x', 'y']])
+            busmap.loc[buses] = b[np.r_[0, m.argmin(axis=0), 1]]
+            all_links = [i for _, i in sum(links, [])]
+
+            s_max_pu = snakemake.config['links']['s_max_pu']
+            name = n.links.loc[all_links, 'length'].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),
+                p_max_pu=s_max_pu,
+                p_min_pu=-s_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))
+
+            n.mremove("Link", all_links)
+
+            static_attrs = n.components["Link"]["attrs"].loc[lambda df: df.static]
+            for attr, default in static_attrs.default.iteritems(): params.setdefault(attr, default)
+            n.links.loc[name] = pd.Series(params)
+
+            # n.add("Link", **params)
+
+    logger.debug("Collecting all components using the busmap")
+
+    _aggregate_and_move_components(n, busmap)
+    return n, busmap
+
+def remove_stubs(n):
+    logger.info("Removing stubs")
+
+    busmap = busmap_by_stubs(n, ['country'])
+    _aggregate_and_move_components(n, busmap)
+
+    return n, busmap
+
+
+if __name__ == "__main__":
+    # Detect running outside of snakemake and mock snakemake for testing
+    if 'snakemake' not in globals():
+        from vresutils.snakemake import MockSnakemake, Dict
+        snakemake = MockSnakemake(
+            path='..',
+            wildcards=Dict(simpl=''),
+            input=Dict(
+                network='networks/elec.nc',
+                regions_onshore='resources/regions_onshore.geojson',
+                regions_offshore='resources/regions_offshore.geojson'
+            ),
+            output=Dict(
+                network='networks/elec_s{simpl}.nc',
+                regions_onshore='resources/regions_onshore_s{simpl}.geojson',
+                regions_offshore='resources/regions_offshore_s{simpl}.geojson'
+            )
+        )
+
+    logger = logging.getLogger()
+    logger.setLevel(snakemake.config['logging_level'])
+
+    n = pypsa.Network(snakemake.input.network)
+
+    n, trafo_map = simplify_network_to_380(n)
+
+    n, simplify_links_map = simplify_links(n)
+
+    n, stub_map = remove_stubs(n)
+
+    busmaps = [trafo_map, simplify_links_map, stub_map]
+
+    if snakemake.wildcards.simpl:
+        n_clusters = int(snakemake.wildcards.simpl)
+        clustering = clustering_for_n_clusters(n_clusters)
+
+        n = clustering.network
+        busmaps.append(clustering.busmap)
+
+    n.export_to_netcdf(snakemake.output.network)
+
+    cluster_regions(busmaps, snakemake.input, snakemake.output)