Show balance slacks owned by buses and areas (#94)

src/call_plots.jl

@@ -1006,7 +1006,7 @@ plot = plot_fuel(res)
 
 # Accepted Key Words
 - `generator_mapping_file` = "file_path" : file path to yaml defining generator category by fuel and primemover
-- `slacks::Bool = true` : display slack variables
+- `slacks::Bool = true` : display the system balance slack variables as "Unserved Energy" and "Over Generation". Nodal and area formulations attach one slack per bus/area; those are summed into a single series per direction. Reactive slacks (the `"Q"` meta of full AC formulations) are excluded, since this is an active-power plot.
 - `load::Bool = true` : display load line
 - `curtailment::Bool = true`: To plot the curtailment in the stack plot
 - `storage::Bool = true`: include storage components (as "<category> In"/"<category> Out" traces)
@@ -1041,8 +1041,8 @@ _report_plot_fuel(backend::PlottingBackend, result; kwargs...) =
 # contributing component must vanish instead of producing all-zero columns.
 
 # TODO upstream: PowerAnalytics has no built-in metrics for these entry types
-# (it should export `calc_system_slack_down` and forecast metrics for the
-# storage/source time-series parameters); build them locally until then. See
+# (it should export forecast metrics for the storage/source time-series
+# parameters); build them locally until then. See
 # https://github.com/PabloBotin/PowerAnalytics.jl/issues/4.
 const _CALC_POWER_OUTPUT =
     PA.make_component_metric_from_entry("PowerOutput", PA.PSI.PowerOutput)
@@ -1054,8 +1054,6 @@ const _CALC_ACTIVE_POWER_OUT_FORECAST = PA.make_component_metric_from_entry(
     "ActivePowerOutForecast",
     PA.PSI.ActivePowerOutTimeSeriesParameter,
 )
-const _CALC_SYSTEM_SLACK_DOWN =
-    PA.make_system_metric_from_entry("SystemSlackDown", PA.PSI.SystemBalanceSlackDown)
 
 # Fallback chain for generators: dispatch power if the component was modeled
 # with a variable, otherwise its forecast parameter (e.g. `FixedOutput`
@@ -1081,12 +1079,10 @@ const _SOURCE_IN_METRICS =
     ((PA.Metrics.calc_active_power_in, -1.0), (_CALC_ACTIVE_POWER_IN_FORECAST, 1.0))
 const _SOURCE_OUT_METRICS =
     ((PA.Metrics.calc_active_power_out, 1.0), (_CALC_ACTIVE_POWER_OUT_FORECAST, 1.0))
-# System balance slacks and their fixed display names, taken from
-# `PA.BALANCE_SLACKVARS` so the naming has a single source of truth.
-const _SLACK_METRICS = (
-    (PA.BALANCE_SLACKVARS[PA.PSI.SystemBalanceSlackUp], PA.Metrics.calc_system_slack_up),
-    (PA.BALANCE_SLACKVARS[PA.PSI.SystemBalanceSlackDown], _CALC_SYSTEM_SLACK_DOWN),
-)
+# System balance slack entry types, in the display order they stack. The names
+# come from `PA.BALANCE_SLACKVARS` so the naming has a single source of truth.
+const _SLACK_ENTRY_TYPES =
+    (PA.PSI.SystemBalanceSlackUp, PA.PSI.SystemBalanceSlackDown)
 
 # Catch-all category for components matched by no rule in the generator
 # mapping; matches the `Other` key in the default mapping and the color
@@ -1219,21 +1215,47 @@ function _accumulate_curtailment!(
     return acc
 end
 
+"""
+Every variable key holding a slack of entry type `T`, whatever component type
+owns it. PowerSimulations attaches the balance slacks to `PSY.System` under
+`CopperPlatePowerModel`/`PTDFPowerModel`, to `PSY.Area` under the area models,
+and to `PSY.ACBus` — one column per bus — under the power-flow models, so
+looking only for the `PSY.System` variant (as PowerAnalytics'
+`calc_system_slack_up` does) makes nodal and area slacks disappear. Full AC
+models split the bus slacks into `"P"` and `"Q"` metas; only `"P"` is returned,
+because the fuel stack is an active-power plot.
+"""
+function _slack_keys(result::IS.Results, ::Type{T}) where {T <: PA.PSI.VariableType}
+    return filter(PA.PSI.list_variable_keys(result)) do key
+        PA.PSI.get_entry_type(key) === T && key.meta != "Q"
+    end
+end
+
 function _accumulate_slacks!(acc::_FuelAccumulator, result::IS.Results)
-    for (name, metric) in _SLACK_METRICS
-        df = try
-            PA.compute(metric, result)
-        catch e
-            # Results without slack variables simply skip the category.
-            _is_missing_result_error(e) || rethrow()
-            continue
-        end
-        _add_fuel_values!(
-            acc,
-            name,
-            Vector{Dates.DateTime}(PA.get_time_vec(df)),
-            Vector{Float64}(PA.get_data_vec(df)),
+    for entry in _SLACK_ENTRY_TYPES
+        slack_keys = _slack_keys(result, entry)
+        # Results without slack variables simply skip the category.
+        isempty(slack_keys) && continue
+        name = PA.BALANCE_SLACKVARS[entry]
+        dfs = PA.PSI.read_results_with_keys(
+            result,
+            slack_keys;
+            table_format = IS.TableFormat.WIDE,
         )
+        for df in values(dfs)
+            # One column per bus/area (exactly one for the `PSY.System` case);
+            # the whole direction stacks as a single series.
+            vals = zeros(Float64, DataFrames.nrow(df))
+            for col in DataFrames.names(df, DataFrames.Not(PA.DATETIME_COL))
+                vals .+= df[!, col]
+            end
+            _add_fuel_values!(
+                acc,
+                name,
+                Vector{Dates.DateTime}(df[!, PA.DATETIME_COL]),
+                vals,
+            )
+        end
     end
     return acc
 end
@@ -1623,7 +1645,7 @@ renderer.
 
 # Accepted Key Words
 - `generator_mapping_file` = "file_path" : file path to yaml defining generator category by fuel and primemover
-- `slacks::Bool = true` : display slack variables
+- `slacks::Bool = true` : display the system balance slack variables as "Unserved Energy" and "Over Generation". Nodal and area formulations attach one slack per bus/area; those are summed into a single series per direction. Reactive slacks (the `"Q"` meta of full AC formulations) are excluded, since this is an active-power plot.
 - `load::Bool = true` : display load line
 - `curtailment::Bool = true`: To plot the curtailment in the stack plot
 - `storage::Bool = true`: include storage components (as "<category> In"/"<category> Out" traces)

test/test_fuel_stack_behavior.jl

@@ -320,3 +320,104 @@ end
     @info("removing test files")
     rm(save_root; recursive = true)
 end
+
+# --- Balance slacks owned by something other than `PSY.System` (issue #94) ---
+
+# PowerSimulations attaches the balance slacks to the component type implied by
+# the network formulation, so a nodal formulation stores one slack column per
+# `ACBus`. The load is scaled up so the slacks are actually nonzero and the
+# aggregation assertions have teeth.
+function run_nodal_slack_model()
+    sys = deepcopy(PSB.build_system(PSB.PSITestSystems, "c_sys5_uc"))
+    for load in get_components(PowerLoad, sys)
+        set_max_active_power!(load, 3 * get_max_active_power(load))
+    end
+    template = ProblemTemplate(NetworkModel(DCPPowerModel; use_slacks = true))
+    set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment)
+    set_device_model!(template, PowerLoad, StaticPowerLoad)
+    set_device_model!(template, Line, StaticBranch)
+    model = DecisionModel(
+        template,
+        sys;
+        optimizer = optimizer_with_attributes(HiGHS.Optimizer),
+        horizon = Hour(6),
+    )
+    build!(model; output_dir = mktempdir())
+    solve!(model)
+    return OptimizationProblemResults(model)
+end
+
+@testset "bus-level balance slacks appear in the fuel stack" begin
+    res = run_nodal_slack_model()
+
+    # The regression of #94: the slacks are keyed on `ACBus`, not `System`, so
+    # looking only for the `System` variant made them vanish from the plot.
+    @test Set(
+        PSI.encode_key_as_string(k) for k in PSI.list_variable_keys(res) if
+        PSI.get_entry_type(k) in keys(PA.BALANCE_SLACKVARS)
+    ) == Set(["SystemBalanceSlackUp__ACBus", "SystemBalanceSlackDown__ACBus"])
+
+    p = plot_fuel_plotly(res; set_display = false, stack = true, auto_units = false)
+    trace_names = [t.name for t in p.data]
+    @test "Unserved Energy" in trace_names
+    @test "Over Generation" in trace_names
+
+    # Each direction is the row-wise sum over the per-bus columns.
+    for (name, entry) in (
+        ("Unserved Energy", PSI.SystemBalanceSlackUp),
+        ("Over Generation", PSI.SystemBalanceSlackDown),
+    )
+        entry_keys =
+            [k for k in PSI.list_variable_keys(res) if PSI.get_entry_type(k) == entry]
+        df = only(
+            values(
+                PSI.read_results_with_keys(
+                    res,
+                    entry_keys;
+                    table_format = IS.TableFormat.WIDE,
+                ),
+            ),
+        )
+        # More than DateTime plus one column, i.e. genuinely nodal.
+        @test ncol(df) > 2
+        expected = vec(sum(Matrix(no_datetime(df)); dims = 2))
+        @test collect(only([t for t in p.data if t.name == name]).y) ≈ expected
+    end
+    # The scaled-up load leaves energy unserved, so the sum above is not
+    # trivially zero.
+    @test sum(only([t for t in p.data if t.name == "Unserved Energy"]).y) > 0
+
+    p_cm = plot_fuel(res; set_display = false, stack = true)
+    @test p_cm.series_count == length(p.data)
+end
+
+@testset "system-level balance slacks are unchanged" begin
+    (_, results_ed) = run_test_sim(TEST_RESULT_DIR, TEST_SIM_NAME)
+
+    # The ED template is CopperPlate with `use_slacks = true`, so the slacks are
+    # keyed on `PSY.System` — the only case PowerAnalytics' own system metrics
+    # handle. Those values are the pre-fix reference and must be reproduced.
+    calc_slack_down =
+        PA.make_system_metric_from_entry("SystemSlackDown", PSI.SystemBalanceSlackDown)
+    p = plot_fuel_plotly(results_ed; set_display = false, stack = true, auto_units = false)
+    for (name, metric) in (
+        ("Unserved Energy", PA.Metrics.calc_system_slack_up),
+        ("Over Generation", calc_slack_down),
+    )
+        expected = PA.get_data_vec(PA.compute(metric, results_ed))
+        @test collect(only([t for t in p.data if t.name == name]).y) ≈ expected
+    end
+end
+
+@testset "results without balance slacks skip the slack categories" begin
+    (results_uc, _) = run_test_sim(TEST_RESULT_DIR, TEST_SIM_NAME)
+
+    # The UC template runs with `use_slacks = false`: no slack variable is
+    # stored, so the categories must be absent rather than raising.
+    @test !any(
+        PSI.get_entry_type(k) in keys(PA.BALANCE_SLACKVARS) for
+        k in PSI.list_variable_keys(results_uc)
+    )
+    p = plot_fuel_plotly(results_uc; set_display = false, stack = true)
+    @test isdisjoint([t.name for t in p.data], ["Unserved Energy", "Over Generation"])
+end