Accurate and Efﬁcient Derivative-Free Three-Phase Power Flow Method for Unbalanced Distribution Networks

: The power ﬂow problem in three-phase unbalanced distribution networks is addressed in this research using a derivative-free numerical method based on the upper-triangular matrix. The upper-triangular matrix is obtained from the topological connection among nodes of the network (i.e., through a graph-based method). The main advantage of the proposed three-phase power ﬂow method is the possibility of working with single-, two-, and three-phase loads, including ∆ and Y -connections. The Banach ﬁxed-point theorem for loads with Y -connection helps ensure the convergence of the upper-triangular power ﬂow method based an impedance-like equivalent matrix. Numerical results in three-phase systems with 8, 25, and 37 nodes demonstrate the effectiveness and computational efﬁciency of the proposed three-phase power ﬂow formulation compared to the classical three-phase backward/forward method and the implementation of the power ﬂow problem in the DigSILENT software. Comparisons with the backward/forward method demonstrate that the proposed approach is 47.01%, 47.98%, and 36.96% faster in terms of processing times by employing the same number of iterations as when evaluated in the 8-, 25-, and 37-bus systems, re-spectively. An application of the Chu-Beasley genetic algorithm using a leader–follower optimization approach is applied to the phase-balancing problem utilizing the proposed power ﬂow in the follower stage. Numerical results present optimal solutions with processing times lower than 5 s, which conﬁrms its applicability in large-scale optimization problems employing embedding master–slave optimization structures.


General Context
Electrical distribution networks are responsible for interconnecting transmission and sub-transmission grids with final electricity users living in urban and rural areas through substations, conforming the well-known medium-and low-voltage networks [1,2]. The main parameter to analyze these networks is power flow. Thus, the power flow problem constitutes a steady-state study where electrical variables (voltages and currents) are determined for a particular load and generation condition [3,4]. The main characteristic of the power flow problem is the conformation of a non-linear system of equations as a result of the presence of products of variables and trigonometric functions, which demands the use of numerical methods to find a solution [5,6]. Typically, the power flow problem in electrical distribution grids has been addressed only with single-phase equivalents in the literature [7]. However, some electrical distribution networks cannot be reduced to single-phase equivalents due to the following characteristics: (i) no transposition in distribution lines, i.e., unbalanced impedances among phases, (ii) unbalanced loads with ∆and Y-connections, (iii) single-phase or two-phase laterals [8][9][10]. Thus, distribution networks must be analyzed directly in their three-phase form in order to capture the actual effect of the imbalances in the electrical variables, including voltages, currents, and powers.

Motivation
Powerful methodologies are always required for the analysis of electrical grids in planning and operational aspects to estimate (calculate) the electrical variables for particular load conditions or variable load cases considering real-time operation scenarios [11]. In the case of three-phase grids, the unbalanced nature of the loads and impedances makes it necessary to propose efficient numerical methods in order to deal with the power flow problem by considering single-, two-, and three-phase loads, including ∆and Y-connections [12]. Due to the importance of developing efficient tools for three-phase networks, this research is motivated by the possibility of improving the upper-triangular power flow method for unbalanced distribution networks with the inclusion of the convergence analysis via the Banach fixed-point theorem, with its main advantage being that it is derivative-free. This is important since this will facilitate its implementation using any programming language with simple operations and without variant matrices, which assists in reducing the required processing times for solving the three-phase power flow problem in unbalanced distribution networks. Additionally, we present the general algorithm to handle single-, two-, and three-phase loads with ∆and Y-connections but without modifying the power flow formulation, as the main advantage of our proposed approach is that it does not depends on the reactance/resistance ratio.

Literature Review
Multiple power flow methods are proposed in the existing literature for the analysis of radial three-phase distribution networks. Authors in [13] presented a matrix reformulation of the classical three-phase backward/forward power flow method for balanced grids. Even if this method is classical [14], the authors demonstrated its convergence using the Banach fixed-point theorem. An improved version for unbalanced three-phase grids that consider loads connected in ∆and Yhas been proposed in [15] to solve the phasebalancing problem by employing a leader-follower optimization approach. The power flow numerical results were compared with the DigSILENT software (three-phase Newton-Raphson method), demonstrating its effectiveness in calculating power losses. Authors in [8] proposed a decoupled power flow method using a quasi-symmetric matrix formed as a result of the inductive and resistive effects of the distribution lines. This method works in real domains by separating the real and imaginary parts of the voltages. Numerical results demonstrate its computational efficiency when compared with the classical yet complex backward/forward formulation in different IEEE test feeders. Garces in [4] proposed a linear formulation for the power flow problem in three-phase grids. This formulation was based on Laurent's series expansion of the hyperbolic relation among voltage and power balance equations. Numerical results contain errors lower than 3% for different IEEE test feeders when compared with the backward/forward power flow method. Authors in [16] proposed a holomorphic embedding power flow approach for determining the electrical variables in medium-and low-voltage three-phase grids. Numerical results demonstrated the efficiency of this approach when compared with the classical three-phase Newton-Raphson approach. In Sereeter, [17], extended versions of the classical Newton-Raphson power flow method have been compared with the backward/forward load flow approach. The algorithms' convergence for different loading conditions, resistance/reactance ratios, and load models was studied by performing numerical experiments on balanced and unbalanced distribution grids. Authors in [18] presented an iterative methodology for the power flow solution in three-phase distribution networks based on decoupled circuit equivalents in the complex domain. Numerical comparisons with classical tools (such as Newton-Raphson in specialized power system software) demonstrated its effectiveness in large-scale distribution grids. Marini et al. [19] proposed a graph-based power flow method for radial and meshed distribution networks using an upper-triangular formulation. This paper's main contribution is to the three-phase modeling of regulators, transformers, and induction machines, among others. Numerical results demonstrate the computational efficiency of the proposed power flow method for different reactance/resistance ratios. Authors of [20] presented a three-phase power flow approach for islanded microgrids that is based on the Newton-Raphson current injection method in order to deal with the lack of a slack node. The proposed algorithm selects system frequency and voltage magnitude of the reference bus as additional variables in power-flow formulation. Numerical simulations confirm the efficiency of this method for different reactance/resistance ratios, including the best processing time performance, compared to classical Newton-Raphson methods and specialized tools available in the PSCAD software.
In order to summarize the three-phase power flow approaches reported in the literature, Table 1 presents some approaches applicable to power distribution grids with unbalanced load conditions.

Solution Methodology
Year Ref.

Contribution and Scope
This paper introduces a modified version of the upper-triangular-based power flow method originally reported in Marini et al. [19]. The proposed formulation considers three-phase unbalanced distribution systems with loads connected to ∆ and Y using a matrix representation. The modified power flow is tested on three different networks with 8, 25, and 37 nodes to demonstrate the algorithm's accuracy and efficiency. The main advantages of the proposed upper-triangular power flow method are as follows: (i) the possibility of working with single-, two-, and three-phase loads with ∆and Y-connections; (ii) independence in the reactance/resistance relations in the distribution lines; and (iii) the processing time performance when compared with the classical backward/forward power flow method.
The computational performance of the triangular-based power flow is assessed within an optimization problem solved through metaheuristics. In this paper, power flow is evaluated using the classical Chu-Beasley genetic algorithm (CBGA) to solve the phasebalancing problem by employing a leader-follower approach. The leader stage is entrusted with defining the configuration of the phases, while the triangular-based power flow method determines the total power losses for each phase configuration defined in the leader stage. The main contributions of this paper to the current literature are as follows: • The analysis of the upper-triangular-based power flow under different load connections, i.e., ∆and Y-connections. • A convergence analysis based on the Banach fixed-point theorem for three-phase networks with pure Y load connections. • Application of the proposed upper-triangular three-phase power flow method to solve the phase-balancing problem using a CBGA.
Note that these contributions are concentrated in the general power flow formulation for three-phase distribution networks with unbalanced loads and radial configurations; however, the analysis of convergence is only developed for Y-connected loads due to the possibility of generalizing its recursive formula with diagonal matrices and vectors, which is not possible in the case of ∆-connected loads. In this sense, the study of the convergence for ∆ loads could be explored in future research.

Structure of the Document
The remainder of this research is structured as follows: Section 2 presents the singlephase formulation of the triangular-based power flow method using a small example to explain the main concepts of the non-derivative power flow formulation. Section 3 presents the three-phase extension of the triangular-based power flow method with a consideration of different load configurations, i.e., Yand ∆-connections. Moreover, Section 4 presents the convergence test using the Banach fixed-point theorem for three-phase loads with solidly grounded Y-connection. Section 5 presents the information pertaining to test systems; these test feeders are composed of 8, 25, and 37 nodes. Section 6 presents the numerical results for the test feeder under different load connections, highlighting comparisons with the classical Newton-Raphson and backward/forward methods. Furthermore, the proposed variable of power flow is tested by employing a CBGA leader-follower optimization structure in order to solve the phase-balancing problem. Finally, Section 7 presents the main conclusions obtained from this research and discusses the possible areas for future work.

Single-Phase Power Flow Formulation
The three-phase power flow problem corresponds to the steady-state formulation of electrical distribution networks with unbalanced loads and different feeder configurations. The idea behind solving this problem is to determine the network state variables, i.e., voltage magnitudes and phase angles. In this paper, we adopted a graph-based representation of electrical distribution grids based on the upper-triangular matrix formulation (it is recommended that the grid is ordered starting from the slack bus). Consider the schematic representation of the 7-node, 6-branch distribution network depicted in Figure 1 with a single-phase representation. Note that branch currents J j can be expressed as a function of the demanded currents I i using the following expression: Here, J j represents the current through the branch j, and I i corresponds to the net injected current at node i, which can be zero in the case of the step-nodes. A general expression can be compacted into Here, J ∈ C b×1 is the vector that contains all the branch currents in the complex domain (b branches); I ∈ C (n−1)×1 represents the vector that contains all the demanded currents in the nodes (n) of the network, except in the slack source; T ∈ R b×b represents the upper-triangular matrix.
Consider that V i represents the voltages at load buses, V 1 is the voltage at the slack bus, and E j is the branch voltage drop. Applying Kirchhoffs' second law to each trajectory that contains the slack source, i.e., V 1 , and all the nodes of the network, the following expression is obtained: This can be expressed in a compact form as: Here, E ∈ C b×1 is a vector that contains all branch voltage drops, V ∈ C (n−1)×1 represents the vector containing the voltages in demand nodes, and V 1 is the voltage at the substation bus. Note that 1 is a vector filled by ones.
Branch voltage drops can be expressed as a function of the branch current after applying Ohm's law. This is performed by using the primitive impedance matrix Z, which has been presented below: Here, Z ∈ C b×b is a complex diagonal matrix containing the branches' impedance, i.e., Z = diag([Z 1 , Z 2 , · · · , Z 6 ]).
To obtain a general formula for nodal voltages V and the nodal currents I, (3) and (1) are substituted into (2), producing the following expression: Note that (4) can be applied in single-phase distribution grids, where the nodal injected current is calculated as follows: Here, S ∈ C (n−1)×1 is a complex vector that contains all constant power demands, and stands for the vector's component conjugate. Now, if (5) is substituted into (4), then the following recursive formula is yielded: Remark 1. Equation (4) relates nodal voltages and injected currents to demand nodes. However, note that constant power demands need to be expressed as current injections. Hence, there is a need to obtain a recursive formula for determining the voltages at demand nodes as a function of constant power demands.

Three-Phase Power Flow Formulation
The general formulation of the triangular-based power flow method in three-phase systems needs to consider three possible cases, i.e., power consumption connected in (i) Y-connection, (ii) ∆-connection, and (iii) both load configurations.

Three-Phase Grid with Demands in Y-Connection
Equation (4) is transformed into a three-phase expression as follows: Here, V 3ϕ ∈ C 3(n−1)×1 is a vector that contains all the phase voltages ordered per node; 1 3ϕ ∈ R 3(n−1)×3 is a rectangular matrix filled by 3 × 3 identity matrices; V 13ϕ ∈ C 3×1 is a vector that contains the voltages of the substation node (fixed); T 3ϕ ∈ R 3b×3b represents the three-phase equivalent of the upper-triangular matrix T, where each connection (1) is converted into a 3 × 3 identity matrix (note that, in case the grid has two-phase or singlephase feeders, the corresponding element in this matrix takes the value zero), and each zero becomes a 3 × 3 zeros matrix; Z 3ϕ ∈ C 3b×3b is the primitive three-phase impedance matrix, which has a three-diagonal structure (note that three-diagonal matrix form corresponds to a matrix composed of 3 × 3 matrices diagonally, with the remainder of the terms being equal to zero); I 3ϕ ∈ C 3(n−1)×1 is a vector that contains all the phase currents ordered per node.
For the purpose of illustration, consider the case where all nodes have constant power loads in the Y-connection and are solidly grounded. Take the generic Y-connection for a load at node k as depicted in Figure 2. In the three-phase load with Y-connection, we assume that the voltage experiencing each load is the phase-to-neutral; with this assumption, each line-current is calculated as follows: This can be compacted as follows: Considering the diagonal structure of the three-phase current at node k, it is possible to generalize the triangular three-phase power flow formulation for three-phase loads connected in Y as follows: Remark 2. To solve the upper-triangular three-phase power flow formulation, (9) is added as an iterative counter t, starting from the initial voltages defined as V 0 3ϕ = 1 3ϕ V 13ϕ . The recursive formulation of the three-phase power flow problem takes the following form: Note that Equation (10) represents a non-derivative power flow method for three-phase distribution grids with balanced or unbalanced structures. The main advantage of this formulation is its simplicity since the impedance-like matrix (i.e., Z bus 3ϕ = T T 3ϕ Z 3ϕ T 3ϕ ) is calculated one time and stored to improve the processing time performance of the power flow solution.

Three-Phase Grid with Demands in ∆-Connection
To obtain the line currents for a three-phase network with loads connected in ∆, consider the load diagram presented in Figure 3 for the k node.
This can be compacted into the following: Here, the matrices M and H take the following form: To solve the problem of the power flow in three-phase network with ∆ loads, the following steps of the iterative procedure were developed: 1. Define the test feeder characteristics and calculate the three-phase impedance-like matrix, i.e., Z bus 3ϕ ; 2. Define the starting point of the voltages using the substation voltage; 3. Calculate the current for each node k using (9) for the iteration t as follows: 4. Calculate all the voltages in the demand nodes using (7) as presented below: 5. Determine whether the convergence error has been fulfilled-if yes, the power flow problem is solved; otherwise, update the starting voltage point and return to step 3.

Three-Phase Grid with Demands in ∆and Y-Connections
In the case that a three-phase distribution grid has loads connected with Y and ∆ structures, the calculation of the demanded currents must consider these connections as presented in Equations (8) and (11), respectively. Algorithm 1 presents the general implementation of the three-phase triangular power flow for unbalanced distribution networks with loads connected in Y and ∆.

Power Losses Calculation
One of the main objectives of the power solution in single-and three-phase grids is computing the amount of power losses in all the conductors of the network after determining all the voltages in the demand nodes. To perform this task, since we know that V 3ϕ and I 3ϕ represent the final voltages and the final demanded currents reported by the power flow solution, which fulfills the convergence criteria, from Equations (1) and (3) in their three-phase forms, we know the following: When these are combined, the following is produced: Here, S loss represents the apparent power losses in the network.

Algorithm 1:
Pseudo code for the solution of the three-phase power flow problem in unbalanced distribution networks with Y and ∆ loads. Data: Define the three-phase grid under study.
Transform the system into a per-unit equivalent; Calculate the three-phase upper-triangular matrix T 3ϕ ; Calculate the primitive branch impedance matrix Z 3ϕ ; Compute the three-phase impedance-like matrix Z bus 3ϕ = T T 3ϕ Z 3ϕ T 3ϕ ; Define the maximum number of iterations t max ; Define the convergence error ; Define the substation voltages: if Load in node k connected to Y then Compute the demanded current I k3ϕ using Equation (8) ; else Compute the demanded current I k3ϕ using Equation (11); end end Calculate the new voltages V t+1 3 f using Equation (7); Report the nodal voltages as V = V 13ϕ ; V t+1 3ϕ ; Report the final three-phase current, i.e., I 3ϕ ; break; else Make V t 3ϕ = V t+1 3ϕ ; end end

Convergence Analysis for Y-Connected Loads
The main advantage of the upper-triangular three-phase power flow formulation is the possibility of ensuring convergence when all the three-phase loads have a Y-connection. Here, we present the convergence test based on the Banach fixed-point theorem. The convergence analysis presented in this section is based on the convergence test for the successive approximation method used in [7] for single-phase distribution networks. The convergence test is based on the following assumptions [27]: Assumption 1. The total consumption of the active active power in the distribution does not cause a voltage collapse, which implies that the power flow equations can be solved.

Assumption 2.
The lower bound for all the voltages is positive, i.e., V min > 0, which is fulfilled by the grid since regulatory policies have made it mandatory. Assumption 3. The impedance-like matrix Z 3ϕ is diagonal dominant, suggesting that Z 3ϕ jj ≥ Z 3ϕ jk , ∀j = k is always ensured.
To present the convergence properties of the triangular-based power flow method for three-phase networks with the structure defined in (10), let us consider the general definition of the Banach fixed-point theorem as follows [13,27]: Fixed-point theorem). The recursive formula defined by Equation (10) is stable and defines a contraction map with the following form: for some V 3ϕ that fulfills Assumption 1 and independent at the starting point V 0 3ϕ , such that: Here, U defines the solution of the three-phase power flow problem, i.e., a vector that contains all the voltages of the demand nodes. In addition, η corresponds to a real number between [0, 1].
Proof. Note that the iterative expression that defines the triangular-based power method showed in (10) can be rewritten as follows: i3ϕ T i∈D (15) with D being the set that contains all the demand buses ordered per node and phase.
Owing to the structure of the three-phase power flow problem, it is possible to say that the solution U complies with U = f (U), and it defines a unique solution, if and only if f (U) corresponds to a contraction mapping on V 3ϕ .
Now, if we consider Assumption 3, relating the nature of the impedance-like matrix, Equation (17) can be transformed into the following: Note that, considering the structure of (18) and taking into account that Z 3ϕjj defines the equivalence impedance at the node j (i.e., the Thévenin equivalent impedance), the following relation is met: Here, 0 ≤ η ≤ 1 can be ensured since the denominator of (19) can be interpreted as the lowest short-circuit current, while the numerator denotes the maximum loading current. At the same time, it is known that the maximum loading current is lower than or equal to the lowest short-circuit current for any demand scenario. This characteristic of the η-coefficient allows confirmation of the recursive power flow formula (10) in order to make sure the three-phase power flow problem is stable and converges to the solution U, which completes the proof.

Test Feeders
To validate the proposed three-phase power flow method, let us consider three test feeders composed of 8, 25, and 37 nodes, which have been reported in [15] for phasebalancing studies.

8-Bus Test Feeder
The eight-bus test feeder corresponds to a radial distribution network composed of eight nodes and seven lines, where the slack node is located at bus 1 with a line-to-line voltage set at 11 kV. The electrical configuration of this test system is depicted in Figure 4. In this eight-bus system, the total active and reactive power demands are 1005 kW and 485 kvar for phase a, 785 kW and 381 kvar for phase b, and 1696 kW and 821 kvar for phase c. The electrical information of the impedances and loads for this test system are presented in Tables 2 and 3, respectively.

Table 2. Loads and lines parameters of the 8-node test system (all powers in kW and kvar).
Line Node i Node j Cond. Length (mi) P ja Q ja P jb Q jb P jc Q jc

25-Bus Test System
The 25-bus system is a radial distribution network with 25 nodes, 24 lines, and 22 constant power loads. The voltage controlled source is located at bus 1, which operates with line-to-line voltage of 4.16 kV. This is a radial unbalanced distribution system with 25 nodes, 24 lines, and 22 loads. The source is located at node 1, and the nominal voltage is 4.16 kV [28]. The grid configuration of this test system is presented in Figure 5. The total loads for phases a, b, and c are 1073 kW and 792 kvar, 1083.3 kW and 801 kvar, and 1083.3 kW and 800 kvar, respectively. The electrical information of the impedances and loads for this test system are presented in Tables 4 and 5, respectively. This parametric information was taken from [29].

IEEE 37-Bus System
This test feeder corresponds to a underground system located in the state of California, which is composed of 37 nodes and 36 lines (see its electrical configuration in Figure 6). The voltage controlled source is connected at node 1 operating with a line-to-line voltage of 4.8 kV. In this system, there are 25 constant power loads, with a total and reactive power consumption in phases a, b, and c of 727 kW and 357 kvar, 639 kW and 314 kvar, and 1091 kW and 530 kvar, respectively. The IEEE 37-bus test system used in this research corresponds to the adaptation proposed in [15], where the voltage regulator connected between nodes 1 and 2 is replaced by a single three-phase line with a length of 1850 feet, and the transformer located at nodes 10 and 24 is removed, including node 24, as recommended in [30].  Tables 6 and 7. These data were taken from [15]. Table 6. Loads and lines parameters of the IEEE 37-bus test system (all powers in kW and kvar).

Computational Validation
In this section, we present the numerical validation of the proposed non-derivative power flow method for three-phase unbalanced distribution grids with loads connected in ∆ and Y. The following cases were considered for simulation purposes: (i) evaluation of the power flow assuming that all the loads are connected to Y, comparing its results with the backward/forward approach reported in [15] and the classical Newton-Raphson method available in the DigSILENT software; (ii) solution of the three-phase power flow considering all the loads connected in ∆; (iii) solution of the phase-balancing problem using a leader-follower approach based on the CBGA in the leader stage and the proposed non-derivative triangular-based power flow method in the follower stage.

Power Flow Solution with Y-Connected Loads
To validate the effectiveness of the non-derivative triangular-based power flow method for solving the power flow problem in three-phase grids, we consider the three test feeders described in the previous section. In addition, the backward/forward and proposed methods are evaluated 10,000 consecutive times each to determine their respective average processing times. Table 8 presents the numerical values indicating the performance of each power flow method. It is important to mention that the processing time of the Newton-Raphson is not reported as this method is implemented in the DigSILENT software (with a different programming environment than the MATLAB used for the backward/forward and the proposed power flow approaches), and the minimum convergence error available is about 1 × 10 −6 . Table 8. Behavior of the power flow methods considering pure-Y load connections.

Method
Mean (ms) Losses (kVA) Iter. Note that the results reported in Table 8 allow observation of the following:

8-bus test feeder
The Newton-Raphson method in the DigSILENT software generates a solution that is pretty similar to the backward/forward and the triangular-based power flow methods, albeit with a few small differences in the decimals. These differences are attributable to the convergence error in that, in DigSiLENT, 1 × 10 −6 can be set as the minimum value, while that for the backward/forward and the triangular methods are set at 1 × 10 −10 .
The number of iterations of the Newton-Raphson approach will always be lower than the backward/forward and triangular-based method. This is an expected since the Newton-Raphson method uses the information of the derivatives in the Jacobian matrix, which provides additional information about the direction of the solution, while the backward/forward and the triangular-based methods are non-derivative and only make some additional iterations if they are necessary to solve the power flow problem. It can be found that the proposed triangular-based power flow method is faster than the backward/forward power flow method for all the test feeders, with the former showing improved processing times of 47.01%, 47.98%, and 36.96% for the 8-, 25and IEEE 37-bus systems, respectively. This general improvement in the processing times results from the fact that the triangular-based power flow method does not use any inverse in the power flow solution, while the backward/forward, based on the incidence matrix, needs inverse matricial operations to determine the nodal admittance matrix, requiring additional processing time.

Power Flow Solution with ∆-Connected Loads
To validate the applicability of the triangular-based power flow method to distribution grids with loads connected in ∆, here, we evaluate a simulation case in which all the loads of the three test feeders are connected in this form (i.e., ∆-connection). The numerical results of this simulation are reported in Table 9. The main findings from this table are as follows: (i) the triangular-based power flow method is superior in terms of processing times with improved performances of 43.35%, 60.14%, and 49.63% for the 8-, 25-, and IEEE 37-bus systems, respectively; (ii) the estimation of power losses by both methods is equal in all the test feeders as well as the number of iterations. It is important to highlight that, compared to the simulation case with loads connected in Y, the total power losses reduce when loads are connected in ∆. These reductions are 21.10%, 2.65%, and 14.40%. In addition, note that no comparison results with the Newton-Raphson are presented, since the DigSILENT software does not allow constant power loads with ∆-connection.

Solution of the Phase-Balancing Problem
In this section, a possible application of the studied upper-triangular three-phase power flow embedded in the metaheuristic optimizer is presented, thereby generating a leader-follower optimization approach. Here, the classical phase-balancing problem in unbalanced distribution networks using the well known CBGA is explored to determine the best phase configuration in the leader stage, while the three-phase power flow is used in the follower stage to determine the total grid power losses.
Note that the main idea of presenting the phase-balancing problem solved with the classical CBGA and the proposed upper-triangular-based power flow method is to show the possibility of using this in master-slave optimization approaches to formulate the power flow solution with short processing times, especially when the three-phase power flow method has solved the problem multiple times. However, it is important to clarify that the power flow approach does not improve the optimal solution of the phase-balancing itself since this solution is dependent on the efficiency of the master optimization stage and not the power flow method used in the slave stage.
The codification adopted for the CBGA is composed of a vector of integral values with dimension 1 × n, i.e., [1, 5, 6, 2, · · · , 3], where the numbers between one and six represent a phase connection. In Table 10, the meaning of the proposed codification is described. Note that, initially, we assumed that all the loads are connected in ABC form, i.e., coded by one. The main characteristics of the implementation of the CBGA to solve the problem of the phase-balancing are presented in Algorithm 2. Note that, in this simulation case, we consider that all the loads have an Y-connection.

Algorithm 2:
Application of the CBGA to the phase-balancing problem.
Data: Define the three-phase grid under study. Determine the population size; Define the number of iterations, i.e., m max ; Generate the initial population; Verify the diversity criteria in the initial population; Evaluate the initial population using Algorithm 1; Arrange the initial population in ascending order; Identify the worst individual in the population, worst m ; for m ≥ m max do Select two individuals randomly from the population; Recombine these individuals; Apply the mutation operator for both of these individuals; Evaluate these individuals using Algorithm 1; Select the winner individual, i.e., winner m ; if p loss (winner m ) < p loss (worst m ) then Verify the diversity criteria ; If diversity criteria is met, then replace the worst individual with the suitable one; end end Report the best solution in the final population; The parameters implemented in the CBGA are as follows: (i) 50 individuals in the population; (ii) 1000 iterations per evaluation; (iii) 100 consecutive evaluations. Numerical results with the three test feeders are summarized in Table 11. The interpretation of the columns in Table 11 from left to right is as follows: test feeders, minimum, mean, maximum, and standard deviation of the power losses of the network; and the total processing time after solving the phase-balancing problem with the classical CBGA. From results in Table 11, the following can be noticed: The processing times to solve the phase-balancing problem in three-phase distribution networks in all the simulation cases is less than 5 seconds, which can be considered a really fast processing time considering that the solution spaces have dimensions of 279,936, 2.8430 × 10 19 , and 6.1887 × 10 28 for each test feeder. The optimal solution obtained with the CBGA on each test feeder is better than that reported in [15] for the 25-and 37-node test feeders and is the same for the 8-bus test system. Notice that the solutions obtained (power losses) in the current paper are 10.5869 kW, 72.2889 kW, and 61.4923 kW for the 8-, 25-, and 37-node systems, respectively, while solutions reported in [15] are 10.5869 kW, 72.2919 kW, 61.5785 kW, respectively. It is important to mention that the improvements reported in this approach can be attributed to the number of iterations and the population size considered in our proposal as compared to the CBGA, reported according to [15]. The standard deviations for the CBGA in Table 11 show that the solution space and the differences between the maximums and minimums have a direct relationship. However, power losses lower than 0.13 kW in the IEEE 37-bus test feeder demonstrate the repeatability properties of the CBGA to solve the phase-balancing problem even considering the dimensions of the solution space, i.e., greater than 1 × 10 28 . When compared to the base cases of each test feeder (see Table 8), reductions in power losses resulting from the CBGA approach in the 8-, 25-, and IEEE 37-bus test systems were 24.34%, 4.15%, and 19.23%, respectively. Reduction in technical losses can be converted into cost reduction by the utility company with minimum inversion efforts since phase balancing does not require special devices, and it only requires a working crew to go over through the feeders to modify the load connections. However, to calculate the net savings that the utility can reach when applied to the phase-balancing problem in its grids, it is necessary to subtract the cost associated with the crew workers of the objective function.

Conclusions
The three-phase power flow problem with Yand ∆-loads considering imbalances is addressed in this research with the derivative-free upper-triangular power flow formulation in the complex domain. The proposed approach also has the ability of working in the three-phase distribution networks with single-, two-, and three-phase loads, offering the main advantage that no assumptions regarding reactance/resistance ratio are required. Numerical comparisons demonstrate that the proposed power flow method reaches the same numerical solution, i.e., power losses, compared to the classical Newton-Raphson and backward/forward power flow methods. Moreover, was found that the proposed power flow converges in the same number of iterations as those in the backward/forward power flow method. Regarding processing times, it was found that the proposed upper-triangularbased power flow method is between 36% and 40% faster than the backward/forward power flow for Yand ∆-loads. Due to its computational efficiency, the proposed power flow method can be used in problems in which the three-phase power flow must be solved recursively; namely, the phase-balancing problem and a plethora of other power systems problems, for example, probabilistic analyses, planning problems, contingency analysis, among others.
A CBGA is used to solve the phase-balancing problem in electric distribution networks. The application was made to demonstrate the efficiency of the proposed three-phase power flow approach and tested on radial systems composed of 8, 25, and 37 nodes. Numerical results show that the computational times required to solve the problem in the three test feeders was 0.7612 s, 3.5526 s and 4.7481 s, which are outstanding, considering that the solution space dimensions for the test feeders are 27.9936 × 10 4 , 2.8430 × 10 19 , and 6.1887 × 10 28 , respectively. Moreover, each test feeder's final objective function value is better than those previously reported in the literature where a CBGA and the backward/forward power flow method were used.
For future research, the paper makes many contributions, such that it will be possible to: (i) implement the studied three-phase power flow method in a multi-period environment to locate and determine the size of distributed generators and capacitors banks embedded in a metaheuristic optimization approach; (ii) extend the three-phase power flow problem to distribution networks with a neutral conductor and meshed configurations, and (iii) include the three-phase power flow formulation voltage controlled nodes.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.