A Model-Free Approach for Maximizing Power Production of Wind Farm Using Multi-Resolution Simultaneous Perturbation Stochastic Approximation

This paper provides a model-free approach based on the Multi-Resolution Simultaneous Perturbation Stochastic Approximation (MR-SPSA) for maximizing power production of wind farms. The main advantage is that the method based on MR-SPSA can achieve fast controller tuning without any plant model by exploiting the information of the wind farm configuration such as turbines location and wind direction. In order to simulate the performance of the model-free scheme, a wind farm model with dynamic characterization of wake interaction between turbines is used and then the proposed method is applied to the Horns Rev wind farm. Simulation results illustrate that the method based on MR-SPSA achieves the maximum total power production with faster convergence compared with other existing model-free methods.


Introduction
Nowadays, wind energy has been one of the most cost-efficient sources of renewable energy. Therefore, control algorithms to maximize energy production of wind farms have been actively studied in the control community.
Over the past two decades, most of the researches in wind turbine control focus on improving the control algorithm of a standalone turbine, which has been well reported in [1][2][3]. On the other hand, the control of an array of turbines in a wind farm is very challenging, because the aerodynamic interactions among the turbines are complex and difficult to model. Hence, the control algorithm for the standalone turbine needs to be improved when it is placed in a wind farm. Many researchers have invested their effort to model the dynamics of the wind farms and synthesize the control law of the wind farms based on the developed models. In [4], a centralized power control of wind farm has been proposed. In this approach, the central control level decides the power reference for each individual turbine at the local control level. In [5,6], the centralized model predictive control technique has been proposed to achieve a certain level of wind farm power production. There, the controllers predict the wake propagation of the derived wind farm models. Similar controllers have been reported by [7,8] in a distributed way, where the power production reference for each turbine is given in advance. For example, in [8], if the wind speed is higher than the nominal one, the references for both the power production and the blade pitch angle are provided to the individual turbine, otherwise rotor speed references are given to each turbine to extract the available power. Next, a stationary wind farm model has been introduced by [9], which consists of the ambient wind speed and interactions among turbines sub-models. Here, the controller adjusts the blade pitch angle of the individual turbine according to its characteristics regarding wind speed and power. In [10], an optimal control method has been proposed for different classes of wake interaction model (e.g., wake propagation under densely and sparsely spaced turbine arrays).
However, the above model-based control strategies must be difficult to apply in practice. This is because their control methods rely on simple linear time invariant models with static wind conditions, which do not accurately describe the chaotic nature of wind and the complex aerodynamic interaction among turbines. Therefore, a model-free approach will be more attractive. So far, several model-free approaches for wind farm optimization problems have been investigated. In [11] and its extended version [12], Game Theory (GT) based approaches have been proposed for optimizing power production in a wind farm. There, distributed learning algorithms have been applied with both local and global knowledge, which are available to each turbine in a wind farm configuration. Another type of GT based approach has been proposed in [13] by considering two control variables of wind turbines. However, these approaches are validated only for a static model of the wind farm, and rely on off-line computation. Therefore, the time efficiency of their learning algorithms, especially online performance, is not guaranteed. To obtain an alternative solution, a dynamic model of the wind farm has been proposed in [14] by approximating a delay structure when a wake travels from one turbine to another, where a model-free approach has been proposed based on the Fixed-Step Maximum Power Point Tracking (FS-MPPT), in which local knowledge is assumed to be available to each turbine. The results show much faster convergence with acceptable degradation in the maximum power production compared with the GT-based approach in [11].
On the other hand, a Simultaneous Perturbation Stochastic Approximation (SPSA) would provide us a promising model-free approach for maximizing the power production of the wind farms. This is because the SPSA method is known to be effective for a variety of model-free optimization problems even for high-dimensional parameter tuning [15]. However, it is not clear whether it works for wind farm problems, since there are few literatures to discuss the application of the SPSA to the problems.
In this paper, we develop a Multi-Resolution Simultaneous Perturbation Stochastic Approximation (MR-SPSA) algorithm and apply it for maximizing power production of the wind farms. Note that the application of the standard SPSA-based method for solving the wind farm optimization problem is still not enough for providing an acceptable convergence speed. For example, in the case of the Horns Rev wind farm with 80 turbines, as we show in this paper, the standard SPSA requires more than 100 h to maximize its total power production, which is heavy computation time. The detail will be shown in Section 5, where it turns out that the exhaustive local search of high-dimensional design variables in the SPSA-based method may not be a practical solution from this point of view. Therefore, we develop the MR-SPSA-based method, which results in a much faster convergence by exploiting the wind farm configuration information. In order to evaluate the performance of the model-free approach, the wind farm model in [14], which has dynamic characterization of wake interaction between turbines, is adopted in this study. The proposed method is tested in an 80-turbine wind farm that replicates the layout of an actual offshore Horns Rev wind farm. Then, the performance of the proposed method is analyzed in terms of the maximum total power production of wind farm and the convergence speed. Finally, a comparative assessment between the MR-SPSA-based method, the standard SPSA-based method, the FS-MPPT-based method [14] and the GT-based method [11] is presented in detail.
The rest of the paper is organized as follows. Section 2 formulates the problem of maximizing power production of wind farms. In Section 3, the methodology and benefit of the MR-SPSA-based method algorithm are shown. Next, the model-free design method for maximizing the total power production of wind farm is presented in Section 4. The wind farm model is briefly described in Section 5. Then, the effectiveness of the MR-SPSA-based method is demonstrated through simulation of the Horns Rev wind farm. Further, the statistical analysis and performance comparison between the MR-SPSA-based method and other existing methods for various wind farm configurations are also presented. Finally, some concluding remarks are given in Section 6.

Problem Formulation
The wind farm to be considered consists of p wind turbines. The turbines are placed arbitrarily in a region. Let a i (i = 1, 2, ..., p) be the control parameter of turbine i and P i (a 1 , a 2 , ..., a p ) (i = 1, 2, ..., p) be the power production of turbine i. Here, a i corresponds to the axial induction factor, which is a generalized form of the control parameters, i.e., blade pitch and rotor speed of turbine i [16]. Note that a wind field with an arbitrarily wind direction and non-static speed (but slow change) is considered in this study. Hence, it is expected that the axial induction factors other than turbine i, which are a 1 , a 2 , ..., a i−1 , a i+1 , ..., a p , would also influence the power production P i due to wake interaction between turbines. Similarly, any changes of axial induction factor a i influence not only P i but also power productions from other turbines, which are P 1 , P 2 , ..., P i−1 , P i+1 , ..., P p . In particular, we can say that P i strongly depends on a i and more or less depends on a 1 , a 2 , ..., a i−1 , a i+1 , ..., a p . Moreover, it is difficult to accurately model the relation among P i and a 1 , a 2 , ..., a p because of the complicated interactions between the wake aerodynamics and the turbine dynamics. However, the total power production of the wind farm is assumed to be measured, which is formally expressed as: Then, our problem can be described as follows: Problem 2.1 Suppose that the explicit forms of the functions P i (i = 1, 2, ..., p) are unknown. Then, find axial induction factors a i (i = 1, 2, ..., p) maximizing P (a 1 , a 2 , ..., a p ).

Multi-Resolution Simultaneous Perturbation Stochastic Approximation
This section presents the main idea to solve Problem 2.1. First, we briefly review the standard SPSA algorithm. Based on this, we propose the MR-SPSA algorithm as a new tool to solve a class of model-free optimization problems.

Standard Simultaneous Perturbation Stochastic Approximation
SPSA is a stochastic approximation algorithm to optimize the design parameters for a pre-specified objective function. The essential feature of the SPSA is to approximate the gradient of the objective function based on the measured value of the objective function and random perturbation. Hence it does not require any explicit form of the objective function, and can be a tool for a variety of model-free optimization problems. Now, we explain the concrete algorithm of the SPSA. Consider the optimization problem given by: where f : R n → R is the objective function and θ ∈ R n is the design parameter. The SPSA algorithm [17] iteratively updates the design parameter to search a local optimal solution θ * ∈ R n of Equation (2). The update law is given by: for k = 0, 1, ..., where d(k) is the gain and g(θ(k)) is the update vector given by: . . .
In Equation (4), c(k) is another gain, and (k) is the n-dimensional random perturbation vector. For example, the gains d(k) and c(k) are given by d(k) =ã/(Ã + k + 1) α and c(k) = c/(k + 1) γ , respectively, for non-negative numbersã, c,Ã, α and γ. Meanwhile, (k) is, for example, drawn from the element-wise Bernoulli distribution: where i (k) is its i-th component. Note that, the selection of non-negative coefficientsã, c,Ã, α and γ will be performed by the guidance reported in [17]. The standard SPSA algorithm consists of the following steps: Step I: Select the non-negative coefficientsã, c,Ã, α and γ for the SPSA gain sequences d(k) =ã/(Ã + k + 1) α and c(k) = c/(k + 1) γ . Set the initial conditions of the design parameter θ(0) and set k = 0.
Step II: Generate an n-dimensional random perturbation vector (k).
In general, it is known that the standard SPSA algorithm has the advantage of solving high-dimensional optimization problem. However, the SPSA algorithm often does not offer an acceptable convergence speed. This is due to the fact that the number of experiments to achieve the convergence state is proportional to the dimension of the design parameter [18]. Consequently, a large computation time may be required for solving Problem 2.1.

Multi-Resolution Simultaneous Perturbation Stochastic Approximation
Consider again the optimization problem in Equation (2). As a solution, we propose a new SPSA algorithm called the MR-SPSA algorithm. The idea of the MR-SPSA algorithm is that the tuning process is divided into several stages with variation in the size of the design parameter's group. This is formalized as follows.
Let SPSA θ (f, θ(0)) denote the output of the above standard SPSA algorithm. Here, the termination criterion in Step VI of the standard SPSA algorithm is given by |f (θ(k + 1)) − f (θ(k))| < ε, where ε is a small number. Let q be a positive integer called the resolution step and let σ(j)(j = 1, 2, ..., q) be positive integers such that Also, we assume that for any ζ j , there exists a unique ζ j+1 such that h (j+1) (ζ j+1 ) = h (j) (ζ j ). Next, the MR-SPSA algorithm is given by: where ζ 1 (0) is given arbitrarily and ζ j (0) is a vector satisfying h (j) (ζ j (0)) = h (j−1) (ζ * j−1 ) for j = 2, 3, ..., q. The output of the MR-SPSA algorithm is denoted by ζ * q . Now, we illustrate the MR-SPSA through a simple example. Consider the design parameter θ ∈ R 9 , which are represented by the black dot in Figure 1. Here we maximize f (θ) using the MR-SPSA algorithm for q = 3. Then, the solution of the optimization problem for each resolution is described as follows: Figure 1a and ζ 1 (0) ∈ R is the given initial condition. Note that the design parameters, which are grouped in the box with a dashed line, have the same value.
Step II: Obtain the solution ζ * j in Equation (6).
Step III: If j = q, the algorithm terminates with the solution ζ * q . Otherwise, set j = j + 1 and go to Step II.
Step 2: Apply (6) to the objective function P and the design parameter a := (a 1 , a 2 , ..., a p ), i.e., by regarding P and a as f and θ, respectively.
Step 3: If j = q, the algorithm terminates with the solution a * = ζ * q . Otherwise, set j = j + 1 and go to Step 2.
In this study, we assume that the initial condition ζ 1 (0) in Step 1 is based on the existing optimal control of the standalone turbine, which produces an axial induction factor of 1/3 [16].
In the above procedure, it is important to adequately select q, σ(j) (j = 1, 2, ..., q) and h (j) (j = 1, 2, ..., q). Usually, they are chosen based on the information of the wind farm layout, wind direction and grouping policy. For example, it is performed in the following way. It would be reasonable to set q in 3 to 5 steps for the wind farm optimization problem. Note that if q is too large, it will increase the computation time while having the same maximum power production. Since most of the existing wind farm layouts locate turbines in a symmetrical way such that the distance between rows and between columns is almost identical, the suggested q (3-5 steps) is more preferable to produce reasonable computation time. Next, a policy on forming the wind turbine group is discussed. Note that the selection of σ(j) (j = 1, 2, ..., q) is based on this grouping policy. First, at each resolution j, define the selected group of the wind turbine as G jk := {j k1 , j k2 , ..., j kn k } (k = 1, 2, ..., σ(j)), where n k =| G jk |, k G jk = {1, 2, ..., n}, and G jl ∩ G jm = ∅ for l = m. Next, from the given design parameter a i (i = 1, 2, ..., n) and group selection G jk , we can construct h (j) : R σ(j) → R n (j = 1, 2, ..., q) in such a way that a j k1 = a j k2 = ... = a j kn k = ζ jk . Here, the basic policy of grouping in MR-SPSA is to cluster the wind turbines whose optimal design parameters are expected to be similar. Specifically, if one wind turbine does not have any other turbine in its downstream, we can expect its optimal parameter to be the same as the standalone turbine (which is known to be 1/3). This is because this wind turbine is allowed to produce its maximum power without considering the wake effects in the downstream. Therefore, it makes sense to form one group in which each member has no wind turbines in the downstream. Further, one reasonable way to make other groups would be to cluster the turbines, depending on the number of their downstream turbines. The detailed policy to group the turbines is stated as follows: Stage 1: Form two groups (σ(1) = 2). The turbines that may affect its nearest downstream turbines through current wind direction information are merged into one group, e.g., G 11 , while the turbines that do not have any turbines in its downstream are clustered into the other group, e.g., G 12 . Then, set the initial condition of the design parameters to be a i (0) = 1/3, which is known to be the optimal parameter of single turbine.

