{simplify,cluster}_network: Fix two-step clustering and introduce m suffix

- a network name like elec_s1000_ now clusters the network down to 1000 buses in
  the simplify step.
- an 'm' after the number of clusters as in elec_s1000_181m_ skips the
  aggregation of renewable generators and just moves them to the new clustered
  bus in the second clustering step.
- to distribute the number of clusters to countries a small quadratic
  optimization is now performed to minimize \sum_i (n_i - L_i/\sum_j L_j N)**2,
  where n_i >= 1.

Snakefile

@@ -1,11 +1,11 @@
 configfile: "config.yaml"
 
-localrules: all, prepare_links_p_nom, base_network, build_powerplants, add_electricity, add_sectors, prepare_network, extract_summaries, plot_network, scenario_comparions
+localrules: all, prepare_links_p_nom, base_network, build_renewable_potentials, build_powerplants, add_electricity, add_sectors, prepare_network, extract_summaries, plot_network, scenario_comparions
 
 wildcard_constraints:
     lv="[0-9\.]+",
     simpl="[a-zA-Z0-9]*",
-    clusters="[0-9]+",
+    clusters="[0-9]+m?",
     sectors="[+a-zA-Z0-9]+",
     opts="[-+a-zA-Z0-9]+"
 
@@ -97,7 +97,7 @@ rule simplify_network:
         regions_offshore="resources/regions_offshore_{network}_s{simpl}.geojson"
     benchmark: "benchmarks/simplify_network/{network}_s{simpl}"
     threads: 1
-    resources: mem_mb=3000
+    resources: mem_mb=4000
     script: "scripts/simplify_network.py"
 
 rule cluster_network:
@@ -133,7 +133,14 @@ rule prepare_network:
     script: "scripts/prepare_network.py"
 
 def partition(w):
-    return 'vres' if int(w.clusters) >= 256 else 'x-men'
+    n_clusters = int(w.clusters[:-1] if w.clusters.endswith('m') else w.clusters)
+    return 'vres' if n_clusters >= 256 else 'x-men'
+
+def memory(w):
+    if w.clusters.endswith('m'):
+        return 61000
+    else:
+        return 4890+310 * int(w.clusters)
 
 rule solve_network:
     input: "networks/{network}_s{simpl}_{clusters}_lv{lv}_{opts}.nc"
@@ -147,7 +154,7 @@ rule solve_network:
     benchmark: "benchmarks/solve_network/{network}_s{simpl}_{clusters}_lv{lv}_{opts}"
     threads: 4
     resources:
-        mem_mb=lambda w: 4890+310 * int(w.clusters), # without 5000 too small for 256
+        mem_mb=memory,
         x_men=lambda w: 1 if partition(w) == 'x-men' else 0,
         vres=lambda w: 1 if partition(w) == 'vres' else 0
     script: "scripts/solve_network.py"

scripts/cluster_network.py

@@ -18,6 +18,8 @@ import networkx as nx
 from six import iteritems
 from six.moves import reduce
 
+import pyomo.environ as po
+
 import pypsa
 from pypsa.io import import_components_from_dataframe, import_series_from_dataframe
 from pypsa.networkclustering import (busmap_by_stubs, busmap_by_kmeans,
@@ -84,10 +86,37 @@ def distribute_clusters_exactly(n, n_clusters):
     else:
         return distribute_clusters(n, n_clusters)
 
+def distribute_clusters_optim(n, n_clusters, solver_name='gurobi'):
+    L = (n.loads_t.p_set.mean()
+         .groupby(n.loads.bus).sum()
+         .groupby([n.buses.country, n.buses.sub_network]).sum()
+         .pipe(normed))
+
+    m = po.ConcreteModel()
+    m.n = po.Var(list(L.index), bounds=(1, None), domain=po.Integers)
+    m.tot = po.Constraint(expr=(po.summation(m.n) == n_clusters))
+    m.objective = po.Objective(expr=po.sum((m.n[i] - L.loc[i]*n_clusters)**2
+                                           for i in L.index),
+                               sense=po.minimize)
+
+    opt = po.SolverFactory(solver_name)
+    if isinstance(opt, pypsa.opf.PersistentSolver):
+        opt.set_instance(m)
+    results = opt.solve(m)
+    assert results['Solver'][0]['Status'].key == 'ok', "Solver returned non-optimally: {}".format(results)
+
+    return pd.Series(m.n.get_values(), index=L.index).astype(int)
+
 def busmap_for_n_clusters(n, n_clusters):
     n.determine_network_topology()
 
-    n_clusters = distribute_clusters_exactly(n, n_clusters)
+    if 'snakemake' in globals():
+        solver_name = snakemake.config['solving']['solver']['name']
+    else:
+        solver_name = "gurobi"
+
+    n_clusters = distribute_clusters_optim(n, n_clusters, solver_name=solver_name)
+
     def busmap_for_country(x):
         prefix = x.name[0] + x.name[1] + ' '
         if len(x) == 1:
@@ -103,11 +132,15 @@ def plot_busmap_for_n_clusters(n, n_clusters=50):
     n.plot(bus_colors=busmap.map(dict(zip(cs, cr))))
     del cs, cr
 
-def clustering_for_n_clusters(n, n_clusters):
+def clustering_for_n_clusters(n, n_clusters, aggregate_renewables=True):
+    aggregate_generators_carriers = (None if aggregate_renewables
+                                     else (pd.Index(n.generators.carrier.unique())
+                                           .difference(['onwind', 'offwind', 'solar'])))
     clustering = get_clustering_from_busmap(
         n, busmap_for_n_clusters(n, n_clusters),
         bus_strategies=dict(country=_make_consense("Bus", "country")),
         aggregate_generators_weighted=True,
+        aggregate_generators_carriers=aggregate_generators_carriers,
         aggregate_one_ports=["Load", "StorageUnit"]
     )
 
@@ -153,8 +186,14 @@ if __name__ == "__main__":
 
     n = pypsa.Network(snakemake.input.network)
 
+    if snakemake.wildcards.clusters.endswith('m'):
+        n_clusters = int(snakemake.wildcards.clusters[:-1])
+        aggregate_renewables = False
+    else:
         n_clusters = int(snakemake.wildcards.clusters)
-    clustering = clustering_for_n_clusters(n, n_clusters)
+        aggregate_renewables = True
+
+    clustering = clustering_for_n_clusters(n, n_clusters, aggregate_renewables)
 
     clustering.network.export_to_netcdf(snakemake.output.network)
 

scripts/simplify_network.py

@@ -208,7 +208,7 @@ if __name__ == "__main__":
 
     if snakemake.wildcards.simpl:
         n_clusters = int(snakemake.wildcards.simpl)
-        clustering = clustering_for_n_clusters(n_clusters)
+        clustering = clustering_for_n_clusters(n, n_clusters)
 
         n = clustering.network
         busmaps.append(clustering.busmap)