Merge branch 'windcosts'

Snakefile

@@ -150,6 +150,7 @@ rule add_electricity:
 rule simplify_network:
     input:
         network='networks/{network}.nc',
+        tech_costs=COSTS,
         regions_onshore="resources/regions_onshore.geojson",
         regions_offshore="resources/regions_offshore.geojson"
     output:

cluster.yaml

@@ -5,7 +5,7 @@ feedin_preparation:
   walltime: "12:00:00"
 
 solve_network:
-  walltime: "02:00:00:00"
+  walltime: "05:00:00:00"
 
 solve:
   walltime: "05:00:00:00"

data/costs.csv

@@ -1,7 +1,7 @@
 technology,year,parameter,value,unit,source
 solar-rooftop,2030,discount rate,0.04,per unit,standard for decentral
-onwind,2030,lifetime,25,years,IEA2010
-offwind,2030,lifetime,25,years,IEA2010
+onwind,2030,lifetime,30,years,DEA https://ens.dk/en/our-services/projections-and-models/technology-data
+offwind,2030,lifetime,30,years,DEA https://ens.dk/en/our-services/projections-and-models/technology-data
 solar,2030,lifetime,25,years,IEA2010
 solar-rooftop,2030,lifetime,25,years,IEA2010
 solar-utility,2030,lifetime,25,years,IEA2010
@@ -16,8 +16,10 @@ lignite,2030,lifetime,40,years,IEA2010
 geothermal,2030,lifetime,40,years,IEA2010
 biomass,2030,lifetime,30,years,ECF2010 in DIW DataDoc http://hdl.handle.net/10419/80348
 oil,2030,lifetime,30,years,ECF2010 in DIW DataDoc http://hdl.handle.net/10419/80348
-onwind,2030,investment,1182,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348
-offwind,2030,investment,2506,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348
+onwind,2030,investment,910,EUR/kWel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data
+offwind,2030,investment,1640,EUR/kWel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data
+offwind-grid,2030,investment,255,EUR/kWel,Haertel 2017; assuming one onshore and one offshore node
+offwind-grid-perlength,2030,investment,0.97,EUR/kWel/km,Haertel 2017
 solar,2030,investment,600,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348
 biomass,2030,investment,2209,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348
 geothermal,2030,investment,3392,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348
@@ -32,8 +34,8 @@ OCGT,2030,investment,400,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348
 nuclear,2030,investment,6000,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348
 CCGT,2030,investment,800,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348
 oil,2030,investment,400,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348
-onwind,2030,FOM,2.961083,%/year,DIW DataDoc http://hdl.handle.net/10419/80348
-offwind,2030,FOM,3.192338,%/year,DIW DataDoc http://hdl.handle.net/10419/80348
+onwind,2030,FOM,2.450549,%/year,DEA https://ens.dk/en/our-services/projections-and-models/technology-data
+offwind,2030,FOM,2.304878,%/year,DEA https://ens.dk/en/our-services/projections-and-models/technology-data
 solar,2030,FOM,4.166667,%/year,DIW DataDoc http://hdl.handle.net/10419/80348
 solar-rooftop,2030,FOM,2,%/year,ETIP PV
 solar-utility,2030,FOM,3,%/year,ETIP PV
@@ -47,8 +49,8 @@ hydro,2030,FOM,1,%/year,DIW DataDoc http://hdl.handle.net/10419/80348
 ror,2030,FOM,2,%/year,DIW DataDoc http://hdl.handle.net/10419/80348
 CCGT,2030,FOM,2.5,%/year,DIW DataDoc http://hdl.handle.net/10419/80348
 OCGT,2030,FOM,3.75,%/year,DIW DataDoc http://hdl.handle.net/10419/80348
-onwind,2030,VOM,0.015,EUR/MWhel,RES costs made up to fix curtailment order
-offwind,2030,VOM,0.02,EUR/MWhel,RES costs made up to fix curtailment order
+onwind,2030,VOM,2.3,EUR/MWhel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data
+offwind,2030,VOM,2.7,EUR/MWhel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data
 solar,2030,VOM,0.01,EUR/MWhel,RES costs made up to fix curtailment order
 coal,2030,VOM,6,EUR/MWhel,DIW DataDoc http://hdl.handle.net/10419/80348 PC (Advanced/SuperC)
 lignite,2030,VOM,7,EUR/MWhel,DIW DataDoc http://hdl.handle.net/10419/80348

scripts/add_electricity.py

@@ -161,6 +161,15 @@ def attach_wind_and_solar(n, costs):
 
         n.add("Carrier", name=tech)
         with xr.open_dataset(getattr(snakemake.input, 'profile_' + tech)) as ds:
+            capital_cost = costs.at[tech, 'capital_cost']
+            if tech + "-grid" in costs.index:
+                if tech + "-grid-perlength" in costs.index:
+                    grid_cost = costs.at[tech + "-grid", "capital_cost"] + costs.at[tech + "-grid-perlength", 'capital_cost'] * ds['average_distance'].to_pandas()
+                    logger.info("Added connection cost of {:0.0f}-{:0.0f} Eur/MW/a to {}".format(grid_cost.min(), grid_cost.max(), tech))
+                else:
+                    grid_cost = costs.at[tech + "-grid", "capital_cost"]
+                    logger.info("Added connection cost of {:0.0f} Eur/MW/a to {}".format(grid_cost, tech))
+                capital_cost = capital_cost + grid_cost
 
             n.madd("Generator", ds.indexes['bus'], ' ' + tech,
                    bus=ds.indexes['bus'],
@@ -169,7 +178,7 @@ def attach_wind_and_solar(n, costs):
                    p_nom_max=ds['p_nom_max'].to_pandas(),
                    weight=ds['weight'].to_pandas(),
                    marginal_cost=costs.at[tech, 'marginal_cost'],
-                   capital_cost=costs.at[tech, 'capital_cost'],
+                   capital_cost=capital_cost,
                    efficiency=costs.at[tech, 'efficiency'],
                    p_max_pu=ds['profile'].transpose('time', 'bus').to_pandas())
 

scripts/build_bus_regions.py

@@ -24,6 +24,8 @@ for country in countries:
     onshore_shape = country_shapes[country]
     onshore_locs = n.buses.loc[c_b & n.buses.substation_lv, ["x", "y"]]
     onshore_regions.append(gpd.GeoDataFrame({
+            'x': onshore_locs['x'],
+            'y': onshore_locs['y'],
             'geometry': voronoi_partition_pts(onshore_locs.values, onshore_shape),
             'country': country
         }, index=onshore_locs.index))
@@ -32,6 +34,8 @@ for country in countries:
     offshore_shape = offshore_shapes[country]
     offshore_locs = n.buses.loc[c_b & n.buses.substation_off, ["x", "y"]]
     offshore_regions_c = gpd.GeoDataFrame({
+            'x': offshore_locs['x'],
+            'y': offshore_locs['y'],
             'geometry': voronoi_partition_pts(offshore_locs.values, offshore_shape),
             'country': country
         }, index=offshore_locs.index)

scripts/build_renewable_profiles.py

@@ -12,6 +12,7 @@ import geokit as gk
 from osgeo import gdal
 from scipy.sparse import csr_matrix, vstack
 
+from pypsa.geo import haversine
 from vresutils import landuse as vlanduse
 from vresutils.array import spdiag
 
@@ -98,8 +99,8 @@ if __name__ == '__main__':
 
     with Pool(initializer=init_globals, initargs=(bounds, dx, dy),
               maxtasksperchild=20, processes=snakemake.config['atlite'].get('nprocesses', 2)) as pool:
-        features = gk.vector.extractFeatures(snakemake.input.regions, onlyAttr=True) #.iloc[:10]
-        buses = pd.Index(features['name'], name="bus")
+        regions = gk.vector.extractFeatures(snakemake.input.regions, onlyAttr=True) #.iloc[:10]
+        buses = pd.Index(regions['name'], name="bus")
         widgets = [
             pgb.widgets.Percentage(),
             ' ', pgb.widgets.SimpleProgress(format='(%s)' % pgb.widgets.SimpleProgress.DEFAULT_FORMAT),
@@ -154,9 +155,20 @@ if __name__ == '__main__':
     layout = xr.DataArray(np.asarray(potmatrix.sum(axis=0)).reshape(cutout.shape),
                           [cutout.meta.indexes[ax] for ax in ['y', 'x']])
 
+    # Determine weighted average distance from substation
+    cell_coords = cutout.grid_coordinates()
+
+    average_distance = []
+    for i in regions.index:
+        row = layoutmatrix[i]
+        distances = haversine(regions.loc[i, ['x', 'y']], cell_coords[row.indices])[0]
+        average_distance.append((distances * (row.data / row.data.sum())).sum())
+    average_distance = xr.DataArray(average_distance, [buses])
+
     ds = xr.merge([(correction_factor * profile).rename('profile'),
-                   capacities.rename('weight'),
-                   p_nom_max.rename('p_nom_max'),
-                   layout.rename('potential')])
+                capacities.rename('weight'),
+                p_nom_max.rename('p_nom_max'),
+                layout.rename('potential'),
+                average_distance.rename('average_distance')])
     (ds.sel(bus=(ds['profile'].mean('time') > config.get('min_p_max_pu', 0.)) & (ds['p_nom_max'] > 0.))
      .to_netcdf(snakemake.output.profile))

scripts/simplify_network.py

@@ -10,7 +10,7 @@ import os
 import re
 import numpy as np
 import scipy as sp
-from scipy.sparse.csgraph import connected_components
+from scipy.sparse.csgraph import connected_components, dijkstra
 import xarray as xr
 import geopandas as gpd
 import shapely
@@ -26,6 +26,7 @@ from pypsa.networkclustering import (busmap_by_stubs, busmap_by_kmeans,
                                      aggregategenerators, aggregateoneport)
 
 from cluster_network import clustering_for_n_clusters, cluster_regions
+from add_electricity import load_costs
 
 def simplify_network_to_380(n):
     ## All goes to v_nom == 380
@@ -62,7 +63,22 @@ def simplify_network_to_380(n):
 
     return n, trafo_map
 
-def _aggregate_and_move_components(n, busmap, aggregate_one_ports={"Load", "StorageUnit"}):
+def _adjust_costs_using_distance(n, distance):
+    costs = load_costs(n.snapshot_weightings.sum() / 8760, snakemake.input.tech_costs,
+                       snakemake.config['costs'], snakemake.config['electricity'])
+
+    for tech in snakemake.config['renewable']:
+        if tech + "-grid-perlength" in costs.index:
+            cost_perlength = costs.at[tech + "-grid-perlength", "capital_cost"]
+            tech_b = n.generators.carrier == tech
+            generator_distance = n.generators.loc[tech_b, "bus"].map(distance).loc[lambda s: s>0]
+            if not generator_distance.empty:
+                n.generators.loc[generator_distance.index, "capital_cost"] += cost_perlength * generator_distance
+                logger.info("Displacing generator(s) {}; capital_cost is adjusted accordingly"
+                            .format(", ".join("`{}` by {:.0f}km".format(b, d) for b, d in generator_distance.iteritems())))
+
+
+def _aggregate_and_move_components(n, busmap, distance, aggregate_one_ports={"Load", "StorageUnit"}):
     def replace_components(n, c, df, pnl):
         n.mremove(c, n.df(c).index)
 
@@ -71,6 +87,8 @@ def _aggregate_and_move_components(n, busmap, aggregate_one_ports={"Load", "Stor
             if not df.empty:
                 import_series_from_dataframe(n, df, c, attr)
 
+    _adjust_costs_using_distance(n, distance)
+
     generators, generators_pnl = aggregategenerators(n, busmap)
     replace_components(n, "Generator", generators, generators_pnl)
 
@@ -84,6 +102,16 @@ def _aggregate_and_move_components(n, busmap, aggregate_one_ports={"Load", "Stor
         df = n.df(c)
         n.mremove(c, df.index[df.bus0.isin(buses_to_del) | df.bus1.isin(buses_to_del)])
 
+def _compute_distance(n, busmap, buses=None, adjacency_matrix=None):
+    if buses is None:
+        buses = busmap.index[busmap.index != busmap.values]
+
+    if adjacency_matrix is None:
+        adjacency_matrix = n.adjacency_matrix(weights=pd.concat(dict(Link=n.links.length, Line=pd.Series(0., n.lines.index))))
+
+    dist = dijkstra(adjacency_matrix, directed=False, indices=n.buses.index.get_indexer(buses))
+    return pd.Series(dist[np.arange(len(buses)), n.buses.index.get_indexer(busmap.loc[buses])], buses)
+
 def simplify_links(n):
     ## Complex multi-node links are folded into end-points
     logger.info("Simplifying connected link components")
@@ -127,6 +155,8 @@ def simplify_links(n):
             seen.add(u)
 
     busmap = n.buses.index.to_series()
+    distance = pd.Series(0., n.buses.index)
+    adjacency_matrix = n.adjacency_matrix(weights=pd.concat(dict(Link=n.links.length, Line=pd.Series(0., n.lines.index))))
 
     for lbl in labels.value_counts().loc[lambda s: s > 2].index:
 
@@ -139,6 +169,8 @@ def simplify_links(n):
             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]]
+            distance.loc[buses] += _compute_distance(n, busmap, buses)
+
             all_links = [i for _, i in sum(links, [])]
 
             p_max_pu = snakemake.config['links'].get('p_max_pu', 1.)
@@ -168,14 +200,17 @@ def simplify_links(n):
 
     logger.debug("Collecting all components using the busmap")
 
-    _aggregate_and_move_components(n, busmap)
+    _aggregate_and_move_components(n, busmap, distance)
     return n, busmap
 
 def remove_stubs(n):
     logger.info("Removing stubs")
 
     busmap = busmap_by_stubs(n) #  ['country'])
-    _aggregate_and_move_components(n, busmap)
+
+    distance = _compute_distance(n, busmap)
+
+    _aggregate_and_move_components(n, busmap, distance)
 
     return n, busmap
 
@@ -199,8 +234,7 @@ if __name__ == "__main__":
             )
         )
 
-    logger = logging.getLogger()
-    logger.setLevel(snakemake.config['logging_level'])
+    logging.basicConfig(level=snakemake.config['logging_level'])
 
     n = pypsa.Network(snakemake.input.network)