Add CVAR and NCVAR resource adequacy reporting

cases.csv

@@ -273,13 +273,14 @@ GSw_PRM_StressOutages,Whether to apply the availability factor (forced + schedul
 GSw_PRM_StressSeedLoadLevel,Region hierarchy level at which to include peak coincident load days as seeded stress periods (or False to ignore peaks),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,transgrp,
 GSw_PRM_StressSeedMinRElevel,Region hierarchy level at which to include minimum wind and solar capacity factor days as seeded stress periods (or False to ignore min-wind/solar CF days),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,interconnect,
 GSw_PRM_StressStorageCutoff,How to select shoulder stress periods giving storage time to recharge before/after high-unserved-energy periods. Two-part switch separated by _. The first argument is 'EUE' or 'capacity' or 'absolute'. If 'EUE' the second argument specifies a PRAS storage headspace / EUE threshold [fraction]; if 'cap' it specifies a headspace / storage capacity threshold [fraction]; if 'abs' it specifies absolute number of periods before/after [integer]. Turned off if set to 'off'.,N/A,EUE_0.1,
-GSw_PRM_StressThresholdMetrics,/-delimited list of metrics for identifying stress periods; supported metrics (case-insensitive) are: LOLD | LOLE | LOLH | NEUE | Depth | Duration,N/A,NEUE,
+GSw_PRM_StressThresholdMetrics,/-delimited list of resource adequacy metrics. LOLD | LOLE | LOLH | NEUE | Depth | Duration identify stress periods; CVAR and NCVAR are report-only metrics and do not add stress periods or update PRM.,N/A,NEUE,
 GSw_PRM_StressThresholdDepth,Outage depth threshold [MW_EUE/MW_peak_load] (fraction); formulated as HierarchyLevel_Depth where HierarchyLevel is a column in hierarchy.csv; Depth is the max outage magnitude in [MW EUE] / [MW peak demand],N/A,transgrp_0.1,
 GSw_PRM_StressThresholdDuration,Outage duration threshold [hours]; formulated as HierarchyLevel_Duration where HierarchyLevel is a column in hierarchy.csv; Duration is the max outage duration in hours,N/A,transgrp_12,
 GSw_PRM_StressThresholdLOLD,LOLD threshold [event-days/year]; formulated as HierarchyLevel_LOLD where HierarchyLevel is a column in hierarchy.csv; LOLD is loss-of-load event-days per year,N/A,transgrp_0.1,
 GSw_PRM_StressThresholdLOLE,LOLE threshold [events/year]; formulated as HierarchyLevel_LOLE where HierarchyLevel is a column in hierarchy.csv; events is loss-of-load events per year,N/A,transgrp_0.1,
 GSw_PRM_StressThresholdLOLH,LOLH threshold [event-hours/year]; formulated as HierarchyLevel_LOLH where HierarchyLevel is a column in hierarchy.csv; LOLH is loss-of-load event-hours per year,N/A,transgrp_2.4,
 GSw_PRM_StressThresholdNEUE,NEUE threshold [ppm]; formulated as HierarchyLevel_NEUE where HierarchyLevel is a column in hierarchy.csv; NEUEppm is normalized expected unserved energy in parts per million,N/A,transgrp_1,
+GSw_PRM_CVARAlpha,Alpha value used for CVAR/NCVAR calculations,float,0.95,
 GSw_PRM_UpdateFraction,Fraction to add to the PRM if a region fails RA threshold (only used if GSw_PRM_UpdateMethod = 1),float,0.02,
 GSw_PRM_UpdateMethod,Option to update PRM: (0) no update; (1) static update set by GSw_PRM_UpdateFraction; (2) dynamic update informed by PRAS; (3) dynamic update but only after all new stress periods have been added,0; 1; 2; 3,0,
 GSw_PRMTRADE_level,hierarchy level within which to allow PRM trading,r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region,country,

cases_test.csv

@@ -3,7 +3,7 @@ ignore,1,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
 GSw_Region,cendiv/Pacific,,country/USA,country/USA,country/USA,,st/ID.WY.NE.IA.IL,st/MA,,,,,,,st/WY,interconnect/western,transreg/PJM,st/NY.VT,st/OR,st/NE.NY.PA,st/NE.NY.PA,st/ID.WY.NE.IA.IL,st/KS,country/USA,country/USA,st/NY.NJ,,,,st/MA.RI.CT.NY.NJ.PA.OH,
 endyear,2032,,2050,2050,2050,2029,2060,2026,,,,,,,,,,,2035,2030,2030,2060,2035,2050,2050,2050,,,,,
 yearset,,,,,,,2010..2060..10,,,,,,,,,,,,,2010..2050..5,2010..2050..5,2010..2060..10,,,2010_2025_2050,2010..2050..5,,,,,
-GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,z134,z54,z48,,,,,,
+GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,z134,z54,z48,z132,,,,,
 GSw_GasCurve,2,,1,1,,,,,,,,,,,,,,,,,,,,1,1,,,,,,
 GSw_Geothermal,,,,2,,,,,,,,,,,,,,,,,,,0,,0,,,,,,
 GSw_GrowthPenalties,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,,
@@ -51,7 +51,7 @@ GSw_StartCost,,,,,,,,,,,,,,,,,,,,,,,0,0,0,,,,,,
 GSw_H2,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,,
 GSw_H2_PTC,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,,
 GSw_H2Combustion,,,,,,,,,,,,,,,,,,,,,,,,,0,,,,,,
-GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth,,,,,
+GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth/CVAR/NCVAR,,,,,
 GSw_DRShed,,,,,,,,,,,,,,,,,,,,,,,,,,,1,,,,
 GSw_MGA_CostDelta,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.01,,,
 GSw_LoadSiteCF,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.95,,

reeds/resource_adequacy/ra_calcs.py

@@ -24,6 +24,7 @@ def run_pras(
         write_surplus=False,
         write_energy=False,
         write_shortfall_samples=False,
+        write_shortfall_samples_totals=False,
         write_availability_samples=False,
         **kwargs,
     ):
@@ -68,7 +69,9 @@ def run_pras(
         f"--write_surplus={int(write_surplus)}",
         f"--write_energy={int(write_energy)}",
         f"--write_shortfall_samples={int(write_shortfall_samples)}",
+        f"--write_shortfall_samples_totals={int(write_shortfall_samples_totals)}",
         f"--write_availability_samples={int(write_availability_samples)}",
+        f"--cvar_alpha={float(sw['GSw_PRM_CVARAlpha'])}",
         f"--iteration={iteration}",
         f"--samples={sw['pras_samples']}",
         f"--overwrite={int(overwrite)}",
@@ -153,13 +156,25 @@ def main(t, tnext, casedir, iteration=0):
         2: True,
     }[int(sw['pras'])]
     if pras_this_solve_year or int(sw.GSw_PRM_StressIterateMax):
+        stress_metrics = [
+            m.strip().upper()
+            for m in sw.GSw_PRM_StressThresholdMetrics.split('/')
+            if m.strip()
+        ]
+        write_shortfall_samples_totals = any(
+            metric in {'CVAR', 'NCVAR'}
+            for metric in stress_metrics
+        )
+
         result = run_pras(
             casedir, t,
             iteration=iteration,
             write_flow=(True if t == max(solveyears) else False),
             write_energy=True,
             write_shortfall_samples=(True if int(sw.GSw_PRM_UpdateMethod) > 1 else False),
+            write_shortfall_samples_totals=write_shortfall_samples_totals,
         )
+
         if result.returncode:
             raise Exception(
                 f"run_pras.jl returned code {result.returncode}. Check gamslog.txt for error trace."

reeds/resource_adequacy/run_pras.jl

@@ -70,10 +70,20 @@ function parse_commandline()
             default = 0
             required = false
         "--write_shortfall_samples"
-            help = "Write the sample-level shortfall"
+            help = "Write per-sample hourly shortfall by region"
             arg_type = Int
             default = 0
             required = false
+        "--write_shortfall_samples_totals"
+            help = "Write per-sample total shortfall by region over the full PRAS time period"
+            arg_type = Int
+            default = 0
+            required = false
+        "--cvar_alpha"
+            help = "Alpha for CVaR (e.g., 0.95)"
+            arg_type = Float64
+            default = 0.95
+            required = false
         "--write_availability_samples"
             help = "Write the sample-level generator and storage availability"
             arg_type = Int
@@ -182,7 +192,7 @@ function run_pras(pras_system_path::String, args::Dict)
     if args["write_energy"] == 1
         resultspec["energy"] = PRAS.StorageEnergy()
     end
-    if args["write_shortfall_samples"] == 1
+    if args["write_shortfall_samples"] == 1 || args["write_shortfall_samples_totals"] == 1
         resultspec["short_samples"] = PRAS.ShortfallSamples()
     end
     if args["write_availability_samples"] == 1
@@ -208,6 +218,38 @@ function run_pras(pras_system_path::String, args::Dict)
     @info "$(PRAS.EUE(results["short"]))"
     @info "$(PRAS.NEUE(results["short"]))"
 
+    #%% Print CVAR and NCVAR for the entire modeled region
+    if args["write_shortfall_samples_totals"] == 1 && haskey(results, "short_samples")
+        _, _, _, _, energyunit = PRAS.get_params(sys)
+        alpha = Float64(args["cvar_alpha"])
+
+        cvar_result = PRAS.CVAR(
+            energyunit,
+            results["short_samples"],
+            alpha,
+        )
+
+        cvar_value = PRAS.val(cvar_result.cvar)
+        cvar_stderr = PRAS.stderror(cvar_result.cvar)
+        cvar_var = cvar_result.var
+
+        ### Normalize CVaR by total load over the full PRAS time period.
+        total_load = sum(sys.regions.load)
+
+        ncvar_value = cvar_value / total_load * 1e6
+        ncvar_stderr = cvar_stderr / total_load * 1e6
+        ncvar_var = cvar_var / total_load * 1e6
+
+        @info(
+            "CVAR = $(cvar_value)±$(cvar_stderr) MWh; " *
+            "VaR = $(cvar_var) MWh; alpha = $(alpha)"
+        )
+        @info(
+            "NCVAR = $(ncvar_value)±$(ncvar_stderr) ppm; " *
+            "VaR = $(ncvar_var) ppm; alpha = $(alpha)"
+        )
+    end
+
     ## Filter out DC regions used for VSC HVDC transmission
     regions = [r for r in sys.regions.names if !(occursin("|", r))]
 
@@ -285,6 +327,7 @@ function run_pras(pras_system_path::String, args::Dict)
         end
         @info("Wrote PRAS surplus to $(surplusfile)")
     end
+
     ### Storage energy
     if args["write_energy"] == 1
         dfenergy = DF.DataFrame()
@@ -302,31 +345,41 @@ function run_pras(pras_system_path::String, args::Dict)
         @info("Wrote PRAS storage energy to $(energyfile)")
     end
 
-    ### Sample-level shortfall
+    ### Per-sample hourly shortfall by region
     if args["write_shortfall_samples"] == 1
-        dictshort = Dict(s => DF.DataFrame() for s = 1:args["samples"])
-        for s in range(1, args["samples"])
-            dictshort[s] = DF.DataFrame(
-                transpose(getindex.(results["short_samples"][:, :], s)),
-                sys.regions.names
-            )
-            # subset to regions (filter out DC regions) 
-            dictshort[s] = dictshort[s][:,findall(regions .∈ Ref(sys.regions.names))]
-        end
-        ## Write it
+        sf = results["short_samples"]
+        ## Use filtered system regions to avoid DC converter pseudo-regions without load.
+        region_names = sf.regions.names
+        idx = [findfirst(==(r), region_names) for r in regions]
+
         shortfile = replace(outfile, ".h5"=>"-shortfall_samples.h5")
         HDF5.h5open(shortfile, "w") do f
-            ## Create a group for each sample. Within each group, write an array for each region.
-            for s in range(1, args["samples"])
+            for s in 1:args["samples"]
                 HDF5.create_group(f, "$s")
-                for column in DF._names(dictshort[s])
-                    f["$s"]["$column", compress=4] = convert(Array, dictshort[s][!, column])
+                for (r, i_r) in zip(regions, idx)
+                    arr = Float64.(sf.shortfall[i_r, :, s])
+                    f["$s"]["$r", compress=4] = arr
                 end
             end
         end
         @info("Wrote PRAS shortfall by sample to $(shortfile)")
     end
 
+    ### Per-sample total shortfall by region over the full PRAS time period, needed for CVaR
+    if args["write_shortfall_samples_totals"] == 1
+        sf = results["short_samples"]
+        ## Use filtered system regions to avoid DC converter pseudo-regions without load.
+        totalsfile = replace(outfile, ".h5"=>"-shortfall_totals_by_sample.h5")
+        HDF5.h5open(totalsfile, "w") do f
+            f["sample", compress=4] = collect(1:args["samples"])
+            f["USA", compress=4] = Float64.(sf[])
+            for r in regions
+                f["$r", compress=4] = Float64.(sf[r])
+            end
+        end
+        @info("Wrote PRAS shortfall totals by sample to $(totalsfile)")
+    end
+
     ### Sample-level generator and storage availability
     if args["write_availability_samples"] == 1
         dictavail = Dict(s => DF.DataFrame() for s = 1:args["samples"])
@@ -457,6 +510,7 @@ if abspath(PROGRAM_FILE) == @__FILE__
     #     "write_surplus" => 0,
     #     "write_energy" => 0,
     #     "write_shortfall_samples" => 1,
+    #     "write_shortfall_samples_totals" => 1,
     #     "write_availability_samples" => 0,
     #     "overwrite" => 1,
     #     "debug" => 0,
@@ -466,6 +520,7 @@ if abspath(PROGRAM_FILE) == @__FILE__
     #     "pras_existing_unit_size" => 1,
     #     "pras_max_unitsize_prm" => 1,
     #     "pras_seed" => 1,
+    #     "cvar_alpha" => 0.95,
     # )
     # reedscase = args["reedscase"]
     # solve_year = args["solve_year"]
@@ -485,4 +540,4 @@ if abspath(PROGRAM_FILE) == @__FILE__
     main(args)
 
     #%%
-end
+end