Add option to calculate AEP as described in the European Wind Atlas, page 95....

Add option to calculate AEP as described in the European Wind Atlas, page 95. This method calculate the power x weibull_weight using linear instead of constant power segments

py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py

@@ -4,7 +4,7 @@ from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines, IEA37Site
 from py_wake import NOJ, Fuga
 from py_wake.site._site import UniformSite
 from py_wake.tests import npt
-from py_wake.examples.data.hornsrev1 import HornsrevV80, Hornsrev1Site, wt_x, wt_y
+from py_wake.examples.data.hornsrev1 import HornsrevV80, Hornsrev1Site, wt_x, wt_y, wt9_x, wt9_y
 from py_wake.tests.test_files.fuga import LUT_path_2MW_z0_0_03
 from py_wake.flow_map import HorizontalGrid
 from py_wake.wind_farm_models.engineering_models import All2AllIterative, PropagateDownwind
@@ -269,3 +269,23 @@ def test_huge_farm():
     peak /= 1024**2
     assert peak < 800
     tracemalloc.stop()
+
+
+def test_aep_wind_atlas_method():
+    site = Hornsrev1Site()
+
+    wt = IEA37_WindTurbines()
+    wfm = IEA37SimpleBastankhahGaussian(site, wt)
+    x, y = [0], [0]
+    wd = np.arange(360)
+    aep_lps = wfm(x, y, wd=wd, ws=np.arange(3, 27)).aep(linear_power_segments=True)
+    aep = wfm(x, y, wd=wd, ws=np.r_[3, np.arange(3.5, 27)]).aep()
+    if 0:
+        plt.plot(aep_lps.ws, np.cumsum(aep_lps.sum(['wt', 'wd'])), '.-', label='Linear power segments')
+        plt.plot(aep.ws, np.cumsum(aep.sum(['wt', 'wd'])), '.-', label='Constant power segments')
+        plt.ylabel('Cumulated AEP [GWh]')
+        plt.xlabel('Wind speed [m/s]')
+        plt.legend()
+        plt.show()
+    npt.assert_almost_equal(aep_lps.sum(), 16.73490444)
+    npt.assert_almost_equal(aep.sum(), 16.69320343)

py_wake/utils/weibull.py

@@ -5,38 +5,38 @@ def cdf(ws, A, k):
     return 1 - np.exp(-(1 / A * ws) ** k)
 
 
-# def WeightedPower(u, power, A, k):
-#     """Calculate the weibull weighted power
-#
-#     Parameters
-#     ----------
-#     Power : xarray DataArray
-#         Power
-#     Returns
-#     -------
-#     y : array_like
-#         y is
-#
-#     Notes
-#     ------
-#     bla bla
-#     """
-#
-#     # see https://orbit.dtu.dk/en/publications/european-wind-atlas, page 95
-#     def G(alpha, k):
-#         # 1/k times the incomplete gamma function of the two arguments 1/k and alpha^k
-#         # Note, the scipy incomplete gamma function, gammainc, must be multiplied with gamma(k) to match the
-#         # the G function used in the European Wind Atlas
-#         import scipy.special as sc
-#         return 1 / k * sc.gamma(1 / k) * sc.gammainc(1 / k, alpha**k)
-#
-#     u0, u1 = u[:-1], u[1:]
-#     alpha0, alpha1 = (u0 / A), (u1 / A)
-#     p0, p1 = power[..., :-1], power[..., 1:]
-#
-#     res = (p0 * np.exp(-alpha0**k) +  # eq 6.5, p0 * cdf(u0:)
-#            (p1 - p0) / (alpha1 - alpha0) * (G(alpha1, k) - G(alpha0, k)) -  # eq 6.4 linear change p0 to p1
-#            p1 * np.exp(-alpha1**k)
-#            )  # eq 6.5, - p1 * cdf(u1:)
-#
-#     return res
+def WeightedPower(u, power, A, k):
+    """Calculate the weibull weighted power
+
+    Parameters
+    ----------
+    Power : xarray DataArray
+        Power
+    Returns
+    -------
+    y : array_like
+        y is
+
+    Notes
+    ------
+    bla bla
+    """
+
+    # see https://orbit.dtu.dk/en/publications/european-wind-atlas, page 95
+    def G(alpha, k):
+        # 1/k times the incomplete gamma function of the two arguments 1/k and alpha^k
+        # Note, the scipy incomplete gamma function, gammainc, must be multiplied with gamma(k) to match the
+        # the G function used in the European Wind Atlas
+        import scipy.special as sc
+        return 1 / k * sc.gamma(1 / k) * sc.gammainc(1 / k, alpha**k)
+
+    u0, u1 = u[:-1], u[1:]
+    alpha0, alpha1 = (u0 / A), (u1 / A)
+    p0, p1 = power[..., :-1], power[..., 1:]
+
+    res = (p0 * np.exp(-alpha0**k) +  # eq 6.5, p0 * cdf(u0:)
+           (p1 - p0) / (alpha1 - alpha0) * (G(alpha1, k) - G(alpha0, k)) -  # eq 6.4 linear change p0 to p1
+           p1 * np.exp(-alpha1**k)
+           )  # eq 6.5, - p1 * cdf(u1:)
+
+    return res

py_wake/wind_farm_models/wind_farm_model.py

@@ -4,7 +4,7 @@ from py_wake.wind_turbines import WindTurbines
 import numpy as np
 from py_wake.flow_map import FlowMap, HorizontalGrid, FlowBox, YZGrid, Grid
 import xarray as xr
-from py_wake.utils import xarray_utils  # register ilk function @UnusedImport
+from py_wake.utils import xarray_utils, weibull  # register ilk function @UnusedImport
 from numpy import newaxis as na
 
 
@@ -159,10 +159,10 @@ class SimulationResult(xr.Dataset):
                                 ('CT', ct_ilk, 'Thrust coefficient'),
                             ]},
                             coords=coords)
-        self['P'] = localWind.P
-        self['WD'] = localWind.WD
-        self['WS'] = localWind.WS
-        self['TI'] = localWind.TI
+        for n in localWind:
+            if n not in ['ws_lower', 'ws_upper']:
+                self[n] = localWind[n]
+
         if yaw_ilk is None:
             self['Yaw'] = self.Power * 0
         else:
@@ -194,7 +194,8 @@ class SimulationResult(xr.Dataset):
          """
         return self.aep(normalize_probabilities=normalize_probabilities, with_wake_loss=with_wake_loss).ilk()
 
-    def aep(self, normalize_probabilities=False, with_wake_loss=True):
+    def aep(self, normalize_probabilities=False, with_wake_loss=True,
+            hours_pr_year=24 * 365, linear_power_segments=False):
         """Anual Energy Production (sum of all wind turbines, directions and speeds) in GWh.
 
         See aep_ilk
@@ -209,7 +210,28 @@ class SimulationResult(xr.Dataset):
         else:
             power_ilk = self.windFarmModel.windTurbines.power(self.WS.ilk(self.Power.shape), self.type)
 
-        return xr.DataArray(power_ilk * self.P.ilk() / norm * 24 * 365 * 1e-9,
+        if linear_power_segments:
+            s = "The linear_power_segments method "
+            assert all([n in self for n in ['Weibull_A', 'Weibull_k', 'Sector_frequency']]),\
+                s + "requires a site with weibull information"
+            assert normalize_probabilities is False, \
+                s + "cannot be combined with normalize_probabilities"
+            assert np.all(self.Power.isel(ws=0) == 0) and np.all(self.Power.isel(ws=-1) == 0),\
+                s + "requires first wind speed to have no power (just below cut-in)"
+            assert np.all(self.Power.isel(ws=-1) == 0),\
+                s + "requires last wind speed to have no power (just above cut-out)"
+            weighted_power = weibull.WeightedPower(
+                self.ws.values,
+                self.Power.ilk(),
+                self.Weibull_A.ilk(),
+                self.Weibull_k.ilk())
+            aep = weighted_power * self.Sector_frequency.ilk() * hours_pr_year * 1e-9
+            ws = (self.ws.values[1:] + self.ws.values[:-1]) / 2
+            return xr.DataArray(aep, [('wt', self.wt), ('wd', self.wd), ('ws', ws)])
+        else:
+            weighted_power = power_ilk * self.P.ilk() / norm
+
+        return xr.DataArray(weighted_power * hours_pr_year * 1e-9,
                             self.Power.coords,
                             name='AEP [GWh]',
                             attrs={'Description': 'Annual energy production [GWh]'})