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

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 current

 

Case 2: minimizing the rated power of PVDGUs under the constraints of load bus voltage, branch current, and harmonic distortions

 

Case 3: improving the voltage profile under the constraints of harmonic distortions and branch current

 

Case 4: reducing the harmonic distortions under the constraint of voltage profile and branch current

In 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 N_g small coyote groups, and there are N_co coyotes in each group. COA method produces two new solution generations per iteration in which the first generation produces N_co new solutions for each group and the second generation produces N_g new solutions for the whole coyote community. Thus, total number of new solutions generated in the COA method is (N_co × N_g + N_g) solutions in which N_co × N_g is equal to population (N_pop). Thus, (N_co × N_g + N_g) is equal to (N_pop + N_g). 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 N_co times new solutions of the second generation. However, COA is coping with low performance of the first generation because it uses a center solution for producing updated step size; meanwhile, the strong point of the center solution is only to produce diversity for solutions, but it is not potential for generating promising updated step size. In the second generation, COA produces one new solution for each group by using a randomization factor. 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 Section 2. 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]:$$\begin{array}{}\text{(1)}& {\text{TPL}}_{\text{DG}}=\sum _{n=1}^{N_{\text{br}}}{I}_{n\text{DG}}^2{R}_n.\end{array}$$

By using the same method, total power losses in all branches before installing PVDGUs, TPL, is determined by means of the following equation [41]:$$\begin{array}{}\text{(2)}& \text{TPL}=\sum _{n=1}^{N_{\text{br}}}{I}_n^2{R}_n.\end{array}$$

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:$$\begin{array}{}\text{(3)}& \text{minimize\hspace{0.17em}}{F}_1=\frac{{\text{TPL}}_{\text{DG}}}{\text{TPL}}.\end{array}$$

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 (RP_DG) 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]:$$\begin{array}{}\text{(4)}& \text{minimize\hspace{0.17em}}{F}_2=\frac{{\sum }_{i=1}^{N_{\text{DG}}}{\text{RP}}_{\text{DG},i}}{{\sum }_{i=1}^{N_{\text{load}}}{P}_{\text{Load}i}}.\end{array}$$

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 the i^th bus can be defined as$$\begin{array}{}\text{(5)}& {\text{VP}}_i=\frac{\left(V_i-V^{\text{min}}\right)\cdot \left(V^{\max }-V_i\right)}{\left(V^{\text{nom}}-V^{\min }\right)\cdot \left(V^{\max }-V^{\text{nom}}\right)},\end{array}$$where V^nom is the nominal or desired bus voltage value and typically taken as 1.0 pu. It is clear that VP_i can get a maximum value (1.0 pu) if the bus voltage equals to V^nom and VP_i can get a negative value if the bus voltage is higher than the maximum limit or lower than the minimum limit. The overall voltage profile index for the distribution system can be defined as$$\begin{array}{}\text{(6)}& \text{maximum\hspace{0.17em}VPI}=\frac{1}{N_{\text{bus}}}\cdot \sum _{i=1}^{N_{\text{bus}}}{\text{VP}}_i.\end{array}$$

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:$$\begin{array}{}\text{(7)}& \text{minimize\hspace{0.17em}}{F}_3=1-\text{VPI}.\end{array}$$

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]$$\begin{array}{}\text{(8)}& {\text{THD}}_{V,i}\left(\%\right)=\frac{\sqrt{{\sum }_{h\ne 1}^H{\left|V_i^h\right|}^2}}{\left|V_i^1\right|}\cdot 100.\end{array}$$

Individual harmonic distortion at each bus can be defined by [44]$$\begin{array}{}\text{(9)}& {\text{IHD}}_{V,i}^h\left(\%\right)=\frac{\left|V_i^h\right|}{\left|V_i^1\right|}\cdot 100.\end{array}$$

According to IEEE Std. 519 [44], ${\text{THD}}_{V,i}$ and ${\text{IHD}}_{V,i}^h$ should not be greater than 5% and 3%, respectively.

The average of ${\text{THD}}_{V,i}$ can be defined by [25]$$\begin{array}{}\text{(10)}& {\text{THD}}_{\text{Vaver}}=\frac{{\sum }_{i=1}^{N_{\text{bus}}}{\text{THD}}_{V,i}}{N_{\text{bus}}}.\end{array}$$

Similarly, the average of ${\text{IHD}}_{V,i}^h$ can be defined by$$\begin{array}{}\text{(11)}& {\text{IHD}}_{\text{Vaver}}^h=\frac{{\sum }_{i=1}^{N_{\text{bus}}}{\text{IHD}}_{V,i}^h}{N_{\text{bus}}}.\end{array}$$

Here, the objective function should be shown as follows:$$\begin{array}{}\text{(12)}& \text{minimize\hspace{0.17em}}{F}_4={\omega }_1\cdot {\text{THD}}_{\text{Vaver}}+{\omega }_2\cdot \max \left({\text{IHD}}_{\text{Vaver}}^h\right),\end{array}$$where ${\omega }_1$ and ${\omega }_2$ are, respectively, weighted factors with ${\text{THD}}_{\text{Vaver}}$ and ${\text{IHD}}_{\text{Vaver}}^h$.${\omega }_1$ and ${\omega }_2$ are constrained by$$\begin{array}{c}\text{(13)}& {\omega }_1+{\omega }_2=1,0\le {\omega }_1,\\ {\omega }_2\le 1.\end{array}$$

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:$$\begin{array}{}\text{(14)}& \sum _{i=1}^{N_{\text{load}}}{P}_{\text{Load},i}+\sum _{n=1}^{N_{\text{br}}}{P}_{\text{Loss},n}-\sum _{i=1}^{N_{\text{DG}}}{\text{RP}}_{\text{DG},i}-P_{\text{Grid}}=0.\end{array}$$

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:$$\begin{array}{}\text{(15)}& V^{\min }\le \sqrt{{\left|V_i^1\right|}^2+\sum _{h=2}^H{\left|V_i^h\right|}^2}\le V^{\max },\hspace{1em}i=1,\dots ,N_{\text{bus}}.\end{array}$$

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]:$$\begin{array}{c}\text{(16)}& {\text{THD}}_{V,i}\le {\text{THD}}^{\max },{\text{IHD}}_{V,i}^h\le {\text{IHD}}^{\max },\end{array}$$where ${\text{THD}}^{\max }$ and ${\text{IHD}}^{\max }$ 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:$$\begin{array}{c}\text{(17)}& {\text{RP}}_{\text{DG}}^{\min }\le {\text{RP}}_{\text{DG},i}\le {\text{RP}}_{\text{DG}}^{\max },\sum _{i=1}^{N_{\text{DG}}}{\text{RP}}_{\text{DG},i}\le \sum _{i=1}^{N_{\text{load}}}{\text{RP}}_{\text{Load},i},\end{array}$$where ${\text{RP}}_{\text{DG}}^{\min }$ and ${\text{RP}}_{\text{DG}}^{\max }$ 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:$$\begin{array}{}\text{(18)}& \left|I_{n\text{DG}}\right|\le I_n^{\max },\end{array}$$where $I_{n\text{DG}}$ and $I_n^{\max }$ are the current magnitude of n^th 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 into N_g groups where each group has N_Co coyotes. Thus, population size of COA is N_pop (where N_pop = N_Co × N_g).

