Stability of Power Grids: State-of-the-art and Future Trends

The state of the art of transient stability and steady-state (small signal) stability in power grids are reviewed. Transient stability concepts are illustrated with simple examples; in particular, we consider two machine learning-based methods for computing region of attraction: ROA produced by Neural Network Lyapunov Function; estimation of the ROA of IEEE 39-bus system using Gaussian process and Converse Lyapunov function. We discuss steady state stability in power systems, and using Prony’s modal analysis for evaluating small signal stability for the 7 Bus Test system and real French power system.


Introduction
Power grids are large, highly complex, ever-changing structures with an objective to deliver electric power reliably and economically from generators to loads (Wood et al., 2012). A critical design consideration for a transmission system is that of maintaining the "stability" of the grid, which is maintaining the balance between the consumption and generation (Willems et al., 1971).
Depending upon the nature of disturbance in power system, the stability problems are classified into two categories: steady state stability and transient stability. Although stability of a power system is essentially a single problem, the various forms of instabilities in the system cannot be properly understood and effectively dealt with by treating it as the single problem. Analysis of stability, including identifying key factors that contribute to instability and devising methods of improving stable operation, is greatly facilitated by classification of stability into appropriate categories. The following classification of the power system stability problem is suggested in (Kundur et al., 2004): 1. Rotor Angle Stability (short-term stability)  Small-disturbance or small-signal rotor angle stability  Large-disturbance rotor angle stability or transient stability 2. Voltage Stability (short-term or long-term stability)  Small-disturbance voltage stability  Large-disturbance voltage stability 3. Frequency Stability (short-term or long-term stability). The objective of this paper is to overview recent progress in transient and steady-state stability of power grids. We discuss transient stability in Section 2. In particular, brief overview of two methods for transient stability analysis, time-domain and direct methods, is presented in Sections 2.1, 2.2. Section 2.3 addresses the problem of computing regions of attraction (ROA). In Section 2.4 we present two examples; three methods are described in more details: ROA produced by Neural Network Lyapunov Function, ROA produced by sum of squares optimization (SOS) method and ROA produced by linear quadratic regulator (LQR) for wind turbine dynamics evaluation; estimation of the ROA of IEEE 39 bus system using Gaussian process and Converse Lyapunov function. In Sect. 3 we discuss steady state stability in power systems, and present an example of using Prony's modal analysis for evaluating small signal stability for the 7 Bus Test system and real French power system.

Transient stability 2.1 Problem Statement
The transient stability of a power system is its ability to maintain synchronous operation of the machines when subjected to a large disturbance. The occurrence of such a disturbance may result in large excursions of the system machine rotor angles; as a consequence, loss of synchronism may result among machines, which, as a rule, develops in very few seconds after the disturbance inception. The nonlinear character of transient stability, its fast evolution, and its disastrous practical implications make it one of the most important and at the same time most problematic issues to assess and even more to control (Gajduk et al.,2014).

Methods
Modern analysis of the stability of electric power systems is based, as a rule, on the numerical integration of nonlinear models of these systems. In the process of integration, the system trajectories are calculated for different fault scenarios, various operating conditions and different network topologies. Special software packages are used to simulate power systems models in the off-line (not real-time) mode. For example, EUROSTAG and RUSTAB software packages are used by the System Operator of United Power Systems of Russia for planning and managing electrical regimes. However, time-domain methods do not provide stability margins, which would give indication "how far" from instability the system is, which would yield suitable sensitivity analysis tools (Ribbens-Pavella et al., 2000). At the same time, for 60 years already, an alternative approach, called "direct methods", is used to control power systems. These methods are based on the use of the so-called energy functions (or energy functionals), and generalize the Lyapunov function method (Chiang, 2011;Chiang, 2013). The main advantage of this approach is the consideration of non-linear and non-stationary characteristics of the power system model. This allows one to calculate power system states after large disturbances that occur during short circuits on power transmission lines or when high power generators and consumers are disconnected from power system (Machowski et al., 2012).
The main idea of direct methods is in the numerical construction of energy functions, which in structure are similar to the Lyapunov functions, but they may not be such. The problem is to numerically estimate the function in the range of admissible values which is also called the region of attraction (ROA). The existence of the Lyapunov function guarantees the stability of a nonlinear system. In the case of direct methods application, the stability must be proved by numerical simulation of the critical power system states in which instability can occur (Frank et al., 2012).

Region of attraction and its computation
The ROA of a general dynamical system normally refers to a region where each state can converge to the stable equilibrium point as time goes to infinity. Likewise, the ROA of a power system is viewed as a set of operating states such as rotor angles and frequencies able to converge to the stable equilibrium, which corresponds to the solution of power flow problem (Mehta et al., 2016), after being subject to a disturbance. The corresponding convergent trajectory is regarded as a stable state trajectory. The estimation of ROA is basically dependent on the construction of Lyapunov function and its level set.
A Lyapunov function by definition is a continuously differentiable locally positive definite function ( ) whose time derivative is locally negative semi definite with the exception of the equilibrium point where the derivative is 0 (see Figure 1) (Zhai et al., 2019). Without loss of generality we shift the equilibrium point to the origin so that = 0. The red ellipse in Figure 1 describes the level set { ∈ 2 | ( ) ≤ 1 } , while the blue one denotes the level set { ∈ 2 | ( ) ≤ 2 } with 1 < 2 and = ( 1 , 2 ) . The black dot represents a stable equilibrium point of dynamical system. Each state in the level set can converge to the equilibrium point on condition of ( )/ < 0.
Finding a good Lyapunov function and calculating that maximizes the region in is not a trivial task. A lot of research effort is devoted to this problem and many proposals have been made. A Lyapunov function can be identified efficiently via a semi-definite program (SDP, (Boyd et al., 2009)) when the dynamics are polynomial and the Lyapunov function is restricted to be a SOS polynomial (Parrilo, 2000). Other methods to compute ROAs include maximization of a measure of ROA volume over system trajectories (Henrion et al., 2014), and sampling-based approaches that generalize information about stability at discrete points to a continuous region (Bobiti et al., 2016).  Most existing works on estimating ROA rely on analytical Lyapunov functions, which are subject to two limitations: the analytic Lyapunov functions may not be always readily available, and the resulting ROA may be overly conservative. A method proposed in (Zhai et al., 2019) overcomes these two limitations by leveraging the converse Lyapunov theorem in control theory to eliminate the need for an analytic Lyapunov function and learning the unknown Lyapunov function with the Gaussian Process (GP) approach. By treating the unknown Lyapunov function by a GP, the authors propose to assess state stability of power systems with the aid of the converse Lyapunov function. For an existing Lyapunov function, such an approach allows for extending the preexisting ROA with a provable confidence level.

