Application of the Multiverse Optimization Method to Solve the Optimal Power Flow Problem in Alternating Current Networks

: In this paper, we solve the optimal power ﬂow problem in alternating current networks to reduce power losses. For that purpose, we propose a master–slave methodology that combines the multiverse optimization algorithm (master stage) and the power ﬂow method for alternating current networks based on successive approximation (slave stage). The master stage determines the level of active power to be injected by each distributed generator in the network, and the slave stage evaluates the impact of the proposed solution on each distributed generator in terms of the objective function and the constraints. For the simulations, we used the 10-, 33-, and 69-node radial test systems and the 10-node mesh test system with three levels of distributed generation penetration: 20%, 40%, and 60% of the power provided by the slack generator in a scenario without DGs. In order to validate the robustness and convergence of the proposed optimization algorithm, we compared it with four other optimization methods that have been reported in the specialized literature to solve the problem addressed here: Particle Swarm Optimization, the Continuous Genetic Algorithm, the Black Hole Optimization algorithm, and the Ant Lion Optimization algorithm. The results obtained demonstrate that the proposed master–slave methodology can ﬁnd the best solution (in terms of power loss reduction, repeatability, and technical conditions) for networks of any size while offering excellent performance in terms of computation time.


Introduction
To understand the importance of this study, it is necessary to analyze the general context of distributed generation and its current importance for electrical networks, as well as the optimal power flow problem in AC grids. This section reviews the main contributions of other authors, presents the contributions of this study, and outlines the structure of this paper.

General Context
Electrical energy has become the backbone of society as we know it today because, when used properly and wisely, it contributes to improving the quality of life of people all around the world [1][2][3]. However, the comfort provided by electrical energy has also led to an increase in its consumption. As a result, power stations have been forced to meet the demand by producing electricity from fossil fuels and renewable energy sources and to implement large-scale electric power systems, which experience significant power losses and unfavorable operating conditions, such as voltage drops and current limit violations [4][5][6][7]. Thus, the ways in which electrical energy is generated and distributed should be expanded while keeping in mind critical aspects, e.g., reduction of greenhouse gas emissions, minimization of power losses associated with energy transport and distribution, and the provision of a reliable and high-quality service to end users. To satisfy these needs, both researchers and network operators have studied and proposed different power generation and distribution technologies and promoted the development and implementation of distributed energy resources such as distributed generators, batteries, and static reactive compensators [8][9][10]. In particular, Distributed Generators (GDs) are the distributed energy resource most widely studied and installed around the world; therefore, they are the focus of this paper.
DGs can be integrated into both Alternating (AC) and Direct Current (DC) networks. Nevertheless, the integration of these devices into AC networks has been more widely studied and applied because AC grids are more common around the world than their DC counterparts [10,11]. An adequate operation of DGs in an electrical system could improve its chargeability, reduce the power losses, and enhance the voltage profiles of the grid, among other aspects [12,13]. Considering these advantages and the widespread use of AC power distribution systems, this study focuses on optimal power flow applied to DGs in AC networks.
To ensure the correct operation of DGs in AC networks, first, we should identify the main components of these networks, which include the slack generator, distributed generators, energy distribution branches, and electrical loads [10]. Second, we need to establish the location of the DGs in the network and their technical constraints: number of generators and nodes, minimum and maximum distributed power levels, branch data, and voltage and current bounds, among others. Finally, we should devise operating strategies to control the power injected by the distributed energy resources into the network to obtain as many economic, technical-operating, and environmental benefits as possible [14]. This problem is referred to as the Optimal Power Flow (OPF) problem in the specialized literature. The OPF is a nonlinear nonconvex problem that must be solved using optimization techniques and numerical methods to find the power configuration of DGs that contributes to the improvement of technical, economical, and environmental indicators [15]. To solve this problem, a master-slave methodology is usually implemented. In this kind of methodology, the master stage is responsible for determining the power to be injected by each distributed energy resource in order to satisfy the objective function defined by the network operator. The slave stage is in charge of solving the load flow problem in order to establish the impact of each solution proposed in the master stage [16]. Multiple authors in the specialized literature have proposed different optimization techniques for the master and slave stages in order to improve the impact of DGs in an electrical grid, as can be observed in the following subsection.

State of the Art
In recent years, several authors interested in the energy sector have looked for and proposed solutions to the various problems that occur in AC networks, including power losses associated with energy distribution [17,18], optimal sizing and location of DGs in the network [19,20], proper energy storage using batteries [15], and reduction of CO 2 emissions [7]. Regarding the optimal power flow in AC networks that include DGs, different solution methodologies have been proposed in the literature. Some of them have used commercial software and solutions based on sequential programming implementing open-source software. These two types of solution methodologies help to find the optimal power levels to be injected by the DGs in a network.
Studies such as [21,22] employed commercial optimization software to solve the OPF problem, satisfying the limits and constraints associated with this problem. In [21], the authors proposed a discrete-continuous solution methodology, in which the Chu-Beasley genetic algorithm was implemented in the DigSILENT programming language. In that study, the tests were only conducted in 6-, 14-, and 39-node test systems, thus ignoring the performance of the proposed methodology in larger systems. Additionally, although the computation times and the voltage profiles were taken into account, the currents flowing through the lines and the standard deviation of the implemented methodology were not considered. As a result, the maximum current levels of the branches that compose the electrical system cannot be evaluated, which allows the method to propose infeasible solutions to the problem under analysis and eliminates the possibility of evaluating the repeatability of the solution each time the algorithm is executed. In [22], the authors presented a multi-objective optimization methodology based on the mixed-integer nonlinear programming model. The problem was modeled in GAMS language and later solved using the DICOPT solver. In that study, the OPF problem was tested in 14-and 30-node systems, thus ignoring its behavior in larger networks. However, it took into account computation times that the study (the same as the previous one) did not consider, the standard deviation of the solution methodologies, the currents flowing through the lines, or the nodal voltages.
Studies such as [23][24][25] employed algorithms based on sequential programming to solve the OPF problem in AC networks, thus avoiding the need for commercial software, which is costly and complex [12]. In [23], the authors proposed a biogeography-based optimization algorithm and used, as the objective function, the minimization of P loss in the IEEE 30-and 57-node test systems. In that study, both the computation times and the nodal voltages were taken into account; however, the standard deviation was not considered to evaluate the repeatability of the solution provided by the proposed methodology. In [24], an artificial bee colony optimization algorithm was employed, the objective function was the minimization of active power losses in AC networks, and the solution was tested in the IEEE 30-and 118-node systems. In the study mentioned, the authors considered the analysis of the nodal voltage constraints but not the computation times or the standard deviation of the techniques they used. In [25], a particle swarm optimization algorithm was proposed to solve the OPF problem in AC networks and tested using the IEEE 57-and 118-node systems. That study took into account the nodal voltages but not the computation times or the standard deviation of the implemented algorithms. Importantly, none of these studies considered the currents flowing through the lines, which does not ensure that the solution obtained with the implemented methodologies satisfies the technical-operating constraints of the test systems used in each paper. In [26], the authors proposed an efficient analytical (EA) method for solving the OPF in an electrical distribution system in order to minimize power losses. Their optimization method was tested in the 33-and 69-node systems, obtaining excellent results. They used multiple methods to make comparisons, but they did not consider the analysis of the processing times or the repeatability of the solutions obtained by the EA. In [27], the authors presented a methodology implementing the Biogeography-Based Optimization (BBO) method to solve the OPF problem in a power transmission system; their objective function was the reduction of power losses. All their simulations were carried out in Matlab using the 114-node test system. However, in that paper, they did not compare the their results with those reported in other studies in the literature. The authors of [28] proposed a master-slave methodology that combines the black hole (BH) optimization algorithm and the triangular-based power flow method to solve the OPF problem in an AC network; their objective function was the reduction of power losses. The numerical results they obtained in the 33-and the 69-node test systems demonstrated the effectiveness and robustness of the proposed approach compared to those of others in the literature. However, they did not analyze the repeatability of the algorithm or the processing times it required; furthermore, they neglected the current limits of the branches in their mathematical model. In [29], the authors used a particle swarm optimization algorithm to solve the optimal power flow problem by sizing a photovoltaic DG, wind DG, and a battery inside an electrical network. Their objective function was the reduction of the investment and management costs of these devices using a mono-nodal electrical network as the test system. Nevertheless, they did not compare their method to others and did not analyze the processing time it required or its repeatability. In [30], the authors addressed the OPF problem in radial distribution networks using a metaheuristic optimization technique known as the sine-cosine algorithm; their objective function was the reduction of power losses. They validated the effectiveness of the proposed solution in the 33-and 69-node test systems by comparing its results with those reported in other papers in the literature. Nevertheless, they did not consider the effects of the solution methods under analysis on processing times. Furthermore, they did not include the current limits of the branches in the constraints that compose the mathematical model.
It is important to highlight that none of the optimization algorithms employed in the studies above were tuned, making it impossible to guarantee the same conditions for the proposed solution techniques. Furthermore, all these methodologies proposed to solve the OPF were evaluated in AC networks with radial topologies, without analyzing their effect in mesh topologies, which are very common in conventional electrical systems.
After reviewing the state of the art, we identified the need for a new methodology based on sequential programming to solve the OPF problem in AC networks that also avoids the use of commercial software and can be replicated using open-source software. The methodology should consider different topologies and sizes of AC networks to obtain high quality solutions in short processing times. It should also analyze the operating constraints of AC network components that were not considered in the state of the art reviewed above, such as nodal voltage, current flowing through the transmission lines in the systems, and power generated by the DGs. To evaluate the repeatability of the solution method, the results should be assessed using a statistical approach that analyzes the standard deviation of each methodology. In order to satisfy this need, this article presents a master-slave methodology to solve the OPF problem that can be replicated using open-source software. The master stage uses the Multiverse Optimization (MVO) algorithm proposed by [31]. It was selected due to its high performance in solving research problems and the fact that it has been applied to issues associated with renewable energies, power generation, and distribution systems [32][33][34][35]. The slave stage uses the Successive Approximation (SA) method proposed in [36] because of its high performance in terms of convergence and processing time to solve the load flow problem, as reported by the authors. To validate this methodology, the 10-, 33-, and 69-node radial systems and a 10-node mesh system were used in this study. Each test system considers a maximum power injection allowed for the DGs of 20%, 40%, and 60% of the power provided by the slack generator in a scenario without distributed generation. Each scenario was tested 100 times to determine the repeatability of the proposed solution methodology.

Scope and Main Contributions
Considering the importance of minimizing power losses associated with energy transport in AC networks, as well as the need to develop new methodologies that help to find reliable and good quality solutions without using commercial software, this study proposes a master-slave methodology to solve the OPF problem in radial and mesh electrical topologies. In the master stage, the proposed methodology uses the MVO [31] to determine the power to be injected by each DG located in the network. In the slave stage, it applies the SA numerical method [36] to evaluate both the impact of the solution provided by the master stage and the set of constraints that compose the problem using the load flow. The following are the main contributions of this study to the state of the art: • A new solution methodology based on sequential programming to solve the OPF problem in AC networks that considers radial and mesh topologies.

•
A robust statistical method that can be used to identify the solution methodology with the best balance of solution quality, processing times, and repeatability.
• An optimization method selected because it yields the best results in terms of solution quality, repeatability, and processing time to solve the OPF problem in AC networks with radial and mesh topologies.
This study also makes the following contributions to the energy sector: • A new methodology that can be used to determine the optimal power level to be injected by each DG in order to reduce power losses and thus better use existing energy resources.

•
A methodology with short processing times to solve the problem of sizing DGs in electrical networks.

Structure of the Paper
This article is organized as follows: Section 2 presents the mathematical formulation of the OPF problem, as well as the constraints that compose it, where the objective function is the reduction of power losses associated with energy distribution. Section 3 describes the proposed solution methodology, which consists of a master stage that uses the MVO and a slave stage that employs the SA numerical method. Section 4 presents and describes the algorithms used here for comparison purposes. Section 5 details the radial and mesh test scenarios employed in this study to perform the simulations. Section 6 reports and analyzes the simulation results obtained by each optimization algorithm. Finally, Section 7 draws the conclusions and proposes future research directions.

Mathematical Formulation
A variety of objective functions can be used in AC networks: the minimization of the operating costs of the network, the optimal selection of power conductors in distribution networks, the reduction of CO 2 emissions, the minimization of power losses in the network, the optimal sizing and location of the DGs, and the improvement of voltage profiles. The selection of the objective function depends on the needs of the network operator. Due to its widespread use to evaluate the efficiency of the methodologies proposed to solve the OPF problem in AC networks [21][22][23][24][25], the objective function selected in this study is the minimization of the active power losses associated with energy distribution in AC networks [25].

