Active Power Dispatch for Supporting Grid Frequency Regulation in Wind Farms Considering Fatigue Load

: This paper proposes an active power control method for supporting grid frequency regulation in wind farms (WF) considering improved fatigue load sensitivity of wind turbines (WT). The control method is concluded into two parts: frequency adjustment control (FAC) and power reference dispatch (PRD). On one hand, the proposed Fuzzy-PID control method can actively maintain the balance between power generation and grid load, by which the grid frequency is regulated when plenty of winds are available. The fast power response can be provided and frequency error can be reduced by the proposed method. On the other hand, the sensitivity of the WT fatigue loads to the power references is improved. The explicit analytical equations of the fatigue load sensitivity are re-derived to improve calculation accuracy. In the process of the optimization dispatch, the re-deﬁned fatigue load sensitivity will be used to minimize fatigue load. Case studies were conducted with a WF under di ﬀ erent grid loads and turbulent wind with di ﬀ erent intensities. By comparing the frequency response of the WF, rainﬂow cycle, and Damage Equivalent Load (DEL) of the WT, the e ﬃ cacy of the proposed method is veriﬁed.


Introduction
Wind energy is one of the rapidly growing renewable energy sources [1][2][3][4]. By far largest wind power market, China installed an additional capacity of 19 Gigawatts, and continues its undisputed position as the world's wind power leader [5][6][7][8]. With the rapid growth of installed wind power capacity, its proportion in the power grid continues to increase. Wind power generally does not participate in frequency regulation due to the decoupling of the WT rotor and the grid frequency, so large-scale wind power access to the grid will significantly weaken the frequency regulation capability of the grid [9,10]. With the continuous increase of the wind power penetration rate, the influence of wind power volatility and uncertainty on the frequency regulation of power systems is also increasing [11][12][13][14][15]. In order to meet grid frequency requirements, active power control (APC) method for wind power is critical to actively maintaining a balance between power generation. However, the power reference is frequently changed by the FAC to adapt to changes in the grid frequency when the WF participates in frequency regulation. The frequent action of the pitch angle of the WT and the torque of the generator is caused by this change, resulting in an increase in the fatigue loads of WTs [16][17][18][19]. Therefore, it is necessary to study which control methods support grid frequency regulation without affecting or even reducing the fatigue loads of WTs.
In recent research, some methods for supporting grid frequency regulation have been proposed. An approach for participation of doubly fed induction generator (DFIG) based wind farms in power system short-term frequency regulation to simplify the high-order average frequency models of bulk power systems was proposed in [20,21]. In [22,23], the inertia constant and primary power reserve for a variable speed wind turbine that operates at derated conditions were formulated in a wind farm to support short-term frequency control in power systems. Reference [24] presents a model-based control method based on Model Predictive Control (MPC) method and on a Kalman-like estimation algorithm to improve the contribution of wind power generators to short-time primary frequency regulation in electric power systems. A frequency control support function responding proportionally to frequency deviation is proposed to take out the kinetic energy of wind turbine for improving the frequency response of the system in reference [25]. However, the grid frequency is only stabilized in a short period of time using these methods. Reference [26] presented two APC methods that are developed based on adaptive pole placement control (APPC) and fuzzy proportional-integral (PI) control approaches to provide rapid power response. The grid frequency can be stabilized for a long time by WF power reduction. However, it can be seen from the DELs results that the fatigue load of WT has been increased using this method. For wind farm power distribution, as long as the power scheduling requirements are met, the fatigue load can be minimized by coordination between the WTs [27]. In [28,29], a distributed MPC based APC method of WFs was presented to reduce the fatigue loads of WTs. In [30], a load sensitivity based optimal active power dispatch algorithm is proposed for wind farms. The explicit analytical equations of the load sensitivity are derived to improve the computational efficiency of controller. Two fatigue loads experienced by shaft torque and tower bending moment are considered.
An APC method is proposed in this paper. This method is used to provide fast power response while minimizing fatigue loading on the WT, which is used to overcome the above problems. The frequency adjustment control of WF was developed based on the Fuzzy-PID control method. The method is designed to track various forms of load while maintaining grid frequency stability when plenty of wind is available. The total active power obtained by the FAC is proportionally distributed to the wind farms. The fast power response can be provided and frequency error can be reduced by the proposed method compared with traditional methods. Based on the resulting power reference value, the PRD adjusts the power reference of WT within the WF to track power. The model of the WT is improved and the analytical equation for fatigue load sensitivity is re-derived. Compared with the original method, the calculation accuracy is improved by the improved method. Meanwhile, PRD minimizes fatigue loads of tower thrust and shaft torque variations based on the improved fatigue load sensitivity model. Following this, a complete system simulation model is built in MATLAB to verify the effectiveness of the proposed method.
The main contributions of this work are described as follows: A control structure support frequency adjustment for large-scale WF is proposed in Section 2. Section 3 designs a Fuzzy-PID controller to respond and recover grid frequencies more quickly. The fatigue load sensitivity model of the wind turbine is improved and the explicit analytical equations of the fatigue load sensitivity are re-derived in Section 4. The improved model can adapt to wind conditions with different turbulence intensity to minimize the shaft torque and tower bending moment fatigue loads. The simulation results demonstrate the effectiveness of the method in Section 5. Section 6 concludes the paper.

