Optimal power sharing of wind farms for frequency response

This paper presents a uniform optimal power sharing strategy to coordinate the wind turbines (WTs) in a wind farm (WF) to provide occasional and continuous frequency response (FR). The coordination of WTs is formulated as an optimisation problem that takes into account the WT dynamics and tries to reduce the long term loss of energy yield caused by the provision of FR. This is achieved by maximising the total kinetic energy of the WF over time while reducing wear and tear of WTs. The proposed optimal power sharing strategy relies on periodic communication between each WT and a WF controller. Local linear approximations are employed to predict the system behaviour and the solution of the optimisation problem is obtained using the proposed centralised and/or distributed algorithm. The distributed algorithm only requires one-way communication between the WF controller and local WTs, reducing the communication overheads. Simulation studies are carried out on a WF model to demonstrate the effectiveness of the proposed strategy. The results show the strategy enables reduction of yield loss over previous methods while avoiding over torque operation during FR provision.

their current operating point. By coordinating individual WTs in a wind farm (WF), FR provision could be shared by WTs to provide an aggregate response that would facilitate the task of the transmission system operator (TSO). Several strategies that implement this idea for coordinating WTs' operation can be found in literature. In these strategies, the power output of the WTs are coordinated to fulfil the total demanded power determined by the WF operator or the TSO, while achieving different goals.
Some publications have addressed the specific problems of reducing the power output of WFs (deloading). The authors in ref. [11] proposed an optimal power dispatch to distribute a command of reducing the total power while minimising the active power losses in a WF. In refs. [12][13][14] authors proposed methods to equalise the utilisation level of the WTs in applications where they took part of the regulation of the frequency. In ref. [15], the authors showed that the drop in energy yield during power sharing could be reduced by adjusting the utilisation level according to each WT's rotor speed. In ref. [16], the authors further developed an optimal power sharing scheme to maximise the total rotational kinetic energy stored in WTs and fulfil the demand order. In ref. [17], a game theory-based distributed algorithm is devised to allocate the deloading demand IET Renew. Power Gener. 2021;15:1005-1018.
wileyonlinelibrary.com/iet-rpg Regarding the case where WFs are expected to generate more power (overload), the authors in ref. [18] proposed an optimal power dispatch to minimise the loss of wind energy caused by FR provision while satisfying the pre-calculated overload capacity constraint.
Some work has considered quantifying the adequacy of WTs to overload or deload depending on their relative position in the WF and the direction of the wind to account for the wake effect [19]. For deloading, the authors argued that the back row of WTs would be given priority whereas the front row would be priorised for overloading.
The aforementioned works [11][12][13][14][15][16][17][18][19] considered steady-state (i.e. long-term constant speed) models. This model has been challenged in some recent publications. The authors of ref. [20] developed an optimal power sharing method for overloading to maximise the rotational kinetic energy of WTs considering the dynamic evolution of the WT during the process. Similarly, the authors in refs. [21][22][23] considered the WTs' transient mechanical loads, and devised an optimal dispatch to minimise them during deloading events.
The purpose of occasional FR is to contribute to the regulation of the frequency under rare large-scale events. In an underfrequency event, typically caused by a large loss of generation, the WF is expected to deliver additional power (i.e. overload) while in an overfrequency event (e.g. large loss of load) the WF is required to reduce the power output (i.e. deload). On the other hand, the purpose of continuous FR is to maintain system frequency closer to 50Hz (or 60Hz) under normal operation and requires the WF to modulate its active power continuously.
In this paper, the power sharing problem is formulated as an optimisation problem that takes into account the dynamics of WTs in order to ensure their proper operation at all times while delivering the required FR. Similar to refs. [16,17,20], the objective function is selected as maximising the total kinetic energy of the WF. This choice enables WTs maintaining a larger reserve of kinetic energy, which could be later released to the grid when the WTs are asked to recover to their original optimal operating point, thus reducing the loss of revenue caused by the continuous FR provision.
During operation of the WTs, three constraints are considered: (1) pitch angle action should be avoided (soft constraint), (2) the generator torque must be kept below maximum rated value and (3) the rotor speed must be kept above the minimum rotor speed at all times. The aim of (1) is to reduce the long term loss of generated energy caused by the provision of FR, and reduce the wear and tear of pitch systems [25,26], whereas (2) and (3) ensure safe operation of WTs and reduce wear and tear, extending the lifespan of WTs [27].
The optimisation problem is solved by a novel power sharing strategy that uses periodic communication with each WT. In the proposed strategy, the dynamics of the rotor speed of each WT are approximated linearly around their current operating point. Then, the three additional constraints are evaluated and the solution of the optimisation problem is obtained by the proposed centralised and/or distributed algorithm.
The main advantages of the proposed strategy are: 1. It is able to handle both the deloading and overloading of WF, which makes it suitable for FR provision in both underfrequency and overfrequency events, and continuous FR under normal operation. 2. It reduces the likelihood of pitch angle action to reduce the loss of yield caused by FR and it avoids exceeding a given WT torque limit. 3. It can be implemented in a distributed way through the proposed distributed algorithm, wherein only one-way communication is required. While for most existing methods, they require two-way communication between the wind farm controller (WFC) and the local WTs, or among WTs.
The paper is structured as follows: Section 2 illustrates the dynamics of WT. Section 3 formulates the optimisation problem and derives the three additional constraints. Section 4 develops the proposed strategy including the centralised and distributed algorithm of solving the optimisation problem. Section 5 verifies the proposed strategy.

