Newer
Older
import numpy as np
from numpy import newaxis as na
"""Initialize AEPCalculator
Parameters
----------
site : py_wake.site.Site
windTurbines : WindTurbines
wake_model : WakeModel
"""
self.site = site
self.wake_model = wake_model
self.windTurbines = windTurbines
def _get_defaults(self, x_i, h_i, type_i, wd, ws):
if type_i is None:
type_i = np.zeros_like(x_i, dtype=np.int)
if h_i is None:
h_i = self.windTurbines.hub_height(type_i)
if wd is None:
return h_i, type_i, wd, ws
def _run_wake_model(self, x_i, y_i, h_i=None, type_i=None, wd=None, ws=None):
h_i, type_i, wd, ws = self._get_defaults(x_i, h_i, type_i, wd, ws)
self.WS_eff_ilk, self.TI_eff_ilk, self.power_ilk, self.ct_ilk, self.WD_ilk, self.WS_ilk, self.TI_ilk, self.P_ilk =\
self.wake_model.calc_wake(x_i, y_i, h_i, type_i, wd, ws, self.site)
def calculate_AEP(self, x_i, y_i, h_i=None, type_i=None, wd=None, ws=None):
"""Calculate AEP
In addition effective wind speed, turbulence intensity, and the
power, ct and probability is calculated
Parameters
----------
x_i : array_like
X position of wind turbines
y_i : array_like
Y position of wind turbines
h_i : array_like or None, optional
Hub height of wind turbines\n
If None, default, the standard hub height is used
type_i array_like or None, optional
Wind turbine types\n
If None, default, the first type is used (type=0)
wd : int, float, array_like or None
Wind directions(s)\n
If None, default, the wake is calculated for site.default_wd
ws : int, float, array_like or None
Wind speed(s)\n
If None, default, the wake is calculated for site.default_ws
Returns
-------
AEP_GWh_ilk : array_like
AEP in GWh
"""
self._run_wake_model(x_i=x_i, y_i=y_i, h_i=h_i, type_i=type_i, wd=wd, ws=ws)
AEP_GWh_ilk = self.power_ilk * self.P_ilk * 24 * 365 * 1e-9
def calculate_AEP_no_wake_loss(self, x_i, y_i, h_i=None, type_i=None, wd=None, ws=None):
"""Calculate AEP without wake loss
In addition local wind speed, turbulence intensity, and the
power, ct and probability is calculated
Parameters
----------
x_i : array_like
X position of wind turbines
y_i : array_like
Y position of wind turbines
h_i : array_like or None, optional
Hub height of wind turbines\n
If None, default, the standard hub height is used
Wind turbine types\n
If None, default, the first type is used (type=0)
wd : int, float, array_like or None
Wind directions(s)\n
If None, default, the wake is calculated for site.default_wd
ws : int, float, array_like or None
Wind speed(s)\n
If None, default, the wake is calculated for site.default_ws
Returns
-------
AEP_GWh_ilk : array_like
AEP (without wake loss) in GWh
"""
h_i, type_i, wd, ws = self._get_defaults(x_i, h_i, type_i, wd=wd, ws=ws)
# Find local wind speed, wind direction, turbulence intensity and probability
self.WD_ilk, self.WS_ilk, self.TI_ilk, self.P_ilk = self.site.local_wind(x_i=x_i, y_i=y_i, wd=wd, ws=ws)
type_ilk = np.zeros(self.WS_ilk.shape, dtype=np.int) + type_i[:, np.newaxis, np.newaxis]
self.power_ilk = self.windTurbines.power(self.WS_ilk, type_ilk)
AEP_GWh_ilk = self.power_ilk * self.P_ilk * 24 * 365 * 1e-9
def _eff_map(self, type, x_j, y_j, h, x_i, y_i, type_i, h_i, wd=None, ws=None):
h_i, type_i, wd, ws = self._get_defaults(x_i, h_i, type_i, wd, ws)
def f(x, N=500, ext=.2):
ext *= (max(x) - min(x))
return np.linspace(min(x) - ext, max(x) + ext, N)
if x_j is None:
x_j = f(x_i)
if y_j is None:
y_j = f(y_i)
if h is None:
h = np.mean(h_i)
X_j, Y_j = np.meshgrid(x_j, y_j)
x_j, y_j = X_j.flatten(), Y_j.flatten()
if len(x_i) == 0:
_, WS_jlk, _, P_ilk = self.site.local_wind(x_i=x_j, y_i=y_j, wd=wd, ws=ws, wd_bin_size=1)
return X_j, Y_j, WS_jlk, P_ilk
self._run_wake_model(x_i, y_i, h_i, type_i, wd, ws)
h_j = np.zeros_like(x_j) + h
_, WS_jlk, TI_jlk, P_ilk = self.site.local_wind(x_i=x_j, y_i=y_j, wd=wd, ws=ws)
dw_ijl, cw_ijl, dh_ijl, _ = self.site.distances(x_i, y_i, h_i, x_j, y_j, h_j, self.WD_ilk.mean(2))
if type == 'WS':
_eff_jlk = self.wake_model.ws_map(self.WS_ilk, self.WS_eff_ilk, self.TI_ilk, self.TI_ilk, dw_ijl,
cw_ijl, dh_ijl, self.ct_ilk, type_i, WS_jlk)
elif type == 'TI':
_eff_jlk = self.wake_model.ti_map(self.WS_ilk, self.WS_eff_ilk, self.TI_ilk, self.TI_ilk, dw_ijl,
cw_ijl, dh_ijl, self.ct_ilk, type_i, TI_jlk)
return X_j, Y_j, _eff_jlk, P_ilk
def wake_map(self, x_j=None, y_j=None, h=None, wt_x=[], wt_y=[], wt_type=None, wt_height=None, wd=None, ws=None):
"""Calculate wake(effective wind speed) map
In addition local wind speed, turbulence intensity, and the
power, ct and probability is calculated
Parameters
----------
x_j : array_like or None, optional
X position map points
y_j : array_like
Y position of map points
h : int, float or None, optional
Height of wake map\n
If None, default, the mean hub height is used
wt_x : array_like, optional
X position of wind turbines
wt_y : array_like, optional
Y position of wind turbines
wt_type : array_like or None, optional
Type of the wind turbines
wt_height : array_like or None, optional
Hub height of the wind turbines\n
If None, default, the standard hub height is used
wd : int, float, array_like or None
Wind directions(s)\n
If None, default, the wake is calculated for site.default_wd
ws : int, float, array_like or None
Wind speed(s)\n
If None, default, the wake is calculated for site.default_ws
Returns
-------
X_j : array_like
2d array of map x positions
Y_j : array_like
2d array of map y positions
WS_eff_avg : array_like
2d array of average effective local wind speed taking into account
the probability of wind direction and speed
See Also
--------
plot_wake_map
X_j, Y_j, WS_eff_jlk, P_ilk = self._eff_map(
'WS', x_j, y_j, h, wt_x, wt_y, wt_type, wt_height, wd, ws)
if P_ilk.sum() > 0:
WS_eff_jlk = WS_eff_jlk * (P_ilk / P_ilk.sum((1, 2)))
return X_j, Y_j, WS_eff_jlk.sum((1, 2)).reshape(X_j.shape)
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
def ti_map(self, x_j=None, y_j=None, h=None, wt_x=[], wt_y=[], wt_type=None, wt_height=None, wd=None, ws=None):
"""Calculate turbulence intensity map
Parameters
----------
x_j : array_like or None, optional
X position map points
y_j : array_like
Y position of map points
h : int, float or None, optional
Height of wake map\n
If None, default, the mean hub height is used
wt_x : array_like, optional
X position of wind turbines
wt_y : array_like, optional
Y position of wind turbines
wt_type : array_like or None, optional
Type of the wind turbines
wt_height : array_like or None, optional
Hub height of the wind turbines\n
If None, default, the standard hub height is used
wd : int, float, array_like or None
Wind directions(s)\n
If None, default, the wake is calculated for site.default_wd
ws : int, float, array_like or None
Wind speed(s)\n
If None, default, the wake is calculated for site.default_ws
Returns
-------
X_j : array_like
2d array of map x positions
Y_j : array_like
2d array of map y positions
WS_eff_avg : array_like
2d array of average effective local wind speed taking into account
the probability of wind direction and speed
See Also
--------
plot_wake_map
"""
X_j, Y_j, TI_eff_jlk, P_ilk = self._eff_map(
'TI', x_j, y_j, h, wt_x, wt_y, wt_type, wt_height, wd, ws)
if P_ilk.sum() > 0:
TI_eff_jlk = TI_eff_jlk * (P_ilk / P_ilk.sum((1, 2)))
return X_j, Y_j, TI_eff_jlk.sum((1, 2)).reshape(X_j.shape)
def plot_wake_map(self, x_j=None, y_j=None, h=None, wt_x=[], wt_y=[], wt_type=None, wt_height=None,
wd=None, ws=None, ax=None, levels=100):
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
"""Plot wake(effective wind speed) map
Parameters
----------
x_j : array_like or None, optional
X position map points
y_j : array_like
Y position of map points
h : int, float or None, optional
Height of wake map\n
If None, default, the mean hub height is used
wt_x : array_like, optional
X position of wind turbines
wt_y : array_like, optional
Y position of wind turbines
wt_type : array_like or None, optional
Type of the wind turbines
wt_height : array_like or None, optional
Hub height of the wind turbines\n
If None, default, the standard hub height is used
wd : int, float, array_like or None
Wind directions(s)\n
If None, default, the wake is calculated for site.default_wd
ws : int, float, array_like or None
Wind speed(s)\n
If None, default, the wake is calculated for site.default_ws
ax : pyplot or matplotlib axes object, default None
Axes to plot on
levels : int or array_like
levels for pyplot.contourf
"""
import matplotlib.pyplot as plt
if ax is None:
ax = plt.gca()
X, Y, Z = self.wake_map(x_j, y_j, h, wt_x, wt_y, wt_type, wt_height, wd, ws)
c = ax.contourf(X, Y, Z, levels, cmap='Blues_r')
plt.colorbar(c, label='wind speed [m/s]')
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
def aep_map(self, x_j=None, y_j=None, type_j=None, wt_x=[], wt_y=[], wt_type=None, wt_height=None, wd=None, ws=None):
"""Calculate AEP map
The map represents the of AEP produced by a new turbine at the specified positions
Parameters
----------
x_j : array_like or None, optional
X position map points (potential turbine positions)
y_j : array_like
Y position of map points (potential turbine positions)
type_j : int, float or None, optional
Type of potential turbine positions\n
If None, default, first turbine type(0) is used
wt_x : array_like, optional
X position of the current wind turbines
wt_y : array_like, optional
Y position of the current wind turbines
wt_type : array_like or None, optional
Type of the current wind turbines
wt_height : array_like or None, optional
Hub height of the current wind turbines\n
If None, default, the standard hub height is used
wd : int, float, array_like or None
Wind directions(s)\n
If None, default, the wake is calculated for site.default_wd
ws : int, float, array_like or None
Wind speed(s)\n
If None, default, the wake is calculated for site.default_ws
Returns
-------
X_j : array_like
2d array of map x positions
Y_j : array_like
2d array of map y positions
WS_eff_avg : array_like
2d array of average effective local wind speed taking into account
the probability of wind direction and speed
"""
h_j = self.windTurbines.hub_height(type_j)
X_j, Y_j, WS_eff_jlk, P_jlk = self._eff_map('WS', x_j, y_j, h_j, wt_x, wt_y, wt_type, wt_height, wd, ws)
P_jlk /= P_jlk.sum((1, 2)) # AEP if wind only comes from specified wd and ws
# power_jlk = self.windTurbines.power_func(WS_eff_jlk, type_j)
# aep_jlk = power_jlk * P_jlk * 24 * 365 * 1e-9
# return X_j, Y_j, aep_jlk.sum((1, 2)).reshape(X_j.shape)
# same as above but requires less memory
return X_j, Y_j, ((self.windTurbines.power(WS_eff_jlk, type_j) * P_jlk).sum((1, 2)) * 24 * 365 * 1e-9).reshape(X_j.shape)
from py_wake.examples.data.iea37 import iea37_path
from py_wake.examples.data.iea37._iea37 import IEA37Site
from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines
from py_wake.wake_models import NOJ
# setup site, turbines and wakemodel
site = IEA37Site(16)
x, y = site.initial_position.T
windTurbines = IEA37_WindTurbines(iea37_path + 'iea37-335mw.yaml')
wake_model = NOJ(windTurbines)
# calculate AEP
aep_calculator = AEPCalculator(site, windTurbines, wake_model)
aep = aep_calculator.calculate_AEP(x, y)[0].sum()
print(aep_calculator.WS_eff_ilk.shape)
# plot wake map
import matplotlib.pyplot as plt
aep_calculator.plot_wake_map(wt_x=x, wt_y=y, wd=[0], ws=[9])
plt.title('AEP: %.2f GWh' % aep)
windTurbines.plot(x, y)
plt.show()
if __name__ == '__main__':
main()