# An Accurate Third-Order Normal Form Approximation for Power System Nonlinear Analysis

Tian Tian, Xavier Kestelyn, Olivier Thomas, Hiroyuki Amano, Arturo Roman Messina

► To cite this version:

Tian Tian, Xavier Kestelyn, Olivier Thomas, Hiroyuki Amano, Arturo Roman Messina. An Accurate Third-Order Normal Form Approximation for Power System Nonlinear Analysis. IEEE Transactions on Power Systems, 2018, 33 (2), pp.2128-2139. 10.1109/TPWRS.2017.2737462. hal-02147596

HAL Id: hal-02147596

https://hal.science/hal-02147596

Submitted on 4 Jun 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Science Arts & Métiers (SAM)

is an open access repository that collects the work of Arts et Métiers ParisTech
researchers and makes it freely available over the web where possible.

This is an author-deposited version published in: https://sam.ensam.eu
Handle ID: http://hdl.handle.net/null

To cite this version :

Tian TIAN, Xavier KESTELYN, Olivier THOMAS, Hiroyuki AMANO, Arturo Roman MESSINA - An Accurate Third-Order Normal Form Approximation for Power System Nonlinear Analysis - IEEE Transactions on Power Systems - Vol. 33, n°2, p.2128-2139 - 2018

<figure><img src="image_2.png" /></figure>
# An Accurate Third-Order Normal Form Approximation for Power System Nonlinear Analysis

Tian Tian, *Student Member, IEEE*, Xavier Kestelyn ID, *Member, IEEE*, Olivier Thomas, Hiroyuki Amano, and Arturo Roman Messina, *Fellow, IEEE*

**Abstract**—The inclusion of higher-order terms in small-signal (modal) analysis has been an intensive research topic in nonlinear power system analysis. Inclusion of second-order terms with the method of normal forms (MNF) has been well developed and investigated, overcoming the linear conventional small-signal methods used in the power system control and stability analysis. However, application of the MNF has not yet been extended to include third-order terms in a mathematically accurate form to account for nonlinear dynamic stability and dynamic modal interactions. Due to the emergence of larger networks and long transmission line with high impedance, modern grids exhibit predominant nonlinear oscillations and existing tools have to be upgraded to cope with this new situation. In this paper, first, fundamentals of normal form theory along with a review of existing tools based on this theory is presented. Second, a new formulation of MNF based on a third-order transformation of the system's dynamic approximation is proposed and nonlinear indexes are proposed to make possible to give information on the contribution of nonlinearities to the system stability and on the presence of significant third-order modal interactions. The induced benefits of the proposed method are compared to those afforded by existing MNFs. Finally, the proposed method is applied to a standard test system, the IEEE 2-area 4-generator system, and results given by the conventional linear small signal and existing MNFs are compared to the proposed approach. The applicability of the proposed MNF to larger networks with more complex models has been evaluated on the New England–New York 16-machine 5-area system.

***Index Terms***—Interconnected power system, methods of normal forms, nonlinear modal interaction, power system dynamic, stability.

## I. INTRODUCTION

TODAY'S standard electrical grids are composed of several generators working in parallel to supply a common load.

Manuscript received January 31, 2017; revised June 23, 2017; accepted August 1, 2017. Date of publication August 11, 2017; date of current version February 16, 2018. Paper no. TPWRS-00106-2017. (Corresponding author: Xavier Kestelyn.)

T. Tian and X. Kestelyn are with University Lille, Centrale Lille, Arts et Metiers Paris Tech, HEI, EA, 2697-L2EP Laboratoire d'Electrotechnique et d'Electronique de Puissance, Lille F-59000, France (e-mail: Tian.Tian@ensam.eu; Xavier.Kestelyn@ensam.eu).

O. Thomas is with Arts et Metiers ParisTech LSIS UMR CNRS 7296 Laboratoire des Sciences de l'Information et des Systemes, Lille 59000, France (e-mail: Olivier.Thomas@ensam.eu).

H. Amano is with the System Engineering Research Laboratory, CRIEPI, Yokosuka 240-0196, Japan (e-mail: h-amano@criepi.denken.or.jp).

A. R. Messina is with the Department of Electrical Engineering, CINVESTAV, Guadalajara, 45070, Mexico (e-mail: aroman@gdl.cinvestav.mx).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TPWRS.2017.2737462

An important problem associated with interconnected power systems is the presence of oscillations that could have dangerous effects on the system. The multiplication of distributed generation units, usually composed of renewable-energy-based-generators, and the increase of energy exchanges through long distance lead to highly stressed power systems. Due to the large amount of power flowing through the lines, the low-frequency oscillations, called in classical power system studies electromechanical oscillations, exhibit predominant nonlinear behaviors. Since these oscillations are essentially caused by modal interactions between the system components after small or large disturbances, they are called nonlinear modal oscillations, higher order modes or higher order modal interactions, inaccurately modeled by the linear analysis based on a linearized model.

Although intensive research has been conducted on the analytical analysis of nonlinear modal oscillations based on the Normal Form Theory with inclusion of 2nd order terms in the system's dynamics, this paper proposes to show that in certain stressed conditions, as modern grids experience more and more, inclusion of 3rd order terms offer some indubitable advantages over existing methods.

The Method of Normal Forms (MNF) being based on successive transformations of increasing orders, the proposed 3rd order-based method inherits the benefits of the linear and the 2nd order-based methods, i.e.:

1) Analytical expressions of decoupled (or invariant) normal dynamics;

2) Physical insights keeping the use of modes to study the contribution of system components to inter-area oscillations;

3) Stability analysis based on the evaluation of the system parameters.

The paper is organized as follows. Section II introduces the need for including higher-order terms in the modal analysis of power systems. A literature review on the major applications of Normal Form Theory for the study of power grids is then conducted along with the proposal of a new formulation at the third order in Section III. Based on this new approximation, nonlinear indexes are proposed in Section IV to make possible to quantify the modal interactions and to give information on the effects of the nonlinearities on the system stability. Section V is dedicated to comparisons of the proposed third-order-based method to the existing methods in the literature. The different methods
reviewed and proposed in this work are applied to IEEE standard test systems (The Kundur's 4 machine 2 area system and the New England New York 16 machine 5 area system) to emphasize the advantages of taking into account high-order terms in the small-signal analysis. Section VI discusses the factors influencing the Normal Form analysis. Section VII proposes some conclusions and possible applications of the proposed method are suggested.

## II. NEED FOR INCLUSION OF HIGHER-ORDER TERMS

Small-signal analysis is the conventional analysis tool for studying electromechanical oscillations that appear in interconnected power systems. It linearizes the power system's equations around an operating point by including only the first-order terms of the Taylor's series expansion of the system's dynamic. The eigenanalysis is made to obtain analytical expressions of the system's dynamic performances and the stability analysis is realized on the basis of the first Lyapunov Method (Analysis of the real parts of the poles). Besides, modal analysis uses the eigenvectors to give an insight of the modal structure of a power system, showing how the components of the power system interact. Thanks to modal analysis, power system stabilizers can be placed at the optimal location in order to stabilize the whole system, ensuring then a small-signal stability [1].

Later, researchers suggested that in certain cases, such as when the system is severely stressed, linear analysis techniques might not provide an accurate picture of the power system modal characteristics. From 1996 to 2005, numerous papers [2]–[11] have been published proving that higher order modal interactions must be studied in case of certain stressed conditions. MNF with the inclusion of 2nd order terms shows its great potential in power system stability analysis and control design. Those achievements are well summarized in the Task-force committee report [12].

The existing 2nd-order-based method gives a better picture of the dynamic performance and the mode interactions than the classical linear modal analysis. However, it fails to take benefit of the system's nonlinearities for studying the stability where the conventional small-signal stability analysis fails. Based on this, [13] proposed to keep a second order transformation but with including some of the third-order terms in order to improve the system stability analysis.

Finally, excepted in [14] and [15], some 3rd-order-based MNFs have been proposed [16], [17] but have not been fully developed yet, not leading to a more useful tool than the ones using linear-based and second-order-based methods. For nonlinear mechanical systems, that often include lightly damped oscillatory modes and possible internal resonances, Normal Forms up to third order are widely used, either to classify the generic families of bifurcations in dynamical systems [18], [19] or to define Nonlinear Modes of vibration and to build reduced-order models [20]–[22].

## III. LITERATURE REVIEW ON THE EXISTING MNFS AND PROPOSAL OF A NEW METHOD FOR THE STUDY OF STRESSED POWER GRIDS

The Methods of Normal Forms (MNF) was initially developed by Poincaré [23] to simplify the system dynamics of nonlinear systems by successive use of near-identity changes

of coordinates. The transformations are chosen in such a way as to eliminate the nonresonant terms of a corresponding order.

The procedure is well documented in [24], [25] and can be adapted to the power system analysis. It consists of eight major steps:

1) Building the differential algebraic equations (DAEs) of the power system: differential equations and power flow constraints;

2) Solving the power flow to obtain the stable equilibrium point (SEP) for the post-fault system, i.e the operating point;

3) Transforming the DAEs on an equivalent system of differential equations and expanding the system of equations around the SEP into Taylor's series up to third-order;

4) Simplifying the linear part of the system by the use of a linear transformation;

5) Simplifying the non-resonant terms of higher-order terms by successive Normal Form (NF) transformations. This paper will use 2nd and 3rd order NF transformations;

6) Simplifying the Normal Forms' dynamics by neglecting (if possible) some resonant terms that can not be annihilated by NF transformations;

7) Reconstructing the original system's dynamic from the Normal Forms' dynamics in order to determine the order of the Taylor's series expansion and the NF transformations to be selected according to the expected accuracy;

8) Using the chosen Normal Forms' approximation for dynamic and stability analysis.

### A. Class of Systems that can be Studied by Methods of Normal Forms

The class of systems that can be studied by MNF are usually modeled using Differential Algebraic Equations (DAEs) [1]. By substituting the algebraic equations into the differential ones, one transforms those DAEs in a dynamical system, which can be written:

$$ \dot{x} = f(x, u), \quad (1) $$

where $x$ is the state-variables vector, $u$ is the system's inputs vector and $f$ is a nonlinear vector field. Expanding this system in Taylor series around a stable equilibrium point, one obtains:

$$ \Delta \dot{x} = H1(\Delta x) + \frac{1}{2!} H2(\Delta x) + \frac{1}{3!} H3(\Delta x) + \mathcal{O}(4) \quad (2) $$

where $Hq$ gathers the $q$-th-order partial derivatives of $f$, i.e., for $j = 1, 2, \dots, n$, $H1_k^j = \partial f_j / \partial x_k$, $H2_{kl}^j = [\partial^2 f_j / \partial x_k \partial x_l]$, $H3_{klm}^j = [\partial^3 f_j / \partial x_k \partial x_l \partial x_m]$ and $\mathcal{O}(4)$ are terms of order 4 and higher.

### B. Simplifying the Linear Terms

The linear part of (2) is simplified using its Jordan form:

$$ \dot{y} = \Lambda y + F2(y) + F3(y) + \mathcal{O}(4) \quad (3) $$
supposed here to be diagonal, where the jth equation of (3) is:

$$ \dot{y}_j = \lambda_j y_j + \sum_{k=1}^{n} \sum_{l=1}^{n} F2_{kl}^j y_k y_l + \sum_{p=1}^{n} \sum_{q=1}^{n} \sum_{r=1}^{n} F3_{pqr}^j y_p y_q y_r + \dots \quad (4) $$

$\lambda_j$ is the jth eigenvalue of matrix **H1** and $j = 1, 2, \dots, n$. **U** and **V** are the matrices collecting the right and left eigenvectors of **H1**. **UV** = **I**, $F2^j = \frac{1}{2} \sum_{i=1}^{n} v_{ji} [U^T H_2^i U]$ and $F3_{pqr}^j = \frac{1}{6} \sum_{i=1}^{n} v_{ji} \sum_{k=1}^{n} \sum_{l=1}^{n} \sum_{m=1}^{n} H_3^i_{klm} u_p^l u_q^m u_r^n$ where $v_{ji}$ is the element at j-th row and i-th column of **V** and $u_p^l$ the element at the p-th row and l-th column of **U**. The fundamentals of linear small-signal analysis use the sign of the real parts $\sigma_j$ of the eigenvalues $\lambda_j = \sigma_j + j\omega_j$ to estimate the system stability. **U** is used to indicate how each mode $y_j$ contribute to the state-variables of **x**, and **V** indicates how the state-variables of **x** are associated to each mode $y_j$.

### C. Simplification of a Maximum Number of Nonlinear Terms

Based on the Normal Form Theory, (3) can be simplified by successive NF transformations.

1) Elimination of Quadratic Terms: When only 2nd order terms are kept in the Taylor's Expansion Series of the system's dynamics, (3) becomes:

$$ \dot{y} = \Lambda y + F2(y) \quad (5) $$

where the jth equation of (5) writes:

$$ \dot{y}_j = \lambda_j y_j + \sum_{k=1}^{n} \sum_{l=1}^{n} F2_{kl}^j y_k y_l \quad (6) $$

To eliminate the second-order terms in (5), the following second-order NF transformation is applied:

$$ y = z + h2(z) \quad (7) $$

or, in component form:

$$ y_j = z_j + \sum_{k=1}^{n} \sum_{l=1}^{n} h2_{kl}^j z_k z_l \quad (8) $$

Use of (7) in (5) leads to:

$$ \dot{z} = \Lambda z - Dh2(z)\Lambda z + \Lambda h2(z) + F2(z) + \mathcal{O}(3) \quad (9) $$

where $Dh2$ is the Jacobian matrix of $h2$ and the jth term of (9) takes the form of:

$$ \dot{z}_j = \lambda_j z_j - \sum_{k=1}^{n} \sum_{l=1}^{n} [ (\lambda_k + \lambda_l - \lambda_j)h2_{kl}^j - F2_{kl}^j ] z_k z_l + \mathcal{O}(3) \quad (10) $$

If one wants to eliminate all second-order terms from (9), transformation (7) must satisfy the following equation, for all $z$:

$$ Dh2(z)\Lambda z - \Lambda h2(z) = F2(z) \quad (11) $$

which gives, for all $j, k, l$:

$$ (\lambda_k + \lambda_l - \lambda_j)h2_{kl}^j - F2_{kl}^j = 0 \quad (12) $$

Equation (12) shows that the computations of coefficients $h2_{kl}^j$ of the change of variables are not possible (or numerically difficult) if the conditions $\lambda_k + \lambda_l - \lambda_j \approx 0$ holds. This condition is only met for undamped systems (or weakly damped systems) and corresponds to the so-called internal resonances (when frequencies are such that $\omega_k + \omega_l - \omega_j = 0$). It is important to note that this condition appears only when the eigenvalues share a commensurability relationship, as for example when $\omega_1 = 2\omega_2$. Although it may be a rare case, [26] proved that it can be the reason of instability in the Japanese power system model.

If it is supposed that no internal resonance occurs, $h2_{kl}^j$ is given by [12]:

$$ h2_{kl}^j = \frac{F2_{kl}^j}{\lambda_k + \lambda_l - \lambda_j} \quad (13) $$

and, by neglecting terms of order greater than 2, the normal dynamics are modeled by a set of decoupled first-order linear differential equations:

$$ \dot{z} = \Lambda z \quad (14) $$

It is very important to note that the normal dynamics given by (14) uses the linear modes of the linearized system. The nonlinearities are then taken into account only through the nonlinear 2nd-order transform (7) and it results in a quadratic combination of linear modes, leading to 2nd-order modal interactions. Then, the stability analyses conducted using this linear normal dynamics give the same conclusion than the analysis that can be conducted using the linear dynamics.

Gradually established and advocated by investigators from Iowa State University, the period 1996 to 2001 [2]–[11] opens the era to apply Normal Forms analysis in studying nonlinear dynamics in power system. Its effectiveness has been shown in many examples [2], [5]–[8], [10], [11], [27], [28], to approximate the stability boundary [5], to investigate the strength of the interaction between oscillation modes [2], [6], [27], to deal with a control design [7], [8], to analyze a vulnerable region over parameter space and resonance conditions [10], [11] and to optimally place controllers [28].

This method is well summarized in [12], which demonstrated the importance of nonlinear modal interactions in the dynamic response of a power system and the utility of including 2nd order terms. Reference [29] proposes an extension of the method to deal with resonant cases. It has become a mature computational tool and other applications have been developed [30], [31], or are still emerging [32].

2) Consideration of Cubic Terms: Although the second-order-based method gives a more accurate picture of the system's dynamics than the linear small-signal analysis, keeping the 3rd order terms in the Taylor's expansion series adds some important features to the method. As examples, improvements in the stability analysis have been proved in [13] and new criteria to design the system controllers in order to improve the transfer capacity of transmissions lines have been established in [33].

Including third-order terms, (3) becomes:

$$ \dot{y} = \Lambda y + F2(y) + F3(y) \quad (15) $$
The j-th state equation being:

$$ \dot{y}_j = \lambda_j y_j + \sum_{k=1}^{n} \sum_{l=1}^{n} F 2_{kl}^j y_k y_l + \sum_{p=1}^{n} \sum_{q=1}^{n} \sum_{r=1}^{n} F 3_{pqr}^j y_p y_q y_r \quad (16) $$

Applying the 2nd Order NF transformation $y = z + h2(z)$, the normal dynamics reads:

$$ \dot{z} = \Lambda z + DF2(z)h2(z) + F3(z) + \mathcal{O}(4) \quad (17) $$

where:

1) $DF2(z)h2(z)$ are third order terms coming from the 2nd Order NF transformation used to cancel the 2nd order terms;

2) $F3(z)$ are the original 3rd order terms from the system (3).

Neglecting terms of order higher than 3, for the jth variable (17) results in:

$$ \dot{z}_j = \lambda_j z_j + \sum_{p=1}^{n} \sum_{q=1}^{n} \sum_{r=1}^{n} C_{pqr}^{j} z_p z_q z_r \quad (18) $$

where $C_{pqr}^j = \sum_{l=1}^n (F2_{pl}^j + F2_{lp}^j)h2_{qr}^p + F3_{pqr}^j$

To simplify as much as possible the third order terms, the following third-order NF transformation is applied to (17) [16]:

$$ z = w + h3(w) \quad (19) $$

where $h3$ is a polynomial of $w$ containing only 3rd-order terms. Applying transformation (19) to (17), it leads to:

$$ \dot{w} = \Lambda w + DF2(w)h2(w) + \Lambda h3(w) - Dh3(w)\Lambda w + F3(w) + \mathcal{O}(5) \quad (20) $$

It should be emphasized that there are neither terms of order 2 nor terms of order 4 in (20), and terms of order 3 are:

1) $DF2(w)h2(w)$ that comes from the 2nd-order NF transformation used to cancel the 2nd-order terms;

2) $\Lambda h3(w) - Dh3(w)\Lambda w$ that comes from the use of transformation (19) in order to cancel the 3rd-order terms;

3) $F3(w)$, the original 3rd order terms of system (16).

For the j-th variable, (20) can be written as:

$$ \dot{w}_j = \lambda_j w_j + \sum_{p=1}^{n} \sum_{q=1}^{n} \sum_{r=1}^{n} [C_{pqr}^{j} - (\lambda_p + \lambda_q + \lambda_r - \lambda_j)h_{pqr}^{3}] \times w_p w_q w_r + \mathcal{O}(5) \quad (21) $$

3) Elimination of Third-order Terms and Effects of the Resonant Terms: The next step is dedicated to the elimination of a maximum number of 3rd-order terms. This elimination is based on the same procedure as the one used to eliminate the second order terms. Equation (21) shows that third order terms can be eliminated with setting $h3$ computed by:

$$ h3_{pqr}^{j} = \frac{F3_{pqr}^{j} + \sum_{l=1}^{n} (F2_{pl}^{j} + F2_{lp}^{j}) h2_{qr}^{p}}{\lambda_{p} + \lambda_{q} + \lambda_{r} - \lambda_{j}} \quad (22) $$

Computation of $h3$ are impossible if the condition $\lambda_p + \lambda_q + \lambda_r - \lambda_j \approx 0$ holds. A careful attention is then mandatory since

resonant terms are systematically present when the system exhibits undamped (or weakly damped) oscillating modes (pairs of complex eigenvalues close the real axis).

If we consider a weak-damped oscillating mode, composed of two conjugated poles $\lambda_{2l}$ and $\lambda_{2l-1}$, the necessary condition for eliminating the associated third-order term $w_j w_{2l} w_{2l-1}$ will not be met for some j leading to the impossibility of computing the coefficient $h3_{j2l2l-1}^j$. As a consequence, term $w_j w_{2l} w_{2l-1}$ cannot then be eliminated and must be kept in the normal dynamics [25], [34].

If we consider that the system possesses $M$ weakly-damped oscillatory modes, $M$ third-order terms can thus not be eliminated from the normal dynamics of the jth variable ($j \in M$). Apart from internal resonances due to commensurability relationships between eigenvalues, as for example when $\omega_1 = 3\omega_2$, $\omega_1 = 2\omega_2 + \omega_3$ or $\omega_1 = \omega_2 + \omega_3 + \omega_4$, all the other third-order terms are eliminated by third-order transformation (19).

Neglecting terms with an order higher than three, the normal dynamics of a system composed of $M$ weakly-damped oscillatory modes are then:

$$ \dot{w}_j = \lambda_j w_j + \sum_{l=1}^{M} c_{2l}^{j} w_j w_{2l} w_{2l-1}, \quad j \in M \quad (23) $$

$$ \dot{w}_j = \lambda_j w_j, \quad j \notin M \quad (24) $$

where $w_{2l-1}$ is the complex conjugate of $w_{2l}$ and coefficients $c_{2l}^j$ are defined as:

$$ \begin{aligned} c_{2l}^{j=2k} = C_{(2k)(2l)(2l-1)}^{2k} &+ C_{(2k)(2l-1)(2l)}^{2k} + C_{(2l)(2k)(2l-1)}^{2k} \\ &+ C_{(2l)(2l-1)(2k)}^{2k} + C_{(2l-1)(2l)(2k)}^{2k} + C_{(2l-1)(2k)(2l)}^{2k} \end{aligned} \\ c_{2l}^{2l} = C_{(2l)(2l)(2l-1)}^{2l} + C_{(2l)(2l-1)(2l)}^{2l} + C_{(2l-1)(j)(j)}^{2l} \quad (25) $$

Equation (23) shows that considering third-order terms in the Normal Dynamic changes the way of making the stability analysis. As shown in [13], [33], the third order terms can have a stabilizing or a destabilizing effect and the inspection of the sign of the eigenvalue real parts is not sufficient to predict the stability of the system.

### D. Comparison with Other Works where Third Order Terms were Considered

In [13], [33], third-order terms are present in the Taylor's series but only the second-order transformation is used. Some 3rd order terms are kept on the normal dynamics and leads for the oscillatory modes to:

$$ \dot{z}_j = \lambda_j z_j + \sum_{l=1}^{M} c_{2l}^j z_j z_{2l} z_{2l-1}, \quad j \in M \quad (26) $$

For the non-oscillatory modes, the normal dynamics are the same as expressed by (14).

In [16], [17], a 3rd-order NF transformation is proposed where the terms $DF2(w)h2(w)$ of (20) are not taken into account, leading to coefficients given by:

$$ \overline{h3_{pqr}^{j}} = \frac{F3_{pqr}^{j}}{\lambda_{p} + \lambda_{q} + \lambda_{r} - \lambda_{j}} \quad (27) $$
Since some 3rd order terms are omitted, the accuracy of the proposed normal dynamics is worse than the one based on co-efficients given by (22).

Moreover, the considered normal dynamics are:

$$ \dot{w}_j = \lambda_j w_j \quad (28) $$

and since all resonant terms in (28) are neglected, it fails to give more information about the system stability than the linearized small-signal stability analysis.

#### IV. NONLINEAR ANALYSIS BASED ON THE 3RD ORDER NORMAL FORMS

NF methods can be used to quantify modal interactions and to give information on the system stability. To evaluate the gain obtained by the 3rd order method over the 2nd order method, normal dynamics (23) and (24) have to be reconstructed with the z coordinates as:

$$ z_j = w_j + \sum_{p=1}^{n} \sum_{q=1}^{n} \sum_{r=1}^{n} h_3^{j}_{pqr} w_p w_q w_r \quad (29) $$

It can be seen from the reconstruction that $z_j$ is composed of two sets. A first set with coefficients $h_3^{j}_{pqr}$ that can be used to quantify the 3rd order modal interactions [17] and a second set where $w_j$ is influenced by $c_{2l}^j$, making then possible to use it for stability analysis [13], [33].

Based on the Normal Forms, $h_3^{j}_{pqr}$ and $c_{2l}^j$ can be used to deduce nonlinear indexes for quantifying the 3rd order modal interactions and to give information on some nonlinear stability margins.

##### A. 3rd Order Modal Interaction Index

As observed from (29), 3rd order oscillations are modeled by $h_3^{j}_{pqr} w_p w_q w_r$. Then, a third order Modal Interaction index $MI_3^{j}_{pqr}$ can be defined as:

$$ MI_3^{j}_{pqr} = \frac{|h_3^{j}_{pqr} w_p^0 w_q^0 w_r^0|}{|w_j^0|} \quad (30) $$

which indicates the participation of the 3rd order modal interaction with the frequency $\omega_p + \omega_q + \omega_r$ in the mode j. $\omega^0$ are the initial values of the variables at the moment of interest. Introducing $MI_3^{j}_{pqr}$ leads to a clearer idea of how the system variables interact with each others and a more precise identification of the source of the oscillations. Those additional informations are crucial, especially when some fundamental, 2nd and 3rd order modal interactions exhibit the same oscillatory frequency. In that case, it is difficult to identify the order of the modal interactions using a FFT analysis. NF analysis offers then the same services as time-domain simulations with signal processing, in the sense that FFT can only identify the oscillatory frequencies while NF analysis can also identify the source of these oscillatory frequencies.

##### B. Stability Index

Since the stability of oscillatory modes in (2) is consistent with (23) and (24), information on the stability can be obtained without performing a time-domain simulation assessment.

For the non-oscillatory modes, $\dot{w}_j = \lambda_j w_j$, the stability is determined by inspecting $\lambda_j$. For the oscillatory modes where $j \in M$, the normal dynamics are:

$$ \dot{w}_j = \lambda_j w_j + \sum_{l=1}^{M} c_{2l}^{j} w_j w_{2l} w_{2l-1}, \quad (31) $$

$$ = \left( \lambda_j + \sum_{l=1}^{M} c_{2l}^{j} |w_{2l}|^2 \right) w_j \quad (32) $$

with $|w_{2l}|$ the magnitude of $w_{2l}$.

A stability interaction index can then be defined as:

$$ SII_{2l}^{j} = c_{2l}^{j} |w_{2l}^{0}|^2, \quad j \in M \quad (33) $$

and, remarking that

$$ \lambda_j + \sum_{l=1}^{M} c_{2l}^{j} |w_{2l}^{0}|^2 = \sigma_j + \text{real} \left( \sum_{l=1}^{M} SII_{2l}^{j} \right) + j(\omega_j + \text{imag}(SII_{2l}^{j})) \quad (34) $$

it is shown that the real part of $SII$ gives an indication on the stability of the system and the imaginary part of $SII$ is in relation with the frequency shift of the mode.

##### C. Nonlinear Modal Persistence Index

The nonlinear indexes $MI_3$ proposed in a previous section indicate the nonlinear interactions at the moment when the disturbance is cleared. To quantify the 3rd order modal interaction in the overall dynamics, a persistence index is defined to indicate how long a nonlinear interaction will influence the dominant modes. Similar to (23) and (24) in [12], $Tr_3$ is the ratio of the time constant associated to a combination of modes over a dominant mode.

$$ Tr_3 = \frac{\text{Time constant for combination of modes } (\lambda_p + \lambda_q + \lambda_r)}{\text{Time constant for dominant mode } (\lambda_j)} \quad (35) $$

A small $Tr_3$ indicates a significant presence of a mode combination. For example, if $Tr_3 = 1$, it means that a nonlinear interaction decays at the same speed as a dominant mode. If $Tr_3$ is very low, it means that the influence of a mode combination $(p, q, r)$ quickly vanishes compared to the dominant mode. A relatively high value of the product $SII \times Tr_3, MI_3 \times Tr_3$ tends to reveal the presence of a persistent modal interaction.

##### D. Stability Assessment

As indicated by (31), it is both the eigenvalues $\lambda_j$ and the stability indexes $SII_{2l}^j$ that give information concerning the system stability. Although $\lambda_j$ keeps constant, $\sum_{l=1}^{M} c_{2l}^j |w_{2l}|^2 w_j$ decays as time runs. For example, if $SII_{2l}^j$ is large, but $SII_{2l}^j \times Tr_3^j (2l)(2l-1)$ is small, the 3rd-order terms may stabilize the system at the beginning of the transient but a long-term instability can exist. Therefore, the time constant associated to $SII_{2l}^j$ must be taken into account when assessing the overall stability. An index of stability assessment can then be
<table><caption>TABLE I
LIST OF THE MNFS PRESENTED IN THIS WORK</caption><thead><tr><th>Method</th><th>Taylor's series</th><th>NF Trans.</th><th>Normal Form</th><th>n<sup>th</sup> order NF Terms</th><th>Related Equation</th></tr></thead><tbody><tr><td>2-2-1</td><td>2nd order</td><td>order 2</td><td>order 1</td><td>/</td><td>(14)</td></tr><tr><td>3-2-3S</td><td>3rd order</td><td>order 2</td><td>order 3</td><td>S</td><td>(26)</td></tr><tr><td>3-3-1</td><td>3rd order</td><td>order 3</td><td>order 1</td><td>/</td><td>(28)</td></tr><tr><td>3-3-3</td><td>3rd order</td><td>order 3</td><td>order 3</td><td>/</td><td>(23)(24)</td></tr></tbody></table>

defined as:

$$SI_j = \sigma_j + \sum_{l=1}^{M} \text{real}(SII_{2l}^j) \times \text{Tr}3_{j(2l-1)(2l)}^j, j \in M \quad (36)$$

with $\sigma$ being the real value of an eigenvalue ($\lambda_j = \sigma_j + j\omega_j$). It immediately follows from this definition that, if $SI_j > 0$ third order terms will have a destabilizing effect and, if $SI_j < 0$ third order terms will have a stabilizing effect. If $SI_j = 0$, the system may stay in limit cycles, or switch between the stable and the unstable phase. Different cases should be discussed and further investigations should be made.

Compared to [13], this approach has a larger range as it takes into account the persistence time of a mode interaction and it is appropriate for cases where $\sigma_j = 0$.

## V. CASE-STUDY

### A. Summary Analysis of the Different Methods Reviewed and Proposed in the Work

The methods presented in this work can be compared on four basic features. Each method is then labeled using three digits and one optional letter where:

1) The first digit gives the order of the Taylor's expansion of the system's dynamic;

2) The second digit gives the order of the Normal Form transformation used;

3) The third digit gives the order of the considered normal dynamics;

4) The optional letter indicates the fact that some terms have not been taken into account in the normal dynamics (S for Simplified)

Table I lists the four MNFs presented and compared in this work. In all the chosen test systems, the variables are in per units, while the rotor angles are measured in rad, considering the scaling problem [35].

### B. The IEEE 4-Machine Test System

The chosen test system to assess the advantages of the proposed method is a well known IEEE standard system, the Kundur's 2-area 4-machine system shown in Fig. 1. It is a classical system suitable for the analysis of modal oscillations for the validation of small-signal analysis [1] and for the validation of 2nd order Normal Form analysis [12].

The generators are modeled using a two-axis fourth-order model and a thyristor exciter with a Transient Gain Reduction. Loads L1 and L2 are modeled as constant impedances and no

<figure><img src="image_3.png" /><figcaption>Fig. 1. IEEE 4 machine test system: 2-area 4-generator.</figcaption></figure>

Power System Stabilizer (PSS) are used. To make the system robust, each area is equipped with large capacitor banks to avoid a voltage collapse. The data for the system and the selected case are provided in the Appendix section. The full numerical time domain simulation is used to assess the performance of the MNFs in approximating the nonlinear system dynamics, based on the well validated demo power_PSS in Matlab 2015a.

The oscillatory modes for the selected cases are listed in Tables II and III with the associated pseudo frequency, time-constant and dominate states. It gives a clear picture of the physical property of the system dynamics by small-signal analysis. Some stable real modes are not listed for the sake of compactness.

The cases analyzed in this section were selected to highlight information provided by the 3rd order Normal Form analysis. The selected system operating conditions for the study is a highly stressed case where the system is close to the voltage collapse, characterized by a tie line flow of 420MW from Area1 to Area2. To consider the emergence of renewable energy based generators, some modifications are made compared to the conventional small-signal and 2nd-order-based NF analysis that have been already conducted on the same test case. The powers generated by generators G1 and G2 in Area1 are unbalanced to consider the production of energy from distant renewable energy based generators in large areas.

Two cases are considered:

Case 1: This case represents a poorly damped situation, where the damping ratio of the inter-area mode is only 5.9% and the system is at its limit of stability according to a linear analysis. A three phase short-circuited fault is applied at Bus 7 and after 0.41 s, line B and line C are tripped. The exciter gain $K_a$ is set as 150 for all generators.

Case 2: This case represents a situation where the damping ratio of two oscillatory modes is negative (modes (5, 6) and (7, 8), as indicated in Table III). The negative damping is introduced by changing the gain of the 4 thyristor exciters by a higher value ($K_a = 240$). A three phase short-circuited fault is applied at Bus 7 and after 0.10s, line B and line C are tripped.

The transient analysis of the overall system based on the NF analysis is presented in the next section, where G1 and G2 exhibit predominant electromechanical oscillations.
<table><caption>TABLE II
OSCILLATORY MODES: CASE 1</caption><thead><tr><th>Mode #</th><th>Eigenvalue</th><th>Pseudo-Freq (Hz)</th><th>Damping Ratio (%)</th><th>Time Constant t = 1/ζω<sub>n</sub></th><th>Dominant States</th></tr></thead><tbody><tr><td>5, 6</td><td>-3.46 ± j93.2</td><td>14.8</td><td>3.17</td><td>0.289</td><td>E'<sub>q1</sub>, E'<sub>q2</sub>, E'<sub>q3</sub>, E'<sub>q4</sub></td></tr><tr><td>7, 8</td><td>-5.94 ± j88.3</td><td>14.0</td><td>6.71</td><td>0.168</td><td>Control unit G2</td></tr><tr><td>9, 10</td><td>-18.1 ± j62.8</td><td>9.99</td><td>27.7</td><td>0.055</td><td>Control unit (G3, G4)</td></tr><tr><td>11, 12</td><td>-17.3 ± j64.5</td><td>10.3</td><td>25.9</td><td>0.057</td><td>Control unit G1</td></tr><tr><td>13, 14</td><td>-1.51 ± j6.73</td><td>1.07</td><td>21.9</td><td>0.662</td><td>Local, Area 2(δ<sub>1</sub>, δ<sub>2</sub>, ω<sub>1</sub>, ω<sub>2</sub>)</td></tr><tr><td>15, 16</td><td>-1.69 ± j6.51</td><td>1.03</td><td>25.1</td><td>0.593</td><td>Local, Area 1(δ<sub>1</sub>, δ<sub>2</sub>, ω<sub>1</sub>, ω<sub>2</sub>)</td></tr><tr><td>17, 18</td><td>-0.135 ± j2.28</td><td>0.36</td><td>5.9</td><td>7.41</td><td>Inter-area (δ<sub>1</sub>, δ<sub>2</sub>, δ<sub>3</sub>, δ<sub>4</sub>)</td></tr></tbody></table>

<table><caption>TABLE III
OSCILLATORY MODES: CASE 2</caption><thead><tr><th>Mode #</th><th>Eigenvalue</th><th>Pseudo-Freq (Hz)</th><th>Damping Ratio (%)</th><th>Time Constant t = 1/ζω<sub>n</sub></th><th>Dominant States</th></tr></thead><tbody><tr><td>5, 6</td><td>3.61 ± j107.9</td><td>17.2</td><td>-3.34</td><td>-0.277</td><td>E'<sub>q1</sub>, E'<sub>q2</sub>, E'<sub>q3</sub>, E'<sub>q4</sub></td></tr><tr><td>7, 8</td><td>1.73 ± j104.1</td><td>16.6</td><td>-1.66</td><td>-5.79</td><td>Control unit G2</td></tr><tr><td>9, 10</td><td>-11.4 ± j77.8</td><td>12.4</td><td>14.5</td><td>0.088</td><td>Control unit (G3, G4)</td></tr><tr><td>11, 12</td><td>-10.3 ± j79.9</td><td>12.7</td><td>12.8</td><td>0.097</td><td>Control unit G1</td></tr><tr><td>13, 14</td><td>-1.56 ± j6.64</td><td>1.05</td><td>22.9</td><td>0.641</td><td>Local, Area 2(δ<sub>1</sub>, δ<sub>2</sub>, ω<sub>1</sub>, ω<sub>2</sub>)</td></tr><tr><td>15, 16</td><td>-1.72 ± j6.44</td><td>1.02</td><td>25.8</td><td>0.581</td><td>Local, Area 1(δ<sub>1</sub>, δ<sub>2</sub>, ω<sub>1</sub>, ω<sub>2</sub>)</td></tr><tr><td>17, 18</td><td>-0.132 ± j2.23</td><td>0.36</td><td>5.92</td><td>7.55</td><td>Inter-area (δ<sub>1</sub>, δ<sub>2</sub>, δ<sub>3</sub>, δ<sub>4</sub>)</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Comparison of different NF approximations: Case 1.</figcaption></figure>

### C. Case 1-2: Benefits of Using 3rd Order Approximation for Modal Interactions and Stability Analysis

When the system is close to its limits of stability as depicted by Case 1, Fig. 2 shows that the linear analysis gives wrong predictions concerning the dynamic behavior of the system. 2nd-order and 3rd-order Normal Form approximations can both better model the system dynamic response than the Linear Method.

The proposed third-order approximation, called 3-3-3, is superior over the other methods since it makes possible to model the interactions of oscillatory of non-oscillatory modes up to

<figure><img src="image_5.png" /><figcaption>Fig. 3. FFT Analysis of different NF approximations: Case 1.</figcaption></figure>

order 3, contributing to a better modeling of the frequency variation of the oscillations as recently formulated in [36] and [37].

Concerning Case 2, the system is such that the linear analysis leads to the computation of two oscillatory modes with a negative damping (See modes (5,6) and (7,8) given by Table III). According to the conventional small-signal stability analysis, the overall system is then considered as unstable. However, as shown in Fig. 4, the system's nonlinearities contribute to the overall stability of the system [13]. The comparison of the different NF approximations shows that only methods keeping 3rd-order terms in the normal dynamics are able to predict the stability of the system (3-2-3S and proposed 3-3-3 approximations).
<figure><img src="image_6.png" /><figcaption>Fig. 4. Comparison of different NF approximations: Case 2.</figcaption></figure>

<table><caption>TABLE IV
INITIAL CONDITIONS</caption><thead><tr><th>j</th><th>y0</th><th>z0221</th><th>z0331</th><th>z0333</th></tr></thead><tbody><tr><td>5</td><td>5.14 + j2.93</td><td>3.89 + j2.38</td><td>3.88 + j2.38</td><td>0.3105 + j0.3263</td></tr><tr><td>7</td><td>-6.36 - j2.73</td><td>-3.24 - j0.95i</td><td>-3.24 - j0.95</td><td>-0.36 - j2.49</td></tr><tr><td>9</td><td>-1.35 - j1.39</td><td>0.48 - j3.14</td><td>0.58 - j2.3</td><td>0.75 - j1.59</td></tr><tr><td>11</td><td>-9.84 - j4.25</td><td>-6.07 - j4.58</td><td>-6.27 - j4.58</td><td>-0.54 - j3.84</td></tr><tr><td>13</td><td>-0.25 - j0.21</td><td>-0.27 - j0.08</td><td>-0.226 - j0.085</td><td>-0.20 - 0.22i</td></tr><tr><td>15</td><td>0.32 + j0.10</td><td>0.04 - 0.0012i</td><td>0.05 - j0.0012</td><td>0.072 - j0.05</td></tr><tr><td>17</td><td>0.87 + j2.52</td><td>0.82 + j2.16</td><td>0.84 + j2.10</td><td>1.70 + j1.31i</td></tr></tbody></table>

### D. Nonlinear Analysis Based on Normal Forms

It has been validated by the time-domain simulation in the previous section that the normal dynamics well approximate the system dynamics. In this section, the MNF is used to make quantitative analyses of the system dynamics with using the nonlinear indexes.

1) **Initial Condition and Magnitude at Fundamental Frequency:** The initial conditions in the Jordan form ($y^0$) and for methods 2-2-1 ($z0221$), 3-3-1 ($z0331$) and 3-3-3 ($z0333$) are listed in Table IV. It has to be noticed that method 3-2-3S has the same initial conditions as method 2-2-1, i.e. $z0323 = z0221$.

From the initial condition, it can be seen that the magnitude of the ratio of $y_{17}^0$, $z0221_{17}$, $z0331_{17}$, $z0333_{17}$ over $z0333_{17}$ is respectively equals to 124%, 107.8%, 107.7% and 100%, which approximately matches the ratios at the fundamental frequency of the different curves in Fig. 3.

2) **Distribution of the Frequency Spectrum:** Comparing results of Table V with data extracted on Fig. 3, it is seen that: 1) $MI3_{17,17,17}^{17} = 0.11 = 11\%$ approximately matches the magnitude ratio of the 3rd order component at approximately 1 Hz; 2) ignoring $DF2h2$, $h3_{pqr}^j$ is small and leads to a too modest prediction of the 3rd order interaction $MI3_{pqr}^j = 0.056 = 5.6\%$.

A deeper analysis can be made using NF methods compared to FFT analysis. For example, $MI3_{17,26,27}^{17} = 0.45$ indicates

<table><caption>TABLE V
3RD ORDER COEFFICIENTS FOR MODE (17, 18)</caption><thead><tr><th>p, q, r</th><th>(DF2h2)<sub>pqr</sub><sup>j</sup></th><th>F3<sub>pqr</sub><sup>j</sup></th><th>h3<sub>pqr</sub><sup>j</sup></th></tr></thead><tbody><tr><td>(17, 17, 17)</td><td>0.046 - j0.118</td><td>0.029 + j0.034</td><td>0.017 - j0.019</td></tr><tr><td>MI3<sub>pqr</sub><sup>j</sup></td><td><u>h3<sub>17,17,17</sub><sup>17</sup></u></td><td><u>MI3<sub>pqr</sub><sup>j</sup></u></td><td>Tr3</td></tr><tr><td>0.1220</td><td>0.0068 + j0.008</td><td>0.0560</td><td>0.33</td></tr></tbody></table>

<table><caption>TABLE VI
RESONANT TERMS ASSOCIATED WITH MODE (17, 18)</caption><thead><tr><th>j</th><th>SII<sub>18</sub><sup>j</sup></th><th>&lambda;<sub>j</sub> + &lambda;<sub>2l</sub> + &lambda;<sub>2l-1</sub></th><th>Tr3</th></tr></thead><tbody><tr><td>5</td><td>-1.30 - j2.59</td><td>-3.73 + j93.2</td><td>0.93</td></tr><tr><td>7</td><td>0.48 + j0.96</td><td>-6.21 + j88.3</td><td>0.96</td></tr><tr><td>9</td><td>0.0117 + j0.17</td><td>-18.37 + j62.8</td><td>0.98</td></tr><tr><td>11</td><td>-0.51 - j1.52</td><td>-17.57 + j64.5</td><td>0.98</td></tr><tr><td>13</td><td>0.074 - j0.25</td><td>-1.78 + j6.73</td><td>0.85</td></tr><tr><td>15</td><td>0.89 + j0.532</td><td>-1.96 + j6.51</td><td>0.86</td></tr><tr><td>17</td><td>0.12 - j0.678</td><td>-0.405 + j2.28</td><td>0.33</td></tr></tbody></table>

<table><caption>TABLE VII
STABILITY INDEXES OF MODE (5, 6), (7, 8)</caption><thead><tr><th>2l - 1</th><th>SII<sub>2l</sub><sup>5</sup></th><th>Tr3</th><th>real(SII) × Tr3<sup>5</sup></th></tr></thead><tbody><tr><td>5</td><td>-53.54 - j109.77</td><td>0.33</td><td>-17.6682</td></tr><tr><td>9</td><td>-3.37 - j6.34</td><td>-0.18</td><td>-0.63</td></tr><tr><td>13</td><td>-3.06 - j7.56</td><td>7.50</td><td>-22.95</td></tr><tr><td>15</td><td>-6.5 - j11.7</td><td>22.50</td><td>-146.25</td></tr><tr><td>17</td><td>-1.60 - j4.15</td><td>1.08</td><td>-1.72</td></tr><tr><td colspan="4">SI<sup>5</sup> = 3.61-17.6682 - 0.63 + 22.95 - 146.25 - 1.72 = -185.6082</td></tr></tbody></table><table><thead><tr><th>2l - 1</th><th>SII<sub>2l</sub><sup>7</sup></th><th>Tr3</th><th>real(SII) × Tr3<sup>7</sup></th></tr></thead><tbody><tr><td>5</td><td>-105.3 - j221.7</td><td>0.2</td><td>-21.06</td></tr><tr><td>9</td><td>-5.0 - j8.94</td><td>0.0821</td><td>-0.41</td></tr><tr><td>13</td><td>-3.21 - j8.81</td><td>1.24</td><td>-3.99</td></tr><tr><td>15</td><td>-4.57 - j8.67</td><td>1.01</td><td>-4.62</td></tr><tr><td>17</td><td>0.29 + j1.31</td><td>1.1801</td><td>0.34</td></tr></tbody></table><p>SI<sup>7</sup> = 1.73 - 21.06 - 0.41 - 3.99 - 4.62 + 0.34 = -28.01</p>

a strong nonlinear interaction, however, since $Tr = 6.7341 \times 10^{-4}$ and $MI3 \times Tr = 3.03 \times 10^{-4}$ such a short duration will be difficult to be captured by FFT analysis.

3) **Frequency Shift of the Fundamental Component:** It can be also noted that there is a shift in the fundamental frequency. As already mentioned, it is indicated by the imaginary part of index $SII_{2l}^j$. Among all those coefficients, coefficients $SII_{18}^j$ are predominant and $Tr3$ indicates a long-time influence as listed in Table VI. Although the real part of $SII_{18}^j$ is too small to contribute to the stability compared to $\lambda_j$, its imaginary part indicates that there is a frequency shift added to the eigen-frequency.

This analysis corroborates with the analysis of the frequency spectrum in Fig. 3. If there are no resonant terms in the normal forms, the fundamental frequency of the system oscillatory dynamics will be exactly as $\omega_j$.

4) **Stability Assessment:** Seen from Table VII, the nonlinear interactions can enhance (e.g $MI3_{5,5,6}^j$) or weaken the stability (e.g $MI3_{7,17,18}^j$). Only taking into account of one specific $MI3$ without considering its time pertinence will lead to a wrong prediction. In this sense, the stability assessment proposed by
<figure><img src="image_7.png" /><figcaption>Fig. 5. 16-Machine 5-area test system [38].</figcaption></figure>

method 3-3-3 is more rigorous and comprehensive than method 3-2-3. As $SI^5 \ll 0, SI^7 \ll 0$, modes (5, 6) and (7, 8) are essentially stable, as validated by the time domain analysis shown in Fig. 4.

The proposed nonlinear stability indicators make possible to better use the power transfer capability of the power grid. For example, in the studied case, the high-gain exciter can improve the system response to a fault, which is not shown using the conventional eigen-analysis.

### E. Applicability to Larger Networks With More Complex Power System Models

For larger networks with more complex models, method 3-3-3 is more demanding since the modeling of modal interactions is more complex, such as for the case of the New England New York 16 machine 5 area system [38], whose typology is shown in Fig. 5. It is composed of five geographical regions out of which NETS and NYPS are represented by a group of generators whereas, the power imported from each of the three other neighboring areas are approximated by equivalent generator models (G14 to G16). G13 also represents a small sub-area within NYPS. Generators G1 to G8 and G10 to G12 have DC excitation systems (DC4B); G9 has a fast static excitation (ST1A), while the rest of the generators (G13 to G16) have manual excitation as they are area equivalents instead of being physical generators [38]. The realistic parameters and well validated Simulink models can be found in [38], [39], where the generators are modeled with the sub-transient models with four equivalent rotor coils. There are 15 pairs of electromechanical modes, among which there are 4 inter-area modes. This system is unstable when no PSS or only one PSS is installed [38].

When placing PSSs on G1 to G12 (the maximum number of possible PSSs) all the local modes are damped, and 3 inter-modes are poorly damped. No more information is available from the linear analysis to damp the inter-area modes [38]. When a three-phase fault is applied to G13 and cleared in 0.25s, the machines exhibit nonlinear inter-area oscillations and finally

<figure><img src="image_8.png" /><figcaption>Fig. 6. FFT analysis of generator rotor angles when G1 to G12 are equipped with PSSs.</figcaption></figure>

<table><caption>TABLE VIII
INTER-AREA MODES OF THE 16 MACHINE 5 AREA SYSTEM WITH PSSS PRESENT ON G6 AND G9</caption><thead><tr><th>Mode</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>Damping ratio</td><td>3.35%</td><td>0.55%</td><td>1.58%</td><td>2.76%</td></tr><tr><td>Frequency (Hz)</td><td>0.7788</td><td>0.6073</td><td>0.5226</td><td>0.3929</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 7. FFT analysis of generator rotor angles when only on G6 and G9 are equipped with PSSs.</figcaption></figure>

damp out to a steady-state equilibrium in 50s, as shown by the frequency spectrum in Fig. 6. As all the local-modes are well damped, the components at frequency higher than 1 Hz are expected as nonlinear interactions.

When placing PSSs only on G6 and G9 (the minimum number of PSSs to ensure the stability of the system) all inter-area modes are poorly damped as listed in Table VIII. When there is a three-phase short circuit fault applied near the end of G13, generators in all the 5 areas exhibit poorly damped electromechanical oscillations and finally damp out to a steady-state equilibrium in a large time. Since the linear participation factor of G13 is [0.0075 0.7997 0.0214 0.6871] on the inter-area modes 1-4, only Mode 2 and Mode 4 will be effectively excited while Mode 1 and Mode 3 are trivially excited. However, as observed from the frequency spectrum in Fig. 7, there are significant components at 0.5 Hz (G14, the fundamental frequency of Mode 3),
<table><caption>TABLE IX
SEARCH FOR INITIAL CONDITIONS IN THE NF COORDINATES</caption><thead><tr><th rowspan="2">Method</th><th colspan="2">Case 1</th><th colspan="2">Case 2</th></tr><tr><th>Iterations</th><th>Resolution</th><th>Iterations</th><th>Resolutions</th></tr></thead><tbody><tr><td>3-3-3</td><td>12</td><td>1.294e-12</td><td>17</td><td>1.954e-6</td></tr><tr><td>2-2-1/3-2-3</td><td>3</td><td>1.799e-11</td><td>6</td><td>1.838e-06</td></tr><tr><td>3-3-1</td><td>7</td><td>1.187e-09</td><td>11</td><td>7.786e-09</td></tr></tbody></table>

