# A Survey on the Electric Vehicle Routing Problem: Variants and Solution Approaches

Authors: Tomislav Erdelić; Tonči Carić
Journal: Journal of Advanced Transportation
 (2019)

---

## Abstract

In order to ensure high-quality and on-time delivery in logistic distribution processes, it is necessary to efficiently manage the delivery fleet. Nowadays, due to the new policies and regulations related to greenhouse gas emission in the transport sector, logistic companies are paying higher penalties for each emission gram ofCO2/km. With electric vehicle market penetration, many companies are evaluating the integration of electric vehicles in their fleet, as they do not have local greenhouse gas emissions, produce minimal noise, and are independent of the fluctuating oil price. The well-researched vehicle routing problem (VRP) is extended to the electric vehicle routing problem (E-VRP), which takes into account specific characteristics of electric vehicles. In this paper, a literature review on recent developments regarding the E-VRP is presented. The challenges that emerged with the integration of electric vehicles in the delivery processes are described, together with electric vehicle characteristics and recent energy consumption models. Several variants of the E-VRP and related problems are observed. To cope with the new routing challenges in E-VRP, efficient VRP heuristics and metaheuristics had to be adapted. An overview of the state-of-the-art procedures for solving the E-VRP and related problems is presented.


---

## Body

## 1. Introduction

The vehicle routing problem (VRP) is an NP-hard optimization problem that aims to determine a set of least-cost delivery routes from a depot to a set of geographically scattered customers, subject to side constraints [1]. The problem was first defined by Dantzig and Ramser [2] as the Truck Dispatching Problem. VRP is a generalization of the well-known traveling salesman problem (TSP), which aims to design one least-cost route to visit all the customers. The problem has applications in several real-life optimization problems, which has led to the definition of many problem variants over the years: limited vehicle load capacity (capacitated VRP, CVRP), customer time windows (VRP with time windows, VRPTW), multiple depots (multidepot VRP, MDVRP), pickup and delivery (VRP with pickup and delivery, VRPPD), time-dependent travel time (time-dependent VRP, TD-VRP), heterogeneous fleet (mixed fleet VRP, MFVRP), etc. [3, 4]. Due to the complexity of the problem, exact procedures are only capable of optimally solving small-sized problems: up to 360 customers for CVRP [5] and 50-100 customers for VRPTW [6]. Over the years, a vast number of heuristics, metaheuristics, and hybrid procedures were proposed for solving different VRP problems.In the past decade, the European Union (EU) has announced many new actions and regulations related to greenhouse gas (GHG) emissions in the transport sector [7]. External factors and the rise of social and ecological awareness have prompted green initiatives in many companies. Conventional internal combustion engine vehicles (ICEVs), which are dependent on limited fossil fuels, severely pollute the environment, especially in congested urban areas. According to WEEA [8], the EU intends to decrease GHG emissions by 20% and 40% by 2020 and 2030, respectively. Sbihi and Eglese [9] introduced the research field of green logistics, which deals with the sustainability of delivery processes by taking into account environmental and social factors. With the electric vehicle (EV) market penetration, many logistic companies evaluated the use of the EVs in their vehicle fleet in order to decrease GHG emissions and, therefore, reduce the charges for every emission gram of CO2/km. EVs have several advantages compared to ICEVs: (i) they do not have local GHG emissions; (ii) they produce minimal noise; (iii) they can be powered from renewable energy sources; and (iv) they are independent of the fluctuating oil price [10, 11]. There are two basic configurations of EVs: the battery electric vehicle (BEV), which is exclusively powered from batteries mounted inside the vehicle, and the hybrid electric vehicle (HEV), which can be powered from batteries inside the vehicle or by other energy sources, most commonly internal combustion engine. The plug-in HEV (PHEV) can be recharged by connecting the plug to the electric power source. In this paper mostly BEVs for logistic purposes are considered, where two main problems come to the fore: limited driving range and the need for additional recharging infrastructure.Due to the limited battery capacity, the range that delivery BEVs can achieve with a fully charged battery is 160-240 km [12], which is much lower than the 480-650 km range of ICEVs [13]. To achieve a similar driving range as ICEVs, BEVs have to visit charging stations (CSs) more frequently. Today, there is still a lack of CSs in the road network infrastructure, and their locations and energy demand should be planned in future infrastructural plans. For an empty BEV to become operable again, battery energy has to be renewed at a CS. This can be performed in two ways: (i) by swapping empty batteries with fully charged ones at Battery Swapping Station (BSS) or (ii) by charging at CS [14, 15]. The former process can be performed in a time comparable to the refueling time of ICEVs. In the latter process BEVs recharge their batteries at CSs by plugging into the electric power source. The recharge time depends on the state of charge (SoC) when entering the CS, the desired SoC level when leaving CS, and the charging function.


### 1.1. Recent Literature Reviews and Scientific Contribution

Here we present, to our knowledge, the most recent literature reviews on the E-VRP and related problems.Juan et al. [16] presented a review regarding the environmental, strategic, and operational challenges of EV integration in logistics and transportation activities. The authors performed a comprehensive analysis of environmental challenges including the transportation impact on the pollution and EVs’ possible contribution to the reduction of carbon emissions. Regarding the strategic challenges, the authors presented issues related to CSs: battery swap technology, different charging technologies, CS location problem, limited number of charges, and charging network; issues related to the EVs: mixed fleet of ICEVs and EVs, economic challenges when integrating EVs in a fleet, and routing constraints. The solution approaches for the E-VRP were presented in general, as a class of VRP solving procedures. Similarly, Margaritis et al. [17] presented practical and research challenges of EVs focusing on battery development, lack of charger compatibility, systematic energy management, the lack of optimization procedures that could minimize the EV routing and scheduling decisions, cooling/heating usage, financial sources, novel policies, measures for EV deployment in transport services, etc.Montoya [18] researched several variants of the E-VRP: green VRP (GVRP), E-VRP with partial recharging and nonlinear charging functions, and the technician routing problem with a mixed fleet of ICEVs and EVs. For each problem, effective solving procedures were proposed: multispace sampling heuristic, iterated local search enhanced with heuristic concentration, and two-phase parallel metaheuristic based on solving a set of subproblems and extended set-covering formulation. The authors also formulated a fixed route vehicle charging problem (FRVCP) with and without time windows to optimize the charging decisions for a route with a fixed customer sequence. The authors did not focus as much on reviewing the E-VRP literature, especially not on the procedures for solving the problem. Compared to Montoya [18], this paper did not focus as much on the FRVCP and technician routing problem or on detailed procedures used to solve those problems.The most recent survey on the E-VRP is presented by Pelletier et al. [19] and it includes technical background on EV types and batteries, EV market penetration, EV competitiveness and incentives, and an overview of the existing research regarding EVs in transportation science. The authors provided a comprehensive review on, at the time, the latest solving procedures for the E-VRP with mixed fleet and optimal paths and covered several key papers regarding partial recharges, hybrid vehicles, and different charging technologies.In this paper, a survey on the E-VRP is presented, which includes approaches for solving the E-VRP and related problems that emerged with BEVs integration in the logistic processes. The focus is not on the economic and environmental challenges related to BEVs. Most of the latest literature reviews were published in 2016; hence, to the best of our knowledge, there is no published research that summarizes the state-of-the-art research in the E-VRP field. In this paper we outline the following contributions:(i)
a review of the recent energy consumption models that could be used in BEV routing models;(ii)
an updated literature review and a concise table summary of already reviewed E-VRP variants such as GVRP, mixed fleet, BSSs, partial recharges, and different charging technologies;(iii)
a review of the additional emerged E-VRP variants, which include hybrid vehicles, CS siting, nonlinear charging function, dynamic traffic conditions and charging schedule optimization;(iv)
a comprehensive analysis of operation research procedures in the E-VRP, which includes an overview of the procedures employed for solving various E-VRP variants, highlighting state-of-the-art procedures, and a concise table summary of the applied procedures.


### 1.2. Organization of This Paper

The remainder of this paper is organized as follows. In Section2, the E-VRP is described and the literature review of the basic problem formulation is presented. In Section 3, basic characteristics and recent energy consumption models of BEVs are described together with the BEV’s application and evaluation in the delivery processes. In Section 4, variants of the E-VRP are presented with some related problems from the literature. In Section 5, approaches for solving E-VRP are presented, which include state-of-the-art exact, heuristic, metaheuristic, and hybrid procedures. The conclusion and future research directions are given in Section 6.


## 1.1. Recent Literature Reviews and Scientific Contribution

Here we present, to our knowledge, the most recent literature reviews on the E-VRP and related problems.Juan et al. [16] presented a review regarding the environmental, strategic, and operational challenges of EV integration in logistics and transportation activities. The authors performed a comprehensive analysis of environmental challenges including the transportation impact on the pollution and EVs’ possible contribution to the reduction of carbon emissions. Regarding the strategic challenges, the authors presented issues related to CSs: battery swap technology, different charging technologies, CS location problem, limited number of charges, and charging network; issues related to the EVs: mixed fleet of ICEVs and EVs, economic challenges when integrating EVs in a fleet, and routing constraints. The solution approaches for the E-VRP were presented in general, as a class of VRP solving procedures. Similarly, Margaritis et al. [17] presented practical and research challenges of EVs focusing on battery development, lack of charger compatibility, systematic energy management, the lack of optimization procedures that could minimize the EV routing and scheduling decisions, cooling/heating usage, financial sources, novel policies, measures for EV deployment in transport services, etc.Montoya [18] researched several variants of the E-VRP: green VRP (GVRP), E-VRP with partial recharging and nonlinear charging functions, and the technician routing problem with a mixed fleet of ICEVs and EVs. For each problem, effective solving procedures were proposed: multispace sampling heuristic, iterated local search enhanced with heuristic concentration, and two-phase parallel metaheuristic based on solving a set of subproblems and extended set-covering formulation. The authors also formulated a fixed route vehicle charging problem (FRVCP) with and without time windows to optimize the charging decisions for a route with a fixed customer sequence. The authors did not focus as much on reviewing the E-VRP literature, especially not on the procedures for solving the problem. Compared to Montoya [18], this paper did not focus as much on the FRVCP and technician routing problem or on detailed procedures used to solve those problems.The most recent survey on the E-VRP is presented by Pelletier et al. [19] and it includes technical background on EV types and batteries, EV market penetration, EV competitiveness and incentives, and an overview of the existing research regarding EVs in transportation science. The authors provided a comprehensive review on, at the time, the latest solving procedures for the E-VRP with mixed fleet and optimal paths and covered several key papers regarding partial recharges, hybrid vehicles, and different charging technologies.In this paper, a survey on the E-VRP is presented, which includes approaches for solving the E-VRP and related problems that emerged with BEVs integration in the logistic processes. The focus is not on the economic and environmental challenges related to BEVs. Most of the latest literature reviews were published in 2016; hence, to the best of our knowledge, there is no published research that summarizes the state-of-the-art research in the E-VRP field. In this paper we outline the following contributions:(i)
a review of the recent energy consumption models that could be used in BEV routing models;(ii)
an updated literature review and a concise table summary of already reviewed E-VRP variants such as GVRP, mixed fleet, BSSs, partial recharges, and different charging technologies;(iii)
a review of the additional emerged E-VRP variants, which include hybrid vehicles, CS siting, nonlinear charging function, dynamic traffic conditions and charging schedule optimization;(iv)
a comprehensive analysis of operation research procedures in the E-VRP, which includes an overview of the procedures employed for solving various E-VRP variants, highlighting state-of-the-art procedures, and a concise table summary of the applied procedures.


## 1.2. Organization of This Paper

The remainder of this paper is organized as follows. In Section2, the E-VRP is described and the literature review of the basic problem formulation is presented. In Section 3, basic characteristics and recent energy consumption models of BEVs are described together with the BEV’s application and evaluation in the delivery processes. In Section 4, variants of the E-VRP are presented with some related problems from the literature. In Section 5, approaches for solving E-VRP are presented, which include state-of-the-art exact, heuristic, metaheuristic, and hybrid procedures. The conclusion and future research directions are given in Section 6.


## 2. Electric Vehicle Routing Problem

With BEV penetration in logistic distribution processes, a problem of routing a fleet of BEVs has emerged: the E-VRP. The E-VRP aims to design least-cost BEV routes in order to serve a set of customers by taking into account often used constraints: vehicle load capacity, customer time windows, working hours, etc. [3, 20]. Additionally, BEVs have the limited driving range which directly corresponds to more frequent recharging events at CSs. CSs can be built at separate locations as public CSs or mounted at customers’ locations as private CSs. The time needed to travel to a CS and the recharging time are important aspects of fleet routing, especially if customer time windows are taken into account.To the best of our knowledge, the first research regarding the routing of an electric fleet was published by Gonçalves et al. [21], with the authors observing VRPPD using a mixed fleet of EVs and ICEVs. Refueling of the vehicle is performed at the location where the need for the refueling occurred and the total refuel time is computed, based on the total distance traveled by the EV. Conrad and Figliozzi [22] formulated recharging VRP in which vehicles with limited driving range are allowed to refuel at customers’ locations during the route to up to 80% of the vehicle’s battery capacity. The authors described an application of the model and analyzed the impact of different driving ranges, fixed charging times, and time windows on EV routes and solution quality. The results indicated that customer time windows greatly limit route distance when recharging time is long and vehicle range is constrained. Erdoğan and Miller-Hooks [23] formulated GVRP in which a fleet of vehicles is powered by alternative fuels (alternative fuel vehicle, AFV): biodiesel, ethanol, hydrogen, methanol, natural gas, electricity, etc. AFVs can refuel at separately located stations, with fixed refueling time. The authors did not consider customer time windows and vehicle load capacity constraints. New problem-specific instances were developed and two heuristics for solving the problem were applied. The results showed that the limitation of driving range severely increased the number of refueling stations and the total traveled distance in the solution. A similar problem was researched by Omidvar and Tavakkoli-Moghaddam [24], in which the authors added vehicle load constraint, customer time windows, and congestion management, as the vehicle can stay at the customer’s location during the congestion hours. The authors minimized emission and pollution costs by applying commercial software for solving small instances and two metaheuristics for solving larger instances. Schneider et al. [25] published the first research on routing a BEV fleet by taking into account possible visits to CSs and charging time dependent on the SoC level when entering a CS. The problem was formulated as E-VRP with time windows (E-VRPTW) as they took into account load, battery, and time window constraints. The authors formulated E-VRPTW as the mixed integer linear program (MILP) on the complete directed graph G, where customers are modeled as graph vertices and paths between customers are modeled as graph arcs. Here, the basic model formulation is given. Let V={1,…,N} be a set of geographically scattered customers who need to be served, and let F be a set of CSs for BEVs. In order to allow multiple visits to the same CS, a virtual set of CSs F′ is defined. Vertices 0 and N+1 denote the depot, and every route begins with vertex 0 and ends with vertex N+1V0,N+1=V∪{0}∪{N+1}. Graph G is defined as G=(V0,N+1∪F′,A), where A is set of arcs, A={i,j∣i,j∈V0,N+1∪F′,i≠j}. Depending on the real-life constraints, different (i,j) arc values can be interpreted as distance dij, travel time tij, energy consumption eij, speed vij, cost cij, etc. The binary variable xij={0,1} is equal to 1 if arc (i,j) is traversed in the solution, and 0 otherwise. The whole MILP program for the E-VRPTW with equations for load, battery, time windows, flow, and subtour constraints is presented by Schneider et al. [25].In the VRP, it is customary for the primary objective to minimize the total number of vehicles used (1) and then to minimize the total distance traveled (2) or some other objective functions [26]. Total vehicle number is a primary objective as generally greater savings can be achieved with fewer vehicles (vehicle fixed costs, labor cost, etc.). Such an objective is contradictory as with fewer vehicles total traveled distance increases and vice versa. By taking into account the high purchase cost of BEVs, such a hierarchical objective seems justifiable in BEV routing applications [25, 27].(1)min∑j∈V∪F′x0j(2)min∑i∈V0∪F′,j∈VN+1∪F′,i≠jdijxijObjective functions can be complex with simultaneous minimization of vehicle number, total traveled distance [28], total travel times [29], total routing cost and planning horizon [11, 27, 30, 31], GHG emission [32, 33], energy consumption [34–36], etc. Total routing costs of BEVs usually consist of BEV acquisition cost, circulation tax, maintenance, costs related to the energy consumption (electric energy price), cost of battery pack renewal after its lifetime, labor costs, etc. Instead of a single-objective function, some authors use multiobjective function, i.e., fuel consumption and total driving time [37], fuel consumption and route cost [38], battery swapping and charge scheduling [39], etc. An overview of different objectives in E-VRP is presented in Table 1 in the columnObjective.Table 1
Overview of the E-VRP variants and related problems.
Reference Problem name Charging Other Objective C TW MIX LRP F L NL PR DFC BS TD H Bektaş and Laporte [40] PRP X X           Emission model, speed limitation Total costs: labor, fuel and emission as a function of load and speed  Conrad and Figliozzi [22] Recharging VRP X X   X         Vehicle number and total traveling costs: distance, service time, recharging time  Gonçalves et al. [21] VRPPD with mixed fleet X  X  X        No CS location, time constraint Total traveling costs: fixed and variable  Demir et al. [32] PRP X X           Emission model Total traveling costs: labor, fuel and emission as a function of load and speed; speed optimization  Erdoğan and Miller-Hooks [23] GVRP     X        AFVs, limited route duration Vehicle number and total traveled distance  Omidvar and Tavakkoli-Moghaddam [24] GVRP X X   X      X  AFVs, limited fuel capacity and route duration, congestion management Total costs: vehicle fixed costs, distance, time and emission  Abdallah [41] PHEVRPTW X X    X  X    X Electric charge cost is neglected Routing costs: time run on the fossil oil  Barco et al. [34, 42] E-VRP and charge scheduling X X    X  X   X  Energy consumption and battery degradation model, private and public CS, time-dependent energy rates Energy consumption  Davis and Figliozzi [10]  X          X  Energy consumption model with speed profiles, limit route duration and energy consumption (battery capacity), no recharging during the route Total costs: vehicle purchase, energy, maintenance, tax incentive, battery replacement, routing  Van Duin et al. [43] FSMVRPTW, EVFSMVRPTW X X X          No recharge, range constrained by battery, relaxed time window constraints, emission Total costs: vehicle fixed costs, time and distance  Adler and Mirchandani [44] Online routing of BEVs          X   Battery reservations, waiting for fully charged battery The average vehicle delay time  Alesiani and Maslekar [45] BEV routing     X        Waiting time cost at CS is proportional to the number of EVs in CS, limiting the number of CSs in route and number of vehicles in CS, energy consumption model Traveling, charging and energy consumption costs  Demir et al. [37] PRP X X           Fuel consumption and emission model Bi-objective minimization of(1) fuel consumption and (2) total driving time  Felipe et al. [46] GVRP-MTPR X     X  X X     Total recharging costs: fixed and variable  Preis et al. [35] Energy-optimized routing of BEVs X X   X        Energy consumption model Energy consumption  Sassi et al. [47] Sassi et al. [48]Sassi et al. [49] HEVRP-TDMF VRP-HFCC VRP-MFHEV X  X   X  X X    Time-dependent charging costs, operating windows and power limitation of CS, compatibility of BEVs with chargers, electricity grid capacity [47, 49] Vehicle number and total costs: fixed, routing, charging and waiting costs  Schneider et al. [25] E-VRPTW X X    X        Vehicle number and total traveled distance  Zündorf [50] EVRC       X X X X   Battery constrained SPP, different CS types: regular, superchargers and BSS, energy consumption model Travel time  Bruglieri et al. [29, 51] E-VRPTW X X    X  X      Vehicle number and total travel, recharging and waiting time  Goeke and Schneider [30] E-VRPTWMF X X X   X       Energy consumption model: varying BEV load, road slope Vehicle number and different objectives: (i) distance, (ii) costs: labor, driver wage, vehicle propulsion - electric energy and diesel costs, (iii) (ii) + battery replacement cost  Lebeau et al. [52] FSMVRPTW-EV X X X   X       Energy consumption model based on the collected data, recharge only at the depot Total costs: vehicles fixed and operating costs, labor costs  Moghaddam [53] E-VRPTWPR X X    X  X     Capacitated CSs Vehicle number and number of CSs  Pourazarm et al. [54] Single/Multi BEV routing      X  X   X  Homogeneous and in-homogeneous CSs Total time  Schneider et al. [55] VRPIS (EVRPRF) X     X        Total travel and fixed vehicle costs  Yang and Sun [56] BSS-EV-LRP X   X      X    Total routing and construction costs  Desaulniers et al. [57] E-VRPTW -SF/MF/SP/MP X X    X  X      Vehicle number and total routing costs  Doppstadt et al. [58] HEV-TSP      X      X Four modes of travel: combustion, electric, charging and boost mode, no CSs visits Total costs as long as maximal route duration is not overrun  Hiermann et al. [31] E-FSMFTW X X X   X        Vehicle number and total costs: vehicle fixed and routing costs  Keskin and Çatay [28] E-VRPTWPR X X    X  X      Vehicle number and total traveled distance  Koç and Karaoglan [59] GVRP     X        AFVs, limited route duration Vehicle number and total traveled distance  Lin et al. [60] E-VRP X     X       Energy consumption model and load effect Total costs: battery charging, travel time and waiting costs  Masliakova [61] Routing and charging of electric buses X     X   X    Energy consumption model, two types of buses depending on the charging event: en-route or at depot, homogeneous and in-homogeneous CSs Investment and operations costs, a travel time of passengers  Mirmohammadi et al. [62] Periodic green VRP X X X        X  Periodic routing, primary and secondary time windows, static traffic conditions within a period Total emissions, total service time and penalties  Montoya et al. [63] GVRP     X        AFVs, limited route duration Vehicle number and total traveled distance  Roberti and Wen [64] E-TSPTW  X    X  X      Total traveled distance  Schiffer et al. [11, 27] E-LRPTWPR X X  X  X  X     CSs at customers’ locations, pickup and delivery, multiple driver shifts and multiple planning periods, emissions Total costs during the planning period: CS and BEV investment costs, fixed costs (tax, maintenance), distance dependent costs (energy)  Wen et al. [65] E-VSP      X  X     Timetable bus trips, multiple depots, time windows of depots and CSs Total cost: vehicles and traveling costs  Andelmin and Bartolini [66] GVRP     X        AFVs, limited route duration Vehicle number and traveled distance  Çatay and Keskin [67] E-VRPTWPR X X    X  X X    Normal and fast charger at CSs Vehicle number and total recharging costs  Froger et al. [68] E-VRP-NL-C       X X X    Capacitated CSs Total travel, service, charging and waiting time  Hof et al. [69] BSS-EV-LRP X   X      X    Total routing and construction costs  Leggieri and Haouari [70] GVRP     X        AFVs, limited route duration Vehicle number and total traveled distance  Mancini [71] HVRP            X Full instant recharge Total traveled distance with penalties for using internal combustion engine  Montoya et al. [72] E-VRP-NL       X X X    CS types: slow, moderate and rapid Total travel and recharging time  Schiffer and Walther [73] E-LRPTWPR X X  X  X  X      Total distance, number of vehicles and CSs used, total costs: investment costs of BEVs and CSs, and operational costs  Shao et al. [74] EVRP-CTVTT X X   X      X   Total costs: travel, charging, penalty, and fixed vehicle costs  Sweda et al. [75] Adaptive routing and recharging policies of EVs      X  X     Heterogeneous CSs - the probability of being available and expected waiting time, origin-destination pairs Traveling, waiting and recharging costs  Vincent et al. [76] HVRP     X       X Vertex demand in CVRP is associated with travel time for each arc in HVRP Total costs  Amiri et al. [39] BSS location & scheduling      X    X    Multi-objective - minimization of(1) battery charging and power loss costs; (2) deviation from nominal voltage; and (3) network capacity releasing  Bruglieri et al. [77] E-VReP  X    X       One way car sharing service, workers with bicycles go to the EVs locations and relocate them, battery level demand request Multi-objective:(1) minimization of the workers employed; (2) minimization of the duration of the longest route; and (3) maximization of the number of served relocation requests  Joo and Lim [78] EV routing             Energy SPP, no recharging, energy consumption model Minimize energy consumption and average speed on the path  Keskin and Çatay [79] E-VRPTW-FC X X    X  X X     Vehicle number and total recharging costs  Keskin et al. [80] E-VRPTW-FC X X     X  X    M/M/1 queuing system at capacitated CSs, battery capacity restriction, four planning intervals in a day, partial recharge not evident in paper Total cost: energy cost, routing, labor and penalties for late arrivals  Kullman et al. [81] E-VRP-PP       X X X    Public and private CSs, capacitated CS, single charging technology per CS Expected time to visit all the customers  Li et al. [82] MBFM & recharging problem   X  X        Electric, diesel, compressed natural gas and hybrid-diesel buses Total network benefit of replacing old vehicles with new ones within the planning horizon and budget constraints  Lu et al. [83] MTFSP X  X  X      X  Travel request are known a priori in time-varying origin-destination tables, service and deadheaded trips with different consumption rate, BEVs have higher priority than ICEVs Total operating cost of a taxi company  Masmoudi et al. [84] DARP-EV X X X       X   Energy consumption model of Genikomsakis and Mitrentsis [85] with constant speed, acceleration and road slope; different vehicle resources: accompanying person seat, handicapped person seat, stretcher and/or a wheelchair, limited user ride time Total routing costs (distance)  Paz et al. [86] MDEVLRPTW - BS/PR/BSPR X X  X  X  X  X   Three MIP models depending on the partial recharge and BSS Total traveled distance  Pelletier et al. [87] EFV-CSP       X X X    Preemptive charging with a limited number of chargers and charging events at the depot, time-dependent energy costs, FRD charge, grid restriction, cyclic and calendar battery degradation Total charging costs  Poonthalir and Nadarajan [38] F-GVRP     X        Fuel consumption model, varying speed, AFVs, limited route duration Bi-objective:(1) routing costs; and (2) fuel consumption  Schiffer and Walther [88] LRPIF X X  X  X  X     Loading and refueling facilities Total costs: investment costs of vehicles and facilities, routing costs  Schiffer and Walther [89] RELRPTWPR X X  X  X  X     Uncertain customer pattern scenarios over working days regarding the spatial customer distribution, demand and service time windows Total costs: investment costs of vehicles and facilities, routing costs  Shao et al. [90] E-VRP X    X        Energy consumption model: cargo load, uncertain travel speed Total costs: travel, charging and fixed vehicle costs  Wang et al. [91] BEV routing      X       Parking fee, capacitated CSs - queuing time Multi-objective minimization of(1) travel time: driving, queuing and charging time; (2) charging costs: electricity, service and parking fee; (3) energy consumption  Zhang et al. [36] E-VRP X            Energy consumption model, emissions, static speed, charging time unknown Energy consumption  Basso et al. [92] 2sEVRP X X    X   X    Each CSs can have different charging rate, energy consumption model for road segments: speed profile, road slope, acceleration Total energy consumption  Breunig et al. [93] E2EVRP X            Two echelons - first ICEVs and second BEVs, full recharge when visiting CSs, charging time unknown Total routing costs  Bruglieri et al. [94] GVRP     X        AFVs, limited route duration Vehicle number and total traveled distance  Froger et al. [95] E-VRP-NL       X X X     Total travel and charging time  Hiermann et al. [96] H 2E-FTW X X X   X  X    X Propulsion mode decision Total costs: fixed and variable  Jie et al. [97] 2E-EVRP-BSS X         X   Routing in two echelons, sensitivity analysis of battery driving range and vehicle emissions Total routing costs, the battery swapping costs and the handling costs at the satellites  Koyuncu and Yavuz [98] MGVRP  X X  X X  X X    ICEVs (fixed recharge time) and AFVs, node- and -arc MILP formulation, single recharge technology per CS, additional modeling: customer demands, customer vehicle restrictions, subscription or pay-as-you-go refueling costs, completely heterogeneous fleet, last-mile delivery and closed time windows Total traveling cost  Macrina et al. [99] GMFVRP-PRTW X X X   X  X X  X  Single recharge technology per CS, but different charging technologies between CSs, energy consumption model for road segments with time-dependent speeds Cost of energy recharged during the route and at the depot, fuel costs, and cost related to traveled distance  Macrina et al. [33] GMFVRP-PRTW X X X   X  X X    Single recharge technology per CS, but different charging technologies between CSs Recharging, routing and activation costs, limit emissions  Normasari et al. [100] CGVRP X    X        AFVs, limited route duration Vehicle number and total traveled distance Reference: referenced paper; Problem name: the name of the analyzed problem; 13 columns representing characteristics of the problem in the following order: C: vehicle load (cargo) capacity, TW: customer time windows, MIX: heterogeneous (mixed) fleet, LRP: location routing problem, F: fixed (constant) refuel (recharge) time, L: linear charging process, NL: nonlinear charging process, PR: partial recharge strategy, DFC: different charging technologies, BS: battery swap strategy, TD: time-dependent travel times, H: hybrid vehicles, Other: some special characteristic of the problem; Objective: the objective function for the optimization.


## 3. Battery Electric Vehicles in Delivery Processes

The major problem that BEVs in delivery processes are facing is the limited driving range. Grunditz and Thiringer [101] analyzed over 40 globally available BEVs, which can be categorized into small, medium-large, high-performing, and sports cars. All of the BEV models utilize lithium-based batteries, especially lithium-ion [102] with battery capacity and distance varying within the ranges 12-90 kWh and 85-528 km, respectively. An average medium-sized personal BEV has a battery capacity of 30 kWh, which is enough to travel 250 km. In the delivery process, mostly light vans and freight BEVs are used, which have a shorter driving range (160-240 km) compared to the driving range of ICEVs (480-650 km) [13, 43, 103]. The reason is that the battery has lower specific energy (130 Wh/kg) than the fossil oil (1233 Wh/kg), and the amount of energy that can be stored in the battery is much lower than in the fossil oil. Batteries mounted in BEVs are mostly the main cause of high acquisition costs and technical limitations as the battery degrades over time, resulting in decreased maximal capacity. Pelletier et al. [104] concluded that the battery should be replaced after five to ten years or after 1,000 to 2,000 cycles with large SoC variations. The authors also described factors that influence such battery degradation: overcharging, overdischarging, high and low temperatures, high SoC during storage, large depth of discharge, etc., and they presented battery degradation models that can be used in goods distribution with BEVs.


### 3.1. BEV Application

BEVs are more likely to be used on short distances and/or in urban areas where they are more effective than ICEVs due to the low driving speed, low noise production, frequent stops, and financial incentives. In cases when the average route length is short, such as the average FedEx route length in the USA, which is 68 km [12], BEVs can be applied directly and recharging can be performed on return to the depot. BEVs are already being applied in such occasions: DHL, UPS, FedEx, and Coca-Cola [105, 106] use BEVs mostly for last-mile deliveries as distances are shorter and vehicle loads are lower. Many companies are performing case studies of integrating BEVs in their delivery fleet. Lin et al. [60] were debating the use of BEVs in time-precise deliveries as the long recharging time at CS causes hard completion of an on-time delivery, which then significantly increases the overall routing costs. Davis and Figliozzi [10] and van Duin et al. [43] evaluated a wide range of scenarios to compare the routing costs of ICEVs and BEVs. Van Duin et al. [43] concluded that BEVs have the ability to efficiently perform urban freight transport and meanwhile to reduce the GHG emissions and noise nuisance. Davis and Figliozzi [10] reported that BEVs are not competitive if the solution to the same problem results in a higher number of BEVs than the number of ICEVs. For BEVs to be competitive, the authors pointed out a combination of several key elements: high daily distance (as much as maximum BEV driving range), low speeds and congestions, frequent customer stops, the reduction of a BEV’s purchase cost by tax incentives or technology development, and long planning horizon. On the other hand, Schiffer et al. [27] conducted a case study using freight BEVs for deliveries with a planning horizon of five years and compared the results to the delivery done by conventional freight trucks. Several characteristics of the performed case study are important: (i) delivery radius up to 190 km from the depot; (ii) CSs located at customers’ locations, which allows simultaneous charging while unloading goods; (iii) already existing strong current at CSs; and (iv) vehicles returning to the depot at least once a day. Overall, the authors concluded that there is no operational limitation when using BEVs compared to ICEVs, as the used number of vehicles and total traveled distance are competitive, and at the same time, the overall costs are lower with almost 25% less CO2 emission. A year later, Schiffer et al. [11] repeated the case study with more realistic costs when a strong current at the CS is not available and concluded that, with the higher CS investment costs, BEVs are no longer competitive. To fully assess the integration of BEVs in logistic processes, the authors pointed out three key elements: different network structures, future CO2 emission policies, and future technology development (battery capacity and charging infrastructure).


### 3.2. Energy Consumption

Due to the low specific energy, energy consumption should be precisely estimated in order to achieve a BEV’s maximal driving range and to reduce the overall routing costs. The energy consumption can be estimated by simulation models but due to the complexity of the E-VRP and unknown driving cycles in advance, mostly macroscopic models with several real-world approximations are applied in the BEV routing models. In the available literature, energy consumption is often estimated using longitudinal dynamics model (LDM). Here, the LDM of Asamer et al. [107] is presented. Force F needed to accelerate and to overcome resistances (grade, rolling, and air) is given by (3), where m is vehicle mass (mostly empty vehicle), a acceleration, v vehicle speed, g gravitational constant, f the inertia force of vehicle rotating parts (up to 5% of the total vehicle mass), α road slope, cr rolling friction coefficient, cd air drag coefficient, ρ air density, and A vehicle frontal air surface. If F≥0, the vehicle is accelerating and power is needed for the movement of BEV (motor mode); otherwise, if F<0, deceleration (braking) or driving downhill is occurring and energy is returned into the BEV’s battery as the electric engine has the ability to return the energy (recuperating mode). By process of recuperation, up to 15% of totally consumed energy can be returned [51, 108]. Electric power that comes from the battery is divided into the auxiliary power P0 and mechanical power Pm=Fv. Auxiliary power is spent on the electronic devices in the vehicle: heating, ventilation, light, etc., which can shorten the BEV’s range up to 30% [109]. Battery power Pb can be computed by (4), where μm is the transmission coefficient between the electric motor and drivetrain, μe is the conversion ratio from chemical energy in the battery to electric energy, and μg is the conversion ratio from mechanical energy on wheels to chemical energy stored in the battery. Energy is returned into the battery only if the force F is lower than zero and speed is higher than the experimentally determined value vmin [107]. Energy consumption can be computed by the time integration of (4).(3)F=mgsin⁡α︸Grade+crmgcos⁡α︸Rolling+0.5cdρAv2︸Air+fma︸Acc.(4)Pb=μeμmFv+P0,ifF≥00,ifv≤vminμgFv+P0,else,ifF<0Goeke and Schneider [30] extended the energy consumption model by taking into account variable vehicle load mass when delivering goods without acceleration and braking processes. In the CVRP variants, as customers are being served, vehicle load is decreasing. The authors concluded that actual load strongly improves the solution quality, as a large number of solutions generated without taking into account load distribution tend to be infeasible due to the violation of battery capacity or time windows.Genikomsakis and Mitrentsis [85] presented a more realistic electric engine model: the load-efficiency curve is approximated by piecewise function and the normalization factor is added to take into account the motor size. The authors used the recuperation energy factor dependent on the vehicle speed as follows: below vmin there is no energy recuperation, beyond vmax maximum energy is recuperated, and in between linear interpolation of energy recuperation is assumed. Therefore, regarding the energy recuperation the authors observed two cases: (i) when recuperated energy exceeds the consumption of the auxiliary devices and the excess of energy is stored into the battery; (ii) when recuperated energy is not sufficient to cover the consumption of the auxiliary devices and thus the energy is drawn from the battery. The authors compared their model on nine characteristic driving cycles (short/long, congested/uncongested, highway, etc.) to the FASTSim simulation tool [110], and this resulted in a relative energy consumption error of up to 4%. The developed energy consumption model was used to create a database of coefficients for several key characteristics: motor type, motor power, battery type, road type, road slope, road speed limit, etc.Macrina et al. [99] developed an energy model for mixed GVRP with partial recharging and time windows. The authors modeled vehicle speed on arc (i,j) with three phases h: acceleration (h=1), constant speed (h=2), and deceleration (h=3), resulting with either triangular function (1→3) or trapezoidal function (1→2→3). The fuel consumption of ICEV is based on the similar LDM model, presented by (5) and (6), where ui is the load (cargo) weight, ζ fuel-to-air mass ratio, κ heating value of typical diesel fuel, ψ conversion factor, k engine friction coefficient, Ne radial engine speed, De engine displacement, μd diesel engine efficiency, μdt drivetrain efficiency, and t¯h travel spent in phase h [40, 111]. For the consumption of BEV, the authors proposed a similar model to that presented by (7), where ηh+ and ηh- are the efficiency of the electric engine in motor and recuperating mode. The authors compared the proposed energy consumption model to the basic model, where energy consumed is proportional to the distance traveled, and the energy consumption model of Goeke and Schneider [30], where acceleration and braking processes are not taken into account. The results showed that, compared to the proposed energy consumption model, basic energy consumption model and the model of Goeke and Schneider [30] produce an average relative error of 70% and 4%, respectively.(5)Fijhui=gsin⁡α+crgcos⁡α+atijhm+ui+0.5cdρAvtijh2(6)FijMui=∑h=1,2,3ζκψkNeDe+Fijhuivtijh︸PijhMui/μdμdtt¯h(7)EijEui=∑h=1,2,3Fijhuivtijh/ηht¯h︸EijhEui,ηh=ηh+≤1,EijhEui>0&0≤PijhMui≤100kWh=1,2ηh-≥1,EijhEui<0&-100≤PijhMui≤0kWh=3Basso et al. [92] developed an energy consumption model for two-stage E-VRP (2sEVRP) that includes detailed topography and speed profiles. The similar LDM model of Asamer et al. [107] is applied for computing the energy consumption on a road segment (link), but with time-dependent speed, varying mass, constant slope, and link distance. For each link, similar to Macrina et al. [99], the three or two characteristic phases of speed curve can be distinguished. In such a way, acceleration and breaking processes before and after the intersection are taken into account. As mass changes during the delivery, the energy consumption is expressed as a linear function of mass. Then, the shortest energy paths between all vertices pairs (customers, depot, CSs) without charging constraint are determined by applying the Bellman-Ford algorithm [112]. A nominal mass value between the mass of full and empty vehicle is used in the computation. The experiments showed an average energy estimation error of 2.28% and improved energy feasibility compared to some of the previous consumption models. The authors also reported that the use of average speed consumed at least 24% less energy.Asamer et al. [107] analyzed the data collected from the BEV trips in order to determine the dependency between BEV energy share and average trip speed. Results showed that on lower speeds v≤30 km/h, most of the energy is spent on the acceleration and auxiliary devices, while on the higher speeds v≥80 km/h, approx. 70% of total energy is spent on overcoming the air drag force. The slope of the terrain in total energy consumption has the lowest share, and its value is even decreasing with the increase of average trip speed. The rolling resistance energy consumption share does not depend on the average trip speed.Preis et al. [35] analyzed the energy optimal routing of BEV, where the routes were designed in different terrain slopes and battery capacity scenarios. The authors concluded that energy savings grow linearly with the maximum altitude difference and that decreasing the battery capacity of BEVs does not affect the recharging schedule but increases the total energy consumption. Fiori et al. [113] compared the power-based consumption model of BEV and ICEV [114]. The authors concluded that BEVs and ICEVs have different fuel/energy-optimized assignment. The faster routes increase the BEV’s energy consumption while the congested and low-speed arterial routes consume less energy. Masmoudi et al. [84] compared the realistic energy consumption model of Genikomsakis and Mitrentsis [85] to the constant consumption model (237.5 W/km) on instances for the dial-a-ride problem with EVs. The authors concluded that using the realistic model is more efficient with 0.14% difference between realistic and constant energy consumption from the best-known solutions (BKS). To emphasize the importance of the energy consumption model and energy minimization, Zhang et al. [36] compared the distance and energy minimization and concluded that the distance-minimizing objective consumes 16.44% more energy than the energy-minimizing objective.In real-life conditions, speeds on the roads are time-dependent and can be described as the speed profile over the observed time period. Speed profile depends on the road type, driver behavior, traffic (accidents, recurrent congestions), weather conditions, etc. [115, 116]. A large number of parameters make it hard to predict a BEV’s energy consumption in different traffic scenarios. Therefore, some researchers apply data-driven approaches to predict the energy consumption of a BEV. De Cauwer et al. [117] developed a model for BEV energy consumption by applying multiple linear regression to real-world measured BEV data. The proposed multiple linear regression model (MLR) has seven features: distance, speed, energy consumed by auxiliary devices, positive elevation, negative elevation, temperature, and kinetic energy change per unit distance. Results showed that the prediction error of trip energy consumption is within 25%. De Cauwer et al. [115] applied a similar MLR model for predicting energy consumption on road segments, suited for BEV routing. The authors used a neural network based on the road-traffic and weather-related features to predict the speed profile of a road segment. The error prediction of the trip’s energy consumption is 12-14%. Lebeau et al. [52] used real-life data to model the energy consumed and recuperated on the trip through least square analysis. The used function depends on the trip duration, measured energy consumption, temperature, and correction parameter. The R2 for the energy consumption model is 0.93 and 0.77 for the recuperated energy model. Wu et al. [118] presented an empirical and analytical method that can estimate a BEV’s instantaneous power in real time and overall trip energy consumption, with the average prediction error up to 15.6%. Fiori and Marzano [119] proposed a backward microscopic power-based LDM energy consumption model based on the known driving cycle. The acceleration, speed, and regenerative braking efficiency were modeled as functions of time. Based on the real driving BEV data, for each vehicle, six parameters were optimized: three parameters regarding the drivetrain and battery efficiency, the parameter for determining when the regenerative braking is occurring, and two traction power thresholds. The authors concluded that the proposed model is flexible enough to be applied to any kind of driving cycle and BEV type with given kinematic profile and characteristics.Pelletier et al. [104] presented battery degradation models that could be used in BEV routing applications, as batteries are a crucial part of the economic BEV routing. Several lithium-ion battery degradation mechanisms with storage and operating conditions that affect battery lifespan were described. Barco et al. [42] applied a battery degradation model based on the temperature, SoC, and depth of discharge in their airport shuttle service with BEVs. It was observed that the consideration of the battery degradation model affects the charging patterns.


## 3.1. BEV Application

BEVs are more likely to be used on short distances and/or in urban areas where they are more effective than ICEVs due to the low driving speed, low noise production, frequent stops, and financial incentives. In cases when the average route length is short, such as the average FedEx route length in the USA, which is 68 km [12], BEVs can be applied directly and recharging can be performed on return to the depot. BEVs are already being applied in such occasions: DHL, UPS, FedEx, and Coca-Cola [105, 106] use BEVs mostly for last-mile deliveries as distances are shorter and vehicle loads are lower. Many companies are performing case studies of integrating BEVs in their delivery fleet. Lin et al. [60] were debating the use of BEVs in time-precise deliveries as the long recharging time at CS causes hard completion of an on-time delivery, which then significantly increases the overall routing costs. Davis and Figliozzi [10] and van Duin et al. [43] evaluated a wide range of scenarios to compare the routing costs of ICEVs and BEVs. Van Duin et al. [43] concluded that BEVs have the ability to efficiently perform urban freight transport and meanwhile to reduce the GHG emissions and noise nuisance. Davis and Figliozzi [10] reported that BEVs are not competitive if the solution to the same problem results in a higher number of BEVs than the number of ICEVs. For BEVs to be competitive, the authors pointed out a combination of several key elements: high daily distance (as much as maximum BEV driving range), low speeds and congestions, frequent customer stops, the reduction of a BEV’s purchase cost by tax incentives or technology development, and long planning horizon. On the other hand, Schiffer et al. [27] conducted a case study using freight BEVs for deliveries with a planning horizon of five years and compared the results to the delivery done by conventional freight trucks. Several characteristics of the performed case study are important: (i) delivery radius up to 190 km from the depot; (ii) CSs located at customers’ locations, which allows simultaneous charging while unloading goods; (iii) already existing strong current at CSs; and (iv) vehicles returning to the depot at least once a day. Overall, the authors concluded that there is no operational limitation when using BEVs compared to ICEVs, as the used number of vehicles and total traveled distance are competitive, and at the same time, the overall costs are lower with almost 25% less CO2 emission. A year later, Schiffer et al. [11] repeated the case study with more realistic costs when a strong current at the CS is not available and concluded that, with the higher CS investment costs, BEVs are no longer competitive. To fully assess the integration of BEVs in logistic processes, the authors pointed out three key elements: different network structures, future CO2 emission policies, and future technology development (battery capacity and charging infrastructure).


## 3.2. Energy Consumption

Due to the low specific energy, energy consumption should be precisely estimated in order to achieve a BEV’s maximal driving range and to reduce the overall routing costs. The energy consumption can be estimated by simulation models but due to the complexity of the E-VRP and unknown driving cycles in advance, mostly macroscopic models with several real-world approximations are applied in the BEV routing models. In the available literature, energy consumption is often estimated using longitudinal dynamics model (LDM). Here, the LDM of Asamer et al. [107] is presented. Force F needed to accelerate and to overcome resistances (grade, rolling, and air) is given by (3), where m is vehicle mass (mostly empty vehicle), a acceleration, v vehicle speed, g gravitational constant, f the inertia force of vehicle rotating parts (up to 5% of the total vehicle mass), α road slope, cr rolling friction coefficient, cd air drag coefficient, ρ air density, and A vehicle frontal air surface. If F≥0, the vehicle is accelerating and power is needed for the movement of BEV (motor mode); otherwise, if F<0, deceleration (braking) or driving downhill is occurring and energy is returned into the BEV’s battery as the electric engine has the ability to return the energy (recuperating mode). By process of recuperation, up to 15% of totally consumed energy can be returned [51, 108]. Electric power that comes from the battery is divided into the auxiliary power P0 and mechanical power Pm=Fv. Auxiliary power is spent on the electronic devices in the vehicle: heating, ventilation, light, etc., which can shorten the BEV’s range up to 30% [109]. Battery power Pb can be computed by (4), where μm is the transmission coefficient between the electric motor and drivetrain, μe is the conversion ratio from chemical energy in the battery to electric energy, and μg is the conversion ratio from mechanical energy on wheels to chemical energy stored in the battery. Energy is returned into the battery only if the force F is lower than zero and speed is higher than the experimentally determined value vmin [107]. Energy consumption can be computed by the time integration of (4).(3)F=mgsin⁡α︸Grade+crmgcos⁡α︸Rolling+0.5cdρAv2︸Air+fma︸Acc.(4)Pb=μeμmFv+P0,ifF≥00,ifv≤vminμgFv+P0,else,ifF<0Goeke and Schneider [30] extended the energy consumption model by taking into account variable vehicle load mass when delivering goods without acceleration and braking processes. In the CVRP variants, as customers are being served, vehicle load is decreasing. The authors concluded that actual load strongly improves the solution quality, as a large number of solutions generated without taking into account load distribution tend to be infeasible due to the violation of battery capacity or time windows.Genikomsakis and Mitrentsis [85] presented a more realistic electric engine model: the load-efficiency curve is approximated by piecewise function and the normalization factor is added to take into account the motor size. The authors used the recuperation energy factor dependent on the vehicle speed as follows: below vmin there is no energy recuperation, beyond vmax maximum energy is recuperated, and in between linear interpolation of energy recuperation is assumed. Therefore, regarding the energy recuperation the authors observed two cases: (i) when recuperated energy exceeds the consumption of the auxiliary devices and the excess of energy is stored into the battery; (ii) when recuperated energy is not sufficient to cover the consumption of the auxiliary devices and thus the energy is drawn from the battery. The authors compared their model on nine characteristic driving cycles (short/long, congested/uncongested, highway, etc.) to the FASTSim simulation tool [110], and this resulted in a relative energy consumption error of up to 4%. The developed energy consumption model was used to create a database of coefficients for several key characteristics: motor type, motor power, battery type, road type, road slope, road speed limit, etc.Macrina et al. [99] developed an energy model for mixed GVRP with partial recharging and time windows. The authors modeled vehicle speed on arc (i,j) with three phases h: acceleration (h=1), constant speed (h=2), and deceleration (h=3), resulting with either triangular function (1→3) or trapezoidal function (1→2→3). The fuel consumption of ICEV is based on the similar LDM model, presented by (5) and (6), where ui is the load (cargo) weight, ζ fuel-to-air mass ratio, κ heating value of typical diesel fuel, ψ conversion factor, k engine friction coefficient, Ne radial engine speed, De engine displacement, μd diesel engine efficiency, μdt drivetrain efficiency, and t¯h travel spent in phase h [40, 111]. For the consumption of BEV, the authors proposed a similar model to that presented by (7), where ηh+ and ηh- are the efficiency of the electric engine in motor and recuperating mode. The authors compared the proposed energy consumption model to the basic model, where energy consumed is proportional to the distance traveled, and the energy consumption model of Goeke and Schneider [30], where acceleration and braking processes are not taken into account. The results showed that, compared to the proposed energy consumption model, basic energy consumption model and the model of Goeke and Schneider [30] produce an average relative error of 70% and 4%, respectively.(5)Fijhui=gsin⁡α+crgcos⁡α+atijhm+ui+0.5cdρAvtijh2(6)FijMui=∑h=1,2,3ζκψkNeDe+Fijhuivtijh︸PijhMui/μdμdtt¯h(7)EijEui=∑h=1,2,3Fijhuivtijh/ηht¯h︸EijhEui,ηh=ηh+≤1,EijhEui>0&0≤PijhMui≤100kWh=1,2ηh-≥1,EijhEui<0&-100≤PijhMui≤0kWh=3Basso et al. [92] developed an energy consumption model for two-stage E-VRP (2sEVRP) that includes detailed topography and speed profiles. The similar LDM model of Asamer et al. [107] is applied for computing the energy consumption on a road segment (link), but with time-dependent speed, varying mass, constant slope, and link distance. For each link, similar to Macrina et al. [99], the three or two characteristic phases of speed curve can be distinguished. In such a way, acceleration and breaking processes before and after the intersection are taken into account. As mass changes during the delivery, the energy consumption is expressed as a linear function of mass. Then, the shortest energy paths between all vertices pairs (customers, depot, CSs) without charging constraint are determined by applying the Bellman-Ford algorithm [112]. A nominal mass value between the mass of full and empty vehicle is used in the computation. The experiments showed an average energy estimation error of 2.28% and improved energy feasibility compared to some of the previous consumption models. The authors also reported that the use of average speed consumed at least 24% less energy.Asamer et al. [107] analyzed the data collected from the BEV trips in order to determine the dependency between BEV energy share and average trip speed. Results showed that on lower speeds v≤30 km/h, most of the energy is spent on the acceleration and auxiliary devices, while on the higher speeds v≥80 km/h, approx. 70% of total energy is spent on overcoming the air drag force. The slope of the terrain in total energy consumption has the lowest share, and its value is even decreasing with the increase of average trip speed. The rolling resistance energy consumption share does not depend on the average trip speed.Preis et al. [35] analyzed the energy optimal routing of BEV, where the routes were designed in different terrain slopes and battery capacity scenarios. The authors concluded that energy savings grow linearly with the maximum altitude difference and that decreasing the battery capacity of BEVs does not affect the recharging schedule but increases the total energy consumption. Fiori et al. [113] compared the power-based consumption model of BEV and ICEV [114]. The authors concluded that BEVs and ICEVs have different fuel/energy-optimized assignment. The faster routes increase the BEV’s energy consumption while the congested and low-speed arterial routes consume less energy. Masmoudi et al. [84] compared the realistic energy consumption model of Genikomsakis and Mitrentsis [85] to the constant consumption model (237.5 W/km) on instances for the dial-a-ride problem with EVs. The authors concluded that using the realistic model is more efficient with 0.14% difference between realistic and constant energy consumption from the best-known solutions (BKS). To emphasize the importance of the energy consumption model and energy minimization, Zhang et al. [36] compared the distance and energy minimization and concluded that the distance-minimizing objective consumes 16.44% more energy than the energy-minimizing objective.In real-life conditions, speeds on the roads are time-dependent and can be described as the speed profile over the observed time period. Speed profile depends on the road type, driver behavior, traffic (accidents, recurrent congestions), weather conditions, etc. [115, 116]. A large number of parameters make it hard to predict a BEV’s energy consumption in different traffic scenarios. Therefore, some researchers apply data-driven approaches to predict the energy consumption of a BEV. De Cauwer et al. [117] developed a model for BEV energy consumption by applying multiple linear regression to real-world measured BEV data. The proposed multiple linear regression model (MLR) has seven features: distance, speed, energy consumed by auxiliary devices, positive elevation, negative elevation, temperature, and kinetic energy change per unit distance. Results showed that the prediction error of trip energy consumption is within 25%. De Cauwer et al. [115] applied a similar MLR model for predicting energy consumption on road segments, suited for BEV routing. The authors used a neural network based on the road-traffic and weather-related features to predict the speed profile of a road segment. The error prediction of the trip’s energy consumption is 12-14%. Lebeau et al. [52] used real-life data to model the energy consumed and recuperated on the trip through least square analysis. The used function depends on the trip duration, measured energy consumption, temperature, and correction parameter. The R2 for the energy consumption model is 0.93 and 0.77 for the recuperated energy model. Wu et al. [118] presented an empirical and analytical method that can estimate a BEV’s instantaneous power in real time and overall trip energy consumption, with the average prediction error up to 15.6%. Fiori and Marzano [119] proposed a backward microscopic power-based LDM energy consumption model based on the known driving cycle. The acceleration, speed, and regenerative braking efficiency were modeled as functions of time. Based on the real driving BEV data, for each vehicle, six parameters were optimized: three parameters regarding the drivetrain and battery efficiency, the parameter for determining when the regenerative braking is occurring, and two traction power thresholds. The authors concluded that the proposed model is flexible enough to be applied to any kind of driving cycle and BEV type with given kinematic profile and characteristics.Pelletier et al. [104] presented battery degradation models that could be used in BEV routing applications, as batteries are a crucial part of the economic BEV routing. Several lithium-ion battery degradation mechanisms with storage and operating conditions that affect battery lifespan were described. Barco et al. [42] applied a battery degradation model based on the temperature, SoC, and depth of discharge in their airport shuttle service with BEVs. It was observed that the consideration of the battery degradation model affects the charging patterns.


## 4. Variants of the Electric Vehicle Routing Problem

Many different VRP variants were researched over the years. With BEV appearance, researchers started to adapt them to the E-VRP context. Due to the specific characteristics of BEV routing, some new problem-specific variants emerged. Here, we present some of the most relevant E-VRP variants and related problems.


### 4.1. Energy Shortest Path Problem and the Electric Traveling Salesman Problem

The energy shortest path problem (E-SPP) and electric TSP (E-TSP) can be considered as two of the simplest forms of the E-VRP. In most of the VRP variants, graph arcs have positive weight values that can represent distance, travel time, cost, etc. To compute the shortest path between customers, the most commonly applied algorithms are Dijkstra, Bellman-Ford,A∗, contraction hierarchies, etc. [108, 120]. By taking into account the recuperated energy of BEV, some arc weights could have a negative value, which makes most of the shortest path algorithms inapplicable. To overcome this problem, Artmeier et al. [108] formalized energy-efficient routing with rechargeable batteries as a special case of the constrained SPP (CSPP) in which the battery charge is limited with its capacity, and there is no recharging at CS. The authors adapted the Bellman-Ford shortest path algorithm with time complexity O(n3) in order to solve the energy CSPP. To overcome negative graph edges, Eisner et al. [121] used Johnson’s shifting technique [122] to transform negative edge cost functions into nonnegative ones and applied Dijsktra’s algorithm with time complexity of O(nlog(n)+m) [123]. Storandt [124] introduced CSPP for EVs with battery swapping, while Sweda and Klabjan [125] added recharging events at CSs with a differentiable charging function. Zündorf [50] solved the CSPP for BEV routing with battery constraints, different types of CSs, and nonlinear charging process. The author developed a charging function propagating algorithm in order to minimize travel time and applied CH and A∗ algorithms to solve the problem. Liao et al. [103] considered the BEV’s shortest travel path problem and presented a dynamic programming algorithm for solving the problem that runs in O(kn2), where k is the upper bound on the number of BSSs. More recently, Strehler et al. [126] observed recharging events in the graph vertices and arcs for the energy-efficient shortest routes of BEVs and HEVs, while Zhang et al. [127] dealt with the electric vehicle route planning with recharging problem (EVRC), which minimizes the overall travel and charging time and takes into account CSs’ locations, partial recharging, nonlinear charging functions, service time duration, and service frequency at CS.As VRP is the generalization of the TSP, the E-VRP is closely related to the E-TSP, in which a set of customers has to be served by only one BEV. Roberti and Wen [64] formulated the E-TSP with time windows (E-TSPTW) as a compact MILP program with binary variables for recharging paths with intermediate stops. The authors observed full and partial recharge policies and applied three-phase heuristic for solving the problem. The authors solved multiple instances, among them the 13 E-VRPTW instances of Schneider et al. [25], which were solved optimally with only one vehicle. Doppstadt et al. [58] formulated the HEV-TSP with four working modes of HEV and provided new test instances. Liao et al. [103] introduced the EV touring problem, as the generalization of TSP, with the minimization of the total required time. Two scenarios with battery swap scheme were observed: the on-site station model, in which each city has a BSS; the off-site station model, in which BSS is located at an acceptable distance from the city. The authors proposed efficient polynomial time algorithms for the problem.


### 4.2. Heterogeneous or Mixed Vehicle Fleet

In today’s vehicle fleets, mostly ICEVs are present. Transition to an almost wholly electric fleet is a very challenging economic task. Therefore, most companies are gradually integrating BEVs into their existing ICEV fleet. Routing algorithms have to be upgraded and adapted for electric fleet characteristics, as on-time priority is harder to achieve.Fleet size and mix VRP (FSM-VRP) was first introduced by Golden et al. [128], where the authors considered routing a fleet of vehicles with different acquisition costs and the routing cost is dependent on the vehicle type. Goeke and Schneider [30] formulated E-VRPTW and mixed fleet (E-VRPTWMF) with equal vehicle load capacities of BEVs and ICEVs, and equal BEV battery capacities. The authors compared the percentage share of BEVs to the total traveled distance obtained with three objective functions: total traveled distance, total costs with battery costs, and total costs without battery costs. The traveled distance on one battery pack was assumed to be 241,350 km (150,000 miles) with the payment of an additional $600 per kWh in case of replacement. The results showed that the BEV’s share in total traveled distance increased significantly when battery costs are not considered, while in the other two scenarios ICEVs performed most of the deliveries. Hiermann et al. [31] analyzed a similar problem, the electric fleet size and mix VRPTW and recharging stations (E-FSMFTW), with different vehicle load and battery capacities. The authors showed the overall positive influence of the heterogeneous vehicle fleet on the generalized total cost function. In most of the instances, three to four vehicle types are used in the final solution. A similar analysis on small-sized instances with different BEVs, ICEVs, and HEVs was undertaken by Lebeau et al. [52]. The authors defined seven groups of vehicle types that could be used for the delivery, from small vans and quadricycles through diesel-only or electric-only groups to a group of all vehicle types. The results showed the following aspects of routing: (i) the fleet with different vehicle types reduced the total routing costs; (ii) in the large van group, ICEVs outperformed BEVs; and (iii) HEVs showed a great application in deliveries solely made by trucks. Sassi et al. [47–49] also observed a heterogeneous fleet of ICEVs and BEVs with different load and battery capacities, different operating costs, time-dependent charges, and the compatibility of BEVs and chargers at CSs. The authors focused on the procedures for solving the problem, so they did not offer any comparison between the ICEV and BEV solutions. Hiermann et al. [96] introduced the Hybrid heterogeneous electric fleet routing problem with time windows and recharging stations (H2E-FTW), in which the fleet consists of ICEVs, BEVs, and PHEVs. The authors compared the overall costs of the solutions obtained with a homogeneous fleet of ICEVs, BEVs, and PHEVs to the optimized solution with a mixed fleet. The gap for the ICEV fleet is the largest, with an average value of 60% in the extreme case when the fuel cost is the highest. In contrast, when the electric cost is the highest, the BEV fleet on average produces a 40% gap, while the PHEV fleet shows the lowest gap value, up to 25%. The authors pointed out that, in mixed solutions, BEVs are preferred for clustered instances, ICEVs for randomly distributed instances, and PHEVs for randomly clustered and distributed instances. Overall, the results showed that operational costs can be 7% lower in a mixed case compared to the homogeneous case.


### 4.3. Hybrid Vehicles

As compensation for the limited driving range of BEVs, HEVs, which have both an internal combustion engine and an electric engine, have been developed. Two main types of HEVs are present on the market: theseries hybrid, in which only the electric motor drives the train and the internal combustion engine is used to recharge the battery pack, and theparallel hybrid, which uses both internal combustion engine and electric engine to drive the train, where the electric engine is more efficient in stop-and-go activities and the internal combustion engine is more efficient at high speeds. We focus on the PHEV as a version of aparallel hybrid in which batteries can be recharged by connecting a plug to the electric power source. The PHEVs have an option to decide during the route to run on either electric energy or fossil oil. This enables the visiting of customers far from the depot with almost no refuel/recharge during the route. As PHEVs have two engines, their load is heavier compared to the BEVs and ICEVs, and therefore they have a higher energy consumption rate. The time spent on the deliveries is shorter, which makes it easier to achieve time-precise deliveries and to reduce the costs of recharging, at the expense of higher costs due to fossil oil consumption.One of the first papers that dealt with the PHEVs routing problem was published by Abdallah [41], where the author defined the problem as the plug-in hybrid electric VRPTW (PHEVRPTW). The objective of the proposed problem is to minimize the routing costs on the internal combustion engine while satisfying the demand and time window constraints. At each customer, the driver can either recharge the vehicle battery or go to the next open time window using an internal combustion engine when the electric energy has been depleted. Doppstadt et al. [58] defined the HEV-TSP in which the authors observed HEVs that are not plug-in and can only be charged while driving. Four working modes were observed, combustion-only mode, electric-only mode, charging mode, and boost mode, when the combined internal combustion engine and electric engine are used. As the authors assumed that there was not a charge left from the day before, the initial charging of the battery was set to zero. The authors presented the positive effects of using HEVs: (i) compared to the ICEV routes, the overall costs for the HEV routes in the tested instances reduced up to 13% and driving time increased up to 11%; (ii) savings depend highly on the depot location, and not on the limited battery capacity, which was an a priori assumption. Mancini [71] introduced hybrid VRP (HVRP), in which PHEVs can change propulsion mode at any time and the electric engine, rather than the internal combustion engine, is promoted. Vincent et al. [76] introduced a similar HVRP problem with PHEVs and mild-HEVs, which do not have an electric-only mode of propulsion. The results showed that PHEVs have a lower average cost per mile in all scenarios compared to the mild-HEVs and ICEVs. Hiermann et al. [96], in their H2E-FTW, used the electric motor of the PHEV as priority mode choice. When a time window of a PHEV route is violated, the mode could be changed at any time during the route by exchanging recharging time for the corresponding amount of fuel at no additional costs. The authors showed the positive impact of a PHEV fleet compared to ICEV and BEV fleets, especially for randomly clustered instances. PHEVs are usually represented with only 20% of the overall number of vehicles used in the solution due to their higher consumption and utility costs but, still, they constitute an important part of the fleet configuration due to their flexibility.


### 4.4. Partial Recharging

In the beginning, in most of the E-VRP problems, full recharge was considered when a BEV visited a CS [25]. This can be time-consuming because, depending on the SoC level, available charging technology, and battery capacity, the vehicle can charge from five minutes to eight hours [15]. Therefore, in real-life applications, partial recharging should be taken into account. The battery should be charged enough to complete the whole route or to surpass a fear that vehicle range will not be enough to perform designated tasks, the so-called range anxiety [75]. This particularly has an effect on customers with narrow time windows, where efficient charge scheduling can enable the feasibility of the route. On the economic side, significant savings can be achieved by applying partial recharging as a minimal amount of energy could be recharged during the day when electricity cost and energy network load are higher, and the rest of the energy could be replenished during the night [46, 47]. In some cases, it is natural to maintain the energy reserve. This can be done by SoC range limitation, i.e., 20,95 [42, 47, 80]. Having an energy reserve seems even more important if energy consumption and range anxiety are taken into account because up to 30% of the consumed energy can be spent on BEV’s auxiliary devices [109]. Limiting SoC value also helps to preserve the battery as battery capacity decreases by overcharging and overdischarging [104].Several papers have analyzed strategies of partial recharging and formulated the problem as E-VRPTW with partial recharging (E-VRPTWPR). Bruglieri et al. [29, 51] and Moghaddam [53] modeled the concept of partial charging in E-VRP. Felipe et al. [46], in GRVP with multiple technologies and partial recharges (GVRP-MTPR), and Keskin and Çatay [28], in E-VRPTWPR, ensured that after every change in the route configuration, at the previous CS, the BEV is charged only with the amount of energy sufficient to finish the segment of the route until the next CS or the depot. This results in the removal of certain CSs in the route, compared to full recharge strategy, and the arrival of a vehicle at the depot with an empty battery. The authors presented the positive impact of partial recharges on the total costs and energy savings. A similar procedure was applied by Sassi et al. [47–49] but with multiple charging technologies, different charging periods, and BEV chargers compatibility checks. Schiffer and Walther [73] compared full and partial recharge by solving small E-VRPTWPR test instances using commercial software. The authors concluded that, in some cases, a partial recharging strategy reduces the total traveled distance and number of visits to CSs. Montoya [18] and Montoya et al. [72] formulated the FRVCP for BEVs in which, for a fixed sequence of customers in the route, CS position and charging amount are optimized. The results on the newly derived test instances showed that good solutions tend to exploit partial recharges. A similar procedure was applied by Keskin and Çatay [79], Hiermann et al. [96], Froger et al. [68], and Schiffer and Walther [88, 89] to enhance the incumbent best solution by optimizing the charging decisions along the BEV route with a fixed customer sequence. Desaulniers et al. [57] presented the effects of partial recharging by optimally solving E-VRPTW instances containing up to 100 customers. In a case with a single recharge per route, the partial recharge reduced the routing costs by 0.97% and the number of vehicles by 2.25%, while in a case with multiple recharges per route these values are 1.91% and 3.80%, respectively, with a significant increase in the average number of recharges per route.


### 4.5. Different Charging Technologies

Today, multiple charging technologies are present: (i) slow, 3 kW (6-8 h); (ii) fast, 7-43 kW (1-2 h); and (iii) rapid, 50-250 kW (5-30 min) [15, 129]. To better control charging time in the E-VRP context, the selection of possible charging technology could also be optimized. This could make some customers who have narrow time windows more accessible by fast charging at previous CSs, or if the time windows are long, an economically better approach could be slow charging. Such a problem could be extended by taking into account CS working hours, time-dependent charging costs, the number of available chargers and their compatibility with BEVs, the power grid load, the charger power, etc. [46–49, 68, 80, 81, 87]. Felipe et al. [46] analyzed the effect of different charging technologies on the recharge cost. The authors concluded that none of the technologies dominated the others. The rapid and fast charging options provided slightly better results than the slow charging, but the best results were obtained when joint technologies were used as the most appropriate one could be chosen in each case. Çatay and Keskin [67] solved small-sized instances to present the insights of the quick charge option. The results showed that quick charging might reduce the fleet size and decrease the cost of energy needed to operate BEVs. Keskin and Çatay [79] formulated E-VRPTW and fast charging (E-VRPTW-FC) with partial recharges and three different charging technologies. The authors provided MILP formulation and minimized total recharging costs while operating a minimum number of vehicles. The authors compared two mixed integer programming (MIP) models: (i) model 1, in which binary variables are used as decision variables to choose which charging technology to use; (ii) model 2, based on the model of Schneider et al. [25] and Keskin and Çatay [28], in which the authors copied each station vertex three times depending on the charging technology. The results showed that in all cases model 1 produced better results with shorter computation time. Testing on larger instances showed that the fast charging option is more beneficial when customers have narrow time windows and that multiple charging technologies improved the overall results (in 28 of 29 instances), while in the instances with larger time windows the average improvement was 0.17%.


### 4.6. Nonlinear Charging Function

In most of the E-VRP related literature, either linear or constant charging time is considered. Most of the BEVs have lithium-ion batteries installed, which are often charged in constant-current constant-voltage (CC-CV) phases: first by constant current until approx. 80% of the SoC value and then by a constant voltage. In the CC phase, SoC increases linearly, and in the CV phase, the current drops exponentially and SoC increases nonlinearly in time, which prolongs charging time [102, 104]. In only a few reviewed papers the author considered nonlinear charging process and solved the E-VRP either by linearization per segments or by estimating charging time by a data-driven approach [50, 68, 72, 87, 95]. Zündorf [50] developed a charging function propagation algorithm for determining an EV battery-constrained route by taking into account piecewise linear and concave functions of the recharging process. Montoya et al. [72] formulated the E-VRP with nonlinear charging functions (E-VRP-NL) and presented the MILP formulation. The authors fitted piecewise linear functions for 11, 22, and 44 kW CSs to the real measured data. The results showed that neglecting the nonlinear charge can lead to infeasible or overly expensive solutions: 12% of the routes in good solutions recharged the battery in the nonlinear part, after 80% of the SoC value. Froger et al. [95] proposed two new formulations for the E-VRP-NL: (i) the improved MILP version of Montoya et al. [72] by arc-based tracking of the SoC value; (ii) a path-based model without CS replication, where a sequence of vertices (a path) between a pair of customers or depots is determined. The second one yielded much better results with shorter computing time. Froger et al. [68] explored similar formulations for the E-VRP-NL with capacitated CSs (E-VRP-NL-C): (i) the CS replication-based formulation with flow and event-based formulations for CS capacity constraints; (ii) recharging path-based formulation.


### 4.7. BEV Routing and CS Location

Due to the currently low BEV market share, the number of CSs installed in the road infrastructure is also relatively low. Therefore, great potential lies in the simultaneous decision-making of CS locations and BEV routes. The classic location routing problem (LRP) consists of determining the locations of the depots and vehicle routes supplying customers from these depots [130]. A modification of the LRP that deals with CS facilities is formulated as electric LRP (E-LRP) [11, 27, 73, 88].Schiffer et al. [11, 27] presented a case study for solving ELRPTWPR. It was assumed that CSs can be located at customers’ locations and that service time can be used for recharging. Overall, results indicated the viability of combining CSs siting and BEV routing for specific cases when a delivery range is not far from the depot. Schiffer and Walther [73] compared ELRPTWPR and E-VRPTWPR solutions on down-scaled test instances and concluded that the ELRPTWPR gives a better solution in all instances. The reason comes from the fact that BEV can be charged while serving, which reduces the overall route time, as there is no travel time wasted on traveling to and from CS. Hof et al. [69] and Yang and Sun [56] addressed the problem of E-LRP but with BSSs. Results indicated that decreasing the construction cost of BSS led to the expected increase in their number in the solution. Schiffer and Walther [88] formulated the LRP with intraroute facilities (LRPIF) as a more general problem, in which intraroute facilities can be used for refueling, loading, or unloading goods. Sun et al. [131] proposed a location model for CSs without routing based on the travel demands of urban residents. For short-distance travelers, slow CSs are utilized, while for long-distance travelers, the rapid CSs are considered. The CSs’ locations are determined to maximize coverage and flow according to the concept of vehicle refueling. In the case study, the authors pointed out two critical aspects: the BEV driving range and the budget constraint.


### 4.8. Battery Swap

Instead of charging at CS, at specially designed BSS, empty or nearly empty batteries can be replaced with fully charged ones [14]. The main advantages of such procedure are the time in which it can be performed and the ability to recharge when energy network load and electricity costs are lower, i.e., during the night. A whole replacement procedure could last less than ten minutes, which is competitive to the refueling time of ICEVs and much faster than one of the fastest charging BEV technologies. The drawbacks of such procedure are the nonstandardized batteries and their installation in BEVs, which makes it hard to swap empty batteries with fully charged ones. Adler and Mirchandani [44] observed the E-VRP with swappable batteries, already determined BSS’s locations, a fixed number of batteries per BSS, full recharge of four hours, and a fixed swapping time of two minutes. Full battery recharge is considered so it is possible that when the vehicle arrives at the station, there is no fully charged battery available and the vehicle has to wait. First, the total routing cost is minimized and then the battery reservations are made, so that the vehicle could avoid BSSs without available batteries. Yang and Sun [56] and Hof et al. [69] simultaneously determined the locations of BSSs and the vehicle routing plan and provided different metaheuristics for solving the problem. Yang and Sun [56] formulated the problem as BSS-EV-LRP and compared the solutions of the BSS-EV-LRP instances to the corresponding best-known CVRP solutions. The total routing costs for 150 km and 300 km driving ranges increased to 7% and 3.5%, respectively, and in some cases, there was no gap between BSS-EV-LRP and CVRP solutions. Hof et al. [69] on the same problem reported that if construction costs of BSSs are equal to zero, the BSSs would be constructed in 83% of total available candidate sites.


### 4.9. Two-Echelon Routing Problem

Breunig et al. [93] proposed the electric two-echelon VRP (E2EVRP), in which goods are transported in two echelons: (i) in the first echelon, goods are transported by conventional freight vehicles from the depot to the satellite facilities; (ii) in the second echelon, goods are transported from the satellite facilities to the customers by light BEVs. Two vehicle types are observed in the problem: ICEVs with higher load capacity located at the depot and BEVs with lower load capacity located at satellite facilities. BEVs are used for the last-mile deliveries due to their lower pollution, noise, and size. The authors formulated the problem on the multigraph and applied exact procedures on smaller instances to solve the problem optimally and a metaheuristic procedure on larger newly developed instances. The results showed that the increase in CS density ρ decreased the detour costs following the expression 1/ρ1.24. Also, the vehicle range has a great impact on the solution quality: a range below 70 km produced infeasible solutions while a range higher than 150 km decreased recharging visits to almost zero. Jie et al. [97] analyzed a similar problem, the two-echelon capacitated E-VRP with BSS (2E-EVRP-BSS), in which, in both echelons, BEVs are used. The BEVs in the first echelon have higher load and battery capacities than the BEVs in the second echelon. The authors reported that the battery driving range is the most important aspect of routing and that it slightly depends on the number of BSSs used.


### 4.10. Charging Schedule

Many companies that use BEVs prefer charging the vehicles at their own facilities in order to charge the vehicles between the delivery routes and during specific periods of the day. In such occasions, there are usually a limited number of chargers at the depot, typically fewer than the fleet size; therefore, the efficient charging schedule at the depot has to be determined.Pelletier et al. [87] considered the electric freight vehicles charge scheduling problem (EFV-CSP) at the depot for the BEVs that deliver goods to a set of customers over multiple days. The authors included multiple charging technologies at the depot, realistic charging process (piecewise linearization of nonlinear charging function), time-dependent charging costs, grid power restriction, battery degradation costs (cyclic and calendar aging), and facility-related demand (FRD) charges representing the maximal demand registered over the billing period. The fixed vehicle routes are known in advance and the number of charging events in the depot is limited to avoid impractical solutions of constantly moving a vehicle from one charger to another. The authors presented the following aspects of charge scheduling tests conducted in summer (higher electricity costs) and winter (lower electricity costs) periods: (i) model tries to keep the SoC lower when battery degradation costs are included; (ii) in summer, vehicles are rarely charged in peak hours, which results in more vehicle charging simultaneously or the use of fast chargers that retrieve more power from the grid in nonpeak hours but incur higher FRD charges; (iii) to avoid cycling the battery in high SoC, it is preferable to split the long routes into smaller ones; (iv) fast chargers are heavily used in a high BEV utilization context; (v) grid power restriction increases overall energy costs, especially in summer months, and leads to infeasible solutions, limiting the number of vehicles simultaneously charging; and (vi) total costs are always lower with larger batteries, as smaller batteries require larger discharge cycles.Barco et al. [34, 42] also analyzed charge scheduling for assigned routes in an airport shuttle scenario. The authors defined the set of charging actions (charge profile) over the determined programming horizon and minimized overall operating costs. Sassi et al. [47, 49] also dealt with determining the charging schedule of BEVs at the depot and included time-dependent charging costs and chargers compatibility checks with BEVs. Adler and Mirchandani [44] provided an online routing model of BEVs by taking into account battery reservations to minimize the average delay of all vehicles by occasionally detouring them. Wen et al. [65] observed the service time of CSs, meaning that a CS can be visited only in some specific time period, usually working hours. The authors proposed a battery reservation model in order to charge the reserved battery at the defined time period. Sweda et al. [75] dealt with the possibility that a CS might not be available at some point in time and rerouting should be performed. The problem is formulated as the adaptive routing and recharging policies for EVs. First, optimal a priori routing and recharging policies were determined and then heuristic procedures were applied for the adaptive routing.


### 4.11. Dynamic Traffic Conditions

Most of the E-VRP research considers static conditions on the road network. The traffic states change recurrently, depending on the time of the day, day of the week, and season, or nonrecurrently when a traffic incident occurs, such as an accident [132–134]. TD-VRP routes a fleet of vehicles by taking into account variable travel time on the road network [135, 136]. Shao et al. [74] observed BEVs in such time-dependent context in their E-VRP with charging time and variable travel time (EVRP-CTVTT). The recharging time is fixed to 30 minutes and batteries are always charged to full capacity. The authors discretized one day into two-minute intervals and applied the dynamic Dijkstra algorithm to find the shortest travel time path when weights in the graph are not constant. The authors presented a real-life problem and solved it by applying a genetic algorithm with a running time of three hours. Omidvar and Tavakkoli-Moghaddam [24] presented a model for routing AFVs that aims to avoid routing during congestion hours, when the pollution costs are high, by waiting at a customer’s locations. These costs were computed based on the road speed profile. Mirmohammadi et al. [62] presented MILP formulation for the periodic green heterogeneous VRP with time-dependent urban traffic and time windows, which is green in terms of the emission minimization and not the utilization of AFVs. The planning horizon is divided into several periods, during which static traffic conditions are considered.


### 4.12. Related Problems

#### 4.12.1. Route Scheduling and Other Electric Variants

Many researchers are dealing with the scheduling of bus/taxi routes that are fixed or could be slightly altered. Li et al. [82] addressed the mixed bus fleet management problem (MBFM) composed of electric, diesel, compressed natural gas and hybrid-diesel buses. The authors maximized the total benefit of replacement of old vehicles with new ones under budget constraints while optimizing the route assignment for each bus during the planning period. Two routing procedures were developed to solve the recharging problem: single-period routing and routing across multiple periods of a day. Analysis of the results on the case study of Hong-Kong showed a cost-effective scheme of fleet configuration in the following order: mixed fleet, electric, natural gas, diesel, and then hybrid-diesel buses. Wen et al. [65] also dealt with the scheduling of electric buses (electric vehicle scheduling problem, E-VSP) in order to efficiently serve a set of timetabled bus trips. The linear partial and full recharge were considered, with the minimization of bus numbers and the total traveled distance. Lu et al. [83] addressed the problem of optimal scheduling of a taxi fleet with mixed BEVs and ICEVs to service advance reservations. The authors presented a multilayer taxi-flow time-space network in order to minimize the total operating fleet cost.Bruglieri et al. [77] formulated the MILP for a one-way electric car-sharing problem as the electric vehicle relocation problem (E-VReP). Relocation of BEVs is performed by workers who come with bicycles to the pickup points, put the bicycle into the BEV’s trunk, and drive the BEV to one of the delivery points. The authors balance a trade-off among the customers’ satisfaction, the workers’ workload balance, and the car-sharing provider’s objective.Masmoudi et al. [84] formulated the dial-a-ride problem with EVs and BSSs (DARP-EV) for customers with special needs and disabilities. The authors observed a fleet of heterogeneous BEVs with equal battery capacity but different available resources for the customers: handicapped person’s seat, stretcher, wheelchair, or accompanying person’s seat. The authors adapted the energy consumption model of Genikomsakis and Mitrentsis [85] with constant speed, acceleration, road slope, and varying mass while loading/unloading passengers.Paz et al. [86] addressed the multidepot electric vehicle location routing problem with time windows (MDEVLRPTW) with three variants: (i) with BSS (MDEVLRPTW-BS); (ii) with partial recharging (MDEVLRPTW-PR) at CSs that could also be located at customers’ locations; and (iii) with BSS and partial recharging (MDEVLRPTW-BSPR) where conventional charging is considered at customers’ locations and battery swap is performed at a separately located BSSs. Three MIP models were presented for each case. The authors determined the number and location of BSSs and/or CSs (without investment costs), the number and location of the depots, the number of vehicles, and the customers’ sequence in route. The proposed model allows minimizing the number of vehicles given a fixed number of CSs or BSSs, and vice versa. For the CS, linear recharge is considered, while the battery swap time is set to 10% of the corresponding linear recharge time. Results were presented by solving small instances with a commercial solver. The results showed that BSPR and PR variants provide better results than the BS variant and that BSPR tends to act as PR in most of the instances.Schiffer and Walther [89] extended the ELRPTWPR model of Schiffer and Walther [73] with uncertain customer patterns related to customers spatial distribution, demands, and service time windows. The authors formulated robust ELRPTWPR (RELRPTWPR) as a MILP program. Five scenarios were observed, which represent customer profile for five working days. The authors used set-based uncertainty representation, which is covered within scenarios that consider customer patterns in daily deliveries. To achieve a computationally tractable model of robust counterpart formulation, the adversarial approach was used, which uses a finite set of scenarios. The robust model was compared to the deterministic model (selected customer pattern), median model (average customer pattern), and worst-case model (worst-case customer pattern), where first the E-LRP was solved, and then the E-VRP. The deterministic models produced infeasible solutions in 35-63% of instances and incurred, on average, 5% higher costs than the robust approach. Also, the robust approach showed the most homogeneous configuration of CSs.Keskin et al. [80] extended E-VRPTW by considering waiting times at the CSs and proposed a solution method for solving small-sized instances. The authors assume single charger at CS, Poisson arrivals, exponential distribution of service times, first-in-first-out strategy, and penalties for late arrivals. The planning horizon is split into morning, afternoon, evening, and night period. The routing decisions are determined based on time-dependent waiting times at CSs.Kullman et al. [81] introduced E-VRP with public-private recharge strategy (E-VRP-PP), where demand at public CSs is unknown. The queuing at CS is modeled as M/M/c with first-in-first-out strategy, and c represents the number of identical chargers at CS. The problem is modeled as a Markov decision process, and an approximate dynamic programming solution is proposed. The authors extended the FRVCP formulation of Montoya et al. [72] with time-dependent waiting times and discrete charging decisions. Dynamic and static routing policies are proposed, where the static policies are incorporated as base policies for dynamic routing. Using dual bound on the optimal policy, the authors concluded that proposed routing policies are within 4.7%. The authors pointed out that the use of a public-private strategy outperformed the use of just private CSs, as overall savings were 20% higher.


#### 4.12.2. Green and Pollution Routing Problem

The green vehicle routing problem focuses on the reduction of routing pollution on the environment. The key idea is to promote the use of sustainable energy sources and minimize overall emissions. Regarding the emission of BEV, Álvarez Fernández [137] indicated that emission savings of BEVs can be significant, up to 80% in some scenarios compared to ICEVs. The authors presented a model for prediction of a BEV’s GHG emissions linked to a route for each country, according to empirical measurements and prediction of routes’ energy consumption.Erdoğan and Miller-Hooks [23] first defined the GVRP with AFVs and full refuel, where also BEVs were included. Therefore, the E-VRP can be observed as a special case of the GVRP. The authors provided MILP formulation of the problem with replication of CSs vertices. Koç and Karaoglan [59] improved the MILP formulation of Erdoğan and Miller-Hooks [23] by removing dummy vertices that represent multiple visits to refueling stations. Instead, a binary variable related to traveling between two vertices by using the station in between is introduced (arc-duplicating). This formulation improved the average gap to the optimal integer solution on all scenarios compared to the model of Erdoğan and Miller-Hooks [23] by 25%. The authors proposed a combination of exact branch-and-cut procedure and simulated annealing with several moving strategies to solve the problem. Leggieri and Haouari [70] formulated the problem as a nonlinear MIP, which was linearized to derive MILP formulation by applying a reformulation-linearization technique. The authors used similar arc-duplicating to that of Koç and Karaoglan [59] but defined the binary variable as the sum of binary variables for the direct arc between the customers and binary variables for the same direct arc but with stations inserted. The proposed formulation without additional constraints was able to solve to optimality 97.5% of test instances, which substantially outperformed the formulations of Erdoğan and Miller-Hooks [23] and Koç and Karaoglan [59]. In addition, the same formulation with prepossessing consistently outperformed the exact branch-and-cut procedure of Koç and Karaoglan [59]. Andelmin and Bartolini [66] took a different approach for the exact solving of GVRP based on a set partitioning formulation in which columns correspond to feasible routes. The problem is modeled using a multigraph in which vertices correspond to customers and each arc represents a possible sequence of stations visited when traveling between customers. The authors added weak subset row inequalities, subset row inequalities, and k-path cuts, with the latter reported as the crucial part of the algorithm efficiency. The results showed tight lower bounds and high efficiency of the algorithm as instances with up to 110 customers were solved to optimality. Bruglieri et al. [94] presented a path-based approach for solving the GVRP, based on two phases. In the first phase, a set of feasible paths is generated, removing from it all those paths that are dominated. In the second phase, through MILP formulation, paths are combined to create routes for the final solution. The authors compared their results to previously observed research papers of Erdoğan and Miller-Hooks [23], Koç and Karaoglan [59], Leggieri and Haouari [70], and Andelmin and Bartolini [66]. On small-sized instances, the proposed algorithm produced optimal solutions with the lowest running time while on the larger instances it produced 2.3% worse solutions compared to Andelmin and Bartolini [66], but in much less time.Poonthalir and Nadarajan [38] formulated biobjective fuel efficient GVRP (F-GVRP) with the minimization of route costs and fuel consumption using goal programming. The authors analyzed the impact of varying speed on fuel consumption where speed is modeled using triangular distribution. To efficiently solve the problem, the authors applied a population-based metaheuristic. The tests were performed by constraining the route cost of the solution to the BKS of GVRP and then the fuel minimization was performed. The authors pointed out that the use of varying speed instead of constant speed reduced the fuel consumption on average by 20.48%, and in half of the instances it reduced the number of vehicles. The lower number of vehicles is a consequence of the lesser amount of time spent at the refueling stations as less fuel is consumed, so the vehicle can visit more customers. The authors pointed out that if the route cost constraint was relaxed, most of the instances might have better route costs than the BKSs of GVRP. The authors also analyzed the impact of speed interval size on the fuel consumption and concluded that fuel consumption and average speed increased with larger interval size.Normasari et al. [100] formulated capacitated GVRP (CGVRP) as a MILP program and solved the problem by applying a simulated annealing heuristic. Compared to the GVRP, the results indicated that the number of tours and total costs for the CGVRP increased by 38.8% and 32.5%, respectively.Macrina et al. [33] investigated the green mixed fleet VRP with partial recharges and time windows (GMFVRP-PRTW), in which polluting emissions are modeled through the functions of both traveled distance and vehicle load. The problem was formulated as a MILP and solved by a commercial solver for small instances and by an iterative local search metaheuristic for large instances. The overall pollution emissions are maintained below the designated threshold value, to simulate emission limits of freight vehicles. Macrina et al. [99] presented MILP formulation for a similar GMFVRP-PRTW problem but added more realistic consumption model presented in Section 3.2, where the authors considered different charging technologies (single technology per CS), cost of energy recharged at the depot, and the fuel cost. A large neighborhood search metaheuristic is applied to solve the problem.Koyuncu and Yavuz [98] compared the node- and arc-duplicating MILP formulation of mixed fleet GVRP (MGVRP) by taking into account the following: (i) a mixed fleet of AFVs and ICEVs; (ii) load and time window constraints; (iii) internal station visits (at customer) and external station visit (at separately located stations); (iv) fixed full, variable full, or variable partial refuel at the station; (v) site-dependent refueling rate; and (vi) the visit of two or more consecutive stations being prohibited. In the node-duplicating formulations [23, 25], each refueling vertex is replicated at most N times, where N is the number of customers. In the arc-duplicating formulation, two subsets of arcs are contained: arcs between customers (direct arcs) and refueling arcs which contain inserted station between customers [70]. In the preprocessing step for both formulations, unfeasible arcs are removed and tighter lower and upper bounds for auxiliary variables and number of vehicles are presented. On the GVRP test instances of Erdoğan and Miller-Hooks [23], the authors reported that the node- and arc-duplicating formulations yielded an average optimality gap of 20.84% and 9.07%, respectively. On the GVRP instances of Yavuz [138], significant conclusions could not be drawn.As BEVs have no local emissions, the E-VRP is closely related to the minimization of GHG emissions, where a problem-specific GVRP variant called the pollution routing problem (PRP) was introduced [32, 40]. Bektaş and Laporte [40] formulated the PRP as a nonlinear MIP where the main objective is the optimization of vehicle speed and GHG emissions. Demir et al. [32] concluded that speed optimization improves the PRP solution and minimizes fuel consumption and driver costs.It can be noted that several papers are dealing with fuel consumption and pollution emissions and not AFVs, but the termgreen VRP is used to indicate a more ecologically aware routing problem [62, 139].


#### 4.12.3. More General Formulation of the E-VRP

The E-VRP can be formulated more generally as the VRP with intermediate stops (VRPIS), in which the vehicle visits intermediate or intraroute facilities to replenish/unload the goods or to refuel [55]. The LRPIF extends the problem in terms of the facilities location optimization [88]. The service time depends on the load and fuel level at the arrival at the facility. Schneider et al. [55] presented the E-VRP with recharging facilities (EVRPRF) as a special case of VRPIS, where BEVs have a fixed recharging time and limited vehicle load capacity without time windows constraints.


#### 4.12.4. Multiobjective Variants

Most of the researchers are using a single-objective function that represents travel distance, total costs, energy consumption, total time, etc. The multiobjective variants are still relatively scarcely applied. Generally, the solution of the multiobjective problem is not the optimal solution for all of the objectives, but rather it is satisfactory in those terms.Demir et al. [37] introduced the biobjective PRP with the minimization of fuel consumption and total driving time. The fuel consumption is based on a model similar to the one presented by (6) but with constant speed. An adaptive large neighborhood search was used to generate nondominated/Pareto optimal solutions while an additional hybrid procedure was used to obtain the solution to the biobjective model. The experiments showed that there is no need to prolong driving time to achieve a significant reduction in fuel consumption.The already mentioned F-GVRP of Poonthalir and Nadarajan [38] (Section 4.12.2) used biobjective minimization of route costs and fuel consumption. Wang et al. [91] considered multiobjective function to determine optimal BEV routes between origin-destination pairs. Three objectives are minimized: travel time, energy consumption, and charging cost. Travel time includes driving time, queuing time, and charging time, while charging cost includes electricity cost, service cost, and parking fee. The multiobjective model is transformed into a single objective by applying a fuzzy programming approach for objectives’ membership functions and fuzzy preference relations to obtain weighting coefficients of each objective. A genetic algorithm is designed to obtain an optimal solution. The influences of various driver trade-offs (weighting conditions) among different objectives are explored using Pareto curves. It was observed that driver trade-offs have an effect on the travel times, and almost no effect on the energy consumption and charging costs.Amiri et al. [39] optimized the location of the battery swapping in the BSSs network and the charging schedule of depleted batteries. The multiobjective model minimizes three objectives: battery charging and power loss costs, deviation from the nominal voltages, and network capacity releasing. The nondominated sorting genetic algorithm is applied to solve the problem. The results on the test case showed that an optimal charging schedule in terms of cost and network constraints can be achieved. Besides, the electricity load profile is more leveled.Bruglieri et al. [77] focused on the E-VReP with three objectives: minimization of workers used to relocate the BEVs (service provider objective), maximization of relocation requests served (users’ objective), and minimization of the longest routes’ duration (workload balance). The problem is solved by approximating the Pareto optimal front based on a two-phase approach: first, feasible solutions are generated by randomized search heuristics, and in the second phase, to find a set of nondominated solutions, the MILP formulation is solved by epsilon constraint programming.


### 4.13. Summary of the E-VRP Variants

Table1 presents a summary of the E-VRP variants and related problems in the available literature. In total, 79 papers were presented in the table. The problem characteristics are observed in the following order: vehicle load (cargo) capacity (48), linear charging/refueling function at station (34), customer time windows (34), partial recharging strategy (30), fixed refueling/recharging time (19), different charging technologies at CSs (17), heterogeneous (mixed) fleet (16), energy consumption model (15), BSS (9), nonlinear charging function (8), dynamic traffic conditions (8), location routing (7), capacitated stations (6), and HEVs (5). We can point out that only a few papers are dealing with (i) simultaneous CS siting and BEV routing (E-LRP) where great savings can be achieved in companies whose business plans include CS investments, as E-LRP generally gives better routing configuration than the E-VRP [73, 88]; (ii) nonlinear charging function, negligence of which can lead to infeasible solutions and additional penalty costs [68, 72]; (iii) HEVs (PHEVs) that have an advantage in several routing scenarios but are being researched only recently due to the higher complexity of the problem as HEVs can change their propulsion mode at any time during the route [96]; (iv) dynamic traffic conditions, which significantly influence BEVs’ energy consumption [74]; (v) different charging technologies at CSs, which could reduce the overall costs by selecting the best possible charger option in each occasion [46, 79]; (vi) CS related challenges: capacitated CS [68, 75, 80], CS reservations [44], and public-private recharge strategy [81]; and (vii) robust E-VRP variant resistant to uncertainties in travel time, service time, and demand [89].


## 4.1. Energy Shortest Path Problem and the Electric Traveling Salesman Problem

The energy shortest path problem (E-SPP) and electric TSP (E-TSP) can be considered as two of the simplest forms of the E-VRP. In most of the VRP variants, graph arcs have positive weight values that can represent distance, travel time, cost, etc. To compute the shortest path between customers, the most commonly applied algorithms are Dijkstra, Bellman-Ford,A∗, contraction hierarchies, etc. [108, 120]. By taking into account the recuperated energy of BEV, some arc weights could have a negative value, which makes most of the shortest path algorithms inapplicable. To overcome this problem, Artmeier et al. [108] formalized energy-efficient routing with rechargeable batteries as a special case of the constrained SPP (CSPP) in which the battery charge is limited with its capacity, and there is no recharging at CS. The authors adapted the Bellman-Ford shortest path algorithm with time complexity O(n3) in order to solve the energy CSPP. To overcome negative graph edges, Eisner et al. [121] used Johnson’s shifting technique [122] to transform negative edge cost functions into nonnegative ones and applied Dijsktra’s algorithm with time complexity of O(nlog(n)+m) [123]. Storandt [124] introduced CSPP for EVs with battery swapping, while Sweda and Klabjan [125] added recharging events at CSs with a differentiable charging function. Zündorf [50] solved the CSPP for BEV routing with battery constraints, different types of CSs, and nonlinear charging process. The author developed a charging function propagating algorithm in order to minimize travel time and applied CH and A∗ algorithms to solve the problem. Liao et al. [103] considered the BEV’s shortest travel path problem and presented a dynamic programming algorithm for solving the problem that runs in O(kn2), where k is the upper bound on the number of BSSs. More recently, Strehler et al. [126] observed recharging events in the graph vertices and arcs for the energy-efficient shortest routes of BEVs and HEVs, while Zhang et al. [127] dealt with the electric vehicle route planning with recharging problem (EVRC), which minimizes the overall travel and charging time and takes into account CSs’ locations, partial recharging, nonlinear charging functions, service time duration, and service frequency at CS.As VRP is the generalization of the TSP, the E-VRP is closely related to the E-TSP, in which a set of customers has to be served by only one BEV. Roberti and Wen [64] formulated the E-TSP with time windows (E-TSPTW) as a compact MILP program with binary variables for recharging paths with intermediate stops. The authors observed full and partial recharge policies and applied three-phase heuristic for solving the problem. The authors solved multiple instances, among them the 13 E-VRPTW instances of Schneider et al. [25], which were solved optimally with only one vehicle. Doppstadt et al. [58] formulated the HEV-TSP with four working modes of HEV and provided new test instances. Liao et al. [103] introduced the EV touring problem, as the generalization of TSP, with the minimization of the total required time. Two scenarios with battery swap scheme were observed: the on-site station model, in which each city has a BSS; the off-site station model, in which BSS is located at an acceptable distance from the city. The authors proposed efficient polynomial time algorithms for the problem.


## 4.2. Heterogeneous or Mixed Vehicle Fleet

In today’s vehicle fleets, mostly ICEVs are present. Transition to an almost wholly electric fleet is a very challenging economic task. Therefore, most companies are gradually integrating BEVs into their existing ICEV fleet. Routing algorithms have to be upgraded and adapted for electric fleet characteristics, as on-time priority is harder to achieve.Fleet size and mix VRP (FSM-VRP) was first introduced by Golden et al. [128], where the authors considered routing a fleet of vehicles with different acquisition costs and the routing cost is dependent on the vehicle type. Goeke and Schneider [30] formulated E-VRPTW and mixed fleet (E-VRPTWMF) with equal vehicle load capacities of BEVs and ICEVs, and equal BEV battery capacities. The authors compared the percentage share of BEVs to the total traveled distance obtained with three objective functions: total traveled distance, total costs with battery costs, and total costs without battery costs. The traveled distance on one battery pack was assumed to be 241,350 km (150,000 miles) with the payment of an additional $600 per kWh in case of replacement. The results showed that the BEV’s share in total traveled distance increased significantly when battery costs are not considered, while in the other two scenarios ICEVs performed most of the deliveries. Hiermann et al. [31] analyzed a similar problem, the electric fleet size and mix VRPTW and recharging stations (E-FSMFTW), with different vehicle load and battery capacities. The authors showed the overall positive influence of the heterogeneous vehicle fleet on the generalized total cost function. In most of the instances, three to four vehicle types are used in the final solution. A similar analysis on small-sized instances with different BEVs, ICEVs, and HEVs was undertaken by Lebeau et al. [52]. The authors defined seven groups of vehicle types that could be used for the delivery, from small vans and quadricycles through diesel-only or electric-only groups to a group of all vehicle types. The results showed the following aspects of routing: (i) the fleet with different vehicle types reduced the total routing costs; (ii) in the large van group, ICEVs outperformed BEVs; and (iii) HEVs showed a great application in deliveries solely made by trucks. Sassi et al. [47–49] also observed a heterogeneous fleet of ICEVs and BEVs with different load and battery capacities, different operating costs, time-dependent charges, and the compatibility of BEVs and chargers at CSs. The authors focused on the procedures for solving the problem, so they did not offer any comparison between the ICEV and BEV solutions. Hiermann et al. [96] introduced the Hybrid heterogeneous electric fleet routing problem with time windows and recharging stations (H2E-FTW), in which the fleet consists of ICEVs, BEVs, and PHEVs. The authors compared the overall costs of the solutions obtained with a homogeneous fleet of ICEVs, BEVs, and PHEVs to the optimized solution with a mixed fleet. The gap for the ICEV fleet is the largest, with an average value of 60% in the extreme case when the fuel cost is the highest. In contrast, when the electric cost is the highest, the BEV fleet on average produces a 40% gap, while the PHEV fleet shows the lowest gap value, up to 25%. The authors pointed out that, in mixed solutions, BEVs are preferred for clustered instances, ICEVs for randomly distributed instances, and PHEVs for randomly clustered and distributed instances. Overall, the results showed that operational costs can be 7% lower in a mixed case compared to the homogeneous case.


## 4.3. Hybrid Vehicles

As compensation for the limited driving range of BEVs, HEVs, which have both an internal combustion engine and an electric engine, have been developed. Two main types of HEVs are present on the market: theseries hybrid, in which only the electric motor drives the train and the internal combustion engine is used to recharge the battery pack, and theparallel hybrid, which uses both internal combustion engine and electric engine to drive the train, where the electric engine is more efficient in stop-and-go activities and the internal combustion engine is more efficient at high speeds. We focus on the PHEV as a version of aparallel hybrid in which batteries can be recharged by connecting a plug to the electric power source. The PHEVs have an option to decide during the route to run on either electric energy or fossil oil. This enables the visiting of customers far from the depot with almost no refuel/recharge during the route. As PHEVs have two engines, their load is heavier compared to the BEVs and ICEVs, and therefore they have a higher energy consumption rate. The time spent on the deliveries is shorter, which makes it easier to achieve time-precise deliveries and to reduce the costs of recharging, at the expense of higher costs due to fossil oil consumption.One of the first papers that dealt with the PHEVs routing problem was published by Abdallah [41], where the author defined the problem as the plug-in hybrid electric VRPTW (PHEVRPTW). The objective of the proposed problem is to minimize the routing costs on the internal combustion engine while satisfying the demand and time window constraints. At each customer, the driver can either recharge the vehicle battery or go to the next open time window using an internal combustion engine when the electric energy has been depleted. Doppstadt et al. [58] defined the HEV-TSP in which the authors observed HEVs that are not plug-in and can only be charged while driving. Four working modes were observed, combustion-only mode, electric-only mode, charging mode, and boost mode, when the combined internal combustion engine and electric engine are used. As the authors assumed that there was not a charge left from the day before, the initial charging of the battery was set to zero. The authors presented the positive effects of using HEVs: (i) compared to the ICEV routes, the overall costs for the HEV routes in the tested instances reduced up to 13% and driving time increased up to 11%; (ii) savings depend highly on the depot location, and not on the limited battery capacity, which was an a priori assumption. Mancini [71] introduced hybrid VRP (HVRP), in which PHEVs can change propulsion mode at any time and the electric engine, rather than the internal combustion engine, is promoted. Vincent et al. [76] introduced a similar HVRP problem with PHEVs and mild-HEVs, which do not have an electric-only mode of propulsion. The results showed that PHEVs have a lower average cost per mile in all scenarios compared to the mild-HEVs and ICEVs. Hiermann et al. [96], in their H2E-FTW, used the electric motor of the PHEV as priority mode choice. When a time window of a PHEV route is violated, the mode could be changed at any time during the route by exchanging recharging time for the corresponding amount of fuel at no additional costs. The authors showed the positive impact of a PHEV fleet compared to ICEV and BEV fleets, especially for randomly clustered instances. PHEVs are usually represented with only 20% of the overall number of vehicles used in the solution due to their higher consumption and utility costs but, still, they constitute an important part of the fleet configuration due to their flexibility.


## 4.4. Partial Recharging

In the beginning, in most of the E-VRP problems, full recharge was considered when a BEV visited a CS [25]. This can be time-consuming because, depending on the SoC level, available charging technology, and battery capacity, the vehicle can charge from five minutes to eight hours [15]. Therefore, in real-life applications, partial recharging should be taken into account. The battery should be charged enough to complete the whole route or to surpass a fear that vehicle range will not be enough to perform designated tasks, the so-called range anxiety [75]. This particularly has an effect on customers with narrow time windows, where efficient charge scheduling can enable the feasibility of the route. On the economic side, significant savings can be achieved by applying partial recharging as a minimal amount of energy could be recharged during the day when electricity cost and energy network load are higher, and the rest of the energy could be replenished during the night [46, 47]. In some cases, it is natural to maintain the energy reserve. This can be done by SoC range limitation, i.e., 20,95 [42, 47, 80]. Having an energy reserve seems even more important if energy consumption and range anxiety are taken into account because up to 30% of the consumed energy can be spent on BEV’s auxiliary devices [109]. Limiting SoC value also helps to preserve the battery as battery capacity decreases by overcharging and overdischarging [104].Several papers have analyzed strategies of partial recharging and formulated the problem as E-VRPTW with partial recharging (E-VRPTWPR). Bruglieri et al. [29, 51] and Moghaddam [53] modeled the concept of partial charging in E-VRP. Felipe et al. [46], in GRVP with multiple technologies and partial recharges (GVRP-MTPR), and Keskin and Çatay [28], in E-VRPTWPR, ensured that after every change in the route configuration, at the previous CS, the BEV is charged only with the amount of energy sufficient to finish the segment of the route until the next CS or the depot. This results in the removal of certain CSs in the route, compared to full recharge strategy, and the arrival of a vehicle at the depot with an empty battery. The authors presented the positive impact of partial recharges on the total costs and energy savings. A similar procedure was applied by Sassi et al. [47–49] but with multiple charging technologies, different charging periods, and BEV chargers compatibility checks. Schiffer and Walther [73] compared full and partial recharge by solving small E-VRPTWPR test instances using commercial software. The authors concluded that, in some cases, a partial recharging strategy reduces the total traveled distance and number of visits to CSs. Montoya [18] and Montoya et al. [72] formulated the FRVCP for BEVs in which, for a fixed sequence of customers in the route, CS position and charging amount are optimized. The results on the newly derived test instances showed that good solutions tend to exploit partial recharges. A similar procedure was applied by Keskin and Çatay [79], Hiermann et al. [96], Froger et al. [68], and Schiffer and Walther [88, 89] to enhance the incumbent best solution by optimizing the charging decisions along the BEV route with a fixed customer sequence. Desaulniers et al. [57] presented the effects of partial recharging by optimally solving E-VRPTW instances containing up to 100 customers. In a case with a single recharge per route, the partial recharge reduced the routing costs by 0.97% and the number of vehicles by 2.25%, while in a case with multiple recharges per route these values are 1.91% and 3.80%, respectively, with a significant increase in the average number of recharges per route.


## 4.5. Different Charging Technologies

Today, multiple charging technologies are present: (i) slow, 3 kW (6-8 h); (ii) fast, 7-43 kW (1-2 h); and (iii) rapid, 50-250 kW (5-30 min) [15, 129]. To better control charging time in the E-VRP context, the selection of possible charging technology could also be optimized. This could make some customers who have narrow time windows more accessible by fast charging at previous CSs, or if the time windows are long, an economically better approach could be slow charging. Such a problem could be extended by taking into account CS working hours, time-dependent charging costs, the number of available chargers and their compatibility with BEVs, the power grid load, the charger power, etc. [46–49, 68, 80, 81, 87]. Felipe et al. [46] analyzed the effect of different charging technologies on the recharge cost. The authors concluded that none of the technologies dominated the others. The rapid and fast charging options provided slightly better results than the slow charging, but the best results were obtained when joint technologies were used as the most appropriate one could be chosen in each case. Çatay and Keskin [67] solved small-sized instances to present the insights of the quick charge option. The results showed that quick charging might reduce the fleet size and decrease the cost of energy needed to operate BEVs. Keskin and Çatay [79] formulated E-VRPTW and fast charging (E-VRPTW-FC) with partial recharges and three different charging technologies. The authors provided MILP formulation and minimized total recharging costs while operating a minimum number of vehicles. The authors compared two mixed integer programming (MIP) models: (i) model 1, in which binary variables are used as decision variables to choose which charging technology to use; (ii) model 2, based on the model of Schneider et al. [25] and Keskin and Çatay [28], in which the authors copied each station vertex three times depending on the charging technology. The results showed that in all cases model 1 produced better results with shorter computation time. Testing on larger instances showed that the fast charging option is more beneficial when customers have narrow time windows and that multiple charging technologies improved the overall results (in 28 of 29 instances), while in the instances with larger time windows the average improvement was 0.17%.


## 4.6. Nonlinear Charging Function

In most of the E-VRP related literature, either linear or constant charging time is considered. Most of the BEVs have lithium-ion batteries installed, which are often charged in constant-current constant-voltage (CC-CV) phases: first by constant current until approx. 80% of the SoC value and then by a constant voltage. In the CC phase, SoC increases linearly, and in the CV phase, the current drops exponentially and SoC increases nonlinearly in time, which prolongs charging time [102, 104]. In only a few reviewed papers the author considered nonlinear charging process and solved the E-VRP either by linearization per segments or by estimating charging time by a data-driven approach [50, 68, 72, 87, 95]. Zündorf [50] developed a charging function propagation algorithm for determining an EV battery-constrained route by taking into account piecewise linear and concave functions of the recharging process. Montoya et al. [72] formulated the E-VRP with nonlinear charging functions (E-VRP-NL) and presented the MILP formulation. The authors fitted piecewise linear functions for 11, 22, and 44 kW CSs to the real measured data. The results showed that neglecting the nonlinear charge can lead to infeasible or overly expensive solutions: 12% of the routes in good solutions recharged the battery in the nonlinear part, after 80% of the SoC value. Froger et al. [95] proposed two new formulations for the E-VRP-NL: (i) the improved MILP version of Montoya et al. [72] by arc-based tracking of the SoC value; (ii) a path-based model without CS replication, where a sequence of vertices (a path) between a pair of customers or depots is determined. The second one yielded much better results with shorter computing time. Froger et al. [68] explored similar formulations for the E-VRP-NL with capacitated CSs (E-VRP-NL-C): (i) the CS replication-based formulation with flow and event-based formulations for CS capacity constraints; (ii) recharging path-based formulation.


## 4.7. BEV Routing and CS Location

Due to the currently low BEV market share, the number of CSs installed in the road infrastructure is also relatively low. Therefore, great potential lies in the simultaneous decision-making of CS locations and BEV routes. The classic location routing problem (LRP) consists of determining the locations of the depots and vehicle routes supplying customers from these depots [130]. A modification of the LRP that deals with CS facilities is formulated as electric LRP (E-LRP) [11, 27, 73, 88].Schiffer et al. [11, 27] presented a case study for solving ELRPTWPR. It was assumed that CSs can be located at customers’ locations and that service time can be used for recharging. Overall, results indicated the viability of combining CSs siting and BEV routing for specific cases when a delivery range is not far from the depot. Schiffer and Walther [73] compared ELRPTWPR and E-VRPTWPR solutions on down-scaled test instances and concluded that the ELRPTWPR gives a better solution in all instances. The reason comes from the fact that BEV can be charged while serving, which reduces the overall route time, as there is no travel time wasted on traveling to and from CS. Hof et al. [69] and Yang and Sun [56] addressed the problem of E-LRP but with BSSs. Results indicated that decreasing the construction cost of BSS led to the expected increase in their number in the solution. Schiffer and Walther [88] formulated the LRP with intraroute facilities (LRPIF) as a more general problem, in which intraroute facilities can be used for refueling, loading, or unloading goods. Sun et al. [131] proposed a location model for CSs without routing based on the travel demands of urban residents. For short-distance travelers, slow CSs are utilized, while for long-distance travelers, the rapid CSs are considered. The CSs’ locations are determined to maximize coverage and flow according to the concept of vehicle refueling. In the case study, the authors pointed out two critical aspects: the BEV driving range and the budget constraint.


## 4.8. Battery Swap

Instead of charging at CS, at specially designed BSS, empty or nearly empty batteries can be replaced with fully charged ones [14]. The main advantages of such procedure are the time in which it can be performed and the ability to recharge when energy network load and electricity costs are lower, i.e., during the night. A whole replacement procedure could last less than ten minutes, which is competitive to the refueling time of ICEVs and much faster than one of the fastest charging BEV technologies. The drawbacks of such procedure are the nonstandardized batteries and their installation in BEVs, which makes it hard to swap empty batteries with fully charged ones. Adler and Mirchandani [44] observed the E-VRP with swappable batteries, already determined BSS’s locations, a fixed number of batteries per BSS, full recharge of four hours, and a fixed swapping time of two minutes. Full battery recharge is considered so it is possible that when the vehicle arrives at the station, there is no fully charged battery available and the vehicle has to wait. First, the total routing cost is minimized and then the battery reservations are made, so that the vehicle could avoid BSSs without available batteries. Yang and Sun [56] and Hof et al. [69] simultaneously determined the locations of BSSs and the vehicle routing plan and provided different metaheuristics for solving the problem. Yang and Sun [56] formulated the problem as BSS-EV-LRP and compared the solutions of the BSS-EV-LRP instances to the corresponding best-known CVRP solutions. The total routing costs for 150 km and 300 km driving ranges increased to 7% and 3.5%, respectively, and in some cases, there was no gap between BSS-EV-LRP and CVRP solutions. Hof et al. [69] on the same problem reported that if construction costs of BSSs are equal to zero, the BSSs would be constructed in 83% of total available candidate sites.


## 4.9. Two-Echelon Routing Problem

Breunig et al. [93] proposed the electric two-echelon VRP (E2EVRP), in which goods are transported in two echelons: (i) in the first echelon, goods are transported by conventional freight vehicles from the depot to the satellite facilities; (ii) in the second echelon, goods are transported from the satellite facilities to the customers by light BEVs. Two vehicle types are observed in the problem: ICEVs with higher load capacity located at the depot and BEVs with lower load capacity located at satellite facilities. BEVs are used for the last-mile deliveries due to their lower pollution, noise, and size. The authors formulated the problem on the multigraph and applied exact procedures on smaller instances to solve the problem optimally and a metaheuristic procedure on larger newly developed instances. The results showed that the increase in CS density ρ decreased the detour costs following the expression 1/ρ1.24. Also, the vehicle range has a great impact on the solution quality: a range below 70 km produced infeasible solutions while a range higher than 150 km decreased recharging visits to almost zero. Jie et al. [97] analyzed a similar problem, the two-echelon capacitated E-VRP with BSS (2E-EVRP-BSS), in which, in both echelons, BEVs are used. The BEVs in the first echelon have higher load and battery capacities than the BEVs in the second echelon. The authors reported that the battery driving range is the most important aspect of routing and that it slightly depends on the number of BSSs used.


## 4.10. Charging Schedule

Many companies that use BEVs prefer charging the vehicles at their own facilities in order to charge the vehicles between the delivery routes and during specific periods of the day. In such occasions, there are usually a limited number of chargers at the depot, typically fewer than the fleet size; therefore, the efficient charging schedule at the depot has to be determined.Pelletier et al. [87] considered the electric freight vehicles charge scheduling problem (EFV-CSP) at the depot for the BEVs that deliver goods to a set of customers over multiple days. The authors included multiple charging technologies at the depot, realistic charging process (piecewise linearization of nonlinear charging function), time-dependent charging costs, grid power restriction, battery degradation costs (cyclic and calendar aging), and facility-related demand (FRD) charges representing the maximal demand registered over the billing period. The fixed vehicle routes are known in advance and the number of charging events in the depot is limited to avoid impractical solutions of constantly moving a vehicle from one charger to another. The authors presented the following aspects of charge scheduling tests conducted in summer (higher electricity costs) and winter (lower electricity costs) periods: (i) model tries to keep the SoC lower when battery degradation costs are included; (ii) in summer, vehicles are rarely charged in peak hours, which results in more vehicle charging simultaneously or the use of fast chargers that retrieve more power from the grid in nonpeak hours but incur higher FRD charges; (iii) to avoid cycling the battery in high SoC, it is preferable to split the long routes into smaller ones; (iv) fast chargers are heavily used in a high BEV utilization context; (v) grid power restriction increases overall energy costs, especially in summer months, and leads to infeasible solutions, limiting the number of vehicles simultaneously charging; and (vi) total costs are always lower with larger batteries, as smaller batteries require larger discharge cycles.Barco et al. [34, 42] also analyzed charge scheduling for assigned routes in an airport shuttle scenario. The authors defined the set of charging actions (charge profile) over the determined programming horizon and minimized overall operating costs. Sassi et al. [47, 49] also dealt with determining the charging schedule of BEVs at the depot and included time-dependent charging costs and chargers compatibility checks with BEVs. Adler and Mirchandani [44] provided an online routing model of BEVs by taking into account battery reservations to minimize the average delay of all vehicles by occasionally detouring them. Wen et al. [65] observed the service time of CSs, meaning that a CS can be visited only in some specific time period, usually working hours. The authors proposed a battery reservation model in order to charge the reserved battery at the defined time period. Sweda et al. [75] dealt with the possibility that a CS might not be available at some point in time and rerouting should be performed. The problem is formulated as the adaptive routing and recharging policies for EVs. First, optimal a priori routing and recharging policies were determined and then heuristic procedures were applied for the adaptive routing.


## 4.11. Dynamic Traffic Conditions

Most of the E-VRP research considers static conditions on the road network. The traffic states change recurrently, depending on the time of the day, day of the week, and season, or nonrecurrently when a traffic incident occurs, such as an accident [132–134]. TD-VRP routes a fleet of vehicles by taking into account variable travel time on the road network [135, 136]. Shao et al. [74] observed BEVs in such time-dependent context in their E-VRP with charging time and variable travel time (EVRP-CTVTT). The recharging time is fixed to 30 minutes and batteries are always charged to full capacity. The authors discretized one day into two-minute intervals and applied the dynamic Dijkstra algorithm to find the shortest travel time path when weights in the graph are not constant. The authors presented a real-life problem and solved it by applying a genetic algorithm with a running time of three hours. Omidvar and Tavakkoli-Moghaddam [24] presented a model for routing AFVs that aims to avoid routing during congestion hours, when the pollution costs are high, by waiting at a customer’s locations. These costs were computed based on the road speed profile. Mirmohammadi et al. [62] presented MILP formulation for the periodic green heterogeneous VRP with time-dependent urban traffic and time windows, which is green in terms of the emission minimization and not the utilization of AFVs. The planning horizon is divided into several periods, during which static traffic conditions are considered.


## 4.12. Related Problems

### 4.12.1. Route Scheduling and Other Electric Variants

Many researchers are dealing with the scheduling of bus/taxi routes that are fixed or could be slightly altered. Li et al. [82] addressed the mixed bus fleet management problem (MBFM) composed of electric, diesel, compressed natural gas and hybrid-diesel buses. The authors maximized the total benefit of replacement of old vehicles with new ones under budget constraints while optimizing the route assignment for each bus during the planning period. Two routing procedures were developed to solve the recharging problem: single-period routing and routing across multiple periods of a day. Analysis of the results on the case study of Hong-Kong showed a cost-effective scheme of fleet configuration in the following order: mixed fleet, electric, natural gas, diesel, and then hybrid-diesel buses. Wen et al. [65] also dealt with the scheduling of electric buses (electric vehicle scheduling problem, E-VSP) in order to efficiently serve a set of timetabled bus trips. The linear partial and full recharge were considered, with the minimization of bus numbers and the total traveled distance. Lu et al. [83] addressed the problem of optimal scheduling of a taxi fleet with mixed BEVs and ICEVs to service advance reservations. The authors presented a multilayer taxi-flow time-space network in order to minimize the total operating fleet cost.Bruglieri et al. [77] formulated the MILP for a one-way electric car-sharing problem as the electric vehicle relocation problem (E-VReP). Relocation of BEVs is performed by workers who come with bicycles to the pickup points, put the bicycle into the BEV’s trunk, and drive the BEV to one of the delivery points. The authors balance a trade-off among the customers’ satisfaction, the workers’ workload balance, and the car-sharing provider’s objective.Masmoudi et al. [84] formulated the dial-a-ride problem with EVs and BSSs (DARP-EV) for customers with special needs and disabilities. The authors observed a fleet of heterogeneous BEVs with equal battery capacity but different available resources for the customers: handicapped person’s seat, stretcher, wheelchair, or accompanying person’s seat. The authors adapted the energy consumption model of Genikomsakis and Mitrentsis [85] with constant speed, acceleration, road slope, and varying mass while loading/unloading passengers.Paz et al. [86] addressed the multidepot electric vehicle location routing problem with time windows (MDEVLRPTW) with three variants: (i) with BSS (MDEVLRPTW-BS); (ii) with partial recharging (MDEVLRPTW-PR) at CSs that could also be located at customers’ locations; and (iii) with BSS and partial recharging (MDEVLRPTW-BSPR) where conventional charging is considered at customers’ locations and battery swap is performed at a separately located BSSs. Three MIP models were presented for each case. The authors determined the number and location of BSSs and/or CSs (without investment costs), the number and location of the depots, the number of vehicles, and the customers’ sequence in route. The proposed model allows minimizing the number of vehicles given a fixed number of CSs or BSSs, and vice versa. For the CS, linear recharge is considered, while the battery swap time is set to 10% of the corresponding linear recharge time. Results were presented by solving small instances with a commercial solver. The results showed that BSPR and PR variants provide better results than the BS variant and that BSPR tends to act as PR in most of the instances.Schiffer and Walther [89] extended the ELRPTWPR model of Schiffer and Walther [73] with uncertain customer patterns related to customers spatial distribution, demands, and service time windows. The authors formulated robust ELRPTWPR (RELRPTWPR) as a MILP program. Five scenarios were observed, which represent customer profile for five working days. The authors used set-based uncertainty representation, which is covered within scenarios that consider customer patterns in daily deliveries. To achieve a computationally tractable model of robust counterpart formulation, the adversarial approach was used, which uses a finite set of scenarios. The robust model was compared to the deterministic model (selected customer pattern), median model (average customer pattern), and worst-case model (worst-case customer pattern), where first the E-LRP was solved, and then the E-VRP. The deterministic models produced infeasible solutions in 35-63% of instances and incurred, on average, 5% higher costs than the robust approach. Also, the robust approach showed the most homogeneous configuration of CSs.Keskin et al. [80] extended E-VRPTW by considering waiting times at the CSs and proposed a solution method for solving small-sized instances. The authors assume single charger at CS, Poisson arrivals, exponential distribution of service times, first-in-first-out strategy, and penalties for late arrivals. The planning horizon is split into morning, afternoon, evening, and night period. The routing decisions are determined based on time-dependent waiting times at CSs.Kullman et al. [81] introduced E-VRP with public-private recharge strategy (E-VRP-PP), where demand at public CSs is unknown. The queuing at CS is modeled as M/M/c with first-in-first-out strategy, and c represents the number of identical chargers at CS. The problem is modeled as a Markov decision process, and an approximate dynamic programming solution is proposed. The authors extended the FRVCP formulation of Montoya et al. [72] with time-dependent waiting times and discrete charging decisions. Dynamic and static routing policies are proposed, where the static policies are incorporated as base policies for dynamic routing. Using dual bound on the optimal policy, the authors concluded that proposed routing policies are within 4.7%. The authors pointed out that the use of a public-private strategy outperformed the use of just private CSs, as overall savings were 20% higher.


### 4.12.2. Green and Pollution Routing Problem

The green vehicle routing problem focuses on the reduction of routing pollution on the environment. The key idea is to promote the use of sustainable energy sources and minimize overall emissions. Regarding the emission of BEV, Álvarez Fernández [137] indicated that emission savings of BEVs can be significant, up to 80% in some scenarios compared to ICEVs. The authors presented a model for prediction of a BEV’s GHG emissions linked to a route for each country, according to empirical measurements and prediction of routes’ energy consumption.Erdoğan and Miller-Hooks [23] first defined the GVRP with AFVs and full refuel, where also BEVs were included. Therefore, the E-VRP can be observed as a special case of the GVRP. The authors provided MILP formulation of the problem with replication of CSs vertices. Koç and Karaoglan [59] improved the MILP formulation of Erdoğan and Miller-Hooks [23] by removing dummy vertices that represent multiple visits to refueling stations. Instead, a binary variable related to traveling between two vertices by using the station in between is introduced (arc-duplicating). This formulation improved the average gap to the optimal integer solution on all scenarios compared to the model of Erdoğan and Miller-Hooks [23] by 25%. The authors proposed a combination of exact branch-and-cut procedure and simulated annealing with several moving strategies to solve the problem. Leggieri and Haouari [70] formulated the problem as a nonlinear MIP, which was linearized to derive MILP formulation by applying a reformulation-linearization technique. The authors used similar arc-duplicating to that of Koç and Karaoglan [59] but defined the binary variable as the sum of binary variables for the direct arc between the customers and binary variables for the same direct arc but with stations inserted. The proposed formulation without additional constraints was able to solve to optimality 97.5% of test instances, which substantially outperformed the formulations of Erdoğan and Miller-Hooks [23] and Koç and Karaoglan [59]. In addition, the same formulation with prepossessing consistently outperformed the exact branch-and-cut procedure of Koç and Karaoglan [59]. Andelmin and Bartolini [66] took a different approach for the exact solving of GVRP based on a set partitioning formulation in which columns correspond to feasible routes. The problem is modeled using a multigraph in which vertices correspond to customers and each arc represents a possible sequence of stations visited when traveling between customers. The authors added weak subset row inequalities, subset row inequalities, and k-path cuts, with the latter reported as the crucial part of the algorithm efficiency. The results showed tight lower bounds and high efficiency of the algorithm as instances with up to 110 customers were solved to optimality. Bruglieri et al. [94] presented a path-based approach for solving the GVRP, based on two phases. In the first phase, a set of feasible paths is generated, removing from it all those paths that are dominated. In the second phase, through MILP formulation, paths are combined to create routes for the final solution. The authors compared their results to previously observed research papers of Erdoğan and Miller-Hooks [23], Koç and Karaoglan [59], Leggieri and Haouari [70], and Andelmin and Bartolini [66]. On small-sized instances, the proposed algorithm produced optimal solutions with the lowest running time while on the larger instances it produced 2.3% worse solutions compared to Andelmin and Bartolini [66], but in much less time.Poonthalir and Nadarajan [38] formulated biobjective fuel efficient GVRP (F-GVRP) with the minimization of route costs and fuel consumption using goal programming. The authors analyzed the impact of varying speed on fuel consumption where speed is modeled using triangular distribution. To efficiently solve the problem, the authors applied a population-based metaheuristic. The tests were performed by constraining the route cost of the solution to the BKS of GVRP and then the fuel minimization was performed. The authors pointed out that the use of varying speed instead of constant speed reduced the fuel consumption on average by 20.48%, and in half of the instances it reduced the number of vehicles. The lower number of vehicles is a consequence of the lesser amount of time spent at the refueling stations as less fuel is consumed, so the vehicle can visit more customers. The authors pointed out that if the route cost constraint was relaxed, most of the instances might have better route costs than the BKSs of GVRP. The authors also analyzed the impact of speed interval size on the fuel consumption and concluded that fuel consumption and average speed increased with larger interval size.Normasari et al. [100] formulated capacitated GVRP (CGVRP) as a MILP program and solved the problem by applying a simulated annealing heuristic. Compared to the GVRP, the results indicated that the number of tours and total costs for the CGVRP increased by 38.8% and 32.5%, respectively.Macrina et al. [33] investigated the green mixed fleet VRP with partial recharges and time windows (GMFVRP-PRTW), in which polluting emissions are modeled through the functions of both traveled distance and vehicle load. The problem was formulated as a MILP and solved by a commercial solver for small instances and by an iterative local search metaheuristic for large instances. The overall pollution emissions are maintained below the designated threshold value, to simulate emission limits of freight vehicles. Macrina et al. [99] presented MILP formulation for a similar GMFVRP-PRTW problem but added more realistic consumption model presented in Section 3.2, where the authors considered different charging technologies (single technology per CS), cost of energy recharged at the depot, and the fuel cost. A large neighborhood search metaheuristic is applied to solve the problem.Koyuncu and Yavuz [98] compared the node- and arc-duplicating MILP formulation of mixed fleet GVRP (MGVRP) by taking into account the following: (i) a mixed fleet of AFVs and ICEVs; (ii) load and time window constraints; (iii) internal station visits (at customer) and external station visit (at separately located stations); (iv) fixed full, variable full, or variable partial refuel at the station; (v) site-dependent refueling rate; and (vi) the visit of two or more consecutive stations being prohibited. In the node-duplicating formulations [23, 25], each refueling vertex is replicated at most N times, where N is the number of customers. In the arc-duplicating formulation, two subsets of arcs are contained: arcs between customers (direct arcs) and refueling arcs which contain inserted station between customers [70]. In the preprocessing step for both formulations, unfeasible arcs are removed and tighter lower and upper bounds for auxiliary variables and number of vehicles are presented. On the GVRP test instances of Erdoğan and Miller-Hooks [23], the authors reported that the node- and arc-duplicating formulations yielded an average optimality gap of 20.84% and 9.07%, respectively. On the GVRP instances of Yavuz [138], significant conclusions could not be drawn.As BEVs have no local emissions, the E-VRP is closely related to the minimization of GHG emissions, where a problem-specific GVRP variant called the pollution routing problem (PRP) was introduced [32, 40]. Bektaş and Laporte [40] formulated the PRP as a nonlinear MIP where the main objective is the optimization of vehicle speed and GHG emissions. Demir et al. [32] concluded that speed optimization improves the PRP solution and minimizes fuel consumption and driver costs.It can be noted that several papers are dealing with fuel consumption and pollution emissions and not AFVs, but the termgreen VRP is used to indicate a more ecologically aware routing problem [62, 139].


### 4.12.3. More General Formulation of the E-VRP

The E-VRP can be formulated more generally as the VRP with intermediate stops (VRPIS), in which the vehicle visits intermediate or intraroute facilities to replenish/unload the goods or to refuel [55]. The LRPIF extends the problem in terms of the facilities location optimization [88]. The service time depends on the load and fuel level at the arrival at the facility. Schneider et al. [55] presented the E-VRP with recharging facilities (EVRPRF) as a special case of VRPIS, where BEVs have a fixed recharging time and limited vehicle load capacity without time windows constraints.


### 4.12.4. Multiobjective Variants

Most of the researchers are using a single-objective function that represents travel distance, total costs, energy consumption, total time, etc. The multiobjective variants are still relatively scarcely applied. Generally, the solution of the multiobjective problem is not the optimal solution for all of the objectives, but rather it is satisfactory in those terms.Demir et al. [37] introduced the biobjective PRP with the minimization of fuel consumption and total driving time. The fuel consumption is based on a model similar to the one presented by (6) but with constant speed. An adaptive large neighborhood search was used to generate nondominated/Pareto optimal solutions while an additional hybrid procedure was used to obtain the solution to the biobjective model. The experiments showed that there is no need to prolong driving time to achieve a significant reduction in fuel consumption.The already mentioned F-GVRP of Poonthalir and Nadarajan [38] (Section 4.12.2) used biobjective minimization of route costs and fuel consumption. Wang et al. [91] considered multiobjective function to determine optimal BEV routes between origin-destination pairs. Three objectives are minimized: travel time, energy consumption, and charging cost. Travel time includes driving time, queuing time, and charging time, while charging cost includes electricity cost, service cost, and parking fee. The multiobjective model is transformed into a single objective by applying a fuzzy programming approach for objectives’ membership functions and fuzzy preference relations to obtain weighting coefficients of each objective. A genetic algorithm is designed to obtain an optimal solution. The influences of various driver trade-offs (weighting conditions) among different objectives are explored using Pareto curves. It was observed that driver trade-offs have an effect on the travel times, and almost no effect on the energy consumption and charging costs.Amiri et al. [39] optimized the location of the battery swapping in the BSSs network and the charging schedule of depleted batteries. The multiobjective model minimizes three objectives: battery charging and power loss costs, deviation from the nominal voltages, and network capacity releasing. The nondominated sorting genetic algorithm is applied to solve the problem. The results on the test case showed that an optimal charging schedule in terms of cost and network constraints can be achieved. Besides, the electricity load profile is more leveled.Bruglieri et al. [77] focused on the E-VReP with three objectives: minimization of workers used to relocate the BEVs (service provider objective), maximization of relocation requests served (users’ objective), and minimization of the longest routes’ duration (workload balance). The problem is solved by approximating the Pareto optimal front based on a two-phase approach: first, feasible solutions are generated by randomized search heuristics, and in the second phase, to find a set of nondominated solutions, the MILP formulation is solved by epsilon constraint programming.


## 4.12.1. Route Scheduling and Other Electric Variants

Many researchers are dealing with the scheduling of bus/taxi routes that are fixed or could be slightly altered. Li et al. [82] addressed the mixed bus fleet management problem (MBFM) composed of electric, diesel, compressed natural gas and hybrid-diesel buses. The authors maximized the total benefit of replacement of old vehicles with new ones under budget constraints while optimizing the route assignment for each bus during the planning period. Two routing procedures were developed to solve the recharging problem: single-period routing and routing across multiple periods of a day. Analysis of the results on the case study of Hong-Kong showed a cost-effective scheme of fleet configuration in the following order: mixed fleet, electric, natural gas, diesel, and then hybrid-diesel buses. Wen et al. [65] also dealt with the scheduling of electric buses (electric vehicle scheduling problem, E-VSP) in order to efficiently serve a set of timetabled bus trips. The linear partial and full recharge were considered, with the minimization of bus numbers and the total traveled distance. Lu et al. [83] addressed the problem of optimal scheduling of a taxi fleet with mixed BEVs and ICEVs to service advance reservations. The authors presented a multilayer taxi-flow time-space network in order to minimize the total operating fleet cost.Bruglieri et al. [77] formulated the MILP for a one-way electric car-sharing problem as the electric vehicle relocation problem (E-VReP). Relocation of BEVs is performed by workers who come with bicycles to the pickup points, put the bicycle into the BEV’s trunk, and drive the BEV to one of the delivery points. The authors balance a trade-off among the customers’ satisfaction, the workers’ workload balance, and the car-sharing provider’s objective.Masmoudi et al. [84] formulated the dial-a-ride problem with EVs and BSSs (DARP-EV) for customers with special needs and disabilities. The authors observed a fleet of heterogeneous BEVs with equal battery capacity but different available resources for the customers: handicapped person’s seat, stretcher, wheelchair, or accompanying person’s seat. The authors adapted the energy consumption model of Genikomsakis and Mitrentsis [85] with constant speed, acceleration, road slope, and varying mass while loading/unloading passengers.Paz et al. [86] addressed the multidepot electric vehicle location routing problem with time windows (MDEVLRPTW) with three variants: (i) with BSS (MDEVLRPTW-BS); (ii) with partial recharging (MDEVLRPTW-PR) at CSs that could also be located at customers’ locations; and (iii) with BSS and partial recharging (MDEVLRPTW-BSPR) where conventional charging is considered at customers’ locations and battery swap is performed at a separately located BSSs. Three MIP models were presented for each case. The authors determined the number and location of BSSs and/or CSs (without investment costs), the number and location of the depots, the number of vehicles, and the customers’ sequence in route. The proposed model allows minimizing the number of vehicles given a fixed number of CSs or BSSs, and vice versa. For the CS, linear recharge is considered, while the battery swap time is set to 10% of the corresponding linear recharge time. Results were presented by solving small instances with a commercial solver. The results showed that BSPR and PR variants provide better results than the BS variant and that BSPR tends to act as PR in most of the instances.Schiffer and Walther [89] extended the ELRPTWPR model of Schiffer and Walther [73] with uncertain customer patterns related to customers spatial distribution, demands, and service time windows. The authors formulated robust ELRPTWPR (RELRPTWPR) as a MILP program. Five scenarios were observed, which represent customer profile for five working days. The authors used set-based uncertainty representation, which is covered within scenarios that consider customer patterns in daily deliveries. To achieve a computationally tractable model of robust counterpart formulation, the adversarial approach was used, which uses a finite set of scenarios. The robust model was compared to the deterministic model (selected customer pattern), median model (average customer pattern), and worst-case model (worst-case customer pattern), where first the E-LRP was solved, and then the E-VRP. The deterministic models produced infeasible solutions in 35-63% of instances and incurred, on average, 5% higher costs than the robust approach. Also, the robust approach showed the most homogeneous configuration of CSs.Keskin et al. [80] extended E-VRPTW by considering waiting times at the CSs and proposed a solution method for solving small-sized instances. The authors assume single charger at CS, Poisson arrivals, exponential distribution of service times, first-in-first-out strategy, and penalties for late arrivals. The planning horizon is split into morning, afternoon, evening, and night period. The routing decisions are determined based on time-dependent waiting times at CSs.Kullman et al. [81] introduced E-VRP with public-private recharge strategy (E-VRP-PP), where demand at public CSs is unknown. The queuing at CS is modeled as M/M/c with first-in-first-out strategy, and c represents the number of identical chargers at CS. The problem is modeled as a Markov decision process, and an approximate dynamic programming solution is proposed. The authors extended the FRVCP formulation of Montoya et al. [72] with time-dependent waiting times and discrete charging decisions. Dynamic and static routing policies are proposed, where the static policies are incorporated as base policies for dynamic routing. Using dual bound on the optimal policy, the authors concluded that proposed routing policies are within 4.7%. The authors pointed out that the use of a public-private strategy outperformed the use of just private CSs, as overall savings were 20% higher.


## 4.12.2. Green and Pollution Routing Problem

The green vehicle routing problem focuses on the reduction of routing pollution on the environment. The key idea is to promote the use of sustainable energy sources and minimize overall emissions. Regarding the emission of BEV, Álvarez Fernández [137] indicated that emission savings of BEVs can be significant, up to 80% in some scenarios compared to ICEVs. The authors presented a model for prediction of a BEV’s GHG emissions linked to a route for each country, according to empirical measurements and prediction of routes’ energy consumption.Erdoğan and Miller-Hooks [23] first defined the GVRP with AFVs and full refuel, where also BEVs were included. Therefore, the E-VRP can be observed as a special case of the GVRP. The authors provided MILP formulation of the problem with replication of CSs vertices. Koç and Karaoglan [59] improved the MILP formulation of Erdoğan and Miller-Hooks [23] by removing dummy vertices that represent multiple visits to refueling stations. Instead, a binary variable related to traveling between two vertices by using the station in between is introduced (arc-duplicating). This formulation improved the average gap to the optimal integer solution on all scenarios compared to the model of Erdoğan and Miller-Hooks [23] by 25%. The authors proposed a combination of exact branch-and-cut procedure and simulated annealing with several moving strategies to solve the problem. Leggieri and Haouari [70] formulated the problem as a nonlinear MIP, which was linearized to derive MILP formulation by applying a reformulation-linearization technique. The authors used similar arc-duplicating to that of Koç and Karaoglan [59] but defined the binary variable as the sum of binary variables for the direct arc between the customers and binary variables for the same direct arc but with stations inserted. The proposed formulation without additional constraints was able to solve to optimality 97.5% of test instances, which substantially outperformed the formulations of Erdoğan and Miller-Hooks [23] and Koç and Karaoglan [59]. In addition, the same formulation with prepossessing consistently outperformed the exact branch-and-cut procedure of Koç and Karaoglan [59]. Andelmin and Bartolini [66] took a different approach for the exact solving of GVRP based on a set partitioning formulation in which columns correspond to feasible routes. The problem is modeled using a multigraph in which vertices correspond to customers and each arc represents a possible sequence of stations visited when traveling between customers. The authors added weak subset row inequalities, subset row inequalities, and k-path cuts, with the latter reported as the crucial part of the algorithm efficiency. The results showed tight lower bounds and high efficiency of the algorithm as instances with up to 110 customers were solved to optimality. Bruglieri et al. [94] presented a path-based approach for solving the GVRP, based on two phases. In the first phase, a set of feasible paths is generated, removing from it all those paths that are dominated. In the second phase, through MILP formulation, paths are combined to create routes for the final solution. The authors compared their results to previously observed research papers of Erdoğan and Miller-Hooks [23], Koç and Karaoglan [59], Leggieri and Haouari [70], and Andelmin and Bartolini [66]. On small-sized instances, the proposed algorithm produced optimal solutions with the lowest running time while on the larger instances it produced 2.3% worse solutions compared to Andelmin and Bartolini [66], but in much less time.Poonthalir and Nadarajan [38] formulated biobjective fuel efficient GVRP (F-GVRP) with the minimization of route costs and fuel consumption using goal programming. The authors analyzed the impact of varying speed on fuel consumption where speed is modeled using triangular distribution. To efficiently solve the problem, the authors applied a population-based metaheuristic. The tests were performed by constraining the route cost of the solution to the BKS of GVRP and then the fuel minimization was performed. The authors pointed out that the use of varying speed instead of constant speed reduced the fuel consumption on average by 20.48%, and in half of the instances it reduced the number of vehicles. The lower number of vehicles is a consequence of the lesser amount of time spent at the refueling stations as less fuel is consumed, so the vehicle can visit more customers. The authors pointed out that if the route cost constraint was relaxed, most of the instances might have better route costs than the BKSs of GVRP. The authors also analyzed the impact of speed interval size on the fuel consumption and concluded that fuel consumption and average speed increased with larger interval size.Normasari et al. [100] formulated capacitated GVRP (CGVRP) as a MILP program and solved the problem by applying a simulated annealing heuristic. Compared to the GVRP, the results indicated that the number of tours and total costs for the CGVRP increased by 38.8% and 32.5%, respectively.Macrina et al. [33] investigated the green mixed fleet VRP with partial recharges and time windows (GMFVRP-PRTW), in which polluting emissions are modeled through the functions of both traveled distance and vehicle load. The problem was formulated as a MILP and solved by a commercial solver for small instances and by an iterative local search metaheuristic for large instances. The overall pollution emissions are maintained below the designated threshold value, to simulate emission limits of freight vehicles. Macrina et al. [99] presented MILP formulation for a similar GMFVRP-PRTW problem but added more realistic consumption model presented in Section 3.2, where the authors considered different charging technologies (single technology per CS), cost of energy recharged at the depot, and the fuel cost. A large neighborhood search metaheuristic is applied to solve the problem.Koyuncu and Yavuz [98] compared the node- and arc-duplicating MILP formulation of mixed fleet GVRP (MGVRP) by taking into account the following: (i) a mixed fleet of AFVs and ICEVs; (ii) load and time window constraints; (iii) internal station visits (at customer) and external station visit (at separately located stations); (iv) fixed full, variable full, or variable partial refuel at the station; (v) site-dependent refueling rate; and (vi) the visit of two or more consecutive stations being prohibited. In the node-duplicating formulations [23, 25], each refueling vertex is replicated at most N times, where N is the number of customers. In the arc-duplicating formulation, two subsets of arcs are contained: arcs between customers (direct arcs) and refueling arcs which contain inserted station between customers [70]. In the preprocessing step for both formulations, unfeasible arcs are removed and tighter lower and upper bounds for auxiliary variables and number of vehicles are presented. On the GVRP test instances of Erdoğan and Miller-Hooks [23], the authors reported that the node- and arc-duplicating formulations yielded an average optimality gap of 20.84% and 9.07%, respectively. On the GVRP instances of Yavuz [138], significant conclusions could not be drawn.As BEVs have no local emissions, the E-VRP is closely related to the minimization of GHG emissions, where a problem-specific GVRP variant called the pollution routing problem (PRP) was introduced [32, 40]. Bektaş and Laporte [40] formulated the PRP as a nonlinear MIP where the main objective is the optimization of vehicle speed and GHG emissions. Demir et al. [32] concluded that speed optimization improves the PRP solution and minimizes fuel consumption and driver costs.It can be noted that several papers are dealing with fuel consumption and pollution emissions and not AFVs, but the termgreen VRP is used to indicate a more ecologically aware routing problem [62, 139].


## 4.12.3. More General Formulation of the E-VRP

The E-VRP can be formulated more generally as the VRP with intermediate stops (VRPIS), in which the vehicle visits intermediate or intraroute facilities to replenish/unload the goods or to refuel [55]. The LRPIF extends the problem in terms of the facilities location optimization [88]. The service time depends on the load and fuel level at the arrival at the facility. Schneider et al. [55] presented the E-VRP with recharging facilities (EVRPRF) as a special case of VRPIS, where BEVs have a fixed recharging time and limited vehicle load capacity without time windows constraints.


## 4.12.4. Multiobjective Variants

Most of the researchers are using a single-objective function that represents travel distance, total costs, energy consumption, total time, etc. The multiobjective variants are still relatively scarcely applied. Generally, the solution of the multiobjective problem is not the optimal solution for all of the objectives, but rather it is satisfactory in those terms.Demir et al. [37] introduced the biobjective PRP with the minimization of fuel consumption and total driving time. The fuel consumption is based on a model similar to the one presented by (6) but with constant speed. An adaptive large neighborhood search was used to generate nondominated/Pareto optimal solutions while an additional hybrid procedure was used to obtain the solution to the biobjective model. The experiments showed that there is no need to prolong driving time to achieve a significant reduction in fuel consumption.The already mentioned F-GVRP of Poonthalir and Nadarajan [38] (Section 4.12.2) used biobjective minimization of route costs and fuel consumption. Wang et al. [91] considered multiobjective function to determine optimal BEV routes between origin-destination pairs. Three objectives are minimized: travel time, energy consumption, and charging cost. Travel time includes driving time, queuing time, and charging time, while charging cost includes electricity cost, service cost, and parking fee. The multiobjective model is transformed into a single objective by applying a fuzzy programming approach for objectives’ membership functions and fuzzy preference relations to obtain weighting coefficients of each objective. A genetic algorithm is designed to obtain an optimal solution. The influences of various driver trade-offs (weighting conditions) among different objectives are explored using Pareto curves. It was observed that driver trade-offs have an effect on the travel times, and almost no effect on the energy consumption and charging costs.Amiri et al. [39] optimized the location of the battery swapping in the BSSs network and the charging schedule of depleted batteries. The multiobjective model minimizes three objectives: battery charging and power loss costs, deviation from the nominal voltages, and network capacity releasing. The nondominated sorting genetic algorithm is applied to solve the problem. The results on the test case showed that an optimal charging schedule in terms of cost and network constraints can be achieved. Besides, the electricity load profile is more leveled.Bruglieri et al. [77] focused on the E-VReP with three objectives: minimization of workers used to relocate the BEVs (service provider objective), maximization of relocation requests served (users’ objective), and minimization of the longest routes’ duration (workload balance). The problem is solved by approximating the Pareto optimal front based on a two-phase approach: first, feasible solutions are generated by randomized search heuristics, and in the second phase, to find a set of nondominated solutions, the MILP formulation is solved by epsilon constraint programming.


## 4.13. Summary of the E-VRP Variants

Table1 presents a summary of the E-VRP variants and related problems in the available literature. In total, 79 papers were presented in the table. The problem characteristics are observed in the following order: vehicle load (cargo) capacity (48), linear charging/refueling function at station (34), customer time windows (34), partial recharging strategy (30), fixed refueling/recharging time (19), different charging technologies at CSs (17), heterogeneous (mixed) fleet (16), energy consumption model (15), BSS (9), nonlinear charging function (8), dynamic traffic conditions (8), location routing (7), capacitated stations (6), and HEVs (5). We can point out that only a few papers are dealing with (i) simultaneous CS siting and BEV routing (E-LRP) where great savings can be achieved in companies whose business plans include CS investments, as E-LRP generally gives better routing configuration than the E-VRP [73, 88]; (ii) nonlinear charging function, negligence of which can lead to infeasible solutions and additional penalty costs [68, 72]; (iii) HEVs (PHEVs) that have an advantage in several routing scenarios but are being researched only recently due to the higher complexity of the problem as HEVs can change their propulsion mode at any time during the route [96]; (iv) dynamic traffic conditions, which significantly influence BEVs’ energy consumption [74]; (v) different charging technologies at CSs, which could reduce the overall costs by selecting the best possible charger option in each occasion [46, 79]; (vi) CS related challenges: capacitated CS [68, 75, 80], CS reservations [44], and public-private recharge strategy [81]; and (vii) robust E-VRP variant resistant to uncertainties in travel time, service time, and demand [89].


## 5. Problem Solving Methods

Since VRP is a well-researched problem, a large number of procedures for solving the problem have been proposed. Due to the NP-hardness of the problem and a large number of customers in real-life problems, most of the procedures used in real-life applications are heuristics, metaheuristics, and hybrid combinations. For small-sized problems, a great number of exact procedures have been proposed. Many of the VRP procedures in the available literature are with an adaptation applicable for solving the E-VRP. An overview of the applied procedures for solving the E-VRP and related problems is presented in Table2.Table 2
Overview of the applied procedures.
Reference Problem name NF Initial Hybrid Metaheuristics Heuristics A Exact & Software Instances Description Bektaş and Laporte [40] PRP     CPLEX - small Generated  Conrad and Figliozzi [22] Recharging VRP  Route construction Iterative route construction and improvement   Solomon - CVRP-TW ( 1 ) Sequential construction of routes(2) Route improvement - improves a group or set of routes  Gonçalves et al. [21] VRPPD with mixed fleet     CPLEX - small  Demir et al. [32] PRP  CWS ALNS + speed opt.   PRP ALNS destroy: Random, WorstDistance, WorstTime, Route, Shaw, ProximityBased, TimeBased, DemandBased, HistoricalKnowledge, Neighborhood, Zone, NodeNeighborhoodALNS repair: Greedy, Regret, GreedyWithNoise, RegretWithNoise, ZoneSpeed opt.: first optimal speeds are computed to find violations if any, and then in the second part the current speeds are revised to optimize the fuel consumption  Erdoğan and Miller-Hooks [23] GVRP  ( 1 ) MCWS(2) DBCA LS  CPLEX - small GVRP LS operator: Exchange  Omidvar and Tavakkoli-Moghaddam [24] GVRP   ( 1 ) SA(2) GA SA CPLEX - small Generated - Solomon procedure No description of the implementation  Abdallah [41] PHEVRPTW  Push forward insertion heuristic + LP optimization of charges Heuristics for feasible upper bound, TS SA MATLAB, CPLEX - small, DP Adapted Solomon Compare three procedures for small instances:(1) Master LP: three different relaxations by applying cutting plane algorithm + CPLEX for MIP subproblem;(2) Best Lagrangian bound + heuristic(3) CPLEX TS neighborhood operators: λ - interchange (λ=1,2), DP for determining optimal route cost  Barco et al. [34, 42] E-VRP and charge scheduling  Random Differential evolution, GA   Real (1) Determine the set of optimal routes(2) Assignment of routes over the planning horizon(3) Charge scheduling + rerouting to CSsMutation: OR, XOR between 3 individuals Crossover: random 2-point  Van Duin et al. [43] FSMVRPTW, EVFSMVRPTW  Sequential insertion heuristics LS    Shortrec softwareLS operators: Bräysy and Gendreau [26]  Adler and Mirchandani [44] Online routing of BEVs     Approximate DP Real Markov chance-decision process  Alesiani and Maslekar [45] BEV routing X Random GA   Generated Solution: a two-dimensional binary array of routes and CSsAfter assigning CS to routes the TSP route is determinedInitial: no charging, battery constrained SPPCrossover: 2-point, the second parent is chosen to reduce the concurrency at CSsSelection: random or based on the fitnessMutation: random - remove/insert CS, ExchangeLS: removing/inserting vertices (Relocate)  Demir et al. [37] PRP   ALNS + speed opt. + heuristic   Generated ALNS of Demir et al. [32] + speed optimization - to generate non-dominated Pareto optimal solutionsHeuristic: adaptive weighting with ϵ-constraint method for bi-objective solutions  Felipe et al. [46] GVRP-MTPR  Constructive-k -PseudoGreedy LS(48A) SA  GVRP-MTPR, GVRP, E-VRPTW LS (48A): RechargeRelocation, 2-Opt, reinsertion with four options: first or best and yes/no update of savings  Preis et al. [35] Energy-optimized routing of BEVs X CWS TS  MIP solver - small Generated Neighborhood operator: Relocate  Sassi et al. [47] HEVRP-TDMF  CRH ILS+LNS   Real CRH -(1) Feasible charging scheme at the depot, (2) Insertion heuristicsLNS destroy: RandomLNS repair: Random, Regret, Score-based  Schneider et al. [25] E-VRPTW X Sweep VNS+TS SA CPLEX - small E-VRPTW, MDVRPI, GVRP VNS - 15 neighborhood structures based on the cyclic exchange moveTS - operators: 2-opt∗, Relocate, Exchange, StationInReDiversification mechanism based on the attributes frequency  Zündorf [50] EVRC   Heuristics to speed up exact algorithms  Charging function propagating algorithm, CH,A∗ Real Battery CSPPMulti-objective search using Pareto-sets  Bruglieri et al. [29, 51] E-VRPTW   VNS with branching  CPLEX E-VRPTW - small VNS with Hamming distance of three: removing CS, inserting CS, merging  Goeke and Schneider [30] E-VRPTWMF X Modified insertion heuristics ALNS+LS SA  E-VRPTWMF, VRPTW, E-VRPTW Initial - long BEV routes are converted to ICEV routesALNS destroy: Random, Worst, Shaw, Cluster, StationVicinityALNS repair: Greedy, Regret, GRASPLS operators: Relocate, Exchange, 2-opt∗, StationInRe  Lebeau et al. [52] FSMVRPTW-EV   CWS + tree branching   FSMVRPTW - EV The tree is used for different combinations of vehicles - in each branch, one vehicle is selected and vertices are inserted  Moghaddam [53] E-VRPTWPR     Gurobi Generated  Pourazarm et al. [54] Single/Multi BEV routing     DP Generated Single - electric vehicle CSPPGrouping subsets of vehicles into sub-flows  Sassi et al. [48] VRP-MFHEV  CRH ITS+LNS RR  Generalized E-VRPTW CRH - without charge scheduling at the depotLNS destroy: RandomLNS repair: RandomFirst and best improvement  Sassi et al. [49] VRP-HFCC  CRH ILS+LNS RR  Real, generalized E-VRPTW LNS destroy: RandomLNS repair: Random, RegretPerturbation: LNS random destroy and repair  Schneider et al. [55] VRPIS (EVRPRF) X MCWS AVNS+LS SA CPLEX - small EVRPRF, GVRP, VRPIRF AVNS shaking phase -24 neighborhood structures: 6 sequence relocation, 18 cyclic exchange(1) First route selection operators: Random, Length, LengthPerDemandUnit, FacilityDensity, FacilityDetour, (2) Next route is selected from a set of spatially close routes, (3) Vertex selection: Random, DistanceToNextRoute, DistanceToNeighboringVertices LS: FacilityInsertion, 2-opt, Or-Opt, Relocate, FacilityReplacement, Exchange  Yang and Sun [56] BSS-EV-LRP  Radius covering + MCWS TS-MCWS+LS  CPLEX - small BSS-EV-LRP LS: Relocate, ExchangeNeighborhood search: replace located and unlocated BSS X Modified sweep + greedy SIGALNS SA Remove all BSSs and then apply the iterated greedy procedure to insert best BSSsALNS destroy: Random, Worst, Related, RequestGraph, SinglePoint, TwoPoint, Binary, StationBasedALNS repair: Greedy, Regret-k, AdvanceGreedy, AdvanceRegret-kImprovement: try to remove BSSs  Desaulniers et al. [57] E-VRPTW -SF/MF/SP/MP     Branch-price-and-cut GENCOL + CPLEX E-VRPTW with modif. for SF/SP  Doppstadt et al. [58] HEV-TSP  ( 1) Random(2) NNH(3) Concorde - TSP ITS  CPLEX HEV-TSP Initial route is driven in combustion mode + LS - 2-Opt hill-climbing 4:1 ratio - four arcs have to be driven in charging mode in order to drive the last one in electric modeLS: λ-mode change (λ=2,5) - changes only the mode, λ-Opt with mode change (λ=2,3) - changes the route configuration and the mode ITS: (1) Search: 2-opt with mode change + arcs in tabu list, (2) Hill-climbing λ- mode change (first or best), (3) Diversification: 3-Opt with the mode change  Hiermann et al. [31] E-FSMFTW X Iterative route construction ALNS+LS SA Branch-and-price+DP E-FSMFTW, E-VRPTW Labeling algorithm for CSs and new extension functionsALNS destroy: Random, RandomAndRelated, WorstAndRelated, InefficientRouteAndNeighbour, TargetALNS repair: SequentialNode, ParallelRegret, SemiParallelConstruction, SemiParallelInsertionLS (best of first 50): 2-opt, 2-opt∗, Relocate, Exchange, InsertRemoveIF (CS), Resize (type), RelocateAndResize (type)  Keskin and Çatay [28] E-VRPTWPR  Insertion heuristics ALNS SA CPLEX - small E-VRPTW ALNS destroy options: remove customer only or with preceding/succeeding CSCustomer: Random, WorstDistance, WorstTime, Shaw, ProximityBased, DemandBased, TimeBased, ZoneRemovalRoute: Random, GreedyStation: Random, WorstDistance, WorstChargeUsage, FullChargeALNS repair: Customer: Greedy, Regret-2/3, TimeBased, Zone;Station: Greedy, GreedyWithComparison, Best  Koç and Karaoglan [59] GVRP  MCWS SA SA Branch-and- cut GVRP Improved MIP formulation of GVRPSA improvement operators: Merge routes, Swap, station Add/DropBranch-and-cut is applied to the initial solutionOn optimal fractional solution, the heuristic is applied to obtain a feasible solution which is further improved by SA  Lin et al. [60] E-VRP     MATLAB Real The linearized model formulation including energy costs and battery charging time  Masliakova [61] Routing and charging of electric buses  Random CSs positions + route construction GA+AC   Real GA: random 2-point crossover, mutation: move the CS to the nearest vertexAC - the ranking of individuals in GA  Mirmohammadi et al. [62] Periodic green VRP     CPLEX Generated  Montoya et al. [63] GVRP  ((1) random NNH, (2) random nearest insertion, (3) random best insertion) + split Modified multi-space sampling heuristic + set partitioning  Gurobi - set partitioning GVRP ( 1 ) Build the TSP route constructed by applying one of the construction procedures(2) Split TSP route and repair (CSPP) - pulse algorithm(3) Set partitioning over a set of routes  Roberti and Wen [64] E-TSPTW X Random Three-phase heuristic: VND + LS + DP  CPLEX, DP E-TSPTW, TSPTW, E-VRPTW - single vehicle (1) Initial - TSP(2) Perturbation: random customer relocation - TSP(3) TW feasible - TSPTW + (VND+LS: Relocate backward and forward, 2-Opt, Exchange)(4) Energy feasible - DP - labeling algorithm  Schiffer et al. [11, 27] E-LRPTWPR X CWS ALNS+DP+LS  DP Real ALNS destroy: large (Add, Drop, SwapPerfect, SwapPerfectOut) - change CS configuration; small (Worst, Related, Route, Shaw, StationVicinity) - change route configurationALNS repair: SequentialNode, SequentialPerturbedLS (best of 100): 2-Opt∗, Or-opt, Exchange, Relocate, StationInsertionDP - optimal determination of recharging visits and charging amount  Wen et al. [65] E-VSP  Greedy constructive ALNS SA CPLEX - small and post opt. E-VSP ALNS destroy: Random, TimeRelated, NeighboringSchedule (problem specific)ALNS repair: regret-2,3,4 (deterministic and stochastic version)Post-optimization: a set of feasible solutions are used in set partitioning model  Andelmin and Bartolini [66] GVRP     Exact procedure GVRP, generated Multigraph - arc represents possible paths between customers with inserted CS Set partitioning formulation + subset row inequalities andk-path cuts (Chvatal–Gomory (CG) cuts of rank one) - MILP formulation of separation + cut-and-column algorithm  Çatay and Keskin [67] E-VRPTWPR     CPLEX - small E-VRPTW  Froger et al. [68] E-VRP-NL-C  NNH-TSP Two-stage heuristic  Gurobi - MILP models Modified E-VRP-NL First stage without CS capacity constraint: ILS - giant TSP route + perturbation (random double bridge) + split procedure on acyclic graph with heuristic for CSs + LS: 2-Opt, Relocate and MILP model for charging decisions; local optimum solutions are stored in a pool of solutionsSecond stage: solutions are assembled using Benders’ like decomposition into route selection master problem, and CS capacity management sub-problemThree strategies for CS capacity management: (1) no revision (only feasibility check), (2) delay the charging operations to satisfy CS capacity constraints, and (3) revise the charging amountThe set partitioning model is solved using the branch-and-bound algorithm and cuts for CS capacity management  Hof et al. [69] BSS-EV-LRP X MCWS AVNS+LS SA CPLEX - small BSS-EV-LRP - Yang and Sun [56] and new AVNS - Schneider et al. [55] - 32 neighborhood structures - (24) + (8) new one based on the facility-removal (4) and facility-replacement (4) LS: 2-Opt, Or-Opt, Relocate, Exchange, FacilityInsertion, FacilityReplacement, FacilityExchange  Leggieri and Haouari [70] GVRP     CPLEX GVRP  Mancini [71] HVRP  Run MILP for5-60 s LNS  Xpress HVRP, GVRP Two routes are destroyed and the MILP mathematical model is run again for10 s only on the customers that were previously involved in the selected routes  Montoya et al. [72] E-VRP-NL  NNH-TSP + split ILS (VND)+ heuristic concentration  FRVCP-Gurobi E-VRP-NL Sequencing first, charging (FRVCP) secondVND: Relocate, 2-Opt and GlobalChargingImprovementPerturbation: build TSP route, perturb and splitHeuristic concentration: set partitioning over a set of routesFRVCP- Gurobi and a greedy heuristic  Schiffer and Walther [73] E-LRPTWPR     Gurobi E-VRPTW - small  Shao et al. [74] EVRP-CTVTT  Encoded mode GA   Real Dynamic Dijkstra algorithm for time-dependent SPP, single string encoding, each individual after the crossover or mutation needs to satisfy all the constraints - if not run the procedure again, elitism  Sweda et al. [75] Adaptive routing and recharging policies of EVs   Two specially designed heuristics   Generated Origin-destination,(1) Optimal a priori policy - CSPP + searching for policies with a lower costs(2.1) Two-stage heuristic: select the EV path to follow and select adaptive recharging policy for the fixed path(2.2) Adjusting vehicle trajectory - find optimal policy between a priori stop locations and permit adaptive decision making for both routing and recharging  Vincent et al. [76] HVRP X Modified NNH SA SA  Various CVRP Neighborhood operators: Insert, Swap, ReverseRestart strategy  Amiri et al. [39] BSS location & scheduling   GA   Generated Non-dominated sorting GA  Bruglieri et al. [77] E-VReP  Construction heuristics Two-phase heuristic  CPLEX Real ( 1 ) First phase: generate a set of feasible solutions: three random constructive heuristics(2) Second phase: ϵ-constraint programming  Joo and Lim [78] EV routing   AC   Generated Energy SPPAC - each ant selects next vertex by a stochastic greedy search based on the energy consumption and speed  Keskin and Çatay [79] E-VRPTW-FC  ( 1 ) BKS(2) E-VRPTW - PR sol. FC(3) Greedy ALNS + opt. SA CPLEX - small and opt. GVRP-MTPR, E-VRPTW ALNS of Keskin and Çatay [28] with modf.: Station destroy procedures: Random, WorstDistance, LeastUsedStation, ExpensiveStation; Customer insertion procedures - only fastest recharging option; modified station insertions algorithms which start from the cheapest and go to the more expensive charging optionsPost optimization of charging decisions (location, amount and type)  Keskin et al. [80] E-VRPTW with waiting times   ALNS + opt. SA CPLEX - small and opt. E-VRPTW ALNS of Keskin and Çatay [28]Post optimization of charging decisions (location, amount and type)  Kullman et al. [81] E-VRP-PP   Routing policies  Exact, FRVCP Generated Markov decision process - approximated DP solutionThe Benders-based branch-and-cut algorithm to solve the decompositionIntegration of static policies into dynamic lookahead policies: pre-decision, post-decision, and one-step rolloutsExact labeling for FRVCP with time-dependent waiting times and discrete charging decisionsNonlinear information penalties that tighten the perfect information bound of the optimal policy  Li et al. [82] MBFM & recharging problem     CPLEX Real Case study, new life additional benefit costs method for the replacement of vehicles (CPLEX); two routing methods: single/multi-period approach for solving the recharging problem  Lu et al. [83] MTFSP   Network decomposition by vehicles and time  CPLEX - sub-models Real Multi-layer taxi-flow time-space networkNetwork decomposition by vehicles: first MTFSP sub-model for k EVs is solved and then the MTFSP for the rest of ICEVsNetwork composition by time: in each slice of the time-sliced network to each request the nearest vehicles are assigned and then the recharging of EV is performed  Masmoudi et al. [84] DARP-EV X Modified insertion heuristics EVO-VNS CA, RR CPLEX - small HDARP, DARP-EV ( 1 ) EVO - chromosome encoded as a sequence of vehicles and pickup/drop off vertices:(1a) population is divided in m groups(1b)s solutions are selected from each group (roulette wheel), and if the best solution was not improved for the last three iterations, merge crossover (MX1) is applied between one solution from s and one from the rest of the solutions in the groups that were not selected(2) VNS - neighborhood operators: Swap, RemoveSequence, CrossExchange; three versions:VNS1: in each iteration roulettewheel selection of routes in neighborhood structures and order of LS operators with the first improvement strategy, CA acceptance functionVNS2: after new solution additional nloc iterations of LS with only one operator in each iteration, RR acceptance functionVNS3: nneig for applying each neighborhood structure, LS, acceptance of only better solutionsLS: 2-Opt∗, RemoveTwoInsertOne, Relocate, 4-Opt, Remove/Insert BatteryStation(3) Insert s new solutions in the population, the rest of the solutions are generated based on the destroy and repair heuristic (2-regret)  Paz et al. [86] MDEVLRPTW - BS/PR/BSPR     CPLEX Modified E-LRPTW - small  Pelletier et al. [87] EFV-CSP     CPLEX Generated The focus is on the deriving meaningful insights and not the computation prowess  Poonthalir and Nadarajan [38] F-GVRP  NNH + random TVa-PSOGMO   GVRP Particle encoding: customer sequence only which is converted to solution by inserting depots and refueling stationsPersonal and global best particles + velocity updateTime-varying inertia and acceleration coefficientsDecoding: rank of value to convert the route to the particle positionGMO and archive of Pareto optimal solutions  Schiffer and Walther [88] LRPIF X MCWS ALNS+LS+DP  DP E-LRPTWPR (old and new), BSS-EV-LRP, E-VRPTW, E-VRPTWPR Penalty terms evaluation - corridor-based approach - a range of possible refueling times + concatenation operatorsALNS and LS: Schiffer et al. [11]DP - optimal place of recharging visits - solve resource constrained SPP for each route by extending the given routes only for facility vertices + resource extension functions for feasibility and dominance checks  Schiffer and Walther [89] RELRPTWPR X MCWS ALNS+LS+DP  DP, Gurobi - small Real Parallelized ALNS of Schiffer and Walther [88]: LS, DP and customer configuration performed in parallel for fixed CS configuration, while CSs configuration, adaptive evaluations and stopping criterion are executed within a critical section of the code Adversarial approach for robust formulation  Shao et al. [90] E-VRP   GA+LS   Real General GA + elitismLS: random Exchange  Wang et al. [91] BEV routing    GA + fuzzy prog.  Generated Fuzzy programming and fuzzy preference relations are applied to transform the objective function into a single objective functionGA: chromosome genes - CSs  Zhang et al. [36] E-VRP X Construction + insertion heuristics AC+ILS, ALNS   Generated Initial:(1) TSP - construction based on the savings and intensity trail value, (2) VRP - insertion of depots, (3) E-VRP - insertion of CSsILS perturbation: random removal of customers and CSs, and then the best insertion procedureLS: 2-Opt∗, Exchange, Relocate, StationInRe ElitismALNS of Goeke and Schneider [30] with modf.: no cluster removal operator  Basso et al. [92] 2sEVRP     CPLEX Generated Bellman-Ford algorithm for SPP  Breunig et al. [93] E2EVRP  Greedy insertion LNS+LS+DP  Exact - route and candidate solution enumerations E2EVRP LNS destroy: RelatedRemoval, Random, CloseSatellite, OpenSatellite, RemoveSingleCustomerRoutesLNS repair: (i) the reinsertion of customers in the second echelon with Greedy (random or demand quantity) operators; (ii) greedy insertion on the first echelon routes(iii) DP for CSs insertionsLS on second echelon: 2-Opt, 2-Opt∗, Relocate, Swap, Swap2-1  Bruglieri et al. [94] GVRP     Exact + Xpress GVRP Path-based approach - the route is a composition of paths for a subset of customers without station visits:(1) All feasible paths are generated, removing from set of the dominated ones(2) MILP formulation is used to combine paths in routes  Froger et al. [95] E-VRP-NL   Heuristic for FRVCP  Exact labeling - FRVCP, GUROBI - small E-VRP-NL Heuristic for FRVCP: at each CS, EV is charged minimum energy required to reach next CS or depot  Hiermann et al. [96] H 2E-FTW X  GA+LS+LNS+DP  CPLEX - set partitioning H 2E-FTW, E-FSMFTW, E-VRPTW, E-VRPTWPR The sequence of customers and vehicle types without CS(1) GA - assignment and sequencing of customer visits - 2 populations: feasible and infeasible(2) DP - sequencing of CS: resource CSPP - labeling algorithm for CSs insertions(3) Greedy extension policy - known sequence, partial recharging, mode selection Route evaluation: (2)&(3), extension functionsGA crossover: (1) binary tournament selection, (2) each parent forms giant TSP tour, (3) OX crossover, (4) split procedureGA mutation: LNS - random selection of destroy operators: Random, RandomRoutes, Similar, Target; uniform selection of repair operators: Greedy, Regret-2LS (only promising edges): 2-opt, 2-opt∗, Swap, Relocate (single and pairs), VehicleSwapSet partitioning over a pool of routes during the searchThe cache memory of move evaluations for unchanged route sequences  Jie et al. [97] 2E-EVRP-BSS X Extended sweep CG-ALNS SA CPLEX - small 2E-EVRP-BSS, 2E-VRP ( 1 ) Second echelon - ALNS: destroy: Random, Related, SatelliteClosing, SatelliteOpening, SatelliteSwapping, RandomRoute, BSS-CustomerRoute, BSS-ProportionRoute, BSS-CostBased, BSS-BestUse; repair: basic and advance Greedy, basic and advance Regret-(2,3)(2) First echelon: column generation (CG) for split delivery EV routing problem, DP and branching scheme for the SPP with resource constraint sub-problem(3) Label-based BSS optimization  Koyuncu and Yavuz [98] MGVRP     CPLEX GVRP Upper and lower bound for energy state (labeling algorithms) and arrival times Lower bound on the number of vehicles  Macrina et al. [99] GMFVRP-PRTW  CPLEX Hybrid LNS   GMFVRP-PRTW Initial: CPLEX to find an initial feasible solutionRandom order of destroy and repair operators and tabu list to forbid the removal of recently removed customers and reinsertion into routes from which they were removedDestroy: Random, WorstDistance, WorstTime, RouteInsertion: Greedy, GreedyAndCS, GreedyNewRoute, GRASP  Macrina et al. [33] GMFVRP-PRTW  Sequential insertion heuristics ILS  CPLEX - small E-VRPTW Initial: two clusters of customers for EV and ICEV routes -(1) ICEV routes, (2) EV routes - all customers from EV cluster and all un-routed customers from ICEV cluster LS with/without penalty function and perturbation operators: relocate-only EV, ICEV or both  Normasari et al. [100] CGVRP  NNH SA+LS SA CPLEX CGVRP, GVRP Random neighborhood structures: Swap, Insertion, station Insert/Delete Reference: referenced paper; Problem name: the name of the analyzed problem; NF: allowing infeasible solutions during the search process; Initial: procedure used for the creation of the initial solution if any; Hybrid heuristics: heuristics, metaheuristics, and hybrid combinations applied for solving the problem; A: the acceptance criteria; Exact & Software: the exact procedure or software applied on the analyzed problem or subproblem; Instances: instances used for testing of the applied procedures where Real denotes the real case-study instances, Generated denotes the artificially created instances which are not benchmark instances, and other denotes the benchmark instances of the problem; Description: short description of procedure used to solve the analyzed problem.


### 5.1. Feasibility

When creating or modifying a VRP solution, two types of search procedures regarding the feasibility of the solution can be observed: allowing only feasible solutions or acceptance of infeasible solutions. The infeasible solution means that some customers are served without satisfying all of the problem constraints. Related, feasible procedures are searching in the space of feasible solutions, while infeasible procedures allow searching the space of infeasible solutions, which broadens the search. In the E-VRP, the feasibility mostly refers to vehicle load, customer time window, and battery (energy) constraints. In infeasible procedures, objective function is often defined with penalization coefficients, which are updated during the search process. At the beginning of the search, infeasible solutions are allowed in order to search a larger solution space. As the search process comes to an end, penalties for infeasible solutions increase. An example of such an objective function is given by (8) where fdist(S) is the total traveled distance in the solution S; penalties coefficients are α, β, and γ; Pdiv(S) is a diversification penalty; and Pcap(S), Ptw(S), and Pbatt(S) are values of constraint violation, respectively, load capacity, time windows, and battery capacity [25]. An important part of such a function is the efficient computation of constraint violations [25, 31, 96, 140].(8)fS=fdistS+αPcapS+βPtwS+γPbattS+PdivSSeveral researchers in the E-VRP field are allowing infeasible solutions during the search [11, 25, 30, 31, 69, 84, 88, 96] while some others are only allowing feasible solutions [23, 28, 46–49, 52, 68, 72, 79]. For example, in the E-VRPTWPR of Keskin and Çatay [28], the customer is inserted in the solution if it satisfies the load and partial time window constraints. If the solution is energy infeasible the CS is inserted to make the solution completely feasible. If the former one was not able to produce the energy feasible solution, the procedure returns to the latest feasible solution. Similarly, Montoya et al. [63] and Froger et al. [68] address the fuel infeasible route by solving a CSPP on the auxiliary graph. A list of papers allowing infeasible or feasible-only solutions is presented in Table 2 in columnNF.


### 5.2. Exact Procedures

Exact procedures are able to find an optimal solution for a smaller number of customers: up to 360 customers for CVRP and up to 100 customers for VRPTW [5, 6]. Some of the most used exact algorithms are branch-bound-cut-and-price, dynamic programming (DP), MIP, set partitioning (SP), etc. [1]. In the E-VRP, many researchers formulate the problem as a MIP and solve the small instances with commercially available software. Koç and Karaoglan [59] applied the branch-and-cut algorithm coupled with simulated annealing and four improvement strategies to solve the GVRP. The authors reported that 22 out of 40 instances were solved to optimality. Hiermann et al. [31] solved the heterogeneous E-VRP by applying a branch-and-price algorithm on small instances up to 15 customers where the pricing problem was solved using a bidirectional labeling algorithm. The algorithm was able to find an optimal solution for all of the instances in a couple of minutes, with only three instances having computation time longer than five minutes. Desaulniers et al. [57] presented branch-price-and-cut algorithms for four E-VRPTW variants depending on the full or partial recharge and on the single or multiple recharges per route. On the 696 instances with 25, 50, and 100 customers, 98%, 90%, and 27%, respectively, were solved optimally. Schiffer et al. [11, 27], Schiffer and Walther [88], and Hiermann et al. [96] applied DP on resource CSPP (RCSPP) to find the optimal facility route configuration in order to enhance the solutions produced by metaheuristic, while Pourazarm et al. [54] applied DP to find the fastest BEV routes in the single and multivehicle environment. Similarly, Montoya et al. [63] used pulse algorithm to determine the optimal placement of refueling facilities in a route with fixed customers. Similar to the FRVCP formulation of Montoya [18], Keskin and Çatay [79] in their E-VRPTW-FC formulated the MILP program to optimize charging related decisions in BEV route with the fixed customer positions. For every Ω iteration, the authors tried to solve this problem exactly by applying CPLEX solver. Liao et al. [103] applied a branch-and-bound algorithm to solve the EV touring problem. Montoya et al. [63, 72] and Hiermann et al. [96] used set partitioning on the pool of routes to select the best subset of routes guaranteeing no overlapping in customers assignments. Froger et al. [68] used a branch-and-bound algorithm to solve the set partitioning model and to discard the selection of routes that are infeasible or for which the total time is underestimated.


### 5.3. Heuristic Procedures

Heuristic procedures seek to solve the problem based on the specific knowledge of the problem, usually suboptimal or close enough to a satisfactory solution. In the research field of VRP, heuristic procedures can be split into constructive and improvement heuristics.


#### 5.3.1. Constructive Heuristics

Constructive heuristics are often used to generate an initial solution by either serial or parallel route construction. Solutions are constructed in a greedy way, which often produces solutions of the VRP that are 10-15% far from an optimal solution [4]. In the E-VRP, constructive heuristics are modified and adapted to the BEV characteristics and feasibility checks. Some well-known general constructive heuristics used in E-VRP are presented.Themodified savings method of Clarke and Wright (MCWS) [141] with the insertion of CSs is used to generate an initial solution of several E-VRP and related problems [23, 32, 52, 69]. The algorithm starts with the creation of back-and-forth routes for each customer. Next, CSs are inserted in a least-cost manner, in either feasible or infeasible routes. Then, the routes that increase the savings the most are merged together and the procedure is repeated until all customers are served.Thesweep algorithm [142] inserts customers in the active route in a circular manner, resulting in efficient space division. The customers are sorted based on the value of a polar angle between the depot and randomly chosen point. Then, customers are inserted in the active route at the position causing minimal increase until a violation of route constraint, when usually a new route is opened [25, 56, 97].Nearest neighbor heuristic (NNH) starts from the depot, and in each iteration the customer with the least-cost increase from the previously selected customers is added to the route. The route is terminated when some/any constraints are violated and a new route is opened [76].Route-first cluster-second [143] approach constructs a giant infeasible TSP tour, usually with NNH, which is then split into several routes. The split procedure can be solved as the SPP on the acyclic graph. The procedure generates energy infeasible solutions, so the CS placement and charging schedule need to be determined to make the solution energy feasible [63, 68, 72].Thecluster-first route-second [144] approach divides the customers into clusters and then each cluster is solved as a subrouting problem. One example used in GVRP is thedensity based clustering algorithm (DBCA) proposed by Erdoğan and Miller-Hooks [23], which determines customer clusters based on their density, by ensuring that each customer has at least some predefined number of customers in the ϵ radial neighborhood. In the second step to obtain the routes, the MCWS is performed.Many researchers are adopting some form of the insertion heuristics, which have some characteristics of already reviewed heuristics. The general idea of the insertion heuristics is to iteratively insert or add customers in available routes. In each step of the algorithm, an unserved customer together with a route and its position in the route are determined to least increase the objective function. This selection can be deterministic but is often stochastic to select at random one fromk-best inserts to diversify the solution, especially in population metaheuristics. In most feasibility ensured cases when the energy constraint is violated, a CS is added to the route to make it energy feasible, and when the other constraints are violated such as vehicle load capacity or time windows, the new route is opened, and the procedure is repeated. Examples of adaptations to the E-VRP arepush forward insertion heuristics (PIH) [41, 145],constructive-k-PseudoGreedy [46], [28, 30, 31, 33, 36] modified insertion heuristics, andcharging routing heuristic (CRH) [47–49], in which first the charging scheme at the depot is determined, and then the least-cost feasible customer insertions are performed.


#### 5.3.2. Improvement Heuristics

Improvement heuristics or local search (LS) explores the neighborhood of the incumbent solution, searching for a better solution. The neighborhood is explored by applying perturbation moves based on the composite neighborhood operators. The local search stops when no improving solution can be found in the neighborhood of the incumbent solution, which is then called local optima. A great variety of researchers have performed LS procedures to intensify the search, coupled together with perturbation moves to escape the local optima. Often, the perturbation moves are similar to the neighborhood operators used in LS phase. Most of the classical VRP neighborhood operators [3, 4] are used in the E-VRP, but some additional neighborhood operators have been developed. Depending on whether the operators perform on only one route or between the routes, they can be divided intointraroute andinterroute operators. Here, we present some of the most used LS neighborhood operators in the E-VRP literature:(i)
2-Opt [146] replaces two arcs with the two new ones, with a possibility of route direction reversal.(ii)
2-Opt∗ [147], unlike the 2-Opt operator, avoids the reversal of a route direction.(iii)
4-Opt [84, 148] replaces four consecutive arcs with four new ones.(iv)
Or-Opt [149] replaces three arcs with three new ones such that a sequence of three vertices is relocated.(v)
Exchange (Swap) [150] exchanges the two vertices in the solution.(vi)
FacilityInsertion [55] andStationInsertion [11] operators insert refueling or loading facility at best position in between two consecutive facility visits where fuel or load violation occurred. The similar operatorInsertRemoveIf [31] inserts CS into the part of the route where the violation occurred to make the route energy feasible and repair violations, but it also removes redundant CS visits. Schneider et al. [55] also proposed the facility-related removal and exchange operatorsFacilityReplacement andExchange.(vii)
RechargeRelocation [46] removes all visits to the CSs and then tries to insert CSs at positions leading to a better solution. Similarly, theGlobalChargingImprovement [72] operator removes all CS visits in the route, and if the route is energy feasible, it stops; otherwise, it solves FRVCP to optimize CSs positions and charging amount.(viii)
Relocate [150] removes one vertex from the solution and inserts it into another position in the solution. Felipe et al. [46] proposed a modification of the relocation operator calledReinsertion,Rbu, for solving GVRP-MTPR. Binary parameters b and u denote whether the first-improve or best-improve strategy is used and whether the values of savings are updated or not. TheRemoveTwoInsertOne [84] operator removes two vertices at random and inserts them one by one into another vehicle at best possible positions.(ix)
CrossExchange [151] removes n1 consecutive vertices from route r1 to route r2, and n2 consecutive vertices from route r2 to route r1 [84]. TheRemoveSequence [84] operator is similar as it removes consecutive vertices from only one route and inserts it into another one.(x)
Resize [31],RelocateAndResize [31], χ-OptWithModeChange [58], andVehicleSwap [96] change the vehicle type assignment or the propulsion mode in the heterogeneous E-VRP and HEV routing problems.(xi)
StationInRe [25] inserts and removes CSs from the solution in so that if an arc to the CS is already in the solution, it is removed; otherwise it is inserted.The order of the operators affects both the solution quality and the execution time. Often, there is a question whether to make afirst better orthe best improvement move in the LS phase [152]. Some E-VRP researchers perform one or the other, and some combine the two approaches, i.e., best move of first 100 [11] or first 50 [31].Neighborhood operators explore only the space of the immediate vicinity of the current solution, which often leads to local optima. In order to search a larger solution space, Shaw [153] proposed thelarge neighborhood search (LNS) heuristic, which is based on destroying and repairing the solution. Heuristic efficiency depends on the implemented destroy and, especially, repair operators. The main drawback of the heuristic is the repeated use of the same destroy and repair operators, which can lead to local optima. A modification of LNS by Christiaens and Vanden Berghe [154] is giving state-of-the-art results on a wide range of VRP variants. In the E-VRP, LNS is successfully implemented by Sassi et al. [47–49] and Zhang et al. [36], but most often it still needs some other guiding heuristic to escape the local optima. Macrina et al. [99] applied an LNS scheme to solve GMFVRP-PRTW. A list of used destroy and repair operators in the E-VRP is presented as a part of the adaptive LNS version in the appendix.


#### 5.3.3. Route or Move Evaluation

Operator evaluation is an important part that highly affects execution time and has to be efficiently implemented. For example, in VRPTW, for most of the basic neighborhood operators, the load capacity check and time window check can be computed in constant timeO(1), if for each customer endogenous variables, i.e., remaining load capacity, latest arrival time, etc., are tracked. In the E-VRP, there is also a constraint regarding the vehicle battery capacity (battery feasibility), which also influences the time window feasibility. Therefore, time window evaluation is more complex. For example, when the customer is trying to be inserted into the route, the load and energy check can be computed in constant time O(1) for most of the basic neighborhood operators. The endogenous variables for vertices between the insertion position and next CS have to be recomputed in order to know the charging time at next CS. The complexity of such approach is O(B) where B is the number of customers between the position of the insertion and next CS. If there is no CS in the second part of the route, no recalculation is needed. By storing additional information about the recharging position and its amount and by the use of general concatenation operators presented by Vidal et al. [155], Hiermann et al. [31] reported that all of the basic neighborhood operators can be evaluated in O(1). Goeke and Schneider [30], instead of storing additional information about recharging, used a surrogate cost function whose evaluation is less demanding. In the search process, the authors store ϵ best moves in the composite neighborhood and evaluated either exact or surrogate cost, and at the end if there is any move with the surrogate cost, its exact cost is evaluated and the best move is chosen. Hiermann et al. [96] presented a different approach for H2E-FTW in which customer sequence without CSs is determined, and then in the route evaluation part time-consuming CSs insertions and charging schedule are performed. Whenever the route needs to be evaluated, the heuristic finds the routing cost by applying a labeling algorithm to determine CSs placement and a greedy algorithm to determine charging amount and running mode of PHEV. As such an evaluation is time-consuming, the authors performed the following: (i) use neighborhood restrictions to consider only promising arcs for each customer vertex; (ii) keep track of move evaluations in cache memory; and (iii) lower bounds to first evaluate the move as ICEV in O(1), and then if it improves the objective function, perform the route evaluation.


### 5.4. Hybrid Metaheuristic Procedures

Many researchers employ metaheuristics to continue the exploration after the first local optima occurrence. Metaheuristics can be defined as heuristics guiding other heuristics and can be divided into neighborhood-oriented metaheuristics and population metaheuristics. Most often, multiple metaheuristics and heuristics are combined together and adapted to the problem; therefore, the term hybrid is used in such occasions.


#### 5.4.1. Neighborhood-Oriented Metaheuristics

Neighborhood-oriented heuristics iteratively explore the neighborhood of the incumbent solution. Here, we present the most used neighborhood heuristics in the E-VRP literature.Simulated annealing (SA) [156] controls the search by the temperature parameter, which mimics the physical cooling process of the material. Most often, it is used as a control heuristic for accepting solutions produced by other heuristics. The probability of accepting solution S2, which is worse than the current best solution S1, is computed as e(f(S1)-f(S2))/T (Boltzmann function [157]), where T is the temperature parameter. The higher the temperature value, the higher is the probability of accepting a deteriorating move to escape the local optima. The procedure starts with high temperatures at the beginning and then dynamically lowers the temperature as the search process comes to an end. In the E-VRP, most often initial temperature is selected such that a μ% worse solution than the initial solution is accepted with a probability of 0.5 [25, 28, 30, 46].Tabu search (TS) [158, 159] is a well-known metaheuristic, which uses memory to forbid the exploration of the space that has already been explored. It escapes the local optima by accepting a worse solution if it is the best in the explored space. Arcs that are deleted from the solution are stored in thetabu-list, which prohibits the reinsertion of the deleted arcs into a specific part of the solution for a designated number of iterations,tabu-tenure. In some implementation cases, when the new best solution is found, some of the moves and arcs are lifted and removed from thetabu-list,aspire criterion [25]. In the E-VRP, TS with neighborhood operators is used to intensify the search [25, 35, 48].Variable neighborhood search (VNS) [160] is a metaheuristic that systematically changes the neighborhood each time there is no improvement in LS phase. The changes are based on predefined neighborhood structures. In each step, the new neighborhood operator is used, which increases the vicinity from the incumbent solution. The selection of neighborhood operators and customers can be deterministic (variable neighborhood decent (VND)) or probabilistic. Schneider et al. [25] successfully adapted it in the E-VRPTW for diversification of the solution based on the cyclic exchange moves, resulting in 15 different neighborhood structures. The VNS in the E-VRP was also combined with exact algorithms like branch-and-price [51], evolutionary algorithms [84], deterministic LS operators [64, 72], and adaptive mechanisms (AVNS) [11, 27, 55, 69].Iterated local search (ILS) [161] is based on successively repeating LS on the incumbent solution. When LS ends up in local optima, perturbation move is applied to escape from local optima. The perturbation move can be scaled dynamically to overcome local optima. The effectiveness of the procedure highly depends on the local search procedures used. Montoya et al. [72] applied ILS together with VND for solving E-VRP-NL with perturbation performed by splitting perturbed giant TSP tour. Macrina et al. [33] used typical ILS with some of the LS operators used in the perturbation phase to solve GMFVRP-PRTW. Zhang et al. [36] used ILS and LNS in an ant colony framework for solving GVRP. The LNS is used as a perturbation mechanism while the LS is used to further improve the solution. Sassi et al. [47, 49] also used a multistart combination of ILS and LNS but without the LS phase. Sassi et al. [48] applied similar ILS paradigm, but they called ititerated tabu search (ITS), as they used TS for the intensification and LNS for the perturbation. A similar ITS procedure is applied by Doppstadt et al. [58] for the HEV-TSP.Adaptive large neighborhood search (ALNS) [162, 163] is the extension of LNS heuristic in which, throughout the search process, different procedures for destroying and repairing the solution are selected in an adaptive manner based on their past performance. The idea is to have versatile destroy and repair operators in order to escape the local optima. For each destroy and repair operator, the adaptive weight and score are increased in the following order (scores σ1≥σ3≥σ2): (i) σ1 if a new best solution is found; (ii) σ2 if the solution is better than the previous one; (iii) σ3 if the solution is worse than the previous one but it is accepted due to the acceptance criteria; and (iv) no increase if the solution is not accepted. The score of worse solution σ3 is greater than the score of better solution σ2, to reward unimproved solutions that diversify the search [28]. For each N iteration, the weights are updated based on the previous weight values and new score values. As a great number of the researchers used ALNS for solving the E-VRP and related problems (Table 2), a vast number of destroy and repair operators were applied. Some of them originated from Pisinger and Ropke [162] but also some new problem-specific ones regarding the CS have been developed [28, 30–32, 79, 88]. The description of the applied destroy and repair operators in E-VRP is presented in the appendix. Schiffer and Walther [89] applied a parallel ALNS version of Schiffer and Walther [88] to solve the RELRPTWPR. The OpenMP was used to solve multiple scenarios in parallel, where the customer destroy and repair operators, as well as the LS and DP, are performed in parallel. The change of facility configuration, adaptive evaluations, and stopping criterion are not parallelized.SA is one of the most used metaheuristics for acceptance criteria, in order to escape the local optima, but some others can also be applied. Tiwari et al. [164] proposed the Cauchy function (CA) T/((f(S1)-f(S2))2+T2) for the acceptance of the solution as it gives more opportunities to escape from local optima. The record-to-record (RR) principle [165], which adds negative value to the objective function of the best solution so far, can also be used to escape the local optima. Several papers in the E-VRP field applied such heuristics for solution acceptance [48, 49, 84].


#### 5.4.2. Population Metaheuristics

Population metaheuristics are based on the natural selection to evolve a population and survival of the fittest. They have been widely applied in the VRP field: genetic algorithm [166], path relinking [167], scatter search [168], ant colony [169], bee colony [170], particle swarm optimization [171, 172], etc.Genetic algorithm (GA) consists of a pool of individuals (VRP solutions) called a population, which goes through the process of evolution. Such evolution process includes evaluation of individuals, selection of parents, crossover, mutation, and a replacement of the old population by a new one. The important aspects of the algorithm are the representation of the VRP solution, efficient management of population diversity, and fitness function to compare the individuals in the population. For further algorithm description and application, we refer to Potvin [173]. Theant colony (AC) algorithm is based on ant behavior when they search for food. As they search, they leave pheromones, which accumulate over the most used path, leading to the optimal path between the food and the nest. AC algorithm consists of three procedures: pheromone initialization, probability adjustment, and pheromone update.Particle swarm optimization (PSO) mimics the behavior of organisms such as fish schooling or bird flocking. Physical movement of the individual (particle) in the swarm is tracked and then well-balanced moves are applied to adapt the individual to the global and local exploration. Each particle has a position in the swarm and corresponding velocity in order to optimize the problem following personal best solution and global best solution [172].Barco et al. [34, 42] applieddifferential evolution to solve the E-VRP and charge scheduling. The initial individual represents the valid route assignment and charging profile. Next, mutation based on OR and XOR operators is performed on three randomly selected individuals. The crossover is a random 2-point while the selection chooses between the target individual and new individual, which is generated by the crossover. Alesiani and Maslekar [45] applied GA for solving a single BEV routing problem. The individual is represented as a two-dimensional binary array, where the rows represent the routes and the columns represent CSs. The initial population contains feasible individuals, individuals without CSs, and individuals containing battery CSPP routes. The authors select parents at random or based on the probability and apply 2-point crossover. Mutation is applied at random to remove, insert, or exchange CS, while the LS is performed to improve the solution. Shao et al. [74] applied GA without detailed description for solving real-life EVRP-CTVTT instances. The authors encoded the individual with all visited vertices including customers, depots, and CSs. The crossover operator tries to create a new individual that satisfies all the constraints. Also, elitism is applied to preserve the elite individuals. The authors reported the average computational time of three hours. A similar algorithm with the LS after the crossover and mutation is applied for solving E-VRP by Shao et al. [90] to obtain an optimal operation scheme including routes, charging plan, and driving paths. Masliakova [61] applied the GA to find the optimal placement of CSs for electric bus routes. The individual is represented as the sequence of CSs and bus stops visited, and a crossover is implemented as the swap between the subroutes. Mutation moves the CS to the nearest bus stop or creates a new bus stop with CS. For the fitness function, the author used an AC approach, which ranks the individuals in the population. Joo and Lim [78] used the AC approach to determine the energy consumption and speed in energy-efficient BEV routing. Each ant selects the next vertex with a greedy stochastic search, which combines heuristic information about energy consumption, speed, and pheromones. The authors concluded that AC strategy is suitable for routing BEVs in terms of energy efficiency. An efficient version of evolutionary algorithm and VNS (EVO-VNS) is applied by Masmoudi et al. [84] for solving the DARP-EV. Each population is divided into m groups in such a way that each group contains good and bad solutions. Further on, s solutions from each group are selected in the VNS phase. If the best solution in the current group is not improved for three iterations, the crossover between the selected solutions and the rest of solutions in the group is performed. The rest of the individuals in the population are generated by applying destroy-repair procedure with regret heuristic. Poonthalir and Nadarajan [38] applied particle swarm optimization with greedy mutation operator and time-varying acceleration coefficient (TVa-PSOGMO) to solve the F-GVRP. Each particle is encoded as the sequence of customers who are converted to feasible routes by inserting depots and refueling stations. Time-varying inertia and acceleration coefficients are employed to better manage diversification and intensification exploration in the multiobjective environment. Additionally, to prevent premature convergence, a greedy mutation operator (GMO) is performed. Hiermann et al. [96] used hybrid GA (HGA) with LS to obtain good assignment and routing decisions for H2E-FTW. The authors represented the solution as a set of vehicle routes and used a biased fitness function. During the search, the authors maintained two populations of feasible and infeasible individuals. Parents are selected using binary tournament scheme from those two populations. The recombination is performed on the parent routes combined in a giant single tour, with an OX crossover and a split algorithm to obtain a complete solution. The mutation phase is performed as LNS with destroy and repair operators similar to ones of Hiermann et al. [31]. In the end, the solution is improved by the LS procedure.


### 5.5. Summary of the Applied Procedures

The overview of the applied procedures in the E-VRP field is presented in Table2, where 78 papers are presented. Searching in the space of infeasible solutions is applied in 17 papers, while the LS for the intensification is applied in 25 papers. Most of the researchers do not use any special criteria for accepting the solution, but 17 papers use SA as acceptance criteria, three RR and one CA. For the initial solution, mostly some form of the universal insertion heuristic is applied (18). Then follow MCWS (11), NNH (7), random (5), split procedure (4), and sweep (3). Regarding the applied heuristics for solving the problem, the order is the following: ALNS (14), GA (10), LNS (7), TS and ITS (6), ILS (6), AC (3), VNS (3), AVNS (2), VND (2), PSO (1), and differential evolution (1). From all of the exact procedures, the DP is the most applied one (12) and then there is the branch-bound-cut-and-price (5). Forty-three papers use commercially available software like CPLEX (32), Gurobi (7), MATLAB (2), and XPRESS (2) for solving either the whole problem or some subproblem during the search, mostly charging decisions.We can note that a vast number of researchers use hybrid heuristics procedures with only a few focusing on the exact procedures [31, 57]. It can be noted that most of the procedures that are giving state-of-the-art results are using exact procedures at some point during the search, usually for (i) optimal determination of charging visits and charging amount [79, 88, 96] and (ii) set partitioning over the set of routes at the end of the search [65, 72, 95, 96].The ALNS is one of the most applied procedures and it is giving the BKSs on various problem variants. Among them, we can point out the following: ALNS of Hiermann et al. [31] (HPH) and ALNS of Goeke and Schneider [30] (GS) for heterogeneous fleet, ALNS of Keskin and Çatay [28] (KC1) and ANLS of Keskin and Çatay [79] (KC2) for partial recharges and different charging technologies, and ALNS of Schiffer and Walther [88] (SW) originally designed for the LRPIF (ELRPTWPR) which is also giving BKSs on the E-VRPTW and E-VRPTWPR. Schiffer and Walther [88] compared their algorithm on the E-VRPTW instances to the HPH, GS, KC1, and VNS-TS of Schneider et al. [25] (SSG). The authors found nine BKSs with an average gap of 0.27% to the previous BKSs and average running time of 3.77 minutes. On the E-VRPTWPR, compared to the KC1, the authors found 41 new BKSs. The proposed ALNS of SW is efficiently parallelized to solve the RELRPTWPR [89].Most of the population-based metaheuristics are not efficiently implemented so they are not giving high-quality solutions. But, the HGA of Hiermann et al. [96] (HHPV) gives some of the BKSs on the E-VRPTW, E-VRPTWPR, and H2E-FTW. The authors compared their results on the E-VRPTW and E-VRPTWPR to the SSG, HPH, GS, and KC1. On the E-VRPTW instances, the HHPV found 11 BKSs with 1.36% increase in the vehicle number, 0.46% distance gap from the previous BKSs, and an average running time of 11.59 minutes. On the E-VRPTWPR, the authors found 28 new BKSs. By comparing the results of HHPV and SW on E-VRPTW, we can conclude that in 20 instances SW outperformed HHPV, while in the opposite case the HHPV was better in 17 instances.


## 5.1. Feasibility

When creating or modifying a VRP solution, two types of search procedures regarding the feasibility of the solution can be observed: allowing only feasible solutions or acceptance of infeasible solutions. The infeasible solution means that some customers are served without satisfying all of the problem constraints. Related, feasible procedures are searching in the space of feasible solutions, while infeasible procedures allow searching the space of infeasible solutions, which broadens the search. In the E-VRP, the feasibility mostly refers to vehicle load, customer time window, and battery (energy) constraints. In infeasible procedures, objective function is often defined with penalization coefficients, which are updated during the search process. At the beginning of the search, infeasible solutions are allowed in order to search a larger solution space. As the search process comes to an end, penalties for infeasible solutions increase. An example of such an objective function is given by (8) where fdist(S) is the total traveled distance in the solution S; penalties coefficients are α, β, and γ; Pdiv(S) is a diversification penalty; and Pcap(S), Ptw(S), and Pbatt(S) are values of constraint violation, respectively, load capacity, time windows, and battery capacity [25]. An important part of such a function is the efficient computation of constraint violations [25, 31, 96, 140].(8)fS=fdistS+αPcapS+βPtwS+γPbattS+PdivSSeveral researchers in the E-VRP field are allowing infeasible solutions during the search [11, 25, 30, 31, 69, 84, 88, 96] while some others are only allowing feasible solutions [23, 28, 46–49, 52, 68, 72, 79]. For example, in the E-VRPTWPR of Keskin and Çatay [28], the customer is inserted in the solution if it satisfies the load and partial time window constraints. If the solution is energy infeasible the CS is inserted to make the solution completely feasible. If the former one was not able to produce the energy feasible solution, the procedure returns to the latest feasible solution. Similarly, Montoya et al. [63] and Froger et al. [68] address the fuel infeasible route by solving a CSPP on the auxiliary graph. A list of papers allowing infeasible or feasible-only solutions is presented in Table 2 in columnNF.


## 5.2. Exact Procedures

Exact procedures are able to find an optimal solution for a smaller number of customers: up to 360 customers for CVRP and up to 100 customers for VRPTW [5, 6]. Some of the most used exact algorithms are branch-bound-cut-and-price, dynamic programming (DP), MIP, set partitioning (SP), etc. [1]. In the E-VRP, many researchers formulate the problem as a MIP and solve the small instances with commercially available software. Koç and Karaoglan [59] applied the branch-and-cut algorithm coupled with simulated annealing and four improvement strategies to solve the GVRP. The authors reported that 22 out of 40 instances were solved to optimality. Hiermann et al. [31] solved the heterogeneous E-VRP by applying a branch-and-price algorithm on small instances up to 15 customers where the pricing problem was solved using a bidirectional labeling algorithm. The algorithm was able to find an optimal solution for all of the instances in a couple of minutes, with only three instances having computation time longer than five minutes. Desaulniers et al. [57] presented branch-price-and-cut algorithms for four E-VRPTW variants depending on the full or partial recharge and on the single or multiple recharges per route. On the 696 instances with 25, 50, and 100 customers, 98%, 90%, and 27%, respectively, were solved optimally. Schiffer et al. [11, 27], Schiffer and Walther [88], and Hiermann et al. [96] applied DP on resource CSPP (RCSPP) to find the optimal facility route configuration in order to enhance the solutions produced by metaheuristic, while Pourazarm et al. [54] applied DP to find the fastest BEV routes in the single and multivehicle environment. Similarly, Montoya et al. [63] used pulse algorithm to determine the optimal placement of refueling facilities in a route with fixed customers. Similar to the FRVCP formulation of Montoya [18], Keskin and Çatay [79] in their E-VRPTW-FC formulated the MILP program to optimize charging related decisions in BEV route with the fixed customer positions. For every Ω iteration, the authors tried to solve this problem exactly by applying CPLEX solver. Liao et al. [103] applied a branch-and-bound algorithm to solve the EV touring problem. Montoya et al. [63, 72] and Hiermann et al. [96] used set partitioning on the pool of routes to select the best subset of routes guaranteeing no overlapping in customers assignments. Froger et al. [68] used a branch-and-bound algorithm to solve the set partitioning model and to discard the selection of routes that are infeasible or for which the total time is underestimated.


## 5.3. Heuristic Procedures

Heuristic procedures seek to solve the problem based on the specific knowledge of the problem, usually suboptimal or close enough to a satisfactory solution. In the research field of VRP, heuristic procedures can be split into constructive and improvement heuristics.


### 5.3.1. Constructive Heuristics

Constructive heuristics are often used to generate an initial solution by either serial or parallel route construction. Solutions are constructed in a greedy way, which often produces solutions of the VRP that are 10-15% far from an optimal solution [4]. In the E-VRP, constructive heuristics are modified and adapted to the BEV characteristics and feasibility checks. Some well-known general constructive heuristics used in E-VRP are presented.Themodified savings method of Clarke and Wright (MCWS) [141] with the insertion of CSs is used to generate an initial solution of several E-VRP and related problems [23, 32, 52, 69]. The algorithm starts with the creation of back-and-forth routes for each customer. Next, CSs are inserted in a least-cost manner, in either feasible or infeasible routes. Then, the routes that increase the savings the most are merged together and the procedure is repeated until all customers are served.Thesweep algorithm [142] inserts customers in the active route in a circular manner, resulting in efficient space division. The customers are sorted based on the value of a polar angle between the depot and randomly chosen point. Then, customers are inserted in the active route at the position causing minimal increase until a violation of route constraint, when usually a new route is opened [25, 56, 97].Nearest neighbor heuristic (NNH) starts from the depot, and in each iteration the customer with the least-cost increase from the previously selected customers is added to the route. The route is terminated when some/any constraints are violated and a new route is opened [76].Route-first cluster-second [143] approach constructs a giant infeasible TSP tour, usually with NNH, which is then split into several routes. The split procedure can be solved as the SPP on the acyclic graph. The procedure generates energy infeasible solutions, so the CS placement and charging schedule need to be determined to make the solution energy feasible [63, 68, 72].Thecluster-first route-second [144] approach divides the customers into clusters and then each cluster is solved as a subrouting problem. One example used in GVRP is thedensity based clustering algorithm (DBCA) proposed by Erdoğan and Miller-Hooks [23], which determines customer clusters based on their density, by ensuring that each customer has at least some predefined number of customers in the ϵ radial neighborhood. In the second step to obtain the routes, the MCWS is performed.Many researchers are adopting some form of the insertion heuristics, which have some characteristics of already reviewed heuristics. The general idea of the insertion heuristics is to iteratively insert or add customers in available routes. In each step of the algorithm, an unserved customer together with a route and its position in the route are determined to least increase the objective function. This selection can be deterministic but is often stochastic to select at random one fromk-best inserts to diversify the solution, especially in population metaheuristics. In most feasibility ensured cases when the energy constraint is violated, a CS is added to the route to make it energy feasible, and when the other constraints are violated such as vehicle load capacity or time windows, the new route is opened, and the procedure is repeated. Examples of adaptations to the E-VRP arepush forward insertion heuristics (PIH) [41, 145],constructive-k-PseudoGreedy [46], [28, 30, 31, 33, 36] modified insertion heuristics, andcharging routing heuristic (CRH) [47–49], in which first the charging scheme at the depot is determined, and then the least-cost feasible customer insertions are performed.


### 5.3.2. Improvement Heuristics

Improvement heuristics or local search (LS) explores the neighborhood of the incumbent solution, searching for a better solution. The neighborhood is explored by applying perturbation moves based on the composite neighborhood operators. The local search stops when no improving solution can be found in the neighborhood of the incumbent solution, which is then called local optima. A great variety of researchers have performed LS procedures to intensify the search, coupled together with perturbation moves to escape the local optima. Often, the perturbation moves are similar to the neighborhood operators used in LS phase. Most of the classical VRP neighborhood operators [3, 4] are used in the E-VRP, but some additional neighborhood operators have been developed. Depending on whether the operators perform on only one route or between the routes, they can be divided intointraroute andinterroute operators. Here, we present some of the most used LS neighborhood operators in the E-VRP literature:(i)
2-Opt [146] replaces two arcs with the two new ones, with a possibility of route direction reversal.(ii)
2-Opt∗ [147], unlike the 2-Opt operator, avoids the reversal of a route direction.(iii)
4-Opt [84, 148] replaces four consecutive arcs with four new ones.(iv)
Or-Opt [149] replaces three arcs with three new ones such that a sequence of three vertices is relocated.(v)
Exchange (Swap) [150] exchanges the two vertices in the solution.(vi)
FacilityInsertion [55] andStationInsertion [11] operators insert refueling or loading facility at best position in between two consecutive facility visits where fuel or load violation occurred. The similar operatorInsertRemoveIf [31] inserts CS into the part of the route where the violation occurred to make the route energy feasible and repair violations, but it also removes redundant CS visits. Schneider et al. [55] also proposed the facility-related removal and exchange operatorsFacilityReplacement andExchange.(vii)
RechargeRelocation [46] removes all visits to the CSs and then tries to insert CSs at positions leading to a better solution. Similarly, theGlobalChargingImprovement [72] operator removes all CS visits in the route, and if the route is energy feasible, it stops; otherwise, it solves FRVCP to optimize CSs positions and charging amount.(viii)
Relocate [150] removes one vertex from the solution and inserts it into another position in the solution. Felipe et al. [46] proposed a modification of the relocation operator calledReinsertion,Rbu, for solving GVRP-MTPR. Binary parameters b and u denote whether the first-improve or best-improve strategy is used and whether the values of savings are updated or not. TheRemoveTwoInsertOne [84] operator removes two vertices at random and inserts them one by one into another vehicle at best possible positions.(ix)
CrossExchange [151] removes n1 consecutive vertices from route r1 to route r2, and n2 consecutive vertices from route r2 to route r1 [84]. TheRemoveSequence [84] operator is similar as it removes consecutive vertices from only one route and inserts it into another one.(x)
Resize [31],RelocateAndResize [31], χ-OptWithModeChange [58], andVehicleSwap [96] change the vehicle type assignment or the propulsion mode in the heterogeneous E-VRP and HEV routing problems.(xi)
StationInRe [25] inserts and removes CSs from the solution in so that if an arc to the CS is already in the solution, it is removed; otherwise it is inserted.The order of the operators affects both the solution quality and the execution time. Often, there is a question whether to make afirst better orthe best improvement move in the LS phase [152]. Some E-VRP researchers perform one or the other, and some combine the two approaches, i.e., best move of first 100 [11] or first 50 [31].Neighborhood operators explore only the space of the immediate vicinity of the current solution, which often leads to local optima. In order to search a larger solution space, Shaw [153] proposed thelarge neighborhood search (LNS) heuristic, which is based on destroying and repairing the solution. Heuristic efficiency depends on the implemented destroy and, especially, repair operators. The main drawback of the heuristic is the repeated use of the same destroy and repair operators, which can lead to local optima. A modification of LNS by Christiaens and Vanden Berghe [154] is giving state-of-the-art results on a wide range of VRP variants. In the E-VRP, LNS is successfully implemented by Sassi et al. [47–49] and Zhang et al. [36], but most often it still needs some other guiding heuristic to escape the local optima. Macrina et al. [99] applied an LNS scheme to solve GMFVRP-PRTW. A list of used destroy and repair operators in the E-VRP is presented as a part of the adaptive LNS version in the appendix.


### 5.3.3. Route or Move Evaluation

Operator evaluation is an important part that highly affects execution time and has to be efficiently implemented. For example, in VRPTW, for most of the basic neighborhood operators, the load capacity check and time window check can be computed in constant timeO(1), if for each customer endogenous variables, i.e., remaining load capacity, latest arrival time, etc., are tracked. In the E-VRP, there is also a constraint regarding the vehicle battery capacity (battery feasibility), which also influences the time window feasibility. Therefore, time window evaluation is more complex. For example, when the customer is trying to be inserted into the route, the load and energy check can be computed in constant time O(1) for most of the basic neighborhood operators. The endogenous variables for vertices between the insertion position and next CS have to be recomputed in order to know the charging time at next CS. The complexity of such approach is O(B) where B is the number of customers between the position of the insertion and next CS. If there is no CS in the second part of the route, no recalculation is needed. By storing additional information about the recharging position and its amount and by the use of general concatenation operators presented by Vidal et al. [155], Hiermann et al. [31] reported that all of the basic neighborhood operators can be evaluated in O(1). Goeke and Schneider [30], instead of storing additional information about recharging, used a surrogate cost function whose evaluation is less demanding. In the search process, the authors store ϵ best moves in the composite neighborhood and evaluated either exact or surrogate cost, and at the end if there is any move with the surrogate cost, its exact cost is evaluated and the best move is chosen. Hiermann et al. [96] presented a different approach for H2E-FTW in which customer sequence without CSs is determined, and then in the route evaluation part time-consuming CSs insertions and charging schedule are performed. Whenever the route needs to be evaluated, the heuristic finds the routing cost by applying a labeling algorithm to determine CSs placement and a greedy algorithm to determine charging amount and running mode of PHEV. As such an evaluation is time-consuming, the authors performed the following: (i) use neighborhood restrictions to consider only promising arcs for each customer vertex; (ii) keep track of move evaluations in cache memory; and (iii) lower bounds to first evaluate the move as ICEV in O(1), and then if it improves the objective function, perform the route evaluation.


## 5.3.1. Constructive Heuristics

Constructive heuristics are often used to generate an initial solution by either serial or parallel route construction. Solutions are constructed in a greedy way, which often produces solutions of the VRP that are 10-15% far from an optimal solution [4]. In the E-VRP, constructive heuristics are modified and adapted to the BEV characteristics and feasibility checks. Some well-known general constructive heuristics used in E-VRP are presented.Themodified savings method of Clarke and Wright (MCWS) [141] with the insertion of CSs is used to generate an initial solution of several E-VRP and related problems [23, 32, 52, 69]. The algorithm starts with the creation of back-and-forth routes for each customer. Next, CSs are inserted in a least-cost manner, in either feasible or infeasible routes. Then, the routes that increase the savings the most are merged together and the procedure is repeated until all customers are served.Thesweep algorithm [142] inserts customers in the active route in a circular manner, resulting in efficient space division. The customers are sorted based on the value of a polar angle between the depot and randomly chosen point. Then, customers are inserted in the active route at the position causing minimal increase until a violation of route constraint, when usually a new route is opened [25, 56, 97].Nearest neighbor heuristic (NNH) starts from the depot, and in each iteration the customer with the least-cost increase from the previously selected customers is added to the route. The route is terminated when some/any constraints are violated and a new route is opened [76].Route-first cluster-second [143] approach constructs a giant infeasible TSP tour, usually with NNH, which is then split into several routes. The split procedure can be solved as the SPP on the acyclic graph. The procedure generates energy infeasible solutions, so the CS placement and charging schedule need to be determined to make the solution energy feasible [63, 68, 72].Thecluster-first route-second [144] approach divides the customers into clusters and then each cluster is solved as a subrouting problem. One example used in GVRP is thedensity based clustering algorithm (DBCA) proposed by Erdoğan and Miller-Hooks [23], which determines customer clusters based on their density, by ensuring that each customer has at least some predefined number of customers in the ϵ radial neighborhood. In the second step to obtain the routes, the MCWS is performed.Many researchers are adopting some form of the insertion heuristics, which have some characteristics of already reviewed heuristics. The general idea of the insertion heuristics is to iteratively insert or add customers in available routes. In each step of the algorithm, an unserved customer together with a route and its position in the route are determined to least increase the objective function. This selection can be deterministic but is often stochastic to select at random one fromk-best inserts to diversify the solution, especially in population metaheuristics. In most feasibility ensured cases when the energy constraint is violated, a CS is added to the route to make it energy feasible, and when the other constraints are violated such as vehicle load capacity or time windows, the new route is opened, and the procedure is repeated. Examples of adaptations to the E-VRP arepush forward insertion heuristics (PIH) [41, 145],constructive-k-PseudoGreedy [46], [28, 30, 31, 33, 36] modified insertion heuristics, andcharging routing heuristic (CRH) [47–49], in which first the charging scheme at the depot is determined, and then the least-cost feasible customer insertions are performed.


## 5.3.2. Improvement Heuristics

