Derate thermal capacity reserve margin contribution with FOR

CONTRIBUTING.md

@@ -6,4 +6,4 @@ For questions or general discussion, please open a [discussion](https://github.c
 
 To contribute code, fork the repository, make your changes, and submit a pull request. Please review our [Developer Guide](https://reeds-model.github.io/ReEDS/developer_best_practices.html) before getting started.
 
-**Note:** if you're considering making a change that might lead to meaningful differences in model capability, workflow, or outputs, we recommend initiating a [discussion](https://github.com/ReEDS-Model/ReEDS/discussions).
+**Note:** if you're considering making a change that might lead to meaningful differences in model capability, workflow, or outputs, we recommend initiating a [discussion](https://github.com/ReEDS-Model/ReEDS/discussions).

cases_test.csv

@@ -19,6 +19,7 @@ GSw_RetirePenalty,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,
 GSw_FakeData,,,,,,1,1,1,,,,,,,,,,,,,,,,,,,,,,
 GSw_PRM_CapCredit,,,,,,,,1,1,,,,,,,,1,,,,,,,,,,,,,
 GSw_PRM_scenario,,,,,,,,,static,,,,,,,,static,,,,,,,,,,,,,
+GSw_PRM_UpdateMethod,,,,,,,,,1,,,,,,,,,,,,,,,,,,,,,
 GSw_HourlyType,,,,,,,,,,wek,year,,,,,,,,,,,,,,,,,,,
 GSw_InterDayLinkage,,,,,,,,,,,,1,,,,,,,,,,,,,,,,,,
 GSw_HourlyWeatherYears,,,,,,,2012_2013,,,,,,2020,2007_2008_2009_2010_2011_2012_2013,,,,,,,,2012_2013,,,,2018,,,,

docs/source/model_documentation.md

@@ -2460,11 +2460,12 @@ Adapted from {cite}`maiIncorporatingStressfulGrid2024`.
 
 The calculation of capacity credit for VRE is described in the [VRE capacity credit](#vre-capacity-credit) section;
 the method for storage is described in the [storage capacity credit](#storage-capacity-credit) section.
-Thermal generators are given a capacity credit of 100% in ReEDS;
-technology-specific [outage rates](#outage-rates) for thermal generators are not considered in this method
-and are instead assumed to be factored into the [planning reserve margin](#planning-reserve-margins).
 
 
+#### Thermal generation capacity credit
+
+For thermal generators (i.e. combined cycle, combustion turbine, nuclear (conventional and SMR), and steam (coal)), ReEDS estimates a seasonal capacity credit for each region/technology combination based on temperature-dependent forced outage rates (described in detail in [outage rates](#outage-rates)). To calculate each technology's contribution to the seasonal reserve margin, its nameplate capacity for each technology is multiplied by $(1 - \bar{FOR})$, where $\bar{FOR}$ is the mean forced outage rate during the top 20 net load hours cross all modeled [weather years](#weather-years) for each season. 
+
 
 #### VRE capacity credit
 

reeds/core/setup/b_inputs.gms

@@ -6118,7 +6118,8 @@ Parameter
     cc_excess(i,r,ccseason,t)              "--MW-- this is the excess capacity credit when assuming marginal capacity credit in intertemporal solve"
     vre_gen_last_year(r,allh,t)            "--MW-- generation from VRE generators in the prior solve year"
     hybrid_cc_derate(i,r,ccseason,sdbin,t) "--fraction-- derate factor for hybrid PV+battery storage capacity credit"
-    m_cc_mar(i,r,ccseason,t)               "--fraction-- marginal capacity credit",
+    m_cc_mar(i,r,ccseason,t)               "--fraction-- marginal capacity credit"
+    mean_forced_outage_rate(i,r,ccseason,t)"--fraction-- mean forced outage rate for each technology, region, and ccseason - used to derate thermal generator capacity"
 * Heuristic climate impacts
     trans_cap_delta(allh,allt)             "--fraction-- fractional adjustment to transmission capacity from climate heuristics"
 * Emissions and policies
@@ -6140,6 +6141,7 @@ cc_excess(i,r,ccseason,t) = 0 ;
 cc_old(i,r,ccseason,t) = 0 ;
 m_cc_mar(i,r,ccseason,t) = 0 ;
 hybrid_cc_derate(i,r,ccseason,sdbin,t)$[pvb(i)$valcap_irt(i,r,t)] = 1 ;
+mean_forced_outage_rate(i,r,ccseason,t) = 0 ;
 
 * Trim some of the largest matrices to reduce file sizes
 cost_vom(i,v,r,t)$[not valgen(i,v,r,t)] = 0 ;

reeds/core/setup/c_model.gms

@@ -1718,6 +1718,7 @@ eq_reserve_margin(r,ccseason,t)
     + sum{(i,v)$[valcap(i,v,r,t)$(not vre(i))$(not hydro(i))$(not storage(i))$(not consume(i))$(not forced_retire(i,r,t))],
           CAP(i,v,r,t)
           * (1 + ccseason_cap_frac_delta(i,v,r,ccseason,t))
+          * (1 - mean_forced_outage_rate(i,r,ccseason,t))
          }
 
 *[plus] firm capacity from existing VRE or CSP

reeds/core/setup/e_solveprep.gms

@@ -40,6 +40,7 @@ parameter
     sdbin_size_load(ccreg,ccseason,sdbin,t)  "--MW-- bin_size power loading in from the cc_out gdx file"
     cc_mar_load(i,r,ccreg,ccseason,t)        "--fraction-- cc_mar loading in from the cc_out gdx file"
     cc_evmc_load(i,r,ccseason,t)               "--fraction--  cc_evmc loading in from the cc_out gdx file"
+    mean_forced_outage_rate_load(i,r,ccseason,t)  "--fraction-- mean_forced_outage_rate loading in from the cc_out gdx file"
 ;
 
 

reeds/core/solve/3_solve_oneyear.gms

@@ -147,6 +147,7 @@ $loaddcr cc_old_load = cc_old
 $loaddcr cc_mar_load = cc_mar
 $loaddcr cc_evmc_load = cc_evmc
 $loaddcr sdbin_size_load = sdbin_size
+$loaddcr mean_forced_outage_rate_load = mean_forced_outage_rate
 $gdxin
 
 *Note: these values are rounded before they are written to the gdx file, so no need to round them here
@@ -181,6 +182,7 @@ cc_old(i,r,ccseason,t)$[tload(t)$pvb(i)$sum{(v,tt)$tprev(t,tt), CAP.l(i,v,r,tt)}
     min{ cc_old(i,r,ccseason,t) / sum{(v,tt)$tprev(t,tt), CAP.l(i,v,r,tt)}, 1 / ilr(i) - bcr(i) }
     * sum{(v,tt)$tprev(t,tt), CAP.l(i,v,r,tt)};
 
+mean_forced_outage_rate(i,r,ccseason,t)$tload(t) = mean_forced_outage_rate_load(i,r,ccseason,t) ;
 $endif.tcheck
 
 

reeds/core/terminus/report.gms

@@ -934,7 +934,10 @@ cc_new(i,r,ccseason,t)$[valcap_irt(i,r,t)$cap_new_cc(i,r,ccseason,t)] = sum{v$iv
 
 cap_firm(i,r,ccseason,t)$[valcap_irt(i,r,t)$[not consume(i)]$tmodel_new(t)$Sw_PRM_CapCredit] =
       sum{v$[(not vre(i))$(not hydro(i))$(not storage(i))$(not storage_hybrid(i)$(not csp(i)))$valcap(i,v,r,t)],
-          CAP.l(i,v,r,t) * (1 + ccseason_cap_frac_delta(i,v,r,ccseason,t)) }
+          CAP.l(i,v,r,t) 
+          * (1 + ccseason_cap_frac_delta(i,v,r,ccseason,t)) 
+          * (1 - mean_forced_outage_rate(i,r,ccseason,t))
+          }
     + cc_old(i,r,ccseason,t)
     + sum{v$[(vre(i) or csp(i) or storage_hybrid(i)$(not csp(i)))$valinv(i,v,r,t)],
          m_cc_mar(i,r,ccseason,t) * (INV.l(i,v,r,t) + INV_REFURB.l(i,v,r,t)$[refurbtech(i)$Sw_Refurb]) }

reeds/resource_adequacy/capacity_credit.py

@@ -152,13 +152,35 @@ def reeds_cc(t, tnext, casedir):
     ### Get the non-duplicated profiles
     resource_profiles = resources.drop_duplicates('resource')
 
+    ### Get forced outage profiles, flatten, and melt
+    forced_outage_rate = reeds.io.get_outage_hourly(inputs_case,'forced')
+
+    # drop techs that won't use this capacity credit before stacking
+    techs_to_drop = ['battery_li','can-imports','distpv', 'electrolyzer',
+                    'hydd', 'hyded', 'hydend', 'hydnd', 'hydnpd', 'hydnpnd',
+                    'hydro', 'hydsd', 'hydsnd', 'hydud', 'hydund', 'pumped-hydro']
+
+    forced_outage_rate = forced_outage_rate.drop(columns=techs_to_drop, errors='ignore')
+
+    forced_outage_rate = (forced_outage_rate
+                            .stack(level=1)
+                            .reset_index()
+                            .rename(columns={"level_0": "timestamp"})
+                            )
+
+    forced_outage_rate = pd.melt(
+        forced_outage_rate,
+        id_vars=['timestamp','r'], 
+    )
+
     # Remove the "8760" safety valve bin
     safety_bin = max(sdb['bin'].values)
     sdb = sdb[sdb['bin'] != max(sdb['bin'])]
     sdb = [int(x) for x in sdb['bin']]
 
     # Temporal definitions
     h_dt_szn = pd.read_csv(os.path.join('inputs_case', 'rep', 'h_dt_szn.csv'))
+    hmap_allyrs = pd.read_csv(os.path.join('inputs_case','rep','hmap_allyrs.csv'))
 
     ccseasons = []
     if sw['cc_calc_annual']:
@@ -464,6 +486,39 @@ def pivot_melt_data(df):
     ### Reorder to match ReEDS convention
     net_load_2012 = net_load_2012.reindex(['ccreg','ccseason','year','h','hour','t','value'], axis=1)
 
+    # export mean forced outage rate for top load hours in each region and ccseason
+    top_net_load_hours = net_load.groupby(['ccreg','ccseason']).head(int(sw['GSw_PRM_CapCreditHours']))
+    # map back to regions from ccreg for ReEDS input
+    top_net_load_hours = top_net_load_hours.merge(
+        hierarchy[['r','ccreg']], on='ccreg', how='left'
+    ).drop(['ccreg','year','h','t','value'], axis=1)
+
+    # map hour to timestamp
+    top_net_load_hours = top_net_load_hours.merge(
+        hmap_allyrs[['hour','*timestamp','actual_h']], on='hour', how='left'
+    ).drop('hour', axis=1)
+    top_net_load_hours['*timestamp'] = pd.to_datetime(top_net_load_hours['*timestamp'])
+    top_net_load_hours['h'] = top_net_load_hours['actual_h']
+
+    # join top hours to outage rates based on timestamp and region
+    forced_outage_rate_top_hours = top_net_load_hours.merge(
+        forced_outage_rate,
+        left_on=['*timestamp','r'],
+        right_on=['timestamp','r'],
+        how='left'
+    ).drop('*timestamp', axis=1)
+
+    # Add t index
+    forced_outage_rate_top_hours['t'] = str(tnext)
+
+    # Find mean forced outage rate across top hours for each region and ccseason
+    mean_forced_outage_rate = (
+        forced_outage_rate_top_hours
+        .groupby(['i','r','ccseason','t'], as_index=False).mean('value')
+        .sort_values(['i','r','ccseason','t'])
+        .reset_index(drop=True)
+    )
+
     if int(sw['GSw_EVMC']):
         cc_evmc = (
             pd.concat(dict_cc_evmc, axis=0)
@@ -484,6 +539,7 @@ def pivot_melt_data(df):
         'sdbin_size': sdbin_size,
         'net_load': net_load,
         'net_load_2012': net_load_2012,
+        'mean_forced_outage_rate' : mean_forced_outage_rate,
     }
 
     return cc_results

reeds/resource_adequacy/ra_calcs.py

@@ -138,6 +138,7 @@ def main(t, tnext, casedir, iteration=0):
             'cc_old': pd.DataFrame(columns=['i','r','ccreg','szn','t','Value']),
             'cc_evmc': pd.DataFrame(columns=['i','r','szn','t','Value']),
             'sdbin_size': pd.DataFrame(columns=['ccreg','szn','bin','t','Value']),
+            'mean_forced_outage_rate': pd.DataFrame(columns=['i','r','szn','t','Value']),
         }
 
     reeds.log.toc(tic=tic, year=t, process='ra/capacity_credit.py')
@@ -165,8 +166,18 @@ def main(t, tnext, casedir, iteration=0):
             )
 
     #%% Identify stress periods
-    print('identifying new stress periods...')
-    tic = datetime.datetime.now()
+    if (
+        int(sw['GSw_PRM_CapCredit']) == 0 
+        and ('user' not in sw['GSw_PRM_StressModel'].lower())
+        ):
+        print('identifying new stress periods...')
+        tic = datetime.datetime.now()
+    elif (
+        int(sw['GSw_PRM_CapCredit']) and (int(sw['GSw_PRM_StressIterateMax']) > 0)
+    ):
+        print('updating PRM values...')
+        tic = datetime.datetime.now()
+
     if (
         ('user' not in sw['GSw_PRM_StressModel'].lower())
         or ((int(sw.GSw_PRM_StressIterateMax)) and int(sw['GSw_PRM_CapCredit']))

runreeds.py

@@ -386,6 +386,16 @@ def check_compatibility(sw):
             except ValueError:
                 raise ValueError(err)
 
+    if int(sw['GSw_PRM_UpdateMethod']) == 0 and int(sw['GSw_PRM_CapCredit']) == 1 and int(sw['GSw_PRM_StressIterateMax']) > 0:
+        raise ValueError(
+            "The combination of GSw_PRM_UpdateMethod=0, GSw_PRM_CapCredit=1, "
+            "and GSw_PRM_StressIterateMax>0 is not supported."
+            "To iteratively update the PRM, set GSw_PRM_UpdateMethod to one of: "
+            "(1) static update set by GSw_PRM_UpdateFraction; "
+            "(2) dynamic update informed by PRAS; "
+            "(3) dynamic update but only after all new stress periods have been added"
+        )
+
     for bir in sw['GSw_PVB_BIR'].split('_'):
         if not (float(bir) >= 0):
             raise ValueError("Fix GSw_PVB_BIR")