- Email: [email protected]
Wind farm layout optimization under uncertainty
Renewable Energy 107 (2017) 288e297
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
Wind farm layout optimization under uncertainty S.A. MirHassani*, A. Yarahmadi Department of Applied Mathematics, Faculty of Mathematics and Computer Science, Amirkabir University of Technology, Tehran, Iran
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 February 2016 Received in revised form 2 January 2017 Accepted 30 January 2017 Available online 1 February 2017
The optimal design of a wind farm with the aim of increasing the exploitation rate and power production is one of the challenging problems in the field of renewable energies. In this paper, the effect of using different hub height wind turbines in a farm on total power is studied. An exact mathematical formulation is presented in terms of the interaction matrix for multi-turbine wake effects considering different hub height wind turbines, and a new mixed-integer quadratic optimization model is developed. Then, the model is generalized to include the changes in wind characteristics and uncertain parameters. The model is solved using the linearization technique and an iterative method. The performance of the model and solution algorithm is evaluated by solving different problems borrowed from the literature. The computation results show the possibility of having a high-quality solution in a reasonable time in almost all cases. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Wind farm layout optimization MIP Wake effect Interaction matrix Uncertainty
1. Introduction Renewable energy is generally defined as energy that comes from natural resources such as sunlight, wind, rain, and waves. This kind of energy is endorsed more than fossil fuel sources for its sustainability and compatibility with climate change and global warming. Wind energy which is used in power production is one of the oldest types of renewable energy. Electrical power is generated by the wind flow in an appropriately designed wind farm. Therefore, determining the optimal location of wind turbines has a significant impact on the performance of the wind farm and the total produced energy. The optimization of the layout of a wind farm is achieved by considering two main factors: the expected power output and the wake effect. The wake effect should be minimized in order to maximize the power output. It has two important features: (i) a decrease in the wind speed; (ii) an increase in the level of turbulence of wind flow. We consider this problem as the wind farm layout optimization (WFLO) problem. 2. Literature review In the past decade, the wind turbine location problem was studied by many researchers. Mosetti et al. [1] solved the WFLO
* Corresponding author. E-mail addresses: [email protected] (S.A. MirHassani), aminyarahmadi@ aut.ac.ir (A. Yarahmadi). http://dx.doi.org/10.1016/j.renene.2017.01.063 0960-1481/© 2017 Elsevier Ltd. All rights reserved.
using a genetic algorithm to find the optimal number and location of wind turbines. They assumed a 10 10 grid of possible turbine locations used Jensen’s wake decay model to analyze the wake effect under various wind speeds and directions. Improving the genetic algorithm and considering the Mosetti objective function, Grady et al. [10] achieved better results. Marmidis et al. [2] used Mossetti and Grady’s objective function in a Monte Carlo simulation approach. Chowdhury et al. [3] determined the location of wind turbines considering the wind farm cost analysis based on the turbine rotor diameters and the number of turbines in the wind farm. They solved the wind farm model using the constrained Particle Swarm Optimization (PSO). Using a nested genetic algorithm, Chen et al. [4] optimized the layout of a 500m 500m wind farm with different hub height wind turbines. They proved that the objective value has improved slightly more than the case with identical hub wind turbines. They also considered different cost models applied to WFLO with different hub height wind turbines, and demonstrated that the result can improve the cost per unit power of a wind farm in this case. A few mathematical programming approaches have been discussed in the literature. Donovan [5] proposed multiple mixed-integer linear models based on vertex packing problem between a couple of wind turbines. Truner et al. [6] developed a new mathematical framework to optimize the layout of a 10 10 grid of possible turbine location in a wind farm. They introduced interaction matrices which modelled the effect of wake interactions between the turbines. They used quadratic integer programming as well as mixed-integer linear programming
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
approaches. In this paper, we propose a mathematical programming framework which is a precise tool for modeling the WFLO problem in order to optimize the layout of a wind farm with different hub height wind turbines. In addition, we formulate an exact wake effect structure and compare the result with the layout in Refs. [6,4,1]. The main factors considered in this study are: (i) the number of wind turbines, (ii) wind directions and speeds, (iii) hub height of wind turbines, (iv) power output, and (v) the exact wake effect interaction. A quadratic integer programming and mixed-integer linear program are developed for different hub height wind turbines. Finally we present a heuristic algorithm with realistic results to solve the proposed models.
3. Wake model When the wind hits a turbine, a cone of slower and more turbulent air is created behind it. This phenomenon, which is known as the wake effect in fluid-aerodynamics (see Fig. 1), has been investigated by researchers [7].
2
a¼
0:5 ln zz0
(3.3)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi a ¼ 0:5 1 1 CT
7 7 7 2 d 5
2a 1þ
aDij
(3.4)
As depicted in Fig. 2, it is preferred to place the turbines in an area with two main characteristics: (i) the highest wind speeds, (ii) outside the wake of other turbines. With the assumption that the terrain is flat, and considering the information in Section 3.1, we use Jensen’s wake decay model as follows: In this study, we assume that the total kinetic energy deficit due to the multiple wake is equal to the sum of individual kinetic energy deficits generated by each wake [7,8]. For more details on the wake model, see Jensen [8]. With this assumption, we can reform Equation (3.1) to calculate the wind speed at position j after multiple sequences of turbines:
2 vj ¼ v0 41
3 sffiffiffiffiffiffiffiffiffiffiffiffi X 2ffi dij 5 Where dij ¼ i2J
3
6 6 v ¼ v0 61 4
289
2a 1þ
aDdij
!2
(3.5)
aDdij þ rr
(3.1)
r1
3.1. The nomenclature
r1 ¼
aDdij
þ rr
(3.2)
All of the parameters involved in the following proposed models and wake model are defined in this section.
J ¼ 1; …; J L ¼ 1; …; L ℕ ¼ 1; …; N S ¼ 1; …; S
Set Set Set Set
K ¼ 1; …; K v0 vj
Set of types of wind turbine, k2K Velocity of free stream (m/s) Wind velocity at position j (m/s)
v r0 rr r1
Ddij
Downstream wind speed (m/s) Turbine radius ðmÞ Downstream rotor radius ðmÞ Wake radius ðmÞ Velocity loss in the wake at position j due to a turbine placed at position i Distance along the wind direction between turbines i and j ðmÞ
Dvij DEij
uij
of of of of
possible turbine locations i2J; j2J wind directions l2L different velocities of free stream n2ℕ scenarios, s2S
NT zk z0 pl pn ps ml a CT
Total number of wind turbines Hub height of wind turbine of type k ðmÞ Surface roughness length ðmÞ Probability of wind directions l Probability of velocity of free stream n Probability of scenario s Slope of line passing through the turbine for each direction
A AO P
Axial induction factor Trust coefficient Entrainment constant Swept area of wind blades ðm2 Þ Overlapping area between wake and a turbine ðm2 Þ Total power output ðkWÞ
Perpendicular distance along the wind direction between turbines i and j ðmÞ
ðXi ; Yj Þ
Coordinates of each turbine ðmÞ
Euclidean distance between turbines i and j ðmÞ
Wij
Interaction matrix
a
4. Interaction matrix calculation The exact structure of wake effect in a wind farm that determines which turbine is in the wake of the upstream turbine plays an important role in optimum wind turbine location. In addition, we need to take into account the wind direction in model formulation. In the present study, we propose a new method to define the exact mathematical structure of the wake effect in a wind farm by means of the interaction matrix concept. The dimensions of interaction matrices are related to the size of grid ðJ JÞ. Therefore, each element of the matrix shows the interaction between itself and J 1 other cells. We attribute a matrix to each direction of incoming wind using the following algorithm: Fig. 1. Schematic representation of the wake effect.
290
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
Fig. 2. The wake effect and a wind farm model.
Table 1 Classification of wind direction and turbines’ coordinates. p q < 3p 4 4 3p q < 5p
Xj < Xi
q < 74p
Xj > Xi
4 5p 4 7p 4
! ! 2 2 2 h2 þ r02 r12 2 1 h þ r1 r0 AO ¼ þ r1 cos 2hr0 2hr1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð h þ r0 þ r1 Þðh r0 þ r1 Þðh þ r0 r1 Þðh þ r0 þ r1 Þ 2 (4.5) r02 cos1
Yj > Yi
4
q < 94p
Yj < Yi
The parameter h indicates the distance between the center of wake effect and the center of turbine blades ðh ¼ Dvij Þ.
4.1. The algorithm i. First, we divide the terrain into four regions. Each region is determined with the direction of incoming wind and the turbine’s coordinates (Table 1). If both conditions in each row of Table 1 are satisfied, then a turbine at position j can be in the wake of a turbine at position i. ii. Formalize the passing line equation of each turbine depending on the wind direction as follows:
y ¼ ml ðx Xi Þ þ Yi
(4.1)
5. Framework with different hub height wind turbines Previous studies proved that using different hub height wind turbines in a wind farm can increase the wind farm’s output. In this section, we present an extended mathematical model for the WFLO problem. In the first part, the Log Law [9] is explained. The details of the generalized model are discussed in the second part. Finally, we appropriately modify the interaction matrices.
iii. According to Fig. 3, compute Dvij for the turbines which meet the first-step conditions, and then calculate Ddij for these turbines considering the Euclidean distance. iv. The value of uij is calculated as follows:
If
If
Dvij > r1 þ rr
uij ¼ 0
Then
r1 rr < Dvij < r1 þ rr
Then
(4.2)
2a
uij ¼ 1þ
aDdij
!2
AO A
aDdij þ rr (4.3)
If Dvij < r1 rr
Then
2a
uij ¼ 1þ
where:
aDdij
!2
(4.4)
aDdij þ rr Fig. 3. Determination of the wake radius at point.
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
5.1. The log law The log law is used to evaluate wind speed from the reference height zr to another level z (see Refs. [4,9]). In this study, the reference height belongs to the turbines which have higher hub heights. Based on the Log Law, one can realize that the wind speed goes up as much as the height increases:
ln zz0 vðzÞ ¼ vðzr Þ ln zr
(5.1)
z0
J K P P
maximize
k¼1 j¼1
291
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 31 u K J uX X ’2 k’ 5A ykj @v0 41 t dkk ij yi 0 2
(5.8)
k’ ¼1 i¼1
P PJ The term Kk¼1 j¼1 ykj v0 is equal to the constant value NT $v0 in which v0 represents the sum of free stream velocities before wind turbines. Therefore, we can replace the maximization model with the minimization one. Moreover considering
’
’
kk ukk ij ¼ dij
2
;
L 2 X ’ ’ ykj ¼ ykj and Wijkk ¼ pl uijkk l
(5.9)
l¼1
vðzr Þ is the free stream velocity before the wind turbines with higher hub heights.
the ultimate quadratic model is formulated as follows:
Minimize 5.2. Mathematical model
J P J P K P K P i¼1 j¼1 k¼1
’
k’ ¼1
’
Wijkk yki ykj
(5.10)
ykj 2Y As already mentioned, to compute the modified constant, we use Equation (3.3) to consider its changes related to turbines with different hub heights. The value of a will change for each turbine and is calculated by the following equation:
ak ¼
0:5 ln zz0k
(5.2)
2
3 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X X k2ffi vj ¼ v0 41 dij 5 where dkij ¼ k
i2J
2a 1þ
ak Ddij
!2 ck2K
ak Ddij þ rr (5.3)
In this case, Equation (3.5) results in the following form: Now, we introduce a new binary variable ykj equal to 1 if a turbine of type k is placed at position j, and to 0 otherwise. The model’s objective function is the generated power which is proportional to the cube of wind speed [6,10]. Thus, our proposed extended model is as follows:
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 313 u K J J uX X K P P ’2 k’ 5A ykj @v0 41 t dkk ij yi 0 2
maximize
k¼1 j¼1
Subject to
J P j¼1
K X
ykj 1
ykj ¼ N k
cj2J
(5.4)
k’ ¼1 i¼1
ck2K
(5.5)
(5.6)
Y shows the feasible region made by constraints (5.5) to (5.7). 5.3. Modified interaction matrix in case of different turbines ’
To evaluate ukk ij , it is essential to determine the type of the upstream turbine. In the case with identical turbines, an interaction matrix was developed for any situation with a different wind direction. In this case, we have jkj2 numbers of interaction matrices for each wind direction. As a result, the modified constant used in Equations (4.2)-(4.4) is specified by the type of turbine at position i, and the parameter h in Equation (4.5) will change as follows (see Fig. 4.):
h¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 Dvij þ ðzk zk’ Þ2
(5.11)
6. Methodology As mentioned in Section 4, we form an interaction matrix for each wind direction. For a given wind direction, the interaction matrix is a sparse and lower triangular matrix. Thus, as more wind directions are considered, the interaction matrix becomes dense. Consequently, solving the model with Formulation (5.10) is very difficult [11]. Two methods are proposed to solve this problem: (i) linearization technique, (ii) direct solution of the quadratic model [6]. We first discuss the linearization technique and generalize it for Formulation (5.10). Then, a heuristic method is introduced to solve the linearized model. 6.1. The linearization technique
k¼1
ykj 2f0; 1g;
cj2J
(5.7)
where the binary variable yki ¼ 1 if a turbine of type k is placed at position i2J The term N k refers to the number of turbines of type k: Thus, the constraint (5.5) requires that N k turbines be placed. The constraint (5.6) enforces that maximum of one turbine be placed at any position. This model is hard to solve due to the nonlinearity of the objective function, and is approximated by applying various transformations. Similarly, we can maximize the wind speed at each turbine instead of maximizing the total generated power [6]. Consequently, the objective function (5.4) will transfer to the following formulation:
minimize
J P J P K P K P i¼1 j¼1 k¼1
Subject to
J P j¼1
K P k¼1
ykj 1
k’ ¼1
ykj ¼ Nk
cj2J
’
’
Wijkk xkk ij
ck2K
(6.1)
(6.2)
(6.3)
292
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
Wake Turbine Overlap Wake
Turbine
Fig. 4. Schematic representation of a wake effect in a wind farm consisting of wind turbines with different hub heights.
’
k xkk ij yi ’
k xkk ij yj
’
’
ci; j2J; isj; ck; k’ 2K
(6.4)
ci; j2J; isj; ck; k’ 2K
(6.5)
’
k k xkk ij yi þ yj 1
ykj 2f0; 1g ’
ci; j2J; isj; ck; k’ 2K
cj2J; ck2K
ci; j2J; isj; ck; k’ 2K
xkk ij 0
(6.6) (6.7) (6.8)
Using this technique leads to the integer quadratic model’s change into a mixed-integer linear (MILP) one. More precisely, we can replace any term of the form yi yj with a new nonnegative variable xij with constraints xij yi , xij yj and xij yi þ yj 1 [12]. By generalizing the definition of the variable xij , we refor’ mulate (5.10) into a mixed-integer linear model where xkk ij ¼ 1 if ’ turbine types k and k are located in positions i and j; respectively, ’ and xkk ij ¼ 0 otherwise. Furthermore, in computational experiments, we add a constraint to accelerate the solution process which states if a turbine of type k and a turbine of type k’ are placed in positions i and j respectively, or j ; i. Therefore: ’
’
kk xkk ij ¼ xji
ci; j2J; isj; ck; k’ 2K
(6.9)
To solve this problem, we consider an iterative method to locate turbines gradually. This technique has two main advantages: (i) the algorithm runs very fast, (ii) the bounds generated by this method are more reliable than those produced by CPLEX through a suboptimal solution. In this technique, the model is first solved by a group (e.g. one, two, or more) of wind turbines, then their location will be fixed. Consequently, the model will locate more turbines in each iteration. The algorithm is terminated when the number of located turbines is equal to the total number of turbines in a wind farm. This procedure works for the MILP model with different hub height wind turbines as well as the case with identical wind turbines. In the first case, the algorithm is allowed to select the type and location of each wind turbine. Moreover, it is possible to determine the type of wind turbines before running the algorithm. To evaluate the efficiency of the iterative method, we considered 36 wind directions with similar probabilities, making the interaction matrix as dense as possible. Furthermore, we increased the size of grid and the number of turbines installed at each stage. The results are shown in Table 2. In addition, the number of turbines installed h i has been determined by 2J in order to challenge the MILP model for moderate values of NT . Table 2 shows the results obtained from the direct solution of the MILP model and iterative method. The first and second columns refer to the values of objective function and running time (s) of the direct method. We placed an upper limit of 24 h on the running time. The next two columns are similarly related to the iterative method, and the Gap column refers to the percentage of difference between these two methods. 7. Computational results
6.2. A heuristic method The computational complexity of the considered MILP models depends on the sparsity of interaction matrix ðWij Þ, size of grid ðJÞ; and the value of NT . We observed that the MILP problem is solvable for a small NT . For a moderate NT , it took too much time and delivered a low-quality solution. By taking into account more wind directions, the density of interaction matrix increases. Therefore, the procedure of solving the MILP problem will be longer and more complicated. As we mentioned, the size and computational complexity are more in the MILP problem with different hub height wind turbines than the case with identical wind turbines. Thus, solving the MILP problem directly is not an appropriate solution.
In this study, two types of wind turbines have been used, k2f1; 2g and, therefore, each wind direction was represented by four interaction matrices. The turbines characteristics and input parameters for WFLO cases are shown in Table 3 (for the cases with identical turbines) and Table 4 (for the cases with different hub height wind turbines). The three cases that have been discussed and reported in the literature have been selected to compare the layout and generated power, see Refs. [1,6,4]. In the first case (I), a dominated wind direction from north to south with the speed of 12 m=s is assumed. The second case (II) consists of 36 wind directions with equal probabilities of occurrence and the speed of 12m=s, and the third
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
293
Table 2 The performance comparison of direct solution and iterative method for linearized formulation (see Ref. [6]). J NT
Direct Method
Iterative Method
Gap
Obj. value
Sol. Time (s)
Obj. value
Sol. Time (s)
94 16 8 25 12 36 18
0:02373819 0:11450618 0:25562266 0:52404598
0:03 0:22 18:52 3278:72
0:023738198 0:114506187 0:256505218 0:524045990
0:00 0:83 2:67 7:25
49 24 64 32 81 40 100 50
0:87460352 1:54624036 2:21679568 3:38630941
23237:74 1005 86400 86400
0:880142860 1:538130101 2:233621949 3:393013736
25:3 55:33 137:35 206:84
0 0 0:0035 ^ 3:05 10ð9Þ 0:006 0:005 0:0076 0:002
Table 3 Parameters used in cases with identical wind turbines. r0 z z0 J Cell size Wind farm size CT P
20 m 60 m 0:3 100 10r0 ð200 m 200 mÞ 2 km 2 km 0:88 PJ 0:3v3i i
Table 4 Parameters used in the case with different hub height turbines. r0 z1 z2 z0 J Cell size Wind farm size CT P
20 m 50 m 78 m 0:3 100 10r0 ð200 m 200 mÞ 2 km 2 km 0:88 PJ;K 0:3vki i;k
Fig. 6. Subdivided wind farm.
subsection, we compare our layouts and generated power value obtained from the WFLO model with identical wind turbines solved by the iterative method, with results from Ref. [1] (MILP) and [6]. The second subsection consists of comparison of results for the case of different hub height wind turbines (6.1)-(6.8) solved by the iterative method and [4]. In the last subsection, we show that the use of wind turbines with different hub heights in a wind farm increases the output with respect to the same case with identical wind turbines. In Section 6.1, we mention that the velocity of incoming wind at the position of wind turbine with a larger hub height is more than that of the smaller one. In case (III), the values of wind speed are assumed to be 8 m=s, 12 m=s and 17 m=s. We suppose that these states represent the velocity of incoming wind at the position of turbines with 78 m of height. Using the Log Law, the velocities of incoming wind at the position of turbines with 50 m of height are 6:44 m=s, 11:04 m=s and 15:64 m=s; respectively. In addition, each state in case (III) is shown as an ordered couple. The first component has two parts which represent the incoming velocity at the
Fig. 5. Case (III), probabilities of occurrence of different wind directions and wind speeds.
one (III) is the case of 36 wind directions with unequal fractions of occurrence and multiple speeds of 8 m=s, 12 m=s and 17 m=s. The probabilities of occurrence for the third case are shown in Fig. 5.
Table 5 8.1.1 Comparison of results. NT
Model
Power (kW)
26
Mosetti et al. [1] Turner et al. [6] Present study Mosetti et al. [1] Turner et al. [6] Present study
9026 9087 9065 9835 9835 9777
8. Case study 30
In this part, we present our results and discuss the performance of our models and solution algorithm in different cases. In the first
294
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
Turner et al. (26 Turbines)
Present study (26 Turbines)
Turner et al. (30 Turbines)
Present study (30 Turbines)
Fig. 7. 8.1.1 layouts.
Table 6 8.1.2 Comparison of results.
Identical wind turbines
ðN 1 ; N 2 Þ
Model
Power (kW)
ð5; 20Þ
Chen et al. [4] Present study
8690 8923
Different hub height wind turbines
heights of 78 m and 50 m; respectively, and the second component demonstrates the direction of incoming wind. Thus, the objective function of MIQP minimization model in case (III) is as follows:
minimize
J P J P 2 P 2 P N P L P k’ ¼1
Subject to
k¼1 ykj 2Y
i¼1 j¼1 n¼1 l¼1
’
’
pn vn0 pl uijkk l yki ykj
(8.1)
50 m hub height wind turbine ’
’
s ¼ vn :ukk l where s corresponds By arranging ps ¼ pn :pl and ukk 0 ij ij to the pair ðn; lÞ, we will have NT L scenarios with the probability P S s of ps . It is clear that s¼1 p ¼ 1. Therefore, Formulation (8.1) is changed to the following MIQP formulation which we can easily transfer to the MILP model:
minimize
J P J P 2 P 2 P S P k¼1 k’ ¼1 i¼1 j¼1 s¼1
p
s
’ ’ uijkk s yki ykj
(8.2)
8.1. Case (I): constant wind speed and direction For this case, we consider a dominated wind direction from the
Chen et al. [4]
Present study
78 m hub height wind turbine Fig. 9. 8.1.3. layouts.
north with a constant wind speed of 12 m=s. The subdivided wind farm and our reference coordinates are shown in Fig. 6. The interaction matrix should be sparse, so the solution time and computational complexity are not too high. 8.1.1 In this section, we compare the WFLO MILP model with identical wind turbines in Ref. [6] solved by the iterative method with the results reported in Refs. [1,6] (see Table 5). In this case, due to the precise calculation of the wake effect (no estimation), there are some differences in the layouts and values of output power in comparison with the results in the literature (Fig. 7). 8.1.2 Table 6 compares the layouts and generated powers of our WFLO model (6.1)-(6.8) solved by the iterative method with the results reported in Ref. [4]. The model was quickly solved for N1ðz1 Þ ¼ 5; N 2ðz2 Þ ¼ 20 and the layout was as good as the one determined in Ref. [4], while the generated power was improved (Fig. 8). 8.1.3 In this part, we investigate the effect of using wind turbines with different hub heights in a wind farm that can improve power generation. For this purpose, we solved the WFLO MILP model in
Table 7 8.2.1 Comparison of results.
50 m hub height wind turbine 78 m hub height wind turbine Fig. 8. 8.1.2 layouts.
NT
Model
Power (kW)
19
Mosetti et al. [1] Turner et al. [6] Present study Turner et al. [6] Present study
6018 6081 6117 9040 9398
39
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
Turner et al. (19 Turbines)
Present study (19 Turbines)
Turner et al. (39 Turbines)
295
Present study (39 Turbines)
Fig. 10. 8.2.1. layouts.
8.2. Case (II): 36 wind directions and constant wind speed
Table 8 8.2.2 Comparison of results. ðN 1 ; N 2 Þ
Model
Power (kW)
ð7; 18Þ
Chen et al. [4] Present study
7189 7412
Chen et al. [4]
Present study
50 m hub height wind turbine
In this case, we assume that the incoming wind is multidirectional with the constant speed of 12 m=s. Each 36 angle has the same possibility of occurrence. The interaction matrices are as dense as possible and, therefore, the solution time and computational complexity are notable. 8.2.1 We fixed the optimal place of 19 and 39 turbines which were taken from Refs. [1,6] in order to calculate the output power in our proposed model (Table 7). The result shows that the output of the generated power is improved in both NT ¼ 19 and NT ¼ 39 cases. Solving the MILP model in Ref. [6] by the iterative method reached the optimal layouts more quickly than the direct method (Fig. 10). 8.2.2 For case (II), we considered the WFLO MILP model (6.1)(6.8) with N1ðz1 Þ ¼ 7; N2ðz2 Þ ¼ 18. Table 8 compares the power output obtained from our optimal layout with those of the layouts created by Chen et al. [4]. The results of the MILP model in this case lead to a layout which is as good as that of Chen et al. (Fig. 11). 8.2.3 For case (II), we considered the WFLO MILP model in Ref. [6] with NT ¼ 25 and z ¼ 78 m. The maximum power output is 7336 kW, while the maximum generated power was 7412 kW for the MILP model with different hub height wind turbines. Fig. 12 demonstrates the resulting layouts.
78 m hub height wind turbine Fig. 11. 8.2.2. layout.
Ref. [6] with NT ¼ 25 and z ¼ 78 m. The output power is 8916 kW while it was 8923 kW for the similar MILP model with different hub height wind turbines. Meanwhile, the layouts in the two cases are almost similar (Fig. 9).
Identical wind turbines
Different hub height wind turbines
8.3. Case (III): 36 wind directions and various wind speeds For this case, Fig. 5 shows 36 wind directions with three speeds of 8 m=s, 12 m=s and 17 m=s with unequal fractions of occurrence. In the case of MILP model with different hub height wind turbines, we suppose that these values of speed represent the velocity of incoming wind at the position of turbines with 78 m of height. Using the Log Law, the velocities of incoming wind at the position of turbines with 50 m of height are 6:44 m=s, 11:04 m=s and 15:64 m=s; respectively. The density of interaction matrices is the same as in case (II). 8.3.1 Table 9 presents the power output obtained from the optimal layouts of 15 and 39 wind turbines by Messetii et al. [1] and Turner et al. [6] compared with the power output of our MILP. As the results demonstrate, our proposed layouts generate a better bound of power output. The resulting layouts are shown in Fig. 13. 8.3.2 For this section of case (III), we considered the WFLO Table 9 8.3.1 Comparison of results.
50 m hub height wind turbine 78 m hub height wind turbine Fig. 12. 8.2.3. layouts.
NT
Model
Power (kW)
15
Mosetti et al. [1] Turner et al. [6] Present study Turner et al. [6] Present study
23883 24168 24587 49007 50639
39
296
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
Turner et al. (19 Turbines)
Present study (19 Turbines)
Turner et al. (39 Turbines)
Present study (39 Turbines)
Fig. 13. 8.3.1. layouts.
Table 10 8.3.2 Comparison of results. ðN 1 ; N 2 Þ
Model
Power (kW)
ð8; 17Þ
Chen et al. [4] Present study
33284 34132
Chen et al. [4]
Present study
optimal layout is more than that obtained by the layouts eventuated by Chen et al. [4] (Table 10). The related layouts are demonstrated in Fig. 14. 8.3.3 As the last part of case (III), we solved the WFLO MILP model in Ref. [6] with NT ¼ 25 and z ¼ 78 m. The output power is 36832 kW that has been increased noticeably in comparison with that of the MILP model with different hub height wind turbines ð34132 kWÞ. It shows that there is an optimal state for selecting the number of small hub height wind turbines in a wind farm. For instance, in this case, if we consider N1ðz1 Þ ¼ 5; N 2ðz2 Þ ¼ 20 the power output bound reaches 35060 kW. Fig. 15 shows the difference between the layouts. In all cases, the AIMMS 3.14 [13] modeling language software with the CPLEX [14] solver is used to run the MILP models. 9. Conclusion
50 m hub height wind turbine 78 m hub height wind turbine Fig. 14. 8.3.2. layouts.
model (6.1)-(6.8) with N 1ðz1 Þ ¼ 8; N2ðz2 Þ ¼ 17. Similar to the results in Sections 8.1.2 and 8.2.2, the power output obtained by our
Identical wind turbines
Different hub height wind turbines
In the present study, a new mathematical algorithm has been introduced to model the Jensen’s wake decay in terms of interaction matrix. The WFLO model was extended with different hub height wind turbines. Using the linearization technique, we transformed model to the mixed-integer linear formulation. Using the iterative method of solving the MILP problem which generates a reliable bound in a reasonable time, both MILP models with identical wind turbines and different hub height wind turbines have been investigated. The obtained layouts were as good as the layouts in previous studies. In the case of using wind turbines with different hub heights, our layouts were almost always better than the optimal layouts found in the literature by the genetic algorithm. In addition, we have shown that using different hub height wind turbines in a wind farm leads to an increase in the power output. Nevertheless, there is an optimal state for selecting between various combinations of short and high hub height wind turbines. As a future study, it is worth investigating the relationship between the power output and wind turbine with different hub height selections. Moreover, when different types of wind turbine are considered, the size of WFLO model will be larger. The iterative algorithm provides an opportunity to accommodate more types of wind turbines with a larger size of grid. In addition, the investigation of the optimum location of turbines and placement of the collector in a wind farm is valuable. References
50 m hub height wind turbine 78 m hub height wind turbine Fig. 15. 8.3.3. layouts.
[1] G. Mosetti, C. Poloni, B. Diviacco, Optimization of wind turbine positioning in large windfarms by means of a genetic algorithm, Wind Eng. Ind. Aerodyn. (1994) 105e116. [2] G. Marmidis, S. Lazarou, E. Pyrgioti, Optimal placement of wind turbines in a wind park using monte carlo simulation, Renew. Energy (2008) 1455e1460. [3] Chowdhury Souma, Zhang Jie, Messac Achille, Castillo Luciano, Unrestricted wind farm layout optimization (UWFLO): investigating key factors influencing
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297 the maximum power generation, Renew. Energy (2012) 16e30. [4] Ying Chen, Hua Li, Kai Jin, Qing Song, Wind farm layout optimization using genetic algorithm with different hub height wind turbines, Energy Convers. Manag. (2013) 56e65. [5] S. Donovan, An improved mixed integer programming model for wind farm layout optimization, in: Proceedings of the 41st Annual Conference of the Operations Research Society, Wellington, 2006. [6] S.D.O. Turner, D.A. Romero, P.Y. Zhang, C.H. Amon, T.C.Y. Chan, A new mathematical programming approach to optimize wind farm layouts, Renew. Energy (2014) 674e680. [7] Michele Samorani, The wind farm layout optimization problem, in: Handbook of Wind Power Systems, 2010. [8] N. Jensen, A note on Wind Generator Interaction, RISO National Laboratory,
297
1983. [9] J.F. Manwell, J.G. Mcgowan, A.L. Rogers, Wind Energy Explained, John Wiley & Sons, United Kingdom, 2009. [10] S.A. Grady, M.Y. Hussaini, M.M. Abdullah, Placement of wind turbines using genetric algorithms, Renew. Energy (2005) 259e270. [11] A. Billionnet, E. Soutif, Using a mixed integer programming tool for solving the 0-1 quadratic knapsack problem, INFORMS J. Comput. (2004) 188e197. [12] Samuel Burer, Adam N. Letchford, Non-convex mixed-integer nonlinear programming: a survey, Surv. Oper. Res. Manag. Sci. (2012) 97e106. [13] Johannes Bisschop, AIMMS-optimization Modeling, Paragon Decision Technology, Harlem, 2012. [14] ILOG, ILOG CPLEX 12.4 User’s Manual, 2011. http://www.ilog.com/products/ cplex.
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
Wind farm layout optimization under uncertainty S.A. MirHassani*, A. Yarahmadi Department of Applied Mathematics, Faculty of Mathematics and Computer Science, Amirkabir University of Technology, Tehran, Iran
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 February 2016 Received in revised form 2 January 2017 Accepted 30 January 2017 Available online 1 February 2017
The optimal design of a wind farm with the aim of increasing the exploitation rate and power production is one of the challenging problems in the field of renewable energies. In this paper, the effect of using different hub height wind turbines in a farm on total power is studied. An exact mathematical formulation is presented in terms of the interaction matrix for multi-turbine wake effects considering different hub height wind turbines, and a new mixed-integer quadratic optimization model is developed. Then, the model is generalized to include the changes in wind characteristics and uncertain parameters. The model is solved using the linearization technique and an iterative method. The performance of the model and solution algorithm is evaluated by solving different problems borrowed from the literature. The computation results show the possibility of having a high-quality solution in a reasonable time in almost all cases. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Wind farm layout optimization MIP Wake effect Interaction matrix Uncertainty
1. Introduction Renewable energy is generally defined as energy that comes from natural resources such as sunlight, wind, rain, and waves. This kind of energy is endorsed more than fossil fuel sources for its sustainability and compatibility with climate change and global warming. Wind energy which is used in power production is one of the oldest types of renewable energy. Electrical power is generated by the wind flow in an appropriately designed wind farm. Therefore, determining the optimal location of wind turbines has a significant impact on the performance of the wind farm and the total produced energy. The optimization of the layout of a wind farm is achieved by considering two main factors: the expected power output and the wake effect. The wake effect should be minimized in order to maximize the power output. It has two important features: (i) a decrease in the wind speed; (ii) an increase in the level of turbulence of wind flow. We consider this problem as the wind farm layout optimization (WFLO) problem. 2. Literature review In the past decade, the wind turbine location problem was studied by many researchers. Mosetti et al. [1] solved the WFLO
* Corresponding author. E-mail addresses: [email protected] (S.A. MirHassani), aminyarahmadi@ aut.ac.ir (A. Yarahmadi). http://dx.doi.org/10.1016/j.renene.2017.01.063 0960-1481/© 2017 Elsevier Ltd. All rights reserved.
using a genetic algorithm to find the optimal number and location of wind turbines. They assumed a 10 10 grid of possible turbine locations used Jensen’s wake decay model to analyze the wake effect under various wind speeds and directions. Improving the genetic algorithm and considering the Mosetti objective function, Grady et al. [10] achieved better results. Marmidis et al. [2] used Mossetti and Grady’s objective function in a Monte Carlo simulation approach. Chowdhury et al. [3] determined the location of wind turbines considering the wind farm cost analysis based on the turbine rotor diameters and the number of turbines in the wind farm. They solved the wind farm model using the constrained Particle Swarm Optimization (PSO). Using a nested genetic algorithm, Chen et al. [4] optimized the layout of a 500m 500m wind farm with different hub height wind turbines. They proved that the objective value has improved slightly more than the case with identical hub wind turbines. They also considered different cost models applied to WFLO with different hub height wind turbines, and demonstrated that the result can improve the cost per unit power of a wind farm in this case. A few mathematical programming approaches have been discussed in the literature. Donovan [5] proposed multiple mixed-integer linear models based on vertex packing problem between a couple of wind turbines. Truner et al. [6] developed a new mathematical framework to optimize the layout of a 10 10 grid of possible turbine location in a wind farm. They introduced interaction matrices which modelled the effect of wake interactions between the turbines. They used quadratic integer programming as well as mixed-integer linear programming
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
approaches. In this paper, we propose a mathematical programming framework which is a precise tool for modeling the WFLO problem in order to optimize the layout of a wind farm with different hub height wind turbines. In addition, we formulate an exact wake effect structure and compare the result with the layout in Refs. [6,4,1]. The main factors considered in this study are: (i) the number of wind turbines, (ii) wind directions and speeds, (iii) hub height of wind turbines, (iv) power output, and (v) the exact wake effect interaction. A quadratic integer programming and mixed-integer linear program are developed for different hub height wind turbines. Finally we present a heuristic algorithm with realistic results to solve the proposed models.
3. Wake model When the wind hits a turbine, a cone of slower and more turbulent air is created behind it. This phenomenon, which is known as the wake effect in fluid-aerodynamics (see Fig. 1), has been investigated by researchers [7].
2
a¼
0:5 ln zz0
(3.3)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi a ¼ 0:5 1 1 CT
7 7 7 2 d 5
2a 1þ
aDij
(3.4)
As depicted in Fig. 2, it is preferred to place the turbines in an area with two main characteristics: (i) the highest wind speeds, (ii) outside the wake of other turbines. With the assumption that the terrain is flat, and considering the information in Section 3.1, we use Jensen’s wake decay model as follows: In this study, we assume that the total kinetic energy deficit due to the multiple wake is equal to the sum of individual kinetic energy deficits generated by each wake [7,8]. For more details on the wake model, see Jensen [8]. With this assumption, we can reform Equation (3.1) to calculate the wind speed at position j after multiple sequences of turbines:
2 vj ¼ v0 41
3 sffiffiffiffiffiffiffiffiffiffiffiffi X 2ffi dij 5 Where dij ¼ i2J
3
6 6 v ¼ v0 61 4
289
2a 1þ
aDdij
!2
(3.5)
aDdij þ rr
(3.1)
r1
3.1. The nomenclature
r1 ¼
aDdij
þ rr
(3.2)
All of the parameters involved in the following proposed models and wake model are defined in this section.
J ¼ 1; …; J L ¼ 1; …; L ℕ ¼ 1; …; N S ¼ 1; …; S
Set Set Set Set
K ¼ 1; …; K v0 vj
Set of types of wind turbine, k2K Velocity of free stream (m/s) Wind velocity at position j (m/s)
v r0 rr r1
Ddij
Downstream wind speed (m/s) Turbine radius ðmÞ Downstream rotor radius ðmÞ Wake radius ðmÞ Velocity loss in the wake at position j due to a turbine placed at position i Distance along the wind direction between turbines i and j ðmÞ
Dvij DEij
uij
of of of of
possible turbine locations i2J; j2J wind directions l2L different velocities of free stream n2ℕ scenarios, s2S
NT zk z0 pl pn ps ml a CT
Total number of wind turbines Hub height of wind turbine of type k ðmÞ Surface roughness length ðmÞ Probability of wind directions l Probability of velocity of free stream n Probability of scenario s Slope of line passing through the turbine for each direction
A AO P
Axial induction factor Trust coefficient Entrainment constant Swept area of wind blades ðm2 Þ Overlapping area between wake and a turbine ðm2 Þ Total power output ðkWÞ
Perpendicular distance along the wind direction between turbines i and j ðmÞ
ðXi ; Yj Þ
Coordinates of each turbine ðmÞ
Euclidean distance between turbines i and j ðmÞ
Wij
Interaction matrix
a
4. Interaction matrix calculation The exact structure of wake effect in a wind farm that determines which turbine is in the wake of the upstream turbine plays an important role in optimum wind turbine location. In addition, we need to take into account the wind direction in model formulation. In the present study, we propose a new method to define the exact mathematical structure of the wake effect in a wind farm by means of the interaction matrix concept. The dimensions of interaction matrices are related to the size of grid ðJ JÞ. Therefore, each element of the matrix shows the interaction between itself and J 1 other cells. We attribute a matrix to each direction of incoming wind using the following algorithm: Fig. 1. Schematic representation of the wake effect.
290
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
Fig. 2. The wake effect and a wind farm model.
Table 1 Classification of wind direction and turbines’ coordinates. p q < 3p 4 4 3p q < 5p
Xj < Xi
q < 74p
Xj > Xi
4 5p 4 7p 4
! ! 2 2 2 h2 þ r02 r12 2 1 h þ r1 r0 AO ¼ þ r1 cos 2hr0 2hr1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð h þ r0 þ r1 Þðh r0 þ r1 Þðh þ r0 r1 Þðh þ r0 þ r1 Þ 2 (4.5) r02 cos1
Yj > Yi
4
q < 94p
Yj < Yi
The parameter h indicates the distance between the center of wake effect and the center of turbine blades ðh ¼ Dvij Þ.
4.1. The algorithm i. First, we divide the terrain into four regions. Each region is determined with the direction of incoming wind and the turbine’s coordinates (Table 1). If both conditions in each row of Table 1 are satisfied, then a turbine at position j can be in the wake of a turbine at position i. ii. Formalize the passing line equation of each turbine depending on the wind direction as follows:
y ¼ ml ðx Xi Þ þ Yi
(4.1)
5. Framework with different hub height wind turbines Previous studies proved that using different hub height wind turbines in a wind farm can increase the wind farm’s output. In this section, we present an extended mathematical model for the WFLO problem. In the first part, the Log Law [9] is explained. The details of the generalized model are discussed in the second part. Finally, we appropriately modify the interaction matrices.
iii. According to Fig. 3, compute Dvij for the turbines which meet the first-step conditions, and then calculate Ddij for these turbines considering the Euclidean distance. iv. The value of uij is calculated as follows:
If
If
Dvij > r1 þ rr
uij ¼ 0
Then
r1 rr < Dvij < r1 þ rr
Then
(4.2)
2a
uij ¼ 1þ
aDdij
!2
AO A
aDdij þ rr (4.3)
If Dvij < r1 rr
Then
2a
uij ¼ 1þ
where:
aDdij
!2
(4.4)
aDdij þ rr Fig. 3. Determination of the wake radius at point.
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
5.1. The log law The log law is used to evaluate wind speed from the reference height zr to another level z (see Refs. [4,9]). In this study, the reference height belongs to the turbines which have higher hub heights. Based on the Log Law, one can realize that the wind speed goes up as much as the height increases:
ln zz0 vðzÞ ¼ vðzr Þ ln zr
(5.1)
z0
J K P P
maximize
k¼1 j¼1
291
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 31 u K J uX X ’2 k’ 5A ykj @v0 41 t dkk ij yi 0 2
(5.8)
k’ ¼1 i¼1
P PJ The term Kk¼1 j¼1 ykj v0 is equal to the constant value NT $v0 in which v0 represents the sum of free stream velocities before wind turbines. Therefore, we can replace the maximization model with the minimization one. Moreover considering
’
’
kk ukk ij ¼ dij
2
;
L 2 X ’ ’ ykj ¼ ykj and Wijkk ¼ pl uijkk l
(5.9)
l¼1
vðzr Þ is the free stream velocity before the wind turbines with higher hub heights.
the ultimate quadratic model is formulated as follows:
Minimize 5.2. Mathematical model
J P J P K P K P i¼1 j¼1 k¼1
’
k’ ¼1
’
Wijkk yki ykj
(5.10)
ykj 2Y As already mentioned, to compute the modified constant, we use Equation (3.3) to consider its changes related to turbines with different hub heights. The value of a will change for each turbine and is calculated by the following equation:
ak ¼
0:5 ln zz0k
(5.2)
2
3 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X X k2ffi vj ¼ v0 41 dij 5 where dkij ¼ k
i2J
2a 1þ
ak Ddij
!2 ck2K
ak Ddij þ rr (5.3)
In this case, Equation (3.5) results in the following form: Now, we introduce a new binary variable ykj equal to 1 if a turbine of type k is placed at position j, and to 0 otherwise. The model’s objective function is the generated power which is proportional to the cube of wind speed [6,10]. Thus, our proposed extended model is as follows:
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 313 u K J J uX X K P P ’2 k’ 5A ykj @v0 41 t dkk ij yi 0 2
maximize
k¼1 j¼1
Subject to
J P j¼1
K X
ykj 1
ykj ¼ N k
cj2J
(5.4)
k’ ¼1 i¼1
ck2K
(5.5)
(5.6)
Y shows the feasible region made by constraints (5.5) to (5.7). 5.3. Modified interaction matrix in case of different turbines ’
To evaluate ukk ij , it is essential to determine the type of the upstream turbine. In the case with identical turbines, an interaction matrix was developed for any situation with a different wind direction. In this case, we have jkj2 numbers of interaction matrices for each wind direction. As a result, the modified constant used in Equations (4.2)-(4.4) is specified by the type of turbine at position i, and the parameter h in Equation (4.5) will change as follows (see Fig. 4.):
h¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 Dvij þ ðzk zk’ Þ2
(5.11)
6. Methodology As mentioned in Section 4, we form an interaction matrix for each wind direction. For a given wind direction, the interaction matrix is a sparse and lower triangular matrix. Thus, as more wind directions are considered, the interaction matrix becomes dense. Consequently, solving the model with Formulation (5.10) is very difficult [11]. Two methods are proposed to solve this problem: (i) linearization technique, (ii) direct solution of the quadratic model [6]. We first discuss the linearization technique and generalize it for Formulation (5.10). Then, a heuristic method is introduced to solve the linearized model. 6.1. The linearization technique
k¼1
ykj 2f0; 1g;
cj2J
(5.7)
where the binary variable yki ¼ 1 if a turbine of type k is placed at position i2J The term N k refers to the number of turbines of type k: Thus, the constraint (5.5) requires that N k turbines be placed. The constraint (5.6) enforces that maximum of one turbine be placed at any position. This model is hard to solve due to the nonlinearity of the objective function, and is approximated by applying various transformations. Similarly, we can maximize the wind speed at each turbine instead of maximizing the total generated power [6]. Consequently, the objective function (5.4) will transfer to the following formulation:
minimize
J P J P K P K P i¼1 j¼1 k¼1
Subject to
J P j¼1
K P k¼1
ykj 1
k’ ¼1
ykj ¼ Nk
cj2J
’
’
Wijkk xkk ij
ck2K
(6.1)
(6.2)
(6.3)
292
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
Wake Turbine Overlap Wake
Turbine
Fig. 4. Schematic representation of a wake effect in a wind farm consisting of wind turbines with different hub heights.
’
k xkk ij yi ’
k xkk ij yj
’
’
ci; j2J; isj; ck; k’ 2K
(6.4)
ci; j2J; isj; ck; k’ 2K
(6.5)
’
k k xkk ij yi þ yj 1
ykj 2f0; 1g ’
ci; j2J; isj; ck; k’ 2K
cj2J; ck2K
ci; j2J; isj; ck; k’ 2K
xkk ij 0
(6.6) (6.7) (6.8)
Using this technique leads to the integer quadratic model’s change into a mixed-integer linear (MILP) one. More precisely, we can replace any term of the form yi yj with a new nonnegative variable xij with constraints xij yi , xij yj and xij yi þ yj 1 [12]. By generalizing the definition of the variable xij , we refor’ mulate (5.10) into a mixed-integer linear model where xkk ij ¼ 1 if ’ turbine types k and k are located in positions i and j; respectively, ’ and xkk ij ¼ 0 otherwise. Furthermore, in computational experiments, we add a constraint to accelerate the solution process which states if a turbine of type k and a turbine of type k’ are placed in positions i and j respectively, or j ; i. Therefore: ’
’
kk xkk ij ¼ xji
ci; j2J; isj; ck; k’ 2K
(6.9)
To solve this problem, we consider an iterative method to locate turbines gradually. This technique has two main advantages: (i) the algorithm runs very fast, (ii) the bounds generated by this method are more reliable than those produced by CPLEX through a suboptimal solution. In this technique, the model is first solved by a group (e.g. one, two, or more) of wind turbines, then their location will be fixed. Consequently, the model will locate more turbines in each iteration. The algorithm is terminated when the number of located turbines is equal to the total number of turbines in a wind farm. This procedure works for the MILP model with different hub height wind turbines as well as the case with identical wind turbines. In the first case, the algorithm is allowed to select the type and location of each wind turbine. Moreover, it is possible to determine the type of wind turbines before running the algorithm. To evaluate the efficiency of the iterative method, we considered 36 wind directions with similar probabilities, making the interaction matrix as dense as possible. Furthermore, we increased the size of grid and the number of turbines installed at each stage. The results are shown in Table 2. In addition, the number of turbines installed h i has been determined by 2J in order to challenge the MILP model for moderate values of NT . Table 2 shows the results obtained from the direct solution of the MILP model and iterative method. The first and second columns refer to the values of objective function and running time (s) of the direct method. We placed an upper limit of 24 h on the running time. The next two columns are similarly related to the iterative method, and the Gap column refers to the percentage of difference between these two methods. 7. Computational results
6.2. A heuristic method The computational complexity of the considered MILP models depends on the sparsity of interaction matrix ðWij Þ, size of grid ðJÞ; and the value of NT . We observed that the MILP problem is solvable for a small NT . For a moderate NT , it took too much time and delivered a low-quality solution. By taking into account more wind directions, the density of interaction matrix increases. Therefore, the procedure of solving the MILP problem will be longer and more complicated. As we mentioned, the size and computational complexity are more in the MILP problem with different hub height wind turbines than the case with identical wind turbines. Thus, solving the MILP problem directly is not an appropriate solution.
In this study, two types of wind turbines have been used, k2f1; 2g and, therefore, each wind direction was represented by four interaction matrices. The turbines characteristics and input parameters for WFLO cases are shown in Table 3 (for the cases with identical turbines) and Table 4 (for the cases with different hub height wind turbines). The three cases that have been discussed and reported in the literature have been selected to compare the layout and generated power, see Refs. [1,6,4]. In the first case (I), a dominated wind direction from north to south with the speed of 12 m=s is assumed. The second case (II) consists of 36 wind directions with equal probabilities of occurrence and the speed of 12m=s, and the third
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
293
Table 2 The performance comparison of direct solution and iterative method for linearized formulation (see Ref. [6]). J NT
Direct Method
Iterative Method
Gap
Obj. value
Sol. Time (s)
Obj. value
Sol. Time (s)
94 16 8 25 12 36 18
0:02373819 0:11450618 0:25562266 0:52404598
0:03 0:22 18:52 3278:72
0:023738198 0:114506187 0:256505218 0:524045990
0:00 0:83 2:67 7:25
49 24 64 32 81 40 100 50
0:87460352 1:54624036 2:21679568 3:38630941
23237:74 1005 86400 86400
0:880142860 1:538130101 2:233621949 3:393013736
25:3 55:33 137:35 206:84
0 0 0:0035 ^ 3:05 10ð9Þ 0:006 0:005 0:0076 0:002
Table 3 Parameters used in cases with identical wind turbines. r0 z z0 J Cell size Wind farm size CT P
20 m 60 m 0:3 100 10r0 ð200 m 200 mÞ 2 km 2 km 0:88 PJ 0:3v3i i
Table 4 Parameters used in the case with different hub height turbines. r0 z1 z2 z0 J Cell size Wind farm size CT P
20 m 50 m 78 m 0:3 100 10r0 ð200 m 200 mÞ 2 km 2 km 0:88 PJ;K 0:3vki i;k
Fig. 6. Subdivided wind farm.
subsection, we compare our layouts and generated power value obtained from the WFLO model with identical wind turbines solved by the iterative method, with results from Ref. [1] (MILP) and [6]. The second subsection consists of comparison of results for the case of different hub height wind turbines (6.1)-(6.8) solved by the iterative method and [4]. In the last subsection, we show that the use of wind turbines with different hub heights in a wind farm increases the output with respect to the same case with identical wind turbines. In Section 6.1, we mention that the velocity of incoming wind at the position of wind turbine with a larger hub height is more than that of the smaller one. In case (III), the values of wind speed are assumed to be 8 m=s, 12 m=s and 17 m=s. We suppose that these states represent the velocity of incoming wind at the position of turbines with 78 m of height. Using the Log Law, the velocities of incoming wind at the position of turbines with 50 m of height are 6:44 m=s, 11:04 m=s and 15:64 m=s; respectively. In addition, each state in case (III) is shown as an ordered couple. The first component has two parts which represent the incoming velocity at the
Fig. 5. Case (III), probabilities of occurrence of different wind directions and wind speeds.
one (III) is the case of 36 wind directions with unequal fractions of occurrence and multiple speeds of 8 m=s, 12 m=s and 17 m=s. The probabilities of occurrence for the third case are shown in Fig. 5.
Table 5 8.1.1 Comparison of results. NT
Model
Power (kW)
26
Mosetti et al. [1] Turner et al. [6] Present study Mosetti et al. [1] Turner et al. [6] Present study
9026 9087 9065 9835 9835 9777
8. Case study 30
In this part, we present our results and discuss the performance of our models and solution algorithm in different cases. In the first
294
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
Turner et al. (26 Turbines)
Present study (26 Turbines)
Turner et al. (30 Turbines)
Present study (30 Turbines)
Fig. 7. 8.1.1 layouts.
Table 6 8.1.2 Comparison of results.
Identical wind turbines
ðN 1 ; N 2 Þ
Model
Power (kW)
ð5; 20Þ
Chen et al. [4] Present study
8690 8923
Different hub height wind turbines
heights of 78 m and 50 m; respectively, and the second component demonstrates the direction of incoming wind. Thus, the objective function of MIQP minimization model in case (III) is as follows:
minimize
J P J P 2 P 2 P N P L P k’ ¼1
Subject to
k¼1 ykj 2Y
i¼1 j¼1 n¼1 l¼1
’
’
pn vn0 pl uijkk l yki ykj
(8.1)
50 m hub height wind turbine ’
’
s ¼ vn :ukk l where s corresponds By arranging ps ¼ pn :pl and ukk 0 ij ij to the pair ðn; lÞ, we will have NT L scenarios with the probability P S s of ps . It is clear that s¼1 p ¼ 1. Therefore, Formulation (8.1) is changed to the following MIQP formulation which we can easily transfer to the MILP model:
minimize
J P J P 2 P 2 P S P k¼1 k’ ¼1 i¼1 j¼1 s¼1
p
s
’ ’ uijkk s yki ykj
(8.2)
8.1. Case (I): constant wind speed and direction For this case, we consider a dominated wind direction from the
Chen et al. [4]
Present study
78 m hub height wind turbine Fig. 9. 8.1.3. layouts.
north with a constant wind speed of 12 m=s. The subdivided wind farm and our reference coordinates are shown in Fig. 6. The interaction matrix should be sparse, so the solution time and computational complexity are not too high. 8.1.1 In this section, we compare the WFLO MILP model with identical wind turbines in Ref. [6] solved by the iterative method with the results reported in Refs. [1,6] (see Table 5). In this case, due to the precise calculation of the wake effect (no estimation), there are some differences in the layouts and values of output power in comparison with the results in the literature (Fig. 7). 8.1.2 Table 6 compares the layouts and generated powers of our WFLO model (6.1)-(6.8) solved by the iterative method with the results reported in Ref. [4]. The model was quickly solved for N1ðz1 Þ ¼ 5; N 2ðz2 Þ ¼ 20 and the layout was as good as the one determined in Ref. [4], while the generated power was improved (Fig. 8). 8.1.3 In this part, we investigate the effect of using wind turbines with different hub heights in a wind farm that can improve power generation. For this purpose, we solved the WFLO MILP model in
Table 7 8.2.1 Comparison of results.
50 m hub height wind turbine 78 m hub height wind turbine Fig. 8. 8.1.2 layouts.
NT
Model
Power (kW)
19
Mosetti et al. [1] Turner et al. [6] Present study Turner et al. [6] Present study
6018 6081 6117 9040 9398
39
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
Turner et al. (19 Turbines)
Present study (19 Turbines)
Turner et al. (39 Turbines)
295
Present study (39 Turbines)
Fig. 10. 8.2.1. layouts.
8.2. Case (II): 36 wind directions and constant wind speed
Table 8 8.2.2 Comparison of results. ðN 1 ; N 2 Þ
Model
Power (kW)
ð7; 18Þ
Chen et al. [4] Present study
7189 7412
Chen et al. [4]
Present study
50 m hub height wind turbine
In this case, we assume that the incoming wind is multidirectional with the constant speed of 12 m=s. Each 36 angle has the same possibility of occurrence. The interaction matrices are as dense as possible and, therefore, the solution time and computational complexity are notable. 8.2.1 We fixed the optimal place of 19 and 39 turbines which were taken from Refs. [1,6] in order to calculate the output power in our proposed model (Table 7). The result shows that the output of the generated power is improved in both NT ¼ 19 and NT ¼ 39 cases. Solving the MILP model in Ref. [6] by the iterative method reached the optimal layouts more quickly than the direct method (Fig. 10). 8.2.2 For case (II), we considered the WFLO MILP model (6.1)(6.8) with N1ðz1 Þ ¼ 7; N2ðz2 Þ ¼ 18. Table 8 compares the power output obtained from our optimal layout with those of the layouts created by Chen et al. [4]. The results of the MILP model in this case lead to a layout which is as good as that of Chen et al. (Fig. 11). 8.2.3 For case (II), we considered the WFLO MILP model in Ref. [6] with NT ¼ 25 and z ¼ 78 m. The maximum power output is 7336 kW, while the maximum generated power was 7412 kW for the MILP model with different hub height wind turbines. Fig. 12 demonstrates the resulting layouts.
78 m hub height wind turbine Fig. 11. 8.2.2. layout.
Ref. [6] with NT ¼ 25 and z ¼ 78 m. The output power is 8916 kW while it was 8923 kW for the similar MILP model with different hub height wind turbines. Meanwhile, the layouts in the two cases are almost similar (Fig. 9).
Identical wind turbines
Different hub height wind turbines
8.3. Case (III): 36 wind directions and various wind speeds For this case, Fig. 5 shows 36 wind directions with three speeds of 8 m=s, 12 m=s and 17 m=s with unequal fractions of occurrence. In the case of MILP model with different hub height wind turbines, we suppose that these values of speed represent the velocity of incoming wind at the position of turbines with 78 m of height. Using the Log Law, the velocities of incoming wind at the position of turbines with 50 m of height are 6:44 m=s, 11:04 m=s and 15:64 m=s; respectively. The density of interaction matrices is the same as in case (II). 8.3.1 Table 9 presents the power output obtained from the optimal layouts of 15 and 39 wind turbines by Messetii et al. [1] and Turner et al. [6] compared with the power output of our MILP. As the results demonstrate, our proposed layouts generate a better bound of power output. The resulting layouts are shown in Fig. 13. 8.3.2 For this section of case (III), we considered the WFLO Table 9 8.3.1 Comparison of results.
50 m hub height wind turbine 78 m hub height wind turbine Fig. 12. 8.2.3. layouts.
NT
Model
Power (kW)
15
Mosetti et al. [1] Turner et al. [6] Present study Turner et al. [6] Present study
23883 24168 24587 49007 50639
39
296
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297
Turner et al. (19 Turbines)
Present study (19 Turbines)
Turner et al. (39 Turbines)
Present study (39 Turbines)
Fig. 13. 8.3.1. layouts.
Table 10 8.3.2 Comparison of results. ðN 1 ; N 2 Þ
Model
Power (kW)
ð8; 17Þ
Chen et al. [4] Present study
33284 34132
Chen et al. [4]
Present study
optimal layout is more than that obtained by the layouts eventuated by Chen et al. [4] (Table 10). The related layouts are demonstrated in Fig. 14. 8.3.3 As the last part of case (III), we solved the WFLO MILP model in Ref. [6] with NT ¼ 25 and z ¼ 78 m. The output power is 36832 kW that has been increased noticeably in comparison with that of the MILP model with different hub height wind turbines ð34132 kWÞ. It shows that there is an optimal state for selecting the number of small hub height wind turbines in a wind farm. For instance, in this case, if we consider N1ðz1 Þ ¼ 5; N 2ðz2 Þ ¼ 20 the power output bound reaches 35060 kW. Fig. 15 shows the difference between the layouts. In all cases, the AIMMS 3.14 [13] modeling language software with the CPLEX [14] solver is used to run the MILP models. 9. Conclusion
50 m hub height wind turbine 78 m hub height wind turbine Fig. 14. 8.3.2. layouts.
model (6.1)-(6.8) with N 1ðz1 Þ ¼ 8; N2ðz2 Þ ¼ 17. Similar to the results in Sections 8.1.2 and 8.2.2, the power output obtained by our
Identical wind turbines
Different hub height wind turbines
In the present study, a new mathematical algorithm has been introduced to model the Jensen’s wake decay in terms of interaction matrix. The WFLO model was extended with different hub height wind turbines. Using the linearization technique, we transformed model to the mixed-integer linear formulation. Using the iterative method of solving the MILP problem which generates a reliable bound in a reasonable time, both MILP models with identical wind turbines and different hub height wind turbines have been investigated. The obtained layouts were as good as the layouts in previous studies. In the case of using wind turbines with different hub heights, our layouts were almost always better than the optimal layouts found in the literature by the genetic algorithm. In addition, we have shown that using different hub height wind turbines in a wind farm leads to an increase in the power output. Nevertheless, there is an optimal state for selecting between various combinations of short and high hub height wind turbines. As a future study, it is worth investigating the relationship between the power output and wind turbine with different hub height selections. Moreover, when different types of wind turbine are considered, the size of WFLO model will be larger. The iterative algorithm provides an opportunity to accommodate more types of wind turbines with a larger size of grid. In addition, the investigation of the optimum location of turbines and placement of the collector in a wind farm is valuable. References
50 m hub height wind turbine 78 m hub height wind turbine Fig. 15. 8.3.3. layouts.
[1] G. Mosetti, C. Poloni, B. Diviacco, Optimization of wind turbine positioning in large windfarms by means of a genetic algorithm, Wind Eng. Ind. Aerodyn. (1994) 105e116. [2] G. Marmidis, S. Lazarou, E. Pyrgioti, Optimal placement of wind turbines in a wind park using monte carlo simulation, Renew. Energy (2008) 1455e1460. [3] Chowdhury Souma, Zhang Jie, Messac Achille, Castillo Luciano, Unrestricted wind farm layout optimization (UWFLO): investigating key factors influencing
S.A. MirHassani, A. Yarahmadi / Renewable Energy 107 (2017) 288e297 the maximum power generation, Renew. Energy (2012) 16e30. [4] Ying Chen, Hua Li, Kai Jin, Qing Song, Wind farm layout optimization using genetic algorithm with different hub height wind turbines, Energy Convers. Manag. (2013) 56e65. [5] S. Donovan, An improved mixed integer programming model for wind farm layout optimization, in: Proceedings of the 41st Annual Conference of the Operations Research Society, Wellington, 2006. [6] S.D.O. Turner, D.A. Romero, P.Y. Zhang, C.H. Amon, T.C.Y. Chan, A new mathematical programming approach to optimize wind farm layouts, Renew. Energy (2014) 674e680. [7] Michele Samorani, The wind farm layout optimization problem, in: Handbook of Wind Power Systems, 2010. [8] N. Jensen, A note on Wind Generator Interaction, RISO National Laboratory,
297
1983. [9] J.F. Manwell, J.G. Mcgowan, A.L. Rogers, Wind Energy Explained, John Wiley & Sons, United Kingdom, 2009. [10] S.A. Grady, M.Y. Hussaini, M.M. Abdullah, Placement of wind turbines using genetric algorithms, Renew. Energy (2005) 259e270. [11] A. Billionnet, E. Soutif, Using a mixed integer programming tool for solving the 0-1 quadratic knapsack problem, INFORMS J. Comput. (2004) 188e197. [12] Samuel Burer, Adam N. Letchford, Non-convex mixed-integer nonlinear programming: a survey, Surv. Oper. Res. Manag. Sci. (2012) 97e106. [13] Johannes Bisschop, AIMMS-optimization Modeling, Paragon Decision Technology, Harlem, 2012. [14] ILOG, ILOG CPLEX 12.4 User’s Manual, 2011. http://www.ilog.com/products/ cplex.