Synchronization of Power Systems
under Stochastic Disturbances

Zhen Wang Kaihua Xi KXI@sdu.edu.cn Aijie Cheng Hai Xiang Lin André C. M. Ran Jan H. van Schuppen Chenghui Zhang School of Mathematics, Shandong University, Jinan, 250100, P. R. China Delft Institute of Applied Mathematics, Delft University of Technology, Delft, 2628 CD, The Netherlands Institute of Environmental Sciences (CML), Leiden University, Leiden 2333 CC, The Netherlands Department of Mathematics, Faculty of Science, Vrije Universiteit, Amsterdam, 1081 HV, The Netherlands Research Focus: Pure and Applied Analytics, North-West University, Potchefstroom, South Africa School of Control Science and Engineering, Shandong University, Jinan, 250061, P. R. China

Abstract

The synchronization of power generators is an important condition for the proper functioning of a power system, in which the fluctuations in frequency and the phase angle differences between the generators are sufficiently small when subjected to stochastic disturbances. Serious fluctuations can prompt desynchronization, which may lead to widespread power outages. Here, we model the stochastic disturbance by a Brownian motion process in the linearized system of the non-linear power systems and characterize the fluctuations by the variances of the frequency and the phase angle differences in the invariant probability distribution. We propose a method to calculate the variances of the frequency and the phase angle differences. For the system with uniform disturbance-damping ratio, we derive explicit formulas for the variance matrices of the frequency and the phase angle differences. It is shown that the fluctuation of the frequency at a node depends on the disturbance-damping ratio and the inertia at this node only, and the fluctuations of the phase angle differences in the lines are independent of the inertia. In particular, the synchronization stability is related to the cycle space of the network. We reveal the influences of constructing new lines and increasing capacities of lines on the fluctuations in the phase angle differences in the existing lines. The results are illustrated for the transmission system of Shandong Province of China. For the system with non-uniform disturbance-damping ratio, we further obtain bounds of the variance matrices.

keywords:

invariant probability distribution, variances, network topology, graph theory, system stability, cycle space.

^†^†thanks: This paper was not presented at any IFAC meeting. This work is financially supported by the Foundation for Innovative Research Groups of the National Natural Science Foundation of China with grant No. 61821004, the National Natural Science Foundation of China with grant No. 62103235 and the Natural Science Foundation of Shandong Province with grant No. ZR2020QF118. Corresponding author: Kaihua Xi.

, , , , , ,

1 Introduction

Power grids deliver a growing share of the energy consumed in the world and are undergoing an unprecedented revolution because of the increasing integration of intermittent power sources such as solar and wind energy and the commercialization of plug-in electric automobiles. These developments will change the structure of power sources and decrease carbon emissions dramatically, but they will also lead to new disturbances associated with fluctuations in energy production and load. These disturbances not only deteriorate the quality of the power supply but may trigger loss of synchronization, which can result in serious blackouts (Marris, 2008). This indicates the necessity to study synchronization under stochastic disturbances.

Here, we focus on the synchronization of power systems under stochastic disturbances. We explore the role of system parameters in a framework of stochastic systems that can be extended to other real complex networks with synchronization. In a synchronous state of a power system, the frequencies of the synchronous machines (e.g., rotor-generators driven by steam or gas turbines) should all be equal or close to the nominal frequency (e.g., 50 Hz or 60 Hz). Here, the frequency is the derivative of the rotational phase angle and is equal to the rotational speed of the synchronous machine in units of $rad/s$ . The synchronization stability is defined as the ability to maintain synchronization under disturbances, which is also called transient stability Kundur (1994).The parameters that determine synchronization include the power flows, inertia (Poolla et al., 2017) and damping coefficients (Motter et al., 2013; Nishikawa et al., 2015) of the synchronous machines as well as the coupling strength (Fazlyab et al., 2017) between the synchronous machines and the network topology, which can be optimized to enhance stability by load-frequency control or by constructing new power generators, virtual inertia and transmission lines. In the analysis of the existence condition of a synchronous state (Dörfler and Bullo, 2012) and the linear (Motter et al., 2013; Nishikawa et al., 2015) and nonlinear stability of that state (Zaborszky et al., 1988; Menck et al., 2013; Chiang et al., 1988), the focus is on the synchronous state, on the local convergence or on the basin of attraction. However, in practice, the state of the power system never stays at the synchronous state and is always fluctuating due to various disturbances. If both the fluctuations of the frequency and the phase angle difference are so large that the system cannot return to the synchronous state, then the synchronization is lost. Hence, the impact of the disturbances cannot be neglected and the size of the fluctuations directly characterizes the stability of the system.

Robust control methods in load frequency control may be used to improve the stability, where the disturbances are considered, see Trip et al. (2020, 2019); Xi et al. (2020). By these methods, the power generation are controlled to balance the disturbances. However, besides the power generations, the stability of the system also depends on the network topology, line capacities, inertia of the synchronous machines and so on, for which the values cannot be specified by the robust control methods. By modelling the disturbances as inputs to the associated linearized system, the fluctuations are evaluated by the $\mathcal{H}_{2}$ norm of the input-output linear system (Tegling et al., 2015; Poolla et al., 2017). However, because the $\mathcal{H}_{2}$ norm equals to the trace of a matrix (Doyle et al., 1989), which is a global metric for the synchronization stability, the fluctuations of the frequency at each node, the phase angle difference in each line and their correlation can hardly be explicitly characterized. Clearly, the nodes with serious fluctuations in the frequencies and the lines with serious fluctuations in the phase angle differences are vulnerable to disturbances. In physics, the focus is on the propagation of the disturbances (Haehne et al., 2019; Kettemann, 2016; Zhang et al., 2020; Auer et al., 2017; Zhang et al., 2019) and the network susceptibility (Manik et al., 2017). For example, the statistics of the fluctuations at the nodes, e.g., the variance of the increment of the frequency distribution, can be calculated via simulations by modelling the disturbances by either Gaussian or non-Gaussian noise (Haehne et al., 2019). With perturbations added to the system parameters, the disturbance arrival time and the vertex and edge susceptibility are estimated in (Zhang et al., 2020; Manik et al., 2017) respectively. The amplitude of perturbation responses of the nodes is used to study the emergent complex response patterns across the network (Zhang et al., 2019). By these investigations on fluctuations, intuitive insights on the impact of the system parameters, e.g., the network topology and the inertia of synchronous machines, on the spread of the disturbances are provided, which may help to develop practical guiding principles for real network design and control.

In this paper, we investigate the fluctuations of the frequency at each node and the phase angle difference in each line in a linear stochastic system. By modelling the disturbances by Gaussian noise, we use the variance in the invariant probability distribution to characterize the fluctuations and propose an efficient method for the calculation of the variance by solving a Lyapunov equation instead of statistics with a large amount of simulations. Under assumptions of uniform disturbance-damping ratio at the nodes, explicit formulas for the variances of the fluctuations in the frequencies and phase angle differences are derived, which can be used to tune the system parameters to improve the synchronization stability. With these explicit formulas, the impact of the network topology on the synchronization stability is considerably clarified.

The contribution of this paper include:

(i)

a new metric, that is the variance in the invariant probability distribution of the frequencies at the nodes and the phase angle differences in the lines, is proposed for the analysis of the synchronization stability. With this metric, the vulnerable nodes and lines can be identified effectively based on solving a matrix Lyapunov equation;
(ii)

under the assumption that the disturbance-damping ratio is uniform, we derive an explicit formula of the variance matrix of the frequency, which reveals the impact of the inertia and the disturbance-damping ratio, and an explicit formula of the variance matrix of the phase angle differences, which reveals the impact of the network topology and the disturbance-damping ratio;
(iii)

for non-uniform disturbance-damping ratio, an upper and a lower bound of the variance matrices of the frequency and the phase angle differences are deduced;
(iv)

the impact of constructing new lines and increasing the capacity of lines on the variance are investigated.

The findings of this paper provide directions for the optimization of the droop control coefficients, the placement of virtual inertia and energy storage, changes in the network topology, and changes in the capacity of lines in the power systems. The framework of this paper for the investigation of synchronization stability may be extended to other networks with stochastic disturbances and problems of synchronization.

This paper is organized as follows. The mathematical model of the power system and the problem formulation are introduced in Section 2. We propose a method to calculate the variance matrices of the invariant probability distribution of the frequency and phase angle differences in Section 3. We derive explicit formulas for the two matrices in Section 4. Based on the explicit formulas, we deduce bounds of the variance matrices for the networks with non-uniform disturbance-damping ratio in Section 5. We study the role of the network topology in Section 6 with verification using a real network in Section 7. We conclude this paper with perspectives in Section 8.

2 Models and problem formulation

The power grid can be modelled by a graph $\mathcal{G}(\mathcal{V},\mathcal{E})$ with nodes $\mathcal{V}$ and edges $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$ , where a node represents a bus and an edge $(i,j)$ represents the transmission line between nodes $i$ and $j$ . We focus on the transmission network and assume the lines are lossless. We denote the number of nodes in $\mathcal{V}$ and edges in $\mathcal{E}$ by $n$ and $m$ , respectively. The dynamics of the power system are described by the following swing equations ((Zaborszky et al., 1988; Menck et al., 2014; Chiang et al., 1988)):


$\displaystyle\dot{\delta}_{i}$	$\displaystyle=\omega_{i},$	(1a)
$\displaystyle m_{i}\dot{\omega}_{i}$	$\displaystyle=P_{i}-d_{i}\omega_{i}-\sum_{j=1}^{n}l_{ij}\sin{(\delta_{i}-\delta_{j})},$	(1b)

where $\delta_{i}$ and $\omega_{i}$ denote the phase angle and the frequency deviation of the synchronous machine at node $i$ ; $m_{i}>0$ describes the inertia of the synchronous generators; $P_{i}$ denotes power generation if $P_{i}>0$ and denotes power load otherwise; $l_{ij}=\hat{b}_{ij}V_{i}V_{j}$ is the effective susceptance, where $V_{i}$ is the voltage; $d_{i}>0$ is the damping coefficient with droop control. Since the dynamics of the voltage and the frequency can be decoupled (Simpson-Porco et al., 2016),we restrict attention to modelling only the dynamics of the frequency and assume that the voltage of each node is a constant. In practice, the voltage can be well controlled by an Automatic Voltage Regulator (Kundur, 1994). This model is often applied to study transient stability and rotor angle stability (Nishikawa and Motter, 2015; Menck et al., 2014; Dörfler and Bullo, 2012). In this paper, we focus on the networks with the following assumption.

Assumption 2.1.

Assume the network $\mathcal{G}(\mathcal{V},\mathcal{E})$ is connected.

With Assumption 2.1, we easily obtain $m\geq n-1$ . In special case of $m=n-1$ , the network is acyclic and when $m\geq n$ , there are cycles in the network.

2.1 The synchronous state

The stable region of system (1) is analysed by Chiang et al. (1988) et al and Zaborszky et al. (1988). The stability analysis of a power system makes use of the concept of the synchronous state which satisfies, for $i=1,2,\cdots,n$ ,

\displaystyle\omega_{i}(t)=\omega_{syn},~\text{and}~\delta_{i}(t)=\omega_{syn}t+\delta_{i}^{*},

where $\delta_{i}^{*}$ and the synchronized frequency $\omega_{syn}$ satisfy,


	$\displaystyle P_{i}-D_{i}\omega_{syn}-\sum_{j=1}^{n}l_{ij}\sin{\delta_{ij}^{*}}$	$\displaystyle=0~~\text{for}~~i=1,\cdots,n,$
	$\displaystyle\omega_{syn}-\frac{\sum_{i=1}^{n}{P_{i}}}{\sum_{i=1}^{n}D_{i}}$	$\displaystyle=0,$

where $\delta_{ij}^{*}=\delta_{i}^{*}-\delta_{j}^{*}$ is the phase angle difference between nodes $i$ and $j$ , which are directly connected by transmission line $(i,j)$ . The power flow in line $(i,j)$ is $l_{ij}\sin{\delta_{ij}^{*}}$ , which is determined by the load frequency control (Kundur, 1994). $\omega_{syn}$ is the deviation of the synchronized frequency from the nominal value of frequency. There are three forms of load frequency control distinguished from fast to slow time-scales, i.e., primary, secondary and tertiary frequency control. Primary control maintains the synchronous state by droop control on a small time-scale. However, this synchronized frequency may deviate from its nominal value in a medium time-scale, which leads to $\omega_{syn}\neq 0$ . Secondary control restores the synchronized frequency to its nominal value such that $\omega_{syn}=0$ on a medium time-scale. With a prediction of power demand, tertiary control calculates the operating point stabilized by primary and secondary control on a large time-scale, which concerns the security and economy of the power system. In the control design for frequency synchronization, the power input $P_{i}$ is determined in the secondary and tertiary control. Thus, it is practical to assume that the power generation and load are balanced in the study of frequency synchronization. Thus, $\sum_{i=1}^{n}{P_{i}}=0$ , which leads to $\omega_{syn}=0$ .

Due to low line capacities, the synchronous state might not exist. For the condition of the existence of the synchronous state, we refer to Dörfler and Bullo (2014). For the number of the synchronous states, see Luxemburg and Huang (1987); Baillieul and Byrnes (1982).

2.2 The linearized model

Assume that there exists a synchronous state $(\bm{\delta}^{*},\bm{0})$ for system (1), which can be linearized as

\displaystyle\left(\begin{array}[]{c}{\bm{\dot{\delta}}}\\ {\bm{\dot{\omega}}}\end{array}\right)

\displaystyle=\left(\begin{array}[]{cc}\bm{0}&\bm{I}_{n}\\ -\bm{M}^{-1}\bm{L}_{c}&-\bm{M}^{-1}\bm{D}\end{array}\right)\left(\begin{array}[]{c}{\bm{\delta}}\\ {\bm{\omega}}\end{array}\right),

(9)

where $\bm{\delta}=\text{col}{(\delta_{i})}\in\mathbb{R}^{n}$ , $\bm{I}_{n}\in\mathbb{R}^{n\times n}$ is the identity matrix, $\bm{\omega}=\text{col}{(\omega_{i})}\in\mathbb{R}^{n}$ , $\bm{M}=\text{diag}(m_{i})\in\mathbb{R}^{n\times n}$ , $\bm{D}=\text{diag}{(d_{i})}\in\mathbb{R}^{n\times n}$ , and $\bm{L}_{c}=(\widetilde{l}_{c_{ij}})\in\mathbb{R}^{n\times n}$ is the Laplacian matrix of the network with weight $l_{ij}\cos{\delta_{ij}^{*}}$ generated by $(\bm{\delta}^{*},\bm{0})$ , which satisfies

\displaystyle\widetilde{l}_{c_{ij}}=\begin{cases}~~-l_{ij}\cos{\delta_{ij}^{*}},~~i\neq j,\\ ~~\displaystyle-\sum_{k\neq i}\widetilde{l}_{c_{ik}},~~i=j.\end{cases}

(10)

Note that the state variables in (9) are the deviations of the phase angles and frequencies from the synchronous state $(\bm{\delta}^{*},\bm{0})$ . By the second Lyapunov method, the stability of $(\bm{\delta}^{*},\bm{0})$ can be determined by the sign of the real part of the eigenvalues of the system matrix of (9). The analysis of the eigenvalue of the system matrix is also called small-signal stability analysis. It has been proven that if $l_{ij}\cos{\delta_{ij}^{*}}\geq 0$ , then the system is stable at the synchronous state $(\bm{\delta}^{*},\bm{0})$ (Zaborsky et al., 1985), which leads to the security condition

\displaystyle\bm{\Theta}=\big{\{}\bm{\delta}\in\mathbb{R}^{n}\big{|}~|\delta_{ij}|<\frac{\pi}{2},\forall(i,j)\in\mathcal{E})\big{\}}.

(11)

It has been proven by (Skar, 1980) that for the power network with a general network topology, the synchronous state in this security range is unique. For the identification of the subset of the $n$ -torus where there exists a synchronous state, we refer to (Jafarpour et al., 2022).

2.3 Problem formulation

In real networks, the state of the power system always fluctuates around the synchronous state due to various disturbances. If the fluctuations are very large, the state may exit the stability region of the synchronous state and lead to instability of the system. A sign of desynchronization is that both the fluctuations of the frequency and the phase angle difference are so large that the system cannot return to the synchronous state. Many factors influence the fluctuations, which include the parameters of the transmission lines, the synchronous machines, the network topology and the disturbances. The source of the disturbances are also various, e.g., the renewable power generation, fault of the devices in the network, etc. We focus on the following problem in this paper.

Problem 2.2.

How do the fluctuations of the frequency and the phase angle differences depend on the parameters of the system and the disturbances?

The solution of this problem provides insights for suppressing the fluctuations by scientific parameter assignments. The choice of a model for the fluctuations in a power system should be based on the criteria that the model is realistic and that the subsequent analysis is not too complex.

A realistic model of the actual disturbances affecting a power system at each node requires an extensive system identification procedure, including the collection of a large amount of data on the fluctuations of the power system. The disturbances come from both the loads and the various power sources, such as wind parks and photovoltaic units. It has been shown that the probability distributions of disturbances are not Gaussian in several real power grids, e.g., grids in North America and Europe, which leads to non-Gaussian distribution of the frequency and is crucial to induce desynchronization in the system, see Haehne et al. (2019); Schäfer et al. (2018); Schmietendorf et al. (2017); Wolff et al. (2019); Xie et al. (2011) etc.. A model could then be a nonlinear stochastic differential equation of the power system driven by either Brownian motion or another process with independent increments. However, the performance evaluation of such nonlinear stochastic system requires requires either the numerical approximation of the solution of a partial differential equation (Wang and Crow, 2013) or a large amount of simulations for the statistics of the frequencies (Haehne et al., 2019). This model is too complicated to obtain an analytic probability distribution of the state of the power system consisting of a large number of synchronous machines.

