diff --git a/Snakefile b/Snakefile index 0c9d22a9..a37f6b32 100644 --- a/Snakefile +++ b/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: 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" 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 04565035..ac4362f9 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 5ee7c7c0..9cc0974e 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 7abd4202..3ffc0f1c 100644 --- a/scripts/build_renewable_profiles.py +++ b/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)) diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index cf807b7f..14811a64 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"].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)