Two Algorithms for Load Flow Analysis in Balanced Radial Distribution Systems

This paper was developed with the purpose of knowing the behavior of radial distribution systems that operate in balanced conditions. Two algorithms were used to calculate the voltage at the network nodes and with them the load flow in each line . With the application of Kirchhoff's laws of electrical circuits, a one-dimensional arrangement of order N and a two-dimensional arrangement of order Nx5 are generated, both depend on the impedances and the shunt admittances. With these two arrangements, a non-linear function f (V1, V2 ,,, VN,, Sc) was generated, which depends on the voltages Vi and the complex power Sc demanded at each node of the distribution system. These equations were solved iteratively, two methods were applied; In the first method, the Halley formula with an acceleration factor α was used to solve the non-linear function f (V1, V2 ,,, VN,, Sc), in the second method, the Gaussian method was applied in the same way that it is used for load flows in interconnected systems. The two algorithms were applied to test systems of 40, 70 and 80 nodes, obtaining results similar to those reported in the specialized literature where they make use of other solution methods. In the author’s opinion, the proposed algorithms represent alternatives for the study of energy distribution systems. An advantage of the methods is that they are simple code and there are differences between them, the second method has a faster convergence.


Introduction
Distribution systems hold a very significant position in the power system since it is the main point of link between bulk power and consumers. Effective planning of radial distribution network is required to meet the present growing domestic, industrial and commercial load day by day. Distribution networks have gained an overwhelming research interest in the academics as well as in the industries community nearly from last three decades. In the literature, there are a number of efficient and reliable load flow solution techniques, such as; Gauss-Seidel, Newton-Raphson and Fast Decoupled Load Flow [1][2][3]. Hitherto they are successfully and widely used for power system operation, control and planning. However, it has repeatedly been shown that these methods may become inefficient in the analysis of distribution systems with high R/X ratios or special network structures [4][5][6]. Accordingly, a number of methods proposed in the literature [7][8][9] specially designed for the solution of power flow problem in radial distribution networks. The methods developed for the solution of ill-conditioned radial distribution systems may be divided into two categories. The first type of methods is based on forward-backward sweep processes using Kirchooff1s Laws or making use of the well-known bi-quadratic equation which, for every branch, relates the voltage magnitude at the receiving end to the voltage at the sending end and the branch power flow for solution of ladder networks [7][8][9]. On the other hand, the second group of methods is utilized by proper modification of the existing methods such as, Newton-Raphson [2][3][4]. The radial part is solved by a straightforward two step procedure in which the branch currents are first computed (backward sweep) and then the bus voltages are updated (forward sweep). In [11,12] use a method to analyse the distribution system by modelling the load in exponential and polynomial form. In [13], they do a review of the forward and backward methods to do the study of power flows in distribution systems. Possibly, the study of a distribution system is more complicated than the studies that are carried out in interconnected power systems. Harmonics, low industrial power factors, unbalanced operation that always occurs in residential feeders, have a direct impact on the distribution system, making your analysis very complicated. Due to this, many researchers and companies dedicate their effort to carry out many studies and to know the behaviour of the distribution system under different scenarios [14,15,9,8,16,17]. Figure 1 shows a section of a distribution system. Considering that the initial values of the voltage in all the nodes are known and that each one demands a complex power S i = P i + jQ i , the proposed methods consist of applying Kirchhoff's current laws and establishing a second-order nonlinear equation to calculate iteratively the voltages of each node [19].  Figure 1; V i , V j , V k and V m , represent the voltages at the nodes, i, j, k and m, in general, is represented as; V i and has units are in pu, For any line; Z ij represents the impedance between node i and node j, and their units are in pu Y shij is the parallel admittance in pu of the line ij, The complex powers demanded at each node have been represented by; S ci , Scj, S ck and S cm . In general, they can be represented as; S i = P i + jQ i and its value can be in MVA or in pu.

Proposed Methods
The symbol '*', as a superscript, means; the conjugate of: In phasor analysis the series impedance and the parallel admittance of the links are represented by their real part and imaginary part as; Z ij = R ij + jX ij and Yshij = jyωC ij , the complex power demanded at each node by; S i = P i + jQ i . Figure 1 shows in general the type of nodes present in an electrical power distribution system. The node (j) that is forming part of the main branch and at the same time is part of another branch to which it is delivering electrical power. If Kirchhoff's current law [19] is applied to this node, we have the equation (1).
Rearranged and defining terms, we have the equation (2).
To calculate the voltage at each node with the proposed methods, the one-dimensional array Y D(N) is generated that contains the series impedances and the parallel admittance of the links, a rectangular array Y nat (N, 5) , which is formed by the series impedances of the links that interconnect the corresponding node. Thus, for the node (j), we have equation (3).
To solve equation (3) with a Halley-type method, function (4) is generated, which depends on the nodal voltages ( 1 , , , ) = * + * + It is derivative with respect to the voltage of node j, and its second derivative must exist to apply the method.
Equation (5) is used to iteratively calculate the voltage at any node in the network. For node i, using Halley's formula [20].
Where; N is the number of nodes and α is an acceleration factor between 2.5 and 10. If equation (2) is solved by applying an iterative Gaussian method, we have equation (6).
The number of elements in equation (4) depends on the node that is being solved. For node i, the equation will consist of only three elements, and for node k, for two, node j and node k. In the nodes placed at the end of the branch, such as node k, applying the law of currents, we have equation (7).
Or also; The admittance Y D for the nodes of the section shown in figure 1 has the structure represented by equation (9).
The admittances, which could be called mutual, are stored in the matrix arrangement represented by equation (10).

Results
Two programs were developed in FORTRAN to iteratively solve equations (4) and (6). For the solution of equation (4), an iterative Halley-type method (Eq. (5)) is used. For the solution of equation (6), a Gauss-type method (Eq. (6)) is used. The free software FORCE-2.0 [24] is used on a 2.0GHz i7 computer and 8 Gb in RAM. It is assumed that there is only one source in the network and the voltage is kept constant throughout the iterative process, corresponding to node 1. A flat profile of the network voltages is initially considered. The stopping criterion in the iterative process considered is the relative error of the voltages and is calculated in each iteration with equation (9).
The tolerance (tol) used in the two algorithms was 1.0e-05. The two algorithms were used to calculate the voltage in each node of a system of 70 nodes (case 1) [22], in a system of 40 nodes, (case 2) [22], and in a system of 80 nodes, (case 3) [22].
A flow chart to solve equation (4) is shown in figure 2; it shows the main steps to solve equation (4) by means of the Halley-type method [20].   (6) In each block of the flow diagrams of figures 2 and 3, the activities carried out are as follows; 1. Reading data, such as; number of nodes, number of links, tolerance, base power and number of branches, elements of each branch, reactive power of each node and initial voltages, 2. Formation of the rectangular matrix Y nat and the one-dimensional arrangement Y D , 3. Calculation of the voltage with equation (5) or with equation (6), 4. Verification of the convergence, if the convergence is fulfilled, voltages are printed, it is not fulfilled, the iterative process is followed. 5. Printing voltages, 6. If desired, power flows are calculated. If not, the program ends.
In the flow charts show in Figures 2 and 3, the Y D admittances and Y nat , are calculated in the third step. If desired, once the voltage is calculated, the power flow can be calculated.

Case 1
Figure 4 corresponds to a system of 70 nodes and 6 branches, the data of impedances and demanded powers are in [21], [22], as well as the graph that shows the magnitude of the voltages and the value of the angle of each voltage.  Figure 5 shows the voltage obtained with the two algorithms used in this work, with equation (5), with equation (6), and the results obtained with a method used for interconnected power systems [23].  Figure 6 shows the angle of the voltages obtained with the two algorithms used in this work, with equation (5), with equation (6), and the results obtained with a method used for interconnected power systems [23].

Case 2
A system of 40 nodes and 20 branches is simulated with the algorithms proposed in this work. The one-line diagram, the system data, and the graphical representation of the magnitude of the voltages and their angle can be seen in [21,22,25]. Figure 7 shows the magnitude of the voltages of each node obtained with equation (5) represented in the figure with (1) and those obtained with equation (6), represented in the figure with (2).  The characteristic of a distribution system is observed in equation (5), for this reason, the one-dimensional arrangement Y D is generated that contains the impedances and admittances of the links of the node that is being solved. The impedances of the adjacent nodes to the one being solved are stored in a rectangular array that has the dimension (N, x), where; N is the size of the system and x can vary depending on the number of interconnections of a node. For the test systems, it was considered that x = 5. Solving the 70-node system with equation (6) requires a running time that ranges from 177 to 249 milliseconds. Figure 9 show the execution time with respect to the acceleration factor α for the two test systems when equation (4) is solved with the Halley-type method.
For α values greater than 2.5 in systems with 70 nodes or more, the simulation time is long. For smaller systems, the execution times are small as can be seen in the same figure 9.

Case 3
Case 3 corresponds to a system of 80 nodes and 6 branches that is derived from the system of 70 nodes [21,22]. Figure 10 shows the magnitude of the voltages in pu.   (5) is used, the execution time is 1345 milliseconds, while with equation (6) it is 374.4 milliseconds. In a distributed system, the branch end nodes are modeled by the equation (8) with the proposed algorithms. As can be seen in the same equation, there are only two nodes that generate its non-linear equation. The table 1 corresponds to the 80 node system. The table was generated to observe the effect that the acceleration factors have on the execution time. We speak of acceleration factors because in order to improve the characteristics of the Halley type method, an acceleration factor α 2 was applied to the end of branch equation and an acceleration factor α 1 to the equations corresponding to any other node Table 1 show that it is necessary to apply the acceleration factor to the end-of-branch nodes to reduce execution time when using the Halley-like method. It is observed that the application of the slowed Newton method (NR 0.9) to the end-of-branch nodes can reduce the execution time. The possibility arises here of applying other iterative methods to solve the non-linear equations that are generated. However, the Gaussian method is still the better of the two algorithms because it requires only 390 milliseconds to solve the 80-node system.

Discussion
Kirchhoff's current law was applied to a distribution system operating under balanced conditions to establish a non-linear equation at each node of the network. Once the nonlinear equation is generated, two iterative methods are used to solve the generated nonlinear equation. Due to the characteristics of the equation, two methods are applied for the solution; the Halley-type method (5) and a Gauss-type approximation (6). An advantage of the methods used is that their formulation is simple, possibly not very robust, which means that several iterations are required to reach the solution, however, the execution times, especially for the method that makes use of the equation (6) could be compared with other methods. To accelerate convergence, the technique of multistep methods was applied [23,25]. The calculation of power flows is a consequence of the calculation of the voltages, because of this, only the voltages obtained with the algorithms proposed here are shown.

Conclusions
The methods or algorithms proposed in this work were devised with the objective of having two more options for the analysis of balanced distribution systems. The application of circuit laws to an electrical energy distribution system, allowed us to establish a non-linear equation in relation to the voltage of each node of the network. This non-linear equation was solved by applying two different methods. The use of a Halley-type method was due to the fact that it is more robust than the Newton-Raphson method when systems greater than 40 nodes are analyzed, even so, the execution times are greater in relation to the Gauss-type method. The results show that, despite their differences, they can be two more alternatives to the calculation of power flows in distribution systems operating under balanced conditions.