Physics-guided Residual Learning for
Probabilistic Power Flow Analysis

PF analysis aims to analyze the steady-state operating points of an electrical grid. The operational data includes power generation, load demands, voltage phasors, and branch flows. There are three types of buses: PQ buses, PV buses, and one slack bus. A PQ bus (a.k.a. load bus) has no generator attached, where its active and reactive power injections are fixed. A PV bus (generator bus) has generators connected, and its active power injection and voltage magnitude are known. Lastly, the slack bus has the given voltage angle and voltage magnitude. Therefore, the unknown variables include the voltage angles and voltage magnitudes of PQ buses and the voltage angles of PV buses.

Consider a power grid with $N$ buses, where there are one slack bus, $N_{g}$ PV buses (denoted by set $\mathcal{N}_{g}$), and $(N-N_{g}-1)$ PQ buses (denoted by set $\mathcal{N}_{l}$). The number of unknown variables is $2\times(N-N_{g}-1)+N_{g}$, which is the same as the number of power balance equations.

Let $P_{i}$ and $Q_{i}$ represent the active and reactive power injections at bus $i$ while $V_{i}$ and $\theta_{i}$ are the voltage magnitude and angle. The AC-PF equations are the forward mappings from voltage phasors to power injections, which are given as:

where $\theta_{ij}:=\theta_{i}-\theta_{j}$ is the voltage angle difference between bus $i$ and bus $j$. $G_{ij}$ and $B_{ij}$ are the real and imaginary parts of the $(i,j)$-th element of the nodal admittance matrix $\mathbf{Y}\in\mathbb{C}^{N\times N}$. Moreover, the active and reactive power branch flows between two connected buses $i$ and $j$ are:

where $b_{ij}^{c}$ is the total line-charging susceptance.

Clearly, the power injections and branch flows are determined by all the voltage phasors, which are defined as the state of the system:

where subscripts $(\cdot)_{s}$, $(\cdot)_{g}$ and $(\cdot)_{l}$ denote the quantities corresponding to the slack bus, generator buses, and load buses, respectively. Let $\mathbf{P}_{g}$, $\mathbf{P}_{l}$, and $\mathbf{Q}_{l}$ denote the active power of generator buses and active/reactive power injections of load buses. Partitioning and reorganizing the admittance matrix $\mathbf{Y}$ in the same manner yielding

where $\mathbf{Y}_{gs}$ is formed from the rows of $\mathbf{Y}$ that correspond to generator buses and the column for the slack bus. Similarly, we can find all other blocks. Define

Hence, the inverse mappings of the AC-PF equations can be compactly expressed as:

The PPF analysis is closely related to the aforementioned AC-PF problem. Consider an electricity grid with stochastic renewable power generation and load demands. PPF studies aim to characterize the uncertainties of voltage phasors induced by the fluctuations of power injections. As a new effort in PPF studies, NNs are employed to learn the end-to-end mapping (7) from historical data pairs $(\mathbf{x},\mathbf{y})$. NNs take in power injections $\mathbf{x}$ and predict the corresponding voltage phasors $\mathbf{y}$. The PDFs of voltage phasors can be further inferred from the output samples. By shifting the time-consuming training process offline, NNs can rapidly predict the corresponding voltage phasors of new power injection samples in the testing stage. Even with lots of new input samples, NNs are still computational efficiently because the forward propagations during the test stage are often very fast.

Unlike the building blocks in the MLP, residual blocks introduce identity mapping as a shortcut connection that skips one or more layers [27]. Fig. 1a illustrates the structure of one residual building block. In Fig. 1a, let $\mathbf{G}(\cdot)$ denote the mapping from $\mathbf{u}$ to $\mathbf{\bar{v}}$. The motivation for formulating the residual building block is as follows. MLP consists of a few stacked layers and shows universal function approximation capabilities. Thus, it is reasonable to assume that the MLP can approximate the residual function $\mathbf{R}(\cdot)$. Then, if those stacked layers are only used to approximate $\mathbf{R}(\cdot)$, the original function $\mathbf{G}(\mathbf{u}):=\mathbf{R}(\mathbf{u})+\mathbf{u}$ can be obtained by introducing an extra identity shortcut connection. It has been shown that ResNets are also universal function approximators [28].

However, the benefits of ResNet over MLP in PF analysis are not obvious. The reason is as follows. The NNs applied in the image recognition domain are usually composed of dozens or hundreds of layers. Therefore, deep ResNets will have apparent benefits over deep MLPs because they do not have degradation issues. However, the advantages are not evident because very deep NNs may not be necessary for approximating the inverse AC-PF equations.

In addition, [27] claimed that if the desired function is close to the identity function, it will be easier for stacked layers to learn the perturbations regarding the identity instead of approximating that desired function directly. Inspired by their work, our proposition is that compared with MLPs that directly approximate the complicated inverse AC-PF equations, pushing the residual correction regarding the linear relationship to zero may be more accessible. This motivates us to design a novel architecture with some initialization methods tailored for PF analysis.

In practical power systems, the voltage magnitudes of buses are approximately 1 per unit, and voltage angle differences are small values [29]. Many existing works use these two properties to linearize the AC-PF equations, and numerical results have shown that these linear models perform well. Therefore, we make novel modifications to the MLP structure to better use this linear property. As shown in Fig. 1b, we replace the identity mapping with a fully connected linear layer. Compared with the MLP that tries to approximate the mapping from $\mathbf{x}$ to $\mathbf{y}$ directly, our designed framework uses a set of stacked layers to only approximate the latent residual functions $\mathbf{R}(\cdot)$. The stacked layers are composed of fully connected linear layers and rectified linear unit (ReLU) activation function [30]. Let $\mathbf{W}_{s}$ and $\mathbf{b}_{s}$ denote the weight matrix and bias of the shortcut connection linear layer, respectively. $\mathbf{W}_{i}$ and $\mathbf{b}_{i}$ are the weight matrix and bias of the $i$-th fully connected linear layer. $\bm{\sigma}$ is the activation function. The residual output $\mathbf{y}_{r}$ and the final output $\mathbf{y}$ can be written as:

In NN training, learnable weights are randomly initialized which can be critical to the NN performance. We develop two physics-guided initialization methods for the weights of the shortcut connection layer. These two methods significantly accelerate the training process and improve the NN performance than random initialization. The goal is to get $\mathbf{y}_{r}$ close to zero in the initial stage of the training process by initializing $\mathbf{W}_{s}$ and $\mathbf{b}_{s}$ properly. If $\mathbf{W}_{s}$ and $\mathbf{b}_{s}$ are randomly initialized, the value of $\mathbf{y}_{r}=\mathbf{y}-(\mathbf{W}_{s}\mathbf{x}+\mathbf{b}_{s})$ cannot be zero. Thus, to drive $\mathbf{y}_{r}$ to zero, $\mathbf{W}_{s}\mathbf{x}+\mathbf{b}_{s}$ should be close to the output $\mathbf{y}$. In other words, the shortcut connection linear layer should serve as an excellent linear approximation from $\mathbf{x}$ to $\mathbf{y}$ after initializing its weights. Therefore, we modify two linear PF models and apply them to initialize $\mathbf{W}_{s}$ and $\mathbf{b}_{s}$ in the proposed framework.

The line parameter profiles are only available from the grid planning files, which are likely outdated [22]. If accurate information on line parameters is unavailable, the aforementioned physics-guided initialization methods are not applicable. Therefore, we propose a data-driven initialization scheme by leveraging ridge regression to deal with this situation. The motivation for using ridge regression is as follows. First, [32] shows the load demands and power generation data have similar rise and fall patterns in adjacent areas. The high correlations among the input variables lead to multicollinearity, weakening the performance of the simple linear regression model. Ridge regression is an effective method for analyzing data that suffer from multicollinearity [33]. It introduces the $\ell_{2}$ regularization, which shrinks the regression coefficients to reduce model variance and avoid overfitting. Secondly, the initial weight parameters of NNs are typically random numbers following certain distributions. Small random values are generally considered better options than larger ones in preventing gradient from exploding or vanishing [34]. Due to the penalty term, ridge regression tends to give smaller random weights than simple linear regression models. Therefore, we employ the coefficients obtained from ridge regression to initialize the shortcut connection layer, which provides better initial values in stabilizing and speeding up the training process.

With $n$ training samples, let $\mathbf{X}\in\mathbb{R}^{n\times(2N-N_{g}-2)}$ and $\mathbf{y}^{o}_{i}\in\mathbb{R}^{n\times 1}$ denote the input data and the $i$-th dimension of the output data, respectively. Ridge regression finds the coefficient vector $\mathbf{w}_{i}\in\mathbb{R}^{(2N-N_{g}-2)\times 1}$ and bias $b_{i}\in\mathbb{R}$ by solving the problem:

where $\mathbf{1}_{n}\in\mathbb{R}^{n}$ is the all-ones column vector. $\mathbf{w}_{i}^{\top}$ and $b_{i}$ will be used to initialized the $i$-th row of $\mathbf{W}_{s}$ and $\mathbf{b}_{s}$. The first term in the loss function (15) calculates fitting errors of data pairs, while the second term shrinks the estimated coefficients. $\lambda$ is a tuning parameter that balances the relative strength of the least square error and the penalty term. When $\lambda=0$, it degrades to the least-square based linear regression.

The closed-form solution to (15) is derived as follows. We can rewrite the objective function $L$ as:

where $\mathbf{I}\in\mathbb{R}^{(2N-N_{g}-2)\times(2N-N_{g}-2)}$ is the identity matrix. The gradients of the convex objective function are

Set $\nabla_{\mathbf{w}_{i}}L=\mathbf{0}$ and $\nabla_{b_{i}}L=0$, we get:

After plugging (20) into (19), we eliminate $b_{i}$ and obtain the closed-form solution of optimal $\mathbf{w}_{i}$ as:

with $\mathbf{H}=\frac{1}{n}\mathbf{1}_{n}\mathbf{1}_{n}^{\top}\in\mathbb{R}^{n\times n}$.

Table II shows our proposed approaches outperform the other methods in predicting the voltage magnitudes and angles for the PF analysis. Besides, we have the following observations:

• •

  Classical non-linear regressors in the machine learning field, including KNN, RR, and SVR, have significant estimation errors in medium- and large-scale bus systems. These solvers cannot effectively extract the data features and recover the complicated mapping relationship in the inverse AC-PF equations.

• •

  The NN-based approaches achieve more accurate predictions than the linear LPF method, the non-linear SML method, and classical non-linear regressors, especially in the voltage angle estimates.

• •

  The performance of the FC method is comparable to that of the ResNet and Random approaches. It shows that only using the residual framework does not have apparent advantages over the vanilla MLP structure.

• •

  Three designed initialization schemes seem straightforward but significantly outperform the FC and random initialization methods in improving the accuracy of voltage phasor estimates.

In addition, Fig. 2 and Fig. 3 show the MAPEs of voltage angles and magnitudes on the IEEE-118 bus system. The average MAPEs of voltage angles and magnitudes of the data-driven method are 0.028% and 0.0042%. Similarly, 0.023% and 0.0065% for the Linearized PF method, and 0.023% and 0.0042% for the Jacobian method.

Table III shows the errors of branch flow calculations. Data-driven, Linearized PF and Jacobian schemes achieve the best performance among all the NN-based solvers. However, on the IEEE-300 and PEGASE-1354 bus systems, the LPF and RR methods obtain accurate branch flow calculations despite their poor performance in voltage phasor estimates. The reason is as follows. The relationship between voltage angles and voltage angle differences is not bijective. Branch flows between two buses are more related to voltage angle differences than voltage angles. However, the outputs of our proposed models are voltage angles instead of voltage angle differences. Therefore, our methods may predict $\theta_{i}$ and $\theta_{j}$ accurately, while there is no guarantee of a similar level of accuracy for $\theta_{ij}$ estimate. This property further affects the calculation accuracy of branch flows [21].

The performance of the proposed approaches is evaluated based on risk assessment, which is critical to power system operation. Tables VII and VIII show the probabilities of exceeding the operational limit of the voltage magnitude and apparent branch flow on the SouthCarolina-500 bus system [48]. Compared with the MCS method, our proposed approaches have achieved promising results in capturing the violation possibility values in the order of $10^{-3}$. Furthermore, by employing the violation metrics for voltage magnitude and apparent branch flow, Tables IX and X present the 95% confidence intervals derived from the MCS, alongside the average violation values yielded by our innovative methods. Notably, we observe that the mean violation degrees estimated through our proposed approach consistently fall within the 95% confidence intervals established by the MCS technique.

In addition, one of the stopping criteria of the MCS method is the variance coefficient for the violation possibility to be smaller than 1% [49]. Within the simulations, accurate estimates of voltage magnitude violations and apparent branch flow violations necessitate a minimum of 2700 and 5100 samples, respectively. Thus, the cumulative computational time required to perform Monte Carlo simulations for power flow analysis across thousands of instances can be substantial. By contrast, our proposed methods can rapidly predict voltage phasors for thousands of data samples.

The simulations are implemented on an iMac with i7-8007 CPU and 32GB RAM and a Linux server with NVIDIA Tesla K20-5 GB GPU. The NN training is implemented via PyTorch 1.7.1 in Python 3.7. Table XI shows the NN training time of each epoch. Residual building blocks need a slightly longer time due to the extra propagation of the shortcut connection layer. Table XII shows our proposed approaches significantly reduce the total computational time compared with the MCS and QMC (with 3000 sampling points) algorithms.

The learning rate affects the loss convergence rate of the NN training process. A significant learning rate helps fast convergence but may lead to the NN weights converging to a suboptimal solution. Hence, we adopt a relatively small learning rate of $10^{-4}$. Note that this small learning rate is only used in this part to indicate the convergence properties. The MSE evolution of the training loss in the first 500 epochs is shown in Fig. 8 and Fig. 9. After training the first epoch, the MSEs of our proposed approaches are more than two orders less than that of others. Besides, even after training 500 epochs, other NN-based methods cannot achieve the same loss level as ours. Therefore, our proposed approaches have significant advantages over others regarding the convergence rate, which can be attractive in the face of a limited training time.

In addition, three designed initialization schemes converge faster than the random initialization. Therefore, we show how designed initialization methods influence the NN weights update during the training process. Fig. 10a shows the randomly distributed parameters of the shortcut connection layer. After training 500 epochs, the pattern of parameters is still very random, as shown in Fig. 10b. In contrast, Fig. 11 shows that the pattern of the updated parameters after training is still quite similar to that of the initial parameters. This phenomenon indicates that these learnable parameters are finely updated based on the initial weights during the training process. Therefore, we conclude that our well-designed initial weights play a critical role in NN training.