Harmonic Solution of Higher-Dimensional Second Order Kuramoto Oscillator Network

Obtaining the solution of Kuramoto oscillator model is an effective way to understand the mechanism behind the dynamical behaviors it represents. However, it is not an easy work to get, even the approximation solution for the second order Kuramoto oscillator network because of its complex coupling interaction and strong nonlinearity. Hence, this paper first establishes a new polynomial mathematical model for the second order Kuramoto oscillator network, which is partly decoupled in a polynomial fashion. Then, a quadratic harmonic balance method (HBM) is developed to get the approximate solution for quadratic differential algebraic equation (QDAE) system. An equivalent quadratic technique is introduced to convert the polynomial Kuramoto oscillator model into a QDAE system aiming to obtain its harmonic solution. Finally, the suggested harmonic solution is utilized to explain the morphological characteristics of the stable motion of oscillators, and the relative motions for the oscillators in Kuramoto network is concluded by a set of nonlinear state output equations. It is also found that the distance between the maximum phase angle and the controlling unstable equilibrium point (UEP) can be used as a stability margin for the oscillator network. A single-machine infinite bus (SMIB) power system, an IEEE 3-machine 9-bus power system and an IEEE 10-machine 39-bus power system are employed as different dimension second order Kuramoto oscillator network cases to demonstrate the validation and accuracy of the proposed method and analysis.


I. INTRODUCTION
The Kuramoto model has been introduced in 1975 to describe the motion of coupled particles [1].In the last forty years, the second order Kuramoto oscillator model has been extended to study problems originated from many disciplines [2], [3], representing the dynamics of electron oscillators [4] in electric engineering, Josephson junction [5] and pendulum equation [6] in physics, synchronous machines [7] in power system community, and so forth.Different from the classical first order Kuramoto oscillator model, the second order Kuramoto model exhibits much more complicated behaviors because of the way of its complicated coupling interaction, posing a challenge for developing method to study the dynamical characteristic of this system it represents.
The associate editor coordinating the review of this manuscript and approving it for publication was Yifan Zhou.
As is well-known, it is very difficult to find closed-form solutions for second order Kuramoto models [8], [9], [10], [11].Hence, it is natural to seek the approximate solution in order to understand the dynamic behavior of the second order Kuramoto oscillator network.In history, several attempts have been reported to get the analytical characteristic of the second order Kuramoto oscillator network.Some early results using asymptotic expansion methods were introduced in [12], [13], and [14] to approximate the stable manifolds of the UEP.Some approximate formulas were introduced in [15], [16], and [17] to represent stability region boundaries which are essentially composed of stable manifold of unstable equilibrium point (UEP) of the second order Kuramoto oscillator model.However, due to the inherent nonlinearity of second order Kuramoto oscillator network, the standard perturbation methods such as normal form method [13], [14], Krylov-Bogoli-Mitropolsky method [18], Lindstedt-Poincaré method [19], multi-time scales method [20], do not perform satisfactorily.The results mentioned above focus on the behavior in the neighborhood of an UEP.Little attention has been paid to the behavior of Kuramoto oscillators inside the stability region of a stable equilibrium point (SEP).Furthermore, it has been recognized that the swing behavior of pendulum systems is typically quasi-periodical [21].This indicates that the second order Kuramoto dynamic system has a natural low-pass filtering characteristic, revealing the potential of harmonic balance method (HBM) in studying the motion of Kuramoto oscillators.
Recently, the incremental HBM has been applied to estimate the limit cycle of swing equation in power system community [22], predicting the homoclinic bifurcation curve effectively, showing the potential of the HBM in dealing with strongly nonlinear systems.In the community of mechanical engineering, a numerical approach to obtain the quasi-periodic vibrations by HBM is developed in [23], demonstrating the satisfactory performance of HBM when the models are experiencing nonlinearities induced by strong engine order excitations.The periodic response of the systems subjected to harmonic forces is obtained by employing HMB method in [24], revealing the dynamic characteristic and vibration transmission behavior of oscillators system.HBM and Galerkin methods are used together to camfollower system to estimate the steady-state response [25], which can predict the fundamental and sub-harmonic solutions accurately and efficiently.These results show the application potential of HBM in dealing with different nonlinear system, and have been extensively studied in various field.For more detailed technique and progress of HBM, the readers are referred to reference [26], [27], and [28], respectively.
Different from previous works focusing the analytical characteristic around the UEP of Kuramoto model mentioned before, this paper investigates the harmonic solution of higher-dimensional second order Kuramoto oscillator model within stability region of a SEP.In other words, the system trajectories reside inside stability region of a SEP.Firstly, a new polynomial model of Kuramoto oscillators network is established via Taylor expansion formula.Then, a multi-frequency quadratic harmonic balance method (QHBM) inspired by the result in [29] and [30], is developed.The method can find the harmonic solution by providing only some coefficient matrices related to the quadratic differential algebraic equation (QDAE) system.Thereby, based on the suggested harmonic solution, the relationships between the eigenvalues and dynamics of the Kuramoto oscillator model are examined, the relative motion equation system between different oscillators is established, and a new phase-domain stability margin index for higher-dimensional oscillator network is proposed based on the celebrated stability region theory.Finally, some power systems data are employed to show the effectiveness of the proposed method.
To sum up, the original contributions of this paper are threefold: 1) Model.A partly decoupled polynomial equation system for the second order Kuramoto model is derived in detail, revealing the morphological characteristics of stable motion of oscillators.
2) Methodology.A new QHBM based on QDAE system is developed systematically to obtain the steady-state or periodic response for zero-damping systems.
3) Analysis.The analysis of solution structure, relative motion of oscillators and large disturbance stability is presented based the obtained harmonic solution.
The rest of the paper is organized as follows.In Section II, a new polynomial model for the second order Kuramoto oscillator network is established, and the structural properties of the derived polynomial equations are presented.Section III gives the fundamental basis of the QHBM rigorously, which provide the details of forming harmonic algebraic equation (HAE).In Section IV, an equivalent quadratic conversion technique is introduced to obtain the harmonic solution of second order Kuramoto oscillator network.Section V investigates the analytical characteristic of second order Kuramoto oscillator network based on the harmonic solution.Section VI shows the effectiveness of the proposed method and analysis using several power system test cases as examples.Section VII concludes this paper.

II. A NEW POLYNOMIAL MODEL FOR SECOND-ORDER KURAMOTO OSCILLATOR NETWORK
In this section, a polynomial model for second order Kuramoto oscillator network is derived firstly.The model is formulated in matrix form so it can accommodate Kuramoto model with arbitrary number of oscillators.An eigen transformation is introduced later to decouple the polynomial model into two subsystems.

A. A TAYLOR EXPANSION FORMULATION FOR SECOND ORDER KURAMOTO NETWORK WITH ANY NUMBER OF OSCILLATORS
Consider the second order Kuramoto oscillator network on a graph G ≜ (N, E) as follows, (1) where N = {1, 2, 3, . . ., n, n + 1} and E = N ⊗ N are the node set and edge set of G, respectively.The symbol ⊗ represents Cartesian product.Here δ i is the phase angle of oscillator i. a ij is the stiffness of elastic spring between oscillator i and oscillator j. γ ij is the phase shift between oscillator i and oscillator j.M i and D i are the inertial time constant and damping coefficient of oscillator i,respectively.i is the driving force of oscillator i.From physical perspective, the second order Kuramoto network is shown in FIGURE 1.The Taylor expansion formula is used to approximate the sine functions in the right-hand side (RHS) of (1) as follows, where ϵ 3 = −1/6, ϵ 5 = 1/120, ϵ 7 = −1/5040.Furthermore, the equilibrium point is moving to the origin by introducing coordinate transformation y = δ − δ, here δ indicates a SEP of (1).Substituting the (2) into (1), the RHS of (1) can be written as follows, RHS Plugging (3) into (1), the polynomial Kuramoto model in y coordinate system is obtained as follows, where coefficients in (4) are shown as follows, In order to write RHS of (4) into matrix form, the following coefficient matrix P is defined as follows, where y ij = y i − y j , 1 ≤ i < j ≤ n + 1 is the angle difference in y coordinate system.Hence, the matrix form of (4) can be represented as follows, where the inertial and damping matrix, respectively.The coefficient matrixes A {2,3,•••,7} ∈ R (n+1)×[(n+1)n / 2] are corresponding to (4), separately.Specifically, A 1 can be written as follows, where A 1 ∈ R (n+1)×(n+1) is the Laplacian matrix of graph G Therefore, A 1 can be formulated via algebraic graph operation equivalently, i.e., n+1) is the adjacency matrices of graph G, and 1 n is vector with 1 being the entries with appropriate dimension.
Notice that the matrix P has structural characteristic because of the constraint imposed by (6).Without loss of generality, an example with n = 3 is illustrated as follows, It can be seen that P • 1 = 0 from (9), which means that the summation of any row of P is equal to 0. Denoting the symbol L indicating the set of zero row sum matrix, then we have P, A 1 ∈ L.

An eigen transformation y
where n are the n nonzero eigenvalues of the matrix M −1 A 1 because of its zero-row sum property.Considering P ∈ L, one obtains, where α ∈ R n α ×n , n α = (n + 1)n/2.Plugging transformation y = U z into (7) with the aid of ( 10) and ( 11), a partly decoupled polynomial model is obtained as follows, where . Suppose that there is uniform damping in Kuramoto network [31], i.e., D i /M i = µ, then ( 12) can be written as follows, where g nl : R n → R n , f nl : R n → R are all polynomial functions.It can be observed that ( 13)-(I) is completely decoupled from variable z n+1 , and ( 13)-(II) completely depends on the main subsystem variable z 1 , revealing the partly decoupled characteristic.

III. A QUADRATIC HARMONIC BALANCE METHOD
In this section, a QHBM for QDAE systems is developed systematically.The principle of QHBM is reviewed in Section III-A, and the detailed procedure on how forming HAE of QDAE system to obtain harmonic solution is presented from Section III-B to III-D.
A. THE PRINCIPLE OF QHBM [29], [30] Consider the following initial value problem of QDAE system, where X ∈ R n e is unknowns.m : R n e → R n e is a linear vector function acting on the derivative of X. l : R n e → R n e and q : R n e × R n e → R n e are linear vector function and quadratic bilinear function acting on X, respectively.Suppose that the system ( 14) has a periodic solution with prior knowledge, and the form of this periodic solution is assumed as follows, where a ∈ R n e ×1 is a constant term, and b k , c k ∈ R n e ×n e are the harmonic coefficients, and ω T = [ω 1 , ω 2 , . . ., ω n e ] T is the frequency vector of harmonic solution, and the symbol H is the maximum order of harmonic.Denote the unknown vector appeared in (15) as follows, where the symbol row(b) indicates the column vector sorted by row vector of matrix b consecutively, and the unknown variable is U ∈ R (n e ×n e )×2H +n e .Inserting the period solution ( 15) into ( 14), each harmonic is collected consecutively according to ascending order of harmonics.After some manipulations, it is not difficult to observe that the unknown vector satisfies n e × n e × 2H + n e HAE as follows, where symbols M(•), C, L(•), Q(•, •) are all operators acting on the variable U, which only depends on the details of the QDAE system (14) and harmonic order H . Adding equations generated by the initial value in (14), one obtains n e equations as follows, It can be observed that the dimension of [U T , ω T ] is equal to 2H ×n 2 e +2n e , which is same with the number of equations simultaneous (17) and (18).Hence, some standard iterative method can be used to solve the algebra equations.

B. LINEAR OPERATOR M (•), L (•) AND CONSTANT OPERATOR C
The linear operator m(•) in ( 14) can be written as follows, where D ∈ R n e ×n e is a specific constant matrix depending on the detail of operator m(•) in ( 14).Specifically, ( 19) can be formulated as follows, Hence, M(U) in ( 17) can be formulated according to the order of harmonic as follows, where M(U) ∈ R (n e ×n e )×2H +n e .The constant operator C in ( 17) is the zeros order harmonic corresponding to the constant term in (14).Hence, we have, where C ∈ R (n e ×n e )×2H +n e .The linear operator l(•) in ( 14) can be written into matrix form as follows, Lc k sin(kωt), (23) where L ∈ R n e ×n e is a specific constant matrix corresponding to the detail of operator l(•) in (14).Thus, L(U) in ( 17) can be obtained as follows, where L(U) ∈ R (n e ×n e )×2H +n e .
The quadratic operator q(•, •) in ( 14) can be written into the matrix form as follows, where A i ∈ R n e ×n e , i ∈ {1, 2, . . ., n e } indicate specific constant matrixes corresponding to the detail of operator q i in ( 14).Substituting the harmonic solution ( 15) into (25), one obtains, where q k , k ∈ {1, 2, . . ., n e } has general expression as follows, In the four parts of ( 27), the third item in (27) can only generate constant bias as follows, The first and second item in (27) can only generate integer harmonic as follows, The general expression of the fourth term in ( 27) can be formulated as follows, where the harmonic crossing term (HCT) mainly produces non-integer sub-harmonics.
Because of the natural low-pass filtering characteristic of second order Kuramoto model, non-integer sub-harmonics only have limited influence to the approximate solution.
Thus, the forms by ignoring the HCT of the fourth term in (27) according to (30) can be obtained as follows, where q 0,k ∈ R, q cj,k ∈ R n e , and q sj,k ∈ R n e are expressed as follows, where the symbol diag(•) represent to take matrix diagonal operation.Combining ( 28), (29), and (31), q(X, X) in ( 25) can be formulated as follows, where each coefficient in (35) of is demonstrated as follows, Hence, the expression of Q(U, U) in (17) according to (25) can be formulated as follows,

IV. HARMONIC SOLUTION OF SECOND ORDER KURAMOTO OSCILLATOR MODEL
In this section, the pure differential equations described in section II is converted into a QDAE system, technically, by introducing auxiliary variables in Section IV-A and the coefficients matrices related to QDAE is derived later.

A. THE QDAE OF SECOND ORDER KURAMOTO OSCILLATOR MODEL
Considering the weak damping property of power system [32], the damping coefficient in ( 13) is set to zero, i.e., µ = 0.A set of auxiliary variables is introduced as follows, Denoting z ′ 1 = ϖ , ( 13)-(I) can be formulated as a QDAE system as follows, where , n e = 2n + 3N α , and the symbol • denotes the Hadamard product of vectors.
So far we get all the coefficient matrices A k , L, D, C of the (40), thus the HAE can be formulated according to (38) aiming to obtain the harmonic solution of ( 13)-(I).

V. ANALYSIS OF THE SECOND ORDER KURAMOTO OSCILLATOR NETWORK
This section analyses the solution structure of Kuramoto oscillator network via obtained harmonic solution.Then the relative motion equation of oscillators in Kuramoto network is concluded by a set of nonlinear state output equations.Finally, a phase-domain margin index is proposed to assess the stability of oscillator system quantitatively.

A. THE SOLUTION STRUCTURE OF THE KURAMOTO OSCILLATOR NETWORK
For ease of analysis, supposing that the order of harmonics is equal to 1, i.e., H = 1, and assuming that the HAE (38) can be solved numerically.The harmonic solution of the second order Kuramoto oscillator network is simplified as follows, where a is constant vector, M bc is the amplitude matrix, and ϕ is the initial phase vector.Higher harmonic terms (HHT) represent those harmonics with order greater than 1.Plugging (52) into ( 13)-(II) with zero-damping assumption, one can obtain as follows, where k 2 , a 2 , ϕ 20 , and ω are the constant determined by the nonlinear function f nl (•).
Neglecting HHT and integrating (53) the first-order harmonic approximate solution of (53) can be obtained as follows, where C 0 , C 1 are the integral constants determined by the initial value condition in (53).The analytical solution in y coordinate system can be formulated by using the inverse transformation as follows, (55) It can be seen from ( 55) that the solution of second order Kuramoto oscillator network can be decomposed into two parts.One part is the multi-harmonic oscillation motion y o , that is non-zero eigenvalue part, which makes the phase angles of the whole network oscillate with each other, and the other is the acceleration motion component y s , i.e., zero eigenvalue part, which pushes the phase angles of the whole network forward.
In particular, if one sets ϵ 3 = ϵ 5 = ϵ 7 = 0, the nonlinear system ( 13) degenerates into a linear system as follows, And the approximate solution of (55) becomes exact solution of (56) as follows, Thus, previous discussion reveals that the zero eigenvalue of A 1 is related to the uniform motion of all oscillators, while the non-zero eigenvalues basically dedicate to the vibration among oscillators.

B. THE RELATIVE MOTION EQUATION OF OSCILLATORS IN KURAMOTO NETWORK
Denoting the transformation matrix U mentioned in (10) as , the phase angle differences can be written as, Without loss of generality, we can set j = n + 1, where i = 1, 2, 3, • • •, n.The angle difference δ can be formulated as follows, where, δ = [y 1,n+1 , y 2,n+1 , . . ., y n,n+1 ] represents relative phase angle, and is the linear mapping matrix.Hence, it can be observed that there exists a linear mapping between the relative angle δ and the state variable z 1 from mathematical perspective.Combining ( 13)-(I) and (59), the relative motion equation of oscillators in Kuramoto network can be written as follows, where

C. A PHASE-DOMAIN STABILITY INDEX
The following celebrated result in the theory of stability region is relevant.
Theorem 1 ( [32]): (Characterization of stability boundary): Consider an autonomous dynamic system ẋ = f (x), let x i , i = 1, 2, 3, • • • be the equilibrium points on the stability boundary ∂A(x s ) of the SEP x s .Under some mild assumptions, there is, where W s (x i ) is the stable manifolds of x i , where φ(t, x) is solution of ẋ = f (x).□ Consider a single oscillator model as follows, The stability region of (63) can be shown in FIGURE 2. Theorem 1 above suggests that the boundary of the single oscillator model is composed of stable manifolds of W s (δ 1 ) and W s (δ 2 ).Hence, the phase-domain distance between the maximum of phase angle differences and the controlling UEP can be viewed as a numerical index to assess the stability margin of the second order Kuramoto oscillator network.
Therefore, a stability margin index can be defined as follows, where

VI. CASE STUDY
In this section, several power system examples are employed as different second order Kuramoto oscillator network cases to show the effectiveness of QHBM and analysis.

