Observability Analysis for Large-Scale Power Systems Using Factor Graphs

TL;DR

This paper introduces a novel observability analysis method for large-scale power systems using factor graphs and Gaussian belief propagation, enabling efficient and distributed determination of observable islands and measurement placement.

Contribution

It pioneers the use of factor graphs and belief propagation for power system observability analysis, offering a scalable and distributed approach for large networks.

Findings

Method scales linearly with system size due to sparsity.

Effective identification of observable islands and measurement placement.

Applicable to power systems with up to 70,000 buses.

Abstract

The state estimation algorithm estimates the values of the state variables based on the measurement model described as the system of equations. Prior to applying the state estimation algorithm, the existence and uniqueness of the solution of the underlying system of equations is determined through the observability analysis. If a unique solution does not exist, the observability analysis defines observable islands and further defines an additional set of equations (measurements) needed to determine a unique solution. For the first time, we utilise factor graphs and Gaussian belief propagation algorithm to define a novel observability analysis approach. The observable islands and placement of measurements to restore observability are identified by following the evolution of variances across the iterations of the Gaussian belief propagation algorithm over the factor graph. Due to sparsity of the underlying power network, the resulting method has the linear computational complexity (assuming a constant number of iterations) making it particularly suitable for solving large-scale systems. The method can be flexibly matched to distributed computational resources, allowing for determination of observable islands and observability restoration in a distributed fashion. Finally, we discuss performances of the proposed observability analysis using power systems whose size ranges between 1354 and 70000 buses.

Equations

$$\mathbf{Hx}=\mathbf{z},$$

$$\text{rank}(\mathbf{H})=n-1,$$

$$\mathbf{H}=\blockarray{*{6}{c}l}\block{*{6}{>{\footnotesize}c<{}}l}x_{1}&x_{2}x_{3}x_{4}x_{5}x_{6}\\ \block{[*{6}{r}]>{\footnotesize}l<{}}1&-10000\bigstrut[t]M_{P_{12}}\\ 0001-10M_{P_{45}}\\ -1-13-100M_{P_{3}}\\ 000-12-1M_{P_{5}}\\$$

$$v_{f_j\to x_i}^{(\tau )}=\sum _{x_b\in {\mathcal{X}}_j\setminus x_i}{v}_{x_b\to f_j}^{(\tau -1)}$$

$${\left(v_{x_i\to f_j}^{(\tau )}\right)}^{-1}=\sum _{f_a\in {\mathcal{F}}_i\setminus f_j}{\left(v_{f_a\to x_i}^{(\tau -1)}\right)}^{-1}$$

$$\begin{array}{c}{\left({\hat{v}}_{x_i}\right)}^{-1}=\sum _{f_a\in {\mathcal{F}}_i}{\left(v_{f_a\to x_i}^{({\tau }_{\max })}\right)}^{-1}\end{array}$$

$$v_{f_1\to x_1}$$

$$v_{f_i\to x_i}$$

$$v_{x_i\to f_i}^{(0)}=1,x_i\in \mathcal{X},f_i\in \mathcal{F}.$$

$$v_{f_9\to x_1}^{(1)}=v_{x_2\to f_9}^{(0)}+v_{x_3\to f_9}^{(0)}+v_{x_4\to f_9}^{(0)}=3.$$

$$\frac{1}{v_{x_1\to f_7}^{(1)}}=\frac{1}{v_{f_9\to x_1}^{(1)}}+\frac{1}{v_{f_1\to x_1}};v_{x_1\to f_7}^{(1)}\to 0.$$

$$\frac{1}{{\hat{v}}_{x_1}}=\frac{1}{v_{f_7\to x_1}^{({\tau }_{\max })}}+\frac{1}{v_{f_9\to x_1}^{({\tau }_{\max })}}+\frac{1}{v_{f_1\to x_1}};{\hat{v}}_{x_1}\to 0$$

$$\frac{1}{{\hat{v}}_{x_2}}=\frac{1}{v_{f_7\to x_2}^{({\tau }_{\max })}}+\frac{1}{v_{f_9\to x_2}^{({\tau }_{\max })}}+\frac{1}{v_{f_2\to x_2}};{\hat{v}}_{x_2}\to 0,$$

$${\mathbf{H}}_{\mathrm{u}}\mathbf{x}={\mathbf{z}}_{\mathrm{u}},$$

$$\left[\begin{array}{c}{\mathbf{H}}_{\text{c}}\\ {\mathbf{H}}_{\text{p}}\end{array}\right]\mathbf{x}=\left[\begin{array}{c}{\mathbf{z}}_{\text{c}}\\ {\mathbf{z}}_{\text{p}}\end{array}\right],$$

$$\left[\begin{array}{c}{\mathbf{H}}_{\text{u}}\\ {\mathbf{H}}_{\text{c}}\\ {\mathbf{H}}_{\text{p}}\end{array}\right]\mathbf{x}+\left[\begin{array}{c}\text{upalpha}\\ \text{upbeta}\\ \text{upgamma}\end{array}\right]=\left[\begin{array}{c}{\mathbf{z}}_{\text{u}}\\ {\mathbf{z}}_{\text{c}}\\ {\mathbf{z}}_{\text{p}}\end{array}\right].$$

$$\left[\begin{array}{cc}\left[\begin{array}{c}{\mathbf{H}}_{\mathrm{u}}\\ {\mathbf{H}}_{\text{c}}\\ {\mathbf{H}}_{\text{p}}\end{array}\right]& {\mathbf{I}}_{d\times d}\\ 0_{d\times n}& {\mathbf{I}}_{d\times d}\end{array}\right]\left[\begin{array}{c}\mathbf{x}\\ \mathbf{r}\end{array}\right]+\left[\begin{array}{c}\text{upalpha}\\ \text{upbeta}\\ \text{upgamma}\\ \text{updelta}\end{array}\right]=\left[\begin{array}{c}{\mathbf{z}}_{\text{u}}\\ {\mathbf{z}}_{\text{c}}\\ {\mathbf{z}}_{\text{p}}\\ 0_{d\times 1}\end{array}\right],$$

$$v_{f_j\to s_i}^{(\tau )}=v_j+\sum _{s_b\in {\mathcal{S}}_j\setminus s_i}{v}_{s_b\to f_j}^{(\tau -1)}$$

$${\left(v_{s_i\to f_j}^{(\tau )}\right)}^{-1}=\sum _{f_a\in {\mathcal{F}}_i\setminus f_j}{\left(v_{f_a\to s_i}^{(\tau -1)}\right)}^{-1}$$

$$v_{f_1\to s_1}$$

$$v_{f_i\to s_i}$$

$$v_{f_{12}\to r_{12}}^{({\tau }_{\max })}=v_{s_1\to f_{12}}^{({\tau }_{\max })}+v_{s_3\to f_{12}}^{({\tau }_{\max })};v_{f_{12}\to r_{12}}^{({\tau }_{\max })}\to \infty .$$

Full text

I Introduction

Motivation: Deregulation of the electric power industry has led to the development of distributed control centres that operate over segments of a large interconnected power system. Such evolution of power systems calls for the distributed monitoring and control operations [1]. As a part of the energy management system, power system monitoring relies on the distributed, scalable and computationally efficient state estimation (SE) routines that include the network topology processor, observability analysis, SE algorithm, and bad data processing.

Typically, the observability analysis enacts before the SE algorithm with bad data processing. According to the available measurements, observability analysis determines whether the network is observable or unobservable. Within an unobservable network, the observability analysis identifies all observable islands that can be solved independently, and obtains a minimal set of additional measurements providing the observability of the entire system, which ensures a unique estimate of the power system state [2]. Over the past decades, comprehensive studies have been published that predominantly rely on distributed SE algorithms with bad data processing. In contrast, computationally efficient and distributed observability analysis received significantly less attention from the research community.

Literature Review: The identification of observable islands and observability restoration is traditionally carried out by applying numerical or topological methods. Numerical methods for the identification of observable islands and the observability restoration are mostly based on the matrix factorisation. For the identification of observable islands, there are several methods depending on the implementation approach: i) iterative methods that use the Jacobian [2], the gain [3, 4], and the Gram matrix [5]; ii) non-iterative methods that use the gain matrix with the inverse of its triangular matrix [6], and the Jacobian matrix based on the branch variable formulation [7]. Similar categorisation can be made for the numerical methods applied to the observability restoration: i) iterative methods that use the gain [6, 3], the Gram [5], and reduced network matrix [8]; ii) non-iterative methods that use the Jacobian matrix [2], the gain matrix with the inverse of its triangular matrix [9], and the Gram matrix of the reduced Jacobian [10, 11]. These approaches depend on the accurate computation of zero pivots. The major drawback of these methods is the sensitivity to the zero pivot threshold value, where the incorrect value may deteriorate observability analysis. Also, non-iterative methods usually require definition and the calculation of additional matrices as compared to the iterative methods. Methods that use the Gram matrices, even for the reduced system, are prone to the risk of losing sparsity.

Unlike numerical methods, topological methods are less analysed with the most of the research efforts focused on the identification of the observable islands. These methods do not use any floating point calculations in the analysis. The decision is strictly based on the logical operations, and therefore requires information about the network connectivity, measurement types and their locations [12, Sec. 4.7]. Typically, topological methods first process all flow measurements to define initial islands, moving to processing injection measurements in an iterative fashion, followed by the identification of observable islands [13, 8, 14, 15]. In the paper [16], authors propose a method based on the symbolic Jacobian matrix which considers only non-redundant measurements. Starting from a random bus, the proposed method utilise search technique to reveal observable islands through repetitive process.

Although topological methods can be adopted in a distributed framework, little work has been done in this area. However, numerical methods are also appealing for the observability analysis in the distributed framework, with many novel methods developed over the years. Several studies [17, 18] proposed the distributed observability analysis based on the numerical methods that require coordination among areas, where observability restoration is accomplished by a non-iterative approach. Authors in [19] analysed the distributed observability in terms of a consensus algorithm, providing conditions which enable solving the system of equations.

Belief Propagation Approach and Contributions: In this paper, we present a novel observability analysis approach based on a probabilistic graphical model, complementing our previous work on Gaussian belief propagation (GBP) based SE [20]. The proposed method can be understood as both numerical and topological, with several advantages over the current state-of-the-art approaches. Our contributions are summarised as follows:

• •

Compared to numerical methods, the proposed method does not utilise direct factorisation or inversion of matrices, avoiding inaccurate computation of zero pivots and incorrect choice of a zero threshold value.

• •

Due to the sparsity of the underlying factor graph, the proposed method has linear computational complexity per GBP iteration, making it particularly suitable for large-scale systems.

