Bc cvar updates

cases.csv

@@ -270,9 +270,12 @@ 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 options are: NEUE, LOLH)",N/A,NEUE,
-GSw_PRM_StressThresholdLOLH,LOLH threshold [hours/year] above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_LOLH_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; LOLH is loss of load event-hours per year [event-h/year]; PeriodAggMethod is 'sum' or 'max' over the hours in each period (only used in period selection) (see README.md for detailed notes),N/A,transgrp_2.4_sum,  
-GSw_PRM_StressThresholdNEUE,Annual NEUE level [ppm] threshold above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_NEUEppm_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; NEUEppm is normalized expected unserved energy in parts per million [ppm]; PeriodAggMethod is 'sum' or 'max' over the hours in each period (only used in period selection) (see README.md for detailed notes),N/A,transgrp_1_sum,  
+GSw_PRM_StressThresholdMetrics,"/-delimited list of resource adequacy metrics. NEUE and LOLH identify stress periods; CVAR and NCVAR are check-only metrics and do not add stress periods or update PRM.",N/A,NEUE,
+GSw_PRM_StressThresholdLOLH,LOLH threshold [hours/year] above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_LOLH_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; LOLH is loss of load event-hours per year [event-h/year]; PeriodAggMethod is 'sum' or 'max' over the hours in each period (only used in period selection) (see README.md for detailed notes),N/A,transgrp_2.4_sum,
+GSw_PRM_StressThresholdNEUE,Annual NEUE level [ppm] threshold above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_NEUEppm_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; NEUEppm is normalized expected unserved energy in parts per million [ppm]; PeriodAggMethod is 'sum' or 'max' over the hours in each period (only used in period selection) (see README.md for detailed notes),N/A,transgrp_1_sum,
+GSw_PRM_CVARAlpha,Alpha value used for CVAR/NCVAR calculations,float,0.95,
+GSw_PRM_StressThresholdNCVAR,Annual NCVAR level [ppm] used only for NCVAR target checks; formulated as HierarchyLevel_NCVARppm_NCVAR_cvar where HierarchyLevel is a column in hierarchy.csv,N/A,transgrp_10_NCVAR_cvar,
+GSw_PRM_StressThresholdCVAR,Annual CVAR level [MWh/year] used only for CVAR target checks; formulated as HierarchyLevel_CVARmwh_CVAR_cvar where HierarchyLevel is a column in hierarchy.csv,N/A,country_50000_CVAR_cvar,
 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=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(m in {'CVAR', 'NCVAR'} for m in stress_metrics)
+            or int(sw.GSw_PRM_UpdateMethod) > 1
+        )
+
         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

@@ -1,5 +1,6 @@
 #%% Imports
 import ArgParse
+import CSV
 import DataFrames
 import Logging
 import LoggingExtras
@@ -70,10 +71,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 +193,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 +219,41 @@ function run_pras(pras_system_path::String, args::Dict)
     @info "$(PRAS.EUE(results["short"]))"
     @info "$(PRAS.NEUE(results["short"]))"
 
+    # CVAR results
+    if haskey(results, "short_samples")
+        _,_,_,_,energyunit = PRAS.get_params(sys)
+        alpha = Float64(args["cvar_alpha"])
+        cvar_obj = PRAS.CVAR(energyunit, results["short_samples"], alpha)
+
+        cvar_value = PRAS.val(cvar_obj.cvar)
+        cvar_stderr = PRAS.stderror(cvar_obj.cvar)
+        cvar_var = cvar_obj.var
+
+        ### Normalize CVaR by total load over the full PRAS time period and report NCVAR in ppm.
+        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)"
+
+        ### Write risk metrics to CSV
+        base = replace(pras_system_path, ".pras"=>"")
+        riskfile = base * "-risk_metrics.csv"
+        dfrisk = DF.DataFrame(
+            metric = ["CVAR",       "NCVAR"],
+            alpha  = [alpha,        alpha],
+            region = ["USA",        "USA"],
+            unit   = ["MWh",        "ppm"],
+            value  = [cvar_value,   ncvar_value],
+            stderr = [cvar_stderr,  ncvar_stderr],
+            var    = [cvar_var,     ncvar_var],
+        )
+        CSV.write(riskfile, dfrisk)
+        @info("Wrote risk metrics to $(riskfile)")
+    end
+
     ## Filter out DC regions used for VSC HVDC transmission
     regions = [r for r in sys.regions.names if !(occursin("|", r))]
 
@@ -285,6 +331,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 +349,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 +514,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 +524,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"]
@@ -478,11 +537,11 @@ if abspath(PROGRAM_FILE) == @__FILE__
 
     #%% Include ReEDS2PRAS
     include(joinpath(
-        args["reedscase"], "reeds", "resource_adequacy", "reeds2pras", "src", "ReEDS2PRAS.jl"
+        args["reeds_path"], "reeds", "resource_adequacy", "reeds2pras", "src", "ReEDS2PRAS.jl"
     ))
 
     #%% Run it
     main(args)
 
     #%%
-end
+end