WIND TURBINE DYNAMICS
In this paper, we consider two most-widely-used variable speed WTs, namely, doubly-fed induction generator (DFIG) WTs and permanent magnet synchronous generator (PMSG) WTs. We use the power characteristic of a WT given by [28] Here, P t is the wind turbine power; is the air density; R is the turbine radius; v w is the wind speed; w t is the rotor speed; C p FIGURE 2 Pitch angle control of WT is the turbine efficiency, which is related to the tip-speed ratio and pitch angle . The approximation of turbine efficiency C p ( , ) is [28] C p ( , ) = 0.22 According to refs. [29,30], the WT one-mass model is adequate for investigating the WT dynamics during FR provision. Thus, the dynamics of rotor speed w t are modelled aṡ where J t is the total mechanical inertia of the rotating shafts; T t is the turbine torque; T g is the generator torque; P t is the turbine power in Equation (1); P g is the output power generated by the generator. Through proper control of the converters of DFIG and PMSG WT [31], P g can be regulated to track its power reference with adequately fast response. Thus, P g can be regarded the same as its power reference. The speed control of the turbine is achieved through its electrical power set-point as well as through its blade pitch angle. In this paper, we consider two modes of operation: maximum power production (MPP) and FR.

MPP mode
In MPP mode, the objective is to maximise the power production. To achieve this, the output power is controlled to optimise the power extraction of WT by regulating w t , i.e. maximise C p . The control law is [32] where C opt p is the optimal power coefficient at = 0; opt is the optimal tip speed ratio which corresponds to the optimal power coefficient C opt p . The pitch angle of WT in MPP mode is controlled to prevent exceeding the maximum rated speed, i.e. to make w t ≤ w t,rated at all times as in ref. [13]. Under normal operation, the pitch angle is kept at zero in order to maximise the yield and it is only increased at high wind speeds in order to keep the rotor speed from exceeding its rated value. The control law is as shown in Figure 2.

FR mode
In FR mode, WT deviates from its MPP operating point and the objective becomes to provide FR service. P g is set as Here, P g,ref is the output power reference, which is determined by the power sharing algorithm; P min and P max are the minimum and the maximum output power respectively, which are determined by the physical limitations of the WT. The pitch angle of WT in FR mode is controlled in the same manner as in MPP mode, shown in Figure 2.

The Power Sharing Problem
In this paper, we consider a WF consisting of N WTs. The scenario of the considered problem is expressed as follows: 1. Initially, the WF is operating in MPP mode, i.e. each WT in the WF is running in MPP mode. The total output power of WF is P 0 . 2. Then, the WF is asked to provide FR with total demanded power P total = P 0 + ΔP for a period of time T .
Here, the value of P total or ΔP is given by the WF operator or the TSO, and in this paper we focusing on the coordination of WTs in the WF to fulfil this power order (i.e. how to optimally distribute this P total or ΔP among the WTs in the WF). Through this setting, the WF provides an aggregate FR to the grid, which will facilitate the task of the TSO. We consider the WF providing two FR services: occasional and continuous FR shown in Figure 1. Thus, the setting of ΔP will be different under different FR cases.
Generally, for occasional FR, ΔP is either positive (corresponds to overload in an underfrequency event) or negative (corresponds to deload in an overfrequency event), and the duration T is within seconds. While for continuous FR, ΔP is time-varying and its modulus and sign are determined by the system frequency. The duration T could be minutes or longer. The value of ΔP could be set as , or as a constant value during underfrequency or overfrequency events as discussed in refs. [5,6].

Effects of FR provision
In order to find the criterion for distributing P total , we first look at the effect of overloading/deloading on a single WT in terms of its operating revenue. The operating revenue is related to the power yield of the WT and it is the key performance indicator for the owners of the WT. A series of tests is performed on a WT transient simulation model described in ref. [33]. The rated wind speed is v w,rated = 10 m/s , all relevant characteristics of the WT are summarised in Appendix. The simulation starts from an MPP operating point at different given wind speeds. The WT is firstly asked to deviate from MPP by ΔP for a period of time T , and then it is asked to recover to original MPP point. We consider the variation of total delivered power from the WT over the whole period (overload or deload period T + recovery period T r ), ΔP * , as We note that ΔP * represents the effect of overloading or deloading on the long term energy delivered to the grid. In Equation (6), ∫ T +T r 0 P g (t )dt is the energy delivered to the grid during the whole period, and it is equal to the captured wind energy ∫ T +T r 0 P t (t )dt because the WT is finally recovered to the original MPP point. Thus, ΔP * can also be computed as

Effects of FR provision: Deload
The WT dynamics during deloading operation are illustrated in Figure 3(a) and Figure 3(b). For low wind speed (v w = 9 m/s < v w,rated ), as shown in Figure 3(a), the rotor speed increases and P t decreases as the WT is shifted away from the MPP point. After a certain amount of time, when the rotor speed reaches the maximum rated speed w t,rated (marked by the red dotted arrow), the pitch angle control kicks in to avoid over-speeding. This leads to the great decline of P t (indicated by the blue dashed arrow), which will waste a great amount of energy according to Equation (7). Thus, ΔP * should be less than zero when deloading WTs with v w < v w,rated . For high wind speed (v w = 11 m/s > v w,rated ), the WT dynamics are shown in Figure 3(b). As the rotor speed is initially at the maximum rated speed w t,rated , the pitch angle increases as soon as the deloading begins. As can be seen, the increase of pitch angle causes the sharp reduction of P t , leading to a great loss of energy. Thus, ΔP * < 0 when deloading WTs with v w ≥ v w,rated .
The quantitative relationship between ΔP * and ΔP under different wind speeds is shown in Figure 3(c). As expected, the ΔP * becomes more negative (meaning waste more energy) the more negative ΔP is. Further, we observe two different behaviours (i.e. slopes) between ΔP * and ΔP. The reason behind this is the pitch angle action. When operating at v w ≥ v w,rated (the grey area), the WT runs at its rated speed; therefore, pitch angle action will happen as soon as deloading starts in order to prevent over-speeding. This leads to great loss of energy (as shown in Figure 3(b)), and increases the wear and tear of pitch systems [25,26]. In contrast, for v w < v w,rated (the red area), the rotor speeds are initially below the rated speed. Thus, WTs can have some headroom to store the extracted wind Relationship between ΔP * and ΔP under different wind speeds energy they can't deliver during deloading, which is later released to the grid in the recovery period, leading to the less loss of energy. Therefore, ΔP * only drops significantly if the deloading event is severe enough for them to reach their maximum rotor speed.
When comparing the WTs with v w < v w,rated under very small ΔP (marked by the red circle in Figure 3(c)), we observe that the loss of energy is lesser than those WTs at higher wind speeds (see Figure 3(c)).
The above observations indicate that (1) deloading WTs will always cause loss of energy. (2) in order to reduce the loss of energy and wear and tear of the pitch systems, it is best to deload WTs operating at higher wind speeds while keeping their rotor speed low enough to avoid their pitch angle action.

