A method to extract instantaneous features of low frequency oscillation based on trajectory section eigenvalues

Abstract

Affected by the nonlinear time-varying factors due to fault scenarios, protection relaying, and control measures, the dynamic behaviors of a power system may be significantly different from the results of previous methods. In order to analyze the oscillation characteristics of complex power systems more accurately and suppress the low frequency oscillation more effectively, this paper improves the trajectory section eigenvalue method. Firstly, the time response of a system is obtained by numerical simulation in a given fault scenario. Secondly, the algebraic variables are substituted to the differential equations along the trajectory. Thus, the original time-varying differential-algebraic equations are approximated by a set of linear ordinary differential equations, which can be updated along the trajectory. On this basis, this paper proposes a method to extract instantaneous features of the oscillation from the micro perspective. The non-equilibrium points with strong nonlinearity or critical eigenmodes are identified by the proposed method. The simulation test results of the IEEE 3-machine 9-bus system and the New England system illustrate the validity of the proposed method.

1 Introduction

In recent years, the low frequency oscillation (LFO) has become increasingly intricate, with the enlargement of power grids [1], the integration of renewable energy and the development of high voltage direct-current transmission [2, 3]. The universality and accuracy of previous analysis approaches for small disturbance stability can hardly satisfy the requirements of modern power systems. Influenced by the essential nonlinear factors and time-varying factors due to fault scenarios, protection relaying, and control measures, the results obtained by traditional methods may be significantly different with the dynamic behaviors of a power system [4]. Hence, it is urgent to put forward an LFO analysis method that can describe the dynamic behaviors of a power system more accurately with a complete consideration of the nonlinear time-varying factors.

Generally, LFOs can be divided into three forms based on the causes: ① steady increase in generator rotor angle due to lack of synchronizing torque; ② rotor oscillation with increasing amplitude due to lack of sufficient damping torque [5]; ③ sustained oscillation due to periodic disturbance resulting from failures or malfunctioning equipment [6]. The eigen-analysis method is most widely used in practice. It describes the dynamic system behaviors near the equilibrium point by a linear time-invariant system. However, the nonlinearity of a power system may lead to bifurcation and blow-up of solutions [7, 8]. Reference [9] addresses the nonlinear singularity phenomenon by Hopf bifurcation theory, which demonstrates that the dynamic behaviors of a simple power system could be different from the results of eigen-analysis. However, limited by the heavy computation in large-scale systems, this method is still under investigation as well in methodological aspects as in concrete applications. In order to consider the effects of nonlinear factors, the normal form theory [10] and the modal series method [11] are proposed. More detailed information of eigenmodes is extracted by solving the high-order approximate equations. One major drawback of these approaches is that the solutions are obtained at the equilibrium point, neglecting time-varying factors. For the sake of considering the effects of random factors, probability models are used to improve the universality of traditional methods [12]. All of these methods reviewed thus far, however, suffer from the limitation of disturbance scenarios. As a result, some critical factors of a power system such as the non-differentiable components and control measures are ignored [13].

To entirely analyze the influence of nonlinear time-varying factors on power system oscillation, it is quite necessary to extract information from trajectories. Some signal analysis techniques are used to identify the oscillation modals. Common approaches like time-frequency transformations take the disturbed trajectory as research object. The frequency features are identified by decoupling frequency components, and the damping features are extracted by differentiating mode amplitudes with respect to time. As a classical method, the Prony analysis replaces the trajectory by a set of exponential components [14]. This method can describe the dynamic behaviors in an interval approximately, though the influences of time-varying factors are averaged. Based on wide-area measurement signals and the proper orthogonal decomposition (POD) method, [15] uses dynamic mode decoupling (DMD) method to achieve oscillation modes decoupling and features extracting. Nevertheless, these methods can hardly describe the instantaneous characteristics, which are time-varying. The empirical mode decoupling (EMD) method extracts instantaneous features by decomposing the original signal into locally symmetric components on the basis of envelope characteristics and the Hilbert transform [16]. However, the precision of EMD method is strongly affected by the frequencies of adjacent components. To extract instantaneous features more accurately, window Prony method [17], and window fast Fourier transformation (ridges) method [18] use slide windows to divide the entire time to a set of overlapping periods along the trajectory. The feature sequence is composed of the modals of all windows. Despite of that, the accuracy of the methods is limited by the width of the window. Moreover, the truncated error caused by the edges on both sides of the signal is inevitable. To improve the adaptive ability of the time-frequency transformation approaches, wavelet (ridge) method is introduced to assist analyzing LFOs [19]. Researches show that the method can identify the modals of an actual trajectory with noises effectively [20], although the choice of wavelet basis relies on empirical parameters. Reference [21] approximates the observed process as a linear system for more information of models. The stability margins are obtained as posterior probabilities that the poles of the estimated system are unstable.

A considerable number of researches have been done to achieve the extraction of instantaneous features by signal analysis techniques during the last decades. To be objective, these approaches are instrumental in identifying modals online. However, all these mentioned methods suffer from a serious limitation. The essence of signal analysis techniques is identification of the average characteristics in a time window. The time scale of the oscillation analysis is limited by the resolution of the methods. On the other hand, the major goal of modal identification of LFO is to suppress the oscillation. To locate critical generators or parameters, the signal analysis techniques have to be unified with the models at the equilibrium point, which can hardly reflect the real dynamic behaviors of a nonlinear time-varying system.

Attempts to resolve this dilemma have resulted in the invention of trajectory section eigenvalue (TSE) theory, which extends the concept of eigenvalue from the equilibrium point to other non-equilibrium points of the trajectory. The pioneer work can be traced to [22]. Reference [23] attempts to use TSE method to analyze transient stability. The method is further developed to estimate the energy of a system before the visual far end point or the visual dynamic saddle point [24]. Reference [25] concludes that the global stability of linear time-varying systems cannot be determined by TSEs. The major reason is that the global dynamics of time-varying systems are affected by the kinetic energy at non-equilibrium points, which could not be considered by TSE method. However, [26] has proved that the effects of the kinetic energy can be ignored when analyzing the instantaneous characteristics in one integration step. Meanwhile, the mathematical and physical foundations of TSE method have been given in this literature.

The contributions of this paper are threefold. Firstly, it improves the mathematics of TSE method to give a better explanation of the physical meaning of the eigenvalues at non-equilibrium points. Based on the improved theory, this paper proposes a method to extract the instantaneous features of LFOs. Secondly, the paper proposes the criterion to identify the non-equilibrium points with strong nonlinearity, and illustrates the mechanism of nonlinear oscillations in power system. Thirdly, it proposes the criterion to identify the dangerous non-equilibrium points and the critical eigenmodes. Compared with previous methods, the critical generators and parameters identified at above non-equilibrium points can achieve better inhibition effects.

2 Mathematics of trajectory section eigenvalues

2.1 Piecewise-linearized system models and solution in time domain

The behaviors of a dynamic system, such as a power system, may be described by a set of n first-order nonlinear differential-algebraic equations (DAEs) as the following form (the derivatives of the state variables are not explicit functions of time) [5]:

$$ \left\{ {\begin{array}{*{20}c}{\dot{\varvec{X}}} = \varvec{f (X ,Y )}\\ {\mathbf{0}} = \varvec{g (X ,Y )}\end{array} } \right.$$

(1)

The column vector X is referred as the state vector, and its entries ${x}_{i}$ as the state variables, such as rotor angles and angular velocities of the rotors. The column vector Y is referred as the algebraic vector, and its entries ${y}_{{i}}$ as the algebraic variables, such as node voltages and node current injections. The solutions of X and Y in the time domain can be obtained by numerical simulation in given fault scenarios. By substituting Y into the differential equations ${\dot{\varvec{X}}} = \varvec{f} (\varvec{X}, \varvec{Y} )$ and linearizing the nonlinear components of the functions along the trajectory, the DAEs are replaced by a set of linear ordinary differential equations (ODEs). For an illustrative purpose, the linear ODE at ${t}_{{k}}$ is taken as an example to introduce TSE method.

$$\Delta \dot{\varvec{X}}_{k} (\Delta {t} )=\varvec{A}_{k}\Delta \varvec{X}_{k} (\Delta {t} )+ \varvec{B}_{k}$$

(2)

where $\varvec{A}_{k} = \left. {\frac{\partial \varvec{f}}{\partial \varvec{X}}} \right|_{{\varvec{X} ({t}_{{k}} )}}$ is the Jacobi matrix of the nonlinear function f; $\varvec{B_{k}} = \varvec{f}\left( {\varvec{X} ({t}_{{k}} ), \varvec{Y} ({t}_{{k}} )} \right)$ is the initial values of the right side of the differential equation at ${t}_{{k}}$. Supposing that U_k is the right modal matrix of A_k. New state vector is defined by $\varvec{Z}_{k} (\Delta {t} )= \varvec{U}_{k}^{ - 1}\Delta \varvec{X}_{k} (\Delta {t} )$. Both sides of (2) are left-multiplied by $\varvec{U}_{k}^{ - 1}$:

$$\dot{\varvec{Z}}_{{k} } (\Delta {t} )=\varvec{ U}_{k}^{ - 1} \varvec{A}_{{k} } \varvec{U}_{{k} } Z_{{k} } (\Delta {t} )+\varvec{ U}_{{k} }^{ - 1} \varvec{B}_{{k} }$$

(3)

Solve (3) with the initial values of the elements of $\varvec{Z}_{{k} } (\Delta {t} )$ respect to $\Delta {t}$ ($\varvec{Z}_{{k} } (0 )= {\mathbf{0}}$), which means the incremental time in $\left[ {{t}_{{k}} ,t_{{{k + }1}} } \right)$.

$$\varvec{Z}_{{k} } (\Delta {t} )=\varvec{ C}_{{k} } \varvec{E}_{{k} } (\Delta {t} )-\varvec{\varLambda}_{{k} }^{ - 1} \varvec{U}_{{k} }^{ - 1} \varvec{B}_{{k} }$$

(4)

where $\varvec{E}_{{k} } (\Delta {t} )= \left[ {{\text{e}}^{{{\lambda }_{{{k,}1}}\Delta {t}}} , {\text{e}}^{{{\lambda }_{{{k,}2}}\Delta {t}}} ,\ldots , {\text{e}}^{{{\lambda }_{{{k,n}}}\Delta {t}}} } \right]^{\text{T}}$; $\varvec{\varLambda}_{{k} }= \varvec{ U}_{{k} }^{ - 1} \varvec{A}_{{k} } \cdot\varvec{U}_{{k} } = {\text{diag (}}{\lambda }_{{{k,}1}} ,{\lambda }_{{{k,}2}} ,\cdots ,{\lambda }_{{{k,n}}} )$ is the eigenvalues of the state matrix A_k at t_k. Two procedure parameter matrices D_k and C_k are defined as follows:

$$\varvec{D}_{{k} } = \left[ {\begin{array}{*{20}c} {{d}_{{{k} , 1}} } \\ {{d}_{{{k} , 2}} } \\ \vdots \\ {{d}_{{{k} ,{{n}} }} } \\ \end{array} } \right] = - \left[ \begin{aligned} \frac{1}{{{\lambda}_{{{k},1}} }}\left( {{{v}}_{{{{k}} ,11}}{ {f}}_{1} + {{v}}_{{{{k}} ,12}} {{f}}_{2} + \cdots + {{v}}_{{{{k}} ,1n}} {{f}}_{{n} } } \right) \hfill \\ \begin{array}{*{20}c} {\frac{1}{{{\lambda }_{{{k},2}} }}\left( {{{v}}_{{{{k}} ,21}} {{f}}_{1} + {{v}}_{{{{k}} ,22}} {{f}}_{2} + \cdots + {{v}}_{{{{k}} ,2n}} {{f}}_{n} } \right)} \\ \vdots \\ \end{array} \hfill \\ \frac{1}{{{\lambda}_{{{k,n}}} }}\left( {{{v}}_{{{{k}} ,n1}} {{f}}_{1} + {{v}}_{{{{k}} ,n2}} {{f}}_{2} + \cdots + {{v}}_{{{{k}} ,n{n}}} {{f}}_{n} } \right) \hfill \\ \end{aligned} \right]$$

(5)

$$\varvec{C}_{{k} } = - \left[ {\begin{array}{*{20}c} {{d}_{{{k},1}} } & 0 & \cdots & 0 \\ 0 & {{d}_{{{k},2}} } & \cdots & 0 \\ \vdots & \vdots & {} & \vdots \\ 0 & 0 & \cdots & {{d}_{{{k},n}} } \\ \end{array} } \right]$$

(6)

where d_k,i is the i^th element of D_k; ${v}_{{{k,ij}}}$ is the element in the i^th row and the j^th column of $\varvec{U}_{{k} }^{ - 1}$; f_i is the i^th element of B_k. Both sides of (4) are left-multiplied by U_k:

$$\Delta \varvec{X}_{{k} } (\Delta {t} )= \varvec{U}_{{k} } \varvec{C}_{{k} } \varvec{E}_{{k} } (\Delta {t} )- \varvec{A}_{{k} }^{ - 1} \varvec{B}_{{k} }$$

(7)

Obviously, (7) describes the dynamic behaviors of the changes of state variables in $[t_{{k} } ,t_{{{k} + 1}} )$. The values of $\varvec{U}_{{k} } \varvec{C}_{{k} }$ indicate the degrees of eigenmodes excitation [27].

2.2 Physical meaning of trajectory section eigenvalues

The changes of the state variables respect to $\Delta {t}$ have been obtained by (7) in the time domain. Considering that the initial values of Z_k(Δt) are zero, we have $\varvec{Z}_{{k} } (0 )=\varvec{ C}_{{k} } \varvec{E}_{{k} } (0 )-\varvec{\varLambda}_{{k} }^{ - 1} \varvec{U}_{{k} }^{ - 1} \varvec{B}_{{k} } \varvec = {\mathbf{0}}$. Equation (7) can be rearranged to:

$$\Delta \varvec{X}_{{k} } (\Delta {t} )= \varvec{U}_{{k} } \varvec{C}_{{k} } \varvec{E}_{{k} } (\Delta {t} )- \varvec{U}_{{k} } \varvec{C}_{{k} } \varvec{E}_{{k} } (0 )$$

(8)

Assuming the existence of a new system, ${t} = {t}_{{k}} +\Delta {t}$, the state vector of this system X_k,new(t) satisfies:

$$\varvec{X}_{{{k} ,new}} ({t} )= \varvec{U}_{{k} } \varvec{C}_{{k} } \varvec{E}_{{k} } ({t} - {t}_{{k} } )$$

(9)

Therefore, the changes of the new state vector in $[t_{k} ,t_{k + 1} )$ can be expressed by:

$$\begin{aligned}\Delta \varvec{X}_{k,new} (\Delta {t} )= &\, \varvec{X}_{k,new} ({t} )- \varvec{X}_{k,new} ({t}_{k} )\\ = &\, \varvec{U}_{{k} } \varvec{C}_{{k} } \varvec{E}_{{k} } ({t} - {t}_{{k} } )- \varvec{U}_{{k} } \varvec{C}_{{k} } \varvec{E}_{{k} } (0 )=\Delta \varvec{X}_{{k} } (\Delta {t} )\\ \end{aligned}$$

(10)

According to (10), the dynamic behavior of the new system is the same with the original system, ignoring the unbalanced power in $[t_{k} ,t_{k + 1} )$. The state equation of the new system can be derived from (9).

$$\begin{aligned} \dot{\varvec{X}}_{k,new} (t )= &\, \varvec{U}_{{k} } \varvec{C}_{{k} }\varvec{\varLambda}_{{k} } \varvec{E}_{{k} } (t - t_{{k} } )\\ = & \,\varvec{U}_{{k} } \varvec{C}_{{k} }\varvec{\varLambda}_{{k} } (\varvec{U}_{{k} } \varvec{C}_{{k} } )^{ - 1} \varvec{X}_{k,new} (t )\\ \end{aligned}$$

(11)

C_k and Λ_k are both diagonal matrices, and the expression in (11) must satisfy $\varvec{U}_{{k} } \varvec{C}_{{k} }\varvec{\varLambda}_{{k} } (\varvec{U}_{{k} } \varvec{C}_{{k} } )^{ - 1} =\varvec{ U}_{{k} }\varvec{\varLambda}_{{k} } \varvec{U}_{{k} }^{ - 1} = \varvec{A}_{{k} }$. Equation (11) can be rearranged as:

$$\dot{\varvec{X}}_{{{k} ,new}} (t )= \varvec{A}_{{k} } \varvec{X}_{{{k} ,new}} (t )$$

(12)

Thus, the original system without regard to unbalanced power is approximated by a linear system in $[t_{k} ,t_{k + 1} )$. Although X_k,new(t) may be different from X_k(t), the changes of their state vectors are the same in $[t_{k} ,t_{k + 1} )$. The section eigenvalues $\left( {{\lambda }_{{{k} , 1}} {,\lambda }_{{{k,} 2}} , \cdots {,\lambda }_{{{k,n}}} } \right)$ describe the dynamic behavior of the imaginary system in $[t_{k} ,t_{k + 1} )$.

Besides, the parameter matrices U_k and C_k are updated along the trajectory according to the results of numerical simulation. The effect of the unbalanced power is reflected in following time sections. Furthermore, the accumulate error is avoided by updating the parameters of (7) and (12) along the trajectory. Consequently, the effects of nonlinear factors and time-varying factors are considered by the section state matrix changing with time.

3 Extraction method of instantaneous features of LFO

To extract the instantaneous features of LFO, the eigenvalues of the section state matrix are calculated. In general, an n-order matrix has n eigenvalues, which may consist of couples of complex eigenmodes and real eigenmodes. A couple of complex eigenmodes represent a sinusoidal component, while a real eigenmode represents an exponential component. In this paper ${f}_{{{k} ,i,ins}}$ is used to represent the imaginary part of the i^th eigenvalue at t_k, divided by a constant 2π. ${\sigma }_{{{k} ,i,ins}}$ is used to represent the real part of the i^th eigenvalue at t_k. As for a complex root, ${f}_{{{k} ,i,ins}}$ represents the instantaneous frequency characteristic of the i^th eigenmode at t_k. ${\sigma }_{{{k} ,i,ins}}$ describes the trend of the amplitude of the i^th eigenvalue at t_k. As for a real root, ${f}_{{{k} ,i,ins}}$ is zero. ${\sigma }_{{{k} ,i,ins}}$ represents the rate of a trend change within one integration step.

3.1 Identification of points with strong nonlinearity

In an oscillation, a couple of complex eigenmodes may reversibly transform into two real eigenmodes occasionally. A one-machine infinite bus (OMIB) system is used to illuminate the mechanism of the variation of TSEs.

The structure of the OMIB system is shown in Fig. 1. V represents the voltage of the infinite bus. P_e is the input power. The summation of reactance is denoted as X_Σ. ${E'}$ and δ are referred as the module and angle of the transient electromotance after the transient reactance of the generator, respectively. The classical model is used to represent the generator. The dynamic equations of the generator can be found in Appendix A (A1). When a disturbance occurs, the TSEs can be obtained along the trajectory:

$${\lambda }_{{{k,}1,2}} = - \frac{{D}}{{2{M}}} \pm \frac{{\sqrt {{D}^{2} - 4{MK}_{{k}} } }}{{2{M}}}$$

(13)

where D is the damping torque coefficient of the generator; ${K}_{{k}} =(E'V{ \cos \delta }_{k} )/{X}_{{\varSigma }}$ represents the slope of the electromagnetic power-rotor angle curve at t_k, geometrically; ${\delta }_{{k}}$ is the value of the rotor angle at ${t}_{{k}}$ with a value range of $\left( {0,{\pi }} \right]$. As can be expected, if ${\delta }_{{k}} \in \left( {0,{ \arccos }\left( {{D}^{2} {X}_{{\varSigma }}/(4{ME'V} )} \right)} \right)$, ${\lambda }_{{{k,}1,2}}$ would be a couple of complex roots. With the increase of ${\delta }_{{k}}$, the imaginary parts of ${\lambda }_{{{k,}1,2}}$ decrease. If ${\delta }_{{k}} \in \left( {{ \arccos }\left( {{D}^{2} {X}_{{\varSigma }} {/} (4{ME'V} )} \right),\uppi} \right]$, ${\lambda }_{{{k,}1,2}}$ would be two different real roots. With the increase of ${\delta }_{{k}}$, one of them decreases, and the other increases. Therefore, the non-equilibrium points with strong nonlinearity can be identified briefly below: ① aggregating the complex eigenmodes at each time section; ② comparing the numbers of the complex eigenmodes at non-equilibrium points to the equilibrium point; ③ those non-equilibrium points are signed, the number of which is inconsistent with the equilibrium point.

Fig. 1

One-machine infinite bus system

The concept of the non-equilibrium points with strong nonlinearity can be used to explain some complex phenomena. For instance, the distortion of a rotor angle trajectory is induced near far end points (FEPs) as shown in Fig. 2a, c. The angular acceleration and velocity approach to zero when FEPs are close to the dynamic saddle point (DSP). Figure 2b depicts the evolution of the imaginary parts of the TSEs, the eigenmodes of which transfer to the exponential components from the oscillation components.

Fig. 2

Distortion phenomenon of rotor angle near FEPs

3.2 Identification of critical points

A two-machine system is presented in Fig. 3 to introduce the critical non-equilibrium points in general system.

Fig. 3

Two-machine system

The TSEs of this two-machine system after being disturbed can be of two forms: ① a couple of complex roots, a real root (can be zero), and a zero root; ② three real roots (can be zero) and a zero root. Equation (12) suggests that the positive signs of the real parts of TSEs represent the injection of dangerous energy. Thus, the critical non-equilibrium points are identified by the signs of the real parts of the eigenvalues along the trajectory.

For the sake of illustrating the relations between instantaneous features and global dynamics, Fig. 4 analyzes a typical convergent oscillation and a typical divergent oscillation in the two-machine system. As shown in Fig. 4a, in a convergent oscillation, the real parts of mode 3 could be positive periodically. As shown in Fig. 4b, in a divergent oscillation, the real parts of mode 1 and mode 2 are not always positive.

Fig. 4

Real parts of TSEs in a two-machine system

In the neighborhood of the equilibrium point, the TSEs are consistent with the eigenvalues at the equilibrium point. However, power system is an essential nonlinear dynamic system. Once the operating point of the system runs away from the equilibrium point, the dynamic behaviors may be different with the eigenvalues at the equilibrium point.

In a linear system, the accelerated power gradually increases with the decreasing relative angular velocity from DCP to FEP. However, affected by the nonlinearity of the P_e–δ curve, the accelerated power of a power system is lower than the linear system before the FEP. As a result, the damping characteristics near FEPs would be worse than the neighborhood of the equilibrium point. The closer the FEP is to the DSP, the lower the stability margin of the system suffers, and the trajectories of the real parts of TSEs change correspondingly. Figure 4a, b depicts the evolution of ${\sigma }_{{{k} ,i,ins}}$ in convergent and divergent oscillation, respectively. The values of ${\sigma }_{{{k} ,i,ins}}$ have obviously positive correlation with the damping characteristics.

Therefore, the trajectories of the real parts of TSEs indicate the instantaneous damping characteristics. Importance should be attached to the critical non-equilibrium points with positive real parts of TSEs.

To suppress power system oscillation, it is necessary to optimize the parameters of critical generators or controllers, which can be located by the eigenmodes with positive real parts. The participation factor and the electromechanical mode correlation ratio can also be extended from the equilibrium point to non-equilibrium points. Section 4 demonstrates the effectiveness of the instantaneous features in locating critical generator and guiding the optimization.

In addition, the relative positions of eigenvalues between adjacent time sections are random, which limits the application of TSE theory in identifying critical eigenmodes. For an illustrative purpose, assuming that the eigenvalues at ${t}_{{a}}$ and ${t}_{{b}}$ are ${\lambda }_{{{a},1}} {,\lambda }_{{{a,}2}} , \cdots {,\lambda }_{{{a,n}}}$ and ${\lambda }_{{{b,}1}} {,\lambda }_{{{b,}2}} , \cdots {,\lambda }_{{{b,n}}}$, respectively. ${\lambda }_{{{a,}1}}$ may relate to ${\lambda }_{{{b,}2}}$ or even ${\lambda }_{{{b,n}}}$. The eigen-polynomial is a continuous function of the parameters of state matrix, and the zero points of the eigen-polynomial are continuous functions of the parameters of eigen-polynomial. Consequently, the eigenvalues of a matrix are continuous functions of the parameters of the matrix. The TSEs between adjacent time sections can match with the similarity of their values.

3.3 Process of instantaneous features extraction

Figure 5 shows the procedure to extract instantaneous features. The method includes six steps as follows.

Fig. 5

Flow diagram of instantaneous features extraction

Step 1: the time-domain responses of state variables are obtained by numerical simulation or phase measurement units.

Step 2: the algebraic variables are substituted to the differential equations along the trajectory and the section state matrices are obtained by linearizing the nonlinear components of the differential equations. In this paper, the size of analysis step is equal to the integration step.

Step 3: the TSEs are calculated from the state matrix and the instantaneous features are extracted along the trajectory.

Step 4: the eigenvalues between adjacent time sections are matched.

Step 5: the points with strong nonlinearity are identified by corresponding criterion. Rⁿ and Iⁿ represent the fields of real number and imaginary number, respectively.

Step 6: critical points are identified by corresponding criterion.

4 Case study

4.1 IEEE 3-machine 9-bus system

In this section, a 3-machine 9-bus system is used to illustrate the limitation of the traditional eigen-analysis method.

For the structure shown in Fig. 6, the generators are all equipped with governors. The classical model is used to represent the generator, the transient electromotance ${E} '$ after ${X}_{d} '$ is constant. The dynamic equations of generator and governor are listed in Appendix A. The constant impedance model is used to represent the load. The line data, the generator data, and the governor data are listed in Appendix B.

Fig. 6

IEEE 3-machine 9-bus system

The eigenvalues at the equilibrium point are listed in Table 1 and the real parts of all eigenvalues are negative. In the view of the traditional eigen-analysis method, the trajectories of state variables should be convergent oscillation when the system is disturbed.

Table 1 Eigenvalues of IEEE 3-machine 9-bus system at equilibrium point

On the other side, the most dangerous eigenmode identified by the traditional eigen-analysis method should be the 3^rd and 4^th eigenvalues with the minimum real parts. Moreover, the corresponding cluster of generators is {G2, G3} vs. {G1}.

The oscillation characteristics of this system are analyzed in different fault scenarios. The numerical integration is used to assess the validity of the proposed method.

4.1.1 Scenario 1: weak nonlinear time-varying scenario

The disturbance is a three-phase transient fault in bus 6, and the fault duration is 0.09 s. The trajectories of the rotor angles are convergent oscillations with weak damping, as shown in Fig. 7.

Fig. 7

Trajectories of rotor angles of system in scenario 1

To demonstrate that the instantaneous features of the oscillation are consistent with the traditional eigen-analysis method in weak nonlinear time-varying scenarios, the instantaneous features are obtained as follows.

As shown in Fig. 8a, after the fault disappears (later than 0.09 s), in most of the time, there are four couples of complex eigenmodes, which represent four oscillation components with different frequency. According to the criterion for points with strong nonlinearity, there are no nonlinear points in the trajectories of ${\lambda }_{1,2}$, ${\lambda }_{3,4}$, and ${\lambda }_{5,6}$. Instead, the nonlinear points appear in the trajectories of ${\lambda }_{7,8}$, which are signed by purple blocks in Fig. 8c. The strong nonlinear eigenmodes and the dangerous eigenmodes may be different. As shown in Fig. 8b, the signs of the real parts of ${\lambda }_{7,8}$ are always negative, and the signs of the real parts of ${\lambda }_{3,4}$ are positive in some time sections, which are signed by red blocks in Fig. 8c. As can be expected, the critical points are related to the points with strong nonlinearity.

Fig. 8

Trajectories of TSEs in scenario 1

4.1.2 Scenario 2: strong nonlinear time-varying scenario

The disturbance is a three-phase transient fault in bus 6, and the fault duration is 0.15 s. The trajectories of the rotor angles are divergent oscillations, as shown in Fig. 9.

Fig. 9

Trajectories of rotor angles of system in scenario 2

The inconsistency between the trajectories of rotor angles and the eigenvalues at the equilibrium point reveals the limitation of traditional eigen-analysis method in strong nonlinear time-varying scenarios. To illustrate the mechanism of the oscillation from the micro point of view, the instantaneous features are obtained as follows.

As shown in Fig. 10a, b, after the fault disappears (later than 0.15 s), points with strong nonlinearity appear in the trajectories of ${\lambda }_{3,4}$, ${\lambda }_{5,6}$ and ${\lambda }_{7,8}$. The critical time sections appear in the trajectories of ${\lambda }_{3,4}$ and ${\lambda }_{5,6}$. The points with strong nonlinearity and critical eigenmodes are signed in Fig. 10c by purple and red blocks, respectively. Obviously, the most dangerous eigenmode is ${\lambda }_{5,6}$ with the maximum positive values in the oscillation, as shown in Fig. 10b. The cluster of the generators is the same as the actual instability mode: {G2}, {G1, G3}.

Fig. 10

Trajectories of TSEs in scenario 2

Therefore, the characteristics of instantaneous features can be summarized as follows: ① in weak nonlinear time-varying scenarios, the operating point of the system runs near the equilibrium point, and the instantaneous features are similar to the traditional eigen-analysis method; ② in strong nonlinear time-varying scenario, the operating point of the system runs far away from the equilibrium point periodically, and the instantaneous features may be totally different from the traditional eigen-analysis method.

4.2 IEEE 10-machine 39-bus New England system

In the view of the traditional eigen-analysis method, the eigenmodes with real parts greater than zero are dangerous. The eigenmodes with real parts approximate to zero are defined as weak damping eigenmodes. For the sake of suppressing LFO, these eigenmodes should be controlled. However, Section 4.1 has revealed that the traditional eigen-analysis method cannot describe the dynamic behaviors accurately. Moreover, the actually dangerous or weak damping eigenmodes may be ignored.

The New England system is used to illuminate the superiority of the proposed method in locating critical eigenmodes.

The structure of the system, the line data and the generator data can be found in [28]. These ten generators are all equipped with governors. The classical model is used to represent the generator, the transient electromotance ${E} '$ after ${X}_{d} '$ are constant. The dynamic equations of generators and governors are shown in Appendix A. The constant impedance model is used to represent the load. The governor data are listed in Appendix C. As a comparison, the eigenvalues at the equilibrium point are listed in Appendix D. It is obvious that the real parts of all eigenvalues are negative. The disturbance is a three-phase transient fault in bus 15 and the fault duration is 0.39 s. The trajectories of the rotor angles are shown in Fig. 11.

Fig. 11

Trajectories of rotor angles of IEEE 10-machine 39-bus system

Before the 3^rd second, the amplitudes of the rotor angles are decreasing and the dominant oscillation mode is {G39} to the set of the other generators. After the 3^rd second, the amplitudes of the rotor angles are increasing and finally the system loses stability. The critical mode is {G34} to the set of the other generators.

4.2.1 Real parts of trajectory section eigenvalues

Instantaneous features are extracted and the real parts of the TSEs are shown in Fig. 12. In order to highlight the critical eigenmodes, the other eigenmodes are hidden, the signs of the real parts of which are negative in all time. In intervals 1 and 2, the critical points mainly appear in eigenmode ${\lambda }_{29,30}$, the same with the eigenvalues at the equilibrium point. In intervals 3 to 6, the critical points mainly appear in the trajectories of a couple of new eigenmodes. Thus, the critical eigenmodes that mainly affect the dynamic behaviors of the oscillation should be the new couple.

Fig. 12

Trajecoties of real parts of critical eigenmodes in IEEE 10-machine 39-bus system

4.2.2 Location of critical generator

In order to locate the critical generator, it is necessary to calculate the participation factors of critical eigenmodes at critical time sections. The generator with highest participation factor in one interval is invariant. For an illustrative purpose, the results of eigenmodes analyses in typical time sections are shown in Table 2.

Table 2 Analyses of eigenmodes at typical time sections from intervals 1 to 6

In intervals 1 and 2, λ₂₉ and λ₃₀ are the most dominant eigenmodes with the highest participation factors and the critical generator is G39. In intervals 3 to 6, the new eigenmodes are the most dominant eigenmodes with the highest participation factor, and the critical generator is G34. The descriptions of the instantaneous features are the same as the actual trajectories of the rotor angles, as shown in Fig. 11. Considering that the electromechanical mode correlation ratios of these eigenmodes are much larger than 1, the state variables participating in the eigenmodes are mainly rotor angles and angular velocities.

To suppress the oscillation, the damping torque coefficient of G34 is up-regulated by 10%, and the numerical simulation shows that the trajectories of the rotor angles turn to be convergent, as shown in the left half of Fig. 13. By comparison, the damping torque coefficient of G39 is up-regulated by 50% and trajectories are still divergent after the optimization, as shown in the right half of Fig. 13.

Fig. 13

Trajectories of rotor angles of IEEE 10G39N system after parameter adjustments of G34 and G39, respectively

5 Conclusion

Based on TSEs, this paper proposes a method to extract instantaneous oscillation features from the micro perspective. To better illustrate the mechanism of power system oscillation and effectively suppress the LFO, the criteria for non-equilibrium points with strong nonlinearity and critical eigenmodes are proposed. At non-equilibrium points with strong nonlinearity, the dynamic behaviors may be totally different from the eigenvalues at the equilibrium point. The critical points rarely appear in weak nonlinear time-varying scenarios. Even though the oscillation is convergent, the real parts of TSEs may be positive at some non-equilibrium points with strong nonlinearity. The instability of oscillations may be caused by the accumulation effects of the critical non-equilibrium points, which massively appear in strong nonlinear time-varying scenarios. Moreover, simulation results clearly indicate the enhanced accuracy of critical generators location.

6 Discussion

Over the last 3 decades, the rotor angle stability has been characterized individually into two subcategories: small-disturbance stability and transient stability. As a result, it is impossible to analyze the relationship between oscillation stability and synchronous stability within the current framework.

The disturbances in some case studies of this paper belong to the traditional transient stability category, for three reasons:

• 1) It is difficult to distinguish small disturbance and large disturbance. In a negative damping oscillation defined by the traditional method, a complex dynamic process is inevitable before synchronous instability. Thus, this study begins with the aim of studying the mechanism and characteristics of the evolution of LFO, no matter the disturbance is large or small. Although TSE method cannot be used to analyze the global stability of a power system, simulation test results indicate a strong correlation between the oscillation characteristics and the proposed features.

• 2) In order to study the stability of a power system, a unified framework has been proposed by [29] to unify the analysis approaches of the oscillation stability and the transient stability in terms of energy. The improved TSE method and the achievements of this paper contribute to the unified framework. This paper is a tentative study on unifying small disturbance analysis and large disturbance analysis through TSE method.

• 3) The eigenvalues can only describe the dynamic behaviors in the neighboring region of one point, whether equilibrium point or non-equilibrium point it is. To analyze the oscillation characteristics in strong nonlinear time-varying scenarios, the state matrix and the variables are updated along the trajectory, according to the numerical results.

Finally, the proposed method would be used in engineering only when the traditional eigen-analysis cannot describe the dynamic behaviors accurately, or detailed analysis is required in some complex scenarios. Thus, compared with the traditional eigen-analysis, the increasing computation is acceptable. The quantity of information and computation can be balanced by determining the step size more intelligently. Moreover, mature methods have been used to reduce the dimension curse when calculating eigenvalues at the equilibrium point. The theorem of the inverse operator can be used to determine the circle where the eigenvalues lie. The radius of the circle can be estimated by Gerschgorin theory. Many results have been made in this field [30]. These achievements can also be used to estimate the distribution of the eigenvalues at the non-equilibrium points.

Appendices

Appendix A

$$\left\{ \begin{aligned} \frac{{{\text{d}}\delta }}{{{\text{d}}{t}}} = \omega \, \\ \frac{{{\text{d}}\omega }}{{{\text{d}}{t}}} = \frac{1}{{M} }\left( {P_{m} - P_{e} - D\omega } \right) \\ \end{aligned} \right.$$

(A1)

where δ represents the angle of the rotors, ω represents the angle velocity of the rotor; M represents the inertial of generator; P_m represents the mechanical torque; and P_e represents the electromagnetic torque. In the data lists of Appendixes B and C, the direct-axis transient reactance of generator is denoted by ${X}_{d} '$, the quadrature axis synchronous reactance is denoted by X_q, the damping torque coefficient is denoted by D.

$$\left\{ {\begin{array}{*{20}l} {\frac{{{\text{d}}{\mu }}}{{{\text{d}}{t}}}{ = \sigma } \, } \\ {\frac{{{\text{d}}{\sigma }}}{{{\text{d}}{t}}} = \frac{1}{{{T}_{{G}} {T}_{{p}} }} ({\mu } - {\xi } - {T}_{{G}} {\sigma } ) { }} \\ {\frac{{{\text{d[}}{\xi } - ({R + D}_{{d}} ){\mu } ]}}{{{\text{d}}{t}}} = \frac{1}{{{T}_{{d}} }} ({R\mu } - {\xi } )} \\ {\frac{{{\text{d(}}{P}_{{m}} + 2{\mu } )}}{{{\text{d}}{t}}} = \frac{2}{{{T}_{{w}} }} ({\mu } - {P}_{{m}} ) { }} \\ \end{array} } \right.$$

(A2)

where μ, σ, ξ represent the state variables of governors and water turbines, as shown in Fig. A1; K_p represents the gain of governor; R represents the regulation coefficient of governor; T_G and T_p represent the servo time constant and the pilot valve time constant of governor; T_d and D_d represents the soft feedback time constant and coefficient of governors; and T_w represents the water starting time of the turbine,

Fig. A1

Diagrams of governor and water turbine

Appendix B

The line data, generator data and governor data for IEEE 3-machine 9-bus system are given in Tables A1, A2 and A3, where V_limit represents the limit speeds of water gate.

Table A1 Line data for IEEE 3-machine 9-bus system

Table A2 Generator data for IEEE 3-machine 9-bus system

Table A3 Governor data for IEEE 3-machine 9-bus system

Appendix C

The governor data for IEEE 10-machine 39-bus system are given in Table A4.

Table A4 Governor data for IEEE 10-machine 39-bus system

Appendix D

The eigenvalues of IEEE 10-machine 39-bus system at equilibrium point are given in Table A5.

Table A5 Eigenvalues of IEEE 10-machine 39-bus system at equilibrium point