Implement daily temperature-based gas price adjustments

[SoLaraS2]

Summary

This PR calculates daily gas price multipliers based on a linear regression fitting heating and cooling degree days on daily deviations from annual regional gas prices (the parameters of which are in inputs/fuelprices/gasreg_degree_day_price_mult_regression_params.csv). The user can choose between these daily regional multipliers, the existing national wintertime markup, and no adjustment.

Technical details

• Added gasreg hierarchy level, which represents the regions at which the degree day / price deviation regressions were fitted. These are mostly just census divisions, with the exceptions being "California" and "Northwest", which are components of the "Pacific" cendiv, and "Southwest" and "Mountain", which are components of the "Mountain" cendiv.
• Daily gasreg-level gas price multipliers are calculated in reeds/input_processing/fuelcostprep.py by calculating historical daily degree days for the weather years of the given run and then rescaling them to match annual degree days projected across the run's solve years. To create multipliers at the r level, the gasreg-level multipliers are copied to their constituent zones. To create multipliers at the cendiv level, the gasreg-level multipliers are aggregated via population-weighted average.
• Added get_daily_gas_price_multipliers function to hourly_writetimeseries.py to get a GAMS compatible csv with representative hourly NG price multipliers. In this function, the multipliers for each region and year are renormalized so that the average across the representative periods is 1, ensuring the year-round average gas price doesn't change.

Additional changes

