Analytical system frequency response model with virtual synchronous wind turbines

The present paper devises a second-order system frequency response (SFR) model for power systems integrating virtual synchronous variable speed wind turbines (VSWTs). The inertial response model of full converter VSWTs, in which grid-side converters are controlled using a virtual synchronous generator concept, is developed. The derived model is incorporated into the traditional SFR model deployed to estimate the frequency behaviour of the synchronous generators following infeed loss contingencies. The proposed SFR model formulates interaction between the VSWTs rotor speed reduction and the system frequency drop. The system loads inertia is also included in the SFR model. Addition-ally, an analytical approach to derive the suggested SFR model for the large power system including VSWTs with different capacities and operating points is presented. Solving the proposed model gives the analytical functions for estimating deviations in the system frequency as well as VSWTs rotor speed variations. Finally, the accuracy of the devised SFR model is veriﬁed through the time-domain simulations of a modiﬁed IEEE 39-bus system. Particularly, impact of the VSWTs virtual inertia on the efﬁciency of the proposed SFR model is studied.


Literature reviews
For traditional power systems, a low-order SFR model is presented in [9]. Whereas this is not the most detailed model, it provides a good enough estimate of the frequency variations of the power systems incorporating SGs. Following the extended SFR model, an analytical method is proposed for aggregating multi-machine network into a single-machine SFR model [10]. On the other hand, several researches have been conducted to integrate dynamics of the VSWTs emulated inertial response into SFR models. For example, a modified SFR model considering inertial responses from de-loaded VSWTs is developed in [11]. Additionally, a simplified linearized model for Type 3 wind turbines (WTs) with a typical emulated inertia [12] is designed in [13] to unify with power systems SFR model. In some studies, the frequency response model of VSWTs is derived by taking different WTs operating conditions into consideration [14,15]. Further, a system identification technique is proposed to aggregate multi-machine VSWTs networks [16]. A generic SFR model is also suggested in [17] to incorporate the dynamics of hydro and wind generation units into the traditional SFR model. In this generic SFR model, a uniform transfer function is employed to describe the dynamics of the aggregated governors [17]. A new methodology is presented in [18] in order to integrate the doubly fed induction generator's (DFIG's) frequency response into the SFR model comprising of hydro and thermal power plants.
In [19], an analytical frequency model for power system including VSWTs is used for tuning typical inertia emulators. However, the frequency model is derived for a wind farm consisting of identical WTs capacities as well as similar operating points [19]. Other researches derived the SFR models of the power systems including CIGs to incorporate into frequency constrained unit commitment problems [20][21][22]. Specifically, the deviations in the system frequency [20,21] along with that of the VSWTs rotor speed [22] are formulated.
In the SFR models introduced in [11] and [13][14][15][16][17][18][19], VSWTs are not controlled using the VSG scheme. However, an emerging trend is to control VSWTs by the VSG concept. In this regard, several strategies have been proposed to control gridside converters (GSCs) in Type 4 WTs using the VSG scheme [23][24][25]. The SFR models for the power systems incorporating VSG-controlled CIGs are developed in [26,27]. However, these studies did not model the impact of the system frequency variations on dynamics of the CIGs energy sources. In other words, it is assumed that the CIGs have infinite energy resources. In the case of VSWTs, a detailed model of VSG-based Type 4 WT is linearized to derive its frequency response model [28]. In this approach, the mechanical power of the VSG is generated by a proportional-integrator (PI) controller. The input of the PI controller is the WT rotor speed deviation which in turn is generated by a maximum power point tracking (MPPT) function [28]. Whilst the VSG virtual inertia and damping coefficient are expressed analytically in a frequency domain, the derived VSWT linearized model cannot be easily incorporated in the large power systems SFR model. Further, a study is conducted in [29] to evaluate the frequency response of a VSG-controlled VSWT, where the VSG mechanical power is directly emanated from an MPPT function. It is shown that these strategies ensure the secure operation for WTs rotor speed during the grid frequency support time interval. However, it is not illustrated how the proposed techniques can be integrated into the traditional SFR model. Furthermore, there is no analytical solution for the SFR model derived in [29].

Procedures and contributions
In order to satisfy the frequency requirements, there has been some literature focusing on incorporating the system frequency constraints in the modern power system operation and planning studies [20][21][22]27]. On the other hand, the size and type of WTs installed in the real-world large-scale power systems varied The control scheme of the virtual synchronous controlled VSWT significantly [30]. In this study, a novel analytical SFR model for the modern power systems including virtual synchronous VSWTs is presented. To this end, a control scheme is devised for the full converter WTs in which GSCs are controlled using the VSG concept. The proposed SFR model formulations are set forth considering WTs with different nominal power and operating points. Therefore, the following contributions are made in this study: • A second-order SFR model is developed by formulating the interaction between VSWTs rotor speed and grid frequency following infeed loss contingencies. • The analytical functions are obtained for estimating deviations in VSWTs rotor speed in addition to the system frequency variations. • The system loads inertia is incorporated in the proposed SFR model. • The devised SFR model is verified investigating various installed capacities, different operating points and virtual inertia for VSG-based VSWTs.
Accuracy of the suggested model in estimating the system frequency dynamics is investigated by time-domain simulations on the IEEE 39-bus system. The VSWTs size and distribution on the test system are selected according to actual wind energy data at Europe [30]. In this regard, the VSWTs installed capacity ranged from 10% to 70% of the system demand considered. Additionally, the VSWTs operating points from 20% to 100% of their nominal power are analysed. This paper is organized as follows: The developed model for VSG-based WTs is concisely described in Section 2. In Section 3, the proposed SFR model is explained in detail. The simulation results are provided in Section 4. Finally, the paper is concluded in Section 5.

VIRTUAL SYNCHRONOUS CONTROLLED VARIABLE SPEED WIND TURBINES MODELLING
Assume that there are m wind farms in the studied power system. In practice, a farm may be constructed by a number of VSWTs. However, each of the farm can be considered to be an aggregation of multiple identical VSWTs [31]. The control scheme of the simulated model for the jth VSG-based VSWT is shown in Figure 1. It is helpful to clarify that in the following expressions, all of the power variables are represented in per unit using the system base power P b . The frequency and rotor speed quantities are also in per unit. First, the WT power P w, j is translated to the WT generator (WTG) shaft power P sh, j by a drive train. Next, the shaft power is converted to the WTG electrical power P g, j . In this study, the WTG is a permanent magnet synchronous generator (PMSG). The DC link voltage is regulated using the machine-side converter (MSC), while the VSG power P ve, j is controlled via the VSG-based GSC. The mechanical power P vm, j for the VSG is determined by a MPPT function [25]. The WT power P w, j , is given by [32]: where V w, j is wind speed in meter per second. The WT rotor speed is denoted by w, j . The parameter C P and the variable are the VSWT power coefficient and tip speed ratio, respectively. The parameters k r and k v are [32]: wherein V wb, j denotes the rated wind speed. Parameter C Pmax is maximum of the power coefficient. Parameter cpmax denotes the tip speed ratio pertinent to C Pmax . The base values for the VSWT power and rotor speed are denoted with P wb, j and wb, j , respectively. The maximum power point WT power determines the VSG mechanical power [32]: in which, The MSC controller regulates the WTG q-axis stator current to stabilize DC link voltage [33]. It also deploys a damping controller to damp the WTG rotor speed oscillations [34]. In this study, the GSC controller is a fourth-order VSG model. Specifically, the two orders simulate the VSG rotor dynamics and the two other orders model dynamics of the VSG virtual field and damper windings. It is worth noting that the energy corresponding to the VSGs virtual inertia comes from the kinetic energy stored in the VSWTs rotor. Therefore, these VSWTs experience temporary reductions in their rotor as well as in the VSGs mechanical power during first seconds following an infeed loss event in the power system. This momentary reduction in the VSGs mechanical power increases burden on the PFR service providers. Accordingly, formulating these interactions into the existing SFR models is necessary in order to help the power system survive to frequency instability problems.

SYSTEM FREQUENCY RESPONSE MODEL WITH VIRTUAL SYNCHRONOUS CONTROLLED VSWTS
A SFR model is established in this section for large power systems including the VSG-based WTs.

3.1
System frequency response model's formulations Consider a power system with n + 1 power stations equipped with SGs and m wind farms constructed by the VSG-based VSWTs. Assume that (n + 1)th SG unit is suddenly tripped off from the system. With these assumptions, the SGs and VSGs electromechanical equations can be expressed as: where P sm,i and P se,i denote the mechanical and electrical powers of the ith SG, respectively. The variable s,i represents the SG rotor speed. H s,i is the inertia constant of the SG, in seconds. The base power of the SG is denoted by P sb,i . These definitions are also applicable on (6) representing dynamics of the VSGs. Similar equations for the VSWTs can be written as: 2H w, j P wb, j w, j d dt w, j = P w, j − P g, j j = 1, 2, … , m, where P w, j , P g, j denote the VSWT power and WTG power pertinent to jth farm, respectively. The variable w, j shows the VSWT rotor speed. Further, the parameter P wb, j is the VSWT base power. Additionally, the parameter H w, j is combined with the inertia constant of the VSWT and WTG (see Appendix B). Summing the individual equations of (5) and (6) gives: where H s and H v denote the aggregated inertia of the SGs and VSGs, respectively. Their detailed derivations are provided in Appendix A. The variables s and v represent the center-ofinertia (COI) frequencies of the SGs and VSGs, correspondingly [35]. Accordingly, the system COI frequency sv is defined as follows: with as the system total inertia. Neglecting the inter-area oscillations between the SGs and VSGs, combining (8) and (9) yields: in which the power variables are the aggregated powers (see Appendix A). On the other hand, summation of the individual equations in (7) results in: where H w denotes the combined inertia for all of the VSWTs. The variables w is the COI VSWTs rotor speed as follows: Now, we assume that all VSWTs rotor speeds w, j are equal to w . This is not a good assumption because of variable speed nature of VSWTs. However, it does not affect the proficiency of the proposed SFR model. This is because the SFR model estimates deviation in the VSWTs rotor speed and not their absolute values. With this assumption, (13) may be rewritten as: where the power variables are the aggregated powers of the VSWTs (see Appendix A).

System frequency response model's linearization
Let us presume that the power system is initially operating in an equilibrium point with sv = 1.0. Further, assume that the (n + 1)th SG generates power P lg . Following tripping off this unit at t = t lg , the cumulative variations in powers of the other SGs and VSGs supply the lost power P lg , the system power loss changes and the loads power variations ΔP d . In particular, the loads demand have a voltage-dependent power and a frequencydependent power. However, the former is usually ignored in the SFR models [9]. With these assumptions and by neglecting power loss, (12) can be linearized as: where the prefix Δ denotes a small deviation. The last term in the above equation models the loads power variation with D as the system load damping constant [36].
To solve (16), the deviation in the mechanical power of the SG units ΔP sm , and that of the VSG units ΔP vm , must be known. The former, namely the system primary frequency response (PFR), can be expressed as the following function [37]: wherein R p fr and T p fr denote PFR ramp rate and time constant, respectively. The maximum available PFR reserve power should satisfy the following equality: where Δ svq denotes the allowable change in the quasi-steadystate system frequency sv following the event [37].
In the case of the VSGs, the mechanical power ΔP vm , is a function of deviation in the COI VSWTs rotor speed w . For calculation, the deviation in w , (15) can be linearized as: with w0 as the pre-event value of w . Note that the power changes in DC link capacitors of the VSWTs are ignored in deriving (19). Furthermore, from (6), ΔP ve can be derived as follows: Finally, by replacing ΔP ve in (19) with right-hand side of (20), the second equation of the proposed SFR model is yielded as follows: The block diagram shown in Figure 2 illustrates the secondorder proposed SFR model where its state variables are Δ sv and Δ w . The frequency dynamics of the SGs and VSGs are modelled in its top section. While, the bottom part of Figure 2 simulates the VSWTs rotor dynamics. In the proposed scheme, dynamics between the VSWTs with those of the SGs and VSGs are simulated through the power deviation ΔP vm . This quantity is in turn proportional to the state variable Δ̄w through the gain k m . How this parameter is calculated is explained in Appendix C. Following the power deficit event, ΔP w can be ignored in comparison to other power deviations, assuming VSWTs operate below their rated power. Note that the suggested SFR model will reduce to the traditional SFR model proposed in [9] when the switch position is changed from 1 to 2.

