Improved Genetic Algorithm for Phase-Balancing in Three-Phase Distribution Networks: A Master-Slave Optimization Approach

: This paper addresses the phase-balancing problem in three-phase power grids with the radial conﬁguration from the perspective of master–slave optimization. The master stage corresponds to an improved version of the Chu and Beasley genetic algorithm, which is based on the multi-point mutation operator and the generation of solutions using a Gaussian normal distribution based on the exploration and exploitation schemes of the vortex search algorithm. The master stage is entrusted with determining the conﬁguration of the phases by using an integer codiﬁcation. In the slave stage, a power ﬂow for imbalanced distribution grids based on the three-phase version of the successive approximation method was used to determine the costs of daily energy losses. The objective of the optimization model is to minimize the annual operative costs of the network by considering the daily active and reactive power curves. Numerical results from a modiﬁed version of the IEEE 37-node test feeder demonstrate that it is possible to reduce the annual operative costs of the network by approximately 20% by using optimal load balancing. In addition, numerical results demonstrated that the improved version of the CBGA is at least three times faster than the classical CBGA, this was obtained in the peak load case for a test feeder composed of 15 nodes; also, the improved version of the CBGA was nineteen times faster than the vortex search algorithm. Other comparisons with the sine–cosine algorithm and the black hole optimizer conﬁrmed the efﬁciency of the proposed optimization method regarding running time and objective function values. than 2 kW; this shows that the phase-balancing approach using the improved CBGA effectively redistributes the loads in all the phases as uniformly as possible.


Introduction
Three-phase distribution grids are responsible for providing electrical energy to residential, commercial, and industrial areas at medium-and low-voltage levels [1,2]. These grids are typically unbalanced owing to the following factors: (i) line configurations are non-symmetric because the transposition criteria are not applicable owing to the short length of the lines (few tens of kilometers) [3]; and (ii) the nature of loads can be single-, two-, and three-phase, which results in intrinsic imbalances in the currents through the lines and voltages in the nodes [4]. The main problem arising from the general unbalanced nature of three-phase networks is the large power losses in the distribution stage (distribution lines and transformers); these losses cause increments in the final billing of the electricity service for all the end-users, which is permitted by regulatory entities within an acceptable range [2]. In Colombia, the power losses of distribution systems range from 6% to 18%, depending on the geographic location of the distribution systems (i.e., rural or urban areas) as well as the level of investment in maintenance activities on the physical layer of the distribution network [5]. The regulatory entity of the electrical distribution service in Colombia, that is, the Energy and Gas Regulation Commission (CREG, based on its acronym in Spanish) allows transferring part of the energy losses to the end-users through the final billing for the electricity service; however, this value cannot exceed 8% of the total energy costs. This implies that, if the distribution company reduces the amount of energy losses to less than 8%, the following benefits can be achieved: (i) an improvement in the potential of the network to supply additional users with the same distribution grid and minimum investment costs (expansion in coverage capacity); (ii) additional economic reedits because lower energy input is needed for supplying the same group of users; and (iii) the possibility to win by billing the difference between the upper regulation bound (i.e., 8%) and the current power losses [6].
Owing to the importance of strategies for reducing the amount of energy loss in electrical distribution networks, multiple well-known approaches have been established. Some of these include (i) optimal sizing and installation of distributed generation [7]; (ii) optimal installation and sizing of reactive power compensators [5]; (iii) optimal grid reconfiguration [8]; and (iv) optimal phase-balancing strategies [2]. Each of the aforementioned methodologies can help distribution companies in significantly reducing the amount of power loss. However, the first two strategies entail large amounts of inversions due to the necessity of integrating new devices into the distribution grid, which implies that the return rate is in the order of years (5 to 15 years) [9]. The third approach, which is based on reconfiguration, requires less investment because few tie lines need to be constructed for the optimal reconfiguration of the grid [10]. The fourth approach is the most economical strategy for reducing power and energy losses in distribution grids, because it only requires a few work crews to reconfigure the load phases without investing in additional devices [11]. Here, considering the low costs associated with the phase-balancing strategy for reducing the amount of energy loss in electrical distribution networks, we propose a new master-slave optimization strategy to address this problem.
As indicated by existing literature, the problem of phase balancing has been solved using multiple optimization approaches, including the classical Chu and Beasley genetic algorithms [12][13][14][15]; tabu search algorithm [16,17]; ant colony optimization [18,19]; immune optimization algorithm [20]; branch and bound and convex methods [4,21]; bat optimization algorithm [22]; vortex search algorithm [2]; particle swarm optimization methods [8,23,24]; differential evolution algorithm [25]; and simulated annealing optimization method [26]. The main characteristics of these approaches are the hybridization of metaheuristic discrete optimization methods (with binary and integer codifications) with three-phase power flow methods, which are typically based on sweep iterative backward/forward power flow methods, and the minimization of the amount of power loss for a particular load condition, which typically corresponds to the peak load condition.
Based on the aforementioned state-of-the-art strategies, in this paper, we propose a master-slave optimization approach based on an improved version of the classical Chu and Beasley genetic algorithm (CBGA) in the master stage, in conjunction with using the three-phase version of the successive approximation power flow method in the slave stage. The master stage is entrusted by defining the phase load configuration using an integer codification between 1 and 6, which represents the six possible methods for connecting a three-phase load. The slave stage is entrusted with the evaluation of the multiperiod power flow problem, in order to determine the amount of power loss in each nodal phase configuration. Improvements in the classical CBGA are realized in the selection, mutation, and recombination procedures based on a probability criterion. If the probability is less than 50%, two individuals are selected in the tournament stage, which are then recombined using a single point. Subsequently, an adaptive multi-point mutation strategy is employed to increase the rate of convergence of the algorithm; however, if the probability exceeds 50%, an offspring is generated using a normal Gaussian distribution probability, which is employed in the evolution process of the vortex search algorithm approach [27]. The main advantage of the vortex search algorithm is the exploration and exploitation of a solution space by using adaptive hyper-ellipses based on one individual solution selected from the current population. It should be noted that this evolution criterion of the proposed CBGA is applied at each iteration; this implies that the evolution process of the algorithm is adaptive.
The main contributions of this research are as follows: (i) we propose an improved version of the CBGA based on adaptive evolution criteria by combining the classical selection, recombination, and mutation stages with the hyper-ellipses used in the vortex search algorithm, in order to explore and exploit the solution space; (ii) the solution of the phase-balancing problem in distribution grids is determined considering daily and reactive power load curves; and (iii) a three-phase version of the successive approximation power flow method, which is suitable for loads connected in Y and ∆, is formulated.
It is important to highlight that the proposed improved version of the classical CBGA with the vortex search algorithm (VSA) increases the possibility of enhancing the exploration and exploitation of the solution space. This is made by introducing normal Gaussian distributions to sweep some promissory solution zones. Hence, this new version of the CBGA with the VSA method can be considered as a powered metaheuristic optimization technique that can be extended to multiple MINLP problems in future researches [27]. Furthermore, an additional important contribution presented in this study corresponds to the extension of the recently developed successive approximation power flow method [28] to three-phase unbalanced distribution networks. The resulting method could be applied to loads with ∆− and Y−connections including possible meshed grid configurations with the possibility of demonstrating its convergence by applying the Banach fixed-point theorem as recently proposed in [29] for the upper-triangular power flow method.
The remainder of the paper is structured as follows: Section 2 presents the mathematical formulation of the phase-balancing problem in three-phase distribution grids by highlighting the mixed-integer nonlinear structure of the optimization problem. Section 3 presents the master-slave optimization approach based on the improved version of the CBGA and the three-phase version of the successive approximation power flow method. Section 4 describes the main characteristics of the test feeders comprising 15 and 37 nodes. Section 5 presents the numerical validation of the proposed optimization approach for the IEEE 37-node system considering the daily load curves. The peak load scenario operative scenario is also compared with the classical CBGA, the black hole optimizer, the sine-cosine algorithm and the vortex search algorithm. Section 6 presents the main concluding remarks derived from this research as well as some possible future works.

Formulation of the Phase-Balancing Problem
The phase-balancing problem in three-phase imbalanced distribution networks can be formulated as a mixed-integer nonlinear programming (MINLP) model, in which the integer variables are related to a combination of the load phases. Moreover, the nonlinearities are produced by the power balance equations, which yield products among voltage magnitudes and trigonometric functions. Here, we present the objective function and a set of constraints that represent the phase-balancing problem, considering multiperiod analyses.

Objective Function
The objective function considered for the phase-balancing problem corresponds to the minimization of costs for the daily energy loss at all the conductors in the grid. This objective function can be expressed as follows: where f 1 represents the annual operative costs associated with the total energy losses in all the branches of the network; C kWh is the average cost of energy; Y k f mg represents the magnitude of the admittance that relates node k at phase f with node m at phase g; V k f t (V m f t ) corresponds to the voltage magnitude at node k (m) in phase f (g) at time period t; δ k f t (δ mgt ) is the angle of the voltage at node k (m) in phase f (g) at time period t; θ k f mg represents the angle of the admittance that relates node k at phase f with node m at phase g; and ∆ t is the time period during which the power demands remain constant. It should be noted that F , N , and T are the sets that contain all the phases, nodes, and time periods, respectively.

Remark 1.
Notably, the objective function is nonlinear and non-convex owing to the product among voltages and trigonometric functions; this implies that its minimization requires advanced optimization strategies that can efficiently deal with nonlinearities [30]. This justifies the use of the master-slave optimization strategy, as in the case of the improved version of the CBGA proposed herein, owing to its simplicity in terms of programming and its high efficiency in relation to the processing times.

Set of Constraints
The set of constraints for the phase-balancing problem is typically associated with the operative conditions of the network, that is, the power balance equations and the voltage regulation bounds [2]. These equations are as follows: where P s k f t represents the active power generated by source s connected at node k in phase f at time period t; Q s k f t represents the reactive power generated by source s connected at node k in phase f at time period t; P d kgt represents the active power demand at node k in phase g at time period t; Q d kgt is the reactive power demand at node k in phase g at time period t; x k f g is a binary variable that defines the connection of the demand at node k at f in phase g; and V min and V max correspond to the minimum and maximum voltage regulation bounds allowed for all nodes of the network at each period of time, respectively. Remark 2. Expressions (2) and (3), associated with the active and reactive power balance equations at each node, phase, and period of time, evidence the complication of the three-phase power flow problem in electrical networks, even when the connection of the loads at these nodes are well-known. This is due to the nonlinear, non-affine equations that must be solved numerically.

Mathematical Model Interpretation
The following interpretations can be drawn from the optimization model described in (1) to (6). Equation (1) corresponds to the objective function of the phase-balancing problem, which is formulated as the minimization of the daily costs of energy loss in all branches of the network. Equations (2) and (3) are the active and reactive power balance constraints that must hold at every phase, node, and time period in the threephase distribution network. These constraints originate from the three-phase power flow problem in electrical grids, which has been solved numerically using numerical methods such as the Newton-Raphson [31], backward/forward [32], and graph-based methods [33]. Equations (4) and (5) guarantee that the loads are connected to the phases in a unique form by using a matrix connection at each node (i.e., bus k), formed by the variables x k f g , in which each row and each column must be equal to 1. Finally, the box-type constraint (6) defines the upper and lower voltage regulation bounds of the network, which are typically ±10% for medium-voltage levels [2].

Solution Methodology
The problem of optimal phase balancing in electrical distribution systems can be solved using a master-slave metaheuristic optimization approach, which represents binary variables of the optimization model (1)-(6) with an integer codification that represents the possible phase connections per node. The phase connection is defined in the master stage, which corresponds to an improved CBGA; in the slave stage, each phase configuration is evaluated to determine the total costs of daily energy losses by using the successive approximation power flow method. Details of the proposed master-slave optimization approach are presented in the following subsections.

Slave Stage: Three-Phase Successive Approximation Power Flow Method
The successive approximation power flow method is a recently developed power flow method for single-phase and direct current distribution grids [28]; however, thus far, its extension to three-phase unbalanced distribution grids with Y and ∆ loads has not been reported. The main characteristics of the successive approximation power flow method are that its recursive formula does not use derivatives and the matrices that intervene in the iteration process are constant, implying that the required processing times to solve the power flow problem are in the order of milliseconds [2]; this power flow method can be formulated directly in a complex domain, which simplifies its computational implementation.
To formulate the three-phase successive approximation power flow method, we present the relation between the nodal voltages and the injected currents using the admittance matrix, as follows: where V s3ϕ,t is the vector that contains all the complex voltages in the slack sources (note that these voltages are clearly known for power flow purposes); V d3ϕ,t is the vector that contains all voltages in the demand nodes (unknown variables of interest); I s3ϕ,t is the vector of the net injected current in the slack nodes; I d3ϕ,t represents the vector that contains all currents in the demand nodes; Y gs3ϕ is a component of the admittance matrix that relates the slack sources among them; Y gd3ϕ = Y T ds3ϕ is the component of the admittance matrix that relates the slack and demand nodes and vice versa; and Y dd3ϕ is the component of the admittance matrix that relates the demand nodes among them. It should be noted that the voltages and currents in Equation (6) are defined under the three-phase condition and ordered by the nodes and phases, respectively; this equation is applied for each time period.
From (6), it is evident that the second row contains the demand nodal voltages (i.e., V s3ϕ,t ), which are the variables of interest in power flow studies; for simplicity, this equation can be rewritten as follows: Equation (8) enables the determination of all the voltage profiles at all the demand nodes; however, it is necessary to determine the form of calculating the demand current because this depends on the voltage profile and the amount of power demanded, including whether the loads have Y− and ∆−connections.
The case where node k has a load with a Y-connection (we assume that this is solidly grounded) is depicted in Figure 1. In the three-phase load with a Y-connection, we assume that the voltage subjected to each load is phase-to-neutral, and each line current is calculated as follows: which can be compacted as presented below: To obtain the line currents for a three-phase network with loads connected in ∆, the load diagram presented in Figure 2 is considered for the k node.
, ∀t ∈ T which can be compacted as follows: where the matrices M and H take the following form: Considering a case where a three-phase distribution grid has loads connected with the Y and ∆ structures, for calculating the demanded currents, these connections must be considered, as presented in Equations (9) and (10), respectively. When solving the power flow problem, the main goal is to determine the daily energy losses for each configuration of the phases provided in the master stage. To this end, we employ the admittance nodal matrix, as follows: where S loss represents the apparent power losses of the network. Algorithm 1 shows the general implementation of the successive approximation power flow for unbalanced distribution networks with loads connected in Y and ∆.
Note that the variable S loss,t contains information regarding the energy losses of the system for the phase configuration provided by the master stage at each time period. To calculate the objective function defined in Equation (1), we use the value of the daily energy losses provided by the three-phase successive approximation power flow in Algorithm 1, as follows: Remark 3. Note that expression (12) presents the connection of the master and slave stages because the slave stage is necessary to determine the objective function value of each phase configuration for the loads provided by the proposed improved CBGA.
Algorithm 1: Solution of the three-phase power flow problem for unbalanced distributions 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 admittance matrix Y bus 3ϕ ; Extract the components Y ds3ϕ and Y dd3ϕ ; Compute and store the three-phase impedance-like matrix Z dd3ϕ = Y −1 dd3ϕ ; Define the maximum number of iterations m max ; Define the convergence error ; Define the number of periods of analysis t max ; Define the substation voltages: Compute the demanded current I m dk3ϕ,t using Equation (9) ; else Compute the demanded current I m dk3ϕ,t using Equation (10) ; Calculate the new nodal voltages V t+1 d3ϕ using Equation (7), as follows Calculate the energy losses using Equation (11); break; else Make V m d3ϕ,t = V m+1 d3ϕ,t ;

Slave Stage: Improved CBGA
The master stage is entrusted with providing the phase configuration to the successive approximation power flow approach defined in Algorithm 1. Here, we propose an improved version of the classical CBGA in the selection and mutation stages by introducing a normal Gaussian distribution, which is used in the vortex search algorithm [34]. To explain the improvements in the CBGA, we begin with the codification of the phase-balancing problem, which is based on the possible connections reported in Table 1. Table 1. Possible connected loads in a three-phase distribution grid [13].

Connection Type
Phases Sequence The proposed codification takes the following structure: where x s i is the individual i in iteration s, which has a dimension of 1 × (n − 1). The main advantage of this codification is that the 9 variables associated with the binary variables x k f g per node are represented solely by an integer between 1 and 6. The improved version of the CBGA is developed using probability criteria at each iteration, which defines the methodology used to generate the offspring. If the probability criterion is lower than 50% in iteration s, the classical CBGA with multi-point mutation is implemented; otherwise, the vortex search algorithm is used to generate the offspring.

Classical Approach
The classical approach is based on the selection of two individuals from the population, that is, x s i and x s j being i = j from the current population. These two solutions are recombined using an arbitrary point randomly generated as a number between 1 and n − 2.
To illustrate this stage, suppose that a test feeder is composed of 11 nodes (1 slack and 10 demands), in which the individuals x s i and x s j take the following form: For each descending individual, the multi-point mutation criteria are applied as follows: (i) the number of mutation points is randomly selected from 1 to 1 + 20%(n − 1)rand(1) of the population size (i.e., 1 and 3). Let us suppose that, for individual y s i , the number of mutation points is 2, whereas for individual y s j , the number of mutation points is 3. Thus, the selection of mutation points is randomly generated for each descending individual; here, we consider that these points are 1 and 8 for y s i and 2, 5, and 7 for y s j . A random number between 1 and 6 is generated at each of these points to define the phase connection at this node. This procedure yields the following structure in the current offspring: The remaining steps involve the evaluation of the objective function (12) to select which among the individuals y s i and y s j has the opportunity to replace the current population, considering the case where the objective function is better than the worst individual in the population and is different from all the other ones (diversity criteria).

Improved Approach: Vortex Search for Offspring Generation
In this stage, the main goal is to add to the classical CBGA the possibility of generating an offspring population based on the evolution criteria of the VSA [2]. The main characteristic of the VSA is the usage of a normal Gaussian distribution to generate a neighborhood around the current solution, called the center of the hyper-ellipse µ [35]. Here, we select the center of the hyper-ellipse as a random individual in the current population during iteration s as µ S = x s i . This center is generated using the following distribution probability: where d = (n − 1) represents the dimension, y is a random vector featuring values between zero and one with dimension d × 1, and Σ is the covariance matrix. If the values of the diagonal elements of Σ are equal and those outside the diagonal are zero, then the resulting shape of the distribution will be a hyper-ellipse; therefore, equal variances with zero covariance can be calculated as the value of Σ, as presented in Equation (15).
where σ represents the variance of the distribution, and I represents the identity matrix of dimension d × d. Equation (16) can be used to calculate the current standard departure σ s of the distribution for iteration s.
where σ s is also considered as the radius r s of the current hyper-ellipse. As a first step, by multiplying the radius, the region of possible solutions for the offspring is generated. The set of neighboring solutions obtained and contained in C s i (y), which are verified using the algorithm and by considering the lower and upper limits set in µ s , is defined by Equation (17).
where rand is a function that generates a random number between 1 and 0 with a normal distribution. After verifying the limits, at this point, the exploiter stage is initiated. It selects the optimal response within the solution space of C s i obtained in (17), which is selected as the potential individual to be part of the current population of the CBGA, provided it is better than the worst solution in the current population and fulfills the diversity criteria.
It should be noted that the radius of the hyper-ellipse is gradually reduced as the optimization process progresses. Here, we selected a linear rule to decrease the radius as follows: where s max is the maximum number of iterations of the CBGA.

Remark 4.
The size of the offspring in C s i (y) is selected as 20% of the population size in order to reduce the number of evaluations required in the slave stage for determining the costs of daily energy losses.

General Flow Diagram of the Proposed Master-Slave Improved CBGA
The implementation of the proposed improved CBGA considering multi-point mutation and offspring generation with the VSA is presented in Algorithm 2. else Select an arbitrary individual x s i from the population and make µ s = x s i ; Define the current radius of the hyper-ellipse r s ; Define the number of individuals that will generate 20% of the population size p s ; Create the neighborhood C s i (y) using (14); Revise the upper and lower bounds for each individual i in C s i (y) using (17); Evaluate the objective function for each individual i in C s i (y) using Algorithm 1; Define the winning individual in C s i (y) as the individual with the lower objective function value and term it as C s best ; if C s best ≤ f 1 (x s worst ) then Verify the diversity criteria and include the winning individual C s best in the current population by replacing x s worst with C s best ; Order the population in the ascending form of the objective function value;

Three-Phase Test Feeders
Two test feeders are considered to validate the proposed optimization approach based on the improved version of the CBGA with the multi-point mutation strategy and the VSA for generating offspring. The main characteristics of the test feeders are discussed below.

Computational Validation
This section presents the numerical validation of the proposed approach in the 15and IEEE 37-bus test feeders, considering the following scenarios: (i) application of the improved CBGA and comparative methodologies to the 15-bus test system under the peak load condition, and (ii) usage of active and reactive demand load curves in the IEEE 37-bus test feeder to minimize the annual operative costs of the network.

Analysis of the Peak Load Condition
In this simulation case, we present the positive effects of the power loss reduction after applying a phase-balancing approach in the 15-bus test feeder, considering the peak load scenario. Table 6 presents the optimal solutions obtained via the classical CBGA [27], the sine-cosine algorithm (SCA) [36], the black hole optimizer (BHO) [37], the VSA [2], and the improved CBGA. Notably, the proposed and comparative algorithms were set with 10 individuals in the population, 1000 iterations, and 100 consecutive evaluations to determine the average processing time. The results in Table 6 show that (i) all the methodologies enable a reduction of more than 18% in the amount of power loss in the peak load case. The worst result was yielded by the BHO, with a reduction of 18.06%, whereas the best result was achieved by the improved CBGA, with a reduction of 18.66%. All the compared methods, including the proposed approach, require less than 40 s to identify an adequate solution by using the phase-balancing configuration. Nevertheless, the proposed approach is the fastest, with an average processing time of 2.0760 s. (iii) The classical CBGA and the VSA afford the second and third best results, respectively. Thus, the combination of the best characteristics of both these methods in the improved CBGA enables excellent numerical results for the phase-balancing problem.
Note the following in Table 6: our proposal-the proposed improved version of the CBGA-is at least three times faster than the classical CBGA, which is the second faster method; in addition, the VSA method proposed in [2] is more than nineteen times slower than our proposed optimization approach. Those comparisons confirms the effectiveness of the approach here presented to solve complex MINLP problems with solution spaces with large dimensions.
To verify the effectiveness of the redistribution of all the loads in the distribution network, Figure 5 presents the amount of power loss per phase for the benchmark case, as well as for the solution obtained using the improved CBGA.
The variations in the power losses per phase, as shown in Figure 5, indicate that, after applying phase balancing with the improved version of the CBGA, the losses of phases b and c increase by 9.1521 kW and 21.8832 kW, respectively. By contrast, the power loss of phase c decreases by approximately 56.0845 kW, which clearly compensates for the increments experienced in phases a and b. In addition, the power losses per phase are balanced at approximately 36 kW, with small differences lower than 2 kW; this shows that the phase-balancing approach using the improved CBGA effectively redistributes the loads in all the phases as uniformly as possible.  In the case where the 15-bus system has all the loads connected in ∆, the power loss in the benchmark case is approximately 110.7841 kW, which can be reduced to 107.0944 kW by using the proposed improved version of the CBGA. The codification that presents this solution is {5, 6, 6, 2, 3, 2, 3, 5, 5, 2, 4, 3, 5, 6}. The difference between the loads connected in Y and ∆ is notable in terms of the amount of power loss in the benchmark case, because the difference between both connections is approximately 23.4631 kW. In addition, the percentage of power loss reduction in the case of ∆−connections is just 3.33%, which is constrained by the Y−connection case, where this reduction exceeds 18%. The average processing time in the case of ∆−connections is 1.6199 s, which is 0.4563 s less than that for the Y−connections.

Minimization of Annual Energy Loss Costs
To demonstrate the effectiveness of the improved version of the CBGA for performing the phase-balancing task in three-phase distribution networks, the IEEE-37 node test feeder was evaluated under a daily power flow environment, considering active and reactive power demand curves. The demand curves are presented in Figure 6 and also reported in Table 7, for a comparison of our results with future research (note that the scaling factor of the data reported in Table 7 for active and reactive power demands is 2).  To evaluate the objective function in terms of the annual operative costs, the cost of the energy is assumed in US$/kWh 0.1390 (this value is the average cost of the energy in Bogotá, Colombia in May 2019 [38]); the number of days is considered as 365, and the length of the power flow period, ∆ t , is 0.5 h. Table 8 presents the numerical performance of the improved version of the CBGA when 10 individuals, 1000 iterations, and 100 consecutive evaluations are adopted for the IEEE 37-bus system. Based on the results reported in Table 8, it is evident that (i) the best numerical solution (see solution 1) enables a reduction in the annual operative costs by approximately US$/year 8121.7220, which corresponds to a reduction of 18.79% with respect to the benchmark case; and (ii) the difference between solution 1 and 10 is just US$/year 75.1586, which implies that the 10 solutions reported in the final population of the improved CBGA are adequate to solve the phase-balancing problem of three-phase networks. The implementation of any one of these depends exclusively on the practices of the utility.
To demonstrate the effectiveness of the proposed optimization approach, Figure 7 presents the percentages of annual reduction in the energy costs for each of the ten solutions reported in Table 8. Note that the differences among all solutions are lower than 0.18%, thereby confirming that the proposed approach can efficiently solve the phase-balancing problem of unbalanced electrical three-phase grids. In addition, it is important to mention that the average processing time for solving the phase-balancing problem in the IEEE 37-bus system is 136.3901 s. This is a short time, considering that, for each combination of the phases provided by the CBGA, 48 power flows are evaluated, which implies that, for a population of 10 individuals, 1000 iterations with two offspring individuals per iteration are evaluated as approximately 96,480 power flows.  To elucidate the effect of phase balancing on the amount of daily energy loss, shown in Figure 8, we present the daily energy losses per phase before and after applying the phase-balancing strategy (solution 1 in Table 8  The results presented in Figure 8 indicate the following: (i) the benchmark case shows that the difference in the daily energy losses between phases b and c is approximately 562.2398 kWh/day, which confirms the higher imbalance in the IEEE 37-bus system; (ii) after applying the phase-balancing strategy, the difference between phases b and c is just 127.1357 kWh/day, which represents a reduction of 77.3876% with respect to the benchmark case; and (iii) the phase-balancing approach leads to better load distribution per phase, which balances the power losses among them with differences lower than 130 kWh/day in the worst case.

Conclusions
In this research, the problem of phase balancing in three-phase distribution grids was addressed from the master-slave optimization perspective. To this end, an improved version of the CBGA using the VSA was proposed. The master stage provides the configuration of the load phases through the CBGA by using an integer codification between 1 and 6, which represents all the possible load connections in a three-phase grid. Alternatively, the slave stage evaluates each of the load configurations in the three-phase power flow tool, corresponding to an extended version of the successive approximation power flow method for three-phase unbalanced networks.
Numerical results for the 15-bus test system demonstrate that the improved version of the CBGA proposed herein attains the highest quality solution with a reduction of 18.66%, as compared with the classical CBGA (18.64%), BHO (18.06%), SCA (18.51%), and VSA (18.57%). In addition, the proposed approach exhibited the fastest processing time (2.0762 s) for solving the phase-balancing problem in the 15-bus system, followed by the classical CBGA (6.9435 s). These processing times are important because the dimension space of the phase-balancing problem potentially increases, that is, 6 (n−1) , implying that the 15-bus system corresponds to 6 14 = 78,364,164,096, that is, more than 78,000 million combinations.
In the IEEE 37-bus system, the improved version of the CBGA reduces the annual energy loss cost by approximately 18.79%, with an average processing time of 136.3901 s. In addition, the difference between the ten best solutions reported by the proposed CBGA is less than 0.18%, which is less than 76 dollars per year of operation. This confirms the effectiveness of the proposed approach in solving complex MINLP models, such as the phase-balancing problem in three-phase networks.
In future works, it will be possible to (i) combine the proposed three-phase-balancing approach with the optimal location of capacitor banks to reduce the annual operating costs of three-phase grids; (ii) study the optimal integration of static distribution compensators in three-phase networks by using a hybrid master-slave approach employing the successive approximation power flow method in the slave stage; and (iii) use the improved version of the CBGA to solve the problem of the optimal selection of conductors in three-phase networks.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.
Acknowledgments: This work was supported in part by the Centro de Investigación y Desarrollo Científico de la Universidad Distrital Francisco José de Caldas under grant 1643-12-2020 associated with the project: "Desarrollo de una metodología de optimización para la gestión óptima de recursos energéticos distribuidos en redes de distribución de energía eléctrica." and in part by the Dirección de Investigaciones de la Universidad Tecnológica de Bolívar under grant PS2020002 associated with the project: "Ubicación óptima de bancos de capacitores de paso fijo en redes eléctricas de distribución para reducción de costos y pérdidas de energía: Aplicación de métodos exactos y metaheurísticos."

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
Here is presented a list with the main abbreviations used along within this document. Angle of the voltage at node k in phase f at time period t (rad). δ mgt Angle of the voltage at node m in phase g at time period t (rad). I d3ϕ,t Complex three-phase current at demand node d in time period t (A). I dk3ϕ,t Complex demanded three-phase current at k in time period t (A). I dka,t Complex demanded current at node k in phase a at time period t (A). I dkb,t Complex demanded current at node k in phase b at time period t (A). I dkc,t Complex demanded current at node k in phase c at time period t (A). I s3ϕ,t Complex three-phase current at source node s in time period t (A). S loss Apparent power losses of the network (VA) S dk3ϕ,t Complex demanded three-phase power at k in time period t (VA). S dka,t Complex demanded power at node k in phase a at time period t (VA). S dkab,t Complex demanded power at node k between phases a and b at time period t (VA). S dkb,t Complex demanded power at node k in phase b at time period t (VA). S dkbc,t Complex demanded power at node k between phases b and c at time period t (VA). S dkc,t Complex demanded power at node k in phase c at time period t (VA). S dkca,t Complex demanded power at node k between phases c and a at time period t (VA). V d3ϕ,t Complex three-phase voltage at demand node d in time period t (V). V dk3ϕ,t Complex demanded three-phase voltage at k in time period t (V). V dka,t Complex demanded voltage at node k in phase a at time period t (V). V dkb,t Complex demanded voltage at node k in phase b at time period t (V). V dkc,t Complex demanded voltage at node k in phase c at time period t (V). V s3ϕ,t Complex three-phase voltage at source node s in time period t (V). Y dd3ϕ Three-phase submatrix of the admittance nodal matrix that relates demand nodes among them (S). Y ds3ϕ Three-phase submatrix of the admittance nodal matrix that relates demand and source nodes among them (S). Y sd3ϕ Three-phase submatrix of the admittance nodal matrix that relates source and demand nodes among them (S). Y ss3ϕ Three-phase submatrix of the admittance nodal matrix that relates source nodes among them (S). θ k f mg Angle of the admittance that relates node k at phase f with node m at phase g (rad). C kWh Average cost of energy (US$/kWh). f 1 Cost of daily energy losses (US$/day). P s k f t Active power generated by source s connected at node k in phase f at time period t (W). P d kgt Active power demanded at node k in phase g at time period t (W). Q s k f t Reactive power generated by source s connected at node k in phase f at time period t (var). Q d kgt Reactive power demanded at node k in phase g at time period t (var). V max Maximum voltage regulation bound (V). V min Minimum voltage regulation bound (V). V k f t Voltage magnitude at node k in phase f at time period t (V).

V mgt
Voltage magnitude at node m in phase g at time period t (V). x k f g Binary variable that defines the connection of the demand in node k at f in phase g. Y k f mg Magnitude of admittance that relates node k at phase f with node m at phase g (S).