diff --git a/reeds/core/terminus/report.gms b/reeds/core/terminus/report.gms index 23d6e7fd..40bc37cc 100644 --- a/reeds/core/terminus/report.gms +++ b/reeds/core/terminus/report.gms @@ -122,34 +122,53 @@ $include autocode%ds%report_params.gms h(h)$[not h_rep(h)] = no ; szn(szn)$[not szn_rep(szn)] = no ; -*================================================= -* -- CAPACITY ABOVE INTERCONNECTION QUEUE LIMIT -- -*================================================= - -cap_above_limit(tg,r,t)$tmodel_new(t) = CAP_ABOVE_LIM.l(tg,r,t) ; - -*===================== -* -- CO2 Reporting -- -*===================== - -CO2_CAPTURED_out(r,h,t)$tmodel_new(t) = CO2_CAPTURED.l(r,h,t) ; -CO2_CAPTURED_out_ann(r,t)$tmodel_new(t) = sum(h,hours(h) * CO2_CAPTURED.l(r,h,t) ); -CO2_STORED_out(r,cs,h,t)$[tmodel_new(t)$csfeas(cs)] = CO2_STORED.l(r,cs,h,t) ; -CO2_STORED_out_ann(r,cs,t)$[tmodel_new(t)$csfeas(cs)] = sum(h,hours(h) * CO2_STORED.l(r,cs,h,t) ); -CO2_TRANSPORT_INV_out(r,rr,t)$tmodel_new(t) = CO2_TRANSPORT_INV.l(r,rr,t) ; -CO2_SPURLINE_INV_out(r,cs,t)$[tmodel_new(t)$csfeas(cs)] = CO2_SPURLINE_INV.l(r,cs,t) ; - -CO2_FLOW_out(r,rr,h,t)$[(ord(r) < ord(rr))$tmodel_new(t)] = CO2_FLOW.l(r,rr,h,t) + CO2_FLOW.l(rr,r,h,t) ; -CO2_FLOW_out_ann(r,rr,t)$[(ord(r) < ord(rr))$tmodel_new(t)] = sum{h, hours(h) * (CO2_FLOW.l(r,rr,h,t) + CO2_FLOW.l(rr,r,h,t)) } ; - -CO2_FLOW_pos_out(r,rr,h,t)$[(ord(r) < ord(rr))$tmodel_new(t)] = CO2_FLOW.l(r,rr,h,t) ; -CO2_FLOW_pos_out_ann(r,rr,t)$[(ord(r) < ord(rr))$tmodel_new(t)] = sum{h, hours(h) * CO2_FLOW.l(r,rr,h,t) } ; - -CO2_FLOW_neg_out(r,rr,h,t)$[(ord(r) < ord(rr))$tmodel_new(t)] = -1 * CO2_FLOW.l(rr,r,h,t) ; -CO2_FLOW_neg_out_ann(r,rr,t)$[(ord(r) < ord(rr))$tmodel_new(t)] = -1 * sum{h, hours(h) * CO2_FLOW.l(rr,r,h,t) } ; - -CO2_FLOW_net_out(r,rr,h,t)$[(ord(r) < ord(rr))$tmodel_new(t)] = CO2_FLOW.l(r,rr,h,t) - CO2_FLOW.l(rr,r,h,t) ; -CO2_FLOW_net_out_ann(r,rr,t)$[(ord(r) < ord(rr))$tmodel_new(t)] = sum{h, hours(h) * (CO2_FLOW.l(r,rr,h,t) - CO2_FLOW.l(rr,r,h,t)) } ; +*===================================== +* -- Parameters moved to Python -- +*===================================== + +$ontext +The calculation of the following output parameters has been moved to report_calcs.py: + +Interconnection queue + cap_above_limit + +CO2 storage and flows + CO2_CAPTURED_out + CO2_CAPTURED_out_ann + CO2_FLOW_neg_out + CO2_FLOW_neg_out_ann + CO2_FLOW_net_out + CO2_FLOW_net_out_ann + CO2_FLOW_out + CO2_FLOW_out_ann + CO2_FLOW_pos_out + CO2_FLOW_pos_out_ann + CO2_SPURLINE_INV_out + CO2_STORED_out + CO2_STORED_out_ann + CO2_TRANSPORT_INV_out + +Transmission + cap_converter_out + invtran_out + tran_cap_energy + tran_cap_grp + tran_cap_prm + tran_flow_all_rep + tran_flow_all_stress + tran_flow_rep + tran_flow_rep_ann + tran_flow_stress + tran_mi_out + tran_mi_out_detail + tran_out + tran_prm_mi_out + tran_prm_out + tran_util_ann_rep + tran_util_ann_stress + tran_util_h_rep + tran_util_h_stress +$offtext *========================= * LCOE @@ -1774,75 +1793,6 @@ excess_load(r,h,t) = EXCESS.l(r,h,t) ; *====================== * Transmission *====================== - -invtran_out(r,rr,trtype,t)$routes_inv(r,rr,trtype,t) = INVTRAN.l(r,rr,trtype,t) ; - -tran_cap_energy(r,rr,trtype,t)$routes(r,rr,trtype,t) = CAPTRAN_ENERGY.l(r,rr,trtype,t) ; -tran_cap_prm(r,rr,trtype,t)$routes(r,rr,trtype,t) = CAPTRAN_PRM.l(r,rr,trtype,t) ; -tran_cap_grp(transgrp,transgrpp,t)$trancap_init_transgroup(transgrp,transgrpp,"AC") - = CAPTRAN_GRP.l(transgrp,transgrpp,t) ; - -tran_out(r,rr,trtype,t)$[(ord(r) df.index.get_level_values('rr') + if isinstance(df, pd.Series): + _out = df.loc[mask] + elif isinstance(df, pd.DataFrame): + _out = df.loc[mask, value] + out = _out.rename_axis(['rr', 'r'] + other_indices).reorder_levels(r_indices + other_indices) + return out + + +def combine_forward_reverse(df, agg:Literal['net','simult']='net', value='Level'): + """Combine forward (r < rr) and reverse (r > rr) into one +/- series""" + r_indices = ['r', 'rr'] + other_indices = [i for i in df.index.names if i not in r_indices] + forward = get_flow(df, 'forward') + reverse = (-1 if agg == 'net' else 1) * get_flow(df, 'reverse') + return pd.concat([forward, reverse]).groupby(r_indices + other_indices).sum() + + +#%% Results calculations +def calc_iq(g): + """Capacity above interconnection queue limit""" + dfs = {} + ## (tg,r,t) + dfs['cap_above_limit'] = g['CAP_ABOVE_LIM'].Level + return dfs + + +def calc_co2_stor(g): + """CO2 capture, transport, and storage""" + dfs = {} + ## (r,h,t) + dfs['CO2_CAPTURED_out'] = g['CO2_CAPTURED'].Level + ## (r,t) + dfs['CO2_CAPTURED_out_ann'] = (g['CO2_CAPTURED'].Level * g['hours']).groupby(['r','t']).sum() + ## (r,cs,h,t) + dfs['CO2_STORED_out'] = g['CO2_STORED'].Level + ## (r,cs,t) + dfs['CO2_STORED_out_ann'] = (g['CO2_STORED'].Level * g['hours']).groupby(['r','cs','t']).sum() + ## (r,rr,t) + dfs['CO2_TRANSPORT_INV_out'] = g['CO2_TRANSPORT_INV'].Level + ## (r,cs,t) + dfs['CO2_SPURLINE_INV_out'] = g['CO2_SPURLINE_INV'].Level + ## (r,rr,h,t) + dfs['CO2_FLOW_out'] = combine_forward_reverse(g['CO2_FLOW'], agg='simult') + ## (r,rr,t) + dfs['CO2_FLOW_out_ann'] = (dfs['CO2_FLOW_out'] * g['hours']).groupby(['r','rr','t']).sum() + ## (r,rr,h,t) + dfs['CO2_FLOW_pos_out'] = get_flow(g['CO2_FLOW'], 'forward') + ## (r,rr,t) + dfs['CO2_FLOW_pos_out_ann'] = (dfs['CO2_FLOW_pos_out'] * g['hours']).groupby(['r','rr','t']).sum() + ## (r,rr,h,t) + dfs['CO2_FLOW_neg_out'] = -get_flow(g['CO2_FLOW'], 'reverse') + ## (r,rr,t) + dfs['CO2_FLOW_neg_out_ann'] = (dfs['CO2_FLOW_neg_out'] * g['hours']).groupby(['r','rr','t']).sum() + ## (r,rr,h,t) + dfs['CO2_FLOW_net_out'] = combine_forward_reverse(g['CO2_FLOW'], agg='net') + ## (r,rr,t) + dfs['CO2_FLOW_net_out_ann'] = (dfs['CO2_FLOW_net_out'] * g['hours']).groupby(['r','rr','t']).sum() + return dfs + + +def calc_transmission(g): + """Transmission capacity and flow""" + dfs = {} + ## (r,rr,trtype,t) + dfs['invtran_out'] = g['INVTRAN'].Level + ## (r,rr,trtype,t) + dfs['tran_cap_energy'] = g['CAPTRAN_ENERGY'].Level + ## (r,rr,trtype,t) + dfs['tran_cap_prm'] = g['CAPTRAN_PRM'].Level + ## (transgrp,transgrpp,t) + dfs['tran_cap_grp'] = g['CAPTRAN_GRP'].Level + ## (r,rr,trtype,t) + dfs['tran_out'] = combine_forward_reverse(dfs['tran_cap_energy'], agg='simult') / 2 + ## (r,rr,trtype,t) + dfs['tran_prm_out'] = combine_forward_reverse(dfs['tran_cap_prm'], agg='simult') / 2 + ## (r,rr,trtype,t) + dfs['tran_mi_out_detail'] = dfs['tran_out'] * g['distance'] + ## (trtype,t) + dfs['tran_mi_out'] = dfs['tran_mi_out_detail'].groupby(['trtype','t']).sum() + ## (trtype,t) + dfs['tran_prm_mi_out'] = (dfs['tran_prm_out'] * g['distance']).groupby(['trtype','t']).sum() + ## (r,t) + dfs['cap_converter_out'] = g['CAP_CONVERTER'].Level + ## (r,rr,h,trtype,t) (r,rr,allh,t,trtype) + dfs['tran_flow_all_rep'] = g['FLOW'].Level.loc[:,:,g['h_rep'].index] + ## (r,rr,allh,trtype,t) + dfs['tran_flow_all_stress'] = ( + g['FLOW'].Level.reset_index() + .merge(g['h_stress_t'], left_on=['allh','t'], right_on=['allh','allt']) + .set_index(['r','rr','allh','trtype','t']).Level + ) + ## (r,rr,h,trtype,t) (r,rr,allh,t,trtype) + dfs['tran_flow_rep'] = combine_forward_reverse(g['FLOW'], agg='net').loc[:,:,g['h_rep'].index] + ## (r,rr,allh,trtype,t) + dfs['tran_flow_stress'] = ( + combine_forward_reverse(g['FLOW']).reset_index() + .merge(g['h_stress_t'], left_on=['allh','t'], right_on=['allh','allt']) + .set_index(['r','rr','allh','trtype','t']).Level + ) + ## (r,rr,trtype,t) + dfs['tran_flow_rep_ann'] = (dfs['tran_flow_rep'] * g['hours']).groupby(['r','rr','trtype','t']).sum() + ## (r,rr,h,trtype,t) + dfs['tran_util_h_rep'] = dfs['tran_flow_all_rep'] / dfs['tran_cap_energy'] + ## (r,rr,allh,trtype,t) + dfs['tran_util_h_stress'] = dfs['tran_flow_all_stress'] / dfs['tran_cap_prm'] + ## (r,rr,trtype,t) + dfs['tran_util_ann_rep'] = ( + (dfs['tran_flow_all_rep'] * g['hours'] / dfs['tran_cap_energy']).groupby(['r','rr','trtype','t']).sum() + / g['hours'].sum() + ) + ## NOTE: We here assume that all solve years use the same total weighting for stress + ## periods and weight all stress hours the same, so we don't weight by the number of hours. + ## If we switch to different weightings for different stress periods, + ## should add an hours(allh,t) parameter. + ## (r,rr,trtype,t) + dfs['tran_util_ann_stress'] = dfs['tran_util_h_stress'].groupby(['r','rr','trtype','t']).mean() + return dfs + + +#%% Procedure +def main(case): + ## NOTE: If calculations slow down for large runs, consider dropping zeros upfront + ## in get_gams_results() to speed up processing + dictin = get_gams_results(case) + dictout = { + **calc_iq(dictin), + **calc_co2_stor(dictin), + **calc_transmission(dictin), + } + ## Drop zeros to reduce file size and match GAMS convention + for key, df in dictout.items(): + _df = df.rename('Value').reset_index() + dictout[key] = _df.loc[_df.Value != 0].dropna().copy() + return dictout diff --git a/reeds/core/terminus/report_data.csv b/reeds/core/terminus/report_data.csv new file mode 100644 index 00000000..b546c322 --- /dev/null +++ b/reeds/core/terminus/report_data.csv @@ -0,0 +1,24 @@ +*** +*** Sets and parameters +*** +distance +h_rep +h_stress_t +hours +r +routes +*** +*** Variables +*** +CAP_ABOVE_LIM +CAP_CONVERTER +CAPTRAN_ENERGY +CAPTRAN_GRP +CAPTRAN_PRM +CO2_CAPTURED +CO2_FLOW +CO2_SPURLINE_INV +CO2_STORED +CO2_TRANSPORT_INV +FLOW +INVTRAN diff --git a/reeds/core/terminus/report_dump.py b/reeds/core/terminus/report_dump.py index 7e3ef747..5169e539 100644 --- a/reeds/core/terminus/report_dump.py +++ b/reeds/core/terminus/report_dump.py @@ -12,6 +12,7 @@ from pathlib import Path sys.path.append(str(Path(__file__).parent.parent.parent.parent)) import reeds +from reeds.core.terminus import report_calcs #%% Generic functions @@ -70,7 +71,7 @@ def dfdict_to_h5( _symbol_list = dfdict.keys() if symbol_list is None else symbol_list ### Check for existing file - _filepath = filepath if filepath.endswith(".h5") else filepath + ".h5" + _filepath = filepath if Path(filepath).suffix == '.h5' else str(filepath) + '.h5' if os.path.exists(_filepath): if overwrite: os.remove(_filepath) @@ -270,9 +271,10 @@ def postprocess_outputs(case, outputs_path=None, verbose=0): # %%### Write results for each gdx file ### outputs gdx print("Loading outputs gdx") - dict_out = gdxpds.to_dataframes( - os.path.join(outputs_path, f"rep_{os.path.basename(case)}.gdx") - ) + dict_out = { + **gdxpds.to_dataframes(Path(outputs_path, f"rep_{os.path.basename(case)}.gdx")), + **report_calcs.main(case), + } print("Finished loading outputs gdx") write_dfdict( @@ -305,6 +307,6 @@ def postprocess_outputs(case, outputs_path=None, verbose=0): #%% All done print("Completed report_dump.py") try: - toc(tic=tic, year=0, path=case, process="report_dump.py") + reeds.log.toc(tic=tic, year=0, path=case, process="report_dump.py") except NameError: print("reeds/log.py not found, so not logging output") diff --git a/reeds/io.py b/reeds/io.py index a2730ea2..5b148b7e 100644 --- a/reeds/io.py +++ b/reeds/io.py @@ -496,7 +496,7 @@ def read_output( Returns: pd.DataFrame """ - if case.endswith('.h5'): + if Path(case).suffix == '.h5': h5path = case else: h5path = os.path.join(case, 'outputs', 'outputs.h5') diff --git a/reeds/results.py b/reeds/results.py index 6522feac..a67c1e8e 100644 --- a/reeds/results.py +++ b/reeds/results.py @@ -9,7 +9,10 @@ ### =========================================================================== import os import sys +import h5py +import numpy as np import pandas as pd +from pathlib import Path from itertools import product sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) @@ -710,4 +713,95 @@ def get_level_map(): return level_map -# %% + +def diff_outputs( + casebase:str|Path, + casecomp:str|Path, + outpath:str|Path|None|bool=None, + threshold_abs:float=0, + threshold_rel:float=0, + verbose:int=0, +) -> dict: + """ + Diff two {case}/outputs/outputs.h5 files and save the result to outpath. + + Args: + casebase: Absolute path to base ReEDS case + casecomp: Absolute path to comparison ReEDS case + outpath: Absolute path to resulting difference .h5 file + If None or True, difference file is saved to + {casebase}/comparisons/diff_{casebase.stem}_{casecomp.stem}.h5. + If False, no difference file is written. + threshold_abs: Absolute cutoff for differences to include + (i.e., to ignore differences of 0.1 MW in an output parameter in units of MW, + set to 0.1) + threshold_rel: [fraction] Relative cutoff for differences to include, relative to base + (i.e., to ignore differences <1%, set to 0.01) + verbose (int): If nonzero, print descriptive logs + + Returns: + dict of pd.DataFrames (keys = entries in outputs.h5 with nonzero diff) + + Inputs for testing: + casebase = Path(reeds.io.reeds_path, 'runs', 'v20260306_itlM0_Pacific') + casecomp = Path(reeds.io.reeds_path, 'runs', 'v20260306_itlM0_Pacific_CC') + outpath = None + """ + ## Check inputs + casebase = Path(casebase) + casecomp = Path(casecomp) + fpaths = { + 'base': Path(casebase, 'outputs', 'outputs.h5'), + 'comp': Path(casecomp, 'outputs', 'outputs.h5'), + } + for key, fpath in fpaths.items(): + if not fpath.is_file(): + raise FileNotFoundError(fpath) + ## Make output directory + if outpath in [None, True]: + outpath = Path( + casebase, 'outputs', 'comparisons', f'diff_{casebase.stem}_{casecomp.stem}.h5' + ) + elif outpath is False: + pass + else: + outpath = Path(outpath) + if outpath: + outpath.parent.mkdir(exist_ok=True, parents=True) + ## Get results + keys = {} + dictin = {} + for case, fpath in fpaths.items(): + _dictin = {} + with h5py.File(fpath, 'r') as f: + keys[case] = list(f) + for key in keys[case]: + df = reeds.io.read_output(fpath, key) + indices = [i for i in df if i.lower() != 'value'] + if len(indices): + df = df.set_index(indices) + _dictin[key] = df.squeeze(1).astype(np.float32) + dictin[case] = _dictin + ## Take diff + allkeys = sorted(set(keys['base'] + keys['comp'])) + dictout = {} + for key in allkeys: + df = pd.concat({case:dictin[case].get(key, None) for case in dictin}, axis=1) + ## Sets have no value column; ignore them + if df.shape[1] == 0: + continue + df['diff'] = df.get('comp', 0) - df.get('base', 0) + df['frac'] = df.get('comp', 0) / df.get('base', 0) - 1 + df = df.loc[(df['diff'].abs() > threshold_abs) & (df['frac'].abs() > threshold_rel)].copy() + if len(df): + dictout[key] = df + if (verbose > 1) or (((verbose > 0) and len(df))): + print(f'{key}: {len(df)} differences') + ## Write it + if outpath: + if outpath.is_file(): + outpath.unlink() + for key, df in dictout.items(): + reeds.io.write_to_h5(df.reset_index(), key, outpath) + print(f'Difference written to {outpath}') + return dictout