From 29b6be0fceff61b5111ad255adcd747b82488b09 Mon Sep 17 00:00:00 2001 From: Burcin Cakir Erdener Date: Wed, 20 May 2026 05:38:13 -0600 Subject: [PATCH 1/2] add CVaR reporting to run_pras.jl --- reeds/resource_adequacy/run_pras.jl | 89 ++++++++++++++++++++++------- 1 file changed, 67 insertions(+), 22 deletions(-) diff --git a/reeds/resource_adequacy/run_pras.jl b/reeds/resource_adequacy/run_pras.jl index a100ac74..7b5ebbdd 100644 --- a/reeds/resource_adequacy/run_pras.jl +++ b/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 @@ -209,9 +220,32 @@ function run_pras(pras_system_path::String, args::Dict) @info "NEUE = $(1e6 * PRAS.EUE(results["short"]).eue.estimate / sum(sys.regions.load)) ppm" # CVAR results - _,_,_,_,energyunit = PRAS.get_params(sys) - alpha = 0.95 - @info (PRAS.CVAR(energyunit, results["short_samples"], alpha)) + 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))] @@ -290,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() @@ -307,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"]) @@ -462,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, @@ -471,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"] @@ -483,11 +528,11 @@ if abspath(PROGRAM_FILE) == @__FILE__ #%% Include ReEDS2PRAS include(joinpath( - args["reedscase"], "reeds", "resource_adequacy", "reeds2pras", "src", "ReEDS2PRAS.jl" + args["reedscase"], "reeds2pras", "src", "ReEDS2PRAS.jl" )) #%% Run it main(args) #%% -end +end \ No newline at end of file From 33d2a02d53aa71c9311a5b2079f9c8860825bc8a Mon Sep 17 00:00:00 2001 From: Burcin Cakir Erdener Date: Mon, 25 May 2026 16:38:48 -0600 Subject: [PATCH 2/2] fix fix --- cases_test.csv | 1 + load_atm.sh | 21 +++++++++++++++++++++ reeds/hpc/srun_template.sh | 4 ++-- reeds/resource_adequacy/run_pras.jl | 2 +- runreeds.py | 5 +++-- 5 files changed, 28 insertions(+), 5 deletions(-) create mode 100755 load_atm.sh diff --git a/cases_test.csv b/cases_test.csv index e97c4207..90a6773a 100644 --- a/cases_test.csv +++ b/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,,,,,,,,,,,,,,,,,,,,,,,,,,,, diff --git a/load_atm.sh b/load_atm.sh new file mode 100755 index 00000000..c3bb649c --- /dev/null +++ b/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 diff --git a/reeds/hpc/srun_template.sh b/reeds/hpc/srun_template.sh index 4bebd071..cd4957cf 100644 --- a/reeds/hpc/srun_template.sh +++ b/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 \ No newline at end of file diff --git a/reeds/resource_adequacy/run_pras.jl b/reeds/resource_adequacy/run_pras.jl index 7b5ebbdd..ac694d55 100644 --- a/reeds/resource_adequacy/run_pras.jl +++ b/reeds/resource_adequacy/run_pras.jl @@ -528,7 +528,7 @@ if abspath(PROGRAM_FILE) == @__FILE__ #%% Include ReEDS2PRAS include(joinpath( - args["reedscase"], "reeds2pras", "src", "ReEDS2PRAS.jl" + args["reeds_path"], "reeds", "resource_adequacy", "reeds2pras", "src", "ReEDS2PRAS.jl" )) #%% Run it diff --git a/runreeds.py b/runreeds.py index 397886af..d67721bc 100644 --- a/runreeds.py +++ b/runreeds.py @@ -859,7 +859,7 @@ def setupEnvironment( "To build the environment for the first time, run:\n" " `conda env create -f environment.yml`\n" "To activate the created environment, run:\n" - " `conda activate reeds2` (or `activate reeds2` on Windows)\n" + " `conda activate /kfs2/projects/atm/Bcakire/conda_envs/reeds2_atm` (or `activate reeds2` on Windows)\n" "Do you want to continue without activating the environment?" ) confirm_env = str(input("Continue? y/[n]: ") or 'n') @@ -1292,8 +1292,9 @@ def write_batch_script( OPATH.writelines("module load conda \n") OPATH.writelines("module load gams \n") - OPATH.writelines("conda activate reeds2 \n") + OPATH.writelines("conda activate /kfs2/projects/atm/Bcakire/conda_envs/reeds2_atm \n") OPATH.writelines('export R_LIBS_USER="$HOME/rlib" \n\n\n') + OPATH.writelines("export JULIA_DEPOT_PATH=/kfs2/projects/atm/Bcakire/julia_depot \n") #%% Write the input_processing script calls big_comment('Input processing', OPATH)