Wake-effect aware optimal online control of wind farms: An explicit solution

Wake effects impose signiﬁcant aerodynamic interactions among wind turbines. To improve the wind farm operating performance, practical wind farm online control considering wake effects becomes very important. To achieve online optimal wind farm control while responding to grid demands, this paper proposes a novel optimal wind farm supervisory control (SC) model and its explicit solutions. From the controller modelling perspective, the two major wind farm operating modes, the maximum power point tracking mode and the set-point tracking mode, are ﬁrst analysed and uniﬁed in one optimisation model while considering wake effects. In this way, wind farm power production and rotor kinetic energy reserve can be simultaneously considered to conveniently modify the operation mode in response to different grid demands. Aside from controller modelling, the collocation method is ﬁrst introduced to address the online application problem of such wake-effect aware optimal WF control. Although a few optimisation algorithms have been proposed to ﬁnd the optimum ofﬂine, online optimal control is still challenging because of the computational complexity brought by wake model non-linearity and non-convexity. The proposed collocation method explicitly approximates the optimal solutions to the proposed supervisory control model, through which only a direct algebraic operation is required for online optimal control instead of repeated optimisations. Case studies are carried out on different wind farms under various wind conditions, showing that the wind farm power production potential and releasable power reserve are improved compared to traditional greedy control in both modes. The accuracy of the collocation method is veriﬁed. A detailed analysis of the wind farm production capacity under different wind speeds and directions is also provided.


INTRODUCTION
Alleviating the wake effect and satisfactorily integrating wind farms (WFs) into the electricity grid are two considerable challenges in the science of wind energy [1]. To address these two challenges, practical WF controllers that can flexibly modify the WF operation modes in response to the grid demands while considering wake effects are of great importance [2]. The typical WF control framework has a hierarchical structure composed of wind turbine (WT) level control and WF level control [3], [4]. Farm-level control, also called WF supervisory control (WF SC) [3], [5], derives the operating points for each WT to deal with various ambient wind conditions while responding to the grid demands. Turbine level control focuses on the operation of a single WT to ensure that references from the SC are reached despite local wind turbulence. Single turbine control has been investigated extensively [6]. However, to handle the wake effect while reacting to different grid demands, controller modelling and online application have not been fully addressed for optimal WF level control.
Optimal WF level control is challenging because of the aerodynamic interactions among WTs via wake effects [7]. The operation of a WT causes a reduced wind speed in downstream wind flows, which is known as the wake effect [8]. The wake effect imposes significant aerodynamic interactions among WTs [9]. In a WF, downstream WTs often operate in the wake of upstream WTs, resulting in lower power production Therefore, coordinated WF control considering the wake effect is essential to increase the overall performance of a WF [9] [10]. Compared with traditional greedy control, where WTs are optimised individually [11], coordinated WF control has recently attracted greater attention.
From the controller modelling perspective, the existing WF control studies mainly focus on two major operation modes separately: the maximum power point tracking (MPPT) mode and the set-point tracking mode. For the MPPT mode, the WF tries to maximise its power production. Reference [12] found that by reducing upstream WT production, it is possible to increase the overall WF production. Reference [13] coordinately optimises WT axis induction factors and successfully increased WF production compared with the greedy control. Reference [14] further extended this model and treated the WF power maximisation problem while accounting for the WT power constraints. However, to simplify the model, the axis induction factor is used as the decision variable. As a result, the WT rotor speed, which stores important kinetic energy reserves, cannot be considered in their models. Apart from power production, power grid frequency regulation is already a requirement of WFs in some jurisdictions to increase the penetration of wind generation in power systems [15]. To this end, the set-point tracking mode was developed, which manipulates rotor kinetic energy. Through simulations and measurement data, [16] found that the wake effect has a significant influence on the WF inertial response capacity, which heavily depends on the overall rotor speeds of WTs in the WF. References [10], [17] proposed an optimal WF SC model for the set-point tracking mode, which explicitly integrating the wake equations in rotor kinetic energy optimisation. However, the wake model in their studies is a simplified one-dimensional model. Only one-column small WFs under one wind direction could be studied.
Generally, from the controller modelling perspective, although the existing literature has investigated both WF operating modes separately, the combination and coordination of the two modes have not yet been addressed. However, to provide different services in response to grid demands, both power production and rotor kinetic energy should be simultaneously considered in the WF supervisory control model [16], [18]. To this end, a unified control model is desired to flexibly modify the operating modes in response to the grid demands.
From the application perspective, the nonlinearity and nonconvexity of the wake effect increase the difficulty of solving such optimal control models integrating wake equations. Some methods have been used and achieved great performance in finding the optimum for off-line calculation. Reference [9], [13] used gradient-based programming methods. Sequential convex programming was used in [9] to maximise the WF power production and obtain quality results. Reference [19] used dynamic programming method. Model-free methods such as game theory [20], adaptive filtering [21], Bayesian ascent [22] and extremum seeking [23] have also been used and achieved good results. However, these optimisation algorithms are not suitable for online applications, which require fast reactions and low computational complexity. To obtain an optimum under one certain wind condition, gradient-based methods need to carefully select the initial point and many iterations [13], while the model-free methods suffer from slow convergence and sensitivity to uncertainties [20]. With the expansion of the farm size, the computational difficulty is increasingly severe [24]. Consequently, to obtain the optimal results in real time addressing various wind conditions, repeated optimisation is not acceptable for large WF online control.
From the application perspective, although some optimisation algorithms have been proposed to solve the WF control problem offline, their computational complexity is still one of the major issues for online control, especially for large WFs.
In this paper, we propose a novel optimal WF SC model and its explicit solutions. From the modelling perspective, both the MPPT and set-point tracking modes are unified into a consociated model to enable flexible switching of operating modes in response to the grid demands. The WT operating limits are also constrained. A two-dimensional wake model is utilised so that large WFs can be analysed under different wind directions. From the application perspective, the collocation method (CM) is introduced [25], [26]. By explicitly mapping the wind condition to the optimal control laws, CM avoids repeatedly solving the computationally expensive optimization problem and gives explicit approximations of optimal control laws of the whole WF. The contributions of this paper are twofold.
(1) A unified optimal WF supervisory control (WF SC) scheme considering the wake effect is proposed. The two distinct WF operating modes, MPPT and set-point tracking mode, are first optimised in one unified optimization model. In this way, both power production and rotor kinetic energy are simultaneously considered and coordinated and the WF can optimise the operating modes in response to different grid demands. Compared with traditional control schemes that do not consider the wake effect, increased WF power production potential and releasable power reserves can be achieved.
(2) CM is first introduced to the WF control field to enable online application. The non-convex WF SC optimisation model is solved offline under only a few specifically selected wind conditions. Then the optimal control signals for different wind conditions can be approximated by the CM. In this way, only direct explicit algebraic operation is required for the online optimal WF control, alleviating the online computational burden. The data efficiency and accuracy of the proposed control methods are also demonstrated.
The remainder of this paper is organised as follows. Section 2 briefly introduces the WT mechanical model and the wake model, formulating the detailed WF model considering WT aerodynamic interactions. Section 3 constructs the proposed optimal WF SC model. Both the MPPT mode and setpoint tracking mode are analysed. Section 4 introduces the CM and its application in the WF control framework. Case studies and conclusions are presented in Sections 5 and 6.

