Transient stability analysis based on the state-plane trajectory and the relative kinetic energy change rate

Due to the widespread application of the wide area measurement system in power systems, it has become feasible to identify the transient stability of a power system on-line. A real-time transient stability identification method combining stateplane trajectory and relative kinetic energy is proposed. First, the correlation between the minimum value of the state-plane trajectory and the system stability is analysed in a single-machine infinite system. On this basis, the relationship between the concavity–convexity of the trajectory and the transient stability is discussed in detail. Second, the calculation of the transient energy is introduced to analyse the energy transformation and a sub-criterion is presented. The comprehensive application of trajectory characteristics and energy criteria can accurately determine the transient stability of all of the units, which is of great practical significance for accurately cutting off the unstable units when system faults occur. A simulation analysis with the IEEE New England 10-machine 39-bus system verified the correctness and effectiveness of the proposed method.


Introduction
Modern power systems have evolved into large-scale regional grid interconnection systems.Although the interconnection of power grids brings huge economic benefits, it also makes the transient stability problem more complex.Transient instability is still one of the greatest threats to a modern power system; thus, effective realtime transient stability prediction and emergency control are required.
The traditional transient stability control strategy mainly adopts the method of off-line control [1].However, this method has limitations.At present, common transient stability analyses include time-domain simulation methods [2,3], direct methods based on the Lyapunov energy function [4,5], a mixed method combining a numerical simulation method with the direct method [6] and the artificial intelligence algorithm [7].
Due to the widespread application of the wide area measurement system (WAMS) [8,9], the identification of transient instability in a real-time power system based on a phasor measurement unit (PMU)/WAMS has become a hot topic.In [10], the validity of the method for stable identification using the stateplane trajectory concavity-convexity proposed in [11] was verified in the one-machine-infinite-bus (OMIB) system.Cen et al. [12] used the least-squares method to fit the first derivative of the stateplane trajectory and judged the transient stability of the system based on the change trend of the first derivative.Zheng et al. [13] proposed a state-plane trajectory-stabilisation method based on the dimensional reduction of the power angle space and reduced the dimension of the power angle trajectory of the system.Li et al. [14] suggested a new transient stability assessment approach based on a characteristic of the key ac branch's disturbed state-plane trajectory.In [15], the trajectory on the plane consisting of the perturbation energy and the voltage phase angle has similar geometric characteristics with the state-plane trajectory, and a transient instability identification method based on the stable boundary of the measured response trajectory was proposed.Zhang et al. [5] proposed a practical dynamic model with good accuracy for transient stability analysis.The critical energy is calculated and the stability of the system can be judged by comparing the transient energy with it.
The transient stability detection method based on the trajectory geometric characteristics of WAMS has no dynamic network equivalence and is easy to implement.It only needs information on the two latest time intervals to detect real-time transient instability.Since the criterion is derived from a two-dimensional autonomous system, the judgment may be erroneous due to the non-Hamiltonian of a complex modern power system.Therefore, in this paper, the aim is to enhance the robustness of the trajectory geometric characteristics method to cope with a complex power system.The correlation between the minimum value of the stateplane trajectory and the system stability is analysed.On this basis, the relationship between the concavity-convexity of the trajectory and the transient stability is discussed.Introducing the transient energy function, the relationship between the variation rate of the transient unbalanced energy and the trajectory characteristics is deduced.A combination of the state-plane trajectory concavityconvexity and the relative kinetic energy change rate is formed to create a comprehensive criterion, and a new method for identifying transient stability is obtained.Finally, the proposed method is tested on the IEEE 39-bus system.The results indicate that the proposed method is effective and reliable in detecting the loss of synchronism.

Relationship between the state-plane trajectory minimum and stability
The basic structure of a simple OMIB system is shown in Fig. 1.
When a line fault occurs in the Hamiltonian OMIB system, by ignoring the damping and excluding the regulator and governor roles, the mathematical expression of the transient process of the generator is as follows: where δ is the generator angle, Δω is the angular velocity, ω 0 is the synchronous angular velocity -the power system generally takes ω 0 = 100π -M is the generator inertia, P M is the mechanical power of the generator and P max is the maximum of the electromagnetic power of the generator.
Assuming that angle δ is the abscissa and the angular velocity Δω is the ordinate, the generator's time-varying motion state can be represented on a plane and the power angle stability can be determined based on the state-plane trajectory characteristics.
When the fault is cleared, the generator may have the following two states: i.The generator does not have a deceleration area and its rotor is always in an accelerated state.The state-plane trajectory moves up and there is no minimum value, thereby resulting in system angle instability.ii.The generator has a deceleration area and the system angle stability depends on its size.After entering the deceleration area and before reaching the unstable equilibrium point, the angular velocity decreases and the state-plane trajectory follows a declining trend.If the system is stable, the stateplane trajectory will cross the horizontal axis.When the angular velocity reaches zero, the angle will not increase and reaches its maximum value and the state-plane trajectory will start to swing back.If the system is unstable, the trajectory cannot reach the horizontal axis and the system reaches an unstable equilibrium point.When the angular velocity reaches a minimum value, the state-plane trajectory starts to move upwards and the angle difference increases.Fig. 2 shows the state-plane trajectories for different fault clearing times when the generator has a decelerating area.
It is obvious that the stability of the system can be identified by whether there is a minimum value of the state-plane trajectory or not.
Proof: From (1), the first derivative of the angular velocity Δω on the power angle δ is Transforming (2) by integrating the two sides gives ∫ Mω 0 ΔωdΔω = ∫ (P M − P max sin δ)dδ which can be further simplified to where C is a constant and is related to the state at the moment of the system state change.
When operating normally, Δω(δ) = 0, where δ is a stable operating point at this time.In a post-fault system, the maximum value of the electromagnetic power is reduced to P 1max , and according to the continuity of the state-plane trajectory before and after the fault, the state-plane trajectory of the system after the fault becomes where C 1 is the value obtained from the continuity of the stateplane trajectory before and after the fault and is related to the motion state at the time of the fault.When the fault is cleared, Δω 1 and δ are known and the maximum value of the electromagnetic power is increased to P 2max .According to its continuity before and after the fault clearing, the state-plane trajectory can be written as where C 2 is the value obtained from the continuity of the stateplane trajectory before and after the fault clearing.The slower the fault clearing time, the greater the value of C 2 .Let f 1 (δ) = 2(P M δ + P 2max cos δ) /Mω 0 + C 2 ; when d f 1 (δ) /dδ = 0, we can obtain the following relationship in the state-plane trajectory: If P M > P 2max , the generator does not have a deceleration area.Equation ( 7) has no real solution; that is to say, ( 6) does not have a minimum value and Δω 2 will increase monotonically.
If P M < P 2max , then there is a decelerating area in the generator.The solution δ = sin −1 (P M /P 2max ) can be solved according to (7).Let δ satisfy (7) as an unstable equilibrium point (δ UEP ), then δ UEP = sin −1 (P M /P 2max ) can be obtained.Therefore, the state-plane trajectory has its minimum value at δ UEP .If f 1 (δ UEP ) > 0, the minimum value meets Δω 2 (δ UEP ) > 0. As the time increases, the power angle δ also increases, the state-plane trajectory starts to move upwards and the system is unstable.If f 1 (δ UEP ) < 0, (6) corresponds to Δω 2 (δ UEP ), i.e. there is no real solution.The stateplane trajectory has not yet reached an unstable equilibrium point and so begins to swing back through the horizontal axis, and the system is stable.In the case of f 1 (δ UEP ) = 0, then Δω 2 (δ UEP ) = 0 and the system reaches its stability limit.□

Identifying the stability based on the concavity-convexity of the state-plane trajectory
For a real-time online system, the aforementioned method requires that the stability of the system can be identified when the power angle reaches an unstable equilibrium point, which is very inefficient in terms of time.Therefore, it is necessary to improve the method, and it is hoped that the stability of the system can be identified earlier in the time domain.By referring to the existing literature and combining the method of stability discrimination in the previous section, the relationship between the concavity and convexity of the state-plane trajectory and the stability of the power angle is analysed.For the trajectory of the state-plane after the fault cleared, the second derivative of the angular velocity Δω 2 to the angle δ is given by When ( 8) is equal to 0, the solution becomes A stable state-plane trajectory does not always intersect with the concave-convex boundary within the concave region and an unstable state-plane trajectory is just the opposite.It will surely pass through the concave-convex boundary within the convex region.Therefore, as long as the state-plane trajectory curve intersects the concave-convex boundary, the system is unstable.The proof is as follows.
Let f 2 (δ) = (P m − P 2max sin δ) 2 / −Mω 0 P 2max cos δ satisfy f 2 (δ UEP ) = 0 when the concavity boundary intersects the horizontal axis.For any δ that satisfies f 1 ′(δ) − f 2 ′(δ) > 0, it means that the descent rate of the state-plane trajectory curve is less than the descent rate of the concave-convex boundary line.Thus, when f 1 (δ UEP ) > 0, according to the theory in Section 2.1, the system is unstable.Therefore, the key to the problem is simply to prove this for π/2 < δ < δ UEP and satisfy For any δ, f 1 ′(δ) and f 2 ′(δ) satisfying ( 10) Therefore, According to (11), when π/2 < δ < π, f 1 ′(δ) − f 2 ′(δ) > 0. Therefore, the tangent slope of the state-plane trajectory is always greater than that of the concave-convex boundary line, so the state-plane trajectory crosses the boundary line and enters the convex region.The decrease rate is slower than that of the state-plane trajectory.
Upon reaching the unstable equilibrium point δ UEP , the state-plane trajectory reaches a minimum value and the generator electromagnetic power is less than the mechanical power.At this time, the angular velocity does not decrease to zero and the system cannot maintain synchronous operation, as shown in Fig. 4.
Through this analysis, it can be concluded that as long as the state-plane trajectory passes through the concavity-convexity demarcation line, its minimum value after the system fault has cleared is always greater than zero, which means that the system power angle is unstable.If the state-plane trajectory of the system does not intersect with the concave-convex boundary, it will not cross the unstable equilibrium point and a backswing will occur, leading to stabilise the system power angle.
Therefore, when μ = Δω ⋅ l < 0, the state-plane trajectory is in the concave region and the system is transiently stable, but when μ = Δω ⋅ l > 0, the state-plane trajectory is in the convex region and the system is transiently destabilised.When μ = Δω ⋅ l = 0, the state-plane trajectory is located on the concave-convex boundary and the system is critically stable.

Application of the proposed approach in multi-machine systems
This section is focused on the identification of transient instability in the two-group first-rate instability mode of a multi-machine system.When large disturbances occur, the generators in the power system can be divided into a leading group S and a lagging group A, and are further equivalents to the OMIB [16].
The basic process of real-time equivalence is mainly divided into the following three steps: i.The generator power angle collected by the system is sorted from the largest to the smallest.ii.The difference between the angles of the adjacent two units in this sequence is calculated and the largest angle gap is selected as the splitting line.The group above the gap is leading group S and the one below the gap is lagging group A. iii.The leading group and the lagging group are further equivalents to the OMIB.
Taking into account the non-linearity of the actual power system, the direction is no longer unique after the state-plane trajectory crosses the boundary of the concavity-convexity.After entering the convex area, the equivalent state-plane trajectory may return to the concave area again due to the existence of nonautonomous factors.Therefore, if the angle stability of the system is identified only by whether the state-plane trajectory passes through the concavity-convexity demarcation line, a misjudgement may occur when the stability of a non-autonomous system with an equivalent value is identified.In [17], the authors proposed an auxiliary index, but it was pointed out in [18] that the method in [17] with incorrect clustering and the time-varying of equivalent single-machine system parameters would cause the method to

Transient energy supplement criteria
The generator state-plane trajectory is only a surface representation of transient instability, and the energy and its changes comprise the essence of power system-stability determination.In this section, we give energy supplemental criteria.Define the transient relative kinetic energy of the generator as follows: The first derivative of the transient relative kinetic energy is The second derivative of the relative kinetic energy is According to the theory in Sections 2.2 and 2.3, when the stateplane trajectory passes through the concave-convex boundary, it will always be in the convex region, thereby satisfying μ = Δω ⋅ l > 0, i.e. −Mcos δP 2max ω 0 Δω 2 > (P m − P 2max sin δ) 2 .Therefore, when the state-plane trajectory enters the convex region, the second derivative of the transient relative kinetic energy is always greater than zero; that is to say, the first derivative of the transient relative kinetic energy will remain monotonously increasing.Therefore, there must be a moment t 0 that satisfies the following formula: Equation (15) indicates that the power angle difference of the power system will be larger, and the system will be unstable.
For OMIB, the relative kinetic energy represented in ( 15) is used as a supplementary criterion and is combined with the stateplane trajectory method to form a comprehensive criterion.When applied to a multi-machine system, the clusters are first equated according to Section 2.3 and then the comprehensive criteria are used for judgments.

Transient instability judgment process
In this article, we assume that the system-wide unit has real-time measurement and communication functions.Subsequently, the comprehensive real-time transient stability judgment process of the multi-machine system based on state-plane trajectory and transient energy change only needs to collect data at two-time points.The system can be monitored by scrolling through the read WAMS data.The specific process is as follows.The stability criteria flowchart is shown in Fig. 5.
For real power systems, the real-time data on main generators are obtained based on PMU/WAMS and includes angles, angular velocity and electromagnetic power.The two indexes in Fig. 5 can be calculated for every synchronous sample.When they are all positive at the same time, the system is considered to be transiently unstable.Real-time detection is accomplished as follows: 1. Collect the data on the system in real time.If the system has a large disturbance, go to Step 2. 2. According to the data of each generator at the current moment, including power angle, angular speed, electromagnetic power and mechanical power, group the generators and make them equivalent to the OMIB.3. Calculate the state-plane trajectory concave-convexity criterion according to the equivalent system.If μ > 0, go to Step 4. Otherwise, go to Step 5. 4. Calculate the relative kinetic energy change rate of the generator according to the equivalent system.If V ke ′ > 0, the system can be judged as unstable.Otherwise, go to Step 5. 5. Determine whether the power angle of the generator reaches the transient instability equilibrium point according to the equivalent system.If yes, judge the system as unstable.
Otherwise, enter Step 6. 6. Wait for the next time that the data are read and return to Step 2.

Simulation analysis
In order to verify the integrated real-time transient stability identification method proposed in this paper, the IEEE 10-machine 39-bus system was taken as an example.The bonneville power administration (BPA) program was used to simulate various faults occurring at different locations and the same faults lasting for different time periods.The calculation step was taken as 0.01 s, the generator was a third-order model and the load was a constant impedance model.Parameters such as the angle, mechanical power, angular velocity and electromagnetic power of the generator were obtained through simulation calculations.Since the BPA program cannot directly verify the correctness of the criteria, it was necessary to write another MATLAB program to analyse and calculate the data.

Critical stability after fault calculations
Fault 1 was set as a three-phase short circuit grounding between buses 4 and 14.After 0.2 s, the fault was cut off and the system was in continuous oscillation.Taking generator 30 as a reference, the relative generators angles changed over time, as shown in Fig. 6.
Using the group equivalent theory in Section 2.3, the equivalent generator angle is shown in Fig. 7.
Through the generator angle, it can be seen that the system was critically stable.The change rate index of the transient relative kinetic energy over time is shown in Fig. 8.When the fault occurred, μ was >0 at 0.62 and 0.98 s, and the state-plane trajectory reached the convex region at these two moments in time.Only relying on the state-plane trajectory concave-convexity criterion would have led to a misjudgement.According to the transient energy supplementation criterion, the transient relative kinetic energy change rate is a necessary condition for system instability.As shown in Fig. 9, at 0.62 and 0.98 s, the transient relative kinetic energy change rate was negative, so the system instability could not be judged at these two moments in time.
The criterion of combining the state-plane trajectory and the transient relative kinetic energy change rate could not satisfy the condition of all positive values at the same time.Therefore, the system was determined as stable, which was consistent with the actual situation.

Example of instability after failure
Fault 2 was set as a three-phase short-circuit grounding between buses 16 and 19.After 0.15 s, the fault was cut off.Each generator angle changed over time, as shown in Fig. 10.Generators 34 and 35 were advanced units, whereas the others were lagging ones.
Using the group equivalent theory in Section 2.3, the equivalent generator angle is shown in Fig. 11.
The equivalent generator angle monotonically increased, so the system was unstable under this fault.
The state-plane trajectory criterion asperity is shown in Fig. 12, where we can see that at 0.3 s, μ was greater than 0 and the system would have had signs of instability.Nevertheless, we needed to analyse the transient energy criteria further.
At 0.42 s, the change rate of the relative kinetic energy was positive and the state-plane trajectory reached the convex region.Therefore, system instability was determined at 0.42 s (Fig. 13).
For currently used engineering methods, it has been determined that the time of instability is 0.76 or 1.23 s if the maximum generator angle before or after the equivalent is used exceeds 180 ∘ as a criterion for instability, respectively.Both methods have established a destabilisation time greater than the method proposed in this paper.

Other fault conditions when changing the operation mode
In order to verify the necessity and accuracy of the proposed method further, we changed the operation mode of the system and simulated more faults.The calculation step was 0.01 s, the generator was a second-order model and there were 4 buses with 80% constant impedance and a 20% constant current load, 7 buses with 70% constant impedance and a 30% constant current load and 15 buses with a 50% constant impedance, a 30% constant current load and a 20% constant power load.The specific situations are shown in Table 1.
Fault 3 is a single-phase broken line provided between buses 17 and 18, and the circuit was cut off at a different cutting time T c .
Fault 4 is a single-phase short-circuit grounding provided between buses 16 and 17, and the cut-off time T c cuts the circuit.
It can be seen from the results that if only the state-plane trajectory concavity-convexity identification method is used, it might misjudge near the critical cutting time.The method proposed in this paper could accurately distinguish the transient stability and had a shorter recognition time than the currently used methods.

Conclusions
In this paper, a comprehensive real-time transient stability identification method based on state-plane trajectory and relative kinetic energy change rate, which solves the problem of judging the stability by trajectory concavity and convexity only and avoids the problem of misjudgement in the condition of critical stability, is proposed.
When analysing the relationship between the minimum value of the trajectory on the state-plane and the stability of the system in a single-machine infinite-bus system, if there is a minimum value of the state-plane trajectory at an unstable equilibrium point, the system is unstable.On this basis, the relationship between the state-plane trajectory concavity-convexity and the system stability was analysed.We conclude that if the state-plane trajectory enters the convex region from the concave region, the system becomes unstable.
The relationship between the state-plane trajectory passing through the concave-convex boundary and the relative kinetic energy change rate was analysed.The transient relative kinetic energy change rate is used as a supplemental criterion for temporary stability identification and is combined with the stateplane trajectory method to form an integrated real-time transient stability criterion, which improves the accuracy of the identification.
Taking the IEEE10-machine 39-bus system as an example, a simulation analysis model using BPA and MATLAB was established to verify the effectiveness of the proposed algorithm.

Fig. 2
Fig. 2 State-plane trajectory at different fault clearing times

Fig. 5
Fig. 5 Stability criteria based on the state-plane trajectory and the transient relative energy change rate

Fig. 13
Fig. 13 Change rate index of the transient relative kinetic energy in fault 2

Table 1
Comparison of the system transient stability identification schemes Fault T c , s System simulation State-plane trajectory concavity-convexity method Proposed method Currently used methods