Effects of FR provision: overload
Regarding the overloading operation of WT, the WT dynamics are shown in Figure 4(a, b). For low wind speed (v w = 9 m/s < v w,rated ), as shown in Figure 4(a), the rotor speed decreases due to the overloading operation. As a result, P t decreases and T g increases. After some time, T g reaches the rated value T g,rated (marked by the red dotted arrow), which will increase the wear and tear of WT. During the whole process (both FR and recovery period), P t is always less than P mpp as the WT is not running at the MPP point, meaning that the long-term power yield of WT is reduced compare with operating the WT at MPP. Thus, according to (7), ΔP * < 0 when overloading WTs with v w ≤ v w,rated .
For high wind speed (v w = 11 m/s > v w,rated ), the WT dynamics are shown in Figure 4(b). T g exceeds the rated value immediately when the overloading starts. Meanwhile, the pitch angle is also decreased to avoid the decline of rotor speed. This leads to the increase of P t , causing an additional production of energy. Thus, ΔP * > 0 when overloading WTs with v w ≥ v w,rated .
The quantitative relationship between ΔP * and ΔP under different wind speeds is shown in Figure 4(c). The results show that there can be a positive correlation between ΔP * and ΔP when v w > v w,rated (the green area). The reason for this is that WTs operating at high wind speeds are likely to have a blade pitch angle greater than zero and will be able to increase the power they extract from the wind during an overloading event to prevent reducing their rotor speed. This causes the increment of P t , leading to additional production of energy as explained in Figure 4(b). However, for these high wind speed WTs, their torque are already at the rated value, overloading them will make the torque exceeding the rated value T g,rated , causing potential damage and increasing wear and tear which will shorten the lifetime of WT [27]. Thus, these WTs should not be overloaded unless the turbine is oversized by design.
On the contrary, for v w ≤ v w,rated (the red area in Figure 4(c), WTs can have some torque headroom for overloading. They slow down and move away from their MPP which drops their efficiency at extracting energy from the wind. This causes the reduction of P t , leading to the waste of energy. These WTs also need to stop overloading when reaching their minimum speed w t,min or their maximum torque T g,rated . We also notice that in the red area, the loss of energy of WTs with higher wind speed is lesser.
The above observations show that (1) overloading WTs should pay great attention to avoid exceeding the WT ratings to ensure safe operation of WTs; (2) higher WTs are more suited for overloading.  In summary, in the distribution of total demanded power, higher WTs should overload/deload more while satisfying the three additional constraints in Table 1. We note that previous works, such as refs. [11][12][13][14][15][16][17][18][19][20][21][22][23][24], did not consider the Constrains 1) and 2) in Table 1 during the calculation of the power sharing of the WTs.

Formulation of the optimisation problem
In this paper, the power sharing problem is formulated as an optimisation problem with its objective function as maximising the total kinetic energy E k,total after the FR is being provided. The purpose of this choice is twofold [16,17,20]. On the WF side, this choice aims to save more kinetic energy after FR provision given the same total demanded power. The additional reserved kinetic energy could be later released when asked to recover to the MPP, ultimately reducing the loss of revenue caused by FR provision, which is of great interest for continuous FR provision. On the TSO side, when WF is asked to overload in an underfrequency event, this choice takes into account the efficient consumption of WT kinetic energy, which makes the WF could sustain for a longer period of time given the same demanded power. This allows more time for the secondary response devices to kick in to regulate the system frequency. Consequently, this enables the TSO to use the cheaper frequency response devices with lower power ramping rates and achieve economic benefits.
The formulated optimisation problem together with the three additional constraints in Table 1 is defined as: where P ti is determined by Equations (1) and (2) and the pitch angle is controlled as shown in Figure 2.

