Multi-Objective Optimal Allocation of Electric Vehicle Charging Stations and Distributed Generators in Radial Distribution Systems using Metaheuristic Optimization Algorithms

The acceptance rate of Electric Vehicles (EVs) in the transport industry has increased substantially due to the augmented interest towards sustainable transportation initiatives. However, their impact in terms of increased power demand on the electric power market can increase real power losses, decrease voltage profile, and consequently decrease voltage stability margins. It is necessary to install Electric Vehicle Charging Stations (EVCSs) and Distributed Generators (DGs) at optimal locations to decrease the EV load effect in the Radial Distribution System (RDS). This paper addresses a multiobjective optimization technique to obtain simultaneous EVCS & DG placement and sizing. The problem is formulated to optimize real power losses, Average Voltage Deviation Index (AVDI), and Voltage Stability Index (VSI) of the electrical distribution system. Simulation studies were performed on the standard IEEE 33-bus and 69-bus test systems. Harries Hawk Optimization (HHO) and Teaching-Learning Based Optimization (TLBO) algorithms were selected to minimize the system objectives. The simulation outcomes reveal that the proposed approach improved system performance in all aspects. Among HHO and TLBO, HHO is reasonably successful in accomplishing the desired goals. Keywords-Electrical Vehicles (EVs); Electric Vehicle Charging Stations (EVCSs); renewable energy sources; Distributed Generators (DGs); Average Voltage Deviation Index (AVDI); Voltage Stability Index (VSI); Harris Hawks Optimization (HHO); Teaching-Learning Based Optimization (TLBO)


INTRODUCTION
Ever-increasing demands for energy, the limited availability of fossil fuels, and the global climate change are some major concerns of the 21 st century. The CO 2 emissions from the transportation sector are one of the major causes of global climate change. Researchers highlighted the significant impact of shifting from conventional vehicles to EVs to reduce the transport sector's greenhouse gas contribution. Developing countries are also attracted to EVs due to the low electricity costs as compared to fossil fuel. One of the disadvantages of EVs is the driving range limitation. After running a certain distance, EVs need to recharge the battery at Charging Stations (CSs). Thus, focus must be put on the construction of charging infrastructure for large-scale installation of EVs. Coordinated EVCS planning is of prevailing importance for the sustainable growth of the EV industry.
In [1], a multi-objective function is built that takes into account the performance indexes of the traffic and electrical distribution network, including traffic conditions, charging wait time, power loss, and node voltage, to plan the charging behavior of EVs and achieve optimum efficiency of the entire system. Nevertheless, for optimum location, the effect on the voltage profile and branch power losses needs to be examined. The low sensitivity voltage bus [2] could be the best location for the charging station. Nevertheless, this approach fails to provide the optimum EVCS on the higher radial bus systems. Heuristic algorithms have been developed to optimize the placement of EVCS not only for servicing EVs but also to minimize losses and thus to increase voltage profile. To examine the effect of charging EVs on the electrical distribution network the Fast Charging Station (FCS) load is designed in [3] by taking into account the relationship between the power consumption of EVs and grid voltage. The suggested model is also compared to the constant ZIP model and constant power load model. In [4], the Genetic Algorithm is used mainly to optimally place the EVCSs without considering losses and conductor thermal limits. This method of positioning EVCSs often impacts the reliability of the grid-connected distribution system. By reconfiguring the given RDS [5], the effect of the charging stations could be minimized. In [6], the location of the EVCSs is identified in two stages. The EVCSs' operation spectrum is calculated in the 1 st stage using the trip success ratio with the account of the volatility of trip distance and variability in the remaining electrical charge of the EVs. In the 2 nd stage, the CS spectrum of operation is calculated for optimal CS location. In [7], a tri-objective optimization function that minimizes the charging station's levelized cost of energy, emissions, and losses on the electrical grid by proposing the use of photovoltaic cells and of a battery storage system. Most studies concentrate on the optimum location and size of EVCSs [8][9][10][11][12]. Improper positioning and sizing of EVCSs hurt the distribution system by increasing the network power losses and cause more depreciation in voltage profile. A detailed review of EV technologies and the impact of electrical demand from Plug-in Electric Vehicles (PEVs) on load profiles have been provided in [13]. Authors in [14,15] addressed the problem of dynamic economic dispatch by incorporating PEVs that charge patterns into a 24-hour load demand in economic and environmental dispatch. Authors in [16] suggested a control mechanism to increase the efficiency of the energy shared between a photovoltaic generator and an unbalanced system. In [17], the detrimental impact on the distribution network of resulting EVCS loads is analyzed for voltage stability, reliability indices, power losses, and economic loss of the IEEE 33 bus system. In [18], the simultaneous placement of EVCSs and DGs in the distribution system is proposed in order to reduce EV user charge and Network Power Losses (NPL) cost for the same station development cost and DG investment. Non-dominated sorting Genetic Algorithm II was utilized for achieving the optimal placement of EVCS and DGs. The accuracy of the proposed methodology is checked on the 118 bus system.
Authors in [19] proposed a multi-objective function by considering the daily real power losses and voltage deviations under a 24-hour load pattern, including residential, industrial, and commercial loads. Particle Swarm Optimization (PSO) and Butterfly Optimization (BO) algorithms were applied to minimize the desired objectives. A repeated forward-backward load flow was used to determine losses and voltages. Compared to PSO, BO does well in achieving the desired goals. In [20], a multi-objective hybrid Shuffled Frog Leap Teaching Learning Based Optimization (SFL-TLBO) algorithm was proposed for better planning of EVCSs and DGs in the combined electrical distribution and transportation network, considering the objectives of voltage deviation, power loss, DG power cost, and energy consumption. The proposed model has been assessed on the IEEE 118 bus system and was confirmed to be accurate and stable for different EV population levels. When considering a large number of PEVs charging and discharging during peak and off-peak grid hours, the optimum positioning of EVs on the IEEE-33 radial distribution system using the PSO-CFA optimization algorithm was suggested in [21] to reduce power losses and enhance voltage profile. In [22], the Artificial Bee Colony (ABC) algorithm was proposed for optimal placement of Battery Swapping Stations (BSS) and DG in the RDS to mitigate the power loss increment due to the EV load. If both the DG and BSS are configured concurrently, the system's power losses and reliability are significantly improved. Accordingly, the authors concluded that the simultaneous placement of DG and BSS systems is the best choice to get the maximum benefit from the deployment of DG and BSS in the network. PSO algorithm [23] was used to reduce real power losses and enhance the voltage profile of the system with optimal placement of charging stations. This paper's main contributions are: • Identifying the optimal locations of EVCSs by considering the distribution system losses, AVDI, and VSI of RDS.
• Finding optimum DG size and location by considering optimum EVCS load to boost system bus voltage profile and reducing RDS power losses.
• Placing different categories of DGs and EVCSs simultaneously while taking into account system constraints for minimizing system power losses and better system bus voltage regulation.
• An efficient multi-objective function formed using TLBO [24] and HHO [25] algorithms is proposed for finding the optimal solution of the optimal number of EVCSs, DGs, and their simultaneous positioning and sizing for minimizing distribution system losses, AVDI and enhance VSI of RDS.

A. Objective Function
The size and placement of EVCSs and DGs should be designed to minimize the impact of the increased EV load on distribution system's performance. Radial Distribution System (RDS) has a large R/X ratio and because of this basic load flow tools such as Newton Raphson or fast decoupled approaches do not provide accurate results. An efficient load flow method is presented in [26] based on the forward-backward method for solving the power flows of RDS. The objectives of this study are minimizing power losses, AVDI, and enhancing VSI. The objective function for power losses of the system is given by: where R i is the resistance of the i th branch, I i is the current flowing in the i th branch, i is the branch number and br is the total number of branches. AVDI is defined in terms of all bus voltage magnitudes are deviation w.r.t. reference voltage magnitude 1.0p.u. and is given in: where V k is the voltage at the k th bus, k is the bus number, and b is the total number of buses. The VSI of a receiving end bus of a branch can be determined by [27]: where V k is the voltage at k th bus, P k is the total real power demand at the k th bus, Q k is the total reactive power demand at the k th bus, r jk is the resistance of the branch 'jk', and x jk is the reactance of the branch 'jk'.
VSI k should be greater than 0 for stable RDS operation. The lowest value of VSI can be considered as the overall stability of the given system. The ultimate objective function is developed to minimize power losses, average voltage deviation, and optimize the index of voltage stability. Mathematically: where PG k and QG k are the active and reactive powers injected by DG at the k th bus, P L and Q L are the active and reactive power losses at the k th bus, P d and Q d are the active and reactive power demands at the k th bus, and NG is the total number of buses.
Inequality Constraints: The acceptable limitation of distribution systems should not extend the highest possible allowable power generated from DGs, voltage deviation should be within 5%, the number of charging points (nCP) and the number of Charging Stations (nCSs) should be between minimum and maximum allowable charging points and charging stations respectively, i.e.: ‫ܩܲ‬

PLANNING
Harris Hawks Optimization (HHO) was proposed in 2019. The harris hawk is a predator bird that lives at steady communities found in the southern side of Arizona, United States. Harris hawk's main tactic to catch prey is "surprise pounce" and it can also be recognized as the strategy of "seven kills". In this strategy, many hawks try to strike together from different directions and concentrate on one rabbit. The assault can be easily accomplished by catching the surprised victim in a couple of seconds, although sometimes, with regards to the prey's escape capability and habits, the "seven kills" can include several, short length, rapid dives near the victim for several minutes. Harris hawks can display a diverse range of hunting styles depending on complex nature and prey's escape patterns. When the best hawk (leader) bows down a prey and departs, a switching tactic occurs, and one of the party members can continue the chase. In different circumstances, such switching activities could be examined as they are helpful for annoying the escaping rabbit. The main benefit of such collaborative tactics is that the harris hawks can pursue to exhaustion the detected rabbit which cannot restore its defensive capabilities by confusing the predators and usually one of the hawks, which is often the most effective and skilled, catches the exhausted rabbit quickly and shares it with the others. Generally, HHO is modeled into exploration and exploitation stages inspired by harris hawks exploration of a prey, surprise pounce, and various attacking strategies.

A. Exploration Phase
The position of the hawk is changed in the exploration phase by random location and the other hawks and is given by (12) [25]: where L is the hawk's position, L k is the position of a randomly chosen hawk, L r is the location of prey, t is the present iteration, u lim and l lim are the threshold limits of the search area, d 1 , d 2 , d 3 , d 4 , and r are five discrete random numbers in [0, 1], and L m is the average location of the present hawks' community, which can be calculated by: where L p is the p th hawk in the community, and H is the number of hawks.

B. Transition from Exploration to Exploitation
The HHO algorithm can be transposed from exploration to exploitation. The behavior of exploration is modified depending on the prey's escaping strength. Mathematically, prey's strength to escape can be measured as where T is the max number of iterations, S 0 is the initial strength produced arbitrary in [-1, 1], and d is a random number in [0, 1]. When the prey's escaping strength S is greater than 1, HHO grants the hawks the possibility to hunt on various regions globally. On the other hand, HHO boosts the local search of the best alternatives around the neighborhood when the prey's escaping strength is less than 1.

C. Exploitation Phase
During this phase, the position of the hawk is changed according to the scenarios described below. Its behavior is controlled based on the prey's escape strength (S) and the prey's chance to escape successfully (q<0.5) or not (q>0.5) before a surprise bounce.

1) Soft Besiege
This phase appears when q≥0.5 and |S|≥0.5, and by using (16), the hawk adjusts its position: where S is the prey's escaping strength, L q is the position of prey, ∆L is the difference between the prey's position and present hawk position, and j is the dive strength. The j and ∆L are calculated by [25]: where d 5 is a random number in [0, 1] that is modified in each iteration.

2) Hard Besiege
This phase appears when q≥0.5 and |S|<0.5. In this scenario, the hawk updates its position using [25]:

3) Soft Besiege with Progressive Rapid Dives
This phase occurs when q<0.5 and |S|≥0.5. The hawk is slowly selecting the best possible dive for capturing the prey. In this situation the hawk's new position is created as [25]: where Y and Q are two recently produced hawks, j is the leaping power, N is the number of dimensions, a is an N dimension random vector, and Levy is the Levy flight function that can be computed as: where µ, ϑ are two separate random numbers created from the normal distribution, and σ is described as: where γ is a constant set to 1.5, and the location of the hawk is modified in this phase by: where F is the fitness function and Y and Q are two solutions derived from (20) and (21).

4) Hard Besiege with Progressive Rapid Dives
This is the last scenario and appears when q<0.5 and |S|<0.5. In this condition, the following two new solutions are created [25]: ܻ = ‫ܮ‬ ሺ‫ݐ‬ሻ − ‫ܮ݆|ܵ‬ ሺ‫ݐ‬ሻ − ‫ܮ‬ ሺ‫ݐ‬ሻ| (25) ܳ = ܻ + ܽ * ‫‪ሺܰሻ‬ݕݒ݁ܮ‬ (26) where L m is the average location of the hawks in the current community. The position of the hawk is subsequently modified as:

FUNCTION
In this section, the sequential steps involved in solving the problem of optimal allocation of EVCSs and DGs using the HHO algorithm are given.
Step 1) Define the objective function subjected to equality and inequality constraints of the design variables and read the line and load data.
Step 2) Initialize the number of hawk searches within the upper and lower limits of DG sizes in kW, locations of EVCSs and DGs, HHO parameters and maximum number of iterations T.
Step 4) Using forward-backward load flow, determine the objective function value using (4) for each search hawk.
Step 5) Identify the best hawk (L r ), which has given the best objective function value.
Step 7) If S≥1, update the hawk's position by using (12), otherwise go to Step 8.
Step 8) If S≥0.5, go to Step 9 otherwise go to Step 10.
Step 9) If q≥0.5, update the hawk's position by using (16) otherwise update by using (24) and go to Step 11.
Step 10) If q≥0.5, update the hawk's position by using (20) otherwise update by using (25) and go to Step 11.
Step 11) Increment the iteration count k by 1 and if (number of iterations<maximum number of iterations) then go to Step 4 otherwise go to Step 12.
Step 12) Return the final best solution stored (DG sizes and EVCSs-DGs locations).
Step 13) Run the load flow and determine the power losses, average voltage deviation index, and voltage stability index.

V. RESULTS AND OBSERVATIONS
This section summarizes the simulation results of the 33bus, 69-bus [28] IEEE test systems using the prescribed methods for optimizing the location and capacity of DGs along with the best placement of EVCS. In this paper, EVCSs and DGs locations are categorically studied with three EVCSs and Type I, II, III, and IV DGs [29]. In addition to the EV models (Chevrolet Volt, Chang An Yidong, Tesla Model X and BMW i3) given in Table I, AC/DC level-2 type charging ports (CPs) are also considered in the CS design. According to SAE J1772 standard, this type of CPs can suit both Battery Electric Vehicles (BEVs) and Plug-in Hybrid Electric Vehicles (PHEVs), and has a maximum power rating of 7kW [30]. Some design features of CPs are the EVs types which can charge at a time in a particular CS and their power ratings in kW, the minimum and maximum number of CPs for different types of EVs and correspondingly the minimum and maximum power rating of CS (Table I). For instance, if all the CSs are designed for only minimum CPs as given in Table I, the power rating of 1 CS is 975kW. Similarly, for maximum number CPs, the demand may increase to 1675.5kW. the DGs for minimizing distribution system losses, AVDI, and enhance VSI of RDS. In this article, the following summaries are presented to address the effect of type-I, II, III & IV DGs on the system.
• Scenario 2: RDS with minimum and maximum EVCS load, without integrating DGs.
• Scenario 3: RDS with three optimal locations of EVCSs with minimum EVCS load without integrating DGs.
• Scenario 4: RDS with three optimal locations of EVCSs with minimum EVCS load and Type-I, II, III, and IV DG integration.
• Scenario 5: RDS with three optimal locations of EVCSs with maximum EVCS load without integrating DGs.
• Scenario 6: RDS with three optimal locations of EVCSs with maximum EVCS load and Type-I, II, III, and IV DG integration.

1) 33 Bus System
The cumulative real and reactive demand on the system is 3715kW and 2300kVar. The minimum voltage of 0.9038p.u is observed at bus number 18, which is outside of the prescribed limit of ±5% after executing the distribution power flow algorithm by using the forward-backward method [26]. Moreover, the cumulative active power loss of the system is 210.9897kW, and the cumulative reactive loss of the system is 143.1284kVar. The AVDI of the system is 0.0041 and VSI min is 0.6661. It is recognized that 5.68% of the power is absorbed in the form of active power loss in the system. The regulation of voltage in the system is 9.63%.