Stage 2:
The intermediate resolutions (1 < j < q) of MR-SPSA are applied in this second stage. Here, the group G 11 is divided into several more groups depending on the number of downstream turbines, and the group G 12 is preserved as before in this stage. Here, σ(j) can be varied between σ(2) and σ(q − 1).

Stage 3:
This stage directly corresponds to the final resolution (j = q) of the MR-SPSA-based method.
Here, the policy is to set the number of groups to be equal to the number of turbines in the wind farm. In other words, σ(q) is set to n, which is similar to the standard n-dimensional tuning problem.
In this resolution step, h (1) (ζ 1 ) = [ζ 12 ζ 12 ζ 12 ζ 12 ζ 11 ζ 11 ζ 11 ζ 12 ζ 11 ζ 11 ζ 11 ζ 12 ζ 11 ζ 11 ζ 11 ζ 12 ] ∈ R 16 . Next, for j = 2, we optimize four design parameters ζ 21 , ζ 22 , ζ 23 and ζ 24 , which correspond to the design parameter of each turbine in groups  16 . The selection of the design parameter in each resolution is illustrated in Figure 2. Remark 4.1 In this study, we choose α < 1 and γ < 1 to obtain better finite-sample performance through maintaining larger step size. Next,Ã is set to be equal to 10% of the maximum number of allowable iterations, whileã is chosen such that d(k) times the magnitude of elements in g(θ(0)) is approximately equal to the smallest value of the desired change magnitudes among the elements of θ in early iterations. Then, we appropriately choose c larger thanã to handle any fluctuation in the measurement of power production.

Remark 4.2
In this study, the objective function is limited to maximize the total power production of the wind farm in below-rated wind conditions, although our MR-SPSA-based method may be extended to a more complicated situation as long as the objective function can be measured. For example, our proposed method may also handle the objective function with mixed value of power production and fatigue load if both of them are measurable. Moreover, the MR-SPSA-based algorithm is also capable of dealing with constraint in the optimization by introducing the penalty term with the original objective function.

Remark 4.3
Since the MR-SPSA-based algorithm seeks the local optimum based on the stochastic approximation of the gradient, its performance strongly depends on the choice of initial condition. However, in the wind farm optimization case, it turns out that the MR-SPSA-based method works well if we set the initial value of the design parameter as the optimal control of the standalone turbine. Also, note that our proposed method is applicable for both smooth and non-smooth objective functions as long as the measurements of the objective function are available. Theoretically, the convergence is guaranteed only for a class of smooth function, which has been reported in [17]. In our wind farm case, the objective function may be non-smooth. However, extensive simulation results demonstrate the effectiveness of our proposed method.

Simulation Results
The MR-SPSA-based method for the wind farm optimization problem is illustrated in this section. Here, a wind farm model that represents a real commercial wind farm is used in order to simulate our model-free approach. We firstly describe the wind farm model proposed by [14], and then the MR-SPSA-based method is tested to the model of commercial offshore Horns Rev wind farm.

