An Investigation of Wind Direction and Speed in a Featured Wind Farm Using Joint Probability Distribution Methods

Wind direction and speed are both crucial factors for wind farm layout; however, the relationship between the two factors has not been well addressed. To optimize wind farm layout, this study aims to statistically explore wind speed characteristics under different wind directions and wind direction characteristics. For this purpose, the angular–linear model for approximating wind direction and speed characteristics were adopted and constructed with specified marginal distributions. Specifically, Weibull–Weibull distribution, lognormal–lognormal distribution and Weibull–lognormal distribution were applied to represent the marginal distribution of wind speed. Moreover, the finite mixture of von Mises function (FVMF) model was used to investigate the marginal distribution of wind direction. The parameters of those models were estimated by the expectation–maximum method. The optimal model was obtained by comparing the coefficient of determination value (R2) and Akaike’s information criteria (AIC). In the numerical study, wind data measured at a featured wind farm in north China was adopted. Results showed that the proposed joint distribution function could accurately represent the actual wind data at different heights, with the coefficient of determination value (R2) of 0.99.


Introduction
With the gradually increasing consumption of nonrenewable energy, such as oil, coal, and natural gas, renewable energy has become the hope of humankind to solve the current energy shortage and future energy crisis [1,2].Wind energy is a well-known renewable energy source and has been used worldwide due to its advantages of clean, renewable, and large reserves [3][4][5].People capture wind energy through wind turbines, and wind speed and wind direction characteristics are therefore crucial for calculating wind power and the design and arrangement of wind turbines [6].The wind speed probability density distribution can be used to estimate wind power density [7,8] and can also be applied to determine the most dominant wind direction [9][10][11].In order to optimize wind farm layout in complex terrain, the joint probability distribution of wind direction and wind speed have gained great attention [12][13][14] as it can be used to obtain wind speed characteristics under different wind directions in complex terrain.Wind speed characteristics under different wind directions play a significant role in arrangement of wind turbines.It should be noted that wind speed with the

Joint Distribution
In this paper, we used the angular-linear model proposed by Johnson and Wehrly [12] to construct the joint probability distribution of wind direction and speed.The joint probability distribution of wind direction and speed is as follows (Equations ( 1) and ( 2)): where f V (v) is the probability distribution of the wind speed; f Θ (θ) is the probability distribution of the wind direction; g(•) is the probability distribution of the angular variable ζ.The angular variable ζ is as follows: where F V (v) is the cumulative probability distribution of the wind speed; F Θ (θ) is the cumulative probability distribution of the wind direction.

Wind Speed Distribution
For the probability distribution of wind speed, we used W-W, given by Equation ( 3), L-L, given by Equation ( 4), and W-L, given by Equation ( 5): where w m are weight coefficients that add to one; the parameter α m is a shape parameter; β m is a scale parameter with the same units as the wind speed.
where w m are weight coefficients that add to one; the parameter µ m is the mean of the associated normal distribution; σ m is standard deviation of the associated normal distribution.
where w 1 and w 2 are weight coefficients that add to one.

Circular Variable Distribution
For the probability distribution of wind direction, we used a finite mixture of von Mises function.The FVMF is a mixture distribution that consists of Mth single von Mises distribution.The finite mixture von Mises function is as follows (Equations ( 6) and ( 7)): where w m are weight coefficients that add to one; the parameter µ m is the mean wind direction; and the parameter k m is the concentration parameter.Here, I 0 (k m ) is the modified Bessel function of the first kind and order zero [33].I 0 (k m ) is as follows: Sustainability 2018, 10, 4338 For the probability distribution of the angular variable ζ, the two component mixtures von Mises function was used [6].

Cumulative Distribution Function
The cumulative probability distribution of the wind speed can be derived through Equation (8).The cumulative probability distribution of the wind direction can be derived through Equation (9).

Parameter Estimation
In this paper, we chose the expectation-maximization method to estimate the parameters of W-W, L-L, W-L, and FVMF [34].

