structures can be efficiently represented as hypernode-based nested graphs
[145], which are composed of hypernodes, which are themselves graphs. Learning
over these types of graphs is done using multiple levels of GNNs, which
operate at different hierarchy levels and exchange information at particular
message passing iterations.

Because data in real-world problems is typically sparsely labelled due to high
labelling costs, techniques that make the best use of both unlabelled and
labelled data samples are required. Semi-supervised learning techniques fall
somewhere between supervised and unsupervised learning, with the goal of
making predictions on unseen data based on small sets of labelled data while
utilizing large amounts of unlabelled data. The self-training method [146],
also known as self-labelling, is a simple form of semi-supervised learning for
classification problems in which an initial classifier is trained on a small
labelled dataset and then used to label the rest of the data using the most
confident predictions. Then, using the larger labelled dataset, a new
classifier is trained, and the process is repeated several times. More
advanced approaches employ generative models to learn the distribution of
real-world data from unlabelled samples [147]. Generative models are then used
to enrich the labelled dataset by generating new data samples with the desired
label. GNNs can adapt well to the semi-supervised concept when dealing with
node-level prediction problems in graph-structured data, where labels are only
available to a small subset of nodes [133]. In this scenario, the GNN model is
trained by backpropagating the supervised learning loss function calculated
using only the labelled nodes. However, because all nodes and their incident
edges participate in the message passing process, the node input features and
connectivity data from unlabelled nodes also contribute to prediction
generation.

Self-supervised learning [148] is another type of methods that falls somewhere
between supervised and unsupervised learning. It uses unlabelled data samples
to create a supervised learning objective using pseudo-labels, resulting in
useful data representations (i.e., embeddings) that can be used as inputs for
other, so-called downstream tasks. Similarly to semi-supervised learning, it
became popular because of its ability to avoid annotating large amounts of
data. Because no manual data labelling is required, some classifications
consider self-supervised learning to be a subset of unsupervised learning. One
of the most common forms of self-supervised learning is contrastive learning
[149], in which data samples are associated with positive and negative pseudo-
labels automatically. A subset of data samples, for example, can be
intentionally corrupted by adding some noise and labelled as negative, whereas
uncorrupted data samples are labelled as positive. Models defined in this
manner are trained using the binary classification loss (e.g., cross-entropy)
as a learning objective and produce data representations that can be reused
for downstream tasks as an intermediate result.

To make use of unlabelled graph-structured data, GNNs can also be used in
self-supervised setting. Link prediction can be thought of as a simple form of
self-supervised contrastive learning based on GNNs. It is essentially a binary
classification task in which pairs of nodes are sampled, and classification
labels are generated based on the presence of direct connections between the
nodes. This task can be generalized so that positive pseudo-labels indicate
that a pair of nodes is close in terms of some neighbourhood metric (e.g.,
$k$-hop or random walk-based), and negative pseudo-labels indicate that they
are not. Node embeddings obtained using these positive and negative label
definitions can be useful for tasks which rely on local graph information. To
obtain more relevant node embeddings for global tasks, contrastive learning
methods that distinguish between original and corrupted graphs can be used
[138].

We conclude the overview of the practical aspects of GNNs by discussing the
augmentation methods that are commonly used with graph-structured data. In
real-world problems, the original, raw graph data is rarely used as an input
for GNNs without being augmented in some way. Since nodes in a graph can lack
input features, feature augmentation in the form of expanding the node input
feature vectors is often used. Let’s consider a case of a graph in which nodes
do not have any input features. A trivial form of feature augmentation is
adding a constant feature (e.g., a scalar with a value of $1.0$) to every
node. As a result, the GNN is able to learn the structural information of a
graph using the aggregation function, which is not possible when aggregating
zeroes. More advanced types of feature augmentation, such as one-hot vector
encoding, cycle counts, node centralities, and clustering coefficients, can
reveal more information about the graph structure and increase the expressive
power of GNN. Aside from feature augmentation, graph augmentation is used when
the input graph is either too sparse, making the message passing process
inefficient, or has long-range dependencies between nodes. Some common graph
augmentation methods include connecting pairs of nodes using virtual edges and
adding virtual nodes that connect all the nodes in a certain subset.

###  Chapter 4 Graph Neural Network-based State Estimation

In this chapter, we explain how GNNs can be applied to both linear and
nonlinear SE, using the power system factor graph-like structures. Since we
used similar methodologies for both linear and nonlinear formulations, they
are presented simultaneously, with specific differences highlighted as
necessary. First, we present the augmentation techniques for the power
system’s factor graph, then the details of the proposed GNN architecture, and
analyse the computational complexity and the distributed implementation of the
GNN model’s inference. Finally, we demonstrate the effectiveness of our
proposed method through numerical evaluations on various test cases.

####  4.1 Power System Factor Graph Augmentation

Inspired by recent work on using probabilistic graphical models for power
system SE [150], we first create a GNN over a graph with a factor graph
topology. This bipartite graph consists of factor and variable nodes connected
by edges, and it is established in accordance with different SE problem
formulations:

  * $\bullet$

Linear SE: Variable nodes are used to calculate a $s$-dimensional node
embedding for all real and imaginary parts of the bus voltages,
$\Re(\mathscr{V}_{i})$ and $\Im(\mathscr{V}_{i})$, $i=1,\dots,n$, using which
state variable predictions can be generated. Factor nodes, two per each
measurement phasor, serve as inputs for the measurement values, variances, and
covariances, also given in rectangular coordinates. These values are then
transformed and sent to variable nodes via GNN message passing. Unlike the
approximative WLS SE defined in 2.2, which neglects measurement covariances,
the GNN includes them, leading to accurate solutions without increasing the
computational complexity. The upper subfigure of Fig. 4.1 illustrates a two-
bus power system, with a PMU on the first bus, containing one voltage and one
current phasor measurement. The nodes connected by full lines represent the
corresponding factor graph in the lower subfigure.

  * $\bullet$

Nonlinear SE: In this case, pairs of variable nodes generate $s$-dimensional
node embeddings for magnitudes ${V}_{i}$ and phase angles ${\theta}_{i}$ of
complex bus voltages $\mathscr{V}_{i}=V_{i}\mathrm{e}^{\mathrm{j}\theta_{i}}$.
The inputs to the factor nodes are the values and variances of both phasor and
legacy measurements. Phasor measurements are expressed in polar coordinates,
which eliminates any correlation between measurement errors. Therefore, in
contrast to linear SE, in nonlinear SE measurement covariances do not need to
be included as inputs into factor nodes. When creating the factor graph from
the bus-branch power system model, each phasor measurement generates two
factor nodes, while each legacy measurement generates one factor node. As an
example, in the upper subfigure of Fig. 4.2 we consider a simple two-bus power
system, in which we placed one voltage phasor measurement on the first bus and
one legacy voltage magnitude measurement on the second bus. Additionally, we
placed one current phasor measurement and one legacy active power flow
measurement on the branch connecting the two nodes. The factor graph of this
simple power system consists of the generated factor and variable nodes,
connected by full-line edges, as shown in the lower subfigure of Fig. 4.2.

$\mathscr{V}_{m1}\mathrm{=}\Re({\mathscr{V}}_{m1})+\mathrm{j}\Im({\mathscr{V}}_{m1})$${\mathscr{V}}_{1}\mathrm{=}\Re{(\mathscr{V}_{1})}+\mathrm{j}\Im{(\mathscr{V}_{1})}$Bus
1$\mathscr{I}_{12}\mathrm{=}\Re{(\mathscr{I}_{12})}+\mathrm{j}\Im{(\mathscr{I}_{12})}$${\mathscr{V}}_{2}\mathrm{=}\Re{(\mathscr{V}_{2})}+\mathrm{j}\Im{(\mathscr{V}_{2})}$Bus
2 (a)

$f_{\Re({\mathscr{V}}_{1})}$$f_{\Im({\mathscr{V}}_{1})}$$\Re({\mathscr{V}}_{1})$$\Im({\mathscr{V}}_{1})$$f_{\Re({\mathscr{I}}_{12})}$$f_{\Im({\mathscr{I}}_{12})}$$\Re({\mathscr{V}}_{2})$$\Im({\mathscr{V}}_{2})$
(b)

Figure 4.1: Subfigure (a) shows a simple two-bus power system with two phasor
measurements from a PMU placed at the bus 1. Subfigure (b) displays the
corresponding factor graph (full-line edges) and augmented factor graph (all
edges). Variable nodes are depicted as circles, and factor nodes are as
squares.

${\mathscr{V}}_{m2}\mathrm{=}V_{m1}\text{e}^{\text{j}\theta_{m1}}$$\mathscr{V}_{1}=V_{1}\text{e}^{\text{j}\theta_{1}}$Bus
1$\mathscr{I}_{12}\mathrm{=}I_{12}\text{e}^{\text{j}\theta_{I_{12}}}$$P_{12}$$\mathscr{V}_{2}=V_{2}\text{e}^{\text{j}\theta_{2}}$Bus
2$V_{m2}$ (a)

$f_{V_{m1}}$$f_{\theta_{m1}}$$V_{1}$$\theta_{1}$$f_{I_{12}}$$f_{\theta_{I_{12}}}$$V_{2}$$\theta_{2}$$f_{V_{m2}}$$f_{P_{12}}$
(b)

Figure 4.2: Subfigure (a) shows a simple two-bus power system containing a PMU
at the bus 1, one legacy active power flow measurement, and one legacy voltage
magnitude measurement at the bus 2. Subfigure (b) displays the corresponding
factor graph (full-line edges) and augmented factor graph (all edges).
Variable nodes are represented as circles, and factor nodes are depicted as
squares, coloured differently to distinguish between phasor and legacy
measurements.

Unlike approaches that apply GNNs over the bus-branch power system model, such
as in [55, 57], the using GNNs over factor-graph-like topology allows for the
incorporation of various types and quantities of measurements on both power
system buses and branches. The ability to simulate the addition or removal of
various measurements can be easily achieved by adding or removing factor nodes
from any location in the graph. In contrast, using a GNN over the bus-branch
power system model would require allocating a single input vector to each bus
that includes all potential measurement data for that bus and its neighbouring
branches. This can cause problems, such as having to fill elements in the
input vector with zeros when not all possible measurements are available, and
making the output sensitive to the order of measurements in the input vector.
This can make it difficult to accurately model the system and generate
reliable results.

Augmenting the factor graph topology by connecting the variable nodes in the
$2$-hop neighbourhood significantly improves the model’s prediction quality in
unobservable scenarios. This is because the graph should remain connected even
when we remove factor nodes to simulate measurement loss. This will allow the
messages to be still propagated in the whole $K$-hop neighbourhood of the
variable node. In other words, a factor node corresponding to a branch current
measurement can be removed while still preserving the physical connection
between the power system buses. This requires an additional set of trainable
parameters for the variable-to-variable message function. Although the
augmented factor graph displayed with both full and dashed lines in Figs. 4.1
and 4.2 is not a factor graph because it is no longer bipartite, we will still
refer to the nodes as factor and variable nodes for simplicity.

To achieve better representation of node’s neighbourhood structure, we perform
variable node feature augmentation using binary index encoding. Since variable
nodes have no additional input features, this encoding allows the GNN to
better capture the relationships between nodes. Compared to one-hot encoding
used in [30], binary index encoding significantly reduces the number of input
neurons and trainable parameters, as well as training and inference time.

####  4.2 Proposed GNN Architecture

Since the proposed GNN operates on a heterogeneous graph, we use two different
types of GNN layers: one for aggregation in factor nodes, and one for variable
nodes. These layers, denoted as
$\operatorname{Layer^{f}}(\cdot|\theta^{\operatorname{Layer^{f}}}):\mathbb{R}^{\textrm{deg}(f)\cdot
s}\mapsto\mathbb{R}^{s}$ and
$\operatorname{Layer^{v}}(\cdot|\theta^{\operatorname{Layer^{v}}}):\mathbb{R}^{\textrm{deg}(v)\cdot
s}\mapsto\mathbb{R}^{s}$, have their own sets of trainable parameters
$\theta^{\operatorname{Layer^{f}}}$ and $\theta^{\operatorname{Layer^{v}}}$,
allowing their message, aggregation, and update functions to be learned
separately. Additionally, we use different sets of trainable parameters for
variable-to-variable and factor-to-variable node messages,
$\textrm{Message}\textsuperscript{f$\rightarrow$v}(\cdot|\theta^{\textrm{Message}\textsuperscript{f$\rightarrow$v}}):\mathbb{R}^{2s}\mapsto\mathbb{R}^{u}$
and
$\textrm{Message}\textsuperscript{v$\rightarrow$v}(\cdot|\theta^{\textrm{Message}\textsuperscript{v$\rightarrow$v}}):\mathbb{R}^{2s}\mapsto\mathbb{R}^{u}$,
in the $\operatorname{Layer^{v}}(\cdot|\theta^{\operatorname{Layer^{v}}})$
layer. In both GNN layers, we use two-layer feed-forward neural networks as
message functions, single layer neural networks as update functions and the
attention mechanism in the aggregation function. Furthermore, we apply a two-
layer neural network
$\operatorname{Pred}(\cdot|\theta^{\operatorname{Pred}}):\mathbb{R}^{s}\mapsto\mathbb{R}$
to the final node embeddings $\mathbf{h}^{K}$ of variable nodes only, to
create state variable predictions $\mathbf{x^{pred}}$. For factor and variable
nodes with indices $f$ and $v$, neighbourhood aggregation and state variable
prediction can be described as:

$\begin{gathered}\mathbf{h_{v}}^{k}=\operatorname{Layer^{v}}(\\{\mathbf{h_{i}}^{k-1}|i\in\mathcal{N}_{v}\\}|\theta^{\operatorname{Layer^{v}}})\\\
\mathbf{h_{f}}^{k}=\operatorname{Layer^{f}}(\\{\mathbf{h_{i}}^{k-1}|i\in\mathcal{N}_{f}\\}|\theta^{\operatorname{Layer^{f}}})\\\
{x_{v}}^{pred}=\operatorname{Pred}(\mathbf{h_{v}}^{K}|\theta^{\operatorname{Pred}})\\\
k\in\\{1,\dots,K\\},\end{gathered}$ (4.1)

where $\mathcal{N}_{v}$ and $\mathcal{N}_{f}$ denote the $1$-hop
neighbourhoods of the nodes $v$ and $f$. All the trainable parameters $\theta$
of the GNN are updated by applying gradient descent, using backpropagation, to
a loss function calculated over a mini-batch of graphs. This loss function is
the mean-squared difference between the predicted state variables and their
corresponding ground-truth values:

$\begin{gathered}L(\theta)=\frac{1}{2nB}\sum_{i=1}^{2nB}({{x_{i}}^{pred}}-{{x_{i}}^{label}})^{2}\\\
\theta=\\{\theta^{\operatorname{Layer^{v}}}\mathop{\cup}\theta^{\operatorname{Layer^{f}}}\mathop{\cup}\theta^{\operatorname{Pred}}\\}\\\
\theta^{\operatorname{Layer^{v}}}=\\{\theta^{\textrm{Message}\textsuperscript{f$\rightarrow$v}}\mathop{\cup}\theta^{\textrm{Message}\textsuperscript{v$\rightarrow$v}}\mathop{\cup}\theta^{\textrm{Aggregate}\textsuperscript{v}}\mathop{\cup}\theta^{\textrm{Update}\textsuperscript{v}}\\}\\\
\theta^{\operatorname{Layer^{f}}}=\\{\theta^{\textrm{Message}\textsuperscript{v$\rightarrow$f}}\mathop{\cup}\theta^{\textrm{Aggregate}\textsuperscript{f}}\mathop{\cup}\theta^{\textrm{Update}\textsuperscript{f}}\\},\end{gathered}$
(4.2)

where the total number of variable nodes in a graph is $2n$, and the number of
graphs in the mini-batch is $B$. In this work, we chose the loss function for
training the GNN based on the fundamental SE problem, where state variables
are obtained from available measurement information. However, if there is a
requirement to include additional constraints in the SE calculation, it is
possible to achieve this by adding new terms to the loss function defined in
(4.2). For example, the loss function can be augmented with the power balance
error at each bus where the constraints are imposed, in addition to minimizing
the prediction error from the labels. A similar approach has been proposed in
[54], where the power flow problem is solved using GNN by minimizing the power
balance errors at each bus. Adding additional constraints (e.g., zero
injection constraints) to the GNN SE loss function can improve the SE results,
especially in distributed power systems with limited measurement coverage.

Fig. 4.3 shows the high-level computational graph for the output of a variable
node from the augmented factor graph given in Fig. 4.1. For simplicity, only
one unrolling of the neighbourhood aggregation is shown, as well as only the
details of the parameters $\theta^{\operatorname{Layer^{v}}}$.

$\operatorname{Layer^{f}}$$\operatorname{Layer^{v}}$$\mathbf{h_{f}}^{K-1}$$\mathbf{h_{v2}}^{K-1}$$\mathbf{h_{v}}^{K}$$\operatorname{Layer^{v}}$$\operatorname{Pred}$Loss$...$$...$$...$$...$$output$$label$
(a)

Messagef$\rightarrow$vMessagev$\rightarrow$vGATvUpdatev$\mathbf{h_{v2}}^{K-1};\mathbf{h_{v}}^{K-1}$$\mathbf{h_{f}}^{K-1};\mathbf{h_{v}}^{K-1}$$...$$...$$\mathbf{h_{v}}^{K}$$\mathbf{h_{v}}^{K-1}$
(b)

Figure 4.3: Subfigure (a) shows a high-level computational graph that starts
with the loss function for the output of a variable node $v$. Subfigure (b)
depicts the detailed structure of a single GNN $\operatorname{Layer^{v}}$.
Functions with trainable parameters are highlighted in yellow.

#####  4.2.1 Computational Complexity and Distributed Inference

Because the node degree in the power system graph does not increase with the
total number of nodes, the same is true for the node degrees in the augmented
factor graph. This means that the inference time per variable node remains
constant, as it only requires information from the node’s $K$-hop
neighbourhood, whose size also does not increase with the total number of
nodes. This implies that the computational complexity of inferring all state
variables is $\mathcal{O}(n)$. To avoid the over-smoothing problem in GNNs
[134], a small value is assigned to $K$, thus not affecting the overall
computational complexity of the inference.

To make the best use of the proposed approach for large-scale power systems,
the inference should be performed in a computationally and geographically
distributed manner. This is necessary because the communication delays between
the PMUs and the central processing unit can hinder the full utilization of
the PMUs’ high sampling rates. The distributed implementation is possible as
long as all the measurements within a node’s $K$-hop neighbourhood in the
augmented factor graph are fed into the computational module that generates
the predictions. For any arbitrary $K$, the GNN inference method only requires
measurements that are physically located within the $\lceil K/2\rceil$-hop
neighbourhood of the power system bus.

####  4.3 Numerical results

In this section, we conduct comprehensive numerical tests to evaluate the
effectiveness of proposed augmented factor graph-based GNN approaches for
linear and nonlinear SE problems. We used the IGNNITION framework [151] for
building and utilizing GNN models, with the hyperparameters presented in Table
4.1, the first three of which were obtained with the grid search
hyperparameter optimization using the Tune tool [152]. All the results
presented in this section are normalized using the corresponding nominal
voltages in the test power systems and a base power of 100 MVA.

Table 4.1: List of GNN hyperparameters. Hyperparameters | Values
---|---
Node embedding size $s$ | $64$
Learning rate | $4\times 10^{-4}$
Minibatch size $B$ | $32$
Number of GNN layers $K$ | $4$
Activation functions | ReLU
Gradient clipping value | $0.5$
Optimizer | Adam
Batch normalization type | Mean

#####  4.3.1 Linear State Estimation

To evaluate the proposed GNN-based linear SE, we create various scenarios
using the IEEE 30-bus system, the IEEE 300-bus system, and the ACTIVSg
2000-bus system [153], on which the GNN model is trained and tested. Training,
validation, and test sets are obtained using WLS solutions of the system
described in (2.9) to label various samples of input measurements.
Measurements are obtained by adding Gaussian noise to the exact power flow
solutions, with each calculation performed using a different, randomly sampled
load profile to capture a wide range of power system states. Due to the strong
interpolation and extrapolation abilities of GNNs [154], our method of
randomly sampling from a wide diversity of loads for training examples is
effective for generalizing the GNN algorithm for state estimation under
varying load conditions. GNN models are tested in three different situations:
i) optimal number of PMUs (minimal measurement redundancy, for which the WLS
SE offers a solution); ii) underdetermined scenarios; iii) scenarios with
maximal measurement redundancy. We also compare the proposed approach with
more conventional deep neural network (DNN)-based SE algorithms and assess its
scalability, sample efficiency, and robustness to outliers.

###### Power System With Optimally Placed PMUs

In this subsection, we conduct a series of experiments on the IEEE 30-bus
power system, using measurement variances of $10^{-5}$, $10^{-3}$, and
$10^{-1}$ for the creation of the training, validation, and test sets. The
number and positions of PMUs are fixed and determined using the optimal PMU
placement algorithm [155], which finds the smallest set of PMUs that make the
system observable. This algorithm has resulted in a total of $10$ PMUs and
$50$ measurement phasors, $10$ of which are voltage phasors and the rest are
current phasors.

Table 4.2 shows the $100$-sample test set results for all the experiments on
the IEEE 30-bus power system, in the form of average mean square errors (MSEs)
between the GNN predictions and the test set labels. These results are
compared with the average MSE between the labels and the approximate WLS SE
solutions defined in Sec. 2.2. The results show that for systems with
optimally placed PMUs and low measurement variances, GNN predictions have very
small deviations from the exact WLS SE, although they are outperformed by the
approximate WLS SE. For higher measurement variances, GNN has a lower
estimation error than the approximate WLS SE, while also having lower
computational complexity in all cases.

Table 4.2: Comparison of GNN and approximative SE test set MSEs for various measurement variances. Variances | GNN | Approx. SE
---|---|---
$10^{-5}$ | $2.48\times 10^{-6}$ | $1.87\times 10^{-8}$
$10^{-3}$ | $8.21\times 10^{-6}$ | $2.25\times 10^{-6}$
$10^{-1}$ | $7.47\times 10^{-4}$ | $2.27\times 10^{-3}$

In Fig. 4.4, we present the predictions and labels for each of the variable
nodes for one of the samples from the test set created with measurement
variance $10^{-5}$. The results for the real and imaginary parts of the
complex node voltages (shown in the upper and lower plots, respectively)
indicate that GNNs can be used as accurate SE approximations.

0.9511.05$\displaystyle\Re({\mathscr{V}}_{i})$
[p.u.]PredictionsLabels051015202530-0.3-0.2-0.10Bus index
($i$)$\displaystyle\Im({\mathscr{V}}_{i})$ [p.u.] Figure 4.4: GNN predictions
and labels for one test example with optimally placed PMUs.

###### Performance in a Partially Observable Scenario

To further assess the robustness of the proposed model, we test it by
excluding several measurement phasors from the previously used test samples
with optimally placed PMUs, resulting in an underdetermined system of
equations that describes the SE problem. These scenarios are relevant even at
higher levels of system redundancy, where partial grid observability can occur
due to multiple component (PMU and communication link) failures caused by
natural disasters or cyberattacks. To exclude a measurement phasor from the
test sample, we remove its real and imaginary parts from the input data, which
is equivalent to removing two factor nodes from the augmented factor graph. We
use the previously used 100-sample test set to create a new test set by
removing selected measurement phasors from each sample while preserving the
same labels obtained as SE solutions of the system with all the measurements
present. As an example, we consider a scenario where two neighbouring PMUs
fail to deliver measurements to the state estimator. In this case, all eight
measurement phasors associated with the removed PMUs are excluded from the GNN
inputs. The average MSE for the test set of 100 samples created by removing
these measurements from the original test set used in this section is
$3.45\times 10^{-3}$. The predictions and labels for a single test set sample,
per variable node index, are shown in Figure 4.5. The figure includes vertical
black dashed lines that indicate the indices of unobserved buses 17 and 18.
These buses have higher prediction errors due to the lack of input measurement
data. Neighbouring buses that are not unobserved, but are affected by
measurement loss, are indicated with vertical green lines and have lower
prediction errors. It can be observed that significant deviations from the
labels occur for some of the neighbouring buses, while the GNN predictions are
a decent fit for the remaining node labels. This demonstrates the proposed
model’s ability to sustain error in the neighbourhood of the malfunctioning
PMU, as well as its robustness in scenarios that cannot be solved using
standard WLS approaches. A possible explanation for the higher susceptibility
to errors in the imaginary parts of the voltage is related to their variance
in the training set. The variance of real parts of voltages is $6.6\times
10^{-4}$, while the variance of imaginary parts of voltages is $4.8\times
10^{-3}$. This indicates that the imaginary parts of voltages have a higher
variability and may therefore be more prone to errors in the prediction model.

0.9511.05$\displaystyle\Re({V}_{i})$ [p.u.]PredictionsLabelsUnobserved
busesAffected buses051015202530-0.3-0.2-0.10Bus index
($i$)$\displaystyle\Im({V}_{i})$ [p.u.] Figure 4.5: GNN predictions and labels
for one test example with phasors from two neighbouring PMUs removed. Vertical
black lines indicate unobserved buses, while green lines represent buses that
are affected by the loss of measurement data.

###### Comparison With the Feed-Forward Deep Neural Network-Based State
Estimation

The main goal of this subsection is to compare the performance of the proposed
GNN-based SE approach with a state-of-the-art deep learning-based approach on
a variety of power systems. We used a 6-layer feed-forward DNN model, proposed
by [48], with the same number of neurons in each layer as the number of input
measurement scalars. This DNN architecture has similar performance as the best
architecture proposed in the same work obtained by unrolling the iterative
nonlinear SE solver, which cannot be applied directly to the linear SE problem
we are considering.

We tested both approaches on the IEEE 30-bus, IEEE 118-bus, IEEE 300-bus, and
the ACTIVSg 2000-bus power systems [153], with measurement variances set to
$10^{-5}$. In contrast to previous examples, we used maximal measurement
redundancies, ranging from $3.73$ to $4.21$. We provide a comparison of the
number of trainable parameters for both GNN and DNN models for various power
system sizes, which is often left out in similar analyses. To compare sample
efficiencies between the GNN and DNN approaches, separate models for each of
the test power systems were trained using smaller and larger datasets,
containing 10 and 10000 samples. The results for all test power systems are
presented in Table 4.3. The first four rows show the 100-sample test set MSE
for GNN and DNN models trained using smaller and larger datasets. The last two
rows of the table show the number of trainable parameters for both approaches,
depending on the power system size.

Table 4.3: A comparison of the performance of GNN and DNN models trained on
different training set sizes, as measured by test set MSE and the number of
trainable parameters.

| IEEE 30 | IEEE 300 | ACTIVSg 2000
---|---|---|---
Small training set GNN | $4.73\times 10^{-6}$ | $5.94\times 10^{-5}$ | $5.08\times 10^{-4}$
Small training set DNN | $9.29\times 10^{-4}$ | $5.92\times 10^{-3}$ | $4.77\times 10^{-3}$
Large training set GNN | $2.48\times 10^{-6}$ | $6.62\times 10^{-6}$ | $3.91\times 10^{-4}$
Large training set DNN | $6.28\times 10^{-6}$ | $2.91\times 10^{-3}$ | $2.61\times 10^{-3}$
GNN parameters | $4.99\times 10^{4}$ | $4.99\times 10^{4}$ | $4.99\times 10^{4}$
DNN parameters | $3.16\times 10^{5}$ | $3.15\times 10^{7}$ | $1.77\times 10^{9}$

The results show that the GNN approaches result in higher overall accuracy
compared to the corresponding DNN approaches for all the power system and
training set sizes. Furthermore, the number of trainable parameters (i.e., the
model size) is constant111More precisely, the number of trainable parameters
in the proposed GNN model remains nearly constant as the number of buses in
the power system increases. This effect would only be noticeable for larger
power systems. The only exception is the number of input neurons for the
binary index encoding of the variable nodes, which grows logarithmically with
the number of variable nodes. However, this increase is insignificant compared
to the total number of GNN parameters. and relatively low for GNN models,
because the number of neurons in a layer is constant regardless of power
system size. In contrast, the number of parameters grows quadratically for DNN
models, because the number of neurons in a layer grows linearly with the input
size, resulting in quadratic growth of the trainable parameter matrices. When
expressed in computer memory units, the GNN models we used had a significantly
smaller memory footprint, taking up only 0.19 MB. In comparison, the DNN model
used for the ACTIVSg 2000 power system required a much larger 6.58 GB of
memory, resulting in more challenging training and inference processes. The
high number of trainable parameters required by DNN models increases their
storage requirements, increases the dimensionality of the training process,
and directly affects the inference speed and computational complexity. Since
the proposed GNNs have a linear computational complexity in the prediction
process, one training iteration of GNN also has a linear computational
complexity. In contrast, one training iteration of DNN-based SE would have at
least quadratic computational complexity per training iteration, making the
overall training process significantly slower. To recall, the reason why the
GNN has a constant number of parameters and generates predictions with linear
computational complexity is that it takes measurements from a limited $K$-hop
neighbourhood for every node, regardless of the size of the power system.

The results indicate that the quality of GNN and DNN model predictions
improves with more training data. However, compared to GNNs, DNN models
performed significantly worse on smaller datasets, suggesting that they are
less sample efficient and more prone to overfitting due to their larger number
of hyperparameters. While we used randomly generated training sets in the
experiments, narrowing the learning space by selecting training samples based
on historical load consumption data could potentially result in even better
performance with small datasets.

The use of GNNs for power systems analysis has several additional advantages
over using DNNs. One advantage is flexibility: spatial GNNs can produce
results even if the input power system topology changes, whereas conventional
DNN methods are trained and tested on the same topology of the power system.
For example, if some measurements are removed from the inputs (as discussed in
Subsection 4.3.1), a DNN would require retraining from scratch with the new
topology. GNNs also have some theoretical advantages over other deep learning
methods in that they are permutation invariant and equivariant. This means
that the output of the GNN is the same regardless of the order in which the
inputs are presented, and that the GNN output changes in a predictable and
consistent way when the inputs are transformed. This property is useful for
problems like SE, where the order of the nodes and edges is not important and
the system can undergo topological changes. In addition, GNNs incorporate
topology information into the learning process by design, whereas many other
deep learning methods in power systems use node-level data as inputs while
ignoring connectivity information. Finally, unlike most deep learning methods,
spatial GNNs can be distributed for evaluation among edge devices.

Similar conclusions can be drawn when comparing GNNs with recurrent and
convolutional neural networks for similar power systems’ analysis algorithms,
as they also require information from the entire power system as input.
Overall, this comparison highlights the potential advantages of using GNNs for
power system modelling and analysis.

###### Robustness to Outliers

To assess how well the proposed model can handle outliers in input data, we
carried out experiments on two separate test sets, each containing samples
with different degrees of outlier intensity. The experiments followed the
setup described in Sec. 4.3.1, with the test samples initially generated using
a measurement variance of $10^{-5}$. We replaced one of the existing
measurements in each test sample with a randomly generated value, using a
variance of either $\sigma_{1}=1.6$ or $\sigma_{2}=1.6\cdot 10^{2}$. WLS SE
solutions without outliers in the inputs are used as ground-truth values. The
results, shown in Table 4.4, indicate the performance of four different
approaches, as well as the WLS SE and the approximative WLS SE algorithm on
the same test sets.

The first approach, which uses the already trained model from Sec. 4.3.1,
results in the highest prediction error on tests set with outliers. The
primary factor contributing to the MSE, particularly in the test set with
outliers generated using larger variances, are significant mismatches from the
ground-truth values in the $K$-hop neighbourhood of the outlier. This occurs
because the ReLU activation function does not constrain its inputs during
neighbourhood aggregation. To address this problem, we propose a second
approach where we train a GNN model with the same architecture as the previous
one, but which uses the tanh (hyperbolic tangent) activation function instead
of ReLU. As presented in Table 4.4, this approach results in a significantly
lower test set MSE for outliers generated using larger variances compared to
the proposed GNN with ReLU activations, WLS SE, and the approximative WLS SE.
The saturation effect of tanh prevents high values stemming from outliers from
propagating through the GNN, but also reduces the training quality due to the
vanishing gradient problem. Specifically, all the experiments we conducted
under the same conditions for GNN with tanh activations required more epochs
to converge to a solution with lower prediction quality compared to the GNN
with ReLU activations. As a third approach, we propose training a GNN model
with ReLU activations on a dataset in which half of the samples contain
outliers, which are generated in the same manner as the test samples used in
this subsection. This approach turned out to be the most effective for both
test cases because the GNN learns to neutralize the effect of unexpected
inputs from the dataset examples while maintaining accurate predictions in the
absence of outliers in the input data. To confirm the validity of this
methodology, we trained the DNN introduced in the subsection 4.3.1 using the
same datasets containing outliers. DNN was able to neutralize the effect of
unexpected inputs because the input power system is small, resulting in the
second-best approach in terms of robustness to outliers, trailing the GNN
trained on the dataset containing outliers.

In summary, as expected, all methods produced better results for the test set
containing outliers with lower variances, while the GNN trained with outliers
demonstrated the best performance for both higher and lower variance outliers.
We note that these are only preliminary efforts to make the GNN model robust
to outliers, and that future research could combine ideas from standard bad
data processing methods in SE with the proposed GNN approach.

Table 4.4:  A comparison of the results of various approaches for two test sets with different degrees of outlier intensity. Approach |  | Test set MSE
---
$\mathbf{\sigma_{1}=1.6}$
| Test set MSE
---
$\mathbf{\sigma_{2}=1.6\cdot 10^{2}}$
GNN | $9.48\times 10^{-3}$ | $1.60\times 10^{3}$
GNN with tanh | $6.87\times 10^{-3}$ | $2.39\times 10^{-2}$
GNN trained with outliers | $4.44\times 10^{-6}$ | $7.99\times 10^{-6}$
DNN trained with outliers | $1.06\times 10^{-5}$ | $4.21\times 10^{-5}$
WLS SE | $1.43\times 10^{-3}$ | $1.41\times 10^{-1}$
Approximative WLS SE | $1.39\times 10^{-3}$ | $1.35\times 10^{-1}$

#####  4.3.2 Scalability and Sample Efficiency Analysis of Linear State
Estimation

In the previous subsection, the GNN model for linear SE demonstrated good
approximation capabilities under normal operating conditions and performed
well in unobservable and underdetermined scenarios. This subsection extends
the previous one in the following ways:

  * $\bullet$

We conduct an empirical analysis to investigate how the same GNN architecture
could be used for power systems of various sizes. We use the IEEE 30-bus
system, the IEEE 118-bus system, the IEEE 300-bus system, and the ACTIVSg
2000-bus system [153], with measurements placed so that measurement redundancy
is maximal. Our main assumption is that the local properties of the graphs in
these systems are similar, leading to local neighbourhoods with similar
structures which can be represented using the same embedding space size and
the same number of GNN layers.

  * $\bullet$

To evaluate the sample efficiency of the GNN model, we run multiple training
experiments on different sizes of training sets. Additionally, we assess the
scalability of the model by training it on various power system sizes and
evaluating its accuracy, training convergence properties, inference time, and
memory requirements. For this purpose we create training sets containing 10,
100, 1000, and 10000 samples for each of the mentioned power systems, while
the GNN models are tested on sets comprising 100 samples.

  * $\bullet$

As a side contribution, the proposed GNN model is tested in scenarios with
high measurement variances, using which we simulate phasor misalignments due
to communication delays [156]. While this is usually simulated by using
variance values that increase over time, as an extreme scenario we fix the
measurement variances to a high value of $5\times 10^{-1}$.

It is important to note that the conclusions that will be made apply to GNN
based-nonlinear SE as well.

###### Properties of Power System Augmented Factor Graphs

For all four test power systems, we create augmented factor graphs using the
methodology described in Section 4.2. Fig. 4.6 illustrates how the properties
of the augmented factor graphs, such as average node degree, average path
length, average clustering coefficient, along with the system’s maximal
measurement redundancy, vary across different test power systems.

30118300200002468101214Number of busesPower systemgraph
propertiesRedundancyAvg. DegreeAvg. Path LengthAvg. Cluster Coeff. Figure 4.6:
Properties of augmented factor graphs along with the system’s measurement
redundancy for different test power systems, labelled with their corresponding
number of buses.

The average path length is a property that characterizes the global graph
structure, and it tends to increase as the size of the system grows. However,
as a design property of high-voltage networks, the other graph properties such
as the average node degree, average clustering coefficient, as well as maximal
measurement redundancy do not exhibit a clear trend of change with respect to
the size of the power system. This suggests that the structures of local,
$K$-hop neighbourhoods within the graph are similar across different power
systems, and that they contain similar factor-to-variable node ratio.
Consequently, it is reasonable to use the same GNN architecture (most
importantly, the number of GNN layers and the node embedding size) for all
test power systems, regardless of their size. In this way, the proposed model
achieves scalability, as it applies the same set of operations to the local,
$K$-hop neighbourhoods of augmented factor graphs of varying sizes without
having to adapt to each individual case.

###### Training Convergence Analysis

First, we analyse the training process for the IEEE 30-bus system with four
different sizes of the training set. As mentioned in 4.2, the training loss is
a measure of the error between the predictions and the ground-truth values for
data samples used in the training process. The validation loss, on the other
hand, is a measure of the error between the predictions and the ground-truth
values on a separate validation set. In this analysis, we used a validation
set of 100 samples.

130609012015000.050.1EpochValidation loss$10$$100$$1000$$10000$ Figure 4.7:
Validation losses for trainings on four different training set sizes.

The training losses for all the training processes converged smoothly, so we
do not plot them for the sake of clarity. Figure 4.7 shows the validation
losses for 150 epochs of training on four different training sets. For smaller
training sets, the validation loss decreases initially but then begins to
increase, which is a sign of overfitting. In these cases, a common practice in
machine learning is to select the model with the lowest validation loss value.
As it will be shown in 4.3.2, the separate test set results for models created
using small training sets are still satisfactory. As the number of samples in
the training set increases, the training process becomes more stable. This is
because the model has more data to learn from and is therefore less prone to
overfitting.

Next, in Table 4.5, we present the training results for the other power
systems and training sets of various sizes. The numbers in the table represent
the number of epochs after which either the validation loss stopped changing
or began to increase. Similarly to the experiments on the IEEE 30-bus system,
the trainings on smaller training sets exhibited overfitting, while others
converged smoothly. For the former, the number in the table indicates the
epoch at which the validation loss reached its minimum and stopped improving.
For the latter, the number in the table represents the epoch when there were
five consecutive validation loss changes less than $10^{-5}$.

Table 4.5: Epoch until validation loss minimum for various power systems and
training set sizes.

Power system | IEEE 118 | IEEE 300 | ACTIVSg 2000
---|---|---|---
10 samples | $61$ | $400$ | $166$
100 samples | $38$ | $84$ | $200$
1000 samples | $24$ | $82$ | $49$
10000 samples | $12$ | $30$ | $15$

Increasing the size of the training set generally results in a lower number of
epochs until the validation loss reaches its minimum. However, the epochs
until the validation loss reaches its minimum vary significantly between the
different power systems. This could be due to differences in the complexity of
the systems or the quality of the data used for training.

###### Accuracy Assessment

Fig. 4.8 reports the MSEs between the predictions and the ground-truth values
on 100-sample sized test sets for all trained models and the approximate WLS
SE. These results indicate that even the GNN models trained on small datasets
outperform the approximate WLS SE, except for the models trained on the IEEE
30-bus system with 10 and 100 samples. These results suggest that the quality
of the GNN model’s predictions and the generalization capabilities improve as
the amount of training data increases, and the models with the best results
(highlighted in bold) have significantly smaller MSEs compared to the
approximate WLS SE. We assume that using carefully selected training samples
based on historical load consumption data instead of randomly generated ones
could potentially lead to even better results with small datasets.

Approx. SE (baseline)GNN SE

$10$$10^{2}$$10^{3}$$10^{4}$$0$$0.02$$0.04$Number of training samplesTest set
MSE(a) IEEE 30

$10$$10^{2}$$10^{3}$$10^{4}$$0$$0.02$$0.04$Number of training samplesTest set
MSE(b) IEEE 118

$10$$10^{2}$$10^{3}$$10^{4}$$0$$0.02$$0.04$Number of training samplesTest set
MSE(c) IEEE 300

$10$$10^{2}$$10^{3}$$10^{4}$$0$$0.02$$0.04$Number of training samplesTest set
MSE(d) ACTIVSg 2000

Figure 4.8: Test set results for various power systems and training set sizes.

###### Inference Time

The plot in Fig. 4.9 shows the ratio of execution times between WLS SE and GNN
SE inference as a function of the number of buses in the system. These times
are measured on a test set of 100 samples. The difference in computational
complexity between GNN, with its linear complexity, and WLS, with more than
quadratic complexity, becomes more apparent on the plot as the number of buses
increases. From the results, it is clear that GNN significantly outperforms
WLS in terms of inference time on larger power systems.

Unlike with GNNs, the hardware implementation of matrix operations in WLS is a
well-established field. However, the hardware implementation of GNNs is an
active area of research, and it is possible that inference times may improve
even further in the future [157].

301183002000121064Number of busesInference time ratio Figure 4.9: A ratio of
the execution times for WLS SE and GNN SE inference on a test set of 100
samples, as a function of the power system size.

#####  4.3.3 Nonlinear State Estimation

Finally, in this subsection, we present the numerical results of the proposed
approach for the nonlinear SE problem formulation. We describe the GNN model’s
training process and test the trained model on various examples to validate
its accuracy, and its robustness under measurement data loss due to
communication failure and cyberattacks in the form of malicious data
injections.

###### Demonstration of prediction accuracy

We conducted separate training experiments for IEEE 30 and IEEE 118-bus test
cases, for which we generated a training set containing 10000 samples and
validation and test sets containing 100 samples each. Similarly to the linear
SE approach, each sample is created by randomly sampling the active and
reactive power injections and solving the power flow problem. Measurement
values are created by adding Gaussian noise to the power flow solutions, and
the nonlinear SE is solved by GN to label the input measurement set in each
sample. We used a Gaussian noise variance of $10^{-5}$ for phasor
measurements, $10^{-3}$ for bus voltage magnitude and active and reactive
power flow legacy measurements, and $10^{-1}$ for active and reactive
injection legacy measurements.

For the IEEE 30-bus test case, we placed $100$ legacy measurements and three
PMUs (i.e., three bus voltage phasors and eight branch current phasors) in
each sample, resulting in $2.03$ measurement redundancy. The trained model
performed well on the test set, with the average test set mean square error of
$1.233\times 10^{-5}$ between predictions and ground truth labels; the average
test set MSE for voltage magnitudes of $5.221\times 10^{-6}$; the average test
set MSE for voltage angles of $1.944\times 10^{-5}$. Fig. 4.10 shows the
average test MSE per each bus, where the upper plot corresponds to voltage
magnitudes and the lower one to voltage angles.

For the IEEE 118-bus test case, we placed $500$ legacy measurements and seven
PMUs (i.e., seven bus voltage phasors and $26$ branch current phasors) in each
sample, resulting in $2.39$ measurement redundancy. The average test set mean
square error equals $2.038\times 10^{-5}$, with the average test set MSE for
voltage magnitudes of $1.572\times 10^{-5}$ and the average test set MSE for
voltage angles of $2.505\times 10^{-5}$. Based on the insights from both
experiments, we can conclude that the proposed GNN model is a good
approximator of the nonlinear SE solved by GN.

$10^{-12}$$10^{-9}$$10^{-6}$$10^{-3}$MSE of
$V_{i}$151015202530$10^{-12}$$10^{-9}$$10^{-6}$$10^{-3}$Bus Index $i$MSE of
${\theta}_{i}$ Figure 4.10: The test set MSE between the predictions and the
labels per each bus for voltage magnitudes and angles in the IEEE 30-bus test
case.

###### Robustness to Loss of Input Data

