Optimal Allocation and Sizing of PV Generation Units in Distribution Networks via the Generalized Normal Distribution Optimization Approach

: The problem of optimal siting and dimensioning of photovoltaic (PV) generators in medium-voltage distribution networks is addressed in this research from the perspective of combinatorial optimization. The exact mixed-integer programming (MINLP) model is solved using a master–slave (MS) optimization approach. In the master stage, the generalized normal distribution optimization (GNDO) with a discrete–continuous codiﬁcation is used to represent the locations and sizes of the PV generators. In the slave stage, the generalization of the backward/forward power method, known as the successive approximation power ﬂow method, is adopted. Numerical simulations in the IEEE 33-bus and 69-bus systems demonstrated that the GNDO approach is the most efﬁcient method for solving the exact MINLP model, as it obtained better results than the genetic algorithm, vortex-search algorithm, Newton-metaheuristic optimizer, and exact solution using the General Algebraic Modeling System (GAMS) software with the BONMIN solver. Simulations showed that, on average, the proposed MS optimizer reduced the total annual operative costs by approximately 27% for both test feeders when compared with the reference case. In addition, variations in renewable generation availability showed that from 30% ahead, positive reductions with respect to the reference case were obtained.


Introduction
Medium-voltage distribution networks cover big areas, with hundreds of kilometers in both urban and rural zones, to attend to the electricity end-users [1]. A distribution grid essentially interfaces the power systems in transformation nodes (i.e., substations) with consumers at medium-and low-voltage magnitudes to provide a quality, reliable, and secure service [2]. The main feature of an electrical network is mainly associated with its topology, as these grids are normally structured in radial form to simplify the protection schemes and reduce investment costs in all infrastructures related to conductors [3]. The radial characteristics of these networks result in considerable energy losses when compared with the transmission/subtransmission grids [4].
On the other hand, the harmful effects of global warming have significantly increased the interest in the reduction of greenhouse gas emissions from fossil fuels [5], that is combustion of coal, diesel, natural gas, and in general, petroleum-derived products [6]. Renewable energy resources are essential elements in the evolution of the electrical sector to achieve this goal, as this is the third emitter of greenhouse gas emissions after beef and pork production [7] and transportation systems [8], respectively. In the context of renewable generation, wind turbines and solar photovoltaic (PV) sources have been widely developed in recent decades, which have allowed essential reductions in the investment and maintenance costs of small-scale installations, making these technologies suitable for installation in medium-and low-voltage distribution grids [9].
For tropical countries, that is countries located near the Equatorial Line (as in the case of Colombia), the most suitable technology to produce clean energy and replace thermal generation from fossil fuels is PV generation since solar radiation suffers minor variations [10]. Instead, wind power is only suitable in coastal areas [11]; thus, large continental territories are relegated to other types of renewable energy resources, that is mainly PV sources [12]. An important aspect that must be considered in the optimal integration of PV generation units is the high variability of the energy production due to the cloud-induced fluctuations, which can produce important changes in the total power generation in intervals on the order of seconds [13]. A complete analysis regarding these fast energy generation fluctuations (due to partially shaded conditions) and their forms was reported in [14].
On the other hand, the low cost and high expansion of PV sources along the distribution grids makes the optimal design and inclusion of this technology not an easy task. This is due to bad planning, which can cause over-voltages and over-currents in nodes and distribution lines, energy loss increments, misoperation of protective devices, and deterioration of the quality of service, among other problems [15] (a complete review regarding problems derived from renewable energies in power systems can be found in [16]). To avoid those problems, distribution companies need to efficiently plan these grids, which requires an optimal integration of PV generation sources. Such a process must be performed by solving the mixed-integer nonlinear programming (MINLP) model that represents the problem of the optimal siting and sizing of the generation sources in distribution grids [17].
The current literature has addressed the problem of the optimally sizing and integration of renewable energy resources from two different perspectives. The first approach considers only the technical improvement of the grid, that is power loss reduction [18], voltage profile improvement [19], and voltage stability improvement [20]; however, the main problem of those approaches is related to the economic feasibility of those projects. The second approach deals with the optimal integration of renewable generation technologies, which considers a planning period based on an economic indicator as the objective function. This approach has the main advantage of ensuring that the final solution is technically and economically feasible [17].
Some approaches, recently reported in the scientific literature, consider economical objective functions in the problem of the optimal siting and dimensioning of renewable energy resources. For example, the application of the Chu and Beasley genetic algorithm (CBGA) using discrete-continuous codification was proposed in [21], which addresses the problem of the optimal siting and dimensioning of PV sources in medium-voltage distribution grids. Computational validations in test feeders formed by 33 and 69 buses demonstrated the efficiency when numerical results were compared with the exact solution of the MINLP model in the general algebraic modeling system (GAMS) and the BONMIN solver. Valencia et al. [9] proposed a linear approximated model to size and site batteries and renewable energy resources in distribution networks. The optimal location of those devices (integer problem) was left to a classical simulated annealing algorithm, and the operation problem (sizing problem) was solved using a linear programming model. Test feeders from 11 to 230 nodes were used to validate the effectiveness of the proposed model; nevertheless, no comparisons with other optimization methods were provided, which is the main problem of that work since it is impossible to ensure that the solution provides global optimization. In [22], the application of a new metaheuristic optimizer called the Newtonmetaheuristic algorithm (NMA) for the same problem was presented. Numerical results in test feeders with 34 and 85 nodes demonstrated the effectiveness of this optimization method through comparisons with the GAMS and CBGA, respectively. In [23], the authors reported the application of the particle swarm optimizer to locate and size PV sources and energy storage systems simultaneously. The main contribution of this research was the complete economic analysis made by the authors, which included the installation, maintenance, and operating costs of the devices. However, the authors simplified the optimization problem by considering a single nodal model for the distribution grid, which reduced the exact MINLP model to an operative nonlinear programming model. Thus, the component associated with the optimal location of the renewables and batteries remains unsolved. Cortes-Caicedo et al. [17] presented a discrete-continuous version of the vortexsearch algorithm (DCVSA) for the PV location and sizing problem. Simulations in the IEEE 33-bus and IEEE 69-bus grids demonstrated the effectiveness of the proposed approach through comparisons with the exact solution in the GAMS software and the application of the CBGA, respectively. A two-stage methodology based on the combination of a mixedinteger quadratic convex model to decide the location of the PV sources and the optimal power flow (PF) solution via the interior point method to determine the PV sizes was proposed in [12]. Numerical results obtained in this work demonstrated that the method reached similar results to the CBGA and the NMA approaches for the IEEE 33-bus and IEEE 85-bus grids.
In addition, optimization methodologies for the locating and sizing of PV generation units in distribution networks, considering technical or economic objective functions, include the Jaya optimization algorithm [24], the heuristic-based local search algorithm [25], the modified gradient-based metaheuristic optimizer [26], the mixed-integer linear approximation [27], the multi-criteria decision system [28], and recursive simulations using multiple PF evaluations [29], among other methods. The main characteristic of those approaches is that they overcome the problem of location from the problem of sizing. The former problem is solved with sensitivity analyses or heuristic algorithms, and the latter problem is solved using multiple PF evaluations.
Considering the previous revision of the literature, the main contribution of the article is the following: the application of a recently developed generalized normal distribution optimizer (i.e., GNDO), with a discrete-continuous codification, to solve the problem of the PV location (selection of nodes) and optimal sizing in the master stage. This solution allows transforming the MINLP model into a simple PF problem for distribution networks, which is solved using the successive approximation power flow (SAPF) method in the slave stage. Simulations in the IEEE 33-bus and IEEE 69-bus grids confirmed the effectiveness of the proposed optimization method, since the final objective function values were better than those provided by current literature approaches. In this way, the satisfactory results reported by the vortex-search algorithm in [17] were improved by approximately USD 89.95 and USD 341.18 per year of operation, respectively.
It is important to remark that the GNDO has not been previously applied to problems of distribution system planning; thus, it is as an opportunity for research addressed in this work. The GNDO algorithm was proposed in 2020 by Zhang et al. [30] to extract the parameters of PV modules with different numbers of diodes. However, that is an optimization problem defined in the domain of the real variables (i.e., continuous variables). Therefore, the work reported in this paper adapted the GNDO method to deal with a mixed codification, which combines integer variables related to the buses where the PV sources will be assigned and continuous variables describing their optimal sizes. The rest of the paper is organized as follows: Section 2 shows the exact MINLP model for the studied problem using a mathematical representation in the complex domain. Then, Section 3 describes the main features of the proposed MS optimization approach, which combines the GNDO in the master stage with the PF method based on the successive approximations in the slave stage. Section 4 presents the test feeder characteristics, including their peak load behavior, demand, and generation curves; moreover, this section also describes the parametrization of the objective function. Section 5 presents the complete analysis and discussion of the results obtained by the proposed MS optimizer and provides comparisons with approaches recently reported in the scientific literature. Finally, Section 6 lists the main conclusions of this work.

Mathematical Formulation
The optimal allocation and sizing of PV generation units in distribution networks correspond to a mixed-integer nonlinear programming (MINLP) model. The integer component of the optimization model is associated with selecting the nodes where the PV generation units are located, which is indeed a binary decision. The remaining variables are continuous, and those are related to the optimal sizes of the PV generation units, the solution of which is equivalent to the PF problem (i.e., voltage, power, and current calculations) [31]. The MINLP model, which represents the problem of the optimal allocation and sizing of PV generation units in distribution grids, is presented below. This work adopted the complex domain representation to reduce the extension of the optimization model with respect to the classical rectangular representation.

Objective Function Representation
The studied problem is a time-dependent optimization problem, which includes the demand and generation curves. Those curves help to improve the estimation of the final grid operating costs with respect to approaches based only on peak demand. The main interest in the integration of PV sources in electrical distribution grids is associated with the minimization of the annualized operative costs of the grid, which includes the investment and maintenance costs of the PV generators, and the energy acquisition costs in the substation node. This objective function is described by Equations (1)- (3). (2) being: where the value of the objective function is A cost . Such a value is related to the annual grid operative costs, where f 1 determines the annualized energy purchasing costs in the substation node, and f 2 measures the total investment and maintenance costs of the PV generation sources installed along with the distribution network. C kWh is the expected energy purchasing cost at the terminals of the substation, which is taken from the spot market. The number of years is defined as T. f a is the cost annualization factor, which depends on the expected internal return rate, that is t a and N t define the length of the planning period in years. s cg k,h represents the complex power output at the terminals of the substation bus connected at node k for each period of time h. ∆ h represents the period of time in which all the electrical variables are constants, that is 1 h for a daily operative scenario. t e is the parameter that measures the expected increment of the energy costs in the substation node for the duration of the planning project (this value is provided by the distribution company). C pv is a parameter that defines the average cost of installing 1 kWp using PV generation. C O&M represents the operating and maintenance costs of the PV generation units, and s pv k,h represents the complex power output of the PV generation unit assigned to k in period h. Finally, the sets that contain all the nodes, periods of time, and years of the planning project are defined as N , H, and T , respectively.

