Feature Extraction-Based Real-Time Transient Stability Analysis

With increasing amounts of asynchronous generation being deployed to meet system energy demands, many transmission corridors may become constrained by angular stability criteria rather than steady-state thermal limitations. In such a system, it is paramount to have the capability to rapidly evaluate the stability margins of the system, particularly in a threatened post-fault state. The use of single machine equivalents (SIME) has been shown to be a powerful and flexible hybrid stability analysis method which can be computed directly from measured PMU data. However, due to the nature of the system reduction employed by SIME, as well as the method of extrapolation to estimate stability margins, there are many cases where the swing of the system is not accurately modelled by the traditional methods until a significant amount of data has been collected, at which point it may be too late to respond to the threat. In this paper, we address some of the limitations imposed by the traditional methods of reduction and prediction. We propose a method where rather than identifying a single critical SIME model for stability prediction, a spectrum of SIME representations of the post-fault system is developed, yielding a more timely and accurate estimate of post-fault stability conditions, through observation and feature extraction.


Introduction
Climate change and greenhouse gas emissions have become recognized as problems of international significance in recent years. Though renewable energy sources are considered a potential solution to these problems, their integration into existing power grids brings a new set of challenges. The unpredictable and uncontrollable nature of these sources of generation will require the development of advanced analytics and distributed control algorithms in order to manage such a highly distributed generation base. As integration of renewable resources continues, though it assists in decarbonisation of the power industry, further challenges arise when the system is sourcing more power from renewables. Since fewer conventional power plants will remain in operation under this scenario, methods of establishing system security will need to be updated to cope with the loss of system inertia.
On the other hand, the deployment of Phasor Measurement Units (PMUs) on large scale in power systems creates the basis for near-real-time and accurate system state observations. With high sampling rate and GPS time synchronisation, PMUs enable a new stream of applications for stability analysis. The data points generated by PMUs can be used to estimate the stability of the network and even the margin.
Stability of power systems has constituted a concern for system operators and several analysis techniques have been developed, such as time-domain simulation methods [1,2]. These offer good accuracy, but they are computationally intensive as they require solutions of ordinary differential equations. They do not provide any information on the margin. Direct methods are faster but can lead to narrow margins [3].
Several strands of hybrid stability methods have been utilized to ensure accuracy and low computational burden. The one used in this paper is the single machine equivalent (SIME) [4][5][6][7][8]. Hybrid methods resemble direct methods by using energy functions, but the input data consists of measurements taken during and immediately after the fault. Determining the generators that are likely to go out of step during a fault is an important task in SIME methods [5,9,10,11], with approaches such as Preventive SIME and Emergency SIME being applied [6]. Several papers have focused on energy functions for power systems property description [12][13][14][15]. The methodology described in this paper can help evolve power systems towards smart grids through advancement of realtime monitoring and control algorithms [16][17][18][19][20][21][22][23]. This paper gives a brief overview of single machine equivalents and describes the time-domain evolution of a SIME model's kinetic energy trajectory, and then analyses the features of post-fault swings. The results presented here prove the usefulness of the method application for large power systems.

Single Machine Equivalents
The aim of SIME is to generate a simplified model of a power system in the post-fault state, employing a hybrid stability analysis approach where measured or simulated data of the swing of the generators serves as a set of inputs for generating the model. This technique creates a simplified power system model that is useful for extrapolation and prediction of system behaviour, and that can be updated with real-time data. The importance of SIME models in the proposed method warrants a brief treatment of the technique. The SIME model formulation begins by separating the generators into two distinct groups in the post-fault state, where C is the set of critical generators, NC is the set of noncritical generators, G k indicates the kth generator in the system, and N G is the total number of generators in the system. While a number of methods could be used to determine which generators belong to the critical set for a given fault, the most often applied approach is through observation of the post-fault angular separation of the generators [5]. Once the generators are split into these two sets, the inertia of the equivalent critical and non-critical machines can be calculated respectively as, and the SIME inertia is then calculated as This renders the SIME representation equivalent to a single machine infinite-bus (SMIB) system.
Similarly, the rotor centre-of-angle for the equivalent critical and non-critical machines is and the frequencies of the critical and non-critical machines as determined respectively as The rotor angle of the SIME model is then found as and the frequency as The mechanical input power to the SIME model is calculated from the generator data as and the electrical output power is Subsequently, the accelerating power of the SIME model is found to be P a;S t ð Þ ¼ P m;S t ð Þ−P e;S t ð Þ: ð10Þ For any given SIME model of a system in the post-fault state, the system can be approximated as where any S contains the set of data which describes a SIME system representation. In prior usage, the critical SMIB refers to the SIME model which first goes unstable, thus having the minimal stability margin of all possible SIME representations.

