Add CVAR and NCVAR resource adequacy reporting

[bcakire]

Summary

This PR adds Conditional Value at Risk (CVAR) and Normalized Conditional Value at Risk (NCVAR) reporting to the ReEDS–PRAS resource adequacy workflow.

CVAR and NCVAR quantify the severity of the highest-shortfall Monte Carlo samples and complement the existing resource adequacy metrics, including NEUE, LOLH, LOLE, LOLD, Duration, and Depth.

The implementation:

• obtains sample-level shortfall results from PRAS,
• calculates CVAR and NCVAR by ReEDS hierarchy level,
• writes the results to the standard ra_metrics_{t}i{iteration}.csv output,
• treats CVAR and NCVAR as report-only metrics.

CVAR and NCVAR do not select stress periods and do not modify the existing PRM update workflow.

Technical details

run_pras.jl

Added the following command-line arguments:

• --write_shortfall_samples_totals
• --cvar_alpha

When sample-total output is requested, PRAS retains PRAS.ShortfallSamples() and writes:

PRAS_{t}i{iteration}-shortfall_totals_by_sample.h5

The H5 file contains:

• the PRAS sample index,
• total shortfall for the entire modeled system,
• total shortfall for each PRAS region.

Each value represents cumulative shortfall over all hours and resource adequacy weather years in one Monte Carlo sample.

For consistency with the existing LOLE, EUE, and NEUE diagnostics, run_pras.jl also logs whole-system:

• CVAR,
• NCVAR,
• alpha,
• VaR,
• standard error.

The CVAR value printed by Julia is the cumulative MWh value over the complete PRAS simulation horizon.

No separate CVAR/NCVAR summary CSV is written by run_pras.jl.

ra_calcs.py

Added support for requesting the new sample-total PRAS output.

When CVAR or NCVAR is included in GSw_PRM_StressThresholdMetrics, ra_calcs.py:

• enables write_shortfall_samples_totals,
• passes GSw_PRM_CVARAlpha to run_pras.jl.

The existing hourly sample output remains unchanged:

• write_shortfall_samples continues to support the existing PRAS-informed PRM update workflow,
• write_shortfall_samples_totals supplies the data needed for CVAR/NCVAR reporting.

stress_periods.py

Added the following CVAR/NCVAR reporting functions:

• get_cvar_alpha()
• get_shortfall_totals_by_sample()
• _sample_cvar()
• calc_cvar()
• calc_ncvar()

CVAR and NCVAR are calculated at the same hierarchy levels used for the other resource adequacy metrics:

• country
• interconnect
• nercr
• transreg
• transgrp
• st
• r

For each hierarchy level, shortfall is first aggregated by Monte Carlo sample. CVAR is then calculated from the highest-shortfall tail of the resulting sample distribution.

The following metrics are added to the standard resource adequacy output:

cvar_mwh_peryear
ncvar_ppm

They are written to:

outputs/ra_metrics_{t}i{iteration}.csv

using the existing long-format structure:

level,metric,region,value
transgrp,cvar_mwh_peryear,NYISO,...
transgrp,ncvar_ppm,NYISO,...
r,cvar_mwh_peryear,p127,...
r,ncvar_ppm,p127,...

CVAR and NCVAR are excluded from the stress-period selection loop. The existing prm_increment_pras() and update_prm() implementations are unchanged.

runreeds.py

Updated the resource adequacy metric compatibility check so that CVAR and NCVAR are accepted in:

GSw_PRM_StressThresholdMetrics

They are classified as report-only metrics and are skipped by the threshold parser.

Compatibility checks for the standard stress-period metrics remain unchanged. The existing requirement that PRM updates use NEUE-based failed regions is also preserved.

Implementation notes

CVAR calculation

CVAR is calculated from sample-level total shortfall, not from hourly EUE.

For each Monte Carlo sample, PRAS produces an hourly shortfall trajectory. The hourly values are summed over the complete PRAS horizon:

sample total shortfall
    = sum of hourly shortfall over all RA weather years

The sample totals are then sorted from highest to lowest. CVAR is the average shortfall of the worst (1 - alpha) fraction of samples.

For example:

pras_samples = 1000
GSw_PRM_CVARAlpha = 0.95
tail size = 50 samples

The tail-size calculation rounds the floating-point result before applying ceil. This prevents a value such as 50.00000000000004 from incorrectly selecting 51 samples.

Annualization

Each PRAS sample spans all configured resource adequacy weather years.

CVAR is first calculated from cumulative shortfall over the complete PRAS horizon. It is then divided by the number of RA weather years and reported as:

cvar_mwh_peryear

Therefore, the reported value is the annualized CVAR of cumulative shortfall over the full RA horizon. It is not the average of independently calculated weather-year CVAR values.

NCVAR calculation

NCVAR is calculated as:

NCVAR = full-horizon CVAR / full-horizon total load × 1e6

and is reported in ppm.

No additional annualization is required because both the CVAR numerator and load denominator cover the same simulation horizon.

Regional aggregation

CVAR is calculated after shortfall is aggregated by sample to the requested hierarchy level.

Regional CVAR values should not generally be added together to reproduce system-level CVAR because different regions may have different samples in their highest-shortfall tails.

Additional changes

• Updated the MultiMetricRA test case to request CVAR and NCVAR reporting.
• Updated the MultiMetricRA test configuration to use z132.
• Preserved the standard sample count in the committed test case; the 1000-sample configuration was used only for dedicated validation.
• CVAR and NCVAR use the existing consolidated ra_metrics output rather than separate reporting files.

Switches added/changed

Added: GSw_PRM_CVARAlpha

Alpha value used for CVAR and NCVAR calculations.

Type: float
Default: 0.95

The valid range is:

0 <= alpha < 1

Changed: GSw_PRM_StressThresholdMetrics

The metric list can now include:

Standard stress-period metrics:

• LOLD
• LOLE
• LOLH
• NEUE
• Depth
• Duration

Report-only metrics:

• CVAR
• NCVAR

Example:

NEUE/LOLH/LOLE/LOLD/Duration/Depth/CVAR/NCVAR

Despite the existing switch name, CVAR and NCVAR are not interpreted as threshold-based stress-period metrics.

No CVAR- or NCVAR-specific threshold switches are added in this PR.

Known incompatibilities and limitations

• CVAR and NCVAR require PRAS.ShortfallSamples(), which increases PRAS memory use when these metrics are enabled.
• Enabling the metrics produces an additional sample-total H5 file for each PRAS year and iteration.
• CVAR and NCVAR are report-only and cannot currently add stress periods.
• CVAR and NCVAR do not trigger PRM updates.
• Higher alpha values use fewer tail samples and may require more Monte Carlo samples for stable estimates.
• Whole-system CVAR printed in the Julia log is a cumulative full-horizon MWh value, while cvar_mwh_peryear in ra_metrics is annualized.

Validation, testing, and comparison report(s)

The implementation was tested using an end-to-end NY MultiMetricRA case with:

• z132,
• 15 resource adequacy weather years,
• 1000 PRAS Monte Carlo samples,
• GSw_PRM_CVARAlpha = 0.95,
• all standard MultiMetricRA metrics plus CVAR/NCVAR.

The test completed all modeled solve years and resource adequacy iterations.

Validation confirmed that:

• PRAS_{t}i{iteration}-shortfall_totals_by_sample.h5 contains one cumulative shortfall value per sample and PRAS region.
• CVAR and NCVAR are written to ra_metrics_{t}i{iteration}.csv.
• Results are produced at all supported hierarchy levels.
• CVAR is reported as annualized MWh/year.
• NCVAR is reported in ppm.
• CVAR and NCVAR do not enter the stress-period selection loop.
• The existing PRM update functions retain their existing hourly shortfall-sample workflow.
• For 1000 samples and alpha = 0.95, the corrected tail contains exactly 50 samples.
• The default workflow is unchanged when CVAR and NCVAR are not requested.

Python syntax checks passed for:

• runreeds.py
• reeds/resource_adequacy/ra_calcs.py
• reeds/resource_adequacy/stress_periods.py

The Julia environment and run_pras.jl imports were also tested successfully.

Additional diagnostic analysis included:

• CVAR by resource adequacy iteration,
• CVAR relative to iteration zero,
• final-iteration CVAR by solve year,
• final-iteration NCVAR by solve year,
• positive-shortfall sample distribution,
• shortfall exceedance curve,
• alpha sensitivity,
• regional CVAR/NCVAR comparisons.

A formal compare_cases.py report is not included because this PR does not modify ReEDS equations, model input data, or default-case results.

Checklist for author

Details to double-check

• Charge code provided to reviewers
• Included comparison reports for appropriate test cases
• Documentation updated if necessary
• Switch descriptions updated in cases.csv
• Code formatting standardized
• Reusable functions used where possible instead of copy/pasted code

General information to guide review

• Zero impact on results of default case
• No large data file(s) added or modified
• No substantive impact on runtime for the full-US reference case when CVAR/NCVAR are not enabled
• No substantive impact on folder size for the full-US reference case when CVAR/NCVAR are not enabled
• No change to the ReEDS solve flow; the runreeds.py change is limited to metric compatibility validation
• No change to code organization
• No change to package requirements (environment.yml or Project.toml)

Did you use LLM tools (chatbot or copilot) in the preparation of this PR? If so, describe how

ChatGPT was used to help draft the PR description and troubleshoot implementation, testing, and Git workflow issues. The code changes, final diffs, calculations, and test results were reviewed by the author.

Suggested reviewers

• Stress periods, resource adequacy, and ReEDS2PRAS/Julia: @patrickbrown4
• Monte Carlo and PRAS sampling: @bsergi

[patrickbrown4]

Thanks Burcin. Would you mind moving this fork to a branch on the ReEDS repo so the github runner can run the tests? (There might be a way to do it for forks but I'm not sure how.)

[patrickbrown4]

The rest of the switch testing is in runreeds.check_compatibility(), so can you move this check there, and then just use float(sw.GSw_PRM_CVARalpha) here?

[patrickbrown4]

What is the 12 here?

[patrickbrown4]

Revert so it works with the current environment

Suggested change

	dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max()
	dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum()
	dfload_agg = dfload.rename(columns=rmap).groupby(axis=1, level=0).sum()
	dflole_agg = dflole.rename(columns=rmap).T.groupby(level=0).max().T
	dfeue_agg = dfeue.rename(columns=rmap).T.groupby(level=0).sum().T
	dfload_agg = dfload.rename(columns=rmap).T.groupby(level=0).sum().T

[patrickbrown4]

For compatibility with pandas 3; see #103

Suggested change

	.groupby(axis=1, level=0)
	.sum()
	.T.groupby(level=0).T.sum()

[patrickbrown4]

Does calculating the new metrics add much runtime? If not, I would instead just skip the if block and always calculate them.

For a bit cleaner organization, I would move shortfall_samples_agg up and then just put the ra_metrics[level, 'cvar_mwh_peryear'] = cvar / numyears and ra_metrics[level, 'ncvar_ppm'] = calc_ncvar(cvar, dfload_agg) one-liners after the rest of the ra_metrics one-liners above.

[patrickbrown4]

All the metrics are always calculated even if they are not compared against the threshold or used to pick stress periods (sorry, I know this changed since the version of the branch you started from), so since CVAR and NCVAR are not yet used to select stress periods, I think this switch can remain unchanged

Suggested change

	GSw_PRM_StressThresholdMetrics,/-delimited list of resource adequacy metrics. LOLD | LOLE | LOLH | NEUE | Depth | Duration identify stress periods; CVAR and NCVAR are report-only metrics and do not add stress periods or update PRM.,N/A,NEUE,
	GSw_PRM_StressThresholdMetrics,/-delimited list of metrics for identifying stress periods; supported metrics (case-insensitive) are: LOLD | LOLE | LOLH | NEUE | Depth | Duration,N/A,NEUE,

[patrickbrown4]

(minor nitpick) personally I think it's a bit more readable to switch to lower case after an acronym, so GSw_PRM_CVARalpha instead of GSw_PRM_CVARAlpha, but not a big deal and probably debatable

[patrickbrown4]

Can this be reverted?

Suggested change

	GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,z134,z54,z48,z132,,,,,
	GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,z134,z54,z48,,,,,,

[patrickbrown4]

I think these can all be reverted

[patrickbrown4]

As noted above, this switch should only indicate metrics that are compared against thresholds and used to select stress periods. As long as the CVAR/NCVAR metrics don't add much runtime/filesize, I think we should always calculate them. We can do a followup PR to add the thresholds and stress-period selection logic.

Suggested change

	GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth/CVAR/NCVAR,,,,,
	GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth,,,,,