diff --git a/cases.csv b/cases.csv index e813a404..9ea18ea7 100644 --- a/cases.csv +++ b/cases.csv @@ -261,14 +261,17 @@ GSw_PRM_NetImportLimit,Turn on (1) or off (0) the eq_firm_transfer_limit constra GSw_PRM_NetImportLimitScen,/-delimited list of {year}_{max percent import OR 'hist' for historical percent OR 'histmax' for max historical percent across all regions} values to use in eq_firm_transfer_limit,N/A,2031_hist/2050_100, GSw_PRM_scenario,"column of inputs/reserves/prm_annual.csv from which to draw the PRM (""nerc"" means the annual values from the 2024 NERC LTRA; ""static"" means the 2024 values from the same source); or a float (like 0.12) which gets applied as a fraction to all regions in all years (0.12 means 12%)",(static|nerc|ramp2025_20by50|^\d*\.?\d+$),0.12, GSw_PRM_StressIncrement,How many stress periods to add per iteration,int,2, -GSw_PRM_StressIterateMax,Max number of times to iterate on a given solve year to achieve GSw_PRM_StressThreshold when using stress periods (0 means don't iterate; 1 means 1 extra ReEDS run; etc),int,5, +GSw_PRM_StressIterateMax,"Max number of times to iterate on a given solve year to achieve adequacy stress levels (e.g., GSw_PRM_StressThresholdNEUE) when using stress periods (0 means don't iterate; 1 means 1 extra ReEDS run; etc)",int,5, GSw_PRM_StressLoadAggMethod,How to aggregate load for stress periods within the chunks specified by GSw_HourlyChunkLengthStress (only used if GSw_HourlyChunkAggMethod=mean; otherwise GSw_HourlyChunkAggMethod is used),mean; max,max, GSw_PRM_StressModel,Model used to identify stress periods: pras or a string starting with 'user' which specifies a file at inputs/temporal/stressperiods_{GSw_PRM_StressModel}.csv,N/A,pras, 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_StressThreshold,/-delimited list of annual NEUE level [ppm] above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_NEUEppm_StressMetric_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; NEUEppm is normalized expected unserved energy in parts per million; StressMetric is EUE or NEUE (only used in period selection); 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_EUE_sum, +GSw_PRM_StressThresholdMetrics,/-delimited list of metrics for identifying stress periods,N/A,NEUE, +GSw_PRM_StressThresholdLOLE,LOLE threshold above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_nLOLE_StressMetric_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; nLOLE is number of loss of load events; StressMetric is LOLE; 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_36_LOLE_sum, +GSw_PRM_StressThresholdEUE,/Annual EUE level [MWh] threshold above which to re-solve the latest model year with new stress periods; formulated as HierarchyLevel_EUE_StressMetric_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; EUE is expected unserved energy; StressMetric EUE (only used in period selection); 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_10000_EUE_sum, +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_StressMetric_PeriodAggMethod where HierarchyLevel is a column in hierarchy.csv; NEUEppm is normalized expected unserved energy in parts per million; StressMetric is EUE or NEUE (only used in period selection); 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_NEUE_sum, 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/postprocessing/compare_cases.py b/postprocessing/compare_cases.py index 57e21efb..722eed96 100644 --- a/postprocessing/compare_cases.py +++ b/postprocessing/compare_cases.py @@ -1436,7 +1436,7 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): #%%### SCOE, NEUE try: - neue_threshold = float(sw.GSw_PRM_StressThreshold.split('_')[1]) + neue_threshold = float(sw.GSw_PRM_StressThresholdNEUE.split('_')[1]) width = 9 + len(cases)*0.5 plt.close() f,ax = plt.subplots( @@ -1549,9 +1549,9 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): #%% Regional NEUE try: - neue_threshold = float(sw.GSw_PRM_StressThreshold.split('_')[1]) + neue_threshold = float(sw.GSw_PRM_StressThresholdNEUE.split('_')[1]) dfmap = reeds.io.get_dfmap(cases[basecase]) - level = sw.GSw_PRM_StressThreshold.split('_')[0] + level = sw.GSw_PRM_StressThresholdNEUE.split('_')[0] regions = dfmap[level].loc[hierarchy[basecase][level].unique()].bounds.minx.sort_values().index _nrows, _ncols, _coords = reeds.plots.get_coordinates(regions, ncols=6) labelcoords = { @@ -2500,7 +2500,10 @@ def two_bars(dfplot, basecase, colors, ax, col=0, ypad=0.02): #%%### Copy some premade single-case plots -level = dictin_sw[basecase]['GSw_PRM_StressThreshold'].split('_')[0] +# 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 metrics = [ 'cap', diff --git a/postprocessing/single_case_plots.py b/postprocessing/single_case_plots.py index 79a31000..58b4dfec 100644 --- a/postprocessing/single_case_plots.py +++ b/postprocessing/single_case_plots.py @@ -701,7 +701,9 @@ print(traceback.format_exc()) try: - level, threshold, _, metric = sw['GSw_PRM_StressThreshold'].split('/')[0].split('_') + _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] plt.close() f,ax = reedsplots.plot_stressperiod_evolution( case=case, level=level, metric=metric) diff --git a/reeds/reedsplots.py b/reeds/reedsplots.py index dae1aa92..851bacca 100644 --- a/reeds/reedsplots.py +++ b/reeds/reedsplots.py @@ -484,7 +484,7 @@ def plot_diff( ymax = max(ax[0].get_ylim()[1], ax[1].get_ylim()[1]) ymin = min(ax[0].get_ylim()[0], ax[1].get_ylim()[0]) if val == 'NEUE (ppm)': - neue_threshold = float(sw.GSw_PRM_StressThreshold.split('_')[1]) + neue_threshold = float(sw.GSw_PRM_StressThresholdNEUE.split('_')[1]) ymax = max(ymax, 10, neue_threshold*1.05) ax[0].axhline(neue_threshold, c='C7', ls='--', lw=0.75) ax[1].axhline(neue_threshold, c='C7', ls='--', lw=0.75) @@ -4137,7 +4137,9 @@ def plot_stressperiod_evolution( """Plot NEUE by year and stress period iteration""" ### Parse inputs sw = reeds.io.get_switches(case) - _level, _threshold, _, _metric = sw['GSw_PRM_StressThreshold'].split('/')[0].split('_') + _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 = _level if level is None else level threshold = float(_threshold) if threshold is None else threshold metric = _metric if metric is None else metric @@ -4261,8 +4263,9 @@ def plot_neue_bylevel( norm = {'sum':1, 'max':1e-4} ylabel = {'sum': 'Sum of NEUE [ppm]', 'max':'Max NEUE [%]'} thresholds = { - i.split('_')[0]: float(i.split('_')[1]) - for i in sw.GSw_PRM_StressThreshold.split('/') + 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('/') } ### Plot it plt.close() @@ -4326,8 +4329,8 @@ def map_neue( neue = reeds.io.read_output(case, f'neue_{year}i{_iteration}.csv') neue = neue.loc[neue.metric==metric].set_index(['level','region']).NEUE_ppm sw = reeds.io.get_switches(case) - neue_threshold = float(sw.GSw_PRM_StressThreshold.split('_')[1]) - neue_threshold_level = sw.GSw_PRM_StressThreshold.split('_')[0] + neue_threshold = float(sw.GSw_PRM_StressThresholdNEUE.split('_')[1]) + neue_threshold_level = sw.GSw_PRM_StressThresholdNEUE.split('_')[0] ### Set up plot levels = ['interconnect','nercr','transreg','transgrp','st','r'] @@ -6195,7 +6198,8 @@ def map_stressors( dfmap[k] = v.to_crs(crs) ### Derived inputs - criterion = sw.GSw_PRM_StressThreshold.split('/')[0] + _first_metric = sw.GSw_PRM_StressThresholdMetrics.split('/')[0].upper() + criterion = sw[f'GSw_PRM_StressThreshold{_first_metric}'] level = criterion.split('_')[0] regions = hierarchy[level].unique() region2rs = { @@ -6656,7 +6660,8 @@ def map_prm(case, tmin=2023, cmap=cmocean.cm.rain, scale=3, fontsize=7, vmax=Non ### Plot it nrows, ncols, coords = plots.get_coordinates(year2iteration.index, aspect=1.8) - level = sw.GSw_PRM_StressThreshold.split('_')[0] + _first_metric = sw.GSw_PRM_StressThresholdMetrics.split('/')[0].upper() + level = sw[f'GSw_PRM_StressThreshold{_first_metric}'].split('_')[0] plt.close() f,ax = plt.subplots( nrows, ncols, figsize=(ncols*scale, nrows*scale*0.8), sharex=True, sharey=True, diff --git a/reeds/resource_adequacy/run_pras.jl b/reeds/resource_adequacy/run_pras.jl index 10f59e0c..ec516737 100644 --- a/reeds/resource_adequacy/run_pras.jl +++ b/reeds/resource_adequacy/run_pras.jl @@ -204,9 +204,9 @@ function run_pras(pras_system_path::String, args::Dict) results = Dict{String,Any}(zip(keys(resultspec), results_tuple)) #%% Print some results for the entire modeled region to show it worked - @info "$(PRAS.LOLE(results["short"])) event-h" - @info "$(PRAS.EUE(results["short"])) MWh" - @info "NEUE = $(1e6 * PRAS.EUE(results["short"]).eue.estimate / sum(sys.regions.load)) ppm" + @info "$(PRAS.LOLE(results["short"]))" + @info "$(PRAS.EUE(results["short"]))" + @info "$(PRAS.NEUE(results["short"]))" ## Filter out DC regions used for VSC HVDC transmission regions = [r for r in sys.regions.names if !(occursin("|", r))] diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 8d159deb..e51c5f7f 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -18,10 +18,15 @@ #%%### Functions -def get_pras_eue(case, t, iteration=0): +def get_pras_stress_metric(case, t, iteration=0, stress_metric='EUE'): """ """ ### Get PRAS outputs + # NEUE load division takes place in the caller function + # get_stress_metric_periods() or get_annual_stress_metric() based on agg_period + if stress_metric.upper() == 'NEUE': + stress_metric = 'EUE' + dfpras = reeds.io.read_pras_results( os.path.join(case, 'handoff', 'PRAS', f"PRAS_{t}i{iteration}.h5") ) @@ -29,20 +34,19 @@ def get_pras_eue(case, t, iteration=0): sw = reeds.io.get_switches(case) dfpras.index = reeds.timeseries.get_timeindex(sw['resource_adequacy_years']) - ### Keep the EUE columns by zone - eue_tail = '_EUE' - dfeue = dfpras[[ + ### Keep the metric columns by zone + metric_tail = '_' + stress_metric.upper() + dfmetric = dfpras[[ c for c in dfpras - if (c.endswith(eue_tail) and not c.startswith('USA')) + if (c.endswith(metric_tail) and not c.startswith('USA')) ]].copy() - ## Drop the tailing _EUE - dfeue = dfeue.rename( - columns=dict(zip(dfeue.columns, [c[:-len(eue_tail)] for c in dfeue]))) - - return dfeue + ## Drop the tailing metric tail + dfmetric = dfmetric.rename( + columns=dict(zip(dfmetric.columns, [c[:-len(metric_tail)] for c in dfmetric]))) + return dfmetric -def get_eue_periods( +def get_stress_metric_periods( case, t, iteration=0, hierarchy_level='transgrp', stress_metric='EUE', @@ -70,37 +74,37 @@ def get_eue_periods( ### Get the region aggregator rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) - ### Get EUE from PRAS - dfeue = get_pras_eue(case=case, t=t, iteration=iteration) + ### Get stress metric from PRAS + dfmetric = get_pras_stress_metric(case=case, t=t, iteration=iteration, stress_metric=stress_metric) ## Aggregate to hierarchy_level - dfeue = ( - dfeue + 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 - if stress_metric.upper() == 'EUE': + if stress_metric.upper() != 'NEUE': ### Aggregate according to period_agg_method dfmetric_period = ( - dfeue - .groupby([dfeue.index.year, dfeue.index.month, dfeue.index.day]) + dfmetric + .groupby([dfmetric.index.year, dfmetric.index.month, dfmetric.index.day]) .agg(period_agg_method) .rename_axis(['y','m','d']) ) - elif stress_metric.upper() == 'NEUE': + else: ### Get load at hierarchy_level dfload = reeds.io.read_h5py_file( os.path.join( case,'handoff','reeds_data',f'pras_load_{t}.h5') ).rename(columns=rmap).groupby(level=0, axis=1).sum() - dfload.index = dfeue.index + dfload.index = dfmetric.index ### Recalculate NEUE [ppm] and aggregate appropriately if period_agg_method == 'sum': dfmetric_period = ( - dfeue - .groupby([dfeue.index.year, dfeue.index.month, dfeue.index.day]) + dfmetric + .groupby([dfmetric.index.year, dfmetric.index.month, dfmetric.index.day]) .agg(period_agg_method) .rename_axis(['y','m','d']) ) / ( @@ -111,8 +115,8 @@ def get_eue_periods( ) * 1e6 elif period_agg_method == 'max': dfmetric_period = ( - (dfeue / dfload) - .groupby([dfeue.index.year, dfeue.index.month, dfeue.index.day]) + (dfmetric / dfload) + .groupby([dfmetric.index.year, dfmetric.index.month, dfmetric.index.day]) .agg(period_agg_method) .rename_axis(['y','m','d']) ) * 1e6 @@ -135,10 +139,10 @@ def get_eue_periods( return dfmetric_top -def plot_eue_diagnostics(sw, t, iteration, high_eue_periods): +def plot_stress_diagnostics(sw, t, iteration, high_stress_periods): try: dates = ( - pd.concat(high_eue_periods) + pd.concat(high_stress_periods) .reset_index().actual_period.map(reeds.timeseries.h2timestamp) .dt.strftime('%Y-%m-%d') .tolist() @@ -200,41 +204,47 @@ def get_and_write_neue(sw, write=True): eue.to_csv(os.path.join(sw['casedir'],'outputs','eue.csv')) return neue - -def get_annual_neue(case, t, iteration=0): +def get_annual_stress_metric(case, t, stress_metric, iteration=0): """ """ - ### Get EUE from PRAS - dfeue = get_pras_eue(case=case, t=t, iteration=iteration) + ### Get values from PRAS + dfmetric = get_pras_stress_metric(case=case, t=t, iteration=iteration, stress_metric=stress_metric) ### Get load (for calculating NEUE) - dfload = reeds.io.read_h5py_file( - os.path.join( - case,'handoff','reeds_data',f'pras_load_{t}.h5') - ) - dfload.index = dfeue.index + if stress_metric.upper() == 'NEUE': + 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'] - _neue = {} + _metric = {} for hierarchy_level in levels: ### Get the region aggregator rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) - ### Get NEUE summed over year - _neue[hierarchy_level,'sum'] = ( - dfeue.rename(columns=rmap).groupby(axis=1, level=0).sum().sum() + + ### Get stress metric summed over year + _metric[hierarchy_level,'sum'] = dfmetric.rename(columns=rmap).groupby(axis=1, level=0).sum().sum() + + ### Get max stress metric hour + _metric[hierarchy_level,'max'] = dfmetric.rename(columns=rmap).groupby(axis=1, level=0).sum().max() + + 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 - _neue[hierarchy_level,'max'] = ( - dfeue.rename(columns=rmap).groupby(axis=1, level=0).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 + ).max() * 1e6 ### Combine it - neue = pd.concat(_neue, names=['level','metric','region']).rename('NEUE_ppm') - - return neue + metric = pd.concat(_metric, names=['level','metric','region']).rename(f'{stress_metric}') + return metric def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods): ## Stop if not needed @@ -256,7 +266,7 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods): 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, ppm, stress_metric, period_agg_method) = criterion.split('_') + (hierarchy_level, stress_level, stress_metric, period_agg_method) = criterion.split('_') ## Aggregate storage energy to hierarchy_level dfenergy_agg = ( @@ -273,13 +283,9 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods): day = pd.Timestamp('-'.join(row[['y','m','d']].astype(str).tolist())) - 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_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] - day_eue = high_eue_periods[criterion, f'high_{stress_metric}'].loc[i,'EUE'] day_index = np.where( timeindex == dfenergy_agg.loc[day.strftime('%Y-%m-%d')].iloc[0].name )[0][0] @@ -288,8 +294,7 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods): day_after = timeindex[(day_index + periodhours) % len(timeindex)] if ( - ((cutofftype == 'eue') and (end_headspace_MWh / day_eue >= float(cutoff))) - or ((cutofftype[:3] == 'cap') and (end_headspace_frac >= float(cutoff))) + (cutofftype[:3] == 'cap') and (end_headspace_frac >= float(cutoff)) or (cutofftype[:3] == 'abs') ): shoulder_periods[criterion, f'after_{row.name}'] = pd.Series({ @@ -299,8 +304,7 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods): print(f"Added {day_after} as shoulder stress period after {day}") if ( - ((cutofftype == 'eue') and (start_headspace_MWh / day_eue >= float(cutoff))) - or ((cutofftype[:3] == 'cap') and (start_headspace_frac >= float(cutoff))) + (cutofftype[:3] == 'cap') and (start_headspace_frac >= float(cutoff)) or (cutofftype[:3] == 'abs') ): shoulder_periods[criterion, f'before_{row.name}'] = pd.Series({ @@ -311,8 +315,84 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods): return shoulder_periods +def _evaluate_stress_threshold_criterion(stress_criteria, criterion, sw, t, iteration, dfenergy_r, stressperiods_this_iteration): + _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'] + + ## Example: criterion = 'transgrp_10_EUE_sum' + (hierarchy_level, stress_level, stress_metric, period_agg_method) = criterion.split('_') + + ### Get stored stress metric + stress_vals = pd.read_csv( + os.path.join(sw.casedir, 'outputs', f'{stress_metric}_{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') + ) -def get_eue_sorted_periods(sw, t, iteration): + ### 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_level)).any(): + _failed[criterion] = this_test.loc[this_test > float(stress_level)] + 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)) + ) + for period, row in _high_stress_periods[criterion, f'high_{stress_metric}'].iterrows(): + print( + f"Added {period} " + f"({reeds.timeseries.h2timestamp(period).strftime('%Y-%m-%d')}) " + f"as stress period for {row.r} " + f"({stress_metric} = {row[stress_metric]})" + ) + + ### Include "shoulder periods" before or after each period + ### if the storage state of charge is low + _shoulder_periods = { + **_shoulder_periods, + **get_shoulder_periods(sw, criterion, dfenergy_r, _high_stress_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 dfenergy = reeds.io.read_pras_results( os.path.join(sw['casedir'], 'handoff', 'PRAS', f"PRAS_{t}i{iteration}-energy.h5") @@ -326,12 +406,6 @@ def get_eue_sorted_periods(sw, t, iteration): .groupby(axis=1, level=0).sum() ) - ### Get NEUE - neue = pd.read_csv( - os.path.join(sw.casedir, 'outputs', f'neue_{t}i{iteration}.csv'), - index_col=['level', 'metric', 'region'], - ).squeeze(1) - ### Load this year's stress periods so we don't duplicate stressperiods_this_iteration = pd.read_csv( os.path.join( @@ -339,78 +413,42 @@ def get_eue_sorted_periods(sw, t, iteration): ) ### Check all stress criteria; for regions that fail, add new stress periods - _eue_sorted_periods = {} - failed = {} - high_eue_periods = {} - shoulder_periods = {} - for criterion in sw.GSw_PRM_StressThreshold.split('/'): - ## Example: criterion = 'transgrp_10_EUE_sum' - (hierarchy_level, ppm, stress_metric, period_agg_method) = criterion.split('_') - - eue_periods = get_eue_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 - _eue_sorted_periods[criterion] = ( - eue_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 = neue[hierarchy_level][period_agg_method] - - if (this_test > float(ppm)).any(): - failed[criterion] = this_test.loc[this_test > float(ppm)] - print(f"GSw_PRM_StressThreshold = {criterion} failed for:") - print(failed[criterion]) - ###### Add GSw_PRM_StressIncrement periods to the list for the next iteration - high_eue_periods[criterion, f'high_{stress_metric}'] = ( - _eue_sorted_periods[criterion].loc[ - ## Only include new stress periods for the region(s) that failed - _eue_sorted_periods[criterion].r.isin(failed[criterion].index) - ## Don't repeat existing stress periods - & ~(_eue_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)) - ) - for period, row in high_eue_periods[criterion, f'high_{stress_metric}'].iterrows(): - print( - f"Added {period} " - f"({reeds.timeseries.h2timestamp(period).strftime('%Y-%m-%d')}) " - f"as stress period for {row.r} " - f"({stress_metric} = {row[stress_metric]})" - ) - - ### Include "shoulder periods" before or after each period - ### if the storage state of charge is low - shoulder_periods = { - **shoulder_periods, - **get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods) - } - - ### Dealing with earlier criteria may also address later criteria, so stop here - break - - else: - print(f"GSw_PRM_StressThreshold = {criterion} passed") - - eue_sorted_periods = pd.concat(_eue_sorted_periods, names=['criterion']) + 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_metrics_units = {'EUE':'MWh', 'NEUE':'ppm', 'LOLE':'days'} + stress_metrics_col_names = {m:f'{m}_{stress_metrics_units[m]}' for m in + stress_metric_switches} + + for metric in stress_metric_switches: + for criterion in sw[f'GSw_PRM_StressThreshold{metric.upper()}'].split('/'): + print(f"Evaluating GSw_PRM_StressThreshold {metric.upper()} with criterion: {criterion}") + stress_criteria = _evaluate_stress_threshold_criterion(stress_criteria, + criterion, + sw, + t, + iteration, + dfenergy_r, + stressperiods_this_iteration, + ) + + 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']) ### 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_eue_periods, **shoulder_periods}, names=['criterion','periodtype'], + {**high_stress_periods, **shoulder_periods}, names=['criterion','periodtype'], ).reset_index().drop_duplicates(subset='actual_period', keep='first') else: return failed, None, None @@ -445,18 +483,17 @@ def get_eue_sorted_periods(sw, t, iteration): combined_periods_write.to_csv(outpath, index=False) ### Tables and plots for debugging - eue_sorted_periods.round(2).rename(columns={'EUE':'EUE_MWh','NEUE':'NEUE_ppm'}).to_csv( - os.path.join(sw.casedir, 'inputs_case', newstresspath, 'eue_sorted_periods.csv') + stress_sorted_periods.round(2).rename(columns=stress_metrics_col_names).to_csv( + os.path.join(sw.casedir, 'inputs_case', newstresspath, 'stress_metrics_sorted_periods.csv') ) - new_stress_periods.round(2).rename(columns={'EUE':'EUE_MWh','NEUE':'NEUE_ppm'}).to_csv( + new_stress_periods.round(2).rename(columns=stress_metrics_col_names).to_csv( os.path.join(sw.casedir, 'inputs_case', newstresspath, 'new_stress_periods.csv'), index=False, ) - plot_eue_diagnostics(sw, t, iteration, high_eue_periods) + plot_stress_diagnostics(sw, t, iteration, high_stress_periods) return failed, new_stressperiods_write, combined_periods_write - def prm_increment_pras(sw, t, iteration, combined_periods_write, failed_regions): try: hmap = pd.read_csv( @@ -579,7 +616,7 @@ def prm_increment_pras(sw, t, iteration, combined_periods_write, failed_regions) return prm_increment -def update_prm(sw, t, iteration, failed, combined_periods_write): +def update_prm(sw, t, iteration, failed, combined_periods_write, stress_metrics): """Update the energy reserve margin by region r for stress periods, either using a static increment (GSw_PRM_UpdateMethod=1) or based on the estimated surplus needed by PRAS to recover the desired reliabiliaty criteria (GSw_PRM_UpdateMethod>1). @@ -590,7 +627,9 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): iteration (int): ReEDS-PRAS iteration failed (dict): Dictionary of regions with unserved energy at the hierarchy_level and their criterion evaluations - stress_hours (pd.DataFrame): data frame of stress periods + combined_periods_write (pd.DataFrame): Data frame of combined stress periods + stress_metrics (list): List of stress metrics. If EUE is not included, + we can only use Fixed-increment update Returns: pd.DataFrame: Table of prm levels for the next PRAS iteration @@ -616,7 +655,13 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): ) ## Fixed-increment update - if int(sw.GSw_PRM_UpdateMethod) == 1: + if int(sw.GSw_PRM_UpdateMethod) == 1 or 'EUE' not in stress_metrics: + if int(sw.GSw_PRM_UpdateMethod) > 1: + print( + 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." + ) prm_increment = failed_regions.copy() prm_increment['fraction'] = float(sw['GSw_PRM_UpdateFraction']) ## PRAS-informed PRM update @@ -647,13 +692,18 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): def main(sw, t, iteration=0, logging=True): """ """ - #%% Write consolidated NEUE so far + #%% Write consolidated stress metrics so far try: - _neue_simple = get_and_write_neue(sw, write=True) - neue = get_annual_neue(sw.casedir, t, iteration=iteration) - neue.round(2).to_csv( - os.path.join(sw.casedir, 'outputs', f"neue_{t}i{iteration}.csv") - ) + ## TODO: Check whether to keep _neue_simple or not. + # _neue_simple = get_and_write_neue(sw, write=True) + # Check if EUE is added, if not add it since it's needed for selection of shoulder stress periods + stress_metrics = sw.GSw_PRM_StressThresholdMetrics.split('/') + for stress_metric in stress_metrics: + print(f"Calculating and writing annual {stress_metric} for iteration {iteration}") + dfmetric = get_annual_stress_metric(sw.casedir, t, stress_metric, iteration=iteration) + dfmetric.round(2).to_csv( + os.path.join(sw.casedir, 'outputs', f"{stress_metric}_{t}i{iteration}.csv") + ) except Exception as err: if int(sw['pras']) == 2: @@ -666,7 +716,7 @@ def main(sw, t, iteration=0, logging=True): return #%% Identify and write new stress periods - failed, new_stressperiods_write, combined_periods_write = get_eue_sorted_periods( + failed, new_stressperiods_write, combined_periods_write = get_stress_metrics_sorted_periods( sw=sw, t=t, iteration=iteration, ) @@ -699,7 +749,7 @@ def main(sw, t, iteration=0, logging=True): index_col='*r', ) else: - prm_next_iteration = update_prm(sw, t, iteration, failed, combined_periods_write) + prm_next_iteration = update_prm(sw, t, iteration, failed, combined_periods_write, stress_metrics) prm_next_iteration.to_csv( os.path.join(sw.casedir, 'inputs_case', newstresspath, 'prm.csv'), diff --git a/runreeds.py b/runreeds.py index 6f15a290..3a88e531 100644 --- a/runreeds.py +++ b/runreeds.py @@ -319,33 +319,53 @@ def check_compatibility(sw): reeds_path = os.path.dirname(__file__) hierarchy = reeds.io.get_hierarchy(GSw_ZoneSet=sw['GSw_ZoneSet']).reset_index() - for threshold in sw['GSw_PRM_StressThreshold'].split('/'): - ## Example: threshold = 'transgrp_10_EUE_sum' - allowed_levels = ['country','interconnect','nercr','transreg','transgrp','st','r'] - (hierarchy_level, ppm, stress_metric, period_agg_method) = threshold.split('_') - if hierarchy_level not in allowed_levels: - raise ValueError( - f"GSw_PRM_StressThreshold: level={hierarchy_level} but must be in:\n" - + '\n'.join(allowed_levels) - ) - if period_agg_method.lower() not in ['sum','max']: - raise ValueError("Fix period agg method in GSw_PRM_StressThreshold") - if not (float(ppm) >= 0): - raise ValueError( - "ppm in GSw_PRM_StressThreshold must be a positive number " - f"but '{ppm}' was provided" - ) - if stress_metric.upper() not in ['EUE','NEUE']: - raise ValueError( - "stress metric in GSw_PRM_StressThreshold must be 'EUE' or 'NEUE' " - f"but '{stress_metric}' was provided" - ) - if (sw['GSw_PRM_StressModel'].lower() != 'pras') and (stress_metric.upper() != 'EUE'): - err = ( - f"The combination of GSw_PRM_StressModel={sw['GSw_PRM_StressModel']} and " - f"stress_metric={stress_metric} is not supported." - ) - raise NotImplementedError(err) + ### Check that the stress metrics specified in GSw_PRM_StressThresholdMetrics + # have corresponding GSw_PRM_StressThreshold{metric} switches + stressThresholdMetricSwitches = [s.split('GSw_PRM_StressThreshold')[1] + for s in sw.keys() if s.startswith('GSw_PRM_StressThreshold') + and not s.endswith('Metrics') + ] + stressTresholdMetrics = sw['GSw_PRM_StressThresholdMetrics'].split('/') + + # Threshold metric added but not specified as a switch + for metric in stressTresholdMetrics: + if metric.upper() not in stressThresholdMetricSwitches: + raise NotImplementedError(f"GSw_PRM_StressThreshold{metric} is not specified as a switch.") + + # # Threshold metric is not added but specified as a switch + # for metric in stressThresholdMetricSwitches: + # if metric.upper() not in stressTresholdMetrics: + # raise NotImplementedError(f"GSw_PRM_StressThreshold{metric} is specified as a switch," + # f"but not added to GSw_PRM_StressThresholdMetrics list {stressTresholdMetrics}") + + for metric in stressThresholdMetricSwitches: + for threshold in sw[f'GSw_PRM_StressThreshold{metric}'].split('/'): + ## Example: threshold = 'transgrp_10_EUE_sum' + allowed_levels = ['country','interconnect','nercr','transreg','transgrp','st','r'] + (hierarchy_level, ppm, stress_metric, period_agg_method) = 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" + + '\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(ppm) >= 0): + raise ValueError( + f"ppm in GSw_PRM_StressThreshold{stress_metric} must be a positive number " + f"but '{ppm}' was provided" + ) + if stress_metric.upper() not in ['EUE','NEUE','LOLE']: + raise ValueError( + f"stress metric in GSw_PRM_StressThreshold{stress_metric} must be 'EUE', 'NEUE', or 'LOLE' " + f"but '{stress_metric}' was provided" + ) + if (sw['GSw_PRM_StressModel'].lower() != 'pras') and (stress_metric.upper() != 'EUE'): + err = ( + f"The combination of GSw_PRM_StressModel={sw['GSw_PRM_StressModel']} and " + f"stress_metric={stress_metric} is not supported." + ) + raise NotImplementedError(err) if sw['GSw_PRM_StressStorageCutoff'].lower() not in ['off','0','false']: metric, value = sw['GSw_PRM_StressStorageCutoff'].split('_')