Optimal Location and Sizing of Distributed Generators in DC Networks Using a Hybrid Method Based on Parallel PBIL and PSO

: This paper addresses the problem of the locating and sizing of distributed generators (DGs) in direct current (DC) grids and proposes a hybrid methodology based on a parallel version of the Population-Based Incremental Learning (PPBIL) algorithm and the Particle Swarm Optimization (PSO) method. The objective function of the method is based on the reduction of the power loss by using a master-slave structure and the consideration of the set of restrictions associated with DC grids in a distributed generation environment. In such a structure, the master stage (PPBIL) ﬁnds the location of the generators and the slave stage (PSO) ﬁnds the corresponding sizes. For the purpose of comparison, eight additional hybrid methods were formed by using two additional location methods and two additional sizing methods, and this helped in the evaluation of the effectiveness of the proposed solution. Such an evaluation is illustrated with the electrical test systems composed of 10, 21 and 69 buses and simulated on the software, MATLAB. Finally, the results of the simulation demonstrated that the PPBIL–PSO method obtains the best balance between the reduction of power loss and the processing time. Future improvements to this study could consider the inclusion of renewable generation curves and hourly power demand into the mathematical model, which could provide a more accurate solution for a particular geographical region. Moreover, the solution can be further developed to integrate energy storage devices into the electrical network, which could be used for improving the electrical behavior and proﬁtability of the grid. In addition, it is possible to propose an optimal tuning strategy for the parametrization of the continuous optimization techniques with the aim of balancing the exploration and exploitation of the solution space, which will enable the identiﬁcation of the best compromise between the number of iterations and the population sizes.


General Context
Research and implementation of the DC networks have been increased considerably Track changes is off Everyone Track changes for everyone You Track changes for You Guests Track changes for guests Current file Overview 3 due to the advantages of DC networks when compared to AC networks in terms of reliability, efficiency, and simplicity of control [1][2][3]. The distributed generators in DC networks are usually interconnected within the electrical system through electronics interfaces, which are intended to control the operation and condition of the generator. Those power electronics interfaces could use batteries to guarantee a given power flow from non-predictable renewable generators (e.g., photovoltaic and wind based) [4,5], or they could be controlled so that they dispatch the required power from predictable sources, such as fuel cells (by using a DC/DC converter) [6] and diesel generators (by using a rectifier) [7].
Moreover, DC networks exhibit a simpler process of analysis and operation, compared to AC networks, as the frequency and the reactive power are not present in the mathematical analysis [8,9]. In this sense, the problem of power flow on DC networks has been studied by using different numerical methods, such as Gauss-Jacobi, Gauss-Seidel and Newton-Raphson [10], linear approximations [11], semi-definite and second order cone programming formulations [12,13] as well as convex quadratic models [14,15].

Motivation
The optimal power flow analysis is one of the main tools used to determine the steady-state conditions of the electrical networks in the process of the analysis of electrical networks operated under the DC paradigm. In this context, the authors of [16] proposed a hybrid methodology based on a Genetic Algorithm (GA) and the Gauss-Seidel method. Similarly, the authors of [17] used a Particle Swarm Optimization (PSO) as the solution method, which considered the integration of Distributed Generators (DGs) and batteries. The solution proposed in that work is analyzed under different test cases; nevertheless, the processing times required in each case are not analyzed. With respect to the optimal integration of DGs in DC grids in [18], the authors have proposed a methodology for integrating those DGs into DC grids, which used a heuristic approach to locate the generators and a mathematical method for determining the size of the photovoltaic generators. Such a solution considers the reduction of the installation and operation costs as its objective function, and this is optimized by using the linear programming tool linprog available for MATLAB [19]. However, no comparison with other methods is provided, and the processing times are not analyzed. Similarly, in Reference [20], the PSO has been used to calculate the optimal sizing and operation of DGs and batteries in DC grids. In this work, the objective function was the reduction of the investment cost and the operation cost of the energy resources by considering a unique nodal representation of the grid that does not consider the multiple buses and branches. In addition, the impact of the distributed energy resources (DERs) on the DC grid was analyzed without the use of the power flow methods or tools and the comparison with other methods reported in the literature.
Therefore, it is necessary to propose mathematical formulations that represent the problem associated with the integration of DGs in DC networks, while considering multiple buses, loads and DGs. Such a model can be used to develop optimization techniques that are able to identify the configuration of the generators with the best impact on the technical aspects of the network with short processing times, which enable them to guarantee the time restrictions of energy projects with short horizon time. In this sense, a correct selection of the method to be used for solving the power flow problem is critical. In particular, the processing time required by such a method is very important, since the power flow must be solved for each configuration evaluated by the optimization algorithm. Based on the state-of-the-art review in this paper, it is observed that the optimal integration of DGs in DC grids by using parallel processing tools have not been addressed; therefore, such a topic must be studied to ensure faster implementation.

Literature Review
The problems associated with the optimal location and sizing of DGs in DC networks have been few explored in the literature, as this topic has become increasingly relevant in the last decade due to the important developments in power electronic interfaces, energy storage systems and renewable generation [10]. Some of the recent work related to this are summarized below: the authors of [21] have proposed a mixed-integer quadratic programming (MIQP) model to address the problem of the optimal location and sizing of dispatchable DGs in DC grids where the proposed model is solved by using optimization solvers available in the general algebraic modeling system (GAMS) and the numerical results demonstrate the effectiveness and robustness of the exact mixed-integer non-linear programming (MINLP) model; however, the global optimal finding is not possible due to the non-convexity of the solution space. Exact mixed integer non-linear programming methods have been proposed in [22,23] to represent the problem of the integration of DGs in DC grids while taking into consideration the reduction of the power loss as the objective function. To solve that problem, a convex relaxation of the mathematical formulation is performed by using a Taylor's series expansion to implement a transformation of variables. The main problem of that approach is the high possibility of converging into a local optima. This issue has been few explored even in DC networks the studied problem, but this is not the case of the conventional AC grids where this topic has been largely investigated, as seen in [24][25][26][27]. This implies that it is possible to evaluate some of the solutions applied in AC networks in the context of their application in DC networks. For example, in Reference [28], the authors proposed a hybrid method for optimal integration of DGs in electrical systems based on a GA and the PSO has been proposed. In that study, the reduction of the active power loss was used as the objective function, thereby validating the proposed method in the test systems of 33 and 69 buses. In this paper, the authors do not compare the performance of the proposed method with other solutions reported in the literature, and the processing times associated with the optimization methods are not discussed. Similarly, in Reference [29], a hybrid methodology that uses the sine cosine algorithm and chaos map theory has been proposed for solving the problem that is analyzed here. This methodology was validated by using the 33 and 69 bus test systems and by comparing the results obtained with other optimization techniques reported in the literature; however, the processing times were not discussed. Another approach has been reported in [30] where a master-slave strategy has been implemented by using a modified version of the PSO to locate the DGs and the capacitors, and the Monte-Carlo algorithm has been used for defining the size of the generators. In that study, the main objectives used were to reduce the power loss and improve the voltage profiles, but the authors have not provided a comparison with other techniques or an analysis of the processing time.
Another approach for AC networks has been presented in [31] where the main objective was the maximization of the energy cost associated with the generation of the distributed energy resources. In Reference [32], a metaheuristic approach based on the population based incremental algorithm and the PSO has been proposed with the aim of locating DGs by using a parallel processing system to improve the calculation speed. In that study, the objective function is based on the active power loss, and the authors provide comparisons with other techniques only for the step of location by using the same sizing method for all the tests, which does not allow them to evaluate the effectiveness of the location methods in comparison with other sizing methods proposed in the literature. Finally, an example of a parallel processing system that can be applied to the problem of the integration of DGs in AC grids has been addressed in [33]. That particular study uses a parallel version of the Monte-Carlo algorithm (PMC) by using the processing time as the performance indicator. However, a comparison with other methods is not provided in it. In conclusion, the previous references provide evidence on the effectiveness of the integration (location and sizing) of DGs in AC networks based on the reduction of power loss while also using the processing times as an efficiency indicator.

Contributions and Scope
This paper proposes a methodology for the optimal location and sizing of DGs in DC grids, which is performed under peak load conditions in order to determine the optimal process of the integration of the DGs, thereby enabling the reduction of the power loss associated with the transport of energy on the grid in the worst operation scenario. The aforementioned scenario was selected in this paper, as it is widely used in the literature for evaluating the proposed methodologies that attempt to solve the problem associated with the optimal integration of DGs in electrical systems [28,34,35]. The solution proposed in this paper is based on a hybrid methodology between the PPBIL and the PSO algorithms that enables the improvement of the quality of the solution and the reduction of the processing time. The selection of the PBIL method for a parallel implementation is motivated by the low memory consumption and low computational complexity of that algorithm [36]. Furthermore, the computational characteristics of the PBIL simplifies the process of performing the parallel implementation aimed at reducing the processing time. The PSO has been widely adopted in the literature to calculate the size of DGs in both AC and DC grids [17,28,30], providing satisfactory results in terms of the solution quality and processing time. Therefore, based on the aforementioned results, the PSO was selected in this paper for calculating the size of the DGs. The main contributions of this work are summarized as follows: • The implementation of a mathematical formulation that considers the characteristics of DC grids, including the set of constraints associated with the representation of DC electrical networks with distributed generation, which is needed to perform an optimal integration of the distributed generators for the minimization of power loss. • A methodology for the optimal integration of DGs in DC grids that uses a parallel master-slave strategy to obtain the satisfactory results in terms of the quality of the solution and the processing time.

•
The adoption of parallel processing tools for the planning of DC grids in order to improve the performance of the methodologies for sizing and location.

