{'fraction_non_alphanumeric': 0.028671575082064502, 'fraction_numerical': 0.011090298770423397, 'mean_word_length': 4.487685731732139, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Power System Transient Stability Assessment Using Couple Machines Method

Analyzing the stability of the power system by using a few machines is promising for transient stability assessment. A hybrid direct-time-domain method that is fully based on the thinking of partial energy function is proposed in this paper. During post-fault period, a pair of machines with high rotor speed difference is defined as couple machines, and the stability analysis of the system is transformed into that of several pairs of couple machines. Based on the prediction of power-angle curve of couple machines within a sampling window after fault clearing, the proposed method avoids the definition of Center of Inertia (COI) and it can also evaluate the stability margin of the system by using the predicted power-angle curve. Simulation results demonstrate its effectiveness in transient stability assessment.


I. INTRODUCTION
Power systems nowadays operate in a complicated state and may be subjected to stressed operating conditions, which demands great accountability of the TSA. Since computational complexity of time-domain simulation becomes intolerable with the increase S of system scale [1], transient energy methods are always attractive to system operators for the online transient stability assessment. In regards to the analysis based on transient energy function, the RUEP method and sustained fault method have received considerable attention and also achieved early advance [2], [3]. Based on these works, Fouad and Stanton used a series of simulations of a practical power system to illustrate the instability phenomenon of the power system, and their work merited two conjectures [4]: (i) Not all the excess kinetic energy at fault clearing contributes directly to the separation of the critical machines form the rest of the system, part of kinetic energy should be corrected.(ii) If more than one machine tends to lose synchronism, the instability is determined by the motion of some unstable critical machines.
Both conjectures marked a milestone to the transient stability studies because they originated two different direct methods. Stimulated by (i), Xue and Pavella developed the EEAC and SIME that were based on kinetic energy correction and group separation [5]- [7]. EEAC and SIME are both proved to be highly efficient, and they had achieved great success in academic research as well as industrial applications. In the meaning time, enlightened by (ii), some work attempted to observe the stability of the system from a non-global perspective. In [8], [9], Vittal and Fouad stated that the instability of a single machine would be identical to the instability of the system and IMEF was firstly proposed therein. Later, Stanton performed a detailed machine by machine analysis of the energy of a multi-machine instability [10]. Stanton also defined PEF to quantify the energy transactions that make a machine transfer from stable to unstable. Then in [11], [12], PEF is utilized to quantify the energy of a local transient control action which is in contrast with the global view. Later, Rastgoufard et al. pointed out that the instability of the system is decided by a few unstable critical machines rather than noncritical machines [13]. The authors also utilized EAC of critical machines to determine the transient stability of a multi-machine system. Haque proposed a strategy to compute the CCT using PEF of critical machines in [14] whereas Ando presented a potential energy ridge which can be used to predict the single machine stability [15]. Recently, Lu and Yu initiated a method using a pair of machines to evaluate the stability index of the system [16], yet the foundation as well as the mechanism of this method has not been fully investigated. Although non-global methods are proved to be less conservative in stability analysis [9], these distinctive methods are at a standstill in recent years, making them a minority and fall into disuse. If one explores deep into these methods, the key problem is that the computation of the critical energy for these non-global methods is rather complicated. To be specific, in IMEF method, the evaluation of critical energy of the critical machines is fully based on the sustained fault trajectory, which means that the computation complexity of IMEF is the same as that of the sustained fault method. In PEF method, the virtual linear trajectory is applied to compute the deceleration area from fault clearing point to the DLP, however its computation is based on the time-consuming simulation of the actual simulated trajectory [11], [12]. In this paper a hybrid direct-time-domain method is proposed. This method can be seen as a continuity of the PEF method. We first clarify that the motion of a critical machine in COI reference in IMEF method and PEF method is identical to that of a two-machine subsystem consisting of a critical machine and a virtual COI machine. Then we replace the virtual COI machine in PEF with a real machine in the system, and define couple machines as the pair of real machines with high rotor-speed difference after fault clearing. Similar to the PEF method that only focuses on the stability state of the critical machine in COI reference, the proposed method only targets on analyzing the stability state of couple machines, neglecting other pairs. We prove that EAC strictly holds for couple machines and the power-angel curve (also known as Kimbark curve [11], [12] ) of couple machines has a distinctive quasi-sinusoidal feature. In light of this feature, we aim to adopt EAC to assess stability index of couple machines based on the prediction of the Kimbark curve of couple machines, then compare the method with timedomain simulations and other classical transient energy methods. Contributions of this paper are summarized as follows: (i) The motion of a critical machine in COI reference is proved to be identical to that of a two-machine system consisting of a critical machine and a virtual COI machine.
(ii) Definition of COI is avoided in the proposed method as the virtual COI machine in PEF is replaced with a real machine in the system.
(iii) Following the thinking of PEF, the method only focuses on analyzing the stability state of couple machines, neglecting other pairs.
(iv) Strict EAC for couple machines is validated.
(v) The problem regarding critical energy computation in PEF method is overcome by using the predicted Kimbark curve of couple machines.
(vi) The method does not rely on equivalence, simplification or aggregation of the whole network or any other machines during stability analysis.
The remaining paper is organized as follows. In Section II, classic PEF method is revisited and discussed. In Section III, EAC of couple machines is analyzed. In Section IV, types of actual Kimbark curves of couple machines are defined and also predictions of these types of Kimbark curves are presented. In Section V, stability measures based on the EAC of couple machines are given. In Section VI, the selection of couple machines and supplementary details are discussed, and the the procedures of the stability analysis by using the proposed method are provided. In Section VII, the stability assessment with the proposed method is compared with time-domain simulation and classical transient energy methods. Conclusions and discussions are provided in Section VIII. Since the proposed method is genuinely derived from PEF method, clarifications about PEF are given in advance in this paper as some concepts in [11], [12] are expressed in an unsystematic and tutorial form. It is noted that "first swing" in this paper corresponds to the swing of couple machines rather than the global first-swing, also only first-swing stability of couple machines is discussed in this paper.

A. SVCS in PEF Method
Trajectory stability theory states that power system transient stability is a trajectory stability problem rather than Lyapunov problem. In other words, transient instability of the system is identical to the "separation" among machines in the system. To describe such separation, COI is designed to depict aggregated motion of all machines in the system or groups. Both the group separation method and PEF method are based on trajectory stability theory. However, the difference between them is that the group separation method states that the instability of the system is identical to the separation between the critical group and noncritical group [5], while PEF method believes that the instability of the system is identical to the separation of a single machine with respect to the system, i.e. COI of the system [11], [12]. In PEF, "a critical machine goes unstable" should be exactly expressed as "a critical machine in COI reference goes unstable". If we take a deep insight into the depiction of the critical machine in PEF method and IMEF method, one can find that the motion of a critical machine in COI reference is identical to the relative motion of a critical machine with respect to the virtual COI machine. In the following analysis, the mathematical model is based on the classical model [9]. For a n-machine system with rotor angle δ i and inertia constant M i , the motion of a single machine i (no matter it is a critical machine or not) in the synchronous reference is governed by differential equations: where voltage behind transient reactance of machine i B ij , G ij transfer susceptance (conductance) in the reduced bus admittance matrix.
Position of the COI of the system is defined by: From (2), the motion of COI is determined by: Eqn.
(3) reveals that COI can also be seen as a virtual "machine" with its own equation of motion being described as the aggregated motion of all machines in the system. Following (1) and (3), since machine i and COI are two "single" machines with interactions, a two-machine subsystem can be formed by using these two machines, which is defined as a SVCS, as shown in Fig. 1.
Substituting (3) into (1), the relative motion between the single machine and virtual COI machine of the SVCS can be given as: where From (4), one can see that the motion of each single machine (no matter it is a critical machine or a non-critical machine) in COI reference is fully identical to the relative motion of a single machine with respect to the virtual COI machine in a SVCS.
Since SVCS is a two-machine subsystem, EAC strictly holds in the Kimbark curve of a critical machine in COI reference without any K.E. or P.E. corrections [11]- [14].

B. Split of the Stability of the System in PEF Method
In transient stability analysis, critical machines are the machines that are severely disturbed by a fault and they tend to accelerate or decelerate from the rest of the machines in the system, depending on the nature of the fault. PEF analysts believe that transient behaviors of critical machines have dominant effects to the stability of the system, viz, only critical machines might go unstable and cause system to go unstable, non-critical machines are always stable and have no contribution to the instability of the system [8]- [14]. Thus, PEF only focuses on analyzing the stability state of critical machines neglecting that of non-critical machines, which makes this distinctive method non-globally. If we observe the motion of a critical machine in COI reference from the angle of SVCS, the transient stability of a multi-machine system with n machines is identical to that of a system consisting of n SVCSs. Since SVCSs with non-critical machines have no contribution to the stability of the system, the stability of the system is fully decided by SVCSs with critical machines, which  implies that the stability analysis of the system in the multidimensional space could be fully split to low-dimensional space problems. Stability mechanism of the PEF method and IMEF method is shown in Fig. 2. From Fig. 2, "a critical machine in COI reference goes unstable" in PEF method and IMEF method should be further expressed as "a pair of machines consisting of a critical machine and the virtual COI machine goes unstable", as shown in Fig.  3 (a). Although transient interactions among all machines in the system are quite complicated, the multi-machine interactions are fully implicated in the relative motion of each SVCS. Thus, the stability of the system can be judged by that of each SVCS with a critical machine independently, which is the essence of PEF method. Following the thinking of Stanton [11], [12], transient stability principle of PEF method can be described as below: (i) The system can be judged as stable when all SVCSs with critical machines are stable.
(ii)The system can be judged as unstable as long as only one SVCS with a critical machine is found to go unstable.

C. Relative Motion of a Pair of Machines
Assume two machines, i and j, are brought out from the multimachine system. A two-machine subsystem consisting of machine i and machine j is shown in Fig. 4.
The swing equations of the pair of machines are: Rewrite (5) in an angle-difference form, we have: where Mi+Mn P ej Eqn. (6) is the swing equation of the two machines. Similar to SVCS, machines i and j also form a two-machine subsystem. It is noted that (6) merely is a mathematical transformation of (5) without any simplification, aggregation or equivalence of machines or network. From the viewpoint of PEF, the swing equation of two machines underlies the stability information of the whole system during the post-fault period as both P e ij and ω ij are impacted by dynamic interactions of all machines in the system. In addition, all parameters in (6) are defined in a synchronous reference without definition of COI.

D. Strict EAC of Couple Machines
In the following analysis, we prove that EAC strictly holds in the two real-machine subsystem along the actual system trajectory.
From (6), we have: Integrating the actual fault-on trajectory till fault clearing, we have where P (F ) eij corresponds to P e ij during fault-on period. Eqn. (8) can be further expressed as: Integrating the actual post-fault trajectory after fault clearing, we have: where P ( eij P F ) corresponds to P e ij during post-fault period. Eqn. (10) can be further expressed as: Substituting (9) into (11) yields: Eqn. (12) validates that EAC strictly holds in the two-machine subsystem that is formed by machine i and j, no matter these two machines are couple machines or not. Following the actual simulated system trajectory, EAC of an unstable couple machines is shown in Fig. 5.
From Fig. 5, both integral parts in (12) can be seen as "areas", which is quite similar to the EAC of the classic OMIB system. However, the EAC of couple machines is a "precise" EAC that is based on the actual system trajectory without any K.E. or P.E. corrections, rather than the simplified EAC of the OMIB system or equivalent EAC of the EEAC system.
In following sections, EAC characteristics of couple machines are explicitly analyzed as the Kimbark curve of couple machines has a distinctive "accelerating-decelerating" characteristic compared with that of non-couple machines. Notice that the analysis of Kimbark curves of couple machines is very simliar to that of critical machines in COI reference in the PEF method.

E. Kimbark Curve of Unstable Couple Machines
The separation of couple machines after fault is cleared can fully be reflected in the Kimbark curve. From numerous simulations, a typical Kimbark curve of unstable couple machines is shown in Fig. 6.
From Fig. 6, couple machines firstly accelerate from P1 to P2 during fault-on period. After fault clearing, couple machines decelerate from P2 to P3 during post-fault period. Once the system trajectory crosses P3, machine i would accelerate and separate with respect to machine j, causing couple machines go unstable. Hence, P3 can be defined as dynamic liberation point of couple machines (CMDLP) [12], [12], and Peij would intersect with Pmij at CMDLP if couple machines go unstable. The Kimbark curve of unstable couple machines has a clear "accelerating-decelerating-accelerating" characteristic.
The unstable couple machines can be characterized by the occurrence of the CMDLP while velocity is positive.
Following EAC, when couple machines go unstable, we have: where CMDLP is an n-dimensional point that lies in the actual post-fault system trajectory although CMDLP is described in an angle-difference form. The acceleration area, deceleration area and CMDLP of couple machines would all vary with different faults. CMDLP is used to describe first-swing instability of couple machines rather than that of the system. At CM DLP ij , the mismatch of P e ij of the couple machines formed by machine i and j is zero while mismatch of other pairs of machines is not zero. If more than one pair of couple machines go unstable after fault clearing, each couple would correspond to a unique CMDLP.

F. Kimbark Curve of Stable or Critical Stable Couple Machines
When fault is not severe and couple machines are stable, the angle difference of the couple machines would initially advance after fault clearing and finally reach its maximum. In this case the velocity of the couple machines would be zero at the instant of maximum angle difference and Kimbark curve of the couple machines would not cross CMDLP. In other words, P e ij would not intersect with P m ij . Instead, it might turn upward or turn downward because the oscillation of other machines impedes backtracking of trajectory. Under this circumstance, the couple machines would be "stable". A typical Kimbark curve of stable couple machines is shown in Fig. 7.
As shown in Fig.7, after fault is cleared, couple machines decelerate from P2 to P3 and the velocity of couple machines would be zero at P3, which means that couple machines is firstswing stable. Hence, P3 with zero velocity can be defined as dynamic stationary point of couple machines (CMDSP) [12], [12]. From the figure, it can be seen that the Kimbark curve of stable couple machines has a clear "accelerating-decelerating" characteristic.
The stable couple machines can be characterized by the occurrence of the CMDSP while deceleration power is positive.
w ij = 0, P eij − P mij > 0 (15) Following EAC, when couple machines go unstable, we have: where )dδ ij In the Kimbark curve, CMDSP is the inflection point where Peij turns upward or downward. Similar to CMDLP, CMDSP is also an n-dimensional point that lies on actual post-fault system trajectory and varies with faults. CMDSP is used to describe first swing stability of a pair of couple machines rather than that of the system. At CMDSPij, the velocity of the couple machines is zero while velocities of other couple machines are not zero. For more than one pair of couple machines that are stable after fault clearing, each pair of stable couple machines would correspond to its unique CMDSP.
After fault is cleared, some couple machines might go unstable while other couples might be stable, which means that CMDSPs and CMDLPs may occur one after another along the actual postfault system trajectory.
Following stable and unstable characterization of couple machines, the critical stable couple machines are characterized by: The Kimbark curve of critical stable couple machines is shown in Fig. 8. It is noted that Eqn. (17) describes an ideal critical stable state of the couple and may not hold in actual simulated cases due to the integral errors in time-domain simulations.
From analysis above, one can merit that the Kimbark curves of stable couple machines and unstable couple machines both have a clear "accelerating-decelerating" characteristic after fault occurs. Since the stability of the system is decided by that of couple machines, the analysis of the Kimbark curve of couple machines is of key importance in transient stability analysis. In fact, the most distinctive characteristic of the Kimbark curve of the couple machines is that it is quasi-sinusoidal and is predictable, which will be analyzed in Section IV.

A. Quasi-Sinusoidal Characteristic of the Kimbark Curve of Couple Machines
Assume Ω cr and Ω ncr are sets of critical machines and noncritical machines after fault clearing, respectively. Machine i is a critical machine that lies in Ω cr and machine j is a non-critical machine that lies in in Ω ncr , then machine i and j would form couple machines as their ω ij is much larger than those pairs of machines brought out from the same set. The power output of P ei and P ej is: Then we have: In (19), δ in and δ jm can be expressed as: (20), δ jn and δ im are angle differences to describe motions inside each set as i, m ∈ Ω cr and j, n ∈ Ω ncr . Substituting (20) into (19) yields: Features of K sin , K cos and K tail in (21) are analyzed as follows: (i) K sin , K cos and K tail signify the impacts forced on the couple machines that come from all the machines in the system although P eij is derived from only two machines.
(ii) As majorities of post-fault trajectories of machines are still in synchronism, ω ij of couple machines is much larger than ω im and ω jn . Such fact implies that ∆δ ij would be much larger than ∆δ im and ∆δ jn in a short time interval, thus most components in K sin , K cos and K tail barely change compared with violent variations of ∆δ ij during the post-fault period. In term of dynamic interactions among all machines of the system, the quasi-stationary features of K sin , K cos and K tail are maintained due to the strong stationary effect of the majority of non-critical machines during the post-fault period. For power system transient stability analysis, a commonly used simplification for post-fault trajectory calculation is that δ im and δ jn are all constants. By applying this simplification, K sin , K cos and K tail would be constant for the couple machines, thereby Kimbark curve of couple machines becomes an ideal sinusoidal curve. However, such simplification is based on the ideal assumption when fault is severe. For most actual trajectories when the fault is not severe and the system still keeps stable, the rotor motions of machines within each set would vary (δ im and δ jn cannot be assumed constant) irrespective of the dynamic response of machines and loads, which means that Ksin, Kcos and Ktail would be effected by the "non-sinusoidal" feature incurred by the dynamic response of the system. Thus the shape of Kimbark curve of couple machines can be described as "quasi-sinusoidal".
(iii) From numerous simulations, the quasi-sinusoidal feature of Kimbark curve of couple machines still holds when the separation pattern of machines in the system is complicated.

B. Parameters Identification
According to the distinctive quasi-sinusoidal feature, Kimbark curve of couple machines can be predicted by using the parameter identification of a formula combined with sine and quadratic functions. The sine function and quadratic function are used to describe the sinusoidal and non-sinusoidal features of Kimbark curve of couple machines, respectively, and can be written as: (22) are unknown parameters to be identified. If H q1 and H q2 is set to be zero, the expression in (22) would be an ideal sinusoidal curve P (sin) eij : Parameters in (22) can be identified by using the least-square approximation method, with data in a short sampling window after fault is cleared. The equation for parameter identification of P (qr) eij is: where M in (24) should be at least five dimensions as H iden has five parameters to be identified. Generally the starting time point for sampling can be set as the fault clearing point, and 10 to 15 samples within 10 ms is appropriate for identification for both small and large scale systems (more than 100 machines). The identification of P  In (23), only the coefficient of sinδ ij is defined by quadratic form. Theoretically, both H cos and H cst can also be replaced with quadratic forms. However, since the primary characteristic of Kimbark curve is sinusoidal, one quadratic form is efficient to describe non-sinusoidal feature of Kimbark curve.
Once H iden are identified, the Kimbark curve being expressed with P (qr) eij can be predicted. It is noted that all parameters for the predictions of Kimbark curve, i.e. P eij in P mea and δ ij in M are all local information of couple machines as COI is not used to describe couple machines, which means that the prediction of Kimbark curves of couple machines is independent of the information of other machines.
In actual simulation cases, although the quadratic form in (22) can well catch the non-sinusoidal feature of Kimbark curve within the sampling window, the identification error of the quadratic form might be enlarged with the increase of simulation time. To ensure robustness of the identification, the effect of quadratic form in (22) should be slightly weakened, thus following weighted formula is utilized.
Step 3 Set H q1 = and H q2 as constants and re-identify P (qr) eij .P (qr) eij , then H q3 , H cos and H cst can be obtained.

C. Typical Types and Predictions of Actual Kimbark Curves of Couple Machines
According to actual simulated Kimbark curves, if couple machines are stable during the first swing, the Kimbark curve might "bend down" or "bend up" in the first swing. After that, the curve might turn upward or turn downward during the second swing, as shown in Fig. 9. e fault clearing time the corresponding Kimbark curves can be classified into four types as shown in Fig. 9. The detailed illustrations are shown in Figs. 10-13. The actual simulated Kimbark curves are drawn with solid lines with arrows directing actual power trajectories.
From Fig. 10 to Fig. 13, common features of all types of actual Kimbark curves are described as below.
(I) For all types of curves, the non-sinusoidal feature would be reflected when couple machines are stable.
(II) For all types of curves, the non-sinusoidal feature fades and the sinusoidal feature gradually dominates with the increase of (III) All types of curves intersect with P mij when couple machines go unstable. The intersection point is the CMDLP as analyzed in Section III D.
(IV) All types of curves do not interest with P mij once the system trajectories reach CMDSP as analyzed in Section III E.
Predictions of the actual Kimbark curve by using P  are uncertain when couple machines are critical stable or critical unstable. However the approximation error is very slight.
From simulations above, P eij has a better approximation than P (sin) eij in most cases. Although P (sin) eij might be slightly better than P (sin) eij in a few cases, the difference generally has no effect to stability assessment. Thus P (qr) eij is more closer to the actual Kimbark curve than P (sin) eij in most cases. In fact, the predictability of Kimbark curve of couple machines is not a coincidence. Following the thinking of PEF, complicated interactions among all machines in n-dimensional space is split into that in m two-dimensional spaces after fault clearing if m pairs of couple machines exist after fault clearing. Thus, the complexity of the stability analysis can be greatly simplified by reducing the dimensions of problem solving, making Kimbark curve of couple machines sine-regular and predictable.

A. Predicted Deceleration Area of Couple Machines
The predictability of the Kimbark curve means that the utilization of EAC of couple machines becomes feasible. Typical predictions of deceleration area of stable and unstable couple machines are shown in Figs. 14 (a) and (b), respectively.
From Fig. 14, the corresponding predicted CMDLP (PCMDLP) is given as the intersection point between the predicted Kimbark curve and the horizontal mechanical power curve. Thus the predicted deceleration area can be expressed as: Stability margin of the couple machines is defined as: In (26), η ij > 0 means the couple is stable, η ij < 0 means the couple is unstable, and η ij => 0 means the couple is critical stable.
The definition of η ij can be validated by Kimbark graphics. From Fig. 25 (a), if couple machines are stable, the PCMDLP would be far from CMDSP and A (P red) ACCij would be larger than the actual deceleration area A DECij , which means that η ij would be positive when couple machines are stable. From Fig.25(b), since P (qr) eij would almost coincide with the actual Kimbark curve when couple machines go unstable as analyzed in Section IV C, the PCMDLP would be very close to the actual one and A (pred) DECij would also be close to the actual A DECij under this circumstance. That is to say, the predicted η ij would be negative and be close to the actual stability margin when couple machines go unstable. Thus, the definition of η ij in (26) is fully in accord with actual simulated stable and unstable cases of couple machines.
Since many types of actual Kimbark curves of couple machines may exist as analyzed in Section IV C, the predictions of Kimbark curves and PCMDLP of couple machines are analyzed in categories in the following sections.
Following different types of predicted Kimbark curves, the PCMDLP of each type is defined as follows.
PCMDLP of A-1: PCMDLP is defined as the intersection point between P mij and P

C. Stability Measure of Couple Machines 1) Stability Measures of A-1 and A-2:
From the Kimbark curve as in Fig. 14, η ij can be computed by using a simple procedure as follows: Step 1 Compute A ACCij by the trapezoidal integration.

The computation of A (pred)
DECij is shown in Fig. 15.

2) Stability Measures of A-3:
Stability measure of A-3 is very similar to that of A-1, the only difference is that PCMDLP is computed by using the intersection point between P

3) Stability Measures of A-4:
Computation of stability index of A-4 is special as PCMDLP does not exist. Theoretically, the predicted deceleration area can be seen as infinity in this case and the couple machines are "infinite" stable. However, as angle difference between two machines cannot be too large, an imperial PCMDLP should be given in advance. In this paper, PCMDLP of A-4 is set as π.
Stability measure of A-4 is defined as: In (28), the minimum margin of couple machines should be set as zero to prevent a negative η ij because couple machines are predicted to be stable in this case.

D. Stability Measure of the System
Multi-pairs of couple machines may exist simultaneously after fault clearing. Since stability of the system is observed by multicouples with the proposed method, the margin of the system is also defined as a multi-dimensional vector that is formed by margins of all couples.
In (29), Ω c is the set of all couple machines. All η ij > 0 means the system is stable. One or a few η ij = 0 while the rest of η ij > 0 means system is critical stable. One or more η ij < 0 means system goes unstable.

A. Identification of Couples for a Certain Fault
For the online TSA or emergency control, the system operators may wish to grasp all severely disturbed machines in the system, thus it would be necessary to identify all couples at fault clearing point under this circumstance, and following key characteristics of couple machines can be utilized: (i) ω of couple machines is much higher than that of noncouples at fault clearing point.
(ii) The Kimbark curve of couple machines has a clear "accelerating-decelerating" characteristic after fault occurs. The identification of couple machines in a small scale system may differ from that in a large scale system. Since the small system has a very weak capability to accommodate fault disturbance, all machines in a small system might be severely disturbed by faults and accelerate simultaneously in synchronous reference. Compared with small scale system, a large scale system preserves a much stronger capability to accommodate fault disturbance. To be specific, only a few machines around fault location might be disturbed while most non-critical machines that are remote from fault location may oscillate slightly after fault occurs. Thus, the identification of couple machines in a large scale system would be simpler than that in a small scale system.
From analysis above, the identification strategies of couple machines in a small scale system and a large scale system are given separately.

a) Identification of couple machines in a small scale system
Step 1 Sort rotor speed of all machines in the system in a descending order at fault clearing point.
Step 2 Assume the sorted index is S = [s 1 , s 2 · · · s N ]. Form the set of test pairs [s 1 s N , s 2 s N −1 , · · · s 1+q s N −q ]. q satisfies following: ω s1+q+1 sN−q−1 < ω cthr , ω s1+q sN−q > ω cthr (30) Step 3 Combine first q machines and last q machines to form the set of candidate couple machines Ω cc . All combinations of candidates are q 2 pairs in Ω cc .
Step 4 Verify candidate couples in Ω cc before sampling window is over. The candidates with clear "accelerating-decelerating" Kimbark curves are defined as couple machines and be preserved in a final set Ω c .
In (30), ω cthr is set as a threshold of selecting candidate couples. ω cthr is given as 0.002 p.u.-0.005 p.u. for small scale systems based on the experiences of numerous simulations.

b) Identification of couple machines in a large scale system
Step 1 Select machines that are close to the fault location and form the set of severely disturbed machines ω se at fault clearing point. Machines in Ω se satisfy: ω i > ω sethr i∈Ω sr (31) Step 2 Select a machine that is remote from the fault location as machine j. ω j should be close to zero. Machine j can be seen as the representative machine of all stationary non-critical machines in the system.
Step 3 Combine machine j with machines in Ω se to form the set of candidate couples Ω cc at fault clearing point.
Step 4 The candidates with clear "accelerating-decelerating" Kimbark curves in Ω cc are defined as couple machines and are preserved in Ω c after sampling window is over.

B. Identification of Couples for Computation of CCT
The identification of couple machines when computing CCT is different from that for a certain fault. The reason is that system operators do not care about all couples under this circumstance but only target the identification of critical stability of the system. Following the analysis in Section III A, the critical stability of the system is identical to the case that the lead couple is critical stable while other couples are marginal stable, thus the system operators can only focus on the lead couple neglecting other stable couples with large margins. In actual simulations, system operators may observe those a few most severely disturbed couples to ensure that the lead couple is involved. Under this circumstance, Ω c can be formed with only a few elements, which greatly simplifies the identification of couples when computing CCT.

C. General Procedures of the Proposed method
General procedures by using the proposed method in TSA are outlined as below.
Step 1 Select couple machines and form Ω c after fault clearing.
Step 2 Kimbark curve of each couple machines is predicted by using P (qr) eij and P (sin) eij simultaneously in a short sampling window.
Step 3 For each pair of couple-machines in Ω c , compute A ACCij and A DECij .
Step 4 η ij of each couple machines is calculated using EAC of couple machines.
Step 5 η sys is calculated with η ij of all couple machines, then the stability of the system is determined.

A. Test Network
The proposed method is firstly tested on a small scale system, i.e., a modified 10-unit New England system (the inertia constant of Unit 1 is modified to 200 pu from 1000 pu). Three simulation cases are described below.
Case-1: Three phase short circuit event occurs at bus 34 at 0.00s and is cleared at 0.23s (the system goes unstable).
Case-2: Three phase short circuit event occurs at bus 4 at 0.00s and is cleared at 0.50s (the system goes unstable).
Case-3: Three phase short circuit event occurs at bus 4 at 0.00s and is cleared at 0.03s (the system is well stable).
Swing curves of machines in Case-1 and Case-2 are shown in Fig. 16. All parameters are set in synchronous reference.

B. Case-1
Rotor speed of machines at fault clearing point is sorted in descending order, as shown in Table I.
According to the procedure of identification of couple machines in a small scale system as in Section VI, q is set as 3. Ω cc and Ω c are shown in Table II. Sampling window for approximating Kimbark curve of couples in Ω c by using both P   types and the margin of couples are shown in Table III. Kimbark curves of couples with machine 36 are very similar to those couples with machine 33 and thus are not shown here. From Table III, couple 34 30 is the lead couple among all couples. The instability of the system is decided by unstable couples 34 30, 34 31 and 34 39 because system can be judged as unstable by finding only one unstable couple. However, the severity of the system is only decided by the lead couple 34 30 and η sys is given as -0.79. As shown in Table II, there are total 9 pairs of couples in this case. However, if classical PEF method is used in Case-1, then only 3 critical machines, i.e., machines 34, 33 and 36 in COI reference could be identified as critical machines, thus the number of critical machines with classical PEF is smaller than the number of couples with the proposed method. The reason is that the virtual COI machine is set as the only stability reference in the classical PEF method, while 3 real machines in the system (machine 30, 31, 39) are set as stability references of couples for the proposed method, which increases the combinations of couples.

C. Case-2
To further illustrate the proposed method, Case-2 with a complicated machine-separation pattern is given. Ω cc and Ω c are provided in Table II. From simulation, the Kimbark curves of all couples in Ω c are also very close to ideal sinusoidal although the identification of critical machines is complicated in this case. Stability judgement of unstable couple machines in Case-2 is also shown in Table III. From Table III, although the separation pattern of machines in Case-2 is more complicated than that of Case-1, the number of unstable couples is even larger than that of Case-1, and the system can be easily determined as unstable.

D. Case-3
Case-3 is designed to test an extreme well stable case using the proposed method. The maximum rotor speed ω 32 at fault clearing point is only 0.0022 p.u. in this case, which indicates that system is slightly disturbed by faults. Ω c with only one couple is shown in Table II. Margin of the couple is shown in Table III. Kimbark curve of couple 32 39 is shown in Fig. 17.
From Fig. 17, the quasi-sinusoidal feature and "acceleratingdecelerating" characteristic of actual Kimbark curve of couple machines are both preserved in the well stable case. The reason is that ω 32 39 is still several times larger than the rotor speed difference of non-couple machines even though ω 32 39 is only 0.0022 p.u in this case. Such results also validate that the proposed method can be utilized in well-stable cases effectively without over conservative judgements.

E. Computation of CCT
Different fault locations in the modified New England system are simulated to evaluate the accuracy of the proposed method, and the test results are shown in Table IV. From Table IV, the proposed method has a high accuracy in identifying CCT. The tiny error is generally incurred by approximation errors of actual Kimbark curve.
The computation of CCT at bus 34 is taken as an example to demonstrate the difference between the selection of couple machines when computing CCT and that for a certain fault in Case-1 and Case-2. With the increase of the number of iterations, it is becoming evident that couples with machine 34 are most severely disturbed couples and have leading effects to the instability of the system. Thus, as t c approaching the real CCT, the proposed method only focuses on couples 34 30, 34 31 and 34 39 by neglecting other stable couples with large margins. By using the proposed method, couples 34 30, 34 31 and 34 39 are all stable lead couples when t c is 0.18s and all go unstable when t c is 0.19s, thus the CCT is set as 0.18s, which fully reflects the distinctive non-global characteristic of the proposed method.

1) Comparison between CUEP Method and Proposed Method:
Two critical unstable cases of the standard IEEE 118-bus 20-unit test system in [17] are re-simulated to demonstrate the difference between the proposed method and CUEP method. The simulated critical-unstable fault clearing time is slightly different from that in [17].
Case-4: Three phase short circuit event occurs at the terminal of Gen. #3 at 0.00s and is cleared at 0.51s.
Case-5: Three phase short circuit event occurs at the terminal of Gen. #2 at 0.00s and is cleared at 0.20s. a) CUEP method For Case-4, possible MODs [Gen. #2, Gen. #3, Gen. #2 and Gen. #3] are given first, then UEP of each MOD is computed by solving a non-linear optimization problem. The MOD with Gen. #2 is finally defined as the dominant MOD because its corresponding normalized ∆V P E is the smallest among all possible MODs. The CUEP in the dominant MOD is (θ 2,CUEP = 2.183rad, and global critical energy is 3.468 p.u. For Case-5, the computation procedure is the same with that in Case-4. It is worthy to point out that both Case-4 and Case-5 have the same CUEP and same amount of global critical energy although both the fault location and fault clearing time in Case-5 are different from that in Case-4 [17]. b) Proposed method For Case-4, Ω c is given in Table II. Couples 2 18, 2 19 and 2 20 go critical unstable while couples 3 18, 3 19 and 3 20 are stable, then system is judged to go unstable. Kimbark curves of couples 2 20 and 3 20 are shown in Figs. 18(a) and (b), respectively. Notice that the predictions of Kimbark curve of couples are both steady and robust, which technically avoids nonconvergence problems of UEP computation in the CUEP method [17].
From Fig. 18, Kimbark curves of couples 2 20 and 3 20 indicate that the critical instability of the system is fully decided by the critical instability of lead couples 2 18, 2 19 and 2 20 (2 18 and 2 19 are not shown in the figure) rather than all couples. In fact, if the system operators only want to know system is stable or not in this case, the computation of the proposed method can be terminated immediately as long as only one couple with machine 2 is found to go unstable and analysis regarding to couples with machine 3 can be neglected. Consequently, the identification of possible MODs required in CUEP method can be avoided when judging stability of the system as it can be judged in a non-global manner by using the proposed method. However, if system operators want to obtain both stability and severity of the system, the identification of the lead couple will be needed under this circumstance and stability of all couples (couples with machine 2 and machine 3) should be judged, which becomes similar to the identification of MOD that needs to judge stability state of all critical machines.
For Case-5, the identified couples are same with that in Case-4. However, the computed deceleration areas and PCMDLPs of couples in both cases are all different, which shows that the proposed method is more sensitive to the fault changes compared with CUEP method.
Since the proposed method is a hybrid method that relies on actual simulated Kimbark curves, the computation efficiency of the proposed might be lower than that of CUEP method in some small scale systems. However, the computation reliability of the proposed method would be revealed with the increase of system scales, as the computation of CUEP would become unbearable in large scale systems especially during online TSA.

2) Comparison between IEEAC Method and Proposed Method:
To test the robustness of the proposed solution, three abnormal operating conditions during optimization are tested. The abnormal operating conditions include loss of communication link, disconnection of load, and loss of agent, which would produce detrimental effects if a centralized method were used.
A practical 2766-bus, 146-unit interconnected system is given to compare the difference between the proposed method and IEEAC method. SYSTEM LC is a regional system consisting of 8 units while SYSTM SD is a main system consisting of 138 units. SYSTEM LC and SYSTM SD are connected by two 500KV The simulated case is given as below. Case-6: The fault location is set at 90% of Line LIAOC TANZ that is close to TANZ. Three phase short circuit event occurs at 0.00s and is cleared at 0.22s (CCT is 0.16s in this case).
Swing curves of machines in the interconnected system is shown in Fig. 20. a) IEEAC method To simplify expression in IEEAC method,Ω n LC is the set with [Gen. #XW, Gen. #HY, Gen. #CP, Gen. # XYRD], Ω s LC is the set with [Gen. #LC1, Gen. #LC2, Gen.#LCR1, Gen. #LCR2], and Ω SD is the set with all machines in SYSTEM SD. From  Fig. 20, all machines in SYSTEM LC accelerate with respect to SYSTEM SD after fault is cleared, thus the separation of all machines in the whole interconnected system is a clear and favorable inter-area oscillation mode for IEEAC method.
Using IEEAC method, all machines in the interconnected system after fault clearing are separated as the critical group Ω A and the non-critical group Ω S . From the swing curves of machines in Fig. 20, possible group separation patterns are given as below.
Pattern-1: Ω A is set as Ω n LC and Ω S is set as Ω s LC ∪ Ω SD . Pattern -2: Ω A is set as Ω n LC ∪ Ω s LC and Ω S is set as Ω SD . Pattern-2 is finally set as the dominated group separation pattern because stability margin of the OMIB system in Pattern-2 is much smaller than that in Pattern-1. For Pattern-2, 8 machines in Ω n LC and Ω S and 138 machines in Ω SD are aggregated as Machine-A and Machine-S, respectively. The Kimbark curve of aggregated OMIB system in Pattern-2 is shown in Fig. 21. The equivalent P m in the figure is not a horizontal line as governor is deployed.
From Fig. 21, one can see that the Kimbark curve of the OMIB system goes across dynamic saddle point [6] which lies in actual post-fault trajectory at 0.37s, thus the interconnected system is judged to go unstable and inter-area oscillation is identified. b) Proposed method Couple identification strategy in the large scale system is utilized first using proposed method. Targeting machines that are close to fault location, Ω se is formed by [XW, HY, CP, XYRD, LC1, LC2, LCR1, LCR2]. Gen. #LZ in SYSTEM SD (not shown in Fig. 19) that is quite remote from the fault location is selected as machine j. Gen. #LZ can be seen as the representative machine of the main SYSTEM SD. Ω cc and Ω c are shown in Table II.
From Table II, unlike IEEAC method that separates all machines in the interconnected system into two groups, the proposed method only focuses on the stability of 8 couples in Ω c that are formed with generators in SYSTEM LC and Gen. #LZ in SYSTEM SD. Using proposed method, all couples in Ω c are judged to go unstable, which reveals that all machines in SYS-TEM LC separate with respect to the same machine Gen. #LZ in SYSTEM SD, thus the inter-area oscillation is substantially identified. Kimbark curves of two unstable couples HY LZ and LCR1 LZ in Ω c are shown in Fig. 22 (a) and (b), respectively.
From Fig. 22, the quasi-sinusoidal characteristic is also fully reflected by Kimbark curves of couples in the large scale system and is well captured by the proposed method, which validates the analysis in Section IV. Since all machines in regional SYSTEM LC are severely disturbed, it is also worthy for system operators to obtain the local-oscillation situation inside SYSTEM LC although inter-area oscillation has already been identified. Using Gen. #LC1 as the reference machine of SYSTEM LC, swing curves of all machines in regional SYSTEM LC are shown in Fig. 23.
From Fig. 23, local oscillation also occurs in regional SYS-TEM LC. Using the proposed method, 16 couples are formed inside SYSTEM LC following the couple identification strategy in the small scale system, as shown in Fig. 23. Kimbark curves of two representative unstable couples HY LCR1 and HY LCR2 inside SYSTEM LC are shown in Fig. 24 (a) and (b), respectively.
From simulations, all couples inside SYSTEM LC go unstable and the local oscillation inside regional SYSTEM LC is also identified by the proposed method.
An extreme case is also given to further demonstrate the non-global characteristic and local-oscillation identification of the proposed method in the online TSA. Assume the anglemeasurements of all machines in SYSTEM LC except Gen. #HY and Gen. #LCR1 malfunction, then the rotor angles of most machines in SYSTEM LC become "unobservable" for the system operators after fault occurs, leaving only Gen. #HY and Gen. #LCR1 "observable". Under such extreme circumstance, the interarea oscillation could not be identified by the proposed method (or any other methods) as rotor angles of most machines in SYSTEM LC are unknown. Meanwhile, η sys is also unknown as system operators do not know whether couple HY LCR1 is the lead couple or not. However, the local oscillation inside SYSTEM LC still can be identified by the proposed method as the only "observable" couple HY LCR1 can be easily judged to go unstable. c) Computation efficiency analysis Compared with IEEAC method that can only identify interarea oscillation in this case, both inter-area oscillation and local oscillation can be identified using the proposed method. In addition, two parallel computations are deployed simultaneously in the proposed method.
(i) Inter-area oscillation mode and local oscillation mode are identified in parallel.
(ii) In each mode, the stability of each couple is judged in parallel.
Since the selection of couples and parameters identification of Kimbark curves in both oscillation modes are simple, the corresponding time consumed can be neglected and the total time consumed is just the length of sampling window for each couple within 150 ms after fault is cleared using the proposed method.
In addition, compared with IEEAC method that relies on transient information of all machines in the interconnected system, local information of only 9 machines is used for online TSA using proposed method without any aggregation or equivalence. To be specific, the local information of each machine for TSA is δ i and P ei from the pre-fault point to the end point of the sampling window, and ω c i at fault clearing point. From simulations above, one can find that different oscillation modes require different identification strategies of couples using the proposed method. Then one question emerges: Could a couple inside SYSTEM LC also be defined as a couple in the interconnected system? In fact, the formation of couple machines is an intrinsic nature of the transient stability and is fully independent of definition of regional or interconnected systems. For instance, if couple HY LCR1 is a couple in SYSTEM LC, then it is also a couple in the interconnected system as SYSTEM LC is a regional system inside the interconnected system. That is to say, SYS-TEM LC and the interconnected system can both be judged to go unstable once couple HY LCR1 goes unstable. However, for the TSA of the interconnected system, the system operator would only focus on the inter-area oscillation to grasp the whole vision of the instability of the interconnected system, thus HY LCR1 is theoretically a real couple but is methodologically neglected in the couple identification strategy of the interconnected system because the instability of HY LCR1 has no contribution to the identification of the inter-area oscillation. d) Discussion of the relationship between IEEAC method and proposed method IEEAC method is fully based on trajectory stability theory as the Kimbark curve of the equivalent OMIB system and the dynamic saddle point are fully obtained from actual simulated trajectory, which means that the stability margin and dynamic saddle point of the OMIB system would vary immediately with the slight change in fault location or fault clearing time. Compared with other classical global methods, IEEAC method is more sensitive to the fault change and could give a more accurate stability judgement, which is also the major advantage of IEEAC method. Similar to the theoretical basis of IEEAC method, the proposed method is also based on trajectory stability theory and CMDLP also lies in the actual simulated post-fault trajectory. The difference between the IEEAC method and the proposed method is the angle of observing the transient instability of the system. For the proposed method, the instability of the system is described as the separation of a few pairs of couple machines, rather than separation of groups.

VII. CONCLUSION AND DISCUSSION
This paper proposes a hybrid direct-time-domain method. The stability of a critical machine in IMEF and PEF is clarified as the stability of a pair of machines consisting of a critical machine and the virtual COI machine. Following the thinking of PEF method, the proposed method transforms the stability analysis of the system into that of a few pairs of couple machines, which avoids the equivalence or aggregation of network. The stability index of couple machines is evaluated by EAC individually and instability of the system is determined by unstable couples, which greatly reduces complexity of stability analysis. The method developed is able to efficiently evaluate the system stability and simulation results obtained are practical in the sense that they are sufficiently accurate compared with time domain simulation and classic global methods.
Further work is to focus on fundamental theories and multiswing stability judgement with the proposed method.