Skip to content

Bc cvar updates#107

Draft
bcakire wants to merge 4 commits into
ReEDS-Model:aa/multi_metricsfrom
bcakire:bc_cvar_updates
Draft

Bc cvar updates#107
bcakire wants to merge 4 commits into
ReEDS-Model:aa/multi_metricsfrom
bcakire:bc_cvar_updates

Conversation

@bcakire

@bcakire bcakire commented May 26, 2026

Copy link
Copy Markdown

Summary

Add CVaR/NCVAR resource adequacy reporting and target checks to the PRAS and stress-period workflow.

This PR extends the existing PRAS risk-metric functionality so that CVaR is no longer only calculated internally during a PRAS run. The PRAS workflow can now also write CVaR/NCVAR reporting outputs and sample-level total shortfall files that can be reused by the Python stress-period workflow.

This PR adds:

CVAR/NCVAR reporting from run_pras.jl
A compact PRAS output file with per-sample total shortfall by region
Annual CVAR/NCVAR calculation in stress_periods.py
Threshold-based CVAR/NCVAR target checks
Check-only behavior for CVAR/NCVAR so risk targets can be evaluated without adding stress periods or changing PRM

Technical details

Implementation notes

run_pras.jl updates:

Added --cvar_alpha argument to control the CVaR confidence level
Added --write_shortfall_samples_totals argument
Added CVAR/NCVAR reporting when PRAS.ShortfallSamples() is available
Writes:
PRAS_{t}i{iteration}-risk_metrics.csv
PRAS_{t}i{iteration}-shortfall_totals_by_sample.h5
risk_metrics.csv reports:
CVAR in MWh
NCVAR in ppm
alpha
estimate
standard error
VaR
shortfall_totals_by_sample.h5 stores:
sample index
USA total shortfall by sample
regional total shortfall by sample

stress_periods.py updates:

Added CVAR_METRICS = {'CVAR', 'NCVAR'}
Added support for reading:
PRAS_{t}i{iteration}-shortfall_totals_by_sample.h5
Added:
get_cvar_alpha()
get_shortfall_totals_by_sample()
_sample_cvar()
get_annual_cvar_stress_metric()
evaluate_cvar_target_check()
CVAR/NCVAR are calculated annually by hierarchy level
NCVAR is normalized by annual load and reported in ppm
CVAR/NCVAR are intentionally excluded from stress-period selection

Additional changes

Replaced memory-intensive hourly shortfall sample aggregation with direct use of PRAS_*shortfall_totals_by_sample.h5 where sample-level total shortfall is sufficient
This avoids needing to read PRAS_*shortfall_samples.h5 only to aggregate it again by sample and region
Existing EUE/NEUE/LOLE stress-period selection behavior is unchanged

Switches added/removed/changed

Added:

GSw_PRM_CVARAlpha
GSw_PRM_StressThresholdNCVAR
GSw_PRM_StressThresholdCVAR

Example usage:

GSw_PRM_StressThresholdMetrics = EUE/NCVAR
GSw_PRM_StressThresholdEUE = transgrp_1000_EUE_sum
GSw_PRM_StressThresholdNCVAR = transgrp_10_NCVAR_cvar
GSw_PRM_CVARAlpha = 0.95

For PRAS reporting, run_pras.jl can be called with:

--write_shortfall_samples_totals 1
--cvar_alpha 0.95

Validation, testing, and comparison report(s)

Validation completed on the Pacific test case.

Tests performed:

Ran run_pras.jl with:
--write_shortfall_samples_totals 1
--cvar_alpha 0.95
Verified generation of:
PRAS_2032i0-risk_metrics.csv
PRAS_2032i0-shortfall_totals_by_sample.h5
Verified that shortfall_totals_by_sample.h5 contains sample-level total shortfall by USA and region
Verified that risk_metrics.csv reports CVAR and NCVAR
Verified:
get_shortfall_totals_by_sample()
get_annual_cvar_stress_metric()
evaluate_cvar_target_check()
Confirmed NCVAR threshold failures do NOT:
add stress periods
trigger PRM updates
Confirmed standard EUE stress-period workflow remains unchanged

Example validation output:

GSw_PRM_StressThreshold = country_-1_NCVAR_cvar failed for:
region
USA 0.0
Name: NCVAR, dtype: float64
NCVAR is check-only: no stress periods and no PRM increment will be added.

Checklist for author

Details to double-check

  • Charge code provided to reviewers
  • Included comparison reports for appropriate test cases
  • Documentation updated if necessary
  • If input data added/modified:
    • Dollar year recorded and converted to 2004$ for GAMS
    • Timeseries are in Central Time
    • Units are specified
    • Preprocessing steps have been documented and committed to ReEDS_Input_Processing
    • New large data files handled with .h5 instead of .csv
    • If spatially resolved inputs are modified, the following visualizations for each file are included in the PR description (time-averaged if the inputs are time-resolved):
      • Map of absolute values before
      • Map of absolute values after
      • Map of differences: (after - before) or (after / before)
    • If entries are added/removed/changed in the EIA-NEMS unit database:
      • Changes have been committed to ReEDS_Input_Processing
      • hourlize/resource.py was rerun to regenerate the existing/prescribed VRE capacity data
  • Code formatting standardized
  • Reusable functions used where possible instead of copy/pasted code

