Merge pull request #57 from open-energy-transition/DSR

table_DSR

Snakefile

@@ -68,6 +68,7 @@ rule plot_electricity_bill:
     output:
         figure_bills=RESULTS+"plot_bill_per_household_{clusters}_{planning_horizon}.png",
         figure_price=RESULTS+"plot_prices_per_MWh_{clusters}_{planning_horizon}.png",
+        figure_opex=RESULTS+"plot_opex_industry_{clusters}_{planning_horizon}.png",
     resources:
         mem_mb=20000,
     script:
@@ -263,6 +264,18 @@ rule plot_co2_level:
         "plots/plot_co2_level.py"
 
 
+rule plot_DSR:
+    params:
+        clusters=config["plotting"]["clusters"],
+        planning_horizon=config["plotting"]["planning_horizon"],
+    output:
+        table=RESULTS+"table_DSR_{clusters}.csv",
+    resources:
+        mem_mb=20000,
+    script:
+        "plots/table_DSR.py"
+
+
 rule plot_co2_levels:
     input:
         expand(
@@ -308,6 +321,7 @@ rule get_infra_saving:
     output:
         table_cap=RESULTS+"table_infra_savings_caps_{clusters}.csv",
         table_costs=RESULTS+"table_infra_savings_costs_{clusters}.csv",
+        table_land=RESULTS+"table_infra_savings_land_{clusters}.csv",
     resources:
         mem_mb=20000,
     script:
@@ -326,6 +340,11 @@ rule get_infra_savings:
             + "table_infra_savings_costs_{clusters}.csv",
             **config["plotting"],
         ),
+        expand(
+            RESULTS
+            + "table_infra_savings_land_{clusters}.csv",
+            **config["plotting"],
+        ),
 
 
 rule plot_transmission_congestion:
@@ -484,6 +503,11 @@ rule plot_all:
             + "table_infra_savings_costs_{clusters}.csv",
             **config["plotting"],
         ),
+        expand(
+            RESULTS
+            + "table_infra_savings_land_{clusters}.csv",
+            **config["plotting"],
+        ),
         expand(
             RESULTS
             + "table_res_share_{clusters}.csv",

plots/_helpers.py

@@ -178,6 +178,11 @@
                    "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}
 
+# land usage for 1 MW nameplate capacity of wind (34.5 ha = 0.345 km^2) URL: https://www.nrel.gov/docs/fy09osti/45834.pdf
+LAND_FOR_WIND = 0.345
+
+# land usage for 1 MW nameplate capacity of solar (4-5 acres = 0.02 km^2) URL: https://amplussolar.com/blogs/1mw-solar-power-plant
+LAND_FOR_SOLAR = 0.02
 
 def rename_techs(label):
 

plots/plot_electricity_bills.py

@@ -3,12 +3,9 @@
 #
 # 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 numpy as np
 import pandas as pd
 import warnings
 warnings.filterwarnings("ignore")
@@ -184,7 +181,36 @@ def plot_electricity_cost(df_prices, name):
         plt.savefig(snakemake.output.figure_price, bbox_inches='tight', dpi=600)
 
 
-def electricity_prices(network, households):
+def plot_industry_opex(df):
+
+    # color codes f
+    color_codes = {
+        "Coal":"#545454",
+        "Methanol":"#468c8b",
+        "Biomass":"#baa741",
+        "Hydrogen":"#f073da",
+        "Oil products": "#aaaaaa",
+        "Gas": "#e05b09",
+        "Heat": "#cc1f1f",
+        "Electricity": "#110d63",
+    }
+
+    fig, ax = plt.subplots(figsize=(3,3))
+    df.T.plot.bar(ax=ax, width=0.7, color=color_codes, stacked=True)
+    ax.set_facecolor("white")
+    handles, labels = ax.get_legend_handles_labels()
+    if planning_horizon in ["2040", "2050"]:
+        labels = ["LREH" if label == "LROH" else label for label in labels]
+    ax.legend(handles[::-1], labels[::-1], loc=[1.05,0], ncol=1, facecolor="white", fontsize='x-small')
+    ax.set_title("")
+    ax.spines['left'].set_color('black')
+    ax.spines['bottom'].set_color('black')
+
+    ylabel = ax.set_ylabel("operating expenses [bn EUR]")
+    ax.set_ylim([0, 550])
+    plt.savefig(snakemake.output.figure_opex, bbox_inches='tight', dpi=600)
+
+def electricity_prices(network):
     n = network
 
     rh_techs_elec = ['residential rural ground heat pump',
@@ -260,6 +286,62 @@ def electricity_prices(network, households):
 
     return energy_price_MWh
 
+def calc_opex(network, scenarioname):
+    mapping = network.loads.bus.map(network.buses.country)
+
+    query = "carrier == 'industry electricity'"
+    load = network.loads.query(query).groupby(mapping).sum().p_set
+    per_MWh = network.buses_t.marginal_price[network.buses.query("carrier == 'AC'").index].T.groupby(network.buses.country).mean().T
+    elec = (load * per_MWh).sum().sum()
+
+    query = "carrier == 'solid biomass for industry'"
+    load = network.loads.query(query).groupby(mapping).sum().p_set
+    per_MWh = network.buses_t.marginal_price[network.buses.query("carrier == 'solid biomass'").index].T.groupby(network.buses.country).mean().T
+    solid_biomass = (load * per_MWh).sum().sum()
+
+    query = "carrier == 'gas for industry'"
+    load = network.loads.query(query).groupby(mapping).sum().p_set
+    per_MWh = network.buses_t.marginal_price[network.buses.query("carrier == 'gas for industry'").index].T.groupby(network.buses.country).mean().T
+    gas = (load * per_MWh).sum().values[0]
+
+    query = "carrier == 'coal for industry'"
+    load = network.loads.query(query).groupby(mapping).sum().p_set.sum()
+    per_MWh = network.buses_t.marginal_price[network.buses.query("carrier == 'coal'").index].T.groupby(network.buses.country).mean().T
+    coal = (load * per_MWh).sum().values[0]
+
+    query = "carrier == 'H2 for industry'"
+    load = network.loads.query(query).groupby(mapping).sum().p_set
+    per_MWh = network.buses_t.marginal_price[network.buses.query("carrier == 'H2'").index].T.groupby(network.buses.country).mean().T
+    H2 = (load * per_MWh).sum().sum()
+
+    query = "carrier == 'H2 for industry'"
+    load = network.loads.query(query).groupby(mapping).sum().p_set
+    per_MWh = network.buses_t.marginal_price[network.buses.query("carrier == 'H2'").index].T.groupby(network.buses.country).mean().T
+    H2 = (load * per_MWh).sum().sum()
+
+    query = "carrier in ['industry methanol', 'shipping methanol']"
+    load = network.loads.query(query).p_set
+    per_MWh = network.buses_t.marginal_price[network.buses.query("carrier in ['industry methanol', 'shipping methanol']").index]
+    methanol = (load * per_MWh).sum().sum()
+
+    query = "carrier in ['shipping oil', 'naphtha for industry']"
+    load = network.loads.query(query).p_set
+    per_MWh = network.buses_t.marginal_price[network.buses.query("carrier in ['shipping oil', 'naphtha for industry']").index]
+    oil = (load * per_MWh).sum().sum()
+
+    query = "carrier == 'low-temperature heat for industry'"
+    load = network.loads.query(query).set_index("bus").p_set
+    per_MWh = network.buses_t.marginal_price[network.buses.query("carrier in ['urban central heat', 'services urban decentral heat']").index].clip(lower=0)
+    heat = (load * per_MWh).fillna(0).sum().sum()
+
+    opex = pd.DataFrame(
+        index=["Electricity", "Biomass", "Gas", "Coal", "Hydrogen", "Methanol", "Oil products", "Heat"],
+        columns=[scenarioname],
+        data=[elec, solid_biomass, gas, coal, H2, methanol, oil, heat]
+    )/1e9
+    return (
+        opex.loc[["Coal", "Methanol", "Hydrogen", "Oil products", "Biomass", "Heat", "Gas", "Electricity"]]
+    )
 
 
 if __name__ == "__main__":
@@ -287,11 +369,16 @@ def electricity_prices(network, households):
     # define scenario namings
     if planning_horizon == BAU_HORIZON:
         scenarios = {"BAU": "BAU"}
+        scenario_abbrev = scenarios
     else:
-        scenarios = {"flexible": "Optimal Renovation and Cost-Optimal Heating", 
-                     "retro_tes": "Optimal Renovation and Electric Heating", 
-                     "flexible-moderate": "Limited Renovation and Cost-Optimal Heating", 
+        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"}
+        scenario_abbrev = {"Optimal Renovation and Cost-Optimal Heating": "OROH",
+                     "Optimal Renovation and Electric Heating": "OREH",
+                     "Limited Renovation and Cost-Optimal Heating": "LROH",
+                     "No Renovation and Electric Heating": "NREH"}
 
     # load networks
     networks = {}
@@ -308,21 +395,25 @@ def electricity_prices(network, households):
     # calculate electricity bills per household for each network and electricity prices
     total_elec_bills = pd.DataFrame()
     total_elec_prices = pd.DataFrame()
+    industry_opex = pd.DataFrame()
     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
         # get electricity bills
-        elec_bills_household = electricity_bills(network, households)
+        elec_bills_household = electricity_bills(network, households).rename(name)
         # get electricity prices
-        elec_bills_MWh = electricity_prices(network, households)
+        elec_bills_MWh = electricity_prices(network).rename(name)
+        # OPEX for industry
+        opex = calc_opex(network, scenario_abbrev[name])
         # rename series name to scenario name
         elec_bills_household.name = name
         elec_bills_MWh.name = name
         # concatenate current results
         total_elec_bills = pd.concat([total_elec_bills, elec_bills_household.to_frame().T], axis=0)
         total_elec_prices = pd.concat([total_elec_prices, elec_bills_MWh.to_frame().T], axis=0)
+        industry_opex = pd.concat([industry_opex, opex], axis=1)
 
     # plot and store electricity bills
     if not total_elec_bills.empty:
@@ -332,3 +423,7 @@ def electricity_prices(network, households):
     if not total_elec_prices.empty:
         plot_electricity_cost(total_elec_prices, "prices")
 
+    # plot and store industry opex
+    if not industry_opex.empty:
+        plot_industry_opex(industry_opex)
+

plots/plot_heat_savings.py

@@ -74,6 +74,9 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
         ax = axes[i]
         ax.set_facecolor("whitesmoke")
 
+        # heat saved ratio
+        heat_saved_ratio = heat_demand.sum()["Heat savings by renovation"] / heat_demand.sum().sum()
+
         print("Hard coded coordinates on x-axis, selected for 3H granularity")
         if not full_year:
             where = [359, 695]
@@ -92,8 +95,46 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
         cumulative = (heat_demand / 1e3).cumsum(axis=1)
         ax.plot(cumulative.iloc[where[0]:where[1]], color='black', linewidth=0.5)
 
+        # x position for arrow
+        if not full_year:
+            x_loc = 436
+            y_shift_arrow = 0
+            x_shift_arrow = -24
+            x_shift_after = 33
+            x_shift_before = -65
+            x_shift_mid = -12
+            y_after = 1550
+            y_before = 2000
+            y_mid = 1800
+        else:
+            x_loc = 90
+            y_shift_arrow = -50
+            x_shift_arrow = 0
+            x_shift_after = 33
+            x_shift_before = -15
+            x_shift_mid = 20
+            y_after = 746
+            y_before = 1600
+            y_mid = 1150
+
+        if name == "Optimal Renovation and Cost-Optimal Heating":
+            # Heating demand before renovation
+            ax.annotate('heat demand before renovation', xy=(x_loc+x_shift_arrow, cumulative['Heat savings by renovation'][x_loc]+y_shift_arrow), xytext=(x_loc+x_shift_before, y_before),
+                arrowprops=dict(facecolor='black', edgecolor='black', arrowstyle='->', linewidth=0.75), fontsize=7, color="black")
+            # Heating demand after renovation
+            ax.annotate('heat demand after renovation', xy=(x_loc, cumulative['Net heat demand'][x_loc]+y_shift_arrow), xytext=(x_loc+x_shift_after, y_after),
+                arrowprops=dict(facecolor='black', edgecolor='black', arrowstyle='->', linewidth=0.75), fontsize=7, color="black")
+        if not name == "BAU":
+            # Heating savings
+            mid_point = cumulative.sum(axis=1) / 2
+            ax.annotate(f'saved heating demand ({100*heat_saved_ratio:.1f}%)', xy=(x_loc, mid_point[x_loc]+y_shift_arrow), xytext=(x_loc+x_shift_mid, y_mid),
+                arrowprops=dict(facecolor='black', edgecolor='black', arrowstyle='->', linewidth=0.75), fontsize=7, color="#a63f3f")
+
         ax.set_xlabel("", fontsize=12)
-        ax.set_ylim([1,1700])
+        if full_year:
+            ax.set_ylim([1,2000])
+        else:
+            ax.set_ylim([1,2200])
 
         if i < 2:
             ax.set_xticks([])
@@ -126,15 +167,15 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
         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
-    )
+    # 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)
@@ -275,6 +316,6 @@ def plot_flexibility(dict_elec):
     # 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_elec_consumption_for_heat(total_heat_data_with_flex, full_year=True)
     # plot heat flexibility
     plot_flexibility(total_flexibility)