Control Structure of WF Participates in Frequency Regulation
The proposed WF control architecture proposed in the paper is shown in Figure 1. The control architecture is mainly divided into two parts.
The first part is FAC of the WF. The measured grid frequency f meas is used as a feedback signal to set up active power control in real-time and maintain the balance between power generation and grid loads. The demanded power of WF P WF demand is calculated by FAC and delivered to the PRD. The grid frequency is regulated by FAC to its rated value despite a changing grid load. The second part is PRD of the WF. This part takes all the WTs as the unit and performs power tracking on the power value assigned by the superior to realize the active power adjustment of the WT in the area under his

Fuzzy-PID Control Method of Supporting Grid Frequency Regulation for WF
The frequency fluctuation of the power system is caused by the imbalance of the generation and consumption of active power. In order to maintain the frequency stability of the power system, the system frequency must be maintained by active power control. The greater the balance between power generation and consumption, the smaller the frequency fluctuations, and the higher the electricity quality. This paper considers active power control at an entire WF level within the general structure shown in Figure 1. A typical large-scale WF including N wind turbines is included.
When the penetration of wind power in the power grid is relatively low, the impact of wind power participation in frequency regulation is minimal, so the wind farm frequency regulation method uses open-loop control. However, with the increase of wind power penetration rate, this open-loop control method cannot adapt to this situation and will cause the frequency to produce steady-state error. Traditional proportion-integral-derivative (PID) controllers are one of the most widely used controllers in industrial applications, and they can eliminate the steady-state error of the frequency [31][32][33]. However, it is difficult to adapt to a wind power conversion system with strong nonlinear dynamic characteristics, which needs to be improved. The fuzzy control system establishes a fuzzy rule table suitable for the actual production process by summarizing expert knowledge and operational experience. It not only reduces the rise time and overshoot of the output response, but also reduces the sensitivity of the system to interference and increases the stability of the system [34]. Fuzzy-PID can adjust the size of the parameters for P, I, and D online, and it can adapt well to dynamic systems [35,36]. In this paper, the Fuzzy-PID controller is used to regulate the grid frequency. In the APC architecture proposed in this paper, the input of the fuzzy controller is the frequency deviation ferror of the grid, and the output is the demanded power of WF P demand WF . The demanded power is then sent to PRD. As it is shown in Figure 2, grid frequency error ferror is calculated by

Fuzzy-PID Control Method of Supporting Grid Frequency Regulation for WF
The frequency fluctuation of the power system is caused by the imbalance of the generation and consumption of active power. In order to maintain the frequency stability of the power system, the system frequency must be maintained by active power control. The greater the balance between power generation and consumption, the smaller the frequency fluctuations, and the higher the electricity quality. This paper considers active power control at an entire WF level within the general structure shown in Figure 1. A typical large-scale WF including N wind turbines is included.
When the penetration of wind power in the power grid is relatively low, the impact of wind power participation in frequency regulation is minimal, so the wind farm frequency regulation method uses open-loop control. However, with the increase of wind power penetration rate, this open-loop control method cannot adapt to this situation and will cause the frequency to produce steady-state error. Traditional proportion-integral-derivative (PID) controllers are one of the most widely used controllers in industrial applications, and they can eliminate the steady-state error of the frequency [31][32][33]. However, it is difficult to adapt to a wind power conversion system with strong nonlinear dynamic characteristics, which needs to be improved. The fuzzy control system establishes a fuzzy rule table suitable for the actual production process by summarizing expert knowledge and operational experience. It not only reduces the rise time and overshoot of the output response, but also reduces the sensitivity of the system to interference and increases the stability of the system [34]. Fuzzy-PID can adjust the size of the parameters for P, I, and D online, and it can adapt well to dynamic systems [35,36]. In this paper, the Fuzzy-PID controller is used to regulate the grid frequency. In the APC architecture proposed in this paper, the input of the fuzzy controller is the frequency deviation f error of the grid, and the output is the demanded power of WF P WF demand . The demanded power is then sent to PRD. As it is shown in Figure 2, grid frequency error f error is calculated by The APC control method determines the WF power demand P demand WF for adjusting the grid frequency through input error . error is a nonlinear differential approximation. KP, KI and KD are calculated by The constants bP, bI, and bD represent a conventional PID controller that provides a good but not optimum system response. These constants can be obtained using any conventional methods for PID controllers. In (2), (3), and (4), the coefficients aP, aI, and aD are obtained according to simulation results to determine the relevant ranges of variations for bP, bI, and bD, respectively. The parameters of P, I, and D onlines are determined by the fuzzy rules in Figure 3. Membership functions for inputs and outputs are shown in Figure 3. The vertical axes represent the degree of membership that is within the range of [0, +1]. The meaning of the linguistic variables is listed in Table 1. The input parameters are blurred and then defuzzified, the explicit control signal obtained after defuzzification is shown in Figure 4. Note that the inputs error and error of fuzzy rules need to be The APC control method determines the WF power demand P WF demand for adjusting the grid frequency through input f error .
. f error is a nonlinear differential approximation. K P , K I and K D are calculated by K P (P) = a P P + b P (2) The constants b P , b I , and b D represent a conventional PID controller that provides a good but not optimum system response. These constants can be obtained using any conventional methods for PID controllers. In (2), (3), and (4), the coefficients a P , a I , and a D are obtained according to simulation results to determine the relevant ranges of variations for b P , b I , and b D , respectively. The parameters of P, I, and D onlines are determined by the fuzzy rules in Figure 3. Membership functions for inputs and outputs are shown in Figure 3. The vertical axes represent the degree of membership that is within the range of [0, +1]. The APC control method determines the WF power demand P demand WF for adjusting the grid frequency through input error . error is a nonlinear differential approximation. KP, KI and KD are calculated by The constants bP, bI, and bD represent a conventional PID controller that provides a good but not optimum system response. These constants can be obtained using any conventional methods for PID controllers. In (2), (3), and (4), the coefficients aP, aI, and aD are obtained according to simulation results to determine the relevant ranges of variations for bP, bI, and bD, respectively. The parameters of P, I, and D onlines are determined by the fuzzy rules in Figure 3. Membership functions for inputs and outputs are shown in Figure 3. The vertical axes represent the degree of membership that is within the range of [0, +1]. The meaning of the linguistic variables is listed in Table 1. The input parameters are blurred and then defuzzified, the explicit control signal obtained after defuzzification is shown in Figure 4. Note that the inputs error and error of fuzzy rules need to be The meaning of the linguistic variables is listed in Table 1. The input parameters are blurred and then defuzzified, the explicit control signal obtained after defuzzification is shown in Figure 4. Note that the inputs f error and . f error of fuzzy rules need to be normalized within the range [−1, +1] before being processing in Figure 4, where e represents f error ; de represents . f error . normalized within the range [−1, +1] before being processing in Figure 4, where e represents error ; de represents error . (a)