OPTIMAL POWER SHARING STRATEGY
In this section, a novel strategy is presented to calculate the power curves of each WT, i.e. P gi,ref (t ), for the optimisation problem (8).
The proposed strategy is employed in a system which is being updated periodically every T c seconds. The illustration of the proposed strategy is shown in are approximated by their linearised model; then the three additional constraints on the output power of each WT are calculated; after that the optimal solution when ignoring the output power constraints is computed; finally, the solution considering the output power constraints is obtained by the proposed centralised and/or distributed algorithm.

Rotor Speed Dynamics Approximation
The rotor speed dynamics of the ith WT is determined by Equations (1), (2) and (3) and the pitch control in Figure 2. During t ∈ [kT c , (k + 1)T c ], it can be approximated using its linearised model, shown below. Here, ΔP ti,ref is defined as the output power deviation with respect to turbine power P ti (kT c ). To evaluate the accuracy of the approximation of the rotor speed, we simulated different step changes of ΔP starting from the MPP operating point with different T c on a WT with a two-mass drivetrain [33] with parameters shown in Appendix and initial wind speed v w = 7 m/s (see Figure 6). The results show that the approximated rotor speed w t,approx given by Equation (9) is very close (about 99%) to the real value w t,real , even considering the effect of the transient speed oscillation caused by the sudden torque change, which is visible over the first few T c intervals.

Constraints Derivation
The three constraints in Table 1 are considered in the optimisation problem (8). Accordingly, at every t = kT c , k = 0, 1, … , P lmt gi,ref (kT c ) is the output power which makes w ti ((k + 1)T c ) = w t,min ; P min gi,ref (kT c ) is the output power which makes w ti ((k + 1)T c ) = w t,rated . In order to meet the three constraints in Table 1, The three output power constraints are calculated by the linearised rotor speed model, show below. Using Equation (9), the generator torque T gi (t ) during t ∈ [kT c , (k + 1)T c ] is approximated as , t ∈ [kT c , (k + 1)T c ] (11) Then, w ti and T gi at t = (k + 1)T c are predicted as Then, by setting T gi ((k + 1)T c ) = T g,rated , w ti ((k + 1)T c ) = w t,min and w ti ((k + 1)T c ) = w t,rated , the three output power constraints at t = kT c are calculated as: T g,rated J ti w 2 ti (kT c ) − P ti (kT c )J ti w ti (kT c ) T g,rated T c + J ti w ti (kT c ) .
Then, together with the physical limitations of the output power, the output power constraints at t = kT c are

Reduced Problem
We note that P total and the wind speed v w may vary over time. Thus, it is not practical to directly obtain the solution for Equation (8) considering these uncertainties. Here, to handle this, the optimisation goal at every t = kT c is reduced to maximise the total kinetic energy at t = (k + 1)T c rather than at t = T . After this treatment, it is now reasonable to assume that the wind speed and the total demanded power are constant during t ∈ [kT c , (k + 1)T c ]. The reduced optimisation problem is Ignoring the inequality constraints on ΔP ti,ref , the solution of Equation (16) can be obtained using the Lagrange multipliers [34]. Its solution is given by: where ΔP t,total (kT c ) = P total (kT c ) − Here, ΔP t,total (kT c ) is the total demanded power deviation with respect to the current total turbine power, and it is varying during all period of FR provision because w t and C p are changing. E ki is the kinetic energy of the ith WT. The non-constraints solution (18) suggests that, to maximise the total kinetic energy, the total demanded power should be weighted-shared based on each WT's kinetic energy. This means that the WTs with higher rotor speed should have more allocated ΔP ti,ref . 1     The final solution of Equation (16) should be determined by both the non-constraints solution (18) and the inequality constraints on ΔP ti,ref .

Centralised Optimal Power Sharing Algorithm
The optimisation problem defined in Equation (16) is solved using the optimal sharing algorithm proposed in [35] modified to consider both upper and lower bounds on ΔP ti,ref . This is deployed in the WFC which communicates with local WTs periodically every T c seconds. As shown in Figure 7(a), each WT sends P ti and E ki to WFC. The WFC then calculates the three output power constraints according to Equation (13). And then it computes the maximum and minimum power constraint ΔP ti,ref and ΔP ti,ref according to Equation (15). After that, the WFC optimally allocates P total under these constraints using Algorithm 1.

