Skip to content

Commit 3919b35

Browse files
staadeckerstaadecker
committed
Create LMP figure
1 parent 165adee commit 3919b35

File tree

3 files changed

+226
-98
lines changed

3 files changed

+226
-98
lines changed

‎papers/Martin_Staadecker_Value_of_LDES_and_Factors/LDES_paper_graphs/Figure Costs.py‎

Lines changed: 184 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,184 @@
1+
# %% GET TOOLS
2+
3+
# Imports
4+
import matplotlib.pyplot as plt
5+
import numpy as np
6+
import pandas as pd
7+
from matplotlib import gridspec
8+
9+
from papers.Martin_Staadecker_Value_of_LDES_and_Factors.LDES_paper_graphs.util import (
10+
set_style,
11+
get_set_e_scenarios,
12+
)
13+
from switch_model.tools.graph.main import GraphTools
14+
15+
# Prepare graph tools
16+
tools = GraphTools(scenarios=get_set_e_scenarios())
17+
tools.pre_graphing(multi_scenario=True)
18+
raw_load_balance = tools.get_dataframe(
19+
"load_balance.csv",
20+
usecols=[
21+
"load_zone",
22+
"timestamp",
23+
"normalized_energy_balance_duals_dollar_per_mwh",
24+
"scenario_name",
25+
],
26+
).rename(columns={"normalized_energy_balance_duals_dollar_per_mwh": "value"})
27+
raw_load_balance.value *= 0.1 # Convert from $/MWh to c/kWh
28+
raw_load_balance = tools.transform.load_zone(raw_load_balance)
29+
30+
# %% CREATE FIGURE
31+
set_style()
32+
plt.close()
33+
fig = plt.figure()
34+
fig.set_size_inches(12, 12)
35+
gs = gridspec.GridSpec(2, 2, figure=fig)
36+
ax1 = fig.add_subplot(gs[0, 0])
37+
ax2 = fig.add_subplot(gs[0, 1])
38+
ax3 = fig.add_subplot(gs[1, 0])
39+
ax4 = fig.add_subplot(gs[1, 1])
40+
41+
y_label = u"Mean estimated LMP (\xa2/kWh)"
42+
x_label = "WECC-wide storage capacity (TWh)"
43+
44+
# %% Variability
45+
46+
ax = ax1
47+
ax.clear()
48+
duals = raw_load_balance.copy()
49+
# duals = duals[duals.region == "WA"] # uncomment to filter by region
50+
duals = duals.groupby("scenario_name").value
51+
variability = pd.DataFrame()
52+
quantiles = [0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99]
53+
for quantile in quantiles:
54+
variability.loc[:, quantile] = duals.quantile(quantile)
55+
56+
cmap = plt.cm.get_cmap("viridis")
57+
58+
ax.fill_between(
59+
variability.index,
60+
variability[0.05],
61+
variability[0.95],
62+
alpha=0.3,
63+
color="dimgray",
64+
# marker=".",
65+
label="5th-95th quantile",
66+
)
67+
ax.fill_between(
68+
variability.index,
69+
variability[0.25],
70+
variability[0.75],
71+
alpha=1,
72+
color="dimgray",
73+
# marker=".",
74+
label="25th-75th quantile",
75+
)
76+
ax.plot(
77+
variability.index,
78+
variability[0.01],
79+
marker=".",
80+
color="dimgray",
81+
label="1st & 99th quantile",
82+
)
83+
ax.plot(variability.index, variability[0.99], marker=".", color="dimgray")
84+
ax.plot(variability.index, variability[0.5], marker=".", color="red", label="Median")
85+
86+
ax.set_xlabel(x_label)
87+
88+
ax.set_ylabel(y_label)
89+
ax.legend()
90+
ax.set_title("A. Distribution of estimated LMPs")
91+
92+
# %% Daily LMP
93+
94+
ax = ax3
95+
ax.clear()
96+
time_mapping = {0: "Midnight", 4: "4am", 8: "8am", 12: "Noon", 16: "4pm", 20: "8pm"}
97+
daily_lmp = raw_load_balance.copy()
98+
daily_lmp = tools.transform.timestamp(daily_lmp)
99+
daily_lmp = daily_lmp.groupby(["scenario_name", "hour"], as_index=False).value.mean()
100+
daily_lmp = daily_lmp.pivot(index="scenario_name", columns="hour", values="value")
101+
daily_lmp = daily_lmp.rename(time_mapping, axis=1).rename_axis(
102+
"Time of day (PST)", axis=1
103+
)
104+
daily_lmp = daily_lmp.sort_values(by=1.94, axis=1, ascending=False)
105+
daily_lmp.plot(
106+
ax=ax,
107+
xlabel=x_label,
108+
marker=".",
109+
ylabel=y_label,
110+
# cmap="tab10"
111+
)
112+
ax.set_title("C. Estimated LMP by time of day")
113+
# %% YEARLY LMP
114+
ax = ax4
115+
ax.clear()
116+
months_map = {
117+
1: "Jan",
118+
2: "Feb-May",
119+
3: "Feb-May",
120+
4: "Feb-May",
121+
5: "Feb-May",
122+
6: "Jun",
123+
7: "Jul",
124+
8: "Aug",
125+
9: "Sep-Oct",
126+
10: "Sep-Oct",
127+
11: "Nov",
128+
12: "Dec",
129+
}
130+
131+
cap = raw_load_balance
132+
cap = cap.groupby(["scenario_name", "timestamp"], as_index=False).value.mean()
133+
cap = tools.transform.timestamp(cap)
134+
cap = cap.set_index("datetime")
135+
cap["month"] = cap.index.month
136+
cap["month"] = cap.month.map(months_map)
137+
cap = cap.groupby(["scenario_name", "month"]).value.mean()
138+
cap = cap.unstack("month").rename_axis("Months", axis=1)
139+
cap = cap.sort_values(by=1.94, ascending=False, axis=1)
140+
# cap = cap[["Jan", "Feb-May", "Jun", "Jul", "Aug", "Sep-Oct", "Nov", "Dec"]]
141+
142+
cap.plot(
143+
ax=ax,
144+
xlabel=x_label,
145+
marker=".",
146+
ylabel=y_label,
147+
# cmap="tab10"
148+
)
149+
ax.set_title("D. Estimated LMP by time of year")
150+
# %% GEOGRAPHICAL LMP
151+
ax = ax2
152+
ax.clear()
153+
154+
geo = raw_load_balance.copy()
155+
region_map = {
156+
"CAN": "Canada",
157+
"WA": "OR, WA",
158+
"OR": "OR, WA",
159+
"ID": "Idaho",
160+
"MT": "Montana",
161+
"WY": "CO, UT, WY",
162+
"CO": "CO, UT, WY",
163+
"UT": "CO, UT, WY",
164+
"NV": "Nevada",
165+
"CA": "California",
166+
"AZ": "AZ, NM, MEX",
167+
"NM": "AZ, NM, MEX",
168+
"MEX": "AZ, NM, MEX",
169+
}
170+
geo.region = geo.region.map(region_map)
171+
geo = geo.groupby(["scenario_name", "region"]).value.mean()
172+
geo = geo.unstack("region")
173+
geo = geo.sort_values(by=1.94, axis=1, ascending=False)
174+
geo = geo.rename_axis("Regions", axis=1)
175+
geo.plot(
176+
ax=ax,
177+
xlabel=x_label,
178+
marker=".",
179+
ylabel=y_label,
180+
# cmap="tab10"
181+
)
182+
ax.set_title("B. Estimated LMP by region")
183+
# %%
184+
fig.tight_layout()

‎papers/Martin_Staadecker_Value_of_LDES_and_Factors/LDES_paper_graphs/Figure impact of LDES on grid.py‎

Lines changed: 24 additions & 98 deletions
Original file line numberDiff line numberDiff line change
@@ -6,31 +6,14 @@
66
from matplotlib.ticker import PercentFormatter
77

88
from papers.Martin_Staadecker_Value_of_LDES_and_Factors.LDES_paper_graphs.util import (
9-
get_scenario,
109
set_style,
10+
get_set_e_scenarios,
1111
)
1212
from switch_model.tools.graph.main import GraphTools
1313

1414

1515
# Prepare graph tools
16-
tools = GraphTools(
17-
scenarios=[
18-
get_scenario("1342", name=1.94),
19-
get_scenario("M7", name=2),
20-
get_scenario("M10", name=2.5),
21-
get_scenario("M9", name=3),
22-
get_scenario("M6", name=4),
23-
get_scenario("M5", name=8),
24-
get_scenario("M11", name=12),
25-
get_scenario("M4", name=16),
26-
get_scenario("M14", name=18),
27-
get_scenario("M13", name=20),
28-
get_scenario("M8", name=24),
29-
get_scenario("M3", name=32),
30-
get_scenario("M12", name=48),
31-
get_scenario("M2", name=64),
32-
]
33-
)
16+
tools = GraphTools(scenarios=get_set_e_scenarios())
3417
tools.pre_graphing(multi_scenario=True)
3518

3619
# %% CREATE FIGURE
@@ -43,76 +26,6 @@
4326
ax3 = fig.add_subplot(2, 2, 3)
4427
ax4 = fig.add_subplot(2, 2, 4)
4528

46-
# %% Daily LMP
47-
48-
ax = ax3
49-
ax.clear()
50-
time_mapping = {
51-
0: "Midnight",
52-
4: "4am",
53-
8: "8am",
54-
12: "Noon",
55-
16: "4pm",
56-
20: "8pm"
57-
}
58-
daily_lmp = tools.get_dataframe(
59-
"load_balance.csv",
60-
usecols=[
61-
"timestamp",
62-
"normalized_energy_balance_duals_dollar_per_mwh",
63-
"scenario_name",
64-
],
65-
).rename(columns={"normalized_energy_balance_duals_dollar_per_mwh": "value"})
66-
daily_lmp = tools.transform.timestamp(daily_lmp)
67-
daily_lmp = daily_lmp.groupby(["scenario_name", "hour"], as_index=False)["value"].mean()
68-
daily_lmp = daily_lmp.pivot(index="scenario_name", columns="hour", values="value")
69-
daily_lmp *= 0.1 # Convert from $/MWh to cents/kWh
70-
daily_lmp = daily_lmp.rename(time_mapping, axis=1).rename_axis("Time of day (PST)", axis=1)
71-
daily_lmp.plot(
72-
ax=ax,
73-
xlabel="WECC-wide storage capacity (TWh)",
74-
marker=".",
75-
ylabel=u"Estimated LMP (\xa2/kWh)",
76-
# cmap="tab10"
77-
)
78-
ax.set_title("C. Average estimated LMP by time of day")
79-
# %% YEARLY LMP
80-
ax = ax4
81-
ax.clear()
82-
scenarios_for_yearly_lmp = [5, 7, 8, 11, 12]
83-
84-
cap = tools.get_dataframe(
85-
"load_balance.csv",
86-
usecols=[
87-
"timestamp",
88-
"normalized_energy_balance_duals_dollar_per_mwh",
89-
"scenario_name",
90-
],
91-
).rename(columns={"normalized_energy_balance_duals_dollar_per_mwh": "value"})
92-
cap = cap.groupby(["scenario_name", "timestamp"], as_index=False).mean()
93-
cap = tools.transform.timestamp(cap)
94-
cap = cap.set_index("datetime")
95-
cap["month"] = cap.index.month
96-
cap = cap.groupby(["scenario_name", "month"]).value.mean()
97-
cap = cap.unstack("month").rename_axis("Storage Capacity (TWh)", axis=1)
98-
# Convert from $/MWh to cents/kWh
99-
cap *= 0.1
100-
cap = cap[scenarios_for_yearly_lmp]
101-
cap = cap.rename({
102-
5: "May",
103-
7: "July",
104-
8: "August",
105-
11: "November",
106-
12: "December",
107-
}, axis=1)
108-
cap.plot(
109-
ax=ax,
110-
xlabel="WECC-wide storage capacity (TWh)",
111-
marker=".",
112-
ylabel=u"Monthly-average estimated LMP (\xa2/kWh)",
113-
# cmap="tab10"
114-
)
115-
ax.set_title("D. Average estimated LMP during key months")
11629
# %% IMPACT ON TX AND GEN
11730

11831
ax = ax2
@@ -149,7 +62,9 @@
14962
# Filter out storage since it's not considered generation
15063
buildout = buildout[buildout["gen_type"] != "Storage"]
15164
# Sum across the entire scenario
152-
buildout = buildout.groupby("scenario_name")["BuildGen"].sum().rename("Built Generation")
65+
buildout = (
66+
buildout.groupby("scenario_name")["BuildGen"].sum().rename("Built Generation")
67+
)
15368

