power gap year-to-year distributions

run.py

@@ -1,34 +1,47 @@
-import pandas as pd 
+import pandas as pd
 import numpy as np
 
-from matplotlib import pyplot as plt 
+from matplotlib import pyplot as plt
 import matplotlib.dates as mdates
 
 import argparse
 import yaml
 
-from mix_simul.scenarios import Scenario 
+from mix_simul.scenarios import Scenario
 
 parser = argparse.ArgumentParser()
-parser.add_argument("--begin", help="begin date (YYYY-MM-DD), between 1985-01-01 and 2015-01-01", required=True)
-parser.add_argument("--end", help="end date (YYYY-MM-DD), between 1985-01-01 and 2015-01-01", required=True)
-parser.add_argument("--flexibility", help="enable load flexibility modeling", action="store_true")
-parser.add_argument("--scenarios", help="path to scenarios parameters yml file", default="scenarios/rte_2050.yml")
+parser.add_argument(
+    "--begin",
+    help="begin date (YYYY-MM-DD), between 1985-01-01 and 2015-01-01",
+    required=True,
+)
+parser.add_argument(
+    "--end",
+    help="end date (YYYY-MM-DD), between 1985-01-01 and 2015-01-01",
+    required=True,
+)
+parser.add_argument(
+    "--flexibility", help="enable load flexibility modeling", action="store_true"
+)
+parser.add_argument(
+    "--scenarios",
+    help="path to scenarios parameters yml file",
+    default="scenarios/rte_2050.yml",
+)
 args = parser.parse_args()
 
-with open(args.scenarios, 'r') as f:
+with open(args.scenarios, "r") as f:
     scenarios = yaml.load(f, Loader=yaml.FullLoader)
 
 potential = pd.read_parquet("data/potential.parquet")
-potential = potential.loc['1985-01-01 00:00:00':'2015-01-01 00:00:00', :]
+potential = potential.loc["1985-01-01 00:00:00":"2015-01-01 00:00:00", :]
 potential.fillna(0, inplace=True)
 
 begin = args.begin
 end = args.end
 flexibility = args.flexibility
 
-potential = potential.loc[(
-    slice(f'{begin} 00:00:00', f'{end} 00:00:00'), 'FR'), :]
+potential = potential.loc[(slice(f"{begin} 00:00:00", f"{end} 00:00:00"), "FR"), :]
 
 # intermittent sources potential
 p = potential[["onshore", "offshore", "solar"]].to_xarray().to_array()
@@ -40,46 +53,41 @@ n_scenarios = len(scenarios)
 
 fig, axes = plt.subplots(nrows=n_scenarios, ncols=2, sharex="col", sharey=True)
 w, h = fig.get_size_inches()
-fig.set_figwidth(w*1.5)
-fig.set_figheight(h*1.5)
+fig.set_figwidth(w * 1.5)
+fig.set_figheight(h * 1.5)
 
-fig_storage, axes_storage = plt.subplots(nrows=n_scenarios, ncols=2, sharex="col", sharey=True)
-fig_storage.set_figwidth(w*1.5)
-fig_storage.set_figheight(h*1.5)
+fig_storage, axes_storage = plt.subplots(
+    nrows=n_scenarios, ncols=2, sharex="col", sharey=True
+)
+fig_storage.set_figwidth(w * 1.5)
+fig_storage.set_figheight(h * 1.5)
 
-fig_dispatch, axes_dispatch = plt.subplots(nrows=n_scenarios, ncols=2, sharex="col", sharey=True)
-fig_dispatch.set_figwidth(w*1.5)
-fig_dispatch.set_figheight(h*1.5)
+fig_dispatch, axes_dispatch = plt.subplots(
+    nrows=n_scenarios, ncols=2, sharex="col", sharey=True
+)
+fig_dispatch.set_figwidth(w * 1.5)
+fig_dispatch.set_figheight(h * 1.5)
 
-fig_gap_distribution, axes_gap_distribution = plt.subplots(nrows=int(np.ceil(n_scenarios/2)), ncols=2, sharex="col", sharey=True)
-fig_gap_distribution.set_figwidth(w*1.5)
-fig_gap_distribution.set_figheight(h*1.5)
+fig_gap_distribution, axes_gap_distribution = plt.subplots(
+    nrows=int(np.ceil(n_scenarios / 2)), ncols=2, sharex="col", sharey=True
+)
+fig_gap_distribution.set_figwidth(w * 1.5)
+fig_gap_distribution.set_figheight(h * 1.5)
 
-months = [
-    "Février",
-    "Juin"
-]
+months = ["Février", "Juin"]
 
 labels = [
     "load (GW)",
     "production (GW)",
     "available power (production-storage) (GW)",
-    "power deficit"
+    "power deficit",
 ]
 
-labels_storage = [
-    "Batteries (TWh)",
-    "STEP (TWh)",
-    "P2G (TWh)"
-]
+labels_storage = ["Batteries (TWh)", "STEP (TWh)", "P2G (TWh)"]
 
-labels_dispatch = [
-    "Hydro (GW)",
-    "Biomass (GW)",
-    "Thermal (GW)"
-]
+labels_dispatch = ["Hydro (GW)", "Biomass (GW)", "Thermal (GW)"]
 
-date_fmt = mdates.DateFormatter('%d/%m')
+date_fmt = mdates.DateFormatter("%d/%m")
 
 # for step in np.linspace(start, stop, 2050-2022, True)[::-1]:
 row = 0
@@ -90,50 +98,74 @@ for scenario in scenarios:
     S, load, production, gap, storage, dp = scenario_model.run(times, p)
 
     print(f"{scenario}:", S, gap.max(), np.quantile(gap, 0.95))
-    print(f"exports: {np.minimum(np.maximum(-gap, 0), 39).sum()/1000} TWh; imports: {np.minimum(np.maximum(gap, 0), 39).sum()/1000} TWh")
-    print(f"dispatchable: " + ", ".join([f"{dp[i].sum()/1000:.2f} TWh" for i in range(dp.shape[0])]))
-
-    potential['load'] = load
-    potential['production'] = production
-    potential['available'] = production-np.diff(storage.sum(axis=0), append=0)
+    print(
+        f"exports: {np.minimum(np.maximum(-gap, 0), 39).sum()/1000} TWh; imports: {np.minimum(np.maximum(gap, 0), 39).sum()/1000} TWh"
+    )
+    print(
+        f"dispatchable: "
+        + ", ".join([f"{dp[i].sum()/1000:.2f} TWh" for i in range(dp.shape[0])])
+    )
+
+    potential["load"] = load
+    potential["production"] = production
+    potential["available"] = production - np.diff(storage.sum(axis=0), append=0)
+    potential["gap"] = gap
 
     for i in range(3):
-        potential[f"storage_{i}"] = np.diff(storage[i,:], append=0)#storage[i,:]/1000
-        potential[f"storage_{i}"] = storage[i,:]/1000
+        potential[f"storage_{i}"] = np.diff(
+            storage[i, :], append=0
+        )  # storage[i,:]/1000
+        potential[f"storage_{i}"] = storage[i, :] / 1000
 
     for i in range(dp.shape[0]):
-        potential[f"dispatch_{i}"] = dp[i,:]
+        potential[f"dispatch_{i}"] = dp[i, :]
 
     potential["dispatch"] = dp.sum(axis=0)
 
     data = [
-        potential.loc[(slice('2013-02-01 00:00:00',
-                       '2013-03-01 00:00:00'), 'FR'), :],
-        potential.loc[(slice('2013-06-01 00:00:00',
-                       '2013-07-01 00:00:00'), 'FR'), :]
+        potential.loc[(slice("2013-02-01 00:00:00", "2013-03-01 00:00:00"), "FR"), :],
+        potential.loc[(slice("2013-06-01 00:00:00", "2013-07-01 00:00:00"), "FR"), :],
     ]
 
     for col in range(2):
         ax = axes[row, col] if axes.ndim > 1 else axes[col]
-        ax.plot(data[col].index.get_level_values(0), data[col]
-                ["load"], label="adjusted load (GW)", lw=1)
-        ax.plot(data[col].index.get_level_values(0), data[col]
-                ["production"], label="production (GW)", ls="dotted", lw=1)
-        ax.plot(data[col].index.get_level_values(0), data[col]["available"],
-                label="available power (production-d(storage)/dt) (GW)", lw=1)
+        ax.plot(
+            data[col].index.get_level_values(0),
+            data[col]["load"],
+            label="adjusted load (GW)",
+            lw=1,
+        )
+        ax.plot(
+            data[col].index.get_level_values(0),
+            data[col]["production"],
+            label="production (GW)",
+            ls="dotted",
+            lw=1,
+        )
+        ax.plot(
+            data[col].index.get_level_values(0),
+            data[col]["available"],
+            label="available power (production-d(storage)/dt) (GW)",
+            lw=1,
+        )
 
         ax.fill_between(
             data[col].index.get_level_values(0),
             data[col]["available"],
             data[col]["load"],
             where=data[col]["load"] > data[col]["available"],
-            color='red',
-            alpha=0.15
+            color="red",
+            alpha=0.15,
         )
 
         ax.xaxis.set_major_formatter(date_fmt)
         ax.text(
-            0.5, 0.87, f"Scénario {scenario} ({months[col]})", ha='center', transform=ax.transAxes)
+            0.5,
+            0.87,
+            f"Scénario {scenario} ({months[col]})",
+            ha="center",
+            transform=ax.transAxes,
+        )
 
         ax.set_ylim(10, 210)
 
