A Comparative Study on Power Flow Methods Applied to AC Distribution Networks with Single-Phase Representation

Abstract

This paper presents a comparative analysis of six different iterative power flow methods applied to AC distribution networks, which have been recently reported in the scientific literature. These power flow methods are (i) successive approximations, (ii) matricial backward/forward method, (iii) triangular-based approach, (iv) product linearization method, (v) hyperbolic linearization method, and (vi) diagonal approximation method. The first three methods and the last one are formulated without recurring derivatives, and they can be directly formulated in the complex domain; the fourth and fifth methods are based on the linear approximation of the power balance equations that are also formulated in the complex domain. The numerical comparison involves three main aspects: the convergence rate, processing time, and the number of iterations calculated using the classical Newton–Raphson method as the reference case. Numerical results from two test feeders composed of 34 and 85 nodes demonstrate that the derivative-free methods have linear convergence, and the methods that use derivatives in their formulation have quadratic convergence.

1. Introduction

The estimation of the state variables in electrical systems corresponds to an essential study in electrical engineering [1], since the variables, i.e., the magnitudes and angles of the voltages in all the nodes of the network, allow us to determine the operative conditions of the grid [2], such as active and reactive power losses, grid efficiency, conductors’ chargeability, etc. The determination of these variables is widely-known in the literature as the power flow problem [3]. The mathematical representation and the solution of the power flow problem in power systems are typically developed using the Newton–Raphson method or its derived methods due to its efficiency and robustness regarding the number of iterations and the quadratic convergence rate [4]. In the case of electrical distribution grids, because of their configurations (radial or weakly meshed), power flow approaches that exploit the grid configuration to reduce the processing time required in the power solution are recommended [5]. For distribution systems, multiple power flow reformulations based on the graph structure of the grid can be found in the literature [6]. The most classical methods correspond to the backward/forward method and its derivations.

Table 1. Power flow methods for distribution networks.

Solution Methodology	Year	Type	Refs.
Power flow solution using multi-port compensation technique for radial and meshed grids	1988	DF	[7]
Backward/forward power flow method considering voltage-controlled nodes	1995	DF	[8]
Current injection power flow method for radial and meshed distribution networks	1996	DF	[9]
Improved Gauss–Seidel power flow method for distribution systems	2002	DF	[10]
Direct power flow solution using LU matrix decomposition	2003	DF	[11]
Improved Newton–Raphson method with Broyden’s method for distribution grids	2008	DB	[12]
Improved Newton–Raphson method with Levenberg–Marquardt method for distribution grids	2008	DB	[13]
Backward/forward power flow solution for radial and weakly meshed distribution grids	2010	DF	[14]
Fast decoupled power flow method for emerging distribution grids	2010	DB	[15]
Triangular formulation based on a real quasi-symmetry matrix	2013	DF	[16]
Fast decoupled power flow method for distribution grids	2015	DB	[17]
Axis rotation fast decoupled load flow on distribution systems	2016	DB	[18]
Linear power flow approximation with hyperbolic linearization	2016	DB	[19]
Linear power flow approximation based on the admittance matrix	2016	DB	[20]
Graph-based power flow using an incidence matrix	2018	DF	[21]
Graph-based power flow using an upper-triangular matrix	2019	DF	[22]
Successive approximations power flow method that guarantees convergence	2020	DF	[23]
Hyperbolic recursive linearization power flow method	2020	DB	[24]
Matricial backward/forward power flow method that guarantees convergence and includes voltage-controlled nodes	2020	DF	[25]
Triangular-based power flow method that guarantees convergence	2021	DF	[26,27]
Product linearization power flow method	2021	DB	[28]
Linearized power flow approach for transmission and distribution networks	2021	DB	[29]

DB: Derivative-based; DF: Derivative-free.

summarizes some references concerning power flow solutions in distribution grids in the last decades.

The power approaches listed in Table 1 can be classified in two main groups. The first group includes the power flow approaches that incorporate the distribution system graph by proposing recursive solutions without having recurring derivatives in their formulations [21,30], implying that these formulations only rearrange the power flow equations to obtain a recursive formula [23]. The second group of methodologies are based on Taylor’s series expansion of the power flow equations in real and complex domains that generate iterative solution methods based on linear approximations [24]. In the case of derivative-free methods, their convergence rate is typically linear, while the derivative-based methods have a quadratic convergence rate [31,32].

Based on the aforementioned revision of the state of the art regarding power flow methods applied to electrical AC distribution networks, in this research, we present a comparative study of the recently developed methodologies (from 2018 to 2021) that include four derivative methods and two derivative-based methods. To evaluate the effectiveness of these methods, the conventional Newton–Raphson method is considered the reference approach in two IEEE test feeders composed of 34 and 85 nodes, respectively, both with radial structure [33]. The numerical assessment considers the average processing time, the number of iterations, and the convergence rate analysis. The main contribution of this research is the summarization of the main recently developed power flow methods for AC distribution networks in a unique work, whereby researchers (also readers) can have multiple options to solve the power flow problem by selecting either derivative-free or derivative-based methods.

The main contributions of this research are as follows:

✓

A complete comparison between recently developed power flow methods comprising four derivative-free and two derivative-based methods. These comparisons present convergence properties that are linear for derivative-free methods and quadratic for derivative-based approaches.

✓

A new iterative method for power flow solution from the family of non-derivative approaches that present linear convergence and use the information of the nodal admittance matrix to obtain its recursive formula.

It is worth mentioning that, in the current literature, there exist some additional power flow approaches based on convex reformulations [34,35], such as semidefinite programming [36,37] and second-order cone programming [38,39]; however, these convex reformulations are typically solved with interior point methods since they are used in optimal power flow analyses, i.e., optimization problems wherein an objective performance (technical, economical, or environmental performance, or a combination of them) is considered. In this regard, these methods are outside the scope of this study since the main objective is to provide a comparison of the numerical methods for power flow solution.

This paper is structured as follows: Section 2 describes the general formulation of the power flow problem, and its subsections present each one of the studied methodologies, with all of them formulated in the complex domain. Section 3 mentions the main characteristics of the test feeders, which include two IEEE radial test feeders having 34 and 85 nodes and all the electrical information to implement the studied power flow methods. Section 4 presents the numerical validation of the studied methods using multiple consecutive evaluations that help determine the average processing time, number of iterations, and convergence rate properties. Section 5 discusses the main conclusions drawn from this research as well as possible future works.

2. Power Flow Formulations

This section presents the mathematical derivation of six power flow methodologies, which are applicable to single-phase distribution networks with radial structure. The methods under study are (i) successive approximations [23], (ii) matricial backward/forward method [21], (iii) triangular-based approach [22], (iv) product linearization method [24], (v) hyperbolic linearization method [28], and (vi) diagonal approximation approach. The main advantage of these power flow approaches is that they are formulated in the complex domain, thereby making their implementation simple in programming environments such as MATLAB or Python.

2.1. Formulation of the General Power Flow Problem

The power flow problem arises due to constant power loads in electrical power systems, which cause a hyperbolic relation between voltages and currents [19]. Let us define the net injected current in an arbitrary node as a function of the nodal admittance matrix and the voltage profile in all nodes of the network [1]:

$$\begin{array}{c}\hfill {\mathbb{I}}_k=\sum _{m\in \mathcal{N}}{\mathbb{Y}}_{km}{\mathbb{V}}_m,\forall k\in \mathcal{N}⟺\mathbb{I}={\mathbb{Y}}_{\mathrm{bus}}\mathbb{V},\end{array}$$

(1)

where ${\mathbb{I}}_k$ is the complex current injected at node k, ${\mathbb{Y}}_{km}$ is the complex component of the nodal admittance matrix that connects nodes k and m, ${\mathbb{V}}_m$ is the complex value of the voltage at node m, $\mathbb{I}$ is the vector that contains all the nodal currents, ${\mathbb{Y}}_{\mathrm{bus}}$ is the square complex nodal admittance matrix, and $\mathbb{V}$ is the vector that contains all the nodal voltages. Note that $\mathcal{N}$ is the set that includes all the nodes of the network.

Now, if we refer to the definition of the average complex power, i.e., Tellegen’s theorem [40], we find the general definition of the power flow problem as follows:

$$\begin{array}{c}{\mathbb{S}}_k^{\star }={\mathbb{V}}_k^{\star }\sum _{m=1}^n{\mathbb{Y}}_{km}{\mathbb{V}}_m,\end{array}$$

(2)

where ${\mathbb{S}}_k^{\star }$ represents the conjugate of the net injected complex power at node k. Note that Equation (2) can be expressed in a compacted form as follows [19]:

$$\begin{array}{c}\hfill {\mathbb{S}}^{\star }=\mathbf{diag}\left({\mathbb{V}}^{\star }\right){\mathbb{Y}}_{\mathrm{bus}}\mathbb{V}.\end{array}$$

(3)

2.2. Successive Approximation Power Flow Method

The successive approximation power flow method is a recursive power flow formulation of AC distribution networks initially proposed in [23], which deals with radial and meshed distribution networks. The main characteristic of this power flow method is its derivative-free formulation, which does not require the inverse of non-diagonal matrices at each iteration [6]. To express this power flow method, let us rearrange the power flow formula (3) as follows:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\mathbb{S}}_s^{\star }\\ -{\mathbb{S}}_d^{\star }\end{array}\right]=\left[\begin{array}{cc}\mathbf{diag}\left({\mathbb{V}}_s^{\star }\right)& 0\\ 0& \mathbf{diag}\left({\mathbb{V}}_d^{\star }\right)\end{array}\right]\left[\begin{array}{cc}{\mathbb{Y}}_{ss}& {\mathbb{Y}}_{sd}\\ {\mathbb{Y}}_{ds}& {\mathbb{Y}}_{dd}\end{array}\right]\left[\begin{array}{c}{\mathbb{V}}_s\\ {\mathbb{V}}_d\end{array}\right].\end{array}$$

(4)

where ${\mathbb{S}}_s^{\star }$ and ${\mathbb{S}}_d^{\star }$ represent the complex conjugate injected power in the slack node and the complex conjugate demanded power in the load nodes, respectively; ${\mathbb{V}}_s^{\star }$ and ${\mathbb{V}}_s^{\star }$ are the complex conjugate voltages in the slack and demand nodes, respectively; ${\mathbb{Y}}_{ss}$, ${\mathbb{Y}}_{sd}$, ${\mathbb{Y}}_{ds}$, and ${\mathbb{Y}}_{dd}$ are the sub-components of the admittance matrix that connect the slack and demand nodes among them.

Note that the first row in (4) is linear since the voltage in the slack node is perfectly known; this implies that if we know all the voltage values in the demand nodes, the power generation in the slack can be determined without performing any iterative process. In addition, note that the vector that contains all values of constant power consumption is preceded by a minus sign, which indicates that the demands in the nodes leave the distribution system.

From considering the second row of (4), one can observe a nonlinear relation among voltage demands, which implies that it is necessary to solve this formula using a numerical method. For doing so, let us obtain a recursive formula for solving ${\mathcal{V}}_d$ from the second row of (4), which produces the following result.

$$\begin{array}{c}\hfill {\mathbb{V}}_d=-{\mathbb{Y}}_{dd}^{-1}\left[{\mathbb{Y}}_{ds}{\mathbb{V}}_s+{\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{\star }\right){\mathbb{S}}_d^{\star }\right].\end{array}$$

(5)

To solve the recursive formula (5), it is necessary to add an iterative counter t to determine the final voltage value by starting from an initial point, which is typically assigned as the slack voltage, i.e., ${\mathbb{V}}_d^t=\mathbf{1}{\mathbb{V}}_s$, with $\mathbf{1}$ a vector constituting by ones with appropriate dimensions. For simplicity, we reassign the inverse of the admittance matrix ${\mathbb{Y}}_{dd}^{-1}$ as ${\mathbb{Z}}_{dd}$, which is calculated once and stored for speeding up the iterative process, computationally speaking. The final recursive formula for the successive approximation method is the following.

$$\begin{array}{c}\hfill {\mathbb{V}}_d^{t+1}=-{\mathbb{Z}}_{dd}\left[{\mathbb{Y}}_{ds}{\mathbb{V}}_s+{\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{t,\star }\right){\mathbb{S}}_d^{\star }\right].\end{array}$$

(6)

The iterative process presented in (6) ends when the convergence criteria are met, which is selected as the difference of the voltage magnitudes between two consecutive iterations [23], i.e.,

$$\begin{array}{c}\hfill max\left\{\left|\left|{\mathbb{V}}_d^{t+1}\right|-\left|{\mathbb{V}}_d^t\right|\right|\right\}\le \epsilon ,\end{array}$$

(7)

which is tested at each iteration. Note that $\epsilon$ is the tolerance error, assigned as $1\times {10}^{-10}$ [6].

Remark 1.

The convergence test for the successive approximation power flow method was provided in [23] using the Banach fixed-point theorem. This proves that the voltage profile always converges to the final solution independent of the starting point if and only if the operative point of the system is located far away from the voltage collapse point.

2.3. Matricial Backward/Forward Power Flow Method

The matricial backward/forward power flow method involves the generalization of the iterative sweep method for distribution grids using a nodal incidence matrix [21]. This reformulation method was initially proposed by the authors of [21,25]. To formulate this power flow method, let us define the branch-to-node incidence matrix as follows:

Branch-to-node incidence matrix: The matrix $\mathbf{A}\in {\mathcal{R}}^{b\times n}$ is the branch-to-node incidence matrix that represents an electrical distribution network composed of b branches and n nodes if their currents are selected with arbitrary senses, and

✓

${\mathbf{A}}_{jk}=1$, if the current of the line j is leaving from the node k.

✓

${\mathbf{A}}_{jk}=-1$, if the current of the line j is arriving at the node k.

✓

${\mathbf{A}}_{jk}=0$, if the line j is not connected to the node k.

With the incidence matrix, it is possible to link the voltages in all the branches, i.e., ${\mathcal{E}}_b$, with the voltage nodes by applying Kirchhoff’s Second Law to each branch, as follows:

$$\begin{array}{c}\hfill {\mathbb{E}}_b={\mathbf{A}}_s{\mathbb{V}}_s+{\mathbf{A}}_d{\mathbb{V}}_d,\end{array}$$

(8)

where ${\mathbf{A}}_s$ is a column vector of the incidence matrix associated with the slack node, and ${\mathbf{A}}_d$ is the remainder columns of the incidence matrix associated with the demand nodes.

Now, when Kirchhoff’s First Law is applied to each node of the network, assuming that all the demanded currents are leaving the system, i.e., they are negative, the following relation between the branch currents ${\mathcal{I}}_b$ and the demanded currents is the following:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\mathbb{I}}_s\\ -{\mathbb{I}}_d\end{array}\right]=\left[\begin{array}{c}{\mathbf{A}}_s^{\top }\\ {\mathbf{A}}_d^{\top }\end{array}\right]{\mathbb{I}}_b.\end{array}$$

(9)

An additional important relation between the branch voltages and currents are defined by Ohm’s law, i.e., ${\mathbb{I}}_b={\mathbb{Y}}_{bb}{\mathbb{E}}_b$, where ${\mathbb{Y}}_{bb}$ is the primitive admittance matrix that contains in its diagonal the inverse of the line impedance of each line. Using this relation in (8) and considering the second row of (9), the following formula is obtained.

$$\begin{array}{c}\hfill -{\mathbb{I}}_d={\mathbf{A}}_d^{\top }{\mathbb{Y}}_{bb}\left[{\mathbf{A}}_s{\mathbb{V}}_s+{\mathbf{A}}_d{\mathbb{V}}_d\right].\end{array}$$

(10)

In (10), it is possible to replace the nonlinear relation between demand voltages and powers, i.e., ${\mathbb{I}}_d={\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{\star }\right){\mathbb{S}}_d^{\star }$, and obtain a formula to solve ${\mathbb{V}}_d$ as follows:

$$\begin{array}{c}\hfill {\mathbb{V}}_d=-{\left[{\mathbf{A}}_d^{\top }{\mathbb{Y}}_{bb}{\mathbf{A}}_d\right]}^{-1}\left[{\mathbf{A}}_d^{\top }{\mathbb{Y}}_{bb}{\mathbf{A}}_s{\mathbb{V}}_s+{\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{\star }\right){\mathbb{S}}_d^{\star }\right].\end{array}$$

(11)

It is worth mentioning that if we redefine ${\mathbb{Y}}_{dd}$ as ${\mathbf{A}}_d^{\top }{\mathbb{Y}}_{bb}{\mathbf{A}}_d$ and ${\mathbb{Y}}_{ds}$ as ${\mathbf{A}}_d^{\top }{\mathbb{Y}}_{bb}{\mathbf{A}}_s$, the recursive formula (10) associated with the matricial backward/forward is equivalent to the successive approximation power flow method as demonstrated by Briñez et al. in [6]. Finally, the iterative formula to solve the power flow problem using the matricial backward/forward approach takes the following form.

$$\begin{array}{c}\hfill {\mathbb{V}}_d^{t+1}=-{\left[{\mathbf{A}}_d^{\top }{\mathbb{Y}}_{bb}{\mathbf{A}}_d\right]}^{-1}\left[{\mathbf{A}}_d^{\top }{\mathbb{Y}}_{bb}{\mathbf{A}}_s{\mathbb{V}}_s+{\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{t,\star }\right){\mathbb{S}}_d^{\star }\right].\end{array}$$

(12)

Remark 2.

In the case of the matricial backward/forward power flow method as demonstrated in [21,41] and [25], its convergence independent of the starting voltage point through the application of the Banach-fixed point theorem.

2.4. Triangular-Based Power Flow Method

The triangular-based power flow method is a graph-based method incorporating an upper-triangular matrix that connects the demand and the branches’ currents. This method was initially proposed by [22]. To formulate the triangular-based power flow method, let us define the general structure of the upper-triangular matrix.

Upper-triangular matrix [42]: The matrix $\mathbf{T}\in {\mathcal{R}}^{b\times (n-1)}$ is an upper-triangular matrix that represents an electrical distribution network composed of b branches and n nodes, which has the following structure:

✓

${\mathbf{T}}_{jk}=1$, if the current of the line j support the current consumption at node k.

✓

${\mathbf{T}}_{jk}=0$, if the current of the line j does not support the current consumption at node k.

Now, if we apply Kirchhoff’s First Law at each node and rewrite the current consumption at each line as a function of the demanded currents, the following relation is obtained.

$$\begin{array}{c}\hfill {\mathbb{I}}_b=\mathbf{T}{\mathbb{I}}_d.\end{array}$$

(13)

Note that in (13), Kirchhoff’s First Law was only applied to the demand nodes, and the slack nodes were not analyzed.

On the other hand, applying Kirchhoff’s Second Law to the closed-loop trajectory containing each demand node and the slack node yields the following results [27].

$$\begin{array}{c}\hfill {\mathbb{V}}_d=\mathbf{1}{\mathbb{V}}_s-{\mathbf{T}}^{\top }{\mathbb{E}}_b.\end{array}$$

(14)

In Equation (14), if we apply the Ohm’s law to each line voltage drop, i.e., ${\mathbb{E}}_b={\mathbb{Z}}_{bb}{\mathbb{I}}_b$, and also it is considered (13) and the hyperbolic relation between voltages and powers, we get the following result.

$$\begin{array}{c}\hfill {\mathbb{V}}_d=\mathbf{1}{\mathbb{V}}_s-{\mathbf{T}}^{\top }{\mathbb{Z}}_{bb}\mathbf{T}{\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{\star }\right){\mathbb{S}}_d^{\star }.\end{array}$$

(15)

Once again, to determine the voltage profile at the demand nodes, we add an iterative counter t into the recursive formula (15), which produces the following iterative equation.

$$\begin{array}{c}\hfill {\mathbb{V}}_d^{t+1}=\mathbf{1}{\mathbb{V}}_s-{\mathbf{T}}^{\top }{\mathbb{Z}}_{bb}\mathbf{T}{\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{t,\star }\right){\mathbb{S}}_d^{\star }.\end{array}$$

(16)

Remark 3.

To complete the iterative evaluation, the same detention criterion defined in (7) is applied. In addition, like the previously mentioned power flow methods, the triangular-based approach ensures convergence to the solution by applying the Banach fixed-point theorem [26].

2.5. Power Flow Approach Based on Voltage Product Linearization

The main characteristic of this power flow approach is that the Taylor series expansion is applied to obtain the equivalent linearized function that defines the product of two continuous variables around a desired operative point. This approach was recently proposed by Montoya et al. in [28]. This method linearizes the product among voltages in the second row of (5), i.e.,

$$\begin{array}{c}\hfill \mathbf{diag}\left({\mathbb{V}}_d^{\star }\right){\mathbb{Y}}_{ds}{\mathbb{V}}_s+\mathbf{diag}\left({\mathbb{V}}_d^{\star }\right){\mathbb{Y}}_{dd}{\mathbb{V}}_d+{\mathbb{S}}_d^{\star }=0,\end{array}$$

(17)

where the main interest of linearization is regarding the product $\mathbf{diag}\left({\mathbb{V}}_d^{\star }\right){\mathbb{Y}}_{dd}{\mathbb{V}}_d$, which can be represented mathematically as follows [28]:

$$\begin{array}{c}\hfill \mathbf{diag}\left({\mathbb{V}}_d^{\star }\right){\mathbb{Y}}_{dd}{\mathbb{V}}_d=\mathbf{diag}\left({\mathbb{V}}_d^{0,\star }\right){\mathbb{Y}}_{dd}{\mathbb{V}}_d+\mathbf{diag}\left({\mathbb{V}}_d^{\star }\right){\mathbb{Y}}_{dd}{\mathbb{V}}_d^0-\mathbf{diag}\left({\mathbb{V}}_d^{0,\star }\right){\mathbb{Y}}_{dd}{\mathbb{V}}_d^0,\end{array}$$

(18)

Now, if we substitute (18) in (17) and rearrange some terms, we get the following result:

$$\begin{array}{c}\hfill \mathbf{diag}\left({\mathbb{V}}_d^{\star }\right)\left[{\mathbb{Y}}_{ds}{\mathbb{V}}_s+{\mathbb{Y}}_{dd}{\mathbb{V}}_d^0\right]+\mathbf{diag}\left({\mathbb{V}}_d^{0,\star }\right){\mathbb{Y}}_{dd}{\mathbb{V}}_d+{\mathbb{S}}_d^{\star }-\mathbf{diag}\left({\mathbb{V}}_d^{0,\star }\right){\mathbb{Y}}_{dd}{\mathbb{V}}_d^0=0.\end{array}$$

(19)

In Equation (19), it is possible to use the properties of diagonal matrices, which allows rewriting the product $\mathbf{diag}\left({\mathbb{V}}_d^{\star }\right)\left[{\mathbb{Y}}_{ds}{\mathbb{V}}_s+{\mathbb{Y}}_{dd}{\mathbb{V}}_d^0\right]$ as $\left[\mathbf{diag}\left({\mathbb{Y}}_{ds}{\mathbb{V}}_s\right)+\mathbf{diag}\left({\mathbb{Y}}_{dd}{\mathbb{V}}_d^0\right)\right]{\mathbb{V}}_d^{\star }$; this implies that Equation (19) takes the following form:

$$\begin{array}{c}\hfill \mathbb{A}{\mathbb{V}}_d^{\star }+\mathbb{B}{\mathbb{V}}_d=\mathbb{C},\end{array}$$

(20)

where

$$\begin{array}{c}\hfill \mathbb{A}=\mathbf{diag}\left({\mathbb{Y}}_{ds}{\mathbb{V}}_s\right)+\mathbf{diag}\left({\mathbb{Y}}_{dd}{\mathbb{V}}_d^0\right),\\ \hfill \mathbb{B}=\mathbf{diag}\left({\mathbb{V}}_d^{0,\star }\right){\mathbb{Y}}_{dd},\\ \hfill \mathbb{C}=\mathbf{diag}\left({\mathbb{V}}_d^{0,\star }\right){\mathbb{Y}}_{dd}{\mathbb{V}}_d^0-{\mathbb{S}}_d^{\star }.\end{array}$$

Note that the solution of (20) requires to write it, into its real and imaginary components, i.e., ${\mathbb{V}}_d={\mathbf{V}}_{dr}+j{\mathbf{V}}_{di}$, $\mathbb{A}={\mathbf{A}}_r+j{\mathbf{A}}_i$, $\mathbb{B}={\mathbf{B}}_r+j{\mathbf{B}}_i$, and $\mathbb{C}={\mathbf{C}}_r+j{\mathbf{C}}_i$, which produce the following matricial relation:

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\mathbf{A}}_r+{\mathbf{B}}_r& {\mathbf{A}}_i-{\mathbf{B}}_i\\ {\mathbf{A}}_i+{\mathbf{B}}_i& {\mathbf{B}}_r-{\mathbf{A}}_r\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}_r\\ {\mathbf{V}}_i\end{array}\right]=\left[\begin{array}{c}{\mathbf{C}}_r\\ {\mathbf{C}}_i\end{array}\right].\end{array}$$

(21)

Note that the solution of Equation (21) is found by inverting the matrix that depends on the real and imaginary parts of the matrices $\mathbb{B}$ and $\mathbb{B}$, which are also dependent on the linearizing point, i.e., ${\mathbb{V}}_d^0$, in general, ${\mathbb{V}}_d^t$, implying that the general solution of (21) takes the following form:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\mathbf{V}}_r^{t+1}\\ {\mathbf{V}}_i^{t+1}\end{array}\right]={\left[\begin{array}{cc}{\mathbf{A}}_r^t+{\mathbf{B}}_r^t& {\mathbf{A}}_i^t-{\mathbf{B}}_i^t\\ {\mathbf{A}}_i^t+{\mathbf{B}}_i^t& {\mathbf{B}}_r^t-{\mathbf{A}}_r^t\end{array}\right]}^{-1}\left[\begin{array}{c}{\mathbf{C}}_r^t\\ {\mathbf{C}}_i^t\end{array}\right].\end{array}$$

(22)

Remark 4.

The iterative process in (22) ends when the tolerance criterion is met as defined in Equation (7). In addition, note that the convergence criteria of the power flow method associated with the linearization of the voltage product has not yet been addressed in the scientific literature, which can be considered as a research opportunity in future works.

2.6. Power Flow Approach Based on Hyperbolic Voltage Relation Linearization

This power flow methodology is based on the linearization of the hyperbolic relation between voltages and powers in the power balance equation. This approach was initially proposed by Garces et. al in [19] and subsequently improved by Bocanegra et. al in [24]. This power formulation rewrites the power flow equation in the demand nodes, i.e., (17), as follows:

$$\begin{array}{c}\hfill {\mathbb{Y}}_{ds}{\mathbb{V}}_s+{\mathbb{Y}}_{dd}{\mathbb{V}}_d+{\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{\star }\right){\mathbb{S}}_d^{\star }=0,\end{array}$$

(23)

where the main idea is to linearize the third component of the left-hand side, i.e., ${\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{\star }\right)$ ${\mathbb{S}}_d^{\star }$. This linearization is achieved through the application of Taylor series expansion, which produces, as presented in [24], the following result:

$$\begin{array}{c}\hfill {\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{\star }\right){\mathbb{S}}_d^{\star }=2{\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{0,\star }\right){\mathbb{S}}_d^{\star }-{\mathbf{diag}}^{-2}\left({\mathbb{V}}_d^{0,\star }\right)\mathbf{diag}\left({\mathbb{V}}_d^{\star }\right){\mathbb{S}}_d^{\star }.\end{array}$$

(24)

Now, if we substitute (24) in (23) and use the properties of the diagonal matrices, the following result is reached.

$$\begin{array}{c}\hfill {\mathbf{diag}}^{-2}\left({\mathbb{V}}_d^{0,\star }\right)\mathbf{diag}\left({\mathbb{S}}_d^{\star }\right){\mathbb{V}}_d^{\star }-{\mathbb{Y}}_{dd}{\mathbb{V}}_d-2{\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{0,\star }\right){\mathbb{S}}_d^{\star }-{\mathbb{Y}}_{ds}{\mathbb{V}}_s=0.\end{array}$$

(25)

Note that the solution of (25) can be reached with the same representation of (20) if the matrices $\mathbb{A}$ and $\mathbb{B}$ and the vector $\mathbb{C}$ take the following structure:

$$\begin{array}{c}\hfill \mathbb{A}={\mathbf{diag}}^{-2}\left({\mathbb{V}}_d^{0,\star }\right)\mathbf{diag}\left({\mathbb{S}}_d^{\star }\right),\\ \hfill \mathbb{B}=-{\mathbb{Y}}_{dd},\\ \hfill \mathbb{C}=2{\mathbf{diag}}^{-1}\left({\mathbb{V}}_d^{0,\star }\right){\mathbb{S}}_d^{\star }+{\mathbb{Y}}_{ds}{\mathbb{V}}_s.\end{array}$$

Remark 5.

The iterative solution of (25) is reached with the recursive formula (22) using the aforementioned matrices $\mathbb{A}$ and $\mathbb{B}$ and the vector $\mathbb{C}$. Its iterative search ends when the stopping criterion is met as defined in (7). Note that this method has not been reported in the literature regarding the convergence test, which offers another research opportunity to be addressed in future works.

2.7. Diagonal Approximation Power Flow Approach

A new iterative approach for solving power flow problem in electrical distribution networks is proposed in this section. This power flow formulation is a derivative-free approach that is created by rearranging Equation (17). Let us organize Equation (17) as follows:

$$\begin{array}{c}\hfill {\mathbb{V}}_d=-{\left[\mathbf{diag}\left({\mathbb{V}}_d^{\star }\right){\mathbb{Y}}_{dd}\right]}^{-1}\left[{\mathbb{S}}_d^{\star }+\mathbf{diag}\left({\mathbb{Y}}_{ds}{\mathbb{V}}_s\right){\mathbb{V}}_d^{\star }\right],\end{array}$$

(26)

Note that Equation (26) is recursive, which implies that it is necessary to add an iterative counter to solve it. The iterative formula for the diagonal approximation power flow approach takes the following form:

$$\begin{array}{c}\hfill {\mathbb{V}}_d^{t+1}=-{\left[\mathbf{diag}\left({\mathbb{V}}_d^{t,\star }\right){\mathbb{Y}}_{dd}\right]}^{-1}\left[{\mathbb{S}}_d^{\star }+\mathbf{diag}\left({\mathbb{Y}}_{ds}{\mathbb{V}}_s\right){\mathbb{V}}_d^{t,\star }\right].\end{array}$$

(27)

As observed in the previous cases, the detention criterion applicable to the recursive formula corresponds to the difference between both consecutive iterations. However, it is important to mention that the inverse of $\mathbf{diag}\left({\mathbb{V}}_d^{t,\star }\right){\mathbb{Y}}_{dd}$ exists as ${\mathbb{Y}}_{dd}$ is a dominant diagonal and $\mathbf{diag}\left({\mathbb{V}}_d^{t,\star }\right)$ is a diagonal whose component is always non-zero.

Remark 6.

The main characteristic of the diagonal approximation method is its similarity to the successive approximation power flow method (see Equation (6)) since both includes the inverse of the demand-to-demand admittance matrix; however, the diagonal approximation method considers the variations of voltage profiles in its numerical performance, which can affect the total processing time as well as the number of iterations.

3. Test Feeders

The assessment of each aforementioned power flow approach is carried out in two test feeders composed of 34 and 85 nodes, respectively. The main characteristics of the test feeders are mentioned below.

3.1. IEEE 34-Bus System

The IEEE 34-bus system is a radial distribution network comprising 34 nodes and 33 lines, which is operated at the substation bus in medium voltage levels with a nominal voltage rate of 11 kV. The electrical configuration of this test feeder is depicted in Figure 1. This test feeder has a total power consumption of $4636.50+j2873.50$ kVA, which produces a total active and reactive power loss of $221.75$ kW and $65.12$ kvar, respectively.

The parametric information concerning this test system is reported in

Table 2. Parametric information regarding the IEEE 34-bus system.

k	m	${\mathit{R}}_{\mathit{km}}$ ($\mathbf{\Omega }$)	${\mathit{x}}_{\mathit{km}}$ ($\mathbf{\Omega }$)	${\mathit{P}}_{\mathit{k}}$ (kW)	${\mathit{Q}}_{\mathit{k}}$ (kW)	k	m	${\mathit{R}}_{\mathit{km}}$ ($\mathbf{\Omega }$)	${\mathit{x}}_{\mathit{km}}$ ($\mathbf{\Omega }$)	${\mathit{P}}_{\mathit{k}}$ (kW)	${\mathit{Q}}_{\mathit{k}}$ (kW)
1	2	0.1170	0.0480	230	142.5	18	19	0.2079	0.0473	230	142.5
2	3	0.1073	0.0440	0	0	19	20	0.1890	0.0430	230	142.5
3	4	0.1645	0.0457	230	142.5	20	21	0.1890	0.0430	230	142.5
4	5	0.1495	0.0415	230	142.5	21	22	0.2620	0.0450	230	142.5
5	6	0.1495	0.0415	0	0	22	23	0.2620	0.0450	230	142.5
6	7	0.3144	0.0540	0	0	23	24	0.3144	0.0540	230	142.5
7	8	0.2096	0.0360	230	142.5	24	25	0.2096	0.0360	230	142.5
8	9	0.3144	0.0540	230	142.5	25	26	0.1310	0.0225	230	142.5
9	10	0.2096	0.0360	0	0	26	27	0.1048	0.0180	137	85
10	11	0.1310	0.0225	230	142.5	7	28	0.1572	0.0270	75	48
11	12	0.1048	0.0180	137	84	28	29	0.1572	0.0270	75	48
3	13	0.1572	0.0270	72	45	29	30	0.1572	0.0270	75	48
13	14	0.2096	0.0360	72	45	10	31	0.1572	0.0270	57	34.5
14	15	0.1048	0.0180	72	45	31	32	0.2096	0.0360	57	34.5
15	16	0.0524	0.0090	13.5	7.5	32	33	0.1572	0.0270	57	34.5
6	17	0.1794	0.0498	230	142.5	33	34	0.1048	0.0180	57	34.5
17	18	0.1645	0.0457	230	142.5	—	—	—	—	—	—

. Note that for all numerical simulations, we consider the voltage and power bases for this test system as 11 kV and 1000 kVA, respectively.

3.2. IEEE 85-Bus System

The IEEE 85-bus system has a radial medium-voltage distribution network composed of 85 nodes and 85 lines, operated with 11 kV. The total active and reactive power demand for this system is $2570.28+j2622.20$ kVA. The electrical configuration of this system is provided in Figure 2, and all its parametric information has been taken from [43], which is given in

Table 3. Parametric information regarding the IEEE 85-bus system.

k	m	${\mathit{R}}_{\mathit{km}}$ ($\mathbf{\Omega }$)	${\mathit{x}}_{\mathit{km}}$ ($\mathbf{\Omega }$)	${\mathit{P}}_{\mathit{k}}$ (kW)	${\mathit{Q}}_{\mathit{k}}$ (kW)	k	m	${\mathit{R}}_{\mathit{km}}$ ($\mathbf{\Omega }$)	${\mathit{x}}_{\mathit{km}}$ ($\mathbf{\Omega }$)	${\mathit{P}}_{\mathit{k}}$ (kW)	${\mathit{Q}}_{\mathit{k}}$ (kW)
1	2	0.108	0.075	0	0	34	44	1.002	0.416	35.28	35.99
2	3	0.163	0.112	0	0	44	45	0.911	0.378	35.28	35.99
3	4	0.217	0.149	56	57.13	45	46	0.911	0.378	35.28	35.99
4	5	0.108	0.074	0	0	46	47	0.546	0.226	14	14.28
5	6	0.435	0.298	35.28	35.99	35	48	0.637	0.264	0	0
6	7	0.272	0.186	0	0	48	49	0.182	0.075	0	0
7	8	1.197	0.820	35.28	35.99	49	50	0.364	0.151	36.28	37.01
8	9	0.108	0.074	0	0	50	51	0.455	0.189	56	57.13
9	10	0.598	0.410	0	0	48	52	1.366	0.567	0	0
10	11	0.544	0.373	56	57.13	52	53	0.455	0.189	35.28	35.99
11	12	0.544	0.373	0	0	53	54	0.546	0.226	56	57.13
12	13	0.598	0.410	0	0	52	55	0.546	0.226	56	57.13
13	14	0.272	0.186	35.28	35.99	49	56	0.546	0.226	14	14.28
14	15	0.326	0.223	35.28	35.99	9	57	0.273	0.113	56	57.13
2	16	0.728	0.302	35.28	35.99	57	58	0.819	0.340	0	0
3	17	0.455	0.189	112	114.26	58	59	0.182	0.075	56	57.13
5	18	0.820	0.340	56	57.13	58	60	0.546	0.226	56	57.13
18	19	0.637	0.264	56	57.13	60	61	0.728	0.302	56	57.13
19	20	0.455	0.189	35.28	35.99	61	62	1.002	0.415	56	57.13
20	21	0.819	0.340	35.28	35.99	60	63	0.182	0.075	14	14.28
21	22	1.548	0.642	35.28	35.99	63	64	0.728	0.302	0	0
19	23	0.182	0.075	56	57.13	64	65	0.182	0.075	0	0
7	24	0.910	0.378	35.28	35.99	65	66	0.182	0.075	56	57.13
8	25	0.455	0.189	35.28	35.99	64	67	0.455	0.189	0	0
25	26	0.364	0.151	56	57.13	67	68	0.910	0.378	0	0
26	27	0.546	0.226	0	0	68	69	1.092	0.453	56	57.13
27	28	0.273	0.113	56	57.13	69	70	0.455	0.189	0	0
28	29	0.546	0.226	0	0	70	71	0.546	0.226	35.28	35.99
29	30	0.546	0.226	35.28	35.99	67	72	0.182	0.075	56	57.13
30	31	0.273	0.113	35.28	35.99	68	73	1.184	0.491	0	0
31	32	0.182	0.075	0	0	73	74	0.273	0.113	56	57.13
32	33	0.182	0.075	14	14.28	73	75	1.002	0.416	35.28	35.99
33	34	0.819	0.340	0	0	70	76	0.546	0.226	56	57.13
34	35	0.637	0.264	0	0	65	77	0.091	0.037	14	14.28
35	36	0.182	0.075	35.28	35.99	10	78	0.637	0.264	56	57.13
26	37	0.364	0.151	56	57.13	67	79	0.546	0.226	35.28	35.99
27	38	1.002	0.416	56	57.13	12	80	0.728	0.302	56	57.13
29	39	0.546	0.226	56	57.13	80	81	0.364	0.151	0	0
32	40	0.455	0.189	35.28	35.99	81	82	0.091	0.037	56	57.13
40	41	1.002	0.416	0	0	81	83	1.092	0.453	35.28	35.99
41	42	0.273	0.113	35.28	35.99	83	84	1.002	0.416	14	14.28
41	43	0.455	0.189	35.28	35.99	13	85	0.819	0.340	35.28	35.99

.

For all numerical simulations, we consider the voltage and power bases for this test system as 11 kV and 1000 kVA, respectively. In reference, [33] reported the amount of active and reactive power losses for this test system: 316.12 kW and 198.60 kvar, respectively. These data are important since they form the losses reference for all the power flow methods studied.

4. Computational Implementation

To solve the power flow problem in the IEEE 34- and IEEE 85-bus systems, all the studied methods as well as the classical Newton–Raphson approach are considered [23]. All the numerical implementations are made in the MATLAB programming environment on a desk computer with INTEL(R) Core(TM) i5-3550, 3.50 GHz, and 8 GB RAM with 64-bit Windows 7 Professional.

To determine the average processing time of each power flow method as well as that of the Newton–Raphson approach, we evaluate each method 100,000 consecutive times. In addition, to determine the total power losses of the network, the following formula is used.

$$\begin{array}{c}\hfill P_{\mathrm{loss}}=\mathrm{real}\left\{{\mathbb{V}}^{\top }{\mathbb{Y}}_{\mathrm{bus}}\mathbb{V}\right\}.\end{array}$$

(28)

4.1. Power Flow Results in the IEEE 34-Bus System

The power flow solution for the IEEE 34-bus system is provided in

Table 4. Power flow performance in the IEEE 34-bus system.

Method	Power Losses (kW)	No. of Iterations	Proc. Time (ms)
Newton–Raphson (NR)		4	2.9929
Successive approximations (SA)		8	0.1893
Matricial backward/forward (MBF)		8	0.6405
Triangular-based (TB)	221.752357	8	0.1662
Product linearization (PL)		4	0.9057
Hyperbolic linearization (HL)		4	0.8568
Diagonal approximation (DA)		8	0.8696

. Numerical results given in this table indicate that (i) all numerical methods lead to the same power losses value for the IEEE 34-bus system with the first six decimals being equal, i.e., $221.752357$ kW; this confirms that all these methods are adequate to solve the power flow problem in radial distribution systems; (ii) the triangular-based power flow approach is the faster method to solve the power flow problem, taking about $0.1662$ ms, which is followed by the successive approximation method with $0.1893$ ms, while the slowest method is the classical Newton–Raphson method that takes almost 3 ms to solve the power flow problem; and (iii) the number of iterations of the derivative-free methods is higher than that of the methods using derivatives in their formulation; this is associated with the convergence rate, as will be discussed in the explanation of Figure 3.

Figure 3 presents the convergence rate of each power flow method studied and its comparison with the classical Newton–Raphson method. It is possible to observe from the figure that the derivative methods have a quadratic convergence, since the logarithm of the error $log\left(ϵ\right)$ decreases with the square of the number of iterations, i.e., $log\left(ϵ\right)=-\alpha t^2$, where $\alpha$ is a positive constant, while all the derivative-free methods converge linearly, i.e., $log\left(ϵ\right)=-\alpha t$, which implies that the number of iterations will be higher when compared with the derivative-based methods.

On the other hand, Figure 4 presents the percentages of the processing time of each studied power flow method when the Newton–Raphson is considered the reference method.