Kinetic Energy Measures
In prior research, SIME methods of estimating stability margins in the post-fault state have focused on fitting and extrapolation of the P a -δ curve of the SIME model [23]. Because the SIME representation is analogous to an SMIB structure, analysis of the P a -δ characteristics is an application of the equal area criterion, which is also equivalent to the Lyupanov criteria for a single machine system [22,23].
While forward-prediction of system stability using the SIME P a -δ curve is effective for cases where the curve is nearly sinusoidal, these cases only occur when a given fault drives the generators to separate into two nearly-coherent sets. For other faults where there is weaker coherency within the critical and non-critical sets, the SIME approximation's P a -δ curve is unlikely to have a sinusoidal shape, thus requiring a more sophisticated extrapolation technique.
In the proposed method, rather than observing the P a -δ curve to calculate a stability margin estimate, the time-domain evolution of a SIME model's kinetic energy trajectory forms the basis for the analysis and prediction of system stability.
Initially, the kinetic energy of a chosen SIME representation S is known to be Calculating the derivative of the SIME kinetic energy with respect to time yields From the swing equation, the relationship between change in frequency and accelerating power is and substitution of this relationship into the time-derivative of kinetic energy yields Finally, the time-derivative of the kinetic energy is inverse weighted by the SIME model's inertia, yielding For any SIME model, Y S (t) forms the characteristic feature which is utilized in the proposed method for model extraction and stability prediction.

Features of Post-Fault Swings
When assessing stability estimates of a system in the post-fault state, additional insight into the behaviour of the system may be obtained by observing not only a single SIME representation of the system, but also through the observation of the behaviour of several SIME representations. Some prior work exists that observes, in addition to the critical SIME model, the P a -δ characteristics of other candidate SIME models that are similar but not identical to the critical SIME model [6].
The approach taken in the proposed method involves the intentional simultaneous generation of a set of SIME models, which are meant to provide additional insight into the postfault behaviour of the system. The set of SIME representations chosen in the proposed method are deliberately selected to observe the system behaviour for several post-fault conditions. Figure 1 demonstrates the concept of the multiple SIME method as proposed in this paper. The spectrum of SIME representations selected for this study are formed by first observing the angular frequencies ω of the generators at the time of the fault clearing t fc . The generators are then sorted into an ordered set according to their angular frequency as sorted from greatest to least. Utilizing this ordered set L, the set of N G -1 SIME representations are formed, where for k = 1… N G -1, the critical and non-critical generator sets are defined as That is, for k = 1, the first generator in L is treated as belonging to the critical set, and the remaining N G -1 generators form the non-critical set. For k = 2, the first two generators in L are treated as belonging to the critical set, and the remaining N G -2 generators form the non-critical set. This is repeated for all k = 1 … N G -1.
With the critical and non-critical sets defined, the SIME equations can be applied to create an SMIB representation of the system based on the kth critical and non-critical generator sets, which are then used to form system properties of this representation as S k . Application of the SIME equations thus yields the spectrum of SIME representations S k , k = 1 … N G -1.
Each of the SIME system representations S k contains the set of data necessary for calculating the characteristic feature Y Sk (t), and thus N G -1separate Y trajectories are calculated from the various generator set groupings. An interesting pattern emerges in the behaviour of these postfault characteristics Y for the spectrum of SIME models S k . As Fig. 1 The concept of Multiple SIME representations seen in Fig. 2, while each of the N G -1 trajectories are unique, they tend to follow a post-fault group trend.
The mean value of the group trend is calculated as the inertially-weighted average, and if this mean trend is removed from the Y trajectories, the remaining mean-adjusted trajectories behave as shown in Fig. 3.
As seen by the post-fault behaviour of the mean-adjusted trajectories Y Sk (t) -〈 Y S (t) 〉, these tend to follow a damped oscillatory pattern prior to the system reaching instability and going out-of-step. The method proposed in this paper seeks to extract the features of the post-fault characteristic trajectories Y Sk (t), by calculating a second-order state space model which encapsulates the group trend as well as the damped oscillatory nature of these signals.

Feature Extraction Technique
In the event of a system fault, if PMUs are deployed within a system, then post-fault estimation of system stability may be achieved through analysis of the discrete-time data set obtained from measurement. Within this study, it is assumed that the measured data obtained from the PMUs is received in discrete time intervals, spaced by a time step of Δt. For each time step, a new set of measurements is obtained from the PMUs deployed at each generator, and this data is subsequently used to form the spectrum of SIME models presented in Section 4. Utilizing the newly obtained generator data for the current time step tc, the spectrum of SIME models S k , k = 1 … NG -1 is formed, and the corresponding Y Sk (t) are calculated, up to the most recently available time interval t c . The sequence of data from the Y Sk (t) are used to compute a least-squares fitting of a second-order state space model, which can then be used for extrapolation and post-fault stability estimation.
The procedure begins with the computation of a set of finite differences, with the first-order finite difference defined as, and the second-order finite difference defined as Computing the relevant finite differences for the most recently available data at time t c over the spectrum of computed SIME models S k yields a coefficient matrix, and the corresponding second-order differences are used to form a column vector as The least-squares solution to the linear problem yields a feature extraction vector f(t c ) which is used to form the second-order state space approximation. This approach can be extended to not only include samples from the most recently available data at time t c , but can be extended to incorporate further past data, increasing the dimension of the least-squares problem. This is achieved through vertical concatenation of previously calculated A and b, forming the expanded leastsquares problem Fig. 2 The calculated characteristic feature Y S (t) for fault number F14 which can be generalized as Note that, in order to ensure that the problem is not underdetermined, there is a constraint on the first time instance after fault clearing where this model can be derived, where given that m T is equal to the maximum order derivative used in the formation of b, plus the number of previous time samples m p used in the formation of A and b.
Calculating the feature extraction vector proceeds by computing where superscript † indicates the pseudo-inverse, though any suitable method of least-squares solution, such as QR decomposition or SVD decomposition, is acceptable.
The desired system behavioural model for approximating the time-domain evolution of the characteristic signals Y Sk (t) is a second-order state space model, thus being composed of four elements in a two-by-two matrix. This behavioural model, F(t c ), contains elements which describe the change in one state variable per unit time with respect to a given state variable. In this instance, the four elements of the model are where the subscript υ indicates the state variable corresponding to the current value of the signal Y Sk (t), and the subscript s indicates the state variable corresponding to the current slope of Y Sk (t). Thus, for example, δ sv would be the change in slope per unit time dependent on the current value of Y Sk (t).
The state space behavioural model F(t c ) is constructed as where the state space coefficients determining the change in value and slope of Y Sk (t) with respect to the current value are set to zero and unity respectively. The coefficients which are dependent on the slope, which form the second column of the state space model, are formed from the previously calculated values in the column vector f(t c ). This forms the state space behavioural model for the current time step, which can then be used for extrapolation and prediction, and which can also be recalculated for each new time step where new data is received.

Extrapolation and Prediction
Having extracted a suitable behavioural model from the characteristic features Y Sk (t) of the multiple SIME representations S k , the method proceeds by utilizing the behavioural model F(t c ) for extrapolation to predict the system trajectory, and subsequently generate an estimate of the stability margin η under the relevant fault conditions. A framework that permits trajectory sensitivity analysis was proposed in [20], but the method proposed in this paper uses parallel simulations of SIME models and feature extraction to speed up the transient analysis.
The first stage in the extrapolation procedure is identical to a stage that is also present in the ESIME method [6], where the goal is to determine the critical SIME model S*. At each time step of the analysis, an estimate of the composition of the critical machine is determined through observation of the rotor angle trajectories of the system generators. Specifically, a Taylor Series extrapolation of the generator rotor angles is calculated, and the system generators are classified into two sets: critical and non-critical, according to the maximum separation between the extrapolated rotor angles.
In this stage, the critical SIME model S* is a separate system representation which will be extrapolated to estimate system stability. While a set of SIME representations S k were utilized to extract the behavioural model F(t c ), the critical SIME model S* determined from the generator angular separation is not necessarily equivalent to one of those in the set S k .
With the critical SIME model S* chosen for the current time step t c , the behavioural model F(t c ) can then be used for extrapolation of the Y S* (t) trajectories. The vector of cur- and for each discrete time step Δt, a forward prediction of the values is calculated as yielding extrapolated estimates of Ŷ S* and Ŷ This forward extrapolation process is continued through repeated incrementation by Δt, until the estimated value of Ŷ S* crosses the zero threshold going from negative to positive. The time which this crossing is predicted to occur is defined as time t s , which also represents the stopping criteria of the extrapolation. From the extrapolation, a discrete time-series estimate of the critical tra-jectoryŶ S* is formed, starting from the most current received measurement time t c until the estimated zero-crossing time t s , This yields a time series of values over which a numerical integration can be performed to estimate the stability margin of the system during the currently analysed fault. An example of the calculated Y S* (t) for fault number F15 is shown in which can thus be seen as the estimation of the amount of accumulated kinetic energy which will be converted to potential energy during the system swing. From the equal area criterion, it is a well-known constraint for post-fault system stability that the accumulated kinetic energy during a fault condition must be sufficiently low, such that the excess kinetic energy can be fully absorbed as potential energy during the post-fault swing. Because the SIME representation yields an SMIB structure, application of the equal area criterion is used to estimate, at each time step, the margin of stability, which is recomputed at each time step by observing the difference between the remaining excess kinetic energy and the estimated stable deceleration area. The stability margin estimate η is thus a time-varying estimate of system stability, which is updated at each discrete time step when a new set of PMU data is collected. Figure 5 is the flowchart of the aforementioned steps for calculation of the stability margin estimate η.
The proposed approach for obtaining these stability estimates is thus an enhancement to the established ESIME method [6], and will be shown in the subsequent section to yield excellent results, matching the performance of the ESIME technique and enhancing the accuracy of the stability estimate in many cases.

Results
The results of testing the proposed method for post-fault system stability estimation are presented in this section, verifying the accuracy and validity of the method. The proposed approach is applied to a number of post-fault test cases for both the IEEE 9-bus and IEEE 39-bus test networks, and a set of performance metrics for the approach are collected. For comparison, an identical set of performance metrics are collected Fig. 6 Comparison of the calculated estimated stability margins for the same fault set, alternatively applying the ESIME method for stability estimation. This section compares the results of the proposed method to the predecessor ESIME method described in the literature, displaying that the proposed technique yields comparable results, while enhancing prediction accuracy.
For both methods, the periodic update of information from the PMUs located at each of the generator buses is assumed to happen at a rate of once per cycle. In this case, the newly received data permits a discrete update of the post-fault stability estimate every 16.6 ms. Every time a periodic update of PMU data is received, a stability estimation algorithm iteration is performed, and the time-varying behaviour of the stability estimate η is recorded.
Comparison of the performance of the proposed method against the ESIME method is observed through the computation of a set of metrics, derived directly from the time varying stability estimates. The first metric μ E is the mean error between all time-varying samples of η and the final converged value of the stability margin estimate η f . As seen in Fig. 6, both methods converge towards an identical η f as time advances and more information is received, but the progression for each method's η may vary.
The second metric N R is the number of samples, prior to the system going unstable, for which the stability margin estimate η is within a ± 10% range of η f . The third metric indicates the number of time samples N + for which the method yields a positive estimate for η, despite that all tested faults are unstable, and should yield a negative stability margin estimate.

9-Bus System
The proposed method was tested by simulating a number of fault conditions on the IEEE 9-Bus system as shown in Fig. 7, then collecting the power, frequency, and angular characteristics of the fault event from each generator bus once per cycle, yielding a data set identical to that which would be available if a PMU were deployed at each generator bus. Utilizing this discrete-time data set, the proposed method was applied to estimate the stability margin once per cycle with each new data update, as well as applying the ESIME method for stability analysis on the same data. A detailed list of created fault scenarios is provided in Table 1. As seen in Table 2, the results of the proposed method yield comparable results to those obtained via the use of ESIME. For this particular system configuration, the ESIME approach tends to have a slightly reduced average error for the stability margin estimate, enabled by the nearly sinusoidal post-fault power angle swings in this simplistic network. Because ESIME applies a quadratic approximation, it is particularly well suited for extrapolative predictions in this simple network. The proposed method obtains stability margin estimates which are nearly as accurate, but also avoids the significant false-positive stability margin estimate during Fault 2. As such, the proposed method need not be seen as a replacement of the ESIME method, but certainly provides stability estimation information which can enhance system insight.

39-Bus System
To verify the performance of the proposed method on a more realistic and larger network, IEEE 39 bus system as shown in Fig. 8 was simulated and a list of twenty different faults were created on generators and randomly selected busbars.
For a busbar fault, both self-clearing fault and the outage of the connecting lines were considered. A detailed list of created fault scenarios has been provided in Table 3. A comparison of results using proposed performance matrices for all the above mentioned fault scenarios are shown in Table 4.
The strengths of the proposed method are clearly observable for the 39-bus system test case, where for nearly every one of the twenty fault scenarios studied, the proposed method yields a more accurate estimate of the stability margin in the post-fault period. Furthermore, it provides an accurate estimate more rapidly, as evidenced by N R , giving a supervisory control system more time to react in case some corrective action is required to prevent an out-of-step condition. It can be up to 100% faster in some cases. On the other hand, the falsely estimated points N + happen for only two samples in the F5 scenario whereas this number is as high as 40 samples for four different fault scenarios.

Conclusions
The increasing share of renewable energy sources within the generation portfolio, along with continual increase in demand, has created conditions which may drive transmission corridors  towards their stability limits. In such systems, rapid evaluation of the stability margins is becoming more vital to provide adequate for corrective measures to be effective in maintaining system stability. On-going advances in continuous realtime monitoring and control of power system stability have provided the enabling infrastructure for online and rapid system stability evaluation. Utilizing the enhancements brought by the PMUs in online system monitoring, this paper presents a novel approach as a direct enhancement to SIME, for online and rapid assessment of system stability. The proposed method provided a derivation creating not only a SIME representation of a post-fault system, but rather a spectrum of SIME representations which are used simultaneously to extract and extrapolate stability features of the system.
Comparative results on a variety of fault conditions for the test networks shows that it is capable in many conditions of generating more accurate estimates of stability margins, as compared to previously utilized methods. In addition to the achieved higher accuracy, in many cases the proposed method also more rapidly achieves an accurate estimate of the post fault stability margin, allowing a supervisory control system more time to react where corrective action is required.
The proposed method underperforms in cases where the number of buses is small, but it is highly efficient for large systems. As future work, the method could be improved by using time series analysis based on the PMU data to estimate the parameters related to each mode of system response to faults.
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.