Wind Farm Model
The well-known wind farm model studied in several literatures [19][20][21] is based on the Park model. This model presents a characterization of wake by estimating the velocity profile of a single turbine. In order to evaluate the time efficiency of the wind farm optimization process, a delay structure was added to the Park model to illustrate the wake travelling time from one turbine to the next [14]. This delay structure is a dynamic version of the Park model and it is practical to explore the model-free approach in a real-time environment.
First, the static wake interaction between the two turbines is briefly presented. Let X = {1, 2, ..., p} be the set of p wind turbines in the wind farm, V ω be the incoming wind speed, D i be the rotor diameter of the turbine i, A j be the rotor swept area of turbine j, A ov i→j be the overlap area between the wake generated by an upstream turbine i and rotor swept area of turbine j, and φ is a roughness coefficient that represents the slope of wake expansion. Let also (x, r) be a point in the wake of the turbine, where x is the distance to the rotor disk plane of the turbine and r is the distance to the centerline of the wind turbine rotor axis. Then, the aggregate wind velocity is given by: for: where x i and x j are the distances to the rotor disk plane of the turbines i and j, respectively. The wake interaction between the two turbines is illustrated in Figure 3. Note that the aggregated wind velocity V j for j ∈ X is evaluated based on the aggregation of the wind velocity deficit created by each upstream turbine. We also assume that the diameter of the wake has a circular cross-section and expands proportionally to the distance x. Further, the power of each turbine can be represented as: Next, the dynamics of the wake interaction is illustrated based on the estimation of the wake travel time from one turbine to another, as studied in [14]. Let m(i) be the index of the nearest neighbor downstream turbine of turbine i, where it is directly influenced by turbine i. Let also R(i) = {i, m(i), m(m(i)), ...} be the set that includes turbine i and the other downstream turbines in a row that are affected by turbine i. Then, the time interval for the wake to travel to the whole wind farm can be approximated as:

Horns Rev Example
In this section, the performance of the MR-SPSA-based method is evaluated for the Horns Rev wind farm layout by using the wind farm model [14] in Section 5.1. This wind farm, which is located in the North Sea off the coast of Denmark, consists of 80 turbines (Vestas V80 2MW) of 80 m diameter. Figure 4 shows the layout of the wind farm, which is in an oblique rectangle with a turbine spacing of 7 turbine diameters, i.e., 560 m, in both the x and y directions. The roughness coefficient is φ = 0.04 and the air density is ρ = 1.225 kg/m 3 . The comparative study of the MR-SPSA, the SPSA and the existing model-free approaches for the cases of different wind directions, non-static incoming winds and several wind turbine failures are discussed in Sections 5.2.1, 5.2.2, and 5.2.3, respectively.

