Application of the Vortex Search Algorithm to the Phase-Balancing Problem in Distribution Systems

: This article discusses the problem of minimizing power loss in unbalanced distribution systems through phase-balancing. This problem is represented by a mixed-integer nonlinear-programming mathematical model, which is solved by applying a discretely encoded Vortex Search Algorithm (DVSA). The numerical results of simulations performed in IEEE 8-, 25-, and 37-node test systems demonstrate the applicability of the proposed methodology when compared with the classical Cuh & Beasley genetic algorithm. In addition, the computation times required by the algorithm to ﬁnd the optimal solution are in the order of seconds, which makes the proposed DVSA a robust, reliable, and efﬁcient tool. All computational implementations have been developed in the MATLAB ® programming environment, and all the results have been evaluated in DigSILENT© software to verify the effectiveness and the proposed three-phase unbalanced power-ﬂow method.


General Context
Energy has become a fundamental base for society, covering the basic needs of people, contributing to technological development and making a change to how people work and live in their houses [1]. Initially, it was necessary to look for a method to deliver energy from generation to consumption locations (end users) in an efficient form using equipment and devices to generate, transport, distribute, and commercialize energy, all of these encompassing the well-known electric power system [1]. One of the fundamental components for power systems is distribution systems, which are defined as, according to [2], the set of elements that transport electric energy from a power substation to an end user. Generally, distribution systems are formed of a power substation, subtransmission substation, distribution substation, primary feeders, distribution transformers, and secondary and auxiliary transformers [3]. The distribution system is entrusted to supplying safe and reliable high-quality energy economically to end users [1].
With technological development and population growth, unbalanced power systems have taken up a prominent role within studies of power distribution systems, given that these generally have unbalanced operation due to the load nature and network configuration, with one-phase distribution transformers connected arbitrarily at the power-system phases [4]. This is consistent since it is more economic to implement branches from one or two phases of the trunk circuit than installing three-phase circuits directly to the end

Motivation
When power loss is reduced, power-supply service provided by utilities will be more efficient [11,12], complying with energy, power and service-quality parameters, meeting end user needs, according to guidelines of regulatory entities [1]. For the Colombian case, the Energy and Gas Regulation Commission (CREG) states that voltage profiles should be within ±10% of nominal operation voltage [13] and frequency within ±3% of nominal operation frequency of the system [14]. If these variables are not in compliance with operation limits, network operators incur financial penalties. In this manner, distribution companies will benefit from system-balancing, with a consequent reduction to operation and investment cost, since less energy has to be purchased to meet user demand and there will be no investment to implement other higher-cost solutions, such as capacitors and DG installation, which can be more expensive than phase-balancing [15]. For this reason, the present research shows the minimization of power loss caused by the unbalance of three-phase power distribution systems. For this purpose, the optimization algorithm known as Discrete Vortex Search Algorithm (DVSA) is used, combined with a three-phase version of the backward-forward sweep iterative method, with which the set of optimal connections for the loads that are part of the system will be determined, with an optimal point of operation, decreasing the phase-unbalancing.

Review of the State of the Art
In the specialized literature, several methods and proposals have been proposed for phase-balancing (PB). This problem was introduced for the first time in 1998 by [16], who formulated PB via mixed-integer programming with linear constraints, seeking to reduce current unbalancing in distribution system lines. Likewise, power loss is reduced, and voltage profiles and power quality are improved, in such a way that the service provided to end users is efficient and of good quality. In [17] PB was developed through a metaheuristic algorithm known as simulated annealing, which minimized the total cost of operation. Voltage drops on the network segments were reduced, therefore reducing energy loss in the system, and the power-supply service to end users less expensive, since less energy was required to meet demand.
In [18], PB is proposed using a Genetic Algorithm (GA), designed to minimize phaseunbalancing of the system and reduce voltage drops on the network segments, power loss, and current through the neutral wire. It is determined that with correct phase reorganization, service quality is significantly improved. In [11], PB was developed via GA, in which power loss is reduced in distribution systems, with an improvement to the line transportation capacity, and decrease to both the voltage drop on the network segments and energy consumption price, in such a way that network operators and end users can be benefit from an economic and reliable power-supply service.
The authors in [19] formulated PB via the backward recurrence search heuristic algorithm, which minimized phase-unbalancing of the system such that neutral current is reduced. It was concluded that a more reliable system was obtained, since the likelihood of overcurrent protection devices actuation was decreased in the neutral. In [20], PB was developed through ant colony optimization, where energy loss was minimized over a 24 h interval for distribution systems, which reduced the neutral current and energy consumption price on the side of the end user. With this methodology, it was concluded that the proposed algorithm is beneficial for power systems with limited economic resources, since the configurations for phase reorganization are at the optimal during the period under study.
In [21], a particle swarm optimization PSO algorithm for the PB was proposed, in which power loss is decreased, and voltage profiles are improved. The right phase reorganization significantly improves the quality of the power-supply service to end users. In [22] PB was formulated through an expert system based on heuristic rules, which had the purpose of minimizing operative costs of the system. The objective function was formulated considering the cost of power loss and cost of service interruption due to the work for phase-swapping, which reduces both system-unbalancing and current flowing through the neutral and energy losses over a 24 h period. The authors concluded that the current of the neutral wire was significantly reduced, enhancing power-system reliability; moreover, transportation capacity of the network segments was improved, decreasing power loss.
In [23], PB was developed via an immune algorithm where the operative costs of the system were reduced over a 24 h interval. The objective function was formulated considering current flowing through the neutral (adding a penalization cost if this exceeds a predetermined value), service interruption cost, and the cost due to phase reorganization, in such a way that service quality provided to end users can be improved, since the reduction of system-unbalancing means a reduction of power loss. In [24], the problem of PB is addressed with automatic phase-swapping via information system-mapping, taking into account the technical characteristics and the load flow, to decide the best phase reorganization that provides the lowest power loss. In this sense, the algorithm developed is a useful tool for power distribution engineers in the direction and maintenance of systems, since the automatic feature allows operation of the system in an efficient and effective manner.
The authors of [25,26] formulated PB using a fuzzy geek hybrid heuristic algorithm and an auto-adaptive hybrid differential evolution algorithm to minimize the system phase-unbalancing based on power-flow equations and capacity and voltage constraints, seeking to balance the current flowing through the system phases and decreasing the neutral current. With the correct phase reorganization, the service quality is considerably improved, and the energy invoice is decreased for end users, since when system is balanced, power loss is also reduced.
In [27], a bacterial foraging optimization algorithm oriented by a PSO algorithm was proposed, which had the objective of minimizing the operative costs of the system, whether it was radial or meshed. The objective function was formulated considering the current flowing through the neutral wire, voltage drop on the network segments, system power losses (adding a penalization cost if this exceeded a set value), and the cost due to phaseswapping, in such a way that the system balance is improved, making as similar as possible both the voltage between phases and the current flowing through the lines, in this way reducing the neutral current and the system power loss.
In [3], PB through a GA was formulated with the objective of minimizing energy loss in distribution systems, and consequently, voltage profiles were improved and the congestion of the system lines was reduced, which had direct effects on investment costs, with optimal system operation points and reduced unbalance factor, power loss, and neutral current. In [28], the PB problem was addressed with a heuristic method based on phase-swapping for load-balancing, minimizing the degree of system-unbalancing. With this, voltage profiles were improved, and energy losses decreased over 7-day period for the implemented test systems. In this manner, the systems were taken to an optimal operation point.
In [8], PB was formulated in the 13-node IEEE test system with a balancing algorithm based on verifying at each candidate node the possible combinations for phase-swapping until the option that balanced the system was obtained. In this manner, a configuration was obtained that reduced the neutral current and power loss, and improved the voltage profiles in the test system. In [29], PB was developed via PSO, which sought to minimize the phase-unbalancing of the system. Likewise, power loss was reduced, and voltage profile and energy efficiency were improved for a good-quality and economic power-supply service to end users.
To address the PB problem in [30], the minimization of the unbalancing in distribution systems was proposed, using a search heuristic algorithm that defined how to perform phase-swapping with load contactors. This method decreased power loss and provided better results depending on the number of contactors used; furthermore, this increased the operation cost of the system. The authors of [31] proposed PB with a GA that used a codification based on group theory that minimized power losses in distribution systems and reduced the existing unbalancing between the phases of the system.
In [32], PB was developed through an analytic approach, taking into consideration that phase-swapping must be minimal. (This was done via the implementation of a search vector that selected the priority nodes considering the load connected at each phase and voltage unbalancing), to reduce the distribution system-unbalancing and therefore the power losses. Finally, in [33], PB implementing a GA was developed which had the objective of reducing energy loss over a 24 h interval and decreased the phase-unbalancing in the system.
As observed in the above review of the state of the art, PB, regardless of the objective and the solution techniques mentioned, seeks to improve the quality of power and energy so that the service provided by network operators to end users is efficient and economic, resulting in a beneficial situation for customers and utilities. Tables 1 and 2 show a summary of the works consulted for the state of the art, with the solution method and objective function employed.
For the development of this research paper, an optimization metaheuristic technique will be used called Vortex Search Algorithm (VSA), as it has not been applied for the PB problem, as observed in Table 1. This technique is used for the solution of problems in the continuous domain [34]. In this case, a discrete version of the VSA, known as DVSA, is proposed, since the optimal combination for load connection must be selected. The objective function to minimize is the power loss in the distribution system, because this has been the least researched variable within the references consulted, as shown in Table 2, to the extent that only the 14% of the authors had it as the main purpose of their research.

Contributions
This study addresses the solution of the optimal PB problem in unbalanced distribution networks. After doing a review of the state of the art, a master-slave optimization is proposed as follows: in the master stage of the DVSA, the set of connections for loads presented in the system is defined; in the slave stage, the solution for the power flow is established via the three-phase sweep iterative method, which will determine the power loss for the set of solutions defined in the master stage. Thus far, this solution method for PB has not been reported in the specialized literature. The contributions of this research work are presented as follows: It is possible to obtain a global optimal solution for the PB problem via the parameterization solution method with a sensitivity analysis between the number of candidate solutions and the number of iterations. The solution developed from the interaction of the DVSA and the three-phase iterative sweep method can be implemented for any radial or meshed unbalanced distribution system. Likewise, the set of optimal connections for loads in the system is not limited by the number of nodes.
Finally, it is important to mention that, while the PB problem has been widely studied and reported in the specialized literature, as shown in the state-of-the-art review, a discrete codification for the VSA that uses integer decision variables will simplify the binary vector dimensions commonly used for the solution of this problem, decreasing the processing time [35]. In addition, to present the effectiveness and robustness of the proposed discrete VSA approach, numerical comparisons are included with the classical Chu & Beasley genetic algorithm to determine the final phase configuration and the total amount of power-loss reduction.

Document Organization
The rest of this document is organized as follows: Section 2 presents the mathematical formulation of the optimal phase-balancing problem in unbalanced distribution systems. Section 3 shows the proposed master-slave methodology, where the DVSA is presented integrated with the three-phase iterative sweep power flow. Section 4 describes the test systems used in this study with their main characteristics. Section 5 shows the results of the operative state of the system (nodal voltage per phase) and the active power loss for a given demand condition, using the three-phase iterative sweep power flow. Section 6 shows the results obtained for the set of connections and the power losses in each of the cases proposed by the master-slave methodology. Section 7 gives conclusions and proposes future works.

Mathematical Formulation
The optimal phase-balancing problem in unbalanced distribution systems can be represented via a mixed-integer non-linear programming (MINLP) model. The integer variables (binary type) represent the decision variables, which are associated with the selection of connection type for each load presented in the system. On the other hand, the continuous variables are associated with the classic power-flow formulation, which uses magnitude and angle voltage per node [36]. Finally, the non-linear nature is due to the trigonometric functions in the power-balance equations and the products between the magnitudes of the nodes [35].
The MINLP formulation is presented below, as an adaptation of the models in [3,37].

Objective Function
The objective function corresponds to total active power loss of the unbalanced distribution system, as presented in (1).
where z represents the value of the objective function. Ω P is the set associated with the system phases. Ω N is the set associated with the system nodes.

Constraints
The set of constraints corresponds to the operative limitations that can be found in any unbalanced distribution system. These are shown in (2) to (11).
where Expressions (2)-(4) represent the active power balance per phase at each node of the system. P gA i , P gB i and P gC i are the active power due to a conventional generator in the phases A, B, and C, respectively, at the node i. P dA ih , P dB ih and P dC ih are the active powers demanded, as per the type of connection h, in the phases A, B, and C, respectively, in the node i. X ih represents the binary decision variable, which takes the value of 1 if connection h is selected for a node i, and 0 otherwise.
The binary variable X ih will have to select among h possible connections for the loads of the nodes of demand i presented in the distribution system, as shown in Table 3. It must be considered that in the case of having inductive loads (motors), the phase sequence should not be swapped, to prevent damage to the equipment [3]. An example of this situation is that if there is a three-phase motor connected in phases A, B, and C, this can only be connected as follows: BCA or CAB, making sure it has the same sequence. Table 3. Possible types of connections for the loads presented in the system [3].

Type of Connection (h)
Connection Sequence As a result, it is proposed that each load connected to the nodes of demand i is associated with a type of connection, which will represent a possible phase-swapping [3]. The type of connections vary from 1 to 6. Type 1 indicates that there was no phase-swapping (default setting); the remaining five types of connection indicate a phase-swapping that may or may not involve a sequence change [3], i.e., if there is a load initially connected in ABC sequence and a Type 6 is required, phase connection will change to BAC [3], as shown in Figure 1. This research work is focused on three-phase loads; however, if there are two-phase or one-phase loads, the type of connections of Table 3 are applied accordingly. In the first case, it is assumed that the power of the phase where the load is not connected is equal to zero, and in the second case it is assumed that the power of the two phases where the load is not connected is equal to zero [3].
{∀i ∈ Ω N } where Expressions (5)-(7) represent the reactive power balance per phase at each node of the system. Q gA i , Q gB i , and Q gC i are the reactive power generated by a conventional generator, in the phases A, B, and C, respectively, in the node i. Q dA ih , Q dB ih , and Q dC ih , are the reactive powers demanded, as per the type of connection h, in the phases A, B, and C, respectively, in the node i.
where Expressions (8)- (10) correspond to the constraints of voltage regulation per phase in the node i. Voltage per phase in the node i can take any value as long as this is within V min and V max .
where the expression (11) is limited to one, the number of possible connections presented in each demand node i. The model described in (1) to (11) represents the optimal phase-balancing problem in unbalanced distribution systems. This model is based on minimizing the power losses in unbalanced distribution systems, where the binary variable X ih will determine the set of connections for the demand nodes that allows the system to find an operation point to reduce phase-unbalancing. The solution for the model described in (1) to (11) requires numeric optimization techniques, since the active and reactive power-balance equations, defined in (2) to (7), are non-linear and non-convex [34]. For this reason, in Section 3 a master-slave methodology is proposed to solve the optimal phase-balancing problem. This calls for combining the DVSA and the three-phase version of the power flow, which is known as the iterative sweep method.

Methodology Proposed
For the solution of the optimal phase-balancing problem in unbalanced distribution systems, to minimize power loss for a given demand condition, in this document the VSA is proposed in its discrete version [35] to be used in the master stage, and in the three-phase iterative sweep power-flow method in the slave stage. In the master stage, the phase connections will be defined in each demand node of the system, which allows balancing, while in the slave stage the constraints will be assessed for the power flow defined in (2) to (7). Each component of the master-slave methodology proposed is described as follows.

Master Stage: Vortex Search Algorithm
VSA is a metaheuristic technique used to solve problems in the continuous domain, and it is based on the behavior of vortex created in agitated fluids [38]. A significant advantage of this optimization method is that candidate solutions are generated around the best current solution via a Gaussian distribution at each iteration [38] and have variable radii that allows the exploration and exploitation of the solution space [34]. When this method is applied to power-flow problems, it is warrantied to reach an optimal global solution with the minimum of standard deviation [39].
As mentioned before, the VSA works with continuous variables. In this paper, a discrete version is employed when these variables are rounded to the closest integer, as proposed in [35], which will allow the slave stage to evaluate a set of connections for the loads that are part of the system, exploring and exploiting the solution space in a proper manner.

Initial Solution
VSA works with non-concentric hyperspheres, where the radius of the external hypersphere represents the size of the solution space, and the center represents the best current solution [34]. The initial center of the hypersphere is defined in (12): where x max and x min are d × 1 vectors that define the upper and lower limits of the decision variables of the optimization problem in a dimensional space d.

Candidate Solutions
The set of candidate solutions is represented as C t (s), where t is the number of iterations. Initially t = 0 and the set of solutions C 0 (s) is generated through a random process using a Gaussian distribution in the dimensional space d, where C 0 (s) = {s 1 , s 2 , . . . , s n }, n being the number of candidate solutions. To generate the Gaussian distribution in a multidimensional space, Expression (13) is used: where d is the dimension of the solution space, x R d×1 correspond to the vector of random variables, µ R d×1 represent the center vector that contains the best solution and Σ andR d×d corresponds to the covariance matrix. If in Σ the diagonal elements (variance) are equal, and if the elements out of the diagonal (covariance) are zero, then the Gaussian distribution will create hyperspheres in the dimensional space d [38]. An easy way to calculate Σ, considering the covariance and the variances are the same, is shown in (14): where σ is the variance of the Gaussian distribution and I d×d is an identity matrix. Notice how the standard deviation of the Gaussian distribution can be defined as shown in (15): where σ 0 can also be considered to be the initial radius r 0 of the hypersphere in the dimensional space d [38]. To reach a proper exploration of the solution space, initially, σ 0 is the largest possible hypersphere. During the search process, the radius decreases to zero and contains the optimal solution [35].

Replacement of the Current Solution
In the solution stage, the best solution of s 0 C 0 (s) is selected and memorized from C 0 (s) to replace the current center of the hypersphere µ 0 . Furthermore, before the selection stage, it must be checked that the candidate solutions can be found within the limits of the solution space. For this it is employed (16), which will adjust within the limits the solutions s l k that are not found in the solution space.
where k = 1, 2, . . . , n, l = 1, 2, . . . , d and r l represent a uniform random distribution of an integer number from 0 to 1. Later, the best solution s t C t (s) is taken as the best solution, if and only if it is better than the previous solution, which will create an update in the center of the hypersphere µ t and in its radius r t , this latter being reduced. Finally, a new set of solutions C t (s) is generated around the new center.

Process of Radius Reduction
The most important process of the VSA is the adaptive adjustment of the radius of the hypersphere, which is done using a variable step approach [38]. For this, the incomplete gamma function is used, as recommended by [38]. The incomplete gamma function can be appreciated in (17).
where a > 0 is known as a shape parameter and x ≥ 0 as a random parameter. However, to determine the adaptive size of the radius in each iteration the inverse incomplete gamma function will be used, as reported in [34], shown in (18).
In MATLAB ® ,(2020a, Natick, MA, USA) the inverse incomplete gamma function for the calculation of the variable radius can be computed as observed in (19) [38].
where a t is a parameter defined as a t = 1 − t t max , t max being the total number of iterations. Additionally, the parameter x is selected as 0.1, as reported in [38].
To exemplify the proposed codification that represents the problem of optimal phasebalancing in unbalanced distribution systems, the following solution vector s k taken from the candidate set of solutions C t (s) is supposed.
This vector is formed by integer numbers that represent the load connection at each demand node (i.e., from demand node 1 to demand node m), in which the values vary from 1 to u, with u = 6, the maximum number of available connections for the loads, as shown in Section 2.
In Algorithm 1, a summary of the VSA implementation is presented, to solve the PB problem in unbalanced distribution systems.
For the solution of the PB problem, using Algorithm 1, it is recommended to set x max in 6.5 and x min in 0.5. Likewise, d is the number of demand nodes presented in the system.

Evaluation of Fitness Function
For the correct exploration and exploitation of the solution space during the iterative process of the DVSA, it is necessary to rewrite the objective function shown in (1). It is proposed that the fitness function takes the following form, shown in (20): where α is the penalization factor and v i is the nodal voltage obtained in (55) when the three-phase power flow is solved in the slave stage.
When a candidate solution is not in compliance with the voltage constraints, established in (8) to (10), the penalization factor will affect the objective function, discarding this solution as a possible best solution.

Algorithm 1 Master-slave optimization proposed.
Data: Read the information of the AC distribution network. Determine the initial center µ 0 from (12); Determine the initial radius r 0 (or the standard deviation σ 0 ) from (15); Establish the best adaptive initial function as z f (s best ) = ∞ (Minimization problem); Generate the candidate solutions using a Gaussian distribution of the center µ t with a standard deviation (radius) r t as defined in (13) to obtain the set of candidate solutions C t (s) with d columns and n rows; If C t (s) exceeds any limit, whether it is upper or lower, adjust the set of solutions within the limits using (16); Round the values of C t (s) to the closest integer; Evaluate the iterative sweep three-phase power flow method (see (55)) for each s k in C t (s) and compute the corresponding adaptive function (z f ) as shown in (20); Choose the best solution s * , contained in C t (s), that produces the lowest Keep the best solution obtained until now as s best ; end Let the center µ t+1 be equal to the best solution s best ; Update the radius r t+1 as shown in (19); t = t + 1; end Result: The best solution can be found for s best and its adaptive function is z f (s best ).

Slave Stage: Three-Phase Iterative Sweep Method
The iterative seep method, according to [40], is based on graph theory where an incidence matrix represents the network topology relating to the nodes with the branches of the system. Kirchhoff laws are also used, to calculate the currents of the system nodes, beginning from the terminal nodes towards the slack node (backward sweep), and compute the voltage drops on the network segments from the slack node towards the terminal nodes (forward sweep) [41]. This solution method is commonly known as backward-forward sweep method [42].
Let us suppose an unbalanced distribution system as shown in Figure 2, where nodes and branches are defined.  The direction of the currents flowing through the lines and the net injected currents of the system are arbitrarily defined. To represent the topology, the relation between the nodes and the phases of the system branches via the three-phase incidence matrix is used, which is shown in (21): If the current enters node i phase p 0, If branch j phase k is not connected to node i phase p (21) where i corresponds to the system nodes, p the phase of the nodes, j the system branches, and k the branches phases, for which the incidence matrix of the example shown in Figure 2 is defined in (22): The three-phase incidence matrix is divided in terms of the generation and demand as shown in (22); the voltage drop per phase through the branches is defined as the difference between the voltage nodes per phase that are connected to each other. The equations that represent the voltage drop in the branches per phase are written as follows: Writing from (23) to (28) in matrix form, (29) is obtained: Notice how (30) can be written in terms of generation and demand, as shown in (31): On the other hand, summing currents per phase in the nodes, we obtain: Node 1: Node 2: Writing (32) to (37) in matrix form, (38) is obtained: Likewise, the voltage drop per phase in the branches can be defined via the three-phase impedance matrix for distribution lines and Ohm's law [43], as shown in (40) and (41): Writing (40) and (41) in matrix form, (42) is obtained: From the theoretical development, results (31), (39), and (43) show the equations that will permit the solution of the three-phase power flow. The objective is to reach an equation of the form V d3OE = f (I d3OE ) that relates to the three-phase voltage demanded with the three-phase current demanded; moreover, it has to be taken into account that matrix A d3OE cannot be inverted, while Z r3OE can be inverted. From (43), I r3OE can be obtained and replaced in (39), reaching (44): Equation (31) is replaced in (44) obtaining (47): In addition, if we organize for V d3OE from (47), then we find (48): Moreover, isolating V d3OE from (47), Expression (48) is obtained: Furthermore, from (48) I d3OE , the three-phase demanded current, is not known. This current will depend on the type of load that is connected, i.e., Wye (Y) or Delta (∆) connection.
where i is the number of the node. For the ∆ load, Figure 4 is taken as reference, deducting (52)-(54).
Finally, using an iterative counter resolves (48), as shown in (55): where t is the iteration counter. The solution of (55) is given if and only if it is fulfilled that |V d3OE t+1 − V d3OE t | ≤ , with being the maximum acceptable error, which is recommended in the specialized literature as 1 × 10 −10 [35].
Finally, Figure 5 presents the flowchart of the three-phase iterative sweep method with which the electrical variables required to evaluate the adaptive function z f will be determined.

Test Systems
In this section, the information of the three-phase radial unbalanced distribution test systems is used to show the efficiency of the optimization method proposed in Section 3. Test systems correspond to the IEEE 8-, 25-, and 37-node radial distribution networks.

8-Node Test System
This is an unbalanced radial distribution system with 8 nodes, with the slack node located at node 1, 7 lines, and 7 demand nodes, with a nominal voltage of 11 kV [44]. The electrical configuration of this system can be appreciated in Figure 6. The total active and reactive power consumption is 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 [44]. For the purpose of this research paper, impedance matrices are modified per type of conductor, given by [44], as shown in Table A2, to properly apply PB. Load consumption per phase is shown in Appendix A in Table A1.

25-Node Test System
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 [45]. The electrical configuration of this system can be seen in Figure 7. The total active and reactive power consumption is 1073 kW and 792 kvar for phase A, 1083.3 kW and 801 kvar for phase B, and 1083.3 kW and 800 kvar for phase C. For this research, the information of load consumption per phase adopted was that proposed by the authors in [46]. The load consumption per phase and impedance matrix information per type of conductor are shown in Tables A3 and A4 from Appendix A, respectively.

IEEE 37-Node Test System
The IEEE 37-node test system is a real radial unbalanced distribution system located in California, totally composed of underground lines [47]. It has 37 nodes, node 1 being the source, 35 lines, and 25 loads. It has a transformer that operates in an unbalanced form and a voltage regulator [47]. The nominal voltage is 4.8 kV [47]. The electric configuration of this system can be appreciated in Figure 8, the nodes being numbered the same as that considered in [3]. The total active and reactive power consumption is 727 kW and 357 kvar for phase A, 639 kW and 314 kvar for phase B, and 1091 kW and 530 kvar for phase C. Initially, this system has total active power losses of 60.563 kW [47]. It is important to mention that in this research neither the transformer with unbalanced operation nor the voltage regulator is considered, as proposed by [3]. The transformer located at nodes 10 and 24 is removed, including node 24 [3]. The voltage regulator located between nodes 1 and 2 is replaced by a type 1 conductor with a length of 1850 feet, as proposed by [47]. Likewise, for the development of this work, it is established that the system loads correspond to a constant power model and the capacitive effect of the lines is neglected, due to the length and voltage level [48].
Load consumption per phase and impedance matrix information per type of conductor are shown in Tables A1-A6.

Evaluation of the Initial Configurations
Due to the modifications posed in Section 4 for the IEEE 8-, 25-, and 37-node test systems, it is necessary to determine the new operative state of the systems (nodal voltage per phase) and the active power losses for the demand conditions given in Tables A1, A3, and A5, respectively. To do this, the power flow presented in Section 3 is evaluated and the results are compared with the software DigSILENT© (Version 15.1, Gomaringen, Germany).

Computational Validation
To validate the developed three-phase iterative sweep power flow, the results obtained will be compared with DigSILENT©, the power-system analysis software. With the purpose of providing a solution to the three-phase power-flow method, MATLAB ® version 2020a is employed in a laptop with an Intel(R) Core(TM) i5-7200U CPU @2.50 GHz processor, with 8.00 GB of RAM memory running Windows 10 Home Single Language operating system of 64 bits. Notice that to determine the run-time, a maximum of 1000 iterations and a maximum tolerance of 1 × 10 −10 was considered.

8-Node Test System
In Table 4, nodal voltages are shown, obtained with the proposed three-phase iterative sweep method, using the software DigSILENT© in relation to the unbalanced three-phase power-flow solution. Moreover, in Figure 9 voltage profiles are presented for each phase, which have been obtained with the three-phase iterative sweep method using the software DigSILENT©. Table 4. Nodal voltage per phase in the 8-node test system.

MATLAB ® DigSILENT©
Node  In Table 5 the numerical results of this test system related to time processing, number of iterations, and power loss per phase of the three-phase iterative sweep method proposed and those obtained with the DigSILENT© simulation can be appreciated. The proposed three-phase iterative sweep method with an error of 1 × 10 −10 takes 5 iterations and a run-time of 0.0905 s to obtain similar results to those reached by DigSILENT©. In Figure 9 it is shown that the behavior in terms of voltage per phase is the same for both solution methods. From Table 4 it is determined that the maximum error between the data acquired by the three-phase iterative sweep method and DigSILENT© for the voltage profiles at each phase are 0.0073% for phase A, 0.0114% for phase B, and 0.0144% for phase C.
On the other hand, from Table 5 the maximum errors for active power loss are 0.0007% for phase A, 0.0436% for phase B, 0.0250% for phase C, and 0.0104% for total loss, which demonstrates that the results obtained with the three-phase iterative sweep method are reliable, since these are comparable with specialized software, such as DigSILENT©. Table 6 shows the nodal voltages per phase obtained with the three-phase iterative sweep method proposed with software DigSILENT©. In Figure 10 the voltage profiles are presented for phase A, phase B, and phase C, respectively, obtained with the three-phase iterative sweep method and DigSILENT©.   In Table 7, the numerical results are appreciated in terms of processing time, the number of iterations, and the power losses per phase of the three-phase iterative sweep method proposed, and those obtained with DigSILENT© simulation.  Tables 6 and 7, it can be confirmed that:

25-Node Test System
The three-phase iterative sweep method proposed with an error of 1 × 10 −10 takes 9 iterations and 0.1882 s of processing time to obtain similar results to those reached with the power-system software DigSILENT©. In Figure 10, it is shown that the behavior in terms of voltage per phase is the same for both solution techniques. From Table 6 it is determined that the maximum error of the results obtained between the three-phase sweep iterative method and DigSILENT© for the voltage profiles at each phase are 0.0098% for phase A, 0.0104% for phase B, and 0.0051% for phase C, respectively. By the other side from Table 7 the maximum error for the active power loss per phase are 0.0115% for phase A, 0.0015% for phase B, 0.0071% for phase C, and 0.0032% for the total loss, which demonstrates that the results obtained with the developed three-phase iterative sweep method are reliable, since they can be compared with the specialized software for power systems DigSILENT©.

IEEE 37-Node Test System
In Table 8 the nodal voltage can be appreciated per phase obtained with the threephase iterative sweep method and software DigSILENT©.  In Figure 11, the voltage profiles are presented for phases A, B, and C, obtained with the three-phase iterative sweep method and DigSILENT©.
In Table 9, the numerical results can be appreciated, with the processing time, number of iterations and power losses per phase of the proposed three-phase iterative sweep method and the specialized software DigSILENT©.  From the results obtained in Tables 8 and 9 it is observed that: The proposed three-phase iterative sweep method with an error of 1 × 10 −10 takes 9 iterations and a processing time of 0.2919 s to obtain similar results to those reached with the specialized software DigSILENT©. In Figure 11, it is shown that the behavior in terms of voltage per phase is the same for both solution methods. From Table 8 it is determined that the maximum error for the voltage profiles obtained between the three-phase iterative sweep method and DigSILENT© at each phase are 0.0068% for phase A, 0.091% for phase B, and 0.0110% for phase C. On the other hand, from Table 8 the maximum error for the active power losses are 0.0003% for phase A, 0.0185% for phase B, 0.0131% for phase C, and 0.0038% for the total loss. This demonstrates that the results obtained with the developed three-phase iterative sweep method are reliable, since these are comparable with the specialized power-system software DigSILENT©.

Results Obtained: Methodology Proposed
With knowledge of the new operative state and the power losses per phase and total losses of the IEEE 8-, 25-, and 37-node test systems, the master-slave optimization is performed to determine the new set of connections for the demand nodes, the reduction of active power losses and reduction on the unbalanced percentage of the test systems proposed.

Computational Validation
To solve the MINLP optimization model shown in (1) to (11) that represents the optimal balancing problem of phases for distribution systems, MATLAB ® V2020a is employed in a laptop with an Intel(R) Core(TM) i5-7200U CPU @2.50 GHz processor, with 8.00 GB of RAM memory running Windows 10 Home Single Language operating system of 64 bits.
To implement the proposed master-slave optimization for the PB problem in unbalanced distribution systems, the information presented in Table 10 is used.

Number of evaluations 100
Notice that for the selection of the parameters shown in Table 10 a sensitivity analysis was done between the number of candidate solutions and the number of iterations employed by the DVSA for the proposed cases in Section 4. For that, the solution space was explored between 8 and 15 candidate solutions, and from 500 to 2000 iterations, which confirmed the proposal of [35]: (i) for a reduced number of iterations and candidate solutions, the optimization process is faster but the chance of reaching an optimal global solution significantly decreases; (ii) for a large number of candidate solutions and iterations, the optimization process is slower, but there is greater chance of obtaining an optimal global solution; and (iii) several iterations between 500 and 1000 with several candidate solutions between 8 and 12 allows a balance between processing time and the possibility of finding the optimal global solution. Therefore, the number of candidate solutions is selected as 10, and the number of iterations is set as 800, for the parameters in the proposed DVSA.
On the other hand, from this sensitivity analysis it was observed that for the number of candidate solutions and selected iterations, after reaching a certain number of iterations, the possible global optimal solution does not change in value. Due to this, it is proposed that after 250 iterations the search process be terminated, when the objective function value is not modified, and the best solution of the current population is reported. This will allow a decrease of the processing time of the algorithm.
Due to the candidate solutions being generated through a random process using a normal distribution (i.e., Gaussian distribution) , it is recommended to do 100 of consecutive evaluations to determine the average processing time algorithm and the number of times the optimal solution of the problem is reached.
Finally, for comparative purposes, we contrasted the proposed discrete VSA approach for the solution of the phase-balancing problem with the classical Chu & Beasley genetic algorithm (CBGA) [3] using the same size of the initial populations and the number of generations assigned for the VSA for all the test feeders.

