Finite-time correlations boost large voltage angle fluctuations in electric power grids

The operational state of an AC power grid is determined by complex voltages $V_i = |V_i| \exp[i \phi_i]$ at each of the $i = 1, \ldots N$ nodes of the grid. In normal operation, voltage amplitudes are fixed not far from their rated value, and voltage angles rotate close to the rated frequency, $\phi_i(t) = \Omega_0 t + \theta_i(t)$, with $\Omega_0/2 \pi = 50$ or 60 Hz. Over time intervals ranging roughly from seconds to several tens of seconds, the transient dynamics of high-voltage power grids is given by the swing equations [5]. They govern the time-evolution of voltage angles in a frame rotating at the rated frequency. In high-voltage power grids, a standard approximation is the lossless line approximation, which neglects Ohmic losses. The swing equations then read

$$\begin{equation} m_i \ddot{\theta}_i + d_i \dot{\theta}_i = P_i - \sum_{j}B_{ij}\sin(\theta_i-\theta_j) \, , \end{equation} \tag{ 1 }$$

with the inertia (m_i ) and damping (d_i ) parameters. The active power P_i is positive for generators and negative for loads, and $B_{ij} = b_{ij} |V_i| |V_j|$ denotes the product of the voltage magnitudes at nodes i and j with the line susceptance. If there is no line between i and j , then $b_{ij} = B_{ij} = 0$ . In the lossless line approximation, line conductances are neglected. The approximation is justified when dealing with very high-voltage power grids, which typically have $g/b \lt 0.1$ for the ratio of conductance over susceptance (see supplemental material).

At equilibrium, electric power grids are synchronized network systems [29]. They lie close to an operational synchronous state where voltage angles are solutions to a set of transcendental equations called the power flow equations. Under our assumptions, they read

$$\begin{equation} P_i = \sum_{j}B_{ij}\sin(\theta_i-\theta_j) \, . \end{equation} \tag{ 2 }$$

Their solution $\{\theta^{(0)}_i\}_{i = 1, \ldots N}$ corresponds to the instantaneous, synchronous operational state of the power grid.

We want to investigate how a local perturbation about the solution to equation (2) propagates across the system and influences voltage angles far away from it. Such an occurrence is illustrated for the PanTaGruEl [30] model of the synchronous grid of continental Europe in figure 1. Initially, the system is in a steady-state solution of equation (2). An abrupt power loss $P_i \rightarrow 0$, corresponding to the disconnection of a large power plant in Spain, brings the system out of equilibrium. Following that perturbation, a voltage angle wave propagates across the grid, which is represented in five consecutive color-coded snapshots in figure 1.

Zoom In Zoom Out Reset image size

In large power grids, even the loss of large power plants is a relatively weak perturbation in a mathematical sense. For instance the European Network of Transmission System Operators for Electricity (ENTSO-E) reference incident considers the simultaneous tripping of two of the largest power plants, connected to the same bus [31]. This corresponds to less than one percent of the total power injected in the synchronous grid of continental Europe. For the case plotted in figure 1 of a power loss of $\Delta P = 900$ MW, corresponding to a large power plant, frequency deviations never exceed $\Delta \omega/2 \pi = 0.12$ Hz, i.e. a fraction of a percent of the rated frequency [17]. Power feed-in noises being by nature smaller than the rated power on which they are superimposed, it is therefore legitimate to investigate them through the linearization of equation (1) about the operational synchronous state. With $\theta_i = \theta^{(0)}_i + \delta\theta_i$ and $P_i = P_i^{(0)} + \delta P_i$, one gets

$$\begin{equation} \mathbf{M} \delta \ddot{\boldsymbol{\theta}} + \mathbf{D} \delta \dot{\boldsymbol{\theta}} = \delta \mathbf{P} - \mathbf{L} \delta {\boldsymbol{\theta}} \, , \end{equation} \tag{ 3 }$$

where we grouped the voltage angle deviations into a vector $\delta {\boldsymbol{\theta}}$, and introduced the diagonal inertia and damping matrices, $\mathbf{M} = {\textrm{diag}}(m_i)$ (with $m_i = 0$ on load nodes), $\mathbf{D} = {\textrm{diag}}(d_i)$ as well as the weighted network Laplacian matrix L,

$$\begin{equation} \mathbf{L}_{ij} = \begin{cases} -B_{ij} \cos(\theta_i^{(0)} - \theta_j^{(0)}) \, , & \text{for } i \neq j \, , \\ \sum_{k}B_{ik}\cos(\theta_i^{(0)}-\theta_k^{(0)}) \, , & \text{for } i = j \, . \end{cases} \end{equation} \tag{ 4 }$$

The perturbation generating a wave of voltage angle and frequency disturbances is encoded in the source term vector $\delta \mathbf{P}$, whose components are non-zero at nodes where the perturbation is active. Below we consider cases of (a) a single noisy perturbation and (b) a collection of independent, geographically distributed noisy perturbations.

Power grids have two types of nodes, corresponding to power plants and loads. They have very different dynamical inertia and damping parameters. Most loads as well as inverter-connected, new renewable sources of energy have no inertia, $m_i = 0$, while traditional power plants have an inertia roughly proportional to their rated power output [5]. Furthermore, loads have a damping parameter significantly smaller than generators [32]. While it is crucial to incorporate these dynamic inhomogeneities in any analysis of realistic power grids, they render analytical approaches intractable. Recent works took a perturbation theory approach to incorporate small deviations about the homogeneous case [33, 34], however most are based on homogeneity assumptions [35–38]. To justify it, one often invokes a Kron-reduction [39] into an effective network with modified line susceptances connecting only inertiaful, generator nodes. This transformation is based on Schur's complement formula [40], and since the reduced load nodes have no inertia and a much smaller damping term, this reduction modifies the dynamics on the generators only marginally. Once the reduction is performed, one furthermore argues that considering uniform damping and inertia, $d_i = d$, $m_i = m$ is justified, because, only large plants, all with large rated power are connected to the high-voltage grids we are focusing on here. Additionally, machine measurements indicate that the ratio of damping over inertia does not vary by much from one machine to another [41]. Hence, in our analytical treatment, we consider noise propagation from equation (3) for a Kron-reduced network with homogeneous dynamic parameters, $d_i = d$, $m_i = m$. However, our numerical simulations are entirely free from this assumption and are based on a realistic, inhomogeneous power grid model. Our numerical results corroborate our analytical results in the particularly relevant regime of long noise correlation time.

Equation (3) is a damped wave equation with a source term. It is defined on a discrete, meshed complex network encoded in the Laplacian matrix L, which accordingly replaces the Laplace operator of continuous wave equations. Given a source term $\delta \mathbf{P}(t)$, we compute the moments $\mu_p \equiv \langle \delta \tilde\theta_i^p(t \rightarrow \infty) \rangle$, $p \leqslant 4$ of the distribution of angle deviations at any node i on the network. Because we are interested in local fluctuations about the instantaneous average response, we measure angle deviations from their geographical average, $\delta \tilde\theta_i(t\rightarrow \infty) = \delta \theta_i(t\rightarrow \infty) - \Delta$ where $\Delta = N^{-1} \sum_j \delta \theta_j(t \rightarrow \infty)$. Here, $t\rightarrow \infty$ means that the observation takes place long after the onset of the noisy perturbation, to avoid transient behaviors. To shorten notations, we do not explicitly write it from now on.

To calculate $\mu_p \equiv \langle \delta \tilde\theta_i^p \rangle$, we use a modal expansion of equation (3) over the set of eigenmodes $\{\mathbf{u}_{\alpha} \}$ of the Laplacian matrix L. We first write the total angle deviation as $\delta\theta_i(t) = \sum_\alpha c_\alpha(t) \, u_{\alpha,i}$ . Equation (3) then gives

$$\begin{eqnarray} &&m\,\ddot{c}_\alpha + d\,\dot{c}_\alpha + \lambda_\alpha c_\alpha = \delta \mathbf{P}(t) \cdot \mathbf{u}_{\alpha} \,, \end{eqnarray} \tag{ 5 }$$

where $\mathbf{L} \textbf{u}_{\alpha} = \lambda_\alpha \textbf{u}_{\alpha}$, with $\lambda_\alpha\geqslant 0$, $\alpha = 1, \ldots N$. Because angle deviations $\delta \tilde\theta_i(t) = \delta \theta_i(t)-\Delta(t)$ are defined relative to their geographical average $\Delta(t)$, the bulk voltage angle shift along the zero-mode $u_{1,i} = 1/\sqrt{N}$ of the Laplacian matrix is subtracted, and the rest of this article considers the dynamics orthogonal to that zero-mode. This in particular removes contributions from time-dependent excursions of angle/frequency averages, which may be another source of non-Gaussianity [12].

Equation (5) is the differential equation for a damped, driven harmonic oscillator. It is easily solved by means of Laplace transforms. The general solution reads [30]

$$\begin{eqnarray} &&c_{\alpha}(t)= m^{-1}{\textrm{e}}^{-(\gamma+\Gamma_{\alpha})t/2}\int_0^{t}{\textrm{e}}^{{\Gamma_{\alpha}}t_2}\int_{0}^{t_2} \, {\textrm{e}}^{(\gamma-\Gamma_{\alpha}) t_1 /2} \, \delta \mathbf{P}(t_1) \cdot \mathbf{u}_{\alpha} \, \, {\textrm{d}}t_1{\textrm{d}}t_2 \;, \end{eqnarray} \tag{ 6 }$$

with $\Gamma_\alpha = \sqrt{\gamma^2-4\lambda_\alpha/m}$ and $\gamma = d/m$ . Moments µ_p of voltage angle deviations are calculated as averages over the noise distribution. From equation (6), µ_p contains an ensemble-average $\langle \delta P_{i_1} (t_1)\delta P_{i_2} (t_2) \ldots \delta P_{i_p}$ $(t_p)\rangle$, over the product of p sources of noise inside exponential integrals. One therefore needs to specify the moments of the noise distribution. We start from a geographically uncorrelated feed-in noise on nodes labeled i₀, whose first two moments are given by

$$\begin{eqnarray} &&\langle \delta P_{i_0}(t_1) \rangle =0 \, , \end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray} &&\langle \delta P_{i_0}(t_1) \delta P_{i_0}(t_2) \rangle = \sigma^2 \, {\textrm{e}}^{- |t_1-t_2|/\tau_0} \, , \end{eqnarray} \tag{ 7b }$$

to which we add non-Gaussianities in the form of finite skewness and kurtosis of the noise distribution,

$$\begin{eqnarray} &&\langle \delta P_{i_0}(t_1) \delta P_{i_0}(t_2) \delta P_{i_0}(t_3) \rangle = a_3 \, \sigma^3 \, \prod_{m < n} {\textrm{e}}^{-|t_{i_m}-t_{i_n}|/\tau_0} \, , \end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray} &&\langle \delta P_{i_0}(t_1) \delta P_{i_0}(t_2) \delta P_{i_0}(t_3) \delta P_{i_0}(t_4) \rangle_c = a_4 \, \sigma^4 \, \, \prod_{m < n} {\textrm{e}}^{-|t_{i_m}-t_{i_n}|/\tau_0} \, , \\ &&\nonumber \end{eqnarray} \tag{ 8b }$$

where $\langle \ldots \rangle_c$ explicitely refers to a cumulant. This in particular substracts all disconnected averages such as $\langle \delta P_{i_0}(t_1) \delta P_{i_0}(t_2)\rangle\langle \delta P_{i_0}(t_3) \delta P_{i_0}(t_4) \rangle$. The parameters $a_{3,4}$ characterize non-Gaussianities in the noise distribution. They correspond to skewed distributions ($a_3 \ne 0$), with tails longer ($a_4 \gt0$) or shorter ($a_4 \lt0$) than the normal distribution.

The moments µ_p are given by exponential integrals and are straightforwardly calculated. However, their exact expressions are somewhat complicated. We give them for the variance only and, for the third and fourth cumulants, discuss limiting cases of long and short correlation time.

Equation (5) makes it clear that, beside τ₀, the other time scales are the damping time $\gamma^{-1} = m/d$, the $\alpha\textrm{th}$ oscillator period $T_\alpha = \sqrt{m/\lambda_\alpha}$ and the combination $\gamma T_\alpha^2 = d/\lambda_\alpha$ of the two [18]. For the high voltage synchronous grid of continental Europe, a detailed analysis based on realistic line admittances and dynamic parameters gave estimates for these time scales as $\gamma^{-1} \simeq 2.5$ s, $T_\alpha \lt 1$ s and $\gamma T_\alpha^2 \lt 0.4$ s $\forall \alpha$. Therefore the regime of long noise correlation time is already reached for $\tau_0 \gtrsim 5-10$ s, while the short correlation time regime requires $\tau_0 \lesssim 1\ \mu{\textrm{s}}$ [17, 18]. In between lies a hybrid regime, where disturbance time scales overlap with spectral time scales [20].

While circuit breakers and other switches may disconnect power lines and put a power plant off-line in a fraction of a second, disconnection-reconnection sequences may occur at most two to three times consecutively by design. It is hard to think of a significant noise perturbation fluctuating persistently on a time scale shorter than a few seconds. Moreover, in the supplemental material we show several examples of feed-in power fluctuations from renewables, which are all characterized by correlation times $\tau_0 \gtrsim 1$ min in the sense of the two-point correlator of equation (7b ). Hence, we conclude that the long correlation time regime is the relevant one for our investigations.

We first consider the case of a single source of noise and discuss multiple noisy nodes in paragraph 3.6.

In the limit of large observation time, the voltage angle variance is given in equation (S10) of the supplemental material. The two limiting cases of long and short noise correlation time τ₀ give key insights on noise propagation.

First, when τ₀ is the largest time scale, the voltage angle variance at node i reads

$$\begin{eqnarray} &&\lim_{\tau_0 \rightarrow \infty} \langle \delta \tilde\theta_i^2 \rangle = \left( \sigma \sum_{\alpha \geqslant 2}\frac{u_{\alpha,i_0}u_{\alpha,i}}{\lambda_\alpha}\right)^2\, . \end{eqnarray} \tag{ 9 }$$

Because we consider deviations from the geographically averaged response, the sum excludes the zero-mode of the network Laplacian matrix. The quantity squared inside the parenthesis in equation (9) is the Green's function for the linear operator L, from the noise source to the observation node i. For optical or electronic waves propagating through disordered mesoscopic systems, quantities similar to $\langle \delta \tilde\theta_i^2 \rangle$ in equation (9) decay as power laws with the distance between i₀ and i, when averaged over a relatively narrow but high-lying spectral interval [42]. Equation (9) instead corresponds to a 'zero-energy' Green's function, indicating that fluctuations with long correlation times are transmitted by a few low-lying, long-wavelength eigenmodes of L, for which the perturbative approaches of [42] cannot be directly applied.

Second, when τ₀ is the shortest time scale, one obtains

$$\begin{eqnarray} &&\lim_{\tau_0 \rightarrow 0} \langle \delta \tilde\theta_i^2 \rangle = 2 \sigma^2 \tau_0 \sum_{\alpha,\beta \geqslant 2}\frac{u_{\alpha,i_0}u_{\beta,i_0}u_{\alpha,i}u_{\beta,i}}{d(\lambda_\alpha + \lambda_\beta) + \frac{m}{2d}(\lambda_\alpha - \lambda_\beta)^2} \, . \qquad \end{eqnarray} \tag{ 10 }$$

In the inertialess limit, m = 0, the variance is given by a two-particle Green's function in this case.

In the limit of long correlation time, it is straightforward to show that higher cumulants behave similarly to the variance, equation (9), namely (see supplemental material section II, C.)

$$\begin{eqnarray} &&\lim_{\tau_0 \rightarrow \infty} \langle \delta \tilde\theta_i^p \rangle_c = a_p \left( \sigma \sum_{\alpha\geqslant 2}\frac{u_{\alpha,i_0}u_{\alpha,i}}{\lambda_\alpha}\right)^p , \qquad \end{eqnarray} \tag{ 11 }$$

with the parameters a_p giving the deviation from Gaussianity in the feed-in fluctuations, equation (8). The most remarkable thing is that, from equations (9) and (11), standardized higher cumulants are given by $\langle \delta \tilde\theta_i^p \rangle_c / \langle \delta \tilde\theta_i^2 \rangle^{p/2} = a_p$, regardless of the distance between the measured node and the noise source. This is one of the main results of this paper: long correlation times enable non-Gaussian fluctuations from a single noise source to propagate over the whole network and persist at their initial relative value compared to the variance, i.e. to Gaussian fluctuations. This result is in particular independent of inertia, suggesting, as previously found in [43], that disturbances with long characteristic times are affected only marginally by inertia, even when the latter is inhomogeneously distributed. This independence of inertia makes sense since in the long-correlation time limit the system response simply tracks the equilibrium point of equation (2) with perturbed power inputs. Large voltage angle fluctuations are therefore boosted by finite-time correlated disturbances. We stress that, as defined in section 3.1, equation (11) refers to voltage angle deviations from their geographical average. Therefore, the network-wide propagation of non-Gaussianities predicted by equation (11) goes beyond the bulk response reported in [7, 20].

In the limit of short τ₀, the third moment of voltage angles is given in equation (S11) in the supplemental material. In the limit of vanishing inertia, m = 0, it gives,

$$\begin{eqnarray} &&\lim_{\tau_0 \rightarrow 0} \langle \delta \tilde\theta_i^3 \rangle = \sigma^3 \tau_0^2 \sum_{\alpha,\beta,\gamma\geqslant 2}\frac{u_{\alpha,i_0}u_{\beta,i_0}u_{\gamma,i_0}u_{\alpha,j}u_{\beta,j}u_{\gamma,j}}{d^2 (\lambda_\alpha+\lambda_\beta+\lambda_\gamma )} \, , \qquad \end{eqnarray} \tag{ 12 }$$

which, together with equation (10), reflects the fact that, for Kuramoto, i.e. inertialess oscillators, non-Gaussianities in the $p\textrm{th}$ cumulants propagate as a p-particle Green's function in the white-noise limit.

As mentioned previously, Green's functions decay exponentially with distance in disordered mesoscopic systems [42], however it is not clear whether this behavior applies to the "zero-energy" case considered here, nor to p-particle Green's functions. To understand better the propagation of non-Gaussianities in the short correlation time regime, we therefore numerically evaluate the expressions in equations (10) and (12). First, figure 2 shows the theoretically predicted standardized skewness and kurtosis (normalized by a₃ and a₄ respectively) for inertialess networks, as a function of the geodesic distance to the source of noise, for various power grid models. In contrast to the the long correlation time prediction, both skewness and kurtosis decay fast away from the noisy node (seemingly exponentially fast) before they saturate at a constant, small value. This means that, for inertialess networks subjected to fast-varying noise, skewness and kurtosis exhibit an additional decay, on top of the exponential decay of the angle variance reported in [9, 20].

Zoom In Zoom Out Reset image size

It is well known that inertia plays an important role in absorbing fast fluctuating disturbances with short correlation time. Therefore, we further show in figure 3 the standardized kurtosis as predicted by the formulae of equations (S10) and (S11), for various values of inertia. As inertia increases, the kurtosis is further suppressed compared to the variance. The influence is especially strong a short distance away from the noise source. This agrees with numerical data reported in [9] and corroborates earlier findings that inertia impacts the dynamics mainly locally and at short times [43].

Zoom In Zoom Out Reset image size

The data presented in figures 2 and 3 show that short noise correlation times suppress the propagation of non-Gaussian voltage angle fluctuations through meshed networks over a network-dependent distance. This is an additional effect, further suppressing non-Gaussianities relative to the already decaying voltage angle variance. The effect becomes stronger at larger inertia. Note also that the grid models used in figures 2 and 3 fairly capture the parameters of actual high-voltage power grids.

We finally consider the case with M distinct, independently but identically distributed sources of power feed-in fluctuations. In that case, there are M contributions similar to that in equation (11), but $M! (M-p/2)!/(p/2)!$ pairings of the noise sources for the moment of even order p. These latter contributions are much more numerous and they result in a Gaussian $p\textrm{th}$ moment—this is the standard mechanism behind central limit and Berry–Esseen theorems [45]. For instance for the p = 4 moment in the large correlation time limit, one obtains

$$\begin{eqnarray} &&\lim_{\tau_0 \rightarrow \infty} \langle \delta \tilde\theta_i^4 \rangle =\sum_{i_0=1}^M \Big(\sigma \sum_{\alpha\geqslant 2}\frac{u_{\alpha,i_0}u_{\alpha,i}}{\lambda_\alpha}\Big)^4 + 3\sum_{i_0 \neq j_0} \Big(\sigma \sum_{\alpha \geqslant 2}\frac{u_{\alpha,i_0}u_{\alpha,i}}{\lambda_\alpha}\Big)^2 \Big(\sigma \sum_{\beta \geqslant 2}\frac{u_{\beta,j_0}u_{\beta,i}}{\lambda_\beta}\Big)^2 , \qquad \end{eqnarray} \tag{ 13 }$$

where the factor 3 in the second line accounts for all possible pairings between the product of four noise sources, $\delta P_{i_l}$, $l = 1, \ldots 4$. The pairing mechanism giving the second term on the right-hand side of equation (13) leads to the convergence of the voltage angle distribution to a Gaussian distribution, with $\langle \delta \tilde\theta_i^4 \rangle/\langle \delta \tilde\theta_i^2 \rangle^2 \rightarrow 3$ [to see this, sum over the noisy nodes $i_0 = 1, \ldots M$ in equation (11) and compare the result with (13)]. The convergence is the same as in the Berry–Esseen theorem [45]. When M is large, this second term dominates over the first one by a factor $(M-1)/2$, so that the ratio of the fourth cumulant—a measure of non-Gaussianity—to the fourth moment becomes $\propto M^{-1}$. With the standard definition [45], non-Gaussianities disappear at a rate $\propto M^{-1/2}$.