• •

In the multi-area scenario, the proposed method can be implemented over the non-overlapping multi-area scenario without the central coordinator. The GBP algorithm ensures measurement data privacy in the distributed architecture, where areas exchange only “beliefs” about specific state variables. The GBP algorithm within local areas can be implemented using synchronous scheduling, where all messages in a given iteration are updated using the output of the previous iteration as an input. Message updates between local areas are not necessarily performed in full synchronism with the scheduling applied within local areas, as this would incur excessive communication delay to each iteration. Thus, the GBP framework allows the algorithm in which asynchronous information exchange is possible [21].

• •

The method is straightforward to distribute and parallelize. In particular, even if implemented in the framework of centralised SE, it can be flexibly mapped to distributed computation resources (e.g., parallel processing on graphical-processing units).

The novelty of our work follows from the fact that, for the first time, we adopt the GBP algorithm to carry out the observability analysis. Moreover, we emphasise a particularly simplifying feature of the proposed method, where within the GBP framework, it is sufficient to track the evolution of the variance values only. In particular, the GBP algorithm convergence depends on the convergence of the means, while variances always converge [22]. To leverage this property, we propose observability analysis based only on variances resulting in an algorithm that always converges. The resulting algorithm relies only on the graph structure and connections between nodes, and in contrast to numerical methods, it is immune to numerical issues.

II Observability Analysis

This section provides a brief review of key concepts used for observability analysis in the power systems. Observability analysis is commonly performed on the linear decoupled measurement model [23, Ch. 7], which we adopt in this paper. In summary, the observability analysis considers the following linearised measurement model:

[TABLE]

where $\mathbf{z}\in\mathbb{R}^{m}$ is the vector of known values, $\mathbf{x}\in\mathbb{R}^{n}$ is the solution, while $\mathbf{H}\in\mathbb{R}^{m\times n}$ is a known matrix. The system is observable if and only if [9]:

[TABLE]

where the state variable corresponding to the slack bus is already included in (1). Ensuring the observability by condition (2) does not imply achieving a good estimation of the state variables. This is because numerical ill-conditioning may deteriorate the SE algorithm, but in the most cases, if the system is observable, it can be efficiently estimated [24]. Note that the requirement posed by the condition (2) implicitly assumes using the least-squares approach in solving (1). In theory, requirement (2) is not sufficient for the system (1) to have a solution (the existence of the solution depends on the vector $\mathbf{z}$). However, (2) guarantees invertibility of the symmetric matrix $\mathbf{H}^{{T}}\mathbf{H}\in\mathbb{R}^{n\times n}$, that is at the heart of the least-squares approach, therefore making (1) solvable.

When the given set of measurements is not sufficient (i.e., $\text{rank}(\mathbf{H})<n-1$), the observability analysis must identify all the possible observable islands that can be independently solved [2], where an observable island is defined as follows:

Definition 1**.**

An observable island is a part of the power system for which the flows across all branches of the observable island can be calculated from the set of available measurements, independent of the values adopted for angular reference [23, Sec. 7.1.1].

In other words, the observability analysis divides the system into the set of islands $\mathcal{S}=\{s_{1},\dots,s_{k}\}$, where $k$ is the number of islands. Formally, the set of state variables $\mathcal{X}=\{x_{1},\dots,x_{n}\}$ is divided into disjoint subsets $\mathcal{X}_{i}\subset\mathcal{X}$, $i=1,\dots,k$. Each island encloses $n_{i}=|\mathcal{X}_{i}|$ state variables, where the number of linearly independent equations $l_{i}$ corresponding to the island $s_{i}$, is equal to $l_{i}=n_{i}-1$. Once the islands are determined, the observability analysis merges these islands in a way to protect previously-determined observable states from being altered by the new set of equations defined by the additional measurements. In general, this can be achieved by ensuring that the set of new measurements is a non-redundant set [23, Sec. 7.3.2], i.e., the set of equations must be linearly independent with regard to the global system. The aim of the observability restoration is to find this non-redundant set. The observability restoration finds exactly $k-1$ linearly independent equations, where the equation that defines the slack bus is known in advance.

III Belief Propagation Based Observability Analysis

We consider a factor graph approach for modelling the system of equations [25]. This approach starts with the construction of the factor graph with the set of variable nodes representing unknown variables, and the set of factor nodes corresponding to the set of equations. Note that a factor node connects to a variable node if and only if the unknown variable associated with the variable node is an argument of the corresponding equation associated with the factor node. Alongside the factor nodes that represent system equations (measurements), we introduce an additional set of virtual factor nodes defined as follows:

Definition 2**.**

The virtual factor node is a singly-connected factor node attached to the variable node, with the variance set111With a slight abuse of terminology, hereinafter we use the term “set”. However, one could argue that “tend” could be deemed to be more appropriate. either to zero or to infinity.

III-A Determination of Observable Islands

In our setup, from the set of equations, we retain a subset that corresponds to the set of the linearly independent equations defined by the power flow and the power injection measurement functions. We start with the construction of the factor graph with the set of variable nodes $\mathcal{X}$ representing state variables, and the set of factor nodes $\mathcal{F}$ corresponding to the linearly independent set of equations. Utilising the property that the GBP is robust to ill-conditioned scenarios caused by significant differences between measurement variances [20, 26], we augment the set of factor nodes $\mathcal{F}$ with the set of virtual factor nodes $\mathcal{F}_{\text{x}}$. Note that virtual factor nodes are equivalent to the pseudo-measurements of the voltage angles at the buses. Furthermore, if its variance tends to the infinity, the factor node has no influence on the rest of the graph.

