Prerpints.org logo
Preprint
Article

This version is not peer-reviewed.

System-Inherent Grid Synchronous Stability Boundary and Spontaneous Synchronization

Yu Yuan  *

Submitted:

31 October 2023

Posted:

02 November 2023

Read the latest preprint version here

Abstract
Synchronization is a prerequisite for the stability of electrical energy, which is the cornerstone of a functioning society. Identifying whether generators are synchronized or not has been attracting immense research interest. Utilizing synchronized stability boundaries can strongly simplify stability discrimination. However, fundamental issues such as determining accurate stability boundaries remain challenging. In this study, an accurate synchronized stability boundary equation that is independent of the system and disturbances is presented, and it is applicable to an arbitrarily coupled grid. The special structure induced by spontaneous synchronization appears to be only near the stability boundary. The results of this study show that: the synchronous stability boundary is an intrinsic property of the system, independent of its structure and parameters. Synchronization on the network may be independent of the network. These conclusions will challenge the traditional perception of synchronization on networks, and the synchronization stability of the grid can be analyzed in a uniform way.
Keywords: 
Power system
;  
spontaneous synchronization
;  
synchronous stability boundary
Subject: 
Engineering  -   Electrical and Electronic Engineering

Introduction

Synchronous stability is a prerequisite for the normal operation of a power grid4. The loss of stability in a power system can cause power supply disruptions for a wide range of customers. Large power systems are complex coupled systems where uncertainties5 and nonlinearities coexist. In the stability analyses of power grids, a high level of real-time performance is required6. To analyze the stability of these systems, researchers have developed numerous insightful stability discriminatory methods7–12, including finding stability regions and their boundaries2,3,13.
Finding the stability boundary of a system is an important fundamental problem3,14,15. The boundary, which can be used to strongly simplify stability discrimination, is the union of unstable equilibrium points (UEPs)16. When an operating point is outside the boundary, the corresponding system is unstable13. The synchronous stability boundary, which is a core concept of grid stability, is closely related to many issues6,17,18. Therefore, studying synchronous stability boundaries has a significant impact on the development of power systems. For decades, scientists and engineers have wanted to find an analytical equation to describe the boundary2,3,19. However, researchers continue to face challenges, such as high computational complexity, difficulty of visualization, and difficulty in adapting to system uncertainty.
In addition, network concepts have been utilized to explain the synchronization of generators by studying the spontaneous synchronization conditions of various networks20. However, the reasonableness of this assumption, which directly equates the spontaneous synchronization with the synchronization stability, has been questioned9,21 and challenges remain in applying it to real networks8. One way to clarify this debate is to obtain data through traditional power system research methods and to obtain evidence of spontaneous synchronization from them. However, such results are rarely reported9. Electrical engineering communities need an explanation with clearer physical meaning.
Here, a simple equation based on a coordinate transformation and new concepts was derived and visualized to describe the stability boundary of a power system. The derivation of the equation is based entirely on physically based concepts. The form of the equation is independent of the network topology, component parameters and fault disturbances and is therefore very generalizable. At the same time, the synchronization stability boundary represented by the equation is independent of the network and the perturbations. This equation can be used to analyze the synchronization stability of the grid in a uniform way. The critical clearing time (CCT) and the generator critical rate were calculated by fitting the disturbed trajectory using the same variable system as the boundary. Additionally, simulation results on fine scales from the IEEE standard arithmetic models demonstrated evidence of spontaneous synchronization1, which occurred near the synchronization stability boundary. This result may also hint that synchronization on the network may be independent of the network. This will challenge the traditional perception of synchronization in networks. Directly discriminating synchronization on the network will be simplified to the relationship between the network and the operating point coordinates. Determining the stability of a coupled grid containing n generators requires only n pairs of variables. These variables have a clear physical meaning and are easily accessible. Spontaneous synchronization and synchronous stability have been widely discussed in many fields8,20. In this case, for potential scenarios where it is difficult to construct a network model, the study demonstrates the success of directly determining synchronization on the network.

Materials and Methods

Extensive interconnections between generators would make stability analysis very difficult(see Extended Figure 1). To solve this problem, the concept of a meta-generator is introduced here. At moment t, the instantaneous values of the n generators system （ u ′ i ( t ) , δ ′ i ( t ) ) , i ∈ ( 1 , 2 , ... , n ) are arranged in descending order by δ ω , relabelled, and then reconstituted as the n meta-generator system （ u i ( t ) , δ i ( t ) ) , i ∈ ( 1 , 2 , ... , n ) .
Preprints 89296 g001
Figure 1. Schematic diagram of the power system operating coordinate system （ U , δ ω ) . The generator variables are selected from the port bus voltage per unit U and the rotor speed per unit ω , U ∈ [ 0 ， + ∞ ） ， ω ∈ （ − ∞ ， + ∞ ） . The coordinates of ω are transformed: δ ω = 2 arctan ( ω ) , δ ω ∈ ( − π , π ) , where δ ω is the angle of rotation rate. The angle of rotation rate of the ith generator δ ω i is abbreviated as δ i in the following. Note: δ ω is not the power angle. a). Synchronised system operation before disturbance, ω 1 = ω 2 = ⋯ = ω n = 1 and δ i = 2 arctan ( ω i ) = π 2 , i = 1 , 2 , ⋯ , n . The magenta dot indicates the operating point of the generator: （ U , δ ω ) . b) The angle of rotation rate of the subsystems differ after the disturbance. u K , u L are the per unit voltage of the port bus of the Kth and Lth generators, respectively, δ K L = δ K − δ L ≥ 0 is the difference in the angle of rotation rate between the Kth and Lth generators. Δ u is defined as the coupling potential difference between generators Kth and Lth (yellow dotted line between Cyan dots). Correspondingly, u ′ is constructed to describe the synchronous potential difference between generators Kth and Lth (solid blue line between magenta square dots). The set of points where u ′ = Δ u is the synchronous stability boundary (see Methods for details: Derivation of the boundary equations).
Figure 1. Schematic diagram of the power system operating coordinate system （ U , δ ω ) . The generator variables are selected from the port bus voltage per unit U and the rotor speed per unit ω , U ∈ [ 0 ， + ∞ ） ， ω ∈ （ − ∞ ， + ∞ ） . The coordinates of ω are transformed: δ ω = 2 arctan ( ω ) , δ ω ∈ ( − π , π ) , where δ ω is the angle of rotation rate. The angle of rotation rate of the ith generator δ ω i is abbreviated as δ i in the following. Note: δ ω is not the power angle. a). Synchronised system operation before disturbance, ω 1 = ω 2 = ⋯ = ω n = 1 and δ i = 2 arctan ( ω i ) = π 2 , i = 1 , 2 , ⋯ , n . The magenta dot indicates the operating point of the generator: （ U , δ ω ) . b) The angle of rotation rate of the subsystems differ after the disturbance. u K , u L are the per unit voltage of the port bus of the Kth and Lth generators, respectively, δ K L = δ K − δ L ≥ 0 is the difference in the angle of rotation rate between the Kth and Lth generators. Δ u is defined as the coupling potential difference between generators Kth and Lth (yellow dotted line between Cyan dots). Correspondingly, u ′ is constructed to describe the synchronous potential difference between generators Kth and Lth (solid blue line between magenta square dots). The set of points where u ′ = Δ u is the synchronous stability boundary (see Methods for details: Derivation of the boundary equations).
Preprints 89296 g001
The essence of the connection between the meta-generators and the original generators is the permutation: S = δ ′ 1 δ ′ 2 ... δ ′ i δ ′ i + 1 ... δ ′ n δ 1 δ 2 ... δ i + 1 δ i ... δ n . When the condition δ ′ i + 1 ( t ) > δ ′ i ( t ) holds, （ u i + 1 ( t ) , δ i + 1 ( t ) ) = （ u ′ i ( t ) , δ ′ i ( t ) ) , （ u i ( t ) , δ i ( t ) ) = （ u ′ i + 1 ( t ) , δ ′ i + 1 ( t ) ) .
Define the distribution matrix H, where H i ′ j is the share of the i'th generator in the jth meta-generator in period T. H can be a potential choice for generator coherent group prediction techniques.
Data sources and experimental procedures
In this study, the IEEE 39-BUS and IEEE 9-BUS models (Extended Figure 4, Figure 5 and Figure 6) were used. The two models are simulated separately using a simulation software package. Here, the fault was set as a three-phase short circuit to ground. The disturbed operating point of each generator was calculated. To observe the movement pattern of the disturbed operating points, no parameters were set for the control elements .
In advance, the fault location was fixed, and the fault duration Δ t was set. This experiment simulated the rotation rate ω ′ i ( t ) and port bus voltage u ′ i ( t ) of the ith generator after different disturbances. Then, the angle of rotation rate of the ith generator δ i ( t ) was calculated. Δ t was increased in a fixed step length and u ′ i ( t ) , δ i ( t ) were calculated again until the system was destabilised. The faulty position was replaced, and the above steps were repeated.
Subsequently, （ u ′ i ( t ) , δ ′ i ( t ) ) , i ∈ ( 1 , 2 , ... , n ) was arranged and relabelled as （ u i ( t ) , δ i ( t ) ) , i ∈ ( 1 , 2 , ... , n ) . This was then averaged as follows:
The mean of （ u i ( t ) , δ i ( t ) ) , i ∈ ( 1 , 2 , ... , n ) over 0 , T was found: u i = 1 T ∫ 0 T u i ( τ ) d τ and δ K L = 1 T ∫ 0 T δ K L ( τ ) d τ = 1 T ∫ 0 T δ K ( τ ) − δ L ( τ ) d τ = 1 T ∫ 0 T δ K ( τ ) − δ L ( τ ) d τ = δ K − δ L .
The mean of （ u i ( t ) , δ i ( t ) ) , i ∈ ( 1 , 2 , ... , n ) over T m , T m + d T was also found: δ i ( t m ) = 1 d T ∫ T m T m + d T δ i ( τ ) d τ , U i ( t m ) = 1 d T ∫ T m T m + d T U i ( τ ) d τ , δ k l ( t m ) = δ k ( t m ) − δ l ( t m ) . There are several definitions of mean, and the simplest, i.e., the arithmetic mean, was used here.
This work added adjacent meta-generator data ( u K , u L , δ K L ) and ( u K ( t m ) , u L ( t m ) , δ K L ( t m ) ) to the coordinate system u K − u L − δ K L to assess the system stability (Figure 2b) and time intervals of instability T m , T m + d T (Figure 2d).
An expression was fitted with Δ t as the independent variable and u K , u L , δ K , δ L , δ K L as the dependent variable (Figure 3). CCT and UEP were then calculated.
Near the boundary, σ δ i , δ i , u i was calculated at a finer scale.
Derivation of the boundary equation
As shown in Figure 1, Δ u = u k 2 + u l 2 − 2 u k u l cos δ k l , u ′ = 2 u k 2 ( 1 − cos δ k l ) .
Following the form of power in electricity P = U 2 R , the coupling power P Δ u = Δ u 2 Z K L is defined to characterise the coupling between the meta-generators. To describe the energy required for the generator to maintain synchronous stability, the synchronous power is constructed: P u ′ = u ′ 2 Z K L . Z K L is the impedance between the Kth and Lth meta-generators.
When the system is synchronized, the meta-generators are not in balance and are still coupled ( δ K L = 0 , u ′ = 0 , Δ u = u K − u L ≠ 0 ). When the system is disturbed, δ K L increases from 0, u ′ increases from 0, and Δ u changes. When the coupling power between the two meta-generators is sufficient to provide synchronous power, i.e., P u ′ ≤ P Δ u , the system is synchronous and stable. Conversely, when P u ′ > P Δ u , there is not enough coupling power to maintain synchronization, and the system is unstable. It is observed that P u ′ P Δ u = 1 ⇔ u ′ Δ u = 1 .
In summary, f u K , u L , δ K L = u ′ Δ u = 1 is the system stability boundary equation. When u ′ Δ u < 1 , the system is stable. When u ′ Δ u > 1 , the system is unstable. Geometrically, f u K , u L , δ K L = 1 describes a curved surface that, together with 0 , u L , δ K L , u K , 0 , δ K L , u K , u L , 0 , encloses a stable domain. In summary, the boundary equation u K = u L ∪ u L u K = 2 cos δ K L − 1 can be found, where u K ≥ u L ≥ 0 , δ K L ∈ [ 0 , π 3 ] . The coordinate system u K − u L − δ K L is established, and the boundary is visualized(Figure 2).
Fitting of operating points to trajectories
The intersection of the disturbed trajectory with the stability boundary is the unstable equilibrium point (UEP), and the failure time of the operating point along the disturbed trajectory to reach the UEP is the CCT. To calculate these important results, it is necessary to fit the disturbed trajectory to the kinematic expression in the （ u , δ ω ) coordinate system. The variables of the operating point obtained from the simulation are fed into commercial software to fit the expression of the disturbed operating point.
Operating point behaviour on the boundary
d δ K L d ( Δ t ) is the derivative of δ K L of the meta-generators with respect to Δ t . Near the boundary, the derivative of d δ K L d ( Δ t ) of the partial meta-generators changes from positive to negative (Figure 3). For this unusual phenomenon, on a finer scale, （ u i ( t ) , δ i ( t ) ) , i ∈ ( 1 , 2 , ... , n ) is calculated sequentially for different Δ t . The standard deviation of δ i is calculated separately (Figure 4): σ δ i = ∑ i = 1 n ( δ i − μ δ i ) 2 n , μ δ i = ∑ i = 1 n δ i n .

Result and Discuss

Stability boundary
Preprints 89296 g002
Figure 2. Stability boundary. a). Visualisation of the stable boundary. The stability boundary (blue surface) and 0 , u l , δ k l , u k , 0 , δ k l , u k , u l , 0 (pink surface) together enclose the stability domain. The boundary equation is u k = u l ∪ u l u k = 2 cos δ k l − 1 (see Methods for details: Derivation of the boundary equations). For visualisation purposes, the radian system is used for δ k l in the figure. b) and c) are the results of the 12BUS three-phase short circuit to ground simulation, respectively. The system is stable at Δ t = 0.28 s , unstable at Δ t = 0.3 s , and unstable in 6 s , 7 s . d). Boundary stability results for the 12BUS three-phase short circuit to ground, where Δ t is the duration of the fault. The positions of all operating points at Δ t = 0.28 s and Δ t = 0.3 s are indicated by cyan and magenta dots, respectively. The plane u l = u k perpendicular to u k − u l is not shown. At Δ t = 0.28 s , all the operating points are clustered together. However, Δ t = 0.3 s , the operating points are far away from the other points and outside the boundary. e). The instability process of Multiswing for the 12BUS three-phase ground fault, where Δ t = 0.3 s . Numbers 1 to 9 indicate period T m after the fault, and d T = 1 s (see Methods for details). The system is destabilised in the time interval 6 s , 7 s (Figure 2). The diagram shows the operating points of the meta-generators 1 and 2. The operating point crosses the boundary at the time of instability, and δ k l rapidly increases after the voltage changes.
Figure 2. Stability boundary. a). Visualisation of the stable boundary. The stability boundary (blue surface) and 0 , u l , δ k l , u k , 0 , δ k l , u k , u l , 0 (pink surface) together enclose the stability domain. The boundary equation is u k = u l ∪ u l u k = 2 cos δ k l − 1 (see Methods for details: Derivation of the boundary equations). For visualisation purposes, the radian system is used for δ k l in the figure. b) and c) are the results of the 12BUS three-phase short circuit to ground simulation, respectively. The system is stable at Δ t = 0.28 s , unstable at Δ t = 0.3 s , and unstable in 6 s , 7 s . d). Boundary stability results for the 12BUS three-phase short circuit to ground, where Δ t is the duration of the fault. The positions of all operating points at Δ t = 0.28 s and Δ t = 0.3 s are indicated by cyan and magenta dots, respectively. The plane u l = u k perpendicular to u k − u l is not shown. At Δ t = 0.28 s , all the operating points are clustered together. However, Δ t = 0.3 s , the operating points are far away from the other points and outside the boundary. e). The instability process of Multiswing for the 12BUS three-phase ground fault, where Δ t = 0.3 s . Numbers 1 to 9 indicate period T m after the fault, and d T = 1 s (see Methods for details). The system is destabilised in the time interval 6 s , 7 s (Figure 2). The diagram shows the operating points of the meta-generators 1 and 2. The operating point crosses the boundary at the time of instability, and δ k l rapidly increases after the voltage changes.
Preprints 89296 g002
Regarding the synchronization stability, one different assumption,which appears parallel to the Kuramoto model8, than before 18: the system requires the synchronous power P u ′ = u ′ 2 Z K L to maintain synchronous stability, and the synchronous power is provided by the coupling power P Δ u = Δ u 2 Z K L within the system, Z K L is the impedance between the Kth and Lth meta-generators. (see Supplementary Materials for details).
Preprints 89296 i001
Equation (1) is the analytical equation for the synchronous stability boundary. Specifically, when δ K L > 0 , the set where u L = u K is the isolated stability domain. There has been previous research on the isolated domain2. Since it is very difficult to always maintain u L = u K after a disturbance of the system, this situation is not discussed in this paper.
In contrast to previous reports2,3,20, here, The synchronous stability boundary described by Equation (1) is independent of the network topology, system parameters, perturbations and the number of the subsystem. It is an inherent property of the grid systems. Regardless of the changes in the network structure and component parameters, it is only the operating points of the system that change rather than the stability boundary. This feature increases the applicability of the boundary equation to different grid systems and allows the stability of multimachine systems to be analyzed even in scenarios where the topology and parameters are not clear22, such as in a black box situation. Thus, the synchronization stability of the grid can be analyzed in a uniform way. This has been validated on different standard arithmetic models (Figure. 2 and Extended Figure 4). Equation (1) is also suitable for identifying the stability of a power system with multiple swings23 and showing the time interval of instability(as in Figure 2c,e).
Equation (1) has a variant as follows:
δ K L c r = arccos ( 1 − u K − u L 2 u L )
When u K , u L are sufficiently close24, the stability margin of the Kth and Lth generator angle rate difference δ K L c r also tends to 0. This indicates that the system may already be in a critical state during normal operation. In this case, even if the difference in the values between δ K and δ L is small, it takes only a very small perturbation to make the system unstable25,26. This may explain why some systems with very small speed differences are more likely to lose synchronization than those with much larger speed differences.
Placing the calculated or measured operating points in the coordinate system u K − u L − δ K L can help directly determine whether the system is stable. A system is considered unstable as long as at least one point ( u K , u L , δ K L ) is outside the boundary (Figure 2b) and far from the other operating points ( σ Δ t δ i > > σ Δ t − ε δ i , 0 < ε < Δ t , where Δ t is the fault duration) (Figure 2b and Figure 4a). To determine the stability of a power system of n generators, only n pairs of variables u i , δ i , i ∈ ( 1 , 2 , ... , n ) are needed, which are physically meaningful and easily obtainable.
As one of the concepts closely related to stability, intentional isolation is an effective way of avoiding widespread blackouts following instability27. Identifying coherent generators is a prerequisite for building intentional islands28. Information regarding the coherence meta-generator groups can be obtained directly from the graph. That is, the partial synchronization of multiple generators, such as chimeric states29–31, can also be identified with a stable boundary. m operating points outside the boundary indicate that the n meta-generators are sequentially divided into m+1 coherent groups.
Disturbed trajectory and parameters
To understand the behavior of the power system after a disturbance and to calculate the CCT and UEP, it is necessary to study the perturbation trajectories consisting of running points.
Preprints 89296 g003
Figure 3. Preprints 89296 i004
Figure 3. Preprints 89296 i004
Preprints 89296 g003
Figure 3a,b show that the operating points move with a uniformly variable speed before the system becomes unstable.
Contrary to intuition32, δ 12 suddenly dropped at Δ t = 0.28 s (Figure 3c). This anomaly suggested that the system appears to have a tendency to maintain its own stability.
Figure 3d and E show that the perturbed trajectories of the subsystems of the coupled system are linearly correlated in the stability domain. This indicates that for a determined power grid, each perturbed trajectory has the same T ′ i , T ′ i being the global invariant of the system. The effects of perturbations are global, reflecting the challenges of controlling the stability of complex systems33–35.
When a high degree of accuracy of the results is not needed, a K L = 0 . The following expression can be derived:
Preprints 89296 i002
This can be used to easily and quickly check the stability margin of the system after a disturbance36. By approximating Δ T ′ as the CCT19,36, the coordinates of the critical stable operating point ( u K c r , u L c r , δ K L c r ) and critical rate of the meta-generator l in the current system can also be estimated37.
Preprints 89296 i003
In summary, the CCT and UEP can be calculated using only information about the rotation rate12, but considering only a single information source may result in more errors. Theoretically, using Δ T ′ directly as the CCT would also lead to a conservative result.
The current power system is receiving an increasing number of renewable energy sources. These sources are connected to the power grid via inverters, which may change the inertia of the system38–40, complicating the coefficients. This issue should be further studied.
On a finer scale, the operating points near the stability boundary exhibit unusual behavior (Figure 4).
Preprints 89296 g004
Figure 4. Boundary effects of the 18bus three-phase short circuit to ground fault, with increased fault duration Δ t from 0.140s to 0.154s ( d ( Δ t ) = 0.001 s ), and the trajectory of the disturbed operating point near the boundary. The arrow shows the direction of increase of Δ t . The simulation result showed instability at Δ t =0.155s. a). σ δ i is the standard deviation of δ. σ δ i started to fall at 0.147s and rose by 1300% at 0.155s when the system became unstable. b). In the δ 1 − δ 2 plane, operating points appeared to cross the barrier before they reached the boundary. The elliptical area marks the position of the barrier. From 0.147s onwards the interval between operating points decreased in the direction of increasing Δ t (shaded area). c). In the u 1 − u 2 plane, the graph is a critical state local attractor, which appears simultaneously with the synchronous barrier. The ellipse indicates the position of the attractor. The graph of the trajectory of the operating point from 0.147s onwards is shown as an attractor (in the shaded area).
Figure 4. Boundary effects of the 18bus three-phase short circuit to ground fault, with increased fault duration Δ t from 0.140s to 0.154s ( d ( Δ t ) = 0.001 s ), and the trajectory of the disturbed operating point near the boundary. The arrow shows the direction of increase of Δ t . The simulation result showed instability at Δ t =0.155s. a). σ δ i is the standard deviation of δ. σ δ i started to fall at 0.147s and rose by 1300% at 0.155s when the system became unstable. b). In the δ 1 − δ 2 plane, operating points appeared to cross the barrier before they reached the boundary. The elliptical area marks the position of the barrier. From 0.147s onwards the interval between operating points decreased in the direction of increasing Δ t (shaded area). c). In the u 1 − u 2 plane, the graph is a critical state local attractor, which appears simultaneously with the synchronous barrier. The ellipse indicates the position of the attractor. The graph of the trajectory of the operating point from 0.147s onwards is shown as an attractor (in the shaded area).
Preprints 89296 g004
As the stability boundary is approached, the perturbed trajectories of the operating points become interesting. The system self-organized and moved toward synchronous evolution1.
The decrease in σ δ i from the highest point indicated a tendency for the system to maintain its own stability before destabilizing and to spontaneously lead the velocities of the subsystems to the mean value. This may be the spontaneous synchronization effect caused by the coupling of the systems (Figure 4). Spontaneous synchronization is precisely a self-organizing behavior, i.e., a phenomenon whereby initially unsynchronized coupled subsystems evolve toward synchronization7,41. Spontaneous synchronization seems to occur near the synchronization stability boundary (Extended Figure 7). For coupled network systems, this strong correlation indicates that the location where spontaneous synchronization occurs is also determined in u K − u L − δ K L by Equation (1)1,20. This may indicate that the mechanism of spontaneous synchronization is not necessarily related to the network topology, i.e., synchronization on the network may be independent of the network42–44. This will challenge the traditional perception of synchronization in networks. At the same time, the correlation may also indicate that physically, the synchronous stability boundary may originate from spontaneous synchronisation effects.
The generators spontaneously exchanged energy through coupling to synchronize and stabilize the system. This caused the coefficients in the fitted representation to change, which was reflected in the disturbed trajectory. When the conditions were correct ( d 2 δ L d ( Δ t ) 2 = a ω ( i ) < 0 and a δ is constant), a wonderful structure emerges from the trajectory of the operating points. Due to the constraint δ K = a δ * δ L + b δ , Δ t ∈ 0 , Δ T ′ and the same Δ T ′ , the perturbed trajectories of all meta-generators simultaneously exhibit this structure. Although the phenomenon of self-organization of synchronization is often used directly to explain the synchronous operation of generators, this structure has rarely been reported in the past. The synchronous barrier and the cycle were the results of the self-organizing behavior in the δ K − δ L and u K − u L planes, respectively (Figure 4). The perturbed trajectory of the generator shows the same result at the same time (Extended Figure 3), proving that this is not caused by a substitution effect but by the emergent nature of the system at the boundary
After the operating point crosses the potential barrier, the energy is concentrated in a very short time in some of the generators, and σ δ i rises sharply (Extended Figure 4). The operating point then crossed the boundary and destabilized the system.
This self-organization phenomenon on the boundary increases the error in the CCT and UEP. This will be explored further in the future.
The results of the self-organizing behavior are complex and do not always result in the formation of a synchronous potential barrier (Extended Figure 6). The reasons for this difference, or rather, the specific conditions for the formation of this particular structure of synchronized barrier, have yet to be investigated.

Conclusions

In this study, a graceful stabilizing boundary equation that accurately describes an ideal synchronous stabilizing boundary is derived. Since the physical quantities in the equation are independent of the network, it is universally applicable to almost any power network. Therefore, the synchronization stability of the grid can be analyzed in a uniform way. It is also shown in a seminal way that the mechanism of synchronization may not be linked to the network topology and parameters. The experimental data are derived from simulations of the IEEE standard arithmetic models. The self-organizing behavior at the operating point demonstrates the existence of spontaneous synchronization on the boundary of the synchronization stability domain, which helps to confirm the argument that spontaneous synchronization is directly equivalent to synchronization stability. Additionally, the self-organizing behavior suggests the existence of a new explanation for the origin of the grid synchronous stability boundary. The concise mathematical tools, simple and universal methods, and ability to assess synchronization stability by simply monitoring the voltage and angular velocity provide great convenience for engineering applications. It is also demonstrated that synchronous stability studies of other coupled systems without clear network details can be performed by finding the correct parameters to directly derive the stability boundary while eliminating the need to construct elaborate network models. Thus, synchronization stability analyses in other disciplines may be able to directly apply the assumption made in the manuscript, with the exact form of Equation (1) depending on the form of the synchronous power and the coupling power in each discipline. Finally, the behavior of the real system's operating points near the boundary is still very complex. It is still difficult to accurately predict these results, such as the specific conditions and duration of the occurrence of synchronous barrier structures, which require further exploration.

Data availability

All the data that support the findings of this study are available at Figshare (DOI:10.6084/m9.figshare.23585961).

Conflicts of Interest

The author declare no conflicts of interest.
Extended Figure 1. Stability boundary for generators
Preprints 89296 i005
a) and b) correspond to IEEE 39-BUS system, generator disturbed instability process Schematic diagram of the disturbed trajectory of the operating point crossing the boundary after a three-phase short-circuit ground fault at 12bus and 18bus in node IEEE39 with gradually increasing ∆t, respectively. (∆t increases from 0.04s, until destabilization)
Extended Figure 2. Fitting results for generator perturbed trajectories of 12bus after a three-phase short circuit to ground fault
Preprints 89296 i006
Preprints 89296 i007
Extended Figure 3. IEEE 39-BUS system, self-organizing behavior of generators near the boundary
Preprints 89296 i008
The 18bus three-phase short-circuit ground fault, increase in fault time Δ t from 0.140s to 0.154s, and the trajectory of disturbed operating point near the boundary. The arrow shows the direction of increase of Δ t . The generator has the same barrier and attractor as the meta-generator. a). In the δ 32 − δ 34 plane, running points appear to cross the barrier before they reach the boundary, and the elliptical area marks the position of the barrier. From 0.147s onwards the interval between running points decreases in the direction of increasing Δ t . b). In the u 32 − u 34 plane, the graph is presented as a critical state local attractor, which appears simultaneously with the synchronous barrier. The ellipse indicates the position of the attractor. The graph of the trajectory of the run point from 0.147s onwards is shown as an attractor (in the shaded area).
Extended Figure 4. Stability boundary for meta-generators for IEEE 9-BUS system
Preprints 89296 i009
A). The same boundary equation applise to the 3-machine 9-node system. 4bus sets the location of the operating point after a three-phase short-circuit ground fault.
B). Results of the multiswing instability for Δt=0.25s.
C) and d are the results of the 4BUS three-phase ground fault simulation, respectively. It is stable at Δ t = 0.24 s , unstable at Δ t = 0.25 s , and unstable in 3 s , 4 s
Extended Figure 5. Fitting results for meta-generator perturbed trajectories of 4bus after a three-phase short circuit to ground fault
Preprints 89296 i010
Fitting results for perturbed trajectories for the 4bus trajectory of disturbed operating points after a three-phase short-circuit ground fault. a). The projection of the disturbed trajectory in the δ 1 − Δ t plane. The result of the fit is a ω 1 2 = 263 . 16677 ± 27 . 27487 , b ω 1 = - 12 . 16449 ± 5 . 45497 , c ω 1 = 90 . The adjusted R-squared value is 0.99997. b). The projection of the disturbed trajectory in the u 1 − Δ t plane. The result of the fit is a u ( 1 ) 2 = - 4 . 32887 ± 0 . 74225 , b u ( 1 ) = 0 . 63292 ± 0 . 21231 , c u ( 1 ) = 1 . 00743 ± 0 . 01298 , and the adjusted R-squared value is 0.97679. c). The projection of the disturbed trajectory in the δ 12 − Δ t plane. The result of the fit is a 12 = - 0 . 01956 ± 0 . 03364 , b 12 = 0 . 06242 ± 0 . 00673 , and the adjusted R-squared value is 0.99519. d). The projection of the disturbed trajectory in the δ 1 − δ 2 plane. The result of the fit is a δ = 1 . 05168   ± 0 . 01287 , b δ = - 4 . 40363 ± 1 . 21501 , and the adjusted R-squared value is 0.99925. e). The projection of the disturbed trajectory in the u 1 − u 2 plane. The result of the fit is a ν = 0 . 98731 ± 0 . 00722 , b ν = 0 . 01791 ± 0 . 00713 , and the adjusted R-squared value is 0.99973.
Extended Figure 6. Spontaneous synchronization behaviour of running points near the boundary
Preprints 89296 i011
a). The standard deviation of δ decreases from ∆t=0.227s. It rises by 2200% at ∆t=0.245s, where there is an instability. b). Derive d δ k l d ( Δ t ) = d ( ( a δ − 1 ) δ l + b δ ) d ( Δ t ) = ( a δ − 1 ) * d δ l d ( Δ t ) from the fitted equation (e). Due to the monotonicity of δ i with respect to Δ t , i.e. d δ l d ( Δ t ) > 0 , Δ t ∈ ( 0 , C C T ) , a δ − 1 changes from positive to negative and a δ changes from greater than 1 to less than 1.
Extended Figure 7. IEEE 39bus system, Spontaneous synchronisation always occurs near the synchronisation stability boundary
Preprints 89296 i012
A three-phase short-circuit ground fault is set at the corresponding node. a), b), c), d), e) correspond to the critical behaviour of the operating points after the failure of 6bus, 12bus, 24bus, 30bus & 36bus respectively. All the results show that near the boundary, there is always a significant decrease in σ δ i . These results demonstrate a strong correlation between synchronisation stability boundary and spontaneous synchronisation.
f) and g) show that the behaviour of the operating points also show the phenomenon of synchronous barrier when faults are set at 6bus and 24bus (shaded area).

References

  1. Motter, A.E.; Myers, S.A.; Anghel, M.; Nishikawa, T. Spontaneous synchrony in power-grid networks. Nat. Phys. 2013, 9, 191–197. [Google Scholar] [CrossRef]
  2. Yu, Y.; Liu, Y.; Qin, C.; Yang, T. Theory and Method of Power System Integrated Security Region Irrelevant to Operation States: An Introduction. Engineering 2020, 6, 754–777. [Google Scholar] [CrossRef]
  3. Yang, P.; Liu, F.; Wei, W.; Wang, Z. Approaching the Transient Stability Boundary of a Power System: Theory and Applications. IEEE Trans. Autom. Sci. Eng. 2022, 1–12. [Google Scholar] [CrossRef]
  4. Molnar, F.; Nishikawa, T.; Motter, A.E. Asymmetry underlies stability in power grids. Nat. Commun. 2021, 12, 1–9. [Google Scholar] [CrossRef] [PubMed]
  5. Martínez, I.; Messina, A.R.; Vittal, V. Normal form analysis of complex system models: A structure-preserving approach. IEEE Trans. Power Syst. 2007, 22, 1908–1915. [Google Scholar] [CrossRef]
  6. Bhui, P.; Senroy, N. Real-Time Prediction and Control of Transient Stability Using Transient Energy Function. IEEE Trans. Power Syst. 2017, 32, 923–934. [Google Scholar] [CrossRef]
  7. Zhu, L.; Hill, D.J. Synchronization of Kuramoto Oscillators: A Regional Stability Framework. IEEE Trans. Automat. Contr. 2020, 65, 5070–5082. [Google Scholar] [CrossRef]
  8. Koronovskii, A.A. , Moskalenko, O.I. & Hramov, A.E. synchronization in complex networks. Tech. Phys. Lett. 2012, 38, 924–927. [Google Scholar]
  9. Casals, M.R. Knowing power grids and understanding complexity science. Int. J. Crit. Infrastructures 2015, 11, 4. [Google Scholar] [CrossRef]
  10. Gurrala, G.; Dimitrovski, A.; Pannala, S.; Simunovic, S.; Starke, M. Parareal in Time for Fast Power System Dynamic Simulations. IEEE Trans. Power Syst. 2016, 31, 1820–1830. [Google Scholar] [CrossRef]
  11. Gurrala, G. Large Multi-Machine Power System Simulations Using Multi-Stage Adomian Decomposition. IEEE Trans. Power Syst. 2017, 32, 3594–3606. [Google Scholar] [CrossRef]
  12. Wang, B.; Fang, B.; Wang, Y.; Liu, H.; Liu, Y. Power System Transient Stability Assessment Based on Big Data and the Core Vector Machine. IEEE Trans. Smart Grid 2016, 7, 2561–2570. [Google Scholar] [CrossRef]
  13. Al-Ammar, E.A.; El-Kady, M.A. Application of operating security regions in power systems. IEEE PES Transm. Distrib. Conf. Expo. Smart Solut. a Chang. World. 2010. [CrossRef]
  14. Kundur, P. Definition and classification of power system stability. IEEE Trans. Power Syst. 2004, 19, 1387–1401. [Google Scholar]
  15. Student Member, B.B.; Senior Member, G.A. On the nature of unstable equilibrium points in power systems. IEEE Trans. Power Syst. 1993, 8, 738–745. [Google Scholar]
  16. Chiang, H.D.; Wu, F.F.; Varaiya, P.P. A BCU Method for Direct Analysis of Power System Transient Stability. IEEE Trans. Power Syst. 1994, 9, 1194–1208. [Google Scholar] [CrossRef]
  17. Shubhanga, K.N.; Kulkarni, A.M. Application of Structure Preserving Energy Margin Sensitivity to Determime the Effectiveness of Shunt and Serles FACTS Devices. IEEE Power Eng. Rev. 2002, 22, 57. [Google Scholar] [CrossRef]
  18. Al Marhoon, H.H.; Leevongwat, I.; Rastgoufard, P. A fast search algorithm for Critical Clearing Time for power systems transient stability analysis. 2014 Clemson Univ. Power Syst. Conf. PSC 2014, 2014. [Google Scholar] [CrossRef]
  19. Rimorov, D.; Wang, X.; Kamwa, I.; Joos, G. An approach to constructing analytical energy function for synchronous generator models with subtransient dynamics. IEEE Trans. Power Syst. 2018, 33, 5958–5967. [Google Scholar] [CrossRef]
  20. Dörfler, F.; Chertkov, M.; Bullo, F. Synchronization in complex oscillator networks and smart grids. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 2005–2010. [Google Scholar] [CrossRef]
  21. Cuadra, L.; Salcedo-Sanz, S.; Del Ser, J.; Jiménez-Fernández, S.; Geem, Z.W. A critical review of robustness in power grids using complex networks concepts. Energies 2015, 8, 9211–9265. [Google Scholar] [CrossRef]
  22. Zhou, J. Large-Scale Power System Robust Stability Analysis Based on Value Set Approach. IEEE Trans. Power Syst. 2017, 32, 4012–4023. [Google Scholar] [CrossRef]
  23. Ajala, O.; Dominguez-Garcia, A.; Sauer, P.; Liberzon, D. A Second-Order Synchronous Machine Model for Multi-swing Stability Analysis. 51st North Am. Power Symp. NAPS 2019, 2019. [Google Scholar] [CrossRef]
  24. Karatekin, C.Z.; Uçak, C. Sensitivity analysis based on transmission line susceptances for congestion management. Electr. Power Syst. Res. 2008, 78, 1485–1493. [Google Scholar] [CrossRef]
  25. Mei, S.; Ni, Y.; Wang, G.; Wu, S. A study of self-organized criticality of power system under cascading failures based on AC-OPF with voltage stability margin. IEEE Trans. Power Syst. 2008, 23, 1719–1726. [Google Scholar]
  26. Dobson, I.; Carreras, B.; Lynch, V.; Newman, D. An initial model for complex dynamics in electric power system blackouts. Proc. Hawaii Int. Conf. Syst. Sci. 2001, 51. [Google Scholar] [CrossRef]
  27. Ding, L.; Gonzalez-Longatt, F.M.; Wall, P.; Terzija, V. Two-step spectral clustering controlled islanding algorithm. IEEE Trans. Power Syst. 2013, 28, 75–84. [Google Scholar] [CrossRef]
  28. Znidi, F.; Davarikia, H.; Rathore, H. Power Systems Transient Stability Indices: Hierarchical Clustering Based Detection of Coherent Groups Of Generators. 2021.
  29. Kuramoto, Y.; Battogtokh, D. Coexistence of Coherence and Incoherence in Nonlocally Coupled Phase Oscillators. Physics (College. Park. Md). 2002, 4, 380–385. [Google Scholar]
  30. Martens, E.A.; Thutupalli, S.; Fourrière, A.; Hallatschek, O. Chimera states in mechanical oscillator networks. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 10563–10567. [Google Scholar] [CrossRef]
  31. Panaggio, M.J.; Abrams, D.M. Chimera states: Coexistence of coherence and incoherence in networks of coupled oscillators. Nonlinearity 2015, 28, R67–R87. [Google Scholar] [CrossRef]
  32. Amirthalingam, K.M.; Ramachandran, R.P. Improvement of transient stability of power system using solid state circuit breaker. Am. J. Appl. Sci. 2013, 10, 563–569. [Google Scholar] [CrossRef]
  33. Liu, X.; Shahidehpour, M.; Cao, Y.; Li, Z.; Tian, W. Risk assessment in extreme events considering the reliability of protection systems. IEEE Trans. Smart Grid 2015, 6, 1073–1081. [Google Scholar] [CrossRef]
  34. Huang, R. Learning and Fast Adaptation for Grid Emergency Control via Deep Meta Reinforcement Learning. IEEE Trans. Power Syst. 2022, 37, 4168–4178. [Google Scholar] [CrossRef]
  35. Guo, M.; Xu, D.; Liu, L. Design of Cooperative Output Regulators for Heterogeneous Uncertain Nonlinear Multiagent Systems. IEEE Trans. Cybern. 2022, 52, 5174–5183. [Google Scholar] [CrossRef] [PubMed]
  36. Roberts, L.G.W.; Champneys, A.R.; Bell, K.R.W.; Di Bernardo, M. Analytical Approximations of Critical Clearing Time for Parametric Analysis of Power System Transient Stability. IEEE J. Emerg. Sel. Top. Circuits Syst. 2015, 5, 465–476. [Google Scholar] [CrossRef]
  37. Owusu-Mireku, R.; Chiang, H.D.; Hin, M. A Dynamic Theory-Based Method for Computing Unstable Equilibrium Points of Power Systems. IEEE Trans. Power Syst. 2020, 35, 1946–1955. [Google Scholar] [CrossRef]
  38. Sajadi, A.; Kenyon, R.W.; Hodge, B.M. Synchronization in electric power networks with inherent heterogeneity up to 100% inverter-based renewable generation. Nat. Commun. 2022, 13, 1–12. [Google Scholar]
  39. Sun, M. On-line power system inertia calculation using wide area measurements. Int. J. Electr. Power Energy Syst. 2019, 109, 325–331. [Google Scholar] [CrossRef]
  40. Zhang, Y.; Bank, J.; Muljadi, E.; Wan, Y.H.; Corbus, D. Angle instability detection in power systems with high-wind penetration using synchrophasor measurements. IEEE J. Emerg. Sel. Top. Power Electron. 2013, 1, 306–314. [Google Scholar] [CrossRef]
  41. Dörfler, F.; Bullo, F. Synchronization in complex networks of phase oscillators: A survey. Automatica 2014, 50, 1539–1564. [Google Scholar] [CrossRef]
  42. Chen, G. Searching for Best Network Topologies with Optimal Synchronizability: A Brief Review. IEEE/CAA J. Autom. Sin. 2022, 9, 573–577. [Google Scholar] [CrossRef]
  43. Li, X.; Wei, W.; Zheng, Z. Promoting synchrony of power grids by restructuring network topologies. Chaos An Interdiscip. J. Nonlinear Sci. 2023, 33, 63149. [Google Scholar] [CrossRef] [PubMed]
  44. Zhang, Y.; Motter, A.E. Symmetry-Independent Stability Analysis of Synchronization Patterns. SIAM Rev. 2020, 62, 817–836. [Google Scholar] [CrossRef]
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
bsky
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

Privacy Settings

© 2026 MDPI (Basel, Switzerland) unless otherwise stated