Frequency Stability Enhancement of Low-Inertia Large-Scale Power System Based on Grey Wolf Optimization

The high penetration of converter-based distributed generations (DGs) to power system can give rise to the lack of rotational inertia while replacing the conventional synchronous generators (SGs), which provide the primary frequency reserve (PFR) in power systems. As the result, the decrease in PFR aggravates the frequency stability. To overcome this problem, the droop coefficients of governors in the remaining conventional SGs must be re-determined newly and properly. This paper proposes a new solution based on the grey wolf optimization (GWO) method to optimally select the droop coefficients of SG governors in the low-inertia large-scale power system due to the high penetration of renewables. The proposed solution is very effective for reducing the computational effort significantly, and it is able to recover not only the steady-state but also the transient frequency stability. To verify the effectiveness of proposed optimization solution based on the GWO method, several case studies are carried out on the practical Korea electric power system with the penetration of wind power plants of 4 GW.


A. MOTIVATION AND INCITEMENT
Globally, many developed countries are increasing the hosting capacity of distributed generations (DGs) based on renewable energies such as photovoltaic (PV) and wind. This is because they give many benefits in terms of environment, energy efficiency, and sustainability, etc. Nevertheless, many studies have reported that this great energy transition can cause the serious problem of frequency stability. For example, the transient frequency stability after a sudden loss of generation has been secured by the inertial response and primary frequency reserve (PFR) provided from synchronous generators (SGs) in the power system so far. If these SGs are replaced with many converter-based DGs, and they always provide their maximum powers without inertia controls, it will cause a serious problem for the transient frequency stability in this low-inertia power system [1]- [5].
When a large disturbance occurs in a power system, several indices are used to evaluate the frequency stability. In particular, the so-called frequency nadir (FN), which is the lowest frequency point during a transient condition, is considered for the transient frequency stability analysis. On the other hand, the settling frequency (SF), which is the value of frequency stabilized in the steady-state after a disturbance occurs, is used to consider the steady-state frequency stability. To keep the system frequency stable by the PFR, the FN needs to be maintained above its required value before the under-frequency load shedding relay starts to operate. Also, the SF must be higher than the maximum steady-state frequency deviation, for which the automatic generation control (AGC) operates to restore the system frequency to its nominal value. The relation between the frequency response in the steady-state and the amount of PFR can be defined with the composite frequency response characteristic, which is so-called β [6]. Thus, it is typically measured in MW/0.1 Hz, and it characterizes the relationship between the active power outputs from conventional SGs and SF after a disturbance occurs.
When the DGs are replaced with the conventional SGs participating in the PFR, the value of β is reduced while causing the decrease in both FN and SF. To figure this out, β must be recovered to its original value by re-determining the droop coefficients of governors in the remaining conventional SGs within their operating margin. In contrast, the FN during a transient state is strongly affected by the ramp rate of the governor as well as β. Note that the ramp rates of governors are also different depending on the type of turbine-governor. Thereby, even when the value of β is remained constant by the same total amount of PFR, the dynamic frequency response of the power system can become totally different according to how the droop coefficients of governors are determined.

B. LITERATURE REVIEW
To solve the frequency stability problem caused by DGs penetration, several studies have been carried out by using the linear and convex optimization methods [7], [8]. However, a large-scale power system has a highly nonlinear behavior, and its optimization requires large computational efforts. To deal with this issue, many heuristic optimization algorithms such as genetic algorithm [9] and particle swarm optimization (PSO) method [10]- [12] have been studied. For example, the PSO method was applied to establish an optimal schedule for renewable energies [13]. Also, the fusion firefly method was applied for PV to improve the maximum power point tracking (MPPT) performance [14]. However, these do not study the optimization in a large-scale power system where significant computational effort is required. After the grey wolf optimization (GWO) method [15] was recently developed, it has also been applied to various optimization problems [16]- [18] of the field of power engineering. In particular, this method has the advantage of reducing the computational effort. This is because after modeling this method for the power system, it does not require additional parameter tuning. In addition, this method has a great advantage on local optimum avoidance. Therefore, the GWO method can provide an effective solution on a large-scale power system where high computational effort is required.
Moreover, various control methods have been developed for DGs to provide additional frequency responses [19], [20]. However, due to the variability of the renewable energy sources, DGs may not provide the appropriate PFR whenever it is required. Also, even though the penetration level of DGs is increasing, most of the PFR is still provided by SGs. Therefore, the operation of the SGs must be considered.

C. CONTRIBUTION AND PAPER ORGANIZATION
This paper makes the new contribution by applying the heuristic optimization algorithm to optimize the droop coefficients of governors in remaining conventional SGs when the inertia of a large-scale power system becomes low due to the high penetration of wind power plants (WPPs). As the result, it can effectively recover both steady-state and transient frequency stability of a large-scale power system while keeping β constant. Also, this paper focuses on the effective and efficient co-simulation framework in the electromagnetic transient program (EMTP) level for the use of DIgSILENT Pow-erFactory ® and MATLAB ® software to solve the nonlinear optimization problem on a large-scale power system. Moreover, the GWO method is used as a heuristic optimization algorithm considering its advantages, such as local optimum avoidance and low parameter tuning requirements [18]. The major contributions of this paper are summarized as follows:  The relationships between factors of SGs and the transient and steady-state frequency stabilities are determined based on mathematical analysis.  An optimization problem for improving the frequency response of the SGs through heuristic optimization algorithms is formulated with the theoretical analysis.  The high computational effort required for optimization in a large-scale power system is dramatically reduced with the proposed DIgSILENT PowerFactory ® and MATLAB ® based co-simulation framework applied with the GWO method.  The proposed method re-determines the droop coefficients of SG governors to provide optimal frequency response for a large-scale power system, improving the system frequency stability.
This paper is organized as follows. Section II describes the frequency stability of the power system with the mathematical analysis. Then, Section III formulates the optimization problem for the application of the GWO method to a largescale practical Korea electric power system. In Section Ⅳ, the PSO and GWO methods are briefly introduced, and thereafter how to implement the GWO method is explained. Section V verifies the effectiveness of the proposed method with several case studies under the DIgSILENT PowerFactory ® and MATLAB ® co-simulation environment. Finally, the conclusion and future work are given in Section VI.

A. ANALYSIS OF FREQUENCY RESPONSE
When a disturbance occurs in a power system, the dynamic frequency response can be analyzed by the swing equation [21] given as where Hsys is the system inertia constant in s, and f(t) is the system frequency in pu. Pm(t) and Pe(t) are the system mechanical and electrical powers in MW, respectively. Also, where f(t0) is the system frequency before a disturbance oc-curs. When a disturbance occurs in the power system, the frequency of each bus has a slightly different dynamic characteristic. However, since there are thousands of buses in a large-scale power system, the frequency of every bus cannot be evaluated for frequency stability evaluation. Therefore, system frequency based on the center of inertia (CoI) is used to assess the frequency stability of the power system [22].
When the frequency deviation exceeds the dead-band (DB) of the governor, the frequency stability is secured by the PFR. As mentioned previous, the frequency response is governed by β, which is calculated as 0 PF 100 1 10 where Si, PFi, and Ri are the power capacity, power factor, and droop coefficient of governor in i-th SG, respectively. Also, f0 is the nominal value of frequency. In particular, as the PFR is provided from SGs, the dynamic frequency response is governed by the ramp rate of governor measured in MW/s. After a large power generator is suddenly disconnected, f(t) in (3) can be modified as where Ploss is the size of generation loss, and Crr-sys is the system ramp rate for all governors in the power system. When the frequency response is beyond the DB, the governors of SGs start to activate, as shown in Fig. 1. Then, the value of the FN can be obtained by considering the time period from t0 to tNAD as where tDB is the time when the frequency response exceeds the DB of the governor, and tNAD is the time when the FN occurs. Also, KEsys is the total system kinetic energy provided from conventional SGs [23]. The first and second integral terms in (6) are represented as A1 and A2, which are yellow and purple areas in Fig. 1, respectively. They become negative because Pe(t) is higher than Pm(t) during the time period from t0 to tNAD. Therefore, the value of f 2 (t0), which is close to 1, will be canceled out to the small deviation value. It is known that the value of f(tNAD) in (6) depends on Hsys, Ssys, KEsys, and Crr-sys. For example, if Hsys and Ssys are large (therefore, KEsys is also large), the FN increases. Also, when Crr-sys is high, the slope of Pm(t) becomes steeper. This causes to decrease A2. As the result, the FN is also increased.
Also, the settling frequency deviation ∆SF at the steadystate after a disturbance shown in Fig. 1 can be calculated by considering the DB of the governor as

B. EFFECT OF HIGH PENETRATION OF RENEWABLES
When the power system has a high penetration of renewables by replacing the conventional SGs with DGs, the frequency stability is directly affected. Firstly, KEsys is calculated as where Hi is the inertia constant of i-th SG. As KEsys is decreased due to the high penetration of DGs, two integral terms in (6) are increased. This causes to make the FN lower. Secondly, when the PFR is provided by only the conventional SGs, replacing the conventional SGs with DGs causes to reduce both β and Crr-sys. Therefore, the FN and SF are decreased by (6) and (7), respectively. In summary, the SF is affected by β depending on the values of Si, Ri, and PFR, as shown in Fig. 2

III. PROBLEM FORMULATION TO APPLY GWO METHOD
As mentioned before, to improve the steady-state frequency stability, β is required to recover by re-determining Ri of governors in the remaining conventional SGs. However, this does not guarantee the enhancement of the transient frequency stability. This is because the value of KEsys is still kept low due to the high penetration of renewables, and Crrsys can be varied depending on how Ri of governors is re-determined. Therefore, this section formulates the optimization problem based on the GWO method, which is recently used in many applications, to recover β while maximizing Crr-sys.

A. ASSUMPTIONS AND CONSTRAINTS
Although Crr-sys is not a constant value in practice [24], it is assumed to be fixed as the slope of Pm(t) from t0 to tNAD in Fig.  1 during the optimization process. Also, Ri of hydro and pump turbine-governors are maintained to their initial values while re-determining Ri of gas and steam turbine-governors. This is because the ramp rate for SGs with hydro and pump turbinegovernors is much lower than that with gas and steam turbinegovernors [25]. Thereby, these turbine-governors have a relatively small effect on Crr-sys.
The built-in models of IEEEG1, IEESGO, TGOV1, GAST, GAST2A, GGOV1, HYGOV, and PIDGOV on the DIgSILENT PowerFactory ® software [26] are used as gas, steam, and hydro and pump turbine-governors in simulation. The practical ranges of Ri (which are defined in the operation regulation of the Korea electric power system) are given in Table I, and they are carefully considered as the constraints of the optimization problem.

B. DEFINITION OF OBJECTIVE FUNCTION
During the optimization process, the objective function, J aims to maximize Crr-sys as where the vector, = [R1, R2, …, Rn] t is formed with the droop coefficients of governors, and n is the total number of the remaining conventional SGs. Also, − indicates the system ramp rate for all governors at k-th iteration during the optimization process, and it is calculated as where is the re-determined droop coefficient of governor in the i-th conventional SG at k-th iteration. Also, Pm,i( ,tNAD) and Pm,i( ,t0) represent the mechanical power of i-th conventional SG with the droop coefficient of at tNAD and t0, respectively. In order to recover the SF effectively, the β at k-th iteration is required to be kept at its original value before the penetration of DGs.
One of the most common ways to adapt the equality constraint is to use the penalty factor [27]. Then, the modified objective function, Jmod for this constrained optimization problem, can be defined as where β k and β 0 are β at k-th iteration and before the penetration of DG, respectively. And, p(β k ) is the penalty factor based on β k . Also, θ[q(β k )] in (12) is defined as As the value of q(β k ) in (13) is increased, a larger penalty factor is applied to Jmod in (11). As the result, Jmod is maximized while finding the optimal droop coefficients of governors in the remaining conventional SGs. This is able to improve both steady-state and transient frequency stability while particularly providing strong support for the FN.

