Improved Coyote Optimization Algorithm for Optimally Installing Solar Photovoltaic Distribution Generation Units in Radial Distribution Power Systems

Power System Optimization Research Group, Faculty of Electrical and Electronics Engineering, Ton Duc !ang University, Ho Chi Minh City, Vietnam Faculty of Automation Technology, !u Duc College of Technology, Ho Chi Minh City, Vietnam Faculty of Electrical and Electronics Engineering, Ho Chi Minh City University of Technology and Education, 01 Vo Van Ngan, !u Duc District, Ho Chi Minh City, Vietnam Institute of Research and Development, Duy Tan University, Da Nang, Vietnam


Importance of Solar Photovoltaic for Distribution Power
Networks. Nowadays, the production of electricity from fossil fuels has been causing many negative impacts such as air pollution, water pollution, and increasing temperature. So, reducing the use of fossil fuels is extremely essential. In fact, in many countries around the world, solar photovoltaic distribution generation units (PVDGUs) are being popularly used as a main renewable energy source due to its economic, technical, and environmental benefits [1]. For example, in California, solar power systems with the capacity of more than 500 MW have been connected to highvoltage power lines, reaching high benefits. In general, during the last decades, the annual growth rate of distributed generation based on solar photovoltaic increased higher than 40% [2]. Many studies have shown that the received benefits will depend on the site and rated power of distributed generation units (DGUs) in general and PVDGUs in particular [2][3][4][5][6]. Proper installation of DGUs can reduce total power losses and improve voltage profile and power quality [7,8]. However, with the inappropriate site and rated power, DGUs will likely cause increase of many undesirable problems in power systems such as voltage flicker, voltage sags, fault current, harmonic distortion, and power losses [9]. erefore, the task of determining the most suitable site and the most optimal rated power of DGUs in the distribution systems plays a very important role [7].

Related Work.
Most of the approaches that have been used for searching the optimal location and sizing of DGUs so far consider power loss reduction and voltage profile enhancement as the main objectives. Genetic algorithm (GA) is one of the most common evolutionary methods and has been applied for the considered problem [1,10,11]. ese authors applied GA to find the suitable location and sizing of DGUs for improving the voltage profile as well as power losses reduction, but total harmonic distortion (THD) and individual harmonic distortion (IHD) have not been considered. e voltage stability and loss reduction were really effective after properly installed DGUs in the distribution system. GA is very common in solving the optimal problems, but the biggest disadvantage is not to guarantee for finding global solutions due to its premature convergence to local optimization zone. In addition to GA, authors in [12] also proposed loss sensitivity factor (LSF) for sizing optimization and applied simulated annealing (SA) to determine the optimal location of DGUs. By combining these two proposed methods, the optimal solution of DGUs was successfully found with the goal of reducing branch losses. is research was implemented at different system power factors, and the obtained results were compared with other methods. Simulation results have shown that this method was better than other methods at different power factors.
is is a good method, but the biggest disadvantage of this method is the complicated calculation process. us, it has not been widely applied for other problems regarding distribution networks. With the same objective of minimizing power losses on branches, the authors in [13,14] have also contributed new algorithms for finding optimal solutions in the power distribution network. To address the suitable position and capacity of DGUs in the distribution system, Plant Growth Simulation Algorithm (PGSA) and Modified Plant Growth Simulation Algorithm (MPGSA) were proposed in [13,14], respectively. Based on the obtained results from the simulation, the proposed solution from the MPGSA method was better than that of the PGSA as compared in the same objective function, the same constraints, and the same system data. is showed that the improvement from the PGSA method (MPGSA) had significant effectiveness. Both PGSA and MPGSA methods depend too much on setting the value of the various parameters in the implementation process of evolutionary. erefore, the quality of the optimal solution largely depends on this initial parameter setting. In the same field of DGU optimization research, different author groups have different solution methods. While the authors in [15] chose selective particle swarm optimization (SPSO), the authors in [16] chose Bat Optimization Algorithm (BAT). Authors that used the SPSO method have done their research on two distribution systems. e result was very satisfactory by installing 3 DGUs in the IEEE 33bus distribution network and 2 DGUs in the IEEE 69-bus distribution network. In addition, BAT also reached very effective results through the installation of various types of DGUs in the IEEE 33-bus distribution network. e results from BAT were superior to other methods with the higher quality solutions. BAT has a fairly fast convergence rate in the first stage, but in the next stage, the convergence rate slows down. In this algorithm, the mathematical analysis to link the parameters with convergence rates is limited. erefore, the accuracy is not high. In addition to BAT, particle swarm optimization (PSO) has been used to deal with the study that DGUs were integrated into the power grid [17]. PSO is one of the popular methods, but its application for the problem is not complicated due to the very simple characteristic. Like GA, one of the most disadvantages of PSO is that it easily falls into local search areas if its search area is too wide. erefore, its convergence to a global optimum is hardly ever guaranteed. ere was not enough evidence to conclude the performance of PSO because only the power loss objective of the IEEE 33bus distribution network has been considered for comparison. Multiple-objective particle swarm optimization (MOPSO) in [18] could reach power loss reduction and increase purchasing power. In [19], human opinion dynamics (HOD) has successfully dealt with the IEEE 14-bus and 30-bus distribution systems.
e algorithm used to build in this article is complex but effective. e main component of the algorithm consists of four pillars as social structure, opinion space, social influence, and update rule. e obtained solution from the proposed method (HOD) shows superiority over other methods. Moreover, author [20] also proposed a method called salp swarm algorithm (SSA) for finding the location and capacity of renewable distributed generators in the 33-bus radial distribution system. Although SSA could reduce losses and operation cost, it was not a good method for the complex objective function. In fact, as observing from convergence characteristic, it suffered premature convergence and could not jump out local zone to more promising search space. Big Bang-Big Crunch method (BB-BC) [21] and Artificial Bee Colony (ABC) Algorithm [22] had the advantage of fairly fast convergence and few parameters, but the biggest disadvantage of them was the early convergence in the searching stage and the accuracy of optimal value cannot be satisfied. In [23], Biogeography-based optimization (BBO) was applied for multiobjective problem considering power losses, operation cost, and the lowest voltage. Different values have been set to weight factors of the three objectives, and the three objectives from BBO could not be better than those from other compared methods. In general, studies from [12][13][14][15][16][17][18][19][20][21][22][23] have tried to reduce power loss and operation cost and improve voltage, but they have ignored one important 2 Complexity factor, harmonic distortions. Two components to evaluate harmonics, THD and IHD, have not been considered as constraints or objective functions in the studies. Actually, as installing DGUs in the distribution system, total harmonic distortion (THD) and individual harmonic distortion (IHD) would be changed depending on the location and capacity of DGUs [24,25]. Study [24] could be a valuable contribution to the distribution system since some types of distributed generation unit (DGU) have been installed in the radial distribution system with many nonlinear loads for testing improvement levels. Four types of DGU that have been proposed to be installed in distribution networks consisted of (1) injection of pure active power, (2) injection of pure reactive power, (3) injection of both active and reactive power, and (4) injection of active power but consuming reactive power. e best type of DGU was similar to the usual generator, i.e., the third DGU type supplying both active and reactive power to the distribution network. In addition, the study has proved that power losses could reduce to 51.16%. However, there was no enough evidence for concluding that the type and applied GA were the best options because only one distribution system with 10 buses was employed as study case. Likewise, the study [25] applied sensitivity analysis and Gravitational Search Algorithm (GSA) for optimizing average harmonic distortion and power loss of the IEEE 69-bus distribution network. In addition to the demonstration of the contribution of DGUs to harmonic distortion reduction, GSA has also been compared to GA and evolution program (EP). GSA has shown better THD and power loss, but the reduction was not significant. In [26], four methods consisting of GA, PSO, ABC, and BBO have been implemented for locating PVDGUs and searching its optimal capacity under constraints of distribution networks and harmonic standards. A biobjective function (power loss minimization and harmonic distortion minimization) has been concerned in which harmonic distortion was the average of THD and IHD. e study has considered power loss minimization to be a priority while THD and IHD have been constrained not higher than 5% and 3%, respectively. e authors have demonstrated that using these applied methods for determining PVDGUs could reduce power losses significantly as compared to the case without PVDGUs. In addition, the authors have recommended BBO should be used rather than three other ones.
Generally, all concerned studies about determining the optimal location and capacity of DG units in radial distribution networks have three main groups with three shortcomings as follows: (1) Neglecting the negative impacts of harmonic currents. THD and IHD have not been considered to be either objective function or constraints.
(2) Considering single objective and reducing the negative impacts of harmonics satisfying IEEE Std. 519 but not demonstrating the performance improvement of the applied method clearly. Improvement level of the applied method over compared methods has not been shown and comparison of convergence speed has not been considered. (3) Considering THD, IHD, and other objectives such as  power loss, voltage profile, and operation cost in  multiobjective function. ere was a trade-off between the obtained objectives. It means some objectives were better, but others were worse. So, it cannot lead to a clear evaluation of the real performance of applied methods and compared methods.

Proposals, Novelties, and Contributions.
Realistically, studies in [24][25][26] have considered THD and IHD and the two factors have been reduced not higher than 5% for THD and 3% for IHD. e two factors always satisfied IEEE Std. 519 [26]. However, the studies in [24][25][26] have not shown the clear performance of applied methods since the results obtained by their applied methods could not be improved better than other methods significantly. Furthermore, the studies in [25,26] were about multiobjective problem with trade-off between two objectives. In this regard, we have found four studies [27][28][29][30] that proposed single objectives for the task of determining the location and sizing of PVDGUs. As pointed in [27], the power losses occur naturally and depend on the type of conductor used, the capacity of the transformer, and other components used in the distribution grid.
us, these power losses cannot be removed but should be minimized for reducing energy losses. Study [28] has focused on the economic goal associated with the initial investment capital for buying PVDGUs. us, reducing the penetration level of PVDGUs is taken into consideration as an objective function. On the contrary, study [29] has considered technical issue as a main goal and improving voltage profile has been proposed to be an objective function. Study [30] has pointed out serious damages of nonlinear loads on distribution grids since they have caused voltage distortion. Hence, the main goal should be the minimization of harmonic distortions. Although the four studies have proposed different single-objective functions, they have had the same viewpoint about technical issues such as pure voltage wave, standard voltage magnitude, and stable operation of girds. For assuring the technical factors such as voltage limits of loads, current limits of conductors and harmonic distortion limits have been constrained seriously. e four single objectives together with constraints of technical factors are considered in the study by the four following cases: In the above four cases, the different harmonic sources are injected into linear loads in two considered systems and make them become nonlinear loads. en, forward/backward sweep technique (FW/BWST) [31] is used to solve the power flow and the harmonic flow. Solving the power flow can obtain the voltage of load buses and branch currents at the fundamental frequency (the first-order harmonic) while solving harmonic flow can reach the two parameters at a frequency of considered harmonics. e voltage and current are directly related to total power losses, voltage profile, and harmonic distortion. e objectives can be highly effective if the location and site of PVDGUs are appropriate.
In recent years, COA has been widely applied for optimization problems in many different fields. e authors in [32] have applied COA to minimize the gas consumption of turbines in combined cycle power plants in Brazil. e proposed solution fully satisfied the physical limits of the turbine and pollution emission regulations. COA method showed superiority to other methods such as Artificial Bee Colony (ABC), Binary Switching Algorithm (BSA), Selfadaptive Differential Evolution (SaDE), Genetic Whale Optimization Algorithm (GWOA), Symbiotic Organism Search (SOS), and PSO. Authors in [33] were successful in applying COA to find the main parameters of three diodes in photovoltaic modules. But the obtained results have not been compared to any results from other methods. In [34], the economic dispatch problem with thermal power plants and wind turbines was solved by COA, GA, and PSO. e obtained results from two studied systems showed that the COA reached a better solution than GA and PSO, but there was no demonstration about the faster speed of COA since settings of population and iterations were ignored. In [35], COA was applied for more complicated problem than previous applications, which was the problem of optimizing the reactive power flow for transmission power systems. However, the solution proposed by COA was not as good as other methods. ereby, the effectiveness of COA was not really convincing in many different problems. is method needs to be reviewed and modified to improve the operational efficiency. In fact, in order to find a solution for any optimal problem, it is necessary to apply the optimal method with high efficiency. In other words, the effectiveness of solutions primarily depends on the applied method. Another demonstration of the inefficiency of COA needs to be improved in [36]. In that study, COA has been applied for dealing with block matching problem and the results of COA were compared with Enhanced Grey Wolf Algorithm (EGWA) and other methods. Discussion of results in [36] indicated that COA was less effective than EGWA and some of the other methods such as Black Hole Algorithm (BHA), Grey Wolf Algorithm (GWA), and PSO. e review on previous applications of COA indicated that COA was not really highly effective for different optimization problems in the engineering field. Clearly, COA is currently owning main disadvantages leading to the restriction of convergence to global optimum. To overcome the disadvantages of COA, we propose an improved coyote optimization algorithm (ICOA) for optimizing the placement of PVDGUs and determination of rated power of these PVDGUs for different objectives of two power systems, IEEE 33-bus distribution network and IEEE 69-bus distribution network. ICOA is an improved metaheuristic algorithm by performing two main modifications on the original coyote optimization algorithm (COA). e two proposed modifications aim at improving optimal solution quality found by the first generation and the second generation in the COA method. COA method was developed in 2018 [37] based on the natural behaviors of coyotes. Each coyote is characterized by two main factors, social condition and quality of the social condition in which social condition is corresponding to optimal solution and quality of the social condition is corresponding to fitness function of the solution. e coyote community is divided into N g small coyote groups, and there are N co coyotes in each group. COA method produces two new solution generations per iteration in which the first generation produces N co new solutions for each group and the second generation produces N g new solutions for the whole coyote community. us, total number of new solutions generated in the COA method is (N co × N g + N g ) solutions in which N co × N g is equal to population (N pop ). us, (N co × N g + N g ) is equal to (N pop + N g ). e number of new solutions for each generation can indicate that the first generation has a more significant impact on the final solution quality because it produces N co times new solutions of the second generation. However, COA is coping with low performance of the first generation because it uses a center solution for producing updated step size; meanwhile, the strong point of the center solution is only to produce diversity for solutions, but it is not potential for generating promising updated step size. In the second generation, COA produces one new solution for each group by using a randomization factor. e new solutions are formed by randomization and control variables in the solution can be either taken from the current solutions or randomly produced within lower and upper bounds. However, it must use one out of three different selections for each new control variable in the solution and two control parameters need to be determined for the selections. Clearly, the second generation is much dependent on randomization and it is time-consuming due to the determination of two control parameters. Consequently, in ICOA method, we propose two modifications on the first and the second new solution generations. In the first generation, center solution is replaced with the so-far best solution with the intent to improve the quality of newly produced solutions and reduce simulation time. In the second modification, each group produces one new solution around the best solution by using either two or three updated steps dependent on the number of pair solutions, which are closed together. e second modification can enhance the effectiveness of local search and find one promising solution for each group. e proposed ICOA method together with six other published methods consisting of COA, BBO, PSO, SFO, SSA, and GA is employed to find the optimal location as well as the optimal capacity of PVDGUs connected in the IEEE 33-bus and IEEE 69-bus distribution networks. e obtained results from ICOA in four cases are compared with those of other methods for performance evaluation. 4 Complexity In summary, the main contributions of this study are as follows: (1) e contributions to the issue of installing PVDGUs in radial distribution networks: (i) Radial distribution networks with the installation of PVDGUs can reduce the negative impacts of harmonic flow, improve the voltage profile of load buses, and reduce total power losses. (ii) Proposing four single-objective functions and considering all constraints for each one. (iii) Indicating particular benefits for each single objective. e benefits help the readers to easily identify the most appropriate single objective of installing PVDGUs for their purpose.
(2) e contributions to the performance improvement of the proposed method: (i) Pointing out the advantage of each one between two proposed modifications on the proposed ICOA (ii) ICOA can have some advantages over COA such as (i) reducing computational time, (ii) improving quality of solutions, and (iii) reaching better stability in the search process. (iii) e results show that ICOA is more suitable than BBO, PSO, GA, SFO, and SSA for installing optimal PVDGUs in distribution networks. e reader should try other methods for the problem.
e remaining parts of the paper can be divided into 6 sections as follows: description, assumption, objective functions, and constraints for four cases are shown in Section 2. e whole search procedure of the coyote optimization algorithm (COA) and improved coyote optimization algorithm (ICOA) is described in detail in Section 4. e whole computation process of using ICOA for the considered problem is shown in Section 5. Simulation results from the systems of 33 buses and 69 buses are compared for evaluation in Section 6. Advantages and disadvantages as well as selected control parameters for setting in the proposed method are also discussed in Section 7. Finally, the conclusions of the work are stated in Section 8.

Problem Formulation
To achieve the goal and fully satisfy all constraints, PVDGUs must be connected in the power system with optimized place and size. However, the installation of PVDGUs always requires power converters and this can cause a negative impact on the distribution system similar to nonlinear loads, i.e., the presence of harmonics. In order to prevent the presence of harmonics from power converters, many manufacturing technologies of power converters have been innovated with the consideration of harmonic filters seriously [38,39]. erefore, in the study, all PVDGUs connected in the systems have been checked for pure sinusoidal voltage waveforms by applying harmonic filters.

Assumption.
e optimal problem for a solar photovoltaic integrated radial distribution system consists of seven components such as (1) main power supply (fossil energy source); (2) solar photovoltaic distribution generation units (PVDGUs); (3) distribution lines; (4) loads (liner loads and nonlinear loads); (5) buses; (6) power conversion technology; and (7) monitoring and control system. e information for all of these components is detailed as follows: (1) Main power supply: work as a slack bus and provide sufficient power capacity for the system to operate stably. (2) Solar photovoltaic distribution generation units: provide the active power for the distribution system through the power converter integrated the energy storage system. Active power and reactive power remain unchanged during the considered period. (6) Buses: voltage magnitude and phase angle are accurately measured. (7) Power converter system: this includes DC-DC converters and DC-AC inverters which are supposed to be working with the efficiency of 100% (without any losses). (8) Harmonic of PVDGUs: by using high-technology power converters with the harmonic filters inside [38,39], it is supposed that there is no harmonic produced from converters and PVDGUs do not produce any harmonic in radial distribution networks. (9) Monitoring and control system: in addition to converters, the monitoring and control system is also essential for the distribution system. is management system can control the operation of energy store system (ESS) and the charging and discharging process [40].
By applying suitable algorithms, the optimal problem for this integrated distribution system is effectively addressed, satisfying technical criteria and economic aspect.

Objective Function.
In considered systems, we focus on four objective functions including total active power losses (TPL), capacity of PVDGUs (RP), voltage profile (VP), and Complexity 5 harmonic distortion (THD & IHD). e objective function along with the constraints is explained by using mathematical formulas as follows.

Total Active Power Losses.
Reducing total active power losses (TPL) due to PVDGU installation can reduce the cost of transmission economically and contribute to improving the stability of the power system. erefore, the reduction of total active power losses is becoming an important objective, which is shown in the following mathematical equation [41]: By using the same method, total power losses in all branches before installing PVDGUs, TPL, is determined by means of the following equation [41]: It shows that if PVDGUs are properly installed, the losses of power system after installation are always smaller than before installation. erefore, for convenience of consideration, another objective function is modified to keep the result in the restricted range between zero and one as shown in the following equation:

e PVDGUs' Penetration Level.
Integration of large distributed generation into electric system can increase investment and operating costs. erefore, minimizing the capacity of PVDGUs (RP DG ) that still satisfies the technical criteria should be considered as the target. e proposed mathematical equation of the objective function is shown as follows [42]: Since the total capacity of PVDGUs is always less than or equal to the total power of the load demand, so the value of the equation (4) always varies from zero to one.

Voltage Profile Index.
One of the important criteria for evaluating the stability of the operating grid is the voltage profile index (VPI). us, the voltage profile index must be considered and the calculation process for VPI should follow the step below [43]. e voltage profile for the i th bus can be defined as where V nom is the nominal or desired bus voltage value and typically taken as 1.0 pu. It is clear that VP i can get a maximum value (1.0 pu) if the bus voltage equals to V nom and VP i can get a negative value if the bus voltage is higher than the maximum limit or lower than the minimum limit. e overall voltage profile index for the distribution system can be defined as e maximum value of VPI is equal to 1.0 pu. So, to keep the value of the objective function varying in the range from zero to one, the mathematical equation is modified and shown as follows: 2.6. Harmonic Distortions. In this study, linear loads are injected by the harmonic current. is creates nonlinear loads in the two considered systems and it makes the current drawn unlike the supply voltage waveform. e flow of harmonic currents creates voltage harmonics when running through the line impedance and distorts the supply voltage of the system. Harmonic distortions which include THD and IHD must be considered in the objective function.
Total voltage harmonic distortion at each bus can be defined by [44] Individual harmonic distortion at each bus can be defined by [44] According to IEEE Std. 519 [44], THD V,i and IHD h V,i should not be greater than 5% and 3%, respectively. e average of THD V,i can be defined by [25] THD Vaver � Similarly, the average of IHD h V,i can be defined by Here, the objective function should be shown as follows: where ω 1 and ω 2 are, respectively, weighted factors with THD Vaver and IHD h Vaver .ω 1 and ω 2 are constrained by 6 Complexity

e Power Balance Constraints.
In this case, the power losses at the higher order frequency are too small and the power losses are mostly at the fundamental frequency in the distribution system [26]. us, power balance constraint is considered at the fundamental frequency and has the following form: 3.2. e Voltage Limits. e rms bus voltage magnitude and voltage limits are calculated and should be maintained in the voltage limits [44,45] as follows: According to IEC Std. 50160, for low-voltage and medium-voltage electricity distribution systems, the lower bound and the upper bound of bus voltage are 0.9 and 1.1 in pu, respectively. However, the range from 0.95 to 1.05 in pu can be considered as the best range for voltage stability [46]. is range is selected to become the constraint in the 33-bus distribution system and 69-bus distribution system.

THD and IHD Limits.
Both 33-bus distribution network and 69-bus distribution network systems have many nonlinear loads, so determining the limits to accept harmonic distortions in that system is essential. e limits of total voltage harmonic distortion and individual voltage harmonic distortion should follow IEEE Std. 519 [44]: where THD max and IHD max are 5% and 3%, respectively.

PVDGUs' Capacity Limit.
e penetration rate of PVDGUs should be surveyed to make preselection because it depends on investment capital, geography, environment, etc. In this study, PVDGUs are considered as a type of distributed generator which only generates active power.
us, the active power of each PVDGU will be set within predefined limits and the sum of the active power of all PVDGUs must not exceed the load demand [26]. e PVDGUs' capacity is constrained by the following inequalities: where RP min DG and RP max DG are the minimum and maximum rated powers of PVDGUs, respectively.

e Branch Current Limits.
e branch current must be kept in the allowable limit of conductor [47] as indicated in the model below: where I nDG and I max n are the current magnitude of n th branch with PVDGUs and the maximum current magnitude of each branch, respectively.

Basic Description of Coyote Optimization Algorithm.
In order to understand coyote optimization algorithm (COA), the following descriptions should be concerned before approaching the main steps of COA method: (1) Each coyote in the community is characterized by social condition and quality of social condition in which social condition is corresponding to a candidate solution and quality of the social condition is corresponding to the quality of solution reflected by fitness function value. (2) Coyote community is divided into N g groups where each group has N Co coyotes. us, population size of 3) e k th coyote in the g th group has two important factors consisting of social condition Co k,g and quality FF k,g in which the former is a solution and the latter is the fitness function of the solution.
e COA method consists of the following steps: Step 1. Create initial coyote community e initial coyote community is corresponding to a set of initial solutions in which each solution is a social condition of each coyote. However, the initialization step is performed for each coyote group as follows: where Co min and Co max are, respectively, lower bound solution and upper bound solution obtained by where var min x and var max x are the minimum value and maximum value of the x th control variable, respectively.
Step 2. Determine fitness function value Social conditions of coyote community are evaluated. is action is similar to the task of determining the fitness function of solutions. e effectiveness level of each new solution is ranked based on a fitness function, which is Complexity normally obtained by the sum of the objective function and penalty function. A typical fitness function for a general optimization problem is expressed as follows: where OF k,g is the objective function of the k th solution in the g th group; PF is a penalty factor, which is selected by experience; PT k,g,m is the penalty term for violating the m th constraint of the k th solution in the g th group; and M is the number of constraints that are taken into account in the considered problem.
In each group, the solution with the best fitness function is called the local best solution and symbolized by Co best,g .
Step 3. Update solutions for each group Updating new solutions is corresponding to updating social condition for coyotes in each group. In COA, new solutions are found around old solutions with two distances based on the local best solution, two randomly picked solutions and a middle solution.
e detail of the newly updated solution strategy can be expressed by the following model: In equation (25), Co 1,g and Co 2,g are two randomly picked solutions in the g th considered group while Co mid,g is the middle solution obtained from the g th group. Each decision variable in the middle solution is obtained by one out of two ways shown in the following equation: where var x,mid1 and var x,mid2 are the x th center decision variables corresponding to two different cases of N Co . If N Co is an odd number, var x,mid1 is taken from the x th decision variable at the position [(N Co + 1)/2]. On the contrary, var x,mid2 is the average of two x th decision variables at two positions, N Co /2 and [(N Co + 1)/2]. However, before determining the two variables in equation (23), each decision variable of all available solutions is sorted in descending order.
Step 4. Evaluate the effectiveness level of new solutions e new solutions are evaluated by calculating fitness function (21).
Step 5. Perform selection technique Now, each coyote k has two social conditions, old social condition corresponding to old solution, Co k,g , and new social condition corresponding to new solution, Co new k,g . e quality of the two social conditions is the fitness function of the two solutions, FF k,g and FF new k,g . So, it should retain one social condition for each k th coyote by using the following rules: Step 6. Create one new social condition in each group In the step, each group g th produces one new solution Co new g by using a randomization mechanism. e new solution and the technique are shown in the following models: var new where var new x,1,g and var new x,2,g are two x th variables randomly picked from two available solutions; var new x,g is the randomly produced variable within lower and upper boundaries; and β is a random number in the range of [0, 1]. en, the new solution is evaluated by using (21) and its quality is symbolized by FF new g .
Step 7. Determine the worst solution and replace it Each new solution in each group obtained by Step 6 will be compared to the worst solution in each group Co worst g for deciding if the worst solution should be replaced with the new one. e decision is taken over by the comparison below: Step 8. Exchange solutions among N g groups To diversify solutions and avoid falling into local optimal zones with premature convergence, COA carries out solution exchange action among N g groups. Two randomly selected groups will provide two randomly picked solutions and their position is then exchanged. e task is dependent on the randomization mechanism below: where ϕ is a random number within the range from 0 to 1. And solution exchange action will take place if the condition in equation (28) is true. Clearly, the number of coyotes in each group is very important for exchanging solutions. e higher N Co is, the higher the probability of exchanged solutions is.
Step 9. Determine the best solution Co Gbest Before terminating one computation iteration, the best solution among N g groups is determined by comparing 8 Complexity fitness function. If the current iteration is the last one, the best solution Co Gbest is one candidate solution of the current run.

Improved Coyote Optimization Algorithm.
As observing equations (22) and (26), COA has produced two generations in each iteration. All solutions in all groups are newly updated in the first generation whereas only one solution is newly produced for each group in the second generation. So, the total number of new solutions is generated in each iteration is (N Co × N g + 1 × N g ) equaling to (N pop + N g ). us, the new solution generation strategy for the first generation has higher contribution to the final solution quality; however, the original technique of COA for the first generation is coping with many disadvantages. Although the second generation produces only one solution for each group, the sole solution has an extremely important role in preventing COA from falling into local optimal zones fast. e sole solution is produced based on the randomization technique, but it must need many parameters and computation steps for randomly producing or randomly choosing each variable in the sole solution. Consequently, we propose two improvements that should replace the two original generations of COA.

e First Improvement in the First Generation.
In the first improvement, the center solution Co mid,g in equation (22) is selected to be replaced with the best solution in the current population. In COA, the impact of the center solution is significant in finding optimal solution as solving benchmark optimization functions. e results shown in the study [37] could demonstrate the role of the solution. However, it should be noted that many of optimal solutions in considered benchmark functions have "0" values while "0" values are the middle points between the lower bound and upper bound. e functions that have zero optimal solution are Sphere, Rosenbrock, Rastrigin, Sum square, etc. So, center solutions in groups tend to contain variables close to 0 and these values are very useful in generating an updated step size. In the considered problem, the optimal solutions are not formed from middle points of lower bound and upper bound. So, the role of the center solution is no longer significant resulting in good quality solutions. Consequently, we propose to replace the center solutions with the best solution in the current population. As a result, the new solution generation strategy is modified as follows: By using the new strategy, the proposed method can reduce several shortcomings: (1) Reducing the task of sorting each variable from all solutions: as pointed in the COA method, each solution has N dv decision variables and the x th decision variable of all solutions must be grouped and then sorted. Finally, equation (23) is applied for determining the most appropriate x th decision variable. e task becomes time-consuming if N dv is high. Furthermore, if the number of groups is high, the number of center solutions is also high, resulting in long execution time for the center solutions.
(2) Finding higher quality solutions: the best solution has a high contribution to the level of solution quality improvement, so the proposed method can help to reduce the number of iterations and shorten execution time.
e performance of the first modification will be investigated by comparing the original COA and COA with the first modification, called ICOA1.

4.4.
e Second Improvement in the Second Generation. In the second improvement, we suggest to modify the strategy of producing new solutions for the second generation. As observing equation (26), it can be indicated that new variables in the sole solution are derived from randomization. e number of decision variables N dv is used in order to establish three selection conditions in which the condition is random and selected variables are also random. Furthermore, random variables within upper bound and lower bound are produced for the third selection condition. In the first two conditions, two variables are taken from two random solutions in the current population while random variables within predetermined boundaries are used in the third condition. Clearly, the combination of the variables cannot form good solutions with high quality. So, in the second improvement, we propose two models for updating new solutions based on the following equations: As seen from the two equations above, (31) can generate a larger step and new solution found by using (31) has a clearer difference. However, it is not sure that equation (31) can produce solutions with better fitness if there is no measurement technique. Let suppose that if the first step size between Co Gbest and Co best,1 and the second step size between Co Gbest and Co best,2 are small, the new solution obtained by equation (30) is so close to the old solution and local optimal solution search strategy may be ineffective. Furthermore, if the phenomenon is repeated in equation (31), the third step size will be high and the new solution is fallen into other zones, which are very far from the old zone. So, if the use is appropriate, equation (31) can prevent the effectiveness of the local search strategy. e shortcoming normally happens once the distances of considered solutions are high at the first some iterations. On the contrary, as considered solutions are very close together at high iterations or at the last iterations, equation (31) is really more effective. So, it should have a proper condition for determining the use of (31) or (30). In order to determine a more Complexity appropriate equation, we suggest that all solutions in each group should be evaluated based on fitness function. If many solutions are close together, a large step is selected. Otherwise, smaller step can be a better choice. At the beginning, we count the number of pairs of solutions that they are close together and then the number divides the maximum number of pairs of solutions. If the ratio is large enough, it means that many solutions are close together and the jumping step should be obtained by using (31). e second improvement can be seen in Algorithm 1. Here, the number of pairs of close solutions is symbolized by Count and determined by using a pseudocode while the maximum number of pairs of close solutions and the ratio can be, respectively, symbolized by Count max and R, which are obtained by equations (32) and (33). As indicated by equation (33), R is within 0 and 1, and its value is directly proportional to the number of pairs of close solutions. If R is close to zero, it means almost all solutions are not close but if R is around 1, approximately the considered solutions are very close together. So, it needs to determine a probability P a for each specific optimization problem. As P a is large enough, equation (31) providing large step size is applied. Otherwise, equation (30) is preferred. In numerical result, we will test the sensitiveness of P a on the obtained result and the most appropriate value of P a will be recommended.

Selection of Decision Variables.
In order to determine the most suitable decision variables for solutions, it needs to analyze data of the problem as well as obtained results. In addition to the use of optimization tools for determining site and rated power of PVDGUs, another technique, FW/ BWST, is also used to compute power flows. In this study, produced harmonic sources have a spectrum as shown in Table 1. Harmonic spectrum contains parameters of 5 harmonic orders with magnitudes (%) and phase angles (°). ese five harmonic flows are injected into the linear loads and treat these loads as harmonic generation sources. e linear loads that are injected by harmonics will be considered as nonlinear loads. Specifically, there are 5 and 10 nonlinear loads in the IEEE 33-bus distribution network and IEEE 69bus distribution network, respectively. Due to the characteristics of harmonics, only the odd-order harmonics are considered because the even-order harmonics are ignored because of having very small amplitude. For each different order, magnitude (%) and phase angle (°) are different in harmonic distortions. After harmonic current with amplitude and phase angle are injected into loads in the system, it is necessary to recalculate the power flow because the harmonic current affects the entire system. In fact, the current source of harmonics is highly dependent on the presence of transformers and line impedances in the power distribution system [48]. e value of harmonic voltage and harmonic current will be found by solving power flow at each harmonic order. Results of voltage and current will be compared with the IEEE Std. 519 for the evaluation of power quality. In this paper, FW/BWST is selected as a technique for calculating the power flow at both fundamental frequency and higher order frequency. Basically, input data of FW/BWST consist of line parameters (resistance, reactance, and maximum current), load parameters (active power demand and reactive power demand), PVDGU parameters (site and rated power), and the harmonic parameters (order, current magnitude, and current phase angle). e obtained parameters after running FW/ BWST are line currents and bus voltage for each harmonic order. As a result, total power losses in all lines together with THD and IHD are calculated by using equations (1), (8), and (9). Among input data, only site and rated power of PVDGUs are unknown parameters. So, the two factors are included in solutions in which the site and the rated power of the j th solar photovoltaic distribution generation unit (PVDGU) corresponding to the k th solution in the g th group are, respectively, represented by Site j,k,g and RP j,k,g . Each solution in each group is represented by the following equation: where N DG is the number of installed PVDGUs in the distribution network.
In order to produce the initial population, lower bound solution Co min and upper bound solution Co max should be defined to make sure that randomly generated solutions can satisfy the limit constraints: where Site min j and Site max j are, respectively, bus 2 and bus N bus and RP min j and RP max j are the minimum rated power and maximum rated power of the j th PVDG unit. However, the minimum power of all units and the maximum power of all units are the same and corresponding to RP min and RP max . e limits of site are easily selected based on the system dimension, but the limits of rated power must be selected thoroughly and exactly. e limits will be discussed in the numerical results. Each solution in the initial population is formed by the following way: Co k,g � Co min + ε · Co max − Co min ; k � 1, . . . , N Co ; g � 1, . . . , N g . (36) It is noted that between site and rated power, site should an integer value but rated power can be a continuous value. So, site is approximated by using command "round" while rated power remains unchanged. e task continues to be done after producing new solutions for the first and the second generations. 10 Complexity

Penalty Terms for Violating of Dependent Variables.
After having decision variables, FW/BWST technique has been run for obtaining dependent variables including branch currents I h nDG and load bus voltage V h i where h is from 1 to H. In the study, we suppose there are 5 harmonic currents consisting of orders 5, 7, 11, 13, and 17. So, h is, respectively, equal to 1, 5, 7, 11, 13, and 17. en, current in branch n, I nDG , and voltage at bus i, V i , are calculated by using the two following equations: As a result, the two parameters are checked and penalized by using the following formulas: In addition to taking voltage and current into consideration, THD and IHD are also constrained for better power quality. THD and IHD are calculated using equations (8) and (9). en, they are checked and penalized as follows: 5.3. Fitness Function. As described in Section 2, we focus on four single objectives consisting of total power loss minimization F 1 , rated power minimization of PVDGUs F 2 , load bus voltage enhancement F 3 , and harmonic distortion minimization F 4 . According to the considered objectives, fitness function can be established as follows: (i) Fitness function for total power loss minimization: (ii) Fitness function for rated power minimization of PVDGUs: (iv) Fitness function for harmonic distortion minimization:

Producing New Solutions and Fixing Violated
Variables. e proposed ICOA experiences two new solution generations. e first generation employs equation (29) to produce N Co new solutions for each group while the second generation applies Algorithm 1 for producing one new solution for each group. After producing new solutions (i.e., new site and new rated power for each PVDGU), the site and the rated power are checked and adjusted by the function of the two following equations: Between site and rated power, site is approximated for getting integer value by using command "round" while rated power remains unchanged.

Termination of Running ICOA and Evaluation of Obtained Results.
In order to evaluate the performance of the proposed ICOA method, we execute 50 independent trial runs for each study case and the best fitness of each run is recorded. Some main factors of fifty fitness functions such as minimum, average, and maximum are calculated and reported for performance comparison with other methods. For executing one trial run, the proposed ICOA carries out IT max iterations in which IT max is the maximum number of iterations. So, the best run and the fifty runs are recorded and plotted in figures.

5.6.
e Computation Steps of Using the Proposed ICOA Method for Solving the Considered Problem. Main stages of the search process by using ICOA for solving PVDGU installation problem have been explained in the previous section and the detail of the implementation is shown in Figure 1.

Simulation Results
In order to investigate the real performance of the proposed ICOA method in finding the optimal location and sizing of PVDGUs, IEEE 33-bus distribution power network and IEEE 69-bus distribution power network [26,49] are used as two main studied tests. e IEEE 33-bus distribution power network shown in Figure 2 has a total load demand of 3715 kW and 2300 kVAR and operates at a nominal voltage of 12.66 kV. e IEEE 69-bus distribution power network shown in Figure 3 operates at a nominal voltage of 12.66 kV and supplies power to total load demand of 3802 kW and 2694 kVAR. e two systems are considering linear loads. So, in order to produce negative impacts of harmonics, we inject the harmonic currents with five harmonic spectrums simultaneously to the IEEE 33-bus distribution power network at 6 load buses 10, 15, 20, 24, 27, and 32 and to the IEEE 69-bus distribution power network at 10 load buses 10,13,18,19,22,25,34,46, 56, and 65. e detailed information of the five harmonics is shown in Table 1 [50].  e capacity for each PVDGU can be defined as a variable varying from 0 to 2 MW [51]. Two limit points 0 and 2 MW are treated as lower and upper bounds of rated power for each PVDGU, respectively. Besides the capacity, the position that we can locate each PVDGU is from the smallest bus number (bus 2) to the highest bus number (bus 33 for the first system and bus 69 for the second system). Due to the total load demand as well as the system dimension, three PVDGUs are installed in the IEEE 33-bus distribution network whereas four PVDGUs are installed in the IEEE 69bus distribution network. e optimal location and optimal capacity of PVDGUs in the two systems aim at reaching four single-objective functions considered in the four following cases. In addition to ICOA, eight other methods consisting of BBO, GA, PSO, SFO, SSA, COA, ICOA1, and ICOA2 are also implemented for dealing with the four cases above for the two power systems. For each case, each method is run 50 trial times and the fitness functions of the trial runs are recorded for reporting the best fitness, average fitness, and maximum fitness in aim to support evaluation of performance. MATLAB program language and a personal computer with a processor of 2.0 GHz and 2.0 GB RAM are the main tools for executing simulation results. Before running these methods, their control parameters are analyzed and selected based on experiments as well as referred from previous studies such as [1,10,11] for GA, [17,18,44] for PSO, [23,52] for BBO, [53] for SFO, [54] for SSA, and [37] for COA. e process of tuning the parameters is explained as follows: (1) For BBO, its five control parameters are selected as follows: maximum immigration (I) and migration rates (E) for each habitat are tried in the range from 0.2 to 1 with a step of 0.2; maximum mutation rate (m max ) is set to five values from 0.1 to 0.5; habitat modification probability (P mod ) is equal to 1; and combining ratio of the selected habitat from the previous generation and current generation for next generation for each iteration is set to the range 20%-80% [26]. (2) For running GA, the optimal parameters are used as follows: the crossover probability (Cr) is set to three   are set to 2 [44]. (4) To implement SFO, the day (d) and the sun (s) are, respectively, fixed at 100 and 1 while the pollination (p) is ranged from 0.6 to 1 with a step of 0.2 [53]. (5) For SSA, c 1 is a function of 2e − (4IT/IT max ) 2 ; c 2 and c 3 are randomly produced numbers in the interval of [0, 1] [54]. (6) In order to run COA, ICOA1, ICOA2, and ICOA methods, the number of coyotes in each group N Co and the number of groups are set to 4. According to the settings, the population N pop of the four methods is 20. In addition, ICOA2 and ICOA have more than one advanced control parameter P a and the parameter is set to five values from 0.2 to 1.0 with a step of 0.2. e impact of P a on simulation results will be surveyed and discussed.
On the other hand, the applied methods have two common control parameters: they are population N pop and the maximum number of iteration IT max . N pop of PSO, GA, BBO, SSA, and SFO is also set to that of COA methods, i.e., 20, while IT max has the same values for all methods. IT max is 75 for the IEEE 33-bus distribution network and 100 for the IEEE 69-bus distribution network.

Survey on Results Obtained by Setting Different Values of P a for ICOA Method.
For finding the most suitable P a value, ICOA is in turn run 50 times corresponding to each P a value by setting P a to 0.2, 0.4, 0.6, 0.8, and 1.0. e second case of minimizing the rated power of PVDGUs is selected for both IEEE 33-bus distribution network and IEEE 69-bus distribution network. e obtained results for the two systems are shown in Tables 2 and 3, respectively. As observed from Table 2, the best fitness corresponding to P a � 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0 is 0.2860, 0.2836, 0.2978, 0.2981, 0.2900, and 0.2978, respectively. rough the results, it points out that ICOA method can reach the best solution only at P a � 0.2. Additionally, the average values of fitness function also confirm the most effective impact of P a � 0.2 on its average fitness is the smallest among five values. e best fitness can conclude P a � 0.2 is the most suitable value for finding the best optimal solution while the best average fitness can confirm the most stable search ability of ICOA if setting P a � 0.2. One more time, the significance of setting P a � 0.2 can be adopted since the best fitness and the average fitness at P a � 0.2 continue to be the smallest values among five values for the second system shown in Table 3. e best fitness and the best average are, respectively, 0.2425 and 0.2939. Consequently, we accept the best P a value is 0.2 and it is used for other remaining cases with intent to shorten simulation time.

Comparison Criteria.
In order to demonstrate the effective search performance of the proposed method with the other methods, the following five criteria have been applied: (1) e best solution of the proposed method and compared methods was compared to show the strong search ability of the proposed ICOA. (2) e average fitness function of 50 trial runs and fitness function of 50 trial runs were compared to evaluate the solution search stability and the oscillation of the proposed ICOA. (3) e number of better runs (N betterrun ) that the proposed ICOA method can reach was compared to the best run of compared methods. (4) Improvement level (IL) of the best fitness, average fitness, and worst fitness from the proposed method over other compared methods is exactly computed. e three factors are obtained by using the following equations: IL of the best fitness(%) � the best fitness of the compared method − the best fitness of ICOA the best fitness of the compared method × 100, IL of average fitness(%) � average fitness of the compared method − average fitness of ICOA average fitness of the compared method × 100, IL of the worst fitness(%) � the worst fitness of the compared method − the worst fitness of ICOA the worst fitness of the compared method × 100.  e implemented methods in this study are based on metaheuristic algorithms. erefore, determining the exact number to show the efficiency of the proposed method and other methods is difficult to accomplish the evaluation [44,45]. However, the evaluation of methods based on criteria which are related to setting parameters and solution quality is always appreciated. Details for evaluating the proposed method and other methods have been analyzed in the section below.

Test System 1: IEEE 33-Bus Distribution Network
Case 1: Minimizing the power losses of the system In this case, total power losses in all branches are considered to be minimized while load voltage profile, branch currents, and total harmonic distortion and individual harmonic distortion are constrained. Fitness function FF 1 is used to determine the quality of solutions from the applied method. Summary of results from 50 trial runs consisting of the best, average, and worst fitness values is reported in Table 4.
As the shown results, the best and average fitness values found by COA, ICOA1, ICOA2, and ICOA are 0.3462, 0.3461, 0.3462, and 0.3461, and 0.3727, 0.3674, 0.3779, and 0.3636. e best values indicate that ICOA1 and ICOA can find the same best fitness less than that of COA and ICOA2. e average value of ICOA1 is better than that of COA but that of ICOA2 is worse than that of COA.
us, it can conclude that the first improvement is really effective for ICOA1, but the second improvement is not effective for ICOA2. Clearly, the individual application of the second improvement does not reach high performance; however, the combination of the two improvements is much more effective for ICOA. And ICOA is the best one among four methods as a result. Figure 4 showing fifty fitness values of the four COA methods can confirm the superiority of ICOA. Because the ICOA has many searched fitness values that are significantly lower than the COA.
As compared to BBO, GA, PSO, SFO, and SSA, the proposed ICOA shows its potential search ability since its best fitness and average fitness are the lowest values among six methods. e improvement level of the best fitness, average fitness, and the worst fitness from ICOA over other methods is significant. In fact, the value is from 0.115% to 11.028% for the best fitness, from 3.860% to 27.671% for average fitness, and 4.239% to 35.260% for the worst fitness.
ere is a coincidence that the second best method and the worst method for the three comparisons are BBO and GA, respectively. In addition, the improvement of the performance of the proposed ICOA over three other methods, PSO, SFO, and SSA, is also high with the improvement level around 10%. Figures 4 and 5 are plotted for illustrating 50 fitness function values of 50 runs obtained by these methods and ICOA. As counting from the figure, ICOA has approximately twenty points with the same fitness functions and below all points of other methods. Furthermore, the oscillations of ICOA are very narrow and the deviation between the highest point and the lowest point is small whereas the oscillations and the deviation of other methods are much higher, especially those from GA. If considering the number of runs between ICOA and other ones, ICOA has 13 runs with better fitness than the best fitness of COA, 0.3462, while the number is much higher as compared to other ones. In fact, it is 15,45,27,38, and 27 runs with better fitness than the best solution of BBO, GA, PSO, SFO, and SSA, respectively. us, it can conclude that the use of the proposed ICOA method for finding the optimal site and capacity of PVDGUs is the most effective among nine methods. Optimal solutions together with IHD, THD, voltage profile, and total power losses obtained by these methods are shown in Table 5.

Case 2: Minimizing PVDG penetration level in the system
In the second case, the sum of the rated power of all PVDGUs is used to be an objective function whereas load voltage profile, branch currents, and total harmonic distortion and individual harmonic distortion are taken into consideration. Fitness function FF 2 in equation (41) is used to evaluate the effectiveness of solutions from implemented approaches. e best, average, and worst fitness values from 50 runs obtained by nine methods together with the improvement level of the ICOA over other methods are reported in Table 6. For clear observation, 50 runs from all methods are separated into two figures in which Figure 6 shows all runs from four COA methods and Figure 7 shows all runs from GA, BBO, PSO, SSA, SFO, and ICOA. For comparison among four COA methods, reported values are good evidences for demonstrating the effectiveness of ICOA over three remaining ones since its improvement levels are much higher than those in Case 1. ICOA can improve the best fitness, average fitness, and the worst fitness up to 5.718%, 6.138%, and 7.621%. ICOA1 continues to be the second best method, but ICOA2 is still the worst method. Clearly, the first proposed improvement still has better contribution to the performance of ICOA. e second proposed improvement is not effective if it is carried out independently, but it can enable ICOA to find better solutions than ICOA1. Figure 6 can support better view of the outstanding performance of ICOA over COA. Consequently, it can conclude that the two proposed improvements should be combined in aim to improve the performance of the COA method.
Comparison results between ICOA and other methods are also similar to those in Case 1 when BBO is still the second best method and GA is still the worst method for average and the worst fitness comparison. SFO is the worst method for the best fitness comparison. e improvement level of ICOA over these methods is also high. It is from about 6% up to 20% for the best fitness function improvement level. e average and the worst fitness improvement levels are, respectively, up to about 36% and 41%. Figure 7 shows that about over 50% points of ICOA are below the points of other methods. Besides, as shown in the Complexity 15 last column of Table 6, ICOA has many runs with better fitness than the best one of other methods such as 5 runs more than BBO, 41 runs more than GA, 16 runs more than PSO, 43 runs with better SFO, and 39 runs more than SSA. Hence, the conclusion based on the numerical results as well as graphics results in Figure 7 is that ICOA is superior to BBO, GA, PSO, SFO, and SSA methods for the second case of the IEEE 33-bus distribution network. Optimal solutions together with IHD, THD, voltage profile, and total rated power of all PVDGUs obtained by these methods are shown in Table 7.
Case 3: Improving the voltage profile In the part, voltage of all loads is considered to be improved by using fitness function FF 3 in equation (42) as a comparison criterion for evaluating applied methods. In the fitness function, objective function F 3 in equation (7) together with branch currents, voltage limits, and total harmonic distortion and individual harmonic distortion is taken into consideration.
Results for comparison are reported in Table 8. e best fitness function of 0.0937 can point out that the proposed ICOA is the best method among compared methods in finding the optimal solution. e lowest average fitness function and the lowest worst fitness function confirm that the proposed method is the most stable tool with the lowest fluctuations. e improvement level of the best fitness is from 0.531% to 17.590%. e improvement levels of the average fitness and the worst fitness are also high, up to 49.739% and 53.211%, respectively. Besides, the number of runs from ICOA with better fitness than the best one of other methods reported in the last column of Table 8 indicates that ICOA has many runs with better fitness than the best one of other methods such as 4 runs more than BBO, 28 runs more than GA, 5 runs more than PSO, 5 runs with better SFO, and 15 runs more than SSA. Figure 8 shows that many points of ICOA are below those of other compared methods and the oscillations of ICOA are much lower than those of other ones. So, the probability that the proposed method finds the optimal solution is much higher than that of other methods.
Optimal solutions together with IHD, THD, voltage profile, and total rated power of all PVDGUs obtained by these methods are shown in Table 9.

Case 4: Reducing the harmonic distortions
In the part, THD and IHD of all buses are improved by using fitness function FF 4 in equation (43) to compare the quality of obtained solutions from implemented methods. In the fitness function, objective function F 4 in equation (12) and branch currents and voltage limits are taken into consideration. From comparing the relation of two equations (8) and (9), THD value depends on IHD value. When     IHD value changes, THD value will also change proportionately. So, to balance between these two values, ω 1 and ω 2 are set to 0.5 and 0.5, respectively. Results reported in Table 10 show that the best optimal solution of the proposed method has approximately equal fitness or not much better fitness than that of other ones since the improvement level is approximately 0% to 4.066%. ICOA1 and ICOA have the same best solution, but all solutions on average of the proposed ICOA are better than those from ICOA1 through average fitness improvement level. In addition, the fluctuations of the proposed method are also smaller. e improvement levels of the average fitness and the best fitness is much better and can be from 0.431% to 21.294% and from 3.503% to 17.133%. e proposed method has 4 runs and 45 runs with better fitness than the best fitness of COA and GA, respectively. Figure 9 shows that many of points of ICOA have lower fitness than those of other methods and many points from ICOA have the same fitness as the best point while points from other methods are far from the best point. So, it can conclude that the proposed method is very effective for the case.

Complexity
Optimal solutions together with IHD, THD, voltage profile, and total rated power of all PVDGUs obtained by these methods are shown in Table 11. In general, the proposed ICOA method outperforms the original COA method for all cases with high improvement level. e improvement level of the best fitness is 0.264% for the first case, 0.082% for the second case, 15.575% for the third case, and 0.293% for the last case. In spite of only the third case with significant improvement, the proposed method has a stronger search ability of finding optimal solutions. In fact, the best solutions are found by the proposed method, but COA does not approach them with the same number of iterations. Furthermore, the probability of finding the best solutions or ones close to the best solutions from the proposed method is very high, but that of COA is very low or nearly zero. For example, as we count the runs    18 Complexity with a fitness function of 0.3027 for the first case, COA has only one run. Similarly, COA has four runs with a fitness function of 0.2427 for the second case, two runs with fitness function of 0.0565 for the third case, and one run with fitness of 0.9547 for the last case. In this regard, the proposed method has nine runs and three runs with better fitness than the best fitness of COA for cases 1 and 2, respectively. For cases 3 and 4, ICOA has four runs with better fitness than the best run of COA. In addition, the stability and oscillation over 50 runs obtained by ICOA and COA can be evaluated via the improvement levels of the average fitness and the worst fitness. It is clear that the values are much higher than the improvement level of the best fitness. e two values are 0.583% and 1.843% for Case 1, 5.711% and 11.674% for Case 2, 10.329% and 11.531% for Case 3, and 0.309% and 1.178% for Case 4. As compared to other remaining methods excluding COA methods, the outstanding performance of ICOA is clearer since the compared values are much higher. For example, the improvement level of the best fitness is from 0.264% to 3. Clearly, most runs of ICOA have lower fitness and smaller oscillations than those from other methods. From the above evidences, it is clearly pointing that the proposed method has a better performance than the other methods in both optimal solution and stability. IHD, THD, voltage profile, and total rated power of all PVDGUs obtained from optimal solutions by implemented methods of the IEEE 69-bus distribution network for four cases are, respectively, shown in Tables 16-19.

Discussions
As analyzed, previous studies mostly focused on considering the power balance constraints, voltage limit, branch current limit, capacity limit of each PVDGU, and limit of total installed capacity of all PVDGs. However, those studies did not pay attention to another important factor, which is harmonic distortion. erefore, this study considered all constraints as in the constraints section, including THD and IHD limits. In addition, in order to prove the superiority in terms of the best solution, the stability of search ability, and fast search speed, the proposed method has been compared with different methods such as BBO, GA, PSO, SFO, SSA, COA, ICOA1, and ICOA2. To get an accurate evaluation, all methods have been implemented in different scenarios with four different objective functions for both the IEEE 33-bus distribution network and IEEE 69-bus distribution network. ese computational efforts make it have much objectivity in evaluating the obtained results from the methods. Results from four cases of two systems and comments are as follows: (i) ICOA1 has found less power loss, smaller capacity of all PVDGUs, better voltage profile, and smaller THD and IHD than COA. Recall that ICOA1 has been developed replacing the central solution Co mid,g in equation (22) with the best solution   shown in (29). us, it is suggested the modification should be applied for COA. (ii) ICOA2 could not reach less power loss, smaller capacity of all PVDGUs, better voltage profile, and    smaller THD and IHD than COA even its performance was worse than COA for some cases. Recall that ICOA2 was developed by replacing the combining three randomizations of COA in equation (26) by Algorithm 1. us, the modification should not be implemented for COA if there is not another modification with it. In fact, the algorithm was really useful for ICOA, which combined ICOA1 and ICOA2.
(iii) e proposed ICOA was more effective than ICOA and much more effective than BBO, GA, PSO, SFO, and SSA. As compared to COA, ICOA has reached less values of minimum fitness, average fitness, and worst fitness. In addition, almost all runs of ICOA had less fitness than those of COA and ICOA has found high number of solutions with better quality than the best solution of COA. Furthermore, the                 improvement of ICOA over other ones was more significant than that over COA.
However, to reach the high performance of ICOA, there must be concerns during the implementation process. Like other methods such as GA, PSO, BBO, SFO, and SSA, the quality of the solutions from the proposed method depends heavily on control parameters. e ICOA method has three basic parameters including N Co , N g , and IT Max . In this method, the population size (N pop ) and the number of new solutions for each iteration are equal to (N Co × N g ) and (N Co × N g + N g ), respectively. During the simulation process, we surveyed and selected two parameters as N Co and N g to find the most appropriate values for two considered systems. e values of N Co and N g are adjusted to find the most suitable values. As observing obtained results, we have found that 4 is suitable for both N Co and N g whereas 75 and 100 are, respectively, reasonable for IT Max for the first system and the second system. erefore, it is not easy to choose the right parameters for each specific network. is takes a long time to test and evaluate the quality of the solution along with the different settings. In addition, in order to obtain objective and accurate results in comparing the proposed solution with other methods, 50 trial runs with the randomly generated initial solution were performed for all four cases including reducing the power losses of the system, minimizing the PVDG penetration level in the system, improving the voltage profile, and mitigating the harmonic distortions. In all four cases, the capacity of PVDGUs variables is in the range from 0 MW to 2 MW and the position variables are from the smallest bus number (except the slack bus) to the highest bus number. As obtained result, the proposed method always gives the optimal solution better than compared methods in 4 cases with two different distribution systems. Finally, the optimal solution (location and sizing of PVDGUs) which is found by the proposed method can bring more economic and technical benefits than other methods as compared.

Conclusions
In this paper, a proposed coyote optimization algorithm has been developed to find the optimal location and sizing of PVDGUs in the radial distribution systems by considering four optimization cases consisting of total power losses, capacity of all PVDGUs, load voltage profile, and harmonic distortions. In addition, constraints from distribution systems such as current limits of conductors, voltage limits of all loads, and limits of total harmonic distortion and individual harmonic distortion have been always taken into consideration seriously. Two standard study cases, IEEE 33bus distribution network and IEEE 69-bus distribution network, have been employed and nonlinear loads have been created by simultaneously injecting five harmonic currents into five loads in the first system and ten loads in the second system. e obtained results show that the proposed method (ICOA) can produce higher quality solutions with better stability than the other methods while the technical criteria are still satisfied. e result comparisons are as follows: (i) As compared to COA, the proposed ICOA could reach the improvement of performance up to 5.718% and 15.575%, and the improvement of the search process stability up to 6.138% and 10.329% for the first system and the second system, respectively (ii) As compared to other methods, the proposed ICOA could reach the improvement of performance up to 23.393% and 57.487%, and the improvement of the search process stability up to 27.671% and 72.908% for the first system and the second system, respectively (iii) For the two systems, the proposed ICOA could find 9 and 13 solutions better than the best solution of COA, and 45 and 48 solutions better than the best solutions of other ones Connecting appropriate DGPVUs in the distribution system can bring both economic and technical benefits. As compared to other methods, the proposed method can (i) Reduce the power loss by 9.002 kW and 2.547 kW and minimize the penetration level of PVDGUs by 0.3398 MW and 0.4311 MW for the first and second systems, respectively (ii) Improve voltage and reduce harmonic distortion by 0.0027 pu and 0.289% for the first system and by 0.019 pu and 0.190% for the second system, respectively Consequently, it can lead to a suggestion that the proposed method should be applied for finding the optimal location and rated power of PVDGUs in distribution power systems for the purposes of economic and technical issues.

Co new
g : e new solution in the g th pack Co best,g and Co worst,g : e best solution and the worst solution in the g th group Co best,1 , Co best,2 , Co best,3 , and Co best,4 : e best solutions in all packs picked up randomly Co Gbest : e best solution in the population Co k : e k th solution in the g th group Co 1,g and Co 2,g : Two randomly picked solutions in the g th pack ΔI n : Current difference at the n th branch ΔIHD i : Individual harmonic distortion difference at the i th bus ΔTHD i : Total harmonic distortion difference at the i th bus ΔV i : Voltage difference at the i th bus F 1 : Objective function of total power losses F 2 : Objective function of the PVDGU penetration level F 3 : Objective function of voltage profile index F 4 : Objective function of harmonic distortion FF new k,g and FF k,g : Fitness function of the k th new solution and the k th old solution in the g th group FF 1 : Fitness function for total power loss minimization FF 2 : Fitness function for rated power minimization FF 3 : Fitness function for improving bus voltage FF 4 : Fitness function for minimizing harmonic distortions FF best,g : Fitness function of the best solution in the g th group FF i,g and FF j,g : Fitness function of the i th and the j th solutions in the g th group FF mean,g : Average fitness function of all solutions in the g th group H: Maximum order frequency I h nDG : Maximum current magnitude in the n th branch at the h th order frequency IHD h V,i : Voltage individual harmonic distortion of the h th order harmonic at the i th bus I n : Current magnitude in the n th branch without PVDGUs I nDG : Current magnitude in the n th branch with PVDGUs I max n : Maximum current magnitude in the n th branch IT: Current iteration IT max : Maximum of iteration N betterruns : Number Active power of the i th load P Loss,n : Active power losses at the n th branch R n : Resistance of the n th branch THD V,i : Voltage harmonic distortion at the i th bus THD max : Maximum total harmonic distortion THD Vaver : Total harmonic distortion average value of the system TPL: Total power losses of all branches without PVDGUs TPL DG : Total power losses of all branches with PVDGUs μ I : Penalty factor for violated current μ V : Penalty factor for violated voltage μ THD : Penalty factor for violated total harmonic distortion μ IHD : Penalty factor for violated individual harmonic distortion V h i : Voltage magnitude of the i th bus at the h th order frequency V min and V max : Lower and upper limitations of load bus voltage magnitude VPI: Voltage profile index of the system.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.