A PMU-based algorithm of synchronous generator stability prediction during a disturbance

A PMU-based algorithm of synchronous generator stability prediction is proposed and tested in the paper. The algorithm is based on the application of the equal area criterion for the following spaces: active power – load angle, active power – time. The deceleration area is formed using the least squares curve fitting by the second order polynomial. The algorithm was tested on the single synchronous machine power system. A number of transients was simulated with the power angle curve polynomial accuracy of 50 ms. The results have shown high adaptability and efficiency of the proposed algorithm. It can be used in adaptive protection and control devices to form control actions aimed at ensuring the stability of a power system operation.


Introduction
Currently, there are significant changes both in the principles of power systems (PS) control and in the composition of theirs equipment, namely one that is generating and consuming electricity. The PS changes are based on digital transformation and the introduction of new mechanisms for the operation of the electricity and capacity markets. The introduction of a significant number of renewable energy sources (RES) and energy storage systems (ESS) [1], as well as technology of the demand response, lead to an increase in irregular components of active power flows. In addition, the introduction of low-carbon generation leads to a decrease in the total inertia of a PS and an increase in the requirements for the response speed of emergency control to ensure transient stability (TS).
Traditionally, the problem of analyzing TS of a complex PS is solved on the basis of equivalent circuits for the PS under protection, which are mathematical models, the parameters of which are most often set according to rated data or data tests and may greatly differ from the actual ones [2]. The irrelevant, rigidly specified parameters of PS models lead to a decrease in the accuracy of the TS analysis.
In the changing conditions of PS operation, increased requirements for speed and accuracy are imposed on the emergency control systems. The purpose of the paper is to develop and test an algorithm for predicting a TS based on the equal area criterion in the time domain and in the area of load angles of a synchronous generator (SG).

Review of methods for transient stability analysis
Historically, the method of numerical integration was one of the first methods for analyzing TS [3], the trajectories of changes in load angles of the SG or changes in their velocities are analyzed in order to determine the TS region.
Almost simultaneously with the methods of numerical integration for the analysis of TS the method of energy ratios analysis began to be applied, later called the equal areas method [4], which is used now [5][6][7][8]. The first application of the equal area method for the analysis of the complex PS TS was described in [9], 10 years later -in [10], and then it was considered as a special case of Lyapunov's theory [11].
As computer technology developed, the pattern recognition began to be applied to the analysis of TS [12]. However, after a certain time it became clear that pattern recognition methods, also known as artificial intelligence (AI) methods, at the time of publication of [12] could not be used for calculations in real time due to significant requirements for computing power.
For modern PSs, the analysis of TS is a nonlinear problem of large dimension. The following methods are applied for its solution [13]: • Numerical integration of the differential equations system that describes the dynamic model of a PS, • Lyapunov's Direct method [14], • AI methods [15], • Equal area method [16][17][18]. Numerical integration methods have found wide application in problems of analysis, planning and emergency control of PSs. This class of methods implies a step-by-step acquisition of the trajectory of change in the PS operation parameters for a given disturbance on the basis of a previously prepared mathematical model. For the TS analysis, two time intervals are consideredperiods of emergency and post-emergency mode. From the point of view of applying the methods of numerical integration for the analysis of TS in real time, there are the following difficulties: • It is necessary to perform a numerical solution of the differential equations system for a significant time interval to determine the TS region, • The methods are based on using the models of PSs, parameters of which can deviate from the real ones.
In connection with the need to analyze TS on the electromechanical transient time scale, the direct Lyapunov method was used [14]. It, in contrast to the methods of numerical integration, does not allow describing the trajectory of power flow parameters change. In [13], the following difficulties of applying the direct Lyapunov method are highlighted: • Ambiguity in forming a Lyapunov function, • The complexity of determining the TS region due to the presence of several stable equilibrium points, • In general case it is impossible to obtain all TS regions.
The Lyapunov's method of TS analysis finds application in hybrid methods that combine the direct Lyapunov method and that of numerical integration [13].
It was proposed to use AI methods for a qualitative analysis of dynamic stability in [15]. These methods are aimed at identifying dependencies in the data of the training sample and translating them to a new data set. The following advantages of using AI methods for the problem of control actions selection are distinguished [13]: • High performance of a correctly trained algorithm, • No necessity to use models of PS elements, • The ability to provide a predetermined accuracy. Application of the AI algorithms to solve differential equations in real time has the following setbacks [13]: • A significant difficulty of forming the training dataset that is sufficient to train an algorithm, • An algorithm must be re-trained if a PS structure was changed either as the result of a power system development or changes of generation capacity.
At present, the equal area method can be considered as one of the main real-time methods of TS analysis [13]. The whole SG group can presumably be divided into two subgroups in case of an instability. In general terms the algorithm of TS analysis, which is based on the equal area method, can be represented as such [13]: • Clustering of SGs into two equal subgroups, • TS analysis based on the criterion of system energy integral.
Most often, the SG clustering is implemented through out-of-step detection. Some authors propose other clustering methods, such as ones based on the energy criterion [19] or monitoring for correlation of SG load angle changes [20]. The single (or two-) machine equivalent is formed after SG division. Then TS is estimated using the equal area method -energies of acceleration and deceleration, namely those of rotor and turbine.
The problem of TS analysis and consequent selection of a control action (CA) on a transient time scale have high response time requirements both in terms of software and hardware. The total time delay is about 210-350 ms. The realtime TS analysis methods are mainly limited by power flow calculation times, which can be reduced up to 5-10 ms by using fast algorithms [21].

The algorithm of synchronous generator stability prediction
The purpose of TS analysis is a synthesis of the emergency control law, which requires the value of energy difference in dimension [MW•s]. The algorithm of a SG TS prediction during a electromechanical transient is based on application of the equal area criterion [22], which can be carried out for two spaces: active power -load angle (P-δ) and active power -time (P-t).
First version allows to estimate values of deceleration and acceleration areas, which are proportional to changes in energies of rotor deceleration and acceleration respectively. Applying this method in the space P-t allows to directly calculate these mentioned values of energy changes.
The block-diagram of the SG TS prediction algorithm during a transient is shown on the Fig.1.
The first stage of the algorithm is a searching of an event inception point by measuring instantaneous currents and voltages [23]. Further, during identification of an event recovery point, the calculation of the acceleration energy is made. This energy is spent on increasing the kinetic energy of rotating SG and turbine masses during a fault. The acceleration energy value is calculated in the space P-t. The deceleration area is estimated for the space P-δ with the transition into the space P-t for a given model of the SG power angle curve. The deceleration area is the one located between curves of the SG active power and turbine power. It is limited by the angular positions of the SG rotor, which correspond to the end of a fault and unstable equilibrium [22]. Predicting the deceleration area in the space P-δ allows to use a simple model for changing the SG active power -power angle curve. At the last stage of the algorithm, the difference between the values of the acceleration energy and the predicted value of the deceleration energy is calculated. A positive energy difference (∆W) means the TS loss and a negative one means stable state.
The following dataset is required for the algorithm of SG TS prediction during an electromechanical transient: • Instantaneous values of currents and voltages at the terminals of the SG stator winding (i, u) are measured using a phasor measurement unit (PMU), • The angular speed of the SG rotor (Ω) measured using a tachometer sensor installed on the SG rotor, • SG load angle (δ) is determined using a load angle sensor installed on the SG rotor (for example, a tachometer sensor) or by means of calculation, • The moment of inertia of the SG rotor (J) can be determined from the SG rated data, during SG testing or operation, for example, using an adaptive model of a synchronous machine [24].
• Energy spent on damping the oscillations of rotationg masses of SG and turbine (Wd).

The acceleration energy estimation
Then the transient time is the time interval between fault inception and termination. Assuming that change in the energy supplied from the turbine during a fault, ∆Wm = 0, (1) can be represented as follows: Hence, difference between the kinetic energy and that of a SG output equals to energy loss due to damping.
According to the Newton's second law the kinetic energy of rotor and turbine can be found as follows: where J is inertia moment of a SG rotor, Ω(t) is angular velocity of a SG rotor at the time t. Due to inertia the load angle cannot change instantly, therefore a fault in proximity to a SG is followed by decrease of its active power output. Therefore, angular velocity of rotor increases with constant mechanical power. The acceleration energy of a SG within the fault time is: ( )

