# Resiliency Improvement of Distribution Network Considering the Charge/Discharge Management of Electric Vehicles in Parking Lots through a Bilevel Optimization Approach

Authors: Mohammad Alizadeh; Meysam Jafari-Nokandi; Majid Shahabi
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

Due to the growing use of Plug-in Electric vehicles (PEVs) in transportation networks, the charge/discharge scheduling of PEVs in Electric Vehicles Parking Lots (EVPLs) can be effective on the distribution network’s (DN) resiliency. This paper presents a bilevel optimization model to improve the resiliency of the DN taking into account the interaction between the DN islanding problem and the charge/discharge scheduling of PEVs in the energized EVPLs. In the Upper-Level (UL) problem, regarding the electrical loads and managing the charge/discharge of PEVs, the islands’ boundaries are determined with the aim of maximizing the amount of restored load. Knowing the islands’ boundaries and the energized EVPLs from the UL problem, the changes in travels characteristic including destination EVPLs are determined in the Lower-Level (LL) problem to identify the nearest energized EVPL to the out-of-service destination EVPL. The number of PEV drivers that change their deenergized destination depends on the distance between the nearest energized EVPL to the destination. A combination of mathematical programming and evolutionary algorithm is applied to reach the final solution. The proposed model is implemented by applying several concurrent faults to the 118-bus active DN, which is coupled with a 25-node traffic network. The results confirm the efficiency of the proposed model for improving the resiliency of DNs with managing the charge/discharge of PEVs in the restored EVPLs.


---

## Body

## 1. Introduction

Recently, environmental concerns, the security of oil supply, and the increasing penetration of intermittent renewable energy resources in the electrical network are increasing the utilization of PEVs in the transportation sector. On the other hand, the charging pattern of PEVs may affect the resiliency of the DN, which is the ability to withstand the low-probability high-impact events, ensuring the minimum load shedding (LS) and enabling a fast recovery to the normal operation state [1]. Some indicators have been presented in several articles to evaluate the power system resiliency. Reference [2] has introduced some indices to assess the resiliency of the power system using the hierarchical analysis and percolation theory. In [3, 4], by proposing a quantitative framework, the resilience of a MicroGrid (MG) against windstorms has been evaluated using the interrupted load recovery index. A comprehensive review has been presented in [5] concerning the resources for resilience enhancement, the mathematical model of operation and planning algorithms, mathematical formulation, and solution algorithm.Islanding and reconfiguration of the DN are the most important plans for improving the resiliency of the DN. A comprehensive literature review including the methods for improving the power system resiliency during severe events has been conducted in [6]. A bilevel model has been proposed in [7] during a high wind event solved by the bilevel genetic algorithm to improve the distribution network resiliency using an operational network reconfiguration strategy. A three-stage self-healing algorithm has been proposed in [8] to enhance distribution network resiliency by maximizing the restored EVPLs via network reconfiguration without intentional islanding with multiple line faults. An efficient technique for optimum scheduling of MGs with multiperiod islanding restrictions has been proposed in [9] utilizing the Cuttle Fish Algorithm (CFA) and Crow Search Algorithm (CSA) to minimize the MG operation cost, the dispatch-able units operation cost, and the power transmission cost. In [10], chance-constrained stochastic programming has been used to model the resiliency-oriented islanding plan for critical load restoration. The proposed model assumes that the critical loads are provided by only one island, while in practice, it is possible to recover the desired loads by more than one island. A hierarchical outage management structure has been introduced in [11] to ameliorate the resiliency of a DN consisting of several MGs. However, the power system topology and its operational constraints have not been considered. Reference [12] has improved the DN resiliency considering the possibility of several DGs on each island; each of them can control the frequency and voltage. The dynamic islanding plan and optimal management of different smart grid technologies have been investigated in [13] to restore the prioritized loads proposing a mixed-integer linear programming (MILP) multiobjective model. However, the PEVs have been neglected in [13]. A hierarchical and stochastic MILP model has been proposed in [14] for simultaneous modeling of an MG scheduling and virtual power plant (VPP) energy management problems to determine the optimal islands’ boundaries for improving the load restoration and reducing the load supply costs in each MG. However, none of the aforementioned articles has investigated the impact of charging/discharging of PEVs in EVPLs on the DN resiliency.Optimal charging/discharging management of PEVs in EVPLs can offer benefits for the DN. Technically, PEVs can operate in two modes: (1) charging or grid-to-vehicle (G2V) mode as a load, and (2) discharging or vehicle-to-grid (V2G) mode as an energy storage unit [15, 16]. A comprehensive review of control structures of PEVs in charging stations, objectives of EV management in power systems, and optimization methods have been presented in [17]. A detailed comparison between different charging/discharging strategies in terms of complexity, economics, power losses, ability to provide ancillary services for integrity of the power grid, and operation aspects has been done in [18]. Reference [1] has proposed a bilevel model to investigate the impact of charging demand of fast charging stations (FCSs) on DN resiliency. EVPLs are reasonable sites for implementing the V2G because PEVs spend many hours of the day in EVPLs [16]. A multiobjective management model for charging/discharging of PEVs in a smart DN to minimize the total operating costs and emissions has been presented in [19], but the traffic network and the characteristics of the PEV trips have not been considered. Reference [20] has investigated the coordinated charging/discharging strategies of PEVs to smooth load and renewable power fluctuations while ensuring the quality of logistics services. In [21], a stochastic charging/discharging management approach has been proposed for lots of EVs parked in an intelligent EVPL, where intelligent EVPLs have been presented as aggregators allowing PEVs to interact with the electric utilities; however, the traffic network has been neglected in this article. In [22], an optimization-based problem has been proposed to optimally manage the charge/discharge of PEVs in EVPLs for minimizing the operation costs of the DN including the PEV charging cost in EVPLs. Reference [23] has presented a stochastic programming framework for optimal energy management of EVPLs considering demand response plans and uncertain behavior of PEVs. In [24], the traffic and grid-based EVPL, allocation and charging scheduling of PEV fleets have been studied in the planning and operation plans considering the driving pattern of PEV drivers. A two-stage model is presented in [25] for optimal energy management in EVPLs. In the first step, a new scheduling method is proposed for the charge and discharge of PEVs. In the second step, an innovative approach is presented governing EVPL and implementing the encouragement and punishment policies. An energy management model for the EVPL community is presented in [26] for the operational scheduling of several EVPLs, which trade energy with each other besides energy exchange with the power distribution grid. To the best of the authors’ knowledge, none of the above articles has examined the effect of charging/discharging management of PEVs in EVPLs on the resilience of the DN.To solve bilevel models, several methods have been proposed and studied in various articles. If the LL problem has no binary variables, converting the bilevel model to a single-level using the Karush–Kuhn–Tucker (KKT) conditions is the most common approach [27]. But if the LL model includes binary variables, the bilevel models can be solved through iterative calculations between the two levels using different solvers [28, 29, 30]. In [28], the optimal sectionalizing problem in the UL is addressed in CPLEX, while minimizing the outage durations of critical loads in LL is implemented in PSO. In [29], a two-layer algorithm has been solved by relying on a combination of GA and mixed-integer second-order cone programming (calling solver CPLEX). The DIgSILENT and MATLAB are linked together for optimal placement, sizing, and daily charge/discharge the of battery energy storage system in [30].To the best of the authors’ knowledge, an extremely few articles have been introduced on the issue of distribution network resilience in the presence of electric vehicles. However, managing the charge/discharge of PEVs in energized EVPLs may result in better decisions for the DN islanding plan and improving the resiliency. Due to the need for some PEVs to be recharged in EVPLs and to earn the potential profit by selling the energy stored in batteries, some PEV drivers tend to park their PEVs in energized EVPLs. The number of PEV drivers who change their destination to the nearest energized EVPL is inversely related to the distance between the out-of-service destination EVPL and the nearest energized one, so that the number of PEVs decreases with increasing the distance. So, reconnecting some EVPLs to the electrical network can play important role in travel characteristics, identifying the destination EVPLs, and finally, the hourly number/features of PEVs in EVPLs. Therefore, in this paper, a bilevel model is proposed to investigate the interaction between the dynamic features of PEV trips (i.e., the dependency of the trip on the availability of EVPLs) in transportation networks and the DN islanding problem (i.e., the optimal islanding formation to maximize the expected amount of restored load). Briefly, the main contributions of this paper are as follows:(i)
A new MILP model has been proposed to determine the nearest energized EVPL to the out-of-service destination EVPLs considering the traffic network.(ii)
This article examines the effect of dynamic PEV travel characteristics and charge/discharge management of PEVs in energized EVPLs on DN resiliency.(iii)
To investigate the interaction between the distribution and transportation networks, a bilevel optimization model is proposed.The rest of this paper is organized as follows: the bilevel model including the DN resiliency problem and identifying the dynamic trip characteristics is presented in Section2. Simulation results and sensitivity analysis are presented in Section 3, and Section 4 concludes the paper results and future works.


## 2. Bilevel Optimization Approach

Islanding and reconfiguration of DNs following the severe disturbances require the amount of load demands and potential power supplies to be predicted. By increasing the penetration rate of PEVs, the charging/discharging power of PEVs can have a significant impact on the DN islanding and reconfiguring problem. Because most of the PEVs’ daily mileages are less than the PEVs’ drive range, destination charging including home and workplace charging is the main method of energy supply for PEVss [31]. So, managing the charge/discharge of PEVs in places such as EVPLs is vital for improving the performance of the DN.In several studies, the movement of PEVs has been modeled as some fleets, moving from a defined origin to the destination and returning from the same path after doing daily work [32, 33, 34]. PEV fleets in large numbers can be regarded as considerable stochastic loads or energy storage units given the electrical grid [16]. To investigate the impact of charging/discharging management of PEV fleets in EVPLs on DN resiliency, knowledge of the characteristics of the fleet (such as the number of PEVs in each fleet) as well as knowledge of the characteristics of the trip (such as departure time/location and arrival time/location) is essential. Having these characteristics, it is possible to determine the hourly features of PEV fleets in restored EVPLs [33].Given the electrical loads and hourly features of the PEV fleets in EVPLs, the optimal islanding plan is determined to maximize the amount of restored load as a resiliency index. Depending on the DSO’s decision on island formation, some EVPLs may be reconnected to the electrical network. Due to the requirement of PEVs to be recharged in destination EVPLs and also due to the incentives considered by DSO at critical conditions for connecting the PEVs to the electrical network, some PEV drivers tend to change their out-of-service destination EVPL to the nearest in-service one. The possible changes in trip destinations may result in changing the DSO’s decisions for the islanding plan. So, in this paper, the interaction between the DN islanding problem and determining the dynamic trip characteristics of PEV fleets is defined as a bilevel problem. This flowchart is presented in Figure1.Figure 1
The proposed bilevel approach.It is worth mentioning that the location of the restored EVPLs is defined as decision variables in the UL problem, whereas they are considered in the LL problem as parameters. Also, the nearest restored EVPLs to the out-of-service destinations are defined as decision variables in the LL problem, but they are fed to the UL problem as parameters.


### 2.1. Upper-Level Problem

By knowing the hourly number/features of PEVs in each energized EVPL, the DN islanding and the reconfiguring problem is modeled as follows:(1)Obj=max∑ω∈Ω∑t∈T∑i∈Iπω⋅pri,t⋅Pi,t,wload,(2)0≤Pi,t,wload≤αi,t⋅Ploadi,t,(3)Pi,t,ωload−Ploadi,t−Pi,t,ωLC≤1−αi,t⋅M,(4)Pi,t,ωload−Ploadi,t−Pi,t,ωLC≥−1−αi,t⋅M,(5)0≤Qi,t,ωload≤αi,t⋅Qloadi,t,(6)Qi,t,ωload−Qloadi,t−tanφi⋅Pi,t,ωLC≤1−αi,t⋅M,(7)Qi,t,ωload−Qloadi,t−tanφi⋅Pi,t,ωLC≥−1−αi,t⋅M,(8)0≤Pm,t,ωDG≤αi,t⋅PmDG,max,i=NCDGm,(9)αi,t⋅QmDG,min≤Qm,t,ωDG≤αi,t⋅QmDG,max,i=NCDGm,(10)αi,t⋅Vimin≤Vi,t,ω≤αi,t⋅Vimax,(11)−αi,t⋅δimax≤δi,t,w≤αi,t⋅δimax(12)Vi,t,ω=αi,t⋅Viset−point,i∈NCDG,(13)δi,t,ω=0,i∈NCDG,(14)Pm,t,ωDG+Pn,t,ωwind+∑e∈NeviPe,t,ωEV,disch=Pi,tload+∑e∈NeviPe,t,ωEV,ch+∑l∈θLiPl,t,ωflow,∀i,m=NADGi,N=NWTi,(15)Qm,t,ωDG−Qi,tload=∑l∈θLiQl,t,ωflow,m=NADGi,∀l∈θLi,(16)0≤Pn,t,wwind≤αi,t⋅Pn,t,ww,prod,i=NWTn,(17)Pl,t,wflow=Pl,t,wz−flow+F2l×δi,t,w−δj,t,w+F1l×Vi,t,w−Vj,t,w,(18)Ql,t,wflow=Ql,t,wz−flow+F1l×δj,t,w−δi,t,w+F2l×Vi,t,w−Vj,t,w,(19)F1l=rlrl2+xl2,F2l=xlrl2+xl2,(20)−1−βl,t⋅M≤Pl,t,wz−flow≤1−βl,t⋅M,(21)−1−βl,t⋅M≤Ql,t,wz−flow≤1−βl,t⋅M,(22)−βl,t⋅M≤Pl,t,wflow≤βl,t⋅M,(23)−βl,t⋅M≤Ql,t,wflow≤βl,t⋅M,(24)Pl,t,ωflow2+Ql,t,ωflow2≤Slmax2,(25)Ue,i,t,wch+Ue,i,t,wdisch+Ue,i,t,widle=Ue,i,t,(26)−bigM⋅Se,i,telec≤Ue,i,t≤bigM⋅Se,i,telec,(27)Se,i,tbigM≤Ue,i,t≤1+Se,i,telecbigM,(28)Se,i,telec=EVe,b,tnum_elec,(29)Ue,i,t,ωch⋅Pech,min≤Pe,i,t,ωch≤Ue,i,t,ωch⋅Pech,max⋅Se,i,telec,(30)Ue,i,t,ωdisch⋅Pedisch,min≤Pe,i,t,ωdisch≤Ue,i,t,ωdisch⋅Pedisch,max⋅Se,i,telec,(31)Ce⋅Se,i,telec⋅1−DODmax≤SOCe,i,t,welec≤Ce⋅Se,i,telec,(32)SOCe,i,t,w=SOCeinitial⋅Ce⋅EVenum,i=Oe,trip1,t=1,(33)SOCe,i,t,welec=se,t,ielecEVenumSOCe,i,t=τe,trip1exit−pecon×se,i,telec,t=τe,trip1enter,(34)SOCe,i,t,wunelec=sev,t,iunelecEVenumSOCe,i,t=τev,trip1exit−Pecon×se,i,tunelec,t=τev,trip1exit(35)SOCe,t,i,welec=αe,t,i×SOCe,t,i,welec+βe,t,i×SOCe,t,i,wunelec|t=τe,trip2exit−Pecon×Se,t,ielec,t=τe,trip2exit(36)ifse,t=τe,trip2enter,i=De,trip2elec≥se,t=τev,trip2exit,i=Oev,trip2elec,thenαe,t,i=1,βe,t,i=se,t=τe,trip2enter,i=De,trip2elec−se,t=τe,trip2exit,i=Oe,trip2elecse,t=τe,trip2exit,i=Oe,trip2unelec,elecαe,t,i=se,t=τe,trip2enter,i=De,trip2elecse,t=τe,trip2exit,i=Oe,trip2elec,βe,t,i=0,(37)SOCe,i,t+1,welec=SOCe,i,t,welec+Pe,i,t,wch⋅η−Pe,i,t,wdischη.The objective function of the UL problem is defined in (1), to maximize the expected amount of the restored active load in all scenarios, taking into account the priority of loads. The amounts of restored active and reactive loads in each bus are determined by applying load shedding in (2) to (7). According to (6) and (7), it is assumed that the power factor is known and constant for all buses. The active and reactive power generations of dispatchable DGs are limited in (8) and (9). The upper and lower limits of the bus voltage magnitude and angle are stated in (10) and (11), respectively. According to (12) and (13), if a bus that has a DG with voltage and frequency control capability is placed inside an island, the voltage magnitude of that bus will be equal to a predetermined value, and the voltage angle will be zero. The balance between active and reactive power production and consumption in each bus is guaranteed by (14) and (15) considering the charging/discharging of EVs in corresponding EVPLs. As stated in (16), the generating power of a wind farm can be delivered to an island if the relevant bus is located within that island, and also the generating power in each scenario will be less than the maximum predictable generating capacity. The linearized power flow equations are expressed by (17)–(19) [35]. Equations (20) and (21) state the limitations of the auxiliary variable utilized in linearized power flow equations. As stated in (22) and (23), in case of destruction or opening of a line by disaster and DSO, respectively, the active and reactive power flowing that line will be zero. The maximum apparent power flowing each line is specified by (24). Equations (25)–(36), model the charge/discharge scheduling of PEV fleets in restored EVPLs. According to (25), a PEV fleet can be in one of the three modes (i.e., charging, discharging, or idle) while stopping in an energized EVPL. The binary variable U is defined by (26) to (27), which is equal to 1 when some PEVs in fleet e are located in a restored EVPL supplied by electrical bus i. In (28), the number of PEVs in the EVPL located in traffic node b is the same as the number of PEVs in the electrical bus i according to the electrical and traffic nodes connected. The maximum and minimum charging/discharging power of each PEV fleet located in an energized EVPL is limited by the number of PEVs and charger spots according to (29) to (30). The minimum allowable SOC of each PEV fleet in an energized EVPL is limited by the maximum permissible Depth of Discharge (DoD) for batteries in (31). The initial SOC of PEV fleets at the start of the daily trip is specified by (32). The SOC of PEVs when entering an energized/deenergized EVPL at the end of the first trip is determined by (33) and (34), respectively. The SOC of PEVs when entering an energized EVPL at the destination of the second trip is determined by (35), where α and β are the percentages of PEVs that enter a EVPL from energized/deenergized origin EVPLs, respectively. Variables α and β are defined in (36). A graphical explanation of (36) is shown in Figure 2. The hourly SOC of PEVs during the stopping time in an energized EVPL is calculated by (37).Figure 2
Graphical description for parametersα and β in 2nd trip.(38)RI=∑ω∈Ω∑t∈T∑i∈Iπω⋅pri,t⋅Pi,t,wload∑t∈T∑i∈Ipri,t⋅Ploadi,t.To evaluate the efficiency of the proposed bilevel model, the ratio of the expected amount of recovered active loads to the total interrupted active loads considering the loads priority has been introduced as the resilience index (RI) as follows:


#### 2.1.1. Ensuring the Radiality of DN

In this paper, the graph theory is used to ensure the radiality of DN [36]. A DN is radial if successive removal of the first-order nodes (i.e., the nodes connected only to one branch) results in segregation of all nodes. The flowchart of the radiality rule is illustrated in Figure 3.Figure 3
The flowchart of the radiality rule.


### 2.2. Lower-Level Problem

The PEV travels are specified by departure time/location and arrival time/location. Considering the trips characteristic, the hourly features of PEVs such as SOC will be determined. Due to the incentives considered by the DSO, some drivers change their out-of-service destination EVPL to the nearest energized EVPL. Knowing the location of the restored EVPLs specified in the UL problem, the nearest energized EVPLs to the out-of-service destinations are determined in the LL problem that is formulated as follows:(39)Obj=Min∑b,b′∈PLDb,b′,(40)Pb,b′,cDG,tr=Pb,b′,cload,tr+∑s∈ScPb,b′,sflow,tr,∀b∈PL,∀b′∈PL,∀c,(41)Pb,b′,cload,tr=Zb,b′,∀b∈PL,∀b′∈PL,c=b′,(42)Pb,b′,cDG,tr=0,∀b∈PL,∀b′∈PL,∀c≠b,(43)Pb,b′,cDG,tr≥0,∀b∈PL,∀b′∈PL,∀c=b,(44)∑b′∈PLZb,b′=1,∀b∈PL,(45)Zb,b′≤elecb′,t,∀b∈PL,∀b′∈PL,(46)Pb,b′,sflow,trM≤Xb,b′,s≤M⋅Pb,b′,sflow,tr,∀b∈PL,∀b′∈PL,∀s,(47)Db,b′=∑s∈SXb,b′,s⋅adjs,∀b∈PL,∀b′∈PL.Determining the nearest energized EVPL to the destination is specified as the objective function of the LL problem in (39). To ensure the continuity of the path between two EVPLs, the method of forming an equivalent electrical network is used in (40)–(43), in which an energy source is placed in the main destination, and an electric load is placed in the new destination [28, 29]. Performing the load flow equations will guarantee the continuity of the route. Logical constraint (44) refers to the fact that each PEV fleet can only choose one of the energized EVPLs as its destination. According to (45), the drivers only choose the EVPLs as their destinations that have been connected to the electrical network by DSO. It should be noted that the parameter elecb,t that shows the energizing state of a EVPL is determined in the UL problem from αi,t considering the islanding formations and restored electrical buses. According to (46), arcs whose fictitious power flow is non-zero will place on the shortest route between each destination and energized EVPL. The length of the shortest route between the nearest in-service EVPL to each destination EVPL is determined in (47). If the primary destination in traffic node b is energized after islanding, no PEV will change its destination, and the variable Db,b′ will be set to zero for all b′. In this paper, the formation of the islands is assumed to be constant during the fault so α, β, and elec variables are time-independent.


#### 2.2.1. Determining the Hourly Number of PEVs in Destination EVPLs

Because, at the time of occurring the faults, PEVs are stationed in a EVPL, so even in the case of deenergizing that EVPL, no driver will change the EVPL. The hourly number of PEVs in an energized EVPL at the origin of the first trip is obtained by the following equation:(48)EVe,b,tnum−en=EVenum×elecb,t.With identifying the nearest energized EVPL to the out-of-service destination, some drivers change their out-of-service destination EVPL to the nearest energized one considering the minimum distance between EVPLs. This assumption is described in Figure4. The number of PEVs that change their destination decreases as the distance increases. If the distance between the out-of-service destination and the nearest energized EVPL is more than a certain amount (reference distance), none of the drivers change their destinations. The hourly number of PEVs in energized/deenergized EVPLs at the end of the trips is obtained from (49)–(51).(49)EVe,b′,tnum−en=EVenum×Gb,b′,∀e,∀b∈PL,∀b′∈PL,(50)EVe,b,tnum_de=EVenum−∑b′∈PLEVe,b′,tnum_en,∀e,∀b∈PL.Figure 4
Ratio of PEV drivers that change their out-of-service destination to the nearest energized PL.(51)Gb,b′=1−Db,b′DrefDb,b′≤Dref,Gb,b′=0Db,b′≥Dref.whereGb,b′ is obtained by the linear descending function as stated in the following equation [37]:Determining the nearest energized EVPL to the destinations in the LL problem will determine the hourly number and features of PEV fleets in EVPLs. It should be noted that the features of PEV fleets in each EVPL obtained in the LL problem appear as parameters in the UL problem.


### 2.3. Linearization of the Proposed Model

To linearize the proposed bilevel model, the nonlinear limitation of the apparent power flowing each line expressed in (24) is replaced by [38](52)sin2πhH−sin2πh−1HPl,t,wflow−cos2πhH−cos2πh−1HQl,t,wflow≤sin2πHSlmax.The absolute value of the fictitious power flowing each traffic arc for ensuring the continuity of traffic route in (46) is expressed by linear equations as shown by (53)–(55).(53)Pb,b′,sflow,tr=Pb,b′,sflow,tr+−Pb,b′,sflow,tr−,∀b∈PL,∀b′∈PL,∀s,(54)Pb,b′,sflow,tr=Pb,b′′sflow,tr++Pb,b′,sflow,tr−,∀b∈PL,∀b′∈PL,∀s,(55)Pb,b′,sflow,tr+≥0,Pb,b′,sflow,tr−≥0,∀b∈PL,∀b′∈PL,∀s.The following relations have been used to linearize the (51).(56)Db,b′−DrefbigM≤kb,b′≤1+Db,b′−DrefbigM,∀b∈PL,∀b′∈PL(57)−bigM⋅kb,b′≤Gb,b′−1−Db,b′Dref≤bigM⋅kb,b′,∀b∈PL,∀b′∈PL,∀s,(58)−bigM⋅1−kb,b′≤Gb,b′≤bigM⋅1−kb,b′,∀b∈PL,∀b′∈PL,∀s.


### 2.4. Flowchart for Solving the Proposed Bilevel Model

In the proposed bilevel model, the lower-level problem includes binary variables. Therefore, it is not possible to convert the bilevel model to a single-level mathematical programming with equilibrium constraints (MPEC) using the Karush–Kuhn–Tucker (KKT) conditions. Hence, a combination of Genetic Algorithm (GA) and mathematical programming has been utilized to solve the proposed bilevel model [29]. The reason for choosing the GA is the simplicity of use and the ability to search for a large space with high quality. The solving flowchart of the presented bilevel model is illustrated in Figure 5.Figure 5
Bilevel optimization flowchart without the possibility of load control.


## 2.1. Upper-Level Problem

By knowing the hourly number/features of PEVs in each energized EVPL, the DN islanding and the reconfiguring problem is modeled as follows:(1)Obj=max∑ω∈Ω∑t∈T∑i∈Iπω⋅pri,t⋅Pi,t,wload,(2)0≤Pi,t,wload≤αi,t⋅Ploadi,t,(3)Pi,t,ωload−Ploadi,t−Pi,t,ωLC≤1−αi,t⋅M,(4)Pi,t,ωload−Ploadi,t−Pi,t,ωLC≥−1−αi,t⋅M,(5)0≤Qi,t,ωload≤αi,t⋅Qloadi,t,(6)Qi,t,ωload−Qloadi,t−tanφi⋅Pi,t,ωLC≤1−αi,t⋅M,(7)Qi,t,ωload−Qloadi,t−tanφi⋅Pi,t,ωLC≥−1−αi,t⋅M,(8)0≤Pm,t,ωDG≤αi,t⋅PmDG,max,i=NCDGm,(9)αi,t⋅QmDG,min≤Qm,t,ωDG≤αi,t⋅QmDG,max,i=NCDGm,(10)αi,t⋅Vimin≤Vi,t,ω≤αi,t⋅Vimax,(11)−αi,t⋅δimax≤δi,t,w≤αi,t⋅δimax(12)Vi,t,ω=αi,t⋅Viset−point,i∈NCDG,(13)δi,t,ω=0,i∈NCDG,(14)Pm,t,ωDG+Pn,t,ωwind+∑e∈NeviPe,t,ωEV,disch=Pi,tload+∑e∈NeviPe,t,ωEV,ch+∑l∈θLiPl,t,ωflow,∀i,m=NADGi,N=NWTi,(15)Qm,t,ωDG−Qi,tload=∑l∈θLiQl,t,ωflow,m=NADGi,∀l∈θLi,(16)0≤Pn,t,wwind≤αi,t⋅Pn,t,ww,prod,i=NWTn,(17)Pl,t,wflow=Pl,t,wz−flow+F2l×δi,t,w−δj,t,w+F1l×Vi,t,w−Vj,t,w,(18)Ql,t,wflow=Ql,t,wz−flow+F1l×δj,t,w−δi,t,w+F2l×Vi,t,w−Vj,t,w,(19)F1l=rlrl2+xl2,F2l=xlrl2+xl2,(20)−1−βl,t⋅M≤Pl,t,wz−flow≤1−βl,t⋅M,(21)−1−βl,t⋅M≤Ql,t,wz−flow≤1−βl,t⋅M,(22)−βl,t⋅M≤Pl,t,wflow≤βl,t⋅M,(23)−βl,t⋅M≤Ql,t,wflow≤βl,t⋅M,(24)Pl,t,ωflow2+Ql,t,ωflow2≤Slmax2,(25)Ue,i,t,wch+Ue,i,t,wdisch+Ue,i,t,widle=Ue,i,t,(26)−bigM⋅Se,i,telec≤Ue,i,t≤bigM⋅Se,i,telec,(27)Se,i,tbigM≤Ue,i,t≤1+Se,i,telecbigM,(28)Se,i,telec=EVe,b,tnum_elec,(29)Ue,i,t,ωch⋅Pech,min≤Pe,i,t,ωch≤Ue,i,t,ωch⋅Pech,max⋅Se,i,telec,(30)Ue,i,t,ωdisch⋅Pedisch,min≤Pe,i,t,ωdisch≤Ue,i,t,ωdisch⋅Pedisch,max⋅Se,i,telec,(31)Ce⋅Se,i,telec⋅1−DODmax≤SOCe,i,t,welec≤Ce⋅Se,i,telec,(32)SOCe,i,t,w=SOCeinitial⋅Ce⋅EVenum,i=Oe,trip1,t=1,(33)SOCe,i,t,welec=se,t,ielecEVenumSOCe,i,t=τe,trip1exit−pecon×se,i,telec,t=τe,trip1enter,(34)SOCe,i,t,wunelec=sev,t,iunelecEVenumSOCe,i,t=τev,trip1exit−Pecon×se,i,tunelec,t=τev,trip1exit(35)SOCe,t,i,welec=αe,t,i×SOCe,t,i,welec+βe,t,i×SOCe,t,i,wunelec|t=τe,trip2exit−Pecon×Se,t,ielec,t=τe,trip2exit(36)ifse,t=τe,trip2enter,i=De,trip2elec≥se,t=τev,trip2exit,i=Oev,trip2elec,thenαe,t,i=1,βe,t,i=se,t=τe,trip2enter,i=De,trip2elec−se,t=τe,trip2exit,i=Oe,trip2elecse,t=τe,trip2exit,i=Oe,trip2unelec,elecαe,t,i=se,t=τe,trip2enter,i=De,trip2elecse,t=τe,trip2exit,i=Oe,trip2elec,βe,t,i=0,(37)SOCe,i,t+1,welec=SOCe,i,t,welec+Pe,i,t,wch⋅η−Pe,i,t,wdischη.The objective function of the UL problem is defined in (1), to maximize the expected amount of the restored active load in all scenarios, taking into account the priority of loads. The amounts of restored active and reactive loads in each bus are determined by applying load shedding in (2) to (7). According to (6) and (7), it is assumed that the power factor is known and constant for all buses. The active and reactive power generations of dispatchable DGs are limited in (8) and (9). The upper and lower limits of the bus voltage magnitude and angle are stated in (10) and (11), respectively. According to (12) and (13), if a bus that has a DG with voltage and frequency control capability is placed inside an island, the voltage magnitude of that bus will be equal to a predetermined value, and the voltage angle will be zero. The balance between active and reactive power production and consumption in each bus is guaranteed by (14) and (15) considering the charging/discharging of EVs in corresponding EVPLs. As stated in (16), the generating power of a wind farm can be delivered to an island if the relevant bus is located within that island, and also the generating power in each scenario will be less than the maximum predictable generating capacity. The linearized power flow equations are expressed by (17)–(19) [35]. Equations (20) and (21) state the limitations of the auxiliary variable utilized in linearized power flow equations. As stated in (22) and (23), in case of destruction or opening of a line by disaster and DSO, respectively, the active and reactive power flowing that line will be zero. The maximum apparent power flowing each line is specified by (24). Equations (25)–(36), model the charge/discharge scheduling of PEV fleets in restored EVPLs. According to (25), a PEV fleet can be in one of the three modes (i.e., charging, discharging, or idle) while stopping in an energized EVPL. The binary variable U is defined by (26) to (27), which is equal to 1 when some PEVs in fleet e are located in a restored EVPL supplied by electrical bus i. In (28), the number of PEVs in the EVPL located in traffic node b is the same as the number of PEVs in the electrical bus i according to the electrical and traffic nodes connected. The maximum and minimum charging/discharging power of each PEV fleet located in an energized EVPL is limited by the number of PEVs and charger spots according to (29) to (30). The minimum allowable SOC of each PEV fleet in an energized EVPL is limited by the maximum permissible Depth of Discharge (DoD) for batteries in (31). The initial SOC of PEV fleets at the start of the daily trip is specified by (32). The SOC of PEVs when entering an energized/deenergized EVPL at the end of the first trip is determined by (33) and (34), respectively. The SOC of PEVs when entering an energized EVPL at the destination of the second trip is determined by (35), where α and β are the percentages of PEVs that enter a EVPL from energized/deenergized origin EVPLs, respectively. Variables α and β are defined in (36). A graphical explanation of (36) is shown in Figure 2. The hourly SOC of PEVs during the stopping time in an energized EVPL is calculated by (37).Figure 2
Graphical description for parametersα and β in 2nd trip.(38)RI=∑ω∈Ω∑t∈T∑i∈Iπω⋅pri,t⋅Pi,t,wload∑t∈T∑i∈Ipri,t⋅Ploadi,t.To evaluate the efficiency of the proposed bilevel model, the ratio of the expected amount of recovered active loads to the total interrupted active loads considering the loads priority has been introduced as the resilience index (RI) as follows:


### 2.1.1. Ensuring the Radiality of DN

In this paper, the graph theory is used to ensure the radiality of DN [36]. A DN is radial if successive removal of the first-order nodes (i.e., the nodes connected only to one branch) results in segregation of all nodes. The flowchart of the radiality rule is illustrated in Figure 3.Figure 3
The flowchart of the radiality rule.


## 2.1.1. Ensuring the Radiality of DN

In this paper, the graph theory is used to ensure the radiality of DN [36]. A DN is radial if successive removal of the first-order nodes (i.e., the nodes connected only to one branch) results in segregation of all nodes. The flowchart of the radiality rule is illustrated in Figure 3.Figure 3
The flowchart of the radiality rule.


## 2.2. Lower-Level Problem

The PEV travels are specified by departure time/location and arrival time/location. Considering the trips characteristic, the hourly features of PEVs such as SOC will be determined. Due to the incentives considered by the DSO, some drivers change their out-of-service destination EVPL to the nearest energized EVPL. Knowing the location of the restored EVPLs specified in the UL problem, the nearest energized EVPLs to the out-of-service destinations are determined in the LL problem that is formulated as follows:(39)Obj=Min∑b,b′∈PLDb,b′,(40)Pb,b′,cDG,tr=Pb,b′,cload,tr+∑s∈ScPb,b′,sflow,tr,∀b∈PL,∀b′∈PL,∀c,(41)Pb,b′,cload,tr=Zb,b′,∀b∈PL,∀b′∈PL,c=b′,(42)Pb,b′,cDG,tr=0,∀b∈PL,∀b′∈PL,∀c≠b,(43)Pb,b′,cDG,tr≥0,∀b∈PL,∀b′∈PL,∀c=b,(44)∑b′∈PLZb,b′=1,∀b∈PL,(45)Zb,b′≤elecb′,t,∀b∈PL,∀b′∈PL,(46)Pb,b′,sflow,trM≤Xb,b′,s≤M⋅Pb,b′,sflow,tr,∀b∈PL,∀b′∈PL,∀s,(47)Db,b′=∑s∈SXb,b′,s⋅adjs,∀b∈PL,∀b′∈PL.Determining the nearest energized EVPL to the destination is specified as the objective function of the LL problem in (39). To ensure the continuity of the path between two EVPLs, the method of forming an equivalent electrical network is used in (40)–(43), in which an energy source is placed in the main destination, and an electric load is placed in the new destination [28, 29]. Performing the load flow equations will guarantee the continuity of the route. Logical constraint (44) refers to the fact that each PEV fleet can only choose one of the energized EVPLs as its destination. According to (45), the drivers only choose the EVPLs as their destinations that have been connected to the electrical network by DSO. It should be noted that the parameter elecb,t that shows the energizing state of a EVPL is determined in the UL problem from αi,t considering the islanding formations and restored electrical buses. According to (46), arcs whose fictitious power flow is non-zero will place on the shortest route between each destination and energized EVPL. The length of the shortest route between the nearest in-service EVPL to each destination EVPL is determined in (47). If the primary destination in traffic node b is energized after islanding, no PEV will change its destination, and the variable Db,b′ will be set to zero for all b′. In this paper, the formation of the islands is assumed to be constant during the fault so α, β, and elec variables are time-independent.


### 2.2.1. Determining the Hourly Number of PEVs in Destination EVPLs

Because, at the time of occurring the faults, PEVs are stationed in a EVPL, so even in the case of deenergizing that EVPL, no driver will change the EVPL. The hourly number of PEVs in an energized EVPL at the origin of the first trip is obtained by the following equation:(48)EVe,b,tnum−en=EVenum×elecb,t.With identifying the nearest energized EVPL to the out-of-service destination, some drivers change their out-of-service destination EVPL to the nearest energized one considering the minimum distance between EVPLs. This assumption is described in Figure4. The number of PEVs that change their destination decreases as the distance increases. If the distance between the out-of-service destination and the nearest energized EVPL is more than a certain amount (reference distance), none of the drivers change their destinations. The hourly number of PEVs in energized/deenergized EVPLs at the end of the trips is obtained from (49)–(51).(49)EVe,b′,tnum−en=EVenum×Gb,b′,∀e,∀b∈PL,∀b′∈PL,(50)EVe,b,tnum_de=EVenum−∑b′∈PLEVe,b′,tnum_en,∀e,∀b∈PL.Figure 4
Ratio of PEV drivers that change their out-of-service destination to the nearest energized PL.(51)Gb,b′=1−Db,b′DrefDb,b′≤Dref,Gb,b′=0Db,b′≥Dref.whereGb,b′ is obtained by the linear descending function as stated in the following equation [37]:Determining the nearest energized EVPL to the destinations in the LL problem will determine the hourly number and features of PEV fleets in EVPLs. It should be noted that the features of PEV fleets in each EVPL obtained in the LL problem appear as parameters in the UL problem.


## 2.2.1. Determining the Hourly Number of PEVs in Destination EVPLs

Because, at the time of occurring the faults, PEVs are stationed in a EVPL, so even in the case of deenergizing that EVPL, no driver will change the EVPL. The hourly number of PEVs in an energized EVPL at the origin of the first trip is obtained by the following equation:(48)EVe,b,tnum−en=EVenum×elecb,t.With identifying the nearest energized EVPL to the out-of-service destination, some drivers change their out-of-service destination EVPL to the nearest energized one considering the minimum distance between EVPLs. This assumption is described in Figure4. The number of PEVs that change their destination decreases as the distance increases. If the distance between the out-of-service destination and the nearest energized EVPL is more than a certain amount (reference distance), none of the drivers change their destinations. The hourly number of PEVs in energized/deenergized EVPLs at the end of the trips is obtained from (49)–(51).(49)EVe,b′,tnum−en=EVenum×Gb,b′,∀e,∀b∈PL,∀b′∈PL,(50)EVe,b,tnum_de=EVenum−∑b′∈PLEVe,b′,tnum_en,∀e,∀b∈PL.Figure 4
Ratio of PEV drivers that change their out-of-service destination to the nearest energized PL.(51)Gb,b′=1−Db,b′DrefDb,b′≤Dref,Gb,b′=0Db,b′≥Dref.whereGb,b′ is obtained by the linear descending function as stated in the following equation [37]:Determining the nearest energized EVPL to the destinations in the LL problem will determine the hourly number and features of PEV fleets in EVPLs. It should be noted that the features of PEV fleets in each EVPL obtained in the LL problem appear as parameters in the UL problem.


## 2.3. Linearization of the Proposed Model

To linearize the proposed bilevel model, the nonlinear limitation of the apparent power flowing each line expressed in (24) is replaced by [38](52)sin2πhH−sin2πh−1HPl,t,wflow−cos2πhH−cos2πh−1HQl,t,wflow≤sin2πHSlmax.The absolute value of the fictitious power flowing each traffic arc for ensuring the continuity of traffic route in (46) is expressed by linear equations as shown by (53)–(55).(53)Pb,b′,sflow,tr=Pb,b′,sflow,tr+−Pb,b′,sflow,tr−,∀b∈PL,∀b′∈PL,∀s,(54)Pb,b′,sflow,tr=Pb,b′′sflow,tr++Pb,b′,sflow,tr−,∀b∈PL,∀b′∈PL,∀s,(55)Pb,b′,sflow,tr+≥0,Pb,b′,sflow,tr−≥0,∀b∈PL,∀b′∈PL,∀s.The following relations have been used to linearize the (51).(56)Db,b′−DrefbigM≤kb,b′≤1+Db,b′−DrefbigM,∀b∈PL,∀b′∈PL(57)−bigM⋅kb,b′≤Gb,b′−1−Db,b′Dref≤bigM⋅kb,b′,∀b∈PL,∀b′∈PL,∀s,(58)−bigM⋅1−kb,b′≤Gb,b′≤bigM⋅1−kb,b′,∀b∈PL,∀b′∈PL,∀s.


## 2.4. Flowchart for Solving the Proposed Bilevel Model

In the proposed bilevel model, the lower-level problem includes binary variables. Therefore, it is not possible to convert the bilevel model to a single-level mathematical programming with equilibrium constraints (MPEC) using the Karush–Kuhn–Tucker (KKT) conditions. Hence, a combination of Genetic Algorithm (GA) and mathematical programming has been utilized to solve the proposed bilevel model [29]. The reason for choosing the GA is the simplicity of use and the ability to search for a large space with high quality. The solving flowchart of the presented bilevel model is illustrated in Figure 5.Figure 5
Bilevel optimization flowchart without the possibility of load control.


## 3. Case Studies

