M Boukheddimi

F Bailly

P Soueres

B Watier

Human gait simulation from a reduced set of low-dimensional tasks using hierarchical control

Keywords: locomotion, hierarchical control, humanoid, robotics

Introduction

Human walking relies on a complex cyclic whole-body organization that involves many degrees of freedom, numerous muscles, and is intrinsically unstable. Many studies have described this motion at the kinematic, kinetic or dynamic levels. However, if researchers have accurately detailed the different sequences of the normal gait cycle, the criteria that underlie its coordination are still poorly understood. Some studies have simulated the gait by focusing on the motion of the lower limbs. Other works have proposed to simulate whole-body motions during gait using direct kinematic [START_REF] Maldonado | Full-body musculoskeletal model for muscle-driven simulation of human gait[END_REF], optimal control [START_REF] Felis | Synthesis of full-body 3-d human gait using optimal control methods[END_REF] or hierarchical control (De Lasa et al., 2016). However, few of them have tried to identify simple organizational principles to generate the walking motion and compare the obtained simulated results with the reference human data. In this paper, we propose to use the hierarchical tasks framework used in humanoid robotics for simulating the human gait from a reduced set of ordered tasks. We simulated human gait with a whole-body model using hierarchical control, and compared the simulated motion with experimental data obtained from a motion capture system. We hypothesized that the control of a reduced set of low-dimensional key tasks could lead to closely reproduce the human gait.

Methods

Subject and data

One healthy man volunteered for this study (age: 24 yr, height: 1.75 m, mass: 78 kg), the data was extracted from [START_REF] Carpentier | On the centre of mass motion in human walking[END_REF] database. Kinematic data was collected using 12 infrared cameras sampling at 200 Hz (Vicon, Oxford Metrics, Oxford, UK) recording 48 reflective markers placed on the subject, based on the International Society of Biomechanics (ISB) recommendations. Data was extracted and filtered with a 4th order butterworth filter after a residual analysis. Scaling parameters were computed according to Dumas anthropomorphic table [START_REF] Dumas | Adjustments to mcconville et al. and young et al. Body segment inertial parameters[END_REF]). The reference data that was used to benchmark the simulation results were extracted from the computation of the inverse kinematics with Opensim.

Model

A 3D whole-body skeletal model, involving 42 degrees of freedom (Dofs) (Table 1) was used to simulate and analyze the gait cycle. The model was frictionless as proposed in [START_REF] Maldonado | Full-body musculoskeletal model for muscle-driven simulation of human gait[END_REF]) (Figure 1). Opposite initial contact: left heel contact.

Gait simulation

As for the human body, the proposed anthropomorphic model is for the majority of tasks a highly redundant system. To manage this redundancy we used a Hierarchical Quadratic Programming solver (HQP) [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF] to simulate the gait. Given selected trajectory tasks as input, we applied the HQP secondorder inverse kinematics, which computes the accelerations of the segments to generate the expected motion. The gait cycle was split into four phases: two single support phases, and two double support phases. This decomposition is induced by the main support phases identified in the normalized walking cycle. The human gait, was simulated by using the following data as input: -The three-dimensional centre of mass trajectory of the whole-body.

-The projection of the heel path of the swinging foot in the transverse plane.

-Final joint positions, (i.e.: matrix shape 1*42) at the end of each phase. Considering that human subjects preferentially controlled the successive positions of their feet when walking, we have chosen to use as a reference the projection of the heel in the transverse plane, without considering its height. On the other hand, we have chosen to keep the three-dimensional centre of mass trajectory as a reference because it is highly related to the whole-body configuration during the walking. Thus, the hierarchy applied throughout our simulation was defined by the following tasks priority order : 1 st : Trajectory of the projection of the heel of the swinging foot in the transverse plane + fixed position of the projection of the heel of the stance foot in the transverse plane. 2 nd : Three dimensional trajectory of the whole-body centre of mass. 3 rd : Final joint positions. Therefore, the height of the swinging foot was not constrained and resulted from the control of the height of the body centre of mass. The simulated gait cycle obtained with this hierarchical motion generation approach was then compared to the subject experimental data measured during one cycle. For this comparison we used several metrics, which are: the Root Mean Square Error (RMSE) that provides the overall error between the reference and the simulation, the correlation coefficient of Pearson (r), which presents a measure of the intensity and the direction of the relationship between the reference and the simulation, and the range of motion (ROM).

Results and discussion

As a preliminary result, we were able to obtain a humanlike kinematic walking motion using a subsentially low number of constraints, regarding to previous works in the literature (see Figure 1). We analyzed the "hip, ankle, knee and shoulder" joint rotations. The simulation results were consistent with experimental reference data.

As an example, we present here the flexion-extension of the knee, which constitutes a major feature of bipedal walking. Figure 2 provides a comparison of the simulated and the reference flexion-extension of the right knees, during one walking cycle. During the flexion phase, the simulation correctly follows the reference. However, the maximal extension of the knee observed for the reference is not entirely reproduced by the simulation. The computed RMSE is 5 degrees. Also, the computed (r) between the reference and the simulation is 0.97, which shows a strong correlation between the two gaits. The ROM of the knee was: 60 degrees for the reference and 48 degrees for the simulated motion. The difference of 12 degrees is induced by a lack of maximal knee extension during the simulation. These results show that controlling only a few degrees of freedom, we can generate a gait close to the human one. It is a well-known idea in human motor control that complex motions should be generated by a reduced set of primitive elements [START_REF] Ijspeert | Dynamical Movement Primitives: Learning Attractor Models for Motor Behaviors[END_REF]. Several representations such as, central pattern generators, motion primitives or synergies have been developed in this direction. We believe that the concept of tasks hierarchy used in robotics to cope with redundancy, notably in humanoid robotics, is closely related to these concepts. The work presented in this paper provides an efficient application of this approach to human waking.

Conclusions

In this paper, we generated human gait using the hierarchical task framework usually developed for humanoid robots control. A hierarchy of tasks was applied to a whole-body model for gait motion generation.

It was possible to reproduce a trajectory close to the reference human walking by controlling a small number of low dimensional ordered tasks. The results of the comparison of the simulated rotation of the joint with the human experimental data confirm that the control proposed in our method is closely related to the one exerted by the central nervous system during human walking. In order to study the robustness and the reliability of our assumption in a future study, we will apply this task hierarchy to a larger panel of subjects, and study in more details the 3D joint rotations.

Figure 1

 1 Figure 1 Snapshots of the simulated walking: (a): Initial contact: right heel contact, (b): Left toe-off, (c):Opposite initial contact: left heel contact.2.3 Gait simulationAs for the human body, the proposed anthropomorphic model is for the majority of tasks a highly redundant system. To manage this redundancy we used a Hierarchical Quadratic Programming solver (HQP)[START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF] to simulate the gait. Given selected trajectory tasks as input, we applied the HQP secondorder inverse kinematics, which computes the accelerations of the segments to generate the expected motion. The gait cycle was split into four phases: two single support phases, and two double support phases. This decomposition is induced by the main support phases identified in the normalized walking cycle. The human gait, was simulated by using the following data as input: -The three-dimensional centre of mass trajectory of the whole-body.-The projection of the heel path of the swinging foot in the transverse plane.-Final joint positions, (i.e.: matrix shape 1*42) at the end of each phase. Considering that human subjects preferentially controlled the successive positions of their feet when walking, we have chosen to use as a reference the projection of the heel in the transverse plane, without considering its height. On the other hand, we have chosen to keep the three-dimensional centre of mass

Figure 2

 2 Figure 2 Joint rotations of the right knee: flexion (-), extension (+)

  Degrees of freedom in the whole-body model

		Lower limbs	Upper limbs	Pelvis back	Neck and head
	Dofs	7 + 7	8 + 8	6	3	3
	Table 1


Student Member, IEEE Tian Tian
email: tian.tian@ensam.eu

Member, IEEE Xavier Kestelyn
email: xavier.kestelyn@ensam.eu

Olivier Thomas
email: olivier.thomas@ensam.eu

Hiroyuki Amano
email: h-amano@criepi.denken.or.jp

Fellow, IEEE Arturo Roman Messina

Xavier Kestelyn Tian

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

Keywords: Interconnected power system, methods of normal forms, nonlinear modal interaction, power system dynamic, stability. I

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.

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-basedgenerators, 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 exiting 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 smallsignal 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 [START_REF] Kundur | Power System Stability and Control[END_REF].

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 [START_REF] Lin | Investigation of modal interaction and its effects on control performance in stressed power systems using normal forms of vector fields[END_REF]- [START_REF] Zhu | Analyzing dynamic performance of power systems over parameter space using normal forms of vector fields part II: Comparison of the system structure[END_REF] 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 [START_REF] Sanchez-Gasca | 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[END_REF].

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, [START_REF] Amano | Nonlinear stability indexes of power swing oscillation using normal form analysis[END_REF] 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 [START_REF] Tian | Analytical investigation of nonlinear interactions between voltage source converters interconnected to a transmission grid[END_REF] and [START_REF] Tian | Third-order based analytical investigation of nonlinear interactions between voltage source converters interconnected to a transmission grid[END_REF], some 3rd-order-based MNFs have been proposed [START_REF] Martínez | Perturbation analysis of power systems: Effects of second-and third-order nonlinear terms on system dynamic behavior[END_REF], [START_REF] Huang | Evaluation of the effect of modal interaction higher than 2nd order in small-signal analysis[END_REF] 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 [START_REF] Iooss | Topics in Bifurcation Theory and Applications[END_REF], [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF] or to define Nonlinear Modes of vibration and to build reduced-order models [START_REF] Jzquel | Analysis of non-linear dynamical systems by the normal form theory[END_REF]- [START_REF] Touzé | Non-linear normal modes for damped geometrically non-linear systems: Application to reduced-order modeling of harmonically forced structures[END_REF].

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é [START_REF] Poincar | Les Mthodes Nouvelles de la Mcanique Cleste[END_REF] 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 [START_REF] Nayfeh | A Method of Normal Forms[END_REF], [START_REF] Lamarque | An upper bound for validity limits of asymptotic analytical approaches based on normal form theory[END_REF] 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) [START_REF] Kundur | Power System Stability and Control[END_REF]. By subsituting the algebraic equations into the differential ones, one transforms those DAEs in a dynamical system, which can be written:

ẋ = f (x, u), ( 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:

ẋ = H1( x) + 1 2! H2( x) + 1 3! H3( x) + O(4)
(2) where Hq gathers the q th-order partial derivatives of f , i.e., for j = 1, 2, • • • , n, H1 j k = ∂f j /∂x k , H2 j kl = ∂ 2 f j /∂x k ∂x l , H3 j klm = ∂ 3 f j /∂x k ∂x l ∂x m and 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:

ẏ = Λy + F 2(y) + F 3(y) + O(4) (3)
supposed here to be diagonal, where the jth equation of ( 3) is:

ẏj = λ j y j + n k =1 n l=1 F 2 j kl y k y l + n p=1 n q =1 n r =1 F 3 j pq r y p y q y r +• • • (4)
λ j is the jth eigenvalue of matrix H1 and j = 1, 2, • • • , n. U and V are the matrices collecting the right and left eigenvectors of H1. U V = I, F

2 j = 1 2 n i=1 v j i U T H i 2 U and F 3 j pq r = 1 6 n i=1 v j i n k =1 n k =1 n m =1 H3 i klm u l p u m q u n r
where v j i is the element at j-th row and i-th column of V and u l p 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 σ j of the eigenvalues λ j = σ j + jω 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:

ẏ = Λy + F 2(y) (5)
where the jth equation of ( 5) writes :

ẏj = λ j y j + n k =1 n l=1 F 2 j kl y k y l (6)
To eliminate the second-order terms in [START_REF] Saha | Stability boundary approximation of a power system using the real normal form of vector fields[END_REF], the following secondorder NF transformation is applied: [START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. i. general theory and procedure[END_REF] or, in component form:

y = z + h2(z)
y j = z j + n k =1 n l=1 h2 j kl z k z l (8)
Use of ( 7) in (5) leads to:

ż = Λz -Dh2(z)Λz + Λh2(z) + F 2(z) + O(3) (9)
where Dh2 is the Jacobian matrix of h2 and the jth term of (9) takes the form of:

żj = λ j z j - n k =1 n l=1 [(λ k + λ l -λ j )h2 j kl -F 2 j kl ]z k z l + O(3) (10)
If one wants to eliminate all second-order terms from (9), transformation [START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. i. general theory and procedure[END_REF] must satisfy the following equation, for all z:

Dh2(z)Λz -Λh2(z) = F 2(z) (11)
which gives, for all j, k, l:

(λ k + λ l -λ j )h2 j kl -F 2 j kl = 0 (12)
Equation [START_REF] Sanchez-Gasca | 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[END_REF] shows that the computations of coefficients h2 j kl of the change of variables are not possible (or numerically difficult) if the conditions λ k + λ l -λ j ≈ 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 ω k + ω lω j = 0). It is important to note that this condition appears only when the eigenvalues share a commensurability relationship, as for example when ω 1 = 2ω 2 . Although it may be a rare case, [START_REF] Kakimoto | Instability of interarea oscillation mode by autoparametric resonance[END_REF] 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 j kl is given by [START_REF] Sanchez-Gasca | 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[END_REF]:

h2 j kl = F 2 j kl λ k + λ l -λ j ( 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:

ż = Λz (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 [START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. i. general theory and procedure[END_REF] 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 [START_REF] Lin | Investigation of modal interaction and its effects on control performance in stressed power systems using normal forms of vector fields[END_REF]- [START_REF] Zhu | Analyzing dynamic performance of power systems over parameter space using normal forms of vector fields part II: Comparison of the system structure[END_REF] opens the era to apply Normal Forms analysis in studying nonlinear dynamics in power system. Its effectiveness has been shown in many examples [START_REF] Lin | Investigation of modal interaction and its effects on control performance in stressed power systems using normal forms of vector fields[END_REF], [START_REF] Saha | Stability boundary approximation of a power system using the real normal form of vector fields[END_REF]- [START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. ii. case studies[END_REF], [START_REF] Zhu | Analyzing dynamic performance of power systems over parameter space using normal forms of vector fieldspart I: Identification of vulnerable regions[END_REF], [START_REF] Zhu | Analyzing dynamic performance of power systems over parameter space using normal forms of vector fields part II: Comparison of the system structure[END_REF], [START_REF] Barocio | Analysis of factors affecting power system normal form results[END_REF], [START_REF] Liu | Assessing placement of controllers and nonlinear behavior using normal form analysis[END_REF], to approximate the stability boundary [START_REF] Saha | Stability boundary approximation of a power system using the real normal form of vector fields[END_REF], to investigate the strength of the interaction between oscillation modes [START_REF] Lin | Investigation of modal interaction and its effects on control performance in stressed power systems using normal forms of vector fields[END_REF], [START_REF] Starrett | Nonlinear measures of mode-machine participation [transmission system stability[END_REF], [START_REF] Barocio | Analysis of factors affecting power system normal form results[END_REF], to dealt with a control design [START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. i. general theory and procedure[END_REF], [START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. ii. case studies[END_REF], to analyze a vulnerable region over parameter space and resonance conditions [START_REF] Zhu | Analyzing dynamic performance of power systems over parameter space using normal forms of vector fieldspart I: Identification of vulnerable regions[END_REF], [START_REF] Zhu | Analyzing dynamic performance of power systems over parameter space using normal forms of vector fields part II: Comparison of the system structure[END_REF] and to optimally place controllers [START_REF] Liu | Assessing placement of controllers and nonlinear behavior using normal form analysis[END_REF].

This method is well summarized in [START_REF] Sanchez-Gasca | 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[END_REF], 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 [START_REF] Betancourt | A real normal form approach to the study of resonant power systems[END_REF] proposes an extension of the method to deal with resonant cases. It has become a mature computational tool and other applications have been developed [START_REF] Martínez | Normal form analysis of complex system models: A structure-preserving approach[END_REF], [START_REF] Davarani | Nonlinear modal analysis of interaction between torsional modes and SVC controllers[END_REF], or are still emerging [START_REF] Lomei | 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[END_REF].

2) Consideration of Cubic Terms: Although the secondorder-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 [START_REF] Amano | Nonlinear stability indexes of power swing oscillation using normal form analysis[END_REF] and new criteria to design the system controllers in order to improve the transfer capacity of transmissions lines have been established in [START_REF] Amano | A new PSS parameter design using nonlinear stability analysis[END_REF].

Including third-order terms, (3) becomes:

ẏ = Λy + F 2(y) + F 3(y) (15)
The j-th state equation being:

ẏj = λ j y j + n k =1 n l=1 F 2 j kl y k y l + n p=1 n q =1 n r =1
F 3 j pq r y p y q y r (16) Applying the 2nd Order NF transformation y = z + h2(z), the normal dynamics reads:

ż = Λz + DF 2(z)h2(z) + F 3(z) + O(4) (17)
where: 1) DF 2(z)h2(z) are third order terms coming from the 2nd Order NF transformation used to cancel the 2nd order terms; 2) F 3(z) are the original 3rd order terms from the system (3). Neglecting terms of order higher than 3, for the jth variable [START_REF] Huang | Evaluation of the effect of modal interaction higher than 2nd order in small-signal analysis[END_REF] results in:

żj = λ j z j + n p=1 n q =1 n r =1 C j pq r z p z q z r ( 18
)
where

C j pq r = n l=1 (F 2 j pl + F 2 j lp )h2 p qr + F 3 j pq r
To simplify as much as possible the third order terms, the following third-order NF transformation is applied to [START_REF] Huang | Evaluation of the effect of modal interaction higher than 2nd order in small-signal analysis[END_REF] [START_REF] Martínez | Perturbation analysis of power systems: Effects of second-and third-order nonlinear terms on system dynamic behavior[END_REF]:

z = w + h3(w) ( 19
)
where h3 is a polynomial of w containing only 3rd-order terms.

Applying transformation [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF] to [START_REF] Huang | Evaluation of the effect of modal interaction higher than 2nd order in small-signal analysis[END_REF], it leads to:

ẇ = Λw + DF 2(w)h2(w) + Λh3(w) -Dh3(w)Λw + F 3(w) + O(5) (20)
It should be emphasized that there are neither terms of order 2 nor terms of order 4 in [START_REF] Jzquel | Analysis of non-linear dynamical systems by the normal form theory[END_REF], and terms of order 3 are: 1) DF 2(w)h2(w) that comes from the 2nd-order NF transformation used to cancel the 2nd-order terms; 2) Λh3(w) -Dh3(w)Λw that comes from the use of transformation [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF] in order to cancel the 3rd-order terms; 3) F 3(w), the original 3rd order terms of system [START_REF] Martínez | Perturbation analysis of power systems: Effects of second-and third-order nonlinear terms on system dynamic behavior[END_REF]. For the j-th variable, (20) can be written as:

ẇj = λ j w j + n p=1 n q =1 n r =1 [C j pq r -(λ p + λ q + λ r -λ j )h3 j pq r ] × w p w q w r + O(5) (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 [START_REF] Touzé | Hardening/softening behaviour in non-linear oscillations of structural systems using non-linear normal modes[END_REF] shows that third order terms can be eliminated with setting h3 computed by:

h3 j pq r = F 3 j pq r + n l=1 (F 2 j pl + F 2 j lp )h2 p qr λ p + λ q + λ r -λ j ( 22
)
Computation of h3 are impossible if the condition λ p + λ q + λ r -λ j ≈ 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 the real axis).

If we consider a weak-damped oscillating mode, composed of two conjugated poles λ 2l and λ 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 j j 2l2l-1 . As a consequence, term w j w 2l w 2l-1 cannot then be eliminated and must be kept in the normal dynamics [START_REF] Lamarque | An upper bound for validity limits of asymptotic analytical approaches based on normal form theory[END_REF], [START_REF] Dulac | Solutions d'un systéme d'équations différentielles dans le voisinage de valeurs singuliéres[END_REF].

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 ∈ M ).

Apart from internal resonances due to commensurability relationships between eigenvalues, as for example when ω 1 = 3ω 2 , ω 1 = 2ω 2 + ω 3 or ω 1 = ω 2 + ω 3 + ω 4 , all the other third-order terms are eliminated by third-order transformation [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF].

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

ẇj = λ j w j + M l=1 c j 2l w j w 2l w 2l-1 , j ∈ M (23) ẇj = λ j w j , j / ∈ M ( 24
)
where w 2l-1 is the complex conjugate of w 2l and coefficients c j 2l are defined as:

c j =2k 2l = C 2k (2k )(2l)(2l-1) + C 2k (2k )(2l-1)(2l) + C 2k (2l)(2k )(2l-1) + C 2k (2l)(2l-1)(2k ) + C 2k (2l-1)(2l)(2k ) + C 2k (2l-1)(2k )(2l) c 2l 2l = C 2l (2l)(2l)(2l-1) + C 2j (2l)(2l-1)(2l) + C 2l (2l-1)(j )(j ) (25
) Equation [START_REF] Poincar | Les Mthodes Nouvelles de la Mcanique Cleste[END_REF] shows that considering third-order terms in the Normal Dynamic changes the way of making the stability analysis. As shown in [START_REF] Amano | Nonlinear stability indexes of power swing oscillation using normal form analysis[END_REF], [START_REF] Amano | A new PSS parameter design using nonlinear stability analysis[END_REF], 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 [START_REF] Amano | Nonlinear stability indexes of power swing oscillation using normal form analysis[END_REF], [START_REF] Amano | A new PSS parameter design using nonlinear stability analysis[END_REF], 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 :

żj = λ j z j + M l=1 c j 2l z j z 2l z 2l-1 , j ∈ M ( 26
)
For the non-oscillatory modes, the normal dynamics are the same as expressed by [START_REF] Tian | Analytical investigation of nonlinear interactions between voltage source converters interconnected to a transmission grid[END_REF].

In [START_REF] Martínez | Perturbation analysis of power systems: Effects of second-and third-order nonlinear terms on system dynamic behavior[END_REF], [START_REF] Huang | Evaluation of the effect of modal interaction higher than 2nd order in small-signal analysis[END_REF], a 3rd-order NF transformation is proposed where the terms DF 2(w)h2(w) of ( 20) are not taken into account, leading to coefficients given by:

h3 j pq r = F 3 j pq r λ p + λ q + λ r -λ j ( 27
)
Since some 3rd order terms are omitted, the accuracy of the proposed normal dynamics is worse than the one based on coefficients given by [START_REF] Touzé | Non-linear normal modes for damped geometrically non-linear systems: Application to reduced-order modeling of harmonically forced structures[END_REF].

Moreover, the considered normal dynamics are: ẇj = λ j w j [START_REF] Liu | Assessing placement of controllers and nonlinear behavior using normal form analysis[END_REF] 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 [START_REF] Poincar | Les Mthodes Nouvelles de la Mcanique Cleste[END_REF] and [START_REF] Nayfeh | A Method of Normal Forms[END_REF] have to be reconstructed with the z coordinates as:

z j = w j + n p=1 n q =1 n r =1 h3 j pq r w p w q w r ( 29
)
It can be seen from the reconstruction that z j is composed of two sets. A first set with coefficients h3 j pq r that can be used to quantify the 3rd order modal interactions [START_REF] Huang | Evaluation of the effect of modal interaction higher than 2nd order in small-signal analysis[END_REF] and a second set where w j is influenced by c j 2l , making then possible to use it for stability analysis [START_REF] Amano | Nonlinear stability indexes of power swing oscillation using normal form analysis[END_REF], [START_REF] Amano | A new PSS parameter design using nonlinear stability analysis[END_REF].

Based on the Normal Forms, h3 j pq r and c j 2l 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 [START_REF] Betancourt | A real normal form approach to the study of resonant power systems[END_REF], 3rd order oscillations are modeled by h3 j pq r w p w q w r . Then, a third order Modal Interaction index MI3 j pq r can be defined as:

MI3 j pq r = |h3 j pq r w 0 p w 0 q w 0 r | |w 0 j | (30)
which indicates the participation of the 3rd order modal interaction with the frequency ω p + ω q + ω r in the mode j. ω 0 are the initial values of the variables at the moment of interest. Introducing MI3 j pq r 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, ẇj = λ j w j , the stability is determined by inspecting λ j . For the oscillatory modes where j ∈ M , the normal dynamics are:

ẇj = λ j w j + M l=1 c j 2l w j w 2l w 2l-1 , (31)
= λ j + M l=1 c j 2l |w 2l | 2 w j ( 32
)
with |w 2l | the magnitude of w 2l .

A stability interaction index can then be defined as:

SII j 2l = c j 2l |w 0 2l | 2 , j ∈ M ( 33
)
and, remarking that

λ j + M l=1 c j 2l |w 0 2l | 2 = σ j + real M l=1 SII j 2l +j(ω j +imag(SII j 2l )) (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 MI3 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 [START_REF] Poincar | Les Mthodes Nouvelles de la Mcanique Cleste[END_REF] and [START_REF] Nayfeh | A Method of Normal Forms[END_REF] in [START_REF] Sanchez-Gasca | 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[END_REF], T r3 is the ratio of the time constant associated to a combination of modes over a dominant mode.

T r3 =

Time constant for combination of modes (λ p + λ q + λ r ) Time constant for dominant mode (λ j )

A small T r3 indicates a significant presence of a mode combination. For example, if T r3 = 1, it means that a nonlinear interaction decays at the same speed as a dominant mode. If T r3 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 × T r3, MI3 × T r3 tends to reveal the presence of a persistent modal interaction.

D. Stability Assessment

As indicated by [START_REF] Davarani | Nonlinear modal analysis of interaction between torsional modes and SVC controllers[END_REF], it is both the eigenvalues λ j and the stability indexes SII j 2l that give information concerning the system stability. Although λ j keeps constant, M l=1 c j 2l |w 2l | 2 w j decays as time runs. For example, if SII j 2l is large, but SII j 2l × T r3 j j (2l)(2l-1) is small, the 3rd-order terms may stabilize the system at the beginning of the transient but a longterm instability can exist. Therefore, the time constant associated to SII j 2l must be taken into account when assessing the overall stability. An index of stability assessment can then be defined as:

SI j = σ j + M l=1 real(SII j 2l ) × T r3 j j (2l-1)(2l) , j ∈ M (36)
with σ being the real value of an eigenvalue (λ j = σ j + jω 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 [START_REF] Amano | Nonlinear stability indexes of power swing oscillation using normal form analysis[END_REF], 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 σ 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 [START_REF] Dobson | Scaling of normal form analysis coefficients under coordinate change[END_REF].

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 [START_REF] Kundur | Power System Stability and Control[END_REF] and for the validation of 2nd order Normal Form analysis [START_REF] Sanchez-Gasca | 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[END_REF].

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 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 P SS 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 Ka 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 [START_REF] Saha | Stability boundary approximation of a power system using the real normal form of vector fields[END_REF][START_REF] Starrett | Nonlinear measures of mode-machine participation [transmission system stability[END_REF] and [START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. i. general theory and procedure[END_REF][START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. ii. case studies[END_REF], 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.

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 order 3, contributing to a better modeling of the frequency variation of the oscillations as recently formulated in [START_REF] Wang | Formulation and characterization of power system electromechanical oscillations[END_REF] and [START_REF] Wang | Properties of the frequency-amplitude curve[END_REF].

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 [START_REF] Saha | Stability boundary approximation of a power system using the real normal form of vector fields[END_REF][START_REF] Starrett | Nonlinear measures of mode-machine participation [transmission system stability[END_REF] and [START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. i. general theory and procedure[END_REF][START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. ii. case studies[END_REF] 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 [START_REF] Amano | Nonlinear stability indexes of power swing oscillation using normal form analysis[END_REF]. 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).

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 0 17 , 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 DF 2h2, h3 j pq r is small and leads to a too modest prediction of the 3rd order interaction MI3 j pq r = 0.056 = 5.6%. A deeper analysis can be made using NF methods compared to FFT analysis. For example, MI3 17 17,26,27 = 0.45 indicates a strong nonlinear interaction, however, since T r = 6.7341 × 10 -4 and MI3 × T r = 3.03 × 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 j 2l . Among all those coefficients, coefficients SII j 18 are predominant and T r3 indicates a long-time influence as listed in Table VI. Although the real part of SII j 18 is too small to contribute to the stability compared to λ j , its imaginary part indicates that the there is a frequency shift added to the eigenfrequency.

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 ω j .

4) Stability Assessment: Seen from Table VII, the nonlinear interactions can enhance (e.g MI3 5 5,5,6 ) or weaken the stability (e.g MI3 j 7,17,18 ). 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 Fig. 5. 16-Machine 5-area test system [START_REF] Singh | IEEE PES task force on benchmark systems for stability controls report on the 68-bus, 16-machine, 5-area system[END_REF]. method 3-3-3 is more rigorous and comprehensive than method 3-2-3. As SI 5 0, SI 7 0, modes (5, 6) and [START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. i. general theory and procedure[END_REF][START_REF] Jang | Effect of nonlinear modal interaction on control performance: Use of normal forms technique in control design. ii. case studies[END_REF] 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 [START_REF] Singh | IEEE PES task force on benchmark systems for stability controls report on the 68-bus, 16-machine, 5-area system[END_REF], 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 [START_REF] Singh | IEEE PES task force on benchmark systems for stability controls report on the 68-bus, 16-machine, 5-area system[END_REF]. The realistic parameters and well validated Simulink models can be found in [START_REF] Singh | IEEE PES task force on benchmark systems for stability controls report on the 68-bus, 16-machine, 5-area system[END_REF], [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 [START_REF] Singh | IEEE PES task force on benchmark systems for stability controls report on the 68-bus, 16-machine, 5-area system[END_REF].

When placing PSSs on G1 to G12 (the maximum number of possible PSSs) all the local modes are damped, and 3 intermodes are poorly damped. No more information is available from the linear analysis to damp the inter-area modes [START_REF] Singh | IEEE PES task force on benchmark systems for stability controls report on the 68-bus, 16-machine, 5-area system[END_REF]. When a three-phase fault is applied to G13 and cleared in 0.25s, the machines exhibit nonlinear inter-area oscillations and finally 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), Non-linear All 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 [START_REF] Liu | Assessing placement of controllers and nonlinear behavior using normal form analysis[END_REF].

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 λ and matrices F 2 and F 3, 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 [START_REF] Sanchez-Gasca | 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[END_REF], [START_REF] Barocio | Analysis of factors affecting power system normal form results[END_REF] 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 [START_REF] Mikofski | Newton-Raphson solver[END_REF], 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 [START_REF] Barocio | Analysis of factors affecting power system normal form results[END_REF] 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 F 2 and F 3 can be extremely tedious. The problem is that matrices F 2 and F 3 are composed of complex values. Using the same Matlab function to calculate matrices A, H2, H3, Λ, F 2 and F 3 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 [START_REF] Martínez | Normal form analysis of complex system models: A structure-preserving approach[END_REF], and 3rd order generator as method 3-2-3 used in [START_REF] Amano | Nonlinear stability indexes of power swing oscillation using normal form analysis[END_REF].

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 σ j =0 . In this case, method 3-3-1 fails to be applied, since h3 j j 2l-12l = ∞, (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 [START_REF] Amano | Nonlinear stability indexes of power swing oscillation using normal form analysis[END_REF], [START_REF] Amano | A new PSS parameter design using nonlinear stability analysis[END_REF] [see [START_REF] Amano | Nonlinear stability indexes of power swing oscillation using normal form analysis[END_REF] and [START_REF] Touzé | Non-linear normal modes for damped geometrically non-linear systems: Application to reduced-order modeling of harmonically forced structures[END_REF]] will be wrongly predicted, as the stability boundary will be calculated as R j = -real(λ j ) real(c j = 2 k 2 k ) = 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

Fig. 1 .

 1 Fig. 1. IEEE 4 machine test system: 2-area 4-generator.

Fig. 3 .

 3 Fig. 3. FFT Analysis of different NF approximations: Case 1.

Fig. 4 .

 4 Fig. 4. Comparison of different NF approximations: Case 2.

Fig. 6 .

 6 Fig. 6. FFT analysis of generator rotor angles when G1 to G12 are equipped with PSSs.

TABLE I LIST

 I OF THE MNFS PRESENTED IN THIS WORK

	Method	Taylor's	NF	Normal	nth order	Related
		series	Trans.	Form	NF Terms	Equation
	2-2-1	2nd order order 2	order 1	/	(14)
	3-2-3S	3rd order	order 2	order 3	S	(26)
	3-3-1	3rd order	order 3	order 1	/	(28)
	3-3-3	3rd order	order 3	order 3	/	(23)(24)

TABLE VIII INTER

 VIII -AREA MODES OF THE 16 MACHINE 5 AREA SYSTEM WITH PSSS PRESENT ON G6 AND G9

	M o d e	1	2	3	4
	Damping ratio	3.35%	0.55%	1.58%	2.76%
	Frequency (Hz)	0.7788 0.6073 0.5226 0.3929
	Fig. 7. FFT analysis of generator rotor angles when only on G6 and G9 are
	equipped with PSSs.

TABLE IX SEARCH

 IX FOR INITIAL CONDITIONS IN THE NF COORDINATES

	Method		Case 1	Case 2
		Iterations	Resolution Iterations	Resolutions
	3-3-3		12	1.294e-12	17	1.954e-6
	2-2-1/3-2-3	3	1.799e-11	6	1.838e-06
	3-3-1		7	1.187e-09	11	7.786e-09
				TABLE X
	PERFORMANCE EVALUATION OF THE STUDIED NF APPROXIMATIONS
	Method Resolution	Order of Modal	Transient	Applicability to
			Interaction	Stability	3rd Order Modes
	Linear 2-2-1 3-2-3S 3-3-1 3-3-3	O(2) O(3) O(3) O(3) O(5)		1 st 2 nd 2 nd 3 rd 3 rd	Linear Linear Non-linear Linear	None None Oscillatory Non-oscillatory

TABLE XI GENERATOR

 XI DATA IN PU ON MACHINE BASE

	Parameter	Value	Parameter	Value
	R a	0.0025	x l	0.002
	x d	1.80	τ d 0	8.0
	x q	1.70	τ q 0	0.40
	x d	0.3	τ d 0	0.0
	x q	0.3	τ d 0	0.0
	x d	0.0	M V A b ase	900
	x q	0.0	H	6.5

TABLE XII POWER

 XII FLOW DATA FOR CASE 1, CASE 2 IN STEADY STATE AFTER THE FAULT IS CLEARED: LOAD (CALCULATED BY MATLAB SIMPOWER)

	Load	Voltage	Z shu nt
	Bus	(pu/230kv)
	B5 B6	1.029∠19.75 • 0.987∠ -47.37 •	P = 665.84M W , Q = -350.82M var P = 1443.38M W , Q = -350.02M var
	Generator	Bus	Terminal Voltage	P	Q
		Type	(pu/20kv)	(MW)	(Mvar)
	G1 G2 G3 G4	P V swingbus P V P V	1.00∠14.02 • 1.000∠0.0 • 1.00∠ -58.17 • 1.00∠ -66.11 •	747.50 406.50 539.25 525.00	49.94 140.74 113.99 143.14