Set of Constraints
The classical power equilibrium equations, maximum and minimum device capabilities, and voltage regulation restrictions constrain the PV location and sizing problem. Those restrictions are listed in Equations (4)- (12).
In the previous expressions, S d k,h is the power demand in the complex domain at node k during period h. v k,h and v j,h are variables in the complex domain associated with the voltages at nodes k and j for each period of time h, respectively. Y kj represents the complex admittance that relates nodes j and k. G are the complex power lower and upper generation bounds, respectively, related to the PV source installed at node k. i kj,h represents the complex current flow through the distribution lines interconnecting nodes k and j at each period h. y k represents the decision variable regarding the installation (y k = 1) or not (y k = 0) of a PV source in the bus k. v min k and v max k represent the lower and upper bounds for the voltage regulation assigned to all the distribution grid nodes at any period of time. Finally, N ava pv corresponds to the PV generation units available for installation.

Model Characterization and Interpretation
The mathematical MINLP formulation presented in Equations (1)- (12) can be interpreted as follows: the objective function is defined in Equation (1), which is formed by the total energy purchasing costs in terminals of the substation bus defined by f 1 in Equation (2) and the investment and maintenance costs of the installed PV sources defined by f 2 in Equation (3). In Equation (4), the complex power equilibrium constraint at each node of the network is defined, which is the most challenging constraint because it is a set of nonlinear non-affine constraints that require numerical methods to be solved efficiently [32]. Equation (5) shows that the generation of the complex power in the PV generation units is a factor, that is their peak sizes, multiplied by the expected generation curve in the zone of influence of the medium-voltage distribution grid. Equation (6) shows that the generation in the PV generation units has a unity power factor, which implies that the reactive power generation is set to zero for all the periods under analysis. The box-type constraint (7) determines the lower and upper apparent power generation bounds admissible for the secure operation of the conventional generation (i.e., substation bus) connected at node k for each period of time. The inequality constraint (8) determines the admissible bounds for power generation in the PV source connected at node k when the binary variable y k takes a value of 1. The box-type constraint (9) bounds the voltage magnitudes in all nodes of the network for any time [33]. The inequality constraint (10) ensures that the current variables in the distribution lines remain between their thermal limits. The inequality (11) fixes the maximum number of PV generators available for integration in the distribution grid. Finally, Constraint (12) shows the binary nature of the decision variables regarding the optimal location of the PV sources in the distribution grid.
The solution of the model (1)- (12) can be addressed using MS optimization methodologies from the family of metaheuristics. Those methodologies allow decoupling the discrete optimization problem from a continuous one. With such a decoupling, it is easy to assess the objective function with a simple PF methodology for distribution grids if a metaheuristic algorithm provides the location and sizes of the PV sources. These MS optimization approaches were applied successfully to this problem by the authors of [21] using a CBGA, the authors of [22] with the application of the Newton-metaheuristic optimizer, by the authors of [9] with the application of the simulated annealing optimization algorithm, as well as by the authors of [17] via the DCVSA.
Based on those recent works, this work proposes the application of the recently developed GNDO in the slave stage to determine the optimal location and sizes of PV sources using a discrete-continuous codification [30], while in the slave stage, the SAPF method is used to evaluate the objective function value [34]. The following sections show that this novel solution provides improved results for the IEEE 33-bus and IEEE 69-node grids in comparison with the solutions recently published in the literature.