A. SINGL MACHINE INFINITE BUSE POWER SYSTEM
Consider a single-machine infinite bus (SMIB) power system [34] with a diagram of topology as shown in FIGURE 3.   The dynamics equation corresponding to SMIB system is shown as follows, where the parameters in (65) are summarized in Table 1.
A three-phase short-circuit fault occurred at the generator terminal bus with fault durations time t c equal to 0.06 s, and the response of SMIB power system with different damping coefficients are shown in FIGURE 4.
It can be observed that the harmonic solution fits the actual response of the SMIB power system well when the damping coefficient is set as zero, but this is not true when damping is not equal to 0. This demonstrates that the QHBM can only be applied to a zero-damping system, whose response is periodic.Moreover, the harmonic solutions with different selections of H are shown in FIGURE 5.
We can see that the performance of QHBM varies with different selection of H and t c .Roughly speaking, depending on the harmonic composition of the actual solution and nonlinear degree of the SMIB system, the higher of order of harmonic, the more precise of harmonic solution.FIGURE 5 (a) shows that the performance of harmonic solution with H equal to 2 is much better than harmonic solution with H equal to 1 because of the consideration of a higher harmonic component in the Furthermore, taking H equal to 2 as an example, the analytical solution expression obtained by QHBM with t c equal to 0.06 s is shown as follows, y = 0.0395 + 0.0347 cos (7. where t ≥ 0.56.It can also be seen that the coefficient of second order sub-harmonics is quite smaller than the fundamental harmonics in (66), which means that the composition of the actual solution of SMIB system with t c equal to 0.06 s includes virtually fundamental harmonics only.Moreover, the trend of the stability margin index with different fault duration times are shown in FIGURE 6.From the FIGURE 6, we can infer that the critical fault time duration of SMIB system is 0.14 s around because of the nearness of stability margin index to zero as 0.134 s.In order to verify this inference, the response of SMIB power system with different fault time duration are shown in FIGURE 7.
Apparently, the time domain simulation indicates that the critical fault time duration is equal to 0.14 s from FIGURE 7 which coincides with the trend in FIGURE 6, demonstrating the effectiveness and accuracy of the phase domain index proposed in Section V-C.Furthermore, the phase portrait of SMIB system can be depicted through time domain simulation and harmonic solution as shown in FIGURE 8.
It can be seen that the harmonic solution fit the numerical solution very well even with fault time duration equal to 0.12 s from FIGURE 8.Moreover, the harmonic solution with t c equal to 0.12 s is shown as follows, y = 0.2603 − 0.0044 cos (6.1883 × (t − 0.62)) + 0.9371 sin (6.1883 × (t − 0.62)) + 0.0926 cos (2 × 6.1883 × (t − 0.62)) − 0.0013 sin (2 × 6.1883 where t ≥ 0.62.Comparing (67) and (66), we can see that the constant harmonic a and fundamental harmonic amplitude increase as the nonlinear degree increase, indicating that the fundamental harmonic still is the main part of the approximate solution.This reflects, as previously asserted in (31), the low-pass filtering characteristics of second order Kuramoto model.

B. IEEE 3-MACHINE-9 BUS POWER SYSTEM
The diagram of the IEEE 3-machine 9-bus power system is shown in FIGURE 9, and the system parameters are obtained from reference [35].
The dynamic equation of the IEEE 3-machine 9-bus power system is a set of second order Kuramoto equations by setting N = 2 in (1), and the governed equation is shown as follows, We set a three-phase short-circuit fault that occurred at bus 7 as a disturbance to (68), and the fault durations time is equal to 0.08 s.The angle responses of this Kuramoto oscillators network with zero-damping are shown in FIGURE 10.
Applying coordinate transformation y = U z, the response in Z coordinate system of this power system can be obtained from FIGURE 10 (b).It can be seen that the phase angle responses in Z space can be divided into three parts, i.e., z 1 , z 2  and z 3 , where the main subsystem variable z 1 and z 2 oscillates against each other, while the secondary subsystem variable z 3 shows the uniform motion characteristic, which validates the feasibility in Section V-A, explaining the morphological characteristics of oscillator's stable motions.
As can be seen from FIGURE 11, the projections of z 1 and z 2 in different directions can generate different oscillatory behaviors.Denoting , the trajectories of the variables Vz 1 = [y v2 , y v3 ] are shown in FIGURE 12.We see that the relative angle responses coincide with the new variables Vz 1 from FIGURE 12, demonstrating that (60) is capable of describing the relative motion of oscillator in Kuramoto network.The response of this power system with different damping coefficients are shown in FIGURE 13.
It can be concluded that QHBM can only be applied to dynamic system with zero-damping system from FIGURE 13 as illustrated in the SMIB system case.Besides, the comparison between numerical solution and harmonic solution with different selection of H is demonstrated in FIGURE 14.
It can be observed that the harmonic solution with H equal to 2 and 3 fits better than H equal to 1 both in the scenarios with t c = 0.08 s and 0.14 s, respectively, from FIGURE 14, showing that the performance of solution with higher order of harmonic is better than lower order of harmonic.
FIGURE 15 demonstrates that the amplitude of constant harmonic a, fundamental harmonic and second order subharmonic is increased as fault time duration time increases, because of the harmonic distortion.This is consistent with the phenomenon in SMIB power system, further indicating  low-pass filtering characteristic of second order Kuramoto model asserted in (31).
Compared with the smaller amplitude component z 2 , QHBM fit well in the larger amplitude component z 1 as fault time duration increases from FIGURE 16.Furthermore, the phase portrait of component z 1 is shown in FIGURE 17.It can be seen that the phase portrait of z 1 is quite similar to the SMIB system portrait in FIGURE 17.We can infer that the critical fault time duration is between 0.22 s and 0.24 s.Hence, the distance between z 1 and the corresponding UEP can be viewed as a stability margin index according to (64).The variation of I (z(t F )) versus fault time duration is shown in FIGURE 18.
From FIGURE 18, one can infer that the critical fault duration time is greater than 0.22 s, this observation is consistent with FIGURE 17, verifying the validity of (64).Moreover, the assertion can also be verified by the following time domain simulation in FIGURE 19.The result in FIGURE 19 demonstrates that this power system keeps stable with t c equal to 0.20 s while it loses stability with t c equal to 0.24 s, indicating actual critical fault time duration time is between 0.22 s and 0.24 s, which is consistent with the assertion revealed by FIGURE 17 and 18.

C. IEEE 10-MACHINE 39-BUS POWER SYSTEM
The diagram of the IEEE 10-machine 39-bus system (a reduced-order system of the actual New England power  system) is shown in FIGURE 20, and the system parameters are obtained from reference [36].
A three-phase short-circuit fault occurred, and the fault durations time is equal to 0.02 s.The responses of the rotor angle of the power system with non-damping are shown in FIGURE 21.
It can be seen that the rotor angle is moving forward after the clearance of fault in δ coordinate system from  FIGURE 21 (a), and those motions can be transformed into ten components in Z coordinate system, as uncovered by (13).The comparison between the numerical solution and harmonic solution with different damping conditions are shown in FIGURE 22.
FIGURE 22 shows that the harmonic solution in the zero-damping case fits better than damping case, which is consistent with the observation aforementioned.Furthermore, the comparison between numerical solution and harmonic solution with different order of harmonic are shown as follows.
FIGURE 23 demonstrates the performance of QHBM with different harmonic order and fault duration time.We can see that the performance of solution with H equal to 2 is better than the harmonic solution than H = 1 in general both in case with fault duration time equal to 0.06 s and 0.08, respectively.Furthermore, the variation of I (z(t F )) versus fault time duration is shown in FIGURE 24.
FIGURE 24 demonstrates that the critical fault duration time is greater than 0.12 s, this inference can be examined by the time domain simulation as shown in FIGURE 25.
The result in FIGURE 25 manifests that the power system is stable with t c equal to 0.12 s while it loses stability with t c equal to 0.14 s, demonstrating that actual critical fault time duration time is between 0.12 s and 0.14 s, which is consistent with the observation revealed by FIGURE 24.This result verifies the validity of the analysis in Section V.

VII. CONCLUSION
In this paper, a new polynomial model for the second order Kuramoto oscillator network is established.It is shown that the model can be partly decoupled by applying an eigen transformation.A new QHBM is proposed and derived in detail, and the harmonic solution by QHBM is efficient to estimate the actually response of the second order Kuramoto oscillator model with zero-damping.Based on the harmonic solution, the solution of Kuramoto oscillator model can be decomposed into two parts, i.e., multi-harmonic oscillation motion component and acceleration motion component, respectively, explaining the morphological characteristics of stable motion of oscillators in mathematic.The relative motion between different oscillators can be depicted by the proposed nonlinear state output equation system and large disturbance stability margin can be measure by phase domain stability index in this paper.Several power system cases used as different second order Kuramoto oscillator networks demonstrate the effectiveness and accuracy of the proposed approach and analysis.

FIGURE 1 .
FIGURE 1. Illustration of the Kuramoto oscillator network.

FIGURE 2 .
FIGURE 2. Stability boundary of single oscillator model.

FIGURE 3 .
FIGURE 3. Diagram of the SMIB power system.

FIGURE 4 .
FIGURE 4. Simulation of result with different damping coefficient and same fault duration time, i.e., t c = 0.06 s.

FIGURE 5 .
FIGURE 5. Simulation of result with different selection of order of harmonic and fault duration time.

FIGURE 6 .FIGURE 7 .
FIGURE 6. Trend of stability margin index with different t c .

FIGURE 10 .
FIGURE 10.Responses of phase angle in different coordinate system with tc equal to 0.08 s.

FIGURE 11 .
FIGURE 11.Response of nonlinear state variables z 1 and z 2 .

FIGURE 12 .
FIGURE 12. Diagram of relative angles.The solid curve is the subtraction of different phase angle, while the dashed curve is generated by (59).

FIGURE 14 .
FIGURE 14. Simulation of result with different H and t c .

FIGURE 16 .
FIGURE 16.Phase angle responses with different tc.The solid line is numerical solution and the dash curve is harmonic solution with = 2.

FIGURE 17 .
FIGURE 17. Phase portrait of IEEE 3-machine 9-bus system.The solid curve is obtained by QHBM and the dashed curve is numerical solution.

FIGURE 18 .
FIGURE 18. Stability margin index with different fault durations time.

FIGURE 19 .
FIGURE 19.Simulation of result for power system with different t c .

FIGURE 21 .
FIGURE 21.Responses of rotor angle in different coordinate system.

FIGURE 22 .
FIGURE 22. Simulation of result with different damping and same fault duration time and order of harmonic., i.e., t c = 0.06 s, H= 1.

FIGURE 23 .
FIGURE 23.Comparison of numerical solution and harmonic solution.

FIGURE 24 .
FIGURE 24.Stability margin index with different fault durations time.

FIGURE 25 .
FIGURE 25.Time domain simulation for power system with different t c .

TABLE 1 .
The parameters of the SMIB power system.