add ddm calculation

Author: prakaa

plotting/plot_ddm.jl

@@ -0,0 +1,110 @@
+using Dates, DataFrames
+using NEMStorageUnderUncertainty
+using CairoMakie
+using JLD2
+
+function add_legend!(
+    fig::Figure,
+    x_loc,
+    y_loc,
+    value_colors::Vector{PolyElement},
+    value_labels::Vector{String},
+    value_title::String,
+)
+    Legend(
+        fig[y_loc, x_loc],
+        [value_colors],
+        [value_labels],
+        [value_title];
+        framevisible=false,
+    )
+    return nothing
+end
+
+function _makie_plot_across_formulation(
+    fig::Figure,
+    plot_data::DataFrame,
+    position::Int64,
+    mw_capacity::Int64,
+    formulations::Vector{String},
+    colors,
+)
+    title = string(round(Int, mw_capacity)) * " MW"
+    ylabel = L"\frac{\textrm{NegRev}_{\textrm{Actual}} -\textrm{NegRev}_{\textrm{Forecast}}}{\textrm{Revenue}_{\textrm{Actual}}-\textrm{Revenue}_{\textrm{Forecast}}}"
+    lookaheads = unique(plot_data.lookahead)
+    ax = Axis(
+        fig[position, 1];
+        xticks=(1:length(lookaheads), string.(lookaheads) .* " min"),
+        yticks=(range(0, 100; step=20), string.(range(0, 100; step=20))),
+        title,
+        ylabel=ylabel,
+        xlabel="Lookahead (minutes)",
+    )
+    xs = [findfirst(x -> x == lk, lookaheads) for lk in plot_data.lookahead]
+    groups = [findfirst(x -> x == f, formulations) for f in plot_data.label]
+    barplot!(
+        ax,
+        xs,
+        plot_data.ddm;
+        dodge=groups,
+        color=colors[groups],
+        fillto=0.0,
+        strokecolor=:gray,
+        strokewidth=1,
+    )
+    return ylims!(ax, 0.0, 100.0)
+end
+
+function plot_ddm_across_formulations(
+    data_path::String, save_path::String
+)
+    data_file = [f for f in readdir(data_path) if f == "ddm_vpl_vpi.jld2"][]
+    all_data = load(joinpath(data_path, data_file))
+    for (state, value) in pairs(all_data)
+        energy = round(Int, unique(value.energy_capacity)[])
+        state_data = all_data[state]
+        plot_mw_capacities = (25, 100, 400)
+        plot_lookaheads = ("5", "60", "240", "480", "900")
+        fig = Figure(; resolution=(800, 1000))
+        state_data.label = map(
+            x -> NEMStorageUnderUncertainty.formulation_label_map[x], state_data.formulation
+        )
+        state_data.param = replace(state_data.param, missing => "")
+        non_param_formulations = unique(state_data[state_data.param.=="", :formulation])
+        param_formulations = unique(state_data[state_data.param.!="", :formulation])
+        state_data[state_data.param.!="", :label] =
+            state_data[state_data.param.!="", :label] .* " [" .*
+            uppercasefirst.(string.(state_data[state_data.param.!="", :param])) .* "]"
+        sims = sort(unique(state_data.label))
+        colors = NEMStorageUnderUncertainty.generate_formulation_colors(
+            param_formulations, non_param_formulations, sims
+        )
+        for (i, mw_capacity) in enumerate(plot_mw_capacities)
+            energy = round(Int, unique(state_data.energy_capacity)[])
+            mw_df = state_data[state_data.power_capacity.==mw_capacity, :]
+            mw_df = mw_df[mw_df.lookahead.∈(plot_lookaheads,), :]
+            mw_df.lookahead = parse.(Int64, mw_df.lookahead)
+            _makie_plot_across_formulation(fig, mw_df, i, mw_capacity, sims, colors)
+        end
+        f_lg = [PolyElement(; polycolor=colors[i]) for i in 1:length(sims)]
+        add_legend!(fig, 2, :, f_lg, sims, "Simulated formulation")
+        title = "$energy MWh BESS - DDM - $state Prices, 2021"
+        Label(fig[0, :]; text=title, fontsize=22, font="Source Sans Pro")
+        filename = "$(state)_$(energy)_allformulations_ddm.pdf"
+        save(joinpath(save_path, filename), fig; pt_per_unit=1)
+    end
+end
+
+
+data_path = joinpath("results", "data")
+plot_path = joinpath("results", "plots", "neg_rev")
+if !ispath(plot_path)
+    mkpath(plot_path)
+end
+
+@assert ispath(data_path) "Results data not compiled. Run 'make compile_results'."
+
+NEMStorageUnderUncertainty.set_project_plot_theme!()
+plot_ddm_across_formulations(
+    joinpath("results", "data"), plot_path
+)

results/data/vpl_vpi.jld2 results/data/ddm_vpl_vpi.jld2