Solution Methodology
The GNDO method is defined as the master algorithm to determine the optimal location and sizing of the PV generators in distribution grids. This is based on a discretecontinuous codification with the structure defined in Equation (13).
x t i = 2, k, · · ·, 10 |1.3846, s pv k , · · ·, 0.6259 (13) In such an expression, x t i is the ith individual in population X for the current iteration t. The main characteristic of this codification is that the size is 1 × 2N ava pv , where the first N ava pv positions are associated with the buses where the PV generators will be assigned, whereas the positions N ava pv + 1 to 2N ava pv are defined as their optimal sizes. The main characteristic of the codification given in Equation (13) is that it solves the problem of the location and sizing of the PV sources in one stage, which allows the transformation of the MINLP model defined in (1)-(12) into a nonlinear programming problem that is easily solvable with a PF method for distribution grids.

Slave Stage: SAPF Method
This PF method is a generalization of the classical backward/forward PF solution proposed by Montoya and Gil-González in [34]. This PF approach deals with the power balance equations in the complex domain using a recursive PF formula that is derivativefree. This recursive formula decouples the slack power balance equation from the demand nodes to obtain a matrix relation as a function of the admittance submatrices. The recursive PF formula is given by Equation (14).
where m is the counter of iterations; V d corresponds to the vector that contains all the complex demanded voltage variables for each period of time h; S pv,h is a vector in the complex domain that contains the power outputs in the PV generators at each period of time h; S d,h is a vector in the complex domain that contains all the active and reactive power consumptions in the demand nodes, which are represented in rectangular coordinates for each period of time h; V s,h corresponds to the complex voltage output in the substation bus at any period of time (this is a known input for PF problems [34]); Y dd contains all the admittance relations among demand nodes, which is a square invertible matrix; Y ds is a rectangular matrix that associates the substation bus with the demand nodes; finally, diag(z) is a diagonal matrix composed of the elements of the z vector, and z is the conjugate version of the z vector.
One of the main features of the SAPF method, defined by the recursive Formula (14), is that the Banach fixed-point theorem can demonstrate its convergence [35]. In this study, the maximum difference of the voltage magnitudes in the demand nodes, in two consecutive iterations, was used to determine the convergence of the recursive Formula (14) to the PF solution. This convergence condition was evaluated using Equation (15), where ζ is defined as the convergence error, which was set as suggested by authors of [34] When the PF solution is reached, the following step in the slave stage corresponds to the calculation of the power generation in the slack bus, as presented below.
Then, with the solution of Equation (16), it is possible to obtain the value of f 1 . In addition, with the codification provided in Equation (13), the value of f 2 is also obtained. This notwithstanding, to deal with the possible unfeasible solutions in the solution space, the objective function (1) is typically replaced by a fitness function [36]. In this study, the proposed fitness (A f ) is presented below: In A f , the variables θ 1 , θ 2 , θ 3 , and θ 4 represent the penalization factors of the objective function. These factors are activated in the case of the voltage regulation bounds, the thermal current limits of the conductors, or the substation active power generation bounds being violated. In this work, the magnitude of those penalty factors is defined as 100 × 10 3 , each one with the adequate units.
The most important characteristic of the fitness function in (17) is that if all the constraints of the optimization model are satisfied, i.e., Restrictions (4)- (12), then the final value of A f is equal to the original objective function A cost .

Master Stage: GNDO
The GNDO is inspired by the Gaussian distribution theory, which is a powerful tool for describing multiple natural phenomena [30]. The general Gaussian distribution takes the form presented in Equation (18), where x is a random variable that corresponds to the distribution probability with the location parameter µ and the scaling parameter δ [37,38].
From Equation (18), the main variables of a normal distribution are the location (µ) and the scaling (δ) parameters. These parameters define the mean value and standard deviation of random variables x. Figure 1 depicts the behavior of a Gaussian distribution for different values of µ and δ [30].  Behavior of a Gaussian distribution with different µ and δ parameters: (a) variation in the µ parameter by fixing the δ parameter; note that the variation of the µ parameter moves the density probability function from left to right with the same amplitude as this parameter increases; (b) variation in the δ parameter by fixing the µ parameter; observe that the increment in the δ parameter makes the density probability increase its amplitude and also reduce its bandwidth.
The behavior of the normal distribution shown in Figure 1 exhibits the following characteristics [30]: When the δ parameter remains constant (see Figure 1a), the Gaussian distribution is moved to the direction where the µ parameter increases; In the case of the µ parameter holding constant (see Figure 1b), the probability density is expanded in the direction of the δ parameter's change.
Then, the previous characteristics of the Gaussian distribution function in (18) were used in [30] to propose an efficient combinatorial optimization method.
Moreover, it is widely known that metaheuristic optimizers based on populations have three stages in their implementation [39]: The initial population is generated with a normal distribution scattered along the solution space. This population evolves by exploring and exploiting the solution space, searching for the global optimum guided by specific movement rules. At the beginning of the searching process, the variance of the positions of all the solution individuals is significant, and the random position of the decision variables around the global optimum can be considered as random variables subject to a normal distribution behavior; ii. Then, the distance between the mean position and the global optimum is continuously reduced, and the variance of the positions among all individuals is also gradually decreased; iii. Finally, the distance between the mean position and the global optimum, as well as the variance of the positions among individuals reach a minimum value.
Based on the previous characteristics of the general evolution process of a populationbased optimizer, the general characteristics of the GNDO approach are described below [30].

Local Exploration
The local exploration process involves finding better optimal solutions around the current locations of all solution individuals. Considering the relationship between all the solution individuals in the current population and the normal distribution, a generalized model can be constructed as follows: where v t i is the trail vector of the ith individual in the current iteration t, µ i represents the generalized mean position of the ith individual, δ i corresponds to the generalized standard variance, and η is a penalty factor. N i represents the number of individuals in the population, and each solution individual has the general structure presented in Equation (13). Then, parameters µ i , δ i , and η are calculated using Equations (20) and (22).
In the previous equations, λ 1 , λ 2 , a, and b are numbers between 0 and 1 with a uniform distribution, x t best is the best solution found so far, and M is the vector with the mean position of the current individuals in the population. Finally, M can be obtained as defined in (23).

Global Exploration
Global exploration in metaheuristics attempts to determine the location of the promising solution regions through the solution space [40]. The authors of [30] proposed a global exploration using three main components, as presented in Equation (24).
where β × (|λ 3 | × v 1 ) contains information about the local solution region and (1 − β) × (|λ 4 | × v 2 ) shares the global information of the solution space. In addition, λ 3 and λ 4 are random numbers with a uniform distribution, β represents an adjusting parameter that is selected randomly between 0 and 1, and v 1 and v 2 correspond to two trail vectors.
In (25) and (26), j, k, and m are three integer numbers between 1 and N i , that is three individuals in the current population, which must be different among each other and also different from the current ith individual.
Prior to deciding whether the potential solution v t i will make part of the next population, each gen in it is revised to ensure the feasibility of the solution space, that is, where x min j and x max j are the minimum and maximum admissible bounds of the decision variable j, respectively. In addition, because of the integer component of the codification, the first N ava pv is rounded to its nearest integer values. Finally, the selection of the next individual is made using the following rule:

General MS Implementation
The application of the GNDO method combined with the SAPF method for the problem under study is presented in Algorithm 1.

Algorithm 1: General implementation of the proposed MS optimizer.
Data: Choose the AC network under study Find the per-unit network equivalent; Select the values of N i and t max ; Select the µ parameter as 1 2 and t = 0; Generate the initial population with the structure of (13); Evaluate the fitness function value (slave stage) for each individual solution x t i , that is A f x t i ; Select x t best as the individual with the minimum value of the fitness function; for t ≤ t max do for i = 1 : N i do Generate a random value γ between 0 and 1; if γ > 1 2 then /*Local exploration search*/; Select the best current solution x t best ; Compute the value of the vector M using Equation (23); Calculate the generalized mean value µ i using Equation (20); Obtain the generalized standard variance δ i through Equation (21); Compute the penalty factor η with Equation (22); Make the local exploration using Equation (19)

First Test Feeder
This is a grid with 32 lines and 33 nodes operated at medium-voltage values with a magnitude of 12.66 kV in the slack source. The grid topology is shown in Figure 2. During the peak load condition, this system absorbed 3715 + j2300 kVA, and with this load, the total power loss increased to 210.9876 kW. The electrical data for this system are reported in Table 1.

Second Test Feeder
This is a grid with 68 lines and 69 nodes operated at medium-voltage values with a magnitude of 12.66 kV in the slack bus. The grid topology is shown in Figure 3. During the peak load condition, this system absorbed 3890.7 + j2693.6 kVA, and with this load, the total power loss increased to 225.0718 kW. The electrical parameters for this system are reported in Table 2.

Objective Function Evaluation
To obtain the objective function values, the information in Table 3 was used.

3
- In order to emulate the daily generation profile of the PV units and the load profile behavior, typical generation and demand curves in Medellín city (Colombia) were used, which are both depicted in Figure 4. This information was provided by Grisales et al. in [41].

Computational Validation
The proposed GNDO was implemented in the MATLAB programming environment to solve the problem under study. For comparison purposes, the following methodologies were considered: the BONMIN solver in the GAMS environment (solution of the exact MINLP formulation), the discrete-continuous version of the CBGA (i.e., DC-CBGA) [21], the Newton-metaheuristic algorithm (i.e., DCNMA) [22], and the vortexsearch algorithm (i.e., DCVSA) [17]. Moreover, to evaluate the combinatorial optimizers, a population of 10 individuals and 100 consecutive repetitions of each algorithm, with 1000 iterations at each running, were considered as the settings. Note that the proposed GNDO works as a comparative methodology with the discrete-continuous codification defined in Equation (13); for this reason, it was renamed as DCGNDO. Finally, for both test feeders, the possibility of installing three PV generation units was considered, where each one of them had a maximum size of 2400 kW.

First Test Feeder
After the execution of the proposed and comparative optimization methods, the optimal solutions reported in Table 4 were obtained. The numerical information in Table 4 reveals that: (i) all the metaheuristic optimizers improved the solution reached by the BONMIN solver, which confirmed that the MINLP model is non-convex in nature; moreover, the presence of binary variables made conven-tional optimization methods stuck in locally optimal solutions; (ii) the proposed DCGNDO method reached the best optimal solution for the first test feeder with a reduction of US-D/year 1,000,783.62, which improved the solution reported in [17] for the DCVSA by about USD/year 89.95; (iii) all the optimization methods allowed reaching reductions between 26.99% in the case of the BONMIN solver and 27.04% in the case of the proposed DCGNDO, with respect to the benchmark case.
An additional important characteristic of the optimal solutions for the proposed and comparative methods, reported in Table 4, is that all the solutions were near certain regions in the distribution network. Those areas were between Nodes 8 and 11, 14 and 18, and 30 and 31, respectively. Those regions will be useful for the distribution company to have multiple possibilities, regarding the final installation of the PV sources, to ensure reductions with respect to the benchmark case above 26.99%.
The average performance of each metaheuristic optimizer after 100 consecutive repetitions is reported in Figure 5, which shows the average reduction with respect to the benchmark case. Reduction (%) The most promising result, shown in Figure 5, is that the proposed DCGNDO had the best average performance with respect to the mean objective function reduction, i.e., 27.01%, which implies that after each 1 of the 100 executions, the expected reduction in the objective function value would be upper higher than 27%.
To demonstrate that the proposed DCGNDO ensures that the voltage profiles are between the minimum and maximum voltage regulation bounds, that is ±10%, the maximum and minimum voltage magnitudes for the first test feeder during its daily operation are depicted in Figure 6.  The voltage behavior shown in Figure 6 depicts that: (i) the maximum voltage magnitude in the first test feeder was presented in Time Period 14, i.e., when the PV generation was maximum, as depicted in Figure 4, with a value of 1.0321 pu; (ii) the minimum voltage magnitude coincided with the benchmark case during Time Periods 20 and 21, with a magnitude of 0.9038 pu, which was expected since during those periods, the demand is maximum (i.e., peak load condition) and the PV generation is null.
However, to confirm that the active power generation in the substation bus is always positive or zero, Figure 7 depicts the output active power in the slack terminals for both the benchmark case and the solution provided by the proposed DCGNDO. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23  The output power behavior in the substation bus, presented in Figure 7, confirmed that when the PV generation starts to increase, the slack active power generation begins to decrease until a value of zero is reached, which corresponds to the maximum PV generation (Time Period 14). In addition, as previously demonstrated in the literature, the slack power generation experienced a duck curve behavior, which is natural when large amounts of PV generation are integrated into the distribution grid [42]; however, the slack generation curve in Figure 7 and the voltage profile behaviors in Figure 6 demonstrate the effectiveness of the fitness function proposed in Equation (17) to ensure the feasibility of the final solution obtained by the proposed DCGNDO approach.

Second Test Feeder
In the second test feeder, it is essential to highlight that the BONMIN solver did not converge to any feasible solution due to the increment in the solution space compared to the first test feeder. For this reason, Table 5 only presents the numerical results obtained by the comparative metaheuristics, as well as the proposed optimization method.
The numerical information reported in Table 5 reveals that: (i) the proposed optimization method also reached a better solution for the second test feeder, with a reduction of 27.16% (USD/year 1,053,276.55); this result improved the solution found by the DCVSA by about USD 341.18 per year of operation; (ii) the margin of the objective function improvement with respect to the benchmark case was between 27.12% for the DCNMA and 27.16% for the proposed DCGNDO, respectively; (iii) the location of the PV sources was distributed in certain regions of the network, which were formed by Nodes 12-16, 22-24, and 60-64. Similar to the first case, those areas are essential because they will help the distribution company select the final implementation of the PV sources, based on additional studies such as the feasibility of the terrain or the size of the vegetation nearest to these nodes, among other factors.
The average performance of each metaheuristic method after 100 consecutive executions (global exploration of the solution space) is depicted in Figure 5, which reports the average reduction with respect to the benchmark case. The behavior of the average objective function reduction, shown in Figure 8, confirms that the proposed DCGNDO approach had the most stable performance since it reached reductions near 27.06%, being only followed by the DCVSA proposed in [17]. Reduction (%) To verify that the solution obtained with the DCGNDO fulfills the minimum and maximum voltage regulation constraints, that is ±10%, Figure 9 depicts the maximum and minimum voltage magnitudes for the second test feeder during its daily operation.  The results reported in Figure 9 show that the maximum voltage value occurred in Period 14, when the PV generation is maximum, with a magnitude of 1.0378 pu. However, when the demand is maximum (Periods 20 and 21 in Figure 4) and the PV generation is null, the minimum voltage profile fully matched the benchmark case with a magnitude of 0.9092 pu.
To ensure that the power generation in the slack source is at least zero or positive when PV generation is installed, the active power output at the substation terminals for both the benchmark case and the DCGNDO solution are compared in Figure 10. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23  The behavior of the active power output, shown in Figure 10, reveals that without the presence of PV generation in the distribution system, the active power in the substation followed the behavior of the demand load, including power losses (see Figure 4). However, when the PV generation units were installed, the power output in the substation decreased as the PV generation increased the power injection into the grid, being zero when the PV generation is maximum at Time Period 14.
As in the previous test, the power behavior shown in Figure 10 in conjunction with the voltage profiles depicted in Figure 9 confirm that the fitness function ensured the practical feasibility of the solution space since the final solution was indeed between the voltage regulation bounds. Finally, the active power was always positive or zero for both test feeders.