Distributed Optimal Power Sharing Algorithm
The centralised method in Section 4.4 would need to collect the operational data from all WTs, which will increase the communication overheads. A distributed approach is hence here proposed.
During t ∈ [kT c , (k + 1)T c ], the power of each WT is ) .

(21)
Here, P gi,ref (kT c ) and P gi,ref (kT c ) are the power constraints during t ∈ [kT c , (k + 1)T c ], which are computed by WTs through Equations (13)- (15). In Equation (21), (t ) is computed in the WFC, which is proportional to the error between the demand and actual delivered power from WF which is measured at the collector bus of WF. The magnitude of should be set large to achieve fast convergence and the maximum value of should be bounded by the maximum power ramp rate of WT. Equation (20b) originates from the non-constraints solution (18) and (t ) is the same for all WTs, which can be interpreted as the power deviation factor from the current P ti (kT c ). In Equation (20c), (t ) is set as the integral of (t ), making that P gi,ref (t ) is gradually moving to the solution of Equation (16).
In the implementation of this distributed method, the WFC measures the power at the collector bus of WF, i.e. ∑ N i=1 P gi (t ), and computes (t ) in Equation (21). Then, it broadcasts (t ) to all WTs. For the WTs in the WF, they compute the power constraints P gi,ref (kT c ) and P gi,ref (kT c ) at every t = kT c , and control P gi,ref (t ) according to Equation (20a-c). The exchange of variables is shown in Figure 7b. Comparing to the existing methods [11][12][13][14][15][16][17][18][19][20][21][22]24] and the proposed centralised algorithm, this distributed approach removes the data collection from all WTs, and it only requires one-way communication between WFC and local WTs, thus reducing the communication overheads.

CASE STUDY
In this section, the proposed strategy is tested on a WF with 10 WTs with parameters shown in Appendix in providing both occasional and continuous FR.
In the simulation, the proposed strategy is employed with T c = 0.5s, which is a realistic value for the application in real WFs. Each WT is running in MPP mode firstly with maximum output power . At s, the WF is asked to deliver total power for seconds.
Two existing methods, the variable utilisation level (VUL) [15] and the optimal overload a (OOS) [20], are tested along with the proposed method for comparison. The VUL only considers the deloading region and adjusts the utilisation level according to WT rotor speed to ensure less reduction of energy production. The distribution law of VUL is }P mpp,i The solution of OOS is to control the P gi,ref to maintain the partial derivative of P ti with respect to E ki at the same value.