Wm Wе
Wd Wl Wk where ∆Ω is change of rotor angular velocity at the time of fault termination, Ω0 is steady-state angular velocity of SG rotor before disturbance.
The energy change ∆We can be found by integrating difference of SC active power and that of turbine: where tstart and tfinish are respectively times of energy estimation start and finish during a fault, Pm(t) is mechanical power of turbine, Pe(t) is active power of SG. Area that is proportional to the acceleration energy can be found by integrating excess power of turbine over space P-δ: where δstart and δfinish are respectively SG load angles related to start and finish of energy estimation Pm(δ) is active power of turbine, Pe(δ) is active power of SG. Area that is proportional to the deceleration energy can be found by integrating excess power of turbine over space P-δ: where δstart and δfinish are respectively SG load angles related to start and finish of energy estimation, Pm(δ) is active power of turbine, Pe(δ) is active power of SG. Based on (2) the energy loss due to damping can be found as follows:

The deceleration energy prediction
One of the following models of curve fitting for active power and load angle can be used to predict an area that is proportional to the deceleration energy of non salient pole SG (Adec): where Pmax is maximal power value of power angle curve.
• Polynomial fit to operation of the AVR-PSS excitation control system (ECS).
Using of this polynomial allows taking into account influence of the ECS operation on a SG stability. On the other hand, consideration of the AVR-PSS ECS implies a priori setting a number of system parameters for each disturbance under consideration, which is impossible on a transient time scale. Equation (9) or (10) can be used to make the procedure of post-contingency power angle curve fitting an adaptive one. The fast algorithm can be used in order to find model coefficients (10) thereby reducing response time [21].
The deceleration energy (Wdec) is estimated using the calculated area value Adec and the procedure of transition between spaces "active power -time" and "active power -load angle".

Transition between spaces "active powertime" and "active power -load angle"
The transition between spaces P-t and P-δ can be carried out by means of the following expressions: where W (P-t) is energy of the space P-t, W (P-δ) is energy of the space P-δ, ∆w(t) is deviation of the SG magnetic field velocity from the synchronous one, Pm(δ) and Pm(t) is active power of turbine, Pe(δ) and Pe(t) is active power of SG, tstart and tfinish respectively times of energy estimation start and finish, δstart and δfinish are respectively SG load angles related to start and finish of energy estimation.
Projections of the curve δ(P,t) on planes P-t, P-δ and δ-t are shown on the Fig.3. The grey plane, which is parallel to the plane δ-t, represents the turbine power. The deceleration energy is coloured orange, while the energy of the space P-δ is coloured blue. The process of transformations (11) and (12) is marked with a red arrow.
The results of checking the expressions (11) and (12) are shown on the Fig.4 and Table 2. It was conducted for the single machine power system [25] with 400 MW being value of SG active power. Headings of plots on the Fig.4 correspond to the column "№" of the Table 2. Three three-phase faults were simulated as disturbances for one of two parallel transmission lines, their durations are 0.14, 0.25 and 0.27 s.  Results of energy change ∆We calculation for the case of fault (№ 1) and post-contingency state (№2, 3) are shown in the Table 2. The calculation is carried out using (5). Also the Table includes the following: area proportional to the acceleration energy Aa (№ 1) is calculated using (6) and that proportional to the deceleration energy Adec (№ 2, 3) using (7); results of transforming these energies to the space P-t using (10). Areas Aa and Adec are coloured grey on the Fig.4. For all cases under consideration the difference between ∆We and W (P-t) is within the range from -0.96 to 0.06 MW•s, which is from -0.64 to 0.04% in relation to ∆We.

Case study
A number of numerical experiments was carried out using Matlab Simulink software in order to test the proposed algorithm. The accepted assumptions are: • SG is a non-salient-pole one, • Load angle of SG is found by means of a shaftangle encoder [26].

Estimation of the acceleration energy
The following transient was simulated in order to estimate energies of acceleration and deceleration: right after a fault at time of 5 s an out-of-step occurs, with duration of 0.27 s. SG parameters, which correspond to acceleration time of the SG, are shown on the Fig.6. Different fragments of a single transient are shown on Fig.6 and Fig.7. ∆Wk is found using (4) Fig. 6. SG operating parameters.
Aa (P-t) means transformation of the acceleration area belonging to the space P-δ to the acceleration energy of the space P-t using (11). Difference between W (P-t) and ∆We is 0,15 MW•s, which is 0,26% in relation to ∆We.  Post-contingency power angle curve prediction was carried out by means of model (10) for training times from 25 to 100 ms. ∆Ω (change in the SG rotor angular speed) prediction was obtained by model (10) the coefficients of which were given for training times from 25 to 100 ms as well. ∆Wd prediction was obtained by means of (8).

Deceleration energy prediction
Results of power angle curve prediction and prediction error distribution of deceleration area Adec are shown in Fig. 8 and Fig. 9.    Fig. 9 shows the prediction error distribution in deceleration area with regard to the reference value, when the training time is varied from 25 to 100 ms. The reference value of the braking area was obtained by expression (7). The minimum prediction error in power angle curve is achieved at a training time of 70 ms.
By means of equation (8), energy loss due to damping the oscillations of the SG rotor during the deceleration can be predicted. To do this, it is necessary to predict changes in the angular velocity of the SG rotor ∆Ω and the value of the SG active power Pe(t). What is more, Pe(t) can be predicted using the power angle diagram that was obtained to predict the deceleration area Adec. Prediction of changes in the SG rotor angular velocity is carried out using a quadratic polynomial. Fig.  10 shows the results of ∆Ω prediction.    Taking the energy spent on damping during the deceleration of the SG rotor into account, difference between the acceleration energy and the predicted deceleration energy for the considered transient is determined as follows: Adec (P-t) means the transformation of the deceleration area from the P-δ space to the acceleration energy in the P-t space using the equation (11). A negative value of ∆W indicates the stability loss of the SG, which is confirmed by the calculation of the electromechanical transient.  Table 4 and Table 5. For case №4 in Sections 4.2 and 4.3, a detailed calculation was given. The headings of the plots (1) -(4) in Fig. 12 correspond to the "№" column in Table 4.

Synchronous generator stability prediction
Loss of stability is shown in the figures as a periodic process, accompanied by a change in the sign of active power when a synchronous machine switches from a generator mode to a motor mode [22]. Sequential operation of a synchronous machine in a generator and motor mode represents a cycle of out-of-step condition.
For the considered transients, the training time 70 ms, based on Fig. 9 and Fig. 11.
"Stability loss (calculation)" column in Table 5 corresponds to the result of the TS analysis by means of the electromechanical transient calculation. "Stability loss (prediction)" column in Table 5 corresponds to the TS prediction using the proposed algorithm.  For the considered transients, the results given by proposed algorithm coincide with the results of electromechanical transients calculations.
The algorithm allowing to predict TS according to PMU data is proposed and tested in the paper. It is based on equal area criterion both in the load angle and the time domains. In order to carry out the transition between these domains the equations (11) and (12) are used. The algorithm uses measurements of instantaneous currents and voltages, the load angle of the SG and the mechanical angular velocity. The inertia constant of the SG rotor is used as a initial data, which can be updated by using the SG adaptive model [25].
As a result of numerical simulation in Matlab Simulink software, the following results were obtained: • For the mathematical model under consideration, the minimum of the predict error for deceleration area corresponds to the training time of 70 ms; • When the training interval value reaches 40 ms, the ∆Wd values differ by less than 0.8 MW•s; • For the four considered numerical experiments of transients, the TS prediction was performed on the training time of 70 ms.
The future work will be aimed at testing the proposed algorithm on the base of physical modeling data. In addition, further investigation is needed for its testing as the main algorithm in the emergency control in real time.
The reported study was funded by RFBR, project number 20-38-90140.