System-Inherent Grid Synchronous Stability Boundary and Spontaneous Synchronization

Preprint

Article

System-Inherent Grid Synchronous Stability Boundary and Spontaneous Synchronization

Altmetrics

Downloads

266

Views

Comments

Yu Yuan^*

This version is not peer-reviewed

Submitted:

31 October 2023

Posted:

02 November 2023

Read the latest preprint version here

Alerts

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:

Subject: Engineering - Electrical and Electronic Engineering

Introduction

Synchronous stability is a prerequisite for the normal operation of a power grid⁴. 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 uncertainties⁵ and nonlinearities coexist. In the stability analyses of power grids, a high level of real-time performance is required⁶. To analyze the stability of these systems, researchers have developed numerous insightful stability discriminatory methods^7–12, including finding stability regions and their boundaries^2,3,13.

Finding the stability boundary of a system is an important fundamental problem^3,14,15. The boundary, which can be used to strongly simplify stability discrimination, is the union of unstable equilibrium points (UEPs)¹⁶. When an operating point is outside the boundary, the corresponding system is unstable¹³. The synchronous stability boundary, which is a core concept of grid stability, is closely related to many issues^6,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 boundary^2,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 networks²⁰. However, the reasonableness of this assumption, which directly equates the spontaneous synchronization with the synchronization stability, has been questioned^9,21 and challenges remain in applying it to real networks⁸. 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 reported⁹. 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 synchronization¹, 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 fields^8,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 \in (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 \in (1, 2, ..., n)

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 \in [0 ， + \infty ） ， ω \in （ - \infty ， + \infty ）

. The coordinates of

ω

are transformed:

δ_{ω} = 2 \arctan (ω), δ_{ω} \in (- π, π)

, 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} = \dots = ω_{n} = 1

and

δ_{i} = 2 \arctan (ω_{i}) = \frac{π}{2}, i = 1, 2, \dots, 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}| \geq 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 \in [0 ， + \infty ） ， ω \in （ - \infty ， + \infty ）

. The coordinates of

ω

are transformed:

δ_{ω} = 2 \arctan (ω), δ_{ω} \in (- π, π)

, 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} = \dots = ω_{n} = 1

and

δ_{i} = 2 \arctan (ω_{i}) = \frac{π}{2}, i = 1, 2, \dots, 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}| \geq 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).

The essence of the connection between the meta-generators and the original generators is the permutation:

S = (\begin{matrix} {δ^{'}}_{1} & {δ^{'}}_{2} & ... & {δ^{'}}_{i} & {δ^{'}}_{i + 1} & ... & {δ^{'}}_{n} \\ δ_{1} & δ_{2} & ... & δ_{i + 1} & δ_{i} & ... & δ_{n} \end{matrix})

. 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 \in (1, 2, ..., n)

was arranged and relabelled as

（ u_{i} (t), δ_{i} (t)), i \in (1, 2, ..., n)

. This was then averaged as follows:

The mean of

（ u_{i} (t), δ_{i} (t)), i \in (1, 2, ..., n)

over

[0, T]

was found:

u_{i} = \frac{1}{T} \int_{0}^{T} u_{i} (τ) d τ

and

δ_{K L} = \frac{1}{T} \int_{0}^{T} δ_{K L} (τ) d τ = \frac{1}{T} \int_{0}^{T} |δ_{K} (τ) - δ_{L} (τ)| d τ = \frac{1}{T} \int_{0}^{T} δ_{K} (τ) - δ_{L} (τ) d τ = δ_{K} - δ_{L}

The mean of

（ u_{i} (t), δ_{i} (t)), i \in (1, 2, ..., n)

over

[T_{m}, T_{m} + d T]

was also found:

δ_{i} (t_{m}) = \frac{1}{d T} \int_{T_{m}}^{T_{m} + d T} δ_{i} (τ) d τ, U_{i} (t_{m}) = \frac{1}{d T} \int_{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 = \sqrt{u_{k}^{2} + u_{l}^{2} - 2 u_{k} u_{l} \cos δ_{k l}}

u^{'} = \sqrt{2 u_{k}^{2} (1 - \cos δ_{k l})}

Following the form of power in electricity

P = \frac{U^{2}}{R}

, the coupling power

P_{Δ u} = \frac{{|Δ 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^{'}} = \frac{{|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} \neq 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^{'}}| \leq |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

\frac{|P_{u^{'}}|}{|P_{Δ u}|} = 1 \Leftrightarrow \frac{|u^{'}|}{|Δ u|} = 1

In summary,

f (u_{K}, u_{L}, δ_{K L}) = \frac{|u^{'}|}{|Δ u|} = 1

is the system stability boundary equation. When

\frac{|u^{'}|}{|Δ u|} < 1

, the system is stable. When

\frac{|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}| \cup \frac{|u_{L}|}{|u_{K}|} = 2 \cos δ_{K L} - 1

can be found, where

|u_{K}| \geq |u_{L}| \geq 0, δ_{K L} \in [0, \frac{π}{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

\frac{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

\frac{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 \in (1, 2, ..., n)

is calculated sequentially for different

Δ t

. The standard deviation of

δ_{i}

is calculated separately (Figure 4):

σ (δ_{i}) = \sqrt{\frac{\sum_{i = 1}^{n} {(δ_{i} - μ (δ_{i}))}^{2}}{n}}, μ (δ_{i}) = \frac{\sum_{i = 1}^{n} δ_{i}}{n}

Result and Discuss

Stability boundary

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}| \cup \frac{|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}| \cup \frac{|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.

Regarding the synchronization stability, one different assumption,which appears parallel to the Kuramoto model⁸, than before ¹⁸: the system requires the synchronous power

P_{u^{'}} = \frac{{|u^{'}|}^{2}}{|Z_{K L}|}

to maintain synchronous stability, and the synchronous power is provided by the coupling power

P_{Δ u} = \frac{{|Δ 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).

(1)

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 domain². 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 reports^2,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 clear²², 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 swings²³ 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 - \frac{|u_{K}| - |u_{L}|}{2 |u_{L}|})

(2)

When

|u_{K}|, |u_{L}|

are sufficiently close²⁴, 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 unstable^25,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 \in (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 instability²⁷. Identifying coherent generators is a prerequisite for building intentional islands²⁸. 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 states^29–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.

Figure 3. Preprints 89296 i004

Figure 3a,b show that the operating points move with a uniformly variable speed before the system becomes unstable.

Contrary to intuition³²,

δ_{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 systems^33–35.

When a high degree of accuracy of the results is not needed,

a_{K L} = 0

. The following expression can be derived:

(3)

This can be used to easily and quickly check the stability margin of the system after a disturbance³⁶. By approximating

Δ T^{'}

as the CCT^19,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 estimated³⁷.

(4)

In summary, the CCT and UEP can be calculated using only information about the rotation rate¹², 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 system^38–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).

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}

Δ t

(shaded area). c). In the

u_{1} - u_{2}

As the stability boundary is approached, the perturbed trajectories of the operating points become interesting. The system self-organized and moved toward synchronous evolution¹.

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 synchronization^7,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 network^42–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 (

\frac{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 \in (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

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

Extended Figure 3. IEEE 39-BUS system, self-organizing behavior of generators near the boundary

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

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

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