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/4] 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/4] 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) From 65fd43aec616d3a46d47299728585df32eb0b7ac Mon Sep 17 00:00:00 2001 From: Burcin Cakir Erdener Date: Fri, 26 Jun 2026 08:45:41 -0600 Subject: [PATCH 3/4] Revert unrelated PR changes --- cases_test.csv | 1 - load_atm.sh | 21 --------------------- reeds/hpc/srun_template.sh | 4 ++-- runreeds.py | 5 ++--- 4 files changed, 4 insertions(+), 27 deletions(-) delete mode 100755 load_atm.sh diff --git a/cases_test.csv b/cases_test.csv index 10f92416..a87de967 100644 --- a/cases_test.csv +++ b/cases_test.csv @@ -63,4 +63,3 @@ pras_scheduled_outage,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0 pras_unitsize_source,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,r2x pras_vre_combine,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1 pras_samples,,,,,,10,10,10,,,,,,,,,,,,,,,,10,10,,,,, -keep_resource_adequacy_files,1,,,,,,,,,,,,,,,,,,,,,,,,,,,, \ No newline at end of file diff --git a/load_atm.sh b/load_atm.sh deleted file mode 100755 index c3bb649c..00000000 --- a/load_atm.sh +++ /dev/null @@ -1,21 +0,0 @@ -#!/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 cd4957cf..4bebd071 100644 --- a/reeds/hpc/srun_template.sh +++ b/reeds/hpc/srun_template.sh @@ -1,9 +1,9 @@ #!/bin/bash -#SBATCH --account=atm +#SBATCH --account=[your HPC allocation] #SBATCH --time=2-00:00:00 #SBATCH --nodes=1 #SBATCH --ntasks-per-node=1 -#SBATCH --mail-user=Burcin.CakirErdener@nlr.gov +#SBATCH --mail-user=[your email address] #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/runreeds.py b/runreeds.py index 4a1f05a6..a58166be 100644 --- a/runreeds.py +++ b/runreeds.py @@ -875,7 +875,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 /kfs2/projects/atm/Bcakire/conda_envs/reeds2_atm` (or `activate reeds2` on Windows)\n" + " `conda activate reeds2` (or `activate reeds2` on Windows)\n" "Do you want to continue without activating the environment?" ) confirm_env = str(input("Continue? y/[n]: ") or 'n') @@ -1308,9 +1308,8 @@ def write_batch_script( OPATH.writelines("module load conda \n") OPATH.writelines("module load gams \n") - OPATH.writelines("conda activate /kfs2/projects/atm/Bcakire/conda_envs/reeds2_atm \n") + OPATH.writelines("conda activate reeds2 \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) From e41316f78204f83033843122cc602776dd1b011e Mon Sep 17 00:00:00 2001 From: Burcin Cakir Erdener Date: Fri, 26 Jun 2026 10:17:16 -0600 Subject: [PATCH 4/4] Address CVAR/NCVAR review comments --- reeds/resource_adequacy/ra_calcs.py | 15 +++++++++ reeds/resource_adequacy/run_pras.jl | 39 ++++++++++++++--------- reeds/resource_adequacy/stress_periods.py | 39 +++++++++++++++-------- 3 files changed, 64 insertions(+), 29 deletions(-) diff --git a/reeds/resource_adequacy/ra_calcs.py b/reeds/resource_adequacy/ra_calcs.py index 214f26ab..33452980 100644 --- a/reeds/resource_adequacy/ra_calcs.py +++ b/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." diff --git a/reeds/resource_adequacy/run_pras.jl b/reeds/resource_adequacy/run_pras.jl index d3690314..80ac0387 100644 --- a/reeds/resource_adequacy/run_pras.jl +++ b/reeds/resource_adequacy/run_pras.jl @@ -71,12 +71,12 @@ function parse_commandline() default = 0 required = false "--write_shortfall_samples" - help = "Write the sample-level shortfall (hourly, large)" + 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 (small)" + help = "Write per-sample total shortfall by region over the full PRAS time period" arg_type = Int default = 0 required = false @@ -224,24 +224,31 @@ function run_pras(pras_system_path::String, args::Dict) _,_,_,_,energyunit = PRAS.get_params(sys) alpha = Float64(args["cvar_alpha"]) cvar_obj = PRAS.CVAR(energyunit, results["short_samples"], alpha) - @info(cvar_obj) + 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 = 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 + 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 = [PRAS.val(cvar_obj.cvar), ncvar_value], - stderr = [PRAS.stderror(cvar_obj.cvar), ncvar_stderr], - var = [cvar_obj.var, ncvar_var], + 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)") @@ -342,9 +349,10 @@ function run_pras(pras_system_path::String, args::Dict) @info("Wrote PRAS storage energy to $(energyfile)") end - ### Sample-level shortfall (hourly, large) + ### Per-sample hourly shortfall by region if args["write_shortfall_samples"] == 1 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] @@ -361,9 +369,10 @@ function run_pras(pras_system_path::String, args::Dict) @info("Wrote PRAS shortfall by sample to $(shortfile)") end - ### Per-sample total shortfall by region (small, needed for CVaR) + ### 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"]) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 40c7fffb..1e12b853 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -183,10 +183,10 @@ def get_and_write_neue(sw, write=True): return neue def get_cvar_alpha(sw): - val = sw.get('GSw_PRM_CVARAlpha', 0.95) - if pd.isna(val) or str(val).strip().lower() in ['', 'none', 'nan']: - return 0.95 - return float(val) + alpha = float(sw.GSw_PRM_CVARAlpha) + if not (0 <= alpha < 1): + raise ValueError(f"GSw_PRM_CVARAlpha must be in [0, 1). Got {alpha}") + return alpha def get_shortfall_totals_by_sample(case, t, iteration=0): @@ -206,8 +206,6 @@ def _sample_cvar(samples, alpha=0.95): x = pd.Series(samples).dropna().astype(float) if x.empty: return np.nan - if not (0 <= alpha < 1): - raise ValueError(f"CVaR alpha must be in [0, 1). Got alpha={alpha}") n_tail = max(1, int(np.ceil((1 - alpha) * len(x)))) return x.sort_values(ascending=False).iloc[:n_tail].mean() @@ -217,21 +215,29 @@ def get_annual_cvar_stress_metric(case, t, stress_metric='NCVAR', iteration=0, a if stress_metric not in CVAR_METRICS: raise NotImplementedError(f"get_annual_cvar_stress_metric only supports {CVAR_METRICS}. Got {stress_metric}.") + ### Read per-sample total shortfall by region from PRAS. shortfall_samples = get_shortfall_totals_by_sample(case=case, t=t, iteration=iteration).drop(columns=['USA'], errors='ignore') + + ### Read hourly PRAS load. This is used to normalize CVAR into NCVAR. dfload = reeds.io.read_h5py_file(os.path.join(case, 'handoff', 'reeds_data', f'pras_load_{t}.h5')) + ### Calculate CVAR/NCVAR at each supported hierarchy level. levels = ['country','interconnect','nercr','transreg','transgrp','st','r'] _metric = {} for hierarchy_level in levels: + ### Aggregate sample-level regional shortfall to the requested hierarchy level. rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) regions = [c for c in shortfall_samples.columns if c in rmap.index] if not regions: continue shortfall_agg = shortfall_samples[regions].rename(columns=rmap).groupby(axis=1, level=0).sum() + + ### CVAR is the average of the highest-shortfall samples in the alpha tail. cvar = shortfall_agg.apply(lambda s: _sample_cvar(s, alpha=alpha), axis=0) if stress_metric == 'NCVAR': + ### NCVAR normalizes CVAR by total load over the PRAS time period and reports ppm. load_regions = [c for c in dfload.columns if c in rmap.index] load_agg = dfload[load_regions].rename(columns=rmap).groupby(axis=1, level=0).sum().sum() cvar = cvar / load_agg.reindex(cvar.index) * 1e6 @@ -240,11 +246,11 @@ def get_annual_cvar_stress_metric(case, t, stress_metric='NCVAR', iteration=0, a return pd.concat(_metric, names=['level','metric','region']).rename(stress_metric) - def evaluate_cvar_target_check(sw, t, iteration=0, stress_metrics=None): if stress_metrics is None: stress_metrics = sw.GSw_PRM_StressThresholdMetrics.split('/') + ### Only evaluate CVAR/NCVAR here. Standard metrics are handled by _evaluate_stress_threshold_criterion(), which can add stress periods. stress_metrics = [m.upper() for m in stress_metrics if str(m).strip()] cvar_metrics = [m for m in stress_metrics if m in CVAR_METRICS] if not cvar_metrics: @@ -262,6 +268,7 @@ def evaluate_cvar_target_check(sw, t, iteration=0, stress_metrics=None): if not criterion: continue + ### CVAR/NCVAR criteria are check-only and use: HierarchyLevel_Threshold_Metric_cvar hierarchy_level, stress_level, stress_metric, metric_name = criterion.split('_') stress_metric = stress_metric.upper() metric_name = metric_name.lower() @@ -271,7 +278,13 @@ def evaluate_cvar_target_check(sw, t, iteration=0, stress_metrics=None): if metric_name != 'cvar': raise ValueError(f"Invalid CVAR criterion: {criterion}. The fourth field must be 'cvar'.") - stress_vals = pd.read_csv(os.path.join(sw.casedir,'outputs',f'{stress_metric.lower()}_{t}i{iteration}.csv',),index_col=['level', 'metric', 'region'],).squeeze(1) + ### Read the annual CVAR/NCVAR metric written earlier in main(). + stress_vals = pd.read_csv( + os.path.join(sw.casedir, 'outputs', f'{stress_metric.lower()}_{t}i{iteration}.csv'), + index_col=['level','metric','region'], + ).squeeze(1) + + ### Compare each region against the check-only threshold. this_test = stress_vals.xs((hierarchy_level, 'cvar'), level=['level','metric']) threshold = float(stress_level) @@ -815,12 +828,11 @@ def main(sw, t, iteration=0, logging=True): os.path.join(sw.casedir, 'outputs', f"{stress_metric.lower()}_{t}i{iteration}.csv") ) - # CVAR / NCVAR are additional check-only metrics. - # They do not replace or modify the standard EUE / NEUE / LOLH calculations above. - cvar_metrics = [m.upper()for m in sw.GSw_PRM_StressThresholdMetrics.split('/')if str(m).strip() and m.upper() in CVAR_METRICS] + # CVAR / NCVAR are additional check-only metrics. They do not replace or modify the standard EUE / NEUE / LOLH calculations above. + cvar_metrics = [m.upper() for m in sw.GSw_PRM_StressThresholdMetrics.split('/') if str(m).strip() and m.upper() in CVAR_METRICS] for stress_metric in cvar_metrics: - print(f"Calculating and writing annual "f"{stress_metric} for iteration {iteration}") + print(f"Calculating and writing annual {stress_metric} for iteration {iteration}") dfmetric = get_annual_cvar_stress_metric(case=sw.casedir, t=t, stress_metric=stress_metric, iteration=iteration, alpha=get_cvar_alpha(sw),) @@ -831,8 +843,7 @@ def main(sw, t, iteration=0, logging=True): dfmetric.round(2).to_csv(os.path.join(sw.casedir,'outputs', f"{stress_metric.lower()}_{t}i{iteration}.csv",)) - # CVAR / NCVAR are check-only: - # no stress periods and no PRM increment. + # CVAR / NCVAR are check-only: no stress periods and no PRM increment. evaluate_cvar_target_check(sw=sw, t=t, iteration=iteration, stress_metrics=cvar_metrics,) except Exception as err: