Enhancing synchronization stability in a multi-area power grid

Maintaining a synchronous state of generators is of central importance to the normal operation of power grids, in which many networks are generally interconnected. In order to understand the condition under which the stability can be optimized, it is important to relate network stability with feedback control strategies as well as network structure. Here, we present a stability analysis on a multi-area power grid by relating it with several control strategies and topological design of network structure. We clarify the minimal feedback gain in the self-feedback control, and build the optimal communication network for the local and global control strategies. Finally, we consider relationship between the interconnection pattern and the synchronization stability; by optimizing the network interlinks, the obtained network shows better synchronization stability than the original network does, in particular, at a high power demand. Our analysis shows that interlinks between spatially distant nodes will improve the synchronization stability. The results seem unfeasible to be implemented in real systems but provide a potential guide for the design of stable power systems.


S.1. NETWORK REDUCTION
When analyzing the synchronization stability, the network has to be reduced to be composed only of generators; load nodes and branching points have to be removed from the network. This process is often named as Kron reduction [1]. According to Kirchhoff's law, we have I = Y V , where I and V are the vectors of currents and voltages; Y is the admittance matrix. Since the non-generator nodes are regarded as constant impedances, their injection currents are zeros. By rearranging Y such that the first n indices correspond to the generator nodes and the remaining r corresponding to the non-generator nodes, we By removing V r , the effective admittance matrix is given as

S.2. THE POWER FLOW EQUATIONS
The swing equations for the entire system composed of two subnetworks 'a' and 'b' are given byθ for i ∈ G a and j ∈ G b .
Since generators are controlled by governors, we have to consider the dynamics of generators and loads separately. The consumed electric power at load i in subnetwork 'a' is represented by P a e,i and Q a e,i , which are the active and reactive power, respectively. For each load i ∈ L, where L denotes the set of loads in the network, the power-balanced condition should be satisfied as follows: for i ∈ L a and j ∈ L b . The sets of nodes in subnetwork 'a' and subnetwork 'b' are denoted Equation (S4) contains 2n a +2n b variables. Given the admittance matrix Y ik = G ik +jB ik with j 2 = −1 and the parameters P a e,i , Q a e,i , P b e,i , and Q b e,i for i ∈ L, the voltage |V i | and the phase θ i , for ∀i ∈ L, can be calculated by the Newton-Raphson method, where the voltages and the phases of the generators are calculated with equation (S3) by the Runge-Kutta method. The voltages of generators are assumed to be constant since constant voltage magnitude is usually maintained by voltage regulators.

S.3. THE EQUILIBRIUM SOLUTION OF THE SWING EQUATIONS FOR THE TWO INTERCONNECTED SUBNETWORKS
We denote the equilibrium solution of equation (S3) as (θ a * ,i , ω a * ,i , θ b * ,j , ω b * ,j ) for i ∈ G a and j ∈ G b , and (θ a i , ω a i , θ b j , ω b j ) is the state obtained by the perturbation around the equilibrium as follows: For simplicity, we assume that the transmission line is lossless, i.e., G aa ik = G bb ik = G ab ik = G ba ik = 0. The linearized equation around the equilibrium is given bẏ The differential equations for generators j in subnetwork 'b' can be expressed in a similar way as follows:

S.4. THE STEADY-STATE STABILITY WITH SELF-FEEDBACK CONTROL
The self-feedback control matrix K and the damping matrix M are defined as The eigenvalue of the matrix K, λ K,i is given as follows: where the control strength for each subarea is assumed to be the same as γ a,i = γ a and γ b,i = γ b , respectively. In order to conduct the analysis, we further assume that λ K,i can be adjusted to a common value as λ K . By changing the scale of variables, the elements of the matrices M a and M b can also be adjusted to a common value λ M =D H .
Denote the matrix L as The synchronization stability is then determined by the following eigenvalues : Our goal is to find a matrixÃ aa such that the synchronization stability can be improved most.
The linearized swing equation for generator i in subnetwork 'a' is given by for ∀i ∈ G a .
By doing similar analysis as we did in the self-feedback control, we obtain an equation for variables X 1 and X 2 as follows:  where the form of the control matrix K can be expressed as: By diagonalizing the matrix C + K, we obtain J CK = Q −1 CK (C + K)Q CK , where Q CK is composed of the eigenvectors. The eigenvalues of J CK are ordered as 0 = λ CK,1 ≤ · · · ≤ λ CK,na+n b . By setting Z 1 = Q −1 CK X 1 and Z 2 = Q −1 CK X 2 , respectively, we have  We denote the matrix L as Under the interaction of the network topology and the communication structure, the synchronization stability is determined by the eigenvalues of L, The expression of the matrix K relies on different local strategies. For instance, if a communication network is built in subnetwork 'b', denoted byÃ bb , the local feedback control matrix K is given by where K bb is defined in a similar way as K aa , given by The general expression of the matrix K can be expressed as: If a global control of generators is implemented, a global communication matrix that connects generators in subnetwork 'a' with generators in subnetwork 'b' is established. Then, the matrixÃ ab na×n b is defined as: ab ik for i ∈ G a , k ∈ G b , where γ ab represents the control strength of generators locating in subnetworks 'a' and 'b'.
The matrix K aa na×na is a diagonal matrix composed of the element K aa ii = − γ ab H a,i k∈G bÃ ab ik for i ∈ G a . The matrixÃ ba n b ×na and the matrix K ba n b ×na can be defined in a similar way. Here, we only present the algorithm for the optimal communication network in subnetwork 'a',Ã aa . The algorithm can be applied to the optimal communication network in subnetwork 'b' and the global communication network in a similar way. Assume that there are totally m communication links that will be built in subnetwork 'a', the process of establishing the communication network is described as follows: 1. Initially, start with the subnetwork 'a', which is only composed of generators. SetÃ aa as a zero matrix and calculate Λ max of the matrix L (equation (16) The ascending order of the real parts of the eigenvalues is given as 0 =λ 1 <λ 2 ≤ · · · ≤λ n .
The largerλ 2 is, the more synchronously stable the network is.
In the following, we apply perturbation analysis to improveλ 2 by adding interlinks appropriately in the two interconnected networks. Assume that the weight of a link connecting nodes i and j is w ij > 0. When the two subnetworks 'a' and 'b' are isolated, the weighted Laplacian matrix is given by where W a and W b represent the weighted Laplacian matrices of the subnetworks 'a' and 'b', respectively. Since W is a real symmetric matrix, it has n a + n b real eigenvalues, which are ordered as 0 =λ 1 =λ 2 ≤λ 3 ≤ · · · ≤λ na+n b , where 0 is the eigenvalue with multiplication 2 due to the two isolated subnetworks with the eigenvectors all of whose components are 1.
When adding a new interlink between the two subnetworks, the new Laplacian matrix of the network W is a perturbed matrix of W , expressed as where represents the perturbed matrix; the matrix Ψ ij = w ij if nodes i and j are interconnected; the matrix D Ψ is a diagonal matrix with element (D Ψ ) ii = j w ij . For instance, by adding a new interlink, the perturbed Laplacian matrix ∆W is The second nonzero eigenvalue of W is perturbed aroundλ 2 , i.e.,λ 2 ( ) =λ 2 + ∆λ 2 + O( ), where is the coupling strength. By setting the eigenvector ofλ 2 as 1 , . . . , u na , . . . , u na+n b ), where the superscription denotes the Fiedler vector while the subscription denotes the node index. We have where < . > is the Euclidean inner product. Since ∆W is semi-definite, we have ∆λ 2 ≥ 0.
The larger ∆λ 2 is, the largerλ 2 is; hence, the more synchronizable the entire network is.
Therefore, we can add such an interlink that maximizes ∆λ 2 , that is, Initially, the two subnetworks are isolated, then 0 is the eigenvalue with multiplicity 2.
However, if we fix u (1) = (1, 1, . . . , 1), then u i − u (2) j is constant for arbitrary i ∈ G a and j ∈ G b . Therefore, we can choose the interlink with the largest value w ij (u j ) 2 . The algorithm for finding the optimal interlinks between two subnetwork is described as follows: 1. Start with two isolated subnetworks and calculate the eigenvalueλ 2 of the Laplacian matrix W of the entire network.
2. Add one interlink that connect with node i and k that maximizes w ik (u