Bc cvar updates

cases.csv

@@ -272,6 +272,9 @@ GSw_PRM_StressThresholdMetrics,/-delimited list of metrics for identifying stres
 GSw_PRM_StressThresholdLOLE,LOLE threshold above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_nLOLE_StressMetric_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; nLOLE is number of loss of load events; StressMetric is LOLE; 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_10_LOLE_sum,
 GSw_PRM_StressThresholdEUE,/Annual EUE level [MWh] threshold above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_EUE_StressMetric_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; EUE is expected unserved energy; StressMetric EUE (only used in period selection); 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_1000_EUE_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_StressMetric_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; NEUEppm is normalized expected unserved energy in parts per million; StressMetric is EUE or NEUE (only used in period selection); 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_NEUE_sum,
+GSw_PRM_CVARAlpha,Alpha value used for CVAR/NCVAR calculations,float,0.95,
+GSw_PRM_StressThresholdNCVAR,/Annual NCVAR level [ppm] threshold used only for CVAR/NCVAR target checks; formulated as HierarchyLevel_NCVARppm_NCVAR_cvar,N/A,transgrp_10_NCVAR_cvar,
+GSw_PRM_StressThresholdCVAR,/Annual CVAR level [MWh] threshold used only for CVAR target checks; formulated as HierarchyLevel_CVAR_CVAR_cvar,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,

cases_test.csv

@@ -60,3 +60,4 @@ pras_scheduled_outage,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0
 pras_unitsize_source,,,,,,,,,,,,,,,,,,,,,,,,,,,,,r2x
 pras_vre_combine,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1
 pras_samples,,,,,,10,10,10,,,,,,,,,,,,,,,10,10,,,,,
+keep_resource_adequacy_files,1,,,,,,,,,,,,,,,,,,,,,,,,,,,,

load_atm.sh

@@ -0,0 +1,21 @@
+#!/bin/bash
+
+cd /kfs2/projects/atm/Bcakire/ReEDS || exit 1
+
+module load anaconda3
+conda activate /kfs2/projects/atm/Bcakire/conda_envs/reeds2_atm
+
+module load julia/1.12.1
+export JULIA_DEPOT_PATH=/kfs2/projects/atm/Bcakire/julia_depot
+
+echo "CONDA_PREFIX=$CONDA_PREFIX"
+echo "Python:"
+which python
+python --version
+
+echo "Julia:"
+which julia
+julia --version
+
+echo "JULIA_DEPOT_PATH=$JULIA_DEPOT_PATH"
+module load gams/51.3.0

reeds/hpc/srun_template.sh

@@ -1,9 +1,9 @@
 #!/bin/bash
-#SBATCH --account=[your HPC allocation]
+#SBATCH --account=atm
 #SBATCH --time=2-00:00:00
 #SBATCH --nodes=1
 #SBATCH --ntasks-per-node=1
-#SBATCH --mail-user=[your email address]
+#SBATCH --mail-user=Burcin.CakirErdener@nlr.gov
 #SBATCH --mail-type=BEGIN,END,FAIL
 #SBATCH --mem=246000        # RAM in MB; up to 246000 for normal or 2000000 for bigmem on kestrel
 # add >>> #SBATCH --qos=high <<< above for quicker launch at double AU cost

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 the sample-level shortfall (hourly, large)"
             arg_type = Int
             default = 0
             required = false
+        "--write_shortfall_samples_totals"
+            help = "Write per-sample total shortfall by region (small)"
+            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
@@ -204,9 +215,37 @@ function run_pras(pras_system_path::String, args::Dict)
     results = Dict{String,Any}(zip(keys(resultspec), results_tuple))
 
     #%% Print some results for the entire modeled region to show it worked
-    @info "$(PRAS.LOLE(results["short"]))"
-    @info "$( PRAS.EUE(results["short"]))"
-    @info "$(PRAS.NEUE(results["short"]))"
+    @info "$(PRAS.LOLE(results["short"])) event-h"
+    @info "$(PRAS.EUE(results["short"])) MWh"
+    @info "NEUE = $(1e6 * PRAS.EUE(results["short"]).eue.estimate / sum(sys.regions.load)) ppm"
+
+    # 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)
+        @info(cvar_obj)
+
+        total_load = sum(sys.regions.load)
+        ncvar_value  = PRAS.val(cvar_obj.cvar)      / total_load * 1e6
+        ncvar_stderr = PRAS.stderror(cvar_obj.cvar) / total_load * 1e6
+        ncvar_var    = cvar_obj.var                  / total_load * 1e6
+
+        ### 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  = [PRAS.val(cvar_obj.cvar),      ncvar_value],
+            stderr = [PRAS.stderror(cvar_obj.cvar), ncvar_stderr],
+            var    = [cvar_obj.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 +324,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 +342,39 @@ function run_pras(pras_system_path::String, args::Dict)
         @info("Wrote PRAS storage energy to $(energyfile)")
     end
 
-    ### Sample-level shortfall
+    ### Sample-level shortfall (hourly, large)
     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"]
+        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 (small, needed for CVaR)
+    if args["write_shortfall_samples_totals"] == 1
+        sf = results["short_samples"]
+        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 +505,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 +515,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 +528,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