From 295ef8ea24d09e2d646c33437db8d0f699870b06 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Tue, 16 Jul 2024 19:55:47 +0500
Subject: [PATCH 01/21] reorder co2 emitting technologies

---
 plots/plot_co2_level.py | 23 +++++++----------------
 1 file changed, 7 insertions(+), 16 deletions(-)

diff --git a/plots/plot_co2_level.py b/plots/plot_co2_level.py
index b4f6495..f0a00e3 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",

From 4b015ab6096d2e104915da8a37ff514d540c30f4 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Wed, 17 Jul 2024 17:02:49 +0500
Subject: [PATCH 02/21] update capacity expansion plot

---
 plots/plot_capacity_expansion.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/plots/plot_capacity_expansion.py b/plots/plot_capacity_expansion.py
index 03002db..4c0d4d4 100644
--- a/plots/plot_capacity_expansion.py
+++ b/plots/plot_capacity_expansion.py
@@ -52,6 +52,16 @@ def plot_capacities(caps_df, clusters, planning_horizon, plot_width=7):
 
     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(

From 6fc3df50d01f4fd06a0420a7849edd68eac54a04 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Fri, 19 Jul 2024 18:57:00 +0500
Subject: [PATCH 03/21] plot for 2 weeks

---
 plots/plot_electricity_for_heat.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/plots/plot_electricity_for_heat.py b/plots/plot_electricity_for_heat.py
index 7839850..28f6a47 100644
--- a/plots/plot_electricity_for_heat.py
+++ b/plots/plot_electricity_for_heat.py
@@ -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"]
@@ -96,19 +96,19 @@ def plot_elec_consumption_for_heat(dict_elec):
         )
 
         ax.set_xlabel("", fontsize=12)
-        ax.set_ylim([1,2000])
-        ax.set_yscale("log")  # Set the y-axis to logarithmic scale
+        ax.set_ylim([1,1600])
+        # ax.set_yscale("log")  # Set the y-axis to logarithmic scale
 
-        # 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())
+        # # 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())
 
         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]

From 46715ed3033b382a973e8c6fd9b79cbb6f57744f Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Fri, 19 Jul 2024 22:21:11 +0500
Subject: [PATCH 04/21] add black border lines for separate technologies

---
 plots/plot_electricity_for_heat.py | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/plots/plot_electricity_for_heat.py b/plots/plot_electricity_for_heat.py
index 28f6a47..fc40030 100644
--- a/plots/plot_electricity_for_heat.py
+++ b/plots/plot_electricity_for_heat.py
@@ -95,14 +95,11 @@ 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,1600])
-        # ax.set_yscale("log")  # Set the y-axis to logarithmic scale
-
-        # # 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())
 
         if i < 3:
             ax.set_xticks([])

From 74610a26073a8f21aacd0c0bb01c8f4d887f3c27 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Sun, 21 Jul 2024 14:23:07 +0500
Subject: [PATCH 05/21] add heat saving plot

---
 Snakefile                  |  37 ++++++++
 plots/plot_heat_savings.py | 184 +++++++++++++++++++++++++++++++++++++
 2 files changed, 221 insertions(+)
 create mode 100644 plots/plot_heat_savings.py

diff --git a/Snakefile b/Snakefile
index 9722223..7fbc3c6 100644
--- a/Snakefile
+++ b/Snakefile
@@ -109,6 +109,33 @@ 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",
+    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"],
+        ),
+
+
 rule plot_electricity_generation:
     params:
         clusters=config["plotting"]["clusters"],
@@ -366,6 +393,16 @@ 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
             + "table_heat_pumps_{clusters}.csv",
diff --git a/plots/plot_heat_savings.py b/plots/plot_heat_savings.py
new file mode 100644
index 0000000..fc8bef4
--- /dev/null
+++ b/plots/plot_heat_savings.py
@@ -0,0 +1,184 @@
+# -*- 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):
+    # 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)
+    # 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 plot_elec_consumption_for_heat(dict_elec, full_year=False):
+    # set heights for each subplots
+    if "BAU" in dict_elec.keys():
+        heights = [1.4]
+    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 = ["#2fb537", "#48f74f", "#d8f9b8", "#db6a25"]
+
+        (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 = "Heating Demands [GW]"
+        axes[0].set_ylabel(ylabel, fontsize=10)
+    else:
+        ylabel = "Heating Demands [GW]"
+        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)
+
+
+
+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 = {}
+    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
+    
+    # plot heat demand data for short period
+    plot_elec_consumption_for_heat(total_heat_data)
+    # plot heat demand data for full year
+    plot_elec_consumption_for_heat(total_heat_data, full_year=True)

