From 475a94af0c565dc4c0c44ab59b7cee43ea061893 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Tue, 23 Jun 2026 16:34:38 -0600 Subject: [PATCH 01/24] stress_periods.py: combine ra metrics into one csv; refactor metric calculations and differentiate LOLE and LOLD --- reeds/resource_adequacy/stress_periods.py | 332 +++++++++++----------- 1 file changed, 171 insertions(+), 161 deletions(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 9e50b5cd..77812b3a 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -6,6 +6,7 @@ from glob import glob import re import matplotlib.pyplot as plt +from pathlib import Path from typing import Literal import reeds @@ -19,11 +20,11 @@ #%%### Functions -def get_pras_stress_metric(case, t, iteration=0, stress_metric:Literal['EUE','LOLE']='EUE'): +def get_pras_shortfall(case, t, iteration=0): """ + Returns: dict of timeseries-indexed dataframes with two keys: 'EUE' and 'LOLE' """ ### Get PRAS outputs - dfpras = reeds.io.read_pras_results( os.path.join(case, 'handoff', 'PRAS', f"PRAS_{t}i{iteration}.h5") ) @@ -32,16 +33,22 @@ def get_pras_stress_metric(case, t, iteration=0, stress_metric:Literal['EUE','LO dfpras.index = reeds.timeseries.get_timeindex(sw['resource_adequacy_years']) ### Keep the metric columns by zone - metric_tail = '_' + stress_metric.upper() - dfmetric = dfpras[[ - c for c in dfpras - if (c.endswith(metric_tail) and not c.startswith('USA')) - ]].copy() - ## Drop the tailing metric tail - dfmetric = dfmetric.rename( - columns=dict(zip(dfmetric.columns, [c[:-len(metric_tail)] for c in dfmetric]))) + dictout = {} + for metric in ['EUE', 'LOLE']: + metric_tail = '_' + metric.upper() + dfmetric = dfpras[[ + c for c in dfpras + if (c.endswith(metric_tail) and not c.startswith('USA')) + ]].copy() + ## Drop the tailing metric tail + dictout[metric] = dfmetric.rename( + columns=dict(zip( + dfmetric.columns, + [c[:-len(metric_tail)] for c in dfmetric] + )) + ) - return dfmetric + return dictout def get_stress_metric_periods( @@ -76,14 +83,17 @@ def get_stress_metric_periods( # Use EUE metric for both EUE and NEUE calculations, since the load division to get NEUE is # peformed afterwards based on agg_period. use_metric_for_pras = { - 'EUE': 'EUE', 'NEUE': 'EUE', 'LOLH': 'LOLE', 'LOLE': 'LOLE', - 'OutageDuration': 'EUE', 'OutageMagnitude': 'EUE', 'NormalizedOutageMagnitude': 'EUE', + 'EUE': 'EUE', + 'NEUE': 'EUE', + 'LOLH': 'LOLE', + 'LOLE': 'LOLE', + 'OutageDuration': 'EUE', + 'OutageMagnitude': 'EUE', + 'NormalizedOutageMagnitude': 'EUE', } - dfmetric = get_pras_stress_metric( - case=case, - t=t, - iteration=iteration, - stress_metric=use_metric_for_pras.get(stress_metric) + dfmetric = ( + get_pras_shortfall(case=case, t=t, iteration=iteration) + [use_metric_for_pras['stress_metric']] ) ## Aggregate to hierarchy_level dfmetric = ( @@ -184,85 +194,157 @@ def get_and_write_neue(sw, write=True): return neue -def get_annual_stress_metric(case, t, stress_metric, iteration=0): - """ - """ - ### Get values from PRAS - # Use EUE for outages duration calculation - use_metric_for_pras = {'EUE':'EUE', - 'NEUE':'EUE', - 'LOLH':'LOLE', - 'LOLE':'LOLE', - 'OutageDuration':'EUE', - 'OutageMagnitude':'EUE', - 'NormalizedOutageMagnitude':'EUE', - } - dfmetric = get_pras_stress_metric( - case=case, - t=t, - iteration=iteration, - stress_metric=use_metric_for_pras[stress_metric], +def get_events(ds:pd.Series, threshold:float=0) -> pd.DataFrame: + """Return a dataframe of events with max and duration""" + starts = ( + ## Convert values > threshold to 1 + (ds > threshold).astype(int) + ## +1 if changes from 0->1, -1 if changes from 1->0 + .diff() + ## If the first value is > threshold, count it as a start + .fillna((ds > threshold).astype(int)) + ## Only keep the beginnings + > 0 + ) + starts = starts.loc[starts > 0].index + ## Same idea for ends but reverse + ends = ( + (ds > threshold).astype(int) + .diff(-1) + .fillna((ds > threshold).astype(int)) + > 0 ) + ends = ends.loc[ends > 0].index + assert len(starts) == len(ends), "Error in event start/end calculation" + ## Get some metrics for each event + dfout = [] + for start, end in zip(starts, ends): + event = ds.loc[start:end] + dfout.append({ + 'start': start, + 'end': end, + 'timesteps': len(event), + 'max': event.max(), + 'sum': event.sum(), + }) + return pd.DataFrame(dfout) + + +def calc_lold(dflole, rmap, threshold=0): + """Count a day as an event-day if at least one hour has LOLE > threshold""" + ## If multiple zones in one level and hour have LOLE, count that as one event, + ## so take the max LOLE across the zones + dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max() + ## Similarly, take the max for each day + ## (That's not quite right if the events are independent) + daily_max = dflole_agg.groupby( + [dflole_agg.index.year, dflole_agg.index.month, dflole_agg.index.day] + ).max() + ## Sum the probability across days + lold = daily_max.sum() + return lold + + +def calc_lole(dflole, rmap, threshold=0): + """Number of events, where an event is >threshold LOLE in contiguous hours""" + ## If multiple zones in one level and hour have LOLE, count that as one event, + ## so take the max LOLE across the zones + dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max() + ## Get the loss-of-load events, keeping the max hourly probability for each + ## (not really right probabilistically) + lole = pd.Series({ + r: get_events(dflole_agg[r], threshold)['max'].sum() for r in dflole_agg + }) + return lole - ### Get load (for calculating NEUE) - if stress_metric.upper() == 'NEUE' or stress_metric == 'NormalizedOutageMagnitude': - dfload = reeds.io.read_h5py_file( - os.path.join( - case,'handoff','reeds_data',f'pras_load_{t}.h5') - ) - dfload.index = dfmetric.index - levels = ['country','interconnect','nercr','transreg','transgrp','st','r'] - _metric = {} - for hierarchy_level in levels: - ### Get the region aggregator - rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) +def calc_max_duration(dfeue, rmap, threshold=0): + """Max event duration, where an event is >threshold EUE [MW] in contiguous hours""" + dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum() + max_duration = pd.Series({ + r: get_events(dfeue_agg[r], threshold)['timesteps'].max() for r in dfeue_agg + }) + return max_duration - if stress_metric == 'LOLE': - # Count a day as an event-day if at least one hour has LOLE > 0 - dfmetric_agg = dfmetric.rename(columns=rmap).groupby(axis=1, level=0).sum() - daily_max = dfmetric_agg.groupby( - [dfmetric_agg.index.year, dfmetric_agg.index.month, dfmetric_agg.index.day] - ).max() - event_days = (daily_max > 0).sum() - _metric[hierarchy_level, 'sum'] = event_days - _metric[hierarchy_level, 'max'] = event_days - continue - if stress_metric in ['OutageDuration','OutageMagnitude','NormalizedOutageMagnitude']: - _metric[hierarchy_level, 'sum'] = get_sum_duration_outage(dfmetric) +def calc_neue(dfeue, dfload, rmap): + """NEUE (sum of EUE / sum of load) in units of ppm""" + dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum() + dfload_agg = dfload.rename(columns=rmap).groupby(axis=1, level=0).sum() + neue = dfeue_agg.sum() / dfload_agg.sum() * 1e6 + return neue - ## TODO: for 'sum', OutageMagnitude returns the same value as `sum` for OutageDuration - if stress_metric == 'OutageMagnitude': - _metric[hierarchy_level, 'max'] = get_max_magnitude_outage(dfmetric) - if stress_metric == 'NormalizedOutageMagnitude': - _metric[hierarchy_level, 'max'] = get_max_magnitude_outage(dfmetric) / dfload.max() - if stress_metric == 'OutageDuration': - _metric[hierarchy_level, 'max'] = get_max_duration_outage(dfmetric) - continue +def calc_peak_eue(dfeue, dfload, rmap, norm:Literal['peak','hourly','absolute']='peak'): + """ + Get the peak hourly outage magnitude - ### Get stress metric summed over year - _metric[hierarchy_level,'sum'] = dfmetric.rename(columns=rmap).groupby(axis=1, level=0).sum().sum() + Args: + norm: How to normalize the hourly EUE [MW] + - 'peak': Divide hourly EUE by peak load -> returns fraction + - 'hourly': Divide hourly EUE by hourly load -> returns fraction + - 'absolute': Do not normalize -> returns MW + """ + dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum() + dfload_agg = dfload.rename(columns=rmap).groupby(axis=1, level=0).sum() + match norm: + case 'peak': + peak_eue = dfeue_agg.max() / dfload_agg.max() + case 'hourly': + peak_eue = (dfeue_agg / dfload_agg).max() + case 'absolute': + peak_eue = dfeue.max() + return peak_eue + + +def calc_ra_metrics( + case:str|Path, + t:int, + iteration:int=0, + levels=['country', 'interconnect', 'nercr', 'transreg', 'transgrp', 'st', 'r'], +): + """ + Calculate all resource adequacy metrics for the specified year/iteration + and regional aggregation levels. + """ + ### Validate inputs + hierarchy = reeds.io.get_hierarchy(case).reset_index() + wrong = [i for i in levels if i not in hierarchy] + if len(wrong): + raise ValueError(f'Invalid levels: {wrong}') - ### Get max stress metric hour - _metric[hierarchy_level,'max'] = dfmetric.rename(columns=rmap).groupby(axis=1, level=0).sum().max() + ### Get values from PRAS + pras_shortfall = get_pras_shortfall(case, t, iteration) + dfeue = pras_shortfall['EUE'] + dflole = pras_shortfall['LOLE'] - if stress_metric.upper() == 'NEUE': - _metric[hierarchy_level,'sum'] = ( - dfmetric.rename(columns=rmap).groupby(axis=1, level=0).sum().sum() - / dfload.rename(columns=rmap).groupby(axis=1, level=0).sum().sum() - ) * 1e6 - ### Get max NEUE hour - _metric[hierarchy_level,'max'] = ( - dfmetric.rename(columns=rmap).groupby(axis=1, level=0).sum() - / dfload.rename(columns=rmap).groupby(axis=1, level=0).sum() - ).max() * 1e6 + ### Get load and number of years (for normalization) + dfload = reeds.io.read_h5py_file( + Path(case,'handoff','reeds_data',f'pras_load_{t}.h5') + ) + dfload.index = dfeue.index + sw = reeds.io.get_switches(case) + numyears = len(sw.resource_adequacy_years_list) + + ### Loop over aggregation levels and calculate all metrics + ra_metrics = {} + for hierarchy_level in levels: + print(f'Calculating RA metrics at {hierarchy_level} level') + ## Get the region aggregator + rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) + ## Calculate the full-timeseries metrics for each region + ra_metrics[hierarchy_level, 'lold_peryear'] = calc_lold(dflole, rmap) / numyears + ra_metrics[hierarchy_level, 'lole_peryear'] = calc_lole(dflole, rmap) / numyears + ra_metrics[hierarchy_level, 'max_duration'] = calc_max_duration(dfeue, rmap) + ra_metrics[hierarchy_level, 'neue_ppm'] = calc_neue(dfeue, dfload, rmap) + ra_metrics[hierarchy_level, 'euemax_peakloadfrac'] = calc_peak_eue(dfeue, dfload, rmap, 'peak') + ra_metrics[hierarchy_level, 'euemax_hourlyloadfrac'] = calc_peak_eue(dfeue, dfload, rmap, 'hourly') + ra_metrics[hierarchy_level, 'euemax_mw'] = calc_peak_eue(dfeue, dfload, rmap, 'absolute') ### Combine it - metric = pd.concat(_metric, names=['level','metric','region']).rename(f'{stress_metric}') + dfout = pd.concat(ra_metrics, names=['level','metric','region']).rename('value') - return metric + return dfout def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods, stress_metric): @@ -755,91 +837,19 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): return prm_next_iteration -def get_max_duration_outage(dfmetric, eue_threshold = 0): - """Return the longest consecutive outage event in hours per region. - - Check regions hourly EUE timeseries for runs of contiguous hours where - EUE exceeds ``eue_threshold``, and returns the length of the longest such run. - A single isolated hour counts as 1. - - Args: - dfmetric (pd.DataFrame): Hourly EUE values (MW) indexed by timestamp, one - column per region. Typically obtained from :func:`get_pras_stress_metric`. - eue_threshold (float, optional): Minimum EUE (MW) to qualify as an outage - hour. Hours at or below this value are treated as no-outage. Defaults to - 0 (any positive shortfall counts). - - Returns: - pd.Series: Length in hours of the longest consecutive outage event, indexed - by region. Regions with no outage hours return 0. - """ - def _max_run(series): - is_outage = series > eue_threshold - groups = (is_outage != is_outage.shift()).cumsum()[is_outage] - return series[is_outage].groupby(groups).size().max() if is_outage.any() else 0 - return dfmetric.apply(_max_run) - - -def get_sum_duration_outage(dfmetric, eue_threshold = 0): - """Return the total number of outage hours per region. - - Counts every hour in the timeseries where EUE exceeds ``eue_threshold`` - - Args: - dfmetric (pd.DataFrame): Hourly EUE values (MW) indexed by timestamp, one - column per region. Typically obtained from :func:`get_pras_stress_metric`. - eue_threshold (float, optional): Minimum EUE (MW) to qualify as an outage - hour. Hours at or below this value are excluded. Defaults to 0. - - Returns: - pd.Series: Total outage hours over the full timeseries, indexed by region. - Regions with no outage hours return 0. - """ - return (dfmetric > eue_threshold).sum() - - -def get_max_magnitude_outage(dfmetric): - """Return the peak single-hour EUE per region (MW). - - Finds the highest magnitude of shortfall for each region. Hours at or - below ``eue_threshold`` are masked before taking the maximum, so they cannot - inflate the result. Regions with no outage hours return 0. - - Args: - dfmetric (pd.DataFrame): Hourly EUE values (MW) indexed by timestamp, one - column per region. Typically obtained from :func:`get_pras_stress_metric`. - eue_threshold (float, optional): Minimum EUE (MW) to qualify as an outage - hour. Hours at or below this value are excluded. Defaults to 0. - - Returns: - pd.Series: Peak hourly EUE (MW) over the full timeseries, indexed by region. - Regions with no outage hours return 0. - """ - return dfmetric.max() - - #%%### Procedure def main(sw, t, iteration=0, logging=True): """ """ #%% Write consolidated stress metrics so far try: - ## TODO: Check if get_and_write_neue() is still needed + ## TODO: Remove get_and_write_neue() _neue_simple = get_and_write_neue(sw, write=True) - ## TODO: check if need to refactor or remove - - for stress_metric in ['EUE', 'NEUE', 'LOLH', 'LOLE', 'OutageDuration', 'OutageMagnitude', 'NormalizedOutageMagnitude']: - print(f"Calculating and writing annual {stress_metric} for iteration {iteration}") - dfmetric = get_annual_stress_metric(sw.casedir, t, stress_metric, iteration=iteration) - # Metric thresholds are defined on a per-year basis, but PRAS reports the total over all resource adequacy years, - # For all metrics except NEUE, divide by the number of resource adequacy years to get the average per year - if stress_metric != 'NEUE': - dfmetric /= len(sw['resource_adequacy_years']) - - dfmetric.round(2).to_csv( - os.path.join(sw.casedir, 'outputs', f"{stress_metric.lower()}_{t}i{iteration}.csv") - ) + ra_metrics = calc_ra_metrics(case=sw.casedir, t=t, iteration=iteration) + ra_metrics.round(3).to_csv( + os.path.join(sw.casedir, 'outputs', f'ra_metrics_{t}i{iteration}.csv') + ) except Exception as err: if int(sw['pras']) == 2: From d7f5dae9cf6c2ff51c2b86f2c3eac5fad1b8ef5a Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Tue, 23 Jun 2026 16:39:00 -0600 Subject: [PATCH 02/24] stress_periods.py: whitespace --- reeds/resource_adequacy/stress_periods.py | 75 +++++++++++++---------- 1 file changed, 43 insertions(+), 32 deletions(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 77812b3a..065b9f6c 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -81,7 +81,7 @@ def get_stress_metric_periods( ### Get stress metric from PRAS # Use EUE metric for both EUE and NEUE calculations, since the load division to get NEUE is - # peformed afterwards based on agg_period. + # peformed afterwards based on agg_period. use_metric_for_pras = { 'EUE': 'EUE', 'NEUE': 'EUE', @@ -108,7 +108,7 @@ def get_stress_metric_periods( .groupby([dfmetric.index.year, dfmetric.index.month, dfmetric.index.day]) .agg(period_agg_method) .rename_axis(['y','m','d']) - ) + ) ### Sort and drop zeros and duplicates dfmetric_top = ( @@ -423,22 +423,22 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods, stress_met return shoulder_periods -def _evaluate_stress_threshold_criterion( - stress_criteria, - criterion, - sw, - t, - iteration, - dfenergy_r, - stressperiods_this_iteration, - stress_metric, -): +def _evaluate_stress_threshold_criterion( + stress_criteria, + criterion, + sw, + t, + iteration, + dfenergy_r, + stressperiods_this_iteration, + stress_metric, +): _stress_sorted_periods = stress_criteria['stress_sorted_periods'] _failed = stress_criteria['failed'] _high_stress_periods = stress_criteria['high_stress_periods'] _shoulder_periods = stress_criteria['shoulder_periods'] - - ## NEUE Example: criterion = 'transgrp_1_sum' + + ## NEUE Example: criterion = 'transgrp_1_sum' (hierarchy_level, stress_value, period_agg_method) = criterion.split('_') ### Get stored stress metric @@ -521,23 +521,24 @@ def _evaluate_stress_threshold_criterion( _shoulder_periods = { **_shoulder_periods, **get_shoulder_periods( - sw, - criterion, - dfenergy_r, - _high_stress_periods, - stress_metric=stress_metric_for_shoulder_periods) + sw, + criterion, + dfenergy_r, + _high_stress_periods, + stress_metric=stress_metric_for_shoulder_periods, + ) } stress_criteria['stress_sorted_periods'] = _stress_sorted_periods stress_criteria['failed'] = _failed stress_criteria['high_stress_periods'] = _high_stress_periods stress_criteria['shoulder_periods'] = _shoulder_periods - + else: print(f"GSw_PRM_StressThreshold = {criterion} passed") return stress_criteria - + def get_stress_metrics_sorted_periods(sw, t, iteration): ### Get storage state of charge (SOC) to use in selection of "shoulder" stress periods @@ -560,20 +561,31 @@ def get_stress_metrics_sorted_periods(sw, t, iteration): ) ### Check all stress criteria; for regions that fail, add new stress periods - stress_criteria = {'stress_sorted_periods': {}, 'failed': {}, 'high_stress_periods': {}, 'shoulder_periods': {}} + stress_criteria = { + 'stress_sorted_periods': {}, + 'failed': {}, + 'high_stress_periods': {}, + 'shoulder_periods': {}, + } # Validation check: Display any GSw_PRM_StressThreshold{metric} # that is not specified in GSw_PRM_StressThresholdMetrics stress_metric_switches = sw.GSw_PRM_StressThresholdMetrics.split('/') - + # stress periods column names for writing outputs stress_metric_units = { - 'EUE':'MWh/year', 'NEUE':'ppm', 'LOLH':'event-h/year', 'LOLE':'event-day/year', - 'OutageDuration':'h', 'OutageMagnitude':'MW', 'NormalizedOutageMagnitude':'p.u. of load', + 'EUE':'MWh/year', + 'NEUE':'ppm', + 'LOLH':'event-h/year', + 'LOLE':'event-day/year', + 'OutageDuration':'h', + 'OutageMagnitude':'MW', + 'NormalizedOutageMagnitude':'p.u. of load', + } + stress_metrics_col_names = { + m: f'{m}_{stress_metric_units[m]}' for m in stress_metric_switches } - stress_metrics_col_names = {m:f'{m}_{stress_metric_units[m]}' for m in - stress_metric_switches} - + for stress_metric in stress_metric_switches: for criterion in sw[f'GSw_PRM_StressThreshold{stress_metric}'].split('/'): print(f"Evaluating GSw_PRM_StressThreshold {stress_metric.upper()} with criterion: {criterion}") @@ -587,7 +599,7 @@ def get_stress_metrics_sorted_periods(sw, t, iteration): stressperiods_this_iteration, stress_metric ) - + stress_sorted_periods = stress_criteria['stress_sorted_periods'] failed = stress_criteria['failed'] high_stress_periods = stress_criteria['high_stress_periods'] @@ -597,7 +609,7 @@ def get_stress_metrics_sorted_periods(sw, t, iteration): ### Get lists of stress periods: new (added this iteration) and all ##TODO: What if same high stress period is added for multiple criteria? - # Currently the duplicates are dropped. + # Currently the duplicates are dropped. if len(failed): new_stress_periods = pd.concat( {**high_stress_periods, **shoulder_periods}, names=['criterion','periodtype'], @@ -781,7 +793,6 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): failed (dict): Dictionary of regions with unserved energy at the hierarchy_level and their criterion evaluations combined_periods_write (pd.DataFrame): Data frame of combined stress periods - Returns: pd.DataFrame: Table of prm levels for the next PRAS iteration @@ -807,7 +818,7 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): pd.concat(_failed_regions) .sort_values(by=['stress_value']) .drop_duplicates(subset='r', keep='first') - ) + ) ## Fixed-increment update if int(sw.GSw_PRM_UpdateMethod) == 1: From 60ed855f24934fe08bb54051996fcd1b6e067e9a Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Tue, 23 Jun 2026 17:14:09 -0600 Subject: [PATCH 03/24] remove get_and_write_neue() and simplified national neue.csv/eue.csv outputs --- postprocessing/bokehpivot/reeds2.py | 13 ---- .../reeds2/standard_report_expanded.py | 1 - postprocessing/compare_cases.py | 1 - reeds/resource_adequacy/stress_periods.py | 59 ++----------------- 4 files changed, 4 insertions(+), 70 deletions(-) diff --git a/postprocessing/bokehpivot/reeds2.py b/postprocessing/bokehpivot/reeds2.py index dcb86cfc..d2a6269f 100644 --- a/postprocessing/bokehpivot/reeds2.py +++ b/postprocessing/bokehpivot/reeds2.py @@ -2761,19 +2761,6 @@ def pre_spur(dfs, **kw): } ), - ('Resource adequacy', - {'file':'neue.csv', - 'columns': ['year', 'iteration', 'neue'], - 'index': ['year', 'iteration'], - 'preprocess': [ - {'func': pre_neue, 'args': {}}, - ], - 'presets': collections.OrderedDict(( - ('Normalized expected unserved energy [%]', {'x':'year', 'y':'neue', 'chart_type':'Bar', 'explode':'scenario', 'bar_width':'1.75'}), - )), - } - ), - ('Average VRE Curtailment', {'file':'curt_rate', 'columns': ['year', 'Curt Rate'], diff --git a/postprocessing/bokehpivot/reports/templates/reeds2/standard_report_expanded.py b/postprocessing/bokehpivot/reports/templates/reeds2/standard_report_expanded.py index 601e4d86..b19e06ae 100644 --- a/postprocessing/bokehpivot/reports/templates/reeds2/standard_report_expanded.py +++ b/postprocessing/bokehpivot/reports/templates/reeds2/standard_report_expanded.py @@ -54,7 +54,6 @@ {'name': 'Total mortality through 2050 (lives)', 'result': 'Health Damages from Emissions', 'preset': 'Total mortality through 2050'}, {'name': 'System cost + health damages: ACS ($/MWh)', 'result': 'Total Social Costs', 'preset': 'System cost + health damages: ACS CR'}, {'name': 'System cost + health damages: H6C ($/MWh)', 'result': 'Total Social Costs', 'preset': 'System cost + health damages: H6C CR'}, - {'name': 'NEUE (ppm)', 'result': 'Resource adequacy', 'preset': 'Normalized expected unserved energy [%]'}, {'name': 'Runtime (hours)', 'result': 'Runtime', 'preset': 'Stacked Bars'}, {'name': 'Runtime by year (hours)', 'result': 'Runtime', 'preset': 'Stacked Bars by Year'}, ] diff --git a/postprocessing/compare_cases.py b/postprocessing/compare_cases.py index 56bada2c..56c91d8d 100644 --- a/postprocessing/compare_cases.py +++ b/postprocessing/compare_cases.py @@ -193,7 +193,6 @@ 'Bulk System Electricity Pric', 'National Average Electricity', 'Present Value of System Cost', - 'NEUE (ppm)', 'Runtime (hours)', 'Runtime by year (hours)', ] diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 065b9f6c..563f59e3 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -1,10 +1,7 @@ #%%### General imports import os -import traceback import pandas as pd import numpy as np -from glob import glob -import re import matplotlib.pyplot as plt from pathlib import Path from typing import Literal @@ -156,44 +153,6 @@ def plot_stress_diagnostics(sw, t, iteration, high_stress_periods): print(err) -def get_and_write_neue(sw, write=True): - """ - Write dropped load across all completed years to outputs - so it can be plotted alongside other ReEDS outputs. - - Notes - ----- - * The denominator of NEUE is exogenous electricity demand; it does not include - endogenous load from losses or H2 production or exogenous H2 demand. - """ - infiles = [ - i for i in sorted(glob( - os.path.join(sw['casedir'], 'handoff', 'PRAS', 'PRAS_*.h5'))) - if re.match(r"PRAS_[0-9]+i[0-9]+.h5", os.path.basename(i)) - ] - eue = {} - for infile in infiles: - year_iteration = os.path.basename(infile)[len('PRAS_'):-len('.h5')].split('i') - year = int(year_iteration[0]) - iteration = int(year_iteration[1]) - eue[year,iteration] = reeds.io.read_pras_results(infile)['USA_EUE'].sum() - eue = pd.Series(eue).rename('MWh') - eue.index = eue.index.rename(['year','iteration']) - - load = reeds.io.read_file(os.path.join(sw['casedir'],'inputs_case','load.h5')) - loadyear = load.sum(axis=1).groupby('year').sum() - - neue = ( - (eue / loadyear * 1e6).rename('NEUE [ppm]') - .rename_axis(['t','iteration']).sort_index() - ) - - if write: - neue.to_csv(os.path.join(sw['casedir'],'outputs','neue.csv')) - eue.to_csv(os.path.join(sw['casedir'],'outputs','eue.csv')) - return neue - - def get_events(ds:pd.Series, threshold:float=0) -> pd.DataFrame: """Return a dataframe of events with max and duration""" starts = ( @@ -853,20 +812,10 @@ def main(sw, t, iteration=0, logging=True): """ """ #%% Write consolidated stress metrics so far - try: - ## TODO: Remove get_and_write_neue() - _neue_simple = get_and_write_neue(sw, write=True) - - ra_metrics = calc_ra_metrics(case=sw.casedir, t=t, iteration=iteration) - ra_metrics.round(3).to_csv( - os.path.join(sw.casedir, 'outputs', f'ra_metrics_{t}i{iteration}.csv') - ) - - except Exception as err: - if int(sw['pras']) == 2: - print(traceback.format_exc()) - if int(sw.GSw_PRM_StressIterateMax): - raise Exception(err) + ra_metrics = calc_ra_metrics(case=sw.casedir, t=t, iteration=iteration) + ra_metrics.round(3).to_csv( + os.path.join(sw.casedir, 'outputs', f'ra_metrics_{t}i{iteration}.csv') + ) #%% Stop here if not iterating or if before ReEDS can build new capacity if (not int(sw.GSw_PRM_StressIterateMax)) or (t < int(sw['GSw_StartMarkets'])): From c3c4400e0a50589915c53e2db8b79082290c8cc6 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Tue, 23 Jun 2026 17:22:19 -0600 Subject: [PATCH 04/24] calc_ra_metrics(): aggregate before the ra calculation functions --- reeds/resource_adequacy/stress_periods.py | 46 ++++++++++------------- 1 file changed, 20 insertions(+), 26 deletions(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 563f59e3..f1ea5b65 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -189,12 +189,9 @@ def get_events(ds:pd.Series, threshold:float=0) -> pd.DataFrame: return pd.DataFrame(dfout) -def calc_lold(dflole, rmap, threshold=0): +def calc_lold(dflole_agg, threshold=0): """Count a day as an event-day if at least one hour has LOLE > threshold""" - ## If multiple zones in one level and hour have LOLE, count that as one event, - ## so take the max LOLE across the zones - dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max() - ## Similarly, take the max for each day + ## Take the max for each day ## (That's not quite right if the events are independent) daily_max = dflole_agg.groupby( [dflole_agg.index.year, dflole_agg.index.month, dflole_agg.index.day] @@ -204,11 +201,8 @@ def calc_lold(dflole, rmap, threshold=0): return lold -def calc_lole(dflole, rmap, threshold=0): +def calc_lole(dflole_agg, threshold=0): """Number of events, where an event is >threshold LOLE in contiguous hours""" - ## If multiple zones in one level and hour have LOLE, count that as one event, - ## so take the max LOLE across the zones - dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max() ## Get the loss-of-load events, keeping the max hourly probability for each ## (not really right probabilistically) lole = pd.Series({ @@ -217,24 +211,21 @@ def calc_lole(dflole, rmap, threshold=0): return lole -def calc_max_duration(dfeue, rmap, threshold=0): +def calc_max_duration(dfeue_agg, threshold=0): """Max event duration, where an event is >threshold EUE [MW] in contiguous hours""" - dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum() max_duration = pd.Series({ r: get_events(dfeue_agg[r], threshold)['timesteps'].max() for r in dfeue_agg }) return max_duration -def calc_neue(dfeue, dfload, rmap): +def calc_neue(dfeue_agg, dfload_agg): """NEUE (sum of EUE / sum of load) in units of ppm""" - dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum() - dfload_agg = dfload.rename(columns=rmap).groupby(axis=1, level=0).sum() neue = dfeue_agg.sum() / dfload_agg.sum() * 1e6 return neue -def calc_peak_eue(dfeue, dfload, rmap, norm:Literal['peak','hourly','absolute']='peak'): +def calc_peak_eue(dfeue_agg, dfload_agg, norm:Literal['peak','hourly','absolute']='peak'): """ Get the peak hourly outage magnitude @@ -244,15 +235,13 @@ def calc_peak_eue(dfeue, dfload, rmap, norm:Literal['peak','hourly','absolute']= - 'hourly': Divide hourly EUE by hourly load -> returns fraction - 'absolute': Do not normalize -> returns MW """ - dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum() - dfload_agg = dfload.rename(columns=rmap).groupby(axis=1, level=0).sum() match norm: case 'peak': peak_eue = dfeue_agg.max() / dfload_agg.max() case 'hourly': peak_eue = (dfeue_agg / dfload_agg).max() case 'absolute': - peak_eue = dfeue.max() + peak_eue = dfeue_agg.max() return peak_eue @@ -289,16 +278,21 @@ def calc_ra_metrics( ra_metrics = {} for hierarchy_level in levels: print(f'Calculating RA metrics at {hierarchy_level} level') - ## Get the region aggregator + ### Aggregate the shortfall and load to this hierarchy level rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) + ## If multiple zones in one level and hour have LOLE, count that as one event, + ## so take the max LOLE across the zones + 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() ## Calculate the full-timeseries metrics for each region - ra_metrics[hierarchy_level, 'lold_peryear'] = calc_lold(dflole, rmap) / numyears - ra_metrics[hierarchy_level, 'lole_peryear'] = calc_lole(dflole, rmap) / numyears - ra_metrics[hierarchy_level, 'max_duration'] = calc_max_duration(dfeue, rmap) - ra_metrics[hierarchy_level, 'neue_ppm'] = calc_neue(dfeue, dfload, rmap) - ra_metrics[hierarchy_level, 'euemax_peakloadfrac'] = calc_peak_eue(dfeue, dfload, rmap, 'peak') - ra_metrics[hierarchy_level, 'euemax_hourlyloadfrac'] = calc_peak_eue(dfeue, dfload, rmap, 'hourly') - ra_metrics[hierarchy_level, 'euemax_mw'] = calc_peak_eue(dfeue, dfload, rmap, 'absolute') + ra_metrics[hierarchy_level, 'lold_peryear'] = calc_lold(dflole_agg) / numyears + ra_metrics[hierarchy_level, 'lole_peryear'] = calc_lole(dflole_agg) / numyears + ra_metrics[hierarchy_level, 'max_duration'] = calc_max_duration(dfeue_agg) + ra_metrics[hierarchy_level, 'neue_ppm'] = calc_neue(dfeue_agg, dfload_agg) + ra_metrics[hierarchy_level, 'euemax_peakloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'peak') + ra_metrics[hierarchy_level, 'euemax_hourlyloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'hourly') + ra_metrics[hierarchy_level, 'euemax_mw'] = calc_peak_eue(dfeue_agg, dfload_agg, 'absolute') ### Combine it dfout = pd.concat(ra_metrics, names=['level','metric','region']).rename('value') From dc8c59e2dfa27b481fc79044c1702480e6f19566 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Tue, 23 Jun 2026 17:52:50 -0600 Subject: [PATCH 05/24] reformat RA threshold switches --- cases.csv | 14 +++++++------- runreeds.py | 52 ++++++++++++++++++---------------------------------- 2 files changed, 25 insertions(+), 41 deletions(-) diff --git a/cases.csv b/cases.csv index 39cdf9e2..4cc03fca 100644 --- a/cases.csv +++ b/cases.csv @@ -270,13 +270,13 @@ GSw_PRM_StressOutages,Whether to apply the availability factor (forced + schedul GSw_PRM_StressSeedLoadLevel,Region hierarchy level at which to include peak coincident load days as seeded stress periods (or False to ignore peaks),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,transgrp, GSw_PRM_StressSeedMinRElevel,Region hierarchy level at which to include minimum wind and solar capacity factor days as seeded stress periods (or False to ignore min-wind/solar CF days),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,interconnect, GSw_PRM_StressStorageCutoff,"How to select ""shoulder"" stress periods giving storage time to recharge before/after high-unserved-energy periods. Two-part switch separated by _. The first argument is ""EUE"" or ""capacity"" or ""absolute"". If ""EUE"" the second argument specifies a PRAS storage headspace / EUE threshold [fraction]; if ""cap"" it specifies a headspace / storage capacity threshold [fraction]; if ""abs"" it specifies absolute number of periods before/after [integer]. Turned off if set to ""off"".",N/A,EUE_0.1, -GSw_PRM_StressThresholdMetrics,"/-delimited list of metrics for identifying stress periods (supported options are: NEUE, LOLH)",N/A,NEUE, -GSw_PRM_StressThresholdLOLH,LOLH threshold [hours/year] above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_LOLH_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; LOLH is loss of load event-hours per year [event-h/year]; 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_2.4_sum, -GSw_PRM_StressThresholdLOLE,LOLE threshold [event-days/year] above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_LOLEdays_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; LOLEdays is loss of load event-days per year (a day is counted if at least one hour has a shortfall) [event-day/year]; 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_0.1_max, -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_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; NEUEppm is normalized expected unserved energy in parts per million [ppm]; 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_sum, -GSw_PRM_StressThresholdOutageDuration,Outage duration; formulated as HierarchyLevel_OutageDurationHours_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; OutageDurationHours is the max outage duration in hours; PeriodAggMethod is 'max' over the hours in each period (only used in period selection) (see README.md for detailed notes),N/A,transgrp_10000_max, -GSw_PRM_StressThresholdOutageMagnitude,Outage magnitude; formulated as HierarchyLevel_OutageMagnitudeMW_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; OutageMagnitudeMW is the max outage magnitude in MW; PeriodAggMethod 'max' over the hours in each period (only used in period selection) (see README.md for detailed notes),N/A,transgrp_0.1_max, -GSw_PRM_StressThresholdNormalizedOutageMagnitude,Normalized outage magnitude; formulated as HierarchyLevel_NormalizedOutageMagnitude_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; NormalizedOutageMagnitude is the max outage normalized magnitude (ratio of the maxmimum load) ; PeriodAggMethod is 'max' over the hours in each period (only used in period selection) (see README.md for detailed notes),N/A,transgrp_10_max, +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, +GSw_PRM_StressThresholdDepth,Outage depth threshold [MW_EUE/MW_peak_load] (fraction); formulated as HierarchyLevel_Depth where HierarchyLevel is a column in hierarchy.csv; Depth is the max outage magnitude in [MW EUE] / [MW peak demand],N/A,transgrp_0.1, +GSw_PRM_StressThresholdDuration,Outage duration threshold [hours]; formulated as HierarchyLevel_Duration where HierarchyLevel is a column in hierarchy.csv; Duration is the max outage duration in hours,N/A,transgrp_12, +GSw_PRM_StressThresholdLOLD,LOLD threshold [event-days/year]; formulated as HierarchyLevel_LOLD where HierarchyLevel is a column in hierarchy.csv; LOLD is loss-of-load event-days per year,N/A,transgrp_0.1, +GSw_PRM_StressThresholdLOLE,LOLE threshold [events/year]; formulated as HierarchyLevel_LOLE where HierarchyLevel is a column in hierarchy.csv; events is loss-of-load events per year,N/A,transgrp_0.1, +GSw_PRM_StressThresholdLOLH,LOLH threshold [event-hours/year]; formulated as HierarchyLevel_LOLH where HierarchyLevel is a column in hierarchy.csv; LOLH is loss-of-load event-hours per year,N/A,transgrp_2.4, +GSw_PRM_StressThresholdNEUE,NEUE threshold [ppm]; formulated as HierarchyLevel_NEUE where HierarchyLevel is a column in hierarchy.csv; NEUEppm is normalized expected unserved energy in parts per million,N/A,transgrp_1, 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, diff --git a/runreeds.py b/runreeds.py index 15be4348..41e84821 100644 --- a/runreeds.py +++ b/runreeds.py @@ -318,48 +318,32 @@ def check_compatibility(sw): hierarchy = reeds.io.get_hierarchy(GSw_ZoneSet=sw['GSw_ZoneSet']).reset_index() ### Check that the stress metrics specified in GSw_PRM_StressThresholdMetrics - # have corresponding GSw_PRM_StressThreshold{metric} switches - ## TODO: Check a regex way to list the switches - stressThresholdMetricSwitches = ['NEUE','LOLH','LOLE', 'OutageDuration', 'OutageMagnitude', 'NormalizedOutageMagnitude'] - UpperstressThresholdMetricSwitches = [s.upper() for s in stressThresholdMetricSwitches] - stressThresholdMetrics = sw['GSw_PRM_StressThresholdMetrics'].split('/') - - # if int(sw['GSw_PRM_UpdateMethod']) == 1 or 'EUE' not in stressThresholdMetrics: - # if int(sw['GSw_PRM_UpdateMethod']) > 1: - # raise NotImplementedError( - # f"Warning: GSw_PRM_UpdateMethod is set to {sw['GSw_PRM_UpdateMethod']}, " - # "but EUE is not included in GSw_PRM_StressThresholdMetrics, so defaulting to fixed increment. " - # "Add EUE to GSw_PRM_StressThresholdMetrics to use PRAS-informed update method." - # ) - - # Threshold metric added but not specified as a switch - for metric in stressThresholdMetrics: - if metric.upper() not in UpperstressThresholdMetricSwitches: - raise NotImplementedError(f"GSw_PRM_StressThreshold{metric} is not specified as a switch.") - - for stress_metric in stressThresholdMetricSwitches: - for threshold in sw[f'GSw_PRM_StressThreshold{stress_metric}'].split('/'): - ## Example: NEUE threshold = 'transgrp_1_sum' - allowed_levels = ['country','interconnect','nercr','transreg','transgrp','st','r'] - (hierarchy_level, stress_value, period_agg_method) = threshold.split('_') + ### are allowed and have well-formed GSw_PRM_StressThreshold{metric} entries + ra_switches = { + i.lower(): f'GSw_PRM_StressThreshold{i}' + for i in ['Depth', 'Duration', 'LOLD', 'LOLE', 'LOLH', 'NEUE'] + } + used_metrics = [i.lower() for i in sw['GSw_PRM_StressThresholdMetrics'].split('/')] + allowed_levels = ['country','interconnect','nercr','transreg','transgrp','st','r'] + + for metric in used_metrics: + if metric not in ra_switches: + raise NotImplementedError(f"GSw_PRM_StressThresholdMetrics = {metric} is not supported") + + for threshold in sw[ra_switches[metric]].split('/'): + ## Example: GSw_PRM_StressThresholdNEUE = 'transgrp_1' + (hierarchy_level, stress_value) = threshold.split('_') if hierarchy_level not in allowed_levels: raise ValueError( - f"GSw_PRM_StressThreshold{stress_metric}: level={hierarchy_level} but must be in:\n" + f"{ra_switches[metric]}: level={hierarchy_level} but must be in:\n" + '\n'.join(allowed_levels) ) - if period_agg_method.lower() not in ['sum','max']: - raise ValueError(f"Fix period agg method in GSw_PRM_StressThreshold{stress_metric}") if not (float(stress_value) >= 0): raise ValueError( - f"stress value in GSw_PRM_StressThreshold{stress_metric} must be a positive number " + f"stress value in {ra_switches[metric]} must be a positive number " f"but '{stress_value}' was provided" ) - # if stress_metric.upper() not in UpperstressThresholdMetricSwitches: - # raise ValueError( - # f"stress metric in GSw_PRM_StressThreshold{stress_metric} must be in {stressThresholdMetricSwitches}" - # f"but '{stress_metric}' was provided" - # ) - + if sw['GSw_PRM_StressStorageCutoff'].lower() not in ['off','0','false']: metric, value = sw['GSw_PRM_StressStorageCutoff'].split('_') if metric.lower()[:3] not in ['eue', 'cap', 'abs']: From 7f4e779f41a7adc5483d5e179d28a2da379d364a Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Tue, 23 Jun 2026 18:11:15 -0600 Subject: [PATCH 06/24] stress_periods.py: write EUE event list --- reeds/resource_adequacy/stress_periods.py | 64 +++++++++++++++-------- 1 file changed, 43 insertions(+), 21 deletions(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index f1ea5b65..e8743e39 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -276,23 +276,23 @@ def calc_ra_metrics( ### Loop over aggregation levels and calculate all metrics ra_metrics = {} - for hierarchy_level in levels: - print(f'Calculating RA metrics at {hierarchy_level} level') + for level in levels: + print(f'Calculating RA metrics at {level} level') ### Aggregate the shortfall and load to this hierarchy level - rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) + rmap = reeds.io.get_rmap(case=case, hierarchy_level=level) ## If multiple zones in one level and hour have LOLE, count that as one event, ## so take the max LOLE across the zones 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() ## Calculate the full-timeseries metrics for each region - ra_metrics[hierarchy_level, 'lold_peryear'] = calc_lold(dflole_agg) / numyears - ra_metrics[hierarchy_level, 'lole_peryear'] = calc_lole(dflole_agg) / numyears - ra_metrics[hierarchy_level, 'max_duration'] = calc_max_duration(dfeue_agg) - ra_metrics[hierarchy_level, 'neue_ppm'] = calc_neue(dfeue_agg, dfload_agg) - ra_metrics[hierarchy_level, 'euemax_peakloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'peak') - ra_metrics[hierarchy_level, 'euemax_hourlyloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'hourly') - ra_metrics[hierarchy_level, 'euemax_mw'] = calc_peak_eue(dfeue_agg, dfload_agg, 'absolute') + ra_metrics[level, 'lold_peryear'] = calc_lold(dflole_agg) / numyears + ra_metrics[level, 'lole_peryear'] = calc_lole(dflole_agg) / numyears + ra_metrics[level, 'max_duration'] = calc_max_duration(dfeue_agg) + ra_metrics[level, 'neue_ppm'] = calc_neue(dfeue_agg, dfload_agg) + ra_metrics[level, 'euemax_peakloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'peak') + ra_metrics[level, 'euemax_hourlyloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'hourly') + ra_metrics[level, 'euemax_mw'] = calc_peak_eue(dfeue_agg, dfload_agg, 'absolute') ### Combine it dfout = pd.concat(ra_metrics, names=['level','metric','region']).rename('value') @@ -300,6 +300,25 @@ def calc_ra_metrics( return dfout +def get_eue_events( + case:str|Path, + t:int, + iteration:int=0, + levels=['country', 'interconnect', 'nercr', 'transreg', 'transgrp', 'st', 'r'], +): + ### Get values from PRAS + dfeue = get_pras_shortfall(case, t, iteration)['EUE'] + + ### Get the list of events at all hierarchy levels + events = {} + for level in levels: + rmap = reeds.io.get_rmap(case=case, hierarchy_level=level) + dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum() + events[level] = pd.concat({r: get_events(dfeue_agg[r]) for r in dfeue_agg}) + dfout = pd.concat(events, names=['level','region','number']) + return dfout + + def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods, stress_metric): ## Stop if not needed if sw.GSw_PRM_StressStorageCutoff.lower() in ['off', '0', 'false']: @@ -493,7 +512,7 @@ def _evaluate_stress_threshold_criterion( return stress_criteria -def get_stress_metrics_sorted_periods(sw, t, iteration): +def get_stress_periods(sw, t, iteration): ### Get storage state of charge (SOC) to use in selection of "shoulder" stress periods dfenergy = reeds.io.read_pras_results( os.path.join(sw['casedir'], 'handoff', 'PRAS', f"PRAS_{t}i{iteration}-energy.h5") @@ -527,13 +546,12 @@ def get_stress_metrics_sorted_periods(sw, t, iteration): # stress periods column names for writing outputs stress_metric_units = { - 'EUE':'MWh/year', - 'NEUE':'ppm', - 'LOLH':'event-h/year', - 'LOLE':'event-day/year', - 'OutageDuration':'h', - 'OutageMagnitude':'MW', - 'NormalizedOutageMagnitude':'p.u. of load', + 'NEUE': 'ppm', + 'LOLD': 'event-days/year', + 'LOLH': 'event-hours/year', + 'LOLE': 'events/year', + 'Duration': 'hours', + 'Depth': 'fraction', } stress_metrics_col_names = { m: f'{m}_{stress_metric_units[m]}' for m in stress_metric_switches @@ -561,8 +579,6 @@ def get_stress_metrics_sorted_periods(sw, t, iteration): stress_sorted_periods = pd.concat(stress_sorted_periods, names=['criterion']) ### Get lists of stress periods: new (added this iteration) and all - ##TODO: What if same high stress period is added for multiple criteria? - # Currently the duplicates are dropped. if len(failed): new_stress_periods = pd.concat( {**high_stress_periods, **shoulder_periods}, names=['criterion','periodtype'], @@ -811,12 +827,18 @@ def main(sw, t, iteration=0, logging=True): os.path.join(sw.casedir, 'outputs', f'ra_metrics_{t}i{iteration}.csv') ) + #%% Write EUE events + eue_events = get_eue_events(case=sw.casedir, t=t, iteration=iteration) + eue_events.round(3).to_csv( + os.path.join(sw.casedir, 'outputs', f'eue_events_{t}i{iteration}.csv') + ) + #%% Stop here if not iterating or if before ReEDS can build new capacity if (not int(sw.GSw_PRM_StressIterateMax)) or (t < int(sw['GSw_StartMarkets'])): return #%% Identify and write new stress periods - failed, new_stressperiods_write, combined_periods_write = get_stress_metrics_sorted_periods( + failed, new_stressperiods_write, combined_periods_write = get_stress_periods( sw=sw, t=t, iteration=iteration, ) From 7d6a1a29d039a4dd7ca21332a4a926770f2c90c0 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Tue, 23 Jun 2026 19:48:28 -0600 Subject: [PATCH 07/24] add MultiMetricRA case to cases_test.csv --- cases_test.csv | 131 +++++++++++++++++++++++++------------------------ 1 file changed, 66 insertions(+), 65 deletions(-) diff --git a/cases_test.csv b/cases_test.csv index a87de967..03af5156 100644 --- a/cases_test.csv +++ b/cases_test.csv @@ -1,65 +1,66 @@ -,Default Value,Pacific,USA_defaults,Mid_Case,USA_decarb,github_Pacific,github_Everything,github_MA_county_CC,Pacific_CC,Pacific_weks,Pacific_full_year,Interday_storage,Pacific_2020,Pacific_rep15,WY_county,WECC_county,PJM_county_CC,NYVT_mixed,OR_water,MonteCarlo_Random,MonteCarlo_LHS,Everything,Simple,USA_fast,USA_faster,Pacific_DR,Pacific_MGA,Pacific_LoadSite95,MARICTNYNJPAOH_Offshore,R2P -ignore,1,0,,,,,,,,,,,,,,,,,,,,,,,,,,,, -GSw_Region,cendiv/Pacific,,country/USA,country/USA,country/USA,,st/ID.WY.NE.IA.IL,st/MA,,,,,,,st/WY,interconnect/western,transreg/PJM,st/NY.VT,st/OR,st/NE.NY.PA,st/NE.NY.PA,st/ID.WY.NE.IA.IL,st/KS,country/USA,country/USA,,,,st/MA.RI.CT.NY.NJ.PA.OH, -endyear,2032,,2050,2050,2050,2029,2060,2026,,,,,,,,,,,2035,2030,2030,2060,2035,2050,2050,,,,, -yearset,,,,,,,2010..2060..10,,,,,,,,,,,,,2010..2050..5,2010..2050..5,2010..2060..10,,,2010_2025_2050,,,,, -GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,,z54,z48,,,,, -GSw_GasCurve,2,,1,1,,,,,,,,,,,,,,,,,,,,1,1,,,,, -GSw_Geothermal,,,,2,,,,,,,,,,,,,,,,,,,0,,0,,,,, -GSw_GrowthPenalties,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,, -GSw_Upstream,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,, -GSw_TransHurdleRate,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,, -distpvscen,,,,,stscen2023_mid_case_95_by_2035,,,,,,,,,,,,,,,,,,,,,,,,, -GSw_AnnualCap,,,,,2,,1,,,,,,,,,,,,,,,1,,,,,,,, -GSw_AnnualCapScen,,,,,start2024_90pct2035_100pct2045,,start2027_95pct2035,,,,,,,,,,,,,,,start2027_95pct2035,,,,,,,, -GSw_LoadProfiles,,,,,EER2025_100by2050,EER2025_IRAlow,EER2025_IRAlow,EER2025_IRAlow,,,,,historic,,,,,,,,,EER2025_100by2050,,,,historic,,,, -GSw_NG_CRF_penalty,,,,,ramp_2045,,ramp_2023_2035,,,,,,,,,,,,,,,ramp_2023_2035,,,,,,,, -GSw_PRM_NetImportLimit,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,, -GSw_RetirePenalty,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,, -GSw_FakeData,,,,,,1,1,1,,,,,,,,,,,,,,,,,,,,,, -GSw_PRM_CapCredit,,,,,,,,1,1,,,,,,,,1,,,,,,,,,,,,, -GSw_PRM_scenario,,,,,,,,,static,,,,,,,,static,,,,,,,,,,,,, -GSw_PRM_UpdateMethod,,,,,,,,,1,,,,,,,,,,,,,,,,,,,,, -GSw_HourlyType,,,,,,,,,,wek,year,,,,,,,,,,,,,,,,,,, -GSw_InterDayLinkage,,,,,,,,,,,,1,,,,,,,,,,,,,,,,,, -GSw_HourlyWeatherYears,,,,,,,2012_2013,,,,,,2020,2007_2008_2009_2010_2011_2012_2013_2016_2017_2018_2019_2020_2021_2022_2023,,,,,,,,2012_2013,,,,2018,,,, -GSw_HourlyClusterMapMethod,,,,,,,,,,,,,,bestfirst,,,,,,,,,,,,,,,, -GSw_WaterCapacity,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,, -GSw_WaterMain,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,, -GSw_WaterUse,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,, -resource_adequacy_years,,,,,,,2011_2012_2013_2021_2022_2023,,,,,,,,,,,,,,,2011_2012_2013_2021_2022_2023,,,,,,,, -GSw_HourlyClusterAlgorithm,,,,,,,,,,,,,,,,,,,,user,user,,,,,,,,, -MCS_runs,,,,,,,,,,,,,,,,,,,,2,2,,,,,,,,, -MCS_dist_groups,,,,,,,,,,,,,,,,,,,,tech.hydro.nuclear.gas.coal.load_country,upv_tri.nuclear_tri.ng_fuel_price_tri.load_country_unif,,,,,,,,, -MCS_lhs,,,,,,,,,,,,,,,,,,,,,1,,,,,,,,, -GSw_PRM_StressIterateMax,,,,,,,,0,,,,,,,,0,0,,,1,1,,,,,,,,, -GSw_ReducedResource,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,, -GSw_SitingUPV,,,,,,limited,limited,limited,,,,,,,,,,,,,,limited,,,,,,,, -GSw_SitingWindOfs,,,,,,limited,limited,limited,,,,,,,,,,,,,,open,,,,,,,, -GSw_SitingWindOns,,,,,,limited,limited,limited,,,,,,,,,,,,,,limited,,,,,,,, -GSw_TransScen,,,,,,,NTP_MT,,,,,,,,,,,,,,,NTP_MT,,,,,,,, -GSw_CO2_Detail,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,, -GSw_DAC,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,, -GSw_NoFossilOffsetCDR,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,, -GSw_Biopower,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,, -GSw_HourlyChunkLengthRep,,,,,,,,,,,,,,,,,,,,,,,6,4,4,,,,, -GSw_HourlyChunkLengthStress,,,,,,,,,,,,,,,,,,,,,,,6,4,4,,,,, -GSw_LfillGas,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,, -GSw_Nuclear,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,, -GSw_OpRes,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,, -GSw_StartCost,,,,,,,,,,,,,,,,,,,,,,,0,0,0,,,,, -GSw_H2,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,, -GSw_H2_PTC,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,, -GSw_H2Combustion,,,,,,,,,,,,,,,,,,,,,,,,,0,,,,, -GSw_DRShed,,,,,,,,,,,,,,,,,,,,,,,,,,1,,,, -GSw_MGA_CostDelta,,,,,,,,,,,,,,,,,,,,,,,,,,,0.01,,, -GSw_LoadSiteCF,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.95,, -GSw_OffshoreZones,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, -GSw_OffshoreBackbone,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, -GSw_OffshoreBackflow,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, -pras_agg_ogs_lfillgas,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1 -pras_existing_unit_size,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0 -pras_scheduled_outage,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0 -pras_unitsize_source,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,r2x -pras_vre_combine,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1 -pras_samples,,,,,,10,10,10,,,,,,,,,,,,,,,,10,10,,,,, +,Default Value,Pacific,USA_defaults,Mid_Case,USA_decarb,github_Pacific,github_Everything,github_MA_county_CC,Pacific_CC,Pacific_weks,Pacific_full_year,Interday_storage,Pacific_2020,Pacific_rep15,WY_county,WECC_county,PJM_county_CC,NYVT_mixed,OR_water,MonteCarlo_Random,MonteCarlo_LHS,Everything,Simple,USA_fast,USA_faster,MultiMetricRA,Pacific_DR,Pacific_MGA,Pacific_LoadSite95,MARICTNYNJPAOH_Offshore,R2P +ignore,1,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,, +GSw_Region,cendiv/Pacific,,country/USA,country/USA,country/USA,,st/ID.WY.NE.IA.IL,st/MA,,,,,,,st/WY,interconnect/western,transreg/PJM,st/NY.VT,st/OR,st/NE.NY.PA,st/NE.NY.PA,st/ID.WY.NE.IA.IL,st/KS,country/USA,country/USA,st/NY,,,,st/MA.RI.CT.NY.NJ.PA.OH, +endyear,2032,,2050,2050,2050,2029,2060,2026,,,,,,,,,,,2035,2030,2030,2060,2035,2050,2050,2050,,,,, +yearset,,,,,,,2010..2060..10,,,,,,,,,,,,,2010..2050..5,2010..2050..5,2010..2060..10,,,2010_2025_2050,2010..2050..5,,,,, +GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,,z54,z48,,,,,, +GSw_GasCurve,2,,1,1,,,,,,,,,,,,,,,,,,,,1,1,,,,,, +GSw_Geothermal,,,,2,,,,,,,,,,,,,,,,,,,0,,0,,,,,, +GSw_GrowthPenalties,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,, +GSw_Upstream,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,, +GSw_TransHurdleRate,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,, +distpvscen,,,,,stscen2023_mid_case_95_by_2035,,,,,,,,,,,,,,,,,,,,,,,,,, +GSw_AnnualCap,,,,,2,,1,,,,,,,,,,,,,,,1,,,,,,,,, +GSw_AnnualCapScen,,,,,start2024_90pct2035_100pct2045,,start2027_95pct2035,,,,,,,,,,,,,,,start2027_95pct2035,,,,,,,,, +GSw_LoadProfiles,,,,,EER2025_100by2050,EER2025_IRAlow,EER2025_IRAlow,EER2025_IRAlow,,,,,historic,,,,,,,,,EER2025_100by2050,,,,,historic,,,, +GSw_NG_CRF_penalty,,,,,ramp_2045,,ramp_2023_2035,,,,,,,,,,,,,,,ramp_2023_2035,,,,,,,,, +GSw_PRM_NetImportLimit,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,, +GSw_RetirePenalty,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,, +GSw_FakeData,,,,,,1,1,1,,,,,,,,,,,,,,,,,,,,,,, +GSw_PRM_CapCredit,,,,,,,,1,1,,,,,,,,1,,,,,,,,,,,,,, +GSw_PRM_scenario,,,,,,,,,static,,,,,,,,static,,,,,,,,,,,,,, +GSw_PRM_UpdateMethod,,,,,,,,,1,,,,,,,,,,,,,,,,,,,,,, +GSw_HourlyType,,,,,,,,,,wek,year,,,,,,,,,,,,,,,,,,,, +GSw_InterDayLinkage,,,,,,,,,,,,1,,,,,,,,,,,,,,,,,,, +GSw_HourlyWeatherYears,,,,,,,2012_2013,,,,,,2020,2007_2008_2009_2010_2011_2012_2013_2016_2017_2018_2019_2020_2021_2022_2023,,,,,,,,2012_2013,,,,,2018,,,, +GSw_HourlyClusterMapMethod,,,,,,,,,,,,,,bestfirst,,,,,,,,,,,,,,,,, +GSw_WaterCapacity,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,,, +GSw_WaterMain,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,,, +GSw_WaterUse,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,,, +resource_adequacy_years,,,,,,,2011_2012_2013_2021_2022_2023,,,,,,,,,,,,,,,2011_2012_2013_2021_2022_2023,,,,,,,,, +GSw_HourlyClusterAlgorithm,,,,,,,,,,,,,,,,,,,,user,user,,,,,,,,,, +MCS_runs,,,,,,,,,,,,,,,,,,,,2,2,,,,,,,,,, +MCS_dist_groups,,,,,,,,,,,,,,,,,,,,tech.hydro.nuclear.gas.coal.load_country,upv_tri.nuclear_tri.ng_fuel_price_tri.load_country_unif,,,,,,,,,, +MCS_lhs,,,,,,,,,,,,,,,,,,,,,1,,,,,,,,,, +GSw_PRM_StressIterateMax,,,,,,,,0,,,,,,,,0,0,,,1,1,,,,,,,,,, +GSw_ReducedResource,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,, +GSw_SitingUPV,,,,,,limited,limited,limited,,,,,,,,,,,,,,limited,,,,,,,,, +GSw_SitingWindOfs,,,,,,limited,limited,limited,,,,,,,,,,,,,,open,,,,,,,,, +GSw_SitingWindOns,,,,,,limited,limited,limited,,,,,,,,,,,,,,limited,,,,,,,,, +GSw_TransScen,,,,,,,NTP_MT,,,,,,,,,,,,,,,NTP_MT,,,,,,,,, +GSw_CO2_Detail,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,, +GSw_DAC,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,, +GSw_NoFossilOffsetCDR,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,, +GSw_Biopower,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,, +GSw_HourlyChunkLengthRep,,,,,,,,,,,,,,,,,,,,,,,6,4,4,,,,,, +GSw_HourlyChunkLengthStress,,,,,,,,,,,,,,,,,,,,,,,6,4,4,,,,,, +GSw_LfillGas,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,, +GSw_Nuclear,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,, +GSw_OpRes,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,, +GSw_StartCost,,,,,,,,,,,,,,,,,,,,,,,0,0,0,,,,,, +GSw_H2,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,, +GSw_H2_PTC,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,, +GSw_H2Combustion,,,,,,,,,,,,,,,,,,,,,,,,,0,,,,,, +GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth,,,,, +GSw_DRShed,,,,,,,,,,,,,,,,,,,,,,,,,,,1,,,, +GSw_MGA_CostDelta,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.01,,, +GSw_LoadSiteCF,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.95,, +GSw_OffshoreZones,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, +GSw_OffshoreBackbone,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, +GSw_OffshoreBackflow,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, +pras_agg_ogs_lfillgas,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1 +pras_existing_unit_size,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0 +pras_scheduled_outage,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0 +pras_unitsize_source,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,r2x +pras_vre_combine,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1 +pras_samples,,,,,,10,10,10,,,,,,,,,,,,,,,,10,10,,,,,, From f9da3338532d194af425e6c3d095e63cc6153451 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Tue, 23 Jun 2026 20:10:10 -0600 Subject: [PATCH 08/24] stress_periods.py: fix get_events() and calc_max_duration() for no events --- reeds/resource_adequacy/stress_periods.py | 36 ++++++++++++++++------- 1 file changed, 26 insertions(+), 10 deletions(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index e8743e39..c7e75912 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -16,6 +16,13 @@ # importlib.reload(functions) +#%%### Constants +RA_SWITCHES = { + i.lower(): f'GSw_PRM_StressThreshold{i}' + for i in ['Depth', 'Duration', 'LOLD', 'LOLE', 'LOLH', 'NEUE'] +} + + #%%### Functions def get_pras_shortfall(case, t, iteration=0): """ @@ -176,17 +183,21 @@ def get_events(ds:pd.Series, threshold:float=0) -> pd.DataFrame: ends = ends.loc[ends > 0].index assert len(starts) == len(ends), "Error in event start/end calculation" ## Get some metrics for each event - dfout = [] + events = [] for start, end in zip(starts, ends): event = ds.loc[start:end] - dfout.append({ + events.append({ 'start': start, 'end': end, 'timesteps': len(event), 'max': event.max(), 'sum': event.sum(), }) - return pd.DataFrame(dfout) + if len(events): + dfout = pd.DataFrame(events) + else: + dfout = pd.DataFrame(columns=['start','end','timesteps','max','sum']) + return dfout def calc_lold(dflole_agg, threshold=0): @@ -215,7 +226,7 @@ def calc_max_duration(dfeue_agg, threshold=0): """Max event duration, where an event is >threshold EUE [MW] in contiguous hours""" max_duration = pd.Series({ r: get_events(dfeue_agg[r], threshold)['timesteps'].max() for r in dfeue_agg - }) + }).fillna(0).astype(int) return max_duration @@ -616,10 +627,18 @@ def get_stress_periods(sw, t, iteration): combined_periods_write.to_csv(outpath, index=False) ### Tables and plots for debugging - stress_sorted_periods.round(2).rename(columns=stress_metrics_col_names).to_csv( + stress_metric_labels = { + 'NEUE': 'neue_ppm', + 'LOLD': 'lold_event-days/year', + 'LOLH': 'lolh_event-hours/year', + 'LOLE': 'lole_events/year', + 'Duration': 'duration_hours', + 'Depth': 'depth_fraction', + } + stress_sorted_periods.round(2).rename(columns=stress_metric_labels).to_csv( os.path.join(sw.casedir, 'inputs_case', newstresspath, 'stress_metrics_sorted_periods.csv') ) - new_stress_periods.round(2).rename(columns=stress_metrics_col_names).to_csv( + new_stress_periods.round(2).rename(columns=stress_metric_labels).to_csv( os.path.join(sw.casedir, 'inputs_case', newstresspath, 'new_stress_periods.csv'), index=False, ) @@ -877,12 +896,9 @@ def main(sw, t, iteration=0, logging=True): os.path.join(sw.casedir, 'inputs_case', newstresspath, 'prm.csv'), ) - #%% Done - return - +# #%%### Option to run script directly for debugging # if __name__ == '__main__': -# #%%### option to run script directly for debugging # casedir = "/path/to/ReEDS/runs/runname" # t = 2030 # previous solve year # iteration = 0 From cf7b5c84220ca1355defa5ca1b8079935fd81f1a Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Wed, 24 Jun 2026 09:41:07 -0600 Subject: [PATCH 09/24] replace get_stress_metric_periods() with get_shortfall_periods() and get_longest_events() --- reeds/resource_adequacy/stress_periods.py | 450 ++++++++++------------ 1 file changed, 202 insertions(+), 248 deletions(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index c7e75912..1f92fe82 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -21,6 +21,14 @@ i.lower(): f'GSw_PRM_StressThreshold{i}' for i in ['Depth', 'Duration', 'LOLD', 'LOLE', 'LOLH', 'NEUE'] } +SWITCH_METRIC = { + 'depth': 'euemax_peakloadfrac', + 'duration': 'max_duration', + 'lold': 'lold_peryear', + 'lole': 'lole_peryear', + 'lolh': 'lolh_peryear', + 'neue': 'neue_ppm', +} #%%### Functions @@ -55,83 +63,6 @@ def get_pras_shortfall(case, t, iteration=0): return dictout -def get_stress_metric_periods( - case, t, iteration=0, - hierarchy_level='transgrp', - stress_metric='EUE', - period_agg_method='sum', - ): - """_summary_ - - Args: - sw (pd.series): ReEDS switches for this run. - t (int): Model solve year. - iteration (int, optional): Iteration number of this solve year. Defaults to 0. - hierarchy_level (str, optional): column of hierarchy.csv specifying the spatial - level over which to calculate stress_metric. Defaults to 'country'. - stress_metric (str, optional): 'EUE' or 'NEUE'. Defaults to 'EUE'. - period_agg_method (str, optional): 'sum' or 'max', indicating how to aggregate - over the hours in each period. Defaults to 'sum'. - - Raises: - NotImplementedError: if invalid value for stress_metric or GSw_PRM_StressModel - - Returns: - pd.DataFrame: Table of periods sorted in descending order by stress metric. - """ - sw = reeds.io.get_switches(case) - ### Get the region aggregator - rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) - - ### Get stress metric from PRAS - # Use EUE metric for both EUE and NEUE calculations, since the load division to get NEUE is - # peformed afterwards based on agg_period. - use_metric_for_pras = { - 'EUE': 'EUE', - 'NEUE': 'EUE', - 'LOLH': 'LOLE', - 'LOLE': 'LOLE', - 'OutageDuration': 'EUE', - 'OutageMagnitude': 'EUE', - 'NormalizedOutageMagnitude': 'EUE', - } - dfmetric = ( - get_pras_shortfall(case=case, t=t, iteration=iteration) - [use_metric_for_pras['stress_metric']] - ) - ## Aggregate to hierarchy_level - dfmetric = ( - dfmetric - .rename_axis('r', axis=1).rename_axis('h', axis=0) - .rename(columns=rmap).groupby(axis=1, level=0).sum() - ) - - ###### Calculate the stress metric by period - dfmetric_period = ( - dfmetric - .groupby([dfmetric.index.year, dfmetric.index.month, dfmetric.index.day]) - .agg(period_agg_method) - .rename_axis(['y','m','d']) - ) - - ### Sort and drop zeros and duplicates - dfmetric_top = ( - dfmetric_period.stack('r') - .sort_values(ascending=False) - .replace(0,np.nan).dropna() - .reset_index().drop_duplicates(['y','m','d'], keep='first') - .set_index(['y','m','d','r']).squeeze(1).rename(stress_metric) - .reset_index('r') - ) - ## Convert to timestamp, then to ReEDS period - dfmetric_top['actual_period'] = [ - reeds.timeseries.timestamp2h(pd.Timestamp(*d), sw['GSw_HourlyType']).split('h')[0] - for d in dfmetric_top.index.values - ] - - return dfmetric_top - - def plot_stress_diagnostics(sw, t, iteration, high_stress_periods): try: dates = ( @@ -204,9 +135,7 @@ def calc_lold(dflole_agg, threshold=0): """Count a day as an event-day if at least one hour has LOLE > threshold""" ## Take the max for each day ## (That's not quite right if the events are independent) - daily_max = dflole_agg.groupby( - [dflole_agg.index.year, dflole_agg.index.month, dflole_agg.index.day] - ).max() + daily_max = dflole_agg.resample('D').max() ## Sum the probability across days lold = daily_max.sum() return lold @@ -222,6 +151,14 @@ def calc_lole(dflole_agg, threshold=0): return lole +def calc_lolh(dflole_agg): + """Event-hours, where an event-hour is >threshold LOLE""" + ## Get the loss-of-load events, keeping the max hourly probability for each + ## (not really right probabilistically) + lolh = dflole_agg.sum() + return lolh + + def calc_max_duration(dfeue_agg, threshold=0): """Max event duration, where an event is >threshold EUE [MW] in contiguous hours""" max_duration = pd.Series({ @@ -299,6 +236,7 @@ def calc_ra_metrics( ## Calculate the full-timeseries metrics for each region ra_metrics[level, 'lold_peryear'] = calc_lold(dflole_agg) / numyears ra_metrics[level, 'lole_peryear'] = calc_lole(dflole_agg) / numyears + ra_metrics[level, 'lolh_peryear'] = calc_lolh(dflole_agg) / numyears ra_metrics[level, 'max_duration'] = calc_max_duration(dfeue_agg) ra_metrics[level, 'neue_ppm'] = calc_neue(dfeue_agg, dfload_agg) ra_metrics[level, 'euemax_peakloadfrac'] = calc_peak_eue(dfeue_agg, dfload_agg, 'peak') @@ -330,7 +268,53 @@ def get_eue_events( return dfout -def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods, stress_metric): +def get_shortfall_periods( + dsmetric:pd.Series, + aggmethod:Literal['max','sum']='sum', + GSw_HourlyType:Literal['day','wek','year']='day', +): + ## Keep days with nonzero metric + metric_day = dsmetric.resample('D').agg(aggmethod).replace(0, np.nan).dropna() + ## Convert to ReEDS period + metric_day.index = metric_day.index.map( + lambda x: reeds.timeseries.timestamp2h(x, GSw_HourlyType).split('h')[0] + ) + ## Sort by value + metric_period = metric_day.groupby(metric_day.index).sum().sort_values(ascending=False) + return metric_period + + +def get_longest_events( + sw, + t:int, + iteration:int, + hierarchy_level:str, + region:str, + num_events:int=1, +): + ## Get the already-identified events + fpath = Path(sw.casedir, 'outputs', f'eue_events_{t}i{iteration}.csv') + eue_events = ( + pd.read_csv(fpath, index_col=['level','region','number']) + .loc[hierarchy_level].loc[region] + .sort_values('timesteps', ascending=False) + .head(num_events) + ) + ## Get the datetimes in each event and keep all unique + dates = [] + for i, row in eue_events.iterrows(): + dates.append( + pd.Series(index=pd.date_range(row.start, row.end, freq='H'), data=1) + .resample('D').count() + ) + metric_period = pd.concat(dates) + metric_period.index = metric_period.index.map( + lambda x: reeds.timeseries.timestamp2h(x, sw.GSw_HourlyType).split('h')[0] + ) + return metric_period.groupby(level=0).sum() + + +def get_shoulder_periods(sw, criterion, dfenergy_agg, high_stress_periods, stress_metric): ## Stop if not needed if sw.GSw_PRM_StressStorageCutoff.lower() in ['off', '0', 'false']: print( @@ -338,7 +322,7 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods, stress_met "so not adding shoulder stress periods based on storage level" ) return {} - if dfenergy_r.empty: + if dfenergy_agg.empty: print( "No storage capacity, so no shoulder stress periods will be added " "based on storage level" @@ -346,34 +330,27 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods, stress_met return {} ## Parse inputs - hierarchy = reeds.io.get_hierarchy(sw.casedir) timeindex = reeds.timeseries.get_timeindex(sw['resource_adequacy_years']) cutofftype, cutoff = sw.GSw_PRM_StressStorageCutoff.lower().split('_') periodhours = {'day':24, 'wek':24*5, 'year':24}[sw.GSw_HourlyType] - (hierarchy_level, stress_value, period_agg_method) = criterion.split('_') - ## Aggregate storage energy to hierarchy_level - dfenergy_agg = ( - dfenergy_r.rename(columns=hierarchy[hierarchy_level]) - .groupby(axis=1, level=0).sum() - ) dfheadspace_MWh = dfenergy_agg.max() - dfenergy_agg dfheadspace_frac = dfheadspace_MWh / dfenergy_agg.max() shoulder_periods = {} - for i, row in high_eue_periods[criterion, f'high_{stress_metric}'].iterrows(): - if row.r not in dfheadspace_MWh: + for i, row in high_stress_periods.iterrows(): + if row.region not in dfheadspace_MWh: continue - day = pd.Timestamp('-'.join(row[['y','m','d']].astype(str).tolist())) + day = reeds.timeseries.h2timestamp(row.period) - start_headspace_MWh = dfheadspace_MWh.loc[day.strftime('%Y-%m-%d'),row.r].iloc[0] - end_headspace_MWh = dfheadspace_MWh.loc[day.strftime('%Y-%m-%d'),row.r].iloc[-1] + start_headspace_MWh = dfheadspace_MWh.loc[day.strftime('%Y-%m-%d'), row.region].iloc[0] + end_headspace_MWh = dfheadspace_MWh.loc[day.strftime('%Y-%m-%d'), row.region].iloc[-1] - start_headspace_frac = dfheadspace_frac.loc[day.strftime('%Y-%m-%d'),row.r].iloc[0] - end_headspace_frac = dfheadspace_frac.loc[day.strftime('%Y-%m-%d'),row.r].iloc[-1] + start_headspace_frac = dfheadspace_frac.loc[day.strftime('%Y-%m-%d'), row.region].iloc[0] + end_headspace_frac = dfheadspace_frac.loc[day.strftime('%Y-%m-%d'), row.region].iloc[-1] - day_eue = high_eue_periods[criterion, f'high_{stress_metric}'].loc[i, stress_metric] + day_eue = high_stress_periods.loc[i, stress_metric] day_index = np.where( timeindex == dfenergy_agg.loc[day.strftime('%Y-%m-%d')].iloc[0].name )[0][0] @@ -388,7 +365,7 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods, stress_met ): shoulder_periods[criterion, f'after_{row.name}'] = pd.Series({ 'actual_period':day_after.strftime('y%Yd%j'), - 'y':day_after.year, 'm':day_after.month, 'd':day_after.day, 'r':row.r, + 'region':row.region, }).to_frame().T.set_index('actual_period') print(f"Added {day_after} as shoulder stress period after {day}") @@ -399,141 +376,134 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods, stress_met ): shoulder_periods[criterion, f'before_{row.name}'] = pd.Series({ 'actual_period':day_before.strftime('y%Yd%j'), - 'y':day_before.year, 'm':day_before.month, 'd':day_before.day, 'r':row.r, + 'region':row.region, }).to_frame().T.set_index('actual_period') print(f"Added {day_before} as shoulder stress period before {day}") return shoulder_periods -def _evaluate_stress_threshold_criterion( - stress_criteria, - criterion, +def check_threshold_and_choose_periods( + stress_metric:str, + criterion:str, sw, - t, - iteration, - dfenergy_r, + t:int, + iteration:int, + dfeue_agg, + dflole_agg, + dfenergy_agg, stressperiods_this_iteration, - stress_metric, ): - _stress_sorted_periods = stress_criteria['stress_sorted_periods'] - _failed = stress_criteria['failed'] - _high_stress_periods = stress_criteria['high_stress_periods'] - _shoulder_periods = stress_criteria['shoulder_periods'] - - ## NEUE Example: criterion = 'transgrp_1_sum' - (hierarchy_level, stress_value, period_agg_method) = criterion.split('_') + ## NEUE Example: criterion = 'transgrp_1' + hierarchy_level, metric_threshold = criterion.split('_') + metric_threshold = float(metric_threshold) + GSw_HourlyType = sw.GSw_HourlyType ### Get stored stress metric - stress_vals = pd.read_csv( - os.path.join(sw.casedir, 'outputs', f'{stress_metric.lower()}_{t}i{iteration}.csv'), + ra_metrics = pd.read_csv( + os.path.join(sw.casedir, 'outputs', f'ra_metrics_{t}i{iteration}.csv'), index_col=['level', 'metric', 'region'], ).squeeze(1) - stress_periods = get_stress_metric_periods( - case=sw.casedir, t=t, iteration=iteration, - hierarchy_level=hierarchy_level, - stress_metric=stress_metric, - period_agg_method=period_agg_method, - ) - - ### Sort in descending stress_metric order - _stress_sorted_periods[criterion] = ( - stress_periods - .sort_values(stress_metric, ascending=False) - .reset_index().set_index('actual_period') - ) - ### Get the threshold(s) and see if any of them failed - this_test = stress_vals[hierarchy_level][period_agg_method] - - if (this_test > float(stress_value)).any(): - _failed[criterion] = this_test.loc[this_test > float(stress_value)] - print(f"GSw_PRM_StressThreshold = {criterion} failed for:") - print(_failed[criterion]) - ###### Add GSw_PRM_StressIncrement periods to the list for the next iteration - _high_stress_periods[criterion, f'high_{stress_metric}'] = ( - _stress_sorted_periods[criterion].loc[ - ## Only include new stress periods for the region(s) that failed - _stress_sorted_periods[criterion].r.isin(_failed[criterion].index) - ## Don't repeat existing stress periods - & ~(_stress_sorted_periods[criterion].index.isin( - stressperiods_this_iteration.actual_period)) - ] - ## Don't add dates more than once - .drop_duplicates(subset=['y','m','d']) - ## Keep the GSw_PRM_StressIncrement worst periods for each region. - ## If you instead want to keep the GSw_PRM_StressIncrement worst periods - ## overall, use .nlargest(int(sw.GSw_PRM_StressIncrement), stress_metric) - .groupby('r').head(int(sw.GSw_PRM_StressIncrement)) + this_test = ra_metrics[hierarchy_level][SWITCH_METRIC[stress_metric]] + failed = this_test.loc[this_test > metric_threshold] + if not len(failed): + print(f"{RA_SWITCHES[stress_metric]} = {criterion} passed") + else: + print(f"{RA_SWITCHES[stress_metric]} = {criterion} failed for:") + print(failed) + match stress_metric: + case 'depth': + metric_periods = { + region: ( + get_shortfall_periods(dfeue_agg[region], 'max', GSw_HourlyType) + .head(int(sw.GSw_PRM_StressIncrement)) + ) + for region in failed.index + } + case 'duration': + ## TODO: Keep whole events (including when they span days) + metric_periods = { + region: get_longest_events( + sw=sw, t=t, iteration=iteration, + hierarchy_level=hierarchy_level, region=region, + ) + for region in failed.index + } + case 'lold' | 'lole': + ## TODO: Double check this approach for LOLD and LOLE + metric_periods = { + region: ( + get_shortfall_periods(dflole_agg[region], 'max', GSw_HourlyType) + .head(int(sw.GSw_PRM_StressIncrement)) + ) + for region in failed.index + } + case 'lolh': + metric_periods = { + region: ( + get_shortfall_periods(dflole_agg[region], 'sum', GSw_HourlyType) + .head(int(sw.GSw_PRM_StressIncrement)) + ) + for region in failed.index + } + case 'neue': + metric_periods = { + region: ( + get_shortfall_periods(dfeue_agg[region], 'sum', GSw_HourlyType) + .head(int(sw.GSw_PRM_StressIncrement)) + ) + for region in failed.index + } + high_stress_periods = ( + pd.concat(metric_periods, names=['region','period']) + .rename(stress_metric) + .reset_index() ) - for period, row in _high_stress_periods[criterion, f'high_{stress_metric}'].iterrows(): + for i, row in high_stress_periods.iterrows(): print( - f"Added {period} " - f"({reeds.timeseries.h2timestamp(period).strftime('%Y-%m-%d')}) " - f"as stress period for {row.r} " + f"Added {row.period} " + f"({reeds.timeseries.h2timestamp(row.period).strftime('%Y-%m-%d')}) " + f"as stress period for {row.region} " f"({stress_metric} = {row[stress_metric]})" ) ### Include "shoulder periods" before or after each period ### if the storage state of charge is low - # Maintain EUE to be used for determining shoulder periods - stress_metric_for_shoulder_periods = 'EUE' - eue_periods = get_stress_metric_periods( - case=sw.casedir, t=t, iteration=iteration, - hierarchy_level=hierarchy_level, - stress_metric=stress_metric_for_shoulder_periods, - period_agg_method=period_agg_method, - ) - - _eue_sorted = ( - eue_periods - .sort_values(stress_metric_for_shoulder_periods, ascending=False) - .reset_index().set_index('actual_period') - ) - - _high_stress_periods[criterion, f'high_{stress_metric_for_shoulder_periods}'] = ( - _eue_sorted.loc[ - _eue_sorted.r.isin(_failed[criterion].index) - & ~(_eue_sorted.index.isin(stressperiods_this_iteration.actual_period)) - ] - .drop_duplicates(subset=['y','m','d']) - .groupby('r').head(int(sw.GSw_PRM_StressIncrement)) - ) - - _shoulder_periods = { - **_shoulder_periods, - **get_shoulder_periods( + if stress_metric.lower() == 'neue': + shoulder_periods = get_shoulder_periods( sw, criterion, - dfenergy_r, - _high_stress_periods, - stress_metric=stress_metric_for_shoulder_periods, + dfenergy_agg, + high_stress_periods, + stress_metric=stress_metric, ) - } - - stress_criteria['stress_sorted_periods'] = _stress_sorted_periods - stress_criteria['failed'] = _failed - stress_criteria['high_stress_periods'] = _high_stress_periods - stress_criteria['shoulder_periods'] = _shoulder_periods + else: + shoulder_periods = {} - else: - print(f"GSw_PRM_StressThreshold = {criterion} passed") - - return stress_criteria + return { + 'failed': failed, + 'high_stress_periods': high_stress_periods, + 'shoulder_periods': shoulder_periods, + } -def get_stress_periods(sw, t, iteration): - ### Get storage state of charge (SOC) to use in selection of "shoulder" stress periods - dfenergy = reeds.io.read_pras_results( +def get_stress_periods(case, sw, t, iteration): + ### Get values from PRAS + pras_shortfall = get_pras_shortfall(case, t, iteration) + dfeue = pras_shortfall['EUE'] + dflole = pras_shortfall['LOLE'] + ## Storage state of charge (SOC) to use in selection of "shoulder" stress periods + dfenergy_unit = reeds.io.read_pras_results( os.path.join(sw['casedir'], 'handoff', 'PRAS', f"PRAS_{t}i{iteration}-energy.h5") ) timeindex = reeds.timeseries.get_timeindex(sw['resource_adequacy_years']) - dfenergy.index = timeindex - ## Sum by region - dfenergy_r = ( - dfenergy - .rename(columns={c: c.split('|')[1] for c in dfenergy.columns}) + dfenergy_unit.index = timeindex + ## Sum over units + dfenergy = ( + dfenergy_unit + .rename(columns={c: c.split('|')[1] for c in dfenergy_unit.columns}) .groupby(axis=1, level=0).sum() ) @@ -544,50 +514,37 @@ def get_stress_periods(sw, t, iteration): ) ### Check all stress criteria; for regions that fail, add new stress periods - stress_criteria = { - 'stress_sorted_periods': {}, - 'failed': {}, - 'high_stress_periods': {}, - 'shoulder_periods': {}, - } - - # Validation check: Display any GSw_PRM_StressThreshold{metric} - # that is not specified in GSw_PRM_StressThresholdMetrics - stress_metric_switches = sw.GSw_PRM_StressThresholdMetrics.split('/') - - # stress periods column names for writing outputs - stress_metric_units = { - 'NEUE': 'ppm', - 'LOLD': 'event-days/year', - 'LOLH': 'event-hours/year', - 'LOLE': 'events/year', - 'Duration': 'hours', - 'Depth': 'fraction', - } - stress_metrics_col_names = { - m: f'{m}_{stress_metric_units[m]}' for m in stress_metric_switches - } + failed = {} + high_stress_periods = {} + shoulder_periods = {} - for stress_metric in stress_metric_switches: - for criterion in sw[f'GSw_PRM_StressThreshold{stress_metric}'].split('/'): - print(f"Evaluating GSw_PRM_StressThreshold {stress_metric.upper()} with criterion: {criterion}") - stress_criteria = _evaluate_stress_threshold_criterion( - stress_criteria, + stress_metrics = [i.lower() for i in sw.GSw_PRM_StressThresholdMetrics.split('/')] + for stress_metric in stress_metrics: + switch = RA_SWITCHES[stress_metric] + for criterion in sw[switch].split('/'): + ### Aggregate the shortfall and load to this hierarchy level + ## Example: criterion = 'transgrp_1' + hierarchy_level, metric_threshold = criterion.split('_') + rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) + dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum() + dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max() + dfenergy_agg = dfenergy.rename(columns=rmap).groupby(axis=1, level=0).sum() + ## Get the stress periods + dictout = check_threshold_and_choose_periods( + stress_metric, criterion, sw, t, iteration, - dfenergy_r, + dfeue_agg, + dflole_agg, + dfenergy_agg, stressperiods_this_iteration, - stress_metric ) - - stress_sorted_periods = stress_criteria['stress_sorted_periods'] - failed = stress_criteria['failed'] - high_stress_periods = stress_criteria['high_stress_periods'] - shoulder_periods = stress_criteria['shoulder_periods'] - - stress_sorted_periods = pd.concat(stress_sorted_periods, names=['criterion']) + if dictout is not None: + failed.update(dictout['failed']) + high_stress_periods.update(dictout['high_stress_periods']) + shoulder_periods.update(dictout['shoulder_periods']) ### Get lists of stress periods: new (added this iteration) and all if len(failed): @@ -635,9 +592,6 @@ def get_stress_periods(sw, t, iteration): 'Duration': 'duration_hours', 'Depth': 'depth_fraction', } - stress_sorted_periods.round(2).rename(columns=stress_metric_labels).to_csv( - os.path.join(sw.casedir, 'inputs_case', newstresspath, 'stress_metrics_sorted_periods.csv') - ) new_stress_periods.round(2).rename(columns=stress_metric_labels).to_csv( os.path.join(sw.casedir, 'inputs_case', newstresspath, 'new_stress_periods.csv'), index=False, @@ -791,20 +745,20 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): for criterion in failed: if not failed[criterion].name == 'NEUE': continue - # Example: criterion = 'transgrp_10_sum' - (hierarchy_level, stress_value, __) = criterion.split('_') + # Example: criterion = 'transgrp_10' + hierarchy_level, metric_threshold = criterion.split('_') # Recover regions where the PRM criterion failed rmap = reeds.io.get_rmap(sw['casedir'], hierarchy_level=hierarchy_level).reset_index() df = rmap.loc[ rmap[hierarchy_level].isin(failed[criterion].index) ].rename(columns={hierarchy_level:'region'}) df['hierarchy_level'] = hierarchy_level - df['stress_value'] = float(stress_value) + df['metric_threshold'] = float(metric_threshold) _failed_regions.append(df) # For zones that failed multiple criteria, use the most stringent (lowest EUE target) failed_regions = ( pd.concat(_failed_regions) - .sort_values(by=['stress_value']) + .sort_values(by=['metric_threshold']) .drop_duplicates(subset='r', keep='first') ) @@ -840,7 +794,7 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): def main(sw, t, iteration=0, logging=True): """ """ - #%% Write consolidated stress metrics so far + #%% Write consolidated stress metrics ra_metrics = calc_ra_metrics(case=sw.casedir, t=t, iteration=iteration) ra_metrics.round(3).to_csv( os.path.join(sw.casedir, 'outputs', f'ra_metrics_{t}i{iteration}.csv') @@ -858,7 +812,7 @@ def main(sw, t, iteration=0, logging=True): #%% Identify and write new stress periods failed, new_stressperiods_write, combined_periods_write = get_stress_periods( - sw=sw, t=t, iteration=iteration, + case=sw.casedir, sw=sw, t=t, iteration=iteration, ) #%% Stop here if all thresholds pass or if there are no new stress periods From 3d525f254db7d643a1eae770dc295aa8742f4185 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Wed, 24 Jun 2026 10:24:13 -0600 Subject: [PATCH 10/24] stress_periods.py: fix stress/shoulder format --- reeds/resource_adequacy/stress_periods.py | 107 ++++++++++++---------- 1 file changed, 59 insertions(+), 48 deletions(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 1f92fe82..caed4e93 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -1,5 +1,6 @@ #%%### General imports import os +import traceback import pandas as pd import numpy as np import matplotlib.pyplot as plt @@ -63,10 +64,10 @@ def get_pras_shortfall(case, t, iteration=0): return dictout -def plot_stress_diagnostics(sw, t, iteration, high_stress_periods): +def plot_stress_diagnostics(sw, t, iteration, new_stressperiods_write): try: dates = ( - pd.concat(high_stress_periods) + new_stressperiods_write .reset_index().actual_period.map(reeds.timeseries.h2timestamp) .dt.strftime('%Y-%m-%d') .tolist() @@ -87,8 +88,8 @@ def plot_stress_diagnostics(sw, t, iteration, high_stress_periods): os.path.join(sw.casedir, 'outputs', 'figures', 'resource_adequacy', savename) ) plt.close() - except Exception as err: - print(err) + except Exception: + print(traceback.format_exc()) def get_events(ds:pd.Series, threshold:float=0) -> pd.DataFrame: @@ -314,7 +315,7 @@ def get_longest_events( return metric_period.groupby(level=0).sum() -def get_shoulder_periods(sw, criterion, dfenergy_agg, high_stress_periods, stress_metric): +def get_shoulder_periods(sw, criterion, dfenergy_agg, high_stress_periods): ## Stop if not needed if sw.GSw_PRM_StressStorageCutoff.lower() in ['off', '0', 'false']: print( @@ -333,26 +334,27 @@ def get_shoulder_periods(sw, criterion, dfenergy_agg, high_stress_periods, stres timeindex = reeds.timeseries.get_timeindex(sw['resource_adequacy_years']) cutofftype, cutoff = sw.GSw_PRM_StressStorageCutoff.lower().split('_') periodhours = {'day':24, 'wek':24*5, 'year':24}[sw.GSw_HourlyType] + fmt = '%Y-%m-%d' dfheadspace_MWh = dfenergy_agg.max() - dfenergy_agg dfheadspace_frac = dfheadspace_MWh / dfenergy_agg.max() - shoulder_periods = {} + _shoulder_periods = {} for i, row in high_stress_periods.iterrows(): if row.region not in dfheadspace_MWh: continue day = reeds.timeseries.h2timestamp(row.period) - start_headspace_MWh = dfheadspace_MWh.loc[day.strftime('%Y-%m-%d'), row.region].iloc[0] - end_headspace_MWh = dfheadspace_MWh.loc[day.strftime('%Y-%m-%d'), row.region].iloc[-1] + start_headspace_MWh = dfheadspace_MWh.loc[day.strftime(fmt), row.region].iloc[0] + end_headspace_MWh = dfheadspace_MWh.loc[day.strftime(fmt), row.region].iloc[-1] - start_headspace_frac = dfheadspace_frac.loc[day.strftime('%Y-%m-%d'), row.region].iloc[0] - end_headspace_frac = dfheadspace_frac.loc[day.strftime('%Y-%m-%d'), row.region].iloc[-1] + start_headspace_frac = dfheadspace_frac.loc[day.strftime(fmt), row.region].iloc[0] + end_headspace_frac = dfheadspace_frac.loc[day.strftime(fmt), row.region].iloc[-1] - day_eue = high_stress_periods.loc[i, stress_metric] + day_eue = high_stress_periods.loc[i, 'value'] day_index = np.where( - timeindex == dfenergy_agg.loc[day.strftime('%Y-%m-%d')].iloc[0].name + timeindex == dfenergy_agg.loc[day.strftime(fmt)].iloc[0].name )[0][0] day_before = timeindex[day_index - periodhours] @@ -363,23 +365,33 @@ def get_shoulder_periods(sw, criterion, dfenergy_agg, high_stress_periods, stres or ((cutofftype[:3] == 'cap') and (end_headspace_frac >= float(cutoff))) or (cutofftype[:3] == 'abs') ): - shoulder_periods[criterion, f'after_{row.name}'] = pd.Series({ - 'actual_period':day_after.strftime('y%Yd%j'), + _shoulder_periods[f'after_{row.period}'] = pd.Series({ 'region':row.region, - }).to_frame().T.set_index('actual_period') - print(f"Added {day_after} as shoulder stress period after {day}") + 'period':day_after.strftime('y%Yd%j'), + }) + print( + f"Added {day_after.strftime(fmt)} as shoulder stress period " + f"after {day.strftime(fmt)}" + ) if ( ((cutofftype == 'eue') and (start_headspace_MWh / day_eue >= float(cutoff))) or ((cutofftype[:3] == 'cap') and (start_headspace_frac >= float(cutoff))) or (cutofftype[:3] == 'abs') ): - shoulder_periods[criterion, f'before_{row.name}'] = pd.Series({ - 'actual_period':day_before.strftime('y%Yd%j'), + _shoulder_periods[f'before_{row.period}'] = pd.Series({ 'region':row.region, - }).to_frame().T.set_index('actual_period') - print(f"Added {day_before} as shoulder stress period before {day}") + 'period':day_before.strftime('y%Yd%j'), + }) + print( + f"Added {day_before.strftime(fmt)} as shoulder stress period " + f"before {day.strftime(fmt)}" + ) + shoulder_periods = ( + pd.concat(_shoulder_periods).unstack(level=1).reset_index() + .rename(columns={'index':'value'})[['region','period','value']] + ) return shoulder_periods @@ -423,7 +435,6 @@ def check_threshold_and_choose_periods( for region in failed.index } case 'duration': - ## TODO: Keep whole events (including when they span days) metric_periods = { region: get_longest_events( sw=sw, t=t, iteration=iteration, @@ -458,7 +469,7 @@ def check_threshold_and_choose_periods( } high_stress_periods = ( pd.concat(metric_periods, names=['region','period']) - .rename(stress_metric) + .rename('value') .reset_index() ) for i, row in high_stress_periods.iterrows(): @@ -466,7 +477,7 @@ def check_threshold_and_choose_periods( f"Added {row.period} " f"({reeds.timeseries.h2timestamp(row.period).strftime('%Y-%m-%d')}) " f"as stress period for {row.region} " - f"({stress_metric} = {row[stress_metric]})" + f"({stress_metric} = {np.around(row.value, 3)})" ) ### Include "shoulder periods" before or after each period @@ -477,10 +488,9 @@ def check_threshold_and_choose_periods( criterion, dfenergy_agg, high_stress_periods, - stress_metric=stress_metric, ) else: - shoulder_periods = {} + shoulder_periods = pd.DataFrame() return { 'failed': failed, @@ -514,9 +524,9 @@ def get_stress_periods(case, sw, t, iteration): ) ### Check all stress criteria; for regions that fail, add new stress periods - failed = {} - high_stress_periods = {} - shoulder_periods = {} + _failed = {} + _high_stress_periods = {} + _shoulder_periods = {} stress_metrics = [i.lower() for i in sw.GSw_PRM_StressThresholdMetrics.split('/')] for stress_metric in stress_metrics: @@ -542,30 +552,39 @@ def get_stress_periods(case, sw, t, iteration): stressperiods_this_iteration, ) if dictout is not None: - failed.update(dictout['failed']) - high_stress_periods.update(dictout['high_stress_periods']) - shoulder_periods.update(dictout['shoulder_periods']) + _failed[stress_metric, criterion] = dictout['failed'] + _high_stress_periods[stress_metric, criterion] = dictout['high_stress_periods'] + _shoulder_periods[stress_metric, criterion] = dictout['shoulder_periods'] + + failed = pd.concat(_failed) + high_stress_periods = pd.concat(_high_stress_periods) + shoulder_periods = pd.concat(_shoulder_periods) ### Get lists of stress periods: new (added this iteration) and all if len(failed): new_stress_periods = pd.concat( - {**high_stress_periods, **shoulder_periods}, names=['criterion','periodtype'], - ).reset_index().drop_duplicates(subset='actual_period', keep='first') + {'stress':high_stress_periods, 'shoulder':shoulder_periods}, + names=['periodtype','metric','criterion','num'], + ).reset_index() + print('All identified stress periods:') + print(new_stress_periods) + new_stress_periods = new_stress_periods.drop_duplicates('period') else: return failed, None, None ## Reproduce the format of inputs_case/stress_period_szn.csv p = 'w' if sw.GSw_HourlyType == 'wek' else 'd' new_stressperiods_write = pd.DataFrame({ - 'rep_period': new_stress_periods.actual_period, - 'year': new_stress_periods.actual_period.map( + 'rep_period': new_stress_periods.period, + 'year': new_stress_periods.period.map( lambda x: int(x.strip('sy').split(p)[0])), - 'yperiod': new_stress_periods.actual_period.map( + 'yperiod': new_stress_periods.period.map( lambda x: int(x.strip('sy').split(p)[1])), - 'actual_period': new_stress_periods.actual_period, + 'actual_period': new_stress_periods.period, }) - ### Add new stress periods to the stress periods used for this year/iteration, then write + ### Add new stress periods to the stress periods used for this year/iteration, + ### drop duplicates, then write newstresspath = f'stress{t}i{iteration+1}' os.makedirs(os.path.join(sw['casedir'], 'inputs_case', newstresspath), exist_ok=True) outpath = os.path.join(sw['casedir'], 'inputs_case', newstresspath, 'period_szn.csv') @@ -584,19 +603,11 @@ def get_stress_periods(case, sw, t, iteration): combined_periods_write.to_csv(outpath, index=False) ### Tables and plots for debugging - stress_metric_labels = { - 'NEUE': 'neue_ppm', - 'LOLD': 'lold_event-days/year', - 'LOLH': 'lolh_event-hours/year', - 'LOLE': 'lole_events/year', - 'Duration': 'duration_hours', - 'Depth': 'depth_fraction', - } - new_stress_periods.round(2).rename(columns=stress_metric_labels).to_csv( + new_stress_periods.to_csv( os.path.join(sw.casedir, 'inputs_case', newstresspath, 'new_stress_periods.csv'), index=False, ) - plot_stress_diagnostics(sw, t, iteration, high_stress_periods) + plot_stress_diagnostics(sw, t, iteration, new_stressperiods_write) return failed, new_stressperiods_write, combined_periods_write From fc09c9aeb2ffb255322af2b3533b293fc2fa3c8c Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Wed, 24 Jun 2026 10:33:19 -0600 Subject: [PATCH 11/24] move plot_stress_diagnostics() to diagnostic_plots.map_outagerate_new_stressperiods() --- reeds/resource_adequacy/diagnostic_plots.py | 38 +++++++++++++++++++++ reeds/resource_adequacy/stress_periods.py | 35 ++----------------- 2 files changed, 40 insertions(+), 33 deletions(-) diff --git a/reeds/resource_adequacy/diagnostic_plots.py b/reeds/resource_adequacy/diagnostic_plots.py index 36145664..79685f7f 100644 --- a/reeds/resource_adequacy/diagnostic_plots.py +++ b/reeds/resource_adequacy/diagnostic_plots.py @@ -85,6 +85,10 @@ def get_inputs(sw): resources['tech'] = reeds.reedsplots.simplify_techs(resources.i, display_level = 'diagnostics') resources['rb'] = resources.r + new_stress_periods = pd.read_csv( + Path(sw.casedir, 'inputs_case', f'stress{sw.t}i{sw.iteration+1}', 'new_stress_periods.csv') + ) + ##### Hourly dispatch by month ### Load and aggregate the VRE generation profiles by tech group try: @@ -242,6 +246,7 @@ def get_inputs(sw): dfs['tech_style'] = tech_style dfs['vre_gen_usa'] = vre_gen_usa dfs['vre_gen'] = vre_gen + dfs['new_stress_periods'] = new_stress_periods return dfs @@ -1023,6 +1028,34 @@ def map_pras_failure_rate(sw, dfs, aggfunc='mean', repair=False): plt.close() +def map_outagerate_new_stressperiods(sw, dfs): + new_stress_periods = dfs['new_stress_periods'] + dates = ( + new_stress_periods + .period.map(reeds.timeseries.h2timestamp) + .dt.strftime('%Y-%m-%d') + .tolist() + ) + vmax = {'forced': 40, 'scheduled': 25, 'both': 50} + aggfunc = 'max' + for outage_type in vmax: + savename = f'map-outage_{outage_type}_{aggfunc}-{t}i{iteration}.png' + plt.close() + f, ax, _ = reeds.reedsplots.map_outage_days( + sw.casedir, + dates=dates, + outage_type=outage_type, + aggfunc=aggfunc, + vmax=vmax[outage_type], + ) + ## Save it + if savefig: + plt.savefig(os.path.join(sw['savepath'],savename)) + if interactive: + plt.show() + plt.close() + + def plot_cc_mar(sw, dfs): """ Marginal capacity credit @@ -1260,6 +1293,11 @@ def main(sw, debug=False): except Exception: print('map_dropped_load() failed:', traceback.format_exc()) + try: + map_outagerate_new_stressperiods(sw, dfs) + except Exception: + print('map_outagerate_new_stressperiods() failed:', traceback.format_exc()) + if int(sw['GSw_PRM_CapCredit']): try: plot_cc_mar(sw, dfs) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index caed4e93..e9a702e2 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -1,9 +1,7 @@ #%%### General imports import os -import traceback -import pandas as pd import numpy as np -import matplotlib.pyplot as plt +import pandas as pd from pathlib import Path from typing import Literal @@ -64,34 +62,6 @@ def get_pras_shortfall(case, t, iteration=0): return dictout -def plot_stress_diagnostics(sw, t, iteration, new_stressperiods_write): - try: - dates = ( - new_stressperiods_write - .reset_index().actual_period.map(reeds.timeseries.h2timestamp) - .dt.strftime('%Y-%m-%d') - .tolist() - ) - vmax = {'forced': 40, 'scheduled': 25, 'both': 50} - aggfunc = 'max' - for outage_type in vmax: - savename = f'map-outage_{outage_type}_{aggfunc}-{t}i{iteration}.png' - plt.close() - f, ax, _ = reeds.reedsplots.map_outage_days( - sw.casedir, - dates=dates, - outage_type=outage_type, - aggfunc=aggfunc, - vmax=vmax[outage_type], - ) - plt.savefig( - os.path.join(sw.casedir, 'outputs', 'figures', 'resource_adequacy', savename) - ) - plt.close() - except Exception: - print(traceback.format_exc()) - - def get_events(ds:pd.Series, threshold:float=0) -> pd.DataFrame: """Return a dataframe of events with max and duration""" starts = ( @@ -602,12 +572,11 @@ def get_stress_periods(case, sw, t, iteration): else: combined_periods_write.to_csv(outpath, index=False) - ### Tables and plots for debugging + ### Tables for debugging new_stress_periods.to_csv( os.path.join(sw.casedir, 'inputs_case', newstresspath, 'new_stress_periods.csv'), index=False, ) - plot_stress_diagnostics(sw, t, iteration, new_stressperiods_write) return failed, new_stressperiods_write, combined_periods_write From ff29cc55d5806820a499f69b3e75bc0736b7dc16 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Wed, 24 Jun 2026 10:42:31 -0600 Subject: [PATCH 12/24] stress_periods.py: fix handling of no new stress periods --- reeds/resource_adequacy/stress_periods.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index e9a702e2..a5897416 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -526,12 +526,11 @@ def get_stress_periods(case, sw, t, iteration): _high_stress_periods[stress_metric, criterion] = dictout['high_stress_periods'] _shoulder_periods[stress_metric, criterion] = dictout['shoulder_periods'] - failed = pd.concat(_failed) - high_stress_periods = pd.concat(_high_stress_periods) - shoulder_periods = pd.concat(_shoulder_periods) - ### Get lists of stress periods: new (added this iteration) and all - if len(failed): + if len(_failed): + failed = pd.concat(_failed) + high_stress_periods = pd.concat(_high_stress_periods) + shoulder_periods = pd.concat(_shoulder_periods) new_stress_periods = pd.concat( {'stress':high_stress_periods, 'shoulder':shoulder_periods}, names=['periodtype','metric','criterion','num'], @@ -540,7 +539,7 @@ def get_stress_periods(case, sw, t, iteration): print(new_stress_periods) new_stress_periods = new_stress_periods.drop_duplicates('period') else: - return failed, None, None + return {}, {}, {} ## Reproduce the format of inputs_case/stress_period_szn.csv p = 'w' if sw.GSw_HourlyType == 'wek' else 'd' From 5388a3d637d6971a6ce83c02014b43c0a660db79 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Wed, 24 Jun 2026 12:34:00 -0600 Subject: [PATCH 13/24] stress_periods.py: streamline logs; diagnostic_plots.py: fix new_stress_periods plot --- reeds/resource_adequacy/diagnostic_plots.py | 12 ++++++------ reeds/resource_adequacy/ra_calcs.py | 14 ++++++-------- reeds/resource_adequacy/stress_periods.py | 16 ++++++---------- 3 files changed, 18 insertions(+), 24 deletions(-) diff --git a/reeds/resource_adequacy/diagnostic_plots.py b/reeds/resource_adequacy/diagnostic_plots.py index 79685f7f..1b4e6b01 100644 --- a/reeds/resource_adequacy/diagnostic_plots.py +++ b/reeds/resource_adequacy/diagnostic_plots.py @@ -85,9 +85,11 @@ def get_inputs(sw): resources['tech'] = reeds.reedsplots.simplify_techs(resources.i, display_level = 'diagnostics') resources['rb'] = resources.r - new_stress_periods = pd.read_csv( - Path(sw.casedir, 'inputs_case', f'stress{sw.t}i{sw.iteration+1}', 'new_stress_periods.csv') - ) + fpath = Path(sw.casedir, 'inputs_case', f'stress{sw.t}i{sw.iteration+1}', 'new_stress_periods.csv') + if fpath.is_file(): + new_stress_periods = pd.read_csv(fpath) + else: + new_stress_periods = pd.DataFrame(columns=['period']) ##### Hourly dispatch by month ### Load and aggregate the VRE generation profiles by tech group @@ -549,9 +551,7 @@ def plot_reeds_pras_capacity(sw, dfs): ### Get coordinates zones = dfs['hierarchy'].index - ncols = int(np.around(np.sqrt(len(zones)) * 1.618, 0)) - nrows = len(zones) // ncols + int(bool(len(zones) % ncols)) - coords = dict(zip(zones, [(row,col) for row in range(nrows) for col in range(ncols)])) + nrows, ncols, coords = reeds.plots.get_coordinates(zones) ### Plot the capacities plt.close() diff --git a/reeds/resource_adequacy/ra_calcs.py b/reeds/resource_adequacy/ra_calcs.py index 214f26ab..cd325d75 100644 --- a/reeds/resource_adequacy/ra_calcs.py +++ b/reeds/resource_adequacy/ra_calcs.py @@ -99,11 +99,11 @@ def run_pras( def main(t, tnext, casedir, iteration=0): # #%% To debug, uncomment these lines and update the run path - # t = 2020 - # tnext = 2023 + # t = 2030 + # tnext = 2035 # reeds_path = reeds.io.reeds_path # casedir = os.path.join( - # reeds_path,'runs','v20260417_reorgM1_Pacific') + # reeds_path,'runs','v20260624_raM0_MultiMetricRA') # iteration = 0 # assert tnext >= t # os.chdir(casedir) @@ -191,11 +191,9 @@ def main(t, tnext, casedir, iteration=0): # #%% Uncomment to run diagnostic_plots # ### (typically run from call_{}.sh script for parallelization) - # try: - # reeds.resource_adequacy.diagnostic_plots.main(sw) - # except Exception as err: - # print('diagnostic_plots.py failed with the following exception:') - # print(err) + # from reeds.resource_adequacy import diagnostic_plots + # sw['iteration'] = iteration + # diagnostic_plots.main(sw) #%% Procedure diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index a5897416..b1774876 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -391,10 +391,14 @@ def check_threshold_and_choose_periods( this_test = ra_metrics[hierarchy_level][SWITCH_METRIC[stress_metric]] failed = this_test.loc[this_test > metric_threshold] if not len(failed): - print(f"{RA_SWITCHES[stress_metric]} = {criterion} passed") + print(f"{RA_SWITCHES[stress_metric]} = {criterion} passed:") + for i, val in this_test.items(): + print(f'{i}: {val} {stress_metric}') else: print(f"{RA_SWITCHES[stress_metric]} = {criterion} failed for:") - print(failed) + for i, val in this_test.items(): + print(f'{i}: {val} {stress_metric}') + ## Get new stress periods since the metric failed match stress_metric: case 'depth': metric_periods = { @@ -442,14 +446,6 @@ def check_threshold_and_choose_periods( .rename('value') .reset_index() ) - for i, row in high_stress_periods.iterrows(): - print( - f"Added {row.period} " - f"({reeds.timeseries.h2timestamp(row.period).strftime('%Y-%m-%d')}) " - f"as stress period for {row.region} " - f"({stress_metric} = {np.around(row.value, 3)})" - ) - ### Include "shoulder periods" before or after each period ### if the storage state of charge is low if stress_metric.lower() == 'neue': From 9bbb98e8a842bcd19dffaa56b3c850d73f3d0374 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Wed, 24 Jun 2026 14:39:22 -0600 Subject: [PATCH 14/24] stress_periods.py: for shoulder periods on first day, loop around to end of timeseries --- reeds/resource_adequacy/stress_periods.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index b1774876..9320cca2 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -327,7 +327,8 @@ def get_shoulder_periods(sw, criterion, dfenergy_agg, high_stress_periods): timeindex == dfenergy_agg.loc[day.strftime(fmt)].iloc[0].name )[0][0] - day_before = timeindex[day_index - periodhours] + ## Loop around with % len(timeindex) if on the first or last day + day_before = timeindex[(day_index - periodhours) % len(timeindex)] day_after = timeindex[(day_index + periodhours) % len(timeindex)] if ( From b718a2209ec80f3d1131852f16ba5b3c2dd8af85 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Wed, 24 Jun 2026 14:39:41 -0600 Subject: [PATCH 15/24] stress_periods.py: fix empty shoulder periods --- cases.csv | 2 +- reeds/resource_adequacy/stress_periods.py | 14 ++++++++------ 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/cases.csv b/cases.csv index 4cc03fca..a7163e24 100644 --- a/cases.csv +++ b/cases.csv @@ -269,7 +269,7 @@ GSw_PRM_StressModel,Model used to identify stress periods: pras or a string star GSw_PRM_StressOutages,Whether to apply the availability factor (forced + scheduled outages) during stress periods,0; 1,1, GSw_PRM_StressSeedLoadLevel,Region hierarchy level at which to include peak coincident load days as seeded stress periods (or False to ignore peaks),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,transgrp, GSw_PRM_StressSeedMinRElevel,Region hierarchy level at which to include minimum wind and solar capacity factor days as seeded stress periods (or False to ignore min-wind/solar CF days),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,interconnect, -GSw_PRM_StressStorageCutoff,"How to select ""shoulder"" stress periods giving storage time to recharge before/after high-unserved-energy periods. Two-part switch separated by _. The first argument is ""EUE"" or ""capacity"" or ""absolute"". If ""EUE"" the second argument specifies a PRAS storage headspace / EUE threshold [fraction]; if ""cap"" it specifies a headspace / storage capacity threshold [fraction]; if ""abs"" it specifies absolute number of periods before/after [integer]. Turned off if set to ""off"".",N/A,EUE_0.1, +GSw_PRM_StressStorageCutoff,How to select shoulder stress periods giving storage time to recharge before/after high-unserved-energy periods. Two-part switch separated by _. The first argument is 'EUE' or 'capacity' or 'absolute'. If 'EUE' the second argument specifies a PRAS storage headspace / EUE threshold [fraction]; if 'cap' it specifies a headspace / storage capacity threshold [fraction]; if 'abs' it specifies absolute number of periods before/after [integer]. Turned off if set to 'off'.,N/A,EUE_0.1, 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, GSw_PRM_StressThresholdDepth,Outage depth threshold [MW_EUE/MW_peak_load] (fraction); formulated as HierarchyLevel_Depth where HierarchyLevel is a column in hierarchy.csv; Depth is the max outage magnitude in [MW EUE] / [MW peak demand],N/A,transgrp_0.1, GSw_PRM_StressThresholdDuration,Outage duration threshold [hours]; formulated as HierarchyLevel_Duration where HierarchyLevel is a column in hierarchy.csv; Duration is the max outage duration in hours,N/A,transgrp_12, diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 9320cca2..1c6c8728 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -358,11 +358,13 @@ def get_shoulder_periods(sw, criterion, dfenergy_agg, high_stress_periods): f"Added {day_before.strftime(fmt)} as shoulder stress period " f"before {day.strftime(fmt)}" ) - - shoulder_periods = ( - pd.concat(_shoulder_periods).unstack(level=1).reset_index() - .rename(columns={'index':'value'})[['region','period','value']] - ) + if len(_shoulder_periods): + shoulder_periods = ( + pd.concat(_shoulder_periods).unstack(level=1).reset_index() + .rename(columns={'index':'value'})[['region','period','value']] + ) + else: + shoulder_periods = pd.DataFrame(columns=['region','period','value']) return shoulder_periods @@ -397,7 +399,7 @@ def check_threshold_and_choose_periods( print(f'{i}: {val} {stress_metric}') else: print(f"{RA_SWITCHES[stress_metric]} = {criterion} failed for:") - for i, val in this_test.items(): + for i, val in failed.items(): print(f'{i}: {val} {stress_metric}') ## Get new stress periods since the metric failed match stress_metric: From 50f0dc174dae3221d70ca47afda57b13c6130b45 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Wed, 24 Jun 2026 15:50:42 -0600 Subject: [PATCH 16/24] fix some plots --- postprocessing/single_case_plots.py | 13 ++-- reeds/reedsplots.py | 101 ++++++++++++++++------------ 2 files changed, 64 insertions(+), 50 deletions(-) diff --git a/postprocessing/single_case_plots.py b/postprocessing/single_case_plots.py index 8ab4306d..b8a5e6da 100644 --- a/postprocessing/single_case_plots.py +++ b/postprocessing/single_case_plots.py @@ -71,7 +71,7 @@ year = args.year # #%% Inputs for testing -# case = os.path.join(reeds_path,'runs','v20251111_15M0_Pacific') +# case = os.path.join(reeds_path,'runs','v20260624_raM1_MultiMetricRA') # year = 0 # interactive = True # write = False @@ -580,7 +580,7 @@ try: plt.close() levels = ['country', 'interconnect', 'transreg', 'transgrp'] - f, ax, _ = reedsplots.plot_neue_bylevel(case=case, levels=levels) + f, ax, _ = reedsplots.plot_ra_metrics_bylevel(case=case, levels=levels) savename = f"plot_stressperiod_neue-{','.join(levels)}.png" if write: plt.savefig(os.path.join(savepath, savename)) @@ -717,13 +717,10 @@ print(traceback.format_exc()) try: - _first_metric = sw['GSw_PRM_StressThresholdMetrics'].split('/')[0].upper() - _parts = sw[f'GSw_PRM_StressThreshold{_first_metric}'].split('_') - level, threshold, metric = _parts[0], _parts[1], _parts[2] + level, threshold = sw['GSw_PRM_StressThresholdNEUE'].split('/')[0].split('_') plt.close() - f,ax = reedsplots.plot_stressperiod_evolution( - case=case, level=level, metric=metric) - savename = f'plot_stressperiod_evolution-{metric}-{level}.png' + f,ax = reedsplots.plot_stressperiod_evolution(case=case, level=level) + savename = f'plot_stressperiod_evolution-neue-{level}.png' if write: plt.savefig(os.path.join(savepath, savename)) if interactive: diff --git a/reeds/reedsplots.py b/reeds/reedsplots.py index f6ecb27e..6d0edd06 100644 --- a/reeds/reedsplots.py +++ b/reeds/reedsplots.py @@ -4121,22 +4121,19 @@ def plot_stressperiod_days(case, repcolor='k', sharey=False, figsize=(10,5)): def plot_stressperiod_evolution( - case, level=None, metric=None, threshold=None, + case, level=None, threshold=None, figsize=None, scale_widths=False, ): """Plot NEUE by year and stress period iteration""" ### Parse inputs sw = reeds.io.get_switches(case) - _first_metric = sw['GSw_PRM_StressThresholdMetrics'].split('/')[0].upper() - _parts = sw[f'GSw_PRM_StressThreshold{_first_metric}'].split('_') - _level, _threshold, _metric = _parts[0], _parts[1], _parts[2] + _level, _threshold = sw['GSw_PRM_StressThresholdNEUE'].split('/')[0].split('_') level = _level if level is None else level threshold = float(_threshold) if threshold is None else threshold - metric = _metric if metric is None else metric ### Load NEUE results - infiles = sorted(glob(os.path.join(case,'outputs','neue_*.csv'))) + infiles = sorted(glob(os.path.join(case,'outputs','ra_metrics_*.csv'))) dictin_neue = { - tuple([int(x) for x in os.path.basename(f)[len('neue_'):-len('.csv')].split('i')]): + tuple([int(x) for x in os.path.basename(f)[len('ra_metrics_'):-len('.csv')].split('i')]): pd.read_csv(f, index_col=['level','metric','region']) for f in infiles } @@ -4144,8 +4141,8 @@ def plot_stressperiod_evolution( dfplot = ( pd.concat(dictin_neue, names=['year','iteration']) .xs(level,0,'level') - .xs(metric,0,'metric') - .NEUE.unstack('region') + .xs('neue_ppm',0,'metric') + .squeeze(1).unstack('region') ) ### Load stress periods for labels dfstress = get_stressperiods(case) @@ -4206,18 +4203,26 @@ def plot_stressperiod_evolution( return f,ax -def plot_neue_bylevel( - case, tmin=2023, - levels=['country','interconnect','transreg','transgrp'], - metrics=['sum','max'], - onlydata=False, - ): +def plot_ra_metrics_bylevel( + case, tmin=2023, + levels=['country','interconnect','transreg','transgrp'], + metrics=[ + 'neue', + 'lolh', + 'lole', + 'lold', + 'duration', + 'depth', + ], + onlydata=False, +): """Plot regional NEUE over time""" + from reeds.resource_adequacy import stress_periods ### Get final iterations year2iteration = ( pd.DataFrame([ - os.path.basename(i).strip('neue_.csv').split('i') - for i in sorted(glob(os.path.join(case, 'outputs', 'neue_*.csv'))) + os.path.basename(i).strip('ra_metrics.csv').split('i') + for i in sorted(glob(os.path.join(case, 'outputs', 'ra_metrics_*.csv'))) ], columns=['year','iteration']).astype(int) .drop_duplicates(subset='year', keep='last') .set_index('year').iteration @@ -4226,58 +4231,70 @@ def plot_neue_bylevel( ### Get NEUE sw = reeds.io.get_switches(case) sw['casedir'] = case - dictin_neue = {} + ra_metrics = {} for t, iteration in year2iteration.items(): try: - dictin_neue[t] = ( - reeds.io.read_output(case, f'neue_{t}i{iteration}.csv') + ra_metrics[t] = ( + reeds.io.read_output(case, f'ra_metrics_{t}i{iteration}.csv') .set_index(['level','metric','region']).squeeze(1) ) except FileNotFoundError: - dictin_neue[t] = ( - reeds.io.read_output(case, f'neue_{t}i{iteration-1}.csv') + ra_metrics[t] = ( + reeds.io.read_output(case, f'ra_metrics_{t}i{iteration-1}.csv') .set_index(['level','metric','region']).squeeze(1) ) - dfin_neue = pd.concat(dictin_neue, axis=0, names=['year']).unstack('year') - dfin_neue = dfin_neue[[c for c in dfin_neue if int(c) >= 2025]].copy() + dfin_ra = pd.concat(ra_metrics, axis=0, names=['year']).unstack('year') + dfin_ra = dfin_ra[[c for c in dfin_ra if int(c) >= 2025]].copy() if onlydata: - return dfin_neue + return dfin_ra ### Plot settings ncols = len(levels) nrows = len(metrics) colors = { level: plots.rainbowmapper( - dfin_neue.xs(level,0,'level').reset_index().region.unique()) + dfin_ra.xs(level,0,'level').reset_index().region.unique()) for level in levels } - norm = {'sum':1, 'max':1e-4} - ylabel = {'sum': 'Sum of NEUE [ppm]', 'max':'Max NEUE [%]'} - thresholds = { - sw[f'GSw_PRM_StressThreshold{m.upper()}'].split('_')[0]: - float(sw[f'GSw_PRM_StressThreshold{m.upper()}'].split('_')[1]) - for m in sw.GSw_PRM_StressThresholdMetrics.split('/') + ylabel = { + 'neue': 'NEUE\n[ppm]', + 'lolh': 'LOLH\n[event-hours/year]', + 'lole': 'LOLE\n[events/year]', + 'lold': 'LOLD\n[event-days/year]', + 'duration': 'Duration\n[hours]', + 'depth': 'Depth\n[% of peak]', } + scale = {'depth':100} ### Plot it plt.close() f,ax = plt.subplots( - nrows, ncols, figsize=(2*ncols, 3*nrows), + nrows, ncols, figsize=(2*ncols, 1.5*nrows), sharex=True, sharey='row', ) for row, metric in enumerate(metrics): + switch = stress_periods.RA_SWITCHES[metric.lower()] + switch_metric = stress_periods.SWITCH_METRIC[metric.lower()] + thresholds = { + i.split('_')[0]: float(i.split('_')[1]) + for i in sw[switch].split('/') + } for col, level in enumerate(levels): for region, c in colors[level].items(): ax[row,col].plot( - dfin_neue.columns, - dfin_neue.loc[(level,metric,region)].values * norm[metric], + dfin_ra.columns, + dfin_ra.loc[(level,switch_metric,region)].values * scale.get(metric,1), c=c, label=region, marker='o', markersize=2, ) + if ( + (level in thresholds) + and (metric in sw.GSw_PRM_StressThresholdMetrics.lower().split('/')) + and (not int(sw.GSw_PRM_CapCredit)) + ): + ax[row,col].axhline(thresholds[level] * scale.get(metric,1), c='k', ls=':', lw=0.75) ## Formatting - ax[row,0].set_ylabel(ylabel[metric]) + ax[row,0].set_ylabel(ylabel.get(metric, metric), ha='right', va='center', rotation=0) ## Formatting for col, level in enumerate(levels): ax[0,col].set_title(level) - if (level in thresholds) and (not int(sw.GSw_PRM_CapCredit)): - ax[0,col].axhline(thresholds[level], c='k', ls=':', lw=0.75) leg = ax[0,col].legend( loc='upper left', frameon=False, fontsize=8, handletextpad=0.3, handlelength=0.7, columnspacing=0.5, labelspacing=0.2, @@ -4292,7 +4309,7 @@ def plot_neue_bylevel( ax[0,0].set_ylim(0, min(ax[0,0].get_ylim()[1], 1000)) plots.despine(ax) - return f, ax, dfin_neue + return f, ax, dfin_ra def map_neue( @@ -4316,7 +4333,7 @@ def map_neue( case=case, year=year, samples=samples) else: _iteration = iteration - neue = reeds.io.read_output(case, f'neue_{year}i{_iteration}.csv') + neue = reeds.io.read_output(case, f'ra_metrics_{year}i{_iteration}.csv') neue = neue.loc[neue.metric==metric].set_index(['level','region']).NEUE sw = reeds.io.get_switches(case) neue_threshold = float(sw.GSw_PRM_StressThresholdNEUE.split('_')[1]) @@ -6632,8 +6649,8 @@ def map_prm(case, tmin=2023, cmap=cmocean.cm.rain, scale=3, fontsize=7, vmax=Non ### Get final iterations year2iteration = ( pd.DataFrame([ - os.path.basename(i).strip('neue_.csv').split('i') - for i in sorted(glob(os.path.join(case, 'outputs', 'neue_*.csv'))) + os.path.basename(i).strip('ra_metrics.csv').split('i') + for i in sorted(glob(os.path.join(case, 'outputs', 'ra_metrics_*.csv'))) ], columns=['year','iteration']).astype(int) .drop_duplicates(subset='year', keep='last') .set_index('year').iteration From b125b1d1254ab6c6073d25b18f0c9eb06d254ca6 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Thu, 25 Jun 2026 11:12:04 -0600 Subject: [PATCH 17/24] docs: adapt ra-flowcharts.png for multimetric RA --- docs/source/figs/docs/ra-flowcharts.png | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/figs/docs/ra-flowcharts.png b/docs/source/figs/docs/ra-flowcharts.png index 52096314..c7b95f0e 100644 --- a/docs/source/figs/docs/ra-flowcharts.png +++ b/docs/source/figs/docs/ra-flowcharts.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:0258efe954e94b093c95fc84487b38b803fe55303aec53af626f3a14500aa5bb -size 306828 +oid sha256:f897397797604afa5ab64d2b1f9af88bd6455f3f24dbbc3d198370bb2234669f +size 288248 From 0a936a4072186bf76220146383b2ca5b1b386308 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Thu, 25 Jun 2026 15:44:31 -0600 Subject: [PATCH 18/24] add profile to eue_events output --- postprocessing/compare_cases.py | 6 +++--- postprocessing/single_case_plots.py | 2 +- reeds/resource_adequacy/stress_periods.py | 5 ++++- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/postprocessing/compare_cases.py b/postprocessing/compare_cases.py index 56c91d8d..a005ce65 100644 --- a/postprocessing/compare_cases.py +++ b/postprocessing/compare_cases.py @@ -472,12 +472,12 @@ def plot_bars_abs_stacked( dictin_neue = {} dictin_neue_all = {} for case in tqdm(cases, desc='NEUE'): - infiles = sorted(glob(os.path.join(cases[case],'outputs','neue_*.csv'))) + infiles = sorted(glob(os.path.join(cases[case],'outputs','ra_metrics_*.csv'))) if not len(infiles): continue df = {} for f in infiles: - y, i = [int(s) for s in os.path.basename(f).strip('neue_.csv').split('i')] + y, i = [int(s) for s in os.path.basename(f).strip('ra_metrics_.csv').split('i')] df[y,i] = pd.read_csv(f, index_col=['level', 'metric', 'region']).squeeze(1) dictin_neue_all[case] = pd.concat(df, names=('t', 'iteration')) indices = ['t', 'level', 'metric', 'region'] @@ -2517,7 +2517,7 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): ] for figname, width, height in [ (f'map_gencap_transcap-{lastyear}', None, SLIDE_HEIGHT), - (f'plot_stressperiod_evolution-sum-{level}', SLIDE_WIDTH, None), + (f'plot_stressperiod_evolution-neue-{level}', SLIDE_WIDTH, None), (f'plot_dispatch-yearbymonth-1-{lastyear}', SLIDE_WIDTH, None), ## Include both versions for backwards compatibility (f'plot_dispatch-yearbymonth-1-{lastyear}-w{weatheryear}', SLIDE_WIDTH, None), diff --git a/postprocessing/single_case_plots.py b/postprocessing/single_case_plots.py index b8a5e6da..d5664336 100644 --- a/postprocessing/single_case_plots.py +++ b/postprocessing/single_case_plots.py @@ -581,7 +581,7 @@ plt.close() levels = ['country', 'interconnect', 'transreg', 'transgrp'] f, ax, _ = reedsplots.plot_ra_metrics_bylevel(case=case, levels=levels) - savename = f"plot_stressperiod_neue-{','.join(levels)}.png" + savename = f"plot_ra_metrics-{','.join(levels)}.png" if write: plt.savefig(os.path.join(savepath, savename)) if interactive: diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 1c6c8728..929e9ec8 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -94,11 +94,12 @@ def get_events(ds:pd.Series, threshold:float=0) -> pd.DataFrame: 'timesteps': len(event), 'max': event.max(), 'sum': event.sum(), + 'profile': event.round(1).values, }) if len(events): dfout = pd.DataFrame(events) else: - dfout = pd.DataFrame(columns=['start','end','timesteps','max','sum']) + dfout = pd.DataFrame(columns=['start','end','timesteps','max','sum','profile']) return dfout @@ -780,6 +781,8 @@ def main(sw, t, iteration=0, logging=True): #%% Write EUE events eue_events = get_eue_events(case=sw.casedir, t=t, iteration=iteration) + ## Store the profile as a |-delimited string + eue_events.profile = eue_events.profile.map(lambda x: '|'.join(str(i) for i in x)) eue_events.round(3).to_csv( os.path.join(sw.casedir, 'outputs', f'eue_events_{t}i{iteration}.csv') ) From dfa9d277537055043c901d3c0fb861b0a7a41725 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Thu, 25 Jun 2026 16:49:39 -0600 Subject: [PATCH 19/24] drop already-modeled stress hours from eue/lole profiles used to identify new stress periods --- reeds/resource_adequacy/stress_periods.py | 52 ++++++++++++++++------- 1 file changed, 37 insertions(+), 15 deletions(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 929e9ec8..908e4291 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -257,18 +257,12 @@ def get_shortfall_periods( def get_longest_events( + dsmetric, sw, - t:int, - iteration:int, - hierarchy_level:str, - region:str, num_events:int=1, ): - ## Get the already-identified events - fpath = Path(sw.casedir, 'outputs', f'eue_events_{t}i{iteration}.csv') eue_events = ( - pd.read_csv(fpath, index_col=['level','region','number']) - .loc[hierarchy_level].loc[region] + get_events(dsmetric) .sort_values('timesteps', ascending=False) .head(num_events) ) @@ -378,7 +372,6 @@ def check_threshold_and_choose_periods( dfeue_agg, dflole_agg, dfenergy_agg, - stressperiods_this_iteration, ): ## NEUE Example: criterion = 'transgrp_1' hierarchy_level, metric_threshold = criterion.split('_') @@ -413,10 +406,17 @@ def check_threshold_and_choose_periods( for region in failed.index } case 'duration': + ## NOTE: This approach adds the next-longest event if the longest event is + ## already included as a stress period (since we've dropped the already- + ## included stress periods from dfeue_agg above). That's less likely to + ## help meet a max-duration threshold than adding the next-worst EUE day + ## is to help meet a NEUE threshold. So consider if there's something else + ## we should do if the longest-event day is already a stress period but + ## we're still not meeting the duration threshold. metric_periods = { region: get_longest_events( - sw=sw, t=t, iteration=iteration, - hierarchy_level=hierarchy_level, region=region, + dfeue_agg[region], sw=sw, + num_events=int(sw.GSw_PRM_StressIncrement), ) for region in failed.index } @@ -492,6 +492,23 @@ def get_stress_periods(case, sw, t, iteration): os.path.join( sw['casedir'], 'inputs_case', f'stress{t}i{iteration}', 'period_szn.csv') ) + stressperiods_this_iteration['start'] = ( + stressperiods_this_iteration.actual_period.map(reeds.timeseries.h2timestamp) + ) + stressperiods_this_iteration['end'] = ( + stressperiods_this_iteration['start'] + + ( + (pd.Timedelta('5D') if sw.GSw_HourlyType == 'wek' else pd.Timedelta('1D')) + - pd.Timedelta('1H') + ) + ) + ## Get already-modeled stress hours so we can exclude them from the hourly + ## EUE and LOLE profiles used to determine new stress periods + covered_hours = [ + pd.date_range(row.start, row.end, freq='1H') + for i,row in stressperiods_this_iteration.iterrows() + ] + covered_hours = [i for sublist in covered_hours for i in sublist] ### Check all stress criteria; for regions that fail, add new stress periods _failed = {} @@ -506,9 +523,9 @@ def get_stress_periods(case, sw, t, iteration): ## Example: criterion = 'transgrp_1' hierarchy_level, metric_threshold = criterion.split('_') rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) - dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum() - dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max() - dfenergy_agg = dfenergy.rename(columns=rmap).groupby(axis=1, level=0).sum() + dfeue_agg = dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum().drop(covered_hours) + dflole_agg = dflole.rename(columns=rmap).groupby(axis=1, level=0).max().drop(covered_hours) + dfenergy_agg = dfenergy.rename(columns=rmap).groupby(axis=1, level=0).sum().drop(covered_hours) ## Get the stress periods dictout = check_threshold_and_choose_periods( stress_metric, @@ -519,7 +536,6 @@ def get_stress_periods(case, sw, t, iteration): dfeue_agg, dflole_agg, dfenergy_agg, - stressperiods_this_iteration, ) if dictout is not None: _failed[stress_metric, criterion] = dictout['failed'] @@ -537,7 +553,13 @@ def get_stress_periods(case, sw, t, iteration): ).reset_index() print('All identified stress periods:') print(new_stress_periods) + ## Remove the existing stress periods and duplicates across regions + new_stress_periods = new_stress_periods.loc[ + ~new_stress_periods.period.isin(stressperiods_this_iteration.actual_period) + ].copy() new_stress_periods = new_stress_periods.drop_duplicates('period') + print('New stress periods after dropping duplicates:') + print(new_stress_periods) else: return {}, {}, {} From d88d298a9c3365c2dd4f5652edf55a5d44731723 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Thu, 25 Jun 2026 17:11:31 -0600 Subject: [PATCH 20/24] add reedsplots.plot_eue_events() --- postprocessing/single_case_plots.py | 24 ++++++++++++- reeds/reedsplots.py | 52 +++++++++++++++++++++++++++++ 2 files changed, 75 insertions(+), 1 deletion(-) diff --git a/postprocessing/single_case_plots.py b/postprocessing/single_case_plots.py index d5664336..dcc38b02 100644 --- a/postprocessing/single_case_plots.py +++ b/postprocessing/single_case_plots.py @@ -5,6 +5,7 @@ import sys import argparse import traceback +import itertools import cmocean sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) import reeds @@ -589,7 +590,28 @@ plt.close() print(savename) except Exception: - print('plot_stressperiod_neue failed:') + print('plot_ra_metrics_bylevel failed:') + print(traceback.format_exc()) + +try: + plt.close() + level = 'transgrp' + xvals = ['timesteps'] + yvals = ['mean','max'] + for xval, yval in itertools.product(xvals, yvals): + f, ax, _ = reedsplots.plot_eue_events( + case=case, year=year, level=level, + xval=xval, yval=yval, + ) + savename = f"plot_eue_events-{yval}_{xval}-{level}-{year}.png" + if write: + plt.savefig(os.path.join(savepath, savename)) + if interactive: + plt.show() + plt.close() + print(savename) +except Exception: + print('plot_eue_events failed:') print(traceback.format_exc()) try: diff --git a/reeds/reedsplots.py b/reeds/reedsplots.py index 6d0edd06..e3bd5658 100644 --- a/reeds/reedsplots.py +++ b/reeds/reedsplots.py @@ -6538,6 +6538,58 @@ def map_stressors( yield f, ax, dfout, plotlabel +def plot_eue_events( + case, year=None, level='transgrp', + xval:Literal['timesteps','mean','max','sum']='timesteps', + yval:Literal['timesteps','mean','max','sum']='mean', + scale=1, alpha=0.7, +): + """ + Plot event metrics against each other for all regions in specified level (columns) + and for each iteration in specified year (rows) + """ + sw = reeds.io.get_switches(case) + dflevel = reeds.io.get_dfmap(case)[level] + regions = dflevel.bounds.minx.sort_values().index + year = int(sw.endyear) if not year else year + iterations = get_stressperiods(case).loc[year].index.get_level_values('iteration').unique() + units = {'timesteps':'hours', 'mean':'MW', 'max':'MW', 'sum':'MWh'} + + nrows, ncols, coords = layout_subplots( + row_list=iterations, col_list=regions, + oneaxis=('columns' if len(regions) > 1 else 'rows'), + ) + dictout = {} + + plt.close() + f,ax = plt.subplots( + nrows, ncols, figsize=(scale*ncols, scale*nrows), + sharex=True, sharey='col', + ) + for iteration in iterations: + ## Get events + fpath = Path(case, 'outputs', f'eue_events_{year}i{iteration}.csv') + events = pd.read_csv(fpath, index_col=['level','region','number']).loc[level] + events['mean'] = events['sum'] / events['timesteps'] + dictout[iteration] = events + ## Plot for each region and iteration + for region in regions: + _ax = ax[coords[iteration, region]] if nrows + ncols > 2 else ax + _ax.plot( + events[xval], events[yval], lw=0, alpha=alpha, + marker='o', markersize=3, markeredgewidth=0, color='C3', + ) + ## Formatting + if iteration == 0: + _ax.set_title(region, weight='bold') + if iteration == max(iterations): + if region == regions[0]: + _ax.set_ylabel(f'EUE {yval} ({units[yval]})', y=0, ha='left') + _ax.set_xlabel(f'EUE {xval} ({units[xval]})', x=0, ha='left') + reeds.plots.despine(ax) + return f, ax, dictout + + def layout_subplots(row_list, col_list, oneaxis='columns'): """ Lay out series of row_list (e.g. cases) and col_list (e.g. years) into array of subplots, From 85acca53d1cb955e8a3f25a34764e2bf3d9ea27e Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Fri, 26 Jun 2026 11:05:34 -0600 Subject: [PATCH 21/24] adapt plot_stressperiod_evolution() for different RA metrics; fix plot_eue_events(); silence map_outagerate_new_stressperiods() if empty --- postprocessing/compare_cases.py | 907 ++++++++++---------- postprocessing/single_case_plots.py | 20 +- reeds/reedsplots.py | 111 ++- reeds/resource_adequacy/diagnostic_plots.py | 3 + 4 files changed, 549 insertions(+), 492 deletions(-) diff --git a/postprocessing/compare_cases.py b/postprocessing/compare_cases.py index a005ce65..a96fc8b2 100644 --- a/postprocessing/compare_cases.py +++ b/postprocessing/compare_cases.py @@ -2149,49 +2149,49 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): except Exception: print(traceback.format_exc()) -#%%### Interregional transfer capability to peak demand ratio -try: - f, ax, dfplot = reedsplots.plot_interreg_transfer_cap_ratio( - case=list(cases.values()), - colors={v: colors[k] for k,v in cases.items()}, - casenames={v:k for k,v in cases.items()}, - level='transreg', tstart=startyear, - ymax=None, - ) - ### Save it - slide = reeds.report_utils.add_to_pptx( - 'Interregional Transmission / Peak Demand', - prs=prs, - height=(SLIDE_HEIGHT if ax.shape[1] <= 8 else None), - width=(SLIDE_WIDTH if ax.shape[1] > 8 else None), - ) - if interactive: - plt.show() - -except Exception: - print(traceback.format_exc()) - - -#%%### Max firm imports -try: - f, ax, dfplot = reedsplots.plot_max_imports( - case=list(cases.values()), - colors={v: colors[k] for k,v in cases.items()}, - casenames={v:k for k,v in cases.items()}, - level='nercr', tstart=startyear, - ) - ### Save it - slide = reeds.report_utils.add_to_pptx( - 'Max Net Stress Imports / Peak Demand', - prs=prs, - height=(SLIDE_HEIGHT if ax.shape[1] <= 8 else None), - width=(SLIDE_WIDTH if ax.shape[1] > 8 else None), - ) - if interactive: - plt.show() - -except Exception: - print(traceback.format_exc()) +# #%%### Interregional transfer capability to peak demand ratio +# try: +# f, ax, dfplot = reedsplots.plot_interreg_transfer_cap_ratio( +# case=list(cases.values()), +# colors={v: colors[k] for k,v in cases.items()}, +# casenames={v:k for k,v in cases.items()}, +# level='transreg', tstart=startyear, +# ymax=None, +# ) +# ### Save it +# slide = reeds.report_utils.add_to_pptx( +# 'Interregional Transmission / Peak Demand', +# prs=prs, +# height=(SLIDE_HEIGHT if ax.shape[1] <= 8 else None), +# width=(SLIDE_WIDTH if ax.shape[1] > 8 else None), +# ) +# if interactive: +# plt.show() + +# except Exception: +# print(traceback.format_exc()) + + +# #%%### Max firm imports +# try: +# f, ax, dfplot = reedsplots.plot_max_imports( +# case=list(cases.values()), +# colors={v: colors[k] for k,v in cases.items()}, +# casenames={v:k for k,v in cases.items()}, +# level='nercr', tstart=startyear, +# ) +# ### Save it +# slide = reeds.report_utils.add_to_pptx( +# 'Max Net Stress Imports / Peak Demand', +# prs=prs, +# height=(SLIDE_HEIGHT if ax.shape[1] <= 8 else None), +# width=(SLIDE_WIDTH if ax.shape[1] > 8 else None), +# ) +# if interactive: +# plt.show() + +# except Exception: +# print(traceback.format_exc()) #%%### Transmission maps @@ -2303,206 +2303,202 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): except Exception: print(traceback.format_exc()) -#%% Flexibly sited load -try: - if any([float(dictin_sw[c].get('GSw_LoadSiteCF', 0)) for c in cases]): - loadsite_cases = { - k:v for k,v in cases.items() if float(dictin_sw[k].get('GSw_LoadSiteCF', 0)) - } - f, ax, dictplot = reeds.reedsplots.map_output_byyear( - case=loadsite_cases, - param='loadsite_cap', - years=[lastyear], - vscale=1e-3, - vmin=0, - title='Sited demand [GW]', - ) - ## Save it - slide = reeds.results.add_to_pptx('Flexibly Sited Demand', prs=prs, width=SLIDE_WIDTH) - if interactive: - plt.show() -except Exception: - print(traceback.format_exc()) - -#%%### RA sharing -if detailed: - try: - ralevel = 'transreg' - scale = 10 - wscale = 7e3 - dfmap = reeds.io.get_dfmap(cases[basecase]) - regions = dfmap[ralevel].bounds.minx.sort_values().index - - ### Calculate aggregated load - dictin_load_stress_agg = {} - for case in cases: - if int(dictin_sw[case].GSw_PRM_CapCredit): - hcol = 'ccseason' - df = dictin_peak_ccseason[case].copy() - else: - hcol = 'h' - df = dictin_load_stress[case].copy() - df = ( - df.assign(region=df.r.map(hierarchy[case][ralevel])) - .groupby(['t','region',hcol]).GW.sum() - .loc[lastyear].unstack('region') - ) - dictin_load_stress_agg[case] = df.sum() - - ### Calculate aggregated stress period flows - tran_flow_stress_agg = {} - for case in cases: - if int(dictin_sw[case].GSw_PRM_CapCredit): - df = dictin_prmtrade[case].copy() - hcol = 'ccseason' - else: - df = dictin_tran_flow_stress[case].copy() - hcol = 'h' - df['aggreg'] = df.r.map(hierarchy[case][ralevel]) - df['aggregg'] = df.rr.map(hierarchy[case][ralevel]) - df['interface'] = df.aggreg + '|' + df.aggregg - - df = ( - df - .loc[df.aggreg != df.aggregg] - .groupby(['t',hcol,'interface']).GW.sum().unstack('interface').fillna(0) - ) - if df.empty: - continue - else: - df = df.loc[lastyear].copy() - ## Order interfaces alphabetically - rename = {} - for interface in df: - r, rr = interface.split('|') - if r > rr: - rename[interface] = f'{rr}|{r}' - df[interface] *= -1 - df = df.rename(columns=rename).groupby(axis=1, level=0).sum() - ## Now reorder interfaces by flow - rename = {} - for interface in df: - r, rr = interface.split('|') - if df[interface].clip(lower=0).sum() < df[interface].clip(upper=0).abs().sum(): - rename[interface] = f'{rr}|{r}' - df[interface] *= -1 - tran_flow_stress_agg[case] = df.rename(columns=rename).copy() - - ### Calculate regional imports/exports - dfimportexport = {} - for case in cases: - df = {} - for region in regions: - df[region] = reedsplots.get_import_export( - region=region, df=tran_flow_stress_agg[case] - ) - dfimportexport[case] = pd.concat(df).sum(axis=1).unstack(level=0) - dfimportexport[case].columns = dfimportexport[case].columns.rename('region') - - ### Plot it - whiteout = dict(zip( - [f'C{i}' for i in range(10)], - [plt.cm.tab20(i*2+1) for i in range(10)] - )) - if any([v not in whiteout for v in list(colors.values())]): - whiteout = {v: (v[0], v[1], v[2], v[3]*0.7) for v in list(colors.values())} - - for label in ['max','average']: - plt.close() - f,ax = plt.subplots( - nrows, ncols, figsize=(SLIDE_WIDTH, SLIDE_HEIGHT), - gridspec_kw={'wspace':0.0,'hspace':-0.1}, - ) - - for case in cases: - ### Formatting - dfmap[ralevel].plot(ax=ax[coords[case]], facecolor='none', edgecolor='C7', lw=0.5) - dfmap['interconnect'].plot(ax=ax[coords[case]], facecolor='none', edgecolor='k', lw=1) - ax[coords[case]].set_title( - case, y=0.95, weight='bold', color=colors[case], fontsize=14) - ax[coords[case]].axis('off') - - ### RA flows - if case not in tran_flow_stress_agg: - continue - - if label == 'max': - ## Max flow - forwardwidth = tran_flow_stress_agg[case].clip(lower=0).max() - reversewidth = abs(tran_flow_stress_agg[case].clip(upper=0).min()) - else: - ## Average when it's flowing - forwardwidth = ( - tran_flow_stress_agg[case].clip(lower=0).sum() - / tran_flow_stress_agg[case].clip(lower=0).astype(bool).sum() - ) - reversewidth = ( - tran_flow_stress_agg[case].clip(upper=0).abs().sum() - / tran_flow_stress_agg[case].clip(upper=0).abs().astype(bool).sum() - ) - - interfaces = tran_flow_stress_agg[case].columns - numdays = ( - len(tran_flow_stress_agg[case]) - * int(dictin_sw[case].GSw_HourlyChunkLengthStress) - // 24 - ) - - ### Head/tail length: - gwh_forward = tran_flow_stress_agg[case].clip(lower=0).sum() - gwh_reverse = abs(tran_flow_stress_agg[case].clip(upper=0).sum()) - - reversefrac = gwh_reverse / (gwh_reverse + gwh_forward) - forwardfrac = gwh_forward / (gwh_reverse + gwh_forward) - - ### Plot it - for interface in interfaces: - r, rr = interface.split('|') - startx, starty = dfmap[ralevel].loc[r, ['x', 'y']] - endx, endy = dfmap[ralevel].loc[rr, ['x', 'y']] - - plots.plot_segmented_arrow( - ax[coords[case]], - reversefrac=reversefrac[interface], - forwardfrac=forwardfrac[interface], - reversewidth=reversewidth[interface]*wscale, - forwardwidth=forwardwidth[interface]*wscale, - startx=startx, starty=starty, endx=endx, endy=endy, - forwardcolor=colors[case], reversecolor=whiteout[colors[case]], - alpha=0.8, headwidthfrac=1.5, - ) - ### Scale - if scale: - (startx, starty, endx, endy) = (-2.0e6, -1.2e6, -1.5e6, -1.2e6) - yspan = ax[coords[case]].get_ylim() - yspan = yspan[1] - yspan[0] - plots.plot_segmented_arrow( - ax[coords[case]], - reversefrac=0, forwardfrac=1, - reversewidth=0, forwardwidth=scale*wscale, - startx=startx, starty=starty, endx=endx, endy=endy, - forwardcolor=colors[case], reversecolor=whiteout[colors[case]], - alpha=0.8, headwidthfrac=1.5, - ) - ax[coords[case]].annotate( - f"{scale} GW\n{label}", ((startx+endx)/2, starty-(scale/2*wscale)-yspan*0.02), - ha='center', va='top', fontsize=14, - ) - - ### Save it - title = f'{ralevel} {label} RA flows' - slide = reeds.report_utils.add_to_pptx(title, prs=prs) - if interactive: - plt.show() - except Exception: - print(traceback.format_exc()) +# #%% Flexibly sited load +# try: +# if any([float(dictin_sw[c].get('GSw_LoadSiteCF', 0)) for c in cases]): +# loadsite_cases = { +# k:v for k,v in cases.items() if float(dictin_sw[k].get('GSw_LoadSiteCF', 0)) +# } +# f, ax, dictplot = reeds.reedsplots.map_output_byyear( +# case=loadsite_cases, +# param='loadsite_cap', +# years=[lastyear], +# vscale=1e-3, +# vmin=0, +# title='Sited demand [GW]', +# ) +# ## Save it +# slide = reeds.results.add_to_pptx('Flexibly Sited Demand', prs=prs, width=SLIDE_WIDTH) +# if interactive: +# plt.show() +# except Exception: +# print(traceback.format_exc()) + +# #%%### RA sharing +# if detailed: +# try: +# ralevel = 'transreg' +# scale = 10 +# wscale = 7e3 +# dfmap = reeds.io.get_dfmap(cases[basecase]) +# regions = dfmap[ralevel].bounds.minx.sort_values().index + +# ### Calculate aggregated load +# dictin_load_stress_agg = {} +# for case in cases: +# if int(dictin_sw[case].GSw_PRM_CapCredit): +# hcol = 'ccseason' +# df = dictin_peak_ccseason[case].copy() +# else: +# hcol = 'h' +# df = dictin_load_stress[case].copy() +# df = ( +# df.assign(region=df.r.map(hierarchy[case][ralevel])) +# .groupby(['t','region',hcol]).GW.sum() +# .loc[lastyear].unstack('region') +# ) +# dictin_load_stress_agg[case] = df.sum() + +# ### Calculate aggregated stress period flows +# tran_flow_stress_agg = {} +# for case in cases: +# if int(dictin_sw[case].GSw_PRM_CapCredit): +# df = dictin_prmtrade[case].copy() +# hcol = 'ccseason' +# else: +# df = dictin_tran_flow_stress[case].copy() +# hcol = 'h' +# df['aggreg'] = df.r.map(hierarchy[case][ralevel]) +# df['aggregg'] = df.rr.map(hierarchy[case][ralevel]) +# df['interface'] = df.aggreg + '|' + df.aggregg + +# df = ( +# df +# .loc[df.aggreg != df.aggregg] +# .groupby(['t',hcol,'interface']).GW.sum().unstack('interface').fillna(0) +# ) +# if df.empty: +# continue +# else: +# df = df.loc[lastyear].copy() +# ## Order interfaces alphabetically +# rename = {} +# for interface in df: +# r, rr = interface.split('|') +# if r > rr: +# rename[interface] = f'{rr}|{r}' +# df[interface] *= -1 +# df = df.rename(columns=rename).groupby(axis=1, level=0).sum() +# ## Now reorder interfaces by flow +# rename = {} +# for interface in df: +# r, rr = interface.split('|') +# if df[interface].clip(lower=0).sum() < df[interface].clip(upper=0).abs().sum(): +# rename[interface] = f'{rr}|{r}' +# df[interface] *= -1 +# tran_flow_stress_agg[case] = df.rename(columns=rename).copy() + +# ### Calculate regional imports/exports +# dfimportexport = {} +# for case in cases: +# df = {} +# for region in regions: +# df[region] = reedsplots.get_import_export( +# region=region, df=tran_flow_stress_agg[case] +# ) +# dfimportexport[case] = pd.concat(df).sum(axis=1).unstack(level=0) +# dfimportexport[case].columns = dfimportexport[case].columns.rename('region') + +# ### Plot it +# whiteout = dict(zip( +# [f'C{i}' for i in range(10)], +# [plt.cm.tab20(i*2+1) for i in range(10)] +# )) +# if any([v not in whiteout for v in list(colors.values())]): +# whiteout = {v: (v[0], v[1], v[2], v[3]*0.7) for v in list(colors.values())} + +# for label in ['max','average']: +# plt.close() +# f,ax = plt.subplots( +# nrows, ncols, figsize=(SLIDE_WIDTH, SLIDE_HEIGHT), +# gridspec_kw={'wspace':0.0,'hspace':-0.1}, +# ) + +# for case in cases: +# ### Formatting +# dfmap[ralevel].plot(ax=ax[coords[case]], facecolor='none', edgecolor='C7', lw=0.5) +# dfmap['interconnect'].plot(ax=ax[coords[case]], facecolor='none', edgecolor='k', lw=1) +# ax[coords[case]].set_title( +# case, y=0.95, weight='bold', color=colors[case], fontsize=14) +# ax[coords[case]].axis('off') + +# ### RA flows +# if case not in tran_flow_stress_agg: +# continue + +# if label == 'max': +# ## Max flow +# forwardwidth = tran_flow_stress_agg[case].clip(lower=0).max() +# reversewidth = abs(tran_flow_stress_agg[case].clip(upper=0).min()) +# else: +# ## Average when it's flowing +# forwardwidth = ( +# tran_flow_stress_agg[case].clip(lower=0).sum() +# / tran_flow_stress_agg[case].clip(lower=0).astype(bool).sum() +# ) +# reversewidth = ( +# tran_flow_stress_agg[case].clip(upper=0).abs().sum() +# / tran_flow_stress_agg[case].clip(upper=0).abs().astype(bool).sum() +# ) + +# interfaces = tran_flow_stress_agg[case].columns +# numdays = ( +# len(tran_flow_stress_agg[case]) +# * int(dictin_sw[case].GSw_HourlyChunkLengthStress) +# // 24 +# ) + +# ### Head/tail length: +# gwh_forward = tran_flow_stress_agg[case].clip(lower=0).sum() +# gwh_reverse = abs(tran_flow_stress_agg[case].clip(upper=0).sum()) + +# reversefrac = gwh_reverse / (gwh_reverse + gwh_forward) +# forwardfrac = gwh_forward / (gwh_reverse + gwh_forward) + +# ### Plot it +# for interface in interfaces: +# r, rr = interface.split('|') +# startx, starty = dfmap[ralevel].loc[r, ['x', 'y']] +# endx, endy = dfmap[ralevel].loc[rr, ['x', 'y']] + +# plots.plot_segmented_arrow( +# ax[coords[case]], +# reversefrac=reversefrac[interface], +# forwardfrac=forwardfrac[interface], +# reversewidth=reversewidth[interface]*wscale, +# forwardwidth=forwardwidth[interface]*wscale, +# startx=startx, starty=starty, endx=endx, endy=endy, +# forwardcolor=colors[case], reversecolor=whiteout[colors[case]], +# alpha=0.8, headwidthfrac=1.5, +# ) +# ### Scale +# if scale: +# (startx, starty, endx, endy) = (-2.0e6, -1.2e6, -1.5e6, -1.2e6) +# yspan = ax[coords[case]].get_ylim() +# yspan = yspan[1] - yspan[0] +# plots.plot_segmented_arrow( +# ax[coords[case]], +# reversefrac=0, forwardfrac=1, +# reversewidth=0, forwardwidth=scale*wscale, +# startx=startx, starty=starty, endx=endx, endy=endy, +# forwardcolor=colors[case], reversecolor=whiteout[colors[case]], +# alpha=0.8, headwidthfrac=1.5, +# ) +# ax[coords[case]].annotate( +# f"{scale} GW\n{label}", ((startx+endx)/2, starty-(scale/2*wscale)-yspan*0.02), +# ha='center', va='top', fontsize=14, +# ) + +# ### Save it +# title = f'{ralevel} {label} RA flows' +# slide = reeds.report_utils.add_to_pptx(title, prs=prs) +# if interactive: +# plt.show() +# except Exception: +# print(traceback.format_exc()) #%%### Copy some premade single-case plots -# Use first stress metric level -## TODO: add a check for choosing level if there are multiple stress metrics -stress_metrics = dictin_sw[basecase]['GSw_PRM_StressThresholdMetrics'].split('/') -level = dictin_sw[basecase][f'GSw_PRM_StressThreshold{stress_metrics[0]}'].split('_')[0] wide = 1 if len(hierarchy[basecase]['transreg'].unique()) > 6 else 0 weatheryear = sw.GSw_HourlyWeatherYears.split('_')[0] metrics = [ @@ -2517,9 +2513,10 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): ] for figname, width, height in [ (f'map_gencap_transcap-{lastyear}', None, SLIDE_HEIGHT), - (f'plot_stressperiod_evolution-neue-{level}', SLIDE_WIDTH, None), + ('plot_stressperiod_evolution-neue', SLIDE_WIDTH, None), (f'plot_dispatch-yearbymonth-1-{lastyear}', SLIDE_WIDTH, None), ## Include both versions for backwards compatibility + ('plot_stressperiod_evolution-sum-transgrp', SLIDE_WIDTH, None), (f'plot_dispatch-yearbymonth-1-{lastyear}-w{weatheryear}', SLIDE_WIDTH, None), ] + [ ( @@ -2541,220 +2538,220 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): print(f'No outputs/figures/{figname}.png for {os.path.basename(cases[case])}') -#%%### Generation capacity maps - -### Shared data -try: - base = cases[list(cases.keys())[0]] - val_r = dictin_cap_r[basecase].r.unique() - dfmap = reeds.io.get_dfmap(base) - dfba = dfmap['r'] - dfstates = dfmap['st'] - if (len(cases) == 2) and (not forcemulti): - for i_plot in maptechs.keys(): - plt.close() - f,ax=plt.subplots( - 1, 3, sharex=True, sharey=True, figsize=(14,8), - gridspec_kw={'wspace':-0.05, 'hspace':0.05}, - dpi=150, - ) - - _,_,dfplot = reedsplots.plot_diff_maps( - val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], - year=lastyear, casebase=casebase, casecomp=casecomp, - plot='base', f=f, ax=ax[0], - cmap=cmocean.cm.rain, - ) - ax[0].annotate( - casebase_name, - (0.1,1), xycoords='axes fraction', fontsize=10) - - _,_,dfplot = reedsplots.plot_diff_maps( - val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], - year=lastyear, casebase=casebase, casecomp=casecomp, - plot='comp', f=f, ax=ax[1], - cmap=cmocean.cm.rain, - ) - ax[1].annotate( - casecomp_name, - (0.1,1), xycoords='axes fraction', fontsize=10) - - _,_,dfplot = reedsplots.plot_diff_maps( - val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], - year=lastyear, casebase=casebase, casecomp=casecomp, - plot='absdiff', f=f, ax=ax[2], - cmap=plt.cm.RdBu_r, - ) - ax[2].annotate( - '{}\n– {}'.format( - casecomp_name, - casebase_name), - (0.1,1), xycoords='axes fraction', fontsize=10) - - reeds.report_utils.add_to_pptx(f'{i_plot} Capacity {lastyear} [GW]', prs=prs) - if interactive: - plt.show() - else: - figwidth = SLIDE_WIDTH - #### Absolute maps - if (nrows == 1) or (ncols == 1): - legendcoords = max(nrows, ncols) - 1 - elif (nrows-1, ncols-1) in coords.values(): - legendcoords = (nrows-1, ncols-1) - else: - legendcoords = (nrows-2, ncols-1) - - ### Set up plot - for tech_type in maptechs.keys(): - ### Get limits - vmin = 0. - vmax = float(pd.concat({ - case: dictin_cap_r[case].loc[ - (dictin_cap_r[case].i.isin(maptechs[tech_type])) - & (dictin_cap_r[case].t.astype(int)==lastyear) - ].groupby('r').MW.sum() - for case in cases - }).max()) / 1e3 - if np.isnan(vmax): - vmax = 0. - if not vmax: - print(f'{tech_type} has zero capacity in {lastyear}, so skipping maps') - continue - ### Set up plot - plt.close() - f,ax = plt.subplots( - nrows, ncols, figsize=(figwidth, SLIDE_HEIGHT), - gridspec_kw={'wspace':0.0,'hspace':-0.1}, - ) - ### Plot it - for case in cases: - dfval = dictin_cap_r[case].loc[ - (dictin_cap_r[case].i.isin(maptechs[tech_type])) - & (dictin_cap_r[case].t.astype(int)==lastyear) - ].groupby('r').MW.sum() - dfplot = dfba.copy() - dfplot['GW'] = (dfval / 1e3).fillna(0) - - ax[coords[case]].set_title( - case if nowrap else plots.wraptext(case, width=figwidth/ncols*0.9, fontsize=14) - ) - dfba.plot( - ax=ax[coords[case]], - facecolor='none', edgecolor='k', lw=0.1, zorder=10000) - dfstates.plot( - ax=ax[coords[case]], - facecolor='none', edgecolor='k', lw=0.2, zorder=10001) - dfplot.plot( - ax=ax[coords[case]], column='GW', cmap=cmap, vmin=vmin, vmax=vmax, - legend=False, - ) - ## Legend - if coords[case] == legendcoords: - plots.addcolorbarhist( - f=f, ax0=ax[coords[case]], data=dfplot.GW.values, - title=f'{tech_type} {lastyear}\ncapacity [GW]', cmap=cmap, vmin=vmin, vmax=vmax, - orientation='horizontal', labelpad=2.25, histratio=0., - cbarwidth=0.05, cbarheight=0.85, - cbarbottom=-0.05, cbarhoffset=0., - ) - - for row in range(nrows): - for col in range(ncols): - if nrows == 1: - ax[col].axis('off') - elif ncols == 1: - ax[row].axis('off') - else: - ax[row,col].axis('off') - ### Save it - slide = reeds.report_utils.add_to_pptx(f'{tech_type} Capacity {lastyear} [GW]', prs=prs) - if interactive: - plt.show() - - #### Difference maps - ### Set up plot - for tech_type in maptechs.keys(): - ### Get limits - dfval = pd.concat({ - case: dictin_cap_r[case].loc[ - (dictin_cap_r[case].i.isin(maptechs[tech_type])) - & (dictin_cap_r[case].t.astype(int)==lastyear) - ].groupby('r').MW.sum() - for case in cases - }, axis=1).fillna(0) / 1e3 - dfdiff = dfval.subtract(dfval[basecase], axis=0) - ### Get colorbar limits - absmax = dfval.stack().max() - diffmax = dfdiff.unstack().abs().max() - - if np.isnan(absmax): - absmax = 0. - if not absmax: - print(f'{tech_type} has zero capacity in {lastyear}, so skipping maps') - continue - ### Set up plot - plt.close() - f,ax = plt.subplots( - nrows, ncols, figsize=(figwidth, SLIDE_HEIGHT), - gridspec_kw={'wspace':0.0,'hspace':-0.1}, - ) - ### Plot it - for case in cases: - dfplot = dfba.copy() - dfplot['GW'] = dfval[case] if case == basecase else dfdiff[case] - - ax[coords[case]].set_title( - case if nowrap else plots.wraptext(case, width=figwidth/ncols*0.9, fontsize=14) - ) - dfba.plot( - ax=ax[coords[case]], - facecolor='none', edgecolor='k', lw=0.1, zorder=10000) - dfstates.plot( - ax=ax[coords[case]], - facecolor='none', edgecolor='k', lw=0.2, zorder=10001) - dfplot.plot( - ax=ax[coords[case]], column='GW', - cmap=(cmap if case == basecase else cmap_diff), - vmin=(0 if case == basecase else -diffmax), - vmax=(absmax if case == basecase else diffmax), - legend=False, - ) - ## Difference legend - if coords[case] == legendcoords: - plots.addcolorbarhist( - f=f, ax0=ax[coords[case]], data=dfplot.GW.values, - title=f'{tech_type} {lastyear}\ncapacity, difference\nfrom {basecase} [GW]', - cmap=(cmap if case == basecase else cmap_diff), - vmin=(0 if case == basecase else -diffmax), - vmax=(absmax if case == basecase else diffmax), - orientation='horizontal', labelpad=2.25, histratio=0., - cbarwidth=0.05, cbarheight=0.85, - cbarbottom=-0.05, cbarhoffset=0., - ) - ## Absolute legend - plots.addcolorbarhist( - f=f, ax0=ax[coords[basecase]], data=dfval[basecase].values, - title=f'{tech_type} {lastyear}\ncapacity [GW]', - cmap=cmap, vmin=0, vmax=absmax, - orientation='horizontal', labelpad=2.25, histratio=0., - cbarwidth=0.05, cbarheight=0.85, - cbarbottom=-0.05, cbarhoffset=0., - ) - - for row in range(nrows): - for col in range(ncols): - if nrows == 1: - ax[col].axis('off') - elif ncols == 1: - ax[row].axis('off') - else: - ax[row,col].axis('off') - ### Save it - slide = reeds.report_utils.add_to_pptx(f'Difference: {tech_type} Capacity {lastyear} [GW]', prs=prs) - if interactive: - plt.show() -except Exception: - print(traceback.format_exc()) +# #%%### Generation capacity maps + +# ### Shared data +# try: +# base = cases[list(cases.keys())[0]] +# val_r = dictin_cap_r[basecase].r.unique() +# dfmap = reeds.io.get_dfmap(base) +# dfba = dfmap['r'] +# dfstates = dfmap['st'] +# if (len(cases) == 2) and (not forcemulti): +# for i_plot in maptechs.keys(): +# plt.close() +# f,ax=plt.subplots( +# 1, 3, sharex=True, sharey=True, figsize=(14,8), +# gridspec_kw={'wspace':-0.05, 'hspace':0.05}, +# dpi=150, +# ) + +# _,_,dfplot = reedsplots.plot_diff_maps( +# val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], +# year=lastyear, casebase=casebase, casecomp=casecomp, +# plot='base', f=f, ax=ax[0], +# cmap=cmocean.cm.rain, +# ) +# ax[0].annotate( +# casebase_name, +# (0.1,1), xycoords='axes fraction', fontsize=10) + +# _,_,dfplot = reedsplots.plot_diff_maps( +# val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], +# year=lastyear, casebase=casebase, casecomp=casecomp, +# plot='comp', f=f, ax=ax[1], +# cmap=cmocean.cm.rain, +# ) +# ax[1].annotate( +# casecomp_name, +# (0.1,1), xycoords='axes fraction', fontsize=10) + +# _,_,dfplot = reedsplots.plot_diff_maps( +# val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], +# year=lastyear, casebase=casebase, casecomp=casecomp, +# plot='absdiff', f=f, ax=ax[2], +# cmap=plt.cm.RdBu_r, +# ) +# ax[2].annotate( +# '{}\n– {}'.format( +# casecomp_name, +# casebase_name), +# (0.1,1), xycoords='axes fraction', fontsize=10) + +# reeds.report_utils.add_to_pptx(f'{i_plot} Capacity {lastyear} [GW]', prs=prs) +# if interactive: +# plt.show() +# else: +# figwidth = SLIDE_WIDTH +# #### Absolute maps +# if (nrows == 1) or (ncols == 1): +# legendcoords = max(nrows, ncols) - 1 +# elif (nrows-1, ncols-1) in coords.values(): +# legendcoords = (nrows-1, ncols-1) +# else: +# legendcoords = (nrows-2, ncols-1) + +# ### Set up plot +# for tech_type in maptechs.keys(): +# ### Get limits +# vmin = 0. +# vmax = float(pd.concat({ +# case: dictin_cap_r[case].loc[ +# (dictin_cap_r[case].i.isin(maptechs[tech_type])) +# & (dictin_cap_r[case].t.astype(int)==lastyear) +# ].groupby('r').MW.sum() +# for case in cases +# }).max()) / 1e3 +# if np.isnan(vmax): +# vmax = 0. +# if not vmax: +# print(f'{tech_type} has zero capacity in {lastyear}, so skipping maps') +# continue +# ### Set up plot +# plt.close() +# f,ax = plt.subplots( +# nrows, ncols, figsize=(figwidth, SLIDE_HEIGHT), +# gridspec_kw={'wspace':0.0,'hspace':-0.1}, +# ) +# ### Plot it +# for case in cases: +# dfval = dictin_cap_r[case].loc[ +# (dictin_cap_r[case].i.isin(maptechs[tech_type])) +# & (dictin_cap_r[case].t.astype(int)==lastyear) +# ].groupby('r').MW.sum() +# dfplot = dfba.copy() +# dfplot['GW'] = (dfval / 1e3).fillna(0) + +# ax[coords[case]].set_title( +# case if nowrap else plots.wraptext(case, width=figwidth/ncols*0.9, fontsize=14) +# ) +# dfba.plot( +# ax=ax[coords[case]], +# facecolor='none', edgecolor='k', lw=0.1, zorder=10000) +# dfstates.plot( +# ax=ax[coords[case]], +# facecolor='none', edgecolor='k', lw=0.2, zorder=10001) +# dfplot.plot( +# ax=ax[coords[case]], column='GW', cmap=cmap, vmin=vmin, vmax=vmax, +# legend=False, +# ) +# ## Legend +# if coords[case] == legendcoords: +# plots.addcolorbarhist( +# f=f, ax0=ax[coords[case]], data=dfplot.GW.values, +# title=f'{tech_type} {lastyear}\ncapacity [GW]', cmap=cmap, vmin=vmin, vmax=vmax, +# orientation='horizontal', labelpad=2.25, histratio=0., +# cbarwidth=0.05, cbarheight=0.85, +# cbarbottom=-0.05, cbarhoffset=0., +# ) + +# for row in range(nrows): +# for col in range(ncols): +# if nrows == 1: +# ax[col].axis('off') +# elif ncols == 1: +# ax[row].axis('off') +# else: +# ax[row,col].axis('off') +# ### Save it +# slide = reeds.report_utils.add_to_pptx(f'{tech_type} Capacity {lastyear} [GW]', prs=prs) +# if interactive: +# plt.show() + +# #### Difference maps +# ### Set up plot +# for tech_type in maptechs.keys(): +# ### Get limits +# dfval = pd.concat({ +# case: dictin_cap_r[case].loc[ +# (dictin_cap_r[case].i.isin(maptechs[tech_type])) +# & (dictin_cap_r[case].t.astype(int)==lastyear) +# ].groupby('r').MW.sum() +# for case in cases +# }, axis=1).fillna(0) / 1e3 +# dfdiff = dfval.subtract(dfval[basecase], axis=0) +# ### Get colorbar limits +# absmax = dfval.stack().max() +# diffmax = dfdiff.unstack().abs().max() + +# if np.isnan(absmax): +# absmax = 0. +# if not absmax: +# print(f'{tech_type} has zero capacity in {lastyear}, so skipping maps') +# continue +# ### Set up plot +# plt.close() +# f,ax = plt.subplots( +# nrows, ncols, figsize=(figwidth, SLIDE_HEIGHT), +# gridspec_kw={'wspace':0.0,'hspace':-0.1}, +# ) +# ### Plot it +# for case in cases: +# dfplot = dfba.copy() +# dfplot['GW'] = dfval[case] if case == basecase else dfdiff[case] + +# ax[coords[case]].set_title( +# case if nowrap else plots.wraptext(case, width=figwidth/ncols*0.9, fontsize=14) +# ) +# dfba.plot( +# ax=ax[coords[case]], +# facecolor='none', edgecolor='k', lw=0.1, zorder=10000) +# dfstates.plot( +# ax=ax[coords[case]], +# facecolor='none', edgecolor='k', lw=0.2, zorder=10001) +# dfplot.plot( +# ax=ax[coords[case]], column='GW', +# cmap=(cmap if case == basecase else cmap_diff), +# vmin=(0 if case == basecase else -diffmax), +# vmax=(absmax if case == basecase else diffmax), +# legend=False, +# ) +# ## Difference legend +# if coords[case] == legendcoords: +# plots.addcolorbarhist( +# f=f, ax0=ax[coords[case]], data=dfplot.GW.values, +# title=f'{tech_type} {lastyear}\ncapacity, difference\nfrom {basecase} [GW]', +# cmap=(cmap if case == basecase else cmap_diff), +# vmin=(0 if case == basecase else -diffmax), +# vmax=(absmax if case == basecase else diffmax), +# orientation='horizontal', labelpad=2.25, histratio=0., +# cbarwidth=0.05, cbarheight=0.85, +# cbarbottom=-0.05, cbarhoffset=0., +# ) +# ## Absolute legend +# plots.addcolorbarhist( +# f=f, ax0=ax[coords[basecase]], data=dfval[basecase].values, +# title=f'{tech_type} {lastyear}\ncapacity [GW]', +# cmap=cmap, vmin=0, vmax=absmax, +# orientation='horizontal', labelpad=2.25, histratio=0., +# cbarwidth=0.05, cbarheight=0.85, +# cbarbottom=-0.05, cbarhoffset=0., +# ) + +# for row in range(nrows): +# for col in range(ncols): +# if nrows == 1: +# ax[col].axis('off') +# elif ncols == 1: +# ax[row].axis('off') +# else: +# ax[row,col].axis('off') +# ### Save it +# slide = reeds.report_utils.add_to_pptx(f'Difference: {tech_type} Capacity {lastyear} [GW]', prs=prs) +# if interactive: +# plt.show() +# except Exception: +# print(traceback.format_exc()) #%% Save the powerpoint file prs.save(savename) diff --git a/postprocessing/single_case_plots.py b/postprocessing/single_case_plots.py index dcc38b02..b64a5b59 100644 --- a/postprocessing/single_case_plots.py +++ b/postprocessing/single_case_plots.py @@ -739,16 +739,16 @@ print(traceback.format_exc()) try: - level, threshold = sw['GSw_PRM_StressThresholdNEUE'].split('/')[0].split('_') - plt.close() - f,ax = reedsplots.plot_stressperiod_evolution(case=case, level=level) - savename = f'plot_stressperiod_evolution-neue-{level}.png' - if write: - plt.savefig(os.path.join(savepath, savename)) - if interactive: - plt.show() - plt.close() - print(savename) + for metric in ['neue','depth','duration','lolh','lole','lold']: + plt.close() + f,ax = reedsplots.plot_stressperiod_evolution(case=case, metric=metric) + savename = f'plot_stressperiod_evolution-{metric}.png' + if write: + plt.savefig(os.path.join(savepath, savename)) + if interactive: + plt.show() + plt.close() + print(savename) except Exception: print('plot_stressperiod_evolution failed:') print(traceback.format_exc()) diff --git a/reeds/reedsplots.py b/reeds/reedsplots.py index e3bd5658..816fb19b 100644 --- a/reeds/reedsplots.py +++ b/reeds/reedsplots.py @@ -4121,32 +4121,48 @@ def plot_stressperiod_days(case, repcolor='k', sharey=False, figsize=(10,5)): def plot_stressperiod_evolution( - case, level=None, threshold=None, + case, + metric:Literal['neue','depth','duration','lolh','lole','lold']='neue', figsize=None, scale_widths=False, ): - """Plot NEUE by year and stress period iteration""" + """ + Plot RA metric by year and stress period iteration. + RA metric thresholds are included in the plot even if they are not applied. + Only the first /-delimited level_threshold pair is used. + """ + from reeds.resource_adequacy import stress_periods ### Parse inputs sw = reeds.io.get_switches(case) - _level, _threshold = sw['GSw_PRM_StressThresholdNEUE'].split('/')[0].split('_') - level = _level if level is None else level - threshold = float(_threshold) if threshold is None else threshold - ### Load NEUE results + switch = stress_periods.RA_SWITCHES[metric.lower()] + switch_metric = stress_periods.SWITCH_METRIC[metric.lower()] + level, threshold = sw[switch].split('/')[0].split('_') + threshold = float(threshold) + ### Load RA results infiles = sorted(glob(os.path.join(case,'outputs','ra_metrics_*.csv'))) - dictin_neue = { + dictin_ra = { tuple([int(x) for x in os.path.basename(f)[len('ra_metrics_'):-len('.csv')].split('i')]): pd.read_csv(f, index_col=['level','metric','region']) for f in infiles } ## Reshape to (year,iteration) x (region) dfplot = ( - pd.concat(dictin_neue, names=['year','iteration']) + pd.concat(dictin_ra, names=['year','iteration']) .xs(level,0,'level') - .xs('neue_ppm',0,'metric') + .xs(switch_metric,0,'metric') .squeeze(1).unstack('region') ) ### Load stress periods for labels dfstress = get_stressperiods(case) ### Plot setup + ylabel = { + 'neue': 'NEUE [ppm]', + 'lolh': 'LOLH [event-hours/year]', + 'lole': 'LOLE [events/year]', + 'lold': 'LOLD [event-days/year]', + 'duration': 'Duration [hours]', + 'depth': 'Depth [% of peak]', + } + scale = {'depth':100} years = [ y for y in dfplot.index.get_level_values('year').unique() if y >= int(sw.GSw_StartMarkets) @@ -4165,7 +4181,7 @@ def plot_stressperiod_evolution( gridspec_kw=gridspec_kw, ) for col, year in enumerate(years): - df = dfplot.loc[year] + df = dfplot.loc[year] * scale.get(metric, 1) for region in regions: ax[col].plot( df.index, df[region], @@ -4196,7 +4212,7 @@ def plot_stressperiod_evolution( loc='upper left', bbox_to_anchor=(1,1), frameon=False, handletextpad=0.3, handlelength=0.7, ) - ax[0].set_ylabel('NEUE [ppm]') + ax[0].set_ylabel(ylabel[metric]) ax[0].set_ylim(0) plots.despine(ax) @@ -4204,7 +4220,7 @@ def plot_stressperiod_evolution( def plot_ra_metrics_bylevel( - case, tmin=2023, + case, tmin=2026, levels=['country','interconnect','transreg','transgrp'], metrics=[ 'neue', @@ -4244,7 +4260,7 @@ def plot_ra_metrics_bylevel( .set_index(['level','metric','region']).squeeze(1) ) dfin_ra = pd.concat(ra_metrics, axis=0, names=['year']).unstack('year') - dfin_ra = dfin_ra[[c for c in dfin_ra if int(c) >= 2025]].copy() + dfin_ra = dfin_ra[[c for c in dfin_ra if int(c) >= tmin]].copy() if onlydata: return dfin_ra ### Plot settings @@ -6542,7 +6558,8 @@ def plot_eue_events( case, year=None, level='transgrp', xval:Literal['timesteps','mean','max','sum']='timesteps', yval:Literal['timesteps','mean','max','sum']='mean', - scale=1, alpha=0.7, + scale=1, alpha=0.9, + showhull=True, ): """ Plot event metrics against each other for all regions in specified level (columns) @@ -6553,18 +6570,30 @@ def plot_eue_events( regions = dflevel.bounds.minx.sort_values().index year = int(sw.endyear) if not year else year iterations = get_stressperiods(case).loc[year].index.get_level_values('iteration').unique() - units = {'timesteps':'hours', 'mean':'MW', 'max':'MW', 'sum':'MWh'} - + units = {'timesteps':'hours', 'mean':'%', 'max':'%', 'sum':'MWh'} + if showhull: + import scipy.spatial + ## Get load if we need it for normalization + peakload = ( + pd.read_csv(Path(case,'inputs_case','peakload.csv'), index_col=['level','region']) + .loc[level][str(year)] + ) + ## Get thresholds + thresholds = {} + if level in sw.GSw_PRM_StressThresholdDuration: + thresholds['timesteps'] = int(sw.GSw_PRM_StressThresholdDuration.split('/')[0].split('_')[-1]) + if level in sw.GSw_PRM_StressThresholdDepth: + thresholds['max'] = float(sw.GSw_PRM_StressThresholdDepth.split('/')[0].split('_')[-1]) * 100 + ## Plot it nrows, ncols, coords = layout_subplots( row_list=iterations, col_list=regions, oneaxis=('columns' if len(regions) > 1 else 'rows'), ) dictout = {} - plt.close() f,ax = plt.subplots( nrows, ncols, figsize=(scale*ncols, scale*nrows), - sharex=True, sharey='col', + sharex=True, sharey=True, ) for iteration in iterations: ## Get events @@ -6575,17 +6604,45 @@ def plot_eue_events( ## Plot for each region and iteration for region in regions: _ax = ax[coords[iteration, region]] if nrows + ncols > 2 else ax - _ax.plot( - events[xval], events[yval], lw=0, alpha=alpha, - marker='o', markersize=3, markeredgewidth=0, color='C3', - ) + if region in events.index.get_level_values('region'): + points = {'x':events.loc[region][xval], 'y':events.loc[region][yval]} + if xval in ['max', 'mean']: + points['x'] = points['x'] / peakload[region] * 100 + if yval in ['max', 'mean']: + points['y'] = points['y'] / peakload[region] * 100 + _ax.plot( + points['x'], points['y'], lw=0, alpha=alpha, + marker='o', markersize=3, markeredgewidth=0, color='C3', + ) + if showhull and len(points['x']) >= 3: + dfpoints = pd.DataFrame(points) + hull = scipy.spatial.ConvexHull(dfpoints) + dfhull = dfpoints.loc[hull.vertices] + _ax.fill(dfhull.x, dfhull.y, color='C3', lw=0, alpha=0.3, zorder=-1) ## Formatting if iteration == 0: - _ax.set_title(region, weight='bold') - if iteration == max(iterations): - if region == regions[0]: - _ax.set_ylabel(f'EUE {yval} ({units[yval]})', y=0, ha='left') - _ax.set_xlabel(f'EUE {xval} ({units[xval]})', x=0, ha='left') + _ax.set_title(region.replace('_','\n'), weight='bold') + if region == regions[0]: + _ax.annotate( + f'{year}i{iteration}', + (0, 1), xycoords='axes fraction', + xytext=(3,-3), textcoords='offset points', + ha='left', va='top', fontsize=12, + path_effects=[pe.withStroke(linewidth=1.5, foreground='w', alpha=0.8)], + ) + if iteration == max(iterations): + _ax.set_ylabel(f'EUE {yval} [{units[yval]}]', y=0, ha='left') + _ax.set_xlabel(f'EUE {xval} [{units[xval]}]', x=0, ha='left') + ## Formatting + _ax.set_xlim(0, _ax.get_xlim()[1]) + _ax.set_ylim(0, _ax.get_ylim()[1]) + for iteration in iterations: + for region in regions: + _ax = ax[coords[iteration, region]] if nrows + ncols > 2 else ax + if xval in thresholds: + _ax.axvline(thresholds[xval], c='k', ls=':', lw=0.75) + if yval in thresholds: + _ax.axhline(thresholds[yval], c='k', ls=':', lw=0.75) reeds.plots.despine(ax) return f, ax, dictout diff --git a/reeds/resource_adequacy/diagnostic_plots.py b/reeds/resource_adequacy/diagnostic_plots.py index 1b4e6b01..a4487d39 100644 --- a/reeds/resource_adequacy/diagnostic_plots.py +++ b/reeds/resource_adequacy/diagnostic_plots.py @@ -1030,6 +1030,9 @@ def map_pras_failure_rate(sw, dfs, aggfunc='mean', repair=False): def map_outagerate_new_stressperiods(sw, dfs): new_stress_periods = dfs['new_stress_periods'] + if not len(new_stress_periods): + print('No new stress periods to plot') + return dates = ( new_stress_periods .period.map(reeds.timeseries.h2timestamp) From 77803852be6d86e666b25d2ecc22a1e5a9d24ee7 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Fri, 26 Jun 2026 11:07:18 -0600 Subject: [PATCH 22/24] revert compare_cases.py bypass --- postprocessing/compare_cases.py | 900 ++++++++++++++++---------------- 1 file changed, 450 insertions(+), 450 deletions(-) diff --git a/postprocessing/compare_cases.py b/postprocessing/compare_cases.py index a96fc8b2..df648cce 100644 --- a/postprocessing/compare_cases.py +++ b/postprocessing/compare_cases.py @@ -2149,49 +2149,49 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): except Exception: print(traceback.format_exc()) -# #%%### Interregional transfer capability to peak demand ratio -# try: -# f, ax, dfplot = reedsplots.plot_interreg_transfer_cap_ratio( -# case=list(cases.values()), -# colors={v: colors[k] for k,v in cases.items()}, -# casenames={v:k for k,v in cases.items()}, -# level='transreg', tstart=startyear, -# ymax=None, -# ) -# ### Save it -# slide = reeds.report_utils.add_to_pptx( -# 'Interregional Transmission / Peak Demand', -# prs=prs, -# height=(SLIDE_HEIGHT if ax.shape[1] <= 8 else None), -# width=(SLIDE_WIDTH if ax.shape[1] > 8 else None), -# ) -# if interactive: -# plt.show() - -# except Exception: -# print(traceback.format_exc()) - - -# #%%### Max firm imports -# try: -# f, ax, dfplot = reedsplots.plot_max_imports( -# case=list(cases.values()), -# colors={v: colors[k] for k,v in cases.items()}, -# casenames={v:k for k,v in cases.items()}, -# level='nercr', tstart=startyear, -# ) -# ### Save it -# slide = reeds.report_utils.add_to_pptx( -# 'Max Net Stress Imports / Peak Demand', -# prs=prs, -# height=(SLIDE_HEIGHT if ax.shape[1] <= 8 else None), -# width=(SLIDE_WIDTH if ax.shape[1] > 8 else None), -# ) -# if interactive: -# plt.show() - -# except Exception: -# print(traceback.format_exc()) +#%%### Interregional transfer capability to peak demand ratio +try: + f, ax, dfplot = reedsplots.plot_interreg_transfer_cap_ratio( + case=list(cases.values()), + colors={v: colors[k] for k,v in cases.items()}, + casenames={v:k for k,v in cases.items()}, + level='transreg', tstart=startyear, + ymax=None, + ) + ### Save it + slide = reeds.report_utils.add_to_pptx( + 'Interregional Transmission / Peak Demand', + prs=prs, + height=(SLIDE_HEIGHT if ax.shape[1] <= 8 else None), + width=(SLIDE_WIDTH if ax.shape[1] > 8 else None), + ) + if interactive: + plt.show() + +except Exception: + print(traceback.format_exc()) + + +#%%### Max firm imports +try: + f, ax, dfplot = reedsplots.plot_max_imports( + case=list(cases.values()), + colors={v: colors[k] for k,v in cases.items()}, + casenames={v:k for k,v in cases.items()}, + level='nercr', tstart=startyear, + ) + ### Save it + slide = reeds.report_utils.add_to_pptx( + 'Max Net Stress Imports / Peak Demand', + prs=prs, + height=(SLIDE_HEIGHT if ax.shape[1] <= 8 else None), + width=(SLIDE_WIDTH if ax.shape[1] > 8 else None), + ) + if interactive: + plt.show() + +except Exception: + print(traceback.format_exc()) #%%### Transmission maps @@ -2303,199 +2303,199 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): except Exception: print(traceback.format_exc()) -# #%% Flexibly sited load -# try: -# if any([float(dictin_sw[c].get('GSw_LoadSiteCF', 0)) for c in cases]): -# loadsite_cases = { -# k:v for k,v in cases.items() if float(dictin_sw[k].get('GSw_LoadSiteCF', 0)) -# } -# f, ax, dictplot = reeds.reedsplots.map_output_byyear( -# case=loadsite_cases, -# param='loadsite_cap', -# years=[lastyear], -# vscale=1e-3, -# vmin=0, -# title='Sited demand [GW]', -# ) -# ## Save it -# slide = reeds.results.add_to_pptx('Flexibly Sited Demand', prs=prs, width=SLIDE_WIDTH) -# if interactive: -# plt.show() -# except Exception: -# print(traceback.format_exc()) - -# #%%### RA sharing -# if detailed: -# try: -# ralevel = 'transreg' -# scale = 10 -# wscale = 7e3 -# dfmap = reeds.io.get_dfmap(cases[basecase]) -# regions = dfmap[ralevel].bounds.minx.sort_values().index - -# ### Calculate aggregated load -# dictin_load_stress_agg = {} -# for case in cases: -# if int(dictin_sw[case].GSw_PRM_CapCredit): -# hcol = 'ccseason' -# df = dictin_peak_ccseason[case].copy() -# else: -# hcol = 'h' -# df = dictin_load_stress[case].copy() -# df = ( -# df.assign(region=df.r.map(hierarchy[case][ralevel])) -# .groupby(['t','region',hcol]).GW.sum() -# .loc[lastyear].unstack('region') -# ) -# dictin_load_stress_agg[case] = df.sum() - -# ### Calculate aggregated stress period flows -# tran_flow_stress_agg = {} -# for case in cases: -# if int(dictin_sw[case].GSw_PRM_CapCredit): -# df = dictin_prmtrade[case].copy() -# hcol = 'ccseason' -# else: -# df = dictin_tran_flow_stress[case].copy() -# hcol = 'h' -# df['aggreg'] = df.r.map(hierarchy[case][ralevel]) -# df['aggregg'] = df.rr.map(hierarchy[case][ralevel]) -# df['interface'] = df.aggreg + '|' + df.aggregg - -# df = ( -# df -# .loc[df.aggreg != df.aggregg] -# .groupby(['t',hcol,'interface']).GW.sum().unstack('interface').fillna(0) -# ) -# if df.empty: -# continue -# else: -# df = df.loc[lastyear].copy() -# ## Order interfaces alphabetically -# rename = {} -# for interface in df: -# r, rr = interface.split('|') -# if r > rr: -# rename[interface] = f'{rr}|{r}' -# df[interface] *= -1 -# df = df.rename(columns=rename).groupby(axis=1, level=0).sum() -# ## Now reorder interfaces by flow -# rename = {} -# for interface in df: -# r, rr = interface.split('|') -# if df[interface].clip(lower=0).sum() < df[interface].clip(upper=0).abs().sum(): -# rename[interface] = f'{rr}|{r}' -# df[interface] *= -1 -# tran_flow_stress_agg[case] = df.rename(columns=rename).copy() - -# ### Calculate regional imports/exports -# dfimportexport = {} -# for case in cases: -# df = {} -# for region in regions: -# df[region] = reedsplots.get_import_export( -# region=region, df=tran_flow_stress_agg[case] -# ) -# dfimportexport[case] = pd.concat(df).sum(axis=1).unstack(level=0) -# dfimportexport[case].columns = dfimportexport[case].columns.rename('region') - -# ### Plot it -# whiteout = dict(zip( -# [f'C{i}' for i in range(10)], -# [plt.cm.tab20(i*2+1) for i in range(10)] -# )) -# if any([v not in whiteout for v in list(colors.values())]): -# whiteout = {v: (v[0], v[1], v[2], v[3]*0.7) for v in list(colors.values())} - -# for label in ['max','average']: -# plt.close() -# f,ax = plt.subplots( -# nrows, ncols, figsize=(SLIDE_WIDTH, SLIDE_HEIGHT), -# gridspec_kw={'wspace':0.0,'hspace':-0.1}, -# ) - -# for case in cases: -# ### Formatting -# dfmap[ralevel].plot(ax=ax[coords[case]], facecolor='none', edgecolor='C7', lw=0.5) -# dfmap['interconnect'].plot(ax=ax[coords[case]], facecolor='none', edgecolor='k', lw=1) -# ax[coords[case]].set_title( -# case, y=0.95, weight='bold', color=colors[case], fontsize=14) -# ax[coords[case]].axis('off') - -# ### RA flows -# if case not in tran_flow_stress_agg: -# continue - -# if label == 'max': -# ## Max flow -# forwardwidth = tran_flow_stress_agg[case].clip(lower=0).max() -# reversewidth = abs(tran_flow_stress_agg[case].clip(upper=0).min()) -# else: -# ## Average when it's flowing -# forwardwidth = ( -# tran_flow_stress_agg[case].clip(lower=0).sum() -# / tran_flow_stress_agg[case].clip(lower=0).astype(bool).sum() -# ) -# reversewidth = ( -# tran_flow_stress_agg[case].clip(upper=0).abs().sum() -# / tran_flow_stress_agg[case].clip(upper=0).abs().astype(bool).sum() -# ) - -# interfaces = tran_flow_stress_agg[case].columns -# numdays = ( -# len(tran_flow_stress_agg[case]) -# * int(dictin_sw[case].GSw_HourlyChunkLengthStress) -# // 24 -# ) - -# ### Head/tail length: -# gwh_forward = tran_flow_stress_agg[case].clip(lower=0).sum() -# gwh_reverse = abs(tran_flow_stress_agg[case].clip(upper=0).sum()) - -# reversefrac = gwh_reverse / (gwh_reverse + gwh_forward) -# forwardfrac = gwh_forward / (gwh_reverse + gwh_forward) - -# ### Plot it -# for interface in interfaces: -# r, rr = interface.split('|') -# startx, starty = dfmap[ralevel].loc[r, ['x', 'y']] -# endx, endy = dfmap[ralevel].loc[rr, ['x', 'y']] - -# plots.plot_segmented_arrow( -# ax[coords[case]], -# reversefrac=reversefrac[interface], -# forwardfrac=forwardfrac[interface], -# reversewidth=reversewidth[interface]*wscale, -# forwardwidth=forwardwidth[interface]*wscale, -# startx=startx, starty=starty, endx=endx, endy=endy, -# forwardcolor=colors[case], reversecolor=whiteout[colors[case]], -# alpha=0.8, headwidthfrac=1.5, -# ) -# ### Scale -# if scale: -# (startx, starty, endx, endy) = (-2.0e6, -1.2e6, -1.5e6, -1.2e6) -# yspan = ax[coords[case]].get_ylim() -# yspan = yspan[1] - yspan[0] -# plots.plot_segmented_arrow( -# ax[coords[case]], -# reversefrac=0, forwardfrac=1, -# reversewidth=0, forwardwidth=scale*wscale, -# startx=startx, starty=starty, endx=endx, endy=endy, -# forwardcolor=colors[case], reversecolor=whiteout[colors[case]], -# alpha=0.8, headwidthfrac=1.5, -# ) -# ax[coords[case]].annotate( -# f"{scale} GW\n{label}", ((startx+endx)/2, starty-(scale/2*wscale)-yspan*0.02), -# ha='center', va='top', fontsize=14, -# ) - -# ### Save it -# title = f'{ralevel} {label} RA flows' -# slide = reeds.report_utils.add_to_pptx(title, prs=prs) -# if interactive: -# plt.show() -# except Exception: -# print(traceback.format_exc()) +#%% Flexibly sited load +try: + if any([float(dictin_sw[c].get('GSw_LoadSiteCF', 0)) for c in cases]): + loadsite_cases = { + k:v for k,v in cases.items() if float(dictin_sw[k].get('GSw_LoadSiteCF', 0)) + } + f, ax, dictplot = reeds.reedsplots.map_output_byyear( + case=loadsite_cases, + param='loadsite_cap', + years=[lastyear], + vscale=1e-3, + vmin=0, + title='Sited demand [GW]', + ) + ## Save it + slide = reeds.results.add_to_pptx('Flexibly Sited Demand', prs=prs, width=SLIDE_WIDTH) + if interactive: + plt.show() +except Exception: + print(traceback.format_exc()) + +#%%### RA sharing +if detailed: + try: + ralevel = 'transreg' + scale = 10 + wscale = 7e3 + dfmap = reeds.io.get_dfmap(cases[basecase]) + regions = dfmap[ralevel].bounds.minx.sort_values().index + + ### Calculate aggregated load + dictin_load_stress_agg = {} + for case in cases: + if int(dictin_sw[case].GSw_PRM_CapCredit): + hcol = 'ccseason' + df = dictin_peak_ccseason[case].copy() + else: + hcol = 'h' + df = dictin_load_stress[case].copy() + df = ( + df.assign(region=df.r.map(hierarchy[case][ralevel])) + .groupby(['t','region',hcol]).GW.sum() + .loc[lastyear].unstack('region') + ) + dictin_load_stress_agg[case] = df.sum() + + ### Calculate aggregated stress period flows + tran_flow_stress_agg = {} + for case in cases: + if int(dictin_sw[case].GSw_PRM_CapCredit): + df = dictin_prmtrade[case].copy() + hcol = 'ccseason' + else: + df = dictin_tran_flow_stress[case].copy() + hcol = 'h' + df['aggreg'] = df.r.map(hierarchy[case][ralevel]) + df['aggregg'] = df.rr.map(hierarchy[case][ralevel]) + df['interface'] = df.aggreg + '|' + df.aggregg + + df = ( + df + .loc[df.aggreg != df.aggregg] + .groupby(['t',hcol,'interface']).GW.sum().unstack('interface').fillna(0) + ) + if df.empty: + continue + else: + df = df.loc[lastyear].copy() + ## Order interfaces alphabetically + rename = {} + for interface in df: + r, rr = interface.split('|') + if r > rr: + rename[interface] = f'{rr}|{r}' + df[interface] *= -1 + df = df.rename(columns=rename).groupby(axis=1, level=0).sum() + ## Now reorder interfaces by flow + rename = {} + for interface in df: + r, rr = interface.split('|') + if df[interface].clip(lower=0).sum() < df[interface].clip(upper=0).abs().sum(): + rename[interface] = f'{rr}|{r}' + df[interface] *= -1 + tran_flow_stress_agg[case] = df.rename(columns=rename).copy() + + ### Calculate regional imports/exports + dfimportexport = {} + for case in cases: + df = {} + for region in regions: + df[region] = reedsplots.get_import_export( + region=region, df=tran_flow_stress_agg[case] + ) + dfimportexport[case] = pd.concat(df).sum(axis=1).unstack(level=0) + dfimportexport[case].columns = dfimportexport[case].columns.rename('region') + + ### Plot it + whiteout = dict(zip( + [f'C{i}' for i in range(10)], + [plt.cm.tab20(i*2+1) for i in range(10)] + )) + if any([v not in whiteout for v in list(colors.values())]): + whiteout = {v: (v[0], v[1], v[2], v[3]*0.7) for v in list(colors.values())} + + for label in ['max','average']: + plt.close() + f,ax = plt.subplots( + nrows, ncols, figsize=(SLIDE_WIDTH, SLIDE_HEIGHT), + gridspec_kw={'wspace':0.0,'hspace':-0.1}, + ) + + for case in cases: + ### Formatting + dfmap[ralevel].plot(ax=ax[coords[case]], facecolor='none', edgecolor='C7', lw=0.5) + dfmap['interconnect'].plot(ax=ax[coords[case]], facecolor='none', edgecolor='k', lw=1) + ax[coords[case]].set_title( + case, y=0.95, weight='bold', color=colors[case], fontsize=14) + ax[coords[case]].axis('off') + + ### RA flows + if case not in tran_flow_stress_agg: + continue + + if label == 'max': + ## Max flow + forwardwidth = tran_flow_stress_agg[case].clip(lower=0).max() + reversewidth = abs(tran_flow_stress_agg[case].clip(upper=0).min()) + else: + ## Average when it's flowing + forwardwidth = ( + tran_flow_stress_agg[case].clip(lower=0).sum() + / tran_flow_stress_agg[case].clip(lower=0).astype(bool).sum() + ) + reversewidth = ( + tran_flow_stress_agg[case].clip(upper=0).abs().sum() + / tran_flow_stress_agg[case].clip(upper=0).abs().astype(bool).sum() + ) + + interfaces = tran_flow_stress_agg[case].columns + numdays = ( + len(tran_flow_stress_agg[case]) + * int(dictin_sw[case].GSw_HourlyChunkLengthStress) + // 24 + ) + + ### Head/tail length: + gwh_forward = tran_flow_stress_agg[case].clip(lower=0).sum() + gwh_reverse = abs(tran_flow_stress_agg[case].clip(upper=0).sum()) + + reversefrac = gwh_reverse / (gwh_reverse + gwh_forward) + forwardfrac = gwh_forward / (gwh_reverse + gwh_forward) + + ### Plot it + for interface in interfaces: + r, rr = interface.split('|') + startx, starty = dfmap[ralevel].loc[r, ['x', 'y']] + endx, endy = dfmap[ralevel].loc[rr, ['x', 'y']] + + plots.plot_segmented_arrow( + ax[coords[case]], + reversefrac=reversefrac[interface], + forwardfrac=forwardfrac[interface], + reversewidth=reversewidth[interface]*wscale, + forwardwidth=forwardwidth[interface]*wscale, + startx=startx, starty=starty, endx=endx, endy=endy, + forwardcolor=colors[case], reversecolor=whiteout[colors[case]], + alpha=0.8, headwidthfrac=1.5, + ) + ### Scale + if scale: + (startx, starty, endx, endy) = (-2.0e6, -1.2e6, -1.5e6, -1.2e6) + yspan = ax[coords[case]].get_ylim() + yspan = yspan[1] - yspan[0] + plots.plot_segmented_arrow( + ax[coords[case]], + reversefrac=0, forwardfrac=1, + reversewidth=0, forwardwidth=scale*wscale, + startx=startx, starty=starty, endx=endx, endy=endy, + forwardcolor=colors[case], reversecolor=whiteout[colors[case]], + alpha=0.8, headwidthfrac=1.5, + ) + ax[coords[case]].annotate( + f"{scale} GW\n{label}", ((startx+endx)/2, starty-(scale/2*wscale)-yspan*0.02), + ha='center', va='top', fontsize=14, + ) + + ### Save it + title = f'{ralevel} {label} RA flows' + slide = reeds.report_utils.add_to_pptx(title, prs=prs) + if interactive: + plt.show() + except Exception: + print(traceback.format_exc()) #%%### Copy some premade single-case plots @@ -2538,220 +2538,220 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): print(f'No outputs/figures/{figname}.png for {os.path.basename(cases[case])}') -# #%%### Generation capacity maps - -# ### Shared data -# try: -# base = cases[list(cases.keys())[0]] -# val_r = dictin_cap_r[basecase].r.unique() -# dfmap = reeds.io.get_dfmap(base) -# dfba = dfmap['r'] -# dfstates = dfmap['st'] -# if (len(cases) == 2) and (not forcemulti): -# for i_plot in maptechs.keys(): -# plt.close() -# f,ax=plt.subplots( -# 1, 3, sharex=True, sharey=True, figsize=(14,8), -# gridspec_kw={'wspace':-0.05, 'hspace':0.05}, -# dpi=150, -# ) - -# _,_,dfplot = reedsplots.plot_diff_maps( -# val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], -# year=lastyear, casebase=casebase, casecomp=casecomp, -# plot='base', f=f, ax=ax[0], -# cmap=cmocean.cm.rain, -# ) -# ax[0].annotate( -# casebase_name, -# (0.1,1), xycoords='axes fraction', fontsize=10) - -# _,_,dfplot = reedsplots.plot_diff_maps( -# val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], -# year=lastyear, casebase=casebase, casecomp=casecomp, -# plot='comp', f=f, ax=ax[1], -# cmap=cmocean.cm.rain, -# ) -# ax[1].annotate( -# casecomp_name, -# (0.1,1), xycoords='axes fraction', fontsize=10) - -# _,_,dfplot = reedsplots.plot_diff_maps( -# val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], -# year=lastyear, casebase=casebase, casecomp=casecomp, -# plot='absdiff', f=f, ax=ax[2], -# cmap=plt.cm.RdBu_r, -# ) -# ax[2].annotate( -# '{}\n– {}'.format( -# casecomp_name, -# casebase_name), -# (0.1,1), xycoords='axes fraction', fontsize=10) - -# reeds.report_utils.add_to_pptx(f'{i_plot} Capacity {lastyear} [GW]', prs=prs) -# if interactive: -# plt.show() -# else: -# figwidth = SLIDE_WIDTH -# #### Absolute maps -# if (nrows == 1) or (ncols == 1): -# legendcoords = max(nrows, ncols) - 1 -# elif (nrows-1, ncols-1) in coords.values(): -# legendcoords = (nrows-1, ncols-1) -# else: -# legendcoords = (nrows-2, ncols-1) - -# ### Set up plot -# for tech_type in maptechs.keys(): -# ### Get limits -# vmin = 0. -# vmax = float(pd.concat({ -# case: dictin_cap_r[case].loc[ -# (dictin_cap_r[case].i.isin(maptechs[tech_type])) -# & (dictin_cap_r[case].t.astype(int)==lastyear) -# ].groupby('r').MW.sum() -# for case in cases -# }).max()) / 1e3 -# if np.isnan(vmax): -# vmax = 0. -# if not vmax: -# print(f'{tech_type} has zero capacity in {lastyear}, so skipping maps') -# continue -# ### Set up plot -# plt.close() -# f,ax = plt.subplots( -# nrows, ncols, figsize=(figwidth, SLIDE_HEIGHT), -# gridspec_kw={'wspace':0.0,'hspace':-0.1}, -# ) -# ### Plot it -# for case in cases: -# dfval = dictin_cap_r[case].loc[ -# (dictin_cap_r[case].i.isin(maptechs[tech_type])) -# & (dictin_cap_r[case].t.astype(int)==lastyear) -# ].groupby('r').MW.sum() -# dfplot = dfba.copy() -# dfplot['GW'] = (dfval / 1e3).fillna(0) - -# ax[coords[case]].set_title( -# case if nowrap else plots.wraptext(case, width=figwidth/ncols*0.9, fontsize=14) -# ) -# dfba.plot( -# ax=ax[coords[case]], -# facecolor='none', edgecolor='k', lw=0.1, zorder=10000) -# dfstates.plot( -# ax=ax[coords[case]], -# facecolor='none', edgecolor='k', lw=0.2, zorder=10001) -# dfplot.plot( -# ax=ax[coords[case]], column='GW', cmap=cmap, vmin=vmin, vmax=vmax, -# legend=False, -# ) -# ## Legend -# if coords[case] == legendcoords: -# plots.addcolorbarhist( -# f=f, ax0=ax[coords[case]], data=dfplot.GW.values, -# title=f'{tech_type} {lastyear}\ncapacity [GW]', cmap=cmap, vmin=vmin, vmax=vmax, -# orientation='horizontal', labelpad=2.25, histratio=0., -# cbarwidth=0.05, cbarheight=0.85, -# cbarbottom=-0.05, cbarhoffset=0., -# ) - -# for row in range(nrows): -# for col in range(ncols): -# if nrows == 1: -# ax[col].axis('off') -# elif ncols == 1: -# ax[row].axis('off') -# else: -# ax[row,col].axis('off') -# ### Save it -# slide = reeds.report_utils.add_to_pptx(f'{tech_type} Capacity {lastyear} [GW]', prs=prs) -# if interactive: -# plt.show() - -# #### Difference maps -# ### Set up plot -# for tech_type in maptechs.keys(): -# ### Get limits -# dfval = pd.concat({ -# case: dictin_cap_r[case].loc[ -# (dictin_cap_r[case].i.isin(maptechs[tech_type])) -# & (dictin_cap_r[case].t.astype(int)==lastyear) -# ].groupby('r').MW.sum() -# for case in cases -# }, axis=1).fillna(0) / 1e3 -# dfdiff = dfval.subtract(dfval[basecase], axis=0) -# ### Get colorbar limits -# absmax = dfval.stack().max() -# diffmax = dfdiff.unstack().abs().max() - -# if np.isnan(absmax): -# absmax = 0. -# if not absmax: -# print(f'{tech_type} has zero capacity in {lastyear}, so skipping maps') -# continue -# ### Set up plot -# plt.close() -# f,ax = plt.subplots( -# nrows, ncols, figsize=(figwidth, SLIDE_HEIGHT), -# gridspec_kw={'wspace':0.0,'hspace':-0.1}, -# ) -# ### Plot it -# for case in cases: -# dfplot = dfba.copy() -# dfplot['GW'] = dfval[case] if case == basecase else dfdiff[case] - -# ax[coords[case]].set_title( -# case if nowrap else plots.wraptext(case, width=figwidth/ncols*0.9, fontsize=14) -# ) -# dfba.plot( -# ax=ax[coords[case]], -# facecolor='none', edgecolor='k', lw=0.1, zorder=10000) -# dfstates.plot( -# ax=ax[coords[case]], -# facecolor='none', edgecolor='k', lw=0.2, zorder=10001) -# dfplot.plot( -# ax=ax[coords[case]], column='GW', -# cmap=(cmap if case == basecase else cmap_diff), -# vmin=(0 if case == basecase else -diffmax), -# vmax=(absmax if case == basecase else diffmax), -# legend=False, -# ) -# ## Difference legend -# if coords[case] == legendcoords: -# plots.addcolorbarhist( -# f=f, ax0=ax[coords[case]], data=dfplot.GW.values, -# title=f'{tech_type} {lastyear}\ncapacity, difference\nfrom {basecase} [GW]', -# cmap=(cmap if case == basecase else cmap_diff), -# vmin=(0 if case == basecase else -diffmax), -# vmax=(absmax if case == basecase else diffmax), -# orientation='horizontal', labelpad=2.25, histratio=0., -# cbarwidth=0.05, cbarheight=0.85, -# cbarbottom=-0.05, cbarhoffset=0., -# ) -# ## Absolute legend -# plots.addcolorbarhist( -# f=f, ax0=ax[coords[basecase]], data=dfval[basecase].values, -# title=f'{tech_type} {lastyear}\ncapacity [GW]', -# cmap=cmap, vmin=0, vmax=absmax, -# orientation='horizontal', labelpad=2.25, histratio=0., -# cbarwidth=0.05, cbarheight=0.85, -# cbarbottom=-0.05, cbarhoffset=0., -# ) - -# for row in range(nrows): -# for col in range(ncols): -# if nrows == 1: -# ax[col].axis('off') -# elif ncols == 1: -# ax[row].axis('off') -# else: -# ax[row,col].axis('off') -# ### Save it -# slide = reeds.report_utils.add_to_pptx(f'Difference: {tech_type} Capacity {lastyear} [GW]', prs=prs) -# if interactive: -# plt.show() -# except Exception: -# print(traceback.format_exc()) +#%%### Generation capacity maps + +### Shared data +try: + base = cases[list(cases.keys())[0]] + val_r = dictin_cap_r[basecase].r.unique() + dfmap = reeds.io.get_dfmap(base) + dfba = dfmap['r'] + dfstates = dfmap['st'] + if (len(cases) == 2) and (not forcemulti): + for i_plot in maptechs.keys(): + plt.close() + f,ax=plt.subplots( + 1, 3, sharex=True, sharey=True, figsize=(14,8), + gridspec_kw={'wspace':-0.05, 'hspace':0.05}, + dpi=150, + ) + + _,_,dfplot = reedsplots.plot_diff_maps( + val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], + year=lastyear, casebase=casebase, casecomp=casecomp, + plot='base', f=f, ax=ax[0], + cmap=cmocean.cm.rain, + ) + ax[0].annotate( + casebase_name, + (0.1,1), xycoords='axes fraction', fontsize=10) + + _,_,dfplot = reedsplots.plot_diff_maps( + val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], + year=lastyear, casebase=casebase, casecomp=casecomp, + plot='comp', f=f, ax=ax[1], + cmap=cmocean.cm.rain, + ) + ax[1].annotate( + casecomp_name, + (0.1,1), xycoords='axes fraction', fontsize=10) + + _,_,dfplot = reedsplots.plot_diff_maps( + val=mapdiff, i_plot=i_plot, titles = maptechs[i_plot], + year=lastyear, casebase=casebase, casecomp=casecomp, + plot='absdiff', f=f, ax=ax[2], + cmap=plt.cm.RdBu_r, + ) + ax[2].annotate( + '{}\n– {}'.format( + casecomp_name, + casebase_name), + (0.1,1), xycoords='axes fraction', fontsize=10) + + reeds.report_utils.add_to_pptx(f'{i_plot} Capacity {lastyear} [GW]', prs=prs) + if interactive: + plt.show() + else: + figwidth = SLIDE_WIDTH + #### Absolute maps + if (nrows == 1) or (ncols == 1): + legendcoords = max(nrows, ncols) - 1 + elif (nrows-1, ncols-1) in coords.values(): + legendcoords = (nrows-1, ncols-1) + else: + legendcoords = (nrows-2, ncols-1) + + ### Set up plot + for tech_type in maptechs.keys(): + ### Get limits + vmin = 0. + vmax = float(pd.concat({ + case: dictin_cap_r[case].loc[ + (dictin_cap_r[case].i.isin(maptechs[tech_type])) + & (dictin_cap_r[case].t.astype(int)==lastyear) + ].groupby('r').MW.sum() + for case in cases + }).max()) / 1e3 + if np.isnan(vmax): + vmax = 0. + if not vmax: + print(f'{tech_type} has zero capacity in {lastyear}, so skipping maps') + continue + ### Set up plot + plt.close() + f,ax = plt.subplots( + nrows, ncols, figsize=(figwidth, SLIDE_HEIGHT), + gridspec_kw={'wspace':0.0,'hspace':-0.1}, + ) + ### Plot it + for case in cases: + dfval = dictin_cap_r[case].loc[ + (dictin_cap_r[case].i.isin(maptechs[tech_type])) + & (dictin_cap_r[case].t.astype(int)==lastyear) + ].groupby('r').MW.sum() + dfplot = dfba.copy() + dfplot['GW'] = (dfval / 1e3).fillna(0) + + ax[coords[case]].set_title( + case if nowrap else plots.wraptext(case, width=figwidth/ncols*0.9, fontsize=14) + ) + dfba.plot( + ax=ax[coords[case]], + facecolor='none', edgecolor='k', lw=0.1, zorder=10000) + dfstates.plot( + ax=ax[coords[case]], + facecolor='none', edgecolor='k', lw=0.2, zorder=10001) + dfplot.plot( + ax=ax[coords[case]], column='GW', cmap=cmap, vmin=vmin, vmax=vmax, + legend=False, + ) + ## Legend + if coords[case] == legendcoords: + plots.addcolorbarhist( + f=f, ax0=ax[coords[case]], data=dfplot.GW.values, + title=f'{tech_type} {lastyear}\ncapacity [GW]', cmap=cmap, vmin=vmin, vmax=vmax, + orientation='horizontal', labelpad=2.25, histratio=0., + cbarwidth=0.05, cbarheight=0.85, + cbarbottom=-0.05, cbarhoffset=0., + ) + + for row in range(nrows): + for col in range(ncols): + if nrows == 1: + ax[col].axis('off') + elif ncols == 1: + ax[row].axis('off') + else: + ax[row,col].axis('off') + ### Save it + slide = reeds.report_utils.add_to_pptx(f'{tech_type} Capacity {lastyear} [GW]', prs=prs) + if interactive: + plt.show() + + #### Difference maps + ### Set up plot + for tech_type in maptechs.keys(): + ### Get limits + dfval = pd.concat({ + case: dictin_cap_r[case].loc[ + (dictin_cap_r[case].i.isin(maptechs[tech_type])) + & (dictin_cap_r[case].t.astype(int)==lastyear) + ].groupby('r').MW.sum() + for case in cases + }, axis=1).fillna(0) / 1e3 + dfdiff = dfval.subtract(dfval[basecase], axis=0) + ### Get colorbar limits + absmax = dfval.stack().max() + diffmax = dfdiff.unstack().abs().max() + + if np.isnan(absmax): + absmax = 0. + if not absmax: + print(f'{tech_type} has zero capacity in {lastyear}, so skipping maps') + continue + ### Set up plot + plt.close() + f,ax = plt.subplots( + nrows, ncols, figsize=(figwidth, SLIDE_HEIGHT), + gridspec_kw={'wspace':0.0,'hspace':-0.1}, + ) + ### Plot it + for case in cases: + dfplot = dfba.copy() + dfplot['GW'] = dfval[case] if case == basecase else dfdiff[case] + + ax[coords[case]].set_title( + case if nowrap else plots.wraptext(case, width=figwidth/ncols*0.9, fontsize=14) + ) + dfba.plot( + ax=ax[coords[case]], + facecolor='none', edgecolor='k', lw=0.1, zorder=10000) + dfstates.plot( + ax=ax[coords[case]], + facecolor='none', edgecolor='k', lw=0.2, zorder=10001) + dfplot.plot( + ax=ax[coords[case]], column='GW', + cmap=(cmap if case == basecase else cmap_diff), + vmin=(0 if case == basecase else -diffmax), + vmax=(absmax if case == basecase else diffmax), + legend=False, + ) + ## Difference legend + if coords[case] == legendcoords: + plots.addcolorbarhist( + f=f, ax0=ax[coords[case]], data=dfplot.GW.values, + title=f'{tech_type} {lastyear}\ncapacity, difference\nfrom {basecase} [GW]', + cmap=(cmap if case == basecase else cmap_diff), + vmin=(0 if case == basecase else -diffmax), + vmax=(absmax if case == basecase else diffmax), + orientation='horizontal', labelpad=2.25, histratio=0., + cbarwidth=0.05, cbarheight=0.85, + cbarbottom=-0.05, cbarhoffset=0., + ) + ## Absolute legend + plots.addcolorbarhist( + f=f, ax0=ax[coords[basecase]], data=dfval[basecase].values, + title=f'{tech_type} {lastyear}\ncapacity [GW]', + cmap=cmap, vmin=0, vmax=absmax, + orientation='horizontal', labelpad=2.25, histratio=0., + cbarwidth=0.05, cbarheight=0.85, + cbarbottom=-0.05, cbarhoffset=0., + ) + + for row in range(nrows): + for col in range(ncols): + if nrows == 1: + ax[col].axis('off') + elif ncols == 1: + ax[row].axis('off') + else: + ax[row,col].axis('off') + ### Save it + slide = reeds.report_utils.add_to_pptx(f'Difference: {tech_type} Capacity {lastyear} [GW]', prs=prs) + if interactive: + plt.show() +except Exception: + print(traceback.format_exc()) #%% Save the powerpoint file prs.save(savename) From 2bfecac5d8407ed6a23e56c534f09088f8f7bded Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Fri, 26 Jun 2026 12:50:38 -0600 Subject: [PATCH 23/24] reedsplots.py: fix scaling in plot_stressperiod_evolution(); bypass hull failures in plot_eue_events() --- reeds/reedsplots.py | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/reeds/reedsplots.py b/reeds/reedsplots.py index 816fb19b..2f1b7b9c 100644 --- a/reeds/reedsplots.py +++ b/reeds/reedsplots.py @@ -4131,12 +4131,22 @@ def plot_stressperiod_evolution( Only the first /-delimited level_threshold pair is used. """ from reeds.resource_adequacy import stress_periods + ### Fixed inputs + ylabel = { + 'neue': 'NEUE [ppm]', + 'lolh': 'LOLH [event-hours/year]', + 'lole': 'LOLE [events/year]', + 'lold': 'LOLD [event-days/year]', + 'duration': 'Duration [hours]', + 'depth': 'Depth [% of peak]', + } + scale = {'depth':100} ### Parse inputs sw = reeds.io.get_switches(case) switch = stress_periods.RA_SWITCHES[metric.lower()] switch_metric = stress_periods.SWITCH_METRIC[metric.lower()] level, threshold = sw[switch].split('/')[0].split('_') - threshold = float(threshold) + threshold = float(threshold) * scale.get(metric, 1) ### Load RA results infiles = sorted(glob(os.path.join(case,'outputs','ra_metrics_*.csv'))) dictin_ra = { @@ -4150,19 +4160,10 @@ def plot_stressperiod_evolution( .xs(level,0,'level') .xs(switch_metric,0,'metric') .squeeze(1).unstack('region') - ) + ) * scale.get(metric, 1) ### Load stress periods for labels dfstress = get_stressperiods(case) ### Plot setup - ylabel = { - 'neue': 'NEUE [ppm]', - 'lolh': 'LOLH [event-hours/year]', - 'lole': 'LOLE [events/year]', - 'lold': 'LOLD [event-days/year]', - 'duration': 'Duration [hours]', - 'depth': 'Depth [% of peak]', - } - scale = {'depth':100} years = [ y for y in dfplot.index.get_level_values('year').unique() if y >= int(sw.GSw_StartMarkets) @@ -4181,7 +4182,7 @@ def plot_stressperiod_evolution( gridspec_kw=gridspec_kw, ) for col, year in enumerate(years): - df = dfplot.loc[year] * scale.get(metric, 1) + df = dfplot.loc[year] for region in regions: ax[col].plot( df.index, df[region], @@ -4213,7 +4214,7 @@ def plot_stressperiod_evolution( handletextpad=0.3, handlelength=0.7, ) ax[0].set_ylabel(ylabel[metric]) - ax[0].set_ylim(0) + ax[0].set_ylim(0, min(ax[0].get_ylim()[1], threshold*50)) plots.despine(ax) return f,ax @@ -6616,9 +6617,12 @@ def plot_eue_events( ) if showhull and len(points['x']) >= 3: dfpoints = pd.DataFrame(points) - hull = scipy.spatial.ConvexHull(dfpoints) - dfhull = dfpoints.loc[hull.vertices] - _ax.fill(dfhull.x, dfhull.y, color='C3', lw=0, alpha=0.3, zorder=-1) + try: + hull = scipy.spatial.ConvexHull(dfpoints) + dfhull = dfpoints.loc[hull.vertices] + _ax.fill(dfhull.x, dfhull.y, color='C3', lw=0, alpha=0.3, zorder=-1) + except Exception as err: + print(err) ## Formatting if iteration == 0: _ax.set_title(region.replace('_','\n'), weight='bold') From 3eaf1fd46693d944a40bc8de35f7051e34b8e848 Mon Sep 17 00:00:00 2001 From: Patrick Brown <25125211+patrickbrown4@users.noreply.github.com> Date: Sat, 27 Jun 2026 09:12:59 -0600 Subject: [PATCH 24/24] fix get_longest_events() when there are no events --- reeds/resource_adequacy/stress_periods.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 908e4291..fe6bbd23 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -273,11 +273,15 @@ def get_longest_events( pd.Series(index=pd.date_range(row.start, row.end, freq='H'), data=1) .resample('D').count() ) - metric_period = pd.concat(dates) - metric_period.index = metric_period.index.map( - lambda x: reeds.timeseries.timestamp2h(x, sw.GSw_HourlyType).split('h')[0] - ) - return metric_period.groupby(level=0).sum() + if len(dates): + metric_period = pd.concat(dates) + metric_period.index = metric_period.index.map( + lambda x: reeds.timeseries.timestamp2h(x, sw.GSw_HourlyType).split('h')[0] + ) + metric_period = metric_period.groupby(level=0).sum() + else: + metric_period = pd.Series(dtype=int) + return metric_period def get_shoulder_periods(sw, criterion, dfenergy_agg, high_stress_periods):