Next, we observe predictions of the GNN models previously trained on IEEE 30
and IEEE 118-bus test data when exposed to the loss of input data caused by
communication failures or measurement device malfunctions. We simulate the
described cases by randomly removing a percentage of all input measurements,
ranging from $0\%$ to $95\%$ with a step of $5\%$. We create 20 test sets per
IEEE test case, each containing samples with the same percentage of excluded
measurements, and show the average test set MSEs in Fig. 4.11. Proposed GNN
models yields predictions in all examples, with an expected growing trend in
MSE as the number of excluded measurements increases. In comparison, the GN
method could not provide a solution for many examples due to underdetermined
and ill-conditioned systems of nonlinear SE equations. A possible explanation
for significantly lower MSEs for the IEEE 118-test case in these scenarios is
that it contains a greater variety of subgraphs for GNN training. To
investigate the GNN predictions further, we create a test set by excluding
five measurements connected to the two directly connected power system buses
from each test sample, resulting in the average test set MSE of $1.488\times
10^{-4}$. Fig. 4.12 shows the results for one test sample, where vertical
dashed lines correspond to the buses in the $1$-hop neighbourhood of the
excluded measurements. We can observe that the deviation from the ground truth
values manifests mainly in the vicinity of the excluded measurements, not
affecting the prediction accuracy in the rest of the power system.

020406080100$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Percentage of
Excluded MeasurementsAverage MSEIEEE 30IEEE 118 Figure 4.11: Average MSEs of
test sets created by randomly excluding measurements. 0.9711.031.061.09Voltage
Magnitude ${V}_{i}$PredictionsLabels051015202530-0.4-0.3-0.2-0.10Bus Index
$i$Voltage Angle ${\theta}_{i}$ Figure 4.12: GNN predictions and labels for
one test example, with all measurements connected to two neighbouring buses
removed. Dashed lines indicate the buses in the $1$-hop neighbourhood of the
excluded measurements.

###### Behaviour Under Malicious Data Injections

We examine the robustness of the proposed GNN model to malicious data
injection type of cyberattacks by randomly altering the values of five
neighbouring measurements in each test sample. We compare the proposed GNN
model’s predictions with the solutions of the GN method and the ground truth
values obtained using the GN method applied on the uncorrupted measurement
data. The GNN model demonstrated an order of magnitude better performance than
the GN method, with the average test set MSEs $1.281\times 10^{-4}$ and
$1.034\times 10^{-3}$, respectively. Fig. 4.13 depicts the comparison of the
state variable predictions under corrupted input data for one example from the
test set.

0.9611.041.081.12Voltage Magnitude ${V}_{i}$PredictionsGN based SEGround
truth051015202530-0.4-0.3-0.2-0.10Bus Index $i$Voltage Angle ${\theta}_{i}$
Figure 4.13: GNN predictions and GN based SE solutions for one test example
with corrupted input data.

####  4.4 Summary and future work

In this chapter, we introduced methods for linear and nonlinear SE based on
the GNN model specialized for operating on augmented power system factor
graphs. The method avoids the problems that traditional SE solvers face, such
as sensitivity to ill-conditioned cases, numerical instabilities and
convergence time depending on the state variable initialization. By testing
the GNN on power systems of various sizes, we observed the prediction accuracy
in the normal operating states of the power system and the sensitivity when
encountering false data injection cyberattacks and input data loss due to
communication irregularities.

The results showed that the proposed approach provides good results for large
power systems, and is an effective approximation method for traditional SE
solutions even with a relatively small number of training samples, indicating
its sample efficiency. The GNN model used in this approach is also fast and
maintains constant memory usage, regardless of the size of the power system.
More specifically, the computational complexity of the proposed GNN model
regarding the number of state variables is linear during the inference phase,
and it is possible to distribute the inference computation across multiple
processing units. Given these characteristics, the approach is worthy of
further consideration for real-world applications.

Since the proposed GNN model generates predictions even for underdetermined
systems of equations describing the SE problem, it could be applied to highly
unobservable distribution power systems. Another application of the proposed
model for nonlinear SE could be the fast and accurate initialization of the
traditional nonlinear SE solver, resulting in a hybrid approach that is both
model-based and data-driven.

While our work shows promising results, an important limitation is the
inability to quantify the uncertainty of the GNN predictions. However, we are
encouraged by ongoing research efforts to address this issue, as quantifying
uncertainty for GNN regression remains an open problem. For instance, [158]
proposes a Bayesian framework that uses assumed density filtering to quantify
aleatoric uncertainty and Monte Carlo dropout to capture epistemic uncertainty
in GNN predictions. In light of this, we believe implementing a similar
approach represents a promising future research direction.

##  Part II Dynamic Distribution Network Reconfiguration and Reinforcement
Learning

###  Chapter 5 Dynamic Distribution Network Reconfiguration

In this chapter, we introduce the foundations of static and dynamic
distribution network reconfiguration, stating their importance in distribution
management software. Furthermore, we provide a mathematical formulation of the
dynamic distribution network reconfiguration problem, which will be
transformed into the equivalent reinforcement learning formulation in Chapter
7111Chapters 5, 6, and 7 introduce a new nomenclature..

####  5.1 Distribution Network Reconfiguration

The electrical distribution network is the part of the electrical power system
which delivers electric power from the transmission system to individual
consumers. Traditional distribution networks consist of passive elements, with
power flows directed only from the network supply point to the customers. Due
to growing power demand, modern distribution networks exhibit changes such as
energy deregulation, increased installation of the distributed generation
coming mainly from renewable energy resources, and the deployment of
controllable loads. These changes bring numerous challenges to the operation
of modern distribution networks, such as bidirectional power flows, increased
system dynamics, transient instabilities, short-circuit conditions supplied by
multiple sources, overall operation inefficiency in terms of increased power
losses, decreased reliability, etc.

To overcome these challenges and improve response time in unforeseen
situations, modern distribution networks increase the level of network
automation by using remotely controlled equipment and employing domain-
specific software solutions such as DMS. DMS, which is usually tightly coupled
with the SCADA system, is used for distribution network monitoring, analysis,
optimization, and planning [159]. It includes functionalities like network
model management and topology processing, load flow, state estimation, short
circuit analysis, relay protection-related functionalities, Volt-VAR
optimization, distributed energy resources monitoring and control, etc.

Distribution network reconfiguration is one of the most important DMS
functionalities used for the optimization of distribution network operation.
In general form, DNR minimizes the objective function, which usually includes
power loss and voltage deviation, by changing the network topology using
manipulations on the switching devices [63]. During topology changes, it is
necessary to satisfy multiple constraints, such as not exceeding the bus
voltage, the apparent power of the branch, and the number of switching
manipulation limits. Additional constraints related to the network topology
are those that enforce the radiality of the network and ensure that all
customer buses are connected to the supply point. Some of the reasons for
enforcing network radiality are lower short circuit currents and simpler relay
protection setup than in the distribution networks containing loops. Fig. 5.1
represents an example of DNR on a simple 15-bus distribution network
consisting of three feeders, with black squares representing closed switches
and white squares representing open switches. The figure displays two
distribution network topologies, before and after the DNR, both of which are
radial and supply all customer buses.

In the scope of the DMS software, DNR functionality can have secondary goals
such as load balancing and Volt-VAR optimization [64], which are achieved by
expanding the DNR objective function. In emergency conditions, DNR can be used
to isolate the part of the distribution network where the fault occurred and
restore supply to the rest of the affected customers. In these cases, since
all the customers can not be supplied, DNR is used to minimize the number of
the disconnected customers, or the amount of energy not supplied. However, in
this thesis, we will not consider these emergency scenarios.

(a)

(b)

Figure 5.1: An example of distribution network before (subfigure a) and after
(subfigure b) the reconfiguration.

Static DNR is defined as DNR performed at a predefined time point, with fixed
load and generation values. In traditional distribution networks with low
network automation levels, where the customer load patterns change only
seasonally, static DNR was usually performed a few times a year. Due to
increasing load and generation dynamics, caused by controllable loads and
renewable energy resources, a need for more frequent and more flexible DNR
arises. Dynamic DNR optimizes the DNR objective function over the specified
time interval by finding the optimal time points when DNR should be performed.
DDNR is enabled by increased levels of network automation, as frequent changes
to the network topology cannot be performed quickly in a manual way. In a
usual formulation, the optimization interval is divided into time intervals,
in each of which a network topology change can be performed. However, since
switching manipulations have their own costs and can cause instabilities
during the topology changes, their number is usually also a subject of
minimization, in addition to being limited in the optimization problem
constraints. In other words, as the switching costs can be larger than the
network reconfiguration benefits, DDNR solves the trade-off problem between
finding the optimal topology in each time interval and performing less
frequent network topology changes to reduce the number of switching
manipulations. Operation planning using DDNR can be performed on a daily or
even hourly bases; therefore, it is necessary to develop fast algorithms that
produce high-quality DDNR solutions.

####  5.2 Mathematical Formulation of the DDNR Problem

The mathematical model of DDNR, formulated as a mixed-integer nonlinear
programming problem like in [72, 105] consists of a multi-objective function
and constraints. The multi-objective function defined as the total cost of
active energy loss and manipulation of switching devices is minimized, subject
to the following constraints: active and reactive power injection constraints,
bus voltage constraints, branch capacity constraints, switching operation
constraints, and a network radiality constraint. Decision variables of the
optimization problem consist of switch statuses in each of the time intervals,
which define the distribution network topology in that time interval.

#####  5.2.1 Objective Function

The DDNR problem is defined using the following multi-objective function:

$\min_{x_{sw}^{t}}\sum_{t=1}^{T}(C_{Loss}T_{int}^{t}P_{Loss}^{t}+C_{SWs}\,SW^{t}).$
(5.1)

In (5.1), $t\in 1\dots T$ denotes the index of a time interval, where $T$ is
the total number of time intervals and $T_{int}^{t}$ is the duration of a time
interval in hours. The decision variables of the optimization problem are the
status of the switches per time interval, where the status of the switch $sw$
in the time interval $t$ is defined as follows:

$x_{sw}^{t}=\begin{cases}1,&\text{if switch $sw$ is closed in time interval
$t$;}\\\ 0,&\text{if switch $sw$ is opened in time interval $t$.}\end{cases}$
(5.2)

$C_{Loss}$ is the cost of energy losses, in $\$$ per $kWh$, $P_{Loss}^{t}$
denotes the active power losses in $kW$, and $C_{SWs}$ is the cost in $\$$ of
the switching action for the $s^{th}$ switch. The total cost of switching
actions in the time interval $t$ is calculated using the number of switches
that had their status changed:

$SW^{t}=\sum_{sw=1}^{N_{sw}}y_{sw}^{t},$ (5.3)

where $sw\in 1\dots N_{sw}$ denotes the switch index and $N_{sw}$ is the total
number of switches in the distribution network. $y_{sw}^{t}$ indicates the
status change of a single switch $sw$ in time interval $t$:

$y_{sw}^{t}=\begin{cases}1,&\text{if the switch $sw$ changed its status in
time interval $t$;}\\\ 0,&\text{otherwise.}\end{cases}$ (5.4)

The active and reactive powers of the branches in the time interval $t$ are
defined as:

$\begin{gathered}P_{b}^{t}=g_{jk}(V_{j}^{t})^{2}-V_{j}^{t}V_{k}^{t}[g_{jk}cos(\theta_{j}^{t}-\theta_{k}^{t})+b_{jk}sin(\theta_{j}^{t}-\theta_{k}^{t})],\\\
Q_{b}^{t}=-b_{jk}(V_{j}^{t})^{2}+V_{j}^{t}V_{k}^{t}[b_{jk}cos(\theta_{j}^{t}-\theta_{k}^{t})-g_{jk}sin(\theta_{j}^{t}-\theta_{k}^{t})],\end{gathered}$
(5.5)

where $b\in 1\dots N_{b}$ denotes the branch index, $N_{b}$ is the total
number of branches, while $j$ and $k$ denote indices of buses at the ends of
the branch $b$. $V_{j}^{t}$ and $\theta_{j}^{t}$ denote the voltage magnitude
and the phase angle at the bus $j$ in the time interval $t$, while $g_{jk}$
and $b_{jk}$ represent elements of the nodal conductance and susceptance
matrices, respectively. Active power losses in the time interval $t$ are
defined as follows:

$P_{Loss}^{t}=\sum_{b=1}^{N_{b}}\alpha_{b}^{t}R_{b}\frac{(P_{b}^{t})^{2}+(Q_{b}^{t})^{2}}{(V_{j}^{t})^{2}},$
(5.6)

where $R_{b}$ denotes the resistance of the branch, and $\alpha_{b}^{t}$
combines information about the switch status and the existence of a switch on
a branch:

$\alpha_{b}^{t}=\begin{cases}1,&\text{if branch $b$ does not have a
switch;}\\\ x_{sw}^{t},&\text{if branch $b$ has the switch with index
$sw$.}\end{cases}$ (5.7)

#####  5.2.2 Constraints

The formulation of the DDNR problem considers multiple constraints listed
below.

  * $\bullet$

Active and reactive injection constraints, defined by the bus power balances
per time interval, as in the classical load flow model [160]:

$\begin{gathered}P_{j}^{t}=V_{j}^{t}\sum_{k=1}^{N_{n}}V_{k}^{t}[g_{jk}cos(\theta_{j}^{t}-\theta_{k}^{t})+b_{jk}sin(\theta_{j}^{t}-\theta_{k}^{t})],\\\
Q_{j}^{t}=V_{j}^{t}\sum_{k=1}^{N_{n}}V_{k}^{t}[g_{jk}sin(\theta_{j}^{t}-\theta_{k}^{t})-b_{jk}cos(\theta_{j}^{t}-\theta_{k}^{t})],\\\
j=1,\dots,N_{n},\end{gathered}$ (5.8)

where $N_{n}$ represents the number of buses in the network. Active and
reactive power injections are equal to the difference between the
corresponding load and generation in bus $j$, and they are assumed constant
during one time interval $t$.

  * $\bullet$

Slack bus constraints, which specify the voltage magnitude and the phase angle
in the root bus (i.e., the supply point) of the distribution network:

$\begin{gathered}V_{0}^{t}=V_{spec}^{t},\\\ \theta_{0}^{t}=0,\end{gathered}$
(5.9)

where $V_{spec}^{t}$ represents the specified slack bus voltage magnitude
value at time interval $t$. The slack bus provides an angular reference for
all other buses and balances the system’s active and reactive power [160].

  * $\bullet$

Bus voltage constraints:

$\begin{gathered}V_{j}^{min}\leq V_{j}^{t}\leq V_{j}^{max},\\\
j=1,\dots,N_{n},\end{gathered}$ (5.10)

where $V_{j}^{min}$ and $V_{j}^{max}$ denote the minimum and maximum voltage
magnitude allowed at the bus $j$.

  * $\bullet$

Branch capacity constraints:

$\begin{gathered}(P_{b}^{t})^{2}+(Q_{b}^{t})^{2}\leq(S_{b}^{max})^{2},\\\
b=1,\dots,N_{b},\end{gathered}$ (5.11)

where $S_{b}^{max}$ represents the maximum apparent power allowed in the
$b^{th}$ branch.

  * $\bullet$

Switching operation constraints:

$\begin{gathered}\sum_{t=1}^{T}|x_{sw}^{t}-x_{sw}^{t-1}|\leq N_{sw}^{max},\\\
sw=1,\dots,N_{sw},\end{gathered}$ (5.12)

where $N_{sw}^{max}$ represents the maximum number of allowed operations for
the $sw^{th}$ switch during the optimization time interval, which depends on
the type and the lifetime of the switch. $x_{sw}^{0}$ are the initial switch
statuses, and they do not belong to the decision variables. Constraints
containing absolute value operators can not be used directly in classical
mixed-integer algorithms, but require reformulation by introducing additional
variables. In the proposed RL approach, we will consider these constraints
directly by adding a penalty term to the reward function, resulting in a
simpler formulation of the problem.

  * $\bullet$

Network radiality constraint, which assures there are no loops in the
distribution network:

$\sum_{b=1}^{N_{b}}\alpha_{b}^{t}=N_{n}-1.$ (5.13)

The DDNR problem formulated in this way is NP-hard, since it has $2^{N_{sw}T}$
possible solutions, and it can not be solved in polynomial time. In the
forthcoming chapters, we will introduce RL algorithms which are trained to
search the solution space based on the agent-environment interaction concept,
and yield quality solutions with low computational effort during the
algorithm’s evaluation.

The proposed multi-objective formulation of the DDNR problem aims to minimize
the total cost of active energy loss and manipulation of switching devices.
However, in many practical applications the DDNR, apart from minimizing the
cost of active energy loss and manipulation of switching devices, the DDNR is
used for voltage deviation minimization [62, 63], load balancing, Volt/Var
optimization, supply restoration [64], the distribution network reliability
maximization [105], limiting the budget [93], etc. Extension of the DDNR
problem formulation assumes incorporating additional criteria into the
objective function. That is achieved by adding the corresponding terms to the
objective function while preserving the constraints (5.8) - (5.13).

###  Chapter 6 Reinforcement Learning

Reinforcement learning (RL), as a machine learning technique, deals with how
software agents learn to take actions in an environment through experience and
exploration, with the goal of finding the optimal strategy that maximizes the
long-term reward obtained [33]. In the RL framework, it is assumed that the
agent interacts with the generally stochastic environment in discrete time
steps. At the beginning of each time step, the agent observes the environment,
that is, it receives the state variables from the environment. Based on the
state variables, the agent takes an action and sends it to the environment.
The environment then changes its state due to the action received, as well as
due to its internal processes. After that, the environment sends the immediate
reward signal for that time step and the state variables for the next time
step to the agent. The goal of an RL algorithm is to find the (close to)
optimal policy, i.e., the action selection that maximizes the long-term
reward, while receiving feedback about its immediate performance. This chapter
presents the theoretical background of RL in Sections 6.1 and 6.2 and the deep
Q-learning algorithm in Section 6.3, which is applied to DDNR in Chapter 7.

####  6.1 Finite Markov Decision Processes

Finite Markov decision processes are discrete-time stochastic control
processes that model decision-making in situations in which the outcome is
partially stochastic and partially under the control of the decision-maker.
The result of the solved MDP is the optimal sequence of actions, that is, the
optimal policy. The RL problem can be mapped onto the MDP, which is defined as
the following tuple:

  * $\bullet$

$\mathcal{S}$ \- finite set of states,

  * $\bullet$

$\mathcal{A}$ \- finite set of actions,

  * $\bullet$

$\mathcal{R}$ \- finite set of immediate rewards,

  * $\bullet$

$p\left(s^{\prime}\,\lvert\,s,a\right)=\Pr\left(S_{t}=s^{\prime}\,\lvert\,S_{t-1}=s,A_{t-1}=a\right)$
\- transition probability function.

Random variables $S_{t}$, $S_{t-1}$, $R_{t}$, $A_{t-1}$ represent the new
state, the previous state, received immediate reward and the action being
taken, respectively, while $s^{\prime},s\in\mathcal{S}$, $r\in\mathcal{R}$,
$a\in\mathcal{A}$ denote the values of these random variables. The transition
function represents the probability of being in the state $s^{\prime}$ on the
condition of being in the state $s$ previously and executing the action $a$.
The agent and the environment interact in discrete time steps, as shown in
Fig. 6.1. In each time step the agent observes the state $s$, makes an action
$a$, upon which the environment changes, and sends the feedback to the agent
in the form of reward $r$, and the next state $s^{\prime}$.

Figure 6.1: The agent-environment interaction process.

Additionally, a problem defined as MDP must satisfy the Markov property,
meaning that the next state $s^{\prime}$ depends only on the current state $s$
and the action $a$, and is independent of all previous states and actions:

$\Pr\left(S_{t}\,\lvert\,S_{0},A_{0},S_{1},A_{1},\dots,S_{t-1},A_{t-1}\right)=\Pr\left(S_{t}\,\lvert\,S_{t-1},A_{t-1}\right)$
(6.1)

Problems that do not satisfy the Markov property can be solved exactly using
traditional MDP methods by expanding the state space with data from previous
states or approximately using the RL methods.

Return in time step $t$, which represents the long-term reward starting from
that time step, is the subject of optimization:

$G_{t}=R_{t+1}+\gamma
R_{t+2}+\gamma^{2}R_{t+3}+\dots=\sum_{i=0}^{\infty}\gamma^{i}r_{t+i+1},$ (6.2)

where $\gamma\in[0,1]$ represents the discount factor. The case of $\gamma=0$
corresponds to the greedy maximization of the immediate reward, while using
$\gamma=1$ implies equal weight on all rewards in the optimization horizon.
Due to convergence problems in the case of long optimization horizons, the
most often used value of the discount factor is $\gamma\in[0.9,0.99]$.

Policy $\pi$ models the action selection in various states. When optimizing
the policy, the long-term reward is optimized. The policy is in the general
case stochastic, i.e., it maps the probability distribution of actions to
states $\pi(a|s)=\Pr(A_{t}=a|S_{t}=s)$, but it can be also defined in a
deterministic way $a=\pi(s)$. The quality of a policy $\pi$ in a state $s$ is
usually expressed as an expectation of the discounted long-term reward $G_{t}$
given the state $s$, enabling the way to compare different policies and
optimize them. This quantity is called the state value function, and it
assigns higher values to the more desirable states in terms of the long-term
reward, if following the policy $\pi$111Superscript $\pi$ denotes that the
quantity is calculated with the assumption that the agent takes actions
according to the policy $\pi$.:

$v^{\pi}(s)=\mathop{\mathbb{E}}_{\pi}\,[G_{t}\,\lvert\,S_{t}=s]=\mathop{\mathbb{E}}_{\pi}\,[\sum_{i=0}^{\infty}\gamma^{i}R_{t+i+1}\,\lvert\,S_{t}=s].$
(6.3)

Similarly, the action value function (Q-function) is a mapping of states-
action pairs to the real numbers, where the value of the state-action pair
represents the expected discounted long-term reward starting from that state,
taking that action, and following a concrete policy $\pi$ afterward:

$q^{\pi}(s,a)=\mathop{\mathbb{E}}_{\pi}\,[G_{t}\,\lvert\,S_{t}=s,A_{t}=a]=\mathop{\mathbb{E}}_{\pi}\,[\sum_{i=0}^{\infty}\gamma^{i}R_{t+i+1}\,\lvert\,S_{t}=s,A_{t}=a].$
(6.4)

The state value function has a more compact representation, while the
Q-function provides a simpler way to determine the action that leads to the
most desirable future state. The state value function and the Q-function can
be recursively expressed using the Bellman equations:

$\begin{split}v^{\pi}(s)&=\mathop{\mathbb{E}}_{\pi}\,[R_{t+1}+\gamma
v^{\pi}(S_{t+1})\,\lvert\,S_{t}=s]\\\
&=\sum_{a\in\mathcal{A}}\pi\left(a\,\lvert\,s\right)(r+\gamma\sum_{s^{\prime}\in\mathcal{S}}p(s^{\prime}|s,a)v^{\pi}(s^{\prime})),\end{split}$
(6.5) $\begin{split}q^{\pi}(s,a)&=\mathop{\mathbb{E}}_{\pi}\,[R_{t+1}+\gamma
q^{\pi}(S_{t+1},A_{t+1})\,\lvert\,S_{t}=s]\\\
&=r+\gamma\sum_{s^{\prime}\in\mathcal{S}}p\left(s^{\prime}\,\lvert\,s,a\right)\sum_{a^{\prime}\in\mathcal{A}}\pi\left(a^{\prime}\,\lvert\,s^{\prime}\right)q^{\pi}(s^{\prime},a^{\prime}).\end{split}$
(6.6)

Finding the optimal state value function or the optimal Q-function results in
finding the optimal policy. The optimal policy can be generated from e.g. the
optimal Q-function by selecting the action with the largest Q-function value
in each state. If the problem is formulated as an MDP, then at least one
optimal solution exists, and an iterative procedure based on dynamic
programming and the Bellman equations that converges to one of those solutions
can be established. Some of the commonly used algorithms are the value
iteration algorithm [161], and the policy iteration algorithm [162].

To find the exact solution to the MDP, it has to be fully defined, i.e., all
transition probabilities and immediate rewards have to be known, and it has to
satisfy the Markov property. Additionally, for large state and action spaces,
solving MDPs exactly could be computationally infeasible. Partial
observability of the environment state additionally increases the
computational time of the traditional algorithms that solve MDPs [163]. The
main idea of RL is to overcome these problems by learning (close to) optimal
policies based on the history of interactions of the agent with the
environment.

####  6.2 Q-Learning

Our work considers model-free off-policy RL algorithms, where the optimal
policy is learned directly from the accumulated experience, i.e., the history
of the interaction process between the agent and the environment. On the other
hand, in the model-based RL, the MDP, or the transition probabilities and
immediate rewards for all state-action-next state triplets are learned from
the accumulated experience and solved to obtain the policy.

Q-learning is a basic model-free RL algorithm, where the values of the
Q-function for each state-action pair are stored in the Q-table and updated
during algorithm training [164]. During the evaluation of the algorithm, for
each state the agent receives from the environment, the action with the
greatest Q-function value in the table is selected. During the algorithm
training, actions are selected randomly with the probability $\epsilon$, and
the actions with the largest Q-function values for the corresponding states
are selected with the probability $1-\epsilon$, where $\epsilon\in[0,1]$ is
the exploration hyperparameter. This way the agent searches the state-action
space and avoids the local optima problem. Algorithm training is performed by
repeating the predefined number of episodes, which consist of time steps. One
time step contains information about one interaction of the agent with the
environment: the current state, the action selected by the agent, the received
reward, and the next state. The length of episodes, i.e., the number of time
steps in them, is generally variable.

Upon one interaction of the agent with the environment, the values in the
Q-table are updated using the following rule, obtained using the idea from the
Bellman equations:

$q(s,a):=(1-\alpha)q(s,a)+\alpha[r+\gamma\;\underset{a^{\prime}}{\max}\;q(s^{\prime},a^{\prime})],$
(6.7)

where $\alpha$ is the learning rate hyperparameter. As well as MDPs, the
Q-learning algorithm assumes discrete state and action spaces and for large
state and action spaces it may be infeasible to learn the Q-function value for
all state-action pairs.

####  6.3 Deep Q-learning

The deep Q-Learning algorithm is one of the basic DRL algorithms, that
utilizes the advances in the deep learning field to improve the traditional RL
algorithms. The idea of the algorithm is to use a deep neural network, also
called the deep Q-network (DQN), as an approximator of the Q-function [165].
Inputs to the DQN are state variables, while output neurons provide the
approximation of Q-function values for each of the actions and for the input
state. Therefore, the state variables can be continuous, which makes the
learning feasible for large continuous state spaces which would have to be
discretized when using the Q-learning algorithm. In deep Q-learning, the
action space must be discrete and finite, since the number of output neurons
is limited. An example of a DQN is shown in Fig. 6.2, where $Q(s,a_{i})$
denotes the DQN output when the agent takes action $a_{i},i=1,\dots,z$, while
being in the state $s$, where $z=|\mathcal{A}|$ denotes the number of possible
actions.

Figure 6.2: An example of a deep Q-network.

Deep Q-learning algorithm introduces target DQN $Q_{target}$ which has the
same model architecture as DQN and in which the parameters of DQN are copied
at the predefined period during the training process. Target DQN is used for
determining the labels for the DQN training, as defined in (6.8), which
significantly improves the training process stability [165]. It reduces the
oscillations of the training by fixing the label generation process during
multiple training steps, allowing the DQN network weights to be updated
steadily. Target networks are a widely used technique in RL, and recent RL
algorithms propose more advanced variations, such as continuously updating the
time-delayed target network [166].

On-policy RL algorithms [33] are trained by updating the same policy using
which the sequence of actions is generated, which results in unstable and
sample inefficient training due to correlation between the actions in the
sequence. Deep Q-learning is an example of an off-policy RL algorithm, which
stores the history of the agent’s interaction with the environment in the
experience replay memory [167], and samples data randomly from it to perform
the DQN training in a supervised learning manner. At each time step, tuple
$(s,a,r,s^{\prime})$ is stored in the replay memory, from which i.i.d. mini-
batch data for DQN training is sampled.

Labels for DQN training are calculated in the following way, using the idea
from the Bellman equations, similarly to the Q-learning algorithm:

$Q_{label}(s,a)=r+\gamma\;\underset{a^{\prime}}{\max}\;Q_{target}(s^{\prime},a^{\prime}|\theta^{Q_{target}}).$
(6.8)

$\theta^{Q}$ and $\theta^{Q_{target}}$ denote the parameters (weights and
biases) of the DQN and the target DQN. DQN is trained using the mini-batch
gradient descent algorithm [168], which minimizes the squared error loss
function that expresses the distance between the labels and the DQN output
during training:

$L(\theta^{Q})=\frac{1}{N_{mb}}\sum_{i=1}^{N_{mb}}(Q_{label}(s_{i},a_{i})-Q(s_{i},a_{i}\,\lvert\,\theta^{Q}))^{2}.$
(6.9)

A trained DQN is evaluated by forwarding the input state $s$ through the
network layers, obtaining the Q-function approximates $Q(s,a_{i})$ for all
actions $a_{i},i=1,\dots,|\mathcal{A}|$, and selecting the action with the
largest Q-function value.

###  Chapter 7 Reinforcement Learning based Dynamic Distribution Network
Reconfiguration

This chapter describes the way DDNR is expressed as an RL problem, how the
objective function and constraints are considered, and the training and
evaluation algorithms of the proposed DQN-based method. Finally, we evaluate
the performance of the proposed approach on three distribution networks:
15-bus test benchmark, real-life large-scale distribution network, and the
IEEE 33-bus network.

####  7.1 Modelling Dynamic Distribution Network Reconfiguration as a Markov
Decision Process

The information flow between the DDNR agent and the environment during their
interaction is presented in Fig. 7.1. Episodes consist of $T$ time steps,
where each episode corresponds to a separate instance of the DDNR applied to
the one-day interval, and each RL agent’s time step corresponds to the
beginning of one time interval in DDNR problem formulation (5.1). At each time
step, active and reactive power consumption data for the next hour are loaded.
Then, the power flow calculation is executed to create the state variables,
which contain the time interval index $t$ and the apparent powers of all
switches in the network $S_{sw}^{t}$, $sw\in 1,\dots,N_{sw}$, as shown in
(7.1):

$s_{t}=(t,S_{1}^{t},\dots,S_{N_{sw}}^{t}).$ (7.1)

This choice of state variables, motivated by the fact that switch statuses,
and hence the network topology can be reconstructed using the apparent powers
of the switches, reduces the state space dimensionality and DQN size. The
current network configuration and power flow results are compressed into a
single set of variables, from which the agent can make its own representation
of the environment, and use it as an input in the decision-making process.

Agent: DDNR controller Environment: consumption data + power flow calculation
Action: switch combination Reward: \- DDNR objective \- Voltage constraints \-
Branch power constraints State: \- Time step index \- Apparent powers of
switches Figure 7.1: The agent-environment interaction process for DDNR.

The action space contains all the switch combinations that lead the network in
a feasible radial configuration, in which all the buses are energized. These
switch combinations are enumerated uniquely so that one output neuron
corresponds to one feasible radial combination. This action space definition
implies that the radiality network constraint stated in (5.13) is always
satisfied, which accelerates the learning process.

The reward value for each time interval $t$ is equal to the negative sum of
the following terms:

  * $\bullet$

Price of active energy losses $C_{Loss}T_{int}^{t}P_{Loss}^{t}$, used to model
the first term in the DDNR objective function stated in (5.1);

  * $\bullet$

Price of the switching manipulations needed to conduct the network from the
previous configuration to the current one $C_{SWs}\,SW^{t}$, used to model the
second term in the DDNR objective function stated in (5.1);

  * $\bullet$

Predefined penalty value $C_{U}$ if the bus voltage constraint is violated in
any of the buses, used to model the bus voltage constraints in (5.10);

  * $\bullet$

Predefined penalty value $C_{S}$ if the branch capacity constraint is violated
for any of the branches, used to model the branch capacity constraints in
(5.11).

As an alternative to adding the predefined penalty to the reward function if
the number of switch manipulations exceeds the predefined limit for any of the
switches, we propose the following way to consider the switching operation
constraints in (5.12):

  * $\bullet$

The subset of available actions is created at the beginning of each episode
and is initially equal to the action set (i.e., all the switch combinations
that lead the network into a feasible radial configuration, in which all the
buses are energized);

  * $\bullet$

During the episode, the number of operations of each switch in that episode is
updated. Prior to action selection, actions that would violate the switching
operation constraints if selected are removed from the subset of available
actions;

  * $\bullet$

In each time interval, the action with the largest Q-value is selected from
the subset of available actions, instead of from the action set.

This approach improves the training efficiency compared to the approach that
would penalize the exceeded number of switch operations. Since the actions
that violate the switching operation constraints cannot be selected in the
first place, the computational effort needed to learn Q-function values for
those state-action pairs is eliminated. It is convenient to consider the
switching operation constraints using this approach since only the selected
action (switch combination) needs to be known, to conclude if the constraint
is violated. The subset of available actions can be determined without
executing the actions. Bus voltage and branch capacity constraints cannot be
modelled using this approach generally, since the action must be executed and
feedback from the environment is required for any conclusions about the
constraint violations. Active and reactive injection constraints in (5.8), as
well as slack bus constraints in (5.9) are satisfied by the design of the
power flow calculation [160] and, therefore, are not considered in the reward
function. Values of the state variables and rewards are normalized, for the
purpose of improving the numerical stability of neural network training.

####  7.2 Training and Evaluation Algorithms

During the algorithm training $N$ episodes are repeated, with each episode
consisting of a predefined number of time steps $T$, where in each time step
the interaction between the agent and the environment takes place, as
described in 7.1. The variety of training scenarios is created by randomly
sampling the daily load curves from some predefined distribution. The way the
exploration hyperparameter $\epsilon$ is updated during the training can also
be tuned. In this work, we used the linear decrease of $\epsilon$ until the
$0.8N^{\textrm{th}}$ episode, and the constant value afterward. A detailed
representation of the agent training procedure is displayed in Algorithm 1.

Initialize the deep Q-network’s $Q(s\mid\theta^{Q})$ parameters randomly

Initialize the target deep Q-network’s $Q^{\prime}(s\mid\theta^{Q_{target}})$
parameters using the original network’s parameters

Initialize the experience replay buffer

for  _$episode=1,2,\ldots,N$_ do

Sample daily load curves randomly

Initialize $\epsilon$

Initialize the subset of available actions to the action set

Run the initial power flow calculation

Send the initial state $s_{1}$ to the agent

for  _$t=1,2,\ldots,T$_ do

$rand$ = random number between $0$ and $1$

if  _$rand >\epsilon$_ then

Based on the current state $s_{t}$ select the action $a_{t}$ with the largest
Q-value in the subset of available actions

else

Select random action $a_{t}$ from the subset of available actions

end if

Update the subset of available actions

Update the network configuration according to $a_{t}$

Run the power flow calculation

Collect the immediate reward $r_{t}$ and the next state $s_{t+1}$ data

Store tuple $(s_{t},a_{t},r_{t},s_{t+1})$ in the experience replay buffer

Sample the mini-batch of tuples from the experience replay buffer

Create labels for deep Q-network training using (6.8)

Update deep Q-network parameters by minimizing the loss function given in
(6.9)

Set loads for the next time interval for each bus

end for

if  _update target network period_ then

$\theta^{Q_{target}}=\theta^{Q}$

Update $\epsilon$

end for

Algorithm 1 Deep Q-network training

Once trained, the deep Q-network model can be evaluated multiple times by
loading the saved trainable parameters, as shown in Algorithm 2. Note that
during the algorithm evaluation target deep Q-network, experience replay
buffer, neural network parameter updates, and the random exploration strategy
are not used, reducing the memory storage and computational requirements per
episode, when compared to the training algorithm. The trained algorithm’s
evaluation reduces to $T+1$ power flow calculations and $T$ neural network
evaluations, which are almost instantaneous, resulting in a computationally
efficient control procedure, which can be used either standalone or as a part
of the more complex power systems’ application.

Load previously saved parameters of the trained deep Q-network
$Q(s\mid\theta^{Q})$

Generate daily load curves for the evaluation example

Initialize the subset of available actions to the action set

Run the initial power flow calculation

Send the initial state $s_{1}$ to the agent

for  _$t=1,2,\ldots,T$_ do

Based on the current state $s_{t}$ select the action $a_{t}$ with the largest
Q-value in the subset of available actions

Update the subset of available actions

Update the network configuration according to $a_{t}$

Run the power flow calculation

Retrieve the next state $s_{t+1}$

Set loads for the next time interval for each bus

end for

Algorithm 2 Deep Q-network evaluation

####  7.3 Numerical Results

In this section, results and discussion are presented for benchmark test
examples – 7.3.1, real-life large-scale distribution network – 7.3.2, and the
IEEE 33-bus radial system – 7.3.3, along with the choice of RL and deep
learning hyperparameters. The proposed algorithms were implemented in Python,
deep neural networks were modelled and trained using the PyTorch deep learning
framework, and power system-related modelling and calculations were completed
using OpenDSSDirect, the Python interface to OpenDSS distribution system
simulation software [169]. The algorithms were executed on a 64-bit Windows 10
with the following hardware configuration: AMD A8-6410 APU with AMD Radeon R5
Graphics 2.00 GHz, 4 cores, and 8 GB of RAM.

#####  7.3.1 Benchmark Test Examples

Fig. 7.2 illustrates a $15$-bus test benchmark where a slack bus is the bus
with the marker $0$ and the other $14$ buses are of the PQ type. Loads are
defined by the chronological daily diagrams, which are sampled uniformly from
the intervals defined by dashed lines in Fig. 7.3. The length of all branches
is $4.5$km. All branches are balanced with the direct sequence impedance
$r+jx=(0.224+j0.109)\Omega/$km. All branches have switching devices, and the
number of switch manipulations is not constrained.

Figure 7.2: Single-line diagram for $15$-bus test benchmark.
00.20.40.60.811.21.4Load [MW]Feeder 1 load Limits00.20.40.60.811.21.4Load
[MW]Feeder 2 loadLimits024681012141618202200.20.40.60.811.21.4Hour [h]Load
[MW]Feeder 3 loadLimits Figure 7.3: Daily load profiles for three feeders.
Full lines represent average load values, and dashed lines represent limits
between which training set loads are sampled.

Penalty values used to model the bus voltage constraints in (5.10) and the
branch capacity constraints in (5.11) are: $C_{U}=C_{S}=10$. The costs of
energy losses and switching operations in the time interval are [100, 170]:

  * $\bullet$

Cost of energy losses ($C_{Loss}$): $6.5625 cents/kWh;

  * $\bullet$

Cost of switching manipulations ($C_{sw}$): $1 per manipulation.

DQN used for this test example consists of the input layer, four hidden
layers, and the output layer. The input layer has $15$ neurons, one for the
time step index variable, and $14$ for the apparent powers of each switch.
Each hidden layer has $256$ neurons, and the output layer consists of $186$
neurons, one for each switch combination that leads to a feasible radial
configuration. Rectified Linear Unit (ReLU) activation function is applied to
each of the hidden layers and the output layer. The neural network is trained
using the Adam optimizer, with a learning rate of $10^{-5}$ and a mini-batch
size of $128$. We experimented with adding batch normalization on several
hidden layers, but it neither improved nor deteriorated the training results.

The number of training episodes we used is $60000$. During the training, loads
for each hour were uniformly sampled from the intervals defined by dashed
curves in Fig. 7.3. The target DQN update frequency is $10$ episodes and the
mini-batches for DQN training were sampled from experience replay memory,
which has the capacity of $10^{6}$ samples. The initial value of the
exploration parameter $\epsilon$ is $1$ and decreases linearly to the episode
index until it reaches the value $0.1$ in $48000^{\textrm{th}}$ episode. The
value of the discount factor $\gamma$ is $0.99$, as in [165].

Fig. 7.4. presents the average value of the DQN loss function, defined in
(6.9), over the episodes. As the training advances, the DQN loss decreases,
which implies that the Q-function is being approximated successfully.
Additionally, the testing performance of the RL algorithm with increasing
training episodes is presented in Fig. 7.4 by displaying the amount of the
normalized received reward per episode and its moving average. These rewards
were obtained by executing the episodes with the exploration parameter
$\epsilon$ equal to zero, after each training episode. The training
demonstrated asymptotic convergence within $20000$ episodes.

Figure 7.4: Average DQN loss per episode (top) and total reward per episode
along with its moving average (bottom).

Table 7.1 presents load, losses, and switch status changes for the proposed
approach in $24$-hour optimization horizon, as well as the comparison of the
same results with the state-of-the-art method from [105]. The method from
[105] minimizes the costs of energy losses, switching manipulations, and
outages. For this comparison, the method from [105] is executed with the cost
of outages set to zero. A graphical representation of switch statuses is
presented in Fig. 7.5. The results are additionally compared with the method
from [75], as is presented in Table 7.2. By comparing the results of the
proposed algorithm and the method from [105], it can be concluded that the
proposed switching actions are not the same, but the total costs only differ
slightly. The execution time of the proposed algorithm is $0.148$s, which is
two orders of magnitude smaller than the execution time of the method from
[105].

Table 7.1: Total load, active power losses and switch status changes in the 24-hour time optimization period for the $15$-bus test benchmark (O–open; C–close). |  | Proposed approach | Method from [105]
---|---|---|---
Hour | Load [kW] | Losses [kW] |  | Switch status
---
changes
Losses [kW] |  | Switch status
---
changes
0 | 4554.0 | 114.06 | 4(O), 14(C) | 131.73 | No changes
1 | 4305.0 | 101.79 | No changes | 115.32 | No changes
2 | 3876.0 | 82.49 | No changes | 89.97 | No changes
3 | 3326.0 | 61.20 | No changes | 62.95 | No changes
4 | 3205.0 | 57.03 | No changes | 57.81 | No changes
5 | 3693.0 | 75.01 | No changes | 80.35 | No changes
6 | 7542.0 | 316.63 |  | 10(O), 14(O),
---
4(C), 13(C)
316.63 | 10(O), 13(C)
7 | 7909.0 | 348.71 | No changes | 348.71 | No changes
8 | 8279.0 | 384.08 | No changes | 384.08 | No changes
9 | 6067.0 | 213.94 | 10(C), 13(O) | 220.01 | No changes
10 | 9369.0 | 465.82 | No changes | 465.82 | 10(C), 13(O)
11 | 9054.0 | 440.61 | No changes | 440.61 | No changes
12 | 8823.0 | 423.20 | No changes | 423.20 | No changes
13 | 8823.0 | 423.20 | No changes | 423.20 | No changes
14 | 7670.0 | 375.02 | No changes | 324.18 | 10(O), 13(C)
15 | 9801.0 | 493.87 | No changes | 493.87 | 10(C), 13(O)
16 | 7440.0 | 278.42 | 4(O), 14(C) | 278.42 | 4(O), 14(C)
17 | 7049.0 | 252.96 | No changes | 252.96 | No changes
18 | 7317.0 | 276.55 | No changes | 276.55 | No changes
19 | 6643.0 | 227.18 | No changes | 227.18 | No changes
20 | 5651.0 | 178.27 | No changes | 178.27 | No changes
21 | 5163.0 | 147.68 | No changes | 147.68 | No changes
22 | 4793.0 | 126.72 | No changes | 126.72 | No changes
23 | 4554.0 | 114.06 | No changes | 114.06 | No changes
Table 7.2: Total losses, number of switch status changes, and total cost in the $24$-hour time optimization period. | Proposed approach | Method from [105] | Method from [75]
---|---|---|---
Total losses [kW] | 5978.48 | 5980.28 | 5967.56
| Number of switch
---
status changes
10 | 10 | 22
Total cost [$] | 402.3 | 402.4 | 413.6
01Other01Switch 1401Switch statusSwitch 1301Switch 1201Switch
100510152001Switch 4 Figure 7.5: Switch status changes during the 24-hour
period.

Based on Table 7.2, it can be concluded that the proposed approach and the
method from [105] provide better solutions than the method from [75].
Solutions from the proposed approach and the method from [105] provide greater
cost savings, and they achieve it with a lower number of switch manipulations.

Total losses in the time period of $24$ hours for $15$-bus test benchmark
without DNR are $6477.81$kW, with the total cost of $\$425.1$. With the
proposed DDNR, losses are reduced to $5978.48$kW and the total cost is reduced
to $\$402.3$. Fig. 7.6 compares loss reduction per hour for the proposed
approach and for the method from [105]. In Fig. 7.6 it can be seen that there
is no reduction of losses in the interval from $10$ to $16$ hours. In Table
7.1 it can be seen that the network configuration from 10 to 16 hours is the
same as the basic network configuration, as in Fig. 7.2 (in the first hour
switch 4 was opened and switch 14 was closed; in the seventh hour switch 4 was
closed and switch 14 was opened and switch 10 was opened and switch 13 was
closed; in the tenth hour switch 10 was closed and switch 13 was opened).
Consequently, there is no reduction of losses from 10 to 16 hours. For the
method from [105] in the $10^{\textrm{th}}$ hour, there is a negative loss
reduction (see Table 7.1). From 17 to $24$ hours, loss reduction is the same
for both methods, because network configuration in that interval is the same,
too.

012345678910111213141516171819202122230102030405060Hour [h]Loss reduction
[kW]Proposed approachMethod from [105] Figure 7.6: Loss reduction using DDNR.

Fig. 7.7 shows the voltage profile before and after reconfiguration. It is
noted that the voltage profile is improved after applying the configuration
found by the proposed algorithm. Bus loads are taken from the daily load
profile presented in Fig. 7.3 at $17^{\textrm{th}}$ hour.

Next, we present the results of the proposed approach when the number of
switch manipulations is limited to two. A graphical representation of the
switch status changes in this case is given in Fig. 7.8. The cost of losses is
equal to $403.5, while the switching manipulation cost is $12, summing up to
the total cost of $415.5, which means that this solution is more expensive
than the solution presented in Table 7.1.

024681012140.920.930.940.950.960.970.980.991Bus indexBus voltages [p.u.]Base
caseAfter reconfiguration Figure 7.7: Voltage profile for 15-bus test
benchmark. 01Other01Switch 1401Switch statusSwitch 1301Switch 1201Switch
1101Switch 1001Switch 70510152001Switch 4 Figure 7.8: Switch status changes
during the 24-hour period when the maximal number of switch manipulations is
two.

#####  7.3.2 Real-Life Large-Scale Distribution Network

To assess the scalability of the proposed approach, we evaluate its results
using a real-life large-scale radial distribution network, which consists of
four feeders, 1015 branches, and 1008 loads. The distribution network is
equipped with 31 remotely controlled switches, where 24 of them are normally
closed and seven are normally open. This large number of remotely controlled
switches is used to demonstrate how well does the developed algorithm scale
with the increase in the number of decision variables. DQN used for this test
example consists of the input layer, four hidden layers, and the output layer.
The input layer has 32 neurons, one for the time step index variable and 31
for the apparent powers of each switch. Each hidden layer has 1024 neurons,
and the output layer consists of 3567 neurons, one for each switch combination
that leads to a feasible radial configuration. The algorithm was trained on
100000 episodes. The initial value of the exploration parameter $\epsilon$ is
1 and decreases linearly to the episode index, until it reaches the value 0.1
in $80000^{\mathrm{th}}$ episode. The rest of the parameters are the same as
in the case of the 15-bus test benchmark network. Costs of energy losses and
switching manipulations are set to the same values as in Section 7.3.1.

Table 7.3 presents the resulting losses and switch manipulations of the
proposed approach for the 24-hour optimization horizon. The execution time of
the proposed algorithm is 3.635s, which is two orders of magnitude smaller
than the execution time of the method from [105] when used for networks with
31 remotely controlled switches. The training of the RL algorithm demonstrated
convergence after 30000 episodes, as evidenced by the analysis of its test
time performance after each training episode.

Table 7.3: Active power losses and switch status changes in the 24-hour time optimization period for the large-scale radial distribution network (O–open; C–close). Hour | Losses [kW] | Switch status changes
---|---|---
0 | 182.30 | 238(O), 900 (O), 994(O), 1011 (C), 1013(C), 1014 (C)
1 | 160.98 | No changes
2 | 127.19 | No changes
3 | 92.46 | No changes
4 | 86.08 | No changes
5 | 196.43 | No changes
6 | 742.05 | 695(O), 742(O), 900 (C), 1015 (C)
7 | 783.53 | No changes
8 | 826.54 | No changes
9 | 368.21 | 947 (O), 1014 (O), 695 (C), 994 (C)
10 | 1008.27 | No changes
11 | 965.66 | No changes
12 | 932.39 | No changes
13 | 891.17 | 1011 (O), 1013 (O), 238 (C), 947 (C)
14 | 746.81 | No changes
15 | 1044.24 | No changes
16 | 786.83 | 947 (O), 1013 (C)
17 | 553.33 | No changes
18 | 744.87 | No changes
19 | 645.53 | 1015 (O), 1014 (C)
20 | 509.13 | 191 (O), 1010 (C)
21 | 455.83 | No changes
22 | 313.87 | No changes
23 | 178.78 | No changes

The total active power losses are equal to 13342.49 kW, with a cost of $875.6.
The switching manipulation cost is equal to $24, which sums up to the total
cost of $899.6. The results demonstrate that the proposed DDNR is applicable
to real-life large-scale distribution networks.

#####  7.3.3 IEEE 33-bus Radial System

In this section, a numerical test is performed by using the IEEE 33‐bus radial
system [171], displayed in Fig. 7.9. We used the same branch resistances,
reactances, and peak loads as in [171]. The loads between buses 2 and 18 are
scaled with the blue full line in Fig. 7.3, loads between buses 19 and 25 are
scaled with the black full line in Fig. 7.3, and the loads between buses 26
and 33 are scaled with the red full line in Fig. 7.3. This network is modified
to match the developed DDNR. The IEEE33 test system is equipped with 20
remotely controlled switches, where 15 of them are normally closed and 5 are
normally open, as shown in Fig. 7.9. Remotely controlled switches are marked
with ”s” and the unique index. If the remotely controlled switch is placed on
the full line, it is closed, otherwise, it is open. Five lines were added to
the original scheme and are presented with dashed lines in Fig. 7.9. Costs of
energy losses and switching manipulations are the same as in Section 7.3.1.

Figure 7.9: IEEE 33‐bus radial system.

The algorithm was trained with the same training hyperparameters as in the
case of the 15-bus test benchmark network. Table 7.4 presents the losses and
switch manipulations for the proposed approach in the 24-hour time
optimization period. The execution time of the algorithm is 0.272s.

Table 7.4: Active power losses and switch status changes in the 24-hour time optimization period for the IEEE 33-bus radial system (O–open; C–close). Hour | Losses [kW] | Switch status changes
---|---|---
0 | 34.87 | 6 (O), 35 (C)
1 | 30.41 | No changes
2 | 23.50 | No changes
3 | 16.03 | No changes
4 | 14.60 | No changes
5 | 20.85 | No changes
6 | 56.49 | No changes
7 | 62.11 | 8 (O), 33 (C)
8 | 70.44 | No changes
9 | 50.05 | No changes
10 | 95.42 | 14 (O), 34 (C)
11 | 87.05 | No changes
12 | 81.22 | No changes
13 | 81.22 | No changes
14 | 57.06 | No changes
15 | 108.47 | No changes
16 | 86.87 | 27 (O), 37 (C)
17 | 75.20 | No changes
18 | 81.43 | No changes
19 | 65.80 | No changes
20 | 47.51 | No changes
21 | 38.90 | No changes
22 | 33.00 | No changes
23 | 29.45 | No changes

Total active power losses are equal to 1347.95 kW, with a cost of $88.5. The
cost of switching manipulation is equal to $8, which sums up to the total cost
of $96.5. The RL algorithm was found to converge after 15000 training
episodes, as determined by an analysis of the test time performance with
increasing episodes. Fig. 7.10 shows the voltage profile before and after the
reconfiguration. It is noted that in most of the buses, the voltage profile is
enhanced after applying the configuration found by the proposed algorithm. Bus
loads are taken from the daily load profile presented in Fig. 7.3 at
$17^{\mathrm{th}}$ hour.

0510152025300.910.920.930.940.950.960.970.980.991Bus indexBus voltages
[p.u.]Base caseAfter reconfiguration Figure 7.10: Voltage profile for IEEE
33-bus radial system.

####  7.4 Summary and future work

In this chapter, we present a multi-objective formulation of the DDNR aimed at
reducing the total cost of energy losses and switching operations. Our
solution is based on DRL, and we demonstrate that the proposed definition of
state variables, which demands lower observability, enables successful
treatment of the DDNR problem within the RL framework. Additionally, the
proposed computationally efficient method of addressing switching operation
constraints by creating a subset of available actions and updating it during
the episode was proven to be effective for this problem. The ideas from the
suggested way of modelling the DDNR problem as MDP can be used for solving
similar power system control and optimization problems as well. Once trained,
the RL algorithm demonstrates faster execution compared to the state-of-the-
art method, while yielding approximately equal total cost savings. The
presented results indicate that the developed algorithm is scalable, i.e., the
computation time does not increase exponentially with the problem dimension.
Therefore, the developed DDNR algorithm can handle large-scale real-life
distribution networks well.

The major drawback of applying typical RL algorithms in power systems is the
need for the usage of power system simulators that enable the exploration
process of an RL agent. Training an RL agent on a real-world power system is
not a common practice due to the high cost of exploration actions. On the
other hand, power system simulators do not represent the real-world power
system ideally, which can cause the lower performance of a trained RL model
when deployed in reality. A possible solution and a topic for future work
could be using offline reinforcement learning [172], which eliminates the need
for using the action exploration, by optimizing the policy based on the
variety of historical actions, collected during the real-world distribution
network operation. Another promising area of research is using safe
reinforcement learning [173] and worst-case reinforcement learning [174],
which enables safe exploration during the training while guaranteeing the
constraint satisfaction, therefore causing no harm to the power system.

###  Chapter 8 Conclusions

Deep learning has demonstrated great potential to improve various tasks in
power systems, including monitoring tasks such as stability assessment and
fault detection, as well as for optimization tasks like Volt-Var optimization,
optimal power flow, etc. One of the current trends in the field is the use of
GNNs and DRL. In this thesis, we applied these methods to SE and DDNR problems
and shown that these methods exhibit high levels of accuracy and improved
performance when compared to traditionally used techniques. As the field
continues to evolve, it is expected that more research and development will be
conducted in these areas, with a focus on implementing these techniques in
real-world power systems to demonstrate their practical potential.

This thesis presents two main contributions. As the first one, we investigate
how GNN can be used as fast and accurate solvers of linear and nonlinear SE.
The proposed graph attention network-based model, specialized for the newly
introduced heterogeneous augmented factor graphs, recursively propagates the
input measurements from the factor nodes to generate predictions in the
variable nodes. Evaluating the trained model on unseen data samples confirms
that the proposed GNN approach can be used as a highly accurate approximator
of the traditional SE solutions, with the added benefit of linear
computational complexity at inference time. The model is robust in
unobservable scenarios that are not solvable using traditional SE and deep
learning methods, such as when individual measurements or entire measurement
units fail to deliver measurement data to the proposed SE solver. Furthermore,
the GNN model performs well when measurement variances are high or outliers
are present in the input data. The proposed approach demonstrates scalability
and sample efficiency when tested on power systems of various sizes, as it
makes good predictions even when trained on a small number of randomly
generated samples. Finally, the proposed GNN model outperforms the more
conventional deep learning-based SE approach in terms of prediction accuracy
and significantly lower number of trainable parameters, especially as the size
of the power system grows. In this work, we focused on using GNNs to solve a
linear and nonlinear transmission system SE model. However, the proposed
learning framework, graph augmentation techniques, and conclusions can be
applied to a wide range of SE formulations. For example, the GNN’s ability to
provide relevant solutions in underdetermined scenarios suggests that it could
be useful for GNN-based SE in highly unobservable distribution systems.

As a second main contribution, this thesis has explored a multi-objective
formulation for DDNR, aimed at minimizing the total cost of energy losses and
switching operations. The proposed solution, based on DRL, demonstrated
successful treatment of DDNR as an MDP through the use of state variables with
reduced observability requirements and a computationally efficient approach
for handling switching operation constraints. The results showed that the
developed algorithm is scalable, performs faster than the state-of-the-art
method while yielding comparable cost savings, and is capable of handling
large-scale real-world distribution networks. This work contributes to the
field of power system control and optimization by providing a novel and
effective solution for DDNR.

### References

  * [1] J. R. Aguero, E. Takayesu, D. Novosel, and R. Masiello, “Modernizing the grid: Challenges and opportunities for a sustainable future,” _IEEE Power Energy Mag._ , vol. 15, no. 3, pp. 74–83, 2017.
  * [2] S. Rusitschka, K. Eger, and C. Gerdes, “Smart grid data cloud: A model for utilizing cloud computing in the smart grid domain,” in _Proc. SmartGridComm_. IEEE, 2010, pp. 483–488.
  * [3] I. J. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_. Cambridge, MA, USA: MIT Press, 2016.
  * [4] I. H. Sarker, “Deep learning: A comprehensive overview on techniques, taxonomy, applications and research directions,” _SN comput. sci._ , vol. 2, 2021.
  * [5] K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” _Neural Netw._ , vol. 2, no. 5, p. 359–366, july 1989.
  * [6] J. Quionero-Candela, M. Sugiyama, A. Schwaighofer, and N. D. Lawrence, _Dataset Shift in Machine Learning_. The MIT Press, 2009.
  * [7] G. D’Antona, “Power system static-state estimation with uncertain network parameters as input data,” _IEEE Trans. Instrum. Meas._ , vol. 65, no. 11, pp. 2485–2494, 2016.
  * [8] W. Zhang, D. Yang, and H. Wang, “Data-driven methods for predictive maintenance of industrial equipment: A survey,” _IEEE Syst. J._ , vol. 13, no. 3, pp. 2213–2227, 2019.
  * [9] O. Kundacina, G. Gojic, M. Mitrovic, D. Miskovic, and D. Vukobratovic, “Supporting future electrical utilities: Using deep learning methods in ems and dms algorithms,” in _Proc. Infoteh-Jahorina_ , 2023, pp. 1–6.
  * [10] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _Proc. ICLR_ , Y. Bengio and Y. LeCun, Eds., 2015.
  * [11] J. Bergstra, D. Yamins, and D. D. Cox, “Making a science of model search: Hyperparameter optimization in hundreds of dimensions for vision architectures,” in _Proc. ICML_ , 2013.
  * [12] A. Poudyal, R. S. Fourney, R. Tonkoski, T. M. Hansen, U. Tamrakar, and R. D. Trevizan, “Convolutional neural network-based inertia estimation using local frequency measurements,” _Proc. NAPS_ , pp. 1–6, 2021.
  * [13] M. Alqudah, M. Pavlovski, T. Dokic, M. Kezunovic, Y. Hu, and Z. Obradovic, “Fault detection utilizing convolution neural network on timeseries synchrophasor data from phasor measurement units,” _IEEE Trans. Power Syst._ , vol. 37, no. 5, pp. 3434–3442, 2022.
  * [14] S. Wang, P. Dehghanian, and L. Li, “Power grid online surveillance through pmu-embedded convolutional neural networks,” _IEEE Trans. Ind. Appl._ , vol. 56, no. 2, pp. 1146–1155, 2020.
  * [15] H. Wu and Z. Xu, “Fast dc optimal power flow based on deep convolutional neural network,” in _Proc. CIEEC_. IEEE, 2022, pp. 2508–2512.
  * [16] R. Jozefowicz, W. Zaremba, and I. Sutskever, “An empirical exploration of recurrent network architectures,” in _Proc. ICML_ , 2015, p. 2342–2350.
  * [17] Q. Li, Y. Xu, B. S. H. Chew, H. Ding, and G. Zhao, “An integrated missing-data tolerant model for probabilistic pv power generation forecasting,” _IEEE Trans. Power Syst._ , vol. 37, no. 6, pp. 4447–4459, 2022.
  * [18] G. Ruan, D. S. Kirschen, H. Zhong, Q. Xia, and C. Kang, “Estimating demand flexibility using siamese lstm neural networks,” _IEEE Trans. Power Syst._ , vol. 37, no. 3, pp. 2360–2370, 2022.
  * [19] Y. Seyedi, H. Karimi, and J. Mahseredjian, “A data-driven method for prediction of post-fault voltage stability in hybrid ac/dc microgrids,” _IEEE Trans. Power Syst._ , vol. 37, no. 5, pp. 3758–3768, 2022.
  * [20] D. Fellner, T. I. Strasser, and W. Kastner, “Applying deep learning-based concepts for the detection of device misconfigurations in power systems,” _Sustain. Energy Grids Netw._ , vol. 32, p. 100851, 2022.
  * [21] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Proc. NIPS_ , Red Hook, NY, USA, 2017, p. 6000–6010.
  * [22] Z. Wang, Y. Ma, Z. Liu, and J. Tang, “R-transformer: Recurrent neural network enhanced transformer,” _arXiv preprint arXiv:1907.05572_ , 2019.
  * [23] W. L. Hamilton, “Graph representation learning,” _Synthesis Lectures on Artificial Intelligence and Machine Learning_ , vol. 14, no. 3, pp. 1–159, 2020\.
  * [24] O. Kundacina, M. Forcan, M. Cosovic, D. Raca, M. Dzaferagic, D. Miskovic, M. Maksimovic, and D. Vukobratovic, “Near real-time distributed state estimation via ai/ml-empowered 5g networks,” in _Proc. SmartGridComm_. IEEE, 2022, pp. 284–289.
  * [25] K. Chen, J. Hu, Y. Zhang, Z. Yu, and J. He, “Fault location in power distribution systems via deep graph convolutional networks,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 1, pp. 119–131, 2020.
  * [26] R. Zhang, W. Yao, Z. Shi, L. Zeng, Y. Tang, and J. Wen, “A graph attention networks-based model to distinguish the transient rotor angle instability and short-term voltage instability in power systems,” _Int. J. Electr. Power Energy Syst._ , vol. 137, p. 107783, 2022.
  * [27] S. Arastehfar, M. Matinkia, and M. R. Jabbarpour, “Short-term residential load forecasting using graph convolutional recurrent neural networks,” _Eng. Appl. Artif. Intell._ , vol. 116, p. 105358, 2022.
  * [28] T. Zhao, M. Yue, and J. Wang, “Structure-informed graph learning of networked dependencies for online prediction of power system transient dynamics,” _IEEE Trans. Power Syst._ , vol. 37, no. 6, pp. 4885–4895, 2022.
  * [29] A. Takiddin, R. Atat, M. Ismail, O. Boyaci, K. R. Davis, and E. Serpedin, “Generalized graph neural network-based detection of false data injection attacks in smart grids,” _IEEE Trans. Emerg. Top. Comput. Intell._ , pp. 1–13, 2023.
  * [30] O. Kundacina, M. Cosovic, and D. Vukobratovic, “State estimation in electric power systems leveraging graph neural networks,” in _Proc. PMAPS_ , 2022, pp. 1–6.
  * [31] O. Kundacina, M. Cosovic, D. Miskovic, and D. Vukobratovic, “Distributed nonlinear state estimation in electric power systems using graph neural networks,” in _Proc. SmartGridComm_. IEEE, 2022, pp. 8–13.
  * [32] Y. Yuan, Z. Wang, and Y. Wang, “Learning latent interactions for event classification via graph neural networks and pmu data,” _IEEE Trans. Power Syst._ , vol. 38, no. 1, pp. 617–629, 2023.
  * [33] R. S. Sutton and A. G. Barto, _Reinforcement Learning: An Introduction_ , 2nd ed. The MIT Press, 2018.
  * [34] M. Glavic, “(deep) reinforcement learning for electric power system control and related problems: A short review and perspectives,” _Annu. Rev. Control_ , 2019.
  * [35] Y. Gao, J. Shi, W. Wang, and N. Yu, “Dynamic distribution network reconfiguration using reinforcement learning,” in _Proc. SmartGridComm_. IEEE, 2019, pp. 1–7.
  * [36] W. Wang, N. Yu, Y. Gao, and J. Shi, “Safe off-policy deep reinforcement learning algorithm for volt-var control in power distribution systems,” _IEEE Trans. Smart Grid_ , vol. 11, no. 4, pp. 3008–3018, 2020.
  * [37] O. Stanojev, O. Kundacina, U. Markovic, E. Vrettos, P. Aristidou, and G. Hug, “A reinforcement learning approach for fast frequency control in low-inertia power systems,” in _Proc. NAPS_ , 2021, pp. 1–6.
  * [38] L. Busoniu, R. Babuska, and B. De Schutter, “A comprehensive survey of multiagent reinforcement learning,” _IEEE Trans. Syst. Man Cybern._ , vol. 38, no. 2, pp. 156–172, 2008.
  * [39] H. Liu and W. Wu, “Online multi-agent reinforcement learning for decentralized inverter-based volt-var control,” _IEEE Trans. Smart Grid_ , vol. 12, no. 4, pp. 2980–2990, 2021.
  * [40] D. Chen, K. Chen, Z. Li, T. Chu, R. Yao, F. Qiu, and K. Lin, “Powernet: Multi-agent deep reinforcement learning for scalable powergrid control,” _IEEE Trans. Power Syst._ , vol. 37, no. 2, pp. 1007–1017, 2022.
  * [41] L. Ding, Z. Lin, X. Shi, and G. Yan, “Target-value-competition-based multi-agent deep reinforcement learning algorithm for distributed nonconvex economic dispatch,” _IEEE Trans. Power Syst._ , vol. 38, no. 1, pp. 204–217, 2023.
  * [42] T. Qian, C. Shao, X. Li, X. Wang, Z. Chen, and M. Shahidehpour, “Multi-agent deep reinforcement learning method for ev charging station game,” _IEEE Trans. Power Syst._ , vol. 37, no. 3, pp. 1682–1694, 2022.
  * [43] A. Monticelli, “Electric power system state estimation,” _Proc. IEEE_ , vol. 88, no. 2, pp. 262–282, 2000.
  * [44] A. Gomez-Exposito, A. Abur, P. Rousseaux, A. de la Villa Jaen, and C. Gomez-Quiles, “On the use of PMUs in power system state estimation,” in _Proc. PSCC_ , 2011.
  * [45] G. N. Korres and N. M. Manousakis, “State estimation and observability analysis for phasor measurement unit measured systems,” _IET Gener. Transm. Dis._ , vol. 6, no. 9, pp. 902–913, September 2012.
  * [46] A. Abur and A. Expósito, _Power System State Estimation: Theory and Implementation_ , ser. Power Engineering. Taylor & Francis, 2004.
  * [47] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl, “Neural message passing for quantum chemistry,” in _Proc. ICML_ , 06–11 Aug 2017, pp. 1263–1272.
  * [48] L. Zhang, G. Wang, and G. B. Giannakis, “Real-time power system state estimation and forecasting via deep unrolled neural networks,” _IEEE Trans. Signal Process._ , vol. 67, no. 15, pp. 4069–4077, 2019.
  * [49] A. S. Zamzam, X. Fu, and N. D. Sidiropoulos, “Data-driven learning-based optimization for distribution system state estimation,” _IEEE Trans. Power Syst._ , vol. 34, no. 6, pp. 4796–4805, 2019.
  * [50] B. Donon, B. Donnot, I. Guyon, and A. Marot, “Graph neural solver for power systems,” in _Proc. IJCNN_ , 2019, pp. 1–8.
  * [51] V. Bolz, J. Rueß, and A. Zell, “Power flow approximation based on graph convolutional networks,” in _Proc. ICMLA_. IEEE, 2019, pp. 1679–1686.
  * [52] D. Wang, K. Zheng, Q. Chen, G. Luo, and X. Zhang, “Probabilistic power flow solution with graph convolutional network,” in _Proc. ISGT Europe_. IEEE, 2020, pp. 650–654.
  * [53] B. Donon, R. Clément, B. Donnot, A. Marot, I. Guyon, and M. Schoenauer, “Neural networks for power flow: Graph neural solver,” _Electr. Power Syst. Res._ , vol. 189, p. 106547, 2020.
  * [54] T. B. Lopez-Garcia and J. A. Domínguez-Navarro, “Power flow analysis via typed graph neural networks,” _Eng. Appl. Artif. Intell._ , vol. 117, p. 105567, 2023.
  * [55] L. Pagnier and M. Chertkov, “Physics-informed graphical neural network for parameter & state estimations in power systems,” _arXiv_ , 2021.
  * [56] Q. Yang, A. Sadeghi, and G. Wang, “Data-driven priors for robust psse via gauss-newton unrolled neural networks,” _IEEE J. Emerg. Sel. Top. Circuits Syst._ , vol. PP, pp. 1–1, 01 2022.
  * [57] M. J. Hossain and M. Rahnamay–Naeini, “State estimation in smart grids using temporal graph convolution networks,” in _Proc. NAPS_ , 2021, pp. 01–05.
  * [58] O. Kundacina, M. Cosovic, D. Miskovic, and D. Vukobratovic, “Graph neural networks on factor graphs for robust, fast, and scalable linear state estimation with PMUs,” _Sustain. Energy Grids Netw._ , vol. 34, p. 101056, 2023.
  * [59] V. G. Satorras and M. Welling, “Neural enhanced belief propagation on factor graphs,” in _Proc. AISTATS_ , 2021.
  * [60] F. Kschischang, B. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” _IEEE Trans. Inf. Theory_ , vol. 47, no. 2, pp. 498–519, 2001.
  * [61] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio, “Graph Attention Networks,” in _Proc. ICLR_ , 2018.
  * [62] M. A. Samman, H. Mokhlis, N. N. Mansor, H. Mohamad, H. Suyono, and N. M. Sapari, “Fast optimal network reconfiguration with guided initialization based on a simplified network approach,” _IEEE Access_ , vol. 8, pp. 11 948–11 963, 2020.
  * [63] V. Fathi, H. Seyedi, and B. M. Ivatloo, “Reconfiguration of distribution systems in the presence of distributed generation considering protective constraints and uncertainties,” _Int. Trans. Electr. Energy Syst._ , vol. 30, no. 5, p. e12346, 2020.
  * [64] S. Mishra, D. Das, and S. Paul, “A comprehensive review on power distribution network reconfiguration,” _Energy Syst._ , vol. 8, pp. 227–284, 2017.
  * [65] S. Civanlar, J. Grainger, H. Yin, and S. Lee, “Distribution feeder reconfiguration for loss reduction,” _IEEE Trans. Power Deliv._ , vol. 3, no. 3, pp. 1217–1223, 1988.
  * [66] D. Shirmohammadi and H. Hong, “Reconfiguration of electric distribution networks for resistive line losses reduction,” _IEEE Trans. Power Deliv._ , vol. 4, no. 2, pp. 1492–1498, 1989.
  * [67] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution systems for loss reduction and load balancing,” _IEEE Power Eng. Rev._ , vol. 9, no. 4, pp. 101–102, 1989.
  * [68] R. Taleski and D. Rajicic, “Distribution network reconfiguration for energy loss reduction,” _IEEE Trans. Power Syst._ , vol. 12, no. 1, pp. 398–406, 1997.
  * [69] Q. Zhou, D. Shirmohammadi, and W.-H. Liu, “Distribution feeder reconfiguration for operation cost reduction,” _IEEE Trans. Power Syst._ , vol. 12, no. 2, pp. 730–735, 1997.
  * [70] A. Borghetti, “A mixed-integer linear programming approach for the computation of the minimum-losses radial configuration of electrical distribution networks,” _IEEE Trans. Power Syst._ , vol. 27, no. 3, pp. 1264–1273, 2012\.
  * [71] H. Ahmadi and J. R. Martí, “Distribution system optimization based on a linear power-flow formulation,” _IEEE Trans. Power Deliv._ , vol. 30, no. 1, pp. 25–33, 2015.
  * [72] M. Lavorato, J. F. Franco, M. J. Rider, and R. Romero, “Imposing radiality constraints in distribution system optimization problems,” _IEEE Trans. Power Syst._ , vol. 27, no. 1, pp. 172–180, 2012.
  * [73] R. A. Jabr, R. Singh, and B. C. Pal, “Minimum loss network reconfiguration using mixed-integer convex programming,” _IEEE Trans. Power Syst._ , vol. 27, no. 2, pp. 1106–1115, 2012.
  * [74] H. Haghighat and B. Zeng, “Distribution system reconfiguration under uncertain load and renewable generation,” _IEEE Trans. Power Syst._ , vol. 31, no. 4, pp. 2666–2675, 2016.
  * [75] E. Ramos, A. Exposito, J. Santos, and F. Iborra, “Path-based distribution network modeling: application to reconfiguration for loss reduction,” _IEEE Trans. Power Syst._ , vol. 20, no. 2, pp. 556–564, 2005.
  * [76] H. M. Khodr, J. Martinez-Crespo, M. A. Matos, and J. Pereira, “Distribution systems reconfiguration based on opf using benders decomposition,” _IEEE Trans. Power Deliv._ , vol. 24, no. 4, pp. 2166–2176, 2009.
  * [77] A. Merlin and H. Back, “Search for a minimal-loss operating spanning tree configuration in an urban power distribution system,” in _Proc. PSCC_ , vol. 1, 09 1975.
  * [78] V. Borozan, D. Rajicic, and R. Ackovski, “Improved method for loss minimization in distribution networks,” _IEEE Trans. Power Syst._ , vol. 10, no. 3, pp. 1420–1425, 1995.
  * [79] C. T. Huddleston, R. P. Broadwater, and A. Chandrasekaran, “Reconfiguration algorithm for minimizing losses in radial electric distribution systems,” _Electr. Power Syst. Res._ , vol. 18, no. 1, pp. 57–66, 1990.
  * [80] I. Roytelman, V. Melnik, S. Lee, and R. Lugtu, “Multi-objective feeder reconfiguration by distribution management system,” _IEEE Trans. Power Syst._ , vol. 11, no. 2, pp. 661–667, 1996.
  * [81] N. Gupta, A. Swarnkar, and K. R. Niazi, “Distribution network reconfiguration for power quality and reliability improvement using genetic algorithms,” _Int. J. Electr. Power Energy Syst._ , vol. 54, pp. 664–671, 2014.