Add CVAR and NCVAR resource adequacy reporting

Project.toml

@@ -23,8 +23,8 @@ HDF5 = "0.16.15"
 InlineStrings = "1.4.0"
 JSON = "0.21.4"
 LoggingExtras = "1.0.3"
-PRAS = "~0.8.0"
-PRASCore = "~0.8.0"
-PRASFiles = "~0.8.0"
+PRAS = "0.9"
+PRASCore = "0.9"
+PRASFiles = "0.9"
 TimeZones = "^1.13.0"
 julia = "=1.12.1"

cases.csv

@@ -283,6 +283,7 @@ GSw_PRM_StressThresholdLOLD,LOLD threshold [event-days/year]; formulated as Hier
 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,

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=True,
         write_availability_samples=False,
         **kwargs,
     ):
@@ -68,6 +69,7 @@ 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"--iteration={iteration}",
         f"--samples={sw['pras_samples']}",
@@ -160,6 +162,7 @@ def main(t, tnext, casedir, iteration=0):
             write_energy=True,
             write_shortfall_samples=(True if int(sw.GSw_PRM_UpdateMethod) > 1 else False),
         )
+
         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,7 +70,12 @@ 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
@@ -182,6 +187,9 @@ function run_pras(pras_system_path::String, args::Dict)
     if args["write_energy"] == 1
         resultspec["energy"] = PRAS.StorageEnergy()
     end
+    ## Hourly per-sample shortfall (region x hour x sample) is memory-heavy, so only
+    ## request it when the hourly samples are written. Per-sample totals (used for CVaR)
+    ## come from the regular Shortfall result (PRAS >= 0.9).
     if args["write_shortfall_samples"] == 1
         resultspec["short_samples"] = PRAS.ShortfallSamples()
     end
@@ -285,6 +293,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 +311,44 @@ 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
+        ## Read the per-sample totals (MWh) from the regular Shortfall result (PRAS >= 0.9).
+        sf = results["short"]
+        ## Use filtered system regions to avoid DC converter pseudo-regions without load.
+        region_names = sf.regions.names
+        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.shortfall_samples)
+            for r in regions
+                i_r = findfirst(==(r), region_names)
+                f["$r", compress=4] = Float64.(sf.shortfall_region_samples[i_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 +479,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,
@@ -485,4 +508,4 @@ if abspath(PROGRAM_FILE) == @__FILE__
     main(args)
 
     #%%
-end
+end

reeds/resource_adequacy/stress_periods.py

@@ -143,6 +143,51 @@ def calc_neue(dfeue_agg, dfload_agg):
     return neue
 
 
+def get_shortfall_totals_by_sample(case, t, iteration=0):
+    filepath = os.path.join(case, 'handoff', 'PRAS', f'PRAS_{t}i{iteration}-shortfall_totals_by_sample.h5')
+    if not os.path.isfile(filepath):
+        raise FileNotFoundError(f"{filepath} not found. Re-run PRAS with --write_shortfall_samples_totals 1.")
+    df = reeds.io.read_pras_results(filepath)
+    df.columns = df.columns.astype(str)
+    if 'sample' in df.columns:
+        df = df.set_index('sample')
+    elif df.index.name != 'sample':
+        df.index = pd.RangeIndex(1, len(df) + 1, name='sample')
+    return df.apply(pd.to_numeric, errors='coerce').clip(lower=0)
+
+
+def _sample_cvar(samples, alpha=0.95):
+    x = pd.Series(samples).dropna().astype(float)
+    if x.empty:
+        return np.nan
+
+    # Round before applying ceil to remove floating-point noise.
+    # For example, a mathematically exact tail size of 50 may be
+    # represented as 50.00000000000004, which would otherwise select 51 samples.
+    tail_size = round((1 - alpha) * len(x), 12)
+    n_tail = max(1, int(np.ceil(tail_size)))
+
+    return x.nlargest(n_tail).mean()
+
+def calc_cvar(shortfall_samples_agg, alpha=0.95):
+    """
+    CVAR from total shortfall by PRAS sample.
+    """
+    cvar = shortfall_samples_agg.apply(
+        lambda s: _sample_cvar(s, alpha=alpha),
+        axis=0,
+    )
+    return cvar
+
+
+def calc_ncvar(cvar, dfload_agg):
+    """
+    NCVAR = CVAR / total load, in ppm.
+    """
+    ncvar = cvar / dfload_agg.sum().reindex(cvar.index) * 1e6
+    return ncvar
+
+
 def calc_peak_eue(dfeue_agg, dfload_agg, norm:Literal['peak','hourly','absolute']='peak'):
     """
     Get the peak hourly outage magnitude
@@ -192,6 +237,12 @@ def calc_ra_metrics(
     sw = reeds.io.get_switches(case)
     numyears = len(sw.resource_adequacy_years_list)
 
+    ### Get total shortfall samples for CVAR calculation
+    shortfall_samples = (
+        get_shortfall_totals_by_sample(case=case, t=t, iteration=iteration)
+        .drop(columns=['USA'], errors='ignore')
+    )
+
     ### Loop over aggregation levels and calculate all metrics
     ra_metrics = {}
     for level in levels:
@@ -212,6 +263,18 @@ def calc_ra_metrics(
         ra_metrics[level, 'euemax_peakloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'peak')
         ra_metrics[level, 'euemax_hourlyloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'hourly')
         ra_metrics[level, 'euemax_mw'] = calc_peak_eue(dfeue_agg, dfload_agg, 'absolute')
+        ## Calculate tail-based metrics (CVAR and NCVAR)
+        regions = [c for c in shortfall_samples.columns if c in rmap.index]
+        if len(regions):
+            shortfall_samples_agg = (
+                shortfall_samples[regions]
+                .rename(columns=rmap)
+                .T.groupby(level=0)
+                .sum().T
+            )
+            cvar = calc_cvar(shortfall_samples_agg, alpha=float(sw.GSw_PRM_CVARalpha))
+            ra_metrics[level, 'cvar_mwh_peryear'] = cvar / numyears
+            ra_metrics[level, 'ncvar_ppm'] = calc_ncvar(cvar, dfload_agg)
 
     ### Combine it
     dfout = pd.concat(ra_metrics, names=['level','metric','region']).rename('value')
@@ -511,13 +574,12 @@ def get_stress_periods(case, sw, t, iteration):
         for i,row in stressperiods_this_iteration.iterrows()
     ]
     covered_hours = [i for sublist in covered_hours for i in sublist]
-
-    ### Check all stress criteria; for regions that fail, add new stress periods
     _failed = {}
     _high_stress_periods = {}
     _shoulder_periods = {}
 
     stress_metrics = [i.lower() for i in sw.GSw_PRM_StressThresholdMetrics.split('/')]
+
     for stress_metric in stress_metrics:
         switch = RA_SWITCHES[stress_metric]
         for criterion in sw[switch].split('/'):

runreeds.py

@@ -330,6 +330,11 @@ def check_compatibility(sw):
                     f"but '{stress_value}' was provided"
                 )
 
+    ## CVAR value in [0,1)
+    alpha = float(sw['GSw_PRM_CVARalpha'])
+    if not (0 <= alpha < 1):
+        raise ValueError(f"GSw_PRM_CVARalpha must be in [0, 1). Got {alpha}")
+
     ### GSw_PRM_UpdateMethod 1-3 (static or PRAS-informed PRM update) is computed from the
     ### NEUE-based shortfall, so it requires NEUE to be an active stress metric
     if int(sw['GSw_PRM_UpdateMethod']) in [1, 2, 3] and 'neue' not in used_metrics: