Coyote Optimization Algorithm-Based Approach for Strategic Planning of Photovoltaic Distributed Generation

The optimal planning for distributed generations (DGs) associated with photovoltaics (PVs) in the utility-owned distribution system is crucial for increasing high penetration of renewables while against practical system operation constraints. Such PV-DG planning is categorized as a complicated mixed-integer nonlinear programming (MINLP) problem and is extremely difficult to solve by using conventional methods. In recent years, several bio-inspired metaheuristic algorithms have been proposed to tackle various complicated real-parameter optimization problems. This paper proposes a two-stage approach including a new bio-inspired algorithm, Coyote Optimization Algorithm (COA), to solve the large-scale MINLP PV-DG sizing problem considering different load levels. The objective function terms under consideration include the total system power loss and voltage regulator tap changes at different load levels while against limits of rms bus voltages, tap changes, and PV-DG constraints at each candidate bus. The proposed method is tested using the IEEE 123-bus unbalanced benchmark system and an actual utility distribution network. Results obtained are then compared with those obtained by a classic MINLP solver-based and four other bio-inspired methods. Moreover, results also show that the proposed method leads to lower loss, a minimum number of regulator tap changes, and higher PV penetration capacity among the compared methods and is suitable for solving the large-scale PV-DG planning problem in distribution systems.


I. INTRODUCTION
The growing awareness of carbon emissions associated with the drastic fossil-fuel consumption in electrical energy production has led to an urgent need for mitigating global warming. Therefore, the photovoltaic (PV) solar energy has become more popular to serve as an alternate resource of fossil-fueled electricity generation [1]. Also, because of the significantly reduced cost of PV modules, MW-scale PV distributed generation (DG) installations at distribution networks are increasing worldwide. Therefore, the hosting capacity of PV generation in a distribution feeder that satisfies power system operation constraints becomes an imminent goal to promote a higher penetration of renewable energy The associate editor coordinating the review of this manuscript and approving it for publication was Tossapon Boongoen .
resources. A literature survey shows that approaches for PV-DG planning of the distribution network fall into two major categories: with and without (i.e. deterministic) including uncertainties associated with either PV generation output and/or load variations [2]- [5]. Both planning categories consider either single-or multi-objective functions while against the system and PV-DG operation constraints. The purpose of studies typically involves decisions of locations, sizes, or both locations and sizes of PV-DG units in the system [2]- [17], [19]- [25].
When considering generation output uncertainties, several methods have been presented for PV-DG planning. A review of recent publications indicates that most proposed methods assume that PV irradiance is modeled by beta probability distribution function and the load variation is modeled by Gaussian distribution function. Most of the test systems are in small or medium scale. Approaches for long-term (one or more years) planning considering uncertainties are not commonly seen because of the problem complexity and the requirement of historical data and producing scenarios for modeling purposes. The methods adopted generally are analytical [6], Monte-Carlo simulation-based [7], metaheuristic-based [8], and hybrid methods [9]- [11]. When considering uncertainties in PV planning, the number of PVs to be placed at candidate busses in the system is limited due to the enormous solution search space. A larger scale of the system with multiple PVs planning while considering load and/or PV output uncertainties remains a challenging task.
On the other hand, the deterministic approaches generally do not consider PV generation or load uncertainties. These methods include analytic [12], classical [13], metaheuristic-based [14]- [17], [19]- [22] and hybrid types [23]- [25]. The analytical type is more suitable for smaller systems with a limited number of PV-DG units in the planning. The classical type of approaches usually formulates the problems by mixed-integer nonlinear programming and solved by commercial off-the-shelf software. A metaheuristic approach is an iterative method that guides a subordinate heuristic by combining different intelligent concepts for exploring the solution search space. Such an approach is often inspired by observing the phenomena occurring in nature. The hybrid type combines two or more of the above approaches.
In this paper, the authors firstly review several metaheuristic approaches not considering uncertainties for finding optimal locations and sizes of PV-DG units while managing the system loss, PV hosting capacity, regulator tap changes, and network voltage profile. Like many other resource allocation problems, the DG planning problem formulation is classified as a mixed-integer nonlinear programming (MINLP) problem and is extremely difficult to solve by using conventional methods because of its highly non-convex, discrete, and constrained nature [13]. Over the past two decades, the bio-inspired metaheuristic methods have gained great interest in applications to the described problems. Not like commonly seen algorithms such as classic Newton-type methods that are easy to trap in the local optimum during the hill-climbing or gradient search, the bio-inspired methods adopt algorithms including certain randomness in the solution procedure and is considered as a higher-level method using specified selection mechanisms and information sharing for finding the global optimum solution. The paradigm of bio-inspired algorithms for global search and optimization can be broadly divided into three classes: swarm intelligence, evolutionary, and ecology-based classes [26]. For instance, the particle swarm optimization (PSO), ant bee colony (ABC), and grey wolf optimizer (GWO) fall into the first class. For the second class, genetic algorithm (GA) and differential evolution (DE) are commonly seen methods. The biogeography-based optimization (BBO) and invasive weed optimization (IWO) methods belong to the third class. The following gives an overview of applying bio-inspired methods to solve the DG planning problems.
In [15], the authors adopted a GA-based method to determine the location and capacity of PV-DG for the area development plan of a distribution network. The objective is to minimize power loss and improve the system voltage profile. Reference [16] presented a PSO-based method for PV-DG planning of a master-slave controlled microgrid, where the master mode is in non-unity power factor operation and the slave mode is in unity power factor operation. The proposed model is formulated as a MINLP problem and is incorporated into an optimal power flow framework considering a variable load profile. However, one of the disadvantages of the GAand PSO-based methods is that there is no guarantee of finding the global optimal solution due to the early trap in the search space. In [17], the imperialist competitive algorithm (ICA) was used to solve the DG distribution planning problem. It is noted that ICA is only applied to some of the standard optimization problems [18]. Reference [19] used the big bang-big crunch method for planning DGs. The method tries to minimize power loss in an unbalanced distribution system. An improved non-dominated sorting GA was used to solve the optimal planning of multiple DG units in [20], which is to minimize load consumption in the network and maintain the bus voltage within the acceptable range. Reference [21] presented a comprehensive teaching learning-based optimization technique for the optimal allocation of DGs in radial distribution systems to improve network loss reduction, voltage profile and annual energy savings. The proposed method possesses immunity to local extrema trappings. However, the selection of the optimal number of DGs in the distribution networks is only limited to three units. In [22], the adaptive quantum-inspired evolutionary algorithm was proposed for placing and sizing DGs and capacitors. Nevertheless, the planning study only considers a given load demand at a specific time instant and only up to three DGs can be placed in the system.
In addition to the aforementioned metaheuristic algorithms for PV-DG planning, several hybrid methods have been proposed to improve the solution. For instance, [23] proposed a GA-based Tabu search method to investigate and analyze the optimal locations of multiple types of DG units with certain capacities for optimizing net present worth subject to economic and technical constraints. The cost terms include capital, replacement, operation and maintenance, and reliability improvement costs. In [24], several performance evaluation indices such as power loss, voltage deviation, reliability, and shift factor are used to develop the multi-objective function while considering different load models. The combined GA and PSO-based solution algorithm is then applied to find the optimal sizing and placement of DGs. A hybrid grey wolf optimization (HGWO) method was proposed for optimal allocation of DGs [25]. The locations, sizes, and the total number of DGs to be placed are under consideration. VOLUME 8, 2020 This paper studies deterministic planning for PV-DG placement and sizing in distribution systems. To assess the effectiveness of the proposed planning method, the electric utility-owned PV-DG planning for annual loss reduction and load balancing while against network and PV-DG operation constraints are considered. In the study, the authors propose a two-stage planning for PV-DGs. The placement problem is firstly solved by calculating the loss reduction sensitivity factor (i.e. LRSF) without and with DG installation at each network bus. Those with top LRSF values are selected as the candidate DG busses. Then, the newly proposed Coyote Optimization Algorithm (COA) is adopted to find the optimal size of each candidate DG bus [27]. The COA is classified as both swarm intelligence and evolutionary heuristic inspired on canis latrans species. It considers social structure and experience exchange during hunting the prey while each coyote is a potential solution and its social condition is the cost of the objective function. In [27], the COA is proposed to solve small-and mid-scale unconstrained and constrained real-parameter single-objective optimization problems through testing 40 benchmark functions with 92 cases not including any engineering problems [28], [29]. The number of solution variables is only up to 100 without any integer variables in the test cases. To test the usefulness of COA for solving practical and larger scale of global constrained optimization problem including both real and integer solution variables and constraints, this paper applies COA to solve PV-DGs planning problem and finds the optimal size of each DG bus through minimizing the total real power loss and the number of tap changes of voltage regulators while the rms voltage at each bus is controlled to satisfy the system operation limits.
To evaluate the multiple objective terms, a weighted sum method is applied for determining the fitness of the multi-objective function and obtain the best solution. The weighted factor depends on the level of importance between the components of the objective function. In this study, the EPRI OpenDSS R distribution system simulation tool and Matlab R are adopted for solving power flow problems [30]. The OpenDSS is to perform sequential-time power flow simulations including PV-DGs over a long time period when the generation alters the load profile. Traditional distribution power flow solvers are formulated with a radial circuit and the forward-backward sweep ladder methods are commonly seen. The default power flow solution method is based on a fixed-point iteration method to solve a set of nonlinear equations, which is computationally efficient for sequential time solutions [31]. When a power flow analysis is completed, the power losses, bus voltages, and branch flows are calculated.
In the study, the IEEE 123-bus benchmark system and an actual 137-bus distribution network are under test [32]. Results obtained by the proposed method are also compared with a conventional MINLP method [13] and four bio-inspired methods including GA [15], PSO [33], biogeography-based optimization (BBO) [34], grey wolf optimizer (GWO) [35]. It shows that the proposed method is superior in both cost minimization and convergence.
The organization of the paper is as follows. In section II the problem formulation for placement and sizing of PV-DG units considering voltage, and tap changes constraints are introduced. Section III illustrates the proposed LRSF-based placement and COA-based sizing procedures. Section IV then reports test results and Section V provides the conclusion.

II. PROBLEM FORMULATION
The problem of optimal placement and sizing of PV-DGs considering multiple objective functions is challenging due to its highly non-convex nature. In the study, the objective is to minimize the total power loss and the tap changes of the voltage regulators while maintaining rms voltage at each network bus, regulator tap positions, and PV-DG capacity and power factor constraints. Listed below describes the problem formulation for the optimal planning of PV-DGs in a distribution network.

A. OBJECTIVE FUNCTION
The objective function to be minimized includes two components: system power loss reduction rate and the number of tap changes of the voltage regulators.

1) TOTAL POWER LOSS
Real power loss is an important index for the economic and technical assessment of PV-DG placements. The total power loss at each load level after the PV-DG installations is expressed by is the power loss of branch b with admittance y b in the distribution system at the l-th load level, V b is the rms voltage across the branch b, B is the set of all network branches, and is the set of all load levels. The objective function term of the total power loss reduction rate is given below.
The smaller the value of (2), the greater the power loss reduction with PV-DG installations.

2) VOLTAGE REGULATOR TAP CHANGES
The reduction of operation cost of voltage regulators means to have the number of tap changes of regulators as low as possible. Equation (3) shows the number of tap changes of regulators at the l-th load level.
where τ l,ϕ is the number of tap changes of the ϕ-th regulator control after connecting PV-DGs at the l-th load level, τ l−1,ϕ is the number of tap changes of the ϕ-th regulator control after connecting PV-DGs at the (l-1)-th load level. is the set of all voltage regulators. The objective function to be minimized then becomes In the study, a method of weighting sum for multi-objective optimization is used to decide the fitness value of the multi-objective function and to obtain the best solution. The weights of (4) are defined according to the degree of importance of each component of the objective function [36].

B. CONSTRAINTS
The constraints of the PV-DG planning problem include the limits of rms voltage, the regulator tap positions, and PV-DG constraints, as described below.

1) RMS VOLTAGE
The rms voltage at each bus must be maintained within an acceptable range, as given in (5) where V n,l is the rms voltage of the n-th bus at the l-th load level, V min and V max are the lower and upper limits of the system voltage profile, respectively. ℵ is the set of all system busses.

2) TAP POSITIONS OF VOLTAGE REGULATOR
In the study, the ϕ-th voltage regulator is assumed to have 2N tap taps for its regulated voltage, V tap ϕ , ranging from -N tap to N tap , as shown in the integer constraint of (6) [37], [38]. V min tap ϕ and V max tap ϕ are the minimum and maximum regulator voltages, respectively, as shown in (7).
C. PV-DG CONSTRAINTS High penetration of PVs can affect the operation of voltage regulation devices and severe voltage fluctuations. The smart inverter control of PVs, which provides different functions such as volt-var and fixed power factor, can help PVs provide reactive power support in response to dynamic variations in voltage at the point of connection [39]. For instance, Fig. 1 shows the fixed power factor function of the PV-DG. PVs can inject a constant real power (kW) and various reactive power (kVar) at a specified power factor range. As shown in Fig. 1, the PV can be operated at cosϕ 1 , cosϕ 2 , or cosϕ 3 corresponding to the output of (P 1 , Q 1 ), (P 2 , Q 2 ), and (P 3 , Q 3 ), respectively. In this study, the fixed power factor function is modeled to control the reactive power of PVs. The reactive power support function introduced by the operation of PVs is given in (8)- (10). The real power output of each PV-DG unit and the total generation of all PV-DG units for a specific load level l must be kept within the maximum generation limit, as shown in (8) and (8). The constraint of (8) is to avoid over-generation of DGs during the lowest load level, P L . Equation (10) indicates the power factor constraint used to control the reactive power of PV-DGs.
where P max PV ,j is the maximum generation limit of the j-th PV-DG unit.
is the set of available PV-DG units for placement.

III. PROPOSED METHOD AND SOLUTION PROCEDURE
The proposed two-stage PV-DG planning problem is tested under the IEEE 123-bus benchmark system and an actual utility distribution feeder. The LRSF-based placement procedure is firstly implemented to efficiently determine PV-DG candidate busses and the proposed COA-based algorithm is then applied to solve the optimal sizing problem.

A. POWER LOSS REDUCTION SENSITIVITY FACTOR-BASED PV-DG PLACEMENT
The key component in the objective function to be minimized is the real power loss reduction after sitting DGs. Assume that a given number of available PV-DG units are planned for placement in the system and each unit is assigned with a maximum kW capacity. The search space of the PV-DG candidate busses usually is enormous. For instance, it is approximately 1.731 × 10 13 combinations of candidate locations for 10 different PV-DG units to be installed in a 100-bus system. Therefore, the brute-force approach to find optimal PV-DG locations is impractical. To reduce the search space for the PV-DG placement problem, the loss reduction sensitivity factor of (11) for a specified bus i at the l-th load level is assessed by placing a small testing PV-DG unit with VOLUME 8, 2020 a capacity of (kW Cap ) PV at the specified bus one at a time.
where TPL l,i(w/oPV ) and TPL l,i(PV ) are the total real power loss of the system without and with the placement of the selected PV-DG unit at the l-th load level and at the i-th bus, respectively. After assessing LRSFs of all busses, the busses with top priorities (i.e. the highest LRSF values) are chosen as candidate busses for siting the available PV-DG units. The size at each candidate bus is then determined by the proposed COA-based algorithm, as described below.

B. OVERVIEW OF COYOTE OPTIMIZATION ALGORITHM
The COA algorithm applies a mathematical model that is described as the birth and death of ''coyote'' inside a pack.
In the COA algorithm, the population of coyotes is separated into N p packs with N c coyotes in each pack. The total population in the algorithm is obtained by the multiplication of N p and N c , each coyote is a possible solution for the optimization problem and its social conditions set, soc, include all decision variables [27]. For the c-th coyote of the p-th pack at the t-th time instant, the social condition is expressed by The first step in the COA is to initialize the global population of coyotes. For the c-th coyote of the p-th pack of the j-th dimension, soc p,t c,j = lb j + r j (ub j − lb j ), j = 1, 2, . . . , J where lb j and ub j are the lower and upper bounds; r j is a real random number uniformly generated in the range of 0 and 1. The next step is to evaluate the current social conditions by (13) representing the cost of the objective function, fit According to [14], to diversify the interactions between all coyotes in the population, the coyotes sometimes leave their packs and become solitary or join a pack instead with probability, P e .
To avoid P e is greater than 1, N c is limited below 14. Then, the new social condition of the coyote is updated by (15).
where r 1 and r 2 are uniformly distributed random numbers in the range of 0 and 1. δ 1 and δ 2 are given by (16) and (17), respectively. where Equation (16) illustrates a cultural difference from a random coyote of the pack (cr 1 ) to the α coyote and (17) shows a cultural difference from a random coyote (cr 2 ) to the cultural tendency of the pack, β. In (19), S p,t is the ranked social conditions of all coyotes of the p-th pack at the t-th time instant for each j, j = 1, 2, . . . , J . Equation (19) implies that the cultural tendency of the pack is calculated as the median social conditions of all coyotes in that pack.
The new social conditions is evaluated by (20), and the new social condition is decided by (21).
Also, the birth and the death of a coyote are considered in COA. The birth of a new coyote is written as a combination of the social conditions of two random parents plus an environmental factor, as shown in (22).
where P s = 1/J (i.e. the number of variables), P a = (1 − P s )/2, k 1 and k 2 are random coyotes from the p-th pack, j 1 and j 2 are two random dimensions of the problem, P s is the scatter probability, P a is the association probability, R j is a random number inside the decision variable bound of the j-th dimension, and rnd j is a random number in the range of 0 to 1.
The pub will survive if the fitness value with the pup smaller than the older; otherwise, the pup will die. Finally, the social condition of the coyote that best adapted itself to the environment is selected and is used as the global solution of the problem.

C. SOLUTION PROCEDURE
The following describes the details of the two-stage solution procedure including the loss reduction sensitivity-based placement and optimal sizing by using the COA algorithm for the PV-DG planning problem.

1) DISTRIBUTION SYSTEM POWER FLOW ANALYSIS
In this study, OpenDSS simulation tool is used for solving power flow problems. Because the load models have been modified in OpenDSS so that power flow solution nearly always converges for very low voltages. While the power flow solutions of other algorithms are difficult to maintain a converged solution over a wide range of voltages [31]. In OpenDSS each power deliver or conversion element is represented by a nodal admittance network model to perform the power flow solution. The power-deliver elements including lines and transformers are represented by the primary admittance matrix, Y p . A power conversion element is typically represented by its Norton equivalent with a constant Y p in parallel with an injection (or compensation) current that compensates for the nonlinear portion. The nodal admittance matrix of each element is then used to construct the system admittance matrix, Y s , where Y s is usually maintained constant for computational efficiency.
An initial guess at the voltages, V, is obtained by performing a direct solution of I = YV, where generators and loads are modeled by their linear equivalents with no injection currents. The power flow iteration starts with obtaining the injection currents from all the power conversion elements in the system and updating them in the injection current vector, I inj . The solution is focused on solving the nonlinear system admittance equation of the form of I inj (V) = Y s V, where I inj (V) is a function of voltage and represents the nonlinear part of the currents from loads, generators, and PV-DGs in the circuit. To solve the nonlinear equations set, a fixed point method shown in (23) is adopted.
The iteration continues until the convergence criterion for the voltage vector is satisfied. This simple iterative solution has been shown to converge well for most distribution systems that have adequate capacity to serve the load demand. When performing yearly simulations such as in our study, the solution at the present time step is used as the starting point for the solution at the next time step. The solution typically converges in two iterations. Therefore, the OpenDSS efficiently performs the power flow calculations [31].

2) PROCEDURE FOR PLACEMENT OF PV-DGs
The placement procedure considers all load levels based on the installation of a testing DG unit by injecting real power at each bus one at a time. Then, the loss reduction sensitivity factor, LRSF, for each bus is calculated. The busses with the top rank of LRSF values obtained by (11)  If yes, proceed to the next step. Otherwise, return to step 2 for the next load level.
7. Select the top-priority K (K <M ) busses from the M busses of the highest load level which are also in the top M busses of each of the other load levels obtained at step 5.
Since the top K buses at the highest load level have higher LRSFs compared to lower load levels, they are selected as candidate placement busses for all load levels.

3) PROCEDURE FOR SIZING THE PV-DGS AT CANDIDATE BUSSES
The following procedure summarizes major steps for sizing the DGs at the selected candidate busses by using the COA algorithm. The sizes and power factors of PV-DGs at candidate busses and tap positions of voltage regulators are defined as the social conditions of the coyote (i.e. a solution).
1. Parameters initialization. Assign the number of packs N p , the number of coyotes in a pack N c , the number of social conditions for each coyote (i.e. a possible solution), and the number of load levels L. Specify the lower and upper bounds of each social condition (i.e. capacity limits of the PV-DG unit at a candidate bus). 2. The initial social conditions are randomly set for each coyote for all load levels. Check bounds for each social condition in the coyote. If bounds are violated, initialize social conditions again until no bound violation. Calculate the fitness value of (4) using the social conditions of the initial coyotes. 3. Specify the maximum number of iteration N i and perform the COA procedure as follows. for each iteration for each load level for each pack -Define the α coyote of the pack.
-Compute the social tendency of the pack using (18). for each coyote of a pack -Update the social condition of coyote using (15).
-Check bounds condition. If the bounds of the social condition are violated, update again. -Calculate the fitness value using the new social conditions of the initial coyotes by (20). -Adapt the social condition using (21). end Perform birth and death inside the pack using (22).
-Check social-condition bounds. If bounds are violated, give birth to a new pup again. -Calculate the fitness value with the pup. If the fitness value with the pup is smaller than the older, the pup will survive. Otherwise, the pup will die. end end Select the coyote with the best fitness value. Check if a coyote can leave the pack and enter another pack according to (14). Then, update the pack information. end VOLUME 8, 2020 4. Select the best adapted coyote with its social conditions (i.e. the optimal sizes, power factors of PV-DGs and tap positions of voltage regulators) from Step 3. 5. Output the fitness value with the best adapted coyote's social conditions (including the optimal size and the power factor of PV-DG at each candidate bus, tap positions of voltage regulators, total power loss, and voltage profile). It is noted that, in the above solution procedure, the initial randomly assigned PV-DG sizes, power factors, and regulator tap positions in Step 2 will be updated in Step 3 for each load level and each iteration. Fig. 2 depicts the flowchart of the proposed procedure for sizing the PV-DGs at candidate busses.

IV. TEST RESULTS
In the study, there are two fairly sizable feeders are under test to show the usefulness of the proposed method. One is the ieee 123-bus benchmark distribution network and another one is an actual Taipower 137-bus distribution system. Assume that 10 PV-DGs are available and each DG capacity is in the range of 50 to 1000 kW and is an integral multiple of 50 kW. the EPRI OpenDSS is used to perform power flow analysis for assessment of the loss reduction sensitivity factor, LRSF,  of (11) at each bus. The priority list of the LRSF values is then determined for the top 10 candidate busses for PV-DG placements. Then each solution (i.e. coyote) is input to OpenDSS through the common object model interface to perform power flow at Steps 2 and 3 of the sizing procedure for social condition assessments and solve the planning problem. Fig. 3 illustrates the diagram of co-simulation between OpenDSS and Matlab, where the COA-based approach is implemented using Matlab. Results include rms voltage at each bus, regulator tap positions, and the system power loss. The fitness value of (4) is then calculated. Results obtained by COA are compared with those obtained by mixed-integer nonlinear programming (MINLP), genetic algorithm (GA), particle swarm optimization (PSO), biogeography-based optimization (BBO), and grey wolf optimizer (GWO).
In the two test cases, the number of annual (i.e. 8760 hrs) load levels, L = 3, and the number of continuous and integer solution variables for case 1 is 450 and for case 2 is 471. the parameter settings of all compared algorithms are listed in Table 1. The parameters are chosen for each algorithm based on those proposed methods in the indicated references.

A. CASE 1: IEEE 123-BUS BENCHMARK DISTRIBUTION SYSTEM
In the study, three load variation levels are included. the highest annual load level is with a peak value of 3490 kW. The medium and low annual load levels are 80% and 60% of the highest load level, respectively. To evaluate the multi-objective function of (4), the weighting factors are γ 1 = 0.7 and γ 2 = 0.3, depending on the importance of each objective component. After applying the COA algorithm for sizing the PV-DGs, the voltage at each bus is limited within the range of 0.95 to 1.05 pu. Table 2 shows the priority list of candidate DG busses for all load levels calculated by (11) after the placement procedure. Table 3 lists the top 10 busses which have the highest LRSFs for all load level selected as the candidate PV-DG busses. The power factor of each PV-DG unit is to be maintained within the range of 0.85 to 1. Each of the seven single-phase voltage regulators has ±16 taps and the corresponding voltage ranges from 0.9 to 1.1 pu. The results shown in Table 4 are the best solutions among 25 independent runs of all methods at all load levels, which include each of the top 10 candidate PV-DG busses with its optimal size in kW and the range of PV power factor, as well as the total installed kW capacity obtained by each compared method. It also indicates that the proposed method leads to the highest penetration capacity among the compared methods. Table 5 lists the loss reductions, bus rms voltages, computational time, and fitness values of (4) obtained by all compared   methods at all load levels before and after PV-DG planning. Table 6 shows the regulator tap positions and the number of tap changes obtained in Step 3 of the sizing procedure.
In table 5, it is seen that the fitness value obtained by the COA method is the lowest among the compared methods and leads to the best solution with a loss reduction ratio of 63.898%. The power loss for all load levels before PV-DG planning is 191.913 kW. Fig. 4 depicts the convergence trend of each method. It is observed that the COA method yields superior convergence than the other methods.    The rms voltage profile is improved and is within the range of 0.99 to 1.05 pu after PV-DG planning obtained by the proposed COA method, as shown in Fig. 5, where the minimum and maximum rms voltages become 0.9912 pu and 1.0489 pu, respectively. Fig. 6 illustrates that tap positions before and after the PV-DG planning obtained by the proposed method. It is observed that the tap positions of all voltage regulators remain unchanged.

B. CASE 2: TAIPOWER 137-BUS DISTRIBUTION SYSTEM
In this case, a Taipower distribution system (69 kV/11.4 kV) is used to test the usefulness of the proposed method. The number of PV-DG units and their sizes are the same as those given in Case 1. Results for comparisons include rms voltages at each bus and the system power loss. The multi-objective   functions of (4) become a single objective function. Since there are no voltage regulators in the system, γ 2 = 0.
There are four feeders in the study system: XD21, XH22, XH21, and XO32. Fig. 7 depicts the single-line diagram of this taipower distribution network. In the system, the peak value of the highest annual load level is 6000 kW. The medium and low annual load levels are 80% and 60% of the highest load level, respectively. Table 7 shows the best solutions among 25 independent runs of all methods at all load levels, which include the top ()10 highest-LRSF candidate busses with optimal kW sizes and the ranges of power factors, as well as the total installed PV-DG capacities obtained by the methods under comparison. Table 8 lists the loss reductions, bus rms voltages, elapsed time, and fitness   values of (4) obtained by all compared methods at all load levels. Results in Table 8 also shows that the fitness value obtained by the COA method is still the lowest among the compared methods and leads to the best solution with a loss reduction ratio of 77.01%.
Simulations based on each compared method are performed at 500 iterations. Fig. 8 represents the convergence trend of each method. In this case, the COA method also performs better than other methods. Fig. 9 depicts the three-phase rms voltage profile corresponding to the number of buses for each phase after installed PV-DGs at all load levels. It is observed that the minimum and maximum rms voltages are 0.9939 pu and 1.0023 pu, respectively, obtained by COA. Table 9 summarizes the maximum PV-DG bus voltages for both cases and it is observed that the bus voltages are well controlled within the allowed upper limit.

V. CONCLUSION
This paper has proposed a two-stage solution algorithm by applying the COA method for PV-DG planning in distribution feeders considering system loss reduction, rms voltage profile improvement, and regulator tap controls under different load levels. The proposed method has been tested on two practical distribution systems. Test results show that the rms voltage at each bus can be improved, and the total number of tap changes of voltage regulators is well controlled. By comparing with one conventional and four other bio-inspired metaheuristic optimization approaches, the results confirm that the proposed method is superior to the compared methods in both power loss reduction and hosting capacity for PV-DG planning.
Though the proposed COA-based approach is not the most computationally efficient compared to the other methods. The solutions obtained by testing both cases show that the maximum PV-DG hosting capacity, the minimum power loss, and the minimum number of tap changes can be achieved by the proposed solution algorithm. For Case 1, the proposed method can install an additional 50 to 300 kW of PV-DG capacity and reduce an additional 1 to 8% of system power loss with no regulator tap changes. For Case 2, the proposed method installs an additional 150 to 700 kW of PV-DG capacity and an additional 1.5 to 11% of system power loss reduction. The study has demonstrated that COA is better than the compared methods in solving the PV-DG planning problem for practical distribution feeders. The proposed method also can be further applied to solve other power system planning problems such as economic dispatch, optimal power flow, and unit commitment.
At the present phase of the study, the PV-DG planning does not model uncertainties of PV generation intermittency associated with seasonal output fluctuations and annual load variations. It is important to appreciate that there is no conflict between the planning either with or without including uncertainties. Both planning processes can be combined and lay a solid foundation for efficient PV-DG planning. Considering PV-DG generation and load uncertainties can be an enhancement of deterministic planning and provide results closer to reality if the uncertainty models are satisfactory accurate. Future work will take the uncertainties into account and include energy storage in the planning problem.