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.

../_images/yield_workflow.png

The yield assessment workflow, from [3]

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,
   ....: );
   ....: 
../_images/turbine_positions.png

To estimate the AEP, we will fist estimate wind climate at the hub height of each turbine:

In [18]: pwc = pw.wasp.generalize_and_downscale(output_locs, bwc, topo_map)

In [19]: print(pwc)
<xarray.Dataset>
Dimensions:        (sector: 12, point: 15)
Coordinates:
  * sector         (sector) float64 0.0 30.0 60.0 90.0 ... 270.0 300.0 330.0
    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_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.434 7.405 7.059 ... 7.77 8.327
    k              (sector, point) float32 2.041 2.045 2.053 ... 2.146 2.143
    wdfreq         (sector, point) float32 0.06445 0.06485 ... 0.07584 0.0792
    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.438 7.261 6.982 7.082 ... 7.555 7.825 8.11
    power_density  (point) float32 413.6 384.6 340.3 368.3 ... 425.6 481.4 531.7
Attributes:
    Conventions:      CF-1.8
    history:          2023-07-26T08:47:42:\twindkit==0.6.3\t wk.spatial.creat...
    title:            WAsP site effects
    Package name:     windkit
    Package version:  0.6.3
    Creation date:    2023-07-26T08:47:43
    Object type:      Met fields
    author:           Neil Davis
    author_email:     neda@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.6.3
    Creation date:    2023-07-26T08:47:43
    author:           Neil Davis
    author_email:     neda@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(pwc, 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.906 2.768 2.552 2.628 ... 3.0 3.194 3.412
Attributes:
    Conventions:      CF-1.8
    history:          2023-07-26T08:47:42:\twindkit==0.6.3\t wk.spatial.creat...
    title:            WAsP site effects
    Package name:     windkit
    Package version:  0.6.3
    Creation date:    2023-07-26T08:47:43
    Object type:      Anual Energy Production
    author:           Neil Davis
    author_email:     neda@dtu.dk
    institution:      DTU Wind

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 [24]: aep_wakes = pw.wasp.potential_aep(pwc, wtg, wind_farm_model="PARK2_onshore")

In [25]: print(aep_wakes)
<xarray.Dataset>
Dimensions:               (sector: 12, point: 15)
Coordinates:
  * sector                (sector) float64 0.0 30.0 60.0 ... 270.0 300.0 330.0
  * point                 (point) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
    height                (point) int64 50 50 50 50 50 50 ... 50 50 50 50 50 50
    south_north           (point) int64 4622313 4622199 ... 4624252 4624142
    west_east             (point) int64 513914 514161 514425 ... 516060 516295
    crs                   int8 0
Data variables:
    potential_AEP_sector  (point, sector) float64 0.1491 0.08913 ... 0.2294
Attributes:
    Conventions:      CF-1.8
    history:          2023-07-26T08:47:42:\twindkit==0.6.3\t wk.spatial.creat...
    title:            WAsP site effects
    Package name:     windkit
    Package version:  0.6.3
    Creation date:    2023-07-26T08:47:43
    Object type:      Anual Energy Production
    author:           Neil Davis
    author_email:     neda@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 [26]: 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 [27]: pwc_grid = pw.wasp.downscale(gwc, topo_map, output_grid, genwc_interp="nearest")

In [28]: flow_map = pw.wasp.wind_farm_flow_map(
   ....:     pwc_grid,
   ....:     wtg,
   ....:     output_locs,
   ....:     wind_farm_model="PARK2_onshore",
   ....: )
   ....: 

In [29]: 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
  * mode                             (mode) int64 0
    crs                              int8 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:          2023-07-26T08:47:43:\twindkit==0.6.3\twk.spatial.create...
    Package name:     windkit
    Package version:  0.6.3
    Creation date:    2023-07-26T08:48:07
    Object type:      Wind Farm Flow Map
    author:           Neil Davis
    author_email:     neda@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 [30]: flow_map["wspd_deficit_sector"].sel(sector=210).plot(vmax=0.1)
../_images/wspd_deficit.png

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 [31]: from functools import partial

In [32]: import py_wake

In [33]: wind_farm_model = partial(
   ....:     py_wake.wind_farm_models.engineering_models.All2AllIterative,
   ....:     wake_deficitModel=py_wake.deficit_models.noj.NOJDeficit(
   ....:         k=0.09, use_effective_ws=True
   ....:     ),
   ....:     rotorAvgModel=py_wake.rotor_avg_models.rotor_avg_model.RotorCenter(),
   ....:     superpositionModel=py_wake.superposition_models.LinearSum(),
   ....:     blockage_deficitModel=py_wake.deficit_models.Rathmann(),
   ....:     deflectionModel=None,
   ....:     turbulenceModel=py_wake.turbulence_models.GCLTurbulence(),
   ....: )
   ....: 

In [34]: aep_wakes = pw.wasp.potential_aep(pwc, wtg, wind_farm_model=wind_farm_model)

In [35]: print(aep_wakes)
<xarray.Dataset>
Dimensions:               (sector: 12, point: 15)
Coordinates:
  * sector                (sector) float64 0.0 30.0 60.0 ... 270.0 300.0 330.0
  * point                 (point) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
    height                (point) int64 50 50 50 50 50 50 ... 50 50 50 50 50 50
    south_north           (point) int64 4622313 4622199 ... 4624252 4624142
    west_east             (point) int64 513914 514161 514425 ... 516060 516295
    crs                   int8 0
Data variables:
    potential_AEP_sector  (point, sector) float64 0.149 0.08926 ... 0.234 0.2299
Attributes:
    Conventions:      CF-1.8
    history:          2023-07-26T08:47:42:\twindkit==0.6.3\t wk.spatial.creat...
    title:            WAsP site effects
    Package name:     windkit
    Package version:  0.6.3
    Creation date:    2023-07-26T08:48:08
    Object type:      Anual Energy Production
    author:           Neil Davis
    author_email:     neda@dtu.dk
    institution:      DTU Wind