Parameter Estimation of W-W Distribution
If we let b m = β m and c m = α m β m , the W-W in Equation (3) can also be written as follows (Equation ( 10)): The W-W parameters (w m , b m , and c m ) can be estimated through the actual wind speed data.The W-W parameters can be derived by Equations ( 11)-( 14): where N is the number of measured wind speed data; p n = p(m|v n ; w m , µ m , σ m ) and p n is the distribution of hidden variables in the expectation-maximization method [11].

Parameter Estimation of L-L Distribution
The L-L parameters (w m , µ m , and σ m ) can be determined through the actual wind speed data.The L-L parameters can be derived by Equations ( 15)-( 18): Sustainability 2018, 10, 4338 where N is the number of measured wind speed data; p n = p( m|v n ; w m , µ m , σ m ) and p n is the distribution of hidden variables in the expectation-maximization method [11].

Parameter Estimation of W-L Distribution
If we let b = β 1 and c = α 1 β 1 , the W-W in Equation ( 5) can also be written as follows (Equation ( 19)): The W-L parameters (w 1 , w 2 , b, c, µ, and σ) can be determined through the actual wind speed data.The W-L parameters can be derived by Equations ( 20)- (26): where N is the number of measured wind speed data; Θ = {w 1 , w 2 , b, c, µ, σ} and p( m|v n ; Θ) is the distribution of hidden variables in the expectation-maximization method [11].
The FVMF parameters {w m , sin(µ m ), cos(µ m ), k m } can be determined through the actual wind direction data [25].The FVMF parameters can be derived by Equations ( 28)-(32): where N is the number of measured wind direction data; p n is the distribution of hidden variables in the expectation-maximization method [11]; I 1 ( km ) is the modified Bessel function of the first kind and order one.Here, the derivative of I 0 (k m ) with respect to k m is I 1 (k m ).Besides, to evaluate the goodness of fit for the theoretical distribution [33], the coefficient of determination (R 2 ) and Akaike's information criteria (AIC) are given by Equations ( 33) and ( 34): where Fn are estimated cumulative probabilities; F = ∑ N n=1 Fn /N; F n are empirical cumulative probabilities.
where L is the likelihood; k is the number of parameter in the fitted model.

Wind Power Density
Wind power density can be used to evaluate the potential of wind energy in a region.The mean wind power density [35][36][37] can be expressed as follows (Equation (35)): where E represent the mean wind power density (W/m 2 ); N is the number of records in the set period; v i is the wind speed of the ith record (m/s); ρ is the density of air (kg/m 3 ).In this paper, the value of ρ was 1.225 kg/m 3 .The probability distribution of wind power density at different wind speed can be determined using Equation (36): In addition, the wind power density can be estimated by Equation ( 37): The joint probability distribution of wind direction and speed can be used to evaluate the wind speed distribution at different wind directions, given by Equation (38): Moreover, two useful wind speeds parameters, namely the most probable wind speed and the wind speed carrying maximum energy, can be derived by the wind speed distribution and wind power density distribution [38,39].The most probable wind speed is the wind speed of maximum point in the wind speed distribution.The wind speed carrying maximum energy is the wind speed of maximum point in the wind power density distribution [40].

