Analytical Analysis of Power Network Stability: Necessary and Sufficient Conditions

Abstract

The investigation of the synchronization of Kuramoto oscillators is a crucial applied model for studying harmonization in oscillating phenomena across physical, biological, and engineering networks. This chapter builds on previous studies by exploring the synchronization of Kuramoto oscillators while also conforming to more realistic models. Using the LaSalle Invariance Principle and contraction property, we introduce the necessary and sufficient conditions for frequency synchronization and phase cohesiveness. The novelty of this chapter’s contents lies in three key areas: First, we consider a heterogeneous second-order model with non-uniformity in coupling topology. Second, we apply a non-zero and non-uniform phase shift in coupling function. Third, we introduce a new Lyapunov-based stability analysis technique. Our findings demonstrate that heterogeneity in the network and the phase shift in the coupling function are key factors in network synchronization. We present the synchronization conditions based on network graph-theoretical characteristics and the oscillators’ parameters. Analysis of the results reveals that an increase in the phase shift and heterogeneity of oscillators will complicate the synchronization conditions. Numerical simulations confirm the validity of our theoretical results. One of the main applications of this study is the development of stability conditions for smart grids with Lossy-Power Network.

1. Introduction

Synchronization in interconnected networks of heterogeneous or homogeneous oscillators is a pervasive phenomenon in various fields such as biology, physics, chemistry, engineering, and social networks. The electrical network, as a prime example, is one of the most important and complex oscillatory networks constructed by engineering science. In these networks, two main factors, namely oscillator dynamics and network communication graph, lead to different dynamic behaviors. Therefore, the overall network response in different regimes can be controlled by the following two approaches:

• Designing and tuning an appropriate control system on the controllable parameters of the oscillatory agents.

• Adjusting the network communication graph based on the parameters of the oscillatory agents.

The future of the electrical network is moving towards increased use of renewable energy sources in the form of distributed generation. On the other hand, the penetration of renewable energy sources leads to an increase in imbalances between electricity production and consumption, along with increased complexity in the dynamic behavior and greater heterogeneity of network factors. Therefore, providing a suitable mathematical model for the studied network [1], understanding the network stability status, and presenting control theories aimed at increasing its stability and robustness against various regimes are among the main concerns of control engineers for the network.

Employing a suitable mathematical model for the studied network is a prerequisite for implementing each of the above approaches [1]. One of the most widely used models in oscillatory networks is the model proposed by Kuramoto. This model describes the dynamic behavior of oscillators whose interactions between agents are proportional to the sine of the phase difference between them.

The electrical network is one of the most important and complex oscillatory networks constructed by engineering science. Due to the oscillatory nature of the dynamic equations of this network, there is a close relationship between its mathematical model and the non-uniform Kuramoto model. Based on this, and considering the limitations of classical methods in studying the stability of electrical networks based on analytical theories, in this thesis, by mapping the oscillation equations of the electrical network to the non-uniform Kuramoto model and presenting analytical theories, a stability condition in the form of an explicit and precise relationship based on the parameters and topology of the network graph is obtained. Furthermore, based on the results obtained from the proposed theories, suggestions for improving network stability are presented.

1.1 Notations

Consider the vector $x={\left(x_1,\dots ,x_n\right)}^T$. The values ${}_{\min }$ and $x_{\max }$ are respectively the smallest and largest elements of this vector. The 2-norm of vector $x$ is represented by ${‖x‖}_2$, and $diag\left(x\right)\in {\mathbb{R}}^{n\times n}$ represents the associated diagonal matrix. The notations $0_{\mathrm{n}}$ and $1_{\mathrm{n}}$ are defined as column vectors with zeros and ones in all entries, respectively. The sign $\circ$ represents the entrywise product (Hadamard Product) of the two vectors, and for a complex number $S$, the terms $\mathfrak{R}\left(S\right)$ and $\mathfrak{I}\left(S\right)$ are the real and imaginary parts of $S$.

A triple set $G=\left(v,\epsilon ,A\right)$ is a weighted directed graph where $A\in {\mathbb{R}}^{n\times n}$ is the adjacency matrix, $v=\left\{1,\dots ,n\right\}$ is the set of vertices, and $\epsilon \subset v\times v$ is the set of edges. In the adjacency matrix, $a_{\mathit{ij}}>0$ if there is a directed edge from vertex $i$ to $j$, otherwise $a_{\mathit{ij}}=0$. The Laplacian matrix of $G$, denoted by $L$, is defined as $L:=diag\left({\sum }_{j=1}^n{a}_{\mathit{ij}}\right)-A$. For the incidence matrix $B=H^T\in {\mathbb{R}}^{n\times \left|\epsilon \right|}$, the entry $B_{\mathit{kl}}=-1$ if edge $l$ is directed to vertex $K$, and $B_{\mathit{kl}}=1$ if edge $l$ is directed from vertex $K$, otherwise $B_{\mathit{kl}}=0$. Notation $ker\left(H\right)$ is relevant to the null space of matrix $H$. The complete graph is defined as a graph where all vertices are connected. If $n$ is equal to the number of vertices of the complete graph, the number of edges is given by $n\left(n-1\right)/2$. If $G$ is a connected graph, then $ker\left(H\right)=ker\left(L\right)=\text{span}\left({\mathbf{1}}_n\right)$, all non-zero eigenvalues of the Laplacian matrix are positive, and ${\lambda }_2\left(L\right)$, the second-smallest eigenvalue of $G$, is named the algebraic connectivity. ${\lambda }_{\min }\left(L\right)$ and ${\lambda }_{\max }\left(L\right)$ correspond to the minimum and maximum eigenvalues of matrix $L$, respectively. Moreover, ${\lambda }_2\left(L\right)=n$ in a complete and uniformly weighted graph ($a_{\mathit{ij}}=1$ for all $i\ne j$).

The torus${\mathbb{T}}^1\mathit{=}\left(-\pi ,\pi ]\right.$ is a set, and angle $\theta$ is a point of the torus where $\theta \in {\mathbb{T}}^1$, and the $\mathit{arc}$ corresponds to the subset of ${\mathbb{T}}^1$. The distance $\left|{\theta }_i-{\theta }_j\right|$ is defined as the minimum length between two angles ${\theta }_i\in {\mathbb{T}}^1$ and ${\theta }_j\in {\mathbb{T}}^1$. For $\rho \in \left[0,\pi /2\right]$, the $∆\left(\rho \right)\subset {\mathbb{T}}^1$ is the set of surrounded angles $\left({\theta }_1,\dots ,{\theta }_n\right)\in {\mathbb{T}}^n$, where the $\mathit{arc}$ with length $\rho$ comprises all the angles ${\theta }_1,\dots ,{\theta }_n$. Accordingly, for every angle of the set $∆\left(\rho \right)$, the inequality ${max}_{\left\{i,j\right\}\in \epsilon }\left|{\theta }_i-{\theta }_j\right|<\rho$ is satisfied. For $\in {\mathbb{R}}^n$, we define the vector-valued function $\mathbf{sin}\left(x\right)=\left(sin\left(x_1\right),\dots ,sin\left(x_n\right)\right)$ and the sinc function $sinc:\mathbb{R}\to \mathbb{R}$ by $sinc\left(x\right)=sin\left(x\right)/x$.

1.2 Classical methods for stability analysis

The classical methods used to determine the stability of a network in the early years of its development were based on available computational tools, such as the ability to maintain network stability during numerical simulations of worst-case scenarios in the design or based on the direct Lyapunov method [2, 3, 4, 5, 6]. However, gradually, with the application of these methods in real networks, some inefficiencies in the performance of classical methods became apparent such as the following items [7, 8]:

1. The need for complex and time-consuming numerical computations in numerical simulations.

2. Inefficiency of the direct Lyapunov method in networks with transmission losses. Considering the energy function as the equation $E={\sum }_{i=1}^n\left(M_i{\dot{\theta }}_i^2/2\right)-{\sum }_{\left\{i,j\right\}\in \epsilon }^n{P}_{\mathit{ij}}cos{\theta }_{\mathit{ij}}-{\sum }_{i=1}^n{\omega }_i{\theta }_i$ (where $\theta ={\left({\theta }_1,\dots ,{\theta }_n\right)}^T$and ${\theta }_{\mathit{ij}}={\theta }_i-{\theta }_j$), and the oscillation equations for lossless electrical networks as $M_i{\ddot{\theta }}_i+D_i{\dot{\theta }}_i={\omega }_i-{\sum }_{j=1}^n{P}_{\mathit{ij}}sin\left({\theta }_i-{\theta }_j\right)\forall i\in \left\{1,\dots ,n\right\}$, the derivative of the energy function is obtained as $\dot{E}=-{\sum }_{k=1}^n{D}_i{\dot{\theta }}_i^2\le 0$. In this case, the network is locally stable [9]. This function can be extended to networks with transmission losses as $E_{\mathit{\text{loss}}}={\sum }_{i=1}^n\left(M_i{\dot{\theta }}_i^2/2\right)-{\sum }_{\left\{i,j\right\}\in \epsilon }^n{P}_{\mathit{ij}}cos\left({\theta }_{\mathit{ij}}+{\varphi }_{\mathit{ij}}\right)-{\sum }_{i=1}^n{\omega }_i{\theta }_i$, which is obtained for $\dot{E}=-{\sum }_{i=1}^n{D}_i{\dot{\theta }}_i^2-2{\sum }_{\left\{i,j\right\}\in \epsilon }^n{P}_{\mathit{ij}}{\dot{\theta }}_{j}cos{\theta }_{\mathit{ij}}cos{\varphi }_{\mathit{ij}}$. As observed, this function is not always decreasing and cannot guarantee network stability.

3. The inability to employ network controllers and topological information for feedback control transmission in stability studies, and the excessive conservatism in these studies.

1.3 Transient stability study through Kuramoto model

In recent years, the increasing random disturbances in the network due to the growing demand for electricity consumption and the penetration of renewable energy sources have highlighted the limitations of classical methods for analyzing network stability. Alongside the consequences of instability in electrical networks and the need to describe network stability through explicit mathematical relationships based on network parameters, the main motivation for using structure-oriented analytical methods as alternatives to classical analytical methods has been the development of measurement and computational technologies that provide a suitable platform for implementing their results in the network [10]. In this regard, a review of the history of studies shows that the oscillatory nature of electrical networks and the similarity between the model of these networks and the Kuramoto oscillatory model have led to a strong tendency to use this model for analyzing the synchrony of electrical networks.

In [11] Jadbabaie demonstrated that the region $\theta \in ∆\left(\pi /2\right)$ for $P>0$ is an invariant set for a first-order oscillatory network with a complete graph and equal natural frequencies, ${\dot{\theta }}_i={\omega }_i+\left(P/n\right){\sum }_{j=1}^{n}sin\left({\theta }_j-{\theta }_i\right)$, state trajectories of each agents simultaneously converge at rate $\left(2P/\mathit{\pi n}\right){\lambda }_2\left(L\right)$. Additionally, for a network with agents having asynchronous natural frequencies and a complete graph, the phase coherence condition of agents is obtained as $P\ge P_L:=2\sqrt{n}{‖\omega ‖}_2/{\lambda }_2\left(L\right)$ and $P\ge P_L:=n{\pi }^2{\lambda }_{\max }\left(L\right){‖\omega ‖}_2/4{\lambda }_2{\left(L\right)}^2$, respectively, in $\theta \in ∆\left(\pi /2\right)$. Further studies by Chopra and Spong obtained the necessary condition for synchrony in terms of $P\ge P_{\mathit{ess}}=n\left({\omega }_{\max }-{\omega }_{\min }\right)/E_{\max }$ (with $E_{\max }$ being a computational parameter of the network model) by studying the dynamics of rotational frequency difference of oscillators [12]. This condition becomes the sufficient condition for synchrony under extreme conditions $P_{\mathit{ess}}=\left({\omega }_{\max }-{\omega }_{\min }\right)/2$. Furthermore, in this reference, the phase coherence condition of agents in $\theta \in ∆\left(\rho \right)\forall \rho \in \pi /2$ is obtained as $P\ge P_{\mathit{suf}}=n\left({\omega }_{\max }-{\omega }_{\min }\right)/cos\left(\pi /2-\rho \right)$. To complement the research [12], Choi et al. used Dini derivatives to derive necessary conditions for stability in phase-locked states, taking into account the initial conditions of the model in the form of $P\ge \left({\omega }_{\max }-{\omega }_{\min }\right)/sin\left({\theta }_{\max }-{\theta }_{\min }\right)$. In [13], they presented these conditions as a function of the initial phase of the model. In [13] sufficient conditions for frequency synchronization were also obtained. In [14], Ha et al. built upon the results from [13] and investigated the effects of dynamic frequency natural factors and phase shifts on the synchronization of oscillators in a limited set of oscillators. They derived stability conditions as functions of the initial phase of the oscillators, the amount of phase shift, and the intensity of network communications. It should be noted that the impact of losses in the dynamic model of transmission lines appears as phase delay. Previously, studies on the stability of oscillatory networks with phase shifts have been evaluated in [15, 16, 17, 18, 19, 20, 21]. In [15], by introducing the parameter link frustration and minimizing the values of dynamic states, the dynamic behavior of networks with different topologies was investigated. In [16], simulations showed how changing the value of phase shift can affect agent asynchrony. However, it is important to note that the results obtained from these studies rely on numerical approaches.

In [21, 22], as one of the first steps in utilizing oscillatory models in the stability analysis of electrical networks, the threshold level of disturbances for each node that leads to loss of synchrony was calculated through numerical analyses, and the relationship between this threshold level and the topological characteristics of the network graph and the time of loss of synchrony was determined. The results obtained from these analyses indicate a clear relationship between the dynamic and topological parameters of the network, such that an increase in the degree of graph nodes increases the minimum resolvable disturbances in the network. These results can be used in the design or upgrading of networks with the aim of improving network stability. Additionally, Dörfler and his colleagues in [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32], analyzed the stability of electrical networks using oscillatory models. Building upon the results of studies [33, 34, 35, 36], Dörfler and Bullo in [26, 27] derived sufficient conditions for transient stability in electrical networks based on the topological characteristics of the network graph, equipment parameter values, and initial conditions, by considering the assumption of a small ratio of inertia parameter to damping coefficient of network components, $M_{\max }/D_{\min }$, and employing techniques based on Singular Perturbation and Lyapunov stability.

In [28], Dörfler and his colleagues presented a more explicit condition for achieving phase coherence among agents by analyzing the first-order oscillatory model on synchronous manifolds. Based on the results obtained for the first-order model, it was shown that by satisfying the condition ${‖B^T{L}^{†}\omega ‖}_{\infty }\le sin\left(\rho \right)$, a non-trivial region for the phase difference of the agents $\mid {\theta }_i-{\theta }_j\mid \le \rho \le \pi /2$ can be established.

Given the importance of frequency control and synchronization in microgrids, the adjustment of network control systems is of great importance. In [25, 26, 28, 31, 32, 37, 38, 39, 40, 41, 42, 43], stability studies of microgrids have been conducted using oscillatory models. In [29] and [30], the frequency stability of an inverter-based microgrid with Power-Frequency Droop Controllers is investigated assuming $M_i=0$, and the necessary and sufficient conditions for synchronization in the microgrid are obtained. In [44], the impact of network decentralization on stability is studied numerically in a power network with distributed energy resources. In [45, 46], using model [47], the necessary and sufficient condition for transient stability in an electric network with controllable injection power to oscillators is presented as ${\dot{\omega }}_i=K_i{\sum }_{j=1}^n{b}_{\mathit{ij}}\left({\dot{\theta }}_j-{\dot{\theta }}_i\right)$. In [38], by employing theories from [27, 30, 48, 49, 50], and designing a distributed adaptive droop controller to enhance the algebraic connectivity of the network, the necessary and sufficient condition for stability of low-inertia microgrids is developed. In [51], by using the Kron reduction process and adjusting parameters including generator droop and network damping coefficient, the eigenvalues of the system are controlled for network stability. In [39], a distributed proportional-integral controller is employed to analyze the synchronization and phase tracking problem.

In [52], using Gronwall’s inequality, the convergence rate is expressed as a function of inertia, coupling strength, and natural frequency for phase-shift-free communications. In [53], the same studies are conducted for a second-order model with phase shift in the communication function. However, the results of these studies are valid for networks with homogeneous agents, uniform communications, and equal phase shifts. Meanwhile, Choi et al. calculate the synchronization conditions of this model in [54] using energy functions as a function of the initial topology, distribution of natural frequencies, and agent inertia for networks with homogeneous agents. In [55], the same studies are conducted for networks with heterogeneous oscillatory agents under certain assumptions. Thanh Long et al. extend their previous stability analysis methods in [56, 57, 58, 59] and classical stability analysis methods in [60, 61, 62, 63, 64, 65] for electric networks with transmission line losses in [8]. They develop a Lyapunov function for networks with lossy transmission lines using an iterative approach based on initial conditions. Bo Li et al. in [66] present an optimization algorithm for the coupling strength and input power to agents in order to increase the value of ${\lambda }_2\left(L\left({\theta }^{\ast }\right)\right)$ in a linear model of a power network. In [67], Grzybowsk et al. provide an estimation of the synchronization basin region in a power network based on the second-order Kuramoto model, considering assumptions such as equal damping of oscillatory agents and neglecting transmission line resistances.

In [68, 69], an estimation of the stable synchronization region of an electric network is provided using the second-order Kuramoto model, based on the Lojasiewicz and Gronwall’s inequalities and the assumption $m_i/d_i=m_j/d_j$ for $i\ne j$. In [70], the same studies are conducted without assumptions and limitations.

The authors in [71] continue their studies from [10, 27, 32, 72], and their previous work in [40, 41], by introducing a Lyapunov function and using the second-order Kuramoto model to obtain a sufficient condition for synchronization in an electric network with lossless transmission lines. In [72], local stability analysis for the electric network model introduced in [47] is also conducted. In [73], Yufeng et al. calculate the critical outage time and return the network to a stable state using the estimated synchronization region for network frequency synchronization in reference to [68] through numerical analysis. In [74, 75, 76, 77, 78, 79], analysis of multi-zone or cluster oscillatory networks is conducted.

Based on the above, the study of stability in oscillatory network has shown that by performing the following transformations in the used model, the results of these studies have become closer to reality over the past two decades:

• First-order models to second-order models;

• Homogeneous agents to heterogeneous agents;

• Uniform coupling graph to non-uniform coupling graph;

• Phase-shift-free coupling to the coupling topology with phase shift.

1.4 Objectives and scope of the chapter

As mentioned earlier, it is not possible to determine the stability condition of a network explicitly based on the characteristics of the agents and the Laplacian matrix of the network using classical stability analysis methods. Therefore, one of the open issues in power system conferences and forums is to determine explicit and accurate conditions based on the network characteristics that guarantee transient stability.

Reviewing the conducted activities shows that the developed theories in recent studies have usually involved assumptions about the oscillatory network agents or specific conditions in the network graph. As a result, the results of each reference have been developed based on previous studies, and the assumptions have gradually decreased. Generally, these assumptions can be categorized as follows:

• Considering homogeneous parameters of network agents such as damping coefficient and inertia in [13, 14, 53].

• Assuming uniform coupling strengths between agents in [11, 12, 13, 14, 53].

• Assuming overdamping conditions in the network and transforming the oscillation equations from second-order models to first-order models using singularity perturbation techniques in [26, 27, 28, 30, 45, 46].

• Local stability studies based on linearized network models in [51, 72, 79, 80].

• Assuming equal inertia-to-damping ratio for network agents in [55, 68, 69].

• Developing results based on analyzing the numerical simulation results in [16, 22, 44, 66, 73, 74].

• Considering small inertia of oscillatory network agents in [25, 26, 27, 28, 31, 32, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46] stability analysis of microgrids that cannot be applied to electric networks with conventional power plants.

• Analyzing lossless transmission line network models except for references [8, 14, 16, 17, 27, 52, 53].

• Analyzing network stability without considering the effect of control signals in the network except for references [25, 29, 30, 39, 45, 46, 51].

In more comprehensive models, such as [8, 52, 53, 67, 68, 69, 70, 71, 73], assumptions and constraints have been introduced. Some examples include:

• Calculation of the Lyapunov function based on iterative algorithms [8].

• Studies on lossless networks [52].

• Assuming equal inertia-to-damping ratio for network agents in [55].

• Limited number of oscillators in [68, 69].

• Employing the assumption of a complete graph in a network with an incomplete graph [71].

Therefore, it can be claimed that many of the theories presented in recent references have provided valuable results in complementing previous studies. In the following sections, the necessary and sufficient conditions for the stability of an electrical network with significant losses, heterogeneous oscillatory agents, and non-uniform communications are presented based on explicit and clear relationships between the characteristics of the network graph, initial conditions, and properties of oscillatory agents.

2. Fundamentals of power network transmission

2.1 Energy balance in generators

In the electric network oscillating model, the dynamics of the network’s generators are usually described using state variables such as phase angle ${\delta }_i$ and rotor angular velocity $d{\delta }_i/\mathit{dt}$. The value of $\mathrm{\Omega }$ corresponds to the nominal frequency of the network and is equal to $\mathrm{\Omega }=2\pi \times 50\mathit{Hz}$ Hz or $\mathrm{\Omega }=2\pi \times 60\mathit{Hz}$ Hz. Therefore, the phase angle of the generators will be:

$${\delta }_i\left(t\right)=\mathrm{\Omega }t+{\Theta }_i\left(t\right)$$E1

