Optimal Power Flow in Direct-Current Power Grids via Black Hole Optimization

This paper addresses the Optimal Power Flow (OPF) problem in DC power microgrids via a combinatorial optimization technique known as Black Hole Optimization (BHO). Such optimization method allows to solve OPF problems via algorithmic strategies trough a master-slave formulation. In the master stage, the total power generated by each Distributed Generator (DG) is determined by the BHO, while the slave strategy is entrusted with solving the resulting conventional power flow problem via a classical Gauss-Seidel (GS) numerical method. For comparison purposes, this work uses nonlinear optimization methods available in General Algebraic Modeling System (GAMS) as well as continuous metaheuristic optimization techniques. Two test feeders with 21 and 69 nodes were considered for validating the proposed hybrid BHO-GS optimization method, which enables to demonstrate its applicability, robustness and efficiency compared to conventional approaches. The results of all the simulations were obtained via MATLAB 2017a.


Introduction
Under the microgrids paradigm, Direct Current networks become a possible reliable and economic alternative for providing electrical services to millions of end-users around the world [1].The main advantage of using DC networks lies in the possibility of integrating multiple distributed energy resources, such as energy storage systems [2] (batteries, supercapacitors or fuel-cells) and renewable generation [3] (solar, wind technologies, among others) with DC-DC converters or AC-DC inverters, thus reducing the back-to-back topologies required to integrate them into conventional AC power grids [4].Those reductions must be reflected in the total cost of the distributed energy resources, as well as lower power losses and high electrical efficiency DC networks in comparison with their AC counterparts [5].Another important driver of the popularization of DC networks is that they are easier to analyze and operate because key concepts of AC grids, such as frequency and reactive power, disappear under this paradigm [5].
Power flow and Optimal Power Flow (OPF) problems have drawn strong attention in the field of DC networks because they represent essential tools to plan, operate, and control DC grids [5], [6] and [7].Aforementioned problems are interesting for the research community because they are nonlinear and non-convex, and solving them requires to propose new and efficient methodologies [8].Regarding the power flow problem for DC grids, [5] proposes a conventional Gauss-Seidel iterative method via successive approximation and demonstrates the Gauss-Seidel convergence via Banach space theorems.Besides, [7] demonstrates the convergence of the Newton-Raphson method through Kantorovitch's theorem for DC grids operating under a master-slave control mode.The authors of [9] present a linear approximation via Taylor's series expansion to solve the power flow problem in low-voltage DC microgrids with high performance and acceptable results in comparison to Gauss-Seidel and Newton-Raphson methods.In the case of OPF, the objective function of the mathematical models employed in the specialized literature is power loss minimization; in addition, conventional power flow, voltage regulation, and distributed generation capabilities are their constraints.Therefore, in order to analyse and solve this problem, [6] presents a second-order cone programming formulation, while [10] proposes a convex model for OPF problems in DC grids with higher penetration of renewable energy resources and energy storage systems, which is solved with the CVX solver.The main disadvantage of the aforementioned convex optimization methods is the quadratic increment in the number of variables required to solve the OPF problem, which can cause computational inefficiencies.On the other hand, [11] presents an approximated OPF model without increasing the number of optimization variables via Taylor's series approximation; notwithstanding, that solution generates approximated results that differ in about 2 % from exact models.In terms of combinatorial methods for DC OPF problems, a continuous genetic algorithm was recently proposed in [1] as a hybrid methodology in conjunction with a Gauss-Seidel numerical method.Such approach was only validated on a small 10-node test feeder, which does not allow to confirm its applicability to large test feeders.
According to the review of the state-of-the-art above, only one combinatorial optimization method has been applied to solve OPF problems in DC power grids, as proposed in [1].In that sense, this paper identifies a research gap in the field and proposes a Black Hole Optimization (BHO) approach [12], in conjunction with a classical Gauss-Seidel (GS) power flow method, to analyze and solve this problem by implementing a masterslave optimization strategy.The main advantage of the method proposed in this work, in comparison with convex approaches, lies in the fact that the number of variables of the power flow problem remains constant and no eigenvalue decomposition is required to recover power flow variables.Additionally, the proposed BHO-GS optimization method does not require specialized software, and it can be implemented in any programming language since its solution methodology is entirely algorithmic.
The remainder of this paper is organized as follows.Section 2. presents the conventional OPF problem for DC power grids with a focus on the possibility of including distributed generation in the grid via controlled percentages of penetration.Section 3. introduces the main characteristics of the BHO method as well as its evolution process and its application to the master problem.Section 4. shows the conventional GS formulation for solving power flow problems in DC power grids and its application to the slave problem.Section 5. details the complete BHO-GS method applied in this document to solve OPF problems in DC power grids as well as its pseudo-code.Section 6. describes the test system, the numerical results, and the General Algebraic Modeling System (GAMS) in conjunction with a large-scale nonlinear solver SCIP, a Continuous Genetic Algorithm (GA), Particle Swarm Optimization (PSO), and an Elephant Swarm Water Search Algorithm (ESWSA) as comparative methods.Finally, Sec. 7. draws the main conclusions and possible future research derived from this work, followed by the acknowledgments and references.

Optimal Power Flow Modeling
The optimal power flow problem in DC power grids corresponds to a nonlinear non-convex optimization problem [6], in which the power balance at all nodes in the grid is defined by hyperbolic constraints [8].The complete formulation of this problem is presented below: 2.1.Objective Function where v ∈ R n×1 corresponds to the vector that contains all the voltage profiles of the network, G L ∈ R n×n contains all the conductive effects associated with all the branches, and p loss is a scalar function that quantifies total power losses in the DC grid.Notice that n corresponds to the total number of nodes in the DC grid.

Set of Constraints
where p g ∈ R n×1 , p dg ∈ R n×1 , and p d ∈ R n×1 correspond to power generation at slack nodes (ideal voltage-controlled sources), distributed generators, and power demanded in constant power load nodes, respectively; p min g ∈ R n×1 p max g ∈ R n×1 represent minimum and maximum power capacities at slack nodes; p min dg ∈ R n×1 p max dg ∈ R n×1 denote the lower and upper limits of distributed generation nodes; v min ∈ R n×1 v max ∈ R n×1 are the minimum and maximum voltage regulation limits on the grid; α is a scalar parameter that represents the maximum power injection allowed via distributed generation, i.e., 0 ≤ α ≤ 1; finally, 1 T ∈ R 1×n and G L ∈ R n×n contain all the conductive effects associated with constant resistive loads.
The mathematical model detailed from Eq. (1) to Eq. ( 6) is interpreted as follows: Expression Eq. ( 1) corresponds to the objective function associated with power losses reduction in the electrical network.Power balance is presented in Eq. ( 2).Slack generation capabilities as well distributed generation capacities are established in Eq. ( 3) and Eq. ( 4), respectively.Voltage regulation bounds are presented in Eq. ( 5), while maximum distributed generation penetration is defined by Eq. ( 6).
Note that, for practical purposes, a constraint associated with power generation capacity at slack nodes is not considered for OPF problems since said nodes are considered ideal, i.e., they can maintain a constant voltage profile regardless of the power requirements of the grid.
Optimization methods are required to solve the OPF problem in DC power grids.Consequently, this paper explores the possibility of solving such problem via a master-slave strategy that decomposes the OPF problem into two sub-problems.The master subproblem determines the power generation of all the DG nodes via a BHO approach, while the slave counterpart resolves the resulting conventional power flow problem via a conventional Gauss-Seidel method.Those processes will be explained in detail in the following sections.

Master Problem: BHO Approach
Black hole optimization is a nature-inspired optimization technique based on the dynamical interaction between stars and black holes at the center of galaxies in the universe, as depicted in Fig. 1 [13].In general terms, this technique allows to solve nonlinear largescale optimization problems via heuristic exploration of the solution spaces [14].
Fig. 1: Typical behavior of a black hole and stars in the center of the Milky Way.Taken from [13].
The BHO method is briefly introduced below by highlighting the main steps in its computational implementation.

Stars Are Born
The BHO approach is a population-based optimization algorithm derived from conventional particle swarm optimization [15]; in that sense, the initial population corresponds to the first set of stars randomly distributed over the solution space, i.e., as a cumulus of stars in the universe [16].During the generation of this set of stars, the number of stars (possible solutions) n i is the number of rows, while number of DGs n dg is the number of columns, as formulated in: where P t represents the star population matrix, o(n i , n dg ) is a rectangular matrix filled with ones, and r(n i , n dg ) corresponds to a rectangular matrix filled with random numbers from zero to one with normal distribution properties.Note that P t corresponds to the current population and each individual (star) inside it represents the total power generation at all nodes that contain DGs.Furthermore, the best solution (lower fitness function for minimization problems) inside the initial population P t is selected as a black hole location.

Movement of Stars
The dynamic behavior of stars in the proximity of a black hole is highly influenced by the intense gravitational force of the latter (see Fig. 1).In that sense, the movement of any star may have a particular behavior as a function of its location with respect to the position of the black hole [14].Such behavior is emulated by the mathematical relationship where p BH t represents the black hole in population t, and P i t+1 denotes the i th individual after its movement.Remarkably, after this process the location of the black hole remains unaltered.

Black Hole Updating
After generating the descending population of stars P t+1 , the location of the black hole must be changed if the i th individual among the descending population exhibits a lower fitness function value than the current black hole, i.e., p BH t = p i t+1 ; otherwise, the location of the black hole remains constant, i.e., p BH t+1 = p BH t .

Star Replacement
The survival of a star in the nearby neighborhood of a black hole depends on its current position with respect to black hole's location.In theoretical physics, any object that crosses the event horizon around a black hole is destined to absolute destruction.Nevertheless, this catastrophic scenario generates stellar material that enables the formation of new stars-a marvelous event in the universe.To emulate the possibility that an arbitrary star in the descending population is absorbed by the black hole, the event horizon radius is defined as: where f (P BH t+1 ) represents the best fitness function value of all individuals contained in the current population (black hole objective function value), while the denominator of Eq. ( 9) corresponds to the sum of the fitness function of all individuals in the same iteration.
To determine if any star crosses the event horizon, the euclidean distance of such star with respect to the black hole's location is defined as: If R EH < D BH−i , a new star is randomly generated to replace the one absorbed by the black hole; otherwise, the star continues in the current population.Notice that the birth of new stars increases the possibility of expanding the exploration of the algorithm over the solution space, which would be considered global exploration [16].

Stopping Criterion
To halt the exploration of the BHO over the solution space, one of the following stopping conditions must be satisfied: • The maximum preset number of iterations has been reached.
• After k consecutive iterations the black hole's location has not been updated.

Slave Problem: Gauss-Seidel Method
Gauss-Seidel is a widely-known iterative method for solving power flow problems adopting the successive approximation theory [5].Such method allows to solve nonlinear nonaffine constraints associated with the power balance at all the nodes in the DC network as given in Eq. ( 2).For this purpose, said equation can be decomposed into two subsets.The first subset of equations is associated with slack nodes (voltagecontrolled nodes), and the second subset corresponds to the constant power loads including distributed generation nodes.This decomposition is expressed as: Note that Eq. ( 11) corresponds to a set of linear equations, since v i ∈ S denotes the voltage profiles at slack nodes, which are perfectly known.Nevertheless, an iterative process is required to solve Eq. ( 12).In that sense, a conventional Gauss-Seidel method is presented below ∀k ∈ D.
where m represents the iteration index.Note that the Gauss-Seidel method always converges for power flow studies, since |g kk | ≥ |g kj |∀k, j ∈ D, which implies that the conductance matrix is diagonally dominant [5].Moreover, the Gauss-Seidel method reaches numerical convergence if the following stopping criterion is fulfilled.
where corresponds to the maximum convergence error.

Solution Strategy: OPF Solution Via BHO-GS
When metaheuristic techniques are used for solving constrained optimization problems, the fitness function value is typically composed of the original objective function value as well as constrains added via penalties.The fitness function value proposed in this paper to guide the master optimization problem is written as: where β 1 to β 5 correspond to penalization factors, which are typically grater than zero.In this paper, each penalization factor equals 1000 in order to force the BHO to fulfill all the conditions imposed over the OPF problem as formulated from Eq. (1) to Eq. ( 6).
When all the constrains are fulfilled, all the penalization factors must be annulled by max{•} and min{•} functions, which turns the fitness function value into an objective function, since in that case z = p loss .
Algorithm 1 shows the application of the masterslave solution method proposed in this study to solve OPF problems in DC power grids via a hybrid BHO-GS approach.
Algorithm 1 Proposed pseudo-code for the hybrid BHO-GS approach.

Test Systems
Two different test feeders were used to evaluate and validate the proposed master-slave BHO-GS for solving the OPF problem in DC networks.First, a 21node test feeder described in the specialized literature is employed [9] and [11]; the second test feeder corresponds to a DC adaptation of the Baran & Wu 69-node test feeder implemented in [17] for optimal location and sizing of distributed generators in AC distribution networks.The information of both test systems is detailed below.

1) 21-node Test Feeder
This test system comprises 21 nodes and 20 lines with multiple constant power loads, which is an adaptation of the system originally presented in [7] and [9].In addition, said system includes two ideal voltage generators at nodes 1 and 21.The electrical configuration of the test system is illustrated in Fig. 2, and information about its power consumption and branches can be consulted in [11].The voltage and power bases in this test system are 1 kV and 100 kW, respectively.Note that the generator located at node 1 is operated with a voltage of 1.0 p.u, while the generator located at node 21 is dispatched with 1.05 p.u of voltage.Moreover, the power losses without DGs reached 0.211 p.u.

2) 69-node Test Feeder
This test system constitutes an adaptation of the AC 69-node test system employed for power losses reduction via distributed generation integration [17].To transform this system into a DC network, we used 12.66 kV and 100 kW as voltage and power bases; in addition, the reactance component in all the branches as well as reactive power consumption in all the nodes were neglected.Figure 3 presents the 69-node test feeder configuration; branch parameters and load information can be consulted in [17].For simulation purposes, 12.66 kV and 100 kVA were used as voltage and power bases for this test system.This study analyzes the possibility of installing from one to three generators considering penetration percentages from 20 % to 60 % of the total power consumption.Note that, via heuristic search methods, nodes 9, 12 and 16 were selected for GD location in the 21-node test feeder; and nodes 26, 61 and 66, in the 69-node test feeder.

Numerical Results
The computational analysis was carried out in a desktop computer with an INTEL(R) Core(TM) i5 − 3550 processor (3.50 GHz) and 8 GB of RAM, running a 64bit Windows 7 Professional operating system.We used the programming environment MATLAB 2017a and the GAMS optimization toolbox [18].

1) 21-node Test Feeder
Table 1 presents the total generation at each DG considering different percentages of power penetration, total power generation, and objective function values in the DC grid.
Note that, in each simulation case, the maximum power injection allowed into the DC grid reaches 1.1080 p.u, 2.2160 p.u and 3.3240 p.u for DG penetration scenarios of 20 %, 40 % and 60 %, respectively.In that sense, from the fifth column in Tab. 1 it can be seen that the proposed BHO-GS method as well as the GAMS optimizer use more than 99 % of the maximum generation allowed to reduce total power losses in the DC network.In terms of objective function minimization, note that the method proposed in this work presents an estimation error lower than 0.30 % (compared to the SCIP solver) for each simulation case.Therefore, such hybrid BHO-GS method offers adequate numerical convergence compared to large-scale nonlinear packages.
It is important to highlight that the power delivered by each distributed generator exhibits small variations (compared to the GAMS optimizer) in the solutions obtained by means of the proposed approach.These differences can be attributed to the nonlinear continuous nature of the OPF problem, which implies that multiple DG power combinations with identical objective functions could exist.Finally, Fig. 4 shows the reduction of power losses in the DC network achieved by the proposed and comparative methods.The proposed BHO-GS method and the comparative method via GAMS implementation exhibit a similar performance in terms of power losses reduction, which clearly validates the proposed approach in terms of numerical convergence.

2) 69-node Test Feeder
For validating the proposed BHO-GS methodology in this test feeder, we employed three different optimization techniques usually implemented in the specialized literature for solving continuous optimization problems: Particle Swarm Optimization (PSO) [17], an Elephant Swarm Water Search Algorithm (ESWSA) [20], and a continuous genetic algorithm [1].To compare all the metaheuristic techniques fairly, we considered a population size of 30 individuals, a maximum number of iterations of 200, and a stopping criterion of 50 (continuous iterations without improving the objective function) to parameterize those techniques.In this test feeder, we assumed that the location of the distributed generators corresponded to nodes 26, 61, and 66 as recommended in [17].Likewise, the total power injection allowed into the system is defined as 60 % of the total power consumption, and each generator's capacity is considered free.
Table 2 presents the comparison of OPF solutions obtained employing the aforementioned metaheuristic techniques for the 69-node test feeder.Note that the PSO algorithm has the best objective function, with 0.05556 p.u. of power losses, followed by the proposed BHO algorithm, with 0.05571 of power losses reduction.In that sense, the PSO algorithm reduced Tab.1: Optimal power dispatch of DGs using the BHO-GS approach and GAMS optimization package [19].It is also important to stress that the proposed BHO reached the lowest power injection into the grid, 22.0233 p.u., followed by the ESWSA with 22.0388 p.u. Furthermore, the AG remains the worst technique in relation to this indicator, followed by the PSO algorithm.

Solutions provided by the BHO-GS method proposed in this work
Regarding power generation at each node, note that node 61 concentrates the highest percentage of power injection according to all the optimization algorithms.Nevertheless, we can say the BHO algorithm outperforms its rivals because it injects 15.6109 p.u., while the ESWSA is in the second position with 15.8613 p.u.
The results above indicate that the proposed BHO-GS algorithm presents an adequate performance in comparison to other continuous optimization techniques for solving the optimal power flow problem in DC networks.

Conclusions
This paper presents a combination of the black hole optimization approach in conjunction with the conventional Gauss-Seidel numerical method for solving the OPF problem in DC power grids.The BHO method is a soft variant of the widely-known particle swarm optimization methods, which offer an adequate performance for OPF problems, as confirmed by the numerical results reported in this work.
The large-scale nonlinear optimization package SCIP solver (available for GAMS) and continuous optimiza-tion techniques (PSO, GA and ESWSA) were used as comparative methods, which confirmed the excellent numerical convergence of the proposed hybrid BHO-GS method.
In future works, the BHO-GS method presented in this paper could be used for solving the optimal location and sizing of distributed generation in DC power grids via binary metaheuristic techniques, such as genetic algorithms or tabu search methods, among others.

About Authors
Orfilio Sebastian VELASQUEZ was born in Cartagena, Bolivar, Colombia.He is currently enrolled in the undergraduate Electrical and Electronic Engineering program at Universidad Tecnollogica de Bolivar, Cartagena.His research interests include power systems analysis and optimization techniques applied to the planning, operation and control of electrical networks.

Fig. 4 :
Fig. 4: Power losses behavior with different levels of DG penetration.

1 :
Data: Initialization parameters.2: for t = 1 : t max do Fernando GRISALES was born in Cartago, Valle, Colombia.He received his B.Sc. and M.Sc.degrees in Electrical Engineering from Universidad Tecnologica de Pereira, Colombia, in 2013 and 2015 respectively.He is currently pursuing a Ph.D. in Automation Engineering at Universidad Nacional, Manizales Colombia, working as a professor at the Department of Electromechanics and Mechatronics at Instituto Tecnologico Metropolitano de Medellin, and a member of the research group MATyER.His research interests include mathematical optimization, planning and control of power systems, renewable energy, energy storage, power electronics and smartgrids.Oscar Danilo MONTOYA was born in Obando, Valle, Colombia.He received his B.Sc. and M.Sc.degrees in Electrical Engineering from Universidad Tecnologica de Pereira, Colombia, in 2012 and 2014 respectively.He is currently pursuing a Ph.D. in Electrical Engineering at Universidad Tecnologica de Pereira, Colombia, and working as an assistant professor at the Electrical Engineering program at Universidad Tecnolgica de Bolivar.His research interests protective devices, passivity-based control, and dynamical analysis.Victor Manuel GARRIDO was born in Sincelejo, Sucre, Colombia.He obtained his B.Sc. in Electrical Engineering in 2009 from Universidad Tecnolgica de Bolivar, Cartagena, Colombia.He received his M.Sc. in Control at Universidad de Pamplona in Colombia.His M.Sc.thesis was about power quality and intelligent techniques.His research interests include power systems, transmission and distribution networks, and power quality.
c 2019 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING