diff --git a/cases.csv b/cases.csv index 525a78bc..11b984fa 100644 --- a/cases.csv +++ b/cases.csv @@ -225,6 +225,8 @@ GSw_MGA_CostDelta,MGA: Fraction by which to allow objective function to increase GSw_MGA_Direction,MGA: Directionality of second optimization,(min|max),min, GSw_MGA_Objective,MGA: Objective for MGA (uses GSw_MGA_SubObjective to specify technology subset if set to capacity),(capacity|generation|transmission|rasharing|co2|employment),capacity, GSw_MGA_SubObjective,MGA: Technology subset to minimize or maximize the capacity or generation of (only used for GSw_MGA_Objective=(capacity or generation)),(battery|ccs|coal|dac|fossil|gas|gentech|geo|h2_combustion|hydro|nuclear|ofswind|onswind|pv|re|storage|upv|vre|wind),gentech, +GSw_MGA_RV_runs,"MGA: Number of samples drawn for the random vector weights that are applied to variables specified in GSw_MGA_Objective and GSw_MGA_SubObjective. Only active when GSw_MGA_CostDelta != 0; set to an integer N>0 to run N ReEDS simulations with sampling (0 disables sampling).",int,0, +GSw_MGA_RV_region,"MGA: regionality to use with random vector method (active when GSw_MGA_RV_runs>=1)",r; nercr; transreg; transgrp; cendiv; st; interconnect; country; usda_region,r, GSw_MinCF,Turn on/off regional min CF constraint (applied at i/r level),0; 1,1, GSw_Mingen,Turn on/off min-gen constraints by r/h/szn,0; 1,0, GSw_MingenFixed,Turn on/off fixed min-gen constraints,0; 1,1, diff --git a/cases_test.csv b/cases_test.csv index f3a122ad..6afd94a0 100644 --- a/cases_test.csv +++ b/cases_test.csv @@ -1,66 +1,67 @@ -,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_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, -endyear,2032,,2050,2050,2050,2029,2060,2026,,,,,,,,,,,2035,2030,2030,2060,2035,2050,2050,2050,,,,, -yearset,,,,,,,2010..2060..10,,,,,,,,,,,,,2010..2050..5,2010..2050..5,2010..2060..10,,2010..2050..5,2010_2025_2050,2010..2050..5,,,,, -GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,z134,z54,z48,,,,,, -GSw_GasCurve,2,,1,1,,,,,,,,,,,,,,,,,,,,1,1,,,,,, -GSw_Geothermal,,,,2,,,,,,,,,,,,,,,,,,,0,,0,,,,,, -GSw_GrowthPenalties,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,, -GSw_Upstream,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,, -GSw_TransHurdleRate,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,, -distpvscen,,,,,stscen2023_mid_case_95_by_2035,,,,,,,,,,,,,,,,,,,,,,,,,, -GSw_AnnualCap,,,,,2,,1,,,,,,,,,,,,,,,1,,,,,,,,, -GSw_AnnualCapScen,,,,,start2024_90pct2035_100pct2045,,start2027_95pct2035,,,,,,,,,,,,,,,start2027_95pct2035,,,,,,,,, -GSw_LoadProfiles,,,,,EER2025_100by2050,EER2025_IRAlow,EER2025_IRAlow,EER2025_IRAlow,,,,,historic,,,,,,,,,EER2025_100by2050,,,,,historic,,,, -GSw_NG_CRF_penalty,,,,,ramp_2045,,ramp_2023_2035,,,,,,,,,,,,,,,ramp_2023_2035,,,,,,,,, -GSw_PRM_NetImportLimit,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,, -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_2016_2017_2018_2019_2020_2021_2022_2023,,,,,,,,2012_2013,,,,,2018,,,, -GSw_HourlyClusterMapMethod,,,,,,,,,,,,,,bestfirst,,,,,,,,,,,,,,,,, -GSw_WaterCapacity,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,,, -GSw_WaterMain,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,,, -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_lhs,,,,,,,,,,,,,,,,,,,,,1,,,,,,,,,, -GSw_PRM_StressIterateMax,,,,,,,,0,,,,,,,,0,0,,,1,1,,,,,,,,,, -GSw_ReducedResource,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,, -GSw_SitingUPV,,,,,,limited,limited,limited,,,,,,,,,,,,,,limited,,,,,,,,, -GSw_SitingWindOfs,,,,,,limited,limited,limited,,,,,,,,,,,,,,open,,,,,,,,, -GSw_SitingWindOns,,,,,,limited,limited,limited,,,,,,,,,,,,,,limited,,,,,,,,, -GSw_TransScen,,,,,,,NTP_MT,,,,,,,,,,,,,,,NTP_MT,,,,,,,,, -GSw_CO2_Detail,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,, -GSw_DAC,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,, -GSw_NoFossilOffsetCDR,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,, -GSw_Biopower,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,, -GSw_HourlyChunkLengthRep,,,,,,,,,,,,,,,,,,,,,,,6,4,4,,,,,, -GSw_HourlyChunkLengthStress,,,,,,,,,,,,,,,,,,,,,,,6,4,4,,,,,, -GSw_LfillGas,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,, -GSw_Nuclear,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,, -GSw_OpRes,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,, -GSw_StartCost,,,,,,,,,,,,,,,,,,,,,,,0,0,0,,,,,, -GSw_H2,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,, -GSw_H2_PTC,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,, -GSw_H2Combustion,,,,,,,,,,,,,,,,,,,,,,,,,0,,,,,, -GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth,,,,, -GSw_DRShed,,,,,,,,,,,,,,,,,,,,,,,,,,,1,,,, -GSw_MGA_CostDelta,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.01,,, -GSw_LoadSiteCF,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.95,, -GSw_OffshoreZones,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, -GSw_OffshoreBackbone,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, -GSw_OffshoreBackflow,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, -pras_agg_ogs_lfillgas,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1 -pras_existing_unit_size,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0 -pras_scheduled_outage,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0 -pras_unitsize_source,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,r2x -pras_vre_combine,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1 -pras_samples,,,,,,10,10,10,,,,,,,,,,,,,,,,10,10,,,,,, +,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, +endyear,2032,,2050,2050,2050,2029,2060,2026,,,,,,,,,,,2035,2030,2030,2060,2035,2050,2050,2050,,,,,, +yearset,,,,,,,2010..2060..10,,,,,,,,,,,,,2010..2050..5,2010..2050..5,2010..2060..10,,2010..2050..5,2010_2025_2050,2010..2050..5,,,,,, +GSw_ZoneSet,,,,,,,z54,z3109,,,,,,,z3109,z3109,z3109,PJMcounty,,,,z54,z134,z54,z48,,,,,,, +GSw_GasCurve,2,,1,1,,,,,,,,,,,,,,,,,,,,1,1,,,,,,, +GSw_Geothermal,,,,2,,,,,,,,,,,,,,,,,,,0,,0,,,,,,, +GSw_GrowthPenalties,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,,, +GSw_Upstream,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,,, +GSw_TransHurdleRate,,,,1,,,,,,,,,,,,,,,,,,,,,,,,,,,, +distpvscen,,,,,stscen2023_mid_case_95_by_2035,,,,,,,,,,,,,,,,,,,,,,,,,,, +GSw_AnnualCap,,,,,2,,1,,,,,,,,,,,,,,,1,,,,,,,,,, +GSw_AnnualCapScen,,,,,start2024_90pct2035_100pct2045,,start2027_95pct2035,,,,,,,,,,,,,,,start2027_95pct2035,,,,,,,,,, +GSw_LoadProfiles,,,,,EER2025_100by2050,EER2025_IRAlow,EER2025_IRAlow,EER2025_IRAlow,,,,,historic,,,,,,,,,EER2025_100by2050,,,,,historic,,,,, +GSw_NG_CRF_penalty,,,,,ramp_2045,,ramp_2023_2035,,,,,,,,,,,,,,,ramp_2023_2035,,,,,,,,,, +GSw_PRM_NetImportLimit,,,,,0,,,,,,,,,,,,,,,,,,,,,,,,,,, +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_2016_2017_2018_2019_2020_2021_2022_2023,,,,,,,,2012_2013,,,,,2018,,,,, +GSw_HourlyClusterMapMethod,,,,,,,,,,,,,,bestfirst,,,,,,,,,,,,,,,,,, +GSw_WaterCapacity,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,,,, +GSw_WaterMain,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,,,, +GSw_WaterUse,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,,,, +resource_adequacy_years,,,,,,,2011_2012_2013_2021_2022_2023,,,,,,,,,,,,,,,2011_2012_2013_2021_2022_2023,,,,,,,,,, +GSw_HourlyClusterAlgorithm,,,,,,,,,,,,,,,,,,,,user,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_lhs,,,,,,,,,,,,,,,,,,,,,1,,,,,,,,,,, +GSw_PRM_StressIterateMax,,,,,,,,0,,,,,,,,0,0,,,1,1,,,,,,,,,,, +GSw_ReducedResource,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,,, +GSw_SitingUPV,,,,,,limited,limited,limited,,,,,,,,,,,,,,limited,,,,,,,,,, +GSw_SitingWindOfs,,,,,,limited,limited,limited,,,,,,,,,,,,,,open,,,,,,,,,, +GSw_SitingWindOns,,,,,,limited,limited,limited,,,,,,,,,,,,,,limited,,,,,,,,,, +GSw_TransScen,,,,,,,NTP_MT,,,,,,,,,,,,,,,NTP_MT,,,,,,,,,, +GSw_CO2_Detail,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,,, +GSw_DAC,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,,, +GSw_NoFossilOffsetCDR,,,,,,,1,,,,,,,,,,,,,,,1,,,,,,,,,, +GSw_Biopower,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,,, +GSw_HourlyChunkLengthRep,,,,,,,,,,,,,,,,,,,,,,,6,4,4,,,,,,, +GSw_HourlyChunkLengthStress,,,,,,,,,,,,,,,,,,,,,,,6,4,4,,,,,,, +GSw_LfillGas,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,,, +GSw_Nuclear,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,,, +GSw_OpRes,,,,,,,,,,,,,,,,,,,,,,,0,,0,,,,,,, +GSw_StartCost,,,,,,,,,,,,,,,,,,,,,,,0,0,0,,,,,,, +GSw_H2,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,,, +GSw_H2_PTC,,,,,,,,,,,,,,,,,,,,,,,,0,0,,,,,,, +GSw_H2Combustion,,,,,,,,,,,,,,,,,,,,,,,,,0,,,,,,, +GSw_PRM_StressThresholdMetrics,,,,,,,,,,,,,,,,,,,,,,,,,,NEUE/LOLH/LOLE/LOLD/duration/depth,,,,,, +GSw_DRShed,,,,,,,,,,,,,,,,,,,,,,,,,,,1,,,,, +GSw_MGA_CostDelta,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.01,0.01,,, +GSw_MGA_RV_runs,,,,,,,,,,,,,,,,,,,,,,,,,,,,,2,,, +GSw_LoadSiteCF,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0.95,, +GSw_OffshoreZones,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, +GSw_OffshoreBackbone,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, +GSw_OffshoreBackflow,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1, +pras_agg_ogs_lfillgas,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1 +pras_existing_unit_size,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0 +pras_scheduled_outage,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,0 +pras_unitsize_source,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,r2x +pras_vre_combine,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1 +pras_samples,,,,,,10,10,10,,,,,,,,,,,,,,,,10,10,,,,,,, \ No newline at end of file diff --git a/docs/source/user_guide.md b/docs/source/user_guide.md index 9de22098..e856602b 100644 --- a/docs/source/user_guide.md +++ b/docs/source/user_guide.md @@ -631,8 +631,17 @@ Options are the column names in the `inputs/tech-subset-table.csv` file. Users familiar with GAMS can add alternative objective functions to the `d_mga.gms` file and associated options to the `GSw_MGA_Objective` switch in `cases.csv`. +By default the MGA min/max is applied to the sum of the variable across all regions being modeled. +The MGA approach also supports an option to randomly sample of a vector of weights to apply to the regional values of the variable being optimized. +This method can be useful to characterizing the uncertainty in the regional distribution of the results. +Weights are sampled as discrete values from a support of {-1,1} to allow for simultaneous minimization and maximization. +The MGA random vector capability is controlled by the following switches: +- `GSw_MGA_RV_runs` (default `0`): Number of random samples of weight vectors to draw; corresponds to the number of runs. +- `GSw_MGA_RV_region` (default `r`): Regionality level (specified by hierarhcy file) over which to sample the random weights. +Note that this capability is currently only supported when `GSw_MGA_Objective = (capacity or generation)`. +Weights for each run are stored in the `mga_weights` parameter. ## Uncertainty Plots diff --git a/reeds/core/setup/b_inputs.gms b/reeds/core/setup/b_inputs.gms index a008fbff..195409dc 100644 --- a/reeds/core/setup/b_inputs.gms +++ b/reeds/core/setup/b_inputs.gms @@ -93,6 +93,8 @@ $gdxin set land(r) "land-based (not offshore) zones" ; land(r)$[not offshore(r)] = yes ; + + sets *The following two sets: *ban - will remove the technology from being considered, anywhere @@ -6025,6 +6027,35 @@ employment_factor_plant(i,"construction") = employment_factor_plant(i,"construction") * upgrade_ratio(i) ; $endif.upgrade_ef +*==================================== +* --- MGA Random Vector Weights --- +*==================================== + +$ifthene.mgaobj ((sameas(%GSw_MGA_Objective%,capacity))or(sameas(%GSw_MGA_Objective%,generation))) + +parameter mga_weights(r,i_subtech) "--unitless-- weight to assign to given MGA subobjective by region" ; + +$ifthene.mga_rv (%GSw_MGA_RV_runs%>=1) +parameter mga_weights_in(r,i_subtech) +/ +$offlisting +$ondelim +$include inputs_case%ds%mga_weights.csv +$offdelim +$onlisting +/ ; +mga_weights(r,i_subtech) = mga_weights_in(r,i_subtech) ; +$else.mga_rv +mga_weights(r,i_subtech) = 1 ; +$endif.mga_rv + +$endif.mgaobj + + +*TODO: add different handling for other subobjectives with different dimensions +* should we also include error checking? + + *================================================================================================ *== h- and szn-dependent sets and parameters (declared here, populated in 2_temporal_params) === *================================================================================================ diff --git a/reeds/core/setup/d_mga.gms b/reeds/core/setup/d_mga.gms index c4802857..e9874ab4 100644 --- a/reeds/core/setup/d_mga.gms +++ b/reeds/core/setup/d_mga.gms @@ -21,7 +21,25 @@ eq_MGA_Objective$Sw_MGA.. $[tmodel(t) $valcap(i,v,r,t) $%GSw_MGA_SubObjective%(i)], - CAP(i,v,r,t) + CAP(i,v,r,t) + * sum{i_subtech$i_subsets(i,i_subtech), mga_weights(r,i_subtech)} + } +; + +* --------------------------------------------------------------------------- + +$elseif.mgaobj %GSw_MGA_Objective% == 'generation' +Equation eq_MGA_Objective "--MW-- Defines generation for MGA" ; +Variable MGA_OBJ "--MWh-- Generation of technology to be minimized/maximied" ; +eq_MGA_Objective$Sw_MGA.. + MGA_OBJ + =e= + sum{(i,v,r,h,t) + $[tmodel(t) + $valgen(i,v,r,t) + $%GSw_MGA_SubObjective%(i)], + GEN(i,v,r,h,t) * hours(h) + * sum{i_subtech$i_subsets(i,i_subtech), mga_weights(r,i_subtech)} } ; diff --git a/reeds/input_processing/mcs_sampler.py b/reeds/input_processing/mcs_sampler.py index 626eb624..0a1dcfb9 100644 --- a/reeds/input_processing/mcs_sampler.py +++ b/reeds/input_processing/mcs_sampler.py @@ -14,6 +14,7 @@ import pandas as pd import scipy.stats import sys +import re import yaml from typing import Tuple, List from collections import defaultdict @@ -1807,7 +1808,7 @@ def write_samples( #%% =========================================================================== ### --- MAIN PROCEDURE --- ### =========================================================================== -def main( +def main_mcs( reeds_path: str, inputs_case: str, n_samples: int = 1, @@ -1872,6 +1873,88 @@ def main( # Write Samples write_samples(sample_group, samples_dict, aux_files) +def main_mga_rv( + reeds_path: str, + inputs_case: str, + sw: dict, + n_samples: int = 1, + lhs_sampling: int = 1, + seed: int = 0, + discrete: bool = True, +): + + # get dimensions based on number of regions and subojectives + + ## regions + # get list of valid regions (val_r generated in copy_files.py) + val_r = list(reeds.io.read_input(inputs_case, 'r')['*']) + + # option to draw samples based on aggregated regions (samples will be mapped back to r regions) + hierarchy = reeds.io.get_hierarchy(GSw_ZoneSet=sw['GSw_ZoneSet']).reset_index() + hierarchy_val_r = hierarchy.loc[hierarchy.r.isin(val_r)] + sampling_regions = hierarchy_val_r[sw['GSw_MGA_RV_region']].unique() + + ## objective (assumes sw.GSw_MGA_Objective in ['capacity', 'generation'] based on + ## check in runreeds.check_compatibility() + ## if the subojective is an aggregated category, break it up into smaller groups + ## otherwise just use the subobjective as the group + mapped_categories ={ + 'gentech': ['coal', 'gas', 'nuclear', 'h2_combustion', 'geo', 'hydro', 'ofswind', 'onswind', 'storage', 'upv'], + 'fossil' : ['coal', 'gas'], + 're': ['geo', 'hydro', 'ofswind', 'onswind', 'pv', 'vre', 'wind'], + 'vre': ['ofswind', 'onswind', 'upv'], + } + subsets = mapped_categories.get(sw.GSw_MGA_SubObjective, [sw.GSw_MGA_SubObjective]) + + ## sample weights for each subojective group and sampling region + dimensions = len(subsets) * len(sampling_regions) + + # setup output + runs_folder_name = os.path.basename(os.path.dirname(inputs_case.rstrip(os.path.sep))) + mga_run_number = int((runs_folder_name.split('_')[-1]).replace('R', '')) + region_labels = np.repeat(sampling_regions, len(subsets)) + subset_labels = np.tile(subsets, len(sampling_regions)) + + # sample using LHS or random approach + if lhs_sampling: + # lhs requires drawing all samples simultaneously, so rather than using + # a run-specific seed we draw for all runs at once using the global seed value + lhs_sampler = scipy.stats.qmc.LatinHypercube(d=dimensions, seed=seed) + # lhs_samples are arranged n x d (n = samples, d = dimensions) + lhs_samples_cdf = lhs_sampler.random(n=n_samples) + if discrete: + # bin CDF samples in discrete choices (-1 or 1 with equal probability) + lhs_samples = np.where(lhs_samples_cdf < 0.5, -1, 1) + else: + # translate CDF samples into weights using uniform distribution (-1 to 1 to allow for simultaneous min/max) + lhs_samples = scipy.stats.uniform.ppf(lhs_samples_cdf, loc=-1, scale=2) + + # record the lhs sampling matrix in each run folder + lhs_samples_out = pd.DataFrame(lhs_samples.round(6)).T + lhs_samples_out.columns = [f"R{i:0>4}" for i in range(1, n_samples + 1)] + lhs_samples_out.index = [f"{region_labels[i]}_{subset_labels[i]}" for i in range(len(region_labels))] + lhs_samples_out.index.name = 'dimension' + lhs_samples_out.to_csv(os.path.join(inputs_case, "mga_rv_latin_hypercube_samples.csv")) + + # get the weights for this specific run (-1 to adjust for zero indexing) + mga_weights_raw = lhs_samples[mga_run_number - 1] + else: + # set random seed using the global seed + MGA run number to allow reproducibility + np.random.seed(seed + mga_run_number) + # get the weights for this specific run (-1 to 1 to allow for simultaneous min/max) + if discrete: + mga_weights_raw = np.random.choice([-1,1], dimensions) + else: + mga_weights_raw = np.random.uniform(-1, 1, dimensions) + + # save vector of weights for this run, mapped back to r regions and rounded to 6 decimal places + mga_weights = pd.DataFrame({sw['GSw_MGA_RV_region']: region_labels, 'i_subtech': subset_labels, 'weight': mga_weights_raw.round(6)}) + if sw['GSw_MGA_RV_region'] != 'r': + mga_weights = pd.merge(mga_weights, hierarchy_val_r[[sw['GSw_MGA_RV_region'], 'r']], on=sw['GSw_MGA_RV_region']) + 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) + if __name__ == '__main__' and not hasattr(sys, 'ps1'): parser = argparse.ArgumentParser(description='Copy files needed for this run') @@ -1902,17 +1985,24 @@ def main( sw = reeds.io.get_switches(inputs_case) MCS_runs = int(sw.get('MCS_runs', 0)) MCS_lhs = int(sw.get('MCS_lhs', 0)) + GSw_MGA_RV_runs = int(sw.get('GSw_MGA_RV_runs', 0)) # get global seed from scalars (used to set the seed for a batch of runs) scalars = reeds.io.get_scalars() seed = int(scalars['MCS_seed']) if MCS_runs >= 1: - print('Starting mcs_sampler.py') - main(reeds_path, inputs_case, n_samples=MCS_runs, lhs_sampling=MCS_lhs, seed=seed) + print('Starting Monte Carlo sampling with mcs_sampler.py') + main_mcs(reeds_path, inputs_case, n_samples=MCS_runs, lhs_sampling=MCS_lhs, seed=seed) else: print('MCS_runs switch is set to 0 or not found. No Monte Carlo sampling will be performed') + if GSw_MGA_RV_runs >= 1: + print('Starting random vector sampling for MGA with mcs_sampler.py') + main_mga_rv(reeds_path, inputs_case, sw, n_samples=GSw_MGA_RV_runs, lhs_sampling=MCS_lhs, seed=seed) + else: + print('GSw_MGA_RV_runs switch is set to 0 or not found. No MGA random vector sampling will be performed') + # Final log/timing update. reeds.log.toc( tic=tic, @@ -1920,3 +2010,5 @@ def main( process='input_processing/mcs_sampler.py', path=os.path.join(os.path.dirname(inputs_case)) ) + +# %% diff --git a/reeds/input_processing/runfiles.csv b/reeds/input_processing/runfiles.csv index a4e55de5..ceda4e6f 100644 --- a/reeds/input_processing/runfiles.csv +++ b/reeds/input_processing/runfiles.csv @@ -132,7 +132,7 @@ i_p.csv,inputs/sets/i_p.csv,1,ignore,ignore,,,,,0,,,set,i_p,mapping from technol i_water_nocooling.csv,inputs/sets/i_water_nocooling.csv,1,ignore,ignore,,,,,0,,,set,i_water_nocooling,technologies that use water but are not differentiated by cooling tech and water source, inflation.csv,inputs/financials/inflation_{inflation_suffix}.csv,1,ignore,ignore,,,,0,0,,,,,, interconnection_queues.csv,inputs/capacity_exogenous/interconnection_queues.csv,1,ignore,ignore,r,"tg,r",,1,0,,,,,, -jtype.csv,inputs/sets/jtype.csv,1,ignore,ignore,,,,,,,,set,jtype,job types used in model (construction and om), +jtype.csv,inputs/sets/jtype.csv,1,ignore,ignore,,,,,0,,,set,jtype,job types used in model (construction and om), lcclike.csv,inputs/sets/lcclike.csv,1,ignore,ignore,,,,,0,,,set,lcclike,transmission capacity types where lines are bundled with AC/DC converters, load_multiplier.csv,inputs/load/demand_{demandscen}.csv,1,ignore,ignore,,,,,0,,,,,, loadsite_annual.csv,inputs/load/loadsite_{GSw_LoadSiteTrajectory}.csv,float(sw.GSw_LoadSiteCF) > 0,ignore,ignore,*loadsitereg,t,,,0,,,,,, @@ -269,4 +269,4 @@ wst_climate.csv,inputs/sets/wst_climate.csv,1,ignore,ignore,,,,,0,,,set,wst_clim yearafter.csv,inputs/sets/yearafter.csv,1,ignore,ignore,,,,,,,,set,yearafter,set to loop over for the final year calculation, Project.toml,Project.toml,1,ignore,ignore,,,,,,,,,,, gamslice.txt,gamslice.txt,0,ignore,ignore,,,,,,,,,,, -runreeds.py,runreeds.py,1,ignore,ignore,,,,,,,,,,, +runreeds.py,runreeds.py,1,ignore,ignore,,,,,,,,,,, \ No newline at end of file diff --git a/reeds/inputs.py b/reeds/inputs.py index c4de0abd..eb62b1fa 100644 --- a/reeds/inputs.py +++ b/reeds/inputs.py @@ -246,14 +246,14 @@ def parse_cases( print("Please change the delimeter in the GSw_Region switch from ',' to '.'") quit() - # If doing a Monte Carlo run, modify dfcases by adding new columns - # for each scenario run. Also validate the distribution file. + # If doing a Monte Carlo or MGA Random Vector run, modify dfcases by adding new columns + # for each scenario run. For Monte Carlo we also validate the distribution file. warned_about_cluster_alg = False - if 'MCS_runs' in dfcases.index: + if 'MCS_runs' in dfcases.index or 'GSw_MGA_RV_runs' in dfcases.index: for c in dfcases.columns: if ( c not in ['Description','Default Value','Choices'] - and (int(dfcases.loc['MCS_runs',c]) > 0) + and ((int(dfcases.loc['MCS_runs',c]) > 0) or (int(dfcases.loc['GSw_MGA_RV_runs',c]) > 0)) and (not int(dfcases.loc['ignore',c])) ): # Warn user if the hourly clustering algorithm is not fixed for Monte Carlo runs @@ -263,7 +263,7 @@ def parse_cases( ): print(f"\n[Warning] Case Column: '{c}'") print( - "You are attempting to run a Monte Carlo simulation with " + "You are attempting to run a Monte Carlo or MGA Random Vector simulation with " "`GSw_HourlyClusterAlgorithm` set to a value other than 'user'.\n" "This may result in inconsistent representative days across MCS runs.\n\n" "To ensure consistency, we strongly recommend setting " @@ -283,23 +283,29 @@ def parse_cases( reeds.io.reeds_path, 'inputs', 'userinput', 'mcs_distributions_{}.yaml'.format(sw.MCS_dist) ) - mcs_sampler.general_mcs_dist_validation(reeds.io.reeds_path, mcs_dist_path, sw) - - # c (column) is a case with monte carlo runs. - # replicate this column N (NumMonteCarloRuns) times - NumMonteCarloRuns = int(dfcases.loc['MCS_runs',c]) + if int(dfcases[c].MCS_runs) > 0: + mcs_sampler.general_mcs_dist_validation(reeds.io.reeds_path, mcs_dist_path, sw) + numruns = int(dfcases[c].MCS_runs) + run_type = 'MC' + else: + numruns = int(dfcases[c].GSw_MGA_RV_runs) + run_type = 'R' + + # c (column) is a case with monte carlo or MGA random vector runs. + # replicate this column N times NewColumnNames = [ - f"{c}_MC{i:0>4}" - for i in range(1, NumMonteCarloRuns + 1) + f"{c}_{run_type}{i:0>4}" + for i in range(1, numruns + 1) ] - # Each new column is a copy of the original column with name c_{MC1,MC2,...} - dfcases_MC = pd.DataFrame( - data=np.array([dfcases[c].values]*NumMonteCarloRuns).T, + # Each new column is a copy of the original column with name + # c_{MCS1,MCS2,...} or c_{R1,R2,...} + dfcases_all = pd.DataFrame( + data=np.array([dfcases[c].values]*numruns).T, index=dfcases.index, columns=NewColumnNames, ) - dfcases = pd.concat([dfcases, dfcases_MC], axis=1) + dfcases = pd.concat([dfcases, dfcases_all], axis=1) # drop the original column dfcases.drop(c, axis=1, inplace=True) @@ -352,6 +358,7 @@ def solvestring_sequential( 'GSw_HourlyWrapLevel', 'GSw_MGA_CostDelta', 'GSw_MGA_Direction', + 'GSw_MGA_RV_runs', 'GSw_PVB_Dur', 'GSw_SkipRAyear', 'GSw_StateCO2ImportLevel', diff --git a/runreeds.py b/runreeds.py index 250eea96..78c0fd6b 100644 --- a/runreeds.py +++ b/runreeds.py @@ -487,6 +487,11 @@ def check_compatibility(sw): ) raise ValueError(err) + if (sw['GSw_MGA_Objective'] not in ['capacity', 'generation']) and int(sw['GSw_MGA_RV_runs']) > 0: + raise NotImplementedError( + f"GSw_MGA_Objective='{sw['GSw_MGA_Objective']}' is not yet supported for MGA random vector sampling." + ) + ### Dependent model availability if ( ((int(sw['pras']) == 2) or int(sw['GSw_PRM_StressIterateMax']))