@@ -105,8 +105,8 @@ function _summarise_simulations(
 end
 
 @doc raw"""
-Calculates values of perfect lookahead and information as absolute values (in AUD) and as
-a percentage of perfect foresight revenue.
+Calculates values of perfect lookahead and information (as absolute values in AUD and as
+a percentage of perfect foresight revenue), and the detrimental decision metric
 
 **Value of perfect lookahead**: What is the additional benefit (revenue) that a participant
 could gain if they were to know exactly what the market prices will be in the *lookahead
@@ -118,6 +118,10 @@ participant could gain if they were to know exactly what the market prices will
 *over the entire year*
   * ``VPI = \textrm{Revenue}_\textrm{Perfect Foresight} -  \textrm{Revenue}_\textrm{Forecast Data Simulation}``
 
+**Detrimental decision metric**: What is the additional negative revenue (i.e. losses)
+incurred when using forecast market prices in the lookahead horizon as a percentage?
+  * ``DDM = \frac{\textrm{NegRev}_{\textrm{Actual Data Simulation}} -\textrm{NegRev}_{\textrm{Forecast Data Simulation}}}{\textrm{Revenue}_{\textrm{Actual Data Simulation}}-\textrm{Revenue}_{\textrm{Forecast Data Simulation}}}``
+
 N.B. This function assumes that the input `df` only has data that corresponds to a device
 of a particular `energy_capacity`.
 
@@ -133,6 +137,7 @@ as a percentage of perfect foresight revenue.
 function calculate_vpl_vpi(df::DataFrame)
     (v_pl_abs, v_pi_abs) = (Float64[], Float64[])
     (v_pl_percentage, v_pi_percentage) = (Float64[], Float64[])
+    ddms = Float64[]
     (power_caps, data) = (Float64[], String[])
     actual_caps = unique(df[df.data_type.=="actual", :power_capacity])
     forecast_caps = unique(df[df.data_type.=="forecast", :power_capacity])
@@ -145,21 +150,25 @@ function calculate_vpl_vpi(df::DataFrame)
             actual_mask = df.data_type .== "actual"
             forecast_mask = df.data_type .== "forecast"
             forecast_rev = df[cap_mask.&forecast_mask.&lk_mask, :revenue]
+            forecast_neg_rev = df[cap_mask.&forecast_mask.&lk_mask, :neg_revenue]
             pf_rev = df[
                 cap_mask.&forecast_mask.&(df.lookahead.=="Perfect Foresight"),
                 :revenue,
             ]
             pi_rev = df[cap_mask.&actual_mask.&lk_mask, :revenue]
+            pi_neg_rev = df[cap_mask.&actual_mask.&lk_mask, :neg_revenue]
             v_pl = pi_rev - forecast_rev
             v_pi = pf_rev - forecast_rev
             v_pl_percentage_pf = @. v_pl / pf_rev * 100
             v_pi_percentage_pf = @. v_pi / pf_rev * 100
+            ddm = @. (pi_neg_rev - forecast_neg_rev) / (pi_rev - forecast_rev) * 100
             push!(power_caps, cap)
             push!(data, lookahead)
             push!(v_pl_abs, v_pl[])
             push!(v_pi_abs, v_pi[])
             push!(v_pl_percentage, v_pl_percentage_pf[])
             push!(v_pi_percentage, v_pi_percentage_pf[])
+            push!(ddms, ddm[])
         end
     end
     return DataFrame(
@@ -171,6 +180,7 @@ function calculate_vpl_vpi(df::DataFrame)
         :vpi_abs => v_pi_abs,
         :vpl_per => v_pl_percentage,
         :vpi_per => v_pi_percentage,
+        :ddm => ddms,
     )
 end
 
@@ -201,7 +211,7 @@ function calculate_vpl_vpi_across_scenarios(summary_folder::String)
             match(r"([A-Z]{2,3})_summary_results.jld2", file).captures[]
         )
 
-        @info "Calculating VPL and VPI for $state"
+        @info "Calculating DDM, VPL and VPI for $state"
         summary_data = load(joinpath(summary_folder, file))
         vpl_vpi_data = DataFrame[]
         for (formulation, summary) in pairs(summary_data)
@@ -215,8 +225,8 @@ function calculate_vpl_vpi_across_scenarios(summary_folder::String)
             push!(vpl_vpi_data, vpl_vpi)
         end
         state_vpl_vpi = vcat(vpl_vpi_data...)
-        @info "Saving VPL and VPI data for $state"
-        jldopen(joinpath(summary_folder, "vpl_vpi.jld2"), "w"; compress=true) do f
+        @info "Saving DDM, VPL and VPI data for $state"
+        jldopen(joinpath(summary_folder, "ddm_vpl_vpi.jld2"), "w"; compress=true) do f
             f["$(state)"] = state_vpl_vpi
         end
     end
@@ -288,9 +298,9 @@ function calculate_summaries_and_vpl_vpi_across_scenarios(sim_folder::String)
             end
         end
     end
-    vpl_vpi_file_name = joinpath(save_path, "vpl_vpi.jld2")
-    if !isfile(vpl_vpi_file_name)
-        @info "Calculating VPL and VPI across scenarios"
+    ddm_vpl_vpi_file_name = joinpath(save_path, "ddm_vpl_vpi.jld2")
+    if !isfile(ddm_vpl_vpi_file_name)
+        @info "Calculating DDM, VPL and VPI across scenarios"
         calculate_vpl_vpi_across_scenarios(save_path)
     end
     return nothing