Example 1** (Constructing a factor graph).**

In this toy example, using 6-bus model presented in Figure 1, we demonstrate the conversion from a bus/branch model with a given measurement configuration to the corresponding factor graph.

The linearised measurement Jacobian matrix can be formed as:

[TABLE]

The corresponding factor graph is given in Figure 2, where the set of variable nodes $\mathcal{X}=\{x_{1},x_{2},x_{3},x_{4},x_{5},x_{6}\}$ is defined according to the vector of state variables $\mathbf{x}$. The set of factor nodes $\mathcal{F}$ is defined by corresponding measurements, where in our example, active power flows $M_{P_{12}}$, $M_{P_{45}}$, and active power injections $M_{P_{3}}$, $M_{P_{5}}$ measurements are mapped into factor nodes $\mathcal{F}=\{f_{7},f_{8},f_{9},f_{10}\}$ (red boxes), respectively. The factor graph reflects the structure of the matrix $\mathbf{H}$, where each row of the matrix corresponds to one factor node, while columns define variable nodes. Finally, we augment the factor graph with the set of virtual factor nodes $\mathcal{F}_{\text{x}}=\{f_{1},f_{2},f_{3},f_{4},f_{5},f_{6}\}$ (orange and green boxes), where the factor node $f_{1}$ is attached to the slack bus.

Gaussian Belief-Propagation is a well-known algorithm for probabilistic inference on factor graphs [27]. It is used to efficiently calculate marginal distributions of random variables in a large-scale, sparse, linear systems of random variables corrupted by Gaussian noise. In this paper, we adapt GBP in order to solve the problem of observability analysis in power systems. We present GBP-based island detection in this subsection and extend it to GBP-based observability restoration in the next subsection.

The GBP-based island detection algorithm tracks the evolution of marginal variances of the state variables, establishing observability criteria in the relation to whether or not these variances converge to a unique fixed point. Let us assume, for the time being, that the system is a fully observable, and the Jacobian matrix is defined according to power flows and injections measurement functions. According to Definition 1, the values of the state variables corresponding to the buses are uniquely determined if at least one of the state variables in this set is declared known (i.e., introducing voltage angle pseudo-measurement). As a result, the Jacobian matrix has a full rank [23, Sec. 7.1.3], while the marginal GBP variances converge to a unique fixed point. Indeed, according to Lemma 1 of [28] and Theorem 2 of [29], variances always converge to a unique fixed point if the fixed point exists. However, when the system is unobservable, using GBP-based approach it is still possible to identify subsets of state variables that converge to a unique point [29]. We explore this feature for the identification of observable islands.

We formalise the proposed setup used to detect observable islands in Algorithm 1. We consider a factor graph with the set of variable nodes $\mathcal{X}$, the set of factor nodes $\mathcal{F}$, and the set of virtual factor nodes $\mathcal{F}_{\text{x}}$ constructed as follows. The set $\mathcal{F}_{\text{x}}$ contains $n$ factor nodes, out of which a single one has a variance set to zero and is attached to an arbitrary (randomly selected) variable node (we call this factor node a “probe”), and the remaining $n-1$ factor nodes have variances set to infinity, and are attached to the remaining $n-1$ variable nodes. After the initialisation (lines 1-9), the identification of observable islands procedure (lines 10-40) starts by selecting a random variable node from the set $\mathcal{X}^{(\nu)}$ to which a probe factor node is attached, where $\nu$ denotes the number of passes through the procedure (line 12). The message-passing starts by setting initial GBP variances from variable nodes $x_{i}\in\mathcal{X}^{(\nu)}$ and factor nodes $f_{j}\in\mathcal{F}^{(\nu)}$ to finite positive values (line 13). The message-passing framework with marginal inference (lines 14-28) represents the main algorithm routine which determines the observable islands (lines 29-34), where the set $\mathcal{F}_{i}\setminus f_{j}$ represents factor nodes incident to the variable node $x_{i}$, excluding the factor node $f_{j}$, while the set $\mathcal{X}_{j}\setminus x_{i}$ represents variable nodes incident to the factor node $f_{j}$, excluding the variable node $x_{i}$. Upon the calculation of the marginal variances for all the variable nodes in $\mathcal{X}^{(\nu)}$, a node $x_{i}\in\mathcal{X}^{(\nu)}$ is proclaimed observable, i.e., $x_{i}\in\mathcal{X}_{\text{o}}^{(\nu)}$, if and only if its marginal variance $\hat{v}_{x_{i}}$ is finite. Note that the resulting observable island depends on the position of the probe factor node. After the discovery of observable variable nodes, we consider a subgraph of the factor graph that contains the set of unobservable variable nodes $\mathcal{X}^{(\nu+1)}=$ $\mathcal{X}^{(\nu)}\setminus$ $\mathcal{X}_{\text{o}}^{(\nu)}$ and the set of all factor nodes $\mathcal{F}^{(\nu+1)}$ and $\mathcal{F}_{\text{x}}^{(\nu+1)}$ connected exclusively to variable nodes in $\mathcal{X}^{(\nu+1)}$ (lines 35-37). The procedure now repeats from the line 11. The next observable island is revealed by introducing a new randomly selected probe, and setting its virtual factor node variance to zero. This causes a change in marginal variances, in a way that all unobservable nodes within the next revealed observable island switch over into finite marginal variance values. The procedure is repeated until all observable islands are revealed.

Example 2** (Identification of Observable Islands).**

We demonstrate the GBP-based island detection algorithm, using the factor graph given in Example 1. We start by selecting the probe factor node $f_{1}$ and setting its variance to zero, while variances of the rest of virtual factor nodes are set to infinity:

[TABLE]

The algorithm starts with initial finite (e.g., equal to one) variance values:

[TABLE]

The GBP iteration $\tau=1$ starts with computing variances from factor nodes $\mathcal{F}$ to variable nodes $\mathcal{X}$, using incoming variances from variable nodes $\mathcal{X}$ to factor nodes $\mathcal{F}$. For example, computation of the variance from the factor node $f_{9}$ to the variable node $x_{1}$ is computed as follows:

[TABLE]

Next, the algorithm proceeds with computing variances from variable nodes $\mathcal{X}$ to factor nodes $\mathcal{F}$, using incoming variances from factor nodes $\mathcal{F}\cup\mathcal{F}_{\text{x}}$ to variable nodes $\mathcal{X}$, for example:

[TABLE]

When the algorithm converges, all variances from factor nodes $\mathcal{F}$ to variable nodes $\mathcal{X}$ are equal to infinity, except the variance from the factor node $f_{7}$ to the variable node $x_{2}$ that tends to zero $v_{f_{7}\to x_{2}}^{(\tau_{\max})}\to 0$. The marginal variance of variable nodes $\mathcal{X}$ can be obtained using variances from factor nodes $\mathcal{F}\cup\mathcal{F}_{\text{x}}$ to variable nodes $\mathcal{X}$:

[TABLE]

[TABLE]

while other marginals are equal to infinity. Hence, the first island $s_{1}$ is defined by state variables $\mathcal{X}_{1}=\{x_{1},x_{2}\}$.

Next, we observe the modified network and the measurement configuration shown in Fig. 3, where the new probe factor node is $f_{4}$. We repeat the process revealing the second island $s_{2}$ defined by state variables $\mathcal{X}_{2}=\{x_{4},x_{5},x_{6}\}$. Finally, the third island $s_{3}$ encapsulates the state variable $\mathcal{X}_{3}=\{x_{3}\}$.

Astonishingly, the GBP-based method for the identification of observable islands is similar to the decoding of low-density parity-check codes over erasure channels [30], which may inspire further work in this domain.

III-B Observability Restoration

For the observability restoration, we consider the unobservable system divided into $k$ islands, where each island has $l_{i}=n_{i}-1$ linearly independent equations. Putting the set of linear equations from all $k$ observable islands together, we have:

[TABLE]

where $\mathbf{H}_{\mathrm{u}}\in\mathbb{R}^{h\times n}$, $\mathbf{x}\in\mathbb{R}^{n}$, $\mathbf{z}_{\mathrm{u}}\in\mathbb{R}^{h}$, and $h=n-k$. In addition, we introduce a new set of equations containing all the candidate measurements that are able to merge islands:

[TABLE]

where $\mathbf{H}_{\mathrm{c}}\in\mathbb{R}^{q\times n}$ and $\mathbf{z}_{\mathrm{c}}\in\mathbb{R}^{q}$ represent measurements which were rejected during observable island identification, while $\mathbf{H}_{\mathrm{p}}\in\mathbb{R}^{b\times n}$ and $\mathbf{z}_{\mathrm{p}}\in\mathbb{R}^{b}$ describe pseudo-measurements. The system of equations described with $\mathbf{H}_{\mathrm{c}}$ contains only linearly independent equations. Furthermore, we note that additional $w=k-q-1$ equations among equations defined by pseudo-measurements are needed to achieve observability, as the equation that defines the slack bus is known in advance.

Let us introduce the vectors $\bm{\upalpha}$, $\bm{\upbeta}$ and $\bm{\upgamma}$ representing an uncorrelated noise with the independent and identically distributed entries, where each element follows a zero-mean Gaussian distribution with associated variances $\alpha_{i}$, $\beta_{i}$ and $\gamma_{i}$ Hence, systems (10) and (11) can be merged together as:

[TABLE]

Next, we apply a simple trick using the values of variances assigned to different measurements. Namely, we first set $\alpha_{i}\to 0$, $i=1,\dots,h$, $\beta_{i}\to 0$, $i=1,\dots,q$, and $\gamma_{i}\to\infty$, $i=1,\dots,b$. Then, we select exactly $w$ linearly independent equations (with regards to the global system) from the set of pseudo-measurements. Without loss of generality, we choose the first $w$ equations and we set $\gamma_{i}\to 0$, $i=1,\dots,w$, i.e., we change the variance of $w$ linearly independent measurements from infinity to zero. After additional $w$ equations are “denoised”, the solution of (12) will be determined only by the equations associated to the variances that are set to zero. Note that the corresponding residuals are also equal to zero. In other words, each of the noiseless equations is a critical equation [12, Sec. 5.2]. However, the described scenario cannot be solved using the traditional weighted least-squares method, because we defined the ill-conditioned scenario caused by significant differences between variances. In contrast, the GBP approach is robust to ill-conditioned scenarios, and above-mentioned scenario can be solved [26].

Introducing the additional vector of measurement residuals $\mathbf{r}\in\mathbb{R}^{d}$, $d=h+q+b$, we can expand (12) as follows:

[TABLE]

where $\mathbf{I}$ is the identity matrix, $\mathbf{0}_{d\times n}$ is the all-zero matrix, $\mathbf{0}_{d\times 1}$ is the vector of zeros, and the vector $\bm{\updelta}$ represents an uncorrelated noise defined in the same way as vectors $\bm{\upalpha}$, $\bm{\upbeta}$ and $\bm{\upgamma}$. By applying the weighted least-squares method to (13) and inspecting the variance-covariance matrix, it is easy to show that the variances of the state variables $\mathbf{r}$ depend only on the variances related to the noise vector $\bm{\updelta}$ in the case when the system (12) contains the set of linearly independent equations, and therefore (12) does not influence the residual values. In other words, if an equation described by a factor node $f_{i}$, sends infinity variance to the corresponding residual described by a variable node $r_{i}$, the equation is linearly independent.

For example, consider a part of the factor graph shown in Figure 4, and focus on a single measurement that defines the factor node $f_{i}$. Note that the message from a factor node to a variable node represents a “belief” about the variable node, and each message exchanged in GBP is completely represented using only two values: mean and variance. Messages from variable nodes $\{x_{l},\dots,x_{L}\}$ to the factor node $f_{i}$ carry a collective evidence of the rest of the factor graph, and according to these messages, the GBP algorithm computes the message $\mu_{f_{i}\to r_{i}}(r_{i})$ from the factor node $f_{i}$ to the variable node $r_{i}$. Next, we set the mean value of the message $\mu_{f_{a}\to r_{i}}(r_{i})$ from the singly-connected factor node $f_{a}$ to the variable node $r_{i}$ to zero, and its variance to a finite positive value. The marginal distribution of the variable node $r_{i}$ is obtained as the product of the incoming messages $\mu_{f_{i}\to r_{i}}(r_{i})$ and $\mu_{f_{a}\to r_{i}}(r_{i})$. The mean of the marginal will remain zero only if the variance of the message $\mu_{f_{i}\to r_{i}}(r_{i})$ tends to infinity. When this is the case, the corresponding equation is linearly independent and does not affect the residual value.

To summarize, the GBP-based algorithm decomposes the contribution of each factor node (equation) to a group of variable nodes (state variables and a measurement residual). This decomposition allows insight into the structure of the measurement residual, where the impact of each measurement can be observed. Hence, if one factor node $f_{i}$ sends finite variance value to the residual $r_{i}$, this means that the equation has an impact on the measurement residual value, i.e., $r_{i}\neq 0$, while infinite variance corresponds to the case when the equation does not have impact to the measurement residual, i.e., $r_{i}=0$.

In order to reduce numerical complexity, it is possible to avoid running GBP inference across the entire factor graph defined by (13). Instead, we can transform the original factor graph into a considerably more compact version which is equivalent from the restoration viewpoint. The variable nodes from sets $\mathcal{X}_{i}$, $i=1,\dots,k$ are connected by factor nodes $\mathcal{F}_{\mathrm{c}}$ and $\mathcal{F}_{\mathrm{p}}$ defined according to matrices $\mathbf{H}_{\mathrm{c}}$ and $\mathbf{H}_{\mathrm{p}}$, respectively. Also, each island $s_{i}$ has $l_{i}=n_{i}-1$ linearly independent equations, where revealing any of the state variables from the set $\mathcal{X}_{i}$ reveals the entire observable island $s_{i}$. This allows the construction of the factor graph, where all state variables $\mathcal{X}_{i}$ enclosed by the island $s_{i}$ are merged into a single variable node i.e., each island becomes the variable node. Thus, we observe the factor graph with a set of variable nodes $\mathcal{S}=\{s_{1},\dots,s_{k}\}$ and the set of factor nodes $\mathcal{F}_{\mathrm{c}}\cup\mathcal{F}_{\mathrm{p}}$. Same as before, we introduce the set of virtual factor nods $\mathcal{F}_{\text{s}}$ attached to the variable nodes from the set $\mathcal{S}$. This strategy reveals an exciting finding: the factor graph defined by the variable nodes $\mathcal{S}$, and factor nodes $\mathcal{F}_{\mathrm{c}}$ and $\mathcal{F}_{\mathrm{p}}$ equates to the matrix $\mathbf{W}_{\text{bc}}$ given in [10]. However, in contrast to [10] where the Gram matrix is needed, we are solving the observability restoration problem directly over the matrix $\mathbf{W}_{\text{bc}}$ using the GBP algorithm.

We formalise the proposed method for the observability restoration in Algorithm 2. The initialisation step (lines 1-13) creates the factor graph, where each factor node $f_{i}$ from the set $\mathcal{F}_{\mathrm{c}}$ has the same variance equal to $v_{i}$, i.e., variance $v_{i}$ corresponds to the equation associated with the factor node $f_{i}$. Each variable node $s_{i}\in\mathcal{S}$ has the virtual factor node with variances set to infinity (lines 8-10). The algorithm continues observability restoration procedure (lines 14-44) until $w=k-q-1$ linearly independent equations are found (line 11). The observability restoration procedure (lines 14-44) starts by selecting a factor node $f_{\nu}$ from the set $\mathcal{F}_{\mathrm{p}}$ followed by changing its variance from the infinity to $v_{i}$. In other words, we select one of the equations defined by pseudo-measurements that can merge islands (line 16). Then, we select one of the variable nodes incident to the factor node $f_{\nu}$ and change the variance of its virtual factor node to zero. This virtual factor node has a role of a slack variable node that exists in the system (lines 17-18). The message-passing procedure (lines 20-29) commence after the initial GBP variances from variable nodes $s_{i}\in{\mathcal{S}}$ to factor nodes $f_{j}\in\mathcal{F}_{\text{c}}\cup\mathcal{F}_{\text{p}}$ are set to finite positive values (line 19), where the set $\mathcal{F}_{i}\setminus f_{j}$ represents factor nodes incident to the variable node $s_{i}$, excluding the factor node $f_{j}$, while the set ${\mathcal{S}}_{j}\setminus s_{i}$ represents variable nodes incident to the factor node $f_{j}$, excluding the variable node $s_{i}$. The algorithm proceeds with the linear independence procedure (lines 30-42). The equation, described by the factor node $f_{\nu}$, is linearly dependent or independent, based on the variance value from factor node $f_{\nu}$ to the residual variable node $r_{\nu}$. The linear independence procedure marks an equation as linearly independent if the variance $v_{f_{\nu}\to r_{\nu}}$ has an infinite value (line 34). Thus, the number of required equations to provide observability is reduced by one (line 35). In contrast, the linearly dependent equation is excluded from the graph by changing its variance from $v_{i}$ to infinity (lines 37-38), and the algorithm proceeds by selecting a new factor node (line 16). The procedure is repeated until observability of the entire system is achieved.

Example 3** (Observability Restoration).**

We demonstrate the GBP-based observability restoration method by leveraging Example 2, where we identified 3 islands $\mathcal{S}=\{s_{1},s_{2},s_{3}\}$, as depicted in Figure 5.

We assume pseudo-measurements can merge islands only if injected at the boundary buses of observable islands. Thus, we consider pseudo-measurements $M_{P_{1}}$, $M_{P_{2}}$, and $M_{P_{4}}$, which define the set $\mathcal{F}_{\text{p}}=\{f_{12},f_{13},f_{14}\}$, respectively. The injection measurement $M_{P_{3}}$ which was rejected during observable island identification define the set $\mathcal{F}_{\text{c}}=\{f_{11}\}$ (Figure 6, red box), with finite positive value variance (e.g., equal to one). This measurement represents a linearly independent equation and it is added at the beginning by default.

The resulting factor graph is shown in Figure 6, where we first select the factor node $f_{12}$ (yellow box) and set its variance to a finite positive value, while the variances of factor nodes $f_{13}$ and $f_{14}$ (grey boxes) are set to infinity. We begin by setting variances from the virtual factor nodes as follows:

[TABLE]

The GBP iterations are performed in a similar way as in Example 2, except that here we use variance $v_{i}$ to compute variances from factor nodes $f_{i}\in\mathcal{F}_{\text{c}}\cup\mathcal{F}_{\text{p}}$ to variable nodes $s_{i}\in{\mathcal{S}}$. When the algorithm converges, we include the residual variable node $r_{12}$ incident to the factor node $f_{12}$ and compute variance $v_{f_{12}\to r_{12}}$ using variances into the factor node $f_{12}$:

[TABLE]

Thus, the factor node $f_{12}$ represents a linearly independent equation, and the system is fully observable using measurements $M_{P_{1}}$ and $M_{P_{3}}$ described by factor nodes $f_{11}$ and $f_{12}$. respectively.

IV Numerical Results

In this section, we assess our proposed method against the numerical and topological state-of-the-art methods available in the literature. In our evaluation, we explore a large set of power systems with different bus configurations. We evaluate the performance of the proposed method using power systems with 1354, 2000, 10000, 25000 and 70000 buses [31, 32]. For each power system, we generate 100000 random measurement configurations in order to get statistically significant results. The measurement set includes active power flow on both branch sides and active power injection measurements. For each power system size, we use different redundancies to control for the resulting number of observable islands, where we aim at producing approximately the same number of measurement configurations per each number of islands ranging between 2 and 33.

In all scenarios, we use box plots where we normalise execution time of each of the different state-of-the-art methods $t_{\text{sa}}$ by the execution time of the GBP-based approach $t_{\text{bp}}$. Hence, the case $t_{\text{sa}}/t_{\text{bp}}>1$ corresponds to the case where the GBP-based method shows better performance (i.e., lower execution time), while $t_{\text{sa}}/t_{\text{bp}}<1$ corresponds to the case where the GBP-based method shows worse performance (i.e., higher execution time) in comparison with the state-of-the-art method. The methods have been implemented using Julia programming language, while the source code of the GBP-based observability analysis can be found as the part of the JuliaGrid package222https://github.com/mcosovic/JuliaGrid.jl.

IV-A Determination of Observable Islands

To evaluate the GBP-based island detection algorithm, we compare the execution time of the GBP-based method $t_{\text{bp}}$ against the numerical method based on the nodal variable formulation $t_{\text{nm}}$ [12, Sec. 4.6.2] and the topological method based on the multi-stage procedure $t_{\text{tm}}$ [13].

