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

Authors: Thang Trung Nguyen; Thai Dinh Pham; Le Chi Kien; Le Van Dai
Journal: Complexity
 (2020)

---

## Abstract

This paper proposes an improved coyote optimization algorithm (ICOA) for optimizing the location and sizing of solar photovoltaic distribution generation units (PVDGUs) in radial distribution systems. In the considered problem, four single objectives consisting of total power losses, capacity of all PVDGUs, voltage profile index, and harmonic distortions are minimized independently while satisfying branch current limits, voltage limits, and harmonic distortion limits exactly and simultaneously. The performance of the proposed ICOA method has been improved significantly since two improvements were carried out on the two new solution generations of the conventional coyote optimization algorithm (COA). By finding four single objectives from two IEEE distribution power systems with 33 buses and 69 buses, the impact of each proposed improvement and two proposed improvements on the real performance of ICOA has been investigated. ICOA was superior to COA in terms of capability of finding higher quality solutions, more stable search ability, and faster convergence speed. Furthermore, we have also applied five other metaheuristic algorithms consisting of biogeography-based optimization (BBO), genetic algorithm (GA), particle swarm optimization algorithm (PSO), sunflower optimization (SFO), and salp swarm algorithm (SSA) for dealing with the same problem and evaluating further performance of ICOA. The result comparisons have also indicated the outstanding performance of ICOA because it could find much better results than these methods, especially SFO, SSA, and GA. Consequently, the proposed ICOA is a very effective method for finding the optimal location and capacity of PVDGUs in radial distribution power systems.


---

## Body

## 1. Introduction

### 1.1. 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 high-voltage 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–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]. Therefore, 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].


### 1.2. 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]. These 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. The 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. This 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. This is a good method, but the biggest disadvantage of this method is the complicated calculation process. Thus, 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. This 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. Therefore, 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. The result was very satisfactory by installing 3 DGUs in the IEEE 33-bus 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. The 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. Therefore, 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. Therefore, its convergence to a global optimum is hardly ever guaranteed. There was not enough evidence to conclude the performance of PSO because only the power loss objective of the IEEE 33-bus 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. The algorithm used to build in this article is complex but effective. The main component of the algorithm consists of four pillars as social structure, opinion space, social influence, and update rule. The 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–23] have tried to reduce power loss and operation cost and improve voltage, but they have ignored one important 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. The 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. The study has considered power loss minimization to be a priority while THD and IHD have been constrained not higher than 5% and 3%, respectively. The 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. There 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.


### 1.3. Proposals, Novelties, and Contributions

Realistically, studies in [24–26] have considered THD and IHD and the two factors have been reduced not higher than 5% for THD and 3% for IHD. The two factors always satisfied IEEE Std. 519 [26]. However, the studies in [24–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–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. Thus, 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. Thus, 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. The four single objectives together with constraints of technical factors are considered in the study by the four following cases:Case 1: minimizing the power losses of the system under the constraints of harmonic distortions, voltage profile, and branch currentCase 2: minimizing the rated power of PVDGUs under the constraints of load bus voltage, branch current, and harmonic distortionsCase 3: improving the voltage profile under the constraints of harmonic distortions and branch currentCase 4: reducing the harmonic distortions under the constraint of voltage profile and branch currentIn the above four cases, the different harmonic sources are injected into linear loads in two considered systems and make them become nonlinear loads. Then, 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. The voltage and current are directly related to total power losses, voltage profile, and harmonic distortion. The 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. The authors in [32] have applied COA to minimize the gas consumption of turbines in combined cycle power plants in Brazil. The 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), Self-adaptive 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. The 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. Thereby, the effectiveness of COA was not really convincing in many different problems. This 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. The 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). The 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. The coyote community is divided into Ng small coyote groups, and there are Nco coyotes in each group. COA method produces two new solution generations per iteration in which the first generation produces Nco new solutions for each group and the second generation produces Ng new solutions for the whole coyote community. Thus, total number of new solutions generated in the COA method is (Nco × Ng + Ng) solutions in which Nco × Ng is equal to population (Npop). Thus, (Nco × Ng + Ng) is equal to (Npop + Ng). The 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 Nco 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. The 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. The second modification can enhance the effectiveness of local search and find one promising solution for each group. The 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. The obtained results from ICOA in four cases are compared with those of other methods for performance evaluation.In summary, the main contributions of this study are as follows:(1)
The 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. The benefits help the readers to easily identify the most appropriate single objective of installing PVDGUs for their purpose.(2)
The 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)
The results show that ICOA is more suitable than BBO, PSO, GA, SFO, and SSA for installing optimal PVDGUs in distribution networks. The reader should try other methods for the problem.The 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 Section2. The whole search procedure of the coyote optimization algorithm (COA) and improved coyote optimization algorithm (ICOA) is described in detail in Section 4. The 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.


## 1.1. 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 high-voltage 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–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]. Therefore, 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].


## 1.2. 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]. These 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. The 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. This 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. This is a good method, but the biggest disadvantage of this method is the complicated calculation process. Thus, 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. This 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. Therefore, 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. The result was very satisfactory by installing 3 DGUs in the IEEE 33-bus 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. The 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. Therefore, 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. Therefore, its convergence to a global optimum is hardly ever guaranteed. There was not enough evidence to conclude the performance of PSO because only the power loss objective of the IEEE 33-bus 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. The algorithm used to build in this article is complex but effective. The main component of the algorithm consists of four pillars as social structure, opinion space, social influence, and update rule. The 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–23] have tried to reduce power loss and operation cost and improve voltage, but they have ignored one important 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. The 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. The study has considered power loss minimization to be a priority while THD and IHD have been constrained not higher than 5% and 3%, respectively. The 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. There 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.


## 1.3. Proposals, Novelties, and Contributions

Realistically, studies in [24–26] have considered THD and IHD and the two factors have been reduced not higher than 5% for THD and 3% for IHD. The two factors always satisfied IEEE Std. 519 [26]. However, the studies in [24–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–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. Thus, 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. Thus, 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. The four single objectives together with constraints of technical factors are considered in the study by the four following cases:Case 1: minimizing the power losses of the system under the constraints of harmonic distortions, voltage profile, and branch currentCase 2: minimizing the rated power of PVDGUs under the constraints of load bus voltage, branch current, and harmonic distortionsCase 3: improving the voltage profile under the constraints of harmonic distortions and branch currentCase 4: reducing the harmonic distortions under the constraint of voltage profile and branch currentIn the above four cases, the different harmonic sources are injected into linear loads in two considered systems and make them become nonlinear loads. Then, 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. The voltage and current are directly related to total power losses, voltage profile, and harmonic distortion. The 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. The authors in [32] have applied COA to minimize the gas consumption of turbines in combined cycle power plants in Brazil. The 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), Self-adaptive 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. The 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. Thereby, the effectiveness of COA was not really convincing in many different problems. This 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. The 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). The 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. The coyote community is divided into Ng small coyote groups, and there are Nco coyotes in each group. COA method produces two new solution generations per iteration in which the first generation produces Nco new solutions for each group and the second generation produces Ng new solutions for the whole coyote community. Thus, total number of new solutions generated in the COA method is (Nco × Ng + Ng) solutions in which Nco × Ng is equal to population (Npop). Thus, (Nco × Ng + Ng) is equal to (Npop + Ng). The 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 Nco 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. The 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. The second modification can enhance the effectiveness of local search and find one promising solution for each group. The 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. The obtained results from ICOA in four cases are compared with those of other methods for performance evaluation.In summary, the main contributions of this study are as follows:(1)
The 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. The benefits help the readers to easily identify the most appropriate single objective of installing PVDGUs for their purpose.(2)
The 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)
The results show that ICOA is more suitable than BBO, PSO, GA, SFO, and SSA for installing optimal PVDGUs in distribution networks. The reader should try other methods for the problem.The 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 Section2. The whole search procedure of the coyote optimization algorithm (COA) and improved coyote optimization algorithm (ICOA) is described in detail in Section 4. The 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.


## 2. 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]. Therefore, in the study, all PVDGUs connected in the systems have been checked for pure sinusoidal voltage waveforms by applying harmonic filters.


### 2.1. Assumption

The 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. The 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.(3)
Placement of solar photovoltaic distribution generation units: suppose that all buses (except for slack bus) in the considered system can be selected for installation.(4)
Lines: transmission capacity limits are predetermined for each branch and remain unchanged during working time. Parameters such as impedance represented the characteristics of the lines and are not affected by over temperature.(5)
Loads: this includes nonlinear loads and linear loads. 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. This 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.


### 2.2. 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 harmonic distortion (THD & IHD). The objective function along with the constraints is explained by using mathematical formulas as follows.


### 2.3. 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. Therefore, the reduction of total active power losses is becoming an important objective, which is shown in the following mathematical equation [41]:(1)TPLDG=∑n=1NbrInDG2Rn.By using the same method, total power losses in all branches before installing PVDGUs, TPL, is determined by means of the following equation [41]:(2)TPL=∑n=1NbrIn2Rn.It shows that if PVDGUs are properly installed, the losses of power system after installation are always smaller than before installation. Therefore, 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:(3)minimizeF1=TPLDGTPL.


### 2.4. The PVDGUs’ Penetration Level

Integration of large distributed generation into electric system can increase investment and operating costs. Therefore, minimizing the capacity of PVDGUs (RPDG) that still satisfies the technical criteria should be considered as the target. The proposed mathematical equation of the objective function is shown as follows [42]:(4)minimizeF2=∑i=1NDGRPDG,i∑i=1NloadPLoadi.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.


### 2.5. Voltage Profile Index

One of the important criteria for evaluating the stability of the operating grid is the voltage profile index (VPI). Thus, the voltage profile index must be considered and the calculation process for VPI should follow the step below [43].The voltage profile for theith bus can be defined as(5)VPi=Vi−Vmin⋅Vmax−ViVnom−Vmin⋅Vmax−Vnom,where Vnom is the nominal or desired bus voltage value and typically taken as 1.0 pu. It is clear that VPi can get a maximum value (1.0 pu) if the bus voltage equals to Vnom and VPi can get a negative value if the bus voltage is higher than the maximum limit or lower than the minimum limit. The overall voltage profile index for the distribution system can be defined as(6)maximum VPI=1Nbus⋅∑i=1NbusVPi.The 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:(7)minimizeF3=1−VPI.


### 2.6. Harmonic Distortions

In this study, linear loads are injected by the harmonic current. This creates nonlinear loads in the two considered systems and it makes the current drawn unlike the supply voltage waveform. The 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](8)THDV,i%=∑h≠1HVih2Vi1⋅100.Individual harmonic distortion at each bus can be defined by [44](9)IHDV,ih%=VihVi1⋅100.According to IEEE Std. 519 [44], THDV,i and IHDV,ih should not be greater than 5% and 3%, respectively.The average ofTHDV,i can be defined by [25](10)THDVaver=∑i=1NbusTHDV,iNbus.Similarly, the average ofIHDV,ih can be defined by(11)IHDVaverh=∑i=1NbusIHDV,ihNbus.Here, the objective function should be shown as follows:(12)minimizeF4=ω1⋅THDVaver+ω2⋅maxIHDVaverh,where ω1 and ω2 are, respectively, weighted factors with THDVaver and IHDVaverh.ω1 and ω2 are constrained by(13)ω1+ω2=1,0≤ω1,ω2≤1.


## 2.1. Assumption

The 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. The 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.(3)
Placement of solar photovoltaic distribution generation units: suppose that all buses (except for slack bus) in the considered system can be selected for installation.(4)
Lines: transmission capacity limits are predetermined for each branch and remain unchanged during working time. Parameters such as impedance represented the characteristics of the lines and are not affected by over temperature.(5)
Loads: this includes nonlinear loads and linear loads. 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. This 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.


## 2.2. 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 harmonic distortion (THD & IHD). The objective function along with the constraints is explained by using mathematical formulas as follows.


## 2.3. 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. Therefore, the reduction of total active power losses is becoming an important objective, which is shown in the following mathematical equation [41]:(1)TPLDG=∑n=1NbrInDG2Rn.By using the same method, total power losses in all branches before installing PVDGUs, TPL, is determined by means of the following equation [41]:(2)TPL=∑n=1NbrIn2Rn.It shows that if PVDGUs are properly installed, the losses of power system after installation are always smaller than before installation. Therefore, 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:(3)minimizeF1=TPLDGTPL.


## 2.4. The PVDGUs’ Penetration Level

Integration of large distributed generation into electric system can increase investment and operating costs. Therefore, minimizing the capacity of PVDGUs (RPDG) that still satisfies the technical criteria should be considered as the target. The proposed mathematical equation of the objective function is shown as follows [42]:(4)minimizeF2=∑i=1NDGRPDG,i∑i=1NloadPLoadi.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.


## 2.5. Voltage Profile Index

One of the important criteria for evaluating the stability of the operating grid is the voltage profile index (VPI). Thus, the voltage profile index must be considered and the calculation process for VPI should follow the step below [43].The voltage profile for theith bus can be defined as(5)VPi=Vi−Vmin⋅Vmax−ViVnom−Vmin⋅Vmax−Vnom,where Vnom is the nominal or desired bus voltage value and typically taken as 1.0 pu. It is clear that VPi can get a maximum value (1.0 pu) if the bus voltage equals to Vnom and VPi can get a negative value if the bus voltage is higher than the maximum limit or lower than the minimum limit. The overall voltage profile index for the distribution system can be defined as(6)maximum VPI=1Nbus⋅∑i=1NbusVPi.The 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:(7)minimizeF3=1−VPI.


## 2.6. Harmonic Distortions

In this study, linear loads are injected by the harmonic current. This creates nonlinear loads in the two considered systems and it makes the current drawn unlike the supply voltage waveform. The 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](8)THDV,i%=∑h≠1HVih2Vi1⋅100.Individual harmonic distortion at each bus can be defined by [44](9)IHDV,ih%=VihVi1⋅100.According to IEEE Std. 519 [44], THDV,i and IHDV,ih should not be greater than 5% and 3%, respectively.The average ofTHDV,i can be defined by [25](10)THDVaver=∑i=1NbusTHDV,iNbus.Similarly, the average ofIHDV,ih can be defined by(11)IHDVaverh=∑i=1NbusIHDV,ihNbus.Here, the objective function should be shown as follows:(12)minimizeF4=ω1⋅THDVaver+ω2⋅maxIHDVaverh,where ω1 and ω2 are, respectively, weighted factors with THDVaver and IHDVaverh.ω1 and ω2 are constrained by(13)ω1+ω2=1,0≤ω1,ω2≤1.


## 3. Constraints

### 3.1. The 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]. Thus, power balance constraint is considered at the fundamental frequency and has the following form:(14)∑i=1NloadPLoad,i+∑n=1NbrPLoss,n−∑i=1NDGRPDG,i−PGrid=0.


### 3.2. The Voltage Limits

The rms bus voltage magnitude and voltage limits are calculated and should be maintained in the voltage limits [44, 45] as follows:(15)Vmin≤Vi12+∑h=2HVih2≤Vmax,i=1,…,Nbus.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]. This range is selected to become the constraint in the 33-bus distribution system and 69-bus distribution system.


### 3.3. 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. The limits of total voltage harmonic distortion and individual voltage harmonic distortion should follow IEEE Std. 519 [44]:(16)THDV,i≤THDmax,IHDV,ih≤IHDmax,where THDmax and IHDmax are 5% and 3%, respectively.


### 3.4. PVDGUs’ Capacity Limit

The 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. Thus, 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]. The PVDGUs’ capacity is constrained by the following inequalities:(17)RPDGmin≤RPDG,i≤RPDGmax,∑i=1NDGRPDG,i≤∑i=1NloadRPLoad,i,where RPDGmin and RPDGmax are the minimum and maximum rated powers of PVDGUs, respectively.


### 3.5. The Branch Current Limits

The branch current must be kept in the allowable limit of conductor [47] as indicated in the model below:(18)InDG≤Inmax,where InDG and Inmax are the current magnitude of nth branch with PVDGUs and the maximum current magnitude of each branch, respectively.


## 3.1. The 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]. Thus, power balance constraint is considered at the fundamental frequency and has the following form:(14)∑i=1NloadPLoad,i+∑n=1NbrPLoss,n−∑i=1NDGRPDG,i−PGrid=0.


## 3.2. The Voltage Limits

The rms bus voltage magnitude and voltage limits are calculated and should be maintained in the voltage limits [44, 45] as follows:(15)Vmin≤Vi12+∑h=2HVih2≤Vmax,i=1,…,Nbus.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]. This range is selected to become the constraint in the 33-bus distribution system and 69-bus distribution system.


## 3.3. 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. The limits of total voltage harmonic distortion and individual voltage harmonic distortion should follow IEEE Std. 519 [44]:(16)THDV,i≤THDmax,IHDV,ih≤IHDmax,where THDmax and IHDmax are 5% and 3%, respectively.


## 3.4. PVDGUs’ Capacity Limit

The 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. Thus, 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]. The PVDGUs’ capacity is constrained by the following inequalities:(17)RPDGmin≤RPDG,i≤RPDGmax,∑i=1NDGRPDG,i≤∑i=1NloadRPLoad,i,where RPDGmin and RPDGmax are the minimum and maximum rated powers of PVDGUs, respectively.


## 3.5. The Branch Current Limits

The branch current must be kept in the allowable limit of conductor [47] as indicated in the model below:(18)InDG≤Inmax,where InDG and Inmax are the current magnitude of nth branch with PVDGUs and the maximum current magnitude of each branch, respectively.


## 4. Improved Coyote Optimization Algorithm (ICOA)

### 4.1. 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 intoNg groups where each group has NCo coyotes. Thus, population size of COA is Npop (where Npop = NCo × Ng).(3)
Thekth coyote in the gth group has two important factors consisting of social condition Cok,g and quality FFk,g in which the former is a solution and the latter is the fitness function of the solution.The COA method consists of the following steps:Step 1.
Create initial coyote community
The 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:(19)Cok,g=Comin+γ⋅Comax−Comin;k=1,…,NCo,g=1,…,Ng,where Comin and Comax are, respectively, lower bound solution and upper bound solution obtained by(20)Comin=varxmin;x=1,…,Ndv,Comax=varxmax;x=1,…,Ndv,where varxmin and varxmax are the minimum value and maximum value of the xth control variable, respectively.Step 2.
Determine fitness function value
Social conditions of coyote community are evaluated. This action is similar to the task of determining the fitness function of solutions. The effectiveness level of each new solution is ranked based on a fitness function, which is 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:(21)FFk,g=OFk,g+PF⋅∑m=1MPTk,g,m2,where OFk,g is the objective function of the kth solution in the gth group; PF is a penalty factor, which is selected by experience; PTk,g,m is the penalty term for violating the mth constraint of the kth solution in the gth 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 byCobest,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. The detail of the newly updated solution strategy can be expressed by the following model:(22)Cok,gnew=Cok,g+γCobest,g−Co1.g+γComid,g−Co2,g.
In equation (25), Co1,g and Co2,g are two randomly picked solutions in the gth considered group while Comid,g is the middle solution obtained from the gth group. Each decision variable in the middle solution is obtained by one out of two ways shown in the following equation:(23)Comid,g=varx,mid1,ifNCois odd,varx,mid2,else,x=1,…Ndv,where varx,mid1 and varx,mid2 are the xth center decision variables corresponding to two different cases of NCo. If NCo is an odd number, varx,mid1 is taken from the xth decision variable at the position [(NCo + 1)/2]. On the contrary, varx,mid2 is the average of two xth decision variables at two positions, NCo/2 and [(NCo + 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
The new solutions are evaluated by calculating fitness function (21).Step 5.
Perform selection technique
Now, each coyotek has two social conditions, old social condition corresponding to old solution, Cok,g, and new social condition corresponding to new solution, Cok,gnew. The quality of the two social conditions is the fitness function of the two solutions, FFk,g and FFk,gnew. So, it should retain one social condition for each kth coyote by using the following rules:(24)Cok,g=Cok,gnew,ifFFk,gnew≤FFk,g,Cok,g,else,FFk,g=FFk,gnew,ifFFk,gnew≤FFk,g,FFk,g,else.Step 6.
Create one new social condition in each group
In the step, each groupgth produces one new solution Cognew by using a randomization mechanism. The new solution and the technique are shown in the following models:(25)Cognew=varx,gnew;x=1,…,Ndv,(26)varx,gnew=varx,1,gnew,ifβ<1Ndv,varx,2,gnew,if1Ndv≤β<1Ndv+0.5,varx,gnew,else,where varx,1,gnew and varx,2,gnew are two xth variables randomly picked from two available solutions; varx,gnew is the randomly produced variable within lower and upper boundaries; and β is a random number in the range of [0, 1]. Then, the new solution is evaluated by using (21) and its quality is symbolized by FFgnew.Step 7.
Determine the worst solution and replace it
Each new solution in each group obtained by Step6 will be compared to the worst solution in each group Cogworst for deciding if the worst solution should be replaced with the new one. The decision is taken over by the comparison below:(27)Cogworst=Cogworst,ifFFgworst<FFgnew,Cognew,else.Step 8.
Exchange solutions amongNg groups
To diversify solutions and avoid falling into local optimal zones with premature convergence, COA carries out solution exchange action amongNg groups. Two randomly selected groups will provide two randomly picked solutions and their position is then exchanged. The task is dependent on the randomization mechanism below:(28)φ<10−22NCo2,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. The higher NCo is, the higher the probability of exchanged solutions is.Step 9.
Determine the best solution CoGbest
Before terminating one computation iteration, the best solution amongNg groups is determined by comparing fitness function. If the current iteration is the last one, the best solution CoGbest is one candidate solution of the current run.


### 4.2. 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 (NCo × Ng + 1 × Ng) equaling to (Npop + Ng). Thus, 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. The 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.


### 4.3. The First Improvement in the First Generation

In the first improvement, the center solutionComid,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. The 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. The 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:(29)Cok,gnew=Cok,g+γ⋅Cobest,g−Co1,g+γ⋅CoGbest−Co2,g.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 hasNdv decision variables and the xth decision variable of all solutions must be grouped and then sorted. Finally, equation (23) is applied for determining the most appropriate xth decision variable. The task becomes time-consuming if Ndv 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.The performance of the first modification will be investigated by comparing the original COA and COA with the first modification, called ICOA1.


### 4.4. The 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. The number of decision variables Ndv 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:(30)CEg=CoGbest+γ∗CoGbest−Cobest,1+γ∗CoGbest−Cobest,2,(31)CEg=CoGbest+γ∗CoGbest−Cobest,1+γ∗CoGbest−Cobest,2+γ∗CoGbest−Cobest,3.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 CoGbest and Cobest,1 and the second step size between CoGbest and Cobest,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. The 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 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). The 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 Countmax 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 Pa for each specific optimization problem. As Pa 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 Pa on the obtained result and the most appropriate value of Pa will be recommended.(32)Countmax=NCo×NCo+12,(33)R=CountCountmax.Algorithm 1: The proposed technique for the second new solution generation.
Count = 0;
Calculate Countmax by using equation (32)
fori = 1 to (NCo − 1)
forj = i + 1 to NCo
if
FF
j
,
g
−
FF
i
,
g
<
FF
mean
,
g
−
FF
best
,
g
Count = Count + 1
end
end
end
CalculateR by using equation (33)
IfR > Pa
Generate new solutions by using equation (31)
else
Generate new solutions by using equation (30)
end


## 4.1. 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 intoNg groups where each group has NCo coyotes. Thus, population size of COA is Npop (where Npop = NCo × Ng).(3)
Thekth coyote in the gth group has two important factors consisting of social condition Cok,g and quality FFk,g in which the former is a solution and the latter is the fitness function of the solution.The COA method consists of the following steps:Step 1.
Create initial coyote community
The 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:(19)Cok,g=Comin+γ⋅Comax−Comin;k=1,…,NCo,g=1,…,Ng,where Comin and Comax are, respectively, lower bound solution and upper bound solution obtained by(20)Comin=varxmin;x=1,…,Ndv,Comax=varxmax;x=1,…,Ndv,where varxmin and varxmax are the minimum value and maximum value of the xth control variable, respectively.Step 2.
Determine fitness function value
Social conditions of coyote community are evaluated. This action is similar to the task of determining the fitness function of solutions. The effectiveness level of each new solution is ranked based on a fitness function, which is 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:(21)FFk,g=OFk,g+PF⋅∑m=1MPTk,g,m2,where OFk,g is the objective function of the kth solution in the gth group; PF is a penalty factor, which is selected by experience; PTk,g,m is the penalty term for violating the mth constraint of the kth solution in the gth 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 byCobest,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. The detail of the newly updated solution strategy can be expressed by the following model:(22)Cok,gnew=Cok,g+γCobest,g−Co1.g+γComid,g−Co2,g.
In equation (25), Co1,g and Co2,g are two randomly picked solutions in the gth considered group while Comid,g is the middle solution obtained from the gth group. Each decision variable in the middle solution is obtained by one out of two ways shown in the following equation:(23)Comid,g=varx,mid1,ifNCois odd,varx,mid2,else,x=1,…Ndv,where varx,mid1 and varx,mid2 are the xth center decision variables corresponding to two different cases of NCo. If NCo is an odd number, varx,mid1 is taken from the xth decision variable at the position [(NCo + 1)/2]. On the contrary, varx,mid2 is the average of two xth decision variables at two positions, NCo/2 and [(NCo + 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
The new solutions are evaluated by calculating fitness function (21).Step 5.
Perform selection technique
Now, each coyotek has two social conditions, old social condition corresponding to old solution, Cok,g, and new social condition corresponding to new solution, Cok,gnew. The quality of the two social conditions is the fitness function of the two solutions, FFk,g and FFk,gnew. So, it should retain one social condition for each kth coyote by using the following rules:(24)Cok,g=Cok,gnew,ifFFk,gnew≤FFk,g,Cok,g,else,FFk,g=FFk,gnew,ifFFk,gnew≤FFk,g,FFk,g,else.Step 6.
Create one new social condition in each group
In the step, each groupgth produces one new solution Cognew by using a randomization mechanism. The new solution and the technique are shown in the following models:(25)Cognew=varx,gnew;x=1,…,Ndv,(26)varx,gnew=varx,1,gnew,ifβ<1Ndv,varx,2,gnew,if1Ndv≤β<1Ndv+0.5,varx,gnew,else,where varx,1,gnew and varx,2,gnew are two xth variables randomly picked from two available solutions; varx,gnew is the randomly produced variable within lower and upper boundaries; and β is a random number in the range of [0, 1]. Then, the new solution is evaluated by using (21) and its quality is symbolized by FFgnew.Step 7.
Determine the worst solution and replace it
Each new solution in each group obtained by Step6 will be compared to the worst solution in each group Cogworst for deciding if the worst solution should be replaced with the new one. The decision is taken over by the comparison below:(27)Cogworst=Cogworst,ifFFgworst<FFgnew,Cognew,else.Step 8.
Exchange solutions amongNg groups
To diversify solutions and avoid falling into local optimal zones with premature convergence, COA carries out solution exchange action amongNg groups. Two randomly selected groups will provide two randomly picked solutions and their position is then exchanged. The task is dependent on the randomization mechanism below:(28)φ<10−22NCo2,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. The higher NCo is, the higher the probability of exchanged solutions is.Step 9.
Determine the best solution CoGbest
Before terminating one computation iteration, the best solution amongNg groups is determined by comparing fitness function. If the current iteration is the last one, the best solution CoGbest is one candidate solution of the current run.


## 4.2. 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 (NCo × Ng + 1 × Ng) equaling to (Npop + Ng). Thus, 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. The 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.


## 4.3. The First Improvement in the First Generation

In the first improvement, the center solutionComid,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. The 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. The 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:(29)Cok,gnew=Cok,g+γ⋅Cobest,g−Co1,g+γ⋅CoGbest−Co2,g.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 hasNdv decision variables and the xth decision variable of all solutions must be grouped and then sorted. Finally, equation (23) is applied for determining the most appropriate xth decision variable. The task becomes time-consuming if Ndv 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.The performance of the first modification will be investigated by comparing the original COA and COA with the first modification, called ICOA1.


## 4.4. The 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. The number of decision variables Ndv 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:(30)CEg=CoGbest+γ∗CoGbest−Cobest,1+γ∗CoGbest−Cobest,2,(31)CEg=CoGbest+γ∗CoGbest−Cobest,1+γ∗CoGbest−Cobest,2+γ∗CoGbest−Cobest,3.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 CoGbest and Cobest,1 and the second step size between CoGbest and Cobest,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. The 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 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). The 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 Countmax 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 Pa for each specific optimization problem. As Pa 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 Pa on the obtained result and the most appropriate value of Pa will be recommended.(32)Countmax=NCo×NCo+12,(33)R=CountCountmax.Algorithm 1: The proposed technique for the second new solution generation.
Count = 0;
Calculate Countmax by using equation (32)
fori = 1 to (NCo − 1)
forj = i + 1 to NCo
if
FF
j
,
g
−
FF
i
,
g
<
FF
mean
,
g
−
FF
best
,
g
Count = Count + 1
end
end
end
CalculateR by using equation (33)
IfR > Pa
Generate new solutions by using equation (31)
else
Generate new solutions by using equation (30)
end


## 5. The Implementation of the Proposed Method for Determining Site and Rated Power of PVDGUs

### 5.1. 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 Table1. Harmonic spectrum contains parameters of 5 harmonic orders with magnitudes (%) and phase angles (°). These five harmonic flows are injected into the linear loads and treat these loads as harmonic generation sources. The 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 69-bus 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]. The 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). The 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 jth solar photovoltaic distribution generation unit (PVDGU) corresponding to the kth solution in the gth group are, respectively, represented by Sitej,k,g and RPj,k,g. Each solution in each group is represented by the following equation:(34)Cog,k=Sitej,g,k,RPj,g,k;j=1,…,NDG,where NDG is the number of installed PVDGUs in the distribution network.Table 1
Five harmonic flows injected into loads in distribution power networks.
Harmonic order
Magnitude (%)
Angle (°)
5; 7; 11; 13; 17
0.765; 0.627; 0.248; 0.127; 0.071
28; −180; −59; 79; −253In order to produce the initial population, lower bound solution Comin and upper bound solution Comax should be defined to make sure that randomly generated solutions can satisfy the limit constraints:(35)Comin=Sitejmin,RPjmin;j=1,…,NDG,Comax=Sitejmax,RPjmax;j=1,…,NDG,where Sitejmin and Sitejmax are, respectively, bus 2 and bus Nbus and RPjmin and RPjmax are the minimum rated power and maximum rated power of the jth PVDG unit. However, the minimum power of all units and the maximum power of all units are the same and corresponding to RPmin and RPmax. The limits of site are easily selected based on the system dimension, but the limits of rated power must be selected thoroughly and exactly. The limits will be discussed in the numerical results. Each solution in the initial population is formed by the following way:(36)Cok,g=Comin+ε⋅Comax−Comin;k=1,…,NCo;g=1,…,Ng.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. The task continues to be done after producing new solutions for the first and the second generations.


### 5.2. Penalty Terms for Violating of Dependent Variables

After having decision variables, FW/BWST technique has been run for obtaining dependent variables including branch currentsInDGh and load bus voltage Vih 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. Then, current in branch n, InDG, and voltage at bus i, Vi, are calculated by using the two following equations:(37)InDG=InDG12+∑h=2HInDGh2,Vi=Vi12+∑h=2HVih2.As a result, the two parameters are checked and penalized by using the following formulas:(38)ΔVi=Vi−Vimax,ifVi>Vmax,Vimin−Vi,ifVi<Vmin,0,else,ΔIn=InDG−Inmax,ifInDG>Inmax,0,else.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). Then, they are checked and penalized as follows:(39)ΔTHDi=THDi−THDmax,ifTHDi>THDmax,0,else,ΔIHDi=IHDi−IHDmax,ifIHDi>IHDmax,0,else.


### 5.3. Fitness Function

As described in Section2, we focus on four single objectives consisting of total power loss minimization F1, rated power minimization of PVDGUs F2, load bus voltage enhancement F3, and harmonic distortion minimization F4. According to the considered objectives, fitness function can be established as follows:(i)
Fitness function for total power loss minimization:(40)FF1=F1+μV∑i=1NbusΔVi2+μI∑i=1NbrΔIn2+μTHD∑i=1NbusΔTHDi2+μIHD∑i=1NbusΔIHDi2.(ii)
Fitness function for rated power minimization of PVDGUs:(41)FF2=F2+μV∑i=1NbusΔVi2+μI∑i=1NbrΔIn2+μTHD∑i=1NbusΔTHDi2+μIHD∑i=1NbusΔIHDi2.(iii)
Fitness function for improving load bus voltage:(42)FF3=F3+μV∑i=1NbusΔVi2+μI∑i=1NbrΔIn2+μTHD∑i=1NbusΔTHDi2+μIHD∑i=1NbusΔIHDi2.(iv)
Fitness function for harmonic distortion minimization:(43)FF4=F4+μV∑i=1NbusΔVi2+μI∑i=1NbrΔIn2.


### 5.4. Producing New Solutions and Fixing Violated Variables

The proposed ICOA experiences two new solution generations. The first generation employs equation (29) to produce NCo 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:(44)Sitej=2,ifSitej<2,Nbus,ifSitej>Nbus,Sitej,else,RPj=RPmin,ifRPj<RPmin,RPmax,ifRPj>RPmax,RPj,else.Between site and rated power, site is approximated for getting integer value by using command “round” while rated power remains unchanged.


### 5.5. 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 ITmax iterations in which ITmax is the maximum number of iterations. So, the best run and the fifty runs are recorded and plotted in figures.


### 5.6. The 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 Figure1.Figure 1
Flowchart of the computation process.


## 5.1. 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 Table1. Harmonic spectrum contains parameters of 5 harmonic orders with magnitudes (%) and phase angles (°). These five harmonic flows are injected into the linear loads and treat these loads as harmonic generation sources. The 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 69-bus 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]. The 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). The 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 jth solar photovoltaic distribution generation unit (PVDGU) corresponding to the kth solution in the gth group are, respectively, represented by Sitej,k,g and RPj,k,g. Each solution in each group is represented by the following equation:(34)Cog,k=Sitej,g,k,RPj,g,k;j=1,…,NDG,where NDG is the number of installed PVDGUs in the distribution network.Table 1
Five harmonic flows injected into loads in distribution power networks.
Harmonic order
Magnitude (%)
Angle (°)
5; 7; 11; 13; 17
0.765; 0.627; 0.248; 0.127; 0.071
28; −180; −59; 79; −253In order to produce the initial population, lower bound solution Comin and upper bound solution Comax should be defined to make sure that randomly generated solutions can satisfy the limit constraints:(35)Comin=Sitejmin,RPjmin;j=1,…,NDG,Comax=Sitejmax,RPjmax;j=1,…,NDG,where Sitejmin and Sitejmax are, respectively, bus 2 and bus Nbus and RPjmin and RPjmax are the minimum rated power and maximum rated power of the jth PVDG unit. However, the minimum power of all units and the maximum power of all units are the same and corresponding to RPmin and RPmax. The limits of site are easily selected based on the system dimension, but the limits of rated power must be selected thoroughly and exactly. The limits will be discussed in the numerical results. Each solution in the initial population is formed by the following way:(36)Cok,g=Comin+ε⋅Comax−Comin;k=1,…,NCo;g=1,…,Ng.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. The task continues to be done after producing new solutions for the first and the second generations.


## 5.2. Penalty Terms for Violating of Dependent Variables

After having decision variables, FW/BWST technique has been run for obtaining dependent variables including branch currentsInDGh and load bus voltage Vih 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. Then, current in branch n, InDG, and voltage at bus i, Vi, are calculated by using the two following equations:(37)InDG=InDG12+∑h=2HInDGh2,Vi=Vi12+∑h=2HVih2.As a result, the two parameters are checked and penalized by using the following formulas:(38)ΔVi=Vi−Vimax,ifVi>Vmax,Vimin−Vi,ifVi<Vmin,0,else,ΔIn=InDG−Inmax,ifInDG>Inmax,0,else.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). Then, they are checked and penalized as follows:(39)ΔTHDi=THDi−THDmax,ifTHDi>THDmax,0,else,ΔIHDi=IHDi−IHDmax,ifIHDi>IHDmax,0,else.


## 5.3. Fitness Function

As described in Section2, we focus on four single objectives consisting of total power loss minimization F1, rated power minimization of PVDGUs F2, load bus voltage enhancement F3, and harmonic distortion minimization F4. According to the considered objectives, fitness function can be established as follows:(i)
Fitness function for total power loss minimization:(40)FF1=F1+μV∑i=1NbusΔVi2+μI∑i=1NbrΔIn2+μTHD∑i=1NbusΔTHDi2+μIHD∑i=1NbusΔIHDi2.(ii)
Fitness function for rated power minimization of PVDGUs:(41)FF2=F2+μV∑i=1NbusΔVi2+μI∑i=1NbrΔIn2+μTHD∑i=1NbusΔTHDi2+μIHD∑i=1NbusΔIHDi2.(iii)
Fitness function for improving load bus voltage:(42)FF3=F3+μV∑i=1NbusΔVi2+μI∑i=1NbrΔIn2+μTHD∑i=1NbusΔTHDi2+μIHD∑i=1NbusΔIHDi2.(iv)
Fitness function for harmonic distortion minimization:(43)FF4=F4+μV∑i=1NbusΔVi2+μI∑i=1NbrΔIn2.


## 5.4. Producing New Solutions and Fixing Violated Variables

The proposed ICOA experiences two new solution generations. The first generation employs equation (29) to produce NCo 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:(44)Sitej=2,ifSitej<2,Nbus,ifSitej>Nbus,Sitej,else,RPj=RPmin,ifRPj<RPmin,RPmax,ifRPj>RPmax,RPj,else.Between site and rated power, site is approximated for getting integer value by using command “round” while rated power remains unchanged.


## 5.5. 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 ITmax iterations in which ITmax is the maximum number of iterations. So, the best run and the fifty runs are recorded and plotted in figures.


## 5.6. The 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 Figure1.Figure 1
Flowchart of the computation process.


## 6. 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. The 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. The 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. The 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. The detailed information of the five harmonics is shown in Table 1 [50].Figure 2
IEEE 33-bus distribution network.Figure 3
IEEE 69-bus distribution network.The 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 69-bus distribution network. The optimal location and optimal capacity of PVDGUs in the two systems aim at reaching four single-objective functions considered in the four following cases.Case 1: Minimizing the total power lossCase 2: Minimizing the rated power of all PVDGUsCase 3: Improving the voltage profile of all loadsCase 4: Minimizing the harmonic distortionIn 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.The 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 (mmax) is set to five values from 0.1 to 0.5; habitat modification probability (Pmod) 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 values, 0.6, 0.8, and 1.0, and mutation probability (Mu) is from 0.1 to 0.5 with a step of 0.1 [10].(3)
For running PSO, two acceleration factorsc1 and c2 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,c1 is a function of 2e−4IT/ITmax2; c2 and c3 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 groupNCo and the number of groups are set to 4. According to the settings, the population Npop of the four methods is 20. In addition, ICOA2 and ICOA have more than one advanced control parameter Pa and the parameter is set to five values from 0.2 to 1.0 with a step of 0.2. The impact of Pa on simulation results will be surveyed and discussed.On the other hand, the applied methods have two common control parameters: they are populationNpop and the maximum number of iteration ITmax. Npop of PSO, GA, BBO, SSA, and SFO is also set to that of COA methods, i.e., 20, while ITmax has the same values for all methods. ITmax is 75 for the IEEE 33-bus distribution network and 100 for the IEEE 69-bus distribution network.


### 6.1. Survey on Results Obtained by Setting Different Values ofPa for ICOA Method

For finding the most suitablePa value, ICOA is in turn run 50 times corresponding to each Pa value by setting Pa to 0.2, 0.4, 0.6, 0.8, and 1.0. The second case of minimizing the rated power of PVDGUs is selected for both IEEE 33-bus distribution network and IEEE 69-bus distribution network. The obtained results for the two systems are shown in Tables 2 and 3, respectively. As observed from Table 2, the best fitness corresponding to Pa = 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. Through the results, it points out that ICOA method can reach the best solution only at Pa = 0.2. Additionally, the average values of fitness function also confirm the most effective impact of Pa = 0.2 on its average fitness is the smallest among five values. The best fitness can conclude Pa = 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 Pa = 0.2. One more time, the significance of setting Pa = 0.2 can be adopted since the best fitness and the average fitness at Pa = 0.2 continue to be the smallest values among five values for the second system shown in Table 3. The best fitness and the best average are, respectively, 0.2425 and 0.2939. Consequently, we accept the best Pa value is 0.2 and it is used for other remaining cases with intent to shorten simulation time.Table 2
The obtained results by ICOA for Case2 of the IEEE 33-bus distribution network.
P
a
0.0
0.2
0.4
0.6
0.8
1.0
Best FF2
0.2860
0.2836
0.2978
0.2981
0.2900
0.2978
Average FF2
0.3881
0.3410
0.3474
0.3485
0.3451
0.3450
Worst FF2
0.5831
0.4109
0.4282
0.4326
0.4429
0.4448Table 3
The obtained results by ICOA for Case2 of the IEEE 69-bus distribution network.
P
a
0.0
0.2
0.4
0.6
0.8
1.0
Best FF2
0.2426
0.2425
0.2425
0.2426
0.2426
0.2426
Average FF2
0.3285
0.2939
0.2969
0.3130
0.3165
0.3271
Worst FF2
0.5143
0.4290
0.4716
0.5002
0.5104
0.5144


### 6.2. 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)
The best solution of the proposed method and compared methods was compared to show the strong search ability of the proposed ICOA.(2)
The 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)
The number of better runs (Nbetterrun) 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. The three factors are obtained by using the following equations:(45)IL of the best fitness%=the best fitness of the compared method−the best fitness of ICOAthe best fitness of the compared method×100,(46)IL of average fitness%=average fitness of the compared method−average fitness of ICOAaverage fitness of the compared method×100,(47)IL of the worst fitness%=the worst fitness of the compared method−the worst fitness of ICOAthe worst fitness of the compared method×100.(5)
Fair control parameter setting for the proposed ICOA and other compared methods was considered during the process of selecting population size and the maximum number of iterations.The implemented methods in this study are based on metaheuristic algorithms. Therefore, 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.


### 6.3. 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 FF1 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. The best values indicate that ICOA1 and ICOA can find the same best fitness less than that of COA and ICOA2. The average value of ICOA1 is better than that of COA but that of ICOA2 is worse than that of COA. Thus, 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. Figure4 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. The 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. There 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%. Figures4 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. Thus, 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.Table 4
The comparison among nine methods of the IEEE 33-bus distribution network for the first case.
Method
Best FF1
Average FF1
Worst FF1
IL of the best FF1
IL of average FF1
IL of the worst FF1
N
betterruns
BBO
0.3465
0.3782
0.4128
0.115
3.860
4.239
15
GA
0.3890
0.5027
0.6106
11.028
27.671
35.260
45
PSO
0.3697
0.4184
0.4492
6.384
13.098
11.999
27
SFO
0.3793
0.4184
0.446
8.753
13.098
11.368
38
SSA
0.3617
0.4054
0.4589
4.313
10.311
13.859
27
COA
0.3462
0.3727
0.4086
0.029
2.442
3.255
13
ICOA1
0.3461
0.3674
0.3963
0.000
1.034
0.252
8
ICOA2
0.3462
0.3779
0.4148
0.029
3.784
4.701
13
ICOA
0.3461
0.3636
0.3953
—
—
—
—Figure 4
The fitness values of four methods in 50 trial times for the first case.Figure 5
The fitness values of six methods in 50 trial times for the first case.Table 5
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case1.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 24; size: 1.0289 MWBus: 14; size: 0.8435 MWBus: 30; size: 1.0669 MW
4.612/2.994
0.969/1.000
73.108
GA
Bus: 13; size: 0.7711 MWBus: 06; size: 1.6119 MWBus: 31; size: 0.4229 MW
4.404/2.859
0.968/1.000
82.039
PSO
Bus: 14; size: 1.0746 MWBus: 25; size: 0.9004 MWBus: 30; size: 0.9302 MW
4.241/2.752
0.968/1.000
77.999
SFO
Bus: 24; size: 0.9188 MWBus: 09; size: 1.4072 MWBus: 30; size: 0.7580 MW
4.617/2.996
0.967/1.000
80.018
SSA
Bus: 14; size: 0.8744 MWBus: 30; size: 0.8465 MWBus: 24; size: 0.7361 MW
4.595/2.982
0.962/1.000
76.309
COA
Bus: 30; size: 1.0341 MWBus: 14; size: 0.8397 MWBus: 24; size: 1.1157 MW
4.620/2.998
0.969/1.000
73.038
ICOA1
Bus: 30; size: 1.0453 MWBus: 14; size: 0.8400 MWBus: 24; size: 1.0903 MW
4.619/2.998
0.969/1.000
73.031
ICOA2
Bus: 30; size: 1.0345 MWBus: 14; size: 0.8389 MWBus: 24; size: 1.1275 MW
4.621/2.999
0.969/1.000
73.043
ICOA
Bus: 24; size: 1.0885 MWBus: 30; size: 1.0467 MWBus: 14; size: 0.8381 MW
4.622/3.000
0.969/1.000
73.017Case 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 FF2 in equation (41) is used to evaluate the effectiveness of solutions from implemented approaches.
The 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 Table6. 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. The 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 Case1 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. The 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. The 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 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.Table 6
The comparison among nine methods of the IEEE 33-bus distribution network for the second case.
Method
Best FF2
Average FF2
Worst FF2
IL of the best FF2
IL of average FF2
IL of the worst FF2
N
betterruns
BBO
0.3024
0.3724
0.4553
6.217
8.432
9.752
5
GA
0.3582
0.5401
0.6977
20.826
36.864
41.106
41
PSO
0.3321
0.4027
0.5327
14.604
15.322
22.865
16
SFO
0.3702
0.4933
0.6259
23.393
30.874
34.351
43
SSA
0.3467
0.4474
0.5879
18.200
23.782
30.107
39
COA
0.3008
0.3633
0.4448
5.718
6.138
7.621
3
ICOA1
0.2975
0.3444
0.4228
4.672
0.987
2.815
2
ICOA2
0.3008
0.3746
0.4466
5.718
8.970
7.994
3
ICOA
0.2836
0.3410
0.4109
—
—
—
—Figure 6
The fitness values of four methods in 50 trial times for the second case.Figure 7
The fitness values of six methods in 50 trial times for the second case.Table 7
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case2.
Method
The best solution
THD/IHD (max) with PVDGUs
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 32; size: 0.5377 MWBus: 14; size: 0.6255 MWBus: 26; size: 0.0240 MW
4.531/2.941%
0.950/1.000
104.587
GA
Bus: 32; size: 0.2652 MWBus: 30; size: 0.2320 MWBus: 14; size: 0.9090 MW
4.604/2.988%
0.950/1.000
100.328
PSO
Bus: 14; size: 0.7796 MWBus: 31; size: 0.0495 MWBus: 32; size: 0.4747 MW
4.592/2.980%
0.951/1.000
101.133
SFO
Bus: 29; size: 0.2121 MWBus: 15; size: 0.8769 MWBus: 31; size: 0.3643 MW
4.269/2.771 pu
0.952/1.000
98.797
SSA
Bus: 32; size: 0.4525 MWBus: 2; size: 0.0004 MWBus: 14; size: 0.9021 MW
4.620/2.999 pu
0.950/1.000
102.365
COA
Bus: 31; size: 0.1411 MWBus: 14; size: 0.6197 MWBus: 33; size: 0.4200 MW
4.439/2.881 pu
0.950/1.000
104.925
ICOA1
Bus: 33; size: 0.2049 MWBus: 15; size: 0.5772 MWBus: 31; size: 0.3859 MW
4.264/2.767 pu
0.950/1.000
105.181
ICOA2
Bus: 31; size: 0.1505 MWBus: 14; size: 0.6201 MWBus: 33; size: 0.4105 MW
4.393/2.851 pu
0.950/1.000
104.887
ICOA
Bus: 33; size: 0.2781 MWBus: 31; size: 0.3144 MWBus: 15; size: 0.5210 MW
4.617/2.996 pu
0.950/1.000
109.299Case 3.
: Improving the voltage profile
In the part, voltage of all loads is considered to be improved by using fitness function FF3 in equation (42) as a comparison criterion for evaluating applied methods. In the fitness function, objective function F3 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 Table8. The 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. The lowest average fitness function and the lowest worst fitness function confirm that the proposed method is the most stable tool with the lowest fluctuations. The improvement level of the best fitness is from 0.531% to 17.590%. The 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 Table9.Table 8
The comparison among nine methods of the IEEE 33-bus distribution network for the third case.
Method
Best FF3
Average FF3
Worst FF3
IL of the best FF3
IL of average FF3
IL of the worst FF3
N
betterruns
BBO
0.0979
0.1273
0.159
4.290
9.112
8.805
4
GA
0.1137
0.2302
0.3099
17.590
49.739
53.211
28
PSO
0.0989
0.134
0.163
5.258
13.657
11.043
5
SFO
0.0990
0.1886
0.3064
5.354
38.653
52.676
5
SSA
0.1035
0.1634
0.2913
9.469
29.192
50.223
15
COA
0.0951
0.1222
0.1473
1.472
5.319
1.561
3
ICOA1
0.0942
0.1178
0.1451
0.531
1.783
0.069
3
ICOA2
0.0953
0.1228
0.1493
1.679
5.782
2.880
3
ICOA
0.0937
0.1157
0.145
—
—
—
—Figure 8
The fitness values of seven methods in 50 trial times for the third case.Table 9
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case3.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 07; size: 1.0925 MWBus: 30; size: 1.6785 MWBus: 14; size: 0.9046 MW
4.3764/2.840
0.983/1.005
114.7856
GA
Bus: 29; size: 1.9804 MWBus: 15; size: 0.3772 MWBus: 09; size: 1.2137 MW
4.625/3.000
0.983/1.007
124.438
PSO
Bus: 14; size: 0.8685 MWBus: 30; size: 1.4705 MWBus: 07; size: 1.3455 MW
4.437/2.880
0.983/1.002
109.115
SFO
Bus: 30; size: 1.6126 MWBus: 06; size: 1.1535 MWBus: 14; size: 0.9324 MW
4.3308/2.811
0.9835/1.004
113.053
SSA
Bus: 07; size: 0.7343 MWBus: 30; size: 1.7797 MWBus: 14; size: 1.0209 MW
4.208/2.731
0.9829/1.005
116.778
COA
Bus: 14; size: 0.9694 MWBus: 07; size: 1.0156 MWBus: 30; size: 1.6609 MW
4.279/2.777
0.9833/1.004
114.890
ICOA1
Bus: 14; size: 0.9464 MWBus: 07; size: 1.1930 MWBus: 30; size: 1.5831 MW
4.309/2.797
0.9836/1.005
114.680
ICOA2
Bus: 14; size: 0.9662 MWBus: 07; size: 1.0155 MWBus: 30; size: 1.6707 MW
4.283/2.779
0.9833/1.004
115.289
ICOA
Bus: 14; size: 0.9564 MWBus: 07; size: 1.1461 MWBus: 30; size: 1.5730 MW
4.298/2.789
0.9834/1.0043
113.132Case 4.
: Reducing the harmonic distortions
In the part, THD and IHD of all buses are improved by using fitness function FF4 in equation (43) to compare the quality of obtained solutions from implemented methods. In the fitness function, objective function F4 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 Table10 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. The 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%. The 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.
Optimal solutions together with IHD, THD, voltage profile, and total rated power of all PVDGUs obtained by these methods are shown in Table11.Table 10
The comparison among nine methods of the IEEE 33-bus distribution network for the fourth case.
Method
Best FF4
Average FF4
Worst FF4
IL of the best FF4
IL of average FF4
IL of the worst FF4
N
betterruns
BBO
1.1368
1.2524
1.428
0.176
7.833
5.469
26
GA
1.1829
1.4666
1.629
4.066
21.294
17.133
45
PSO
1.1392
1.2516
1.4708
0.386
7.774
8.220
33
SFO
1.1781
1.289
1.5417
3.675
10.450
12.441
45
SSA
1.1403
1.2571
1.4695
0.482
8.178
8.139
36
COA
1.1350
1.202
1.4074
0.018
3.968
4.086
4
ICOA1
1.1348
1.1593
1.3989
0.000
0.431
3.503
1
ICOA2
1.1350
1.2159
1.4619
0.018
5.066
7.661
4
ICOA
1.1348
1.1543
1.3499
—
—
—
—Figure 9
The fitness values of seven methods in 50 trial times for the fourth case.Table 11
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case4.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 19; size: 0.5830 MWBus: 15; size: 1.1293 MWBus: 33; size: 2.0001 MW
2.371/1.538
0.981/1.016
150.419
GA
Bus: 19; size: 0.5807 MWBus: 08; size: 1.6284 MWBus: 11; size: 1.5059 MW
2.699/1.751
0.961/1.022
165.954
PSO
Bus: 15; size: 1.8908 MWBus: 33; size: 1.2287 MWBus: 19; size: 0.5933 MW
2.444/1.585
0.981/1.035
164.398
SFO
Bus: 12; size: 1.7359 MWBus: 19; size: 0.5882 MWBus: 17; size: 1.0495 MW
2.709/1.757
0.955/1.050
218.820
SSA
Bus: 15; size: 1.2310 MWBus: 19; size: 0.5802 MWBus: 32; size: 1.8416 MW
2.404/1.559
0.9814/1.009
136.360
COA
Bus: 17; size: 1.1842 MWBus: 19; size: 0.5739 MWBus: 32; size: 1.9566 MW
2.403/1.558
0.982/1.018
154.611
ICOA1
Bus: 19; size: 0.5741 MWBus: 17; size: 1.1839 MWBus: 32; size: 1.9588 MW
2.402/1.558
0.982/1.018
154.760
ICOA2
Bus: 17; size: 1.1847 MWBus: 19; size: 0.5744 MWBus: 32; size: 1.9548 MW
2.401/1.558
0.982/1.018
154.512
ICOA
Bus: 19; size: 0.5779 MWBus: 17; size: 1.3284 MWBus: 32; size: 1.8086 MW
2.410/1.563
0.982/1.025
156.902


### 6.4. Test System 2: IEEE 69-Bus Distribution Network

In this section, an IEEE 69-bus distribution network is utilized as a main tool for running applied methods and the proposed ICOA method. Four single objectives consisting of total power losses, capacity of all PVDGUs, load voltage profile, and harmonic distortions similar to those in the previous section are in turn independently optimized. For better comparison of the performance of methods, values of fitness functions shown in equations (40)–(43) obtained from 50 runs are summarized in Tables 12–15 and reported in detail in Figures 10–13. Numerical tables can support enough comparison criteria such as quality of the best optimal solution, stability of all runs, and oscillations of runs while figures can enable a clear view of all runs.Table 12
The comparison among nine methods of the IEEE 69-bus distribution network for the first case.
Method
Best FF1
Average FF1
Worst FF1
IL of the best FF1
IL of average FF1
IL of the worst FF1
N
betterruns
BBO
0.3027
0.311
0.3402
0.264
1.383
7.613
8
GA
0.3132
0.3415
0.3968
3.608
10.190
20.791
48
PSO
0.3103
0.3267
0.3468
2.707
6.122
9.371
45
SFO
0.3131
0.3311
0.3632
3.577
7.369
13.464
48
SSA
0.3045
0.3219
0.3891
0.854
4.722
19.224
15
COA
0.3027
0.3085
0.3202
0.264
0.583
1.843
9
ICOA1
0.3024
0.3079
0.3184
0.165
0.390
1.288
3
ICOA2
0.3027
0.3091
0.3221
0.264
0.776
2.422
9
ICOA
0.3019
0.3067
0.3143
—
—
—
—Table 13
The comparison among nine methods of the IEEE 69-bus distribution network for the second case.
Method
Best FF2
Average FF2
Worst FF2
IL of the best FF2
IL of average FF2
IL of the worst FF2
N
betterruns
BBO
0.2427
0.3302
0.5152
0.082
10.993
16.731
1
GA
0.2752
0.4619
0.6518
11.882
36.372
34.182
17
PSO
0.2717
0.4102
0.592
10.747
28.352
27.534
17
SFO
0.3363
0.5192
0.6574
27.892
43.394
34.743
43
SSA
0.2563
0.4344
0.6854
5.384
32.343
37.409
12
COA
0.2427
0.3117
0.4857
0.082
5.711
11.674
3
ICOA1
0.2425
0.2966
0.4527
0.000
0.910
5.235
1
ICOA2
0.2427
0.3165
0.4991
0.082
7.141
14.045
3
ICOA
0.2425
0.2939
0.429
—
—
—
—Table 14
The comparison among nine methods of the IEEE 69-bus distribution network for the third case.
Method
Best FF3
Average FF3
Worst FF3
IL of the best FF3
IL of average FF3
IL of the worst FF3
N
betterruns
BBO
0.0571
0.0838
0.1631
16.462
5.728
20.969
5
GA
0.1122
0.2916
0.3915
57.487
72.908
67.075
44
PSO
0.1020
0.1991
0.3371
53.235
60.321
61.762
40
SFO
0.0829
0.2544
0.4031
42.461
68.947
68.023
34
SSA
0.0658
0.1584
0.3345
27.508
50.126
61.465
15
COA
0.0565
0.0881
0.1457
15.575
10.329
11.531
4
ICOA1
0.0517
0.0795
0.1302
7.737
0.629
0.998
2
ICOA2
0.0568
0.0895
0.1461
16.021
11.732
11.773
4
ICOA
0.0477
0.079
0.1289
—
—
—
—Table 15
The comparison among nine methods of the IEEE 69-bus distribution network for the fourth case.
Method
Best FF4
AverageFF4
Worst FF4
IL of the best FF4
IL of average FF4
IL of the worst FF4
N
betterruns
BBO
0.9547
1.0063
1.1189
0.293
0.646
5.523
4
GA
1.0282
1.0846
1.1964
7.421
7.819
11.643
44
PSO
1.0265
1.0767
1.1364
7.267
7.142
6.978
43
SFO
1.0064
1.0751
1.1269
5.415
7.004
6.194
29
SSA
0.9812
1.0705
1.1496
2.986
6.604
8.046
13
COA
0.9547
1.0029
1.0697
0.293
0.309
1.178
4
ICOA1
0.9526
1.0000
1.0572
0.073
0.020
0.009
2
ICOA2
0.9550
1.0048
1.0808
0.325
0.498
2.193
5
ICOA
0.9519
0.9998
1.0571
—
—
—
—Figure 10
The fitness values of seven methods in 50 trial times for the first case.Figure 11
The fitness values of seven methods in 50 trial times for the second case.Figure 12
The fitness values of seven methods in 50 trial times for the third case.Figure 13
The fitness values of seven methods in 50 trial times for the fourth case.In general, the proposed ICOA method outperforms the original COA method for all cases with high improvement level. The 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 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. The two values are 0.583% and 1.843% for Case1, 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.608% for Case 1, from 0.082% to 27.892% for Case 2, from 16.462% to 57.487% for Case 3, and from 0.293% to 7.421% for Case 4. Furthermore, the improvement levels of the average fitness and the worst fitness are much higher and can be up to 10.190% and 20.791% for Case 1, 43.394% and 37.409% for Case 2, 68.947% and 68.023% for Case 3, and 7.819% and 11.643% for Case 4. The number of runs from ICOA with better fitness than the best fitness of these compared methods is also reported. For almost cases, BBO is the second best method while GA, PSO, and SFO are the worst methods, especially GA. ICOA can find 8, 1, 5, and 4 solutions with better fitness than the best solutions of BBO for Cases 1, 2, 3, and 4, respectively. As compared to GA, PSO, and SFO, ICOA can reach much higher number of runs with better fitness. They are, respectively, 48, 45, and 48 for Case 1, 17, 17, and 43 for Case 2, 44, 40, and 34 for Case 3, and 44, 43, and 29 for Case 4. One more time, the quality of runs from ICOA can be confirmed by observing figures. 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 Tables16–19.Table 16
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case1.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 21; size: 0.3057 MWBus: 49; size: 0.7413 MWBus: 10; size: 0.5958 MWBus: 61; size: 1.7118 MW
4.260/2.753
0.979/1.000
67.951
GA
Bus: 39; size: 0.3089 MWBus: 17; size: 0.5140 MWBus: 61; size: 1.6919 MWBus: 53; size: 0.6702 MW
3.887/2.512
0.980/1.000
70.312
PSO
Bus: 08; size: 0.7504 MWBus: 49; size: 0.7036 MWBus: 61; size: 1.6883 MWBus: 17; size: 0.2579 MW
4.351/2.812
0.978/1.000
69.668
SFO
Bus: 61; size: 1.6826 MWBus: 39; size: 0.3095 MWBus: 17; size: 0.5122 MWBus: 53; size: 0.6785 MW
3.893/2.515
0.980/1.000
70.299
SSA
Bus: 17; size: 0.3835 MWBus: 49; size: 0.9538 MWBus: 61; size: 1.6617 MWBus: 66; size: 0.3690 MW
4.144/2.678
0.976/1.000
68.360
COA
Bus: 61; size: 1.7076 MWBus: 17; size: 0.3894 MWBus: 49; size: 0.9110 MWBus: 66; size: 0.4886 MW
4.115/2.659
0.978/1.000
67.954
ICOA1
Bus: 49; size: 0.7485 MWBus: 10; size: 0.5841 MWBus: 61; size: 1.7367 MWBus: 21; size: 0.3572 MW
4.160/2.688
0.980/1.000
67.889
ICOA2
Bus: 49; size: 0.7493 MWBus: 10; size: 0.5938 MWBus: 21; size: 0.3069 MWBus: 61; size: 1.7068 MW
4.259/2.752
0.979/1.000
67.949
ICOA
Bus: 49; size: 0.7183 MWBus: 21; size: 0.3231 MWBus: 61; size: 1.6881 MWBus: 11; size: 0.6243 MW
4.233/2.735
0.978/1.000
67.765Table 17
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case2.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 24; size: 0.0005 MWBus: 61; size: 0.4310 MWBus: 64; size: 0.5456 MWBus: 03; size: 0.0003 MW
4.453/2.875
0.950/1.000
112.543
GA
Bus: 17; size: 0.0947 MWBus: 64; size: 0.4732 MWBus: 63; size: 0.5328 MWBus: 42; size: 0.0075 MW
3.901/2.519
0.952/1.000
106.372
PSO
Bus: 68; size: 0.0159 MWBus: 69; size: 0.0477 MWBus: 61; size: 0.5273 MWBus: 64; size: 0.5028 MW
4.503/2.902
0.952/1.000
106.844
SFO
Bus: 21; size: 0.3091 MWBus: 61; size: 0.7208 MWBus: 63; size: 0.3180 MWBus: 64; size: 0.0602 MW
3.897/2.519
0.954/1.000
92.108
SSA
Bus: 64; size: 0.1936 MWBus: 24; size: 0.0361 MWBus: 63; size: 0.0633 MWBus: 62; size: 0.7390 MW
4.713/3.000
0.950/1.000
109.137
COA
Bus: 02; size: 0.0003 MWBus: 61; size: 0.4312 MWBus: 64; size: 0.5455 MWBus: 24; size: 0.0004 MW
4.457/2.878
0.950/1.000
112.332
ICOA1
Bus: 15; size: 0.0002 MWBus: 62; size: 0.4975 MWBus: 64; size: 0.4792 MWBus: 02; size: 0.0001 MW
4.557/2.937
0.950/1.000
112.644
ICOA2
Bus: 24; size: 0.0005 MWBus: 61; size: 0.4311 MWBus: 63; size: 0.5456 MWBus: 03; size: 0.0002 MW
4.452/2.875
0.950/1.000
112.543
ICOA
Bus: 15; size: 0.0002 MWBus: 62; size: 0.4976 MWBus: 64; size: 0.4791 MWBus: 02; size: 0.0001 MW
4.557/2.937
0.950/1.000
112.643Table 18
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case3.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 62; size: 2.0022 MWBus: 15; size: 0.5342 MWBus: 55; size: 0.7809 MWBus: 68; size: 0.4592 MW
3.725/2.399
0.994/1.003
85.109
GA
Bus: 11; size: 0.2395 MWBus: 65; size: 1.9581 MWBus: 56; size: 0.9560 MWBus: 17; size: 0.6186 MW
3.710/2.390
0.994/1.020
124.062
PSO
Bus: 65; size: 1.9993 MWBus: 09; size: 0.2583 MWBus: 17; size: 0.5867 MWBus: 56; size: 0.9538 MW
3.823/2.464
0.994/1.022
126.773
SFO
Bus: 55; size: 0.0799 MWBus: 53; size: 1.1452 MWBus: 18; size: 0.5801 MWBus: 61; size: 1.9702 MW
3.569/2.299
0.994/1.001
80.482
SSA
Bus: 60; size: 1.3378 MWBus: 66; size: 0.7844 MWBus: 17; size: 0.4479 MWBus: 63; size: 0.9821 MW
3.860/2.494
0.994/1.003
85.481
COA
Bus: 62; size: 2.0000 MWBus: 55; size: 0.7803 MWBus: 15; size: 0.5341 MWBus: 68; size: 0.4552 MW
3.726/2.400
0.994/1.002
84.943
ICOA1
Bus: 10; size: 0.7891 MWBus: 61; size: 2.0311 MWBus: 21; size: 0.3902 MWBus: 56; size: 0.5777 MW
3.960/2.559
0.994/1.002
82.884
ICOA2
Bus: 62; size: 2.0061 MWBus: 55; size: 0.7815 MWBus: 68; size: 0.4548 MWBus: 15; size: 0.5339 MW
3.724/2.399
0.994/1.003
85.177
ICOA
Bus: 10; size: 0.7786 MWBus: 21; size: 0.4094 MWBus: 56; size: 0.5689 MWBus: 61; size: 2.0001 MW
3.936/2.543
0.994/1.001
81.781Table 19
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case4.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 18; size: 0.2705 MWBus: 23; size: 1.4079 MWBus: 63; size: 1.5097 MWBus: 55; size: 0.6023 MW
2.778/1.795
0.984/1.044
135.120
GA
Bus: 66; size: 0.0256 MWBus: 09; size: 1.0095 MWBus: 60; size: 1.9590 MWBus: 17; size: 0.1785 MW
2.899/1.873
0.982/1.000
82.736
PSO
Bus: 17; size: 0.1639 MWBus: 09; size: 0.9271 MWBus: 69; size: 0.0089 MWBus: 60; size: 2.0000 MW
3.011/1.945
0.982/1.000
83.398
SFO
Bus: 63; size: 1.7029 MWBus: 21; size: 0.0497 MWBus: 10; size: 1.3583 MWBus: 19; size: 0.4482 MW
4.648/3.000
0.986/1.005
81.975
SSA
Bus: 24; size: 0.6126 MWBus: 61; size: 0.8376 MWBus: 60; size: 1.6438 MWBus: 02; size: 0.0067 MW
4.599/2.961
0.994/1.005
89.803
COA
Bus: 18; size: 0.2706 MWBus: 23; size: 1.4078 MWBus: 63; size: 1.5098 MWBus: 55; size: 0.6024 MW
2.777/1.795
0.984/1.044
135.121
ICOA1
Bus: 65; size: 1.6407 MWBus: 55; size: 0.6677 MWBus: 23; size: 1.2326 MWBus: 17; size: 0.2436 MW
2.633/1.701
0.988/1.036
144.027
ICOA2
Bus: 23; size: 1.4105 MWBus: 18; size: 0.2700 MWBus: 63; size: 1.5370 MWBus: 55; size: 0.6570 MW
2.685/1.735
0.986/1.044
136.356
ICOA
Bus: 18; size: 0.2475 MWBus: 65; size: 1.5442 MWBus: 55; size: 0.6293 MWBus: 25; size: 1.3124 MW
2.709/1.750
0.985/1.047
156.052


## 6.1. Survey on Results Obtained by Setting Different Values ofPa for ICOA Method

For finding the most suitablePa value, ICOA is in turn run 50 times corresponding to each Pa value by setting Pa to 0.2, 0.4, 0.6, 0.8, and 1.0. The second case of minimizing the rated power of PVDGUs is selected for both IEEE 33-bus distribution network and IEEE 69-bus distribution network. The obtained results for the two systems are shown in Tables 2 and 3, respectively. As observed from Table 2, the best fitness corresponding to Pa = 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. Through the results, it points out that ICOA method can reach the best solution only at Pa = 0.2. Additionally, the average values of fitness function also confirm the most effective impact of Pa = 0.2 on its average fitness is the smallest among five values. The best fitness can conclude Pa = 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 Pa = 0.2. One more time, the significance of setting Pa = 0.2 can be adopted since the best fitness and the average fitness at Pa = 0.2 continue to be the smallest values among five values for the second system shown in Table 3. The best fitness and the best average are, respectively, 0.2425 and 0.2939. Consequently, we accept the best Pa value is 0.2 and it is used for other remaining cases with intent to shorten simulation time.Table 2
The obtained results by ICOA for Case2 of the IEEE 33-bus distribution network.
P
a
0.0
0.2
0.4
0.6
0.8
1.0
Best FF2
0.2860
0.2836
0.2978
0.2981
0.2900
0.2978
Average FF2
0.3881
0.3410
0.3474
0.3485
0.3451
0.3450
Worst FF2
0.5831
0.4109
0.4282
0.4326
0.4429
0.4448Table 3
The obtained results by ICOA for Case2 of the IEEE 69-bus distribution network.
P
a
0.0
0.2
0.4
0.6
0.8
1.0
Best FF2
0.2426
0.2425
0.2425
0.2426
0.2426
0.2426
Average FF2
0.3285
0.2939
0.2969
0.3130
0.3165
0.3271
Worst FF2
0.5143
0.4290
0.4716
0.5002
0.5104
0.5144


## 6.2. 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)
The best solution of the proposed method and compared methods was compared to show the strong search ability of the proposed ICOA.(2)
The 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)
The number of better runs (Nbetterrun) 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. The three factors are obtained by using the following equations:(45)IL of the best fitness%=the best fitness of the compared method−the best fitness of ICOAthe best fitness of the compared method×100,(46)IL of average fitness%=average fitness of the compared method−average fitness of ICOAaverage fitness of the compared method×100,(47)IL of the worst fitness%=the worst fitness of the compared method−the worst fitness of ICOAthe worst fitness of the compared method×100.(5)
Fair control parameter setting for the proposed ICOA and other compared methods was considered during the process of selecting population size and the maximum number of iterations.The implemented methods in this study are based on metaheuristic algorithms. Therefore, 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.


## 6.3. 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 FF1 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. The best values indicate that ICOA1 and ICOA can find the same best fitness less than that of COA and ICOA2. The average value of ICOA1 is better than that of COA but that of ICOA2 is worse than that of COA. Thus, 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. Figure4 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. The 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. There 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%. Figures4 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. Thus, 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.Table 4
The comparison among nine methods of the IEEE 33-bus distribution network for the first case.
Method
Best FF1
Average FF1
Worst FF1
IL of the best FF1
IL of average FF1
IL of the worst FF1
N
betterruns
BBO
0.3465
0.3782
0.4128
0.115
3.860
4.239
15
GA
0.3890
0.5027
0.6106
11.028
27.671
35.260
45
PSO
0.3697
0.4184
0.4492
6.384
13.098
11.999
27
SFO
0.3793
0.4184
0.446
8.753
13.098
11.368
38
SSA
0.3617
0.4054
0.4589
4.313
10.311
13.859
27
COA
0.3462
0.3727
0.4086
0.029
2.442
3.255
13
ICOA1
0.3461
0.3674
0.3963
0.000
1.034
0.252
8
ICOA2
0.3462
0.3779
0.4148
0.029
3.784
4.701
13
ICOA
0.3461
0.3636
0.3953
—
—
—
—Figure 4
The fitness values of four methods in 50 trial times for the first case.Figure 5
The fitness values of six methods in 50 trial times for the first case.Table 5
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case1.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 24; size: 1.0289 MWBus: 14; size: 0.8435 MWBus: 30; size: 1.0669 MW
4.612/2.994
0.969/1.000
73.108
GA
Bus: 13; size: 0.7711 MWBus: 06; size: 1.6119 MWBus: 31; size: 0.4229 MW
4.404/2.859
0.968/1.000
82.039
PSO
Bus: 14; size: 1.0746 MWBus: 25; size: 0.9004 MWBus: 30; size: 0.9302 MW
4.241/2.752
0.968/1.000
77.999
SFO
Bus: 24; size: 0.9188 MWBus: 09; size: 1.4072 MWBus: 30; size: 0.7580 MW
4.617/2.996
0.967/1.000
80.018
SSA
Bus: 14; size: 0.8744 MWBus: 30; size: 0.8465 MWBus: 24; size: 0.7361 MW
4.595/2.982
0.962/1.000
76.309
COA
Bus: 30; size: 1.0341 MWBus: 14; size: 0.8397 MWBus: 24; size: 1.1157 MW
4.620/2.998
0.969/1.000
73.038
ICOA1
Bus: 30; size: 1.0453 MWBus: 14; size: 0.8400 MWBus: 24; size: 1.0903 MW
4.619/2.998
0.969/1.000
73.031
ICOA2
Bus: 30; size: 1.0345 MWBus: 14; size: 0.8389 MWBus: 24; size: 1.1275 MW
4.621/2.999
0.969/1.000
73.043
ICOA
Bus: 24; size: 1.0885 MWBus: 30; size: 1.0467 MWBus: 14; size: 0.8381 MW
4.622/3.000
0.969/1.000
73.017Case 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 FF2 in equation (41) is used to evaluate the effectiveness of solutions from implemented approaches.
The 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 Table6. 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. The 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 Case1 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. The 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. The 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 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.Table 6
The comparison among nine methods of the IEEE 33-bus distribution network for the second case.
Method
Best FF2
Average FF2
Worst FF2
IL of the best FF2
IL of average FF2
IL of the worst FF2
N
betterruns
BBO
0.3024
0.3724
0.4553
6.217
8.432
9.752
5
GA
0.3582
0.5401
0.6977
20.826
36.864
41.106
41
PSO
0.3321
0.4027
0.5327
14.604
15.322
22.865
16
SFO
0.3702
0.4933
0.6259
23.393
30.874
34.351
43
SSA
0.3467
0.4474
0.5879
18.200
23.782
30.107
39
COA
0.3008
0.3633
0.4448
5.718
6.138
7.621
3
ICOA1
0.2975
0.3444
0.4228
4.672
0.987
2.815
2
ICOA2
0.3008
0.3746
0.4466
5.718
8.970
7.994
3
ICOA
0.2836
0.3410
0.4109
—
—
—
—Figure 6
The fitness values of four methods in 50 trial times for the second case.Figure 7
The fitness values of six methods in 50 trial times for the second case.Table 7
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case2.
Method
The best solution
THD/IHD (max) with PVDGUs
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 32; size: 0.5377 MWBus: 14; size: 0.6255 MWBus: 26; size: 0.0240 MW
4.531/2.941%
0.950/1.000
104.587
GA
Bus: 32; size: 0.2652 MWBus: 30; size: 0.2320 MWBus: 14; size: 0.9090 MW
4.604/2.988%
0.950/1.000
100.328
PSO
Bus: 14; size: 0.7796 MWBus: 31; size: 0.0495 MWBus: 32; size: 0.4747 MW
4.592/2.980%
0.951/1.000
101.133
SFO
Bus: 29; size: 0.2121 MWBus: 15; size: 0.8769 MWBus: 31; size: 0.3643 MW
4.269/2.771 pu
0.952/1.000
98.797
SSA
Bus: 32; size: 0.4525 MWBus: 2; size: 0.0004 MWBus: 14; size: 0.9021 MW
4.620/2.999 pu
0.950/1.000
102.365
COA
Bus: 31; size: 0.1411 MWBus: 14; size: 0.6197 MWBus: 33; size: 0.4200 MW
4.439/2.881 pu
0.950/1.000
104.925
ICOA1
Bus: 33; size: 0.2049 MWBus: 15; size: 0.5772 MWBus: 31; size: 0.3859 MW
4.264/2.767 pu
0.950/1.000
105.181
ICOA2
Bus: 31; size: 0.1505 MWBus: 14; size: 0.6201 MWBus: 33; size: 0.4105 MW
4.393/2.851 pu
0.950/1.000
104.887
ICOA
Bus: 33; size: 0.2781 MWBus: 31; size: 0.3144 MWBus: 15; size: 0.5210 MW
4.617/2.996 pu
0.950/1.000
109.299Case 3.
: Improving the voltage profile
In the part, voltage of all loads is considered to be improved by using fitness function FF3 in equation (42) as a comparison criterion for evaluating applied methods. In the fitness function, objective function F3 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 Table8. The 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. The lowest average fitness function and the lowest worst fitness function confirm that the proposed method is the most stable tool with the lowest fluctuations. The improvement level of the best fitness is from 0.531% to 17.590%. The 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 Table9.Table 8
The comparison among nine methods of the IEEE 33-bus distribution network for the third case.
Method
Best FF3
Average FF3
Worst FF3
IL of the best FF3
IL of average FF3
IL of the worst FF3
N
betterruns
BBO
0.0979
0.1273
0.159
4.290
9.112
8.805
4
GA
0.1137
0.2302
0.3099
17.590
49.739
53.211
28
PSO
0.0989
0.134
0.163
5.258
13.657
11.043
5
SFO
0.0990
0.1886
0.3064
5.354
38.653
52.676
5
SSA
0.1035
0.1634
0.2913
9.469
29.192
50.223
15
COA
0.0951
0.1222
0.1473
1.472
5.319
1.561
3
ICOA1
0.0942
0.1178
0.1451
0.531
1.783
0.069
3
ICOA2
0.0953
0.1228
0.1493
1.679
5.782
2.880
3
ICOA
0.0937
0.1157
0.145
—
—
—
—Figure 8
The fitness values of seven methods in 50 trial times for the third case.Table 9
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case3.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 07; size: 1.0925 MWBus: 30; size: 1.6785 MWBus: 14; size: 0.9046 MW
4.3764/2.840
0.983/1.005
114.7856
GA
Bus: 29; size: 1.9804 MWBus: 15; size: 0.3772 MWBus: 09; size: 1.2137 MW
4.625/3.000
0.983/1.007
124.438
PSO
Bus: 14; size: 0.8685 MWBus: 30; size: 1.4705 MWBus: 07; size: 1.3455 MW
4.437/2.880
0.983/1.002
109.115
SFO
Bus: 30; size: 1.6126 MWBus: 06; size: 1.1535 MWBus: 14; size: 0.9324 MW
4.3308/2.811
0.9835/1.004
113.053
SSA
Bus: 07; size: 0.7343 MWBus: 30; size: 1.7797 MWBus: 14; size: 1.0209 MW
4.208/2.731
0.9829/1.005
116.778
COA
Bus: 14; size: 0.9694 MWBus: 07; size: 1.0156 MWBus: 30; size: 1.6609 MW
4.279/2.777
0.9833/1.004
114.890
ICOA1
Bus: 14; size: 0.9464 MWBus: 07; size: 1.1930 MWBus: 30; size: 1.5831 MW
4.309/2.797
0.9836/1.005
114.680
ICOA2
Bus: 14; size: 0.9662 MWBus: 07; size: 1.0155 MWBus: 30; size: 1.6707 MW
4.283/2.779
0.9833/1.004
115.289
ICOA
Bus: 14; size: 0.9564 MWBus: 07; size: 1.1461 MWBus: 30; size: 1.5730 MW
4.298/2.789
0.9834/1.0043
113.132Case 4.
: Reducing the harmonic distortions
In the part, THD and IHD of all buses are improved by using fitness function FF4 in equation (43) to compare the quality of obtained solutions from implemented methods. In the fitness function, objective function F4 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 Table10 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. The 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%. The 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.
Optimal solutions together with IHD, THD, voltage profile, and total rated power of all PVDGUs obtained by these methods are shown in Table11.Table 10
The comparison among nine methods of the IEEE 33-bus distribution network for the fourth case.
Method
Best FF4
Average FF4
Worst FF4
IL of the best FF4
IL of average FF4
IL of the worst FF4
N
betterruns
BBO
1.1368
1.2524
1.428
0.176
7.833
5.469
26
GA
1.1829
1.4666
1.629
4.066
21.294
17.133
45
PSO
1.1392
1.2516
1.4708
0.386
7.774
8.220
33
SFO
1.1781
1.289
1.5417
3.675
10.450
12.441
45
SSA
1.1403
1.2571
1.4695
0.482
8.178
8.139
36
COA
1.1350
1.202
1.4074
0.018
3.968
4.086
4
ICOA1
1.1348
1.1593
1.3989
0.000
0.431
3.503
1
ICOA2
1.1350
1.2159
1.4619
0.018
5.066
7.661
4
ICOA
1.1348
1.1543
1.3499
—
—
—
—Figure 9
The fitness values of seven methods in 50 trial times for the fourth case.Table 11
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case4.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 19; size: 0.5830 MWBus: 15; size: 1.1293 MWBus: 33; size: 2.0001 MW
2.371/1.538
0.981/1.016
150.419
GA
Bus: 19; size: 0.5807 MWBus: 08; size: 1.6284 MWBus: 11; size: 1.5059 MW
2.699/1.751
0.961/1.022
165.954
PSO
Bus: 15; size: 1.8908 MWBus: 33; size: 1.2287 MWBus: 19; size: 0.5933 MW
2.444/1.585
0.981/1.035
164.398
SFO
Bus: 12; size: 1.7359 MWBus: 19; size: 0.5882 MWBus: 17; size: 1.0495 MW
2.709/1.757
0.955/1.050
218.820
SSA
Bus: 15; size: 1.2310 MWBus: 19; size: 0.5802 MWBus: 32; size: 1.8416 MW
2.404/1.559
0.9814/1.009
136.360
COA
Bus: 17; size: 1.1842 MWBus: 19; size: 0.5739 MWBus: 32; size: 1.9566 MW
2.403/1.558
0.982/1.018
154.611
ICOA1
Bus: 19; size: 0.5741 MWBus: 17; size: 1.1839 MWBus: 32; size: 1.9588 MW
2.402/1.558
0.982/1.018
154.760
ICOA2
Bus: 17; size: 1.1847 MWBus: 19; size: 0.5744 MWBus: 32; size: 1.9548 MW
2.401/1.558
0.982/1.018
154.512
ICOA
Bus: 19; size: 0.5779 MWBus: 17; size: 1.3284 MWBus: 32; size: 1.8086 MW
2.410/1.563
0.982/1.025
156.902


## 6.4. Test System 2: IEEE 69-Bus Distribution Network

In this section, an IEEE 69-bus distribution network is utilized as a main tool for running applied methods and the proposed ICOA method. Four single objectives consisting of total power losses, capacity of all PVDGUs, load voltage profile, and harmonic distortions similar to those in the previous section are in turn independently optimized. For better comparison of the performance of methods, values of fitness functions shown in equations (40)–(43) obtained from 50 runs are summarized in Tables 12–15 and reported in detail in Figures 10–13. Numerical tables can support enough comparison criteria such as quality of the best optimal solution, stability of all runs, and oscillations of runs while figures can enable a clear view of all runs.Table 12
The comparison among nine methods of the IEEE 69-bus distribution network for the first case.
Method
Best FF1
Average FF1
Worst FF1
IL of the best FF1
IL of average FF1
IL of the worst FF1
N
betterruns
BBO
0.3027
0.311
0.3402
0.264
1.383
7.613
8
GA
0.3132
0.3415
0.3968
3.608
10.190
20.791
48
PSO
0.3103
0.3267
0.3468
2.707
6.122
9.371
45
SFO
0.3131
0.3311
0.3632
3.577
7.369
13.464
48
SSA
0.3045
0.3219
0.3891
0.854
4.722
19.224
15
COA
0.3027
0.3085
0.3202
0.264
0.583
1.843
9
ICOA1
0.3024
0.3079
0.3184
0.165
0.390
1.288
3
ICOA2
0.3027
0.3091
0.3221
0.264
0.776
2.422
9
ICOA
0.3019
0.3067
0.3143
—
—
—
—Table 13
The comparison among nine methods of the IEEE 69-bus distribution network for the second case.
Method
Best FF2
Average FF2
Worst FF2
IL of the best FF2
IL of average FF2
IL of the worst FF2
N
betterruns
BBO
0.2427
0.3302
0.5152
0.082
10.993
16.731
1
GA
0.2752
0.4619
0.6518
11.882
36.372
34.182
17
PSO
0.2717
0.4102
0.592
10.747
28.352
27.534
17
SFO
0.3363
0.5192
0.6574
27.892
43.394
34.743
43
SSA
0.2563
0.4344
0.6854
5.384
32.343
37.409
12
COA
0.2427
0.3117
0.4857
0.082
5.711
11.674
3
ICOA1
0.2425
0.2966
0.4527
0.000
0.910
5.235
1
ICOA2
0.2427
0.3165
0.4991
0.082
7.141
14.045
3
ICOA
0.2425
0.2939
0.429
—
—
—
—Table 14
The comparison among nine methods of the IEEE 69-bus distribution network for the third case.
Method
Best FF3
Average FF3
Worst FF3
IL of the best FF3
IL of average FF3
IL of the worst FF3
N
betterruns
BBO
0.0571
0.0838
0.1631
16.462
5.728
20.969
5
GA
0.1122
0.2916
0.3915
57.487
72.908
67.075
44
PSO
0.1020
0.1991
0.3371
53.235
60.321
61.762
40
SFO
0.0829
0.2544
0.4031
42.461
68.947
68.023
34
SSA
0.0658
0.1584
0.3345
27.508
50.126
61.465
15
COA
0.0565
0.0881
0.1457
15.575
10.329
11.531
4
ICOA1
0.0517
0.0795
0.1302
7.737
0.629
0.998
2
ICOA2
0.0568
0.0895
0.1461
16.021
11.732
11.773
4
ICOA
0.0477
0.079
0.1289
—
—
—
—Table 15
The comparison among nine methods of the IEEE 69-bus distribution network for the fourth case.
Method
Best FF4
AverageFF4
Worst FF4
IL of the best FF4
IL of average FF4
IL of the worst FF4
N
betterruns
BBO
0.9547
1.0063
1.1189
0.293
0.646
5.523
4
GA
1.0282
1.0846
1.1964
7.421
7.819
11.643
44
PSO
1.0265
1.0767
1.1364
7.267
7.142
6.978
43
SFO
1.0064
1.0751
1.1269
5.415
7.004
6.194
29
SSA
0.9812
1.0705
1.1496
2.986
6.604
8.046
13
COA
0.9547
1.0029
1.0697
0.293
0.309
1.178
4
ICOA1
0.9526
1.0000
1.0572
0.073
0.020
0.009
2
ICOA2
0.9550
1.0048
1.0808
0.325
0.498
2.193
5
ICOA
0.9519
0.9998
1.0571
—
—
—
—Figure 10
The fitness values of seven methods in 50 trial times for the first case.Figure 11
The fitness values of seven methods in 50 trial times for the second case.Figure 12
The fitness values of seven methods in 50 trial times for the third case.Figure 13
The fitness values of seven methods in 50 trial times for the fourth case.In general, the proposed ICOA method outperforms the original COA method for all cases with high improvement level. The 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 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. The two values are 0.583% and 1.843% for Case1, 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.608% for Case 1, from 0.082% to 27.892% for Case 2, from 16.462% to 57.487% for Case 3, and from 0.293% to 7.421% for Case 4. Furthermore, the improvement levels of the average fitness and the worst fitness are much higher and can be up to 10.190% and 20.791% for Case 1, 43.394% and 37.409% for Case 2, 68.947% and 68.023% for Case 3, and 7.819% and 11.643% for Case 4. The number of runs from ICOA with better fitness than the best fitness of these compared methods is also reported. For almost cases, BBO is the second best method while GA, PSO, and SFO are the worst methods, especially GA. ICOA can find 8, 1, 5, and 4 solutions with better fitness than the best solutions of BBO for Cases 1, 2, 3, and 4, respectively. As compared to GA, PSO, and SFO, ICOA can reach much higher number of runs with better fitness. They are, respectively, 48, 45, and 48 for Case 1, 17, 17, and 43 for Case 2, 44, 40, and 34 for Case 3, and 44, 43, and 29 for Case 4. One more time, the quality of runs from ICOA can be confirmed by observing figures. 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 Tables16–19.Table 16
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case1.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 21; size: 0.3057 MWBus: 49; size: 0.7413 MWBus: 10; size: 0.5958 MWBus: 61; size: 1.7118 MW
4.260/2.753
0.979/1.000
67.951
GA
Bus: 39; size: 0.3089 MWBus: 17; size: 0.5140 MWBus: 61; size: 1.6919 MWBus: 53; size: 0.6702 MW
3.887/2.512
0.980/1.000
70.312
PSO
Bus: 08; size: 0.7504 MWBus: 49; size: 0.7036 MWBus: 61; size: 1.6883 MWBus: 17; size: 0.2579 MW
4.351/2.812
0.978/1.000
69.668
SFO
Bus: 61; size: 1.6826 MWBus: 39; size: 0.3095 MWBus: 17; size: 0.5122 MWBus: 53; size: 0.6785 MW
3.893/2.515
0.980/1.000
70.299
SSA
Bus: 17; size: 0.3835 MWBus: 49; size: 0.9538 MWBus: 61; size: 1.6617 MWBus: 66; size: 0.3690 MW
4.144/2.678
0.976/1.000
68.360
COA
Bus: 61; size: 1.7076 MWBus: 17; size: 0.3894 MWBus: 49; size: 0.9110 MWBus: 66; size: 0.4886 MW
4.115/2.659
0.978/1.000
67.954
ICOA1
Bus: 49; size: 0.7485 MWBus: 10; size: 0.5841 MWBus: 61; size: 1.7367 MWBus: 21; size: 0.3572 MW
4.160/2.688
0.980/1.000
67.889
ICOA2
Bus: 49; size: 0.7493 MWBus: 10; size: 0.5938 MWBus: 21; size: 0.3069 MWBus: 61; size: 1.7068 MW
4.259/2.752
0.979/1.000
67.949
ICOA
Bus: 49; size: 0.7183 MWBus: 21; size: 0.3231 MWBus: 61; size: 1.6881 MWBus: 11; size: 0.6243 MW
4.233/2.735
0.978/1.000
67.765Table 17
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case2.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 24; size: 0.0005 MWBus: 61; size: 0.4310 MWBus: 64; size: 0.5456 MWBus: 03; size: 0.0003 MW
4.453/2.875
0.950/1.000
112.543
GA
Bus: 17; size: 0.0947 MWBus: 64; size: 0.4732 MWBus: 63; size: 0.5328 MWBus: 42; size: 0.0075 MW
3.901/2.519
0.952/1.000
106.372
PSO
Bus: 68; size: 0.0159 MWBus: 69; size: 0.0477 MWBus: 61; size: 0.5273 MWBus: 64; size: 0.5028 MW
4.503/2.902
0.952/1.000
106.844
SFO
Bus: 21; size: 0.3091 MWBus: 61; size: 0.7208 MWBus: 63; size: 0.3180 MWBus: 64; size: 0.0602 MW
3.897/2.519
0.954/1.000
92.108
SSA
Bus: 64; size: 0.1936 MWBus: 24; size: 0.0361 MWBus: 63; size: 0.0633 MWBus: 62; size: 0.7390 MW
4.713/3.000
0.950/1.000
109.137
COA
Bus: 02; size: 0.0003 MWBus: 61; size: 0.4312 MWBus: 64; size: 0.5455 MWBus: 24; size: 0.0004 MW
4.457/2.878
0.950/1.000
112.332
ICOA1
Bus: 15; size: 0.0002 MWBus: 62; size: 0.4975 MWBus: 64; size: 0.4792 MWBus: 02; size: 0.0001 MW
4.557/2.937
0.950/1.000
112.644
ICOA2
Bus: 24; size: 0.0005 MWBus: 61; size: 0.4311 MWBus: 63; size: 0.5456 MWBus: 03; size: 0.0002 MW
4.452/2.875
0.950/1.000
112.543
ICOA
Bus: 15; size: 0.0002 MWBus: 62; size: 0.4976 MWBus: 64; size: 0.4791 MWBus: 02; size: 0.0001 MW
4.557/2.937
0.950/1.000
112.643Table 18
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case3.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 62; size: 2.0022 MWBus: 15; size: 0.5342 MWBus: 55; size: 0.7809 MWBus: 68; size: 0.4592 MW
3.725/2.399
0.994/1.003
85.109
GA
Bus: 11; size: 0.2395 MWBus: 65; size: 1.9581 MWBus: 56; size: 0.9560 MWBus: 17; size: 0.6186 MW
3.710/2.390
0.994/1.020
124.062
PSO
Bus: 65; size: 1.9993 MWBus: 09; size: 0.2583 MWBus: 17; size: 0.5867 MWBus: 56; size: 0.9538 MW
3.823/2.464
0.994/1.022
126.773
SFO
Bus: 55; size: 0.0799 MWBus: 53; size: 1.1452 MWBus: 18; size: 0.5801 MWBus: 61; size: 1.9702 MW
3.569/2.299
0.994/1.001
80.482
SSA
Bus: 60; size: 1.3378 MWBus: 66; size: 0.7844 MWBus: 17; size: 0.4479 MWBus: 63; size: 0.9821 MW
3.860/2.494
0.994/1.003
85.481
COA
Bus: 62; size: 2.0000 MWBus: 55; size: 0.7803 MWBus: 15; size: 0.5341 MWBus: 68; size: 0.4552 MW
3.726/2.400
0.994/1.002
84.943
ICOA1
Bus: 10; size: 0.7891 MWBus: 61; size: 2.0311 MWBus: 21; size: 0.3902 MWBus: 56; size: 0.5777 MW
3.960/2.559
0.994/1.002
82.884
ICOA2
Bus: 62; size: 2.0061 MWBus: 55; size: 0.7815 MWBus: 68; size: 0.4548 MWBus: 15; size: 0.5339 MW
3.724/2.399
0.994/1.003
85.177
ICOA
Bus: 10; size: 0.7786 MWBus: 21; size: 0.4094 MWBus: 56; size: 0.5689 MWBus: 61; size: 2.0001 MW
3.936/2.543
0.994/1.001
81.781Table 19
The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case4.
Method
The best solution
THD/IHD (max) with PVDGUs (%)
Volt (min/max) with PVDGUs (pu)
Total power losses with PVDGUs (kW)
BBO
Bus: 18; size: 0.2705 MWBus: 23; size: 1.4079 MWBus: 63; size: 1.5097 MWBus: 55; size: 0.6023 MW
2.778/1.795
0.984/1.044
135.120
GA
Bus: 66; size: 0.0256 MWBus: 09; size: 1.0095 MWBus: 60; size: 1.9590 MWBus: 17; size: 0.1785 MW
2.899/1.873
0.982/1.000
82.736
PSO
Bus: 17; size: 0.1639 MWBus: 09; size: 0.9271 MWBus: 69; size: 0.0089 MWBus: 60; size: 2.0000 MW
3.011/1.945
0.982/1.000
83.398
SFO
Bus: 63; size: 1.7029 MWBus: 21; size: 0.0497 MWBus: 10; size: 1.3583 MWBus: 19; size: 0.4482 MW
4.648/3.000
0.986/1.005
81.975
SSA
Bus: 24; size: 0.6126 MWBus: 61; size: 0.8376 MWBus: 60; size: 1.6438 MWBus: 02; size: 0.0067 MW
4.599/2.961
0.994/1.005
89.803
COA
Bus: 18; size: 0.2706 MWBus: 23; size: 1.4078 MWBus: 63; size: 1.5098 MWBus: 55; size: 0.6024 MW
2.777/1.795
0.984/1.044
135.121
ICOA1
Bus: 65; size: 1.6407 MWBus: 55; size: 0.6677 MWBus: 23; size: 1.2326 MWBus: 17; size: 0.2436 MW
2.633/1.701
0.988/1.036
144.027
ICOA2
Bus: 23; size: 1.4105 MWBus: 18; size: 0.2700 MWBus: 63; size: 1.5370 MWBus: 55; size: 0.6570 MW
2.685/1.735
0.986/1.044
136.356
ICOA
Bus: 18; size: 0.2475 MWBus: 65; size: 1.5442 MWBus: 55; size: 0.6293 MWBus: 25; size: 1.3124 MW
2.709/1.750
0.985/1.047
156.052


## 7. 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. Therefore, 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. These 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 solutionComid,g in equation (22) with the best solution shown in (29). Thus, 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. Thus, 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)
The 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. The ICOA method has three basic parameters includingNCo, Ng, and ITMax. In this method, the population size (Npop) and the number of new solutions for each iteration are equal to (NCo × Ng) and (NCo × Ng + Ng), respectively. During the simulation process, we surveyed and selected two parameters as NCo and Ng to find the most appropriate values for two considered systems. The values of NCo and Ng are adjusted to find the most suitable values. As observing obtained results, we have found that 4 is suitable for both NCo and Ng whereas 75 and 100 are, respectively, reasonable for ITMax for the first system and the second system. Therefore, it is not easy to choose the right parameters for each specific network. This 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.


## 8. 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 33-bus 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. The 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. The 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 onesConnecting 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, respectivelyConsequently, 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.



---

*Source: 1603802-2020-02-13.xml*