Improve state of charge plot

Author: staadecker

switch_model/generators/extensions/storage.py

@@ -7,6 +7,8 @@
 storage, when to charge, energy accounting, etc.
 """
 import math
+
+import pandas as pd
 from scipy import fft
 
 from pyomo.environ import *
@@ -444,68 +446,72 @@ def post_solve(instance, outdir):
 @graph(
     "state_of_charge",
     title="State of Charge Throughout the Year",
-    supports_multi_scenario=True
+    supports_multi_scenario=True,
+    note="The daily charge/discharge amount is calculated as"
+         " the difference between the maximum and minimum"
+         " state of charge in a 1-day rolling window.\n"
+         "The black line is the 14-day rolling mean of the state of charge."
 )
 def graph_state_of_charge(tools):
+    # Each panel is a period and scenario
+    panel_group = ["period", "scenario_name"]
+    rolling_mean_window_size = '14D'
+
     # Get the total state of charge per timepoint and scenario
-    df = tools.get_dataframe("storage_dispatch")
-    df = df.groupby(["timepoint", "scenario_name"], as_index=False)["StateOfCharge"].sum()
+    soc = tools.get_dataframe("storage_dispatch.csv").rename({"StateOfCharge": "value"}, axis=1)
+    soc = soc.groupby(["timepoint", "scenario_name"], as_index=False).value.sum()
+    # Convert values to TWh
+    soc.value /= 1e6
     # Add datetime information
-    df = tools.transform.timestamp(df, key_col="timepoint")
+    soc = tools.transform.timestamp(soc, key_col="timepoint")[panel_group + ["datetime", "value"]]
     # Count num rows
-    num_periods = len(df["period"].unique())
+    num_periods = len(soc["period"].unique())
 
-    # Get the total capacity per period and scenario
-    capacity = tools.get_dataframe("storage_capacity.csv")
-    capacity = capacity.groupby(["period", "scenario_name"], as_index=False)["OnlineEnergyCapacityMWh"].sum()
+    # Used later
+    grouped_soc = soc.set_index("datetime").groupby(panel_group, as_index=False)
 
-    # Add the capacity to our dataframe
-    df = df.merge(
-        capacity,
-        on=["period", "scenario_name"],
-        validate="many_to_one",
-        how="left"
-    )
+    # Calculate the weekly SOC
+    weekly_soc = grouped_soc.rolling(rolling_mean_window_size, center=True) \
+        .value \
+        .mean() \
+        .reset_index()
 
-    # Convert values to TWh
-    df["StateOfCharge"] /= 1e6
-    df["OnlineEnergyCapacityMWh"] /= 1e6
+    # Get the total capacity per period and scenario
+    capacity = tools.get_dataframe("storage_capacity.csv")
+    capacity = capacity.groupby(panel_group, as_index=False)["OnlineEnergyCapacityMWh"]\
+        .sum()\
+        .rename({"OnlineEnergyCapacityMWh": "value"}, axis=1)
+    capacity.value /= 1e6
+    capacity["type"] = "Total Energy Capacity"
 
     # Add information regarding the diurnal cycle to the dataframe
     # Find the difference between the min and max for every day of the year
-    group = df.groupby(["period", "scenario_name", tools.pd.Grouper(freq="D", key="datetime")])["StateOfCharge"]
-    daily_size = (group.max() - group.min()).reset_index()
+    group = grouped_soc.rolling('D', center=True).value
+    daily_size = (group.max() - group.min()).reset_index().groupby(panel_group, as_index=False)
     # Find the mean between the difference of the min and max
-    daily_size = daily_size.groupby(["period", "scenario_name"], as_index=False).mean().reset_index()
-    # Add the mean to the dataframe under the name DiurnalCycle
-    df = df.merge(
-        daily_size.rename({"StateOfCharge": "DiurnalCycle"}, axis=1),
-        on=["period", "scenario_name"],
-        how="left"
-    )
+    avg_daily_size = daily_size.mean()[panel_group + ["value"]]
+    avg_daily_size["type"] = "Mean Daily Charge/Discharge"
+    max_daily_size = daily_size.max()[panel_group + ["value"]]
+    max_daily_size["type"] = "Maximum Daily Charge/Discharge"
 
     # Determine information for the labels
-    y_axis_max = df["OnlineEnergyCapacityMWh"].max()
-    label_offset = y_axis_max * 0.05
-    label_x_pos = df["datetime"].median()
+    y_axis_max = capacity.value.max()
+    label_x_pos = soc["datetime"].median()
+
+    hlines = pd.concat([capacity, avg_daily_size, max_daily_size])
+
     # For the max label
-    df["label_position"] = df["OnlineEnergyCapacityMWh"] + label_offset
-    df["label"] = df["OnlineEnergyCapacityMWh"].round(decimals=2)
-    # For the diurnal cycle label
-    df["diurnal_label"] = df["DiurnalCycle"].round(decimals=2)
-    df["diurnal_label_pos"] = df["DiurnalCycle"] + label_offset
+    hlines["label_pos"] = hlines.value + y_axis_max * 0.05
+    hlines["label"] = hlines.value.round(decimals=2)
 
     # Plot with plotnine
     pn = tools.pn
-    plot = pn.ggplot(df, pn.aes(x="datetime", y="StateOfCharge")) \
-           + pn.geom_line() \
+    plot = pn.ggplot(soc, pn.aes(x="datetime", y="value")) \
+           + pn.geom_line(color='gray') \
+           + pn.geom_line(data=weekly_soc, color='black') \
            + pn.labs(y="State of Charge (TWh)", x="Time of Year") \
-           + pn.geom_hline(pn.aes(yintercept="OnlineEnergyCapacityMWh"), linetype="dashed", color='blue') \
-           + pn.geom_text(pn.aes(label="label", x=label_x_pos, y="label_position"), fontweight="light", size="10") \
-           + pn.geom_hline(pn.aes(yintercept="DiurnalCycle"), linetype="dashed", color='red') \
-           + pn.geom_text(pn.aes(label="diurnal_label", x=label_x_pos, y="diurnal_label_pos"), fontweight="light", size="10", color="red")
-
-
+           + pn.geom_hline(pn.aes(yintercept="value", label="label", color="type"), data=hlines, linetype="dashed") \
+           + pn.geom_text(pn.aes(label="label", x=label_x_pos, y="label_pos"), data=hlines, fontweight="light", size="10")
     tools.save_figure(by_scenario_and_period(tools, plot, num_periods).draw())