The proposed system frequency response model
The time-domain solutions for the system frequency deviation Δ sv as well as the WTs rotor speed deviation Δ w , can be obtained by solving (16) and (21). Considering (17), the variables Δ sv and Δ w have two components as illustrated in Figure 3. Applying Laplace transform on (16) and (21), and using inverse Laplace transform yields: in which: Similarly, the components of the Δ w can be derived as: where Further, the terms Δ 1 to Δ 4 are exponential functions as follows: The time constants 1 and 2 can be calculated using the following characteristic polynomial: Therefore, time constants 1 and 2 may be expressed as: These time constants may be complex values depending on the amount of wind generation. The coefficients a 1 to a 4 , which are introduced in (26) and (27), are provided in Appendix A.
As observed, the derived functions for estimating the system frequency consists of the exponential functions with two different time constants. While the traditional SFR model has only a time constant [9]. It is worth noting that when the VSWTs installed capacity tends to zero, it can be proved that the proposed SFR model reduces to the traditional model.
We can now derive the analytical expressions for maximum deviation in the system frequency as well as the VSWTs rotor speed. In Figure 3, these quantities are denoted with Δ svm and Δ wm , respectively. When the system frequency deviation reaches its maximum before t ′ reaches T p fr , the time interval t svm can be estimated as: where t 0 is the initial guess for t svm1 . Actually, t 0 is the exact value of the time interval t svm pertinent to the traditional SFR model that can be calculated as follows: In this case, the maximum deviation in the system frequency is given by substituting t ′ with t svm1 into (23) as follows: On the other hand, if the Δ sv reaches its maximum after t ′ reaches T p fr , taking the derivative of the second equation in (23) with respect to t ′′ and setting it equal to zero yields: .
In this case, the maximum deviation in the system frequency is given by substituting t ′′ with t svm2 into (23) as follows: Accordingly, the time interval t svm may be determined as: and the maximum deviation in the system frequency is: Following a similar procedure, when the VSWTs rotor speed deviation reaches its maximum before t ′ reaches T p fr , taking the derivative of the first equation in (24) with respect to t ′ and setting it equal to zero obtains: .
In this case, the maximum absolute value of Δ w is obtained by replacing t ′ with t wm1 into (24) as follows: On the contrary, if the Δ w reaches its maximum after t ′ reaches T p fr , taking the derivative of the second equation in (24) with respect to t ′′ and setting it equal to zero gives: .
In this case, the maximum absolute value of Δ w is obtained by replacing t ′′ with t wm2 into (24) as follows: Accordingly, the time interval t wm may be calculated as: and the maximum absolute value of Δ w is given by:

SIMULATION RESULTS AND DISCUSSIONS
In this section, the studied power system is firstly described. Then, the simulation results are provided to compare the proposed SFR model with the detailed model of the system.

4.1
Description of the studied system The studied system is a modified IEEE 39-bus system with its single line diagram is portrayed in Figure 4 [38]. This system is implemented in DIgSILENT PowerFactory environment. It has 11 SG-based power stations and 9 VSG-based wind farms. All of the SG and VSG units are equipped with IEEE Type 1 exciter without a power system stabilizer [39]. The detailed parameters of the SG units and those of the wind farms are provided in Appendix B. The system demand is 6150 MW that is selected as the system base power P b . Half of the demand is consumed by static loads and the remaining half is used by induction motors [38]. The system load damping constant D is 1.4 in the base of P b . In this study, the largest infeed loss is considered as 5% of P b [40][41][42]. Thus, to contain the quasi-steadystate frequency deviation to 1% of the system base frequency f b ( = 50 Hz) [37], the PFR reserve becomes: P res = P lg + DΔ svq = 0.05 + 1.4 (−0.01) = 0.036p.u. (44) The P res can be delivered within 10 s, that is, T p fr is 10 s [37]. In the PFR controller, the up-rate power limit is R p fr while the down-rate power is set to 10 times of R p fr [39]. The PFR mechanism is conducted by SG10.
The VSGs virtual inertia constant is defined as follows [32]: in which, k h is a factor to control the VSGs inertia. The values of the power P ve0, j the VSWT inertia constant H w, j and the WTG inertia constant H g, j are provided in Appendix B.

System frequency response model verification
In this case, the total wind farms base power P wb is 50% of the system base power. The wind farms individual base power and frequency dynamics are different with each other. Additionally, the various operating points, range from 45% to 100% loadings, are considered for the VSWTs. Further details in this regard are provided in Appendixes B and E.
The estimated SFR to tripping off the generator SG11 is compared to its true value shown in Figures 5-7. The true values are achieved through time-domain simulations of the detailed test system model. The quantities estimated by the traditional [9,20] and the proposed SFR models are labelled by Traditional   Figure 5a, we see that the estimated frequency nadir is 0.1 Hz larger and smaller than its true value using the proposed and the traditional SFR models, correspondingly. Recall that half of the system demand is consumed by the induction motors. Thus, inertial responses of the motors should be incorporated into the SFR model to achieve more accurate estimations. To this end, the cumulative motors inertia constant H m should be added to that of the SGs, that is, to H s . In the studied system, the amount of H m is about 0.625 s, calculated based on the system base power [38]. The updated estimations are shown in Figure 5b. In the proposed SFR model, the The true values of the maximum deviation of the COI frequency versus VSWTs (left) installed capacity and (right) operating point discrepancies between the estimated and true values is significantly vanished by considering the motor loads inertia. However, accuracy of the traditional model is deteriorated when the motors inertia is considered. This observation can also be seen in Figure 6, where the estimated system frequency RoCoF is compared with its true value. The RoCoF values are calculated using the frequency values with 0.1 s time intervals. Finally, it can be observed from Figure 7 that the suggested model gives a superior estimation of the COI VSWTs rotor speed variations. However, the traditional model is not able to estimate this variable. Note that the VSWTs are running in different operating points before the event inception. The detailed simulation results of the VSWTs are provided in Appendix E.
The above simulations results are calculated assuming that only generator SG10 provides PFR service. However, the accuracy of the proposed SFR model is verified in Appendix D considering multiple PFR service providers with distinguished characteristics.

Impact of wind power generation on the SFR model
The true values of the maximum deviation in the system frequency following the tripping off the generator SG11 is illustrated for various wind power generation levels in Figure 8. In particular, the results are provided for the VSWTs installed capacity P wb ranged from 10% to 70% of the system base power. Furthermore, the VSWTs operating points ranged from 20% to 100% of their base powers are considered. In order to gain insight into the relation between Δ svm and the wind generation level, it is plotted in Figure 9a versus the pre-event total wind generations P ve0 . Furthermore, the corresponding estimation errors ignoring and considering the motors inertia are depicted in Figures 9b and 9c, respectively. In both cases, the estimation error of the traditional model increases dramatically with the wind generation level. On the contrary, the proposed model results in 10% estimation error when the loads inertia is neglected regardless of the wind generation. Additionally, the error is reduced by three times by considering the motors inertia in the suggested SFR model. Note that a negative estimation error means the absolute value of the estimated quantity is lower than that of its true value. Further, regarding the time interval t svm estimation, it is observed that the estimation error ranged from −5% to −11% using the traditional model with modelling   Figure 10 illustrates the true values of the maximum deviation in the COI VSWTs rotor speed for the various wind installed capacities and power generations. These values are depicted versus the pre-event wind power generation in Figure 11a. The estimation errors pertinent to Δ wm obtained by the proposed model with considering the loads inertia are shown in   Figure 11b. Accordingly, the estimation errors ranged from −6% to −11% depending on the wind generation level.
The time constants of the exponential functions in the proposed SFR model for various wind installed capacities are presented in Figure 12. For each installed capacity, the time constants are shown for the VSWTs operating points ranged from 20% to 100% of the VSWTs base power. Several observations are in order. First, the time constants are complex and conjugate values for installed capacities greater than 10% of the system demand. Second, it can be seen that the real part of the time constants increases with the wind generation reduction. It is also interesting to point out that 2 tends to sv , that is about 8 s, as the wind generation declines to zero.

Impact of wind turbine virtual inertia on the SFR model
In this case, it is assumed that all VSWTs generate their rating power before the loss of generation incident. The effect of the VSWTs virtual inertia on the accuracy of the estimated maximum frequency deviation is illustrated in Figure 13. In particular, Figure 13a shows the true values of this quantity for virtual inertia gain k h ranged from 1 to 3. The achieved estimation errors considering the motors inertia are depicted in Figure 13b. Interestingly, the accuracy of the proposed model is roughly independent of the VSWTs virtual inertia while this is not the case for the traditional model. In terms of the time t svm , it is also observed that the estimation errors are ranged from −1% to −8% using the suggested model. However, it ranged from −2% to −20% using the traditional one.   Figure 14a shows the true values of the maximum deviation in the COI VSWTs rotor speed for the considered virtual inertia gains. The attained estimation errors to estimate Δ wm deploying the proposed model are presented in Figure 14b. For low VSWTs installed capacities, the errors ranged from −10% to at most −12% which is proportional with the virtual inertia gain. However, the estimation error tends to be independent of the VSWTs virtual inertia gain as the wind farms installed capacity increases.

Impact of SGs inertia on the SFR model
Another quantity which is used in the SFR model is the total inertia of the SGs (excluding the tripped unit). In this regard, the estimation errors achieved by the SFR models considering ±10% errors in the inertia H s are compared with the base case in Figure 15. The base case (the solid traces) is identical with the reported case in the previous sub-section with k h = 1. Figure 15a and b show the estimation errors of the maximum frequency deviation obtained by the traditional and proposed SFR models, respectively. While, the estimation error of the maximum COI VSWTs rotor speed deviation attained by the proposed SFR model is provided in Figure 15c. It can be observed that the estimation errors achieved with 0.9H s and 1.1H s tend to that of the base case when the wind farms installed capacity increases. This is clearly due to the fact that the SGs inertia participation in the system inertia decreases with the wind farms integration level. In addition, it is revealed that the proposed SFR model gives more reliable frequency nadir estimations in comparison to the traditional approach, subjecting to the errors in the SGs inertia.