Performance of the MR-SPSA-Based Algorithm with Different Wind Directions
First, we assume that the wind speed is constant at V ω = 8 m/s, and study the cases when the wind direction is 170 • , 200 • , 220 • , 240 • , 250 • and 270 • . Note that the wind farm experiences a much larger wake effect when the wind direction is 220 • and 270 • [22].
In the case of wind direction 170 • , the axial induction factors of the 80 turbines are selected as design parameters (n = 80) and the time interval for the wake to travel to the whole wind farm is approximated as T w = 980 s. It is assumed that the wind farm is equipped with a single sensor to measure the total power production P after T w . Therefore, at each iteration, the sampling interval for the MR-SPSA-based method is equal to 3 T w since it requires two measurements of the objective function P to update the design parameter and another one measurement to observe the resultant P of the updated design parameter. Next, we set the parameters of the MR-SPSA-based method d(k) = 6.5 × 10 −7 /(k + 109) 0.8 , c(k) = 0.0001/(k + 1) 1/3 , q = 3 with σ(1) = 2, σ(2) = 8 and σ(3) = 80. Here, we group the turbines for each resolution according to the procedures stated in Section 4. The small number ε = 0.01 is chosen as the termination criterion in each resolution. In order to compare the proposed method with the standard SPSA-based method and the existing model-free approaches, which are the FS-MPPT-based method [14] and the GT-based method [11], the parameters for those methods are provided in this paper. Here, the gain sequences d(k) and c(k) of the standard SPSA-based method are set to be similar to the MR-SPSA-based method. Next, the FS-MPPT approach with a scaling factor for the size of the design parameter step K F = 0.0008 and the GT approach with a size of interval for random step on the design parameter K G = 0.03 and the probability of using a new random setting E = 0.3 are used. See [11,14] for the detail of the GT-based and FS-MPPT-based algorithms, respectively. In order to observe the randomization effect, we perform 100 trials for the MR-SPSA, SPSA and GT-based approaches. Figure 5 shows the response of the total power production P (a 1 , a 2 , ..., a 80 ) of each method for the first 10 h (left) and full simulation time 700 h (right). Here, the responses of the SPSA, GT, FS-MPPT and MR-SPSA are represented by the red-dash, green, blue-dot and black-thick lines, respectively. Note that the total power production response of MR-SPSA, SPSA and GT are shown from 100 trials due to the stochastic nature of the algorithms. In general, each method successfully improves the total power production during the 700 h simulation time. The statistical analysis of the final value of the total power production and the convergence time is recorded in Tables 1 and 2. Here, the convergence time is defined as the time it takes for the total power increase to exceed 90% of the final value at 700 h simulated time. Notice that, the MR-SPSA-based method achieves the highest total power production, which is 39.6056 MW. This is followed by the SPSA (39.6010 MW), GT (39.5897 MW) and FS-MPPT-based (38.3153 MW) methods. Similar pattern also holds for the mean and worst values of the total power production. However, in terms of the convergence speed, the MR-SPSA-based and FS-MPPT-based methods converge much faster than the SPSA-based and GT-based methods. This also can be observed from Figure 5 for the first 10 h of the simulation time. On the other hand, since the wake propagation in this case is distributed symmetrically, it is expected that the optimal design parameters of the turbines in the same row are same. Here, the best optimal design parameters of the MR-SPSA-based method for each row are depicted as row {1, 2, 3, 4, 5, 6, 7, 8} = {0.207, 0.162, 0.166, 0.167, 0.173, 0.116, 0.264, 0.333}, where the turbines in the first row receive the incoming wind speed V ω . It is observed that the optimal design parameters of the turbines in the final row are maintained as the initial design parameters. However, the optimal design parameters of the other turbines are decreased since they produce the wake effect to its downstream turbines. Hence, this fact supports the validity of our grouping policy, which yields a faster convergence time. Also, the optimal axial induction factors of the turbines in the first row are higher than the turbines in the intermediate rows. This result is consistent with previous study on wind farm optimization, e.g., [11].
Next, we show the results for the cases of wind direction 200 • , 220 • , 240 • , 250 • and 270 • . For each case, we have a new wind farm configuration and the parameters of each model-free algorithm need to be selected again. However, in order to observe the robustness of each algorithm to the variation of the wind farm configuration, we set the parameters of the GT, SPSA and FS-MPPT-based methods to be similar to the case of wind direction 170 • . Similarly, the gain sequences d(k) and c(k), and q for the MR-SPSA-based method are maintained as in the previous case. However, the group selected in the first and second resolutions is different for each wind direction. Hence, the function h (j) (j = 1, 2) is varied for each wind direction. Note that the function h (j) (j = q) is similar for all wind directions according to the grouping policy in Section 4. For example, the selected groups for wind direction 220 • and 270 • are shown in Figures 6 and 7, respectively. For both figures, the groups in the first and second resolutions are represented in the left and right sub-figures, respectively. For each wind direction, we also perform 100 trials for the MR-SPSA, SPSA and GT-based approaches.    According to the grouping policy in Section 4, the selection of group of turbines in Stages 1 and 3 is uniquely determined, while in Stage 2, there is some flexibility in choosing the groups. Hence, we also evaluate the performance of our MR-SPSA-based method for various group selections in Stage 2. Here, we consider the wind direction 220 • again. Next, the parameters of the MR-SPSA-based algorithm are set to be similar to the previous case with the various group selections in the second resolution. Here, four types of potential groups are considered as shown in Figure 8, e.g., H 1 , H 2 , H 3 , and H 4 . Note that the variation of group selection follows the grouping policy in Section 4. Table 3 shows the performance analysis of the total power production and the convergence time for different group selections in the second resolution for 100 trials. One can see clearly that the statistical values for the total power production and the convergence time are almost similar among selected groups. Moreover, it is also similar to the previous statistical analysis in Tables 1 and 2 for wind direction 220 • . Therefore, the performance of our MR-SPSA-based method is not very sensitive to the choice of group selection in Stage 2.  Next, we consider the incoming wind with varying speeds and directions. Figure 9 shows the response of the wind speed and direction in the top and bottom subplot, respectively, during 10 h simulation time. Note that the parameters of the MR-SPSA, SPSA, GT and FS-MPPT-based methods are set to be similar to the case in Section 5.2.1. Figure 10 shows the responses of the total power improvement by the MR-SPSA, SPSA, GT and FS-MPPT-based methods to the non-static incoming wind in Figure 9. It can be seen that the MR-SPSA-based method yields higher total power improvement than the other model-free algorithms when the incoming wind speed and direction are changed every hour. Here, the objective function is always fluctuating due to the non-static incoming wind speed and direction. Hence, the MR-SPSA always computes the first resolution of its algorithm during non-static incoming wind, which yields a rapid power increment.  Remark 5.1 Note that the performance of our proposed method is limited to slowly varying wind speed and direction conditions only.

