An Improved Equilibrium Optimizer for Optimal Placement of Distributed Generators in Distribution Systems considering Harmonic Distortion Limits

Tis paper proposes an improved equilibrium optimizer (IEO) for determining optimal location and efective size of distributed generation units (DGUs) in the distribution systems in order to minimize the total power loss on distribution branches, investment cost, and operation and maintenance cost. In a good obtained solution, limits of voltage, current, and harmonic fows are also seriously considered, exactly satisfying predetermined ranges. Especially, individual harmonic distortion (IHD) and total harmonic distortion (THD) of bus voltage must fall into IEEE Std. 519. Te proposed IEO is developed from the original equilibrium optimizer (EO), which was motivated by control volume mass balance models. Tis novel algorithm can efectively expand the search area and avoid the premature convergence to low-quality solution spaces. With the determined solutions from IEO, not only are the voltages well improved but also the harmonics are mitigated from the violated values down to the allowable values of IEEE Std. 519. Moreover, the total power loss is signifcantly reduced from 0.2110MW to 0.0815MW, 0.2245MW to 0.07197MW, and 0.3161MW to 0.1515MW for IEEE 33-bus, IEEE 69-bus, and IEEE 85-bus radial distribution systems, respectively. Not only that, the total cost of DGUs is also more economical and consumes only $3.4753 million, $3.2840 million, and $3.0593 million corresponding to the three systems for a 20-year planning period. Te performance of the proposed algorithm is compared to three other implemented methods consisting of artifcial bee colony (ABC) algorithm, salp swarm algorithm (SSA), and EO and eight previously published methods for the three test systems. Te comparisons of results imply that IEO is better than other methods in terms of performance, stability, and convergence characteristics.


Introduction
Nowadays, it is very common to integrate distributed generation sources into distribution systems due to the development of renewable energies. Properly determining the location and capacity of DGUs in the integrated systems can reach great achievement such as minimizing the losses on the lines, improving the voltage quality of the system, mitigating the harmonic distortions, and reducing the operation and maintenance costs [1]. On the contrary, if the installation of DGUs is not appropriately planned, it can cause undesirable efects such as increasing power losses, decreasing voltage quality, increasing operation costs and other costs, and so on [2]. Tus, optimal planning strategy of DGUs in the distribution systems is essential to maximize the economic and technical achievements.
Te achieved benefts derived from the installation of DGUs in distribution systems mainly depend on the placement locations and the rated power of DGUs in which optimization algorithms are utilized for determining the planning of DGUs in the integrated system. In [3][4][5], the authors used GA for determining the two factors of DGUs with intent to reduce losses on transmission sectors and boost the bus voltage to the best stable region in the distribution systems. GA is a simple optimization technique, and it is known as a local search method reaching inefective solutions for the optimization problem. To overcome the disadvantages of GA, other authors in [6][7][8] proposed hybrid methods which divide tasks of each phase in the considered problem to apply the most appropriate algorithm to each phase. Tose authors applied the hybrid methods between GA and other techniques such as fuzzy, PSO, and RNN to search more favorable solutions for the proper integration of DGUs in 33-bus, 51-bus, and 69-bus distribution systems. In addition to GA, PSO is also a commonly selected method for fnding the efective solution in solving the problem of determining placement and capacity of DGUs [9][10][11]. Tese authors optimally determined the penetration of multiple DGUs into distribution systems. As a result, the voltage profle was greatly improved and losses on branches were efectively cut. However, this algorithm has reached a low convergence rate and unguaranteed stability. Tus, many researchers have sought to improve the original PSO. Te studies [12,13] proposed improvements by integrating uniform distribution, exponential attenuation inertia weight, cubic spline interpolation function, and learning factor of enhanced control. Tese studies have signifcantly contributed to a high efciency of the improved algorithms. In addition to the abovementioned methods, ABC is also widely used for optimizing the installation of DGUs [14,15]. Te authors have proved the efects of diferent types of DGUs to the loss reduction target in the distribution system, and they have successfully found a suitable planning strategy for DGU's installation with low total losses.
Likewise, the studies in [16,17] also proposed efcient hybrid methods for optimal installation of DGUs such as VSA-CBGA and BPSO-SLFA, respectively. By properly connecting DGUs, the voltage as well as power loss has changed positively with guaranteed energy balance in the system. Additionally, another strong method is applied in solving the same problem, named BBO [18]. Te authors demonstrated the obtained multi-beneft for the grid connected DGU system. Tis study considered the multi-objective function of minimizing the electricity purchase cost from the distribution company and maximizing the profitability of the DGU's owners. Similarly, in [19,20], with the same applied optimization method, the authors examined harmonic distortions on 33 and 69-bus systems. Te results from simulation have also proved that harmonic distortions could be signifcantly mitigated to IEEE Std. 519 after integrating DGUs into the considered systems. Besides, with the goal of enhancing distribution system capacity without causing adverse efects, revised IEEE Std. 1547-2018 has introduced the smart functions of interfacing photovoltaic inverter [21], and recent studies have also suggested to connect inverters which have smart functions in controlling voltage, frequency, and power fow in the DGU integrated system [22,23]. In addition, a few other studies have analyzed the benefts for voltage regulation, loss reduction, and congestion mitigation through integrating photovoltaic inverter to generate reactive power [24,25]. In the DGU integrated system, the inverter can be used for injecting or absorbing reactive power to maximize the achieved economic and technical welfare [26]. For the optimization problem in generating the active and reactive powers, the power factor of DGUs needs to be considered and parameter should be inserted accordingly [27]. About the environmental aspect, several researchers have looked at reducing annual operation costs and pollutant gas emissions [28]. Tese authors have suggested the multi-objective evolutionary algorithm which is based on the Pareto optimality, called SPEA2, to solve the optimization problem of network reconfguration and connection of DGUs. Not only that, the fuzzy set theory is also used for selecting the best compromise solution among achieved Pareto solution set in the distribution system under the consideration of many electrical stability criteria and variation of loads. Similarly, another positive method has recently been created for fnding the optimal solutions, named SSA [29]. SSA is developed based on the swarming behaviors of salps in oceans. In the mathematical model of the algorithm, the salp population is divided into two groups, the leader and the followers. Te leader is the leading salp and stands in front of the chain. Other remaining salps are considered as the followers. Te leader determines the direction of movement for followers. By comparing the obtained results with some other metaheuristic algorithms, the study [29] showed the superior efectiveness of SSA in solving the single and multiobjective functions. Additionally, SFO is also recently developed for solving the optimization problems [30]. It was known as a nature-inspired optimization algorithm that was inspired by sunfowers' motion toward the sun to get the radiation. Sunfowers that are located near the sun are calmer because greater heat is received in this space. Conversely, fowers that are further away from the sun tend to move closer to the sun because they receive less heat. Tis cycle is repeated daily in the morning [31]. Based on this natural feature, the algorithm has been built and applied in fnding the optimal solutions efectively. Besides, a recently efective method, called EO, was also published in 2020 [32]. Tis algorithm is developed based on control volume mass balance models that are applied to predict equilibrium state and dynamic state. In the mathematical model, each particle with its concentration is responsible for fnding the feasible solution in the space, and they are considered as the search agents. Te concentration of these agents is updated with respect to the current best solution to reach the equilibrium state. EO's strengths are in exploration and local minima avoidance. Its efectiveness has been demonstrated by comparing it with many other existing powerful methods. EO is really a great method in solving optimization problems with high stability. However, the performance of the EO can be enhanced with an improvement in the update equation, and this has been demonstrated in this study.
In this paper, an improvement in the original algorithm (EO) is implemented, called improved equilibrium optimizer (IEO). Te suggested IEO can inherit the strengths of EO and efectively avoid disadvantages of EO to reach a greater performance than EO. Tus, IEO is employed to fnd 2 Complexity the optimal solution in determining the sitting and sizing of solar photovoltaic distributed generation units in three distribution systems of 33 buses, 69 buses, and 85 buses. To satisfy many aspects of the technical and economic problems in the grid integrated DGU system, this study has considered the multi-objective function under various constraints. Te main objective of the study is to reduce the power loss, enhance the voltage stability, and minimize the costs of DGUs while the branch current, bus voltage, and harmonic distortions in the systems operating nonlinear loads are kept within the allowable limits. In addition, to solve harmonic power fow, the backward/forward sweep technique (BW/ FWST) has been used and the detailed description is clearly presented in [33]. Moreover, the weighted sum method is applied to decide the compromise solution for the multiobjective function [31]. In this technique, each weighting coefcient is assigned to each single objective function and the value of the weighting coefcient depends on the importance of each component in the multi-objective function.
Overall, this study includes the following major contributions: (1) For a project of DGU installation in distribution system, two major factors that should be calculated to evaluate the efectiveness are power loss and costs of investment, operation, and maintenance for the installed DGUs. Terefore, this study proposes to consider the multi-objective function under constraints of voltage profle, branch current, and harmonic distortions in diferent distribution systems of IEEE 33-bus, IEEE 69-bus, and IEEE 85-bus radial distribution systems. Te two objectives included in the multi-objective function are as follows: Minimize power loss on distribution lines in systems.
Minimize the sum of investment cost, operation cost, and maintenance cost for DGUs, which are integrated in distribution systems during a 20-year project.
(2) When DGUs are installed at buses in distributions systems, four key parameters including current of branches, power loss on branches, voltage of loads, and harmonic distortions will vary negatively or positively due to the power supply from installed DGUs. Hence, this paper analyzes the four parameters for two cases, before and after the installation of DGUs. As a result, proper determination of the sitting and sizing for DGUs can reach the reduction of power loss, the enhancement of voltage quality, the reduction of the costs of DGUs, and the mitigation of harmonic distortions. (3) Te gained benefts from integrating DGUs in a distribution system greatly depend on the found optimal solutions by applied algorithms. In other words, the selection or development of an efcient and stable optimal method will contribute to the both economic and technical benefts. Terefore, this paper has developed a novel algorithm (called IEO)   based on modifcations on equations updating solutions in the original algorithm (EO). Tereby, IEO  has signifcant improvements in fnding feasible  solutions with high quality while EO is less efective  when testing them on IEEE 33-bus, IEEE 69-bus, and  IEEE 85-bus radial distribution systems. Te rest of the paper is organized as follows. Firstly, section of the problem formulation (Section 2) describes the objective functions and constraints of this research. Secondly, section of the proposed algorithm (Section 3) shows the structure of EO and improvements of IEO on mechanism of EO. Tirdly, section of application of IEO for the DGUs problem (Section 4) presents computation steps to apply IEO for the problem. Next, section of simulation results (Section 5) analyzes the obtained results from integrating DGUs into IEEE 33-bus, IEEE 69-bus, and 85-bus radial distribution systems. Besides, the discussion of the performance of the proposed method and compared methods is also covered in this section. Te shortcomings of the proposed method and research expansion are also discussed in Section 5. Finally, Section 6 presents the summary of results, contributions, and future work of the study.

Objective Function.
Total power loss reduction on conductors is one of the key factors that mainly afect the economic and technical problems of distribution systems. Tus, it is important to take into account the minimization of the loss. Besides, when investing in a new integrated distribution system, the investors are always interested in the optimal strategy which can help to decrease costs on investments, maintenance, and operation. Tus, in order to satisfy the important criteria as mentioned above, the multi-objective function is proposed for use and it includes total active power loss (TAPL) and total cost of DGUs (TCDG). (TAPL). In the distribution system, the total active power loss is increasing more and more due to the development of loads as well as the extension of grid. Terefore, determining the optimal strategy in connecting DGUs to the system is a great way for the power loss minimization. Te reduction of power loss mainly contributes to the efectiveness of the power grid in terms of technical as well as economic factors in saving energy and reducing operation costs. Hence, in this study, TAPL is considered as the frst target in the multi-target problem. When performing the integration of DGUs into the grid, the mathematical equation of TAPL is described as follows [12]:

