An WAMS-Based Online Generators Coherency Identification Approach for Controlled Islanding

In this paper, an online model-free approach based on the largest Lyapunov exponents (LLE) and angular velocity deviation is provided for coherent groups identification (CGI). Firstly, by using a model-free LLE calculation method, the generator coherency identification (GCI) criterions are proposed by establishing the mathematical relationship between the LLE and angular velocity difference. Secondly, an online CGI scheme for power systems is designed. Compared with the existing measurement-based methods, the breakthrough work of this paper is that the proposed GCI criterions combines the stability theory of nonlinear dynamic trajectory based on LLE with the physical mechanism of the out-of-step between generators, such that it can improve both the speed and the accuracy of the GCI. Besides, the proposed approach needs not to find the optimal observation window size and shorten the required observation window size greatly. Moreover, it can be applied to online CGI with small computation requirement. Extensive test results on New England 39-Bus System and practical East China (EC) power grid, as well as the comparisons with existing CGI methods verify high accuracy and efficiency of the proposed approach.


I. INTRODUCTION
The coherence phenomenon in power system refers that there exhibits the same or similar dynamic behaviors among the generators of the system in a dynamic process after fault [1]. When oscillations of the post-fault power system occur during a dynamic process, synchronous generators are probably to behave in forms of several coherent groups. The timely and effective coherent groups identification (CGI) is critical for model reduction, wide-area control and controlled islanding of interconnected power system [2], [3]. Since the advancement of wide-area measurement system (WAMS) based on phasor measurement units (PMUs) has made it possible to obtain the real-time dynamic information of power systems [4]- [6], many online coherency identification approaches are developed for real-time power system stability analysis and control.
Existing methods for coherency identification are mainly classified into two kinds as model-based methods [7], [8] and measurement-based methods [9], [10]. The model-based The associate editor coordinating the review of this manuscript and approving it for publication was Siqi Bu . methods, such as the slow coherency method in [11], [12] and the methods based on graph-theoretic in [13], [14], are used for CGI by time-domain analysis with the linearized dynamic model of power systems. However, such methods depend on the detailed system topology data which may not always be available, and cannot capture the dynamic features of post-fault power systems, so they are inappropriate for online identification especially in large power grids.
The measurement-based CGI methods, employing realtime dynamic information of the power system from PMUs, are independent of system model such that inaccurate identification brought about by model errors can be avoided. The artificial intelligence (AI) approach is a typical modelfree method and has received great interests. The reported AI methods in the technical literatures are mainly based on support vector cluster [15], artificial neural networks [16], K-harmonic means clustering [17] and Fuzzy C-means [18]. However, these methods require excessive samples obtained off-line and massive training must be carried out, which are hard to realize in practical large-scale power system. References [19] and [20] employ the independent component analysis (ICA) method for coherency identification based on PMU signals. In [21], the Koopman modes is defined for power system and the established nonlinear Koopman modes which has the full information of the system is used for coherency identification. In [22], a coherency identification approach based on multi-flock model by transforming generator data from the observation space to an information space is proposed. In [23], a coherency identification method based on frequency deviation signals from WAMS is presented to dynamically identify coherent generators and electrical areas of an interconnected power system. In [24], ten trajectory dissimilarity indices for the generators' rotor angle and rotor speed trajectories are presented for coherency identification. However, the limitations of these approaches include 1) difficulty to predefine an optimal number of coherent groups, 2) sufficiently long observation window requirement, 3) high computing burden.
Since the generators coherency is to describe the similarity among the post-fault rotor angle trajectories, it can be reflected by the convergency of the rotor angle difference trajectories between the pair of generators. Therefore, the generators coherency problem can be analyzed by using the largest Lyapunov exponents (LLE) stability theory because it involves studying whether the trajectories surrounding the post-fault equilibrium point diverge or converge. The main advantages of using LLE stability theory are that 1) a Lyapunov function is an invariant property for nonlinear trajectories [25], [26], 2) the coherency between two generators can be easily identified just according to the signs of LLE values, specifically a negative (positive) LLE indicates that the two generators are coherent (noncoherent), 3) with model-free LLE estimation approaches, the online modelfree CGI can be realized only using the measured rotor angle data, 4) model-free LLE estimation can be realized by using simple algorithm with small computation, and 5) the LLEbased methods need not to predefine the number of clusters of power systems. Two LLE-based coherency identification approaches are presented in [27] and [28]. They can identify the generators coherency according to the final signs of LLEs by only using the PMU measurement data. However, the two methods require a relatively long observation time due to the LLEs may fluctuate between negative and positive values for quite a long time after disturbance. The optimal window size, which is crucial for the accurate coherency identification, is difficult to be determined when applying the two methods. Besides, the sampling rate and the measurement noise have a big negative impact on the performance of the coherence identification with these methods.
Considering above limitations in existing methods, this paper developed a novel measurement-based approach for identifying coherent groups based on two indexes, LLE and angular velocity deviation. Main contributions of this paper include the following. 1) Combining the stability theory of nonlinear dynamic trajectory based on LLE with the physical mechanism of the out-of-step between generators, the proposed method can realize the coherency identification with low computation cost and high accuracy. 2) The observation  window time required for GCI is shortened greatly to enable timely controlled islanding in smart security and stability control system. 3) The robustness to measurement noise and sampling rate is strengthened due to the fact that very small amount of sample data is required by the proposed approach. 4) The proposed approach is applied to the New England 39-bus system and practical east-China (EC) power grid, and the results are compared with existing model-based and measurement-based methods to verify its performance.
The rest of this paper is organized as follows. The coherency identification criterions and the online CGI scheme are proposed in Section II. Section III presents the results obtained from the application of the proposed approach to the New England 39-bus system and practical EC power grid. In Section IV, some key issues of the proposed technique are discussed. Finally, the conclusion is given in Section V.

II. PROPOSED COHERENT GROUPS IDENTIFICATION APPROACH A. PHYSICAL MECHANISM OF GENERATORS COHERENCY BASED ON POST-FAULT PHASE TRAJECTORY
Since the generators coherency is to describe the similarity between the rotor angle curves of generators in a post-fault power system [7], the GP is taken as the research object. Denote G i -G j (i, j = 1, 2, . . . , n, i = j) is one GP of a n-machine power system, and suppose that G i is the relative lead generator of the GP. The network diagram between G i and G j is shown in Fig. 1. The parameters involved in this section are shown in Table 1.
The motion equations for generator G i and G j are separately described as equation (1) and equation (2) The motion equations of GP G i -G j can be written as where Because equation (3) is similar to the motion equation of single machine infinite bus system [33], the dynamic features of one GP can be analyzed using the dynamic response process of the power system. Therefore, the typical dynamic characteristics of post-fault phase trajectories (δ-ω curves) for both GP and single machine infinite bus system are the same. Typical dynamic characteristics of post-fault δ-ω curves of G i -G j including four different stability conditions are shown in Fig.2 [30], in which the black spots denote the points of δ-ω curves at the fault-clearing time.
• Synchronous condition 1: In this condition, the fault is removed as soon as it occurs. Generator i and j are synchronous because the GP shows stable rotor angle. Corresponding to curve 1 in Fig. 2, δ ij tends to be stable after a period of oscillation, and ω ij increases with a reducing rate and then tends to be zero after a period of oscillation.
• Synchronous condition 2: Generator i and j are synchronous as the GP finally tends to be rotor angle stable. Corresponding to curve 2 in Fig. 2, δ ij oscillates first and then tends to be stable, while ω ij decreases, then oscillates and finally tends to be zero.
• Asynchronous condition 1: Generators i and j are asynchronous as the GP exhibits rotor angle instability. Corresponding to curve 3 in Fig. 2, δ ij keeps increasing, while ω ij keeps increasing quickly after decreasing for a short time.
• Asynchronous condition 2: Generators i and j are asynchronous due to the long duration of fault. Corresponding to curve 4 in Fig. 2, δ ij keeps increasing and ω ij accelerates.

B. COHERENCY IDENTIFICATION CRITERIONS BASED ON LLE AND ANGULAR VELOCITY DEVIATION 1) ESTIMATING LLE FROM TIME SERIES
In this paper, the LLE of the rotor angle of a GP is estimated with the trajectory tracking-based method in [27]. The LLE estimation algorithm of a GP by using rotor angle measurements from PMUs is introduced as follows. i) Obtain the rotor angle samples of the generators of a GP from PMUs. δ i (t) and δ j (t) are separately defined as the rotor angle vectors valued time series data for t = 0, t, 2 t, . . . , N t of the two generators of the GP, where i = j, and t is the sampling period. Generator j is taken as the reference, and the relative rotor angle is defined as: ii) Suppose the number of initial conditions is M , the LLE of δ ij (t) is calculated by equation (5) λ where m = 1, 2, . . . , M , M < k. There is one dimension of embedding space in the LLE calculation algorithm, and the algorithm is an easy approach with small computation cost. Accordingly, it is convenient to realize online LLE calculation especially for generators coherency identification (GCI) of large-scale power systems.

2) ESTABLISHING THE MATHEMATICAL RELATIONSHIP OF THE LLE AND ANGULAR VELOCITY DIFFERENCE
According to motion equation (3), the following equation based on the time series can be obtained as Therefore, equation (5) can be written as Denote algebraic expression | ω ij ((k+m) t)| | ω ij (m t)| as H , then the following calculation equation is obtained.

3) DERIVING THE COHERENCY IDENTIFICATION CRITERIONS
For M does not affect the ultimate LLE-t curve trend of GP [27]- [29], M is considered as a small number in the following analysis. According to the typical dynamic characteristics of the post-fault phase trajectories of GP as shown in Fig. 2, five different swing types of post-fault angular velocity deviation curves can be identified [30]. Fig. 3 shows the angular velocity deviation curves of G i -G j after clearing faults ( ω-t curves) and the corresponding LLE-t curves of the relative rotor angle for different conditions as shown in Fig.2. It should be noted that the lagging generator of the GP is considered as the reference in our analysis, since there is no correlation between the estimating results of the LLE of a GP and the reference generator [27], [29]. Type 1: From the left graph ( ω-t curve) in Fig. 3(a), the angular velocity deviation ω decreases to a value less than 0 rapidly, and then it reaches to the minimum value less than −u. ω finally tends to be 0 after a period of damped oscillations. Correspondingly, there is 0 < H < 1, then there are H > 1 and H < 1, alternately. According to the characteristic of the logarithmic function of equation (8), it can be obtained that the corresponding LLE-t curve starts with a negative value, and finally tends to be a negative value after a period of oscillations, as shown in the right graph (LLE-t curve) of Fig. 3

(a).
Type 2: The ω-t curve in Fig. 3(b) shows that ω soon decreases to a value less than 0, and finally tends to be 0 after a period of damped oscillations. Type 1 and Type 2 are two different kinds of ω-t curve which corresponds to Synchronous condition 1 as shown in Fig. 2. However, the distinction between these two types is that there is always ω > −u for Type 2. According to equation (8) and as shown in the right graph of Fig. 3(b), it can be seen that the corresponding LLE-t curve of Type 2 cannot intersect the horizontal axis compared with Type 1. Therefore, the only difference of the LLE-t curves between the two types is that whether they pass through the zero point.
Type 3: This type corresponds to Synchronous condition 2 as shown in Fig. 2. From ω-t curve in Fig. 3(c), ω first increases with a reducing rate for a short time, then decreases to a value less than 0 and specially it can decrease to a value less than −u. ω finally tends to be 0 after a period of damped oscillations. Correspondingly, H changes from the beginning of H > 1 to 0 < H < 1 quickly, then there are H > 1 and 0 < H < 1 alternately, and there is finally 0 < H < 1. According to equation (8), the corresponding LLE-t curve changes from the beginning possitive value to a negative value soon, then passes through the zero line and . Different angular velocity deviation curves after clearing fault and the corresponding LLE-t curves for the five distinct types. t ZCP is the time when LLE-t curve passes through the zero line from negative to positive for the first time, ω ZCP and δ ZCP are the angular velocity deviation and the relative rotor angle at t ZCP respectively; u (u > 0) represents the value of ω ij at the fault-clearing time; point a and b are the intersect points between the straight line ω = −u and ω-t ; LLE t 0 represents the starting point of LLE-t curve; t m is the time when LLE-t curve monotonically increases to the peak value for the first time; LLE tm and ω t m are the LLE and the angular velocity deviation at t m , respectively.
oscillates, and finally evolves with time maintaining negative values, as shown in the right graph of Fig. 3(c).
Type 4: This type corresponds to Asynchronous condition 1 as shown in Fig. 2, in which the phase trajectory does not pass through the zero line, i.e., ω keeps positive. VOLUME 8, 2020 The ω-t curve in Fig. 3(d) presents that ω decreases for a short time, then keeps increasing. Correspondingly, there is 0 < H < 1 for a short time, then there is always H > 1. According to equation (8), the LLEt curve starts with a negative value, then soon passes through the zero line and finally tends to be positive values, as shown in the right graph of Fig. 3(d).
Type 5: This type corresponds to Asynchronous condition 2 as shown in Fig. 2, in which ω keeps positive. The corresponding H is greater than 1 over the entire time. According to equation (8), the corresponding LLE is always positive, as shown in the right graph of Fig. 3(e).
Based on the above analysis, the following criterions can be used to identify the coherent generators.
• Criterion 1: If the LLE-t curve of the relative rotor angle of a GP passes through zero line, ω ZCP < 0 ( ω ZCP > 0) demonstrates the two generators of the GP are coherent (noncoherent), and t ZCP is the coherency identification time. To further confirm the identification time for Criterion 2, the following auxiliary criterions are proposed, in which LLE t0 represents the starting point of LLE-t curve, and t m is the time when LLE-t curve monotonically increases to the peak value for the first time. LLE tm and ω tm are the LLE and the angular velocity deviation at t m , respectively. Because Criterion 2 is proposed for the condition in which LLE-t curve does not pass through the zero line, the following auxiliary criterions are presented for two stability conditions, Synchronous condition 2 and Asynchronous condition 2 depicted in Section II. A.
• Auxiliary criterion 1: If LLE t0 < 0 and LLE tm < 0, the two generators of the GP are considered to be coherent. (Auxiliary criterion 1 corresponds to Synchronous condition 2).
• Auxiliary criterion 2: If LLE t0 > 0 and ω tm > 0, the two generators of the GP are considered to be noncoherent. (Auxiliary criterion 2 corresponds to Asynchronous condition 2).
In the proposed auxiliary criterions, t m is the coherency identification time.
According to the proposed generators coherency identification crterions, the coherency identification time is the time instant when LLE-t curve passes through the zero line from negative to positive for the first time, or the time instant when LLE-t curve monotonically increases to the peak value for the first time, i.e. the coherency identification time is t ZCP or t m as shown in Fig. 3. It can be observed that the observation window required for generator coherency identification is greatly shortened by using the proposed method.
The flow chart of the proposed GCI approach is shown in Fig. 4. It should be noted that the condition of ii) Select the generator which has the maximum rotor angle at a certain moment (t 1 ) after clearing fault as the reference denoted as G r1 , and then calculate the relative rotor angle between G r1 and the rest of , denoted as δ r 1 i (t) = δ r 1 i − δ i , i = 1, . . . , n, i = r 1 . Meanwhile, calulate the angular velocity deviations of all GPs of the post-fasult power system, denoted as set = { ω ij , i, j = 1, 2, . . . , N , i = j}. It should be noted that the rotor angle would not swing obviously after clearing fault becasue the rotor shaft of generator has a big inertia [31], thus t 1 is the moment of several cycles after clearing fault.
iii) Calculate the LLE sequences for all relative rotor angles by using the algorithm in Section II. B1) to obtain LLE-t curves; then identify the coherent generators of G r1 based on the proposed GCI approach as shown in Fig. 4 to get the first coherent group of the post-fault power system denoted as set r1 ; denote the noncoherent generators of G r1 as set non−r1 . Suppose that the number of r1 is l 1 , and the number of non−r1 is n-l 1 .
iv) Select the generator with the maximum rotor angle of non−r1 as the new reference denoted as G r2 , and calculate the LLE sequences for all relative rotor angles between G r2 and the rest of non−r1 ; identify the coherent generators of G r2 to get the second coherent group of the post-fault power system, denoted as set r2 . Denote the noncoherent generators of G r2 as set non−r2 . Suppose that the number of r2 is l 2 , and the number of non−r2 is n-l 1 -l 2 . 89634 VOLUME 8, 2020 v) Repeat step iv), and the algorithm ends when n-l 1l 2 − . . . − l p = 0. Fig. 5 shows the flow chart of the proposed online CGI scheme. It should be noted that the two machines of which the rotor angel difference is zero are identified as the coherent generators in our proposed algorithm as shown in Fig.5.

III. SIMULATION RESULTS
The proposed CGI approach is tested on New England 39-bus system and EC power grid. Be similar with other measurement-based methods, the proposed method requires generator buses can be fully observed. In the test cases, the sampling rate of the WAMS measurements used for LLE calculation is 120 samples per second, and the number of initial conditions M is chosen as 10 according to a great number of testing results. In order to compare the proposed approach with the traditional LLE-based methods in [27], [28], the observation window size is set as 10 s. Besides, the initial time for calculating is the fault-clearing time. The response results to a disturbance extracted from Power System Department Bonneville Power Administration (PSD-BPA) (a Chinaversion BPA software developed by China Electric Power Research Institute) are used as the real-time measurement data of WAMS. All the tests are performed in MATLAB environment (computer specifications: Intel Core i7-5500U CPU at 2.39 GHz, 8 G RAM).

A. TESTS ON NEW ENGLAND 39-BUS SYSTEM
The New England 39-bus system has 10 generators, 39 buses and 19 load buses. The bus node of generator 39 is the slack bus and the bus nodes of generator 30-38 are PQ bus nodes. The diagram and more detailed parameters related to the test system can be found in reference [32]. A three-phase shortcircuit ground fault occurs in the middle of line 6-7 at 0 s, and  then the fault is cleared at 12 cycles and 16 cycles for Case 1 and Case 2, respectively. Then the system begins to oscillate and finally become unstable for both cases.
Case 1: The response trajectories of rotor angle for all generators after fault are given in Fig. 6. According to the online CGI algorithm shown in Fig. 5, generator 38 is selected as G r1 . The simulation results including LLEt curves and phase trajectories for all GPs, which are formed by generator 38 and the other generators respectively, as shown in Fig. 7. Besides, the involved indexes to identify the coherent generators of generator 38 are presented in Table 2. These results show that only generator 39 is noncoherent with generator 38. The program ends after one step in this case. Hence, two coherent groups are identified in Case 1, which are {30, 31, . . . , 38} and {39}.
Moreover, it can be seen from the time indexes in Table 2 that the coherency identification time for all GPs involved in Case 1 ranges from 0.53 s to 0.74 s, i.e., the time required for the coherent groups identification of the system in Case 1 is 0.58 s after clearing fault. Besides, all values of δ shown in Table 2 are less than 180 • . Therefore, the coherent groups can be identified before the instability of the system in Case 1. The coherency identification result of the proposed approach is the same as the result of the proposed methods in reference [27], [28], which identify the coherency just by observing the final signs of LLE-t curves. However, most of the LLE-t curves as shown in Fig. 7 fluctuate between positive and negative values for quite a long time, and the required coherency identification time is more than 8 s by using the method in [27], [28].
Case 2: To further verify that the generators coherency of post-fault system is related to the fault duration, the fault duration is extended to 16 cycles in Case 2. The response trajectories of rotor angle for all generators after fault are depicted in Fig. 8. According to the proposed algorithm, generator 38 is selected as G r1 . The simulation results including VOLUME 8, 2020   LLE-t curves and phase trajectories for all GPs, which are formed by generator 38 and the other generators respectively, are shown in Fig. 9(a). Besides, the corresponding indexes to identify the coherent generators of generator 38 are presented in Table 3. It can be seen from the involved identification indexes shown in Table 3 that generators 31, 32 and 39 are noncoherent with the generator 38. Hence, the first coherent generator group is denoted as r1 = {30, 33, 34, 35, 36, 37, 38}, and non−r1 = {31, 32, 39}. Then, generator 32 is selected as G r2 . Fig. 9(b) shows the simulation results for all GPs which are formed by generator 32 and the other generators of non−r1 respectively. The involved indexes are presented in Table 3. These coherency identification indexes show that generator 31 is coherent with generator 32, while generator 39 is noncoherent with 32. Thus, the two coherent  Furthermore, it can be observed from the time indexes shown in Table 3 that the time needed for the coherency identification of the system in Case 2 is 0.47 s after clearing fault, and all values of δ shown in Table 3 are less than 180 • s. Therefore, the coherent groups can be identified before the instability of the system occurs in case 2. The coherency identification result of the proposed method is in accordance with the result of the proposed methods in [27], [28] for Case 2. However, the coherency identification time for the system in Case 2 is about 10 s by using the traditional LLEbased methods in [27], [28], which is distinctly longer than the proposed approach.

B. TESTS ON EC POWER GRID
The proposed approach is also tested on the EC power grid. This practical power system is a multi-area interconnected and ultra-high voltage (UHV) power grid of China. The EC power grid is composed by five regional power grids: Shanghai (SH) power grid, Anhui (AH) power grid, Fujian (FJ) power grid, Zhejiang (ZJ) power grid and Jiangsu (JS) power grid. There are 451 generators in EC power grid, and the loads are described by different proportions of constant resistance load model and induction motors. The simplified geographical wiring diagram of EC power grid is shown in Fig. 10. The proposed CGI program is conducted in many different instability scenarios of EC power grid, including power plant instability scenarios and regional grid instability scenarios. For the practical power system is large-scale and complex, we select the following instability scenario as a sample case to present the practical applicability and validity of the proposed approach.
Case 3: The three-phase to ground fault occurs on the tieline connected AH power grid and JS power grid, i.e., the UHV transmission line interconnecting bus HN and NJ at time 0 s, and the fault is cleared after a duration of 8 cycles for Case 3. According to the time-domain simulation results, power plant YZ, PW, PS in AH power grid lose stability in this case because of unbalanced power. There are separately 4, 2 and 2 generators in YZ, PW and PS. The rotor angle curves for all generators after fault are given in Fig. 11. The proposed coherency identification algorithm is performed in Case 3. Specifically: Step 1: Generator YZ-G 3 of power plant YZ in AH power grid is selected as G r1 . According to the proposed CGI algorithm, the relative rotor angle between YZ-G 4 and YZ-G 3 are 0. Therefore, YZ-G 4 is the coherent generator of YZ-G 3 . The LLE-t curves for all GPs, which are formed by YZ-G 3 and the other generators except YZ-G 4 respectively, are shown in Fig. 12(a). By using the proposed GCI criterions for GP, all generators of YZ, PW and PS are identified as the coherent generators of YZ-G 3 . Therefore, the first coherent group of the power grid can be denoted as r1 = {YZ-G 1 , YZ-G 2 , YZ-G 3 , YZ-G 4 , PW-G 1 , PW-G 2 , PS-G 1 , PS-G 2 }, as shown in Fig 10. Step 2: Generator PS-G 1 of PS is selected as G r2 . The relative rotor angle between PS-G 1 and PS-G 2 are 0, so PS-G 2 is the coherent generator of PS-G 1 . The LLE-t curves for all GPs, which are formed by PS-G 1 and the other generators except PS-G 2 , are shown in Fig. 12(b). According to the proposed GCI criterions, all the generators except PS-G 2 are noncoherent with PS-G 1 . Therefore, the second coherent group is denoted as r2 = {PS-G 1 , PS-G 2 }, as shown in Fig. 10.
Step 3: Generator HP-G 1 of power plant HP in AH power grid is selected as G r3 , and LLE-t curves for all GPs in this step are shown in Fig. 12(c). According to the proposed GCI criterions, all generators are coherent with HP-G 1 . Hence, the third coherent group is the set of the rest of generators excepting YZ, PW and PS, denoted as r3 , as shown in Fig. 10. Because the number of non−r3 is 0, the program ends after this step. Therefore, three coherent groups for Case 3 are identified by using the proposed CGI program.
The time durations of coherency identification for GPs involved in the above three steps are summarized in Table 4. The coherency identification time for GPs ranges from 0.27 s to 0.98 s. Therefore, the time needed for the coherency identification of the system in Case 3 is 0.82 s after clearing fault.
Moreover, all of the relative rotor angle values of the GPs when the coherent groups are identified are less than 180 • . Therefore, the coherent groups can be identified before the instability of the system in Case 3. In contrast, it can be inferred from the LLE-t curves in Fig.12 that the traditional LLE-based method has the same result of the proposed method. However, most of the LLE-t curves fluctuate between positive and negative values for quite a long time as shown in Fig.12. As a result, it takes a long time to identify the coherent groups of the system with the methods in [27], [28].

IV. DISCUSSIONS A. REQUIRED OBSERVATION WINDOW SIZE AND COMPUTATIONAL TIME
In LLE-based identification approaches, the time window size needs to be prespecified for online LLE observation. A large window size will lead to reliable results, but hardly applicating in real time, while a too small window size will lead to inaccurate results. Hence, LLE-based methods usually need to find the optimal window size. However, as one prominent advantage of the proposed approach, theoretically the optimal window size needs not to be found, because the LLE and ω are combined together so that the physical mechanism of the out-of-step between generators is utilized for coherency identification. Besides, the needed window size for calculation is small as the coherency between generators can be identified when the LLE-t curve passes through the zero line from negative to positive for the first time or when the LLE-t curve monotonically increases to the peak value for the first time.
To further determine the observation window size, extensive simulation tests have been performed in the New England 39-bus system and EC power grid. Specifically, a three-phase short-circuit ground fault occurred on each un-generator bus of the New England 39-bus system and critical tie-lines of EC power grid, and the fault-clearing time (t c ) were 4 cycles, 8 cycles, 12 cycles and 16 cycles respectively. In the extensive tests, the observation window sizes are set as 1.5 s. There are 217 instability cases in all, and the simulation tests were conducted in these cases by using the proposed approach. The test results are summarized in Table 5.
It can be seen from Table 5 that the coherency identification time needed for the New England 39-bus system and EC power grid ranges from 0.41 s to 1.01 s in all the 217 cases, and the average identification time is 0.775 s. The window size for calculating LLE-t curves in most cases are no more than 1 s. More impressively, the identification accuracy is 100% when the observation window sizes are set 1.5 s. Therefore, the window size of only 1.5 s can achieve accurate CGI with a fast performance by using the proposed approach.
Besides, in the extensive tests, the LLE calculation time for the observation window size of 1.5 s is 0.035 ms, and the total computing time for the GCI is less than 0.07 s on a computer with Intel Core 2.39 GHz i7-5500U CPU and 8G RAM, which means the proposed method runs fast enough to meet the online application requirements.

B. COMPARISONS WITH EXISTING METHODS
In order to further demonstrate the effectiveness of the proposed method, comparisons with the model-based method, i.e., slow-coherency method (SC) [12], and two measurement-based methods, i.e., the traditional LLE-based method (TLLE) [27] as well as the frequency deviation signals based method (FDS) in [23] are presented. Take Case 1 and Case 2 of New England 39-bus system as examples, the comparison results obtained by SC, TLLE and FDS methods are given in Table 6.
When applying the model-based SC method to Case 1 and Case 2, the same initial conditions for both cases are used. It can be seen from Table 6 that the coherent groups results obtained by applying the SC method on the two cases are same. The coherent groups results are fixed under the same initial condition for a power system due to the theory of the SC method. Therefore, the SC method dose not capture the dynamic features of the post-fault power systems such that it cannot be used for different operation conditions and is not suitable for online applications.
In order to compare the identification speed of the proposed method with the two measurement-based methods of TLLE and FDS, the observation window sizes are predefined as 1.5 s and 10 s respectively when applying TLLE and   Table 6 show that the coherent groups are different with the two distinct observation window size of 1.5 s and 10 s by using TLLE and FDS. Comparisons among Table 2, III and VI show that the results obtained by the proposed method are same to that obtained by TLLE and FDS methods with the observation window size of 10 s. In other word, the method TLLE and FDS cannot identify the coherent groups correctly in a short observation time, while the correct results can be obtained by the proposed approach rapidly.

C. SAMPLING RATE
The sampling rate is a critical parameter which should be considered when conducting online CGI by using measurementbased methods. However, for the proposed approach, the sampling rate does not have so much effect on the coherency identification results because only a small number of samples is needed. In on-line application, the sampling rate of a commercially available PMU ranges from 1 sample per two cycles to 2 samples per cycle [35], [36]. Take a noncoherent GP (38-39) and a coherent GP (38-30) in Case 1 as two examples, the coherency identification results of different sampling rates are shown in Fig. 13. It can be seen that the coherency identification results with different sample rates for the two GPs are same. Besides, the differences of the identification time t ZCP for different sample rates are very small. Since the coherency identification for power systems is based on the coherency identification of GP, the accuracy and rapidity of the CGI for systems stays unaffected from different sampling rates using the proposed approach.

D. THE ANTI-NOISE PERFOAMANCE
The measurement noise is inevitable in practice, and it has impact on the performance of measurement-based methods.  For commercially available PMUs, signal-to-noise ratio (SNR) is at least 100, i.e., total vector error is less than 1% [33]. To evaluate its anti-noise performance, the proposed approach is tested under the presence of Gaussian noise with SNR of 40 dB. Fig. 14 shows the testing results for identifying coherency of 38-30 and 38-39 in Case 1. It can be seen that both the accuracy and the rapidity of the GCI are not affected by the measurement noise. It demonstrates that the proposed approach can tolerate more than 2.5% noise which is higher than the noise immunity of the other measurement-based methods in [21]- [25]. In theory, measurement noise has not too much effect on the simulation results because low volume of measurement data from PMUs (usually less than 1 s) is required with the proposed method, and the correctness of this theory is validated by the noise testing results.

V. CONCLUSIONS
This paper presents an online GCI approach based on the LLE and angular velocity difference. The mathematical relationship between the LLE and angular velocity difference is established, and based on the two indexes, the coherency identification criterions are proposed, which are used to develop an online CGI algorithm for power systems. Extensive tests are performed on the New England 39-bus system and practical larger system (EC power grid), and the comparisons with the existing coherency identification methods, including the model-based method and measurementbased methods, are presented. The proposed method features several beneficial properties that motivate its application for power systems.
First, combining the stability theory of nonlinear dynamic system with the physical mechanism of synchronization for generators, the proposed approach can achieve coherency identification before the instability of the system occurring effectively. The observation window size can be shortened greatly compared with other measurement-based methods.
Second, the proposed approach is computationally efficient. In the extensive simulations of New England 39-bus system and EC power grid, the coherency can be determined in milliseconds, indicating that the proposed method runs fast enough to meet the online calculation requirements.
Third, as very small amount of sample data is required to realize online coherency identification, high robustness of this method is assured, minimizing the effect of measurement noise and sampling rate of the practical PMU.
These above advantages demonstrate the potential that our proposed method is an effective tool for rapid and robust GCI in controlled islanding. In the next step, the controlled islanding problem will be considered to avoid blackout for further improvement. Criterion 1 corresponds to Type 1, Type 3 and Type 4 as shown in Fig. 3. Both Type 1 and Type 3 are two synchronous conditions for generators, and Type 4 is one asynchronous condition for generators.

