Calculate AEP
PyWAsP can estimate the Annual Energy Production (AEP) of wind turbines.
For now, this includes estimating the gross and potential AEP (as per the workflow below), with only limited support, for now, for estimating other losses and P50 and P90.
We will start by importing needed libraries and some test data from the Serra Santa Luzia site.
In [1]: import numpy as np
In [2]: import pandas as pd
In [3]: import xarray as xr
In [4]: import matplotlib.pyplot as plt
In [5]: import windkit as wk
In [6]: import pywasp as pw
In [7]: bwc = wk.read_bwc("./source/tutorials/data/SerraSantaLuzia.omwc", "EPSG:4326")
In [8]: bwc = wk.spatial.reproject(bwc, "EPSG:32629")
In [9]: elev_map = wk.read_vector_map(
...: "./source/tutorials/data/SerraSantaLuzia.map",
...: crs="EPSG:32629",
...: map_type="elevation"
...: )
...:
In [10]: lc_map, lc_tbl = wk.read_vector_map(
....: "./source/tutorials/data/SerraSantaLuzia.map",
....: crs="EPSG:32629",
....: map_type="roughness"
....: )
....:
In [11]: topo_map = pw.wasp.TopographyMap(elev_map, lc_map, lc_tbl)
We have used the site in the WAsP Flow Model page, but here we will
go further and include some turbines. Their locations are stored in a .csv,
so we will create a xarray.Dataset to hold their spatial information: x, y, and hub height.
We can plot their locations together with the elevation map to get an idea of where they are located:
In [12]: wtg_locs = pd.read_csv('./source/tutorials/data/turbine_positions.csv')
In [13]: output_locs = wk.spatial.create_dataset(
....: wtg_locs.Easting.values,
....: wtg_locs.Northing.values,
....: wtg_locs['Hub height'].values,
....: crs="EPSG:32629"
....: )
....:
In [14]: fig, ax = plt.subplots()
In [15]: ax.set(xlim=(510000, 521500), ylim=(4616000, 4627500))
Out[15]: [(510000.0, 521500.0), (4616000.0, 4627500.0)]
In [16]: elev_map.plot("elev", ax=ax);
In [17]: ax.scatter(
....: output_locs.west_east.values,
....: output_locs.south_north.values,
....: c="k",
....: marker="x",
....: s=15,
....: zorder=2,
....: );
....:
To estimate the AEP, we will fist estimate wind climate at the hub height of each turbine:
In [18]: wwc = pw.wasp.predict_wwc(bwc, topo_map, output_locs)
In [19]: print(wwc)
<xarray.Dataset>
Dimensions: (point: 15, sector: 12)
Coordinates:
height (point) int64 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
south_north (point) int64 4622313 4622199 4622336 ... 4624252 4624142
west_east (point) int64 513914 514161 514425 ... 515808 516060 516295
crs int8 0
* sector (sector) float64 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
sector_ceil (sector) float64 15.0 45.0 75.0 105.0 ... 285.0 315.0 345.0
sector_floor (sector) float64 345.0 15.0 45.0 75.0 ... 255.0 285.0 315.0
Dimensions without coordinates: point
Data variables:
A (sector, point) float32 7.432 7.405 7.063 ... 7.716 8.268
k (sector, point) float32 2.025 2.033 2.041 ... 2.131 2.127
wdfreq (sector, point) float32 0.06444 0.06484 ... 0.07591 0.07929
site_elev (point) float32 459.7 460.0 460.0 427.1 ... 520.0 520.0 520.0
air_density (point) float32 1.163 1.163 1.163 1.167 ... 1.156 1.156 1.156
wspd (point) float32 7.485 7.307 7.027 7.126 ... 7.595 7.861 8.145
power_density (point) float32 422.7 392.7 347.5 375.1 ... 433.9 489.4 540.2
Attributes:
Conventions: CF-1.8
history: 2024-06-11T13:35:50+00:00:\twindkit==0.8.0\t wk.spatial...
title: WAsP site effects
Package name: windkit
Package version: 0.8.0
Creation date: 2024-06-11T13:35:53+00:00
Object type: Met fields
author: DTU CI Config
author_email: pywasp@dtu.dk
institution: DTU Wind
Gross AEP
Before we estiamte the AEP we have to get a Wind Turbine Generator (WTG) corresponding to the turbine models used.
We use windkit.read_wtg() to read the .wtg file.
In [20]: wtg = wk.read_wtg("./source/tutorials/data/Bonus_1_MW.wtg")
In [21]: print(wtg)
<xarray.Dataset>
Dimensions: (mode: 1, wind_speed: 22)
Coordinates:
* wind_speed (wind_speed) float64 4.0 5.0 ... 24.0 25.0
* mode (mode) int64 0
Data variables:
power_output (mode, wind_speed) float64 2.41e+04 ... 1e+06
thrust_coefficient (mode, wind_speed) float64 0.915 ... 0.161
air_density (mode) float64 1.225
stationary_thrust_coefficient (mode) float64 0.161
wind_speed_cutin (mode) float64 4.0
wind_speed_cutout (mode) float64 25.0
rated_power (mode) float64 1e+06
name <U10 'Bonus 1 MW'
rotor_diameter float64 54.2
hub_height float64 50.0
regulation_type int64 2
Attributes:
Conventions: CF-1.8
Package name: windkit
Package version: 0.8.0
Creation date: 2024-06-11T13:35:53+00:00
author: DTU CI Config
author_email: pywasp@dtu.dk
institution: DTU Wind
The wtg holds power and thrust curves for a number of wind speeds and operational modes (in this case just 1).
Now we can estimate the gross AEP using gross_aep():
In [22]: aep = pw.wasp.gross_aep(wwc, wtg)
In [23]: print(aep)
<xarray.Dataset>
Dimensions: (point: 15)
Coordinates:
height (point) int64 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
south_north (point) int64 4622313 4622199 4622336 ... 4624252 4624142
west_east (point) int64 513914 514161 514425 ... 515808 516060 516295
crs int8 0
Dimensions without coordinates: point
Data variables:
gross_AEP (point) float32 2.945 2.807 2.591 2.665 ... 3.032 3.222 3.439
Attributes:
Conventions: CF-1.8
history: 2024-06-11T13:35:50+00:00:\twindkit==0.8.0\t wk.spatial...
title: WAsP site effects
Package name: windkit
Package version: 0.8.0
Creation date: 2024-06-11T13:35:53+00:00
Object type: Anual Energy Production
author: DTU CI Config
author_email: pywasp@dtu.dk
institution: DTU Wind
Passing simply a wwc and a WTG works well when your climate is predicted for exactly the hub heights and locations of the turbines
you want to know the AEP for, but sometimes you have different turbines associated with different WTGs.
To handle this, you can pass a dictionary of WTG’s and a dataset of turbine locations, heights, turbine id’s, group id’s, and wtg keys that
map to the WTG’s in the dictionary. Let’s try it here:
In [24]: wtg_dict = {"Bonus_1_MW": wtg}
In [25]: wind_turbines = wk.create_wind_turbines_from_arrays(
....: west_east=output_locs.west_east.values,
....: south_north=output_locs.south_north.values,
....: height=output_locs.height.values,
....: wtg_keys=["Bonus_1_MW"]*15,
....: group_ids=np.zeros(15, dtype=int),
....: turbine_ids=np.arange(15, dtype=int),
....: )
....:
In [26]: aep = pw.wasp.gross_aep(wwc, wtg_dict, wind_turbines)
In [27]: print(aep)
<xarray.Dataset>
Dimensions: (point: 15)
Coordinates:
height (point) int64 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
south_north (point) int64 4622313 4622199 4622336 ... 4624252 4624142
west_east (point) int64 513914 514161 514425 ... 515808 516060 516295
crs int8 0
Dimensions without coordinates: point
Data variables:
gross_AEP (point) float32 2.945 2.807 2.591 2.665 ... 3.032 3.222 3.439
Attributes:
Conventions: CF-1.8
history: 2024-06-11T13:35:50+00:00:\twindkit==0.8.0\t wk.spatial...
title: WAsP site effects
Package name: windkit
Package version: 0.8.0
Creation date: 2024-06-11T13:35:53+00:00
Object type: Anual Energy Production
author: DTU CI Config
author_email: pywasp@dtu.dk
institution: DTU Wind
Here, a single type of turbine is used, but it is possible to add more turbines to the collection by passing several WTG’s in the dictionary and adding more turbine locations to the dataset. The turbine id’s and group id’s are used to group turbines together for calculating the AEP. The turbine id’s are used to identify individual turbines,
Potential AEP
To add wind farm effects (such as wakes and blockage) to the AEP estimation and obtain the
“potential AEP”, PyWAsP integrates with PyWake. PyWAsP defines a number
of named wind farm models (see potential_aep()), including “PARK2_onshore” and “PARK2_offshore”,
that are indentical to the wake model in the WAsP GUI. Here we will use the onshore version:
In [28]: aep_wakes = pw.wasp.potential_aep(wwc, wtg_dict, wind_turbines, wind_farm_model="PARK2_onshore", turbulence_intensity=0.0)
In [29]: print(aep_wakes)
<xarray.Dataset>
Dimensions: (sector: 12, point: 15, mode: 1)
Coordinates:
* sector (sector) float64 0.0 30.0 ... 300.0 330.0
height (point) int64 50 50 50 50 ... 50 50 50 50
south_north (point) int64 4622313 4622199 ... 4624142
west_east (point) int64 513914 514161 ... 516295
crs int8 0
* point (point) int64 0 1 2 3 4 ... 10 11 12 13 14
* mode (mode) int64 0
Data variables: (12/15)
potential_AEP_sector (point, sector) float64 0.1492 ... 0.2267
gross_AEP_sector (mode, point, sector) float64 0.15 ... 0...
AEP_deficit_sector (point, mode, sector) float64 0.004901 ....
wspd_sector (point, sector) float64 6.585 ... 7.323
wspd_eff_sector (point, sector) float64 6.585 ... 7.323
wspd_deficit_sector (point, sector) float64 0.0 0.00249 ... 0.0
... ...
gross_AEP (mode, point) float64 2.955 2.818 ... 3.446
AEP_deficit (point, mode) float64 0.0462 ... 0.03699
wspd (point) float64 7.485 7.307 ... 7.861 8.145
wspd_eff (point) float64 7.359 7.19 ... 7.653 8.0
wspd_deficit (point) float64 0.01847 0.01814 ... 0.01782
turbulence_intensity_eff (point) float64 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
Conventions: CF-1.8
history: 2024-06-11T13:35:50+00:00:\twindkit==0.8.0\t wk.spatial...
title: WAsP site effects
Package name: windkit
Package version: 0.8.0
Creation date: 2024-06-11T13:35:54+00:00
Object type: Anual Energy Production
author: DTU CI Config
author_email: pywasp@dtu.dk
institution: DTU Wind
Wind farm flow map
To map an area around a wind farm PyWAsP has wind_farm_flow_map(), which
takes a resource grid as input and enriches it with wind farm variables including
effects of nearby turbines:
In [30]: output_grid = wk.spatial.create_dataset(
....: np.linspace(511700.0, 517700.0, 41),
....: np.linspace(4619500.0, 4625500.0, 41),
....: [50.0],
....: crs="EPSG:32629",
....: struct="cuboid",
....: )
....:
In [31]: wwc_grid = pw.wasp.downscale(gwc, topo_map, output_grid, interp_method="nearest")
In [32]: flow_map = pw.wasp.wind_farm_flow_map(
....: wwc_grid,
....: wtg_dict,
....: wind_turbines,
....: output_grid,
....: wind_farm_model="PARK2_onshore",
....: turbulence_intensity=0.0
....: )
....:
In [33]: print(flow_map)
<xarray.Dataset>
Dimensions: (sector: 12, height: 1, south_north: 41,
west_east: 41, mode: 1)
Coordinates:
* sector (sector) float64 0.0 30.0 ... 300.0 330.0
* height (height) float64 50.0
* south_north (south_north) float64 4.62e+06 ... 4.626...
* west_east (west_east) float64 5.117e+05 ... 5.177e+05
crs int8 0
* mode (mode) int64 0
Data variables: (12/15)
potential_AEP_sector (height, south_north, west_east, sector) float64 ...
gross_AEP_sector (height, south_north, west_east, sector, mode) float64 ...
AEP_deficit_sector (height, south_north, west_east, sector, mode) float64 ...
wspd_sector (height, south_north, west_east, sector) float64 ...
wspd_eff_sector (height, south_north, west_east, sector) float64 ...
wspd_deficit_sector (height, south_north, west_east, sector) float64 ...
... ...
gross_AEP (height, south_north, west_east, mode) float64 ...
AEP_deficit (height, south_north, west_east, mode) float64 ...
wspd (height, south_north, west_east) float64 ...
wspd_eff (height, south_north, west_east) float64 ...
wspd_deficit (height, south_north, west_east) float64 ...
turbulence_intensity_eff (height, south_north, west_east) float64 ...
Attributes:
Conventions: CF-1.8
history: 2024-06-11T13:35:54+00:00:\twindkit==0.8.0\twk.spatial....
Package name: windkit
Package version: 0.8.0
Creation date: 2024-06-11T13:36:35+00:00
Object type: Wind Farm Flow Map
author: DTU CI Config
author_email: pywasp@dtu.dk
institution: DTU Wind
title: WAsP site effects
We can plot the wind speed deficit for one direction to show the wake effects:
In [34]: flow_map["wspd_deficit_sector"].sel(sector=210).plot(vmax=0.1)
PyWake integration in PyWAsP
PyWAsP can use PyWake wind farm models to estimate the potential AEP. This is done by passing the name of the PyWake wind farm model to the potential_aep() and wind_farm_flow_map() functions.
To ensure stability and consistency of PyWAsP, the PyWake versions for each PyWAsP release is pinned to a specific version of PyWake. When installing PyWAsP via mamba, the correct version of PyWake is installed automatically.
We periodically update the pinned version of PyWake to the latest version that is compatible with PyWAsP.
Using custom PyWake wind farm models
Instead of using a named wind farm model, it is possible to construct a
PyWake wind farm model manually, using functools.partial(), and use it for adding wind farm effects in PyWAsP:
In [35]: from functools import partial
In [36]: import py_wake
In [37]: wind_farm_model = partial(
....: py_wake.wind_farm_models.engineering_models.All2AllIterative,
....: wake_deficitModel=py_wake.deficit_models.noj.NOJLocalDeficit(
....: a=[0, 0.06],
....: ct2a=py_wake.deficit_models.utils.ct2a_mom1d,
....: use_effective_ws=True,
....: use_effective_ti=False,
....: rotorAvgModel=py_wake.rotor_avg_models.AreaOverlapAvgModel(),
....: ),
....: superpositionModel=py_wake.superposition_models.LinearSum(),
....: blockage_deficitModel=py_wake.deficit_models.Rathmann(),
....: deflectionModel=None,
....: turbulenceModel=None,
....: )
....:
In [38]: aep_wakes = pw.wasp.potential_aep(wwc, wtg_dict, wind_turbines, wind_farm_model=wind_farm_model, turbulence_intensity=0.0)
In [39]: print(aep_wakes)
<xarray.Dataset>
Dimensions: (sector: 12, point: 15, mode: 1)
Coordinates:
* sector (sector) float64 0.0 30.0 ... 300.0 330.0
height (point) int64 50 50 50 50 ... 50 50 50 50
south_north (point) int64 4622313 4622199 ... 4624142
west_east (point) int64 513914 514161 ... 516295
crs int8 0
* point (point) int64 0 1 2 3 4 ... 10 11 12 13 14
* mode (mode) int64 0
Data variables: (12/15)
potential_AEP_sector (point, sector) float64 0.1492 ... 0.2273
gross_AEP_sector (mode, point, sector) float64 0.15 ... 0...
AEP_deficit_sector (point, mode, sector) float64 0.005137 ....
wspd_sector (point, sector) float64 6.585 ... 7.323
wspd_eff_sector (point, sector) float64 6.584 ... 7.331
wspd_deficit_sector (point, sector) float64 9.693e-05 ... -0...
... ...
gross_AEP (mode, point) float64 2.955 2.818 ... 3.446
AEP_deficit (point, mode) float64 0.05318 ... 0.04445
wspd (point) float64 7.485 7.307 ... 7.861 8.145
wspd_eff (point) float64 7.337 7.161 ... 7.617 7.967
wspd_deficit (point) float64 0.0217 0.0225 ... 0.0218
turbulence_intensity_eff (point) float64 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
Conventions: CF-1.8
history: 2024-06-11T13:35:50+00:00:\twindkit==0.8.0\t wk.spatial...
title: WAsP site effects
Package name: windkit
Package version: 0.8.0
Creation date: 2024-06-11T13:36:35+00:00
Object type: Anual Energy Production
author: DTU CI Config
author_email: pywasp@dtu.dk
institution: DTU Wind
Calculating AEP including wake effects
In PyWAsP, the “park” method is used in py_wake.Site.XRSsite.from_pwc() for calculating the speed-up when passing the wind climate to a py_wake site object.
The steps used to calculate the AEP including wake effects are as follows:
Flow cases are defined to ensure all possible/relevant directions and speeds are covered. The wind speed in each flow case is relative to the maximum mean wind speed in the sector.
The wind farm simulation results from py_wake contains the variable
WS, which is inhomogeneous background wind speed used to obtain the probability density from the indivudal sector-wise weibull distributions. The variablesWS_effcontains the effective wind speed including wind farm effects.The AEP is determined by finding the effective power curves relative to the background wind speed
P_eff(WS)and performing an integral over wind speed.