IV. IMPLEMENTATION OF GWO METHOD
When the heuristic optimization algorithm is applied to solve the nonlinear optimization problem on a large-scale power system, it generally requires a large computational effort. For example, to run the Korea electric power system, which has 400 SGs and 1900 buses, for a simulation time of 25 s by the DIgSILENT PowerFactory ® software only, it takes 150 s under the computational environment with Intel Core I7 CPU at 3.2 GHz and RAM of 8 GB. Moreover, population-based heuristic optimization algorithms like the PSO and GWO methods are required to repeatedly run this single-shot simulation for all its particles in every iteration. This results in excessive computational effort to solve the nonlinear optimization problem.
This study significantly reduces the large computational VOLUME XX, 2021 5 effort by proposing the new co-simulation framework with the effective interaction of both DIgSILENT PowerFactory ® and MATLAB ® software while still obtaining the global optimization solution. That is, the information of turbine-governor models, their input, and the initial amount of PFR of all SGs are firstly obtained in the DIgSILENT PowerFactory ® based simulation before the optimization process. Then, only the droop coefficients of governors are updated during the optimization process by the heuristic optimization algorithm with the fixed (initial) input of turbine-governor models in MATLAB ® based simulation. Finally, these optimized values are used in the DIgSILENT PowerFactory ® based simulation.
Various population-based heuristic optimization algorithms require additional parameter tuning to obtain a better solution for the optimization problem. Unfortunately, as high computational effort is required for optimization in a largescale power system, additional parameter tuning can increase even more computational effort. However, unlike other population-based heuristic optimizations, the GWO method does not require parameter tuning after designing the algorithm [16]. Also, this method has a good performance on local optimum avoidance. Thus, the GWO method is adopted for optimization in a large-scale power system.

A. OVERVIEW OF PARTICLE SWARM OPTIMIZATION
To verify the effectiveness of the proposed framework with the GWO method, its performance is compared with that of the PSO method, which is one of the representative heuristic optimization algorithms. In general, the PSO method provides a good performance in solving nonlinear optimization problems. The particles used in the PSO method iteratively change their position in the search space to optimize the fitness value of the given objective function. The optimization process begins with a randomly selected position on each particle. At each iteration, the particles move stochastically with the velocity based on three factors, which are the current particle position, best position experienced on each particle itself, and best position experienced globally by the entire swarm. Then, the velocity and position are updated in every iteration as x v (16) where the vectors, (∈ R n×1 ) and (∈ R n×1 ) are the position and velocity of j-th particle at k-th iteration, respectively. Also, the vectors, , (∈ R n×1 ) and (∈ R n×1 ) are the best position experienced by j-th particle and the global best position experienced by the swarm at the k-th iteration, respectively. ω is the inertia weight. The rand1 and rand2 are the random values of uniform distribution from 0 to 1. The coefficients, b1 and b2 are the cognitive and social parameters for controlling the velocity, respectively. They are both set to 2. However, in order to obtain the optimal solution, b1 and b2 have to be additionally tuned. In other words, additional computational effort is required for parameter tuning. More details about the PSO method are given in [10]- [12].

B. OVERVIEW OF GREY WOLF OPTIMIZATION
The GWO method adopts population-based heuristic optimization. It imitates the group hunting mechanism of grey wolves in nature. The group hierarchy of grey wolves is composed of alpha, delta, gamma, and omega wolves. The alpha wolf is the first level of the group, which is the leader of wolves. It has the role of decision-making during hunting. The delta wolf is the second level of the group, which both helps the alpha wolf and dominates the lower-level wolves. The gamma wolf is the third level of the group, which dominates the omega wolf, which is the last level of the group. To mathematically model the group of grey wolves with the hunting mechanism, the distance of prey and the wolf is regarded as the fitness value of the objective function. Then, the solution with the most optimized fitness value is considered as the alpha wolf. Consequently, the delta and gamma wolves are considered as the second and third best solutions. In summary, the hunting mechanism of grey wolves are modeled as 1 kk mp where the vector, (∈ R n×1 ) is the distance from prey to grey wolf. And, the vectors, and are the prey and m-th grey wolf position vectors at k-th iteration, respectively. The operator, is the element-wise multiplication of vectors. Then, the vectors, and in (17) and (18) (20) where h is linearly decreased from 2 to 0 during the optimization process. And, the vectors, 1 and 2 have random numbers in the range from 0 to 1. Because the alpha, delta, and gamma wolves are assumed to be the wolves closest to the prey in each iteration, the remaining wolves are led by these wolves. Then, +1 in (18) can be defined as 12 3 ,, This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Step 1 Step 2 Step 3 Step 4 Step 5

FIGURE 3. Overall procedure to implement the GWO method by using DIgSILENT PowerFactory ® and MATLAB ® based co-simulation framework.
Note that unlike other heuristic optimization algorithms, including the PSO method, the GWO method does not require parameter tuning after being modeled. Therefore, it can significantly reduce the computational effort when optimization is carried out on a large-scale power system.

C. IMPLEMENTATION OF OPTIMIZATION IN CO-SIMU-LATION FRAMEWORK
The overall procedure to implement the GWO method by using DIgSILENT PowerFactory ® and MATLAB ® based cosimulation framework is shown in Fig. 3. Firstly, after replacing the conventional SGs with DGs, the information on the system frequency, decreased β, and all parameters of remaining turbine-governors (including droop coefficients) in DIgSILENT PowerFactory ® based simulation are transmitted to MATLAB ® based simulation environment, where the optimization by the GWO method is implemented as follows.
Step 1) Before the optimization process starts, all parameters are initialized to apply the GWO method. Note that each wolf represents the group of turbine-governors in the remaining conventional SGs, and its position indicates the droop coefficients. They are randomly selected within the range given in Table I depending on the type of turbine-governor. Also, the population size and the maximum number of iterations, kmax, are set to 50 and 500, respectively.
Step 2) The wolf positions are applied to the turbine-governors modeled in MATLAB ® based simulation. Then, they generate their frequency output responses with the input of frequency deviation information. The computational time required to generate this response is about 6 s. Thereby, it significantly reduces the computational effort compared to running the entire power system to obtain the frequency output response in DIgSILENT PowerFactory ® based simulation. Finally, it calculates the fitness value, which is − of wolves, by (10).
Step 3) The wolves that violate the equality constraint given in (13) are sorted out. By this process, the fitness value of the wolves whose β k exceeds β 0 is canceled out by the penalty factor calculated by (12). Note that β k of each wolf is calculated based on their updated position (re-determined droop coefficients) by (4). Therefore, wolves with high fitness value due to the exceeding β k instead of the improved − are sorted out.
Step 4) By comparing − of each wolf, it determines the alpha, delta, and gamma wolves. Therefore, these wolves are the ones with improved − while having the same value of β k as β 0 . VOLUME XX, 2021 7 Step 5) The optimization process from Step 2) to Step 4) is repeated until the number of iterations reaches to kmax.
After the optimization process, the final position of the alpha wolf represents the optimal droop coefficients of governors in the remaining conventional SGs, while maximizing Crr-sys. It is important to note that the amount of PFR is still kept to β 0 . Then, they are applied to the governors in DIg-SILENT PowerFactory ® based simulation, where the corresponding dynamic frequency response and maximum support for raising the FN are evaluated. The optimized droop coefficients might not recover the FN up to the limit of under-frequency load shedding relay. This means that the amount of PFR needs to be increased. Thereby, if it happens, the above procedure is required to implement with the increased β 0 in Step 3).
As mentioned previously, when the personal computer with Intel Core I7 CPU at 3.2 GHz and RAM of 8 GB is used, the computational efforts by the proposed co-simulation framework and the only DIgSILENT PowerFactory ® software are compared in Table II. It is observed that it takes 6 s and 150 s, respectively, for the simulation run of 25 s to obtain the dynamic frequency response. When the GWO method is applied with 50 wolves and 500 maximum number of iterations, a total of 25,000 simulation runs of 25 s are required to find the optimized solution. As the result, it takes 42 h and 1042 h (estimated) for applying two methods, respectively. This clearly proves that the computational effort is dramatically reduced by the proposed co-simulation framework. Note that the weekly and monthly supply and demand planning for Korea electric power system is carried out every Friday and every 25 th of every month, respectively. Therefore, it is possible to provide the solution on optimal droop coefficients of SG governors for the long-term planning within time.

