[1]:
from pathlib import Path
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("ticks")
SMALL_SIZE = 12
MEDIUM_SIZE = 14
BIGGER_SIZE = 16
plt.rc('font', size=SMALL_SIZE) # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
plt.rc('figure', titlesize=22) # fontsize of the figure title
plt.rc('axes', titlesize=20) # fontsize of the figure title
[2]:
mws = np.linspace(0, 227.05, 10)
m_cost = [0 if mw < (227.05 - 127.05) else 15 for mw in mws]
plt.step(mws, m_cost)
plt.xlabel("MWh")
plt.ylabel("Marginal Cost [$/MWh]")
plt.title("Wind-PEM Marginal Cost Curve at T=t")
[2]:
Text(0.5, 1.0, 'Wind-PEM Marginal Cost Curve at T=t')
[3]:
import glob
import numpy as np
from dispatches.case_studies.renewables_case.double_loop_utils import read_rts_gmlc_wind_inputs
from dispatches_sample_data import rts_gmlc
from dispatches.case_studies.renewables_case.load_parameters import wind_cap_cost, wind_op_cost, pem_cap_cost, pem_op_cost, pem_var_cost, PA
from dispatches_data.api import files
wind_cf = read_rts_gmlc_wind_inputs(rts_gmlc.source_data_path, gen_name="303_WIND_1", agg_func="first")['303_WIND_1-RTCF']
results_file_pattern = f"RE/input_data/sweep_parameters_15_1000.csv"
params_file = files("dynamic_sweep", pattern=results_file_pattern)
if not params_file:
raise LookupError(f"No files found with pattern {results_file_pattern}")
params = pd.read_csv(params_file[0])
records = []
csv_files_to_load = files("dynamic_sweep", pattern=f"RE/results_parameter_sweep_15_1000/*.csv")
for filename in csv_files_to_load:
res = pd.read_csv(filename)
ind = int(Path(filename).stem.split('_')[-1])
pem_bid = params['PEM_bid'][ind]
pem_size = params['PEM_power_capacity'][ind]
e_revenue = (res["Dispatch DA"] * res["LMP DA"] + (res["Dispatch"] - res["Dispatch DA"]) * res["LMP"]).sum()
pem_elec = np.clip(wind_cf.values * 847 - res['Dispatch'], 0, pem_size)
h2_price = 3
h_revenue = pem_elec.sum() / 54.953 * h2_price * 1e3
cap_cost = (847 * wind_cap_cost + pem_size * pem_cap_cost) * 1e3
fixed_op_cost = (847 * wind_op_cost + pem_size * pem_op_cost) * 1e3
var_op_cost = pem_var_cost * pem_elec.sum()
npv = -cap_cost + PA * (e_revenue + h_revenue - fixed_op_cost - var_op_cost)
records.append({
"e_revenue": e_revenue,
"h_revenue": h_revenue,
"pem_bid": pem_bid,
"pem_mw": pem_size,
'NPV': npv})
sweep_results = pd.DataFrame(records)
sweep_results['pem_bid_round'] = sweep_results['pem_bid'].astype('int')
sweep_results['pem_mw_round'] = sweep_results['pem_mw'].astype('int')
sweep_results['e_revenue_mil'] = (sweep_results['e_revenue'] * 1e-6)
sweep_results['h_revenue_mil'] = (sweep_results['h_revenue'] * 1e-6)
sweep_results['NPV_bil'] = (sweep_results['NPV'] * 1e-9)
[4]:
# get results file from RE_surrogate_optimization_steadystate.py
re_case_dir = Path.cwd()
results = pd.read_csv(re_case_dir / "surrogate_results_ss.csv")
results = results.drop(columns=["Unnamed: 0", 'wind_mw'])
results['pem_bid_round'] = results['pem_bid'].astype('int')
results['pem_mw_round'] = results['pem_mw'].astype('int')
results['e_revenue_mil'] = (results['e_revenue'] * 1e-6)
results['h_revenue_mil'] = (results['h_revenue'] * 1e-6)
results['NPV_bil'] = (results['NPV'] * 1e-9)
results
[4]:
| pem_mw | pem_bid | e_revenue | h_revenue | NPV | freq_day_0 | freq_day_1 | freq_day_2 | freq_day_3 | freq_day_4 | ... | freq_day_15 | freq_day_16 | freq_day_17 | freq_day_18 | freq_day_19 | pem_bid_round | pem_mw_round | e_revenue_mil | h_revenue_mil | NPV_bil | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 127.050 | 15.0 | -8.370107e+06 | 2.129433e+07 | -1.552816e+09 | 0.028895 | 0.063540 | 0.427766 | 0.003281 | 0.000458 | ... | 0.123763 | 0.023888 | 0.005713 | 0.034070 | 0.073425 | 15 | 127 | -8.370107 | 21.294334 | -1.552816 |
| 1 | 148.225 | 15.0 | -7.738664e+06 | 2.342615e+07 | -1.561226e+09 | 0.025751 | 0.056873 | 0.427467 | 0.004410 | 0.000105 | ... | 0.124810 | 0.024191 | 0.012067 | 0.034141 | 0.073441 | 15 | 148 | -7.738664 | 23.426149 | -1.561226 |
| 2 | 169.400 | 15.0 | -6.992188e+06 | 2.584096e+07 | -1.563620e+09 | 0.022988 | 0.049106 | 0.426976 | 0.005417 | 0.000059 | ... | 0.125266 | 0.024992 | 0.018411 | 0.034165 | 0.073034 | 15 | 169 | -6.992188 | 25.840958 | -1.563620 |
| 3 | 190.575 | 15.0 | -6.139070e+06 | 2.721534e+07 | -1.579989e+09 | 0.021001 | 0.041341 | 0.426691 | 0.005945 | 0.000062 | ... | 0.125027 | 0.025733 | 0.023153 | 0.033734 | 0.071100 | 15 | 190 | -6.139070 | 27.215340 | -1.579989 |
| 4 | 211.750 | 15.0 | -5.196320e+06 | 2.773554e+07 | -1.607797e+09 | 0.019871 | 0.033633 | 0.426775 | 0.005686 | 0.001641 | ... | 0.124345 | 0.026239 | 0.025169 | 0.032097 | 0.065010 | 15 | 211 | -5.196320 | 27.735539 | -1.607797 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 190 | 338.800 | 45.0 | 5.366495e+06 | 6.490788e+07 | -1.188242e+09 | 0.079991 | 0.004136 | 0.389319 | 0.049944 | 0.047801 | ... | 0.014834 | 0.120509 | 0.066808 | 0.000282 | 0.002565 | 45 | 338 | 5.366495 | 64.907877 | -1.188242 |
| 191 | 359.975 | 45.0 | 6.100078e+06 | 6.749057e+07 | -1.188315e+09 | 0.079942 | 0.004292 | 0.389276 | 0.055595 | 0.026800 | ... | 0.014218 | 0.120158 | 0.067412 | 0.000584 | 0.002580 | 45 | 359 | 6.100078 | 67.490574 | -1.188315 |
| 192 | 381.150 | 45.0 | 6.669175e+06 | 6.906862e+07 | -1.205993e+09 | 0.080035 | 0.004216 | 0.389166 | 0.059080 | 0.013766 | ... | 0.013994 | 0.119731 | 0.067629 | 0.000838 | 0.002450 | 45 | 381 | 6.669175 | 69.068617 | -1.205993 |
| 193 | 402.325 | 45.0 | 7.101820e+06 | 6.970906e+07 | -1.239839e+09 | 0.080240 | 0.004013 | 0.389038 | 0.061097 | 0.006632 | ... | 0.014017 | 0.119420 | 0.067614 | 0.000902 | 0.002229 | 45 | 402 | 7.101820 | 69.709057 | -1.239839 |
| 194 | 423.500 | 45.0 | 7.424840e+06 | 7.002385e+07 | -1.280257e+09 | 0.080448 | 0.003809 | 0.388930 | 0.062328 | 0.002577 | ... | 0.014115 | 0.119213 | 0.067537 | 0.000888 | 0.002023 | 45 | 423 | 7.424840 | 70.023854 | -1.280257 |
195 rows × 30 columns
[5]:
elec_rev_min = min(results['e_revenue_mil'].min(), sweep_results['e_revenue_mil'].min())
elec_rev_max = max(results['e_revenue_mil'].max(), sweep_results['e_revenue_mil'].max())
h2_rev_min = min(results['h_revenue_mil'].min(), sweep_results['h_revenue_mil'].min())
h2_rev_max = max(results['h_revenue_mil'].max(), sweep_results['h_revenue_mil'].max())
npv_min = min(results['NPV_bil'].min(), sweep_results['NPV_bil'].min())
npv_max = max(results['NPV_bil'].max(), sweep_results['NPV_bil'].max())
[6]:
common_bids = set(results['pem_bid'].unique()).intersection(set(sweep_results['pem_bid'].unique()))
common_sizes = set(results['pem_mw'].round(2).unique()).intersection(set(sweep_results['pem_mw'].round(2).unique()))
results = results[(results['pem_bid'].isin(common_bids)) & (results['pem_mw'].round(2).isin(common_sizes))]
fig, axs = plt.subplots(1, 3, figsize=(14, 4))
sns.heatmap(results.pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values='e_revenue_mil',
aggfunc='mean'), annot=True, fmt=".1f", ax=axs[0],
vmin=elec_rev_min, vmax=elec_rev_max)
axs[0].set_title("Elec Revenue [$Mil]")
axs[0].set_xlabel("PEM Size [MW]")
axs[0].set_ylabel("PEM Bid [$/MW]")
sns.heatmap(results.pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values='h_revenue_mil',
aggfunc='mean'), annot=True, fmt=".0f", ax=axs[1],
vmin=h2_rev_min, vmax=h2_rev_max)
axs[1].set_title("H2 Revenue [$Mil]")
axs[1].set_xlabel("PEM Size [MW]")
axs[1].set_ylabel("PEM Bid [$/MW]")
sns.heatmap(results.pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values='NPV_bil',
aggfunc='mean'), annot=True, fmt=".1f", ax=axs[2],
vmin=npv_min, vmax=npv_max)
axs[2].set_title("NPV [$Bil]")
axs[2].set_xlabel("PEM Size [MW]")
axs[2].set_ylabel("PEM Bid [$/MW]")
for ax in axs:
ax.invert_yaxis()
plt.suptitle("Results from Surrogates")
plt.tight_layout()
Compare with Input Dataset from Prescient Sweep
[7]:
fig, axs = plt.subplots(1, 3, figsize=(14, 4))
sns.heatmap(sweep_results.pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values='e_revenue_mil',
aggfunc='mean'), annot=True, fmt=".1f", ax=axs[0],
vmin=elec_rev_min, vmax=elec_rev_max)
axs[0].set_title("Elec Revenue [$Mil]")
axs[0].set_xlabel("PEM Size [MW]")
axs[0].set_ylabel("PEM Bid [$/MW]")
sns.heatmap(sweep_results.pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values='h_revenue_mil',
aggfunc='mean'), annot=True, fmt=".0f", ax=axs[1],
vmin=h2_rev_min, vmax=h2_rev_max)
axs[1].set_title("H2 Revenue [$Mil]")
axs[1].set_xlabel("PEM Size [MW]")
axs[1].set_ylabel("PEM Bid [$/MW]")
sns.heatmap(sweep_results.pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values='NPV_bil',
aggfunc='mean'), annot=True, fmt=".1f", ax=axs[2],
vmin=npv_min, vmax=npv_max)
axs[2].set_title("NPV [$Bil]")
axs[2].set_xlabel("PEM Size [MW]")
axs[2].set_ylabel("PEM Bid [$/MW]")
for ax in axs:
ax.invert_yaxis()
plt.suptitle("Results from Sweep")
plt.tight_layout()
[8]:
def calc_diff_df(value):
sweep_diff = sweep_results[(sweep_results['pem_bid'].isin(common_bids)) & (sweep_results['pem_mw'].isin(common_sizes))].pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values=value,
aggfunc='mean')
diff_df = (results[(results['pem_bid'].isin(common_bids)) & (results['pem_mw'].isin(common_sizes))].pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values=value,
aggfunc='mean') - sweep_diff) / sweep_diff
return diff_df * 1e2
fig, axs = plt.subplots(1, 3, figsize=(14, 4))
sns.heatmap(calc_diff_df('e_revenue_mil'), annot=True, fmt=".0f", ax=axs[0])
axs[0].set_title("Elec Revenue [%]")
axs[0].set_xlabel("PEM Size [MW]")
axs[0].set_ylabel("PEM Bid [$/MW]")
sns.heatmap(calc_diff_df('h_revenue_mil'), annot=True, fmt=".0f", ax=axs[1])
axs[1].set_title("H2 Revenue [%]")
axs[1].set_xlabel("PEM Size [MW]")
axs[1].set_ylabel("PEM Bid [$/MW]")
sns.heatmap(calc_diff_df('NPV_bil'), annot=True, fmt=".0f", ax=axs[2])
axs[2].set_title("NPV [%]")
axs[2].set_xlabel("PEM Size [MW]")
axs[2].set_ylabel("PEM Bid [$/MW]")
for ax in axs:
ax.invert_yaxis()
plt.suptitle("Percent Difference: (Double Loop - Sweep) / Sweep")
plt.tight_layout()
[9]:
def calc_diff_df(value):
sweep_diff = sweep_results[(sweep_results['pem_bid'].isin(common_bids)) & (sweep_results['pem_mw'].isin(common_sizes))].pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values=value,
aggfunc='mean')
diff_df = (results[(results['pem_bid'].isin(common_bids)) & (results['pem_mw'].isin(common_sizes))].pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values=value,
aggfunc='mean') - sweep_diff)
print(min(diff_df.min()), min(diff_df.max()))
return diff_df
fig, axs = plt.subplots(1, 3, figsize=(14, 4))
sns.heatmap(calc_diff_df('e_revenue_mil'), annot=True, fmt=".0f", ax=axs[0])
axs[0].set_title("Elec Revenue")
axs[0].set_xlabel("PEM Size [MW]")
axs[0].set_ylabel("PEM Bid [$/MW]")
sns.heatmap(calc_diff_df('h_revenue_mil'), annot=True, fmt=".0f", ax=axs[1])
axs[1].set_title("H2 Revenue")
axs[1].set_xlabel("PEM Size [MW]")
axs[1].set_ylabel("PEM Bid [$/MW]")
sns.heatmap(calc_diff_df('NPV_bil'), annot=True, fmt=".0f", ax=axs[2])
axs[2].set_title("NPV")
axs[2].set_xlabel("PEM Size [MW]")
axs[2].set_ylabel("PEM Bid [$/MW]")
for ax in axs:
ax.invert_yaxis()
plt.suptitle("Difference: (Double Loop - Sweep)")
plt.tight_layout()
-3.376360794079675 0.8547649441334917
-8.450876253965717 -1.5617455661226742
-0.14683758471369557 -0.03333816706580417
Compare day frequencies
[10]:
for day in range(20):
results[f'percent_day_{day}'] = results[f'freq_day_{day}'] * 1e2
fig, axs = plt.subplots(1, 1, figsize=(8, 5))
sns.heatmap(results.pivot_table(index=f'pem_bid_round',
columns=f"pem_mw_round",
values=f'percent_day_{day}',
aggfunc='mean'), annot=True, fmt=".1f", ax=axs)
axs.set_title(f"Percent of Day {day}")
axs.set_xlabel("PEM Size [MW]")
axs.set_ylabel("PEM Bid [$/MW]")
axs.invert_yaxis()
/var/folders/yh/xvlf1svn0xq6r2l4hzsggz2xt3qfdr/T/ipykernel_73964/1946447195.py:2: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
results[f'percent_day_{day}'] = results[f'freq_day_{day}'] * 1e2