Figure 7 shows distributions of the normalised execution times $t_{\text{nm}}/t_{\text{bp}}$ and $t_{\text{tm}}/t_{\text{bp}}$ for power systems with 1354 and 2000 buses. We split results into two groups corresponding to the particular range of the number of islands. For all experiments, the topological and GBP-based methods achieve significantly lower execution times compared to the numerical method. Both topological and the GBP-based methods result in comparable execution times across all experiments with the topological method achieving marginally lower execution times. Note that median values of the execution time for the topological and GBP-based methods are between $0.78\,\text{ms}$ and $1.37\,\text{ms}$, and between $1.38\,\text{ms}$ and $2.96\,\text{ms}$, respectively333The methods have been tested on a 64 bit Windows 10 with an Intel CoreTM i7-8850H 2.60 GHz CPU, and 32 GB of RAM. It is important to emphasise that the same GBP iteration scheme (i.e., stopping criterion threshold, the maximum number of iterations) was used throughout all scenarios. However, for the specific system, by tuning these parameters, execution times can be further improved..

Using a funnel-based approach, we omit numerical method from further experiments due their large execution times on power systems with a smaller number of buses, making it unpractical for larger systems. Figure 8 shows box plots of the normalised execution time $t_{\text{tm}}/t_{\text{bp}}$ for the power systems with 10000, 25000 and 70000 buses with regard to the number of islands. The GBP-based method outperforms topological method for larger power systems. Also, for the fixed power system size (i.e., the number of buses), increasing the number of islands results in lower execution times for the GBP-based method compared to the topological method. Note that the execution time medians of the GBP-based method for 70000-buses are $92.97\,\text{ms}$ and $123.43\,\text{ms}$ for 2-17 and 18-33 group of islands, respectively.

IV-B Observability Restoration

To evaluate the GBP-based observability restoration algorithm, we compare the execution time of the GBP-based method $t_{\text{bp}}$ against the numerical method based on the Gram matrix associated with a Jacobian matrix corresponding to a reduced network and measurement set $t_{\text{nm}}$ [10].

Figure 9 shows distributions of the normalised execution times $t_{\text{nm}}/t_{\text{bp}}$ for the power system with 70000 buses. We split results into four groups corresponding to the particular range of the number of islands. For all experiments, execution times for both methods are on order of few milliseconds. Still, the GBP-based method achieves lower execution times compared to the numerical method. Median values of the execution time for the GBP-based and numerical methods are between $1.04\,\text{ms}$ and $3.33\,\text{ms}$, and between $5.45\,\text{ms}$ and $5.66\,\text{ms}$, for 2-9 and 26-33 group of islands, respectively. As a sideline, among 100000 random measurement configurations, the numerical method failed to find non-redundant set of equations in $239$ cases due to the incorrect zero pivot threshold, while the GBP-based method succeeded throughout all measurement configurations. Note that the redundancy has a significant effect on the number of failed attempts to find the non-redundant set of equations for the numerical method, and increases as the redundancy decreases.

IV-C Computational complexity analysis

The complexity of the GBP algorithm, considering only the asymptotic behaviour, depends on: i) the scaling of the number of edges in the factor graph, and ii) the scaling of the number of iterations required for the algorithm to converge.

The first point exploits the underlying system sparsity (i.e., measurement factor nodes are typically connected to a very low number of state variable nodes), as the computational effort per iteration is proportional to the number of edges in the factor graph. For each of the $k$ measurements, the degree of the corresponding factor node is limited by a typically small constant. Indeed, for any type of measurements, the corresponding measurement function depends only on a few state variables. As the number of state variables $n$ and the number of measurements $k$ grows large, the number of edges in the factor graph scales as $\mathcal{O}(n)$. Therefore the computational complexity of the GBP scales linearly per iteration. The above analysis is directly related to the sparsity of the underlying factor graph [20, 28].

For the second point, the scaling of the number of the GBP iterations as $n$ grows large is a more challenging problem related to the convergence rate of the GBP algorithm. Based on a discussion in [33] for full matrices, we know that the number of iterations is likely to scale with the condition number of the underlying matrix, which, for well-conditioned matrices may scale as low as $\mathcal{O}(1)$ (constant). In our method, we only track the evolution of the variance values, because variances converge much faster than means. Tracking only variances results in a significant simplification. In our framework, we expect the number of the GBP iterations to scale as a constant.

Another advantage of using the GBP algorithm is the possibility of distributing GBP method across disjoint areas by arbitrarily segmenting the underlying factor graph. In the extreme case of the fully-distributed GBP algorithm, each factor graph node operates locally and independently. In this scenario, the problem is distributed across $\mathcal{O}(n)$ nodes, and, if implemented to run in parallel, it can be $\mathcal{O}(n)$ times faster than the centralised solution.

V Conclusions

We presented a novel computationally efficient GBP-based observability analysis suitable for large-scale data processing in electric power systems that can be easily implemented as part of a fully distributed architecture. In essence, proposed algorithms use only addition operations over the factor graph. Both, determination of observable islands and observability restoration algorithm, are implemented by tracking evolution of the variances, which results in the algorithm that always converges. While algorithms use floating point calculations, most of the time they operate as topological methods. Finally, for large-scale power systems, GBP-based method significantly outperforms numerical and topological state-of-the-art methods, reducing execution time of island detection procedure while being more robust in a search for correct non-redundant set of equations.

VI Acknowledgment

This paper has received funding from the European Union’s Horizon 2020 research and innovation programme under Grant Agreement number 856967.