(Zheng et al., 2021a) defines an sBTTC index that monotonically represents the system stability condition based on wide-area branch response information. Based on this, Eq. 5 represents the simplified sBTTC index. The deterioration of system stability is reflected by setting a stability threshold. sBTTC = U i m U i n ⁡ cos 2 Δ θ 2 (5)

Where: U _im and U _in represent the node voltage amplitude at the two ends of the i branch, respectively. Δ θ signifies the difference in the mode value of the node voltage vector at the two ends of the i branch.

In the destabilized system, the phase difference between the nodes at the two ends of the key branch shows a continuous increase and tends to be unbounded, while the voltage product keeps decreasing. This characteristic is reflected in Eq. 5, which shows that the phase difference between the two ends of the branch is 180° when the sBTTC index drops to 0, indicating system instability. During the initial descent phase of the sBTTC index, the sBTTC index values of all observed branches in the network are ordered from smallest to largest. According to Eq. 6, the key branch k has the smallest index value, and the candidate the key branch clusters have the next smallest index value. Additionally, the branch droop voltage coefficient ξ _v is introduced and combined with the branch droop voltage coefficient located within the interval [1,2] to distinguish the key branches (Zheng et al., 2021b; Zheng et al., 2022). Meanwhile, in order to improve the accuracy of discriminating the key branch, the branch droop voltage coefficient ξ _v is introduced and combined with the branch droop voltage coefficient located in the interval [1,2] for discrimination (Zheng et al., 2021a; Zheng et al., 2022). k = argmin i ∈ 1,2 , ⋯ N sBTTC i (6) 1 ≤ ξ v ≤ 2 (7)

When a system generator destabilisation occurs, at least one branch exists in the network. The phase difference between the nodes at the ends of this branch is consistent with the change in the power angle difference between the overrun and lagging units, and the centre of oscillation is located here. Among the AC lines in the whole network, the key branches have the weakest stability performance. The dynamic response process of sBTTC index of the key branch can directly reflect the transient angle stability of the whole system.

The sBTTC exponential variation curves for the key branch B3-B4 are illustrated in Figure 4 for the single machine infinity system in Section 2 under both the instability and stability cases.

The sBTTC exponential curve of branch B3-B4 is transformed to the sBTTC-d (sBTTC)/dt plane as shown in Figure 5. Calculate d (sBTTC)/dt as in Eq. 8: d sBTTC d t m Δ t , i = sBTT C m Δ t , i − sBTT C m − 1 Δ t , i Δ t (8) where d (sBTTC)/dt and sBTTC_(m−1)Δt,i represent the sBTTC index values of branch i at time instances mΔt and (m−1) Δt, respectively.

As can be seen from Figure 5A, if the system transient angle is stability, the trajectory of sBTTC-d (sBTTC)/dt curve will gradually converge to a stable point, at which time the value of sBTTC-d (sBTTC)/dt is close to zero; from Figure 5B if the system transient angle is instability, the curve of sBTTC-d (sBTTC)/dt will be diverging around the zero axis, and cannot be converted to a stability equilibrium point. Therefore, the system angle stability analysis problem can be transformed into the problem of whether the sBTTC exponential curve can converge to a certain point. With the help of MLE, it can be used to determine the motion behavior of the power system; a negative MLE characterizes the convergence of trajectories, which represents the stability of the power system; a positive MLE characterizes the divergence of trajectories, which represents the instability of the power system.

The selection of the initial time window requires ensuring that the distance change between the original and neighboring trajectories exhibits clear dynamic characteristics. In Section 3.1.2, the sBTTC-d (sBTTC)/dt plane trajectories have obvious dynamic convergence or divergence characteristics at the beginning. Both d (sBTTC)/dt and lgL, as calculated by Eq. 8 and 9 respectively, closely relate to the disparity between two points along the trajectories. Therefore, the beginning stage of original and neighbouring trajectory calculation also has obvious dynamic characteristics. In this paper, the initial value of the initial time window is chosen as the initial moment of the MLE calculation.

To examine the impact of different initial time window lengths on the MLE-t curve, we conducted calculations for two scenarios: M = 5 and M = 50, focusing on the stability and destabilization cases of a two-machine system. The resulting MLE-t curves for the sBTTC indices of key branches B3-B4 are presented in Figures 7A, B.

Under the stability and instability calculations for the stand-alone system, in order to analyse the effect of the selection of the initial time window length on the MLE-t curves, the initial window length M = 5 as well as M = 50 is set for the calculation. Figure 7 demonstrates the MLE-t curves for the sBTTC index of the B3-B4 key branch.

From Figure 7, the MLE-t curve can be divided into four phases:

(i) Initial stage: marked by significant fluctuations occurring in a brief span of time.

In the stability and instability cases, both M = 5 and M = 50 have no effect on the final sign of the MLE-t curve. Regardless of the length of the initial time window, the MLE-t curve eventually tends to be negative in the instability case and positive in the stability case. In this paper, the initial time window length M is taken as shown in Eq. 12, which can be adapted to systems of different sizes. M = 5 int lg ⁡ n bus (12)

Where: n _bus denotes the number of system nodes, and int represents upward rounding.

Maximum Lyapunov Exponent to discriminate the stability of the system needs to judge the final sign characteristics of the MLE. A time observation window of finite window length is usually predefined. The optimal time window depends on many factors and is difficult to determine. Theoretically, the length of the observation window tends to infinity to reflect the final MLE of the system. However, it is impossible to calculate the MLE at infinite moments due to the limitation of practical computational conditions.

If the optimal time observation window is too long, it may lead to the instability system crossing the unstable equilibrium point long time ago, and then the control action is implemented too late; if the optimal time observation window is too short, it may lead to the observation window not reflecting the final situation of the MLE-t curve, and then there is a misjudgment.

The sBTTC index can reflect the system stability level. Therefore, the initial moment of the optimal time window in this paper is set as the beginning moment of the initial decline phase of the sBTTC index, and the ending moment is the end moment of the initial decline phase of the sBTTC index. The dynamic observation time window set in this paper can be adapted to a variety of fault scenarios and can effectively discriminate the final sign of MLE.

It should be noted that, under the complex multivariable system, if the MLE sign is clear within the optimal time window, the observation of the MLE curve does not need to wait until the end of the optimal time window. On the contrary, if the MLE curve within the observation window has not entered the oscillation phase, the MLE sign cannot be determined. In this case, it is necessary to observe the value of the sBTTC index at the end of the dynamic observation window to make a determination.

The aim is to reduce the effect of noise on the signal and improve the curve smoothness. This paper employs the Savitzky-Golay filtering signal processing technique to refine the calculation of the MLE.

Following the fault occurrence, the sBTTC index information of the key branch can be acquired through real-time PMU sampling. Utilizing Eq. 9, the logarithmic Euclidean variation distance of the sBTTC index generates a one-dimensional signal sequence of length N, denoted as lgL lgL = [lgL ₁, lgL ₂, …. LgL _N ]. To apply the filter, a window of length 2d + 1 is selected, resulting in the subsequence lgL_w = [lgL _n-d , …lgL _n, …lgL _n+d ], where lgL _n serves as the central data. The order of polynomial fitting is set as p, and the window’s subsequence is represented by the fitting polynomial, as shown in Eq. 13: lg ⁡ L _ w ≈ a 0 + a 1 X + a 2 X 2 + . . . + a p X p (13)

Where: X is the relative position within the window, i.e., − d , − d + 1 , … d .

The objective function corresponds to the squared discrepancy between the subsequence fitting polynomial and the original signal subsequence, expressed as shown in Eq. 14. min ⁡ = ∑ lg ⁡ L _ w − a 0 + a 1 X + . . . + a p X p 2 (14)

The fitting polynomial coefficients, denoted as a 0 , a 1 , … a p , can be obtained by achieving an objective function min value of 0 as shown in Eq. 14. Utilizing a₀, the smoothed value of the central data can be estimated. Thereafter, employing a sliding window approach to each position of the signal, the entire signal sequence lgL is comprehensively filtered. Subsequently, Eq. 11 is employed to calculate the MLE value of the sBTTC exponential trajectory.

A three-phase grounded short-circuit fault occurs on the front side of the bus16-bus17 branch at 0.2 s. The fault removal moment is set to 0.4 s. The system is transiently stable. The sBTTC index of the wide-area branch circuit based on WAMS measurement is shown in Figure 10.

In Figure 10A, it is evident that the sBTTC index experiences the initial decline phase from t = 0.8 s to t = 2.03 s, which represents the MLE observation time window. Within this window, the bus15-bus16 branch exhibits the lowest sBTTC index value. The droop voltage position coefficient of this branch is illustrated in Figure 10B, and it satisfies the conditions stated in Eq. 7, confirming its key ity.

Next, the MLE curve of the sBTTC index trajectory for the bus15-bus16 branch is determined and depicted in Figure 11. At 1.33 s, the MLE curve exhibits its first swing vertex with a negative sign, while the value of the sBTTC index at that moment is 0.88, which is greater than the threshold value of 0.8 set previously. Based on these observations, we can conclude that the system is in a state of transient stability.

At 0.2 s, a three-phase grounded short-circuit fault emerges on the front side of the bus16-bus17 branch, with the fault removal scheduled for 0.6 s. Consequently, the system experiences a single-swing destabilisation. In Figure 12, we observe the response curves of the generator angle and the entire network node voltage for the node system, with G30 serving as the reference machine.

Figures 12A, B clearly indicates that the generator angle loses synchronisation after the fault, leading to a significant decline in node voltage, followed by continuous periodic oscillations. The system units divide into two groups: the leading group S, comprising 6 generators, and the lagging group A, comprising the remaining 4 generators, represented as S = {G31, G32, G33, G34, G35, G36} and A = {G30, G37, G38, G39}, respectively.

Upon calculating the sBTTC indices of the branches based on the disturbed response information and identifying the key branches, we assess whether their branch droop voltage location coefficients fall within the interval [1, 2], as depicted in Figure 13.

From Figure 13A, we can observe that within the time frame of 0.92 s–4.66 s, the branch sBTTC index consistently decreases to a value of zero. Throughout this process, the Bus14-Bus15 branch maintains the minimum sBTTC index value, and the branch footing voltage is situated atop this branch, confirming it as the key branch of the system.

Subsequently, we obtain the sBTTC index traces of the key branch and compute the MLE curve, as shown in Figure 14. Also, the MLE calculation curves for the traditional method of obtaining the most overrun unit G35 are given in Figure 14 (Dasgupta et al., 2015).

In Figure 14A, the system stability is discriminated using the method proposed in this paper. The initial phase of the MLE-t curve is not given due to the large fluctuation, and the oscillation phase has long entered in the observation window. The first swing vertex of the MLE curve can be detected quickly at point C, and the value of the first swing vertex is 0.21. The value of the sBTTC index at this moment is 0.71, which is smaller than the set threshold value of 0.8. Therefore, this paper determines that the system is unstable for 1.3 s, and the maximum relative angle difference of the generator at this time is 1.91 rad.

In Figure 14B, the traditional method is used to calculate the MLE and determine the stability of the system. If the observation window is set to 3 s, the MLE sign can be determined to be positive at 2.6 s, and then the system is judged to be unstable. If the observation window is set too short, such as 2 s, the system will be misjudged as stable.

In Figure 14C, the time required to discriminate the system stability using the traditional stability engineering criterion (generator rotor angle over 180) is 3.4 s.

In summary, the time required to judge the stability for MLE analysis based on the SBTTC exponential curve is 1.3 s, the time required to judge the stability for MLE analysis based on the generator rotor angle curve is 2.6 s, and the time required to judge the stability of the system based on the traditional engineering stability criterion is 3.4 s. Therefore, the method proposed in this paper is characterised by its rapidity in determining the stability of the system.

Considering 0.1 s communication delay, emergency control measures are implemented at 1.4 s and 2.7 s respectively. Cutting out the most overrun unit G35, the system generator angle curve is obtained as shown in Figures 15A, B. After the implementation of emergency control measures at 1.4 s, the generator group can maintain synchronous operation. After the implementation of emergency control measures in 2.7 s, the generator group loses synchronous stability and is unable to maintain system stability.

The method proposed in this study incorporates the Savitzky-Golay filtering processing technique for MLE calculation, endowing it with a considerable degree of noise immunity. Typically, the signal noise ratio (SNR) of PMU measurements surpasses 100, ensuring that the measurement error remains below 1 percent (Dasgupta and Paramasivam, 2013). A Gaussian noise with a SNR of 60 dB is added to the above instability example to simulate the actual PMU measurement noise. The results of the MLE curves for both algorithms are shown in Figure 16.

From Figure 16A, in the method proposed in this paper, the addition of noise does not have a delayed effect on the time point corresponding to when the MLE curve reaches the apex of the first swingback in the discriminant. There is no effect on the time required for the stability discrimination.

From Figure 16B, in the traditional MLE method, the MLE curve is affected by noise with a larger magnitude in the first 5 s, when it is not appropriate to discriminate the MLE sign. Until the 5 s curve is slightly smooth, then the MLE symbols can be discriminated, but the discriminatory timeliness is poor.

From this, the method proposed in this paper has strong anti-noise ability.

To further validate the reliability of the proposed method, a traversing fault test is performed at the 39-node system. So that all lines excluding the generator node have a three-phase short-circuit fault (N−1) at the moment of 0.2 s. The double return line defaults to disconnect the first return line. The fault removal moments are traversed sequentially from 0.2 to 3 s in steps of 0.2 s. Total 480 groups of faults.

In each set of tests, the key branch was quickly identified after fault removal using the method proposed in Section 4.3. The MLE curve of the sBTTC index of the key branch is analyzed to determine the system stability. Comparison is made with two judgments. These two criteria are the criterion based on the stability characteristic quantity (generator angle greater than 180°) and the criterion based on the traditional MLE calculation method. To increase the conservatism, the time window based on the traditional MLE calculation method is chosen as the end moment of the simulation. The third of these methods uses the batBPA batch program, and the first and second methods are implemented by building a Python interface to BPA. The statistical results are shown in the following table.

From Table 1, the method of this paper can successfully identify 404 sets of destabilizing samples through the traversal failure test. If based on the amount of stabilization features, the recognition rate of the destabilization samples of this paper’s criterion is 100%. The stabilization misjudgment rate is 2.6%. And in the test results, the average time required to judge the system stability is within 1.4 s after failure. The first method in the above table is obviously better than the second.