diff --git a/cases.csv b/cases.csv index c856b1a3..1404fc43 100644 --- a/cases.csv +++ b/cases.csv @@ -263,14 +263,20 @@ 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_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, +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/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,,,,,, 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 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 cda0f96e..df648cce 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)', ] @@ -473,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'] @@ -1436,7 +1435,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 +1548,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 +2499,6 @@ 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] wide = 1 if len(hierarchy[basecase]['transreg'].unique()) > 6 else 0 weatheryear = sw.GSw_HourlyWeatherYears.split('_')[0] metrics = [ @@ -2515,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-sum-{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), ] + [ ( diff --git a/postprocessing/single_case_plots.py b/postprocessing/single_case_plots.py index 5acf8dfa..b64a5b59 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 @@ -71,7 +72,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,8 +581,8 @@ try: plt.close() levels = ['country', 'interconnect', 'transreg', 'transgrp'] - f, ax, _ = reedsplots.plot_neue_bylevel(case=case, levels=levels) - savename = f"plot_stressperiod_neue-{','.join(levels)}.png" + f, ax, _ = reedsplots.plot_ra_metrics_bylevel(case=case, levels=levels) + savename = f"plot_ra_metrics-{','.join(levels)}.png" if write: plt.savefig(os.path.join(savepath, savename)) if interactive: @@ -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: @@ -717,17 +739,16 @@ print(traceback.format_exc()) try: - level, threshold, _, metric = sw['GSw_PRM_StressThreshold'].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' - 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 99764d4e..7840111f 100644 --- a/reeds/reedsplots.py +++ b/reeds/reedsplots.py @@ -465,7 +465,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) @@ -4120,30 +4120,46 @@ def plot_stressperiod_days(case, repcolor='k', sharey=False, figsize=(10,5)): def plot_stressperiod_evolution( - case, level=None, metric=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 + ### 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) - _level, _threshold, _, _metric = sw['GSw_PRM_StressThreshold'].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'))) - dictin_neue = { - tuple([int(x) for x in os.path.basename(f)[len('neue_'):-len('.csv')].split('i')]): + 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) * scale.get(metric, 1) + ### Load RA results + infiles = sorted(glob(os.path.join(case,'outputs','ra_metrics_*.csv'))) + 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(metric,0,'metric') - .NEUE_ppm.unstack('region') - ) + .xs(switch_metric,0,'metric') + .squeeze(1).unstack('region') + ) * scale.get(metric, 1) ### Load stress periods for labels dfstress = get_stressperiods(case) ### Plot setup @@ -4196,25 +4212,33 @@ 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_ylim(0) + ax[0].set_ylabel(ylabel[metric]) + ax[0].set_ylim(0, min(ax[0].get_ylim()[1], threshold*50)) plots.despine(ax) 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=2026, + 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 @@ -4223,57 +4247,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) >= tmin]].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 = { - i.split('_')[0]: float(i.split('_')[1]) - for i in sw.GSw_PRM_StressThreshold.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, @@ -4288,7 +4325,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( @@ -4312,11 +4349,11 @@ 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 = neue.loc[neue.metric==metric].set_index(['level','region']).NEUE_ppm + 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_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'] @@ -6179,7 +6216,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 = { @@ -6510,6 +6548,102 @@ 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.9, + showhull=True, +): + """ + 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':'%', '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=True, + ) + 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 + 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) + 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') + 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 + + 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, @@ -6621,8 +6755,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 @@ -6640,7 +6774,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/diagnostic_plots.py b/reeds/resource_adequacy/diagnostic_plots.py index 36145664..a4487d39 100644 --- a/reeds/resource_adequacy/diagnostic_plots.py +++ b/reeds/resource_adequacy/diagnostic_plots.py @@ -85,6 +85,12 @@ def get_inputs(sw): resources['tech'] = reeds.reedsplots.simplify_techs(resources.i, display_level = 'diagnostics') resources['rb'] = resources.r + 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 try: @@ -242,6 +248,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 @@ -544,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() @@ -1023,6 +1028,37 @@ 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'] + if not len(new_stress_periods): + print('No new stress periods to plot') + return + 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 +1296,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/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/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..ec90b33a 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -1,11 +1,9 @@ #%%### 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 +import pandas as pd +from pathlib import Path +from typing import Literal import reeds from reeds.input_processing import hourly_writetimeseries @@ -17,9 +15,25 @@ # importlib.reload(functions) +#%%### Constants +RA_SWITCHES = { + 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 -def get_pras_eue(case, t, iteration=0): +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( @@ -29,214 +43,248 @@ 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[[ - c for c in dfpras - if (c.endswith(eue_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 - - -def get_eue_periods( - case, t, iteration=0, - hierarchy_level='transgrp', - stress_metric='EUE', - period_agg_method='sum', - ): - """_summary_ + ### Keep the metric columns by zone + 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] + )) + ) - 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'. + return dictout - 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 EUE from PRAS - dfeue = get_pras_eue(case=case, t=t, iteration=iteration) - ## Aggregate to hierarchy_level - dfeue = ( - dfeue - .rename_axis('r', axis=1).rename_axis('h', axis=0) - .rename(columns=rmap).groupby(axis=1, level=0).sum() +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 ) - - ###### Calculate the stress metric by period - if stress_metric.upper() == 'EUE': - ### Aggregate according to period_agg_method - dfmetric_period = ( - dfeue - .groupby([dfeue.index.year, dfeue.index.month, dfeue.index.day]) - .agg(period_agg_method) - .rename_axis(['y','m','d']) - ) - elif stress_metric.upper() == 'NEUE': - ### 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 - - ### Recalculate NEUE [ppm] and aggregate appropriately - if period_agg_method == 'sum': - dfmetric_period = ( - dfeue - .groupby([dfeue.index.year, dfeue.index.month, dfeue.index.day]) - .agg(period_agg_method) - .rename_axis(['y','m','d']) - ) / ( - dfload - .groupby([dfload.index.year, dfload.index.month, dfload.index.day]) - .agg(period_agg_method) - .rename_axis(['y','m','d']) - ) * 1e6 - elif period_agg_method == 'max': - dfmetric_period = ( - (dfeue / dfload) - .groupby([dfeue.index.year, dfeue.index.month, dfeue.index.day]) - .agg(period_agg_method) - .rename_axis(['y','m','d']) - ) * 1e6 - - ### 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') + 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 ) - ## 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 + ends = ends.loc[ends > 0].index + assert len(starts) == len(ends), "Error in event start/end calculation" + ## Get some metrics for each event + events = [] + for start, end in zip(starts, ends): + event = ds.loc[start:end] + events.append({ + 'start': start, + 'end': end, + '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','profile']) + return dfout + + +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.resample('D').max() + ## Sum the probability across days + lold = daily_max.sum() + return lold + + +def calc_lole(dflole_agg, threshold=0): + """Number of events, where an event is >threshold LOLE in contiguous hours""" + ## 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 -def plot_eue_diagnostics(sw, t, iteration, high_eue_periods): - try: - dates = ( - pd.concat(high_eue_periods) - .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 as err: - print(err) +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 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. +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 - 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')) +def calc_neue(dfeue_agg, dfload_agg): + """NEUE (sum of EUE / sum of load) in units of ppm""" + neue = dfeue_agg.sum() / dfload_agg.sum() * 1e6 return neue -def get_annual_neue(case, t, iteration=0): - """ +def calc_peak_eue(dfeue_agg, dfload_agg, norm:Literal['peak','hourly','absolute']='peak'): """ - ### Get EUE from PRAS - dfeue = get_pras_eue(case=case, t=t, iteration=iteration) + Get the peak hourly outage magnitude - ### Get load (for calculating NEUE) + 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 + """ + 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_agg.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 values from PRAS + pras_shortfall = get_pras_shortfall(case, t, iteration) + dfeue = pras_shortfall['EUE'] + dflole = pras_shortfall['LOLE'] + + ### Get load and number of years (for normalization) dfload = reeds.io.read_h5py_file( - os.path.join( - case,'handoff','reeds_data',f'pras_load_{t}.h5') + Path(case,'handoff','reeds_data',f'pras_load_{t}.h5') ) dfload.index = dfeue.index - - levels = ['country','interconnect','nercr','transreg','transgrp','st','r'] - _neue = {} - 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() - / 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() - / dfload.rename(columns=rmap).groupby(axis=1, level=0).sum() - ).max() * 1e6 + sw = reeds.io.get_switches(case) + numyears = len(sw.resource_adequacy_years_list) + + ### Loop over aggregation levels and calculate all metrics + ra_metrics = {} + 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=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[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') + 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 - neue = pd.concat(_neue, names=['level','metric','region']).rename('NEUE_ppm') - - return neue + dfout = pd.concat(ra_metrics, names=['level','metric','region']).rename('value') + + 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_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( + dsmetric, + sw, + num_events:int=1, +): + eue_events = ( + get_events(dsmetric) + .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() + ) + 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_r, high_eue_periods): +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( @@ -244,7 +292,7 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods): "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" @@ -252,39 +300,34 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods): 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, ppm, stress_metric, period_agg_method) = criterion.split('_') + fmt = '%Y-%m-%d' - ## 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: + _shoulder_periods = {} + 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(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.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(fmt), row.region].iloc[0] + end_headspace_frac = dfheadspace_frac.loc[day.strftime(fmt), row.region].iloc[-1] - day_eue = high_eue_periods[criterion, f'high_{stress_metric}'].loc[i,'EUE'] + 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] + ## 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 ( @@ -292,141 +335,251 @@ def get_shoulder_periods(sw, criterion, dfenergy_r, high_eue_periods): 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'), - 'y':day_after.year, 'm':day_after.month, 'd':day_after.day, 'r':row.r, - }).to_frame().T.set_index('actual_period') - print(f"Added {day_after} as shoulder stress period after {day}") + _shoulder_periods[f'after_{row.period}'] = pd.Series({ + 'region':row.region, + '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'), - 'y':day_before.year, 'm':day_before.month, 'd':day_before.day, 'r':row.r, - }).to_frame().T.set_index('actual_period') - print(f"Added {day_before} as shoulder stress period before {day}") - + _shoulder_periods[f'before_{row.period}'] = pd.Series({ + 'region':row.region, + 'period':day_before.strftime('y%Yd%j'), + }) + print( + f"Added {day_before.strftime(fmt)} as shoulder stress period " + f"before {day.strftime(fmt)}" + ) + 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 -def get_eue_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( +def check_threshold_and_choose_periods( + stress_metric:str, + criterion:str, + sw, + t:int, + iteration:int, + dfeue_agg, + dflole_agg, + dfenergy_agg, +): + ## 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 + 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) + + ### Get the threshold(s) and see if any of them failed + 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:") + for i, val in this_test.items(): + print(f'{i}: {val} {stress_metric}') + else: + print(f"{RA_SWITCHES[stress_metric]} = {criterion} failed for:") + for i, val in failed.items(): + print(f'{i}: {val} {stress_metric}') + ## Get new stress periods since the metric 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': + ## 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( + dfeue_agg[region], sw=sw, + num_events=int(sw.GSw_PRM_StressIncrement), + ) + 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('value') + .reset_index() + ) + ### Include "shoulder periods" before or after each period + ### if the storage state of charge is low + if stress_metric.lower() == 'neue': + shoulder_periods = get_shoulder_periods( + sw, + criterion, + dfenergy_agg, + high_stress_periods, + ) + else: + shoulder_periods = pd.DataFrame() + + return { + 'failed': failed, + 'high_stress_periods': high_stress_periods, + 'shoulder_periods': shoulder_periods, + } + + +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() ) - ### 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( sw['casedir'], 'inputs_case', f'stress{t}i{iteration}', 'period_szn.csv') ) - - ### 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') + 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] - ### 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)) + ### Check all stress criteria; for regions that fail, add new stress 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: + 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().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, + criterion, + sw, + t, + iteration, + dfeue_agg, + dflole_agg, + dfenergy_agg, ) - 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']) + if dictout is not None: + _failed[stress_metric, criterion] = dictout['failed'] + _high_stress_periods[stress_metric, criterion] = dictout['high_stress_periods'] + _shoulder_periods[stress_metric, criterion] = dictout['shoulder_periods'] ### Get lists of stress periods: new (added this iteration) and all - if len(failed): + if len(_failed): + failed = pd.concat(_failed, names=['stress_metric', 'criterion']) + high_stress_periods = pd.concat(_high_stress_periods) + shoulder_periods = pd.concat(_shoulder_periods) new_stress_periods = pd.concat( - {**high_eue_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) + ## 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 failed, None, None + return pd.DataFrame(), pd.DataFrame(), pd.DataFrame() ## 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') @@ -444,15 +597,11 @@ def get_eue_sorted_periods(sw, t, iteration): else: 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') - ) - new_stress_periods.round(2).rename(columns={'EUE':'EUE_MWh','NEUE':'NEUE_ppm'}).to_csv( + ### Tables for debugging + new_stress_periods.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) return failed, new_stressperiods_write, combined_periods_write @@ -588,26 +737,44 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): sw (pd.series): ReEDS switches for this run. t (int): Model solve year. 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 + failed (pd.Series): Series of failed regions with unserved energy at the + hierarchy_level, indexed by (stress_metric, criterion, region) + combined_periods_write (pd.DataFrame): Data frame of combined stress periods Returns: pd.DataFrame: Table of prm levels for the next PRAS iteration """ + ## Existing PRM, in case there's nothing to update it with below + prm = pd.read_csv( + os.path.join(sw['casedir'], 'inputs_case', f'stress{t}i{iteration}', 'prm.csv'), + index_col='*r', + ).fraction + # Get regions that failed criteria + # Use NEUE-based failed regions only, since the PRM update is computed from the + # NEUE-based shortfall (in ppm) and doesn't generalize to the other stress metrics. + # NEUE may not appear in `failed` at all this iteration (either it's not a configured + # stress metric, or it simply didn't fail while another metric did) -- in that case + # there's nothing to update the PRM with, so leave it unchanged. _failed_regions = [] - for criterion in failed: - # Example: criterion = 'transgrp_10_EUE_sum' - (hierarchy_level, ppm, __, __) = 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['ppm'] = float(ppm) - _failed_regions.append(df) + if 'neue' in failed.index.get_level_values('stress_metric'): + neue_failed = failed.xs('neue', level='stress_metric') + for criterion in neue_failed.index.get_level_values('criterion').unique(): + sub = neue_failed.xs(criterion, level='criterion') + # 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(sub.index) + ].rename(columns={hierarchy_level:'region'}) + df['hierarchy_level'] = hierarchy_level + df['ppm'] = float(metric_threshold) + _failed_regions.append(df) + + if not len(_failed_regions): + return prm.round(3) + # For zones that failed multiple criteria, use the most stringent (lowest EUE target) failed_regions = ( pd.concat(_failed_regions) @@ -634,10 +801,6 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): ) ## Add the PRM increment to last iteration's PRM - prm = pd.read_csv( - os.path.join(sw['casedir'], 'inputs_case', f'stress{t}i{iteration}', 'prm.csv'), - index_col='*r', - ).fraction prm_next_iteration = prm.add(prm_increment, fill_value=0).round(3) return prm_next_iteration @@ -647,27 +810,27 @@ def update_prm(sw, t, iteration, failed, combined_periods_write): def main(sw, t, iteration=0, logging=True): """ """ - #%% Write consolidated NEUE 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") - ) + #%% 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') + ) - except Exception as err: - if int(sw['pras']) == 2: - print(traceback.format_exc()) - if int(sw.GSw_PRM_StressIterateMax): - raise Exception(err) + #%% 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') + ) #%% 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_eue_sorted_periods( - sw=sw, t=t, iteration=iteration, + failed, new_stressperiods_write, combined_periods_write = get_stress_periods( + case=sw.casedir, sw=sw, t=t, iteration=iteration, ) #%% Stop here if all thresholds pass or if there are no new stress periods @@ -705,12 +868,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 diff --git a/runreeds.py b/runreeds.py index d76a6149..8ab81b84 100644 --- a/runreeds.py +++ b/runreeds.py @@ -294,34 +294,42 @@ 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 + ### 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"{ra_switches[metric]}: level={hierarchy_level} but must be in:\n" + + '\n'.join(allowed_levels) + ) + if not (float(stress_value) >= 0): + raise ValueError( + f"stress value in {ra_switches[metric]} must be a positive number " + f"but '{stress_value}' was provided" + ) + + ### GSw_PRM_UpdateMethod 1-3 (static or PRAS-informed PRM update) is computed from the + ### NEUE-based shortfall, so it requires NEUE to be an active stress metric + if int(sw['GSw_PRM_UpdateMethod']) in [1, 2, 3] and 'neue' not in used_metrics: + raise ValueError( + f"GSw_PRM_UpdateMethod={sw['GSw_PRM_UpdateMethod']} requires 'NEUE' to be included " + f"in GSw_PRM_StressThresholdMetrics (={sw['GSw_PRM_StressThresholdMetrics']}), " + "since PRM updates are computed from the NEUE-based shortfall." + ) + 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']: