From b30f4a2fff71a420eba3b4682359d64627746d3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6rsch?= Date: Fri, 26 Oct 2018 21:32:04 +0200 Subject: [PATCH 1/6] Update technology costs for wind Calculate grid extension costs for offshore based on weighted average distance from weather cell to substation. --- Snakefile | 1 - data/costs.csv | 18 ++++++++++-------- scripts/add_electricity.py | 11 ++++++++++- scripts/build_bus_regions.py | 4 ++++ scripts/build_renewable_profiles.py | 16 +++++++++++++++- 5 files changed, 39 insertions(+), 11 deletions(-) diff --git a/Snakefile b/Snakefile index 97e2aa0a..ccae9908 100644 --- a/Snakefile +++ b/Snakefile @@ -112,7 +112,6 @@ rule build_renewable_potentials: rule build_renewable_profiles: input: - base_network="networks/base.nc", potentials="resources/potentials_{technology}.nc", regions=lambda wildcards: ("resources/regions_onshore.geojson" if wildcards.technology in ('onwind', 'solar') diff --git a/data/costs.csv b/data/costs.csv index 5fdcea72..ee216d36 100644 --- a/data/costs.csv +++ b/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 diff --git a/scripts/add_electricity.py b/scripts/add_electricity.py index 2a460db1..465f6414 100644 --- a/scripts/add_electricity.py +++ b/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()) diff --git a/scripts/build_bus_regions.py b/scripts/build_bus_regions.py index 6ba67b50..7034df61 100644 --- a/scripts/build_bus_regions.py +++ b/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) diff --git a/scripts/build_renewable_profiles.py b/scripts/build_renewable_profiles.py index 1fe6391e..b7e00cba 100644 --- a/scripts/build_renewable_profiles.py +++ b/scripts/build_renewable_profiles.py @@ -6,6 +6,8 @@ import numpy as np import xarray as xr import pandas as pd import geopandas as gpd +from pypsa.geo import haversine +from vresutils.array import spdiag import logging logger = logging.getLogger(__name__) @@ -47,9 +49,21 @@ p_nom_max = xr.DataArray([np.nanmin(relativepotentials[row.nonzero()[1]]) for row in indicatormatrix.tocsr()], [capacities.coords['bus']]) * capacities +# Determine weighted average distance to substation +matrix_weighted = indicatormatrix * spdiag(layout.stack(spatial=('y', 'x'))) +cell_coords = cutout.grid_coordinates() + +average_distance = [] +for i, bus in enumerate(regions.index): + row = matrix_weighted[i] + distances = haversine(regions.loc[bus, ['x', 'y']], cell_coords[row.indices])[0] + average_distance.append((distances * (row.data / row.data.sum())).sum()) +average_distance = xr.DataArray(average_distance, [regions.index.rename("bus")]) + ds = xr.merge([(correction_factor * profile).rename('profile'), capacities.rename('weight'), p_nom_max.rename('p_nom_max'), - layout.rename('potential')]) + layout.rename('potential'), + average_distance.rename('average_distance')]) (ds.sel(bus=ds['profile'].mean('time') > config.get('min_p_max_pu', 0.)) .to_netcdf(snakemake.output.profile)) From f87d7d9f94657f54c125692b28c81f910bd9f9d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6rsch?= Date: Sat, 27 Oct 2018 03:34:40 +0200 Subject: [PATCH 2/6] base_network: Important small fixes - Allow connecting at substations - Fix offshore shape selection --- scripts/base_network.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/base_network.py b/scripts/base_network.py index 4df0e925..44a8ddba 100644 --- a/scripts/base_network.py +++ b/scripts/base_network.py @@ -279,7 +279,7 @@ def _set_countries_and_substations(n): countries = snakemake.config['countries'] country_shapes = gpd.read_file(snakemake.input.country_shapes).set_index('id')['geometry'] offshore_shapes = gpd.read_file(snakemake.input.offshore_shapes).set_index('id')['geometry'] - substation_b = buses['symbol'].str.contains('substation', case=False) + substation_b = buses['symbol'].str.contains('substation|converter station', case=False) def prefer_voltage(x, which): index = x.index @@ -307,7 +307,7 @@ def _set_countries_and_substations(n): buses.loc[onshore_country_b, 'country'] = country - if country not in offshore_shapes: continue + if country not in offshore_shapes.index: continue offshore_country_b = buses_in_shape(offshore_shapes[country]) offshore_b |= offshore_country_b @@ -325,7 +325,7 @@ def _set_countries_and_substations(n): c_nan_b = buses.country.isnull() if c_nan_b.sum() > 0: c_tag = _get_country(buses.loc[c_nan_b]) - c_tag.loc[c_tag.index.difference(countries)] = np.nan + c_tag.loc[~c_tag.isin(countries)] = np.nan n.buses.loc[c_nan_b, 'country'] = c_tag c_tag_nan_b = n.buses.country.isnull() From 9c64034950633a473572e2b378826b3718aa3b4e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6rsch?= Date: Sat, 27 Oct 2018 03:41:38 +0200 Subject: [PATCH 3/6] simplify_network: Add simplification cost updates Depends on pypsa commit https://github.com/PyPSA/PyPSA/commit/9ec406252add1b6 --- Snakefile | 1 + scripts/simplify_network.py | 40 +++++++++++++++++++++++++++++++++---- 2 files changed, 37 insertions(+), 4 deletions(-) diff --git a/Snakefile b/Snakefile index ccae9908..43dc919d 100644 --- a/Snakefile +++ b/Snakefile @@ -155,6 +155,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: diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index 576b3ae5..c7fb1d76 100644 --- a/scripts/simplify_network.py +++ b/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"]].join(distance.rename("distance"), on="bus")['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) @@ -127,6 +145,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 +159,9 @@ 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]] + dist = dijkstra(adjacency_matrix, directed=False, indices=n.buses.index.get_indexer(buses)) + distance.loc[buses] += [dist[i,j] for i, j in enumerate(n.buses.index.get_indexer(busmap.loc[buses]))] + all_links = [i for _, i in sum(links, [])] p_max_pu = snakemake.config['links'].get('p_max_pu', 1.) @@ -168,14 +191,23 @@ 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) + indices, = np.where(busmap.index != busmap.values) + buses = busmap.index[indices] + + adjacency_matrix = n.adjacency_matrix(busorder=busmap.index, weights=pd.concat(dict(Link=n.links.length, Line=pd.Series(0., n.lines.index)))) + dist = dijkstra(adjacency_matrix, directed=False, indices=indices) + distance = pd.Series(dist[np.arange(len(indices)), busmap.index.get_indexer(busmap.iloc[indices])], buses) + logger.info("The following offshore buses are displaced: {}" + .format(", ".join("{} by {:.0f}".format(b, d) for b,d in distance.loc[offwind_buses].iteritems()))) + + _aggregate_and_move_components(n, busmap, distance) return n, busmap From cc7d356db05169de272fbb2415e9221041474b31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6rsch?= Date: Sat, 27 Oct 2018 03:45:28 +0200 Subject: [PATCH 4/6] simplify_network: Set up logger correctly --- scripts/simplify_network.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index c7fb1d76..e1ca1476 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -231,8 +231,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) From 9df54c9530561d2d085d16ac662d43cc1610003c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6rsch?= Date: Sat, 27 Oct 2018 16:42:58 +0200 Subject: [PATCH 5/6] simplify_network: Refactor compute_distance --- scripts/simplify_network.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index e1ca1476..ab118f8a 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -71,7 +71,7 @@ def _adjust_costs_using_distance(n, distance): 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"]].join(distance.rename("distance"), on="bus")['distance'].loc[lambda s: s>0] + 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" @@ -102,6 +102,16 @@ def _aggregate_and_move_components(n, busmap, distance, aggregate_one_ports={"Lo 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") @@ -159,8 +169,7 @@ 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]] - dist = dijkstra(adjacency_matrix, directed=False, indices=n.buses.index.get_indexer(buses)) - distance.loc[buses] += [dist[i,j] for i, j in enumerate(n.buses.index.get_indexer(busmap.loc[buses]))] + distance.loc[buses] += _compute_distance(n, busmap, buses) all_links = [i for _, i in sum(links, [])] @@ -198,14 +207,8 @@ def remove_stubs(n): logger.info("Removing stubs") busmap = busmap_by_stubs(n, ['country']) - indices, = np.where(busmap.index != busmap.values) - buses = busmap.index[indices] - adjacency_matrix = n.adjacency_matrix(busorder=busmap.index, weights=pd.concat(dict(Link=n.links.length, Line=pd.Series(0., n.lines.index)))) - dist = dijkstra(adjacency_matrix, directed=False, indices=indices) - distance = pd.Series(dist[np.arange(len(indices)), busmap.index.get_indexer(busmap.iloc[indices])], buses) - logger.info("The following offshore buses are displaced: {}" - .format(", ".join("{} by {:.0f}".format(b, d) for b,d in distance.loc[offwind_buses].iteritems()))) + distance = _compute_distance(n, busmap) _aggregate_and_move_components(n, busmap, distance) From e262f2264dd2292de150c09633ebbaf6efb3080d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jonas=20H=C3=B6rsch?= Date: Wed, 31 Oct 2018 15:21:14 +0100 Subject: [PATCH 6/6] Increase walltime for long-running models --- cluster.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cluster.yaml b/cluster.yaml index 0f99a28d..02bb07ab 100644 --- a/cluster.yaml +++ b/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"