Performance of the MR-SPSA-based Algorithm with Turbine Failures
In this section, we focus on the case where five turbines as shown in the circle in Figure 11 are not operational. In this case, a wind direction 270 • with a constant incoming wind V ω = 8 m/s is assumed. It is also assumed that the controllers are redesigned for all algorithms when the failure is detected. Since five turbines fail, only n = 75 design parameters are selected. Here, the parameter of the MR-SPSA, SPSA, GT and FS-MPPT-based algorithms are set to be similar to the case of wind direction 170 • in the previous section. The result of the total power production P (a 1 , a 2 , ..., a 75 ) of each method for the remaining 75 turbines is shown in Figure 12, while the statistical analysis of 100 trials is summarized in Table 4. From the comparative assessment, the MR-SPSA-based method still generates the highest total power production with the fastest convergence even when the wind farm configuration is changed due to the five turbine failures. Figure 12. Convergence of the total power production P (a 1 , a 2 , ..., a 75 ) with five turbine failures.

Conclusions
In this paper, a model-free approach based on Multi-Resolution Simultaneous Perturbation Stochastic Approximation (MR-SPSA) has been proposed for a wind farm optimization problem. The method tries to achieve fast convergence of the optimization process even for a large scale wind farm by taking account of the wind farm configuration and the wind directions. The proposed method is tested through extensive simulation using the Horns Rev wind farm layout, which takes account of wake aerodynamic interactions in the wind farm. The simulation results illustrate that the proposed method outperforms other existing model-free approaches from the viewpoints of maximum power production and convergence speed for various wind directions. Furthermore, the simulation has been performed in the presence of wind turbine failures, slow wind changes and communication disruptions, and the MR-SPSA-based approach is shown to be robust in such cases. p Number of wind turbines in the wind farm a i Control parameter of turbine i a Control parameter vector of the turbines in the wind farm P i Power production of turbine i P Total power production of the wind farm n Number of design parameters f Objective function for R n → R θ, θ * Design parameter vector a, c,Ã, α, γ Nonnegative coefficients of SPSA parameter (k) Random perturbation vector at k-th step i (k) i-th component of (k) ε Small number q Resolution step σ(j) (j = 1, 2, ..., q) Number of design parameters at each resolution j h (j) (j = 1, 2, ..., q) Function for R σ(j) → R n at each resolution j ζ j , ζ * j (j = 1, 2, ..., q) Design parameter vector at each resolution j ζ jk (j = 1, 2, ..., q) k-th component of ζ j G jk (k = 1, 2, ..., σ(j)) Groups of the wind turbine at each resolution j n k Number of turbines in the group X Set of p wind turbines in the wind farm V ω Incoming wind speed D i Rotor diameter of the turbine i A j Rotor swept area of turbine j A ov i→j Overlap area between the wake generated by an upstream turbine i and A j φ Roughness coefficient x Distance to the rotor disk plane of the turbine r Distance to the centerline of the wind turbine rotor axis V j Aggregate wind velocity ρ Air density m(i) (i = 1, 2, ..., p) Index of the nearest neighbor downstream turbine of turbine i R(i) (i = 1, 2, ..., p) Set of turbine i and other downstream turbines in a row that are affected by turbine i T w Time interval for the wake to travel to the whole wind farm K F Scaling factor for the size of the design parameter step in the FS-MPPT-based method K G Size of interval for random step on the design parameter in the GT-based method E Probability of using a new random setting in the GT-based method