Estimation of the ROA of the wind turbine using Lyapunov neural network algorithm
The case example deals with the novel converter less wind power generator with a control framework that consists of an excitation synchronous generator, permanent magnet synchronous servo motor, signal sensors, and servo control system. The wind and servo motor powers are integrated with each other and transmitted to the excitation synchronous generator via a coaxial configuration as proposed in (Jyothi et al., 2015).
A method for stabilization of the inverted pendulumbased servo motor control system was developed in (Tomin et al., 2019). The presented control model can be used to improve the wind power generation using the novel converter from (Yunhai et al., 2010) (see Figure 2). The servo motor controller plays an important role in output power and grid voltage phase tracking. A situation in which the controller detects a power increase from the servo motor implies decreasing wind speeds. At this moment, the system regulates the exciter current to reduce the excitation generator output power. A chain reaction subsequently occurs in which the servo motor power returns to a balanced level. During the energy balance periods, the servo motor consumes only a slight amount of energy to stabilize the shaft speed.  Figure 3 shows the diagram of the system produced by a single-rotation-fulcrum inverted pendulum. Using servo motor control system as a rotary actuators, rotary arm is connected with the motor shaft through the slip ring, one end of the installation is disk displacement detection sensor, an inverted pendulum is connected through the bearing and optical code disc which is fixed on the rotating arm.
The inverted pendulum is governed by the differential equation 2 2 2̈2 = 2 2 sin 2 −̇2 + , with state ≔ ( 2 ,̇2), where 2 is the rotation angle for the rotating arm from the upright equilibrium 0 = 0, is the input torque, 2 is the pendulum mass, is the gravitational acceleration, 2 is the pole length, and is the friction coefficient. We discretize the dynamics with a time step of Δ = 0.01 and enforce a saturation constraint ∈ [− , ̅ ̅], such that the pendulum falls over and cannot recover when it goes beyond a certain angle. In this section, we present details for the implementation of Lyapunov Neural Network (NN) algorithm proposed in (Richards et al., 2019) to learn the largest safe (stability) region of a simulated inverted pendulum wind turbine. We also provide experimental results in a comparison to other methods of computing Lyapunov functions. The paper (Richards et al., 2019) proposed the Lyapunov candidate ( ) = θ ( ) ( ) to learn the requisite features, where θ : ℝ → ℝ is a feed-forward neural network with parameter vector θ . The property of neural networks to approximating any continuous function was exploited together with gradientbased parameter training to closely match the true ROA with a level set of the candidate without the need to engineer individual features of . Let and be a feedback policy and a closed-loop dynamics, respectively. Then, we can always use the one-step decrease condition as a provable safety certificate to identify safe level sets that are subsets of the largest safe region .
Thus, we specify a training algorithm that adapts the candidate to the shape of the dynamical system's trajectories via classification of states as safe or unsafe. For simplicity, we use the perceptron loss, which penalizes misclassifications more when they occur far from the decision boundary. The training algorithm was designed to adapt the parameters θ such that the resulting Lyapunov candidate satisfies Δ ( ) < 0 throughout We implement a Lyapunov NN algorithm on the inverted pendulum benchmark with the Python code, which uses an open-source software library TensorFlow. For the neural network Lyapunov candidate , we used three layers of 64 Tanh activation units each. We prescribe ( ) with = 1 as the level set that delineates the safe region . Figure 4 shows the results of training with Lyapunov NN algorithm, and the largest safe level set ( 18 ) with 10 stochastic gradient descent (SGD) iterations per update. The true largest ROA is estimated by forward-simulating all of the states in a state space discretization, and set volume is estimated by counting discrete states. Figure 4 also shows the largest safe sets for a LQR Lyapunov candidate and a SOS Lyapunov candidate. The LQR candidate ( ) = is computed in closed-form for the same discretized, linearized, unconstrained form of the dynamics used to determine the LQR policy   can certify more of the true largest safe region than other common Lyapunov candidates. This has important implications for safe exploration during learning for complex dynamical systems such as power grids.

Estimation of the ROA of IEEE 39-bus system using Gaussian process and Converse Lyapunov function
In this case study, we investigate the problem of estimating the ROA for power systems using GP approach to assess state stability of power systems with the aid of the converse Lyapunov function.
The GP approach for learning the unknown Lyapunov function (i.e., GP regression) was proposed in (Zhai et al.,  2019). Normally, a general GP regression requires a prior distribution of unknown functions specified by a mean function, a covariance function, and the probability of the observations and sampling data to obtain the posterior distribution. Without loss of generality, we consider an unknown function ℎ( ) as a GP, which can be sequentially measured by ( ) = ℎ( ) + , ∈ + , where ( ) refers to the observed function value for the input at the -th sampling step, and the measurement noise is zero-mean, independent and bounded by . With the GP approach, we can obtain the posterior distribution over ℎ( ) by using sampling data in the training set {( 1 , 1 ), ( 2 , 2 ), … , ( , )}.
The unknown Lyapunov function ( ) can be approximated by a GP. The covariance or kernel function ( , 0 ) encodes the smoothness property of ( ) from the GP. If the prior distribution over ( ) is unknown, ( ) can be approximated by the linear combination of kernel functions in the reproducing kernel Hilbert space (RKHS), i.e. ( ) = ∑ ( , ( ) ) . In brief, the function ( ) gets smoother as ‖ ‖ decreases.
In addition, a Gaussian Process Upper Confidence Bound (GP-UCB) based sampling algorithm is designed to reconcile the trade-off between the exploitation for enlarging the ROA and the exploration for reducing the uncertainty of sampling region. For a given value of the parameter, the specified confidence level ∈ (0,1), and the sampling domain ∈ ℝ , which is a subset of the state space, our goal is to maximize the region of attraction, wherein each point converges to the origin with the probability of at least 1−δ.
Numerical simulations were conducted on the IEEE 39-bus 10-machine system ( Figure 5). Bus 31 with a generator is assigned as the swing bus. We then considered the dynamics for the remaining 9 machines. The detailed parameter setting and the algorithm can be found in MATLAB code made available in GitHub (GitHub). Figure 6 shows the assessment result for the ROA of IEEE 39 bus system with 9 machines. In each panel of this figure, the state points inside the red dashed ellipse is guaranteed to converge to the origin, while those in the yellow region are asymptotically convergent with the probability of at least 95%.  By Monte Carlo like estimate, the volume of yellow region is about 2.3×10 4 times larger than that of the certified ROA. This achievement is due to the large number of state dimensions or the size of the considered dynamical system. Actually, the profile of yellow regions largely depends on the distribution and size of sampling points as well as the choice of kernel functions for GP learning.

Problem Statement
Steady state stability analysis is the study of the power system stability when the system is subjected to small perturbations such as incremental changes in system load. The time frame of interest is in the order of several seconds. These disturbances are considered to be sufficiently small so that the linearization of the system equations is permissible for the purpose of analysis. There is extensive literature related to linearized models in power systems (Kundur et al., 1990;Wang et al., 1989). They are obtained numerically from the nonlinear differential algebraic models. Linear models are not only useful for small signal stability analysis of power systems, but also for developing of power system control techniques.
Small signal stability is concerned with the system response to small disturbances, where the system is operating in a linear region. Inter-area modes are characterized by oscillations between groups of generators, or groups of plants, separated by long geographical and electrical distances. As a result, the rotor speeds of synchronous generators can oscillate slowly in response to a disturbance in electrical power. For generators separated by long electrical distances these oscillations tend to be out of phase. The frequency of these electromechanical oscillations, typically in the range of 0.1 to 1.0 Hz, is a function of the inertia and damping in each area as well as the impedance of the transmission system.

Methods
Currently, one of the main tool for the small signal stability calculation is the modal analysis (Wang et al., 2016;Gibbard et al., 2015). The spectral methods are among the most popular methods for estimating the power system stability and for diagnosis of power system technical condition. These methods are based on the calculation of the spectrum of the linearized power system models. The terms "selective modal analysis" and "participation factors" were introduced in (Verhese et al., 1982; Perez-Arriaga et al., 1990) for the analysis of largescale power systems. Participation factors allow to analyze the relationship between the state variables and the eigenvalues of the linearized dynamic model. They were also interpreted in terms of stochastic averaging, mode energy, as well as their observability, controllability and mobility (Abed et al., 2000;Tawalbeh, 2010). In a broader sense, selective modal analysis can be understood as a set of methods for the fast calculation of eigenvalues and eigenvectors in the selected matrix spectrum. The effective methods for isolating critical modes are the QR method, the inverse iteration method, the Lanczoc method (Wang et al., 1989), modified Arnoldi method (Arnoldi, 1951;Kundur et al., 1990), calculation of the spectrum of dominant poles (Martins, 1997) and the method of matrix signum functions (Misrikhanov et al., 2008). The RSE Research Center (Milan) recently developed a new software package for selective modal analysis (Gaglioti, 2014). An alternative approach to conventional modal analysis was proposed in (Yadykin, 2010;Yadykin et al., 2014). In these papers, the solution of the Lyapunov equations was represented as a sum of Hermitian matrices corresponding either to individual eigenvalues of the dynamics matrix or to their pairwise combinations. Each eigen-part was called a sub-Gramian. The closest to the method of sub-Gramians is the approach, in which the spectral properties of the solutions of the Lyapunov equations were used to reduce the dimension of large dynamical systems. In particular, in (Antoulas, 2005) the spectral decompositions of the square H2-norm of the transfer function into its poles were obtained. Further, these decompositions were applied to obtain the approximating dynamic power system model of reduced dimension. Model order reduction techniques are evaluated in connection with their approximation accuracy, computational reliability and efficiency, as well as their ability to preserve the dynamic qualities of the original system (Baur et al., 2014). Some of the most popular methods for reducing the dimensionality of models are based on the analysis of Gramians, in particular, methods of balanced truncation (Moore, 1981) (when the system is transformed to a basis in which the controllability and observability Gramians have the same diagonal form) and the method of using cross-Gramians (Fernando, 1984).
Another approach for performing small signal analysis of large scale power systems is nonlinear simulation, followed by analysis of input-output data. Typical stimuli for the system include resistive brake insertions and modulation of power injection at different nodes in the system. One advantage of this approach is that its results can be compared to synchrophasor measurement unit data of real-world disturbances for model validation (Byrne et al., 2016). Once the data has been obtained from the simulation, common analysis techniques include Eigensystem Realization Algorithm (ERA), Prony analysis, and the Matrix pencil algorithm. The ERA algorithm is described in (Juang et al., 1985)

Case studies
Small-signal stability analysis is the ability of the power system to maintain synchronism under small disturbances. Such disturbances occur continually in the system because of small variations in loads and generation. Instability mat arise in two forms: a) an increase of the rotor angle due to lack of sufficient synchronization torque, or b) the rotor oscillations of increasing amplitude due to lack of sufficient damping torque. The small-disturbance stability problem is equivalent to ensure sufficient damping in system oscillations (Prabha, 1994). This index is computed from an estimate of the eigenvalues of the system, which are determined using time series from dynamic simulations. The eigenvalues allow to calculate the damping ratio of each system oscillation to provide a measure of stability for a given contingency. Estimating the damping of the system time series is a challenging problem. First of all, the appropriate signal processing techniques have to be applied.
The modal analysis problem may be posed as follows: given a set of measurements that vary with time, it is desired to fit a time-varying waveform of pre-specified form (such as (2)) to the actual waveform. Each state variable of the system is given by: where the parameter is the residue (?) of the mode , 0 is obtained from initial conditions, and represents the (possibly complex) eigenvalue of (что такое А?). The estimation of these responses yields modal information about the system that can be used to obtain the prediction of possible unstable behavior, controller design, parametric summaries for damping studies, and modal interaction information. Prony's method is designed to directly estimate the parameters for the exponential terms in (1), by fitting a function to an observed record for ( ) . The Prony method is a "polynomial" method in a sense that it includes the process of finding the roots of a characteristic polynomial. In the first test the 7-bus test system from (Glower et al., 2012) is used to demonstrate the Prony's modal analysis. The diagram of the test system is shown in Fig.  7. The results are shown in Figure 8 and the reconstructed signals are marked with dots. Table 1  The next example presents the validation of the smallsignal stability index using simulations from the transmission system operator of France (RTE) for French power system (RTE, 2012). Figure 9 depicts the simulations results after applying the small-signal index to the simulations of the French network. The active power in lines (118 lines) has been analysed as shown in Figure 10. Only one mode with a frequency between 0.1 and 1 Hz has been identified. This mode has an estimated damping ratio of 1 =1.74% and a frequency of 1 = 1.8 , it is presumably a local mode. After applying Prony's approach, information about the signals such as the identified modes is available and is used to compute the three-layer small-signal stability index. This index provides a measure related to the stability of the system and is based on the damping ratio of the system's modes. The measure is the angular distance (in radians) from each mode to a pre-defined damping ratio. • Single Mode Index (SMI), a matrix. • All Modes Index (AMI), a vector. • Global Modes Index (GMI), a scalar. The SMI is the most detailed source of information, but for a system with a large number of modes, it is not easy to interpret the whole matrix of data. It provides the individual distance of each mode to a pre-defined damping ratio, e.g. 0 = 0%, 5 = 5%, and 10 = 10%,. If any of the elements of SMI is negative, this indicates that the mode corresponding to the specific row has a damping ratio less than the damping ratio for that required column. AMI facilitates the interpretation of SMI, it gives the minimum distance of the modes with respect to each of the damping ratios. If one value of the AMI row is negative, this means that at least one mode has a damping less than the one required. GMI gives a global interpretation of the modes with respect to all pre-defined damping ratios. It provides the minimum distance among all modes respect to all pre-defined ratios. In Table 2, it can be seen that at least one criterion was violated for French power system, as indicated by the negative number in GMI. At least one mode has broken the 5 = 5% and 10 = 10% damping requirement, as indicated by the negative numbers in the second and third column of AMI. Mode 1 (M1) has damping ratio below 5 = 5% and 10 = 10% and greater than 0 = 0% . Since only one mode was identified, SMI is a vector rather than matrix and is equal to AMI (SMI=AMI).

Conclusions
In this paper we have reviewed the recent progress in transient and steady-state stability of power grids. We have considered two methods for transient stability analysis, time-domain and direct method.
In direct methods one computes region of attraction; two learning-based methods for their computations are described in more details: Lyapunov Neural Network, and Gaussian process and Converse Lyapunov function. We have also discussed steady-state stability in power systems; an example of using Prony's modal analysis for evaluating small signal stability for the 7 Bus Test system and real French power system has been presented.