QUASI-STATIONARY WF MODEL
In this section, we introduce the detailed quasi-stationary WF model, which formulates the operating constraints for the optimal control model in Section 3.

WT model
According to aerodynamic theory, the mechanical power generated by WTi can be computed as [27] where subscript i is the index for the ith WT, is the air density, and A i is the blade swept area. V i is the upstream wind speed. In a WF, V i is determined by the control variables of all the upstream WTs, and we will formulate it in Section 2.2. C p is the power coefficient, which is a function of the blade pitch angle and the tip speed ratio . can be represented by where i is the rotor speed and R is the rotor ratio. C p is determined by the aerodynamic design of the WT blades, usually approximated by curve fitting with the manufacturer data [10], and can be written as the polynomial form [28]. The kinetic energy stored in the rotating mass of WTi can be expressed as where J is the moment of inertia of the WT rotating body.

Wake model
The wake model used in this work, the Jensen model, is the most prevalent analytical wake model. Reference [29] compared different wake models and found that the Jensen model provides the best fit with realistic WF measurements. Reference [16] validated the effectiveness of the simulation based on the Jensen model against the Horn rev measurements. Unlike the one-dimensional model in [10], the Jensen model is a good tool for analysing large WFs with multiple rows and columns and has been widely used in various WF control studies [13], [16], [19], [21]. As shown in Figure 1, the wake behind a WT is assumed to be linearly expanded, and the area outside the dashed lines is not affected. In a WF, the wind speed in front of a downstream WT i can be formulated as where U is the ambient free stream wind speed. v i is the aggregated deficit factor quantifying the overall speed reduction caused by all upstream wakes. A single wake deficit produced by one upstream WT j is calculated as where D is the blade diameter. is the decay parameter. It should be calibrated with the measurements or high-fidelity data to modify the model error for different WFs [16]. However, this is outside of the scope of this work. x i , x j are the WT i and WT j position coordinates in the wind direction. C t ( j , j ) is the thrust coefficient of upstream WT j , which is also a non-linear function of and and is usually approximated by curve fitting with the manufacturer data [10]. Based on Equation (5), the aggregated wind deficit considering all the wakes of upstream WTs is given by is the overlapping area of the wake corridor of WT j and the blade disk plane of WT i , as shown in Figure 1. Φ i denotes the set of upstream WTs of WT i . Φ i and Φ i denote the control variables of all WTs in Φ i . Because of the wake effect, V i in front of each WT is influenced by the control variables of all upstream WTs in Φ i .
Note that Φ i is jointly determined by the WF layout and the wind direction . Take a 4-WT small WF as an example. As shown in Figure 2 is empty and the four WTs operate separately. Therefore, Φ i is a discrete function of and determined by the WF layout. We denote it as Φ i ( ).
Then, the detailed WF power model incorporating the wake interactions can be modelled in the following steps.
Step 1: Given the WT type, the C t and C p curve fitting in terms of and can be determined [10].
Step 2: Given the position coordinates of all in-service WTs in the WF and the wind direction , Φ i ( ) for each WTi are simple to determine.

FIGURE 3 Overall WF control framework
Step 3: V i ( Φ i , Φ i ) in front of each WT can then be modelled through the Jensen model (4)-(6).
Step 4: By summing the power of each turbine, the overall WF power production P WF can be formulated as

OPTIMAL WF SUPERVISORY CONTROL
As shown in Figure 3, the typical task of the optimal SC is to find the optimal operating point of each WT in response to the grid demands under various wind conditions. From the grid perspective, the main WF operating modes are the MPPT mode and set-point tracking mode [3]. In MPPT, the WF tries to maximise its power production. On the other hand, with increasing wind energy penetration, WF is preferred for operating as a controllable entity like conventional power plants and for providing frequency support [30], [31]. To this end, set-point tracking control is proposed to constrain the power output to a specified reference. The two distinct modes were investigated separately. However, to provide different services in response to grid demands, both modes should be coordinated in WF supervisory control. To this end, a unified control model is proposed that simultaneously optimises WF power production and rotor kinetic energy so that the WF can flexibly modify the operating modes in response to different grid demands. To optimise the quasi-stationary operating point of the entire WF, the proposed WF SC model unifies both operating modes and considers the WT operating constraints. The quasi-stationary optimisation model of an N -turbine WF is formulated as follows: WT model : Wake interactions : = {0, 1} is a predefined mode parameter in response to the grid demand. = 0 when the WF is required by the power system to operate in the MPPT mode, while = 1 operates in the set-point tracking mode. We introduce to unify the two operating modes in one form. It is pre-determined in response to the grid demand before the optimisation is solved. E rot WF = ∑ N i=1 E i and denotes the total kinetic energy of the WF stored in the rotating mass [10]. E i is defined in Equation (3). E 0 is the predefined WF rotor kinetic energy level, which is also determined by the power system [32]. P WF represents the total WF power generation defined in Equation (7). P max WF (U, ) denotes the maximum WF power generation potential under the specific wind condition [U, ], which is a constant. The grid active power demand is P dem = (1 − )P max WF . is the WF deloading proportion. , P dem and E 0 are all predefined parameters in response to the grid demand. i is the index for the ith WT in the WF.
The rotor speed and pitch angle of individual WTs are constrained in Equations (11) and (12), respectively. Constraint (13) limits the WT output power by the rated capacity P rate . Given the wind condition [U, ] and the farm layout, Equations (14) and (15) formulate the detailed quasi-stationary WF model.
1) = 0, WF maximum power production mode (WF MPPT) MPPT is one of the main issues for WF operation. When is set to 0, the objective function f (U, ) = P WF , so the objective function (8) maximises the total WF power generation. Constraint (10) is inactive. Constraint (9) is used to maintain E rot WF above the predefined rotor kinetic energy level E 0 . Given the wind condition [U, ], P max WF (U, ) is the optimal value of the objective function obtained by solving this mode.
2) = 1, set-point tracking mode This problem focuses on improving the amount of additional power a WF can provide during the frequency excursion and therefore determines the WF frequency support capability. Therefore, in this mode, the optimisation model maximises the releasable kinetic energy in the quasi-stationary state in the WF level. When power production is required to maintain a certain level by the power system, the WF should operate in this mode. The objective function f (U, ) = E rot WF , witch maximises the releasable kinetic energy in the quasi-stationary state in the WF level. Constraint (9) is inactive. Constraint (10) enforces the WF power output following the grid demand and maintains a power reserve level . In this mode, the power reserve and the frequency regulation services can be provided by the WF [10].
According to the grid demands, the WF is deloaded by a predefined margin to provide the power reserve and follow P dem according to the constraint (10). Additionally, by exploiting wake effects, this mode provides a coordinated quasi-stationary distribution of P dem to individual WTs. By exploiting the wake effect and coordinating WTs, the proposed model can possibly increase the WF kinetic energy without reducing the total output power. In this way, the WF frequency regulation capability is improved without any additional cost. We call this distribution scheme coordinated distribution (CD).

Application difficulties of optimisation algorithms in online WF control
For WF control, f (U, ) and the resolved optimal decision variables [ opt , opt ] need to be determined for various wind conditions covering the main operating range. Because of the wake model (4)-(6), the SC optimisation model is non-linear and non-convex. For offline calculations, some methods have been proposed and achieved good performance in finding quality solutions [9], [13], [22]. However, these algorithms are not suitable for online application, which requires fast reactions and low computational complexity. For one optimization under one wind condition, model-free algorithms such as game theory methods suffer from slow convergence and sensitivity to uncertainties [20]. Additionally, the model-based methods such as sequential convex programming (SCP) [9] require many iterations. Because of the non-convexity, the results obtained by repeated optimisations are influenced by the initial point. Although we can modify the solving procedure and try different initial points to find the optimum offline, the complexity is not acceptable for online control applications. With the expansion of the farm size, repeatedly solving the non-convex optimisation problem for various wind conditions is computationally expensive for online control.
Therefore, the lack of explicit control laws hinders the application of optimal coordinated WF control in online control. To this end, we first introduce the collocation method (CM) to the WF control field.

4.2
Explicit WF control based on the collocation method CM is a polynomial approximation method with high data efficiency. The power of the CM lies in its ability to select appropri-ate interpolation points to create a polynomial model that has a higher order approximation accuracy. The optimisation model only needs to be solved under a few specially selected wind conditions offline. Then, explicit approximations of optimal control laws can be derived by the proposed CM. For the online optimal WF control problem, only a direct explicit algebraic operation is required instead of repeated optimisations.
To clarify the theory and application procedure of the proposed CM, we first take the MPPT mode as an example. The theoretical procedure for both modes is the same. In the MPPT mode, with the variation of the wind condition [U, ], the resolved optimal decision variables ≜ [ opt , opt ] and the optimal objective values of the model (8)-(15) also change. Thus, the optimal control laws with respect to different wind conditions [U, ] are defined implicitly through the optimisation model (8)- (15). Our target is to formulate an explicit approximation to relate [U, ] to the optimal control laws with as few interpolation points as possible since each optimisation is computationally expensive. The CM-based approximation of the optimal control law can be expressed as a linear combination of basis functions * where * represents the approximation, j (U, ) is the j th basis of the basis function set { j (U, )}  j =1 , and c k j is the coefficient of the j th basis.  is the basis size. Index k means the kth control variable. For the set-point tracking mode, the parameter space extends to [U, , ]. The WF deloading proportion is also included to modify the WF operation according to the required power level. Given the approximation order n, the basis function set is constructed in the tensor product form as where P belongs to a particular orthogonal polynomial set. According to approximation theory [25], [33] and our previous works [34], the Legendre orthogonal polynomial series provides good approximation accuracy. For the set-point tracking mode, the basis function set is { j (U, , )}  j =1 accordingly. We choose to construct the basis using Legendre orthogonal polynomials here. Readers who are interested in more details of orthogonal polynomials are referred to [25], [26]. P m 1 (U ) denotes the m 1 th degree Legendre polynomial in U formulated as where P m 2 ( ) is the m 2 th degree Legendre polynomial in . ( denotes the combination operation. We omit the proofs, but it can be proven that each P m 1 has exactly m 1 roots [26]. Now, the task is to identify the  coefficients in Equation (16). According to Equation (17),  = n 2 for the MPPT mode, n 2 coefficients in Equation (16) need to be determined. For the set point tracking mode, there are n 3 coefficients because the parameter space is three dimensional. This can be done by solving the optimisation model (8)-(15) many times under many different wind conditions [U, ] and applying a least squares algorithm. However, we are attempting to solve the optimisation model (8)-(15) under as few wind conditions as possible because each optimisation is computationally complex and requires carefully selecting the initial point to obtain the optimal results. The power of CM lies in its ability to judiciously select appropriate situation points to create the polynomial that has a higher order accuracy. To this end, CM judiciously interpolates the polynomial with the optimisation results at [U l , l ] as follows: The specially selected  wind conditions {[U l , l ]}  l =1 , which are the roots of the equation P n (U )P n ( ) = 0, are determined in Equation (19). Because each P n has exactly n roots [26], n 2 wind conditions are selected. This selection algorithm is based on the Gaussian quadrature integration (GQI) theory. For the set point tracking mode, accordingly, the selected interpolation points also include : ] | P n (U )P n ( )P n ( ) = 0}. According to GQI, if the polynomial is interpolated using exactly the right points, the expected value of an (n-1)th degree approximation polynomial is identical to the expected value of any polynomial of degree less or equal to the (2n-1)th degree [26]. Similarly, higher order moments are also well approximated. This means that the polynomial interpolation can obtain a higher order accuracy based on the CM and the points selected by Equation (19) ensures a high data efficiency. The math proof will be further discussed in Section 4.3.
Therefore, the optimisation model (8)-(15) is solved offline under only the specially selected  wind conditions by Equation (19). Then, the following linear equations can be solved to determine the  coefficients c k j where k (U, ) are the optimisation results of the optimal SC model (8)- (15) for the kth control variable at the chosen wind condition [U, ] in {[U l , l ]}  l =1 . k (U, ) can be derived by solving the optimisation model offline [9], [13], [22], [23]. In this literature, we use sequential convex programming (SCP) and carefully adjust the initial point for different wind conditions. The quality of the offline solutions by SCP has been demonstrated in [9] and is shown in Section 5. Readers interested in the details of SCP, which is a commonly used method for non-convex optimisation and has been used to solve such optimisations integrating wake equations, are referred to [9]. For the set-point tracking mode, the parameter space extends to [U, , ]. The theoretical procedure for both modes is all the same.

Discussion on the accuracy of explicit control and its implementation procedure
The CM-based approximation polynomial can obtain the same moments as a higher order model based on the selection of suitable points by Equation (19). The proof is based on orthogonal polynomials and GQI [26] and is provided here. Let w(x) be the probability density function describing the distribution of a system parameter x. A set of polynomials { j }  j =1 is said to be orthonormal if and only if the following relationship holds for all P n in { j }  j =1 : where P n is an orthogonal polynomial in this set { j }  j =1 of order n.
Theorem: For any polynomials of degrees up to 2n − 1, their expected values can be accurately identified by n samples.
Proof. Let ℙ 2n−1 denote all polynomials of (2n − 1) degrees and let {x l } n l =1 be the zeros of n degree orthogonal polynomial P n . Let f * = ∑ n−1 j =0 c j P j (x) be the (n-1) degree interpolation polynomial constructed by the sum of orthogonal polynomials, and c j is the constant coefficient. Then, any f ∈ ℙ 2n−1 can be expressed as where r (x) is a polynomial of degree of at most (n-1) [25]. By orthogonality of P n in { j }  j =1 , the expected value of f (i.e. ∫ fwdx) can be easily derived To identify c 0 , we need to choose interpolation nodes and operate through Equation (22). Only n samples are sufficient if we use the zeros {x l } n l =1 of n degree orthogonal polynomial P n as interpolation nodes because P n (x l ) = 0 in Equation (22). Otherwise, more interpolation points are required to identify r (x). Therefore, only n samples are needed to compute the expected value of the (2n − 1) degree polynomial f . ■ Algorithm 1 Application procedure of CM-based explicit optimal WF control (Take MPPT mode as an example) Require: Approximation order n; power system demands: .
Ensure: Explicit approximation of optimal WF control laws under the input approximation order and power system demands for all wind conditions: Equation ( 4: Solve the optimisation models and store the offline solutions. In this paper, we use SCP [9]. Other proven optimization algorithms can also be used.

5: Formulate (20) based on the optimization results of
Step 4. Solve (20) directly by matrix inversion to obtain the coefficients c k j in Equation (16). The approximation polynomial is derived.
Higher order moments are also well approximated in this way, although a firm theoretical basis is lacking [25]. Readers who are interested in more details of GQI theory are referred to previous studies [25], [26]. Section 5 presents numerical experiments to further demonstrate the accuracy and the data efficiency of the proposed CM-based explicit WF control. This proof gives the math theories of CM advantages over normal interpolation methods. The overall implementation procedure is summarised in Algorithm 1.

Discussion on the polynomial order selection
Typically, the approximation order n is a hyper-parameter of the proposed polynomial interpolation method. Just like hyperparameters in the normal interpolation or approximation methods [26], [34], the higher order n typically results in a higher order accuracy, but more sample points at the same time. So the polynomial order should be chosen based on the required approximation accuracy. For the explicit WF control application in this paper, the parameters space is only two dimensional for the MPPT mode and three dimensional for the set-point tracking mode. The sample space is small. So the hyper-parameter can be chosen based on the standard hyper-parameter choice method, i.e. try-and-error experience to obtain required approximation accuracy [26], [34].
For future WF researches with more complicated grid demands and higher-dimensional parameter space, Smolyak adaptive sparse algorithm can be used alternatively to incrementally add higher order terms to the polynomial in an adaptive way until a given error tolerance is reached [35], [36]. The basic idea is to sequentially add decoupled terms, low-order coupled terms and high-order coupled terms into the polynomial. Each time a new term is added, the polynomial is updated and compared to the error tolerance to determine whether a new term needs to be further added or not. The impact of a higher order term on the polynomial interpolation precision is assessed by the its coefficients. Only the terms with a noticeable impact are added. Since this algorithm is not a contribution nor the application focus of this paper, to save space, we do not discuss its equation details further. Instead, details of the procedure can be found in [35].
Through CM, there are fewer required sample points compared to those of the common methods such as least squares algorithms [26] or lookup tables. Additionally, the choice of sample points is independent of the target control viable. One set of optimisation results is sufficient for all the interest outputs of the model. Therefore, CM is very suitable for WF control. The approximation accuracy and the online control performance are illustrated in Section 5. Via CM, the online optimal WF control problem is switched from repeated optimisation (8)- (15) to explicit algebraic operation (16), which is suitable for online control.

CASE STUDIES
In this section, two examples are presented: a 25-WT WF and an 80-WT WF that replicates the layout of Horns Rev1 in Denmark. The WF simulation is carried out based on Aeolus SimWindFarm [37], which is a widely used fast WF simulation toolbox based on the MATLAB/Simulink environment. WTs are the commonly used NREL 5 MW Type III WT [11], whose main parameters are provided in Appendix A. Other WT types are also applicable.

25-Turbine example
A WF with 25 WTs is first studied, as shown in Figure 4. The WTs are placed in a grid pattern with a distance of 5D. 1) The influence of U and on P max WF For the MPPT mode, the main objective is to maximise the overall power production of the WF. Therefore, the performance measure of the proposed control scheme is the power production improvement. Figures 5 and 7 show the total WF production of 25 turbines derived by the proposed WF SC in WF MPPT mode and traditional single WT MPPT under different wind conditions. Traditional WT MPPT optimises a single WT operation individually [6], and the algorithm is mature. The power production of the proposed WF MPPT is always higher than that with greedy WT MPPT. The improvement is plotted with the percentage value on the right. Figure 5 illustrates the impact of wind direction . We only consider the wind directions ranging from 0 o to 45 o because of the symmetry of the WF layout. The wind speed is set to 9 m/s. When = 0 o , the wake effect is the most severe due to the alignment between the wind direction and the WT array. Therefore, the power production is the lowest. In this case, the proposed WF MPPT can improve the power output by 20.45%. Traditional WT MPPT optimises a single WT operation individually [6], and the algorithm is mature. By considering the aerodynamic interactions, the proposed WF MPPT significantly alleviates the wake effect. Along the direction = 45 o , the relative distance between WTs is larger, and the wake interaction is weaker. When = 13 o and 34 o , the wake effect is minimum since WTs are staggered in this direction. As expected, the power production is highest, and the improvement of the proposed WF MPPT compared to WT MPPT is only 1.32%. The wind direction affects the relative locations of WTs and thus the FIGURE 7 Improvement of WF total power production under various wind speeds FIGURE 8 WF total kinetic energy for different wind speeds and different wake interaction patterns. The wind speed throughout the entire WF under the three directions is depicted in Figure 6. Figure 7 investigated the influence of the free-stream wind speed U on P max WF . After the cut-in speed, the wind speed range can be divided into three sections based on its influence on WF power production. As shown in Figure 7, when U is lower than 10.5m∕s, all WTs in the WF are partially loaded [11], in which and are jointly optimised and Equation (13) is inactive. In this wind speed section, the total power production increases with the wind speed. However, the improvement percentage of the proposed WF MPPT compared with WT MPPT remains at the same level. This is because, in this section, the wake interaction pattern among WTs is independent of U in the MPPT mode [9]. With the increase in U , the upstream WTs are fully loaded [11]. Equation (13) is constrained for these WTs. However, the downstream WTs are still operating in partially loaded regions due to the wind speed deficit. Thus, in this transition region, the improvement percentage decreases because the wind power is spilled for upstream WTs. When all the WTs have reached the rated output as U further increases, we say that the WF reaches its rated generation. The exact wind speed does not affect the WF power output any more before the cut-out speed.
2) Set-point tracking mode In the set-point tracking mode, the main objective is to maximise the releasable kinetic energy. The performance measure is the releasable kinetic energy improvement. Figure 8 gives an overview of the kinetic energy E rot WF under different deloading levels and wind speeds U . To highlight the effectiveness of the proposed method, is 13 o in this case in which direction the wake effect is minimal. As shown in Figure 8, with increasing U , E rot WF also increases until all the WTs reach the rotor speed upper limit max , which is 1.2671 rad/s [11]. In addition, the decrease in will reduce the WF holistic rotor speed and consequently E rot WF . Take U = 6.5m∕s as an example. For = 30% and 25%, E rot WF reaches the maximum WF rotor kinetic energy limited by the rotor speed upper limit. For = 5%, only 56.4% of the maximum is reached. This is because of the nature of the WT C p ( , ) curve [27], which gives rise to this trade-off between E rot WF and P WF . By considering the wake effect, the proposed WF SC in setpoint tracking mode, i.e. coordinated distribution (CD), successfully improves the releasable energy without reducing the total output power compared to the commonly used proportional distribution (PD) [5] [38]. The improvement by CD is shown  Figure 9. The currently commonly-used PD distributes the power demand proportionally to the available power capabilities of each WT, which do not consider the wake effect. We will further introduce it in Section 5.2. The wind conditions of 13 o and 8m∕s are used as an example. With decreasing , the kinetic energy derived by PD decreases much more significantly than that of CD, and the energy gap between the two methods increases. When = 5%, CD can increase E rot WF by 65.03% compared with PD. For other wind directions when the wake is more severe, such improvement is more distinct.
By exploiting the wake effect and coordinating WTs, the proposed WF SC model increases the WF kinetic energy without reducing the total output power. In this way, the WF frequency regulation capability is improved without any additional cost.
3) Explicit polynomial approximations of optimal control laws Figures 10, 11, and 12 compare the results derived by explicit algebraic operation through approximation polynomi- For CM, we only select 100 wind conditions according to the GQI in Equation (19) for optimisations offline. Namely, n in Equations (17) and (19) is set to 10. Then, the polynomials Equation (16) are derived to approximate the relationship from and U to P max WF and the optimal control laws of all WTs. Note that the same 100 points are utilised to determine the polynomial models for all the outputs of interest.
Specifically, the approximate optimal control law for the of WT 1 1 as circled in Figure 4 takes the following form. We omit 88 items to save space. To verify the accuracy of derived explicit approximations, the optimization model (8) As shown in Figures 10, 11 and 12, the explicit approximations are in great agreement with the optimisation results throughout all the wind conditions. For P max WF ( , U ) in Figure 10, the mean percentage error between the explicit approximations and optimisation results in the 1 729 points is only 1.34%; 95.2 percent of the gaps are less than 5%. Figures 11 and 12 illustrate the derived and for WT 1 1 . For , the mean percentage error of the explicit control laws is only 1.03% with a maximum deviation of 6.13%. For optimal , the mean error of the explicit approximations is only 0.09 degrees.
As shown in Figure 12, when WT 1 1 is upstream of other WTs (for example, < 14 o ), the optimal reference value of increases obviously compared to the wind direction where its wake does not influence other WTs and thus it operates alone as WT MPPT. These results reveal that the deloading of the upstream WTs is performed by the WF MPPT for the overall WF power maximisation.
To illustrate the data efficiency and accuracy of the proposed CM, we also compare CM to the standard look up table (LUT) method based on linear interpolation. With the same number of 100 sample points, where each dimension of [U, ] is equally discrete by ten points, the mean percentage error of LUT for CM is also applied in the set-point tracking mode to further verify the advantages. Table 1 shows the accuracy of LUT to approximate the optimal control law of for WT 1 1 under the set-point tracking mode. The sample efficiency of LUT limits its application when the parameter space dimension increases. For the proposed CM, with the help of smolyak to adaptively determine the added high-order terms, only 294 points are required to obtain the accuracy of 2.43%, which are much less than the points required by the LUT.
For the optimisation problem under one wind condition, the SCP takes approximately 24 minutes on average with greedy control as the initial point in a 3.6 GHz Intel Core i7 processor in MATLAB R2016a, which is apparently not acceptable for real-time online control. However, with CM, the only necessary optimisations are the few to identify the polynomials offline. In real-time online control, only an explicit algebraic operation through the derived explicit polynomials is needed to obtain the optimal control signals for different wind conditions, which only takes milliseconds to compute. The online computational complexity is reduced significantly in this way, making the proposed WF SC feasible for large WF online control applications.

Horn rev WF
In this section, we proceed to a realistic WF that replicates the commonly employed offshore WF Horns Rev1 located in Denmark, as shown in Figure 13. Eighty WTs are placed in the corners of a parallelogram layout with a 7.2 o tilt, whose side lengths are approximately 7D. In this case, we illustrate the contribution of the proposed explicit online optimal WF SC to WF perfor-  Figures 14 and 15 present the performance of the WF by the widely used PD [38] and the proposed explicit optimal control, respectively. P dem is the grid power demand. P out is the actual WF power generation. P avail is the maximum available wind power of the entire WF V i eff (t ) is the effective wind speed in front of WT i at time t . Many estimation techniques for V i eff (t ) have been developed [39], but they are outside the scope of this paper. C max p is the maximum power coefficient, which is a constant for a certain WT type. Traditionally, the widely used PD allocates P dem to WTs as follows, which do not consider the wake effect [5] P WTi = P dem P avail P WTi avail .
P WTi is the power command to WT i . When the set points of WTs change, the wind field will change accordingly and in turn affects the WF operation. If the WF control can not modify the operation modes while considering the wake effect, the generation capacity may be evaluated incorrectly and the control signals may not be persistently acceptable in the steady state. As shown in Figure 14, at t = 900 s, the grid demand is increased to 320 MW. According to the current wind speed estimation, P avail is sufficient to support P dem . However, when the reduced wind speed caused by up-WTs successively reaches the down-WTs after a time delay, P avail decreases gradually. As a result, the WF cannot provide P dem after 1 060 s.
However, because the proposed WF SC considers the wake interactions of the whole WF, the WF generation capacity can be correctly evaluated. In addition, because the unified model can modify the operation modes conveniently in response to different grid demands, the power demand signal is well-followed by the proposed methods. Figure 15 shows the WF controlled by the proposed SC with the CM explicit polynomials. P dem is the same as that in Figure 14. By operating the WF as a whole, the wake effect is alleviated, and the wind speed experienced by the downstream WTs is increased; thus, the maximum steady-state P out and P avail are improved. In addition, because the WF capacity is correctly evaluated, we do not need to worry about the failure of following P dem caused by the wake effect.

CONCLUSION
This paper proposes an wake-effect aware online optimal WF SC model. Two major WF operation modes are first optimised in one unified model so that the WF operation can be flexibly modified in response to different grid demands. The WF optimal operating points in the quasi-stationary state are calculated for both MPPT and set-point tracking modes. Through the proposed control scheme, the WF power production potential and releasable power reserve are improved significantly compared to those of traditional single WT control. In addition, from the application perspective, CM is first introduced in the WF control field. By solving the WF SC model offline under only a few specifically selected wind conditions, the optimal SC results can be explicitly approximated by the proposed CM. For the online optimal WF control problem, only a direct explicit algebraic operation is required instead of repeated optimisations, making the proposed WF SC feasible for large WF online control. For future work, an important direction is the modelling and control of the WF in the dynamic state considering the wake effect time delay.