@@ -142,52 +174,103 @@ for scenario in scenarios:
             if i == 2:
                 base = 0
             else:
-                base = np.sum([data[col][f"storage_{j}"] for j in np.arange(i+1,3)], axis=0)
-            
-            ax.fill_between(data[col].index.get_level_values(0), base, base+data[col]
-                [f"storage_{i}"], label=f"storage {i}", alpha=0.5)
-
-            ax.plot(data[col].index.get_level_values(0), base+data[col]
-                [f"storage_{i}"], label=f"storage {i}", lw=0.25)
+                base = np.sum(
+                    [data[col][f"storage_{j}"] for j in np.arange(i + 1, 3)], axis=0
+                )
+
+            ax.fill_between(
+                data[col].index.get_level_values(0),
+                base,
+                base + data[col][f"storage_{i}"],
+                label=f"storage {i}",
+                alpha=0.5,
+            )
+
+            ax.plot(
+                data[col].index.get_level_values(0),
+                base + data[col][f"storage_{i}"],
+                label=f"storage {i}",
+                lw=0.25,
+            )
 
         ax.xaxis.set_major_formatter(date_fmt)
         ax.text(
-            0.5, 0.87, f"Scénario {scenario} ({months[col]})", ha='center', transform=ax.transAxes)
+            0.5,
+            0.87,
+            f"Scénario {scenario} ({months[col]})",
+            ha="center",
+            transform=ax.transAxes,
+        )
 
         ax = axes_dispatch[row, col] if axes.ndim > 1 else axes_dispatch[col]
         for i in range(dp.shape[0]):
             if i == 0:
                 base = 0
             else:
-                base = np.sum([data[col][f"dispatch_{j}"] for j in np.arange(i)], axis=0)
-            
-            ax.fill_between(data[col].index.get_level_values(0), base, base+data[col]
-                [f"dispatch_{i}"], label=f"dispatch {i}", alpha=0.5)
-
-            ax.plot(data[col].index.get_level_values(0), base+data[col]
-                [f"dispatch_{i}"], label=f"dispatch {i}", lw=0.25)
+                base = np.sum(
+                    [data[col][f"dispatch_{j}"] for j in np.arange(i)], axis=0
+                )
+
+            ax.fill_between(
+                data[col].index.get_level_values(0),
+                base,
+                base + data[col][f"dispatch_{i}"],
+                label=f"dispatch {i}",
+                alpha=0.5,
+            )
+
+            ax.plot(
+                data[col].index.get_level_values(0),
+                base + data[col][f"dispatch_{i}"],
+                label=f"dispatch {i}",
+                lw=0.25,
+            )
 
         ax.xaxis.set_major_formatter(date_fmt)
         ax.text(
-            0.5, 0.87, f"Scénario {scenario} ({months[col]})", ha='center', transform=ax.transAxes)
+            0.5,
+            0.87,
+            f"Scénario {scenario} ({months[col]})",
+            ha="center",
+            transform=ax.transAxes,
+        )
+
+        div = int(np.ceil(n_scenarios / 2))
+        ax = (
+            axes_gap_distribution[row % div, row // div]
+            if axes_gap_distribution.ndim > 1
+            else axes_gap_distribution[row % div]
+        )
 
-        div = int(np.ceil(n_scenarios/2))
-        ax = axes_gap_distribution[row%div, row//div] if axes_gap_distribution.ndim > 1 else axes_gap_distribution[row%div]
         hist, bin_edges = np.histogram(-gap, bins=1000)
         hist = np.cumsum(hist)
-        hist = 100*(hist - hist.min()) / hist.ptp()
-
+        hist = 100 * (hist - hist.min()) / hist.ptp()
         keep = np.abs(bin_edges[:-1]) < 50
-        ax.plot(bin_edges[:-1][keep], hist[keep], label="power gap")
-        ax.text(
-            0.5, 0.87, f"Scénario {scenario}", ha='center', transform=ax.transAxes)
-        
+        ax.plot(bin_edges[:-1][keep], hist[keep], lw=1, label="power gap", color="#ff7f00")
+
+        years = pd.date_range(start=begin, end=end, freq="Y")
+        for i in range(len(years) - 1):
+            year_data = potential.loc[
+                (slice(f"{years[i]} 00:00:00", f"{years[i+1]} 00:00:00"), "FR"), "gap"
+            ]
+
+            hist, bin_edges = np.histogram(-year_data, bins=1000)
+            hist = np.cumsum(hist)
+            hist = 100 * (hist - hist.min()) / hist.ptp()
+            keep = np.abs(bin_edges[:-1]) < 50
+            ax.plot(
+                bin_edges[:-1][keep], hist[keep], lw=0.5, alpha=0.2, color="#ff7f00"
+            )
+
+        ax.text(0.5, 0.87, f"Scénario {scenario}", ha="center", transform=ax.transAxes)
 
     row += 1
 
 for axs in [axes, axes_dispatch, axes_storage]:
     if axes.ndim > 1:
-        for label in axs[-1, 0].get_xmajorticklabels() + axs[-1, 1].get_xmajorticklabels():
+        for label in (
+            axs[-1, 0].get_xmajorticklabels() + axs[-1, 1].get_xmajorticklabels()
+        ):
             label.set_rotation(30)
             label.set_horizontalalignment("right")
     else:
@@ -198,28 +281,58 @@ for axs in [axes, axes_dispatch, axes_storage]:
 flex = "With" if flexibility else "Without"
 
 plt.subplots_adjust(wspace=0, hspace=0)
-fig.suptitle(f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
-fig.text(1, 0, 'Lucas Gautheron', ha="right")
-fig.legend(labels, loc='lower right', bbox_to_anchor=(1, -0.1),
-           ncol=len(labels), bbox_transform=fig.transFigure)
+fig.suptitle(
+    f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)"
+)
+fig.text(1, 0, "Lucas Gautheron", ha="right")
+fig.legend(
+    labels,
+    loc="lower right",
+    bbox_to_anchor=(1, -0.1),
+    ncol=len(labels),
+    bbox_transform=fig.transFigure,
+)
 fig.savefig("output/load_supply.png", bbox_inches="tight", dpi=200)
 
-fig_storage.suptitle(f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
-fig_storage.text(1, 0, 'Lucas Gautheron', ha="right")
-fig_storage.legend(labels_storage, loc='lower right', bbox_to_anchor=(1, -0.1),
-           ncol=len(labels_storage), bbox_transform=fig_storage.transFigure)
+fig_storage.suptitle(
+    f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)"
+)
+fig_storage.text(1, 0, "Lucas Gautheron", ha="right")
+fig_storage.legend(
+    labels_storage,
+    loc="lower right",
+    bbox_to_anchor=(1, -0.1),
+    ncol=len(labels_storage),
+    bbox_transform=fig_storage.transFigure,
+)
 fig_storage.savefig("output/storage.png", bbox_inches="tight", dpi=200)
 
-fig_dispatch.suptitle(f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
-fig_dispatch.text(1, 0, 'Lucas Gautheron', ha="right")
-fig_dispatch.legend(labels_dispatch, loc='lower right', bbox_to_anchor=(1, -0.1),
-           ncol=len(labels_dispatch), bbox_transform=fig_dispatch.transFigure)
+fig_dispatch.suptitle(
+    f"Simulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)"
+)
+fig_dispatch.text(1, 0, "Lucas Gautheron", ha="right")
+fig_dispatch.legend(
+    labels_dispatch,
+    loc="lower right",
+    bbox_to_anchor=(1, -0.1),
+    ncol=len(labels_dispatch),
+    bbox_transform=fig_dispatch.transFigure,
+)
 fig_dispatch.savefig("output/dispatch.png", bbox_inches="tight", dpi=200)
 
-fig_gap_distribution.suptitle(f"Power gap cumulative distribution (%)\nSimulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)")
-fig_gap_distribution.legend(["Power gap (available-load) (GW)"], loc='lower right', bbox_to_anchor=(1, -0.1),
-           ncol=1, bbox_transform=fig_dispatch.transFigure)
-fig_gap_distribution.text(1, 0, 'Lucas Gautheron', ha="right")
-fig_gap_distribution.savefig("output/gap_distribution.png", bbox_inches="tight", dpi=200)
+fig_gap_distribution.suptitle(
+    f"Power gap cumulative distribution (%)\nSimulations based on {begin}--{end} weather data.\n{flex} consumption flexibility; no nuclear seasonality (unrealistic)"
+)
+fig_gap_distribution.legend(
+    ["Power gap (available-load) (GW)"],
+    loc="lower right",
+    bbox_to_anchor=(1, -0.1),
+    ncol=1,
+    bbox_transform=fig_dispatch.transFigure,
+)
+fig_gap_distribution.text(1, 0, "Lucas Gautheron", ha="right")
+fig_gap_distribution.savefig(
+    "output/gap_distribution.png", bbox_inches="tight", dpi=200
+)
 
 plt.show()