Multistability and variations in basin of attraction in power-grid systems

A power grid refers to an interconnected network of electric power transmission lines, whose node elements are composed of power plants and substations. In general substations are considered as consumers, since substations are the gateways where electricity flows out to the consumers through distribution networks. The other group is composed of power plants, which plays the role of producers of electricity in power grids. We specify the power producers with positive values of net-power input ${P}_{i}\gt 0$ and the consumers as negative values ${P}_{i}\lt 0$ ($| {P}_{i}|$ amount of power output) in the grid, for each node i, which will be used in the forthcoming equation describing the phase dynamics.

A power-grid node is an oscillator in terms of the phase angle of its alternating current that oscillates sinusoidally with a designated angular frequency. The oscillation dynamics is described by the motion of a damped driven pendulum, the swing equation [32], which is a Kuramoto-type model. Although the model simplifies the detailed parameters of real power grids, it can still effectively simulate the essence of the physical interactions. As a result, it is widely used in power-grid studies [3, 16, 17, 19, 23, 24, 33]. However, in contrast to the conventional (first-order) Kuramoto model where the system is completely governed by the first derivative of the phase angle [33], the swing equation involves the second derivative of the phase angle, or the time derivative of the angular frequency as presented in the following.

Applying the swing equation to power-grid nodes, the equation of motion for the phase angle θ_i of each node i is given by [22]:

$$\begin{eqnarray}&&{\ddot{\theta }}_{i}={\dot{\omega }}_{i}={P}_{i}-{\alpha }_{i}{\dot{\theta }}_{i}-\displaystyle \sum _{i\ne j}{K}_{{ij}}{a}_{{ij}}\sin ({\theta }_{i}-{\theta }_{j}),\end{eqnarray} \tag{ 1 }$$

where the source term ${P}_{i}$ is the net-power input of node i (${P}_{i}\gt 0$ if the node i is a producer and ${P}_{i}\lt 0$ for a consumer, as we introduced before—collectively, we denote the two cases by P_p > 0 and P_c < 0, respectively.); the damping term is the product of the dissipation coefficient α_i and the deviation of the angular frequency to the reference frame ${\dot{\theta }}_{i};$ K_ij is the coupling strength that corresponds to the transmission capacity of the line between nodes i and j; the adjacency matrix element, which describes the interaction or network structure, is given by a_ij = 1 if nodes i and j are connected, and a_ij = 0 otherwise. On top of the interaction substrate given by a_ij, the coupling between the phase angles of nodes i and j is given by the sine function, as in the conventional Kuramoto model [33]. Note that the phase angle of node i denoted by θ_i represents the phase difference, and ${\omega }_{i}={\dot{\theta }}_{i}$ is the difference of the angular frequency of node i, from the desired constant angular frequency called rated frequency; these reference points are set in order to ensure the vanishing average frequency deviation over all of the nodes, i.e. ${\sum }_{i}{\dot{\theta }}_{i}={\sum }_{i}{\omega }_{i}=0$.

When power-grid nodes are synchronized in a steady state, the angular frequency is equal to the rated frequency, such that ω_i = 0 for all i's, which means that all of the nodes keep the rated frequency of the power system (frequency synchronization). To characterize the stability of power-grid systems, we thoroughly investigate the response of nodes against the coupling strength change in terms of synchronization dynamics. In particular, we examine the detailed transition pattern of the basin stability observed in the basic units of power-grid systems composed of simple graph structures, which sometimes entail nonmonotonic pattern of sudden drooping.

For the sake of simplicity, we use the global coupling strength K_ij = K as the control parameter adjusted from zero to the values large enough to guarantee the power grid's recovery to its original synchrony against the given range of perturbations. We also set the dissipation coefficients α_i = α = 0.1 for all of the nodes [16, 17, 24]. For numerical integration of equation (1), we use the fourth-order Runge–Kutta method [34] and the convergence criteria [16] of the deviation of the angular frequency from the rated frequency, $| {\omega }_{i}| \lt 5\times {10}^{-2}$ for every node i, which we believe is reasonable based on the observed fluctuations data.

We use the numerical Monte Carlo approach [24] to estimate the basin stability B, our measure of node's ability to recover its synchrony. More precisely, B is the fraction of volume in the phase space composed of the phase and angular frequency disturbance as the external perturbation, from which the synchronization recovery occurs. In order to measure B, after the system reaches the synchrony, i.e. ${\dot{\theta }}_{i}={\omega }_{i}=0$ for all i, we perturb each node j by imposing a new phase and angular frequency values $({\theta }_{j}^{* },{\omega }_{j}^{* })$ chosen from the phase space of the intervals ${\theta }_{j}^{* }\in [-\pi ,\pi )$ and ${\omega }_{j}^{* }\in [-100,100]$, respectively, uniformly at random [17]. When we apply such perturbation to the system, it may recover the synchronized steady state again or fall into a limit cycle (out of synchrony) after enough time [17]. For a given number of trials of random sampling of perturbation, we measure a node's B value by the number of successful trials of synchronization recovery from the random samples of phase space $({\theta }_{j}^{* },{\omega }_{j}^{* })\in [-\pi ,\pi )\times [-100,100]$ as described above, divided by the total number of trials. We repeat this procedure for various coupling strength K to systematically investigate B as a function of K.

In terms of synchronization stability, larger values of B indicate that the system endures a wider range of phase and angular frequency disturbance, which means that the node is more stable against certain perturbations. On the other hand, a node with B ≃ 0 is almost always unable to recover the synchronized steady state, even from small perturbation. Each node has its own transition pattern of the basin stability B, but in general, B is usually an increasing function of the value of coupling strength K [17, 20]. One can intuitively understand the larger stability at larger K values because larger coupling strength enhances the interaction between nodes and such a tightly coupled system tends to absorb perturbation via their mutual collective robustness. However, it is not always true. The counterintuitive phenomena of suddenly decreased B have been witnessed from seemingly quite stable cases with B ≃ 1 even when K is increased [21]. This phenomenon, therefore, leaves characteristic peaks of B as a function of K, and implies that the power grid may suddenly become unstable unexpectedly, even with larger values of coupling strength.

In order to investigate the stability of power-grid nodes systematically in terms of basic building blocks or motifs constituting power grids, we conduct the basin stability analysis on nodes in small graph structures from simple chains up to the enumerative set of the four graphlets [31]. One may consider a national level of power grid as a single large network. In reality, however, small-size micro grids, loosely isolated from the national grid, are also important, e.g. the local power plants in industrial complexes or highly populated areas supply the regional power demand. A more specific example is the power plants for renewable energy sources such as solar photo-voltaic, which supply electricity to keep the local power grid alive, even when the main power grid becomes shut down for some reasons. These small parts of the power grid operate independently from the main power grid in the form of micro power grids, which is called 'islanding' [35]. Understanding the synchronization behavior and stability of these islands is important even in the perspective of the entire power grid stability, because the islands should maintain the local functioning of the power grid even when the blackout in the global scale occurs [5, 35]. Therefore, studying a set of small graph configurations is a good start to explore large structures [36] and, at the same time, it is also important to manage the local power stability.

We compare the transition patterns of B for nodes located at different topological positions in various structures and scales of power grids. The amount of power input from the producers is assigned as $1\leqslant {P}_{p}\leqslant 2$ depending on their topological properties, which will be presented later. The power output to the consumers, the relative portion of which depends on the consumers' topological properties as well, is set for the total amount of power input and output to be balanced as ${\displaystyle \sum }_{i=1}^{n}{P}_{i}=0$, where n is the total number of nodes.

Let us consider a very large fully-connected network as an ideal power grid, where all of the nodes are well connected so that synchronization between them never fails. In this ideal power grid, any perturbation to a power-grid node is instantly absorbed and the angular frequency of it is fixed at the rated value. Conventionally, a node of this type of ideal power grid is often replicated by an oscillator coupled with a busbar—so-called infinite busbar—and the system that can absorb an infinite amount of perturbation, analogous to the heat reservoir in the context of heat transfer in statistical mechanics [14, 23]. In practice, the phase and frequency deviations of the infinite busbar are set to be constant, such that θ_f = 0, ω_f = 0 during the simulation.

The transition pattern of the infinite busbar is the gradual and monotonic increase when the transmission strength increases, except for minor stochastic fluctuations from random perturbations in the phase space (see figure 1). It is also reported that at K values large enough to rule over the perturbation, the basin stability reaches the maximum value of unity, which means that the system can absorb any perturbation within the given range of phase and frequency space [17, 21].

Zoom In Zoom Out Reset image size

The pattern of monotonically increasing basin stability as K is increased and the appearance of the perfect stability B = 1 above a certain threshold K_c also occurs in the power grid composed of two nodes as shown in figure 2. The two nodes are a power producer and a consumer with P_p > 0 and P_c = −P_p < 0, respectively, for the balanced power input and output, i.e. P_p + P_c = 0. The transition patterns of B of the both producer and consumer are identical in this case.

Zoom In Zoom Out Reset image size

The transition shape of B according to K sustains, although we tune the parameter values. Changing P value in the two-node case just shifts the K_c after where the synchronization stability holds unity. Increasing the amount of input, for example, $| P| =1.5\to 2$ from figures 2(a) to (b), the threshold for K is also increased (K_c ≃ 278 → K_c ≃ 494), but the overall shape of the transition pattern is not altered. Based on this observation, the monotonic increase seems to be a universal characteristics of the basin stability transition, which is not always true as we present in the following.

In contrast to the monotonically increasing transition pattern of the infinite busbar, Kim et al [21] revealed the emergence of nonmonotonic transition patterns of the basin stability as K is increased in small networks. The stability measured by B becomes quite large for a short window of K, followed by a drastic decrease of B with the increased value of K, counterintuitively. In other words, there are local peaks of B as a function of K, for small power-grid networks. The nonmonotonic pattern is observable in small-scale power grids and is also related with the community characteristics [20]. However, the specific underlying mechanism of the emergence of the nonmonotonic behavior is yet to be elucidated.

Our systematic numerical simulations show that the nonmonotonic behavior in B even appears in small power-grid networks with simple topology. For instance, in the case of ${P}_{p}=1$, with additional power consumers to a single power plant, a power producer shows a sudden increase and decrease of B when it is connected to a number of consumers (see figure 3). The peak of B appears only when there are two or three consumers, as shown in figures 3(b) and (c). In these cases, even if B reaches the value close to unity for a certain K value, the value of B can suddenly droop to an intermediate value when we increase K (the cases of nodes denoted by N8, N9, and N11 in figure 3). This peak of B does not appear in the aforementioned infinite busbar and the two-node power grid.

Zoom In Zoom Out Reset image size

The transition pattern with the peak is related not only to the size of the network but also to the location of power-grid components and the structural properties of the network. For a given number of nodes, different network structures induce various transition patterns. For instance, connected networks composed of three nodes can belong to three different types of topology (up to isomorphism) as shown in figure 4, with a power producer (mint) with P_p = 2 and two consumers (red) with P_c = −1. When the power producer is located in the center, the node denoted by N27 in figure 4(a), it shows a similar transition pattern to the infinite busbar or two-nodes power grids, except for the value of the critical coupling strength (K_c ≃ 186). However, the consumers on the left and right sides (N26 and N28) have a sharp peak before the nodes reach B = 1, with much smaller K_c ≃ 46 than that of the central node.

Zoom In Zoom Out Reset image size

In this case, a local change in a topology can affect the synchronization dynamics of the entire network. Switching the positions of the central producer node and the left consumer node in figure 4(a) (N26 $\leftrightarrow$ N27) will result in the network topology in figure 4(b), which shows a completely different transition pattern of B. Interestingly, the level of instability in a certain range of $K\in [100,180]$ seems to become less severe from the left-end (N29, the producer) to the right-end (N31, one of the consumers). Moreover, when an additional edge connects both end nodes to complete the triangle (N32–N34), the peak disappears and the transition pattern become monotonic. These various transition patterns exemplify that predicting the synchronization stability is nontrivial.

Such a nonmonotonic transition pattern with the peak structure becomes more common in four-node networks. Figure 5 shows the various transition patterns of B for four-node networks [31]. In the case of connected four-node networks, 11 types of network topology are possible, when we consider the balanced network consisting of the two producers and the two consumers.

Zoom In Zoom Out Reset image size

The various characteristic transition patterns of B shown in the three-node power grids shown in figure 4 are also observed in the four-node cases. Note that the networks where the producer and consumer sides are polarized (g1, g4, and g5 in figure 5) show significant instability. For example, the nodes located between the producer and consumer side have a low level of B for the entire range of K. In figure 5, the two central nodes in g1 and the nodes with three neighbors in g4 and g5 are such cases. The observation is in good agreement with the report that centralized power grids are more vulnerable than distributed ones [18].

We also find the nonmonotonic behavior in some nodes in topologies that have four nodes and five edges in figure 5, including the fully-connected one (g11). As a specific example, figure 5(g) shows an enlarged view of the transition pattern of the bottom left node of g10 in figure 5(c). The local peak persists with a very high value of the basin stability for $13\lesssim K\lesssim 15$, before it decreases again for $15\lesssim K\lesssim 18$. The basin stability reaches a final transition point to unity at $K\gtrsim 25$. The onset of the local peak is much sharper than its decay, which is also observed for other topologies such as the three-node power gird in figure 4(b) (see the plot for N29).

In the four-node networks denoted by g9, g10, and g11 in figure 5(c), the power producers and consumers show very similar transition patterns in B as K increases. The local peak commonly appears in the range of $13\lesssim K\lesssim 19$ for both power producers and consumers. The overall transition shapes show a similar pattern for all nodes. The global basin stability, defined as the value of B averaged over all of the nodes [28] and denoted by $\langle B\rangle$ in figures 5(e) and (f), shows that such drastically different patterns for individual nodes cannot always be observed by just looking at the global basin stability.

In order to investigate this stability drooping in more detail, we investigate the variations in the instantaneous frequency of each node after the transient period, when we apply the same perturbation to a specific node. Since the basin of attraction in phase space changes as the coupling strength K increases, the same point in phase space may not lead to the same attractor for different values of K. Here, we carefully choose an external perturbation to a target node such that all of the four nodes fail to recover synchronization for various values of K, except between K = 13 and 14. In this range, all nodes show B ≃ 1 implying that they almost always return to the synchronized state.

As a specific example, we focus on the four-node network denoted by g10 in figure 5 and shown in figure 6(a). The producer nodes are located on the left (red and orange) and the consumer nodes on the right (blue and mint). The external perturbation is applied to the red producer node with ω_red = 16.2. After a sufficiently long time (t > 495) to minimize any transient effects, we observe the dynamics of the instantaneous frequency as shown in figure 6(b) for various K values. While for K = 13 and K = 14 we recover a synchronized state of ${\omega }_{i}=\dot{{\theta }_{i}}=0$ as expected, the frequency and phase of the nodes orbit a limit cycle and do not converge to the synchronized state in all other cases. More interestingly, comparing K = 12 and K = 15, we observe a period-doubling phenomenon in the emergent oscillations of two of the nodes—the period of the producer node (orange) and the consumer node (mint) becomes twice the period of the other nodes. When K = 15, just after the stability drooping, these two nodes are also in anti-phase with each other, while they were in-phase and identical for K = 12. However, as the coupling strength K increases, we observe a reversed period-doubling and for K = 18 the dynamics becomes qualitatively identical to the case of a simple limit cycle at K = 12: The orange and mint nodes have identical dynamics and they are out of phase with the other nodes, while satisfying the condition ${\sum }_{i}{\omega }_{i}=0$.

Zoom In Zoom Out Reset image size

In order to understand the observed response of the oscillator network to the chosen perturbation better, we now explicitly explore the emergent dynamics in a more extended part of phase space (θ, ω) for K = 12, 13, ⋯, 19, focusing on initial conditions θ ∈ [−π, π) and ω ∈ [−10, 50] for the red node in figure 6(a), while the initial conditions of all other nodes are (0, 0). Therefore, this corresponds to a 2-dimensional subspace of the full 8-dimensional phase space, but it allows for a better visualization. Figure 7 displays eight plots of the basins of attraction for different K values, where each light green point in phase space corresponds to an initial condition leading to the synchronized state. The fraction of the light green area in each plot in figure 7 approximately corresponds to the basin stability B for the given K.

Zoom In Zoom Out Reset image size

When the coupling strength K = 12 in figure 7, more than half of the area of the shown phase space makes the system fail to recover synchronization, corresponding to the dark blue area for $\omega \gtrsim 10$, which results in $B\lesssim 0.5$ in figure 5(g). However comparing figures 7 and 5(g) further, when K = 13, the basin stability significantly increases and the dark blue area breaks down into small islands. This tendency continues when we increase K to K = 14, which results in B ≃ 1. Then, the basin stability B suddenly starts to decrease when we increase K to K = 15. As shown in the plot for K = 15 in figure 7, the number of the dark blue speckles increases again and in addition and more importantly a solid blue area or channel around ω ≃ 10 emerges in phase space. Interestingly, the speckled area of dark blue, corresponding to the values of perturbation from which the node fails to recover synchronization, forms a moiré-like pattern for K = 16 and K = 17. Eventually, around the top half of the phase space turns into the basin of attraction for the unsynchronized state (see K = 18 and K = 19 in figure 7.) These speckled patterns and the moiré-like fine structures of the basin boundary provide a more detailed insight into the stability droop.

In order to analyze the observed multistability and its relationship to the changes in basin stability in more detail, we now focus on how the stability of the different limit cycles shown in figure 6(b) changes with K. The simple limit cycle observed for K = 12 remains stable over a range of K values. As shown in figure 8(a), at K ≃ 12.555 the limit cycle suddenly becomes unstable and decays to the fixed point corresponding to the synchronized state. This change in stability coincides with the onset of the local peak in basin stability (see figure 5(g), where the first vertical yellow line corresponds to the instability of the simple limit cycle) indicating that the basin of attraction of the formerly stable limit cycle now almost exclusively belongs to the basin of attraction of the synchronized state.

Zoom In Zoom Out Reset image size

In contrast, the period-doubled limit cycle observed for K = 15 in figure 6(b) loses its stability for decreasing K. As shown in figure 8(b), this limit cycle becomes unstable for K ≃ 14.61 and the synchronized state is the new attractor. This bifurcation point of the limit cycles is indicated by the second yellow vertical line in figure 5(g). The instability coincides with the initial decay of the local peak suggesting a direct connection between the basin of attraction of the period-doubled limit cycle and the basin stability. Both instabilities of the two different limit cycles leading to the synchronized state as the new attractor share the property that the frequency amplitudes of the limit cycle remain finite up to the respective instability as shown in figures 8(a) and (b). This suggests that the type of bifurcation is the saddle-node bifurcation of cycles in both cases [37].

Interestingly, as we increase K gradually as shown in figure 8(c) the evolution of the period-doubled limit cycle follows the same pattern as observed in figure 6(b) for $K\geqslant 15$. Specifically, the period-doubled limit cycle first undergoes an inverse period-doubling transition around K ≃ 15.6 and later at K ≃ 17.3 a transition to the simple limit cycle, which is identical to that of K = 12. These two transition points are indicated by red vertical lines in figure 5(g), which are consistent with the more slowly decaying tail of the local peak in B.

Since both the simple limit cycle and the period-doubled limit cycle are unstable between K ≃ 12.555 and K ≃ 14.61—which is consistent with the range of K values in figure 7 for which no extended blue areas are visible—the initial conditions for which synchrony is not recovered for this range of K values must correspond to another attractor. Indeed, there exists another stable limit cycle as shown in figure 9. This limit cycle has different characteristics from the others. Both producers nodes are in-phase with each other and the two consumer nodes are in-phase with each other as well. Yet, producer and consumer nodes are in anti-phase with each other. The period of the limit cycle is almost half of the former stable limit cycle in figure 8(a). In order to investigate its stability, we decrease and increase K values continuously as shown in figure 9. This limit cycle shows very stable behavior for a wide range of K. The period does not change with K, but its amplitude does.

Zoom In Zoom Out Reset image size