Distributed State Estimation of Multi-region Power System based on Consensus Theory

Abstract

Effective state estimation is critical to the security operation of power systems. With the rapid expansion of interconnected power grids, there are limitations of conventional centralized state estimation methods in terms of heavy and unbalanced communication and computation burdens for the control center. To address these limitations, this paper presents a multi-area state estimation model and afterwards proposes a consensus theory based distributed state estimation solution method. Firstly, considering the nonlinearity of state estimation, the original power system is divided into several non-overlapped subsystems. Correspondingly, the Lagrange multiplier method is adopted to decouple the state estimation equations into a multi-area state estimation model. Secondly, a fully distributed state estimation method based on the consensus algorithm is designed to solve the proposed model. The solution method does not need a centralized coordination system operator, but only requires a simple communication network for exchanging the limited data of boundary state variables and consensus variables among adjacent regions, thus it is quite flexible in terms of communication and computation for state estimation. In the end, the proposed method is tested by the IEEE 14-bus system and the IEEE 118-bus system, and the simulation results verify that the proposed multi-area state estimation model and the distributed solution method are effective for the state estimation of multi-area interconnected power systems.

1. Introduction

Power system state estimation has an important impact on power system security assessment and on-line security control. At present, the centralized state estimation method is widely researched, and usually a control center is needed to estimate the state of the whole system based on comprehensive measurement information provided by the Supervisory Control and Data Acquisition (SCADA) system [1]. With the development of interconnected power systems and the expanded scale of power grids, a control center is used to collect huge measurement data for the state estimation, and thus this center usually faces the defects of heavy communication and computation burden. Therefore, it is significant to divide interconnected large power grids into several sub-areas and accordingly design a distributed state estimation models with efficient solution methods.

The distributed state estimation models with proper solution approaches have been investigated by some scholars [2,3,4,5]. A distributed state estimation model was established with a hierarchical framework in [6], where the local estimation results of the first layer need to be globally coordinated at the center of the second layer. However, the reliability of this hierarchical estimation model was slightly poor. In [7], based on the Richeardson iterative and linear equations, a distributed method was proposed for the weighted least square estimation. This method could well adapt to power network topology changes, but the convergence performance was not convincing enough. Authors of [8] presented an effective and accurate method for multi-area power system state estimation, which could improve the bad data detection capability and the convergence speed by using well-tuned measurement weights. In [9], an overlapped sub-region partitioning was proposed for the distributed state estimation problem, and a parallel algorithm based on the principle of virtual boundary measurement was designed to solve the problem. Authors of [10] proposed a multi-area paralleled state estimation method on the basis of sub-regions independently performing the local state estimation and the control center coordinately calculating the global optimal solution with the boundary variables. The distributed state estimation approach presented in [11] was based on a recursive estimation algorithm with time-varying weight parameters, which needed global communication network information to calculate the eigenvalues of Laplacian matrix and afterward determined the parameter range of the time-varying weights, and the number of convergence iterations was usually large. In [12], a state estimation algorithm was designed with optimal weighting factors for power system state estimation and the data packet loss was also considered in the model.

In recent years, popular applications of wide area measurement system (WAMS) and phasor measurement unit (PMU) in power systems have significantly promoted the development of state estimation [13,14]. Since sampling phase errors are inevitable among PMUs in their deployments, authors of [15] proposed a novel power state estimation algorithm with consideration to PMU phase errors. Authors of [16] presented a PMU-based decomposition method for distributed state estimation and evaluated the impact of measurement errors on state estimation accuracy. Authors of [17] divided the whole system into observable and unobservable sub-regions and solved them separately. With the help of WAMS/SCADA mixed measurements, a distributed state estimation algorithm was formed based on matrix-splitting strategy to improve the robustness of state estimation in [18]. In [19], a bi-linear state estimation model with mixed measurements was presented and afterward solved by the alternating direction multiplier method.

This paper proposes a nonlinear multi-area state estimation model and presents a consensus theory based distributed state estimation method. For the proposed method, each sub-system independently conducts the individual state estimation according to the local measurements, and only limited information of consensus variables and boundary-bus state variables is shared among adjacent regions, while the system-level global optimal solution can be obtained within several iterations. The main contributions of this paper are threefold.

(1) A multi-area state estimation nonlinear model is proposed in this paper by using Lagrange multiplier function, and a consensus theory based distributed state estimation method is designed accordingly to solve the model.

(2) The proposed distributed state estimation method does not need a coordination center, but only requires a simple communication network for exchanging limited data of boundary state variables and consensus variables among adjacent regions, thus it is quite flexible in terms of communication and computation for state estimation.

(3) Compared with the centralized method, simulation results of the IEEE 14-bus and the IEEE 118-bus systems demonstrate that the proposed multi-area distributed state estimation method is effective and flexible for state estimation of interconnected power systems.

The main contents of this paper are organized as follows: Section 2 introduces the traditional state estimation model. Section 3 establishes a multi-area state estimation model, while Section 4 covers the consensus theory based distributed state estimation method designed to solve this model. In Section 5, simulations of the IEEE 14-bus and the IEEE 118-bus systems conducted to validate the effectiveness of proposed distributed method are addressed, and conclusions are drawn in the last section.

2. Traditional State Estimation Model

This section generally describes the traditional state estimation model. Assume the measurement vector of power systems are z including the voltage magnitudes, power flows, power injections, etc. The system state vector is x representing the voltage angles and magnitudes of all buses. In a traditional state estimation model, the state variable x is described in measurement equation h(x), which could be used to establish the relationship of system state vector x and the measurement vector z. When further considering the measurement noise, the relationship can be refined as (1):

$$z=h(x)+v$$

(1)

where h(x) and v are the function vector and noise vector of measurement data, respectively.

Since the accuracy of each measurement is different, the weighted least square (WLS) method is widely used in state estimation problems in order to improve the accuracy of state estimation, and each measurement is tuned with a proper weight related to its variance. If the variance of the measurement error is denoted as a diagonal matrix R = diag[σ₁², σ₂²,…, σ_m²], where σ_i² is the variance of the measurement error v_i, the objective function of WLS method for state estimation problems is [20].

$$\min J(x)={(z-h(x))}^T{R}^{-1}(z-h(x))$$

(2)

If the initial value of the state variable is denoted as x⁽⁰⁾ and h(x) is expanded around x⁽⁰⁾ using the Taylor series, the function vector h(x) can be linearized at x⁽⁰⁾ as (3):

$$h(x)=h(x^{(0)})+H(x^{(0)})\Delta x$$

(3)

where H(x⁽⁰⁾) is the Jacobian matrix of the function vector h(x). By substituting (3) into the weighted least squares criterion, we can obtain J(x) = [∆z − H(x⁽⁰⁾)∆x]^TR⁻¹[∆z − H(x⁽⁰⁾)∆x], where ∆z=z − h(x⁽⁰⁾). To minimize the function J(x) in (2), the following equation must be satisfied:

$$\Delta \widehat{x}=G^{-1}(x^{(0)})H^T(x^{(0)})R^{-1}\Delta z$$

(4)

where G(x⁽⁰⁾) = H^T(x⁽⁰⁾)R^-1H(x⁽⁰⁾). With (4), the variable $\widehat{x}$ can be updated from the initial value x⁽⁰⁾ by (5):

$$\widehat{x}=x^{(0)}+\Delta \widehat{x}=x^{(0)}+G^{-1}(x^{(0)})H^T(x^{(0)})R^{-1}[z-h(x^{(0)})]$$

(5)

With the iteration number (indicated as k) increased, the following formula can be adopted to continuously update $\widehat{x}$ until achieving a satisfactory accuracy:

$$\begin{array}{cc}\Delta {\widehat{x}}^{(k)}& =G^{-1}({\widehat{x}}^{(k)})H^T({\widehat{x}}^{(k)})R^{-1}[z-h({\widehat{x}}^{(k)})]\\ & ={[H^T({\widehat{x}}^{(k)})R^{-1}H({\widehat{x}}^{(k)})]}^{-1}{H}^T({\widehat{x}}^{(k)})R^{-1}[z-h({\widehat{x}}^{(k)})]\\ & {\widehat{x}}^{(k+1)}={\widehat{x}}^{(k)}+\Delta {\widehat{x}}^{(k)}\end{array}$$

(6)

3. Multi-Area State Estimation Model

In this section, a nonlinear multi-area state estimation model is presented. Assume that a power network has n buses and is portioned into r non-overlapped sub-system S_i (i = 1, 2, ..., r, each sub-system S_i includes n_i buses) during the system planning phase based on the following principles: (1) each sub-system has a similar network scale with the consideration of geographical location; (2) designed with few tie lines between subsystems, which can reduce the amount of data exchange and communication burden among subsystems for distributed state estimation.

Figure 1 shows a non-overlapped sub-area partition diagram, and the measurements include voltage amplitude, nodal power injections (indicated by the small solid square), and line power flow (indicated by the small filled-dot). In the multi-area state estimation model, the buses with no connection to any neighboring regions are called internal buses, and the buses connected to neighboring regions are called boundary buses. Correspondingly, the measurements related to boundary buses are called boundary measurements, as circled by the large red dot circle, and the measurements inside the sub-areas are called internal measurements, as circled by the large blue dot circle.

Assume that the slack bus of the whole network is included in sub-system i, and its voltage amplitude is measurable. The other sub-area has a local reference bus, and its voltage and phase angle are set as pseudo-measurements. It is also assumed that all sub-areas are observable. Based on these assumptions and the non-overlapped partition scheme in Figure 1, the state estimation measurement equation of the system can be divided into the internal and boundary equations:

$$\begin{array}{l}{z}_i=h_i(x_i)+v_i,i=1,2\cdots r\\ z_c=h_c(x)+v_c\end{array}$$

(7)

where z_i is the internal measurement for sub-system i, and it is a column vector with dimensions m_i × 1, while m_i is the number of internal measurements; x_i is a 2n_i × 1 column vector of voltage magnitudes and phase angles for internal buses; z_c is a m_c × 1 column vector of boundary measurements where m_c is the number of boundary measurements; x = [x₁, x₂,…,x_r]^T is the state variable of the original system; v_i and v_c are the corresponding measurement errors for internal and boundary buses. The proposed multi-area state estimation can be expressed in terms of the WLS minimization problem with boundary equality constraints as (8):

$$\begin{array}{ll}\hfill \min J(x)& =r_c^T{R}_c^{-1}{r}_c+{\displaystyle \sum _{i=1}^r{r}_i^T(x_i)R_i^{-1}r}(x_i)\\ \hfill r_c& =z_c-h_c(x)\end{array}$$

(8)

where the measurement error r_i (x_i) = z_i − h_i (x_i), the internal measurement variance matrix R_i = cov(v_i) = diag[σ₁², σ₂²… σ_mi²], and the boundary measurement variance matrix R_c = cov(v_c) = diag[σ₁², σ₂²…σ_mc²].

By using Lagrange multiplier method, the state estimation problem (8) can be transformed into the Lagrange function L(x, λ) as (9):

$$L(x,\lambda )=\frac{1}{2}{r}_c^T{R}_{\mathrm{c}}^{-1}{r}_c+\frac{1}{2}{\displaystyle \sum _{i=1}^r{r}_i^T(x_i)R_i^{-1}{r}_i(x_i)}+{\lambda }^T(r_c-z_c+h_c(x))$$

(9)

where λ is a m_c × 1 vector of Lagrange multiplier. To minimize L(x, λ) in (9), the following Equations (10)–(12) should be satisfied at the optimal solution:

$$\frac{\partial L}{\partial x}=0\Rightarrow -\left(\begin{array}{ccc}{H}_1^T({\widehat{x}}_1)R_1^{-1}& & \\ & \ddots & \\ & & H_r^T({\widehat{x}}_r)R_r^{-1}\end{array}\right)\times \left(\begin{array}{c}{z}_1-h_1({\widehat{x}}_1)\\ \cdots \\ z_r-h_r({\widehat{x}}_r)\end{array}\right)+H_c^T(\widehat{x})\widehat{\lambda }=0$$

(10)

$$\frac{\partial L}{\partial r_c}=0\Rightarrow R_c^{-1}{r}_c+\widehat{\lambda }=0$$

(11)

$$\frac{\partial L}{\partial \lambda }=0\Rightarrow r_c-z_c+h_c(\widehat{x})=0$$

(12)

From Equations (11) and (12), we have:

$$h_c(\widehat{x})-z_c-R_c\widehat{\lambda }=0$$

(13)

The transformed nonlinear Equations (10) and (13) can be solved by Newton’s iterative method as follows:

$$\left(\begin{array}{ccc}{G}_1(x_1^k)& & \\ & \ddots & \\ & & G_r(x_r^k)\end{array}\right)\times \left(\begin{array}{c}\Delta x_1^k\\ \cdots \\ \Delta x_r^k\end{array}\right)+H_c^T(x^k){\lambda }^{k+1}=\left(\begin{array}{c}{H}_1^T(x_1^k)R_1^{-1}\Delta z_1^k\\ \cdots \\ H_r^T(x_r^k)R_r^{-1}\Delta z_r^k\end{array}\right)$$

(14)

$$H_c(x^k)\Delta x^k-R_c{\lambda }^{k+1}=\Delta z_c^k$$

(15)

where k is the number of iterations, and ∆x^k = [∆x₁^k, ∆x₂^k,…∆x_r^k]^T = [(x₁^k+1 − x₁^k), (x₂^k+1 − x₂^k),…,(x_r^k+1 − x_r^k)]^T; ∆z_i^k = z_i − h_i(x_i^k), I = 1,2…r; ∆z_c^k = z_c − h_c(x^k); G_i^k(x_i^k) = H_i^T(x_i^k)R_i^-1H_i(x_i^k); H_c(x^k) = [H_c1(x₁^k), H_c2(x₂^k),…, H_cr(x_r^k)], and H_ci is the Jacobian matrix (m_c×2n_i) of state variables for the boundary measurement in sub-region S_i.

The Equations of (14) and (15) can be iteratively solved by (16)–(19) with the help of introducing an intermediate variable $\Delta y_i^k$ in Equation (16) [10]:

$$\Delta y_i^k=G_i^{-1}(x_i^k)H_i^T(x_i^k)R_i^{-1}\Delta z_i^k,i=1,2\cdots r$$

(16)

$${\lambda }^{k+1}=G_c^{-1}(x^k)(\Delta z_c^k-{\displaystyle \sum _{i=1}^r{H}_{ci}(x_i^k)}\Delta y_i^k)$$

(17)

$$u_i^{k+1}=G_i^{-1}(x_i^k)H_{ci}^T(x_i^k){\lambda }^{k+1},i=1,2\cdots r$$

(18)

$$\Delta x_i^k=\Delta y_i^k+u_i^{k+1},i=1,2\cdots r$$

(19)

where G_c is a m_c order square matrix calculated by (20):

$$G_c(x^k)={\displaystyle \sum _{i=1}^r(R_{ci}+H_{ci}(x_i^k)G_i^{-1}(x_i^k)H_{ci}^T(x_i^k))}$$

(20)

As can be observed above, though Equations (16) and (19) can be solved independently for each sub-area i, the calculation of G_c(x^k), H_ci(x_i^k) and ∆z_c^k in Equations (17) and (18) depends on the boundary measurements and state variables of boundary buses. Therefore, a control center is needed for the centralized method to collect the necessary information of these boundary bus measurements and boundary bus state variables from the whole system for calculating λ^k+1 and u_i^k+1. However, with the increased scale of the power grid, the communication and calculation burden of the control center for centralized state estimation are very likely to be overloaded, which may jeopardize the reliability. In order to balance computational burden and improve reliability, a fully distributed state estimation method based on consensus theory is proposed in the next section.

4. Consensus Algorithm Based Distributed State Estimation Approach

4.1. Graph Description and Consensus Algorithm

4.1.1. Graph Description

A graph G can be used to model the communication network topology of power networks, which is a pair of sets (V_G, E_G), where V_G is a vertex set containing all vertices of graph G, and E_G is an edge set of unordered binary arrays of all elements in V_G. If there is a path between any two different vertices in graph G, the graph is called a connected graph, and the structure of the connected graph can be reflected by an adjacency matrix A. When G is a finite and simple undirected graph containing n vertices, the adjacency matrix A = [a_ij] of G is a n × n symmetric matrix. The element a_ij is a (0,1)-matrix with zero diagonal [21]. If there is a connection between vertice i and vertice j, a_ij = 1, otherwise a_ij = 0. The row-stochastic matrix D = [d_ij] could be designed based on the adjacency matrix A, and its element d_ij, representing the communication coefficient between agents i and j, could be defined by (21) [22,23]:

$$d_{ij}=\{\begin{array}{cc}1/(b_i+b_j+2)& j\in N_i\\ 1-{\displaystyle {\sum }_{j\in N_i}1/(b_i+b_j+2)}& j\notin N_i\\ 0& otherwise\end{array}$$

(21)

where N_i represents the index of adjacent vertices of vertex i, and b_i is the sum of the number of edges associated with vertex i, defined as:

$$b_i={\displaystyle \sum _{j=1}^n{a}_{ij}}$$

(22)

4.1.2. Consensus Algorithm

Denote the state of vertex i by w_i. According to the consensus theory, if and only if all vertex states in the network are equal, i.e., w₁ = w₂ =…w_n, the vertices of the network reach consensus. A consensus algorithm is given in [24] as:

$${\dot{w}}_i=-{\displaystyle \sum _{j=1}^n{a}_{ij}(w_i-w_j)} i=1,2\cdots n$$

(23)

where a_ij is an element of (i, j) in the adjacency matrix A. When it takes a fixed time for information share between vertices, we need to model the consensus network as a discrete-time system using d_ij and D, which is described as follows [23]:

$$w_i\left[t+1\right]={\displaystyle \sum _{j=1}^n{d}_{ij}{w}_j\left[t\right]} i=1,2\cdots n$$

(24)

$$w\left[t+1\right]=Dw\left[t\right]$$

(25)

The stopping criterion for the consensus procedures is that the value w_i[t] of each node i satisfies $\left|w_i\left[t\right]-w_j\left[t\right]\right|\le \epsilon , \forall j\in N_i$, where the parameter ε is a preset threshold for error tolerance [25].

4.2. Propose Distributed State Estimation Method

For the multi-area state estimation model proposed in Section 3, the key point of realizing the distributed method is to solve λ^k+1 of (17) in a distributed manner. Specifically, understanding how to effectively calculate the two terms G_c(x^k) and ∆z_c^k − ∑H_ci(x_i^k)∆y_i^k) of (17) in a distributed way is essential for the proposed distributed state estimation method, which is elaborated upon in the following section.

Equation (17) can be rewritten in the form of (26):

$${\lambda }^{k+1}={\left(\frac{1}{N}{G}_{\mathrm{c}}(x^k)\right)}^{-1}\left(\frac{1}{N}(\Delta z_c^k-{\displaystyle \sum _{i=1}^r{H}_{ci}(x_i^k)}\Delta y_i^k)\right)$$

(26)

where the error ∆z_c^k and ∆z_ci^k are defined as:

$$\Delta z_c^k={\left[\begin{array}{c}{\left(z_{c1}-h_{c1}(x^k)\right)}_{m_{c1}\times 1}\\ ⋮\\ {\left(z_{ci}-h_{ci}(x^k)\right)}_{m_{ci}\times 1}\\ ⋮\\ {\left(z_{c\mathrm{r}}-h_{cr}(x^k)\right)}_{m_{cr}\times 1}\end{array}\right]}_{m_c\times 1}\begin{array}{c}Sub-system 1\\ ⋮\\ Sub-system i\\ ⋮\\ Sub-system r\end{array}$$

(27)

$$\Delta z_{ci}^k={\left[\begin{array}{c}{\left(\begin{array}{c}0\\ ⋮\\ 0\end{array}\right)}_{m_{c1}\times 1}\\ ⋮\\ {\left(z_{ci}-h_{ci}(x^k)\right)}_{m_{ci}\times 1}\\ ⋮\\ {\left(\begin{array}{c}0\\ ⋮\\ 0\end{array}\right)}_{m_{cr}\times 1}\end{array}\right]}_{m_c\times 1}\begin{array}{c}Sub-system 1\\ ⋮\\ ⋮\\ Sub-system i\\ ⋮\\ ⋮\\ Sub-system r\end{array}$$

(28)

where m_ci is the number of boundary measurements for region S_i, and h_ci(•) is the function vector for the boundary measurements. Thus, ∆z_c^k can be converted into the following form:

$$\Delta z_c^k={\displaystyle \sum _{i=1}^r\Delta z_{ci}^k}$$

(29)

Afterward, Equation (26) can be converted to:

$${\lambda }^{k+1}={\left(\frac{1}{N}{G}_{\mathrm{c}}(x^k)\right)}^{-1}\left(\frac{1}{N}{\displaystyle \sum _{i=1}^r(\Delta z_{ci}^k-H_{ci}(x_i^k)}\Delta y_i^k)\right)$$

(30)

4.2.1. Distributed Solution Process

According to the consensus theory, the average value of consensus variables can be iteratively approached via limited iterations of exchanging local data. Therefore, there is no need to set up a control center to collect the global information from all sub-areas for computing G_c(x^k) and ∑(∆z_ci^k − H_ci(x_i^k)∆y_i^k) of (30) in a centralized way. Instead, the consensus algorithm can be used in each sub-system to solve λ in a distributed manner. In the following section, we demonstrate how to conveniently calculate the terms G_c(x^k) and ∑(∆z_ci^k − H_ci(x_i^k)∆y_i^k) of (30) in a distributed manner.

(1) For distributed calculation of G_c(x^k)

Apply the consensus algorithm to solve Equation (30) for the multi-area state estimation. The term (1/N)G_c(x^k) [refer to Equation (20)] can be calculated in a distributed manner as follows:

$${\left[\begin{array}{c}{R}_{c1}+H_{c1}(x_1^k)G_1^{-1}(x_1^k)H_{c1}^T(x_1^k)\\ R_{c2}+H_{c2}(x_2^k)G_2^{-1}(x_2^k)H_{c2}^T(x_2^k)\\ ⋮\\ R_{cr}+H_{cr}(x_r^k)G_r^{-1}(x_r^k)H_{cr}^T(x_r^k)\end{array}\right]}^{\left[t+1\right]}=D\cdot {\left[\begin{array}{c}{R}_{c1}+H_{c1}(x_1^k)G_1^{-1}(x_1^k)H_{c1}^T(x_1^k)\\ R_{c2}+H_{c2}(x_2^k)G_2^{-1}(x_2^k)H_{c2}^T(x_2^k)\\ ⋮\\ R_{cr}+H_{cr}(x_r^k)G_r^{-1}(x_r^k)H_{cr}^T(x_r^k)\end{array}\right]}^{\left[t\right]}$$

(31)

$$\left[\begin{array}{c}\frac{1}{N}{G}_c(x^k)\\ \frac{1}{N}{G}_c(x^k)\\ ⋮\\ \frac{1}{N}{G}_c(x^k)\end{array}\right]=\underset{t\to \infty }{\lim }{\left[\begin{array}{c}{R}_{c1}+H_{c1}(x_1^k)G_1^{-1}(x_1^k)H_{c1}^T(x_1^k)\\ R_{c2}+H_{c2}(x_2^k)G_2^{-1}(x_2^k)H_{c2}^T(x_2^k)\\ ⋮\\ R_{cr}+H_{cr}(x_r^k)G_r^{-1}(x_r^k)H_{cr}^T(x_r^k)\end{array}\right]}^{\left[t+1\right]}=\underset{t\to \infty }{\lim }{D}^t\cdot {\left[\begin{array}{c}{R}_{c1}+H_{c1}(x_1^k)G_1^{-1}(x_1^k)H_{c1}^T(x_1^k)\\ R_{c2}+H_{c2}(x_2^k)G_2^{-1}(x_2^k)H_{c2}^T(x_2^k)\\ ⋮\\ R_{cr}+H_{cr}(x_r^k)G_r^{-1}(x_r^k)H_{cr}^T(x_r^k)\end{array}\right]}^{\left[0\right]}$$

(32)

where D is the row-stochastic matrix calculated by Equation (21) based on the communication network topology, k is the number of external iteration to update the state variables, and t is the number of internal iteration to update the consensus variables. Theoretically, when t approaches infinity, each sub-area can get the exact mean value of consensus variables in Equation (32). In practice, very few iterations are needed to get the approximated mean value of consensus variables with satisfied engineering accuracy.

Based on Equation (32), each sub-area i independently calculates the local consensus variable R_ci + H_ci(x_i^k)G_i⁻¹(x_i^k)H_ci^T(x_i^k) and then interacts with the neighboring sub-areas through the row-stochastic matrix D, and afterwards obtains the mean value of network variables via few iterations. In this way, the term (1/N)G_c(x^k) can be calculated in a distributed manner.

(2) For distributed calculation of ∑(∆z_ci^k − H_ci(x_i^k)∆y_i^k)

Similarly, the second term (1/N)∑(∆z_ci^k − H_ci(x_i^k)∆y_i^k) in (30) can also be solved in a distributed manner by (33) and (34) based on consensus algorithm:

$${\left[\begin{array}{c}\Delta z_{c1}^k-H_{c1}(x_1^k)\Delta y_1^k\\ \Delta z_{c2}^k-H_{c2}(x_2^k)\Delta y_2^k\\ ⋮\\ \Delta z_{cr}^k-H_{cr}(x_r^k)\Delta y_r^k\end{array}\right]}^{\left[t+1\right]}=D\cdot {\left[\begin{array}{c}\Delta z_{c1}^k-H_{c1}(x_1^k)\Delta y_1^k\\ \Delta z_{c2}^k-H_{c2}(x_2^k)\Delta y_2^k\\ ⋮\\ \Delta z_{cr}^k-H_{cr}(x_r^k)\Delta y_r^k\end{array}\right]}^{\left[t\right]}$$

(33)

$$\left[\begin{array}{c}\frac{1}{N}{\displaystyle \sum _{i=1}^r(\Delta z_{ci}^k-H_{ci}(x_i^k)}\Delta y_i^k)\\ \frac{1}{N}{\displaystyle \sum _{i=1}^r(\Delta z_{ci}^k-H_{ci}(x_i^k)}\Delta y_i^k)\\ ⋮\\ \frac{1}{N}{\displaystyle \sum _{i=1}^r(\Delta z_{ci}^k-H_{ci}(x_i^k)}\Delta y_i^k)\end{array}\right]=\underset{t\to \infty }{\lim }{\left[\begin{array}{c}\Delta z_{c1}^k-H_{c1}(x_1^k)\Delta y_1^k\\ \Delta z_{c2}^k-H_{c2}(x_2^k)\Delta y_2^k\\ ⋮\\ \Delta z_{cr}^k-H_{cr}(x_r^k)\Delta y_r^k\end{array}\right]}^{\left[t+1\right]}=\underset{t\to \infty }{\lim }{D}^t\cdot {\left[\begin{array}{c}\Delta z_{c1}^k-H_{c1}(x_1^k)\Delta y_1^k\\ \Delta z_{c2}^k-H_{c2}(x_2^k)\Delta y_2^k\\ ⋮\\ \Delta z_{cr}^k-H_{cr}(x_r^k)\Delta y_r^k\end{array}\right]}^{\left[0\right]}$$

(34)

With the help of (32) and (34), the Lagrange multiplier λ in (30) or (17) can be obtained in a distributed manner. Considering that G_i⁻¹(x_i^k) and H_ci^T(x_i^k) could be calculated in each sub-system independently, u_i^k of (18) can be computed in a distributed manner with λ obtained by (30). Therefore, the proposed state estimation solution approach essentially including Equations (16)–(19) can be completely calculated in a distributed manner.

Remarks

The essential step of realizing a distributed state estimation algorithm is to calculate (in a distributed manner) the global variable Lagrange multiplier λ in (17), which is actually coupled with the global state variables x^k in terms of G_c(x^k) and ∆z_c^k − ∑H_ci(x_i^k)∆y_i^k. The proposed distributed algorithm transforms (17) into (30), and then introduces the consensus theory to compute the terms G_c(x^k) and ∑(∆z_ci^k − H_ci(x_i^k)∆y_i^k) by Equations (32) and (34) in a distributed manner. For the proposed algorithm, there is no need to set up a control center to collect the measurements of the whole system. Each sub-system only requires the local measurements in its own sub–area and exchanges consensus variables with the adjacent regions to calculate the average value of consensus variables, thereby realizing a distributed state estimation approach.

4.2.2. Main Steps of Proposed Algorithm

The detailed steps of the proposed multi-area distributed state estimation method are as follows:

Step 1 Initialize k = 0 and set initial values for state variables as well as the internal iteration convergence precision ɛ and external iteration convergence precision δ;

Step 2 Collect network structure parameters, and form an adjacency matrix A according to the communication network topology among sub-systems;

Step 3 The following procedures are performed simultaneously in each sub-system S_i:

Step 3.1 Obtain internal measurement z_i and boundary measurement z_ci, as well as R_i (the internal measurement variance matrix) and R_ci (the boundary measurement variance matrix in sub-system i);

Step 3.2 Calculate the internal measurement function h_i(x^k), internal Jacobi matrix H_i(x^k), and G_i(x^k) = H_i^T(x^k)R_i⁻¹H_i(x^k) for each subsystem

Step 3.3 Calculate internal error ∆z_i^k=z_i − h_i(x_i^k) then compute ∆y_i^k according to Equation (16);

Step 4 Each sub-area interacts with neighboring areas to compute h_ci(x^k) and H_ci(x^k) and calculate ∆z_ci^k = z_ci − h_ci(x^k)

Step 5 Each sub-area exchanges consensus variables G_ci(x^k) and ∆z_ci^k − H_ci(x_i^k)∆y_i^k, iteratively calculates (1/N)G_c(x^k) and (1/N)∑(∆z_ci^k − H_ci(x_i^k)∆y_i^k) according to Equations (32) and (34) until the internal convergence precision ɛ is satisfied, and then computes λ^k+1 according to Equation (30);

Step 6 Each sub-area calculates u_i^k according to Equation (18) and then obtains ∆x_i^k and x_i^k+1 according to Equations (19) and (6);

Step 7 Check the convergence of each subsystem. If all subsystems meet the convergence conditions (∆x_i^k < δ), x^k is the global optimal solution; otherwise k = k + 1, return to step 3.2.

5. Simulation Results

5.1. Case 1: IEEE 14-bus System

The effectiveness of the proposed method is verified by the simulation results of the IEEE 14-bus system. The internal state variable convergence precision ɛ is settled as 10⁻⁴, and the external state variables convergence precision δ is settled as 10⁻⁵ for the following simulations. The system partition and measurement configuration are shown in Figure 2. Set bus 1 as the slack bus with zero phase angle and measurable voltage amplitude. The other measurements include the active and reactive power flow and the active and reactive power injections, as shown in Figure 2. The total number of measurements for the whole system is 43. The communication topology of sub-systems is shown in Figure 3. The proposed distributed state estimation method has no special requirements for communication networks, and the mature communication technologies, such as carrier communication and optical fiber communication, are all applicable.

In this case study, the measured data are assumed as the solutions of power flow calculation plus Gaussian white noise with standard deviation σ, where voltage amplitude measurement error σ_v = 0.03, bus injection power measurement error σ_Pi = 0.01, and line transmission power error σ_Pij = 0.004. The boundary measurement may belong to either of the two adjacent areas according to the actual measurement position. In addition to the reference bus located in sub-system i, each of the other sub-systems is provided with a zero phase angle pseudo measurement and 1 p.u. voltage pseudo measurement {in this example [(θ₃, V₃), (θ₆, V₆), (θ₉, V₉)]}. The detailed system measurements are as follows:

$$\begin{gathered} z_1=\left(\begin{array}{l}0.0000\\ 1.5708\\ 0.7594\\ -0.4081\\ 1.0600\\ -0.1748\\ 0.0532\\ -0.0193\end{array}\right)|\begin{array}{l}{\theta }_1\\ P_{1\text{-}2}\\ P_{1\text{-}5}\\ P_{5\text{-}2}\\ V_1\\ Q_{1\text{-}2}\\ Q_{1\text{-}5}\\ Q_{5\text{-}2}\end{array}{z}_2=\left(\begin{array}{l}0.0000\\ -0.2510\\ 0.2707\\ 0\\ 1.0000\\ 0.0473\\ -0.1540\\ 0.2103\end{array}\right)|\begin{array}{l}{\theta }_3\\ P_{3\text{-}4}\\ P_{4\text{-}7}\\ P_8\\ V_3\\ Q_{3\text{-}4}\\ Q_{4\text{-}7}\\ Q_8\end{array}{z}_3=\left(\begin{array}{l}0.0000\\ -0.0816\\ 0.1834\\ 0.0188\\ -0.0610\\ 1.0000\\ -0.0864\\ 0.0998\\ 0.0141\\ -0.0160\end{array}\right)|\begin{array}{l}{\theta }_6\\ P_{11\text{-}6}\\ P_{6\text{-}13}\\ P_{12\text{-}13}\\ P_{12}\\ V_6\\ Q_{11\text{-}6}\\ Q_{6\text{-}13}\\ Q_{12\text{-}13}\\ Q_{12}\end{array}{z}_4=\left(\begin{array}{l}0.0000\\ 0.0745\\ -0.1051\\ 1.0000\\ 0.0665\\ -0.0366\end{array}\right)|\begin{array}{l}{\theta }_9\\ P_{9\text{-}10}\\ P_{14\text{-}9}\\ V_9\\ Q_{9\text{-}10}\\ Q_{14\text{-}9}\end{array} \\ {z}_{c1}=\left(\begin{array}{l}0.1830\\ 0.6006\\ 0.3523\\ -0.1006\end{array}\right)|\begin{array}{l}{P}_2\\ P_{5\text{-}4}\\ Q_2\\ Q_{5\text{-}4}\end{array}{z}_{c2}=\left(\begin{array}{l}-0.9420\\ -0.5427\\ 0.2707\\ 0.0000\\ 0.0876\\ 0.0213\\ 0.1480\\ 1.0000\end{array}\right)|\begin{array}{l}{P}_3\\ P_{4\text{-}2}\\ P_{7\text{-}9}\\ {\theta }_3\\ Q_3\\ Q_{4\text{-}2}\\ Q_{7\text{-}9}\\ V_3\end{array}{z}_{c3}=\left(\begin{array}{l}-0.0350\\ 0.4589\\ 0.0000\\ -0.0180\\ -0.2084\\ 1.0000\end{array}\right)|\begin{array}{l}{P}_{11}\\ P_{5\text{-}6}\\ {\theta }_6\\ Q_{11}\\ Q_{5\text{-}6}\\ V_6\end{array}{z}_{c4}=\left(\begin{array}{l}-0.1490\\ 0.1546\\ 0.0000\\ -0.0500\\ -0.0264\\ 1.0000\end{array}\right)|\begin{array}{l}{P}_{14}\\ P_{4\text{-}9}\\ {\theta }_9\\ Q_{14}\\ Q_{4\text{-}9}\\ V_9\end{array} \end{gathered}$$

Based on the above measurement configurations, the proposed multi-area distributed estimation method is applied to the IEEE 14-bus system. The solutions of proposed distributed estimation method are benchmarked with that of the conventional centralized weighted least square state estimation method, and these comparisons are shown in

Table 1. Solution comparisons of distributed and centralized state estimation method for the IEEE 14-bus system.

State Variable	Centralized Method (°)	Distributed Method (°)	State Variable	Centralized Method (p.u.)	Distributed Method (p.u.)
θ1	0	0	V1	1.0600	1.0605
θ2	−4.9808	−4.9816	V2	1.0450	1.0446
θ3	−12.7176	−12.7860	V3	1.0100	1.0080
θ4	−10.3241	−10.2464	V4	1.0186	1.0140
θ5	−8.7825	−8.7684	V5	1.0203	1.0174
θ6	−14.2223	−14.4500	V6	1.0700	1.0683
θ7	−13.3680	−13.2353	V7	1.0620	1.0480
θ8	−13.3680	−13.2353	V8	1.0900	1.0822
θ9	−14.9462	−14.8114	V9	1.0563	1.0333
θ10	−15.1039	−15.0375	V10	1.0513	1.0256
θ11	−14.7949	−14.8601	V11	1.0571	1.0435
θ12	−15.0771	−15.3064	V12	1.0569	1.0512
θ13	−15.1586	−15.3386	V13	1.0504	1.0444
θ14	−16.0386	−16.0839	V14	1.0358	1.0166

.

As shown in Table 1, the estimation results in terms of voltage amplitudes and phase angles solved by the proposed distributed state estimation method are very similar to those of the centralized state estimation method with the maximal voltage amplitude error 2.4% for the bus 10 and the maximal phase angle error 1.5% for the bus 12. From the vivid bar chart based solution comparisons of two methods in Figure 4, it is clear that the estimation results of the proposed distributed method are very similar to those of the centralized method, and therefore the proposed distributed state estimation method has a good accuracy.

In order to verify the convergence of the proposed method, Figure 5 shows the solution errors of the distributed method benchmarked with the centralized state estimation method in terms of nodal voltage amplitudes and phase angles at each external iteration. For the proposed distributed algorithm, one step of the external iteration is with several internal consensus iteration steps, and it is performed only after the internal iteration loop converges. Therefore, the convergence of algorithm errors can be reflected with regard to the external iterations in Figure 5. It can be seen that the voltage amplitude and phase angle differences of every node decrease rapidly when the number of iterations increases, and finally, they can converge to a satisfactory accuracy via only seven or eight iterations. These curves indicate that the proposed multi-area distributed state estimation method has good convergence performance.

It should also be noted that the adjacency matrix A varies for different communication network topologies. However, a different matrix A only affects the number of internal iterations for consensus variables convergence in Equation (17), while the solution accuracy of the distributed state estimation method is not affected. Based on the investigations of the influence of communication network topologies on convergence speed in [15,21], increasing the communication network connection density enhances the information exchange efficiency and thus accelerates the convergence speed.

5.2. Case 2: IEEE 118-bus System

In order to check the effectiveness and scalability of the distributed state estimation algorithm, an IEEE 118-bus system is adopted to test the proposed method. The IEEE 118-bus system including 179 lines is partitioned into six subsystems, as shown in Figure 6 (black lines represent internal lines in sub-systems and red lines represent boundary lines) with the communication network topology shown in Figure 7. For each sub-area, a local phase angle reference bus is settled as: S₁: {bus 1}, S₂: {bus 17}, S₃: {bus 33}, S₄: {bus 24}, S₅: {bus 74}, S₆: {bus 91}. The total number of measurements is 389 (except for pseudo measurements). The standard deviation σ of Gaussian white noise for various measurement data is the same as that in Case 1.

With the above measurement configurations, the traditional centralized weighted least square state estimation method and the proposed distributed state estimation method are both used to solve the state estimation problem of the IEEE 118-bus system, and Figure 8 clearly demonstrates the solution differences between these two methods at each iteration. It can be observed from Figure 8 that the results of the proposed distributed method gradually approach the state estimation results of the centralized method, and finally, the former matches well with the solutions of the centralized method after 10–12 iterations, which is demonstrated by the magnification part in Figure 8. Furthermore, when comparing Figure 5 for the IEEE 14-bus system with Figure 8 for the IEEE 118-bus system, it can be observed that with the increasing number of system buses, the error of the distributed state estimation results is not deteriorated for a large power network. These comparisons verify that the proposed distributed state estimation method is effective and scalable for the state estimation of large scale power systems.

6. Conclusions

In this paper, a multi-area state estimation model is proposed by using Lagrange multiplier algorithm, and then a novel consensus theory based distributed state estimation method is designed to solve the model. The proposed distributed state estimation method does not need a coordination center and only requires a simple communication network for exchanging the limited data of boundary state variables and consensus variables among adjacent regions. Compared with the conventional centralized method, the proposed distributed state estimation method can solve state estimation problems with satisfactory accuracy and good convergence. Simulation results of the IEEE 14-bus and the IEEE 118-bus systems demonstrate that the proposed distributed method is effective and flexible for state estimation of multi-area interconnected power systems.