General information to guide review

  • [X ] Zero impact on results of default case
  • [X ] No large data file(s) added/modified
  • [X ] No substantive impact on runtime for full-US reference case
  • No substantive impact on folder size for full-US reference case
  • No change to process flow (runreeds.py, reeds/core/solve/solve.py)
  • 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

Tag points of contact here if you would like additional review of the relevant parts of the model

@bcakire bcakire marked this pull request as draft May 28, 2026 06:19
@patrickbrown4 patrickbrown4 self-requested a review June 2, 2026 23:50

@bsergi bsergi left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here are some initial thoughts from a quick look. Happy to touch base when you get back to discuss any of these. Let me know when you've taken another pass here and I can review again.

One high-level comment: can you organize the PR summary text a bit more? The technical details section has a lot of information but is a bit difficult to follow. Also, I didn't quite understand your test output GSw_PRM_StressThreshold = country_-1_NCVAR_cvar failed. What does it mean to have a -1 threshold? Ultimately for the PR we'll also want a comparison of default case from cases.csv from both this branch and main to show there weren't any changes.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

revert the changes to this file before merging

required = false
"--write_shortfall_samples"
help = "Write the sample-level shortfall"
help = "Write the sample-level shortfall (hourly, large)"

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What are the large and small designations? I might suggest just "Write per-sample hourly shortfall by region" for this one and "Write per-sample total shortfall by region" for the one below.

default = 0
required = false
"--write_shortfall_samples_totals"
help = "Write per-sample total shortfall by region (small)"

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is total over the entire PRAS time period and not annual, right?

_,_,_,_,energyunit = PRAS.get_params(sys)
alpha = Float64(args["cvar_alpha"])
cvar_obj = PRAS.CVAR(energyunit, results["short_samples"], alpha)
@info(cvar_obj)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think PRAS is setup so that the above info statements provide context (e.g., @info "$(PRAS.EUE(results["short"]))" results in EUE = 851000±5000 MWh/131400h MWh) being printed). Does this do something similar?

end
## Write it
sf = results["short_samples"]
region_names = sf.regions.names

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the previous code filtered using sys.regions.names. I think this was because there are some pseudo regions for DC converter stations that don't have load and thus don't need to be included here, so we may want to continue to drop them (@patrickbrown4 might be able to weigh in on this one).

return x.sort_values(ascending=False).iloc[:n_tail].mean()


def get_annual_cvar_stress_metric(case, t, stress_metric='NCVAR', iteration=0, alpha=0.95):

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this function and the next one get pretty long, so it would be good to add comments throughout to provide some documentation on what they are doing.

x = pd.Series(samples).dropna().astype(float)
if x.empty:
return np.nan
if not (0 <= alpha < 1):

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Feels like this check might fit better in get_cvar_alpha. An even better approach might be to add a check before the ReEDS run which would avoid having runs get this far and then running into an error (best place for that would probably be here:

def check_compatibility(sw):
)


def get_annual_cvar_stress_metric(case, t, stress_metric='NCVAR', iteration=0, alpha=0.95):
stress_metric = stress_metric.upper()
if stress_metric not in CVAR_METRICS:

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The multi-metrics PR (#99) now drops the metric tag in the switch value; for example, GSw_PRM_StressThresholdNEUE now just has a value of transgrp_1_sum instead of transgrp_1_EUE_sum (see here). I think we could do the same in this PR and then you could drop the CVAR_METRICS and the check and just iterate over whichever switches are activated.

)
for criterion in sw[f'GSw_PRM_StressThreshold{metric}'].split('/'):
print(f"Evaluating GSw_PRM_StressThreshold {metric} with criterion: {criterion}")
stress_criteria = _evaluate_stress_threshold_criterion(stress_criteria, criterion, sw, t, iteration, dfenergy_r, stressperiods_this_iteration,

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line looks like just a formatting change--I'd suggest keeping as the original had it.

return pd.concat(_metric, names=['level','metric','region']).rename(stress_metric)


def evaluate_cvar_target_check(sw, t, iteration=0, stress_metrics=None):

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Much of this function, including the core check of failed = this_test.loc[this_test > threshold], seems to replicate behavior in _evaluate_stress_threshold_criterion. Do you think we could merge the CVAR treatment into that existing function and just have some control statements to turn off adding new stress periods for now when using CVAR metrics? It's possible having a second function makes the most sense here, but I worry a little bit about maintaining two similar functions.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants