Optimization of the distribution of wind speeds using convexly combined Weibull densities

Abstract

This paper presents a new approach for the determination of the wind speed distribution based on wind speed data. This approach is based on the fact that, in general, wind speed distributions restricted to seasons of year or months are different. Therefore, instead of one Weibull density function, a convex combination of Weibull density functions is considered for a calendar year. This model improves the maximum likelihood of the estimated wind speed distribution. Numerical results including a Kolmogorov–Smirnov test are given for a site at Jamaica. Numerical comparisons are carried out for different sites and various known methods for the estimation of the wind speed distribution.

Introduction

For the forecast of the annual revenue of wind power stations, one needs a good estimate of the probability distribution of wind speeds (compare also Wang et al. 2016b; Zhao et al. 2016; Sohoni et al. 2016a). By default, one generally works with a Weibull probability density function (PDF) for wind power potential calculations (e.g. see Hennessey 1977; Bowden et al. 1983; Genc et al. 2005; Sohoni et al. 2016b). Quite often, such an estimated PDF leads to an incorrect prediction of the produced energy so that additional costs may occur (e.g. see Tye et al. 2014). The use of only one Weibull PDF seems to be problematic, and at special sites, e.g. the wind farm Chungtun located at a small island in Taiwan Trait (see Liu et al. 2014 for details), a bimodal mixture Weibull PDF has shown to be more useful (see also Jaramillo and Borja 2004). Other approaches such as the truncated normal-Weibull PDF, the mixture Gamma–Weibull PDF and the mixture truncated normal PDF are known from the special literature (e.g. see Chang 2011; Akpinar and Akpinar 2009; Carta and Mentado 2007; Wang et al. 2016a; Tian Pau 2011; Kollu et al. 2012). Better PDF estimates can be expected, as proposed by Bischoff and Jahn (2016), using convex combinations of different Weibull PDFs. The present paper extends these investigations in such a way that monthly distributions are taken into account. This leads to an improvement of the estimate, which is achieved by a high numerical effort for the solution of a constrained optimization problem with a highly nonlinear objective function.

The goal of this paper is to present this new approach. This method is based on a highly nonlinear optimization problem, which can be solved by standard algorithms of numerical smooth optimization. Since this approach uses much more parameters than the known methods, one gets an improved resulting PDF of wind speeds.

This paper is organized as follows: the next section describes preliminaries, and then, convex combinations of Weibull PDFs are investigated. The algorithmic approach is presented in the fourth section followed by numerical results and a Kolmogorov–Smirnov test. In the last section, numerical comparisons are carried out for known estimation methods applied to different sites.

Preliminaries

Let the random variable V describe the wind speed (in m/s) at an arbitrary site of a wind farm. The PDF of V is very often assumed to be a Weibull density function $f_{k,c}$ given as

$$\begin{aligned} f_{k,c}(v):= {\left\{ \begin{array}{ll} 0 & \quad \text {if } v < 0\\ \frac{k}{c} \bigl ( \frac{v}{c}\bigr ) ^{k-1}e^{-(\frac{v}{c})^k} &\quad \text {if } v \ge 0, \end{array}\right. } \end{aligned}$$

where $k>0$ denotes the so-called Weibull form parameter and $c>0$ denotes the Weibull scale parameter in $\frac{m}{s}$ (compare also Rinne 2008).

In general, the Weibull PDF is estimated on the basis of wind speed forecasts. For instance, for every hour per year one uses a forecast $v_{i}$ with $i\in \{ 1,\ldots ,8760 \}$. Then, these data are used for the determination of the Weibull parameters k and c (e.g. see Gupta et al. 1998). It is outlined by Akdağ and Dinler (2009) that there are different methods for the computation of these parameters. We restrict ourselves to the maximum likelihood estimation, which estimates the parameters k and c in such a way that the data are generated by the corresponding distribution with maximal probability. From a mathematical point of view, one solves the nonlinear optimization problem

$$\begin{aligned} \max _{k,c > 0} \ \prod _{j=1}^n \;f_{k,c}(v_{i_{j}}), \end{aligned}$$

(1)

where one uses only positive values $v_{i_{j}}$ with $j\in \{ 1,\ldots ,n\}$ for some $n\in \{ 1,\ldots ,8760 \}$. So, wind speeds of the type 0 m/s are dropped.

Fig. 1

Histogram of the measured wind speeds with Weibull PDF at site 1 (discretization by 0.3 m/s)

The maximal solutions of the optimization problem (1) are the so-called maximum likelihood estimators of the two Weibull parameters. For simplification, one considers the logarithm of the objective function of problem (1), i.e. one maximizes

$$\begin{aligned} \ln \prod _{j=1}^n \;f_{k,c}(v_{i_{j}})=\sum _{j=1}^n \;\ln f_{k,c}(v_{i_{j}}). \end{aligned}$$

(2)

As an example, Fig. 1 shows the histogram of measured wind speeds and the corresponding Weibull PDF, which is computed using the maximum likelihood estimation for site 1 (see Table 1 for details).

The standard Weibull PDF is certainly not appropriate for site 1. This already shows the known fact that a Weibull PDF is not always the best choice. Wind power potential calculations require a better approximation of the PDF.

The standard Weibull approach has the following disadvantages:

• 1. The data of wind speeds $v_{1},\ldots ,v_{8760}$ are ordered in time. This ordering is not considered in problem (1). Therefore, the structure of the wind profiles is not completely used.

• 2. If there are wind speeds of the form $v_{i}=0$ for some $i\in \{ 1,\ldots ,8760\}$, then this information is unused in problem (1). This leads to an incorrect estimate of the PDF.

These disadvantages may be corrected with convex combinations of Weibull PDFs, which are discussed in the next section.

Fig. 2

Illustration of the monthly mean wind speeds for the Caribbean sites 1 (red curve) and 2 (blue curve)

Table 1 Characteristics of three sites

For sites in the Caribbean, it is well known (compare Wang 2007) that mean wind speeds have two local maxima in summer and winter and two local minima in fall and spring. Figure 2 illustrates monthly mean wind speeds for sites 1 and 2 given in Table 1. Based on these observations, it certainly makes sense to incorporate monthly distributions into an approach with convexly combined Weibull PDFs. This leads to an significant improvement of the PDF for difficult sites.

Convex combinations of Weibull PDFs

Taking the temporal order of the data into account, one can consider (seasonal or) monthly wind speeds. Consequently, for every month of a year the Weibull PDF is estimated with the maximum likelihood method. Then, we consider a convex combination of these 12 Weibull PDFs, i.e. we formulate a PDF $\bar{f}_{\bar{\lambda },\bar{k},\bar{c}}$ for the whole year with

$$\begin{aligned} {\bar{f}_{\bar{\lambda },\bar{k},\bar{c}}(v)}\\& {}:=\, \sum \limits _{j=1}^{12} \lambda _j f_{k_j, c_j}(v)\\& {}= \left\{ \begin{array}{ll} 0 &{}\quad \text {if } v < 0\\ \sum \limits _{j=1}^{12} \lambda _j \frac{k_j}{c_j} \left( \frac{v}{c_j}\right) ^{k_j-1}e^{-(\frac{v}{c_j})^{k_j}} &{}\quad\text {if } v \ge 0. \end{array}\right. \end{aligned}$$

Fig. 3

Example of a convexly combined PDF (3) with $\tilde{\lambda }=(0.1,\,0.5,\, 0.4,\, 0,\, \ldots ,\, 0)$, $\bar{k} = (2,\, 4.5,\, 1,\, \ldots ,\, 1)$ and $\tilde{c}=(5,\, 2.5,\, 8,\, 1,\, \ldots ,\, 1)$

Here we have $\bar{\lambda } := (\lambda _1, \dots , \lambda _{12})$, $\bar{k} := (k_1, \dots , k_{12})$, $\bar{c} := (c_1, \dots , c_{12})$ with $k_j, c_j >0$, $\lambda _j \in [0,1]$ for all $j\in \{ 1,\ldots ,12\}$ and $\sum\nolimits_{j=1}^{12} \lambda _j =1$. The coefficients $\lambda _1,\ldots ,\lambda _{12}$ can be chosen as quotient of the number of days per considered month and the number of days per year. Since wind speeds with 0 m/s are possible, we consider an exponential PDF as a special Weibull PDF $f_{1, c_0}$ with $k_0=1$, $c_0 >0$ and

$$\begin{aligned} f_{1, c_0}(v) = \left\{ \begin{array}{ll} 0 & \quad \text {if } v < 0\\ \frac{1}{c_0}e^{-\frac{v}{c_0}} & \quad \text {if } v \ge 0. \end{array}\right. \end{aligned}$$

This special Weibull PDF is then added to the convex combination of the 12 PDFs so that we investigate the new convex combination $\tilde{f}_{ \tilde{\lambda },\bar{k},\tilde{c} }$ with $\tilde{\lambda } := (\lambda _0, \bar{\lambda })$, $\tilde{c} := (c_0, \bar{c})$ where $c_0>0$, $\lambda _j \in [0,1]$ for all $j\in \{0,\ldots ,12\}$, and $\sum\nolimits_{j=0}^{12} \lambda _j =1$. This new convex combination is then given by

$$\begin{aligned} {\tilde{f}_{\tilde{\lambda },\bar{k},\tilde{c}}(v)}\nonumber \\:=\, & {} \sum \limits _{j=0}^{12} \lambda _j f_{k_j, c_j}(v)\nonumber \\=\, & {} \left\{ \begin{array}{ll} 0 &{}\quad\text {if } v < 0\\ \sum \limits _{j=0}^{12} \lambda _j \frac{k_j}{c_j}\left( \frac{v}{c_j}\right) ^{k_j-1}e^{-(\frac{v}{c_j})^{k_j}} &{}\quad \text {if } v \ge 0. \end{array}\right. \end{aligned}$$

(3)

An example of such a convex combination is illustrated in Fig.  3.

Fig. 4

Illustration of the objective function (2) with the variables k and c with wind speeds at site 1 and the optimal solution at $k=1.46$ and $c=4.04$

If we apply the maximum likelihood method to the convex combination (3) of Weibull PDFs with the logarithmic simplification to Eq. (2), we get the following nonlinear optimization problem

$$\begin{aligned} \begin{array}{l} \displaystyle \max \quad \sum _{i=1}^{8760} \ln \sum _{j=0}^{12} \lambda _j \frac{k_j}{c_j}\Bigl ( \frac{v_i}{c_j}\Bigr )^{k_j-1} e^{-(\frac{v_i}{c_j})^{k_j}} \\ \text {subject to the constraints}\\ \lambda _j \ge 0 \ \forall \ j\in \{ 0,\ldots ,12\}\\ k_j \ge \varepsilon \ \forall \ j\in \{ 1,\ldots ,12\}\\ c_j \ge \delta \ \forall \ j\in \{ 0,\ldots ,12\}\\ \displaystyle \sum _{j=0}^{12} \lambda _j =1\\ \lambda _0 = c_0h_0\\ (\tilde{\lambda },\bar{k}, \tilde{c})\in \mathbb {R}^{13}\times \mathbb {R}^{12}\times \mathbb {R}^{13}, \end{array} \end{aligned}$$

(4)

where $\epsilon , \delta >0$ are given lower bounds, $k_{0}:=1$ and $h_0$ equals the relative frequency of the wind speeds with 0 m/s. In problem (4), the following adaptations are already modelled:

• 1. The original objective function appears in a logarithmic form.

• 2. All observed wind speeds are taken into account including wind speeds with 0 m/s.

• 3. The positivity of the parameters $k_j$ ($j\in \{ 1,\ldots , 12\}$) and $c_j$ ($j\in \{ 0,\ldots ,12\}$) is ensured by the lower bounds $\varepsilon$ and $\delta$.

• 4. The last constraint ensures the right PDF value at 0 m/s.

The optimization problem (4) is a constrained problem with a highly nonlinear objective function. In general, methods of continuous optimization determine at most local optima. Figure 4 illustrates the graph of the logarithmic objective function (2) for the classical maximum likelihood method using wind speeds at Jamaica. This figure already highlights the complexity of this problem.

Fig. 5

Optimized PDF of the wind speeds at site 1 (discretization of the histogram by 0.3 m/s)

Procedure

Based on the remarks of the previous section, we now present a procedure for the optimization of the PDF of wind speeds.

Fig. 6

Optimized PDF of the wind speeds at site 1 (discretization of the histogram by 0.05 m/s)

Instead of the SQP method, one can also choose any numerical method of smooth constrained optimization. Since the objective function in problem (4) is highly nonlinear, one cannot expect that an SQP method finds the global solution of this problem. It is known that the computed solution strongly depends on the choice of the starting point. Therefore, the SQP method, which is not a method of global optimization, is repeatedly applied to different starting points. Among all computed points, one then selects this one with largest objective function value. This leads to more realistic numerical results.

Numerical results

The algorithm in the previous section is now applied to the wind speeds at site 1. At this site, we have $h_0=0$, i.e. there are no wind speeds with 0 m/s.

A first investigation uses the special starting vector with the calculated monthly parameters $\lambda _j$, $k_j$ and $c_j$ ($j\in \{ 1,\ldots ,12\}$) and sets $\ell _{\text {max}}=1$, i.e. the optimization problem (4) is only solved with this special starting vector. The data of this starting vector are given in the columns $\tilde{\lambda }_{\text {start}}$, $\bar{k}_{\text {start}}$ and $\tilde{c}_{\text {start}}$ in Table 2. The parameters of the exponential PDF are chosen as $\lambda _0:=0$ and $c_0:=1$. It is interesting to note that the $k_{j}$ and $c_{j}$ Weibull parameters ($j\in \{ 1,\ldots ,12\}$) vary significantly among the months. This shows that a convex combination of Weibull PDFs certainly makes sense.

Table 2 Starting vector and solution vector for $\ell _{\text {max}}=1$

Table 3 Solution vector for $\ell _{\text {max}}=45,000$

The constrained optimization problem (4) is solved by the SQP method of the optimization toolbox of MATLAB. The components of the obtained solution vector can be found in the columns $\tilde{\lambda }_{\text {opt}}$, $\bar{k}_{\text {opt}}$ and $\tilde{c}_{\text {opt}}$ of Table 2. The parameters of the exponential PDF are unchanged. It is evident from the data in Table 2 that the components of the starting vector are quite different from the components of the solution vector. This optimization leads to an improvement of the value of the objective function by 7.08% in comparison with the objective function value at the starting vector.

In a second investigation, the algorithm is used as given in the previous section. Now the parameter $\ell _{\text {max}}=45,000$ is chosen, i.e. 45,000 constrained optimization problems are to solve. An average CPU time for the execution of the SQP method is 213 s on an 8 core processor workstation. Table 3 presents the solution vector. The objective function value at this solution is improved by 7.38% in comparison with the objective function value at the starting vector given in Table 2. Figures 5 and 6 illustrate the optimized PDF with different discretization. In Figure 5, the histogram of wind speeds is discretized by 0.3 m/s, whereas the discretization of 0.05 m/s is chosen in Fig. 6. The finer discretization in Fig. 6 makes clear why a standard Weibull approach may lead to unacceptable results at difficult sites like the one at Jamaica. For the optimized PDF, one can easily determine the cumulative distribution function (CDF) illustrated in Fig. 7.

Fig. 7

Optimized CDF of the wind speeds at site 1.

Kolmogorov–Smirnov test

In the previous sections, we have concentrated ourselves to a good type of approximation of the CDF of wind speeds at a specific site. But now we test the hypothesis that the wind speed as random variable has the optimized CDF obtained by the presented algorithm. One accepts this hypothesis, if the optimized CDF and the empirical CDF are in a certain sense close together. The well-known Kolmogorov–Smirnov (KS) test (e.g. see D’Agostino and Stephens 1986) can be used for the test of this hypothesis.

For the Kolmogorov–Smirnov test, the wind speeds (8760 numbers) at site 1 are randomly splitted into two data sets with 4380 numbers. The first data set is used for the application of the algorithm of this paper. This leads to an optimized CDF, which is then compared with the empirical CDF of the second data set. Then, the Kolmogorov–Smirnov test is applied to these two CDFs. We get the result that with a level of significance of 5% the hypothesis is accepted that the optimized CDF is the true CDF of the second sample. In fact, the calculated test statistic value 0.017 is less than the critical value 0.021 of the Kolmogorov–Smirnov test. This shows that the approach of this paper is suitable for a good determination of the CDF of wind speeds.

If one works with the whole data set of 8760 wind speeds per year, the critical value in the Kolmogorov–Smirnov test at a level of significance of 5% is given by $1.358/\sqrt{8,760}\approx 0.0145$, i.e. for the supremum of deviations below this value, the hypothesis is accepted that the calculated CDF is the true CDF of the wind speeds as random variable. Assuming the correctness of the hypothetical CDF, there is a maximum probability of 5% observing test statistic values above the critical value, thus rejecting the hypothesis falsely.

If one considers only the classical Weibull CDF and a site with a higher number of hours with wind speeds 0 m/s, then the classical CDF F and the empirical CDF $\hat{F}$ have the derivatives $F'(0)=0$ and $\hat{F}'(0)=h_{0}$ (given in the algorithm). So, the expression $\sup_{v\ge 0}| F(v)-\hat{F}(v)|$ may be greater than the critical value 0.0145 so that the tested hypothesis is rejected. The convex combination presented in this paper tries to avoid this disadvantage.

Fig. 8

Histogram and estimated PDFs of wind speeds at site 1: Weibull PDF (blue curve), bimodal Weibull PDF (green curve), Gamma–Weibull PDF (black curve) and convexly combined Weibull PDF (red curve)

Table 4 Numerical results for different sites and various approaches

Numerical comparisons

The presented new method is now compared with other approaches for an estimation of the CDF of wind speeds for different wind sites. Table 1 gives some characteristics of three sites (see https://mesonet.agron.iastate.edu/request/download.phtml?network=CU__ASOS, https://mesonet.agron.iastate.edu/request/download.phtml?network=JM__ASOS and https://mesonet.agron.iastate.edu/request/download.phtml?network=JP__ASOS).

Fig. 9

Histogram and estimated PDFs of wind speeds at site 2: Weibull PDF (blue curve), bimodal Weibull PDF (green curve), Gamma–Weibull PDF (black curve) and convexly combined Weibull PDF (red curve)

For these sites, the PDF of wind speeds is calculated for various standard approaches. First of all, the (standard) Weibull PDF is determined for the three sites. Moreover, the bimodal Weibull PDF also known as Weibull–Weibull PDF and the mixture Gamma–Weibull PDF are calculated with the wind speed data. Figures 8, 9 and 10 illustrate the histograms of the measured wind speeds together with the PDFs obtained with the standard Weibull approach, the bimodal Weibull method, the Gamma–Weibull approach and the new method of this paper with $\ell _{\text {max}}=1$. All numerical results are listed in Table 4.

Figures 8, 9 and 10 and Table 4 show that there are significant differences between the computed PDFs. It is obvious that the standard Weibull approach is not suitable for difficult international sites.

Furthermore, the new method of this paper seems to be superior in contrast to the other methods. These discrepancies between the PDFs of the considered approaches are certainly smaller, if one investigates wind sites with a more uniform PDF.

The Kolmogorov–Smirnov test is carried out for all four approaches and all three sites. For every site, the KS test statistic value of the convexly combined Weibull PDF is the smallest among all used methods, which means that the new method determines the best approximation of the CDF. But this better performance of the new approach is reached by a higher numerical effort.

Conclusion

This paper modifies the classical Weibull PDF for wind speeds using a convex combination of Weibull PDFs. Optimal parameters can be obtained with the maximum likelihood estimation as an optimal solution of a highly nonlinear constrained optimization problem. By a monthly splitting of wind data, one gets for site 1 at Jamaica an improvement of more than 7% of the objective function, and with a level of significance of 5%, we can accept the hypothesis that the optimized CDF is the true CDF of wind speeds. With such an optimized CDF, we are able to investigate and analyse wind speeds more precisely than with the known techniques as shown by numerical comparisons.

Fig. 10

Histogram and estimated PDFs of wind speeds at site 3: Weibull PDF (blue curve), bimodal Weibull PDF (green curve), Gamma–Weibull PDF (black curve) and convexly combined Weibull PDF (red curve)