*Article* **An Improved Crow Search Algorithm Applied to the Phase Swapping Problem in Asymmetric Distribution Systems**

**Brandon Cortés-Caicedo <sup>1</sup> , Laura Sofía Avellaneda-Gómez <sup>1</sup> , Oscar Danilo Montoya 2,3,\* , Lázaro Alvarado-Barrios <sup>4</sup> and César Álvarez-Arroyo <sup>5</sup>**


**Citation:** Cortés-Caicedo, B.; Avellaneda-Gómez, L.S.; Montoya, O.D.; Alvarado-Barrios, L.; Álvarez-Arroyo, C. An Improved Crow Search Algorithm Applied to the Phase Swapping Problem in Asymmetric Distribution Systems. *Symmetry* **2021**, *13*, 1329. https:// doi.org/10.3390/sym13081329

Academic Editors: Juan Alberto Rodríguez Velázquez and Alejandro Estrada-Moreno

Received: 6 July 2021 Accepted: 21 July 2021 Published: 23 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Abstract:** This paper discusses the power loss minimization problem in asymmetric distribution systems (ADS) based on phase swapping. This problem is presented using a mixed-integer nonlinear programming model, which is resolved by applying a master–slave methodology. The master stage consists of an improved version of the crow search algorithm. This stage is based on the generation of candidate solutions using a normal Gaussian probability distribution. The master stage is responsible for providing the connection settings for the system loads using integer coding. The slave stage uses a power flow for ADSs based on the three-phase version of the iterative sweep method, which is used to determine the network power losses for each load connection supplied by the master stage. Numerical results on the 8-, 25-, and 37-node test systems show the efficiency of the proposed approach when compared to the classical version of the crow search algorithm, the Chu and Beasley genetic algorithm, and the vortex search algorithm. All simulations were obtained using MATLAB and validated in the DigSILENT power system analysis software.

**Keywords:** improved crow search algorithm; normal Gaussian distribution; phase swapping problem; power losses; asymmetric distribution grids; vortex search algorithm

## **1. Introduction**

Due to the economic and population growth, the dependence on electrical systems has equally grown to satisfy humanity's basic needs, changing the habits and customs of how individuals live and work [1]. To ensure this, three-phase distribution networks are used, which are responsible for interconnecting transmission and sub-transmission networks with end-users (i.e., residential, industrial, and commercial areas) requiring medium and low voltage [2,3]. These systems generally operate in an asymmetric manner due to the following factors. (i) The configurations on the distribution lines are asymmetrical since the transposition criterion is not applicable due to the short length of the lines [4,5]. (ii) The nature of the loads may be 1*ϕ*, 2*ϕ*, or 3*ϕ*, which generates unbalances in voltages at the nodes and in the line currents [6]. (iii) The arbitrary location of single-phase transformers on the phases of the system causes an unbalance in the currents through the lines [7]. Load unbalances in distribution systems create undesirable scenarios such as the increase of current in any phase system, which produces an increase in power losses through its constituent elements [8]. These power losses can exceed the capacity required to supply the demand, cause equipment to age, and increase investment and operating costs for network operators [7,9].

The importance of reducing power losses in distribution networks has established multiple approaches, such as (i) optimal placement and sizing of distributed generation [10], (ii) optimal capacitor placement and sizing [11], (iii) optimal network reconfiguration [12], (iv) optimal conductor sizing in distribution networks [13], (v) optimal power system restoration [14], and (vi) optimal phase swapping [15,16]. These strategies can significantly help distribution companies to reduce the number of power losses. However, the first two approaches involve significant investments since they integrate new devices into the distribution network [17]. The third approach requires less investment since few distribution lines need to be constructed to realize the optimal network reconfiguration [18]. The fourth approach also requires high investment as the system conductors need to be renewed [19]. The fifth methodology is more appropriate for operation of the power system after fault isolation. The sixth strategy is the most economical as it requires few teams to reconfigure the system loads without investing in new equipment [20]. Bearing in mind the low phase swapping costs to minimize the power losses in ADSs, a new master–slave optimization strategy is proposed to solve this problem.

In the specialized literature, the balance phase problem, with the minimizing power losses approach, has been solved using different optimization methods, including the Chu and Beasley genetic algorithms [8,16,21–24], particle swarm optimization [9], mixedinteger convex optimization [25], bat optimization algorithm [26], differential evolution algorithm [27], simulated annealing optimizer [28], and vortex search algorithm [15], among others.

The main feature of the optimization methodologies described above is that they employ the master–slave optimization scheme to solve the problem [15]. The master defines the connection of the loads to the nodes. The slave is typically a power flow tool that allows one to revise and exploit the solution space through the power losses calculation [16].

Similar to the metaheuristic optimization methods described above, a master–slave methodology is proposed in this work to solve the phase swapping problem in ADSs. The proposed optimization algorithm corresponds to an improved version of the crow optimization algorithm (CSA) to select the connection of the loads in the master stage, together with the use of the iterative sweep power flow method in its three-phase version in the slave stage. In the master stage, the connection of each load is defined using an integer encoding between 1 and 6, which represents the six possible connection forms for a three-phase charge [8]. The slave stage is responsible for evaluating the power flow to determine the total power losses for the connection set provided in the master stage [15]. Improvements in the classical CSA are carried out in the crow avoidance stage based on a probability criterion [29]. If the probability is higher than the crow knowledge probability (*Ap*), the new crow position *i* is provided using the classical CSA exploration proposed in [29]. Likewise, if the possibility is less than the crow knowledge probability (i.e., *Ap*), the new crow position *i* is generated through a regular Gaussian distribution (GD) used in the process of evolution of the vortex search algorithm (VSA) [30]. The main benefit of the VSA is that the solution space can be explored and exploited through the use of hyper-spheres derived from the selection of an individual from the current population. It also clarifies that the criterion of evolution in our proposal is applied at each iteration, which implies that this process is of the adaptive type.

The main contributions of our proposal are listed below.


It is relevant to mention that, upon analyzing the specialized literature, there was no evidence of the CSA application to the phase swapping problem in distribution systems, which corresponds to a research gap that this work intends to fill. In addition, the numerical results obtained in test systems of 8, 25, and 37 nodes prove the quality of the algorithm when compared with classical metaheuristic optimization methodologies.

The structure of the remainder of the document takes the following form: Section 2 presents the optimal phase swapping problem representation in ADSs, Section 3 presents the ICSA incorporated with the TPPF method, Section 4 portrays the electrical networks used in this research, and Section 5 represents the obtained results for the connections set and the grid power losses. Finally, Section 6 states the conclusions drawn from the development of this article.

#### **2. Mathematical Formulation**

The optimal phase swapping problem in ADSs is represented by a mixed-integer nonlinear programming (MINLP) model [16]. The binary variables are the decision variables, which correspond to the connections set for each load presented in the system [8]. Additionally, the power flow formulation provides the continuous part of the decision variables. The nonlinear nature of the products appears between the different voltage magnitudes at the nodes and the trigonometric functions [25]. Next, the objective function and the set of constraints representing the phase swapping problem are presented.

#### *2.1. Formulation of the Objective Function*

The phase swapping problem has an objective function associated with the minimization of total active power losses of the ADS, as presented in Equation (1):

$$\min \ z = \sum\_{n \in \mathcal{N}} \sum\_{m \in \mathcal{N}} \sum\_{f \in \mathcal{F}} \sum\_{g \in \mathcal{F}} Y\_{nfmg} V\_{nf} V\_{mg} \cos(\delta\_{nf} - \delta\_{mg} - \theta\_{nmfg}), \tag{1}$$

where *z* defines the value of the objective function. Further, *Yn f mg* is the admittance magnitude associated with node *n* in the electrical phase *f* with node *m* in the electrical phase *g*, *Vn f*(*Vmg*) corresponds to the voltage magnitude at node *n*(*m*) in the electrical phase *f*(*g*), *δn f*(*δmg*) represents the angle of the voltage at node *n*(*m*) in the electrical phase *f*(*g*), and *θnm f g* represents the admittance angle associated with node *n* in the electrical phase *f* with node *m* in the electrical phase *g*. It is relevant to mention that F and N are the sets containing all phases and nodes, respectively.

**Remark 1.** *The product between the magnitudes of the voltages and the trigonometric functions makes the objective function nonlinear and nonconvex [16]. The structure of the objective function makes advanced numerical optimization techniques necessary to minimize it efficiently [15]. A master–slave methodology is proposed, as is the case of the improved version of the developed CSA, due to its simplicity in programming terms.*

#### *2.2. Set of Constraints*

The phase swapping problem has a set of constraints that corresponds to the different operating limitations in an ADS [15]. These are shown from Equation (2) to Equation (6):

$$P\_{nf}^{\mathrm{s}} - \sum\_{\mathcal{g} \in \mathcal{F}} \mathbf{x}\_{nf\mathcal{g}} P\_{ng}^{d} = V\_{nf} \sum\_{m \in \mathcal{N}} \sum\_{\mathcal{g} \in \mathcal{F}} Y\_{nf\mathsf{m}\mathsf{g}} V\_{mg} \cos(\delta\_{nf} - \delta\_{\mathsf{m}\mathsf{g}} - \theta\_{nf\mathsf{m}\mathsf{g}}), \begin{cases} \forall f \in \mathcal{F} \\ \forall n \in \mathcal{N} \end{cases} \tag{2}$$