Results and Discussion
Figure 1 was drawn by analyzing the data at a height of 70 m, and it demonstrates the wind speed distributions (Figure 1a) and the probability distribution of wind power density (Figure 1b).The parameters in the wind speed distributions, the corresponding R 2 coefficient, and AIC are shown in Table 1.It can be seen from Figure 1a that W-L distribution was better than W-W and L-L distribution for the histogram of wind speed data, which is also clear from in Table 1.Therefore, we chose W-L distribution as the marginal distribution of wind speed to construct the joint distribution.Figure 1b illustrates the probability density distribution of wind power density.The estimated wind power density was 270.79 W/m 2 , 307.34 W/m 2 , and 271.97 W/m 2 corresponding to W-W, L-L, and W-L, respectively.The mean wind power density was 270.53 W/m 2 .Figure 2 shows that the wind direction rose at a height of 70 m (Figure 2a) and the wind energy rose at a height of 70 m (Figure 2b).As shown in Figure 2, wind direction 1, 8, and 13 (N, SSE, and W) had a large wind direction frequency of 9.28%, 7.66%, and 10.41%, respectively, while wind direction 12 (WSW) had a large wind power density frequency (Figure 2b) of 14.92%.Therefore, in order to investigate their internal correspondence, we needed to study the wind speed characteristics under different wind directions and investigate their characteristics.Figure 3 was obtained by numerical calculation of the circular variable, which was measured at a height 70 m.Figure 3a shows the wind direction distributions, and Figure 3b depicts the probability distribution of the angular variable ζ.It can be seen from the histogram of the measured wind direction data in Figure 3a that there were up to four main wind directions; therefore, the FVMF with M = 2, M = 3, and M = 4 was chosen to fit the wind direction distribution.Figure 3b illustrates that there were two distinct peaks in the histogram of the angular variable ζ; Therefore, the FVMF with Figure 3 was obtained by numerical calculation of the circular variable, which was measured at a height 70 m.Figure 3a shows the wind direction distributions, and Figure 3b depicts the probability distribution of the angular variable ζ.It can be seen from the histogram of the measured wind direction data in Figure 3a that there were up to four main wind directions; therefore, the FVMF with M = 2, M = 3, and M = 4 was chosen to fit the wind direction distribution.Figure 3b illustrates that there were two distinct peaks in the histogram of the angular variable ζ; Therefore, the FVMF with M = 2 was selected to construct the probability distribution of the angular variable ζ.The parameters in the wind direction distributions, the corresponding R 2 coefficient, and AIC are shown in Table 2.The parameters in the angular variable ζ distribution function and the corresponding R 2 coefficient are shown in Table 3. From Figure 3a and Table 2, it can be seen that the FVMF with M = 4 was the best fitted model for the wind direction data.In addition, in order to estimate the parameters of the FVMF model, we redefined sin(µ m ) and cos(µ m ) instead of µ m .Here, the µ m should be calculated by sin(µ m ) and cos(µ m ) to describe the mean wind direction.The FVMF with M = 4 reflected that the most dominant wind direction was 4.69 rad (268.6 • ).As displayed in Figure 3b and Table 3, the FVMF with M = 2 was suitable to describe the angular variable ζ distribution.3. From Figure 3a and Table 2, it can be seen that the FVMF with M = 4 was the best fitted model for the wind direction data.In addition, in order to estimate the parameters of the FVMF model, we redefined sin ( ) and cos ( ) instead of .Here, the should be calculated by sin ( ) and cos ( ) to describe the mean wind direction.The FVMF with M = 4 reflected that the most dominant wind direction was 4.69 rad (268.6°).As displayed in Figure 3b and Table 3, the FVMF with M = 2 was suitable to describe the angular variable ζ distribution.Figure 4 shows the joint probability distribution of actual wind direction and speed at a height of 70 m.Figure 5 shows the joint probability distribution at a height of 70 m, which was obtained with the angular-linear distribution model.According to the proposed joint probability distribution at a height 70 m, the wind speed distribution under 16 wind directions could be determined.The corresponding R 2 value can be seen in Table 4.As can be seen from Figures 4 and 5 and Table 4, the proposed model had a great degree of fit to the sample data at a height 70 m.Moreover, in order to verify the validity of the selected model, the proposed model was used to construct the joint probability distribution at a height of 50 m.It should be borne in mind that the wind direction varied slightly from 70 m to 50 m.Figure 4 shows the joint probability distribution of actual wind direction and speed at a height of 70 m.Figure 5 shows the joint probability distribution at a height of 70 m, which was obtained with the angular-linear distribution model.According to the proposed joint probability distribution at a height 70 m, the wind speed distribution under 16 wind directions could be determined.The corresponding R 2 value can be seen in Table 4.As can be seen from Figures 4 and 5 and Table 4, the proposed model had a great degree of fit to the sample data at a height 70 m.Moreover, in order to verify the validity of the selected model, the proposed model was used to construct the joint probability distribution at a height of 50 m.It should be borne in mind that the wind direction varied slightly from 70 m to 50 m.Figure 6 shows the joint probability distribution at a height of 50 m.In the proposed joint distribution model, W-L distribution was selected to process the wind speed data, and the FVMF with M = 4 was the fitted model for the wind direction data.Moreover, the FVMF with M = 2 was chosen to describe the angular variable ζ distribution.The parameters in the joint probability distribution at a height of 50 m and the corresponding R 2 coefficient are shown in Table 5.According to the proposed joint probability distribution at a height of 50 m, the wind speed distribution under 16 wind directions could be determined.The corresponding R 2 value can be seen in Table 6.It can be seen from the above analysis that the proposed model could well present wind characteristics at different heights.Figure 6 shows the joint probability distribution at a height of 50 m.In the proposed joint distribution model, W-L distribution was selected to process the wind speed data, and the FVMF with M = 4 was the fitted model for the wind direction data.Moreover, the FVMF with M = 2 was chosen to describe the angular variable ζ distribution.The parameters in the joint probability distribution at a height of 50 m and the corresponding R 2 coefficient are shown in Table 5.According to the proposed joint probability distribution at a height of 50 m, the wind speed distribution under 16 wind directions could be determined.The corresponding R 2 value can be seen in Table 6.It can be seen from the above analysis that the proposed model could well present wind characteristics at different heights.The wind direction was recorded in 16 intercardinal directions.Each wind direction of the 16 intercardinal directions had a corresponding wind speed and a wind power density distribution.Due to the limitation of space, we have only given the results of two wind directions-wind direction 12 (E) and wind direction 13 (W)-in Figure 7.The figure was drawn by applying the joint probability density distribution at a height of 70 m.As can be seen from Figure 7, the most probable wind speeds for wind direction 12 and 13 were 6.5 m/s and 6.95 m/s, respectively; the wind speeds carrying maximum energy were 9.34 m/s and 9.63 m/s, respectively.
Sustainability 2018, 10, x; doi: FOR PEER REVIEW 11 of 12 The wind direction was recorded in 16 intercardinal directions.Each wind direction of the 16 intercardinal directions had a corresponding wind speed and a wind power density distribution.Due to the limitation of space, we have only given the results of two wind directions-wind direction 12 (E) and wind direction 13 (W)-in Figure 7.The figure was drawn by applying the joint probability density distribution at a height of 70 m.As can be seen from Figure 7, the most probable wind speeds for wind direction 12 and 13 were 6.5 m/s and 6.95 m/s, respectively; the wind speeds carrying maximum energy were 9.34 m/s and 9.63 m/s, respectively.
The FVMF with M = 4 was a good fitted model for the wind direction data and could be used to investigate the wind direction characteristics for four seasons.The four seasons were defined as follows: (a) spring: March, April, and May; (b) summer: June, July, and August; (c) Autumn: September, October, and November; (d) winter: December, January, and February.Figure 8 depicts the wind direction distributions over four seasons.It was drawn using the wind direction data of 70 m for three years.The parameters in the wind direction model over four seasons and the corresponding R 2 coefficient are shown in Table 7.As can be seen from Figure 8, the most dominant wind direction was 5 rad (286.6°) in the spring season, 6 rad (348.8°) in the summer season, 4.5 rad (260.2°) in the autumn season, and 4.7 rad (267.1°) in the winter season.There were more than two prevailing directions in the four seasons.It was concluded from Figure 8 that the wind direction had seasonal features.The FVMF with M = 4 was a good fitted model for the wind direction data and could be used to investigate the wind direction characteristics for four seasons.The four seasons were defined as follows: (a) spring: March, April, and May; (b) summer: June, July, and August; (c) Autumn: September, October, and November; (d) winter: December, January, and February.Figure 8 depicts the wind direction distributions over four seasons.It was drawn using the wind direction data of 70 m for three years.The parameters in the wind direction model over four seasons and the corresponding R 2 coefficient are shown in Table 7.As can be seen from Figure 8, the most dominant wind direction was 5 rad (286.

12 Figure 1 .
Figure 1.The distribution of wind speed characteristics at a height of 70 m: (a) wind speed distributions, (b) probability distribution of wind power density.

Figure 1 .
Figure 1.The distribution of wind speed characteristics at a height of 70 m: (a) wind speed distributions, (b) probability distribution of wind power density.

Figure 2 .
Figure 2. The rose diagrams at a height of 70 m: (a) wind direction rose; (b) wind energy rose.

Figure 2 .
Figure 2. The rose diagrams at a height of 70 m: (a) wind direction rose; (b) wind energy rose.
Sustainability 2018, 10, x; doi: FOR PEER REVIEW 8 of 12 M = 2 was selected to construct the probability distribution of the angular variable ζ.The parameters in the wind direction distributions, the corresponding R 2 coefficient, and AIC are shown in Table 2.The parameters in the angular variable ζ distribution function and the corresponding R 2 coefficient are shown in Table

Figure 3 .
Figure 3.The fitting of the finite mixture of von Mises function for the circular variable: (a) wind direction distributions; (b) probability distribution of the angular variable ζ.

Figure 3 .
Figure 3.The fitting of the finite mixture of von Mises function for the circular variable: (a) wind direction distributions; (b) probability distribution of the angular variable ζ.

Figure 4 .
Figure 4.The joint probability distribution of actual wind direction and speed at a height of 70 m.

Figure 4 .Figure 4 .
Figure 4.The joint probability distribution of actual wind direction and speed at a height of 70 m.

Figure 5 . 20 WFigure 5 .
Figure 5.The joint probability distribution at a height of 70 m.

Figure 6 .
Figure 6.The joint probability distribution at a height of 50 m.

Figure 6 .
Figure 6.The joint probability distribution at a height of 50 m.

Figure 7 .
Figure 7.The distribution of wind speed characteristics at a height 70 m: (a) wind speed distributions under different wind directions; (b) probability distributions of wind power density under different wind directions.

Figure 7 .
Figure 7.The distribution of wind speed characteristics at a height 70 m: (a) wind speed distributions under different wind directions; (b) probability distributions of wind power density under different wind directions.

6 •
) in the spring season, 6 rad (348.8 • ) in the summer season, 4.5 rad (260.2 • ) in the autumn season, and 4.7 rad (267.1 • ) in the winter season.There were more than two prevailing directions in the four seasons.It was concluded from Figure 8 that the wind direction had seasonal features.Sustainability 2018, 10, x; doi: FOR PEER REVIEW 12 of 12

Figure 8 .
Figure 8.The wind direction distributions over four seasons.

Figure 8 .
Figure 8.The wind direction distributions over four seasons.

Table 1 .
Fitted parameters in wind speed distributions.

Table 1 .
Fitted parameters in wind speed distributions.

Table 2 .
Fitted parameters in wind direction models.

Table 2 .
Fitted parameters in wind direction models.

Table 4 .
The R 2 value of the 70-m high wind speed distribution under 16 wind directions.

Table 4 .
The R 2 value of the 70-m high wind speed distribution under 16 wind directions.

Table 5 .
Fitted parameters in the joint probability distribution at a height of 50 m.

Table 5 .
Fitted parameters in the joint probability distribution at a height of 50 m.

Table 6 .
The R 2 value of the 50-m high wind speed distribution under 16 wind directions.

Table 7 .
Fitted parameters in wind direction models over four seasons.

Table 7 .
Fitted parameters in wind direction models over four seasons.