Occasional Frequency Response: Overfrequency Event
In this section, the WF is asked to provide occasional FR in an overfrequency event. In this test, the wind speed of each WT are v w10 = 11.5 m/s and v wi = v w10 × (0.76 To mitigate the overfrequency event, WF needs to decrease its total power [36,37]. Similar to ref. [15], the total demanded power of WF is set as constant during FR provision. In the test, we set ΔP = −1.0MW, and the duration time is set as T = 10 s. As all WTs are demanded to deload, the rotor speed will increase, thus the generator torque will decrease. Hence, no WT will violate Constraint (2) and (3) in Table 1.
The results of VUL are shown in Figure 8(a). The VUL does not consider the Constraint 1) in Table 1. Instead, ΔP is always weighted-distributed based on P mpp,i (w t,max − w ti ). As a result, the pitch angle of WT 10 are increased immediately when the FR provision starts. As time passes, the pitch angles of WT 9 , WT 8 , WT 7 are also increased. The increasing of pitch angles will lead to additional loss of yield and wear and tear of pitch systems, which would be a concern if the WF deloading happens regularly. The proposed strategy can avoid this pitch angle action and still deliver the same FR service as shown next.
The results of the proposed strategy are shown in Figure 8(b). As the total demanded power is always allocated under Constraint (1), none of the pitch angles are changed during the FR provision. This will increase the stored kinetic energy of WF during the FR provision. According to Equation (18), the total demanded power is weighted-shared based on each WT's kinetic energy under Constraint (1). Consequently, as shown in Figure 8(b), WT 10 does not participate in deloading as its rotor speed is already at the rated value; for the remaining WTs, as w t 9 is greater than the other WTs, WT 9 has more shared ΔP (i.e. deload more) when the FR provision starts. As time passes, w t 9 and P g9,ref increase as well according to Equations (13)- (15).
At t = 201.5s, P g9,ref becomes greater than WT 9 's nonconstraints solution according to (18). Thus, P g9 is allocated to  Figure 8(c) shows the rotor speed variation of the proposed strategy. Δw t 10 is always zero as it does not overload under the proposed strategy. When the FR provision starts, Δw t 9 is the largest as it has more shared ΔP. With time passing by, P g9 is allocated to be P g9,ref at t = 201.5s; therefore,  Figure 9(a) shows the discrepancy between the solution of proposed centralised algorithm P c gi and the solution of proposed distributed algorithm P d gi under different . We observe that when is increased from 5.0 × 10 −7 to 5.0 × 10 −6 , the speed of convergence is increased and the discrepancy is reduced as well. The initial discrepancy at t = 200 s is due to the initial deload setting. The discrepancy is virtually eliminated 0.2 s for = 5.0 × 10 −7 (0.01 s for = 5.0 × 10 −6 ) after the event is detected and remains within 10% of the centralised solution. As shown in Figure 8(b), at t = 201. 5, 203, 205, 206.5, 208, 209 s, the non-constraint solution of WT i reaches P gi,ref , and thus P gi is increased to P gi,ref . This causes the discrepancy at these times. Furthermore, the initial discrepancy can be eliminated if, at the onset of the FR provision, the initial solution is calculated at the WFC and broadcasted to all participating WTs. Figure 9(b) compares the total kinetic energy of VUL and the proposed strategy. As can be seen, due to the consideration of Constraint (1), the total kinetic energy of the proposed strategy is greater than VUL during the whole FR provision period. At the end of FR provision, the stored total kinetic energy of the proposed strategy is 9.73 MJ, which is 18.66% larger than VUL (8.20 MJ). Here, we note that, the advantages of the proposed strategy also include: (1) the VUL is a centralised algorithm, wherein two-way communication between WFC and local WTs is required. While for the proposed strategy, it can be implemented in a distributed way through the proposed distributed algorithm (see Figure 9(a)), wherein only one-way communication is required, thus reducing the communication overheads. (2) the proposed strategy is able to handle both the overloading and deloading of WTs. While the VUL only considers the deloading of WTs.