Objective Function
The objective function is the variable to be optimized while respecting the constraints associated with the problem. The objective function selected for this study is represented by Equation (1) below: In this equation, P loss denotes the active power losses (i.e., the variable to be minimized in this study); v, a vector containing the complex voltages of the system calculated based on the load flow; Y L , the admittance matrix associated with the distribution lines; and v T , a transposed vector of v.

Set of Constraints
The set of constraints refers to the technical-operating limits of the equipment that composes the AC network. These limits, which are detailed below, are directly linked to the objective function selected in this study: Equations (2)-(6) represent the OPF problem. Equation (2) expresses the total (active and reactive) power balance in the network, where S CG is the complex power provided by the slack generator; S DG , the complex power supplied by the DGs; S D , the complex power demand; D(v), a symmetrical matrix that contains the complex voltages of the system in its diagonal; and Y L , a matrix of admittances associated with the nodal effects. Equation (3) presents the apparent power limits allowed for the DGs, where S min DG and S max DG are the minimum and maximum power that each DG can supply to the network, respectively. Importantly, the DGs only inject active power into the network. Equation (4) defines the nodal voltage limits of the system, where v min and v max are the minimum and maximum voltages allowed in each node of the system. Equation (5) represents the maximum current allowed inside the power distribution lines, where I B denotes the current flowing through the transmission lines; and I max B , the maximum current that can pass through the conductor. Equation (6) defines the maximum penetration percentage allowed for the DGs in the network, where 1 T is a row vector filled with ones; and α, the penetration percentage allowed for each DG, which can be 20%, 40%, or 60% of the power provided by the slack generator. Finally, in this equation, Real{X} corresponds to the real part of the complex vector (X) associated with the injection of active power by the DGs: In addition to the equations that represent the set of constraints, we also present Equation (7), which guarantees that the constraints (Equations (2)- (6)) are respected at all times. This equation penalizes the algorithm if the limits presented above are violated. In this equation, β 1 to β 6 are the penalty factors, which, for this study, were heuristically calculated and are equal to 1000. When all the constraints are respected, all the penalty values are canceled out using the min{.} and max{.} functions, which causes z to be equal to P loss .
To facilitate the comprehension of this article, Table 1 presents the notation and nomenclature used for each variable in Equations (1)-(7), as well as the acronyms of the OPF problem.

Proposed Solution Methodology
The equations presented in the previous section represent the OPF problem taking into account the operation of the DGs located in the AC network. Since this is a nonlinear nonconvex problem, it must be solved using optimization techniques and numerical methods. This study, therefore, proposes dividing the OPF problem into two stages: master and slave. In the first stage (master), the MVO [31] is used to determine the active power to be injected by each DG installed in the network. In the second stage (slave), the SA numerical method [36] is employed to solve the power flow problem and calculate the electrical variables provided by the master stage in order to determine their impact on the objective function we selected. The master-slave (MVO-SA) methodology proposed in this paper is described in the next subsections.

Master Stage: Multiverse Optimizer (MVO)
Although there are various theories about the origin of the universe, the Big Bang theory is one of the most well-known [37,38]. According to it, the Big Bang was an explosion that created the known universe [39,40]. Physicists have also proposed several theories suggesting the existence of universes other than the one we live in. One of them is known as the Multiverse theory, which claims that there are other universes and that each began as a result of an explosion (or Big Bang) [41][42][43][44]. The MVO is inspired by this theory, and it employs three of its main concepts: white holes, black holes, and wormholes [31].
Physicists regard the Big Bang as a white hole and the main component for the start of a new universe. White holes create and expel matter and energy [45][46][47][48]. Black holes, whose behavior is the complete opposite, absorb everything that comes near their event horizon, and nothing can escape them due to their colossal gravitational force [49][50][51][52][53]. Finally, wormholes connect two distant points in a universe, allowing objects to pass through them and creating instantaneous travel paths in the vastness of the universe [54][55][56][57].
The MVO mathematically models the behavior of white holes, black holes, and wormholes, and its search process consists of two stages. The first stage is exploration, in which white holes and black holes are used. In this stage, the algorithm explores the search space and suggests the most promising regions where the best local optima can be found. The second stage is exploitation, in which wormholes are employed. In this stage, a search is carried out in the regions where the best local optima have been found in order to find the global optimum.
The following subsections present the steps followed in this study for the computational development and to find a solution to the OPF problem using the MVO.

Generating the Initial Population
Equation (8) is used to generate the initial population of the universes, where each universe constitutes a solution to the problem addressed here (i.e., the power level to be injected by the DGs). In this equation, subscript i denotes the i-th universe, and subscript j represents the j-th object in universe i. In the OPF problem, each row represents a possible solution to the problem; and each column, the number of nodes (other than the slack node) that are allowed to inject power. Each universe (U i,j ) created with Equation (8) is located within the solution space, which is limited by the technical constraints of the OPF problem. In this equation, ub and lb are the maximum and minimum power that each DG can inject, respectively. In order for the algorithm to explore larger regions of the search space, the initial population of U i,j is generated randomly using the rand command, which produces random values between 0 and 1 for each object j within universe i: To generate the entire population of universes, this study proposes a population contained in a matrix of size nxd, where n is the number of universes, and d denotes the number of variables in the problem addressed here (see Equation (9)). In Equation (9), U n is the n-th universe in the M Universes matrix (which, in the case of the OPF problem, represents the possible solutions), and U d denotes the number of DGs located in the network:

Calculating the Objective Function
To assess the impact of every solution proposed in the M Universes matrix, the objective function of each U n is evaluated using the SA numerical method. This method evaluates the power levels injected by the DGs, which makes it possible to determine the power losses associated with each universe (U n ). The obtained values are stored in a matrix of size n × 1 called MO Universes , which is presented in Equation (10): Additionally, the MO Universes matrix is sorted in ascending order, where the first position represents the objective function with the best value, i.e., the lowest P loss for the system. This is selected as the incumbent solution to the OPF problem and stored in Best_Fob (Equation (11)). The variables that compose the best objective function are also stored in Best_U (Equation (12)), and this universe is selected as the one with the best solution found at the current iteration:

Movement Strategy
The MVO technique is divided into two stages: exploration and exploitation. In the exploration stage, this technique uses black and white holes to explore the search space, and the algorithm tries to discover the most promising regions where the global optimum can be potentially found. In the exploitation stage, this technique uses wormholes to exploit the promising regions found in the exploration stage, and the algorithm converges toward the global optimum (moves forward toward the best solution). In the MVO method, each solution is assigned an inflation rate (IR i ) that corresponds to the best P loss obtained so far by each universe. Such inflation rate is proportional to the fitness function provided by the slave stage after evaluating the objects j proposed by each one of the universes (i.e., the power levels to be injected by the DGs). In addition, the following rules are applied to the universes of the MVO: • The higher the IR, the greater the probability of having a white hole.

•
The lower the IR, the greater the probability of having a black hole.

•
Universes with a high IR tend to send objects through white holes.

•
Universes with a low IR tend to receive more objects through black holes. • Any object in any universe can randomly move toward the best universe through wormholes, regardless of the IR.
Furthermore, black or white holes must change their behavior depending on the required solution approach (i.e., a maximization or minimization problem). In this case, it is a minimization problem.
In the MVO method, objects (search particles) move between universes through black/white holes, i.e., objects can be exchanged between universes. When a relationship is established between white/black holes in two universes, the universe with the higher IR is considered to have a white hole, whereas that with the lower IR is considered to have black holes. Thus, information is transferred from the white holes in the source universe to the black holes in the target universe. Moreover, regardless of their IR, all universes are assumed to have wormholes in order to randomly transport objects. Wormhole existence probability (WEP) and travel distance rate (TDR), which are denoted by Equations (13) and (14), respectively, are two fundamental coefficients for the optimal operation of the MVO. Equation (13) defines WEP as the probability of existence of wormholes in the universes, which increases linearly over iterations to emphasize exploitation as progress made in the optimization process: This equation is able to generate a rate of change based on the sum of a minimum value (Min) and a maximum value (Max), both of which depend on the algorithm and the iterative process. In this equation, l is the current iteration, and L denotes the maximum number of iterations. This allows for changes at the exploration and exploitation levels of the optimization technique. Importantly, the Min and Max values were optimized to improve the quality of the solution to the problem addressed here.
Additionally, in Equation (14), TDR defines the distance over which an object can be transported through a wormhole around the best universe obtained thus far. Unlike WEP, TDR increases over iterations to provide a more precise local exploitation/search around the best universe found thus far: In this equation, P denotes the level of precision of the exploitation of the algorithm, i.e., the higher its value, the faster and more precise the search exploitation of the algorithm.
In this study, parameter P was optimized to find the best solution to the OPF problem in AC networks.

Evolution of the Universes in the Iterative Process
The relationship between white and black holes allows objects to be exchanged between universes through white/black tunnels. In addition, a universe can update the value of an object based on the information of an object from another universe. Their interaction occurs as a function of IR. Objects travel from a source white hole to a target black hole. Black holes make objects in a universe disappear, and they are replaced with others coming from white holes. Equation (15) shows how the normalized inflation rate of the i-th universe is obtained, where the IR values calculated for each universe in the population are normalized. When this equation is applied to the set of universes, the vector of size n × 1 (where n is the number of solutions or universes) is obtained with normalized values in the [0-1] range: The white/black tunnel is generated by means of Equation (16), where each object j in universe i can be replaced with an object j from universe k, and the position of k is selected using a roulette wheel that offers the option of arbitrarily selecting any universe in the solution space. The exchange of objects between universes occurs for the same element, and this is achieved by comparing a random number r1 with the normalized inflation rate (NI(IR i )).
In the equation above, U k,j is the object j of universe k selected by the roulette wheel and transported through the white/black tunnel to universe i. An object is transferred to universe i if and only if the random value r1 in the [0-1] range is below the N I(IR i ). Otherwise, there is no exchange of objects through white/black tunnels, i.e., the object of universe U i,j will not be updated and will keep its current value. Equation (17) is used to determine the k-th universe to which object j will be transported. In this equation, k is chosen using the roulette wheel function [58], which selects a random integer in the [1-n] range, where n is the number of solutions proposed for the problem. Algorithm 1 presents the pseudocode to select universe k by means of Equation (17): Importantly, since this is a minimization problem, the selection using the roulette wheel must be performed with −N I. If the problem to be solved is a maximization one, −N I should be replaced with the positive N I. Algorithm 1 Proposed pseudocode of the roulette wheel to select universe k to transport element j to universe i.

Updating Universes Based on Wormholes
Wormholes allow for the instantaneous teleportation of objects j to the best universe found thus far, where IR represents the evaluation of the objective function (in this case, the minimization of power losses in the system). As the algorithm advances, each universe experiences random teleportation of its objects to the best universe (Best_U(1, j)) through wormholes. This process is repeated until a stopping criterion is met.
To guarantee changes in each universe, there must be a high probability of improvement of the IR through wormholes, assuming that the latter are always established between new universes randomly generated depending on parameter r4 and the best universe (Best U (1, j)), which are represented in Equation (18): This equation can be used to generate new objects j for universe i (U i,j ) at each iteration, which depends on WEP and random parameter r2 in the [0-1] range. If r2 < WEP, a new element j will be generated for universe i depending on the value that parameter r3 takes, which is a random number in the [0-1] range. As selection criteria, we have that, if r3 < 0.5, Equation (18) (a) will be used, but, if r3 ≥ 0.5, Equation (18) (b) will be employed. However, if r2 ≥ WEP, a new object j will not be generated in universe i; instead, the current value of object j will be kept until r2 < WEP.

Stopping Criteria
The following two stopping criteria are employed for the master stage: • Maximum number of iterations: The iterative process will finish when the algorithm reaches the maximum number of iterations (L), which are controlled by counter l.

•
Number of non-improvement iterations: The algorithm will stop when the incumbent solution is not updated after n consecutive iterations.

