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

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

---

## Abstract

This paper proposes an improved equilibrium optimizer (IEO) for determining optimal location and effective 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 flows 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. The proposed IEO is developed from the original equilibrium optimizer (EO), which was motivated by control volume mass balance models. This novel algorithm can effectively 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 significantly reduced from 0.2110 MW to 0.0815 MW, 0.2245 MW to 0.07197 MW, and 0.3161 MW to 0.1515 MW 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. The performance of the proposed algorithm is compared to three other implemented methods consisting of artificial bee colony (ABC) algorithm, salp swarm algorithm (SSA), and EO and eight previously published methods for the three test systems. The comparisons of results imply that IEO is better than other methods in terms of performance, stability, and convergence characteristics.


---

## Body

## 1. 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 effects such as increasing power losses, decreasing voltage quality, increasing operation costs and other costs, and so on [2]. Thus, optimal planning strategy of DGUs in the distribution systems is essential to maximize the economic and technical achievements.The achieved benefits 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–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 ineffective solutions for the optimization problem. To overcome the disadvantages of GA, other authors in [6–8] proposed hybrid methods which divide tasks of each phase in the considered problem to apply the most appropriate algorithm to each phase. Those 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 finding the effective solution in solving the problem of determining placement and capacity of DGUs [9–11]. These authors optimally determined the penetration of multiple DGUs into distribution systems. As a result, the voltage profile was greatly improved and losses on branches were effectively cut. However, this algorithm has reached a low convergence rate and unguaranteed stability. Thus, many researchers have sought to improve the original PSO. The studies [12, 13] proposed improvements by integrating uniform distribution, exponential attenuation inertia weight, cubic spline interpolation function, and learning factor of enhanced control. These studies have significantly contributed to a high efficiency of the improved algorithms. In addition to the abovementioned methods, ABC is also widely used for optimizing the installation of DGUs [14, 15]. The authors have proved the effects of different 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 efficient 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]. The authors demonstrated the obtained multi-benefit for the grid connected DGU system. This 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. The results from simulation have also proved that harmonic distortions could be significantly mitigated to IEEE Std. 519 after integrating DGUs into the considered systems. Besides, with the goal of enhancing distribution system capacity without causing adverse effects, 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 flow in the DGU integrated system [22, 23]. In addition, a few other studies have analyzed the benefits 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]. These authors have suggested the multi-objective evolutionary algorithm which is based on the Pareto optimality, called SPEA2, to solve the optimization problem of network reconfiguration 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 finding 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. The leader is the leading salp and stands in front of the chain. Other remaining salps are considered as the followers. The leader determines the direction of movement for followers. By comparing the obtained results with some other metaheuristic algorithms, the study [29] showed the superior effectiveness of SSA in solving the single and multi-objective 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 sunflowers’ motion toward the sun to get the radiation. Sunflowers that are located near the sun are calmer because greater heat is received in this space. Conversely, flowers that are further away from the sun tend to move closer to the sun because they receive less heat. This cycle is repeated daily in the morning [31]. Based on this natural feature, the algorithm has been built and applied in finding the optimal solutions effectively. Besides, a recently effective method, called EO, was also published in 2020 [32]. This 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 finding the feasible solution in the space, and they are considered as the search agents. The 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 effectiveness 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). The suggested IEO can inherit the strengths of EO and effectively avoid disadvantages of EO to reach a greater performance than EO. Thus, IEO is employed to find 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. The 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 flow, 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 multi-objective function [31]. In this technique, each weighting coefficient is assigned to each single objective function and the value of the weighting coefficient 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 effectiveness are power loss and costs of investment, operation, and maintenance for the installed DGUs. Therefore, this study proposes to consider the multi-objective function under constraints of voltage profile, branch current, and harmonic distortions in different distribution systems of IEEE 33-bus, IEEE 69-bus, and IEEE 85-bus radial distribution systems. The 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)
The gained benefits 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 efficient and stable optimal method will contribute to the both economic and technical benefits. Therefore, this paper has developed a novel algorithm (called IEO) based on modifications on equations updating solutions in the original algorithm (EO). Thereby, IEO has significant improvements in finding feasible solutions with high quality while EO is less effective when testing them on IEEE 33-bus, IEEE 69-bus, and IEEE 85-bus radial distribution systems.The rest of the paper is organized as follows.Firstly, section of the problem formulation (Section2) 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. Thirdly, 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. The 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.


## 2. Problem Formulation

### 2.1. Objective Function

Total power loss reduction on conductors is one of the key factors that mainly affect the economic and technical problems of distribution systems. Thus, 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. Thus, 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).


#### 2.1.1. Total Active Power Loss (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. Therefore, determining the optimal strategy in connecting DGUs to the system is a great way for the power loss minimization. The reduction of power loss mainly contributes to the effectiveness 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 first 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]:(1)TAPLDGU=∑bh=1NbhIbhDGU2Rbh,where TAPLDGU is the total power loss after connecting DGUs to the system.However, the compromise solution in the multi-objective function is difficult to determine because the single objective values have large deviations. Therefore, it is necessary to convert the objectives to the same considered range. Specifically, in the first single objective, we suggest the normalization equation and it is presented as equation (2) for determining the best compromise solution [20].(2)MinimizeOFI=TAPLDGUTAPLNoDGU,where TAPLNoDGU is the total power loss of the original system without DGUs and it can be found by using the equation below [17]:(3)TAPLNoDGU=∑bh=1NbhIbh2Rbh.If DGUs are successfully integrated into the distribution system, the value ofTAPLDGU will be smaller than the value of TAPLNoDGU. In other words, the value of OFI will vary in the range from 0 to 1 and the small value of the objective is expected. This demonstrates the positive effect of connecting DGUs in the system in reducing the total power loss.


#### 2.1.2. Total Cost of DGUs (TCDG)

In order to improve technical efficiency and maximize profits 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. The initial investment cost of all DGUsICDGU,dopt$$)) and the total cost of maintenance and operation sector MOCDGU,dopt$$)) are mainly influenced by the penetration level of DGUs APDGU,doptMW. The larger the penetration level is, the higher the above costs are. In addition, ICDGU,dopt$$) and MOCDGU,dopt$$) also depend on the current market on the investment price (CdinvM/MW) as well as maintenance and operation price Cdm&o/MWh for each connected DGU, respectively. ICDGU,dopt$$) and MOCDGU,dopt$$) can be found by [34](4)ICDGU,dopt=∑d=1NdgCdinv.APDGU,dopt,(5)MOCDGU,dopt=8760.∑y=1PP∑dNdgαy.Cdm&o.APDGU,dopt.In equation (5), MOCDGU,dopt$$) is considered for a 20-year planning period (i.e., PP = 20) and αy is the cumulative present value factor and obtained by(6)αy=11+dry,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:(7)TCDGDGUopt=ICDGU,dopt+MOCDGU,dopt.Similar to the first objective, the second objective (OFII) must be also converted into the same range from 0 to 1 for evaluating the multi-objective function correctly. Thus, the second objective can be normalized as follows:(8)MinimizeOFII=TCDGDGUoptTCDGDGUmax.Here,TCDGDGUmax is the total cost of the maximum generated power of DGUs in the distribution system. The value of TCDGDGUmax can be found by using equations (9)–(11) [34]:(9)ICDGU,dmax=∑d=1NdgCdinv.APDGU,dmax,(10)MOCDGU,dmax=8760.∑y=1PP∑dNdgαy.Cdm&o.APDGU,dmax,(11)TCDGDGUmax=ICDGU,dmax+MOCDGU,dmax,where ICDGU,dmax$$) is the maximum initial investment cost of all DGUs and MOCDGU,dmax 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. Therefore, OFII always varies in the range from 0 to 1. The smaller OFII is, the more profitable 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:(12)Minimize OF=ωI.OFI+ωII.OFII,where OFI and ωI are the first normalized single target and the weighting factor associated to the target. OFII 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](13)ωI+ωII=1and0≤ωI,ωII≤1.The values ofωI and ωII are chosen dependent on the importance of the components (single objective) in the multi-objective function [34]. These values can be adjusted by user. The larger the weighting factor is, the greater its significance in the multi-objective function is.


### 2.2. Constraints

#### 2.2.1. Power Balance Constraints

The 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 or decreased, respectively [35]. Thus, it is necessary to take the power balance between power generation and power consumption sides into account. The power generation side includes the supplied power from the grid APgr and the injected power by connected DGUs APDGU,d; meanwhile, the total power consumption side is comprised of the consumed power by loads APld,b and the power losses APls,bh. These two power sides should be balanced as follows [31]:(14)APgr=∑b=1NldAPld,b+∑bh=1NbhAPls,bh−∑d=1NdgAPDGU,d.


#### 2.2.2. Voltage Limits

The voltage profile 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. Thus, 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, 69-bus, and 85-bus radial distribution systems, the best range for the bus voltage is maintained within the minimum limitVmin and the maximum limit Vmax 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. The bus voltage and voltage limits are considered at the fundamental frequency and constrained as follows [17]:(15)Vmin≤Vb≤Vmax,b=1,…,Nb.


#### 2.2.3. 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 defined as the deformation of a waveform by integer multiples of a fundamental frequency. As shown in Figure1, they are the distorted waveforms of the 12-pulse converter, AC voltage regulator, and fluorescent lighting [35].Figure 1
The waveform of some disturbing loads.Harmonic distortions are caused by nonlinear load’s operations (saturated electric machines or rectifiers), and they can be minimized by using passive and active filters. To evaluate the power quality of the systems, the two parameters including total harmonic distortion and individual harmonic distortion must satisfy allowable ranges in Table1 [36, 37].Table 1
The harmonic voltage limits for power producer.
Bus voltage at PCCIndividual harmonic distortion (IHD) (%)Total harmonic distortion (THD) (%)69 kV and below3.05.069 kV to 161 kV1.52.5169 kV and above1.01.5The two factors which indicate the harmonic noise level of the system are total voltage harmonic distortion (THDb%) and individual voltage harmonic distortion IHDbh%. Their values can be obtained by voltage values at the fundamental frequency Vb1 and other higher order frequencies Vbh as follows [20]:(16)THDb%=∑h≠1HVbh2Vb1×100≤THDmax%,(17)IHDbh%=VbhVb1×100≤IHDmax%.By applying the IEEE Std. 519 of harmonic limits as shown in Table1, THDmax% and IHDmax% 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].


#### 2.2.4. DGU’s Capacity Limits

Before photovoltaic strategy planning, for each DGU, the lower and upper bounds of locationLDGUminandLDGUmax as well as capacity APDGUminandAPDGUmax should be predetermined and imposed on DGUs as shown in equations (18) and (19), respectively [4]. In other words, the selection of location LDGU,d and rated power (APDGU,d) of the dth 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 APld as equation (20) for this study [12, 20].(18)LDGUmin≤LDGU,d≤LDGUmax,(19)APDGUmin≤APDGU,d≤APDGUmax,(20)∑d=1NdgAPDGU,d≤80%×APld.


#### 2.2.5. Branch Current Limits

When DGUs are connected to the distribution system, the current in the branches can be significantly increased depending on the location and power of DGUs. This can disrupt the original structure of the grid and cause damage to existing conductors and other electric components. Therefore, the branch current limit should be considered as one of the most serious constraints in systems after integrated DGUs as follows [28]:(21)IbhDGU≤Ibhmax,bh=1,…,Nbh.


## 2.1. Objective Function

Total power loss reduction on conductors is one of the key factors that mainly affect the economic and technical problems of distribution systems. Thus, 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. Thus, 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).


### 2.1.1. Total Active Power Loss (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. Therefore, determining the optimal strategy in connecting DGUs to the system is a great way for the power loss minimization. The reduction of power loss mainly contributes to the effectiveness 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 first 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]:(1)TAPLDGU=∑bh=1NbhIbhDGU2Rbh,where TAPLDGU is the total power loss after connecting DGUs to the system.However, the compromise solution in the multi-objective function is difficult to determine because the single objective values have large deviations. Therefore, it is necessary to convert the objectives to the same considered range. Specifically, in the first single objective, we suggest the normalization equation and it is presented as equation (2) for determining the best compromise solution [20].(2)MinimizeOFI=TAPLDGUTAPLNoDGU,where TAPLNoDGU is the total power loss of the original system without DGUs and it can be found by using the equation below [17]:(3)TAPLNoDGU=∑bh=1NbhIbh2Rbh.If DGUs are successfully integrated into the distribution system, the value ofTAPLDGU will be smaller than the value of TAPLNoDGU. In other words, the value of OFI will vary in the range from 0 to 1 and the small value of the objective is expected. This demonstrates the positive effect of connecting DGUs in the system in reducing the total power loss.


### 2.1.2. Total Cost of DGUs (TCDG)

In order to improve technical efficiency and maximize profits 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. The initial investment cost of all DGUsICDGU,dopt$$)) and the total cost of maintenance and operation sector MOCDGU,dopt$$)) are mainly influenced by the penetration level of DGUs APDGU,doptMW. The larger the penetration level is, the higher the above costs are. In addition, ICDGU,dopt$$) and MOCDGU,dopt$$) also depend on the current market on the investment price (CdinvM/MW) as well as maintenance and operation price Cdm&o/MWh for each connected DGU, respectively. ICDGU,dopt$$) and MOCDGU,dopt$$) can be found by [34](4)ICDGU,dopt=∑d=1NdgCdinv.APDGU,dopt,(5)MOCDGU,dopt=8760.∑y=1PP∑dNdgαy.Cdm&o.APDGU,dopt.In equation (5), MOCDGU,dopt$$) is considered for a 20-year planning period (i.e., PP = 20) and αy is the cumulative present value factor and obtained by(6)αy=11+dry,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:(7)TCDGDGUopt=ICDGU,dopt+MOCDGU,dopt.Similar to the first objective, the second objective (OFII) must be also converted into the same range from 0 to 1 for evaluating the multi-objective function correctly. Thus, the second objective can be normalized as follows:(8)MinimizeOFII=TCDGDGUoptTCDGDGUmax.Here,TCDGDGUmax is the total cost of the maximum generated power of DGUs in the distribution system. The value of TCDGDGUmax can be found by using equations (9)–(11) [34]:(9)ICDGU,dmax=∑d=1NdgCdinv.APDGU,dmax,(10)MOCDGU,dmax=8760.∑y=1PP∑dNdgαy.Cdm&o.APDGU,dmax,(11)TCDGDGUmax=ICDGU,dmax+MOCDGU,dmax,where ICDGU,dmax$$) is the maximum initial investment cost of all DGUs and MOCDGU,dmax 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. Therefore, OFII always varies in the range from 0 to 1. The smaller OFII is, the more profitable 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:(12)Minimize OF=ωI.OFI+ωII.OFII,where OFI and ωI are the first normalized single target and the weighting factor associated to the target. OFII 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](13)ωI+ωII=1and0≤ωI,ωII≤1.The values ofωI and ωII are chosen dependent on the importance of the components (single objective) in the multi-objective function [34]. These values can be adjusted by user. The larger the weighting factor is, the greater its significance in the multi-objective function is.


## 2.1.1. Total Active Power Loss (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. Therefore, determining the optimal strategy in connecting DGUs to the system is a great way for the power loss minimization. The reduction of power loss mainly contributes to the effectiveness 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 first 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]:(1)TAPLDGU=∑bh=1NbhIbhDGU2Rbh,where TAPLDGU is the total power loss after connecting DGUs to the system.However, the compromise solution in the multi-objective function is difficult to determine because the single objective values have large deviations. Therefore, it is necessary to convert the objectives to the same considered range. Specifically, in the first single objective, we suggest the normalization equation and it is presented as equation (2) for determining the best compromise solution [20].(2)MinimizeOFI=TAPLDGUTAPLNoDGU,where TAPLNoDGU is the total power loss of the original system without DGUs and it can be found by using the equation below [17]:(3)TAPLNoDGU=∑bh=1NbhIbh2Rbh.If DGUs are successfully integrated into the distribution system, the value ofTAPLDGU will be smaller than the value of TAPLNoDGU. In other words, the value of OFI will vary in the range from 0 to 1 and the small value of the objective is expected. This demonstrates the positive effect of connecting DGUs in the system in reducing the total power loss.


## 2.1.2. Total Cost of DGUs (TCDG)

In order to improve technical efficiency and maximize profits 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. The initial investment cost of all DGUsICDGU,dopt$$)) and the total cost of maintenance and operation sector MOCDGU,dopt$$)) are mainly influenced by the penetration level of DGUs APDGU,doptMW. The larger the penetration level is, the higher the above costs are. In addition, ICDGU,dopt$$) and MOCDGU,dopt$$) also depend on the current market on the investment price (CdinvM/MW) as well as maintenance and operation price Cdm&o/MWh for each connected DGU, respectively. ICDGU,dopt$$) and MOCDGU,dopt$$) can be found by [34](4)ICDGU,dopt=∑d=1NdgCdinv.APDGU,dopt,(5)MOCDGU,dopt=8760.∑y=1PP∑dNdgαy.Cdm&o.APDGU,dopt.In equation (5), MOCDGU,dopt$$) is considered for a 20-year planning period (i.e., PP = 20) and αy is the cumulative present value factor and obtained by(6)αy=11+dry,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:(7)TCDGDGUopt=ICDGU,dopt+MOCDGU,dopt.Similar to the first objective, the second objective (OFII) must be also converted into the same range from 0 to 1 for evaluating the multi-objective function correctly. Thus, the second objective can be normalized as follows:(8)MinimizeOFII=TCDGDGUoptTCDGDGUmax.Here,TCDGDGUmax is the total cost of the maximum generated power of DGUs in the distribution system. The value of TCDGDGUmax can be found by using equations (9)–(11) [34]:(9)ICDGU,dmax=∑d=1NdgCdinv.APDGU,dmax,(10)MOCDGU,dmax=8760.∑y=1PP∑dNdgαy.Cdm&o.APDGU,dmax,(11)TCDGDGUmax=ICDGU,dmax+MOCDGU,dmax,where ICDGU,dmax$$) is the maximum initial investment cost of all DGUs and MOCDGU,dmax 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. Therefore, OFII always varies in the range from 0 to 1. The smaller OFII is, the more profitable 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:(12)Minimize OF=ωI.OFI+ωII.OFII,where OFI and ωI are the first normalized single target and the weighting factor associated to the target. OFII 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](13)ωI+ωII=1and0≤ωI,ωII≤1.The values ofωI and ωII are chosen dependent on the importance of the components (single objective) in the multi-objective function [34]. These values can be adjusted by user. The larger the weighting factor is, the greater its significance in the multi-objective function is.


## 2.2. Constraints

### 2.2.1. Power Balance Constraints

The 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 or decreased, respectively [35]. Thus, it is necessary to take the power balance between power generation and power consumption sides into account. The power generation side includes the supplied power from the grid APgr and the injected power by connected DGUs APDGU,d; meanwhile, the total power consumption side is comprised of the consumed power by loads APld,b and the power losses APls,bh. These two power sides should be balanced as follows [31]:(14)APgr=∑b=1NldAPld,b+∑bh=1NbhAPls,bh−∑d=1NdgAPDGU,d.


### 2.2.2. Voltage Limits

The voltage profile 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. Thus, 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, 69-bus, and 85-bus radial distribution systems, the best range for the bus voltage is maintained within the minimum limitVmin and the maximum limit Vmax 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. The bus voltage and voltage limits are considered at the fundamental frequency and constrained as follows [17]:(15)Vmin≤Vb≤Vmax,b=1,…,Nb.


### 2.2.3. 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 defined as the deformation of a waveform by integer multiples of a fundamental frequency. As shown in Figure1, they are the distorted waveforms of the 12-pulse converter, AC voltage regulator, and fluorescent lighting [35].Figure 1
The waveform of some disturbing loads.Harmonic distortions are caused by nonlinear load’s operations (saturated electric machines or rectifiers), and they can be minimized by using passive and active filters. To evaluate the power quality of the systems, the two parameters including total harmonic distortion and individual harmonic distortion must satisfy allowable ranges in Table1 [36, 37].Table 1
The harmonic voltage limits for power producer.
Bus voltage at PCCIndividual harmonic distortion (IHD) (%)Total harmonic distortion (THD) (%)69 kV and below3.05.069 kV to 161 kV1.52.5169 kV and above1.01.5The two factors which indicate the harmonic noise level of the system are total voltage harmonic distortion (THDb%) and individual voltage harmonic distortion IHDbh%. Their values can be obtained by voltage values at the fundamental frequency Vb1 and other higher order frequencies Vbh as follows [20]:(16)THDb%=∑h≠1HVbh2Vb1×100≤THDmax%,(17)IHDbh%=VbhVb1×100≤IHDmax%.By applying the IEEE Std. 519 of harmonic limits as shown in Table1, THDmax% and IHDmax% 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].


### 2.2.4. DGU’s Capacity Limits

Before photovoltaic strategy planning, for each DGU, the lower and upper bounds of locationLDGUminandLDGUmax as well as capacity APDGUminandAPDGUmax should be predetermined and imposed on DGUs as shown in equations (18) and (19), respectively [4]. In other words, the selection of location LDGU,d and rated power (APDGU,d) of the dth 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 APld as equation (20) for this study [12, 20].(18)LDGUmin≤LDGU,d≤LDGUmax,(19)APDGUmin≤APDGU,d≤APDGUmax,(20)∑d=1NdgAPDGU,d≤80%×APld.


### 2.2.5. Branch Current Limits

When DGUs are connected to the distribution system, the current in the branches can be significantly increased depending on the location and power of DGUs. This can disrupt the original structure of the grid and cause damage to existing conductors and other electric components. Therefore, the branch current limit should be considered as one of the most serious constraints in systems after integrated DGUs as follows [28]:(21)IbhDGU≤Ibhmax,bh=1,…,Nbh.


## 2.2.1. Power Balance Constraints

The 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 or decreased, respectively [35]. Thus, it is necessary to take the power balance between power generation and power consumption sides into account. The power generation side includes the supplied power from the grid APgr and the injected power by connected DGUs APDGU,d; meanwhile, the total power consumption side is comprised of the consumed power by loads APld,b and the power losses APls,bh. These two power sides should be balanced as follows [31]:(14)APgr=∑b=1NldAPld,b+∑bh=1NbhAPls,bh−∑d=1NdgAPDGU,d.


## 2.2.2. Voltage Limits

The voltage profile 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. Thus, 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, 69-bus, and 85-bus radial distribution systems, the best range for the bus voltage is maintained within the minimum limitVmin and the maximum limit Vmax 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. The bus voltage and voltage limits are considered at the fundamental frequency and constrained as follows [17]:(15)Vmin≤Vb≤Vmax,b=1,…,Nb.


## 2.2.3. 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 defined as the deformation of a waveform by integer multiples of a fundamental frequency. As shown in Figure1, they are the distorted waveforms of the 12-pulse converter, AC voltage regulator, and fluorescent lighting [35].Figure 1
The waveform of some disturbing loads.Harmonic distortions are caused by nonlinear load’s operations (saturated electric machines or rectifiers), and they can be minimized by using passive and active filters. To evaluate the power quality of the systems, the two parameters including total harmonic distortion and individual harmonic distortion must satisfy allowable ranges in Table1 [36, 37].Table 1
The harmonic voltage limits for power producer.
Bus voltage at PCCIndividual harmonic distortion (IHD) (%)Total harmonic distortion (THD) (%)69 kV and below3.05.069 kV to 161 kV1.52.5169 kV and above1.01.5The two factors which indicate the harmonic noise level of the system are total voltage harmonic distortion (THDb%) and individual voltage harmonic distortion IHDbh%. Their values can be obtained by voltage values at the fundamental frequency Vb1 and other higher order frequencies Vbh as follows [20]:(16)THDb%=∑h≠1HVbh2Vb1×100≤THDmax%,(17)IHDbh%=VbhVb1×100≤IHDmax%.By applying the IEEE Std. 519 of harmonic limits as shown in Table1, THDmax% and IHDmax% 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].


## 2.2.4. DGU’s Capacity Limits

Before photovoltaic strategy planning, for each DGU, the lower and upper bounds of locationLDGUminandLDGUmax as well as capacity APDGUminandAPDGUmax should be predetermined and imposed on DGUs as shown in equations (18) and (19), respectively [4]. In other words, the selection of location LDGU,d and rated power (APDGU,d) of the dth 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 APld as equation (20) for this study [12, 20].(18)LDGUmin≤LDGU,d≤LDGUmax,(19)APDGUmin≤APDGU,d≤APDGUmax,(20)∑d=1NdgAPDGU,d≤80%×APld.


## 2.2.5. Branch Current Limits

When DGUs are connected to the distribution system, the current in the branches can be significantly increased depending on the location and power of DGUs. This can disrupt the original structure of the grid and cause damage to existing conductors and other electric components. Therefore, the branch current limit should be considered as one of the most serious constraints in systems after integrated DGUs as follows [28]:(21)IbhDGU≤Ibhmax,bh=1,…,Nbh.


## 3. The Proposed Method

### 3.1. 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. These search agents update their concentration randomly with the desire of finding 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. The quality of the solution entirely depends on the fitness function value under the constraints, and the use of the fitness function result is intended to evaluate the mass-concentration balance within the adjusted volume. Overall, EO is an ideal algorithm for exploration, finding high-quality optimal solutions and avoiding local traps. The algorithm is expressed as follows.


#### 3.1.1. Initial Population Generation

In EO, each solutionFi in the population contains a set of adjustment variables and there must be a predetermined range for the variables. The lower and upper bounds of the variables are called the minimum solution (Fmin) and the maximum solution (Fmax), respectively. During the application of EO, each solution Fi is always kept within these limits. The mentioned solutions and the boundary solutions are formulated as follows:(22)Fi=Vaji;j=1,2,3,…,Nad&i=1,2,3,…,Nc,(23)Fmin=Vajmin&Fmax=Vajmax;j=1,…,Nad.The generation of initial solution set or initial population is implemented by [32](24)Fi=Fmin+rand0,1.Fmax−Fmin;i=1,…,Nc.


#### 3.1.2. Calculation Process for Updating Concentration

Quality of all solutions is calculated to rank the solutions in population in which solution with the smallest fitness value is the best and other solutions with worse fitness values are from the second-best solution to the worst solution. After obtaining the fitness values of solutions, the top four best solutions are collected and assigned toFLbest1, FLbest2, FLbest3, FLbest4 for the calculation of the arithmetic mean (FLmean). The new solution Finew is updated through a combination of three components. The first component (Fs) is randomly selected from the solution set including FLbest1, FLbest2, FLbest3, FLbest4, and FLmean. 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 difference between the selected solution (Fs) and the considered old solution (Fi). This 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. The three components are useful in finding better solutions, which are not close to the current good solutions. Finally, the new solution Finew can be calculated by applying the following equation [32].(25)Finew=Fs+Fi−FsET+GRl1−ET.In equation (25), two important terms, which mainly affect 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].(26)ET=β1signµ−0.5e−tc.l−1.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 coefficient tc is a value that changes with each iteration, and its value depends entirely on the maximum iteration number ItMax and the current iteration number It. The value of tc can be found by applying the following equation [32].(27)tc=1−ItItMaxβ2It/ItMax,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 efficient solutions through enhancing the exploitation phase and GR can be determined by using the following equation [32].(28)GR=ET.Gq.Fs−l.Fi,where(29)Gq=θ12,ifθ2≥q,0,else.In equation (28), Gq is considered as the control coefficient 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. Therefore, 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].


#### 3.1.3. 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:(30)Finew=Fmin,ifFinew<Fmin,Fmax,ifFinew>Fmax,Finew,else.By calculating fitness function, the quality of each concentration is determined. Solutions with better quality (i.e., concentrations with better mass balance) are retained. Therefore, new concentration qualityFitinew should be compared with the current concentration quality (Fiti) for storing better ones. The explanation can be formulated as follows:(31)Fiti=Fitinew,ifFiti≥Fitinew,Fiti,else,(32)Fi=Finew,ifFiti≥Fitinew,Fi,else.


#### 3.1.4. Improved Equilibrium Optimizer Algorithm (IEO)

By substituting equation (28) into equation (25), new solutions can be obtained:(33)Finew=Fs+Fi−Fs.ET+Fs−l.Fi.Gq.ETl.1−ET.Equation (33) has three terms including Fs, Fi−Fs.ET, and Fs−l.FiGq.ET/l.1−ET in which the second term and third term are the updated jumping steps to make the new solution Finew different from the old solution Fs. In this case, the scaling factor of the first 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. The 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. Thus, in this case, to minimize the randomization that can affect 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. ET and ET′ will be found with respect to l values for determining their operating regions. The operating regions of ET and ET′ are plotted in Figures 2 and 3.Figure 2
The maximum and minimum values of ET.Figure 3
The maximum and minimum values of ET′.Through Figure2, the obtained results clearly indicate that the values of ET completely fluctuate in the range of [−0.6, 0.6], and this implies that the scaling factor of the first jumping step is only suitable for finding 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. Therefore, 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 difficulty 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 finding the optimal solutions in the local search space. In this study, equation (33) has been improved to extend the search area.(34)Finew=FLbest1+Fi−FLbest1.ET+Fr1−Fr2.γ.In the first jumping step, the current solutionFi should be moved to a new location near the location of the best current solution FLbest1 to inherit the best current solution’s experience at the previous generation. This will contribute to increase the probability of finding 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 (Fr1 and Fr2) by taking variables in the Fs solution will bring many advantages in approaching a higher quality solution. As a result, Fr1−Fr2 is always greater than zero and Fr1−Fr2 is multiplied by a randomization vector γ which is comprised of terms between 0 and 1. The results of multiplication facilitate the expansion of the search space to avoid local traps. Generally, equation (33) has advantages in finding solutions in the local space, and equation (34) is useful in expanding the search area. Therefore, the great combination of equations (33) and (34) will improve the performance of EO significantly. In this research, the choice of the update equation depends on the quality of the current solution (Fiti) and the mean fitness of all solutions (Fitmean). If the quality value of Fiti is better than the value of Fitmean, equation (33) is used for generating new solutions. For other cases, equation (34) is applied for new solution update.The implementation process of IEO for a typical problem is shown in Figure4.Figure 4
The implementation process for the whole of IEO for the typical problem.


### 3.2. Application of Improved Equilibrium Optimization

#### 3.2.1. Selection of Decision Variables

In this paper, three DGUs are considered for connecting to IEEE 33-bus, IEEE 69-bus, and IEEE 85-bus radial distribution systems. In this paper, five harmonic flows are simultaneously injected into several linear loads and the harmonic currents affect the entire system. In this study, FW/BWST is chosen to compute power flows at the fundamental frequency and other higher order frequencies. After reaching the power flows, the obtained parameters such as the branch current and the bus voltage at the order frequencies are used to compute the fitness function. Among the input data, the location and capacity of the DGUs are the elements to be determined by using IEO. Thus, each solution (Fi) will consist of location (LDGU,d,i) and capacity (APDGU,d,i) as follows:(35)Fi=LDGU,d,i,APDGU,d,i;d=1,…Nd,gandi=1,…,Nc.In equation (35), LDGU,d,i and APDGU,d,i are the location and capacity of the dth DGU of the ith solution with the limitations as follows:(36)LDGUmin≤LDGU,d,i≤LDGUmax,(37)APDGUmin≤APDGU,d,i≤APDGUmax.The two parameters of DGUs will be randomly generated by(38)LDGU,d,i=roundLDGUmin+randLDGUmax−LDGUmin,(39)APDGU,d,i=APDGUmin+randAPdgumax_APDGUmin.


#### 3.2.2. Penalty Term for Violation

After theith new solution is generated, voltage, current, and the harmonic violations are checked as follows:(40)ΔVb,i=Vmin−Vb,i,ifVb,i<Vmin,Vmax−Vb,i,ifVb,i<Vmax,0,else,(41)ΔIbh,i=Ibhmax−IbhDGU,i,ifIbhDGU,i>Ibhmax,0,else,(42)ΔTHDb,i=THDmax−THDb,iifTHDb,i>THDmax0else,(43)ΔIHDb,ih=IHDmax−IHDb,ihifIHDb,ih>IHDmax0else.


#### 3.2.3. Fitness Function

The fitness function is a combination of objective function and penalty functions. It is used to evaluate the quality of each solution. In this paper, the fitness function of theith solution is constructed as follows:(44)Fiti=OF+σv.∑b=1NbΔVb,i2+σ1.∑bh=1NbhΔVbh,i2+σTHD.∑b=1NbΔTHDb,i2+σ1HD∑b=1NbΔIHDb,ih2.


#### 3.2.4. 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:(45)LDGU,d,i=LDGUmaxifLDGU,d,i>LDGUmaxLDGUminifLDGU,d,i>LDGUminLDGU,d,ielse,(46)APDGU,d,i=APDGUmax,ifAPDGU,d,i>APDGUmax,APDGUmin,ifAPDGU,d,i>APDGUmin,APDGU,d,i,else.


#### 3.2.5. Implementation Process of the Whole IEO for DGU’s Problem

The process of searching location and capacity of DGUs in the distribution systems by using IEO is briefly presented in Figure5, and the implemented steps are as follows:Step 1: survey and select parameters of proposed method.The 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.The 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 fitness function as equation (44). The function value is also quality of each solution.Step 3: start loop in the algorithm.The first iteration is started to perform the promising solution search iterative algorithm.Step 4: determine the group of solutions with high quality.Based on their fitness values in Step 2, the four most effective solutions are determined and assigned toFLbest1, FLbest2, FLbest3, FLbest4 for the calculation of the arithmetic mean (FLmean). The selection of good solutions contributes significantly 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 influence 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. They are important terms of the new solution generation equation.Step 6: check the conditions for selecting the method of new solution generation.Compare the fitness value of each solution (Fiti) with the mean fitness value (Fitmean) to select the equation for creating the next generation Finew.Step 7: create the new solution generation.Either equation (33) or equation (34) is employed to calculate the ith new solution. If Fiti ≥ Fitmean occurs, the ith 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 fitness 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). Then, fitness value in equation (44) is obtained for each new solution.Step 9: find the best current solution.Through the evaluation of the fitness values, the best current solution (FLbest1) and its fitness value (FitLbest1) are retained by applying the rules in equations (31) and (32).Step 10: check the termination condition.The condition for stopping the iteration is checked by comparing It and ItMax. If It = ItMax 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.Figure 5
The implementation process of IEO for placing DGUs in radial distribution networks.


## 3.1. 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. These search agents update their concentration randomly with the desire of finding 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. The quality of the solution entirely depends on the fitness function value under the constraints, and the use of the fitness function result is intended to evaluate the mass-concentration balance within the adjusted volume. Overall, EO is an ideal algorithm for exploration, finding high-quality optimal solutions and avoiding local traps. The algorithm is expressed as follows.


### 3.1.1. Initial Population Generation

In EO, each solutionFi in the population contains a set of adjustment variables and there must be a predetermined range for the variables. The lower and upper bounds of the variables are called the minimum solution (Fmin) and the maximum solution (Fmax), respectively. During the application of EO, each solution Fi is always kept within these limits. The mentioned solutions and the boundary solutions are formulated as follows:(22)Fi=Vaji;j=1,2,3,…,Nad&i=1,2,3,…,Nc,(23)Fmin=Vajmin&Fmax=Vajmax;j=1,…,Nad.The generation of initial solution set or initial population is implemented by [32](24)Fi=Fmin+rand0,1.Fmax−Fmin;i=1,…,Nc.


### 3.1.2. Calculation Process for Updating Concentration

Quality of all solutions is calculated to rank the solutions in population in which solution with the smallest fitness value is the best and other solutions with worse fitness values are from the second-best solution to the worst solution. After obtaining the fitness values of solutions, the top four best solutions are collected and assigned toFLbest1, FLbest2, FLbest3, FLbest4 for the calculation of the arithmetic mean (FLmean). The new solution Finew is updated through a combination of three components. The first component (Fs) is randomly selected from the solution set including FLbest1, FLbest2, FLbest3, FLbest4, and FLmean. 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 difference between the selected solution (Fs) and the considered old solution (Fi). This 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. The three components are useful in finding better solutions, which are not close to the current good solutions. Finally, the new solution Finew can be calculated by applying the following equation [32].(25)Finew=Fs+Fi−FsET+GRl1−ET.In equation (25), two important terms, which mainly affect 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].(26)ET=β1signµ−0.5e−tc.l−1.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 coefficient tc is a value that changes with each iteration, and its value depends entirely on the maximum iteration number ItMax and the current iteration number It. The value of tc can be found by applying the following equation [32].(27)tc=1−ItItMaxβ2It/ItMax,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 efficient solutions through enhancing the exploitation phase and GR can be determined by using the following equation [32].(28)GR=ET.Gq.Fs−l.Fi,where(29)Gq=θ12,ifθ2≥q,0,else.In equation (28), Gq is considered as the control coefficient 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. Therefore, 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].


### 3.1.3. 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:(30)Finew=Fmin,ifFinew<Fmin,Fmax,ifFinew>Fmax,Finew,else.By calculating fitness function, the quality of each concentration is determined. Solutions with better quality (i.e., concentrations with better mass balance) are retained. Therefore, new concentration qualityFitinew should be compared with the current concentration quality (Fiti) for storing better ones. The explanation can be formulated as follows:(31)Fiti=Fitinew,ifFiti≥Fitinew,Fiti,else,(32)Fi=Finew,ifFiti≥Fitinew,Fi,else.


### 3.1.4. Improved Equilibrium Optimizer Algorithm (IEO)

By substituting equation (28) into equation (25), new solutions can be obtained:(33)Finew=Fs+Fi−Fs.ET+Fs−l.Fi.Gq.ETl.1−ET.Equation (33) has three terms including Fs, Fi−Fs.ET, and Fs−l.FiGq.ET/l.1−ET in which the second term and third term are the updated jumping steps to make the new solution Finew different from the old solution Fs. In this case, the scaling factor of the first 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. The 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. Thus, in this case, to minimize the randomization that can affect 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. ET and ET′ will be found with respect to l values for determining their operating regions. The operating regions of ET and ET′ are plotted in Figures 2 and 3.Figure 2
The maximum and minimum values of ET.Figure 3
The maximum and minimum values of ET′.Through Figure2, the obtained results clearly indicate that the values of ET completely fluctuate in the range of [−0.6, 0.6], and this implies that the scaling factor of the first jumping step is only suitable for finding 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. Therefore, 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 difficulty 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 finding the optimal solutions in the local search space. In this study, equation (33) has been improved to extend the search area.(34)Finew=FLbest1+Fi−FLbest1.ET+Fr1−Fr2.γ.In the first jumping step, the current solutionFi should be moved to a new location near the location of the best current solution FLbest1 to inherit the best current solution’s experience at the previous generation. This will contribute to increase the probability of finding 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 (Fr1 and Fr2) by taking variables in the Fs solution will bring many advantages in approaching a higher quality solution. As a result, Fr1−Fr2 is always greater than zero and Fr1−Fr2 is multiplied by a randomization vector γ which is comprised of terms between 0 and 1. The results of multiplication facilitate the expansion of the search space to avoid local traps. Generally, equation (33) has advantages in finding solutions in the local space, and equation (34) is useful in expanding the search area. Therefore, the great combination of equations (33) and (34) will improve the performance of EO significantly. In this research, the choice of the update equation depends on the quality of the current solution (Fiti) and the mean fitness of all solutions (Fitmean). If the quality value of Fiti is better than the value of Fitmean, equation (33) is used for generating new solutions. For other cases, equation (34) is applied for new solution update.The implementation process of IEO for a typical problem is shown in Figure4.Figure 4
The implementation process for the whole of IEO for the typical problem.


## 3.1.1. Initial Population Generation

In EO, each solutionFi in the population contains a set of adjustment variables and there must be a predetermined range for the variables. The lower and upper bounds of the variables are called the minimum solution (Fmin) and the maximum solution (Fmax), respectively. During the application of EO, each solution Fi is always kept within these limits. The mentioned solutions and the boundary solutions are formulated as follows:(22)Fi=Vaji;j=1,2,3,…,Nad&i=1,2,3,…,Nc,(23)Fmin=Vajmin&Fmax=Vajmax;j=1,…,Nad.The generation of initial solution set or initial population is implemented by [32](24)Fi=Fmin+rand0,1.Fmax−Fmin;i=1,…,Nc.


## 3.1.2. Calculation Process for Updating Concentration

Quality of all solutions is calculated to rank the solutions in population in which solution with the smallest fitness value is the best and other solutions with worse fitness values are from the second-best solution to the worst solution. After obtaining the fitness values of solutions, the top four best solutions are collected and assigned toFLbest1, FLbest2, FLbest3, FLbest4 for the calculation of the arithmetic mean (FLmean). The new solution Finew is updated through a combination of three components. The first component (Fs) is randomly selected from the solution set including FLbest1, FLbest2, FLbest3, FLbest4, and FLmean. 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 difference between the selected solution (Fs) and the considered old solution (Fi). This 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. The three components are useful in finding better solutions, which are not close to the current good solutions. Finally, the new solution Finew can be calculated by applying the following equation [32].(25)Finew=Fs+Fi−FsET+GRl1−ET.In equation (25), two important terms, which mainly affect 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].(26)ET=β1signµ−0.5e−tc.l−1.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 coefficient tc is a value that changes with each iteration, and its value depends entirely on the maximum iteration number ItMax and the current iteration number It. The value of tc can be found by applying the following equation [32].(27)tc=1−ItItMaxβ2It/ItMax,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 efficient solutions through enhancing the exploitation phase and GR can be determined by using the following equation [32].(28)GR=ET.Gq.Fs−l.Fi,where(29)Gq=θ12,ifθ2≥q,0,else.In equation (28), Gq is considered as the control coefficient 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. Therefore, 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].


## 3.1.3. 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:(30)Finew=Fmin,ifFinew<Fmin,Fmax,ifFinew>Fmax,Finew,else.By calculating fitness function, the quality of each concentration is determined. Solutions with better quality (i.e., concentrations with better mass balance) are retained. Therefore, new concentration qualityFitinew should be compared with the current concentration quality (Fiti) for storing better ones. The explanation can be formulated as follows:(31)Fiti=Fitinew,ifFiti≥Fitinew,Fiti,else,(32)Fi=Finew,ifFiti≥Fitinew,Fi,else.


## 3.1.4. Improved Equilibrium Optimizer Algorithm (IEO)

By substituting equation (28) into equation (25), new solutions can be obtained:(33)Finew=Fs+Fi−Fs.ET+Fs−l.Fi.Gq.ETl.1−ET.Equation (33) has three terms including Fs, Fi−Fs.ET, and Fs−l.FiGq.ET/l.1−ET in which the second term and third term are the updated jumping steps to make the new solution Finew different from the old solution Fs. In this case, the scaling factor of the first 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. The 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. Thus, in this case, to minimize the randomization that can affect 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. ET and ET′ will be found with respect to l values for determining their operating regions. The operating regions of ET and ET′ are plotted in Figures 2 and 3.Figure 2
The maximum and minimum values of ET.Figure 3
The maximum and minimum values of ET′.Through Figure2, the obtained results clearly indicate that the values of ET completely fluctuate in the range of [−0.6, 0.6], and this implies that the scaling factor of the first jumping step is only suitable for finding 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. Therefore, 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 difficulty 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 finding the optimal solutions in the local search space. In this study, equation (33) has been improved to extend the search area.(34)Finew=FLbest1+Fi−FLbest1.ET+Fr1−Fr2.γ.In the first jumping step, the current solutionFi should be moved to a new location near the location of the best current solution FLbest1 to inherit the best current solution’s experience at the previous generation. This will contribute to increase the probability of finding 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 (Fr1 and Fr2) by taking variables in the Fs solution will bring many advantages in approaching a higher quality solution. As a result, Fr1−Fr2 is always greater than zero and Fr1−Fr2 is multiplied by a randomization vector γ which is comprised of terms between 0 and 1. The results of multiplication facilitate the expansion of the search space to avoid local traps. Generally, equation (33) has advantages in finding solutions in the local space, and equation (34) is useful in expanding the search area. Therefore, the great combination of equations (33) and (34) will improve the performance of EO significantly. In this research, the choice of the update equation depends on the quality of the current solution (Fiti) and the mean fitness of all solutions (Fitmean). If the quality value of Fiti is better than the value of Fitmean, equation (33) is used for generating new solutions. For other cases, equation (34) is applied for new solution update.The implementation process of IEO for a typical problem is shown in Figure4.Figure 4
The implementation process for the whole of IEO for the typical problem.


## 3.2. Application of Improved Equilibrium Optimization

### 3.2.1. Selection of Decision Variables

In this paper, three DGUs are considered for connecting to IEEE 33-bus, IEEE 69-bus, and IEEE 85-bus radial distribution systems. In this paper, five harmonic flows are simultaneously injected into several linear loads and the harmonic currents affect the entire system. In this study, FW/BWST is chosen to compute power flows at the fundamental frequency and other higher order frequencies. After reaching the power flows, the obtained parameters such as the branch current and the bus voltage at the order frequencies are used to compute the fitness function. Among the input data, the location and capacity of the DGUs are the elements to be determined by using IEO. Thus, each solution (Fi) will consist of location (LDGU,d,i) and capacity (APDGU,d,i) as follows:(35)Fi=LDGU,d,i,APDGU,d,i;d=1,…Nd,gandi=1,…,Nc.In equation (35), LDGU,d,i and APDGU,d,i are the location and capacity of the dth DGU of the ith solution with the limitations as follows:(36)LDGUmin≤LDGU,d,i≤LDGUmax,(37)APDGUmin≤APDGU,d,i≤APDGUmax.The two parameters of DGUs will be randomly generated by(38)LDGU,d,i=roundLDGUmin+randLDGUmax−LDGUmin,(39)APDGU,d,i=APDGUmin+randAPdgumax_APDGUmin.


### 3.2.2. Penalty Term for Violation

After theith new solution is generated, voltage, current, and the harmonic violations are checked as follows:(40)ΔVb,i=Vmin−Vb,i,ifVb,i<Vmin,Vmax−Vb,i,ifVb,i<Vmax,0,else,(41)ΔIbh,i=Ibhmax−IbhDGU,i,ifIbhDGU,i>Ibhmax,0,else,(42)ΔTHDb,i=THDmax−THDb,iifTHDb,i>THDmax0else,(43)ΔIHDb,ih=IHDmax−IHDb,ihifIHDb,ih>IHDmax0else.


### 3.2.3. Fitness Function

The fitness function is a combination of objective function and penalty functions. It is used to evaluate the quality of each solution. In this paper, the fitness function of theith solution is constructed as follows:(44)Fiti=OF+σv.∑b=1NbΔVb,i2+σ1.∑bh=1NbhΔVbh,i2+σTHD.∑b=1NbΔTHDb,i2+σ1HD∑b=1NbΔIHDb,ih2.


### 3.2.4. 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:(45)LDGU,d,i=LDGUmaxifLDGU,d,i>LDGUmaxLDGUminifLDGU,d,i>LDGUminLDGU,d,ielse,(46)APDGU,d,i=APDGUmax,ifAPDGU,d,i>APDGUmax,APDGUmin,ifAPDGU,d,i>APDGUmin,APDGU,d,i,else.


### 3.2.5. Implementation Process of the Whole IEO for DGU’s Problem

The process of searching location and capacity of DGUs in the distribution systems by using IEO is briefly presented in Figure5, and the implemented steps are as follows:Step 1: survey and select parameters of proposed method.The 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.The 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 fitness function as equation (44). The function value is also quality of each solution.Step 3: start loop in the algorithm.The first iteration is started to perform the promising solution search iterative algorithm.Step 4: determine the group of solutions with high quality.Based on their fitness values in Step 2, the four most effective solutions are determined and assigned toFLbest1, FLbest2, FLbest3, FLbest4 for the calculation of the arithmetic mean (FLmean). The selection of good solutions contributes significantly 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 influence 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. They are important terms of the new solution generation equation.Step 6: check the conditions for selecting the method of new solution generation.Compare the fitness value of each solution (Fiti) with the mean fitness value (Fitmean) to select the equation for creating the next generation Finew.Step 7: create the new solution generation.Either equation (33) or equation (34) is employed to calculate the ith new solution. If Fiti ≥ Fitmean occurs, the ith 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 fitness 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). Then, fitness value in equation (44) is obtained for each new solution.Step 9: find the best current solution.Through the evaluation of the fitness values, the best current solution (FLbest1) and its fitness value (FitLbest1) are retained by applying the rules in equations (31) and (32).Step 10: check the termination condition.The condition for stopping the iteration is checked by comparing It and ItMax. If It = ItMax 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.Figure 5
The implementation process of IEO for placing DGUs in radial distribution networks.


## 3.2.1. Selection of Decision Variables

In this paper, three DGUs are considered for connecting to IEEE 33-bus, IEEE 69-bus, and IEEE 85-bus radial distribution systems. In this paper, five harmonic flows are simultaneously injected into several linear loads and the harmonic currents affect the entire system. In this study, FW/BWST is chosen to compute power flows at the fundamental frequency and other higher order frequencies. After reaching the power flows, the obtained parameters such as the branch current and the bus voltage at the order frequencies are used to compute the fitness function. Among the input data, the location and capacity of the DGUs are the elements to be determined by using IEO. Thus, each solution (Fi) will consist of location (LDGU,d,i) and capacity (APDGU,d,i) as follows:(35)Fi=LDGU,d,i,APDGU,d,i;d=1,…Nd,gandi=1,…,Nc.In equation (35), LDGU,d,i and APDGU,d,i are the location and capacity of the dth DGU of the ith solution with the limitations as follows:(36)LDGUmin≤LDGU,d,i≤LDGUmax,(37)APDGUmin≤APDGU,d,i≤APDGUmax.The two parameters of DGUs will be randomly generated by(38)LDGU,d,i=roundLDGUmin+randLDGUmax−LDGUmin,(39)APDGU,d,i=APDGUmin+randAPdgumax_APDGUmin.


## 3.2.2. Penalty Term for Violation

After theith new solution is generated, voltage, current, and the harmonic violations are checked as follows:(40)ΔVb,i=Vmin−Vb,i,ifVb,i<Vmin,Vmax−Vb,i,ifVb,i<Vmax,0,else,(41)ΔIbh,i=Ibhmax−IbhDGU,i,ifIbhDGU,i>Ibhmax,0,else,(42)ΔTHDb,i=THDmax−THDb,iifTHDb,i>THDmax0else,(43)ΔIHDb,ih=IHDmax−IHDb,ihifIHDb,ih>IHDmax0else.


## 3.2.3. Fitness Function

The fitness function is a combination of objective function and penalty functions. It is used to evaluate the quality of each solution. In this paper, the fitness function of theith solution is constructed as follows:(44)Fiti=OF+σv.∑b=1NbΔVb,i2+σ1.∑bh=1NbhΔVbh,i2+σTHD.∑b=1NbΔTHDb,i2+σ1HD∑b=1NbΔIHDb,ih2.


## 3.2.4. 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:(45)LDGU,d,i=LDGUmaxifLDGU,d,i>LDGUmaxLDGUminifLDGU,d,i>LDGUminLDGU,d,ielse,(46)APDGU,d,i=APDGUmax,ifAPDGU,d,i>APDGUmax,APDGUmin,ifAPDGU,d,i>APDGUmin,APDGU,d,i,else.


## 3.2.5. Implementation Process of the Whole IEO for DGU’s Problem

The process of searching location and capacity of DGUs in the distribution systems by using IEO is briefly presented in Figure5, and the implemented steps are as follows:Step 1: survey and select parameters of proposed method.The 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.The 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 fitness function as equation (44). The function value is also quality of each solution.Step 3: start loop in the algorithm.The first iteration is started to perform the promising solution search iterative algorithm.Step 4: determine the group of solutions with high quality.Based on their fitness values in Step 2, the four most effective solutions are determined and assigned toFLbest1, FLbest2, FLbest3, FLbest4 for the calculation of the arithmetic mean (FLmean). The selection of good solutions contributes significantly 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 influence 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. They are important terms of the new solution generation equation.Step 6: check the conditions for selecting the method of new solution generation.Compare the fitness value of each solution (Fiti) with the mean fitness value (Fitmean) to select the equation for creating the next generation Finew.Step 7: create the new solution generation.Either equation (33) or equation (34) is employed to calculate the ith new solution. If Fiti ≥ Fitmean occurs, the ith 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 fitness 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). Then, fitness value in equation (44) is obtained for each new solution.Step 9: find the best current solution.Through the evaluation of the fitness values, the best current solution (FLbest1) and its fitness value (FitLbest1) are retained by applying the rules in equations (31) and (32).Step 10: check the termination condition.The condition for stopping the iteration is checked by comparing It and ItMax. If It = ItMax 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.Figure 5
The implementation process of IEO for placing DGUs in radial distribution networks.


## 4. Simulation Results

This 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 profile, 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 coefficient of c1 is the obtained result from function of 2.e−4.It/ItMax2, and the coefficients of c2 and c3 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. Thus, the value of each weighting coefficient depends on the importance of each element in the multi-objective function [31]. Its value can be adjusted by user. The larger value of a weighting coefficient, the greater its significance 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. Thus, the total power loss receives significant 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 five harmonic flows are simultaneously injected to linear loads, and the detailed information is presented in Table 2 [20].Table 2
The detailed information of five harmonic flows.
Harmonic orderMagnitude (%)Angle (°)50.7652870.627−180110.248−59130.12779170.071−253


### 4.1. Case 1: IEEE 33-Bus Radial Distribution System

The configuration of the 33-bus distribution system is shown in Figure6. The load and line data of the system are collected from [31]. The 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 different 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 flows are also injected simultaneously into six buses 9, 15, 20, 24, 29, and 32.Figure 6
IEEE 33-bus radial distribution system.After 50 trial runs are implemented, the worst, average, and best fitness values are collected and clearly presented in Table3. According to the presented results in Table 3, the best fitness value of the implemented methods such as ABC, SSA, EO, and IEO is 0.3866, 0.3837, 0.3832, and 0.3828, respectively. The value of the proposed method (IEO) is lower than that of ABC, SSA, and EO. In other words, IEO can find a more optimal solution than ABC, SSA, and EO. Moreover, the mean and worst fitness values of IEO are smaller than those of others. For a clearer view of performance comparison, fitness 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 fluctuation of IEO is very small in the range of [0.3828, 0.3896]. Thus, IEO is the most stable method.Table 3
The comparison of fitness values found by four implemented methods for case 1.
MethodABCSSAEOIEOThe worst fitness value0.41410.41980.39690.3896The mean fitness value0.39800.39510.38610.3856The best fitness value0.38660.38370.38320.3828Figure 7
The fitness values of 50 trial runs in increasing arrangement for the IEEE 33-bus system.Table4 shows the objective values in the multi-objective function for comparison. The first 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]). The 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. The figure 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. The loss reduction on other remaining branches is not significant, but there are still small values. The very small reduction cannot be identified in the figure due to the constraint of figure size. In general, the installation of DGUs can bring high benefits to the system in reducing power loss.Table 4
The comparison of the objective function values for the IEEE 33-bus system.
MethodOFIOFIIOFReal loss (MW)Cost of DGUs (million $)GA [39]0.50280.49900.50210.10614.6999PSO [39]0.49930.49800.49900.10544.6903GA-PSO [39]0.49000.49800.49160.10344.6901LSF-SA [40]0.38890.41130.39330.08203.8735BFOA [41]0.46590.27840.42840.09832.6234BSOA [42]0.42210.27800.39330.08912.6182ABC0.39410.35650.38660.08323.3572SSA0.38630.37330.38370.08153.5161EO0.38670.36920.38320.08163.4774IEO0.38630.36900.38280.08153.4753Figure 8
Power loss on distribution lines in IEEE 33-bus radial distribution system before and after placing DGUs.The 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 effective 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. This is a trade-off in solving multi-objective function problem. This result is due to the selection of weighting coefficients, and it cannot be avoided when comparing a proposed method to many other methods. However, the proposed IEO has reached better fitness function than the three methods; meanwhile, fitness function is also the multi-objective function and used to evaluate performance of methods. As a result, the proposed IEO is also more effective 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 profile 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 modified system with nonlinear loads and without DGU connection has exceeded the limitations of both THD and IHD as shown in Figure9. Specifically, 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.Figure 9
THD (%) and IHD (%) of the IEEE 33-bus system before and after DGUs connection.Similarly, the voltage profile of the modified system before and after placing three DGUs is also plotted in Figure10. In the system without DGUs, many buses have lower voltage than 0.95 pu and the worst voltage is 0.9038 pu. But in the system with three placed DGUs, the voltage is significantly improved and it has the range from 0.9613 pu to 1.0 pu. Clearly, DGUs can improve the voltage of the system effectively.Figure 10
Voltage profile of the IEEE 33-bus system before and after DGUs connection.Optimal solutions of the applied system are reported in TableS1 in Supplementary Materials.


### 4.2. Case 2: IEEE 69-Bus Radial Distribution System

Figure11 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 flows shown in Table 2 are injected simultaneously into nine buses 12, 18, 19, 22, 25, 34, 46, 56, and 65.Figure 11
The IEEE 69-bus radial distribution system.The worst, mean, and best fitness values for the 50 trial runs are clearly presented in Table5. According to the collected data from Table 5, the best fitness values for ABC, SSA, EO, and IEO are 0.3279, 0.3267, 0.3264, and 0.3262, respectively. The best fitness value of IEO is smaller than that of all other methods. Clearly, the modifications carried out on IEO have a positive impact on the effectiveness 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 fifty trial runs and the fifty runs in detail are also reported in Table 5 and Figure 12. The 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 fifty fitness 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 fluctuations while others have significant fluctuations. In summary, IEO can reach the best performance and the most stable search ability among the four run methods.Table 5
The comparison of fitness values found by four implemented methods for the IEEE 69-bus system.
MethodABCSSAEOIEOThe worst fitness value0.33960.34540.32830.3282The mean fitness value0.33170.32980.32690.3267The best fitness value0.32790.32670.32640.3262Figure 12
The fitness values of 50 trial runs in increasing arrangement for the IEEE 69-bus system.Table6 shows a summary of two objectives and corresponding power loss and cost obtained by implemented and compared methods. The first 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. The 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 significantly different. Specifically, the power loss reduction is greater on branches 4 to 9, 11 to 14, and 52 to 60. The 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 benefits.Table 6
The comparison of the objective function values for the IEEE 69-bus system.
MethodOFIOFIIOFReal loss (MW)Cost of DGUs (million $)MOBA [34]0.32680.33330.32810.07343.1393GA [39]0.39620.49960.41680.08894.7049PSO [39]0.37150.49800.39680.08344.6900GA-PSO [39]0.36120.49800.38860.08114.6901LSF-SA [40]0.34340.36350.34740.07713.4239BFOA [41]0.33500.34800.33760.07523.2776HSA [43]0.38660.29550.36840.08682.7833ABC0.32250.34930.32790.07243.2893SSA0.32160.34720.32670.07233.2699EO0.32070.34900.32640.07203.2868IEO0.32060.34870.32620.071973.2840Figure 13
Power loss on distribution lines in IEEE 69-bus radial distribution system before and after placing DGUs.The second objective and the corresponding cost of IEO are 0.3487 and $3.2840 million, respectively. These 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-off in solving the multi-objective function problem. The results are affected 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. Thus, IEO is actually more effective than others in solving the problem of the location of DGUs in the IEEE 69-bus system.Using the optimal solution obtained by the proposed IEO method, THD and IHD of each bus in the base system without DGUs and modified system with DGUs are plotted in Figure14. Data from the figure indicate that the maximum values of THD and IHD of the base system are 5.4287% and 3.5087%, but those from the modified 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 significant improvement of the modified 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. The minimum bus voltage is 0.9092 pu for the base system but is much higher and equal to 0.9761 pu for the modified system with DGU connection. Voltage values at buses 5 to 28 and 51 to 69 in the modified 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 benefit in reducing harmonic distortions and improving voltage profile.Figure 14
THD (%) and IHD (%) of the IEEE 69-bus system before and after DGUs connection.Figure 15
Voltage profile of the IEEE 69-bus system before and after DGUs connection.Optimal solutions of the applied system are reported in TableS2 in Supplementary Materials.


### 4.3. Case 3: IEEE 85-Bus Radial Distribution System

Figure16 shows the configuration of the IEEE 85-bus radial distribution system. The load data and line data of the system are collected from [31]. The system has the total loss of 0.3161 MW and the total loads of 2.5703 MW and 2.6221 MVar. For reducing loss and cost, location and size of three DGUs are selected from bus No. 2 to bus No. 85 and from 0.0 MW to 2.0 MW. To produce distortion of voltage wave, five harmonic flows in Table 2 are also injected to seven buses 5, 12, 25, 34, 56, 65, and 77.Figure 16
The IEEE 85-bus radial distribution system.The best, mean, and worst fitness values of ABC, SSA, EO, and IEO are presented in Table7. IEO is also still the best method with all smallest values excluding the comparison of the worst fitness value with EO. The fifty runs are also presented in detail in Figure 17 by sorting fitness from the smallest to the highest values shown in different 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 7
The comparison of fitness values found by four implemented methods for the IEEE 85-bus system.
MethodABCSSAEOIEOThe worst fitness value0.53810.56080.48190.4842The mean fitness value0.48600.47330.46330.4614The best fitness value0.45460.44920.44890.4483Figure 17
The fitness values of 50 trial runs in increasing arrangement for the IEEE 85-bus system.Table8 presents the results obtained by the four methods in detail. There are no comparisons with other published methods for the system because there were no studies regarding the IEEE 85-bus system before. The first 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 effective in terms of power loss reduction than others. To see the effectiveness 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. There is outstanding power loss reduction on the first 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 effective as on branch 7, the reduction is relatively high, from some kW to 20 kW. The loss reduction on other branches is not as effective or even there is no loss reduction on some branches. The loss reduction on branches 24 to 29, branch 56, and branch 57 is still high while the reduction on other remaining branches is not identified. In general, power loss reduction is effective in all branches, but the total loss reduction is still high, from 0.3161 MW to 0.1515 MW. The high loss reduction is only for one hour, and it indicates that the placement of DGUs in distribution system is very useful.Table 8
The comparison of the objective function values for the IEEE 85-bus system.
MethodOFIOFIIOFReal loss (MW)Cost of DGUs (million $)ABC0.49120.30830.45460.15532.9036SSA0.48220.31760.44920.15242.9915EO0.48120.31980.44890.15213.0115IEO0.47920.32480.44830.15153.0593Figure 18
Power loss on distribution lines in IEEE 85-bus radial distribution system before and after placing DGUs.The 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-off in solving the multi-objective function problem.By applying the best solution from IEO, THD and IHD of each bus and voltage profile in the base system without DGUs and modified system with DGUs are plotted in Figures19 and 20, respectively. The maximum values of THD and IHD in the base system are 6.1504% and 3.9840%, but those from the modified system with the application of IEO are 4.0445% and 2.6233%, respectively. Approximately all buses in the modified system can reach much better THD and IHD than those in the base system. The 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 profile is improved drastically with the lowest voltage value of 0.9500 pu. The view on voltage of all buses in the two cases can indicate that the modified system can have a much better voltage profile 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.Figure 19
THD (%) and IHD (%) of the IEEE 85-bus system before and after DGUs connection.Figure 20
Voltage profile of the IEEE 85-bus system before and after DGUs connection.Optimal solutions of the applied system are reported in TableS3 in Supplementary Materials.


## 4.1. Case 1: IEEE 33-Bus Radial Distribution System

The configuration of the 33-bus distribution system is shown in Figure6. The load and line data of the system are collected from [31]. The 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 different 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 flows are also injected simultaneously into six buses 9, 15, 20, 24, 29, and 32.Figure 6
IEEE 33-bus radial distribution system.After 50 trial runs are implemented, the worst, average, and best fitness values are collected and clearly presented in Table3. According to the presented results in Table 3, the best fitness value of the implemented methods such as ABC, SSA, EO, and IEO is 0.3866, 0.3837, 0.3832, and 0.3828, respectively. The value of the proposed method (IEO) is lower than that of ABC, SSA, and EO. In other words, IEO can find a more optimal solution than ABC, SSA, and EO. Moreover, the mean and worst fitness values of IEO are smaller than those of others. For a clearer view of performance comparison, fitness 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 fluctuation of IEO is very small in the range of [0.3828, 0.3896]. Thus, IEO is the most stable method.Table 3
The comparison of fitness values found by four implemented methods for case 1.
MethodABCSSAEOIEOThe worst fitness value0.41410.41980.39690.3896The mean fitness value0.39800.39510.38610.3856The best fitness value0.38660.38370.38320.3828Figure 7
The fitness values of 50 trial runs in increasing arrangement for the IEEE 33-bus system.Table4 shows the objective values in the multi-objective function for comparison. The first 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]). The 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. The figure 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. The loss reduction on other remaining branches is not significant, but there are still small values. The very small reduction cannot be identified in the figure due to the constraint of figure size. In general, the installation of DGUs can bring high benefits to the system in reducing power loss.Table 4
The comparison of the objective function values for the IEEE 33-bus system.
MethodOFIOFIIOFReal loss (MW)Cost of DGUs (million $)GA [39]0.50280.49900.50210.10614.6999PSO [39]0.49930.49800.49900.10544.6903GA-PSO [39]0.49000.49800.49160.10344.6901LSF-SA [40]0.38890.41130.39330.08203.8735BFOA [41]0.46590.27840.42840.09832.6234BSOA [42]0.42210.27800.39330.08912.6182ABC0.39410.35650.38660.08323.3572SSA0.38630.37330.38370.08153.5161EO0.38670.36920.38320.08163.4774IEO0.38630.36900.38280.08153.4753Figure 8
Power loss on distribution lines in IEEE 33-bus radial distribution system before and after placing DGUs.The 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 effective 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. This is a trade-off in solving multi-objective function problem. This result is due to the selection of weighting coefficients, and it cannot be avoided when comparing a proposed method to many other methods. However, the proposed IEO has reached better fitness function than the three methods; meanwhile, fitness function is also the multi-objective function and used to evaluate performance of methods. As a result, the proposed IEO is also more effective 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 profile 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 modified system with nonlinear loads and without DGU connection has exceeded the limitations of both THD and IHD as shown in Figure9. Specifically, 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.Figure 9
THD (%) and IHD (%) of the IEEE 33-bus system before and after DGUs connection.Similarly, the voltage profile of the modified system before and after placing three DGUs is also plotted in Figure10. In the system without DGUs, many buses have lower voltage than 0.95 pu and the worst voltage is 0.9038 pu. But in the system with three placed DGUs, the voltage is significantly improved and it has the range from 0.9613 pu to 1.0 pu. Clearly, DGUs can improve the voltage of the system effectively.Figure 10
Voltage profile of the IEEE 33-bus system before and after DGUs connection.Optimal solutions of the applied system are reported in TableS1 in Supplementary Materials.


## 4.2. Case 2: IEEE 69-Bus Radial Distribution System

Figure11 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 flows shown in Table 2 are injected simultaneously into nine buses 12, 18, 19, 22, 25, 34, 46, 56, and 65.Figure 11
The IEEE 69-bus radial distribution system.The worst, mean, and best fitness values for the 50 trial runs are clearly presented in Table5. According to the collected data from Table 5, the best fitness values for ABC, SSA, EO, and IEO are 0.3279, 0.3267, 0.3264, and 0.3262, respectively. The best fitness value of IEO is smaller than that of all other methods. Clearly, the modifications carried out on IEO have a positive impact on the effectiveness 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 fifty trial runs and the fifty runs in detail are also reported in Table 5 and Figure 12. The 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 fifty fitness 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 fluctuations while others have significant fluctuations. In summary, IEO can reach the best performance and the most stable search ability among the four run methods.Table 5
The comparison of fitness values found by four implemented methods for the IEEE 69-bus system.
MethodABCSSAEOIEOThe worst fitness value0.33960.34540.32830.3282The mean fitness value0.33170.32980.32690.3267The best fitness value0.32790.32670.32640.3262Figure 12
The fitness values of 50 trial runs in increasing arrangement for the IEEE 69-bus system.Table6 shows a summary of two objectives and corresponding power loss and cost obtained by implemented and compared methods. The first 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. The 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 significantly different. Specifically, the power loss reduction is greater on branches 4 to 9, 11 to 14, and 52 to 60. The 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 benefits.Table 6
The comparison of the objective function values for the IEEE 69-bus system.
MethodOFIOFIIOFReal loss (MW)Cost of DGUs (million $)MOBA [34]0.32680.33330.32810.07343.1393GA [39]0.39620.49960.41680.08894.7049PSO [39]0.37150.49800.39680.08344.6900GA-PSO [39]0.36120.49800.38860.08114.6901LSF-SA [40]0.34340.36350.34740.07713.4239BFOA [41]0.33500.34800.33760.07523.2776HSA [43]0.38660.29550.36840.08682.7833ABC0.32250.34930.32790.07243.2893SSA0.32160.34720.32670.07233.2699EO0.32070.34900.32640.07203.2868IEO0.32060.34870.32620.071973.2840Figure 13
Power loss on distribution lines in IEEE 69-bus radial distribution system before and after placing DGUs.The second objective and the corresponding cost of IEO are 0.3487 and $3.2840 million, respectively. These 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-off in solving the multi-objective function problem. The results are affected 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. Thus, IEO is actually more effective than others in solving the problem of the location of DGUs in the IEEE 69-bus system.Using the optimal solution obtained by the proposed IEO method, THD and IHD of each bus in the base system without DGUs and modified system with DGUs are plotted in Figure14. Data from the figure indicate that the maximum values of THD and IHD of the base system are 5.4287% and 3.5087%, but those from the modified 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 significant improvement of the modified 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. The minimum bus voltage is 0.9092 pu for the base system but is much higher and equal to 0.9761 pu for the modified system with DGU connection. Voltage values at buses 5 to 28 and 51 to 69 in the modified 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 benefit in reducing harmonic distortions and improving voltage profile.Figure 14
THD (%) and IHD (%) of the IEEE 69-bus system before and after DGUs connection.Figure 15
Voltage profile of the IEEE 69-bus system before and after DGUs connection.Optimal solutions of the applied system are reported in TableS2 in Supplementary Materials.


## 4.3. Case 3: IEEE 85-Bus Radial Distribution System

Figure16 shows the configuration of the IEEE 85-bus radial distribution system. The load data and line data of the system are collected from [31]. The system has the total loss of 0.3161 MW and the total loads of 2.5703 MW and 2.6221 MVar. For reducing loss and cost, location and size of three DGUs are selected from bus No. 2 to bus No. 85 and from 0.0 MW to 2.0 MW. To produce distortion of voltage wave, five harmonic flows in Table 2 are also injected to seven buses 5, 12, 25, 34, 56, 65, and 77.Figure 16
The IEEE 85-bus radial distribution system.The best, mean, and worst fitness values of ABC, SSA, EO, and IEO are presented in Table7. IEO is also still the best method with all smallest values excluding the comparison of the worst fitness value with EO. The fifty runs are also presented in detail in Figure 17 by sorting fitness from the smallest to the highest values shown in different 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 7
The comparison of fitness values found by four implemented methods for the IEEE 85-bus system.
MethodABCSSAEOIEOThe worst fitness value0.53810.56080.48190.4842The mean fitness value0.48600.47330.46330.4614The best fitness value0.45460.44920.44890.4483Figure 17
The fitness values of 50 trial runs in increasing arrangement for the IEEE 85-bus system.Table8 presents the results obtained by the four methods in detail. There are no comparisons with other published methods for the system because there were no studies regarding the IEEE 85-bus system before. The first 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 effective in terms of power loss reduction than others. To see the effectiveness 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. There is outstanding power loss reduction on the first 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 effective as on branch 7, the reduction is relatively high, from some kW to 20 kW. The loss reduction on other branches is not as effective or even there is no loss reduction on some branches. The loss reduction on branches 24 to 29, branch 56, and branch 57 is still high while the reduction on other remaining branches is not identified. In general, power loss reduction is effective in all branches, but the total loss reduction is still high, from 0.3161 MW to 0.1515 MW. The high loss reduction is only for one hour, and it indicates that the placement of DGUs in distribution system is very useful.Table 8
The comparison of the objective function values for the IEEE 85-bus system.
MethodOFIOFIIOFReal loss (MW)Cost of DGUs (million $)ABC0.49120.30830.45460.15532.9036SSA0.48220.31760.44920.15242.9915EO0.48120.31980.44890.15213.0115IEO0.47920.32480.44830.15153.0593Figure 18
Power loss on distribution lines in IEEE 85-bus radial distribution system before and after placing DGUs.The 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-off in solving the multi-objective function problem.By applying the best solution from IEO, THD and IHD of each bus and voltage profile in the base system without DGUs and modified system with DGUs are plotted in Figures19 and 20, respectively. The maximum values of THD and IHD in the base system are 6.1504% and 3.9840%, but those from the modified system with the application of IEO are 4.0445% and 2.6233%, respectively. Approximately all buses in the modified system can reach much better THD and IHD than those in the base system. The 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 profile is improved drastically with the lowest voltage value of 0.9500 pu. The view on voltage of all buses in the two cases can indicate that the modified system can have a much better voltage profile 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.Figure 19
THD (%) and IHD (%) of the IEEE 85-bus system before and after DGUs connection.Figure 20
Voltage profile of the IEEE 85-bus system before and after DGUs connection.Optimal solutions of the applied system are reported in TableS3 in Supplementary Materials.


## 5. 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 Figures7, 12, and 17, IEO can reach many better solutions than EO, ABC, and SSA. Furthermore, IEO also has very small fluctuations, but the fluctuations 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 effective, but it cannot avoid shortcomings such as solving larger and more complicated systems. In fact, IEO is successful in finding 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. The location can be from 2 to 33 for the first 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. This 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. The 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 first system, 230 for the second system, and 250 iterations for the last system. Due to the two limitations, IEO can be ineffective 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. This assumption may not be true in practice. So, in the future, the authors will consider 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 reconfiguration, 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]. Then, 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. These independent variables are checked and penalized if they are beyond predetermined boundaries. The objective function and the penalty terms become fitness function to evaluate the performance of IEO as well as other metaheuristic algorithms. In summary, a higher number of installed DGUs lead to a higher number of considered control variables, and the implementation of IEO for optimally placing separately or simultaneously different types of DGUs is the same.


## 6. Conclusions

This 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 benefits. Overall, the main contributions of the study can be briefly summarized as follows:(i)
The proposed method could reach more promising solutions for the problem than conventional EO. Realistically, the average values of the fitness function of IEO for the three systems were 0.3856, 0.3267, and 0.4614 while these values of the second-best 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). This proves that the improvements in IEO have been remarkably effective.(ii)
For the first 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 first 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 effecitvely. The total cost over 20 years was only $3.4753 million, $3.2840 million, and $3.0593 million for the three systems, respectively.(iii)
The proposed solutions from IEO also satisfied 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 profiles were also significantly enhanced with the voltage range from [0.9038, 1.0] pu to [0.9613, 1.0] pu, from [0.9092, 1.0] pu to [0.9761, 1.0] pu, and from [0.8713, 1.0] pu to [0.95, 1.0] pu for three systems, respectively.Clearly, the contributions of the placement of DGUs are significant for distribution systems. However, DGUs will bring more benefits 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 solar- or wind-based DGUs produce higher power than load demand. The optimization operation of the distribution systems with the presence of DGUs, smart inverters, and BESS problem is an upcoming direction of the study.



---

*Source: 3755754-2022-12-13.xml*