Dynamics and Collapse in a Power System Model with Voltage Variation: The Damping Effect

Complex nonlinear phenomena are investigated in a basic power system model of the single-machine-infinite-bus (SMIB) with a synchronous generator modeled by a classical third-order differential equation including both angle dynamics and voltage dynamics, the so-called flux decay equation. In contrast, for the second-order differential equation considering the angle dynamics only, it is the classical swing equation. Similarities and differences of the dynamics generated between the third-order model and the second-order one are studied. We mainly find that, for positive damping, these two models show quite similar behavior, namely, stable fixed point, stable limit cycle, and their coexistence for different parameters. However, for negative damping, the second-order system can only collapse, whereas for the third-order model, more complicated behavior may happen, such as stable fixed point, limit cycle, quasi-periodicity, and chaos. Interesting partial collapse phenomena for angle instability only and not for voltage instability are also found here, including collapse from quasi-periodicity and from chaos etc. These findings not only provide a basic physical picture for power system dynamics in the third-order model incorporating voltage dynamics, but also enable us a deeper understanding of the complex dynamical behavior and even leading to a design of oscillation damping in electric power systems.


Introduction
Electric power system is a very complicated dynamical system [1], due to its intrinsic property of high nonlinearity. Its dynamical behavior is of essential importance for a secure and stable operation of power systems, such as rotor-angle stability, voltage instability (or voltage collapse), and frequency stability. Therefore, it is very natural that the study of dynamical behavior in power systems, namely, the so-called power system dynamics [2], has been a very important topic in power electrical engineering. Meanwhile, it is no doubt that nonlinear dynamics, which has been developed to study general dynamical systems with nonlinearity, has become an inevitable and important mathematical tool in power system dynamics.
The remainder of this paper is organized as follows: In the section of mathematical model, the third-order SMIB model of a power system is introduced. The results for its dynamic behaviors by considering both positive and negative damping coefficients are presented in detail in the section of rich dynamics and system collapse. In the section of theoretic analysis, the local stability of the equilibrium state is studied. Finally, the conclusions and discussions are given. Fig 1(A) schematically shows the SMIB power system, where the electric power is generated from the synchronous generator (G) and transferred to an infinitely large bus through the transformer (X T ) and two parallel transmission lines (X L ). Here the infinite bus means that the voltage magnitude and phase are constant and unchanged with time. Due to its simplicity, this representation has been extensively used in the study of power system dynamics. The models of power systems can be presented at several diverse levels of complexity, relying on the studied problems. (See our discussions in the introduction part.) Here with some typical assumptions, we will employ and derive the classical third-order flux decay model describing also the voltage dynamics [32][33][34].