CONCLUSIONS AND FUTURE WORKS
An SFR model to evaluate the frequency response of the largescale power systems integrating VSWTs was presented. The derived SFR model consists of two differential equations which formulate interaction between VSWTs rotor speed and the system frequency variations following the infeed loss incidents. The system loads inertia was also included in the SFR model. Moreover, a systematic approach was presented to derive the SFR model for the power systems including VSWTs with different capacities as well as various operating points. The developed SFR model was analytically solved for time variations of the system frequency and the VSWTs' rotor speed. It has been revealed that the solutions contain exponential functions with two different time constants, which in turn were dependent upon the wind generation level. Furthermore, the system frequency nadir and maximum reduction in the VSWTs' rotor speed were analytically derived. Finally, the proposed SFR model was verified for estimating the frequency dynamics of a modified IEEE 39-bus system. The results were provided for various VSWTs installed capacities and operating points. It was observed that the proposed SFR model can estimate the frequency nadir with about 3% error regardless of the VSWTs generation and their virtual inertia. In case of the maximum deviation in the VSWTs rotor speed, the estimation error was about 10% for different system parameters. For future work, the proposed SFR model can be used to determine optimal values of the VSGs ensuring the power system stability. On the other hand, the proposed SFR model can be extended to model the frequency response of other VSG-based generations such as photovoltaic systems.

APPENDIX B
For the studied test system, the parameters of the SG power stations and those of the wind farms are listed in Tables B.1 and B.2, respectively. All power quantities are given in percentage of the system base power. The WTs size and distribution on the network are selected according to the wind energy in Europe in 2018 [30]. The WTs models and their nominal wind speed are borrowed from [43]. In Table B.1, the SGs reactive power is denoted by Q s,i . The quantities V st ,i and PF sn,i represent the SGs terminal voltage and their nominal power factor, respectively. In Table B.2, number of the WTs installed in the jth wind farm is given by the parameter N w, j . The quantity P wt , j shows the WTs rated power in MW. In this study, the participation of the wind farms in the total wind installed capacity is specified by introducing the wind farm participation factor, denoted by k wb, j . The summation of this factor for all of the wind farms is: Accordingly, the individual wind farms base power would be given by: where the total wind farms base power is given by: The parameters of the above equation are listed in Table B.3 [44]. Similarly, the WTGs inertia constants can be calculated as follows:

APPENDIX C
The aggregate VSGs mechanical power may be written in a form similar to (3) as follows: where k mpp denotes an optimal power point coefficient for all of the wind farms. On the other hand, assuming the jth wind farm generates P ve0, j before the loss of generation event, the pre-event total wind farms power generation is obtained as: As the system is initially on an equilibrium point, the initial value of P vm equals with P ve0 . Therefore, the coefficient k mpp can be given by: in which, the pre-event COI WTs rotor speed w0 is obtained by substituting the following equation into (14): The flowchart illustrating how this parameter is calculated is presented in Figure C

APPENDIX D
For further validation of the suggested SFR model, a case including two PFR services has been considered as illustrated in Figure D.1. Particularly, the generators SG4 and SG10 provide 40% and 60% of the required frequency response, respectively. The generator SG10 has a 10 s delivery time with no activation delay. While, the generator SG4 has a 5 s delivery time with 1.0 s activation delay [45]. The achieved results for estimation of the COI frequency variations and VSWTs rotor speed deviation with considering the motors inertia are shown in Figure D.2. As observed, even for cases with different PFR services, the proposed SFR model gives better estimations against the traditional model.

APPENDIX E
This appendix provides the detailed simulation results pertinent to Section 4.2. The true results of the wind farms which are achieved through the time-domain simulations of the detailed test system model, are illustrated in Figures E.1 and E.2. In particular, the time evolutions of the VSWTs' rotor speeds are shown in Figure E.1a,b. Additionally, the VSWTs rotor speed deviations with respect to their pre-event values are compared in Figure E.1c. Note that the rotor speeds are represented in percentage of the individual VSWTs base speed. Furthermore, the time evolutions of the VSGs electrical power are illustrated in Figure E.2a,b. These power traces are represented in percentage of the individual VSWTs base power. In addition, the deviation in the aggregated electrical and mechanical powers of the VSGs associated with all the wind farms are given in Figure E.2c, in percentage of the system base power. It can be clearly seen that the VSWTs are running in discrepant operating points during the pre-event condition.