• Added a new switch GSw_GasPriceAdjMethod which controls whether the daily price multipliers are applied, the national wintertime markup is applied, or no adjustment is applied.
• In main, the wintertime markup of 1.054 (szn_adj_gas_winter in scalars.csv) mentioned in the model documentation is only being applied when GSw_GasCurve = 1. A separate wintertime multiplier of 1.04 (gasprice_ref_frac_adder in scalars.csv) is being applied when GSw_GasCurve = 3. No multiplier is applied when GSw_GasCurve = 0 or 2. In this PR, the gasprice_ref_frac_adder multiplier is deleted and the szn_adj_gas multiplier is applied in every GSw_GasCurve scenario (assuming GSw_GasPriceAdjMethod = 1).
• In main, the wintertime markup only increases the wintertime price but doesn't make corresponding adjustments to maintain the pre-markup year-round average price. As a result, the year-round natural gas price is being erroneously inflated. In this PR, the seasonal multipliers are renormalized so that wintertime prices are higher than nonwinter prices while the year-round average price remains unchanged.
• In main, there are four legacy zones (p32, p35, p47, and p59) that are assigned to census divisions that don't match their state's census division (according to https://www.eia.gov/consumption/commercial/maps.php#census:~:text=CBECS%20climate%20zones-,U.S.%20census%20regions,-and%20divisions%3A). This was initially done to remedy a supply curve related issue that is now fixed, so this PR reassigns those zones to their states' census divisions for consistency with those zones' assigned gasregs.

Switches added/removed/changed

• Added GSw_GasPriceAdjMethod (0 = no adjustment, 1 = national wintertime markup, 2 = daily adjustments based on regional temperatures)

Validation, testing, and comparison report(s)

Results comparing main and this branch using different GSw_GasPriceAdjMethod and GSw_GasCurve values are below.

GSw_GasCurve = 1:

• GSw_GasPriceAdjMethod = 1:
  • results-0618_Main_USA_defaults,0618_NatGas_USA_defaults.pptx
  • Slightly more natural gas is built due to the fact that the wintertime markup is now renormalized so that the year-round average gas price remains unchanged (previously we only increased the wintertime gas price without a commensurate decrease in the non-wintertime price, which increased the year-round average price).
• GSw_GasPriceAdjMethod = 2:
  • results-0618_Main_USA_defaults,0618_NatGas_USA_DailyAdj.pptx
  • More natural gas is built than when GSw_GasPriceAdjMethod = 1, but still not a large amount
• GSw_GasPriceAdjMethod = 0:
  • results-0618_Main_USA_defaults,0618_NatGas_USA_NoAdj.pptx
  • Similar to the other cases, no adjustment results in the year-round average price being unchanged when it was previously being inflated, so more natural gas is built.

GSw_GasCurve = 2

• GSw_GasPriceAdjMethod = 1:
  • results-0618_Main_USA_GasCurve2,0618_NatGas_USA_GasCurve2.pptx
  • Slightly less gas
• GSw_GasPriceAdjMethod = 2:
  • results-0618_Main_USA_GasCurve2,0618_NatGas_USA_GasCurve2_DailyAdj.pptx
  • Slightly more gas
• GSw_GasPriceAdjMethod = 0:
  • results-0618_Main_USA_GasCurve2,0618_NatGas_USA_GasCurve2_NoAdj.pptx
  • Slightly less gas

GSw_GasCurve = 3

• GSw_GasPriceAdjMethod = 1:
  • results-0618_Main_USA_GasCurve3,0618_NatGas_USA_GasCurve3.pptx
  • Virtually no change
• GSw_GasPriceAdjMethod = 2:
  • results-0618_Main_USA_GasCurve3,0618_NatGas_USA_GasCurve3_DailyAdj.pptx
  • Very minor changes
• GSw_GasPriceAdjMethod = 0:
  • results-0618_Main_USA_GasCurve3,0618_NatGas_USA_GasCurve3_NoAdj.pptx
  • Very minor changes

GSw_GasCurve = 0

• GSw_GasPriceAdjMethod = 1:
  • results-0618_Main_USA_GasCurve0,0618_NatGas_USA_GasCurve0.pptx
  • Slightly more gas
• GSw_GasPriceAdjMethod = 2:
  • results-0618_Main_USA_GasCurve0,0618_NatGas_USA_GasCurve0_DailyAdj.pptx
  • Similar gas buildout to GSw_GasCurve = 1 with slightly more gas generation, but still not a huge amount
• GSw_GasPriceAdjMethod = 0:
  • results-0618_Main_USA_GasCurve0,0618_NatGas_USA_GasCurve0_NoAdj.pptx
  • Minor changes

The 2050 cendiv-level daily multipliers produced by these runs, compared with the new seasonal adjustment, are below.

Checklist for author

Details to double-check

• Charge code provided to reviewers
• Included comparison reports for appropriate test cases
• Documentation updated if necessary
• If input data added/modified:
  • Units are specified
  • Preprocessing steps have been documented and committed to ReEDS_Input_Processing: https://github.com/ReEDS-Model/ReEDS_Input_Processing/tree/ko/temp_gasprice_regression/aeo_updates/temperature_gas_price_regression
• Code formatting standardized
• Reusable functions used where possible instead of copy/pasted code
• Verify all the temporal tests from cases_test.csv still work
  • Pacific_CC
  • Pacific_weks
  • Pacific_full_year
  • Interday_storage
  • Pacific_rep15

General information to guide review

• Zero impact on results of default case
• No large data file(s) added/modified
• No substantive impact on runtime for full-US reference case
• No substantive impact on folder size for full-US reference case
• No change to process flow (runbatch.py, d_solve_iterate.py)
• No change to code organization
• No change to package requirements (environment.yml or Project.toml)

Did you use LLM tools (chatbot or copilot) in the preparation of this PR? If so, describe how

• LLMs used for timestamp compatibility between different outputs to ensure consistency in the final csv outputted
• LLMs used to run through the files and find relevant predefined functions I (Lara) may have originally missed
• LLMs used to debug h5 file handling

Tag points of contact here if you would like additional review of the relevant parts of the model

[kodiobika]

Thanks for your work on this! It's looking great so far -- mostly just some stylistic suggestions

[wesleyjcole]

I noticed that the z error_check went down when you turned on the daily natural gas price adjustments with both GSw_GasCurve = 1 and GSw_GasCurve = 2. I couldn't see any reasons why that would change. Any thoughts as to why that might be? It moved in a good direction, so if you don't have ideas for why it changed, that's fine.

[kodiobika]

I noticed that the z error_check went down when you turned on the daily natural gas price adjustments with both GSw_GasCurve = 1 and GSw_GasCurve = 2. I couldn't see any reasons why that would change. Any thoughts as to why that might be? It moved in a good direction, so if you don't have ideas for why it changed, that's fine.

Weird, yeah I'm not sure why it changed

[patrickbrown4]

Can you note the default setting?

[patrickbrown4]

Have you shown that these two methods match for the cendiv regions? I worry a bit about using different methods for different regions - personally my preference would be to either use cendiv (if that's the resolution at which you have the AEO data) or use the single NOAA method everywhere.

[kodiobika]

They differ slightly due to AEO using different (paywalled) population projections and AEO's version of the Pacific cendiv including AK and HI. But good point, just using one method would be better for consistency and easier to explain - I'll update that to use just the NOAA method

[patrickbrown4]

To make sure I understand, is the projection step to capture climate change? Is a particular climate scenario assumed (if so, it'd be good to note it)?

[kodiobika]

That's right, specifically we wanted our HDD/CDD baseline for each year to match what's in the AEO. So I guess the scenario would just be AEO 2026. I'll update the documentation to note that that's where the projection method comes from

[kodiobika]

Although now that I think about it, we initially wanted to make sure the HDD/CDD baselines match because we were regressing daily HDD/CDDs on daily prices, and the projected HDD/CDDs are baked into the AEO prices. But now we're regressing daily HDD/CDDs on the deviation of daily prices from an annual average (i.e., on daily multipliers), so I think we can just remove the projections entirely and just apply the regression coefficients to the historical (weather year) daily degree days. @wesleyjcole Does this sound right to you?

[kodiobika]

Actually maybe not, thinking about it more - if 2050 is a hotter year for a region than 2025, we would want its price multipliers to be higher on average in 2050, even if the hourly weather patterns are the same right?

[kodiobika]

Thinking about it more (again), if the HDD/CDD projections are already baked into the year-round average AEO prices, then actually we might be double-counting the hotness or coldness of the model years (maybe the normalization of the multipliers makes it all moot anyway though?). So I think we should just be using the historical degree days

[wesleyjcole]

Thanks for putting your thinking in this thread. This is complicated, and I think I can convince myself that either approach is reasonable. In the end, I don't think it will matter that much, and it would be much simpler (and easier to maintain) if we didn't use the projections of the HDD/CDD, so I'm for taking your suggestion and just use the historical degree days.

[patrickbrown4]

Yeah, I'm not sure I follow everything, but that sounds most in keeping with the current approach of using historical weather for wind/solar/load/outages. Eventually it'd be great to have a unified approach for incorporating defined climate projections into all of these, so I like the idea of keeping the gas price treatment more straightforwardly historical for now, and then thinking more about how to update them all together.

[patrickbrown4]

It might be helpful to note the default treatment (the one used in the USA_defaults case)

[patrickbrown4]

Do you think it'd be worthwhile to add a test case for the new daily adjustment method to cases_test.csv? If there's a lot of processing that only runs when GSw_GasPriceAdjMethod = 2, having a test case might help to keep it working. If the differences are relatively minor and/or you don't think the feature will be used often, it's ok to leave it out.

[kodiobika]

The daily adjustment processing gets run in every case and then the multipliers are overwritten if GSw_GasPriceAdjMethod != 2, so I think we should be okay

[patrickbrown4]

If you use multiple representative years (e.g. GSw_HourlyWeatherYears=2012_2013), it looks like you'll get different average prices in the different weather years (which is what I would want to happen, since it does actually happen) - is that right?

[kodiobika]

Yep that's right - the average price within a model year will be the same, but the average price across different weather years can differ

[patrickbrown4]

I checked the versions of daily_gasprice_multipliers_cendiv.csv and daily_gasprice_multipliers_r.csv saved to the inputs_case/stress{year}i0 folders, and they're all empty. I think that means we multiply gas prices by zero during stress periods, which will mess up the dispatch order during stress periods.

So I think you want resource_adequacy_years here instead of GSw_HourlyWeatherYears. I'm not sure if the downstream code in hourly_writetimeseries.py will also need to be adjusted.

[patrickbrown4]

If I'm reading it right, it seems like gasprice_adj_cendiv won't have an effect here, since it's only used in the conditional on gb. Should it multiply gasprice instead?

Although actually, repgasprice is only used to calculate repgasprice_filt, which used to be used in Augur (probably for the storage arbitrage value calculation) but is no longer used. So this whole block should be able to be deleted. (Or can save for a later cleanup PR, either works.)

[patrickbrown4]

Just to double check (I'm not sure I'm interpreting correctly as I don't think I've dealt with this equation before), I don't think I follow how gasprice_adj_cendiv is having an effect here, as it shouldn't affect the gb conditional. Should there be a GASUSED.l * gasprice_adj_cendiv * gasprice term?

[kodiobika]

Ah yeah sorry, I think I was just looking for all gas price-related terms that summed over h and added the adjustment term there without thinking. Looking at it carefully now I think what this statement is doing is finding the max gasprice among gas bins where the total amount of gas used across all hours is non-zero? So I think that sum statement should remain unchanged (from main) and instead the gasprice term should be replaced by gasprice(cendiv,gb,t) * sum{h, gasprice_adj_cendiv(cendiv,h,t)}. Does that sound right?

[kodiobika]

Oh actually we don't want to multiply by the sum there, I think we'd want to find the hour that maximizes the gas price, so instead it'd be gasprice(cendiv,gb,t) * smax{h, gasprice_adj_cendiv(cendiv,h,t)}?

[kodiobika]

But then again, I guess we'd only want to consider hours where GASUSED.l is non-zero (assuming my interpretation is right). And actually I'm not sure if we want repgasprice to include the adjustment terms if it's just supposed to be "representative", so maybe it should just be reverted to what's in main. @wesleyjcole Any opinion on this one?

[wesleyjcole]

Patrick - thanks for catching this one.

Kodi - I think we need to calculate an annual gasprice_adj_cendiv value, weighted by sum{gb, GASUSED.l(cendiv,gb,h,t) }. This would give a weighted average price adjustment based on when gased in used. I expect that it will be close to one, given that the adjustment was built to not change the annual input price.