Complementary Analysis
The effect of the final objective function value when the PV generation varies from 30% to 100% of its power availability is presented in Figure 11. Such results confirm that the optimal solution was found by the proposed DCGNDO approach under possible generation variations along the planning project.  Figure 11. Variations of the PV generation availability from 30% to 100% during the planning period. Note that the shaded area is a part of the possible uncertainties of the generation availability during the planning horizon.
Considering the variations in the PV generation availability for both test feeders, expected reductions reached by the solution of the DCGNDO approach with respect to the benchmark case are presented in Figure 12. 30   The effect of the renewable generation availability variation, shown in Figure 12, shows that: (i) the minimum reduction with the availability of 30% in the PV power output allowed reducing about 0.42% and 0.54% of the total grid operative costs for the IEEE 33-bus and IEEE 69-bus systems, i.e., USD/year 15,265.74 and USD/year 21,028.33, respectively; (ii) when the availability of the PV generation varied from 50% to 80%, the expected reductions in the first test feeder changed from 8.33% to 19.72% and for the second test feeder between 8.47% and 19.86%, respectively; these values imply that the expected reductions were higher than 8.30% with a variable gain of approximately 11.40% in this range of renewable generation variations; (iii) a the expected additional gain when the PV generation increased by 10%, its power availability was approximately 3.80% for both test feeders, which implies that for the first test feeder, there was an additional annual gain of USD/year 140,617.30 and for the second test feeder of USD/year 147371.60, respectively.
It is important to remark that the proposed DCGNDO required 268.49 s for the first test feeder and 1237.23 s for the second test feeder. These values confirm the effectiveness of the proposed optimization approach for solving the problem, where the solution space regarding the nodal location of the PV sources has 4960 possibilities for the first test feeder and 50,116 combinations for the second test feeder. However, there are infinite possible sizes between 0 kW and 2400 kW regarding the dimensions of the PV generation units for each one of these nodal combinations. In addition, the time required to obtain the solution was less than 4.47 min and 20.62 min in both test feeders, respectively. These values imply that the distribution company can perform multiple simulations to investigate the best PV location and size alternatives prior to the final physical implementation in the selected nodes.
Note that the proposed economic analysis is subject to improvements by considering the high variation in the power generation of the PV sources induced by partial shading due to clouds [13]. In this sense, a detailed analysis in the time scales of seconds or minutes is suggested as potential improvement for the proposed methodology. In addition, it is also important to highlight that the proposed optimization methodology does not depend on the size of the distribution system, i.e., the codification retains the same structure, and the evolution rules of the DCGNDO do not change. However, the number of nodes of the distribution grid affects the processing times required by the slave stage for the solution of the PF problem, which increases the total computational time to reach the optimal solution for the problem.

Conclusions and Future Works
The problem of the optimal placement and sizing of PV generators in radial distribution networks was addressed in this research by applying an MS optimization approach. In the master stage, applying the GNDO approach with a discrete-continuous codification solved the location and sizing problems with a unified codification. An SAPF method was applied to calculate the fitness function value. Numerical comparisons with metaheuristic optimizers applied in the test feeders provided the following conclusions: The expected reductions in the annual operative costs were 27.04% and 27.16% for both test feeders. These values imply annual reductions in the operation costs of USD 1,000,783.62 and USD 1,053,276.55 per year of operation in each test feeder, respectively; The average behavior of the proposed DCGNDO after 100 repetitions showed that the expected reductions in the power losses, for both test feeders, was higher than 27%, which implies that after each execution, the reductions in the objective function would be higher than USD/year 999,122.95 and USD/year 1,047,113.9811 for the test feeders, respectively; When the renewable generation varies from the expected output, between 30% and 100%, we note the following: If during all the 20 years, the PV generation is only 30% of the expected value, then both test feeders experience positive reductions of the total annual operative costs of 0.42% and 0.54%, respectively. However, if the PV availability exceeds 60% of the nominal curve, then in the first test feeder, the expected annual grid operative costs would be reduced to 12.18% or higher. For the second test feeder, this reduction would be 12.33% or higher. Funding: This work was supported by the Universidad Nacional de Colombia under the project "Reconfiguración de paneles fotovoltaicos para la maximización de extracción de potencia bajo condiciones de sombreado basado en procesamiento de imágenes con inteligencia artificial" (Hermes 51503).

Informed Consent Statement: Not applicable.
Data Availability Statement: Data is contained within the article.