From 45b2ffd70854cdcd2bff0214a02d8aba93beb9a2 Mon Sep 17 00:00:00 2001 From: Burcin Cakir Erdener Date: Fri, 24 Jul 2026 13:16:21 -0600 Subject: [PATCH] Add CVAR and NCVAR resource adequacy reporting --- cases.csv | 3 +- cases_test.csv | 4 +- reeds/resource_adequacy/ra_calcs.py | 15 +++ reeds/resource_adequacy/run_pras.jl | 91 ++++++++++++++---- reeds/resource_adequacy/stress_periods.py | 110 +++++++++++++++++++--- runreeds.py | 32 ++++--- 6 files changed, 209 insertions(+), 46 deletions(-) diff --git a/cases.csv b/cases.csv index 525a78bc7..3b302256d 100644 --- a/cases.csv +++ b/cases.csv @@ -273,13 +273,14 @@ GSw_PRM_StressOutages,Whether to apply the availability factor (forced + schedul GSw_PRM_StressSeedLoadLevel,Region hierarchy level at which to include peak coincident load days as seeded stress periods (or False to ignore peaks),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,transgrp, GSw_PRM_StressSeedMinRElevel,Region hierarchy level at which to include minimum wind and solar capacity factor days as seeded stress periods (or False to ignore min-wind/solar CF days),false; False; FALSE; r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region; ccreg,interconnect, GSw_PRM_StressStorageCutoff,How to select shoulder stress periods giving storage time to recharge before/after high-unserved-energy periods. Two-part switch separated by _. The first argument is 'EUE' or 'capacity' or 'absolute'. If 'EUE' the second argument specifies a PRAS storage headspace / EUE threshold [fraction]; if 'cap' it specifies a headspace / storage capacity threshold [fraction]; if 'abs' it specifies absolute number of periods before/after [integer]. Turned off if set to 'off'.,N/A,EUE_0.1, -GSw_PRM_StressThresholdMetrics,/-delimited list of metrics for identifying stress periods; supported metrics (case-insensitive) are: LOLD | LOLE | LOLH | NEUE | Depth | Duration,N/A,NEUE, +GSw_PRM_StressThresholdMetrics,/-delimited list of resource adequacy metrics. LOLD | LOLE | LOLH | NEUE | Depth | Duration identify stress periods; CVAR and NCVAR are report-only metrics and do not add stress periods or update PRM.,N/A,NEUE, GSw_PRM_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_CVARAlpha,Alpha value used for CVAR/NCVAR calculations,float,0.95, 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 f3af6b246..a3c8fa38a 100644 --- a/cases_test.csv +++ b/cases_test.csv @@ -3,7 +3,7 @@ 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.NJ,,,,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,z134,z54,z48,,,,,, +GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,z134,z54,z48,z132,,,,, GSw_GasCurve,2,,1,1,,,,,,,,,,,,,,,,,,,,1,1,,,,,, GSw_Geothermal,,,,2,,,,,,,,,,,,,,,,,,,0,,0,,,,,, GSw_GrowthPenalties,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,, @@ -51,7 +51,7 @@ 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_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth/CVAR/NCVAR,,,,, GSw_DRShed,,,,,,,,,,,,,,,,,,,,,,,,,,,1,,,, GSw_MGA_CostDelta,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.01,,, GSw_LoadSiteCF,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.95,, diff --git a/reeds/resource_adequacy/ra_calcs.py b/reeds/resource_adequacy/ra_calcs.py index cd325d75d..9bc2edff3 100644 --- a/reeds/resource_adequacy/ra_calcs.py +++ b/reeds/resource_adequacy/ra_calcs.py @@ -24,6 +24,7 @@ def run_pras( write_surplus=False, write_energy=False, write_shortfall_samples=False, + write_shortfall_samples_totals=False, write_availability_samples=False, **kwargs, ): @@ -68,7 +69,9 @@ def run_pras( f"--write_surplus={int(write_surplus)}", f"--write_energy={int(write_energy)}", f"--write_shortfall_samples={int(write_shortfall_samples)}", + f"--write_shortfall_samples_totals={int(write_shortfall_samples_totals)}", f"--write_availability_samples={int(write_availability_samples)}", + f"--cvar_alpha={float(sw['GSw_PRM_CVARAlpha'])}", f"--iteration={iteration}", f"--samples={sw['pras_samples']}", f"--overwrite={int(overwrite)}", @@ -153,13 +156,25 @@ def main(t, tnext, casedir, iteration=0): 2: True, }[int(sw['pras'])] if pras_this_solve_year or int(sw.GSw_PRM_StressIterateMax): + stress_metrics = [ + m.strip().upper() + for m in sw.GSw_PRM_StressThresholdMetrics.split('/') + if m.strip() + ] + write_shortfall_samples_totals = any( + metric in {'CVAR', 'NCVAR'} + for metric in stress_metrics + ) + result = run_pras( casedir, t, iteration=iteration, write_flow=(True if t == max(solveyears) else False), write_energy=True, write_shortfall_samples=(True if int(sw.GSw_PRM_UpdateMethod) > 1 else False), + write_shortfall_samples_totals=write_shortfall_samples_totals, ) + if result.returncode: raise Exception( f"run_pras.jl returned code {result.returncode}. Check gamslog.txt for error trace." diff --git a/reeds/resource_adequacy/run_pras.jl b/reeds/resource_adequacy/run_pras.jl index ec516737b..d39c7964d 100644 --- a/reeds/resource_adequacy/run_pras.jl +++ b/reeds/resource_adequacy/run_pras.jl @@ -70,10 +70,20 @@ function parse_commandline() default = 0 required = false "--write_shortfall_samples" - help = "Write the sample-level shortfall" + help = "Write per-sample hourly shortfall by region" arg_type = Int default = 0 required = false + "--write_shortfall_samples_totals" + help = "Write per-sample total shortfall by region over the full PRAS time period" + arg_type = Int + default = 0 + required = false + "--cvar_alpha" + help = "Alpha for CVaR (e.g., 0.95)" + arg_type = Float64 + default = 0.95 + required = false "--write_availability_samples" help = "Write the sample-level generator and storage availability" arg_type = Int @@ -182,7 +192,7 @@ function run_pras(pras_system_path::String, args::Dict) if args["write_energy"] == 1 resultspec["energy"] = PRAS.StorageEnergy() end - if args["write_shortfall_samples"] == 1 + if args["write_shortfall_samples"] == 1 || args["write_shortfall_samples_totals"] == 1 resultspec["short_samples"] = PRAS.ShortfallSamples() end if args["write_availability_samples"] == 1 @@ -208,6 +218,38 @@ function run_pras(pras_system_path::String, args::Dict) @info "$(PRAS.EUE(results["short"]))" @info "$(PRAS.NEUE(results["short"]))" + #%% Print CVAR and NCVAR for the entire modeled region + if args["write_shortfall_samples_totals"] == 1 && haskey(results, "short_samples") + _, _, _, _, energyunit = PRAS.get_params(sys) + alpha = Float64(args["cvar_alpha"]) + + cvar_result = PRAS.CVAR( + energyunit, + results["short_samples"], + alpha, + ) + + cvar_value = PRAS.val(cvar_result.cvar) + cvar_stderr = PRAS.stderror(cvar_result.cvar) + cvar_var = cvar_result.var + + ### Normalize CVaR by total load over the full PRAS time period. + total_load = sum(sys.regions.load) + + ncvar_value = cvar_value / total_load * 1e6 + ncvar_stderr = cvar_stderr / total_load * 1e6 + ncvar_var = cvar_var / total_load * 1e6 + + @info( + "CVAR = $(cvar_value)±$(cvar_stderr) MWh; " * + "VaR = $(cvar_var) MWh; alpha = $(alpha)" + ) + @info( + "NCVAR = $(ncvar_value)±$(ncvar_stderr) ppm; " * + "VaR = $(ncvar_var) ppm; alpha = $(alpha)" + ) + end + ## Filter out DC regions used for VSC HVDC transmission regions = [r for r in sys.regions.names if !(occursin("|", r))] @@ -285,6 +327,7 @@ function run_pras(pras_system_path::String, args::Dict) end @info("Wrote PRAS surplus to $(surplusfile)") end + ### Storage energy if args["write_energy"] == 1 dfenergy = DF.DataFrame() @@ -302,31 +345,41 @@ function run_pras(pras_system_path::String, args::Dict) @info("Wrote PRAS storage energy to $(energyfile)") end - ### Sample-level shortfall + ### Per-sample hourly shortfall by region if args["write_shortfall_samples"] == 1 - dictshort = Dict(s => DF.DataFrame() for s = 1:args["samples"]) - for s in range(1, args["samples"]) - dictshort[s] = DF.DataFrame( - transpose(getindex.(results["short_samples"][:, :], s)), - sys.regions.names - ) - # subset to regions (filter out DC regions) - dictshort[s] = dictshort[s][:,findall(regions .∈ Ref(sys.regions.names))] - end - ## Write it + sf = results["short_samples"] + ## Use filtered system regions to avoid DC converter pseudo-regions without load. + region_names = sf.regions.names + idx = [findfirst(==(r), region_names) for r in regions] + shortfile = replace(outfile, ".h5"=>"-shortfall_samples.h5") HDF5.h5open(shortfile, "w") do f - ## Create a group for each sample. Within each group, write an array for each region. - for s in range(1, args["samples"]) + for s in 1:args["samples"] HDF5.create_group(f, "$s") - for column in DF._names(dictshort[s]) - f["$s"]["$column", compress=4] = convert(Array, dictshort[s][!, column]) + for (r, i_r) in zip(regions, idx) + arr = Float64.(sf.shortfall[i_r, :, s]) + f["$s"]["$r", compress=4] = arr end end end @info("Wrote PRAS shortfall by sample to $(shortfile)") end + ### Per-sample total shortfall by region over the full PRAS time period, needed for CVaR + if args["write_shortfall_samples_totals"] == 1 + sf = results["short_samples"] + ## Use filtered system regions to avoid DC converter pseudo-regions without load. + totalsfile = replace(outfile, ".h5"=>"-shortfall_totals_by_sample.h5") + HDF5.h5open(totalsfile, "w") do f + f["sample", compress=4] = collect(1:args["samples"]) + f["USA", compress=4] = Float64.(sf[]) + for r in regions + f["$r", compress=4] = Float64.(sf[r]) + end + end + @info("Wrote PRAS shortfall totals by sample to $(totalsfile)") + end + ### Sample-level generator and storage availability if args["write_availability_samples"] == 1 dictavail = Dict(s => DF.DataFrame() for s = 1:args["samples"]) @@ -457,6 +510,7 @@ if abspath(PROGRAM_FILE) == @__FILE__ # "write_surplus" => 0, # "write_energy" => 0, # "write_shortfall_samples" => 1, + # "write_shortfall_samples_totals" => 1, # "write_availability_samples" => 0, # "overwrite" => 1, # "debug" => 0, @@ -466,6 +520,7 @@ if abspath(PROGRAM_FILE) == @__FILE__ # "pras_existing_unit_size" => 1, # "pras_max_unitsize_prm" => 1, # "pras_seed" => 1, + # "cvar_alpha" => 0.95, # ) # reedscase = args["reedscase"] # solve_year = args["solve_year"] @@ -485,4 +540,4 @@ if abspath(PROGRAM_FILE) == @__FILE__ main(args) #%% -end +end \ No newline at end of file diff --git a/reeds/resource_adequacy/stress_periods.py b/reeds/resource_adequacy/stress_periods.py index 7f41e0fd0..20150ae7e 100644 --- a/reeds/resource_adequacy/stress_periods.py +++ b/reeds/resource_adequacy/stress_periods.py @@ -14,6 +14,8 @@ # import importlib # importlib.reload(functions) +CVAR_METRICS = {'CVAR', 'NCVAR'} + #%%### Constants RA_SWITCHES = { @@ -142,6 +144,51 @@ def calc_neue(dfeue_agg, dfload_agg): neue = dfeue_agg.sum() / dfload_agg.sum() * 1e6 return neue +def get_cvar_alpha(sw): + alpha = float(sw.GSw_PRM_CVARAlpha) + if not (0 <= alpha < 1): + raise ValueError(f"GSw_PRM_CVARAlpha must be in [0, 1). Got {alpha}") + return alpha + + +def get_shortfall_totals_by_sample(case, t, iteration=0): + filepath = os.path.join(case, 'handoff', 'PRAS', f'PRAS_{t}i{iteration}-shortfall_totals_by_sample.h5') + if not os.path.isfile(filepath): + raise FileNotFoundError(f"{filepath} not found. Re-run PRAS with --write_shortfall_samples_totals 1.") + df = reeds.io.read_pras_results(filepath) + df.columns = df.columns.astype(str) + if 'sample' in df.columns: + df = df.set_index('sample') + elif df.index.name != 'sample': + df.index = pd.RangeIndex(1, len(df) + 1, name='sample') + return df.apply(pd.to_numeric, errors='coerce').clip(lower=0) + + +def _sample_cvar(samples, alpha=0.95): + x = pd.Series(samples).dropna().astype(float) + if x.empty: + return np.nan + n_tail = max(1, int(np.ceil(round((1 - alpha) * len(x), 12))),) + return x.sort_values(ascending=False).iloc[:n_tail].mean() + +def calc_cvar(shortfall_samples_agg, alpha=0.95): + """ + CVAR from total shortfall by PRAS sample. + """ + cvar = shortfall_samples_agg.apply( + lambda s: _sample_cvar(s, alpha=alpha), + axis=0, + ) + return cvar + + +def calc_ncvar(cvar, dfload_agg): + """ + NCVAR = CVAR / total load, in ppm. + """ + ncvar = cvar / dfload_agg.sum().reindex(cvar.index) * 1e6 + return ncvar + def calc_peak_eue(dfeue_agg, dfload_agg, norm:Literal['peak','hourly','absolute']='peak'): """ @@ -192,6 +239,19 @@ def calc_ra_metrics( sw = reeds.io.get_switches(case) numyears = len(sw.resource_adequacy_years_list) + cvar_metrics = { + i.strip().upper() + for i in sw.GSw_PRM_StressThresholdMetrics.split('/') + if i.strip().upper() in CVAR_METRICS + } + + if len(cvar_metrics): + shortfall_samples = ( + get_shortfall_totals_by_sample(case=case, t=t, iteration=iteration) + .drop(columns=['USA'], errors='ignore') + ) + cvar_alpha = get_cvar_alpha(sw) + ### Loop over aggregation levels and calculate all metrics ra_metrics = {} for level in levels: @@ -200,9 +260,9 @@ def calc_ra_metrics( 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).T.groupby(level=0).max().T - dfeue_agg = dfeue.rename(columns=rmap).T.groupby(level=0).sum().T - dfload_agg = dfload.rename(columns=rmap).T.groupby(level=0).sum().T + 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 @@ -212,6 +272,22 @@ def calc_ra_metrics( 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') + if len(cvar_metrics): + regions = [c for c in shortfall_samples.columns if c in rmap.index] + if len(regions): + shortfall_samples_agg = ( + shortfall_samples[regions] + .rename(columns=rmap) + .groupby(axis=1, level=0) + .sum() + ) + cvar = calc_cvar(shortfall_samples_agg, alpha=cvar_alpha) + + if 'CVAR' in cvar_metrics: + ra_metrics[level, 'cvar_mwh_peryear'] = cvar / numyears + + if 'NCVAR' in cvar_metrics: + ra_metrics[level, 'ncvar_ppm'] = calc_ncvar(cvar, dfload_agg) ### Combine it dfout = pd.concat(ra_metrics, names=['level','metric','region']).rename('value') @@ -232,7 +308,7 @@ def get_eue_events( events = {} for level in levels: rmap = reeds.io.get_rmap(case=case, hierarchy_level=level) - dfeue_agg = dfeue.rename(columns=rmap).T.groupby(level=0).sum().T + 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 @@ -268,7 +344,7 @@ def get_longest_events( dates = [] for i, row in eue_events.iterrows(): dates.append( - pd.Series(index=pd.date_range(row.start, row.end, freq='h'), data=1) + pd.Series(index=pd.date_range(row.start, row.end, freq='H'), data=1) .resample('D').count() ) if len(dates): @@ -486,7 +562,7 @@ def get_stress_periods(case, sw, t, iteration): dfenergy = ( dfenergy_unit .rename(columns={c: c.split('|')[1] for c in dfenergy_unit.columns}) - .T.groupby(level=0).sum().T + .groupby(axis=1, level=0).sum() ) ### Load this year's stress periods so we don't duplicate @@ -501,23 +577,26 @@ def get_stress_periods(case, sw, t, iteration): stressperiods_this_iteration['start'] + ( (pd.Timedelta('5D') if sw.GSw_HourlyType == 'wek' else pd.Timedelta('1D')) - - pd.Timedelta('1h') + - 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') + pd.date_range(row.start, row.end, freq='1H') for i,row in stressperiods_this_iteration.iterrows() ] covered_hours = [i for sublist in covered_hours for i in sublist] - - ### Check all stress criteria; for regions that fail, add new stress periods _failed = {} _high_stress_periods = {} _shoulder_periods = {} - stress_metrics = [i.lower() for i in sw.GSw_PRM_StressThresholdMetrics.split('/')] + stress_metrics = [] + for metric in sw.GSw_PRM_StressThresholdMetrics.split('/'): + metric = str(metric).strip() + if metric and metric.upper() not in CVAR_METRICS: + stress_metrics.append(metric.lower()) + for stress_metric in stress_metrics: switch = RA_SWITCHES[stress_metric] for criterion in sw[switch].split('/'): @@ -525,9 +604,9 @@ def get_stress_periods(case, sw, t, iteration): ## Example: criterion = 'transgrp_1' hierarchy_level, metric_threshold = criterion.split('_') rmap = reeds.io.get_rmap(case=case, hierarchy_level=hierarchy_level) - dfeue_agg = dfeue.rename(columns=rmap).T.groupby(level=0).sum().T.drop(covered_hours) - dflole_agg = dflole.rename(columns=rmap).T.groupby(level=0).max().T.drop(covered_hours) - dfenergy_agg = dfenergy.rename(columns=rmap).T.groupby(level=0).sum().T.drop(covered_hours) + 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, @@ -865,6 +944,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__': diff --git a/runreeds.py b/runreeds.py index 2688a713c..2c5b175bd 100644 --- a/runreeds.py +++ b/runreeds.py @@ -300,21 +300,28 @@ def check_compatibility(sw): 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('/')] + report_only_metrics = {'cvar', 'ncvar'} + + used_metrics = [i.strip().lower() for i in sw['GSw_PRM_StressThresholdMetrics'].split('/') if i.strip()] allowed_levels = ['country','interconnect','nercr','transreg','transgrp','st','r'] for metric in used_metrics: + if metric in report_only_metrics: + continue + 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('_') + 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 " @@ -835,19 +842,22 @@ def setupEnvironment( #%% Check whether the ReEDS conda environment is activated if (not skip_checks) and ( - ('reeds' not in os.environ['CONDA_DEFAULT_ENV'].lower()) - or (not pd.__version__.startswith('3')) + ('reeds2' not in os.environ['CONDA_DEFAULT_ENV'].lower()) + or (not pd.__version__.startswith('2')) ): - err = ( + print( f"Your environment is {os.environ['CONDA_DEFAULT_ENV']} and your pandas " - f"version is {pd.__version__}.\nThe supported environment is 'reeds', with\n" - "pandas version 3.x.\n" + f"version is {pd.__version__}.\nThe default environment is 'reeds2', with\n" + "pandas version 2.x, so the python parts of ReEDS are unlikely to work.\n" "To build the environment for the first time, run:\n" " `conda env create -f environment.yml`\n" "To activate the created environment, run:\n" - " `conda activate reeds` (or `activate reeds` on Windows)" + " `conda activate reeds2` (or `activate reeds2` on Windows)\n" + "Do you want to continue without activating the environment?" ) - raise ValueError(err) + confirm_env = str(input("Continue? y/[n]: ") or 'n') + if confirm_env not in ['y','Y','yes','Yes','YES']: + quit() #%% Load specified case file, infer other settings from cases.csv if cases_suffix in ['', 'default']: @@ -1275,7 +1285,7 @@ def write_batch_script( OPATH.writelines("module load conda \n") OPATH.writelines("module load gams \n") - OPATH.writelines("conda activate reeds \n") + OPATH.writelines("conda activate reeds2 \n") OPATH.writelines('export R_LIBS_USER="$HOME/rlib" \n\n\n') #%% Write the input_processing script calls @@ -1283,7 +1293,6 @@ def write_batch_script( for s in [ 'copy_files', 'mcs_sampler', - 'climateprep', 'hydcf', 'h2_storage', 'calc_financial_inputs', @@ -1292,6 +1301,7 @@ def write_batch_script( 'writesupplycurves', 'writedrshift', 'plantcostprep', + 'climateprep', 'hourly_load', 'recf', 'forecast',