Wind farm optimisation#
import os
import cartopy.crs as ccrs
import contextily as cx
import geopandas as gpd
import mapclassify as mc
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
from matplotlib import patheffects
from matplotlib_scalebar.scalebar import ScaleBar
from h2ss import capacity as cap
from h2ss import compare
from h2ss import data as rd
from h2ss import functions as fns
from h2ss import optimisation as opt
# basemap cache directory
cx.set_cache_dir(os.path.join("data", "basemaps"))
Halite data#
ds, extent = rd.kish_basin_data_depth_adjusted(
dat_path=os.path.join("data", "kish-basin"),
bathymetry_path=os.path.join("data", "bathymetry"),
)
Constraints#
# exploration wells
_, wells_b = fns.constraint_exploration_well(
data_path=os.path.join(
"data",
"exploration-wells",
"Exploration_Wells_Irish_Offshore.shapezip.zip",
)
)
# wind farms
wind_farms = fns.constraint_wind_farm(
data_path=os.path.join(
"data", "wind-farms", "marine-area-consent-wind.zip"
)
)
# frequent shipping routes
_, shipping_b = fns.constraint_shipping_routes(
data_path=os.path.join(
"data", "shipping", "shipping_frequently_used_routes.zip"
),
dat_extent=extent,
)
# shipwrecks
_, shipwrecks_b = fns.constraint_shipwrecks(
data_path=os.path.join(
"data", "shipwrecks", "IE_GSI_MI_Shipwrecks_IE_Waters_WGS84_LAT.zip"
),
dat_extent=extent,
)
# subsea cables
_, cables_b = fns.constraint_subsea_cables(
data_path=os.path.join("data", "subsea-cables", "KIS-ORCA.gpkg"),
dat_extent=extent,
)
# distance from salt formation edge
edge_buffer = fns.constraint_halite_edge(dat_xr=ds)
Zones of interest#
zones, zds = fns.zones_of_interest(
dat_xr=ds,
constraints={"net_height": 120, "min_depth": 500, "max_depth": 2000},
)
Generate caverns#
caverns = fns.generate_caverns_hexagonal_grid(
zones_df=zones,
dat_extent=extent,
)
caverns = fns.cavern_dataframe(
dat_zone=zds,
cavern_df=caverns,
depths={"min": 500, "min_opt": 1000, "max_opt": 1500, "max": 2000},
)
# label caverns by depth and heights
caverns = fns.label_caverns(
cavern_df=caverns,
heights=[120],
depths={"min": 500, "min_opt": 1000, "max_opt": 1500, "max": 2000},
)
caverns, _ = fns.generate_caverns_with_constraints(
cavern_df=caverns,
exclusions={
"wells": wells_b,
"wind_farms": wind_farms,
"shipwrecks": shipwrecks_b,
"shipping": shipping_b,
"cables": cables_b,
"edge": edge_buffer,
},
)
Without constraints...
Number of potential caverns: 568
------------------------------------------------------------
Excluding salt formation edges...
Number of potential caverns: 539
------------------------------------------------------------
Exclude shipping...
Number of potential caverns: 261
Caverns excluded: 51.58%
------------------------------------------------------------
Exclude wind farms...
Number of potential caverns: 218
Caverns excluded: 59.55%
------------------------------------------------------------
Exclude cables...
Number of potential caverns: 218
Caverns excluded: 59.55%
------------------------------------------------------------
Exclude wells...
Number of potential caverns: 218
Caverns excluded: 59.55%
------------------------------------------------------------
Exclude shipwrecks...
Number of potential caverns: 218
Caverns excluded: 59.55%
------------------------------------------------------------
Capacity#
caverns["cavern_total_volume"] = cap.cavern_volume(
height=caverns["cavern_height"]
)
caverns["cavern_volume"] = cap.corrected_cavern_volume(
v_cavern=caverns["cavern_total_volume"]
)
caverns["t_mid_point"] = cap.temperature_cavern_mid_point(
height=caverns["cavern_height"], depth_top=caverns["cavern_depth"]
)
(
caverns["p_operating_min"],
caverns["p_operating_max"],
) = cap.pressure_operating(
thickness_overburden=caverns["TopDepthSeabed"],
depth_water=-caverns["Bathymetry"],
)
caverns["rho_min"], caverns["rho_max"] = cap.density_hydrogen_gas(
p_operating_min=caverns["p_operating_min"],
p_operating_max=caverns["p_operating_max"],
t_mid_point=caverns["t_mid_point"],
)
(
caverns["working_mass"],
caverns["mass_operating_min"],
caverns["mass_operating_max"],
) = cap.mass_hydrogen_working(
rho_h2_min=caverns["rho_min"],
rho_h2_max=caverns["rho_max"],
v_cavern=caverns["cavern_volume"],
)
caverns["capacity"] = cap.energy_storage_capacity(
m_working=caverns["working_mass"]
)
Power curve [MW] and Weibull wind speed distribution#
# extract data for wind farms at 150 m
weibull_wf_df = fns.read_weibull_data(
data_path_weibull=os.path.join(
"data", "weibull-parameters-wind-speeds", "Weibull_150m_params_ITM.zip"
),
data_path_wind_farms=os.path.join(
"data", "wind-farms", "marine-area-consent-wind.zip"
),
)
weibull_powercurve = opt.weibull_distribution(weibull_wf_data=weibull_wf_df)
weibull_powercurve.head()
| wind_speed | power_curve | Codling Wind Park | Dublin Array | North Irish Sea Array | |
|---|---|---|---|---|---|
| 0 | 0.00 | 0.0 | 0.000000 | 0.000000 | 0.000000 |
| 1 | 0.01 | 0.0 | 0.000250 | 0.000260 | 0.000070 |
| 2 | 0.02 | 0.0 | 0.000484 | 0.000503 | 0.000154 |
| 3 | 0.03 | 0.0 | 0.000711 | 0.000739 | 0.000243 |
| 4 | 0.04 | 0.0 | 0.000935 | 0.000972 | 0.000337 |
ax = weibull_powercurve.plot(
x="wind_speed",
y="power_curve",
linewidth=3,
color="tab:blue",
figsize=(10, 5.5),
legend=False,
)
ax.set_xlabel("Wind speed [m s\N{SUPERSCRIPT MINUS}\N{SUPERSCRIPT ONE}]")
ax.set_ylabel("Power [MW]")
plt.yticks([3 * n for n in range(6)])
sns.despine()
ax.xaxis.grid(True, linewidth=0.25)
ax.yaxis.grid(True, linewidth=0.25)
plt.tight_layout()
# plt.savefig(
# os.path.join("graphics", "Figure_6.png"), format="png", dpi=600
# )
plt.show()
plt.figure(figsize=(10, 5.5))
ax = sns.lineplot(
data=weibull_powercurve.drop(columns=["power_curve"]).melt(
id_vars="wind_speed"
),
x="wind_speed",
y="value",
hue="variable",
linestyle="dashed",
linewidth=2.25,
alpha=0.85,
palette=sns.color_palette(["tab:red", "tab:gray", "tab:blue"]),
)
ax.set_xlabel("Wind speed [m s\N{SUPERSCRIPT MINUS}\N{SUPERSCRIPT ONE}]")
ax.set_ylabel(
"Weibull probability distribution function "
"[s m\N{SUPERSCRIPT MINUS}\N{SUPERSCRIPT ONE}]"
)
ax.xaxis.grid(True, linewidth=0.25)
ax.yaxis.grid(True, linewidth=0.25)
sns.despine()
ax.legend(title=None, fontsize=10.5)
plt.tight_layout()
# plt.savefig(
# os.path.join("graphics", "fig_weibull.jpg"), format="jpg", dpi=600
# )
plt.show()
Number of reference wind turbines#
# number of 15 MW turbines, rounded down to the nearest integer
weibull_wf_df["n_turbines"] = opt.number_of_turbines(
owf_cap=weibull_wf_df["cap"]
)
Annual energy production [MWh]#
weibull_wf_df = opt.annual_energy_production(weibull_wf_data=weibull_wf_df)
Annual hydrogen production [kg]#
weibull_wf_df["AHP"] = opt.annual_hydrogen_production(aep=weibull_wf_df["AEP"])
AHP as a proportion of the total working mass#
weibull_wf_df["AHP_frac"] = (
weibull_wf_df["AHP"] / caverns[["working_mass"]].sum().iloc[0]
)
AHP converted to storage demand [GWh]#
weibull_wf_df["demand"] = cap.energy_storage_capacity(
m_working=weibull_wf_df["AHP"]
)
Number of caverns required based on cumulative working mass and AHP#
compare.calculate_number_of_caverns(
cavern_df=caverns, weibull_wf_data=weibull_wf_df
)
Codling Wind Park
Working mass [kg]: 1.095672E+08
Number of caverns required: 25–51
Capacity (approx.) [GWh]: 3,735.76
------------------------------------------------------------------------------
Dublin Array
Working mass [kg]: 6.742514E+07
Number of caverns required: 15–33
Capacity (approx.) [GWh]: 2,289.53
------------------------------------------------------------------------------
North Irish Sea Array
Working mass [kg]: 4.515738E+07
Number of caverns required: 10–24
Capacity (approx.) [GWh]: 1,556.50
------------------------------------------------------------------------------
Total number of caverns required: 50–108
------------------------------------------------------------------------------
Number of caverns required as a percentage of all available caverns:
22.94–49.54%
------------------------------------------------------------------------------
Total maximum cavern capacity (approx.): 7,581.79 GWh
Transmission distance [km]#
caverns, injection_point = opt.transmission_distance(
cavern_df=caverns, wf_data=wind_farms
)
Electrolyser capacity [MW]#
weibull_wf_df["E_cap"] = opt.electrolyser_capacity(
n_turbines=weibull_wf_df["n_turbines"]
)
CAPEX for pipeline [€ km⁻¹]#
weibull_wf_df["CAPEX"] = opt.capex_pipeline(e_cap=weibull_wf_df["E_cap"])
weibull_wf_df
| name | cap | (c, min) | (c, max) | (c, mean) | (k, min) | (k, max) | (k, mean) | n_turbines | AEP | integral | abserr | AHP | AHP_frac | demand | E_cap | CAPEX | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Codling Wind Park | 1300 | 10.2 | 10.8 | 10.500000 | 1.9 | 2.0 | 1.950000 | 86 | 5.550255e+06 | 8.335975 | 2.528685e-07 | 1.095672e+08 | 0.154167 | 3651.022453 | 1070 | 1.451786e+06 |
| 1 | Dublin Array | 824 | 9.9 | 10.6 | 10.292857 | 1.9 | 2.0 | 1.950000 | 54 | 3.415500e+06 | 8.169630 | 3.778893e-07 | 6.742514e+07 | 0.094871 | 2246.755615 | 672 | 1.205538e+06 |
| 2 | North Irish Sea Array | 500 | 10.7 | 11.2 | 10.950000 | 2.1 | 2.2 | 2.133333 | 33 | 2.287500e+06 | 8.953423 | 1.516351e-07 | 4.515738e+07 | 0.063539 | 1504.744246 | 410 | 1.032447e+06 |
# totals
weibull_wf_df[
["cap", "n_turbines", "AEP", "AHP", "AHP_frac", "demand", "E_cap", "CAPEX"]
].sum()
cap 2.624000e+03
n_turbines 1.730000e+02
AEP 1.125325e+07
AHP 2.221497e+08
AHP_frac 3.125772e-01
demand 7.402522e+03
E_cap 2.152000e+03
CAPEX 3.689771e+06
dtype: float64
compare.electricity_demand_ie(data=weibull_wf_df["demand"])
Energy capacity as a percentage of Ireland's electricity demand
in 2050 (84–122 TWh electricity), assuming a conversion efficiency
of 50%: 3.03–4.41%
Energy capacity as a percentage of Ireland's hydrogen demand
in 2050, assuming it is 17% of the electricity demand
(84–122 TWh electricity): 35.69–51.84%
compare.hydrogen_demand_ie(data=weibull_wf_df["demand"])
Energy capacity as a percentage of Ireland's domestic hydrogen
demand in 2050 (4.6–39 TWh hydrogen): 18.98–160.92%
Energy capacity as a percentage of Ireland's domestic and
non-domestic hydrogen demand in 2050 (19.8–74.6 TWh hydrogen): 9.92–37.39%
LCOT for pipeline [€ kg⁻¹]#
caverns = opt.lcot_pipeline(weibull_wf_data=weibull_wf_df, cavern_df=caverns)
caverns[
[
"cavern_depth",
"working_mass",
"capacity",
"distance_ip",
]
+ list(caverns.filter(like="dist_"))
+ list(caverns.filter(like="LCOT_"))
].describe()
| cavern_depth | working_mass | capacity | distance_ip | dist_Codling_Wind_Park | dist_Dublin_Array | dist_North_Irish_Sea_Array | LCOT_Codling_Wind_Park | LCOT_Dublin_Array | LCOT_North_Irish_Sea_Array | |
|---|---|---|---|---|---|---|---|---|---|---|
| count | 218.000000 | 2.180000e+02 | 218.000000 | 218.000000 | 218.000000 | 218.000000 | 218.000000 | 218.000000 | 218.000000 | 218.000000 |
| mean | 1154.020877 | 3.260108e+06 | 108.634041 | 29.197826 | 58.408245 | 43.581633 | 44.393167 | 0.078615 | 0.079154 | 0.103102 |
| std | 365.142804 | 7.786651e+05 | 25.946851 | 6.830597 | 11.452952 | 13.055932 | 7.523248 | 0.015415 | 0.023713 | 0.017473 |
| min | 502.745339 | 1.678950e+06 | 55.946359 | 16.263066 | 33.027322 | 17.542931 | 32.980512 | 0.044454 | 0.031862 | 0.076596 |
| 25% | 872.651330 | 2.685380e+06 | 89.482842 | 23.216441 | 51.592318 | 33.503367 | 36.808911 | 0.069441 | 0.060850 | 0.085488 |
| 50% | 1128.329146 | 3.277441e+06 | 109.211632 | 30.947028 | 58.988376 | 47.378240 | 44.619590 | 0.079396 | 0.086050 | 0.103628 |
| 75% | 1433.385847 | 3.886101e+06 | 129.493525 | 34.863247 | 67.109072 | 54.107289 | 50.773606 | 0.090326 | 0.098271 | 0.117920 |
| max | 1980.680682 | 4.765555e+06 | 158.798899 | 42.990158 | 80.806729 | 70.176200 | 67.926218 | 0.108763 | 0.127456 | 0.157757 |
pd.Series(caverns[list(caverns.filter(like="dist_"))].values.flat).describe()
count 654.000000
mean 48.794348
std 12.862247
min 17.542931
25% 36.924385
50% 49.890968
75% 57.488679
max 80.806729
dtype: float64
pd.Series(caverns[list(caverns.filter(like="LCOT_"))].values.flat).describe()
count 654.000000
mean 0.086957
std 0.022312
min 0.031862
25% 0.071538
50% 0.087836
75% 0.101471
max 0.157757
dtype: float64
fig, axes = plt.subplots(1, 2, figsize=(10.5, 4.5))
for n, col, lab, show_legend in zip(
[0, 1],
["dist_", "LCOT_"],
["Transmission distance [km]", "Pipeline LCOT [€ kg⁻¹]"],
[False, True],
):
sns.boxplot(
caverns.filter(like=col)
.set_axis(list(wind_farms["name"]), axis=1)
.melt(),
y="value",
hue="variable",
palette=sns.color_palette(["tab:red", "tab:gray", "tab:blue"]),
width=0.35,
ax=axes[n],
legend=show_legend,
linecolor="black",
linewidth=1.1,
gap=0.15,
showmeans=True,
meanprops={
"marker": "d",
"markeredgecolor": "black",
"markerfacecolor": "none",
},
)
axes[n].set_ylabel(lab)
axes[n].tick_params(axis="x", bottom=False)
axes[n].yaxis.grid(True, linewidth=0.25)
axes[1].legend(loc="lower right", fontsize=11.25)
sns.despine(bottom=True)
plt.tight_layout()
# plt.savefig(
# os.path.join("graphics", "Figure_10.png"),
# format="png",
# dpi=600,
# )
plt.show()
Maps#
shape = rd.halite_shape(dat_xr=ds).buffer(1000).buffer(-1000)
def plot_map_facet(
cavern_df,
classes,
fontsize=11.5,
):
"""Helper function for plotting LCOT facet maps"""
fig1 = plt.figure(figsize=(11, 11.75))
xmin_, ymin_, xmax_, ymax_ = cavern_df.total_bounds
colours = [int(n * 255 / (len(classes) - 1)) for n in range(len(classes))]
legend_handles = []
classes = sorted(classes)
for n1, c in enumerate(colours):
if n1 == 0:
label = f"< {classes[n1]:.2f}"
elif n1 == len(colours) - 1:
label = f"≥ {classes[-2]:.2f}"
else:
label = f"{classes[n1 - 1]:.2f} – < {classes[n1]:.2f}"
legend_handles.append(
mpatches.Patch(
facecolor=sns.color_palette("flare", 256)[c], label=label
)
)
for a, wf1 in enumerate(list(wind_farms["name"])):
ax1 = fig1.add_subplot(2, 2, a + 1, projection=ccrs.epsg(rd.CRS))
gpd.GeoDataFrame(cavern_df, geometry=cavern_df.centroid).plot(
ax=ax1,
scheme="UserDefined",
classification_kwds={"bins": classes},
column=f"LCOT_{wf1.replace(' ', '_')}",
zorder=2,
marker=".",
cmap="flare",
markersize=20,
)
shape.plot(
ax=ax1, color="white", alpha=0.5, edgecolor="slategrey", zorder=1
)
cx.add_basemap(
ax1,
crs=rd.CRS,
source=cx.providers.CartoDB.VoyagerNoLabels,
attribution=False,
)
ax1.gridlines(
draw_labels={"bottom": "x"},
color="lightslategrey",
alpha=0.25,
xlabel_style={"fontsize": fontsize},
xformatter=LongitudeFormatter(auto_hide=False, dms=True),
)
if not a == 1:
ax1.gridlines(
draw_labels={"left": "y"},
color="lightslategrey",
alpha=0.25,
ylabel_style={"fontsize": fontsize, "rotation": 89.9},
yformatter=LatitudeFormatter(auto_hide=False, dms=True),
)
if a == 2:
ax1.add_artist(
ScaleBar(
1,
box_alpha=0,
location="lower right",
color="darkslategrey",
font_properties={"size": fontsize},
)
)
plt.legend(
loc="lower right",
bbox_to_anchor=(1, 0.075),
handles=legend_handles,
title="Pipeline LCOT [€ kg⁻¹]",
fontsize=fontsize,
title_fontsize=fontsize,
)
plt.xlim(xmin_ - 1000, xmax_ + 1000)
plt.ylim(ymin_ - 1000, ymax_ + 1000)
ax1.set_title(list(wind_farms["name"])[a])
plt.tight_layout()
# plt.savefig(
# os.path.join("graphics", "Figure_11.png"),
# format="png",
# dpi=600,
# )
plt.show()
plot_map_facet(
caverns,
list(mc.Quantiles(caverns[list(caverns.filter(like="LCOT_"))], k=6).bins),
)
def plot_map_extent(cavern_df):
"""Helper function for plotting extent map"""
xmin_, _, xmax_, _ = extent.total_bounds
_, ymin_, _, ymax_ = wind_farms.total_bounds
plt.figure(figsize=(6, 6))
ax2 = plt.axes(projection=ccrs.epsg(rd.CRS))
shape.plot(
ax=ax2, color="white", alpha=0.5, edgecolor="slategrey", zorder=1
)
gpd.GeoDataFrame(cavern_df, geometry=cavern_df.centroid).plot(
ax=ax2, marker=".", markersize=2.5, color="firebrick"
)
gpd.GeoDataFrame(cavern_df, geometry=cavern_df.centroid).dissolve().buffer(
1000
).envelope.boundary.plot(ax=ax2, color="darkslategrey")
wind_farms.plot(
ax=ax2, facecolor="none", hatch="///", edgecolor="royalblue"
)
plt.xlim(xmin_ - 25000, xmax_ + 9500)
plt.ylim(ymin_ - 3000, ymax_ + 3000)
injection_point.plot(ax=ax2, marker="*", color="darkslategrey")
basemap = cx.providers.CartoDB.VoyagerNoLabels
cx.add_basemap(
ax2,
crs=rd.CRS,
source=basemap,
attribution=False,
)
ax2.text(xmin_ - 24250, ymin_ - 2150, basemap["attribution"], fontsize=9)
map_labels = zip(
zip(wind_farms.centroid.x, wind_farms.centroid.y), wind_farms["name"]
)
for xy, lab2 in map_labels:
ax2.annotate(
text=lab2,
xy=xy,
path_effects=[
patheffects.withStroke(linewidth=2.5, foreground="white")
],
fontsize=13,
va="center",
fontweight="semibold",
)
ax2.annotate(
text="DUBLIN\nPORT",
xy=(
injection_point.centroid.x[0] - 1500,
injection_point.centroid.y[0],
),
path_effects=[patheffects.withStroke(linewidth=2, foreground="white")],
fontsize=13,
va="center",
ha="right",
font="Montserrat",
fontweight="semibold",
)
ax2.add_artist(
ScaleBar(
1,
box_alpha=0,
location="lower right",
color="darkslategrey",
font_properties={"size": 11.5},
)
)
plt.tight_layout()
# plt.savefig(
# os.path.join("graphics", "Figure_11_inset.png"),
# format="png",
# dpi=600,
# )
plt.show()
plot_map_extent(caverns)
