key: cord-0270406-2cto9hqi authors: Lin, Bo; Ghaddar, Bissan; Nathwani, Jatin title: Electric Vehicle Routing and Charging/Discharging under Time-Variant Electricity Prices date: 2020-12-17 journal: nan DOI: 10.1016/j.trc.2021.103285 sha: 4522774ca06a7f4f930ecc28e3388b7310e16c0e doc_id: 270406 cord_uid: 2cto9hqi The integration of electric vehicles (EVs) with the energy grid has become an important area of research due to the increasing EV penetration in today's transportation systems. Under appropriate management of EV charging and discharging, the grid can currently satisfy the energy requirements of a considerable number of EVs. Furthermore, EVs can help enhance the reliability and stability of the energy grid through ancillary services such as energy storage. This paper proposes the EV routing problem with time windows under time-variant electricity prices (EVRPTW-TP) which optimizes the routing of an EV fleet that are delivering products to customers, jointly with the scheduling of the charging and discharging of the EVs from/to the grid. The proposed model is a multiperiod vehicle routing problem where EVs can stop at charging stations to either recharge their batteries or inject stored energy to the grid. Given the energy costs that vary based on time-of-use, the charging and discharging schedules of the EVs are optimized to benefit from the capability of storing energy by shifting energy demands from peak hours to off-peak hours when the energy price is lower. The vehicles can recover the energy costs and potentially realize profits by injecting energy back to the grid at high price periods. EVRPTW-TP is formulated as an optimization problem. A Lagrangian relaxation approach and a hybrid variable neighborhood search/tabu search heuristic are proposed to obtain high quality lower bounds and feasible solutions, respectively. Numerical experiments on instances from the literature are provided. The proposed heuristic is also evaluated on a case study of an EV fleet providing grocery delivery at the region of Kitchener-Waterloo in Ontario, Canada. Insights on the impacts of energy pricing, service time slots, range reduction in winter as well as fleet size are presented. Over the recent years, sustainability has become a paramount global concern. In the transportation sector, public and private institutions are attempting to increase the penetration of electric vehicles 1 arXiv:2012.09357v4 [math.OC] 10 Jun 2021 (EVs) due to their ability to mitigate greenhouse gas emission and their direct impact on reducing particulate matter pollution (Boulanger et al. 2011 , Waraich et al. 2013 . Companies are also investigating the use of new green technologies due the potential brand benefits given the growing demands for green products (Kleindorfer et al. 2005 , Dekker et al. 2012 . UPS, FedEx, and Walmart are among the leading companies that have deployed EV fleets in their operations (Winston 2018) . As a consequence, the past decade has seen a rapid expansion of EV adoption (Hertzke et al. 2019 ). According to Statistics Canada (2019), the number of new EV registrations has increased from 25,163 in 2014 to 69,010 in 2018. As projected by Canada Energy Board (National Energy Board 2018), EVs will account for over 60% of the new motor vehicle registrations in Canada by 2040. Not surprisingly, the growing penetration of EVs has a significant impact on the power system. Previous studies have shown that, without proper management, the EVs will represent a sizable fraction of the total demand for energy (Triviño-Cabrera et al. 2019 , Dyke et al. 2010 . As a result, the gap between peak and off-peak demands will increase, and the ramp requirements will affect the stability and reliability of the power network (Villar et al. 2012) . With optimized scheduling of the charging, the existing power system infrastructure can accommodate the energy requirements of a considerable number of EVs (Razeghi and Samuelsen 2016 , Kintner-Meyer et al. 2007 , Letendre et al. 2008 . In doing so, the need for installing new capacity, which is expensive, time-consuming, and harmful to the environment, can be minimized (Villar et al. 2012) . Furthermore, EVs can be incorporated into the power grid as a reliable and cost-effective distributed power storage (Kempton and Letendre 1997) . By optimizing the charging and discharging to the grid, an EV fleet connected to the energy grid can assist to level out peaks in the overall electricity consumption and support the utilization of intermittent renewable energy. Although the vehicle-to-grid (V2G) connectivity represents an enticing idea, it nonetheless remains in the pilot stages of development (Sovacool et al. 2018) , and is mainly focused around a centralized architecture where a controller manages the ancillary services of a large group of EVs that are charging and discharging energy from/to the grid (Guille and Gross 2009, Sortomme and El-Sharkawi 2010) . Commercial EV fleet owners such as logistics and e-commerce companies are naturally strong candidates for realizing the benefits of EV integration in the energy grid as they aggregate a large number of EVs. This paper considers a delivery service system that operates a fleet of EVs that are primarily used to deliver products to customers. The EVs can be charged or discharged at the home depot or at charging stations in the network. When charging, a cost is paid according to the energy price at the time of use. If energy is discharged to the grid, a profit is creditted to the EV. We thus propose the EV routing problem with time windows under time-variant electricity prices (EVRPTW-TP) which optimizes the monetary cost of the EV fleet operation while allowing the charging and discharging of EVs at time-of-use electricity prices. EVRPTW-TP extends the EV routing problem with time windows (Schneider et al. 2014 ) by incorporating additional operational constraints, allowing the partial charging and discharging of the EVs, and accounting for the time-varying electricity prices. In order to solve EVRPTW-TP, a hybrid Variable Neighborhood Search/Tabu Search (VNS/TS) heuristic that can generate high quality feasible solutions efficiently is developed. A Lagrangian Relaxation that is solved by a cutting plane approach is also proposed to obtain lower bounds. The results are evaluated using a variation of the widely-used vehicle routing instances of Solomon (1987) . Finally, the model is evaluated using a case study of an EV fleet performing grocery delivery in the Kitchener-Waterloo region in Ontario, Canada. Managerial insights are drawn with respect to electricity pricing, time slots design, winter range reduction, and fleet size. This paper is the first to investigate the joint optimization of routing and charge/discharge scheduling of multiple EVs under time-variant electricity prices. The proposed model provides operational decisions to support commercial EV fleet operators in order to lower the overall energy costs. The proposed model also offers important implications for policy makers and can assist the power system regulators to better understand the impact of EV fleets on energy markets, and to predict and estimate the market reaction to energy price adjustments. The managerial insights extracted can help policy makers in creating more efficient energy pricing schemes to maximize the environmental and economic benefits from the widespread adoption of commercial EV fleets. The rest of this paper is organized as follows. Following this introductory section, Section 2 reviews the related literature. The proposed problem is formulated in Section 3. The Lagrangian relaxation is presented in Section 4 and the proposed VNS/TS heuristic is then discussed in Section 5. Computational results and the case study are presented in Sections 6 and 7, respectively. Finally, Section 8 concludes and highlights future research directions. The vehicle routing problem (VRP) was first proposed by Dantzig and Ramser (1959) as a generalization of the well-known travelling salesman problem. In general, given a set of geographically scattered customers each associated with a demand, the VRP seeks to assign customers to vehicles in such a manner that the demand of each customer is satisfied while the total distance travelled by the fleet of vehicles to serve all the customers is minimized. Since the introduction of the classical VRP, numerous variants were developed and investigated to account for realistic constraints and objectives. One of the most common variants is the VRP with time windows where visits to individual customers are restricted to fixed time intervals (Russell 1977) . Another variation is the green vehicle routing problem introduced by Erdogan and Miller-Hooks (2012) which particularly models alternative fuel vehicles and accounts for the opportunity to extend a vehicle's distance limitation by visiting an en-route station facility to replenish. Schneider et al. (2014) tailors this framework specifically to EVs and proposes the EV routing problem with time windows (EVRPTW). Instead of using a constant replenishment time as in Erdogan and Miller-Hooks (2012) , EVRPTW assumes a linear energy charging time that is associated with the battery level of the electric vehicle upon arrival to a station. Following Erdogan and Miller-Hooks (2012) and Schneider et al. (2014) , various studies have investigated the optimization of EV routing and charging/discharging operations which are summarized in Table 1 . In the uni-directional V2G contexts, Felipe et al. (2014) and Keskin and Ç atay (2016) consider partial charging strategies with multiple types of chargers, each with a different charging speed and static unit cost. Yang et al. (2015) optimize over the monetary cost of an EV providing pickup and delivery services under time-variant charging prices. The problem is extended by Barco et al. (2017) to a multi-vehicle case incorporating the battery degradation cost. Also under time-variant electricity prices, Yu and Lam (2018) models a fleet of autonomous EVs that provide customer delivery and renewable energy storage to the energy grid. A quadraticly-constrained mixed integer program is formulated and three objective functions are proposed to either minimize the total driving distance, to maximize the amount of energy charged from renewable sources, or to minimize the amount of time until the vehicles reach the final destinations. In the bi-directional V2G contexts where EVs are allowed to inject energy back to the grid, Tang et al. (2017) consider a set of EVs travelling from their origins to the corresponding destinations without en-route customers. EVs can detour for en-route charging and discharging to two types of stations, one providing renewable energy at a low price while the other is a normal station with higher charging cost and provide discharging reward. The energy prices do not vary across time, the objective is to minimize the overall monetary cost of the EV fleet. Triviño-Cabrera et al. (2019) extend the problem by incorporating intermediate stops for EVs and considering time-variant electricity prices as well as battery degradation cost. However, the model is not able to coordinate the schedules for different EVs and does not have a customer delivery component. Abdulaal et al. (2017) study a similar problem for an EV fleet considering EV congestion at charging/discharging stations. The customer assignments to different EVs are assumed to be given, the problem thus degenerate to a single EV case. Moreover, the routing and charging decisions are made sequentially, which are likely to be sub-optimal. To the best of our knowledge, no previous research has been conducted to jointly optimize the routing and charging/discharging operations of multiple EVs under time-variant electricity prices. The model proposed in this paper fills the gap between considering energy networks and transportation independently. The proposed model extends the work of Schneider et al. (2014) to include energy discharging to the grid in a multi-period framework that also accounts for the changing energy prices. Thus the proposed model can be seen as a bi-directional V2G system where the fleet of EVs whose primary purpose is to deliver products to customers can also be used to store and redistribute energy from/to the energy grid. The EV routing problem with time windows under time-variant electricity prices is presented next. Felipe et al. (2014) partial Yang et al. (2015) partial Keskin and Ç atay (2016) partial Desaulniers et al. (2016) partial Tang et al. (2017) partial Abdulaal et al. (2017) partial Barco et al. (2017) partial Yu and Lam (2018) partial Triviño-Cabrera et al. (2019) partial EVRPTW-TP partial To formulate EVRPTW-TP, a complete directed graph G(V c,s,od , E) is considered where V c,s,od denotes the set of all nodes and E is the set of all edges. The nodes are partitioned into three distinct categories: customer nodes, station nodes, and a depot. Let V c = {1, 2, . . . , N } denote the set of N customers where each node i is associated with a demand q i , a service time s i , and a time window [e i , l i ] during which an EV should arrive to node i. Let V s = {N + 1, . . . , N + S} be the set of S en-route stations where an EV can charge or discharge energy. The depot is denoted by two nodes 0 and N + S + 1, i.e. V od = {0, N + S + 1}, where node 0 is the start of the vehicle route and N + S + 1 is the end of the route. Each station or depot node i ∈ V s ∪ V od has a time window [e i , l i ] from the start to the end of the planning horizon. Each edge (i, j) is associated with a distance a ij and time t ij that denotes the travel distance and time between nodes i and j, respectively. An EV's instantaneous power depends on its mass, acceleration, as well as aerodynamic, rolling, and grade resistances (Wu et al. 2015) . The change in EV mass is negligible for certain applications, e.g. small-package delivery, while the resistances could be captured by the driving speed and acceleration assuming related parameters that are given in (Schneider et al. 2014) . Additionally, the power consumption does not directly translate to the battery's state of charge due to the nonlinear battery efficiency and other real-world characteristics such as temperature (Hannan et al. 2017 ). Incorporating these realistic features could more accurately describe the energy consumption and recharging/discharging processes, yet they will introduce computational burdens. For the sake of keeping the model tractable, following Schneider et al. (2014) and Desaulniers et al. (2016) , we assume a constant energy consumption rate g and a constant charging speed 1 α . The energy consumption for edge (i, j) is thus given by c ij = ga ij , the time required to recharge the energy consumed by traveling along edge (i, j) is given by f ij = αc ij . We assume a commercial EV fleet consisting of K homogeneous EVs, each with a load capacity Q and a battery capacity C. The time required to fully recharge the battery from 0 is defined as B = αC. At the beginning of the planning horizon, all the EVs are at the depot (node 0) with a full battery. All EVs should return to the depot (node N +S +1) before the end of the planning horizon. During each visit to a station or to the depot, an EV can either pay to charge its battery or make profits by injecting energy back to the grid from its battery. The charging cost and discharging reward vary according to the time period (time-of-use energy pricing). We assume that EVs are allowed to perform either charging or discharging during their station visits, but are not allowed to perform both during the same time period. The planning horizon is formed of |T | consecutive discrete periods, each of length δ. Each time period t refers to a time interval [δ(t − 1), δt) and is associated with a buying (from the grid) energy price p t b and a selling (to the grid) energy price p t s . We assume that if an EV is to charge or discharge during a time period, then it has to do so for the full time period. Due to the constant linear charging/discharging assumption, a fixed amount of energy δ B C would be charged/discharged during each period, we therefore can calculate a charging cost P t re = δ B Cp t b and a discharging reward P t dis = δ B Cp t s for each period. EVs are recharged back to full battery capacity during the night using at a lower energy price p night . For the ease of presentation, we define P night = Table 2 . The problem formulation that is proposed next jointly optimizes the routes and charging/discharging schedule of the K vehicles by minimizing the net electricity cost given that all customer demands are satisfied. To formulate EVRPTW-TP, the following binary variables are introduced: Amount of time required to fully charge the EV battery from empty α The reciprocal of the constant charging speed g Energy consumption rate with respect to distance traveled a ij Travel distance of edge (i, j) t ij Travel time of edge (i, j) c ij Energy consumption along edge (i, j) f ij The amount of time required to charge the energy consumed along edge (i, j) e i Earliest service start time at node i l i Latest service start time at node i s i Required service time at node i q i The demand at node i p night Unit electricity buying (from the grid) price at night (¢/kWh) p t b Unit electricity buying (from the grid) price during period t (¢/kWh) p t s Unit electricity selling (to the grid) price during period t (¢/kWh) P night The cost of charging per unit of time at night P t re Cost of charging during period t P t dis Reward for discharging during period t b ik : remaining energy (in terms of charging time) in vehicle k upon arrival to node i, u ik : remaining cargo in vehicle k upon arrival to node i. The decision variables r itk and d itk are only associated with the station and depot nodes, while the remaining variables are associated with all the nodes. EVRPTW-TP is formulated as i∈Vc,s,o j∈Vc,s,o The objective function (1) minimizes the net cost, i.e. the total cost of charging the batteries minus the total reward earned from discharging the batteries. The first part of the objective function corresponds to the net cost during the planning horizon, while the second part is the cost of fully recharging back all EVs at night. Constraints (2) ensure that every customer is served by exactly one EV. Constraints (3) force all the EVs to return to the depot by the end of the planning horizon. Constraints (4) guarantee that no route ends at a customer or a station node. Constraints (5) ensure time feasibility of the edges leaving the customer and the depot nodes, while constraints (6) deal with the edges originating from charging stations. Constraints (7) ensure that the time windows of all the nodes are not violated. Constraints (8) indicate that every EV is fully charged before leaving the depot. Constraints (9)-(10) track the battery capacity along the route. Constraints (11) indicate that an EV battery cannot be recharged to a level that exceeds its capacity, while Constraints (12) state that an EV battery cannot be discharged to a level below 0. Constraints (13) indicate that an EV is allowed to discharge or recharge its battery at a station or the depot nodes but is not allowed to discharge and recharge during the same time period. Constraints (14) ensure that an EV cannot start charging/discharging at a station/depot before arrival and before the start of the time period. Constraints (15) As shown in Section 6.3, solving EVRPTW-TP to optimality is computationally very challenging. Lagrangian relaxation is a well known algorithm that has been used to address many complex optimization problems. Particularly, in the context of vehicle routing, Lagrangian relaxation has been used to address several variants of VRP (Fisher et al. 1997 , Kallehauge et al. 2006 ). For the EVRPTW-TP, all the constraints other than Constraints (2) are associated with a particular vehicle k. Given this special structure which is common in vehicle routing problems, Constraints (2) are relaxed and the violation is penalized in the objective function using the Lagrangian multipliers λ i . The resulting relaxed problem is Since the EV fleet is homogeneous, problem Z LR (λ) decomposes into K identical sub-problems The value of the Lagrangian bound Z LR (λ) is and the best Lagrangian bound is given by max λ∈R |Vc| Z LR (λ). Given H, the set of feasible solutions of the Lagrangian subproblem, the best Lagrangian bound can be found by solving the following problem where r h it , d h it , and x h ij describe the charging, discharging, and routing schedule associated with solution h ∈ H, respectively, which is equivalent to the Lagrangian master problem Since the set H is not known beforehand, we implemented the stablized cutting-plane approach developed by Kallehauge et al. (2006) to obtain Z M P iteratively starting with an empty set H. The determination of the set of cutting planes require solution of the subproblem. Given fixed values of the Lagrangian multipliers λ i , the Lagrangian subproblem is solved to obtain Z SP (λ) and a new feasible solution h ∈ H. The resulting solution generates a cut that is added to the master problem. The relaxed master problem is solved to obtain new values for the Lagrangian multipliers. The solution of the relaxed master problem provides an upper bound on the optimal Lagrangian bound while the optimal solution of Z LR (λ) provides a lower bound. The algorithm iterates until the gap between the upper and lower bounds is sufficiently small. Alternatively, the Lagrangian bound can also be obtained by solving the Dantzig-Wolfe reformulation of Z M P using column generation. Both methods require solving Z SP (λ) repetitively. The Lagrangian sub-problem has similar structure as the column generation sub-problem introduced by Desaulniers et al. (2016) for EVRPTW. However, it is nontrivial to apply the labeling algorithm proposed by Desaulniers et al. (2016) to solve Z SP (λ) because of the difference between their objective functions. Although the EV's energy consumption along the route is linear with respect to the total distance, the cost of a unit distance in Z SP (λ) depends on the time of charging which introduces extra complexity to the problem. In addition, the decision about discharging is largely independent of the travelling distance and hence we can not formulate the Lagrangian sub-problem as a elementary shortest path problem with resource constraints. For these reasons, we use a standard mixed integer programming (MIP) solver to solve Z SP (λ). The Lagrangian relaxation presented in Section 4 provides lower bounds on the optimal solution of problem (1)-(20). To obtain upper bounds and very importantly to implement in practice, good quality feasible solutions are needed relatively quickly. Following the framework presented in Schneider et al. (2014) for EVRPTW, this section presents a hybrid variable neighborhood search and tabu search (VNS/TS) meta-heuristic with an annealing mechanism to solve EVRPTW-TP. Hybrid VNS/TS meta-heuristics have been previously applied successfully for routing problems (Melechovskỳ et al. 2005 , Tarantilis et al. 2008 . The overall framework of the VNS/TS is shown in if counter ≥ η early then 12: Break 13: else if Accept(S, S ) then S ← S In the initialization step, routes are constructed such that each customer with a positive demand is visited exactly once. For that, the well known sweep heuristic is used to obtain the initial feasible solution without charging and discharging operations (Cordeau et al. 2001) . Customers are first sorted according to an increasing order of the geometric angle using as a reference a randomly selected point. Then, starting from the customer with the smallest angle, customers are inserted to the active route at the position resulting in the minimal increase of the travel distance of the route. Once the battery or cargo constraints of the active route are violated, a new route is initiated until the number of routes used so far is equal to the EV fleet size. Then, all the remaining customers are inserted into the last route. The VNS/TS meta-heuristic considers infeasible solutions during the search. Similar to Schneider et al. (2014), a cost function is used to evaluate the quality of a solution S. The generalized cost function f gen (S) is given by where f elec (S) is the net cost of electricity (charging cost minus discharging reward); Φ tw (S), In order to evaluate the electricity cost and the violations, let r i denote the i th node along route R of length n. The following metrics are then defined for every node along a given route T E i is the earliest service start time at node r i without violating any time window constraints before it. Similarly, T L i is the latest departure time from r i that will not result in any time window violations after it. T F S i is the forward cumulative slack time, i.e. the difference between earliest arrival time and the earliest service start time, from the last station/depot to r i . T BS i is the backward cumulative slack time, i.e. the difference between the latest departure time and the latest service time, from r i to the preceding station/depot. F i is the amount of time required to recharge the energy consumption from the last station/depot to r i . The cargo capacity violation Φ cargo (R) of a route R is calculated as Time Window Violation We calculate the time window violation for route R as the sum of violations at each node along R. For the violation at node r i that follows r i−1 with time window violation, we assume the arrival time at r i−1 is its latest service start time. This assumption is similar to that of Schneider et al. (2014) to prevent time window violations from propagating along the route and avoid penalizing a good customer sequence only because it follows a node with time window violation. The time window violation Φ tw (R) for a route R is Battery Capacity Violation The battery violations Φ batt (R) for a route R is To calculate the electricity cost of a given route, a feasible charging/discharging schedule for each route is computed. Due to the complexity of the problem, the following assumptions are made to limit the number of potential feasible schedules: 1. The maximum number of station visits is limited to two. In practice, it is likely that this assumption is reasonable as the EV fleet needs to primarily service customers, and charging activities are to make sure that the vehicle has enough energy to complete the route and/or to offer the ancillary service of discharging for additional gain. Furthermore, frequent battery charging and discharging have a negative impact on the battery health and as such in practice it is expected that charging and discharging activities are limited to a few station visits. 2. If a vehicle visits a station to recharge, then the vehicle recharges the battery just enough to be able to complete the route to reach the depot. This assumption is applicable to the two-station visit case and excludes the possibility of over-charging at one station to discharge at another station later. 3. A vehicle is allowed to perform either charging or discharging at stations, but is not allowed to do both during a single station visit. This assumption reduces the solution space defined by the MIP formulation. Given a depot or a station node which is the i th node on route R denoted by r i ∈ R, the time window during which a vehicle can charge or discharge is T E r i , T L r i . Given that the planning horizon is discretized, then the set of time periods during which the vehicle can charge/discharge is given , namely the connected periods for r i . Given the i th and the j th node (i < j) of route R where a vehicle can charge/discharge (station or depot node) and T i and T j are the connected periods corresponding to i and j respectively, if a vehicle charges/discharges at time period t ∈ T i , then this might make it impossible for the vehicle to arrive at node j before the start of time period t . Thus M t i is defined as the mutually exclusive set for t ∈ T i which is the set of all the time periods t ∈ T j where the vehicle cannot reach node j if it is charging/discharging during time period t. The mutually exclusive set M t i is given by Given that at most two stations are allowed in each route, three different cases are analyzed to find a feasible charging/discharging schedule corresponding to the three special cases of zero stations, one station, and two stations routes, respectively. Zero stations If there are no stations along the vehicle's route and given that the EV is fully charged at the starting depot, then the only potential activities are to discharge the EV at the starting and/or ending depot. The maximum amount of time periods during which the EV can discharge while ensuring that there is enough energy available to cover the full route is given by . Given the connected periods T 1 and T n at the two depot nodes (the first and the last node along R) and the set of mutually exclusive time periods M t 1 for each time period t ∈ T 1 , the following mixed integer program maximizes the revenue from energy discharge s.t. where d m j is a binary decision variable indicating if the EV discharges at node j during time period m. We note that by enumerating the end period where the EV discharges at r 1 , the MIP simplifies to |T 1 | 0-1 knapsack problems with unit item weights that can be solved very efficiently. Given the optimal solution d * , the net cost of electricity for route R is given by One station Given that an EV is visiting a station (k th node along route R), then one possibility is that the EV has to recharge the battery at the station node to make sure there is enough energy to reach the end-of-route depot (n th node along route R). This case occurs if ∆ = F k + F n − B > 0. Otherwise, the EV can either discharge energy at the station or recharge in order to discharge later at the depot. Given the sets of connected time periods T 1 , T k , and T n , for the 1 st , k th , and last node along route R, respectively, and the mutually exclusive sets M t 1 for each time period t ∈ T 1 , and M t k for each time period t ∈ T k , the following mixed integer program maximizes the revenue from energy charging and discharging where d m j and r m j are binary decision variables indicating if the EV discharges and charges at node j during time period m, respectively. Constraint (37) specifies the amount of charging and discharging given different values of ∆. Constraint (38) ensures the EV battery capacity is not violated after charging/discharging at the en-route station. Constraint (39) guarantee the EV has enough energy to complete the trip. Constraints (40)-(41) force that the EV can perform either charging or discharging at the station. Constraints (42)-(43) describe the mutual exclusiveness among periods at r 1 , r k and r n . Constraint (44)-(45) specify the domains of the decision variables. Given the optimal solution d * and r * , the net cost of electricity for route R is given by Two stations Given an EV route that includes visits to two stations k 1 and k 2 with k 1 < k 2 . If ∆ = F k 1 + F k 2 + F n − B > 0 then the EV needs to recharge en-route in order to complete the trip to the final depot. In order to formulate the problem to minimize the cost of energy recharging, we first define the sets of connected time periods T k 1 and T k 2 corresponding to stations k 1 and k 2 respectively. Furthermore, let M t k 1 be the mutually exclusive sets for each time period t ∈ T k 1 . Given the binary variables r t k j , the problem that minimizes the cost of energy recharging is given by The objective function (47) minimizes the en-route electricity recharging cost. Constraints (48) make sure the EV has enough energy to complete the trip. Constraints (49) ensure that the amount of recharged energy at the first station does not exceed the available battery capacity and will allow the EV to reach the next station. Similarly, constraints (50) enforce the capacity limits at the second station. Finally, constraints (51) enforce the mutual exclusive conditions for the time periods. Given r * , the optimal solution of problem (47)-(51), the route's cost is then given by If ∆ < 0, then the EV can discharge energy at the depots (nodes 1 and n) or the stations (nodes k 1 and k 2 ). To formulate the problem that maximizes the value of the discharged energy, we introduce binary decision variable d t k that takes value 1 if the EV discharges at node k during period t and takes value 0 otherwise. Given the connected time periods T 1 , T k 1 , T k 2 , and T n corresponding to the starting depot, stations k 1 and k 2 , and the ending depot, respectively, and M t 1 , M t k 1 , and M t k 2 the mutually exclusive sets for each time period t ∈ T 1 , t ∈ T k 1 , and t ∈ T k 2 , the problem is formulated as Given the optimal solution d * , the net cost of electricity for route R is given by The variable neighborhood search heuristic was proposed by Mladenović and Hansen (1997) . Given the current solution S, the neighborhood structure is defined by a cyclic-exchange operator which selects N r routes from S to form an exchange cycle. The cyclic-exchange selects in each route R i , a random number υ i of consecutive nodes which form an exchange block. These blocks are then reversed and exchanged between the routes. Each exchange of blocks forms a new neighboring solution S which may be feasible or infeasible. Figure 1 shows an example of the cyclic exchange operator where N r = 3. The selected routes are R i , R j and R k with υ i = 2, υ j = 3, υ k = 2, respectively. The three blocks on the left form an exchange cycle. The blocks are reversed and transferred forming the three new routes R i , R j , and R k on the right. T emp of this annealing phase is initialized to T emp 0 such that a solution with cost f gen (S ) that is κ times worse than the current best solution will be accepted with a probability of 50%. The temperature is then linearly decreased by a factor after each VNS iteration so that in the last 20% of the iterations, the temperature is below 0.0001. These parameters are the same as the ones used by Schneider et al. (2014) . Tabu search (TS) is applied to every solution S generated by VNS. The neighborhood solutions are generated using the three widely-used operators for tabu search: 2-opt * , exchange, and relocate, as well as the stationInRe operator that is discussed in Schneider et al. (2014) . The operators are visualized in Figure 2 where the nodes that are on the same row are travelled by an EV before the operator is applied, the dashed arrows are edges to be removed, and the stripped and shadowed nodes are the nodes that are selected by the algorithm. The operators are applied as follows: • 2-opt*: Select two routes and remove one edge from each of them. Connect the first part of the first route with the second part of the second route and vice versa. • Exchange: Exchange the positions of two nodes. The two nodes could either be in the same route or in different routes. • Relocate: Select one route and remove one node from this route. Reinsert the selected node at another position. The new position could either be in the same route or in another route. • StationInRe: Perform insertion or removal of a station node. At each TS iteration, the four operators are applied on S and the candidate solutions are filtered against the tabu list. The solution that results in the maximal decrease in the generalized cost f gen (S) is selected. The reinsertion of the removed edges which led to the selected solution, is prohibited for a fixed number of tabu iterations called tabu tenure. The tabu tenure for each deleted edge is randomly selected from an interval [v min , v max ]. The procedure is repeated until no improvement can be made or for at most η tabu iterations. The best solution that is obtained is denoted as S . This section presents extensive computational results to evaluate the proposed Lagrangian relaxation approach and the VNS/TS hybrid heuristic. Table 4 ). The instances that start with "R" are random instances where customers are uniformly distributed, and those that start with "C" are clustered instances in which customers are clustered into small groups. The customer distribution of the "RC" instances is a mixture of random and clustered distributions. The EV fleet is homogeneous, each with a cargo capacity Q = 200, a battery capacity B = 270, and a range of 150 km. The travel speed v is set to a constant value of 30 km/h, and the batteries are recharged at a constant speed α = 1.8. The discharging speed is assumed the same as the charging speed. The choice of these values are discussed in more details in the case study. The locations and time windows in the Schneider instances are normalized values which makes it difficult to relate the charging and discharging time to real-world values so as to estimate the potential costs and gains. We adjust the instances in a way similar to Schiffer and Walther (2018) . 12:00 AM 7:00 AM 6.5 6.5 7:00 AM 11:00 AM 9.4 8.0 11:00 AM 5:00 PM 13.4 10.0 5:00 PM 7:00 PM 9.4 8.0 7:00 PM 12:00 AM 6.5 6.5 In particular, we set the maximum distance from the depot to customers/stations as 100 km and convert all other distances to km proportionally for each instances. In addition, since fully charging an EV requires B = 270 minutes (4.5 hours), we set the planning horizon as 5am − 12am ([0, 1140] in minutes) so that an EV has enough time to recharge at night. We scale all time windows proportionally to minutes to fit the planning horizon. We find through computational experiments that the time windows in the scaled instances are relatively tight which prevent vehicles from detouring to perform charging and discharging activities. We thus further relax the time windows to three periods, morning (5am − 12pm), afternoon (12pm − 6pm), and evening (6pm − 12am), based on their service start time. For charging and discharging, the length of each period is set to one hour, i.e. δ = 60 minutes. The energy prices which are shown in Table 3 that an EV can make profits by discharging at peak hours (11:00 AM -5:00 PM) and recharging the battery later. Note that the reward rates are also economically beneficial to the grid because the discharging reward that the grid pays to the EV owner is lower than the corresponding market price. All the tests are performed on a MacBook Pro running OS X 10.13.6, using a single 2.30 GHZ CPU, 16 GB of RAM and CPLEX 20.1.0.0 is used as an optimization solver. The time limit is set to 7200 seconds and the memory limit is set to 10 GB. The algorithms are implemented as single thread codes in Python. For the VNS/TS hybrid heuristic, the penalty parameters β tw , β batt , β cargo are set to 10 and the number of tabu iterations per round is set as η tabu = 30. The number of VNS iterations η vns and the early stopping criterion η early are varied with respect to the number of customers included. For the instances with 5, 10, and 15 customers, η vns is set to 10, 20, and 30 respectively. η early is set to 10 for instances with more than 10 customers and 5 otherwise. For the cyclic operator, N r is set to 2 when the fleet consists of 3 EVs or less and to 3 otherwise. The length of each exchange block is randomly selected from {1, 2, 3}. The upper and lower bounds of the tabu tenure are v min = 5 and v max = 15. Parameter κ for the annealing mechanism is set to 0.5. The results for CPLEX, the VNS/TS heuristic, the Lagrangian relaxation and the LP relaxation are presented in Table 4 . The upper and lower bounds achieved by the four algorithms are presented in the "UB" and "LB" columns, respectively. The "Time" column indicates the computational time in seconds. The "Gap" columns indicate the gap between the corresponding upper bound and the lower bound obtained by the Lagrangian relaxation. Column "Best Iter" indicates the VNS iteration during which the best solution was found. For the small instances with 5 customers, the VNS/TS heuristic outperforms CPLEX in terms of solution time for the majority of the tested instances. In terms of bound quality, the Lagrangian relaxation obtains the optimal solution for every instance, while the LP relaxation yields bounds Table 5 Performance of VNS/TS heuristic on medium instances of poor quality. The VNS/TS heuristic obtains the optimal solution for 9 out of 12 instances with the largest gap being 4.30%. As the number of customers increases to 10, the performance of CPLEX worsens significantly where only 1 out of the 12 instances are solved to optimality within the time and memory limits. CPLEX is able to find feasible solutions for 11 out of 12 instances before termination yet the gaps are on average 34.31%, which is significantly higher than the average gap of the VNS/TS heuristic 0.85%. The Lagrangian relaxation obtains bounds for 9 out of the 12 instances and the gap compared to the upper bound that is found by the heuristic is no more than 4.45%. We note that solving the Lagrangian subproblem is computationally challenging for the instances with 10 customers and no solution is obtained within the time limits for 3 out of the 12 tested instances. The lower bounds obtained by the LP relaxation are far away from the optimal values, which partially explains the low efficiency of CPLEX on this problem. For the instances with 15 customers, both CPLEX and the Lagrangian relaxation fail to solve the majority of the instances. CPLEX finds a feasible solution for only 1 out of the 11 instances. VNS/TS identifies solutions for all the instances within 317 seconds on average. The results for the VNS/TS heuristic on medium instances are shown in Table 5 . The VNS/TS heuristic is able to solve all the instances with less than 35 customers within 2 hours (20 iterations), allowing EV operators to implement the proposed approach for real-world tasks that do not require dynamic EV dispatching. Nevertheless, the proposed VNS/TS heuristic contains a large number of MIPs, which hinders its scalability. Future research efforts could possibly be made to re-formulate the embedded MIPs as quadratic knapsack problems and implement dynamic programming methods to solve them more efficiently. In this section, the proposed EVRPTW-TP model along with the VNS/TS solution heuristic are evaluated on a case study of an EV fleet providing online grocery delivery services, which attracted great attention amid the COVID-19 pandemic in the Kitchener-Waterloo (KW) region in Ontario, Canada. Through computational experiments, we particularly investigate the impact of electricity pricing schemes, time windows, winter range reductions, length of charging/discharging time period and fleet size on the fleet's routing and charging/discharging behaviors based on which we provide managerial insights. We consider a local grocery store that uses EVs to provide online grocery delivery services. This business model has been increasingly popular in recent years (Begley et al. 2020) , especially during the COVID-19 pandemic where people avoided going to supermarkets and local stores. In the KW region, retail giants, such as Walmart, Zehrs, and T&T Supermarket, as well as local small businesses are all providing such services. Customers would browse available items online and place the orders. Delivery providers perform touch-free delivery during the predetermined time slot. In this case study, we only consider orders that are scheduled for delivery on the next day and accordingly apply EVRPTW-TP to optimize the routing and charging/discharging of the EV fleet for the planning horizon of 5am -12am. For the purpose of the case study, we first consider a small fleet consisting of 3 EVs. The fleet size is then increased to 6 in order to investigate the impact on fleet operations (Section 7.6). Service is provided to 30 customers in the designated area, which reflects a real situation for a small business. The EVs in the fleet are homogeneous with identical range and cargo capacity. The EV that is considered has a range of 150 kilometers with a 32.4 kWh battery and a cargo capacity of 200 order units. The travel time between locations is calculated based on Google Maps. The geographical locations of the actual customers, the stations, and the depot are shown in Figure 3 as blue, red, and green dots, respectively. The planning horizon is first divided into 3 time periods: morning (5am -12pm), afternoon (12pm -6pm), and evening (6pm -12am). Results with different time windows are presented in Section 7.3. The service time at each customer location is set to 10 minutes. The 7 selected stations are level-2 stations which are compatible with most EV models. The power supply of a level-2 charger is 240V AC/30A. One hour of charge using a level-2 The solution time for the cases is on average 45 minutes which allows the fleet operator to set the routing, charging and discharging schedule at night. In this section, we study the impact of the electricity pricing on the EV fleet operations. As shown in Table 6 , we consider four electricity pricing schemes: Scheme A: The price varies between on-peak, mid-peak, and off-peak hours. The reward for energy discharging is the same as the price for energy charging. The price is only lower in the evening hours. Scheme D: Only energy charging is allowed (energy discharging to the grid is not allowed). In scheme A, the charging rates are the real hydro rates that were in effect between May 1, 2018 and Nov 1, 2019. We note that the province of Ontario currently does not have a set reward for injecting electricity to the grid. Thus we evaluate the case where the reward for discharging is lower than the price of electricity (Scheme A) and the case where the reward is the same as the price of energy (Scheme B). We also evaluate the case where discharging is not allowed (Scheme D). The electricity price periods are defined by the Ontario Energy Board differently for the summer and the winter (see Figure 4 ). We thus evaluate both summer and winter periods. The results obtained by solving EVRPTW-TP under the different electricity pricing schemes are shown in Table 7 . The first two columns present the total distance traveled and the total electricity cost of the fleet. The total cost is decomposed into charging cost and discharging reward illustrated in the following two columns. The negative numbers in the "Total Cost" column indicate that the EV fleet makes a net profit due to selling energy. The last three columns present the number of hours spent by the fleet on charging and discharging electricity while on route and the number of station visits of the fleet, respectively. The first insight is that allowing energy discharging is economically positive to EV fleet owners. When discharging is not allowed (Scheme D), the fleet's daily energy cost is $2.02, while the gain reaches $2.80 per day (Scheme B, winter) when discharging is allowed. The daily profits due to energy discharging could be as large as $4.82 (the difference between schemes B and D), corresponding to an annual profits of $1, 759.30. According to Fernández et al. (2013) and Millner (2010) the EV battery life span is around 4.6 -5.5 years under 100% depth of discharge (DoD) daily usage in mild temperature. Therefore, considering the 8-year warranty offered by most EV producers, then enabling discharging can result in 2.5 -3.4 years of reduction in battery life span. The total profits the fleet is able to gain from discharging before its end of life is $7,978 -$9,539. Assuming EV battery price of 163 $/kWh (Scott 2020) , the battery degradation cost is $4,951 -$6,734 for a fleet of 3 EVs, each with a battery capacity of 32.4 kWh. The profits gained from discharging can cover the battery degradation cost. Considering that the EV drivers are paid a fixed monthly salary, the additional profits assist to reduce the EV fleet's operational cost. Not surprisingly, electricity prices impact fleet operations. In particular, since the reward rate in scheme B is much higher than in scheme A, EVs are more willing to discharge during peak hours and detour to recharge their batteries later. The travel distance under scheme B is, on average 32 kilometers longer than under scheme A. EVs also spend 1 to 3 more hours discharging energy under scheme B than under scheme A. Moreover, the overall electricity cost varies significantly between these two cases. For schemes A and B, the total electricity cost is lower in winter than in summer. Taking scheme B as an example, in Figure 5 each blue line depicts one EV battery level through the day and the charging/discharging behaviors are highlighted in red. The green dashed lines represent the discharging reward rates. In summer, in order to perform discharging during peak hours (11am -5pm), at least one EV has to detour during the day because some customers should be served before 12pm (the end of the morning period). However, in winter, since the on-peak periods are 7am -11am and 5pm -7pm all the EVs can stay at the depot before 11am to perform discharge and leave to perform delivery services after that. As a consequence, the travel distance in winter is shorter than in summer for both schemes. Accordingly, EV operators can adjust their time slot settings with respect to the season. For instance, reducing the morning service period may decrease the overall electricity cost. Accounting for the changes in time-of-use prices and the corresponding time between the seasons is therefore important when designing the service time slots. Thus, next we evaluate the impact of the time windows on the EV operations. In this section, we consider three variations of the time windows. We first consider 2-hour time windows, then the 3-periods presented in the previous section, and finally eliminate the time windows. Currently, 2-hour time windows are the most commonly used for grocery delivery in the KW area to allow customers to know the time that they would be served more precisely, thus reducing their waiting times and enhancing user experience. This of course comes at the expense of the flexibility of the EV fleet in terms of routing and scheduling. The results presented in Table 8 and Figure 6 reflect the trade-off between timely service and operational flexibility. Under the 2-hours setting, customers that are geographically close might book time slots that are very far away from each other. In that case, an EV cannot always wait at a customer's location until the beginning of the time slot of the nearest customer. Instead, the EV will travel to a customer that is relatively distant but with a time slot starting earlier, leaving the nearest customer to be served by another EV. When a 3-period time window is used, more EVs are able to serve the customers in an order based on their geographical locations. As shown in Figure 6 , two EVs have to travel to Cambridge under the 2-periods setting, while only one has to do so under the other two settings. The increased flexibility leads to a decrease in the total travel distance from over 210 kilometers in the 2-hours setting to around 150 kilometers in the 3-period setting. The driving distance decreases to 131.81 kilometers for the case of no time window constraints. In terms of charging and discharging activities, intuitively the fleet has the greatest flexibility in scheduling discharging in the case where no time windows are used. In both winter and summer, the fleet is able to discharge for 12 hours in total during the planning horizon. For the "3-period" case, the fleet performs 9 hours of discharging, while for the 2-hour time windows, the fleet performs 6 hours of discharging in summer, and 7 hours in winter. As a consequence, as shown in Table 8 , the cost of electricity is 1.85 dollars and 1.64 dollars for summer and winter respectively under the 2-hours setting. A profit of 1.17 dollars and 1.12 dollars is generated when no time windows are enforced in summer and winter, respectively. This difference can be seen as the price of providing "timely service" delivery. In this section, we investigate the impact of range reduction under low ambient temperature. According to the data presented by Yuksel and Michalek (2015) and Lohse-Busch et al. (2013) , an EV's driving range could drop by 50% in severe weather conditions. We thus consider EV ranges from 100% to 50% of its stated range as shown in Table 7 .4. We note that battery capacity is assumed to be fixed, we vary EV range by adjusting the energy consumption rate. For example, increasing energy consumption rate by 100% is equivalent as reducing EV range by 50%. We also assume the charging and discharging processes are not affected by ambient temperature. Not surprisingly, Table 7 .4 shows that the fleet's total energy cost increases as EV range drops. When EVs are operated at close-to-full ranges (90% or more), the fleet can discharge for 10 hours during peak periods and still have the capacity to visit charging stations 3 times to charge the energy back during off-peak hours. Consequently, the discharging rewards cancel out all the charging costs. When EV range drops to 70-80%, an EV can not make its detour for en-route charging to have the same amount of energy being discharged during peak periods. It thus discharges one hour less and cancels the en-route station visit, which leads to higher total costs but less driving distance. If the EV range drops to 50-60% of its stated range, the fleet is still able to discharge for 8 hours, yet the total cost increases to around 2 dollars per day as more charging is needed. In this section, we investigate the impact of the length of each charging/discharging time period. Shortening time periods enhance the fleet's operational flexibility, allowing fleet operators to consider schedules of higher resolution. For instance, consider an EV that visits a station at 11:15am and plan to discharge for one hour during the visit, the EV would have to wait for 45 minutes before it starts discharging if δ = 60, while the waiting time is 15 and 0 minutes when δ = 30 and we observe that some EVs may finish serving assigned customers and return to the depot before the end of evening peak hours. These EVs usually don't have enough energy to discharge for one hour, thus can not make additional profits when δ = 60. Nevertheless, when δ = 15, the fleet is allowed to inject less energy back to the grid each time, hence can jointly discharge for 1.5 and 1.75 additional hours under summer and winter settings, respectively. Consequently, as presented in Table 10 , the total energy cost decreases as the value of δ decreases. Although shorter time periods lead to decreased energy cost, it also introduces extra computational burdens. The average computational time for the cases of δ = 60, 30 and 15 is 45, 110, and 160 minutes, respectively. In practice, fleet operators can choose the length of time period based on the available computational resources and the problem size. In this section, we consider the impact of EV fleet size on the overall electricity cost under the two restrictive time window settings, 2-hours and 3-periods. The number of EVs is varied from 3 to 6 and the corresponding results are shown in Table 11 and Figure 8 . Intuitively, increasing fleet size reduces the overall travel distance and provides more spare time to discharge energy to the grid. However, given the high acquisition cost, using an EV solely for the purpose of trading energy with the grid is currently not economical. Figure 7 shows the battery levels of an EV during summer and winter hours. In summer, the EV discharges its battery during the on-peak rate for 4 hours and recharges under the off-peak rate at night, bringing a profit of 1.00 dollar. Similarly in winter, the EV makes a profit of 1.09 dollars by discharging for 6 hours under on-peak rates and recharging for 2 and 4 hours under mid-peak and off-peak rates respectively. Adding one additional EV increases the profit by 1.01 dollars in summer hours an 1.09 dollars in winter hours. In that case the additional EV is acting as a static battery storage. The impact of fleet size on the travel distance is shown in Figure 9 . Each bar consists of several sub-bars highlighted in different colors each representing the travel distance of an EV in the fleet. As shown, for all the four cases, when the fleet size is 3, the fleets' overall travel distance is the highest. As the fleet size increases, the overall travel distance decreases. For the "2-hours" time window cases (red and blue bars), increasing the fleet size from 3 to 4 results in a relatively large reduction in the overall travel distance, while for other cases, the distance reductions are marginal. This observation suggests that the newly added EVs contribute to reducing the overall electricity cost mainly through grid ancillary services instead of improving the routing. Additional EVs thus are increasingly playing the role of a battery storage which is not economical under current energy prices. This paper presented the EV routing problem with time windows under time-variant electricity prices. The formulated problem jointly optimizes the routing of a fleet of electric vehicles along with the scheduling of the charging and discharging activities. The proposed framework permits the EV fleet to provide the primary service of delivering products to customers while minimizing the energy costs by scheduling energy charging while taking into consideration the varying energy prices throughout the planning horizon. Furthermore, the electric vehicles can store energy at periods with relatively low energy prices and inject it back to the grid at high price periods which helps in cost recovery and potentially realize profits. This framework demonstrates the capabilities of EV fleets in helping energy grids to smooth out the demand and reduce the gap between on-peak and off-peak by acting as energy storage. The proposed problem is formulated as a multiperiod optimization model that is challenging to solve to optimality. A Lagrangian relaxation that provides lower bounds is proposed. To find feasible solutions quickly which is important for implementation in practice, The proposed framework is a step towards realizing the full integration of transportation and energy networks. The joint optimization of both systems promises increasing opportunities to improve reliability and reduce costs. Nevertheless, it is not without limitations. Some assumptions, such as linear charging/discharging, constant energy consumption, are made to have a tractable problem. Moreover, battery degradation cost which could potentially be incurred by frequent charging and discharging is not taken into account as it leads to longer time horizon than the one considered. Incorporating these realistic characteristics presents interesting future research directions. From the perspective of solution method, the proposed Lagrangian relaxation approach and the VNS/TS heuristic both have sub-problems that are computationally demanding, which hinders their scalability. Future research efforts could be made to improve solution efficiency of both approaches. While in the framework that is proposed in this paper, the energy prices in a period are known in advance, a model that integrates dynamic and uncertain prices of energy will be investigated in future work. Furthermore, future work will particularly focus on integrating energy generated from renewable sources to further advance the development of green logistics systems. Solving the multivariant ev routing problem incorporating V2G and G2V options Optimal routing and scheduling of charge for electric vehicles: A case study Digital disruption at the grocery store Vehicle electrification: Status and issues A unified tabu search heuristic for vehicle routing problems with time windows The truck dispatching problem Operations research for green logistics-An overview of aspects, issues, contributions and challenges Exact algorithms for electric vehicle-routing problems with time windows The impact of transport electrification on electrical networks A green vehicle routing problem A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges Capacity fade and aging models for electric batteries and optimal charging strategy for electric vehicles Vehicle routing with time windows: Two optimization algorithms A conceptual framework for the vehicle-to-grid (V2G) implementation. Energy policy A review of lithium-ion battery state of charge estimation and management system in electric vehicle applications: Challenges and recommendations A variable neighborhood search heuristic for periodic routing problems Expanding electricvehicle adoption despite early growing pains Lagrangian duality applied to the vehicle routing problem with time windows Electric vehicles as a new power source for electric utilities Partial recharge strategies for the electric vehicle routing problem with time windows Impacts assessment of plug-in hybrid vehicles on electric utilities and regional us power grids, part 1: Technical analysis Sustainable operations management Plug-in hybrid vehicles and the vermont grid: a scoping analysis Ambient temperature (20°f, 72°f and 95°f) impact on fuel and energy consumption for several conventional vehicles, hybrid and plug-in hybrid electric vehicles and battery electric vehicle A metaheuristic to solve a location-routing problem with non-linear costs Modeling lithium ion battery degradation in electric vehicles Variable neighborhood search Canada's energy future Ontario Ministry of Transportation Impacts of plug-in electric vehicles in a balancing area An effective heuristic for the m-tour traveling salesman problem with some side conditions Strategic planning of electric logistics fleet networks: A robust location-routing approach The electric vehicle-routing problem with time windows and recharging stations Ever-cheaper batteries bring cost of electric cars closer to gas guzzlers Algorithms for the vehicle routing and scheduling problems with time window constraints Optimal charging strategies for unidirectional vehicle-to-grid The neglected social dimensions to a vehicle-to-grid (V2G) transition: a critical and systematic review Table 20-10-0021-01 new motor vehicle registrations An adaptive variable neighborhood search algorithm for a vehicle routing problem arising in small package shipping Joint routing and charging scheduling optimizations for smart-grid enabled electric vehicle networks A hybrid guided local search for the vehicle-routing problem with intermediate replenishment facilities Joint routing and scheduling for electric vehicles in smart grids with V2G. Energy Impact of plug-in-electric vehicles penetration on electricity demand, prices and thermal generation dispatch Plug-in hybrid electric vehicles and smart grids: Investigations based on a microsimulation Inside ups's electric vehicle strategy Electric vehicles' energy consumption measurement and estimation Electric vehicle route optimization considering time-of-use electricity price by learnable partheno-genetic algorithm Autonomous vehicle logistic system: Joint routing and charging strategy Effects of regional temperature on electric vehicle efficiency, range, and emissions in the united states. Environmental science & technology Bo Lin was supported by the Energy Council of Canada energy policy research fellowship and Bissan Ghaddar was supported by NSERC Discovery Grant 2017-04185. We greatly thank the referees, the Associate Editor, and the Editor for their thorough and thoughtful comments that helped us improve the quality of the paper.