From 99299564d62152d226b7f4e0e94c87a0fa8e3fd6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6rsch?= Date: Tue, 30 Jan 2018 23:11:16 +0100 Subject: [PATCH] Split simplify_network from cluster_network --- Snakefile | 39 +++++-- scripts/cluster_network.py | 140 +++++++++-------------- scripts/simplify_network.py | 218 ++++++++++++++++++++++++++++++++++++ 3 files changed, 296 insertions(+), 101 deletions(-) create mode 100644 scripts/simplify_network.py diff --git a/Snakefile b/Snakefile index c658868a..7e22d980 100644 --- a/Snakefile +++ b/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', - regions_onshore="resources/regions_onshore_{network}_{clusters}.geojson", - regions_offshore="resources/regions_offshore_{network}_{clusters}.geojson" - benchmark: "benchmarks/cluster_network/{network}_{clusters}" + network='networks/{network}_s{simpl}.nc', + regions_onshore="resources/regions_onshore_{network}_s{simpl}.geojson", + regions_offshore="resources/regions_offshore_{network}_s{simpl}.geojson" + 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' - output: 'networks/elec_{clusters}_lv{lv}_{opts}.nc' + input: 'networks/{network}_s{simpl}_{clusters}.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" - output: "results/networks/{clusters}_lv{lv}_{opts}.nc" + input: "networks/{network}_s{simpl}_{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", - python="logs/{clusters}_lv{lv}_{opts}_python.log" - benchmark: "benchmarks/solve_network/{clusters}_lv{lv}_{opts}" + gurobi="logs/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_gurobi.log", + python="logs/{network}_s{simpl}_{clusters}_lv{lv}_{opts}_python.log" + 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" diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py index bd621b26..13ee0a5d 100644 --- a/scripts/cluster_network.py +++ b/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): - ## 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): +def weighting_for_country(n, 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): - n_clusters = distribute_clusters_exactly(n_clusters) +def busmap_for_n_clusters(n, 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): - busmap = busmap_for_n_clusters(n_clusters) +def plot_busmap_for_n_clusters(n, n_clusters=50): + 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') - geom_c = regions.geometry.groupby(clustering.busmap).apply(shapely.ops.cascaded_union) + regions = gpd.read_file(getattr(input, which)).set_index('name') + 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 - import yaml - snakemake = Dict() - with open('../config.yaml') as f: - snakemake.config = yaml.load(f) - snakemake.wildcards = Dict(clusters='37') - snakemake.input = Dict(network='../networks/elec.nc', - regions_onshore='../resources/regions_onshore.geojson', - regions_offshore='../resources/regions_offshore.geojson') - snakemake.output = Dict(network='../networks/elec_{clusters}.nc'.format(**snakemake.wildcards), - regions_onshore='../resources/regions_onshore_{clusters}.geojson'.format(**snakemake.wildcards), - regions_offshore='../resources/regions_offshore_{clusters}.geojson'.format(**snakemake.wildcards)) + from vresutils.snakemake import MockSnakemake, Dict + snakemake = MockSnakemake( + path='..', + wildcards=Dict(network='elec', simpl='', clusters='45'), + input=Dict( + 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=Dict( + 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' + ) + ) + + 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,)) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py new file mode 100644 index 00000000..a8502be7 --- /dev/null +++ b/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)