diff --git a/cases.csv b/cases.csv index 02c1b34db..97f358ae3 100644 --- a/cases.csv +++ b/cases.csv @@ -37,7 +37,7 @@ 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, egsnearfieldsupplycurve,EGS (near field) Supply Curve,ATB_2023,ATB_2023, -ngscen,Natural gas price inputs (see inputs\fuelprices),N/A,AEO_2026_baseline, +ngscen,Natural gas price inputs,^AEO_(2025_(HOG|reference|LOG)|2026_(HOG|baseline|LOG))$,AEO_2026_baseline, pshsupplycurve,PSH Supply Curve,(8|10|12)hr_(ref|wEphemeral|wExist|wEph|wExist_wEph)_(apr2025|mar2024),12hr_ref_apr2025, pvbscen,PVB cost-sharing scenario,benchmark2020,benchmark2020, dacscen,Electrical DAC cost and performance inputs (see inputs/consume/[dacscen].csv starting with 'dac_elec_').,BVRE_2021_(high|mid|low),BVRE_2021_mid, diff --git a/cases_test.csv b/cases_test.csv index 805d5829d..eea0ebea7 100644 --- a/cases_test.csv +++ b/cases_test.csv @@ -1,6 +1,6 @@ ,Default Value,Pacific,USA_defaults,Mid_Case,USA_decarb,github_Pacific,github_Everything,github_MA_county_CC,Pacific_CC,Pacific_weks,Pacific_full_year,Interday_storage,Pacific_2020,Pacific_rep15,WY_county,WECC_county,PJM_county_CC,NYVT_mixed,OR_water,MonteCarlo_Random,MonteCarlo_LHS,Everything,Simple,USA_fast,USA_faster,MultiMetricRA,Pacific_DR,Pacific_MGA,Pacific_MGA_RV,Pacific_LoadSite95,MARICTNYNJPAOH_Offshore,R2P ignore,1,0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, -GSw_Region,cendiv/Pacific,,country/USA,country/USA,country/USA,,st/ID.WY.NE.IA.IL,st/MA,,,,,,,st/WY,interconnect/western,transreg/PJM,st/NY.VT,st/OR,st/NE.NY.PA,st/NE.NY.PA,st/ID.WY.NE.IA.IL,st/KS,country/USA,country/USA,st/NY.NJ,,,,,st/MA.RI.CT.NY.NJ.PA.OH, +GSw_Region,cendiv/Pacific,,country/USA,country/USA,country/USA,,st/ID.WY.NE.IA.IL,st/MA,,,,,,,st/WY,interconnect/western,transreg/PJM,st/NY.VT,st/OR,,,st/ID.WY.NE.IA.IL,st/KS,country/USA,country/USA,st/NY.NJ,,,,,st/MA.RI.CT.NY.NJ.PA.OH, endyear,2035,,2050,2050,2050,2030,2060,2030,,,,,,,,,,,2035,2030,2030,2060,2035,2050,2050,2050,,,,,, yearset,,,,,,,2010..2060..10,,,,,,,,,,,,,,,2010..2060..10,2010_2015_2020..2050..3,,2010_2025_2050,,,,,,, GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,z134,z54,z48,,,,,,, @@ -30,7 +30,7 @@ GSw_WaterUse,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,,,, resource_adequacy_years,,,,,,,2011_2012_2013_2021_2022_2023,,,,,,,,,,,,,,,2011_2012_2013_2021_2022_2023,,,,,,,,,, GSw_HourlyClusterAlgorithm,,,,,,,,,,,,,,,,,,,,user,user,,,,,,,,,,, MCS_runs,,,,,,,,,,,,,,,,,,,,2,2,,,,,,,,,,, -MCS_dist_groups,,,,,,,,,,,,,,,,,,,,tech.hydro.nuclear.gas.coal.load_country,upv_tri.nuclear_tri.ng_fuel_price_tri.load_country_unif,,,,,,,,,,, +MCS_dist_groups,,,,,,,,,,,,,,,,,,,,load_flat_state.gasprice_flat.siting_flat_state.upv.onswind.ofswind.battery.nuclear.geo,load_country_unif.gasprice_tri.siting_country_discrete_flat.upv_tri.nuclear_tri,,,,,,,,,,, MCS_lhs,,,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,, GSw_PRM_StressIterateMax,,,,,,,,0,,,,,,,,0,0,,,1,1,,,,,,,,,,, GSw_ReducedResource,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,,, diff --git a/inputs/capacity_exogenous/README.md b/inputs/capacity_exogenous/README.md index ce7af1a7c..744b0ea8f 100644 --- a/inputs/capacity_exogenous/README.md +++ b/inputs/capacity_exogenous/README.md @@ -8,35 +8,6 @@ - Active when `GSw_NuclearDemo = 1` - For more information, refer to the **notes** section of [https://www.energy.gov/oced/advanced-reactor-demonstration-projects-0](https://www.energy.gov/oced/advanced-reactor-demonstration-projects-0) -- `exog_cap_geohydro_allkm_reference.csv`: - -- `exog_cap_geohydro_reference.csv`: - -- `exog_cap_upv_*.csv`: - - `exog_cap_upv_limited.csv`: - - `exog_cap_upv_open.csv`: - - `exog_cap_upv_reference.csv`: - -- `exog_cap_wind-ons_*.csv`: - - `exog_cap_wind-ons_limited.csv`: - - `exog_cap_wind-ons_open.csv`: - - `exog_cap_wind-ons_reference.csv`: - - `interconnection_queues.csv`: -- `prescribed_builds_wind-ofs_meshed_*.csv`: - - `prescribed_builds_wind-ofs_meshed_limited.csv`: - - `prescribed_builds_wind-ofs_meshed_open.csv`: - - `prescribed_builds_wind-ofs_meshed_reference.csv`: - -- `prescribed_builds_wind-ofs_radial_*.csv`: - - `prescribed_builds_wind-ofs_radial_limited.csv`: - - `prescribed_builds_wind-ofs_radial_open.csv`: - - `prescribed_builds_wind-ofs_radial_reference.csv`: - -- `prescribed_builds_wind-ons_*.csv`: - - `prescribed_builds_wind-ons_limited.csv`: - - `prescribed_builds_wind-ons_open.csv`: - - `prescribed_builds_wind-ons_reference.csv`: - - `ReEDS_generator_database_final_EIA-NEMS.csv`: ReEDS generator database derived from the EIA NEMS plant file and EIA 860M. See https://github.com/ReEDS-Model/ReEDS_Input_Processing/tree/main/nems_database_processing. diff --git a/inputs/supply_curve/dollaryear.csv b/inputs/supply_curve/dollaryear.csv index 1540b9a69..468d6bd91 100644 --- a/inputs/supply_curve/dollaryear.csv +++ b/inputs/supply_curve/dollaryear.csv @@ -6,8 +6,8 @@ evmc_rsc_Baseline,2020 geo_rsc_ATB_2023,2015 hydcost,2004 PSHcostn,2004 -supplycurve_egs,2023 -supplycurve_geohydro,2023 -supplycurve_upv,2023 -supplycurve_wind-ofs,2023 -supplycurve_wind-ons,2022 \ No newline at end of file +supplycurve_init_egs,2023 +supplycurve_init_geohydro,2023 +supplycurve_init_upv,2023 +supplycurve_init_wind-ofs,2023 +supplycurve_init_wind-ons,2022 diff --git a/inputs/userinput/futurefiles.csv b/inputs/userinput/futurefiles.csv index abca66c93..7995ad0d3 100644 --- a/inputs/userinput/futurefiles.csv +++ b/inputs/userinput/futurefiles.csv @@ -428,6 +428,11 @@ storinmaxfrac.csv,.csv,1,None,9999,None,9999,0,constant,None,None,done,constant stressperiods_user.csv,.csv,1,None,9999,None,9999,0,constant,None,None,new; check if we need to expand the input file, superpeak_hour_mapper.csv,.csv,1,None,9999,None,9999,0,constant,None,None,done,constant supply_chain_adjust.csv,.csv,0,None,0,None,0,None,constant,None,None,new, +supplycurve_init_egs.csv,.csv,1,None,9999,None,9999,0,constant,None,None,new, +supplycurve_init_geohydro.csv,.csv,1,None,9999,None,9999,0,constant,None,None,new, +supplycurve_init_upv.csv,.csv,1,None,9999,None,9999,0,constant,None,None,new, +supplycurve_init_wind-ofs.csv,.csv,1,None,9999,None,9999,0,constant,None,None,new, +supplycurve_init_wind-ons.csv,.csv,1,None,9999,None,9999,0,constant,None,None,new, supplycurve_egs.csv,.csv,1,None,9999,None,9999,0,constant,None,None,new, supplycurve_geohydro.csv,.csv,1,None,9999,None,9999,0,constant,None,None,new, supplycurve_upv.csv,.csv,1,None,9999,None,9999,0,constant,None,None,new, diff --git a/inputs/userinput/mcs_distributions_default.yaml b/inputs/userinput/mcs_distributions_default.yaml index 530c9cacb..a5a8d38db 100644 --- a/inputs/userinput/mcs_distributions_default.yaml +++ b/inputs/userinput/mcs_distributions_default.yaml @@ -201,52 +201,52 @@ ### ---------------------------------------------------### ### --- Distributions for NG-fuel prices ### ---------------------------------------------------### -- name: ng_fuel_price +- name: gasprice_flat assignments_list: - - ngscen: [AEO_2025_LOG, AEO_2025_reference, AEO_2025_HOG] + - ngscen: [AEO_2026_HOG, AEO_2026_LOG] dist: dirichlet - dist_params: [1, 1, 1] + dist_params: [1, 1] weight_r: country -- name: ng_fuel_price_tri +- name: gasprice_tri assignments_list: - - ngscen: [AEO_2025_HOG, AEO_2025_reference, AEO_2025_LOG] + - ngscen: [AEO_2026_HOG, AEO_2026_baseline, AEO_2026_LOG] dist: triangular weight_r: country ### ---------------------------------------------------### ### --- Distributions for load ### ---------------------------------------------------### -- name: load_country +- name: load_flat_country assignments_list: - - GSw_LoadProfiles: [EER2025_100by2050, EER2025_IRAlow] + - GSw_LoadProfiles: [EER2025_Baseline_AEO2023, EER2025_100by2050] dist: dirichlet dist_params: [1, 1] weight_r: country +- name: load_tri_country + assignments_list: + - GSw_LoadProfiles: [EER2025_Baseline_AEO2023, EER2025_IRAlow, EER2025_100by2050] + dist: dirichlet + dist_params: [1, 1, 1] + weight_r: country + - name: load_country_unif assignments_list: - GSw_LoadProfiles: [EER2025_IRAlow, EER2025_100by2050] dist: uniform weight_r: country -- name: load_st +- name: load_flat_state assignments_list: - - GSw_LoadProfiles: [EER2025_100by2050, EER2025_IRAlow] + - GSw_LoadProfiles: [EER2025_Baseline_AEO2023, EER2025_100by2050] dist: dirichlet dist_params: [1, 1] weight_r: st -- name: load_ba - assignments_list: - - GSw_LoadProfiles: [EER2025_100by2050, EER2025_IRAlow] - dist: dirichlet - dist_params: [1, 1] - weight_r: ba - - name: load_discrete_st assignments_list: - - GSw_LoadProfiles: [EER2025_100by2050, EER2025_IRAlow] + - GSw_LoadProfiles: [EER2025_IRAlow, EER2025_100by2050] dist: discrete dist_params: [1, 1] weight_r: st @@ -289,7 +289,7 @@ dist_params: [1, 1] weight_r: country -- name: wind_solar_siting_country +- name: siting_lowmid_country assignments_list: - GSw_SitingUPV: [limited, reference] - GSw_SitingWindOfs: [limited, reference] @@ -298,25 +298,25 @@ dist_params: [1, 1] weight_r: country -- name: wind_solar_siting_st +- name: siting_tri_country assignments_list: - - GSw_SitingUPV: [limited, reference] - - GSw_SitingWindOfs: [limited, reference] - - GSw_SitingWindOns: [limited, reference] + - GSw_SitingUPV: [limited, reference, open] + - GSw_SitingWindOfs: [limited, reference, open] + - GSw_SitingWindOns: [limited, reference, open] dist: dirichlet - dist_params: [1, 1] - weight_r: st + dist_params: [1, 1, 1] + weight_r: country -- name: wind_solar_siting_ba +- name: siting_flat_state assignments_list: - - GSw_SitingUPV: [limited, reference] - - GSw_SitingWindOfs: [limited, reference] - - GSw_SitingWindOns: [limited, reference] + - GSw_SitingUPV: [limited, open] + - GSw_SitingWindOfs: [limited, open] + - GSw_SitingWindOns: [limited, open] dist: dirichlet dist_params: [1, 1] - weight_r: ba + weight_r: st -- name: wind_solar_siting_discrete_st +- name: siting_discrete_st assignments_list: - GSw_SitingUPV: [limited, reference] - GSw_SitingWindOfs: [limited, reference] @@ -325,6 +325,15 @@ dist_params: [1, 1] weight_r: st +- name: siting_country_discrete_flat + assignments_list: + - GSw_SitingUPV: [limited, reference, open] + - GSw_SitingWindOfs: [limited, reference, open] + - GSw_SitingWindOns: [limited, reference, open] + dist: discrete + dist_params: [1, 1, 1] + weight_r: country + ### ---------------------------------------------------### ### --- Miscellaneous distributions diff --git a/reeds/input_processing/copy_files.py b/reeds/input_processing/copy_files.py index dc3ee4a1a..de0ae8e5f 100644 --- a/reeds/input_processing/copy_files.py +++ b/reeds/input_processing/copy_files.py @@ -419,6 +419,11 @@ def subset_to_valid_regions( case=os.path.dirname(os.path.normpath(inputs_case)), agg=agg, ).reset_index() + keepcols = [i for i in df if i in [ + 'sc_point_gid', 'class', 'capacity', 'capital_adder_per_mw', + 'cf', 'mean_resource_temp', 'region', + ]] + df = df[keepcols].copy() elif filename == 'techs_banned.csv': df, nuclear_ban_regions = read_banned_tech_file( full_path, diff --git a/reeds/input_processing/mcs_sampler.py b/reeds/input_processing/mcs_sampler.py index 7ccf3b626..844c71aa3 100644 --- a/reeds/input_processing/mcs_sampler.py +++ b/reeds/input_processing/mcs_sampler.py @@ -7,7 +7,6 @@ ### --- IMPORTS --- ### =========================================================================== import argparse -import copy import datetime import numpy as np import os @@ -16,7 +15,6 @@ import sys import yaml from typing import Tuple, List -from collections import defaultdict # Local Imports from pathlib import Path @@ -40,18 +38,21 @@ class MCSConstants: ### --- Fixed columns that should not be modified in most cases OTHER_INDICES = ['columns', 'p', '*p'] # 'p' is used in h2_exog_cap.csv - NONMODIFIABLE_FINANCIAL_COLUMNS = ['debt_fraction', 'tax_rate'] - FIXED_COLUMN_NAMES = YEAR_SYNONYMS + TECH_DESCRIPTOR + OTHER_INDICES + NONMODIFIABLE_FINANCIAL_COLUMNS + REGION_SYNONYMS + OTHER_FIXED_COLUMNS = ['debt_fraction', 'tax_rate', 'class', 'cf'] + FIXED_COLUMN_NAMES = ( + YEAR_SYNONYMS + + TECH_DESCRIPTOR + + OTHER_INDICES + + OTHER_FIXED_COLUMNS + + REGION_SYNONYMS + ) ### --- Files that require special treatment SUPPLY_CURVE_FILES = [ - "supplycurve_upv.csv", - "supplycurve_wind-ofs.csv", - "supplycurve_wind-ons.csv", + "supplycurve_init_upv.csv", + "supplycurve_init_wind-ofs.csv", + "supplycurve_init_wind-ons.csv", ] - EXOG_CAP_FILES = ["exog_cap_upv.csv", "exog_cap_wind-ons.csv"] - PRESCRIBED_BUILDS_FILES = ["prescribed_builds_wind-ofs.csv", "prescribed_builds_wind-ons.csv"] - RECF_FILES = ["recf_wind-ons.h5", "recf_wind-ofs.h5", "recf_upv.h5"] ### --- Switch-File(s) combinations hardcoded in copy_files.py # These files are explicitly handled in copy_files.py, bypassing the standard @@ -411,17 +412,8 @@ def general_mcs_dist_validation(reeds_path: str, mcs_dist_path: str, sw: pd.Seri f"The following switches are not valid (check cases.csv): {invalid_switches}" ) - ## sampling using siting switches is currently disabled - siting_set = set(MCSConstants.SITING_SWITCHES) - used_siting = all_switch_names & siting_set - if used_siting: - raise ValueError( - f"Sampling using siting switches {MCSConstants.SITING_SWITCHES} is " - "currently disabled. For details see " - "https://github.com/ReEDS-Model/ReEDS/issues/41." - ) - ## siting switches can only use specific distributions + siting_set = set(MCSConstants.SITING_SWITCHES) for _, row in df_input_dist.iterrows(): if set(row['switch_names']) & siting_set: if row['dist'] not in ['dirichlet', 'discrete']: @@ -499,8 +491,27 @@ def general_mcs_dist_validation(reeds_path: str, mcs_dist_path: str, sw: pd.Seri f"{sample_group['dist']} for {sample_group['name']} requires " f"{assignment_rule} entries for each switch assignment." ) - - + + +def load_mcs_dist(case) -> pd.DataFrame: + """ + Load the full set of possible Monte Carlo distributions, + downselect to the ones used in the specified case, + and return them as a dataframe. + """ + mcs_dist_path = Path(reeds.io.standardize_case(case), 'inputs_case', 'mcs_distributions.yaml') + with open(mcs_dist_path, 'r') as f: + data = yaml.safe_load(f) + all_options = pd.DataFrame(data) + + sw = reeds.io.get_switches(case) + mcs_dist_groups = sw['MCS_dist_groups'].split('.') + + keep = all_options.loc[all_options['name'].isin(mcs_dist_groups)].reset_index(drop=True) + + return keep + + def get_dist_instructions(reeds_path: str, inputs_case: str) -> Tuple[pd.DataFrame, dict]: """ Obtain the instructions to sample the distributions for each switch @@ -516,17 +527,7 @@ def get_dist_instructions(reeds_path: str, inputs_case: str) -> Tuple[pd.DataFra """ print('Reading the input distribution information for Monte Carlo sampling') - # Read yaml file with the input distribution information. - mcs_dist_path = os.path.join(inputs_case, 'mcs_distributions.yaml') - with open(mcs_dist_path, 'r') as f: - data = yaml.safe_load(f) - df_input_dist = pd.DataFrame(data) - - sw = reeds.io.get_switches(inputs_case) - mcs_dist_groups = sw['MCS_dist_groups'].split('.') - - # Ignore all cases not in mcs_dist_groups - df_input_dist = df_input_dist[df_input_dist['name'].isin(mcs_dist_groups)].reset_index(drop=True) + df_input_dist = load_mcs_dist(inputs_case) # Expand df_input_dist with new information to facilitate the Monte Carlo sampling process. # Sample ID here is used to uniquely identify each sample-process. @@ -700,13 +701,6 @@ def __init__( self.r_weights = get_all_region_weights( self.distribution, self.dist_params, self.hierarchy_file, self.sample_hierarchy_lvl) - # Store the weights for the recf files (CF files) - # Those are computed during the the supply curve file sampling - self.recf_weights_map = {} - - # Flag to validate that recf_weights_map was normalized - self.flag_recf_normalization = defaultdict(lambda: False) - def _validate_inputs(self, dist_files: list, sw_name: str, file_name: str) -> None: """ Validate inputs @@ -727,15 +721,15 @@ def _validate_inputs(self, dist_files: list, sw_name: str, file_name: str) -> No unique_sample_levels = self.hierarchy_file[self.sample_hierarchy_lvl].unique() single_r_weight = len(unique_sample_levels) == 1 - # Group files that require special treatment - except_files = MCSConstants.SUPPLY_CURVE_FILES + MCSConstants.EXOG_CAP_FILES + ( - MCSConstants.PRESCRIBED_BUILDS_FILES + MCSConstants.RECF_FILES) - # Return an error if you have multiple weight assignments but the mcs_distributions.yaml object is # pointing to a set of switches that have no region columns # e.g. asking for a region-based sampling for swicthes.csv, or plantchar type files. - if not single_r_weight and not columns_in_hierarchy and not generic_region_columns and ( - file_name not in except_files): + if ( + (not single_r_weight) + and (not columns_in_hierarchy) + and (not generic_region_columns) + and (file_name not in MCSConstants.SUPPLY_CURVE_FILES) + ): raise ValueError( f"Invalid sampling configuration for file: {file_name}\n" f"Switch: {sw_name}\n" @@ -789,10 +783,6 @@ def get_df_weights( if file_name in MCSConstants.SUPPLY_CURVE_FILES: return self._get_weights_supply_curve(dist_files, modifiable_columns, sw_name) - elif file_name in MCSConstants.RECF_FILES: - return self._get_weights_recf(sw_name) - elif file_name in MCSConstants.EXOG_CAP_FILES + MCSConstants.PRESCRIBED_BUILDS_FILES: - return self._get_weights_exog_prescribed(dist_files) else: return self._get_weights_general(dist_files, modifiable_columns, sw_name, file_name) @@ -901,21 +891,17 @@ def _get_weights_supply_curve( # Dictionary to store computed weights for the modifiable columns # file index -> pd.DataFrame dict_df_weights = {} + base_scenario = reeds.io.get_siting_switchval(switchname=sw_name, case=inputs_case) + siting_scenarios = { + key: val for entry in self.sample_group.assignments_list + for key, val in entry.items() + }[sw_name] + base_scenario_index = siting_scenarios.index(base_scenario) + dfbase = dist_files[base_scenario_index].set_index('sc_point_gid') - # Store weights to use later in the recf files (CF files) - self.recf_weights_map[sw_name] = {} - - # Create a new column with the class|region combination (like in the CF file) - dist_files_copy = [copy.deepcopy(df) for df in dist_files] - - for df in dist_files_copy: - df["old c|r"] = ( - df["class"].astype(int).astype(str) + "|" + df["region"].astype(str) - ) - - for f, df in enumerate(dist_files_copy): + for f, df in enumerate(dist_files): # Initial skeleton of the weights DataFrame - w_df_tmp = df[["region", "sc_point_gid", "old c|r"]] + w_df_tmp = df[["region", "sc_point_gid"]].copy() # Create a mapping from each unique region to its corresponding weight region_to_weight = { @@ -936,131 +922,23 @@ def _get_weights_supply_curve( # Join the modifiable columns back into the original DataFrame w_df_tmp = w_df_tmp.join(modifiable_df) - # Save the intermediate weights for the recf files (These are weights multiplied by capacity) - self.recf_weights_map[sw_name][f] = w_df_tmp[["old c|r","class"]].rename(columns={"class": "weight"}) - # Store in dictionary - dict_df_weights[f] = w_df_tmp.drop(columns=["old c|r"]) - - # Normalize the weights to sum to 1 - # Divide the weights by the sum of the weights across all files - sum_weights = sum(dict_df_weights[f][modifiable_columns] for f in range(len(dist_files))) - sum_weights[sum_weights == 0] = 1 - - for f in range(len(dist_files)): - dict_df_weights[f][modifiable_columns] /= sum_weights - # recf_weights_map is not normalized here because it will be normalized later - # values need to be aggregated according to the new c|r column from supply curves - - return dict_df_weights - - def _get_weights_exog_prescribed(self, dist_files: list) -> dict: - """ - Get the weights for exogenous capacity and prescribed builds files. - - Args: - dist_files (list of pd.DataFrame): List of reference dataframes (ajusted to have the same # of rows) - - Returns: - Dict[int, pd.DataFrame]: Dictionary mapping reference file index to - the weight DataFrame for that file. - """ - dict_df_weights = {} - for f, df in enumerate(dist_files): + dict_df_weights[f] = w_df_tmp - region_to_weight = { - r: self.r_weights[r][f] - for r in df["region"].unique() - } + # Normalize the weights to sum to 1 + # Divide the weights by the sum of the weights across all files + sum_weights = sum(dict_df_weights[f][modifiable_columns] for f in range(len(dist_files))) + sum_weights[sum_weights == 0] = 1 - dict_df_weights[f] = pd.DataFrame( - data=df["region"].map(region_to_weight).values, - columns=["capacity"], - index=df.index, - ) + for f in range(len(dist_files)): + dict_df_weights[f][modifiable_columns] /= sum_weights + ## Inherit everything but capacity from the base scenario + for col in ['class','cf','capital_adder_per_mw']: + dict_df_weights[f][col] = 1 if f == base_scenario_index else 0 + dist_files[f][col] = dist_files[f].sc_point_gid.map(dfbase[col]) return dict_df_weights - def _get_weights_recf(self, sw_name: str) -> dict: - """ - Get the weights for the recf files (CF files). This file construction is - dependent on the new supply curve samples and therefore is computed after - the supply curve files are sampled. - - Args: - sw_name (str): Name of the switch we are getting the weights for. - """ - # From get_dist_instructions(.) the supply curve file is deliberaly - # placed before the recf files, so that recf_weights_map is already populated. - - # Check if recf_weights_map is not empty and that it was normalized. - if not self.recf_weights_map[sw_name]: - raise ValueError( - f"The recf_weights_map for switch {sw_name} was not populated" - ) - - if not self.flag_recf_normalization[sw_name] : - raise ValueError( - f"The recf_weights_map for switch {sw_name} was not normalized" - ) - - return self.recf_weights_map[sw_name] - - def normalize_recf_weights_map(self, samples_sw: pd.DataFrame, sw_name: str) -> None: - """ - The recf map is responsible for informing how the old class/region data files - need to be put together (weights) to form the new class/region data. - After creating the new supply curve sample, we normalize the weights to sum to 1. - - Args: - samples_sw (pd.DataFrame): The sampled supply curve DataFrame. - sw_name (str): Name of the switch being sampled. - - Updates: - self.recf_weights_map (dict): Dictionary with the normalized weights for the recf files. - Each element of this dictionary is a pd.DataFrame (for the reference file f) - with the normalized weights, indexed by new and old class|region (c|r). - """ - n_files = len({key[1] for key in self.recf_weights_map[sw_name].keys()}) - - # The normalization can change depending on the sample # - for f in range(n_files): - # Add a new column with the new class|region combination - self.recf_weights_map[sw_name][f]["new c|r"] = ( - samples_sw["class"].astype(str) + "|" + - samples_sw["region"].astype(str) - ) - - # Sum the weights for each new class|region combination - self.recf_weights_map[sw_name][f] = self.recf_weights_map[sw_name][f].groupby( - ["new c|r","old c|r"], as_index=False).sum() - - # Remove cases with 0 weight (e.g old c|r had no capacity -> class 0) - self.recf_weights_map[sw_name][f] = self.recf_weights_map[sw_name][f][ - self.recf_weights_map[sw_name][f]["weight"] > 0 - ] - - # Go over all files and obtain the total sum of weights for each new class|region - sum_weights_recf_map = ( - pd.concat( - [self.recf_weights_map[sw_name][f] for f in range(n_files)] - ) - .groupby("new c|r")["weight"] - .sum() - .to_dict() - ) - - for f in range(n_files): - # Get current DataFrame - df = self.recf_weights_map[sw_name][f] - # Perform division using "new c|r" as the reference - df["weight"] = df["weight"] / df["new c|r"].map(sum_weights_recf_map) - # Assign back to original structure - self.recf_weights_map[sw_name][f] = df.set_index(["new c|r", "old c|r"]) - - # Flag to validade that recf_weights_map was normalized - self.flag_recf_normalization[sw_name] = True - #%% =========================================================================== ### --- MAIN SAMPLING CLASS --- @@ -1133,9 +1011,7 @@ def prepare_ref_data( # (e.g For the supply curves we will make sure that all files are # ajusted to contain all regions and sc_point_gid combinations) map_files2ref_columns = { - **{file: ["region", "sc_point_gid"] for file in MCSConstants.SUPPLY_CURVE_FILES}, - **{file: ["region", "year", "sc_point_gid"] for file in MCSConstants.EXOG_CAP_FILES}, - **{file: ["region", "year"] for file in MCSConstants.PRESCRIBED_BUILDS_FILES}, + file: ["region", "sc_point_gid"] for file in MCSConstants.SUPPLY_CURVE_FILES } if file_name in map_files2ref_columns: @@ -1165,25 +1041,23 @@ def prepare_ref_data( } # Ensure all dist_files have the same set of general columns - if file_name not in MCSConstants.RECF_FILES: - for i, df in enumerate(dist_files[1:], start=1): - current_cols = {col for col in df.keys() if col not in MCSConstants.FIXED_COLUMN_NAMES} - if current_cols != general_mult_columns: - error_msg = ( - f"Column mismatch between dist_files[0] and dist_files[i]:\n" - "This usually happens when you run MCS on a file whose columns " - "vary by switch assignment (e.g. RECF_FILES).\n If you really need to support " - f"'{file_name}' here, add the necessary handling in prepare_ref_data()." - ) - raise ValueError(error_msg) + for i, df in enumerate(dist_files[1:], start=1): + current_cols = {col for col in df.keys() if col not in MCSConstants.FIXED_COLUMN_NAMES} + if current_cols != general_mult_columns: + error_msg = ( + f"Column mismatch between dist_files[0] and dist_files[i]:\n" + "This usually happens when you run MCS on a file whose columns " + "vary by switch assignment.\n If you really need to support " + f"'{file_name}' here, add the necessary handling in prepare_ref_data()." + ) + raise ValueError(error_msg) - exceptions_mult_col = { - **{file: ["class"] + list(general_mult_columns) for file in MCSConstants.SUPPLY_CURVE_FILES}, - **{file: ["capacity"] for file in MCSConstants.EXOG_CAP_FILES}, - **{file: ["capacity"] for file in MCSConstants.PRESCRIBED_BUILDS_FILES}, - **{file: [] for file in MCSConstants.RECF_FILES}, # treated separately - } - modifiable_columns = exceptions_mult_col.get(file_name, list(general_mult_columns)) + ## Supply curve files are special; capacity is sampled, but all other parameters + ## use values from the most-permissive scenario (the same as hourly profiles) + if file_name in MCSConstants.SUPPLY_CURVE_FILES: + modifiable_columns = ['capacity'] + else: + modifiable_columns = list(general_mult_columns) ### =========================================================================== ### --- Map for the number of decimals in each column we will change @@ -1257,51 +1131,9 @@ def _adjust_supply_curve_sample(self, samples_sw: pd.DataFrame, sw_name: str, sa # Convert class to integer samples_sw["class"] = samples_sw["class"].astype(int) - # Update the recf weights map (weight_calc.recf_weights_map) - self.weight_calc.normalize_recf_weights_map(samples_sw, sw_name) - - # Remove samples with no capacity. - # Need to do this after normalizing the recf weights - samples_sw = samples_sw[samples_sw["capacity"] > 0] - - return samples_sw - - def _adjust_exog_cap_samples(self, samples_sw: pd.DataFrame, file_name: str) -> pd.DataFrame: - """ - Adjust samples for exogenous capacity files: - - Remove rows with no capacity. - - Adjust the tech classes based on available classes per sc_point_gid. - - Args: - samples_sw (pd.DataFrame): The sampled exogenous capacity DataFrame. - file_name (str): Name of the file being sampled. - - Returns: - pd.DataFrame: Adjusted exogenous capacity sample. - """ - # Remove samples with no capacity - samples_sw = samples_sw[samples_sw["capacity"] > 0].copy() - - tech_mapping = { - "exog_cap_upv.csv": ("upv", "supplycurve_upv.csv"), - "exog_cap_wind-ons.csv": ("wind-ons", "supplycurve_wind-ons.csv"), - } - tech_name, Sample_ID = tech_mapping[file_name] - - # Get the class available for each sc_point_gid - class_sc_point_map = self.samples[Sample_ID][["sc_point_gid", "class"]] - class_sc_point_map = class_sc_point_map.set_index("sc_point_gid").to_dict()["class"] - - # Remove any rows from samples_sw that cannot be mapped - # These are cases with zero supply in the region - valid_sc_point_gids = samples_sw["sc_point_gid"].isin(class_sc_point_map.keys()) - samples_sw = samples_sw[valid_sc_point_gids].copy() - - # Create a new tech name for each sc_point_gid - new_tech_name = [tech_name + "_" + str(int(c)) for c in - samples_sw["sc_point_gid"].map(class_sc_point_map).values] - - samples_sw["*tech"] = new_tech_name + ## Remove samples with no capacity or class + ## (meaning they are not in the access scenario used for hourly profiles) + samples_sw = samples_sw.loc[(samples_sw["capacity"] > 0) & (samples_sw['class'] > 0)].copy() return samples_sw @@ -1315,7 +1147,7 @@ def _apply_weights_general( ): """ Apply the distribution weights to the reference files. - Applicable to all cases but recf files and switches.csv. + Applicable to all cases but switches.csv. Args: dist_files (List[pd.DataFrame]): List of input DataFrames for sampling. @@ -1344,14 +1176,6 @@ def _apply_weights_general( if file_name in MCSConstants.SUPPLY_CURVE_FILES: adjusted_samples = self._adjust_supply_curve_sample(samples_sw, sw_name, sample_idx) - - elif file_name in MCSConstants.EXOG_CAP_FILES: - adjusted_samples = self._adjust_exog_cap_samples(samples_sw, file_name) - - elif file_name in MCSConstants.PRESCRIBED_BUILDS_FILES: - # Remove samples with no capacity - adjusted_samples = samples_sw[samples_sw["capacity"] > 0] - else: # For all other files we can directly apply the weights adjusted_samples = samples_sw @@ -1359,45 +1183,6 @@ def _apply_weights_general( # Save the adjusted samples. self.samples[Sample_ID] = adjusted_samples - def _apply_weights_recf( - self, - dist_files: list, - sample_idx: int - ): - """ - Apply the distribution weights to the recf files. - This file gets compleatly overwriten so need to be treated separately - - Args: - dist_files (List[pd.DataFrame]): List of input DataFrames for sampling. - sample_idx (int): Index of the Sample_ID in sample_group. - - Update: - self.samples (Dict[str, pd.DataFrame]): Dictionary with the samples for each switch/file_name. - """ - - Sample_ID = self.sample_group['Sample_ID'][sample_idx] - sw_name = self.sample_group['switch_names'][sample_idx] - - # For the recf files we need to apply the weights to the old class|region combinations - weights = self.weight_calc.recf_weights_map[sw_name] - # Index is the same for all files (time) - indexes = dist_files[0].index - - # get initial switch values - sample_sw = defaultdict(int) - - for f, df in enumerate(dist_files): - # Get the old and new class|region combinations from weights[(s, f)] - for (new_c_r, old_c_r) in weights[f].index: - sample_sw[new_c_r] += df[old_c_r] * weights[f].loc[(new_c_r,old_c_r)].values[0] - - # Round numbers to 9 decimal places and allow min/max values of 0 and 1 - for new_c_r in sample_sw.keys(): - sample_sw[new_c_r] = sample_sw[new_c_r].round(9).clip(0,1) - - self.samples[Sample_ID] = pd.DataFrame(sample_sw, index=indexes) - def _apply_weights_switches_csv( self, dist_files: list, @@ -1629,7 +1414,7 @@ def apply_lhs_switches_csv(self, sample_group_num, sample_idx, dist_files, n_dec def apply_lhs_general(self, sample_group_num, Sample_ID, dist_files, aux_files, modifiable_columns): - """Apply LHS sampling to a general (non-switch, non-RECF) file. + """Apply LHS sampling to a general (non-switch) file. For each modifiable column, draws values from the configured distribution using the LHS quantile matrix and stores the resulting DataFrame in self.samples. @@ -1723,8 +1508,6 @@ def get_samples(self, aux_files, sample_group_num): # apply latin hybercube sampling based on file type if file_name == "switches.csv": self.apply_lhs_switches_csv(sample_group_num, sample_idx, dist_files, n_decimals) - elif file_name in MCSConstants.RECF_FILES: - self._apply_weights_recf(dist_files, sample_idx) else: self.apply_lhs_general(sample_group_num, Sample_ID, dist_files, aux_files, modifiable_columns) else: @@ -1732,8 +1515,6 @@ def get_samples(self, aux_files, sample_group_num): # Dispatch weight application based on file type if file_name == "switches.csv": self._apply_weights_switches_csv(dist_files, n_decimals, dict_df_weights, sample_idx) - elif file_name in MCSConstants.RECF_FILES: - self._apply_weights_recf(dist_files, sample_idx) else: self._apply_weights_general(dist_files, modifiable_columns, n_decimals, dict_df_weights, sample_idx) @@ -2009,8 +1790,9 @@ def main_mga_rv( mga_weights = mga_weights.rename(columns={'r':'*r'})[['*r','i_subtech','weight']] mga_weights = mga_weights.sort_values(by=['*r', 'i_subtech'], ascending=True) mga_weights.to_csv(os.path.join(inputs_case, "mga_weights.csv"), index=False) - + +#%% Procedure if __name__ == '__main__' and not hasattr(sys, 'ps1'): parser = argparse.ArgumentParser(description='Copy files needed for this run') parser.add_argument('reeds_path', help='ReEDS directory') @@ -2024,7 +1806,7 @@ def main_mga_rv( # ---- Settings for testing ---- # reeds_path = reeds.io.reeds_path - # inputs_case = os.path.join(reeds_path,'runs','v20250825_revM2_MonteCarlo_MC1','inputs_case') + # inputs_case = os.path.join(reeds_path,'runs','v20261008_mcM0_MC_tri_country_MC0001','inputs_case') # n_samples = 1 # seed = 0 @@ -2065,4 +1847,3 @@ def main_mga_rv( process='input_processing/mcs_sampler.py', path=os.path.join(os.path.dirname(inputs_case)) ) - diff --git a/reeds/input_processing/process_unitdata.py b/reeds/input_processing/process_unitdata.py index f21a8e348..69c6b22af 100644 --- a/reeds/input_processing/process_unitdata.py +++ b/reeds/input_processing/process_unitdata.py @@ -39,12 +39,8 @@ def assign_gids_to_unitdata(sw, df, offland_gdf, land_gdf): df_sub = df[df.tech.isin(tech_sub)] # Read supply curves - if tech == 'geohydro': - # Use egs supply curve for geohydro for now - geo_tech = 'egs' - supply_curve = pd.read_csv(os.path.join(inputs_case,'supplycurve_'+geo_tech+'.csv')) - else: - supply_curve = pd.read_csv(os.path.join(inputs_case,'supplycurve_'+tech+'.csv')) + usetech = 'egs' if tech == 'geohydro' else tech + supply_curve = pd.read_csv(Path(inputs_case, f'supplycurve_init_{usetech}.csv')) # Only consider the sc_point_gids that are in supply curves # (to avoid unmatched units later) diff --git a/reeds/input_processing/recf.py b/reeds/input_processing/recf.py index e1c768fcf..19532cac8 100644 --- a/reeds/input_processing/recf.py +++ b/reeds/input_processing/recf.py @@ -159,7 +159,7 @@ def calculate_class_region_cf_hourly( sw = reeds.io.get_switches(inputs_case) # Get supply curve information df_sc = reeds.io.assemble_supplycurve( - os.path.join(inputs_case, f'supplycurve_{tech}.csv'), + os.path.join(inputs_case, f'supplycurve_init_{tech}.csv'), case=os.path.dirname(inputs_case), agg=True, ) @@ -189,8 +189,23 @@ def calculate_class_region_cf_hourly( year=year, case=inputs_case, ) - # Downselect to relevant sites - weather_year_site_cf_hourly = weather_year_site_cf_hourly[df_sc.index] + # Downselect to relevant sites. + # For scenarios not using Monte Carlo, we require an exact match to make + # sure our datasets stay in sync. + # For Monte Carlo scenarios, we allow for missing values, because the + # limited scenario contains sites not in reference/open, but we only use profiles + # from the most-open scenario in the distribution. + if not int(sw['MCS_runs']): + weather_year_site_cf_hourly = weather_year_site_cf_hourly[df_sc.index] + else: + keep_sites = [i for i in weather_year_site_cf_hourly if i in df_sc.index] + if len(keep_sites) != len(df_sc): + dropped = [i for i in df_sc.index if i not in keep_sites] + print( + f'MCS: Dropped {len(dropped)} {tech} sites: ' + + ', '.join([str(i) for i in dropped]) + ) + weather_year_site_cf_hourly = weather_year_site_cf_hourly[keep_sites] # Calculate the capacity-weighted average CF for each class-region pair weather_year_class_region_cf_hourly = ( weather_year_site_cf_hourly.mul(df_sc['capacity']) @@ -582,7 +597,7 @@ def main(reeds_path, inputs_case): # #%% Settings for testing # reeds_path = reeds.io.reeds_path - # inputs_case = os.path.join(reeds_path,'runs','v20260804_inputsM0_MARICTNYNJPAOH_Offshore','inputs_case') + # inputs_case = os.path.join(reeds_path,'runs','v20261006_mcM1_MonteCarlo_Random_MC0001','inputs_case') log = reeds.log.makelog( scriptname=__file__, diff --git a/reeds/input_processing/runfiles.csv b/reeds/input_processing/runfiles.csv index e22922330..3b5695ded 100644 --- a/reeds/input_processing/runfiles.csv +++ b/reeds/input_processing/runfiles.csv @@ -223,10 +223,10 @@ storage_duration.csv,inputs/storage/storage_duration.csv,int(sw.GSw_Storage) != storage_mandates.csv,inputs/state_policies/storage_mandates.csv,int(sw.GSw_BatteryMandate) != 0,ignore,ignore,*st,t,,0,0,,,,,, stressperiods_user.csv,inputs/temporal/stressperiods_{GSw_PRM_StressModel}.csv,sw.GSw_PRM_StressModel == 'user*',ignore,ignore,,,,,,,,,,, supply_chain_adjust.csv,inputs/financials/supply_chain_adjust.csv,1,ignore,ignore,,,,,0,,,,,, -supplycurve_egs.csv,inputs/supply_curve/supplycurve_egs-reference.csv,int(sw.GSw_Geothermal) != 0,ignore,ignore,region,,,,0,,,,,,No fix_col (only region_col) because this is an intermediate input file -supplycurve_upv.csv,inputs/supply_curve/supplycurve_upv-{GSw_SitingUPV}.csv,1,ignore,ignore,region,,,,0,,,,,,No fix_col (only region_col) because this is an intermediate input file -supplycurve_wind-ofs.csv,inputs/supply_curve/supplycurve_wind-ofs-{GSw_SitingWindOfs}.csv,int(sw.GSw_OfsWind) != 0,ignore,ignore,region,,,,0,,,,,,No fix_col (only region_col) because this is an intermediate input file -supplycurve_wind-ons.csv,inputs/supply_curve/supplycurve_wind-ons-{GSw_SitingWindOns}.csv,1,ignore,ignore,region,,,,0,,,,,,No fix_col (only region_col) because this is an intermediate input file +supplycurve_init_egs.csv,inputs/supply_curve/supplycurve_egs-reference.csv,int(sw.GSw_Geothermal) != 0,ignore,ignore,region,,,,0,,,,,,No fix_col (only region_col) because this is an intermediate input file +supplycurve_init_upv.csv,inputs/supply_curve/supplycurve_upv-{GSw_SitingUPV}.csv,1,ignore,ignore,region,,,,0,,,,,,No fix_col (only region_col) because this is an intermediate input file +supplycurve_init_wind-ofs.csv,inputs/supply_curve/supplycurve_wind-ofs-{GSw_SitingWindOfs}.csv,int(sw.GSw_OfsWind) != 0,ignore,ignore,region,,,,0,,,,,,No fix_col (only region_col) because this is an intermediate input file +supplycurve_init_wind-ons.csv,inputs/supply_curve/supplycurve_wind-ons-{GSw_SitingWindOns}.csv,1,ignore,ignore,region,,,,0,,,,,,No fix_col (only region_col) because this is an intermediate input file wst_surface.csv,inputs/sets/wst_surface.csv,1,ignore,ignore,,,,,0,,,set,wst_surface,surface water types where access is based on consumption not withdrawal, tc_phaseout_schedule.csv,inputs/financials/tc_phaseout_schedule_{GSw_TCPhaseout_schedule}.csv,int(sw.GSw_TCPhaseout) != 0,ignore,ignore,,,,,0,,,,,, tech-subset-table.csv,inputs/tech-subset-table.csv,1,ignore,ignore,,,,,0,,,,,, diff --git a/reeds/input_processing/writecapdat.py b/reeds/input_processing/writecapdat.py index a084c2cfb..3e4e3c264 100644 --- a/reeds/input_processing/writecapdat.py +++ b/reeds/input_processing/writecapdat.py @@ -31,16 +31,135 @@ import sys import pandas as pd from pathlib import Path +from typing import Literal sys.path.append(str(Path(__file__).parent.parent.parent)) import reeds #%%################# ### FIXED INPUTS ### -WINDOFS_FIXED_CLASSES = list(range(5)) -WINDOFS_FLOATING_CLASSES = list(range(6,11)) +WINDOFS_CLASS_TYPE = { + **{i: 'fixed' for i in range(1,6)}, + **{i: 'floating' for i in range(6,11)}, +} ONLINEYEAR_COLUMNS = ['i', 'r', 'StartYear', 'RetireYear', 'summer_power_capacity_MW'] +# Only keep neccessary columns from unitdata to work with +# And rename column names for easier processing +COLNAMES = { + 'capexog_rsc': ( + ['tech','r','RetireYear','StartYear','sc_point_gid','summer_power_capacity_MW'], + ['tech','region','year','onlineyear','sc_point_gid','MW'] + ), + 'capnonrsc': ( + ['tech','coolingwatertech','r','ctt','wst','summer_power_capacity_MW'], + ['i','coolingwatertech','r','ctt','wst','value'] + ), + 'capnonrsc_energy': ( + ['tech','r','energy_capacity_MWh'], + ['i','r','value'] + ), + 'prescribed_nonRSC': ( + ['StartYear','tech','vin','r','coolingwatertech','ctt','wst','summer_power_capacity_MW'], + ['t','i','v','r','coolingwatertech','ctt','wst','value'] + ), + 'prescribed_nonRSC_energy': ( + ['StartYear','tech','vin','r','coolingwatertech','ctt','wst','energy_capacity_MWh'], + ['t','i','v','r','coolingwatertech','ctt','wst','value'] + ), + 'prescribed_RSC': ( + ['StartYear','tech','vin','r','summer_power_capacity_MW'], + ['t','i','v','r','value'] + ), + 'rsc': ( + ['tech','r','v','ctt','wst','summer_power_capacity_MW'], + ['i','r','v','ctt','wst','value'] + ), + 'rsc_wsc': ( + ['r','tech','summer_power_capacity_MW'], + ['r','i','value'] + ), + 'prsc_csp': ( + ['StartYear','r','tech','ctt','wst','summer_power_capacity_MW'], + ['t','r','i','ctt','wst','value'] + ), + 'prsc_geo': ( + ['StartYear','r','tech','summer_power_capacity_MW'], + ['t','r','i','value'] + ), + 'retirements': ( + ['tech','v','r','RetireYear','StartYear','coolingwatertech','ctt','wst','type','summer_power_capacity_MW'], + ['i','v','r','t','tt','coolingwatertech','ctt','wst','type','value'] + ), + 'retirements_energy': ( + ['tech','v','r','RetireYear','StartYear','type','energy_capacity_MWh'], + ['r','i','v','t','tt','type','value'] + ), + 'windret': ( + ['r','tech','RetireYear','summer_power_capacity_MW'], + ['r','i','t','value'] + ), + 'georet': ( + ['r','tech','RetireYear','summer_power_capacity_MW'], + ['r','i','t','value'] + ), +} + +TECH = { + 'capnonrsc': [ + 'battery_li', 'biopower', 'coal-igcc', 'coal-new', + 'coaloldscr','coalolduns','gas-cc', 'gas-ct', + 'lfill-gas','nuclear', 'o-g-s', 'pumped-hydro' + ], + 'capnonrsc_energy': [ + 'battery_li' + ], + 'prescribed_nonRSC': [ + 'battery_li', 'biopower', 'coal-igcc', 'coal-new', + 'coaloldscr', 'coalolduns', 'gas-cc', 'gas-ct', + 'hydED', 'hydEND', 'hydUD', 'hydUND', 'hydND', 'hydNPND', + 'lfill-gas', 'nuclear', 'o-g-s', 'pumped-hydro' + ], + 'prescribed_nonRSC_energy': [ + 'battery_li', + ], + 'storage' : ['battery_li', 'pumped-hydro' + ], + 'rsc_pv_all': ['upv','pvb','pvb_pv','csp-ns'], + 'rsc_upv': ['upv','pvb'], + 'rsc_w': ['wind-ons','wind-ofs'], + 'rsc_csp': ['csp-ns'], + 'rsc_wsc': ['upv','pvb','csp-ns','csp-ws','wind-ons','wind-ofs', + 'geohydro_allkm','egs_allkm'], + 'prsc_csp': ['csp-ns','csp-ws'], + 'prsc_geo': ['geohydro_allkm','egs_allkm'], + 'retirements': [ + 'coalolduns', 'o-g-s', 'hydED', 'hydEND', 'gas-ct', 'lfill-gas', + 'coaloldscr', 'biopower', 'gas-cc', 'coal-new', + 'battery_li','nuclear', 'pumped-hydro', 'coal-igcc' + ], + 'retirements_energy': [ + 'battery_li' + ], + 'windret': ['wind-ons'], + 'georet': ['geohydro_allkm','egs_allkm'], + # This is not all technologies that do not having cooling, but technologies + # that are (or could be) in the plant database. + 'no_cooling': [ + 'upv', 'pvb', 'gas-ct', 'geohydro_allkm','egs_allkm', + 'battery_li', 'pumped-hydro', 'pumped-hydro-flex', + 'hydUD', 'hydUND', 'hydD', 'hydND', 'hydSD', 'hydSND', 'hydNPD', + 'hydNPND', 'hydED', 'hydEND', 'wind-ons', 'wind-ofs', + ], +} + +quartershorten = {'spring':'spri', 'summer':'summ', 'fall':'fall', 'winter':'wint'} + +hotcold_months = { + 'NOV':'cold', 'DEC':'cold', 'JAN':'cold', 'FEB':'cold', + 'JUN':'hot', 'JUL':'hot', 'AUG':'hot' +} + #%% =========================================================================== ### --- FUNCTIONS --- ### =========================================================================== @@ -56,63 +175,104 @@ def create_rsc_wsc(gendb,TECH,startyear): return rsc_wsc -def create_exog_rsc(reeds_path,inputs_case,gendb,TECH,COLNAMES,sw,startyear): - # Mappings to resource class are based on the resource quality of the technology as it comes from reV - # Establish resource classification inputs for technologies (UPV, wind-ons, wind-ofs) - # from supply curves - - - rsc_class = {} - rsc_class["upv"] = get_class_cf_bounds(reeds_path, tech='upv', - access_case=sw.GSw_SitingUPV, subtech='') - rsc_class["wind-ons"] = get_class_cf_bounds(reeds_path, tech='wind-ons', - access_case=sw.GSw_SitingWindOns, subtech='') - - # for offshore wind, specify 'fixed' or 'floating' tech - wind_ofs_subtech_list = ['fixed','floating'] - wind_ofs_class_all = [] - for wind_ofs_subtech in wind_ofs_subtech_list: - wind_ofs_class_subtech = get_class_cf_bounds(reeds_path, tech='wind-ofs', - access_case=sw.GSw_SitingWindOfs,subtech=wind_ofs_subtech) - wind_ofs_class_all.append(wind_ofs_class_subtech) - rsc_class["wind-ofs"] = pd.concat(wind_ofs_class_all, ignore_index=True) - - # Read resource classification inputs for geothermal - rsc_class["geohydro_allkm"] = ( - pd.read_csv(os.path.join(inputs_case, 'classification_geothermal.csv')) - .query(f"access_case == '{sw.GSw_SitingGeo}'") + +def get_class_cf_min( + inputs_case, + tech:Literal['upv','wind-ons','wind-ofs'], + subtech:None|Literal['fixed','floating']=None, +) -> pd.DataFrame: + """ + Establish class cut offs based on capacity factors. + Because the max of one class is the min of the next class, we only need the min. + """ + ## Get the supply curve for this run + dfsc = pd.read_csv(Path(inputs_case, f'supplycurve_init_{tech}.csv')) + + ## Downselect offshore wind to fixed or floating + if tech == 'wind-ofs': + dfsc['subtech'] = dfsc['class'].map(WINDOFS_CLASS_TYPE) + dfsc = dfsc.loc[dfsc['subtech'] == subtech] + + class_cf_min = dfsc.groupby('class')['cf'].min() + return class_cf_min + + +def get_class_cf_min_alltechs(inputs_case) -> dict: + """ + Generate a dictionary of dataframes for supply-curve-based generation technologies + indicating the CF (or temperature for geothermal) range for each resource class + + Returns: + dict: + keys are ['upv', 'wind-ons', 'wind-ofs', 'geohydro_allkm'] + values are dataframes mapping from a resource class to a CF range + """ + sw = reeds.io.get_switches(inputs_case) + tech_class_cutoffs = {} + for tech in ['upv', 'wind-ons']: + tech_class_cutoffs[tech] = get_class_cf_min(inputs_case, tech) + + ## Offshore wind is separated into fixed and floating subtechs + for subtech in ['fixed', 'floating']: + tech_class_cutoffs['wind-ofs', subtech] = get_class_cf_min(inputs_case, 'wind-ofs', subtech) + + ## Geothermal is hardcoded, and the classes are ordered in the opposite direction + ## from wind/solar (the lowest geo class is the highest temperature, while the + ## lowest wind/solar class is the lowest CF) + fpath_geo = Path(inputs_case, 'classification_geothermal.csv') + tech_class_cutoffs["geohydro_allkm"] = ( + pd.read_csv(fpath_geo, index_col='access_case') + .loc[sw.GSw_SitingGeo] + .set_index('class').min_reV_mean_temp.sort_values(ascending=True) ) - # Check if any rsc_wsc tech class in unitdata does not match with a resource class - missing_resource_class(gendb,rsc_class) + return tech_class_cutoffs + + +def assign_class(cf:float, tech_class_cutoff:pd.Series): + """Assign a resource class based on a capacity factor""" + ## Work down from the highest-CF class and keep the first match + for resource_class, value in tech_class_cutoff.iloc[::-1].items(): + if cf >= value: + return resource_class + ## Only end up here if the provided cf is outside the range of the lookup table + err = f'cf={cf} is not in the provided tech_class_cutoff table:\n{tech_class_cutoff}' + raise ValueError(err) + + +def create_exog_rsc(inputs_case, gdb_use_cap_exog, sw, startyear): + tech_class_cutoffs = get_class_cf_min_alltechs(inputs_case) + colname = { + 'geohydro_allkm':'reV_mean_resource_temp', + 'egs_allkm':'reV_mean_resource_temp', + } + colname_default = 'reV_capacity_factor_ac' cap_exog = {} for tech in TECH['rsc_wsc']: - print(tech) + ## All solar techs use UPV classes; all exogenous offshore wind is fixed (not floating) + match tech: + case _ if tech in TECH['rsc_pv_all']: + usetech = 'upv' + case 'wind-ofs': + usetech = ('wind-ofs', 'fixed') + case _: + usetech = tech + print(tech, usetech) # Filter active plants - cap_exog[tech]= gendb.loc[(gendb['tech']==tech) & - (gendb['StartYear'] < startyear) & - (gendb['RetireYear'] > startyear)].copy() + cap_exog[tech]= gdb_use_cap_exog.loc[ + (gdb_use_cap_exog['tech']==tech) + & (gdb_use_cap_exog['StartYear'] < startyear) + & (gdb_use_cap_exog['RetireYear'] > startyear) + ].copy() if len(cap_exog[tech]) > 0: - # Assigning each geothermal unit in unit database to a class based on - # groups' temperatures - if tech in ['geohydro_allkm','egs_allkm']: - cap_exog[tech]["class"] = cap_exog[tech]["reV_mean_resource_temp"].apply( - lambda x: assign_class(x, tech, rsc_class[tech])) - cap_exog[tech]["tech"] = (cap_exog[tech]["tech"].astype(str) + "_" + - cap_exog[tech]["class"].astype(str)) - # Assigning each solar, wind unit in unit database to a class based on - # groups' minimum and maximum capacity factors - elif tech in TECH['rsc_pv_all']: - cap_exog[tech]["class"] = cap_exog[tech]["reV_capacity_factor_ac"].apply( - lambda x: assign_class(x, tech, rsc_class['upv'])) - cap_exog[tech]["tech"] = ('upv' + "_" + - cap_exog[tech]["class"].astype(str)) - else: - cap_exog[tech]["class"] = cap_exog[tech]["reV_capacity_factor_ac"].apply( - lambda x: assign_class(x, tech, rsc_class[tech])) - cap_exog[tech]["tech"] = (cap_exog[tech]["tech"].astype(str) + "_" + - cap_exog[tech]["class"].astype(str)) + cap_exog[tech]["class"] = cap_exog[tech][colname.get(tech,colname_default)].apply( + lambda x: assign_class(x, tech_class_cutoffs[usetech]) + ) + cap_exog[tech]["tech"] = ( + ('upv' if usetech == 'upv' else cap_exog[tech]["tech"].astype(str)) + + "_" + cap_exog[tech]["class"].astype(str) + ) cap_exog[tech] = cap_exog[tech][COLNAMES['capexog_rsc'][0]] cap_exog[tech].columns = COLNAMES['capexog_rsc'][1] @@ -120,69 +280,8 @@ def create_exog_rsc(reeds_path,inputs_case,gendb,TECH,COLNAMES,sw,startyear): cap_exog[tech] = pd.concat([expand_exog_cap(row, startyear) for _, row in cap_exog[tech].iterrows()], ignore_index=True) - return cap_exog, rsc_class - -def get_class_cf_bounds(reeds_path, tech, access_case, subtech): - """Establish class cut offs based on capacity factors""" - class_def_name = 'reV_cf_ac' + return cap_exog - # Load the supply curve raw file produced by reV - df = pd.read_csv(os.path.join( - reeds_path,'inputs','supply_curve', - 'supplycurve_'+tech+'-'+access_case+'.csv')) - - # Aggregate min/max by class and attach access_case - if tech == 'wind-ofs': - df['subtech'] = 'fixed' - df.loc[df['class'].isin(WINDOFS_FLOATING_CLASSES),'subtech'] = 'floating' - df_sub = df[df['subtech']==subtech] - summary_df = df_sub.groupby('class')['cf'].agg(['min', 'max']).reset_index() - summary_df['subtech'] = subtech - summary_df['access_case'] = access_case - summary_df.columns = ['class', f'min_{class_def_name}', - f'max_{class_def_name}', 'subtech', 'access_case'] - else: - summary_df = df.groupby('class')['cf'].agg(['min', 'max']).reset_index() - summary_df['access_case'] = access_case - summary_df.columns = ['class', f'min_{class_def_name}', - f'max_{class_def_name}', 'access_case'] - - # Pin each class's min CF to the max CF of the previous class to avoid gaps - summary_df = summary_df.sort_values(by=['class',f'min_{class_def_name}']) - for c in summary_df['class'].unique().tolist(): - if c > min(summary_df['class'].unique().tolist()): - summary_df.loc[ - summary_df['class']==c,f'min_{class_def_name}' - ] = summary_df.loc[summary_df['class']==c-1][f'max_{class_def_name}'].iloc[0] - - # Round values to 4 decimal places - summary_df[f'min_{class_def_name}'] = summary_df[f'min_{class_def_name}'].round(4) - summary_df[f'max_{class_def_name}'] = summary_df[f'max_{class_def_name}'].round(4) - - return summary_df - -# Assign each wind, solar and geothermal unit in unit database to a class -def assign_class(cf, tech, df_class): - # Each unit is assigned to the class associated with the min and max performance - # its mean performance falls between - value = 'mean_temp' if tech in ['geohydro_allkm', 'egs_allkm'] else 'cf_ac' - row = df_class[(df_class[f'min_reV_{value}'] < cf) & (cf <= df_class[f'max_reV_{value}'])] - # Handle min cutoff point - if cf == df_class[f'min_reV_{value}'].min(): - row = df_class[cf == df_class[f'min_reV_{value}']] - - if len(row) == 1: - return row.iloc[0]['class'] - # If a offshore wind cf matches with both fixed and floating - # resources, assign a fixed resource - elif (len(row) > 1) & (tech == 'wind-ofs'): - row = row[row['subtech']=='fixed'] - return row.iloc[0]['class'] - else: - # If a unit's capacity factor/mean temp does not fall between any two max and min values - # specified in the classificalion file, it is unclassified and gives an error - raise ValueError('Unclassified ' + tech + ' technology with cf= ' + str(cf) + - ', check capacity factor/mean temperature values in unitdata.csv and classification files.') # Expand each row into multiple rows (startyear → retirement_year) def expand_exog_cap(row, start_year): @@ -248,31 +347,35 @@ def get_capacity_weighted_onlineyear(df, start_year, end_year): return pd.concat(output, ignore_index=True) -# Assign each calendar year to its appropriate modeledyear -# (For example: capacity that comes online in 2016 will -# show up in modeled year 2020) + def assign_modeledyear(x,years_list): + """ + Assign each calendar year to its appropriate modeledyear + (For example: capacity that comes online in 2016 will + show up in modeled year 2020) + """ for m in years_list: if x <= m: return m return None -# Check if there are any rsc techs in unitdata without resource classes -def missing_resource_class(gendb,rsc_class): + +def missing_resource_class(gdb_use_cap_exog, tech_class_cutoffs): + """Check if there are any rsc techs in unitdata without resource classes""" # Find the tech classes in unitdata that need to be matched to resource classes - matched_techs = [i for i in gendb['tech'].unique().tolist() if i in TECH['rsc_wsc']] + matched_techs = [i for i in gdb_use_cap_exog['tech'].unique() if i in TECH['rsc_wsc']] # Do not count csp-ns as it is matched to upv resources matched_techs = [i for i in matched_techs if i not in TECH['rsc_csp'] ] # Find tech classes in unitdata that are without assigned resource classes - missing_techs = list(set(matched_techs) - set(rsc_class)) + missing_techs = list(set(matched_techs) - set(tech_class_cutoffs)) if len(missing_techs) > 0: - raise ValueError(f'{missing_techs} are in unitdata but not matched with any resource classes. Exiting program.') + raise ValueError(f'{missing_techs} are in unitdata but not matched with any resource classes') else: print('All rsc/geothermal tech classes in unitdata are matched with available resource classes.') def process_ivt(years, inputs_case): - + """For vintages between model years, use the vintage of the next modeled year""" ivt_df= pd.read_csv(os.path.join(inputs_case,'ivt.csv')) ### modify set of technology name as lower case and convert all columns except the first to string ivt_df.iloc[:, 0] = ivt_df.iloc[:, 0].str.lower() @@ -294,137 +397,11 @@ def process_ivt(years, inputs_case): return ivt_df -# Only keep neccessary columns from unitdata to work with -# And rename column names for easier processing -COLNAMES = { - 'capexog_rsc': ( - ['tech','r','RetireYear','StartYear','sc_point_gid','summer_power_capacity_MW'], - ['tech','region','year','onlineyear','sc_point_gid','MW'] - ), - 'capnonrsc': ( - ['tech','coolingwatertech','r','ctt','wst','summer_power_capacity_MW'], - ['i','coolingwatertech','r','ctt','wst','value'] - ), - 'capnonrsc_energy': ( - ['tech','r','energy_capacity_MWh'], - ['i','r','value'] - ), - 'prescribed_nonRSC': ( - ['StartYear','tech','vin','r','coolingwatertech','ctt','wst','summer_power_capacity_MW'], - ['t','i','v','r','coolingwatertech','ctt','wst','value'] - ), - 'prescribed_nonRSC_energy': ( - ['StartYear','tech','vin','r','coolingwatertech','ctt','wst','energy_capacity_MWh'], - ['t','i','v','r','coolingwatertech','ctt','wst','value'] - ), - 'prescribed_RSC': ( - ['StartYear','tech','vin','r','summer_power_capacity_MW'], - ['t','i','v','r','value'] - ), - 'rsc': ( - ['tech','r','v','ctt','wst','summer_power_capacity_MW'], - ['i','r','v','ctt','wst','value'] - ), - 'rsc_wsc': ( - ['r','tech','summer_power_capacity_MW'], - ['r','i','value'] - ), - 'prsc_csp': ( - ['StartYear','r','tech','ctt','wst','summer_power_capacity_MW'], - ['t','r','i','ctt','wst','value'] - ), - 'prsc_geo': ( - ['StartYear','r','tech','summer_power_capacity_MW'], - ['t','r','i','value'] - ), - 'retirements': ( - ['tech','v','r','RetireYear','StartYear','coolingwatertech','ctt','wst','type','summer_power_capacity_MW'], - ['i','v','r','t','tt','coolingwatertech','ctt','wst','type','value'] - ), - 'retirements_energy': ( - ['tech','v','r','RetireYear','StartYear','type','energy_capacity_MWh'], - ['r','i','v','t','tt','type','value'] - ), - 'windret': ( - ['r','tech','RetireYear','summer_power_capacity_MW'], - ['r','i','t','value'] - ), - 'georet': ( - ['r','tech','RetireYear','summer_power_capacity_MW'], - ['r','i','t','value'] - ), - } -#%% =========================================================================== -### --- SUPPLEMENTAL DATA --- -### =========================================================================== - -######################### -### STATIC DICTIONARY ### -TECH = { - 'capnonrsc': [ - 'battery_li', 'biopower', 'coal-igcc', 'coal-new', - 'coaloldscr','coalolduns','gas-cc', 'gas-ct', - 'lfill-gas','nuclear', 'o-g-s', 'pumped-hydro' - ], - 'capnonrsc_energy': [ - 'battery_li' - ], - 'prescribed_nonRSC': [ - 'battery_li', 'biopower', 'coal-igcc', 'coal-new', - 'coaloldscr', 'coalolduns', 'gas-cc', 'gas-ct', - 'hydED', 'hydEND', 'hydUD', 'hydUND', 'hydND', 'hydNPND', - 'lfill-gas', 'nuclear', 'o-g-s', 'pumped-hydro' - ], - 'prescribed_nonRSC_energy': [ - 'battery_li', - ], - 'storage' : ['battery_li', 'pumped-hydro' - ], - 'rsc_pv_all': ['upv','pvb','pvb_pv','csp-ns'], - 'rsc_upv': ['upv','pvb'], - 'rsc_w': ['wind-ons','wind-ofs'], - 'rsc_csp': ['csp-ns'], - 'rsc_wsc': ['upv','pvb','csp-ns','csp-ws','wind-ons','wind-ofs', - 'geohydro_allkm','egs_allkm'], - 'prsc_csp': ['csp-ns','csp-ws'], - 'prsc_geo': ['geohydro_allkm','egs_allkm'], - 'retirements': [ - 'coalolduns', 'o-g-s', 'hydED', 'hydEND', 'gas-ct', 'lfill-gas', - 'coaloldscr', 'biopower', 'gas-cc', 'coal-new', - 'battery_li','nuclear', 'pumped-hydro', 'coal-igcc' - ], - 'retirements_energy': [ - 'battery_li' - ], - 'windret': ['wind-ons'], - 'georet': ['geohydro_allkm','egs_allkm'], - # This is not all technologies that do not having cooling, but technologies - # that are (or could be) in the plant database. - 'no_cooling': [ - 'upv', 'pvb', 'gas-ct', 'geohydro_allkm','egs_allkm', - 'battery_li', 'pumped-hydro', 'pumped-hydro-flex', - 'hydUD', 'hydUND', 'hydD', 'hydND', 'hydSD', 'hydSND', 'hydNPD', - 'hydNPND', 'hydED', 'hydEND', 'wind-ons', 'wind-ofs', - ], -} - - - #%% =========================================================================== ### --- MAIN FUNCTION --- ### =========================================================================== -def main(reeds_path, inputs_case): - - ######################### - ### SUPPLEMENTAL DATA ### - - quartershorten = {'spring':'spri','summer':'summ','fall':'fall','winter':'wint'} - - hotcold_months = {'NOV':'cold', 'DEC':'cold', 'JAN':'cold', 'FEB':'cold', - 'JUN':'hot', 'JUL':'hot', 'AUG':'hot' - } - +def main(inputs_case): #%% Inputs from switches sw = reeds.io.get_switches(inputs_case) GSw_WaterMain = int(sw.GSw_WaterMain) @@ -440,11 +417,6 @@ def main(reeds_path, inputs_case): regions = reeds.io.read_input(inputs_case, 'r').squeeze(1).values - #################### - ### DICTIONARIES ### - - - #%% print('Importing generator database:') gdb_use = pd.read_csv(os.path.join(inputs_case,'unitdata.csv'), @@ -714,59 +686,54 @@ def main(reeds_path, inputs_case): # Rename generic geothermal tech category to geohydro_allkm_1 geoexist['i'] = 'geohydro_allkm_1' - ###################################### + #%%################################### # -- RSC Exogenous Capacity -- # ###################################### - (cap_exog, rsc_class) = create_exog_rsc(reeds_path, inputs_case, gdb_use_cap_exog, TECH, COLNAMES, sw, startyear) - + tech_class_cutoffs = get_class_cf_min_alltechs(inputs_case) + cap_exog = create_exog_rsc(inputs_case, gdb_use_cap_exog, sw, startyear) #%%#################################### # -- RSC Prescribed Capacity -- # ####################################### print('Gathering RSC Prescribed Capacity...') + colnames = { + 'geohydro_allkm':'reV_mean_resource_temp', + 'egs_allkm':'reV_mean_resource_temp', + } + colname_default = 'reV_capacity_factor_ac' cap_pres = {} for tech in TECH['rsc_wsc']: + match tech: + case _ if tech in TECH['rsc_upv'] + TECH['prsc_csp']: + usetech = 'upv' + case 'wind-ofs': + usetech = ('wind-ofs', 'fixed') + case _: + usetech = tech + print(tech, usetech) + colname = colnames.get(tech,colname_default) + cap_pres[tech]= gdb_use.loc[(gdb_use['tech']==tech) & (gdb_use['StartYear'] >= startyear) & (gdb_use['StartYear'] <= endyear) ].copy() mask = ivt_df['Unnamed: 0'].str.contains(tech, case=False, na=False) if len(cap_pres[tech]) != 0: - # DUPV, PVB and UPV values are collected at the same time here: - if tech in TECH['rsc_upv']: - print(tech) - cap_pres[tech]["class"] = cap_pres[tech]["reV_capacity_factor_ac"].apply( - lambda x: assign_class(x, tech, rsc_class['upv'])) - cap_pres[tech]["tech"] = (cap_pres[tech]["tech"].astype(str) + "_" + - cap_pres[tech]["class"].astype(str)) - # Load in wind builds: - elif tech in TECH['rsc_w']: - print(tech) - cap_pres[tech]["class"] = cap_pres[tech]["reV_capacity_factor_ac"].apply( - lambda x: assign_class(x, tech, rsc_class[tech])) - cap_pres[tech]["tech"] = (cap_pres[tech]["tech"].astype(str) + "_" + - cap_pres[tech]["class"].astype(str)) - # Add prescribed csp builds: - # Note: Since csp is affected by GSw_WaterMain, it must be dealt with separate - # from the other RSC tech (dupv, upv, wind, etc) - elif tech in TECH['prsc_csp']: - print(tech) - cap_pres[tech]["class"] = cap_pres[tech]["reV_capacity_factor_ac"].apply( - lambda x: assign_class(x, tech, rsc_class['upv'])) - cap_pres[tech]["tech"] = (cap_pres[tech]["tech"].astype(str) + "_" + - cap_pres[tech]["class"].astype(str)) - if GSw_WaterMain == 1: - cap_pres[tech]["tech"] = np.where( cap_pres[tech]["tech"]=='csp-ws', - cap_pres[tech]["tech"]+'_'+cap_pres[tech]['ctt']+'_'+cap_pres[tech]['wst'], - 'csp-ws') - # Load in geo builds: - elif tech in TECH['prsc_geo']: - cap_pres[tech]["class"] = cap_pres[tech]["reV_mean_resource_temp"].apply( - lambda x: assign_class(x, tech, rsc_class[tech])) - cap_pres[tech]["tech"] = (cap_pres[tech]["tech"].astype(str) + "_" + - cap_pres[tech]["class"].astype(str)) + cap_pres[tech]["class"] = cap_pres[tech][colname].apply( + lambda x: assign_class(x, tech_class_cutoffs[usetech]) + ) + cap_pres[tech]["tech"] = ( + cap_pres[tech]["tech"].astype(str) + "_" + cap_pres[tech]["class"].astype(str) + ) + ## CSP is the only RSC tech with water cooling + if (tech in TECH['prsc_csp']) and (GSw_WaterMain == 1): + cap_pres[tech]["tech"] = np.where( + cap_pres[tech]["tech"]=='csp-ws', + cap_pres[tech]["tech"]+'_'+cap_pres[tech]['ctt']+'_'+cap_pres[tech]['wst'], + 'csp-ws' + ) # assign vintages based on start year of the unit ivt_df_mask = (ivt_df[mask] # filter rows .iloc[:, 1:] # drop first technology column @@ -1149,8 +1116,7 @@ def main(reeds_path, inputs_case): inputs_case = args.inputs_case # #%% Settings for testing - # reeds_path = reeds.io.reeds_path - # inputs_case = os.path.join(reeds_path,'runs','v20260804_inputsM0_Simple','inputs_case') + # inputs_case = os.path.join(reeds.io.reeds_path,'runs','v20261006_mcM0_github_MA_county_CC','inputs_case') #%% Set up logger log = reeds.log.makelog( @@ -1160,7 +1126,7 @@ def main(reeds_path, inputs_case): print('Starting writecapdat.py') #%% Run procedure - files_out = main(reeds_path, inputs_case) + files_out = main(inputs_case) #%% Write outputs print('Writing out capacity data') diff --git a/reeds/input_processing/writesupplycurves.py b/reeds/input_processing/writesupplycurves.py index 1a6aa4d49..1dfd3abd1 100644 --- a/reeds/input_processing/writesupplycurves.py +++ b/reeds/input_processing/writesupplycurves.py @@ -95,7 +95,7 @@ def agg_supplycurve( dfin = reeds.io.assemble_supplycurve( scfile=scpath, case=os.path.dirname(os.path.normpath(inputs_case)), - ).reset_index().drop(columns=['FIPS','cf'], errors='ignore') + ).reset_index().drop(columns=['FIPS'], errors='ignore') ## Convert dollar year and recalculate total cost transcost_cols = [c for c in dfin if 'cost' in c] dfin.loc[:, transcost_cols] *= deflate['interconnection'] @@ -218,7 +218,7 @@ def main( for s in wind_types: windin[s], wind[s] = agg_supplycurve( - scpath=os.path.join(inputs_case,f'supplycurve_wind-{s}.csv'), + scpath=Path(inputs_case, f'supplycurve_init_wind-{s}.csv'), inputs_case=inputs_case, numbins_tech=numbins[f'wind-{s}'], deflate=deflate, @@ -313,13 +313,14 @@ def main( os.path.join(inputs_case, f"exog_wind_{s}_rsc.csv") ) exog_onlineyear_list.append(exog_onlineyear) + windin[s].round(3).to_csv(Path(inputs_case, f'supplycurve_wind-{s}.csv')) # %%############### # -- PV -- # ################## upvin, upv = agg_supplycurve( - scpath=os.path.join(inputs_case, 'supplycurve_upv.csv'), + scpath=Path(inputs_case, 'supplycurve_init_upv.csv'), inputs_case=inputs_case, numbins_tech=numbins['upv'], deflate=deflate, @@ -355,6 +356,7 @@ def main( ) exog_upv_rsc.round(3).to_csv(os.path.join(inputs_case, "exog_upv_rsc.csv")) exog_onlineyear_list.append(exog_onlineyear) + upvin.round(3).to_csv(Path(inputs_case, 'supplycurve_upv.csv')) ### Normalize formatting upv = upv.reset_index() @@ -401,7 +403,7 @@ def main( if int(sw["GSw_CSP"]): _, csp = agg_supplycurve( - scpath=os.path.join(inputs_case, 'supplycurve_csp.csv'), + scpath=Path(inputs_case, 'supplycurve_init_csp.csv'), inputs_case=inputs_case, numbins_tech=numbins['csp'], deflate=deflate, @@ -480,7 +482,7 @@ def main( geoin[s], geo[s] = agg_supplycurve( scpath=os.path.join( inputs_case, - f'supplycurve_{s}.csv' + f'supplycurve_init_{s}.csv' ), numbins_tech=numbins[s], inputs_case=inputs_case, @@ -493,6 +495,8 @@ def main( .assign(rscbin="bin" + geo[s].reset_index()["bin"].astype(str)) .rename(columns={"region": "r"}) ) + if write: + geoin[s].round(3).to_csv(Path(inputs_case, f'supplycurve_{s}.csv')) geoall = ( pd.concat(geo, axis=0) diff --git a/reeds/inputs.py b/reeds/inputs.py index 256f95114..e1045895e 100644 --- a/reeds/inputs.py +++ b/reeds/inputs.py @@ -272,12 +272,12 @@ def parse_cases( "This may result in inconsistent representative days across MCS runs.\n\n" "To ensure consistency, we strongly recommend setting " "`GSw_HourlyClusterAlgorithm = user` in your switch configuration.\n" - "Do you want to proceed with the current setup?" ) - user_input = input("Type 'yes' to proceed, or 'no' to exit: ").strip().lower() - if user_input not in ['yes', 'y']: - print("\nPlease update the `GSw_HourlyClusterAlgorithm` switch and restart.") - quit() + if not skip_checks: + user_input = input("Do you want to proceed with the current setup? y/[n]") + if user_input.strip().lower() not in ['yes', 'y']: + print("\nPlease update the `GSw_HourlyClusterAlgorithm` switch and restart.") + quit() warned_about_cluster_alg = True print() @@ -441,7 +441,7 @@ def get_bin( bin_ser = bin_ser.rank(method='dense') df[bin_out_col] = bin_ser.values # we need the same index ordering for apply to work - df = df.reindex(index=orig_index) + df = df.reindex(index=orig_index).drop(columns='cum_cap') df[bin_out_col] = df[bin_out_col].astype(int) return df diff --git a/reeds/io.py b/reeds/io.py index 70ae5c611..aba3582a0 100644 --- a/reeds/io.py +++ b/reeds/io.py @@ -1099,7 +1099,48 @@ def get_temperatures(case, tz_in='UTC', tz_out='Etc/GMT+6', subset_years=True): return temperatures -def get_site_cf_hourly(tech, year, case=None, **kwargs): +def get_siting_switchval( + tech:None|Literal['upv','wind-ons','wind-ofs']=None, + switchname:None|Literal['GSw_SitingUPV','GSw_SitingWindOns','GSw_SitingWindOfs']=None, + case=None, + **kwargs +) -> str: + """ + Get the siting scenario to use for the specified tech. + If not using Monte Carlo sampling, simply reads from the case switches. + If using Monte Carlo sampling, returns the most permissive scenario used in the distribution + of siting scenarios for the specified technology (e.g., if the distribution includes + 'open' and 'limited', it returns 'open'). + """ + sw = reeds.io.get_switches(case, **kwargs) + if switchname is None: + switchname = { + 'upv': 'GSw_SitingUPV', + 'wind-ons': 'GSw_SitingWindOns', + 'wind-ofs': 'GSw_SitingWindOfs', + }[tech] + options_openest_first = ['open', 'reference', 'limited'] + if not int(sw['MCS_runs']): + return sw[switchname] + else: + mc_inputs = reeds.input_processing.mcs_sampler.load_mcs_dist(case) + keep = mc_inputs.assignments_list.map(lambda x: any([switchname in entry for entry in x])) + ## Siting should only be specified in one entry + assert keep.sum() <= 1 + if keep.sum() == 0: + return sw[switchname] + dist = mc_inputs.loc[keep].squeeze(0) + siting_scenarios = { + key: val for entry in dist.assignments_list for key, val in entry.items() + }[switchname] + ## Keep the openest scenario (which is first in the list) + for scen in options_openest_first: + if scen in siting_scenarios: + print(f'MCS: Using {scen}-access profiles for {switchname}') + return scen + + +def get_site_cf_hourly(tech:Literal['upv','wind-ons','wind-ofs'], year, case=None, **kwargs): """ Get hourly site-level capacity factor profiles for the given tech and year in UTC. Note that "distpv" is not a valid input to the "tech" parameter for @@ -1115,30 +1156,12 @@ def get_site_cf_hourly(tech, year, case=None, **kwargs): precedence. If {case} is None and a keyword argument is not provided for a switch, the default switch values specified in cases.csv are used. """ - sw = reeds.io.get_switches(case, **kwargs) - match tech: - case 'upv': - fname = f'cf_upv_{sw.GSw_SitingUPV}' - case 'wind-ons': - fname = f'cf_wind-ons_{sw.GSw_SitingWindOns}' - case 'wind-ofs': - fname = f'cf_wind-ofs_{sw.GSw_SitingWindOfs}' - case None: - raise ValueError( - "A technology must be provided if no case is " - "provided or if inputs_case/recf.h5 does not exist." - ) - case _: - raise NotImplementedError( - f"The provided tech '{tech}' does not have CF profiles." - ) + access = get_siting_switchval(tech=tech, case=case, **kwargs) + fname = f'cf_{tech}_{access}.h5' + h5path = Path(reeds_path, 'inputs', 'profiles_cf', fname) + if not h5path.is_file(): + raise FileNotFoundError(h5path) - h5path = os.path.join( - reeds_path, - 'inputs', - 'profiles_cf', - f'{fname}.h5' - ) with h5py.File(h5path, 'r') as f: time_index = pd.to_datetime( pd.Series(f[f'time_index_{year}'][:]) diff --git a/runreeds.py b/runreeds.py index 1572ce3ef..520294e06 100644 --- a/runreeds.py +++ b/runreeds.py @@ -1296,8 +1296,8 @@ def write_batch_script( big_comment('Input processing', OPATH) for s in [ 'copy_files', - 'process_unitdata', 'mcs_sampler', + 'process_unitdata', 'climateprep', 'hydcf', 'h2_storage',