From fd08e18899c8af75b2c3cf68bc69872ea7eff130 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Mon, 22 Jul 2024 19:14:47 +0500
Subject: [PATCH 06/21] change color schemes

---
 plots/plot_heat_savings.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/plots/plot_heat_savings.py b/plots/plot_heat_savings.py
index fc8bef4..4caa2b5 100644
--- a/plots/plot_heat_savings.py
+++ b/plots/plot_heat_savings.py
@@ -56,7 +56,7 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
             heat_demand = heat_demand.resample('24H').max()
         
         heat_demand = heat_demand.reset_index(drop=True)
-        colors = ["#2fb537", "#48f74f", "#d8f9b8", "#db6a25"]
+        colors = ["#a63f3f", "#e39191"]
 
         (heat_demand/1e3).iloc[where[0] : where[1]].plot(
             kind='area', stacked=True, color=colors, legend=False,

From c1f335b8594434be45bd389f3fff9bb8eaeaebd4 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Tue, 23 Jul 2024 16:03:13 +0500
Subject: [PATCH 07/21] add flexibility into 2 week plot

---
 plots/plot_heat_savings.py | 33 ++++++++++++++++++++++++---------
 1 file changed, 24 insertions(+), 9 deletions(-)

diff --git a/plots/plot_heat_savings.py b/plots/plot_heat_savings.py
index 4caa2b5..35725b3 100644
--- a/plots/plot_heat_savings.py
+++ b/plots/plot_heat_savings.py
@@ -19,23 +19,35 @@
                      LINE_LIMITS, CO2L_LIMITS, BAU_HORIZON
 
 
-def get_heat_data(n):
+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)
-    # 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"]
+    # 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 plot_elec_consumption_for_heat(dict_elec, full_year=False):
     # set heights for each subplots
     if "BAU" in dict_elec.keys():
-        heights = [1.4]
+        heights = [1.2]
     else:
         heights = [1.4] * 3
     fig = plt.figure(figsize=(6.4, sum(heights)))
@@ -105,7 +117,7 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
         reversed(handles1[0:7]),
         [
             "Heat savings by renovation",
-            "Net heat demand"
+            "Net heat demand",
         ],
         loc=[1.02, -.2], fontsize=10
     )
@@ -170,6 +182,7 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
     change_path_to_base()
    
     total_heat_data = {}
+    total_heat_data_with_flex = {}
     for name, network in networks.items():
         if network is None:
             # Skip further computation for this scenario if network is not loaded
@@ -177,8 +190,10 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
             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
     
     # plot heat demand data for short period
-    plot_elec_consumption_for_heat(total_heat_data)
+    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)

From 5a2c5c5f1f67ecd622867bf9b978512d5dfbfa85 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Thu, 25 Jul 2024 14:46:14 +0500
Subject: [PATCH 08/21] TW for installed capacity

---
 plots/plot_total_costs.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

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:

From e3e318b92a2ea3272d49e42668189b8e1c91e0e5 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Thu, 25 Jul 2024 15:05:58 +0500
Subject: [PATCH 09/21] TW for capacity expansion

---
 plots/plot_capacity_expansion.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/plots/plot_capacity_expansion.py b/plots/plot_capacity_expansion.py
index 4c0d4d4..6975513 100644
--- a/plots/plot_capacity_expansion.py
+++ b/plots/plot_capacity_expansion.py
@@ -39,14 +39,14 @@ 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])
 
@@ -85,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:

From d4b923f32f1a543f9d1181ee9b57576b1c0ab9cd Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Thu, 25 Jul 2024 16:15:07 +0500
Subject: [PATCH 10/21] plot in Billion tonn_CO2

---
 plots/plot_co2_level.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/plots/plot_co2_level.py b/plots/plot_co2_level.py
index f0a00e3..ebbb33c 100644
--- a/plots/plot_co2_level.py
+++ b/plots/plot_co2_level.py
@@ -186,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)
@@ -220,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

From 578d12f8a7e25f20639664b80ddac013f0a760dc Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Thu, 25 Jul 2024 18:08:19 +0500
Subject: [PATCH 11/21] total heat pump capacity in GWel

---
 plots/table_heat_pumps.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

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

From 5c2ca7de177eedb4d6b2d6fd46311a533cb9c018 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Sat, 27 Jul 2024 23:42:47 +0500
Subject: [PATCH 12/21] convert saving into trillion m3

---
 plots/plot_co2_level.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/plots/plot_co2_level.py b/plots/plot_co2_level.py
index ebbb33c..d1155e2 100644
--- a/plots/plot_co2_level.py
+++ b/plots/plot_co2_level.py
@@ -366,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

From c65d757db6fed41e7ce75e2c68e58a063ea383a2 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Mon, 29 Jul 2024 13:41:43 +0500
Subject: [PATCH 13/21] variable ylim for horizons

---
 plots/plot_electricity_for_heat.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/plots/plot_electricity_for_heat.py b/plots/plot_electricity_for_heat.py
index fc40030..230dd1a 100644
--- a/plots/plot_electricity_for_heat.py
+++ b/plots/plot_electricity_for_heat.py
@@ -66,7 +66,7 @@ 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]
@@ -99,7 +99,10 @@ def plot_elec_consumption_for_heat(dict_elec):
         ax.plot(cumulative.iloc[where[0]:where[1]], color='black', linewidth=0.5)
 
         ax.set_xlabel("", fontsize=12)
-        ax.set_ylim([1,1600])
+        # get in y limits for all horizons
+        y_limits = {'2020':1600, '2030':1100, '2040':900, '2050': 900}
+
+        ax.set_ylim([1, y_limits[horizon]])
 
         if i < 3:
             ax.set_xticks([])
@@ -200,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)
 

From 167377bb7e7f3843df575379416d82fbdb3a6beb Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Mon, 29 Jul 2024 14:03:40 +0500
Subject: [PATCH 14/21] add flexibility plot

---
 Snakefile                  | 11 ++++++
 plots/plot_heat_savings.py | 80 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 91 insertions(+)

diff --git a/Snakefile b/Snakefile
index 7fbc3c6..13b471a 100644
--- a/Snakefile
+++ b/Snakefile
@@ -116,6 +116,7 @@ rule plot_heat_saving:
     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:
@@ -134,6 +135,11 @@ rule plot_heat_savings:
             + "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:
@@ -403,6 +409,11 @@ rule plot_all:
             + "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/plots/plot_heat_savings.py b/plots/plot_heat_savings.py
index 35725b3..25e5158 100644
--- a/plots/plot_heat_savings.py
+++ b/plots/plot_heat_savings.py
@@ -44,6 +44,20 @@ def get_heat_data(n, with_flexibility=False):
     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():
@@ -135,6 +149,66 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
         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)
+        ax.plot((heat_demand["heat"]/1e3).iloc[where[0]:where[1]], color="black", linewidth=0.5)
+
+        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
+
+    if len(dict_elec.keys()) == 1:
+        ylabel = "Heating Demands [GW]"
+        axes[0].set_ylabel(ylabel, fontsize=10)
+    else:
+        ylabel = "Heating Demands [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():
@@ -183,6 +257,7 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
    
     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
@@ -192,8 +267,13 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
         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

From 622815897297b44e83fac9f7d476a2052073aa18 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Mon, 29 Jul 2024 14:48:09 +0500
Subject: [PATCH 15/21] add legend and fix figure size

---
 plots/plot_heat_savings.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/plots/plot_heat_savings.py b/plots/plot_heat_savings.py
index 25e5158..43d8a7d 100644
--- a/plots/plot_heat_savings.py
+++ b/plots/plot_heat_savings.py
@@ -61,7 +61,7 @@ def get_flexibility(network):
 def plot_elec_consumption_for_heat(dict_elec, full_year=False):
     # set heights for each subplots
     if "BAU" in dict_elec.keys():
-        heights = [1.2]
+        heights = [1.3]
     else:
         heights = [1.4] * 3
     fig = plt.figure(figsize=(6.4, sum(heights)))
@@ -170,8 +170,8 @@ def plot_flexibility(dict_elec):
         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)
-        ax.plot((heat_demand["heat"]/1e3).iloc[where[0]:where[1]], color="black", linewidth=0.5)
+        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])
@@ -199,6 +199,7 @@ def plot_flexibility(dict_elec):
         ax.set_title(name, fontsize=10)
         i+= 1
 
+    axes[0].legend(loc=[1.02, -.2], fontsize=10)
     if len(dict_elec.keys()) == 1:
         ylabel = "Heating Demands [GW]"
         axes[0].set_ylabel(ylabel, fontsize=10)

From c5ac1145c3243547d705785c2cdc4c8343778fbb Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Mon, 29 Jul 2024 20:35:32 +0500
Subject: [PATCH 16/21] add historic prices

---
 plots/_helpers.py               |  7 +++++++
 plots/plot_electricity_bills.py | 23 +++++++++++++++++------
 2 files changed, 24 insertions(+), 6 deletions(-)

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_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)

From 4e01e42f16d2f565af4d20acab62ea59e7acc66c Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Mon, 29 Jul 2024 20:41:52 +0500
Subject: [PATCH 17/21] set numeric focus to retro_tes 2050 to prevent numeric
 issue

---
 configs/EEE_study/config.retro_tes-industry_2050.yaml | 1 +
 1 file changed, 1 insertion(+)

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

From 5e49f7f443116f71cc1d72c150959e583a8a37fd Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Mon, 29 Jul 2024 20:52:35 +0500
Subject: [PATCH 18/21] adjust ylabel

---
 plots/plot_heat_savings.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/plots/plot_heat_savings.py b/plots/plot_heat_savings.py
index 43d8a7d..4c40a3f 100644
--- a/plots/plot_heat_savings.py
+++ b/plots/plot_heat_savings.py
@@ -136,10 +136,10 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
         loc=[1.02, -.2], fontsize=10
     )
     if len(dict_elec.keys()) == 1:
-        ylabel = "Heating Demands [GW]"
+        ylabel = "Flexibility [GW]"
         axes[0].set_ylabel(ylabel, fontsize=10)
     else:
-        ylabel = "Heating Demands [GW]"
+        ylabel = "Flexibility [GW]"
         axes[1].set_ylabel(ylabel, fontsize=10)
         # axes[1].yaxis.set_label_coords(-0.1, 1)
     plt.subplots_adjust(hspace=0.3)

From 3b6a2e7a792e47b15e75f63deec5cbfbdd092689 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Mon, 29 Jul 2024 20:53:25 +0500
Subject: [PATCH 19/21] fix ylabel

---
 plots/plot_heat_savings.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/plots/plot_heat_savings.py b/plots/plot_heat_savings.py
index 4c40a3f..abaf602 100644
--- a/plots/plot_heat_savings.py
+++ b/plots/plot_heat_savings.py
@@ -136,10 +136,10 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
         loc=[1.02, -.2], fontsize=10
     )
     if len(dict_elec.keys()) == 1:
-        ylabel = "Flexibility [GW]"
+        ylabel = "Heat demand [GW]"
         axes[0].set_ylabel(ylabel, fontsize=10)
     else:
-        ylabel = "Flexibility [GW]"
+        ylabel = "Heat demand [GW]"
         axes[1].set_ylabel(ylabel, fontsize=10)
         # axes[1].yaxis.set_label_coords(-0.1, 1)
     plt.subplots_adjust(hspace=0.3)
@@ -201,10 +201,10 @@ def plot_flexibility(dict_elec):
 
     axes[0].legend(loc=[1.02, -.2], fontsize=10)
     if len(dict_elec.keys()) == 1:
-        ylabel = "Heating Demands [GW]"
+        ylabel = "Flexibility [GW]"
         axes[0].set_ylabel(ylabel, fontsize=10)
     else:
-        ylabel = "Heating Demands [GW]"
+        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)

From dbabb2af900d4269986f4596b9ee22e75b1b9084 Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Mon, 29 Jul 2024 22:12:04 +0500
Subject: [PATCH 20/21] set GW_th

---
 plots/plot_heat_savings.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/plots/plot_heat_savings.py b/plots/plot_heat_savings.py
index abaf602..d416cd7 100644
--- a/plots/plot_heat_savings.py
+++ b/plots/plot_heat_savings.py
@@ -136,10 +136,10 @@ def plot_elec_consumption_for_heat(dict_elec, full_year=False):
         loc=[1.02, -.2], fontsize=10
     )
     if len(dict_elec.keys()) == 1:
-        ylabel = "Heat demand [GW]"
+        ylabel = "Heat demand [$GW_{th}$]"
         axes[0].set_ylabel(ylabel, fontsize=10)
     else:
-        ylabel = "Heat demand [GW]"
+        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)

From 135bd712d32655acb7217311b5a811d20036e87c Mon Sep 17 00:00:00 2001
From: yerbol-akhmetov <yerbol.akhmetov.93@gmail.com>
Date: Tue, 30 Jul 2024 13:42:51 +0500
Subject: [PATCH 21/21] update infra saving to have installed capacity and
 capital costs

---
 plots/plot_electricity_for_heat.py |  6 ++---
 plots/table_infra_savings.py       | 43 +++++++++---------------------
 2 files changed, 16 insertions(+), 33 deletions(-)

diff --git a/plots/plot_electricity_for_heat.py b/plots/plot_electricity_for_heat.py
index 230dd1a..f6a11bc 100644
--- a/plots/plot_electricity_for_heat.py
+++ b/plots/plot_electricity_for_heat.py
@@ -69,7 +69,7 @@ def get_elec_consumption_for_heat(n):
 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)))
@@ -136,10 +136,10 @@ def plot_elec_consumption_for_heat(dict_elec, horizon):
         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)
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]"