An alternative to the modelling approach described above is to formulate a deterministic linear system obtained by linearization of a nonlinear power system at a synchronous state. The deterministic linear system is then transformed into a linear stochastic differential equation driven by Brownian motion. Such models are often used in control engineering and in mathematical finance, and these models are regarded as reasonable approximations of realistic models. Moreover, these models have a low algebraic complexity. It is well known that for a linear stochastic differential equation with a system matrix that is Hurwitz, there exists an invariant probability distribution of the state that is a Gaussian probability distribution characterized by the mean value and the variance of the state (Kwakernaak and Sivan, 1972, Theorem 1.53)(Karatzas and Shreve, 1988, Theorem 6.7). For power systems, the fluctuations are described by the variance matrices in the invariant probability distribution of the associated linear stochastic system. The dependence on the system parameters is indicated. The complexity of the performance of this model is manageable. Though the analysis of the stochastic linearized system is valid only for comparatively small disturbances, it still provides intuitive insights on the stability of the power system.

When subjected to disturbances, the state of the power system deviates from the synchronous state. Hence, we study the deviation of the frequency and the phase angle difference from the synchronous state, which is the state of the linearized system of the nonlinear power system. We model the disturbances by a Brownian motion process, which is then an input to a linear system, and study the stochastic system


$\displaystyle\text{d}\bm{\delta}(t)$	$\displaystyle=\bm{\omega}(t)\text{d}t,$	(12a)
$\displaystyle\text{d}\bm{\omega}(t)$	$\displaystyle=-\bm{M}^{-1}\big{(}\bm{L}_{c}\bm{\delta}(t)+\bm{D}\bm{\omega}(t)\big{)}\text{d}t+\bm{M}^{-1}\widetilde{\bm{B}}\text{d}\bm{\mu}(t),$	(12b)

with the state variable, system matrix and input matrix,

\displaystyle\bm{x}(t)=\begin{bmatrix}\bm{\delta}(t)\\ \bm{\omega}(t)\end{bmatrix},~~\bm{A}=\begin{bmatrix}\bm{0}&\bm{I}_{n}\\ -\bm{M}^{-1}\bm{L}_{c}&-\bm{M}^{-1}\bm{D}\end{bmatrix},~~\bm{B}

\displaystyle=\begin{bmatrix}\bm{0}\\ \bm{M}^{-1}\widetilde{\bm{B}}\end{bmatrix},

(13)

where the notations $\bm{\delta}(t),\bm{\omega}(t)$ , $\bm{M},\bm{D},\bm{L}_{c}$ are defined as for (9), $\widetilde{\bm{B}}=\text{diag}(b_{i})\in\mathbb{R}^{n\times n}$ where $b_{i}\in\mathbb{R}$ and $b_{i}^{2}$ is used to characterize the strength of the disturbance; $\bm{\mu}(t)=\text{col}(\mu_{i}(t))\in\mathbb{R}^{n}$ is a vector of $n$ independent scalar Brownian motion processes $\mu_{i}$ , which are also all independent of the initial state $\bm{x}(0)$ . A Brownian motion process has increments with a Gaussian probability distribution. Here, we refer to $l_{ij}$ as the line capacity of line $e_{k}$ , which is also called the coupling strength between generators, and refer to $l_{c_{ij}}=l_{ij}\cos{\delta_{ij}^{*}}$ as the weight of line $e_{k}$ . It is obvious that the weights of the lines are determined by the line capacity and the power flows at the synchronous state.

In the model (12), the disturbances denoted by $\mu_{i}(t)$ at node $i$ are assumed to be independent, which is reasonable because the locations of the power generators, including renewable power generators, are usually far from each other. Because the system (12) is linear, at any time, the probability distribution of the state is Gaussian. To address Problem 2.2, we focus on the variance matrices of the frequency and of the phase angle difference in the invariant probability distribution of the linear stochastic system, which reflect the dependence of the fluctuations of the frequency and the phase angle difference on the system parameters. When considering the variance matrix in the invariant probability distribution, we set the output matrix so that

\displaystyle\bm{y}=\bm{C}\bm{x},~~\bm{y}=\begin{bmatrix}\bm{y}_{\delta}\\ \bm{y}_{\omega}\end{bmatrix},~~\bm{C}=\begin{bmatrix}\widetilde{\bm{C}}^{\bm{\top}}&\bm{0}\\ \bm{0}&\bm{I}_{n}\end{bmatrix}\in\mathbb{R}^{(m+n)\times 2n}.

(14)

The $m$ elements in $\bm{y}_{\delta}$ are the phase angle differences in the $m$ lines, and the $n$ elements in $\bm{y}_{\omega}$ are the frequencies at the $n$ nodes. The matrix $\widetilde{\bm{C}}=(C_{ik})\in\mathbb{R}^{n\times m}$ is the incidence matrix of the network, which is defined as

\displaystyle C_{ik}=\begin{cases}+1,&\text{if node $i$ is the begin of line $e_{k}$},\\ -1,&\text{if node $i$ is the end of line $e_{k}$},\\ 0,&\text{otherwise},\end{cases}

where the direction of line $e_{k}$ is specified arbitrarily without influence on the study below. By the complex network theory (Biggs, 1993), the incidence matrix $\widetilde{\bm{C}}$ satisfies

\displaystyle\widetilde{\bm{C}}\bm{R}\widetilde{\bm{C}}^{\bm{\top}}=\bm{L}_{c},

(15)

where $\bm{R}=\text{diag}(R_{k})\in\mathbb{R}^{m\times m}$ is defined such that $R_{k}=l_{c_{ij}}$ is the weight of line $e_{k}$ connecting nodes $i$ and $j$ .

Because $\bm{x}(t)$ is the deviation of the frequency and phase angle difference from the synchronous state, it is natural to assume that $\bm{x}(0)\in G(\bm{0},\bm{Q}_{x_{0}})$ where $\bm{Q}_{x_{0}}\in\mathbb{R}^{2n\times 2n}$ . Problem 2.2 requires the calculation of the invariant probability distribution of the deviations of the frequencies and of the phase angle differences, and requires an analysis of how this distribution depends on the parameters of the power system in particular on the intensities of the stochastic disturbances. It will be shown in Theorem 3.2 that the variance matrix in the invariant probability distribution is independent of the initial distribution. Below we restrict attention to the computation of the invariant probability distribution of the state of the linear stochastic power system. From that distribution, the variances of the outputs can be computed.

3 Derivation of the variance matrices

We denote the variance matrix of the frequencies and the phase angle differences at the invariant probability distribution by

\displaystyle\bm{Q}=\begin{bmatrix}\bm{Q}_{\delta}&\bm{Q}_{\delta\omega}^{\bm{\top}}\\ \bm{Q}_{\delta\omega}&\bm{Q}_{\omega}\end{bmatrix}\in\mathbb{R}^{(m+n)\times(m+n)},

where $\bm{Q}_{\delta}\in\mathbb{R}^{m\times m}$ denotes the variance matrix of the phase angle differences, $\bm{Q}_{\omega}\in\mathbb{R}^{n\times n}$ denotes the variance matrix of the frequencies, and $\bm{Q}_{\delta\omega}\in\mathbb{R}^{n\times m}$ denotes the covariance of the phase angle differences and the frequencies. Based on the theory of linear stochastic Gaussian systems, $\bm{Q}$ is derived by solving a Lyapunov equation, as presented in Definition A.1 in the Appendix. However, for a linear stochastic power system, the system matrix $\bm{A}$ is not Hurwitz. This is due to the singularity of the Laplacian matrix $\bm{L}_{c}$ , which has a zero eigenvalue. Therefore, the variance matrix $\bm{Q}$ cannot be calculated directly from the corresponding Lyapunov equation. A coordinate transformation is required. Before introducing the transformation, we present a lemma for the symmetrizable matrix $\bm{M}^{-1}\bm{L}_{c}$ (Xi et al., 2020, Appendix).

Lemma 3.1.

Consider the Laplacian matrix $\bm{L}_{c}$ and the positive-definite diagonal matrix $\bm{M}^{-1}$ in system (12). The matrix $\bm{L}_{c}$ has a zero eigenvalue with eigenvector $\bm{1}_{n}\in\mathbb{R}^{n}$ which is a vector with all its elements equal to one and there exists an orthogonal matrix $\bm{U}\in\mathbb{R}^{n\times n}$ such that

\displaystyle\bm{U}^{\bm{\top}}\bm{M}^{-1/2}\bm{L}_{c}\bm{M}^{-1/2}\bm{U}=\bm{\Lambda}_{n},

(16)

where $\bm{\Lambda}_{n}=\text{diag}(\lambda_{i})\in\mathbb{R}^{n\times n}$ with $0=\lambda_{1}<\lambda_{2}\cdots\lambda_{n}$ being the eigenvalues of the matrix $\bm{M}^{-1/2}\bm{L}_{c}\bm{M}^{-1/2}$ , $\bm{U}=\begin{bmatrix}\bm{u}_{1}&\bm{u}_{2}&\cdots&\bm{u}_{n}\end{bmatrix}$ with $\bm{u}_{i}\in\mathbb{R}^{n}$ being the eigenvector corresponding to $\lambda_{i}$ for $i=1,\cdots,n$ . In addition, $\bm{u}_{1}=\sigma\bm{M}^{1/2}\bm{1}_{n}$ where $\sigma$ is a constant.

Based on Lemma 3.1, we transform the coordinates of $(\bm{\delta},\bm{\omega})$ into the eigen-space as follows. Let $\bm{x}_{1}=(\bm{M}^{-1/2}\bm{U})^{-1}\bm{\delta}$ , $\bm{x}_{2}=(\bm{M}^{-1/2}\bm{U})^{-1}\bm{\omega}$ and insert (16) into (12), we derive


$\displaystyle\hskip-10.0pt\text{d}\bm{x}_{1}$	$\displaystyle=\bm{x}_{2}\text{d}t,$	(17a)
$\displaystyle\hskip-10.0pt\text{d}\bm{x}_{2}$	$\displaystyle=-\big{(}\bm{\Lambda}_{n}\bm{x}_{1}+\bm{U}^{\bm{\top}}\bm{M}^{-1}\bm{D}\bm{U}\bm{x}_{2}\big{)}\text{d}t+\bm{U}^{\bm{\top}}\bm{M}^{-1/2}\widetilde{\bm{B}}\text{d}{\bm{\mu}(t)},$	(17b)

with the state variable, system matrix and input matrix becoming

\displaystyle\bm{x}_{e}=\begin{bmatrix}\bm{x}_{1}\\ \bm{x}_{2}\end{bmatrix},~\bm{A}_{e}=\begin{bmatrix}\bm{0}&\bm{I}_{n}\\ -\bm{\Lambda}_{n}&-\bm{U}^{\bm{\top}}\bm{M}^{-1}\bm{D}\bm{U}\end{bmatrix},\bm{B}_{e}=\begin{bmatrix}\bm{0}\\ \bm{U}^{\bm{\top}}\bm{M}^{-1/2}\widetilde{\bm{B}}\end{bmatrix},

(18)

and initial distribution $\bm{x}_{e}(0)\in G(\bm{0},\bm{Q}_{x_{e_{0}}})$ such that

	$\displaystyle\bm{Q}_{x_{e_{0}}}$	$\displaystyle=\bm{T}\bm{Q}_{x_{0}}\bm{T}^{\bm{\top}}\in\mathbb{R}^{2n\times 2n},$
	$\displaystyle\bm{T}$	$\displaystyle=\begin{bmatrix}(\bm{M}^{-1/2}\bm{U})^{-1}&\bm{0}\\ \bm{0}&(\bm{M}^{-1/2}\bm{U})^{-1}\end{bmatrix}\in\mathbb{R}^{2n\times 2n}.$

The output (14) becomes

\displaystyle\bm{y}=\bm{C}_{e}\bm{x}_{e},~~\bm{C}_{e}=\begin{bmatrix}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\bm{U}&\bm{0}\\ \bm{0}&\bm{M}^{-1/2}\bm{U}\end{bmatrix}\in\mathbb{R}^{(m+n)\times 2n}.

(19)

Because $\widetilde{\bm{C}}$ is an incidence matrix of the network, it satisfies $\widetilde{\bm{C}}^{\bm{\top}}\bm{1}_{n}=\bm{0}$ . Thus, $\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\bm{u}_{1}=\bm{0}$ since $\bm{u}_{1}=\sigma\bm{M}^{1/2}\bm{1}_{n}$ , which leads to

\displaystyle\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\bm{U}=\begin{bmatrix}\bm{0}&\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\bm{u}_{2}&\cdots&\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\bm{u}_{n}\end{bmatrix},

where the entries in the first column are all zero. So the entries in the first column of $\bm{C}_{e}$ are all zero. Because the diagonal matrix $\bm{\Lambda}_{n}$ has a zero entry at position $(1,1)$ , the entries of the first column of $\bm{A}_{e}$ are also all zero. In addition, the entries of the first row of $\bm{B}_{e}$ are all zero. Hence, we decompose the system matrix $\bm{A}_{e}$ , the input matrix $\bm{B}_{e}$ , and the output matrix $\bm{C}_{e}$ into

\displaystyle\bm{A}_{e}=\begin{bmatrix}0&\bm{A}_{12}\\ \bm{0}&\bm{A}_{2}\end{bmatrix},~~\bm{B}_{e}=\begin{bmatrix}\bm{0}\\ \bm{B}_{2}\end{bmatrix},~~\bm{C}_{e}=\begin{bmatrix}\bm{0}&\bm{C}_{2}\end{bmatrix},

(20)

where $\bm{A}_{12}\in\mathbb{R}^{1\times(2n-1)}$ and $\bm{A}_{2}\in\mathbb{R}^{(2n-1)\times(2n-1)}$ , $\bm{B}_{2}\in\mathbb{R}^{(2n-1)\times n}$ and $\bm{C}_{2}$ is the matrix obtained by removing the first column of the output matrix in (19) so that

\displaystyle\bm{C}_{2}=\begin{bmatrix}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\widehat{\bm{U}}&\bm{0}\\ \bm{0}&\bm{M}^{-1/2}\bm{U}\end{bmatrix}\in\mathbb{R}^{(m+n)\times(2n-1)},

(21)

with $\widehat{\bm{U}}=\begin{bmatrix}\bm{u}_{2}&\bm{u}_{3}&\cdots&\bm{u}_{n}\end{bmatrix}\in\mathbb{R}^{n\times(n-1)}$ . According to these decompositions, the matrix $\bm{Q}_{x_{e_{0}}}$ is further rewritten as

\displaystyle\bm{Q}_{x_{e_{0}}}=\begin{bmatrix}Q_{e_{1}}&\bm{Q}_{e_{12}}^{\bm{\top}}\\ \bm{Q}_{e_{12}}&\bm{Q}_{e_{2}}\end{bmatrix},

(22)

where $Q_{e_{1}}\in\mathbb{R},~\bm{Q}_{e_{12}}\in\mathbb{R}^{2n-1},~\bm{Q}_{e_{2}}\in\mathbb{R}^{(2n-1)\times(2n-1)}$ .

In (20), $\bm{A}_{2}$ is obtained from $\bm{A}_{e}$ by removing the first column and the first row and $\bm{B}_{2}$ is obtained from $\bm{B}_{e}$ by removing the first row. Since the eigenvalues of $\bm{A}_{e}$ all have non-positive real parts and $\text{rank}(\bm{A}_{e})=2n-1$ , $\bm{A}_{2}$ is Hurwitz. With (18) and (20), $\bm{A}_{2},$ and $\bm{B}_{2}$ are further written into block matrices,

\displaystyle\bm{A}_{2}=\begin{bmatrix}\bm{0}&\bm{A}_{22}\\ \bm{A}_{23}&\bm{A}_{24}\end{bmatrix},~~\bm{B}_{2}=\begin{bmatrix}\bm{0}\\ \bm{B}_{22}\end{bmatrix},

(23)

where


	$\displaystyle\bm{A}_{22}=\begin{bmatrix}\bm{0}&\bm{I}_{n-1}\end{bmatrix}\in\mathbb{R}^{(n-1)\times n},~~\bm{A}_{23}^{\bm{\top}}=\begin{bmatrix}\bm{0}&-\bm{\Lambda}_{n-1}\end{bmatrix}\in\mathbb{R}^{(n-1)\times n},$		(24a)
	$\displaystyle\bm{A}_{24}=-\bm{U}^{\bm{\top}}\bm{M}^{-1}\bm{D}\bm{U}\in\mathbb{R}^{n\times n},~~\bm{B}_{22}=\bm{U}^{\bm{\top}}\bm{M}^{-1/2}\widetilde{\bm{B}}\in\mathbb{R}^{n\times n}.$		(24b)

Here, $\bm{\Lambda}_{n-1}=\text{diag}(\lambda_{i},i=2,\cdots,n)\in\mathbb{R}^{(n-1)\times(n-1)}$ is obtained by removing the first column and the first row of the diagonal matrix $\bm{\Lambda}_{n}$ . With the above notations, for the variance matrix of the output of the system (12), we have the following theorem.

Theorem 3.2.

The variance matrix $\bm{Q}$ of the output $\bm{y}$ of the system (12) in the invariant probability distribution satisfies

\displaystyle\bm{Q}=\bm{C}_{2}\bm{Q}_{x}\bm{C}_{2}^{\bm{\top}},

(25)

where $\bm{C}_{2}$ is defined in (21), $\bm{Q}_{x}\in\mathbb{R}^{(2n-1)\times(2n-1)}$ is the unique solution of the following Lyapunov equation

\displaystyle\bm{A}_{2}\bm{Q}_{x}+\bm{Q}_{x}\bm{A}_{2}^{\bm{\top}}+\bm{B}_{2}\bm{B}_{2}^{\bm{\top}}=\bm{0},

(26)

where $\bm{A}_{2},~~\bm{B}_{2}$ are defined in (23) with blocks in (24).

Proof: We decompose the state variable $\bm{x}_{e}=(x_{e_{1}},\bm{x}_{e_{2}}^{\bm{\top}})^{\bm{\top}}$ with $x_{e_{1}}\in\mathbb{R}$ , $\bm{x}_{e_{2}}\in\mathbb{R}^{2n-1}$ . From (17-18) and the decomposition of matrices in (20), we obtain the stochastic process

\displaystyle\text{d}\bm{x}_{e_{2}}(t)=\bm{A}_{2}\bm{x}_{e_{2}}(t)\text{d}t+\bm{B}_{2}\text{d}\bm{\mu}(t),

(27)

where $\bm{A}_{2}$ is Hurwitz. From (20), it is seen that the entries in the first column of $\bm{C}_{e}$ are all zero. Thus, the output $\bm{y}(t)$ satisfies

\displaystyle\bm{y}(t)=\bm{C}_{e}\bm{x}_{e}(t)=\bm{C}_{2}\bm{x}_{e_{2}}(t),

(28)

From (22), we obtain the initial value of $\bm{x}_{e_{2}}(0)\in G(\bm{0},\bm{Q}_{e_{2}})$ . Consider the stochastic process (27) with output in (28). Following Definition A.1 in the Appendix, the variance of the output $\bm{y}(t)$ is

\displaystyle\bm{Q}_{y}(t)=\bm{C}_{2}\text{e}^{\bm{A}_{2}t}\bm{Q}_{e_{2}}\text{e}^{\bm{A}_{2}^{\bm{\top}}t}\bm{C}_{2}^{\bm{\top}}+\int_{0}^{t}{\bm{C}_{2}\text{e}^{\bm{A}_{2}\tau}\bm{B}_{2}\bm{B}_{2}^{\bm{\top}}\text{e}^{\bm{A}_{2}^{\bm{\top}}\tau}\bm{C}_{2}^{\bm{\top}}}\text{d}\tau.

With the Hurwitz condition of $\bm{A}_{2}$ , we obtain the variance matrix of $\bm{y}(t)$ ,

\displaystyle\bm{Q}=\lim_{t\rightarrow}\bm{Q}_{y}(t)=\int_{0}^{+\infty}{\bm{C}_{2}\text{e}^{\bm{A}_{2}\tau}\bm{B}_{2}\bm{B}_{2}^{\bm{\top}}\text{e}^{\bm{A}_{2}^{\bm{\top}}\tau}\bm{C}_{2}^{\bm{\top}}}\text{d}\tau,

which can be solved from (25) with

\displaystyle\bm{Q}_{x}=\int_{0}^{+\infty}{\text{e}^{\bm{A}_{2}\tau}\bm{B}_{2}\bm{B}_{2}^{\bm{\top}}\text{e}^{\bm{A}_{2}^{\bm{\top}}\tau}}\text{d}\tau,

which is the Controllability Gramian of the pair $(\bm{A}_{2},\bm{B}_{2})$ and is the unique solution of the Lyapunov equation (26). $\square$

It is seen that the invariant matrix $\bm{Q}$ is independent of the initial distribution of the original process $\bm{x}(t)$ defined in (12-13). With Theorem 3.2 and the formulation of $\bm{A}_{2}$ and $\bm{B}_{2}$ in (23-24), the variance matrix $\bm{Q}_{x}$ can be obtained by solving the Lyapunov equation using Matlab and the variance matrix $\bm{Q}$ can be further calculated from (25). Clearly, the larger the variances, the more serious the fluctuations in the nodes and lines will be. Thus, from the diagonal elements of $\bm{Q}$ , the vulnerable nodes and lines with large variances can be identified.

Remark 3.3.

If $\bm{Q}_{\delta}$ is needed only, the output is set for the system (17) as

\displaystyle\bm{y}=\bm{C}_{e}\bm{x}_{e},~~\bm{C}_{e}=\begin{bmatrix}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\bm{U}&\bm{0}\end{bmatrix}\in\mathbb{R}^{m\times 2n}.

By removing the first column of $\bm{C}_{e}$ , we obtain

\displaystyle\bm{C}_{2}=\begin{bmatrix}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\widehat{\bm{U}}&\bm{0}\end{bmatrix}\in\mathbb{R}^{m\times(2n-1)}

(29)

for the calculation of $\bm{Q}_{\delta}$ by (25).

If $\bm{Q}_{\omega}$ is needed only, the output is set for the system (17) as

\displaystyle\bm{y}=\bm{C}_{e}\bm{x}_{e},~~\bm{C}_{e}=\begin{bmatrix}\bm{0}&\bm{M}^{-1/2}\bm{U}\end{bmatrix}\in\mathbb{R}^{n\times 2n}.

By removing the first column of $\bm{C}_{e}$ , we obtain

\displaystyle\bm{C}_{2}=\begin{bmatrix}\bm{0}&\bm{M}^{-1/2}\bm{U}\end{bmatrix}\in\mathbb{R}^{n\times(2n-1)}

(30)

for the calculation of $\bm{Q}_{\omega}$ by (25).

The variance of the frequency at a node can also be calculated via the $\mathcal{H}_{2}$ norm of the input-output system with the output being the frequency at this node. However, when considering the variances of the frequencies at all the nodes, $n$ Lyapunov equations need to be solved. Similarly, when considering the variances of the phase angle differences, the solution of $m$ Lyapunov equations are required. These computations have a high computational complexity. Furthermore, the correlation of the outputs cannot be derived in this way.

4 Explicit formulas of the variance matrices for networks with uniform disturbance-damping ratio

Based on the following assumption, we derive the explicit formula of the solution $\bm{Q}$ .

Assumption 4.1.

Consider the stochastic system (12). Assume that the uniform disturbance-damping ratio holds, in which there exists a strictly positive number $\eta\in(0,+\infty)$ such that for all nodes $i\in\mathcal{V}$ , $b_{i}^{2}/d_{i}=\eta$ .

In practice, in order to achieve fair power sharing, the drooping coefficients $d_{i}$ are often scheduled proportionally to the rating of the power source. Thus, it is reasonable to expect that the strength of the disturbance, that is characterized by $b_{i}^{2}$ , is proportional to the rating of the power source. On contrast to this assumption, one says that the non-uniform disturbance-damping ratio holds in the complementary case, or, equivalently, if there exist $i,~j\in\mathcal{V}$ with $i\neq j$ such that $b_{i}^{2}/d_{i}\neq b_{j}^{2}/d_{j}$ .

To compute the variance matrix $\bm{Q}$ one has to first compute the variance matrix $\bm{Q}_{x}$ as stated next.

Lemma 4.2.

We write the matrix $\bm{Q}_{x}$ defined in Theorem 3.2 into a block matrix,

\displaystyle\bm{Q}_{x}=\begin{bmatrix}\bm{Q}_{1}&\bm{Q}_{2}\\ \bm{Q}_{2}^{\bm{\top}}&\bm{Q}_{3}\end{bmatrix},

where $\bm{Q}_{1}\in\mathbb{R}^{(n-1)\times(n-1)}$ , $\bm{Q}_{2}\in\mathbb{R}^{(n-1)\times n}$ and $\bm{Q}_{3}\in\mathbb{R}^{n\times n}$ . If Assumption 4.1 holds and $\bm{Q}_{x}$ satisfies the Lyapunov equation (26), then

\displaystyle\bm{Q}_{1}=\frac{1}{2}\eta\bm{\Lambda}_{n-1}^{-1},~~\bm{Q}_{2}=\bm{0},~~\bm{Q}_{3}

\displaystyle=\frac{1}{2}\eta\bm{I}_{n},

(31)

where $\bm{\Lambda}_{n-1}$ is obtained from the matrix $\bm{\Lambda}_{n}$ by removing the first column and the first row as in (24).

Proof: With the block matrices $\bm{A}_{2}$ and $\bm{B}_{2}$ in (23) and the corresponding blocks $\bm{A}_{22}$ , $\bm{A}_{23}$ , $\bm{A}_{24}$ and $\bm{B}_{22}$ in (24), we derive from the Lyapunov equation (26) that

	$\displaystyle\begin{bmatrix}\bm{0}&\bm{A}_{22}\\ \bm{A}_{23}&\bm{A}_{24}\end{bmatrix}\begin{bmatrix}\bm{Q}_{1}&\bm{Q}_{2}\\ \bm{Q}_{2}^{\bm{\top}}&\bm{Q}_{3}\end{bmatrix}$	$\displaystyle+\begin{bmatrix}\bm{Q}_{1}&\bm{Q}_{2}\\ \bm{Q}_{2}^{\bm{\top}}&\bm{Q}_{3}\end{bmatrix}\begin{bmatrix}\bm{0}&\bm{A}_{22}\\ \bm{A}_{23}&\bm{A}_{24}\end{bmatrix}^{\bm{\top}}$
		$\displaystyle+\begin{bmatrix}\bm{0}\\ \bm{B}_{22}\end{bmatrix}\begin{bmatrix}\bm{0}&\bm{B}_{22}^{\bm{\top}}\end{bmatrix}=\bm{0}$

which yields


$\displaystyle\bm{Q}_{2}\bm{A}_{22}^{\bm{\top}}+\bm{A}_{22}\bm{Q}_{2}^{\bm{\top}}$	$\displaystyle=\bm{0},$	(32a)
$\displaystyle\bm{Q}_{1}\bm{A}_{23}^{\bm{\top}}+\bm{Q}_{2}\bm{A}_{24}^{\bm{\top}}+\bm{A}_{22}\bm{Q}_{3}$	$\displaystyle=\bm{0},$	(32b)
$\displaystyle\bm{Q}_{2}^{\bm{\top}}\bm{A}_{23}^{\bm{\top}}+\bm{Q}_{3}\bm{A}_{24}^{\bm{\top}}+\bm{A}_{23}\bm{Q}_{2}+\bm{A}_{24}\bm{Q}_{3}$	$\displaystyle=-\bm{B}_{22}\bm{B}_{22}^{\bm{\top}}.$	(32c)

The idea to solve the above equations is as follows. We first assume $\bm{Q}_{2}=\bm{0}$ , then solve $\bm{Q}_{3}$ and $\bm{Q}_{1}$ from (32c) and (32b) respectively, finally we check whether these three matrices satisfy all the equations in (32). If that is true, then from the uniqueness of the solution of (26), we have obtained the solution $\bm{Q}_{x}$ for (26). From (32c) with the formula for $\bm{A}_{24}$ and $\bm{B}_{22}$ in (24) and $\bm{Q}_{2}=\bm{0}$ , we derive

\displaystyle\bm{Q}_{3}\bm{U}^{\bm{\top}}\bm{M}^{-1}\bm{D}\bm{U}+\bm{U}^{\bm{\top}}\bm{M}^{-1}\bm{D}\bm{U}\bm{Q}_{3}=\bm{U}^{\bm{\top}}\bm{M}^{-1/2}\widetilde{\bm{B}}\widetilde{\bm{B}}^{\bm{\top}}\bm{M}^{-1/2}\bm{U}

which has a unique solution

\displaystyle\bm{Q}_{3}=\frac{1}{2}\bm{U}^{\bm{\top}}\bm{D}^{-1}\widetilde{\bm{B}}^{2}\bm{U},

where the fact that $\bm{M},~~\bm{D},~~\widetilde{\bm{B}}$ are diagonal matrices and $\widetilde{\bm{B}}^{2}=\widetilde{\bm{B}}\widetilde{\bm{B}}^{\bm{\top}}$ are used. It is obvious that the diagonal entries of $\bm{D}^{-1}\widetilde{\bm{B}}^{2}$ are $b_{i}^{2}/d_{i}=\eta$ for $i=1,\cdots,n$ , which yields $\bm{D}^{-1}\widetilde{\bm{B}}^{2}=\eta\bm{I}_{n}$ . Thus $\bm{Q}_{3}=\frac{1}{2}\eta\bm{I}_{n}$ . From (32b) with the formulas for $\bm{A}_{23}$ and $\bm{A}_{22}$ in (24) and $\bm{Q}_{2}=\bm{0}$ , we derive

\displaystyle\bm{Q}_{1}\begin{bmatrix}\bm{0}&-\bm{\Lambda}_{n-1}\end{bmatrix}+\frac{1}{2}\eta\begin{bmatrix}\bm{0}&\bm{I}_{n-1}\end{bmatrix}\bm{I}_{n}=\bm{0},

which leads to

\displaystyle-\bm{Q}_{1}\bm{\Lambda}_{n-1}+\frac{1}{2}\eta\bm{I}_{n-1}=\bm{0}.

Thus, $\bm{Q}_{1}=\frac{1}{2}\eta\bm{\Lambda}^{-1}_{n-1}$ . In conclusion, by assuming $\bm{Q}_{2}=\bm{0}$ , we have obtained the explicit formulas for $\bm{Q}_{1}$ and $\bm{Q}_{3}$ as presented in (31). Furthermore, it can be verified that $\bm{Q}_{1},\bm{Q}_{2}$ and $\bm{Q}_{3}$ satisfy (32) which is equivalent to the Lyapunov equation (26). $\square$

By Lemma 4.2, we derive the independence of the stochastic process of the frequency to the phase angle differences in the lines. In addition, an explicit formula for the variance matrix $\bm{Q}_{\omega}$ of the frequencies at the nodes is deduced.

Theorem 4.3.

Consider the system (12) with Assumption 4.1.

(i)

The invariant probability distributions of the frequencies and of the phase angle differences are independent, i.e., $\bm{Q}_{\delta\omega}=\bm{0}$ .
(ii)

The variance matrix of the frequencies is

$\displaystyle\bm{Q}_{\omega}=\frac{1}{2}\eta\bm{M}^{-1}.$ (33)

Proof: (i) We take $\bm{C}_{2}$ in (21) as the output matrix for the system (17). By Theorem 3.2, we obtain that the variance matrix $\bm{Q}$ satisfies

	$\displaystyle\bm{Q}$	$\displaystyle=\bm{C}_{2}\bm{Q}_{x}\bm{C}_{2}^{\bm{\top}}$
		$\displaystyle=\begin{bmatrix}\frac{\eta}{2}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\widehat{\bm{U}}\bm{\Lambda}_{n-1}^{-1}\widehat{\bm{U}}^{\bm{\top}}\bm{M}^{-1/2}\widetilde{\bm{C}}&\bm{0}\\ \bm{0}&\frac{\eta}{2}\bm{M}^{-1/2}\bm{U}\bm{U}^{\bm{\top}}\bm{M}^{-1/2}\end{bmatrix}$

from which we obtain that the blocks of $\bm{Q}$ satisfy

	$\displaystyle\bm{Q}_{\delta}=\frac{\eta}{2}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\widehat{\bm{U}}\bm{\Lambda}_{n-1}^{-1}\widehat{\bm{U}}^{\bm{\top}}\bm{M}^{-1/2}\widetilde{\bm{C}},$		(34)
	$\displaystyle\bm{Q}_{\omega}=\frac{\eta}{2}\bm{M}^{-1/2}\bm{U}\bm{U}^{\bm{\top}}\bm{M}^{-1/2},$		(35)
	$\displaystyle\bm{Q}_{\delta\omega}=\bm{0}.$

Since the off-diagonal block matrix $\bm{Q}_{\delta\omega}$ of the variance matrix $\bm{Q}$ of the output (14) is a zero matrix and the stochastic process (13) is Gaussian, the invariant probability distribution of the frequency and the phase angle differences are independent.

(ii) Given the fact that $\bm{U}$ is an orthogonal matrix and $\bm{M}$ is a diagonal matrix, (35) is simply rewritten into (33). $\square$

In the proof, the fact is applied that two Gaussian distributed random variables are independent if and only if their covariance equals to zero. Due to the independence between the frequencies and the phase angle differences, Theorem 4.3(i) indicates that the fluctuations of the frequencies and those of the phase-angle differences have no stochastic relation with each other.

Formula (33) is verified in an example presented in Section 7. This formula shows the dependence of the variances of the frequencies at the nodes on the system parameters. First, the variance matrix $\bm{Q}_{\omega}$ is a diagonal matrix with $\bm{M}=\text{diag}(m_{i})\in\mathbb{R}^{n\times n}$ ; thus, the frequencies in different nodes are independent in the invariant probability distribution. Second, the variance of the frequency at each node increases linearly with the disturbance-damping ratio and is inversely proportional to the inertia of the synchronous machine at this node. This shows the importance of the inertia and the damping coefficient in suppressing the frequency deviation in the power network. However, increasing the inertia at a node suppresses the fluctuations of the frequency only at this node, without any effect on the other nodes. The vulnerable nodes are the ones with small inertia. Those nodes will then have large variances and these variances are not influenced by the disturbances at the other nodes. Finally, the parameters, the power generation and the loads, which determine the synchronous state $(\bm{\delta}^{*},\bm{0})$ and play roles in determining the value $l_{c_{ij}}$ as shown in (10), the line capacity and the network topology are all absent from the formula. It is surprising that these parameters have no impact on the variances of the frequencies. This might be due to the assumption of uniform disturbance-damping ratio in Assumption 4.1. Whether this occurs in the systems without the assumption still needs further study.

The trace of $\bm{Q}_{\omega}$ is derived directly from (33) as presented in the following corollary, which is actually the $\mathcal{H}_{2}$ norm of a linear input-output system (Poolla et al., 2017).

Corollary 4.4.

From Theorem 4.3 one obtains that

\displaystyle\text{trace}(\bm{Q}_{\omega})=\frac{\eta}{2}\sum_{i=1}^{n}\frac{1}{m_{i}}.

