Merge pull request #103 from staadecker/plots

Plot improvements

Author: pesap

papers/Martin_Staadecker_Value_of_LDES_and_Factors/LDES_paper_graphs/plots.py

@@ -251,11 +251,12 @@ def graph_energy_balance(tools):
         load_balance["normalized_energy_balance_duals_dollar_per_mwh"], errors="coerce") / 10
     load_balance = load_balance[["energy_balance_duals", "time_row"]]
     load_balance = load_balance.pivot(columns="time_row", values="energy_balance_duals")
+    percent_of_zeroes = sum(load_balance == 0) / len(load_balance) * 100
     # Don't include the zero-valued duals
     load_balance = load_balance.replace(0, tools.np.nan)
     if load_balance.count().sum() != 0:
         load_balance.plot.box(
-            ax=tools.get_axes(),
+            ax=tools.get_axes(note=f"{percent_of_zeroes:.1f}% of duals are zero"),
             xlabel='Period',
             ylabel='Energy balance duals (cents/kWh)',
             showfliers=False

switch_model/balancing/load_zones.py

@@ -776,6 +776,23 @@ def graph_buildout_per_tech(tools):
     ax = tools.get_axes()
     df.plot(ax=ax, kind='line', color=colors, xlabel='Period', marker="x")
     # Set the y-axis to use percent
-    ax.yaxis.set_major_formatter(tools.mplt.ticker.PercentFormatter(1.0))
+    ax.yaxis.set_major_formatter(tools.plt.ticker.PercentFormatter(1.0))
     # Horizontal line at 100%
     ax.axhline(y=1, linestyle="--", color='b')
+
+@graph(
+    "buildout_map",
+    title="Map of online capacity per load zone."
+)
+def buildout_map(tools):
+    buildout = tools.get_dataframe("gen_cap.csv").rename({"GenCapacity": "value"}, axis=1)
+    buildout = tools.transform.gen_type(buildout)
+    buildout = buildout.groupby(["gen_type", "gen_load_zone"], as_index=False)["value"].sum()
+    ax = tools.maps.graph_pie_chart(buildout)
+    transmission = tools.get_dataframe("transmission.csv", convert_dot_to_na=True).fillna(0)
+    transmission = transmission.rename({"trans_lz1": "from", "trans_lz2": "to", "BuildTx": "value"}, axis=1)
+    transmission = transmission[["from", "to", "value", "PERIOD"]]
+    transmission = transmission.groupby(["from", "to", "PERIOD"], as_index=False).sum().drop("PERIOD", axis=1)
+    # Rename the columns appropriately
+    transmission.value *= 1e-3
+    tools.maps.graph_transmission(transmission, cutoff=0.1, ax=ax, legend=True)

switch_model/generators/core/build.py

@@ -806,8 +806,117 @@ def graph_curtailment_per_tech(tools):
     df.plot(ax=ax, kind='line',  xlabel='Period', marker="x", **kwargs)
 
     # Set the y-axis to use percent
-    ax.yaxis.set_major_formatter(tools.mplt.ticker.PercentFormatter(1.0))
+    ax.yaxis.set_major_formatter(tools.plt.ticker.PercentFormatter(1.0))
     # Horizontal line at 100%
     # ax.axhline(y=1, linestyle="--", color='b')
 
 
+@graph(
+    "energy_balance_2",
+    title="Balance between demand, generation and storage for last period",
+    note="Dashed green and red lines are total generation and total demand (incl. transmission losses),"
+         " respectively.\nDotted line is the total state of charge (scaled for readability)."
+         "\nWe used a 14-day rolling mean to smoothen out values."
+)
+def graph_energy_balance_2(tools):
+    # Get dispatch dataframe
+    dispatch = tools.get_dataframe("dispatch.csv", usecols=[
+        "timestamp", "gen_tech", "gen_energy_source", "DispatchGen_MW", "scenario_name"
+    ]).rename({"DispatchGen_MW": "value"}, axis=1)
+    dispatch = tools.transform.gen_type(dispatch)
+
+    # Sum dispatch across all the projects of the same type and timepoint
+    dispatch = dispatch.groupby(["timestamp", "gen_type"], as_index=False).sum()
+    dispatch = dispatch[dispatch["gen_type"] != "Storage"]
+
+    # Get load dataframe
+    load = tools.get_dataframe("load_balance.csv", usecols=[
+        "timestamp", "zone_demand_mw", "TXPowerNet", "scenario_name"
+    ])
+
+    def process_time(df):
+        df = df.astype({"period": int})
+        df = df[df["period"] == df["period"].max()].drop(columns="period")
+        return df.set_index("datetime")
+
+    # Sum load across all the load zones
+    load = load.groupby(["timestamp"], as_index=False).sum()
+
+    # Include Tx Losses in demand and flip sign
+    load["value"] = (load["zone_demand_mw"] + load["TXPowerNet"]) * -1
+
+    # Rename and convert from wide to long format
+    load = load[["timestamp", "value"]]
+
+    # Add the timestamp information and make period string to ensure it doesn't mess up the graphing
+    dispatch = process_time(tools.transform.timestamp(dispatch))
+    load = process_time(tools.transform.timestamp(load))
+
+    # Convert to TWh (incl. multiply by timepoint duration)
+    dispatch["value"] *= dispatch["tp_duration"] / 1e6
+    load["value"] *= load["tp_duration"] / 1e6
+
+    days = 14
+    freq = str(days) + "D"
+    offset = tools.pd.Timedelta(freq) / 2
+
+    def rolling_sum(df):
+        df = df.rolling(freq, center=True).value.sum().reset_index()
+        df["value"] /= days
+        df = df[(df.datetime.min() + offset < df.datetime) & (df.datetime < df.datetime.max() - offset)]
+        return df
+
+    dispatch = rolling_sum(dispatch.groupby("gen_type", as_index=False))
+    load = rolling_sum(load).set_index("datetime")["value"]
+
+    # Get the state of charge data
+    soc = tools.get_dataframe("StateOfCharge.csv", dtype={"STORAGE_GEN_TPS_1": str}) \
+        .rename(columns={"STORAGE_GEN_TPS_2": "timepoint", "StateOfCharge": "value"})
+    # Sum over all the projects that are in the same scenario with the same timepoint
+    soc = soc.groupby(["timepoint"], as_index=False).sum()
+    soc["value"] /= 1e6  # Convert to TWh
+    max_soc = soc["value"].max()
+
+    # Group by time
+    soc = process_time(tools.transform.timestamp(soc, use_timepoint=True, key_col="timepoint"))
+    soc = soc.rolling(freq, center=True)["value"].mean().reset_index()
+    soc = soc[(soc.datetime.min() + offset < soc.datetime) & (soc.datetime < soc.datetime.max() - offset)]
+    soc = soc.set_index("datetime")["value"]
+
+
+    dispatch = dispatch[dispatch["value"] != 0]
+    dispatch = dispatch.pivot(columns="gen_type", index="datetime", values="value")
+    dispatch = dispatch[dispatch.std().sort_values().index].rename_axis("Technology", axis=1)
+    total_dispatch = dispatch.sum(axis=1)
+
+    max_val = max(total_dispatch.max(), load.max())
+
+    # Scale soc to the graph
+    soc *= max_val / max_soc
+
+    # Plot
+    # Get the colors for the lines
+    # plot
+    ax = tools.get_axes(ylabel="Average Daily Generation (TWh)")
+    ax.set_ylim(0, max_val * 1.05)
+    dispatch.plot(
+        ax=ax,
+        color=tools.get_colors()
+    )
+    soc.plot(ax=ax, color="black", linestyle="dotted")
+    load.plot(ax=ax, color="red", linestyle="dashed")
+    total_dispatch.plot(ax=ax, color="green", linestyle="dashed")
+    ax.fill_between(total_dispatch.index, total_dispatch.values, load.values, alpha=0.2, where=load<total_dispatch, facecolor="green")
+    ax.fill_between(total_dispatch.index, total_dispatch.values, load.values, alpha=0.2, where=load>total_dispatch, facecolor="red")
+
+
+@graph(
+    "dispatch_map",
+    title="Dispatched electricity per load zone"
+)
+def dispatch_map(tools):
+    dispatch = tools.get_dataframe("dispatch_zonal_annual_summary.csv").rename({"Energy_GWh_typical_yr": "value"},
+                                                                               axis=1)
+    dispatch = tools.transform.gen_type(dispatch)
+    dispatch = dispatch.groupby(["gen_type", "gen_load_zone"], as_index=False)["value"].sum()
+    tools.maps.graph_pie_chart(dispatch)

switch_model/generators/core/dispatch.py

@@ -32,6 +32,8 @@
 
 """
 import math
+
+import pandas as pd
 from scipy import fft
 
 from pyomo.environ import *
@@ -422,48 +424,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())
-
-    # 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()
+    num_periods = len(soc["period"].unique())
 
-    # Add the capacity to our dataframe
-    df = df.merge(
-        capacity,
-        on=["period", "scenario_name"],
-        validate="many_to_one",
-        how="left"
-    )
+    # Used later
+    grouped_soc = soc.set_index("datetime").groupby(panel_group, as_index=False)
 
-    # Convert values to TWh
-    df["StateOfCharge"] /= 1e6
-    df["OnlineEnergyCapacityMWh"] /= 1e6
+    # Calculate the weekly SOC
+    weekly_soc = grouped_soc.rolling(rolling_mean_window_size, center=True) \
+        .value \
+        .mean() \
+        .reset_index()
 
-    # Determine information for the label
-    y_axis_lim = df["OnlineEnergyCapacityMWh"].max()
-    offset = y_axis_lim * 0.05
-    df["label_position"] = df["OnlineEnergyCapacityMWh"] + offset
-    df["label"] = df["OnlineEnergyCapacityMWh"].round(decimals=2)
-    label_x_pos = df["datetime"].median()
+    # 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 = 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
+    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 = capacity.value.max()
+    label_x_pos = soc["datetime"].median()
+
+    hlines = pd.concat([capacity, avg_daily_size, max_daily_size])
+
+    # For the max label
+    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="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())
 
 

switch_model/generators/extensions/storage.py

@@ -11,8 +11,8 @@ def add_arguments(parser):
     parser.add_argument("--modules", default=None, nargs='+',
                         help="Modules to load the graphing functions for. "
                              "If not specified reads the modules from modules.txt.")
-    parser.add_argument("-f", "--figure", default=None,
-                        help="Name of the figure to graph. Figure names are the first argument in the @graph() decorator."
+    parser.add_argument("-f", "--figures", default=None, nargs='+',
+                        help="Name of the figures to graph. Figure names are the first argument in the @graph() decorator."
                              " If unspecified graphs all the figures.")
     parser.add_argument("--ignore-modules-txt", default=False, action="store_true",
                         help="When true modules in modules txt are not loaded")
@@ -23,4 +23,4 @@ def graph_scenarios_from_cli(scenarios, args):
         if args.modules is None:
             args.modules = []  # Provide an empty list of modules
 
-    graph_scenarios(scenarios, args.graph_dir, args.overwrite, args.skip_long, args.modules, args.figure)
+    graph_scenarios(scenarios, args.graph_dir, args.overwrite, args.skip_long, args.modules, args.figures)