PRD Method for WT Based on Fatigue Load Sensitivity Using Quadratic Programming Algorithm
In the wind farm power reference dispatch section, the power reference dispatch controller regulates all the WT active outputs of the wind farm. The controller aims to minimize the fatigue load of WT and track the active power allocated by the frequency adjustment controller. Typically, the sampling time of the wind farm controller is in seconds [28]. Therefore, the rapid dynamics of generators and pitch actuators can be ignored [37]. In addition, shaft torsion and tower point oscillations are ignored to reduce the complexity of the model.
The fatigue loads of WTs can be divided into two parts: one is aerodynamic loads and gravity loads (external), and the other is structural loads (internal) [38]. In this paper, the fatigue loads mainly focus on the loads of the drive train due to the torsion of the shaft and the loads of the tower structure due to the tower deflection. Compared with static loads, the dynamic stress causing structural damage of WTs is a much bigger issue. By reducing the fluctuations of low-speed shaft torque Ts and thrust force Ft, the related fatigue loads can be reduced.

PRD Method for WT Based on Fatigue Load Sensitivity Using Quadratic Programming Algorithm
In the wind farm power reference dispatch section, the power reference dispatch controller regulates all the WT active outputs of the wind farm. The controller aims to minimize the fatigue load of WT and track the active power allocated by the frequency adjustment controller. Typically, the sampling time of the wind farm controller is in seconds [28]. Therefore, the rapid dynamics of generators and pitch actuators can be ignored [37]. In addition, shaft torsion and tower point oscillations are ignored to reduce the complexity of the model.
The fatigue loads of WTs can be divided into two parts: one is aerodynamic loads and gravity loads (external), and the other is structural loads (internal) [38]. In this paper, the fatigue loads mainly focus on the loads of the drive train due to the torsion of the shaft and the loads of the tower structure due to the tower deflection. Compared with static loads, the dynamic stress causing structural damage of WTs is a much bigger issue. By reducing the fluctuations of low-speed shaft torque T s and thrust force F t , the related fatigue loads can be reduced.