•
The adoption of sequential programming in order to enable the use of free programming languages for the optimal sizing and location of DGs in DC grids and avoid the use of specialized and costly software.
It is worth mentioning that this hybrid solution has been previously proposed in [32] in the context of AC distribution networks, but it has not been applied to DC systems. Furthermore, the work reported in [32] has only evaluated the performance of the locating technique. Therefore, the PSO algorithm has not yet been demonstrated as the most suitable solution for the development of the hybrid methodology for DC networks. Instead, this study considers two additional location techniques and two additional sizing techniques, which were used for demonstrating the improvement created by the PPBIL-PSO in DC grids. The two additional location techniques are the GA [28] and the PMC [33] algorithms. The selection of those solution methods was based on the satisfactory results reported in the literature for the same problem [34,[37][38][39]. The two additional sizing techniques are based on a continuous GA proposed in [16], and the black-hole (BH) optimization method has been proposed in [40]. As in the previous case, those methods were selected based on the satisfactory results reported in the literature.
The study described in this paper considers the combination of the three location techniques with the three sizing techniques, thereby obtaining nine hybrid solutions, one of them being the proposed PPBIL-PSO method. The evaluation of these solutions is performed in three test systems composed of 10, 21 and 69 nodes, respectively, while considering that a maximum of three DGs are to be located. The simulations were carried out in MATLAB by using the successive approximation method for evaluating the power flows inside the sizing methods [14]. The results obtained demonstrated that the proposed mathematical formulation and the methodology of the solution were able to solve the problems concerning the optimal location and sizing of DGs in DC grids, as described in the previous paragraphs.

Organization of the Document
The remainder of this document is organized as presented below: Section 2 shows the mathematical formulation used for the optimal integration of DGs in DC grids. Then, Section 3 presents the master-slave strategy based on both the PPBIL and the PSO algorithms proposed in this paper. The test scenarios and considerations used are explained in Section 4. Then Section 5 presents the results of the simulations and its analysis. Finally, the conclusions and the possibilities of future research are described in Section 6.

Mathematical Formulation
The problem of the optimal location and sizing of DGs in DC networks was formulated by means of an optimization problem mono objective where the objective function takes into consideration the reduction of the power loss in the whole system, and this is subject to a set of technical constraints associated with the DC networks in a distributed generation environment [11]. The power loss is selected as the objective function, as this technical aspect is widely used in AC and DC networks and are associated with the same problem [16,17,[41][42][43].
The mathematical model that describes the problem of the integration of DGs in DC grids is formed by: Objective function: In such an objective function, P Loss represents the power loss related to the transport of energy by the branches where N denotes the set of buses that form the system, G ij is associated with the ijth component of the conductance matrix and G i0 corresponds to the conductance associated with the resistive load connected at bus i. Finally, v i and v j represent the voltages in the buses i and j, respectively.
Set of constraints: The mathematical formulation defined from (1) to (8) considers the minimization of the objective function defined in (1), and therefore, the minimization of the total power loss of the DC grid. Equation (2) imposes the power balance at each bus where P g i and P d i denote the power generated and demanded at bus i, respectively. Equation (3) imposes the voltage bounds where (v max i ) and (v min i ) correspond to the maximum and the minimum limits for the bus voltages. Expression (4) shows the current bound of each branch installed on the DC system where i ij is the current of the line ij within the set of branches that form the electrical network (B), and I max ij is the maximum current allowed in that branch. The power limits to be injected by the DG connected at bus i are defined in (5) where D denotes the set of buses selected for the installation of DGs. In this equation, x DG i represents a binary variable that takes a value of 1 when a distributed generator is installed at bus i, and it takes a value of 0 otherwise; the nature of x DG i is described in Equation (8). Finally, Equation (6) presents the maximum number of DGs that can be located on a DC grid (NDG max ), while Equation (7) represents the maximum level of penetration (P max DG ) allowed into the electrical network.
Finally, this optimization problem is non-linear and non-convex, and therefore, this paper adopts a master-slave structure to decouple the problem of the location of DGs (discrete optimization) from the problem of the sizing of DGs (continuous optimization).

Proposed Hybrid Method
The integration of DGs in electrical grids is traditionally divided into two sub-problems [28,41,42,44]: the location of DGs on the electrical network and the calculation of the size of those generators. This paper adopts a master-slave structure to connect the proposed PPBIL and PSO techniques. The PPBIL is entrusted with solving the discrete optimization problem of locating the DGs, and the PSO algorithm is responsible for solving the continuous optimization problem of sizing the DGs.

Codification of the Solution
The codification scheme for the location and sizing problems is depicted in Figure 1. The integration of the DGs is performed by using two vectors of (|N | − 1) columns where each element of the vector corresponds to a particular bus of the electrical system while excluding the voltage controlled node, i.e., the slack bus. |N | represents the number of nodes present in the DC grid. Figure 1a shows the binary codification for the location of the DGs in which it is assigned a value of 1 when a bus is a candidate for the locating of a generator and 0 is assigned in the opposite case. For example, the vector in Figure 1a reports that buses 1, 3 and (|N| − 2) were selected for locating DGs.  . The example given in Figure 1b reports that buses 1, 3 and (|N | − 2) will supply 0.85, 1.25 and 0.75 p.u, respectively.
Finally, both vectors are used in the master and slave processes, which are described in the next subsections.

Master Stage: Optimal Location of DGs by Using PPBIL
In the master stage, the location of the generators is performed by using the PPBIL algorithm; this metaheuristic optimization algorithm is selected due to the satisfactory results provided by this in the case of AC networks [32]. The PPBIL includes a modification of the iterative processing of the traditional PBIL [45] in order to enable parallel processing, which reduces the computation time required by the algorithm.
The PBIL algorithm is an estimation of the distribution algorithm [46], and it uses probabilistic methods to find a solution with a given value on the objective function. This algorithm has been widely used to solve combinatorial optimization problems [47,48], as is the case of the location of DGs. The main advantages of using the PBIL algorithm are the low memory consumption, low computational complexity and mainly, the capability of modifying the learning rate, which enables it to increase or decrease the exploration in the solution space in agreement with the application [49]. The PPBIL is implemented by parallelly evaluating the objective function of all the individuals of the population [32]. Similar approaches have been used in the literature in order to reduce the computation times of the optimization and simulation processes by taking advantage of the graphics processing units (GPU) or processors with parallel cores [33,50,51].
Algorithm 1 describes the PPBIL algorithm. In the following section, each step of the algorithm is explained in detail.
Algorithm 1: Pseudo-code proposed for the PPBIL algorithm Data: Step 1 → Assign initial conditions Step 2 → Initialize probability matrix; Step 3 → Generate the population; Step 4 → Evaluate the fitness function; Step 5 → Select the best individual; Step 6 → Update the probability matrix; Step 7 → Update the learning rate; Step 8 → Calculate the entropy; Step 9 → Evaluate the stopping criteria; end Step 10 → Extract the best individual from the probability matrix; Step 11 → Evaluate the fitness function;

•
Step 1. Assign initial conditions: In this step, four main parameters are selected to start the iterative process of the PBIL. First, the population size (PS) is assigned to it. In order to ensure that the evaluation of the complete population is done in one cycle of the parallel units, the PS selected is kept equal to the number of cores of the processor. The second parameter is the initial probability of each option in the probability matrix. For an adequate exploration of the solution space, the PBIL starts the iterative process with the same probability for all the possible options. The initial probability for this application is equal to 0.5 due to which each bus has two options: the option of locating (50%) or not locating (50%) a DG on the bus. The third parameter includes the type of learning rate (LR), for example, linear, exponential, sigmoidal or bell shaped. The LR bounds are within the range of 0 to 1 [49]. In regard to the location of DGs, it was found in a heuristic form that the best type of LR is sigmoidal and the best values for the LR bounds are 0.25 and 0.50, respectively. The fourth parameter is the stopping criterion, which corresponds to the entropy tolerance, and it is selected as 0.1 in accordance with the work reported in [32].

•
Step 2. Initialization of the probability matrix (PM): The PM is formed by (n) columns, which represent the elements to be considered while formulating the solution of the problem, and (m) rows, which represent the possible solutions for each element.
For the location of DGs, n is equal to (|N | − 1), thereby excluding the slack bus, and m = 2 since only two options are available for each node. Figure 2 shows an example of the probability matrix. At the start of the iterative process, all the components of the PM start with the same probability of 1/m = 0.5. Moreover, the PM must satisfy the constraint given in (9) where P(j, h) represents the probability of the option j to be selected on element h. Such a restriction ensures that the accumulative probability of the options for a single element is 100%. • Step 3. Generate the population: The population is generated by using the PM data. The objective is to create individuals by using the information of each option with each element of the PM to avoid identical individuals and guarantee an adequate exploration of the solution space. Therefore, when an individual is repeated, it is replaced by a randomly generated one.

•
Step 4. Evaluate the Fitness Function (FF): The evaluation of the FF of each individual is the step that takes up more time in the PBIL process. Therefore, the solution proposed in this work uses the parallel processing structure reported in [32], which enables the evaluation of the FF of all the individuals at the same time. If the population size is not equal to the number of cores in the processor (W), the time required by the paralleling process (PPT) is equal to PPT = CEIL(PS/W) · MTRP where MTRP is the maximum time required to evaluate the FF of an individual and CEIL (·) returns the integer that is higher than or equal to a real number. In this paper, the evaluation of the FF is performed in the slave stage, which optimizes the size of the distributed generators for each individual of the PPBIL population.

•
Step 5. Select the best individual: The individual with the best FF value is selected. In this case, the individual with the minimum power loss fulfills the set of constraints discussed in Section 2.
• Steps 6 and 7. Update the probability matrix and the learning rate: The values of the probability matrix are updated by using Equation (10) where P(i, j) Old is the non-updated probability of positions (i, j) and P(i, j) Act is the updated probability.
Then, the index k is assigned to the option with the highest fitness value of the element j, and the update of the probability matrix is performed by using Equation (11); the probability of the option with the highest fitness value is increased while the probability of the other options (different to k) are reduced. In such an equation, P(i, j) New denotes the actualized value of the probability at positions (i, j). Finally, all the elements must satisfy the restriction (9) after updating the probability matrix.
Algorithm 1 reports that the learning rate must also be updated. The value of LR is updated by using Equation (12) where LR min and LR max are the bounds assigned in Step 1. Moreover, E n is the entropy of the PM in the first iteration of the algorithm E n = 1; therefore, the probability matrix exhibits a high dispersion of the data, as all the options have the same probability in each element [45,49].
• Steps 8 and 9. Calculate the entropy and evaluate the stopping criteria: The entropy E n enables the quantification of the distribution of the probabilities of the solution.
In that way, the probabilities in PM are completely dispersed when E n = 1, and they converge into a single solution when E n = 0. However, a particular tolerance E TOL is used as an approximation to E n = 0 in order to allow converging into a solution in a finite time. Equation (13) describes the mathematical formulation of the entropy, which is calculated by adding all the probabilities in PM and using the total number of elements n to normalize the entropy into 0≤E n ≤1. The stopping criterion for the PPBIL algorithm is E n < E TOL with E TOL = 0.1. The selection of the E TOL value has been discussed in [32].
• Steps 10 and 11. Extract the best individual from the probability matrix and evaluate the FF: When the PPBIL algorithm reaches the stopping criterion, the best solution is extracted from PM by selecting the options on each element with the highest probability. Finally, the FF of the best configuration is calculated, which enables the evaluation of the impact associated with the best solution on the DC network.

Slave Problem: Optimal Sizing of DGs by Using PSO
In the slave stage, the PSO is used to determine the size of the generators present in each individual generated by the master stage, and it is a bio-inspired algorithm that uses the flocking behavior of birds and fish swarms for solving continuous problems [52]. This optimization algorithm was selected for solving the problem of the sizing of DGs in AC networks due to the satisfactory results reported in the literature [28,30,44,53]. In such a context, a detailed description of the PSO algorithm is given in [52].
The processing procedure of the PSO is summarized in Algorithm 2: the initial conditions of the PSO correspond to the number of particles, the velocity limits, the inertia factor, the cognitive and social components and the stopping criterion. The stopping criteria used for this paper are associated with the maximum number of iterations and a number of iterations in which the FF is not improved. In the first iteration, the particle swarm is generated by assigning the active power to be supplied by DGs as the coordinates of position. Therefore, the number of coordinates for each particle is equal to the number of generators to be installed. Subsequently, the power flow is solved and the FF for each particle is evaluated. Then, such information is used to select the best solution for each particle and the best swarm solution, which corresponds to the solution for the sizing problem. In the subsequent iterations of the algorithm, the velocity vector is calculated in each iteration, which is used to update the position of the particles while also considering the last position, the inertia factor, the cognitive and social components as well as the previous best positions for each particle and the swarm. In the next step of the PSO, the FF is evaluated for all the particles, thereby updating the position and the value of the best solution for each particle and for the swarm. Then, the stopping criteria are evaluated. Finally, the best solution reached by the PSO algorithm is sent back to the master stage.
Algorithm 2: Pseudo-code for the PSO algorithm Data: Assign initial conditions for t = 1 : t max do if t == 1 then Generate the swarm particles; Solve the load power flow for each particle; Evaluate the fitness function; Select the solution and position presented by each particle as the best particle solution and position; Select the best solution presented by the swarm particles and its associated position as the best swarm solution and position; else Calculate the velocity vector; Update the position of the swarm particles; Solve the power flow for each particle; Evaluate the fitness function; Update the best particle solution and position; Update the best swarm solution and position; if Have any stopping criteria been met? then solution found; Finish the optimization process;

Break; else
Continue; end end end Finally, the flowchart reported in Figure 3 summarizes the master-slave method based on both the PPBIL and the PSO algorithms. The figure proves that parallel processing is used to evaluate the objective function of all the PPBIL individuals at the same time. This is performed by the PSO algorithm to obtain the optimal size of each generator in each configuration.

Test Systems
In this manuscript, three different test systems are used for validating the hybrid method proposed for the optimal integration of DGs into DC grids.
The first two test systems are adapted from the 10 and 21 bus systems that are widely used in the literature to study DC networks [2,[10][11][12]. The modifications made to those test systems include the replacement of the distributed generators and batteries by loads with the same power magnitude. Therefore, the resulting test systems are only formed by a main generator and electrical loads; Figures 4 and 5 present the system architectures. The objective of the previous modifications is to increase the technical problems on the grids (increase power loss, reduce voltage profiles, etc.) and to eliminate the presence of distributed elements to provide a base case.
The DC version used in the third bus test system is an adaptation of the classical 69 bus test system, which is widely used in AC networks to evaluate the planning and operation strategies [32,[54][55][56]. The modifications consist of eliminating the reactance associated with the branches and the reactive power generated and demanded in the electrical system, thereby obtaining a DC version of the system. This test system enables the testing of the proposed solution in a large DC network. Furthermore, the 69 bus test system was modified to increase the power loss and to reduce the voltage profiles. The main characteristics of the three test systems are presented in the following subsections.

10 Bus Test System
This test system corresponds to the DC electrical network proposed in [8]. It is composed of 10 buses and 9 branches, and the electrical configuration is illustrated in Figure 4. The information about the branches, power supply and demand are shown in Table 1. The system considers two types of loads: constant power loads (P) and resistive loads (R), which produce an initial power loss of 0.1436 p.u. for a total demand of 3.6 p.u. Moreover, the slack generator (Gen) is set with a voltage of 1.0 p.u [11]. The voltage and power bases used in this system are 1 kV and 100 kW, respectively. The 21 bus test system presented in Figure 5 is consists of 21 nodes and 20 lines, and it considers the constant power loads [14]. Moreover, the generator is operated at a voltage of 1.0 p.u, the total power demanded is 5.54 p.u and the power loss is of 0.2760 p.u. In addition, the same bases used for the 10 bus test feeder are considered for this test system [2]. Finally, Table 2 shows the load, branches and bus data used for the 21 test system [2].  The modification made to the classical 69 bus test system to obtain a DC version consists of eliminating the reactive component in both the branches and the power consumption for all the nodes. On the other hand, the system preserves the number of nodes and branches, the resistance values of branches and the power load and voltage level reported in [57]. Figure 6 presents the line diagram of the grid. The test system considers 12.66 kV and 100 kW voltage and power bases, respectively. Finally, the power loss was found to be 1.5385 p.u for a total power consumption of 38.8925 p.u.

General Considerations
The following considerations and conditions are imposed to provide a fair comparison between the different optimization solutions [43]: • All nodes except the slack node are candidates for locating DGs.

•
It is considered that NDG max = 3, and therefore, a maximum of three generators can be installed. Such a limit is in agreement with the test scenarios adopted in the literature for the same problem [32,41,55].

•
The maximum power level that can be supplied by each generator is 1.2 p.u for the 10 bus system, 1.5 p.u for the 21 bus system and 12 p.u for the 69 bus system. The minimum power level of each generator is 0 p.u. for all the cases [15].

•
The current ratings of the DC lines for all the test systems were calculated by using a power flow analysis without taking the integration of DGs into consideration. The calibers and maximum current levels were defined based on those values of the branch currents and by using In accordance with the recommendation given in [58], all test systems were assigned a maximum distributed power penetration level equal to 40% of the total power supplied by the slack bus without considering the integration of DGs (base case). Such a constraint was considered in order to account for the restriction of the distributed generation present in some national electrical regulations; hence, the proposed mathematical formulation provides a general solution because such a limit could be set to zero when it is not needed. Moreover, this consideration enables the analysis of the impact of the penetration level of the DG on the electrical system.

Fitness Function
The metaheuristic techniques that are used to solve the constrained optimization problems can adopt a Fitness Function (FF) in order to satisfy the set of constraints that represent the problem and explore outside of the feasible solution space [16,35] in search of good quality solutions. The FF adopted in this paper is formed by the original value of the objective function and all the constraints introduced as penalties: In the previous FF, the constants β 1 to β 7 correspond to the penalization factors. In this study, all the penalization factors were set to 1000 by using trial and error in an attempt to force the optimization methods to satisfy all the constraints that compose the mathematical formulation defined in Section 2. When all the constraints are fulfilled, the penalization factors must be nullified by the max{·} and min{·} functions. Therefore, the fitness function becomes FF = min (P loss ).

Comparison Methods and Parameters
Two additional binary optimization methods were used to locate the DGs in the master stage, which enabled the evaluation of the performance of the PPBIL solution. Those binary methods were the GA and the PMC. The GA is widely used in the literature to locate DGs in electrical grids [28,44,59]. This technique uses the selection, recombination and mutation methods to find a combination of generators with lower power loss. The second binary method (PMC) implements a Monte-Carlo algorithm, which performs several random samples of the different possible locations for the DGs, simultaneously evaluating the FF of all the configurations in order to detect the locations with lower power loss [33]. The adoption of the PMC for evaluation purposes enabled the comparison of the proposed solution with another parallel processing method. The parameters of the PPBIL, GA and PMC methods used for the master stage are given in Table 3. The size of the population for both the PPBIL and the PMC were equal to the number of workers available in the PC processor. Otherwise, the processing time was significantly increased since the entire population can not be evaluated at the same time. Instead, the size of the genetic algorithm population was defined by using trial and error and selecting the number of individuals that can provide the best balance between the power losses and the processing time. Finally, the other parameters of the binary methods were obtained by using a heuristic method in accordance with the methodology reported in Reference [32]. For the evaluation of the slave stage as well, two additional optimization methods were considered. The first method is a continuous adaptation of the GA [16], which uses the selection, recombination and mutation methods to obtain a continuous population that represents the sizes of the distributed generators located in the master stage. The second method is the BH, which is a nature-inspired optimization technique. The BH algorithm is based on the dynamical interaction between stars and black holes in the universe [60]. This technique was used in [61] to solve the problem of the sizing of DGs in DC grids by implementing a particle swarm and some application-based criteria for the elimination and generation of stars through a heuristic approach.
The parameters for each continuous method used in the tests are reported in Table 4. For the proposed PSO, the parameters were found in a heuristic form in order to provide a satisfactory balance between the number of iterations, the size of the population and the total processing time required. This is important because the relation between the exploration vs. exploitation (i.e., exploration/exploitation) could be decreased dramatically and the real improvements in the objective function related to the minimization of power loss could be negligible [32]. The stop criterion based on the 50 iterations without improvements corresponds to 25% of the total iterations, which was considered as an adequate criterion for finishing the exploitation process, as it allows the determination of whether the solution is a local or a global optimum. The selection of this percentage was made by using trial and error, which can take different values depending on the problem under analysis [62,63]. The parameters used for the CGA and the BH were taken from the values reported by the authors in the literature [43] in order to solve the problem addressed here. These algorithms were implemented in MATLAB by using the power flow method for DC grids reported in [14], which enabled the identification of the best size of the DGs. In order to validate the proposed PPBIL-PSO method, the four additional algorithms were combined with both PPBIL and PSO algorithms to form eight additional techniques, which can be used to evaluate the performance of the proposed solution. The validation of the effectiveness and robustness of the hybrid methodologies was carried out by evaluating the same test cases with all optimization methods. The simulation results obtained are shown in Tables 5-7, which report from left to right: the hybrid methodology, the location and size of the DGs, the mean power loss, the relative mean reduction of the power loss in relation with the base case, the mean square error of the voltage profile and its relative reduction in relation with the base case, the worst voltage profile and its location in the system, the maximum current on the branches and the averaged processing time. The average values were obtained from 1000 continuous executions for each solution method on each test scenario.

Simulation Results
To carry out the simulation, a Dell Precision T7600 Workstation with an Intel(R) Xeon(R) CPU ES-2670 at 2.50 GHz, 32 GB of RAM memory and 12 processing cores was used. The simulation results obtained from each test system are presented below: Figure 7 present the impact on the power losses and processing times of the different optimization methods in the case of the 10 bus test system. From Figure 7a, it is observed that the proposed method provides the best solution with a reduction of 66.21% of the power loss, thereby exhibiting an additional average reduction of 1.68% over the other techniques. In contrast, the PMC-BH solution produces the minimum reduction of active power loss (60.47%), thereby exhibiting the worst performance. Moreover, from column 6 of Table 5, it is observed that all hybrid solutions satisfy the maximum allowable current in the branches of this test system. If we observe the results obtained by the solution methods with respect to the objective function, we can appreciate that the methods that employ the PPBIL in the master stage achieve the best results in general terms, while the methods that use the PMC and the GA present similar results. This demonstrates the effectiveness of the solution methods in solving the problem addressed here based on the estimation of distribution algorithms in small electrical networks. Furthermore, the standard deviation (STD) of the mean power loss obtained by each hybrid solution is presented in Figure 8 where the PPBIL/PSO presents the minimum STD (0.26%), which proves the precision and repeatability of the solution provided by this method. The PPBIL/GA method also provides a low STD (0.29%), while the other hybrid methods exhibit high STD values.  Figure 7b shows the processing time required by all optimization techniques, and the proposed solution exhibits the second lowest processing time, which is 33.71% longer than the time required by the faster PMC-BH solution. However, though the PMC-BH technique is faster, it produces the worst solution in terms of power loss. Finally, after excluding the PMC-BH method, the proposed PBIL-PSO method exhibits a processing time that 53.61% lower on average in comparison with the other optimization techniques. The processing times required by the different methods in the 10 bus test system show that the methods that work with the GA and CGA in the master and slave stages require long calculation times, which can be attributed to the evolution procedures used by the genetic algorithms when they explore the solution space.

Table 5 and
Furthermore, in order to analyze the impact of each optimization technique on the voltage profile, the square voltage error (SVE) is calculated with respect to the voltage base (V base ) by using expression (15) [28]. Figure 7c shows the reduction of the square error voltage profiles in percentage where the maximum reduction is obtained by the PPBIL-PSO (67.85%) and the minimum reduction is provided by the GA-BH (63.82%). Therefore, the proposed method provides a better voltage profile in comparison with both the base case and the other techniques. In this way, Figure 7d depicts the worst voltage nodes given by each solution, which exhibit a deviation from the nominal voltage (1 p.u.) that is lower than 2%. Such a deviation according to the load and type of network is acceptable; hence, it guarantees a secure and reliable operation of the grid [9,32,64,65]. In any case, the best operation condition in terms of the voltage profiles is provided by the PPBIL-PSO, as the worst voltage node is closer to 1 p.u., in comparison with the other solutions. In contrast, the worst operation condition is provided by the PMC-BH technique.

21 Bus Test System Results
The results obtained with the 21 bus system are presented in Table 6 and illustrated in Figure 9. Figure 9a shows that the higher reduction in power loss is achieved once again by the PPBIL-PSO (78.37%), while and the lower impact is created by the PMC-BH (52.36%). In comparison with the other techniques, the proposed method provides an averaged reduction of 11.76% in the power loss. In this test system, the results obtained show that these methods help the PPBIL improve the impact on the objective function, while the other methods (PMC and GA) present a negative impact in general terms with respect to the results obtained in the 10 bus test system. This demonstrates that these methods can be stuck in local optimal solutions when the dimension of the solution space increases. In addition, we can observe a tendency to find solutions with bad quality when the BH is used in the master slave. Figure 10 presents the STD of the power loss obtained in the case of the 21 bus test system. As in the previous case, the PPBIL/PSO method provides the lower STD value (1.42%), followed by the PPBIL/GA (2.38%), and the other methods exhibit higher STD values.   With respect to the processing time, the proposed method takes second place, requiring 39.58% more time than the PMC-BH method. In any case, the proposed PPBIL-PSO is 50.35% faster on average than the other seven solutions. When all the processing times presented in Figure 9b are analyzed, we can appreciate that they increase with respect to the 10 bus test system, as these methods work with the GA and the CGA solution methodologies to exhibit the highest increment. Figure 9c,d analyze the square error voltage profile where the PPBIL-PSO method provides the best performance with a reduction of 87.61% with respect to the base case, and with an average reduction of 13.80% in comparison with the other techniques. As in the previous case, the solution with the worst reduction in the square voltage error is the PMC-BH (51.22%). Then, Figure 9d reports that the method with the best nodal voltage is the PPBIL-PSO (0.9760 p.u at bus 9), while the worst voltage is obtained by the PMC-BH (0.9383 p.u at bus 17). In addition, the figure also shows that all the methods provide an operation of around ±10% of the nominal voltage, which ensures a safe grid operation [9,64]. Finally, column 6 of Table 6 reports the maximum branch currents obtained by each solution method, which fulfill the maximum allowable current in the 21 bus test system.

69 Bus Test System Results
The results obtained with the 69 bus system are reported in Table 7 and Figure 11. Figure 11a shows the impact of the nine methods on the 69 bus system where the proposed PPBIL-PSO again provides a higher reduction in power loss with respect to the base case (90.99%) by providing an additional average reduction of 24.83% in power loss with respect to the other approaches. By analyzing the results obtained regarding the reduction of the total power loss, it can be stated that the methodologies based on the PPBIL and GA achieve good results, while the solution methods that use the PMC remain trapped in local optimal solutions; this situation can occur because this method is based in a random search and a stochastic procedure [33]. In addition, we can appreciate that the effectiveness of the sizing methods are limited to the performance of the location method while addressing the problem of the optimal integration of DGs, which implies that it is vital to use an adequate relation between the master and slave stages. The aforementioned results prove that in the case of both small and large systems, the proposed method provides the best location and size for the DGs. The high difference between the worst and the best solution presented in Figure 11a is due to the hybrid method based on the PMC/CGA, which becomes trapped in a local optimum; the same problem occurs when the PMC/CGA, GA/PSO, GA/CGA and GA/BH solution methods are used in the 69 bus test system. Therefore, those results show that the aforementioned methods are not adequate for solving the problem of optimal location and sizing of the distributed generation in large DC grids. This is also explained by Figure 12, which demonstrates that the methods that use in the master stage that use the GA and the PMC converge to regular and bad solutions, respectively. However, the PMC/CGA method exhibits a high standard deviation of the power loss. All of these indicate that many of the 1000 executions of the aforementioned algorithms are trapped in different local optimums. To illustrate the repeatability and precision of the solution methods in the 69 bus test system, Figure 12 presents the STD values of the power loss exhibited by these. This picture shows that the PBIL/PSO once again achieved the lowest STD value (5.25%), followed by the PPBIL/GA (5.37%) and PPBIL/BH (5.98%) methods, while the other ones exhibited higher STD values. Finally, based on the previous STD results, it can be concluded that the PBIL/PSO method provides the highest precision and repeatability of the solution independent of the size of the electrical network under study.    The effectiveness of the proposed method in terms of speed is shown in Figure 11b, which reports that the PPBIL-PSO is in third position, while the method PPBIL-BH (39.79% faster) is faster and the PMC-BH (39.23% faster) is in the second position. Finally, PPBIL-PSO is 52.46% faster than the other six methods. In addition, this picture demonstrates that all methods increase the processing times in an exponential form, as they use the GA and CGA methods to present the highest increment when the size of the electrical network increases.
The impact of the nine methods on the voltage profiles is illustrated in both Figure 11c,d. As in the previous cases, the proposed PPBIL-PSO provides the higher reduction in the square error voltages with a reduction of 91.91%. Moreover, the worst voltage profile provided by this method is 0.984 p.u at bus 64. Instead, the worst solution in terms of voltage profiles was obtained by the PMC-CGA with a reduction of 28.68% in the square error voltage, followed closely by the PMC-BH with a reduction of 34.27%. Both of these methods also provide the worst node voltages: 0.9333 p.u for PMC-CGA and 0.9339 p.u for PMC-BH. Finally, all the solution methods fulfill the maximum allowable current in the branches of the 69 bus test system, which is reported in column 6 of Table 7. Figure 13 shows the impact of the different hybrid methods used to find the optimal location and size of the DGs. Such a figure enables us to observe the trade-offs presented by all the hybrid methods in terms of power loss and processing time for any network size. The vertical axis presents the average value of power loss [%] with respect to the base cases for all test systems, while the horizontal axis shows the average processing time used by each method with respect to the GA/CGA in the form of percentage. This is the hybrid method that requires the longest processing time in all the test systems. This figure shows that the PMC/BH is the faster method with an average processing time of 9.08% but with the worst average power loss (55.64% of power loss in the systems without any DG). Similarly, the worst solution in terms of processing time is provided by the GA/CGA method, which has an average processing time of 791.73 s. The PPBIL/PSO is the second fastest method with an average value of 14.46%, but it obtained the best results in terms of power loss with an average value 21.50%. In this figure, the trade-off between power loss and processing time corresponds to the distance from the origin (0,0), which represents the ideal solution for the problem addressed in this study: the power loss and the processing time is equal to zero. The figure reports that the proposed PPBIL/PSO solution provides the best trade-off between power loss and speed, as it exhibits a shorter distance from the origin. The previous result demonstrated that the proposed solution provides the lowest power loss with an acceptable processing time.

Conclusions
In this paper, a hybrid method, formed with the PPBIL and the PSO algorithms was proposed for the optimal integration of DGs in DC grids. The method uses a mathematical formulation to analyze the impact of the distributed generation on the grid. In particular, the objective function is based on the reduction of the power loss while accounting for the set of restrictions required in this type of system. Finally, the performance of the proposed solution was contrasted with the eight other options, based on the optimization algorithms widely adopted in the literature.
The proposed PPBIL-PSO hybrid method provides the location and size of the DGs that produce the highest reduction in the power loss for small, medium and large DC grids. In terms of location, the GA provides similar results, but it requires a higher processing time. With respect to the sizing problem, the fastest method is the BH, but it provides solutions that are much less efficient, i.e., it leads to higher power loss than the PSO. Moreover, the processing time and power loss of all the methods were normalized for the three test cases, which demonstrated that the proposed PPBIL-PSO solution obtains the best trade-off between speed and power loss for any grid size. In addition, the PBIL/PSO exhibits the lower standard deviation for power loss in all the test scenarios; hence, the proposed PBIL/PSO provides the highest precision and repeatability for the solution when it is compared with the other hybrid approaches.
Future improvements to this study could consider the inclusion of renewable generation curves and hourly power demand into the mathematical model, which could provide a more accurate solution for a particular geographical region. Moreover, the solution can be further developed to integrate energy storage devices into the electrical network, which could be used for improving the electrical behavior and profitability of the grid. In addition, it is possible to propose an optimal tuning strategy for the parametrization of the continuous optimization techniques with the aim of balancing the exploration and exploitation of the solution space, which will enable the identification of the best compromise between the number of iterations and the population sizes. Funding: The development of this work was supported by the Instituto Tecnológico Metropolitano, Universidad Tecnologica de Bolívar, Universidad Nacional de Colombia and Minciencias (Fondo nacional de financiamiento para ciencia, la tecnología y la innovación Francisco José de Caldas) under the National Scholarship Program (call for applications 727-2015) and the projects P17211, PS2020002 and "Estrategia de transformación del sector energético Colombiano en el horizonte de 2030-Energética 2030"-"Generación distribuida de energía eléctrica en Colombia a partir de energía solar y eólica" (Code: 58838, Hermes: 38945).