Three main results can be observed in this graphic: (i) even if the matricial backward/forward method and the successive approximation power flow method are mathematically equivalent, the former method takes additional processing time to solve the optimization problem, which is attributed to the need to construct an incidence matrix to generate the ${\mathbb{Y}}_{\mathrm{bus}}$ components; in contrast, the latter method constructs this matrix directly, avoiding additional matricial calculations that help reduce the time required to solve the power flow problem; (ii) the proposed diagonal approximation power flow method, even though similar to the successive approximation approach, takes about 4.6 times more time to solve the power flow problem, attributable to the fact that the diagonal method requires the calculation of the inverse matrix $\mathbf{diag}\left({\mathbb{V}}_d^{t,\star }\right){\mathbb{Y}}_{dd}$ at each iteration, which extends the total processing time required to solve the problem; and (iii) the second slowest method is the product linearization method since it takes about $30.26\%$ of the time required in the Newton–Raphson method, followed by the diagonal approximation approach with $29.06\%$; however, numerical results concerning the IEEE 34-bus system show that all the methods improve the Newton–Raphson time by more than $69\%$, which confirm their effectiveness in solving the studied problem.

4.2. Power Flow Results in the IEEE 85-Bus System

The power flow solution for the IEEE 85-bus system is given in

Table 5. Power flow performance in the IEEE 85-bus system.

Method	Power Losses (kW)	No. of Iterations	Proc. Time (ms)
Newton-Raphson (NR)		5	20.8949
Successive approximations (SA)		11	1.0791
Matricial backward/forward (MBF)		11	2.7924
Triangular-based (TB)	316.117496	11	0.5450
Product linearization (PL)		5	4.6976
Hyperbolic linearization (HL)		4	3.3564
Diagonal approximation (DA)		11	5.3632

, wherein the numerical results suggest that: (i) all methods calculate the total grid power losses with the same level of efficiency, i.e., up to six decimals, with the value being $316.117496$ kW; (ii) the Newton–Raphson method takes more than $20$ms to solve the power flow problem, being the slowest method, while the triangular-based approach requires a total processing time of less than $0.6$ ms, followed only by the successive approximation method with less than $1.1$ms; and (iii) the number of iterations for the Newton–Raphson method and the product linearization method increases to five iterations in comparison with the IEEE 34-bus system, while for the hyperbolic linearization approach, it remains at four iterations, and for the derivative-free methods number of iterations increase from 8 to 11 for this test feeder. These increments are mainly associated with the increment in the number of nodes of the system; however, in the case of the hyperbolic approximation, the rate of convergence is better, which demonstrates that it is possible to reach the desired convergence error with fewer iterations.

On the other hand, the results in Figure 5 confirm that the convergence rate of the derivative-based methods is quadratic, while that of the derivative-free approaches is linear, with the hyperbolic linearization approach showing the most efficient convergence rate with a slope higher than the Newton–Raphson method and the product linearization approach.

Finally, Figure 6 shows the comparison in percentage of the processing times of all the studied methods compared with the Newton–Raphson approach.

These processing times confirm that in the case of radial grids, the best power flow approach is the triangular-based power flow method, followed by the successive approximation method, with less than $2.7$% and $5.2$%, respectively, of the total processing time taken by the Newton–Raphson approach. In addition, all the methods, when compared with the Newton–Raphson approach, are possibly at least $74\%$ faster in the IEEE 85-bus system.

4.3. Additional Comments

Some important facts were generally observed with the power flow solutions for the IEEE 34- and IEEE 85-bus system when all the four derivative-free methods and the two derivative-based approaches are analyzed.

• All derivative-free methods exhibit a linear convergence since their formulation is merely based on the reorganization of the power flow equations in an iterative manner, which implies that information regarding the gradient direction is not included to accelerate their performance regarding the number of iterations, while the derivative-based methods (including the Newton–Raphson approach) have in their formulation variable gains (marices) that help with the reduction of the voltage error after each iteration, which causes these methods to have a quadratic convergence rate.
• The diagonal approximation method, as well as the successive approximation method, is derived from the same power flow formula (see the second row of (4)); therefore, their numerical performance is very similar regarding the number of iterations and the convergence rate. However, the main advantage of the successive approximation method (6) over the diagonal approximation method (27) is the lesser total processing time required to solve the power flow problem since the former method uses a constant matrix that is once time inverted and stored, while the latter method requires inverse calculation at each time, leading to an additional computational effort.
• All the studied methods provide the power flow solution regarding the final power losses estimation and voltage calculations; in fact, any one of them can be selected as a power flow tool for specialized optimization algorithms; however, the importance of knowing the total processing times deals with the selection of the most adequate method for specialized algorithms that evaluate thousands of power flows since, in these recursive optimization algorithms, small differences in the processing times of the power flow can produce important differences in the total execution time of the complete optimization algorithm.

5. Conclusions and Future Works

The classical problem of power flow calculation in AC distribution networks was addressed in this research with six different iterative methods that include four derivative-free and two derivative-based methods; the former methods exhibited linear convergence rate, while the latter methods present a quadratic convergence rate (including the Newton–Raphson approach). The convergence rate produced a direct effect on the total number of iterations to reach the desired tolerance error, with the hyperbolic linearization method with the highest slope that allows finding the solution of the power flow method in four iterations for the IEEE 85-bus system, while the remainder methods take five or more numbers of iterations to solve the power flow problem.

In terms of processing time, the triangular-based approach is the better methodology for radial distribution networks with processing times of $0.1662$ ms and $0.5450$ ms in the IEEE 34- and IEEE 85-bus systems, respectively, which is followed by the successive approximation method with $0.1893$ ms and $1.0791$ ms, respectively. However, all the derivative-methods starting with the present Newton–Raphson, product linearization, and hyperbolic linearization methods take higher processing times since, in each iteration, some matrices are updated to be inverted. The same happened with the diagonal approximation power flow method as was evidenced in both test feeders. Notwithstanding, for the IEEE 34- and IEEE 85-bus system, it was demonstrated that all the numerical methods improve the processing times of the Newton–Raphson approach more than $69\%$ and $74\%$, respectively.

Future works could possibly deal with the following topics: (i) to extend the six studied methods to three-phase electrical networks with $\Delta -$ and $Y-$connected loads considering grid imbalances and also extend these methods for bipolar DC networks; and (ii) to compare the first iteration of these power flow methods (linear power flow solution) with different linear power approaches in the current literature to determine the best power flow approximation in AC distribution grids.