Extended Equal Area Criterion Revisited: A Direct Method for Fast Transient Stability Analysis

: For transient stability analysis of a multi-machine power system, the Extended Equal Area Criterion (EEAC) method applies the classic Equal Area Criterion (EAC) concept to an approximate One Machine Inﬁnite Bus (OMIB) equivalent of the system to ﬁnd the critical clearing angle. The system-critical clearing time can then be obtained by numerical integration of OMIB equations. The EEAC method was proposed in the 1980s and 1990s as a substitute for time-domain simulation for Transmission System Operators (TSOs) to provide fast, transient stability analysis with the limited computational power available those days. To ensure the secure operation of the power system, TSOs have to identify and prevent potential critical scenarios through ofﬂine analyses of a few dangerous ones. These days, due to increased uncertainties in electrical power systems, the number of these critical scenarios is increasing, substantially, calling for fast, transient stability analysis techniques once more. Among them, the EEAC is a unique approach that provides not only valuable information, but also a graphical representation of system dynamics. This paper revisits the EEAC but from a modern, functional point of view. First, the deﬁnition of the OMIB model of a multi-machine power system is redrawn in its general form. To achieve fast, transient stability analysis, EEAC relies on approximate models of the true OMIB model. These approximations are clariﬁed, and the EAC concept is redeﬁned with a general deﬁnition for instability, and its conditions. Based on the deﬁned conditions and deﬁnitions, functions are developed for each EEAC building block, which are later put out together to provide a full-resolution, functional scheme. This functional scheme not only covers the previous literature on the subject, but also allows to introduce several possible new EEAC approaches and provides a detailed description of their implementation procedure. A number of approaches are applied to the French EHV network, and the approximations are examined.


Introduction
To ensure the secure operation of the power system, Transmission System Operators (TSOs) perform offline Transient Stability Analysis (TSA) for a few dangerous scenarios and design remedial actions for the critical ones, i.e., the ones with lower Critical Clearing Time (CCT) (CCT is the maximum fault elimination time without the system losing its capability to recover a normal operating condition [1]). The reference technique for TSA is time-domain simulation using numerical integration of the nonlinear differential equations representing the system dynamic model. Time-domain simulation is flexible, and it can consider a detailed model for almost any component of the power system. However, this detail comes at a cost, a high computational time. Moreover, the time-domain simulation cannot provide a direct indication of the CCT, thus the system equations should be solved for different fault elimination times to search for the critical time. In the early 20th century, stable (in the first-swing) if after the fault the rotor angle does not increase continuously; therefore, it reaches a maximum value and thereafter decreases. It can be shown that the variations of the SMIB system rotor angle δ, is linked to the area between its input mechanical power P m and output electrical power P e in δ − P plane. As shown in Figure 1, the area has a positive portion for which P m > P e and a negative portion for which P m < P e . The first portion is linked to the energy gained during rotor acceleration, and the second portion is linked to the energy dissipated during rotor deceleration. If there is a δ m after which the sum of the areas between P m and P e becomes negative, the system is stable and at δ m the rotor angle starts to decrease. By increasing the fault elimination angle δ e , we reach an angle after which clearing the fault cannot maintain the system's stability, i.e., the sum of the different areas is always positive and the rotor angle increases continuously. This critical angle is the CCA of the SMIB system. Once the CCA has been calculated, the CCT can be obtained by numerical integration of SMIB differential equations up to CCA. Equal area criterion for a synchronous generator connected to an infinite bus. There is a δ m after which the area becomes negative and δ starts to decrease, i.e., stable case For a multi-machine power system, the idea of EEAC was proposed in the late 1980s. It belongs to a class of transient stability analysis methods that rely on the OMIB equivalent model of the power system. The OMIB-based methods are based on the observation that the loss of synchronism of a multi-machine power system originates from the irrevocable separation of its synchronous generators into two groups: the Critical Cluster of generators (CC) which push the system towards instability, and the remaining Non-critical Cluster of generators (NC). These methods replace these two groups with a two-machine system and then with an OMIB equivalent [1]. This section presents the OMIB-equivalent general equations and discusses the approximations that let us estimate the CCA of the OMIB model quickly and without a formal solution of its equations.

One-Machine Infinite Bus Concept and General Formulation
The OMIB can be considered as an approximate transformation of the multidimensional state-space of multi-machine power system dynamics to a lower dimension space. Based on the transformation approach, the resultant OMIB model can be a 'time-invariant' where P m denotes the mechanical power, P e is the electrical output power, and ω 0 = 2π f 0 .
With summation over set CR equations: Considering Equation (2), we can rewrite this equation in the following form: Similarly, with summation over set NC swing equations we have: Dividing both sides of Equations (5) and (6) by M cr and M nc , respectively, and subtracting them we get: To derive the OMIB model, we define the OMIB angle δ as follows: Multiplying both sides of Equation (7) where: P m = M nc ∑ k∈CC P m k − M cr ∑ l∈NC P m l M T (10) P e = M nc ∑ k∈CC P e k − M cr ∑ l∈NC P e l M T The OMIB stability can be inferred from Conjecture 2, given below: An OMIB is first-swing unstable, if after fault inception its angle increases continuously with time. The OMIB CCA, if it exists, is the smallest fault elimination angle after which clearing the fault cannot maintain system stability.

Extended Equal Area Criterion Concept
Let us consider a simple four-machine power system shown in Figure 2. For a shortcircuit fault in one of the transmission lines, the variations of generator angles with time can be obtained using time-domain simulation, based on which one can judge the system transient stability. For example, Figure 3 shows the angle variations of the four machines for two different fault elimination times. The solid lines show the trajectories when the fault is successfully cleared at fault elimination time t 1 e , and the generator angles recover towards stability. The dashed lines show the trajectories for a slightly larger fault elimination time t 2 e , after which the generator angles diverge continuously, and the system losses its transient stability.  The main idea behind the EEAC is to apply the classic EAC concept to the OMIB equivalent of a multi-machine power system. Multiplying both sides of the swing equation of Equation (9) by 2dδ/dt and integrating we get: (12) where δ 0 is the OMIB equivalent initial pre-fault angle. The above equation shows that the OMIB angle variation is linked to the area between P m and P e . Consider Figure 4, which shows the δ − P and t − δ curves for the OMIB equivalent of detailed model of the four-machine power system, when the fault is initiated at t = 0 and corresponding δ 0 , and cleared at t 1 e . The curves are obtained from time-domain simulation results using OMIB equations in their general form. After the fault, δ starts to increase. Based on Equation (9), as at the initial instant of the fault P m > P e , the change of δ in t − δ plane is concave-up (the rate of change increases). At t 1 e and corresponding δ 1 e , the fault is eliminated. A short moment after fault elimination, P m < P e and the curve in t − δ plane becomes concave-down (δ is still increasing, but the rate of change is decreasing). The OMIB angle δ will continue increasing until dδ/dt in Equation (12) becomes zero, i.e., the area between P m and P e becomes zero.
As shown, the area can be divided to two portions. The positive area is where P m > P e and the OMIB acceleration in Equation (9) is positive. The negative area is where P m < P e and the acceleration is negative. As shown, the negative area does not necessarily start after fault elimination. If, as in Figure 4, there is a δ m after which the sum of the areas between P m and P e becomes negative, the system is first-swing stable. The δ m at which the sum of the areas becomes zero and δ starts to decrease is the return angle δ r . Figure 5 shows the δ − P and t − δ curves for the same scenario, but with a slightly larger fault elimination time t 2 e . Similar to the previous case, a short moment after fault elimination, P m < P e and the variations in t − δ plane becomes concave-down, i.e., OMIB decelerates. However, a δ m cannot be found after which the area becomes negative. Therefore, δ increases and reaches δ u at time t u . After this point, P m > P e , the OMIB starts to accelerate again, δ increases continuously, and the system loses its first-swing stability.  It might occur that at the initial instant of the fault P m is lower than P e . In this case, δ will start to decrease. Based on Equation (9), as P m < P e following the fault, the change of δ in t − δ plane will be concave-down. At the point where P e becomes lower than P m , due to fault elimination, decreased electrical power, or increased mechanical power, the variations of δ become concave-up, but it will decrease until an angle δ m is reached, at which dδ/dt in Equation (12) becomes zero. If there cannot be a δ m which satisfies this condition, δ will decrease continuously and the system will be 'backward-swing unstable'.
The EEAC method relies on the general swing equation of Equation (9) in conjunction with the traditional EAC concept to provide fast transient stability assessment. In this equation, P m can be calculated by substituting the mechanical powers of synchronous generators of critical and non-critical sets in Equation (10). Without considering the governor controls, P m is a known parameter equal to its pre-fault value. However, P e in Equation (11) depends on the electrical power of synchronous generators, which are functions of their variable angles.
Let us consider Figure 3 which shows the angle trajectories of the four-machine power system subjected to a short-circuit fault and cleared at two different elimination times. It can be seen that for any time instant of the during-fault period, trajectories of generator angles are similar for different t e . However, for the post-fault period, they depend on t e . In other words, with an exact simulator of the system, for each generator, and for the during-fault period, we can calculate the trajectory of its angle with respect to time. However, for a post-fault period the trajectory depends on δ e , and thus on t e .
Let us consider the OMIB δ − P curves in Figure 6 plotted for the same fault scenario and the two different fault elimination times. For the during-fault period, as the variation of each generator angle with respect to time can be defined uniquely, the variations of the electrical power of each generator with respect to time can be uniquely defined. Therefore, as shown, the δ − P curves for different t e coincide in the during-fault period. However, for a post-fault period, the variations of each generator angle with respect to time depend on t e . Therefore, the variations of the OMIB P e with δ depend on the fault elimination time. Consequently, as shown, for the post-fault period the δ − P curves are different and dependent on δ e .
The problem is that we intend to apply the EAC on the OMIB model to find the critical clearing angle, as the smallest δ e after which the area between P m and P e is always negative. However, the post-fault P e of the OMIB model is uncertain and dependent on δ e . Therefore, the area itself depends on δ e . The scientific question here is how to find the critical clearing angle using uncertain δ − P curves which are dependent on the clearing angle.
In response to this question, an approach could be to obtain δ c iteratively. Using hybrid EAC and time-domain techniques, such as Single-Machine Equivalent (SIME), δ c can be obtained iteratively by examining different clearing times. SIME uses a generalized OMIB model and infers its parameters from the multi-machine temporal data. These data can be obtained either from time-domain transient stability simulations (Preventive SIME) or from real-time measurements (Emergency SIME) [1].
However accurate, the iterative approach may not satisfy the computation speed requirements for fast TSA. Therefore, some approximation in OMIB equations is proposed which enable the direct and rapid estimation of CCA.

Approximations for Rapid Estimation of CCA
All OMIB-based transient stability analysis methods rely on the observation that the loss of synchronism involves the irrevocable separation of generators into two groups: critical generators and non-critical generators [1]. Among the general class, the ones which are able to provide a direct estimation of the CCA are all based on some assumptions and approximations. As the first and common assumption, these approaches model the synchronous generators as a constant voltage behind the transient reactance, as shown in Figure 7.
To discuss the further assumptions, let us consider that for each critical generator δ i = δ cr + ξ i , and for each non-critical generator δ j = δ nc + ξ j , where ξ i and ξ j show the angular offset of each generator in set CR and NC from their respective PCOA δ cr and δ nc . As described in Appendix A, by considering the classical model of synchronous generators, the equation for OMIB electrical power becomes: where g ij = E i E j G ij and b ij = E i E j B ij , E i is the constant voltage behind direct axis transient reactance in the classical model of generator i, and B ij and G ij are the real and imaginary parts of the element of row i and column j of the admittance matrix reduced to synchronous generator internal nodes. The first approximate technique for rapid estimation of CCA assumes that for all generators ξ k and ξ l are zero, which means that the δ i of each generator is equal to its corresponding PCOA [6,7,9]. In this paper, this approach is referred to as Zero Offset OMIB (ZOOMIB). As described in Appendix A, with this assumption, the post-fault electrical power in Equation (13) becomes purely sinusoidal and is no longer dependent on δ e : where, as discussed in Appendix A, P C , P max , and v are constant values for the ZOOMIB model. As a second approximation approach, instead of using PCOA for all generators, it is possible to assume that the angle offsets are constant and equal to their pre-fault values. In this paper, this approach is referred to as Constant Offset OMIB (COOMIB). As described in Appendix A, similar to considering zero offsets, also with this assumption, the post-fault electrical power in Equation (13) becomes purely sinusoidal, independent of δ e and in the form of Equation (14).
ZOOMIB and COOMIB are both time-invariant models which freeze the generator angle offsets at the fault inception instant. However, the fact is that the angle offsets are not constant. In response to this limitation, a piecewise time-variant method is proposed in [10,11]. This approach, usually referred to as Dynamic OMIB (DOMIB), first makes an initial estimation of δ c and δ u . It then specifies some points between δ 0 and δ u , and simulates the generator angle trajectories by individual Taylor series described in Appendix C.2. Having the generator angles at each point, DOMIB updates the angle offsets and it considers an updated curve with constant offsets between the points. As shown in Figure 8, this approach leads to a piecewise sinusoidal δ − P curve. The most interesting point about the time-invariant and piecewise time-variant approaches is that the assumption of zero or constant ξ i values makes the post-fault P e of the OMIB model independent of δ e . Therefore, the variations of post-fault P e with δ can be defined uniquely. The merit is that EAC can be applied on the unique δ − P curves to obtain the critical clearing angle. Table 1 compares the assumptions and approximations of different methods based on the OMIB equivalent. With the sinusoidal P e in Equation (14), the area between P m and P e among an initial angle δ i and a final angle δ f can be obtained as follows:

Definition and Conditions of OMIB Stability
This section presents definitions and conditions for OMIB stability, that are based on three main assumptions: (1) generators' separation to critical and non-critical groups, (2) independence of the post-fault OMIB equivalent electrical power of the clearing time, and (3) for any δ there is only one possible value for P e and P m . These definitions and conditions provide a basis for defining a general function for CCA calculation.
Definition 1 (first-swing instability). An OMIB is first-swing unstable if and only if after fault inception its angle increases continuously with time.
Definition 2 (conditions for critical clearing angle for first-swing instability). Let δ be the angle of the OMIB equivalent of a power system. After a fault, the critical clearing angle δ c is the angle that satisfies the following conditions: where A D and A P are the area between P m and during or post-fault P e , and ε δ denotes a very small positive angle increment.
The condition in Equation (16) checks the attainability of δ c . The condition in Equation (17) states that by clearing the fault at δ c + ε δ , dδ/dt will remain positive for any δ above δ c , i.e., continuously increasing angle. The condition of Equation (18) states that for clearing angles less than or equal to δ c , the system remains stable. Definition 3 (backward-swing instability). An OMIB is backward-swing unstable if and only if after fault inception its angle decreases continuously with time.
Definition 4 (conditions for critical clearing angle for backward-swing instability). Let δ be the angle of the OMIB equivalent of a power system. After a fault, the critical clearing angle δ c is the angle that satisfies the following conditions: In this case, the direction of angle variations is backward. The condition in Equation (19) checks the attainability of δ c . The condition in Equation (20) states that by clearing the fault at δ c − ε δ , dδ/dt will remain negative for any δ below δ c , i.e., continuously decreasing angle. The condition in Equation (21) states that for clearing angles greater than or equal to δ c the system remains stable.
Theorem 1 (equivalency of first-swing instability and negative backward-swing instability). Let OMIB denote a one machine infinite node equivalent of a power system and let OMIB − denote an OMIB at which the signs of electrical power, mechanical power, and angles are negated. The OMIB model of the power system is backward-swing unstable, if and only if the OMIB − model is first-swing unstable. The critical clearing angle for the OMIB is equal to the negative of the equivalent OMIB − critical clearing angle.

Proof of Theorem.
Considering that A(δ a , δ b ) = −A(δ b , δ a ), by negating the signs of electrical power, mechanical power and angles in the conditions of Definition 2 for the critical clearing angle of a first-swing unstable OMIB we get: These equations represent the conditions for the critical clearing angle of a first-swing unstable OMIB − . Considering δ min = −δ max , the above conditions are equivalent to the conditions of Definition 4 for the critical clearing angle of a backward-swing unstable OMIB. Therefore, the critical clearing angle for the first-swing stability of OMIB − is equal to the negative of the critical clearing angle for the backward-swing stability of OMIB.

The EEAC Functions
In previous sections, we discussed the concept of EEAC, OMIB general equations, and approximations that enable one to have uniquely defined P e and P m curves for all δ. We also presented general definitions and conditions for OMIB stability.
Two basic functions are required for OMIB stability evaluation. The first function, as shown in Figure 9, takes the CC, NC, the reduced system admittance matrix, and the synchronous generators' data to calculate the OMIB model based on the specified type (ZOOMIB, COOMIB, or DOMIB). The output of this function is the OMIB inertia coefficient, P m and P e defined by P c , P max , and v as in Equation (14). For DOMIB, the constants defining P e will be different for each interval. The function should be called to form the pre-fault, during-fault, and post-fault OMIB equivalents when required. A pseudocode is presented in Appendix B, Algorithm A1, which details the process of the OMIB function.
-CC: set of names of critical machines -NC: set of names of non-critical machines -type of OMIB (ZOOMIB, COOMIB or DOMIB) -reduced system admittance matrix -synchronous generators data -synchronous generators angles for DOMIB

OMIB Inputs Outputs
-OMIB electrical power for each interval -OMIB mechanical power -OMIB inertia coefficient The second main function is a function using the conditions presented in the previous section to find the CCA of the developed OMIB models. As shown in Figure 10, this function takes the during-fault and pre-fault OMIB models as inputs. It also requires two input parameters: the angle step size ∆δ and the maximum integration limit δ max . The idea is to start at the OMIB initial angle δ 0 and to increase the fault elimination angle by an angle increment ∆δ to find δ c as the smallest fault elimination angle after which for any δ m ≤ δ max the area between P e and P m is positive.
-OMIB during-fault electrical power for each interval -OMIB post-fault electrical power for each interval -OMIB mechanical power -OMIB inertia constant

CCA Inputs Outputs
Parameters -critical clearing angle -return angle -direction of angular deviations (`forward swing' or `backward swing') -the type of the case detected (`always stable', `always unstable' or `potentially stable') : angle step size : maximum integration limit Figure 10. Schematic representation of CCA function to find the CCA of the equivalent OMIB model.
There are some exceptional cases that should be considered to avoid unreasonable results. The first case might happen when the maximum post-fault electrical power is less than the OMIB mechanical power P m . In such a case, the area between P m and P e is always positive and the system is unstable. Even if the maximum post-fault electrical power is more than P m , there might be other situations where the system is always unstable. As shown in Figure 11a, the area between P m and P e might be such that even for δ c = δ 0 it is always positive. Another exceptional case might happen for less severe disturbances, where the maximum during-fault electrical power is much more than P m . In such cases, shown in Figure 11b, the system will remain stable even without removing the fault. The function should be able to handle such cases and also backward-swing instability. It outputs the direction of angular deviations, 'first-swing' or 'backward-swing'; the type of the case detected, 'always stable case', 'always unstable case', or 'potentially stable case'; and the critical clearing angle. A pseudocode is presented in Appendix B, Algorithm A2, which details the process of CCA function.

Critical Machines Identification and Critical Cluster Formation
In a multi-machine power system, transient stability phenomena are governed by the critical machines, i.e., the set of machines responsible for the loss of synchronism following a large disturbance. Up to now, we have assumed that the critical cluster CC and the non-critical cluster NC are known. However, identification of the CC is one of the first steps of the EEAC algorithm and a prerequisite of OMIB equivalent model formation. Different OMIB equivalents can be formed for different possible sets of CC and NC. The true sets of CC and NC will be the ones with the smallest CCT. The reasoning behind this is that adding any critical machine of the true CC to the true NC, or adding any non-critical machine of the true NC to the true CC, will lead to slower OMIB dynamics, i.e., higher CCT. For a power system with n generators, the true CC can be identified by examining all possible combinations of n generators, i.e., 2 n − 1 candidates to find the ones with the smallest CCT. This exhaustive process would, however, be computationally demanding. The other solution is to find of a limited list of candidate critical generators in a CMI process. Then, in a CCF process, different pairs of CC and NC can be formed. An OMIB equivalent should be formed for each pair. The OMIB with the smallest CCT corresponds to the true pair of CC and NC. The next subsections present different methods for CCI and CCF.

Critical Machines Identification
All the techniques proposed for CMI are designed to provide a ranked list of critical machines to limit the number of possible combinations. Some are based on indices which rank the list of machines based on a criterion calculated for the fault inception time. Some others are based on a pre-estimation of CCT, to obtain the generators' t − δ trajectories, and to rank them based on their estimated rotor angles at an appropriate time after fault inception.

Acceleration Criterion
In the earlier stages of the EEAC development, the first approach for CMI was based on the initial accelerations the generators acquire at the disturbance inception [6,7,9]. According to this so-called "accelerations criterion", generators likely to be critical were considered to be those with the largest initial accelerations. For a given contingency, this approach first ranks the generators in a decreasing order of their initial accelerations calculated using Equation (A32) immediately following the fault inception.
Despite the encouraging results of this approach, the studies revealed difficulties of two types [9]. First, cm needs to be limited to avoid computational intractability. This may lead to unacceptable results for stability cases involving several critical generators. Second, it may happen that some generators not appearing at the top of the initial acceleration list experience considerable variations in their rotor angles after clearing the fault and eventually become unstable. In such cases, the initial acceleration criterion is not valid.

Composite Criterion
To improve the acceleration criterion, the "composite criterion" is proposed in [9]. The "composite criterion" relies on the initial accelerations together with the generators' pre-fault electrical distance to the fault to better identify the critical generators. It also considers the post-fault electrical distance of the generators to the fault to obtain a sense of the post-fault network.
To define it in the form of a criterion to rank the generators, for each generator k we can write: where γ k | t=0 + is generator k initial acceleration, and dist pre k and dist post k denote the preand post-fault electrical distances to the fault bus f , calculated as follows: where z ij is the magnitude of the element of row i and column j of non-reduced system impedance matrixẐ. The composite criterion was shown to perform better than the acceleration criterion in ranking the generators [9]. Nevertheless, it requires inversion of the bus admittance matrix to findẐ and to calculate the electrical distances. Moreover, calculation of electrical distance would be problematic when network splitting happens after fault clearance.

Trajectory Criterion
The trajectory criterion is proposed in [10,11] to rank the generators in order of their criticality. It is conjectured that the degree of criticality of a given generator is directly proportional to the magnitude of its rotor angle observed at an appropriate instant of time, in its evolution along an appropriate trajectory. The appropriate trajectory is a nearcritically cleared one, i.e., cleared at a time nearly above the actual CCT, and the appropriate observation time is defined as the time to reach the unstable equilibrium point of the OMIB equivalent of the power system.
The estimation of the generators' appropriate trajectory can be obtained by numerical integration. In [10,11], however, the Taylor series is employed as a quick substitute. Having an initial estimation of the CCT and the observation time, the trajectories can be obtained using individual Taylor series detailed in the Appendix C.

Critical Cluster Formation
The CMI gives a ranked list of critical machines. The aim of CCF is to form different combinations of CC and NC. The combinations should be later evaluated to identify the true combination of the clusters. Different techniques are proposed to form the clusters. A simple approach is to consider all possible combinations of critical machines as possible CCs [6,8,9]. A more efficient technique is presented in [10]. This technique selects cm candidate CCs composed of the first from the top, the first two from the top, . . ., up to all cm machines in the CC set.

CMI and CCF Functions
Three methods are discussed for CMI. As shown in Figure 12, for the acceleration criterion, the CMI function inputs the synchronous generators' data, their initial angle, and the during-fault system admittance matrix to provide a ranked list of generators based on calculated initial accelerations. For the composite criterion, the function also needs the pre-fault and post-fault system admittance matrix and the index of the fault bus (for line faults a virtual node should be added at the fault location) in the matrices to calculate the distances in Equation (26). For trajectory criterion, the reduced post-fault system admittance matrix, the fault elimination time, and the observation time are required. The individual Taylor series is employed to obtain the generators' angle trajectory and to rank them based on their angles at the observation time. In this paper, the observation time is defined as the time to reach the OMIB return angle. After ranking the generators with any of the criteria, the generators which are close to the top generator based on a predefined threshold are selected as critical ones and are outputted as a ranked list. A pseudocode is presented in Appendix B, Algorithm A6, which details the process of the CMI function.
The CCF function, as shown in Figure 13, receives the ranked list of the critical generators and forms different candidate pairs of CC and NC. A pseudocode is presented in Appendix B, Algorithm A7, presenting one simple method, among others, for CCF.

CMI Inputs Outputs
Parameters -ranked set of names of synchronous machines identified as critical CMI threshold : system base frequency -synchronous generators data and initial angle -type of CMI (acceleration, composite or trajectory criterion) -reduced during-fault system admittance matrix -reduced post-fault system admittance matrix (for trajectory criterion) -observation time and fault elimination time (for trajectory criterion) -pre-fault and post-fault system admittance matrix (for composite criterion) -index of the faulted bus in system admittance matrix (for composite criterion)

Integration
The EAC, despite all the information it provides, cannot directly give an indication of CCT which is of interest in transient stability studies. The CCT may be assessed by integrating the dynamics of the OMIB up to the point where it reaches CCA. In principle any numerical integration algorithm can be used. In [5][6][7][8][9][10], however, the Taylor series is employed as a handy and quick substitute for numerical integration. In the context of EAC-based methods, the Taylor series expansion can be applied to the OMIB equivalent of a power system, or to an individual generator to obtain its rotor angle evolution with time. The equations and the process of the Taylor series is presented in Appendix C.
As shown in Figures 14 and 15, despite the integration techniques employed, two main functions are required. The angle-to-time function inputs the OMIB equivalent model, the initial angle, the initial angular speed, and a desired angle (e.g., the CCA). It integrates the OMIB equations interval by interval to find the time and angular speed associated to the desired angle. The pseudocodes presented in Appendix B, Algorithms A8 and A9, present the details of this process with the Taylor series. As discussed in [9], for large departures of the desired angle from the initial angle, the Taylor series estimation may be inaccurate or may fail to give a result.
The Trajectory function, inputs the generators' data, their initial angle, during-fault and post-fault reduced system admittance matrices, the fault elimination time, and the desired final time of an individual generator's angle trajectory. It also receives the desired number of during-fault and post-fault intervals as input parameters. The function specifies some time instants based on the number of intervals and the time spans. It then calculates each individual generator's angle and angular speed interval by interval to reach the desired final time. The outputs will be the time instants, and the generator's angle and angular speed at each time. The pseudocodes presented in Appendix B, Algorithms A10 and A11, detail this process with Taylor series.

angle-to-time Inputs Outputs
-desired time at the given desired angle -desired angular speed at the given desired angle -OMIB electrical power for each interval -OMIB mechanical power -OMIB inertia constant -initial angle -initial angular speed -desired angle

Combining Algorithms for a Full-Resolution Scheme
The functions presented in the previous sections can be combined in different ways to provide an estimation of the CCT. The main functions include CCF and CMI, which can be based on acceleration, composite or trajectory criteria; OMIB which can be of type ZOOMIB, COOMIB, or DOMIB; CCA to estimate the critical clearing angle of the OMIB; angle-to-time to find the time corresponding to an OMIB angle; and Trajectory to find the trajectory of generator angles in time. These functions provide an insight to rethink the schemes proposed in previous literature, and to think of new schemes for direct CCT estimation.
The first step for any TSA technique is the preparation of the synchronous generators and the network data, and the formation and reduction of admittance matrices for pre-fault, during-fault, and post-fault states. EEAC relies on the classical model of the power system. Synchronous generators will be modeled with the classical model, and the admittance matrices include the loads and the generators' direct axis transient reactance. Admittance matrix reduction can be done by Kron method considering that all the nodes have zero injection currents except the internal nodes of synchronous generators.
The first and the simplest EEAC scheme that was proposed in [5,6,8] is as presented in Figure 16. The types are mentioned within green parenthesis, while the variable parameters are shown within red parenthesis. The pseudocodes presented in Appendix B, Algorithm A12, present the details of the basic-eeac scheme.
The scheme starts by identifying the critical machines (Here, the CMI is done with acceleration criterion, but it can also be done with composite criterion.) and forming a set of CCs and a set of NCs. It then evaluates each pair of CC and NC. For each pair, it first forms the pre-fault, during-fault, and post-fault OMIB equivalents (Here, the OMIB models are derived with ZOOMIB assumptions, but they can also be derived with COOMIB assumptions.). Having the OMIB equivalents defined, the CCA function is applied to find the CCA of the pair under consideration. Then, the CCT is calculated as the time to reach the CCA. After repeating these steps for each pair, the true CC and true NC are identified as the ones with minimum CCT. The algorithm finally returns the CCT, the identified clusters, the CCA and the angular speed, and the observation time as the time to reach the return angle δ r from δ 0 and may later be used for subsequent calculations. When an estimation of the CCT is made with the basic-eeac, the estimation can be improved in many ways. One approach can be as shown in Figure 17. Having a first estimation of CCT and observation time, the Trajectory function can be used to estimate the individual generator's angles for d during-fault and p post-fault intervals within δ 0 to δ max . With these estimated angles, the OMIB function can be recalled to make a better estimation of the OMIB equivalent model with DOMIB model assumptions. Then, the CCA function can be applied to the updated OMIB model to estimate the CCA, and the angle-to-time function can be employed to calculate the refined CCT. The above refinement process just re-estimates the CCT and does not update the estimate of the CC and NC. The second refinement process shown in Figure 18 uses the calculated CCT to find the individual generator angle trajectories and to rank them based on their angles at the estimated observation time. It then runs the CCF function to form a set of CCs and a set of NCs, and evaluates each pair of CC and NC. For each pair, it first forms the pre-fault, during-fault, and post-fault DOMIB equivalents. The CCA and angle-to-time functions are then applied to estimate the CCT for each pair. Finally, the pair with the minimum CCT is identified as the true pair (Here, the OMIB models are derived with DOMIB assumptions, but they can also be derived with ZOOMIB or COOMIB assumptions.).  Figure 19 shows a more sophisticated scheme. This scheme was proposed in [10]. Similar to the second refinement scheme, this scheme relies on the 'trajectory' CCI. However, for each pair of the CC and NC, it first forms a ZOOMIB equivalent, applies the angleto-time function, and obtains the generator angle trajectories. It then forms a DOMIB equivalent model using the obtained trajectories, applies the CCA function on the model, and calculates the CCT corresponding to each obtained CCA. Finally, the pair with the smallest CCT is identified as true CC and NC. The pseudocodes presented in Appendix B, Algorithm A13, present the details of the third refinement scheme. The pseudocodes of the other refinement schemes are simplified versions of this scheme and are not presented. The interesting point about this functional point of view of the EEAC is that the output of the basic-eeac is identical to the inputs and outputs of all three refinement schemes discussed. Therefore, as shown in Figure 20, these schemes can be repeated after each other to achieve the desired accuracy. As shown, after running the basic-scheme or each of the refinement functions, it is possible to terminate the calculations and output the estimated CCT and CC. We define each path from the input to the outputs as a branch. As ZOOMIB and COOMIB equivalents, and acceleration and composite CCIs can be used interchangeably, there are four possible combinations and thus four variants for each branch in Figure 20. Each branch has a different computational time which is a function of the system scale, the short-circuit scenario, and the chosen values for parameters. The parameters of each branch can be optimized to achieve the best performance in terms of CCT estimation accuracy and the computational time.

Simulation Studies and Discussions
This section presents the results for the application of the EEAC method on two test systems. The first is the four-machine system discussed before. The second is the French EHV power system with more than 400 synchronous machines, 2900 transmission lines, and 8800 transformers. Applying the EEAC on the four-machine system helps to investigate its approximations in detail, while studies of the French network helps to evaluate its performance for a real-life, large-scale network. For both test systems, the EEAC results are compared against time-domain simulations.

Four-Machine System
The considered scenario for the first test system, shown in Figure 2, is a three-phase short-circuit fault at one of the transmission lines, which is cleared by opening the line circuit breakers. As shown in Figure 3, for this case study, there is a clear separation between the NC generators and CC generators with increasing angles. A correct estimation of CC and NC allows us to evaluate the effect of assumptions of OMIB equivalent models, and also the effect of considering the classical model for synchronous generators.
To evaluate the assumptions, Figure 21 represents the δ − P curves obtained with ZOOMIB-and DOMIB-equivalent models. The curves are compared against the δ − P curves obtained from the time-domain simulation results, with the classical model for synchronous machines, and with the detailed model with regulators, i.e., speed governor and automatic voltage regulator. For time-domain simulation results, the angles and powers are obtained using Equations (8), (10), and (11).
As can be seen, the OMIB δ − P curve obtained with ZOOMIB assumptions is close to simulation results with the classical model. However, the modified piecewise curves obtained with DOMIB assumptions are much closer. As all the approximate OMIB equivalent models are based on the classical model of synchronous generators, the results they provide do not necessarily match the result obtained with the detailed model with regulators. However, the estimations of CCA are still close.
To better highlight the differences in estimations, Table 2 compares the CCT and CCA values obtained with different OMIB model assumptions against the values obtained from time-domain simulations. The Taylor series is employed to estimate the CCTs with the OMIB equivalents. In comparison with the results obtained with the classical model, the results obtained with ZOOMIB and COOMIB assumptions are acceptable. There is a clear improvement in the results obtained with DOMIB assumptions, and they are close to the classical model results. By considering a detailed model in simulations, the estimation errors increase; however, the DOMID still performs better.

French Network
This section briefly discusses the results of several fault scenarios considered in the French network. The considered scenarios are three-phase short-circuit faults on transmission lines, on bus-bars or on transformers. The French network covers several voltage levels. The faults are applied at different locations of the network, some close to large power plants on 400 kV, others on 225 kV portions of the network. The model considered for synchronous machines is a detailed model with regulators.
As discussed in Section 5, the functions presented in the previous sections can be combined in different ways to provide an estimation of the CCT. For each fault scenario, four of the possible schemes are examined, the basic scheme shown in Figure 16, and the basic scheme followed by the refinement schemes shown in Figures 17-19. The pseudocodes of these schemes are presented in Appendix B. For the basic scheme, the acceleration criterion is used for CMI. The scheme is evaluated with ZOOMIB and COOMIB assumptions.
The results are obtained by considering certain default values for parameters. However, the parameters can be optimized to find more accurate results. The considered values for angle step size ∆δ and maximum integration limit δ max are 0.1 and 360 degrees, respectively. Five intervals are considered for each of the during-fault and post-fault periods, and 50% is considered as the threshold for CMI.
The error in CCT calculation is calculated as follows: where CCT e and CCT a are the estimated CCT with EEAC and Taylor series, and the actual CCT obtained from the time-domain simulations, respectively. Table 3 compares the CCT values obtained with time-domain simulations against the basic scheme with ZOOMIB and COOMIB assumptions. Table 4 compares the CCT values obtained with time-domain simulations against the refined schemes for the same fault scenarios. As discussed in Section 2.5, there are some exceptional cases that should be considered to avoid unrealistic results. In both tables, 'stable' denotes the detection of an 'always stable' case, and 'unstable' denotes the detection of an 'always unstable' case (see Figure 11). For each of the considered case studies, the error percentage (Error) is calculated using Equation (27). Figures 22 and 23 compare the number of cases in each Error interval.
It is better for TSO to be more conservative and have an estimated CCT lower than the actual CCT, than having higher values. In other words, the TSO prefers to have a positive Error rather than a negative Error. As can be seen in the tables and the figures, with the basic scheme and ZOOMIB model assumptions, the errors are mainly positive, while with COOMIB assumptions fewer cases have a high negative Error. The refined schemes produce more accurate results for some fault scenarios. However, as discussed, the use of DOMIB assumptions modifies the estimation towards the classical model and the results are not necessarily close to the results obtained with the detailed model. On average, the basic scheme with ZOOMIB assumptions has better estimations than the other schemes. The first and third refined scheme decrease the maximum error, but all refined schemes have a larger error on average. Moreover, for more cases they detect an exception and do not provide a result.
The results show that a more complicated scheme does not necessarily provide more accurate results, while it might involve more computational time. For the case studies considered, the average computation time for basic schemes with ZOOMIB or COOMIB assumptions was around 30 s, while it was around 45, 70, and 100 s for the first, second, and third refined schemes, respectively. Moreover, the basic scheme does not identify any scenario as stable, thus less risk for the TSO.

Conclusions
These days, due to increased uncertainties in electrical power systems, there are an increasing number of transient stability scenarios that are of concern. Therefore, TSOs need fast TSA techniques to filter the scenarios and to identify the critical ones with lower CCT for detailed analysis. EEAC was proposed in the late 1980s as a promising and fast TSA method. However, despite the encouraging results and approaches presented through several papers, it was difficult to obtain a synthetic view of the key building blocks upon which the EEAC was built. This paper has revisited the EEAC from scratch. It has presented its very basic concept, the detailed equations, and the idea behind the approximations for fast TSA. New definitions and conditions have been defined for approximate models forward swing and backward swing stabilities. Based on these definitions and conditions, functions were developed for each EEAC building block, together with detailed pseudocodes. The idea was to propose a general full-resolution functional scheme that not only covers all the previous literature on the subject, but also introduces interesting possibilities for several new approaches.
Our studies show that the accuracy of the EEAC, though acceptable, depends on the selection of the sequence of functions and parameters. Once the optimal sequence of functions and parameters has been identified, the EEAC can serve as an effective tool for contingency filtering. It had reduced the time required for the analysis of a fault scenario in the French network from around 15 minutes for time-domain simulation to just a few seconds. However, further studies are required to design result quality indicators to tag cases where EEAC may not perform well. Moreover, the EEAC equations should be developed to consider non-synchronous generation units (e.g., windfarms) and HVDC links.

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. OMIB Electrical Power with the Classical Model
The classical model of a power system considers a constant-voltage-behind-transientreactance model for synchronous generators. In this model of the system, by dividing the network nodes to n synchronous generator internal nodes and r remaining nodes, the relationship between the bus voltages, nodal current injections, and the network admittance matrix is given by: whereẼ n denotes the synchronous generators' internal voltage behind their transient reactance,Ĩ n is the generators' current, andṼ r denotes the voltages of the remaining network nodes.Ŷ is the network admittance matrix which includes the load impedances and generator transient reactances. This matrix can be partitioned as follows: To obtain the electrical power, first we find the reduced admittance matrices by eliminating all the nodes except for the internal nodes of the synchronous generators. The reduction can be achieved through matrix operations considering that all the nodes have zero injection currents except for the source nodes. By eliminatingṼ r we have: In a power system with n synchronous generators, to calculate the output electrical power of each one we can write: whereẼ i = E i ∠δ i is the voltage behind direct axis transient reactance of the synchronous generator i, andŷ ij = y ij ∠θ ij is the element of row i and column j of the reduced admittance matrix. Therefore, for each generator we have: where G ij and B ij are conductance and susceptance parts of the admittance element of row i and column j. Expanding Equation (A4) and separating the real and imaginary parts we get: For each generator k of set CR we can rewrite Equation (A5) in the following form: We can consider that for each critical generator δ i = δ cr + ξ i , and for each non-critical generator δ j = δ nc + ξ j , where ξ i and ξ j show the angular deviation of each generator in sets CR and NC from their respective PSOA δ cr and δ nc . Therefore, for each generator k of set CR we can write: Similarly, for each generator l of set NC we can write: Considering δ i = δ cr + ξ k and δ j = δ nc + ξ j , we have: Substituting Equations (A7) and (A9) in Equation (11), we have: To have a simpler form of equation, by considering g ij = E i E j G ij and b ij = E i E j B ij , we have: Considering that M nc +M cr M T = 1, Equation (A11) can be written as follows: In this section, we simplify the P e for the OMIB model by considering zero rotor angle offsets. The assumptions are: With these assumptions, we can simplify Equation (A12) to the following form: On the other hand, in general we have: The equation for P e becomes: where: where P C , P max , and v are dependent on M nc , M cr , and g ij , i.e., constants.

Appendix A.2. Considering Constant Rotor Angle Offsets with Respect to PCOA
In this section, we simplify the P e for the OMIB model by assuming that ∀i ∈ CC, ∀j ∈ NC, ξ i , and ξ j are not necessarily zero, but that they are constant with respect to δ. With this assumption, we can simplify Equation (A12) to the following form: where: By separating sine and cosine terms we have: The equation for P e becomes: where: where P C , P max , and v are dependent on M nc , M cr , g ij , b ij , and ξ i , i.e., constants.

Appendix B. Pseudocodes
This appendix presents the pseudocodes of all the algorithms discussed. Algorithm A1 pseudocode details the function to compute the OMIB model. The function is designed for the DOMIB model, but it can be applied for COOMIB or ZOOMIB, which are similar to DOMIB but with only one interval. It inputs a data class of s generators including their name, inertia constant, internal voltage, mechanical power, and their initial and final angles for n intervals. It also requires the system admittance matrix to be reduced to generators' internal nodes, and the sets of critical and non-critical generators. The function should be recalled for each pre-fault, during-fault, and post-fault periods. For each period, the pseudocode first estimates OMIB P m and M, which are constant values. Then, for each time interval in the considered period, it estimates the OMIB angle at the beginning and end of the interval, and the terms defining its electrical power in that interval. for k = 1 : s do: end for 9: end for ∆δ. With two loops, for any fault elimination angle angle δ e , the algorithm searches for a return angle δ r as a δ m at which the sum of the areas become negative. If a δ r is found, the algorithm increases the δ e . The search continues until finding δ c , as the first δ e after which for any δ m ≤ δ max the sum of the areas is always positive. A 'potentially stable case' refers to a case where the system stability can be maintained by removing the fault at an angle below the identified δ c . The algorithm can also identify an 'always stable case' at which for δ e = δ 0 system is unstable, and an 'always stable case' at which for δ e = δ max system is stable. Algorithms A3-A5 serve Algorithm A2 as auxiliary functions. dflag: presents the direction of angular deviations, 'first-swing' or 'backward-swing': str 2.

Algorithm A2
tflag: indicates the type of the case detected, 'always stable case', 'always unstable case' or 'potentially stable case': str 3.
δ : for a 'potentially stable case' δ gives the critical clearing angle, for an 'always stable case' it gives δ max , and for an 'always unstable case' it gives δ 0 : float 4. δ r : for a 'potentially stable case' δ gives the return angle, for an 'always stable case' it gives δ max , and for an 'always unstable case' it gives δ 0 : float Algorithm A2 Cont. Parameter 1.
For trajectory criterion, the algorithm uses the Trajectory function to estimate the angles of individual generators at the observation time as the criterion. For the other two criteria, the algorithm first calculates generators' initial acceleration. If the criterion is acceleration, the initial accelerations are taken as the criterion. Otherwise, the pre-fault and post-fault distances to fault are estimated to calculate the composite criterion. Finally, the generators for which the calculated criterion is close to that of the top generator based on a predefined threshold are selected as critical generators, and a sorted list of them is outputted. type: type of CMI technique, 'Acc' for acceleration criterion, 'Comp' for composite criterion and 'Traj' for trajectory criterion: str 3.Ŷ red dur : reduced during-fault system admittance matrix: matrix of complex numbers 4.Ŷ red post : reduced post-fault system admittance matrix: matrix of complex numbers 5.Ŷ pre : pre-fault system admittance matrix: matrix of complex numbers 6.Ŷ post : post-fault system admittance matrix: matrix of complex numbers 7.

Algorithm A6
f : index of the faulted bus in system admittance matrix, set to 1 by default: int 8. end for 21: for j = 1 : s do: 22: if criterion[j] > threshold · max(criterion) then 23: append S[j].name to CM 24: end for 25: sort CM in decreasing order of criterion 26: return CM Algorithm A7 pseudocode details the function to form different candidate CCs and NCs. This pseudocode represents a simple method among others. For cm critical machines, the algorithm selects cm candidate CCs composed of the first from the top, the first two from the top, ..., up to all cm machines in the CC set. It outputs the set of candidate CCs and the set of candidate NCs.
S.name: set of names of synchronous generators: set of str 2.
CM: ranked set of names of synchronous generators identified as critical: set of str Output 1.
SCC: sets of CC: sets of str 2.
SNC: sets of NC: sets of str

SNC[j]=S.name not in SCC[j]
5: end for 6: return SCC , SNC Algorithm A8 pseudocode details the function to estimate the time associated with a desired angle for the OMIB model. The function angle-to-time updates the time and OMIB angular speed interval by interval up to reaching the desired time associated to the desired angle. This function relies on Algorithm A9 pseudocode which presents the GTS function. For each interval, this function employs the Taylor series equations to find the time and angular speed associated to a desired angle using the interval initial values. δ 0 = sin −1 P O [1].P m −P O [1].P c P O [1].P max + P O [1].v δ 0 = sin −1 P O [1].P m −P O [1].P c P O [1].P max + P O [1].v δ 0 = sin −1 P O [1].P m −P O [1].P c P O [1].P max + P O [1].v To obtain the evolution of δ k and ω k with time, Equations (A30) and (A35) should be updated together. The obtained values at each time instant should be employed to initialize the Taylor series for the next time step.