Add price-weighted capacity credit output

[wesleyjcole]

Summary

This pull request builds on #171 to add a price-weighted variant of the capacity-credit-like plots for using stress periods. We have an option that looks at the top 10 hours, but this instead looks at the all hours and weights them by their price. For regions that have less than 10 hours with a meaningful price, this will mean the ~$0 periods won't get included. Here's an example from Alabama for 2050, which only has 7 higher priced hours (the % values show how much of the price weight occurs on that day):

For other regions that have more than 10 hours, this allows all of those hours get captured. For example, here is the CA_LA region for 2050:

It also changes the function's default to this new method because I think it is the most accurate way to estimate a marginal capacity credit.

Technical details

Implementation notes

I added a filter for infinity values in the capacity-credit division (gen / cap). There is sometimes a region with zero capacity but a tiny nonzero residual in generation (~1e-16), which creates an infinity value. This is was present for any metric, and I don't think it impacted them, but filtering it out means it won't cause unintended issues.

I also changed the ymax value of the plots to be dynamic because the capacity credit can go above 100% for some technologies (we measure capacity credit using net summer capacity, and net winter capacity can be meaningfully higher than net summer).

Validation, testing, and comparison report(s)

The output for this new metric looks like this (drawn from the run performed for #176):

Checklist for author

Details to double-check

• Charge code provided to reviewers
• [ ] Included comparison reports for appropriate test cases
• [ ] Documentation updated if necessary
• Code formatting standardized
• Reusable functions used where possible instead of copy/pasted code

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 (runreeds.py, reeds/core/solve/solve.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

Yes — I implemented this with Claude (I used plan mode for design, I asked it to stop for my review after it completed each step, and I used to for tracing for the inf values).

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

[patrickbrown4]

Thanks Wesley, sorry for the slow response.

[patrickbrown4]

How about changing this to stress_weight_price? Then it would match the format of the other inputs, which are of the form stress_{how to select}_{which parameter to use to select}. (Would need to update the calls elsewhere too.)

[patrickbrown4]

Simpler alternative would be to default to the existing ymax if it's over 100:

Suggested change

	finite_vals = np.concatenate([
	capcredit.values[np.isfinite(capcredit.values)],
	repfraction.values[np.isfinite(repfraction.values)],
	])
	datamax = finite_vals.max() if finite_vals.size else 0
	ymax = datamax * 1.02 if datamax > 100 else 100
	_ax.set_ylim(0, ymax)
	_ax.set_ylim(0, max(100, _ax.get_ylim()[1]))

[patrickbrown4]

The default doesn't really matter since we always specify, but could you say more about why you think this method is the most accurate way to estimate a marginal capacity credit? I understand that it reflects the hours that drive capacity revenues, but it seems the most different from how we've calculated capacity credit in the past.

My other concern is about stability; if it's based on fewer hours, it seems like it could jump around more between years/regions/runs and make it hard to tell trends from noise. For example, I'm not sure how "real" the low values for storage in Texas from 2032-2038 are:

None of these appraoches are perfect because we're showing dispatch and not availability. So I guess it's more about whether we want a metric that's more similar to the way we calculate capacity credit (top 10 net load or top 10 load), or one that better captures the timing of the price spikes (price weighted).

[patrickbrown4]

Reshaping big dataframes with r/h indices can be slow (maybe not for the default case but for high-temporal/spatial-resolution cases). You can sum by t first and avoid the big reshape:

Suggested change

	price_long = price_stress.stack('r').rename('price').reset_index()
	price_sum = price_long.groupby(['t','r'], as_index=False).price.sum()
	price_sum = price_stress.groupby('t').sum().stack('r')

[patrickbrown4]

For the suggested change below (avoid reshaping)

Suggested change

	price_stress = price_stress.groupby(['t','r','h']).Value.max()

[patrickbrown4]

For the suggested change below. Import it here to localize the special-case dependence

Suggested change

	from reeds.core.terminus import report_calcs
	### Check inputs

[patrickbrown4]

For the suggested change below

Suggested change

	hours_t = (
	report_calcs.get_gams_results(case, 'hours_t')['hours_t'].reset_index()
	.astype({'t':int}).rename(columns={'allh':'h'}).set_index(['h','t']).squeeze(1)
	)

[patrickbrown4]

Avoid reshaping

Suggested change

	dfindex = price_stress.reorder_levels(['t','r','h'])

[patrickbrown4]

If we change the weighting for stress hours (discussed in #154), summing over hours directly will be wrong. Safer to first multiply by the hour weight in hours_t(allh,t). Here's my suggestion:

Suggested change

	## sum_h(gen*price) / sum_h(price)
	price_long = price_stress.stack('r').rename('price').reset_index()
	price_sum = price_long.groupby(['t','r'], as_index=False).price.sum()
	gen_price = gen_h_stress.merge(price_long, on=['t','r','h'], how='left')
	gen_price['price'] = gen_price['price'].fillna(0)
	gen_price['gen_x_price'] = gen_price.MW * gen_price.price
	numer = gen_price.groupby(['t','i','r'], as_index=False).gen_x_price.sum()
	df = numer.merge(price_sum, on=['t','r'], how='left')
	df['MW'] = df.gen_x_price / df.price
	df = df.set_index(['t','i','r']).MW.unstack('r')
	## Guard against exact-cancellation division by zero (positive/negative hourly
	## prices summing to ~0 while gen_x_price is nonzero); NaN (0/0) is also possible
	## and both should collapse to 0
	df = df.replace([np.inf, -np.inf], np.nan).fillna(0)
	## sum_h(gen * price * hours) / sum_h(price * hours)
	price_weighted = price_stress.multiply(hours_t, axis=0).dropna()
	df = (
	(gen_h_stress.set_index(['t','i','r','h']).squeeze(1) * price_weighted)
	.groupby(['t','i','r']).sum()
	/ price_weighted.groupby(['t','r']).sum()
	).unstack('r').replace([np.inf, -np.inf], np.nan).fillna(0)

I'll note the other changes above/below

[patrickbrown4]

Let's match the same options used below to facilitate comparison (not sure if the suggestion will apply correctly since it's outside the changed region but you get the idea)

Suggested change

	metrics = [
	'cap',
	'rep_mean',
	'stress_mean',
	'stress_top10_load',
	'stress_top10_netload',
	'stress_bottom10_vregen',
	'stress_top10_price',
	'stress_price_weighted',
	]

[patrickbrown4]

Sorry, the one other required change is to get_gams_results() in report_calcs.py; it should now be:

def get_gams_results(case, key=None):
    print('Loading results.gdx')
    if isinstance(key, str):
        dictin = {key: gdxpds.to_dataframe(Path(case, 'outputs', 'results.gdx'), key)}
    else:
        dictin = gdxpds.to_dataframes(Path(case, 'outputs', 'results.gdx'))
    ## Set indices as multiindex
    valcols = ['Value','Level','Marginal','Lower','Upper','Scale']
    for key, df in dictin.items():
        indices = [i for i in df if i not in valcols]
        dictin[key] = df.set_index(indices).squeeze(1)
    print('Finished loading results.gdx')
    return dictin