Slave Stage
The slave stage allows us to determine the impact of each solution proposed in the master stage by evaluating the objective function associated with each solution (U n,d ). In other words, the slave stage assesses the electrical variables provided by the master stage, which are used to estimate the P loss in the system. Hence, a numerical method is necessary to solve the load flow problem in an iterative manner and ensure that the constraints of the problem are respected. Considering that, this paper proposes implementing the SA numerical method presented in [59], which was chosen due to its excellent performance in solving the power flow problem and because it is able to solve the load flow problem for radial as well as mesh networks. This method is based on Equation (19): In this equation, Y dd and Y dg are the components of the nodal admittances matrix of the network that relate load nodes to generator nodes, respectively. The variables v d and v g represent the voltages in the demand and generator nodes, respectively. The nodal voltage profiles in the demand nodes can be calculated by means of a mathematical development applied to Equation (19): In order to estimate the nodal voltages (v d ) of the system with an almost-null convergence error, an iterative process must be implemented. For this purpose, a counter (t) is included in Equation (20). As a result, the following equation can be used to calculate the nodal voltages of the system: Figure 1 shows a flowchart that describes the iterative process of the master-slave (MVO-SA) methodology used here to solve the OPF problem.
START Read the data of the electrical systems and the MVO algorithm.
Generate initial universe population.
Size the DGs fixed for each individual that composes the population by evaluating their objective function and technical constraints.
Select the universe with the best objective function as incumbent.

Methods Used for Comparison Purposes
To evaluate the capacity of the MVO to converge toward the best solution in terms of quality and computation times, it was compared with four other optimization techniques that have been reported in the specialized literature to solve the OPF problem: the Particle Swarm Optimization (PSO) algorithm, the Ant Lion Optimization (ALO) algorithm, the Black Hole (BH) optimization algorithm, and a Continuous Genetic Algorithm (CGA). The PSO algorithm is inspired by the behavior of schools of fish when exploring a region in search of food [60]. The ALO mimics the hunting mechanism of ant lions, which dig a cone-shaped hole in the ground to catch their prey [61]. The BH algorithm is based on the dynamic interaction between stars and black holes [62]. Finally, the CGA is inspired by the genetic crossover, selection, and mutation processes of living beings, through which parents pass on their DNA to their offspring [63]. These four optimization algorithms were selected because of their excellent performance in solving the OPF problem in AC networks [19,21,30,64]. Each one of them was implemented here in the master stage of the master-slave methodology. Thus, the master stage employed one of these optimization algorithms, and the slave stage always used the SA numerical method to solve the load flow problem [36].
In this study, we used the 10-, 33-, and 69-node radial test systems and the 10-node mesh test system (which are described in Section 5) to solve the OPF problem and assess the effectiveness of the proposed methodology in these types of electrical systems.
To ensure a fair comparison of all the optimization algorithms we employed, each one of them was tuned to guarantee the same conditions and that they could find the best solution to the problem addressed in this paper. To carry out such tuning, we used the PSO algorithm [60] with an initial population of 10 individuals and a maximum number of iterations of 300. The tuning parameters were a population size in the [1-100] range, a maximum number of iterations for each algorithm in the [1-1000] range, and a number of non-improvement iterations in the [1-1000] range. Moreover, for the MVO, parameters P, Max, and Min were tuned using [1][2][3][4][5][6][7][8][9][10], [0.001-1], and [0.001-1] ranges, respectively. Table 2 presents the turned parameters that allow each optimization algorithm to find the best solution to the OPF problem. Importantly, the parameters reported in Table 2 can be used to replicate each optimization technique in networks of any size.

Test Scenarios and Considerations
To evaluate the robustness and quality of the solution provided by each method, we used the 10-, 33-, and 69-node radial test systems, but we also employed a 10-node mesh test system to assess the effectiveness of the proposed methodology in mesh topologies. These four test systems were selected in order to assess the ability of each optimization algorithm to find the best solution in networks of different sizes. In each test system, DGs were allowed to inject 20%, 40%, and 60% of the power provided by the slack generator. In addition, each system features a single slack generator in a scenario without DGs (i.e., the base case).

Radial Test Systems
This subsection presents the radial test systems used in this paper. Figure 2 illustrates the 10-node system, which consists of 9 lines and 10 nodes. In a scenario without distributed generation, the slack generator produces 12,591.4181 kW, and the active power losses amount to 223.4181 kW. In this system, the DGs were located at nodes 5, 9, and 10, considering three maximum penetration levels: 20%, 40%, and 60% of the active power provided by the slack generator in the base case. For all the DGs, the minimum injection power was 0 kW at each penetration level. Likewise, the maximum injection powers for the 20%, 40%, and 60% penetration levels were 2518.2836 kW, 5036.5673 kW, and 7554.8509 kW, respectively. The maximum current for this system (i.e., 581.2757 A) was calculated using the load flow. Considering this value, we used a 1250-kcmil conductor operating at 75 • C, which establishes a maximum line current of 590 A. Since this study investigates a non-telescopic electrical network, this current limit applies to all the lines in the system. For all the calculations related to this system, we employed a base power of 100 kVA and a base voltage of 23 kV.  Table 3 reports the electrical parameters of the 10-node system. The first two columns specify the output and input nodes. The third column shows the resistance of the conducting lines. The fourth column details the reactances. Finally, the fifth and sixth columns present the active and reactive powers, respectively, which correspond to the powers demanded by each node in the system.  Figure 3 presents the 33-node system, which consists of 33 nodes and 32 lines. The base case of this system uses a base voltage of 12.66 kV and a base apparent power of 100 kVA. In the base case, the slack generator injects an active power of 3925.9785 kW, and the active power losses amount to 210.9785 kW. The locations of the DGs in this system were defined as in [10], i.e., they were located at nodes 12, 15, and 31. As in the previous system, DGs are allowed to inject 20%, 40%, and 60% of the power supplied by the slack generator in the base case. For all the DGs, the minimum active power was 0 kW, and the maximum active powers were 785.1957 kW, 1570.3914 kW, and 2355.5871 kW for the 20%, 40%, and 60% penetration levels, respectively. When the base case was evaluated, the maximum line current was 365.2518 A. Based on this value, we used a 700-kcmil conductor operating at 60 • C in each segment of the system, which establishes a maximum line current of 385 A.  Table 4, which is organized in the same way as Table 3, reports the electrical parameters of the 33-node system.  Figure 4 shows the 69-node system, which consists of 69 nodes and 68 lines. The scenario without DGs in this system employs a base voltage of 12.66 kV and a base apparent power of 100 kVA. In addition, the slack generator supplies an active power of 4132.8423 kW, and the active power losses amount to 242.1523 kW. As in the other two systems, DGs are allowed to inject 20%, 40%, and 60% of the power provided by the slack generator. The location of the DGs in this system was defined as in [10], i.e., they were located at nodes 26, 61, and 66. For all the DGs, the minimum active power was 0 kW, and the maximum active powers were 826.5685 kW, 1653.1369 kW, and 2479.7054 kW for the 20%, 40%, and 60% penetration levels, respectively. When the load flow was simulated for this system, the maximum current was 394.4489 A. Therefore, a 750-kcmil conductor operating at 60 • C was used in each segment of the system, with a maximum current of 400 A.   Table 5, which is organized in the same way as Tables 3 and 4, presents the parameters of the 69-node system.

Mesh Test System
This subsection describes the mesh test system used in this study. Figure 5 presents the topology of the 10-node mesh network. This 10-node system is the same we used previously, but, in this case, it has 10 nodes and 11 lines. The same voltage and apparent power of the 10-node radial test system are used in the base case. In the case without distributed generation, the slack generator provides a power of 12,558.3237 kW, and the active power losses equal 190.3237 kW. As in its radial counterpart, the DGs are located at nodes 5, 9, and 10, and this feeder uses penetration percentages of 20%, 40%, and 60% of the power produced by the slack generator. This way, each DG can produce a minimum power of 0 kW and maximum powers of 2511.6647 kW, 5023.3295 kW, and 7534.9942 kW, respectively. The maximum voltage flowing through the lines is 579.7276 A. Thus, a 1250-kcmil conductor at 75 • C was selected according to NTC 2050 standard, and a maximum line power of 590 A was established for each section of the network. Using the same order as in Tables 3-10 presents the electrical parameters of this system.     Table 9. Simulation results in the 69-node radial test system.   Table 10. Simulation results in the 10-node mesh test system.

Simulations and Results
This section discusses the results obtained by each optimization algorithm used to solve the OPF problem in the four test systems mentioned above. All the simulations were performed in Matlab R (version 2021b), a numerical computing system, running on a laptop with 4 GB of RAM, an Intel R Core TM i5-8250U @1.60 GHz 1.80 GHz processor, a 225-GB solid-state drive, and Windows 10 PRO. To guarantee the same conditions and evaluate the repeatability and precision of each solution technique, the codes were executed 100 times. The following subsections present the results obtained in each radial and mesh test system.

Results in the Radial Test Systems
This subsection presents and analyzes the results obtained in each radial system used as test scenario. Table 7 reports the results obtained by each optimization algorithm after performing the simulations in the 10-node radial system at three levels of distributed generation penetration: 20%, 40%, and 60% of the power supplied by the slack generator. From left to right, this table details the optimization algorithm implemented; the node where the DGs are located and the power that each one of them injects into the system in kW; the minimum power losses (P loss ) in kW and the percentage of reduction with respect to the base case (%); the average P loss in kW and the percentage of reduction with respect to the base case (%); the computation time required by each technique to obtain the solution (s); the standard deviation of each optimization technique (%); the worst voltage profile (p.u) and the node where it occurs; and the maximum current in the electrical system (A). Additionally, the P loss in the base case (system without DGs), i.e., 223.4181 kW, are shown in the upper part of the table. Table 7 compares the results obtained by the five optimization algorithms employed to solve the OPF problem in the 10-node system. Figures 6 and 7, which were constructed based on such results, illustrate the differences between these algorithms in terms of reduction of minimum and average power losses. Figure 6 compares the reduction of minimum P loss obtained by each algorithm at the three penetration levels: 20%, 40%, and 60%. Considering that the MVO is the base scenario, this figure presents the reduction percentage obtained by the MVO in comparison to the other optimization algorithms. The number in red indicates that the proposed MVO was outperformed by the PSO. At 20% penetration, the MVO achieved a reduction of 47.6667%, which was outperformed by the PSO by 0.0001% but was better than those of the ALO, CGA, and BH by 0.0113%, 0.0535%, and 0.4486%, respectively. At 40% penetration, the MVO and PSO obtained a reduction percentage of 63.8522%, outperforming the CGA by 0.0089%, ALO by 0.0140%, and BH by 0.0965%. At 60% penetration, the MVO and PSO exhibited the same reduction, i.e., 67.7170%, thus outperforming ALO, the CGA, and BH by 0.0021%, 0.0038%, and 0.0107%, respectively.  Figure 7 compares the reduction of average P loss obtained by the optimization techniques at the three penetration levels: 20%, 40%, and 60% of the power provided by the slack generator in the base case. As in Figure 6, this figure presents the MVO as the base scenario. At 20% penetration, the MVO achieved the best reduction, i.e., 47.6654%, thus outperforming PSO by 0.1284%, the CGA by 0.2484%, ALO by 0.4448%, and BH by 2.0591%. At 40% penetration, the MVO obtained the best reduction of average P loss , outperforming PSO, the CGA, ALO, and BH by 0.0969%, 0.1100%, 0.4888%, and 0.7498%, respectively. At 60% penetration, the MVO and PSO exhibited the best reduction of average P loss , i.e., 67.7170%, thus outperforming ALO (in the third place) by 0.0021%, the CGA (fourth place) by 0.0038%, and BH (fifth place) by 0.0107%. To complete the analysis of the 10-node system, Figure 8, which was created based on the results presented in Table 7, compares the standard deviation obtained by the optimization techniques at the three penetration levels: 20%, 40%, and 60%. At 20% penetration, the MVO presented the best standard deviation, i.e., 0.005%, thus outperforming the CGA by 0.1685%, ALO by 0.7162%, PSO by 1.3230%, and BH by 1.7415%. At 40% penetration, the MVO exhibited the best standard deviation, outperforming the other optimization algorithms by an average percentage of 0.9915%. At 60% penetration, the MVO and PSO obtained standard deviations of 9.4 × 10 −7 and 1.22 × 10 −8 , respectively; thus, the MVO outperformed the CGA by 0.0610%, BH by 1.1291%, and ALO by 1.6134%. This demonstrates that, when evaluated in small networks, the MVO produces excellent results in terms of the repeatability of the solution. According to these results, the MVO has excellent precision and repeatability, which guarantees that a good quality solution is obtained every time it is executed. Moreover, it can be concluded that the MVO is the most appropriate technique to solve the OPF problem in small AC networks, as it ensures that a good quality solution is found in a relatively short time in any distributed generation scenario. Table 8, which is organized in the same way as Table 7, reports the results of the simulations performed in the 33-node system after the algorithms were executed 100 times to solve the OPF problem. Importantly, the P loss in the base case were 210.9785 kW, as shown in Table 8. The same as with the 10-node system above, Figures 9-11 compare the minimum P loss , average P loss , and standard deviation, respectively, obtained by each algorithm at the three penetration levels (i.e., 20%, 40%, and 60%). Figure 9 shows the reduction of minimum P loss obtained by each optimization method in the 33-node radial system at the three penetration levels. At 20% penetration, the MVO and PSO achieved the same reduction of minimum P loss with respect to the base case, i.e., 39.5681%, thus outperforming ALO by 0.0021%, the CGA by 0.0098%, and BH by 0.0603%. At 40% penetration, the MVO and PSO exhibited the same reduction percentage, i.e., 57.1629%, thus outperforming ALO by 0.0043%, the CGA by 0.0118%, and BH by 0.582%. At 60% penetration, the MVO and PSO obtained the same reduction percentage once again, i.e., 59.3423%, thus outperforming ALO, the CGA, and BH by 0.0007, 0.0011%, and 0.0121%, respectively. Considering that the MVO is the base scenario, Figure 10 reports the differences in the reductions of average P loss obtained by the MVO and the other optimization algorithms at the three penetration levels. At 20% penetration, the MVO presented a reduction of average P loss of 39.5676%, outperforming the CGA (in the second place) by 0.0496%, ALO (third place) by 0.0605%, PSO (fourth place) by 0.1857%, and BH (fifth place) by 0.4508%. At 40% penetration, the MVO exhibited the best reduction of average P loss , i.e., 57.1626%, thus outperforming the other algorithms by an average percentage of 0.2443%. At 60% penetration, both the MVO and PSO achieved a reduction of 59.3423%, outperforming the CGA by 0.0100%, ALO by 0.1094%, and BH by 0.2806%.

33-Node Radial Test System
To complete the analysis of the 33-node system, Figure 11 illustrates the standard deviation obtained by each technique used to solve the OPF problem in AC networks. As can be seen in this figure, the MVO is superior to the other optimization algorithms in each distributed generation scenario. At 20% penetration, the MVO presented a standard deviation of 0.001%, outperforming the CGA by 0.0430%, ALO by 0.0901%, BH by 0.4034%, and PSO by 0.5231%. At 40% penetration, the MVO exhibited the highest precision, with a standard deviation of 0.001%, outperforming the other optimization algorithms by an average percentage of 0.5510%. At 60% penetration, the MVO obtained a standard deviation of 6.1 × 10 −7 %, thus outperforming PSO by 1 × 10 −5 %, the CGA by 0.0168%, ALO by 0.3471%, and BH by 0.6068%. After a general analysis of these results, the MVO proves to be superior to the other algorithms in terms of standard deviation and precision, which guarantees that a good quality solution is obtained in a short computation time every time this algorithm is executed. Therefore, it can be concluded that the MVO is the technique that offers the best solution, in a short computation time, for medium-sized networks.  Table 9, which is organized in the same way as Tables 7 and 8, presents the results obtained by the optimization algorithms used to solve the OPF problem in the 69-node radial test system.
The results reported in Table 9 were used to construct Figures 12-14, which compare the minimum P loss , average P loss , and standard deviation, respectively, obtained by each optimization algorithm implemented here to solve the OPF problem in the 69-node radial test system. The same distributed generation scenarios considered in the 10-and 33-node radial networks above were employed for this system. Regarding minimum P loss , Figure 12 shows the reduction obtained by the MVO in comparison to the other optimization algorithms. At 20% penetration, the MVO presented a reduction of minimum P loss of 44.8433% with respect to the base case, which was only outperformed by the PSO by 0.0003% but was better than those of ALO, the CGA, and BH by 0.0289%, 0.0533%, and 0.1584%, respectively. At 40% penetration, the MVO and PSO exhibited a reduction in minimum P loss of 64.2963%, thus outperforming the CGA by 0.0040%, ALO by 0.0101%, and BH by 0.2166%. At 60% penetration, the MVO and PSO achieved a reduction of 68.2193%, outperforming ALO by 0.0006%, the CGA by 0.0006%, and BH by 0.0168%. Figure 12. Percentage of reduction of minimum power losses obtained by the MVO in the 69-node radial system compared to other methodologies. Figure 13 illustrates the differences in average P loss obtained by the MVO and the other optimization techniques. At 20% penetration, the MVO achieved the best reduction of average P loss , i.e., 44.8410%, thus outperforming the other optimization algorithms by an average percentage of 0.6703%. At 40% penetration, the MVO exhibited the best reduction of average P loss , i.e., 64.2958%, thus outperforming the CGA, PSO, ALO, and BH by 0.0587%, 0.0788%, 0.2508%, and 1.6601%, respectively. Finally, at 60% penetration, the MVO and PSO presented the same reduction of average P loss , i.e., 68.2193%, thus outperforming the CGA by 0.0116%, ALO by 0.1788%, and BH by 0.8731%.
To complete the analysis of the 69-node radial system, Figure 14 shows the standard deviation obtained by each optimization technique used to solve the OPF problem in in this network at the three penetration percentages. At 20% penetration, the MVO presented a standard deviation of 0.003%, thus outperforming the other optimization algorithms by an average percentage of 0.9329%. At 40% penetration, the MVO exhibited a standard deviation of 0.002%, thus outperforming the CGA, ALO, PSO, and BH by 0.0958%, 0.6241%, 1.6621%, and 1.9223%, respectively. Finally, at 60% penetration, the MVO and PSO obtained standard deviations of 1.3 × 10 −6 and 1.5 × 10 −6 , respectively; therefore, the MVO outperformed the CGA by 0.0237%, ALO by 0.7409%, and BH by 1.8238%.
These results demonstrate once again that the MVO is superior to the other techniques in terms of standard deviation, making it the most precise and appropriate method to solve the OPF problem in large AC networks because it guarantees that good quality solutions are obtained every time it is executed.

