From 436ade06d57c566aa1e90419b33b22492b804b36 Mon Sep 17 00:00:00 2001 From: bsergi Date: Wed, 15 Apr 2026 14:35:39 -0400 Subject: [PATCH 01/15] branch transfer finito linkage --- ReEDS_Augur/prep_data.py | 32 +++- b_inputs.gms | 50 ++++-- c_supplymodel.gms | 88 ++++++++-- c_supplyobjective.gms | 49 +++--- cases.csv | 6 + cases_finito.csv | 46 +++++ createmodel.gms | 14 ++ d1_temporal_params.gms | 13 +- d2_varfix.gms | 6 + d3_data_dump.gms | 18 +- d_solveoneyear.gms | 19 +++ e_report.gms | 141 ++++++++++----- e_report_dump.py | 27 ++- e_report_params.csv | 1 + input_processing/copy_files.py | 5 + input_processing/hourly_load.py | 47 +++++ input_processing/hourly_repperiods.py | 27 ++- input_processing/hourly_writetimeseries.py | 3 + inputs/emission_constraints/co2_tax.csv | 84 ++++----- reeds/reedsplots.py | 1 + restart_runs.py | 14 +- runbatch.py | 189 ++++++++++++++++++++- 22 files changed, 716 insertions(+), 164 deletions(-) create mode 100644 cases_finito.csv diff --git a/ReEDS_Augur/prep_data.py b/ReEDS_Augur/prep_data.py index 8f4febfc..9805767d 100644 --- a/ReEDS_Augur/prep_data.py +++ b/ReEDS_Augur/prep_data.py @@ -147,8 +147,38 @@ def main(t, casedir, iteration=0): h_dt_szn.index.map(hmap_allyrs.set_index(['year', 'hour'])['*timestamp'])) h_dt_szn = h_dt_szn.reset_index().set_index('timestamp') + # load exogenous demand seen by ReEDS load = reeds.io.read_file(os.path.join(inputs_case, 'load.h5'), parse_timestamps=True) + load_year = load.loc[t] + + # when running linked model with FINITO (GSw_FINITO_Link=1), + # we also need to add in the load from FINITO + if int(sw.GSw_FINITO_Link): + load_finito_rt = gdxreeds['load_finito_rt'].rename(columns={'allh':'h', 'Value':'load_MW'}) + + # down select to relevant model year + load_finito = load_finito_rt.loc[load_finito_rt.t==t].drop('t', axis=1).copy() + #TODO: why so much shifting? validate the dynamics here with FINITO team + #TODO: also check against reference quantity + + # map from rep day to actual hour + # since we don't have multi-year profiles for FINITO load + # we assume they repeat across all weather years + # TODO: test that this works with multiple regions + load_finito = pd.merge(load_finito, h_dt_szn.reset_index(), on='h', how='outer')[['timestamp', 'r', 'load_MW']] + load_finito = load_finito.rename(columns={'timestamp':'datetime'}) + # add model year back in + #load_finito = load_finito.assign(year=t) + load_finito = load_finito.pivot(index=['datetime'], columns='r', values='load_MW') + + # convert timezone and fill any missing columns + load_finito.index = load_finito.index.tz_convert(load_year.index.tz) + load_finito = load_finito.reindex(columns=load_year.columns, fill_value=0) + + # add to load + load_new = load_year + load_finito + resources = pd.read_csv(os.path.join(inputs_case, 'resources.csv')) recf = reeds.io.read_file(os.path.join(inputs_case, 'recf.h5'), parse_timestamps=True) recf.columns = pd.MultiIndex.from_tuples([tuple(x.split('|')) for x in recf.columns], @@ -367,7 +397,7 @@ def intify(v): .merge(h_dt_szn_load_years[['h']].reset_index(), left_on='allh', right_on='h') .pivot(index='timestamp', columns='r', values='Value') ) - load_year = load.loc[t].add(can_exports, fill_value=0) + load_year = load_year.add(can_exports, fill_value=0) ### PRAS doesn't yet handle flexible load, so include all H2/DAC load in the ### version we write for PRAS diff --git a/b_inputs.gms b/b_inputs.gms index 0d59e562..3c135f24 100644 --- a/b_inputs.gms +++ b/b_inputs.gms @@ -1341,7 +1341,8 @@ set tmodel(t) "years to include in the model", tfix(t) "years to fix variables over when summing over previous years", tprev(t,tt) "previous modeled tt from year t", stfeas(st) "states to include in the model", - tsolved(t) "years that have solved" ; + tsolved(t) "years that have solved", + tfuel(t) "years that use ReEDS fuel supply curve module (otherwise uses supply curves in FINITO)" ; *following parameters get re-defined when the solve years have been declared parameter mindiff(t) "minimum difference between t and all other tt that are in tmodel(t)" ; @@ -1354,7 +1355,7 @@ tfix(t) = no ; stfeas(st) = no ; tprev(t,tt) = no ; tsolved(t) = no ; - +tfuel(t)=no; *============================== * Year specification @@ -1379,6 +1380,12 @@ tprev(t,tt)$[tmodel_new(t)$tmodel_new(tt)$(tt.valmindiff(t))] = no ; +* If FINITO linkage is on, remove all modeled years from tfuel after first_year_finito +* as the FINITO supply curves will be used instead of those in ReEDS; +* otherwise, all modeled years use ReEDS supply curves and are eligible for tfuel +tfuel(t)$[tmodel_new(t)]=yes; +tfuel(t)$[tmodel_new(t)$Sw_FINITO_Link$(t.val>=%first_year_finito%)] = no ; + * In order to fill all necessary dimensions of upgrade techs parameters, we require * Sw_UpgradeYear in ban(i) to be a modeled year and thus we compute as either * the GSw_UpgradeYear option or the next modeled years after GSw_UpgradeYear @@ -2331,6 +2338,7 @@ noncumulative_prescriptions(pcat,r,t)$tmodel_new(t) ], prescribednonrsc(tt,pcat,r,"value") + prescribedrsc(tt,pcat,r,"value") } ; +noncumulative_prescriptions(pcat,r,t)$noncumulative_prescriptions(pcat,r,t) = round(noncumulative_prescriptions(pcat,r,t),6); parameter noncumulative_prescriptions_energy(pcat,r,t) "--MWh-- prescribed energy capacity that comes online in a given year" ; noncumulative_prescriptions_energy(pcat,r,t)$tmodel_new(t) @@ -2777,6 +2785,7 @@ h2_ptc(i,v,r,t)$valcap(i,v,r,t) = h2_ptc_in(i,v,t) ; * Otherwise, we assume it receives $0/kg because the cleanliness of its carbon cannot be proven h2_ptc("electrolyzer",v,r,t)$[(not Sw_H2_PTC)] = 0; + set h2_ptc_years(t) "years in which the hydrogen production incentive is active"; h2_ptc_years(t) = tmodel_new(t)$[sum{(i,v,r),h2_ptc(i,v,r,t)}]; @@ -3571,6 +3580,7 @@ $include inputs_case%ds%tsc_binwidth.csv $offdelim $onlisting / ; +tsc_binwidth(r,rr,tscbin) = 1e9; parameter tsc_forward(r,rr,tscbin) "--$/MW-- transmission upgrade cost for forward direction" / @@ -6543,22 +6553,25 @@ $offdelim $ondigit $onlisting / , - min_co2_spurline_distance "--mi-- minimum distance for a spur line (used to provide a floor for pipeline distances in r_cs_distance)" + min_co2_spurline_distance "--mi-- minimum distance for a spur line (used to provide a floor for pipeline distances in r_cs_distance)", + min_r_cs_distance(r) "--mi-- mininum euclidean distance between BA transmission endpoints and storage formations" ; $offempty -* Wherever BA centroids fall within formation boundaries, assume some average spur line distance to connect a CCS or DAC plant with an injection site -min_co2_spurline_distance = 20 ; -r_cs_distance(r,cs)$[r_cs_distance(r,cs) < min_co2_spurline_distance] = min_co2_spurline_distance ; - -* Assign spurline costs -cost_co2_spurline_cap(r,cs,t)$[r_cs(r,cs)$tmodel_new(t)] = Sw_CO2_spurline_cost * r_cs_distance(r,cs) ; +* Wherever BA centroids fall within formation boundaries, r_cs_distance will be equal to 0, to ensure these are not dropped from the set, set these elements to a small number +r_cs_distance(r,cs)$[r_cs(r,cs)$(r_cs_distance(r,cs) = 0)] = 1e-3 ; +* find the closest site to each region +min_r_cs_distance(r) = smin(cs$[r_cs(r,cs)], r_cs_distance(r,cs)); -* CO2 pipelines can be build between any two adjacent BAs -cost_co2_pipeline_cap(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_cost * pipeline_distance(r,rr) ; -cost_co2_pipeline_fom(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_fom * pipeline_distance(r,rr) ; +$ifthene.rcslimit %GSw_CO2_LimitStorageSites% == 1 +* remove the region site combinations that are not the closest +r_cs_distance(r,cs)$[r_cs(r,cs)$(r_cs_distance(r,cs) <= min_r_cs_distance(r))] = min_r_cs_distance(r) ; +r_cs_distance(r,cs)$[r_cs(r,cs)$(r_cs_distance(r,cs) > min_r_cs_distance(r))] = 0 ; +$endif.rcslimit -co2_routes(r,rr)$routes_adjacent(r,rr) = yes ; +* Wherever BA centroids fall within formation boundaries, assume some average spur line distance to connect a CCS or DAC plant with an injection site +min_co2_spurline_distance = 20 ; +r_cs_distance(r,cs)$[r_cs_distance(r,cs)$(r_cs_distance(r,cs) <= min_co2_spurline_distance )] = min_co2_spurline_distance ; $onempty table co2_char(cs,*) "co2 site characteristics including injection rate limit, total storage limit, and break even cost" @@ -6577,9 +6590,18 @@ cost_co2_stor_bec(cs,t) = co2_char(cs,"bec_%GSw_CO2_BEC%"); csfeas(cs)$[co2_storage_limit(cs)$co2_injection_limit(cs)] = yes ; * only want to consider r_cs pairs which have available capacity r_cs(r,cs)$[not csfeas(cs)] = no ; +r_cs(r,cs)$[not r_cs_distance(r,cs)] = no ; +* Assign spurline costs +cost_co2_spurline_cap(r,cs,t)$[r_cs(r,cs)$tmodel_new(t)] = Sw_CO2_spurline_cost * r_cs_distance(r,cs) ; cost_co2_spurline_fom(r,cs,t)$[r_cs(r,cs)$tmodel_new(t)] = Sw_CO2_spurline_fom * r_cs_distance(r,cs) ; +* CO2 pipelines can be build between any two adjacent BAs +cost_co2_pipeline_cap(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_cost * pipeline_distance(r,rr) ; +cost_co2_pipeline_fom(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_fom * pipeline_distance(r,rr) ; + +co2_routes(r,rr)$routes_adjacent(r,rr) = yes ; + cost_co2_pipeline_cap(r,rr,t) = %GSw_CO2_CostAdj% * cost_co2_pipeline_cap(r,rr,t); cost_co2_pipeline_fom(r,rr,t) = %GSw_CO2_CostAdj% * cost_co2_pipeline_fom(r,rr,t); cost_co2_stor_bec(cs,t) = %GSw_CO2_CostAdj% * cost_co2_stor_bec(cs,t) ; @@ -6771,4 +6793,4 @@ emit_rate(etype,e,i,v,r,t)$[not valgen(i,v,r,t)] = 0 ; *============================================================ valinv_init(i,v,r,t) = valinv(i,v,r,t) ; -valcap_init(i,v,r,t) = valcap(i,v,r,t) ; +valcap_init(i,v,r,t) = valcap(i,v,r,t) ; \ No newline at end of file diff --git a/c_supplymodel.gms b/c_supplymodel.gms index 029d53f1..71389843 100644 --- a/c_supplymodel.gms +++ b/c_supplymodel.gms @@ -396,6 +396,14 @@ eq_loadcon(r,h,t)$tmodel(t).. * (the effect is the same but avoiding the h-indexed OP_LOADSITE reduces solve time) + OP_LOADSITE(r,h,t)$[Sw_LoadSiteCF$(Sw_LoadSiteCF<1)$val_loadsite(r)] + CAP_LOADSITE(r,t)$[(Sw_LoadSiteCF=1)$val_loadsite(r)] + +* [plus] load for industrial and converted fuel facilities (FINITO), +* including the PRM for stress periods +* TODO: should this include distloss? +* [MWh/hr = MW] +$ifthene.linked_load Sw_FINITO_Link==1 + + USE_ELE_FINITO(r,h,t)$[t_finito(t)] * (1 + prm(r,t)$h_stress(h)) +$endif.linked_load ; * --------------------------------------------------------------------------- @@ -1040,7 +1048,6 @@ eq_rsc_INVlim(r,i,rscbin,t)$[tmodel(t) *plus exogenous (pre-start-year) capacity, using its level in the first year (tfirst) + sum{(ii,v,tt)$[tfirst(tt)$rsc_agg(i,ii)$exog_rsc(i)], capacity_exog_rsc(ii,v,r,rscbin,tt) } - ; * --------------------------------------------------------------------------- @@ -2881,7 +2888,7 @@ eq_national_gen(t)$[tmodel(t)$national_gen_frac(t)$Sw_GenMandate].. * --------------------------------------------------------------------------- *gas used from each bin is the sum of all gas used -eq_gasused(cendiv,h,t)$[tmodel(t)$((Sw_GasCurve=0) or (Sw_GasCurve=3))].. +eq_gasused(cendiv,h,t)$[tmodel(t)$tfuel(t)$((Sw_GasCurve=0) or (Sw_GasCurve=3))].. sum{gb,GASUSED(cendiv,gb,h,t) } @@ -2898,7 +2905,7 @@ eq_gasused(cendiv,h,t)$[tmodel(t)$((Sw_GasCurve=0) or (Sw_GasCurve=3))].. * --------------------------------------------------------------------------- * gas from each bin needs to less than its capacity -eq_gasbinlimit(cendiv,gb,t)$[tmodel(t)$(Sw_GasCurve=0)].. +eq_gasbinlimit(cendiv,gb,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=0)].. gaslimit(cendiv,gb,t) @@ -2909,7 +2916,7 @@ eq_gasbinlimit(cendiv,gb,t)$[tmodel(t)$(Sw_GasCurve=0)].. * --------------------------------------------------------------------------- -eq_gasbinlimit_nat(gb,t)$[tmodel(t)$(Sw_GasCurve=3)].. +eq_gasbinlimit_nat(gb,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=3)].. gaslimit_nat(gb,t) @@ -2922,7 +2929,7 @@ eq_gasbinlimit_nat(gb,t)$[tmodel(t)$(Sw_GasCurve=3)].. * --------------------------------------------------------------------------- -eq_gasaccounting_regional(cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)].. +eq_gasaccounting_regional(cendiv,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)].. sum{fuelbin, VGASBINQ_REGIONAL(fuelbin,cendiv,t) } @@ -2935,7 +2942,7 @@ eq_gasaccounting_regional(cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)].. * --------------------------------------------------------------------------- -eq_gasaccounting_national(t)$[tmodel(t)$(Sw_GasCurve=1)].. +eq_gasaccounting_national(t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)].. sum{fuelbin,VGASBINQ_NATIONAL(fuelbin,t) } @@ -2948,7 +2955,7 @@ eq_gasaccounting_national(t)$[tmodel(t)$(Sw_GasCurve=1)].. * --------------------------------------------------------------------------- -eq_gasbinlimit_regional(fuelbin,cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)].. +eq_gasbinlimit_regional(fuelbin,cendiv,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)].. Gasbinwidth_regional(fuelbin,cendiv,t) @@ -2959,7 +2966,7 @@ eq_gasbinlimit_regional(fuelbin,cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)].. * --------------------------------------------------------------------------- -eq_gasbinlimit_national(fuelbin,t)$[tmodel(t)$(Sw_GasCurve=1)].. +eq_gasbinlimit_national(fuelbin,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)].. Gasbinwidth_national(fuelbin,t) @@ -2973,10 +2980,10 @@ eq_gasbinlimit_national(fuelbin,t)$[tmodel(t)$(Sw_GasCurve=1)].. *============================== * -- Bioenergy Supply Curve -- *============================== +* defer to FINITO representation when models are linked * --------------------------------------------------------------------------- - -eq_bioused(r,t)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }$tmodel(t)].. +eq_bioused(r,t)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }$tmodel(t)$tfuel(t)].. sum{bioclass, BIOUSED(bioclass,r,t) } @@ -2995,7 +3002,7 @@ eq_bioused(r,t)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }$tmodel(t)].. * --------------------------------------------------------------------------- * biomass consumption limit is annual -eq_biousedlimit(bioclass,usda_region,t)$tmodel(t).. +eq_biousedlimit(bioclass,usda_region,t)$[tmodel(t)$tfuel(t)].. biosupply(usda_region,bioclass,"cap") @@ -3565,6 +3572,11 @@ eq_h2_demand(p,t)$[(sameas(p,"H2"))$tmodel(t)$(yeart(t)>=h2_demand_start)$(Sw_H2 + sum{(i,v,r,h)$[valgen(i,v,r,t)$h2_combustion(i)$h_rep(h)], GEN(i,v,r,h,t) * hours(h) * h2_combustion_intensity * heat_rate(i,v,r,t) } + +* hydrogen demand from indusrty (FINITO): FINITO demand [MMBtu/yr] * conversion [metric tons-H2/MMBtu-H2] +$ifthene.linked_h2_nat Sw_FINITO_Link==1 + + [sum{(r,h)$h_rep(h), hours(h) * USE_H2_FINITO(r,h,t) * h2_metric_tons_per_mmbtu }]$t_finito(t) +$endif.linked_h2_nat ; * --------------------------------------------------------------------------- @@ -3596,6 +3608,11 @@ eq_h2_demand_regional(r,h,t) + sum{(i,v)$[valgen(i,v,r,t)$h2_combustion(i)], GEN(i,v,r,h,t) * h2_combustion_intensity * heat_rate(i,v,r,t) } + +* regional hydrogen demand for industry (FINITO) +$ifthene.linked_h2_reg Sw_FINITO_Link==1 + + [ USE_H2_FINITO(r,h,t) * h2_metric_tons_per_mmbtu ]$t_finito(t) +$endif.linked_h2_reg ; * --------------------------------------------------------------------------- @@ -3838,11 +3855,16 @@ eq_co2_capture(r,h,t) * capture from DAC + sum{(i,v)$[dac(i)$valcap(i,v,r,t)$i_p(i,"DAC")], PRODUCE("DAC",i,v,r,h,t) }$Sw_DAC +* (linked) capture from industry [metric_tons-CO2/hr]: +* calculation: hours_per_year [yrs/hr] * capture [scaled_metric_tons-CO2/yr] / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2] +$ifthene.linked_co2_capture Sw_FINITO_Link==1 + + [CAPTURE_CO2EM(r,h,t) / co2_scale]$[t_finito(t)$h_rep(h)] +$endif.linked_co2_capture ; * --------------------------------------------------------------------------- -eq_co2_transport_caplimit(r,rr,h,t)$[co2_routes(r,rr)$Sw_CO2_Detail +eq_co2_transport_caplimit(r,rr,h,t)$[co2_routes(r,rr)$Sw_CO2_Detail$h_rep(h) $tmodel(t)$(yeart(t)>=co2_detail_startyr)].. *capacity computed as cumulative investments of co2 pipelines up to the current year @@ -3858,7 +3880,8 @@ eq_co2_transport_caplimit(r,rr,h,t)$[co2_routes(r,rr)$Sw_CO2_Detail * --------------------------------------------------------------------------- -eq_co2_spurline_caplimit(r,cs,h,t)$[Sw_CO2_Detail$r_cs(r,cs)$tmodel(t)$(yeart(t)>=co2_detail_startyr)].. +eq_co2_spurline_caplimit(r,cs,h,t)$[Sw_CO2_Detail$r_cs(r,cs)$h_rep(h)$tmodel(t) + $(yeart(t)>=co2_detail_startyr)].. *capacity computed as cumulative investments of co2 spurlines up to the current year sum{tt$[(yeart(tt)<=yeart(t))$(tmodel(tt) or tfix(tt))$(yeart(tt)>=co2_detail_startyr)], @@ -3867,11 +3890,16 @@ eq_co2_spurline_caplimit(r,cs,h,t)$[Sw_CO2_Detail$r_cs(r,cs)$tmodel(t)$(yeart(t) =g= CO2_STORED(r,cs,h,t) +* (linked) extraction of CO2 [metric_tons-CO2/hr] +* calculation: hours_per_year [yrs/hr] * (1 / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2]) * use [scaled_metric_tons-CO2/yr] +$ifthene.linked_co2_spurline_caplimit Sw_FINITO_Link==1 + + [ (1 / co2_scale) * EXTRACT_CO2_CS(cs,r,h,t) ]$[h_rep(h)$t_finito(t)] +$endif.linked_co2_spurline_caplimit ; * --------------------------------------------------------------------------- -eq_co2_sink(r,h,t)$[tmodel(t)$Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)].. +eq_co2_sink(r,h,t)$[tmodel(t)$Sw_CO2_Detail$h_rep(h)$(yeart(t)>=co2_detail_startyr)].. *the amount of co2 stored from r in all of its cs sites sum{cs$r_cs(r,cs), CO2_STORED(r,cs,h,t) } @@ -3886,11 +3914,24 @@ eq_co2_sink(r,h,t)$[tmodel(t)$Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)].. * net trade + sum{rr$co2_routes(r,rr), CO2_FLOW(rr,r,h,t) - CO2_FLOW(r,rr,h,t) } + +* (linked) extraction - use of CO2 [metric_tons-CO2/hr] +*calculation: hours_per_year [yrs/hr] * (1 / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2]) * extraction/use [scaled_metric_tons-CO2/yr] +$ifthene.linked_co2_sink Sw_FINITO_Link==1 + + (1 / co2_scale) * [ +* extraction from all cs sites in r + + sum{cs$[csfeas(cs)$r_cs(r,cs)], EXTRACT_CO2_CS(cs,r,h,t)} +* total use of CO2, equivalent to supply + - USE_CO2(r,h,t) + ]$[h_rep(h)$t_finito(t)] +$endif.linked_co2_sink + ; * --------------------------------------------------------------------------- -eq_co2_injection_limit(cs,h,t)$[Sw_CO2_Detail$tmodel(t)$(yeart(t)>=co2_detail_startyr)$csfeas(cs)].. +eq_co2_injection_limit(cs,h,t)$[tmodel(t)$Sw_CO2_Detail$h_rep(h) + $(yeart(t)>=co2_detail_startyr)$csfeas(cs)].. * exogenously defined injection limit co2_injection_limit(cs) @@ -3899,6 +3940,14 @@ eq_co2_injection_limit(cs,h,t)$[Sw_CO2_Detail$tmodel(t)$(yeart(t)>=co2_detail_st * must exceed metric tons per hour entering storage sum{r$r_cs(r,cs), CO2_STORED(r,cs,h,t) } + + + +* (linked) extraction of CO2 for use [metric_tons-CO2/hr] +* calculation: hours_per_year [yrs/hr] * (1 / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2]) * use [scaled_metric_tons-CO2/yr] +$ifthene.linked_co2_injection_limit Sw_FINITO_Link==1 + + (1 / co2_scale) * sum{r$[r_cs(r,cs)], EXTRACT_CO2_CS(cs,r,h,t) }$[h_rep(h)$t_finito(t)] +$endif.linked_co2_injection_limit ; * --------------------------------------------------------------------------- @@ -3915,7 +3964,16 @@ eq_co2_cumul_limit(cs,t)$[tmodel(t)$Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr) $[(yeart(tt)<=yeart(t))$(tmodel(tt) or tfix(tt))$(yeart(tt)>=co2_detail_startyr) $r_cs(r,cs)$h_rep(h)], yearweight(tt) * hours(h) * CO2_STORED(r,cs,h,tt) } + +* (linked) cumulative amount extracted over time +$ifthene.linked_co2_storage_cumul_limit Sw_FINITO_Link==1 + - sum{(r,h,tt) + $[(yeart(tt)<=yeart(t))$(tmodel(tt) or tfix(tt))$(yeart(tt)>=co2_detail_startyr) + $r_cs(r,cs)$(tfuel(tt))], + yearweight(tt) * hours(h) * EXTRACT_CO2_CS(cs,r,h,tt) }$[t_finito(t)] +$endif.linked_co2_storage_cumul_limit ; + * --------------------------------------------------------------------------- *=================== diff --git a/c_supplyobjective.gms b/c_supplyobjective.gms index 38a5b0b3..a5955399 100644 --- a/c_supplyobjective.gms +++ b/c_supplyobjective.gms @@ -19,7 +19,17 @@ Variable Z "--$-- total cost of operations and investment, scale varie * objective function is the sum over modeled years of the investment * and operations components -eq_ObjFn.. Z =e= cost_scale * sum{t$tmodel(t), Z_inv(t) + Z_op(t) } ; +eq_ObjFn.. Z =e= cost_scale * ( +* electricity and H2 costs + sum{t$tmodel(t), Z_inv(t) + Z_op(t) } +* economy-wide costs from FINITO: deflate from $2018 to $2004 +* and remove any FINITO scaling +$ifthene.linked_objective Sw_FINITO_Link==1 + + deflator('2018') / obj_scale + * sum{t$[t_finito(t)], Z_finito(t) } +$endif.linked_objective + ) +; *======================================================= * -- Investment component of the objective function -- @@ -221,7 +231,7 @@ eq_Objfn_op(t)$tmodel(t).. * plus cost of H2 fuel when using fixed price (Sw_H2=0) or during stress periods. * When using endogenous H2 price (Sw_H2=1 or Sw_H2=2), H2 fuel cost is captured elsewhere * via the capex + opex costs of H2 production and its associated electricity demand. - + sum{(i,v,r,h)$[valgen(i,v,r,t)$heat_rate(i,v,r,t) + + sum{(i,v,r,h)$[tfuel(t)$valgen(i,v,r,t)$heat_rate(i,v,r,t) $(not gas(i))$(not bio(i))$(not cofire(i)) $((not h2_combustion(i)) or h2_combustion(i)$[(Sw_H2=0) or h_stress(h)])], hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN(i,v,r,h,t) } @@ -232,50 +242,49 @@ eq_Objfn_op(t)$tmodel(t).. * --cofire coal consumption--- * cofire bio consumption already accounted for in accounting of BIOUSED - + sum{(i,v,r,h)$[valgen(i,v,r,t)$cofire(i)$heat_rate(i,v,r,t)], + + sum{(i,v,r,h)$[tfuel(t)$valgen(i,v,r,t)$cofire(i)$heat_rate(i,v,r,t)], (1-bio_cofire_perc) * hours(h) * heat_rate(i,v,r,t) * fuel_price("coal-new",r,t) * GEN(i,v,r,h,t) } * --- cost of natural gas--- *Sw_GasCurve = 2 (static natural gas prices) *first - gas consumed for electricity generation - + sum{(i,v,r,h)$[valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)$(Sw_GasCurve = 2)], + + sum{(i,v,r,h)$[tfuel(t)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)$(Sw_GasCurve = 2)], hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN(i,v,r,h,t) } *second - gas consumed by gas-powered DAC - + sum{(v,r,h)$[valcap("dac_gas",v,r,t)$(Sw_GasCurve = 2)], - hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE("DAC","dac_gas",v,r,h,t) }$Sw_DAC_Gas + + sum{(v,r,h)$[tfuel(t)$valcap("dac_gas",v,r,t)$(Sw_GasCurve = 2)], + hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE("DAC","dac_gas",v,r,h,t) }$[Sw_DAC_Gas] *Sw_GasCurve = 0 (census division supply curves natural gas prices) - + sum{(cendiv,gb), sum{h, hours(h) * GASUSED(cendiv,gb,h,t) } + + sum{(cendiv,gb)$[tfuel(t)], sum{h, hours(h) * GASUSED(cendiv,gb,h,t) } * gasprice(cendiv,gb,t) - }$(Sw_GasCurve = 0) + }$[(Sw_GasCurve = 0)] *Sw_GasCurve = 3 (national supply curve for natural gas prices with census division multipliers) - + sum{(h,cendiv,gb), hours(h) * GASUSED(cendiv,gb,h,t) + + sum{(h,cendiv,gb)$[tfuel(t)], hours(h) * GASUSED(cendiv,gb,h,t) * gasadder_cd(cendiv,t,h) + gasprice_nat_bin(gb,t) - }$(Sw_GasCurve = 3) + }$[(Sw_GasCurve = 3)] *Sw_GasCurve = 1 (national and census division supply curves for natural gas prices) *first - anticipated costs of gas consumption given last year's amount - + (sum{(i,r,v,cendiv,h)$[valgen(i,v,r,t)$gas(i)], + + (sum{(i,r,v,cendiv,h)$[valgen(i,v,r,t)$gas(i)$tfuel(t)], gasmultterm(cendiv,t) * szn_adj_gas(h) * cendiv_weights(r,cendiv) * hours(h) * heat_rate(i,v,r,t) * GEN(i,v,r,h,t) } *second - adjustments based on changes from last year's consumption at the regional and national level - + sum{(fuelbin,cendiv), + + sum{(fuelbin,cendiv)$tfuel(t), gasbinp_regional(fuelbin,cendiv,t) * VGASBINQ_REGIONAL(fuelbin,cendiv,t) } - + sum{(fuelbin), + + sum{(fuelbin)$tfuel(t), gasbinp_national(fuelbin,t) * VGASBINQ_NATIONAL(fuelbin,t) } - )$[Sw_GasCurve = 1] + )$[(Sw_GasCurve = 1)] * ---cost of biofuel consumption and biomass transport--- - + sum{(r,bioclass)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }], + + sum{(r,bioclass)$[tfuel(t)$sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }], BIOUSED(bioclass,r,t) * (sum{usda_region$r_usda(r,usda_region), biosupply(usda_region, bioclass, "price") } + bio_transport_cost) } - * --- hurdle costs for transmission flow --- + sum{(r,rr,h,trtype)$[routes(r,rr,trtype,t)$cost_hurdle(r,rr,t)], cost_hurdle(r,rr,t) * FLOW(r,rr,h,t,trtype) * hours(h) } @@ -340,15 +349,15 @@ eq_Objfn_op(t)$tmodel(t).. CO2_SPURLINE_INV(r,cs,tt) } }$[Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)] * --- CO2 injection break even costs - + sum{(r,cs,h)$r_cs(r,cs), hours(h) * CO2_STORED(r,cs,h,t) * cost_co2_stor_bec(cs,t) }$[Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)] + + sum{(r,cs,h)$[r_cs(r,cs)$h_rep(h)], hours(h) * CO2_STORED(r,cs,h,t) * cost_co2_stor_bec(cs,t) }$[Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)] * --- Tax credit for CO2 stored --- * note conversion to 12-year CRF given length of CO2 captured incentive payments - - sum{(i,v,r,h)$[valgen(i,v,r,t)$co2_captured_incentive(i,v,r,t)], + - sum{(i,v,r,h)$[valgen(i,v,r,t)$h_rep(h)$co2_captured_incentive(i,v,r,t)], (crf(t) / crf_co2_incentive(t)) * co2_captured_incentive(i,v,r,t) * hours(h) * capture_rate("CO2",i,v,r,t) * GEN(i,v,r,h,t)} * --- Tax credit for CO2 stored for DAC --- - - sum{(p,i,v,r,h)$[dac(i)$valcap(i,v,r,t)$i_p(i,p)$h_rep(h)], + - sum{(p,i,v,r,h)$[dac(i)$valcap(i,v,r,t)$i_p(i,p)$h_rep(h)$co2_captured_incentive(i,v,r,t)], (crf(t) / crf_co2_incentive(t)) * co2_captured_incentive(i,v,r,t) * hours(h) * PRODUCE(p,i,v,r,h,t)} * --- PTC value for electric power generation --- @@ -362,7 +371,7 @@ eq_Objfn_op(t)$tmodel(t).. - sum{(p,v,r,h)$[valcap("electrolyzer",v,r,t)$(sameas(p,"H2"))$h2_ptc("electrolyzer",v,r,t)$h_rep(h)], hours(h) * PRODUCE(p,"electrolyzer",v,r,h,t) * (crf(t) / crf_h2_incentive(t)) * h2_ptc("electrolyzer",v,r,t) * 1e3} - $[Sw_H2_PTC$Sw_H2$h2_ptc_years(t)$(yeart(t) >= h2_demand_start)] + $[(Sw_H2_PTC)$Sw_H2$h2_ptc_years(t)$(yeart(t) >= h2_demand_start)] *end multiplier for pvf_onm ) diff --git a/cases.csv b/cases.csv index 9d0752e7..59f86c65 100644 --- a/cases.csv +++ b/cases.csv @@ -33,6 +33,11 @@ coalscen,"Coal price scenario (inputs\fuelprices\{coalscen}.csv). The options fo demandscen,Demand scenario (see inputs\load),N/A,AEO_2025_reference, dr_shedscen,DR shed scenario,demo_data_January_2025,demo_data_January_2025, evmcscen,EVMC scenario for load and supply curve and flexibility,Baseline,Baseline, +GSw_FINITO_Link,Switch to enable linkage with FINITO,0; 1,0, +finito_dir,Directory with cloned FINITO repository,N/A,/projects/finito/molmezt/FINITO, +finito_cases_file,"extension of the FINITO cases file including all scenarios to be run, e.g., 'linked' for cases_linked.csv",N/A,linked, +finito_case,Case name for FINITO in {finito_dir}/cases_linked.csv; use 'same' if the scenario name in cases_linked.csv matches the scenario name in the cases file,N/A,BAU, +first_year_finito,First year when the FINITO module is enabled,N/A,2018, geodiscov,Annual discover rates for new geothermal sites,BAU; TI,BAU, geohydrosupplycurve,Geohydro Supply Curve,ATB_2023,ATB_2023, egssupplycurve,EGS Supply Curve,ATB_2023; reV,reV, @@ -116,6 +121,7 @@ GSw_CO2_pipeline_fom,Set FOM cost for CO2 Pipeline. 2004$/(tonne-mi)/hr-yr. Sugg GSw_CO2_spurline_cost,Set capex cost for CO2 Spurline. 2004$/(tonne-mi)/hr. Suggested value: 2257,N/A,2257, GSw_CO2_spurline_fom,Set FOM cost for CO2 Spurline. 2004$/(tonne-mi)/hr-yr. Suggested value: 60,N/A,60, GSw_CO2_Detail,Switch to enable detailed representation CO2 transportation and storage,0; 1,0, +GSw_CO2_LimitStorageSites,Switch to limit storage sites to the closest site for each region,0; 1,0, GSw_CO2_CostAdj,multiplier for co2 infrastructure when Sw_CO2_Detail = 1,float,1, GSw_CO2_BEC,"Break even cost capacity factor assignment. Must be a suffix to ""bec_"" in the columns of co2_site_char.csv",N/A,70, GSw_CoalIGCC,Turn on/off COAL IGCC,0; 1,1, diff --git a/cases_finito.csv b/cases_finito.csv new file mode 100644 index 00000000..25fb39db --- /dev/null +++ b/cases_finito.csv @@ -0,0 +1,46 @@ +,Default Value,ReEDSonly,BAUreg,Reference,Reference-AEO,Reference-HOG,TradeReference,TradeReference-AEO,Trade-HOG,Reference-LowBio,Reference-HighBio,Tax200,Tax200-AEO,NetZero,NetZero-LowBio,NetZero-HighBio,NetZero-SlackLim1500,NetZero-SlackLim1500-CO2Stor-Mid,NetZero-SlackLim1500-CO2Stor-Low,NetZero-SlackLim500-CO2Stor-Low,NetZero-SlackCost2500-CO2Stor-Low,aggreg_54_level +ignore,1,1,0,1,1,1,1,1,1,,,,,,,,,,,,, +GSw_FINITO_Link,1,0,,,,,,,,,,,,,,,,,,,, +finito_dir,/path/to/FINITO/directory,0,,,,,,,,,,,,,,,,,,,, +finito_cases_file,linked,,,,,,,,,,,,,,,,,,,,, +finito_case,BAU,,BAUreg,,BAU-AEO,AEO-HOG,TradeReEDS,TradeReEDS-AEO,TradeReEDS-HOG,BAU-LowBio,BAU-HighBio,,BAU-AEO,NetZero,NetZero-LowBio,NetZero-HighBio,NetZero-SlackLim1500,NetZero-SlackLim1500-CO2Stor-Mid,NetZero-SlackLim1500-CO2Stor-Low,NetZero-SlackLim500-CO2Stor-Low,NetZero-SlackCost2500-CO2Stor-Low,aggreg_54_level +coalscen,AEO_2025_reference,,,,,,,,,,,,,,,,,,,,, +ngscen,AEO_2025_reference,,,,,AEO_2025_HOG,,,AEO_2025_HOG,,,,,,,,,,,,, +uraniumscen,AEO_2025_reference,,,,,,,,,,,,,,,,,,,,, +GSw_Region,country/USA,,st/CO,,,,,,,,,,,,,,,,,,, +GSw_GasCurve,1,,2,,,,,,,,,,,,,,,,,,, +GSw_gopt,2,,,,,,,,,,,,,,,,,,,,, +GSw_NewValCapShrink,1,,,,,,,,,,,,,,,,,,,,, +endyear,2050,,2035,,,,,,,,,,,,,,,,,,, +yearset,2010_2015_2020..2050..5,,,,,,,,,,,,,,,,,,,,, +first_year_finito,2018,,,,,,,,,,,,,,,,,,,,, +incentives_suffix,obbba,,,,,,,,,,,,,,,,,,,,, +GSw_Upgrades,0,,,,,,,,,,,,,,,,,,,,, +GSw_AnnualCap,0,,,,,,,,,,,,,,,,,,,,, +GSw_AnnualCapScen,start2023_100pct2035,,,,,,,,,,,,,,,,,,,,, +GSw_Clean_Air_Act,0,,,,,,,,,,,,,,,,,,,,, +GSw_CSP,0,,,,,,,,,,,,,,,,,,,,, +GSw_Geothermal,0,,,,,,,,,,,,,,,,,,,,, +GSw_H2,2,,,,,,,,,,,,,,,,,,,,, +GSw_H2_PTC,0,,,,,,,,,,,,,,,,,,,,, +GSw_H2_Demand_Case,BAU,,,,,,,,,,,,,,,,,,,,, +GSw_H2_SMR,1,,,,,,,,,,,,,,,,,,,,, +GSw_H2_Transport,0,,,,,,,,,,,,,,,,,,,,, +GSw_H2Combustion,1,,,,,,,,,,,,,,,,,,,,, +GSw_H2Combustionupgrade,1,,,,,,,,,,,,,,,,,,,,, +GSw_H2CombinedCycle,1,,,,,,,,,,,,,,,,,,,,, +GSw_DAC,1,,,,,,,,,,,,,,,,,,,,, +GSw_BECCS,1,,,,,,,,,,,,,,,,,,,,, +GSw_RegionResolution,ba,,,,,,,,,,,,,,,,,,,,,aggreg +GSw_TCPhaseout,1,,,,,,,,,,,,,,,,,,,,, +GSw_CO2_Detail,1,,,,,,,,,,,,,,,,,,,,, +GSw_CO2_LimitStorageSites,1,,,,,,,,,,,,,,,,,,,,, +numhintage,4,,,,,,,,,,,,,,,,,,,,, +GSw_CarbTax,0,,,,,,,,,,,1,1,,,,,,,,, +GSw_CarbTaxOption,default,,,,,,,,,,,t200_2050,t200_2050,,,,,,,,, +GSw_PRM_CapCredit,0,,,,,,,,,,,,,,,,,,,,, +GSw_PRM_StressIterateMax,0,,,,,,,,,,,,,,,,,,,,, +GSw_HourlyNumClusters,9,,,,,,,,,,,,,,,,,,,,, +GSw_HourlyChunkLengthRep,4,,,,,,,,,,,,,,,,,,,,, +GSw_HourlyChunkLengthStress,4,,,,,,,,,,,,,,,,,,,,, +GSw_HierarchyFile,default,,,,,,,,,,,,,,,,,,,,,agg54 diff --git a/createmodel.gms b/createmodel.gms index 36ab4603..d8ba1dca 100644 --- a/createmodel.gms +++ b/createmodel.gms @@ -1,5 +1,19 @@ $include b_inputs.gms + +$ifthene.finito_link %GSw_FINITO_Link% == 1 +$include finito/model/finito_input.gms +$endif.finito_link + +$ifthene.finito_link %GSw_FINITO_Link% == 1 +$include finito/model/finito_vars_eqs.gms +$endif.finito_link + $include c_supplymodel.gms + +$ifthene.finito_link %GSw_FINITO_Link% == 1 +$include finito/model/finito_model.gms +$endif.finito_link + $include c_supplyobjective.gms $include c_mga.gms $include d_solveprep.gms diff --git a/d1_temporal_params.gms b/d1_temporal_params.gms index 444c7fc5..fcbf8578 100644 --- a/d1_temporal_params.gms +++ b/d1_temporal_params.gms @@ -299,9 +299,6 @@ can_exports_h(r,h,t)$[hours(h)] = can_exports(r,t) * can_exports_h_frac(h) / hou $endif.Canada -* zero Canada exports when Canada is not modeled -can_imports_szn(r,szn,t)$[Sw_Canada=0] = 0 ; -can_exports_h(r,h,t)$[Sw_Canada=0] = 0 ; $onempty parameter canmexload(r,allh) "load for canadian and mexican regions" @@ -887,12 +884,20 @@ szn_adj_gas(h)$frac_h_quarter_weights(h,"wint") = szn_adj_gas(h) + frac_h_quarter_weights(h,"wint") * szn_adj_gas_winter ; +*============================================= +* ----- ReEDS-FINITO temporal parameters ----- +*============================================= +* remove focused sectors load from load_exog_static +$ifthene.linked_finito_temporal_params %GSw_FINITO_Link% ==1 +$include finito/model/linked_finito_temporal_params.gms +$endif.linked_finito_temporal_params + + *============================================= * -- Round parameters for GAMS -- *============================================= avail(i,r,h)$avail(i,r,h) = round(avail(i,r,h),3) ; -can_imports_szn(r,szn,t)$can_imports_szn(r,szn,t) = round(can_imports_szn(r,szn,t),3) ; can_exports_h(r,h,t)$can_exports_h(r,h,t) = round(can_exports_h(r,h,t),3) ; h_weight_csapr(h)$h_weight_csapr(h) = round(h_weight_csapr(h),3) ; load_exog(r,h,t)$load_exog(r,h,t) = round(load_exog(r,h,t),3) ; diff --git a/d2_varfix.gms b/d2_varfix.gms index 85db2442..3a523732 100644 --- a/d2_varfix.gms +++ b/d2_varfix.gms @@ -9,6 +9,10 @@ if(Sw_RemoveSmallNumbers = 1, INV_ENERGY.l(i,v,r,tfix)$[abs(INV_ENERGY.l(i,v,r,tfix)) < rhs_tolerance] = 0 ; INV_RSC.l(i,v,r,rscbin,tfix)$[abs(INV_RSC.l(i,v,r,rscbin,tfix)) < rhs_tolerance] = 0 ; INV_POI.l(r,tfix)$[abs(INV_POI.l(r,tfix)) < rhs_tolerance] = 0 ; + INV_REFURB.l(i,v,r,tfix)$[abs(INV_REFURB.l(i,v,r,tfix)) < rhs_tolerance] = 0; + CO2_SPURLINE_INV.l(r,cs,tfix)$[abs(CO2_SPURLINE_INV.l(r,cs,tfix)) < rhs_tolerance] = 0; + CO2_TRANSPORT_INV.l(r,rr,tfix)$[abs(CO2_TRANSPORT_INV.l(r,rr,tfix)) < rhs_tolerance] = 0; + H2_TRANSPORT_INV.l(r,rr,tfix)$[abs(H2_TRANSPORT_INV.l(r,rr,tfix)) < rhs_tolerance] = 0 ; H2_STOR_INV.l(h2_stor,r,tfix)$[abs(H2_STOR_INV.l(h2_stor,r,tfix)) < rhs_tolerance] = 0 ; H2_TRANSPORT_INV.l(r,rr,tfix) $[abs(H2_TRANSPORT_INV.l(r,rr,tfix) ) < rhs_tolerance] = 0 ; ); @@ -119,7 +123,9 @@ H2_STOR_LEVEL_SZN.fx(h2_stor,r,actualszn,tfix)$[(h2_stor_r(h2_stor,r))$(Sw_H2=2) *CO2-related variables CO2_CAPTURED.fx(r,h,tfix)$Sw_CO2_Detail = CO2_CAPTURED.l(r,h,tfix) ; +$ifthene.linked_varfix Sw_FINITO_Link==0 CO2_STORED.fx(r,cs,h,tfix)$[Sw_CO2_Detail$r_cs(r,cs)] = CO2_STORED.l(r,cs,h,tfix) ; +$endif.linked_varfix CO2_FLOW.fx(r,rr,h,tfix)$[Sw_CO2_Detail$co2_routes(r,rr)] = CO2_FLOW.l(r,rr,h,tfix) ; CO2_TRANSPORT_INV.fx(r,rr,tfix)$[Sw_CO2_Detail$co2_routes(r,rr)] = CO2_TRANSPORT_INV.l(r,rr,tfix) ; CO2_SPURLINE_INV.fx(r,cs,tfix)$[Sw_CO2_Detail$r_cs(r,cs)] = CO2_SPURLINE_INV.l(r,cs,tfix) ; diff --git a/d3_data_dump.gms b/d3_data_dump.gms index 29c97036..30a8ff0b 100644 --- a/d3_data_dump.gms +++ b/d3_data_dump.gms @@ -8,6 +8,7 @@ This file creates a gdx file with all of the data necessary for the Augur module - availability rates (1 - outage rates) - transmission capacities and loss rates - technology sets + - industrial electricity demand from FINITO when running linked model $offtext $if not set start_year $setglobal start_year %startyear% @@ -62,6 +63,7 @@ h2_usage_regional(r,allh,t) "--metric tons-- H2 usage by region" inv_cond_filt(i,v,t) "--set-- vintage-year mapping for investments by technology" inv_ivrt(i,v,r,t) "--MW-- investments in power generation capacity" inv_energy_ivrt(i,v,r,t) "--MWh-- investments in energy generation capacity" +load_finito_rt(r,allh,t) "--MWh-- total load from industrial sectors and converted fuels modeled by FINITO" m_cf_filt(i,v,r,allh) "--fraction-- capacity factor used in the model" m_cf_szn_filt(i,v,r,allszn) "--fraction-- modelled capacity factors filtered for hydro resources to set seasonal energy constraints" minloadfrac_filt(r,i,allszn) "--fraction-- modelled mingen fraction filtered for hydro resources to set mingen constraints" @@ -161,6 +163,10 @@ h2_usage_regional(r,h,t)$tcur(t) = h2_exogenous_demand_regional(r,'h2',h,t) + sum{(i,v)$[valgen(i,v,r,t)$h2_combustion(i)], GEN.l(i,v,r,h,t) * h2_combustion_intensity * heat_rate(i,v,r,t)} +* regional hydrogen demand for industry (FINITO) +$ifthene.linked_h2_reg_report Sw_FINITO_Link==1 + + [USE_H2_FINITO.l(r,h,t) * h2_metric_tons_per_mmbtu ]$t_finito(t) +$endif.linked_h2_reg_report ) ; @@ -284,15 +290,24 @@ cap_converter_filt(r) = sum{t$tcur(t), CAP_CONVERTER.l(r,t) } ; routes_filt(r,rr,trtype) = sum{t$tcur(t), routes(r,rr,trtype,t) } ; *============================ -* Flexible load data +* Load data *============================ +* flexible load flex_load(r,h) = sum{(flex_type,t)$tcur(t), load_exog_flex(flex_type,r,h,t) } ; flex_load_opt(r,h) = sum{(flex_type,t)$tcur(t), FLEX.l(flex_type,r,h,t) } ; ra_cap_loadsite(r,t)$[Sw_LoadSiteCF$val_loadsite(r)] = CAP_LOADSITE.l(r,t) ; +* FINITO load +$ifthene.linked_load Sw_FINITO_Link==1 + load_finito_rt(r,h,t) = USE_ELE_FINITO.l(r,h,t) ; +$else.linked_load + load_finito_rt(r,h,t) = 0 ; +$endif.linked_load + + *============================ * Extra consumption data *============================ @@ -382,6 +397,7 @@ execute_unload 'ReEDS_Augur%ds%augur_data%ds%reeds_data_%cur_year%.gdx' inv_ivrt inv_energy_ivrt ivt_num + load_finito_rt m_cf_filt m_cf_szn_filt maxage diff --git a/d_solveoneyear.gms b/d_solveoneyear.gms index b3484a80..a62eaaee 100644 --- a/d_solveoneyear.gms +++ b/d_solveoneyear.gms @@ -23,6 +23,13 @@ tload(t) = no ; tmodel(t) = no ; tmodel("%cur_year%") = yes ; +* (linked) also reset t_finito +$ifthene.linked_finito_time %GSw_FINITO_Link% == 1 +t_finito(t) = no ; +t_finito(t)$[tmodel(t)$(t.val>=first_year_finito)$(t.val<=endyear)] = yes ; +$endif.linked_finito_time + + $log 'Solving sequential case for...' $log ' Case: %case%' $log ' Year: %cur_year%' @@ -241,6 +248,7 @@ $endif.diagnose * ------------------------------ * Solve the Model * ------------------------------ +option limrow = 1e6; $ifthen.valstr %GSw_ValStr% == 1 OPTION lp = convert ; ReEDSmodel.optfile = 1 ; @@ -276,11 +284,22 @@ $endif.mga *** Adjust some parameters based on the solution for this solve year $include d2_post_solve_adjustments.gms +*** Do the same for finito if necessary +$ifthene.linked_finito_shrink %GSw_FINITO_Link% == 1 +* Industrial variable fix +$include finito/model/finito_post_solve_adjustments.gms +$endif.linked_finito_shrink *** Fix decision variables to their optimized levels for this solve year tfix("%cur_year%") = yes ; $include d2_varfix.gms +*** Fix FINITO decision variables to their optimized levels for this solve year +$ifthene.linked_finito_varfix %GSw_FINITO_Link% == 1 +* Industrial variable fix +$include finito/model/finito_varfix.gms +$endif.linked_finito_varfix + *** Dump data used in calculations between solve years $include d3_data_dump.gms diff --git a/e_report.gms b/e_report.gms index ebd14d10..f251e657 100644 --- a/e_report.gms +++ b/e_report.gms @@ -107,7 +107,7 @@ rev_cat "categories for renvenue streams" /load, res_marg, oper_res, rps, charge lcoe_cat "categories for LCOE calculation" /capcost, upgradecost, rsccost, fomcost, vomcost, gen / -loadtype "categories for types of load" / end_use, dist_loss, trans_loss, stor_charge, h2_prod, h2_network, dac / +loadtype "categories for types of load" / end_use, dist_loss, trans_loss, stor_charge, h2_prod, h2_network, dac, finito / h2_demand_type / "electricity", "cross-sector"/ @@ -495,92 +495,121 @@ ptc_out(i,v,t)$[tmodel_new(t)$ptc_value_scaled(i,v,t)] = ptc_value_scaled(i,v,t) * Case 2: the resource of one or more biomass classes ARE exhausted, i.e., BIOUSED.l(bioclass) = biosupply(bioclass) * Marginal Biomass Price = maximum difference between eq_bioused.m and eq_biousedlimit.m(bioclass) across all biomass classes in a region -repbioprice(r,t)$tmodel_new(t) = max{0, smax{bioclass$BIOUSED.l(bioclass,r,t), eq_bioused.m(r,t) - - sum{usda_region$r_usda(r,usda_region), eq_biousedlimit.m(bioclass,usda_region,t) } } } / pvf_onm(t) ; +repbioprice(r,t)$[tmodel_new(t)$tfuel(t)] = + [max{0, smax{bioclass$BIOUSED.l(bioclass,r,t), eq_bioused.m(r,t) - + sum{usda_region$r_usda(r,usda_region), eq_biousedlimit.m(bioclass,usda_region,t) } } } / pvf_onm(t) ] ; -* quantity of biomass used (convert from mmBTU to dry tons using biomass energy content) -bioused_out(bioclass,r,t)$tmodel_new(t) = BIOUSED.l(bioclass,r,t) / bio_energy_content ; -bioused_usda(bioclass,usda_region,t)$tmodel_new(t) = sum{r$r_usda(r,usda_region), bioused_out(bioclass,r,t) } ; +* quantity of biomass used in the power sector (convert from mmBTU to dry tons using biomass energy content) +bioused_out(bioclass,r,t)$[tmodel_new(t)$tfuel(t)] = + [BIOUSED.l(bioclass,r,t) / bio_energy_content ]; + +bioused_usda(bioclass,usda_region,t)$[tmodel_new(t)$tfuel(t)] = + [sum{r$r_usda(r,usda_region), bioused_out(bioclass,r,t) } ]; * 1e9 converts from MMBtu to Quads -repgasquant(cendiv,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 3)$tmodel_new(t)] = - sum{(gb,h), GASUSED.l(cendiv,gb,h,t) * hours(h) } * gas_scale/ 1e9 ; +repgasquant(cendiv,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 3)$tmodel_new(t)$tfuel(t)] = + [sum{(gb,h), GASUSED.l(cendiv,gb,h,t) * hours(h) } * gas_scale/ 1e9 ]; -repgasquant(cendiv,t)$[(Sw_GasCurve = 1 or Sw_GasCurve = 2)$tmodel_new(t)] = - ( sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], +repgasquant(cendiv,t)$[(Sw_GasCurve = 1 or Sw_GasCurve = 2 or Sw_FINITO_Link = 1)$tmodel_new(t)] = + [( sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t)} + sum{(v,r,h)$[valcap("dac_gas",v,r,t)$r_cendiv(r,cendiv)], hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE.l("DAC","dac_gas",v,r,h,t) }$Sw_DAC_Gas + sum{(p,i,v,r,h)$[r_cendiv(r,cendiv)$valcap(i,v,r,t)$smr(i)], hours(h) * smr_methane_rate * PRODUCE.l(p,i,v,r,h,t) }$Sw_H2 - ) / 1e9 ; + ) / 1e9 ]; -repgasquant_irt(i,r,t)$tmodel_new(t) = - ( sum{(v,h)$[valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], +repgasquant_irt(i,r,t)$[tmodel_new(t)] = + [( sum{(v,h)$[valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } + sum{(v,h)$[valcap("dac_gas",v,r,t)], hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE.l("DAC","dac_gas",v,r,h,t) }$Sw_DAC_Gas + sum{(p,v,h)$[valcap(i,v,r,t)$smr(i)], hours(h) * smr_methane_rate * PRODUCE.l(p,i,v,r,h,t) }$Sw_H2 - ) / 1e9 ; + ) / 1e9 ]; -repgasquant_nat(t)$tmodel_new(t) = sum{cendiv, repgasquant(cendiv,t) } ; +repgasquant_nat(t)$[tmodel_new(t)] = + [sum{cendiv, repgasquant(cendiv,t) } ]; *for reported gasprice (not that used to compute system costs) *scale back to $ / mmbtu -repgasprice(cendiv,t)$[(Sw_GasCurve = 0)$tmodel_new(t)$repgasquant(cendiv,t)] = - smax{gb$[sum{h, GASUSED.l(cendiv,gb,h,t) }], gasprice(cendiv,gb,t) } / gas_scale ; +repgasprice(cendiv,t)$[(Sw_GasCurve = 0)$tmodel_new(t)$repgasquant(cendiv,t)$tfuel(t)]= + [smax{gb$[sum{h, GASUSED.l(cendiv,gb,h,t) }], gasprice(cendiv,gb,t) } / gas_scale ]; -repgasprice(cendiv,t)$[(Sw_GasCurve = 2)$tmodel_new(t)$repgasquant(cendiv,t)] = - sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], +repgasprice(cendiv,t)$[(Sw_GasCurve = 2)$tmodel_new(t)$repgasquant(cendiv,t)$tfuel(t)] = + [sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], hours(h)*heat_rate(i,v,r,t)*fuel_price(i,r,t)*GEN.l(i,v,r,h,t) - } / (repgasquant(cendiv,t) * 1e9) ; + } / (repgasquant(cendiv,t) * 1e9) ]; -repgasprice_r(r,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 2)$tmodel_new(t)] = sum{cendiv$r_cendiv(r,cendiv), repgasprice(cendiv,t) } ; -repgasprice_r(r,t)$[(Sw_GasCurve = 1)$tmodel_new(t)] = - ( sum{(h,cendiv), - gasmultterm(cendiv,t) * szn_adj_gas(h) * cendiv_weights(r,cendiv) * - hours(h) } / sum{h, hours(h) } +* gas price by timeslice when linked with FINITO (see similar calculation in finito_report.gms) [$2004/MMBtu] +$ifthene.finitogasprice Sw_FINITO_Link == 1 +repgasprice_finito(cendiv,h,t)$[tmodel_new(t)$(not tfuel(t))] = + ( 1/obj_scale * 1/pvf_onm(t) * deflator('2018') * + [ max{ + smax{(eus,usep)$[map_roe_ec_usep("NG",eus,usep)], eq_supplydemand_fe_pool.m("NG",eus,usep,cendiv,h,t)}, + smax{(usep)$[map_roi_ec_usep("NG",usep)],eq_supplydemand_fe_roi.m("NG",usep,cendiv,h,t)}, + smax{(cfp2,cfvin2,fi,r)$[map_ec_ei("NG",fi)$valei_cf(cfp2,cfvin2,fi,r,t)$r_cendiv(r,cendiv)],eq_use_fi_cf.m(cfp2,cfvin2,fi,r,h,t)}, + smax{(tech,vin,fi,r)$[map_ec_ei("NG",fi)$valei_ind(tech,vin,fi,r,t)$r_cendiv(r,cendiv)],eq_use_fi_ind.m(tech,vin,fi,r,h,t)} + } / hours(h) + ] + ) +; +$else.finitogasprice + repgasprice_finito(cendiv,h,t)$[tmodel_new(t)$(not tfuel(t))] = 0 ; +$endif.finitogasprice + +repgasprice(cendiv,t)$[tmodel_new(t)$(not tfuel(t))] = +* average price over hours when linked with FINITO + sum{h, hours(h) * repgasprice_finito(cendiv,h,t) } +; + +repgasprice_r(r,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 2)$tmodel_new(t)$tfuel(t)] = + [sum{cendiv$r_cendiv(r,cendiv), repgasprice(cendiv,t) } ]; + +repgasprice_r(r,t)$[(Sw_GasCurve = 1)$tmodel_new(t)$tfuel(t)] = + sum{(h,cendiv), + gasmultterm(cendiv,t) * szn_adj_gas(h) * cendiv_weights(r,cendiv) * + hours(h) } / sum{h, hours(h) } - + smax((fuelbin,cendiv)$[VGASBINQ_REGIONAL.l(fuelbin,cendiv,t)$r_cendiv(r,cendiv)], gasbinp_regional(fuelbin,cendiv,t) ) + + smax((fuelbin,cendiv)$[VGASBINQ_REGIONAL.l(fuelbin,cendiv,t)$r_cendiv(r,cendiv)], gasbinp_regional(fuelbin,cendiv,t) ) - + smax(fuelbin$VGASBINQ_NATIONAL.l(fuelbin,t), gasbinp_national(fuelbin,t) ) - ) ; + + smax(fuelbin$VGASBINQ_NATIONAL.l(fuelbin,t), gasbinp_national(fuelbin,t) ) + ; -repgasprice(cendiv,t)$[(Sw_GasCurve = 1)$tmodel_new(t)$repgasquant(cendiv,t)] = - sum{(i,r)$r_cendiv(r,cendiv), repgasprice_r(r,t) * repgasquant_irt(i,r,t) } / repgasquant(cendiv,t) ; +repgasprice(cendiv,t)$[((Sw_GasCurve = 1) or (Sw_FINITO_Link = 1))$tmodel_new(t)$repgasquant(cendiv,t)] = + [sum{(i,r)$r_cendiv(r,cendiv), repgasprice_r(r,t) * repgasquant_irt(i,r,t) } / repgasquant(cendiv,t) ] ; repgasprice_nat(t)$[tmodel_new(t)$sum{cendiv, repgasquant(cendiv,t) }] = - sum{cendiv, repgasprice(cendiv,t) * repgasquant(cendiv,t) } - / sum{cendiv, repgasquant(cendiv,t) } ; + [sum{cendiv, repgasprice(cendiv,t) * repgasquant(cendiv,t) } + / sum{cendiv, repgasquant(cendiv,t) } ]; *======================================== * NATURAL GAS FUEL COSTS *======================================== gasshare_ba(r,cendiv,t)$[r_cendiv(r,cendiv)$tmodel_new(t)$repgasquant(cendiv,t)] = - sum{i$[valgen_irt(i,r,t)$gas(i)],repgasquant_irt(i,r,t) / repgasquant(cendiv,t) } ; + [sum{i$[valgen_irt(i,r,t)$gas(i)],repgasquant_irt(i,r,t) / repgasquant(cendiv,t) } ]; gasshare_techba(i,r,cendiv,t)$[r_cendiv(r,cendiv)$tmodel_new(t)$repgasquant(cendiv,t)$gas(i)] = - repgasquant_irt(i,r,t) / repgasquant(cendiv,t) ; + [repgasquant_irt(i,r,t) / repgasquant(cendiv,t) ]; -gasshare_cendiv(cendiv,t)$[sum{cendiv2,repgasquant(cendiv2,t)}] = repgasquant(cendiv,t) / sum{cendiv2,repgasquant(cendiv2,t)} ; +gasshare_cendiv(cendiv,t)$[sum{cendiv2,repgasquant(cendiv2,t)}] = + [repgasquant(cendiv,t) / sum{cendiv2,repgasquant(cendiv2,t)} ]; -gascost_cendiv(cendiv,t)$tmodel_new(t) = +* cost of natural gas - standalone ReEDS +gascost_cendiv(cendiv,t)$[tmodel_new(t)$tfuel(t)] = *cost of natural gas for Sw_GasCurve = 2 (static natural gas prices) - + sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t) - $[not bio(i)]$[not cofire(i)]$[Sw_GasCurve = 2]], - hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN.l(i,v,r,h,t) } + + [ sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t) + $[not bio(i)]$[not cofire(i)]$[Sw_GasCurve = 2]], + hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN.l(i,v,r,h,t) } *cost of natural gas for Sw_GasCurve = 0 (census division supply curves natural gas prices) - + sum{gb, sum{h,hours(h) * GASUSED.l(cendiv,gb,h,t) } * gasprice(cendiv,gb,t) - }$[Sw_GasCurve = 0] + + sum{gb, sum{h,hours(h) * GASUSED.l(cendiv,gb,h,t) } * gasprice(cendiv,gb,t) }$[Sw_GasCurve = 0] *cost of natural gas for Sw_GasCurve = 3 (national supply curve for natural gas prices with census division multipliers) - + sum{(h,gb), hours(h) * GASUSED.l(cendiv,gb,h,t) - * gasadder_cd(cendiv,t,h) + gasprice_nat_bin(gb,t) - }$[Sw_GasCurve = 3] + + sum{(h,gb), hours(h) * GASUSED.l(cendiv,gb,h,t) + * gasadder_cd(cendiv,t,h) + gasprice_nat_bin(gb,t) }$[Sw_GasCurve = 3] *cost of natural gas for Sw_GasCurve = 1 (national and census division supply curves for natural gas prices) *first - anticipated costs of gas consumption given last year's amount + (sum{(i,v,r,h)$[valgen(i,v,r,t)$gas(i)], @@ -593,15 +622,25 @@ gascost_cendiv(cendiv,t)$tmodel_new(t) = + sum{(fuelbin), gasbinp_national(fuelbin,t) * VGASBINQ_NATIONAL.l(fuelbin,t) } * gasshare_cendiv(cendiv,t) - )$[Sw_GasCurve = 1]; + )$[Sw_GasCurve = 1] + ]; + +* cost of natural gas - linked with FINITO ('not tfuel' indicates years using FINITO supply curves) +gascost_cendiv(cendiv,t)$[tmodel_new(t)$(not tfuel(t))] = +* cost = gas price multiplied by gas usage [$ = $/MMBtu * MMBtu/MWh * MW * h] + sum{(i,v,r,h)$[valgen(i,v,r,t)$gas(i)$r_cendiv(r, cendiv)], + repgasprice_finito(cendiv,h,t) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) * hours(h) + } +; + *======================================== * BIOFUEL COSTS *======================================== -bioshare_techba(i,r,t)$[(cofire(i) or bio(i))$tmodel_new(t)] = +bioshare_techba(i,r,t)$[(cofire(i) or bio(i))$tmodel_new(t)$tfuel(t)] = * biofuel-based generation of tech i in the BA (biopower + cofire) - (( sum{(v,h)$[valgen(i,v,r,t)$bio(i)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } + [ (( sum{(v,h)$[valgen(i,v,r,t)$bio(i)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } + sum{(v,h)$[cofire(i)$valgen(i,v,r,t)], bio_cofire_perc * hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } ) / * biofuel-based generation of all techs in the BA (biopower + cofire) @@ -611,6 +650,7 @@ bioshare_techba(i,r,t)$[(cofire(i) or bio(i))$tmodel_new(t)] = )$[ sum{(ii,v,h)$[valgen(ii,v,r,t)$bio(ii)], hours(h) * heat_rate(ii,v,r,t) * GEN.l(ii,v,r,h,t) } + sum{(ii,v,h)$[cofire(ii)$valgen(ii,v,r,t)], bio_cofire_perc * hours(h) * heat_rate(ii,v,r,t) * GEN.l(ii,v,r,h,t) } ] + ] ; *========================= @@ -1740,6 +1780,11 @@ error_check('z') = ( ) * account for penalty paid to deploy capacity beyond interconnection queue limits + sum{(tg,r), cap_penalty(tg) * CAP_ABOVE_LIM.l(tg,r,t) } +* account for costs from FINITO: deflate from $2018 to $2004, +* remove any FINITO scaling, and then apply ReEDS scaling +$ifthene.linked_objective Sw_FINITO_Link==1 + + cost_scale * ( Z_finito.l(t)$t_finito(t) * deflator('2018') / obj_scale ) +$endif.linked_objective } ) / z.l ; @@ -2033,6 +2078,8 @@ load_cat("h2_network",r,t)$tmodel_new(t) = ; load_cat("dac",r,t)$tmodel_new(t) = sum{i$dac(i), prod_load_ann(i,r,t) } ; +load_cat("finito",r,t)$tmodel_new(t) = sum{h, hours(h) * USE_ELE_FINITO.l(r,h,t) } ; + *======================================== * H2 NETWORK *======================================== diff --git a/e_report_dump.py b/e_report_dump.py index faa3a20a..aa3aa78a 100644 --- a/e_report_dump.py +++ b/e_report_dump.py @@ -67,6 +67,8 @@ def dfdict_to_h5( """ Write dictionary of dataframes to one .h5 file """ + print(f"Saving results to {filepath}") + ### unless a subset is specified, iterate over all keys in dict _symbol_list = dfdict.keys() if symbol_list is None else symbol_list @@ -108,6 +110,8 @@ def dfdict_to_excel( """ Write dictionary of dataframes to one .xlsx file """ + print(f"Saving results to {filepath}") + ### unless a subset is specified, iterate over all keys in dict _symbol_list = dfdict.keys() if symbol_list is None else symbol_list @@ -201,10 +205,13 @@ def postprocess_outputs(case, outputs_path=None, verbose=0): _outputs_path = os.path.join(case, 'outputs') if outputs_path is None else outputs_path ## System cost - reeds.output_calc.calc_systemcost(case).to_csv( - os.path.join(_outputs_path, 'post_systemcost_annualized.csv'), - index=False, - ) + try: + reeds.output_calc.calc_systemcost(case).to_csv( + os.path.join(_outputs_path, 'post_systemcost_annualized.csv'), + index=False, + ) + except FileNotFoundError: + print("Skipping system cost calculation, missing files") ## Reinforcement and spur-line reeds.output_calc.calc_reinforcement_spur_capacity_miles(case).to_csv( @@ -276,6 +283,18 @@ def postprocess_outputs(case, outputs_path=None, verbose=0): ) print("Finished loading outputs gdx") + ### FINITO outputs when running linked model + if int(sw.GSw_FINITO_Link): + print("Loading FINITO outputs gdx") + # TODO: this is much slower than loading the ReEDS outputs, + # need to figure out why and address + finito_outputs = gdxpds.to_dataframes( + os.path.join(outputs_path, f"finito_reeds_outputs_{os.path.basename(case)}.gdx") + ) + print("Finished loading FINITO outputs gdx") + dict_out.update(finito_outputs) + + write_dfdict( dfdict=dict_out, outputs_path=outputs_path, diff --git a/e_report_params.csv b/e_report_params.csv index 70b37020..f7792cb3 100644 --- a/e_report_params.csv +++ b/e_report_params.csv @@ -142,6 +142,7 @@ raw_op_cost(t),2004$,sum of operational costs from systemcost,,, RE_gen_price_nat(t),2004$/MWh,marginal cost of the national RE generation constraint,,, "repbioprice(r,t)",2004$/MMBtu,highest marginal bioprice of utilized bins for each region,,, "repgasprice(cendiv,t)",2004$/MMBtu,highest marginal gas price of utilized gas bins for each census division,,, +"repgasprice_finito(cendiv,allh,t)",2004$/MMBtu,highest marginal gas price of utilized gas bins by hour for each census division when using linked model,,, "repgasprice_r(r,t)",2004$/MMBtu,highest marginal gas price of utilized gas bins for each region,1,, repgasprice_nat(t),2004$/MMBtu,weighted-average national natural gas price assuming that plants pay the marginal price,,, "repgasquant(cendiv,t)",Quads,quantity of gas consumed in each census division,,, diff --git a/input_processing/copy_files.py b/input_processing/copy_files.py index 90e22a09..8c53b164 100644 --- a/input_processing/copy_files.py +++ b/input_processing/copy_files.py @@ -1614,6 +1614,11 @@ def main(reeds_path, inputs_case): source_deflator_map = get_source_deflator_map(reeds_path) + ### Rewrite the finito scalar as GAMS-readable definitions + if int(sw['GSw_FINITO_Link']) == 1: + pd.read_csv(os.path.join(os.path.dirname(inputs_case),'finito','inputs','scalars.csv')).iloc[:, :-1].to_csv(os.path.join(os.path.dirname(inputs_case),'finito','inputs','scalars.csv'),index=False) + scalar_csv_to_txt(os.path.join(os.path.dirname(inputs_case),'finito','inputs','scalars.csv')) + # Copy non-region files write_non_region_files(non_region_files, sw, inputs_case, regions_and_agglevel, source_deflator_map) diff --git a/input_processing/hourly_load.py b/input_processing/hourly_load.py index 8b2114fb..a29ee8f5 100644 --- a/input_processing/hourly_load.py +++ b/input_processing/hourly_load.py @@ -596,6 +596,44 @@ def reaggregate_to_model_regions( return regional_load_hourly +#TODO: add docstrings +def get_hourly_finito_load( + inputs_case: str, +) -> pd.DataFrame: + # load reference FINITO load + inputs_case_finito = Path(inputs_case).parent / 'finito' / Path(inputs_case).name + load_finito = pd.read_csv(inputs_case_finito / "load_finito.csv") + + # reshape to match load data. with year in index and region in columns + load_finito = load_finito.melt(id_vars='r', var_name='year', value_name='load_MWh') + load_finito.year = load_finito.year.astype(int) + load_finito = load_finito.pivot(index='year', columns='r', values='load_MWh') + load_finito.columns.name = None + + # allocate annual load to hours, assuming flat demand + # TODO: should this use h_weight_finito? + hours_per_year = 8760 + load_hourly_finito = load_finito / hours_per_year + load_hourly_finito = load_hourly_finito.astype(np.float32) + + return load_hourly_finito + +def remove_finito_load( + load_hourly: pd.DataFrame, + inputs_case: str, +) -> pd.DataFrame: + + # get FINITO reference load + load_hourly_finito = get_hourly_finito_load(inputs_case) + + # subtract FINITO reference load from ReEDS load data, + # aligning by model year (index) and region (columns) + result = load_hourly - load_hourly_finito + + #TODO: add this + # Validation check: confirm that total load - FINITO load > 0 + + return load_hourly #%% =========================================================================== ### --- MAIN FUNCTION --- @@ -697,6 +735,15 @@ def main(reeds_path, inputs_case): peakload = calculate_peak_load(regional_load_hourly, hierarchy) #%%%######################################### + # -- FINITO Load Adjustment -- # + ############################################# + + # note that this step occurs after peakload calculation so that the latter + # includes a baseline estimate of industrial load captured by FINITO + if int(sw.GSw_FINITO_Link): + regional_load_hourly = remove_finito_load(regional_load_hourly, inputs_case) + + ############################################# # -- DR Shed Load Modifications -- # ############################################# diff --git a/input_processing/hourly_repperiods.py b/input_processing/hourly_repperiods.py index 2c57bcc8..7d13b20e 100644 --- a/input_processing/hourly_repperiods.py +++ b/input_processing/hourly_repperiods.py @@ -36,6 +36,7 @@ import hourly_plots sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) import reeds +import hourly_load ## Time the operation of this script tic = datetime.datetime.now() @@ -75,11 +76,18 @@ def szn2period(szn): # -- Load Processing -- # ############################### -def get_load(inputs_case, keep_modelyear=None, keep_weatheryears=[2012]): +def get_load(inputs_case, sw, keep_modelyear=None, keep_weatheryears=[2012]): """ """ ### Subset to modeled regions load = reeds.io.read_file(os.path.join(inputs_case,'load.h5'), parse_timestamps=True) + + ### When running the linked model (GSw_FINITO_Link=1) we add reference load estimates + ### for FINITO load back in to use when selecting representative periods + if int(sw.GSw_FINITO_Link): + load_hourly_finito = hourly_load.get_hourly_finito_load(inputs_case) + load = load + load_hourly_finito + ### Subset to keep_modelyear if provided if keep_modelyear: load = load.loc[keep_modelyear].copy() @@ -578,6 +586,7 @@ def main( print("Collecting 8760 load data") load = get_load( inputs_case=inputs_case, + sw=sw, keep_modelyear=(int(sw['GSw_HourlyClusterYear']) if int(sw['GSw_HourlyClusterYear']) in modelyears else max(modelyears)), @@ -736,7 +745,7 @@ def main( and (sw['GSw_PRM_StressSeedLoadLevel'].lower() not in ['false','none']) ): ## Get load for all model and weather years - load_allyears = get_load(inputs_case, keep_weatheryears='all').loc[modelyears] + load_allyears = get_load(inputs_case, sw, keep_weatheryears='all').loc[modelyears] ## Add descriptive index load_allyears = load_allyears.merge( timestamps[['year', 'yperiod', 'h_of_period']], left_on='datetime', right_index=True) @@ -925,10 +934,14 @@ def main( description='Create the necessary 8760 and capacity factor data for hourly resolution') parser.add_argument('reeds_path', help='ReEDS-2.0 directory') parser.add_argument('inputs_case', help='ReEDS-2.0/runs/{case}/inputs_case directory') + parser.add_argument('--nolog', '-n', default=False, action='store_true', help='turn off logging for debugging') + args = parser.parse_args() reeds_path = args.reeds_path inputs_case = args.inputs_case + logging = not args.nolog + # #%% Settings for testing # reeds_path = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) @@ -938,10 +951,11 @@ def main( # interactive = True #%% Set up logger - log = reeds.log.makelog( - scriptname=__file__, - logpath=os.path.join(inputs_case,'..','gamslog.txt'), - ) + if logging: + log = reeds.log.makelog( + scriptname=__file__, + logpath=os.path.join(inputs_case,'..','gamslog.txt'), + ) print('Starting hourly_repperiods.py') #%% Inputs from switches sw = reeds.io.get_switches(inputs_case) @@ -956,6 +970,7 @@ def main( sw=sw, reeds_path=reeds_path, inputs_case=inputs_case, periodtype='rep', make_plots=1, + logging=logging ) ############################################ diff --git a/input_processing/hourly_writetimeseries.py b/input_processing/hourly_writetimeseries.py index 5fa32608..57c5f436 100644 --- a/input_processing/hourly_writetimeseries.py +++ b/input_processing/hourly_writetimeseries.py @@ -199,6 +199,9 @@ def get_yearly_demand(sw, hmap_myr, hmap_allyrs, inputs_case, periodtype='rep'): reload the raw demand and extract the demand on the modeled days for each year. """ ### Get original demand data, subset to cluster year + ### Note that this does not include FINITO demand when running with the linked model + ### (GSw_FINITO_Link=1) even though the reference estimates for that load is used when + ### identifying the representative days load_in = reeds.io.read_file( os.path.join(inputs_case,'load.h5'), parse_timestamps=True).unstack(level=0) load_in.columns = load_in.columns.rename(['r','t']) diff --git a/inputs/emission_constraints/co2_tax.csv b/inputs/emission_constraints/co2_tax.csv index d9e79551..1e08f624 100644 --- a/inputs/emission_constraints/co2_tax.csv +++ b/inputs/emission_constraints/co2_tax.csv @@ -1,42 +1,42 @@ -t,default,t050_2035,t050_2050,t100_2035,t100_2050,t150_2035,t150_2050,t200_2035,t200_2050 -2010,0,0,0,0,0,0,0,0,0 -2011,0,0,0,0,0,0,0,0,0 -2012,0,0,0,0,0,0,0,0,0 -2013,0,0,0,0,0,0,0,0,0 -2014,0,0,0,0,0,0,0,0,0 -2015,0,0,0,0,0,0,0,0,0 -2016,0,0,0,0,0,0,0,0,0 -2017,0,0,0,0,0,0,0,0,0 -2018,0,0,0,0,0,0,0,0,0 -2019,0,0,0,0,0,0,0,0,0 -2020,0,0,0,0,0,0,0,0,0 -2021,0,0,0,0,0,0,0,0,0 -2022,30,1,1,1,1,1,1,1,1 -2023,31.5,4.77,2.75,8.62,4.54,12.46,6.32,16.31,8.11 -2024,33.075,8.54,4.5,16.23,8.07,23.92,11.64,31.62,15.21 -2025,34.72875,12.31,6.25,23.85,11.61,35.38,16.96,46.92,22.32 -2026,36.4651875,16.08,8,31.46,15.14,46.85,22.29,62.23,29.43 -2027,38.28844688,19.85,9.75,39.08,18.68,58.31,27.61,77.54,36.54 -2028,40.20286922,23.62,11.5,46.69,22.21,69.77,32.93,92.85,43.64 -2029,42.21301268,27.38,13.25,54.31,25.75,81.23,38.25,108.15,50.75 -2030,44.32366331,31.15,15,61.92,29.29,92.69,43.57,123.46,57.86 -2031,46.53984648,34.92,16.75,69.54,32.82,104.15,48.89,138.77,64.96 -2032,48.8668388,38.69,18.5,77.15,36.36,115.62,54.21,154.08,72.07 -2033,51.31018074,42.46,20.25,84.77,39.89,127.08,59.54,169.38,79.18 -2034,53.87568978,46.23,22,92.38,43.43,138.54,64.86,184.69,86.29 -2035,56.56947427,50,23.75,100,46.96,150,70.18,200,93.39 -2036,59.39794798,50,25.5,100,50.5,150,75.5,200,100.5 -2037,62.36784538,50,27.25,100,54.04,150,80.82,200,107.61 -2038,65.48623765,50,29,100,57.57,150,86.14,200,114.71 -2039,68.76054953,50,30.75,100,61.11,150,91.46,200,121.82 -2040,72.19857701,50,32.5,100,64.64,150,96.79,200,128.93 -2041,75.80850586,50,34.25,100,68.18,150,102.11,200,136.04 -2042,79.59893115,50,36,100,71.71,150,107.43,200,143.14 -2043,83.57887771,50,37.75,100,75.25,150,112.75,200,150.25 -2044,87.7578216,50,39.5,100,78.79,150,118.07,200,157.36 -2045,92.14571268,50,41.25,100,82.32,150,123.39,200,164.46 -2046,96.75299831,50,43,100,85.86,150,128.71,200,171.57 -2047,101.5906482,50,44.75,100,89.39,150,134.04,200,178.68 -2048,106.6701806,50,46.5,100,92.93,150,139.36,200,185.79 -2049,112.0036897,50,48.25,100,96.46,150,144.68,200,192.89 -2050,117.6038742,50,50,100,100,150,150,200,200 \ No newline at end of file +t,default,t050_2035,t050_2050,t100_2035,t100_2050,t150_2035,t150_2050,t200_2035,t200_2050,t200alt_2050 +2010,0,0,0,0,0,0,0,0,0,0 +2011,0,0,0,0,0,0,0,0,0,0 +2012,0,0,0,0,0,0,0,0,0,0 +2013,0,0,0,0,0,0,0,0,0,0 +2014,0,0,0,0,0,0,0,0,0,0 +2015,0,0,0,0,0,0,0,0,0,0 +2016,0,0,0,0,0,0,0,0,0,0 +2017,0,0,0,0,0,0,0,0,0,0 +2018,0,0,0,0,0,0,0,0,0,0 +2019,0,0,0,0,0,0,0,0,0,0 +2020,0,0,0,0,0,0,0,0,0,0 +2021,0,0,0,0,0,0,0,0,0,0 +2022,30,1,1,1,1,1,1,1,1,0 +2023,31.5,4.77,2.75,8.62,4.54,12.46,6.32,16.31,8.11,0 +2024,33.075,8.54,4.5,16.23,8.07,23.92,11.64,31.62,15.21,0 +2025,34.72875,12.31,6.25,23.85,11.61,35.38,16.96,46.92,22.32,0 +2026,36.4651875,16.08,8,31.46,15.14,46.85,22.29,62.23,29.43,0 +2027,38.28844688,19.85,9.75,39.08,18.68,58.31,27.61,77.54,36.54,0 +2028,40.20286922,23.62,11.5,46.69,22.21,69.77,32.93,92.85,43.64,0 +2029,42.21301268,27.38,13.25,54.31,25.75,81.23,38.25,108.15,50.75,0 +2030,44.32366331,31.15,15,61.92,29.29,92.69,43.57,123.46,57.86,0 +2031,46.53984648,34.92,16.75,69.54,32.82,104.15,48.89,138.77,64.96,10 +2032,48.8668388,38.69,18.5,77.15,36.36,115.62,54.21,154.08,72.07,20 +2033,51.31018074,42.46,20.25,84.77,39.89,127.08,59.54,169.38,79.18,30 +2034,53.87568978,46.23,22,92.38,43.43,138.54,64.86,184.69,86.29,40 +2035,56.56947427,50,23.75,100,46.96,150,70.18,200,93.39,50 +2036,59.39794798,50,25.5,100,50.5,150,75.5,200,100.5,60 +2037,62.36784538,50,27.25,100,54.04,150,80.82,200,107.61,70 +2038,65.48623765,50,29,100,57.57,150,86.14,200,114.71,80 +2039,68.76054953,50,30.75,100,61.11,150,91.46,200,121.82,90 +2040,72.19857701,50,32.5,100,64.64,150,96.79,200,128.93,100 +2041,75.80850586,50,34.25,100,68.18,150,102.11,200,136.04,110 +2042,79.59893115,50,36,100,71.71,150,107.43,200,143.14,120 +2043,83.57887771,50,37.75,100,75.25,150,112.75,200,150.25,130 +2044,87.7578216,50,39.5,100,78.79,150,118.07,200,157.36,140 +2045,92.14571268,50,41.25,100,82.32,150,123.39,200,164.46,150 +2046,96.75299831,50,43,100,85.86,150,128.71,200,171.57,160 +2047,101.5906482,50,44.75,100,89.39,150,134.04,200,178.68,170 +2048,106.6701806,50,46.5,100,92.93,150,139.36,200,185.79,180 +2049,112.0036897,50,48.25,100,96.46,150,144.68,200,192.89,190 +2050,117.6038742,50,50,100,100,150,150,200,200,200 diff --git a/reeds/reedsplots.py b/reeds/reedsplots.py index da1013b1..1630ca17 100644 --- a/reeds/reedsplots.py +++ b/reeds/reedsplots.py @@ -5020,6 +5020,7 @@ def plot_seed_stressperiods( ### Recalculate peak load days since we dropped duplicates above load_allyears = hourly_repperiods.get_load( os.path.join(case, 'inputs_case'), + sw, keep_weatheryears='all').loc[years] timestamps = pd.read_csv(os.path.join(case, 'inputs_case', 'rep', 'timestamps.csv')) resource_adequacy_years = [int(y) for y in sw.resource_adequacy_years.split('_')] diff --git a/restart_runs.py b/restart_runs.py index dee9a85c..71eb72e7 100644 --- a/restart_runs.py +++ b/restart_runs.py @@ -4,6 +4,8 @@ import subprocess import argparse import pandas as pd +import re +import reeds from glob import glob from runbatch import submit_slurm_parallel_jobs from runstatus import get_run_status @@ -92,7 +94,14 @@ #%% Copy additional files if desired for f in more_copyfiles: - shutil.copy(os.path.join(reeds_path,f), os.path.join(case,f)) + if f.lower().startswith('finito'): + # if file starts with 'finito' append the finito directory + sw = reeds.io.get_switches(case) + if int(sw.GSw_FINITO_Link): + f_copy = re.sub("^finito/", "" , f, flags=re.IGNORECASE) + shutil.copy(os.path.join(sw.finito_dir,f_copy), os.path.join(case,f)) + else: + shutil.copy(os.path.join(reeds_path,f), os.path.join(case,f)) #%% Make a backup copy of the original bash and sbatch scripts callfile = os.path.join(case,f'call_{casename}.sh') @@ -165,7 +174,8 @@ # It is a single case or we are not on HPC if copy_srun_template: writelines_srun_case = writelines_srun.copy() - writelines_srun_case.append(f"\n#SBATCH --job-name={casename}\n") + writelines_srun_case.append(f"#SBATCH --job-name={casename}") + writelines_srun_case.append(f"#SBATCH --output={os.path.join(case, 'slurm-%j.out')}\n") writelines_srun_case.append(f"sh {callfile}") with open(sbatchfile, 'w') as f: for line in writelines_srun_case: diff --git a/runbatch.py b/runbatch.py index 33b32979..106d2735 100644 --- a/runbatch.py +++ b/runbatch.py @@ -2,6 +2,10 @@ ### --- IMPORTS --- ### =========================================================================== +#%% =========================================================================== +### --- IMPORTS --- +### =========================================================================== + import os import git import queue @@ -11,6 +15,8 @@ import csv import importlib import numpy as np +import gdxpds +from gdxpds.gdx import GdxFile,GdxSymbol,GamsDataType import pandas as pd import subprocess import re @@ -18,6 +24,7 @@ import argparse from pathlib import Path import reeds +import sys # Assert core programs are accessible CORE_PROGRAMS = ["gams"] @@ -99,12 +106,23 @@ def create_case_lists(df_cases:pd.DataFrame, BatchName:str, single:str=''): continue # Add switch settings to list of options passed to GAMS shcom = f' --case={BatchName}_{case}' - for i,v in df_cases[case].items(): + # (linked) Combine the cases files for the linked model + if int(df_cases[case].loc['GSw_FINITO_Link']) == 1 : + # add the FINITO switches to the caseSwitches + df_cases_combine=linked_cases(df_cases,case) + for i,v in df_cases_combine.items(): + #exclude certain switches that don't need to be passed to GAMS + if i not in ['file_replacements','keep_run_terminal']: + shcom = shcom + ' --' + i + '=' + v + caseList.append(shcom) + caseSwitches.append(df_cases_combine.to_dict()) + else: + for i,v in df_cases[case].items(): #exclude certain switches that don't need to be passed to GAMS - if i not in ['file_replacements','keep_run_terminal']: - shcom += f' --{i}={v}' - caseList.append(shcom) - caseSwitches.append(df_cases[case].to_dict()) + if i not in ['file_replacements','keep_run_terminal']: + shcom += f' --{i}={v}' + caseList.append(shcom) + caseSwitches.append(df_cases[case].to_dict()) return caseSwitches, casenames, caseList @@ -554,6 +572,7 @@ def solvestring_sequential( 'GSw_StateCO2ImportLevel', 'GSw_StartMarkets', 'GSw_ValStr', + 'GSw_FINITO_Link', 'solver', 'debug', 'startyear', @@ -807,6 +826,71 @@ def setup_window( OPATH.writelines( "python valuestreams.py" + '\n') +# (linked) update 'df_cases' to include FINITO switches +def linked_cases(df_cases,case): + """ + Updates the cases dataframe to include FINITO-specific switches. + When a switch is duplicated in FINITO and ReEDS, then we default to + the ReEDS value. + + For the FINITO switches, the combined cases file defaults to the case-specific + 'Default Value' in cases_linked.csv first, before using the universal FINITO + 'Default Value' in cases.csv for any un-assigned switches. + """ + # check for valid finito_dir + if not os.path.isdir(df_cases[case]['finito_dir']): + raise ValueError(f"finito_dir = {df_cases[case]['finito_dir']} is not a valid path.") + + # define path to and read the FINITO check_inputs function + finito_check_inputs_path = os.path.join(df_cases[case]['finito_dir'], 'input_processing', 'processing') + sys.path.append(finito_check_inputs_path) + from check_inputs import check_inputs + + ### load the default values for all FINITO switches from ~\FINITO\cases.csv + df_cases_finito = pd.read_csv(os.path.join(df_cases[case]['finito_dir'],'cases.csv'), dtype=object, index_col=0) + df_cases_finito = df_cases_finito[['Choices', 'Default Value']] + + ### load the scenario-specific switches from ~\FINITO\cases_linked.csv + cases_linked_path = os.path.join(df_cases[case]['finito_dir'],f"cases_{df_cases[case]['finito_cases_file']}.csv") + df_cases_suf_finito = pd.read_csv(cases_linked_path, dtype=object, index_col=0) + ## check that case names are unique in cases_linked.csv + # grab the scenario names **exactly** as they appear in the csv file + with open(cases_linked_path, 'r', newline='') as csvfile: + reader = csv.reader(csvfile) + header = next(reader) + # find the duplicate column names and raise an error if any are found + duplicate_columns = set([x for x in header if header.count(x) > 1]) + if duplicate_columns: + raise ValueError(f"The FINITO cases_{df_cases[case]['finito_cases_file']}.csv has the following duplicate column names: {duplicate_columns}") + ### identify the FINITO case + if df_cases[case]['finito_case'] == 'same': + finito_case=case + else: + finito_case=df_cases[case]['finito_case'] + # ensures **exact** match of names between the ReEDS cases_{}.csv and the FINITO cases_linked.csv + if finito_case not in (df_cases_suf_finito.columns): + raise ValueError(f"The 'finito_case' input '{finito_case}' in the ReEDS cases file does not exist in FINITO's cases_{df_cases[case]['finito_cases_file']}.csv.") + + ### first use 'Default Value' from the FINITO cases_linked.csv to fill missing switches + if 'Default Value' in df_cases_suf_finito.columns: + df_cases_suf_finito[finito_case] = df_cases_suf_finito[finito_case].fillna(df_cases_suf_finito['Default Value']) + ### then, use 'Default Value' from the FINITO cases.csv to fill un-assigned switches + df_cases_suf_finito.drop(['Choices','Default Value'], axis='columns',inplace=True, errors='ignore') + df_cases_finito = df_cases_finito.join(df_cases_suf_finito, how='outer') + df_cases_finito[finito_case] = df_cases_finito[finito_case].fillna(df_cases_finito['Default Value']) + + #### create new dataframe for the combined ReEDS and FINITO switches + df_cases_combine = pd.concat([df_cases[case],df_cases_finito[finito_case]]) + ### drop duplicated switches, defaulting to reeds + df_cases_combine = df_cases_combine[~df_cases_combine.index.duplicated(keep='first')] + + #%% Check for incompatibility of FINITO switches + model_sectors = df_cases_finito['Default Value']['focus_sectors'].split('.') + check_inputs(case = case, df_case = df_cases_combine, model_sectors=model_sectors) + + return df_cases_combine + + #%% =========================================================================== ### --- PROCEDURE --- ### =========================================================================== @@ -1260,6 +1344,7 @@ def write_batch_script( yearset_augur = os.path.join('inputs_case','modeledyears.csv') toLogGamsString = ' logOption=4 logFile=gamslog.txt appendLog=1 ' + ## Copy code folders for dirname in ['reeds', 'ReEDS_Augur', 'input_processing', 'reeds2pras']: shutil.copytree( @@ -1271,6 +1356,66 @@ def write_batch_script( #make the augur_data folder os.makedirs(os.path.join(casedir,'ReEDS_Augur','augur_data'), exist_ok=True) os.makedirs(os.path.join(casedir,'ReEDS_Augur','PRAS'), exist_ok=True) + + ## (linked) Copy FINITO code folders and inputs into [casedir]/finito + if int(caseSwitches['GSw_FINITO_Link'])==1: + ## Redirecting the FINITO case-specific directories + # ... finito directory within the case directory + casedir_finito = os.path.join(casedir,'finito') + # ... define the inputs case directory for FINITO + inputs_case_finito = os.path.join(casedir,'finito','inputs_case') + + ## copy directories + os.makedirs(casedir_finito, exist_ok=True) + shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'inputs'),os.path.join(casedir,'finito', 'inputs')) + shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'input_processing'),os.path.join(casedir,'finito', 'input_processing')) + shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'model'),os.path.join(casedir,'finito', 'model')) + shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'visualization'),os.path.join(casedir, 'finito', 'visualization')) + shutil.copytree(os.path.join(caseSwitches['finito_dir'],'inputs'), os.path.join(casedir,'inputs'), dirs_exist_ok=True) + ## copy over the FINITO cases files + shutil.copy2(os.path.join(caseSwitches['finito_dir'], 'cases.csv'), os.path.join(casedir, 'finito')) + shutil.copy2(os.path.join(caseSwitches['finito_dir'], f"cases_{caseSwitches['finito_cases_file']}.csv"), os.path.join(casedir, 'finito')) + + + #%% (GSw_Trade_PriceResponse > 0) If doing a price-responsive trade run, retrieve the reference exports/imports prices + if int(caseSwitches['GSw_Trade_PriceResponse']) > 0: + initialize_price_response_path = os.path.join(casedir_finito, 'input_processing', 'processing', 'initialize_price_response.py') + # Collect all arguments for initialize_price_response.py + initialize_price_response_args = f" -c {casedir} -b {BatchName} -cr {caseSwitches['GSw_Trade_PriceResponse_RefScen']} -l {caseSwitches['GSw_FINITO_Link']}" + # Call initialize_price_response.py file before starting the runs + os.system('python ' + initialize_price_response_path + initialize_price_response_args) + + #%% Filter and copy all input files for each scenario + # Call FINITO copy_files.py file before starting the runs + copy_files_run = os.system( + 'python ' + os.path.join(caseSwitches['finito_dir'], 'input_processing', 'processing', 'copy_files.py') + + f" -c {casedir_finito} -d {inputs_case_finito} --link" + ) + + # Print an error message if copy_files.py encounters any issue + if copy_files_run != 0: + print("\n❌ ERROR: copy_files.py encountered an issue and did not complete successfully.") + print("Please check the console output above for details.") + print("The issue could be due to regionality, focus sector filtering, or file reading errors.") + os._exit(1) + + ## Populate sets for each linked run + + # Pass AEO base year as an argument + # NOTE: make this dynamic before merging in? + aeo_year = caseSwitches['aeo_year'] + # Collect all arguments for autopop_set.py + focus_sectors = caseSwitches['focus_sectors'].replace('.', ' ') + autopop_args = f" -c {casedir_finito} -d {inputs_case_finito} -a {aeo_year} -s {focus_sectors} -f {caseSwitches['GSw_FixedCostSupply']} -rwf {caseSwitches['GSw_Trade_Partners']} -e {caseSwitches['GSw_ROE_EndUses']}" + # Call autopop_set.py file before starting the runs + autopop_path = os.path.join(caseSwitches['finito_dir'], 'input_processing','processing','autopop_set.py') + os.system('python ' + autopop_path + autopop_args) + ## Call read_mecs_heat.py to generate heat/nonheat/feedstock ratios for FINITO Rest of Industry (ROI) + mecs_sectors = focus_sectors + read_mecs_path = os.path.join(caseSwitches['finito_dir'], 'input_processing', 'processing', 'mecs', 'read_mecs_heat.py') + # Collect all arguments for read_mecs_heat.py + read_mecs_args = f' -s {mecs_sectors} -d {inputs_case_finito}' + os.system('python ' + read_mecs_path + read_mecs_args) ###### Replace files according to 'file_replacements' in cases. Ignore quotes in input text. # << is used to separate the file that is to be replaced from the file that is used @@ -1309,6 +1454,7 @@ def write_batch_script( OPATH.writelines("module load conda \n") OPATH.writelines("module load gams \n") + OPATH.writelines("conda deactivate \n") OPATH.writelines("conda activate reeds2 \n") OPATH.writelines('export R_LIBS_USER="$HOME/rlib" \n\n\n') @@ -1412,7 +1558,8 @@ def write_batch_script( + ' gdxcompress=1' + toLogGamsString + f"--fname={batch_case}" - + f" --GSw_calc_powfrac={caseSwitches['GSw_calc_powfrac']} \n" + + f" --GSw_calc_powfrac={caseSwitches['GSw_calc_powfrac']}" + + f" --first_year_finito={caseSwitches['first_year_finito']} \n" ) OPATH.writelines(writescripterrorcheck("e_report.gms")) if not LINUXORMAC: @@ -1425,6 +1572,32 @@ def write_batch_script( + f" {os.path.join(reeds_path,'postprocessing','diagnose','diagnose_process.py')}" + f" --casepath {casedir} \n\n" ) + OPATH.writelines(f'python e_report_dump.py {casedir} -c\n\n') + # (linked) calls finito reporting (from [casedir]/finito/visualization) + if int(caseSwitches['GSw_FINITO_Link'])==1: + OPATH.writelines( + 'gams ' + + f"{os.path.join(casedir,'finito', 'model', 'finito_report.gms')}" + + f" o={os.path.join('lstfiles',f'finito_report_{batch_case}.lst')}" + + (' license=gamslice.txt' if hpc else '') + + (' r=$r' if LINUXORMAC else ' r=!r!') + + ' gdxcompress=1' + + toLogGamsString + + f"--case={batch_case}" + + f" --casedir={casedir}" + + f" --GSw_FINITO_Link={caseSwitches['GSw_FINITO_Link']}" + + f" --GSw_RetailAdder={caseSwitches['GSw_RetailAdder']}" + + f" --finito_inputs_dir={os.path.join('finito','inputs')}" + + f" --inputs_case_finito_dir={os.path.join('finito','inputs_case')}" + + f" --linked_report_dir={caseSwitches['linked_report_dir']} \n" + ) + # (linked) calls finito postprocessing (from [casedir]/finito/visualization) + OPATH.writelines( + f"python {os.path.join(casedir, 'finito', 'visualization', 'postprocessing.py')}" + + " -b 0" + + f" -l {caseSwitches['GSw_FINITO_Link']}" + + f' -c {batch_case} \n\n' + ) ### Run the retail rate module OPATH.writelines( @@ -1727,7 +1900,7 @@ def launch_single_case_run( shellscript = subprocess.Popen( [os.path.join(casedir, 'call_' + batch_case + ext)], shell=True) # Wait for it to finish before killing the thread - shellscript.wait() + #shellscript.wait() else: if int(caseSwitches['keep_run_terminal']) == 1: terminal_keep_flag = ' /k ' @@ -1855,7 +2028,7 @@ def main( help="Check inputs but don't start runs") args = parser.parse_args() - + main( BatchName=args.BatchName, cases_suffix=args.cases_suffix, single=args.single, simult_runs=args.simult_runs, forcelocal=args.forcelocal, skip_checks=args.skip_checks, From 0424508d93b090aad91d14afe9e1e064b15c36ec Mon Sep 17 00:00:00 2001 From: bsergi Date: Thu, 16 Apr 2026 12:42:43 -0400 Subject: [PATCH 02/15] branch transfer finito linkage --- ReEDS_Augur/prep_data.py | 32 +++- b_inputs.gms | 50 ++++-- c_supplymodel.gms | 88 ++++++++-- c_supplyobjective.gms | 49 +++--- cases.csv | 6 + cases_finito.csv | 46 +++++ createmodel.gms | 14 ++ d1_temporal_params.gms | 13 +- d2_varfix.gms | 6 + d3_data_dump.gms | 18 +- d_solveoneyear.gms | 19 +++ e_report.gms | 141 ++++++++++----- e_report_dump.py | 27 ++- e_report_params.csv | 1 + input_processing/copy_files.py | 5 + input_processing/hourly_load.py | 47 +++++ input_processing/hourly_repperiods.py | 27 ++- input_processing/hourly_writetimeseries.py | 3 + inputs/emission_constraints/co2_tax.csv | 84 ++++----- reeds/reedsplots.py | 1 + restart_runs.py | 14 +- runbatch.py | 189 ++++++++++++++++++++- 22 files changed, 716 insertions(+), 164 deletions(-) create mode 100644 cases_finito.csv diff --git a/ReEDS_Augur/prep_data.py b/ReEDS_Augur/prep_data.py index 8f4febfc..9805767d 100644 --- a/ReEDS_Augur/prep_data.py +++ b/ReEDS_Augur/prep_data.py @@ -147,8 +147,38 @@ def main(t, casedir, iteration=0): h_dt_szn.index.map(hmap_allyrs.set_index(['year', 'hour'])['*timestamp'])) h_dt_szn = h_dt_szn.reset_index().set_index('timestamp') + # load exogenous demand seen by ReEDS load = reeds.io.read_file(os.path.join(inputs_case, 'load.h5'), parse_timestamps=True) + load_year = load.loc[t] + + # when running linked model with FINITO (GSw_FINITO_Link=1), + # we also need to add in the load from FINITO + if int(sw.GSw_FINITO_Link): + load_finito_rt = gdxreeds['load_finito_rt'].rename(columns={'allh':'h', 'Value':'load_MW'}) + + # down select to relevant model year + load_finito = load_finito_rt.loc[load_finito_rt.t==t].drop('t', axis=1).copy() + #TODO: why so much shifting? validate the dynamics here with FINITO team + #TODO: also check against reference quantity + + # map from rep day to actual hour + # since we don't have multi-year profiles for FINITO load + # we assume they repeat across all weather years + # TODO: test that this works with multiple regions + load_finito = pd.merge(load_finito, h_dt_szn.reset_index(), on='h', how='outer')[['timestamp', 'r', 'load_MW']] + load_finito = load_finito.rename(columns={'timestamp':'datetime'}) + # add model year back in + #load_finito = load_finito.assign(year=t) + load_finito = load_finito.pivot(index=['datetime'], columns='r', values='load_MW') + + # convert timezone and fill any missing columns + load_finito.index = load_finito.index.tz_convert(load_year.index.tz) + load_finito = load_finito.reindex(columns=load_year.columns, fill_value=0) + + # add to load + load_new = load_year + load_finito + resources = pd.read_csv(os.path.join(inputs_case, 'resources.csv')) recf = reeds.io.read_file(os.path.join(inputs_case, 'recf.h5'), parse_timestamps=True) recf.columns = pd.MultiIndex.from_tuples([tuple(x.split('|')) for x in recf.columns], @@ -367,7 +397,7 @@ def intify(v): .merge(h_dt_szn_load_years[['h']].reset_index(), left_on='allh', right_on='h') .pivot(index='timestamp', columns='r', values='Value') ) - load_year = load.loc[t].add(can_exports, fill_value=0) + load_year = load_year.add(can_exports, fill_value=0) ### PRAS doesn't yet handle flexible load, so include all H2/DAC load in the ### version we write for PRAS diff --git a/b_inputs.gms b/b_inputs.gms index 0d59e562..3c135f24 100644 --- a/b_inputs.gms +++ b/b_inputs.gms @@ -1341,7 +1341,8 @@ set tmodel(t) "years to include in the model", tfix(t) "years to fix variables over when summing over previous years", tprev(t,tt) "previous modeled tt from year t", stfeas(st) "states to include in the model", - tsolved(t) "years that have solved" ; + tsolved(t) "years that have solved", + tfuel(t) "years that use ReEDS fuel supply curve module (otherwise uses supply curves in FINITO)" ; *following parameters get re-defined when the solve years have been declared parameter mindiff(t) "minimum difference between t and all other tt that are in tmodel(t)" ; @@ -1354,7 +1355,7 @@ tfix(t) = no ; stfeas(st) = no ; tprev(t,tt) = no ; tsolved(t) = no ; - +tfuel(t)=no; *============================== * Year specification @@ -1379,6 +1380,12 @@ tprev(t,tt)$[tmodel_new(t)$tmodel_new(tt)$(tt.valmindiff(t))] = no ; +* If FINITO linkage is on, remove all modeled years from tfuel after first_year_finito +* as the FINITO supply curves will be used instead of those in ReEDS; +* otherwise, all modeled years use ReEDS supply curves and are eligible for tfuel +tfuel(t)$[tmodel_new(t)]=yes; +tfuel(t)$[tmodel_new(t)$Sw_FINITO_Link$(t.val>=%first_year_finito%)] = no ; + * In order to fill all necessary dimensions of upgrade techs parameters, we require * Sw_UpgradeYear in ban(i) to be a modeled year and thus we compute as either * the GSw_UpgradeYear option or the next modeled years after GSw_UpgradeYear @@ -2331,6 +2338,7 @@ noncumulative_prescriptions(pcat,r,t)$tmodel_new(t) ], prescribednonrsc(tt,pcat,r,"value") + prescribedrsc(tt,pcat,r,"value") } ; +noncumulative_prescriptions(pcat,r,t)$noncumulative_prescriptions(pcat,r,t) = round(noncumulative_prescriptions(pcat,r,t),6); parameter noncumulative_prescriptions_energy(pcat,r,t) "--MWh-- prescribed energy capacity that comes online in a given year" ; noncumulative_prescriptions_energy(pcat,r,t)$tmodel_new(t) @@ -2777,6 +2785,7 @@ h2_ptc(i,v,r,t)$valcap(i,v,r,t) = h2_ptc_in(i,v,t) ; * Otherwise, we assume it receives $0/kg because the cleanliness of its carbon cannot be proven h2_ptc("electrolyzer",v,r,t)$[(not Sw_H2_PTC)] = 0; + set h2_ptc_years(t) "years in which the hydrogen production incentive is active"; h2_ptc_years(t) = tmodel_new(t)$[sum{(i,v,r),h2_ptc(i,v,r,t)}]; @@ -3571,6 +3580,7 @@ $include inputs_case%ds%tsc_binwidth.csv $offdelim $onlisting / ; +tsc_binwidth(r,rr,tscbin) = 1e9; parameter tsc_forward(r,rr,tscbin) "--$/MW-- transmission upgrade cost for forward direction" / @@ -6543,22 +6553,25 @@ $offdelim $ondigit $onlisting / , - min_co2_spurline_distance "--mi-- minimum distance for a spur line (used to provide a floor for pipeline distances in r_cs_distance)" + min_co2_spurline_distance "--mi-- minimum distance for a spur line (used to provide a floor for pipeline distances in r_cs_distance)", + min_r_cs_distance(r) "--mi-- mininum euclidean distance between BA transmission endpoints and storage formations" ; $offempty -* Wherever BA centroids fall within formation boundaries, assume some average spur line distance to connect a CCS or DAC plant with an injection site -min_co2_spurline_distance = 20 ; -r_cs_distance(r,cs)$[r_cs_distance(r,cs) < min_co2_spurline_distance] = min_co2_spurline_distance ; - -* Assign spurline costs -cost_co2_spurline_cap(r,cs,t)$[r_cs(r,cs)$tmodel_new(t)] = Sw_CO2_spurline_cost * r_cs_distance(r,cs) ; +* Wherever BA centroids fall within formation boundaries, r_cs_distance will be equal to 0, to ensure these are not dropped from the set, set these elements to a small number +r_cs_distance(r,cs)$[r_cs(r,cs)$(r_cs_distance(r,cs) = 0)] = 1e-3 ; +* find the closest site to each region +min_r_cs_distance(r) = smin(cs$[r_cs(r,cs)], r_cs_distance(r,cs)); -* CO2 pipelines can be build between any two adjacent BAs -cost_co2_pipeline_cap(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_cost * pipeline_distance(r,rr) ; -cost_co2_pipeline_fom(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_fom * pipeline_distance(r,rr) ; +$ifthene.rcslimit %GSw_CO2_LimitStorageSites% == 1 +* remove the region site combinations that are not the closest +r_cs_distance(r,cs)$[r_cs(r,cs)$(r_cs_distance(r,cs) <= min_r_cs_distance(r))] = min_r_cs_distance(r) ; +r_cs_distance(r,cs)$[r_cs(r,cs)$(r_cs_distance(r,cs) > min_r_cs_distance(r))] = 0 ; +$endif.rcslimit -co2_routes(r,rr)$routes_adjacent(r,rr) = yes ; +* Wherever BA centroids fall within formation boundaries, assume some average spur line distance to connect a CCS or DAC plant with an injection site +min_co2_spurline_distance = 20 ; +r_cs_distance(r,cs)$[r_cs_distance(r,cs)$(r_cs_distance(r,cs) <= min_co2_spurline_distance )] = min_co2_spurline_distance ; $onempty table co2_char(cs,*) "co2 site characteristics including injection rate limit, total storage limit, and break even cost" @@ -6577,9 +6590,18 @@ cost_co2_stor_bec(cs,t) = co2_char(cs,"bec_%GSw_CO2_BEC%"); csfeas(cs)$[co2_storage_limit(cs)$co2_injection_limit(cs)] = yes ; * only want to consider r_cs pairs which have available capacity r_cs(r,cs)$[not csfeas(cs)] = no ; +r_cs(r,cs)$[not r_cs_distance(r,cs)] = no ; +* Assign spurline costs +cost_co2_spurline_cap(r,cs,t)$[r_cs(r,cs)$tmodel_new(t)] = Sw_CO2_spurline_cost * r_cs_distance(r,cs) ; cost_co2_spurline_fom(r,cs,t)$[r_cs(r,cs)$tmodel_new(t)] = Sw_CO2_spurline_fom * r_cs_distance(r,cs) ; +* CO2 pipelines can be build between any two adjacent BAs +cost_co2_pipeline_cap(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_cost * pipeline_distance(r,rr) ; +cost_co2_pipeline_fom(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_fom * pipeline_distance(r,rr) ; + +co2_routes(r,rr)$routes_adjacent(r,rr) = yes ; + cost_co2_pipeline_cap(r,rr,t) = %GSw_CO2_CostAdj% * cost_co2_pipeline_cap(r,rr,t); cost_co2_pipeline_fom(r,rr,t) = %GSw_CO2_CostAdj% * cost_co2_pipeline_fom(r,rr,t); cost_co2_stor_bec(cs,t) = %GSw_CO2_CostAdj% * cost_co2_stor_bec(cs,t) ; @@ -6771,4 +6793,4 @@ emit_rate(etype,e,i,v,r,t)$[not valgen(i,v,r,t)] = 0 ; *============================================================ valinv_init(i,v,r,t) = valinv(i,v,r,t) ; -valcap_init(i,v,r,t) = valcap(i,v,r,t) ; +valcap_init(i,v,r,t) = valcap(i,v,r,t) ; \ No newline at end of file diff --git a/c_supplymodel.gms b/c_supplymodel.gms index 029d53f1..71389843 100644 --- a/c_supplymodel.gms +++ b/c_supplymodel.gms @@ -396,6 +396,14 @@ eq_loadcon(r,h,t)$tmodel(t).. * (the effect is the same but avoiding the h-indexed OP_LOADSITE reduces solve time) + OP_LOADSITE(r,h,t)$[Sw_LoadSiteCF$(Sw_LoadSiteCF<1)$val_loadsite(r)] + CAP_LOADSITE(r,t)$[(Sw_LoadSiteCF=1)$val_loadsite(r)] + +* [plus] load for industrial and converted fuel facilities (FINITO), +* including the PRM for stress periods +* TODO: should this include distloss? +* [MWh/hr = MW] +$ifthene.linked_load Sw_FINITO_Link==1 + + USE_ELE_FINITO(r,h,t)$[t_finito(t)] * (1 + prm(r,t)$h_stress(h)) +$endif.linked_load ; * --------------------------------------------------------------------------- @@ -1040,7 +1048,6 @@ eq_rsc_INVlim(r,i,rscbin,t)$[tmodel(t) *plus exogenous (pre-start-year) capacity, using its level in the first year (tfirst) + sum{(ii,v,tt)$[tfirst(tt)$rsc_agg(i,ii)$exog_rsc(i)], capacity_exog_rsc(ii,v,r,rscbin,tt) } - ; * --------------------------------------------------------------------------- @@ -2881,7 +2888,7 @@ eq_national_gen(t)$[tmodel(t)$national_gen_frac(t)$Sw_GenMandate].. * --------------------------------------------------------------------------- *gas used from each bin is the sum of all gas used -eq_gasused(cendiv,h,t)$[tmodel(t)$((Sw_GasCurve=0) or (Sw_GasCurve=3))].. +eq_gasused(cendiv,h,t)$[tmodel(t)$tfuel(t)$((Sw_GasCurve=0) or (Sw_GasCurve=3))].. sum{gb,GASUSED(cendiv,gb,h,t) } @@ -2898,7 +2905,7 @@ eq_gasused(cendiv,h,t)$[tmodel(t)$((Sw_GasCurve=0) or (Sw_GasCurve=3))].. * --------------------------------------------------------------------------- * gas from each bin needs to less than its capacity -eq_gasbinlimit(cendiv,gb,t)$[tmodel(t)$(Sw_GasCurve=0)].. +eq_gasbinlimit(cendiv,gb,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=0)].. gaslimit(cendiv,gb,t) @@ -2909,7 +2916,7 @@ eq_gasbinlimit(cendiv,gb,t)$[tmodel(t)$(Sw_GasCurve=0)].. * --------------------------------------------------------------------------- -eq_gasbinlimit_nat(gb,t)$[tmodel(t)$(Sw_GasCurve=3)].. +eq_gasbinlimit_nat(gb,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=3)].. gaslimit_nat(gb,t) @@ -2922,7 +2929,7 @@ eq_gasbinlimit_nat(gb,t)$[tmodel(t)$(Sw_GasCurve=3)].. * --------------------------------------------------------------------------- -eq_gasaccounting_regional(cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)].. +eq_gasaccounting_regional(cendiv,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)].. sum{fuelbin, VGASBINQ_REGIONAL(fuelbin,cendiv,t) } @@ -2935,7 +2942,7 @@ eq_gasaccounting_regional(cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)].. * --------------------------------------------------------------------------- -eq_gasaccounting_national(t)$[tmodel(t)$(Sw_GasCurve=1)].. +eq_gasaccounting_national(t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)].. sum{fuelbin,VGASBINQ_NATIONAL(fuelbin,t) } @@ -2948,7 +2955,7 @@ eq_gasaccounting_national(t)$[tmodel(t)$(Sw_GasCurve=1)].. * --------------------------------------------------------------------------- -eq_gasbinlimit_regional(fuelbin,cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)].. +eq_gasbinlimit_regional(fuelbin,cendiv,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)].. Gasbinwidth_regional(fuelbin,cendiv,t) @@ -2959,7 +2966,7 @@ eq_gasbinlimit_regional(fuelbin,cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)].. * --------------------------------------------------------------------------- -eq_gasbinlimit_national(fuelbin,t)$[tmodel(t)$(Sw_GasCurve=1)].. +eq_gasbinlimit_national(fuelbin,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)].. Gasbinwidth_national(fuelbin,t) @@ -2973,10 +2980,10 @@ eq_gasbinlimit_national(fuelbin,t)$[tmodel(t)$(Sw_GasCurve=1)].. *============================== * -- Bioenergy Supply Curve -- *============================== +* defer to FINITO representation when models are linked * --------------------------------------------------------------------------- - -eq_bioused(r,t)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }$tmodel(t)].. +eq_bioused(r,t)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }$tmodel(t)$tfuel(t)].. sum{bioclass, BIOUSED(bioclass,r,t) } @@ -2995,7 +3002,7 @@ eq_bioused(r,t)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }$tmodel(t)].. * --------------------------------------------------------------------------- * biomass consumption limit is annual -eq_biousedlimit(bioclass,usda_region,t)$tmodel(t).. +eq_biousedlimit(bioclass,usda_region,t)$[tmodel(t)$tfuel(t)].. biosupply(usda_region,bioclass,"cap") @@ -3565,6 +3572,11 @@ eq_h2_demand(p,t)$[(sameas(p,"H2"))$tmodel(t)$(yeart(t)>=h2_demand_start)$(Sw_H2 + sum{(i,v,r,h)$[valgen(i,v,r,t)$h2_combustion(i)$h_rep(h)], GEN(i,v,r,h,t) * hours(h) * h2_combustion_intensity * heat_rate(i,v,r,t) } + +* hydrogen demand from indusrty (FINITO): FINITO demand [MMBtu/yr] * conversion [metric tons-H2/MMBtu-H2] +$ifthene.linked_h2_nat Sw_FINITO_Link==1 + + [sum{(r,h)$h_rep(h), hours(h) * USE_H2_FINITO(r,h,t) * h2_metric_tons_per_mmbtu }]$t_finito(t) +$endif.linked_h2_nat ; * --------------------------------------------------------------------------- @@ -3596,6 +3608,11 @@ eq_h2_demand_regional(r,h,t) + sum{(i,v)$[valgen(i,v,r,t)$h2_combustion(i)], GEN(i,v,r,h,t) * h2_combustion_intensity * heat_rate(i,v,r,t) } + +* regional hydrogen demand for industry (FINITO) +$ifthene.linked_h2_reg Sw_FINITO_Link==1 + + [ USE_H2_FINITO(r,h,t) * h2_metric_tons_per_mmbtu ]$t_finito(t) +$endif.linked_h2_reg ; * --------------------------------------------------------------------------- @@ -3838,11 +3855,16 @@ eq_co2_capture(r,h,t) * capture from DAC + sum{(i,v)$[dac(i)$valcap(i,v,r,t)$i_p(i,"DAC")], PRODUCE("DAC",i,v,r,h,t) }$Sw_DAC +* (linked) capture from industry [metric_tons-CO2/hr]: +* calculation: hours_per_year [yrs/hr] * capture [scaled_metric_tons-CO2/yr] / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2] +$ifthene.linked_co2_capture Sw_FINITO_Link==1 + + [CAPTURE_CO2EM(r,h,t) / co2_scale]$[t_finito(t)$h_rep(h)] +$endif.linked_co2_capture ; * --------------------------------------------------------------------------- -eq_co2_transport_caplimit(r,rr,h,t)$[co2_routes(r,rr)$Sw_CO2_Detail +eq_co2_transport_caplimit(r,rr,h,t)$[co2_routes(r,rr)$Sw_CO2_Detail$h_rep(h) $tmodel(t)$(yeart(t)>=co2_detail_startyr)].. *capacity computed as cumulative investments of co2 pipelines up to the current year @@ -3858,7 +3880,8 @@ eq_co2_transport_caplimit(r,rr,h,t)$[co2_routes(r,rr)$Sw_CO2_Detail * --------------------------------------------------------------------------- -eq_co2_spurline_caplimit(r,cs,h,t)$[Sw_CO2_Detail$r_cs(r,cs)$tmodel(t)$(yeart(t)>=co2_detail_startyr)].. +eq_co2_spurline_caplimit(r,cs,h,t)$[Sw_CO2_Detail$r_cs(r,cs)$h_rep(h)$tmodel(t) + $(yeart(t)>=co2_detail_startyr)].. *capacity computed as cumulative investments of co2 spurlines up to the current year sum{tt$[(yeart(tt)<=yeart(t))$(tmodel(tt) or tfix(tt))$(yeart(tt)>=co2_detail_startyr)], @@ -3867,11 +3890,16 @@ eq_co2_spurline_caplimit(r,cs,h,t)$[Sw_CO2_Detail$r_cs(r,cs)$tmodel(t)$(yeart(t) =g= CO2_STORED(r,cs,h,t) +* (linked) extraction of CO2 [metric_tons-CO2/hr] +* calculation: hours_per_year [yrs/hr] * (1 / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2]) * use [scaled_metric_tons-CO2/yr] +$ifthene.linked_co2_spurline_caplimit Sw_FINITO_Link==1 + + [ (1 / co2_scale) * EXTRACT_CO2_CS(cs,r,h,t) ]$[h_rep(h)$t_finito(t)] +$endif.linked_co2_spurline_caplimit ; * --------------------------------------------------------------------------- -eq_co2_sink(r,h,t)$[tmodel(t)$Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)].. +eq_co2_sink(r,h,t)$[tmodel(t)$Sw_CO2_Detail$h_rep(h)$(yeart(t)>=co2_detail_startyr)].. *the amount of co2 stored from r in all of its cs sites sum{cs$r_cs(r,cs), CO2_STORED(r,cs,h,t) } @@ -3886,11 +3914,24 @@ eq_co2_sink(r,h,t)$[tmodel(t)$Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)].. * net trade + sum{rr$co2_routes(r,rr), CO2_FLOW(rr,r,h,t) - CO2_FLOW(r,rr,h,t) } + +* (linked) extraction - use of CO2 [metric_tons-CO2/hr] +*calculation: hours_per_year [yrs/hr] * (1 / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2]) * extraction/use [scaled_metric_tons-CO2/yr] +$ifthene.linked_co2_sink Sw_FINITO_Link==1 + + (1 / co2_scale) * [ +* extraction from all cs sites in r + + sum{cs$[csfeas(cs)$r_cs(r,cs)], EXTRACT_CO2_CS(cs,r,h,t)} +* total use of CO2, equivalent to supply + - USE_CO2(r,h,t) + ]$[h_rep(h)$t_finito(t)] +$endif.linked_co2_sink + ; * --------------------------------------------------------------------------- -eq_co2_injection_limit(cs,h,t)$[Sw_CO2_Detail$tmodel(t)$(yeart(t)>=co2_detail_startyr)$csfeas(cs)].. +eq_co2_injection_limit(cs,h,t)$[tmodel(t)$Sw_CO2_Detail$h_rep(h) + $(yeart(t)>=co2_detail_startyr)$csfeas(cs)].. * exogenously defined injection limit co2_injection_limit(cs) @@ -3899,6 +3940,14 @@ eq_co2_injection_limit(cs,h,t)$[Sw_CO2_Detail$tmodel(t)$(yeart(t)>=co2_detail_st * must exceed metric tons per hour entering storage sum{r$r_cs(r,cs), CO2_STORED(r,cs,h,t) } + + + +* (linked) extraction of CO2 for use [metric_tons-CO2/hr] +* calculation: hours_per_year [yrs/hr] * (1 / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2]) * use [scaled_metric_tons-CO2/yr] +$ifthene.linked_co2_injection_limit Sw_FINITO_Link==1 + + (1 / co2_scale) * sum{r$[r_cs(r,cs)], EXTRACT_CO2_CS(cs,r,h,t) }$[h_rep(h)$t_finito(t)] +$endif.linked_co2_injection_limit ; * --------------------------------------------------------------------------- @@ -3915,7 +3964,16 @@ eq_co2_cumul_limit(cs,t)$[tmodel(t)$Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr) $[(yeart(tt)<=yeart(t))$(tmodel(tt) or tfix(tt))$(yeart(tt)>=co2_detail_startyr) $r_cs(r,cs)$h_rep(h)], yearweight(tt) * hours(h) * CO2_STORED(r,cs,h,tt) } + +* (linked) cumulative amount extracted over time +$ifthene.linked_co2_storage_cumul_limit Sw_FINITO_Link==1 + - sum{(r,h,tt) + $[(yeart(tt)<=yeart(t))$(tmodel(tt) or tfix(tt))$(yeart(tt)>=co2_detail_startyr) + $r_cs(r,cs)$(tfuel(tt))], + yearweight(tt) * hours(h) * EXTRACT_CO2_CS(cs,r,h,tt) }$[t_finito(t)] +$endif.linked_co2_storage_cumul_limit ; + * --------------------------------------------------------------------------- *=================== diff --git a/c_supplyobjective.gms b/c_supplyobjective.gms index 38a5b0b3..a5955399 100644 --- a/c_supplyobjective.gms +++ b/c_supplyobjective.gms @@ -19,7 +19,17 @@ Variable Z "--$-- total cost of operations and investment, scale varie * objective function is the sum over modeled years of the investment * and operations components -eq_ObjFn.. Z =e= cost_scale * sum{t$tmodel(t), Z_inv(t) + Z_op(t) } ; +eq_ObjFn.. Z =e= cost_scale * ( +* electricity and H2 costs + sum{t$tmodel(t), Z_inv(t) + Z_op(t) } +* economy-wide costs from FINITO: deflate from $2018 to $2004 +* and remove any FINITO scaling +$ifthene.linked_objective Sw_FINITO_Link==1 + + deflator('2018') / obj_scale + * sum{t$[t_finito(t)], Z_finito(t) } +$endif.linked_objective + ) +; *======================================================= * -- Investment component of the objective function -- @@ -221,7 +231,7 @@ eq_Objfn_op(t)$tmodel(t).. * plus cost of H2 fuel when using fixed price (Sw_H2=0) or during stress periods. * When using endogenous H2 price (Sw_H2=1 or Sw_H2=2), H2 fuel cost is captured elsewhere * via the capex + opex costs of H2 production and its associated electricity demand. - + sum{(i,v,r,h)$[valgen(i,v,r,t)$heat_rate(i,v,r,t) + + sum{(i,v,r,h)$[tfuel(t)$valgen(i,v,r,t)$heat_rate(i,v,r,t) $(not gas(i))$(not bio(i))$(not cofire(i)) $((not h2_combustion(i)) or h2_combustion(i)$[(Sw_H2=0) or h_stress(h)])], hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN(i,v,r,h,t) } @@ -232,50 +242,49 @@ eq_Objfn_op(t)$tmodel(t).. * --cofire coal consumption--- * cofire bio consumption already accounted for in accounting of BIOUSED - + sum{(i,v,r,h)$[valgen(i,v,r,t)$cofire(i)$heat_rate(i,v,r,t)], + + sum{(i,v,r,h)$[tfuel(t)$valgen(i,v,r,t)$cofire(i)$heat_rate(i,v,r,t)], (1-bio_cofire_perc) * hours(h) * heat_rate(i,v,r,t) * fuel_price("coal-new",r,t) * GEN(i,v,r,h,t) } * --- cost of natural gas--- *Sw_GasCurve = 2 (static natural gas prices) *first - gas consumed for electricity generation - + sum{(i,v,r,h)$[valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)$(Sw_GasCurve = 2)], + + sum{(i,v,r,h)$[tfuel(t)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)$(Sw_GasCurve = 2)], hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN(i,v,r,h,t) } *second - gas consumed by gas-powered DAC - + sum{(v,r,h)$[valcap("dac_gas",v,r,t)$(Sw_GasCurve = 2)], - hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE("DAC","dac_gas",v,r,h,t) }$Sw_DAC_Gas + + sum{(v,r,h)$[tfuel(t)$valcap("dac_gas",v,r,t)$(Sw_GasCurve = 2)], + hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE("DAC","dac_gas",v,r,h,t) }$[Sw_DAC_Gas] *Sw_GasCurve = 0 (census division supply curves natural gas prices) - + sum{(cendiv,gb), sum{h, hours(h) * GASUSED(cendiv,gb,h,t) } + + sum{(cendiv,gb)$[tfuel(t)], sum{h, hours(h) * GASUSED(cendiv,gb,h,t) } * gasprice(cendiv,gb,t) - }$(Sw_GasCurve = 0) + }$[(Sw_GasCurve = 0)] *Sw_GasCurve = 3 (national supply curve for natural gas prices with census division multipliers) - + sum{(h,cendiv,gb), hours(h) * GASUSED(cendiv,gb,h,t) + + sum{(h,cendiv,gb)$[tfuel(t)], hours(h) * GASUSED(cendiv,gb,h,t) * gasadder_cd(cendiv,t,h) + gasprice_nat_bin(gb,t) - }$(Sw_GasCurve = 3) + }$[(Sw_GasCurve = 3)] *Sw_GasCurve = 1 (national and census division supply curves for natural gas prices) *first - anticipated costs of gas consumption given last year's amount - + (sum{(i,r,v,cendiv,h)$[valgen(i,v,r,t)$gas(i)], + + (sum{(i,r,v,cendiv,h)$[valgen(i,v,r,t)$gas(i)$tfuel(t)], gasmultterm(cendiv,t) * szn_adj_gas(h) * cendiv_weights(r,cendiv) * hours(h) * heat_rate(i,v,r,t) * GEN(i,v,r,h,t) } *second - adjustments based on changes from last year's consumption at the regional and national level - + sum{(fuelbin,cendiv), + + sum{(fuelbin,cendiv)$tfuel(t), gasbinp_regional(fuelbin,cendiv,t) * VGASBINQ_REGIONAL(fuelbin,cendiv,t) } - + sum{(fuelbin), + + sum{(fuelbin)$tfuel(t), gasbinp_national(fuelbin,t) * VGASBINQ_NATIONAL(fuelbin,t) } - )$[Sw_GasCurve = 1] + )$[(Sw_GasCurve = 1)] * ---cost of biofuel consumption and biomass transport--- - + sum{(r,bioclass)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }], + + sum{(r,bioclass)$[tfuel(t)$sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }], BIOUSED(bioclass,r,t) * (sum{usda_region$r_usda(r,usda_region), biosupply(usda_region, bioclass, "price") } + bio_transport_cost) } - * --- hurdle costs for transmission flow --- + sum{(r,rr,h,trtype)$[routes(r,rr,trtype,t)$cost_hurdle(r,rr,t)], cost_hurdle(r,rr,t) * FLOW(r,rr,h,t,trtype) * hours(h) } @@ -340,15 +349,15 @@ eq_Objfn_op(t)$tmodel(t).. CO2_SPURLINE_INV(r,cs,tt) } }$[Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)] * --- CO2 injection break even costs - + sum{(r,cs,h)$r_cs(r,cs), hours(h) * CO2_STORED(r,cs,h,t) * cost_co2_stor_bec(cs,t) }$[Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)] + + sum{(r,cs,h)$[r_cs(r,cs)$h_rep(h)], hours(h) * CO2_STORED(r,cs,h,t) * cost_co2_stor_bec(cs,t) }$[Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)] * --- Tax credit for CO2 stored --- * note conversion to 12-year CRF given length of CO2 captured incentive payments - - sum{(i,v,r,h)$[valgen(i,v,r,t)$co2_captured_incentive(i,v,r,t)], + - sum{(i,v,r,h)$[valgen(i,v,r,t)$h_rep(h)$co2_captured_incentive(i,v,r,t)], (crf(t) / crf_co2_incentive(t)) * co2_captured_incentive(i,v,r,t) * hours(h) * capture_rate("CO2",i,v,r,t) * GEN(i,v,r,h,t)} * --- Tax credit for CO2 stored for DAC --- - - sum{(p,i,v,r,h)$[dac(i)$valcap(i,v,r,t)$i_p(i,p)$h_rep(h)], + - sum{(p,i,v,r,h)$[dac(i)$valcap(i,v,r,t)$i_p(i,p)$h_rep(h)$co2_captured_incentive(i,v,r,t)], (crf(t) / crf_co2_incentive(t)) * co2_captured_incentive(i,v,r,t) * hours(h) * PRODUCE(p,i,v,r,h,t)} * --- PTC value for electric power generation --- @@ -362,7 +371,7 @@ eq_Objfn_op(t)$tmodel(t).. - sum{(p,v,r,h)$[valcap("electrolyzer",v,r,t)$(sameas(p,"H2"))$h2_ptc("electrolyzer",v,r,t)$h_rep(h)], hours(h) * PRODUCE(p,"electrolyzer",v,r,h,t) * (crf(t) / crf_h2_incentive(t)) * h2_ptc("electrolyzer",v,r,t) * 1e3} - $[Sw_H2_PTC$Sw_H2$h2_ptc_years(t)$(yeart(t) >= h2_demand_start)] + $[(Sw_H2_PTC)$Sw_H2$h2_ptc_years(t)$(yeart(t) >= h2_demand_start)] *end multiplier for pvf_onm ) diff --git a/cases.csv b/cases.csv index 9d0752e7..59f86c65 100644 --- a/cases.csv +++ b/cases.csv @@ -33,6 +33,11 @@ coalscen,"Coal price scenario (inputs\fuelprices\{coalscen}.csv). The options fo demandscen,Demand scenario (see inputs\load),N/A,AEO_2025_reference, dr_shedscen,DR shed scenario,demo_data_January_2025,demo_data_January_2025, evmcscen,EVMC scenario for load and supply curve and flexibility,Baseline,Baseline, +GSw_FINITO_Link,Switch to enable linkage with FINITO,0; 1,0, +finito_dir,Directory with cloned FINITO repository,N/A,/projects/finito/molmezt/FINITO, +finito_cases_file,"extension of the FINITO cases file including all scenarios to be run, e.g., 'linked' for cases_linked.csv",N/A,linked, +finito_case,Case name for FINITO in {finito_dir}/cases_linked.csv; use 'same' if the scenario name in cases_linked.csv matches the scenario name in the cases file,N/A,BAU, +first_year_finito,First year when the FINITO module is enabled,N/A,2018, geodiscov,Annual discover rates for new geothermal sites,BAU; TI,BAU, geohydrosupplycurve,Geohydro Supply Curve,ATB_2023,ATB_2023, egssupplycurve,EGS Supply Curve,ATB_2023; reV,reV, @@ -116,6 +121,7 @@ GSw_CO2_pipeline_fom,Set FOM cost for CO2 Pipeline. 2004$/(tonne-mi)/hr-yr. Sugg GSw_CO2_spurline_cost,Set capex cost for CO2 Spurline. 2004$/(tonne-mi)/hr. Suggested value: 2257,N/A,2257, GSw_CO2_spurline_fom,Set FOM cost for CO2 Spurline. 2004$/(tonne-mi)/hr-yr. Suggested value: 60,N/A,60, GSw_CO2_Detail,Switch to enable detailed representation CO2 transportation and storage,0; 1,0, +GSw_CO2_LimitStorageSites,Switch to limit storage sites to the closest site for each region,0; 1,0, GSw_CO2_CostAdj,multiplier for co2 infrastructure when Sw_CO2_Detail = 1,float,1, GSw_CO2_BEC,"Break even cost capacity factor assignment. Must be a suffix to ""bec_"" in the columns of co2_site_char.csv",N/A,70, GSw_CoalIGCC,Turn on/off COAL IGCC,0; 1,1, diff --git a/cases_finito.csv b/cases_finito.csv new file mode 100644 index 00000000..25fb39db --- /dev/null +++ b/cases_finito.csv @@ -0,0 +1,46 @@ +,Default Value,ReEDSonly,BAUreg,Reference,Reference-AEO,Reference-HOG,TradeReference,TradeReference-AEO,Trade-HOG,Reference-LowBio,Reference-HighBio,Tax200,Tax200-AEO,NetZero,NetZero-LowBio,NetZero-HighBio,NetZero-SlackLim1500,NetZero-SlackLim1500-CO2Stor-Mid,NetZero-SlackLim1500-CO2Stor-Low,NetZero-SlackLim500-CO2Stor-Low,NetZero-SlackCost2500-CO2Stor-Low,aggreg_54_level +ignore,1,1,0,1,1,1,1,1,1,,,,,,,,,,,,, +GSw_FINITO_Link,1,0,,,,,,,,,,,,,,,,,,,, +finito_dir,/path/to/FINITO/directory,0,,,,,,,,,,,,,,,,,,,, +finito_cases_file,linked,,,,,,,,,,,,,,,,,,,,, +finito_case,BAU,,BAUreg,,BAU-AEO,AEO-HOG,TradeReEDS,TradeReEDS-AEO,TradeReEDS-HOG,BAU-LowBio,BAU-HighBio,,BAU-AEO,NetZero,NetZero-LowBio,NetZero-HighBio,NetZero-SlackLim1500,NetZero-SlackLim1500-CO2Stor-Mid,NetZero-SlackLim1500-CO2Stor-Low,NetZero-SlackLim500-CO2Stor-Low,NetZero-SlackCost2500-CO2Stor-Low,aggreg_54_level +coalscen,AEO_2025_reference,,,,,,,,,,,,,,,,,,,,, +ngscen,AEO_2025_reference,,,,,AEO_2025_HOG,,,AEO_2025_HOG,,,,,,,,,,,,, +uraniumscen,AEO_2025_reference,,,,,,,,,,,,,,,,,,,,, +GSw_Region,country/USA,,st/CO,,,,,,,,,,,,,,,,,,, +GSw_GasCurve,1,,2,,,,,,,,,,,,,,,,,,, +GSw_gopt,2,,,,,,,,,,,,,,,,,,,,, +GSw_NewValCapShrink,1,,,,,,,,,,,,,,,,,,,,, +endyear,2050,,2035,,,,,,,,,,,,,,,,,,, +yearset,2010_2015_2020..2050..5,,,,,,,,,,,,,,,,,,,,, +first_year_finito,2018,,,,,,,,,,,,,,,,,,,,, +incentives_suffix,obbba,,,,,,,,,,,,,,,,,,,,, +GSw_Upgrades,0,,,,,,,,,,,,,,,,,,,,, +GSw_AnnualCap,0,,,,,,,,,,,,,,,,,,,,, +GSw_AnnualCapScen,start2023_100pct2035,,,,,,,,,,,,,,,,,,,,, +GSw_Clean_Air_Act,0,,,,,,,,,,,,,,,,,,,,, +GSw_CSP,0,,,,,,,,,,,,,,,,,,,,, +GSw_Geothermal,0,,,,,,,,,,,,,,,,,,,,, +GSw_H2,2,,,,,,,,,,,,,,,,,,,,, +GSw_H2_PTC,0,,,,,,,,,,,,,,,,,,,,, +GSw_H2_Demand_Case,BAU,,,,,,,,,,,,,,,,,,,,, +GSw_H2_SMR,1,,,,,,,,,,,,,,,,,,,,, +GSw_H2_Transport,0,,,,,,,,,,,,,,,,,,,,, +GSw_H2Combustion,1,,,,,,,,,,,,,,,,,,,,, +GSw_H2Combustionupgrade,1,,,,,,,,,,,,,,,,,,,,, +GSw_H2CombinedCycle,1,,,,,,,,,,,,,,,,,,,,, +GSw_DAC,1,,,,,,,,,,,,,,,,,,,,, +GSw_BECCS,1,,,,,,,,,,,,,,,,,,,,, +GSw_RegionResolution,ba,,,,,,,,,,,,,,,,,,,,,aggreg +GSw_TCPhaseout,1,,,,,,,,,,,,,,,,,,,,, +GSw_CO2_Detail,1,,,,,,,,,,,,,,,,,,,,, +GSw_CO2_LimitStorageSites,1,,,,,,,,,,,,,,,,,,,,, +numhintage,4,,,,,,,,,,,,,,,,,,,,, +GSw_CarbTax,0,,,,,,,,,,,1,1,,,,,,,,, +GSw_CarbTaxOption,default,,,,,,,,,,,t200_2050,t200_2050,,,,,,,,, +GSw_PRM_CapCredit,0,,,,,,,,,,,,,,,,,,,,, +GSw_PRM_StressIterateMax,0,,,,,,,,,,,,,,,,,,,,, +GSw_HourlyNumClusters,9,,,,,,,,,,,,,,,,,,,,, +GSw_HourlyChunkLengthRep,4,,,,,,,,,,,,,,,,,,,,, +GSw_HourlyChunkLengthStress,4,,,,,,,,,,,,,,,,,,,,, +GSw_HierarchyFile,default,,,,,,,,,,,,,,,,,,,,,agg54 diff --git a/createmodel.gms b/createmodel.gms index 36ab4603..d8ba1dca 100644 --- a/createmodel.gms +++ b/createmodel.gms @@ -1,5 +1,19 @@ $include b_inputs.gms + +$ifthene.finito_link %GSw_FINITO_Link% == 1 +$include finito/model/finito_input.gms +$endif.finito_link + +$ifthene.finito_link %GSw_FINITO_Link% == 1 +$include finito/model/finito_vars_eqs.gms +$endif.finito_link + $include c_supplymodel.gms + +$ifthene.finito_link %GSw_FINITO_Link% == 1 +$include finito/model/finito_model.gms +$endif.finito_link + $include c_supplyobjective.gms $include c_mga.gms $include d_solveprep.gms diff --git a/d1_temporal_params.gms b/d1_temporal_params.gms index 444c7fc5..fcbf8578 100644 --- a/d1_temporal_params.gms +++ b/d1_temporal_params.gms @@ -299,9 +299,6 @@ can_exports_h(r,h,t)$[hours(h)] = can_exports(r,t) * can_exports_h_frac(h) / hou $endif.Canada -* zero Canada exports when Canada is not modeled -can_imports_szn(r,szn,t)$[Sw_Canada=0] = 0 ; -can_exports_h(r,h,t)$[Sw_Canada=0] = 0 ; $onempty parameter canmexload(r,allh) "load for canadian and mexican regions" @@ -887,12 +884,20 @@ szn_adj_gas(h)$frac_h_quarter_weights(h,"wint") = szn_adj_gas(h) + frac_h_quarter_weights(h,"wint") * szn_adj_gas_winter ; +*============================================= +* ----- ReEDS-FINITO temporal parameters ----- +*============================================= +* remove focused sectors load from load_exog_static +$ifthene.linked_finito_temporal_params %GSw_FINITO_Link% ==1 +$include finito/model/linked_finito_temporal_params.gms +$endif.linked_finito_temporal_params + + *============================================= * -- Round parameters for GAMS -- *============================================= avail(i,r,h)$avail(i,r,h) = round(avail(i,r,h),3) ; -can_imports_szn(r,szn,t)$can_imports_szn(r,szn,t) = round(can_imports_szn(r,szn,t),3) ; can_exports_h(r,h,t)$can_exports_h(r,h,t) = round(can_exports_h(r,h,t),3) ; h_weight_csapr(h)$h_weight_csapr(h) = round(h_weight_csapr(h),3) ; load_exog(r,h,t)$load_exog(r,h,t) = round(load_exog(r,h,t),3) ; diff --git a/d2_varfix.gms b/d2_varfix.gms index 85db2442..3a523732 100644 --- a/d2_varfix.gms +++ b/d2_varfix.gms @@ -9,6 +9,10 @@ if(Sw_RemoveSmallNumbers = 1, INV_ENERGY.l(i,v,r,tfix)$[abs(INV_ENERGY.l(i,v,r,tfix)) < rhs_tolerance] = 0 ; INV_RSC.l(i,v,r,rscbin,tfix)$[abs(INV_RSC.l(i,v,r,rscbin,tfix)) < rhs_tolerance] = 0 ; INV_POI.l(r,tfix)$[abs(INV_POI.l(r,tfix)) < rhs_tolerance] = 0 ; + INV_REFURB.l(i,v,r,tfix)$[abs(INV_REFURB.l(i,v,r,tfix)) < rhs_tolerance] = 0; + CO2_SPURLINE_INV.l(r,cs,tfix)$[abs(CO2_SPURLINE_INV.l(r,cs,tfix)) < rhs_tolerance] = 0; + CO2_TRANSPORT_INV.l(r,rr,tfix)$[abs(CO2_TRANSPORT_INV.l(r,rr,tfix)) < rhs_tolerance] = 0; + H2_TRANSPORT_INV.l(r,rr,tfix)$[abs(H2_TRANSPORT_INV.l(r,rr,tfix)) < rhs_tolerance] = 0 ; H2_STOR_INV.l(h2_stor,r,tfix)$[abs(H2_STOR_INV.l(h2_stor,r,tfix)) < rhs_tolerance] = 0 ; H2_TRANSPORT_INV.l(r,rr,tfix) $[abs(H2_TRANSPORT_INV.l(r,rr,tfix) ) < rhs_tolerance] = 0 ; ); @@ -119,7 +123,9 @@ H2_STOR_LEVEL_SZN.fx(h2_stor,r,actualszn,tfix)$[(h2_stor_r(h2_stor,r))$(Sw_H2=2) *CO2-related variables CO2_CAPTURED.fx(r,h,tfix)$Sw_CO2_Detail = CO2_CAPTURED.l(r,h,tfix) ; +$ifthene.linked_varfix Sw_FINITO_Link==0 CO2_STORED.fx(r,cs,h,tfix)$[Sw_CO2_Detail$r_cs(r,cs)] = CO2_STORED.l(r,cs,h,tfix) ; +$endif.linked_varfix CO2_FLOW.fx(r,rr,h,tfix)$[Sw_CO2_Detail$co2_routes(r,rr)] = CO2_FLOW.l(r,rr,h,tfix) ; CO2_TRANSPORT_INV.fx(r,rr,tfix)$[Sw_CO2_Detail$co2_routes(r,rr)] = CO2_TRANSPORT_INV.l(r,rr,tfix) ; CO2_SPURLINE_INV.fx(r,cs,tfix)$[Sw_CO2_Detail$r_cs(r,cs)] = CO2_SPURLINE_INV.l(r,cs,tfix) ; diff --git a/d3_data_dump.gms b/d3_data_dump.gms index 29c97036..30a8ff0b 100644 --- a/d3_data_dump.gms +++ b/d3_data_dump.gms @@ -8,6 +8,7 @@ This file creates a gdx file with all of the data necessary for the Augur module - availability rates (1 - outage rates) - transmission capacities and loss rates - technology sets + - industrial electricity demand from FINITO when running linked model $offtext $if not set start_year $setglobal start_year %startyear% @@ -62,6 +63,7 @@ h2_usage_regional(r,allh,t) "--metric tons-- H2 usage by region" inv_cond_filt(i,v,t) "--set-- vintage-year mapping for investments by technology" inv_ivrt(i,v,r,t) "--MW-- investments in power generation capacity" inv_energy_ivrt(i,v,r,t) "--MWh-- investments in energy generation capacity" +load_finito_rt(r,allh,t) "--MWh-- total load from industrial sectors and converted fuels modeled by FINITO" m_cf_filt(i,v,r,allh) "--fraction-- capacity factor used in the model" m_cf_szn_filt(i,v,r,allszn) "--fraction-- modelled capacity factors filtered for hydro resources to set seasonal energy constraints" minloadfrac_filt(r,i,allszn) "--fraction-- modelled mingen fraction filtered for hydro resources to set mingen constraints" @@ -161,6 +163,10 @@ h2_usage_regional(r,h,t)$tcur(t) = h2_exogenous_demand_regional(r,'h2',h,t) + sum{(i,v)$[valgen(i,v,r,t)$h2_combustion(i)], GEN.l(i,v,r,h,t) * h2_combustion_intensity * heat_rate(i,v,r,t)} +* regional hydrogen demand for industry (FINITO) +$ifthene.linked_h2_reg_report Sw_FINITO_Link==1 + + [USE_H2_FINITO.l(r,h,t) * h2_metric_tons_per_mmbtu ]$t_finito(t) +$endif.linked_h2_reg_report ) ; @@ -284,15 +290,24 @@ cap_converter_filt(r) = sum{t$tcur(t), CAP_CONVERTER.l(r,t) } ; routes_filt(r,rr,trtype) = sum{t$tcur(t), routes(r,rr,trtype,t) } ; *============================ -* Flexible load data +* Load data *============================ +* flexible load flex_load(r,h) = sum{(flex_type,t)$tcur(t), load_exog_flex(flex_type,r,h,t) } ; flex_load_opt(r,h) = sum{(flex_type,t)$tcur(t), FLEX.l(flex_type,r,h,t) } ; ra_cap_loadsite(r,t)$[Sw_LoadSiteCF$val_loadsite(r)] = CAP_LOADSITE.l(r,t) ; +* FINITO load +$ifthene.linked_load Sw_FINITO_Link==1 + load_finito_rt(r,h,t) = USE_ELE_FINITO.l(r,h,t) ; +$else.linked_load + load_finito_rt(r,h,t) = 0 ; +$endif.linked_load + + *============================ * Extra consumption data *============================ @@ -382,6 +397,7 @@ execute_unload 'ReEDS_Augur%ds%augur_data%ds%reeds_data_%cur_year%.gdx' inv_ivrt inv_energy_ivrt ivt_num + load_finito_rt m_cf_filt m_cf_szn_filt maxage diff --git a/d_solveoneyear.gms b/d_solveoneyear.gms index b3484a80..a62eaaee 100644 --- a/d_solveoneyear.gms +++ b/d_solveoneyear.gms @@ -23,6 +23,13 @@ tload(t) = no ; tmodel(t) = no ; tmodel("%cur_year%") = yes ; +* (linked) also reset t_finito +$ifthene.linked_finito_time %GSw_FINITO_Link% == 1 +t_finito(t) = no ; +t_finito(t)$[tmodel(t)$(t.val>=first_year_finito)$(t.val<=endyear)] = yes ; +$endif.linked_finito_time + + $log 'Solving sequential case for...' $log ' Case: %case%' $log ' Year: %cur_year%' @@ -241,6 +248,7 @@ $endif.diagnose * ------------------------------ * Solve the Model * ------------------------------ +option limrow = 1e6; $ifthen.valstr %GSw_ValStr% == 1 OPTION lp = convert ; ReEDSmodel.optfile = 1 ; @@ -276,11 +284,22 @@ $endif.mga *** Adjust some parameters based on the solution for this solve year $include d2_post_solve_adjustments.gms +*** Do the same for finito if necessary +$ifthene.linked_finito_shrink %GSw_FINITO_Link% == 1 +* Industrial variable fix +$include finito/model/finito_post_solve_adjustments.gms +$endif.linked_finito_shrink *** Fix decision variables to their optimized levels for this solve year tfix("%cur_year%") = yes ; $include d2_varfix.gms +*** Fix FINITO decision variables to their optimized levels for this solve year +$ifthene.linked_finito_varfix %GSw_FINITO_Link% == 1 +* Industrial variable fix +$include finito/model/finito_varfix.gms +$endif.linked_finito_varfix + *** Dump data used in calculations between solve years $include d3_data_dump.gms diff --git a/e_report.gms b/e_report.gms index ebd14d10..f251e657 100644 --- a/e_report.gms +++ b/e_report.gms @@ -107,7 +107,7 @@ rev_cat "categories for renvenue streams" /load, res_marg, oper_res, rps, charge lcoe_cat "categories for LCOE calculation" /capcost, upgradecost, rsccost, fomcost, vomcost, gen / -loadtype "categories for types of load" / end_use, dist_loss, trans_loss, stor_charge, h2_prod, h2_network, dac / +loadtype "categories for types of load" / end_use, dist_loss, trans_loss, stor_charge, h2_prod, h2_network, dac, finito / h2_demand_type / "electricity", "cross-sector"/ @@ -495,92 +495,121 @@ ptc_out(i,v,t)$[tmodel_new(t)$ptc_value_scaled(i,v,t)] = ptc_value_scaled(i,v,t) * Case 2: the resource of one or more biomass classes ARE exhausted, i.e., BIOUSED.l(bioclass) = biosupply(bioclass) * Marginal Biomass Price = maximum difference between eq_bioused.m and eq_biousedlimit.m(bioclass) across all biomass classes in a region -repbioprice(r,t)$tmodel_new(t) = max{0, smax{bioclass$BIOUSED.l(bioclass,r,t), eq_bioused.m(r,t) - - sum{usda_region$r_usda(r,usda_region), eq_biousedlimit.m(bioclass,usda_region,t) } } } / pvf_onm(t) ; +repbioprice(r,t)$[tmodel_new(t)$tfuel(t)] = + [max{0, smax{bioclass$BIOUSED.l(bioclass,r,t), eq_bioused.m(r,t) - + sum{usda_region$r_usda(r,usda_region), eq_biousedlimit.m(bioclass,usda_region,t) } } } / pvf_onm(t) ] ; -* quantity of biomass used (convert from mmBTU to dry tons using biomass energy content) -bioused_out(bioclass,r,t)$tmodel_new(t) = BIOUSED.l(bioclass,r,t) / bio_energy_content ; -bioused_usda(bioclass,usda_region,t)$tmodel_new(t) = sum{r$r_usda(r,usda_region), bioused_out(bioclass,r,t) } ; +* quantity of biomass used in the power sector (convert from mmBTU to dry tons using biomass energy content) +bioused_out(bioclass,r,t)$[tmodel_new(t)$tfuel(t)] = + [BIOUSED.l(bioclass,r,t) / bio_energy_content ]; + +bioused_usda(bioclass,usda_region,t)$[tmodel_new(t)$tfuel(t)] = + [sum{r$r_usda(r,usda_region), bioused_out(bioclass,r,t) } ]; * 1e9 converts from MMBtu to Quads -repgasquant(cendiv,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 3)$tmodel_new(t)] = - sum{(gb,h), GASUSED.l(cendiv,gb,h,t) * hours(h) } * gas_scale/ 1e9 ; +repgasquant(cendiv,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 3)$tmodel_new(t)$tfuel(t)] = + [sum{(gb,h), GASUSED.l(cendiv,gb,h,t) * hours(h) } * gas_scale/ 1e9 ]; -repgasquant(cendiv,t)$[(Sw_GasCurve = 1 or Sw_GasCurve = 2)$tmodel_new(t)] = - ( sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], +repgasquant(cendiv,t)$[(Sw_GasCurve = 1 or Sw_GasCurve = 2 or Sw_FINITO_Link = 1)$tmodel_new(t)] = + [( sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t)} + sum{(v,r,h)$[valcap("dac_gas",v,r,t)$r_cendiv(r,cendiv)], hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE.l("DAC","dac_gas",v,r,h,t) }$Sw_DAC_Gas + sum{(p,i,v,r,h)$[r_cendiv(r,cendiv)$valcap(i,v,r,t)$smr(i)], hours(h) * smr_methane_rate * PRODUCE.l(p,i,v,r,h,t) }$Sw_H2 - ) / 1e9 ; + ) / 1e9 ]; -repgasquant_irt(i,r,t)$tmodel_new(t) = - ( sum{(v,h)$[valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], +repgasquant_irt(i,r,t)$[tmodel_new(t)] = + [( sum{(v,h)$[valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } + sum{(v,h)$[valcap("dac_gas",v,r,t)], hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE.l("DAC","dac_gas",v,r,h,t) }$Sw_DAC_Gas + sum{(p,v,h)$[valcap(i,v,r,t)$smr(i)], hours(h) * smr_methane_rate * PRODUCE.l(p,i,v,r,h,t) }$Sw_H2 - ) / 1e9 ; + ) / 1e9 ]; -repgasquant_nat(t)$tmodel_new(t) = sum{cendiv, repgasquant(cendiv,t) } ; +repgasquant_nat(t)$[tmodel_new(t)] = + [sum{cendiv, repgasquant(cendiv,t) } ]; *for reported gasprice (not that used to compute system costs) *scale back to $ / mmbtu -repgasprice(cendiv,t)$[(Sw_GasCurve = 0)$tmodel_new(t)$repgasquant(cendiv,t)] = - smax{gb$[sum{h, GASUSED.l(cendiv,gb,h,t) }], gasprice(cendiv,gb,t) } / gas_scale ; +repgasprice(cendiv,t)$[(Sw_GasCurve = 0)$tmodel_new(t)$repgasquant(cendiv,t)$tfuel(t)]= + [smax{gb$[sum{h, GASUSED.l(cendiv,gb,h,t) }], gasprice(cendiv,gb,t) } / gas_scale ]; -repgasprice(cendiv,t)$[(Sw_GasCurve = 2)$tmodel_new(t)$repgasquant(cendiv,t)] = - sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], +repgasprice(cendiv,t)$[(Sw_GasCurve = 2)$tmodel_new(t)$repgasquant(cendiv,t)$tfuel(t)] = + [sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], hours(h)*heat_rate(i,v,r,t)*fuel_price(i,r,t)*GEN.l(i,v,r,h,t) - } / (repgasquant(cendiv,t) * 1e9) ; + } / (repgasquant(cendiv,t) * 1e9) ]; -repgasprice_r(r,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 2)$tmodel_new(t)] = sum{cendiv$r_cendiv(r,cendiv), repgasprice(cendiv,t) } ; -repgasprice_r(r,t)$[(Sw_GasCurve = 1)$tmodel_new(t)] = - ( sum{(h,cendiv), - gasmultterm(cendiv,t) * szn_adj_gas(h) * cendiv_weights(r,cendiv) * - hours(h) } / sum{h, hours(h) } +* gas price by timeslice when linked with FINITO (see similar calculation in finito_report.gms) [$2004/MMBtu] +$ifthene.finitogasprice Sw_FINITO_Link == 1 +repgasprice_finito(cendiv,h,t)$[tmodel_new(t)$(not tfuel(t))] = + ( 1/obj_scale * 1/pvf_onm(t) * deflator('2018') * + [ max{ + smax{(eus,usep)$[map_roe_ec_usep("NG",eus,usep)], eq_supplydemand_fe_pool.m("NG",eus,usep,cendiv,h,t)}, + smax{(usep)$[map_roi_ec_usep("NG",usep)],eq_supplydemand_fe_roi.m("NG",usep,cendiv,h,t)}, + smax{(cfp2,cfvin2,fi,r)$[map_ec_ei("NG",fi)$valei_cf(cfp2,cfvin2,fi,r,t)$r_cendiv(r,cendiv)],eq_use_fi_cf.m(cfp2,cfvin2,fi,r,h,t)}, + smax{(tech,vin,fi,r)$[map_ec_ei("NG",fi)$valei_ind(tech,vin,fi,r,t)$r_cendiv(r,cendiv)],eq_use_fi_ind.m(tech,vin,fi,r,h,t)} + } / hours(h) + ] + ) +; +$else.finitogasprice + repgasprice_finito(cendiv,h,t)$[tmodel_new(t)$(not tfuel(t))] = 0 ; +$endif.finitogasprice + +repgasprice(cendiv,t)$[tmodel_new(t)$(not tfuel(t))] = +* average price over hours when linked with FINITO + sum{h, hours(h) * repgasprice_finito(cendiv,h,t) } +; + +repgasprice_r(r,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 2)$tmodel_new(t)$tfuel(t)] = + [sum{cendiv$r_cendiv(r,cendiv), repgasprice(cendiv,t) } ]; + +repgasprice_r(r,t)$[(Sw_GasCurve = 1)$tmodel_new(t)$tfuel(t)] = + sum{(h,cendiv), + gasmultterm(cendiv,t) * szn_adj_gas(h) * cendiv_weights(r,cendiv) * + hours(h) } / sum{h, hours(h) } - + smax((fuelbin,cendiv)$[VGASBINQ_REGIONAL.l(fuelbin,cendiv,t)$r_cendiv(r,cendiv)], gasbinp_regional(fuelbin,cendiv,t) ) + + smax((fuelbin,cendiv)$[VGASBINQ_REGIONAL.l(fuelbin,cendiv,t)$r_cendiv(r,cendiv)], gasbinp_regional(fuelbin,cendiv,t) ) - + smax(fuelbin$VGASBINQ_NATIONAL.l(fuelbin,t), gasbinp_national(fuelbin,t) ) - ) ; + + smax(fuelbin$VGASBINQ_NATIONAL.l(fuelbin,t), gasbinp_national(fuelbin,t) ) + ; -repgasprice(cendiv,t)$[(Sw_GasCurve = 1)$tmodel_new(t)$repgasquant(cendiv,t)] = - sum{(i,r)$r_cendiv(r,cendiv), repgasprice_r(r,t) * repgasquant_irt(i,r,t) } / repgasquant(cendiv,t) ; +repgasprice(cendiv,t)$[((Sw_GasCurve = 1) or (Sw_FINITO_Link = 1))$tmodel_new(t)$repgasquant(cendiv,t)] = + [sum{(i,r)$r_cendiv(r,cendiv), repgasprice_r(r,t) * repgasquant_irt(i,r,t) } / repgasquant(cendiv,t) ] ; repgasprice_nat(t)$[tmodel_new(t)$sum{cendiv, repgasquant(cendiv,t) }] = - sum{cendiv, repgasprice(cendiv,t) * repgasquant(cendiv,t) } - / sum{cendiv, repgasquant(cendiv,t) } ; + [sum{cendiv, repgasprice(cendiv,t) * repgasquant(cendiv,t) } + / sum{cendiv, repgasquant(cendiv,t) } ]; *======================================== * NATURAL GAS FUEL COSTS *======================================== gasshare_ba(r,cendiv,t)$[r_cendiv(r,cendiv)$tmodel_new(t)$repgasquant(cendiv,t)] = - sum{i$[valgen_irt(i,r,t)$gas(i)],repgasquant_irt(i,r,t) / repgasquant(cendiv,t) } ; + [sum{i$[valgen_irt(i,r,t)$gas(i)],repgasquant_irt(i,r,t) / repgasquant(cendiv,t) } ]; gasshare_techba(i,r,cendiv,t)$[r_cendiv(r,cendiv)$tmodel_new(t)$repgasquant(cendiv,t)$gas(i)] = - repgasquant_irt(i,r,t) / repgasquant(cendiv,t) ; + [repgasquant_irt(i,r,t) / repgasquant(cendiv,t) ]; -gasshare_cendiv(cendiv,t)$[sum{cendiv2,repgasquant(cendiv2,t)}] = repgasquant(cendiv,t) / sum{cendiv2,repgasquant(cendiv2,t)} ; +gasshare_cendiv(cendiv,t)$[sum{cendiv2,repgasquant(cendiv2,t)}] = + [repgasquant(cendiv,t) / sum{cendiv2,repgasquant(cendiv2,t)} ]; -gascost_cendiv(cendiv,t)$tmodel_new(t) = +* cost of natural gas - standalone ReEDS +gascost_cendiv(cendiv,t)$[tmodel_new(t)$tfuel(t)] = *cost of natural gas for Sw_GasCurve = 2 (static natural gas prices) - + sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t) - $[not bio(i)]$[not cofire(i)]$[Sw_GasCurve = 2]], - hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN.l(i,v,r,h,t) } + + [ sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t) + $[not bio(i)]$[not cofire(i)]$[Sw_GasCurve = 2]], + hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN.l(i,v,r,h,t) } *cost of natural gas for Sw_GasCurve = 0 (census division supply curves natural gas prices) - + sum{gb, sum{h,hours(h) * GASUSED.l(cendiv,gb,h,t) } * gasprice(cendiv,gb,t) - }$[Sw_GasCurve = 0] + + sum{gb, sum{h,hours(h) * GASUSED.l(cendiv,gb,h,t) } * gasprice(cendiv,gb,t) }$[Sw_GasCurve = 0] *cost of natural gas for Sw_GasCurve = 3 (national supply curve for natural gas prices with census division multipliers) - + sum{(h,gb), hours(h) * GASUSED.l(cendiv,gb,h,t) - * gasadder_cd(cendiv,t,h) + gasprice_nat_bin(gb,t) - }$[Sw_GasCurve = 3] + + sum{(h,gb), hours(h) * GASUSED.l(cendiv,gb,h,t) + * gasadder_cd(cendiv,t,h) + gasprice_nat_bin(gb,t) }$[Sw_GasCurve = 3] *cost of natural gas for Sw_GasCurve = 1 (national and census division supply curves for natural gas prices) *first - anticipated costs of gas consumption given last year's amount + (sum{(i,v,r,h)$[valgen(i,v,r,t)$gas(i)], @@ -593,15 +622,25 @@ gascost_cendiv(cendiv,t)$tmodel_new(t) = + sum{(fuelbin), gasbinp_national(fuelbin,t) * VGASBINQ_NATIONAL.l(fuelbin,t) } * gasshare_cendiv(cendiv,t) - )$[Sw_GasCurve = 1]; + )$[Sw_GasCurve = 1] + ]; + +* cost of natural gas - linked with FINITO ('not tfuel' indicates years using FINITO supply curves) +gascost_cendiv(cendiv,t)$[tmodel_new(t)$(not tfuel(t))] = +* cost = gas price multiplied by gas usage [$ = $/MMBtu * MMBtu/MWh * MW * h] + sum{(i,v,r,h)$[valgen(i,v,r,t)$gas(i)$r_cendiv(r, cendiv)], + repgasprice_finito(cendiv,h,t) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) * hours(h) + } +; + *======================================== * BIOFUEL COSTS *======================================== -bioshare_techba(i,r,t)$[(cofire(i) or bio(i))$tmodel_new(t)] = +bioshare_techba(i,r,t)$[(cofire(i) or bio(i))$tmodel_new(t)$tfuel(t)] = * biofuel-based generation of tech i in the BA (biopower + cofire) - (( sum{(v,h)$[valgen(i,v,r,t)$bio(i)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } + [ (( sum{(v,h)$[valgen(i,v,r,t)$bio(i)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } + sum{(v,h)$[cofire(i)$valgen(i,v,r,t)], bio_cofire_perc * hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } ) / * biofuel-based generation of all techs in the BA (biopower + cofire) @@ -611,6 +650,7 @@ bioshare_techba(i,r,t)$[(cofire(i) or bio(i))$tmodel_new(t)] = )$[ sum{(ii,v,h)$[valgen(ii,v,r,t)$bio(ii)], hours(h) * heat_rate(ii,v,r,t) * GEN.l(ii,v,r,h,t) } + sum{(ii,v,h)$[cofire(ii)$valgen(ii,v,r,t)], bio_cofire_perc * hours(h) * heat_rate(ii,v,r,t) * GEN.l(ii,v,r,h,t) } ] + ] ; *========================= @@ -1740,6 +1780,11 @@ error_check('z') = ( ) * account for penalty paid to deploy capacity beyond interconnection queue limits + sum{(tg,r), cap_penalty(tg) * CAP_ABOVE_LIM.l(tg,r,t) } +* account for costs from FINITO: deflate from $2018 to $2004, +* remove any FINITO scaling, and then apply ReEDS scaling +$ifthene.linked_objective Sw_FINITO_Link==1 + + cost_scale * ( Z_finito.l(t)$t_finito(t) * deflator('2018') / obj_scale ) +$endif.linked_objective } ) / z.l ; @@ -2033,6 +2078,8 @@ load_cat("h2_network",r,t)$tmodel_new(t) = ; load_cat("dac",r,t)$tmodel_new(t) = sum{i$dac(i), prod_load_ann(i,r,t) } ; +load_cat("finito",r,t)$tmodel_new(t) = sum{h, hours(h) * USE_ELE_FINITO.l(r,h,t) } ; + *======================================== * H2 NETWORK *======================================== diff --git a/e_report_dump.py b/e_report_dump.py index faa3a20a..aa3aa78a 100644 --- a/e_report_dump.py +++ b/e_report_dump.py @@ -67,6 +67,8 @@ def dfdict_to_h5( """ Write dictionary of dataframes to one .h5 file """ + print(f"Saving results to {filepath}") + ### unless a subset is specified, iterate over all keys in dict _symbol_list = dfdict.keys() if symbol_list is None else symbol_list @@ -108,6 +110,8 @@ def dfdict_to_excel( """ Write dictionary of dataframes to one .xlsx file """ + print(f"Saving results to {filepath}") + ### unless a subset is specified, iterate over all keys in dict _symbol_list = dfdict.keys() if symbol_list is None else symbol_list @@ -201,10 +205,13 @@ def postprocess_outputs(case, outputs_path=None, verbose=0): _outputs_path = os.path.join(case, 'outputs') if outputs_path is None else outputs_path ## System cost - reeds.output_calc.calc_systemcost(case).to_csv( - os.path.join(_outputs_path, 'post_systemcost_annualized.csv'), - index=False, - ) + try: + reeds.output_calc.calc_systemcost(case).to_csv( + os.path.join(_outputs_path, 'post_systemcost_annualized.csv'), + index=False, + ) + except FileNotFoundError: + print("Skipping system cost calculation, missing files") ## Reinforcement and spur-line reeds.output_calc.calc_reinforcement_spur_capacity_miles(case).to_csv( @@ -276,6 +283,18 @@ def postprocess_outputs(case, outputs_path=None, verbose=0): ) print("Finished loading outputs gdx") + ### FINITO outputs when running linked model + if int(sw.GSw_FINITO_Link): + print("Loading FINITO outputs gdx") + # TODO: this is much slower than loading the ReEDS outputs, + # need to figure out why and address + finito_outputs = gdxpds.to_dataframes( + os.path.join(outputs_path, f"finito_reeds_outputs_{os.path.basename(case)}.gdx") + ) + print("Finished loading FINITO outputs gdx") + dict_out.update(finito_outputs) + + write_dfdict( dfdict=dict_out, outputs_path=outputs_path, diff --git a/e_report_params.csv b/e_report_params.csv index 70b37020..f7792cb3 100644 --- a/e_report_params.csv +++ b/e_report_params.csv @@ -142,6 +142,7 @@ raw_op_cost(t),2004$,sum of operational costs from systemcost,,, RE_gen_price_nat(t),2004$/MWh,marginal cost of the national RE generation constraint,,, "repbioprice(r,t)",2004$/MMBtu,highest marginal bioprice of utilized bins for each region,,, "repgasprice(cendiv,t)",2004$/MMBtu,highest marginal gas price of utilized gas bins for each census division,,, +"repgasprice_finito(cendiv,allh,t)",2004$/MMBtu,highest marginal gas price of utilized gas bins by hour for each census division when using linked model,,, "repgasprice_r(r,t)",2004$/MMBtu,highest marginal gas price of utilized gas bins for each region,1,, repgasprice_nat(t),2004$/MMBtu,weighted-average national natural gas price assuming that plants pay the marginal price,,, "repgasquant(cendiv,t)",Quads,quantity of gas consumed in each census division,,, diff --git a/input_processing/copy_files.py b/input_processing/copy_files.py index 90e22a09..8c53b164 100644 --- a/input_processing/copy_files.py +++ b/input_processing/copy_files.py @@ -1614,6 +1614,11 @@ def main(reeds_path, inputs_case): source_deflator_map = get_source_deflator_map(reeds_path) + ### Rewrite the finito scalar as GAMS-readable definitions + if int(sw['GSw_FINITO_Link']) == 1: + pd.read_csv(os.path.join(os.path.dirname(inputs_case),'finito','inputs','scalars.csv')).iloc[:, :-1].to_csv(os.path.join(os.path.dirname(inputs_case),'finito','inputs','scalars.csv'),index=False) + scalar_csv_to_txt(os.path.join(os.path.dirname(inputs_case),'finito','inputs','scalars.csv')) + # Copy non-region files write_non_region_files(non_region_files, sw, inputs_case, regions_and_agglevel, source_deflator_map) diff --git a/input_processing/hourly_load.py b/input_processing/hourly_load.py index 8b2114fb..a29ee8f5 100644 --- a/input_processing/hourly_load.py +++ b/input_processing/hourly_load.py @@ -596,6 +596,44 @@ def reaggregate_to_model_regions( return regional_load_hourly +#TODO: add docstrings +def get_hourly_finito_load( + inputs_case: str, +) -> pd.DataFrame: + # load reference FINITO load + inputs_case_finito = Path(inputs_case).parent / 'finito' / Path(inputs_case).name + load_finito = pd.read_csv(inputs_case_finito / "load_finito.csv") + + # reshape to match load data. with year in index and region in columns + load_finito = load_finito.melt(id_vars='r', var_name='year', value_name='load_MWh') + load_finito.year = load_finito.year.astype(int) + load_finito = load_finito.pivot(index='year', columns='r', values='load_MWh') + load_finito.columns.name = None + + # allocate annual load to hours, assuming flat demand + # TODO: should this use h_weight_finito? + hours_per_year = 8760 + load_hourly_finito = load_finito / hours_per_year + load_hourly_finito = load_hourly_finito.astype(np.float32) + + return load_hourly_finito + +def remove_finito_load( + load_hourly: pd.DataFrame, + inputs_case: str, +) -> pd.DataFrame: + + # get FINITO reference load + load_hourly_finito = get_hourly_finito_load(inputs_case) + + # subtract FINITO reference load from ReEDS load data, + # aligning by model year (index) and region (columns) + result = load_hourly - load_hourly_finito + + #TODO: add this + # Validation check: confirm that total load - FINITO load > 0 + + return load_hourly #%% =========================================================================== ### --- MAIN FUNCTION --- @@ -697,6 +735,15 @@ def main(reeds_path, inputs_case): peakload = calculate_peak_load(regional_load_hourly, hierarchy) #%%%######################################### + # -- FINITO Load Adjustment -- # + ############################################# + + # note that this step occurs after peakload calculation so that the latter + # includes a baseline estimate of industrial load captured by FINITO + if int(sw.GSw_FINITO_Link): + regional_load_hourly = remove_finito_load(regional_load_hourly, inputs_case) + + ############################################# # -- DR Shed Load Modifications -- # ############################################# diff --git a/input_processing/hourly_repperiods.py b/input_processing/hourly_repperiods.py index 2c57bcc8..7d13b20e 100644 --- a/input_processing/hourly_repperiods.py +++ b/input_processing/hourly_repperiods.py @@ -36,6 +36,7 @@ import hourly_plots sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) import reeds +import hourly_load ## Time the operation of this script tic = datetime.datetime.now() @@ -75,11 +76,18 @@ def szn2period(szn): # -- Load Processing -- # ############################### -def get_load(inputs_case, keep_modelyear=None, keep_weatheryears=[2012]): +def get_load(inputs_case, sw, keep_modelyear=None, keep_weatheryears=[2012]): """ """ ### Subset to modeled regions load = reeds.io.read_file(os.path.join(inputs_case,'load.h5'), parse_timestamps=True) + + ### When running the linked model (GSw_FINITO_Link=1) we add reference load estimates + ### for FINITO load back in to use when selecting representative periods + if int(sw.GSw_FINITO_Link): + load_hourly_finito = hourly_load.get_hourly_finito_load(inputs_case) + load = load + load_hourly_finito + ### Subset to keep_modelyear if provided if keep_modelyear: load = load.loc[keep_modelyear].copy() @@ -578,6 +586,7 @@ def main( print("Collecting 8760 load data") load = get_load( inputs_case=inputs_case, + sw=sw, keep_modelyear=(int(sw['GSw_HourlyClusterYear']) if int(sw['GSw_HourlyClusterYear']) in modelyears else max(modelyears)), @@ -736,7 +745,7 @@ def main( and (sw['GSw_PRM_StressSeedLoadLevel'].lower() not in ['false','none']) ): ## Get load for all model and weather years - load_allyears = get_load(inputs_case, keep_weatheryears='all').loc[modelyears] + load_allyears = get_load(inputs_case, sw, keep_weatheryears='all').loc[modelyears] ## Add descriptive index load_allyears = load_allyears.merge( timestamps[['year', 'yperiod', 'h_of_period']], left_on='datetime', right_index=True) @@ -925,10 +934,14 @@ def main( description='Create the necessary 8760 and capacity factor data for hourly resolution') parser.add_argument('reeds_path', help='ReEDS-2.0 directory') parser.add_argument('inputs_case', help='ReEDS-2.0/runs/{case}/inputs_case directory') + parser.add_argument('--nolog', '-n', default=False, action='store_true', help='turn off logging for debugging') + args = parser.parse_args() reeds_path = args.reeds_path inputs_case = args.inputs_case + logging = not args.nolog + # #%% Settings for testing # reeds_path = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) @@ -938,10 +951,11 @@ def main( # interactive = True #%% Set up logger - log = reeds.log.makelog( - scriptname=__file__, - logpath=os.path.join(inputs_case,'..','gamslog.txt'), - ) + if logging: + log = reeds.log.makelog( + scriptname=__file__, + logpath=os.path.join(inputs_case,'..','gamslog.txt'), + ) print('Starting hourly_repperiods.py') #%% Inputs from switches sw = reeds.io.get_switches(inputs_case) @@ -956,6 +970,7 @@ def main( sw=sw, reeds_path=reeds_path, inputs_case=inputs_case, periodtype='rep', make_plots=1, + logging=logging ) ############################################ diff --git a/input_processing/hourly_writetimeseries.py b/input_processing/hourly_writetimeseries.py index 5fa32608..57c5f436 100644 --- a/input_processing/hourly_writetimeseries.py +++ b/input_processing/hourly_writetimeseries.py @@ -199,6 +199,9 @@ def get_yearly_demand(sw, hmap_myr, hmap_allyrs, inputs_case, periodtype='rep'): reload the raw demand and extract the demand on the modeled days for each year. """ ### Get original demand data, subset to cluster year + ### Note that this does not include FINITO demand when running with the linked model + ### (GSw_FINITO_Link=1) even though the reference estimates for that load is used when + ### identifying the representative days load_in = reeds.io.read_file( os.path.join(inputs_case,'load.h5'), parse_timestamps=True).unstack(level=0) load_in.columns = load_in.columns.rename(['r','t']) diff --git a/inputs/emission_constraints/co2_tax.csv b/inputs/emission_constraints/co2_tax.csv index d9e79551..1e08f624 100644 --- a/inputs/emission_constraints/co2_tax.csv +++ b/inputs/emission_constraints/co2_tax.csv @@ -1,42 +1,42 @@ -t,default,t050_2035,t050_2050,t100_2035,t100_2050,t150_2035,t150_2050,t200_2035,t200_2050 -2010,0,0,0,0,0,0,0,0,0 -2011,0,0,0,0,0,0,0,0,0 -2012,0,0,0,0,0,0,0,0,0 -2013,0,0,0,0,0,0,0,0,0 -2014,0,0,0,0,0,0,0,0,0 -2015,0,0,0,0,0,0,0,0,0 -2016,0,0,0,0,0,0,0,0,0 -2017,0,0,0,0,0,0,0,0,0 -2018,0,0,0,0,0,0,0,0,0 -2019,0,0,0,0,0,0,0,0,0 -2020,0,0,0,0,0,0,0,0,0 -2021,0,0,0,0,0,0,0,0,0 -2022,30,1,1,1,1,1,1,1,1 -2023,31.5,4.77,2.75,8.62,4.54,12.46,6.32,16.31,8.11 -2024,33.075,8.54,4.5,16.23,8.07,23.92,11.64,31.62,15.21 -2025,34.72875,12.31,6.25,23.85,11.61,35.38,16.96,46.92,22.32 -2026,36.4651875,16.08,8,31.46,15.14,46.85,22.29,62.23,29.43 -2027,38.28844688,19.85,9.75,39.08,18.68,58.31,27.61,77.54,36.54 -2028,40.20286922,23.62,11.5,46.69,22.21,69.77,32.93,92.85,43.64 -2029,42.21301268,27.38,13.25,54.31,25.75,81.23,38.25,108.15,50.75 -2030,44.32366331,31.15,15,61.92,29.29,92.69,43.57,123.46,57.86 -2031,46.53984648,34.92,16.75,69.54,32.82,104.15,48.89,138.77,64.96 -2032,48.8668388,38.69,18.5,77.15,36.36,115.62,54.21,154.08,72.07 -2033,51.31018074,42.46,20.25,84.77,39.89,127.08,59.54,169.38,79.18 -2034,53.87568978,46.23,22,92.38,43.43,138.54,64.86,184.69,86.29 -2035,56.56947427,50,23.75,100,46.96,150,70.18,200,93.39 -2036,59.39794798,50,25.5,100,50.5,150,75.5,200,100.5 -2037,62.36784538,50,27.25,100,54.04,150,80.82,200,107.61 -2038,65.48623765,50,29,100,57.57,150,86.14,200,114.71 -2039,68.76054953,50,30.75,100,61.11,150,91.46,200,121.82 -2040,72.19857701,50,32.5,100,64.64,150,96.79,200,128.93 -2041,75.80850586,50,34.25,100,68.18,150,102.11,200,136.04 -2042,79.59893115,50,36,100,71.71,150,107.43,200,143.14 -2043,83.57887771,50,37.75,100,75.25,150,112.75,200,150.25 -2044,87.7578216,50,39.5,100,78.79,150,118.07,200,157.36 -2045,92.14571268,50,41.25,100,82.32,150,123.39,200,164.46 -2046,96.75299831,50,43,100,85.86,150,128.71,200,171.57 -2047,101.5906482,50,44.75,100,89.39,150,134.04,200,178.68 -2048,106.6701806,50,46.5,100,92.93,150,139.36,200,185.79 -2049,112.0036897,50,48.25,100,96.46,150,144.68,200,192.89 -2050,117.6038742,50,50,100,100,150,150,200,200 \ No newline at end of file +t,default,t050_2035,t050_2050,t100_2035,t100_2050,t150_2035,t150_2050,t200_2035,t200_2050,t200alt_2050 +2010,0,0,0,0,0,0,0,0,0,0 +2011,0,0,0,0,0,0,0,0,0,0 +2012,0,0,0,0,0,0,0,0,0,0 +2013,0,0,0,0,0,0,0,0,0,0 +2014,0,0,0,0,0,0,0,0,0,0 +2015,0,0,0,0,0,0,0,0,0,0 +2016,0,0,0,0,0,0,0,0,0,0 +2017,0,0,0,0,0,0,0,0,0,0 +2018,0,0,0,0,0,0,0,0,0,0 +2019,0,0,0,0,0,0,0,0,0,0 +2020,0,0,0,0,0,0,0,0,0,0 +2021,0,0,0,0,0,0,0,0,0,0 +2022,30,1,1,1,1,1,1,1,1,0 +2023,31.5,4.77,2.75,8.62,4.54,12.46,6.32,16.31,8.11,0 +2024,33.075,8.54,4.5,16.23,8.07,23.92,11.64,31.62,15.21,0 +2025,34.72875,12.31,6.25,23.85,11.61,35.38,16.96,46.92,22.32,0 +2026,36.4651875,16.08,8,31.46,15.14,46.85,22.29,62.23,29.43,0 +2027,38.28844688,19.85,9.75,39.08,18.68,58.31,27.61,77.54,36.54,0 +2028,40.20286922,23.62,11.5,46.69,22.21,69.77,32.93,92.85,43.64,0 +2029,42.21301268,27.38,13.25,54.31,25.75,81.23,38.25,108.15,50.75,0 +2030,44.32366331,31.15,15,61.92,29.29,92.69,43.57,123.46,57.86,0 +2031,46.53984648,34.92,16.75,69.54,32.82,104.15,48.89,138.77,64.96,10 +2032,48.8668388,38.69,18.5,77.15,36.36,115.62,54.21,154.08,72.07,20 +2033,51.31018074,42.46,20.25,84.77,39.89,127.08,59.54,169.38,79.18,30 +2034,53.87568978,46.23,22,92.38,43.43,138.54,64.86,184.69,86.29,40 +2035,56.56947427,50,23.75,100,46.96,150,70.18,200,93.39,50 +2036,59.39794798,50,25.5,100,50.5,150,75.5,200,100.5,60 +2037,62.36784538,50,27.25,100,54.04,150,80.82,200,107.61,70 +2038,65.48623765,50,29,100,57.57,150,86.14,200,114.71,80 +2039,68.76054953,50,30.75,100,61.11,150,91.46,200,121.82,90 +2040,72.19857701,50,32.5,100,64.64,150,96.79,200,128.93,100 +2041,75.80850586,50,34.25,100,68.18,150,102.11,200,136.04,110 +2042,79.59893115,50,36,100,71.71,150,107.43,200,143.14,120 +2043,83.57887771,50,37.75,100,75.25,150,112.75,200,150.25,130 +2044,87.7578216,50,39.5,100,78.79,150,118.07,200,157.36,140 +2045,92.14571268,50,41.25,100,82.32,150,123.39,200,164.46,150 +2046,96.75299831,50,43,100,85.86,150,128.71,200,171.57,160 +2047,101.5906482,50,44.75,100,89.39,150,134.04,200,178.68,170 +2048,106.6701806,50,46.5,100,92.93,150,139.36,200,185.79,180 +2049,112.0036897,50,48.25,100,96.46,150,144.68,200,192.89,190 +2050,117.6038742,50,50,100,100,150,150,200,200,200 diff --git a/reeds/reedsplots.py b/reeds/reedsplots.py index da1013b1..1630ca17 100644 --- a/reeds/reedsplots.py +++ b/reeds/reedsplots.py @@ -5020,6 +5020,7 @@ def plot_seed_stressperiods( ### Recalculate peak load days since we dropped duplicates above load_allyears = hourly_repperiods.get_load( os.path.join(case, 'inputs_case'), + sw, keep_weatheryears='all').loc[years] timestamps = pd.read_csv(os.path.join(case, 'inputs_case', 'rep', 'timestamps.csv')) resource_adequacy_years = [int(y) for y in sw.resource_adequacy_years.split('_')] diff --git a/restart_runs.py b/restart_runs.py index dee9a85c..71eb72e7 100644 --- a/restart_runs.py +++ b/restart_runs.py @@ -4,6 +4,8 @@ import subprocess import argparse import pandas as pd +import re +import reeds from glob import glob from runbatch import submit_slurm_parallel_jobs from runstatus import get_run_status @@ -92,7 +94,14 @@ #%% Copy additional files if desired for f in more_copyfiles: - shutil.copy(os.path.join(reeds_path,f), os.path.join(case,f)) + if f.lower().startswith('finito'): + # if file starts with 'finito' append the finito directory + sw = reeds.io.get_switches(case) + if int(sw.GSw_FINITO_Link): + f_copy = re.sub("^finito/", "" , f, flags=re.IGNORECASE) + shutil.copy(os.path.join(sw.finito_dir,f_copy), os.path.join(case,f)) + else: + shutil.copy(os.path.join(reeds_path,f), os.path.join(case,f)) #%% Make a backup copy of the original bash and sbatch scripts callfile = os.path.join(case,f'call_{casename}.sh') @@ -165,7 +174,8 @@ # It is a single case or we are not on HPC if copy_srun_template: writelines_srun_case = writelines_srun.copy() - writelines_srun_case.append(f"\n#SBATCH --job-name={casename}\n") + writelines_srun_case.append(f"#SBATCH --job-name={casename}") + writelines_srun_case.append(f"#SBATCH --output={os.path.join(case, 'slurm-%j.out')}\n") writelines_srun_case.append(f"sh {callfile}") with open(sbatchfile, 'w') as f: for line in writelines_srun_case: diff --git a/runbatch.py b/runbatch.py index 33b32979..106d2735 100644 --- a/runbatch.py +++ b/runbatch.py @@ -2,6 +2,10 @@ ### --- IMPORTS --- ### =========================================================================== +#%% =========================================================================== +### --- IMPORTS --- +### =========================================================================== + import os import git import queue @@ -11,6 +15,8 @@ import csv import importlib import numpy as np +import gdxpds +from gdxpds.gdx import GdxFile,GdxSymbol,GamsDataType import pandas as pd import subprocess import re @@ -18,6 +24,7 @@ import argparse from pathlib import Path import reeds +import sys # Assert core programs are accessible CORE_PROGRAMS = ["gams"] @@ -99,12 +106,23 @@ def create_case_lists(df_cases:pd.DataFrame, BatchName:str, single:str=''): continue # Add switch settings to list of options passed to GAMS shcom = f' --case={BatchName}_{case}' - for i,v in df_cases[case].items(): + # (linked) Combine the cases files for the linked model + if int(df_cases[case].loc['GSw_FINITO_Link']) == 1 : + # add the FINITO switches to the caseSwitches + df_cases_combine=linked_cases(df_cases,case) + for i,v in df_cases_combine.items(): + #exclude certain switches that don't need to be passed to GAMS + if i not in ['file_replacements','keep_run_terminal']: + shcom = shcom + ' --' + i + '=' + v + caseList.append(shcom) + caseSwitches.append(df_cases_combine.to_dict()) + else: + for i,v in df_cases[case].items(): #exclude certain switches that don't need to be passed to GAMS - if i not in ['file_replacements','keep_run_terminal']: - shcom += f' --{i}={v}' - caseList.append(shcom) - caseSwitches.append(df_cases[case].to_dict()) + if i not in ['file_replacements','keep_run_terminal']: + shcom += f' --{i}={v}' + caseList.append(shcom) + caseSwitches.append(df_cases[case].to_dict()) return caseSwitches, casenames, caseList @@ -554,6 +572,7 @@ def solvestring_sequential( 'GSw_StateCO2ImportLevel', 'GSw_StartMarkets', 'GSw_ValStr', + 'GSw_FINITO_Link', 'solver', 'debug', 'startyear', @@ -807,6 +826,71 @@ def setup_window( OPATH.writelines( "python valuestreams.py" + '\n') +# (linked) update 'df_cases' to include FINITO switches +def linked_cases(df_cases,case): + """ + Updates the cases dataframe to include FINITO-specific switches. + When a switch is duplicated in FINITO and ReEDS, then we default to + the ReEDS value. + + For the FINITO switches, the combined cases file defaults to the case-specific + 'Default Value' in cases_linked.csv first, before using the universal FINITO + 'Default Value' in cases.csv for any un-assigned switches. + """ + # check for valid finito_dir + if not os.path.isdir(df_cases[case]['finito_dir']): + raise ValueError(f"finito_dir = {df_cases[case]['finito_dir']} is not a valid path.") + + # define path to and read the FINITO check_inputs function + finito_check_inputs_path = os.path.join(df_cases[case]['finito_dir'], 'input_processing', 'processing') + sys.path.append(finito_check_inputs_path) + from check_inputs import check_inputs + + ### load the default values for all FINITO switches from ~\FINITO\cases.csv + df_cases_finito = pd.read_csv(os.path.join(df_cases[case]['finito_dir'],'cases.csv'), dtype=object, index_col=0) + df_cases_finito = df_cases_finito[['Choices', 'Default Value']] + + ### load the scenario-specific switches from ~\FINITO\cases_linked.csv + cases_linked_path = os.path.join(df_cases[case]['finito_dir'],f"cases_{df_cases[case]['finito_cases_file']}.csv") + df_cases_suf_finito = pd.read_csv(cases_linked_path, dtype=object, index_col=0) + ## check that case names are unique in cases_linked.csv + # grab the scenario names **exactly** as they appear in the csv file + with open(cases_linked_path, 'r', newline='') as csvfile: + reader = csv.reader(csvfile) + header = next(reader) + # find the duplicate column names and raise an error if any are found + duplicate_columns = set([x for x in header if header.count(x) > 1]) + if duplicate_columns: + raise ValueError(f"The FINITO cases_{df_cases[case]['finito_cases_file']}.csv has the following duplicate column names: {duplicate_columns}") + ### identify the FINITO case + if df_cases[case]['finito_case'] == 'same': + finito_case=case + else: + finito_case=df_cases[case]['finito_case'] + # ensures **exact** match of names between the ReEDS cases_{}.csv and the FINITO cases_linked.csv + if finito_case not in (df_cases_suf_finito.columns): + raise ValueError(f"The 'finito_case' input '{finito_case}' in the ReEDS cases file does not exist in FINITO's cases_{df_cases[case]['finito_cases_file']}.csv.") + + ### first use 'Default Value' from the FINITO cases_linked.csv to fill missing switches + if 'Default Value' in df_cases_suf_finito.columns: + df_cases_suf_finito[finito_case] = df_cases_suf_finito[finito_case].fillna(df_cases_suf_finito['Default Value']) + ### then, use 'Default Value' from the FINITO cases.csv to fill un-assigned switches + df_cases_suf_finito.drop(['Choices','Default Value'], axis='columns',inplace=True, errors='ignore') + df_cases_finito = df_cases_finito.join(df_cases_suf_finito, how='outer') + df_cases_finito[finito_case] = df_cases_finito[finito_case].fillna(df_cases_finito['Default Value']) + + #### create new dataframe for the combined ReEDS and FINITO switches + df_cases_combine = pd.concat([df_cases[case],df_cases_finito[finito_case]]) + ### drop duplicated switches, defaulting to reeds + df_cases_combine = df_cases_combine[~df_cases_combine.index.duplicated(keep='first')] + + #%% Check for incompatibility of FINITO switches + model_sectors = df_cases_finito['Default Value']['focus_sectors'].split('.') + check_inputs(case = case, df_case = df_cases_combine, model_sectors=model_sectors) + + return df_cases_combine + + #%% =========================================================================== ### --- PROCEDURE --- ### =========================================================================== @@ -1260,6 +1344,7 @@ def write_batch_script( yearset_augur = os.path.join('inputs_case','modeledyears.csv') toLogGamsString = ' logOption=4 logFile=gamslog.txt appendLog=1 ' + ## Copy code folders for dirname in ['reeds', 'ReEDS_Augur', 'input_processing', 'reeds2pras']: shutil.copytree( @@ -1271,6 +1356,66 @@ def write_batch_script( #make the augur_data folder os.makedirs(os.path.join(casedir,'ReEDS_Augur','augur_data'), exist_ok=True) os.makedirs(os.path.join(casedir,'ReEDS_Augur','PRAS'), exist_ok=True) + + ## (linked) Copy FINITO code folders and inputs into [casedir]/finito + if int(caseSwitches['GSw_FINITO_Link'])==1: + ## Redirecting the FINITO case-specific directories + # ... finito directory within the case directory + casedir_finito = os.path.join(casedir,'finito') + # ... define the inputs case directory for FINITO + inputs_case_finito = os.path.join(casedir,'finito','inputs_case') + + ## copy directories + os.makedirs(casedir_finito, exist_ok=True) + shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'inputs'),os.path.join(casedir,'finito', 'inputs')) + shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'input_processing'),os.path.join(casedir,'finito', 'input_processing')) + shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'model'),os.path.join(casedir,'finito', 'model')) + shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'visualization'),os.path.join(casedir, 'finito', 'visualization')) + shutil.copytree(os.path.join(caseSwitches['finito_dir'],'inputs'), os.path.join(casedir,'inputs'), dirs_exist_ok=True) + ## copy over the FINITO cases files + shutil.copy2(os.path.join(caseSwitches['finito_dir'], 'cases.csv'), os.path.join(casedir, 'finito')) + shutil.copy2(os.path.join(caseSwitches['finito_dir'], f"cases_{caseSwitches['finito_cases_file']}.csv"), os.path.join(casedir, 'finito')) + + + #%% (GSw_Trade_PriceResponse > 0) If doing a price-responsive trade run, retrieve the reference exports/imports prices + if int(caseSwitches['GSw_Trade_PriceResponse']) > 0: + initialize_price_response_path = os.path.join(casedir_finito, 'input_processing', 'processing', 'initialize_price_response.py') + # Collect all arguments for initialize_price_response.py + initialize_price_response_args = f" -c {casedir} -b {BatchName} -cr {caseSwitches['GSw_Trade_PriceResponse_RefScen']} -l {caseSwitches['GSw_FINITO_Link']}" + # Call initialize_price_response.py file before starting the runs + os.system('python ' + initialize_price_response_path + initialize_price_response_args) + + #%% Filter and copy all input files for each scenario + # Call FINITO copy_files.py file before starting the runs + copy_files_run = os.system( + 'python ' + os.path.join(caseSwitches['finito_dir'], 'input_processing', 'processing', 'copy_files.py') + + f" -c {casedir_finito} -d {inputs_case_finito} --link" + ) + + # Print an error message if copy_files.py encounters any issue + if copy_files_run != 0: + print("\n❌ ERROR: copy_files.py encountered an issue and did not complete successfully.") + print("Please check the console output above for details.") + print("The issue could be due to regionality, focus sector filtering, or file reading errors.") + os._exit(1) + + ## Populate sets for each linked run + + # Pass AEO base year as an argument + # NOTE: make this dynamic before merging in? + aeo_year = caseSwitches['aeo_year'] + # Collect all arguments for autopop_set.py + focus_sectors = caseSwitches['focus_sectors'].replace('.', ' ') + autopop_args = f" -c {casedir_finito} -d {inputs_case_finito} -a {aeo_year} -s {focus_sectors} -f {caseSwitches['GSw_FixedCostSupply']} -rwf {caseSwitches['GSw_Trade_Partners']} -e {caseSwitches['GSw_ROE_EndUses']}" + # Call autopop_set.py file before starting the runs + autopop_path = os.path.join(caseSwitches['finito_dir'], 'input_processing','processing','autopop_set.py') + os.system('python ' + autopop_path + autopop_args) + ## Call read_mecs_heat.py to generate heat/nonheat/feedstock ratios for FINITO Rest of Industry (ROI) + mecs_sectors = focus_sectors + read_mecs_path = os.path.join(caseSwitches['finito_dir'], 'input_processing', 'processing', 'mecs', 'read_mecs_heat.py') + # Collect all arguments for read_mecs_heat.py + read_mecs_args = f' -s {mecs_sectors} -d {inputs_case_finito}' + os.system('python ' + read_mecs_path + read_mecs_args) ###### Replace files according to 'file_replacements' in cases. Ignore quotes in input text. # << is used to separate the file that is to be replaced from the file that is used @@ -1309,6 +1454,7 @@ def write_batch_script( OPATH.writelines("module load conda \n") OPATH.writelines("module load gams \n") + OPATH.writelines("conda deactivate \n") OPATH.writelines("conda activate reeds2 \n") OPATH.writelines('export R_LIBS_USER="$HOME/rlib" \n\n\n') @@ -1412,7 +1558,8 @@ def write_batch_script( + ' gdxcompress=1' + toLogGamsString + f"--fname={batch_case}" - + f" --GSw_calc_powfrac={caseSwitches['GSw_calc_powfrac']} \n" + + f" --GSw_calc_powfrac={caseSwitches['GSw_calc_powfrac']}" + + f" --first_year_finito={caseSwitches['first_year_finito']} \n" ) OPATH.writelines(writescripterrorcheck("e_report.gms")) if not LINUXORMAC: @@ -1425,6 +1572,32 @@ def write_batch_script( + f" {os.path.join(reeds_path,'postprocessing','diagnose','diagnose_process.py')}" + f" --casepath {casedir} \n\n" ) + OPATH.writelines(f'python e_report_dump.py {casedir} -c\n\n') + # (linked) calls finito reporting (from [casedir]/finito/visualization) + if int(caseSwitches['GSw_FINITO_Link'])==1: + OPATH.writelines( + 'gams ' + + f"{os.path.join(casedir,'finito', 'model', 'finito_report.gms')}" + + f" o={os.path.join('lstfiles',f'finito_report_{batch_case}.lst')}" + + (' license=gamslice.txt' if hpc else '') + + (' r=$r' if LINUXORMAC else ' r=!r!') + + ' gdxcompress=1' + + toLogGamsString + + f"--case={batch_case}" + + f" --casedir={casedir}" + + f" --GSw_FINITO_Link={caseSwitches['GSw_FINITO_Link']}" + + f" --GSw_RetailAdder={caseSwitches['GSw_RetailAdder']}" + + f" --finito_inputs_dir={os.path.join('finito','inputs')}" + + f" --inputs_case_finito_dir={os.path.join('finito','inputs_case')}" + + f" --linked_report_dir={caseSwitches['linked_report_dir']} \n" + ) + # (linked) calls finito postprocessing (from [casedir]/finito/visualization) + OPATH.writelines( + f"python {os.path.join(casedir, 'finito', 'visualization', 'postprocessing.py')}" + + " -b 0" + + f" -l {caseSwitches['GSw_FINITO_Link']}" + + f' -c {batch_case} \n\n' + ) ### Run the retail rate module OPATH.writelines( @@ -1727,7 +1900,7 @@ def launch_single_case_run( shellscript = subprocess.Popen( [os.path.join(casedir, 'call_' + batch_case + ext)], shell=True) # Wait for it to finish before killing the thread - shellscript.wait() + #shellscript.wait() else: if int(caseSwitches['keep_run_terminal']) == 1: terminal_keep_flag = ' /k ' @@ -1855,7 +2028,7 @@ def main( help="Check inputs but don't start runs") args = parser.parse_args() - + main( BatchName=args.BatchName, cases_suffix=args.cases_suffix, single=args.single, simult_runs=args.simult_runs, forcelocal=args.forcelocal, skip_checks=args.skip_checks, From 5fd48485826754d2e972fad27e56b32eb2878849 Mon Sep 17 00:00:00 2001 From: mmowers Date: Thu, 23 Apr 2026 17:38:10 -0400 Subject: [PATCH 03/15] Fix windows bug in outputs for linked ReEDS-FINITO runs --- runbatch.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/runbatch.py b/runbatch.py index 106d2735..104b690e 100644 --- a/runbatch.py +++ b/runbatch.py @@ -1562,7 +1562,7 @@ def write_batch_script( + f" --first_year_finito={caseSwitches['first_year_finito']} \n" ) OPATH.writelines(writescripterrorcheck("e_report.gms")) - if not LINUXORMAC: + if not LINUXORMAC and int(caseSwitches['GSw_FINITO_Link']) != 1: OPATH.writelines("endlocal\n") OPATH.writelines(f'python {logger}\n') OPATH.writelines(f'python e_report_dump.py {casedir} -c\n\n') @@ -1598,7 +1598,8 @@ def write_batch_script( + f" -l {caseSwitches['GSw_FINITO_Link']}" + f' -c {batch_case} \n\n' ) - + if not LINUXORMAC: + OPATH.writelines("endlocal\n") ### Run the retail rate module OPATH.writelines( "python" From 99147d174e00466f0e0b0ef6925ac2710f6355cf Mon Sep 17 00:00:00 2001 From: bsergi Date: Wed, 13 May 2026 12:46:33 -0400 Subject: [PATCH 04/15] Update to autopop_set.py call --- runbatch.py | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/runbatch.py b/runbatch.py index 067880d7..b0c937b6 100644 --- a/runbatch.py +++ b/runbatch.py @@ -1399,7 +1399,7 @@ def write_batch_script( shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'input_processing'),os.path.join(casedir,'finito', 'input_processing')) shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'model'),os.path.join(casedir,'finito', 'model')) shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'visualization'),os.path.join(casedir, 'finito', 'visualization')) - shutil.copytree(os.path.join(caseSwitches['finito_dir'],'inputs'), os.path.join(casedir,'inputs'), dirs_exist_ok=True) + ## copy over the FINITO cases files shutil.copy2(os.path.join(caseSwitches['finito_dir'], 'cases.csv'), os.path.join(casedir, 'finito')) shutil.copy2(os.path.join(caseSwitches['finito_dir'], f"cases_{caseSwitches['finito_cases_file']}.csv"), os.path.join(casedir, 'finito')) @@ -1427,19 +1427,15 @@ def write_batch_script( print("The issue could be due to regionality, focus sector filtering, or file reading errors.") os._exit(1) - ## Populate sets for each linked run - - # Pass AEO base year as an argument - # NOTE: make this dynamic before merging in? - aeo_year = caseSwitches['aeo_year'] - # Collect all arguments for autopop_set.py - focus_sectors = caseSwitches['focus_sectors'].replace('.', ' ') - autopop_args = f" -c {casedir_finito} -d {inputs_case_finito} -a {aeo_year} -s {focus_sectors} -f {caseSwitches['GSw_FixedCostSupply']} -rwf {caseSwitches['GSw_Trade_Partners']} -e {caseSwitches['GSw_ROE_EndUses']}" - # Call autopop_set.py file before starting the runs - autopop_path = os.path.join(caseSwitches['finito_dir'], 'input_processing','processing','autopop_set.py') - os.system('python ' + autopop_path + autopop_args) + ## Populate sets for each linked run using autopop_set.py + #autopop_args = f" -c {casedir_finito} -d {inputs_case_finito} -a {aeo_year} -s {focus_sectors} -f {caseSwitches['GSw_FixedCostSupply']} -rwf {caseSwitches['GSw_Trade_Partners']} -e {caseSwitches['GSw_ROE_EndUses']}" + os.system( + 'python ' + os.path.join(caseSwitches['finito_dir'], 'input_processing', 'processing', 'autopop_set.py') + + f" -c {casedir_finito} -d {inputs_case_finito} --link" + ) + ## Call read_mecs_heat.py to generate heat/nonheat/feedstock ratios for FINITO Rest of Industry (ROI) - mecs_sectors = focus_sectors + mecs_sectors = caseSwitches['focus_sectors'].replace('.', ' ') read_mecs_path = os.path.join(caseSwitches['finito_dir'], 'input_processing', 'processing', 'mecs', 'read_mecs_heat.py') # Collect all arguments for read_mecs_heat.py read_mecs_args = f' -s {mecs_sectors} -d {inputs_case_finito}' From 6c98dde896f345d659a0528063360b0b1ee02b0d Mon Sep 17 00:00:00 2001 From: mmowers Date: Sat, 30 May 2026 04:43:29 +0000 Subject: [PATCH 05/15] Update region switches in cases_finito.csv --- cases_finito.csv | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/cases_finito.csv b/cases_finito.csv index 25fb39db..086d15d9 100644 --- a/cases_finito.csv +++ b/cases_finito.csv @@ -8,6 +8,7 @@ coalscen,AEO_2025_reference,,,,,,,,,,,,,,,,,,,,, ngscen,AEO_2025_reference,,,,,AEO_2025_HOG,,,AEO_2025_HOG,,,,,,,,,,,,, uraniumscen,AEO_2025_reference,,,,,,,,,,,,,,,,,,,,, GSw_Region,country/USA,,st/CO,,,,,,,,,,,,,,,,,,, +GSw_ZoneSet,z132,,,,,,,,,,,,,,,,,,,,,z54 GSw_GasCurve,1,,2,,,,,,,,,,,,,,,,,,, GSw_gopt,2,,,,,,,,,,,,,,,,,,,,, GSw_NewValCapShrink,1,,,,,,,,,,,,,,,,,,,,, @@ -31,7 +32,6 @@ GSw_H2Combustionupgrade,1,,,,,,,,,,,,,,,,,,,,, GSw_H2CombinedCycle,1,,,,,,,,,,,,,,,,,,,,, GSw_DAC,1,,,,,,,,,,,,,,,,,,,,, GSw_BECCS,1,,,,,,,,,,,,,,,,,,,,, -GSw_RegionResolution,ba,,,,,,,,,,,,,,,,,,,,,aggreg GSw_TCPhaseout,1,,,,,,,,,,,,,,,,,,,,, GSw_CO2_Detail,1,,,,,,,,,,,,,,,,,,,,, GSw_CO2_LimitStorageSites,1,,,,,,,,,,,,,,,,,,,,, @@ -43,4 +43,3 @@ GSw_PRM_StressIterateMax,0,,,,,,,,,,,,,,,,,,,,, GSw_HourlyNumClusters,9,,,,,,,,,,,,,,,,,,,,, GSw_HourlyChunkLengthRep,4,,,,,,,,,,,,,,,,,,,,, GSw_HourlyChunkLengthStress,4,,,,,,,,,,,,,,,,,,,,, -GSw_HierarchyFile,default,,,,,,,,,,,,,,,,,,,,,agg54 From 046b21de2d271495d705f32264b2739ed73d965f Mon Sep 17 00:00:00 2001 From: mmowers Date: Tue, 2 Jun 2026 19:58:47 -0400 Subject: [PATCH 06/15] Fix bug in report.gms for linked FINITO runs --- reeds/core/terminus/report.gms | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/reeds/core/terminus/report.gms b/reeds/core/terminus/report.gms index e0026dd2..882a2621 100644 --- a/reeds/core/terminus/report.gms +++ b/reeds/core/terminus/report.gms @@ -547,8 +547,8 @@ $ifthene.finitogasprice Sw_FINITO_Link == 1 repgasprice_finito(cendiv,h,t)$[tmodel_new(t)$(not tfuel(t))] = ( 1/obj_scale * 1/pvf_onm(t) * deflator('2018') * [ max{ - smax{(eus,usep)$[map_roe_ec_usep("NG",eus,usep)], eq_supplydemand_fe_pool.m("NG",eus,usep,cendiv,h,t)}, - smax{(usep)$[map_roi_ec_usep("NG",usep)],eq_supplydemand_fe_roi.m("NG",usep,cendiv,h,t)}, + smax{(eus,usep,r)$[map_roe_ec_usep("NG",eus,usep)$r_cendiv(r,cendiv)], eq_supplydemand_fe_pool.m("NG",eus,usep,r,h,t)}, + smax{(usep,r)$[map_roi_ec_usep("NG",usep)$r_cendiv(r,cendiv)],eq_supplydemand_fe_roi.m("NG",usep,r,h,t)}, smax{(cfp2,cfvin2,fi,r)$[map_ec_ei("NG",fi)$valei_cf(cfp2,cfvin2,fi,r,t)$r_cendiv(r,cendiv)],eq_use_fi_cf.m(cfp2,cfvin2,fi,r,h,t)}, smax{(tech,vin,fi,r)$[map_ec_ei("NG",fi)$valei_ind(tech,vin,fi,r,t)$r_cendiv(r,cendiv)],eq_use_fi_ind.m(tech,vin,fi,r,h,t)} } / hours(h) From 8cecf3cdb1064a2eb58745553692ee731ba79640 Mon Sep 17 00:00:00 2001 From: mmowers Date: Tue, 2 Jun 2026 20:05:34 -0400 Subject: [PATCH 07/15] Add comment to report.gms --- reeds/core/terminus/report.gms | 1 + 1 file changed, 1 insertion(+) diff --git a/reeds/core/terminus/report.gms b/reeds/core/terminus/report.gms index 882a2621..8c0d8441 100644 --- a/reeds/core/terminus/report.gms +++ b/reeds/core/terminus/report.gms @@ -543,6 +543,7 @@ repgasprice(cendiv,t)$[(Sw_GasCurve = 2)$tmodel_new(t)$repgasquant(cendiv,t)$tfu * gas price by timeslice when linked with FINITO (see similar calculation in finito_report.gms) [$2004/MMBtu] +* We maybe should be taking weighted averages instead of max across regions and categories (pool vs roi etc.) $ifthene.finitogasprice Sw_FINITO_Link == 1 repgasprice_finito(cendiv,h,t)$[tmodel_new(t)$(not tfuel(t))] = ( 1/obj_scale * 1/pvf_onm(t) * deflator('2018') * From b0720c9b312090383c38149788db23794a5982ef Mon Sep 17 00:00:00 2001 From: mmowers Date: Wed, 3 Jun 2026 11:21:58 -0400 Subject: [PATCH 08/15] Convert USE_ELE_FINITO to busbar load when using in ReEDS --- reeds/core/setup/c_model.gms | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/reeds/core/setup/c_model.gms b/reeds/core/setup/c_model.gms index 3f95c908..a81e5ee4 100644 --- a/reeds/core/setup/c_model.gms +++ b/reeds/core/setup/c_model.gms @@ -399,10 +399,10 @@ eq_loadcon(r,h,t)$tmodel(t).. * [plus] load for industrial and converted fuel facilities (FINITO), * including the PRM for stress periods -* TODO: should this include distloss? +* USE_ELE_FINITO is end-use, so divide by (1-distloss) to convert it to busbar * [MWh/hr = MW] $ifthene.linked_load Sw_FINITO_Link==1 - + USE_ELE_FINITO(r,h,t)$[t_finito(t)] * (1 + prm(r,t)$h_stress(h)) + + (USE_ELE_FINITO(r,h,t) / (1.0 - distloss))$[t_finito(t)] * (1 + prm(r,t)$h_stress(h)) $endif.linked_load ; From bf47a9ba02bf72e004ca3c06527f9a486171945d Mon Sep 17 00:00:00 2001 From: bsergi Date: Mon, 8 Jun 2026 11:04:45 -0400 Subject: [PATCH 09/15] error handling for missing finito gdx file --- reeds/core/terminus/report_dump.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/reeds/core/terminus/report_dump.py b/reeds/core/terminus/report_dump.py index 4b3e5c81..2a18ed4b 100644 --- a/reeds/core/terminus/report_dump.py +++ b/reeds/core/terminus/report_dump.py @@ -282,13 +282,13 @@ def postprocess_outputs(case, outputs_path=None, verbose=0): ### FINITO outputs when running linked model if int(sw.GSw_FINITO_Link): print("Loading FINITO outputs gdx") - # TODO: this is much slower than loading the ReEDS outputs, - # need to figure out why and address - finito_outputs = gdxpds.to_dataframes( - os.path.join(outputs_path, f"finito_reeds_outputs_{os.path.basename(case)}.gdx") - ) - print("Finished loading FINITO outputs gdx") - dict_out.update(finito_outputs) + finito_gdx = os.path.join(outputs_path, f"finito_reeds_outputs_{os.path.basename(case)}.gdx") + if os.path.isfile(finito_gdx): + finito_outputs = gdxpds.to_dataframes(finito_gdx) + dict_out.update(finito_outputs) + print("Finished loading FINITO outputs gdx") + else: + print(f"FINITO outputs in {finito_gdx} could not be found, skipping.") write_dfdict( From eaaaa7e9f5513a499713a008d2862215a1fe9643 Mon Sep 17 00:00:00 2001 From: bsergi Date: Mon, 8 Jun 2026 11:07:51 -0400 Subject: [PATCH 10/15] error handling for finito gdx --- reeds/core/terminus/report_dump.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/reeds/core/terminus/report_dump.py b/reeds/core/terminus/report_dump.py index 2a18ed4b..617433c9 100644 --- a/reeds/core/terminus/report_dump.py +++ b/reeds/core/terminus/report_dump.py @@ -283,12 +283,13 @@ def postprocess_outputs(case, outputs_path=None, verbose=0): if int(sw.GSw_FINITO_Link): print("Loading FINITO outputs gdx") finito_gdx = os.path.join(outputs_path, f"finito_reeds_outputs_{os.path.basename(case)}.gdx") - if os.path.isfile(finito_gdx): + try: finito_outputs = gdxpds.to_dataframes(finito_gdx) dict_out.update(finito_outputs) print("Finished loading FINITO outputs gdx") - else: - print(f"FINITO outputs in {finito_gdx} could not be found, skipping.") + except Exception as err: + print(err) + print(f"Error loading FINITO outputs in {finito_gdx}, skipping.") write_dfdict( From 877b591623ce72443d45cb2c063cae17c553af02 Mon Sep 17 00:00:00 2001 From: mmowers Date: Tue, 9 Jun 2026 01:19:38 -0400 Subject: [PATCH 11/15] Fix bug in standalone ReEDS outputs by putting if statement around FINITO output --- reeds/core/terminus/report.gms | 2 ++ 1 file changed, 2 insertions(+) diff --git a/reeds/core/terminus/report.gms b/reeds/core/terminus/report.gms index 8c0d8441..cdc0b943 100644 --- a/reeds/core/terminus/report.gms +++ b/reeds/core/terminus/report.gms @@ -2079,7 +2079,9 @@ load_cat("h2_network",r,t)$tmodel_new(t) = ; load_cat("dac",r,t)$tmodel_new(t) = sum{i$dac(i), prod_load_ann(i,r,t) } ; +$ifthene.finitoloadcat Sw_FINITO_Link == 1 load_cat("finito",r,t)$tmodel_new(t) = sum{h, hours(h) * USE_ELE_FINITO.l(r,h,t) } ; +$endif.finitoloadcat *======================================== * H2 NETWORK From 86c4093e8baebc020c9b493724730ba1702e0285 Mon Sep 17 00:00:00 2001 From: mmowers Date: Tue, 9 Jun 2026 12:09:19 -0400 Subject: [PATCH 12/15] Update to z54 by default in cases_finito.csv --- cases_finito.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cases_finito.csv b/cases_finito.csv index 086d15d9..ca5459a3 100644 --- a/cases_finito.csv +++ b/cases_finito.csv @@ -8,7 +8,7 @@ coalscen,AEO_2025_reference,,,,,,,,,,,,,,,,,,,,, ngscen,AEO_2025_reference,,,,,AEO_2025_HOG,,,AEO_2025_HOG,,,,,,,,,,,,, uraniumscen,AEO_2025_reference,,,,,,,,,,,,,,,,,,,,, GSw_Region,country/USA,,st/CO,,,,,,,,,,,,,,,,,,, -GSw_ZoneSet,z132,,,,,,,,,,,,,,,,,,,,,z54 +GSw_ZoneSet,z54,,,,,,,,,,,,,,,,,,,,,z54 GSw_GasCurve,1,,2,,,,,,,,,,,,,,,,,,, GSw_gopt,2,,,,,,,,,,,,,,,,,,,,, GSw_NewValCapShrink,1,,,,,,,,,,,,,,,,,,,,, From 74443a6da0ee4e2bc7dd5be63f8fdf949502fb61 Mon Sep 17 00:00:00 2001 From: mmowers Date: Sat, 13 Jun 2026 00:24:20 -0400 Subject: [PATCH 13/15] Remove round statement to align with main --- reeds/core/setup/b_inputs.gms | 1 - 1 file changed, 1 deletion(-) diff --git a/reeds/core/setup/b_inputs.gms b/reeds/core/setup/b_inputs.gms index c8a410a7..c45bda48 100644 --- a/reeds/core/setup/b_inputs.gms +++ b/reeds/core/setup/b_inputs.gms @@ -1989,7 +1989,6 @@ noncumulative_prescriptions(pcat,r,t)$tmodel_new(t) ], prescribednonrsc(tt,pcat,r,"value") + prescribedrsc(tt,pcat,r,"value") } ; -noncumulative_prescriptions(pcat,r,t)$noncumulative_prescriptions(pcat,r,t) = round(noncumulative_prescriptions(pcat,r,t),6); parameter noncumulative_prescriptions_energy(pcat,r,t) "--MWh-- prescribed energy capacity that comes online in a given year" ; noncumulative_prescriptions_energy(pcat,r,t)$tmodel_new(t) From 968806411a267b7abfb362b2d8c709410c548e53 Mon Sep 17 00:00:00 2001 From: mmowers Date: Thu, 18 Jun 2026 21:49:35 -0400 Subject: [PATCH 14/15] Cleaning up unnecessary changes for cleaner comparison to main --- cases.csv | 1 - helpers/restart_runs.py | 3 +- inputs/emission_constraints/co2_tax.csv | 84 ++++++++++++------------- reeds/core/setup/b_inputs.gms | 30 +++------ reeds/core/setup/c_model.gms | 12 ++-- reeds/core/setup/d_objective.gms | 15 ++--- reeds/core/solve/3_solve_oneyear.gms | 1 - reeds/core/solve/5_varfix.gms | 6 -- reeds/core/terminus/report.gms | 69 ++++++++++---------- runreeds.py | 14 ++--- 10 files changed, 102 insertions(+), 133 deletions(-) diff --git a/cases.csv b/cases.csv index b7a343d1..ece24ca2 100644 --- a/cases.csv +++ b/cases.csv @@ -121,7 +121,6 @@ GSw_CO2_pipeline_fom,Set FOM cost for CO2 Pipeline. 2004$/(tonne-mi)/hr-yr. Sugg GSw_CO2_spurline_cost,Set capex cost for CO2 Spurline. 2004$/(tonne-mi)/hr. Suggested value: 2257,N/A,2257, GSw_CO2_spurline_fom,Set FOM cost for CO2 Spurline. 2004$/(tonne-mi)/hr-yr. Suggested value: 60,N/A,60, GSw_CO2_Detail,Switch to enable detailed representation CO2 transportation and storage,0; 1,0, -GSw_CO2_LimitStorageSites,Switch to limit storage sites to the closest site for each region,0; 1,0, GSw_CO2_CostAdj,multiplier for co2 infrastructure when Sw_CO2_Detail = 1,float,1, GSw_CO2_BEC,"Break even cost capacity factor assignment. Must be a suffix to ""bec_"" in the columns of co2_site_char.csv",N/A,70, GSw_CoalIGCC,Turn on/off COAL IGCC,0; 1,1, diff --git a/helpers/restart_runs.py b/helpers/restart_runs.py index 9b02315c..d353d2ae 100644 --- a/helpers/restart_runs.py +++ b/helpers/restart_runs.py @@ -181,8 +181,7 @@ # It is a single case or we are not on HPC if copy_srun_template: writelines_srun_case = writelines_srun.copy() - writelines_srun_case.append(f"#SBATCH --job-name={casename}") - writelines_srun_case.append(f"#SBATCH --output={os.path.join(case, 'slurm-%j.out')}\n") + writelines_srun_case.append(f"\n#SBATCH --job-name={casename}\n") writelines_srun_case.append(f"sh {callfile}") with open(sbatchfile, 'w') as f: for line in writelines_srun_case: diff --git a/inputs/emission_constraints/co2_tax.csv b/inputs/emission_constraints/co2_tax.csv index 1e08f624..d9e79551 100644 --- a/inputs/emission_constraints/co2_tax.csv +++ b/inputs/emission_constraints/co2_tax.csv @@ -1,42 +1,42 @@ -t,default,t050_2035,t050_2050,t100_2035,t100_2050,t150_2035,t150_2050,t200_2035,t200_2050,t200alt_2050 -2010,0,0,0,0,0,0,0,0,0,0 -2011,0,0,0,0,0,0,0,0,0,0 -2012,0,0,0,0,0,0,0,0,0,0 -2013,0,0,0,0,0,0,0,0,0,0 -2014,0,0,0,0,0,0,0,0,0,0 -2015,0,0,0,0,0,0,0,0,0,0 -2016,0,0,0,0,0,0,0,0,0,0 -2017,0,0,0,0,0,0,0,0,0,0 -2018,0,0,0,0,0,0,0,0,0,0 -2019,0,0,0,0,0,0,0,0,0,0 -2020,0,0,0,0,0,0,0,0,0,0 -2021,0,0,0,0,0,0,0,0,0,0 -2022,30,1,1,1,1,1,1,1,1,0 -2023,31.5,4.77,2.75,8.62,4.54,12.46,6.32,16.31,8.11,0 -2024,33.075,8.54,4.5,16.23,8.07,23.92,11.64,31.62,15.21,0 -2025,34.72875,12.31,6.25,23.85,11.61,35.38,16.96,46.92,22.32,0 -2026,36.4651875,16.08,8,31.46,15.14,46.85,22.29,62.23,29.43,0 -2027,38.28844688,19.85,9.75,39.08,18.68,58.31,27.61,77.54,36.54,0 -2028,40.20286922,23.62,11.5,46.69,22.21,69.77,32.93,92.85,43.64,0 -2029,42.21301268,27.38,13.25,54.31,25.75,81.23,38.25,108.15,50.75,0 -2030,44.32366331,31.15,15,61.92,29.29,92.69,43.57,123.46,57.86,0 -2031,46.53984648,34.92,16.75,69.54,32.82,104.15,48.89,138.77,64.96,10 -2032,48.8668388,38.69,18.5,77.15,36.36,115.62,54.21,154.08,72.07,20 -2033,51.31018074,42.46,20.25,84.77,39.89,127.08,59.54,169.38,79.18,30 -2034,53.87568978,46.23,22,92.38,43.43,138.54,64.86,184.69,86.29,40 -2035,56.56947427,50,23.75,100,46.96,150,70.18,200,93.39,50 -2036,59.39794798,50,25.5,100,50.5,150,75.5,200,100.5,60 -2037,62.36784538,50,27.25,100,54.04,150,80.82,200,107.61,70 -2038,65.48623765,50,29,100,57.57,150,86.14,200,114.71,80 -2039,68.76054953,50,30.75,100,61.11,150,91.46,200,121.82,90 -2040,72.19857701,50,32.5,100,64.64,150,96.79,200,128.93,100 -2041,75.80850586,50,34.25,100,68.18,150,102.11,200,136.04,110 -2042,79.59893115,50,36,100,71.71,150,107.43,200,143.14,120 -2043,83.57887771,50,37.75,100,75.25,150,112.75,200,150.25,130 -2044,87.7578216,50,39.5,100,78.79,150,118.07,200,157.36,140 -2045,92.14571268,50,41.25,100,82.32,150,123.39,200,164.46,150 -2046,96.75299831,50,43,100,85.86,150,128.71,200,171.57,160 -2047,101.5906482,50,44.75,100,89.39,150,134.04,200,178.68,170 -2048,106.6701806,50,46.5,100,92.93,150,139.36,200,185.79,180 -2049,112.0036897,50,48.25,100,96.46,150,144.68,200,192.89,190 -2050,117.6038742,50,50,100,100,150,150,200,200,200 +t,default,t050_2035,t050_2050,t100_2035,t100_2050,t150_2035,t150_2050,t200_2035,t200_2050 +2010,0,0,0,0,0,0,0,0,0 +2011,0,0,0,0,0,0,0,0,0 +2012,0,0,0,0,0,0,0,0,0 +2013,0,0,0,0,0,0,0,0,0 +2014,0,0,0,0,0,0,0,0,0 +2015,0,0,0,0,0,0,0,0,0 +2016,0,0,0,0,0,0,0,0,0 +2017,0,0,0,0,0,0,0,0,0 +2018,0,0,0,0,0,0,0,0,0 +2019,0,0,0,0,0,0,0,0,0 +2020,0,0,0,0,0,0,0,0,0 +2021,0,0,0,0,0,0,0,0,0 +2022,30,1,1,1,1,1,1,1,1 +2023,31.5,4.77,2.75,8.62,4.54,12.46,6.32,16.31,8.11 +2024,33.075,8.54,4.5,16.23,8.07,23.92,11.64,31.62,15.21 +2025,34.72875,12.31,6.25,23.85,11.61,35.38,16.96,46.92,22.32 +2026,36.4651875,16.08,8,31.46,15.14,46.85,22.29,62.23,29.43 +2027,38.28844688,19.85,9.75,39.08,18.68,58.31,27.61,77.54,36.54 +2028,40.20286922,23.62,11.5,46.69,22.21,69.77,32.93,92.85,43.64 +2029,42.21301268,27.38,13.25,54.31,25.75,81.23,38.25,108.15,50.75 +2030,44.32366331,31.15,15,61.92,29.29,92.69,43.57,123.46,57.86 +2031,46.53984648,34.92,16.75,69.54,32.82,104.15,48.89,138.77,64.96 +2032,48.8668388,38.69,18.5,77.15,36.36,115.62,54.21,154.08,72.07 +2033,51.31018074,42.46,20.25,84.77,39.89,127.08,59.54,169.38,79.18 +2034,53.87568978,46.23,22,92.38,43.43,138.54,64.86,184.69,86.29 +2035,56.56947427,50,23.75,100,46.96,150,70.18,200,93.39 +2036,59.39794798,50,25.5,100,50.5,150,75.5,200,100.5 +2037,62.36784538,50,27.25,100,54.04,150,80.82,200,107.61 +2038,65.48623765,50,29,100,57.57,150,86.14,200,114.71 +2039,68.76054953,50,30.75,100,61.11,150,91.46,200,121.82 +2040,72.19857701,50,32.5,100,64.64,150,96.79,200,128.93 +2041,75.80850586,50,34.25,100,68.18,150,102.11,200,136.04 +2042,79.59893115,50,36,100,71.71,150,107.43,200,143.14 +2043,83.57887771,50,37.75,100,75.25,150,112.75,200,150.25 +2044,87.7578216,50,39.5,100,78.79,150,118.07,200,157.36 +2045,92.14571268,50,41.25,100,82.32,150,123.39,200,164.46 +2046,96.75299831,50,43,100,85.86,150,128.71,200,171.57 +2047,101.5906482,50,44.75,100,89.39,150,134.04,200,178.68 +2048,106.6701806,50,46.5,100,92.93,150,139.36,200,185.79 +2049,112.0036897,50,48.25,100,96.46,150,144.68,200,192.89 +2050,117.6038742,50,50,100,100,150,150,200,200 \ No newline at end of file diff --git a/reeds/core/setup/b_inputs.gms b/reeds/core/setup/b_inputs.gms index c45bda48..771cc656 100644 --- a/reeds/core/setup/b_inputs.gms +++ b/reeds/core/setup/b_inputs.gms @@ -2432,7 +2432,6 @@ h2_ptc(i,v,r,t)$valcap(i,v,r,t) = h2_ptc_in(i,v,t) ; * Otherwise, we assume it receives $0/kg because the cleanliness of its carbon cannot be proven h2_ptc("electrolyzer",v,r,t)$[(not Sw_H2_PTC)] = 0; - set h2_ptc_years(t) "years in which the hydrogen production incentive is active"; h2_ptc_years(t) = tmodel_new(t)$[sum{(i,v,r),h2_ptc(i,v,r,t)}]; @@ -3130,7 +3129,6 @@ $include inputs_case%ds%tsc_binwidth.csv $offdelim $onlisting / ; -tsc_binwidth(r,rr,tscbin) = 1e9; parameter tsc_forward(r,rr,tscbin) "--$/MW-- transmission upgrade cost for forward direction" / @@ -5950,19 +5948,18 @@ $include inputs_case%ds%r_cs_distance_mi.csv $offdelim $ondigit $onlisting -/ , - min_r_cs_distance(r) "--mi-- minimum euclidean distance between BA transmission endpoints and storage formations" +/ ; $offempty -* find the closest storage site to each region -min_r_cs_distance(r) = smin(cs$[r_cs(r,cs)], r_cs_distance(r,cs)); +* Assign spurline costs +cost_co2_spurline_cap(r,cs,t)$[r_cs(r,cs)$tmodel_new(t)] = Sw_CO2_spurline_cost * r_cs_distance(r,cs) ; + +* CO2 pipelines can be build between any two adjacent BAs +cost_co2_pipeline_cap(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_cost * pipeline_distance(r,rr) ; +cost_co2_pipeline_fom(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_fom * pipeline_distance(r,rr) ; -$ifthene.rcslimit %GSw_CO2_LimitStorageSites% == 1 -* remove the region-site combinations that are not the closest -r_cs_distance(r,cs)$[r_cs(r,cs)$(r_cs_distance(r,cs) <= min_r_cs_distance(r))] = min_r_cs_distance(r) ; -r_cs_distance(r,cs)$[r_cs(r,cs)$(r_cs_distance(r,cs) > min_r_cs_distance(r))] = 0 ; -$endif.rcslimit +co2_routes(r,rr)$[routes_adjacent(r,rr)$pipeline_distance(r,rr)] = yes ; $onempty table co2_char(cs,*) "co2 site characteristics including injection rate limit, total storage limit, and break even cost" @@ -5981,18 +5978,9 @@ cost_co2_stor_bec(cs,t) = co2_char(cs,"bec_%GSw_CO2_BEC%"); csfeas(cs)$[co2_storage_limit(cs)$co2_injection_limit(cs)] = yes ; * only want to consider r_cs pairs which have available capacity r_cs(r,cs)$[not csfeas(cs)] = no ; -r_cs(r,cs)$[not r_cs_distance(r,cs)] = no ; -* Assign spurline costs -cost_co2_spurline_cap(r,cs,t)$[r_cs(r,cs)$tmodel_new(t)] = Sw_CO2_spurline_cost * r_cs_distance(r,cs) ; cost_co2_spurline_fom(r,cs,t)$[r_cs(r,cs)$tmodel_new(t)] = Sw_CO2_spurline_fom * r_cs_distance(r,cs) ; -* CO2 pipelines can be build between any two adjacent BAs -cost_co2_pipeline_cap(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_cost * pipeline_distance(r,rr) ; -cost_co2_pipeline_fom(r,rr,t)$[routes_adjacent(r,rr)$tmodel_new(t)] = Sw_CO2_pipeline_fom * pipeline_distance(r,rr) ; - -co2_routes(r,rr)$[routes_adjacent(r,rr)$pipeline_distance(r,rr)] = yes ; - cost_co2_pipeline_cap(r,rr,t) = %GSw_CO2_CostAdj% * cost_co2_pipeline_cap(r,rr,t); cost_co2_pipeline_fom(r,rr,t) = %GSw_CO2_CostAdj% * cost_co2_pipeline_fom(r,rr,t); cost_co2_stor_bec(cs,t) = %GSw_CO2_CostAdj% * cost_co2_stor_bec(cs,t) ; @@ -6167,4 +6155,4 @@ emit_rate(etype,e,i,v,r,t)$[not valgen(i,v,r,t)] = 0 ; *============================================================ valinv_init(i,v,r,t) = valinv(i,v,r,t) ; -valcap_init(i,v,r,t) = valcap(i,v,r,t) ; \ No newline at end of file +valcap_init(i,v,r,t) = valcap(i,v,r,t) ; diff --git a/reeds/core/setup/c_model.gms b/reeds/core/setup/c_model.gms index 3fa3758c..1f9e28a8 100644 --- a/reeds/core/setup/c_model.gms +++ b/reeds/core/setup/c_model.gms @@ -1048,6 +1048,7 @@ eq_rsc_INVlim(r,i,rscbin,t)$[tmodel(t) *plus exogenous (pre-start-year) capacity, using its level in the first year (tfirst) + sum{(ii,v,tt)$[tfirst(tt)$rsc_agg(i,ii)$exog_rsc(i)], capacity_exog_rsc(ii,v,r,rscbin,tt) } + ; * --------------------------------------------------------------------------- @@ -3864,7 +3865,7 @@ $endif.linked_co2_capture * --------------------------------------------------------------------------- -eq_co2_transport_caplimit(r,rr,h,t)$[co2_routes(r,rr)$Sw_CO2_Detail$h_rep(h) +eq_co2_transport_caplimit(r,rr,h,t)$[co2_routes(r,rr)$Sw_CO2_Detail $tmodel(t)$(yeart(t)>=co2_detail_startyr)].. *capacity computed as cumulative investments of co2 pipelines up to the current year @@ -3880,8 +3881,7 @@ eq_co2_transport_caplimit(r,rr,h,t)$[co2_routes(r,rr)$Sw_CO2_Detail$h_rep(h) * --------------------------------------------------------------------------- -eq_co2_spurline_caplimit(r,cs,h,t)$[Sw_CO2_Detail$r_cs(r,cs)$h_rep(h)$tmodel(t) - $(yeart(t)>=co2_detail_startyr)].. +eq_co2_spurline_caplimit(r,cs,h,t)$[Sw_CO2_Detail$r_cs(r,cs)$tmodel(t)$(yeart(t)>=co2_detail_startyr)].. *capacity computed as cumulative investments of co2 spurlines up to the current year sum{tt$[(yeart(tt)<=yeart(t))$(tmodel(tt) or tfix(tt))$(yeart(tt)>=co2_detail_startyr)], @@ -3899,7 +3899,7 @@ $endif.linked_co2_spurline_caplimit * --------------------------------------------------------------------------- -eq_co2_sink(r,h,t)$[tmodel(t)$Sw_CO2_Detail$h_rep(h)$(yeart(t)>=co2_detail_startyr)].. +eq_co2_sink(r,h,t)$[tmodel(t)$Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)].. *the amount of co2 stored from r in all of its cs sites sum{cs$r_cs(r,cs), CO2_STORED(r,cs,h,t) } @@ -3930,8 +3930,7 @@ $endif.linked_co2_sink * --------------------------------------------------------------------------- -eq_co2_injection_limit(cs,h,t)$[tmodel(t)$Sw_CO2_Detail$h_rep(h) - $(yeart(t)>=co2_detail_startyr)$csfeas(cs)].. +eq_co2_injection_limit(cs,h,t)$[Sw_CO2_Detail$tmodel(t)$(yeart(t)>=co2_detail_startyr)$csfeas(cs)].. * exogenously defined injection limit co2_injection_limit(cs) @@ -3973,7 +3972,6 @@ $ifthene.linked_co2_storage_cumul_limit Sw_FINITO_Link==1 yearweight(tt) * hours(h) * EXTRACT_CO2_CS(cs,r,h,tt) }$[t_finito(t)] $endif.linked_co2_storage_cumul_limit ; - * --------------------------------------------------------------------------- *=================== diff --git a/reeds/core/setup/d_objective.gms b/reeds/core/setup/d_objective.gms index 2400f944..af8672bf 100644 --- a/reeds/core/setup/d_objective.gms +++ b/reeds/core/setup/d_objective.gms @@ -259,12 +259,12 @@ eq_Objfn_op(t)$tmodel(t).. *Sw_GasCurve = 0 (census division supply curves natural gas prices) + sum{(cendiv,gb)$[tfuel(t)], sum{h, hours(h) * GASUSED(cendiv,gb,h,t) } * gasprice(cendiv,gb,t) - }$[(Sw_GasCurve = 0)] + }$(Sw_GasCurve = 0) *Sw_GasCurve = 3 (national supply curve for natural gas prices with census division multipliers) + sum{(h,cendiv,gb)$[tfuel(t)], hours(h) * GASUSED(cendiv,gb,h,t) * gasadder_cd(cendiv,t,h) + gasprice_nat_bin(gb,t) - }$[(Sw_GasCurve = 3)] + }$(Sw_GasCurve = 3) *Sw_GasCurve = 1 (national and census division supply curves for natural gas prices) *first - anticipated costs of gas consumption given last year's amount @@ -279,12 +279,13 @@ eq_Objfn_op(t)$tmodel(t).. + sum{(fuelbin)$tfuel(t), gasbinp_national(fuelbin,t) * VGASBINQ_NATIONAL(fuelbin,t) } - )$[(Sw_GasCurve = 1)] + )$[Sw_GasCurve = 1] * ---cost of biofuel consumption and biomass transport--- + sum{(r,bioclass)$[tfuel(t)$sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }], BIOUSED(bioclass,r,t) * (sum{usda_region$r_usda(r,usda_region), biosupply(usda_region, bioclass, "price") } + bio_transport_cost) } + * --- hurdle costs for transmission flow --- + sum{(r,rr,h,trtype)$[routes(r,rr,trtype,t)$cost_hurdle(r,rr,t)], cost_hurdle(r,rr,t) * FLOW(r,rr,h,t,trtype) * hours(h) } @@ -349,15 +350,15 @@ eq_Objfn_op(t)$tmodel(t).. CO2_SPURLINE_INV(r,cs,tt) } }$[Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)] * --- CO2 injection break even costs - + sum{(r,cs,h)$[r_cs(r,cs)$h_rep(h)], hours(h) * CO2_STORED(r,cs,h,t) * cost_co2_stor_bec(cs,t) }$[Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)] + + sum{(r,cs,h)$r_cs(r,cs), hours(h) * CO2_STORED(r,cs,h,t) * cost_co2_stor_bec(cs,t) }$[Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)] * --- Tax credit for CO2 stored --- * note conversion to 12-year CRF given length of CO2 captured incentive payments - - sum{(i,v,r,h)$[valgen(i,v,r,t)$h_rep(h)$co2_captured_incentive(i,v,r,t)], + - sum{(i,v,r,h)$[valgen(i,v,r,t)$co2_captured_incentive(i,v,r,t)], (crf(t) / crf_co2_incentive(t)) * co2_captured_incentive(i,v,r,t) * hours(h) * capture_rate("CO2",i,v,r,t) * GEN(i,v,r,h,t)} * --- Tax credit for CO2 stored for DAC --- - - sum{(p,i,v,r,h)$[dac(i)$valcap(i,v,r,t)$i_p(i,p)$h_rep(h)$co2_captured_incentive(i,v,r,t)], + - sum{(p,i,v,r,h)$[dac(i)$valcap(i,v,r,t)$i_p(i,p)$h_rep(h)], (crf(t) / crf_co2_incentive(t)) * co2_captured_incentive(i,v,r,t) * hours(h) * PRODUCE(p,i,v,r,h,t)} * --- PTC value for electric power generation --- @@ -371,7 +372,7 @@ eq_Objfn_op(t)$tmodel(t).. - sum{(p,v,r,h)$[valcap("electrolyzer",v,r,t)$(sameas(p,"H2"))$h2_ptc("electrolyzer",v,r,t)$h_rep(h)], hours(h) * PRODUCE(p,"electrolyzer",v,r,h,t) * (crf(t) / crf_h2_incentive(t)) * h2_ptc("electrolyzer",v,r,t) * 1e3} - $[(Sw_H2_PTC)$Sw_H2$h2_ptc_years(t)$(yeart(t) >= h2_demand_start)] + $[Sw_H2_PTC$Sw_H2$h2_ptc_years(t)$(yeart(t) >= h2_demand_start)] *end multiplier for pvf_onm ) diff --git a/reeds/core/solve/3_solve_oneyear.gms b/reeds/core/solve/3_solve_oneyear.gms index 43e6c747..d059ff3a 100644 --- a/reeds/core/solve/3_solve_oneyear.gms +++ b/reeds/core/solve/3_solve_oneyear.gms @@ -247,7 +247,6 @@ $endif.diagnose * ------------------------------ * Solve the Model * ------------------------------ -option limrow = 1e6; $ifthen.valstr %GSw_ValStr% == 1 OPTION lp = convert ; ReEDSmodel.optfile = 1 ; diff --git a/reeds/core/solve/5_varfix.gms b/reeds/core/solve/5_varfix.gms index 6230d69a..9b133772 100644 --- a/reeds/core/solve/5_varfix.gms +++ b/reeds/core/solve/5_varfix.gms @@ -9,10 +9,6 @@ if(Sw_RemoveSmallNumbers = 1, INV_ENERGY.l(i,v,r,tfix)$[abs(INV_ENERGY.l(i,v,r,tfix)) < rhs_tolerance] = 0 ; INV_RSC.l(i,v,r,rscbin,tfix)$[abs(INV_RSC.l(i,v,r,rscbin,tfix)) < rhs_tolerance] = 0 ; INV_POI.l(r,tfix)$[abs(INV_POI.l(r,tfix)) < rhs_tolerance] = 0 ; - INV_REFURB.l(i,v,r,tfix)$[abs(INV_REFURB.l(i,v,r,tfix)) < rhs_tolerance] = 0; - CO2_SPURLINE_INV.l(r,cs,tfix)$[abs(CO2_SPURLINE_INV.l(r,cs,tfix)) < rhs_tolerance] = 0; - CO2_TRANSPORT_INV.l(r,rr,tfix)$[abs(CO2_TRANSPORT_INV.l(r,rr,tfix)) < rhs_tolerance] = 0; - H2_TRANSPORT_INV.l(r,rr,tfix)$[abs(H2_TRANSPORT_INV.l(r,rr,tfix)) < rhs_tolerance] = 0 ; H2_STOR_INV.l(h2_stor,r,tfix)$[abs(H2_STOR_INV.l(h2_stor,r,tfix)) < rhs_tolerance] = 0 ; H2_TRANSPORT_INV.l(r,rr,tfix) $[abs(H2_TRANSPORT_INV.l(r,rr,tfix) ) < rhs_tolerance] = 0 ; ); @@ -123,9 +119,7 @@ H2_STOR_LEVEL_SZN.fx(h2_stor,r,actualszn,tfix)$[(h2_stor_r(h2_stor,r))$(Sw_H2=2) *CO2-related variables CO2_CAPTURED.fx(r,h,tfix)$Sw_CO2_Detail = CO2_CAPTURED.l(r,h,tfix) ; -$ifthene.linked_varfix Sw_FINITO_Link==0 CO2_STORED.fx(r,cs,h,tfix)$[Sw_CO2_Detail$r_cs(r,cs)] = CO2_STORED.l(r,cs,h,tfix) ; -$endif.linked_varfix CO2_FLOW.fx(r,rr,h,tfix)$[Sw_CO2_Detail$co2_routes(r,rr)] = CO2_FLOW.l(r,rr,h,tfix) ; CO2_TRANSPORT_INV.fx(r,rr,tfix)$[Sw_CO2_Detail$co2_routes(r,rr)] = CO2_TRANSPORT_INV.l(r,rr,tfix) ; CO2_SPURLINE_INV.fx(r,cs,tfix)$[Sw_CO2_Detail$r_cs(r,cs)] = CO2_SPURLINE_INV.l(r,cs,tfix) ; diff --git a/reeds/core/terminus/report.gms b/reeds/core/terminus/report.gms index cdc0b943..5a0f5e8d 100644 --- a/reeds/core/terminus/report.gms +++ b/reeds/core/terminus/report.gms @@ -495,51 +495,46 @@ ptc_out(i,v,t)$[tmodel_new(t)$ptc_value_scaled(i,v,t)] = ptc_value_scaled(i,v,t) * Case 2: the resource of one or more biomass classes ARE exhausted, i.e., BIOUSED.l(bioclass) = biosupply(bioclass) * Marginal Biomass Price = maximum difference between eq_bioused.m and eq_biousedlimit.m(bioclass) across all biomass classes in a region -repbioprice(r,t)$[tmodel_new(t)$tfuel(t)] = - [max{0, smax{bioclass$BIOUSED.l(bioclass,r,t), eq_bioused.m(r,t) - - sum{usda_region$r_usda(r,usda_region), eq_biousedlimit.m(bioclass,usda_region,t) } } } / pvf_onm(t) ] ; +repbioprice(r,t)$[tmodel_new(t)$tfuel(t)] = max{0, smax{bioclass$BIOUSED.l(bioclass,r,t), eq_bioused.m(r,t) - + sum{usda_region$r_usda(r,usda_region), eq_biousedlimit.m(bioclass,usda_region,t) } } } / pvf_onm(t) ; -* quantity of biomass used in the power sector (convert from mmBTU to dry tons using biomass energy content) -bioused_out(bioclass,r,t)$[tmodel_new(t)$tfuel(t)] = - [BIOUSED.l(bioclass,r,t) / bio_energy_content ]; - -bioused_usda(bioclass,usda_region,t)$[tmodel_new(t)$tfuel(t)] = - [sum{r$r_usda(r,usda_region), bioused_out(bioclass,r,t) } ]; +* quantity of biomass used (convert from mmBTU to dry tons using biomass energy content) +bioused_out(bioclass,r,t)$[tmodel_new(t)$tfuel(t)] = BIOUSED.l(bioclass,r,t) / bio_energy_content ; +bioused_usda(bioclass,usda_region,t)$[tmodel_new(t)$tfuel(t)] = sum{r$r_usda(r,usda_region), bioused_out(bioclass,r,t) } ; * 1e9 converts from MMBtu to Quads repgasquant(cendiv,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 3)$tmodel_new(t)$tfuel(t)] = - [sum{(gb,h), GASUSED.l(cendiv,gb,h,t) * hours(h) } * gas_scale/ 1e9 ]; + sum{(gb,h), GASUSED.l(cendiv,gb,h,t) * hours(h) } * gas_scale/ 1e9 ; repgasquant(cendiv,t)$[(Sw_GasCurve = 1 or Sw_GasCurve = 2 or Sw_FINITO_Link = 1)$tmodel_new(t)] = - [( sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], + ( sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t)} + sum{(v,r,h)$[valcap("dac_gas",v,r,t)$r_cendiv(r,cendiv)], hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE.l("DAC","dac_gas",v,r,h,t) }$Sw_DAC_Gas + sum{(p,i,v,r,h)$[r_cendiv(r,cendiv)$valcap(i,v,r,t)$smr(i)], hours(h) * smr_methane_rate * PRODUCE.l(p,i,v,r,h,t) }$Sw_H2 - ) / 1e9 ]; + ) / 1e9 ; -repgasquant_irt(i,r,t)$[tmodel_new(t)] = - [( sum{(v,h)$[valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], +repgasquant_irt(i,r,t)$tmodel_new(t) = + ( sum{(v,h)$[valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } + sum{(v,h)$[valcap("dac_gas",v,r,t)], hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE.l("DAC","dac_gas",v,r,h,t) }$Sw_DAC_Gas + sum{(p,v,h)$[valcap(i,v,r,t)$smr(i)], hours(h) * smr_methane_rate * PRODUCE.l(p,i,v,r,h,t) }$Sw_H2 - ) / 1e9 ]; + ) / 1e9 ; -repgasquant_nat(t)$[tmodel_new(t)] = - [sum{cendiv, repgasquant(cendiv,t) } ]; +repgasquant_nat(t)$tmodel_new(t) = sum{cendiv, repgasquant(cendiv,t) } ; *for reported gasprice (not that used to compute system costs) *scale back to $ / mmbtu -repgasprice(cendiv,t)$[(Sw_GasCurve = 0)$tmodel_new(t)$repgasquant(cendiv,t)$tfuel(t)]= - [smax{gb$[sum{h, GASUSED.l(cendiv,gb,h,t) }], gasprice(cendiv,gb,t) } / gas_scale ]; +repgasprice(cendiv,t)$[(Sw_GasCurve = 0)$tmodel_new(t)$repgasquant(cendiv,t)$tfuel(t)] = + smax{gb$[sum{h, GASUSED.l(cendiv,gb,h,t) }], gasprice(cendiv,gb,t) } / gas_scale ; repgasprice(cendiv,t)$[(Sw_GasCurve = 2)$tmodel_new(t)$repgasquant(cendiv,t)$tfuel(t)] = - [sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], + sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)], hours(h)*heat_rate(i,v,r,t)*fuel_price(i,r,t)*GEN.l(i,v,r,h,t) - } / (repgasquant(cendiv,t) * 1e9) ]; + } / (repgasquant(cendiv,t) * 1e9) ; * gas price by timeslice when linked with FINITO (see similar calculation in finito_report.gms) [$2004/MMBtu] @@ -579,38 +574,40 @@ repgasprice_r(r,t)$[(Sw_GasCurve = 1)$tmodel_new(t)$tfuel(t)] = ; repgasprice(cendiv,t)$[((Sw_GasCurve = 1) or (Sw_FINITO_Link = 1))$tmodel_new(t)$repgasquant(cendiv,t)] = - [sum{(i,r)$r_cendiv(r,cendiv), repgasprice_r(r,t) * repgasquant_irt(i,r,t) } / repgasquant(cendiv,t) ] ; + sum{(i,r)$r_cendiv(r,cendiv), repgasprice_r(r,t) * repgasquant_irt(i,r,t) } / repgasquant(cendiv,t) ; repgasprice_nat(t)$[tmodel_new(t)$sum{cendiv, repgasquant(cendiv,t) }] = - [sum{cendiv, repgasprice(cendiv,t) * repgasquant(cendiv,t) } - / sum{cendiv, repgasquant(cendiv,t) } ]; + sum{cendiv, repgasprice(cendiv,t) * repgasquant(cendiv,t) } + / sum{cendiv, repgasquant(cendiv,t) } ; *======================================== * NATURAL GAS FUEL COSTS *======================================== gasshare_ba(r,cendiv,t)$[r_cendiv(r,cendiv)$tmodel_new(t)$repgasquant(cendiv,t)] = - [sum{i$[valgen_irt(i,r,t)$gas(i)],repgasquant_irt(i,r,t) / repgasquant(cendiv,t) } ]; + sum{i$[valgen_irt(i,r,t)$gas(i)],repgasquant_irt(i,r,t) / repgasquant(cendiv,t) } ; gasshare_techba(i,r,cendiv,t)$[r_cendiv(r,cendiv)$tmodel_new(t)$repgasquant(cendiv,t)$gas(i)] = - [repgasquant_irt(i,r,t) / repgasquant(cendiv,t) ]; + repgasquant_irt(i,r,t) / repgasquant(cendiv,t) ; gasshare_cendiv(cendiv,t)$[sum{cendiv2,repgasquant(cendiv2,t)}] = - [repgasquant(cendiv,t) / sum{cendiv2,repgasquant(cendiv2,t)} ]; + repgasquant(cendiv,t) / sum{cendiv2,repgasquant(cendiv2,t)} ; * cost of natural gas - standalone ReEDS gascost_cendiv(cendiv,t)$[tmodel_new(t)$tfuel(t)] = *cost of natural gas for Sw_GasCurve = 2 (static natural gas prices) - + [ sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t) - $[not bio(i)]$[not cofire(i)]$[Sw_GasCurve = 2]], - hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN.l(i,v,r,h,t) } + + sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t) + $[not bio(i)]$[not cofire(i)]$[Sw_GasCurve = 2]], + hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN.l(i,v,r,h,t) } *cost of natural gas for Sw_GasCurve = 0 (census division supply curves natural gas prices) - + sum{gb, sum{h,hours(h) * GASUSED.l(cendiv,gb,h,t) } * gasprice(cendiv,gb,t) }$[Sw_GasCurve = 0] + + sum{gb, sum{h,hours(h) * GASUSED.l(cendiv,gb,h,t) } * gasprice(cendiv,gb,t) + }$[Sw_GasCurve = 0] *cost of natural gas for Sw_GasCurve = 3 (national supply curve for natural gas prices with census division multipliers) - + sum{(h,gb), hours(h) * GASUSED.l(cendiv,gb,h,t) - * gasadder_cd(cendiv,t,h) + gasprice_nat_bin(gb,t) }$[Sw_GasCurve = 3] + + sum{(h,gb), hours(h) * GASUSED.l(cendiv,gb,h,t) + * gasadder_cd(cendiv,t,h) + gasprice_nat_bin(gb,t) + }$[Sw_GasCurve = 3] *cost of natural gas for Sw_GasCurve = 1 (national and census division supply curves for natural gas prices) *first - anticipated costs of gas consumption given last year's amount + (sum{(i,v,r,h)$[valgen(i,v,r,t)$gas(i)], @@ -623,8 +620,7 @@ gascost_cendiv(cendiv,t)$[tmodel_new(t)$tfuel(t)] = + sum{(fuelbin), gasbinp_national(fuelbin,t) * VGASBINQ_NATIONAL.l(fuelbin,t) } * gasshare_cendiv(cendiv,t) - )$[Sw_GasCurve = 1] - ]; + )$[Sw_GasCurve = 1]; * cost of natural gas - linked with FINITO ('not tfuel' indicates years using FINITO supply curves) gascost_cendiv(cendiv,t)$[tmodel_new(t)$(not tfuel(t))] = @@ -641,7 +637,7 @@ gascost_cendiv(cendiv,t)$[tmodel_new(t)$(not tfuel(t))] = bioshare_techba(i,r,t)$[(cofire(i) or bio(i))$tmodel_new(t)$tfuel(t)] = * biofuel-based generation of tech i in the BA (biopower + cofire) - [ (( sum{(v,h)$[valgen(i,v,r,t)$bio(i)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } + (( sum{(v,h)$[valgen(i,v,r,t)$bio(i)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } + sum{(v,h)$[cofire(i)$valgen(i,v,r,t)], bio_cofire_perc * hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) } ) / * biofuel-based generation of all techs in the BA (biopower + cofire) @@ -651,7 +647,6 @@ bioshare_techba(i,r,t)$[(cofire(i) or bio(i))$tmodel_new(t)$tfuel(t)] = )$[ sum{(ii,v,h)$[valgen(ii,v,r,t)$bio(ii)], hours(h) * heat_rate(ii,v,r,t) * GEN.l(ii,v,r,h,t) } + sum{(ii,v,h)$[cofire(ii)$valgen(ii,v,r,t)], bio_cofire_perc * hours(h) * heat_rate(ii,v,r,t) * GEN.l(ii,v,r,h,t) } ] - ] ; *========================= diff --git a/runreeds.py b/runreeds.py index fa7fea75..0b07da7f 100644 --- a/runreeds.py +++ b/runreeds.py @@ -13,8 +13,6 @@ import csv import importlib import numpy as np -import gdxpds -from gdxpds.gdx import GdxFile,GdxSymbol,GamsDataType import pandas as pd import subprocess import re @@ -1358,7 +1356,6 @@ def write_batch_script( OPATH.writelines("module load conda \n") OPATH.writelines("module load gams \n") - OPATH.writelines("conda deactivate \n") OPATH.writelines("conda activate reeds2 \n") OPATH.writelines('export R_LIBS_USER="$HOME/rlib" \n\n\n') @@ -1438,8 +1435,7 @@ def write_batch_script( ################################# # -- OUTPUT PROCESSING -- # ################################# - - ## create reporting files + #create reporting files big_comment('Output processing', OPATH) if not LINUXORMAC: OPATH.writelines("setlocal enabledelayedexpansion\n") @@ -1513,7 +1509,7 @@ def write_batch_script( + f" --casepath {casedir} \n\n" ) - ## Run the retail rate module + ### Run the retail rate module OPATH.writelines( "python" + f" {os.path.join(reeds_path,'postprocessing','retail_rate_module','retail_rate_calculations.py')}" @@ -1527,7 +1523,7 @@ def write_batch_script( f"{casedir}\n\n" ) - ## Make script to unload all data to .gdx file + ### Make script to unload all data to .gdx file command = ( f"gams {Path('reeds','core','terminus','dump_alldata.gms')}" + ' o='+os.path.join('lstfiles','dump_alldata_{}_{}.lst'.format(BatchName,case)) @@ -1816,7 +1812,7 @@ def launch_single_case_run( shellscript = subprocess.Popen( [os.path.join(casedir, 'call_' + batch_case + ext)], shell=True) # Wait for it to finish before killing the thread - #shellscript.wait() + shellscript.wait() else: if int(caseSwitches['keep_run_terminal']) == 1: terminal_keep_flag = ' /k ' @@ -1944,7 +1940,7 @@ def main( help="Check inputs but don't start runs") args = parser.parse_args() - + main( BatchName=args.BatchName, cases_suffix=args.cases_suffix, single=args.single, simult_runs=args.simult_runs, forcelocal=args.forcelocal, skip_checks=args.skip_checks, From 93bd58bc4f2df90e6b22ce350a4ca29c144f99dc Mon Sep 17 00:00:00 2001 From: mmowers Date: Mon, 29 Jun 2026 15:09:08 -0400 Subject: [PATCH 15/15] Remove old switch from cases_finito.csv --- cases_finito.csv | 1 - 1 file changed, 1 deletion(-) diff --git a/cases_finito.csv b/cases_finito.csv index ca5459a3..6e6c732e 100644 --- a/cases_finito.csv +++ b/cases_finito.csv @@ -34,7 +34,6 @@ GSw_DAC,1,,,,,,,,,,,,,,,,,,,,, GSw_BECCS,1,,,,,,,,,,,,,,,,,,,,, GSw_TCPhaseout,1,,,,,,,,,,,,,,,,,,,,, GSw_CO2_Detail,1,,,,,,,,,,,,,,,,,,,,, -GSw_CO2_LimitStorageSites,1,,,,,,,,,,,,,,,,,,,,, numhintage,4,,,,,,,,,,,,,,,,,,,,, GSw_CarbTax,0,,,,,,,,,,,1,1,,,,,,,,, GSw_CarbTaxOption,default,,,,,,,,,,,t200_2050,t200_2050,,,,,,,,,