fixed bug causing coal not to phase out

config/config.yaml

@@ -22,15 +22,13 @@ scenario:
   simpl:
   - ''
   ll:  # allowed transmission line volume expansion, can be any float >= 1.0 with a prefix v|c (today) or "copt"
-  - v1.0
-  # clusters:  # number of nodes in Europe, any integer between 37 (1 node per country-zone) and several hundred
-  # - 37
+  - v1.5
   gb_regions:
   - 'eso' # regions can be 'eso' zones (emphasizing grid conjestion), or 'dno'
   fes: # Future energy scenario considered, can be 
   # - 'CT' # Consumer Transformation
-  - 'FS' # Falling Short
-  # - 'LW' # Leading the Way 
+  # - 'FS' # Falling Short
+  - 'LW' # Leading the Way 
   # - 'ST' # System Transformation
   year: # investment years for myopic and perfect; for overnight, year of cost assumptions can be different and is defined under 'costs'
   - 2023
@@ -41,8 +39,8 @@ scenario:
   - 2045
   - 2050
   flexopts:
-  # - reg-ss-bev-heat
-  - heat
+  - reg-ss-v2g-heat
+  # - heat
   - ''
   # reg  # regular weekly demand flexibility (2 days max), one hour
   # ss   # one hour events, any weekday in winter, Saving Sessions style
@@ -97,6 +95,8 @@ flexibility:
     - 'offwind-ac'
     - 'offwind-dc'
     - 'nuclear'
+
+  coal_phaseout_year_uk: 2026  # assuming no coal capacity in the system after this year
 
 
 # here is a list of all possible countries: 

scripts/build_fes_constraints.py

@@ -13,6 +13,7 @@
 """
 
 import pandas as pd
+import numpy as np
 
 from _fes_helpers import get_data_point
 from _helpers import configure_logging
@@ -76,20 +77,27 @@
         "value": val,
         "sense": "==",})
 
+    coal_phaseout = snakemake.config["flexibility"]["coal_phaseout_uk"]
+    caps.loc[len(caps)] = pd.Series({
+        "carrier": "coal",
+        "attr": "p_nom",
+        "value": np.interp(year, [2023, coal_phaseout], [2_520., 0]),
+        "sense": "==",})
+
     """
     val = get_data_point("bioenergy_capacity", fes, year)
     caps.loc[len(caps)] = pd.Series({
         "carrier": "biomass",
         "attr": "p_nom",
         "value": val,
         "sense": "==",})
-    
+
     val = get_data_point("bioenergy_ccs_capacity", fes, year)
     caps.loc[len(caps)] = pd.Series({
         "carrier": "biomass ccs",
         "attr": "p_nom",
         "value": val,
         "sense": "==",})
     """
-    
+
     caps.to_csv(snakemake.output["capacity_constraints"])

scripts/prepare_network.py

@@ -520,12 +520,6 @@ def add_heat_pump_load(
     )
 
     gb_nodes = heat_demand.columns
-    print("heat demand index")
-    print(gb_nodes)
-
-    #########################################################################################################
-    # add infrastructure for flexibility
-    # grid -> house -> link to thermal inertia -> thermal inertia store -> link to house     
 
     hp_heat_demand = heat_demand / heat_demand.sum().sum() * hp_load_future
 
@@ -559,19 +553,13 @@ def add_heat_pump_load(
         marginal_cost=0.,
     )
 
-    import matplotlib.pyplot as plt
-    fig, ax = plt.subplots(1, 1, figsize=(16, 4))
-
-    hp_heat_demand.iloc[:200].sum(axis=1).plot(ax=ax, label="hp_heat_demand")
-    plt.savefig("hp_heat_demand.png")
-
-    plt.show()
-
     n.loads_t.p_set[gb_regions] -= heat_demand / heat_demand.sum().sum() * hp_load_base
 
     share_smart_tariff = snakemake.config["flexibility"]["heat_share_smart_tariff"]
     heatflex = "heat" in snakemake.wildcards["flexopts"].split("-")
 
+    # add infrastructure for flexibility
+    # grid -> house -> link to thermal inertia -> thermal inertia store -> link to house     
     if share_smart_tariff > 0. and heatflex:
 
         logger.info("Adding heat flexibility")
@@ -621,17 +609,6 @@ def add_heat_pump_load(
         p_max_pu = (daily_p_max / p_nom).mul(charging_window, axis=0)
         p_min_pu = - (daily_p_max / p_nom).mul(discharging_window, axis=0)
 
-        print("p_nom")
-        print(p_nom)        
-        print("p_max_pu")
-        print(p_max_pu)
-        print("p_min_pu")
-        print(p_min_pu)
-
-        # p_nom.index += " thermal ine"
-        # p_max_pu.columns += " lhv heat"
-        # p_min_pu.columns += " lhv heat"
-
         daily_e_max = (
             hp_heat_demand
             .rolling(shift_size)
@@ -645,9 +622,6 @@ def add_heat_pump_load(
         e_nom = daily_e_max.max() * share_smart_tariff
         e_max_pu = (daily_e_max / e_nom).mul(store_use_window, axis=0)
 
-        # e_nom.index += " thermal inertia"
-        # e_max_pu.columns += " thermal inertia"
-
         n.madd(
             "Bus",
             gb_nodes,
@@ -680,41 +654,6 @@ def add_heat_pump_load(
 
 
 
-
-        print("daily max")
-        print(daily_p_max.head())
-
-        print("gb nodes")
-        print(gb_nodes)
-
-        import matplotlib.pyplot as plt
-        fig, axs = plt.subplots(3, 1, figsize=(16, 10))
-
-        steps = 200
-
-        e_max_pu.mean(axis=1).iloc[:steps].plot(ax=axs[0], label="e_nom_max")
-        p_max_pu.mean(axis=1).iloc[:steps].plot(ax=axs[0], label="p_nom_max")
-        p_min_pu.mean(axis=1).iloc[:steps].plot(ax=axs[0], label="p_nom_max")
-
-        p_nom.plot.bar(ax=axs[1])
-        e_nom.plot.bar(ax=axs[2])
-
-        # daily_e_max.mean(axis=1).iloc[:steps].mean(axis=1).plot(ax=axs[1], label="daily max")
-        # daily_p_max.mean(axis=1).iloc[:steps].mean(axis=1).plot(ax=axs[1], label="daily max")
-        # store_use_max.iloc[:steps].mean(axis=1).plot(ax=ax, label="store use max")
-        # discharging_max.iloc[:steps].mean(axis=1).plot(ax=ax, label="discharging max", linestyle=":")
-        # charging_max.iloc[:steps].mean(axis=1).plot(ax=ax, label="charging max", linestyle=":")
-
-        axs[0].legend()
-
-        plt.savefig("heat_flexibility.pdf")
-        plt.show()
-
-
-    #########################################################################################################
-    # n.loads_t.p_set[gb_regions] += heat_demand
-
-
 def add_bev(n, transport_config):
     """Adds BEV load and respective stores units;
     adapted from `add_land_transport` method in `scripts/prepare_sector_network.py`"""
@@ -1296,13 +1235,15 @@ def add_batteries(n):
 
     flexopts = snakemake.wildcards.flexopts.split("-")
 
-    print(flexopts)
-    if not (len(flexopts) == 1 and "" in flexopts):
-        assert sum(list(map(lambda x: x in {"ss", "reg", "bev", "heat"}, flexopts))) == len(flexopts), (
-            "flexopts must be a combination of 'ss', 'reg', 'bev', 'heat', currently is {}".format(
-                snakemake.wildcards.flexopts
-            )
-        )
+
+    logger.warning("Flexopts wildcard constraints commented out currently.")
+    # print(flexopts)
+    # if not (len(flexopts) == 1 and "" in flexopts):
+    #     assert sum(list(map(lambda x: x in {"ss", "reg", "bev", "heat"}, flexopts))) == len(flexopts), (
+    # "flexopts must be a combination of 'ss', 'reg', 'bev', 'heat', currently is {}".format(
+    #             snakemake.wildcards.flexopts
+    #        )
+    #)
 
     logger.info(f"Using Flexibility Options: {flexopts}")
 

scripts/summarize_gb.py

@@ -88,7 +88,7 @@ def to_csv(df):
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
 
-        snakemake = mock_snakemake("make_summary")
+        snakemake = mock_snakemake("summarize_gb")
 
     logging.basicConfig(level=snakemake.config["logging"]["level"])
 
@@ -188,10 +188,9 @@ def sort_mindex(index, networks_dict):
                 print("Setting snapshots")
                 print(snapshots[(fes, year)])
                 print("We took {} percent of snapshots".format(len(snapshots[(fes, year)])/8760))
-
-
-            # inflow = inflow.loc[snapshots[(fes, year)]]
-            # outflow = outflow.loc[snapshots[(fes, year)]]
+
+            inflow = inflow.loc[snapshots[(fes, year)]]
+            outflow = outflow.loc[snapshots[(fes, year)]]
 
             df[output] = df[output].reindex(
                 df[output].index.union(