.. py:currentmodule:: pywasp.wasp

.. _calculate_aep:

=============
Calculate AEP
=============

PyWAsP can estimate the Annual Energy Production (AEP) of wind turbines.

This includes estimating the gross and potential AEP (as shown in the workflow below),
as well as support for estimating net AEP (P50) and P90.

.. figure:: ../_static/yield_workflow.png
    :align: center

    The yield assessment workflow, from :cite:`mortensen2015offshore`

We will start by importing the necessary libraries and some test data from the
Serra Santa Luzia site.

.. ipython:: python
    :okwarning:

    import numpy as np
    import pandas as pd
    import xarray as xr
    import matplotlib.pyplot as plt

    import windkit as wk
    import pywasp as pw

.. ipython:: python
    :okwarning:

    tutorial_data = wk.get_tutorial_data("serra_santa_luzia")

    bwc = tutorial_data.bwc
    elev_map = tutorial_data.elev
    rgh_map = tutorial_data.rgh

    topo_map = pw.wasp.TopographyMap(elev_map, rgh_map)

    gwc = pw.wasp.generalize(bwc, topo_map)

We have used this site in the :ref:`wasp_flow_model` page, but here we will
go further and include some turbines. Their locations are stored in a ``.csv`` file,
so we will create an `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:

.. ipython:: python
    :okwarning:

    wtg_locs = pd.read_csv('./source/tutorials/data/turbine_positions.csv')

    output_locs =  wk.spatial.create_point(
        wtg_locs.Easting.values,
        wtg_locs.Northing.values,
        wtg_locs['Hub height'].values,
        crs="EPSG:32629"
    )

    fig, ax = plt.subplots()
    ax.set(xlim=(510000, 521500), ylim=(4616000, 4627500))
    elev_map.plot("elev", ax=ax);
    @savefig turbine_positions.png align=center width=6in
    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 first estimate the wind climate at the hub height of each turbine:

.. ipython:: python
    :okwarning:

    wwc = pw.wasp.predict_wwc(bwc, topo_map, output_locs)
    print(wwc)

Gross AEP
=========

Before we estimate the AEP, we need to get a Wind Turbine Generator (wtg) corresponding to the turbine models used.
We use :py:func:`windkit.read_wtg` to read the ``.wtg`` file.

.. ipython:: python
    :okwarning:

    wtg = wk.read_wtg("./source/tutorials/data/Bonus_1_MW.wtg")
    print(wtg)

The ``wtg`` holds power and thrust curves for a number of wind speeds and operational modes (in this case, just one operational mode).
The ``wtg`` also contains extra information, which can be either mode-specific, such as ``air_density``, ``cut_in_wind_speed``, ``cut_out_wind_speed``, ``rated_wind_speed``,
or dimensionless, such as ``name``, ``manufacturer``, ``rotor_diameter``, ``hub_height``, and ``regulation_type``.

We can now estimate the gross AEP using :py:func:`pywasp.gross_aep`, passing at least a ``wwc`` and a ``wtg`` object:

.. ipython:: python
    :okwarning:

    aep = pw.gross_aep(wwc, wtg)

    print(aep)
    print(aep["gross_aep"].values.sum(), "GWh")

Passing simply a ``wwc`` and a ``wtg`` works well when your climate is predicted for exactly the hub heights and locations
of the same turbine type for which you want to estimate AEP.
However, it may be the case that a project has different turbines across the same wind farm.
To handle this, a more flexible approach is needed. Instead of a single ``wtg``, a dictionary of wtg's and a dataset of
turbine locations, heights, turbine ids, group ids, and wtg keys that
map to the wtg's in the dictionary should be passed. In the example below, we show how to do this and how to make use of the extra keyword arguments in the :py:func:`pywasp.gross_aep` function.

.. figure:: ../_static/yield_workflow_2_turbines.png
    :align: center

.. ipython:: python
    :okwarning:

    v90_wtg  = wk.read_wtg("./source/tutorials/data/V90-3.0MW_VCRS_60hz.wtg")
    wtg_dict = {"Bonus_1_MW":wtg, "V90-3.0MW_VCRS_60hz":v90_wtg }

    output_locs_updated  = output_locs.copy(deep=True)
    output_locs_updated["height"].loc[dict(point=range(10, 15))] = 80

    wwc_updated  = pw.wasp.predict_wwc(bwc, topo_map, output_locs_updated)

    wind_turbines = wk.create_wind_turbines_from_arrays(
        west_east=output_locs_updated.west_east.values,
        south_north=output_locs_updated.south_north.values,
        height=output_locs_updated.height.values,
        wtg_keys=["Bonus_1_MW"] * 10 + ["V90-3.0MW_VCRS_60hz"] * 5,
        group_ids=np.zeros(15),
        turbine_ids=np.arange(15),
    )

    aep_updated = pw.gross_aep(
        wwc_updated,
        wtg_dict,
        wind_turbines,
        use_sectors=True,
        air_density_correction="infer",
        air_density=None)

    print(aep_updated)
    print(aep_updated["gross_aep"].values.sum(), "GWh")

In the example above, we introduce an additional turbine type, the "V90-3.0MW_VCRS_60hz",
which has a different hub height compared to the existing "Bonus_1_MW" turbines.
Specifically, we effectively replace the last five turbines with the new turbine model.

Therefore, we must first predict the wind climate at the correct hub heights.
Initially, the output locations assume a hub height of 50 m for the first ten turbines (Bonus_1_MW).
However, since the newly introduced turbines operate at 80 m, we update the output locations accordingly to reflect this difference.

Once the wind climate has been downscaled to the revised hub heights,
we proceed with creating a wind turbines dataset that correctly maps turbine types to their respective locations.
This dataset, represented by the ``wind_turbines`` object, contains:
- turbine locations
- hub heights
- wtg keys (to name the turbine model)
- group ids (to categorize the different groups of turbines within a wind farm)
- turbine ids (unique integer identifiers for each turbine)

When working with multiple turbine types, it is essential to
make sure that the ``wtg_keys`` name matches the dictionary key passed in the ``wtg_dict`` object.
Notice that the last 5 turbines in the ``wind_turbines`` object have a value of "V90-3.0MW_VCRS_60hz" in the ``wtg_keys`` attribute.
There is no limit to the number of turbine models that can be used in this approach.

Additionally, turbine IDs must be unique integers, typically assigned based on the order of
west-east and south-north locations for consistency.

Moreover, extra keyword arguments in the :py:func:`pywasp.gross_aep` function can be used to further customize the AEP estimation.
In this example, we use the default values:

    - ``use_sectors``: Defaults to ``True`` to use sector-wise "A" and "k". If set to ``False``, "A_combined" and "k_combined" are used for the AEP calculation.
    - ``air_density_correction``: Defaults to ``"infer"``, which means that air density is inferred from the wind climate and the ``wtg`` power curve gets corrected to account for site-specific air density.
      If set to ``"none"``, no correction is applied to the power curve, and the default air density mode is used to calculate produced power.
    - ``air_density``: Can be passed with a value to override the inferred air density.

Potential AEP
=============================

To add wind farm effects, such as wakes and blockage, to the AEP estimation and obtain the
"potential AEP", PyWAsP integrates with :doc:`PyWake <py_wake:index>`. PyWAsP defines a number
of named wind farm models, see :py:func:`pywasp.potential_aep`, including "PARK2_onshore" and "PARK2_offshore",
that are identical to the wake model in the WAsP GUI. Here we will use the onshore version:

.. ipython:: python
    :okwarning:

    aep_wakes = pw.potential_aep(wwc,
        wtg,
        wind_farm_model="PARK2_onshore",
        turbulence_intensity=0.0
    )

    print(aep_wakes)
    print(aep_wakes["potential_aep"].values.sum(), "GWh")

We can also calculate the potential AEP for the updated turbine models and locations:

.. ipython:: python
    :okwarning:

    aep_wakes_updated = pw.potential_aep(
        wwc_updated,
        wtg_dict,
        wind_turbines,
        wind_farm_model="PARK2_onshore",
        turbulence_intensity=0.0,
        shear = None,
        air_density_correction="infer"
    )

    print(aep_wakes_updated)
    print(aep_wakes_updated["potential_aep"].values.sum(), "GWh")

Notice that :py:func:`pywasp.potential_aep` needs to be passed with extra keyword arguments ``turbulence_intensity`` (mandatory) and ``shear`` (optional).
Moreover, the same characteristics regarding the air density correction arguments in :py:func:`pywasp.gross_aep` function apply for :py:func:`pywasp.potential_aep`,
where again, there is no limit to the number of turbine models that can be used.

From the complex case above we can see how, even though, we achieve higher power with the larger turbines,
these also generate a considerable amount of wake losses.

.. figure:: ../_static/wakes_plot.png
    :align: center

Wind farm flow map
==================

To map an area around a wind farm, PyWAsP has :py:func:`pywasp.wind_farm_flow_map`, which
takes a resource grid as input and enriches it with wind farm variables including
effects of nearby turbines. This function gives a detailed visualization of the wind farm effects.
Currently, the use of different turbine models populating the wind farm is not supported.
Thus, the example below corresponds to the case where the wind farm is fully populated with the "Bonus_1_MW" turbine.

.. ipython:: python
    :okwarning:

    output_grid = wk.spatial.create_dataset(
        np.linspace(511700.0, 517700.0, 31),
        np.linspace(4619500.0, 4625500.0, 31),
        [50.0],
        crs="EPSG:32629",
    )

    wwc_grid = pw.wasp.downscale(gwc, topo_map, output_grid, interp_method="nearest")

    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),
    )

    flow_map = pw.wind_farm_flow_map(
        wwc_grid,
        wtg,
        wind_turbines,
        output_grid,
        wind_farm_model="PARK2_onshore",
        turbulence_intensity=0.0
    )

    print(flow_map)

Notice that in :py:func:`wind_farm_flow_map`, the first argument, namely the wind climate, is a gridded wind climate.
This is because we need to wrap the wind turbines of the wind farm in a large enough cuboid to capture the upstream and downstream effects.
By modifying the ``output_grid`` ``x`` and ``y`` linespace, a finer or coarser resolution in the grid can be obtained.
Remember to set the cuboid height to the hub height of the turbines, in this case 50 m.
Additionally, the ``wind_turbines`` object is also mandatory to be passed, as the function needs to know the exact location of the turbines within the cuboid.

The :py:func:`pywasp.wind_farm_flow_map` function outputs an updated `xarray.Dataset` with the wind speed and AEP deficits. We can later use these to generate plots,
for example, the wind speed deficit at one directional sector to visualize the wake effects:

.. ipython:: python
    :okwarning:

    @savefig wspd_deficit.png align=center width=6in
    flow_map["wspd_deficit_sector"].sel(sector=210).plot(vmax=0.1)

Net AEP
=============================

To add other losses after wake effects towards the AEP estimation and obtain the "net AEP", we can do so by importing a loss table with the `windkit` function :py:func:`wk.get_loss_table`.
For now, only one table is available, that is a "dtu-default" loss table, which can be modified and personalized by the user.
`pywasp` :doc:`Tutorial 6 <../tutorials/tutorial6_nb>`, shows an extended use case of the loss module.
Once the loss table is imported, the :py:func:`pywasp.net_aep` function can be used to estimate the net AEP.

.. ipython:: python
    :okwarning:

    loss_table = wk.get_loss_table()

    print(wk.loss_table_summary(loss_table))

    net_aep = pw.net_aep(loss_table, aep_wakes)

    print(net_aep["net_aep"].sum().values, "GWh")

Px yield
=============================
The P90 yield, that is the AEP bankable number with 90% probability of being achieved or surpassed, can be estimated in a similar way,
but has some extra steps involved in the process. For a more detailed overview, please refer to :doc:`Tutorial 6 <../tutorials/tutorial6_nb>` to see the full process and usage of these functions.

In this simplified approach, first an uncertainty table must be loaded, and if desired, personalized.
Next, we must estimate the sensitivity factor, which tells the model how sensitive the AEP is to wind speed. For this, the :py:func:`pywasp.estimate_sensitivity_factor` function is used.
Note that this function requires the wind climate, the wind turbine generator, and a wind perturbation factor.
Ideally, for the ``wind_perturbation_factor``, the wind-speed dimensional standard deviation should be passed. This can be retrieved using the windkit function :py:func:`windkit.total_uncertainty` and later multiplying by the local wind speed.

Finally, the P90 yield is estimated with the :py:func:`pywasp.px_aep` function, although one could also estimate any other AEP associated with a probability, such as P50, P75, etc., by changing the percentile argument.

.. ipython:: python
    :okwarning:

    uncertainty_table = wk.get_uncertainty_table()

    Uave = wk.mean_wind_speed(wwc, bysector=False)
    _, wind_uncertainty, _ = wk.total_uncertainty(uncertainty_table)
    sigma_wind_dimensional =  wind_uncertainty/100 * Uave.values.mean()
    print(sigma_wind_dimensional, "m/s")

    sf = pw.estimate_sensitivity_factor(wwc, wtg, wind_perturbation_factor=sigma_wind_dimensional)
    print("The sensitivity factor is:", sf)

    p90 = pw.px_aep(uncertainty_table, net_aep, sensitivity_factor=sf, percentile=90)
    print("The P90 is:", p90["P_90_aep"].values.sum(), "GWh")

PyWake integration in PyWAsP
============================

PyWAsP can use :doc:`PyWake <py_wake:index>` wind farm models to estimate the potential AEP. This is done by passing the name of the PyWake wind farm model to the :py:func:`pywasp.potential_aep` and :py:func:`pywasp.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 :py:func:`functools.partial`, and use it for adding wind farm effects in PyWAsP:

.. ipython:: python
    :okwarning:

    from functools import partial
    import py_wake

    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,
    )

    aep_wakes_personalized_model = pw.potential_aep(
        wwc,
        wtg,
        wind_farm_model=wind_farm_model,
        turbulence_intensity=0.0
    )

    print(aep_wakes_personalized_model)
    print(aep_wakes_personalized_model["potential_aep"].values.sum(), "GWh")

Calculating AEP including wake effects
---------------------------------------

In PyWAsP, the "park" method is used in :py:meth:`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:

1. 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.
2. The wind farm simulation results from py_wake contains the variable two different wind speed variables. One is ``WS``, which is inhomogeneous background (free) wind speed used to obtain the probability
   density from the indivudal sector-wise weibull distributions. Another one is variable ``WS_eff``, which contains the effective wind speed after including wind farm effects.
3. 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.