diff --git a/Snakefile b/Snakefile index 9722223..13b471a 100644 --- a/Snakefile +++ b/Snakefile @@ -109,6 +109,39 @@ rule plot_electricity_for_heats: ), +rule plot_heat_saving: + params: + clusters=config["plotting"]["clusters"], + planning_horizon=config["plotting"]["planning_horizon"], + output: + figure=RESULTS+"plot_heat_savings_{clusters}_{planning_horizon}.png", + figure_full=RESULTS+"plot_heat_savings_full_year_{clusters}_{planning_horizon}.png", + figure_flexibility=RESULTS+"plot_heat_flexibility_{clusters}_{planning_horizon}.png", + resources: + mem_mb=10000, + script: + "plots/plot_heat_savings.py" + + +rule plot_heat_savings: + input: + expand( + RESULTS + + "plot_heat_savings_{clusters}_{planning_horizon}.png", + **config["plotting"], + ), + expand( + RESULTS + + "plot_heat_savings_full_year_{clusters}_{planning_horizon}.png", + **config["plotting"], + ), + expand( + RESULTS + + "plot_heat_flexibility_{clusters}_{planning_horizon}.png", + **config["plotting"], + ), + + rule plot_electricity_generation: params: clusters=config["plotting"]["clusters"], @@ -366,6 +399,21 @@ rule plot_all: + "plot_electricity_for_heat_{clusters}_{planning_horizon}.png", **config["plotting"], ), + expand( + RESULTS + + "plot_heat_savings_{clusters}_{planning_horizon}.png", + **config["plotting"], + ), + expand( + RESULTS + + "plot_heat_savings_full_year_{clusters}_{planning_horizon}.png", + **config["plotting"], + ), + expand( + RESULTS + + "plot_heat_flexibility_{clusters}_{planning_horizon}.png", + **config["plotting"], + ), expand( RESULTS + "table_heat_pumps_{clusters}.csv", diff --git a/configs/EEE_study/config.retro_tes-industry_2050.yaml b/configs/EEE_study/config.retro_tes-industry_2050.yaml index 2e05d75..0d108e3 100644 --- a/configs/EEE_study/config.retro_tes-industry_2050.yaml +++ b/configs/EEE_study/config.retro_tes-industry_2050.yaml @@ -121,3 +121,4 @@ solving: threads: 8 BarConvTol: 1.e-4 BarHomogeneous: 1 + NumericFocus: 3 diff --git a/plots/_helpers.py b/plots/_helpers.py index e3d92ba..2fc3c1c 100644 --- a/plots/_helpers.py +++ b/plots/_helpers.py @@ -171,6 +171,13 @@ ] ) +# day-ahead prices for 2021 URL:https://www.ffe.de/en/publications/european-day-ahead-electricity-prices-in-2022/ +HISTORIC_PRICES = {"AT": 106.9, "BE": 104.1, "CH": 114.9, "CZ": 100.7, "DE": 96.8, + "DK": 88.1, "EE": 86.7, "ES": 111.9, "FI": 72.3, "FR": 109.2, + "GR": 116.4, "HU": 113.9, "IT": 125.2, "LT": 90.5, "LV": 88.8, + "NL": 103.0, "NO": 74.7, "PL": 87.0, "PT": 112.0, "RO": 110.9, + "RS": 114.0, "SE": 66.0, "SI": 115.0, "SK": 102.8} + def rename_techs(label): diff --git a/plots/plot_capacity_expansion.py b/plots/plot_capacity_expansion.py index 03002db..6975513 100644 --- a/plots/plot_capacity_expansion.py +++ b/plots/plot_capacity_expansion.py @@ -39,19 +39,29 @@ def plot_capacities(caps_df, clusters, planning_horizon, plot_width=7): df = df[df.index != 'solid biomass transport'] # convert to GW - df = df / 1e3 + df = df / 1e6 df = df.groupby(df.index.map(rename_techs)).sum() - caps_threshold = 10 + caps_threshold = 10e-3 to_drop = df.index[df.max(axis=1) < caps_threshold] #df < logger.info( - f"Dropping technology with capacity below {caps_threshold} GW" + f"Dropping technology with capacity below {caps_threshold} TW" ) logger.debug(df.loc[to_drop]) df = df.drop(to_drop) + # remove extra technologies that are not expanded + drop_techs = ["gas boiler", "process emissions", + "solid biomass", "solid biomass CHP", + "gas storage", + "kerosene for aviation", "land transport oil", "naphtha for industry", + "shipping methanol", "shipping oil", + "coal", "lignite", "oil", + "coal for industry", "gas for industry"] + df = df.drop(set(df.index).intersection(set(drop_techs))) + logger.info(f"Total optimal capacity is {round(df.sum())} GW") new_index = PREFERRED_ORDER.intersection(df.index).append( @@ -75,11 +85,11 @@ def plot_capacities(caps_df, clusters, planning_horizon, plot_width=7): plt.xticks(rotation=0, fontsize=10) - ax.set_ylabel("Capacity expansion [GW]") + ax.set_ylabel("Capacity expansion [TW]") ax.set_xlabel("") - ax.set_ylim([0,16000]) - ax.set_yticks(np.arange(0, 17000, 2000)) + ax.set_ylim([0,16]) + ax.set_yticks(np.arange(0, 17, 2)) x_ticks = list(df.columns) if planning_horizon in ["2040", "2050"] and "Limited \nRenovation &\nCost-Optimal Heating" in x_ticks: diff --git a/plots/plot_co2_level.py b/plots/plot_co2_level.py index b4f6495..d1155e2 100644 --- a/plots/plot_co2_level.py +++ b/plots/plot_co2_level.py @@ -110,18 +110,18 @@ "solid biomass for industry", "gas for industry", "coal for industry", - "methanol", - "oil", - "lignite", - "coal", "shipping oil", "shipping methanol", "naphtha for industry", "land transport oil", "kerosene for aviation", - - "transmission lines", - "distribution lines", + "process emissions", + "SMR", + "coal", + "lignite", + "methanol", + "oil", + "gas pipeline", "H2 pipeline", @@ -132,26 +132,17 @@ "methanation", "BEV charger", "V2G", - "SMR", "methanolisation", "battery storage", "gas storage", "H2 storage", "TES", - - "hydroelectricity", "OCGT", "CCGT", - "onshore wind", - "offshore wind", - "solar PV", - "solar thermal", - "solar rooftop", "co2", "CO2 sequestration", - "process emissions", "gas CHP", "resistive heater", @@ -195,12 +186,14 @@ def get_co2_balance(n, nice_name): def plot_co2_balance(co2_df, clusters, planning_horizon, plot_width=7): + # convert to Billion tCO2_eq + co2_df = co2_df / 1e9 # filter out technologies with very small emission max_emissions = co2_df.abs().sum().max() / 2 co2_threshold = max_emissions / 100 # 1% of max to_drop = co2_df.index[co2_df.abs().max(axis=1) < co2_threshold] logger.info( - f"Dropping technology with co2 balance below {co2_threshold} ton CO2_eq per year" + f"Dropping technology with co2 balance below {co2_threshold} Bton CO2_eq per year" ) logger.debug(co2_df.loc[to_drop]) co2_df = co2_df.drop(to_drop) @@ -229,10 +222,11 @@ def plot_co2_balance(co2_df, clusters, planning_horizon, plot_width=7): labels = labels[:-neg_co2_length] + labels[-neg_co2_length:][::-1] plt.xticks(rotation=0, fontsize=10) - ax.set_ylabel("CO$_2$ emissions [tCO$_{2-eq}$]") + ax.set_ylabel("CO$_2$ emissions [BtCO$_{2-eq}$]") ax.set_xlabel("") - ax.set_ylim([-4e9,4e9]) - ax.set_yticks(np.arange(-4.0e9, 4.0e9, 5e8)) + ax.set_yticks(np.arange(-4.0, 4.0, 0.5)) + ax.set_ylim([-3.6,3.6]) + x_ticks = list(co2_df.columns) if planning_horizon in ["2040", "2050"] and "Limited \nRenovation &\nOptimal Heating" in x_ticks: # replace name for Limited Renovation scenario for 2030 to be LROH @@ -372,7 +366,7 @@ def fill_table_df(df, planning_horizon, scenarios, values): co2_savings_df = co2_emission_BAU - table_emissions_df co2_savings_df.index = ["CO2 emission savings [tCO2_eq]"] # saving in cubic meter of gas (1 m3 natural gas = 1,9 kg CO2) Source: https://www.eeagrants.gov.pt/media/2776/conversion-guidelines.pdf - co2_savings_df.loc["Equivalent gas savings [m3]"] = co2_savings_df.loc["CO2 emission savings [tCO2_eq]", :] * 1000 / 1.9 + co2_savings_df.loc["Equivalent gas savings [trillion m3]"] = co2_savings_df.loc["CO2 emission savings [tCO2_eq]", :] * 1000 / 1.9 / 1e12 # move to base directory diff --git a/plots/plot_electricity_bills.py b/plots/plot_electricity_bills.py index 9bbb054..a944e60 100644 --- a/plots/plot_electricity_bills.py +++ b/plots/plot_electricity_bills.py @@ -14,7 +14,7 @@ warnings.filterwarnings("ignore") from _helpers import mock_snakemake, update_config_from_wildcards, load_network, \ change_path_to_pypsa_eur, change_path_to_base, \ - CO2L_LIMITS, LINE_LIMITS, BAU_HORIZON + CO2L_LIMITS, LINE_LIMITS, BAU_HORIZON, HISTORIC_PRICES def get_households(): @@ -123,6 +123,10 @@ def electricity_bills(network, households): total_cost_country += total_gas_cost_country # electricity bill per household [EUR/household] (households given in thousands) elec_bills_household = total_cost_country / (households*1e3) + + # add EU average bills + average_bills_household = total_cost_country.sum() / (households.sum() * 1e3) + elec_bills_household["EUR"] = average_bills_household return elec_bills_household @@ -141,7 +145,14 @@ def plot_electricity_cost(df_prices, name): "Optimal Renovation and Electric Heating":"limegreen", "Limited Renovation and Cost-Optimal Heating":"royalblue", "No Renovation and Electric Heating":"#f4b609", - "BAU": "grey"} + "BAU": "grey", + "Historic data": '#900C3F'} + + # add historic electricity price for BAU plot + if "BAU" in df_prices.index and name == "prices": + historic_prices = pd.Series(HISTORIC_PRICES).to_frame().T + historic_prices.index = ["Historic data"] + sorted_df_prices = pd.concat([sorted_df_prices, historic_prices], axis=0) # plot as bar plot fig, ax = plt.subplots(figsize=(7,3)) @@ -243,9 +254,9 @@ def electricity_prices(network, households): # electricity bill per MWh [EUR/MWh] energy_price_MWh = prices / total_load_country - # drop EU - if "EU" in energy_price_MWh.index: - energy_price_MWh.drop("EU", axis=0, inplace=True) + # add EU average + average_price_MWh = total_cost.sum() / total_load_country.sum() + energy_price_MWh["EUR"] = average_price_MWh return energy_price_MWh @@ -256,7 +267,7 @@ def electricity_prices(network, households): snakemake = mock_snakemake( "plot_electricity_bill", clusters="48", - planning_horizon="2030", + planning_horizon="2020", ) # update config based on wildcards config = update_config_from_wildcards(snakemake.config, snakemake.wildcards) diff --git a/plots/plot_electricity_for_heat.py b/plots/plot_electricity_for_heat.py index 7839850..f6a11bc 100644 --- a/plots/plot_electricity_for_heat.py +++ b/plots/plot_electricity_for_heat.py @@ -66,10 +66,10 @@ def get_elec_consumption_for_heat(n): return elec_demand_f_heating -def plot_elec_consumption_for_heat(dict_elec): +def plot_elec_consumption_for_heat(dict_elec, horizon): # set heights for each subplots if "BAU" in dict_elec.keys(): - heights = [1.4] + heights = [1.5] else: heights = [1.4] * 4 fig = plt.figure(figsize=(6.4, sum(heights))) @@ -85,7 +85,7 @@ def plot_elec_consumption_for_heat(dict_elec): ax.set_facecolor("whitesmoke") print("Hard coded coordinates on x-axis, selected for 3H granularity") - where = [359, 527] + where = [359, 695] elec_demand_f_heating = elec_demand_f_heating.reset_index()[["ground heat pump", "air heat pump", "resistive heater", "gas boiler"]] colors = ["#2fb537", "#48f74f", "#d8f9b8", "#db6a25"] @@ -95,20 +95,20 @@ def plot_elec_consumption_for_heat(dict_elec): linewidth=0, ax=ax ) + cumulative = (elec_demand_f_heating / 1e3).cumsum(axis=1) + ax.plot(cumulative.iloc[where[0]:where[1]], color='black', linewidth=0.5) + ax.set_xlabel("", fontsize=12) - ax.set_ylim([1,2000]) - ax.set_yscale("log") # Set the y-axis to logarithmic scale + # get in y limits for all horizons + y_limits = {'2020':1600, '2030':1100, '2040':900, '2050': 900} - # Custom y-ticks for logarithmic scale - ax.yaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=10)) - ax.yaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs='auto')) - ax.yaxis.set_minor_formatter(ticker.NullFormatter()) + ax.set_ylim([1, y_limits[horizon]]) if i < 3: ax.set_xticks([]) ax.set_xlabel("") if i == 3 or "BAU" in dict_elec.keys(): - ticks = [i for i in range(where[0], where[1], 24)] + ticks = [i for i in range(where[0], where[1], 48)] ax.set_xticks(ticks) # Set the tick positions ticks = [ str(n.snapshots[i]).split(" ")[0].split("-")[2]+"."+str(n.snapshots[i]).split(" ")[0].split("-")[1] @@ -136,10 +136,10 @@ def plot_elec_consumption_for_heat(dict_elec): loc=[1.02, -.2], fontsize=10 ) if len(dict_elec.keys()) == 1: - ylabel = "Electricity and Gas \nConsumption to Cover \nHeating Demands [GW]" + ylabel = "Electricity and Gas \nConsumption to Cover \nHeating Demands [$GW_{el}$]" axes[0].set_ylabel(ylabel, fontsize=10) else: - ylabel = "Electricity and Gas Consumption to \n Cover Heating Demands [GW]" + ylabel = "Electricity and Gas Consumption to \n Cover Heating Demands [$GW_{el}$]" axes[2].set_ylabel(ylabel, fontsize=10) axes[2].yaxis.set_label_coords(-0.1, 1) plt.subplots_adjust(hspace=0.3) @@ -203,5 +203,5 @@ def plot_elec_consumption_for_heat(dict_elec): total_elec_consumption_for_heat[name] = elec_consumption_for_heat # plot and store electricity prices - plot_elec_consumption_for_heat(total_elec_consumption_for_heat) + plot_elec_consumption_for_heat(total_elec_consumption_for_heat, planning_horizon) diff --git a/plots/plot_heat_savings.py b/plots/plot_heat_savings.py new file mode 100644 index 0000000..d416cd7 --- /dev/null +++ b/plots/plot_heat_savings.py @@ -0,0 +1,280 @@ +# -*- coding: utf-8 -*- +# SPDX-FileCopyrightText: Open Energy Transition gGmbH +# +# SPDX-License-Identifier: AGPL-3.0-or-later + +import os +import sys +sys.path.append("../submodules/pypsa-eur") +import pypsa +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import matplotlib.ticker as ticker +import numpy as np +import pandas as pd +import warnings +warnings.filterwarnings("ignore") +from _helpers import mock_snakemake, update_config_from_wildcards, load_network, \ + change_path_to_pypsa_eur, change_path_to_base, \ + LINE_LIMITS, CO2L_LIMITS, BAU_HORIZON + + +def get_heat_data(n, with_flexibility=False): + # get heat load + heat_load = n.loads_t.p_set.filter(like="heat").multiply(n.snapshot_weightings.objective, axis=0).sum(axis=1) + # get heat savings by retrofitting + heat_savings = n.generators_t.p.filter(like="retrofitting").multiply(n.snapshot_weightings.objective, axis=0).sum(axis=1) + # get heat flexibility + if with_flexibility: + heat_flex_carr = ['residential rural heat', 'residential urban decentral heat', 'urban central heat'] + heat_flexibility_techs = n.stores.query("carrier in @heat_flex_carr").index + heat_flexibility = n.stores_t.p[heat_flexibility_techs].multiply(n.snapshot_weightings.objective, axis=0).sum(axis=1) + + # compute heat load after renovation and flexbility + net_heat_load = heat_load - heat_savings - heat_flexibility + # concatenate net heat load and heat_savings + heat_data = pd.concat([net_heat_load, heat_flexibility + heat_savings], axis=1) + heat_data.columns = ["Net heat demand", "Heat savings by renovation"] + else: + # compute heat load after renovation + net_heat_load = heat_load - heat_savings + # concatenate net heat load and heat_savings + heat_data = pd.concat([net_heat_load, heat_savings], axis=1) + heat_data.columns = ["Net heat demand", "Heat savings by renovation"] + return heat_data + + +def get_flexibility(network): + # heat flexibility + heat_flex_carr = ['residential rural heat', 'residential urban decentral heat', 'urban central heat'] + heat_flexibility_techs = n.stores.query("carrier in @heat_flex_carr").index + heat_flexibility = n.stores_t.p[heat_flexibility_techs].multiply(n.snapshot_weightings.objective, axis=0).sum(axis=1) + # EV flexibility + ev_flex_carr = ['Li ion'] + ev_flexibility_techs = n.stores.query("carrier in @ev_flex_carr").index + ev_flexibility = n.stores_t.p[ev_flexibility_techs].multiply(n.snapshot_weightings.objective, axis=0).sum(axis=1) + # concatenate flexibility + flexibility = pd.concat([ev_flexibility, heat_flexibility], axis=1) + flexibility.columns = ["EV", "heat"] + return flexibility + +def plot_elec_consumption_for_heat(dict_elec, full_year=False): + # set heights for each subplots + if "BAU" in dict_elec.keys(): + heights = [1.3] + else: + heights = [1.4] * 3 + fig = plt.figure(figsize=(6.4, sum(heights))) + gs = gridspec.GridSpec(len(heights), 1, height_ratios=heights) + axes = [fig.add_subplot(gs[i]) for i in range(len(heights))] + if "BAU" in dict_elec.keys(): + axes = axes * 3 + i=0 + for name, heat_demand in dict_elec.items(): + ax = axes[i] + ax.set_facecolor("whitesmoke") + + print("Hard coded coordinates on x-axis, selected for 3H granularity") + if not full_year: + where = [359, 695] + else: + where = [0, 365] + heat_demand = heat_demand.resample('24H').max() + + heat_demand = heat_demand.reset_index(drop=True) + colors = ["#a63f3f", "#e39191"] + + (heat_demand/1e3).iloc[where[0] : where[1]].plot( + kind='area', stacked=True, color=colors, legend=False, + linewidth=0, ax=ax + ) + + cumulative = (heat_demand / 1e3).cumsum(axis=1) + ax.plot(cumulative.iloc[where[0]:where[1]], color='black', linewidth=0.5) + + ax.set_xlabel("", fontsize=12) + ax.set_ylim([1,1700]) + + if i < 2: + ax.set_xticks([]) + ax.set_xlabel("") + if i == 2 or "BAU" in dict_elec.keys(): + if not full_year: + ticks = [i for i in range(where[0], where[1], 48)] + ax.set_xticks(ticks) # Set the tick positions + ticks = [ + str(n.snapshots[i]).split(" ")[0].split("-")[2]+"."+str(n.snapshots[i]).split(" ")[0].split("-")[1] + for i in ticks + ] + ax.set_xticklabels(ticks, fontsize=10) # Set the tick labels + else: + snapshots_series = pd.Series(1, index=n.snapshots) + first_day_each_month = snapshots_series.resample('MS').first().index + days_diff = first_day_each_month.to_series().diff().dt.days + days_diff.iloc[0] = 0 + cumulative_days = days_diff.cumsum() + cumulative_days_list = cumulative_days.tolist() + ax.set_xticks([int(x) for x in cumulative_days_list]) + ticks = [str(i).split(" ")[0].replace("-",".")[5:] for i in first_day_each_month] + ax.set_xticklabels(ticks, fontsize=10, rotation=15) # Set the tick labels + plt.xticks(fontsize=10) + plt.yticks(fontsize=10) + + ax.set_xlim([where[0], where[1]-1]) + # change name to LR for 2040 and 2050 + name = "Limited Renovation and Electric Heating" if name == "Limited Renovation and Cost-Optimal Heating" and planning_horizon in ["2040", "2050"] else name + ax.set_title(name, fontsize=10) + i+= 1 + + handles1, _ = axes[0].get_legend_handles_labels() + axes[0].legend( + reversed(handles1[0:7]), + [ + "Heat savings by renovation", + "Net heat demand", + ], + loc=[1.02, -.2], fontsize=10 + ) + if len(dict_elec.keys()) == 1: + ylabel = "Heat demand [$GW_{th}$]" + axes[0].set_ylabel(ylabel, fontsize=10) + else: + ylabel = "Heat demand [$GW_{th}$]" + axes[1].set_ylabel(ylabel, fontsize=10) + # axes[1].yaxis.set_label_coords(-0.1, 1) + plt.subplots_adjust(hspace=0.3) + if not full_year: + plt.savefig(snakemake.output.figure, bbox_inches='tight', dpi=600) + else: + plt.savefig(snakemake.output.figure_full, bbox_inches='tight', dpi=600) + + +def plot_flexibility(dict_elec): + # set heights for each subplots + if "BAU" in dict_elec.keys(): + heights = [1.2] + else: + heights = [1.4] * 3 + fig = plt.figure(figsize=(6.4, sum(heights))) + gs = gridspec.GridSpec(len(heights), 1, height_ratios=heights) + axes = [fig.add_subplot(gs[i]) for i in range(len(heights))] + if "BAU" in dict_elec.keys(): + axes = axes * 3 + i=0 + colors = ["red", "black"] + for name, heat_demand in dict_elec.items(): + ax = axes[i] + ax.set_facecolor("whitesmoke") + + print("Hard coded coordinates on x-axis, selected for 3H granularity") + where = [0, 8759] + heat_demand = heat_demand.reset_index(drop=True) + + ax.plot((heat_demand["EV"]/1e3).iloc[where[0]:where[1]], color="red", linewidth=0.5, label="EV flexibility") + ax.plot((heat_demand["heat"]/1e3).iloc[where[0]:where[1]], color="black", linewidth=0.5, label="Heat flexibility") + + ax.set_xlabel("", fontsize=12) + ax.set_ylim([-550,550]) + ax.set_yticks(np.arange(-500, 600, 250)) + + if i < 2: + ax.set_xticks([]) + ax.set_xlabel("") + if i == 2 or "BAU" in dict_elec.keys(): + snapshots_series = pd.Series(1, index=n.snapshots) + first_day_each_month = snapshots_series.resample('MS').first().index + hours_diff = first_day_each_month.to_series().diff().dt.total_seconds() / 3600 + hours_diff.iloc[0] = 0 + cumulative_hours = hours_diff.cumsum() + cumulative_hours_list = cumulative_hours.tolist() + ax.set_xticks([int(x) for x in cumulative_hours_list]) + ticks = [str(i).split(" ")[0].replace("-",".")[5:] for i in first_day_each_month] + ax.set_xticklabels(ticks, fontsize=10, rotation=15) # Set the tick labels + plt.xticks(fontsize=10) + plt.yticks(fontsize=10) + + ax.set_xlim([where[0], where[1]-1]) + # change name to LR for 2040 and 2050 + name = "Limited Renovation and Electric Heating" if name == "Limited Renovation and Cost-Optimal Heating" and planning_horizon in ["2040", "2050"] else name + ax.set_title(name, fontsize=10) + i+= 1 + + axes[0].legend(loc=[1.02, -.2], fontsize=10) + if len(dict_elec.keys()) == 1: + ylabel = "Flexibility [GW]" + axes[0].set_ylabel(ylabel, fontsize=10) + else: + ylabel = "Flexibility [GW]" + axes[1].set_ylabel(ylabel, fontsize=10) + plt.subplots_adjust(hspace=0.3) + plt.savefig(snakemake.output.figure_flexibility, bbox_inches='tight', dpi=600) + + + +if __name__ == "__main__": + if "snakemake" not in globals(): + snakemake = mock_snakemake( + "plot_heat_saving", + clusters="48", + planning_horizon="2030", + ) + # update config based on wildcards + config = update_config_from_wildcards(snakemake.config, snakemake.wildcards) + + # network parameters + co2l_limits = CO2L_LIMITS + line_limits = LINE_LIMITS + clusters = config["plotting"]["clusters"] + planning_horizon = config["plotting"]["planning_horizon"] + opts = config["plotting"]["sector_opts"] + lineex = line_limits[planning_horizon] + sector_opts = f"Co2L{co2l_limits[planning_horizon]}-{opts}" + + # move to submodules/pypsa-eur + change_path_to_pypsa_eur() + + # define scenario namings + if planning_horizon == BAU_HORIZON: + scenarios = {"BAU": "BAU"} + else: + scenarios = {"flexible": "Optimal Renovation and Cost-Optimal Heating", + "retro_tes": "Optimal Renovation and Electric Heating", + "flexible-moderate": "Limited Renovation and Cost-Optimal Heating"} #, + #"rigid": "No Renovation and Electric Heating"} + + # load networks + networks = {} + for scenario, nice_name in scenarios.items(): + n = load_network(lineex, clusters, sector_opts, planning_horizon, scenario) + if n is None: + # Skip further computation for this scenario if network is not loaded + print(f"Network is not found for scenario '{scenario}', planning year '{planning_horizon}'. Skipping...") + continue + + networks[nice_name] = n + + # move to base directory + change_path_to_base() + + total_heat_data = {} + total_heat_data_with_flex = {} + total_flexibility = {} + for name, network in networks.items(): + if network is None: + # Skip further computation for this scenario if network is not loaded + print(f"Network is not found for scenario '{scenario}', planning year '{planning_horizon}'. Skipping...") + continue + heat_data= get_heat_data(network) + total_heat_data[name] = heat_data + heat_data_with_flex = get_heat_data(network, with_flexibility=True) + total_heat_data_with_flex[name] = heat_data_with_flex + # get flexibility data + flexibility = get_flexibility(network) + total_flexibility[name] = flexibility + + # plot heat demand data for short period + plot_elec_consumption_for_heat(total_heat_data_with_flex, full_year=False) + # plot heat demand data for full year + plot_elec_consumption_for_heat(total_heat_data, full_year=True) + # plot heat flexibility + plot_flexibility(total_flexibility) \ No newline at end of file diff --git a/plots/plot_total_costs.py b/plots/plot_total_costs.py index 5e7244c..8940789 100644 --- a/plots/plot_total_costs.py +++ b/plots/plot_total_costs.py @@ -156,11 +156,11 @@ def plot_capacities(caps_df, clusters, planning_horizon, plot_width=7): # drop solid biomass transport df = df[df.index != 'solid biomass transport'] - # convert to GW - df = df / 1e3 + # convert to TW + df = df / 1e6 df = df.groupby(df.index.map(rename_techs)).sum() - caps_threshold = 20 + caps_threshold = 20e-3 to_drop = df.index[df.max(axis=1) < caps_threshold] #df < logger.info( @@ -170,7 +170,7 @@ def plot_capacities(caps_df, clusters, planning_horizon, plot_width=7): df = df.drop(to_drop) - logger.info(f"Total optimal capacity is {round(df.sum())} GW") + logger.info(f"Total optimal capacity is {round(df.sum())} TW") new_index = PREFERRED_ORDER.intersection(df.index).append( df.index.difference(PREFERRED_ORDER) @@ -193,11 +193,11 @@ def plot_capacities(caps_df, clusters, planning_horizon, plot_width=7): plt.xticks(rotation=0, fontsize=10) - ax.set_ylabel("Installed capacities [GW]") + ax.set_ylabel("Installed capacities [TW]") ax.set_xlabel("") - ax.set_ylim([0,22000]) - ax.set_yticks(np.arange(0, 23000, 2000)) + ax.set_ylim([0,22]) + ax.set_yticks(np.arange(0, 23, 2)) x_ticks = list(df.columns) if planning_horizon in ["2040", "2050"] and "Limited \nRenovation &\nCost-Optimal Heating" in x_ticks: diff --git a/plots/table_heat_pumps.py b/plots/table_heat_pumps.py index 42163fa..b01d225 100644 --- a/plots/table_heat_pumps.py +++ b/plots/table_heat_pumps.py @@ -19,9 +19,10 @@ def define_heat_pump_dataframe(): # Define index levels - index_level_0 = ['Optimal Heat Pump Capacity [MW_el]', - 'Optimal Heat Pump Capacity [MW_el]', - 'Optimal Heat Pump Capacity [MW_el]', + index_level_0 = ['Optimal Heat Pump Capacity [GW_el]', + 'Optimal Heat Pump Capacity [GW_el]', + 'Optimal Heat Pump Capacity [GW_el]', + 'Optimal Heat Pump Capacity [GW_el]', 'Number of heat pumps [millions]', 'Number of heat pumps [millions]', 'Number of heat pumps [millions]', @@ -29,6 +30,7 @@ def define_heat_pump_dataframe(): index_level_1 = ['Residential decentral', 'Services', 'Central', + 'Total', 'Residential decentral', 'Services', 'Central', @@ -110,11 +112,12 @@ def define_heat_pump_dataframe(): p_nom_opt = n.links.filter(like=h, axis=0).filter(like="heat pump", axis=0).p_nom_opt.sum() cop = n.links_t.efficiency.filter(like=h).filter(like="heat pump").mean().mean() - df_heat_pumps.loc[("Optimal Heat Pump Capacity [MW_el]", h_name), (planning_horizon, nice_name)] = p_nom_opt.round(2) + df_heat_pumps.loc[("Optimal Heat Pump Capacity [GW_el]", h_name), (planning_horizon, nice_name)] = (p_nom_opt / 1e3).round(2) df_heat_pumps.loc[("Number of heat pumps [millions]", h_name), (planning_horizon, nice_name)] = (p_nom_opt * cop / (hp_capacity[h] * 1e3)).round(2) # estimate heat pump df_heat_pumps.loc[("Number of heat pumps [millions]", "Total"), (planning_horizon, nice_name)] = df_heat_pumps.loc[("Number of heat pumps [millions]", heat_pumps.values()), (planning_horizon, nice_name)].sum() + df_heat_pumps.loc[("Optimal Heat Pump Capacity [GW_el]", "Total"), (planning_horizon, nice_name)] = df_heat_pumps.loc[("Optimal Heat Pump Capacity [GW_el]", heat_pumps.values()), (planning_horizon, nice_name)].sum() # set heat plan amount based on EU action plan df_heat_pumps.loc[('Number of heat pumps [millions]', 'Total'), ("2030", "EU action plan (Announced Pledges Scenario) [2]")] = 45 diff --git a/plots/table_infra_savings.py b/plots/table_infra_savings.py index aa8eb94..623632b 100644 --- a/plots/table_infra_savings.py +++ b/plots/table_infra_savings.py @@ -41,7 +41,8 @@ # define scenario namings scenarios = {"flexible": "Optimal Renovation and Cost-Optimal Heating", "retro_tes": "Optimal Renovation and Electric Heating", - "flexible-moderate": "Limited Renovation and Cost-Optimal Heating", + "flexible-moderate": "Limited Renovation and Cost-Optimal Heating", + "rigid": "No Renovation and Electric Heating" } # define dataframe to store infra savings @@ -61,8 +62,8 @@ ("2050", "wind"), ("2050", "solar"), ("2050", "gas") ] ) - df_savings.columns = pd.MultiIndex.from_tuples(df_savings.columns, names=['horizon','tech']) - cost_savings.columns = pd.MultiIndex.from_tuples(cost_savings.columns, names=['horizon','tech']) + df_savings.columns = pd.MultiIndex.from_tuples(df_savings.columns, names=['horizon','Installed capacity [GW]']) + cost_savings.columns = pd.MultiIndex.from_tuples(cost_savings.columns, names=['horizon','Capital cost [BEur]']) for planning_horizon in planning_horizons: lineex = line_limits[planning_horizon] @@ -82,41 +83,23 @@ # estimate upper and lower limits of congestion of grid solar_carriers = ["solar", "solar rooftop"] - solar = ( - n.generators.query("carrier in @solar_carriers").p_nom_opt.sum() - - b.generators.query("carrier in @solar_carriers").p_nom_opt.sum() - )/1e3 + solar = n.generators.query("carrier in @solar_carriers").p_nom_opt.sum() / 1e3 wind_carriers = ["onwind", "offwind-ac", "offwind-dc"] - wind = ( - n.generators.query("carrier in @wind_carriers").p_nom_opt.sum() - - b.generators.query("carrier in @wind_carriers").p_nom_opt.sum() - )/1e3 + wind = n.generators.query("carrier in @wind_carriers").p_nom_opt.sum() / 1e3 CCGT_carriers = ["CCGT"] - gas = ( - n.links.query("carrier in @CCGT_carriers").p_nom_opt.multiply(n.links.efficiency).sum() - - b.links.query("carrier in @CCGT_carriers").p_nom_opt.multiply(b.links.efficiency).sum() - )/1e3 + gas = n.links.query("carrier in @CCGT_carriers").p_nom_opt.multiply(n.links.efficiency).sum() / 1e3 - df_savings.loc[nice_name, (planning_horizon, "solar")] = solar - df_savings.loc[nice_name, (planning_horizon, "wind")] = wind - df_savings.loc[nice_name, (planning_horizon, "gas")] = gas + df_savings.loc[nice_name, (planning_horizon, "solar")] = solar.round(2) + df_savings.loc[nice_name, (planning_horizon, "wind")] = wind.round(2) + df_savings.loc[nice_name, (planning_horizon, "gas")] = gas.round(2) cap_costs = compute_costs(n, nice_name, "Capital") wind_costs_carriers = ["Generator:Offshore Wind (AC)", "Generator:Offshore Wind (DC)", "Generator:Onshore Wind"] - cost_savings.loc[nice_name, (planning_horizon, "wind")] = ( - b_costs.loc[wind_costs_carriers].sum()[0]- - cap_costs.loc[wind_costs_carriers].sum()[0] - )/1e9 + cost_savings.loc[nice_name, (planning_horizon, "wind")] = (cap_costs.loc[wind_costs_carriers].sum()[0] / 1e9).round(2) solar_costs_carriers = ["Generator:Solar", "Generator:solar rooftop"] - cost_savings.loc[nice_name, (planning_horizon, "solar")] = ( - b_costs.loc[solar_costs_carriers].sum()[0]- - cap_costs.loc[solar_costs_carriers].sum()[0] - )/1e9 + cost_savings.loc[nice_name, (planning_horizon, "solar")] = (cap_costs.loc[solar_costs_carriers].sum()[0] / 1e9).round(2) gas_costs_carriers = ["Store:gas", "Link:Open-Cycle Gas"] - cost_savings.loc[nice_name, (planning_horizon, "gas")] = ( - b_costs.loc[gas_costs_carriers].sum()[0]- - cap_costs.loc[gas_costs_carriers].sum()[0] - )/1e9 + cost_savings.loc[nice_name, (planning_horizon, "gas")] = (cap_costs.loc[gas_costs_carriers].sum()[0] / 1e9).round(2) # add name for columns df_savings.index.name = "Scenario [GW]"