Where ${\Theta }_i$ represents the phase difference between the generator and the reference generator $\mathrm{\Omega }t$. During the rotation of the generators in the electric network, the amount of dissipated power is obtained by the expression $P_{\mathit{\text{diss}},i}=K_{D,i}{\dot{\delta }}_i^2$, where $K_{D,i}$ is equal to the overall losses coefficient, friction coefficient or other power loss factors. The amount of kinetic energy and mechanical power in the electric generators are equal to $E_{\mathit{kin},i}=I_i{\dot{\delta }}_i^2/2$ and $P_{\mathit{acc},i}={\mathit{dE}}_{\mathit{kin},i}/\mathit{dt}$ respectively, where $I_i$ is equivalent to the moment of inertia of the $i$-th generator. If there is power flow between two nodes $i$ and $j$, then the transfer power between node $i$ and $j$ is denoted by $P_{\mathit{\text{trans}},\mathit{ij}}$. In order to achieve power balance in a generator, it is necessary for the injected power $P_{\mathit{\text{source}},i}$ to the generator to be equal to the sum of kinetic power, dissipated power, and transferred power from the generator, $P_{\mathit{\text{trans}},i}={\sum }_{j=1}^N{P}_{\mathit{\text{trans}},\mathit{ij}}$. In other words, $P_{\mathit{\text{source}},i}=P_{\mathit{\text{diss}},i}+P_{\mathit{acc},i}+{\sum }_{j=1}^N{P}_{\mathit{\text{trans}},\mathit{ij}}$. Substituting the power balance equation in terms of generator states, we obtain:

$$P_{\mathit{\text{source}},i}=K_{D,i}{\dot{\delta }}_i^2+I_i{\ddot{\Theta }}_i{\dot{\delta }}_i+P_{\mathit{\text{trans}},i}$$E2

Utilizing ${\delta }_i\left(t\right)=\mathrm{\Omega }t+{\Theta }_i\left(t\right)$ in the prior equation results:

$$\begin{array}{rcl}{P}_{\mathit{\text{source,}}i}& =& K_{D,i}{\left(\mathrm{\Omega }+{\dot{\Theta }}_i\right)}^2+I_i{\ddot{\Theta }}_i\left(\mathrm{\Omega }+{\dot{\Theta }}_i\right)+P_{\mathit{\text{trans,}}i}\\ & =& K_{D,i}{\mathrm{\Omega }}^2+2K_{D,i}\mathrm{\Omega }{\dot{\Theta }}_i+K_{D,i}{\dot{\Theta }}_i^2+I_i\mathrm{\Omega }{\ddot{\Theta }}_i+I_i{\ddot{\Theta }}_i{\dot{\Theta }}_i+P_{\mathit{\text{trans,}}i}\end{array}$$E3

Where, $2K_{D,i}\mathrm{\Omega }$ is equal to the damping coefficient. Assuming that the impact of disturbances on the rotational frequency of network generators compared to the nominal frequency of the network ${\dot{\Theta }}_i\ll \mathrm{\Omega }$, then the terms $I_i{\ddot{\Theta }}_i{\dot{\Theta }}_i$ and $K_{D,i}{\dot{\Theta }}_i^2$ can be neglected compared to $I_i\mathrm{\Omega }{\ddot{\Theta }}_i$ and $K_{D,i}{\mathrm{\Omega }}^2$[81]. The power balance equation in network generators will be as follows:

$$I_i\mathrm{\Omega }{\ddot{\Theta }}_i=P_{\mathit{\text{source}},i}-K_{D,i}{\mathrm{\Omega }}^2-2K_{D,i}\mathrm{\Omega }{\dot{\theta }}_i-P_{\mathit{\text{trans}},i}$$E4

To simplify, we define $M_i:=I_i\mathrm{\Omega }$ as the moment of inertia, $P_{m,i}:=P_{\mathit{\text{source}},i}-K_{D,i}{\mathrm{\Omega }}^2$ as the net power provided to the $i$-th generator, and $D_i:=2K_{D,i}\mathrm{\Omega }$ is defined as the damping coefficient for the $i$-th generator. Thus the power balance network generators can be simplified to the following:

$$M_i{\ddot{\Theta }}_i=P_{m,i}-D_i{\dot{\Theta }}_i-P_{\mathit{\text{trans}},i}$$E5

2.2 Electric network model reduction

To calculate $P_{\mathit{\text{trans}},i}$ as Eq. (5), the real part of matrix$\mathfrak{R}\left(V\circ I^{\ast }\right)$, the vectors of node voltages$V$, and the injected current to the nodes $I$ are necessary. The value of the vector $I$ is obtained the reference node form as $I=Y_{\mathit{net}}V$, where $Y_{\mathit{net}}$ is the admittance matrix of the network. With the assumption of static loads in the electric network, the total of $N$ network nodes can be divided into dynamic nodes $N_G$ and static nodes $N_G$. So the network admittance matrix $Y_{\mathit{net}}$ can be divided into four separate sections as follows:

$$Y_{\mathit{net}}=\left[\begin{array}{cc}{Y}_{\mathit{GG}}& Y_{\mathit{LG}}\\ Y_{\mathit{LG}}& Y_{\mathit{LL}}\end{array}\right]$$E6

The injected current values to each of the network nodes will be as follows:

$$\left[\begin{array}{c}{I}_G\\ I_L\end{array}\right]=Y_{\mathit{net}}\left[\begin{array}{c}{V}_G\\ V_L\end{array}\right]\to \begin{array}{c}{I}_G=Y_{\mathit{GG}}{V}_G+Y_{\mathit{LG}}{V}_L\\ I_L=Y_{\mathit{LG}}{V}_G+Y_{\mathit{LL}}{V}_L\end{array}$$E7

Where $V_G={\left[V_1\dots V_{N_G}\right]}^T$ represents the node voltages of the network generators, $V_L={\left[\begin{array}{ccc}{V}_{N_G+1}& \dots & V_N\end{array}\right]}^T$ is the node voltages of the network loads, $I_G={\left[\begin{array}{ccc}{I}_1& \dots & I_{N_G}\end{array}\right]}^T$ is the injected currents to the network generators, and $I_L={\left[\begin{array}{ccc}{I}_{N_G+1}& \dots & I_N\end{array}\right]}^T$ represents the injected currents to the network loads. In the following, the assumption of zero resistance for transmission lines has been disregarded, and an attempt has been made to use a more realistic model of the network in stability analyses.

As mentioned before, assuming the static nature of loads compared to the dynamics of network generators, the network loads can be modeled as admittances shunted to the reference node in an $N_L$ node model. It is worth noting that since transient stability studies are conducted in this chapter and considering the time frame for transient stability studies, the assumption of static loads is reasonable. In this case, by applying the Kron reduction method, we can obtain the voltage $V_L$ in terms of $V_G$ and substitute it into equation $I_G$.

$$V_L=Y_{\mathit{LL}}^{-1}{I}_L-Y_{\mathit{LL}}^{-1}{Y}_{\mathit{LG}}{V}_G$$E8

$$I_G=Y_{\mathit{GG}}{V}_G-Y_{\mathit{LG}}\left(Y_{\mathit{LL}}^{-1}{Y}_{\mathit{LG}}{V}_G\right)+Y_{\mathit{LG}}{Y}_{\mathit{LL}}^{-1}{I}_L=\left(Y_{\mathit{GG}}-Y_{\mathit{LG}}{Y}_{\mathit{LL}}^{-1}{Y}_{\mathit{LG}}\right)V_G+Y_{\mathit{LG}}{Y}_{\mathit{LL}}^{-1}{I}_L$$E9

So the reduced admittance matrix of the network will be $Y_{\mathit{red},g}=Y_{\mathit{GG}}-Y_{\mathit{LG}}{Y}_{\mathit{LL}}^{-1}{Y}_{\mathit{LG}}$. By extracting the shunt admittance due to static loads from the matrix $Y_{\mathit{red},g}$ and inputting it into the generator nodes, Eq. (7) will be $I_G=Y_{\mathit{red}}{V}_G$. In this regard, the net injected power by the generators into the network will be given by the following:

$$S_{\mathit{\text{trans}}}=V_G\circ I_G^{\ast }=V_G\circ \left(Y_{\mathit{red}}^{\ast }{V}_G^{\ast }\right)=V_G\circ Y_{\mathit{red}}^{\ast }{V}_G^{\ast }+S_{L2G}$$E10

Where $S_{L2G}=V_G\circ {\left(Y_{\mathit{LG}}{Y}_{\mathit{LL}}^{-1}{I}_L\right)}^{\ast }$ represents the apparent power effect of network loads on network generator nodes. By substituting the voltage phasors of the generator nodes and the reduced admittance matrix, we have:

$$S_{\mathit{\text{trans}},i}=S_{L2G,i}+\sum _{j=1}^{N_G}\left|V_i\right|\left|V_j\right|\left|Y_{\mathit{red},\mathit{ij}}\right|e^{i\left({\Theta }_i-{\Theta }_j-{\phi }_{\mathit{ij}}\right)}$$E11

Where ${\Theta }_i$ and ${\Theta }_j$ respectively represent the voltage angles of generator $i$ and $j$, and ${\phi }_{\mathit{ij}}$ represents the phase angle of the reduced admittance matrix between two nodes $i,j\in \left\{1,\dots ,N_G\right\}$, which is obtained from the equation ${\phi }_{\mathit{ij}}=\text{atan}\left(\mathfrak{I}\left(Y_{\mathit{red},\mathit{ij}}\right)/\mathfrak{R}\left(Y_{\mathit{red},\mathit{ij}}\right)\right)$. Taking into account the power effect caused by network loads on network generator nodes as $P_{L2G,i}=\mathfrak{R}\left(S_{L2G,i}\right)$, the active transfer power from generator $i$ will be as follows:

$$P_{\mathit{\text{trans}},i}=P_{L2G,i}+\sum _{j=1}^{N_G}\left|V_i\right|\left|V_j\right|\left|Y_{\mathit{red},\mathit{ij}}\right|cos\left({\Theta }_i-{\Theta }_j-{\phi }_{\mathit{ij}}\right)$$E12

Employing the trigonometric identity $sin\left(x\right)=cos\left(x-\pi /2\right)$ and defining the parameter ${\varphi }_{\mathit{ij}}=-{\phi }_{\mathit{ij}}+\pi /2$, the equation for active transfer power from generator $i$ will be transformed to:

$$P_{\mathit{\text{trans}},i}=P_{L2G,i}+\sum _{j=1}^{N_G}\left|V_i\right|\left|V_j\right|\left|Y_{\mathit{red},\mathit{ij}}\right|sin\left({\Theta }_i-{\Theta }_j+{\varphi }_{\mathit{ij}}\right)$$E13

Moreover, applying the identities $sin\left(\alpha \pm \beta \right)=sin\left(\alpha \right)cos\left(\beta \right)\pm cos\left(\alpha \right)sin\left(\beta \right)$ and $cos\left(\alpha \pm \beta \right)=cos\left(\alpha \right)cos\left(\beta \right)\mp sin\left(\alpha \right)sin\left(\beta \right)$), the relationship between the angle ${\varphi }_{\mathit{ij}}$ and the reduced admittance matrix will be as follows:

$$tan\left({\varphi }_{\mathit{ij}}\right)=tan\left(-{\phi }_{\mathit{ij}}+\frac{\pi }{2}\right)=\frac{sin\left(-{\phi }_{\mathit{ij}}+\frac{\pi }{2}\right)}{cos\left(-{\phi }_{\mathit{ij}}+\frac{\pi }{2}\right)}=\frac{sin\left(\frac{\pi }{2}\right)cos\left(-{\phi }_{\mathit{ij}}\right)}{-sin\left(\frac{\pi }{2}\right)sin\left(-{\phi }_{\mathit{ij}}\right)}=\frac{\mathfrak{R}\left(Y_{\mathit{red},\mathit{ij}}\right)}{\mathfrak{I}\left(Y_{\mathit{red},\mathit{ij}}\right)}$$E14

Lemma 1: Applying the Kron reduction method in an electrical network, the off-diagonal and diagonal elements of the reduced admittance matrix $Y_{\mathit{red}}$ will be located in the regions $\left(\pi /2,\pi \right)$ and $\left(-\pi /2,0\right)$ (respectively)

Proof: In an electrical network, the off-diagonal and diagonal elements of the $Y_{\mathit{net}}$ matrix, representing non-diagonal and diagonal elements, are located in the regions $\left(\pi /2,\pi \right)$ and $\left(-\pi /2,0\right)$ respectively. Utilizing the Kron reduction method and eliminating the node corresponding to the load with index $p\in \left\{N_G+1,\dots ,N\right\}$ in the $Y_{\mathit{net}}$ matrix, the reduced element $Y_{\mathit{red},\mathit{jk}}$ can be obtained from the following equation:

$$Y_{\mathit{red},\mathit{jk}}=Y_{\mathit{net},\mathit{jk}}-\frac{Y_{\mathit{net},\mathit{jp}}{Y}_{\mathit{net},\mathit{pk}}}{Y_{\mathit{net},\mathit{pp}}}$$E15

The elements of the network admittance matrix are represented in phasor form as $Y_{\mathit{net},\mathit{jk}}=\mid Y_{\mathit{net},\mathit{jk}}\mid \measuredangle {\Theta }_{\mathit{jk}}$, $Y_{\mathit{net},\mathit{jp}}=\mid Y_{\mathit{net},\mathit{jp}}\mid \measuredangle {\Theta }_{\mathit{jp}}$, $Y_{\mathit{net},\mathit{pk}}=\mid Y_{\mathit{net},\mathit{pk}}\mid \measuredangle {\Theta }_{\mathit{pk}}$, and $Y_{\mathit{net},\mathit{pp}}=\mid Y_{\mathit{net},\mathit{pp}}\mid \measuredangle {\Theta }_{\mathit{pp}}$, where their phases lie in the regions ${\Theta }_{\mathit{pk}},{\Theta }_{\mathit{jp}},{\Theta }_{\mathit{jk}}\in \left(\pi /2,\pi \right)$ and ${\Theta }_{\mathit{pp}}\in \left(-\pi /2,0\right)$. Based on these angles, the phase of the element $\mathit{jk}$-th in the reduced matrix $Y_{\mathit{red}}$ can be expressed as a phasor, $Y_{\mathit{red},\mathit{ij}}=\mid Y_{\mathit{red},\mathit{ij}}\mid \measuredangle {\phi }_{\mathit{red},\mathit{jk}}\forall {\phi }_{\mathit{red},\mathit{jk}}\in \left(\pi /2,\pi \right)$. In other words, $\mathfrak{R}\left(Y_{\mathit{red},\mathit{ij}}\right)$ and $\mathfrak{I}\left(Y_{\mathit{red},\mathit{ij}}\right)$ will have negative and positive values, respectively. Therefore, the phases of the off-diagonal and diagonal elements of the $Y_{\mathit{red}}$ matrix will lie in the regions $\left(\pi /2,\pi \right)$ and $\left(-\pi /2,0\right)$, respectively. This concludes the proof of this Lemma.

As mentioned in Lemma 1, in the reduced network admittance matrix of the electrical network, without considering losses in transmission lines, the angle ${\phi }_{\mathit{ij}}=\text{atan}\left(\mathfrak{I}\left(Y_{\mathit{red},\mathit{ij}}\right)/\mathfrak{R}\left(Y_{\mathit{red},\mathit{ij}}\right)\right)$ will fall within the range ${\phi }_{\mathit{ij}}\in \left(\pi /2,\pi \right)$. Consequently, the angle ${\varphi }_{\mathit{ij}}=-{\phi }_{\mathit{ij}}+\pi /2$ will also lie within the range ${\varphi }_{\mathit{ij}}\in \left(0,\pi /2\right)$. By substituting $P_{\mathit{\text{trans}},i}$ from Eq. (13) into Eq. (5), the oscillation equations of the electrical network generators can be expressed as follows:

$$M_i{\ddot{\Theta }}_i=P_{m,i}-P_{L2G,i}-D_i{\dot{\Theta }}_i-\sum _{j=1}^{N_G}\left|V_i\right|\left|V_j\right|\left|Y_{\mathit{red},\mathit{ij}}\right|sin\left({\Theta }_i-{\Theta }_j+{\varphi }_{\mathit{ij}}\right)$$E16

It is evident that two nodes, denoted as $i$ and $j$, within the reduced network or dynamic model, will be linked by an edge if there is a path from node $i$ to node $j$ in the dynamic-algebraic model. Consequently, the reduced graph of a network with a connected graph will exhibit a complete graph structure. Additionally, node $i$ in the reduced network will have an internal loop if there is a path from node $i$ to a node connected to the reference node in the original network, and the reference node is not included as a separate node in the admittance matrix of the network $Y_{\mathit{net}}$. Hence, when the original graph is connected and one of the removed nodes is connected to the reference node, all nodes in the reduced network will be connected to the reference node.

3. Transient stability analysis

According to the previous section, by reducing the dimension of the dynamic-algebraic model of an electrical network to a dynamic model, the power transfer between nodes of the dynamic network will be a function of the sine of the phase angle difference between node voltages and the phase shift ${\varphi }_{\mathit{ij}}$. Based on Eq. (16), the swing equations of the generators in the electrical network can be presented as follows:

$$M_i{\ddot{\Theta }}_i+D_i{\dot{\Theta }}_i={\omega }_i-\sum _{j=1,j\ne i}^{N_G}{P}_{\mathit{ij}}sin\left({\Theta }_i-{\Theta }_j+{\varphi }_{\mathit{ij}}\right)i\in \left\{1,\dots ,N_G\right\}$$E17

Where, ${\omega }_i=P_{m,i}-P_{L2G,i}$, $P_{\mathit{ij}}=\left|Y_{\mathit{red},\mathit{ij}}\right|\mid V_i\Vert V_j\mid$, ${\varphi }_{\mathit{ij}}=-{\phi }_{\mathit{ij}}+\pi /2,{\varphi }_{\mathit{ij}}\in \left(0,\pi /2\right)$, and ${\phi }_{\mathit{ij}}=\text{atan}\left(\mathfrak{I}\left(Y_{\mathit{red},\mathit{ij}}\right)/\mathfrak{R}\left(Y_{\mathit{red},\mathit{ij}}\right)\right)$, $Y_{\mathit{red}}$ is the reduced admittance matrix of the network, ${\Theta }_i$ is the phase angle of voltage $V_i$, $P_{m,i}$ is the mechanical power injected into generator $i$, $P_{L2G,i}$ is the effect of network loads on generator $i$, and $N_G$ is the number of generators in the network.

The presented model in Eq. (17) corresponds to the second-order non-uniform Kuramoto model with a non-zero phase shift ${\varphi }_{\mathit{ij}}\ne 0$. Considering this model and the results obtained in the previous section, the transient stability of the electrical network is investigated in this chapter. Based on the obtained results, an explicit and direct method is presented to determine the necessary and sufficient condition for the frequency synchronization and phase coherence of the network generators. Dividing both sides of Eq. (17) by the inertia factors, the dynamic of the oscillatory agent $i$ in the electrical network can be obtained as follows:

$${\ddot{\Theta }}_i=\frac{{\omega }_i}{M_i}-\frac{D_i}{M_i}{\dot{\Theta }}_i-\sum _{j=1,j\ne i}^{N_G}\frac{P_{\mathit{ij}}}{M_i}sin\left({\Theta }_i-{\Theta }_j+{\varphi }_{\mathit{ij}}\right)$$E18

It is evident that the variation in the inertia values of network components leads to the asymmetry of coupling strengths in the dynamic model of oscillating agents. According to Eq. (18), for two oscillating agents $i$ and $j$ that are coupled, the coupling strength in the dynamics of agent $i$ is proportional to $P_{\mathit{ij}}/M_i$ while the coupling strength in the dynamics of agent $j$ is proportional to $P_{\mathit{ji}}/M_j$. This complexity arises from the analytical computations performed on the dynamic model of oscillating agent, considering the variation of inertia values. However, in order to analyze the dynamic behavior of the network model and determine sufficient conditions for phase coherence, it is necessary to specify the dynamic model of the frequency differences of oscillating agents in the electrical network. Through employing Eq. (18), the dynamic model of the network can be derived as follows:

$$\frac{d}{\mathit{dt}}\left(H\dot{\Theta }\right)=H{M}^{-1}\omega -H{M}^{-1}D\dot{\Theta }-H{M}^{-1}{H}^{T}diag\left(P_{\mathit{ij}}\right)\mathbf{sin}\left(\mathit{H\Theta }+\varphi \right)$$E19

where $\varphi =diag\left({\varphi }_{\mathit{ij}}\right){\mathbf{1}}_{\left|\epsilon \right|}\in {\mathbb{R}}^{\left|\epsilon \right|\times 1}$ is equivalent to a column vector. As mentioned, except for [27], the cases where the necessary and sufficient condition for frequency synchronization and phase coherence of the first or second-order Kuramoto model is explicitly presented in the form of analytical methods, are related to conditions where symmetry is established in agent communications, and references [12, 33, 52, 53, 54] are among the most important of these references. The theorem presented in [27] also provides a sufficient condition for phase coherence of agents in an electrical network with an analytical approach. However, the results of this reference are reliable only under the condition that the value of $M_i/D_i$ for small network generators is considered.

In this regard, and with the aim of providing a comprehensive theory for phase coherence of agents in an electrical network, first, the necessary condition for frequency synchronization of the network based on physical and topological characteristics is presented in the form of theorem 1. Then, using two approaches, the direct Lyapunov method and the compressibility property, a sufficient condition for boundedness of the phase difference between generators is extracted in sections 3–2 and 3–3, respectively, in the form of theorems 2 and 3.

3.1 Necessary condition for synchronization

In continuation of studies [11, 12, 13, 26, 33],this section provides the “necessary condition” for synchronization in the electrical network with significant losses in transmission lines and an asymmetric communication matrix.

Theorem 1. Consider the power network oscillators’ phase difference dynamic model Eq. (17), where the coupling strength $\mathrm{P}={\mathrm{P}}^{\mathrm{T}}$ and $B=H^T\in {\mathbb{R}}^{N_G\times \left|\epsilon \right|}$ is the incidence matrix of the complete graph. The symmetric phase shift ${\mathrm{\varphi }}_{\mathrm{ij}}$ and the phase ${\mathrm{\Theta }}_{\mathrm{i}}$ are located in the regions ${\varphi }_{\mathit{ij}}\in \left(0,\pi /2\right)$and$\left\{\Theta \left(t\right)\in ∆\left(\rho \right):\rho \in \left(0,\pi /2\right)\right\}$. In the synchronization conditions, the network states are located on the synchronous manifold $A=\left\{\left(\Theta ,f\right)\forall \theta \in ∆\left(\rho \right),f={\omega }_{\mathit{syn}}\right\}$, and the dynamic model will be simplified to the following algebraic equation:

$$D_i{\omega }_{\mathit{syn}}={\omega }_i-\sum _{j=1,j\ne i}^{N_G}{P}_{\mathit{ij}}sin\left({\Theta }_i-{\Theta }_j+{\varphi }_{\mathit{ij}}\right)i\in \left\{1,\dots ,N_G\right\}$$E20

By defining the vector X with element i as $X_i={\sum }_{j=1,j\ne i}^{N_G}{P}_{\mathit{ij}}sin\left({\varphi }_{\mathit{ij}}\right)cos\left({\Theta }_i-{\Theta }_j\right)$ and parameters $S_{\omega }={\sum }_{i=1}^{N_G}{X}_i$and${\varphi }_{\max }=max\left\{{\varphi }_{\mathit{ij}}\right\}$, the following propositions will hold:

Statement 1) The necessary condition for synchronization in an electrical network, assuming that the phase of all agents is limited to the set $\left\{\Theta \left(t\right)\in ∆\left(\rho \right):\rho \in \left(0,\pi /2\right)\right\}$, will be as follows:

$$\begin{array}{c}\left({\mathrm{{\rm Y}}}_{\max }-{\mathrm{{\rm Y}}}_{\min }\right)\le K_{\mathit{suf}}\left({deg}_i+{deg}_j\right)\mathit{\text{for}}\mathit{all}\left\{i,j\right\}\in \epsilon {\mathrm{{\rm Y}}}_{\max }\\ =max\left\{\frac{{\omega }_i}{D_i}\right\},\\ {\mathrm{{\rm Y}}}_{\min }=min\left\{\frac{{\omega }_i}{D_i}\right\},K_{\mathit{suf}}=sin\left(\rho +{\varphi }_{\max }\right)/D_{\min }\mathit{\text{for}}\left(\rho +{\varphi }_{\max }\right)\\ \le \left(\frac{\pi }{2}\right)K_{\mathit{suf}}=1/D_{\min }\mathit{\text{for}}\left(\rho +{\varphi }_{\max }\right)\ge \left(\frac{\pi }{2}\right)\end{array}$$E21

Statement 2) The actual frequency of oscillations ${\omega }_{\mathit{\text{real}}}$, is determined by Eq. (22), where $\omega$ corresponds to the nominal frequency of the network.

$${\omega }_{\mathit{\text{real}}}=\omega +{\omega }_{\mathit{syn}},{\omega }_{\mathit{syn}}=\sum _{i=1}^{N_G}{\omega }_i-S_{\omega }/\sum _{i=1}^{N_G}{D}_i$$E22

Proof of Statement 1: Synchronization in Eq. (20) requires that the maximum power transfer between nodes through inter-agent communications is greater than the maximum difference in input power to the agents. Substituting the maximum value of the term $sin\left({\Theta }_i-{\Theta }_j+{\varphi }_{\mathit{ij}}\right)$ and $\left({\mathrm{{\rm Y}}}_{\max }-{\mathrm{{\rm Y}}}_{\min }\right)$ on the right side of the this model, we will reach to the minimum connectivity strength to satisfy the synchronization condition. Considering the assumption of Theorem 1 regarding the phase of the agents $\left({\Theta }_1,\dots ,{\Theta }_{N_G}\right)\in {\mathbb{T}}^{N_G}$ belonging to the region $∆\left(\rho \right)\subset {\mathbb{T}}^1\forall \rho \in \left[0,\pi /2\right]$, it is possible to conservatively replace the term $sin\left({\Theta }_i-{\Theta }_j+{\varphi }_{\mathit{ij}}\right)$ with $sin\left(\rho +{\varphi }_{\max }\right)$ or $sin\left(\pi /2\right)$. Therefore, by satisfying the inequality $\left({\mathrm{{\rm Y}}}_{\max }-{\mathrm{{\rm Y}}}_{\min }\right)\le K_{\mathit{suf}}\left({deg}_i+{deg}_j\right)$, frequency synchronization will be achieved in the network. End of proof of Statement 1.

Proof of Statement 2: In a synchronized network, Eq. (20) holds for all $N_G$ generators. Applying trigonometric transformations and summing both sides of the equation for $N_G$, we will have ${\omega }_{\mathit{syn}}{\sum }_{i=1}^{N_G}{D}_i$ on the left side. On the right side, due to the oddness of the sine function and the symmetry of phase shifts between agents, the expression ${\sum }_{i=1}^{N_G}{\omega }_i-{\sum }_{i=1}^{N_G}{X}_i$ is obtained. So, the actual frequency of the network will be given by Eq. (22). End of proof of Statement 2.

Result 1: In a network with losses or non-zero phase shift in the dynamic model, the angle ${\varphi }_{\mathit{ij}}$ leads to network oscillation at a frequency other than the nominal frequency $\omega$. Through constant input power to the oscillatory agents while changing the phase shift ${\varphi }_{\mathit{ij}}$, the value of the network oscillation frequency will change as shown in Figure 1.

Figure 1.

Effect of phase shift on the rotational frequency of agents in the second-order non-uniform Kuramoto model.

Result 2: For a network with homogeneous agents $D_i=D_j=D$ and uniform connectivity strength given by $P_{\mathit{ij}}=P/N_G$, the necessary condition for synchronization is derived as $P\ge \frac{N_G\left({\omega }_{\max }-{\omega }_{\min }\right)}{2\left(N_G-1\right)sin\left(\rho +{\varphi }_{\max }\right)}$. For a lossless network ${\varphi }_{\mathit{max}}=0$ and $\rho =\pi /2$, this condition simplifies to $P\ge \frac{N_G\left({\omega }_{\max }-{\omega }_{\min }\right)}{2\left(N_G-1\right)}$, which is equivalent to the condition presented in [33] and Theorem 2.1 in Ref. [12]. This condition further reduces to $P\ge \left({\omega }_{\max }-{\omega }_{\min }\right)$ for $N_G=2$ and $P\ge \frac{\left({\omega }_{\max }-{\omega }_{\min }\right)}{2}$ as $N_G\to \infty$ approaches infinity. Additionally, assuming zero phase shift, the condition becomes $P\ge \frac{N_G\left({\omega }_{\max }-{\omega }_{\min }\right)}{2\left(N_G-1\right)sin\left(\rho \right)}$, which is more conservative compared to the result presented in Theorem 3.1 in Ref. [13].

3.2 Sufficient condition for synchronization - LaSalle invariance principle

Theorem 2: Consider the power network oscillators’ dynamic model (11), where the coupling strength$P=P^T$and$B=H^T\in {\mathbb{R}}^{N_G\times \left|\epsilon \right|}$is the incidence matrix of the complete graph. For$\rho \in \left[0,\pi -\left(H{\overline{\theta }}_e+2H\overline{\vartheta }\right)\right]$, let the coupling topology between oscillators be such that the algebraic connectivity of the network graph${\lambda }_2\left(L\left(P_{\mathit{ij}}cos\left({\vartheta }_{\mathit{ij}}\right)\right)\right)$exceeds the${\lambda }_2^{\mathit{\text{Critical}}1}$or${\lambda }_2^{\mathit{\text{Critical}}2}$, which are defined as:

$${\lambda }_2\left(L\left(P_{\mathit{ij}}cos\left({\vartheta }_{\mathit{ij}}\right)\right)\right)\ge {\lambda }_2^{\mathit{\text{Critical}}1}=\frac{{max}_{i\ne j}\left\{D_i{D}_j\right\}{‖H{D}^{-1}\omega ‖}_2+d{‖diag\left(P_{\mathit{ij}}sin\left(|{\vartheta }_{\mathit{ij}}|\right)\right){\mathbf{1}}_{\mid \epsilon \mid }‖}_2}{\mathit{\kappa \rho }\left(d/N_G\right)sinc\left(\rho \right)}$$E23

$${\lambda }_2\left(L\left(P_{\mathit{ij}}cos\left({\vartheta }_{\mathit{ij}}\right)\right)\right)\ge {\lambda }_2^{\mathit{\text{Critical}}2}=\frac{{max}_{i\ne j}\left\{D_i{D}_j\right\}{‖H{D}^{-1}\omega ‖}_2{K}_{\max }}{\mathit{\kappa \rho }\left(d/N_G\right)sinc\left(\rho \right)}$$E24

Here,$K_{\max }={max}_{i\ne j}\left\{{\mathbf{1}}_{\mid \epsilon \mid }-diag\left(P_{\mathit{ij}}/P_{\mathit{ij}}^0\right)\mathbf{cos}\left(\rho \right)\right\}$, ${\theta }_e$ is equal with network equilibrium point, $H{\overline{\theta }}_e={max}_{i\ne j}\left\{|{\theta }_{e,i}-{\theta }_{e,j}|\right\}$,$H{\overline{\Theta }}^0={max}_{i\ne j}\left\{|{\Theta }_i^0-{\Theta }_j^0|\right\}$, $H\overline{\vartheta }=H{\overline{\Theta }}^0+{\varphi }_{\max }$, $P_{\mathit{ij}}^0$ is equal with network coupling strength before the transient, $0\le \kappa \le Q$, $d:={\sum }_{i=1}^{N_G}{D}_i$, $m:={\sum }_{i=1}^{N_G}{M}_i$, ${\vartheta }_{\mathit{ij}}={\Theta }_{\mathit{ij}}^0+{\varphi }_{\mathit{ij}}$, ${\Theta }_{\mathit{ij}}^0$ represents the initial phase difference of the generators, $Q=\sqrt{{\lambda }_U^{\min }/\left({\lambda }_U^{\max }+{\mathcal{S}}_{\max }\right)}$, ${\mathcal{S}}_{\max }=\left(\alpha +\beta \right)\left(m/2\right){max}_{i\ne j}\left\{P_{\mathit{ij}}\right\}$, ${\varphi }_{\max }={max}_{i\ne j}\left\{{\varphi }_{\mathit{ij}}\right\}$and${\lambda }_U^{\max }$and${\lambda }_U^{\min }$ are respectively the largest and smallest eigenvalues of the matrix $\mathrm{U}$ with the following relationship:

$$U=\left[\begin{array}{cc}diag\left(D_i{D}_j\right)& diag\left(D_i{D}_j\right)H{D}^{-1}M\\ {\mathit{MD}}^{-1}{H}^{T}diag\left(D_i{D}_j\right)& \alpha H^{T}diag\left(M_i{M}_j\right)H+\mathit{\beta mM}\end{array}\right]$$E25

And $\alpha$ and $\beta$ are also satisfied in the $\alpha >{max}_{i\ne j}\left\{D_j{M}_i/D_i{M}_j\right\}$and$\beta \ge \left(d/m\right)max\left\{M_i/D_i\right\}+\left(\alpha -1\right)\left(N_G-2\right)max\left\{M_i\right\}/m-\alpha$. (Then)

Statement 1: Phase Cohesiveness. The set ${\mathrm{\Omega }}_{\varsigma }=\left\{\left(\mathit{H\theta },\dot{\theta }\right)\in {\mathbb{R}}^{\mid \epsilon \mid \times N_G}:{‖\mathit{H\theta }\left(t\right)‖}_2\le \left(\kappa /Q\right)\rho \right\}$ is a positively invariant set for the electrical network, and under the given conditions, if the initial phase difference of the network generators satisfies $‖\mathit{H\theta }\left(t_0\right)‖=\mathit{\kappa \rho }\le \mathit{Q\rho }$, the network states will be uniformly bounded and the phase difference of the generators will remain within the region $‖\mathit{H\theta }\left(t\right)‖\le \left(\kappa /Q\right)\rho$.

Statement 2: Frequency Synchronization. Under the initial conditions mentioned above, the states of the oscillatory components in the dynamic model of the electrical network converge towards the largest positively invariant set given by $S_{\mathit{inv}}=\left(\mathit{H\theta }\left(t\right),\dot{\theta }\right)=\left(H{\theta }_e\left(t\right),{\mathit{\omega }}_{\mathit{syn}}\right)$, which is equivalent to achieving frequency synchronization in the electrical network in the form $\underset{t\to \infty }{lim}\mid {\dot{\theta }}_i\left(t\right)-{\dot{\theta }}_j\left(t\right)\mid =0,\forall i\ne j$, or${‖H\dot{\theta }\left(t\right)‖}_2=0$. In this set,${\mathit{\omega }}_{\mathit{syn}}={\mathbf{1}}_{N_G}{\omega }_{\mathit{syn}}$n, and${\omega }_{\mathit{syn}}$is obtained from Eq.(22). The largest positively invariant set with the equilibrium point of the system$\left(\mathit{H\theta }\left(t\right),\dot{\theta }\right)=\left(H{\theta }_e\left(t\right),{\mathit{\omega }}_{\mathit{syn}}\right)$is unique.

Proof: To prove this theorem, the Lyapunov function is defined as follows. The details of its proof are provided in [82, 83].

$$\begin{array}{c}{E}_1\left(\mathit{H\theta },H\dot{\theta }\right)=\frac{1}{2}\left[\begin{array}{cc}{\left(\mathit{H\theta }\right)}^T& {\dot{\theta }}^T\end{array}\right]\left[\begin{array}{cc}diag\left(D_i{D}_j\right)& diag\left(D_i{D}_j\right)H{D}^{-1}M\\ {\mathit{MD}}^{-1}{H}^{T}diag\left(D_i{D}_j\right)& \alpha H^{T}diag\left(M_i{M}_j\right)H+\mathit{\beta mM}\end{array}\right]\left[\begin{array}{c}\mathit{H\theta }\\ \dot{\theta }\end{array}\right]\\ E_2\left(\mathit{H\theta }\right)=\left(\alpha +\beta \right){\int }_0^{\mathit{H\theta }}diag\left(M_i{M}_j\right)\left[-H{M}^{-1}\omega +H{M}^{-1}{H}^{T}diag\left(P_{\mathit{ij}}\right)\mathbf{sin}\left(\mathit{H\theta }+\mathit{H\vartheta }\right)\right)]d\left(\mathit{H\theta }\right)\end{array}$$E26

Result 3 (The region of Attraction Development with Input Power Control): According to (23) and (24), it is evident that reducing the term$D_i{D}_j\mid {\omega }_i/D_i-{\omega }_j/D_j\mid$in the edges of the electrical network graph leads to a decrease in the term${max}_{i\ne j}\left\{D_i{D}_j\right\}{‖H{D}^{-1}\omega ‖}_2$and ultimately a decrease in${\lambda }_2^{\mathit{\text{Critical}}1}$and${\lambda }_2^{\mathit{\text{Critical}}2}$. The term${\omega }_i$in the electrical network is defined as${\omega }_i=P_{m,i}-P_{L2G,i}$, where$P_{m,i}$represents the mechanical power input to generator$i$and$P_{L2G,i}$represents the influence of network loads on generator$i$. Therefore, applying an appropriate control strategy to regulate parameters$P_{m,i}$and$P_{L2G,i}$in the network can provide a sufficient condition for synchronization.Figure 2 illustrates the variations of ${\lambda }_2^{\mathit{\text{Critical}}1}$ and ${\lambda }_2^{\mathit{\text{Critical}}2}$ with respect to changes in the term ${\mathit{max}}_{i\ne j}\left\{D_i{D}_j\right\}{‖H{D}^{-1}\omega ‖}_2$.

Figure 2.

Algebraic connectivity as a function of ${max}_{i\ne j}\left\{D_i{D}_j\right\}{‖H{D}^{-1}\omega ‖}_2$.

Result 4 (Impact of Oscillators’ Initial Phase on the Critical Coupling): According to the results, as the initial phase difference ${\mathrm{\Theta }}_{\mathrm{ij}}^0$ increases towards $\mathrm{\pi }/2$, the value of the term $cos\left({\mathrm{\vartheta }}_{\mathrm{ij}}\right)$ decreases and the value of the term $sin\left(|{\mathrm{\vartheta }}_{\mathrm{ij}}|\right)$ increases. Consequently, the values of ${\lambda }_2^{\mathit{\text{Critical}}1}$ and ${\lambda }_2^{\mathit{\text{Critical}}2}$λ decrease. Considering the above, in case of reducing the fault clearance time in the network, the stability conditions and the region of the attraction of the network will improve.

Result 5 (The region of attraction Development with Control of Agent Parameters): As the inertia and damping coefficient of generators in the electrical network become closer to a homogeneous state, the value of critical coupling decreases. In the extreme case where the inertia and damping coefficient are the same for all oscillatory agents, we have:

$${\lambda }_2^{\mathit{\text{Critical}}1}=\frac{{‖\mathit{H\omega }‖}_2+N_G{‖diag\left(P_{\mathit{ij}}sin\left(|{\vartheta }_{\mathit{ij}}|\right)\right){\mathbf{1}}_{\mid \epsilon \mid }‖}_2}{\mathit{\kappa \rho }sinc\left(\rho \right)}$$E27

$${\lambda }_2^{\mathit{\text{Critical}}2}=\frac{{‖\mathit{H\omega }‖}_2{K}_{\max }}{\mathit{\kappa \rho }sinc\left(\rho \right)}$$E28

In this case, the value of critical coupling is not dependent on the parameters of the generators in the network. Figure 3 illustrates the variations of ${\mathrm{\lambda }}_2^{\text{Critical}1}$ and ${\mathrm{\lambda }}_2^{\text{Critical}2}$ for two homogeneous and heterogeneous networks with respect to changes in $\mathrm{\rho }$. It is evident that the value of critical coupling in the homogeneous network is smaller than that in the heterogeneous network.

Figure 3.

Algebraic connectivity in two “heterogeneous” and “homogenous” networks (abbreviated to “Het” and “Hom”, respectively).

Result 6 (The effect of transmission line resistance in the dynamic of an electrical network): Employing trigonometric transformation, it is evident that the phase shift${\varphi }_{\mathit{ij}}$in the electrical network model not only affects the magnitude of the connections between network generators from$P_{\mathit{ij}}$to$P_{\mathit{ij}}cos\left({\varphi }_{\mathit{ij}}\right)$, but also adds the term$-{\sum }_{j=1}^{N_G}{P}_{\mathit{ij}}sin\left({\varphi }_{\mathit{ij}}\right)cos\left({\Theta }_i-{\Theta }_j\right)<0,$to the$i$th generator node. Both terms have a negative impact on the stability conditions of the network. Therefore, it can be concluded that increasing the resistance of transmission lines makes it more difficult to establish the sufficient condition for phase coherence among network generators.

Table 1 presents a recursive algorithm for estimating the critical clearing time $t_c$ under conditions where an error occurs in the electrical network topology.

Stage	Function
1	Obtain input information for each bus in the network, including specifications of network generators, transmission line impedances, load specifications, transformer specifications, and convert them into an appropriate format for AC load flow analysis.
2	Perform AC load flow analysis and determine the values of network’s state
3	Apply the Kron reduction technique and develop the dynamic model of the network with a complete graph.
4	Develop network topology in pre-fault, during fault, and post-fault conditions based on the postulated scenario.
5	Determine the time step size $t_s={0.001}^{\ast }\mathit{sec}$ for numerical simulations and set the fault removal time as $t_c=t_s$.
6	Perform simulations for the state equations from the time of fault occurrence until fault removal time $t_c$ to calculate initial values of state variables.
7	Calculate initial conditions for network state variables and calculate parameters $Q$, ${\mu }_{c1}={\lambda }_2^{\mathit{\text{Critical}}1}$, ${\mu }_{c2}={\lambda }_2^{\mathit{\text{Critical}}2}$, ${\lambda }_2\left(L\left(P_{\mathit{ij}}cos\left({\vartheta }_{\mathit{ij}}\right)\right)\right)$, and determine $\rho$ at the fault clearing time.
8	If ${\mu }_{c1}\le {‖\mathit{H\theta }\left(t_0\right)‖}_2\le \mathit{Q\rho }$ or ${\mu }_{c2}\le {‖\mathit{H\theta }\left(t_0\right)‖}_2\le \mathit{Q\rho }$, then the network state variables are in the region of attraction at the fault clearing time. The network maintains its transient stability and proceed to step 9.; If ${‖\mathit{H\theta }\left(t_0\right)‖}_2\ge \mathit{Q\rho }$, the network cannot maintain stability condition at the fault removal time $t_c$. However, the network has maintained transient stability condition in the previous step, and the critical clearing time or CCT is equal to $t_{\mathit{cr}}=t_c-t_s$. Proceed to step 10.
9	Update the fault removal time as $t_c=t_c+t_s\mathit{sec}$ sec and go back to step 6.
10	Determine the value of CCT as $t_{\mathit{cr}}$ in the output.
* Since the period of the first swing angle of the rotor generator is approximately 2–4 seconds, a time step of $t_s=0.001$ has been chosen.

Table 1.

Estimating the region of attraction and critical clearing time $t_c$ based on the initial condition.

3.3 Sufficient condition for synchronization –contraction property

The challenges in finding a Lyapunov function and applying the Lyapunov method to analyze the stability of certain systems, especially nonlinear multi-agent systems with time-varying communications, demonstrate the need for more advanced tools and methods in studying these systems. In [84, 85, 86, 87, 88, 89, 90, 91, 92], several examples of these tools have been developed in the form of practical theories for analyzing the stability properties of nonlinear monotone systems using max-separable Lyapunov functions. Conditions for convergence of states in multi-agent systems with nonlinear dynamics have been developed using analytical theories in the field of non-smooth analysis. In this section, we invoke the contraction property introduced in [93] for consensus protocols in self-organizing systems. According to the mentioned theory, in order to achieve convergence of agent states to a common value, it is necessary for $max\left\{x_1,\dots ,x_n\right\}$ to be a non-increasing function of time and at the same time, $min\left\{x_1,\dots ,x_n\right\}$ should not decrease. Therefore, by combining these two characteristics in [93], the Lyapunov function $V\left(x\right)=max\left\{x_1,\dots ,x_n\right\}-min\left\{x_1,\dots ,x_n\right\}$ or in another form, $V\left(x\right)=max\left\{|x_i-x_j||i,j\in \left\{1,\dots ,n\right\}\right\}$ has been introduced to analyze the stability of the consensus protocol. Since the max-separable Lyapunov function is not necessarily continuously differentiable, the direct Lyapunov method cannot be applied to this class of Lyapunov functions.

In [84], with reference to [94, 95], it has been shown that for a system $\dot{x}=f\left(t,x\right),f:\mathbb{R}\times \mathcal{D}\to {\mathbb{R}}^n$ where $\mathcal{D}\subset {\mathbb{R}}^n$ is a set and $V\left(t,x\right):\mathbb{R}\times \mathcal{D}\to \mathbb{R}$ is a continuous function that satisfies the Lipschitz condition, the Upper Right Dini Derivative bound of the function $V\left(t,x\right)$ along the state trajectories $f\left(t,x\right)$ can be presented as the $D_f^{+}V$ function, as follows:

$$D_f^{+}V\left(t,x\right)=\underset{\tau \to 0^{+}}{limsup}\frac{V\left(t+\tau ,x\left(t,\tau \right)\right)-V\left(t,x\right)}{\tau }=\underset{\tau \to 0^{+}}{limsup}\frac{V\left(t+\tau ,x+\mathit{\tau f}\left(t,x\right)\right)-V\left(t,x\right)}{\tau }$$E29

In this case, considering $V_i\left(t,x\left(t\right)\right):\mathbb{R}\times \mathcal{D}\to \mathbb{R}$ as a class $C^1$ function and defining $V\left(t,x\left(t\right)\right)=\underset{i\in \left\{1,\dots ,n\right\}}{max}{V}_i\left(t,x\left(t\right)\right)$, the equality $D^{+}V\left(t,x\left(t\right)\right)=\underset{i\in \left\{1,\dots ,n\right\}}{max}{\dot{V}}_i\left(t,x\left(t\right)\right)$ holds. The stability and control analysis theories derived from this approach can be applied in areas related to distributed decision-making, coordination, synchronization, and coherence of states in multi-agent systems and high-dimensional systems. The development of Theorem 3 is based on this approach.

Lemma 2: We recall the dynamic model of a lossy power network (19), where the network connectivity matrix $P=P^T$ is a complete graph, the symmetric phase shift ${\varphi }_{\mathit{ij}}$ and phase ${\Theta }_i$ are restricted to the regions ${\varphi }_{\mathit{ij}}\in \left(0,\pi /2\right)$ and $\left\{\Theta \left(t\right)\in ∆\left(\rho \right):\rho \in \left(0,\pi /2\right)\right\}$, and the network coupling matrix is denoted as $B=H^T\in {\mathbb{R}}^{N_G\times \left|\epsilon \right|}$. The upper bound values for the oscillators’ frequency and its rate of change are given by ${\dot{\theta }}_{\max }={\omega }_{ِ\mathit{Dif}}+c$and${\ddot{\theta }}_{\max }=\left(\frac{D_{\min }+D_{\max }}{M_{\min }}\right){\omega }_{ِ\mathit{Dif}}+\frac{D_{\max }}{M_{\min }}c$, where${\omega }_{\mathit{syn}}^0=\left({\sum }_{i=1}^{N_G}{\omega }_i-S_{\omega }^0\right)/{\sum }_{i=1}^{N_G}{D}_i$, $S_{\omega }^0={\sum }_{i=1}^{N_G}{X}_i^0$, $X_i^0={\sum }_{j=1,j\ne i}^{N_G}{P}_{\mathit{ij}}^{0}sin\left({\varphi }_{\mathit{ij}}\right)cos\left({\Theta }_i^0-{\Theta }_j^0\right)$, ${\omega }_{\mathit{Dif}}=\left({\omega }_{\max }+sin\left(H\overline{\vartheta }\right){deg}_{\max }\right)/D_{\min }$, $c\ge \Vert {\omega }_{\mathit{Dif}}\mid -\mid {\omega }_{\mathit{syn}}^0\Vert$, $H\overline{\vartheta }=H{\overline{\Theta }}^0+{\varphi }_{\max }$, and$H{\overline{\Theta }}^0={max}_{i\ne j}\left\{|{\Theta }_i^0-{\Theta }_j^0|\right\}$. Additionally, ${\Theta }_i^0$ and $P_{\mathit{ij}}^0$ represent the initial phase of agent $i$ and the initial strength of communication between agents $i$ and $j$ at time $t=0$, (respectively)

Proof: The proof of this lemma is provided in [96].

In continuation of the studies [12] (Theorem 4.1), [13] (Theorem 3.1), [27], [33] (Theorem 2), [34, 68, 69, 70, 71], this section provides a sufficient condition for achieving synchrony in a lossy electrical network with heterogeneous oscillators, non-uniform inter-agent connectivity matrix. The results of this theorem also hold for networks with uncertainty in the agents’ input power in the range ${\omega }_i\in \left\{{\omega }_{\min }-{\omega }_{\max }\right\}$.

Theorem 3: We consider the dynamic model of an electrical network with non-zero phase shifts Eq.(17). If the following condition is satisfied in the network graph:

$$\begin{array}{ll}{N}_G\underset{i\ne j}{min}\left\{\frac{P_{\mathit{ij}}}{D_i}cos\left({\varphi }_{\mathit{ij}}\right)\right\}sin\left(\rho \right)\ge & \underset{i\ne j}{max}\left|\frac{{\omega }_i}{D_i}-\frac{{\omega }_j}{D_j}\right|+2{\ddot{\Theta }}_{\max }\underset{i\in \left\{1,\dots ,N_G\right\}}{max}\left\{\frac{M_i}{D_i}\right\}\\ & +N_G\underset{i\ne j}{max}\left\{\frac{P_{\mathit{ij}}}{D_i}sin\left({\varphi }_{\mathit{ij}}\right)\right\}-N_G\underset{i\ne j}{min}\left\{\frac{P_{\mathit{ij}}}{D_i}sin\left({\varphi }_{\mathit{ij}}\right)\right\}cos\left(\rho \right)\end{array}$$E30

Then, the set$\Theta \in ∆\left(\rho \right)$is an attractive region for the electrical network with respect to$\rho \in \left[0,\pi /2-{\varphi }_{\max }\right]$, and under the condition that the initial phases of the agents fall within this non-empty set, phase coherence will be established in the network.

Proof: To prove this Theorem, the max-separable Lyapunov function $V\left(\Theta \left(t\right)\right)=\underset{i,j\in \left\{1,\dots ,N_G\right\}}{max}{V}_{\mathit{ij}}\left(\Theta \left(t\right)\right)={\Theta }_p\left(t\right)-{\Theta }_q\left(t\right)$ is used. The phases ${\Theta }_p$ and ${\Theta }_q$ correspond to the oscillators ${\mathit{Osc}}_p$ and ${\mathit{Osc}}_q$ at the boundary positions of the arc $\rho$, as shown in Figure 4. In order for $∆\left(\rho \right)$ to be invariant, it is necessary that the geodesic distance between the two boundary agents decreases. This distance is equal to the maximum $V_{\mathit{ij}}\left(\Theta \left(t\right)\right)$, which is why we define the candidate Lyapunov function $V\left(\Theta \left(t\right)\right)=\underset{i,j\in \left\{1,\dots ,N_G\right\}}{max}{V}_{\mathit{ij}}\left(\Theta \left(t\right)\right)={\Theta }_p\left(t\right)-{\Theta }_q\left(t\right)$ on ${\mathbb{T}}^1$. The function $V\left(\Theta \left(t\right)\right)$ is not continuously differentiable, but it satisfies the Lipschitz continuity conditions. Therefore, according to the Lemma 2.2 in [84], to ensure the decremental nature of $V\left(\Theta \left(t\right)\right)$, it is necessary that the upper bound of the Dini Derivative of $V\left(\Theta \left(t\right)\right)$ along the trajectory of the system Eq. (17) is negative.

Figure 4.

Network agents’ positions in an extreme state.

$$D_f^{+}V\left(\Theta \left(t\right)\right)=\underset{\tau \to 0^{+}}{limsup}\frac{V\left(\Theta \left(t+\tau \right)\right)-V\left(\Theta \left(t\right)\right)}{\tau }={\dot{\Theta }}_p\left(t\right)-{\dot{\Theta }}_q\left(t\right)<0$$E31

The details of the proof of this Theorem are presented in Ref. [96].

Result 7 (The region of attractive expansion through power control): Considering the sufficient condition for phase coherence of the agents given in Eq. (30), it is clear that reducing the differential term $\underset{\mathrm{i}\ne \mathrm{j}}{max}\mid {\mathrm{\omega }}_{\mathrm{i}}/{\mathrm{D}}_{\mathrm{i}}-{\mathrm{\omega }}_{\mathrm{j}}/{\mathrm{D}}_{\mathrm{j}}\mid$ in the edges of the graph associated with the dynamic model of the electrical network leads to a decrease in the right-hand side of the Eq. (30). In other words, the heterogeneity of ${\mathrm{\omega }}_{\mathrm{i}}/{\mathrm{D}}_{\mathrm{i}}$ in each node of the network will have a negative impact on the stability condition. Figure 5 illustrates the effect of ${\mathrm{\omega }}_{\mathrm{i}}/{\mathrm{D}}_{\mathrm{i}}$ heterogeneity on the critical angle ${\mathrm{\varphi }}_{\mathrm{cri}}$ (for a network with a fixed topology, the values of the left-hand side and the right-hand side of Eq. (30) will be equal if the phase shift in the coupling function is equal to ${\mathrm{\varphi }}_{\mathrm{cri}}$) in a sample network. According to Eq. (17), the value of ${\mathrm{\omega }}_{\mathrm{i}}$ in the dynamic model of the electrical network is given by ${\mathrm{\omega }}_{\mathrm{i}}={\mathrm{P}}_{\mathrm{m},\mathrm{i}}-{\mathrm{P}}_{\mathrm{L}2\mathrm{G},\mathrm{i}}$, where ${\mathrm{P}}_{\mathrm{m},\mathrm{i}}$ represents the mechanical power input to the $\mathrm{i}$-th generator and ${\mathrm{P}}_{\mathrm{L}2\mathrm{G},\mathrm{i}}$ represents the impact of network loads on the $\mathrm{i}$-th generator after applying the Kron reduction method. Thanks to the control systems in the electrical network, the value of the mechanical input power and partially the power consumption at the network loads are adjustable. By implementing appropriate control strategies to regulate these parameters, it is possible to make the synchronization condition stated in Eq. (30) more achievable.

Figure 5.

The value of ${\varphi }_{\mathit{cri}}$ as a function of $\underset{i\ne j}{max}\mid {\omega }_i/D_i-{\omega }_j/D_j\mid .$

Result 8: The term$2{\ddot{\Theta }}_{\max }\underset{i\in \left\{1,\dots ,N_G\right\}}{max}\left\{M_i/D_i\right\}$on the right-hand side of (30) indicates that for two networks with the same network graph and generators’ parameters$M_i$and$D_i$being pairwise equal, it will be easier to satisfy the phase coherence condition of the agents in the network with a smaller${\ddot{\Theta }}_{\mathit{max}}$.

Result 9: As the inertia and damping coefficients of the generators approach a homogeneous state, the magnitudes of the terms$\underset{i\ne j}{min}\left\{P_{\mathit{ij}}/D_{i}sin\left({\varphi }_{\mathit{ij}}\right)\right\}cos\left(\rho \right)$and$\underset{i\ne j}{max}\left\{P_{\mathit{ij}}/D_{i}sin\left({\varphi }_{\mathit{ij}}\right)\right\}$become closer to each other, and the value of$\underset{i\ne j}{max}\left|{\omega }_i/D_i-{\omega }_j/D_j\right|$also decreases. Considering this aspects, having a large standard deviation of the generator parameters reduces the right-hand side of the inequality Eq.(30). In the extreme case where the inertia and damping coefficients of all network generators are the same and$\rho =\pi /2-{\varphi }_{\max }$, the condition will simplifies to:

$$\begin{array}{c}{N}_G\underset{i\ne j}{min}\left\{P_{\mathit{ij}}cos\left({\varphi }_{\mathit{ij}}\right)\right\}cos\left({\varphi }_{\max }\right)\ge \underset{i\ne j}{max}\left|{\omega }_i-{\omega }_j\right|+2M{\ddot{\Theta }}_{\max }+N_{G}sin\left({\varphi }_{\max }\right)\\ (\underset{i\ne j}{max}\left\{P_{\mathit{ij}}\right\}-\underset{i\ne j}{min}\left\{P_{\mathit{ij}}sin\left({\varphi }_{\mathit{ij}}\right)\right\}\end{array}$$E32

Result 10: The closer the coupling strength between network generators becomes homogeneous, it will be easier to achieve the sufficient condition for synchronization. In the extreme case, for a weakened network with homogeneous connections and parameters and $\rho =\pi /2-\varphi$, the sufficient condition will be as follows:

$$P\ge \frac{\underset{i\ne j}{max}\left|\frac{{\omega }_i}{D_i}-\frac{{\omega }_j}{D_j}\right|+2{\ddot{\Theta }}_{\max }\underset{i\in \left\{1,\dots ,N_G\right\}}{max}\left\{\frac{M_i}{D_i}\right\}}{\underset{i\ne j}{min}\left\{1/D_i\right\}-\underset{i\ne j}{max}\left\{sin\left(\varphi \right)/D_i\right\}}$$E33

Furthermore, if the desired network also consists of homogeneous oscillators, the condition simplifies to:

$$P\ge \frac{\underset{i\ne j}{max}\left|{\omega }_i-{\omega }_j\right|+2{\ddot{\Theta }}_{\max }M}{1-sin\left(\varphi \right)}$$E34

Result 11: The arc consisting the phase of agents,$\rho \in \left[0,\pi /2-{\varphi }_{\max }\right]$, is introduced in the form of two terms$sin\left(\rho \right)$and$cos\left(\rho \right)$in the inequality Eq.(30). It is evident that as the value of $\mathrm{\rho }$ increases towards $\pi /2-{\varphi }_{\max }$, the right-hand side of the inequality decreases, while the left-hand side increases. In the extreme case of $\rho =\pi /2-{\varphi }_{\max }$, it becomes easier to satisfy the sufficient condition, and it transforms into the following equation:

$$\begin{array}{ll}{N}_G\underset{i\ne j}{min}\left\{\frac{P_{\mathit{ij}}}{D_i}cos\left({\varphi }_{\mathit{ij}}\right)\right\}cos\left({\varphi }_{\max }\right)\ge & \underset{i\ne j}{max}\left|\frac{{\omega }_i}{D_i}-\frac{{\omega }_j}{D_j}\right|+2{\ddot{\Theta }}_{\max }\underset{i\in \left\{1,\dots ,N_G\right\}}{max}\left\{\frac{M_i}{D_i}\right\}\\ & +N_G\underset{i\ne j}{max}\left\{\frac{P_{\mathit{ij}}}{D_i}sin\left({\varphi }_{\mathit{ij}}\right)\right\}\\ & -N_G\underset{i\ne j}{min}\left\{\frac{P_{\mathit{ij}}}{D_i}sin\left({\varphi }_{\mathit{ij}}\right)\right\}sin\left({\varphi }_{\max }\right)\end{array}$$E35

Result 12: According to Eq. (30), which shows the impact of transmission line resistance in the form of phase shift ${\varphi }_{\mathit{ij}}$, since the term $N_G\underset{i\ne j}{max}\left\{P_{\mathit{ij}}/D_{i}sin\left({\varphi }_{\mathit{ij}}\right)\right\}-N_G\underset{i\ne j}{min}\left\{P_{\mathit{ij}}/D_{i}sin\left({\varphi }_{\mathit{ij}}\right)\right\}cos\left(\rho \right)$ is always greater than zero, increasing the resistance of transmission lines leads to an increase in the value of phase shift ${\varphi }_{\mathit{ij}}$, making it more difficult to achieve the sufficient condition for the phase coherence of network generators. In the case where the resistance value is zero, condition in Eq. (30) simplifies to the following equation:

$$N_G\underset{i\ne j}{min}\left\{\frac{P_{\mathit{ij}}}{D_i}\right\}sin\left(\rho \right)\ge \underset{i\ne j}{max}\left|\frac{{\omega }_i}{D_i}-\frac{{\omega }_j}{D_j}\right|+2{\ddot{\Theta }}_{\max }\underset{i\in \left\{1,\dots ,N_G\right\}}{max}\left\{\frac{M_i}{D_i}\right\}$$E36

Result 13: In the networks of oscillators with the first order model $M=0$, homogenous agents $D_i=D_j$, and uniform connections $P_{\mathit{ij}}=P/N_G$, when the phase shift ${\varphi }_{\mathit{ij}}$ is equal to zero, the condition Eq. (30) is equivalent to the condition of Theorem 3.1 in [13] and is expressed as $P\ge \left({\omega }_{\max }-{\omega }_{\min }\right)/sin\left(\rho \right)$.

4. Conclusion

In this chapter, we present the estimation of the region of attraction in an electric network by applying the LaSalle invariance principle and the contraction property for consensus protocols in systems with autonomous agents. We detail this estimation in Theorem 2 and 3. Furthermore, Theorem 1 introduces the necessary condition for achieving synchronization based on the network graph and its parameters.

In Theorem 2, we explicitly express the estimation of the region of attraction of the electric network as a function of the underlying network parameters, the input power intensity of the network agents, the intensity of the coupling topology, and the initial conditions of the model. This estimation is presented in two distinct forms:

• The region of attraction based on the topology of network graph ${\lambda }_2\left(L\left(P_{\mathit{ij}}cos\left({\vartheta }_{\mathit{ij}}\right)\right)\right)\ge {\lambda }_2^{\mathit{\text{Critical}}}$.

• The region of attraction based on the network’s initial conditions ${\mu }_c\le {‖\mathit{H\theta }\left(t_0\right)‖}_2\le \mathit{Q\rho }$.

Furthermore, based on the results of Theorem 2, an algorithm for estimating the critical clearing time in the network is presented.

In Theorem 3, the obtained estimation for the region of attraction of the electric network is presented as an invariant set of agent phases and also an explicit function of the underlying network parameters, input power intensity of the network agents, the intensity of their interactions, and the initial conditions of the model.

Based on the results obtained from theorems 2 and 3, it is necessary to increase the distance between the value of ${\lambda }_2\left(L\left(P_{\mathit{ij}}cos\left({\vartheta }_{\mathit{ij}}\right)\right)\right)$ and the critical value in order to develop the estimated regions of attractions. This can be achieved by keeping ${\lambda }_2\left(L\left(P_{\mathit{ij}}cos\left({\vartheta }_{\mathit{ij}}\right)\right)\right)$ constant and reducing ${\lambda }_2^{\mathit{\text{Critical}}}$, or by keeping ${\lambda }_2^{\mathit{\text{Critical}}}$ constant and increasing ${\lambda }_2\left(L\left(P_{\mathit{ij}}cos\left({\vartheta }_{\mathit{ij}}\right)\right)\right)$. The development of the region of attraction is introduced in the following methods:

1. Increasing the connectivity intensity between the generators of the electric network, which consequently increases the value of the connectivity distance of the network, ${\lambda }_2\left(L\left(P_{\mathit{ij}}cos\left({\varphi }_{\mathit{ij}}\right)\right)\right)$.

2. Increasing the density of connections in the electric network (pre-reduction Kron topology). As a result of this process, the value of the algebraic connectivity of the network, ${\lambda }_2\left(L\left(P_{\mathit{ij}}cos\left({\varphi }_{\mathit{ij}}\right)\right)\right)$ (increases)

3. Having the term $D_i{D}_j\mid {\omega }_i/D_i-{\omega }_j/D_j\mid$ in the network generators $i\in \left\{1,\dots ,N_G\right\}$ close to each other (details in the result of 6–2).

4. The inertia and damping coefficients of the generators should be closer to a state of homogeneity.

Based on the details discussed in this chapter, and the analytical presented theories, the following topics can be suggested for further exploration within the scope of this chapter. These topics are based on providing analytical results for second-order Kuramoto oscillator networks under different operating conditions and improving the capabilities of electric networks to maintain stability.

4.1 Further network analysis with more details

• Investigation and analysis of the impact of uncertainty in the parameters of oscillatory agents (damping, inertia, input power or natural frequency, and the intensity of inter-agent connection) on the stability conditions addressed in this chapter for second-order Kuramoto oscillator networks.

• Examination and analysis of the effect of the delay in the second-order Kuramoto oscillator model on the stability condition and identification of the permissible delay before the model reaches an unstable state (necessary and sufficient conditions for stability in networks with phase delay should be provided).

4.2 Moving from network analysis towards control theory formulation

• Enhancement of the region of attraction and transient response in second-order Kuramoto oscillator networks through the design of optimized controllers.

• Improvement of the estimated the region of attraction and transient response in second-order Kuramoto oscillator networks with time-varying connecting topology through the design of adaptive controllers to control the input power to the oscillators.

• Controller design for electric networks to improve the region of attraction and transient response of the network by mapping the network model to a second-order Kuramoto oscillator model and identifying the critical clearing time based on the derived stability condition.

• Analysis of the impact of delay or pocket loss in the communication feedback on the stability conditions of electric networks.

4.3 Numerical analysis

• Performing numerical analyses and investigating the impact of oscillatory agent parameters and inter-agent connections on the stability conditions of oscillatory networks.