In order to reveal the dependence of the variances of the phase angle differences on the system parameters, we further deduce the formula of $\bm{Q}_{\delta}$ based on Lemma 4.2. Before presenting this explicit formula, we first introduce a lemma for the properties of the matrix $\widehat{\bm{C}}=\bm{R}^{1/2}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}$ .

Lemma 4.5.

Consider the symmetric matrix $\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}$ with $\widehat{\bm{C}}=\bm{R}^{1/2}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}$ . There exists an orthogonal matrix $\bm{W}\in\mathbb{R}^{m\times m}$ for $m\geq n$ such that

\displaystyle\bm{W}^{\bm{\top}}\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}\bm{W}=\bm{\Lambda}_{m},~~\bm{\Lambda}_{m}=\begin{bmatrix}\bm{\Lambda}_{n-1}&\bm{0}\\ \bm{0}&\bm{0}\end{bmatrix}\in\mathbb{R}^{m\times m}

(36)

with $\bm{\Lambda}_{n-1}$ defined in Lemma 4.2. If we denote $\bm{W}=[\bm{w}_{1},\bm{w}_{2},\cdots,\bm{w}_{n-1},\bm{w}_{n},\cdots,\bm{w}_{m}]$ , then the vector $\bm{w}_{i}$ is the orthonormal eigenvector of $\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}$ corresponding to the nonzero eigenvalue $\lambda_{i+1}$ for $i=1,\cdots,n-1$ and $\bm{w}_{i}$ for $i=n,\cdots,m$ are the orthonormal eigenvectors corresponding to the zero eigenvalue. For the case with $m=n-1$ , all the eigenvalues of $\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}$ are non-zero, and

\displaystyle\bm{\Lambda}_{m}=\bm{\Lambda}_{n-1},~~~~\bm{W}=[\bm{w}_{1},\bm{w}_{2},\cdots,\bm{w}_{n-1}].

Proof: For a connected graph, we have $\text{rank}(\widetilde{\bm{C}})=n-1$ , which leads to $\text{rank}(\widehat{\bm{C}})=n-1$ . Since the kernel of $\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}\bm{X}=\bm{0}$ and $\widehat{\bm{C}}^{\bm{\top}}\bm{X}=0$ are identical, $\text{rank}(\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}})=n-1$ . Based on Theorem A.2, we only need to prove that the non-zero diagonal elements of $\bm{\Lambda}_{m}$ are the non-zero eigenvalues of $\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}$ . We obtain from (15,16) that

	$\displaystyle\bm{U}^{\bm{\top}}\widehat{\bm{C}}^{\bm{\top}}\widehat{\bm{C}}\bm{U}$	$\displaystyle=\bm{U}^{\bm{\top}}\bm{M}^{-1/2}\widetilde{\bm{C}}\bm{R}\widetilde{\bm{C}}^{\bm{\top}}M^{-1/2}\bm{U}$
		$\displaystyle=\bm{U}^{\bm{\top}}\bm{M}^{-1/2}\bm{L}_{c}\bm{M}^{-1/2}\bm{U}=\bm{\Lambda}_{n}.$

With the left multiplication of $\widehat{\bm{C}}\bm{U}$ to the above equation, we obtain

\displaystyle\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}\widehat{\bm{C}}\bm{U}=\widehat{\bm{C}}\bm{U}\bm{\Lambda}_{n}.

(37)

We write $\bm{U}$ into the form $\begin{bmatrix}\bm{u}_{1}&\bm{u}_{2}&\bm{u}_{3}&\cdots&\bm{u}_{n}\end{bmatrix}$ . From Lemma 3.1 and $\widetilde{\bm{C}}^{\bm{\top}}\bm{1}_{n}=\bm{0}$ , we obtain $\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\bm{u}_{1}=\bm{0}$ , which leads to $\widehat{\bm{C}}\bm{u}_{1}=\bm{0}$ . Hence, we derive from (37) that

\displaystyle\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}\begin{bmatrix}\widehat{\bm{C}}\bm{u}_{2},\widehat{\bm{C}}\bm{u}_{3},\cdots,\widehat{\bm{C}}\bm{u}_{n}\end{bmatrix}=\begin{bmatrix}\lambda_{2}\widehat{\bm{C}}\bm{u}_{2}&\lambda_{3}\widehat{\bm{C}}\bm{u}_{2}&\cdots&\lambda_{n}\widehat{\bm{C}}\bm{u}_{n}\end{bmatrix},

which indicates that $\lambda_{i}$ and $\widehat{\bm{C}}\bm{u}_{i}$ for $i=2,\cdots,n$ are the eigenvalues and the corresponding eigenvectors of the matrix $\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}$ . $\square$

Based on Lemma 4.5, we present the explicit formula for the variance matrix $\bm{Q}_{\delta}$ in the following theorem.

Theorem 4.6.

Consider the system (12) with Assumption 4.1. The variance matrix of the phase angle differences in the invariant probability distribution satisfies

\displaystyle\bm{Q}_{\delta}=\frac{1}{2}\eta\bm{R}^{-1/2}\big{(}\bm{I}_{m}-\sum_{i=1}^{m-n+1}\bm{X}_{i}\bm{X}_{i}^{\bm{\top}}\big{)}\bm{R}^{-1/2}

(38)

where $\{\bm{X}_{i}\in\mathbb{R}^{m},i=1,2,\cdots,m-n+1\}$ is an orthonormal basis vector of the kernel of the matrix $\widetilde{\bm{C}}\bm{R}^{1/2}$ such that $\widetilde{\bm{C}}\bm{R}^{1/2}\bm{X}_{i}=\bm{0}$ . Clearly, because the inertia values are absent from the formula, they have no impact on the variance of the phase angle difference in each line.

Proof: From (34), we obtain

$\displaystyle\bm{Q}_{\delta}$	$\displaystyle=\frac{\eta}{2}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\widehat{\bm{U}}\bm{\Lambda}_{n-1}^{-1}\widehat{\bm{U}}^{\bm{\top}}\bm{M}^{-1/2}\widetilde{\bm{C}}$
	by Theorem A.2
	$\displaystyle=\frac{\eta}{2}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\big{(}\bm{M}^{-1/2}\bm{L}_{c}\bm{M}^{-1/2}\big{)}^{{\dagger}}\widehat{\bm{U}}\widehat{\bm{U}}^{\bm{\top}}\bm{M}^{-1/2}\widetilde{\bm{C}}$
	by $[\bm{u}_{1},\widehat{\bm{U}}][\bm{u}_{1},\widehat{\bm{U}}]^{\top}=\bm{u}_{1}\bm{u}_{1}^{\top}+\widehat{\bm{U}}\widehat{\bm{U}}^{\top}=\bm{I}_{n}$
	$\displaystyle=\frac{\eta}{2}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\big{(}\bm{M}^{-1/2}\bm{L}_{c}\bm{M}^{-1/2}\big{)}^{{\dagger}}(\bm{I}_{n}-\bm{u}_{1}\bm{u}_{1}^{\bm{\top}})\bm{M}^{-1/2}\widetilde{\bm{C}}$
	by $\widetilde{\bm{C}}^{\top}\bm{M}^{-1/2}\bm{u}_{1}=\bm{0}$ obtained from Lemma 3.1
	$\displaystyle=\frac{\eta}{2}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\big{(}\bm{M}^{-1/2}\bm{L}_{c}\bm{M}^{-1/2}\big{)}^{{\dagger}}\bm{M}^{-1/2}\widetilde{\bm{C}}$	(39)
	by (15)
	$\displaystyle=\frac{\eta}{2}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\big{(}\bm{M}^{-1/2}\widetilde{\bm{C}}\bm{R}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}\big{)}^{{\dagger}}\bm{M}^{-1/2}\widetilde{\bm{C}}$

where $(\cdot)^{{\dagger}}$ denotes the Moore-Penrose pseudo inverse of a matrix. With $\widehat{\bm{C}}=\bm{R}^{1/2}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}$ as in Lemma 4.5, we further obtain

\displaystyle\bm{Q}_{\delta}=\frac{\eta}{2}\bm{R}^{-1/2}\widehat{\bm{C}}\Big{(}\widehat{\bm{C}}^{\bm{\top}}\widehat{\bm{C}}\Big{)}^{{\dagger}}\widehat{\bm{C}}^{\bm{\top}}\bm{R}^{-1/2}.

(40)

By Lemma 4.5 and left multiplying (36) by $\widehat{\bm{C}}^{{\bm{\top}}}\bm{W}$ , we get

\displaystyle\widehat{\bm{C}}^{\bm{\top}}\widehat{\bm{C}}\widehat{\bm{C}}^{{\bm{\top}}}\bm{W}=\widehat{\bm{C}}^{{\bm{\top}}}\bm{W}\bm{\Lambda}_{m},

which indicates that the column vectors of $\widehat{\bm{C}}^{{\bm{\top}}}\bm{W}$ are the eigenvectors of $\widehat{\bm{C}}^{\bm{\top}}\widehat{\bm{C}}$ . We focus on the first $n-1$ eigenvectors $\widehat{\bm{C}}^{\bm{\top}}\bm{w}_{1},\cdots,\widehat{\bm{C}}^{\bm{\top}}\bm{w}_{n-1}$ in matrix $\widehat{\bm{C}}^{{\bm{\top}}}\bm{W}$ , which are orthogonal. The normalization of $\widehat{\bm{C}}^{\bm{\top}}\bm{w}_{i}$ for $i=1,\cdots,n-1$ yields

\displaystyle\lambda_{2}^{-1/2}\widehat{\bm{C}}^{\bm{\top}}\bm{w}_{1},\lambda_{3}^{-1/2}\widehat{\bm{C}}^{\bm{\top}}\bm{w}_{2},\cdots,\lambda_{n}^{-1/2}\widehat{\bm{C}}^{\bm{\top}}\bm{w}_{n-1}.

With these unit vectors, we obtain from Theorem A.2 that the Moore-Penrose pseudo inverse of $\widehat{\bm{C}}^{{\bm{\top}}}\widehat{\bm{C}}$ satisfies

\displaystyle\Big{(}\widehat{\bm{C}}^{{\bm{\top}}}\widehat{\bm{C}}\Big{)}^{{\dagger}}=\sum_{i=2}^{n}{\frac{1}{\lambda_{i}^{2}}(\widehat{\bm{C}}^{\bm{\top}}\bm{w}_{i-1})(\widehat{\bm{C}}^{\bm{\top}}\bm{w}_{i-1})^{{\bm{\top}}}}.

With (36), we further obtain

	$\displaystyle\widehat{\bm{C}}\Big{(}\widehat{\bm{C}}^{{\bm{\top}}}\widehat{\bm{C}}\Big{)}^{{\dagger}}\widehat{\bm{C}}^{{\bm{\top}}}$	$\displaystyle=\sum_{i=2}^{n}{\frac{1}{\lambda_{i}^{2}}\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}\bm{w}_{i-1}\bm{w}_{i-1}^{\bm{\top}}\widehat{\bm{C}}\widehat{\bm{C}}^{\bm{\top}}}=\sum_{i=2}^{n}{\bm{w}_{i-1}\bm{w}_{i-1}^{\bm{\top}}}$
		$\displaystyle=\bm{I}_{m}-\sum_{i=n}^{m}{\bm{w}_{i}\bm{w}_{i}^{\bm{\top}}}.$

By Lemma 4.5, $\bm{w}_{i}$ for $i=n,\cdots,m$ are the orthonormal eigenvectors corresponding to the zero eigenvalue such that $\bm{w}_{i}^{{\bm{\top}}}\widehat{\bm{C}}\widehat{\bm{C}}^{{\bm{\top}}}\bm{w}_{i}=0$ from which we obtain $\widehat{\bm{C}}^{{\bm{\top}}}\bm{w}_{i}=0$ . Since $\widehat{\bm{C}}^{{\bm{\top}}}=\bm{M}^{-1/2}\widetilde{\bm{C}}\bm{R}^{1/2}$ , $\widetilde{\bm{C}}\bm{R}^{1/2}\bm{w}_{i}=0$ which indicates that the vectors $\bm{w}_{i}$ for $i=n,\cdots,m$ form an orthonormal basis of the kernel of $\widetilde{\bm{C}}\bm{R}^{1/2}$ . Define $\bm{X}_{i}=\bm{w}_{i+n-1}$ for $i=1,\cdots,m-n+1$ to complete the proof. $\square$

Corollary 4.7.

If Assumption 4.1 holds and $l_{c_{ij}}=\gamma$ for all the lines, the variance matrix $\bm{Q}_{\delta}$ becomes

\displaystyle\bm{Q}_{\delta}=\frac{\eta}{2\gamma}\big{(}\bm{I}_{m}-\sum_{i=1}^{m-n+1}\bm{X}_{i}\bm{X}_{i}^{\bm{\top}}\big{)}

(41)

where $\bm{X}_{i}$ becomes the orthonormal basis of the kernel of the incidence matrix $\widetilde{\bm{C}}$ . Furthermore, the trace of $\bm{Q}_{\delta}$ is $\displaystyle\frac{\eta}{2\gamma}(n-1)$ .

The proof follows directly from Theorem 4.6 with $\bm{R}=\gamma\bm{I}_{m}$ and

\displaystyle\text{trace}\big{(}\displaystyle\sum_{i=1}^{m-n+1}\bm{X}_{i}\bm{X}_{i}^{{\bm{\top}}}\big{)}=\sum_{i=1}^{m-n+1}\bm{X}_{i}^{\bm{\top}}\bm{X}_{i}=m-n+1.

The trace of $\bm{Q}_{\delta}$ has been obtained by the $\mathcal{H}^{2}$ norm of input-output linear systems as in Poolla et al. (2017); Tegling et al. (2015), which is consistent with the result in the above corollary.

Following the procedure described in Appendix A.3, the vector $\bm{X}_{i}$ can be calculated from the basis vectors of the kernel of $\widetilde{\bm{C}}\bm{R}^{1/2}$ . Due to the non-uniqueness of the basis vectors $\bm{\xi}_{c}$ for $c=1,\cdots,m-n+1$ of the kernel of $\bm{C}$ , the set of the orthonormal basis vectors of the kernel of $\bm{C}\bm{R}^{1/2}$ is also non-unique. However, for the kernel, a set of orthonormal basis vectors can be obtained from any set of basis vectors by a linear transformation consisting of an orthogonal matrix. Such a transformation does not influence the calculation of the multiplication $\bm{X}_{i}\bm{X}_{i}^{\bm{\top}}$ . The explicit formula (38) of $\bm{Q}_{\delta}$ describes the dependence of the variances of the phase angle differences on the system parameters. It is shown that the variances of the phase angle differences increase linearly as the disturbance-damping ratio $\eta$ increases. Because the variance of the phase-angle differences does not depend on the inertia, the control objective of rotor angle stability hardly be improved by changing the virtual inertia. Here, the rotor angle stability is the ability of the phase angles to maintain their coherence.

In particular, formula (38) reveals the role of the network topology with weight $l_{c_{ij}}$ for line $e_{k}$ . In the complex network theory, the kernel of $\widetilde{\bm{C}}$ is the cycle space of the graph $\mathcal{G}$ . Hence, it follows from formula (38) that the stability of the power system is related to the cycle space of the graph. The way that changes in the topology of the power network affect the variances of the phase angle differences and hence stability can be investigated by a study of the cycle space of the graph. In Section 6, we make a further study on the impact of the network topology by studying the cycle space of graphs.

5 Bounds of the variance matrices for networks with non-uniform disturbance-damping ratio

In the previous sections, we discussed the roles of the parameters in systems with a uniform disturbance-damping ratio at the nodes. In this section, we present the findings for a system with non-uniform ratios. We define $\overline{\eta}=\max\{\eta_{i},i=1,\cdots,n\}$ and $\underline{\eta}=\min\{\eta_{i},i=1,\cdots,n\}$ with $\eta_{i}=b_{i}^{2}/d_{i}$ . For $\bm{A},\bm{B}\in\mathbb{R}^{n\times n}$ , we say that $\bm{A}\leq\bm{B}$ if the matrix $\bm{A}-\bm{B}$ is semi-negative-definite.

Lemma 5.1.

Define $\overline{\eta}=\max\{\eta_{i},i=1,\cdots,n\}$ and $\underline{\eta}=\min\{\eta_{i},i=1,\cdots,n\}$ with $\eta_{i}=b_{i}^{2}/d_{i}$ , and define $\overline{\bm{\beta}}=(\overline{\eta}\bm{D})^{1/2}$ and $\underline{\bm{\beta}}=(\underline{\eta}\bm{D})^{1/2}$ . The solution $\bm{Q}_{x}$ of the Lyapunov equation (26) satisfies the following inequalities where the various matrices are also defined

\displaystyle\bm{Q}_{\underline{\beta}}\leq\bm{Q}_{x}\leq\bm{Q}_{\overline{\beta}},

(42)

where

\displaystyle\bm{Q}_{\overline{\beta}}=\int_{0}^{\infty}e^{\bm{A}_{2}t}\overline{\bm{B}}_{2}\overline{\bm{B}}_{2}^{{\bm{\top}}}e^{\bm{A}_{2}^{{\bm{\top}}}t}\emph{d}t,~\bm{Q}_{\underline{\beta}}=\int_{0}^{\infty}e^{\bm{A}_{2}t}\underline{\bm{B}}_{2}\underline{\bm{B}}_{2}^{{\bm{\top}}}e^{\bm{A}_{2}^{{\bm{\top}}}t}\emph{d}t

with $\overline{\bm{B}}_{2},\underline{\bm{B}}_{2}\in\mathbb{R}^{(2n-1)\times n}$ such that

\displaystyle\overline{\bm{B}}_{2}=\begin{bmatrix}\bm{0}\\ \bm{U}^{{\bm{\top}}}\bm{M}^{-1/2}\overline{\bm{\beta}}\end{bmatrix},~~\underline{\bm{B}}_{2}=\begin{bmatrix}\bm{0}\\ \bm{U}^{{\bm{\top}}}\bm{M}^{-1/2}\underline{\bm{\beta}}\end{bmatrix}.

Proof: By the definition of $\overline{\bm{\beta}}$ and $\underline{\bm{\beta}}$ and $\underline{\eta}d_{i}\leq b_{i}^{2}=\eta_{i}d_{i}\leq\overline{\eta}d_{i})$ for all the nodes, we obtain

\displaystyle\underline{\eta}\text{diag}(d_{i})=\underline{\bm{\beta}}\underline{\bm{\beta}}^{\bm{\top}}\leq\widetilde{\bm{B}}\widetilde{\bm{B}}^{\bm{\top}}=\text{diag}(b_{i}^{2})\leq\overline{\bm{\beta}}\overline{\bm{\beta}}^{\bm{\top}}=\overline{\eta}\text{diag}(d_{i}).

Hence, with the definition of $\bm{B}_{2}$ in (23)

\displaystyle\underline{\bm{B}}_{2}\underline{\bm{B}}_{2}^{\bm{\top}}\leq\bm{B}_{2}\bm{B}_{2}^{{\bm{\top}}}\leq\overline{\bm{B}}_{2}\overline{\bm{B}}_{2}^{\bm{\top}}

which leads to (42). $\square$

Based on Lemma 5.1, we deduce bounds for $\bm{Q}_{\omega}$ and $\bm{Q}_{\delta}$ .

Theorem 5.2.

Consider the system (12). The variance matrix $\bm{Q}_{\omega}$ of the frequencies at the nodes satisfies

\displaystyle\frac{1}{2}\underline{\eta}\bm{M}^{-1}\leq\bm{Q}_{\omega}\leq\frac{1}{2}\overline{\eta}\bm{M}^{-1},

(43)

the variance matrix $\bm{Q}_{\delta}$ of the phase angle differences in the lines satisfies

\displaystyle\frac{1}{2}\underline{\eta}\widehat{\bm{Q}}\leq\bm{Q}_{\delta}\leq\frac{1}{2}\overline{\eta}\widehat{\bm{Q}},~~\widehat{\bm{Q}}=\bm{R}^{-1/2}\big{(}\bm{I}_{m}-\sum_{i=1}^{m-n+1}\bm{X}_{i}\bm{X}_{i}^{{\bm{\top}}}\big{)}\bm{R}^{-1/2},

(44)

where $\underline{\eta}$ and $\overline{\eta}$ are defined in Lemma 5.1 and $\bm{X}_{i}$ is as defined in Theorem 4.6.

Proof: In Lemma 5.1, the matrices $\underline{\bm{B}}_{2}$ and $\overline{\bm{B}}_{2}$ are defined such that the disturbance-damping ratio $b_{i}^{2}/d_{i}=\underline{\eta}$ and $b_{i}^{2}/d_{i}=\overline{\eta}$ for all the nodes respectively. Hence, using Lemma 4.2, $\bm{Q}_{\underline{\beta}}$ and $\bm{Q}_{\overline{\beta}}$ are solved explicitly as

\bm{Q}_{\underline{\beta}}=\begin{bmatrix}\frac{1}{2}\underline{\eta}\bm{\Lambda}_{n-1}^{-1}&\bm{0}\\ \bm{0}&\frac{1}{2}\underline{\eta}\bm{I}\end{bmatrix},~~\bm{Q}_{\overline{\beta}}=\begin{bmatrix}\frac{1}{2}\overline{\eta}\bm{\Lambda}_{n-1}^{-1}&\bm{0}\\ \bm{0}&\frac{1}{2}\overline{\eta}\bm{I}\end{bmatrix}.

From (42), we obtain

\displaystyle\bm{C}_{2}\bm{Q}_{\underline{\beta}}\bm{C}_{2}^{\bm{\top}}\leq\bm{Q}=\bm{C}_{2}\bm{Q}_{x}\bm{C}_{2}^{\bm{\top}}\leq\bm{C}_{2}\bm{Q}_{\overline{\beta}}\bm{C}_{2}^{\bm{\top}},

(45)

where $\bm{C}_{2}$ is the one in (21) or in (29) or in (30).

To prove (43), we consider the output as the frequency and take $\bm{C}_{2}$ in (30). Following the procedure to calculate the variances of the frequencies in Theorem 4.3 with $b_{i}^{2}/d_{i}=\overline{\eta}$ and $b_{i}^{2}/d_{i}=\underline{\eta}$ for all the nodes, we get

\displaystyle\bm{C}_{2}\bm{Q}_{\underline{\beta}}\bm{C}_{2}^{\bm{\top}}=\frac{1}{2}\underline{\eta}\bm{M}^{-1}~~\text{and}~~\bm{C}_{2}\bm{Q}_{\overline{\beta}}\bm{C}_{2}^{\bm{\top}}=\frac{1}{2}\overline{\eta}\bm{M}^{-1},

which lead to (43) with (45).

To prove (44), we consider the output as the phase angle differences and insert $\bm{C}_{2}$ of (29) into (45), then obtain the upper bound of $\bm{Q}_{\delta}$ from (39) such that

\displaystyle\hskip-15.0pt\bm{C}_{2}\bm{Q}_{\overline{\beta}}\bm{C}_{2}^{\bm{\top}}=\frac{1}{2}\overline{\eta}\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}(\bm{M}^{-1/2}\bm{L}_{c}\bm{M}^{-1/2})^{{\dagger}}\bm{M}^{-1/2}\widetilde{\bm{C}}.

Following the procedure to deduce the explicit formula in (40), we obtain

\displaystyle\widetilde{\bm{C}}^{\bm{\top}}\bm{M}^{-1/2}(\bm{M}^{-1/2}\bm{L}_{c}\bm{M}^{-1/2})^{{\dagger}}\bm{M}^{-1/2}\widetilde{\bm{C}}=\widehat{\bm{Q}}.

Hence, the upper bound of $\bm{Q}_{\delta}$ satisfies $\bm{C}_{2}\bm{Q}_{\overline{\beta}}\bm{C}_{2}^{\bm{\top}}=\frac{1}{2}\overline{\eta}\widehat{\bm{Q}}$ . Similarly, the lower bound satisfies $\bm{C}_{2}\bm{Q}_{\underline{\beta}}\bm{C}_{2}^{\bm{\top}}=\frac{1}{2}\underline{\eta}\widehat{\bm{Q}}$ . With these two bounds and (45), we obtain (44). $\square$

It is well known that the diagonal elements of a semi-positive definite symmetric matrix are all non-negative. Hence, the bounds of the variances of the frequencies at the nodes and the phase angle differences in the lines are derived directly from (43) and (44).

Formula (43) reveals the factors that impact the variances of the frequencies at nodes in networks with a non-uniform disturbance-damping ratio. First, as in networks with a uniform disturbance-damping ratio, the inertias of the synchronous machines locally impact the variances of the frequencies at the nodes, and the network topology and the parameter $l_{c_{ij}}$ have little impact because they are absent in the formula. Second, in networks with a non-uniform disturbance-damping ratio, the variances of the frequencies will increase as the minimum value $\underline{\eta}$ increases and decrease as the maximum value $\overline{\eta}$ decreases. Hence, by decreasing all the disturbance-damping ratios, the variances of the frequencies will be decreased, which is consistent with the findings in networks with a uniform disturbance-damping ratio. In addition, by decreasing the maximum value $\overline{\eta}$ , there are nodes at which the variances of the frequencies will be decreased.

Formula (44) illustrates the roles played by the system parameters in determining the variances of the phase angle differences in networks with a non-uniform disturbance-damping ratio. First, the roles of the values $\underline{\eta}$ and $\overline{\eta}$ in determining the variances of the phase angle differences is the same as that in determining the variances of the frequencies. Decreasing the largest disturbance-damping ratio can decrease the variances of the phase angle differences at some lines. For example, energy storage in combination with droop control, which affects the parameter $d_{i}$ at the relevant nodes, will directly decrease the disturbance-damping ratios. Second, as in a network with a uniform disturbance-damping ratio, the inertia is absent from the formula, and the role of the network topology is also reflected by the basis of the cycle space. Hence, the inertia has little impact on the variances of the phase angle differences, and by forming small cycles, the variances of the phase angle differences can also be effectively decreased in the network. Third, the impact of constructing new lines to form cycles and increasing the capacities of the lines on the upper and lower bounds are the same as in the networks with a uniform disturbance-damping ratio.

In regard to the impact of the scales of the power systems on the stability, we have the following conclusion. From formulas (33,38, 43,44), we see that, if the scale of the network is increased by constructing nodes that have small effects on the power flows and possess disturbance-damping ratios close to $\eta$ , the fluctuations in the frequency or in the phase angle differences in the network will not be dramatically increased or decreased. Hence, the stability will be changed little by increasing the scale of the network. This follows formula (33) for networks with a uniform disturbance-damping ratio, which states that the newly connected nodes with disturbance-damping ratios equal to $\eta$ will not bring fluctuations to the frequency at the other nodes. Since $\delta_{i}^{*}\approx\delta_{j}^{*}$ for all the nodes, the newly connected nodes have little influence on the phase angle difference in the synchronous state, and it is indicated by formula (38) that the fluctuation of the phase angle difference will not change greatly. Similarly, for networks with a non-uniform disturbance-damping ratio, the newly connected nodes with disturbance-damping ratios in the set $[\underline{\eta},\overline{\eta}]$ will not change the bounds of the variance, as follows from the formulas (43) and (44). This conclusion is different from that obtained by a study of linear stability (Xi et al., 2017), where the linear stability decreases if the scale of the network increases. However, if nodes that consume a large amount of power and have large disturbance-damping ratios $b_{i}^{2}/d_{i}$ are added to the network, the variance of the frequency and the phase angle difference may increase because the weights of lines may decrease and the disturbances may propagate from these nodes to the other nodes in the network.

6 The role of the network topology

To fully explore the role of the network topology from the formula (38), we introduce three concepts for graphs,

Definition 6.1.

Consider a connected and undirected graph $\mathcal{G}$ . (i) A single line is defined as a line that does not belong to any cycle; (ii) Line $e_{1}$ is called a cycle-shared line of line $e_{2}$ if there exists at least one cycle containing both $e_{1}$ and $e_{2}$ ; (iii) A cycle-cluster is a subgraph of $\mathcal{G}$ obtained in the following way. One starts from a subgraph of one cycle and extends it by adding the lines in all the cycles with which the subgraph has at least one line in common, then one obtains a cycle-cluster.

It is deduced that a graph is composed of cycle-clusters and single lines, a line either belongs to a cycle-cluster or is a single line and in a cycle-cluster each pair of lines are cycle-shared lines. In the following example, we explain the definitions and the formulation of the basis vectors of the cycle space.

Refer to caption — Figure 1: A network with two cycle-clusters and a single line.

Example 6.2.

Consider the network show in Fig. 1. There are two cycle-clusters, i.e., $\{e_{1},e_{2},e_{7}\}$ and $\{e_{4},e_{5},e_{6},e_{8},e_{9}\}$ , and a single line $e_{3}$ that does not belong to any cycle. Each pairs of lines in the cycle-cluster $\{e_{1},e_{2},e_{7}\}$ are cycle-shared respectively, similarly for the lines in the cycle-cluster $\{e_{4},e_{5},e_{6},e_{8},e_{9}\}$ . However, two lines belonging to two different cycle-clusters are not cycle-shared, because a cycle containing both of these two lines cannot be found, for example $e_{1}$ and $e_{4}$ . The directions of lines are specified for the formulation of the incidence matrix $\widetilde{\bm{C}}$ and the calculation of the basis vectors of the cycle space. The directions of all the cycles are clock-wise. Following the procedure to calculate the basis vectors of the cycle space in Appendix A.3, we get the basis vectors of the cycle space of this network, $\bm{\xi}_{1}=\begin{bmatrix}-1&1&0&0&0&0&-1&0&0\end{bmatrix}^{\bm{\top}}$ , $\bm{\xi}_{2}=\begin{bmatrix}0&0&0&-1&1&0&0&0&1\end{bmatrix}^{\bm{\top}}$ , and $\bm{\xi}_{3}=\begin{bmatrix}0&0&0&0&0&1&0&-1&-1\end{bmatrix}^{\bm{\top}}$ which are corresponding to the fundamental cycles $\{e_{1},e_{2},e_{7}\}$ , $\{e_{4},e_{5},e_{9}\}$ and $\{e_{6},e_{8},e_{9}\}$ respectively. Obviously, $\bm{\xi}_{1}$ is orthogonal to $\bm{\xi}_{2}$ and $\bm{\xi}_{3}$ . This indicates that the basis vectors corresponding the cycles in different cycle-clusters are orthogonal. Due to the non-uniqueness of the spanning tree selected to form the fundamental cycles, the basis vectors are also non-unique. Thus the basis vectors of the cycle space of the network in Fig. 1 can also be $\bm{\xi}_{1}=\begin{bmatrix}-1&1&0&0&0&0&-1&0&0\end{bmatrix}^{\bm{\top}}$ , $\bm{\xi}_{2}=\begin{bmatrix}0&0&0&-1&1&0&0&0&1\end{bmatrix}^{\bm{\top}}$ , and $\bm{\xi}_{3}=\begin{bmatrix}0&0&0&-1&1&1&0&-1&0\end{bmatrix}^{\bm{\top}}$ which are corresponding to the fundamental cycles $\{e_{1},e_{2},e_{7}\}$ , $\{e_{4},e_{5},e_{9}\}$ and $\{e_{4},e_{6},e_{8},e_{5}\}$ respectively.

The network topology has two effects on the stability of the power system: the power flows at the synchronous state $(\bm{\delta}^{*},\bm{0})$ and the variance of the phase angle differences. Formula (38) indicates that the variance also depends on the power flows because $R_{k}=l_{c_{ij}}$ and $l_{c_{ij}}=l_{ij}\cos{\delta_{ij}^{*}}$ . This demonstrates the nonlinear character of the impacts of the network topology on stability. A network can be constructed mathematically in two steps, i.e., first connecting all the nodes to form a tree network and then constructing new lines or replacing the existing lines by ones with larger capacities. By following these steps, in addition to investigating the tree network, we reveal the role of the network topology by studying the impact of constructing new lines and increasing the capacity of the lines.

For the power flows, we have the following proposition.

Proposition 6.3.

Consider the power system (1) with a synchronous state that satisfies the security condition (11). (i) If the capacity of a single line is increased, then the power flows in all the other lines remain unchanged. (ii) If in a cycle-cluster a new line is constructed or the capacity of a line is increased, the power flows in the lines that are not in this cycle-cluster, remain unchanged.

Proof: Without loss of generality, we assume there three sub-graphs in graph $\mathcal{G}$ , i.e., $\mathcal{G}_{1}(\mathcal{V}_{1},\mathcal{E}_{1})$ , $\mathcal{G}_{2}(\mathcal{V}_{2},\mathcal{E}_{2})$ and $\mathcal{G}_{3}(\mathcal{V}_{3},\mathcal{E}_{3})$ where $\mathcal{G}_{1}(\mathcal{V}_{1},\mathcal{E}_{2})$ is either a cycle-cluster or single-line, $\mathcal{V}_{1}\cup\mathcal{V}_{2}\cup\mathcal{V}_{3}=\mathcal{V}$ , $\mathcal{E}_{1}\cup\mathcal{E}_{2}\cup\mathcal{E}_{3}=\mathcal{E}$ , $\mathcal{E}_{i}\cap\mathcal{E}_{j}=\emptyset$ for $i\neq j$ , $\mathcal{V}_{1}\cap\mathcal{V}_{2}=\{k\}$ , $\mathcal{V}_{1}\cap\mathcal{V}_{3}=\{q\}$ and $\mathcal{V}_{2}\cap\mathcal{V}_{3}=\emptyset$ . We prove that the power flows in the lines in $\mathcal{G}_{2}$ remain unchanged when the capacity of a line is increased or a new line is constructed in $\mathcal{G}_{1}$ . In the power flow calculation, we choose node $k$ as the reference node with $\delta_{k}=0$ . Thus, the power flow in the $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ are decoupled, where the power flows in cycle-cluster $\mathcal{G}_{2}$ satisfy

	$\displaystyle P_{i}-\sum_{j\in\mathcal{V}_{2}}l_{ij}\sin{(\delta_{i}-\delta_{j})}$	$\displaystyle=0,~i\in\mathcal{V}_{2}~~\text{and}~~i\neq k,$
	$\displaystyle\delta_{k}$	$\displaystyle=0.$

that are not changed by adding new lines or increasing the capacities of lines in cycle-cluster $\mathcal{G}_{1}$ . Similarly, it is proven that the power flows in $\mathcal{G}_{3}$ remain unchanged by constructing new lines or increasing the capacity of the lines in $\mathcal{G}_{1}$ . $\square$

From Proposition 6.3, we obtain that the phase angle differences $\delta_{ij}^{*}$ in lines at the synchronous state in a cycle-cluster are independent of the power flows in the other cycle-clusters. Hence, the weights $l_{c_{ij}}=l_{ij}\cos{\delta_{ij}^{*}}$ of lines in the cycle-cluster will not be changed by adding lines or increasing line capacity in the other cycle-clusters.

Based on the theory of the cycle space, we obtain the following Corollary of Theorem 4.6.

Corollary 6.4.

Consider the system (12) with Assumption 4.1.

(i)

The invariant probability distribution of the phase angle difference in a single line connecting nodes $i$ and $j$ is independent of those of the phase angle differences in all the other lines in the network, and the variance of the phase angle differences in this line is $\displaystyle\frac{1}{2}\eta l_{c_{ij}}^{-1}$ .
(ii)

According to the invariant probability distribution, the phase-angle differences of all lines in a particular cycle cluster are independent of the phase-angle differences of all lines which are not in this cycle cluster.
(iii)

Increasing the weight of a line or constructing new lines in a cycle-cluster without changing the weights of all the other lines decreases the variances of the phase angle differences in the lines of this cycle-cluster.

(iv)

For a cycle-cluster with only one cycle with lines in set $\mathcal{E}_{c}$ in the graph, the variance of the phase angle differences in the line connecting nodes $i$ and $j$ in this cycle-cluster is

\displaystyle\frac{\eta}{2}\Big{(}l_{c_{ij}}^{-1}-l_{c_{ij}}^{-2}\Big{(}\sum_{(r,q)\in\mathcal{E}_{c}}l_{c_{rq}}^{-1}\Big{)}^{-1}\Big{)}.

(46)

If $l_{c_{ij}}=\gamma$ for all the lines in this cycle, the variances of the phase angle differences in these lines are $\displaystyle\frac{\eta}{2\gamma}(1-\frac{1}{N})$ where $N$ is the length of the cycle.

Proof: (i) For an acyclic network, it follows from (38) that the variance matrix of the phase angle difference is $\displaystyle\frac{\eta}{2}\bm{R}^{-1}$ because the cycle space of the acyclic network is empty. Thus, the variance in line $e_{k}$ is $\displaystyle\frac{1}{2}\eta l_{c_{ij}}^{-1}$ . For a network with cycles and single lines, without loss of generality, assume line $e_{1}$ is a single-line. Following the method to formulate the basis of the cycle space in Appendix A.3, the base vector has the form $\bm{\xi}_{i}=\begin{bmatrix}0&\xi_{i,2}&\xi_{i,3}&\cdots&\xi_{i,m}\end{bmatrix}^{\bm{\top}}$ where $\xi_{i,j}$ is either $-1$ , $1$ or $0$ , and $\bm{X}_{i}$ has the form $\bm{X}_{i}=\begin{bmatrix}0&x_{i,2}&x_{i,3}&\cdots&x_{i,m}\end{bmatrix}^{\bm{\top}}$ obtained by Gram-Schmidt orthogonalization of $\bm{R}^{-1/2}\bm{\xi}_{i}$ . Because the elements in the first column and the first row of $\bm{X}_{i}\bm{X}_{i}^{\bm{\top}}$ are all zero, we derive the independence of the invariant probability distribution of the phase angle difference in this line to those of the phase angle difference in all the other lines. By (38), we obtain that the variance in this line is $\displaystyle\frac{\eta}{2l_{c_{ij}}}$ .

(ii) We partition the graph $\mathcal{G}$ into two sub-graphs, $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ , where $\mathcal{G}_{1}$ is either a cycle-cluster or a single line. If $\mathcal{G}_{1}$ is a single line, we obtain this conclusion directly from Corollary 6.4(i) directly. We now consider the case where $\mathcal{G}_{1}$ is a cycle-cluster. Denote the number of lines in these two sub-graphs by $N$ and $m-N$ , the number of fundamental cycles by $m_{1}$ and $m_{2}$ , the lines in $\mathcal{G}_{1}$ by $e_{1},\cdots,e_{N}$ and those in $\mathcal{G}_{2}$ by $e_{N+1},\cdots,e_{m}$ respectively. Here, $m_{1}+m_{2}=m-n+1$ . The basis vectors of the cycles in $\mathcal{G}_{1}$ have the form $\bm{\xi}_{i}=\begin{bmatrix}\xi_{i,1}&\xi_{i,2}&\cdots&\xi_{i,N}&0&\cdots&0\end{bmatrix}^{\bm{\top}}$ for $i=1,\cdots,m_{1}$ and those of the cycles in $\mathcal{G}_{2}$ have the form $\bm{\xi}_{i}=\begin{bmatrix}0&0&\cdots&0&\xi_{i,{N+1}}&\cdots&\xi_{i,m}\end{bmatrix}^{\bm{\top}}$ for $i=m_{1}+1,\cdots,m-n+1$ . In these vectors, $\xi_{i,j}$ are either $1$ , $-1$ , or $0$ . By Gram-Schmidt orthogonalization of $\bm{R}^{-1/2}\bm{\xi}_{i}$ , we get the orthonormal vectors $\bm{X}_{i}=\begin{bmatrix}x_{i,1}&\cdots&x_{i,N}&0&\cdots&0\end{bmatrix}^{\bm{\top}}$ for $i=1,\cdots,m_{1}$ and $\bm{X}_{i}=\begin{bmatrix}0&\cdots&0&x_{i,{N+1}}&\cdots&x_{i,m}\end{bmatrix}^{\bm{\top}}$ for $i=m_{1}+1,\cdots,m-n+1$ . It is obvious that the entries in the first $N$ columns and the first $N$ rows of the matrix $\displaystyle\sum_{i=m_{1}+1}^{m-n+1}\bm{X}_{i}\bm{X}_{i}^{\bm{\top}}$ are all $0$ . This indicates that the lines in $\mathcal{G}_{2}$ have no contributions to the first $N$ columns and the first $N$ rows of $\bm{Q}_{\delta}$ . Similarly, the lines in $\mathcal{G}_{1}$ have no contributions to the last $m-N$ columns and the last $m-N$ rows of $\bm{Q}_{\delta}$ . Hence, the invariant probability distribution of the phase angle differences in the lines of $\mathcal{G}_{1}$ are independent of those in the lines of $\mathcal{G}_{2}$ .

(iii) The case in which the weight of a line in a cycle-cluster increases is considered first. Assume the graph is a cycle-cluster, where the weight of line $e_{1}$ increases. Denote the dimension of the kernel of $\widetilde{\bm{C}}\bm{R}^{1/2}$ by $N$ , which equals to $m-n+1$ . Thus, there are $N$ fundamental cycles in the cycle-cluster. The basis vectors are chosen below. The basis vectors corresponding to the $N-1$ fundamental cycles which do not include line $e_{1}$ have the form $\bm{\xi}_{i}=\begin{bmatrix}0&\xi_{i,2}&\xi_{i,3}&\cdots&\xi_{i,m}\end{bmatrix}^{\bm{\top}}$ for $i=1,\cdots,N-1$ , where $\xi_{i,q}=1,-1~\text{or}~0$ for $q=2,\cdots,m$ and that corresponding to the fundamental cycle which includes line $e_{1}$ has the form $\bm{\xi}_{N}=\begin{bmatrix}\xi_{N,1}&\xi_{N,2}&\cdots&\xi_{N,m}\end{bmatrix}^{\bm{\top}}$ where $\xi_{N,1}=1~\text{or}~-1$ and $\xi_{N,q}=1,-1~\text{or}~0$ for $q=2,\cdots,m$ . This can be done by changing the basis vectors of the cycle space properly. By the Gram-Schmidt orthogonalization of $\bm{R}^{-1/2}\bm{\xi}_{i}$ , we obtain $\bm{X}_{i}=\begin{bmatrix}0&x_{i,2}&x_{i,3}&\cdots&x_{i,m}\end{bmatrix}^{\bm{\top}}$ for $i=1,\cdots,N-1$ which is independent of the weight $l_{1}$ of line $e_{1}$ . The last unit vector $\bm{X}_{N}$ can be obtained by the normalization of the vector $\bm{X}_{N}^{\prime}=\bm{R}^{-1/2}\bm{\xi}_{N}-\alpha_{1}\bm{X}_{1}-\cdots-\alpha_{N-1}\bm{X}_{N-1}$ where $\displaystyle\alpha_{i}=\frac{(\bm{R}^{-1/2}\bm{\xi}_{N})^{\bm{\top}}\bm{X}_{i}}{\bm{X}_{i}^{\bm{\top}}\bm{X}_{i}}$ . Because the first element of $\bm{X}_{i}$ is zero for $i=1,\cdots,N-1$ , $\alpha_{i}$ is independent of $l_{1}$ . Hence $\bm{X}_{N}^{\prime}$ has the form $\bm{X}_{N}^{\prime}=\begin{bmatrix}l_{1}^{-1/2}\xi_{N,1}&x^{\prime}_{N,2}&x^{\prime}_{N,3}&\cdots&x^{\prime}_{N,m}\end{bmatrix}^{\bm{\top}}$ where $x^{\prime}_{N,q}$ is independent of $l_{1}$ for $q=2,\cdots,m$ . By the normalization of $\bm{X}_{N}^{\prime}$ , we obtain $\bm{X}_{N}=a\bm{X}_{N}^{\prime}$ where $a=\big{(}l_{1}^{-1}+\displaystyle\sum_{i=2}^{m}{x^{\prime}_{N,i}}^{2}\big{)}^{-1/2}$ . Hence, the diagonal element of $\bm{X}_{N}\bm{X}_{N}^{\top}$ equals to $a^{2}l_{1}^{-1}$ for $i=1$ and equals to $a^{2}{x^{\prime}_{N,i}}^{2}$ for $i=2,\cdots,m$ . Inserting $\bm{X}_{N}\bm{X}_{N}^{\bm{\top}}$ into (41), we obtain the variance of the phase angle difference in line $e_{1}$ which equals to $\displaystyle\frac{1}{2}\eta(l_{1}^{-1}-a^{2}l_{1}^{-2})$ and that in line $e_{q}$ which equals to $\displaystyle\frac{1}{2}\eta l_{q}^{-1}(1-a^{2}{x^{\prime}_{N,q}}^{2}-\sum_{i=1}^{N-1}{x_{i,q}^{2}})$ for $q=2,\cdots,m$ . It is obvious that if $l_{1}$ increases, these variances decrease.

We next consider the case when a new line is constructed in a cycle-cluster without changing the weight of all the other lines. Assume line $e_{1}$ is the new line. Following the above calculation, we obtain that the variance in the line with weight $l_{q}$ equals to $\displaystyle\frac{1}{2}\eta l_{q}^{-1}(1-a^{2}{x^{\prime}_{N,q}}^{2}-\sum_{i=1}^{N-1}{x_{i,q}^{2}})$ . For the variances in lines before constructing line $e_{1}$ , by choosing the basis vector corresponding to the $N-1$ fundamental cycles which do not include line $e_{1}$ and the Gram-Schmidt orthogonalization of these vectors, we obtain the variance in line $e_{q}$ with weight $l_{q}$ is $\displaystyle\frac{1}{2}\eta l_{q}^{-1}(1-\sum_{i=1}^{N-1}{x_{i,q}^{2}})$ for $q=2,\cdots,m$ . Clearly, the variance decreases after adding line $e_{1}$ .

(iv) The lines in the cycle are denoted by $e_{1},e_{2},\cdots,e_{N}$ with weights $l_{1},l_{2},\cdots,l_{N}$ . Assume the direction of these lines are consistent with the direction of the cycle. The vectors corresponding to this cycle and the other cycles are denoted by $\bm{\xi}_{1}$ and $\bm{\xi}_{i}$ with $i=2,\cdots,m-n+1$ respectively. Following Appendix A.3, we obtain $\bm{\xi}_{1}=\begin{bmatrix}1&1&\cdots&1&0&\cdots&0\end{bmatrix}^{\bm{\top}}$ where the first $N$ elements equal to $1$ and the last $m-N$ elements equal to $0$ , and $\bm{\xi}_{i}=\begin{bmatrix}0&0&\cdots&0&\xi_{i,{N+1}}&\cdots&\xi_{i,m}\end{bmatrix}^{\bm{\top}}$ where the first $N$ elements are all $0$ and the last $m-N$ elements equal to either $0$ , $1$ or $-1$ . Obviously, the vector $\bm{R}^{-1/2}\bm{\xi}_{1}$ is orthogonal to the vector $\bm{R}^{-1/2}\bm{\xi}_{i}$ for $i=2,\cdots,m-n+1$ . By Gram-Schmidt orthogonalization, we derive $\bm{X}_{1}=\Big{(}\sum_{k=1}^{N}l_{k}^{-1}\Big{)}^{-1/2}\begin{bmatrix}l_{1}^{-1/2}&l_{2}^{-1/2}&\cdots&l_{N}^{-1/2}&0&\cdots&0\end{bmatrix}^{\bm{\top}}$ from $\bm{R}^{-1/2}\bm{\xi}_{1}$ and $\bm{X}_{i}=\begin{bmatrix}0&0&\cdots&0&x_{i,{N+1}}&\cdots&x_{i,m}\end{bmatrix}^{\bm{\top}}$ for the linear subspace composed of the vectors $\bm{R}^{-1/2}\bm{\xi}_{i}$ with $i=2,\cdots,m-n+1$ . Because the first $N$ elements of $\bm{X}_{i}$ for $i=2,\cdots,m-n+1$ are all $0$ , the matrix $\bm{X}_{i}\bm{X}_{i}^{\bm{\top}}$ has no contributions to the first $N$ columns and the first $N$ rows of $\bm{Q}_{\delta}$ . Hence, the invariant probability distribution of the phase angle differences in the lines of the cycle are independent from those in the other lines. Further more, by (38), we obtain that the $k$ th diagonal element of $\bm{Q}_{\delta}$ for $k=1,\cdots,N$ is

\frac{\eta}{2}\Big{(}l_{k}^{-1}-l_{k}^{-2}\big{(}\sum_{r=1}^{N}l_{r}^{-1}\big{)}^{-1}\Big{)}

from which we obtain (46) by replacing $l_{k}$ by $l_{c_{ij}}$ for line $e_{k}$ . If $l_{k}=\gamma$ for $k=1,\cdots,N$ , we further get the first $N$ diagonal elements of $\bm{Q}_{\delta}$ equal to $\displaystyle\frac{\eta}{2\gamma}(1-\frac{1}{N})$ . $\square$

Remark 6.5.

From Proposition 6.3 and Corollary 6.4, we get the following findings.
(i) The variance of the phase angle difference in a single line connecting nodes $i$ and $j$ is $\displaystyle\frac{\eta}{2}l_{c_{ij}}^{-1}$ , which is not influenced by either constructing a new line without forming a cycle-cluster that includes this line or increasing the capacities of the other lines. Thus, a single line is likely to be a vulnerable line. This is because neither the construction of new lines nor the increase in the capacity of the other lines changes the power flow $l_{ij}\sin{\delta_{ij}^{*}}$ in this line, which is stated in Proposition 6.3, and the invariant probability distribution of the phase angle difference in the single line is independent of those of the phase angle differences in all the other lines, which is obtained from Corollary 6.4-(i).
(ii) Constructing new lines and increasing the capacities of lines in a cycle-cluster have no impact on the variances of the phase angle differences in the lines that are not in this cycle-cluster. This is because constructing new lines or increasing the capacities of lines in a cycle-cluster has no influence on the power flows in other cycle-clusters and single lines, which is indicated by Propostion 6.3, and the invariant probability distribution of the phase angle differences in the lines of a cycle-cluster is independent of those in the lines that are not in this cycle-cluster, which is demonstrated by Corollary 6.4-(ii).
(iii) By either increasing the weights $l_{c_{ij}}$ of lines or constructing new lines without changing the weights of the other lines in a cycle-cluster, the variances of the phase angle differences in this cycle-cluster will decrease.
(iv) For a cycle-cluster with only one cycle with lines in set $\mathcal{E}_{c}$ in the graph, the variance of the phase angle difference in the line connecting nodes $i$ and $j$ can be calculated from (46). In addition, based on (iii), we obtain that formula (46) provides a conservative estimation of the variances in the lines in cycle-clusters, i.e., the variance in a line that is in multiple cycles can be approximated by formula (46) by taking the smallest cycle that includes this line.

These findings provide guidelines on how to reduce the negative effects of vulnerable lines and designing future power networks, which should have low variances in phase angle differences when subjected to stochastic disturbances from power sources and power loads. The term remedy will be used for the reduction of these negative effects. Changing a power network by adding lines to form small cycles or by increasing the capacity of particular lines will suppress the fluctuations in the phase differences in the lines of the corresponding cycle-cluster. The benefit of forming small cycles is that the fluctuations in the phase angle differences decrease by $O(1/N)$ , where $N$ denotes the length of the cycle. This is consistent with the findings obtained by studying the energy barrier of a nonlinear system with a cyclic network in Xi et al (Xi et al., 2017). The fluctuations in the phase angle differences can be decreased by replacing transmission lines with small line capacities by ones with large line capacities. This is the same rule as for the transient stability analysis of the Single Machine Infinite Bus (SMIB) model (Kundur, 1994) by the equal area criterion. Because the variances of the phase angle differences decrease linearly with the parameter $l_{c_{ij}}=l_{ij}\cos{\delta_{ij}^{*}}$ , the control of the power flows to increase the value $\cos{\delta_{ij}^{*}}$ can also decrease the fluctuations of the phase angle differences in the lines. These findings will be further explained in an example in Section 7.

7 Case study

In this section, we verify the formulas (33, 38) for the networks with uniform disturbance-damping ratio, the bounds (43, 44) for the variance matrices for the networks with non-uniform disturbance-damping ratio and the findings presented in Remark 6.5. We take the 500KV transmission network of Shandong Province of China (Ye et al., 2016) as an example.

Example 7.1.

Consider the 500 KV transmission network of Shandong Province as shown in Fig. 2. There are 5 nodes with generators and 18 nodes with loads only. The nodes of squares denote power generators and the nodes of cycles denotes power loads. Line $e_{4}$ does not exist in practice, which is constructed virtually in order to explain our findings. Before constructing line $e_{4}$ , all the red lines are single lines and all the black lines are in a cycle-cluster which has 6 fundamental cycles. After constructing line $e_{4}$ , there is one more cycle-cluster, which is composed of $(e_{1},e_{2},e_{3},e_{4})$ . We set $m_{i}=i$ for the generators and $m_{i}=1$ for all the loads. We study the 7 cases with different settings of $P_{i}$ , $l_{ij}$ and $b_{i}^{2}/d_{i}$ as shown in Table 1.

In cases 1-3, it holds that the phase angle difference $\delta_{ij}^{*}=0$ because of $P_{i}=0$ for all the nodes. Thus, when disturbances occur, the frequencies at the nodes and the power flows in the lines fluctuate around zero. The weights of the lines satisfy $l_{c_{ij}}=l_{ij}$ . The weights of the lines in Cases 4-6 are shown in Table 3, which are calculated by solving the power flow equations. The variances of the frequencies at the nodes and the phase angle differences in the lines are presented in Tables 2 and 4, respectively. The values in the tables are first calculated by formulas (33) and (38) and then verified using Matlab following the procedure in Theorem 3.2. In Cases 1-3, because $\delta_{ij}^{*}=0$ for all the lines in the networks, the power flows are independent of the network topology. In this case, the impact of the network topology alone on the variance of the phase angle difference can be observed. In Cases 4-6, because $P_{i}$ is nonzero, updating the network topology, such as constructing new lines and increasing the line capacities, may change the weight $l_{c_{ij}}$ or the cycle space. Hence, the overall impact of the network topology can be analysed. In case 3 and 6, the capacity of line $e_{23}$ is increased from $10$ to $20$ in order to observe the changes of the variances in the other lines. This line is selected because in Case 4 the variance in this line is the largest one in the cycle-cluster that includes this line, as shown in Table 4.

First, let us focus on the variances of the phase angle differences in the single lines. Lines $e_{1}-e_{8}$ in cases 1, 4, and $e_{5}-e_{8}$ in cases 2-3, 5-6 are single lines. It is verified in Table 4 that the variances of the phase angle differences in these lines equal $\displaystyle\frac{\eta}{2}l_{c_{ij}}^{-1}$ with the weights of the lines shown in Table 3. In particular, the variances in lines $e_{5}-e_{8}$ are affected neither by constructing $e_{4}$ in cases 2 and 5 nor by increasing the capacity of $e_{23}$ in cases 3 and 6. This verifies the finding in Remark 6.5-(i).

Second, by comparing the weights and the variances in the lines in case 4 with those in case 5, it may be noted in Table 3 and 4 that both the weights and the variances in $e_{5}-e_{29}$ are not changed when $e_{4}$ is constructed in case 5. This is because these lines are not in the cycle-cluster that includes $e_{4}$ . Similarly, by comparing the weights and the variances in the lines in case 5 with those in case 6, it is seen in Table 3 and 4 that both the weights and the variances in $e_{1}-e_{8}$ are not influenced by increasing the capacity of $e_{23}$ . This is due to the fact that these lines are not in the cycle-cluster that includes $e_{23}$ . Thus, the findings in Remark 6.5-(ii) is verified .

Third, we evaluate the findings in Remark 6.5-(iii). The effects of constructing $e_{4}$ have already been analyzed, where the variances in the lines in the cycle-cluster of $(e_{1},e_{2},e_{3},e_{4})$ all decrease while those in the other lines are not affected. When comparing the variances in the lines in case 2 with those in case 3 in Table 4, it is found that the variances in $e_{11},e_{14},e_{17}-e_{29}$ all decrease after increasing the capacity of $e_{23}$ from $10$ to $20$ . We remark that those in $e_{9},e_{10},e_{12},e_{13}$ also decrease, which are not explicitly shown in the table because of the limited precision. This indicates that the variances of the lines in a cycle-clusters all decrease if the capacity of a line in this cycle-cluster increases. However, in practice, constructing new lines or increasing the capacity of lines also changes the power flows, which further influence the weight $l_{c_{ij}}$ . For example, when comparing the weights in case 5 with those in case 6, it is shown in Table 3 that after increasing the capacity of $e_{23}$ in case 6, the weights of $e_{24},e_{25},e_{26}$ decrease from $9.7528,9.9266,9.9978$ to $9.6934,9.8933,9.9896$ respectively. We remark that similar as in case 3, the variances in $e_{9},e_{10},e_{12},e_{15},e_{16}$ also decrease, which are not explicitly shown due to the limitation of the precision. Although only some of the weights decrease, as shown in Table 4, the variances in $e_{9}-e_{29}$ all decrease. This is due to the fact that the negative impact brought by the decrease in the weights cannot overcome the positive impact brought by increasing the capacity of $e_{23}$ . However, if the negative impact surpasses the positive impact, then the variance will increase, which may happen in a subset of networks.

Finally, we verify the findings in Remark 6.5-(iv). We focus on Cases 2-3 with $\delta_{ij}^{*}=0$ for all the lines that are not changed by either constructing new lines or increasing the line capacity. The cycle-cluster $\{e_{1},e_{2},e_{3},e_{4}\}$ includes a cycle. The basis vector corresponding to this cycle is $\bm{\xi}_{1}=[-1,-1,-1,1,0,\cdots,0]^{\bm{\top}}$ . By scaling this vector to unit length, we obtain $\bm{X}_{1}=[-1/2,-1/2,-1/2,1/2,0,\cdots,0]^{\bm{\top}}$ . From formula (38), we obtain that the diagonal elements $\bm{Q}_{\delta}$ at positions (1-4) are all $3/80$ , which is consistent with the values shown in Table 4. Hence, the construction of $e_{4}$ decreases the variances of the phase angle differences, and the size of the decrease depends on the length of the cycle. It is verified that the variances in $e_{1}-e_{4}$ in cases 5 and 6 can also be calculated by (46) for simplicity. Let us next focus on the conservative estimation of the variances in the lines in a cycle-cluster by formula (46). For example, the variance in $e_{29}$ in case 2 can be approximated as $0.0333$ for simplicity from formula (46) by taking $\mathcal{E}_{c}=\{e_{22},e_{28},e_{29}\}$ . This value is larger than $0.0326$ as shown in Table 4. Because constructing new lines to form cycles or increasing the capacities of lines changes the power flows, which may decrease the weights of the lines in the cycle-cluster or even destroy the synchronization, it is complicated to analyse how the variances of the lines of this cycle-cluster change. However, in a real network, the phase angle differences are usually small, and the weight $l_{c_{ij}}\approx l_{ij}$ , which is often assumed in the investigation of the synchronization of power systems (Poolla et al., 2017; Tegling et al., 2015). In this case, the negative influences on the weight can be neglected and the variances decrease if new lines are constructed to form small cycles or the capacities of the lines are increased. The reduction in the variances can be approximated using (46).

In regard to the bounds of the variance matrices for the networks with non-uniform disturbance-damping ratio, it is shown for case 7 in Table 2 and 4 that the variances of the frequencies at the nodes and the phase angle differences in the lines are both constrained by the lower bound and the upper bound in (43) and (44) respectively. For the frequency, it is demonstrated that the variance at the node which possesses the largest disturbance-damping ratio is closer to the upper bound and that at the node with the smallest disturbance-damping ratio is closer to the lower bound. For example, the variance at node $5$ , which has the largest disturbance-damping ratio $\eta_{5}=\sqrt{5}$ , is $0.1901$ which is closer to the upper bound $\sqrt{5}/10=0.2236$ . However, those at the nodes $1,6-23$ with the smallest disturbance-damping ratio are all closer to the lower bound $1/2$ . For the phase angle differences, the variance in the lines which connect nodes with larger disturbance-damping ratio are usually larger. For example, the variance in $e_{8}$ is $0.0904$ which becomes closer to the upper bound $0.1198$ compared with its value in case 6. However, the variances in the lines which are far away from the nodes with larger disturbance-damping ratio are closer to the lower bounds. This is seen from the variance in lines $e_{26}-e_{29}$ .

In regard to the vulnerable nodes, it is found in Table 2 that nodes $1,6-23$ in case 1-6 and node $6$ in case 7 are the most vulnerable nodes. The remedy methods include increasing the inertia and decreasing the disturbance-damping ratio at these nodes or their neighbor nodes. With respect to the vulnerable lines, it is seen in Table 4 for cases 1-7 that the single lines are usually the vulnerable lines, for example, lines $e_{5}-e_{8}$ . The remedy method includes increasing the capacities of these lines and constructing new lines to include these lines into cycles.

8 Conclusion

In this paper, based on a stochastic Gaussian system, we have investigated the dependence of the fluctuations in a power system on system parameters when subjected to stochastic disturbances. The dynamics of turbine-governors of the synchronous machines and that of voltage may be considered in the system (Trip et al., 2019). By the method proposed in this paper, the impact of the system parameters on the fluctuations of the frequency and voltage at each node and the phase angle difference in each line can be investigated. In that case, the system parameters includes the ones in the dynamics of the turbine-governor and voltage besides those studied in this paper.

A future investigation will address the deduction of explicit formulas for the variance matrices of the frequencies at the nodes and the phase angle differences in the lines in the network with non-uniform disturbance-damping ratio and lossy transmission lines.

References

(1)
Auer et al. (2017) Auer, S., Hellmann, F., Krause, M. and Kurths, J. (2017). Stability of synchrony against local intermittent fluctuations in tree-like power grids, Chaos 27(12): 127003.
Baillieul and Byrnes (1982) Baillieul, J. and Byrnes, C. (1982). Geometric critical point analysis of lossless power system models, IEEE Trans. Circuits Syst. 29(11): 724–737.
Biggs (1993) Biggs, N. (1993). Algebraic Graph Theory, 2nd edn, Cambridge University Press.
Chiang et al. (1988) Chiang, H. D., Hirsch, M. and Wu, F. (1988). Stability regions of nonlinear autonomous dynamical systems, IEEE Trans. Autom. Control 33(1): 16–27.
Diestel (2000) Diestel, R. (2000). Graph Theory, Springer-Verlag, New York 1997.
Dörfler and Bullo (2012) Dörfler, F. and Bullo, F. (2012). Synchronization and transient stability in power networks and nonuniform Kuramoto oscillators, SIAM J. Control Optim. 50(3): 1616–1642.
Dörfler and Bullo (2014) Dörfler, F. and Bullo, F. (2014). Synchronization in complex networks of phase oscillators: A survey, Automatica 50(6): 1539 – 1564.
Doyle et al. (1989) Doyle, J. C., Glover, K., Khargonekar, P. P. and Francis, B. A. (1989). State-space solutions to standard H2 and H infinty control problems, IEEE Trans. Autom. Control 34(8): 831–847.
Fazlyab et al. (2017) Fazlyab, M., Dörfler, F. and Preciado, V. M. (2017). Optimal network design for synchronization of coupled oscillators, Automatica 84: 181 – 189.
Haehne et al. (2019) Haehne, H., Schmietendorf, K., Tamrakar, S.and Peinke, J. and Kettemann, S. (2019). Propagation of wind-power-induced fluctuations in power grids, Phys. Rev. E 99: 050301.
Horn and Johnson (2013) Horn, R. A. and Johnson, C. R. (2013). Matrix Analysis, 2nd edn, Cambridge University Press, Cambridge, UK.
Jafarpour et al. (2022) Jafarpour, S., Huang, E. Y., Smith, K. D. and Bullo, F. (2022). Flow and elastic networks on the $n$ -torus: Geometry, analysis, and computation, SIAM Review 64(1): 59–104.
Karatzas and Shreve (1988) Karatzas, I. and Shreve, S. (1988). Brownian motion and stochastic calculus, 2nd edn, Springer-Verlag, Berlin.
Kettemann (2016) Kettemann, S. (2016). Delocalization of disturbances and the stability of ac electricity grids, Phys. Rev. E 94: 062311.
Kundur (1994) Kundur, P. (1994). Power system stability and control, McGraw-Hill, New York.
Kwakernaak and Sivan (1972) Kwakernaak, H. and Sivan, R. (1972). Linear optimal control systems, Wiley-Interscience, New York.
Luxemburg and Huang (1987) Luxemburg, L. A. and Huang, G. (1987). On the number of unstable equilibria of a class of nonlinear systems, 26th IEEE Conf. Decision Control, Vol. 20, IEEE, pp. 889–894.
Manik et al. (2017) Manik, D., Rohden, M., Ronellenfitsch, H.and Zhang, X., Hallerberg, S., Witthaut, D. and Timme, M. (2017). Network susceptibilities: Theory and applications, Phys. Rev. E 95: 012319.
Marris (2008) Marris, E. (2008). Energy: Upgrading the Grid, Nature 454: 570–573.
Menck et al. (2014) Menck, P. J., Heitzig, J., Kurths, J. and Joachim Schellnhuber, H. (2014). How dead ends undermine power grid stability, Nat. Commun. 5: 3969.
Menck et al. (2013) Menck, P. J., Heitzig, J., Marwan, N. and Kurths, J. (2013). How basin stability complements the linear-stability paradigm, Nat. Phys. 9(2): 89–92.
Motter et al. (2013) Motter, A. E., Myers, S. A., Anghel, M. and Nishikawa, T. (2013). Spontaneous synchrony in power-grid networks, Nat. Phys. 9(3): 191–197.
Nishikawa et al. (2015) Nishikawa, T., Molnar, F. and Motter, A. E. (2015). Stability landscape of power-grid synchronization, IFAC-PapersOnLine 48(18): 1 – 6.
Nishikawa and Motter (2015) Nishikawa, T. and Motter, A. E. (2015). Comparative analysis of existing models for power-grid synchronization, New J. Phys. 17(1): 015012.
Poolla et al. (2017) Poolla, B. K., Bolognani, S. and Dörfler, F. (2017). Optimal placement of virtual inertia in power grids, IEEE Trans. Autom. Control 62(12): 6209–6220.
Schäfer et al. (2018) Schäfer, B., Beck, C., Aihara, K., Witthaut, D. and Timme, M. (2018). Non-gaussian power grid frequency fluctuations characterized by lévy-stable laws and superstatistics, Nat. Energy 3(2): 119–126.
Schmietendorf et al. (2017) Schmietendorf, K., Peinke, J. and Kamps, O. (2017). The impact of turbulent renewable energy production on power grid stability and quality, Eur. Phys. J. B 90(11): 1–6.
Simpson-Porco et al. (2016) Simpson-Porco, J. W., Dörfler, F. and Bullo, F. (2016). Voltage collapse in complex power grids, Nat. Commun. 7: 10790.
Skar (1980) Skar, S. J. (1980). Stability of multi-machine power systems with nontrivial transfer conductances, SIAM J. Appl. Math. 39(3): 475–491.
Tegling et al. (2015) Tegling, E., Bamieh, B. and Gayme, D. F. (2015). The price of synchrony: Evaluating the resistive losses in synchronizing power networks, IEEE Trans. Control Netw. Syst. 2(3): 254–266.
Toscano (2013) Toscano, R. (2013). Structured controllers for uncertain systems, Springer-verlag, London.
Trip et al. (2020) Trip, S., C., M., Persis, C. D., Ferrara, A. and Scherpen, J. M. A. (2020). Robust load frequency control of nonlinear power networks, Int. J. Control 93(2): 346–359.
Trip et al. (2019) Trip, S., Cucuzzella, M., De Persis, C., van der Schaft, A. and Ferrara, A. (2019). Passivity-based design of sliding modes for optimal load frequency control, IEEE Trans. Control Syst. Technol. 27(5): 1893–1906.
Wang and Crow (2013) Wang, K. and Crow, M. L. (2013). The Fokker-Planck equation for power system stability probability density function evolution, IEEE Trans. Power Syst. 28(3): 2994–3001.
Wolff et al. (2019) Wolff, M. F., Schmietendorf, K., Lind, P. G., Kamps, O., Peinke, J. and Maass, P. (2019). Heterogeneities in electricity grids strongly enhance non-gaussian features of frequency fluctuations under stochastic power input, Chaos 29(10): 103149.
Xi et al. (2017) Xi, K., Dubbeldam, J. L. A. and Lin, H. X. (2017). Synchronization of cyclic power grids: equilibria and stability of the synchronous state, Chaos 27(1): 013109.
Xi et al. (2020) Xi, K., Lin, H. X., Shen, C. and Van Schuppen, J. H. (2020). Multilevel power-imbalance allocation control for secondary frequency control of power systems, IEEE Trans. Autom. Control 65(7): 2913–2928.
Xie et al. (2011) Xie, L., Carvalho, P. M. S., Ferreira, L. A. F. M., Liu, J., Krogh, B. H., Popli, N. and Ilić, M. D. (2011). Wind integration in power systems: Operational challenges and possible solutions, Proceedings of the IEEE 99(1): 214–232.
Ye et al. (2016) Ye, H., Liu, Y. and Zhang, P. (2016). Efficient eigen-analysis for large delayed cyber-physical power system using explicit infinitesimal generator discretization, IEEE Trans. Power Syst. 31(3): 2361–2370.
Zaborsky et al. (1985) Zaborsky, J., Huang, G., Leung, T. C. and Zheng, B. (1985). Stability monitoring on the large electric power system, 24th IEEE Conf. Decision Control, Vol. 24, IEEE, pp. 787–798.
Zaborszky et al. (1988) Zaborszky, J., Huang, G., Zheng, B. and Leung, T. C. (1988). On the phase portrait of a class of large nonlinear dynamic systems such as the power system, IEEE Trans. Autom. Control 33(1): 4–15.
Zhang et al. (2019) Zhang, X., Hallerberg, S., Matthiae, M., Witthaut, D. and Timme, M. (2019). Fluctuation-induced distributed resonances in oscillatory networks, Sci. Adv. 5(7): eaav1027.
Zhang et al. (2020) Zhang, X., Witthaut, D. and Timme, M. (2020). Topological determinants of perturbation spreading in networks, Phys. Rev. Lett. 125: 218301.

Appendix A Appendix

A.1 The variance matrix of a linear stochastic process

Definition A.1.