(3)

The k^th coyote in the $g^{\text{th}}$ group has two important factors consisting of social condition ${\text{Co}}_{k,g}$ and quality ${\text{FF}}_{k,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:$$\begin{array}{}\text{(19)}& {\text{Co}}_{k,g}={\text{Co}}^{\min }+\gamma \cdot \left({\text{Co}}^{\max }-{\text{Co}}^{\min }\right);\hspace{1em}k=1,\dots ,N_{\text{Co}},g=1,\dots ,N_g,\end{array}$$where Co^min and Co^max are, respectively, lower bound solution and upper bound solution obtained by$$\begin{array}{c}\text{(20)}& {\text{Co}}^{\min }=\left[{\mathrm{var}}_x^{\min }\right];\hspace{1em}x=1,\dots ,N_{\text{dv}},{\text{Co}}^{\max }=\left[{\mathrm{var}}_x^{\max }\right];\hspace{1em}x=1,\dots ,N_{\text{dv}},\end{array}$$where ${\mathrm{var}}_x^{\min }$ and ${\mathrm{var}}_x^{\max }$ are the minimum value and maximum value of the x^th control variable, respectively.

Step 2.

Determine fitness function value

Social conditions of coyote community are evaluated. 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:$$\begin{array}{}\text{(21)}& {\text{FF}}_{k,g}={\text{OF}}_{k,g}+\text{PF}\cdot \sum _{m=1}^M{\left({\text{PT}}_{k,g,m}\right)}^2,\end{array}$$where ${\text{OF}}_{k,g}$ is the objective function of the k^th solution in the $g^{\text{th}}$ group; PF is a penalty factor, which is selected by experience; ${\text{PT}}_{k,g,m}$ is the penalty term for violating the m^th constraint of the k^th solution in the $g^{\text{th}}$ group; and M is the number of constraints that are taken into account in the considered problem.

In each group, the solution with the best fitness function is called the local best solution and symbolized by ${\text{Co}}_{\text{best},g}$.

Step 3.

Update solutions for each group

Updating new solutions is corresponding to updating social condition for coyotes in each group. In COA, new solutions are found around old solutions with two distances based on the local best solution, two randomly picked solutions and a middle solution. The detail of the newly updated solution strategy can be expressed by the following model:$$\begin{array}{}\text{(22)}& {\text{Co}}_{k,g}^{\text{new}}={\text{Co}}_{k,g}+\gamma \left({\text{Co}}_{\text{best},g}-{\text{Co}}_{1.g}\right)+\gamma \left({\text{Co}}_{\text{mid},g}-{\text{Co}}_{2,g}\right).\end{array}$$

In equation (25), ${\text{Co}}_{1,g}$ and ${\text{Co}}_{2,g}$ are two randomly picked solutions in the $g^{\text{th}}$ considered group while ${\text{Co}}_{\text{mid},g}$ is the middle solution obtained from the $g^{\text{th}}$ group. Each decision variable in the middle solution is obtained by one out of two ways shown in the following equation:$$\begin{array}{}\text{(23)}& {\text{Co}}_{\text{mid},g}=\left\{\begin{array}{cc}{\mathrm{var}}_{x,\text{mid}1},& \text{if\hspace{0.17em}}{N}_{\text{Co}}\text{\hspace{0.17em}is\hspace{0.17em}odd},\\ {\mathrm{var}}_{x,\text{mid}2},& \text{else,}\\ x=1,\dots N_{\text{dv}},& \end{array}\right.\end{array}$$where var_x,mid1 and var_x,mid2 are the x^th center decision variables corresponding to two different cases of N_Co. If N_Co is an odd number, var_x,mid1 is taken from the x^th decision variable at the position [(N_Co + 1)/2]. On the contrary, var_x,mid2 is the average of two x^th decision variables at two positions, N_Co/2 and [(N_Co + 1)/2]. However, before determining the two variables in equation (23), each decision variable of all available solutions is sorted in descending order.

Step 4.

Evaluate the effectiveness level of new solutions

The new solutions are evaluated by calculating fitness function (21).

Step 5.

Perform selection technique

Now, each coyote k has two social conditions, old social condition corresponding to old solution, ${\text{Co}}_{k,g}$, and new social condition corresponding to new solution, ${\text{Co}}_{k,g}^{\text{new}}$. The quality of the two social conditions is the fitness function of the two solutions, ${\text{FF}}_{k,g}$ and ${\text{FF}}_{k,g}^{\text{new}}$. So, it should retain one social condition for each k^th coyote by using the following rules:$$\begin{array}{c}\text{(24)}& {\text{Co}}_{k,g}=\left\{\begin{array}{cc}{\text{Co}}_{k,g}^{\text{new}},& \text{if\hspace{0.17em}}{\text{FF}}_{k,g}^{\text{new}}\le {\text{FF}}_{k,g},\\ {\text{Co}}_{k,g},& \text{else},\end{array}\right.{\text{FF}}_{k,g}=\left\{\begin{array}{cc}{\text{FF}}_{k,g}^{\text{new}},& \text{if\hspace{0.17em}}{\text{FF}}_{k,g}^{\text{new}}\le {\text{FF}}_{k,g},\\ {\text{FF}}_{k,g},& \text{else}.\end{array}\right.\end{array}$$

Step 6.

Create one new social condition in each group

In the step, each group $g^{\text{th}}$ produces one new solution ${\text{Co}}_g^{\text{new}}$ by using a randomization mechanism. The new solution and the technique are shown in the following models:$$\begin{array}{}\text{(25)}& C{o}_g^{\text{new}}=\left[{\mathrm{var}}_{x,g}^{\text{new}}\right];\hspace{1em}x=1,\dots ,N_{\text{dv}},\text{(26)}& {\mathrm{var}}_{x,g}^{\text{new}}=\left\{\begin{array}{cc}{\mathrm{var}}_{x,1,g}^{\text{new}},& \text{if\hspace{0.17em}}\beta <\frac{1}{N_{\text{dv}}},\\ {\mathrm{var}}_{x,2,g}^{\text{new}},& \text{if\hspace{0.17em}}\frac{1}{N_{\text{dv}}}\le \beta <\frac{1}{N_{\text{dv}}}+0.5,\\ {\mathrm{var}}_{x,g}^{\text{new}},& \text{else},\end{array}\right.\end{array}$$where ${\mathrm{var}}_{x,1,g}^{\text{new}}$ and ${\mathrm{var}}_{x,2,g}^{\text{new}}$ are two x^th variables randomly picked from two available solutions; ${\mathrm{var}}_{x,g}^{\text{new}}$ 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 $F{F}_g^{\text{new}}$.

Step 7.

Determine the worst solution and replace it

Each new solution in each group obtained by Step 6 will be compared to the worst solution in each group ${\text{Co}}_g^{\text{worst}}$ for deciding if the worst solution should be replaced with the new one. The decision is taken over by the comparison below:$$\begin{array}{}\text{(27)}& {\text{Co}}_g^{\text{worst}}=\left\{\begin{array}{cc}{\text{Co}}_g^{\text{worst}},& \text{if\hspace{0.17em}}{\text{FF}}_g^{\text{worst}}<{\text{FF}}_g^{\text{new}},\\ {\text{Co}}_g^{\text{new}},& \text{else}.\end{array}\right.\end{array}$$

Step 8.

Exchange solutions among N_g groups

To diversify solutions and avoid falling into local optimal zones with premature convergence, COA carries out solution exchange action among N_g groups. Two randomly selected groups will provide two randomly picked solutions and their position is then exchanged. The task is dependent on the randomization mechanism below:$$\begin{array}{}\text{(28)}& \phi <\frac{{10}^{-2}}{2}{N}_{\text{Co}}^2,\end{array}$$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 N_Co is, the higher the probability of exchanged solutions is.

Step 9.

Determine the best solution Co_Gbest

Before terminating one computation iteration, the best solution among N_g groups is determined by comparing fitness function. If the current iteration is the last one, the best solution Co_Gbest 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 (N_Co × N_g + 1 × N_g) equaling to (N_pop + N_g). 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 solution ${\text{Co}}_{\text{mid},g}$ in equation (22) is selected to be replaced with the best solution in the current population. In COA, the impact of the center solution is significant in finding optimal solution as solving benchmark optimization functions. 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:$$\begin{array}{}\text{(29)}& {\text{Co}}_{k,g}^{\text{new}}={\text{Co}}_{k,g}+\left[\gamma \cdot \left({\text{Co}}_{\text{best},g}-{\text{Co}}_{1,g}\right)+\gamma \cdot \left({\text{Co}}_{\text{Gbest}}-{\text{Co}}_{2,g}\right)\right].\end{array}$$

By using the new strategy, the proposed method can reduce several shortcomings:

(1)

Reducing the task of sorting each variable from all solutions: as pointed in the COA method, each solution has N_dv decision variables and the x^th decision variable of all solutions must be grouped and then sorted. Finally, equation (23) is applied for determining the most appropriate x^th decision variable. The task becomes time-consuming if N_dv is high. Furthermore, if the number of groups is high, the number of center solutions is also high, resulting in long execution time for the center solutions.

(2)

Finding higher quality solutions: the best solution has a high contribution to the level of solution quality improvement, so the proposed method can help to reduce the number of iterations and shorten execution time.

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 N_dv is used in order to establish three selection conditions in which the condition is random and selected variables are also random. Furthermore, random variables within upper bound and lower bound are produced for the third selection condition. In the first two conditions, two variables are taken from two random solutions in the current population while random variables within predetermined boundaries are used in the third condition. Clearly, the combination of the variables cannot form good solutions with high quality. So, in the second improvement, we propose two models for updating new solutions based on the following equations:$$\begin{array}{}\text{(30)}& {\text{CE}}_g={\text{Co}}_{\text{Gbest}}+{\gamma }^{\ast }\left({\text{Co}}_{\text{Gbest}}-{\text{Co}}_{\text{best},1}\right)+{\gamma }^{\ast }\left({\text{Co}}_{\text{Gbest}}-{\text{Co}}_{\text{best,}2}\right),\text{(31)}& {\text{CE}}_g={\text{Co}}_{\text{Gbest}}+\gamma \ast \left({\text{Co}}_{\text{Gbest}}-{\text{Co}}_{\text{best,}1}\right)+{\gamma }^{\ast }\left({\text{Co}}_{\text{Gbest}}-{\text{Co}}_{\text{best,}2}\right)+{\gamma }^{\ast }\left({\text{Co}}_{\text{Gbest}}-{\text{Co}}_{\text{best,}3}\right).\end{array}$$

As seen from the two equations above, (31) can generate a larger step and new solution found by using (31) has a clearer difference. However, it is not sure that equation (31) can produce solutions with better fitness if there is no measurement technique. Let suppose that if the first step size between Co_Gbest and Co_best,1 and the second step size between Co_Gbest and Co_best,2 are small, the new solution obtained by equation (30) is so close to the old solution and local optimal solution search strategy may be ineffective. Furthermore, if the phenomenon is repeated in equation (31), the third step size will be high and the new solution is fallen into other zones, which are very far from the old zone. So, if the use is appropriate, equation (31) can prevent the effectiveness of the local search strategy. 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 Count_max and R, which are obtained by equations (32) and (33). As indicated by equation (33), R is within 0 and 1, and its value is directly proportional to the number of pairs of close solutions. If R is close to zero, it means almost all solutions are not close but if R is around 1, approximately the considered solutions are very close together. So, it needs to determine a probability P_a for each specific optimization problem. As P_a is large enough, equation (31) providing large step size is applied. Otherwise, equation (30) is preferred. In numerical result, we will test the sensitiveness of P_a on the obtained result and the most appropriate value of P_a will be recommended.$$\begin{array}{}\text{(32)}& {\text{Count}}_{\max }=\frac{N_{\text{Co}}\times \left(N_{\text{Co}}+1\right)}{2},\text{(33)}& R=\frac{\text{Count}}{{\text{Count}}_{\max }}.\end{array}$$

Algorithm 1: The proposed technique for the second new solution generation.

 

Count = 0;

 

Calculate Count_max by using equation (32)

 

for i = 1 to (N_Co − 1)

 

for j = i + 1 to N_Co

 

$\text{if\hspace{0.17em}}\left|{\text{FF}}_{j,g}-{\text{FF}}_{i,g}\right|<\left|{\text{FF}}_{\text{mean},g}-{\text{FF}}_{\text{best},g}\right|$

 

Count = Count + 1

 

end

 

end

 

end

 

Calculate R by using equation (33)

 

If R > P_a

 

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 Table 1. 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 j^th solar photovoltaic distribution generation unit (PVDGU) corresponding to the k^th solution in the $g^{\text{th}}$ group are, respectively, represented by ${\text{Site}}_{j,k,g}$ and ${\text{RP}}_{j,k,g}$. Each solution in each group is represented by the following equation:$$\begin{array}{}\text{(34)}& {\text{Co}}_{g,k}=\left[{\text{Site}}_{j,g,k},{\text{RP}}_{j,g,k}\right];\hspace{1em}j=1,\dots ,N_{\text{DG}},\end{array}$$where N_DG 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; −253

In order to produce the initial population, lower bound solution Co^min and upper bound solution Co^max should be defined to make sure that randomly generated solutions can satisfy the limit constraints:$$\begin{array}{c}\text{(35)}& {\text{Co}}^{\min }=\left[{\text{Site}}_j^{\min },{\text{RP}}_j^{\min }\right];\hspace{1em}j=1,\dots ,N_{\text{DG}},{\text{Co}}^{\max }=\left[{\text{Site}}_j^{\max },{\text{RP}}_j^{\max }\right];\hspace{1em}j=1,\dots ,N_{\text{DG}},\end{array}$$where ${\text{Site}}_j^{\min }$ and ${\text{Site}}_j^{\max }$ are, respectively, bus 2 and bus N_bus and ${\text{RP}}_j^{\min }$ and ${\text{RP}}_j^{\max }$ are the minimum rated power and maximum rated power of the j^th PVDG unit. However, the minimum power of all units and the maximum power of all units are the same and corresponding to ${\text{RP}}^{\min }$ and ${\text{RP}}^{\max }$. 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:$$\begin{array}{}\text{(36)}& {\text{Co}}_{k,g}={\text{Co}}^{\min }+\epsilon \cdot \left({\text{Co}}^{\max }-{\text{Co}}^{\min }\right);\hspace{1em}k=1,\dots ,N_{\text{Co}};g=1,\dots ,N_g.\end{array}$$

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 currents $I_{n\text{DG}}^h$ and load bus voltage $V_i^h$ 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, I_nDG, and voltage at bus i, V_i, are calculated by using the two following equations:$$\begin{array}{c}\text{(37)}& I_{n\text{DG}}=\sqrt{{\left|I_{n\text{DG}}^1\right|}^2+\sum _{h=2}^H{\left|I_{n\text{DG}}^h\right|}^2},V_i=\sqrt{{\left|V_i^1\right|}^2+\sum _{h=2}^H{\left|V_i^h\right|}^2}.\end{array}$$

As a result, the two parameters are checked and penalized by using the following formulas:$$\begin{array}{c}\text{(38)}& \Delta V_i=\left\{\begin{array}{cc}{V}_i-V_i^{\max },& \text{if\hspace{0.17em}}{V}_i>V^{\max },\\ V_i^{\min }-V_i,& \text{if\hspace{0.17em}}{V}_i<V^{\min },\\ 0,& \text{else,}\end{array}\right.\Delta I_n=\left\{\begin{array}{cc}{I}_{n\text{DG}}-I_n^{\max },& \text{if\hspace{0.17em}}{I}_{n\text{DG}}>I_n^{\max },\\ 0,& \text{else}.\end{array}\right.\end{array}$$

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:$$\begin{array}{c}\text{(39)}& \Delta {\text{THD}}_i=\left\{\begin{array}{cc}{\text{THD}}_i-{\text{THD}}^{\max },& \text{if\hspace{0.17em}}{\text{THD}}_i>{\text{THD}}^{\max },\\ 0,& \text{else,}\end{array}\right.\Delta {\text{IHD}}_i=\left\{\begin{array}{cc}{\text{IHD}}_i-{\text{IHD}}^{\max },& \text{if\hspace{0.17em}}{\text{IHD}}_i>{\text{IHD}}^{\max },\\ 0,& \text{else}.\end{array}\right.\end{array}$$

5.3. Fitness Function

As described in Section 2, we focus on four single objectives consisting of total power loss minimization F₁, rated power minimization of PVDGUs F₂, load bus voltage enhancement F₃, and harmonic distortion minimization F₄. According to the considered objectives, fitness function can be established as follows:

(i)

Fitness function for total power loss minimization:

$$\begin{array}{}\text{(40)}& {\text{FF}}_1=F_1+{\mu }_V\sum _{i=1}^{N_{\text{bus}}}\Delta V_i^2+{\mu }_I\sum _{i=1}^{N_{\text{br}}}\Delta I_n^2+{\mu }_{\text{THD}}\sum _{i=1}^{N_{\text{bus}}}\Delta {\text{THD}}_i^2+{\mu }_{\text{IHD}}\sum _{i=1}^{N_{\text{bus}}}\Delta {\text{IHD}}_i^2.\end{array}$$

(ii)

Fitness function for rated power minimization of PVDGUs:

$$\begin{array}{}\text{(41)}& {\text{FF}}_2=F_2+{\mu }_V\sum _{i=1}^{N_{\text{bus}}}\Delta V_i^2+{\mu }_I\sum _{i=1}^{N_{\text{br}}}\Delta I_n^2+{\mu }_{\text{THD}}\sum _{i=1}^{N_{\text{bus}}}\Delta {\text{THD}}_i^2+{\mu }_{\text{IHD}}\sum _{i=1}^{N_{\text{bus}}}\Delta {\text{IHD}}_i^2.\end{array}$$

(iii)

Fitness function for improving load bus voltage:

$$\begin{array}{}\text{(42)}& {\text{FF}}_3=F_3+{\mu }_V\sum _{i=1}^{N_{\text{bus}}}\Delta V_i^2+{\mu }_I\sum _{i=1}^{N_{\text{br}}}\Delta I_n^2+{\mu }_{\text{THD}}\sum _{i=1}^{N_{\text{bus}}}\Delta {\text{THD}}_i^2+{\mu }_{\text{IHD}}\sum _{i=1}^{N_{\text{bus}}}\Delta {\text{IHD}}_i^2.\end{array}$$

(iv)

Fitness function for harmonic distortion minimization:

$$\begin{array}{}\text{(43)}& {\text{FF}}_4=F_4+{\mu }_V\sum _{i=1}^{N_{\text{bus}}}\Delta V_i^2+{\mu }_I\sum _{i=1}^{N_{\text{br}}}\Delta I_n^2.\end{array}$$

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 N_Co new solutions for each group while the second generation applies Algorithm 1 for producing one new solution for each group. After producing new solutions (i.e., new site and new rated power for each PVDGU), the site and the rated power are checked and adjusted by the function of the two following equations:$$\begin{array}{c}\text{(44)}& {\text{Site}}_j=\left\{\begin{array}{cc}2,& \text{if\hspace{0.17em}}{\text{Site}}_j<2,\\ N_{\text{bus}},& \text{if\hspace{0.17em}}{\text{Site}}_j>N_{\text{bus}},\\ {\text{Site}}_j,& \text{else,}\end{array}\right.{\text{RP}}_j=\left\{\begin{array}{cc}{\text{RP}}^{\min },& \text{if\hspace{0.17em}}{\text{RP}}_j<{\text{RP}}^{\min },\\ {\text{RP}}^{\max },& \text{if\hspace{0.17em}}{\text{RP}}_j>{\text{RP}}^{\max },\\ {\text{RP}}_j,& \text{else}.\end{array}\right.\end{array}$$

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

5.6. 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 Figure 1.

[figure omitted; refer to PDF]

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 omitted; refer to PDF]

[figure omitted; refer to PDF]

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 loss

 

Case 2: Minimizing the rated power of all PVDGUs

 

Case 3: Improving the voltage profile of all loads

 

Case 4: Minimizing the harmonic distortion

In addition to ICOA, eight other methods consisting of BBO, GA, PSO, SFO, SSA, COA, ICOA1, and ICOA2 are also implemented for dealing with the four cases above for the two power systems. For each case, each method is run 50 trial times and the fitness functions of the trial runs are recorded for reporting the best fitness, average fitness, and maximum fitness in aim to support evaluation of performance. MATLAB program language and a personal computer with a processor of 2.0 GHz and 2.0 GB RAM are the main tools for executing simulation results. Before running these methods, their control parameters are analyzed and selected based on experiments as well as referred from previous studies such as [1, 10, 11] for GA, [17, 18, 44] for PSO, [23, 52] for BBO, [53] for SFO, [54] for SSA, and [37] for COA.

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 (m_max) is set to five values from 0.1 to 0.5; habitat modification probability (P^mod) is equal to 1; and combining ratio of the selected habitat from the previous generation and current generation for next generation for each iteration is set to the range 20%–80% [26].

(2)

For running GA, the optimal parameters are used as follows: the crossover probability (Cr) is set to three 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 factors c₁ and c₂ are set to 2 [44].

(4)

To implement SFO, the day (d) and the sun (s) are, respectively, fixed at 100 and 1 while the pollination (p) is ranged from 0.6 to 1 with a step of 0.2 [53].

(5)

For SSA, c₁ is a function of $2e^{-{\left(4\text{IT}/{\text{IT}}_{\max }\right)}^2}$; c₂ and c₃ are randomly produced numbers in the interval of [0, 1] [54].

(6)

In order to run COA, ICOA1, ICOA2, and ICOA methods, the number of coyotes in each group N_Co and the number of groups are set to 4. According to the settings, the population N_pop of the four methods is 20. In addition, ICOA2 and ICOA have more than one advanced control parameter P_a and the parameter is set to five values from 0.2 to 1.0 with a step of 0.2. The impact of P_a on simulation results will be surveyed and discussed.

On the other hand, the applied methods have two common control parameters: they are population N_pop and the maximum number of iteration IT_max. N_pop of PSO, GA, BBO, SSA, and SFO is also set to that of COA methods, i.e., 20, while IT_max has the same values for all methods. IT_max is 75 for the IEEE 33-bus distribution network and 100 for the IEEE 69-bus distribution network.

6.1. Survey on Results Obtained by Setting Different Values of P_a for ICOA Method

For finding the most suitable P_a value, ICOA is in turn run 50 times corresponding to each P_a value by setting P_a to 0.2, 0.4, 0.6, 0.8, and 1.0. 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 P_a = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0 is 0.2860, 0.2836, 0.2978, 0.2981, 0.2900, and 0.2978, respectively. Through the results, it points out that ICOA method can reach the best solution only at P_a = 0.2. Additionally, the average values of fitness function also confirm the most effective impact of P_a = 0.2 on its average fitness is the smallest among five values. The best fitness can conclude P_a = 0.2 is the most suitable value for finding the best optimal solution while the best average fitness can confirm the most stable search ability of ICOA if setting P_a = 0.2. One more time, the significance of setting P_a = 0.2 can be adopted since the best fitness and the average fitness at P_a = 0.2 continue to be the smallest values among five values for the second system shown in Table 3. The best fitness and the best average are, respectively, 0.2425 and 0.2939. Consequently, we accept the best P_a value is 0.2 and it is used for other remaining cases with intent to shorten simulation time.

Table 2

The obtained results by ICOA for Case 2 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.4448

Table 3

The obtained results by ICOA for Case 2 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 (N_betterrun) that the proposed ICOA method can reach was compared to the best run of compared methods.

(4)

Improvement level (IL) of the best fitness, average fitness, and worst fitness from the proposed method over other compared methods is exactly computed. The three factors are obtained by using the following equations:

$$\begin{array}{}\text{(45)}& \text{IL\hspace{0.17em}of\hspace{0.17em}the\hspace{0.17em}best\hspace{0.17em}fitness}\left(\%\right)=\frac{\text{the\hspace{0.17em}best\hspace{0.17em}fitness\hspace{0.17em}of\hspace{0.17em}the\hspace{0.17em}compared\hspace{0.17em}method}-\text{the\hspace{0.17em}best\hspace{0.17em}fitness\hspace{0.17em}of\hspace{0.17em}ICOA}}{\text{the\hspace{0.17em}best\hspace{0.17em}fitness\hspace{0.17em}of\hspace{0.17em}the\hspace{0.17em}compared\hspace{0.17em}method}}\times 100,\text{(46)}& \text{IL\hspace{0.17em}of\hspace{0.17em}average\hspace{0.17em}fitness}\left(\%\right)=\frac{\text{average\hspace{0.17em}fitness\hspace{0.17em}of\hspace{0.17em}the\hspace{0.17em}compared\hspace{0.17em}method}-\text{average\hspace{0.17em}fitness\hspace{0.17em}of\hspace{0.17em}ICOA}}{\text{average\hspace{0.17em}fitness\hspace{0.17em}of\hspace{0.17em}the\hspace{0.17em}compared\hspace{0.17em}method}}\times 100,\text{(47)}& \text{IL\hspace{0.17em}of\hspace{0.17em}the\hspace{0.17em}worst\hspace{0.17em}fitnes}s\left(\%\right)=\frac{\text{the\hspace{0.17em}worst\hspace{0.17em}fitness\hspace{0.17em}of\hspace{0.17em}the\hspace{0.17em}compared\hspace{0.17em}method}-\text{the\hspace{0.17em}worst\hspace{0.17em}fitness\hspace{0.17em}of\hspace{0.17em}ICOA}}{\text{the\hspace{0.17em}worst\hspace{0.17em}fitness\hspace{0.17em}of\hspace{0.17em}the\hspace{0.17em}compared\hspace{0.17em}method}}\times 100.\end{array}$$

(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 FF₁ 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. Figure 4 showing fifty fitness values of the four COA methods can confirm the superiority of ICOA. Because the ICOA has many searched fitness values that are significantly lower than the COA.

As compared to BBO, GA, PSO, SFO, and SSA, the proposed ICOA shows its potential search ability since its best fitness and average fitness are the lowest values among six methods. 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%. Figures 4 and 5 are plotted for illustrating 50 fitness function values of 50 runs obtained by these methods and ICOA. As counting from the figure, ICOA has approximately twenty points with the same fitness functions and below all points of other methods. Furthermore, the oscillations of ICOA are very narrow and the deviation between the highest point and the lowest point is small whereas the oscillations and the deviation of other methods are much higher, especially those from GA. If considering the number of runs between ICOA and other ones, ICOA has 13 runs with better fitness than the best fitness of COA, 0.3462, while the number is much higher as compared to other ones. In fact, it is 15, 45, 27, 38, and 27 runs with better fitness than the best solution of BBO, GA, PSO, SFO, and SSA, respectively. 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 omitted; refer to PDF]

[figure omitted; refer to PDF]

Table 5

The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case 1.

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.017

Case 2.

: Minimizing PVDG penetration level in the system

In the second case, the sum of the rated power of all PVDGUs is used to be an objective function whereas load voltage profile, branch currents, and total harmonic distortion and individual harmonic distortion are taken into consideration. Fitness function FF₂ 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 Table 6. For clear observation, 50 runs from all methods are separated into two figures in which Figure 6 shows all runs from four COA methods and Figure 7 shows all runs from GA, BBO, PSO, SSA, SFO, and ICOA. For comparison among four COA methods, reported values are good evidences for demonstrating the effectiveness of ICOA over three remaining ones since its improvement levels are much higher than those in Case 1. ICOA can improve the best fitness, average fitness, and the worst fitness up to 5.718%, 6.138%, and 7.621%. ICOA1 continues to be the second best method, but ICOA2 is still the worst method. Clearly, the first proposed improvement still has better contribution to the performance of ICOA. 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 Case 1 when BBO is still the second best method and GA is still the worst method for average and the worst fitness comparison. SFO is the worst method for the best fitness comparison. 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 omitted; refer to PDF]

[figure omitted; refer to PDF]

Table 7

The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case 2.

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.299

Case 3.

: Improving the voltage profile

In the part, voltage of all loads is considered to be improved by using fitness function FF₃ in equation (42) as a comparison criterion for evaluating applied methods. In the fitness function, objective function F₃ in equation (7) together with branch currents, voltage limits, and total harmonic distortion and individual harmonic distortion is taken into consideration.

Results for comparison are reported in Table 8. 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 Table 9.

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 omitted; refer to PDF]

Table 9

The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case 3.

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.132

Case 4.

: Reducing the harmonic distortions

In the part, THD and IHD of all buses are improved by using fitness function FF₄ in equation (43) to compare the quality of obtained solutions from implemented methods. In the fitness function, objective function F₄ 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, ω₁ and ω₂ are set to 0.5 and 0.5, respectively.

Results reported in Table 10 show that the best optimal solution of the proposed method has approximately equal fitness or not much better fitness than that of other ones since the improvement level is approximately 0% to 4.066%. ICOA1 and ICOA have the same best solution, but all solutions on average of the proposed ICOA are better than those from ICOA1 through average fitness improvement level. In addition, the fluctuations of the proposed method are also smaller. 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 Table 11.

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 omitted; refer to PDF]

Table 11

The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 33-bus distribution network in Case 4.

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 omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

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 Case 1, 5.711% and 11.674% for Case 2, 10.329% and 11.531% for Case 3, and 0.309% and 1.178% for Case 4. As compared to other remaining methods excluding COA methods, the outstanding performance of ICOA is clearer since the compared values are much higher. For example, the improvement level of the best fitness is from 0.264% to 3.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 Tables 16–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 Case 1.

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.765

Table 17

The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case 2.

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.643

Table 18

The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case 3.

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.781

Table 19

The best solutions and the total capacity of PVDGUs found by nine implemented methods for the IEEE 69-bus distribution network in Case 4.

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 solution ${\text{Co}}_{\text{mid},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 including N_Co, N_g, and IT_Max. In this method, the population size (N_pop) and the number of new solutions for each iteration are equal to (N_Co $\times$ N_g) and (N_Co $\times$ N_g + N_g), respectively. During the simulation process, we surveyed and selected two parameters as N_Co and N_g to find the most appropriate values for two considered systems. The values of N_Co and N_g are adjusted to find the most suitable values. As observing obtained results, we have found that 4 is suitable for both N_Co and N_g whereas 75 and 100 are, respectively, reasonable for IT_Max for the first system and the second system. 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 ones

Connecting appropriate DGPVUs in the distribution system can bring both economic and technical benefits. As compared to other methods, the proposed method can

(i)

Reduce the power loss by 9.002 kW and 2.547 kW and minimize the penetration level of PVDGUs by 0.3398 MW and 0.4311 MW for the first and second systems, respectively

(ii)

Improve voltage and reduce harmonic distortion by 0.0027 pu and 0.289% for the first system and by 0.019 pu and 0.190% for the second system, respectively

Consequently, it can lead to a suggestion that the proposed method should be applied for finding the optimal location and rated power of PVDGUs in distribution power systems for the purposes of economic and technical issues.

Glossary

Nomenclature

${\text{Co}}_g^{\text{new}}$:The new solution in the $g^{\text{th}}$ pack

${\text{Co}}_{\text{best},g}$ and ${\text{Co}}_{\text{worst},g}$:The best solution and the worst solution in the $g^{\text{th}}$ group

Co_best,1, Co_best,2, Co_best,3, and Co_best,4:The best solutions in all packs picked up randomly

Co_Gbest:The best solution in the population

Co_k:The k^th solution in the $g^{\text{th}}$ group

${\text{Co}}_{1,g}$ and ${\text{Co}}_{2,g}$:Two randomly picked solutions in the $g^{\text{th}}$ pack

$\Delta I_n$:Current difference at the n^th branch

$\Delta {\text{IHD}}_i$:Individual harmonic distortion difference at the i^th bus

$\Delta {\text{THD}}_i$:Total harmonic distortion difference at the i^th bus

$\Delta V_i$:Voltage difference at the i^th bus

$F_1$:Objective function of total power losses

$F_2$:Objective function of the PVDGU penetration level

$F_3$:Objective function of voltage profile index

$F_4$:Objective function of harmonic distortion

${\text{FF}}_{k,g}^{\text{new}}$ and ${\text{FF}}_{k,g}$:Fitness function of the k^th new solution and the k^th old solution in the $g^{\text{th}}$ group

${\text{FF}}_1$:Fitness function for total power loss minimization

${\text{FF}}_2$:Fitness function for rated power minimization

${\text{FF}}_3$:Fitness function for improving bus voltage

${\text{FF}}_4$:Fitness function for minimizing harmonic distortions

${\text{FF}}_{\text{best},g}$:Fitness function of the best solution in the $g^{\text{th}}$ group

${\text{FF}}_{i,g}$ and ${\text{FF}}_{j,g}$:Fitness function of the i^th and the j^th solutions in the $g^{\text{th}}$ group

${\text{FF}}_{\text{mean},g}$:Average fitness function of all solutions in the $g^{\text{th}}$ group

H:Maximum order frequency

$I_{n\text{DG}}^h$:Maximum current magnitude in the n^th branch at the h^th order frequency

${\text{IHD}}_{V,i}^h$:Voltage individual harmonic distortion of the h^th order harmonic at the i^th bus

$I_n$:Current magnitude in the n^th branch without PVDGUs

$I_{n\text{DG}}$:Current magnitude in the n^th branch with PVDGUs

$I_n^{\text{max}}$:Maximum current magnitude in the n^th branch

$\text{IT}$:Current iteration

${\text{IT}}_{\text{max}}$:Maximum of iteration

N_betterruns:Number of better runs of MCOA as compared to other methods

$N_{\text{br}}$:Number of branches

$N_{\text{bus}}$:Number of buses

N_Co:Number of coyotes in each group

N_DG:Number of PVDGUs

N_dv:Number of decision variables

N_g:Number of groups

N_load:Number of loads

N_pop:Number of coyotes in all groups

$P_{\text{Grid}}$:Total active power supplied by grid

$P_{\text{Load},i}$:Active power of the i^th load

$P_{\text{Loss},n}$:Active power losses at the n^th branch

$R_n$:Resistance of the n^th branch

THD_V,i:Voltage harmonic distortion at the i^th bus

THD^max:Maximum total harmonic distortion

THD_Vaver:Total harmonic distortion average value of the system

TPL:Total power losses of all branches without PVDGUs

TPL_DG:Total power losses of all branches with PVDGUs

${\mu }_I$:Penalty factor for violated current

${\mu }_V$:Penalty factor for violated voltage

${\mu }_{\text{THD}}$:Penalty factor for violated total harmonic distortion

${\mu }_{\text{IHD}}$:Penalty factor for violated individual harmonic distortion

$V_i^h$:Voltage magnitude of the i^th bus at the h^th order frequency

V^min and V^max:Lower and upper limitations of load bus voltage magnitude

VPI:Voltage profile index of the system.

References

[1] M. Kashyap, A. Mittal, S. Kansal, "Optimal placement of distributed generation using genetic algorithm approach," Proceeding of the Second International Conference on Microelectronics, Computing & Communication Systems, pp. 587-597, .

[2] A. Z. Kamaruzzaman, A. Mohamed, H. Shareef, "Effefct of grid-connected photovoltaic systems on static and dynamic voltage stability with analysis techniques—a review," Universiti Kebangsaan Malaysia., vol. 1 no. 6, pp. 134-138, DOI: 10.15199/48.2015.06.27, 2015.

[3] M. A. Al-Sobounchi, J. Gow, M. Al-Akaidi, "Optimal sizing and location of large PV plants on radial distribution feeders for minimum line losses," Proceedings of the 4th International Conference on Electric Power and Energy Conversion Systems (EPECS),DOI: 10.1109/epecs.2015.7368528, .

[4] E. C. Ashigwuike, S. A. Benson, "Optimal location and sizing of distributed generation in distribution network using adaptive neuro-fuzzy logic technique," European Journal of Engineering Research and Science, vol. 4 no. 4, pp. 83-89, DOI: 10.24018/ejers.2019.4.4.1237, 2019.

[5] M. Al-Maghalseh, S. Odeh, A. Saleh, "Optimal sizing and allocation of DGs for real power loss reduction and voltage profile improvement in radial L.V networks," Proceedings of the 14th International Conference on Smart Cities: Improving Quality of Life Using ICT & IoT (HONET-ICT), pp. 21-25, DOI: 10.1109/HONET.2017.8102213, .

[6] R. S. Ghatak, P. Acharjee, "Optimal allocation of DG using exponential PSO with reduced search space," Proceedings of the Second International Conference on Computational Intelligence & Communication Technology (CICT), pp. 489-494, DOI: 10.1109/CICT.2016.103, .

[7] A. Daud, A. F. A. Kadir, K. C. Gan, R. A. Abdullah, "Optimal allocation and sizing of photovoltaic-based distributed generation for voltage dip improvement," Proceedings of the IEEE International Conference on Power and Energy (PECon), pp. 89-94, DOI: 10.1109/PECON.2016.7951539, .

[8] A. R. Gupta, "Effect of optimal allocation of multiple DG and D-STATCOM in radial distribution system for minimizing losses and THD," Proceedings of the 7th International Symposium on embedded Computing and System Design,DOI: 10.1109/ISED.2017.8303936, .

[9] A. M. Eltamaly, Y. S. Mohamed, A. M. El-sayed, A. N. A. Elghaffar, "Impact of distributed generation (DG) on the distribution system network," International Journal of Engineering, vol. 17, pp. 165-170, 2019.

[10] R. Sulistyowati, C. D. Riawan, M. Ashari, "PV farm placement and sizing using GA for area development plan of distribution network," Proceedings of the International Seminar on Intelligent Technology and Its Application, pp. 509-514, DOI: 10.1109/ISITIA.2016.7828712, .

[11] I. I. Ali, Y. S. Saeed, "Optimal reconfiguration and distributed generation placement in Baghdad distribution sector," Engineering and Technology Journal, vol. 36 no. 4, pp. 333-343, DOI: 10.30684/etj.36.3a.13, 2018.

[12] S. K. Injeti, N. Prema Kumar, "A novel approach to identify optimal access point and capacity of multiple DGs in a small, medium and large scale radial distribution systems," International Journal of Electrical Power & Energy Systems, vol. 45 no. 1, pp. 142-151, DOI: 10.1016/j.ijepes.2012.08.043, 2013.

[13] S. A. Heydari, T. Heydarzadeh, N. M. Tabatabaei, "A combined approach for loss reduction and voltage profile improvement in distribution systems," International Journal on Technical and Physical Problems of Engineering, vol. 26 no. 8, pp. 30-35, 2016.

[14] R. Rajaram, K. Sathish Kumar, N. Rajasekar, "Power system reconfiguration in a radial distribution network for reducing losses and to improve voltage profile using modified plant growth simulation algorithm with distributed generation (DG)," Energy Reports, vol. 1, pp. 116-122, DOI: 10.1016/j.egyr.2015.03.002, 2015.

[15] A. B. Sarfaraz, S. Singh, "Optimal allocation and sizing of distributed generation for power loss reduction," Proceedings of the International Conference & Workshop on Electronics and Telecommunication Engineering, pp. 15-20, DOI: 10.1049/cp.2016.1116, .

[16] T. Yuvaraj, K. R. Devabalaji, K. Ravi, "Optimal allocation of DG in the radial distribution network using bat optimization algorithm," Advances in Power Systems and Energy Management, vol. 436, pp. 563-569, DOI: 10.1007/978-981-10-4394-9_55, 2018.

[17] W. Phuangpornpitak, K. Bhumkittipich, "Principle optimal placement and sizing of single distributed generation for power loss reduction using particle swarm optimization," Research Journal of Applied Sciences, Engineering and Technology, vol. 7 no. 6, pp. 1211-1216, DOI: 10.19026/rjaset.7.382, 2014.

[18] A. Ameli, S. Bahrami, F. Khazaeli, M.-R. Haghifam, "A multiobjective particle swarm optimization for sizing and placement of DGs from DG owner’s and distribution company’s viewpoints," IEEE Transactions on Power Delivery, vol. 29 no. 4, pp. 1831-1840, DOI: 10.1109/tpwrd.2014.2300845, 2014.

[19] S. Mahajan, S. Vadhera, "Optimal location and sizing of distributed generation unit using human opinion dynamics optimization technique," Distributed Generation & Alternative Energy Journal, vol. 33 no. 2, pp. 38-57, DOI: 10.1080/21563306.2018.12002410, 2018.

[20] M. Tolba, H. Rezk, A. A. Z. Diab, M. Al-Dhaifallah, "A novel robust methodology bases salp swarm algorithm for allocation and capacity of renewable distributed generators on distribution grids," Energies, vol. 11 no. 10,DOI: 10.3390/en11102556, 2018.

[21] M. M. Othman, W. El-Khattam, Y. G. Hegazy, A. Y. Abdelaziz, "Optimal placement and sizing of distributed generators in unbalanced distribution systems using supervised big bang-big Crunch method," IEEE Transactions on Power Systems, vol. 30 no. 2, pp. 911-919, DOI: 10.1109/tpwrs.2014.2331364, 2015.

[22] M. Dixit, P. Kundu, R. H. Jariwala, "Optimal placement and sizing of DG in distribution system using artificial bee Colony algorithm," Proceedings of the IEEE 6th International Conference on Power Systems (ICPS),DOI: 10.1109/ICPES.2016.7584010, .

[23] A. S. Vizhiy, K. R. Santhi, "Biogeography based optimal placement of distributed generation units in distribution networks: optimal placement of distributed generation units," Proceedings of the International Conference on Electrical, Electronics, and Optimization Techniques (ICEEOT), pp. 2245-2250, DOI: 10.1109/ICEEOT.2016.7755092, .

[24] F. Umar, M. Ashari, O. Penangsang, "Optimal location, size and type of DGs to reduce power losses and voltage deviation considering THD in radial unbalanced distribution systems," Proceedings of the International Seminar on Intelligent Technology and Its Applications (ISITIA), pp. 557-582, DOI: 10.1109/ISITIA.2016.7828724, .

[25] A. F. A. Kadir, H. Shareef, C. Z. M. Wanik, A. A. Ibrahim, Optimal Sizing and Placement of Distributed Generation in Distribution System Considering Losses and THDv Using Gravitational Search Algorithm, 2013.

[26] M. Q. Duong, T. D. Pham, T. T. Nguyen, A. T. Doan, H. V. Tran, "Determination of optimal location and sizing of solar photovoltaic distribution generation units in radial distribution systems," Energies, vol. 12 no. 1, pp. 174-198, DOI: 10.3390/en12010174, 2019.

[27] G. Srinivasan, S. Visalakshi, "Application of AGPSO for power loss minimization in radial distribution network via DG units, capacitors and NR," Energy Procedia, vol. 117, pp. 190-200, DOI: 10.1016/j.egypro.2017.05.122, 2017.

[28] B. Wojszcyk, M. Brandao, "High penetration of distributed generation and its impact on electric grid performance—utility perspective," Proceedings of the IEEE PES Innovative Smart Grid Technologies,DOI: 10.1109/ISGT-Asia.2011.6167146, .

[29] S. J. Rudresha, G. S. Ankaliki, "Voltage stability improvement with integration of DG into the distribution system," International Journal of Emerging Trends & Technology in Computer Science (IJETTCS), vol. 6 no. 5, pp. 196-200, 2017.

[30] A. A. Abdelsalam, A. A. Zidan, E. F. El-Saadany, "Optimal DG allocation in radial distribution systems with high penetration of non-linear loads," Electric Power Components and Systems, vol. 43 no. 13, pp. 1487-1497, DOI: 10.1080/15325008.2015.1043601, 2015.

[31] J.-H. Teng, C.-Y. Chang, "Backward/forward sweep-based harmonic analysis method for distribution systems," IEEE Transactions on Power Delivery, vol. 22 no. 3, pp. 1665-1672, DOI: 10.1109/tpwrd.2007.899523, 2007.

[32] J. Pierezan, G. Maidl, E. M. Yamao, L. D. S. Coelho, C. V. Mariani, "Cultural coyote optimization algorithm applied to a heavy duty gas turbines operation," Energy Conversion and Management, vol. 199,DOI: 10.1016/j.enconman.2019.111932, 2019.

[33] H. M. Qais, M. H. Hasanien, S. Alghuwainem, S. A. Nouh, "Coyote optimization algorithm for parameters extraction of three-diode photovoltaic models of photovoltaic modules," Energy, vol. 187,DOI: 10.1016/j.energy.2019.116001, 2019.

[34] U. Güvenç, E. Kaymaz, "Economic dispatch integrated wind power using coyote optimization algorithm," Proceedings of the 7th International Istanbul Smart Grids and Cities Congress and Fair, pp. 179-183, DOI: 10.1109/SGCF.2019.8782354, .

[35] T. T. Nguyen, D. N. Vo, H. Van Tran, L. Van Dai, "Optimal dispatch of reactive power using modified stochastic fractal search algorithm," Complexity, vol. 2019,DOI: 10.1155/2019/4670820, 2019.

[36] A. Betka, N. Terki, A. Toumi, H. Dahmani, "Grey Wolf optimizer-based learning automata for solving block matching problem," Signal, Image and Video Processing,DOI: 10.1007/s11760-019-01554-w, 2019.

[37] J. Pierezan, L. D. S. Coelho, "Coyote optimization algorithm: a new metaheuristic for global optimization problems," Proceedings of the IEEE Congress on Evolutionary Computation (CEC),DOI: 10.1109/CEC.2018.8477769, .

[38] A. B. Allah, L. Djamel, "Control of power and voltage of solar grid connected," International Journal of Electrical and Computer Engineering (IJECE), vol. 6 no. 1, pp. 26-33, 2016.

[39] H. Hu, Q. Shi, Z. He, J. He, S. Gao, "Potential harmonic resonance impacts of PV inverter filters on distribution systems," IEEE Transactions on Sustainable Energy, vol. 6 no. 1, pp. 151-161, DOI: 10.1109/tste.2014.2352931, 2015.

[40] E. George, Control System Design Guide, 2012.

[41] M. Suresh, B. E. J., "Optimal placement of DG units for loss reduction in distribution systems using one rank cuckoo search algorithm," International Journal of Grid and Distributed Computing, vol. 11 no. 1, pp. 37-44, DOI: 10.14257/ijgdc.2018.11.1.04, 2018.

[42] A. S. Hosseini, M. Karami, K. S. S. Madahi, F. Razavi, A. A. Ghadimi, "Determining the optimal capacity and place of DGs in distribution systems," Applied Mechanics and Materials, vol. 110–116, pp. 5195-5199, DOI: 10.4028/www.scientific.net/amm.110-116.5195, 2012.

[43] H. Iyer, S. Ray, R. Ramakumar, "Voltage profile improvement with distributed generation," Proceedings of the IEEE Power Engineering Society General Meeting,DOI: 10.1109/PES.2005.1489406, .

[44] M. R. Sasirajia, K. Muthulakshmi, S. V. Kumar, T. Abinaya, "PSO based optimal distributed generation placement and capacity by considering harmonic limits," Tehnički Vjesnik, vol. 24, pp. 391-398, 2017.

[45] A. Wazir, N. Arbad, "Analysis and optimization of IEEE 33 bus radial distributed system using optimization algorithm," (JETAE) Journal of Emerging Trends in Applied Engineering, vol. 1 no. 2, pp. 17-21, 2016.

[46] L. Gong, W. Cao, K. Liu, J. Zhao, X. Li, "Spatial and temporal optimization strategy for plug-in electric vehicle charging to mitigate impact on distribution network," Energies, vol. 11 no. 6,DOI: 10.3390/en11061373, 2018.

[47] H. Xing, X. Sun, "Distributed generation locating and sizing in active distribution network considering network reconfiguration," IEEE Access, vol. 5, pp. 14768-14774, DOI: 10.1109/access.2017.2732353, 2017.

[48] L. Manuel, B. Josep, M. Juan, C. Montserrat, G. Eulia, "Modelling harmonics drawn by nonlinear loads," Proceedings of the 2015 9th International Conference on Compatibility and Power Electronics (CPE), pp. 93-97, DOI: 10.1109/cpe.2015.7231055, .

[49] D. R. Prabha, T. Jayabarathi, "Determining the optimal location and sizing of distributed generation unit using plant growth simulation algorithm in a radial distribution network," WSEAS Transactions on Systems, vol. 13, pp. 543-550, 2014.

[50] EPRI, Harmonic Spectrum, 2014. https://github.com/tshort/OpenDSS/blob/master/Test/indmachtest/Spectrum.DSS

[51] V. Vita, "Development of a decision-making algorithm for the optimum size and placement of distributed generation units in distribution networks," Energies, vol. 10 no. 9, pp. 1433-1446, DOI: 10.3390/en10091433, 2017.

[52] W. Guo, M. Chen, L. Wang, Y. Mao, Q. Wu, "A survey of biogeography-based optimization," Neural Computing and Applications, vol. 28 no. 8, pp. 1909-1926, DOI: 10.1007/s00521-016-2179-x, 2017.

[53] G. F. Gomes, S. S. da Cunha, A. C. Ancelotti, "A sunflower optimization (SFO) algorithm applied to damage identification on laminated composite plates," Engineering with Computers, vol. 35 no. 2, pp. 619-626, DOI: 10.1007/s00366-018-0620-8, 2018.

[54] S. Mirjialili, H. A. Gandomi, Z. S. Mirjalili, S. Saremi, H. Faris, M. S. Mirjalili, "Salp swarm algorithm: a bio-inspired optimizer for engineering design problems," Advances in Engineering Software, vol. 114, pp. 163-191, 2017.