Low-frequency oscillations in coupled phase oscillators with inertia

This work considers a second-order Kuramoto oscillator network periodically driven at one node to model low-frequency forced oscillations in power grids. The phase fluctuation magnitude at each node and the disturbance propagation in the network are numerically analyzed. The coupling strengths in this work are sufficiently large to ensure the stability of equilibria in the unforced system. It is found that the phase fluctuation is primarily determined by the network structural properties and forcing parameters, not the parameters specific to individual nodes such as power and damping. A new “resonance” phenomenon is observed in which the phase fluctuation magnitudes peak at certain critical coupling strength in the forced system. In the cases of long chain and ring-shaped networks, the Kuramoto model yields an important but somehow counter-intuitive result that the fluctuation magnitude distribution does not necessarily follow a simple attenuating trend along the propagation path and the fluctuation at nodes far from the disturbance source could be stronger than that at the source. These findings are relevant to low-frequency forced oscillations in power grids and will help advance the understanding of their dynamics and mechanisms and improve the detection and mitigation techniques.

. Phase fluctuation in a two-node network consisting of a generator (P 1 =+ 1) and a load (P 2 = −1) driven by an external periodic forcing A sin (2πft) where f = 0.2Hz. The forcing is added at t = 30 s and removed at t = 120 s. (a) Results when the forcing is applied to the generator node 1. Here k = 2, α = 0.5, A = 0.2. (b) Results when the forcing is applied to the load node 2 under the same conditions of (a). (c) The phase difference between the two nodes under the same conditions of (a). (d) Dependence of generator node phase fluctuation peak-to-peak magnitude Δ 1 on the damping coefficient α and forcing amplitude A with k = 2. (e) Dependence of load node phase fluctuation peak-to-peak magnitude Δ 2 on the damping coefficient α and forcing amplitude A with k = 2. (f) Dependence of phase fluctuation peak-to-peak magnitudes Δ 1,2 on the damping coefficient α with k = 2 and 5 and A = 0.2. In this work, we analyze the response of second-order Kuramoto oscillator network motifs to sinusoidal forcing, the general model of which was first given in 5 . This problem is of specific relevance to power grid reliability and stability. On the one hand, modeling power grid as Kuramoto oscillator networks has driven a new interdisciplinary research thrust in the last decade 29,[35][36][37][38][39] (some studies were based on first-order Kuromoto model 40,41 ). Following the initial model formulation 36 , a study demonstrated that decentralized power sources might be conducive to self-organized synchrony, which would enhance gird robustness 37 . New synchronization conditions were derived in terms of network topology and parameters to improve the practical applicability of the theoretical results to smart grid dynamics and control 38,39 . Despite the simplifying assumptions and limitations 37 , this approach roots in the conceptual framework based on the electromechanical model of rotating electric machines and thus able to generate illustrative and insightful results justifiable for real-world applications 42 .
On the other hand, in electric power engineering, the dynamics of the grid under substantial oscillatory disturbances (desynchronized states that may have catastrophic consequences such as blackouts) has been investigated from two perspectives. The first considers the propagation of disturbances in the grid in the form of low frequency (order of 0.1-1 Hz) electromechanical waves 43,44 . A continuum model was constructed to describe the traveling wave along transmission lines 43 , based on which a control method was proposed to damp the dynamics 45 . With the grid getting more and more extensive and complex, it is imperative to ensure that effective protection schemes are in place against damaging electromechanical waves 46 and advanced techniques are used to monitor them 47 . The second perspective scrutinizes low-frequency (sometimes also called inter-area) oscillations observed at specific nodes [48][49][50] . Classical analysis based on equivalent oscillating circuit model 48 explained how "natural" interaction among grid devices results in free oscillations. In contrast, forced oscillations, i.e., the grid (typically in the regime where free oscillations are inhibited) response to persistent external forcing, have become a growing concern in recent decades 50 , stimulating the development of detection algorithms 51 and mitigation methods 52 .
The Kuramoto model is a promising alternative to the continuum model and the circuit model that not only captures disturbance propagation (like the continuum model) and forced oscillation (like the circuit model) in the grid but also sheds light on the effects of network topology (which continuum model is unable to show) and coupling strength (which circuit model does not clarify) on system dynamics. This work performs a numerical study of the desynchronization dynamics of some representative periodically-driven Kuramoto oscillator networks. In the context of this paper, desynchronization refers to a phase unlocked regime (the phase differences between nodes are no longer constant). As shown below, there is an "ac steady state" in which the phase of each node fluctuates at the same frequency as the external forcing around a state that deviates www.nature.com/scientificreports www.nature.com/scientificreports/ from the synchrony. Special attention is paid to the magnitude of this type of phase fluctuation and the effects of network and forcing parameters on the magnitude. In power grids, the knowledge of phase fluctuation magnitude at each node will be essential to the specification of its hosting capacity as well as necessary control measures.

Model
The model in this work follows that of 36,37 , with the addition of sinusoidal forcing at a selected node. Consider a coupled network of N synchronous generators (production) and motors (load). Each node i(i = 1, …, N) can be described as a phase oscillator with its electromechanical phase θ φ = Ω + t i i , where Ω is the constant grid reference frequency and φ i is the relative phase (or simply called phase here). For each oscillator, the generated (consumed) power should balance the power transmitted to (received from) the grid plus (minus) that for local acceleration and dissipation. Assuming , the governing equation of the system at node i is: where P i and α i are the equivalent power and damping coefficient of the i-th oscillator, [K ij ] N×N is the matrix of coupling strength with each element K ij being the product of the connectivity (1 if connected; 0 otherwise) and the coupling strength between nodes i and j, A and ω = 2πf are the magnitude and angular frequency of the external forcing applied at node l in time interval [t 1 , t 2 ], u(t) is the unit step function, and δ il is the Kronecker symbol. Equation (1) is solved using the 4 th order Runge-Kutta method.
In the un-disturbed case (A = 0), to reach synchrony, it is necessary that ∑ = = P 0 i N i 1 (energy balance). Further, in this work, the coupling strengths are so chosen that ≥ | | ↔ K P min j i ij i (j↔i means the two nodes are connected). This generally guarantees the existence of attractive fixed points (stable synchrony) of the system. It is also true in actual power grid since the capacity of transmission lines are by design higher than local generation or load. On the other hand, for each set of parameters (P i , K ij ) used in the simulations, we numerically solve the un-disturbed steady state equation of Eq. (1) and confirm that the system has only one attractive fixed point. The regime with multiple fixed points is avoided in this study to make the simulation results comparable and without ambiguity. www.nature.com/scientificreports www.nature.com/scientificreports/ Under these conditions, for each numerical case, the un-disturbed oscillators' initial phases are chosen randomly from [0, 2π). Once all φ i 's reach equilibria, we record their equilibrium values as the initial conditions for solving Eq. (1). It is expected that before the application of external forcing, the system will remain in the initially synchronized state, and after the removal of forcing, the system will return to the initial state. Under the periodic forcing, the peak-to-peak magnitudes of all φ i 's fluctuations are measured to see how seriously an oscillator is affected. To describe the disturbance propagation, we record the starting time of each oscillator's forced oscillation practically defined as the time when φ i first crosses a pre-set value between the initial synchronized state and the time average of the ac steady state. By comparing the starting times of the oscillators, one can analyze the speed and path of propagation.

Results
two-node networks. We start with the two-node network due to its simplicity. When P 1 = −P 2 = P 0 , α 1 = α 2 = α, 1 2 , and the +/− means the forcing is applied at node 1/2. This is the equation of a damped nonlinear pendulum driven by an external force that has both dc (2P 0 ) and ac (A sin (ωt)) components. Without the ac component, it can be shown that there exists an attractive fixed point at φ = − P k arcsin( / ) . This is the regime under study in this work, i.e., there is a steady state with constant φ 1−2 . Usually this means both φ 1 and φ 2 are constants. Although, as shown in 37 , for smaller k's and certain initial conditions, φ 1 and φ 2 may show identical oscillations which will be canceled in φ 1−2 , we do not consider this rare case and choose a higher k or an initial condition that leads to constant φ 1 and φ 2 . Now with the ac component included, φ 1−2 will fluctuate around its steady-state value. We numerically solve this problem, and some representative results are shown in Fig. 1.
Comparing Fig. 1(a,b), one can see that due to symmetry, the two cases (the forcing at node 1 and node 2) have the same pattern of phase fluctuation. The phase fluctuation has the same frequency as the external forcing. Figure 1(c) shows that φ 1−2 also reaches an ac steady state with the same fluctuation frequency. This is an interesting result of relevance to the inter-area low-frequency oscillation in the power grid. The two nodes in Fig. 1 represent two grid zones where φ 1−2 indicates the direction and amount of power flowing from node 1 to node 2. Under the conditions of Fig. 1, < 1 P k 0 , which practically guarantees that φ 1−2 does not have sustained free oscillations; however, with the periodic forcing, φ 1−2 fluctuates around its normal operating point, which generally www.nature.com/scientificreports www.nature.com/scientificreports/ undermines the reliability and stability of a power grid. Note that the natural frequency of motion is k 2 (no damping or forcing). The resonant frequency of the damped driven system is lower than this value. The closer the external forcing frequency f is to the resonant frequency, the higher the magnitude of the fluctuation. In the case of Fig. 1, = k 2 2 Hz; so the most interesting results are to be observed in the low frequency range, i.e., f ≲ 2 Hz. Higher frequency fluctuations are less important in terms of magnitude (see Fig. S1).
Further, in Fig. 1(d,e), we present the phase fluctuation peak-to-peak magnitudes (Δ 1,2 ) at the two nodes with different damping coefficients (α) and forcing amplitudes (A). Keeping other conditions the same as Fig. 1(a), Δ 1,2 increase about linearly with A, while increasing α results in higher Δ 1 and lower Δ 2 . This counter-intuitive effect of damping is also plotted in Fig. 1(f), where, for comparison, the case of k = 5 is shown. Since this phenomenon can only be observed when f ≤ 0.2 Hz, its implication is that damping might not be effective for the reduction of very low-frequency phase fluctuations with relatively weak links. The detailed analysis will be reserved for continuing studies.
On the other hand, when f > 0.2 Hz, it is found that not only the dependency of Δ 1,2 on α is similar to the case of k = 5 in Fig. 1(f), but also there is a "optimal" coupling strength k at which Δ 1,2 peak. As shown in Fig. 2(a), when f = 0.5 Hz, both Δ 1 and Δ 2 reaches maxima at k ≈ 6. For smaller k's, Δ 1 > Δ 2 ; for larger k's, Δ 1 < Δ 2 . Similar results are in Fig. 2(b) for f = 1 Hz, with the "optimal" coupling strength k ≈ 21. Note that in both cases, the forcing amplitude A = 0.2, and the peak Δ 1,2 can be several times higher than A. Unlike the majority of literature on synchronization that requires some critical (minimum) coupling strength 38 , it is a new finding that, in a system as simple as the two-node network, there is another critical coupling strength under which the forced phase fluctuations resemble resonance. The difference is that here the variable is the coupling strength instead of the driving forcing frequency. In Fig. 2(c,d), we show that the "optimal" coupling strength k as a function of the forcing frequency f is insensitive to other factors such as α and A, implying that this phenomenon is due to the nonlinear interaction between the two nodes. Interestingly, Fig. 2 bears some similarity to the frequency response of forced Duffing oscillator 53 (more details in SI, Note 1). three-and four-node networks. We now move on to three-node (and four-node) systems. Figure 3 presents the phase fluctuation in a linear three-node network with forcing (f = 1 Hz) applied to node 1. Comparing the four cases, one concludes that 1) the nodal power distribution has an insignificant (<2%) effect on Δ's ((a) versus (c), or (b) versus (d)), and 2) the variation in coupling strengths can lead to distinct patterns of Δ ((a) versus (b), or (c) versus (d)). To confirm the first point, we perform more simulations for confirmation (see Fig. S2). The second point is consistent with what we have seen in the two-node network. It is worth mentioning that these conclusions also hold for circular three-node network and four-node network with a branch (see Figs. S3 and S4). www.nature.com/scientificreports www.nature.com/scientificreports/ Physically, the nodal power serves as a dc driving force, so it is understandable that it has limited effect on the oscillatory motion. It is the interaction between nodes (characterized by the coupling strength) that "spread" the external forcing effects over the network.
The "resonance" phenomena are also found in the linear three-node network. In Fig. 4, keeping k 12 = 30 and varying k 23 , we plot the phase fluctuation peak-to-peak magnitudes Δ 1,2,3 under different forcing frequencies. Similar to the two-node case, when f > 0.2 Hz, the "optimal" coupling strength starts to appear at which at least one node's Δ peaks. The "optimal" coupling strength also increases with increasing f. For 0.4 ≤ f ≤ 0.7 Hz, we have Δ 3 > Δ 1 > Δ 2 at the "optimal" coupling strength; increasing f to 0.9 Hz, Δ 1 becomes the highest and Δ 2 does not show an obvious peak near the "resonance" point. Note that the forcing is added at node 1. In power grids, to detect forced oscillations, it would be desired to monitor a node with maximum Δ. From Fig. 4, one can see that this node is one of the end nodes.

Longer chains and rings.
To study the forced oscillation in more complex systems, we consider the propagation of phase fluctuation as well as the distribution of phase fluctuation magnitudes. Figure 5(a) illustrates the configuration of a linear 20-node network with global coupling strength driven by external periodic forcing at node 1. Figure 5(b) shows the propagation of the phase oscillations along the chain. The higher the coupling strength k, the faster the phase oscillation propagates to the last node in the chain. After k is over 100, this trend comes to saturation. Due to the limitation of the Kuramoto model, this cannot be directly related to the speed of propagation since there is no length for each link. Nevertheless, one can still see that stronger links facilitates the propagation of low-frequency phase oscillations. Additional results in Fig. 5(c,d) show that the forcing frequency and the damping coefficient have negligible impacts on the propagation. In Fig. 5(c), the apparent delay in the f = 0.2 Hz case results from the method used to measure t osc . Figure 6 presents the distributions of Δ along the 20-node chain under various conditions. If Fig. 5 views the phase fluctuation as "traveling wave", Fig. 6 displays something like "standing wave". In Fig. 6(a), we find that the higher the forcing frequency, the more "bumps" there will be in the distribution of Δ. There is a decaying trend from node 1 to node 20, but the decay is not monotonic. Also in general, higher frequency corresponds to lower Δ. Comparing Fig. 6(b) with 6(a), one sees that the higher the coupling strength, the fewer "bumps" there will be in the distribution of Δ. Figure 6(c,d) indicate that the effect of the damping coefficient on the distribution of Δ is more significant under higher coupling strengths. Similarly, we obtain the distributions of ∆ along an N-node chain where 3 ≤ N ≤ 20 (Fig. S5). An interesting observation is that the fluctuation magnitude at node N is always a local maximum and in the cases of N ≥ 7, the value of ∆ at the 7th node counting backwards starting from node N is always a local minimum. In the cases of N ≥ 14, the value of ∆ at the 8th node counting backwards starting from this local minimum becomes a local maximum again. This pattern continues with longer chains and under different parameter settings (e.g., when f = 0.5 Hz, the numbers will be 3 and 4 instead of 7 and 8). The distribution of fluctuation magnitudes can be qualitatively inferred from the node that is the farthest from the forcing node, while the common sense may suggest the other direction. In Fig. S6, we plot the values of ∆ at node N in all the cases above, which display an oscillatory damping trend as the chain elongates.
When a link is added between node 1 and node 20 in Fig. 5(a), the above network becomes a ring, as illustrated in Fig. 7(a). The resulting distributions of Δ are in Fig. 7(b-d). As expected, since there are now two pathways of propagation, the distributions of Δ show a decaying trend along both pathways from node 1 to node 11. The three sets of results point out that the higher the forcing frequency, the more difficult it can "penetrate" the network (most links with coupling strength k = 10). In the case of f = 1 Hz, the phase fluctuations are localized near nodes 1, 2, and 20. This phenomenon could inspire the design of power grids to damp the forced oscillations and improve system stability.
To make it comparable with Fig. 6, we simulate the forced oscillation in the ring now with global coupling strength k. Figure 8(a) shows the propagation along the ring. Figure 8(b) presents the propagation times from node 1 to node 11 (the farthest one from external forcing) under various k's and f's. Again, the forcing frequency does not affect how fast the phase fluctuations propagate. The higher the coupling strength, the faster the propagation. With f = 0.2 Hz and k in a range from 30 to 100, the distributions of Δ in the system are plotted in Fig. 8(c). Different from the chain in Fig. 6, the ring's Δ distributions have similar profiles under various k's (increasing k tends to flatten the profile), and node 11 always has the highest Δ. In addition, as shown in Fig. 8(d), the ring's Δ distributions are symmetric along the two pathways. In Fig. S7, we show the results of a ring with 19 nodes and other conditions the same as Fig. 8(d). Now there are two farthest nodes (10 and 11) and the distribution features are very similar. Therefore the Δ distribution is not determined by even or odd number of nodes in the ring.

conclusions
This work considers a second-order Kuramoto model periodically driven at one node as the model of forced oscillations in power grids. Note that the un-forced system stays within the regime of stable synchrony and the external forcing does not introduce new dynamic regimes or cause regime changes. Motivated by the power grid application, we study the phase fluctuation magnitude at each node as well as the disturbance propagation in the network. There are three major findings. Firstly, given the network topology, coupling strength, the