Physics-Informed Graphical Neural Network
for Parameter & State Estimations in Power Systems

The term “Physics Informed Machine Learning” (PIML) in reference to ML was coined in the series of Los Alamos “Physics Informed Machine Learning” PIML-Workshops which took place in 2016, 2018 and 2020. [1]. The term came in reference to two complementary aspects – (a) Physics providing an intuitive and often constructive input into design, analysis and implementation of various ML schemes, and (b) when ML helps to solve a physical problem. This manuscript is concerned with the later aspect of PIML. The approach often translates into embedding equations, describing related physics for the application of interest, into a ML scheme. These ideas originate from early work, e.g. of [27], on tuning a NN to satisfy the output of an equation (algebraic or differential). Noticeable applications of the methodology came in the context of dynamical equations (ODEs and PDEs), including the so-called Sparse Identification of Nonlinear Dynamics (SINDy) [5], Physics Informed Neural Network (PINN) [44] which was also applied to the PS dynamics in [34], and Neural ODE [8] frameworks. See also most recent review [53], and references there in, e.g. emphasizing different objectives in combining elements PIML with the state-of-the-art ML models to leverage their complementary strengths. Physics informed NN modeling in applications to learning SE under limited observability was also discussed in [35].

Graph (Convolutional) Neural Network (GNN) is a NN which we build making relations between variables in the (hidden) layers based on the known graph associated with an application, e.g. PS. In this regards, GNN is informed, at least in part, about physical laws and controls associated with the power system operations. One of the first GNN references [25] focuses on citation networks, and then the approach spread into many applications [58] including geo-location [42], infrastructure modeling [11], system of interacting particles [45] and PSs (mentioned in the next Subsection).

PS-agnostic techniques have focused on making SE based on historical and (limited) measurement data utilizing various modern tools of ML, such as Feed-Forward Neural Networks (FFNNs) [48, 33, 55, 56], autoencoders [3], and most recently Convolutional Neural Network (CNN) [28] and Graph Neural Network (GNN) [4, 17, 36, 30, 24, 29, 7] (some of these in the context of a related problem of fault detection). It is also worth mentioning most recent NeuroIPS 2020 Learning To Run Power System (L2RPS) competition [2], won by teams utilizing PS-agnostic approaches. Physics-Informed approaches to SE, and related problem of PE in static and dynamic models of PS, actual or reduced, state the problem as a regression based on the PF equations [32, 23, 10, 59, 22, 60, 9, 6, 52, 31, 35], and also extend in the related areas of probing the proximity to instability and helping in design of the corresponding emergency control actions [20, 51], optimization and resource allocation [41, 13]. Some of these approaches also dealt with the partial observability [12, 14, 38, 15, 31], which is especially well pronounced on the lower voltage (distribution) side of the PS. The hybrid approach, attempting to mix the best of the two aforementioned approaches (and inspired by universal methodologies [27, 5, 43] suggesting to blend equations with modern ML) was also explored in [34] to estimate the PS dynamics. Ability to predict SE and learn physical models, i.e. to make PE, based primarily on the measurements represents an attractive feature of the framework we advance in this manuscript. In addition to applications mentioned above, validated SE and PE are significant for many other tasks of importance for the system instability, protection and energy management reviewed in many classical books on the subject, e.g. [26, 46].

Importance of developing reduced order static and dynamic models, providing computationally light equivalent representation of actual PS, was long time recognized as important, see e.g. review [19] and extensive historical references therein. We will highlight in this brief Subsection only a small subset of the waste literature on the subject directly related to the so-called Kron reduction methodology [16, 18] most relevant to the manuscript.

Power system should be balanced on the time scale of seconds. This means that electricity generation ought to match the consumption (i.e. the combined value of electric loads and electric/heating losses) at any time. Electricity demand varies over a range of scales – hourly, daily, weekly or seasonal variations are of a concern in different settings. For example, in an area with a continental climate, winter demand is generally higher than summer demand. This variations may be linked to local specifics, e.g. dependence on the electric heating, as in France. In warmer areas, such in the Southern regions of US, the summer demand exceed the winter demand, it is mainly due to air-conditioning. Shorter time scale fluctuations of loads, observed on the scale from seconds to minutes, are mainly due to activities of many small appliances. Note also that renewables, such as wind and solar, even though they contribute generation, puts an additional pressure on the system injecting uncertainty due to wind and solar which needs to be balanced by other generators or batteries. Operational redistribution of generation between flexible generators is implemented via energy market mechanisms. Basic optimization behind this task, the so-called Optimal Power Flow (OPF), is resolved frequently, in some energy markets as often as once every 3-5 minutes. OPF is a constrained optimization, where in addition to generation constraints (limits on the generation possible at any generator) and thermal constraints (limiting power or current flowing through each line of the power system) we also add the so-called Power Flow (PF) equations, expressing Ohm’s laws, and relating injections/consumptions of power at all the nodes of the system to phases, voltages at the nodes and power flows over the power lines. PF equations are defined and discussed in more details in the following Subsection.

Consider a Power System (PS) which operates over a grid-graph, ${\cal G}=({\cal V},{\cal E})$, where ${\cal V}$ and ${\cal E}$ are the set on nodes (generators or loads) and set of lines, respectively. PF equations, governing steady redistribution of power over the system, are stated in terms of the complex-valued (also called complex) powers, $\forall i\in{\cal V}:\ {S}_{i}\equiv p_{i}+{\bm{i}}q_{i}$, and in terms of the complex-valued (electric) potentials, $\forall i\in{\cal V}:\ {V}_{i}\equiv v_{i}\exp({\bm{i}}\theta_{i})$, where $v_{i}$ and $\theta_{i}$ denote voltage (magnitude) and phase of the potential at the node $i$, and ${\bm{i}}^{2}=-1$. In these notations, the PF equations become, see e.g. [26, 46], $\forall i\in{\cal V}$ ¹¹1We simplify a bit and did not include in the PF equations shunt capacitors, associated with nodes of the systems and representing effects of self-inductance and self-conductance of the PS nodes/buses. :

where $b_{ij}$ and $g_{ij}$ are susceptance and conductance of the power lines, $\{i,j\}\in{\cal E}$; PF Eqs. (1,2) can be pictured as defining implicitly the PF map, ${\bm{\Pi}}:{\bm{S}}\equiv({S}_{i}|i\in{\cal V})\mapsto{\bm{V}}\equiv({V}_{i}|i\in{\cal V})$ (The map is implicit, because it requires solving the PF equations, which may have no or may have many solutions.) It can also be stated as an (explicit) inverse map, ${\bm{\Pi}}^{-1}:{\bm{V}}\mapsto{\bm{S}}$. To emphasize parametric dependence of the PF map on the graph, ${\cal G}$, and also on the matrix of the power line admittances defined over the graph, ${\bm{Y}}\equiv{\bm{g}}+{\bm{i}}{\bm{b}}=(g_{ij}+{\bm{i}}b_{ij}|\{i,j\}\in{\cal E})$, we use the following notations, ${\bm{S}}={\bm{\Pi}}^{-1}_{\bm{Y}}({\bm{V}})$.

SE – main task in the center of this manuscript – consists in, given available observations, to predict other relevant physical characteristics of the system. For example, assuming that the PS network, i.e. ${\cal G}$ and ${\bm{Y}}$, are known, and that voltages and phases are measured by PMUs at all nodes of the system, the SE task may be to estimate injected/consumed active and reactive powers. This task can obviously be resolved by explicit (that is simple) application of the inverse PF map, ${\bm{\Pi}}^{-1}$, set by the PF equations.

In the case of a limited observability, i.e. when some nodes are unobserved the task of SE can be viewed in two complementary ways. We can pose the question of complementing missing power injections/consumptions only at the nodes where voltages and phases are actually measured. Alternatively, we can also pose a more challenging version of the SE under limited observability — in addition to completing measurements at the observed nodes, also to reconstruct injected/consumed powers and also voltages and phases at all nodes of the system. In this manuscript we will focus on the former, less complicated, but still challenging enough version of the SE.

The problem of PE consists in reconstructing from PMU measurements structure of the PS graph and also line characteristics, i.e. ${\bm{Y}}$. This question is relevant because parameters of the power line are known to change due to change in the operational conditions, added or removed vegetation and aging. These changes are generally infrequent but still significant, e.g. for making accurate state estimations under complete or partial observability.

Motivated by the case of partial observability, we are interested to pose the question of finding equivalent models of power systems providing a description similar to the one provided by the PF Eqs. (1,2), however stated over a reduced set of nodes. Given that the PF Eqs. (1,2) are nonlinear the quest of finding an equivalent (reduced) model is nontrivial. One, admittedly phenomenological (empirical) approach may consists in drawing some useful conclusion from the current-voltage version of the problem (and not power-voltage used in the PF setting), where model reduction can be evaluated explicitly.

Specifically, let us discuss the so-called Kron reduction – see [18] and references therein. We divide the nodes of the system into two not-overlapping subsets of the observed (labeled $o$) and unobserved (labeled $u$) nodes. Then injected/consumed currents, ${I}_{i}={P}_{i}/{V}_{i}$ (not powers!) and potentials of the PS are related to each other according to the Ohm’s law describing linear relations between currents and potentials (see Supplementary File for details). Linearity of the relation between the injected/consumed currents and potentials allows to get a closed relation between currents and voltages associated solely with the observed nodes:

It is therefore suggestive to use the graph structure, ${\cal G}^{(r)}\equiv({\cal V}^{(o)},{\cal E}^{(r)})$, associated with the reduced admittance matrix ${\bm{Y}}^{(r)}\doteq(\{i,j\}|Y_{ij}^{(r)}\neq 0)$, for building an equivalent PF model – akin Eqs. (1,2), $\forall i\in{\cal V}^{(o)}$, i.e.

where, ${\bm{Y}}^{(r)}$ denote the reduced (we can also call it “equivalent”) admittance matrix associated with the effective (not necessarily real) power lines, $\{i,j\}\in{\cal E}^{(r)}$. ${\bm{Y}}^{(r)}$ is to be guessed, or as we do it in the following, to be learned.

Eq. (4) sets up the backbone of the PIML scheme we develop for the case of partial observability. Specifically, we aim to solve the following parallel SE and PE learning problem, which we call Power-GNN:

where the $l_{2}$-norm, $\|\cdots\|^{2}$, assumes summation over observed nodes of the PS; $|{\cal V}^{(o)}|$ is the number of the observed nodes; ${\bm{S}}^{(o)}_{n}$ and ${\bm{V}}^{(o)}_{n}$, with $n=1,\cdots,N$, denote $N$ samples of the vectors of the complex power and of the electric potential, respectively, measured at the observed nodes, $i\in{\cal V}^{(o)}$. Here in Eq. (6), $\Sigma_{\bm{\varphi}}\big{(}{\cal V}_{n}^{(o)},S_{n}^{(o)}\big{)}$ is introduced to represent effects of the part of the system which is not observed; ${\bm{\varphi}}$ stands for the vector of parameters representing the hidden part, and $\mathcal{R}({\bm{\varphi}})=\alpha\|{\bm{\varphi}}\|^{2}$ is the regularization term, chosen $l_{2}$.

In the following, and specifically in Section IV.2, we will introduce the Power-GNN scheme implementing optimization (5). The scheme will be training simultaneously physics-loaded parameters (elements of the admittance matrix, ${\bm{Y}}^{(r)}$), and physics-blind parameters representing the $\Sigma_{\bm{\varphi}}\big{(}{\cal V}_{n}^{(o)},S_{n}^{(o)}\big{)}$ function, via a Graphical Neural Networks (GNN) defined over ${\cal G}^{(r)}$, constructed according to the Kron reduction procedure introduced above in Section III.4.

To mimic realistic market operations, the PF physics, as well as variability and uncertainty due to fluctuations of loads and renewables, we build the following procedure for generating data samples, which are used in the following to train the ML models we build for the SE.

We work with three exemplary networks, shown in Fig. (1), representing small, medium and large PSs, and generate synthetic data (to train respective ML models) utilizing a popular OPF solver – the MATPOWER  [61]. For each of the networks we generate five data sets that we label as cases #1-5. (Later on we will also explain an additional case #6.) Our data generation procedure for case #1 is as follows:
$\bullet$ Pick a reference configuration of loads and find respective solution of the OPF problem with the generation cost fixed.
$\bullet$ Devise a family of load configurations each with the same proportion of the loads at different nodes (the same participation factor), but rescaled differently. The rescaling factor is chosen to guarantee the system reside within the pre-defined interval between the reference (base) load and the peak load. For each load configuration from the family solve the OPF (with the same cost of generation per generator as in the reference case).
$\bullet$ We arrive at the family of samples each represented by complex potentials and injected/consumed powers known for all nodes of the system.

To produce case #2 data set we modify procedure injecting random load fluctuation at every node (distributed according to i.i.d. Gaussian distribution with variance equal to a constant fraction of the load at the respective node). For the case #3 data set active and reactive load are allowed to vary independently. For case #4, one third of generators is removed from the pool of “ready to generate” generators (This relates to the unit commitment problem). A new set of excluded generators is randomly selected for each sample. Finally, for case #5 data set, the generation cost of generators is randomly chosen in a predefined range, new costs are picked for each sample. Each of the data sets for each of the models consists of 2000 samples.

We conclude this Subsection with a methodological remark, emphasizing specifics of the datasets representing realistic PS. Actual states of a PS are distributed over respective space in a highly structured way, which is certainly rather far from uniform. This is due to special correlations induced by the OPF, PF and other procedures, inheriting structure and correlations of the underlying physics and controls. This highly nonuniform representation of the phase space makes the PS learning problem more challenging than in other applications, as it requires to go beyond the standard interpolation — here we ought to extrapolate into regimes which were either unseen or rarely visited.

In this brief Subsection we describe two NN schemes – our custom Power-GNN – combining elements of the PF equations and GNN, and – a standard (one may call it “vanilla”) NN introduced for comparison (as a bench-mark for the Power-GNN).

Power-GNN: is the algorithm which solves Eq. (5) by minimizing the Loss Function (6) via the gradient descent (we use the Adam algorithm) over the physical parameters, representing the matrix of admittances, ${\bm{Y}}^{(r)}$ ²²2We also include in this matrix self-admittances of the nodes associated with shunt-capacitors., and also parameters of the Neural Network representing effect of unobserved degrees of freedom on the observed physical characteristics and expressed via $\Sigma_{\bm{\varphi}}\big{(}{\cal V}_{n}^{(o)},S_{n}^{(o)}\big{)}$ term in Eq. (5). The total number of physical parameters which we learn in the process of training of the PIML model is twice the number of nodes plus twice the number of edges in ${\cal E}^{(r)}$. Even though we call parameters of the NN – ${\bm{\varphi}}$ – physics blind, we are still attempting to inject in them some physical meaning by building the NN based on the graph, ${\cal G}^{(r)}$, associated with the Kron reduction (see discussion of Section III.4). Specifically, we build Graphical Neural Network (GNN) where graph-convolutions, guiding relation between neurons in different layers, follow ${\cal G}^{(r)}$. More info on the construction of the GNN is given below.

Vanilla NN: is presented as a benchmark. We choose a NN with four fully-connected (hidden) layers each containing $2\cdot N_{\rm bus}$ neurons. We choose this simple architecture after checking, empirically, that it is the smallest layered fully-connected NN which is trained successfully on all the SE data-sets described in Section IV.1. We tested both RELU and SoftSign activation functions – and observe similar results. We use the following Loss Function

to train the Vanilla NN, mapping samples of the electric potentials to the samples of complex powers at the observed nodes. We call training successful when $L_{\text{NN}}$ reaches $10^{-4}$ on the training set. (Learning rate and number of epochs is chosen to reach the goal. (See supplementary file for additional details on our choice of the Vanilla NN, including selection of the hyper-parameters.)

We train both Power-GNN and Vanilla-NN using only fifth of the respective dataset, and then validate it (also to verify that we do not overfit) on the complete dataset.