Model-Free MLE Estimation for Online Rotor Angle Stability Assessment with PMU Data

TL;DR

This paper introduces a real-time, model-free method for assessing rotor angle stability using PMU data by estimating the maximal Lyapunov exponent with a recursive least squares approach, enabling timely stability evaluation.

Contribution

It proposes a novel online stability assessment method that accurately estimates the MLE from real-time measurements without relying on system models, addressing fluctuation issues.

Findings

Effective stability assessment on test systems

Reliable detection of stability conditions

Timely evaluation using real-time data

Abstract

Recent research has demonstrated that the rotor angle stability can be assessed by identifying the sign of the system maximal Lyapunov exponent (MLE). A positive (negative) MLE implies unstable (stable) rotor angle dynamics. However, because the MLE may fluctuate between positive and negative values for a long time after a severe disturbance, it is difficult to determine the system stability when observing a positive or negative MLE without knowing its further fluctuation trend. In this paper, a new approach for online rotor angle stability assessment is proposed to address this problem. The MLE is estimated by a recursive least square (RLS) based method based on real-time rotor angle measurements, and two critical parameters, the Theiler window and the MLE estimation initial time step, are carefully chosen to make sure the calculated MLE curves present distinct features for different stability conditions. By using the proposed stability assessment criteria, the developed approach can provide timely and reliable assessment of the rotor angle stability. Extensive tests on the New-England 39-bus system and the Northeast Power Coordinating Council 140-bus system verify the effectiveness of the proposed approach.

Tables

Table 1. TABLE I: Summary of the Tests on the New-England 39-Bus System

$t_c/$s	Occurrence Times of the Swing Patterns
	I	II	III	IV	V	VI
1.08	N/A	N/A	4	43	6	3
1.16	N/A	6	10	36	3	1
1.24	4	13	9	27	2	1
1.32	38	12	2	4	N/A	N/A
Success Rate	100%	100%	100%	100%	100%	100%

Table 2. TABLE II: Summary of the Tests on the NPCC 140-Bus System

$t_c/$s	Occurrence Times of the Swing Patterns
	I	II	III	IV	V	VI
0.18	N/A	N/A	6	166	2	5
0.26	N/A	10	26	143	N/A	N/A
0.32	8	57	38	76	N/A	N/A
0.40	60	12	49	58	N/A	N/A
Success Rate	100%	100%	100%	100%	100%	100%

Equations

$${\lambda }_k\approx \frac{1}{k\Delta t}\log \left(\frac{d\left(m(n),n,k\right)}{d\left(m(n),n,0\right)}\right)$$

$$=\frac{1}{k\Delta t}\log \left(\frac{||{\mathbf{X}}_{m(n)+k}-{\mathbf{X}}_{n+k}||}{||{\mathbf{X}}_{m(n)}-{\mathbf{X}}_n||}\right),$$

$${\mathbf{X}}_n,{\mathbf{X}}_{m(n)},{\mathbf{X}}_{n+k},\text{and}{\mathbf{X}}_{m(n)+k}\in \text{Phase II},$$

$${\mathbf{\Theta }}_t={\left[{\theta }_t,{\theta }_{t-\tau },{\theta }_{t-2\tau },\cdots ,{\theta }_{t-2N\tau }\right]}^{\top },$$

$$d(m(n),n,i)=||{\mathbf{\Theta }}_{m(n)+i}-{\mathbf{\Theta }}_{n+i}||,i\in \{0,k\},$$

$$d(m(n),n,i)=|{\theta }_{m(n)+i}-{\theta }_{n+i}|,i\in \{0,k\}.$$

$$L(m(n)+i)$$

$$=\log \left(|{\theta }_{m(n)+i}-{\theta }_{n+i}|\right),i=0,1,\dots ,k.$$

$$L(m(n)+i)={\lambda }_k\cdot (m(n)+i)\Delta t+C_k+{\xi }_k,$$

$$i=0,1,\dots ,k,$$

$$\displaystyle\hat{\boldsymbol{E}}_{k}=\left[\begin{array}[]{c}\lambda_{k}\\ C_{k}\end{array}\right]=(\boldsymbol{X}_{k}^{\top}\boldsymbol{X}_{k})^{-1}\boldsymbol{X}_{k}^{\top}\boldsymbol{Y}_{k},$$

$$\displaystyle\boldsymbol{X}_{k}\!=\!\left[\!\begin{array}[]{c c}m(n)\Delta t&1\\ (m(n)+1)\Delta t&1\\ \cdots&\cdots\\ (m(n)+k)\Delta t&1\end{array}\!\right],$$

$$\displaystyle\boldsymbol{Y}_{k}\!=\!\left[\!\begin{array}[]{c}L(m(n))\\ L(m(n)+1)\\ \cdots\\ L(m(n)+k)\\ \end{array}\!\right].$$

$${\hat{\mathbf{E}}}_{k+1}={\hat{\mathbf{E}}}_k+{\mathbf{G}}_{k+1}\left[y_{k+1}-{\mathbf{x}}_{k+1}^{\top }{\hat{\mathbf{E}}}_k\right],$$

$${\mathbf{G}}_{k+1}=\frac{{\mathbf{P}}_k{\mathbf{x}}_{k+1}}{1+{\mathbf{x}}_{k+1}^{\top }{\mathbf{P}}_k{\mathbf{x}}_{k+1}},$$

$${\mathbf{P}}_{k+1}={\mathbf{P}}_k-{\mathbf{G}}_{k+1}{\mathbf{x}}_{k+1}^{\top }{\mathbf{P}}_k,$$

$$d_j$$

$$\approx \left|\left({\theta }_{0+j}+\sum _{t=j+1}^{j+w}{v}_t\cdot \Delta t\right)-{\theta }_{0+j}\right|=\left|\sum _{t=j+1}^{j+w}{v}_t\cdot \Delta t\right|,$$

Full text

I Introduction

Transient rotor angle stability refers to the ability of synchronous generators of an interconnected power system to remain in synchronism after a severe disturbance [1]. With the development of synchrophasor technologies, utilities are now able to track rotor angle deviations and take actions to respond to emergency events. However, since the dynamics of power systems are complex, online rotor angle stability assessment is still very challenging [2, 3].

In [4], an adaptive out-of-step relay is proposed for the Florida-Georgia system. The equal area criterion is applied to change the settings of the protection system based on phasor measurement unit (PMU) measurements. In [5], the dynamics of the power transfer paths are monitored based on the energy functions of the two-machine equivalent system, and the PMU data are used to identify the parameters of the energy functions. In [6], PMU measurements are used as inputs for estimating the differential/algebraic equation model to predict the post-fault dynamics. In [7], an online dynamic security assessment scheme is proposed based on self-adaptive decision trees, where the PMU data are used for online identification of the system critical attributes. In [8], the rotor angle stability is estimated by using artificial neural networks and the measured voltage and current phasors are used as inputs of the offline trained estimation model. In [9], a systematic scheme for building fuzzy rule-based classifiers for fast stability assessment is proposed. By testing on a large and highly diversified database, it is demonstrated that the analysis of post-fault short-term PMU data can extract useful features satisfying the requirements of stability assessment.

Lyapunov exponents (LEs) those characterize the separation rate of infinitesimally close trajectories are important indices for quantifying the stability of dynamical systems. If the system’s maximal Lyapunov exponent (MLE) is positive, the system is unstable, and vice versa. LEs are first applied to power system stability analysis in [10], in which it is verified that LEs can predict the out-of-step conditions of power systems. In [2], a model-based MLE method is proposed for online prediction of the rotor angle stability with PMU measurements. The work builds solid analytical foundations for the LE-based rotor angle stability assessment. In [11], the LEs are calculated with dynamic component and network models to identify the coherent groups of generators.

Although the model-based MLE estimation approaches have made significant progresses on online rotor angle stability assessment, they are usually computationally expensive especially when applied to large power systems. Therefore, two model-free MLE estimation approaches have been proposed for transient voltage stability assessment [12] and rotor angle stability assessment [3], for which the MLEs can be estimated by only using PMU measurements.

The model-free LE-based stability assessment approaches are attractive, because they can eliminate model errors and simplify the calculation. However, when applying these approaches, a time window has to be pre-specified for the MLE observation. The window size is crucial for obtaining reliable and timely assessment results, i.e., too small window size will lead to unreliable assessment results while too large window size will lead to accurate but untimely results. The window size is difficult to be determined in advance, because the estimated MLEs may fluctuate between positive and negative values for quite a long time after disturbances and the window size should change with different fault scenarios.

In this paper, a LE-based model-free rotor angle stability assessment approach is proposed. The MLEs are estimated by a recursive least square (RLS) based method based on real-time rotor angle measurements. By properly choosing two critical parameters according to the characteristics of the relative rotor angles of the selected generator pairs, the calculated MLE curves will present distinct features for different stability conditions, based on which the stability criteria are correspondingly designed to capture the MLE features and perform online rotor angle stability assessment. Compared with the existing approaches, the proposed approach does not need a pre-specified time window to identify the sign of MLEs. Instead, the proposed approach can always make a reliable and timely assessment as soon as the crucial features are observed from the estimated MLE curves.

The remainder of this paper is organized as follows. Section II introduces the theoretical basis of model-free MLE estimation. Section III proposes a rotor angle stability assessment approach, and discusses the parameter selection principles and stability criteria. In Section IV, simulation results on the New-England 39-bus system and the Northeast Power Coordinating Council (NPCC) 140-bus system are presented to validate the effectiveness of the proposed approach. Finally, the conclusions are drawn in Section V.

II MLE Estimation From Time Series

LEs can reflect the exponential divergence or convergence of neighboring trajectories in the state space of a dynamic system [13]. An $N$-dimensional dynamic system has $N$ LEs and the largest one is defined as the MLE of the system. MLE is a useful indicator of system stability: A positive MLE indicates unstable system dynamics while a negative MLE indicates asymptotically stable dynamics [2]. The MLE can be estimated by using Jacobian matrix based (model-based) methods [13, 2] or direct model-free methods [13, 12, 3]. Compared with the Jacobin matrix based methods, direct methods are more suitable for online stability assessment mainly because they do not need the repeated computing of the Jacobian matrix or even the dynamic model of the system.

According to Oseledec’s multiplicative ergodic theorem [13, 14, 15], for a reference point $\boldsymbol{X}_{0}$ and its neighboring point $\boldsymbol{X}_{m(0)}$ chosen from the state space of a nonlinear dynamic system, the distance between the trajectories emerged from $\boldsymbol{X}_{0}$ and $\boldsymbol{X}_{m(0)}$, i.e., the original trajectory and the neighboring trajectory, will have three different growth phases as illustrated in Fig. 1. In Phase I, the difference vector between the states of the trajectories gradually converges towards the most expanding direction, and the distance between the trajectories will exhibit fluctuations. In Phase II, the distance experiences an exponential growth characterized by the MLE, which corresponds to a linear segment in the semi-logarithmic plot. Finally, in Phase III the separation of the trajectories is saturated and the distance converges to a constant value.

As proved by [13] and [16, 17, 15], the MLE can be estimated from the mean logarithmic separation rate of the trajectories in Phase II as

[TABLE]

where $\lambda_{k}$ is the estimated MLE, $k$ is the lagged time steps for the MLE estimation, $\Delta t$ is the time duration for each time step, $\boldsymbol{X}_{n}$ and $\boldsymbol{X}_{m(n)}$ are the initial points for the MLE estimation on the original and neighboring trajectories, respectively, $\boldsymbol{X}_{n+k}$ is the $k$th point behind $\boldsymbol{X}_{n}$ on the original trajectory, $\boldsymbol{X}_{m(n)+k}$ is the $k$th point behind $\boldsymbol{X}_{m(n)}$ on the neighboring trajectory, $d\left(m(n),n,0\right)$ is the Euclidean distance between the MLE estimation initial points, $d\left(m(n),n,k\right)$ is the Euclidean distance between the $k$th points behind the MLE estimation initial points, and $||A-B||$ denotes the Euclidean distance between points $A$ and $B$.

It should be noted that the original and the neighboring trajectories are usually from the same observed time series with different initial points, as illustrated in Fig. 2. In order to make sure the trajectories are temporally separated and thus can be seen as different trajectories, the trajectory initial points should satisfy $|m(0)-0|>w$, where $w$ is called the Theiler window [13] and should be determined according to the characteristics of the system.

III LE-Based Rotor Angle Stability Assessment

Here we propose a rotor angle stability assessment approach based on MLE estimation. It is data-driven and can perform online stability assessment only using the rotor angle and rotor speed of the generators [2, 3]. Although the rotor angle and rotor speed may not be directly available from PMU measurements, they can be estimated by various dynamic state estimation methods [18, 19, 20, 21, 22].

III-A State Variable Selection

Even a moderate-size power system may still have hundreds of state variables. Due to both the calculation intractability and the insufficiency of measurements, it is impractical to use all these variables to form the state space for MLE estimation. It is more feasible to reconstruct the power system dynamics only with a small number of state variables.

According to Takens’ theorem [23, 24], the dynamics of a nonlinear system can be reconstructed by the observations of a single state variable, and the reconstructed state vector can be expressed as

[TABLE]

where $N$ is the dimension of the original system, $\tau$ is the lag time, and $\theta_{t}$, $\theta_{t-\tau}$, $\theta_{t-2\tau}$, $\cdots$, $\theta_{t-2N\tau}$ are the observations of the observed variable at the corresponding time steps.

In the reconstructed state space, the Euclidean distance in (II) becomes

[TABLE]

where $\boldsymbol{\Theta}_{n+i}$ and $\boldsymbol{\Theta}_{m(n)+i}$ are the observations on the reconstructed trajectory.

Moreover, another simplification can be made based on [13] and [17], in which it is shown that it is sufficient to only consider the first component of the reconstructed state vector to estimate MLE, because all components will grow exponentially at the rate of MLE. Therefore, the $\boldsymbol{\Theta}_{n+i}$ and $\boldsymbol{\Theta}_{m(n)+i}$ in (3) can be replaced by their first components, and the distance in (II) can be calculated by

[TABLE]

In this paper, the relative rotor angle of the severely disturbed generator pair (SDGP) is selected as the observed state variable, because these SDGPs are, in general, responsible for the system dynamics after considerable disturbances [25]. An SDGP should be composed of a severely disturbed generator and the least disturbed generator in order to reflect the dynamics of the severely disturbed generators [26]. In particular, the SDGPs can be identified as follows.

Obtain the rotor speed of all generators at the fault clearing moment, $\omega_{t_{c},n},\ n=1,2,\dots,N_{G}$, from PMU measurements, where $N_{G}$ is the number of generators. 2. 2.

Obtain the maximal absolute value of the rotor speed, i.e., $\omega_{t_{c}}^{*}=\max\limits_{n=1,2,\ldots,N_{G}}\left|\omega_{t_{c},n}\right|$. Define generator $g$ as one of the severely disturbed generators, if $\left|\omega_{t_{c},g}\right|/\omega_{t_{c}}^{*}>\sigma$, where $\sigma$ is a predetermined threshold. In this paper $\sigma$ is chosen as 0.7, as in [27]. 3. 3.

Find the least disturbed generator with the minimal absolute value of the rotor speed. 4. 4.

Form a SDGP by combing one of the severely disturbed generators and the least disturbed generator. Iterate over the severely disturbed generators and form all SDGPs.

III-B RLS-Based MLE Estimation

The MLE can be estimated by calculating the slope of the logarithmic distance curve in Phase II. Considering the influences of measurement errors and nonlinear fluctuations, we adopt the least square algorithm to estimate the MLE.

From (II) and Fig. 2 it is seen that the MLE is estimated starting from time step $m(n)$. From this time step, $k+1$ sequential logarithmic distances can be obtained as

[TABLE]

Then, the MLE estimation model can be expressed as

[TABLE]

where $\lambda_{k}$ is the MLE to be estimated, $C_{k}$ is the constant term, and $\xi_{k}$ is the residual term.

The solution of the least square estimation is

[TABLE]

where $\boldsymbol{X}_{k}$ is the coefficient matrix and $\boldsymbol{Y}_{k}$ is the observation vector. According to the aforementioned definitions, $\boldsymbol{X}_{k}$ and $\boldsymbol{Y}_{k}$ can be expressed as

[TABLE]

Moreover, to avoid the repeated calculation of the inverse matrix in (9), a recursive estimation algorithm is applied [28], which can be formulated as

[TABLE]

where $\boldsymbol{x}_{k+1}$ is $\left[(m(n)+k+1)\Delta t\ \ 1\right]^{\top}$, $y_{k+1}$ is the new observation $L(m(n)+k+1)$, $\hat{\boldsymbol{E}}_{k}$ and $\hat{\boldsymbol{E}}_{k+1}$ are the estimation results before and after obtaining the new observation, $\boldsymbol{P}_{k}$ and $\boldsymbol{P}_{k+1}$ are the covariance matrices, and $\boldsymbol{G}_{k+1}$ is the gain vector.

The algorithm includes the following three steps:

With the first 2 groups of data ($k=1$ in (III-B)), set the initial values of $\hat{\boldsymbol{E}}_{1}$ and $\boldsymbol{P}_{1}$ to be $(\boldsymbol{X}_{1}^{\top}\boldsymbol{X}_{1})^{-1}\boldsymbol{X}_{1}^{\top}\boldsymbol{Y}_{1}$ and $(\boldsymbol{X}_{1}^{\top}\boldsymbol{X}_{1})^{-1}$, respectively. 2. 2.

Obtain a new observation of the logarithmic distance between the trajectories, and then sequentially calculate $\boldsymbol{G}_{k+1}$, $\hat{\boldsymbol{E}}_{k+1}$, and $\boldsymbol{P}_{k+1}$ according to (III-B). 3. 3.

Set $k=k+1$ and return to step 2.

III-C Parameter Setting

Because the Theiler window $w$ and the MLE estimation initial time step $m(n)$ determine the shape of the estimated MLE curve, they are crucial for a quick and reliable rotor angle stability assessment. Here, we discuss how to choose these parameters according to the rotor angle swing features.

(1) Theiler window selection

The Theiler window $w$ determines the temporal separation between the initial points $\theta_{0}$ and $\theta_{m(0)}$. It should be large enough to ensure that $\theta_{0}$ and $\theta_{m(0)}$ are the initial points of different trajectories. However, too large $w$ will cause unnecessary waiting time and delay the stability assessment.

The post-fault rotor angles of the SDGPs have significant swing patterns [29, 30, 31]. According to the features of the relative rotor speed of the SDGPs, six distinct swing patterns can be identified as shown in Fig. 3. Different $w$ are chosen for different patterns as follows.

• •

Pattern I: In this pattern, the relative rotor speed increases after fault clearing. No decelerating area exists for the SDGP and the system will lose stability during the first swing. Since the original and the neighboring trajectories separate rapidly, $w$ is set to 1 (the smallest positive integer) to minimize the estimation waiting time.

• •

Pattern II: In this pattern, the relative rotor speed decreases after fault clearing. However, since the decelerating area is relatively small, the relative rotor speed increases again after a short time period, and the increasing trend continues until the system loses stability. The key feature of this pattern is that the relative rotor speed at the fault clearing moment, denoted by $v_{0}$ in Fig. 3, appears again after the initial decrease. In order to achieve obvious separation between the trajectories, $w$ is set to be the time step lags of the reappearance of $v_{0}$.

• •

Pattern III: In this pattern, the relative rotor speed first decreases to $-v_{0}$, and then exhibits periodic oscillations. The decelerating area is large enough to reduce the relative rotor speed to zero, and the system stability depends on the damping characteristics of the post-fault equilibrium point. In this case, $w$ is set to be the time step lags of the first appearance of $-v_{0}$.

• •

Pattern IV: In this pattern, the relative rotor speed first decreases to some value greater than $-v_{0}$, and then oscillates periodically. The key feature of this pattern is that the relative rotor speed $v_{0}$ and $-v_{0}$ cannot be observed after fault clearing. This pattern is a special case of Pattern III, and $w$ is set to be the time step lags of the first appearance of the local minimal relative rotor speed after fault clearing.

• •

Pattern V: This pattern is similar to Pattern III, and usually appears after very quick fault clearing. The key feature of this pattern is that the relative rotor speed shows decelerated growth immediately after fault clearing, and the relative rotor speed $-v_{0}$ can be observed after that. The $w$ is set in the same way as in Pattern III.

• •

Pattern VI: This pattern is a special case of Pattern V. The key feature of this pattern is that the relative rotor speed shows decelerated growth immediately after the fault clearing, and the relative rotor speed $-v_{0}$ cannot be observed during the oscillations. The $w$ is set in the same way as in Pattern IV.

(2) MLE estimation initial time step selection

The MLE estimation initial time step $m(n)$ should ensure that the slope estimation is performed for Phase II of the logarithmic distance growth. In this phase, the original and neighboring trajectories have been sufficiently separated, and thus the logarithmic distance curve has shown clear development trend [13].

In fact, the distance between the rotor angle trajectories has a close relationship with the relative rotor speed of the SDGP. According to the definition of the Theiler window, the distance between the corresponding points of the original and the neighboring trajectories at any time step $j$ can be reformulated as

[TABLE]

where $v_{t}$ is the relative rotor speed at the relevant time step.

It is revealed in (III-C) that the distance $d_{j}$ is equal to the area enclosed by the time axis and the relative rotor speed curve within the Theiler window, as indicated by the shaded regions in Fig. 3. Therefore, $m(n)$ can be determined as follows.

• •

Patterns I–II: With the selected $w$, it is seen in Fig. 3 that the distance $d_{j}$ of these patterns monotonically increases after fault clearing, which indicates that the logarithmic distance growth is in Phase II as soon as the fault is cleared. Therefore, in these two patterns, $m(n)$ is set to be $w$ to minimize the estimation waiting time.

• •

Patterns III–VI: For these patterns the distance will exhibit periodic fluctuations after fault clearing. To ensure the trajectories have been sufficiently separated and to catch the main trend of the fluctuations, $m(n)$ is set to be $w+j^{*}$, where $j^{*}$ is the time step when $d_{j}$ reaches its first local maximum after fault clearing.

III-D Rotor Angle Stability Assessment Criteria

When the MLE curve of a SDGP is estimated, the following criteria can be used to determine the stability of the SDGP.

• •

Criterion I: If the MLE of the SDGP increases at the beginning, the SDGP is unstable.

• •

Criterion II: If the MLE decreases at the beginning, it will have oscillations. If the first peak point of the oscillation is positive, the SDGP is unstable.

• •

Criterion III: If the MLE decreases at the beginning and the first peak point of the oscillation is negative, the SDGP is stable.

Typical MLE curves corresponding to these criteria are shown in Fig. 4. If the condition in Criterion I is satisfied, Pattern I or II in Fig. 3 will happen, for which the distance will increase immediately after the MLE estimation initial time step and the trend will last until the SDGP loses stability.

Criteria II and III correspond to Patterns III–VI. In these patterns, the relative rotor angle will oscillate after fault clearing. However, by using the selected parameters, the logarithmic distance will exhibit significant development trend as illustrated in Fig. 5. If the relative rotor speed has a undamped oscillation, the SDGP is unstable and the logarithmic distance curve will fluctuate periodically as shown in Fig. 5(a). Since the MLE is the average slope of the logarithmic distance curve from the MLE estimation initial time step, it will first decrease and then increase to a positive peak value, as in Fig. 4. By contrast, if the relative rotor speed has a damped oscillation, the SDGP is stable and the logarithmic distance curve will look like Fig. 5(b). The MLE will first decrease and then increase to a negative peak value, as in Fig. 4.

Because SDGPs are responsible for the system dynamics after disturbances [25], we finally have the following criterion on the angle stability of the system.

• •

Criterion IV: If all SDGPs are stable, the system is stable; otherwise, the system is unstable.

III-E Rotor Angle Stability Assessment Procedure

The proposed online rotor angle stability assessment procedure is shown in Fig. 6, which includes the measurement data preparation module, the parameter setting module, the MLE estimation module, and the stability assessment module. Specifically,

When a fault is detected, the measurement data preparation module will be immediately activated to identify the SDGPs and collect the corresponding relative rotor angle and rotor speed measurements in real time. 2. 2.

The relative rotor speed variation patterns are identified online according to the rules described in Section III-C, and the parameters are assigned accordingly. 3. 3.

The MLE sequences can be calculated by using the RLS-based algorithm in Section III-B. 4. 4.

Finally the stability condition can be assessed by using the criteria provided in Section III-D according to the features of the estimated MLE curves.

IV Case Studies

The proposed approach is tested on the New-England 39-bus system and the NPCC 140-bus system. Simulations are performed with Power System Analysis Toolbox (PSAT) [32] in MATLAB.

IV-A New-England 39-Bus System

The New-England 39-bus system has 10 generators and 46 branches. The parameters can be found in [33]. Unless otherwise specified, all generators in the tests are described by the fourth-order transient model with Type I turbine governor (TG), Type II automatic voltage regulator (AVR), and Type II power system stabilizer (PSS) (see PSAT documentation). All loads are described by the ZIP model and the ratios of the constant impedance, constant current, and constant power loads are $0.4$, $0.5$, and $0.1$, respectively. The sampling rate of the PMU measurements used for MLE estimation is 120 samples/s [3].

To verify the effectiveness of the proposed approach, a three-phase-to-ground fault is applied at bus 2 at $t=1$ s, and the fault is cleared by opening line $2$–$3$ at $t_{c}=1.243$ s and $t_{c}=1.244$ s for Scenarios I and II, respectively. According to the time-domain simulation, the system is stable under Scenario I and unstable under Scenario II.

By using the method in Section III-A, the generator pairs $38$–$39$ and $37$–$39$ are identified as the SDGPs for both scenarios. Figs. 7–8 show the relative rotor angles and the estimated MLEs of the SDGPs.

From Fig. 7(a), it is seen that the relative rotor angles of the SDGPs tend to be stable after a long period of oscillations. In Fig. 7(b), the MLEs of the SDGPs decrease immediately after fault clearing, and their first peak points of the oscillations are both less than 0. Therefore, according to Criteria III and IV, the system is stable under this scenario.

By contrast, under Scenario II the MLE of the generator pair $38$–$39$ increases immediately after fault clearing, as shown in Fig. 8(b). Therefore, according to Criteria I and IV, the system is unstable under this scenario.

For the two scenarios, the system stability can be assessed, respectively, within $2.82$ s and $1.40$ s after fault clearing, which demonstrates that the proposed approach can provide early detection of instability. In fact, according to the proposed approach, the first-swing instability can be identified very quickly, because the instability can be detected at the beginning of the MLE curve. On the other hand, for the multi-swing stability or instability, the assessment time is a little longer since the approach needs to check the first peak point of the MLE curve after fault clearing.

In industrial applications, a predetermined relative rotor angle value, i.e., $\pi$ rad, is usually set as the threshold for determining rotor angle stability [34]. Although this pragmatic criterion is easy to execute, the relative rotor angle corresponding to rotor angle instability may change significantly with the change of topologies, parameters, and operating conditions. For instance, under Scenario I, the relative rotor angle between generators $38$ and $39$ can reach up to $3.207$ rad while the system is still stable.

It should also be noted that there are foundational differences between the proposed approach and that in [3]. In the proposed approach, the system stability can be explicitly determined at latest when the first peak point of the MLE curve is observed. By contrast, for the approach in [3], the MLE curve must be observed for a long period to ensure the sign of MLE. Unfortunately, it is actually difficult for the approach in [3] to predetermine the observing window size to provide reliable and timely assessment results (for instance, see Figs. 1–2 in [3]). From this perspective, the proposed approach is more time efficient and more reliable.

In order to further verify the accuracy of the proposed approach, extensive tests are performed on the New-England 39-bus system. Specifically, a three-phase-to-ground fault is created for each bus (except the generator buses) at $t=1$ s, and is cleared at $1.08$ s, $1.16$ s, $1.24$ s, and $1.32$ s, respectively. According to the test results, the proposed approach successfully determines the rotor angle stability in all $224$ tests. The occurrence frequencies of the fault patterns and the success rate of the proposed stability assessment approach are summarized in Table I.

It is seen in Table I that Patterns I and II those correspond to first-swing instability usually occur with longer fault clearing time. In contrast, Patterns V and IV usually occur with the shorter fault clearing time, and the system will have a good chance to maintain stability in these patterns. Among the tests, the first-swing instability identification time ranges from $1.2$ s to $1.5$ s, and the multi-swing stability assessment time ranges from $2.2$ s to $2.5$ s, which further confirms the efficient of the proposed approach.

Because the system has sufficient damping to ensure the stability of the post-fault equilibrium points, there are no multi-swing instability cases. In order to generate a multi-swing instability case, all of the PSSs are removed and the parameters of the AVRs are tuned. In the modified system, a three-phase-to-ground fault is applied at bus 28 at $t=1$ s, and is cleared by opening line $27$–$28$ at $t_{c}=1.12$ s. In this test, only generator pair $38$–$39$ is identified as the SDGP, whose relative rotor angle curve and MLE curve are shown in Fig. 9.

In Fig. 9(b), the MLE decreases immediately after fault clearing, and the following peak MLE is positive. According to Criteria II and IV, the MLE variation pattern indicates that the system is multi-swing unstable, which is verified by the relative rotor angle curve. In this case, the stability assessment time is $2.36$ s, which is much shorter than the time required for directly looking at the relative rotor angle curve.

IV-B NPCC 140-Bus System

The proposed approach is also tested on the NPCC 140-bus system [35]. The settings are the same as those for the New-England 39-bus system.

A three-phase-to-ground fault is first applied at bus $35$ at $t=0.1$ s, and is cleared by opening line $34$–$35$ at $t_{c}=0.307$ s (Scenario III, stable) and $t_{c}=0.308$ s (Scenario IV, unstable), respectively. Generator pairs $1$–$48$ and $2$–$48$ are identified as the SDGPs for both cases. The relative rotor angle curves and MLE curves under these two scenarios are shown in Figs. 10–11, respectively.

In Fig. 10, it is seen that the curves of Scenario III are similar to those of Scenario I. In this case, the system is determined as stable because of the same reason as Scenario I. As shown in Fig. 11, under Scenario IV the MLE curves of both SDGPs exhibit unstable features, which indicates that both generator pairs will lose stability and thus the system is unstable.

Extensive tests are also executed on the NPCC 140-bus system. A three-phase-to-ground fault is applied at each bus (except the generator buses) at $t=0.1$ s, and is cleared at $0.18$ s, $0.26$ s, $0.32$ s and $0.40$ s, respectively. Table II lists the occurrence frequencies of the fault patterns and the success rate of the proposed stability assessment approach during the tests.

According to the test results, the proposed approach can accurately determine system stability in all $716$ tests. The first-swing instability identification time ranges from $1.1$ s to $1.7$ s, and the multi-swing stability assessment time ranges from $1.8$ s to $2.4$ s.

V Conclusion

In this paper, a model-free approach for online rotor angle stability assessment is proposed based on MLE. By using the proposed MLE estimation algorithm, parameter setting rules and the stability criteria, the approach can online identify the system stability condition with PMU measurements. The approach does not need a predetermined observing window to identify the sign of the MLE, and can provide reliable and timely assessment results by analyzing the features of the estimated MLE curve. To verify the performance of the proposed approach, extensive tests are performed on the New-England 39-bus system and the NPCC 140-bus system. The proposed approach can successfully determine the system stability conditions in all 945 tests. Moreover, among all the tests, the first-swing stability can be assessed within $1.7$ s and the multi-swing stability can be assessed within $2.5$ s.