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.

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]:

By using the same method, total power losses in all branches before installing PVDGUs, TPL, is determined by means of the following equation [41]:

It shows that if PVDGUs are properly installed, the losses of power system after installation are always smaller than before installation. 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:

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]:

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 aswhere 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

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:

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]

Individual harmonic distortion at each bus can be defined by [44]

According to IEEE Std. 519 [44], and should not be greater than 5% and 3%, respectively.

The average of can be defined by [25]

Similarly, the average of can be defined by

Here, the objective function should be shown as follows:where and are, respectively, weighted factors with and . and are constrained by

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:

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:

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]:where and 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:where and 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:where and 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 group has two important factors consisting of social condition and quality 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:where Co^min and Co^max are, respectively, lower bound solution and upper bound solution obtained bywhere and 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:where is the objective function of the k^th solution in the group; PF is a penalty factor, which is selected by experience; is the penalty term for violating the m^th constraint of the k^th solution in the 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 .

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:In equation (25), and are two randomly picked solutions in the considered group while is the middle solution obtained from the group. Each decision variable in the middle solution is obtained by one out of two ways shown in the following equation:where var_x,mid1 and var_x,mid2 are the x^th center decision variables corresponding to two different cases of N_Co. If N_Co is an odd number, var_x,mid1 is taken from the x^th decision variable at the position [(N_Co + 1)/2]. On the contrary, var_x,mid2 is the average of two x^th decision variables at two positions, N_Co/2 and [(N_Co + 1)/2]. However, before determining the two variables in equation (23), each decision variable of all available solutions is sorted in descending order.

Step 4. Evaluate the effectiveness level of new solutions
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, , and new social condition corresponding to new solution, . The quality of the two social conditions is the fitness function of the two solutions, and . So, it should retain one social condition for each k^th coyote by using the following rules:

Step 6. Create one new social condition in each group
In the step, each group produces one new solution by using a randomization mechanism. The new solution and the technique are shown in the following models:where and are two x^th variables randomly picked from two available solutions; 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 .

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 for deciding if the worst solution should be replaced with the new one. The decision is taken over by the comparison below:

Step 8. Exchange solutions among N_g groups
To diversify solutions and avoid falling into local optimal zones with premature convergence, COA carries out solution exchange action among N_g groups. Two randomly selected groups will provide two randomly picked solutions and their position is then exchanged. The task is dependent on the randomization mechanism below:where ϕ is a random number within the range from 0 to 1. And solution exchange action will take place if the condition in equation (28) is true. Clearly, the number of coyotes in each group is very important for exchanging solutions. 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 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:

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:

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.

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 group are, respectively, represented by and . Each solution in each group is represented by the following equation:where N_DG is the number of installed PVDGUs in the distribution network.

In order to produce the initial population, lower bound solution Co^min and upper bound solution Co^max should be defined to make sure that randomly generated solutions can satisfy the limit constraints:where and are, respectively, bus 2 and bus N_bus and and 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 and . 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:

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 and load bus voltage 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:

As a result, the two parameters are checked and penalized by using the following formulas:

In addition to taking voltage and current into consideration, THD and IHD are also constrained for better power quality. THD and IHD are calculated using equations (8) and (9). Then, they are checked and penalized as follows:

5.3. Fitness Function

As described in Section 2, we focus on four single objectives consisting of total power loss minimization F₁, 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:(ii)Fitness function for rated power minimization of PVDGUs:(iii)Fitness function for improving load bus voltage:(iv)Fitness function for harmonic distortion minimization:

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:

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 1 

Flowchart of the computation process.

6. Simulation Results

In order to investigate the real performance of the proposed ICOA method in finding the optimal location and sizing of PVDGUs, IEEE 33-bus distribution power network and IEEE 69-bus distribution power network [26, 49] are used as two main studied tests. The IEEE 33-bus distribution power network shown in Figure 2 has a total load demand of 3715 kW and 2300 kVAR and operates at a nominal voltage of 12.66 kV. The IEEE 69-bus distribution power network shown in Figure 3 operates at a nominal voltage of 12.66 kV and supplies power to total load demand of 3802 kW and 2694 kVAR. The two systems are considering linear loads. So, in order to produce negative impacts of harmonics, we inject the harmonic currents with five harmonic spectrums simultaneously to the IEEE 33-bus distribution power network at 6 load buses 10, 15, 20, 24, 27, and 32 and to the IEEE 69-bus distribution power network at 10 load buses 10, 13, 18, 19, 22, 25, 34, 46, 56, and 65. The detailed information of the five harmonics is shown in Table 1 [50].

Figure 2 

IEEE 33-bus distribution network.

Figure 3 

IEEE 69-bus distribution network.

The capacity for each PVDGU can be defined as a variable varying from 0 to 2 MW [51]. Two limit points 0 and 2 MW are treated as lower and upper bounds of rated power for each PVDGU, respectively. Besides the capacity, the position that we can locate each PVDGU is from the smallest bus number (bus 2) to the highest bus number (bus 33 for the first system and bus 69 for the second system). Due to the total load demand as well as the system dimension, three PVDGUs are installed in the IEEE 33-bus distribution network whereas four PVDGUs are installed in the IEEE 69-bus distribution network. The optimal location and optimal capacity of PVDGUs in the two systems aim at reaching four single-objective functions considered in the four following cases. Case 1: Minimizing the total power 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 ; 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.

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:(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.

Figure 4 

The fitness values of four methods in 50 trial times for the first case.

Figure 5 

The fitness values of six methods in 50 trial times for the first case.

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.

Figure 6 

The fitness values of four methods in 50 trial times for the second case.

Figure 7 

The fitness values of six methods in 50 trial times for the second case.

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.

Figure 8 

The fitness values of seven methods in 50 trial times for the third case.

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.

Figure 9 

The fitness values of seven methods in 50 trial times for the fourth case.

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.