The proposed model is run on 11 kV DN shown in Figure6. This network consists of 3 feeders, 118 buses, 3 breakers, 30 sectionalizers, and 9 tie lines [39]. The peak active and reactive loads of the DN are 22.71 MW and 17.04 Mvar, respectively. The priority of all loads is considered equal. The details of the DN such as electrical loads and line parameters are given in [39]. The required data about the location and capacity of the DGs are shown in Table 1. This network consists of 14 DGs, 7 of which (DG1 to DG7) can control the voltage and frequency. The formation of an island depends on the presence of a DG with the frequency and voltages control ability. It means that, on each island, there should be at least one DG with such ability. So, in this study network, it is possible to form a maximum of 7 islands with the existence of DG1 to DG7. There are also 10 WTs operating at the unity power factor. The power generating capacity of each WT is equal to 500 kW, and their locations are shown in Figure 6.Figure 6
118 bus distribution network.Table 1
DGs capacity and location data.
unitBusPmax (kW)Qmax (kVAr)unitBusPmax (kW)Qmax (kVAr)DG1171000600DG87200150DG2241000800DG933300200DG351900800DG1043300200DG45912001000DG1188200100DG5671000800DG12103500300DG6761100800DG13113300150DG710715001200DG14117300150The uncertainty in load predicting is neglected. The total amount of active and reactive loads of DN during 24 hours is 436,026 kWh and 327,188 kVARh, respectively. In this paper, it is assumed that the electrical loads cannot be individually decreased or removed, but parts of adjacent loads can be disconnected from DN by opening some sectionalizers or tie switches. Scenarios corresponding to the uncertainty of the wind power and the daily load profile are illustrated in Figures7 and 8, respectively [12].Figure 7
Wind power scenarios (p.u.).Figure 8
Daily load profile (p.u.).The transportation network consisting of 25 nodes is represented in Figure9 [31]. The number next to each arc expresses the distance between nodes in miles. Also, the overlap of the electrical and traffic nodes containing EVPLs is shown in Table 2. The model proposed in [32, 33, 34] is used to model the movement of PEVs between EVPLs on the transportation network. Three PEV fleets with three different routes move from a specific origin to the destination, and again after doing their daily work, they return to the same origin from the same route. For specifying the arrival/departure times and arrival/departure EVPLs as shown in Table 3, it is assumed that fleets 1 and 2 will travel from the residential area to the industrial area, and fleet 3 will travel from the residential area to the commercial area. The traffic routes of each fleet are shown in Figure 6. Table 4 specifies the number of PEVs in each fleet and their features. The charge/discharge power of each charger in EVPLs is 7.3/6.2 kW, respectively. The maximum battery DOD, PEV energy consumption, and the charge/discharge efficiency are considered equal to 80%, 0.27 kWh/mile, and 85%, respectively. The reference distance (Dref) is assumed to be 30 miles for both trips.Figure 9
The 25-node transportation network.Table 2
Overlap of traffic and electrical nodes containing PLs.
PL number123456789101112Electrical node11152261334866428692101108Traffic node2578913141517202325Table 3
PEV fleet features.
FleetSOC at origin (%)Trip length (miles)EV numberCapacity (kWh)110022350192901330027.43801925027.4Table 4
Travel characteristics.
FleetFirst tripSecond tripDepartureArrivalDepartureArrivalhourBushourBushourBushourBus1611742174218112733892189219333152216101191012022In this paper, to validate the proposed model, three different cases are considered as follows.(i)
Case 1: neglecting the PEVs.(ii)
Case 2: considering the PEVs in EVPLs, but PEV drivers do not tend to change their destinations.(iii)
Case 3: considering the charge/discharge management of PEVs in EVPLs, some PEV drivers tend to change their out-of-service destination EVPL to the nearest energized one.Subsequently, a sensitivity analysis is performed to evaluate the impact of reference distance on PEV drivers’ behaviors and travel characteristics and finally on DN resiliency.In this paper, the proposed bilevel model is solved utilizing GAMS and MATLAB simultaneously. For an appropriate setting of the GA parameters, different values have been examined. In this paper, by choosing the GA population size equal to 50 or even higher, all populations converge to the same answer, so the value of 50 has been chosen for the GA population size to reduce the computation time. The amount of crossover and mutation rates has been considered equal to 0.8 and 0.5, respectively. The simulation results are analyzed in the following:


### 3.1. Case 1

In this case, the effect of charging/discharging management of PEV fleets on DN resiliency is not considered. In this case, the optimal islanding plan along with specifying the boundaries of the islands is illustrated in Figure10. Also, the characteristics of each island including restored buses, the expected value of restored active/ reactive loads, and the DGs on each island are listed in Table 5.Figure 10
Optimal islanding plan in case 1.Table 5
The features of six formed islands in case 1.
IslandRestored busesRestored active loads (MWh)Restored reactive loads (MVARh)Energy resources111–168,7274,245DG1, WT125–9, 22–2718,2678,674DG2, DG8346–5316,2609,509DG3, WT4441–43, 59–6214,4467,247DG4, WT3, DG105————675–77, 90, 97–9919,9678,731DG6, WT9794–96, 106–10916,12616,125DG7, WT8Total restored load107, 07654,533As shown in Figure10, six islands are formed in case 1 where no islands were created in the presence of DG5 due to the lack of enough generating capacity and a large amount of loads on the smallest possible island that can be created by opening the sectionalizers S19 and S20. The imbalance between production and consumption active power on this unformed island is shown in Figure 11 where the peak of active loads on this island is only 6.2 kW more than the generation capacity of DG5 at 15 o’clock, so this island is not formed.Figure 11
Active power imbalance in island with DG5 in case 1.By solving the proposed model, the amount of 107,076 kWh of interrupted active loads is restored. Therefore, according to (38), the resiliency index is 24.5%. Also, 54,533 kVARh of the reactive loads is restored.


### 3.2. Case 5

In this case, the impact of charging/ discharging management of PEVs in EVPLs is considered, but PEV drivers do not tend to change their destination EVPLs even if the destination EVPLs have not been connected to the electrical network and therefore are out of service. The boundaries and optimal islanding formation in case 2 are obtained the same as in case 1. The reason is that, due to the movement of PEV fleets between EVPLs, every fleet is temporarily in one of the EVPLs; therefore, because of the fixed islanding formation, none of the fleets can serve as a permanent energy source for restoring more loads. Some islands may need temporary energy resources at specific hours (for example, island 5 at 15 o’clock), but in this study network, there is no PEV fleet inside island 5. In this situation, no more load will be added to the formed islands relative to case 1, and by restoring the 107,076 kWh of interrupted loads, only 24.5% of the interrupted loads is recovered. So, the movement of PEVs between EVPLs and charge/discharge management of PEVs have no positive effect on improving the DN resiliency in this study network.


### 3.3. Case 3