15469
cap = tools.get_dataframe(
15570
"gen_cap.csv",
@@ -166,7 +81,7 @@
16681
# Convert to percent against baseline
16782
df = (df / df.iloc[0] - 1) * 100
16883

169-
dotted_tx = df.loc[[1.94,3,20,64], ["Built Transmission"]]
84+
dotted_tx = df.loc[[1.94, 3, 20, 64], ["Built Transmission"]]
17085

17186
# Plot
17287
colors = tools.get_colors()
@@ -200,9 +115,15 @@
200115
# Add the gen_type column
201116
curtailment = tools.transform.gen_type(curtailment)
202117
# Convert to GW
203-
curtailment["value"] = curtailment["Curtailment_MW"] * curtailment["tp_weight_in_year_hrs"] / 1000
204-
curtailment = curtailment.groupby(["scenario_name", "gen_type"], as_index=False).value.sum()
205-
curtailment = curtailment.pivot(index="scenario_name", columns="gen_type", values="value")
118+
curtailment["value"] = (
119+
curtailment["Curtailment_MW"] * curtailment["tp_weight_in_year_hrs"] / 1000
120+
)
121+
curtailment = curtailment.groupby(
122+
["scenario_name", "gen_type"], as_index=False
123+
).value.sum()
124+
curtailment = curtailment.pivot(
125+
index="scenario_name", columns="gen_type", values="value"
126+
)
206127
curtailment /= 1000
207128
curtailment = curtailment.rename_axis("Technology", axis=1)
208129
curtailment.plot(ax=ax, color=tools.get_colors(), marker=".")
@@ -211,7 +132,9 @@
211132
ax.set_title("A. Impact of LDES on curtailment")
212133
ax.tick_params(top=False, bottom=False, right=False, left=False)
213134
# %%
214-
plt.subplots_adjust(hspace=0.2, wspace=0.25, left=0.07, right=0.97, top=0.95, bottom=0.07)
135+
plt.subplots_adjust(
136+
hspace=0.2, wspace=0.25, left=0.07, right=0.97, top=0.95, bottom=0.07
137+
)
215138

216139
# %% CALCULATIONS
217140
cap_total = cap["Solar"] + cap["Wind"]
@@ -224,8 +147,8 @@
224147
1 - cap / cap.iloc[0]
225148
cap - cap.iloc[0]
226149
# %% change from 20 to 64
227-
1 - cap / cap.loc[20] # wind decrease %
228-
cap / cap.loc[20] - 1 # solar increase %
150+
1 - cap / cap.loc[20] # wind decrease %
151+
cap / cap.loc[20] - 1 # solar increase %
229152
(cap - cap.loc[20]) / 1000
230153
# %% transmission
231154
tx / tx.iloc[0] * 100
@@ -237,4 +160,7 @@
237160
# %% average drop in daily duals
238161
df = daily_lmp.mean(axis=1)
239162
df / df.iloc[0] - 1
240-
daily_lmp / daily_lmp.iloc[0] - 1
163+
# daily_lmp / daily_lmp.iloc[0] - 1
164+
# %% night time drop
165+
df = daily_lmp[["Midnight", "8pm", "4am"]].mean(axis=1)
166+
(1 - df.loc[3] / df.iloc[0]) * 100 / ((3 - 1.94) * 10)

‎papers/Martin_Staadecker_Value_of_LDES_and_Factors/LDES_paper_graphs/util.py‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,3 +16,21 @@ def set_style():
1616
plt.rcParams["font.family"] = "sans-serif"
1717
plt.rcParams["ytick.minor.visible"] = False
1818
plt.rcParams["xtick.minor.visible"] = False
19+
20+
def get_set_e_scenarios():
21+
return [
22+
get_scenario("1342", name=1.94),
23+
get_scenario("M7", name=2),
24+
get_scenario("M10", name=2.5),
25+
get_scenario("M9", name=3),
26+
get_scenario("M6", name=4),
27+
get_scenario("M5", name=8),
28+
get_scenario("M11", name=12),
29+
get_scenario("M4", name=16),
30+
get_scenario("M14", name=18),
31+
get_scenario("M13", name=20),
32+
get_scenario("M8", name=24),
33+
get_scenario("M3", name=32),
34+
get_scenario("M12", name=48),
35+
get_scenario("M2", name=64),
36+
]

0 commit comments

Comments
 (0)