1) PROOF OF SYNCHRONOUS CONDITIONS
Necessity: In the synchronous conditions of Type 1 and Type 3, according to the logarithmic function characteristics of equation (8), the LLE at the time when ω-t curve passes through the zero-line for the first time must be less than zero, and the LLE-t curve soon passes through the zero-line from negative to positive after this time because there is ω < −u during the first oscillation period. The corresponding ω at the time when the LLE becomes positive from a negative for the first time must be in the ω-t curve of interval of [a, b] as shown in Fig. 3. Hence, there must be ω < 0 when LLE-t curve passes through the zero-line for the first time, i.e., ω ZCP < 0. Therefore, the coherent GP implies that ω ZCP < 0. Sufficiency: For the phase trajectories do not pass through the zero line for asynchronous conditions, there must be ω ZCP > 0 for the asynchronous conditions. Therefore, if there is ω ZCP < 0, the GP is coherent.

2) PROOF OF ASYNCHRONOUS CONDITIONS
Necessity: In the asynchronous conditions, the phase trajectories do not pass through the zero line, i.e., there must be ω ZCP > 0 for the asynchronous conditions. Therefore, noncoherent GP implies that ω ZCP > 0.
Sufficiency: According to the above proof for synchronous conditions, the necessary and sufficient condition for the synchronous condition is ω ZCP < 0. Therefore, ω ZCP > 0 implies that the GP is noncoherent.

B. PROOFS OF AUXILIARY CRITERION 1) PROOF OF AUXILIARY CRITERIONS 1
Auxiliary criterion 1 corresponds to Synchronous Condition 2.
Necessity: In Synchronous condition 2, ω soon decreases to a value less than 0 and then experiences damped oscillation. In other time, there is always ω > −u. Thus, the corresponding LLE-t curve starts with a negative value, i.e., LLE t0 < 0, and it does not pass through the zeroline. Therefore, there must be LLE tm < 0 for Synchronous condition 2.
Sufficiency: It corresponds to three kinds of conditions as Synchronous condition 1, Synchronous condition 2 and Asynchronous condition 1 for LLE t0 < 0. Because the LLEt curves in both Synchronous condition 1 and Asynchronous condition 1 must pass through the zero-line, there must be LLE tm > 0 in these two conditions. Besides, LLE t0 < 0 and LLE tm > 0 indicate that the LLE-t curve must pass through the zero-line. Moreover, LLE t0 < 0 and LLE tm > 0 are the sufficient and necessary conditions of Synchronous condition 2 and Asynchronous condition 1. Therefore, LLE t0 < 0 and LLE tm < 0 imply that the condition is related to Synchronous condition 2, and hence the two generators of the GP are coherent.
Necessity: In Asynchronous condition 2, there is always H > 1 for the increasing ω, and hence the corresponding LLE-t curve does not pass through the zero-line, i.e., there must be LLE t0 > 0. Moreover, because the phase trajectory does not pass through the zero line for asynchronous condition, there must be ω tm > 0. Therefore, there must be LLE t0 > 0 and ω tm > 0 in Asynchronous condition 2.
Sufficiency: For LLE t0 > 0 implies ω > 0, there corresponds to two kinds of conditions as Synchronous condition 3 and Asynchronous condition 2. There must be LLE tm > 0 because the LLE-t curve must pass through the zero-line in Synchronous condition 3. According to the characteristics of the logarithmic function in equation (8), ω tm in the ω-t curve is between [a, b], so there is ω tm < 0. Furthermore, ω tm < 0 indicates that the GP is synchronous, thus LLE tm > 0 and ω tm < 0 imply that the synchronous condition belongs to Synchronous condition 3. Therefore, if there are LLE tm > 0 and ω tm > 0, the two generators of the GP are noncoherent, and it belongs to Asynchronous condition 2.