8-Node Test System
In Table 11, the results obtained for the DVSA are shown, for the 8-node tests system. Since there were modifications in the benchmark case, the results reported by the set of connections are compared with CBGA and the specialized power-system software DigSILENT©. From Table 11, the following situations can be observed: (i) the proposed DVSA reduces the total power loss by 24.34%, i.e., from 13.9925 kW to 10.5869 kW. Moreover, it is observed that the power loss is increased by 59.08% for phase A and 75.74% for phase B, but is reduced in 62.1795% for phase C, making the loss per phase to be similar to each other, reducing the system-unbalancing; (ii) the proposed DVSA reaches an optimal global solution by connection change only at 3 demand nodes of the system, i.e., nodes 2, 4, and 6, which represent 42.86% of the demand nodes; and (iii) the errors for the active power loss are 0% for phases A, B, and C, respectively, and the total loss in comparison with CBGA; in the same way, the maximum errors for the active power loss are 0.0052% for phase A, 0.0005% for phase B, 0.0051% for phase C, and 0.0003% for total loss, respectively, in comparison with DigSILENT©, which validates the connections for the loads obtained by the DVSA, since the results are comparable with CBGA and the specialized power-system software DigSILENT©.
In Figure 12, the active and reactive power unbalance percentage is presented for the modified benchmark case and those obtained by the DVSA. In Figure 12a, a reduction on the active power unbalance can be appreciated of 5.51% for phase A, 30.03% for phase B, and 40.36% for phase C. In Figure 12b, a reduction of the reactive power unbalance can be appreciated of 5.87% for phase A, 29.99% for phase B, and 40.36% for phase C. Likewise, the results obtained for the active and reactive power unbalance per phase are very similar to each other. Figure 13a shows the voltage profiles per phase for the modified benchmark case, while Figure 13b shows the PB effect with the methodology proposed for the 8-node test system. Graphically, it can be appreciated that when reducing the phase unbalance, the voltage profile behaviors are similar. Additionally, before performing PB, the lowest value of voltage was 0.9923 pu, corresponding to node 4 at phase C; after the PB, this value is 0.9960 pu, which indicates an improvement of 0.3730%.
On the other hand, after doing 100 consecutive evaluations of the DVSA proposed for the 8-node test system, it is found that the best solution is 10.5869 kW, which is repeated in 92% of the evaluations performed, while CBGA got the same solution but in 70% of the evaluations. It is also observed that only one evaluation obtains this solution in one iteration with DVSA proposed and only one evaluation obtains this solution in three iteration with CBGA. Likewise, for the 8-node test system, the standard deviation for the DVSA proposed

25-Node Test System
In Table 12, the results obtained are shown for the DVSA proposed in the 25-node test system, and the respective comparison with CBGA and the software DigSILENT©.  Table 12, the following situations can be observed: (i) the proposed DVSA reduces the total power losses by 4.1526%, i.e., from 75.4207 kW to 72.2888 kW. It is observed that the power loss is reduced by 30.41% for phase A, and 13.8654% for phase C; however, there is an increase of 76.9604% for phase B, making the power losses per phase similar, reducing the system unbalance; (ii) in comparison with the CBGA, the proposed DVSA reduces the total power loss by 0.005%; (iii) the proposed DVSA reaches a global optimal solution via connection changing on 21 demand nodes of the system, which represents 87.50% of the demand nodes; and (iv) it can be appreciated that the maximum errors for the active power loss are 0.0003% for phase A, 0.0110% for phase B, 0.0144% for phase C, and 0.0002% for the total loss, respectively, which validates the load connections for the DVSA, since the results are comparable with the specialized power-system software DigSILENT©.
In Figure 14 the active and reactive power unbalancing percentage is presented for the modified benchmark case and those obtained for the DVSA. In Figure 14a, a reduction of the active power unbalance by 22.86% can be appreciated for phase A and 24.51% for phase B, and an increase of 0.34% for phase C. In Figure 14b a reduction of the reactive power unbalance is shown of 17.60% for phase A, 19.53% for phase B, and an increase of 1.01% for phase C. Although the unbalanced percentage increases for phase C, for both active and reactive power, this increase can neglected in comparison with the unbalanced percentage reduction of the phases A and B for active and reactive power. In the same sense, the results obtained for the active and reactive power unbalance per phase are very similar to each other. Figure 15a shows the voltage profiles per phase for the modified benchmark case, while Figure 15b shows the PB effect with the methodology proposed for the 25-node test system. Graphically it can be noted that when there is a reduction of unbalance for phases A and B, the voltage profiles have similar behavior; conversely, when there is an unbalance for phase C, the voltage profiles diverge. Additionally, before PB, the lowest voltage value was 0.9352 pu, corresponding to node 13 in phase A; after the PB this value is 0.9472 pu, which indicates an improvement of 1.28%.
On the other hand, after doing 100 consecutive evaluations of the proposed DVSA for the 25-node test system, it is found that the best solution is 72.2888 kW, which is repeated only once in the evaluations performed. It is also observed that this value is obtained in 395 iterations. Although CBGA found that the best solution is 72.2919 kW, which is repeated only once in evaluations performed, it is also observed that this value is obtained in 550 iterations. Likewise, it is found that for the 25-node test system the standard deviation of the proposed DVSA is 0.0233 kW for the objective function and 126.71 for the iterations performed; moreover, the standard deviation of the CBGA is 0.