Consider a linear stochastic system

	$\displaystyle\text{\emph{d}}\bm{x}(t)$	$\displaystyle=\bm{A}\bm{x}(t)\text{\emph{d}}t+\bm{B}\text{\emph{d}}\bm{\mu}(t),\bm{x}(0)=\bm{x}_{0},$
	$\displaystyle\bm{y}$	$\displaystyle=\bm{C}\bm{x}(t),$

where $\bm{x}\in\mathbb{R}^{n}$ , $\bm{A}\in\mathbb{R}^{n\times n}$ , $\bm{B}\in\mathbb{R}^{n\times m}$ , $\bm{x}_{0}\in G(\bm{0},\bm{Q}_{x_{0}})$ is a Gaussian random variable where $\bm{Q}_{x_{0}}\in\mathbb{R}^{n\times n}$ is the variance matrix of $\bm{x}_{0}$ , $\bm{C}\in\mathbb{R}^{z\times n}$ , $\bm{\mu}\in\mathbb{R}^{m}$ is standard Brownian motion, $\bm{y}\in\mathbb{R}^{z}$ is the output. It follows from (Kwakernaak and Sivan, 1972, Theorem 1.52) that the state $\bm{x}$ and $\bm{y}$ are Gaussian process, i.e., for all $t>0$ ,

\displaystyle\bm{x}(t)\in G(\bm{0},\bm{Q}_{x,tv}(t)),~~\bm{y}(t)\in G(\bm{0},\bm{Q}_{y,tv}(t))

where the variance matrix $\bm{Q}_{x,tv}(t)\in\mathbb{R}^{n\times n}$ of $\bm{x}(t)$ is

\displaystyle\bm{Q}_{x,tv}(t)=\emph{e}^{\bm{A}t}\bm{Q}_{x_{0}}\emph{e}^{\bm{A}^{\top}t}+\int_{0}^{t}{\emph{e}^{\bm{A}\tau}\bm{B}\bm{B}^{\bm{\top}}\emph{e}^{\bm{A}^{\bm{\top}}\tau}\text{\emph{d}}\tau}

and the variance matrix $\bm{Q}_{y,tv}(t)\in\mathbb{R}^{z\times z}$ of $\bm{y}(t)$ satisfies $\bm{Q}_{y,tv}(t)=\bm{C}\bm{Q}_{x,tv}(t)\bm{C}^{\top}$ . The matrix $\bm{Q}_{x,tv}(t)$ satisfies the matrix differential equation


$\displaystyle\dot{\bm{Q}}_{x,tv}(t)$	$\displaystyle=\bm{A}\bm{Q}_{x,tv}(t)+\bm{Q}_{x,tv}(t)\bm{A}^{\top}+\bm{B}\bm{B}^{\top},$	(47a)
$\displaystyle\bm{Q}_{x,tv}(0)$	$\displaystyle=\bm{Q}_{x_{0}}.$	(47b)

In addition, if $\bm{A}$ is Hurwitz, then there exists an invariant distribution of the stochastic processes $\bm{x}(t)$ and $\bm{y}(t)$ with asymptotic variance matrices

\displaystyle\bm{Q}_{x}=\lim_{t\rightarrow\infty}\bm{Q}_{x,tv}(t)=\int_{0}^{+\infty}{\emph{e}^{\bm{A}\tau}\bm{B}\bm{B}^{\bm{\top}}\emph{e}^{\bm{A}^{\bm{\top}}\tau}\text{\emph{d}}\tau},

and $\bm{Q}_{y}=\bm{C}\bm{Q}_{x}\bm{C}^{\top}$ . The matrix $\bm{Q}_{x}$ , which is called the controllability Gramian of the pair $(\bm{A},\bm{B})$ , is the unique solution of the Lyapunov equation due to the Hurwitz condition (Toscano, 2013; Doyle et al., 1989),

\displaystyle\bm{A}\bm{Q}_{x}+\bm{Q}_{x}\bm{A}^{\bm{\top}}+\bm{B}\bm{B}^{\bm{\top}}=\bm{0}.

(48)

which can be either derived from the limit of the differential equation (47a) or from

	$\displaystyle\bm{A}\bm{Q}_{x}+\bm{Q}_{x}\bm{A}^{\bm{\top}}$	$\displaystyle=\int_{0}^{+\infty}\frac{\emph{d}}{\emph{d}t}(\emph{e}^{\bm{A}t}\bm{B}\bm{B}^{\bm{\top}}\emph{e}^{\bm{A}^{\bm{\top}}t})\emph{d}t$
		$\displaystyle=(\emph{e}^{\bm{A}t}\bm{B}\bm{B}^{\bm{\top}}\emph{e}^{\bm{A}^{\bm{\top}}t})\Big{\|}_{0}^{+\infty}=-\bm{B}\bm{B}^{\top}.$

A.2 The Moore-Penrose Pseudo inverse of real symmetric matrices

Theorem A.2.

Consider a real symmetric matrix $\bm{S}\in\mathbb{R}^{n\times n}$ . There exists an orthogonal matrix $\bm{V}\in\mathbb{R}^{n\times n}$ such that

\displaystyle\bm{V}^{\bm{\top}}\bm{S}\bm{V}=\bm{\Sigma}

where $\bm{\Sigma}=\text{diag}(\sigma_{i})\in\mathbb{R}^{n\times n}$ is a diagonal matrix with the diagonal elements $\sigma_{i}$ being the eigenvalues of $\bm{S}$ , the column vectors $\bm{v}_{i}$ of $\bm{V}$ are orthonormal eigenvectors of $\bm{S}$ corresponding to the eigenvalue $\sigma_{i}$ . In addition, the Moore-Penrose pseudo inverse is defined by the formula

\displaystyle\bm{S}^{{\dagger}}=\bm{V}\Sigma^{{\dagger}}\bm{V}^{\bm{\top}}=\sum_{i=1}^{n}\sigma_{i}^{*}\bm{v}_{i}\bm{v}_{i}^{\bm{\top}}

where $\bm{\Sigma}^{\dagger}=\text{diag}(\sigma_{i}^{*})\in\mathbb{R}^{n\times n}$ with $\sigma_{i}^{*}=1/\sigma_{i}$ if $\sigma_{i}\neq 0$ , otherwise $\sigma_{i}^{*}=0$ (Horn and Johnson, 2013).

A.3 The basis vectors of the kernel of $\widetilde{\bm{C}}\bm{R}^{1/2}$

The cycle space of a graph is defined as the kernel of the incidence matrix $\widetilde{\bm{C}}$ , which is a vector subspace in $\mathbb{R}^{m}$ . By graph theory, we have $\text{rank}(\widetilde{\bm{C}})=n-1$ . Hence, the dimension of the cycle space is $m-n+1$ . It is obvious that the cycle space of an acyclic graph is an empty space. For a graph with cycles, the basis for the cycle space is derived by the following method: Considering a cycle $\mathcal{C}$ with a set $\mathcal{E}_{c}$ of edges in the graph $\mathcal{G}$ , we specify a direction for $\mathcal{C}$ ; then, the vector $\bm{\xi}_{c}=[\xi_{c,1},\xi_{c,2},\cdots,\xi_{c,m}]^{\bm{\top}}\in\mathbb{R}^{m}$ such that

\displaystyle\xi_{c,k}=\begin{cases}+1,&\text{\footnotesize if $e_{k}\in\mathcal{E}_{c}$ with direction $=$ the cycle direction},\\ -1,&\text{\footnotesize if $e_{k}\in\mathcal{E}_{c}$ with direction $\neq$ the cycle direction},\\ 0,&\text{\footnotesize otherwise.}\end{cases}

belongs to the kernel of $\widetilde{\bm{C}}$ such that $\widetilde{\bm{C}}\bm{\xi}_{c}=\bm{0}$ (Biggs, 1993). The basis for the cycle space can be derived by taking the vectors as $\bm{\xi}_{c}$ for $c=1,\cdots,m-n+1$ corresponding to the $(m-n+1)$ fundamental cycles (Diestel, 2000, Theorem 1.9.6) in the graph. Because $\bm{R}$ is non-singular, the vectors $\bm{R}^{-1/2}\bm{\xi}_{c}$ for all the cycles are the basis vectors of the kernel of $\widetilde{\bm{C}}\bm{R}^{1/2}$ . The orthonormal basis vectors $\bm{X}_{i}$ are obtained by Gram-Schmidt orthogonalization of the basis vectors $\bm{R}^{-1/2}\bm{\xi}_{c}$ . The fundamental cycles can be obtained by the following method. Let $\mathcal{T}$ be a spanning tree of the graph $\mathcal{G}$ . Then $\mathcal{T}$ has $n-1$ edges and there are $m-n+1$ edges of $\mathcal{G}$ lie outside of $\mathcal{T}$ . Then for each of these $m-n+1$ edges $e\in\mathcal{E}\setminus\mathcal{E}(\mathcal{T})$ , the graph $\mathcal{T}+e$ contains a cycle, which is a fundamental cycle. Note that the basis vectors of the cycle space may not be unique due to the non-uniqueness of the spanning tree.

Table 1: Table describing the 7 cases in the example; line

e_{4}

is present in the network if the label is

1

and not present if the label is

0

Case	Line	Source	Load	Line $l_{ij}$		$b_{i}^{2}/d_{i}$
	$e_{4}$	$P_{i}$	$P_{i}$	$e_{23}$	others	1-5	others
1	0	0	0	10	10	1	1
2	1	0	0	10	10	1	1
3	1	0	0	20	10	1	1
4	0	3.6	-1	10	10	1	1
5	1	3.6	-1	10	10	1	1
6	1	3.6	-1	20	10	1	1
7	1	3.6	-1	20	10	$\sqrt{i}$	1

Table 2: The variances of the frequencies in the 7 cases in the example; 7L and 7U denotes the lower and upper bounds in case 7;

\sqrt{5}\approx 2.2236

case	1	2	3	4	5	6	7	8	9	10	11	12
1-6, 7L	$1/2$	$1/4$	$1/6$	$1/8$	$1/10$	$1/2$	$1/2$	$1/2$	$1/2$	$1/2$	$1/2$	$1/2$
7	$0.5002$	$0.3012$	$0.2600$	$0.2275$	$0.1901$	$0.5107$	$0.5024$	$0.5054$	$0.5034$	$0.5002$	$0.5006$	$0.5038$
7U	$\sqrt{5}/2$	$\sqrt{5}/4$	$\sqrt{5}/6$	$\sqrt{5}/8$	$\sqrt{5}/10$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$
case	13	14	15	16	17	18	19	20	21	22	23
1-6, 7L	$1/2$	$1/2$	$1/2$	$1/2$	$1/2$	$1/2$	$1/2$	$1/2$	$1/2$	$1/2$	$1/2$
7	$0.5025$	$0.5002$	$0.5000$	$0.5000$	$0.5000$	$0.5000$	$0.5001$	$0.5004$	$0.5042$	$0.5003$	$0.5001$
7U	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$	$\sqrt{5}/2$

Table 3: The weights of the lines of the network in Cases (4-7) in the example.

Case	$e_{1}$	$e_{2}$	$e_{3}$	$e_{4}$	$e_{5}$	$e_{6}$	$e_{7}$	$e_{8}$	$e_{9}$	$e_{10}$	$e_{11}$	$e_{12}$	$e_{13}$	$e_{14}$	$e_{15}$
4	$9.9499$	$7.8460$	$9.3295$	$-$	$9.9499$	$9.6561$	$9.8712$	$9.3295$	$9.9109$	$9.9945$	$9.8215$	$9.9193$	$9.7394$	$9.9549$	$9.9860$
5	$9.8530$	$9.3706$	$9.9602$	$9.6262$	$9.9499$	$9.6561$	$9.8712$	$9.3295$	$9.9109$	$9.9945$	$9.8215$	$9.9193$	$9.7394$	$9.9549$	$9.9860$
6	$9.8530$	$9.3706$	$9.9602$	$9.6262$	$9.9499$	$9.6561$	$9.8712$	$9.3295$	$9.9092$	$9.9941$	$9.8312$	$9.9209$	$9.7423$	$9.9607$	$9.9897$
Case	$e_{16}$	$e_{17}$	$e_{18}$	$e_{19}$	$e_{20}$	$e_{21}$	$e_{22}$	$e_{23}$	$e_{24}$	$e_{25}$	$e_{26}$	$e_{27}$	$e_{28}$	$e_{29}$
4	$9.7337$	$9.9540$	$9.9086$	$9.7745$	$9.6346$	$9.9380$	$9.8829$	$9.4709$	$9.7528$	$9.9266$	$9.9978$	$9.9687$	$9.9965$	$9.9198$
5	$9.7337$	$9.9540$	$9.9086$	$9.7745$	$9.6346$	$9.9380$	$9.8829$	$9.4709$	$9.7528$	$9.9266$	$9.9978$	$9.9687$	$9.9965$	$9.9198$
6-7	$9.7430$	$9.9434$	$9.9271$	$9.7973$	$9.6722$	$9.9495$	$9.9069$	$19.6990$	$9.6934$	$9.8933$	$9.9896$	$9.9852$	$9.9984$	$9.9300$

Table 4: The variances of the phase angle differences in the 7 Cases in the example; 7L and 7U denotes the lower and upper bounds in case 7.

Case	$e_{1}$	$e_{2}$	$e_{3}$	$e_{4}$	$e_{5}$	$e_{6}$	$e_{7}$	$e_{8}$	$e_{9}$	$e_{10}$	$e_{11}$	$e_{12}$	$e_{13}$	$e_{14}$	$e_{15}$
1	$0.0500$	$0.0500$	$0.0500$	$-$	$0.050$	$0.050$	$0.050$	$0.050$	$0.0394$	$0.0394$	$0.0298$	$0.0394$	$0.0394$	$0.0340$	$0.0281$
2	$0.0375$	$0.0375$	$0.0375$	$0.0375$	$0.050$	$0.050$	$0.050$	$0.050$	$0.0394$	$0.0394$	$0.0298$	$0.0394$	$0.0394$	$0.0340$	$0.0281$
3	$0.0375$	$0.0375$	$0.0375$	$0.0375$	$0.050$	$0.050$	$0.050$	$0.050$	$0.0394$	$0.0394$	$0.0297$	$0.0394$	$0.0394$	$0.0339$	$0.0281$
4	$0.0503$	$0.0637$	$0.0536$	$-$	$0.0503$	$0.0518$	$0.0507$	$0.0536$	$0.0397$	$0.0395$	$0.0302$	$0.0397$	$0.0403$	$0.0343$	$0.0283$
5	$0.0383$	$0.0396$	$0.0380$	$0.0389$	$0.0503$	$0.0518$	$0.0507$	$0.0536$	$0.0397$	$0.0395$	$0.0302$	$0.0397$	$0.0403$	$0.0343$	$0.0283$
6, 7L	$0.0383$	$0.0396$	$0.0380$	$0.0389$	$0.0503$	$0.0518$	$0.0507$	$0.0536$	$0.0397$	$0.0395$	$0.0301$	$0.0397$	$0.0402$	$0.0342$	$0.0283$
7	$0.0397$	$0.0444$	$0.0518$	$0.0530$	$0.0551$	$0.0566$	$0.0524$	$0.0904$	$0.0398$	$0.0397$	$0.0304$	$0.0398$	$0.0403$	$0.0354$	$0.0290$
7U	$0.0856$	$0.0884$	0.0849	$0.0869$	$0.1124$	$0.1158$	$0.1133$	$0.1198$	$0.0889$	$0.0883$	$0.0674$	$0.0888$	$0.0900$	$0.0765$	$0.0632$
Case	$e_{16}$	$e_{17}$	$e_{18}$	$e_{19}$	$e_{20}$	$e_{21}$	$e_{22}$	$e_{23}$	$e_{24}$	$e_{25}$	$e_{26}$	$e_{27}$	$e_{28}$	$e_{29}$
1	$0.0262$	$0.0304$	$0.0292$	$0.0352$	$0.0343$	$0.0352$	$0.0302$	$0.0430$	$0.0430$	$0.0430$	$0.0430$	$0.0430$	$0.0326$	$0.0326$
2	$0.0262$	$0.0304$	$0.0292$	$0.0352$	$0.0343$	$0.0352$	$0.0302$	$0.0430$	$0.0430$	$0.0430$	$0.0430$	$0.0430$	$0.0326$	$0.0326$
3	$0.0262$	$0.0303$	$0.0291$	$0.0351$	$0.0341$	$0.0351$	$0.0300$	$0.0231$	$0.0424$	$0.0424$	$0.0424$	$0.0424$	$0.0325$	$0.0325$
4	$0.0267$	$0.0306$	$0.0296$	$0.0359$	$0.0353$	$0.0356$	$0.0305$	$0.0451$	$0.0440$	$0.0434$	$0.0431$	$0.0432$	$0.0327$	$0.0328$
5	$0.0267$	$0.0306$	$0.0296$	$0.0359$	$0.0353$	$0.0356$	$0.0305$	$0.0451$	$0.0440$	$0.0434$	$0.0431$	$0.0432$	$0.0327$	$0.0328$
6, 7L	$0.0267$	$0.0305$	$0.0293$	$0.0357$	$0.0350$	$0.0354$	$0.0302$	$0.0235$	$0.0436$	$0.0429$	$0.0426$	$0.0426$	$0.0326$	$0.0327$
7	$0.0272$	$0.0400$	$0.0383$	$0.0365$	$0.0352$	$0.0356$	$0.0303$	$0.0318$	$0.0456$	$0.0432$	$0.0426$	$0.0426$	$0.0326$	$0.0328$
7U	$0.0597$	$0.0683$	$0.0656$	$0.0799$	$0.0783$	$0.0792$	$0.0676$	$0.0525$	$0.0976$	$0.0960$	$0.0952$	$0.0952$	$0.0729$	$0.0731$

Synchronization of Power Systems under Stochastic Disturbances