A. CHARACTERISTICS OF KOREA ELECTRIC POWER SYSTEM
The practical Korea electric power system has about 400 SGs with a total power generation capacity of 145 GW. Also, 350 SGs among them are equipped with the governor. The peak load demand and corresponding power generation during winter in 2020 are about 82.4 GW and 83.9 GW, respectively, and Hsys is 4.46 s. For the frequency stability maintenance, under-frequency load shedding relay is operated if the system frequency falls beyond 1 Hz, and the maximum steady-state frequency deviation band is 0.2 Hz. According to provinces, the Korea electric power system is divided into six areas, as shown in Fig. 4. The amount of power generation from different types of SGs and the load demand in each area are given in Table III. The load demand of Area 1 including Seoul, which is the capital of South Korea, takes a large portion of the total load, and it is much higher than the amount of power generation from Area 1. Therefore, a large amount of power must be transmitted from the other areas via the high-voltage transmission lines. In particular, the composition of several types of SGs (including types of turbine-governors) is also different in each area. The nuclear power plants, which are mostly located in Areas 5 and 6, play the role of the base power generation in South Korea. In addition, the coal power plants as load-followers are concentrated in Area 4.  Note that many power plants have a turbine-governor to provide the PFR. In particular, the coal power plants have a steam turbine-governor, and the combined cycle power plants are composed of steam and gas turbine-governors. Moreover, the nuclear power plants use a large steam turbine-governor. However, they do not participate in the PFR. As the result, the Korea electric power system has unique regional characteristics depending on several factors including the composition of different types of SGs and the amount of load demand, etc.

B. EFFECT OF HIGH PENETRATION OF WIND POWER PLANTS
For the base case study, the existing 6 WPPs shown in Fig.  4, which are located in areas of Muchang, Pyeongchang, Shinan, Taebaek, Yeonggwang, and Yeongyang, are considered. Their total power capacity is 529 MW. The detailed modeling of the permanent magnet synchronous generator (PMSG) is used for the wind turbine. Also, assume that it operates in the MPPT mode [28]. Note that the government of South Korea is planning to increase its capacity up to 17.7 GW by 2030.
To analyze the effect of high penetration of WPPs on frequency stability of the system, the planned WPPs of 4 GW are additionally applied to 4 areas, which are Seonamhae, Shinan, Yeongdeok, and Yeonggwang, as shown in Fig. 4. Then, the ShinKori nuclear power plant of 1400 MW in Area 6, which is one of the biggest power plants in South Korea, is suddenly disconnected at 2 s. The result is shown in Fig.  5. When the additional WPPs are installed, KEsys is decreased from 465.8 to 452.6 GVAs. Also, β is decreased from 1940 to 1820 MW/0.1 Hz, and Crr-sys is decreased from 210.48 to 186.08 MW/s. As the result, it is clearly observed that the SF is decreased from 59.908 to 59.904 Hz, and the FN is also decreased from 59.866 to 59.857 Hz.

C. EFFECT OF OPTIMIZED SYSTEM RAMP RATE FOR ALL GOVERNORS
For the base case study in Fig. 5, 199 SGs among a total of 236 operating conventional SGs are equipped with a governor. Thereafter, because some conventional SGs are replaced with the WPPs of 4 GW, the droop coefficients of governors in other 124 SGs are re-determined to recover the reduced β of 120 MW/0.1 Hz using the proposed co-simulation framework shown in Fig. 3. As mentioned previously, the SGs of nuclear power plants and the SGs with the hydro and pump turbine-governors are not considered for optimization. After defining the objective function as described in Section III-B, the required information of the system frequency, decreased β, and all parameters of the remaining turbine-governors (including droop coefficients) are transmitted to MATLAB ® based simulation environment. Thereafter, the solution is obtained using the PSO and GWO methods, respectively. The variations of the average droop coefficient in each area during the optimization process using the PSO and GWO methods are shown in Fig. 6. Also, the fitness values of the objective function (which is Crr-sys) are shown in Fig.  7. It is observed that the average droop coefficients are converged by the PSO and GWO methods after 438 and 398 iterations, respectively. Also, when the GWO method is applied, the converged value of Crr-sys is 5.43 MW/s higher than that of the PSO method. This means that it can successfully find the more optimal solution than the PSO method. Note that the solution from the PSO method can be improved by tuning the parameters of b1 and b2 in (15). However, this would significantly increase the computational effort since it has to run the whole power system iteratively. Table IV summarizes the optimized values of the average droop coefficient in each area and the corresponding Crr-sys by the two methods. Therefore, it is clearly shown that Crr-sys may differ depending on how the droop coefficients of governors are determined under the same value of β.

D. CASE STUDIES FOR ENHANCING THE FREQUENCY STABILITY
To verify the performance of the GWO method to improve the frequency stability, three additional case studies with recovering β are carried out on the practical Korea electric power system as given in Table V. Then, their results are compared by those of the base case and the case with the penetration of WPPs of 4 GW without recovering β. The results of system frequency and system mechanical power when the biggest power plant in South Korea is suddenly disconnected at 2 s are shown in Fig. 8. It is observed from the response of system frequency in Fig. 8(a) that the decreased SF is recovered for all cases when compared to the case with the high penetration of WPPs without recovering β. This is because β is recovered from 1820 to 1940 MW/0.1 Hz for all cases. However, the value of Crr-sys varies in each case, as given in Table VI, and therefore the degree of improvement on the FN is different. For example, for cases 1 and 2, the FN is not recovered to that of the base case. This is because the improvement on Crr-sys is not enough to compensate for the reduced KEsys due to the high penetration of WPPs, while recovering β. Therefore, the FN in cases 1 and 2 is lower than that of the base case. On the other hand, the FN in case 3 is recovered to that of the base case, and the FN in case 4 becomes higher than that of the base case. This means that the improvement of Crr-sys in case 4 is most successfully achieved. The values of power capacity of the WPPs, KEsys, Crr-sys, β, FN, and SF for all cases are compared in Table VI.  Correspondingly, the state trajectories between system mechanical power and frequency for all cases are shown in Fig. 8(b). As expected, all SF points of the state trajectories for cases 1 to 4 (including the base case) are the same because β is equally set to 1940 MW/0.1 Hz for these cases. Furthermore, the spiral curves for cases 1 and 2 (see the green and cyan dashed curves) are deeper than that of the base case (see the black solid curve). This means that the FN is arrested lower as shown in Table VI, and the peak value of system mechanical power is higher. As the result, the corresponding oscillation of the mechanical power becomes high. In contrast, the spiral curves of cases 3 and 4 are shallower than those of cases 1 and 2. In particular, it is clearly shown that the spiral curve of case 4 (see the red solid curve) is the shallowest among all cases. This means that the peak value of system mechanical power is the smallest, and the FN is arrested at the highest among all cases.
Additionally, the performance of the proposed method is evaluated when the output power from the DGs is suddenly decreased by 1800 MW at 2 s due to the variability of the renewable energy sources. The results of system frequency and system mechanical power are shown in Fig. 9. It is observed that the SF is decreased by 0.034 Hz in all cases compared to the previous scenario where the power plant was suddenly disconnected. This is because the value of Ploss is increased by 400 MW while decreasing the SF by (7). Also, the FN is arrested at 59.804 Hz, 59.795 Hz, 59.804 Hz, 59.805 Hz, 59.809 Hz, and 59.819 Hz, respectively, for the base case, the case with high penetration of WPPs, case 1, case 2, case 3, and case 4. In addition, it is clearly shown that the spiral curve of case 4 (see the red solid curve in Fig. 9(b)) is the shallowest among all cases as in the previous scenario.
In summary, all results above verify that the proposed optimization solution by the GWO method is able to effectively improve the frequency stability of a large-scale power system with the high penetration of WPPs.

VI. CONCLUSION AND FUTURE WORK
This paper proposed the new solution to optimally select the droop coefficients of governors of many conventional generators in the low-inertia large-scale power system due to the high penetration of renewables by using the grey wolf optimization (GWO) method. Firstly, the impact on frequency stability due to the penetration of converter-based distributed generation (DG) was carefully analyzed. Then, after the problem formulation to apply the GWO method, its implementation on a practical large-scale power system was accomplished in DIgSILENT PowerFactory ® and MATLAB ® based co-simulation framework.
The results from several studies verified that the proposed optimization solution based on the GWO method improves not only the steady-state frequency stability, but also the transient frequency stability most effectively and efficiently. Moreover, the proposed solution significantly reduced the computational effort required for optimizing a large-scale power system. Therefore, the proposed method would be expected to be preferably applied in a large-scale power system with a high penetration of renewables to keep its frequency stability.
Furthermore, the renewable energy-based DGs are required to provide various frequency responses by using improved control methods instead of the conventional maximum power point tracking control method. Therefore, the system operator must consider the integrated operations of both new DGs and conventional synchronous generators (SGs) to keep the power system stability.