In this case, the effect of dynamic travel characteristics considering the tendency of PEV drivers to change their out-of-service destination EVPL to the nearest energized one is analyzed. The optimal islanding formation is shown in Figure12. Unlike cases 1 and 2, where 6 islands were obtained for DN islanding, in this case, 7 islands are designated as the optimal number. Totally, the amount of 126,405 kWh and 63,663 kVARh of active and reactive disconnected load is restored during 24 hours, respectively. So, according to (38), the resiliency index will be equal to 29%, which increased by 4.5% compared to cases 1 and 2. The proposed model has reached the optimal answer of 126,405 kWh after 20 repetitions. The convergence rate of the GA for the proposed model is shown in Figure 13.Figure 12
Optimal islanding formation in case 3.Figure 13
Convergence rate of the genetic algorithm for the proposed model in case 3.Once the optimal islanding is formed by the DSO, the EVPLs located in the traffic nodes 2, 5, 7, 8, 13, 14, 15, and 25 are inserted into the formed islands and become energized, so these EVPLs can provide charging/discharging service to PEVs. Also, the EVPLs located in the traffic nodes 9, 17, 20, and 23 are not restored by DSO and so are out of service. As mentioned in (50), due to the lack of restoration of some EVPLs by DSO, some PEV drivers change their out-of-service destination EVPLs to the nearest in-service EVPL considering the minimum distance between EVPLs. Table 6 shows the energized/deenergized states of the destination EVPLs, the nearest energized EVPL to the out-of-service destination, the distances between EVPLs, and the ratio of PEVs that change their out-of-service destination. Due to the changes in some destinations, the characteristics of the travels will change as shown in Table 7. By comparing Tables 3 and 7, it is deduced that the destination of the first trip for fleets 2 and 3 and the destination of the second trip for fleet 2 are changed. In Table 7, the changed destinations have been highlighted for easier realization.Table 6
Changing the out-of-service destinations PLs in case 3.
DestinationService state of the PLThe nearest energized PLMinimum distance between PLs (mile)GTraffic (electrical) nodeTraffic (electrical) nodes2 (11)In-service2 (11)007 (22)In-service7 (22)009 (33)Out-of-service8 (61)60.815 (42)In-service15 (42)0020 (92)Out-of-service14 (66)40.8723 (101)Out-of-service14 (66)70.77Table 7
Travel characteristics in case 3.
FleetFirst tripSecond tripDepartureArrivalDepartureArrivalhourElectrical bushourElectrical bushourElectrical bushourElectrical bus16117421742181127338661866196131522166619662022Due to the changes in the trip characteristics for fleets 2 and 3, the number of PEVs that are willing to change their out-of-service destination to the nearest energized EVPL is shown in Figure14. The number written next to each line indicates the number of PEVs moving on the route. Energized and deenergized EVPLs are marked with circles and rectangles around the node number, respectively.Figure 14
Travel characteristics considering the restored PLs (a) fleet 2; (b) fleet3.
(a)(b)PEVs in fleet 2 travel from the electric bus 33 to 92 on their first trip, where none of the origin/destination EVPLs are restored by DSO. According to the aforementioned assumption, none of the PEVs change the origin out-of-service EVPL connected to electrical bus 33 on the first trip, but some PEV drivers change their out-of-service destination EVPL connected to bus 92 to the closest energized EVPL that is connected to the electrical bus 66. According to the information in Table6 and considering the distance of 4 miles between these two EVPLs, 87% of the PEVs in fleet 2 (i.e., 260 PEVs) change their destination from electrical bus 92 to 66, but the others do not change their destination. On the second trip, where destination bus 33 is not connected to the electrical network, the closest energized EVPL to the destination is located on bus 61, which is due to the 6-mile distance between these two EVPLs, and 80% of PEVs (i.e., 240 PEVs) change their out-of-service destination from bus 33 to 61. Because the number of these PEVs is less than the number of PEVs that were in energized EVPL at bus 66 at the origin of the second trip, therefore, 92% of these PEVs change their destination again to electrical bus 61 (α = 0.92).PEVs in fleet 3 also move from electrical bus 22 to 101 on the first trip. The destination of the first trip on bus 101 is deenergized. According to the information in Table6, the nearest energized EVPL to the out-of-service destination EVPL at the end of the first trip is the EVPL located at bus 66. Considering the distance of 7 miles between these two EVPLs, 77% of PEVs change their destination from bus 101 to 66 (i.e., 192 PEVs). On the second trip, where destination EVPL at bus 22 is connected to the electrical network, no PEV changes its destination. So, all 250 PEVs (58 PEVs from bus 101 and 192 PEVs from bus 66) go to EVPL connected to electrical bus 22, so the variables α and β are equal to 1.The SOC of PEV fleets in case 3 is shown in Figures15 to 17. As shown in Figure 16, fleet 2 has lost part of its SOC at the first trip destination EVPL connected to bus 66 to supply the load on island 5 (bus 66–70) in overloading hours. So, island 5 is formed by managing the charge/discharge of PEV fleet 2 in the EVPL7 connected to electrical bus 66.Figure 15
SOC of PEV fleet 1 in case 3 (the first scenario).Figure 16
SOC of PEV fleet 2 in case 3 (the first scenario).Figure 17
SOC of PEV fleet 3 in case 3 (the first scenario).As shown in Figure11, without considering the PEVs, the overload of 6.2 kW at hour 15 prevents the formation of island 5. It means that only a discharging of 6.2 kWh by PEV fleet 2 in destination EVPL located in electrical bus 66 at hour 15 was enough to form island 5. The reason for the overdischarging of PEV fleet 2 in hours of 15 to 16 is that operating costs of DGs, WTs, and depletion cost of PEV batteries are not considered in the objective function of the proposed bilevel model. The charge/discharge scheduling of PEV fleets is shown in Figures 18–20, respectively. It should be noted that, in Figures 18–20, the negative and positive values indicate the charging/discharging power of PEV fleets in electrical buses connected to the corresponding EVPLs, respectively. To compare the results of three different simulated cases better, the summary of the obtained results is given in Table 8.Figure 18
Charge/discharge power of PEV fleet 1 in case 3 (the first scenario).Figure 19
Charge/discharge power of PEV fleet 2 in case 3 (the first scenario).Figure 20
Charge/discharge power of PEV fleet 3 in case 3 (the first scenario).Table 8
Comparing the results for three different cases.
Charging/discharging of PEVsTrips characteristicsRestored active load (kWh)Restored reactive load (kVARh)Resiliency index (%)Number of formed islandsCase 1✗✗107,07654,53324.56Case 2✓Static107,07654,53324.56Case 3✓Dynamic126,40563,663297


### 3.4. Sensitivity Analysis

A sensitivity analysis is performed concerning the changes in the reference distances on the PEV drivers’ behavior and DN resiliency, which is shown in Figure21.Figure 21
Sensitivity analysis of objective function with respect to different reference distance (Dref).As is clear from the results, by setting a reference distance of fewer than 4 miles for PEV drivers to change their out-of-service destination EVPL, according to (4), no driver changes their out-of-service destination. Therefore, by not forming island 5 consisting of electrical buses 66–70, only 107,076 kWh of active loads are recovered by DSO. By setting the reference distance to 4.02 miles, for PEVs in fleet 2 that the minimum distance between the out-of-service destinations to the nearest energized EVPL is 4 miles, 0.49% of PEVs change their deenergized destination to EVPL connected to electrical bus 66. For PEVs in fleet 3 that the minimum distance between the out-of-service destination to the nearest energized EVPL is 7 miles, no driver is willing to change destinations to the nearest energized EVPL. Due to the power shortage of 6.2 kW to form island 5, by setting the reference distance to 4.02 miles, the ratio of 0.49% of PEVs in fleet 2 is sufficient to compensate for this shortage; therefore, the expected amount of restored load increases to 126,405 kWh. Further increase in reference distance does not have much effect on objective function because these PEVs are in a specified EVPL for a short time and cannot provide the energy for a permanent state.


## 3.1. Case 1

In this case, the effect of charging/discharging management of PEV fleets on DN resiliency is not considered. In this case, the optimal islanding plan along with specifying the boundaries of the islands is illustrated in Figure10. Also, the characteristics of each island including restored buses, the expected value of restored active/ reactive loads, and the DGs on each island are listed in Table 5.Figure 10
Optimal islanding plan in case 1.Table 5
The features of six formed islands in case 1.
IslandRestored busesRestored active loads (MWh)Restored reactive loads (MVARh)Energy resources111–168,7274,245DG1, WT125–9, 22–2718,2678,674DG2, DG8346–5316,2609,509DG3, WT4441–43, 59–6214,4467,247DG4, WT3, DG105————675–77, 90, 97–9919,9678,731DG6, WT9794–96, 106–10916,12616,125DG7, WT8Total restored load107, 07654,533As shown in Figure10, six islands are formed in case 1 where no islands were created in the presence of DG5 due to the lack of enough generating capacity and a large amount of loads on the smallest possible island that can be created by opening the sectionalizers S19 and S20. The imbalance between production and consumption active power on this unformed island is shown in Figure 11 where the peak of active loads on this island is only 6.2 kW more than the generation capacity of DG5 at 15 o’clock, so this island is not formed.Figure 11
Active power imbalance in island with DG5 in case 1.By solving the proposed model, the amount of 107,076 kWh of interrupted active loads is restored. Therefore, according to (38), the resiliency index is 24.5%. Also, 54,533 kVARh of the reactive loads is restored.


## 3.2. Case 5

In this case, the impact of charging/ discharging management of PEVs in EVPLs is considered, but PEV drivers do not tend to change their destination EVPLs even if the destination EVPLs have not been connected to the electrical network and therefore are out of service. The boundaries and optimal islanding formation in case 2 are obtained the same as in case 1. The reason is that, due to the movement of PEV fleets between EVPLs, every fleet is temporarily in one of the EVPLs; therefore, because of the fixed islanding formation, none of the fleets can serve as a permanent energy source for restoring more loads. Some islands may need temporary energy resources at specific hours (for example, island 5 at 15 o’clock), but in this study network, there is no PEV fleet inside island 5. In this situation, no more load will be added to the formed islands relative to case 1, and by restoring the 107,076 kWh of interrupted loads, only 24.5% of the interrupted loads is recovered. So, the movement of PEVs between EVPLs and charge/discharge management of PEVs have no positive effect on improving the DN resiliency in this study network.


## 3.3. Case 3

In this case, the effect of dynamic travel characteristics considering the tendency of PEV drivers to change their out-of-service destination EVPL to the nearest energized one is analyzed. The optimal islanding formation is shown in Figure12. Unlike cases 1 and 2, where 6 islands were obtained for DN islanding, in this case, 7 islands are designated as the optimal number. Totally, the amount of 126,405 kWh and 63,663 kVARh of active and reactive disconnected load is restored during 24 hours, respectively. So, according to (38), the resiliency index will be equal to 29%, which increased by 4.5% compared to cases 1 and 2. The proposed model has reached the optimal answer of 126,405 kWh after 20 repetitions. The convergence rate of the GA for the proposed model is shown in Figure 13.Figure 12
Optimal islanding formation in case 3.Figure 13
Convergence rate of the genetic algorithm for the proposed model in case 3.Once the optimal islanding is formed by the DSO, the EVPLs located in the traffic nodes 2, 5, 7, 8, 13, 14, 15, and 25 are inserted into the formed islands and become energized, so these EVPLs can provide charging/discharging service to PEVs. Also, the EVPLs located in the traffic nodes 9, 17, 20, and 23 are not restored by DSO and so are out of service. As mentioned in (50), due to the lack of restoration of some EVPLs by DSO, some PEV drivers change their out-of-service destination EVPLs to the nearest in-service EVPL considering the minimum distance between EVPLs. Table 6 shows the energized/deenergized states of the destination EVPLs, the nearest energized EVPL to the out-of-service destination, the distances between EVPLs, and the ratio of PEVs that change their out-of-service destination. Due to the changes in some destinations, the characteristics of the travels will change as shown in Table 7. By comparing Tables 3 and 7, it is deduced that the destination of the first trip for fleets 2 and 3 and the destination of the second trip for fleet 2 are changed. In Table 7, the changed destinations have been highlighted for easier realization.Table 6
Changing the out-of-service destinations PLs in case 3.
DestinationService state of the PLThe nearest energized PLMinimum distance between PLs (mile)GTraffic (electrical) nodeTraffic (electrical) nodes2 (11)In-service2 (11)007 (22)In-service7 (22)009 (33)Out-of-service8 (61)60.815 (42)In-service15 (42)0020 (92)Out-of-service14 (66)40.8723 (101)Out-of-service14 (66)70.77Table 7
Travel characteristics in case 3.
FleetFirst tripSecond tripDepartureArrivalDepartureArrivalhourElectrical bushourElectrical bushourElectrical bushourElectrical bus16117421742181127338661866196131522166619662022Due to the changes in the trip characteristics for fleets 2 and 3, the number of PEVs that are willing to change their out-of-service destination to the nearest energized EVPL is shown in Figure14. The number written next to each line indicates the number of PEVs moving on the route. Energized and deenergized EVPLs are marked with circles and rectangles around the node number, respectively.Figure 14
Travel characteristics considering the restored PLs (a) fleet 2; (b) fleet3.
(a)(b)PEVs in fleet 2 travel from the electric bus 33 to 92 on their first trip, where none of the origin/destination EVPLs are restored by DSO. According to the aforementioned assumption, none of the PEVs change the origin out-of-service EVPL connected to electrical bus 33 on the first trip, but some PEV drivers change their out-of-service destination EVPL connected to bus 92 to the closest energized EVPL that is connected to the electrical bus 66. According to the information in Table6 and considering the distance of 4 miles between these two EVPLs, 87% of the PEVs in fleet 2 (i.e., 260 PEVs) change their destination from electrical bus 92 to 66, but the others do not change their destination. On the second trip, where destination bus 33 is not connected to the electrical network, the closest energized EVPL to the destination is located on bus 61, which is due to the 6-mile distance between these two EVPLs, and 80% of PEVs (i.e., 240 PEVs) change their out-of-service destination from bus 33 to 61. Because the number of these PEVs is less than the number of PEVs that were in energized EVPL at bus 66 at the origin of the second trip, therefore, 92% of these PEVs change their destination again to electrical bus 61 (α = 0.92).PEVs in fleet 3 also move from electrical bus 22 to 101 on the first trip. The destination of the first trip on bus 101 is deenergized. According to the information in Table6, the nearest energized EVPL to the out-of-service destination EVPL at the end of the first trip is the EVPL located at bus 66. Considering the distance of 7 miles between these two EVPLs, 77% of PEVs change their destination from bus 101 to 66 (i.e., 192 PEVs). On the second trip, where destination EVPL at bus 22 is connected to the electrical network, no PEV changes its destination. So, all 250 PEVs (58 PEVs from bus 101 and 192 PEVs from bus 66) go to EVPL connected to electrical bus 22, so the variables α and β are equal to 1.The SOC of PEV fleets in case 3 is shown in Figures15 to 17. As shown in Figure 16, fleet 2 has lost part of its SOC at the first trip destination EVPL connected to bus 66 to supply the load on island 5 (bus 66–70) in overloading hours. So, island 5 is formed by managing the charge/discharge of PEV fleet 2 in the EVPL7 connected to electrical bus 66.Figure 15
SOC of PEV fleet 1 in case 3 (the first scenario).Figure 16
SOC of PEV fleet 2 in case 3 (the first scenario).Figure 17
SOC of PEV fleet 3 in case 3 (the first scenario).As shown in Figure11, without considering the PEVs, the overload of 6.2 kW at hour 15 prevents the formation of island 5. It means that only a discharging of 6.2 kWh by PEV fleet 2 in destination EVPL located in electrical bus 66 at hour 15 was enough to form island 5. The reason for the overdischarging of PEV fleet 2 in hours of 15 to 16 is that operating costs of DGs, WTs, and depletion cost of PEV batteries are not considered in the objective function of the proposed bilevel model. The charge/discharge scheduling of PEV fleets is shown in Figures 18–20, respectively. It should be noted that, in Figures 18–20, the negative and positive values indicate the charging/discharging power of PEV fleets in electrical buses connected to the corresponding EVPLs, respectively. To compare the results of three different simulated cases better, the summary of the obtained results is given in Table 8.Figure 18
Charge/discharge power of PEV fleet 1 in case 3 (the first scenario).Figure 19
Charge/discharge power of PEV fleet 2 in case 3 (the first scenario).Figure 20
Charge/discharge power of PEV fleet 3 in case 3 (the first scenario).Table 8
Comparing the results for three different cases.
Charging/discharging of PEVsTrips characteristicsRestored active load (kWh)Restored reactive load (kVARh)Resiliency index (%)Number of formed islandsCase 1✗✗107,07654,53324.56Case 2✓Static107,07654,53324.56Case 3✓Dynamic126,40563,663297


## 3.4. Sensitivity Analysis

A sensitivity analysis is performed concerning the changes in the reference distances on the PEV drivers’ behavior and DN resiliency, which is shown in Figure21.Figure 21
Sensitivity analysis of objective function with respect to different reference distance (Dref).As is clear from the results, by setting a reference distance of fewer than 4 miles for PEV drivers to change their out-of-service destination EVPL, according to (4), no driver changes their out-of-service destination. Therefore, by not forming island 5 consisting of electrical buses 66–70, only 107,076 kWh of active loads are recovered by DSO. By setting the reference distance to 4.02 miles, for PEVs in fleet 2 that the minimum distance between the out-of-service destinations to the nearest energized EVPL is 4 miles, 0.49% of PEVs change their deenergized destination to EVPL connected to electrical bus 66. For PEVs in fleet 3 that the minimum distance between the out-of-service destination to the nearest energized EVPL is 7 miles, no driver is willing to change destinations to the nearest energized EVPL. Due to the power shortage of 6.2 kW to form island 5, by setting the reference distance to 4.02 miles, the ratio of 0.49% of PEVs in fleet 2 is sufficient to compensate for this shortage; therefore, the expected amount of restored load increases to 126,405 kWh. Further increase in reference distance does not have much effect on objective function because these PEVs are in a specified EVPL for a short time and cannot provide the energy for a permanent state.


## 4. Conclusion and Future Research

The coordinated interaction of islanding scheme and charging management of PEV fleets in electric vehicle parking lots can maximize the amount of restored load and improve the resiliency of the DN. The results of the proposed model verified that it can determine the optimal islanding plan considering the dynamic travel features. The efficiency of the proposed bilevel model was investigated in different cases. According to the results, only 24.5% of the DN active loads are restored, without the presence of PEVs. However, if PEV drivers change their out-of-service destination EVPL to the nearest energized EVPL, the resiliency of distribution network can be improved by 4.5%. In addition, increasing the reference distance increases the motivation of PEV owners to change their destination EVPL and improves network resilience. Therefore, providing incentives for PEV owners in such conditions can be beneficial for both distribution network operator and PEV owners.In this paper, the objective function contains only the expected amount of restored load and does not include any cost term, so it is not possible to decide on the optimal output power of dispatchable/renewable DGs and the charging/discharging power of PEVs to supply the loads during islanding. Therefore, in the future work, the cost of energy can be added to the objective function of the islanding problem to determine the optimal share of each of the energy resources in formed islands.



---

*Source: 3878440-2022-09-23.xml*