Improvement heuristics or local search (LS) explores the neighborhood of the incumbent solution, searching for a better solution. The neighborhood is explored by applying perturbation moves based on the composite neighborhood operators. The local search stops when no improving solution can be found in the neighborhood of the incumbent solution, which is then called local optima. A great variety of researchers have performed LS procedures to intensify the search, coupled together with perturbation moves to escape the local optima. Often, the perturbation moves are similar to the neighborhood operators used in LS phase. Most of the classical VRP neighborhood operators [3, 4] are used in the E-VRP, but some additional neighborhood operators have been developed. Depending on whether the operators perform on only one route or between the routes, they can be divided intointraroute andinterroute operators. Here, we present some of the most used LS neighborhood operators in the E-VRP literature:(i)
2-Opt [146] replaces two arcs with the two new ones, with a possibility of route direction reversal.(ii)
2-Opt∗ [147], unlike the 2-Opt operator, avoids the reversal of a route direction.(iii)
4-Opt [84, 148] replaces four consecutive arcs with four new ones.(iv)
Or-Opt [149] replaces three arcs with three new ones such that a sequence of three vertices is relocated.(v)
Exchange (Swap) [150] exchanges the two vertices in the solution.(vi)
FacilityInsertion [55] andStationInsertion [11] operators insert refueling or loading facility at best position in between two consecutive facility visits where fuel or load violation occurred. The similar operatorInsertRemoveIf [31] inserts CS into the part of the route where the violation occurred to make the route energy feasible and repair violations, but it also removes redundant CS visits. Schneider et al. [55] also proposed the facility-related removal and exchange operatorsFacilityReplacement andExchange.(vii)
RechargeRelocation [46] removes all visits to the CSs and then tries to insert CSs at positions leading to a better solution. Similarly, theGlobalChargingImprovement [72] operator removes all CS visits in the route, and if the route is energy feasible, it stops; otherwise, it solves FRVCP to optimize CSs positions and charging amount.(viii)
Relocate [150] removes one vertex from the solution and inserts it into another position in the solution. Felipe et al. [46] proposed a modification of the relocation operator calledReinsertion,Rbu, for solving GVRP-MTPR. Binary parameters b and u denote whether the first-improve or best-improve strategy is used and whether the values of savings are updated or not. TheRemoveTwoInsertOne [84] operator removes two vertices at random and inserts them one by one into another vehicle at best possible positions.(ix)
CrossExchange [151] removes n1 consecutive vertices from route r1 to route r2, and n2 consecutive vertices from route r2 to route r1 [84]. TheRemoveSequence [84] operator is similar as it removes consecutive vertices from only one route and inserts it into another one.(x)
Resize [31],RelocateAndResize [31], χ-OptWithModeChange [58], andVehicleSwap [96] change the vehicle type assignment or the propulsion mode in the heterogeneous E-VRP and HEV routing problems.(xi)
StationInRe [25] inserts and removes CSs from the solution in so that if an arc to the CS is already in the solution, it is removed; otherwise it is inserted.The order of the operators affects both the solution quality and the execution time. Often, there is a question whether to make afirst better orthe best improvement move in the LS phase [152]. Some E-VRP researchers perform one or the other, and some combine the two approaches, i.e., best move of first 100 [11] or first 50 [31].Neighborhood operators explore only the space of the immediate vicinity of the current solution, which often leads to local optima. In order to search a larger solution space, Shaw [153] proposed thelarge neighborhood search (LNS) heuristic, which is based on destroying and repairing the solution. Heuristic efficiency depends on the implemented destroy and, especially, repair operators. The main drawback of the heuristic is the repeated use of the same destroy and repair operators, which can lead to local optima. A modification of LNS by Christiaens and Vanden Berghe [154] is giving state-of-the-art results on a wide range of VRP variants. In the E-VRP, LNS is successfully implemented by Sassi et al. [47–49] and Zhang et al. [36], but most often it still needs some other guiding heuristic to escape the local optima. Macrina et al. [99] applied an LNS scheme to solve GMFVRP-PRTW. A list of used destroy and repair operators in the E-VRP is presented as a part of the adaptive LNS version in the appendix.


## 5.3.3. Route or Move Evaluation

Operator evaluation is an important part that highly affects execution time and has to be efficiently implemented. For example, in VRPTW, for most of the basic neighborhood operators, the load capacity check and time window check can be computed in constant timeO(1), if for each customer endogenous variables, i.e., remaining load capacity, latest arrival time, etc., are tracked. In the E-VRP, there is also a constraint regarding the vehicle battery capacity (battery feasibility), which also influences the time window feasibility. Therefore, time window evaluation is more complex. For example, when the customer is trying to be inserted into the route, the load and energy check can be computed in constant time O(1) for most of the basic neighborhood operators. The endogenous variables for vertices between the insertion position and next CS have to be recomputed in order to know the charging time at next CS. The complexity of such approach is O(B) where B is the number of customers between the position of the insertion and next CS. If there is no CS in the second part of the route, no recalculation is needed. By storing additional information about the recharging position and its amount and by the use of general concatenation operators presented by Vidal et al. [155], Hiermann et al. [31] reported that all of the basic neighborhood operators can be evaluated in O(1). Goeke and Schneider [30], instead of storing additional information about recharging, used a surrogate cost function whose evaluation is less demanding. In the search process, the authors store ϵ best moves in the composite neighborhood and evaluated either exact or surrogate cost, and at the end if there is any move with the surrogate cost, its exact cost is evaluated and the best move is chosen. Hiermann et al. [96] presented a different approach for H2E-FTW in which customer sequence without CSs is determined, and then in the route evaluation part time-consuming CSs insertions and charging schedule are performed. Whenever the route needs to be evaluated, the heuristic finds the routing cost by applying a labeling algorithm to determine CSs placement and a greedy algorithm to determine charging amount and running mode of PHEV. As such an evaluation is time-consuming, the authors performed the following: (i) use neighborhood restrictions to consider only promising arcs for each customer vertex; (ii) keep track of move evaluations in cache memory; and (iii) lower bounds to first evaluate the move as ICEV in O(1), and then if it improves the objective function, perform the route evaluation.


## 5.4. Hybrid Metaheuristic Procedures

Many researchers employ metaheuristics to continue the exploration after the first local optima occurrence. Metaheuristics can be defined as heuristics guiding other heuristics and can be divided into neighborhood-oriented metaheuristics and population metaheuristics. Most often, multiple metaheuristics and heuristics are combined together and adapted to the problem; therefore, the term hybrid is used in such occasions.


### 5.4.1. Neighborhood-Oriented Metaheuristics

Neighborhood-oriented heuristics iteratively explore the neighborhood of the incumbent solution. Here, we present the most used neighborhood heuristics in the E-VRP literature.Simulated annealing (SA) [156] controls the search by the temperature parameter, which mimics the physical cooling process of the material. Most often, it is used as a control heuristic for accepting solutions produced by other heuristics. The probability of accepting solution S2, which is worse than the current best solution S1, is computed as e(f(S1)-f(S2))/T (Boltzmann function [157]), where T is the temperature parameter. The higher the temperature value, the higher is the probability of accepting a deteriorating move to escape the local optima. The procedure starts with high temperatures at the beginning and then dynamically lowers the temperature as the search process comes to an end. In the E-VRP, most often initial temperature is selected such that a μ% worse solution than the initial solution is accepted with a probability of 0.5 [25, 28, 30, 46].Tabu search (TS) [158, 159] is a well-known metaheuristic, which uses memory to forbid the exploration of the space that has already been explored. It escapes the local optima by accepting a worse solution if it is the best in the explored space. Arcs that are deleted from the solution are stored in thetabu-list, which prohibits the reinsertion of the deleted arcs into a specific part of the solution for a designated number of iterations,tabu-tenure. In some implementation cases, when the new best solution is found, some of the moves and arcs are lifted and removed from thetabu-list,aspire criterion [25]. In the E-VRP, TS with neighborhood operators is used to intensify the search [25, 35, 48].Variable neighborhood search (VNS) [160] is a metaheuristic that systematically changes the neighborhood each time there is no improvement in LS phase. The changes are based on predefined neighborhood structures. In each step, the new neighborhood operator is used, which increases the vicinity from the incumbent solution. The selection of neighborhood operators and customers can be deterministic (variable neighborhood decent (VND)) or probabilistic. Schneider et al. [25] successfully adapted it in the E-VRPTW for diversification of the solution based on the cyclic exchange moves, resulting in 15 different neighborhood structures. The VNS in the E-VRP was also combined with exact algorithms like branch-and-price [51], evolutionary algorithms [84], deterministic LS operators [64, 72], and adaptive mechanisms (AVNS) [11, 27, 55, 69].Iterated local search (ILS) [161] is based on successively repeating LS on the incumbent solution. When LS ends up in local optima, perturbation move is applied to escape from local optima. The perturbation move can be scaled dynamically to overcome local optima. The effectiveness of the procedure highly depends on the local search procedures used. Montoya et al. [72] applied ILS together with VND for solving E-VRP-NL with perturbation performed by splitting perturbed giant TSP tour. Macrina et al. [33] used typical ILS with some of the LS operators used in the perturbation phase to solve GMFVRP-PRTW. Zhang et al. [36] used ILS and LNS in an ant colony framework for solving GVRP. The LNS is used as a perturbation mechanism while the LS is used to further improve the solution. Sassi et al. [47, 49] also used a multistart combination of ILS and LNS but without the LS phase. Sassi et al. [48] applied similar ILS paradigm, but they called ititerated tabu search (ITS), as they used TS for the intensification and LNS for the perturbation. A similar ITS procedure is applied by Doppstadt et al. [58] for the HEV-TSP.Adaptive large neighborhood search (ALNS) [162, 163] is the extension of LNS heuristic in which, throughout the search process, different procedures for destroying and repairing the solution are selected in an adaptive manner based on their past performance. The idea is to have versatile destroy and repair operators in order to escape the local optima. For each destroy and repair operator, the adaptive weight and score are increased in the following order (scores σ1≥σ3≥σ2): (i) σ1 if a new best solution is found; (ii) σ2 if the solution is better than the previous one; (iii) σ3 if the solution is worse than the previous one but it is accepted due to the acceptance criteria; and (iv) no increase if the solution is not accepted. The score of worse solution σ3 is greater than the score of better solution σ2, to reward unimproved solutions that diversify the search [28]. For each N iteration, the weights are updated based on the previous weight values and new score values. As a great number of the researchers used ALNS for solving the E-VRP and related problems (Table 2), a vast number of destroy and repair operators were applied. Some of them originated from Pisinger and Ropke [162] but also some new problem-specific ones regarding the CS have been developed [28, 30–32, 79, 88]. The description of the applied destroy and repair operators in E-VRP is presented in the appendix. Schiffer and Walther [89] applied a parallel ALNS version of Schiffer and Walther [88] to solve the RELRPTWPR. The OpenMP was used to solve multiple scenarios in parallel, where the customer destroy and repair operators, as well as the LS and DP, are performed in parallel. The change of facility configuration, adaptive evaluations, and stopping criterion are not parallelized.SA is one of the most used metaheuristics for acceptance criteria, in order to escape the local optima, but some others can also be applied. Tiwari et al. [164] proposed the Cauchy function (CA) T/((f(S1)-f(S2))2+T2) for the acceptance of the solution as it gives more opportunities to escape from local optima. The record-to-record (RR) principle [165], which adds negative value to the objective function of the best solution so far, can also be used to escape the local optima. Several papers in the E-VRP field applied such heuristics for solution acceptance [48, 49, 84].


### 5.4.2. Population Metaheuristics

Population metaheuristics are based on the natural selection to evolve a population and survival of the fittest. They have been widely applied in the VRP field: genetic algorithm [166], path relinking [167], scatter search [168], ant colony [169], bee colony [170], particle swarm optimization [171, 172], etc.Genetic algorithm (GA) consists of a pool of individuals (VRP solutions) called a population, which goes through the process of evolution. Such evolution process includes evaluation of individuals, selection of parents, crossover, mutation, and a replacement of the old population by a new one. The important aspects of the algorithm are the representation of the VRP solution, efficient management of population diversity, and fitness function to compare the individuals in the population. For further algorithm description and application, we refer to Potvin [173]. Theant colony (AC) algorithm is based on ant behavior when they search for food. As they search, they leave pheromones, which accumulate over the most used path, leading to the optimal path between the food and the nest. AC algorithm consists of three procedures: pheromone initialization, probability adjustment, and pheromone update.Particle swarm optimization (PSO) mimics the behavior of organisms such as fish schooling or bird flocking. Physical movement of the individual (particle) in the swarm is tracked and then well-balanced moves are applied to adapt the individual to the global and local exploration. Each particle has a position in the swarm and corresponding velocity in order to optimize the problem following personal best solution and global best solution [172].Barco et al. [34, 42] applieddifferential evolution to solve the E-VRP and charge scheduling. The initial individual represents the valid route assignment and charging profile. Next, mutation based on OR and XOR operators is performed on three randomly selected individuals. The crossover is a random 2-point while the selection chooses between the target individual and new individual, which is generated by the crossover. Alesiani and Maslekar [45] applied GA for solving a single BEV routing problem. The individual is represented as a two-dimensional binary array, where the rows represent the routes and the columns represent CSs. The initial population contains feasible individuals, individuals without CSs, and individuals containing battery CSPP routes. The authors select parents at random or based on the probability and apply 2-point crossover. Mutation is applied at random to remove, insert, or exchange CS, while the LS is performed to improve the solution. Shao et al. [74] applied GA without detailed description for solving real-life EVRP-CTVTT instances. The authors encoded the individual with all visited vertices including customers, depots, and CSs. The crossover operator tries to create a new individual that satisfies all the constraints. Also, elitism is applied to preserve the elite individuals. The authors reported the average computational time of three hours. A similar algorithm with the LS after the crossover and mutation is applied for solving E-VRP by Shao et al. [90] to obtain an optimal operation scheme including routes, charging plan, and driving paths. Masliakova [61] applied the GA to find the optimal placement of CSs for electric bus routes. The individual is represented as the sequence of CSs and bus stops visited, and a crossover is implemented as the swap between the subroutes. Mutation moves the CS to the nearest bus stop or creates a new bus stop with CS. For the fitness function, the author used an AC approach, which ranks the individuals in the population. Joo and Lim [78] used the AC approach to determine the energy consumption and speed in energy-efficient BEV routing. Each ant selects the next vertex with a greedy stochastic search, which combines heuristic information about energy consumption, speed, and pheromones. The authors concluded that AC strategy is suitable for routing BEVs in terms of energy efficiency. An efficient version of evolutionary algorithm and VNS (EVO-VNS) is applied by Masmoudi et al. [84] for solving the DARP-EV. Each population is divided into m groups in such a way that each group contains good and bad solutions. Further on, s solutions from each group are selected in the VNS phase. If the best solution in the current group is not improved for three iterations, the crossover between the selected solutions and the rest of solutions in the group is performed. The rest of the individuals in the population are generated by applying destroy-repair procedure with regret heuristic. Poonthalir and Nadarajan [38] applied particle swarm optimization with greedy mutation operator and time-varying acceleration coefficient (TVa-PSOGMO) to solve the F-GVRP. Each particle is encoded as the sequence of customers who are converted to feasible routes by inserting depots and refueling stations. Time-varying inertia and acceleration coefficients are employed to better manage diversification and intensification exploration in the multiobjective environment. Additionally, to prevent premature convergence, a greedy mutation operator (GMO) is performed. Hiermann et al. [96] used hybrid GA (HGA) with LS to obtain good assignment and routing decisions for H2E-FTW. The authors represented the solution as a set of vehicle routes and used a biased fitness function. During the search, the authors maintained two populations of feasible and infeasible individuals. Parents are selected using binary tournament scheme from those two populations. The recombination is performed on the parent routes combined in a giant single tour, with an OX crossover and a split algorithm to obtain a complete solution. The mutation phase is performed as LNS with destroy and repair operators similar to ones of Hiermann et al. [31]. In the end, the solution is improved by the LS procedure.


## 5.4.1. Neighborhood-Oriented Metaheuristics

Neighborhood-oriented heuristics iteratively explore the neighborhood of the incumbent solution. Here, we present the most used neighborhood heuristics in the E-VRP literature.Simulated annealing (SA) [156] controls the search by the temperature parameter, which mimics the physical cooling process of the material. Most often, it is used as a control heuristic for accepting solutions produced by other heuristics. The probability of accepting solution S2, which is worse than the current best solution S1, is computed as e(f(S1)-f(S2))/T (Boltzmann function [157]), where T is the temperature parameter. The higher the temperature value, the higher is the probability of accepting a deteriorating move to escape the local optima. The procedure starts with high temperatures at the beginning and then dynamically lowers the temperature as the search process comes to an end. In the E-VRP, most often initial temperature is selected such that a μ% worse solution than the initial solution is accepted with a probability of 0.5 [25, 28, 30, 46].Tabu search (TS) [158, 159] is a well-known metaheuristic, which uses memory to forbid the exploration of the space that has already been explored. It escapes the local optima by accepting a worse solution if it is the best in the explored space. Arcs that are deleted from the solution are stored in thetabu-list, which prohibits the reinsertion of the deleted arcs into a specific part of the solution for a designated number of iterations,tabu-tenure. In some implementation cases, when the new best solution is found, some of the moves and arcs are lifted and removed from thetabu-list,aspire criterion [25]. In the E-VRP, TS with neighborhood operators is used to intensify the search [25, 35, 48].Variable neighborhood search (VNS) [160] is a metaheuristic that systematically changes the neighborhood each time there is no improvement in LS phase. The changes are based on predefined neighborhood structures. In each step, the new neighborhood operator is used, which increases the vicinity from the incumbent solution. The selection of neighborhood operators and customers can be deterministic (variable neighborhood decent (VND)) or probabilistic. Schneider et al. [25] successfully adapted it in the E-VRPTW for diversification of the solution based on the cyclic exchange moves, resulting in 15 different neighborhood structures. The VNS in the E-VRP was also combined with exact algorithms like branch-and-price [51], evolutionary algorithms [84], deterministic LS operators [64, 72], and adaptive mechanisms (AVNS) [11, 27, 55, 69].Iterated local search (ILS) [161] is based on successively repeating LS on the incumbent solution. When LS ends up in local optima, perturbation move is applied to escape from local optima. The perturbation move can be scaled dynamically to overcome local optima. The effectiveness of the procedure highly depends on the local search procedures used. Montoya et al. [72] applied ILS together with VND for solving E-VRP-NL with perturbation performed by splitting perturbed giant TSP tour. Macrina et al. [33] used typical ILS with some of the LS operators used in the perturbation phase to solve GMFVRP-PRTW. Zhang et al. [36] used ILS and LNS in an ant colony framework for solving GVRP. The LNS is used as a perturbation mechanism while the LS is used to further improve the solution. Sassi et al. [47, 49] also used a multistart combination of ILS and LNS but without the LS phase. Sassi et al. [48] applied similar ILS paradigm, but they called ititerated tabu search (ITS), as they used TS for the intensification and LNS for the perturbation. A similar ITS procedure is applied by Doppstadt et al. [58] for the HEV-TSP.Adaptive large neighborhood search (ALNS) [162, 163] is the extension of LNS heuristic in which, throughout the search process, different procedures for destroying and repairing the solution are selected in an adaptive manner based on their past performance. The idea is to have versatile destroy and repair operators in order to escape the local optima. For each destroy and repair operator, the adaptive weight and score are increased in the following order (scores σ1≥σ3≥σ2): (i) σ1 if a new best solution is found; (ii) σ2 if the solution is better than the previous one; (iii) σ3 if the solution is worse than the previous one but it is accepted due to the acceptance criteria; and (iv) no increase if the solution is not accepted. The score of worse solution σ3 is greater than the score of better solution σ2, to reward unimproved solutions that diversify the search [28]. For each N iteration, the weights are updated based on the previous weight values and new score values. As a great number of the researchers used ALNS for solving the E-VRP and related problems (Table 2), a vast number of destroy and repair operators were applied. Some of them originated from Pisinger and Ropke [162] but also some new problem-specific ones regarding the CS have been developed [28, 30–32, 79, 88]. The description of the applied destroy and repair operators in E-VRP is presented in the appendix. Schiffer and Walther [89] applied a parallel ALNS version of Schiffer and Walther [88] to solve the RELRPTWPR. The OpenMP was used to solve multiple scenarios in parallel, where the customer destroy and repair operators, as well as the LS and DP, are performed in parallel. The change of facility configuration, adaptive evaluations, and stopping criterion are not parallelized.SA is one of the most used metaheuristics for acceptance criteria, in order to escape the local optima, but some others can also be applied. Tiwari et al. [164] proposed the Cauchy function (CA) T/((f(S1)-f(S2))2+T2) for the acceptance of the solution as it gives more opportunities to escape from local optima. The record-to-record (RR) principle [165], which adds negative value to the objective function of the best solution so far, can also be used to escape the local optima. Several papers in the E-VRP field applied such heuristics for solution acceptance [48, 49, 84].


## 5.4.2. Population Metaheuristics

Population metaheuristics are based on the natural selection to evolve a population and survival of the fittest. They have been widely applied in the VRP field: genetic algorithm [166], path relinking [167], scatter search [168], ant colony [169], bee colony [170], particle swarm optimization [171, 172], etc.Genetic algorithm (GA) consists of a pool of individuals (VRP solutions) called a population, which goes through the process of evolution. Such evolution process includes evaluation of individuals, selection of parents, crossover, mutation, and a replacement of the old population by a new one. The important aspects of the algorithm are the representation of the VRP solution, efficient management of population diversity, and fitness function to compare the individuals in the population. For further algorithm description and application, we refer to Potvin [173]. Theant colony (AC) algorithm is based on ant behavior when they search for food. As they search, they leave pheromones, which accumulate over the most used path, leading to the optimal path between the food and the nest. AC algorithm consists of three procedures: pheromone initialization, probability adjustment, and pheromone update.Particle swarm optimization (PSO) mimics the behavior of organisms such as fish schooling or bird flocking. Physical movement of the individual (particle) in the swarm is tracked and then well-balanced moves are applied to adapt the individual to the global and local exploration. Each particle has a position in the swarm and corresponding velocity in order to optimize the problem following personal best solution and global best solution [172].Barco et al. [34, 42] applieddifferential evolution to solve the E-VRP and charge scheduling. The initial individual represents the valid route assignment and charging profile. Next, mutation based on OR and XOR operators is performed on three randomly selected individuals. The crossover is a random 2-point while the selection chooses between the target individual and new individual, which is generated by the crossover. Alesiani and Maslekar [45] applied GA for solving a single BEV routing problem. The individual is represented as a two-dimensional binary array, where the rows represent the routes and the columns represent CSs. The initial population contains feasible individuals, individuals without CSs, and individuals containing battery CSPP routes. The authors select parents at random or based on the probability and apply 2-point crossover. Mutation is applied at random to remove, insert, or exchange CS, while the LS is performed to improve the solution. Shao et al. [74] applied GA without detailed description for solving real-life EVRP-CTVTT instances. The authors encoded the individual with all visited vertices including customers, depots, and CSs. The crossover operator tries to create a new individual that satisfies all the constraints. Also, elitism is applied to preserve the elite individuals. The authors reported the average computational time of three hours. A similar algorithm with the LS after the crossover and mutation is applied for solving E-VRP by Shao et al. [90] to obtain an optimal operation scheme including routes, charging plan, and driving paths. Masliakova [61] applied the GA to find the optimal placement of CSs for electric bus routes. The individual is represented as the sequence of CSs and bus stops visited, and a crossover is implemented as the swap between the subroutes. Mutation moves the CS to the nearest bus stop or creates a new bus stop with CS. For the fitness function, the author used an AC approach, which ranks the individuals in the population. Joo and Lim [78] used the AC approach to determine the energy consumption and speed in energy-efficient BEV routing. Each ant selects the next vertex with a greedy stochastic search, which combines heuristic information about energy consumption, speed, and pheromones. The authors concluded that AC strategy is suitable for routing BEVs in terms of energy efficiency. An efficient version of evolutionary algorithm and VNS (EVO-VNS) is applied by Masmoudi et al. [84] for solving the DARP-EV. Each population is divided into m groups in such a way that each group contains good and bad solutions. Further on, s solutions from each group are selected in the VNS phase. If the best solution in the current group is not improved for three iterations, the crossover between the selected solutions and the rest of solutions in the group is performed. The rest of the individuals in the population are generated by applying destroy-repair procedure with regret heuristic. Poonthalir and Nadarajan [38] applied particle swarm optimization with greedy mutation operator and time-varying acceleration coefficient (TVa-PSOGMO) to solve the F-GVRP. Each particle is encoded as the sequence of customers who are converted to feasible routes by inserting depots and refueling stations. Time-varying inertia and acceleration coefficients are employed to better manage diversification and intensification exploration in the multiobjective environment. Additionally, to prevent premature convergence, a greedy mutation operator (GMO) is performed. Hiermann et al. [96] used hybrid GA (HGA) with LS to obtain good assignment and routing decisions for H2E-FTW. The authors represented the solution as a set of vehicle routes and used a biased fitness function. During the search, the authors maintained two populations of feasible and infeasible individuals. Parents are selected using binary tournament scheme from those two populations. The recombination is performed on the parent routes combined in a giant single tour, with an OX crossover and a split algorithm to obtain a complete solution. The mutation phase is performed as LNS with destroy and repair operators similar to ones of Hiermann et al. [31]. In the end, the solution is improved by the LS procedure.


## 5.5. Summary of the Applied Procedures

The overview of the applied procedures in the E-VRP field is presented in Table2, where 78 papers are presented. Searching in the space of infeasible solutions is applied in 17 papers, while the LS for the intensification is applied in 25 papers. Most of the researchers do not use any special criteria for accepting the solution, but 17 papers use SA as acceptance criteria, three RR and one CA. For the initial solution, mostly some form of the universal insertion heuristic is applied (18). Then follow MCWS (11), NNH (7), random (5), split procedure (4), and sweep (3). Regarding the applied heuristics for solving the problem, the order is the following: ALNS (14), GA (10), LNS (7), TS and ITS (6), ILS (6), AC (3), VNS (3), AVNS (2), VND (2), PSO (1), and differential evolution (1). From all of the exact procedures, the DP is the most applied one (12) and then there is the branch-bound-cut-and-price (5). Forty-three papers use commercially available software like CPLEX (32), Gurobi (7), MATLAB (2), and XPRESS (2) for solving either the whole problem or some subproblem during the search, mostly charging decisions.We can note that a vast number of researchers use hybrid heuristics procedures with only a few focusing on the exact procedures [31, 57]. It can be noted that most of the procedures that are giving state-of-the-art results are using exact procedures at some point during the search, usually for (i) optimal determination of charging visits and charging amount [79, 88, 96] and (ii) set partitioning over the set of routes at the end of the search [65, 72, 95, 96].The ALNS is one of the most applied procedures and it is giving the BKSs on various problem variants. Among them, we can point out the following: ALNS of Hiermann et al. [31] (HPH) and ALNS of Goeke and Schneider [30] (GS) for heterogeneous fleet, ALNS of Keskin and Çatay [28] (KC1) and ANLS of Keskin and Çatay [79] (KC2) for partial recharges and different charging technologies, and ALNS of Schiffer and Walther [88] (SW) originally designed for the LRPIF (ELRPTWPR) which is also giving BKSs on the E-VRPTW and E-VRPTWPR. Schiffer and Walther [88] compared their algorithm on the E-VRPTW instances to the HPH, GS, KC1, and VNS-TS of Schneider et al. [25] (SSG). The authors found nine BKSs with an average gap of 0.27% to the previous BKSs and average running time of 3.77 minutes. On the E-VRPTWPR, compared to the KC1, the authors found 41 new BKSs. The proposed ALNS of SW is efficiently parallelized to solve the RELRPTWPR [89].Most of the population-based metaheuristics are not efficiently implemented so they are not giving high-quality solutions. But, the HGA of Hiermann et al. [96] (HHPV) gives some of the BKSs on the E-VRPTW, E-VRPTWPR, and H2E-FTW. The authors compared their results on the E-VRPTW and E-VRPTWPR to the SSG, HPH, GS, and KC1. On the E-VRPTW instances, the HHPV found 11 BKSs with 1.36% increase in the vehicle number, 0.46% distance gap from the previous BKSs, and an average running time of 11.59 minutes. On the E-VRPTWPR, the authors found 28 new BKSs. By comparing the results of HHPV and SW on E-VRPTW, we can conclude that in 20 instances SW outperformed HHPV, while in the opposite case the HHPV was better in 17 instances.


## 6. Conclusion and Future Research Directions

With the increase in popularity of electric-powered vehicles, many logistics companies have started to integrate electric vehicles into their delivery fleets. Routing a fleet of electric vehicles for delivering goods was formulated as the electric vehicle routing problem. Besides the often used load (cargo) capacity and time window constraints, E-VRP routing models have to account for the limited driving range of BEVs, which directly corresponds to the more frequent recharging needs at CS.In this paper, a literature review on recent developments in the E-VRP field is presented. We consider BEVs to be vehicles powered only from batteries mounted inside the vehicle. A general overview of the BEV’s characteristics for goods delivery includes the driving range, battery capacity, application, and case studies. As energy consumption estimation is an important part of BEV routing, we summarized the recent research on theoretical and data-driven energy consumption models.Due to the BEVs specific characteristics, new E-VRP variants have emerged: a mixed fleet of electric and conventional vehicles, partial recharging, simultaneous CS siting and BEV routing, nonlinear charging, different charging technologies, battery swap technology, hybrid vehicles, green routing, etc. The development of efficient heuristics was necessary to find optimal or near-optimal solutions to the new routing problems. We reviewed the state-of-the-art exact, heuristic, and hybrid procedures applied for solving various E-VRP variants. The adaptive large neighborhood search [88] and hybrid genetic algorithm [96] produced high-quality solutions on various E-VRP variants. Most of the state-of-the-art procedures use exact algorithms during the search to determine the optimal placement of CSs and the charging amount.From the literature review, it can be noted that the electric vehicle routing research community has grown rapidly in the last few years and many problem variants have already been explored. Nevertheless, we can highlight possible future directions as follows.By our observation, there is a lack of papers regarding case studies and application cases, where actual E-VRP models could be evaluated, and some meaningful insights could be drawn. Several energy consumption models were reviewed, but only few are predicting realistic energy consumption at the road segment level in the road network. Only recently, have researchers started to integrate a nonlinear charging process, CS location problem, and hybrid electric vehicles into the E-VRP models. A few papers addressed the problem of CS capacity, as a limited number of BEVs can charge simultaneously at a CS. Several variations were observed: waiting times [68, 80, 81], battery reservations [44], and adaptive routing with uncertainties in CS availability [75, 81]. Furthermore, the following real-life characteristics have not been sufficiently studied in the E-VRP field: dynamic traffic conditions [74], uncertainties in demand, travel time, time windows, service time and charging time [75, 81, 89], compatibility of BEVs with chargers in CSs [48, 49], and recharging at public or private CSs [81].For future research regarding the solution procedures, we can highlight the following. Although there have been a couple of exact procedures developed for GVRP, only a few exact procedures have been proposed for the E-VRP and its extensions. A great number of researchers applied a population metaheuristic to solve the problem, but only a handful of them produced high-quality solutions in reasonable computation time. To improve the computation time, parallelized procedures could be used [18, 89]. Furthermore, more general procedures could be investigated as most of the state-of-the-art procedures are using complex problem-specific heuristics.



---

*Source: 5075671-2019-05-09.xml*