On the Matricial Formulation of Iterative Sweep Power Flow for Radial and Meshed Distribution Networks with Guarantee of Convergence

Abstract

This paper presents a general formulation of the classical iterative-sweep power flow, which is widely known as the backward–forward method. This formulation is performed by a branch-to-node incidence matrix with the main advantage that this approach can be used with radial and meshed configurations. The convergence test is performed using the Banach fixed-point theorem while considering the dominant diagonal structure of the demand-to-demand admittance matrix. A numerical example is presented in tutorial form using the MATLAB interface, which aids beginners in understanding the basic concepts of power-flow programming in distribution system analysis. Two classical test feeders comprising 33 and 69 nodes are used to validate the proposed formulation in comparison with conventional methods such as the Gauss–Seidel and Newton–Raphson power-flow formulations.

1. Introduction

1.1. General Context

Electrical distribution networks represent the most significant portion of the power system and are entrusted with providing electricity to end-users at medium- and low-voltage levels by connecting the transmission/sub-transmission system with consumers [1,2]. Typical voltages in these electrical networks oscillate between 4.16 kV and 33 kV [3], and these networks are typically operated with a radial structure to minimize installation, operation, and maintenance costs, including protective devices coordination simplifications [4,5,6]. To determine the behavior of electrical networks under well-defined voltage conditions, power flow methodologies are used, which facilitate the calculation of the voltage profiles at all the nodes of the network for a particular load condition [7,8]. The main challenge in the power-flow analysis of distribution networks is the constant power loads that produce nonlinear relationships between the voltages and powers [9,10], which makes the use of numerical methods for solving the power-flow problem necessary [11]. A typical tendency in the power-flow analysis of electrical distribution networks is the use of graph-based methods to address the power-flow problem based on the radial structure of the grid [12,13]. However, these methodologies are not useful in the case of weakly or strongly meshed distribution networks [14], which can cause major issues in modern power systems where meshed structures can help with grid performance in terms of realizing lower power losses and improvement voltage profiles [15]. In the case of commercial solutions for a distribution system analysis such as the DigSILENT and ETAP software, the Newton–Raphson (NR) method is the most-used approach [16] as it can be used with radial and meshed structures as well as multiple slack and voltage-controlled nodes [17]. In this research, we tackle the power-flow problem in electrical distribution networks while considering a unique slack node and radial and meshed configurations by reformulating the conventional backward–forward power flow using the branch-to-nodal admittance matrix [18].

1.2. Motivation

The power-flow analysis in electrical distribution networks is one of the most classical and largely studied problems in power-system analysis [14] as the power-flow is a necessary calculation in the determination of grid performance [10]. A power-flow method is a tool that is used to determine, via an iterative procedure, the final values of all the voltages in an electrical network with a certain tolerance acceptance in order to determine the complete operative state [14], i.e., currents through lines, voltage regulation in all the nodes, active and reactive power losses, and voltage stability. To solve the power-flow problem, the most classical methods used are the NR and Gauss–Seidel (GS) approaches, which guarantee convergence based on the Kantarovich and Banach fixed-point theorems [19]. In addition, some graph-based methods such as triangular formulation and quasi-symmetric matrices have also been reported for handling radial grid configurations [12,20]. However, this research is motivated by the fact that the classical backward–forward power-flow method is commonly formulated via sequential steps that require that the grid is ordered in layers [21]. It has a radial structure as, all the currents are calculated in the backward stage, while all the voltages are defined in the forward stage [13]. As this operation can be efficient, these stages exclude meshed configurations, which implies that it is not applicable to weakly or strongly meshed distribution networks. Therefore, herein, a reformulation of the backward–forward power flow in distribution networks is presented based on the branch-to-nodal incidence matrix [13], which permits the handling of radial and meshed distribution networks that include voltage-controlled nodes. Furthermore, the convergence of the proposed method is demonstrated by applying the Banach fixed-point theorem [22].

1.3. Literature Review

Power-flow solutions in the literature present a vast universe of possibilities based on iterative procedures, linearization methods, and convex reformulations. Some of these approaches are presented below.

The most classical power-flow methods in power-systems analysis correspond to the Gauss-Seidel (GS) and NR approaches [14]. The GS approach can be easily implemented using any programming language with its main advantage being that it can be used with complex numbers. However, it exhibits the worse performance in terms of processing time and number of iterations. With the purpose of improving the efficiency of the GS method, an accelerating factor $\alpha$ that reduces the total number of iterations and the required overall processing time is used [17]. In the case of the NR approach, it is a methodology that is more commonly used in commercial software and is widely adopted by utilities [23]. The main advantage of the NR approach is that it converges in a few iterations, has low processing times, and can possibly be used with radial and meshed networks comprising multiple voltage-controlled nodes. The main problem with using the NR approach in distribution networks is the high dispersion of the Jacobian matrix, which may cause convergence problems when the matrix is inverted for updating the voltage magnitude and angle variables [19]. The authors of [24] proposed the use of the Levenberg–Marquardt (LM) algorithm for solving power-flow problems in electrical networks; this approach essentially comprises an alternative manner of presenting the NR formulation with the main advantage that a factor is added to the Jacobian in order to reduce the possible singularities in this matrix.

For distribution networks, some of the most recognized approaches for addressing the power-flow problem are graph-based methodologies known as backward–forward and triangular methods [12,20]. The backward–forward method is typically designed with an algorithmic structure using an iterative sweep based on currents in the backward stage to update the voltages in the forward stage [21]. An important improvement in the backward–forward approach was presented by [13], wherein radial and meshed structures were included in a matricial formulation. The triangular approach works for radial and meshed structures [12]; however, the related convergence analysis is not provided. It should be noted that in these methods, the authors do not present an analysis of the inclusion of voltage-controlled nodes, which we proposed in our matricial reformulation. A new approach based on successive approximations was recently reported in [25]. This approach demonstrates the convergence of the recursive formulation for radial and meshed networks using the Banach fixed-point theorem. However, voltage-controlled sources were not considered in this formulation, and the behavior of the $\theta$-coefficient was not determined to confirm the convergence at well-defined voltage conditions. Two additional important power flow contributions have been reported in specialized literature by [26,27] regarding secondary distribution networks to study the problem of the harmonic penetration caused by higher penetration of photovoltaic sources. The authors in these approaches have proposed a harmonic power flow analysis using affine arithmetic which allows reaching better results regarding harmonic distortion estimations with low-computational effort when compared with classical Monte-Carlo simulations; in addition, these methods were validate in real secondary distribution networks with promissory results for utilities.

Other approaches have also been presented in specialized literature in relation to linear and convex approximations for power-flow solutions in distribution networks. In [11], a convex approximation for a power-flow analysis of radial and meshed distribution networks was presented while considering a semidefinite programming approximation by using a rectangular representation of the power-balance equations. In [8], a second-order cone programming approximation was described while considering branch variables related to currents and powers in all the lines. It is essential to mention that the main disadvantage of these approaches is the low rate of convergence when the number of nodes (n) increases as these formulations create $n^2$ variables related to voltages, which increase the required processing times to reach the optimal solution in the iterative procedure [28]. In the case of linear methodologies, the authors of [7] proposed a Taylor’s series expansion of the hyperbolic relation between the voltage and powers to reach a linear formulation that can be used to solve directly without the use of iterative procedures. In [10], an improvement of this linear approximation was presented using the straight equation to include the minimum expected voltage in the approximation and improve the final result of the power-loss calculation. The authors of [29] described a linear approximation based on a logarithmic transformation of the voltage magnitudes added to the Taylor’s series expansion for applications in power systems. It should be noted that these linear approximations have speed responses; however, the final result of the power-loss estimation varies in the range of 1% and 10%.

1.4. Contributions and Scope

The main contributions of this research can be summarized as follows:

• The reformulation of the classical backward–forward power-flow method for application to distribution networks with the ability to handle radial and meshed configurations by rewriting the branch variables into nodal variables using the branch-to-node incidence matrix.
• The parametric independence of the power-flow formulation, as no assumptions about relations reactance/resistance, are required in the proposed matricial formulation.
• The possibility of guaranteeing convergence under well-defined voltage conditions by applying the Banach fixed-point theorem to the recursive solution, which only requires that the short-circuit current be more significant to the load current in all the nodes to ensure the convergence of the algorithm.
• The presentation of the proposed matricial formulation in an intuitive manner to introduce students of electrical engineering to power-flow analysis by providing the MATLAB/OCTAVE algorithm for solving a small test feeder as a numerical example comprising radial and meshed structures.

It is important to mention that a similar approach was previously proposed in [13] that combined branch-to-nodal matrices with triangular matrices; moreover, its convergences were also proved. However, the authors of [13] did not present the behavior of the $\theta$-coefficient as a function of the load current and the short-circuit current relations per node. In addition, a simulation case is presented herein, where the power system contains a voltage-controlled node (VC node) in a small transmission network, which demonstrates that the proposed matricial approach is extensible to power system applications with multiple meshes and VC nodes.

1.5. Organization of the Document

The remainder of this paper is arranged as follows. Section 2 presents the matricial formulation of the classical backward–forward power flow method using a small distribution test feeder with multiple meshes. In Section 3, the convergence analysis is presented based on the Banach fixed-point theorem, which results in the maximum load current divided by the minimum short-circuit current per node. Section 4 presents a small numerical example comprising a 13-node test feeder with mesh and radial topologies and presents the MATLAB/OCTAVE implementation that can be used as a power-flow analysis tutorial for beginners. In Section 4, the configurations of two classical distribution test feeders comprising 33 and 69 nodes are presented. Section 6 presents the numerical results obtained for these test feeders on comparing the proposed matricial formulation with the classical power- flow methods such as the GS, NR, and LM methods; in addition, a small transmission network comprising four nodes that comprise a meshed structure with a voltage-controlled node is presented, which aids in understanding the extension of the proposed formulation to power-system analyses. Section 7 presents the main conclusions derived from this work as well as its possible future improvements.

2. Matricial Power-Flow Formulation

The classical backward–forward power-flow method was developed using the concept of iterative sweeps, which achieves the following [30]:

• Calculates the total current demands at all the loads with assumed known voltages.
• Determines all the currents that flow in all the branches of the network by applying the first Kirchhoff’s law at each node.
• Calculates the voltage drops by starting from the source and using an ordering stage that defines the layers at which nodes are located.
• Updates the voltage profile in all the demand nodes and repeats all the stages until the convergence tolerance is reached.

This structure of the backward–forward method has two main disadvantages: (i) It can only handle radial grid configurations, and (ii) it requires nodal ordering in layers [20]. However, to address these problems, it is possible to reformulate a power-flow problem using branch-to-nodal incidence matrices [18].

Definition 1

(Branch-to-node incidence matrix $\mathcal{A}$). An electrical network that is represented by a connected graph with b links and n vertices can be represented by a rectangular matrix $\mathcal{A}\in {\mathbb{R}}^{b\times n}$ by assuming arbitrary directions in the current flow through the lines as follows:

• ${\mathcal{A}}_{i,j}=1$ if the line i is connected to the node j and its current leaves from this node.
• ${\mathcal{A}}_{i,j}=-1$ if the line i is connected to the node j and its current arrives at this node.
• ${\mathcal{A}}_{i,j}=0$ if the line i is not connected to the node j.

To exemplify the structure of the branch-to-node incidence matrix $\mathcal{A}$, let us consider the electrical network depicted in Figure 1. This grid comprises six links (branches) and five nodes. The names of the lines are listed in

Table 1. Node connections.

Branch Number i	Send Node j	Receiving Node k
1	1	2
2	1	3
3	2	4
4	3	4
5	3	5
6	4	5

.

On considering the branch information reported in Table 1 and the electrical-grid configuration depicted in Figure 1, the branch-to-node matrix $\mathcal{A}$ can take the following form:

$$\begin{array}{c}\hfill \mathcal{A}=\left[\begin{array}{ccccc}1& -1& 0& 0& 0\\ 1& 0& -1& 0& 0\\ 0& 1& 0& -1& 0\\ 0& 0& 1& -1& 0\\ 0& 0& 1& 0& -1\\ 0& 0& 0& 1& -1\end{array}\right]\end{array}$$

(1)

It should be noted that this matrix can be split into two sub-matrices that are related to the slack source and the rest of the nodes as follows (This is made possible by assuming the slack node is located at node 1.):

$$\begin{array}{c}\hfill \mathcal{A}=\left[{\mathcal{A}}_s{\mathcal{A}}_d\right],\end{array}$$

(2)

where ${\mathcal{A}}_s$ is the first column of $\mathcal{A}$ (slack node), and ${\mathcal{A}}_d$ corresponds to the rest of the columns in $\mathcal{A}$, i.e., the demand nodes.

Now, let us define the voltage drop at line i as ${\mathbf{E}}_i$ and the voltages at its terminals as ${\mathbf{V}}_j$. Then, from Figure 1, we know that

$$\begin{array}{c}\hfill {\mathbf{E}}_1={\mathbf{V}}_1-{\mathbf{V}}_2,\\ \hfill {\mathbf{E}}_2={\mathbf{V}}_1-{\mathbf{V}}_3,\\ \hfill {\mathbf{E}}_3={\mathbf{V}}_2-{\mathbf{V}}_3,\\ \hfill {\mathbf{E}}_4={\mathbf{V}}_2-{\mathbf{V}}_4,\\ \hfill {\mathbf{E}}_5={\mathbf{V}}_3-{\mathbf{V}}_4,\\ \hfill {\mathbf{E}}_6={\mathbf{V}}_4-{\mathbf{V}}_5,\end{array}$$

which can be rewritten using (1) and (2) as presented below:

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

(3)

where ${\mathbf{V}}_s={\mathbf{V}}_1$, $\mathbf{E}$ is a vector that comprises all the voltage drops, and ${\mathbf{V}}_d$ is a vector that comprises all the voltages in the demand nodes.

To relate the currents $\mathbf{J}$ in all the branches with the voltage drops, it is possible to apply Ohm’s law as presented below:

$$\begin{array}{c}\hfill \mathbf{E}=\mathbb{Z}\mathbf{J},\end{array}$$

(4)

where $\mathbb{Z}=\mathbf{diag}(Z_1,Z_2,\dots ,Z_b)$.

To determine the relation between the currents injected at the nodes (i.e., I) and all the branch currents, let us apply Kirchhoff’s second law at each node while considering the current directions defined in Table 1. We thus obtain

$$\begin{array}{cc}\hfill {\mathbf{I}}_1& ={\mathbf{J}}_1+{\mathbf{J}}_2,\hfill \\ \hfill {\mathbf{I}}_2& =-{\mathbf{J}}_1+{\mathbf{J}}_3,\hfill \\ \hfill {\mathbf{I}}_3& =-{\mathbf{J}}_2+{\mathbf{J}}_4+{\mathbf{J}}_5,\hfill \\ \hfill {\mathbf{I}}_4& =-{\mathbf{J}}_3-{\mathbf{J}}_4+{\mathbf{J}}_6,\hfill \\ \hfill {\mathbf{I}}_5& =-{\mathbf{J}}_5-{\mathbf{J}}_6.\hfill \end{array}$$

It should be noted that $\mathbf{I}$ can also be split as

$$\begin{array}{c}\hfill \mathbf{I}={\left[{\mathbf{I}}_s{\mathbf{I}}_d\right]}^T,\end{array}$$

(5)

where ${\mathbf{I}}_s$ and ${\mathbf{I}}_d$ are the current of generators and demands, respectively.

Using the relation between the branch and nodal currents presented above, it is possible to find that

$$\begin{array}{c}\hfill {\mathbf{I}}_s={\mathcal{A}}_s^T\mathbf{J},\end{array}$$

(6)

$$\begin{array}{c}\hfill {\mathbf{I}}_d={\mathcal{A}}_d^T\mathbf{J}.\end{array}$$

(7)

Now if expressions (3), (4), and (7) are combined based on $\mathbb{Y}={\mathbb{Z}}^{-1}$, then the following result is obtained:

$$\begin{array}{c}\hfill {\mathbf{I}}_d={\mathcal{A}}_d^T\mathbb{Y}{\mathcal{A}}_s{\mathbf{V}}_s+{\mathcal{A}}_d^T\mathbb{Y}{\mathcal{A}}_d{\mathbf{V}}_d,\end{array}$$

(8)

To find the relation between the injected currents and powers, it is possible to apply the Tellegen’s theorem, which defines apparent power as follows:

$$\begin{array}{c}\hfill {\mathbf{I}}_s={\mathbf{diag}}^{-1}\left({\mathbf{V}}_s^{\star }\right){\mathbf{S}}_s^{\star },\end{array}$$

(9)

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

(10)

where the superscript $\star$ represents the conjugate operator, ${\mathbf{S}}_s^{\star }$ is the apparent power generation in the slack node, and ${\mathbf{S}}_d^{\star }$ are the apparent power consumption in the demand nodes.

Now, if expressions (8) and (10) are combined to obtain ${\mathbf{V}}_d$, the following result is obtained:

$$\begin{array}{c}\hfill {\mathbf{V}}_d=-{\mathbf{Z}}_{dd}{\mathbf{diag}}^{-1}\left({\mathbf{V}}_d^{\star }\right){\mathbf{S}}_d^{\star }-\mathbf{b},\end{array}$$

(11)

where ${\mathbf{Z}}_{dd}$ is defined as ${\left[{\mathcal{A}}_d^T\mathbb{Y}{\mathcal{A}}_d\right]}^{-1}$, and $\mathbf{b}=+{\mathbf{Z}}_{dd}{\mathcal{A}}_d^T\mathbb{Y}{\mathcal{A}}_s{\mathbf{V}}_s$.

It should be noted that to solve expression (11), an iterative procedure is required, as this expression has a recursive structure. For this purpose, let us to add an iterative counter t that starts from the initial point $t=0$ (i.e., ${\mathbf{V}}_d^0$ is perfectly known), which allows the rewriting of (11) as presented below:

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

(12)

Remark 1.

The iterative procedure in (12) is performed until the convergence error ϵ is reached, i.e., $\max \left\{{\mathbf{V}}_d^{t+1}-{\mathbf{V}}_d^t\right\}\le ϵ$ or the total number of iterations $t_{\max }$ is completed.

3. Convergence Analysis

To prove that the recursive formula (12) guarantees the convergence of the power-flow problem in the case of distribution networks comprising radial and mesh structures formulated using branch-to-nodal incidence matrices, let us consider the following.

Assumption 1.

The amount of power consumed by the loads does produce voltage instabilities in the distribution network, i.e., the system operates under well-defined voltage conditions.

Assumption 2.

There exists a positive lower bound for all the voltage profiles, i.e., $V^{\min }>0$, which is typically defined by the utility and regulatory policies.

Assumption 3.

The demand-to-demand impedance matrix, i.e., ${\mathbb{Z}}_{dd}$, is diagonal dominant, which implies that, always, $\left|{\mathbb{Z}}_{d{d}_{jj}}\right|\ge \left|{\mathbb{Z}}_{d{d}_{jk}}\right|,\forall j\ne k$.

To demonstrate the convergence properties of the recursive power flow approach defined by (12), let us present the general definition of the Banach fixed-point theorem as follows.

Theorem 1

(Banach fixed-point theorem). The iterative formula presented in (12) is stable and corresponds to a contraction map with the form

$$\begin{array}{c}\hfill {\mathbf{V}}_d^{t+1}=g\left({\mathbf{V}}_d^t\right),\end{array}$$

(13)

for some $\mathbf{V}$ that fulfills Assumption 1 irrespective of the starting point, i.e., ${\mathbf{V}}^0$, such that

$$\begin{array}{c}\hfill \left|\left|g\left({\mathbf{V}}_d^0\right)-g\left(\mathbf{W}\right)\right|\right|\le \theta \left|\left|{\mathbf{V}}^0-\mathbf{W}\right|\right|\end{array}$$

(14)

where $\mathbf{W}$ is the solution of the power-flow problem, and θ is a real number between 0 and 1.

Proof. 

The recursive formula for the iterative backward–forward power-flow method presented in (12) can be rewritten as follows:

$$\begin{array}{c}\hfill {\mathbf{V}}_d^{t+1}=g\left({\mathbf{V}}_d^t\right)=-\mathbf{b}-{\mathbf{Z}}_{dd}{\left[\frac{{\mathbb{S}}_{d_i}^{\star }}{{\mathbb{V}}_{d_i}^{t,\star }}\right]}_{i\in {\mathsf{\Omega }}_S}^T\end{array}$$

(15)

where ${\mathsf{\Omega }}_S$ is the set that contains all the demand nodes.

In addition, based on the structure of the Banach fixed-point theorem, it is possible to affirm that the solution of the power-flow problem, i.e., $\mathbb{W}$, satisfies $\mathbb{W}=g\left(\mathbb{W}\right)$, and it is unique if and only if $g\left(\mathbb{W}\right)$ is a contraction mapping on ${\mathbb{V}}_d$.

$$\begin{array}{cc}\hfill \left|\left|{\mathbf{V}}_d^{t+1}-\mathbf{W}\right|\right|=& \left|\left|g\left({\mathbf{V}}_d^{t+1}\right)-g\left(\mathbf{W}\right)\right|\right|\hfill \\ \hfill & \left|\left|{\mathbb{Z}}_{dd}{\mathbf{S}}_d^{\star }{\left[\frac{1}{{\mathbb{W}}_i^{\star }}-\frac{1}{{\mathbb{V}}_{d_i}^{t,\star }}\right]}_{i\in {\mathsf{\Omega }}_S}^T\right|\right|\hfill \\ \hfill & \left|\left|{\mathbb{Z}}_{dd}{\mathbf{S}}_d^{\star }\right|\right|\left|\left|{\left[\frac{{\mathbb{V}}_{d_i}^{t,\star }-{\mathbb{W}}_i^{\star }}{{\mathbb{W}}_i^{\star }{\mathbb{V}}_{d_i}^{t,\star }}-\right]}_{i\in {\mathsf{\Omega }}_S}^T\right|\right|\hfill \\ \hfill & \le \theta \left|\left|{\mathbf{V}}_d^t-\mathbf{W}\right|\right|,\hfill \end{array}$$

(16)

where

$$\begin{array}{c}\hfill \theta =\frac{\left|\left|{\mathbb{Z}}_{dd}{\mathbf{S}}_d^{\star }\right|\right|}{{\left({\mathbf{V}}^{\min }\right)}^2}.\end{array}$$

(17)

Now, on considering Assumption 2, which is regarding the positive definiteness of the impedance matrix, the expression (17) can be reduced to

$$\begin{array}{c}\hfill \theta =\underset{i\in {\mathsf{\Omega }}_S}{\max }\left\{\frac{\left|{\mathbb{Z}}_{d{d}_{ii}}\right|\left|{\mathbb{S}}_{d_i}^{\star }\right|}{{\left({\mathbf{V}}^{\min }\right)}^2}\right\}.\end{array}$$

(18)

It should be noted that, considering the mathematical structure presented in (18) and that ${\mathbb{Z}}_{d{d}_{ii}}$ is the equivalent impedance at node i, the following relation can be obtained.

$$\begin{array}{c}\hfill \theta =\underset{i\in {\mathsf{\Omega }}_S}{\max }\left\{\frac{\frac{\left|{\mathbb{S}}_{d_i}^{\star }\right|}{{\mathbf{V}}^{\min }}}{\frac{{\mathbf{V}}^{\min }}{\left|{\mathbb{Z}}_{d{d}_{ii}}\right|}}\right\},\end{array}$$

(19)

where we can ensure that $0\le \theta \le 1$, as the denominator in (19) corresponds to the minimum short-circuit current, and the numerator represents the maximum load current, which is always lower than any short-circuit current during normal operative conditions. This implies that the recursive matricial sweep power flow (12) converges to the power-flow solution, which completes the proof.  □

4. Numerical Example

This section presents a small test feeder comprising 13 nodes that can handle radial or mesh configurations. This facilitates the demonstration of the convergence of the proposed matricial reformulation of the classical backward—forward power flow approach. Figure 2 presents the electrical configuration of the 13-node test feeder.

Table 2. Branch and load parameters of the 13-node test feeder.

Node i	Node j	${\mathit{R}}_{\mathbf{ij}}$ [$\mathbf{\Omega }$]	${\mathit{X}}_{\mathbf{ij}}$ [$\mathbf{\Omega }$]	${\mathit{P}}_{\mathit{j}}$ [kW]	${\mathit{Q}}_{\mathit{j}}$ [kvar]
1	2	0.3968	0.5290	2000	1600
2	3	0.4232	0.5819	3000	400
2	4	0.4761	0.9522	2000	−400
4	5	0.2116	0.2116	1500	1200
1	6	0.5819	0.5819	4000	2700
6	7	0.4232	0.5819	5000	1800
6	8	0.5819	0.5819	1000	900
7	9	0.5819	0.5819	600	−400
1	10	0.5819	0.5819	1000	900
10	11	0.4761	0.6348	1000	−1100
10	12	0.4232	0.5819	1000	900
12	13	0.2116	0.2116	2100	−800
Tie-lines
3	9	0.2116	0.2116	—	—
8	11	0.2116	0.2116	—	—
5	13	0.4761	0.6348	—	—

lists the branch and load information for the 13-node test feeder that has a total demand of active and reactive power of 24,200 kW and 7700 kVAr, respectively, and it operates on a line-to-neutral voltage of 23 kV.

To solve the power-flow problem in the 13-node test feeder, we consider the following: (i) a radial configuration is reached when all the dashed lines in Figure 2 are open, and a meshed configuration is obtained when one or more dashed lines are connected; (ii) the recursive formula (12) was implemented in MATLAB using scripts in a tutorial style.

Figure 3 reports the MATLAB implementation of the proposed matricial backward–forward (MBF) power flow while considering the meshed configuration for the 13-node test feeder depicted in Figure 2.

In the MATLAB implementation presented in Figure 3, we can observe the following:

√

From lines 1 to 36, all the numerical information in the test system is presented and transformed into a per-unit representation.

√

Lines 37 to 48 comprise the branch-to-node incidence matrix, impedance matrix, and all the constant components for evaluating the recursive power-flow formula.

√

From lines 49 to 59, the iterative procedure for solving the power-flow problem while considering a matricial formulation is implemented while considering a convergence tolerance of approximately $ϵ=1\times {10}^{-10}$. It should be noted that lines 53 and 54 present the required calculations for determining the total grid power losses using branch variables, i.e., $\mathbb{J}$ and $\mathbb{E}$, and the line impedance matrix $\mathbb{Z}$.

√

Lines 60 to 65 comprise the $\theta$ coefficients for all the loads for proving that each of them fulfills the Banach fixed-point theorem condition.

It is important to mention that if the last three lines related to the branches are removed in the matrix (see Figure 3), then the radial test system can be evaluated. However, Figure 4 presents the $\theta$-coefficient for the radial and mesh 13-node configurations.

From Figure 4, it can be concluded that expression (19) is fulfilled for all the demand nodes as the $\theta$ values are between 0 and 1. For both the simulation cases, the maximum value of $\theta$ is reached at node 7 and is $0.016$ and $0.009$ in the case of the radial and meshed configurations, respectively. It should be noted that, in the case of the meshed configuration, the minimum voltage is $0.983$ p.u at node 7, and in the radial configuration, this value is $0.976$ p.u at node 9.

5. Test Systems

This section presents two classical distribution test feeders with radial configurations that are typically used in the optimal location of distributed generators [5], the optimal location of capacitor banks [32,33,34], and optimal battery scheduling approaches [35]. The first system corresponds to the 33-node test feeder, and the second system comprises the 69-node test feeder, which is widely known as the Baran & Wu distribution network.

5.1. 33-Three-Node Test Feeder

This test system comprises 33 nodes and 32 branches with an operating voltage of $12.66$ kV. The slack node is located at node 1, and its configuration is illustrated in Figure 5. This feeder has 3715 kW and 2300 kVAr of total active and reactive power demand, respectively. The initial active power losses of this system are equal to $210.988$ kW. Here we considered the voltage and power base values of $12.66$ kV and 1000 kW, respectively.

The information of all the branches as well as the load consumption of the 33-node test feeder are listed in

Table 3. Electrical parameters of the 33-node test feeder.

Node i	Node j	${\mathit{R}}_{\mathbf{ij}}$ [$\mathbf{\Omega }$]	${\mathit{X}}_{\mathbf{ij}}$ [$\mathbf{\Omega }$]	${\mathit{P}}_{\mathit{j}}$ [kW]	${\mathit{Q}}_{\mathit{j}}$ [kW]	Node i	Node j	${\mathit{R}}_{\mathbf{ij}}$ [$\mathbf{\Omega }$]	${\mathit{X}}_{\mathbf{ij}}$ [$\mathbf{\Omega }$]	${\mathit{P}}_{\mathit{j}}$ [kW]	${\mathit{Q}}_{\mathit{j}}$ [kW]
1	2	0.0922	0.0477	100	60	17	18	0.7320	0.5740	90	40
2	3	0.4930	0.2511	90	40	2	19	0.1640	0.1565	90	40
3	4	0.3660	0.1864	120	80	19	20	1.5042	1.3554	90	40
4	5	0.3811	0.1941	60	30	20	21	0.4095	0.4784	90	40
5	6	0.8190	0.7070	60	20	21	22	0.7089	0.9373	90	40
6	7	0.1872	0.6188	200	100	3	23	0.4512	0.3083	90	50
7	8	1.7114	1.2351	200	100	23	24	0.8980	0.7091	420	200
8	9	1.0300	0.7400	60	20	24	25	0.8960	0.7011	420	200
9	10	1.0400	0.7400	60	20	6	26	0.2030	0.1034	60	25
10	11	0.1966	0.0650	45	30	26	27	0.2842	0.1447	60	25
11	12	0.3744	0.1238	60	35	27	28	1.0590	0.9337	60	20
12	13	1.4680	1.1550	60	35	28	29	0.8042	0.7006	120	70
13	14	0.5416	0.7129	120	80	29	30	0.5075	0.2585	200	600
14	15	0.5910	0.5260	60	10	30	31	0.9744	0.9630	150	70
15	16	0.7463	0.5450	60	20	31	32	0.3105	0.3619	210	100
16	17	1.2890	1.7210	60	20	32	33	0.3410	0.5302	60	40

.

5.2. 69-Node Test Feeder

This test system comprises 69 nodes and 68 branches with an operating voltage of $12.66$ kV. The slack node is located at node 1, and its configuration is depicted in Figure 6. This feeder has a total active and reactive power demand of $3890.7$ kW and $2693.6$ kVAr, respectively. The initial active power losses of this system are equal to $225.072$ kW. Here, we assumed that the voltage and power base values were $12.66$ kV and 1000 kW, respectively.

The information of all the branches as well as the load consumption of the 69-node test feeder are listed in

Table 4. Electrical parameters of the 69-node test feeder.

Node i	Node j	${\mathit{R}}_{\mathbf{ij}}$ [$\mathbf{\Omega }$]	${\mathit{X}}_{\mathbf{ij}}$ [$\mathbf{\Omega }$]	${\mathit{P}}_{\mathit{j}}$ [kW]	${\mathit{Q}}_{\mathit{j}}$ [kW]	Node i	Node j	${\mathit{R}}_{\mathbf{ij}}$ [$\mathbf{\Omega }$]	${\mathit{X}}_{\mathbf{ij}}$ [$\mathbf{\Omega }$]	${\mathit{P}}_{\mathit{j}}$ [kW]	${\mathit{Q}}_{\mathit{j}}$ [kW]
1	2	0.0005	0.0012	0	0	3	36	0.0044	0.0108	26	18.55
2	3	0.0005	0.0012	0	0	36	37	0.0640	0.1565	26	18.55
3	4	0.0015	0.0036	0	0	37	38	0.1053	0.1230	0	0
4	5	0.0251	0.0294	0	0	38	39	0.0304	0.0355	24	17
5	6	0.3660	0.1864	2.6	2.2	39	40	0.0018	0.0021	24	17
6	7	0.3810	0.1941	40.4	30	40	41	0.7283	0.8509	1.2	1
7	8	0.0922	0.0470	75	54	41	42	0.3100	0.3623	0	0
8	9	0.0493	0.0251	30	22	42	43	0.0410	0.0475	6	4.3
9	10	0.8190	0.2707	28	19	43	44	0.0092	0.0116	0	0
10	11	0.1872	0.0619	145	104	44	45	0.1089	0.1373	39.22	26.3
11	12	0.7114	0.2351	145	104	45	46	0.0009	0.0012	39.22	26.3
12	13	1.0300	0.3400	8	5	4	47	0.0034	0.0084	0	0
13	14	1.0440	0.3450	8	5.5	47	48	0.0851	0.2083	79	56.4
14	15	1.0580	0.3496	0	0	48	49	0.2898	0.7091	384.7	274.5
15	16	0.1966	0.0650	45.5	30	49	50	0.0822	0.2011	384.7	274.5
16	17	0.3744	0.1238	60	35	8	51	0.0928	0.0473	40.5	28.3
17	18	0.0047	0.0016	60	35	51	52	0.3319	0.1114	3.6	2.7
18	19	0.3276	0.1083	0	0	9	53	0.1740	0.0886	4.35	3.5
19	20	0.2106	0.0690	1	0.6	53	54	0.2030	0.1034	26.4	19
20	21	0.3416	0.1129	114	81	54	55	0.2842	0.1447	24	17.2
21	22	0.0140	0.0046	5	3.5	55	56	0.2813	0.1433	0	0
22	23	0.1591	0.0526	0	0	56	57	1.5900	0.5337	0	0
23	24	0.3460	0.1145	28	20	57	58	0.7837	0.2630	0	0
24	25	0.7488	0.2475	0	0	58	59	0.3042	0.1006	100	72
25	26	0.3089	0.1021	14	10	59	60	0.3861	0.1172	0	0
26	27	0.1732	0.0572	14	10	60	61	0.5075	0.2585	1244	888
23	28	0.0044	0.0108	26	18.6	61	62	0.0974	0.0496	32	23
28	29	0.0640	0.1565	26	18.6	62	63	0.1450	0.0738	0	0
29	30	0.3978	0.1315	0	0	63	64	0.7105	0.3619	227	162
30	31	0.0702	0.0232	0	0	64	65	1.0410	0.5302	59	42
31	32	0.3510	0.1160	0	0	11	66	0.2012	0.0611	18	13
32	33	0.8390	0.2816	14	10	66	67	0.0047	0.0014	18	13
33	34	1.7080	0.5646	19.5	14	12	68	0.7394	0.2444	28	20
34	35	1.4740	0.4873	6	4	68	69	0.0047	0.0016	28	20

.

6. Computational Validation

To validate the proposed MBF power-flow reformulation, we consider four classical iterative power-flow approaches reported in the specialized literature as comparative methodologies: (i) classical GS [14], (ii) accelerated version of GS (AGS) [14], (iii) NR [17], and iv) LM [24]. To perform a fair comparison, all these methods are evaluated 1000 consecutive times to determine their computational costs (processing times) by considering a maximum of 100 iterations and a tolerance error of $1\times {10}^{-10}$.

6.1. Results for the 33- and 69-Node Test Feeders

Table 5. Numerical results for the 33- and 69-node test systems.

33-node Test Feeder
Method	Proc. time [ms]	Iterations	Losses [p.u]
GS	441.974	2313	2.110
AG$(\alpha =1.82)$	38.555	227	2.110
NR	10.751	5	2.110
LM	10.882	5	2.110
MBF	1.323	10	2.110
69-node Test Feeder
GS	31107.756	49031	2.422
AG$(\alpha =1.92)$	1662.691	2455	2.422
NR	38.303	5	2.422
LM	42.719	5	2.422
MBF	5.369	10	2.422

lists the computational efficiency of the proposed MBF method as compared to the classical approaches reported in the literature. These results were reached after performing 1000 consecutive evaluations, which confirms the following:

√

The proposed approach is the most efficient in terms of processing times as compared with the classical approaches as it only takes 1.32 ms for completion in the case of the 33-node test feeder and 5.37 ms for completion in the case of the 69-node test feeder, which implies that it is at least seven times faster than the NR method, which is the most used approach in research and the industry.

√

The classical GS method presents the worse performance in terms of processing times and number of iterations. However, its accelerated version with the $\alpha$-coefficient presents a significant improvement in its performance. In the case of the 33-node test feeder, this method is approximately ten times faster, and in the case of the 69-node test feeder, it is at least 18 times faster, in terms of the total processing time required to solve the power-flow problem.

√

The NR and LM methods belong to the same class as both comprise the use of Jacobian matrices. Therefore, their performance with respect to the number of iterations and the processing times is very similar. Even if the proposed approach requires double the number of iterations, in terms of processing times, we can confirm that the MBF is the best numerical method for distribution power-flow analysis as compared with the NR and LM approaches.

√

In terms of the power losses calculation, it is important to mention that all the numerical methods listed in Table 1 arrive at the same solution with negligible errors of estimation (lower than $1\times {10}^{-10}$). This implies that each of them is suitable for implementation; nevertheless, our proposed method, i.e., the MBF reformulation, is the most desirable approach due to its speediest performance.

Figure 7 presents the voltage profiles for the 33- and 69-node test feeders for the proposed MBF and NR power-flow methods. We can observe that these plots confirm that, numerically speaking, both these methods demonstrate the same performance in terms of voltage determination. However, the main advantage of the MBF when compared to the NR approach (see Table 1) is its faster calculation times as it guarantees convergence on well-defined voltage conditions.

To demonstrate that the proposed MBF approach guarantees convergence in the 33- and 69-node test feeders, Figure 8 presents the behavior of the $\theta$-coefficient. From the results in Figure 8, it is possible to confirm that hypothesis (19) is fulfilled, and the proposed BFM guarantees the convergence of the power-flow problem. In the case of the 33-node test feeder, the maximum coefficient occurs at node 30 with a value of $0.027$, while in the case of the 69-node test feeder, this occurs at node 61 with a value of $0.062$. These values facilitate the verification that $\theta$ is comprehended between 0 and 1 in the case of both the test feeders.

6.2. Inclusion of Voltage-Controlled Nodes

One of the main advantages of the proposed MBF power-flow reformulation is the possibility of working with voltage-controlled nodes, i.e., nodes that control active power injection and the magnitude of the voltage profile. These nodes are common in power-system analyses [14], wherein multiple synchronous machines are interconnected via the transmission network. To realize the VC inclusion in the proposed model, we can apply the classical GS method, where the reactive power generation at the VC source nodes is calculated as follows:

$${\mathbb{S}}_{g_k}^{\mathrm{VC},t}=P_{g_k}^{\mathrm{VC}}+j{Q}_{g_k}^{\mathrm{VC},t},$$

(20a)

$$Q_{g_k}^{\mathrm{VC},t}=-\mathrm{Imag}\left\{{\mathbf{S}}_{d_k}^{\star ,t}-{\mathbf{V}}_{d,k}^{\star ,t}\sum _{j=1}^n{\mathbb{Y}}_{kj}{\mathbf{V}}_{d,j}^{\star ,t}\right\}$$

(20b)

where ${\mathbb{S}}_{g_k}^{\mathrm{VC},t}$ is the apparent power injected at node k by the VC source, which can be split into its real and imaginary powers, i.e., $P_{g_k}^{\mathrm{VC}}$ and $Q_{g_k}^{\mathrm{VC},t}$, which are the active power constant and reactive power variable for maintaining a constant magnitude of the voltage output; ${\mathbb{Y}}_{kj}$ is the admittance that relates nodes k and j, respectively.

The voltage-controlled nodes can be easily included in the proposed MBF method by modifying (12) as follows.

$$\begin{array}{c}\hfill {\mathbf{V}}_d^{t+1}={\mathbf{Z}}_{dd}{\mathbf{diag}}^{-1}\left({\mathbf{V}}_d^{t,\star }\right)\left({\mathbb{S}}_g^{\mathrm{VC},t}-{\mathbf{S}}_d^{\star }\right)-\mathbf{b}.\end{array}$$

(21)

Remark 2.

Once all the voltages are calculated, in the case of the VC source, the magnitude of the voltage is fixed as a constant, and from (21), only the angle at each iteration is updated, i.e.,

$$\begin{array}{c}\hfill {\mathit{V}}_{g_k}^{t+1,VC}=\left|{\mathbf{V}}_{g_k}^{VC}\right|\angle {\mathit{V}}_{g_k}^{t+1,VC}.\end{array}$$

being ${\mathbf{V}}_{g_k}^{VC}$ the specified voltage in the VC source.

To demonstrate the applicability of the proposed MBF power-flow method, we implement a small four-node test feeder with a slack node at bus 1 and a VC source at node 4. This system was initially presented in [14] to introduce the GS and NR methods. The comparison in this section is performed using the widely used power-system software called DigSILENT. Figure 9 presents the configuration of this test feeder and the DigSILENT results obtained after evaluating the AC power flow via the NR method.

The information of the branches and loads for this four-node test feeder are obtained from [14] and are listed in

Table 6. Line parameters of the four-node test feeder.

Line	R [$\mathbf{\Omega }$]	X [$\mathbf{\Omega }$]	${\mathit{Y}}_{\mathbf{shunt}}$ [$\mathsf{\mu }$S]	Current [A]
1–2	5.3323	26.6616	193.7618	200
1–3	3.9358	19.6788	146.5028	400
2–4	3.9358	19.6788	146.5028	500
3–4	6.7289	33.6444	241.0208	400

and

Table 7. Nodal parameters of the four-node test feeder.

Node	${\mathit{P}}_{\mathit{g}}$ [MW]	${\mathit{Q}}_{\mathit{g}}$ [MVAr]	${\mathit{P}}_{\mathit{d}}$ [MW]	${\mathit{Q}}_{\mathit{d}}$ [MVAr]	V [pu]
1 (slack)	—	—	50	30.99	1.00$\angle 0$
2 (load)	0	0	170	105.35	1.00$\angle 0$
3 (load)	0	0	200	123.94	1.00$\angle 0$
4 (VC)	318	—	80	49.58	1.02$\angle 0$

.

Once the power-flow problem is solved using the proposed MBF method, the results reported in

Table 8. Voltages in the four-node test feeder.

Node	DigSILENT	MBF	Node	DigSILENT	MBF
1	$1.000\angle {0.00}^{\circ }$	$1.000\angle {0.00}^{\circ }$	3	$0.969\angle {-1.87}^{\circ }$	$0.969\angle {-1.87}^{\circ }$
2	$0.982\angle {-0.98}^{\circ }$	$0.982\angle {-0.98}^{\circ }$	4	$1.020\angle {-1.52}^{\circ }$	$1.020\angle {-1.52}^{\circ }$

are obtained. It should be noted that these are compared to the DigSILENT solution presented in Figure 9.

The results presented in Table 8 confirm the possibility of using our proposal approach to solve power-flow problems in power systems by using the MBF reformulation with the main advantage that it is comparable with sophisticate software such as DigSILENT. In Appendix A has been illustrated the codification of the proposed MBF method for power systems that comprise voltage-controlled nodes.

Remark 3.

It should be noted that, if the voltage profiles are the same while comparing the NR and MBF methods, as presented in Table 8, then all the active and reactive power flows reported in Figure 9 will demonstrate exactly the same performance, which confirms the possibility of using the proposed approach in radial and meshed networks comprising voltage controlled nodes.

7. Conclusions and Future Works

A matricial version of the classical backward–forward power-flow method has been presented in this paper with the main advantage that it guarantees convergence under well-defined voltage conditions. In addition, numerical validations were performed using 33- and 69-node test feeders to demonstrate that the MBF proposal requires the minimum processing times as compared with those of classical methods such as GS, NR, and LM.

The proof of convergence was performed by applying the Banach fixed-point theorem. This has permitted the conclusion that the proposed MBF method always converges under well-defined load conditions as the $\theta$-coefficient was defined as the relation between the maximum load current and minimum short-circuit current. This implies that $\theta$ is always be comprehended between 0 and 1, thus guaranteeing the numerical convergence of the power-flow problem.

The inclusion of voltage-controlled nodes was also considered by implementing a small four-node test feeder and comparing the numerical results reported using well-known software used for solving power-flow problems called DigSILENT and those of our proposed approach. The voltage profiles obtained using both these methods allowed use to conclude that the proposed MBF method can handle radial and meshed networks comprising VC nodes and will be suitable for multiple power-system applications shortly.

Some additional works that can be derived from this work are listed below: (i) the application of this power-flow method to electrical networks comprising renewable energy sources and battery energy-storage systems in an embedded strategy that facilitates the optimal dispatch of these devices with metaheuristics; (ii) the application of our proposed approach for conducting optimal power-flow analyses in distribution networks via sequential programming approaches based on combinatorial optimization; (iii) the extension of the MBF formulation to electrical systems comprising hybrid structures that contain AC and DC feeders interfaced with power electronic converters.