2) 69 Bus System
The cumulative real and reactive demand on the system is 3801.4kW & 2693.6kVAr. The maximum voltage deviation occurred at bus 65, and it is the worst among all other busses. The voltage at the 65 th bus is 0.9092p.u., the cumulative active power loss of the system is 224.8807 kW, and the cumulative reactive loss of the system is 102.1647 kVAr. The AVDI of the system is 0.0014, and VSI min is 0.6823. It can be noticed that 5.92% of the power is absorbed in the form of active power loss in the system. The regulation of voltage in the system is 9.08%.

B. Scenario 2 1) 33 Bus System
Incorporation of a total of 3 EVCSs load in the system is assumed. By imposing a total power demand of 2925kW (975kW×3) for 3 CSs with the minimum number of CPs at all EVCSs as given in Table I, the new loading condition of the test system is determined. According to this, the real power load is increased from 3715kW to 6640kW. The corresponding test system performance is given in Table II. Due to the increased EVCS load, the real losses are increased to 576.1705kW from 210.9897kW, the average voltage deviation index is increased to 0.0108 from 0.0041, the stability index is decreased to 0.4984 from 0.6661, and the minimum voltage is decreased to 0.8408p.u. from 0.9038p.u. Similarly, for the maximum number of CPs in each EVCS, the load can increase on the network by 5023.5kW (1674.5×3) for the 3 EVCSs. The corresponding system performance is given in Table II. The losses are increased to 1024.3908kW, the voltage deviation index is raised to 0.0187, the stability index is decreased to 0.3854, and the minimum voltage reaches 0.7888p.u.

2) 69 Bus System
As in the case of the 33-bus test system, 3 EVCSs were considered with different CPs. By imposing a total power demand of 2925kW for 3 CSs with minimum CPs, the test loading condition increased from 3801.4kW to 6726.4kW. Due to the increased EV load, the real losses increased to 613.4994kW from 224.8807kW, the average voltage deviation index increased to 0.004 from 0.0014, the stability index decreased to 0.5114 from 0.6823, and the minimum voltage was reduced to 0.8462p.u. from 0.9092p.u. Similarly, for the maximum number of CPs in each EVCS, the load on the network increased by 5023.5kW (1674.5×3) for the 3 EVCSs. The corresponding system performance is given in Table II. The losses increased to 1108.6266kW, the voltage deviation index increased to 0.0072, the stability index decreased to 0.3949, and the minimum voltage reached 0.7935p.u.

C. Scenario 3
It is assumed that 3 EVCSs are integrated at optimal locations in the test systems with a minimum number of CPs at best locations (bus-2, 19, and 25), see Table III. The real losses are decreased to 295.6474kW from 576.1705kW, the average voltage deviation index decreased to 0.0047 from 0.0108, the voltage stability index improved to 0.6499 from 0.4984, and the minimum voltage increased to 0.8982p.u from 0.8408p.u for the 33 bus system. For the 69 bus system, the best locations are bus-2, 28, and 47. The losses decreased to 225.2186kW from 613.4994kW, AVDI decreased to 0.0014 from 0.004, and VSI increased to 0.6822 from 0.51114. Also, the minimum voltage increased to 0.9092p.u from 0.8462p.u.

VI. CONCLUSIONS
Massive EV deployment has adverse impacts on the power grid. A significant amount of EVCSs connected to the grid increases the system power losses and significant voltage variation at remote buses from sources. In this paper, a novel approach for the simultaneous allocation of DGs and EVCSs to increase bus voltage profile and decrease power losses in the distribution network. In addition to the AC/DC Level-2 EV-CSs suitable for BEVs and PHEVs, different EV models (Chevrolet Volt, Chang An Yidong, Tesla Model X and BMW i3) are taken into account while designing the EVCS with multiple CPs and all four types of DGs. The multi-objective function is formulated and optimized using HHO and TLBO algorithms. The simulation results presented on standard IEEE 33-bus and 69-bus test systems have highlighted the technical benefits which can be achieved through the optimal allocation of EVCSs and DGs even with increased EV loading conditions. The main findings of the simulation are summarized as follows: • HHO algorithm has superior performance when compared to TLBO in assessing optimal size and location of EVCSs and DGs.
• The optimum placement and sizing of the three type-III DGs led to significant reduction in power losses and to voltage enhancement.