From 3ae8021cb69f18d7a80efc6d92b082ccf06f469a Mon Sep 17 00:00:00 2001 From: martavp Date: Sat, 26 Dec 2020 17:47:32 +0100 Subject: [PATCH 01/10] Adapt nomenclature from "YearCommissioned" to "DataIn" This was breaking add_existing_baseyear.py and it is now fixed. Column name for commissioning year in powerplantmatching has changed. Now "DataIn" is used as column name, also when renewable capacities per country are added to the power plants dataframe --- scripts/add_existing_baseyear.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py index 852c5bb0..09b47da6 100644 --- a/scripts/add_existing_baseyear.py +++ b/scripts/add_existing_baseyear.py @@ -125,7 +125,7 @@ def add_existing_renewables(df_agg): if capacity > 0.: df_agg.at[name,"Fueltype"] = tech df_agg.at[name,"Capacity"] = capacity - df_agg.at[name,"YearCommissioned"] = year + df_agg.at[name,"DateIn"] = year df_agg.at[name,"cluster_bus"] = node def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, baseyear): @@ -182,7 +182,7 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas add_existing_renewables(df_agg) df_agg["grouping_year"] = np.take(grouping_years, - np.digitize(df_agg.YearCommissioned, + np.digitize(df_agg.DateIn, grouping_years, right=True)) @@ -249,7 +249,7 @@ def add_heating_capacities_installed_before_baseyear(n, baseyear, grouping_years grouping_years : intervals to group existing capacities - linear decomissioning of heating capacities from 2020 to 2045 is + linear decommissioning of heating capacities from 2020 to 2045 is currently assumed heating capacities split between residential and services proportional @@ -408,18 +408,18 @@ if __name__ == "__main__": if 'snakemake' not in globals(): from vresutils.snakemake import MockSnakemake snakemake = MockSnakemake( - wildcards=dict(network='elec', simpl='', clusters='39', lv='1.0', - sector_opts='Co2L0-168H-T-H-B-I-solar3-dist1', - co2_budget_name='b30b3', + wildcards=dict(network='elec', simpl='', clusters='45', lv='1.0', + sector_opts='Co2L0-3H-T-H-B-I-solar3-dist1', planning_horizons='2020'), - input=dict(network='pypsa-eur-sec/results/test/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{co2_budget_name}_{planning_horizons}.nc', + input=dict(network='pypsa-eur-sec/results/version-2/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{planning_horizons}.nc', powerplants='pypsa-eur/resources/powerplants.csv', busmap_s='pypsa-eur/resources/busmap_{network}_s{simpl}.csv', busmap='pypsa-eur/resources/busmap_{network}_s{simpl}_{clusters}.csv', - costs='pypsa-eur-sec/data/costs/costs_{planning_horizons}.csv', + costs='technology_data/outputs/costs_{planning_horizons}.csv', cop_air_total="pypsa-eur-sec/resources/cop_air_total_{network}_s{simpl}_{clusters}.nc", - cop_soil_total="pypsa-eur-sec/resources/cop_soil_total_{network}_s{simpl}_{clusters}.nc"), - output=['pypsa-eur-sec/results/test/prenetworks_brownfield/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{planning_horizons}.nc'], + cop_soil_total="pypsa-eur-sec/resources/cop_soil_total_{network}_s{simpl}_{clusters}.nc", + clustered_pop_layout="pypsa-eur-sec/resources/pop_layout_{network}_s{simpl}_{clusters}.csv",), + output=['pypsa-eur-sec/results/version-2/prenetworks_brownfield/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{planning_horizons}.nc'], ) import yaml with open('config.yaml', encoding='utf8') as f: From ec5efcc93d02ef21df0bda976d9cd642804f20fc Mon Sep 17 00:00:00 2001 From: martavp Date: Mon, 28 Dec 2020 15:39:05 +0100 Subject: [PATCH 02/10] Make explicit solver_dir=tmpdir Running the rule solve_network in the university cluster, I was getting a "No space left on device" error. Making solve_dir=tmpdir by default avoids this error and makes it easier to identify any problem with temp files --- scripts/solve_network.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 795c3327..85251caa 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -272,6 +272,7 @@ def solve_network(n, config=None, solver_log=None, opts=None): solver_name=solver_name, solver_logfile=solver_log, solver_options=solver_options, + solver_dir=tmpdir, extra_functionality=extra_functionality, formulation=solve_opts['formulation']) #extra_postprocessing=extra_postprocessing From 5db196f92ad4c353d86cf349e87ee61fbdbc9990 Mon Sep 17 00:00:00 2001 From: martavp Date: Tue, 29 Dec 2020 11:31:00 +0100 Subject: [PATCH 03/10] build_eea_co2() now reads version23 of UNFCCC inventory. **This requires updating to UNFCCC_v23.csv in the data bundle** This enables that the same function is used to read emissions in 2018, which are assumed to remain constant in 2019 and 2020 and subtracted from carbon budget (estimated from 2018 on). I checked and 1990 emissions calculated with UNFCCC_v23 are very similar to those calculated with UNFCCC_v21 (<1% differences in all values at EU level). Some countries show higher deviations, mainly in domestic aviation and navigation. I guess because those values started to be reported later and v23 should include more accurate values. --- scripts/build_energy_totals.py | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/scripts/build_energy_totals.py b/scripts/build_energy_totals.py index 0cfa709b..4d80abb9 100644 --- a/scripts/build_energy_totals.py +++ b/scripts/build_energy_totals.py @@ -390,12 +390,12 @@ def build_energy_totals(): return clean_df -def build_eea_co2(): +def build_eea_co2(year=1990): # see ../notebooks/compute_1990_Europe_emissions_for_targets.ipynb - #https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-14 - #downloaded 190222 (modified by EEA last on 181130) - fn = "data/eea/UNFCCC_v21.csv" + #https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-16 + #downloaded 201228 (modified by EEA last on 201221) + fn = "data/eea/UNFCCC_v23.csv" df = pd.read_csv(fn, encoding="latin-1") df.loc[df["Year"] == "1985-1987","Year"] = 1986 df["Year"] = df["Year"].astype(int) @@ -418,16 +418,14 @@ def build_eea_co2(): e['waste management'] = '5 - Waste management' e['other'] = '6 - Other Sector' e['indirect'] = 'ind_CO2 - Indirect CO2' - e["total wL"] = "Total (with LULUCF, with indirect CO2)" - e["total woL"] = "Total (without LULUCF, with indirect CO2)" + e["total wL"] = "Total (with LULUCF)" + e["total woL"] = "Total (without LULUCF)" pol = "CO2" #["All greenhouse gases - (CO2 equivalent)","CO2"] cts = ["CH","EUA","NO"] + eu28_eea - year = 1990 - emissions = df.loc[idx[cts,pol,year,e.values],"emissions"].unstack("Sector_name").rename(columns=pd.Series(e.index,e.values)).rename(index={"All greenhouse gases - (CO2 equivalent)" : "GHG"},level=1) #only take level 0, since level 1 (pol) and level 2 (year) are trivial @@ -467,7 +465,7 @@ def build_eurostat_co2(year=1990): return eurostat_co2 -def build_co2_totals(year=1990): +def build_co2_totals(eea_co2, eurostat_co2, year=1990): co2 = eea_co2.reindex(["EU28","NO","CH","BA","RS","AL","ME","MK"] + eu28) @@ -486,10 +484,6 @@ def build_co2_totals(year=1990): #doesn't include non-energy emissions co2.loc[ct,'agriculture'] = eurostat_co2[ct,"+","+","Agriculture / Forestry"].sum() - - - co2.to_csv(snakemake.output.co2_name) - return co2 @@ -547,7 +541,7 @@ if __name__ == "__main__": snakemake.output['transport_name'] = "data/transport_data.csv" snakemake.input = Dict() - snakemake.input['nuts3_shapes'] = 'resources/nuts3_shapes.geojson' + snakemake.input['nuts3_shapes'] = '../pypsa-eur/resources/nuts3_shapes.geojson' nuts3 = gpd.read_file(snakemake.input.nuts3_shapes).set_index('index') population = nuts3['pop'].groupby(nuts3.country).sum() @@ -566,6 +560,7 @@ if __name__ == "__main__": eurostat_co2 = build_eurostat_co2() - build_co2_totals() - + co2=build_co2_totals(eea_co2, eurostat_co2, year) + co2.to_csv(snakemake.output.co2_name) + build_transport_data() From 7479ba0424ad95d531094ee7313d76e83e05c281 Mon Sep 17 00:00:00 2001 From: martavp Date: Wed, 30 Dec 2020 12:14:08 +0100 Subject: [PATCH 04/10] Fix unicode error due to dash before sawdust A quick fix to https://github.com/PyPSA/pypsa-eur-sec/issues/79 --- config.default.yaml | 2 +- scripts/build_biomass_potentials.py | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index d43d0e29..7907859a 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -69,7 +69,7 @@ biomass: year: 2030 scenario: "Med" classes: - solid biomass: ['Primary agricultural residues', 'Forestry energy residue', 'Secondary forestry residues', 'Secondary Forestry residues – sawdust', 'Forestry residues from landscape care biomass', 'Municipal waste'] + solid biomass: ['Primary agricultural residues', 'Forestry energy residue', 'Secondary forestry residues', 'Secondary Forestry residues sawdust', 'Forestry residues from landscape care biomass', 'Municipal waste'] not included: ['Bioethanol sugar beet biomass', 'Rapeseeds for biodiesel', 'sunflower and soya for Biodiesel', 'Starchy crops biomass', 'Grassy crops biomass', 'Willow biomass', 'Poplar biomass potential', 'Roundwood fuelwood', 'Roundwood Chips & Pellets'] biogas: ['Manure biomass potential', 'Sludge biomass'] diff --git a/scripts/build_biomass_potentials.py b/scripts/build_biomass_potentials.py index c959680f..44fd04b5 100644 --- a/scripts/build_biomass_potentials.py +++ b/scripts/build_biomass_potentials.py @@ -1,3 +1,4 @@ +# coding: utf-8 import pandas as pd @@ -57,7 +58,12 @@ if __name__ == "__main__": snakemake.input['jrc_potentials'] = "data/biomass/JRC Biomass Potentials.xlsx" snakemake.output = Dict() snakemake.output['biomass_potentials'] = 'data/biomass_potentials.csv' + snakemake.output['biomass_potentials_all']='resources/biomass_potentials_all.csv' with open('config.yaml', encoding='utf8') as f: snakemake.config = yaml.safe_load(f) - + + if 'Secondary Forestry residues sawdust' in snakemake.config['biomass']['classes']['solid biomass']: + snakemake.config['biomass']['classes']['solid biomass'].remove('Secondary Forestry residues sawdust') + snakemake.config['biomass']['classes']['solid biomass'].append('Secondary Forestry residues – sawdust') + build_biomass_potentials() From 3457d4ee386e5760eca5c1e73e5993ae3d385336 Mon Sep 17 00:00:00 2001 From: martavp Date: Wed, 30 Dec 2020 15:55:08 +0100 Subject: [PATCH 05/10] Add function build_carbon_budget() For the myopic method, based on the carbon budget indicated in the config file (sector_opts), a CO2 limit is calculated for every planning_horizon following an exponential or beta decay. A file with CO2 limit in every planning_horizon and a plot showing historical and planned CO2 emissions are saved in the results --- Snakefile | 1 + config.default.yaml | 4 + scripts/prepare_sector_network.py | 226 +++++++++++++++++++++++++++++- 3 files changed, 226 insertions(+), 5 deletions(-) diff --git a/Snakefile b/Snakefile index ee0fb8cf..0ac8bd99 100644 --- a/Snakefile +++ b/Snakefile @@ -292,6 +292,7 @@ rule build_retro_cost: output: retro_cost="resources/retro_cost_{network}_s{simpl}_{clusters}.csv", floor_area="resources/floor_area_{network}_s{simpl}_{clusters}.csv" + resources: mem_mb=1000 script: "scripts/build_retro_cost.py" diff --git a/config.default.yaml b/config.default.yaml index 7907859a..08cda6c4 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -24,11 +24,15 @@ scenario: # B for biomass supply, I for industry, shipping and aviation # solarx or onwindx changes the available installable potential by factor x # dist{n} includes distribution grids with investment cost of n times cost in data/costs.csv + # for myopic/perfect foresight cb states the carbon budget in GtCO2 (cumulative + # emissions throughout the transition path in the timeframe determined by the + # planning_horizons), be:beta decay; ex:exponential decay planning_horizons : [2030] # investment years for myopic and perfect; or costs year for overnight # for example, set to [2020, 2030, 2040, 2050] for myopic foresight # CO2 budget as a fraction of 1990 emissions # this is over-ridden if CO2Lx is set in sector_opts +# this is also over-ridden if cb is set in sector_opts co2_budget: 2020: 0.7011648746 2025: 0.5241935484 diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 34a02e0f..088b124f 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -20,6 +20,8 @@ import pytz from vresutils.costdata import annuity +from scipy.stats import beta +from build_energy_totals import build_eea_co2, build_eurostat_co2, build_co2_totals #First tell PyPSA that links can have multiple outputs by #overriding the component_attrs. This can be done for @@ -45,6 +47,206 @@ override_component_attrs["Generator"].loc["lifetime"] = ["float","years",np.nan, override_component_attrs["Store"].loc["build_year"] = ["integer","year",np.nan,"build year","Input (optional)"] override_component_attrs["Store"].loc["lifetime"] = ["float","years",np.nan,"lifetime","Input (optional)"] + +def co2_emissions_year(year): + """ + calculate co2 emissions in one specific year (e.g. 1990 or 2018). + """ + eea_co2 = build_eea_co2(year) + + #TODO: read Eurostat data from year>2014, this only affects the estimation of + # CO2 emissions for "BA","RS","AL","ME","MK" + if year > 2014: + eurostat_co2 = build_eurostat_co2(year=2014) + else: + eurostat_co2 = build_eurostat_co2(year) + + co2_totals=build_co2_totals(eea_co2, eurostat_co2, year) + + pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0) + pop_layout["ct"] = pop_layout.index.str[:2] + cts = pop_layout.ct.value_counts().index + + co2_emissions = co2_totals.loc[cts, "electricity"].sum() + + if "T" in opts: + co2_emissions += co2_totals.loc[cts, [i+ " non-elec" for i in ["rail","road"]]].sum().sum() + if "H" in opts: + co2_emissions += co2_totals.loc[cts, [i+ " non-elec" for i in ["residential","services"]]].sum().sum() + if "I" in opts: + co2_emissions += co2_totals.loc[cts, ["industrial non-elec","industrial processes", + "domestic aviation","international aviation", + "domestic navigation","international navigation"]].sum().sum() + co2_emissions *=0.001 #MtCO2 to GtCO2 + return co2_emissions + +def historical_emissions(): + """ + read historical emissions to add them to the carbon budget plot + """ + #https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-16 + #downloaded 201228 (modified by EEA last on 201221) + fn = "data/eea/UNFCCC_v23.csv" + df = pd.read_csv(fn, encoding="latin-1") + df.loc[df["Year"] == "1985-1987","Year"] = 1986 + df["Year"] = df["Year"].astype(int) + df = df.set_index(['Year', 'Sector_name', 'Country_code', 'Pollutant_name']).sort_index() + + e = pd.Series() + e["electricity"] = '1.A.1.a - Public Electricity and Heat Production' + e['residential non-elec'] = '1.A.4.b - Residential' + e['services non-elec'] = '1.A.4.a - Commercial/Institutional' + e['rail non-elec'] = "1.A.3.c - Railways" + e["road non-elec"] = '1.A.3.b - Road Transportation' + e["domestic navigation"] = "1.A.3.d - Domestic Navigation" + e['international navigation'] = '1.D.1.b - International Navigation' + e["domestic aviation"] = '1.A.3.a - Domestic Aviation' + e["international aviation"] = '1.D.1.a - International Aviation' + e['total energy'] = '1 - Energy' + e['industrial processes'] = '2 - Industrial Processes and Product Use' + e['agriculture'] = '3 - Agriculture' + e['LULUCF'] = '4 - Land Use, Land-Use Change and Forestry' + e['waste management'] = '5 - Waste management' + e['other'] = '6 - Other Sector' + e['indirect'] = 'ind_CO2 - Indirect CO2' + e["total wL"] = "Total (with LULUCF)" + e["total woL"] = "Total (without LULUCF)" + + pol = ["CO2"] # ["All greenhouse gases - (CO2 equivalent)"] + + + pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0) + pop_layout["ct"] = pop_layout.index.str[:2] + cts = pop_layout.ct.value_counts().index.to_list() + if "GB" in cts: + cts.remove("GB") + cts.append("UK") + + year = np.arange(1990,2018).tolist() + + idx = pd.IndexSlice + co2_totals = df.loc[idx[year,e.values,cts,pol],"emissions"].unstack("Year").rename(index=pd.Series(e.index,e.values)) + + co2_totals = (1/1e6)*co2_totals.groupby(level=0, axis=0).sum() #Gton CO2 + + co2_totals.loc['industrial non-elec'] = co2_totals.loc['total energy'] - co2_totals.loc[['electricity', 'services non-elec','residential non-elec', 'road non-elec', + 'rail non-elec', 'domestic aviation', 'international aviation', 'domestic navigation', + 'international navigation']].sum() + + emissions = co2_totals.loc["electricity"] + if "T" in opts: + emissions += co2_totals.loc[[i+ " non-elec" for i in ["rail","road"]]].sum() + if "H" in opts: + emissions += co2_totals.loc[[i+ " non-elec" for i in ["residential","services"]]].sum() + if "I" in opts: + emissions += co2_totals.loc[["industrial non-elec","industrial processes", + "domestic aviation","international aviation", + "domestic navigation","international navigation"]].sum() + return emissions + + +def build_carbon_budget(o): + #distribute carbon budget following beta or exponential transition path + if "be" in o: + #beta decay + carbon_budget = float(o[o.find("cb")+2:o.find("be")]) + be=float(o[o.find("be")+2:]) + if "ex" in o: + #exponential decay + carbon_budget = float(o[o.find("cb")+2:o.find("ex")]) + r=float(o[o.find("ex")+2:]) + + e_1990 = co2_emissions_year(year=1990) + + #emissions at the beginning of the path (last year available 2018) + e_0 = co2_emissions_year(year=2018) + #emissions in 2019 and 2020 assumed equal to 2018 and substracted + carbon_budget -= 2*e_0 + planning_horizons = snakemake.config['scenario']['planning_horizons'] + CO2_CAP = pd.DataFrame(index = pd.Series(data=planning_horizons, + name='planning_horizon'), + columns=pd.Series(data=[], + name='paths', + dtype='float')) + t_0 = planning_horizons[0] + if "be" in o: + #beta decay + t_f = t_0 + (2*carbon_budget/e_0).round(0) # final year in the path + #emissions (relative to 1990) + CO2_CAP[o] = [(e_0/e_1990)*(1-beta.cdf((t-t_0)/(t_f-t_0), be, be)) for t in planning_horizons] + + if "ex" in o: + #exponential decay without delay + T=carbon_budget/e_0 + m=(1+np.sqrt(1+r*T))/T + CO2_CAP[o] = [(e_0/e_1990)*(1+(m+r)*(t-t_0))*np.exp(-m*(t-t_0)) for t in planning_horizons] + + CO2_CAP.to_csv(path_cb + 'carbon_budget_distribution.csv', sep=',', + line_terminator='\n', float_format='%.3f') + + """ + Plot historical carbon emissions in the EU and decarbonization path + """ + import matplotlib.pyplot as plt + import matplotlib.gridspec as gridspec + import seaborn as sns; sns.set() + sns.set_style('ticks') + plt.style.use('seaborn-ticks') + plt.rcParams['xtick.direction'] = 'in' + plt.rcParams['ytick.direction'] = 'in' + plt.rcParams['xtick.labelsize'] = 20 + plt.rcParams['ytick.labelsize'] = 20 + + plt.figure(figsize=(10, 7)) + gs1 = gridspec.GridSpec(1, 1) + ax1 = plt.subplot(gs1[0,0]) + ax1.set_ylabel('CO$_2$ emissions (Gt per year)',fontsize=22) + ax1.set_ylim([0,5]) + ax1.set_xlim([1990,planning_horizons[-1]+1]) + ax1.plot(e_1990*CO2_CAP[o],linewidth=3, + color='dodgerblue', label=None) + + emissions = historical_emissions() + + ax1.plot(emissions, color='black', linewidth=3, label=None) + + #plot commited and uder-discussion targets + #(notice that historical emissions include all countries in the + # network, but targets refer to EU) + ax1.plot([2020],[0.8*emissions[1990]], + marker='*', markersize=12, markerfacecolor='black', + markeredgecolor='black') + + ax1.plot([2030],[0.45*emissions[1990]], + marker='*', markersize=12, markerfacecolor='white', + markeredgecolor='black') + + ax1.plot([2030],[0.6*emissions[1990]], + marker='*', markersize=12, markerfacecolor='black', + markeredgecolor='black') + + ax1.plot([2050, 2050],[x*emissions[1990] for x in [0.2, 0.05]], + color='gray', linewidth=2, marker='_', alpha=0.5) + + ax1.plot([2050],[0.01*emissions[1990]], + marker='*', markersize=12, markerfacecolor='white', + linewidth=0, markeredgecolor='black', + label='EU under-discussion target', zorder=10, + clip_on=False) + + ax1.plot([2050],[0.125*emissions[1990]],'ro', + marker='*', markersize=12, markerfacecolor='black', + markeredgecolor='black', label='EU commited target') + + ax1.legend(fancybox=True, fontsize=18, loc=(0.01,0.01), + facecolor='white', frameon=True) + + path_cb_plot = snakemake.config['results_dir'] + snakemake.config['run'] + '/graphs/' + if not os.path.exists(path_cb_plot): + os.makedirs(path_cb_plot) + print('carbon budget distribution saved to ' + path_cb_plot + 'carbon_budget_plot.pdf') + plt.savefig(path_cb_plot+'carbon_budget_plot.pdf', dpi=300) + def add_lifetime_wind_solar(n): """ Add lifetime for solar and wind generators @@ -1775,15 +1977,15 @@ def get_parameter(item): return item -#%% + if __name__ == "__main__": # Detect running outside of snakemake and mock snakemake for testing if 'snakemake' not in globals(): from vresutils.snakemake import MockSnakemake snakemake = MockSnakemake( wildcards=dict(network='elec', simpl='', clusters='37', lv='1.0', - opts='', planning_horizons='2030', co2_budget_name="go", - sector_opts='Co2L0-120H-T-H-B-I-solar3-dist1'), + opts='', planning_horizons='2020', + sector_opts='120H-T-H-B-I-solar3-dist1-cb40ex0'), input=dict( network='../pypsa-eur/networks/{network}_s{simpl}_{clusters}_ec_lv{lv}_{opts}.nc', energy_totals_name='resources/energy_totals.csv', co2_totals_name='resources/co2_totals.csv', @@ -1819,10 +2021,10 @@ if __name__ == "__main__": solar_thermal_total="resources/solar_thermal_total_{network}_s{simpl}_{clusters}.nc", solar_thermal_urban="resources/solar_thermal_urban_{network}_s{simpl}_{clusters}.nc", solar_thermal_rural="resources/solar_thermal_rural_{network}_s{simpl}_{clusters}.nc", - retro_cost_energy = "resources/retro_cost_{network}_s{simpl}_{clusters}.csv", + retro_cost_energy = "resources/retro_cost_{network}_s{simpl}_{clusters}.csv", floor_area = "resources/floor_area_{network}_s{simpl}_{clusters}.csv" ), - output=['pypsa-eur-sec/results/test/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{co2_budget_name}_{planning_horizons}.nc'] + output=['results/version-8/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{planning_horizons}.nc'] ) import yaml with open('config.yaml', encoding='utf8') as f: @@ -1926,6 +2128,20 @@ if __name__ == "__main__": limit = get_parameter(snakemake.config["co2_budget"]) print("CO2 limit set to",limit) + for o in opts: + if "cb" in o: + path_cb = snakemake.config['results_dir'] + snakemake.config['run'] + '/csvs/' + if not os.path.exists(path_cb): + os.makedirs(path_cb) + try: + CO2_CAP=pd.read_csv(path_cb + 'carbon_budget_distribution.csv', index_col=0) + except: + build_carbon_budget(o) + CO2_CAP=pd.read_csv(path_cb + 'carbon_budget_distribution.csv', index_col=0) + + limit=CO2_CAP.loc[investment_year] + print("overriding CO2 limit with scenario limit",limit) + for o in opts: if "Co2L" in o: limit = o[o.find("Co2L")+4:] From 7c464f166b468c3ae64fb725036f83b7f2471a04 Mon Sep 17 00:00:00 2001 From: martavp Date: Wed, 30 Dec 2020 15:56:34 +0100 Subject: [PATCH 06/10] Calculate and save cumulative costs for the myopic approach Creates an additional file in results/csvs including the cumulative costs with different social discount rates --- scripts/make_summary.py | 33 ++++++++++++++++++++++++++++++--- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/scripts/make_summary.py b/scripts/make_summary.py index 4425ad5c..43f64b2c 100644 --- a/scripts/make_summary.py +++ b/scripts/make_summary.py @@ -196,7 +196,25 @@ def calculate_costs(n,label,costs): return costs +def calculate_cumulative_cost(): + planning_horizons = snakemake.config['scenario']['planning_horizons'] + cumulative_cost = pd.DataFrame(index = df["costs"].sum().index, + columns=pd.Series(data=np.arange(0,0.1, 0.01), name='social discount rate')) + + #discount cost and express them in money value of planning_horizons[0] + for r in cumulative_cost.columns: + cumulative_cost[r]=[df["costs"].sum()[index]/((1+r)**(index[-1]-planning_horizons[0])) for index in cumulative_cost.index] + + #integrate cost throughout the transition path + for r in cumulative_cost.columns: + for cluster in cumulative_cost.index.get_level_values(level=0).unique(): + for lv in cumulative_cost.index.get_level_values(level=1).unique(): + for sector_opts in cumulative_cost.index.get_level_values(level=2).unique(): + cumulative_cost.loc[(cluster, lv, sector_opts,'cumulative cost'),r] = np.trapz(cumulative_cost.loc[idx[cluster, lv, sector_opts,planning_horizons],r].values, x=planning_horizons) + + return cumulative_cost + def calculate_nodal_capacities(n,label,nodal_capacities): #Beware this also has extraneous locations for country (e.g. biomass) or continent-wide (e.g. fossil gas/oil) stuff for c in n.iterate_components(n.branch_components|n.controllable_one_port_components^{"Load"}): @@ -564,16 +582,17 @@ if __name__ == "__main__": snakemake.config = yaml.safe_load(f) #overwrite some options - snakemake.config["run"] = "test" + snakemake.config["run"] = "version-8" snakemake.config["scenario"]["lv"] = [1.0] - snakemake.config["scenario"]["sector_opts"] = ["Co2L0-168H-T-H-B-I-solar3-dist1"] + snakemake.config["scenario"]["sector_opts"] = ["3H-T-H-B-I-solar3-dist1"] snakemake.config["planning_horizons"] = ['2020', '2030', '2040', '2050'] snakemake.input = Dict() snakemake.input['heat_demand_name'] = 'data/heating/daily_heat_demand.h5' + snakemake.input['costs'] = snakemake.config['costs_dir'] + "costs_{}.csv".format(snakemake.config['scenario']['planning_horizons'][0]) snakemake.output = Dict() for item in outputs: snakemake.output[item] = snakemake.config['summary_dir'] + '/{name}/csvs/{item}.csv'.format(name=snakemake.config['run'],item=item) - + snakemake.output['cumulative_cost'] = snakemake.config['summary_dir'] + '/{name}/csvs/cumulative_cost.csv'.format(name=snakemake.config['run']) networks_dict = {(cluster, lv, opt+sector_opt, planning_horizon) : snakemake.config['results_dir'] + snakemake.config['run'] + '/postnetworks/elec_s{simpl}_{cluster}_lv{lv}_{opt}_{sector_opt}_{planning_horizon}.nc'\ .format(simpl=simpl, @@ -592,6 +611,7 @@ if __name__ == "__main__": print(networks_dict) Nyears = 1 + costs_db = prepare_costs(snakemake.input.costs, snakemake.config['costs']['USD2013_to_EUR2013'], snakemake.config['costs']['discountrate'], @@ -603,3 +623,10 @@ if __name__ == "__main__": df["metrics"].loc["total costs"] = df["costs"].sum() to_csv(df) + + if snakemake.config["foresight"]=='myopic': + cumulative_cost=calculate_cumulative_cost() + cumulative_cost.to_csv(snakemake.config['summary_dir'] + '/' + snakemake.config['run'] + '/csvs/cumulative_cost.csv') + + + \ No newline at end of file From cf71858bc0df2d52e57f4453abf22becde4d53b7 Mon Sep 17 00:00:00 2001 From: martavp Date: Fri, 1 Jan 2021 17:53:59 +0100 Subject: [PATCH 07/10] update documentation and release notes --- doc/installation.rst | 13 ++++++++++--- doc/myopic.rst | 25 ++++++++++++++----------- doc/release_notes.rst | 3 +++ doc/supply_demand.rst | 4 ++-- 4 files changed, 29 insertions(+), 16 deletions(-) diff --git a/doc/installation.rst b/doc/installation.rst index 997f5221..f76a9104 100644 --- a/doc/installation.rst +++ b/doc/installation.rst @@ -16,7 +16,7 @@ its dependencies. Clone the repository: .. code:: bash - projects % git clone git@github.com:PyPSA/pypsa-eur.git + projects % git clone https://github.com/PyPSA/pypsa-eur.git then download and unpack all the PyPSA-Eur data files by running the following snakemake rule: @@ -32,7 +32,7 @@ Next install the technology assumptions database `technology-data xarray version >= 0.15.1, you will need the latest master branch of atlite version 0.0.2. +You can create an enviroment using the environment.yaml file in pypsa-eur/envs: + +.../pypsa-eur % conda env create -f envs/environment.yaml + +.../pypsa-eur % conda activate pypsa-eur + +See details in `PyPSA-Eur Installation `_ Data requirements ================= diff --git a/doc/myopic.rst b/doc/myopic.rst index 6faaaeff..14f651dd 100644 --- a/doc/myopic.rst +++ b/doc/myopic.rst @@ -6,7 +6,7 @@ Myopic transition path The myopic code can be used to investigate progressive changes in a network, for instance, those taking place throughout a transition path. The capacities installed in a certain time step are maintained in the network until their operational lifetime expires. -The myopic approach was initially developed and used in the paper `Early decarbonisation of the European Energy system pays off (2020) `__ but the current implementation complies with the pypsa-eur-sec standard working flow and is compatible with using the higher resolution electricity transmission model `PyPSA-Eur `__ rather than a one-node-per-country model. +The myopic approach was initially developed and used in the paper `Early decarbonisation of the European Energy system pays off (2020) `__ but the current implementation complies with the pypsa-eur-sec standard working flow and is compatible with using the higher resolution electricity transmission model `PyPSA-Eur `__ rather than a one-node-per-country model. The current code applies the myopic approach to generators, storage technologies and links in the power sector and the space and water heating sector. @@ -61,12 +61,15 @@ Wildcards The {planning_horizons} wildcard indicates the timesteps in which the network is optimized, e.g. planning_horizons: [2020, 2030, 2040, 2050] +Options +============= +The total carbon budget for the entire transition path can be indicated in the ``scenario.sector_opts`` in ``config.yaml``. +The carbon budget can be split among the ``planning_horizons`` following an exponential or beta decay. +E.g. ``'cb40ex0'`` splits the a carbon budget equal to 40 GtCO_2 following an exponential decay whose initial linear growth rate $r$ is zero -**{co2_budget_name} wildcard** +$e(t) = e_0 (1+ (r+m)t) e^(-mt)$ -The {co2_budget_name} wildcard indicates the name of the co2 budget. - -A csv file is used as input including the planning_horizons as index, the name of co2_budget as column name, and the maximum co2 emissions (relative to 1990) as values. +See details in Supplementary Note 1 of the paper `Early decarbonisation of the European Energy system pays off (2020) `__ Rules overview ================= @@ -74,17 +77,17 @@ Rules overview General myopic code structure =============================== -The myopic code solves the network for the time steps included in planning_horizons in a recursive loop, so that: +The myopic code solves the network for the time steps included in ``planning_horizons`` in a recursive loop, so that: 1.The existing capacities (those installed before the base year are added as fixed capacities with p_nom=value, p_nom_extendable=False). E.g. for baseyear=2020, capacities installed before 2020 are added. In addition, the network comprises additional generator, storage, and link capacities with p_nom_extendable=True. The non-solved network is saved in ``results/run_name/networks/prenetworks-brownfield``. -The base year is the first element in planning_horizons. Step 1 is implemented with the rule add_baseyear for the base year and with the rule add_brownfield for the remaining planning_horizons. +The base year is the first element in ``planning_horizons``. Step 1 is implemented with the rule add_baseyear for the base year and with the rule add_brownfield for the remaining planning_horizons. -2.The 2020 network is optimized. The solved network is saved in ‘results/run_name/networks/postnetworks’ +2.The 2020 network is optimized. The solved network is saved in ``results/run_name/networks/postnetworks`` 3.For the next planning horizon, e.g. 2030, the capacities from a previous time step are added if they are still in operation (i.e., if they fulfil planning horizon <= commissioned year + lifetime). In addition, the network comprises additional generator, storage, and link capacities with p_nom_extendable=True. The non-solved network is saved in ``results/run_name/networks/prenetworks-brownfield``. -Steps 2 and 3 are solved recursively for all the planning_horizons included in the configuration file. +Steps 2 and 3 are solved recursively for all the planning_horizons included in ``config.yaml``. rule add_existing baseyear @@ -110,8 +113,8 @@ Then, the resulting network is saved in ``results/run_name/networks/prenetworks- rule add_brownfield =================== -The rule add_brownfield loads the network in ‘results/run_name/networks/prenetworks’ and performs the following operation: +The rule add_brownfield loads the network in ``results/run_name/networks/prenetworks`` and performs the following operation: -1.Read the capacities optimized in the previous time step and add them to the network if they are still in operation (i.e., if they fulfil planning horizon < commissioned year + lifetime) +1.Read the capacities optimized in the previous time step and add them to the network if they are still in operation (i.e., if they fulfill planning horizon < commissioned year + lifetime) Then, the resulting network is saved in ``results/run_name/networks/prenetworks_brownfield``. diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 9d9acc4c..66594614 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -2,6 +2,9 @@ Release Notes ########################################## +Future release +=================== +*For the myopic option, a carbon budget and a type of decay (exponential or beta) can be selected in the config file to distribute the budget across the planning_horizons. PyPSA-Eur-Sec 0.4.0 (11th December 2020) ========================================= diff --git a/doc/supply_demand.rst b/doc/supply_demand.rst index 823bdd93..002bd16c 100644 --- a/doc/supply_demand.rst +++ b/doc/supply_demand.rst @@ -106,7 +106,7 @@ Thermal energy storage using hot water tanks Small for decentral applications. -Big pit storage for district heating. +Big water pit storage for district heating. Hydrogen demand @@ -122,7 +122,7 @@ Industry (ammonia, precursor to hydrocarbons for chemicals and iron/steel). Hydrogen supply ================= -SMR, SMR+CCS, electrolysers. +Steam Methane Reforming (SMR), SMR+CCS, electrolysers. Methane demand From 93360119d941f3e9fe8c0484968b139f64dfdb30 Mon Sep 17 00:00:00 2001 From: martavp Date: Mon, 4 Jan 2021 10:07:29 +0100 Subject: [PATCH 08/10] Allow decimals in sector_opts --- Snakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Snakefile b/Snakefile index 0ac8bd99..6e484aae 100644 --- a/Snakefile +++ b/Snakefile @@ -8,7 +8,7 @@ wildcard_constraints: clusters="[0-9]+m?", sectors="[+a-zA-Z0-9]+", opts="[-+a-zA-Z0-9]*", - sector_opts="[-+a-zA-Z0-9]*" + sector_opts="[-+a-zA-Z0-9\.]*" From b98235c1f1ae4d502a688d3982a72ec711900422 Mon Sep 17 00:00:00 2001 From: martavp Date: Mon, 4 Jan 2021 10:13:30 +0100 Subject: [PATCH 09/10] Allow specifying an option to alter the capital cost of carriers by a factor indicated in the config file, eg: solar+c0.5 This is almost a direct copy PyPSA-Eur #167 https://github.com/PyPSA/pypsa-eur/pull/167 A factor altering the maximum capacity (p_nom_max) can also be specified by e.g. solar+p3 One should be careful when using this for solar because the factor is applied to all the generators whose carrier includes the string 'solar' (i.e., it is applied to both utility and rooftop solar) I would suggest implementing 'solar utility' and 'solar rooftop' as carriers, since this can be useful for other selecting processes. Is there is any reason for keeping 'solar' as a carrier for 'solar utility'? The previous way of increasing maximum capacity via the config file (e.g 'solar3') is still present in the code. --- config.default.yaml | 1 + doc/release_notes.rst | 1 + scripts/prepare_sector_network.py | 33 +++++++++++++++++++++++++------ 3 files changed, 29 insertions(+), 6 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 08cda6c4..946977f9 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -24,6 +24,7 @@ scenario: # B for biomass supply, I for industry, shipping and aviation # solarx or onwindx changes the available installable potential by factor x # dist{n} includes distribution grids with investment cost of n times cost in data/costs.csv + # solar+c0.5 reduces the capital cost of solar to 50\% of reference value # for myopic/perfect foresight cb states the carbon budget in GtCO2 (cumulative # emissions throughout the transition path in the timeframe determined by the # planning_horizons), be:beta decay; ex:exponential decay diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 66594614..1e43b92c 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -5,6 +5,7 @@ Release Notes Future release =================== *For the myopic option, a carbon budget and a type of decay (exponential or beta) can be selected in the config file to distribute the budget across the planning_horizons. +*Added an option to alter the capital cost of carriers by a factor via ``carrier+factor`` in the ``{opts}`` wildcard. This can be useful for exploring uncertain cost parameters. Example: ``solar+0.5`` reduces the capital cost of solar to 50\% of original values. PyPSA-Eur-Sec 0.4.0 (11th December 2020) ========================================= diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 088b124f..a875a635 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1984,8 +1984,8 @@ if __name__ == "__main__": from vresutils.snakemake import MockSnakemake snakemake = MockSnakemake( wildcards=dict(network='elec', simpl='', clusters='37', lv='1.0', - opts='', planning_horizons='2020', - sector_opts='120H-T-H-B-I-solar3-dist1-cb40ex0'), + opts='', planning_horizons='2050', + sector_opts='Co2L0-180H-T-H-B-I-solar3-dist1-solar+c0.5-Sabatier+c0.5'), input=dict( network='../pypsa-eur/networks/{network}_s{simpl}_{clusters}_ec_lv{lv}_{opts}.nc', energy_totals_name='resources/energy_totals.csv', co2_totals_name='resources/co2_totals.csv', @@ -2024,7 +2024,7 @@ if __name__ == "__main__": retro_cost_energy = "resources/retro_cost_{network}_s{simpl}_{clusters}.csv", floor_area = "resources/floor_area_{network}_s{simpl}_{clusters}.csv" ), - output=['results/version-8/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{planning_horizons}.nc'] + output=['results/version-0/prenetworks/{network}_s{simpl}_{clusters}_lv{lv}__{sector_opts}_{planning_horizons}.nc'] ) import yaml with open('config.yaml', encoding='utf8') as f: @@ -2150,11 +2150,10 @@ if __name__ == "__main__": print("adding CO2 budget limit as per unit of 1990 levels of",limit) add_co2limit(n, Nyears, limit) - - + for o in opts: for tech in ["solar","onwind","offwind"]: - if tech in o: + if tech in o and "+" not in o: limit = o[o.find(tech)+len(tech):] limit = float(limit.replace("p",".").replace("m","-")) print("changing potential for",tech,"by factor",limit) @@ -2170,6 +2169,28 @@ if __name__ == "__main__": if snakemake.config["sector"]['electricity_distribution_grid']: insert_electricity_distribution_grid(n) + + for o in opts: + if "+" in o: + oo = o.split("+") + carrier_list=np.hstack((n.generators.carrier.unique(), n.links.carrier.unique(), + n.stores.carrier.unique(), n.storage_units.carrier.unique())) + suptechs = map(lambda c: c.split("-", 2)[0], carrier_list) + if oo[0].startswith(tuple(suptechs)): + carrier = oo[0] + # handles only p_nom_max as stores and lines have no potentials + attr_lookup = {"p": "p_nom_max", "c": "capital_cost"} + attr = attr_lookup[oo[1][0]] + factor = float(oo[1][1:]) + if carrier == "AC": # lines do not have carrier + n.lines[attr] *= factor + else: + comps = {"Generator", "Link", "StorageUnit", "Store"} + for c in n.iterate_components(comps): + sel = c.df.carrier.str.contains(carrier) + c.df.loc[sel,attr] *= factor + print("changing", attr ,"for",carrier,"by factor",factor) + if snakemake.config["sector"]['gas_distribution_grid']: insert_gas_distribution_costs(n) if snakemake.config["sector"]['electricity_grid_connection']: From a6a88e26daefeef1938683a37f3e6f70da9fecf0 Mon Sep 17 00:00:00 2001 From: martavp Date: Mon, 4 Jan 2021 16:22:39 +0100 Subject: [PATCH 10/10] Allow white space in sector_opts, e.g. to investigate sensitivity to cost of 'H2 Electrolysis' --- Snakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Snakefile b/Snakefile index 6e484aae..31370603 100644 --- a/Snakefile +++ b/Snakefile @@ -8,7 +8,7 @@ wildcard_constraints: clusters="[0-9]+m?", sectors="[+a-zA-Z0-9]+", opts="[-+a-zA-Z0-9]*", - sector_opts="[-+a-zA-Z0-9\.]*" + sector_opts="[-+a-zA-Z0-9\.\s]*"