Mathematical Model
According to the basic principle of a synchronous generator, its rotor motion equation can be represented by where δ is the rotor angle of the synchronous generator, ω is the rotor speed with respect to synchronous reference, M is the moment of inertia, D is the damping factor of the synchronous generator, P m indicates the mechanical input power to the synchronous generator which is usually assumed to be constant in short-term power system stability problems, and P e represents the active electrical power delivered by the synchronous generator.
According to the principle of the conservation flux of the closure windings, the rotor windings dynamic equation of the synchronous generator can be modeled by ignoring the impact of damping windings. Fig 1(B) schematically shows the fixed abc coordinate system defined by the aa', bb' and cc' axis of the stator and the dq coordinate system rotating with the rotor in a synchronous generator (d-axis centered in the rotor field's magnetic north pole N and the qaxis perpendicular to the S-N axis). Taking into account the electromagnetical interactions between the involved field, damping and stator windings given in the abc-stator coordinate system (a, b and c denoting the three stator phases), finally we have which is formulated in the dq-rotor coordinate by the abc ! dq0 Park's transformation. Here E is the quadrature-axis transient voltage of the synchronous generator, E q is the quadratureaxis voltage of the synchronous generator, T d0 represents the direct-axis open-circuit transient time constant of the synchronous generator, and E f denotes the field excitation voltage. Based on the corresponding vector relation [35], we obtain the following algebraic equations: Fig 1. (a) and (b) Schematic shows for the single-machine-infinite-bus power system and the coordinate systems in a synchronous generator, respectively. In (a), here G represents the synchronous generator, x T indicates the reactance of the transformer, x L is the reactance of the transmission line, V s denotes the infinite bus voltage, and V t denotes the terminal voltage of the synchronous generator. In (b), the fixed abc coordinate system is defined by the aa', bb' and cc' axis of the stator with the dq coordinate system rotating with the rotor (d-axis centered in the rotor field's magnetic north pole N and the q-axis perpendicular to the S-N axis).
The active power P e of the synchronous generator G is expressed as follows, where x q∑ = x q +x T +x L and x ' d∑ = x ' d +x T +x L , with x q the quadrature-axis synchronous reactance, x ' d the direct-axis transient reactance, x T the reactance of the transformer, and x L the reactance of the transmission line. The symbols I d and I q indicate the direct-axis and quadrature-axis currents of the synchronous generator, respectively, and V s and V t represent the infinite bus voltage and the terminal voltage of the synchronous generator, respectively, as shown in Fig 1. Substituting the electrical Eqs 4-7 into the mechanical and electrical dynamics equations of system Eqs 1-3, we obtain the complete mathematical model of the third-order SMIB power system as follows: In this paper, for convenience, we set B = 1=x 0 Clearly now we have three state variables δ, ω and E in the model, and they are expressed together as the state vector X = [δ, ω, E] T . The other letters including M, B, V s , X, E f , α, P m and γ are all system parameters, which do not change with time. In fact, here we can factor out the inertia constant (M) in the left of Eq 12 completely. However, in order to investigate its effect, we keep it and make M as a system parameter. Recently the low-inertia stability issues including virtual inertia-emulating devices and virtual inertia placement schemes are of key importance in renewables integration to grids and power system dynamics, and thus it would be very interesting to study the effect of the inertia constant for different M's [33]. Compared with the classical swing equation, which includes only Eqs 11 and 12 and treat the voltage magnitude E as constant (E = 1.0 is chosen in the paper, without losing generality), the whole three-order model Eqs 11-13 treats E as a state variable, which can vary with time. In the power system field, this model more precisely characterizes the synchronous generator dynamics with excitation control.
In our simulations, we choose the mechanical input power P m and the damping coefficient γ as our primary operation parameters, but fix all other parameters: B = 1.0, V s = 1.0, E f = 1.0, α = 2.0 and X = 1.0 [25], which are typical in per-unit in power systems. The classical fixed-step fourth-order Runge-Kutta method for solving differential equations has been used with the time step being 0.01, which is small compared to the period. Therefore, below we will mainly investigate the dynamical behaviors with the variation of P m and γ, and especially focus on the difference with different values and signs of the damping coefficient γ, which can be positive or negative. We are also interested in the difference with the classical second-order swing model which ignores the voltage dynamics.

Rich Dynamics and System Collapse
The case for positive damping In this section, we will mainly rely on numerical simulations by using several standard techniques in nonlinear dynamics, such as time series analysis, phase diagram, Poincare section, Lyapunov exponents, etc. First let us consider the case for the unit one inertia constant, M = 1, and the positive damping, γ>0. The effect of the inertia constant will be studied in the end of this section. The second-order SMIB power system [Eqs 11 and 12] with E = 1 fixed in (12) and the third-order SMIB power system [Eqs 11-13] will be treated and compared. Fig 2(A) and 2(B) are the phase diagrams for different dynamic behaviors in the P mγ parameter space for these two systems, respectively. Simply we use I to indicate the area of stable fixed point, III for stable limit cycle, and II for their coexistence. Based on these observations, although the precise values for different types of dynamics are different, their qualitative behaviors are similar [see the three major areas and compare the different locations of the critical curves in the figures]. In the whole parameter region, no more complicated behavior, such as quasi-periodic, can be found. Therefore, we may conclude that if the damping is positive, considering the voltage dynamics may have not significant influence on the corresponding dynamics. Note that these three different types of dynamics in the second-order system have been well recognized and studied in different fields, such as in power system dynamics, in mechanics for forced pendulum, and in physics for Josephson junctions, etc [13].
To show the similarity to the second-order system in detail, we choose the parameter set (P m = 0.5 and γ = 0. Therefore, we know that although the inertia constant does not affect the equilibria, which can be verified in the theoretic analysis part, but their local stability properties and the transient behavior are entirely different. On the other hand, it is well-known that robustness to disturbance is more important than a large parameter region of attraction, in term of proper operation of power systems. Thus further parameter sensitive analysis and even dynamical analysis under either noise or intentional disturbance are imperative for our deeper understanding of inertia effects on such a classical third-order flux decay equation. The case for negative damping Now we will study dynamic behaviors for negative damping (M = 1 fixed first without losing generality). As we have stated in the introduction, many researchers have considered the positive damping effect and investigated the related excitation control. Several decades ago, however, Concordia and Carter have analyzed the negative damping of electrical machinery and proposed that the negative damping coefficient is related to the generator angle and the armature resistance [36]. Based on the fact that these features determine the property of the generator, the negative damping will play a significant role in the stability of the synchronous generator and further of the whole power systems. Thus it becomes very crucial to study the negative damping effect in power system dynamics. Below we will see that with the consideration of voltage dynamics, the system dynamics indeed becomes much more complicated and distinctive, compared with the classical second-order model. Fig 5 illustrates the phase diagram in the P mγ parameter space for the third-order system with γ<0 (M = 1), where it can be roughly divided into three areas again. Now we use I to indicate stable fixed point, II to represent system collapse, and III to show very complicated dynamical behaviors with different colors including periodic orbits (orange-yellow), quasiperiodic orbits (red), and chaos (yellow). The area of quasi-periodic motions is located only in the small region of the left corner of area III. As there are periodic windows in the chaotic regions, the areas for periodic orbits and chaotic orbits are scattered and interlaced within area III. We will go to more details in the subsequent parts.
For the sake of brevity, several different parameter values corresponding to the typical dynamic behaviors are selected within different parameter areas. For example, within area I as the mechanical power P m and the damping coefficient γ equal 1.2 and −0.117, respectively, the system clearly shows a stable fixed point (Fig 6). This finding verifies that we can still find a large parameter region outside of the traditional region for γ>0 as shown in Fig 2. This new region for a stable equilibrium point with negative damping is of great significance for our understanding of dynamical behaviors and proper operation in power systems, and it can be viewed as one key difference to the second-order power system.
Next let us concentrate on the dynamics within area III and see how the stable fixed points within area I become unstable. Without losing generality, we fix P m = 0.5 and change the values of γ. The bifurcation diagram with the extreme value of the time series δ(t) for changing γ is shown in Fig 7, which exhibits two pieces of periodic and chaotic regions but with clearly different nature. From right to left, with a decrease of γ, the first region of periodicity happens only within a very tiny parameter region, as indicated by the arrow, after the immediate instability of the fixed point. The first region of chaos happens from the instability of fixed points and subsequently from the instability of periodic orbits. Our analysis shows that the fixed point within area I becomes unstable via a supercritical Hopf bifurcation and this critical loci from area I to area III can be analyzed, as we will show later. In contrast, the second region for chaos happens due to a period-doubling cascade. Therefore, the whole route becomes: fixed point ! periodic motion ! chaos ! periodic motion ! period-doubling cascade ! chaos.

Fig 2. Comparison of phase diagrams for γ>0 in the second-order power system [Eqs 11-12] (a) and the third-order power system [Eqs 11-13] (b).
Here the symbol I indicates that of stable fixed point, III denotes that of stable limit cycle, and II represents the area for their coexistence (namely, for different initial conditions two different attractors can appear and coexist). This comparison for qualitative similarity indicates that the voltage dynamics has not significant influence on the dynamic behaviors if the damping is positive. In 2(a), E = 1.0 is chosen for the second-order system throughout the whole work. In 2(b), the letter A indicates the parameters: P m = 0.5 and γ = 0.2. All other parameters M = 1.0, B = 1.0, V s = 1.0, E f = 1.0, α = 2.0 and X = 1.0 are fixed. Clearly a constant critical parameter P mSNB = 1.299 exists for the third-order system, compared to P mSNB = 1 for the secondorder one, which comes from the similar saddle-node bifurcation and can be analytically obtained.   Fig 2(B) for the third-order system]. They are simulated with the same parameters but two different initial conditions. This coexistence can also be extensively found in the secondorder system under the condition of positive damping. In (d), the angle is treated as mod 2π. treated as mod 2π, corresponding to a periodic rotation of the pendulum. Obviously chaos appears in 8(c) and 8(f). This means that in the three-order power system by including voltage dynamics, chaos becomes very typical if the negative damping coefficient has been properly chosen. This can also be viewed as one key difference with the second-order power system.
With further decreasing γ, the chaotic motion may also become unstable and the system collapses. As an example, Fig 9 shows such an instability for γ = −0.1304. It is clear that the system collapses after a very long transient time with a chaotic behavior. This pattern for the disappearance of a chaotic attractor and the accompanying system collapse usually appears when a strange chaotic attractor collides with another unstable limit set, which is probably an unstable limit cycle or a saddle. In the literature, it is generally called boundary crisis [37,38]. In addition, we find that the system collapses for the phase variables (δ and ω) only, whereas the voltage magnitude (E) still keeps a nearly constant value after the transient. Later we will see this type of partial collapse is very common in this model.
Next in order to exhibit the generality of the route to chaos by period-doubling bifurcation within the area III, we study the dynamics with the increase of the mechanical power P m but let the damping coefficient unchanged (γ = −0.117). The results are shown in Fig 10, with period-1, period-2, and chaotic motions presented in Fig 10(A), 10(B) and 10(C), respectively. Correspondingly, P m = 0.01, 0.024, and 0.035. These panels prove that the period-doubling scenario transition to chaos is very common in the third-order system for negative damping. In addition, to our surprise, with the increase of P m and still within the area III, we find another unusual behavior; it is still chaotic in each piece of attractor but nearly periodic between these pieces. This phenomenon is exemplified in Fig 10(D) for P m = 0.0375. We may call it periodic chaos. Detailed calculation for the equilibrium points and their stabilities shows that the first piece of chaotic attractor rotates around one equilibrium point, which is unstable, and after several rotations it suddenly collides with the other equilibrium point, which is also unstable, and jumps to the next piece of the chaotic attractor. This phenomenon can extend further with time. Thus this type of behavior is similar to the multi-scroll attractors [39,40], such as the well-known Chua's double scroll, where multiple equilibrium points in the scrollbased chaotic attractor by modifications of nonlinear characteristics are necessary.
Relying on the calculation of the Lyapunov exponents and dynamical analysis, still within the area III we can even find another type of dynamical behavior, quasi-periodic motion, which appears only in a very tiny parameter region close to P m = 1.3 and γ = −0.5. To be clearer, we use the symbols QP in Fig 5 to  In addition, we find that this quasi-periodic motion can become unstable. See, e.g., the pattern in Fig 11(C), 11(D) and 11(E) for P m = 1.298968 and γ = −0.50001. Clearly this is another type of collapse, which is different from the usual collapse of chaos. However, we find again that the phase variables (δ and ω) go infinity immediately and the voltage magnitude (E) remains finite after the transient. This observation may not be a surprise; actually it can be proven by the applications of variations of constant formula in Eq 13.
Furthermore, let us move on the general behavior within the area II. Without losing generality, two different parameter sets are chosen, one within the lower II and the other within the upper II and both are far from the critical curves. The explosive behaviors for a partial system Dynamics and Collapse collapse for the phase variables only are expected there. Indeed we find that now the phase (rotor-angle) δ goes to infinity after a very short transient, whereas the variable E keeps finite but fluctuates, as shown in Fig 12. This phenomenon demonstrates that although the third- Phase diagram for γ<0 in the third-order power system. Here area I indicates stable fixed point, II represents explosive solutions for system collapse, and III shows very rich dynamics, such as periodic orbits (orange-yellow), chaos (yellow), and quasiperiodic orbits (red and at the left corner of III, indicated by the capital letters QP). In addition, the theoretic results for the saddle-node bifurcation (the horizontal line) and the supercritical Hopf bifurcation (the curved line) have been superimposed with two dashed lines, showing perfect fits with the numerical results. Clearly a constant critical parameter P mSNB = 1.299 still exists, the same as in Fig 2(B) for γ>0. In a sharp contrast to all these, the second-order power system is quite simple; the whole parameter space is occupied by the explosive solution. Finally, similar to the study in Fig 4 for positive damping, let us consider and study the inertia effects in the third-order power system for negative damping. The phase diagram on the Mγ parameter space is shown in Fig 13(A); P m = 0.5. As one example, the bifurcation diagram  Fig 13(B). Again area I indicates a stable fixed point, II represents explosive solutions for system collapse, and III shows rich dynamics, such as periodic orbits (orange-yellow), and chaos (yellow). From right to left in Fig 13(B), with decreasing γ, the system appears a variety of dynamic behavior, from fixed point ! periodic motion ! chaos ! periodic motion ! period-doubling cascade ! chaos ! collapse, except that now the system collapses at γ = −0.2087. Comparing these diagrams with those in Figs 5 and 7 for M = 1, we find that the behavior is quite similar. Observing Fig 13(A), however, we find that the region for area I for the stable fixed point becomes larger with the increase of M; this behavior is quite different to that for positive damping in Fig 5, where smaller inertia favors the stable fixed point, in term of the size of stable parameter region of fix points. This difference can be regarded as one key distinction induced by inertia.
In the power system dynamics field, engineers are quite familiar with the technology of modal analysis in studying small-signal stability (or called small-disturbance stability) [2]. Therefore, it would be very interesting to view the behavior of the eigenvalues of the fixed point as the damping moves from positive to negative. The result of this change for different system parameters γ from γ = 2.0 to γ = −0.5 with Δγ = −0.1 in the eigenvalue space is shown in Fig 14, where it is clear that a pair of leading eigenvalues (λ 1 and λ 2 ) comes across the imaginary axis from left to right at the bifurcation point γ = −0.0151, whereas the other eigenvalue (λ 3 ) remains real and negative. This value of bifurcation point coincides with that in the phase diagram in Fig 5, and this type of behavior is indicative of a super-critical Hopf bifurcation. In the theoretic analysis part, we will further obtain more details in a more mathematical manner. In addition, we know that the fixed point is always locally stable under the condition of positive damping (γ>0); this point is also in accordance with the observation in Fig 2(B). It is apparent that with this local linearization method of the fixed point, we can only obtain local stability or instability information but we cannot get further information for parameters within the region of complex behavior in

Theoretic Analysis
So far, based on numerical simulations we have obtained exhausted observations and have found several unusual behaviors, such as the stable fixed point, chaos, periodic motion, quasiperiodic motion, and their collapses when a negative damping parameter is considered. Apparently their appearance fundamentally contradicts with our common sense that negative damping can only disrupt the system, such as making the rotor-angle diverge in either oscillating or monotonic manner in power systems. In particular, the newly born region of Fig 9. System collapse due to boundary crisis. The parameters P m = 0.5 and γ = −0.1304 are chosen quite close to but outside of the stable chaotic region (i.e., the region near the critic curve for dividing areas II and III in Fig 5). Clearly here the system exhibits chaos in the transient process, and E still keeps a finite value after the transient. equilibrium state is also of special interest from the power electric engineering point of view, as it is a novel type of phenomenon when the third-order model incorporating voltage dynamics is considered. Hence it becomes valuable to further understand such an unusual behavior mathematically, which will be the main objective of this section. For convenience of theoretical analysis, we take M = 1.0.

Dynamics and Collapse
The saddle-node bifurcation According to Eqs 11-13, the equilibrium points for the system are the solutions of the following set of algebraic equations: 0 Based on the first two equations [Eqs 14 and 15], we immediately know that all the equilibrium points are determined by the mechanical power P m only and they are unrelated to the damping coefficient γ. Further we can obtain how the two variables E and δ and δ are related to the only adjustable parameter P m , respectively, as follows: The functional forms in the above two equations are very complicated, and therefore, we cannot have the explicit expressions for E and δ as a function of P m . Basically they have two positive solutions for E and δ, respectively, and the system will have two equilibrium points. Their relation can be stated as follows: 1. If 0 P m < P mSNB , then the system has exactly two equilibrium points denoted as X s = (δ s ,0,E s ) and X u = (δ u ,0,E u ), where X s is the stable equilibrium point, and X u is the unstable one with δ s <δ u and E s >E u ; 2. If P m = P mSNB , then the system has only one equilibrium point with X SNB = (δ SNB ,0,E SNB ).
The parameter point P m = P mSNB is the so-called static bifurcation point in power system dynamics; 3. If P m >P mSNB , then these two equilibrium points collide and annihilate, indicative of a saddle-node bifurcation. After several algebraic manipulations, we specifically get For the parameters used in the paper: B = 1.0, V s = 1.0, E f = 1.0, α = 2.0, and X = 1.0, we have the predicted value for P mSNB : P mSNB = 1.299. P mSNB is independent of γ, and this value exactly fits with our observations [see the horizontal lines for both γ>0 and γ<0 in the phase diagrams in Figs 2(B) and 5, respectively]. Correspondingly, we have the explicit expressions for several special values of the rotor angle δ and the transient voltage E for the saddle-node bifurcation, such as δ min , δ max , δ SNB , E min , E max , and E SNB , where δ min and δ max (E min and E max ) are minimal and maximal values of δ (E), coming from the solution in Eqs 17 and 18. With the stable branch δ s 2(δ min ,δ max ), E s 2(E SNB , E max ) and the unstable branch δ u 2(δ min ,δ max ), E u 2(E min , E SNB ), we obtain dynamics, such as periodic orbits (orange-yellow), chaos (yellow). In (b), from right to left: the transition scenario is similar, from fixed point ! periodic motion ! chaos ! periodic motion ! period-doubling cascade ! chaos ! collapse at γ = −0.2087. We can find that opposite to the case for positive damping in Fig 4, and know the solution regions for the stable or unstable equilibrium points.
The super-critical Hopf bifurcation The system Jacobin matrix based on the linearized equation can be stated below: where the values of δ 0 and E 0 for equilibrium points are determined by the algebraic equations [Eqs [14][15][16].
Solving the characteristic equation for the above Jacobin matrix, we yield for the characteristic value (eigenvalue) λ.
Based on the observations of the instability of the equilibrium points and the immediate appearance of limit cycles within the negative damping region in Assuming the characteristic equation above has the following form: we obtain which produces the following equities: Further we have which can be simplified to This is the exact condition for a Hopf bifurcation. Obviously these theoretic results help us to understand the dynamical behaviors after the static bifurcation. To prove them quantitatively, we plot these predicted values as the critical curves with two dashed lines in Fig 5: one for the saddle-node bifurcation with the horizontal line and the other for the supercritical Hopf bifurcation with the curved line; their fit with the numerical results in the original phase diagram is perfect. Further calculation can prove that the Hopf bifurcation belongs to the supercritical type. The Hopf bifurcation showing the physical picture for a pair of characteristic eigenvalues comes across the imaginary axis from left to right accompanying with the critical parameter, γ = −0.0151 in Fig 14, can also be well confirmed by this analysis. In addition, all these bifurcation analysis results including the saddlenode bifurcation and the supercritical Hopf bifurcation have been confirmed by the ordinary differential equation bifurcation analysis software MATCONT in MATLAB [41].

Conclusions and Discussions
In conclusion, nonlinear dynamic properties of the third-order SMIB power system with the classical flux decay generator model have been investigated in detail and in particular damping effects with both positive and negative damping coefficients have been considered. We have mainly found on the one hand, when the damping is positive, the third-order system and the second-order one are quite similar, exhibiting only three distinct types of dynamical behavior: stable fixed point, stable limit cycle, and their bistability. Therefore, even in this third-order differential equation, positive damping cannot induce quasi-periodicity or chaos. On the other hand for the parameter region for negative damping, we have uncovered very rich and complicated dynamics, such as stable fixed point, limit cycle, quasi-periodic motion, chaotic motion, and periodic chaos, which is similar to multi-scroll behavior. In addition, we have uncovered that the system can collapse due to a boundary crisis via chaos, and even via quasi-periodic motion. The partial collapse as an interesting type of system collapse for certain variables only is found as a quite general phenomenon here. From power systems point of view, the finding of partial collapse for angle variable only indicates that even with the third-order model studied here, rotor-angle stability can only be observed. From nonlinear dynamics point of view, it is also very interesting. All these complicated dynamics are fundamentally different from the behavior for the second-order power system, where the system collapses immediately after the damping is switched from positive to negative. Moreover, two different types of local bifurcations from the equilibrium state, the saddle-node bifurcation and the supercritical Hopf bifurcation, are theoretically analyzed. These analytical results are of great significance to understand the negative-damping-induced stable equilibrium state and to take a fresh look at the negative damping effect and oscillation phenomena in power systems. They are indeed well confirmed by the numerical results, as shown in the phase diagram (Fig 5). In addition, the inertia effects are studied for either positive (γ>0) or negative damping (γ<0) in Figs 4 and 13, respectively. The results show that all complicated dynamics for M = 1 can be found for any arbitrarily chosen M. However, for positive γ, smaller inertia enlarges the dominant stable fixed point region, and oppositely, for negative γ, larger inertia enlarges the stable fixed point region. Thus, inertia does show different impacts for either positive or negative damping. Further systematical comparison with transient behavior, parameter sensitivity, and dynamical behavior under either noise or intentional disturbance for small or large inertia constants are important.
Finally it is valuable to give further discussions as follows. (1) As the negative damping coefficient is closely related to the synchronous generator angle and the armature resistance, and the negative damping coefficient will affect the stability of a power system, the study in the present work considering both positive and negative damping coefficients is crucial and meaningful. In the third-order model in the paper, we considered the magnetic flux decay and we did not take the controllers into consideration, such as Automatic Voltage Regulator (AVR) and Power System Stabilizer (PSS). The controllers AVR and PSS are just designed to avoid the appearance of negative damping. It has ever been a very hot topic to study these important controllers in power systems [11], e.g., an adaptive power system stabilizer which can cancel the negative damping torque of a synchronous generator was designed [12]. However, it is necessary to investigate these pure negative damping effect in the absence of these classical controllers. (2) Other nonlinear dynamics results here, such as the fixed point, chaos, quasiperiodicity, and different collapses etc, are of great importance to understand the complicated dynamical process in power systems and might help us to design an oscillation damping device. (3) These results considering the voltage dynamics of synchronous generator are also highly relevant to voltage stability and voltage collapse, which are still unsolved hard problems in power systems. (4) When considering only equilibria and their local stability properties, there are quite a few studies, similar to ours, that show the swing equation itself may deliver very different results from higher-order models. See, e.g., several recent papers [7,8]. (5) Last but not least, for more complicated cases, such as more detailed models for synchronous generator considering damping windings, the interaction of multiple synchronous generators [42][43][44][45][46][47][48], different types of stability in power systems including rotor angle stability [49] and voltage stability [50][51], and interplay with delayed power price [52] so on, further theoretic analysis methods have to be developed. It is of interest to study these problems from nonlinear dynamics point of view in the future.