- Email: [email protected]
Frequency stability in modern power network from complex network viewpoint
Physica A xxx (xxxx) xxx
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Frequency stability in modern power network from complex network viewpoint ∗
Hai-Peng Ren a , , Yuan Gao a , Long Huo a , Ji-hong Song a , Celso Grebogi a,b a
Shaanxi Key Laboratory of Complex System Control and Intelligent Information Processing, Xi’an University of Technology, Xi’an, Shaanxi, 710048, China b Institute for Complex Systems and Mathematical Biology, King’s College, University of Aberdeen, Aberdeen AB24 3UE, UK
article
info
Article history: Received 21 July 2018 Received in revised form 23 October 2019 Available online xxxx Keywords: Renewable Energy and Storage Multi-agent system Frequency synchronization Stability
a b s t r a c t Grid-connected operation of Renewable Energy and Storage (RES) nodes make the dynamics of modern power grid to be more complex. A model of power grid, considering RES nodes, is being proposed to address frequency synchronization and stability analysis. First, a unified dynamical model of four different types of nodes are established according to the swing equation, including the RES nodes with power-frequency droop inverter controllers, and the storage nodes with charging and discharging states. Second, based on the multi-agent linear time-varying protocol, the stability of the frequency synchronization solution of power system is given by a theorem, which is significant for theoretical progress and practical engineering applications. Third, verification and validation are carried out by simulation experiments and case study in local grid such as Western System Coordinating Council and the Shaanxi North Power Grid model. © 2019 Elsevier B.V. All rights reserved.
1. Introduction To solve the problem of energy consumption and carbon dioxide emission caused by the traditional fossil energy generation, the proportion of renewable energy generation is being increased in power grids. Grid-connected operation of intermittent energy, like solar energy or wind energy, leads to fluctuations in energy supplying and demanding, in spite of energy storage systems reduce such fluctuations. The modern power grid, therefore, is a time-varying, nonlinear coupled network, heterogeneously composed of Renewable Energy and Storage (RES) nodes and the conventional generation nodes [1]. Therefore the traditional frequency synchronization methods are difficult to apply. The synchronization concept can be traced back to 1665, when C. Huygens found that coupled pendula would synchronously swing. The Kuramoto model, proposed by biologist Winfree and Kuramoto, has made much contribution to the synchronization research. A review [2] summarized a series of research results about the Kuramoto oscillators, including the stability analysis. Kuramoto model with a bimodal parameter distribution was proposed to describe the synchronous oscillation of the generator and load in [3]. However, there were the following shortcomings: (1) Lack of theoretical proofs about the stability and convergence of the frequency synchronization point; (2) The coupling weights of each edge and damping constant of each node is unnecessarily assumed to be equal; (3) The analysis complexity increases with the presence of higher order derivatives; (4) The model does not take RES nodes into account. The non-uniform Kuramoto model and the multi-agent consensus protocol were used to derive the sufficient condition for the frequency synchronization of power grid in [4]. However, the frequency stability analysis of power grid in [4] only considered the network connectivity and the initial phase of the ∗ Corresponding author. E-mail address: [email protected] (H.-P. Ren). https://doi.org/10.1016/j.physa.2019.123558 0378-4371/© 2019 Elsevier B.V. All rights reserved.
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
2
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
nodes, without taking into account the effect of the Laplacian matrix on the stability. The critical coupling strength for frequency synchronization in power grid was obtained in [5] and [6] in 2012 by assuming that the coupling weights of each edge and the damping constant of each node are equal and did not consider the RES nodes as well [5,6], which was not agreeable with practical grids. The output of the generation nodes were controlled cooperatively according to the information exchanged with the ambient agent and the multi-agent consensus protocol, in order to guarantee the frequency stability of the grid [7]. This article considered the dynamics of renewable energy equipped with droop controller, but did not consider the storage nodes with different working states. A necessary and sufficient condition for the existence of frequency synchronization solution in micro-grid with droop-controlled inverters were derived in [8]. However, all loads were assumed to be constant power loads in [8]. The methods of analyzing power grid frequency synchronization are mainly divided into two approaches: the grid topology and the grid dynamics. From the grid topology approach, the synchronous ability of the various grid sizes and the evolution patterns are primarily compared and analyzed using numerical experiments and statistical physics methods [9,10], which demonstrates the static characteristics in an analytical way, but cannot reflect the dynamical behavior of the power grid. In this paper, to analyze frequency stability from the dynamics approach, we establish an unified dynamical model for the power grid including RES nodes, which simplifies the analysis complexity about the power grid frequency synchronization. We show that the frequency stability of power grid is influenced not only by the network connectivity and the initial phase of the nodes, but also by the time-varying Laplacian matrix of the power grid. Ref. [11] applies a set-valued Lyapunov approach to consider discrete-time consensus algorithms with unidirectional time-dependent communication links. Ref. [4] gives the condition for the existence of frequency synchronized solution, but the result cannot guarantee the stability of the synchronized solution. This paper provides a theorem to decide whether the synchronized solution is stable by estimating the eigenvalues of the time-varying Laplacian matrix. Thereby, we derive the sufficient conditions for frequency synchronization and its stability in power grids that contains RES nodes. Finally, the conditions are verified by the Western System Coordinating Council model and the Shaanxi North Power Grid. The remainder of this paper is organized as follows: in Section 2, some relevant preliminaries and notations are introduced. In Section 3, the uniform mathematical power grid model for four different types of nodes are given. In Section 4, the sufficient condition for frequency synchronization and stability of power grids are derived. The simulation results of both the Western System Coordinating Council and the Shaanxi North Power Grid are given in Section 5. Finally, some conclusions are drawn in Section 6. 2. Preliminaries and notation 2.1. The set of phase differences and positive invariance Geometry on the n-torus: S = (−π, +π] denotes a torus so that any phase angle θ belongs to S. A n-dimensional product set T n is defined as S × S ...S . Let |θ1 − θ2 | denote the geodesic distance between two angles θ1 , θ2 ∈ S. Let
n
∆ (γ ) be the set of angle array (θ1 , θ2 , . . . , θn ) that belongs to T n such that there exists an arc of length γ⏐ containing all ⏐ θ1 , θ2 , . . . , θn in its interior. That means that an array of angles θ ∈ ∆ (γ ) satisfies maxi,j∈{1,2,...,n} ⏐θi − θj ⏐ ⩽ γ , where θi is a function of time t. At each time t, there exists an arc of length γ containing all angles θi (t), if for t > ts , where ts is constant, γ is no longer increasing, then the set of phase differences ∆ (γ ) is defined as positively invariance [4]. 2.2. Algebraic graph theory [12] An arbitrary complex network can be represented by the relevant weighted graph. A weighted directed graph is a triple set G = (v, ε , A), where v = {1, 2, . . . , n} is the set of nodes, ε = {(i, j) |i, j ∈ v} is the set of directed edges, and A is the adjacency matrix. A = [aij ], with aij > 0 for each directed edge (i, j) ∈ ε and otherwise aij =0. a ij is the coupling ∑ n strength between node i and j. Degree matrix De is a diagonal matrix whose main diagonal entries are j=1 aij . Then the Laplacian matrix is defined by L =De −A, and it is positive semidefinite with eigenvalues satisfying 0 = λ1 ≤ λ2 ≤ ...λ|n| , if the non-main diagonal elements of the Laplace matrix are less than or equal to zero [13]. 2.3. Linear time-varying consensus protocol of multi-agent system Before giving⏐ Lemma 1, we ⏐ define two concepts used in Lemma 1. Asymptotic consensus is defined as a state in a network where ⏐ξi (t ) − ξj (t )⏐ → 0 as t → ∞ from any initial condition ξi (0), ξj (0), i, j = 1, . . . , n, n is the node number in the network. Spanning tree is defined as a graph connection where every certex (node) is connected each other but without any cycle. Lemma 1 ([14]). Consider ξi is the state variable of agent i, aij (t) is the (i, j) entry of the adjacency matrix of the associated graph at time t. Assume that each node satisfies the following equation,
ξ˙i (t) = −
n ∑
aij (t)(ξi (t) − ξj (t)),
(1)
j=1
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
3
Eq. (1) can be written in matrix form as
ξ˙ (t) = −L(t)ξ (t),
(2)
where ξ = [ξ1 , . . . ξn ] is the state and L = [lij ] is the Laplacian matrix, with lii = j̸ =i aij (t), and lij = −aij (t), j ̸ = i. Matrix L(t) in (2) is time-varying. The system (2) achieves asymptotic consensus if the collection of directed interaction graphs across some time intervals has a spanning tree frequently enough as the system evolves. T
∑
∑n
Lemma 2 ([14]). Given a matrix B = [bij ] where bii < 0, bij > 0, ∀i ̸ = j, and j=1 bij = 0 for each i, then B has at least one zero eigenvalue and all of the nonzero eigenvalues are in the open left half plane. Furthermore, B has exactly one zero eigenvalue if and only if the directed graph associated with B has a spanning tree. According to Lemmas 1 and 2, we learn that, with the network connectivity varying, if −L(t) has exactly one zero eigenvalue and all of the nonzero eigenvalues are in the open left half plane, then, with the network connectivity varying, the system (2) achieve asymptotic consensus, that is, the asymptotic synchronization of the nodes in the network is achieved. 3. Mathematical model of different nodes in a power grid 3.1. The dynamical model of the generator The dynamical model of the generator is given as follows
⎛ Di θ˙i = Pmi − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(θi − θj )⎠ ,
(3)
j=1
where, in terms of node i, Pmi is the mechanical power of generator, Di is the damping constants, respectively, θi of node i is the angle of the shaft in radians with respect to a synchronously rotating reference, and Ei is the voltage magnitude of the synchronous machine. Yij is the admittance of transmission line between node i and node j, and Gii is the self-conductance of node i. More details about the model parameters and derivation of Eq. (3) could be found in Appendix A. 3.2. The dynamical node of the load Similar to the generator, the dynamics of the load is given as follows
⎛ Di θ˙i = −PLi − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(θi − θj )⎠ .
(4)
j=1
Note that, the main difference between (4) and (A.29) in Appendix A is the mechanical power. Since the load node is the energy consumption node, PLi represents the active power that the load node consumes. 3.3. Model of RES nodes based on droop control For grid-connected RES nodes, it is very necessary to convert the DC or the non-industrial frequency AC power to the industrial frequency AC power by means of inverters. Power grid typically contains a bank of RES nodes connected using converters equipped with power-frequency droop controllers, and operated in parallel. Therefore, we establish the dynamics of RES nodes based on the droop control. Droop control is a basic way of grid connected converters control, and is widely applied in the RES batteries converter technique [15]. The frequency droop control formula is fi = f0i − npi Pi ,
(5)
where the parameter npi = Xi /2π Ei E is referred to as the droop coefficient, Ei is the output voltage of inverter connected to the RES nodes, E is the voltage of DC bus, jXi is the line impedance between inverter and bus, f0i is the output frequency of the unload inverter, and Pi is the active power of the RES node. Proofs about droop control can be found in Appendix B. For simplicity, let Di = 1/npi , Eq. (5) is rewritten as Di fi = Di f0i − Pi .
(6)
In addition, we denote fi by θ˙i , the nominal power Di f0i by PDi , and the output power Pi is treated as done in (A.26) in Appendix A. Then Eq. (6) is rewritten as
⎛ Di θ˙i = PDi − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(θi − θj )⎠ .
(7)
j=1
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
4
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Eq. (7) is the model of RES nodes for discharging state. The storage node works not only in discharging state, but also in charging state that can be treated as a load node, Thus, its model can be described by (4). 4. Frequency synchronization and stability analysis of a power grid based on complex network theory Based on the theory in Section 2.2, we define the power grid including n nodes and m transmission lines denoted as
|v| = n, |ε| = m. Then all the nodes can be denoted as {vG , vL , vD }, where vG represents the set of the generator nodes, vL represents the set of the load nodes, and vD represents the set of the RES nodes. In order to describe the various types of grid nodes using a uniform matrix form, we make the following transformation ∗ Pmi = Pmi − Ei 2 Gii , i ∈ vG ,
PLi∗ = −PLi − Ei 2 Gii , i ∈ vL ,
(8)
∗ PDi = ±PDi − Ei 2 Gii , i ∈ vD ,
when the storage node works in charging state, the sign before PDi should be negative, otherwise it should be positive. Then Eqs. (3), (4) and (7) can be written as ∗ − Di θ˙i = Pmi
n ∑
Ei Ej ⏐Yij ⏐ sin(θi − θj ), i ∈ vG ,
⏐ ⏐
j=1
Di θ˙i = PLi∗ −
n ∑
Ei Ej ⏐Yij ⏐ sin(θi − θj ), i ∈ vL ,
⏐ ⏐
(9)
j=1
∗ Di θ˙i = PDi −
n ∑
Ei Ej ⏐Yij ⏐ sin(θi − θj ), i ∈ vD .
⏐ ⏐
j=1
Define
⏐ ⏐ αij ≜ Ei Ej ⏐Yij ⏐ , (i, j) ∈ ε,
(10)
so we can rewrite the model in Eq. (9) in the following form Di θ˙i = Pi∗ −
n ∑
) ( αij sin θi − θj , i ∈ v.
(11)
j=1
To achieve a stable frequency synchronization solution in a power grid, first, it is necessary for the synchronization solution of Eq. (11) to exist, and second, the synchronization solution should guarantee to be stable. A solution θ of Eq. (11) is said to be synchronized if θ˙i = θ˙j for each node, so the frequency difference between two arbitrary nodes can be represented as follows
θ˙i − θ˙j =
Pi∗ Di
−
Pj∗ Dj
−
n ( ∑ αik
Di
k=1
sin (θi − θk ) −
αjk Dj
sin θj − θk
(
) )
, i, j ∈ v.
(12)
To achieve frequency synchronization in power grid, the coupling strength between nodes i and j has to dominate the network’s non-uniformity for all t, that is Pi∗ Di
−
Pj∗ Dj
<
n ( ∑ αik k=1
Di
By defining Γmin = n mini̸=j
sin (θi − θk ) − αij
{ } Di
αjk Dj
sin θj − θk
(
⏐ ∗ ⏐P
, Γcritical = maxi̸=j ⏐ Di − i
Γmin sin (γ ) ⩾ Γcritical .
) )
, i, j ∈ v.
(13)
Pj∗
⏐ ⏐ , we rewrite Eq. (13) as Dj ⏐ (14)
There exists the frequency synchronized solution if the inequalities in Eq. (14) is satisfied. The synchronization frequency can be obtained by summing all the rows of the matrix equation, Eq. (11), such that n ∑ i=1
Di θ˙i =
n ∑ i=1
Pi∗ +
n n ∑ ∑
( ) αij sin θi − θj .
(15)
i=1 j=1
Assuming that∑ there∑ is the synchronization frequency ωsync , the property of the sine function implies that ( ) ∑non-symmetric ∑n n n n ∗ the last term i=1 j=1 αij sin θi − θj is zero, then, ωsync = i=1 Pi / i=1 Di . Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
5
After obtaining the synchronization solution, we need to analyze the stability of this solution. To prove the stability of the synchronized solution of Eq. (11), we consider the derivative of Eq. (11)
θ¨i = −
n ∑ αij
(
Di
j=1
αij
By defining aij =
θ¨i = −
cos θi − θj
) θ˙i − θ˙j .
(16)
cos θi − θj , we rewrite Eq. (16) as
(
Di
n ∑
)(
)
aij θ˙i − θ˙j ,
(
)
(17)
j=1
where aij is the adjacency matrix. The matrix form of Eq. (17) is
˙ θ¨ = −L(t)θ,
(18)
) ) { } (∑n (∑n a − aij is the Laplacian matrix, diag where θ = [θ1 , . . . θn ] is the nodes status, while L(t) = diag j=1 aij { } ∑j=n 1 ij T
represents a diagonal matrix whose main diagonal entries are j=1 aij , and aij represents a matrix with all entries being aij . The form of Eq. (18) is similar to the form of the consensus protocol of the linear time invariant multi-agent system [14] shown in Eq. (2). So we use Lemmas 1 and 2 to analyze the stability of the frequency synchronization solution of a power grid. Theorem 1. Consider (Γcritical /Γmin ), ⏐ ∗ (11),P ∗ ⏐if θ (0) ∈ ∆ (γ ), γ ∈ [γmin , π/2], where γmin = arcsin { } the power system ∑ ∑ α ⏐ ⏐P Γmin = n mini̸=j Dij , Γcritical = maxi̸=j ⏐ Di − Dj ⏐, then there exists the synchronization frequency ωsync = ni=1 Pi∗ / ni=1 Di . i
i
j
Furthermore, assuming with the change of time, Laplacian matrix - L(t) has ) exactly { } one zero eigenvalue and all of the nonzero (∑ n eigenvalues locate in the open left half plane, where L(t) = diag j=1 aij − aij , then the synchronization frequency is stable, and the frequency θ˙i (t ) exponentially converges to the following common value θ˙∞ = ωsync . Proof. Note that
γmin = arcsin
(
Γcritical Γmin
)
.
(19)
Then, sin (γmin ) =
Γcritical . Γmin
(20)
Since for γ ∈ [γmin , π /2], sin(γ ) is strictly monotonously increasing on the interval [0, π/2), then sin (γ ) ⩾ sin (γmin ) =
Γcritical . Γmin
(21)
Due to Γmin > Γcritical > 0, Eq. (14) holds from the fact that it can be derived by multiplying both sides of Eq. (21) by Γmin , consequently, there exists a synchronization frequency. In the following, we will prove the synchronization frequency is asymptotic stable solution. According to Section 2.2, we have the Laplacian matrix:
−L (t ) =
⎧ n ∑ ( ) ⎪ αij ⎪ ⎪ cos θi − θj , i = j ⎨− j=1,j̸ =i
αij
⎪ ⎪ ⎪ ⎩
Di
Di
(22)
( ) cos θi − θj , i ̸ = j.
If
⎧ n ∑ ( ) ⎪ αij ⎪ ⎪ cos θi − θj ⩽ 0 ⎨− j=1,j̸ =i
⎪ ⎪ ⎪ ⎩
αij Di
Di
(23)
cos θi − θj ⩾ 0,
(
)
according to Section 2.2 and Lemma 2, we have Laplacian matrix - L(t) has at least one zero eigenvalue and all of the nonzero eigenvalues locate in the open left half plane. Furthermore, with the change of time, if - L(t) has exactly one zero eigenvalue, then system (16) achieves asymptotic consensus, that is the asymptotic synchronization of all nodes in the network is achieved. Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
6
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Fig. 1. The grid topology of the Western System Coordinating Council. Table 1 The nodes parameters of the Western System Coordinating Council. Node number
1 2 3 4 5 6
(Bus (Bus (Bus (Bus (Bus (Bus
1) 2) 3) 4) 5) 6)
Node power
Node damping
U(p.u.)
Node voltage
δ (◦ )
P(p.u.)
D(p.u.)
1.0400 1.0250 1.0250 0.9956 1.0127 1.0159
0 9.2800 4.6648 −3.9888 −3.6874 0.7275
16.30 8.00 0.85 −1.00 −0.90 −0.80
0.2 0.1 0.1 0.1 0.1 0.1
Table 2 The line parameters of the Western System Coordinating Council. Source node
Sink node
Branch resistance
Branch reactance
1 1 2 2 3 3
4 5 4 6 5 6
0.0100 0.0170 0.0320 0.0085 0.0390 0.0119
0.2034 0.2104 0.3433 0.2543 0.4099 0.2708
The conditions for the existence of synchronization solution given in [4] is described as follows: γ ∈ [γmin , π/2], γmin = arcsin (Γcritical /Γmin ), then the set ∆ (γ ) is positively invariant. However, Theorem 1 in this work shows that, when the conditions given in [4] are satisfied, the existence of synchronization solution of the system can be proved, but the stability of this solution cannot be guaranteed. We prove that if Laplacian matrix −L(t) associated with the power system has exactly one zero eigenvalue and all of the nonzero eigenvalues locate in the open left half plane, the synchronization frequency is stable. Our Theorem in this paper gives the stability of the frequency synchronization solution. The following simulation examples validate the applicability of this Theorem. 5. Model simulation and analysis The grid topology of the Western System Coordinating Council (WSCC 9-BUS) [16] is shown in Fig. 1, where buses 1, 2 and 3 connected to generators are treated as generator nodes, buses 5, 6 and 8 connected to loads can be viewed as load nodes, buses 4, 7 and 9 are transformer buses, and can be treated as single paths. We thus obtain a 6-node network model. Additionally, we assume that there are three generator nodes, where node 1 is the classical generator node, node 2 is the renewable energy node, and node 3 is the storage node that operates in discharging state. The nodes parameters and lines parameters of the Western System Coordinating Council in the simulations are given in Tables 1 and 2. From Table 2, the branch resistance value is an order of magnitude smaller ⏐ { than } the branch reactance value. ⏐ αij
⏐ P∗ = 598.9573, Γcritical = maxi̸=j ⏐ Di −
Pj∗
⏐ = 91.5, Dj ⏐ ⏐ ⏐ Γcritical γmin = arcsin Γ = 8.7872◦ and γ = maxi̸=j ⏐θi (0) − θj (0)⏐ = max{|0◦ − 9.28◦ | , . . . , |0.7275◦ − (−3.6874◦ )|} = min ◦ ◦ |9.28 − (−3.9888 )| = 13.269◦ . Since γ ∈ (γmin , π /2), according to Theorem 1, there exists the frequency synchronized solution, since Laplace matrix ( ) α L(t) is time-varying, its eigenvalues change as well, where aij = Dij cos θi − θj . Fig. 2(b) illustrates that some eigenvalues i of the −L(t) are far greater than zero after 9 s, which does not satisfy the stability condition given by Theorem 1. Fig. 2(c) We can calculate from the parameters to obtain: Γmin = n mini̸=j
(
)
Di
i
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
7
Fig. 2. Subplots (a), (b), (c) illustrate the nodes frequency variation, the eigenvalues of −L(t) and the node phase variation of the Western System Coordinating Council, respectively. Simulation results using D = [0.2, 0.1, 0.1, 0.1, 0.1, 0.1]T and P = [16.3, 8, 0.85, −1, −0.9, −0.8]T .
illustrates the phase differences divergence. Therefore, although there exists a frequency synchronization solution, the solution is unstable. Thus the system does not obtain the stable synchronization frequency as shown in Fig. 2(a). In Fig. 3 the simulation parameters are changed. The damping vector is D = [2, 1, 1, 1, 1, 1]T . The injection power P= { } α
[1.63, 0.8, 0.85, −1, −0.9, −0.8]T . The simulation results are shown in Fig. 3. In this case we have Γmin = n mini̸=j Dij = i ⏐ ∗ ⏐ ( ) P∗ ⏐ ⏐P Γ 59.986, Γcritical = maxi̸=j ⏐ Di − Dj ⏐ = 1.85, γmin = arcsin Γcritical = 1.77◦ and γ is the same as that in Fig. 2. i j min For γ ∈ (γmin , π /2), there exists a frequency synchronization solution. Moreover, −L(t) has exactly one and only
one zero eigenvalue after 3 s and others are smaller than zero, which satisfies the stability condition given by Theorem 1. Fig. 3(c) illustrates that the phase difference between arbitrary two nodes converges, that is, the synchronization frequency is stable for these parameters. Therefore, the solutions of the system exponentially converge to stable∑synchronization ∑n n ∗ frequency in Fig. 3(a). In another words, the frequency of all nodes converge to the stable value θ˙∞ = i=1 Pi / i=1 Di . The North Local Power Grid of Shaanxi Province with 58 generator nodes and 57 load nodes is used to further verify ⏐ ⏐ { } Theorem 1 in this paper. We figured out that Γmin = n mini̸=j
γmin = arcsin
(
Γcritical Γmin
)
αij Di
⏐ P∗ = 11.0764, Γcritical = maxi̸=j ⏐ Di − i
Pj∗ Dj
⏐ ⏐ = 4.3505,
⏐ ⏐ = 23.1246◦ and γ = maxi̸=j ⏐θi (0) − θj (0)⏐ = 60.2112◦ . Obviously, for γ ∈ (γmin , π/2), there exists a frequency synchronization solution. Moreover, as shown in Fig. 4(b), −L(t) has one and only one zero eigenvalue at each time point and others are smaller than zero, which satisfies the stability condition given in Theorem 1. Fig. 4(c) illustrates that the phase difference between arbitrary two nodes converges, that is the synchronization frequency is stable. Therefore, the solutions of the system converge to stable synchronization frequency as shown ( in Fig. 4(a). ⏐ ⏐ { exponentially } ) αij
⏐ P∗ = 11.0764, Γcritical = maxi̸=j ⏐ Di −
Pj∗
⏐ = 8.8901, γmin = arcsin ΓΓcritical = Dj ⏐ i min ⏐ ⏐ 53.3825◦ and γ = maxi̸=j ⏐θi (0) − θj (0)⏐ = 60(.2112◦ . According ) { } to Theorem 1, there exists the frequency synchronization ∑n solution, since Laplace matrix L(t) = diag is time-varying, its eigenvalues change as well, where j=1 aij − aij ( ) α aij = Dij cos θi − θj . Fig. 5(b) illustrates that some eigenvalues of −L(t) are greater than zero, which does not satisfy i In Fig. 5, we have Γmin = n mini̸=j
Di
the stability condition given in Theorem 1. As shown in Fig. 5(c), the phase differences are divergent. Therefore, although there does exist the frequency synchronization solution, the solution is unstable. Thus the system does not obtain the stable synchronization frequency as shown in Fig. 5(a). Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
8
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Fig. 3. Subplots (a), (b), (c) illustrate the nodes frequency variation, the eigenvalues of −L(t) and the node phase variation of Western System Coordinating Council, respectively. Simulation results using D = [2, 1, 1, 1, 1, 1]T and P = [1.63, 0.8, 0.85, −1, −0.9, −0.8]T .
6. Conclusions The models of generator and load nodes, according to rotor dynamics of the generator and the RES nodes and equipped with power-frequency droop inverter controllers, are given in this paper, where the storage nodes have two operating states. We proved that the stability of synchronization solutions depends on the eigenvalues of −L(t) given in Eq. (18), using the consensus protocol of the linear time-varying multi-agent system. With the change of time, if matrix −L(t) has exactly one zero eigenvalue and all of the nonzero eigenvalues locate in the open left half plane, then the synchronization solutions is stable. Simulation examples validate the applicability of this Theorem using the power grid of the Western System Coordinating Council and the North Local Power Grid of Shaanxi Province. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The work is supported in part by Key Program of Natural Science Foundation of Shaanxi Province (2016ZDJC-01), in part by Shaanxi Provincial Special Support Program for Science and Technology Innovation Leader, China, in part by the Shaanxi Industrial Key Project, China (2018GY-165), and in part by Key Lab Research Program from Department of Education of Shaanxi Provincial Government, China (18JS081). Appendix A The motion equation of the generator rotor is governed by J δ¨ = Ta ,
(A.1)
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
9
Fig. 4. Subplots (a), (b), (c) illustrate the nodes frequency variation, the eigenvalues of −L(t) with enlarged part to show the positive ones and the node phase variation of the North Local Power Grid of Shaanxi Province, respectively.
where δ is the mechanical angle of the shaft in radians with respect to a fixed reference, J is the moment of inertia of the rotor in kg .m2 , and Ta is the sum torque in Newton meters. J δ¨ = Tm − TL J δ¨ = Tm − (Te + Td ) .
(A.2)
The angular reference is chosen with respect to a synchronously rotating reference frame, then
δ = ωR t + θm , δ˙ = ωR + θ˙m = ωR + ωm and δ¨ = θ¨m = ω ˙ m,
(A.3) (A.4)
where θm is the mechanical angle of the shaft in radians with respect to a synchronously rotating reference frame, θ˙m is the angular velocity of the shaft, ωR is synchronously rotating angular velocity. Another form of Eq. (A.2) is obtained by multiplying its both sides with ωm , as shown by Eq. (A.5) J ωm θ¨m = Pm − PL M θ¨m = Pm − PL ,
(A.5)
where Pm = Tm ωm is the mechanical power, PL = TL ωm is the load power. The quantity J ωm is called the inertia constant 2 and is denoted by M, it is related to kinetic energy of the rotating masses Wk , where Wk = 12 J ωm . Then M can be calculated as follows M = J ωm =
2Wk
ωm
.
(A.6)
The electrical angle θe can be used to denote all of the equations above,
θ = θe =
p 2
θm ,
(A.7)
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
10
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Fig. 5. Subplots (a), (b), (c) illustrate the nodes frequency variation, the eigenvalues of −L(t) with enlarged part to show the positive ones and node phase variation of the North Local Power Grid of Shaanxi Province, respectively.
where p is the number of pole pairs of the (generator) motor, which is calculated as follows 120fR
p=
nR
,
(A.8)
where, nR is the rated shaft speed in rad/min, fR is the synchronously rotating frequency in Hz. The Eq. (A.1) can be equivalently written as
(
2J
)
θ¨ =
p
(
2J
)
p
ω˙ = Ta ( N · m).
(A.9)
Now, normalize Eq. (A.9) by dividing both sides by a constant equal to the rated torque at the rated speed SB3
TB =
ωmR
=
60SB3 2π nR
,
(A.10)
where SB3 is the rated three-phase power (VA), ωmR = 2π nR /60 is the rated angular velocity of the shaft. According to Eqs. (A.8)–(A.10), we have J π 2 n2R
(
)
900ωR SB3
ω˙ =
Ta TB
= Tau p.u.,
(A.11) (
746 W R2
Ta TB
where = Tau . The moment of inertia J = 500 unit of J. Then, by the following transformation 746 W R2 π 2 n2R
(
H≜
can be obtained by using the method in [16], and W R2 /g is the
(
)
)
550 g × 900ωR SB3 Wk =
g
)
1
2 Wk SB3
2 J ωm =
,
1 2
ω˙ = Tau p.u. ( 2)
×
746 W R 550 g
(A.12)
×
(2π nR )2 3600
(
=
746 W R
) 2
π 2 n2R
550 g × 1800
(A.13) (A.14)
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
11
we have
(
)
2Wk SB3 ωR
(
2H
)
ωR
ω˙ = Tau p.u..
(A.15)
ω˙ = Tau p.u..
(A.16)
For a classical model of the synchronous machine, recognizing that the angular speed is nearly constant, the accelerating power Pa in p.u. is nearly equal to the accelerating torque Ta [16], then the generator becomes
(
2H
)
ωR
ω˙ ∼ = Pa p.u..
(A.17)
In addition, the damping torque (damping power) includes both the mechanical and electrical damping, which is represented by Dω, then Eq. (A.17) can be rewritten as
(
2H
)
ωR
ω˙ = Pm − Pe − Pd = Pm − Pe − Dω pu,
(A.18)
where Pd = Dω is the damping power. The concrete formula of the electrical output power Pe of the generator is derived as follows. Consider a power system consisting of one machine connected to an infinite bus through a transmission line. A schematic representation of this case is shown in Fig. A.1(a), E is the internal voltage of the synchronous machine, L is the equivalent inductance of the transmission line, and V is the voltage of the infinite bus, which is used as reference. The equivalent electrical circuit of the system is given in Fig. A.1(b), E = E ̸ δ is the voltage vector of the synchronous machine and δ is the voltage phase, and V = V ̸ 0 is the voltage vector of the infinite bus, x′d is the transient reactance of the machine, V¯ t is the terminal voltage of the synchronous machine, which can be eliminated by using Y − ∆ transformation, Z¯s is the equivalent shunt impedance at the machine terminal, and Z¯TL is the series impedance of the transmission line. Another equivalent circuit of Fig. 1(b) is given in Fig. A.1(c) by a Y − ∆ transformation. Assumptions for deriving the model are the following: (1) The mechanical power input remains constant during the transient period. (2) Damping or asynchronous power is negligible. (3) The synchronous machine can be represented by a constant voltage source and its mechanical angle coincides with the electrical phase angle of the voltage. (4) If a local load is fed at the terminal voltage of the machine, it can be represented by a constant impedance to neutral point. In Fig. A.1(a) we define: Nodes 0, 1 and 2 are the reference nodes, the internal voltage node of the synchronous machine and the infinite bus node, respectively; y10 , y12 , y20 are the three admittance elements, respectively, derived by using a Y − ∆ transformation. The symbols with an arrow in the following equations are the vector variables corresponding to Fig. A.1(c). For node 1, according to Kirchhoff’s law, we obtain: I1 = (E − V) y12 + Ey10 = Ey12 − Vy12 + Ey10 = E (y12 + y10 ) − Vy12 .
(A.19)
Define: Y11 = Y11 ̸ θ11 = y12 + y10 , θ11 = arctan Y12 = Y12 ̸ θ12 = −y12 , θ12 = arctan
(
(
real(Y11 ) imag (Y11 )
real(Y12 ) imag (Y12 )
)
) (A.20)
.
Then, I1 = EY11 + VY12 = EY11 ej(δ+θ11 ) + V Y12 ejθ12 .
(A.21)
So the power input Pe at node 1 is Pe = real EI∗1
(
)
(A.22)
EI∗1 = Eejδ EY11 e−j(δ+θ11 ) + V Y12 e−jθ12 = E 2 Y11 e−jθ11 + EV Y12 e−j(θ12 −δ)
(A.23)
Pe = real EI1 = E Y11 cos θ11 + EV Y12 cos (θ12 − δ) .
(A.24)
(
(
)
∗
)
2
When the voltage angle of the infinite bus is different from zero, we can also get the corresponding Pe , given as Pe = real EI∗1 = E 2 Y11 cos θ11 + EV Y12 cos (θ12 + α − δ) .
(
)
Now, we get the compact formula of the electrical output power of the generator as Eq. (A.25). Define β = G11 ≜ Y11 cos θ11 , G11 is the self-conductance, then Eq. (A.25) is written as Pe = E 2 G11 + EV Y12 sin (δ − α + β) .
(A.25) π 2
− θ12 and (A.26)
Consequently, we can get the mathematical model of the generator node in the power system [16] as
(
2H
ωR
)
ω˙ = Pm − Pd − Pe = Pm − Dω − E 2 G11 − EV Y12 sin (δ − α + β) .
(A.27)
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
12
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Fig. A.1. One machine connected to an infinite bus through a transmission line.
In nodes power grid, the dynamic models of the generator nodes can be written as
(
2Hmi
ωR
)
⎛ ω˙ i = Pmi − Di ωi − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(δi − δj + ϕij )⎠ .
(A.28)
j=1
It can be also expressed as
⎛ Mi θ¨i = Pmi − Di θ˙i − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(θi − θj + ϕij )⎠ ,
(A.29)
j=1
where Mi θ¨i is the inertial term. When frequency synchronization analysis of Eq. (A.29) is carried out, the inertia term Mi θ¨i can be omitted due to the fact that it only affects the convergence time to synchronize the state but not its existence [16]. In addition, the resistive component is much smaller than the inductive component of the generator’s output impedance and the line impedance, then ϕij = 0. In consequence, the dynamic models of the generator nodes is given by
⎛ Di θ˙i = Pmi − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(θi − θj + ϕij )⎠ ,
(A.30)
j=1
where, in terms of node i, Di and Mi are the damping and inertia constants, respectively, θi of node i is the angle of the shaft in radians with respect to a synchronously rotating reference, ϕij is the phase shifts caused by transfer conductance, Ei is the voltage magnitude of the synchronous machine, Yij is the admittance of the transmission line between nodes i and j, and Gii ≜ Yii cos θi is the self-conductance. Appendix B The model for the droop controlled inverter connected RES node is the following. Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
13
Fig. B.1. The equivalent circuit of two inverters operating in parallel.
Fig. C.1. Schematic diagram of local power grid in northern Shaanxi.
Power grid typically contains a bank of RES nodes equipped with active power-frequency droop inverter controllers operating in parallel. For the convenience of analysis, using two inverters in parallel as an example, Fig. B.1 is the equivalent circuit for two inverters operating in parallel, where E1 , E2 are the input voltage and I1 , I2 are the input current of inverter 1 and 2, respectively. Ri is the input impedance of inverter i, jXi is the line impedance between inverter and bus, E is the voltage of load bus, I0 , Z0 are the load current and load impedance, respectively, and δ1 , δ2 are the voltage phase difference between inverter 1, 2 and the bus, respectively. We get the inverter output active power from Fig. B.1.
Pi =
Ri Ei E cos δi − Ri E 2 Xi2
+
R2i
+
Xi E i E Xi2 + R2i
sin δi , i = 1, 2.
(B.1)
Assume that the entire inverter system is inductive, the resistive portion of the inverter output impedance and line impedance is much smaller than the inductive portion, and assume that the two inverters output impedance are equal, i.e, R1 = R2 ≈ 0, then Eq. (B.1) can be equivalently transformed as Pi =
Ei E Xi
sin δi , i = 1, 2.
(B.2)
Compared with the load impedance, inverter source output reactance is very small, so the output voltage phase angle difference δi is also very small. Then we have sin δi ≈ δi , therefore Eq. (B.2) can equivalently be given by: Pi =
Ei E Xi
δi , i = 1, 2.
(B.3)
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
14
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Additionally, in a practical inverter system, the variation range of output voltage amplitude is not large, which can be approximately viewed as a constant, then the output power of inverter system is only related to the phase angle difference
δi
∆Pi =
Ei E Xi
∆δi , i = 1, 2.
(B.4)
Since the phase control is achieved by adjusting the output frequency, the output frequency change of the inverter is δi , so we can change the inverter output phase by adjusting the output frequency, thereby, altering the ∆fi = ∆2πωi = 2∆ π ∆t active power output of the inverter. Define Xi npi = , (B.5) 2 π Ei E where the parameter npi is referred to as the droop coefficient, then the frequency droop control formula can be obtained fi = f0i − npi Pi ,
(B.6)
where f0i is the output frequency of the inverter without load. Appendix C The topology of the local power grid in northern Shaanxi is given in Fig. C.1. References [1] C.C. Chu, H.C. Iu, Complex networks theory for modern smart grid spplications: a survey, IEEE J. Emerg. Sel. Top. Circuits Syst. 7 (2) (2017) 177–191. [2] F.A. Rodrigues, T.K.D.M. Peron, P. Ji, J. Kurths, The Kuramoto model in complex networks, Phys. Rep. 610 (2016) 1–98. [3] G. Filatrella, A.H. Nielsen, N.F. Pedersen, Analysis of a power grid using a Kuramoto-like model, Eur. Phys. J. B 61 (4) (2008) 485–491. [4] F. Dörfler, F. Bullo, Synchronization and transient stability in power networks and non-uniform Kuramoto oscillators, SIAM J. Control Optim. 50 (3) (2012) 1616–1642. [5] R. Carareto, M.S. Baptista, C. Grebogi, Natural synchronization in power-grids with anti-correlated units, Commun. Nonlinear Sci. Numer. Simul. 18 (4) (2012) 1035–1046. [6] S. Lozano, L. Buzna, A.D. Guilera, Role of network topology in the synchronization of power systems, Eur. Phys. J. B 231 (85) (2012) 231–238. [7] L.Y. Lu, C.C. Chu, Consensus-based secondary frequency and voltage droop control of virtual synchronous generators for isolated AC micro-grids, IEEE J. Emerg. Sel. Top. Circuits Syst. 5 (3) (2015) 443–455. [8] J.W.S. Porco, F. Dörfler, F. Bullo, Synchronization and power sharing for droop-controlled inverters in islanded microgrids, Automatica 49 (9) (2013) 2603–2611. [9] J. Wang, Y.H. Liu, Y. Jiao, H.Y. Hu, Cascading dynamics in congested complex networks, Eur. Phys. J. B 67 (1) (2009) 95–100. [10] L.X. Yang, J. Jiang, Impacts of link addition and removal on synchronization of an elementary power network, Physica A 479 (2017) 99–107. [11] L. Moreau, Stability of continuous-time distributed consensus algorithms, in: IEEE Conference on Decision and Control, Vol. 4, 2004, pp. 3998-4003. [12] C. Godsil, Algebraic Graph Theory, Springer, 2001. [13] F. Dörfler, F. Bullo, Kron, Reduction of graphs with applications to electrical networks, IEEE Trans. Circuits Syst. I Regul. Pap. 60 (1) (2013) 150–163. [14] Wei. Ren, R.W. Beard, Beard Distributed Consensus in Multi-Vehicle Cooperative Control, Springer, London, 2008. [15] X.F. Sun, Y.C. Hao, Q.F. Wu, X.Q. Guo, B.C. Wang, A multifunctional and wireless droop control for distributed energy storage units in islanded AC microgrid applications, IEEE Trans. Power Electron. 32 (1) (2017) 736–751. [16] P.M. Anderson, A.A. Fouad, Power System Control and Stability, Wiley, 2003.
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Frequency stability in modern power network from complex network viewpoint ∗
Hai-Peng Ren a , , Yuan Gao a , Long Huo a , Ji-hong Song a , Celso Grebogi a,b a
Shaanxi Key Laboratory of Complex System Control and Intelligent Information Processing, Xi’an University of Technology, Xi’an, Shaanxi, 710048, China b Institute for Complex Systems and Mathematical Biology, King’s College, University of Aberdeen, Aberdeen AB24 3UE, UK
article
info
Article history: Received 21 July 2018 Received in revised form 23 October 2019 Available online xxxx Keywords: Renewable Energy and Storage Multi-agent system Frequency synchronization Stability
a b s t r a c t Grid-connected operation of Renewable Energy and Storage (RES) nodes make the dynamics of modern power grid to be more complex. A model of power grid, considering RES nodes, is being proposed to address frequency synchronization and stability analysis. First, a unified dynamical model of four different types of nodes are established according to the swing equation, including the RES nodes with power-frequency droop inverter controllers, and the storage nodes with charging and discharging states. Second, based on the multi-agent linear time-varying protocol, the stability of the frequency synchronization solution of power system is given by a theorem, which is significant for theoretical progress and practical engineering applications. Third, verification and validation are carried out by simulation experiments and case study in local grid such as Western System Coordinating Council and the Shaanxi North Power Grid model. © 2019 Elsevier B.V. All rights reserved.
1. Introduction To solve the problem of energy consumption and carbon dioxide emission caused by the traditional fossil energy generation, the proportion of renewable energy generation is being increased in power grids. Grid-connected operation of intermittent energy, like solar energy or wind energy, leads to fluctuations in energy supplying and demanding, in spite of energy storage systems reduce such fluctuations. The modern power grid, therefore, is a time-varying, nonlinear coupled network, heterogeneously composed of Renewable Energy and Storage (RES) nodes and the conventional generation nodes [1]. Therefore the traditional frequency synchronization methods are difficult to apply. The synchronization concept can be traced back to 1665, when C. Huygens found that coupled pendula would synchronously swing. The Kuramoto model, proposed by biologist Winfree and Kuramoto, has made much contribution to the synchronization research. A review [2] summarized a series of research results about the Kuramoto oscillators, including the stability analysis. Kuramoto model with a bimodal parameter distribution was proposed to describe the synchronous oscillation of the generator and load in [3]. However, there were the following shortcomings: (1) Lack of theoretical proofs about the stability and convergence of the frequency synchronization point; (2) The coupling weights of each edge and damping constant of each node is unnecessarily assumed to be equal; (3) The analysis complexity increases with the presence of higher order derivatives; (4) The model does not take RES nodes into account. The non-uniform Kuramoto model and the multi-agent consensus protocol were used to derive the sufficient condition for the frequency synchronization of power grid in [4]. However, the frequency stability analysis of power grid in [4] only considered the network connectivity and the initial phase of the ∗ Corresponding author. E-mail address: [email protected] (H.-P. Ren). https://doi.org/10.1016/j.physa.2019.123558 0378-4371/© 2019 Elsevier B.V. All rights reserved.
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
2
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
nodes, without taking into account the effect of the Laplacian matrix on the stability. The critical coupling strength for frequency synchronization in power grid was obtained in [5] and [6] in 2012 by assuming that the coupling weights of each edge and the damping constant of each node are equal and did not consider the RES nodes as well [5,6], which was not agreeable with practical grids. The output of the generation nodes were controlled cooperatively according to the information exchanged with the ambient agent and the multi-agent consensus protocol, in order to guarantee the frequency stability of the grid [7]. This article considered the dynamics of renewable energy equipped with droop controller, but did not consider the storage nodes with different working states. A necessary and sufficient condition for the existence of frequency synchronization solution in micro-grid with droop-controlled inverters were derived in [8]. However, all loads were assumed to be constant power loads in [8]. The methods of analyzing power grid frequency synchronization are mainly divided into two approaches: the grid topology and the grid dynamics. From the grid topology approach, the synchronous ability of the various grid sizes and the evolution patterns are primarily compared and analyzed using numerical experiments and statistical physics methods [9,10], which demonstrates the static characteristics in an analytical way, but cannot reflect the dynamical behavior of the power grid. In this paper, to analyze frequency stability from the dynamics approach, we establish an unified dynamical model for the power grid including RES nodes, which simplifies the analysis complexity about the power grid frequency synchronization. We show that the frequency stability of power grid is influenced not only by the network connectivity and the initial phase of the nodes, but also by the time-varying Laplacian matrix of the power grid. Ref. [11] applies a set-valued Lyapunov approach to consider discrete-time consensus algorithms with unidirectional time-dependent communication links. Ref. [4] gives the condition for the existence of frequency synchronized solution, but the result cannot guarantee the stability of the synchronized solution. This paper provides a theorem to decide whether the synchronized solution is stable by estimating the eigenvalues of the time-varying Laplacian matrix. Thereby, we derive the sufficient conditions for frequency synchronization and its stability in power grids that contains RES nodes. Finally, the conditions are verified by the Western System Coordinating Council model and the Shaanxi North Power Grid. The remainder of this paper is organized as follows: in Section 2, some relevant preliminaries and notations are introduced. In Section 3, the uniform mathematical power grid model for four different types of nodes are given. In Section 4, the sufficient condition for frequency synchronization and stability of power grids are derived. The simulation results of both the Western System Coordinating Council and the Shaanxi North Power Grid are given in Section 5. Finally, some conclusions are drawn in Section 6. 2. Preliminaries and notation 2.1. The set of phase differences and positive invariance Geometry on the n-torus: S = (−π, +π] denotes a torus so that any phase angle θ belongs to S. A n-dimensional product set T n is defined as S × S ...S . Let |θ1 − θ2 | denote the geodesic distance between two angles θ1 , θ2 ∈ S. Let
n
∆ (γ ) be the set of angle array (θ1 , θ2 , . . . , θn ) that belongs to T n such that there exists an arc of length γ⏐ containing all ⏐ θ1 , θ2 , . . . , θn in its interior. That means that an array of angles θ ∈ ∆ (γ ) satisfies maxi,j∈{1,2,...,n} ⏐θi − θj ⏐ ⩽ γ , where θi is a function of time t. At each time t, there exists an arc of length γ containing all angles θi (t), if for t > ts , where ts is constant, γ is no longer increasing, then the set of phase differences ∆ (γ ) is defined as positively invariance [4]. 2.2. Algebraic graph theory [12] An arbitrary complex network can be represented by the relevant weighted graph. A weighted directed graph is a triple set G = (v, ε , A), where v = {1, 2, . . . , n} is the set of nodes, ε = {(i, j) |i, j ∈ v} is the set of directed edges, and A is the adjacency matrix. A = [aij ], with aij > 0 for each directed edge (i, j) ∈ ε and otherwise aij =0. a ij is the coupling ∑ n strength between node i and j. Degree matrix De is a diagonal matrix whose main diagonal entries are j=1 aij . Then the Laplacian matrix is defined by L =De −A, and it is positive semidefinite with eigenvalues satisfying 0 = λ1 ≤ λ2 ≤ ...λ|n| , if the non-main diagonal elements of the Laplace matrix are less than or equal to zero [13]. 2.3. Linear time-varying consensus protocol of multi-agent system Before giving⏐ Lemma 1, we ⏐ define two concepts used in Lemma 1. Asymptotic consensus is defined as a state in a network where ⏐ξi (t ) − ξj (t )⏐ → 0 as t → ∞ from any initial condition ξi (0), ξj (0), i, j = 1, . . . , n, n is the node number in the network. Spanning tree is defined as a graph connection where every certex (node) is connected each other but without any cycle. Lemma 1 ([14]). Consider ξi is the state variable of agent i, aij (t) is the (i, j) entry of the adjacency matrix of the associated graph at time t. Assume that each node satisfies the following equation,
ξ˙i (t) = −
n ∑
aij (t)(ξi (t) − ξj (t)),
(1)
j=1
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
3
Eq. (1) can be written in matrix form as
ξ˙ (t) = −L(t)ξ (t),
(2)
where ξ = [ξ1 , . . . ξn ] is the state and L = [lij ] is the Laplacian matrix, with lii = j̸ =i aij (t), and lij = −aij (t), j ̸ = i. Matrix L(t) in (2) is time-varying. The system (2) achieves asymptotic consensus if the collection of directed interaction graphs across some time intervals has a spanning tree frequently enough as the system evolves. T
∑
∑n
Lemma 2 ([14]). Given a matrix B = [bij ] where bii < 0, bij > 0, ∀i ̸ = j, and j=1 bij = 0 for each i, then B has at least one zero eigenvalue and all of the nonzero eigenvalues are in the open left half plane. Furthermore, B has exactly one zero eigenvalue if and only if the directed graph associated with B has a spanning tree. According to Lemmas 1 and 2, we learn that, with the network connectivity varying, if −L(t) has exactly one zero eigenvalue and all of the nonzero eigenvalues are in the open left half plane, then, with the network connectivity varying, the system (2) achieve asymptotic consensus, that is, the asymptotic synchronization of the nodes in the network is achieved. 3. Mathematical model of different nodes in a power grid 3.1. The dynamical model of the generator The dynamical model of the generator is given as follows
⎛ Di θ˙i = Pmi − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(θi − θj )⎠ ,
(3)
j=1
where, in terms of node i, Pmi is the mechanical power of generator, Di is the damping constants, respectively, θi of node i is the angle of the shaft in radians with respect to a synchronously rotating reference, and Ei is the voltage magnitude of the synchronous machine. Yij is the admittance of transmission line between node i and node j, and Gii is the self-conductance of node i. More details about the model parameters and derivation of Eq. (3) could be found in Appendix A. 3.2. The dynamical node of the load Similar to the generator, the dynamics of the load is given as follows
⎛ Di θ˙i = −PLi − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(θi − θj )⎠ .
(4)
j=1
Note that, the main difference between (4) and (A.29) in Appendix A is the mechanical power. Since the load node is the energy consumption node, PLi represents the active power that the load node consumes. 3.3. Model of RES nodes based on droop control For grid-connected RES nodes, it is very necessary to convert the DC or the non-industrial frequency AC power to the industrial frequency AC power by means of inverters. Power grid typically contains a bank of RES nodes connected using converters equipped with power-frequency droop controllers, and operated in parallel. Therefore, we establish the dynamics of RES nodes based on the droop control. Droop control is a basic way of grid connected converters control, and is widely applied in the RES batteries converter technique [15]. The frequency droop control formula is fi = f0i − npi Pi ,
(5)
where the parameter npi = Xi /2π Ei E is referred to as the droop coefficient, Ei is the output voltage of inverter connected to the RES nodes, E is the voltage of DC bus, jXi is the line impedance between inverter and bus, f0i is the output frequency of the unload inverter, and Pi is the active power of the RES node. Proofs about droop control can be found in Appendix B. For simplicity, let Di = 1/npi , Eq. (5) is rewritten as Di fi = Di f0i − Pi .
(6)
In addition, we denote fi by θ˙i , the nominal power Di f0i by PDi , and the output power Pi is treated as done in (A.26) in Appendix A. Then Eq. (6) is rewritten as
⎛ Di θ˙i = PDi − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(θi − θj )⎠ .
(7)
j=1
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
4
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Eq. (7) is the model of RES nodes for discharging state. The storage node works not only in discharging state, but also in charging state that can be treated as a load node, Thus, its model can be described by (4). 4. Frequency synchronization and stability analysis of a power grid based on complex network theory Based on the theory in Section 2.2, we define the power grid including n nodes and m transmission lines denoted as
|v| = n, |ε| = m. Then all the nodes can be denoted as {vG , vL , vD }, where vG represents the set of the generator nodes, vL represents the set of the load nodes, and vD represents the set of the RES nodes. In order to describe the various types of grid nodes using a uniform matrix form, we make the following transformation ∗ Pmi = Pmi − Ei 2 Gii , i ∈ vG ,
PLi∗ = −PLi − Ei 2 Gii , i ∈ vL ,
(8)
∗ PDi = ±PDi − Ei 2 Gii , i ∈ vD ,
when the storage node works in charging state, the sign before PDi should be negative, otherwise it should be positive. Then Eqs. (3), (4) and (7) can be written as ∗ − Di θ˙i = Pmi
n ∑
Ei Ej ⏐Yij ⏐ sin(θi − θj ), i ∈ vG ,
⏐ ⏐
j=1
Di θ˙i = PLi∗ −
n ∑
Ei Ej ⏐Yij ⏐ sin(θi − θj ), i ∈ vL ,
⏐ ⏐
(9)
j=1
∗ Di θ˙i = PDi −
n ∑
Ei Ej ⏐Yij ⏐ sin(θi − θj ), i ∈ vD .
⏐ ⏐
j=1
Define
⏐ ⏐ αij ≜ Ei Ej ⏐Yij ⏐ , (i, j) ∈ ε,
(10)
so we can rewrite the model in Eq. (9) in the following form Di θ˙i = Pi∗ −
n ∑
) ( αij sin θi − θj , i ∈ v.
(11)
j=1
To achieve a stable frequency synchronization solution in a power grid, first, it is necessary for the synchronization solution of Eq. (11) to exist, and second, the synchronization solution should guarantee to be stable. A solution θ of Eq. (11) is said to be synchronized if θ˙i = θ˙j for each node, so the frequency difference between two arbitrary nodes can be represented as follows
θ˙i − θ˙j =
Pi∗ Di
−
Pj∗ Dj
−
n ( ∑ αik
Di
k=1
sin (θi − θk ) −
αjk Dj
sin θj − θk
(
) )
, i, j ∈ v.
(12)
To achieve frequency synchronization in power grid, the coupling strength between nodes i and j has to dominate the network’s non-uniformity for all t, that is Pi∗ Di
−
Pj∗ Dj
<
n ( ∑ αik k=1
Di
By defining Γmin = n mini̸=j
sin (θi − θk ) − αij
{ } Di
αjk Dj
sin θj − θk
(
⏐ ∗ ⏐P
, Γcritical = maxi̸=j ⏐ Di − i
Γmin sin (γ ) ⩾ Γcritical .
) )
, i, j ∈ v.
(13)
Pj∗
⏐ ⏐ , we rewrite Eq. (13) as Dj ⏐ (14)
There exists the frequency synchronized solution if the inequalities in Eq. (14) is satisfied. The synchronization frequency can be obtained by summing all the rows of the matrix equation, Eq. (11), such that n ∑ i=1
Di θ˙i =
n ∑ i=1
Pi∗ +
n n ∑ ∑
( ) αij sin θi − θj .
(15)
i=1 j=1
Assuming that∑ there∑ is the synchronization frequency ωsync , the property of the sine function implies that ( ) ∑non-symmetric ∑n n n n ∗ the last term i=1 j=1 αij sin θi − θj is zero, then, ωsync = i=1 Pi / i=1 Di . Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
5
After obtaining the synchronization solution, we need to analyze the stability of this solution. To prove the stability of the synchronized solution of Eq. (11), we consider the derivative of Eq. (11)
θ¨i = −
n ∑ αij
(
Di
j=1
αij
By defining aij =
θ¨i = −
cos θi − θj
) θ˙i − θ˙j .
(16)
cos θi − θj , we rewrite Eq. (16) as
(
Di
n ∑
)(
)
aij θ˙i − θ˙j ,
(
)
(17)
j=1
where aij is the adjacency matrix. The matrix form of Eq. (17) is
˙ θ¨ = −L(t)θ,
(18)
) ) { } (∑n (∑n a − aij is the Laplacian matrix, diag where θ = [θ1 , . . . θn ] is the nodes status, while L(t) = diag j=1 aij { } ∑j=n 1 ij T
represents a diagonal matrix whose main diagonal entries are j=1 aij , and aij represents a matrix with all entries being aij . The form of Eq. (18) is similar to the form of the consensus protocol of the linear time invariant multi-agent system [14] shown in Eq. (2). So we use Lemmas 1 and 2 to analyze the stability of the frequency synchronization solution of a power grid. Theorem 1. Consider (Γcritical /Γmin ), ⏐ ∗ (11),P ∗ ⏐if θ (0) ∈ ∆ (γ ), γ ∈ [γmin , π/2], where γmin = arcsin { } the power system ∑ ∑ α ⏐ ⏐P Γmin = n mini̸=j Dij , Γcritical = maxi̸=j ⏐ Di − Dj ⏐, then there exists the synchronization frequency ωsync = ni=1 Pi∗ / ni=1 Di . i
i
j
Furthermore, assuming with the change of time, Laplacian matrix - L(t) has ) exactly { } one zero eigenvalue and all of the nonzero (∑ n eigenvalues locate in the open left half plane, where L(t) = diag j=1 aij − aij , then the synchronization frequency is stable, and the frequency θ˙i (t ) exponentially converges to the following common value θ˙∞ = ωsync . Proof. Note that
γmin = arcsin
(
Γcritical Γmin
)
.
(19)
Then, sin (γmin ) =
Γcritical . Γmin
(20)
Since for γ ∈ [γmin , π /2], sin(γ ) is strictly monotonously increasing on the interval [0, π/2), then sin (γ ) ⩾ sin (γmin ) =
Γcritical . Γmin
(21)
Due to Γmin > Γcritical > 0, Eq. (14) holds from the fact that it can be derived by multiplying both sides of Eq. (21) by Γmin , consequently, there exists a synchronization frequency. In the following, we will prove the synchronization frequency is asymptotic stable solution. According to Section 2.2, we have the Laplacian matrix:
−L (t ) =
⎧ n ∑ ( ) ⎪ αij ⎪ ⎪ cos θi − θj , i = j ⎨− j=1,j̸ =i
αij
⎪ ⎪ ⎪ ⎩
Di
Di
(22)
( ) cos θi − θj , i ̸ = j.
If
⎧ n ∑ ( ) ⎪ αij ⎪ ⎪ cos θi − θj ⩽ 0 ⎨− j=1,j̸ =i
⎪ ⎪ ⎪ ⎩
αij Di
Di
(23)
cos θi − θj ⩾ 0,
(
)
according to Section 2.2 and Lemma 2, we have Laplacian matrix - L(t) has at least one zero eigenvalue and all of the nonzero eigenvalues locate in the open left half plane. Furthermore, with the change of time, if - L(t) has exactly one zero eigenvalue, then system (16) achieves asymptotic consensus, that is the asymptotic synchronization of all nodes in the network is achieved. Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
6
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Fig. 1. The grid topology of the Western System Coordinating Council. Table 1 The nodes parameters of the Western System Coordinating Council. Node number
1 2 3 4 5 6
(Bus (Bus (Bus (Bus (Bus (Bus
1) 2) 3) 4) 5) 6)
Node power
Node damping
U(p.u.)
Node voltage
δ (◦ )
P(p.u.)
D(p.u.)
1.0400 1.0250 1.0250 0.9956 1.0127 1.0159
0 9.2800 4.6648 −3.9888 −3.6874 0.7275
16.30 8.00 0.85 −1.00 −0.90 −0.80
0.2 0.1 0.1 0.1 0.1 0.1
Table 2 The line parameters of the Western System Coordinating Council. Source node
Sink node
Branch resistance
Branch reactance
1 1 2 2 3 3
4 5 4 6 5 6
0.0100 0.0170 0.0320 0.0085 0.0390 0.0119
0.2034 0.2104 0.3433 0.2543 0.4099 0.2708
The conditions for the existence of synchronization solution given in [4] is described as follows: γ ∈ [γmin , π/2], γmin = arcsin (Γcritical /Γmin ), then the set ∆ (γ ) is positively invariant. However, Theorem 1 in this work shows that, when the conditions given in [4] are satisfied, the existence of synchronization solution of the system can be proved, but the stability of this solution cannot be guaranteed. We prove that if Laplacian matrix −L(t) associated with the power system has exactly one zero eigenvalue and all of the nonzero eigenvalues locate in the open left half plane, the synchronization frequency is stable. Our Theorem in this paper gives the stability of the frequency synchronization solution. The following simulation examples validate the applicability of this Theorem. 5. Model simulation and analysis The grid topology of the Western System Coordinating Council (WSCC 9-BUS) [16] is shown in Fig. 1, where buses 1, 2 and 3 connected to generators are treated as generator nodes, buses 5, 6 and 8 connected to loads can be viewed as load nodes, buses 4, 7 and 9 are transformer buses, and can be treated as single paths. We thus obtain a 6-node network model. Additionally, we assume that there are three generator nodes, where node 1 is the classical generator node, node 2 is the renewable energy node, and node 3 is the storage node that operates in discharging state. The nodes parameters and lines parameters of the Western System Coordinating Council in the simulations are given in Tables 1 and 2. From Table 2, the branch resistance value is an order of magnitude smaller ⏐ { than } the branch reactance value. ⏐ αij
⏐ P∗ = 598.9573, Γcritical = maxi̸=j ⏐ Di −
Pj∗
⏐ = 91.5, Dj ⏐ ⏐ ⏐ Γcritical γmin = arcsin Γ = 8.7872◦ and γ = maxi̸=j ⏐θi (0) − θj (0)⏐ = max{|0◦ − 9.28◦ | , . . . , |0.7275◦ − (−3.6874◦ )|} = min ◦ ◦ |9.28 − (−3.9888 )| = 13.269◦ . Since γ ∈ (γmin , π /2), according to Theorem 1, there exists the frequency synchronized solution, since Laplace matrix ( ) α L(t) is time-varying, its eigenvalues change as well, where aij = Dij cos θi − θj . Fig. 2(b) illustrates that some eigenvalues i of the −L(t) are far greater than zero after 9 s, which does not satisfy the stability condition given by Theorem 1. Fig. 2(c) We can calculate from the parameters to obtain: Γmin = n mini̸=j
(
)
Di
i
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
7
Fig. 2. Subplots (a), (b), (c) illustrate the nodes frequency variation, the eigenvalues of −L(t) and the node phase variation of the Western System Coordinating Council, respectively. Simulation results using D = [0.2, 0.1, 0.1, 0.1, 0.1, 0.1]T and P = [16.3, 8, 0.85, −1, −0.9, −0.8]T .
illustrates the phase differences divergence. Therefore, although there exists a frequency synchronization solution, the solution is unstable. Thus the system does not obtain the stable synchronization frequency as shown in Fig. 2(a). In Fig. 3 the simulation parameters are changed. The damping vector is D = [2, 1, 1, 1, 1, 1]T . The injection power P= { } α
[1.63, 0.8, 0.85, −1, −0.9, −0.8]T . The simulation results are shown in Fig. 3. In this case we have Γmin = n mini̸=j Dij = i ⏐ ∗ ⏐ ( ) P∗ ⏐ ⏐P Γ 59.986, Γcritical = maxi̸=j ⏐ Di − Dj ⏐ = 1.85, γmin = arcsin Γcritical = 1.77◦ and γ is the same as that in Fig. 2. i j min For γ ∈ (γmin , π /2), there exists a frequency synchronization solution. Moreover, −L(t) has exactly one and only
one zero eigenvalue after 3 s and others are smaller than zero, which satisfies the stability condition given by Theorem 1. Fig. 3(c) illustrates that the phase difference between arbitrary two nodes converges, that is, the synchronization frequency is stable for these parameters. Therefore, the solutions of the system exponentially converge to stable∑synchronization ∑n n ∗ frequency in Fig. 3(a). In another words, the frequency of all nodes converge to the stable value θ˙∞ = i=1 Pi / i=1 Di . The North Local Power Grid of Shaanxi Province with 58 generator nodes and 57 load nodes is used to further verify ⏐ ⏐ { } Theorem 1 in this paper. We figured out that Γmin = n mini̸=j
γmin = arcsin
(
Γcritical Γmin
)
αij Di
⏐ P∗ = 11.0764, Γcritical = maxi̸=j ⏐ Di − i
Pj∗ Dj
⏐ ⏐ = 4.3505,
⏐ ⏐ = 23.1246◦ and γ = maxi̸=j ⏐θi (0) − θj (0)⏐ = 60.2112◦ . Obviously, for γ ∈ (γmin , π/2), there exists a frequency synchronization solution. Moreover, as shown in Fig. 4(b), −L(t) has one and only one zero eigenvalue at each time point and others are smaller than zero, which satisfies the stability condition given in Theorem 1. Fig. 4(c) illustrates that the phase difference between arbitrary two nodes converges, that is the synchronization frequency is stable. Therefore, the solutions of the system converge to stable synchronization frequency as shown ( in Fig. 4(a). ⏐ ⏐ { exponentially } ) αij
⏐ P∗ = 11.0764, Γcritical = maxi̸=j ⏐ Di −
Pj∗
⏐ = 8.8901, γmin = arcsin ΓΓcritical = Dj ⏐ i min ⏐ ⏐ 53.3825◦ and γ = maxi̸=j ⏐θi (0) − θj (0)⏐ = 60(.2112◦ . According ) { } to Theorem 1, there exists the frequency synchronization ∑n solution, since Laplace matrix L(t) = diag is time-varying, its eigenvalues change as well, where j=1 aij − aij ( ) α aij = Dij cos θi − θj . Fig. 5(b) illustrates that some eigenvalues of −L(t) are greater than zero, which does not satisfy i In Fig. 5, we have Γmin = n mini̸=j
Di
the stability condition given in Theorem 1. As shown in Fig. 5(c), the phase differences are divergent. Therefore, although there does exist the frequency synchronization solution, the solution is unstable. Thus the system does not obtain the stable synchronization frequency as shown in Fig. 5(a). Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
8
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Fig. 3. Subplots (a), (b), (c) illustrate the nodes frequency variation, the eigenvalues of −L(t) and the node phase variation of Western System Coordinating Council, respectively. Simulation results using D = [2, 1, 1, 1, 1, 1]T and P = [1.63, 0.8, 0.85, −1, −0.9, −0.8]T .
6. Conclusions The models of generator and load nodes, according to rotor dynamics of the generator and the RES nodes and equipped with power-frequency droop inverter controllers, are given in this paper, where the storage nodes have two operating states. We proved that the stability of synchronization solutions depends on the eigenvalues of −L(t) given in Eq. (18), using the consensus protocol of the linear time-varying multi-agent system. With the change of time, if matrix −L(t) has exactly one zero eigenvalue and all of the nonzero eigenvalues locate in the open left half plane, then the synchronization solutions is stable. Simulation examples validate the applicability of this Theorem using the power grid of the Western System Coordinating Council and the North Local Power Grid of Shaanxi Province. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The work is supported in part by Key Program of Natural Science Foundation of Shaanxi Province (2016ZDJC-01), in part by Shaanxi Provincial Special Support Program for Science and Technology Innovation Leader, China, in part by the Shaanxi Industrial Key Project, China (2018GY-165), and in part by Key Lab Research Program from Department of Education of Shaanxi Provincial Government, China (18JS081). Appendix A The motion equation of the generator rotor is governed by J δ¨ = Ta ,
(A.1)
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
9
Fig. 4. Subplots (a), (b), (c) illustrate the nodes frequency variation, the eigenvalues of −L(t) with enlarged part to show the positive ones and the node phase variation of the North Local Power Grid of Shaanxi Province, respectively.
where δ is the mechanical angle of the shaft in radians with respect to a fixed reference, J is the moment of inertia of the rotor in kg .m2 , and Ta is the sum torque in Newton meters. J δ¨ = Tm − TL J δ¨ = Tm − (Te + Td ) .
(A.2)
The angular reference is chosen with respect to a synchronously rotating reference frame, then
δ = ωR t + θm , δ˙ = ωR + θ˙m = ωR + ωm and δ¨ = θ¨m = ω ˙ m,
(A.3) (A.4)
where θm is the mechanical angle of the shaft in radians with respect to a synchronously rotating reference frame, θ˙m is the angular velocity of the shaft, ωR is synchronously rotating angular velocity. Another form of Eq. (A.2) is obtained by multiplying its both sides with ωm , as shown by Eq. (A.5) J ωm θ¨m = Pm − PL M θ¨m = Pm − PL ,
(A.5)
where Pm = Tm ωm is the mechanical power, PL = TL ωm is the load power. The quantity J ωm is called the inertia constant 2 and is denoted by M, it is related to kinetic energy of the rotating masses Wk , where Wk = 12 J ωm . Then M can be calculated as follows M = J ωm =
2Wk
ωm
.
(A.6)
The electrical angle θe can be used to denote all of the equations above,
θ = θe =
p 2
θm ,
(A.7)
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
10
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Fig. 5. Subplots (a), (b), (c) illustrate the nodes frequency variation, the eigenvalues of −L(t) with enlarged part to show the positive ones and node phase variation of the North Local Power Grid of Shaanxi Province, respectively.
where p is the number of pole pairs of the (generator) motor, which is calculated as follows 120fR
p=
nR
,
(A.8)
where, nR is the rated shaft speed in rad/min, fR is the synchronously rotating frequency in Hz. The Eq. (A.1) can be equivalently written as
(
2J
)
θ¨ =
p
(
2J
)
p
ω˙ = Ta ( N · m).
(A.9)
Now, normalize Eq. (A.9) by dividing both sides by a constant equal to the rated torque at the rated speed SB3
TB =
ωmR
=
60SB3 2π nR
,
(A.10)
where SB3 is the rated three-phase power (VA), ωmR = 2π nR /60 is the rated angular velocity of the shaft. According to Eqs. (A.8)–(A.10), we have J π 2 n2R
(
)
900ωR SB3
ω˙ =
Ta TB
= Tau p.u.,
(A.11) (
746 W R2
Ta TB
where = Tau . The moment of inertia J = 500 unit of J. Then, by the following transformation 746 W R2 π 2 n2R
(
H≜
can be obtained by using the method in [16], and W R2 /g is the
(
)
)
550 g × 900ωR SB3 Wk =
g
)
1
2 Wk SB3
2 J ωm =
,
1 2
ω˙ = Tau p.u. ( 2)
×
746 W R 550 g
(A.12)
×
(2π nR )2 3600
(
=
746 W R
) 2
π 2 n2R
550 g × 1800
(A.13) (A.14)
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
11
we have
(
)
2Wk SB3 ωR
(
2H
)
ωR
ω˙ = Tau p.u..
(A.15)
ω˙ = Tau p.u..
(A.16)
For a classical model of the synchronous machine, recognizing that the angular speed is nearly constant, the accelerating power Pa in p.u. is nearly equal to the accelerating torque Ta [16], then the generator becomes
(
2H
)
ωR
ω˙ ∼ = Pa p.u..
(A.17)
In addition, the damping torque (damping power) includes both the mechanical and electrical damping, which is represented by Dω, then Eq. (A.17) can be rewritten as
(
2H
)
ωR
ω˙ = Pm − Pe − Pd = Pm − Pe − Dω pu,
(A.18)
where Pd = Dω is the damping power. The concrete formula of the electrical output power Pe of the generator is derived as follows. Consider a power system consisting of one machine connected to an infinite bus through a transmission line. A schematic representation of this case is shown in Fig. A.1(a), E is the internal voltage of the synchronous machine, L is the equivalent inductance of the transmission line, and V is the voltage of the infinite bus, which is used as reference. The equivalent electrical circuit of the system is given in Fig. A.1(b), E = E ̸ δ is the voltage vector of the synchronous machine and δ is the voltage phase, and V = V ̸ 0 is the voltage vector of the infinite bus, x′d is the transient reactance of the machine, V¯ t is the terminal voltage of the synchronous machine, which can be eliminated by using Y − ∆ transformation, Z¯s is the equivalent shunt impedance at the machine terminal, and Z¯TL is the series impedance of the transmission line. Another equivalent circuit of Fig. 1(b) is given in Fig. A.1(c) by a Y − ∆ transformation. Assumptions for deriving the model are the following: (1) The mechanical power input remains constant during the transient period. (2) Damping or asynchronous power is negligible. (3) The synchronous machine can be represented by a constant voltage source and its mechanical angle coincides with the electrical phase angle of the voltage. (4) If a local load is fed at the terminal voltage of the machine, it can be represented by a constant impedance to neutral point. In Fig. A.1(a) we define: Nodes 0, 1 and 2 are the reference nodes, the internal voltage node of the synchronous machine and the infinite bus node, respectively; y10 , y12 , y20 are the three admittance elements, respectively, derived by using a Y − ∆ transformation. The symbols with an arrow in the following equations are the vector variables corresponding to Fig. A.1(c). For node 1, according to Kirchhoff’s law, we obtain: I1 = (E − V) y12 + Ey10 = Ey12 − Vy12 + Ey10 = E (y12 + y10 ) − Vy12 .
(A.19)
Define: Y11 = Y11 ̸ θ11 = y12 + y10 , θ11 = arctan Y12 = Y12 ̸ θ12 = −y12 , θ12 = arctan
(
(
real(Y11 ) imag (Y11 )
real(Y12 ) imag (Y12 )
)
) (A.20)
.
Then, I1 = EY11 + VY12 = EY11 ej(δ+θ11 ) + V Y12 ejθ12 .
(A.21)
So the power input Pe at node 1 is Pe = real EI∗1
(
)
(A.22)
EI∗1 = Eejδ EY11 e−j(δ+θ11 ) + V Y12 e−jθ12 = E 2 Y11 e−jθ11 + EV Y12 e−j(θ12 −δ)
(A.23)
Pe = real EI1 = E Y11 cos θ11 + EV Y12 cos (θ12 − δ) .
(A.24)
(
(
)
∗
)
2
When the voltage angle of the infinite bus is different from zero, we can also get the corresponding Pe , given as Pe = real EI∗1 = E 2 Y11 cos θ11 + EV Y12 cos (θ12 + α − δ) .
(
)
Now, we get the compact formula of the electrical output power of the generator as Eq. (A.25). Define β = G11 ≜ Y11 cos θ11 , G11 is the self-conductance, then Eq. (A.25) is written as Pe = E 2 G11 + EV Y12 sin (δ − α + β) .
(A.25) π 2
− θ12 and (A.26)
Consequently, we can get the mathematical model of the generator node in the power system [16] as
(
2H
ωR
)
ω˙ = Pm − Pd − Pe = Pm − Dω − E 2 G11 − EV Y12 sin (δ − α + β) .
(A.27)
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
12
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Fig. A.1. One machine connected to an infinite bus through a transmission line.
In nodes power grid, the dynamic models of the generator nodes can be written as
(
2Hmi
ωR
)
⎛ ω˙ i = Pmi − Di ωi − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(δi − δj + ϕij )⎠ .
(A.28)
j=1
It can be also expressed as
⎛ Mi θ¨i = Pmi − Di θ˙i − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(θi − θj + ϕij )⎠ ,
(A.29)
j=1
where Mi θ¨i is the inertial term. When frequency synchronization analysis of Eq. (A.29) is carried out, the inertia term Mi θ¨i can be omitted due to the fact that it only affects the convergence time to synchronize the state but not its existence [16]. In addition, the resistive component is much smaller than the inductive component of the generator’s output impedance and the line impedance, then ϕij = 0. In consequence, the dynamic models of the generator nodes is given by
⎛ Di θ˙i = Pmi − ⎝Ei2 Gii +
n ∑
⎞ Ei Ej Yij sin(θi − θj + ϕij )⎠ ,
(A.30)
j=1
where, in terms of node i, Di and Mi are the damping and inertia constants, respectively, θi of node i is the angle of the shaft in radians with respect to a synchronously rotating reference, ϕij is the phase shifts caused by transfer conductance, Ei is the voltage magnitude of the synchronous machine, Yij is the admittance of the transmission line between nodes i and j, and Gii ≜ Yii cos θi is the self-conductance. Appendix B The model for the droop controlled inverter connected RES node is the following. Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
13
Fig. B.1. The equivalent circuit of two inverters operating in parallel.
Fig. C.1. Schematic diagram of local power grid in northern Shaanxi.
Power grid typically contains a bank of RES nodes equipped with active power-frequency droop inverter controllers operating in parallel. For the convenience of analysis, using two inverters in parallel as an example, Fig. B.1 is the equivalent circuit for two inverters operating in parallel, where E1 , E2 are the input voltage and I1 , I2 are the input current of inverter 1 and 2, respectively. Ri is the input impedance of inverter i, jXi is the line impedance between inverter and bus, E is the voltage of load bus, I0 , Z0 are the load current and load impedance, respectively, and δ1 , δ2 are the voltage phase difference between inverter 1, 2 and the bus, respectively. We get the inverter output active power from Fig. B.1.
Pi =
Ri Ei E cos δi − Ri E 2 Xi2
+
R2i
+
Xi E i E Xi2 + R2i
sin δi , i = 1, 2.
(B.1)
Assume that the entire inverter system is inductive, the resistive portion of the inverter output impedance and line impedance is much smaller than the inductive portion, and assume that the two inverters output impedance are equal, i.e, R1 = R2 ≈ 0, then Eq. (B.1) can be equivalently transformed as Pi =
Ei E Xi
sin δi , i = 1, 2.
(B.2)
Compared with the load impedance, inverter source output reactance is very small, so the output voltage phase angle difference δi is also very small. Then we have sin δi ≈ δi , therefore Eq. (B.2) can equivalently be given by: Pi =
Ei E Xi
δi , i = 1, 2.
(B.3)
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.
14
H.-P. Ren, Y. Gao, L. Huo et al. / Physica A xxx (xxxx) xxx
Additionally, in a practical inverter system, the variation range of output voltage amplitude is not large, which can be approximately viewed as a constant, then the output power of inverter system is only related to the phase angle difference
δi
∆Pi =
Ei E Xi
∆δi , i = 1, 2.
(B.4)
Since the phase control is achieved by adjusting the output frequency, the output frequency change of the inverter is δi , so we can change the inverter output phase by adjusting the output frequency, thereby, altering the ∆fi = ∆2πωi = 2∆ π ∆t active power output of the inverter. Define Xi npi = , (B.5) 2 π Ei E where the parameter npi is referred to as the droop coefficient, then the frequency droop control formula can be obtained fi = f0i − npi Pi ,
(B.6)
where f0i is the output frequency of the inverter without load. Appendix C The topology of the local power grid in northern Shaanxi is given in Fig. C.1. References [1] C.C. Chu, H.C. Iu, Complex networks theory for modern smart grid spplications: a survey, IEEE J. Emerg. Sel. Top. Circuits Syst. 7 (2) (2017) 177–191. [2] F.A. Rodrigues, T.K.D.M. Peron, P. Ji, J. Kurths, The Kuramoto model in complex networks, Phys. Rep. 610 (2016) 1–98. [3] G. Filatrella, A.H. Nielsen, N.F. Pedersen, Analysis of a power grid using a Kuramoto-like model, Eur. Phys. J. B 61 (4) (2008) 485–491. [4] F. Dörfler, F. Bullo, Synchronization and transient stability in power networks and non-uniform Kuramoto oscillators, SIAM J. Control Optim. 50 (3) (2012) 1616–1642. [5] R. Carareto, M.S. Baptista, C. Grebogi, Natural synchronization in power-grids with anti-correlated units, Commun. Nonlinear Sci. Numer. Simul. 18 (4) (2012) 1035–1046. [6] S. Lozano, L. Buzna, A.D. Guilera, Role of network topology in the synchronization of power systems, Eur. Phys. J. B 231 (85) (2012) 231–238. [7] L.Y. Lu, C.C. Chu, Consensus-based secondary frequency and voltage droop control of virtual synchronous generators for isolated AC micro-grids, IEEE J. Emerg. Sel. Top. Circuits Syst. 5 (3) (2015) 443–455. [8] J.W.S. Porco, F. Dörfler, F. Bullo, Synchronization and power sharing for droop-controlled inverters in islanded microgrids, Automatica 49 (9) (2013) 2603–2611. [9] J. Wang, Y.H. Liu, Y. Jiao, H.Y. Hu, Cascading dynamics in congested complex networks, Eur. Phys. J. B 67 (1) (2009) 95–100. [10] L.X. Yang, J. Jiang, Impacts of link addition and removal on synchronization of an elementary power network, Physica A 479 (2017) 99–107. [11] L. Moreau, Stability of continuous-time distributed consensus algorithms, in: IEEE Conference on Decision and Control, Vol. 4, 2004, pp. 3998-4003. [12] C. Godsil, Algebraic Graph Theory, Springer, 2001. [13] F. Dörfler, F. Bullo, Kron, Reduction of graphs with applications to electrical networks, IEEE Trans. Circuits Syst. I Regul. Pap. 60 (1) (2013) 150–163. [14] Wei. Ren, R.W. Beard, Beard Distributed Consensus in Multi-Vehicle Cooperative Control, Springer, London, 2008. [15] X.F. Sun, Y.C. Hao, Q.F. Wu, X.Q. Guo, B.C. Wang, A multifunctional and wireless droop control for distributed energy storage units in islanded AC microgrid applications, IEEE Trans. Power Electron. 32 (1) (2017) 736–751. [16] P.M. Anderson, A.A. Fouad, Power System Control and Stability, Wiley, 2003.
Please cite this article as: H.-P. Ren, Y. Gao, L. Huo et al., Frequency stability in modern power network from complex network viewpoint, Physica A (2019) 123558, https://doi.org/10.1016/j.physa.2019.123558.