Merge pull request #103 from staadecker/plots

Plot improvements

Author: pesap

papers/Martin_Staadecker_Value_of_LDES_and_Factors/LDES_paper_graphs/plots.py

@@ -287,11 +287,12 @@ def graph_energy_balance(tools):
     )
     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

@@ -915,6 +915,33 @@ 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

@@ -922,6 +922,146 @@ 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 *
@@ -502,58 +504,91 @@ def post_solve(instance, outdir):
     "state_of_charge",
     title="State of Charge Throughout the Year",
     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())
+
+    # Used later
+    grouped_soc = soc.set_index("datetime").groupby(panel_group, as_index=False)
+
+    # Calculate the weekly SOC
+    weekly_soc = (
+        grouped_soc.rolling(rolling_mean_window_size, center=True)
+        .value.mean()
+        .reset_index()
+    )
 
     # 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()
+    capacity = (
+        capacity.groupby(panel_group, as_index=False)["OnlineEnergyCapacityMWh"]
+        .sum()
+        .rename({"OnlineEnergyCapacityMWh": "value"}, axis=1)
+    )
+    capacity.value /= 1e6
+    capacity["type"] = "Total Energy Capacity"
 
-    # Add the capacity to our dataframe
-    df = df.merge(
-        capacity, on=["period", "scenario_name"], validate="many_to_one", how="left"
+    # 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"
 
-    # Convert values to TWh
-    df["StateOfCharge"] /= 1e6
-    df["OnlineEnergyCapacityMWh"] /= 1e6
+    # 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])
 
-    # 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()
+    # 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()
+        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"),
+            pn.aes(yintercept="value", label="label", color="type"),
+            data=hlines,
             linetype="dashed",
-            color="blue",
         )
         + pn.geom_text(
-            pn.aes(label="label", x=label_x_pos, y="label_position"),
+            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

@@ -29,9 +29,10 @@ def add_arguments(parser):
     )
     parser.add_argument(
         "-f",
-        "--figure",
+        "--figures",
         default=None,
-        help="Name of the figure to graph. Figure names are the first argument in the @graph() decorator."
+        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(
@@ -53,5 +54,5 @@ def graph_scenarios_from_cli(scenarios, args):
         args.overwrite,
         args.skip_long,
         args.modules,
-        args.figure,
+        args.figures,
     )