$$\mathbf{Q}\_{\mathrm{nf}}^{\mathrm{s}} - \sum\_{\mathcal{g} \in \mathcal{F}} \mathbf{x}\_{\mathrm{nf}} \mathbf{Q}\_{\mathrm{ng}}^{\mathrm{d}} = V\_{\mathrm{nf}} \sum\_{m \in \mathcal{N}} \sum\_{\mathcal{g} \in \mathcal{F}} Y\_{\mathrm{nfm}g} V\_{\mathrm{mg}} \mathrm{sin}(\delta\_{\mathrm{nf}} - \delta\_{\mathrm{mg}} - \theta\_{\mathrm{nfm}g}), \begin{cases} \forall f \in \mathcal{F} \\ \forall n \in \mathcal{N} \end{cases} \,, \tag{3}$$

$$\sum\_{\mathcal{g}\in\mathcal{F}}\mathfrak{x}\_{nf\mathcal{g}}=\mathfrak{1},\{\forall f\in\mathcal{F},\quad\forall n\in\mathcal{N}\},\tag{4}$$

$$\sum\_{f \in \mathcal{F}} \mathbf{x}\_{nfg} = \mathbf{1}\_{\prime} \{ \forall g \in \mathcal{F}\_{\prime} \quad \forall n \in \mathcal{N} \},\tag{5}$$

$$V\_{\min} \le V\_{nf} \le V\_{\max} \left\{ \forall \mathcal{g} \in \mathcal{F}\_{\prime} \quad \forall n \in \mathcal{N} \right\}\_{\prime} \tag{6}$$

where *P s n f* is the variable associated with the active power produced at generator *s* connected to node *n* in the electrical phase *f* , and *Q<sup>s</sup> n f* is the generated reactive power at

generator *s* sited at node *n* in the electrical phase *f* . *P d ng* indicates the active demanded power at node *n* in the electrical phase *g*, while *Q<sup>d</sup> ng* describes the reactive power needed at node *n* in the electrical phase *g*. *xn f g* is a binary variable defining the configuration of demand node *n* at *f* in the electrical phase *g*. Finally, *V*min and *V*max correspond to the allowable limits of voltage regulation for all the system nodes.

**Remark 2.** *Equations (2) and (3) are nonlinear and nonconvex. These features highlight the complexity of the TPPF problem for electrical systems, making it necessary to use numerical methods, such as the iterative sweep method, to solve it.*

Note that Equation (1) defines the form of the objective function for the phase swapping problem formulated as the minimization of power losses under a given demand condition in all network sections of the system. Equations (2) and (3) define the apparent power balance constraints maintained at each phase and node of the ADS. Equations (4) and (5) ensure that loads take a unique connection form by using a matrix of connections (i.e., *xn f g*) at each node [5]. Finally, the constraint presented in (6) defines the allowable limits of voltage regulation for all nodes of the system [15].

#### **3. Methodology Proposed**

To solve the optimal phase swapping problem in ADS with the objective of minimize power losses for a specific demand condition, this paper proposes to use an ICSA [29] as the master stage in conjunction with the iterative swept TPPF as the slave stage [15]. The master stage defines the phase configurations at each system demand node to achieve the most balanced system possible, while the slave stage evaluates the power flow constraints defined in (2) and (3). The section below will describe each component of the proposed master–slave methodology.

#### *3.1. Slave Stage: TPPF Method*

The iterative sweep power flow method is a numerical method typically used for single-phase distribution networks [31]. Nevertheless, this method has been adapted for three-phase ADSs with wye (i.e., *Y*) and delta (i.e., ∆) loads [15]. This method is derived from graph theory, where the topology of the network is represented by an incidence matrix, which relates the nodes and the links of the system [31]. First, Kirchhoff's first law is used to calculate currents in the system nodes, starting from the terminal nodes and until the source node, which corresponds to the implementation of the backward sweep stage. Then, Kirchhoff's second law, which corresponds to the implementation of the forward sweep stage, is used to calculate the voltage drops in the network sections from the slack node to the terminal nodes [31].

One of the most important aspects of the iterative sweep power flow method is that it is derivative-free. Likewise, the matrices involved in the calculations are constant, which implies that the computing times required to obtain a solution are in milliseconds [15].

In order to expose the iterative swept TPPF method developed by [15], in any *n*−node system, we used the relationship between the nodal voltage and the injected current that is presented using the equivalent between the admittance matrix and the incidence matrix [32], as shown in Equation (7).

$$
\begin{bmatrix}
\mathbb{I}\_{\mathbb{S}^{3\varphi}} \\
\mathbb{I}\_{d\mathbb{S}\varphi}
\end{bmatrix} = \begin{bmatrix}
\mathbf{A}\_{\mathbb{S}^{3\varphi}} \mathbb{Z}\_{r\mathbb{S}\varphi}^{-1} \mathbf{A}\_{\mathbb{S}^{3\varphi}}^{T} & \mathbf{A}\_{\mathbb{S}^{3\varphi}} \mathbb{Z}\_{r\mathbb{S}\varphi}^{-1} \mathbf{A}\_{d\mathbb{S}\varphi}^{T} \\
\mathbf{A}\_{d\mathbb{S}\varphi} \mathbb{Z}\_{r\mathbb{S}\varphi}^{-1} \mathbf{A}\_{\mathbb{S}^{3\varphi}}^{T} & \mathbf{A}\_{d\mathbb{S}\varphi} \mathbb{Z}\_{r\mathbb{S}\varphi}^{-1} \mathbf{A}\_{d\mathbb{S}\varphi}^{T}
\end{bmatrix} \begin{bmatrix}
\mathbb{V}\_{\mathbb{S}^{3\varphi}} \\
\mathbb{V}\_{d\mathbb{S}\varphi}
\end{bmatrix},
\tag{7}
$$

where V*g*3*<sup>ϕ</sup>* is the vector containing all the voltages at the slack node that are known for power flow purposes [21]. V*d*3*<sup>ϕ</sup>* is the vector containing all the unknown variables of interest, i.e., the demand voltages. Further, I*g*3*<sup>ϕ</sup>* represents the vector with the net current injections at the slack node, I*d*3*<sup>ϕ</sup>* represents the vector that involves all the currents at the nodes of consumption, and Z*r*3*<sup>ϕ</sup>* is the matrix that contains all the impedance matrices of the distribution lines present in the system. **A***g*3*<sup>ϕ</sup>* is the component of the incidence matrix that associates the source nodes to each other, while **A***d*3*<sup>ϕ</sup>* is the component of the incidence matrix relating the nodes of consumption to each other.

**Remark 3.** *The voltage variables and current ones in Equation (7) are organized by nodes and phases according to the three-phase condition.*

From Equation (7), it can be observed that the second row has the nodal voltages at the demand nodes (i.e., V*d*3*ϕ*), which are the unknown variables in power flow studies [16]. Equation (7) can then be rewritten as follows: (8).

$$\mathbb{V}\_{d\Im\phi} = -\mathbb{Y}\_{d d\Im\phi}^{-1} \left[ \mathbb{Y}\_{d\lg\Im\phi} \mathbb{V}\_{\mathbb{S}^{\Re}} - \mathbb{I}\_{d\Im\phi} \right],\tag{8}$$

where Y*dg*3*<sup>ϕ</sup>* = **A***d*3*ϕ*Z −1 *r*3*ϕ***A***<sup>T</sup> g*3*ϕ* and Y*dd*3*<sup>ϕ</sup>* = **A***d*3*ϕ*Z −1 *r*3*ϕ***A***<sup>T</sup> d*3*ϕ* .

Equation (8) allows the determination of all the nodal demand voltages per phase. However, it is necessary to consider the type of load, either *Y* or ∆, to establish the demand current (i.e., I*d*3*ϕ*).

In the case where node *m* has a constant power load with a *Y* structure (assuming it is solidly earthed [33]), the demand current can be shown as in Equation (9), as reported in [21].

$$\mathbb{I}\_{dm\Im\rho} = -\mathbf{diag}^{-1}(\mathbb{V}\_{dm\Im\rho}^\*)\mathbb{S}\_{dm\Im\rho}^\* \tag{9}$$

If node *m* has a load with a connection ∆, the demand current can be as shown in Equation (10), as reported in [21].

$$\mathbb{I}\_{dm3q} = - (\mathbf{diag}^{-1}(\mathbf{M}\mathbb{V}\_{dm3q}^\*) - \mathbf{diag}^{-1}(\mathbf{M}^T\mathbb{V}\_{dm3q}^\*)\mathbf{H})\mathbb{S}\_{dm3q\prime}^\* \tag{10}$$

where the **H** and **M** matrices are defined as follows:

$$\mathbf{H} = \begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}, \mathbf{M} = \begin{bmatrix} 1 & -1 & 0 \\ 0 & 1 & -1 \\ -1 & 0 & 1 \end{bmatrix}$$

Through a *t* iteration counter, the solution of Equation (8) is obtained if max<sup>n</sup> ||V *t*+1 *d*3*ϕ* | − |V*<sup>t</sup> d*3*ϕ* ||o <sup>≤</sup> *<sup>e</sup>*, where *<sup>e</sup>* is the maximum tolerance suggested as <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>10</sup> [34].

When solving the TPPF, the main objective is to evaluate the power losses for the phase connection set established in the master stage. For this purpose, Equation (11), as described in [32], is used.

$$\mathbb{P}\_{\text{loss}} = \text{real}\left\{ \left| \mathbf{J}\_{3\varphi} \right|^{T} \mathbb{Z}\_{r \mathbf{3}\varphi} \left| \mathbf{J}\_{3\varphi} \right| \right\} \,, \tag{11}$$

where Ploss describes the total system effective power losses and **J**3*<sup>ϕ</sup>* represents the current per phase flowing through the system branches expressed as shown in Equation (12) through Ohm's Law, as reported in [15].

$$\mathbf{J}\_{3\boldsymbol{\varphi}} = \mathbb{Z}\_{r3\boldsymbol{\varphi}}^{-1} \mathbf{E}\_{3\boldsymbol{\varphi}}.\tag{12}$$

where **E**3*<sup>ϕ</sup>* represents the voltage drop per phase in the system branches, which can be written in terms of the generation and demand using the three-phase incidence matrix as shown in Equation (13), as reported in [15].

$$\mathbf{E}\_{3\varphi} = \mathbf{A}\_{\mathcal{g}\mathfrak{Z}\varphi}^T \mathbb{V}\_{\mathcal{g}\mathfrak{Z}\varphi} + \mathbf{A}\_{d\mathfrak{Z}\varphi}^T \mathbb{V}\_{d\mathfrak{Z}\varphi} \tag{13}$$

Algorithm 1 shows the general implementation of the TPPF method by an iterative sweep for ADSs with connected loads in *Y* and ∆.

```
Algorithm 1: Solution of the TPPF problem for ADSs with Y and ∆ loads
  Define the characteristics of the unbalanced three-phase system under study;
  Obtain the per-unit equivalent of the system;
  Generate the 3ϕ incidence matrix A3ϕ;
  Extract the components Ag3ϕ and Ad3ϕ;
  Define Zr3ϕ;
  Calculate Ydg3ϕ and Ydd3ϕ;
  Select tmax;
  Define e;
  Chose the voltages per phase of the Slack node: Vg3ϕ =
                                                            -

                                                             1∠0, 1∠ − 2π
                                                                           3
                                                                             , 1∠
                                                                                   2π
                                                                                    3
                                                                                      T
                                                                                        ;
  Do t = 0;
  for m ≥ n − 1 do
      Do Vt
           dm3ϕ = Vg3ϕ;
  end
  for t ≤ tmax do
      Define m = 1;
      for m ≥ n − 1 do
         if node load m is connected at Y then
             Calculate I
                        t
                        dm3ϕ
                             using Equation (9);
         else
             Calculate I
                        t
                        dm3ϕ
                             using Equation (10);
         end
      end
      Calculate the new voltages at the demand nodes Vt
                                                          d3ϕ
                                                              using Equation (8);
      if maxn
              ||V
                 t+1
                 d3ϕ
                     | − |Vt
                           d3ϕ
                              ||o
                                  ≤ e then
         Report the nodal voltages as V3ϕ =
                                               -

                                                Vg3ϕ, Vd3ϕ
                                                            T
                                                               ;
         Calculate the voltage drop across the branches of the system using Equation (13);
         Calculate the current flowing in the system branches using Equation (12);
         Calculate the power losses using Equation (11) ;
         break;
      else
         Do Vt
               d3ϕ = V
                       t+1
                       d3ϕ
                          ;
      end
  end
```
**Remark 4.** *The convergence of the matrix iterative sweep method can be demonstrated with the Banach fixed-point theorem, as reported in [31]. So, if the system is far enough from the stress collapse point, it can be guaranteed that the solution of Equations (2) and (3) obtain any combination of nodal loads provided by the master stage.*

#### *3.2. Master Stage: ICSA*

The master stage is responsible for providing the nodal connections set for evaluation in the iterative sweep TPPF presented in Algorithm 2. This paper proposes an improved version of the classical CSA modifying the solution space exploration by introducing a normal GD employed by the VSA optimization method [35]. Before explaining the improvements made to the CSA, the encoding of the phase swapping problem in ADSs according to the possible configurations is shown in Table 1.


**Table 1.** Possible connection types for the system loads [8].

The coding used to represent an individual *i* in iteration *k* is presented as follows:

$$X\_{i,k} = [6, \text{ 2, 4, 5, ..., c, ..., 1}]$$

Here, the dimensions are 1 × (*n* − 1) and *c* is an integer that defines the type of connection (see Table 1) [15]. The ICSA is developed using a probability criterion at each iteration, which will define the methodology to be used to define the new site of crow *i*. If the probability criterion is greater than the crow knowledge probability (i.e., *AP*) at iteration *k*, the traditional classical CSA search is employed [29]; otherwise, it will use the GD of the VSA to generate the new crow position *i*.

#### 3.2.1. Classical Approach: CSA

CSA is a metaheuristic optimization technique inspired by the intelligent performance of crows [29]. Crows are considered to be the smartest birds in nature; they possess a brain much larger in relation to the size of their bodies [36]. In groups, crows show notable traits of intelligence. They can learn and remember faces, use tools, communicate using sophisticated manners, and manage their food throughout the seasons because they hide their excess food in certain places (caches) in their environment and retrieve it when necessary [29].

Crows are ambitious birds as they pursue each other for better food reserves, observe where other birds hide their food, and steal it once they have left [37]. If a crow has stolen something, it takes additional cautions such as relocating its hiding places to prevent becoming a future victim [36]. They use their experience in theft to predict another robber's behavior and to determine the quickest course to protect their caches from being stolen [38]. The main bases of CSA are listed below [29]: (i) crows live in flocks, (ii) crows remember the site of their hiding places, (iii) crows pursue others to commit robberies, and (iv) crows shield their caches from being robbed using probability.

The following is an example of how the CSA mechanism works. Crows explore and exploit their environment, which is the solution space. Each environment cache corresponds to a feasible solution, the quality of the food source is the objective function. The best food source in the environment is the solution to the problem, Thus, the CSA seeks to simulate the intelligent behavior of crows to find the optimal solution [29].

The CSA is also based on a population. The population size (i.e., flock) consists of *N* individuals (number of crows), which belong to a *d* dimensional solution space, where *d* = (*n* − 1). The position *Xi*,*<sup>k</sup>* of crow *i* at iteration *k* is described in Equation (14) and represents a feasible problem solution.

$$X\_{i,k} = [\mathfrak{x}\_{i,k'}^1 \mathfrak{x}\_{i,k'}^2 \dots \mathfrak{x}\_{i,k}^n]\_\prime \tag{14}$$

where *i* = 1, 2, ..., *N* and *k* = 1, 2, ..., *t*max. Further, *t*max is the maximum number of iterations of the exploration and exploitation process of the solution space. Each crow (individual) can memorize the position of its hiding place. At iteration *k*, the position of crow *i* hiding place is represented as *Mi*,*<sup>k</sup>* , being the best position that crow *i* has obtained so far. Of course, the position of its best experience is memorized because crows move in their environment and search for the best food source (i.e., hiding places).

At the start of the optimization process, it is assumed that at iteration *k*, crow *j* wants to visit its hideout, *Mj*,*<sup>k</sup>* . In this iteration, crow *i* decided to follow crow *j* to approach crow *j*'s hideout. In this case, two situations can occur.

#### Situation 1: Search

Crow *j* does not know that crow *i* is following it. As a result, crow *i* approaches crow *j*'s hiding place. In this case, through Equation (15), the new position of crow *i* is obtained.

$$X\_{i,k+1} = X\_{i,k} + r\_i f l\_{i,k} \left( M\_{j,k} - X\_{i,k} \right),\tag{15}$$

where *r<sup>i</sup>* a random number between 0 and 1, and *f li*,*<sup>k</sup>* is the flight length of the crow *i* at each iteration *k*. According to [29], small values of *f l* allow a local exploration of the solution space (close to *Xi*,*<sup>k</sup>* ), while large values of *f l* allow a global solution space exploration (far from *Xi*,*<sup>k</sup>* ).

#### Situation 2: Evasion

Crow *j* knows that crow *i* is following it. As a result, crow *j* tries to deceive crow *i* by heading towards another position in the solution space to protect its caches from being stolen. Either way, situations 1 and 2 can be represented as shown in Equation (16).

$$X\_{i,k+1} = \begin{cases} X\_{i,k} + r\_i f l\_{i,k} (M\_{j,k} - X\_{i,k}) & \text{if } r\_j \ge A\_P \\ \text{random} & \text{otherwise} \end{cases} \tag{16}$$

where *r<sup>j</sup>* is a random number between 0 and 1, and *A<sup>P</sup>* denotes the probability of crow *j*'s knowledge of crow *i*.

Once the crows' positions are modified, the memory of each crow is updated based on the objective function values of the new spots. So, if the objective function of the new location is better than the objective function of the memorized position, the crows update their memory to the new area, as shown in Equation (17).

$$M\_{i,k+1} = \begin{cases} F(X\_{i,k+1}) & \text{if } F(X\_{i,k+1}) < F(M\_{i,k}) \\\ M\_{i,k} & \text{otherwise} \end{cases} \tag{17}$$

where *F*(·) represents the minimized objective function.

#### 3.2.2. Proposed ICSA

The CSA has confirmed its capability to reach the optimal solution for particular solution space configurations [29,39,40]. Nevertheless, its convergence is not ensured due to the inefficient exploration of its search strategy [37]. Therefore, it presents difficulties when facing high-dimensional problems [37,41]. In the original CSA method [29], there are two elements responsible for the search process: knowledge probability (i.e., *AP*) and random motion (i.e., Situation 2: Avoidance) [41].

The value of *A<sup>P</sup>* is entrusted with providing an adequate equilibrium between diversification and intensification [29]. With small *A<sup>P</sup>* values, a local solution space search is obtained, increasing intensification [29]. On the other hand, with large *A<sup>P</sup>* values, a global solution space search is obtained, which increases diversification [29]. Since metaheuristic algorithms require a balance between diversification and intensification to find a globally optimal solution when solving problems with large dimensions [42,43], *A<sup>P</sup>* is taken as an intermediate value, i.e., 50%.

Moreover, the random motion specifically impacts the CSA search mechanism since it resets the candidate solutions, deviating them from the current best solution and delaying the convergence of the problem [37]. In the proposed ICSA, the random motion is reformulated, as shown below.

#### 3.2.3. Improved Approach: VSA for Random Movement

The evasion behavior is simulated by a random motion implementation computed through a uniformly distributed random value [29]. In the proposed ICSA, to have a better solution space diversification, the possibility of generated candidate solutions based on the VSA evolution criteria is added to the classical CSA [15], considering that the main feature of the VSA is the use of a regular GD to generate neighbors around the current best solution named as the center of the hyper-sphere *µ* [44]. In this paper, the hyper-sphere center is selected as the best solution in the current population during iteration *k* as *µ<sup>k</sup>* = *X*best,*<sup>k</sup>* .

The set of candidate solutions *Cv*,*<sup>k</sup>* (*y*) is generated using a random GD in the space *d* around the best solution *Xbest*,*<sup>k</sup>* , as shown in Equation (18).

$$\mathcal{L}\_{\upsilon,k}(y) = p(y|\mu\_k, \Sigma) = \frac{1}{\sqrt{(2\pi)^d |\Sigma|}} \exp\left\{-\frac{1}{2} (y - \mu\_k)^T \Sigma^{-1} (y - \mu\_k) \right\},\tag{18}$$

where *<sup>d</sup>* is the solution space dimension, *<sup>y</sup>* <sup>∈</sup> <sup>R</sup>*d*×<sup>1</sup> corresponds to a vector of random variables with values between zero and one, and <sup>Σ</sup> <sup>∈</sup> <sup>R</sup>*d*×*<sup>d</sup>* is the matrix of covariance. If, in Σ, the on-diagonal elements (variance) are equally defined and if the off-diagonal components (covariance) are chosen as zero, then the GD will generate hyper-spheres in the *d*−dimensional space [35]. Equation (19) displays a simple way to calculate Σ, considering equal variance and zero covariances.

$$
\Sigma = \sigma^2 \mathbf{I}\_{d \times d \prime} \tag{19}
$$

where *<sup>I</sup>d*×*<sup>d</sup>* is an identity matrix and *σ* represents the variance of the GD. Note that the standard deviation of the GD can be defined as shown in Equation (20).

$$
\sigma\_0 = \frac{\max\{y^{\max}\} - \min\{y^{\min}\}}{2},
\tag{20}
$$

where *y* max and *y* min are vectors of dimension *<sup>d</sup>* <sup>×</sup> <sup>1</sup> that define the upper and lower bounds of the decision variables of the optimization problem, respectively. Here, *σ*<sup>0</sup> can also be considered as the initial radius *r*<sup>0</sup> of the hyper-sphere [35]. To achieve a proper exploration of the solution space, initially, *σ*<sup>0</sup> is the largest possible hyper-sphere.

The candidate solutions obtained and contained in the set *Cv*,*<sup>k</sup>* (*y*) must guarantee that the results lie within the bounds of the solution space; so, Equation (21) is employed.

$$\mathcal{C}\_{\upsilon,k}(y) = \begin{cases} \mathcal{C}\_{\upsilon,k}(y) & y^{\min} \le y \le y^{\max} \\ y^{\min} + (y^{\max} - y^{\min}) \text{rand} & otherwise \end{cases} \tag{21}$$

where rand generates random numbers between 0 and 1. Once verified the limits, the best solution obtained in the set *Cv*,*<sup>k</sup>* (*y*) is selected to be the new position of the crow *i*.

It is necessary to mention that the radius of the hyper-sphere decreases as the iteration process progresses using an inverse incomplete gamma function, as reported in [45] and as shown in Equation (22).

$$r\_k = \sigma\_0 \gamma^{-1}(y, a\_k) \tag{22}$$

The inverse incomplete gamma function for the variable radius calculation can be calculated in MATLAB®, as shown in Equation (23) [35].

$$r\_k = \sigma\_0 \frac{1}{y} \text{gammaincinc} \text{inv}(y, a\_k), \tag{23}$$

where *a<sup>k</sup>* is a parameter defined as *<sup>a</sup><sup>k</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>k</sup> t*max . Moreover, the parameter *y* is chosen as 0.1, as recommended in [35].

Algorithm 2 summarizes the implementation of the CSA, considering the VSA evolution criterion for Situation 2, to solve the phase swapping problem in ADSs.

**Remark 5.** *For the phase swapping problem solution, when using Algorithm 2, it is recommended to take y* max *as 6.5 and y* min *as 0.5. Further, d is the number of system demand nodes, except for the slack node.*

**Remark 6.** *The size of the solution set Cv*,*<sup>k</sup>* (*y*) *is chosen as 30% of the CSA population size to minimize the number of evaluations needed in the slave stage to determine the power losses.*

**Algorithm 2:** General ICSA implementation for solving the phase swapping problem in ADSs

Read information from the AC distribution system; Set the initial values *N*,*AP*, *f l*, and *t*max; Initialize the crows' position *Xi*,0 randomly; Calculate the objective function value for each crow *F*(*xi*,0) using Algorithm 1; Select the initial memory value *Mi*,0 for each crow *i*; Select the initial radius *r*<sup>0</sup> (or the standard deviation *σ*0) of (20); Set *k* = 1; **while** *k* ≤ *iter max* **do for** *i = 1: N* **do** Randomly select a crow *j* for tracking.; **if** *r<sup>j</sup>* ≥ *AP* **then** *Xi*,*k*+<sup>1</sup> = *Xi*,*<sup>k</sup>* + *r<sup>i</sup>* · *f li*,*<sup>k</sup>* · (*Mj*,*<sup>k</sup>* − *Xi*,*<sup>k</sup>* ) ; **else** Determine the center of the hyper-sphere *µ<sup>k</sup>* as *Xbest*,*<sup>k</sup>* ; Select the radius of the hyper-sphere *r<sup>k</sup>* ; Define the individuals' number *v* as 30% of *N*; Create the set of candidate solutions *Cv*,*<sup>k</sup>* (*y*) using (18); Check the lower and upper bounds for each *v* en *Cv*,*<sup>k</sup>* (*y*) using (21); Calculate the objective function value for each crow *v* in *Cv*,*<sup>k</sup>* (*y*) using Algorithm 1; Select *Xi*,*k*+<sup>1</sup> = as the individual with the best solution of *Cv*,*<sup>k</sup>* (*y*); **end end** Verify feasibility of the new positions *Xi*,*k*+<sup>1</sup> ; Evaluate the crows' new position *F*(*Xi*,*k*+<sup>1</sup> ); Update the crows' memory *Mi*,*k*+<sup>1</sup> ; Update the radius *rk*+<sup>1</sup> as shown in en (22); *k* = *k* + 1; **end Result:** Report the best solution *Xi*,*t*max , and its obj. func. value, i.e., *F*(*Xi*,*t*max ).

#### **4. Three-Phase Test Feeders**

We consider three test systems to validate the proposed ICSA. These test systems correspond to the 8, 25, and 37 node systems with radial topology reported in [15] for the phase swapping study using the VSA. Their main characteristics are presented below.

#### *4.1. 8-Bus Test Feeder*

The 8-node test system is a three-phase ADS formed by eight nodes and seven distribution lines, which operates with a nominal voltage of 11 kV at the main node. Figure 1 shows the electrical configuration of the test system. Note that the benchmark active power losses for this system take a value of 13.9925 kW. The electrical parameters for this test feeder can be found in [15].

**Figure 1.** Electrical design of the 8-node system.

#### *4.2. 25-Bus Test Feeder*

The 25-node test system is a three-phase ADS formed by twenty five nodes and twenty four distribution lines, which operates with a nominal voltage of 4.16 kV at the main node. Figure 2 displays the grid configuration. Note that the benchmark active power losses for this system take a value of 75.4207 kW. The electrical parameters for this test feeder can be found in [15].

**Figure 2.** Electrical design of the 25-node system.

#### *4.3. 37-Bus Test Feeder*

The 37-node system is part of a current ADS, consisting entirely of subway lines situated in California, USA. It has thirty seven nodes and thirty five distribution lines, and it operates at a nominal voltage of 4.8 kV at the substation. Figure 3 shows the grid electrical design of this test feeder. Note that the benchmark active power losses for this system take a value of 76.1357 kW. The electrical parameters for this test feeder can be found in [15].

**Figure 3.** Electrical design of the 37-node system.

#### **5. Numerical Simulations**

This section contains the numerical validation of the developed methodology to solve the phase swapping problem in 8-, 25-, and 37-node test systems, considering a given demand condition. In that sense, it uses the information provided by [15], which presents two methodologies to solve the proposed problem. Note that the aim is to compare the results obtained by the proposed procedure with each optimization technique reported in [15].

To solve the MINLP, formulated from (1) through (6), that represents the optimal phase swapping problem for ADSs, MATLAB® V2020*a* software is used on a laptop computer with an Intel(R) Core(TM) i5-7200U CPU @2.50 Ghz, a RAM of 8.00 GB, and a Windows 10 Home Single Language 64-bit operating system.

To test the efficiency of the proposed algorithm, the ICSA is compared with the Chu and Beasley genetic algorithm (CBGA) [15], the VSA [15], and the CSA. The parameters used for the CSA and ICSA are as provided in Table 2. Furthermore, the parameters were established with ten individuals in the population, six hundred iterations, and one hundred evaluations to calculate the average processing time. The parameters for the CSA are those recommended by the author of the algorithm in [29].



#### *5.1. Results for the 8-Node System*

Table 3 shows the results obtained by the proposed ICSA for the 8-node test system as follows: (i) all methodologies allow a reduction of more than 24% in total power losses; (ii) the solution obtained with the ICSA for the 8-node system equals those reported in [15], which states that the optimal global solution for this system is 10.5968 kW, albeit with a much longer processing time; and (iii) the standard deviation for the ICSA shown in Table <sup>3</sup> is 1.099 <sup>×</sup> <sup>10</sup>−<sup>5</sup> kW, which is at least ten times lower than in the VSA [15]. These results indicate that the repeatability of the solutions is close to 100% when solving the phase swapping problem in the 8-node test system, bearing in mind that the solution space has dimensions of 279,936.

**Table 3.** Performance of the power losses after implementing the phase-swapping plan in the 8-bus test feeder.


Figure 4 displays the variations of the phase power losses before and after the application of phase swapping with the ICSA, where the phase losses of *a* and *b* increase by 1.025 kW and 1.663 kW, respectively. On the contrary, the power losses in the electrical phase *c* decrease by about 6.10 kW, which increases the offsets seen in the power losses of phases *a* and *b*. Additionally, the power losses per phase are close to the average of the total power losses of approximately 3.50 kW, with differences of less than 0.80 kW. These results indicate that phase swapping by ICSA is an effective way to redistribute the loads in the phases of the system as evenly as possible.

**Figure 4.** Effect of phase swapping on power losses in the 8-bus test feeder.

#### *5.2. Results for the 25-Node System*

Table 4 shows the results obtained by the proposed ICSA for the 25-node test system, where the following is evident: (i) the solution obtained with the ICSA for the 25-node system equals the one reported in [15] for the VSA, which states that the optimal global solution for this system is 72.2888 kW; (ii) the solution obtained with the ICSA for the 25-node system outperforms the solution obtained with the CSA, which shows that the improvement made to the classical CSA explores and exploits the solution space for systems of large dimensions (i.e., 24 nodes in this case) in a better way; and (iii) the standard deviation for the ICSA, shown in Table 4, is 0.0116 kW lower than those reported in [15] and in CSA. This affirms the repeatability properties of the ICSA for solving the phase swapping problem, considering the size of the solution space, i.e., 2.8430 <sup>×</sup> <sup>10</sup><sup>19</sup> .


**Table 4.** Performance of the power losses after implementing the phase-swapping plan in the 25-bus test feeder.

Figure 5 represents the variations of the phase power losses before and after the application of phase swapping with the ICSA, where the phase losses of *b* increase by 11.3776 kW. On the other hand, the phase power losses of *a* and *c* approximately decrease by 11.2156 kW and 3.294 kW, respectively, which offsets the increase seen in the electrical phase *b* power losses. Additionally, the phase power losses are close to the average of the total power losses, approximately 24 kW, with differences of less than 3.60 kW. These results indicate that phase swapping by ICSA is an effective way to redistribute the loads in the phases of the system as evenly as possible.

**Figure 5.** Effect of phase swapping on power losses in the 25-bus test feeder.

Figure 6 shows a comparison of the percentage reductions of the total active power losses of the different methods displayed in Table 4, in contrast to the initial power losses. All methodologies allow a cutback of more than 4%. ICSA obtains a reduction of 4.15%, akin to the VSA reported in [15].

**Figure 6.** Reduction of the power losses in the 25-bus system.

#### *5.3. Results for the 37-Node System*

Table 5 shows the results obtained by the proposed ICSA for the 37-node test system, where the following is evident: (i) the solution obtained with the ICSA for the 37-node system improves the one reported in [15] for the VSA, which indicates that the optimal

solution for this system is 61.4781 kW; (ii) the solution obtained with the ICSA for the 37-node system outperforms the solution obtained with the CSA, which shows that the improvement made to the classical CSA explores and exploits the solution space for systems of large dimensions (i.e., 37-nodes) in a better way; and (iii) the standard deviation for the ICSA shown in Table 5 is 0.1344 kW, which is lower than those reported in [15] and the one obtained by the CSA. This affirms the repeatability properties of the ICSA for solving the phase swapping problem, considering the size of the solution space, i.e., 6.1887 <sup>×</sup> <sup>10</sup><sup>28</sup> .

**Table 5.** Performance of the power losses after implementing the phase-swapping plan in the 37-bus test feeder.


Figure 7 portrays the variations of the phase power losses before and after the application of phase swapping with the ICSA, where the phase *b* losses increase by 9.9184 kW. On the contrary, the phase power losses in *a* and *c* decrease by approximately 6.7771 kW and 17.7989 kW, respectively, offsetting the increase seen in the electrical phase *b* power losses. Additionally, the phase power losses are close to the average total power losses of approximately 20.50 kW, with differences of less than 1.40 kW. These results indicate that phase swapping by ICSA is an effective way to redistribute the loads in the phases of the system as evenly as possible.

**Figure 7.** Effect of phase swapping on power losses in the 37-bus test feeder.

Figure 8 compares the reductions in the total active power losses percentage of the different methods presented in Table 5 to the initial power losses. All methodologies allow a decrease of more than 19%, where ICSA obtains a reduction of 19.252%, a higher loss than VSA, as reported in [15], which reported a decrease of 19.249%.

In Figure 8, it is possible to observe that the proposed ICSA allows an additional improvement about 0.003% when compared to the power losses reduction with the VSA. This reduction implies a difference of 0.1784 kW in the total power losses minimization for the IEEE 37-bus system. Even if this power losses value corresponds to a small power losses reduction for this system, this demonstrates that the proposed ICSA finds an optimal solution for the IEEE 37-bus system, which supports the best current literature report obtained by [15] with the VSA method. Thus, this new solution will serve as a reference

point for future approaches that can be proposed to solve the phase swapping problem in ADSs.

**Figure 8.** Reduction of the power losses in the 37-bus test feeder.

#### *5.4. Complementary Results*

The following can be concluded from the results obtained in Section 5.


X The comparison of the effect in the energy losses redistribution at each phase using the classic and improved CSA methods is presented in Table 6.


**Table 6.** Comparison between the CSA and the proposed ICSA regarding power losses among phases.

With these results, we can make several observations. (i) In the 8-bus system, the phase *c* presents power losses higher than 4 kW, while the ICSA does not support this value; however, both solutions are indeed optimal since the total power losses is the same for both methods. (ii) In the 25-bus system both methods, i.e., the CSA and its improved version, the phase *c* has been identified to present higher power losses surpassing 26 kW; however, regarding the final power losses, the ICSA presents better load redistribution, since the amount of power losses is about 0.0408 kW. Finally, (iii) in the 37-bus system, the CSA method presents a difference between phases *a* and *b* power losses of 2.0450 kW, while the proposed ICSA has a minor difference between both phases with a value of 1.1067 kW. This directly impacts the final grid results with a general improvement of about 0.1784 kW in favor of the proposed ICSA, which is indeed the best global optimum reported in the current literature for the IEEE 37-bus system.

#### **6. Conclusions and Future Works**

This paper presented a master–slave methodology to solve the phase equilibrium problem in ADSs. For this purpose, an ICSA was proposed using the VSA evolution mechanism. The master stage determined the set of phase configurations for the system three-phase charges through the ICSA, using an integer encoding between 1 and 6 representing the load connections in the three-phase system. On the other hand, the slave stage evaluated each of the load connections in the TPPF, which correspond to the extended version of the iterative sweep power flow method for unbalanced three-phase systems.

The numerical results on the 8-, 25-, and 37-node test systems showed that the proposed ICSA, compared to the VSA, achieves identical solutions for the 8-node and 25-node systems. However, for the 37-node system, the ISCA obtains a better optimal solution when compared with the current report employing the VSA method. The solutions obtained by the ICSA are 10.5869 kW, 72.2888 kW, and 61.4781 kW for the 8-, 25-, and 37-node systems, representing a reduction of 24.34%, 4.152%, and 19.252%, respectively. In the same way, the solutions for the VSA are 10.5869 kW, 72.2888 kW, and 61.4801 kW, respectively.

In addition, the proposed methodology has minor standard deviations for solving the phase swapping problem for the 8-, 25-, and IEEE 37-node test systems, which were 1.099 <sup>×</sup> <sup>10</sup>−<sup>5</sup> kW, 0.0116 kW, and 0.1344 kW, respectively. This demonstrates better repeatability properties of the improved algorithm, since all the solutions are contained inside small hyper-spheres around the global optimum. Further, the solution space for the phase equilibrium problem potentially increases as a function of the demand nodes. Thus, for the

IEEE 37-node system, there is a solution space bigger than <sup>1</sup> <sup>×</sup> <sup>10</sup>28, which is by far higher than 100 billion combinations. This result confirmed the effectiveness of the proposed methodology for solving complex MINLP models such as the phase equilibrium problem in ADSs as well as in being better than the optimal solution reported in the current literature using the VSA.

In future work, it is possible to accomplish the following: (i) combine the proposed phase swapping with the optimal placement of distributed generators to reduce power losses in ADSs; (ii) employ typical active and reactive power demand curves to solve the phase swapping problem using the ICSA to reduce power losses; and (iii) use the ICSA to solve the optimal reconfiguration problem in three-phase radial distribution networks.

**Author Contributions:** Conceptualization, methodology, software, and writing—review and editing, B.C.-C., L.S.A.-G., O.D.M., L.A.-B. and C.Á.-A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported in part by the Laboratorio de Simulación Hardware-inthe-loop para Sistemas Ciberfísicos de la Universidad Loyola Andalucía.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**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.

#### **References**


## *Article* **Optimal Demand Reconfiguration in Three-Phase Distribution Grids Using an MI-Convex Model**

**Oscar Danilo Montoya 1,2 , Andres Arias-Londoño <sup>3</sup> , Luis Fernando Grisales-Noreña <sup>3</sup> , José Ángel Barrios 4,\* and Harold R. Chamorro <sup>5</sup>**


**Abstract:** The problem of the optimal load redistribution in electrical three-phase medium-voltage grids is addressed in this research from the point of view of mixed-integer convex optimization. The mathematical formulation of the load redistribution problem is developed in terminals of the distribution node by accumulating all active and reactive power loads per phase. These loads are used to propose an objective function in terms of minimization of the average unbalanced (asymmetry) grade of the network with respect to the ideal mean consumption per-phase. The objective function is defined as the *l*<sup>1</sup> -norm which is a convex function. As the constraints consider the binary nature of the decision variable, each node is conformed by a 3 × 3 matrix where each row and column have to sum 1, and two equations associated with the load redistribution at each phase for each of the network nodes. Numerical results demonstrate the efficiency of the proposed mixed-integer convex model to equilibrate the power consumption per phase in regards with the ideal value in three different test feeders, which are composed of 4, 15, and 37 buses, respectively.

**Keywords:** load redistribution; leveling power consumption per phase; three-phase asymmetric distribution networks; ideal power consumption; mixed-integer convex optimization

## **1. Introduction**

Most of the electricity users are typically connected to medium- and low-voltage levels, corresponding to three-phase distribution system structures [1]. The main characteristics of these networks are: (i) the radial connection among nodes helps to reduce the investment costs in protective schemes [2]; (ii) the existence of multiple single-, two-, and threephase loads produce current unbalances that increases the amount of power losses with respect to the perfectly balanced load scenario [3]; and the high grade of active and reactive power imbalances in terminals of the substation causes deterioration of voltage profile in the nodes located at the end of the feeder [4]. The importance of the threephase distribution networks to supply medium- and low-voltage users shows the need of proposing optimization methodologies to improve their electrical performance when the connection of new loads is required and the consumption at industrial nodes is increased [5]. The most common methodologies to improve the operative performance of the distribution networks are: optimal placement and sizing of reactive power compensators, i.e., capacitor banks and static distribution compensators [6,7]; optimal placement and sizing of disperse generation [8,9]; optimal grid reconfiguration [10]; and optimal phase-balancing [3,11].

**Citation:** Montoya, O.D.; Arias-Londoño, A.; Grisales-Noreña, L.F.; Barrios, J.Á.; Chamorro, H.R. Optimal Demand Reconfiguration in Three-Phase Distribution Grids Using an MI-Convex Model. *Symmetry* **2021**, *13*, 1–15. https://doi.org/ 10.3390/sym13071124

Academic Editors: Juan Alberto Rodríguez Velázquez and Alejandro Estrada-Moreno

Received: 30 May 2021 Accepted: 23 June 2021 Published: 24 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Note that the first two methodologies are based on the connection of new devices (shunt generators and compensators) to the network, which implies large amounts of investment to improve the quality of the grid, hardly recovered in short periods of time, i.e., 5 to 15 years [12–14]; the third methodology based on grid reconfiguration involves moderate investments in tie-lines and reconfiguration of protective devices [15]; whereas the phase-balancing method is the most simple strategy to reduce power losses in threephase distribution networks with minimum investment efforts, since devices are not required to implement the phase-balancing plan and it is only necessary to send few working crews along with the grid infrastructure to interchange the phase connections in the required nodes [3,16,17].

In the current literature can be found multiple optimization strategies, most of them based on evolutionary optimization algorithms to address the problem of the phasebalancing in three-phase networks. Some of these works apply to the following optimization methods: genetic algorithms [18,19]; tabu search algorithm [20]; particle swarm optimization [21]; ant colony optimization [22], and the vortex search algorithm [3], among others. The main characteristic of these evolutionary algorithms is the master–slave optimization strategy, where the master stage is entrusted of defining the connection of the loads using an integer or binary codification, while the slave stage determines the amount of energy losses at each connection provided by the master stage. Even if the master–slave optimization approach is widely accepted in the current literature, its main problem arises with not ensuring the global optimum finding. Recently, authors of [23] have proposed a mixed-integer quadratic programming model that allows to ensure the global optimum finding of the phase-balancing problem in three-phase networks; however, the effectiveness of this methodology was only tested in a small low-voltage microgrid. This fact has reduced the real impact of the convex methodologies in the general operation improvement of the electrical networks.

Based on the aforementioned arguments, this research deals with a problem similar to the phase-balancing approach, which is known as the load redistribution at terminals of the substation, i.e., the problem studied corresponds to leveling the active power demands per phase, making these to be accumulated in the main bus of the network (without considering the effect of the three-phase lines). The main advantage of this optimization problem is that it can be formulated with a mixed-integer convex (MIC) model that allows to ensure the global optimum finding by combining the Branch & and Bound method with the interior point method. The main contributions of this research are listed below:


The remainder of this document is structured as presented below:

Section 2 presents the proposed mixed-integer convex optimization model to represent the problem of the load redistribution in the terminals of the main substation; Section 3 presents the main characteristics of the optimization methodology based on the combination of the branch and bound (B&B) method with linear programming methods, ensuring the global optimum finding for MIC optimization models; Section 4 presents the main characteristics of the test feeders composed of 4, 15, and 37 nodes used to validate the proposed MIC optimization model using the CVX optimization tool with the MOSEK solver in the MATLAB programming environment; Section 5 describes the main numerical achievements in the three test feeders under study regarding the minimization of the general average grade of unbalance in terminals of the substation and the grid power losses. Finally, Section 6 presents the main concluding remarks derived from this work

based on the main numerical results obtained after solving the proposed MIC model with the MOSEK solver.

#### **2. Exact Mathematical Formulation**

In general terms, the problem of the load redistribution in electrical three-phase medium-voltage grids is a mixed-integer non-linear programming model due to the presence of the power balance equations [23]; however, here, we propose a mixed-integer convex (MIC) model that allows redistributing the loads at all the nodes by the accumulation in the main substation (i.e., by neglecting the effect of the electrical distribution grid) [11]. The main objective of the MIC is to find the global optimum by combining the (B&B) method with linear programming search methods. Figure 1 presents a schematic model of the load redistribution in a particular node of the network before and after the solution of the proposed MIC model.

**Figure 1.** Redistribution of the load in a particular node of the network after solving the proposed MIC model.

Note that the main objective of the proposed MIC model is to reduce the level of asymmetry within all the loads of the network, analyzed at terminal of the main substation, i.e., redistribute all the loads in the nodes to reach the maximum level of balanced (symmetry) in the power consumption.

The complete optimization model for load redistribution in electrical asymmetric networks is fully described below.

#### *2.1. Objective Function*

The objective function of the problem is the minimization of general grid unbalance of the network with respect to the ideal consumption per phase. The objective function proposed is defined by Equation (1).

$$\min \mathcal{U}\_{\%} = \left(\frac{100}{3P\_{\text{ave}}}\right) \sum\_{f \in \mathcal{F}} \left| P\_f - P\_{\text{ave}} \right| \tag{1}$$

where *U*% represents the average grid unbalance; *P*ave corresponds to the average active power consumption per phase, which is calculated as the total active power load divided by three; and *P<sup>f</sup>* is the total active power consumption per phase. Note that the F represents the set that contains all the phases of the network.

**Remark 1.** *The main advantage of the objective function defined in (1) comes from the fact that this corresponds to the l*1−*norm (i.e., absolute value) which is a convex function. The convexity property is important since it is possible to ensure the global optimum finding of the optimization problem if, and only if, the set of constraints are also convex, or by using the MIC through the combination of the (B&B) method with the Simplex method.*

#### *2.2. Set of Constraints*

The problem of the load redistribution is subject to linear constraints which correspond to the load reconfiguration at each phase, the binary nature of the decision variable, and the calculation of the total load per phase, among others. The complete list of constraints is defined as follows:

