Newer
Older
import numpy as np
from abc import ABC, abstractmethod
"""
suffixs:
- i: Local point (wind turbines)
- j: Local point (downstream turbines or positions)
- l: Wind directions
- k: Wind speeds
class Site(ABC):
def __init__(self):
self.default_ws = np.arange(3, 26)
self.default_wd = np.arange(360)
@abstractmethod
def local_wind(self, x_i, y_i, h_i=None, wd=None, ws=None, wd_bin_size=None, ws_bin_size=None):
"""Local free flow wind conditions
Parameters
----------
Local x coordinate
y_i : array_like
Local y coordinate
h_i : array_like, optional
Local h coordinate, i.e., heights above ground
Global wind speed(s). Override self.default_ws
wd_bin_size : int or float, optional
Size of wind direction bins. default is size between first and
second element in default_wd
ws_bin_size : int or float, optional
Size of wind speed bins. default is size between first and
second element in default_ws
Returns
-------
WD_ilk : array_like
local free flow wind directions
WS_ilk : array_like
local free flow wind speeds
TI_ilk : array_like
local free flow turbulence intensity
Probability/weight
"""
pass
@abstractmethod
def probability(self, x_i, y_i, h_i, WD_lk, WS_lk, wd_bin_size, ws_bin_size):
"""Probability of wind situation (wind speed and direction)
Parameters
----------
x_i : array_like
Local x coordinate
y_i : array_like
Local y coordinate
h_i : array_like
Local height
Wind speed
wd_bin_size : int or float
size of wind direction sectors
ws_bin_size : int or float
size of wind speed bins
Returns
-------
P_ilk : float or array_like
Probability of wind speed and direction at local positions
pass
@abstractmethod
def distances(self, src_x_i, src_y_i, src_h_i, dst_x_j, dst_y_j, dst_h_j, wd_il):
"""Calculate down/crosswind distance between source and destination points
src_x_i : array_like
Source x position
src_y_i : array_like
Source y position
Source height above ground level
dst_x_j : array_like
Destination x position
dst_y_j : array_like
Destination y position
dst_h_j : array_like
Destination height above ground level
wd_il : array_like, shape (#src, #wd)
Local wind direction at the source points for all global wind directions
hcw_ijl : array_like
horizontal cross wind distances. Positive when the wind is in your
back and the turbine lies on the left.
dh_ijl : array_like
vertical distances
dw_order_indices_l : array_like
indices that gives the downwind order of source points
"""
pass
def wt2wt_distances(self, x_i, y_i, h_i, wd_il):
return self.distances(x_i, y_i, h_i, x_i, y_i, h_i, wd_il)
@abstractmethod
def elevation(self, x_i, y_i):
"""Local terrain elevation (height above mean sea level)
Parameters
----------
x_i : array_like
Local x coordinate
y_i : array_like
Local y coordinate
Returns
-------
elevation : array_like
"""
pass
def wd_bin_size(self, wd, wd_bin_size=None):
if wd_bin_size is not None:
return wd_bin_size
else:
return 360 / len(np.atleast_1d(wd))
def ws_bin_size(self, ws, ws_bin_size=None):
if ws_bin_size is None:
if hasattr(ws, '__len__') and len(ws) > 1:
return ws[1] - ws[0]
else:
return 1
else:
return ws_bin_size
def plot_ws_distribution(self, x=0, y=0, h=70, wd=[0], include_wd_distribution=False, ax=None):
"""Plot wind speed distribution
Parameters
----------
x : int or float
Local x coordinate
y : int or float
Local y coordinate
h : int or float
Local height above ground
wd : int or array_like
Wind direction(s) (one curve pr wind direction)
include_wwd_distributeion : bool, default is False
If true, the wind speed probability distributions are multiplied by
the wind direction probability. The sector size is set to 360 / len(wd).
This only makes sense if the wd array is evenly distributed
ax : pyplot or matplotlib axes object, default None
"""
if ax is None:
import matplotlib.pyplot as plt
ax = plt
ws = np.arange(0.05, 30.05, .1)
x = np.atleast_1d(x)
y = np.atleast_1d(y)
wd = np.atleast_1d(wd)
for wd_ in wd:
wd_bin_size = 360 / len(wd)
if include_wd_distribution:
v = wd_bin_size / 2
wd_l = np.arange(wd_ - v, wd_ + v) % 360
WD_lk, WS_lk = np.meshgrid(wd_l, ws, indexing='ij')
p = self.probability(x, y, h, WD_lk, WS_lk, wd_bin_size=1, ws_bin_size=.1).sum((0, 1))
lbl = r"Wind direction: %d$\pm$%s deg" % (wd_, (int(v), v)[(wd_bin_size % 2) != 0])
else:
# WD_lk = np.array([wd_])[:, na]
# WS_lk = ws
WD_lk, WS_lk = np.meshgrid([wd_], ws, indexing='ij')
p = self.probability(x, y, h, WD_lk, WS_lk, wd_bin_size=wd_bin_size, ws_bin_size=.1)[0, 0]
p /= p.sum()
lbl = "Wind direction: %d deg" % (wd_)
ax.plot(ws, p * 10, label=lbl)
ax.xlabel('Wind speed [m/s]')
ax.ylabel('Probability')
def plot_wd_distribution(self, x=0, y=0, h=70, n_wd=12, ws_bins=None, ax=None):
"""Plot wind direction (and speed) distribution
Parameters
----------
x : int or float
Local x coordinate
y : int or float
Local y coordinate
h : int or float
Local height above ground
n_wd : int
Number of wind direction sectors
ws_bins : None, int or array_like, default is None
Splits the wind direction sector pies into different colors to show
the probability of different wind speeds\n
If int, number of wind speed bins in the range 0-30\n
If array_like, limits of the wind speed bins limited by ws_bins,
e.g. [0,10,20], will show 0-10 m/s and 10-20 m/s
ax : pyplot or matplotlib axes object, default None
"""
if ax is None:
import matplotlib.pyplot as plt
ax = plt
x = np.atleast_1d(x)
y = np.atleast_1d(y)
h = np.atleast_1d(h)
wd = np.linspace(0, 360, n_wd, endpoint=False)
theta = wd / 180 * np.pi
if not ax.__class__.__name__ == 'PolarAxesSubplot':
if hasattr(ax, 'subplot'):
ax = ax.figure.add_subplot(111, projection='polar')
ax.set_theta_direction(-1)
ax.set_theta_offset(np.pi / 2.0)
s = 360 / n_wd
if ws_bins is None:
WD_lk, WS_lk = np.meshgrid(np.arange(-s / 2, s / 2) + 1, [100], indexing='ij')
p = [self.probability(x_i=x, y_i=y, h_i=h, WD_lk=(WD_lk + wd_) % 360,
WS_lk=WS_lk,
wd_bin_size=1, ws_bin_size=200).sum() for wd_ in wd]
ax.bar(theta, p, width=s / 180 * np.pi, bottom=0.0)
else:
if not hasattr(ws_bins, '__len__'):
ws_bins = np.linspace(0, 30, ws_bins)
else:
ws_bins = np.asarray(ws_bins)
ws = ((ws_bins[1:] + ws_bins[:-1]) / 2)
ws_bin_size = ws[1] - ws[0]
WD_lk, WS_lk = np.meshgrid(np.arange(-s / 2, s / 2) + 1, ws, indexing='ij')
p = [self.probability(x_i=x, y_i=y, h_i=h, WD_lk=(WD_lk + wd_) % 360, WS_lk=WS_lk,
wd_bin_size=1, ws_bin_size=ws_bin_size).sum((0, 1)) for wd_ in wd]
cum_p = np.cumsum(p, 1).T
start_p = np.vstack([np.zeros_like(cum_p[:1]), cum_p[:-1]])
for ws1, ws2, p_ws1, p_ws2 in zip(ws_bins[:-1], ws_bins[1:], start_p, cum_p):
ax.bar(theta, p_ws2 - p_ws1, width=s / 180 * np.pi, bottom=p_ws1, label="%s-%s m/s" % (ws1, ws2))
ax.legend(bbox_to_anchor=(1.15, 1.1))
ax.set_rlabel_position(-22.5) # Move radial labels away from plotted line
ax.grid(True)
"""Site with uniform (same wind over all, i.e. flat uniform terrain) and
constant wind speed probability of 1. Only for one fixed wind speed
"""
def __init__(self, p_wd, ti, ws=12, interp_method='piecewise', alpha=None, h_ref=None):
super().__init__()
self.h_ref = h_ref
self.default_ws = ws
self.ti = Sector2Subsector(np.atleast_1d(ti), interp_method=interp_method)
self.alpha = Sector2Subsector(np.atleast_1d(alpha), interp_method=interp_method)
self.p_wd = Sector2Subsector(p_wd / np.sum(p_wd), interp_method=interp_method) / (360 / len(p_wd))
def probability(self, x_i, y_i, h_i, WD_lk, WS_lk, wd_bin_size, ws_bin_size):
P_lk = np.ones_like(WS_lk, dtype=np.float) * self.p_wd[np.round(WD_lk).astype(np.int) % 360] * wd_bin_size
return P_lk[na]
def local_wind(self, x_i=None, y_i=None, h_i=None, wd=None, ws=None, wd_bin_size=None, ws_bin_size=None):
if wd is None:
wd = self.default_wd
if ws is None:
ws = self.default_ws
ws_bin_size = self.ws_bin_size(ws, ws_bin_size)
wd_bin_size = self.wd_bin_size(wd, wd_bin_size)
WD_ilk, WS_ilk = [np.tile(W, (len(x_i), 1, 1)).astype(np.float)
for W in np.meshgrid(wd, ws, indexing='ij')]
WD_index_ilk = np.round(WD_ilk).astype(np.int)
if h_i is not None and self.h_ref is not None and self.alpha is not None:
if not np.all(h_i == self.h_ref):
wind_shear_ratio = (h_i / self.h_ref)[:, na, na] ** self.alpha[WD_index_ilk]
WS_ilk = WS_ilk * wind_shear_ratio
TI_ilk = self.ti[WD_index_ilk]
P_ilk = self.probability(0, 0, 0, WD_ilk[0], WS_ilk[0], wd_bin_size, ws_bin_size)
return WD_ilk, WS_ilk, TI_ilk, P_ilk
def distances(self, src_x_i, src_y_i, src_h_i, dst_x_j, dst_y_j, dst_h_j, wd_il):
wd_l = np.mean(wd_il, 0)
dx_ij, dy_ij, dh_ij = [np.subtract(*np.meshgrid(dst_j, src_i, indexing='ij')).T
for src_i, dst_j in [(src_x_i, dst_x_j),
(src_y_i, dst_y_j),
(src_h_i, dst_h_j)]]
src_x_i, src_y_i = map(np.asarray, [src_x_i, src_y_i])
theta_l = np.deg2rad(90 - wd_l)
cos_l = np.cos(theta_l)
sin_l = np.sin(theta_l)
dw_il = cos_l[na, :] * src_x_i[:, na] + sin_l[na] * src_y_i[:, na]
dw_ijl = -cos_l[na, na, :] * dx_ij[:, :, na] - sin_l[na, na, :] * dy_ij[:, :, na]
hcw_ijl = sin_l[na, na, :] * dx_ij[:, :, na] - cos_l[na, na, :] * dy_ij[:, :, na]
dh_ijl = np.zeros_like(dw_ijl)
dh_ijl[:, :, :] = dh_ij[:, :, na]
# dw_order_indices_l = np.argsort(-dw_il, 0).astype(np.int).T
dw_order_indices_l = np.argsort((dw_ijl > 0).sum(0), 0).T
return dw_ijl, hcw_ijl, dh_ijl, dw_order_indices_l
def elevation(self, x_i, y_i):
return np.zeros_like(x_i)
class UniformWeibullSite(UniformSite):
"""Site with uniform (same wind over all, i.e. flat uniform terrain) and
weibull distributed wind speed
"""
def __init__(self, p_wd, a, k, ti, interp_method='nearest', alpha=None, h_ref=None):
"""Initialize UniformWeibullSite
Parameters
----------
p_wd : array_like
Probability of wind direction sectors
a : array_like
Weilbull scaling parameter of wind direction sectors
k : array_like
Weibull shape parameter
ti : float or array_like
Turbulence intensity
interp_method : 'nearest', 'linear' or 'spline'
p_wd, a, k, ti and alpha are interpolated to 1 deg sectors using this
method
alpha : float or array_like
Power shear profile exponent of the wind direction sectors
h_ref : int or float
Reference height for the power shear profile
Notes
------
The wind direction sectors will be: [0 +/- w/2, w +/- w/2, ...]
where w is 360 / len(p_wd)
"""
super().__init__(p_wd, ti, interp_method=interp_method, alpha=alpha, h_ref=h_ref)
self.a = Sector2Subsector(a, interp_method=interp_method)
self.k = Sector2Subsector(k, interp_method=interp_method)
def cdf(ws, A=A, k=k):
return 1 - np.exp(-(ws / A) ** k)
def probability(self, x_i, y_i, h_i, WD_lk, WS_lk, wd_bin_size, ws_bin_size):
i_lk = np.round(WD_lk).astype(np.int) % 360
p_wd = self.p_wd[i_lk] * wd_bin_size
if wd_bin_size == 360:
p_wd /= p_wd.sum(0)
P_lk = self.weibull_weight(WS_lk, self.a[i_lk], self.k[i_lk], ws_bin_size) * p_wd
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
def Sector2Subsector(para, axis=-1, wd_binned=None, interp_method='piecewise'):
""" Expand para on the wind direction dimension, i.e., increase the nubmer
of sectors (sectors to subsectors), by interpolating between sectors, using
specified method.
Parameters
----------
para : array_like
Parameter to be expand, it can be sector-wise Weibull A, k, frequency.
axis : integer
Denotes which dimension of para corresponds to wind direction.
wd_binned : array_like
Wind direction of subsectors to be expanded to.
inter_method : string
'piecewise'/'linear'/'spline', based on interp1d in scipy.interpolate,
'spline' means cubic spline.
--------------------------------------
Note: the interpolating method for sector-wise Weibull distributions and
joint distribution of wind speed and wind direction is referred to the
following paper:
Feng, J. and Shen, W.Z., 2015. Modelling wind for wind farm layout
optimization using joint distribution of wind speed and wind direction.
Energies, 8(4), pp.3075-3092. [https://doi.org/10.3390/en8043075]
"""
if wd_binned is None:
wd_binned = np.arange(360)
para = np.array(para)
num_sector = para.shape[axis]
wd_sector = np.linspace(0, 360, num_sector, endpoint=False)
try:
interp_index = ['nearest', 'piecewise', 'linear', 'spline'].index(interp_method)
interp_kind = ['nearest', 'nearest', 'linear', 'cubic'][interp_index]
except ValueError:
raise NotImplementedError(
'interp_method={0} not implemeted yet.'.format(interp_method))
wd_sector_extended = np.hstack((wd_sector, 360.0))
para_sector_extended = np.concatenate((para, para.take([0], axis=axis)),
axis=axis)
if interp_kind == 'cubic' and len(wd_sector_extended) < 4:
interp_kind = 'linear'
f_interp = interpolate.interp1d(wd_sector_extended, para_sector_extended,
kind=interp_kind, axis=axis)
para_expanded = f_interp(wd_binned % 360)
return para_expanded
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
f = [0.035972, 0.039487, 0.051674, 0.070002, 0.083645, 0.064348,
0.086432, 0.117705, 0.151576, 0.147379, 0.10012, 0.05166]
A = [9.176929, 9.782334, 9.531809, 9.909545, 10.04269, 9.593921,
9.584007, 10.51499, 11.39895, 11.68746, 11.63732, 10.08803]
k = [2.392578, 2.447266, 2.412109, 2.591797, 2.755859, 2.595703,
2.583984, 2.548828, 2.470703, 2.607422, 2.626953, 2.326172]
ti = .1
h_ref = 100
alpha = .1
site = UniformWeibullSite(f, A, k, ti, alpha=alpha, h_ref=h_ref)
x_i = y_i = np.arange(5)
wdir_lst = np.arange(0, 360, 90)
wsp_lst = np.arange(1, 20)
WD_ilk, WS_ilk, TI_ilk, P_ilk = site.local_wind(x_i=x_i, y_i=y_i, wd=wdir_lst, ws=wsp_lst)
import matplotlib.pyplot as plt
site.plot_ws_distribution(0, 0, wdir_lst)
plt.figure()
z = np.arange(1, 100)
u = [site.local_wind(x_i=[0], y_i=[0], h_i=[z_], wd=0, ws=10)[1][0][0] for z_ in z]
plt.plot(u, z)
plt.xlabel('Wind speed [m/s]')
plt.ylabel('Height [m]')
plt.show()
if __name__ == '__main__':
main()