Improved Model of Fatigue Load Sensitivity
The WT (NREL 5 MW) developed by the National Renewable Energy Laboratory (NREL) is used in this paper [39,40]. The oscillations in the shaft torsion and tower nodding are disregarded, the fluctuations of wind speed are ignored based on the previous report [30]. In order to better optimize the calculation, the equations of fatigue load sensitivity are re-derived. In the process of the optimization dispatch, the redefined fatigue load sensitivity will be used. The active power of the WT output is still controlled by adjusting the pitch angle and torque.
The equivalent mass J t of drive train system is described by [41] where J r is the rotor mass; J g is the generator mass; η g is the gear box ratio. The low-shaft motion equation is described by where ω r is the measured rotor speed; T rot is the aerodynamic torque; T g is the generator torque; The measured generator speed ω g is filtered by a low-pass filter and the filtered speed ω f is where τ f is time constant of the filter of ω g . According to the deviation of ω f from generator rated speed ω g-rated , pitch angle reference θ ref can be obtained by the PI controller.
where k p is the proportional gain; k i is the integral gain; k a is a function of θ ref , defined by k a k a1 + k a2 θ ref .
Where k a1 and k a2 are the constants. By defining β k a θ ref According to the motion equation of shaft torque T s , . ω r , and . ω g can be described by .
where ω g is the measured generator speed; B is the main shaft viscous friction coefficient. According to (11) and (12), . ω g can be expressed by .
The time of the operating point is assumed to be k. The wind speed v is a variable that can be estimated or measured [42]. In this study, v is estimated. The value at t = k is v 0 and is assumed to be constant over a short control period. The measured power output, generator speed, filtered speed, and pitch angle are defined as P g0 , ω g0 , ω f0 , and θ 0 at t = k, respectively. The T rot and T g at t = k can be defined as T rot0 and T g0 , respectively.
Based on (6), (7), (10), and (13), the incremental form can be obtained as ∆ The aerodynamic torque T rot is calculated by where R is the length of the blade; ρ is the air density; v is the wind speed on the rotor; C p is the power coefficient; λ is the tip speed ratio, defined by λ ω r R v ; In order to simplify the expression, P sim is defined by P sim 0.5πR 2 ρv 3 .
According to (18), ∆T rot can be calculated by where C p is described in a lookup table derived from the inputs λ and θ, as is shown in Table 2. Where n and m are the corresponding rows and columns, respectively. The generator torque reference T g_ref is filtered by a low-pass filter and the generator torque T g is derived by, where τ g is the time constant of the filter of T g_ref . T g_ref is calculated by According to (23) and (24), ∆T g and ∆T g_ref can be calculated by (26) where ∆t is the control cycle. According to (14)- (17), the continuous state space model for WT is formulated as where x ≈ [∆ω g , ∆β, ∆ω f , ∆ω r , and the state space matrices are These matrices change every other dispatch cycle. Then, the continuous state space model is discretized with the sampling period t s , which is Substituting (6) and (11), the shaft torque T s can be calculated by Accordingly, Based on (19) and (25), (30) can be transformed into with Based on (28) and (31) ∆T s (k Hence In order to simplify the expressions, a i T s and b i T s are defined by Therefore, the fatigue load sensitivity of the drive train can be expressed as According to [43], The tower dynamics is not included in the simplified WT model. According to, it is assumed the tower base overturning moment M t can be approximately derived by where H is the tower height. F t is thrust force. The thrust force F t is calculated by where C t is the thrust coefficient. Accordingly Similar to (31), Based on (28) and (41) Hence Similar to (34) and (35), Therefore, the fatigue load sensitivity of the tower structure can be expressed as The fatigue load sensitivities considered in the paper are dynamic loads causing structural damage. The equation shown in (33) and (42) are not equations for calculating a fatigue load, but equations for calculating the change in the shaft torque ∆M s and the tower bending moment ∆M t associated with the fatigue load. The parameters in the equation are changing at different times. By reducing ∆M s and ∆M t , the corresponding fatigue load can be reduced. The calculations are carried out by this law of effect. Reductions on the changes of the moments ∆M s and ∆M t are highly corelated to reductions in the damage equivalent fatigue load [44]. ∆M s and ∆M t are related to changes in power reference. ∆M s and ∆M t can be reduced by properly distributing the active power. ∆M s and ∆M t are reduced for each sampling period, so the fatigue loads are reduced.

Cost Function and Constraints
In Figure 1, every individual WT is equipped with an exclusive control system that can follow the power references provided by the frequency adjustment controller.
The controller minimizes the variation of shaft torque T s and thrust force F t to reduce the fatigue load. Accordingly, the cost function is expressed as For the convenience of calculation, the equivalent calculation formula is expressed as where ξ is the weight coefficient.
where P WF demand is the demanded power of WF; P WTi available is the maximum available power of WT-i, it can be estimated by where P i rated is the rated power of WT; v nac is the nacelle wind speed of WT. The corresponding data of C p and C t for this study can be accessed in the wind turbine model of SimWindFarm. The plot of C p (λ, θ) and C t (λ, θ) based on the lookup table is shown in Appendix A.
The presented problems can be expressed as standard quadratic programming (QP) issues [45,46]. It can be effectively solved by a commercial solver. This optimization problem can be solved by different optimization methods. Group intelligent algorithms such as particle swarm optimization (PSO) and genetic algorithm (GA) have good search performance for solving complex problems, but the calculation time is long and is not suitable for real-time online optimization [47][48][49][50]. The programming algorithm has a good ability to search for non-complex solving problems, and the calculation time is short, which is suitable for online optimization. The model proposed in this paper belongs to the non-complex online solution model. Then, the matrix H and the matrix f of QP can be expressed as If the available power of WF does not meet the demand, the following procedure is available.

System Setup
The WF control architecture illustrated in Figure 1 is used to test the performance of control. Figure 1 shows the typical configuration of a WF, which consists of 10, 5 MW WTs. The simulation model is based on SimWindFarm [40] developed in the EU-FP7 project by AEOLUS. The SimWindFarm model allows the real time simulation of flows in the wind field and includes the main aerodynamic effects of the wind farm in MATLAB/Simulink. It consists of four elementary components: wind turbine dynamics, wind field interactions and dynamics, wind farm controller, and electrical network operator. It can represent the main dynamics of a wind farm. The model toolbox includes a wind field generator where mean wind speed, turbulence intensity, and grid resolution can be specified. Furthermore, a number of post-processing tools are included for performance, such as a rain-flow counting algorithm and wake animation. Simulations are conducted for WF over 200 s of run time.
In the test system, the frequency adjustment controller of WF sends instructions every 0.5 s. 0.5 s is to match the wind farm dispatch algorithm and faster to achieve frequency regulation. The parameters of frequency adjustment controller are shown in Table 3. The PRD optimization (OPT) controllers send instructions every 0.5 s. The detailed WT parameters are listed in Table 4. The typical turbulent winds with different intensities generated by the Aeolus toolbox with same average wind speed of 12 m/s are shown in Figure 5. aerodynamic effects of the wind farm in MATLAB/Simulink. It consists of four elementary components: wind turbine dynamics, wind field interactions and dynamics, wind farm controller, and electrical network operator. It can represent the main dynamics of a wind farm. The model toolbox includes a wind field generator where mean wind speed, turbulence intensity, and grid resolution can be specified. Furthermore, a number of post-processing tools are included for performance, such as a rain-flow counting algorithm and wake animation. Simulations are conducted for WF over 200 s of run time.
In the test system, the frequency adjustment controller of WF sends instructions every 0.5 s. 0.5 s is to match the wind farm dispatch algorithm and faster to achieve frequency regulation. The parameters of frequency adjustment controller are shown in Table 3. The PRD optimization (OPT) controllers send instructions every 0.5 s. The detailed WT parameters are listed in Table 4. The typical turbulent winds with different intensities generated by the Aeolus toolbox with same average wind speed of 12 m/s are shown in Figure 5.      Table 4. Parameters of WT.

Parameter Value
Rotor inertia: J r 3.54 × 10 7 (kg·m 2 ) Generator inertia: J g 5.34 × 10 2 (kg·m 2 ) Gear box ratio: Through post-processing, the fatigue cycles based on the rainflow counting method are derived to evaluate the performance of the proposed method. Furthermore, the damage equivalent load (DEL) is based on the Miner's rule and depends on the material properties specified by the slope of the S-N curve for quantizing the load minimization. In this study, the relevant calculations are performed by MCrunch developed by NREL [51].

Performance for the Improved Model of Fatigue Load Sensitivity
In the simple WT model, the drive train can be considered to be rigid if the fatigue load is ignored. However, the drive train is a flexible system. A comparison of measured ω r and ω g η g is shown in Figure 6. As shown in Figure 6, ω r ω g η g . So ∆ω g (k) and ∆ω r (k) are calculated separately in the improved models, and the calculation formula of ω g = η g ω r is not used.
Through post-processing, the fatigue cycles based on the rainflow counting method are derived to evaluate the performance of the proposed method. Furthermore, the damage equivalent load (DEL) is based on the Miner's rule and depends on the material properties specified by the slope of the S-N curve for quantizing the load minimization. In this study, the relevant calculations are performed by MCrunch developed by NREL [51].

Performance for the Improved Model of Fatigue Load Sensitivity
In the simple WT model, the drive train can be considered to be rigid if the fatigue load is ignored. However, the drive train is a flexible system. A comparison of measured ωr and is shown in Figure 6. As shown in Figure 6, ωr ≠ . So Δωg(k) and Δωr(k) are calculated separately in the improved models, and the calculation formula of ω g = η g ω r is not used.  Figures 7-10a, which will not be described in detail in the following paper. The wind speed of the WF is 12 m/s with 0.1 turbulence intensity. The values of ∆ω g (k + 1) and ∆ω r (k + 1) obtained by calculating WT1 are shown in Figures 7 and 8. Note that x(k + 1) is obtained by discretizing x(k). Figures 7, 8, 9 and 10b are enlarged image of the inside of the dashed box in Figures 7, 8, 9 and 10a, which will not be described in detail in the following paper.
Through post-processing, the fatigue cycles based on the rainflow counting method are derived to evaluate the performance of the proposed method. Furthermore, the damage equivalent load (DEL) is based on the Miner's rule and depends on the material properties specified by the slope of the S-N curve for quantizing the load minimization. In this study, the relevant calculations are performed by MCrunch developed by NREL [51].

Performance for the Improved Model of Fatigue Load Sensitivity
In the simple WT model, the drive train can be considered to be rigid if the fatigue load is ignored. However, the drive train is a flexible system. A comparison of measured ωr and is shown in Figure 6. As shown in Figure 6, ωr ≠ . So Δωg(k) and Δωr(k) are calculated separately in the improved models, and the calculation formula of ω g = η g ω r is not used.  Figures 7-10a, which will not be described in detail in the following paper. Through post-processing, the fatigue cycles based on the rainflow counting method are derived to evaluate the performance of the proposed method. Furthermore, the damage equivalent load (DEL) is based on the Miner's rule and depends on the material properties specified by the slope of the S-N curve for quantizing the load minimization. In this study, the relevant calculations are performed by MCrunch developed by NREL [51].

Performance for the Improved Model of Fatigue Load Sensitivity
In the simple WT model, the drive train can be considered to be rigid if the fatigue load is ignored. However, the drive train is a flexible system. A comparison of measured ωr and is shown in Figure 6. As shown in Figure 6, ωr ≠ . So Δωg(k) and Δωr(k) are calculated separately in the improved models, and the calculation formula of ω g = η g ω r is not used.  Figures 7-10a, which will not be described in detail in the following paper. As is shown in Figures 7b and 8b, The values of ∆ω g and ∆ω r after 0.5 s can be accurately calculated using the improved model. A better basis for optimizing results is provided by this precise calculation.
The effect of changes in wind speed is considered when calculating ∆T rot . Comparison of the case of adding ∂T rot ∂v ∆v and not adding is shown in Figure 9.
As is shown in Figures 7b and 8b, The values of Δωg and Δωr after 0.5 s can be accurately calculated using the improved model. A better basis for optimizing results is provided by this precise calculation.
The effect of changes in wind speed is considered when calculating ΔTrot. Comparison of the case of adding ∂T rot ∂v ∆v and not adding is shown in Figure 9. The model is simulated based on SimWindFarm. Since the control time is 0.5 s, sampling is performed every 0.5 s. Because Trot is cube of v, even small changes of v can have a big impact on ΔTrot. As is shown in Figure 9b, in the vicinity of 154 s, there is a big step in the measured value, which is caused by the change in wind speed. Mutations in ΔTrot caused by sudden changes in wind speed can still be accurately calculated by the improved model. However, the original models that did not consider ∂T rot ∂v ∆v cannot adapt to this mutation. The value of ΔTrot after 0.5 s is accurately calculated.
After a detailed improvement, the ΔFt calculation is more accurate. The comparison between the improved value and the measured value is shown in Figure 10. As is shown in Figure 10b. In the vicinity of 82 s, there is a big step in the measured value, which is caused by the change in wind speed. The improved model can better adapt to the ΔFt mutation like ΔTrot. This improvement makes the calculation of the optimization target more accurate and the optimization effect is better. The re-derived fatigue load sensitivity equations can improve the calculation accuracy very well.

Performance for Different Turbulence Intensity
The normal PRD controller uses a proportional distribution method (NORM). The reference power is assigned to the WT in the available power ratio. As is calculated by The model is simulated based on SimWindFarm. Since the control time is 0.5 s, sampling is performed every 0.5 s. Because T rot is cube of v, even small changes of v can have a big impact on ∆T rot . As is shown in Figure 9b, in the vicinity of 154 s, there is a big step in the measured value, which is caused by the change in wind speed. Mutations in ∆T rot caused by sudden changes in wind speed can still be accurately calculated by the improved model. However, the original models that did not consider ∂T rot ∂v ∆v cannot adapt to this mutation. The value of ∆T rot after 0.5 s is accurately calculated. After a detailed improvement, the ∆F t calculation is more accurate. The comparison between the improved value and the measured value is shown in Figure 10. As is shown in Figures 7b and 8b, The values of Δωg and Δωr after 0.5 s can be accurately calculated using the improved model. A better basis for optimizing results is provided by this precise calculation.
The effect of changes in wind speed is considered when calculating ΔTrot. Comparison of the case of adding ∂T rot ∂v ∆v and not adding is shown in Figure 9. The model is simulated based on SimWindFarm. Since the control time is 0.5 s, sampling is performed every 0.5 s. Because Trot is cube of v, even small changes of v can have a big impact on ΔTrot. As is shown in Figure 9b, in the vicinity of 154 s, there is a big step in the measured value, which is caused by the change in wind speed. Mutations in ΔTrot caused by sudden changes in wind speed can still be accurately calculated by the improved model. However, the original models that did not consider ∂T rot ∂v ∆v cannot adapt to this mutation. The value of ΔTrot after 0.5 s is accurately calculated.
After a detailed improvement, the ΔFt calculation is more accurate. The comparison between the improved value and the measured value is shown in Figure 10. As is shown in Figure 10b. In the vicinity of 82 s, there is a big step in the measured value, which is caused by the change in wind speed. The improved model can better adapt to the ΔFt mutation like ΔTrot. This improvement makes the calculation of the optimization target more accurate and the optimization effect is better. The re-derived fatigue load sensitivity equations can improve the calculation accuracy very well.

Performance for Different Turbulence Intensity
The normal PRD controller uses a proportional distribution method (NORM). The reference power is assigned to the WT in the available power ratio. As is calculated by As is shown in Figure 10b. In the vicinity of 82 s, there is a big step in the measured value, which is caused by the change in wind speed. The improved model can better adapt to the ∆F t mutation like ∆T rot . This improvement makes the calculation of the optimization target more accurate and the optimization effect is better. The re-derived fatigue load sensitivity equations can improve the calculation accuracy very well.

Performance for Different Turbulence Intensity
The normal PRD controller uses a proportional distribution method (NORM). The reference power is assigned to the WT in the available power ratio. As is calculated by In order to test the efficacy of the improved model, the operation of wind farm under different turbulence intensity were studied and the FAC is not used. The average wind speed of the test wind is 12 m/s, and the turbulence intensity is 0.1, 0.2, and 0.3 respectively. The results of OPT and NORM are compared under different turbulence intensities. The tower bending moment was mainly considered in this paper, so the weight coefficient ξ to 3000 is selected. Scenario 1: The turbulence intensity is 0.1. The wind farm reference power P WF demand is 35 MW under this scenario. The calculated DELs of T s and M t for all WTs are listed in Table 5. As shown in Table 5, some of the WTs have more growth in DELs for T s , but this growth will be offset as the simulation time increases. The total DELs of wind farms using OPT is still lower than NORM. From the whole WF point of view, the reductions of the total DEL for T s is 2.21%.   Table 6. As shown in Table 6, the DELs of both T s and M t with the OPT are reduced compared with NORM. When the turbulence intensity is 0.2, a more detailed cumulative rainflow cycle is shown in Figure 11. A representative WT (WT09) is used as an example. Compared with the NORM, less cycles are found for OPT, which implies less fatigue loads experienced by the WT.
When the turbulence intensity is 0.3, a more detailed cumulative rainflow cycle is shown in Figure 12. A representative WT02 is used as an example. Compared with the NORM, less cycles are found for OPT, which implies less fatigue loads experienced by the WT.
Fatigue load optimization results with three different turbulence intensities were obtained. It can be seen from the results that the improved model can optimize the fatigue load well under different turbulence intensities.

WF Controller Performance
In order to prove the performance of the frequency control of the developed solution, the wind farm is assumed to be able to get enough power from the wind. When the WF is operating in the frequency regulation mode, the measured grid frequency is used as a feedback signal to establish real-time active power control. The baseline model for electrical network operator in SimWindFarm can function in different power command modes such as: absolute, delta, and frequency regulation modes. Basically, baseline control scheme in frequency regulation mode can maintain the necessary balance between power generation and loads. The demand power of the baseline control is calculated by where t d and t c are two constants (t c > t d ) of dead bands and control bands, defined by t d = 0.002 and t c = 0.02 which can enable the baseline control to function optimally. Moreover, P up and P down are defined by P up = P max + P min (58) where P min and P max are the prescribed minimum and maximum limits for the total power of WF. In order to verify the performance of the control method, different functions of the grid load are defined and used for WF testing. Baseline control is compared to Fuzzy-PID. NORM dispatch method is compared to OPT method. The initial frequency is assumed to be 50 Hz. The test wind is 12 m/s with 0.3 turbulence intensity. Figures 13a and 14a show the performance of the active power control method for the WF in response to the different loads considered, respectively. Figures 13b and 14b show the grid frequency responses with respect to the different loads, respectively. The power and frequency responses of Fuzzy-PID and baseline control are compared in Figures 13 and 14, respectively. The power curve of the wind farm output is the same whether the WF dispatch method uses NORM or OPT.
where td and tc are two constants (tc > td) of dead bands and control bands, defined by td = 0.002 and tc = 0.02 which can enable the baseline control to function optimally. where Pmin and Pmax are the prescribed minimum and maximum limits for the total power of WF. In order to verify the performance of the control method, different functions of the grid load are defined and used for WF testing. Baseline control is compared to Fuzzy-PID. NORM dispatch method is compared to OPT method. The initial frequency is assumed to be 50 Hz. The test wind is 12 m/s with 0.3 turbulence intensity. Figures 13a and 14a show the performance of the active power control method for the WF in response to the different loads considered, respectively. Figures 13b and 14b show the grid frequency responses with respect to the different loads, respectively. The power and frequency responses of Fuzzy-PID and baseline control are compared in Figures 13 and 14, respectively. The power curve of the wind farm output is the same whether the WF dispatch method uses NORM or OPT.  where Pmin and Pmax are the prescribed minimum and maximum limits for the total power of WF. In order to verify the performance of the control method, different functions of the grid load are defined and used for WF testing. Baseline control is compared to Fuzzy-PID. NORM dispatch method is compared to OPT method. The initial frequency is assumed to be 50 Hz. The test wind is 12 m/s with 0.3 turbulence intensity. Figures 13a and 14a show the performance of the active power control method for the WF in response to the different loads considered, respectively. Figures 13b and 14b show the grid frequency responses with respect to the different loads, respectively. The power and frequency responses of Fuzzy-PID and baseline control are compared in Figures 13 and 14, respectively. The power curve of the wind farm output is the same whether the WF dispatch method uses NORM or OPT.  As shown in Figure 13a, the grid load A is a step variable. As is shown in Figure 13b, the baseline control cannot achieve a frequency of 50 Hz, and there is still an error of 0.004 Hz as the control time increases. The frequency error of the Fuzzy-PID control method is smaller than that of the baseline control. As the control time increases, the frequency gradually approaches 50 Hz and the error disappears. The Fuzzy-PID control method can quickly track the frequency so that the error is approximately zero at the 100 s. Therefore, the Fuzzy-PID method has better control effects than the baseline control.
As is shown in Figure 14a, grid load B is the varying load. Baseline control failed to track the grid load during the initial phase of 20-50 s, which also caused a relatively large frequency fluctuation of 20-50 s as shown in Figure 14b. Fuzzy-PID control can track grid load very well by compared to baseline control. As is shown in Figure 14b, Fuzzy-PID can track the frequency very well, and the frequency error of Fuzzy-PID is smaller. As are shown in Figures 13 and 14, less power fluctuations make power tracking better.
To check the tracking accuracy of APC methods, the normalized root-mean-squared error (NRMSE) is defined as To illustrate the ability of Fuzzy-PID methods in the regulation of the grid frequency, Table 7 quantitatively compares the obtained grid frequency responses in terms of standard deviation (SD) and NRMSE of power responses. For a step grid load like A, compared to baseline control, NRMSE for power responses increased by 22.02%, and SD for grid frequency responses is reduced by 88.45%. For a varying grid load like B, compared to baseline control, NRMSE for power responses is reduced by 73.15% and SD for grid frequency responses is reduced by 81.88%.
By testing the control method of the WF with different loads, the Fuzzy-PID control has a good performance regardless of the speed or error of the tracking frequency.

Fatigue Loads Performance
The fatigue load results of the WTs using Baseline Control + NORM, Fuzzy-PID + NORM and Fuzzy-PID + OPT are shown in Figure 15 when load B is used to test the overall controller. baseline control.
As is shown in Figure 14a, grid load B is the varying load. Baseline control failed to track the grid load during the initial phase of 20-50 s, which also caused a relatively large frequency fluctuation of 20-50 s as shown in Figure 14b. Fuzzy-PID control can track grid load very well by compared to baseline control. As is shown in Figure 14b, Fuzzy-PID can track the frequency very well, and the frequency error of Fuzzy-PID is smaller. As are shown in Figure 13 and 14, less power fluctuations make power tracking better.
To check the tracking accuracy of APC methods, the normalized root-mean-squared error (NRMSE) is defined as To illustrate the ability of Fuzzy-PID methods in the regulation of the grid frequency, Table 7 quantitatively compares the obtained grid frequency responses in terms of standard deviation (SD) and NRMSE of power responses. By testing the control method of the WF with different loads, the Fuzzy-PID control has a good performance regardless of the speed or error of the tracking frequency.

Fatigue Loads Performance
The fatigue load results of the WTs using Baseline Control + NORM, Fuzzy-PID + NORM and Fuzzy-PID + OPT are shown in Figure 15 when load B is used to test the overall controller. As is shown in Figure 15, the frequency adjustment effect by the Fuzzy-PID is better, but the DELs of T s is significantly larger than Baseline Control + NORM. The specific results are shown in Table 7, the fatigue load using Fuzzy-PID + NORM increased by 6.47% compared with the fatigue load using Baseline Control + NORM. The fatigue load using Fuzzy-PID + OPT still increased by 2.44% compared with Baseline Control + NORM, but decreased by 3.79% compared with Fuzzy-PID + NORM.
The fatigue load of M t using Fuzzy-PID + NORM was significantly higher than that of Baseline Control + NORM, and the specific results were increased by 8.15% as shown in Table 8. The fatigue load using Fuzzy-PID + OPT is 0.32% higher than Baseline Control + NORM and 7.24% lower than Fuzzy-PID + NORM. The fatigue load of the WT will increase since the WF needs to frequently change the pitch angle and reference power when participating in frequency regulation. The better the frequency adjustment effect, the more frequent this change, the greater the fatigue load caused. The overall fatigue load of the WT is reduced, DELs of M t is hardly affected especially after considering the fatigue load.

Conclusions
In this paper, an APC method for supporting grid frequency regulation in WFs considering fatigue load of WT is proposed. The Fuzzy-FID control is used for a frequency adjustment controller to maintain the balance between power generation and grid load. The equations of fatigue load sensitivity are re-derived. The calculation accuracy is improved by the re-derived equations. The re-derived equations are used for PRD to minimize the fatigue load for all the WTs. Case studies show that the frequency adjustment control method based on Fuzzy-PID can respond to and recover grid frequency deviations more quickly. Less power fluctuations make power tracking better. Compared to baseline control, NRMSE for power responses and SD for grid frequency of Fuzzy-PID are reduced by 81.45% and 81.88%, separately. In addition, the explicit analytical equations of fatigue load sensitivity are re-derived to improve the calculation accuracy of the fatigue loads. DELs for T s reduced by 2.21% to 5.13%, DELs for M t reduced by 5.32% to 7.32%. This improvement minimizes fatigue loads for wind farms under different turbulent wind and different loads. The proposed solution is suitable for real-time control of large-scale WF.
Since this study is mainly focused on the wind farm control level, some of the dynamics of WTs are ignored. Rainflow counting is applied to the time histories of thrust and torque from the simulation results, so the evaluation of fatigue load cannot be fully reflected. Quantitative conclusions are uncertain and can only represent trends in fatigue load changes. Since this study is a simulation based on SimWindFarm, the study cannot consider all practical engineering problems. Verifying this method in engineering is what we will do next.
Author Contributions: The paper was a collaborative effort among the authors. Y.W. carried out relevant theoretical research, performed the simulation, analyzed the data, and wrote the paper. Y.L. and X.W. provided critical comments. J.Z. and W.H.L. revised the paper. Grid frequency error f meas Grid measured frequency f N Normal frequency of the grid J t Equivalent mass of the drive-train J r Rotational inertia of the rotor J g Rotational inertia of the generator η g Gear box ratio ω r Measured rotor speed ω g Measured generator speed ω g-rated Rated generator speed F t Thrust force T rot Aerodynamic torque T g Generator torque T g_ref Generator torque reference Constants of k a B Main shaft viscous friction coefficient P g0 Output power of a turbine at t = k ω g0 Generator speed at t = k ω f0 Filtered speed of the generator speed at t = k θ 0 Pitch angle at t = k T rot0 Aerodynamic torque at t = k T g Generator torque at t = k R Length of the blade H Tower height ρ Air density v Wind speed of hub height C p Power coefficient C t Thrust coefficient λ Tip speed ratio