Total Active Power Loss
where TAPL DGU is the total power loss after connecting DGUs to the system. However, the compromise solution in the multi-objective function is difcult to determine because the single Complexity objective values have large deviations. Terefore, it is necessary to convert the objectives to the same considered range. Specifcally, in the frst single objective, we suggest the normalization equation and it is presented as equation (2) for determining the best compromise solution [20].
where TAPL NoDGU is the total power loss of the original system without DGUs and it can be found by using the equation below [17]: If DGUs are successfully integrated into the distribution system, the value of TAPL DGU will be smaller than the value of TAPL NoDGU . In other words, the value of OF I will vary in the range from 0 to 1 and the small value of the objective is expected. Tis demonstrates the positive efect of connecting DGUs in the system in reducing the total power loss.

Total Cost of DGUs (TCDG).
In order to improve technical efciency and maximize profts during the operation in a long period, it is necessary to minimize the related costs. In this case, the sum of investment cost and the cost of operation and maintenance is considered and becomes the second target in the multi-objective function. Te initial investment cost of all DGUs (IC can be found by [34] In equation (5), MOC opt DGU,d ($) is considered for a 20year planning period (i.e., PP � 20) and α y is the cumulative present value factor and obtained by where dr is the discount rate and y is the considered year varying from 1 to 20 (years). Finally, the total cost of DGUs is the sum of initial investment cost and maintenance and operation cost as the following equation: Similar to the frst objective, the second objective (OF II ) must be also converted into the same range from 0 to 1 for evaluating the multi-objective function correctly. Tus, the second objective can be normalized as follows: Here, TCDG max DGU is the total cost of the maximum generated power of DGUs in the distribution system. Te value of TCDG max DGU can be found by using equations (9)-(11) [34]: where IC max DGU,d ($) is the maximum initial investment cost of all DGUs and MOC max DGU,d is the maximum total cost of maintenance and operation.
In the considered optimization problem, the optimal generation of each DGU is searched within the minimum and maximum limits of the output power, which are given. Terefore, OF II always varies in the range from 0 to 1. Te smaller OF II is, the more proftable the distribution system is.
As mentioned, this study considers the multi-objective function consisting of two objectives as the power loss and the costs of DGUs. To determine the best compromise solution for the multi-objective function, the weighted sum method is used [31]. Lastly, the multi-objective function is established as follows: where OF I and ω I are the frst normalized single target and the weighting factor associated to the target. OF II and ω II are the second normalized single target and the weighting factor associated to the target, respectively. In equation (12), the two weighting factors are constrained by [20] ω I + ω II � 1 and 0 ≤ ω I , ω II ≤ 1.
Te values of ω I and ω II are chosen dependent on the importance of the components (single objective) in the multi-objective function [34]. Tese values can be adjusted by user. Te larger the weighting factor is, the greater its signifcance in the multi-objective function is.

Power Balance Constraints.
Te power balance is a key factor in keeping the stability for system frequency. If the power generation is more or less than the total power consumption, the frequency of the system will be increased 4 Complexity or decreased, respectively [35]. Tus, it is necessary to take the power balance between power generation and power consumption sides into account. Te power generation side includes the supplied power from the grid (AP gr ) and the injected power by connected DGUs (AP DGU,d ); meanwhile, the total power consumption side is comprised of the consumed power by loads (AP ld,b ) and the power losses (AP ls,bh ). Tese two power sides should be balanced as follows [31]:

Voltage Limits.
Te voltage profle plays a major role in evaluating quality power of the operating systems. However, when we consider the penetration of DGUs as another power source in the distribution system, the voltage has variations and its range is expanded. Tus, the voltage should be considered to keep in the limits. According to IEC and European EN 50160 standards, the allowable voltage limits are from 0.90 pu to 1.10 pu. However, in 33-bus, 69bus, and 85-bus radial distribution systems, the best range for the bus voltage is maintained within the minimum limit (V min ) and the maximum limit (V max ) of 0.95 pu and 1.05 pu, respectively. Obviously, the limit can also satisfy the standard of IEC and European EN 50160, but it is more serious than the two organizers. Te bus voltage and voltage limits are considered at the fundamental frequency and constrained as follows [17]:

Total Voltage Harmonic Distortion and Individual
Voltage Harmonic Distortion Limits. Harmonic distortions are considered a prime concern in the electrical power system because it can cause overheating (for capacitors, generators, and transformers), disoperation (for electronics, relays, and switchgear), reduction of life time of the connected devices in the system, and so on. Harmonic distortions can be defned as the deformation of a waveform by integer multiples of a fundamental frequency. As shown in Figure 1, they are the distorted waveforms of the 12-pulse converter, AC voltage regulator, and fuorescent lighting [35]. Harmonic distortions are caused by nonlinear load's operations (saturated electric machines or rectifers), and they can be minimized by using passive and active flters. To evaluate the power quality of the systems, the two parameters including total harmonic distortion and individual harmonic distortion must satisfy allowable ranges in Table 1 [36,37].
Te two factors which indicate the harmonic noise level of the system are total voltage harmonic distortion (THD b (%)) and individual voltage harmonic distortion . Teir values can be obtained by voltage values at the fundamental frequency (V 1 b ) and other higher order frequencies (V h b ) as follows [20]: By applying the IEEE Std. 519 of harmonic limits as shown in Table 1, THD max (%) and IHD max (%) in equations (16) and (17) are set to 5% and 3%, respectively, because the voltage level at the point of common coupling in the considered systems is less than 69 kV [33].

DGU's Capacity Limits.
Before photovoltaic strategy planning, for each DGU, the lower and upper bounds of location (L min DGU and L max DGU ) as well as capacity (AP min DGU and AP max DGU ) should be predetermined and imposed on DGUs as shown in equations (18) and (19), respectively [4]. In other words, the selection of location (L DGU,d ) and rated power (AP DGU,d ) of the d th DGU is performed by applying the predetermined allowed regions. Moreover, the total generating capacity of all DGUs must not exceed 80% of the total load demand (AP ld ) as equation (20) for this study [12,20].
2.2.5. Branch Current Limits. When DGUs are connected to the distribution system, the current in the branches can be signifcantly increased depending on the location and power of DGUs. Tis can disrupt the original structure of the grid and cause damage to existing conductors and other electric components. Terefore, the branch current limit should be Complexity 5 considered as one of the most serious constraints in systems after integrated DGUs as follows [28]:

Equilibrium Optimization (EO).
A strong optimization algorithm based on the mass balance control models was published in 2020 called equilibrium optimizer (EO) [32]. In this algorithm, each particle along with its concentration works as a search agent. Tese search agents update their concentration randomly with the desire of fnding the best quality solution. In other words, each solution represents an internal concentration and the adjustment variables in the solutions are the concentration parameters of EO. Te quality of the solution entirely depends on the ftness function value under the constraints, and the use of the ftness function result is intended to evaluate the massconcentration balance within the adjusted volume. Overall, EO is an ideal algorithm for exploration, fnding highquality optimal solutions and avoiding local traps. Te algorithm is expressed as follows.

Initial Population Generation.
In EO, each solution F i in the population contains a set of adjustment variables and there must be a predetermined range for the variables. Te lower and upper bounds of the variables are called the minimum solution (F min ) and the maximum solution (F max ), respectively. During the application of EO, each solution F i is always kept within these limits. Te mentioned solutions and the boundary solutions are formulated as follows: Te generation of initial solution set or initial population is implemented by [32]

Calculation Process for Updating Concentration.
Quality of all solutions is calculated to rank the solutions in population in which solution with the smallest ftness value is the best and other solutions with worse ftness values are from the second-best solution to the worst solution. After obtaining the ftness values of solutions, the top four best solutions are collected and assigned to updated through a combination of three components. Te frst component (F s ) is randomly selected from the solution set including F Lbest1 , F Lbest2 , F Lbest3 , F Lbest4 , and F Lmean . Clearly, the selection of a solution in the good group can contribute to a positive orientation for reaching the next generation. For the second component, it is the diference between the selected solution (F s ) and the considered old solution (F i ). Tis component is regarded as a part of distance between the current candidate solution and the new found solution. Another component is also another part of distance between old and new solutions, but it is built without the experience of the current solution. Te three components are useful in fnding better solutions, which are not close to the current good solutions. Finally, the new solution (F new i ) can be calculated by applying the following equation [32].
In equation (25), two important terms, which mainly afect the quality of produced new solutions, are the exponential term (ET) and the generation rate (GR). ET and GR values can be found by solving complex equations.
For the exponential term (ET), it can be calculated by applying the following equation [32].
Its main purpose is built to create the suitable balance between exploration and exploitation. Practically, the turnover rate will be changed in the real control volume with time. So, l is added in the equation to represent random variation over time in the range of (0, 1). In equation (26), the constant (β 1 ) is selected to be 2 and the integer number (μ) is randomly taken between 0 and 1. Besides, the time coefcient (tc) is a value that changes with each iteration, and its value depends entirely on the maximum iteration number (It Max ) and the current iteration number (It). Te value of tc can be found by applying the following equation [32].
where β 2 is a constant and it is selected to be 1 for the sake of simplicity.
For the generation rate (GR), it is the most important item in the algorithm of EO to suggest efcient solutions through enhancing the exploitation phase and GR can be determined by using the following equation [32]. where In equation (28), Gq is considered as the control coeffcient of GR. It positively contributes to the process of creating new solutions. Its value depends on random generated variables in the range of (0, 1) such as θ 1 and θ 2 and generation probability (q). Obviously, once Gq takes on zero, the update equation (25) will no longer exist the exploitation rate term. Terefore, the value of q is important in determining the existence of this term and q should be chosen to be 0.5 to achieve a good balance [32].

New Concentration Correction and Update.
Each new concentration is always checked for determining the violation in the predetermined limits, and the found violation will be corrected according to the rule below: By calculating ftness function, the quality of each concentration is determined. Solutions with better quality (i.e., concentrations with better mass balance) are retained. Terefore, new concentration quality (Fit new i ) should be compared with the current concentration quality (Fit i ) for storing better ones. Te explanation can be formulated as follows:

Improved Equilibrium Optimizer Algorithm (IEO).
By substituting equation (28) into equation (25), new solutions can be obtained: Equation (33) has three terms including F s , in which the second term and third term are the updated jumping steps to make the new solution F new i diferent from the old solution F s . In this case, the scaling factor of the frst jumping step is ET. Similarly, the scaling factor of the second jumping step is Gq.ET.(1 − ET)/l and it is set to ET′ for the sake of simplicity. Te values of ET and ET′ are investigated to avoid the incorrect evaluation for ET and ET′ due to the missing values from randomization in the range of 0 to 1. Tus, in this case, to minimize the randomization that can afect the evaluation results, l is run from 0 to 1 with a step size of 0.0001, which is equivalent to 10,000 steps. In other words, 10,000 values of l will be produced for evaluation. ETand ET′ will be found with respect to l values for determining their operating regions. Te operating regions of ET and ET′ are plotted in Figures 2 and 3.
Trough Figure 2, the obtained results clearly indicate that the values of ET completely fuctuate in the range of [− 0.6, 0.6], and this implies that the scaling factor of the frst jumping step is only suitable for fnding local solutions. Besides, as presented in Figure 3, the limits of ET′ are in the range of [− 1.2, 0.9]. However, only 806/10,000 values (corresponding to 8.06%) are out of range [− 0.6 0.6] and 5,014/10,000 values (corresponding to 50.14%) are zero as counted. Terefore, most of the scaling factor values make the second jumping step fall into the local trap, and more than half the chance of existence of ET′ is completely eliminated. In other words, most of operating ranges of ET and ET′ are only concentrated in the narrow margins. So, the use of ET and ET′ leads to a difculty to expand the search space. Clearly, equation (33) restricts the expansion of search spaces by the impact of the two scaling factors and this makes the performance of the algorithm low. Hence, the original algorithm (EO) is only considered appropriate for fnding the optimal solutions in the local search space. In this study, equation (33) has been improved to extend the search area.
In the frst jumping step, the current solution (F i ) should be moved to a new location near the location of the best current solution (F Lbest1 ) to inherit the best current solution's experience at the previous generation. Tis will contribute to increase the probability of fnding better quality solutions. In summary, the position of the best current solution is chosen instead of random selection from a group that contains the top four best solutions and an average solution of them at the previous generation. Additionally, two randomly selected solutions (F r1 and F r2 ) by taking variables in the F s solution will bring many advantages in approaching a higher quality solution. As a result, (F r1 − F r2 ) is always greater than zero and (F r1 − F r2 ) is multiplied by a randomization vector (c) which is comprised of terms between 0 and 1. Te results of multiplication facilitate the expansion of the search space to avoid local traps. Generally, equation (33) has advantages in fnding solutions in the local space, and equation (34) is useful in expanding the search area. Terefore, the great combination of equations (33) and (34) will improve the performance of EO signifcantly. In this research, the choice of the update equation depends on the quality of the current solution (Fit i ) and the mean ftness of all solutions (Fit mean ). If the quality value of Fit i is better than the value of Fit mean , equation (33) is used for generating new solutions. For other cases, equation (34) is applied for new solution update.
Te implementation process of IEO for a typical problem is shown in Figure 4.
In equation (35), L DGU,d,i and AP DGU,d,i are the location and capacity of the d th DGU of the i th solution with the limitations as follows: Te two parameters of DGUs will be randomly generated by

Penalty Term for Violation.
After the i th new solution is generated, voltage, current, and the harmonic violations are checked as follows:

Fitness Function.
Te ftness function is a combination of objective function and penalty functions. It is used to evaluate the quality of each solution. In this paper, the ftness function of the i th solution is constructed as follows: 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 Step No.  Step No.

Correction for the Violated Variables.
When a new solution (i.e., the new location and the new capacity of DGUs) is produced, it has to be checked and adjusted throughout the search for the optimal solution as follows:

Implementation Process of the Whole IEO for DGU's
Problem. Te process of searching location and capacity of DGUs in the distribution systems by using IEO is briefy presented in Figure 5, and the implemented steps are as follows: Step 1: survey and select parameters of proposed method. Te survey is implemented to select the appropriate number of concentrations and the maximum number of iterations. Besides, the parameters of the proposed algorithm are also selected.
Step 2: generate the initial solutions at random. Te initial solutions of the placement and sizing of DGUs are randomly created within the allowable limits, which are predetermined by using equations (38) and (39). For each newly generated solution, all constraints are checked to calculate the penalty terms for bus voltage, branch current, total harmonic distortion, and individual harmonic distortion by applying equations (40)-(43), respectively. After that, the sum of these penalty terms and the objective function are determined to calculate the ftness function as equation (44). Te function value is also quality of each solution.
Step 3: start loop in the algorithm. Te frst iteration is started to perform the promising solution search iterative algorithm.
Step 4: determine the group of solutions with high quality.
Based on their ftness values in Step 2, the four most efective solutions are determined and assigned to F Lbest1 , F Lbest2 , F Lbest3 , F Lbest4 for the calculation of the arithmetic mean (F Lmean ). Te selection of good solutions contributes signifcantly in improving the quality of the next generation.
Step 5: calculate important terms for the new solution generation equation. Two important terms with a high infuence on the performance of the proposed algorithm are the exponential term (ET) and the generation rate (GR). ET and GR can be found by solving complex equations (26) and (28), respectively. Tey are important terms of the new solution generation equation.
Step 6: check the conditions for selecting the method of new solution generation. Compare the ftness value of each solution (Fit i ) with the mean ftness value (Fit mean ) to select the equation for creating the next generation (F new i ).
Step 7: create the new solution generation. Either equation (33) or equation (34) is employed to calculate the i th new solution. If Fit i ≥ Fit mean occurs, the i th new solution is found by utilizing equation (34). For another case, equation (33) is used.
Step 8: check the violation for correction and calculate the ftness value for each new solution. Each created new solution should be checked for limit violations and corrected as the rules in equations (45)- (46). Penalty and objective values of these corrected solutions are calculated by using equations (40)- (43). Ten, ftness value in equation (44) is obtained for each new solution.
Step 9: fnd the best current solution. Trough the evaluation of the ftness values, the best current solution (F Lbest1 ) and its ftness value (Fit Lbest1 ) are retained by applying the rules in equations (31) and (32).
Step 10: check the termination condition. Te condition for stopping the iteration is checked by comparing It and It Max . If It � It Max happens, stop implementing the proposed method and report the best results. For other cases, It is increased to (It + 1) and go back to Step 4.

Simulation Results
Tis section presents the applications of EO, ABC, SSA, and the proposed IEO for minimizing power loss on all distribution branches and costs of DGUs while satisfying all constraints including branch current, voltage profle, and harmonic distortions. In addition to the comparison with EO, ABC, and SSA, the proposed IEO is also compared to other published methods such as GA, PSO, GA/PSO, SA, HSA, BFOA, and BSOA via the simulation for three IEEE distribution systems with 33, 69, and 85 buses.
In this paper, each implemented method is independently run 50 times on a personal computer with a 1.8 GHz processor and 8.0 GB RAM in MATLAB environment. To get a fair evaluation for the implemented methods, the population is surveyed from 20 to 40 with step size of 10, and the suitable population is selected as 30. Besides, ItMax is set to 200, 230, and 250 iterations for 33, 69, and 85 buses systems, respectively. Moreover, the parameters for methods are also taken appropriately. For ABC algorithm, the colony size, which is the sum of employed bees and onlooker bees, is set to the population, and employed bees and onlooker bees have the same quantity [38]. For SSA, the coefcient of c 1 is the obtained result from function of 2.e − (4.It/It Max ) 2 , and the  Figure 5: Te implementation process of IEO for placing DGUs in radial distribution networks. Complexity coefcients of c 2 and c 3 are generated numbers randomly from 0 to 1 [29]. To run the proposed IEO and EO, parameters of α 1 , α 2 , and GP are set to 2, 1, and 0.5, respectively [32]. For calculating the costs of DGUs in the multi-objective function, investment and maintenance and operation costs are, respectively, 770 ($/kW) and 0.01 ($/kWh) [4] while the discount rate is selected to be 9% for a 20-year plan. Besides, this study applied the weighted sum method to decide the best compromise solution. Tus, the value of each weighting coefcient depends on the importance of each element in the multi-objective function [31]. Its value can be adjusted by user. Te larger value of a weighting coefcient, the greater its signifcance in the multi-objective function. In this study, the authors considered minimizing of power loss to be more important than the costs of DGUs in the distribution system. Tus, the total power loss receives signifcant weight (ω I ) of 0.8 and light weight of the cost (ω II ) of 0.2. Furthermore, all test systems are evaluated in the case of harmonic distortions. Hence, the fve harmonic fows are simultaneously injected to linear loads, and the detailed information is presented in Table 2 [20]. Figure 6. Te load and line data of the system are collected from [31]. Te basic system has total power loss of 0.2110 MW and total load demand of 3.715 MW and 2.300 MVar. In this case, three DGUs can be installed at three diferent buses from bus 2 to bus 33, and the capacity of each DGU can be chosen from 0.0 MW to 2.0 MW [34]. Furthermore, harmonic fows are also injected simultaneously into six buses 9, 15, 20, 24, 29, and 32.

Case 1: IEEE 33-Bus Radial Distribution System. Te confguration of the 33-bus distribution system is shown in
After 50 trial runs are implemented, the worst, average, and best ftness values are collected and clearly presented in Table 3. According to the presented results in Table 3, the best ftness value of the implemented methods such as ABC, SSA, EO, and IEO is 0.3866, 0.3837, 0.3832, and 0.3828, respectively. Te value of the proposed method (IEO) is lower than that of ABC, SSA, and EO. In other words, IEO can fnd a more optimal solution than ABC, SSA, and EO. Moreover, the mean and worst ftness values of IEO are smaller than those of others. For a clearer view of performance comparison, ftness values of 50 trial runs of these methods are sorted and plotted in Figure 7. All points on the curve of IEO are located at the lower positions than those of others, and the fuctuation of IEO is very small in the range of [0.3828, 0.3896]. Tus, IEO is the most stable method. Table 4 shows the objective values in the multi-objective function for comparison. Te frst objective and the corresponding power loss of IEO are 0.3863 and 0.0815 MW, which are equal to those of SSA and better than those of others. About the contribution to power loss, the proposed solution of IEO can drastically reduce losses in the system from 0.2110 MW to 0.0815 MW, while others can reduce the loss to 0.0816 MW (the second-best method, EO) and 0.1061 MW (the worst method, GA [39]). Te power loss on branches in the base system and the hybrid system with DGUs by using optimal solution of IEO is plotted in Figure 8. Te fgure indicates that the power loss on branches in the IEEE 33-bus system before and after placing DGUs has a relatively high deviation and approximately all branches of the system with the installation of DGUs have lower power losses than those in the system without DGUs. Especially, the loss reduction is very high for branches close to slack bus 1, including branches 1, 2, 3, 4, 5, 7, and 27. Other branches cannot reach the same high power loss reduction, but the loss reduction is still relatively high, for example, from branches 8 to 12 and from branches 25 to 30. Te loss reduction on other remaining branches is not signifcant, but there are still small values. Te very small reduction cannot be identifed in the fgure due to the constraint of fgure size. In general, the installation of DGUs can bring high benefts to the system in reducing power loss.
Te second objective and the corresponding cost of IEO are 0.3690 and $3.4753 million, which are smaller than those of others excluding BFOA [41], BSOA [42], and ABC. Clearly, the solution from IEO can save $1.2246 million and $0.0021 million in comparison to the worst method (GA [39]) and the original method (EO), respectively. So, it can lead to a conclusion that the proposed IEO is more efective than GA [39], PSO [39], GA-PSO [39], LSF-SA [40], SSA, and EO. For the comparison with BFOA [41], BSOA [42], and ABC, the proposed IEO cannot reach better cost, but it reaches much better power loss. Tis is a trade-of in solving multi-objective function problem. Tis result is due to the selection of weighting coefcients, and it cannot be avoided when comparing a proposed method to many other methods. However, the proposed IEO has reached better ftness function than the three methods; meanwhile, ftness function is also the multi-objective function and used to evaluate performance of methods. As a result, the proposed IEO is also more efective than the compared methods in solving the problem of DGU placement in the IEEE 33-bus distribution system.
In this research, harmonic distortions and voltage profle are considered as the two constraints of the multi-objective function. According to IEEE Std. 519, THD and IHD should not exceed 5% and 3%, respectively. However, the modifed system with nonlinear loads and without DGU connection has exceeded the limitations of both THD and IHD as shown in Figure 9. Specifcally, the largest values of THD and IHD at frequency orders are 7.2504% and 4.7037%, respectively. But the violation of harmonic distortions is solved after placing three DGUs by using IEO, and these factors are reduced to 4.6203% and 3.0%, respectively.
Similarly, the voltage profle of the modifed system before and after placing three DGUs is also plotted in Figure 10. In the system without DGUs, many buses have Optimal solutions of the applied system are reported in Table S1 in Supplementary Materials. Figure 11 shows the structure of the IEEE 69-bus radial distribution system, and data of the system can be taken from [20]. In addition, general information of the system is comprised of total power loss of 0.2245 MW and total loads of 3.8019 MW and 2.6941 MVar. In this system, three DGUs can be installed from buses 2 to 69 and the power of each DGU can be chosen to be from 0.0 MW to 2.0 MW. Since this research is being considered under harmonic conditions, the harmonic fows shown in Table 2 are injected simultaneously into nine buses 12, 18, 19, 22, 25, 34, 46, 56, and 65.

Case 2: IEEE 69-Bus Radial Distribution System.
Te worst, mean, and best ftness values for the 50 trial runs are clearly presented in Table 5. According to the collected data from Table 5, the best ftness values for ABC, SSA, EO, and IEO are 0.3279, 0.3267, 0.3264, and 0.3262, respectively. Te best ftness value of IEO is smaller than that of all other methods. Clearly, the modifcations carried out on IEO have a positive impact on the efectiveness of the method and they support IEO reach the most optimal solution among the four run methods. To prove the high stability of IEO, the mean of ffty trial runs and the ffty runs in detail are also reported in Table 5 and Figure 12. Te mean value of IEO is 0.3267, and it is smaller than 0.3317, 0.3298, and 0.3269 from ABC, SSA, and EO, respectively. In Figure 12, the ffty ftness values of the methods are re-sorted and the curve of IEO is almost below the curves of other methods. Furthermore, this curve has very small fuctuations while others have signifcant fuctuations. In summary, IEO can reach the best performance and the most stable search ability among the four run methods. Table 6 shows a summary of two objectives and corresponding power loss and cost obtained by implemented and compared methods. Te frst objective and the corresponding power loss of IEO are 0.3206 and 0.07197 MW, whereas those of the second-best method (EO) and the worst method (GA [39]) are (0.3207 and 0.0720 MW) and (0.3962 and 0.0889 MW), respectively. Te power loss on branches in the base system and the hybrid system with DGUs by applying the optimal solution from IEO is clearly presented in Figure 13. It has been shown that the power loss in the IEEE 69-bus system before and after the integration of DGUs is signifcantly diferent. Specifcally, the power loss reduction is greater on branches 4 to 9, 11 to 14, and 52 to 60. Te loss reduction on other remaining branches is much smaller than these branches. Overall, thanks to the proper installation of DGUs, branch losses are drastically reduced for more economic and technical benefts.
Te second objective and the corresponding cost of IEO are 0.3487 and $3.2840 million, respectively. Tese values are higher than those of MOBA [34], BFOA [41], HSA [43], and SSA but lower than other methods such as GA [39], PSO [39], GA-PSO [39], LSF-SA [40], ABC, and EO. Excluding the comparison with MOBA [34], BFOA [41], HSA [43], and SSA, IEO can save up to $1.4209 million and $0.0028 million when compared to the worst method (GA [39]) and the second-best method (EO). Comparing IEO with MOBA, BFOA, HSA, and SSA, the optimal solution from the proposed IEO method cannot achieve better cost, but it can get much better power loss reduction. As mentioned, this is a trade-of in solving the multi-objective function problem. Te results are afected by the selected weighting factors, and it is unavoidable when compared with many other methods.  However, the multi-objective function of IEO with 0.3262 is smaller than these methods with 0.3281, 0.3376, 0.3684, and 0.3267. Tus, IEO is actually more efective than others in solving the problem of the location of DGUs in the IEEE 69bus system.
Using the optimal solution obtained by the proposed IEO method, THD and IHD of each bus in the base system without DGUs and modifed system with DGUs are plotted in Figure 14. Data from the fgure indicate that the maximum values of THD and IHD of the base system are 5.4287% and 3.5087%, but those from the modifed system with the application of IEO are 3.0232% and 1.9495%, respectively. In addition, the view on THD and IHD of the whole 69 buses also sees the signifcant improvement of the modifed system from buses 5 to 28 and from buses 51 to 69. In another view, the improvement of voltage is also presented in Figure 15. Te minimum bus voltage is 0.9092 pu for the base system but is much higher and equal to 0.9761 pu for the modifed system with DGU connection. Voltage values at buses 5 to 28 and 51 to 69 in the modifed system with DGUs are much higher than those in the base system without DGUs. Clearly, the installation of DGUs in the distribution system can reach an extra beneft in reducing harmonic distortions and improving voltage profle.
Optimal solutions of the applied system are reported in Table S2 in Supplementary Materials. Figure 16 shows the confguration of the IEEE 85-bus radial distribution system. Te load data and line data of the system     Table 2 are also injected to seven buses 5, 12, 25, 34, 56, 65, and 77. Te best, mean, and worst ftness values of ABC, SSA, EO, and IEO are presented in Table 7. IEO is also still the best method with all smallest values excluding the comparison of the worst ftness value with EO. Te ffty runs are also presented in detail in Figure 17 by sorting ftness from the smallest to the highest values shown in diferent color curves. Almost all points on red curve of IEO are lower than other points on three other curves of ABC, SSA, and EO. So, IEO is the most powerful method among the four applied ones for the system. Table 8 presents the results obtained by the four methods in detail. Tere are no comparisons with other published methods for the system because there were no studies regarding the IEEE 85-bus system before. Te frst objective and the corresponding power loss of IEO are, respectively, 0.4792 and 0.1515 MW, which are better than those of others. In other words, the solution of IEO can reduce the power loss from 0.3161 MW to 0.1515 MW while ABC, SSA, and EO can reduce the loss to a higher value, 0.1553 MW, 0.1524 MW, and 0.1521 MW, respectively. It has been strongly proven that the solution from IEO is more efective in terms of power loss reduction than others. To see the efectiveness of loss reduction thanks to the solution obtained by IEO, the power losses on all branches in the IEEE 85-bus system for two cases, with and without DGUs, are plotted in Figure 18. Tere is outstanding power loss reduction on the frst branches, from branch 1 to branch 7. Especially, the highest loss reduction is reached on branch 7 and the loss reduces from higher than 100 kW to under 50 kW. Although the loss reduction on branches 1 to 6 is not as efective as on branch 7, the reduction is relatively high, from some kW to 20 kW. Te loss reduction on other branches is not as efective or even there is no loss reduction on some branches. Te loss reduction on branches 24 to 29, branch 56, and branch 57 is still high while the reduction on other remaining branches is not identifed. In general, power loss reduction is efective in all branches, but the total loss reduction is still high, from 0.3161 MW to 0.1515 MW. Te high loss reduction is only for one hour, and it indicates that the placement of DGUs in distribution system is very useful.

Case 3: IEEE 85-Bus Radial Distribution System.
Te second objective and the corresponding cost of IEO are 0.3248 and $3.0593 million, which are slightly better than other compared methods. As presented, this is a trade-of in solving the multi-objective function problem.
By applying the best solution from IEO, THD and IHD of each bus and voltage profle in the base system without DGUs and modifed system with DGUs are plotted in Figures 19 and 20, respectively. Te maximum values of THD and IHD in the base system are 6.1504% and 3.9840%, but those from the modifed system with the application of IEO are 4.0445% and 2.6233%, respectively. Approximately all buses in the modifed system can reach much better THD and IHD than those in the base system. Te base system has the smallest voltage with 0.8713 pu, which is much smaller than the permissible voltage limits with [0.95, 1.05] pu. However, after the three DGUs are connected to the system, the voltage profle is improved drastically with the lowest voltage value of 0.9500 pu. Te view on voltage of all buses in the two cases can indicate that the modifed system can have a much better voltage profle than the base system. Clearly, the installation of three DGUs can cut the harmonic distortion and increase the voltage of buses for the IEEE 85-bus system. Finally, these outstanding results show that the proposed method (IEO) is indeed more powerful than others in solving the problem of DGU placement in the considered distribution system.
Optimal solutions of the applied system are reported in Table S3 in Supplementary Materials.

Shortcomings of the Proposed Method and Research Expansion
In this paper, the proposed method (IEO) has outstanding performance over EO, ABC, and SSA. As shown in Figures 7,12, and 17, IEO can reach many better solutions than EO, ABC, and SSA. Furthermore, IEO also has very small fuctuations, but the fuctuations from ABC, SSA, and EO are very high. As compared to previous methods in other studies, Tables 4 and 6 indicate that the proposed IEO also reaches better solutions with smaller cost and loss. For the three test systems, IEO is really efective, but it cannot avoid shortcomings such as solving larger and more complicated systems. In fact, IEO is successful in fnding optimal solutions of the problems because the number of control variables is not high. It is 6 for all three test systems, including three locations and three rated powers. In addition, the variable boundaries are not large. Te location can be from 2 to 33 for the frst system, to 69 for the second system, and to 85 for the last system, and the rated powers can be from 0 to 2.0 MW. IEO is based on four best solutions to update new solutions, and almost all newly obtained solutions are near the four best solutions. So, local search has a huge contribution to the promising results obtained by IEO. Tis feature can be the highest limitation of IEO when it is applied for large-scale problems with a high number of control variables and very large search spaces. Te second limit of IEO is convergence speed to high-quality solutions. Although IEO is tested on the three test systems with only six control variables and small search space, it needs high enough values for control parameters. Population size is set to 30 while iteration number is 200 for the frst system, 230 for the second system, and 250 iterations for the last system. Due to the two limitations, IEO can be inefective for larger and more complicated problems, and it may need more improvements. Additionally, in this study, we supposed the highest power of DGUs is 2.0 MW, and DGUs can change in the range between 0 and 2.0 MW. Tis assumption may not be true in practice. So, in the future, the authors will consider    16 Complexity the uncertainties of solar radiation as well as wind speed for solar photovoltaic units and wind turbines. Realistically, the power of solar photovoltaic units and wind turbines mainly depends on the primary energies [44,45]. For considering these uncertainties, they can be based on recorded historical data and combined with the probability density function to predict the power for solar photovoltaic units and wind turbines [46]. For the cases that systems change topology because of expansion or reconfguration, the application of the proposed IEO or other metaheuristic algorithms to solve the optimization problem is the same as long as the grid data     such as load and line data are given or calculated at each computation iteration. For other cases that the systems need the installation of capacitors and/or wind turbines, the problem will have more control variables, including the location and rated power of the capacitors and/or wind turbines. It is noted that rated power of capacitors is reactive power while that of wind turbines is comprised of active and reactive power [45]. Ten, the additional control variables are put in the data and the forward/backward sweep technique is run to reach dependent variables as those in the current problem, including current, voltage, and individual and total harmonic distortions. Tese independent variables are checked and penalized if they are beyond predetermined boundaries. Te objective function and the penalty terms

Conclusions
Tis research proposed improved equilibrium optimizer (IEO) to reach higher solutions and more stable search performance than other methods in solving the optimization problem of installing DGUs in three distribution systems for maximizing the economic and technical benefts. Overall, the main contributions of the study can be briefy summarized as follows: (i) Te proposed method could reach more promising solutions for the problem than conventional EO. Realistically, the average values of the ftness function of IEO for the three systems were 0.3856, 0.3267, and 0.4614 while these values of the secondbest method (EO) were 0.3861, 0.3269, and 0.4633, respectively. As compared from the obtained results, the proposed method is not only better than the original method but also better than other implemented methods (ABC, SSA, and EO) and previously published methods (MOBA, GA, PSO, GA-PSO, LSF-SA, HSA, BFOA, and BSOA). Tis proves that the improvements in IEO have been remarkably efective. (ii) For the frst objective of minimizing power loss, by applying the optimal solution of IEO, it could reduce the loss from 0.2110 MW to 0.0815 MW for the frst system, from 0.2245 MW to 0.07197 MW for the second system, and from 0.3161 MW to 0.1515 MW for the third system. In addition, the second objective was also optimized efecitvely. Te total cost over 20 years was only $3.4753 million, $3.2840 million, and $3.0593 million for the three systems, respectively. (iii) Te proposed solutions from IEO also satisfed all the constraints. It has given excellent optimal solutions for placing DGUs that can mitigate harmonic distortions and keeps it in IEEE Std. 519. In addition, the voltage profles were also signifcantly enhanced with the voltage range from [0 Clearly, the contributions of the placement of DGUs are signifcant for distribution systems. However, DGUs will bring more benefts to the systems if DGUs are combined with smart inverters and battery energy store system (BESS). Smart inverters can stabilize the operation of DGUs, meanwhile BESS can store energy for the case that the solaror wind-based DGUs produce higher power than load demand. Te optimization operation of the distribution systems with the presence of DGUs, smart inverters, and BESS problem is an upcoming direction of the study.

ABC:
Artifcial bee colony algorithm BBO: Biogeography-based optimization BFOA: Bacterial foraging optimization algorithm BPSO-SLFA: Binary particle swarm optimization and shufed frog leap algorithm BSOA: Backtracking search optimization algorithm DGU(s): Distributed generation unit(s) GA: Genetic algorithm HSA: Harmony search algorithm LSF-SA: Loss sensitivity factor and simulated annealing algorithm MOBA: Multi-objective bat algorithm PCC: Point at common coupling PSO: Particle swarm optimization RNN: Recurrent      Number of loads in the system θ 1 , θ 2 and l: Random numbers in the range of (0, 1) R bh : Resistance of the bh th branch S i : Te i th solution

Data Availability
Data of the employed systems were extracted from [20,31].

Conflicts of Interest
Te authors declare that there are no conficts of interest.