Occasional Frequency Response: Underfrequency Event
In this section, the WF is demanded to provide occasional FR in an underfrequency event under the same wind conditions as before.
In this case, the WF needs to increase its total power to mitigate the underfrequency event [9,10,38]. Similar to [20], the total demanded power of WF is set as a constant value during FR provision. In the test, we set ΔP = +0.7 MW and T = 10 s.
As shown in Figure 10(a), the OOS [20] distributes the total power among low wind speed WTs (i.e. m/s) without considering the torque constraints. As a result, T g9 and T g8 exceed the rated value T g,rated when the FR provision starts. Moreover, with time passing by, the generator torque of all WTs increase and T gi , i = 9, 8, 7, 6 exceed T g,rated by 10.37%, 7.46%, 4.58% and 1.74% respectively after FR provision. This will require sizing the turbine for higher torque and may increase wear and tear over time [27]. The proposed method can prevent this and still provide the same service as shown next.
The result of proposed strategy (centralised version) is shown in Figure 10(b). As can be seen, the generator torque of all WTs are always below T g,rated because of the consideration of the Constraint (2) and (3) in Table 1. As T g10 is already at the rated torque, WT 10 does not overload during the whole FR provision. When the FR provision begins, P gi,ref , i = 9, 8, 7, which are determined by the current torque according to Equation  Figure 10(c) shows the rotor speed variation of the proposed strategy. Δw t 10 is always zer o as it does not overload under proposed strategy. Δw t 9 , Δw t 8 and Δw t 7 are larger than Δw ti , i = 1, 2, … , 6 because P g9 , P g8 and P g7 are saturated by P g9,ref , P g8,ref , and P g7,ref respectively. With time passing by, P g6,ref is allocated FIGURE 10 Results of WF in providing occasional FR in an underfrequency event, (a) OOS, (b) proposed strategy-centralised optimal power sharing algorithm, (c) rotor speed variation of proposed strategy to be P g6,ref after t = 202.5 s. This slows down the decrease rate of Δw t 6 , and thus Δw t 6 becomes larger than Δw ti , i = 1, 2, … , 5 after t = 205.5 s (indicated by the red arrow). A similar pattern can be observed for Δw t 5 . When P g5,ref is saturated to be P g5,ref after t = 205 s, Δw t 5 becomes larger than Δw ti , i = 1, 2, … , 4 after around t = 208.5 s (indicated by the black arrow). Figure 11(a) shows the discrepancy between the distributed solution P d gi and the centralised solution P c gi . As the WF is asked to overload, the initial discrepancy at t = 200 s is expected. The discrepancy is virtually eliminated 0.2 s after the event is detected and remains within 2% of the centralised solution. Meanwhile, we also observe that when is increased from 5.0 × 10 −7 to 5.0 × 10 −6 , the speed of convergence is increased and the discrepancy is reduced as well. Furthermore, the initial FIGURE 11 WF providing occasional FR in an underfrequency event, (a) discrepancy between the centralised solution P c gi and the distributed solution P d gi : |P c gi −P d gi | P c gi × 100%, (b) Comparison of the total kinetic energy of WF discrepancy can be eliminated if, at the onset of the FR provision the WFC calculates and broadcasts the solution to all participating WTs. Figure 11(b) presents the curves of the total kinetic energy of the WF. We notice that the consumed kinetic energy of the proposed strategy is 7.19 MJ and that of the OOS is 7.16 MJ, which is only slightly smaller than proposed strategy. However, in the OOS the torque constraints of WT are violated, which will increase the wear and tear of WT. Here, we note that, the benefits of the proposed strategy also include: (1) it can be implemented in a distributed way through the proposed distributed algorithm (see Figure 11(a)), wherein only one-way communication is needed. While for the OOS, two-way communication between WFC and local WTs is required. (2) the proposed strategy considers both the overloading and deloading of WTs. While the OOS only considers the overloading region.

Continuous Frequency Response
In this section, the WF is required to provide continuous FR to regulate the system frequency in normal operation. In this test, the average wind speed of each WT is set the same as that in

FIGURE 12
Test settings in Continuous FR, (a) wind speed of each WT, (b) total demanded power occasional FR; apart from that, an additional time-varying signal is added, as shown in Figure 12(a). The WF is asked to provide continuous FR for 120 s. In this case, P total is adjusted up and down to regulate the system frequency, and its value is shown Figure 12(b). For the comparison, VUL is used in the deload region, and OOS is used in the overload region. In the following, it is named as VUL+OOS.
As shown in Figure 13(a), in VUL+OOS, WT 9 is running in over torque for 109s (91% of the FR period) in total; while in the proposed strategy, there is almost no over torque operation during the whole FR provision period (shown in Figure 13(b)). As shown in Figure 13(c), the kinetic energy of the proposed strategy at the end of the simulated period is 6.20 MJ greater than OOS+VUL, which is equivalent to 51.67 kW average (0.52% of the set-point of the WF at the beginning of the period).

FINAL REMARKS
WF is expected to operate as a single aggregate entity in providing FR service to facilitate the planning of TSO. In this paper, a novel strategy is presented to allocate the total demanded power such that the total kinetic energy is maximised after providing FR taking into account the WT dynamics. The proposed strategy includes three additional constraints to ensure the safe operation of WT, and reduce wear and tear and loss of energy. Using a dedicated scheme, which updates in a periodic manner, the three constraints are computed by approximating the rotor speed dynamics and the optimal solution is then obtained using proposed centralised and/or distributed algorithm. Case studies are performed on a WF providing both occasional and continuous FR, and have shown the efficacy of proposed strategy and the gained enhancement over existing methods.

APPENDIX A
The parameters of studied WT are shown in Table A.1.