$$P\_{kf}^d = \sum\_{\mathcal{G} \in \mathcal{G}} \mathbf{x}\_{kf\mathcal{G}} P\_{k\mathcal{G}'}^d \ \{\forall k \in \mathcal{N}, \ \forall f \in \mathcal{F}\},\tag{2}$$

$$Q\_{kf}^d = \sum\_{\mathcal{G} \in \mathcal{G}} \mathbf{x}\_{kf\mathcal{g}} Q\_{\mathcal{kg}'}^d \ \{\forall k \in \mathcal{N}, \ \forall f \in \mathcal{F}\} \,\tag{3}$$

$$P\_f = \sum\_{k \in \mathcal{N}} P\_{kf'}^d \ \{ \forall f \in \mathcal{F} \}\_{\prime} \tag{4}$$

$$Q\_f = \sum\_{k \in \mathcal{N}} Q\_{kf'}^d \ \{ \forall f \in \mathcal{F} \}, \tag{5}$$

$$\sum\_{g \in \mathcal{F}} \mathfrak{x}\_{kf\mathcal{g}} = 1,\ \{\forall k \in \mathcal{N},\ \forall f \in \mathcal{F}\},\tag{6}$$

$$\sum\_{f \in \mathcal{F}} \mathfrak{x}\_{kf\mathcal{g}} = 1,\ \{\forall k \in \mathcal{N},\ \forall \mathcal{g} \in \mathcal{F}\},\tag{7}$$

where *P d k f* and *<sup>Q</sup><sup>d</sup> k f* are the active and reactive power connected at node *k* in phase *f* after the redistribution of the loads; *P d kg* and *<sup>Q</sup><sup>d</sup> kg* correspond to the active and reactive power connected at the node *k* in the phase *f* before the redistribution of the loads; *xk f g* is the binary variable that determines if the load connected in phase *g* is reassigned to the phase *f* in the node *k*; *P<sup>f</sup>* represents the total active power consumption of the network in the phase *f* after redistributing all the loads; *Q<sup>f</sup>* defines the total reactive power consumption of the network in the phase *f* after redistributing all the loads. Note that N represents the set that contains all the buses of the network.

The set of constraints defined from (2) to (7) are explained as follows: Equations (2) and (3) determine the amount of active and reactive power consumption at each phase and node after redistributing the loads in all the buses and phases of the network. Equations (4) and (5) determine the final equivalent active and reactive power consumption in the terminals of the substation after redistributing all the loads in order to minimize the average unbalance of the network; finally, Equations (6) and (7) ensure that each load is uniquely connected to one phase of the network.

**Remark 2.** *The general structure of the set of constraints above presented show that the optimization model defined from (1) to (7) is indeed MIC, which implies that its optimal solution is achievable with conventional mixed-integer optimization methods [24].*

**Remark 3.** *The solution of the optimization model presented in (1) to (7) with conventional optimization techniques such as the (B&B) and interior point methods ensures the global optimum finding based on the mixed-integer convex theory [25]; however, it is not possible that the combination of the variables that produce the optimum value is unique. This behavior is observed in the numerical analysis presented in the Results' section.*

To illustrate the effect of the three-dimensional characteristics of the decision variable *xk f g*, let us consider the possible three-phase load connections presented in Table 1.


**Table 1.** Possible load distributions in a three-phase network [3].

Note that the decision variables in Table 1 in the last column represent the six possible combinations for the connection of the three-phase loads [23], which clearly fulfills the requirements in equality constraints (6) and (7) associated with the uniqueness of the loads per phase.

It is worth mentioning that the exact formulation of the load redistribution problem in three-phase asymmetric networks is indeed a mixed-integer non-linear programming (MINLP) problem due to the power balance equations that relates the power injection at each node with the voltage and angle variables [18]; the main difficulty of the exact MINLP model lies in the non-convexity of the power balance equations that makes impossible to find the global optimum with exact or metaheuristic methods. To reduce this model complexity, here, we propose a reformulation of the load balancing problem in two stages which corresponds to the exact MIC model formulated from (1) to (7) in the first stage that helps with finding the optimal redistribution of the loads in all the nodes. The solution obtained in the first stage is evaluated at the second one to determine the final level of power losses of the grid. The two-level methodology proposed in this research is easily implemented at any optimization software that combines the (B&B) and interior point methods with the main advantage of ensuring the optimal solution as demonstrated in [25,26], some solvers that can deal with this type of problems are available in the GAMS and AMPL software. Numerical results that will be reported in the Results' section were corroborated with the CPLEX software in the GAMS software [27].

#### **3. Methodology of Solution**

To efficiently solve the MIC optimization model defined from (1) to (7) it is possible to use any programming language that deals with convex optimization. Here, we adopt the CVX optimization package in the MATLAB programming environment with the MOSEK solver. The main characteristic of the optimization model is that the objective function is defined as the *l*1-norm which is a convex function. The illustration of the objective function in a three-dimensional space is presented in Figure 2.

**Figure 2.** Representation of the objective function in a three-dimensional space *z* = |*x*1| + |*x*2|.

It is important to mention, that as with the most of the integer optimization models, an MIC can be solved with a modification of the (B&B) method as presented in Figure 3 [28]. Note that at each iteration, it is solved a linear programming model which ensures the optimal solution finding at each nodal exploration [29].

**Figure 3.** General application of the B&B method for addressing MIC problems.

**Remark 4.** *Notice that the MIC model defined from (1) to (7) can be rewritten as a mixed-integer quadratic programming problem which is also convex, i.e., it is also possible to find its global optimum by combining the interior point and the (B&B) method [25].*

To verify the efficiency of the optimization model to balance the total power consumption in the substation terminals and its positive effects on the minimization of the power losses, we evaluate the final load reconfiguration in an asymmetrical three-phase power flow method to determine the final power losses and compare with the initial state of the network. The power flow methodology used in this research corresponds to the matrix version of the backward–forward method reported in [3,30].

#### **4. Electric Distribution Grids**

The computational validation of the proposed MIC to redistribute loads in three-phase networks considering a simplified model in the substation terminals is made with three different test feeders composed of 4, 15, and 37 buses, respectively. The information of these test feeders is presented below.

#### *4.1. 4-Bus Test Feeder*

The 4-bus system is a medium-voltage grid with 4 nodes and 3 lines with a nominal line-to-line voltage of 11.4 kV. The information of the loads and branches are listed in Tables 2 and 3, respectively. This information was obtained from [3].

**Table 2.** Parametric information of the 4-bus test system (kW and kvar units are used for all powers).


**Table 3.** Impedances' information for the conductors used in the 4-bus system.


#### *4.2. 15-Bus Test Feeder*

This test feeder is composed by 15 buses and 14 branches, asymmetric three-phase nature network with 13.2 kV of nominal phase voltage at the node of the substation, which corresponds to the typical operating voltage in Colombian power distribution grids. In Figure 4, it is shown the electrical configuration of this test feeder.

**Figure 4.** Nodal connections in the 15-bus test feeder.

Tables 3 and 4 present the parametric data of the 15-bus system.


**Table 4.** Parametric information of the 15-bus test system (kW and kvar units are used for all power values).

#### *4.3. IEEE 37-Bus Test Feeder*

The IEEE 37-bus system is a three-phase unbalanced network that is a portion of a real power grid located in California, USA. This grid has 37 nodes with radial connection among them. The line-to-line voltage assigned to the substation bus is 4.8 kV. Note that the electrical configuration of this test feeder was taken from [18] where some variations to the grid topology were included. The single-phase equivalent diagram of the IEEE 37-bus system is presented in Figure 5.

**Figure 5.** Nodal connection of the IEEE 37-bus system.

The complete parametric information for this test feeder is reported in Tables 5 and 6. Note that this information was taken from [3].


**Table 5.** Parametric information of the IEEE 37-bus test system (kW and kvar units are used for all power values).

**Table 6.** Data of impedance for the conductors used in the IEEE 37-bus system.


method [30].

*5.1. 4-Bus System*

**5. Computational Validation**

#### **5. Computational Validation** The computational validation of the proposed MIC programming model is made in the MATLAB environment with the CVX optimization package and the MOSEK solver.

The computational validation of the proposed MIC programming model is made in the MATLAB environment with the CVX optimization package and the MOSEK solver. In addition, we evaluate the power losses before and after solving the MIC model using the matrix version of the backward–forward asymmetrical three-phase power flow method [30]. In addition, we evaluate the power losses before and after solving the MIC model using the matrix version of the backward/forward asymmetrical three-phase power flow

#### *5.1. 4-Bus System* This test feeder presents an initial power losses of 68.6292 kW with an average grid

*Symmetry* **2021**, *1* 10

This test feeder presents an initial power losses of 68.6292 kW with an average grid unbalance of 22.47%. After solving the load reconfiguration problem, total grid power losses is 62.5449 kW, i.e., a reduction with respect to the benchmark case is about 8.87%; in addition, the general grid unbalance is reduced until 0.74%. In Figure 6 is presented the comparison between the initial and the final grade of unbalance per phase. unbalance of 22.47%. After solving the load reconfiguration problem, total grid power losses is 62.5449 kW, i.e., a reduction with respect to the benchmark case is about 8.87%; in addition, the general grid unbalance is reduced until 0.74%. In Figure 6 is presented the comparison between the initial and the final grade of unbalance per phase.

**Figure 6.** Initial and final unbalances in the 4-bus system. **Figure 6.** Initial and final unbalances in the 4-bus system.

Note that all the phases are effectively balanced with differences lower than 1.20 % with respect to the ideal consumption case, i.e., *P*ave; in addition, the final active power load at phases *a*, *b*, and *c* in terminals of the substation are 1220 kW, 1200 kW, and 1200 kW, respectively; where the variations with respect to the initial case were 30 kW, 370 kW, and Note that all the phases are effectively balanced with differences lower than 1.20% with respect to the ideal consumption case, i.e., *P*ave; in addition, the final active power load at phases *a*, *b*, and *c* in terminals of the substation are 1220 kW, 1200 kW, and 1200 kW, respectively; where the variations with respect to the initial case were 30 kW, 370 kW, and 400 kW, for phases *a*, *b*, and *c*, respectively.

400 kW, for phases *a*, *b*, and *c*, respectively. Table 7 presents the final solution regarding load connections after solving the MIC model to redistribute all the loads. The most important result observed in this table corresponds to the existence at least of two possible solutions for the MIC model in the 4-bus system. This happens in this test feeder since two of the phases ends with 1200 kW of total load, which implies that some rotations in the phase connections will exhibit the Table 7 presents the final solution regarding load connections after solving the MIC model to redistribute all the loads. The most important result observed in this table corresponds to the existence at least of two possible solutions for the MIC model in the 4-bus system. This happens in this test feeder since two of the phases ends with 1200 kW of total load, which implies that some rotations in the phase connections will exhibit the same final active power losses.

same final active power losses. **Table 7.** Optimal solutions reached by the MIC optimization model in the 4-bus system.


Solution 1 {1, 6, 4, 5} 62.7868 8.51 0.74 Solution 2 {1, 2, 1, 3} 62.5449 8.87 0.74 It is important to mention that the solutions obtained with the MOSEK solver in It is important to mention that the solutions obtained with the MOSEK solver in the CVX environment were validated with the GAMS optimization package with the COUENNE solver. In addition, the average processing time in MATLAB including the power flow evaluations was about 1.83 s.

In this test feeder, previous to the application of the MIC model to redistribute all the loads among the phases of the network, we know that the initial power losses is 134.2472 kW, caused by a general unbalance of 20.48 %, which is distributed as 2.68 % for

the CVX environment were validated with the GAMS optimization package with the COUENNE solver. In addition, the average processing time in MATLAB including the

power flow evaluations was about 1.83 s.

*5.2. 15-Bus System*

#### *5.2. 15-Bus System*

In this test feeder, previous to the application of the MIC model to redistribute all the loads among the phases of the network, we know that the initial power losses is 134.2472 kW, caused by a general unbalance of 20.48%, which is distributed as 2.68% for phase *a*, 30.72% pf phase *b*, and 28.04% for phase *c*. Once the MIC model defined from (1) to (7) is executed, we observe that the active power losses is uniformly distributed for all the phases with respect to the average value. In Table 8 is presented the load redistribution in the 15-bus system.

**Table 8.** Comparison between initial and final load distributions per phase (all the values in kW and kVar).


Note that results in Table 8 show that: (i) there are at least two solutions of the optimization model (1)–(7) that present the same objective function performance which in this system is exactly zero; this implies that all the phases have the same active power load consumption per phase, i.e., 9354 kW; (ii) the general unbalance in the case of reactive power for this system in the benchmark case is 26.01% which is reduced to 8.42% after making the load redistribution; this result implies an important effect when redistributing the total consumption per phase in the substation terminals, since the modification of the active power load connection is directly connected with the total reactive power consumption; (iii) the final power losses for solutions 1 and 2 are 117.8982 kW and 115.1107 kW, with reductions respect to the benchmark case of about 12.18%, and 14.25%, respectively; and (iv) note that the main difference between solutions 1 and 2 corresponds to the rotation of the loads connected between phases *a* and *c* in all the nodes; this results important since 4 additional solutions can be obtained making possible the 6 load rotations presented in Table 1 with the same objective function of 0.00%, and power losses between 117.8982 kW and 115.1107 kW, respectively.

In regards with the total processing time the MOSEK solver using the MATLAB/CVX environment takes about 174.45 s; which is a quite small processing time taking into account that there are 6 *<sup>n</sup>*−<sup>1</sup> possible combinations of the loads, where *n* is the total number of nodes, i.e., 78,364,164,096, this is, more than 78,000 million of combinations.

#### *5.3. IEEE 37-Bus System*

The initial average unbalance in the IEEE 37-bus system is 22.14%, which is distributed in 11.23% for the phase *a*, 21.98% for the phase *b*, and 33.21 for the phase *c*, respectively. Once it is solved the proposed MIC model the final average unbalance in the network is 1.71%. Figure 7 depicts the initial and final unbalances per phase.

Numerical results per phase in Figure 7 show that phases *a*, *b*, and *c* are improved in about 9.16%, 19.41%, and 32.72%, respectively, with respect to the benchmark case of active power; which confirms the efficiency of the proposed approach for optimizing the general average grid unbalance, guaranteeing the global minimum of the problem.

that there are 6

*5.3. IEEE 37-Bus System*

and 115.1107 kW, respectively.

in the 15-bus system.

kVar).

phase *a*, 30.72 % pf phase *b*, and 28.04 % for phase *c*. Once the MIC model defined from (1) to (7) is executed, we observe that the active power losses is uniformly distributed for all the phases with respect to the average value. In Table 8 is presented the load redistribution

**Table 8.** Comparison between initial and final load distributions per phase (all the values in kW and

Note that results in Table 8 show that: (i) there are at least two solutions of the optimization model (1)–(7) that present the same objective function performance which in this system is exactly zero; this implies that all the phases have the same active power load consumption per phase, i.e., 9354 kW; (ii) the general unbalance in the case of reactive power for this system in the benchmark case is 26.01% which is reduced to 8.42% after making the load redistribution; this result implies an important effect when redistributing the total consumption per phase in the substation terminals, since the modification of the active power load connection is directly connected with the total reactive power consumption; (iii) the final power losses for solutions 1 and 2 are 117.8982kW and 115.1107kW, with reductions respect to the benchmark case of about 12.18%, and 14.25%, respectively; and (iv) note that the main difference between solutions 1 and 2 corresponds to the rotation of the loads connected between phases *a* and *c* in all the nodes; this results important since 4 additional solutions can be obtained making possible the 6 load rotations presented in Table 1 with the same objective function of 0.00%, and power losses between 117.8982 kW

In regards with the total processing time the MOSEK solver using the MATLAB/CVX environment takes about 174.45 s; which is a quite small processing time taking into account

The initial average unbalance in the IEEE 37-bus system is 22.14 %, which is distributed

i.e., 78364164096, this is, more than 78,000 million of combinations.

1.71 %. Figure 7 depicts the initial and final unbalances per phase.

*<sup>n</sup>*−<sup>1</sup> possible combinations of the loads, where *n* is the total number of nodes,

**Scenario** *P<sup>a</sup> Q<sup>a</sup> P<sup>b</sup> Q<sup>b</sup> P<sup>c</sup> Q<sup>c</sup> U*% (%) Benchmark case 9605 5226 6480 4940 11977 8778 28.04 Solution 1 9354 7112 9354 5794 9354 6038 0.00 Solution 2 9354 6038 9354 5794 9354 7112 0.00

**Figure 7.** Initial and final unbalances in the IEEE 37-bus system. **Figure 7.** Initial and final unbalances in the IEEE 37-bus system.

Numerical results per phase in Figure 7 show that phases *a*, *b* and *c* are improved in about 9.16 %, 19.41 %, and 32.72 % respectively with respect to the benchmark case of In relation with the amount of power losses, the benchmark case presents the initial power losses of 76.1357 kW; however, for this test system after solving the MIC model there are six possible combinations that produce different levels of power losses reduction. Table 9 reports all the possible solutions in regards with power losses obtained by our proposed optimization approach.


**Table 9.** Optimal solutions reached by the MIC optimization model in the 15-bus system.

Results in Table 9 allow concluding that: (i) the first solution obtained by the MIC approach presents the best numerical performance regarding grid power losses with a reduction of 12.55% in comparison to the benchmark case; (ii) the worst solution regarding of power losses corresponds to solution 4, which is obtained by rotating solution 1 from XYZ to XZY in all the nodes since the reduction of the power losses in this case decreases until 11.37%; (iii) all the solutions in Table 9 are indeed the global optimum for the optimization model (1)–(7) since the general grid imbalance is 1.71% for all the solution cases; however, the calculation of the final power losses can be considered a decision criterion to select the most attractive solution from the point of view of the grid operator, which, in this context, is solution 1.

Finally, with respect to total processing time, the MOSEK solver using the MAT-LAB/CVX environment takes about 29,879.44 s; which is an acceptable processing time taking into account that there are 6 *<sup>n</sup>*−<sup>1</sup> possible combination of the loads, where *n* is the total number of nodes, i.e., 1.03144247984905 <sup>×</sup> <sup>10</sup><sup>28</sup> .

#### *5.4. Additional Comments*

It is worth mentioning that all the numerical results reported with the CVX tool in the MATLAB environment for the 4-, 15-, and 37-bus systems were confirmed by the CPLEX solver in the GAMS environment with simulation times that do not get over 10 s [27]. These processing times confirm the effectiveness of the MIC model to solve the problem of the load redistribution in asymmetrical three-phase networks by ensuring the global optimum finding in the first stage of the two-level proposed optimization method [26]. The solutions provided in the first stage were evaluated in the triangular-based three-phase power flow which takes less than 10 ms to solve it and determine the final level of power losses in the grid at the second stage [3].

On the other hand, the proposed two-stage optimization methodology to redistribute the load connections in three-phase networks is suitable to be applied in the improvement of the resilience level of the electricity distribution activity [31]. This is in the context of the physical- or cyber-attacks to the distribution network or any electrical disturbance that implies the reconfiguration of the grid topology; since the proposed optimization model can be applied to each possible grid topology to ensure that after clarifying the disturbance the resulting electrical network has the minimum level of unbalance, in other words, minimum power losses.

#### **6. Conclusions**

The problem of the load redistribution in three-phase distribution networks was addressed in this research from the point of view of the mathematical optimization, by proposing a mixed-integer convex model that ensures the global optimum finding via (B&B) and linear programming methods. Numerical results in the 4-, 15-, and 37-bus systems demonstrate the effectiveness of the proposed optimization model to reduce the general average grid unbalance of the network by reducing from 22.47% to 0.74% for the 4-bus system, 20.48% to 0.00% for the 15-bus system; and 22.14% to 1.71% for the IEEE-37 bus system, respectively.

In addition, for each test feeder it was observed that when the final load reconfiguration is rotated for all the possible phase combinations, the amount of power losses in the final configuration changes with the same objective function value, which confirms the multi-modal behavior of the load redistribution optimization problem in terminals of the main substation. For the 4-, 15-, and IEEE 37-bus systems the maximum power losses reductions with respect to the benchmark case were 8.87%, 14.25%, and 12.55%, respectively. These reductions demonstrate the strong relationship between the load redistribution problem and the total grid power losses, which can be taken as an advantage of the grid owner to improve the quality of the electricity service at the same time that increases the net profit due to the reduction in the costs of the energy losses.

Regarding the processing times it was observed that for each test feeder the MOSEK solver in the CVX package using the MATLAB environment takes 1.73 s, 174.45 s, and 29,879.44 s, for the 4-, 15-, and IEEE 37-bus systems, respectively, which can be considered as short processing times due to the large dimension of the solution space. This latter can be calculated as 6 *n*−1 , being 216, 7.8364164096 <sup>×</sup> <sup>10</sup>10, and 1.031442479849054 <sup>×</sup> <sup>10</sup><sup>28</sup> for the test feeders mentioned previously. As future work it will be possible to analyze the inclusion in the proposed MIC model of the power balance constraints with a second-order cone representation that will ensure the global optimum finding for the problem of the phase-balancing problem in three-phase networks.

**Author Contributions:** Conceptualization, O.D.M., A.A.-L., L.F.G.-N., H.R.C. and J.A.B.; Methodology, O.D.M., A.A.-L., L.F.G.-N., H.R.C. and J.A.B.; Investigation, O.D.M., A.A.-L., L.F.G.-N., H.R.C.; and writing—review and editing, O.D.M., A.A.-L., L.F.G.-N., H.R.C. and J.A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** 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".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

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

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

#### **References**