<table><caption>TABLE X
PERFORMANCE EVALUATION OF THE STUDIED NF APPROXIMATIONS</caption><thead><tr><th>Method</th><th>Resolution</th><th>Order of Modal Interaction</th><th>Transient Stability</th><th>Applicability to 3rd Order Modes</th></tr></thead><tbody><tr><td>Linear</td><td>O(2)</td><td>1<sup>st</sup></td><td>Linear</td><td>None</td></tr><tr><td>2-2-1</td><td>O(3)</td><td>2<sup>nd</sup></td><td>Linear</td><td>None</td></tr><tr><td>3-2-3S</td><td>O(3)</td><td>2<sup>nd</sup></td><td>Non-linear</td><td>Oscillatory</td></tr><tr><td>3-3-1</td><td>O(3)</td><td>3<sup>rd</sup></td><td>Linear</td><td>Non-oscillatory</td></tr><tr><td>3-3-3</td><td>O(5)</td><td>3<sup>rd</sup></td><td>Non-linear</td><td>All</td></tr></tbody></table>

at 0.8 Hz (G15, the fundamental frequency of Mode 1), around 1.0 Hz (G6, G9, G13, G14, G15, the frequency is not corroborated with any inter-area mode). Therefore, *SII* and *MI3* are expected to give additional information to identify the modal interactions.

In this case, method 3-3-3 may provide more information to: 1) reduce the number of PSSs; 2) damp the inter-area modes. The siting of PSSs based on NF analysis is not the issue to be dealt with in this paper, and it can be found in [28].

## VI. FACTORS INFLUENCING NORMAL FORM ANALYSIS

### A. Computational Burden

The essence of NF analysis is to calculate the nonlinear indexes to predict the modal interactions and to give information on the parameters influencing the system stability, which is composed of two phases of computations:

1) the SEP Initialization in order to obtain the eigenmatrix $\lambda$ and matrices $F2$ and $F3$, which depends on the post-fault SEP;

2) the Disturbance Initialization in order to obtain the initial points in the NF coordinates, which depend on the disturbances the system experiences.

The values of nonlinear indexes depend both on the SEP and the disturbances. Item 2 has been discussed in detail in [12], [27] while Item 1 is somewhat neglected in the literature.

1) *Search for the Initial Conditions:* In this paper, the search for the initial conditions is performed using the Newton-Raphson (NR) method [40], the starting search point is $y^0$, and it converges in several iterations, as shown in Table IX. This is because the normal coefficients $h2$ and $h3$ are small in size. When these coefficient are large in size (near strong resonance case), the search for $z^0, w^0$ can be extremely slow and even fails to converge. A more robust algorithm is proposed in [27] to circumvent some disadvantages of NR method, which can be adopted also for the 3rd order NF methods.

2) *SEP Initialization:* If the power load changes, the SEP initialization must be restarted. The calculation of the nonlinear matrices $F2$ and $F3$ can be extremely tedious. The problem is that matrices $F2$ and $F3$ are composed of complex values. Using the same Matlab function to calculate matrices $A, H2, H3, \Lambda, F2$ and $F3$ for the IEEE 4 machine case it leads to differences from 200s to 2000s. Reducing the time needed to multiply complex matrix can be a direction to optimize the program.

Once the perturbation model around the SEP is established in Jordan form, calculation of nonlinear indexes can be computed in a short time.

### B. Strong Resonance and Model Dependence

Validated using 4th order or higher order generator models, the proposed method 3-3-3 also inherits the applicability to systems modeled by 2nd order generator models as method 3-3-1 used in [30], and 3rd order generator as method 3-2-3 used in [13].

The technique works where methods 3-3-1 and 3-3-1 are not applicable. For example, using classical models with zero mechanical dampings, the eigenvalues will be pure imaginary, i.e $\sigma_j=0$. In this case, method 3-3-1 fails to be applied, since $\overline{h_j^3}_{2l-12l} = \infty$, (Mode (2l-1, 2l) being a conjugate pair that leads to a strong resonance).

In addition, as the eigenvalues are purely imaginary, the stability boundary proposed in [13], [33] [see (13) and (22)] will be wrongly predicted, as the stability boundary will be calculated

$$R_j = \sqrt{-\frac{\text{real}(\lambda_j)}{\text{real}(c_{2k}^{j=2k})}} = 0.$$

Method 3-3-3 is more complete, as it can make possible to predict both the importance of modal interactions and to proposes nonlinear stability indexes for a broader ranges of implemented models.

## VII. CONCLUSION

### A. Significance of the Proposed Research

With the nonlinear indexes, the proposed method 3-3-3 makes possible to quantify the third-order modal interactions and offers some pertinent information for stability analysis, providing a better tool compared to the linear or other existing normal forms methods. The indication given by the nonlinear indexes are validated by time-domain simulation and FFT analysis. Besides, this paper gives a good review of existing MNFs along with a performance evaluation of the different NF approximations studied in this work (see Table X).

### B. Possible Applications of the Proposed 3-3-3 Method

Several potential applications of the method can be envisaged. Among them, let cite methods to locate power system stabilizers, sensibility analysis based on the higher-order participation factors to aid in the design of the power system structure and controller parameters tuning and nonlinear stability analysis to better predict the power transfer limits of the system.
# APPENDIX

<table><caption>TABLE XI
GENERATOR DATA IN PU ON MACHINE BASE</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>R<sub>a</sub></td><td>0.0025</td><td>x<sub>l</sub></td><td>0.002</td></tr><tr><td>x<sub>d</sub></td><td>1.80</td><td>τ'<sub>d0</sub></td><td>8.0</td></tr><tr><td>x<sub>q</sub></td><td>1.70</td><td>τ'<sub>q0</sub></td><td>0.40</td></tr><tr><td>x'<sub>d</sub></td><td>0.3</td><td>τ''<sub>d0</sub></td><td>0.0</td></tr><tr><td>x'<sub>q</sub></td><td>0.3</td><td>τ''<sub>d0</sub></td><td>0.0</td></tr><tr><td>x''<sub>d</sub></td><td>0.0</td><td>MVA<sub>base</sub></td><td>900</td></tr><tr><td>x''<sub>q</sub></td><td>0.0</td><td>H</td><td>6.5</td></tr></tbody></table>

<table><caption>TABLE XII
POWER FLOW DATA FOR CASE 1, CASE 2 IN STEADY STATE AFTER THE FAULT IS CLEARED: LOAD (CALCULATED BY MATLAB SIMPOWER)</caption><thead><tr><th>Load Bus</th><th>Voltage (pu/230kv)</th><th colspan="3">Z<sub>shunt</sub></th></tr></thead><tbody><tr><td>B5</td><td>1.029∠19.75°</td><td>P = 665.84MW, Q = -350.82Mvar</td><td></td><td></td></tr><tr><td>B6</td><td>0.987∠-47.37°</td><td>P = 1443.38MW, Q = -350.02Mvar</td><td></td><td></td></tr><tr><th>Generator</th><th>Bus Type</th><th>Terminal Voltage (pu/20kv)</th><th>P (MW)</th><th>Q (Mvar)</th></tr><tr><td>G1</td><td><i>PV</i></td><td>1.00∠14.02°</td><td>747.50</td><td>49.94</td></tr><tr><td>G2</td><td><i>swingbus</i></td><td>1.00∠0.0°</td><td>406.50</td><td>140.74</td></tr><tr><td>G3</td><td><i>PV</i></td><td>1.00∠-58.17°</td><td>539.25</td><td>113.99</td></tr><tr><td>G4</td><td><i>PV</i></td><td>1.00∠-66.11°</td><td>525.00</td><td>143.14</td></tr></tbody></table>

# REFERENCES

[1] P. Kundur, N. J. Balu, and M. G. Lauby, *Power System Stability and Control*, vol. 7. New York, NY, USA: McGraw-Hill, 1994.

[2] C.-M. Lin, V. Vittal, W. Kliemann, and A. Fouad, "Investigation of modal interaction and its effects on control performance in stressed power systems using normal forms of vector fields," *IEEE Trans. Power Syst.*, vol. 11, no. 2, pp. 781-787, May 1996.

[3] Y. X. Ni, V. Vittal, W. Kliemann, and A. A. Fouad, "Nonlinear modal interaction in hvdc/ac power systems with dc power modulation," *IEEE Trans. Power Syst.*, vol. 11, no. 4, pp. 2011-2017, Nov. 1996.

[4] J. Thapar, V. Vittal, W. Kliemann, and A. A. Fouad, "Application of the normal form of vector fields to predict interarea separation in power systems," *IEEE Trans. Power Syst.*, vol. 12, no. 2, pp. 844-850, May 1997.

[5] S. Saha, A. A. Fouad, W. H. Kliemann, and V. Vittal, "Stability boundary approximation of a power system using the real normal form of vector fields," *IEEE Trans. Power Syst.*, vol. 12, no. 2, pp. 797-802, May 1997.

[6] S. K. Starrett and A. A. Fouad, "Nonlinear measures of mode-machine participation [transmission system stability]," *IEEE Trans. Power Syst.*, vol. 13, no. 2, pp. 389-394, May 1998.

[7] G. Jang, V. Vittal, and W. Kliemann, "Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. i. general theory and procedure," *IEEE Trans. Power Syst.*, vol. 13, no. 2, pp. 401-407, May 1998.

[8] G. Jang, V. Vittal, and W. Kliemann, "Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. ii. case studies," *IEEE Trans. Power Syst.*, vol. 13, no. 2, pp. 408-413, May 1998.

[9] V. Vittal, W. Kliemann, Y. X. Ni, D. G. Chapman, A. D. Silk, and D. J. Sobajic, "Determination of generator groupings for an islanding scheme in the manitoba hydro system using the method of normal forms," *IEEE Trans. Power Syst.*, vol. 13, no. 4, pp. 1345-1351, Nov. 1998.

[10] S. Zhu, V. Vittal, and W. Kliemann, "Analyzing dynamic performance of power systems over parameter space using normal forms of vector fields-part I: Identification of vulnerable regions," *IEEE Trans. Power Syst.*, vol. 16, no. 3, pp. 444-450, Aug. 2001.

[11] S. Zhu, V. Vittal, and W. Kliemann, "Analyzing dynamic performance of power systems over parameter space using normal forms of vector

fields part II: Comparison of the system structure," *IEEE Power Eng. Rev.*, vol. 21, no. 6, pp. 71-71, Jun. 2001.

[12] J. Sanchez-Gasca *et al.*, "Inclusion of higher order terms for small-signal (modal) analysis: Committee report-task force on assessing the need to include higher order terms for small-signal (modal) analysis," *IEEE Trans. Power Syst.*, vol. 20, no. 4, pp. 1886-1904, Nov. 2005.

[13] H. Amano, T. Kumano, and T. Inoue, "Nonlinear stability indexes of power swing oscillation using normal form analysis," *IEEE Trans. Power Syst.*, vol. 21, no. 2, pp. 825-834, May 2006.

[14] T. Tian, X. Kestelyn, G. Denis, X. Guillaud, and O. Thomas, "Analytical investigation of nonlinear interactions between voltage source converters interconnected to a transmission grid," in *Proc. IEEE Power Energy Soc. General Meeting*, Jul. 2016, pp. 1-5.

[15] T. Tian, X. Kestelyn, and O. Thomas, "Third-order based analytical investigation of nonlinear interactions between voltage source converters interconnected to a transmission grid," in *Proc. 18th Eur. Conf. Power Electron. Appl.*, Sep. 2016, pp. 1-10.

[16] I. Martínez, A. Messina, and E. Barocio, "Perturbation analysis of power systems: Effects of second- and third-order nonlinear terms on system dynamic behavior," *Elect. Power Syst. Res.*, vol. 71, no. 2, pp. 159-167, 2004. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S037877960400029X

[17] Q. Huang, Z. Wang, and C. Zhang, "Evaluation of the effect of modal interaction higher than 2nd order in small-signal analysis," in *Proc. Power Energy Soc. General Meeting*, Jul. 2009, pp. 1-6.

[18] G. Iooss and M. Adelmeyer, *Topics in Bifurcation Theory and Applications*. Singapore: World Scientific, 1992.

[19] J. Guckenheimer and P. Holmes, *Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields*. New York, NY, USA: Springer-Verlag, 1983.

[20] L. Jzquel and C. H. Lamarque, "Analysis of non-linear dynamical systems by the normal form theory," *J. Sound Vib.*, vol. 149, no. 3, pp. 429-459, 1991.

[21] C. Touzé, O. Thomas, and A. Chaigne, "Hardening/softening behaviour in non-linear oscillations of structural systems using non-linear normal modes," *J. Sound Vib.*, vol. 273, nos. 1/2, pp. 77-101, 2004.

[22] C. Touzé and M. Amabili, "Non-linear normal modes for damped geometrically non-linear systems: Application to reduced-order modeling of harmonically forced structures," *J. Sound Vib.*, vol. 298, nos. 4/5, pp. 958-981, 2006.

[23] H. Poincaré, *Les Méthodes Nouvelles de la Mécanique Céleste*. Paris, France: Gautheur-Villars, 1899.

[24] A. H. Nayfeh, *A Method of Normal Forms*. New York, NY, USA: Wiley, 2011.

[25] C.-H. Lamarque, C. Touzé, and O. Thomas, "An upper bound for validity limits of asymptotic analytical approaches based on normal form theory," *Nonlinear Dyn.*, vol. 70, no. 3, pp. 1931-1949, 2012. [Online]. Available: http://dx.doi.org/10.1007/s11071-012-0584-y

[26] N. Kakimoto, A. Nakanishi, and K. Tomiyama, "Instability of interarea oscillation mode by autoparametric resonance," *IEEE Trans. Power Syst.*, vol. 19, no. 4, pp. 1961-1970, Nov. 2004.

[27] E. Barocio, A. Messina, and J. Arroyo, "Analysis of factors affecting power system normal form results," *Elect. Power Syst. Res.*, vol. 70, no. 3, pp. 223-236, 2004. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0378779603003092

[28] S. Liu, A. R. Messina, and V. Vittal, "Assessing placement of controllers and nonlinear behavior using normal form analysis," *IEEE Trans. Power Syst.*, vol. 20, no. 3, pp. 1486-1495, Aug. 2005.

[29] R. Betancourt, E. Barocio, J. Arroyo, and A. Messina, "A real normal form approach to the study of resonant power systems," *IEEE Trans. Power Syst.*, vol. 21, no. 1, pp. 431-432, Feb. 2006.

[30] I. Martínez, A. R. Messina, and V. Vittal, "Normal form analysis of complex system models: A structure-preserving approach," *IEEE Trans. Power Syst.*, vol. 22, no. 4, pp. 1908-1915, Nov. 2007.

[31] R. Z. Davarani, R. Ghazi, and N. Pariz, "Nonlinear modal analysis of interaction between torsional modes and SVC controllers," *Elect. Power Syst. Res.*, vol. 91, pp. 61-70, 2012. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0378779612001277

[32] H. Lomei, M. Assili, D. Sutanto, and K. M. Muttaqi, "A new approach to reduce the nonlinear characteristics of a stressed power system by using the normal form technique in the control design of the excitation system," *IEEE Trans. Ind. Appl.*, vol. 53, no. 1, pp. 492-500, Jan. 2017.

[33] H. Amano and T. Inoue, "A new PSS parameter design using nonlinear stability analysis," in *Proc. Power Eng. Soc. General Meeting*, Jun. 2007, pp. 1-8.
[34] H. Dulac, "Solutions d'un système d'équations différentielles dans le voisinage de valeurs singulières," *Bull. de la Société Mathématique de France*, vol. 42, pp. 324-383, 1912. [Online]. Available: http://www.numdam.org/item?id=BSMF_1912_40_324_0

[35] I. Dobson and E. Barocio, "Scaling of normal form analysis coefficients under coordinate change," *IEEE Trans. Power Syst.*, vol. 19, no. 3, pp. 1438-1444, Aug. 2004.

[36] B. Wang and K. Sun, "Formulation and characterization of power system electromechanical oscillations," *IEEE Trans. Power Syst.*, vol. 31, no. 6, pp. 5082-5093, Nov. 2016.

[37] B. Wang, X. Su, and K. Sun, "Properties of the frequency-amplitude curve," *IEEE Trans. Power Syst.*, vol. 32, no. 1, pp. 826-827, Jan. 2017.

[38] A. K. Singh and B. C. Pal, "IEEE PES task force on benchmark systems for stability controls report on the 68-bus, 16-machine, 5-area system," *IEEE Power and Energy Soc.*

[39] [Online]. Available: http://www.sel.eesc.usp.br/ieee/

[40] M. Mikofski, "Newton-Raphson solver," 2014. [Online]. Available: https://fr.mathworks.com/matlabcentral/fileexchange/43097-newton-raphso n-solver