37 Nodes Test System
In Table 13, the results obtained for the proposed DVSA are found for the IEEE 37-node test system and the comparison with CBGA and the specialized power-system software DigSILENT©. From Table 13, the following situations can be observed: (i) the proposed DVSA reduces the total power losses by 19.2493%, i.e., from 76.136 kW to 61.4801 kW. Moreover, it is observed that the power losses are reduced by 22.4195% for phase A and 49.5108% for phase C; however, there is an increase of 82.1248% for phase B, making the losses per phase similar to each other, reducing the system unbalance; (ii) in comparison with the CBGA, the proposed DVSA reduces the total power losses by 0.16%; (iii) the proposed DVSA reaches a global optimal solution, making the change of connection of 25 demand nodes of the system, which represents 71.43% of the demand nodes; and (iv) it is appreciated that the largest error for the active power loss is 0.0058% for phase A, 0.0071% for phase B, 0.0008% for phase C, and 0.0001% for the total loss, which validates the load connections obtained by the DVSA, since the results are comparable to the specialized power-system software DigSILENT©.
In Figure 16 the unbalanced percentage of active and reactive power is presented for the modified benchmark case and those obtained by the DVSA. In figure 16a, a reduction of the active power unbalance is appreciated of 4.40% for phase A, 6.11% for phase B, and 24.18% for phase C. In Figure 16b, a reduction of the reactive power unbalance is appreciated of 3.50% for phase A, 5.41% for phase B, and 23.56% for phase C. Figure 17a shows the voltage profiles per phase for the modified benchmark case, while Figure 17b shows the PB effect with the methodology proposed for the IEEE 37-odes test system. Graphically, it is shown that when the phase unbalance is reduced, the voltage profiles have a similar behavior. Additionally, before PB, the lowest voltage value was 0.9369 pu, corresponding to node 22 in phase A; after PB, this value is 0.9593 pu, which shows an improvement of 2.39%.
On the other hand, after doing 100 consecutive evaluations of the proposed DVSA for the IEEE 37-node test system, it is found that the best solution is 61.4801 kW, which is repeated only once in the evaluations performed. It is also observed that this value is obtained in 387 iterations. Although CBGA found that the best solution is 61.5785 kW, which is repeated only once in the evaluations performed, it is also observed that this value is obtained in 427 iterations. Likewise, for this test system, the standard deviation of the proposed DVSA is 0.3286 kW for the objective function, and 149.68 for the iterations performed; by comparison with CBGA, the standard deviation is 0.4274 kW for the objective function and 532.934 for the iterations performed. Regarding the processing time, the proposed DVSA takes 50.0262 s and CBGA takes on average 14.1816 s to reach the global optimum.
All the solutions obtained by the DVSA proposed for the IEEE 8-, 25-, and 37-node test systems require a longer computing time compared to the results obtained with the CBGA. However, when the proposed methodology is applied, an improvement in the results is evident compared with the CBGA, it can be observed that the proposed DVSA requires a smaller number of iterations to reach the best solution.

Conclusions and Future Works
In this research paper, an optimization algorithm is presented to solve the optimal phase-balance problem in radial and meshed unbalanced distribution systems. This algorithm uses a master-slave optimization methodology. In the master stage, a discrete version of the Vortex Search Algorithm is used, which determines the set of phase connections for the demand nodes of the system; in the slave stage, a three-phase version of the classic method of iterative sweep power flow is implemented, which determines the operative state of the system including the active power total losses. The objective function under analysis is focused on active power-loss minimization in unbalanced distribution systems. The numerical results show the applicability and the efficiency of the optimization method developed for the proposed test systems. In the 8-node test system a 24.34% reduction in total power loss is reached; in the 25nodes system the reduction accounts for 4.1526% of the total power loss, and in the IEEE 37-node test system this reduction is 19.2483% of total power loss. In the same manner, in the methodology proposed, a reduced processing time is shown to determine the global optimal solution. For the 8-node system, the average processing time is 6.059 s, for the 25-node test system the average processing time is 39.69 s, while for the IEEE 37-node test system the processing time is 50.026 s, on average.
One of the main characteristics of the developed DVSA is that a discrete codification is employed, using integer numbers as decision variables, which allows the representation of the problem of optimal phase-balancing. This avoids the excessive use of large-size binary vectors, and consequently the processing time is reduced. Additionally, after 100 consecutive evaluations, the proposed optimization algorithm reaches a possible global optimal solution for the 8-node test system in one iteration, for the 25-node system in 395 iterations, and for the 37-node test system in 387 iterations, which confirms that a sensitivity analysis for setting the DVSA allows the reaching of a global optimal solution with a minimum of iterations. Likewise, it is appreciated that all the solutions obtained by the DVSA for the IEEE 8-, 25-, and 37-node test systems are similar to each other, since the standard deviation is 0.4 W, 0.023 kW, and 0.32 kW, respectively. Finally, it is noted that the application for optimal phase-balancing in unbalanced distribution systems reduces the unbalanced percentage, which leads to an improvement of the system voltage profiles and therefore the total power loss is reduced accordingly, in such a way that the system finds its best operation point.
Numerical comparisons with the CBGA demonstrated that the proposed discrete VSA approach can find the best objective function values when the same number of iterations and population sizes are used during the comparison process. It was observed that the CBGA required lower processing times to reach the solution of the problem when compared with the proposed VSA; however, for planning purposes, computational times lower than 60 s reported by the proposed approach can be considered extremely efficient (i.e., negligible) for solving the problem of phase-balancing in distribution grids.
In future works, it will be possible to address the following problems: (i) extend the problem presented in this work with typical demand curves using the master-slave methodology to reduce the energy losses through the determination of the set of phase connections for the demand nodes of the system; (ii) use the proposed master-slave methodology to determine the siting and sizing of AC three-phase distributed generators to balance the system using a discrete codification of integer number; and (iii) employ the proposed master-slave methodology to determine the optimal wire gauge per phase for each network segment of an unbalanced three-phase distribution system and, accordingly, reduce total active power loss.

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.

Appendix A. Electrical Parameters
Tables A1, A3, and A5 show the parameters for the loads and lines of the IEEE 8-, 25-, and 37-node test systems, respectively. Tables A2, A4 and A6 show the impedance matrices for the type of conductor of the IEEE 8-, 25-, and 37-node test systems, respectively.