Wind Farm Yaw Optimization via Random Search Algorithm

One direction in optimizing wind farm production is reducing wake interactions from upstream turbines. This can be done by optimizing turbine layout as well as optimizing turbine yaw and pitch angles. In particular, wake steering by optimizing yaw angles of wind turbines in farms has received significant attention in recent years. One of the challenges in yaw optimization is developing fast optimization algorithms which can find good solutions in real-time. In this work, we developed a random search algorithm to optimize yaw angles. Optimization was performed on a layout of 39 turbines in a 2 km by 2 km domain. Algorithm specific parameters were tuned for highest solution quality and lowest computational cost. Testing showed that this algorithm can find near-optimal (<1% of best known solutions) solutions consistently over multiple runs, and that quality solutions can be found under 200 iterations. Empirical results show that as wind farm density increases, the potential for yaw optimization increases significantly, and that quality solutions are likely to be plentiful and not unique.


Introduction
Wind turbine wakes are undesirable as they reduce energy production and increase mechanical stresses of downstream turbines. These wake interactions can reduce annual energy production up to 10% to 20% [1]. As wind passes through a wind turbine it creates a wake that detrimentally affects the power generated by turbines downstream. One way to reduce wake effect during planning stage is place wind turbines strategically within a domain, and this has been well studied [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. Many of these studies focused on developing novel optimization algorithms to minimize wake effects. Other studies focused on developing optimization models and algorithms to solve multi-objective functions beyond purely minimizing wake effects, such cost, cabling, land constraints, noise.
In addition to optimizing turbine locations, yaw optimization and control of individual turbines can also reduce wake loss and improve the power production [22][23][24][25][26][27][28][29][30][31][32]. By yawing upstream turbines, wakes can be deflected away from downstream turbines. In particular, several studies [28,29,31,32] were conducted to study the potential improvement of controlling turbine yaw angles. Fleming et al. [24] performed an numerical analysis on two turbines and concluded that the overall power production could be increased with optimized yaw angles. There are some existing work in yaw angle optimization [25,28,32] which demonstrated that yawing optimization can reduce wake losses significantly, and that optimizing yaw angles based on wind direction is more important than optimizing to wind speed. Although yawed turbines will produce less energy, the downstream wakes are deflected away from turbines downstream, leading to an overall increase in energy production. However, optimizing yaw angles within a farm is an non-trivial task, as there are infinite combinations of yaw angles in every wind farm.
The aim of wind farm yaw control and optimization is to maximize power production by minimizing wake interactions. One approach is to accurately describe the complex physical interactions within the wind farm, and directly develop control strategies based on this description [23,27,32,33]. However, due to the complexity of turbulent atmospheric flow, the computational cost required to accurately simulate all the conditions is immense, thus real-time wind farm power optimization using high-fidelity flow simulations is not yet possible [29,34].
In recent years, low cost wake models capable of simulating yawed wakes which are suitable for fast yaw optimization have been developed. Jiménez et al. [23] derived a simple wake deflection model of turbines in yaw, but subsequent studies [35][36][37] demonstrated that the model tends to overestimate wake deflection. Bastankhah and Porté-Agel [26] developed a computationally inexpensive analytical wake model to predict wake deflection and velocity distribution in the far wake; however, this model does not adequately account for turbulence intensity, which is important in wake development. Qian and Ishihara [35] developed a mechanistic analytical wake model accounting for turbulence intensity, and validated with experimental results. Then, Lopez et al. [38] developed an accurate and low-cost numerical field model, building on the unyawed wake model by Kuo et al. [39] and the yawed wake model by Qian and Ishihara [35]. Shapiro et al. [40] derived a new model based on the aerodynamic lifting line theory.
In conjunction with these wake models and wind farm data, control algorithms can be used to optimize yaw angles. Several studies have developed optimization-based control approaches to adjust to the uncertainties, using approaches such as Bayesian optimization [41], genetic algorithm [42], game theory [43,44], random search [45], sequential optimization [32], gradient descent [28,46], greedy control [47], particle swarm optimization [48], dynamic programming [49,50]. Most recently, Howland et al. [37] developed an control scheme based on an empirical fitted analytical wake model [26,40]. This analytical wake model was calibrated for a six turbine group within a wind farm in Alberta, Canada based on historical data. Then wake steering is done by solving a yaw optimization model in real-time using a gradient descent strategy called Adam optimization [51]. Over the period tested, the wind farm production increased by 7-13% near average wind speeds and decreased production variability by up to 72%.
To solve the optimization problem in real-time, particularly for higher number of turbines, it is important to reduce computational cost. A number of optimization algorithms have been developed to optimize wind farm layout [6,16,[52][53][54][55]. Keeping in mind of the lessons learned from layout optimization literature, the objective of this work is to develop an optimization algorithm to optimize yaw angles of a wind farm, with emphasis on solution quality and low computational cost. In this work, we propose a random search (RS) algorithm, and we apply this algorithm in conjunction with a steady analytical wake model by Qian and Ishihara [35] to compare its performance with that of a binary quadratic programming model, in terms computational cost and solution quality. The remainder of the paper is structured as follows: Section 2 focuses on the wake model used in this work; Section 3 discusses the proposed RS algorithm and the model to be solved; Section 4 illustrates and discusses the performance of the algorithm; and Section 5 summarizes the findings.

Wake Model
This work is utilizes the wake model by Qian and Ishihara [35] with the proposed optimization algorithm. This yaw wake model states the relationship between the average wake velocity, U w , and velocity deficit, ∆U, and is described as where, U 0 is the freestream velocity of the wind.
As shown in Figure 1, ∆U is related to both x (i.e. the distance measured downstream from the turbine) and y d (i.e. radial distance measured from the center of the yawed wake) in x − y plane, and ∆U can be found by using ∆U = Sφ (2) where S in Equation (2) is the streamwise function, and can be found by using, Here, γ represents the yaw angle, D represents the diameter of turbine (illustrated in Figure 1), I a is the ambient turbulence intensity, and C T is the coefficient of thrust under no yawing conditions. In this form, Equation (3) is applicable in the far wake region. x − y plane of Gaussian wake as described by Qian and Ishihara [35].
To capture the change in velocity deficit (Equation (2)) in the spanwise direction, a spanwise function φ is used, and this is described using a Guassian function, Here, σ is the standard deviation of the wind velocity distribution at each cross-section, and can be calculated through where r in Equation (4) is the radial distance measured from the center of the yawed wake in x − z plane, as illustrated in Figure 2, and can be defined by Here, H is the hub height; y d is the wake deflection, and its expression depends on its location within the wake. In the near wake region, y d is assumed to be linear with distance downstream of turbine and is expressed as and, where θ 0 is the initial skew angle. In the far wake region, y d is expressed as where y d0 is the wake deflection and σ 0 is the standard deviation of the velocity distribution at the end of near wake region, x 0 , and and x − z plane of Gaussian wake as described by Qian and Ishihara [35].
This analytical model can only be used to calculate the downstream velocity caused by a single yawed turbine. If a wake of ith turbine is influenced by n multiple nearby upstream yawed turbines, the downstream velocity, U i , needs to be calculated through where, U ij represents the wind velocity at ith turbine due to the wake of jth turbine. If the turbine is partially immersed in an upstream turbine wake, then A w,j is the rotor area immersed in the upstream wake. Details on A w can be found in the work by González et al. [56].

Optimization Formulation
The objective of the optimization problem is to maximize power production of a wind farm. The problem formulation for the optimization problem is as follows: where the number of turbines in a wind farm is N and n is the number of turbines upstream of turbine i. A single wake model is used to solve for the velocity at turbine i due to the wake of turbine j, denoted as U ij . In this formulation, each turbine i is allowed to yaw freely, at angle γ i . Furthermore, ρ is density of air, A is rotor area of turbine, and C P,i is the power coefficient of the turbine i. The power curve of a turbine is dependent on the model and manufacturer, but an ideal curve is used in this study, and we assume that the coefficient of power of a yawed turbine is proportional to cosine of yaw angle squared, similar to empirically observations [25,57,58]. That is, C P,i = C P,unyawed cos 2 γ at an yaw angle of γ. Note that there is no physical explanation behind this cosine dependence other than it created a good fit with experimental results.

Proposed Random Search Algorithm
The authors initially developed an BQP model [59], which allows turbines to yaw a predefined set of angles. The authors' experience developing and working with mathematical programming models for wind farm layout optimization problems suggested that this BQP model can be quickly solved using state-of-the-art solvers such as Gurobi and CPLEX, and perhaps outperform certain types of metaheuristic algorithms. However, it became clear that the BQP model is not suitable for fast optimization purposes as the size of the problem grows exponentially with increasing size of the set of angles. From our experience developing the BQP model and analyzing results, we realized that near-optimal (<1% of best known solutions) solutions are plentiful and that systematically guessing solutions is likely to do well. This leads to the proposed RS algorithm.
The core idea behind the proposed algorithm is to improve solution quality based on best known solution during the optimization process. Thus, instead of determining the yaw angle, γ, the algorithm determines an incremental yaw angle with respect to the current yaw angle. This incremental yaw angle is randomly selected from a range of −β • to β • , where the value of β is predefined. As shown in Figure 3, the algorithm starts with the determination of the power output of a farm with all non-yawed turbines, which is considered as the initial farm configuration. This initial configuration is considered to be the best known configuration. Then, the yaw configuration of the farm is changed by yawing F% of the turbines to angles ranging from −β • to β • . After that, the power output of the farm with new configuration is calculated and compared to that of the best configuration. If the new configuration produces higher power production, then the best known configuration will be replaced by the new configuration. Based on the best known solution, the turbines will yaw additionally and randomly between the range of −β • and β • and compared with the previous solution. This process is repeated for a specified number of iterations. In other words, the yaw angle of a turbine i, γ i , will be the sum of all the random angles chosen between −β • and β • from previous iterations. For example, suppose β = 20 • , and turbine 1 yaws 9 • (random selected between −20 • to 20 • ) from the initial yaw angle of 0 • relative to the wind in the first iteration, then yaws 8 • (random selected between −20 • to 20 • ) from its current yaw angle in the second iteration, so the yaw angle of the turbine after two iterations is 9 • + 8 • = 17 • , if both configurations produced higher larger power output than previous configurations.
This simple algorithm builds upon best known solutions to find the better solutions. To reduce the computational cost and improve the quality of solutions, determination of the percentage of turbines to be yawed, F, and the range of angle [−β • , β • ] are crucial to achieve those objectives. For instance, if β is large, there is a strong likelihood to find very different solutions than the initial baseline solution. Conversely, if β is small, it is likely that the new solution will be similar to the initial baseline solution. Thus, the value of β can be seen as a mean to control intensification and diversification of the optimization process. The percentage of turbines yawed in a single iteration, F, is important as well in terms of solution quality and computational cost. If F is a small value, the change per iteration will be small and may likely slow down convergence and further increase computational cost. If the value of F is large, the additional computational cost associated with each mathematical operation may not become negligible, but there is no guarantee that the solution quality will increase, especially during the intensification process. The detailed explanation of choice of F and β is illustrated in the Section 4. The termination criteria in the proposed algorithm is dependent on the number of iterations. Of course, more complex convergence criteria could be implemented; however, since computational cost is of importance in real-time yaw control, constraining the computational cost by fixing the maximum number of iterations is also applicable. Clearly, as the value of maximum iteration increases, better solutions will always be found, thus we are interested in determining the best possible solutions within a realistic time frame. The value of maximum iteration was set to be 1000 iterations.

Results
The optimal wind farm layout used in this work is from Grady et al., shown in Figure 4. This optimized layout is based on a wind speed of 12 m/s uniformly from 36 different wind direction sectors separated by 10 • , and each wind direction sector is a 10 degree wind sector. For example, 0 • wind direction sector is a wind sector ranging from −5 • to 5 • . The diameter of the 39 turbines is D = 60 m , hub height of 60 m, and C T = 0.88. These parameters are shown in Table 1. The optimal layout from this simple wind regime shows turbines evenly spread out throughout the domain and that there are no turbine clusters. We will be performing yaw optimization on this layout to study the performance of the RS algorithm.
It is important to note that in this study, we only consider fixed values for C T and C P . However, under experimental condition, these parameters are functions of tip speed ratio and yaw angle. Thus, these simplified assumptions would lead to errors in predicting power production. Nevertheless, as this work is focused on studying the performance of optimization approaches, these simplifications are justified as they are commonly used in the literature. The algorithm in this paper was implemented in MATLAB (i.e. MATLAB 2019) and solved using an HP EliteDesk with a 4 core Intel i5-6500 3.2 GHz processor and 8 GB of RAM. To ensure a high quality of optimal solution and low computational cost, the maximum number of iterations is set to be 1000. Under current MATLAB implementation, it takes approximately 3 s per iteration on the 4-core desktop. It is clear that even at 1000 iterations, the computational cost for RS algorithm is significantly lower than the 16 CPU hours used to solve the BQP model. The optimization problem is solved for each wind direction at 0.1 • degree increment, and then averaged over their corresponding sector, for all 36 different wind direction sectors. The results of this work is compared to that of BQP model, discussed in the following subsections. should be used if we wish to explore the search space (diversification) and the smaller values of β (e.g. 5 • and 2 • ) if we wish to concentrate on promising regions (intensification). During the diversification process, the focus is to find a number of good solutions earlier in the optimization process. Then, later in the process, an intensification process would take place to refine solution quality. Thus, it may not be ideal to consider a single β value throughout the the optimization process. Thus, we consider a decaying function for β, which varies the value of β during the optimization process to balance between solution quality and computational time, and the results are compared with fixed β values. The idea of this decaying function is similar in concept to the probabilistic function in simulated annealing [61]. The value of β should be larger early on in the optimization process to take risks and explore different neighborhoods, and gradually decrease with increasing number of iterations. Then, we define β as β = 7 • e −kg/1000 + m • , where k is an exponential constant to be determined and g is the number of iteration, and has an asymptote of m • . Note that this equation was developed with an termination point of 1000 iterations (hence the exponent −kg/1000), and that this should be modified for different maximum iteration. Here we consider the value of m to be 4, and k to be 5 and 8. We compare the results with the decay function with constant β values, shown in Table 2, which shows the percent improvement in power (solution quality) after 1000 iterations over five trials, with the top five objective values underlined. It can be seen that in terms of solution quality, the results using the decay function are similar to that of fixed β values. Three important observations can be drawn from these experiments. First, most trials reach nearly identical objective values. There are likely numerous local optima with similar objective values, and that all approaches will eventually find one. Second, when β value is small, e.g., 2 • , the number of iterations it takes to find near-optimal (<1% of best known solutions) solution increases. As the β value increases to 5 • , the number of iterations required to find near-optimal solution decreases. However, when the value of β increases to 10 • and 20 • , although good solutions are found early on, solution quality improvement becomes stagnant later in the optimization process. Third, all β values tried with the exception of 2 • all resulted in similar performance, ranging between 6.3% to 6.5%. This means that the algorithm is robust in finding good solutions and that diversification is likely to be more important than intensification in this type of problem.

Identifying F
Identifying the percentage of turbines, F, to be yawed in each iteration is important. It is worth noting that as the value of F increases, the ability to fine tune solutions decreases as more randomness is present in each iteration, especially later on in the optimization process. The value of F is varied to study how it affects solution quality and computational cost. We performed five trials for 13 different values of F, corresponding to different number of turbines which will be yawed in each iteration. These trials were performed with a decay function constants of k = 5 and m = 4. As previously mentioned in Section 3, when the value of F is small, it would take many iterations to find good solutions as it searches through the search space at a slower pace. When the value of F is very large, randomness will dominate as a single turbine could be yawed multiple times within a single iteration, thus producing similar effects as increasing the range of β. Table 3 shows the average solution (of five trials) found at the end of 100, 200, 300, 500, and 1000 iterations. The values are shown in terms of percentage increase from unyawed solution, averaged over five trials. The best solutions after 50 iterations are found with F values in the range of 17.9% and 51.3%. After 200 to 300 iterations, the best range of F lies between 7.7% to 51.3%. However, after 500 iterations, the best F values reduces to a range of 7.7% to 25.6%. After 1000 iterations, the ideal range lies between 5.1% to 25.6%. In other words, higher F value emphasizes on exploring the search space, and will produce better solutions earlier on in the search, while lower F value intensifies the search, which finds better solutions later on in the optimization process. However, after 1000 iterations, the quality of solutions found for all F values tested, with the exception of F = 2.5% and F = 100%, resulted in greater than 6% increase in performance, illustrating the robustness of the algorithm. For the remainder of this work, the values F = 5.1, k = 5, and m = 4 are used. Table 3. Percent improvement in objective value (from unyawed solution) as function of number of iterations with varying F values at D = 60 m with northerly wind (0 • direction sector, 10 degree sector, averaged over 0.1 • increment), averaged over 5 trials, with k = 5 and m = 4. Best objective value at each iteration is underlined.

Solution Reproducibility
Randomness is an inherent part of the proposed RS algorithm. To test the RS algorithm's sensitivity to initial seeds, five trials were performed to demonstrate that the algorithm is able to explore the solution space and that the solution quality is insensitive of initial guesses. The results are shown in Figure 5, and two observations can be made. First, potential improvement is strongly dependent on wind direction. In particular, it is clear that production improvements are greatest when winds are blowing from cardinal and intercardinal directions, when turbines are aligned due to the structured-grid nature of the layout. Second, in each of the five trials, the improvements are nearly identical for all wind direction sectors. At first glance, this was very surprising as this implied that the RS algorithm is able to converge to similar objective values with different initial guesses. Upon closer inspection of the five solutions, it is likely that this type of problem has multiple local optima with similar solution quality and that the RS algorithm is able to adequately explore the solution space.

Solution Quality
To the best of the authors' knowledge, the only existing parametric studies on wake steering is the work by Kuo et al. [59], which illustrated that the potential of wake steering is a function of wind farm density. That work employed a binary quadratic programming (BQP) formulation. The formulation only allows turbines to yaw a discrete number of angles. Although this discrete formulation limited the quality of solutions, the formulation allowed well-developed mathematical programming methods to be applied, as the BQP solution is known to be optimal. This work can be compared to the BQP formulation in terms of solution quality. Table 4 shows the power density and the average percent improvement of the proposed work and the BQP model, and Figure 6 shows the comparison for different wind direction sectors at D = 60 m. Two observations can be drawn from Table 4. First, as power density increases with turbine diameter, production improvement also increased, as previously observed in the work by Kuo et al. Second, the production improvement from RS algorithm significantly outperforms the BQP model. There are two reasons for this, (1) although the RS algorithm restricts yawing within a range of incremental angle β for each iteration, it does not restrict turbines from yawing beyond that range. For example, even if β is restricted to 20 • for each iteration, it is possible for a turbine to have a final yaw angle of 47 • after a number of iterations. However, this is not the case for the BQP formulation, where the range of angle is a global constraint.
(2) The RS algorithm permits a more dynamic range of angles (e.g. varying β values), while the BQP formulation only allows a predefined set of angles throughout the optimization process. Table 4. Power density and percent production improvement of layout from Grady et al. [60] using binary quadratic programming (BQP) model [59] compared with proposed random search (RS) algorithm.

Sensitivity to Wind Direction
It is known that wake effects are very sensitive to changing wind directions, thus in this work, we perform simulations at 0.1 • resolution, then the results are averaged over 10 degree wind sectors. This approach is important as this particular wind farm layout is very structured such that turbines are mostly aligned in certain angles. This means that if wakes from upstream turbines are directly aligned with downstream turbines due to perfect alignment with wind direction, then upstream turbines would need to yaw significantly to reduce wake interactions downstream. This would also reduce overall gains as upstream turbines would produce less power as a result of larger yaw angles. If the wind is slightly misaligned with turbine positions, then improvements from yaw optimization should be larger. We illustrate this sensitivity by presenting the percentage improvement at 0 • wind sector, or wind directions from −5 • to 5 • at 0.1 • increment, shown in Figure 7. As discussed, it is clear that when wind direction moves away from 0 • (perfect alignment due to structured layout), the percent improvement from yaw optimization increases as turbines move out of perfect alignment from incoming wind. Although in this work, the approach of averaging the solutions at 0.1 • increment is performed, it is not possible in practice to change the yaw angles for every 0.1 • change in wind direction. Furthermore, there is no guarantee that a near-optimal solution for one particular wind direction will be a near-optimal solution within the wind sector. Thus, the future research direction is to develop a robust yaw optimization model, similar to that of the robust wind farm layout optimization model developed [17] for varying wind directions, and that the proposed RS algorithm can be used to solve the robust model.

Conclusions
This paper describes a random search algorithm capable of optimizing yaw angles of wind turbines. A series of tests have been conducted to understand the algorithm's performance in terms of computational cost and solution quality. The proposed RS algorithm has been shown to significantly outperform previous work of a BQP formulation solved using a state-of-the-art mathematical programming solver, and the production improvement is tripled compared with BQP algorithm. Through the empirical study, the following observations have been made, (a) The quality of solutions is consistent and is reproducible over multiple runs with minor variations; (b) The solution quality is dependent on the range of yaw angles and the percentage of turbines yawed in each iteration; (c) During optimization process, it is more preferable to have larger range of yaw angle at the beginning and gradually tighten the range as the process continues to renew the results; (d) The optimal range of the percentage of turbines yawed in each iteration is between 5% and 50%; and (e) Near-optimal solutions (<1% of best known solutions) can be found within 500 iterations, and that excellent (<2% of best known solutions) solutions can be found under 300 iterations. While testing the capabilities of the proposed algorithm, four observations related to wind farm yaw optimization were made. First, as power density increases, the potential to improve production through yaw optimization increases. Second, good solutions are likely to be plentiful and are not unique. Third, yaw optimization can further improve power production in an optimized layout. Finally, yaw optimization is sensitive to changes in wind direction, so robust yaw optimization considering varying wind direction should be considered in future work.