Results in the Mesh Test System
This section presents and analyzes the results obtained in the mesh system used as a test scenario. Table 10 presents the results obtained by the optimization algorithms to solve the OPF problem in the 10-node mesh system after they were executed 100 times. The information appears in the same order as in Tables 3-5. Figures 15-17 resulted from analyzing Table 10. These figures compare the results obtained by the optimization algorithms in terms of minimum P loss , average P loss , and standard deviation, respectively. Figure 15 compares the results achieved by the optimization algorithms to obtain the minimum P loss in the 10-node mesh system with 20%, 40%, and 60% penetration. In the figure mentioned, the y-axis denotes the improvement, in terms of minimum power loss reduction, obtained by the MVO compared to the other optimization algorithms; in turn, the x-axis specifies the methods used for comparison. At 20% penetration, the MVO has the best reduction in minimum P loss , outperforming PSO, ALO, the CGA, and BHO by 0.0001%, 0.0250%, 0.0296%, and 0.1150%, respectively. At 40% penetration, the MVO presents a reduction in minimum P loss of 69.2705%, outperforming PSO by 0.0002%, ALO by 0.0068%, the CGA by 0.0179, and BH by 0.0758%. Finally, at 60% penetration, the MVO presents the best reduction in minimum P loss again, i.e., 79.3054% (the same as that of PSO), thus outperforming the CGA by 0.0021%, ALO by 0.0057%, and BH by 0.0704%. Figure 15. Percentage of reduction of minimum power losses obtained by the MVO in the 10-node mesh test system compared to other methodologies. Figure 16 compares the average P loss obtained by the MVO and the other optimization algorithms in the 10-node mesh test system (with the MVO as the base case). This analysis was performed employing the same three penetration percentages as in Figure 15. At 20% penetration, the MVO presents a reduction of average P loss of 44.9601%, outperforming ALO, the CGA, PSO, and BHO by 0.1485%, 0.1639%, 0.2987%, and 0.6524%, respectively. At 40% penetration, the MVO exhibits a reduction of average P loss of 69.2691%, outperforming the other optimization algorithms by 1.0645%. At 60% penetration, the MVO has the best reduction of average P loss again, outperforming the CGA, ALO, BHO, and PSO by 0.0428%, 0.1449%, 0.6585%, and 0.7125%, respectively.

10-Node Mesh Test System
To conduct the analysis of the mesh system, the precision and repeatability of all the algorithms are compared in Figure 17. As with the minimum and average P loss , the standard deviation is assessed with 20%, 40%, and 60% penetration. In the first case, the MVO presents an outstanding standard deviation of 0.002%, outperforming the CGA by 0.1152%, ALO by 0.1775%, BH by 0.5359%, and PSO by 1.8050%. In a 40% penetration scenario, the MVO has a standard deviation of 0.006%, outperforming the other optimization algorithms by a significant average of 6.46%. Finally, at 60% penetration, the MVO presents an outstanding solution precision with a standard deviation of 0.002%, outperforming the CGA by 0.1026%, ALO by 0.6884%, BH by 1.5748, and PSO by 10.3097%.
Analyzing the results obtained and discussed in this section, it is possible to note that, in general, the MVO provides the best solution. It offers an outstanding solution quality, short processing time, and a good standard deviation; as a result, it finds an excellent solution every time the algorithm is run. Therefore, the MVO is the most suitable optimization methodology to solve the OPF problem in AC networks of any type and size.

Conclusions
In this paper, we proposed the implementation of the MVO to solve the OPF problem in AC networks using a master-slave methodology. The master stage uses the MVO to determine the optimal power to be injected by each DG in the AC network. The slave stage employs SA (a numerical method) to assess the impact of every solution provided by the master stage on the objective function and the constraints of the problem. In order to demonstrate the robustness of the proposed algorithm, it was compared with some of the best optimization techniques that have been reported in the specialized literature to solve the OPF problem in AC networks: PSO, ALO, BH, and the CGA. All the algorithms employed in this study were tuned so that they could find the best solution to the OPF problem. To test each algorithm, the simulations were performed in the 10-, 33-, and 69-node radial systems and the 10-node mesh system using three levels of distributed generation penetration: 20%, 40%, and 60% of the power supplied by the slack generator in the base case. Moreover, each algorithm was executed 100 times to evaluate its standard deviation and precision.
According to the results in radial systems, the MVO obtained the best solution regarding the reduction of average P loss in an adequate computation time and presented the lowest standard deviation. Moreover, in most cases, the MVO achieved the best percentage of reduction of minimum P loss , although it was outperformed only by PSO by an almost negligible difference in the 10-and 69-node radial systems at a 20% penetration level. In terms of repeatability, the MVO presented low standard deviations (its highest value was 1 × 10 −3 ) and outperformed most of the other optimization techniques in the different test scenarios. Thus, this technique proves to be efficient in terms of precision and repeatability, obtaining good quality solutions to the OPF problem every time the algorithm is executed.
In the mesh topology, the MVO algorithm obtained the best solution in all the scenarios of distributed generation evaluated in this study, reducing minimum power losses by 0.029% and average power losses by 0.59% with respect to the other methods used here for comparison. The MVO achieved these results in a short processing time with the lowest standard deviation (under 0.006%). This proves its efficiency in terms of solution quality and processing times, as well as the repeatability of its solutions. Therefore, the MVO obtains a good-quality solution to the OPF problem in AC mesh networks every time its algorithm is run. From the paragraphs above, it is possible to conclude that the MVO is the best strategy to solve the OPF problem in mesh or radial AC networks.
Future studies can use this methodology for optimal power management in mesh or radial AC networks by employing power dispatch in energy storage devices and power loads considering their behavior in a 24-h cycle. This methodology can also be implemented in the planning of AC electrical systems to size distributed energy resources (generators, reactive elements, and energy storage) when the objective functions are to reduce the operating cost of the network, minimize CO 2 emissions, and improve the operating conditions. Finally, to reduce processing times and use parallel processing tools, future research can propose a parallel version of this methodology to shorten the processing times of the individuals that compose the populations at each iteration. This would considerably reduce the total processing times of the method. Funding: The authors would like to thank Minciencias (through the Fondo Nacional de Financiamiento para la Ciencia, la Tecnología y la Innovación, Fondo Francisco José de Caldas); the Instituto Tecnológico Metropolitano; the Universidad Nacional de Colombia; and the Universidad del Valle (through the project entitled "Estrategias de dimensionamiento, planeación y gestión inteligente de energía a partir de la integración y la optimización de las fuentes no convencionales, los sistemas de almacenamiento y cargas eléctricas, que permitan la generación de soluciones energéticas confiables para los territorios urbanos y rurales de Colombia," which is part of the research program entitled "Estrategias para el desarrollo de sistemas energéticos sostenibles, confiables, eficientes y accesibles para el futuro de Colombia").

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