A Heuristic Approach for a Real-World Electric Vehicle Routing Problem

To develop a non-polluting and sustainable city, urban administrators encourage logistics companies to use electric vehicles instead of conventional (i.e., fuel-based) vehicles for transportation services. However, electric energy-based limitations pose a new challenge in designing reasonable visiting routes that are essential for the daily operations of companies. Therefore, this paper investigates a real-world electric vehicle routing problem (VRP) raised by a logistics company. The problem combines the features of the capacitated VRP, the VRP with time windows, the heterogeneous fleet VRP, the multi-trip VRP, and the electric VRP with charging stations. To solve such a complicated problem, a heuristic approach based on the adaptive large neighborhood search (ALNS) and integer programming is proposed in this paper. Specifically, a charging station adjustment heuristic and a departure time adjustment heuristic are devised to decrease the total operational cost. Furthermore, the best solution obtained by the ALNS is improved by integer programming. Twenty instances generated from real-world data were used to validate the effectiveness of the proposed algorithm. The results demonstrate that using our algorithm can save 7.52% of operational cost.


Introduction
With industrial development and urbanization, air pollution in cities has become increasingly serious in recent years.For example, in December 2016, Beijing was covered by smog for six days, forcing authorities to announce the highest-level smog alert of 2016.A major cause of air pollution in cities is the large number of conventional motor vehicles with internal combustion engines that run on diesel or gasoline.They emit carbon oxides, nitrogen oxides and particulate matters that cause air pollution.To alleviate this hazardous situation and to reduce the number of conventional vehicles, green vehicles, such as alternative fuel vehicles and electric vehicles, which produce fewer emissions that contribute to air pollution than conventional vehicles, have gained much attention from city managers.Logistics companies are encouraged by city managers to use electric vehicles to construct a non-polluting and sustainable urban logistics system.For example, in January 2018, one of the largest Chinese online retailers, JD.com, used electric vehicles instead of conventional vehicles to deliver parcels in Beijing, and is going to replace all its conventional vehicles in the next two years.Electric vehicles are constrained by the capacity of their battery, the location of charging stations, and long charging times.The routing approaches for conventional vehicles are not applicable to electric vehicles, which leads to a new optimization problem called the electric vehicle routing problem (EVRP).
The basic EVRP is to obtain a set of routes with the minimum operational cost that serve customers' requests while satisfying the driving range limitations and charging requirements of electric vehicles.A real-world EVRP proposed by a logistics company in Wuhan, China is presented in this paper.In the problem, a fleet of heterogeneous electric vehicles departs from the distribution center, delivers or picks up customer parcels within predefined time windows, and finishes at the center.While in transit, a vehicle can return to the center to load parcels to be delivered or unload collected parcels, and it can also go to charging stations as well as the distribution center to fully charge its battery in a fixed time (e.g., 30 min).The objective is to minimize the sum of the acquisition costs of used vehicles, the travel costs of vehicles, the waiting costs at customers and centers, and the charging costs.In Figure 1, two types of electric vehicles, i.e., Vehicle 1 and Vehicle 2, are used to fulfill delivery or pickup requests of customers; the lines having the same color denote a visiting route of a vehicle in the day.The red route does not need to charge; the purple route requires two charges at two different stations, one of which it shares with the black route; and the yellow route with two trips needs to go back to the distribution center.Obviously, the problem combines the features of the capacitated VRP [1], the VRP with time windows [2,3], the heterogeneous fleet VRP [4], the multi-trip VRP [5], and the EVRP with charging stations [6].Considering its complexity, we must resort to meta-heuristics to efficiently solve the proposed problem.
Algorithms 2019, 12, x FOR PEER REVIEW 2 of 18 applicable to electric vehicles, which leads to a new optimization problem called the electric vehicle routing problem (EVRP).The basic EVRP is to obtain a set of routes with the minimum operational cost that serve customers' requests while satisfying the driving range limitations and charging requirements of electric vehicles.A real-world EVRP proposed by a logistics company in Wuhan, China is presented in this paper.In the problem, a fleet of heterogeneous electric vehicles departs from the distribution center, delivers or picks up customer parcels within predefined time windows, and finishes at the center.While in transit, a vehicle can return to the center to load parcels to be delivered or unload collected parcels, and it can also go to charging stations as well as the distribution center to fully charge its battery in a fixed time (e.g., 30 min).The objective is to minimize the sum of the acquisition costs of used vehicles, the travel costs of vehicles, the waiting costs at customers and centers, and the charging costs.In Figure 1, two types of electric vehicles, i.e., Vehicle 1 and Vehicle 2, are used to fulfill delivery or pickup requests of customers; the lines having the same color denote a visiting route of a vehicle in the day.The red route does not need to charge; the purple route requires two charges at two different stations, one of which it shares with the black route; and the yellow route with two trips needs to go back to the distribution center.Obviously, the problem combines the features of the capacitated VRP [1], the VRP with time windows [2,3], the heterogeneous fleet VRP [4], the multi-trip VRP [5], and the EVRP with charging stations [6].Considering its complexity, we must resort to meta-heuristics to efficiently solve the proposed problem.This paper focuses on developing a heuristic approach based on the adaptive large neighborhood search (ALNS) and integer programming (IP) to solve the proposed problem.In the algorithm, we devise a charging station adjustment heuristic based on labeling algorithm to obtain the optimal selection and location of charging stations in a given sequence of customers and a departure time adjustment heuristic to calculate the optimal departure time when vehicles should first leave from the distribution center.After the termination of the ALNS iterative process, the best solution is further improved by set-partitioning-based integer programming.In computational experiments, the effectiveness of the components of the proposed algorithm was validated by 20 instances generated from real-world data.
The rest of this paper is organized as follows.Section 2 provides a brief review of the pertinent literature about the EVRP.Section 3 outlines a definition for the proposed problem.Section 4 presents the solution approach for the problem.Section 5 reports the computational results.Section 6 presents the conclusions and future research directions.This paper focuses on developing a heuristic approach based on the adaptive large neighborhood search (ALNS) and integer programming (IP) to solve the proposed problem.In the algorithm, we devise a charging station adjustment heuristic based on labeling algorithm to obtain the optimal selection and location of charging stations in a given sequence of customers and a departure time adjustment heuristic to calculate the optimal departure time when vehicles should first leave from the distribution center.After the termination of the ALNS iterative process, the best solution is further improved by set-partitioning-based integer programming.In computational experiments, the effectiveness of the components of the proposed algorithm was validated by 20 instances generated from real-world data.

Literature Review
The rest of this paper is organized as follows.Section 2 provides a brief review of the pertinent literature about the EVRP.Section 3 outlines a definition for the proposed problem.Section 4 presents the solution approach for the problem.Section 5 reports the computational results.Section 6 presents the conclusions and future research directions.

Literature Review
Research on the VRP and its variants is increasing tremendously in the literature [7][8][9].Pelletier et al. [10] presented a survey of the existing research in the goods distribution with electric vehicles.This section mainly reviews the very recent developments of the EVRP and its variants.Schneider et al. [11] introduced the EVRP with time windows which considers the possibility of partial recharges at any station.They proposed a hybrid heuristic based on a variable neighborhood search (VNS) algorithm and a tabu search (TS) heuristic.Felipe et al. [12] exploited constructive and local search heuristics within a non-deterministic simulated annealing framework to solve an EVRP with multiple technologies and partial recharges.Sassi et al. [13] proposed an iterated TS for the mix fleet VRP with heterogeneous electric vehicles.Yang and Sun [14] presented an electric vehicle battery swap location routing problem that determines the location of battery swap stations and the routing of electric vehicles under battery driving range limitations.In the problem, vehicles can stop at battery swap stations to swap their battery for a fully charged battery.They proposed a four-phase heuristic called SIGALNS (Sweep heuristic, Iterated Greedy, Adaptive Large Neighborhood Search) and a two-phase TS-modified Clarke and Wright Savings heuristic to solve the problem.Wen et al. [15] addressed the electric vehicle scheduling problem in which a set of timetabled bus trips should be carried out by a set of electric vehicles.The vehicles can be recharged fully or partially at any recharging station.They proposed an ALNS algorithm to solve it.Keskin and Çatay [16] studied the EVRP with time windows and partial recharging and design an ALNS algorithm to solve it.Hiermann et al. [6] introduced the electric fleet size and mix VRP with time windows in which the vehicle is assumed to recharge to full capacity on visit of a recharging station.They solved the problem using the branch-and-price algorithm and the enhanced ALNS algorithm.Desaulniers et al. [17] presented four variants of the EVRP based on the allowable times of recharges per route and partial or full recharges.For each variant, they proposed the branch-price-and-cut algorithm with customized mono-directional and bidirectional labeling algorithms to solve it.Montoya et al. [18] considered the EVRP with a nonlinear charging function, which assumes that the battery-charge level is not a linear function of the charging time.They developed a hybrid meta-heuristic combining an iterated local search (ILS) and a heuristic concentration method to solve the problem.Schiffer and Walther [19] introduced a location-routing problem that considers routing of electric vehicles and positioning decisions for charging stations.Partial recharges are allowed in the problem.They presented an ALNS, which is improved by local search and labeling algorithms, and devised new penalty functions for neighborhood evaluation.Hof et al. [20] designed an adaptive VNS algorithm to solve the battery swap station location-routing problem proposed by Yang and Sun [14].Zhang et al. [21] devised an ant colony algorithm to solve the EVRP with recharging stations to minimize energy consumption.Macrina et al. [22] proposed an ILS heuristic to the mixed fleet vehicle routing problem with partial battery recharging and time windows in which the fleet is composed of electric and conventional vehicles.Froger et al. [23] proposed a new model, a heuristic, and an exact labeling algorithm for the EVRP with nonlinear charging functions.
From the literature review, we can find that differing charge strategies and objectives make the solution approaches not interchangeable and that problem-specific optimization components must be designed to accommodate the features of our problem, which motivates us to design a sophisticated heuristic for solving it.

Problem Description
The problem can be defined on a complete directed graph G = (V, A) with a set of nodes V = {0} ∪ N ∪ F and a set of arcs A = {(i, j) | i, j ∈ V, i = j}.Node 0 is the distribution center (also called the depot), the set N = {1, 2, . . ., |N|} is the set of customers, and the set F = {|N|+1, . . ., |N|+|F|} is the set of charging stations.The travel time and distance on arc (i, j) are denoted by t ij and d ij , respectively.Each customer i ∈ N has a load q i , which is positive for a pickup request and negative for a delivery request, a service time t i , and a time window [a i , b i ], where a i and b i are the earliest and latest allowable start service times, respectively.Each customer should be visited exactly once.A fleet of heterogeneous electric vehicles is available at Node 0. The fleet composes K different types of vehicles.Each vehicle type k = {1, 2, . . ., |K|} has a maximal load capacity Q k , a maximal driving range D k (due to limited battery capacity), a charging time t k c , a charging cost c k c per unit of time, an acquisition cost c k f if it visits any customer, and a travel cost c k v per unit of distance.Each vehicle departs from Node 0 after a given time a 0 , serves some customers, and finally ends at Node 0 before a specified time b 0 .The vehicle capacity Q k must be respected at customer i.While in transit, the vehicle can return to Node 0 to charge, reload parcels or unload parcels (i.e., multi-trips) if necessary.The total travel distance of a vehicle after the last charge should not exceed its maximal driving range D k .Each vehicle should start the service of a customer within the given time windows.If the vehicle arrives earlier than a i , a waiting cost proportional to the waiting time is incurred.Meanwhile, when the vehicle returns to Node 0 to prepare for the next trip (including charging), a given waiting time t w is spent, which also incurs a waiting cost proportional to the waiting time.The parameter c w is the cost per unit of waiting time.Each vehicle should go to one of the charging stations before it runs out of electricity and charge to full capacity in its predefined charging time t k c [24].The problem consists of determining a route plan for satisfying customers' demands while minimizing the sum of the acquisition costs, travel costs, charging costs, and waiting costs at customers and the distribution center.It contains four aspects: (1) the vehicle type assigned to each customer; (2) the visiting sequence of customers; (3) the charging station or distribution center and time chosen by the vehicle if it needs recharging; and (4) the departure time when the vehicle should first leave the distribution center.
For a detailed description and mathematical formulation of the proposed problem, interested readers are referred to similar papers [6,19].

Solution Approach
Considering the complexity of the proposed problem, we devise a heuristic approach based on the ALNS and IP to solve it.The ALNS was first proposed by Ropke and Pisinger [25].Its basic idea is to use a leaning mechanism to bias the selection of a variety of removal and insertion operators that are used to generate new solutions.It has been widely used to solve various variants of VRPs [26][27][28][29][30][31][32][33][34][35].In this paper, we incorporate the features of the proposed problem into the ALNS.Its details are described in Algorithm 1.
An initial solution s 0 is first generated at the beginning of the ALNS.Then, at each iteration, as on Line 5, a removal and insertion heuristic pair is chosen based on their scores and weights in previous iterations.On Line 6, a given number of customers are first removed from current solution s current using the chosen removal heuristic and put into the set of removed customers.The locations of charging stations and depot in the partial solution remain unchanged.Using the corresponding insertion heuristic, the removed customers are reinserted into the partial solution to generate a new solution s .On Lines 7 and 8, solution s is first improved by a local search, and then the charging stations and departure time of all routes are optimized.On Lines 9-17, if s is better than current best infeasible solution s inf * or current solution s current or meets the acceptance criteria, then it replaces s current .Accordingly, the route set R model is updated with new feasible routes in s .On Lines 18-20, the current best feasible solution s fea * is updated.Then, on Lines 21 and 22, the scores and weights of the selected removal and insertion heuristics are updated, and the penalty coefficients of the generalized cost function are also updated.The search is repeated until some stopping criterion is satisfied.The stopping criterion is that the runtime of the ALNS exceeds the given value T ALNS .Then, a set-partitioning mathematical model based on the routes in set R model is constructed and solved by optimization software.Finally, the solution obtained from the model is outputted as the solution of the problem.
Algorithm 1: Solution approach for the proposed problem.

Construction of Initial Solution
To generate an initial feasible solution, a sequential route construction heuristic is implemented.At each iteration, using identical unassigned customers, the routes associated with all vehicle types are created independently.For the construction of a route with respect to a vehicle type, the customer with the least cost increase is added into the partial route while satisfying all constraints.If one insertion violates the battery capacity (or driving range limitation), the algorithm attempts to insert the charging station closest to the customer before or after the customer.The insertion process continues until no more customers can be added into the route.The vehicle type with the least cost is selected.The customers in the corresponding route are then removed from the unassigned customer set, and the next iteration is performed until no unassigned customer is available.The newly generated routes are further optimized using the labeling algorithm and departure time adjustment heuristic in Sections 4.8 and 4.9.

Penalty Functions
Both feasible and infeasible solutions are allowed in the ALNS using a generalized cost function f gen (s) to evaluate a solution s.In addition to a term f obj (s) denoting the objective value, penalty terms for time window violations f tw (s), capacity violations f cap (s), and battery capacity violations f batt (s) are taken into consideration Factors α, β, and γ are used to weight the penalties.They are initialized with (α 0 , β 0 , γ 0 ) and dynamically adjusted between a given lower (α min , β min , γ min ) and upper bound (α max , β max , γ max ).
To achieve a good balance between diversification and intensification, the weights are multiplied by a factor ω if a penalty occurred in the last iteration, and are divided by ω if no penalty occurred in the last iteration.
A time window violation occurs at a customer if the vehicle arrives later than the latest allowable time for starting service which equals the arrival time minus the latest time for starting service.The total time window violation f tw (s) is the sum of violations of all customers.Because multiple trips and customers with delivery or pickup requests are allowed in a route, a capacity violation occurs if the maximal demand loaded by the vehicle exceeds its capacity before returning the depot.The capacity violation during one trip equals the maximal demand minus the capacity of vehicles.The total capacity violation f cap (s) is the sum of violations on all trips.A battery capacity violation occurs at a charging station or the depot if the distance traveled by the vehicle exceeds its battery capacity after the last charge.The total battery capacity violation f batt (s) is the sum of violations of charging stations and the depot.

Removal Heuristics
This section introduces four removal heuristics.These heuristics take a solution and an integer q as input.The output is a partial solution where q customers are removed and the order of other customers and charging stations remains the same.Moreover, the heuristics Shaw removal, worst removal, and historical knowledge node removal have parameters p Shaw , p worst and p historical that introduce some randomness in the selection of customers.

Random Removal Heuristic
This simple removal heuristic removes q customers selected randomly from current solution s.The idea of randomly selecting nodes helps diversify the search.

Shaw Removal Heuristic
The Shaw removal heuristic was first proposed by Shaw.Its general idea is to remove similar customers, as it can be expected that it is reasonably easy to reshuffle similar customers and create new, perhaps better, solutions.We use the distance d ij between customers i and j to define the relatedness R(i, j) between the two customers.The lower R(i, j) is, the more similar are the two customers.The first customer i to be removed from solution s is selected at random.This customer i is added into the set of removed customers D. Thereafter, at each iteration, one customer i * is randomly selected from D. An array L containing all visited customers from solution s not in D is constructed.This array is sorted according to increasing relatedness values R(i * , j).Then, a random number y is chosen from the interval [0, 1) and customer L[y pshaw |L|] is removed and put into D.This process is repeated until q customers are removed from solution s.

Worst Removal Heuristic
The worst removal heuristic attempts to remove customers with high cost and insert them at another position in the solution to obtain a better solution.Given customer i served by some vehicle in solution s, cost(i, s) = f gen (s) − f gen −i (s) where f gen −i (s) is the cost of the partial solution without customer i.At each iteration, an array L containing all visited customers from solution s not in D is constructed.This array is sorted according to descending costs.Then, a random number y is chosen from the interval [0, 1) and customer L[y pworst |L|] is removed and put into D.This process is repeated until q customers are removed from solution s.

Historical Knowledge Node Removal Heuristic
The historical knowledge node removal heuristic is derived from the one used in Demir et al. [26].It records the position cost of every customer i, defined as the sum of the distances between its preceding and following nodes, and computed as d i = d i−1 , i + d i, i+1 .At each iteration, the best position cost d i * is updated to be the minimum of all d i found so far.An array L containing all visited customers from solution s not in D is constructed.This array is sorted according to descending values of the maximum deviation between current and best position cost (i.e., d i − d i * ).Then, a random number y is chosen from the interval [0, 1) and customer L[y phistorical |L|] is removed and put into D.This process is repeated until q customers are removed from solution s.

Insertion Heuristics
To repair a partial solution, three insertion heuristics basic greedy insertion, deep greedy insertion, and regret-k insertion are employed [28].They try to reinsert the customers removed by removal heuristics into the partial solution without changing charging stations.Infeasible solutions are accepted by insertion heuristics.

Basic Greedy Insertion Heuristic
Let ∆ ik be the change in the generalized cost function value incurred by inserting customer i into route k at the position that increases the function value the least.Furthermore, let cost i = min k∈K {∆ ik } be the cost of inserting customer i at its best position overall which is also called the minimum cost position.This basic greedy insertion heuristic inserts customers in the set of removed customers D following their removal sequence.After one customer has been inserted at its best position, cost i is calculated again and the process is repeated until all removed customers are reinserted into the partial solution.

Deep Greedy Insertion Heuristic
The difference between this heuristic and the previous one lies in the insertion sequence of customers in D. In each iteration, the deep greedy insertion heuristic selects customer i having the minimum global cost (i.e., min i∈D {cost i }), and inserts it into its minimum cost position.

Regret-k Insertion Heuristic
Let r ik denote the route for which customer i has the kth lowest insertion cost, that is, ∆ i,r ik ≤ ∆ i,r ik for k ≤ k .In the regret heuristic, the regret value c i * is defined as the difference in the cost of inserting customer i in its best route r i1 and its second-best route r i2 , i.e., cost i * = ∆ i,r i2 − ∆ i,r i1 .In each iteration, the regret heuristic chooses customer i having the maximum global regret value, i.e., min i∈D {cost i * }, and inserts it into its minimum cost position.Ties are broken by selecting the lowest cost insertion.The process is repeated until all removed customers are reinserted into the partial solution.
The above technique can be extended naturally to define a class of regret heuristics by computing the difference in cost of inserting customer i in its best, 2nd best, kth best route, where k is a user-defined parameter.The resulting heuristic is called the regret-k heuristic.The regret-k heuristic selects customer i based on max i∈D {∑ k j=2 ∆ i,r ik − ∆ i,r i1 } and then inserts it in its least cost position.Regret-2, regret-3, regret-4, and regret-all are used in our ALNS.The regret-all heuristic considers all routes when selecting customers to be inserted.

Choosing a Removal and Insertion Heuristic
At each iteration, a pair of removal and insertion heuristics is selected from removal and insertion heuristics to generate new solutions.Each heuristic i has a weight w i and a score π i .The values of w i and π i are initialized as one and zero, respectively, at the beginning of the ALNS.The score of a heuristic is updated as follows.If a pair of removal and insertion heuristics finds a new global best solution, their scores are increased by σ 1 ; if it finds a new solution that has not been accepted before and is better than the current one, their scores are increased by σ 2 ; if it finds a new non-improving solution that has not been accepted before but is accepted in current iteration, their scores are increased by σ 3 .
The ALNS iterative process consists of a number of segments.Each segment contains 100 iterations.At the start of each segment, the scores of all the heuristics are set to zero.When one segment ends, the weight of heuristic i to be used in the next segment is updated as w i = w i (1 − γ) + γπ i /θ i where θ i is the number of times heuristic i is used during the last segment and γ ∈ [0, 1] is a reaction factor that controls how quickly the weight adjustment reacts to changes in the effectiveness of the heuristics.The choice of the removal and insertion heuristics is determined by the well-known roulette-wheel mechanism.Specifically, given k insertion (or removal) heuristics with weights w i , i ∈ {1, 2, . . ., k}, heuristic j is selected with a probability of w j / ∑ k i=1 w i .Note that the insertion heuristic is chosen independently of the removal heuristic (and vice versa).

Acceptance Criteria
Similar to most literature about the ALNS, we use the acceptance criteria from simulated annealing to avoid getting trapped in a local minimum.That is, solution s is always accepted if it is better than s, i.e., f gen (s ) ≤ f gen (s); otherwise, solution s is accepted with a probability of e -(f gen(s ) -f gen(s))/T , where T denotes the current temperature.Following Ropke and Pisinger [25], T is initialized as T start and is decreased every iteration using the equation T = T • c, where c ∈ (0, 1) is the cooling rate.Let f gen (s 0 ) denote the generalized cost function value of the initial solution s 0 .The initial temperature T start equals -w•f gen (s 0 )/ln0.5,where w is the initial temperature control parameter.

Local Search
The local search procedure uses a list of operators, including inter-and intra-route 2-opt, swap, and relocation [36], to improve the solution generated by removal and insertion heuristics.The procedure works in a cyclic manner and uses a first improvement strategy.Each operator is explored until no further improvement is found, after which the next operator is chosen and explored.When the last operator of the list is explored, the search starts again from the first operator.The iterative process continues until a local minimum is attained in all operators.Note that the operators are executed between customers instead of charging stations.

Labeling Algorithm
The removal and insertion heuristics and local search do not alter the charging stations of routes in the solution.To obtain the optimal charging stations or intermediate depots for a given sequence of customers, we first construct an auxiliary graph by adding all charging stations and the depot into the customer sequence.Taking a customer sequence {0, 1, 2, 3, 4, 5, 0} as an example, its auxiliary graph is illustrated in Figure 2, where node E i denotes the charging station and m is the number of charging stations.We need to find a path with minimal cost that starts at Node 0, visits all customers and returns to Node 0. Such a problem can be solved by the labeling algorithm [6,37].

Acceptance Criteria
Similar to most literature about the ALNS, we use the acceptance criteria from simulated annealing to avoid getting trapped in a local minimum.That is, solution s′ is always accepted if it is better than s, i.e., fgen(s′) ≤ fgen(s); otherwise, solution s′ is accepted with a probability of e -(fgen(s′) -fgen(s))/T , where T denotes the current temperature.Following Ropke and Pisinger [25], T is initialized as Tstart and is decreased every iteration using the equation T = T • c, where c ∈ (0, 1) is the cooling rate.Let fgen(s0) denote the generalized cost function value of the initial solution s0.The initial temperature Tstart equals -w•fgen(s0)/ln0.5,where w is the initial temperature control parameter.

Local Search
The local search procedure uses a list of operators, including inter-and intra-route 2-opt, swap, and relocation [36], to improve the solution generated by removal and insertion heuristics.The procedure works in a cyclic manner and uses a first improvement strategy.Each operator is explored until no further improvement is found, after which the next operator is chosen and explored.When the last operator of the list is explored, the search starts again from the first operator.The iterative process continues until a local minimum is attained in all operators.Note that the operators are executed between customers instead of charging stations.

Labeling Algorithm
The removal and insertion heuristics and local search do not alter the charging stations of routes in the solution.To obtain the optimal charging stations or intermediate depots for a given sequence of customers, we first construct an auxiliary graph by adding all charging stations and the depot into the customer sequence.Taking a customer sequence {0, 1, 2, 3, 4, 5, 0} as an example, its auxiliary graph is illustrated in Figure 2, where node Ei denotes the charging station and m is the number of charging stations.We need to find a path with minimal cost that starts at Node 0, visits all customers and returns to Node 0. Such a problem can be solved by the labeling algorithm [6,37].The label algorithm uses labels to construct a partial path from Node 0 to any customer in the sequence.The label  = ( ,  ,  ,  ,  ,  ,  , ) at node i is defined as:  is the arrival time at node i;  and  are the total pickup and delivery demand after the last departure from Node 0, respectively;  is the maximal demand after the last departure from Node 0;  is the total travel distance after the last charge;  is the cost from the starting point to node i;  is the adjustment cost used and calculated in the dominance rule; and i is the last reached node.The label algorithm uses labels to construct a partial path from Node 0 to any customer in the sequence.The label l = t a i , q p i , q d i , q max i , d i , f o i , f a i , i at node i is defined as: t a i is the arrival time at node i; q p i and q d i are the total pickup and delivery demand after the last departure from Node 0, respectively; q max i is the maximal demand after the last departure from Node 0; d i is the total travel distance after the last charge; f o i is the cost from the starting point to node i; f a i is the adjustment cost used and calculated in the dominance rule; and i is the last reached node.
At the starting point 0, t a 0 , q p 0 , q d 0 , q max 0 , d 0 , f o 0 , f a 0 , 0 are initialized to zero.When a label t a i , q p i , q d i , q max i , d i , f o i , f a i , i of customer i is extended to customer j, (m+2) labels t a j , q p j , q d j , q max j , d j , f o j , f a j , j are generated, including one label associated with a path directly traveling from customers i to j using the extension rule in Algorithm 2, m labels associated with each path going though one charging station using the extension rule in Algorithm 3, and one label associated with a path going through the depot using the extension rule in Algorithm 4.
Algorithm 2: Extension rule for directly traveling from customers i to j.
if customer j is a delivery node then 3: q d j := q d i + q j , q p j := q p i 4: if customer j is a pickup node then 6: q d j := q d i , q p j := q p i + q j 7: end if 8: q max j := max{q max i , q d j + q p j } 9: if j is the depot then 11: end if Algorithm 3: Extension rule considering charging stations between a pair of nodes. 1: if customer j is a delivery node then 3: q d j := q d i + q j , q p j := q p i 4: if customer j is a pickup node then 6: q d j := q d i , q p j := q p i + q j 7: end if 8: if j is the depot then 11: end if Algorithm 4: Extension rule considering depot between a pair of nodes.
1: t a j := max{a i , t a i } + t s i + t i0 + t w + t 0j 2: if customer j is a delivery node then 6: q d j := q j , q p j := 0 7: end if 8: if customer j is a pickup node then 9: q d j := 0, q p j := q j 10: end if 11: q max j := max{0, q d j + q p j } 12: An exponential number of labels will be generated during the extension process.Therefore, a dominance rule is essential to remove the labels that cannot lead to optimal paths.The rule is defined as follows.

Definition 1 (Dominance rule):
If two labels l = t a i , and at least one of the inequalities is strict.
The arrival time of label l is earlier than that of label l , which does not guarantee that the cost of the former is smaller than that of the latter during the later extension.Therefore, we adjust the arrival time of label l to the same as that of label l , compute the adjustment cost f a i , and identify the dominance relationship between the two labels.The adjustment cost f a i is computed as follows.
The notation used in the labeling algorithm is listed as: L i is the set of labels associated with customer i; l i is one label associated with customer i; and v i and n s are the ith node and the number of customers in the sequence {v 0 , v 1 , . . .,v i , . . .,v n s , v n s +1 }, where v 0 and v n s +1 denote depot 0. The pseudo code of the labeling algorithm is presented as Algorithm 5.

Departure Time Adjustment Heuristic
The waiting cost in the objective function depends on the first departure time from the depot.We devise a departure time adjustment heuristic to compute the optimal departure time for a given route whose first departure time is zero.The heuristic relies on the difference between the arrival time and time window of each customer.Taking a route {v 0 , v 1 , . . .,v i , . . .,v n s , v n s +1 } as an example, where v i and n s are the ith node and the number of nodes except the starting and ending nodes.The details of the heuristic are described as Algorithm 6.It first sets the departure time from the depot t a 0 to zero on Line 1, and then on Lines 2-11 calculates the difference between the arrival time and time window of each customer v i , i.e., t e v i max{a v i − t a v i , 0} and t l v i max{b v i − t a v i , 0}.The adjustment iterative process is described on Lines 12-35.In each iteration, the first customer that incurs a waiting cost and the maximal allowable adjustment time t min are found on Lines 14-22, and the values of t a 0 , t e v i , and t l v i are updated on Lines 23-33.The iterative process continues until the value of t l v i at customer v i is zero.
1: t a 0 := 0 2: for i = 0 to n s do // initialize t e v i and t l v i 3: else if v i+1 is a charging station then 9: end for 12: repeat 13: t min := +∞, i min := 0 14: for i = 0 to n s do 15: if t e v i > 0 then // find the first customer which incurs waiting cost 16: i The set R model is populated during the iterations of the ALNS using feasible routes of the local optimum at each iteration.When adding the routes in R model , we check their feasibility and remove the routes with the same set of customers and larger route costs.After the iterative process of the ALNS, we use the set R model to construct a set-partitioning model as follows.
where c r is the route cost of route r; α ri is 1 if customer i is visited by route r and otherwise 0; and x r is a 0-1 decision variable that is 1 if route r is selected in the optimal solution and otherwise 0. When the travel time and distance between two customers respect the triangle inequality, the "=" in the constraints in Equation ( 4) can be replaced by "≥" (i.e., the set-covering model).
We solve such a model using optimization software and expect that some better solutions missed by the ALNS can be found.To reduce the runtime of solving the model, the best solution obtained by the ALNS is used as an initial solution of the model.Meanwhile, the optimization software continues until the optimal solution of the model is obtained or the maximal runtime exceeds the predefined time T IP .

Computational Experiments
In this section, we first introduce the construction of test instances from real-world data provided by a logistics company in Wuhan, China.Then, we report the values of major parameters in the proposed algorithm.Finally, we validate the effectiveness of the proposed algorithm by testing the instances.The algorithm was coded in C++ language.The mathematical model was solved by IBM ILOG CPLEX 12.6.3.All experiments were run on a computer with Intel Xeon E5-1620 V4 3.50 GHz CPU and 32 GB RAM under Windows 10 operating system.

Construction of Test Instances
The test instances were obtained from real-life data provided by a logistics company in Wuhan, China.The company uses two types of electric vehicles to provide logistics services for customers.The details of electric vehicles are listed in Table 1.The company also has detailed information about customers and charging stations, such as their locations.To generate the instances, we randomly selected approximately 100 customers and 10 charging stations around a distribution center.The loading or unloading time (t i ) at each customer is 20 min.If a vehicle returns to the distribution center, it must wait 40 min (t w ) to prepare for the next trip.The waiting cost per minute (c w ) is 0.5.Twenty instances were generated to validate the effectiveness of the proposed algorithm.Instances are represented in the following format: R1 indicates the first instance in the set.

Parameter Settings
The parameters of the proposed algorithm are listed in Table 2. To set all parameters, each parameter in turn took several values, while the others were fixed.We ran the proposed algorithm 10 times for each parameter setting, and the parameter setting producing the best average results was selected.Table 2 also shows the best parameter setting for the proposed algorithm.According to a series of preliminary experiments, the number of customers removed by removal heuristics (q) had the largest influence on the performance of the proposed algorithm.To give a sensitivity analysis to this parameter, six different intervals were tested: [0.05|N|, 0.25|N|], [0.1|N|, 0.3|N|], [0.05|N|, 0.35|N|], [0.1|N|, 0.4|N|], [0.15|N|, 0.35|N|], and [0.15|N|, 0.45|N|].At each iteration of the ALNS, q was selected uniformly at random in corresponding interval.Table 3 shows the average of the objective function value of all instances in different intervals.The results prove that [0.1|N|, 0.4|N|] was among the best of all intervals.

Performance Analysis of ALNS
In the proposed algorithm, the best solution obtained by the ALNS is further improved by set-partitioning-based integer programming.To validate the effectiveness of such an algorithm design, we compared the results only using the ALNS with those improved by integer programming (IP).Meanwhile, in the original framework of ALNS, the local search component does not exist.Therefore, we also tested the original ALNS without using the improvement strategies.For each instance, we ran the algorithms ten times.In Table 4, the column "Customers" presents the number of customers included in the instance; the column "Charging Stations" denotes the number of charging stations vehicles, customers are always visited by type-1 vehicles according to our experiments.Taking the best solution of R1 as an example, it uses 10 type-1 vehicles, charges six times, waits 509 min, and travels 949.707 kilometers.Then, the total acquisition cost is 2000 Yuan, the total charging cost is 270 Yuan, the total waiting cost is 254.5 Yuan, and the total travel cost is 9497.07Yuan.As shown in Figure 3, the total travel cost is 79% of the total operational cost, indicating that designing reasonable visiting routes is essential for companies to save operational cost.

Comparison with Company's Algorithm
The logistics company uses optimization software developed by a consulting company to address the proposed EVRP.To further demonstrate the performance of our algorithm, we compared the results obtained by it with that by optimization software.Due to the confidentiality agreement, the core algorithm embedded in the software cannot be described in detail.The results of instances obtained by the software and the ALNS are listed in Table 5.In the table, the column "Customers" presents the number of customers included in the instance; the column "Charging Stations" denotes the number of charging stations that can be used in the instance; the column "Company's Algorithm" shows the results of instances obtained by the software; the column "Proposed Algorithm" lists the best results of instances obtained by our algorithm; and the column "Gap" demonstrates the relative gap between the software and our algorithm and equals (Proposed Algorithm − Company's Algorithm)/Proposed Algorithm × 100%.

Comparison with Company's Algorithm
The logistics company uses optimization software developed by a consulting company to address the proposed EVRP.To further demonstrate the performance of our algorithm, we compared the results obtained by it with that by optimization software.Due to the confidentiality agreement, the core algorithm embedded in the software cannot be described in detail.The results of instances obtained by the software and the ALNS are listed in Table 5.In the table, the column "Customers" presents the number of customers included in the instance; the column "Charging Stations" denotes the number of charging stations that can be used in the instance; the column "Company's Algorithm" shows the results of instances obtained by the software; the column "Proposed Algorithm" lists the best results of instances obtained by our algorithm; and the column "Gap" demonstrates the relative gap between the software and our algorithm and equals (Proposed Algorithm − Company's Algorithm)/Proposed Algorithm × 100%.Table 5 shows that for all instances the proposed algorithm outperformed the optimization software used in the company.The minimal and maximal relative gaps between the two approaches were 4.49% (R6) and 9.91% (R2), respectively, and the average relative gap was 7.52%, indicating that using our algorithm could save 7.52% of operational cost.We showed the solutions obtained by our algorithm to the managers in the company, and they were very satisfied with these solutions, which can be directly used in practice.These results clearly show the effectiveness of our proposed algorithm.

Conclusions and Future Research Works
This paper develops a heuristic approach to solve a real-world EVRP proposed by a logistics company in Wuhan, China.In the approach, the ALNS is followed by solving a set-partitioning mathematical model that strengthens its effectiveness.Meanwhile, to reduce charging costs and waiting costs, two problem-specific heuristics, including the charging station adjustment heuristic and the departure time adjustment heuristic, are devised.The performance of the proposed algorithm was validated by 20 test instances generated from real-world data.The contributions of this paper are summarized as follows.First, this paper investigates a new variant of VRPs that originates from a real-world application and combines the features of different VRPs.Second, an efficient approach based on the ALNS and IP is developed to solve it.In the approach, a charging station adjustment heuristic and a departure time adjustment heuristic are devised to accommodate the features of the proposed problem.Third, the superiority of the proposed approach over the original ALNS is validated through computational experiments.
This work also offers several directions for future research.First, different charge strategies, such as allowing partial charges, will be included in the proposed problem.Second, different charge techniques, such as different charging modes, will be considered in the proposed problem.Finally, different meta-heuristics, such as the genetic algorithm [38] and the VNS [39], will be developed to enrich the solution approaches to the proposed problem.

Figure 1 .
Figure 1.Illustration for the proposed problem.

Figure 1 .
Figure 1.Illustration for the proposed problem.

Figure 2 .
Figure 2. Auxiliary graph in the labeling algorithm.

Figure 2 .
Figure 2. Auxiliary graph in the labeling algorithm.

Figure 3 .
Figure 3. Percentages of costs of the best solution of R1.

Figure 3 .
Figure 3. Percentages of costs of the best solution of R1.

Table 1 .
Parameters of electric vehicles.

Table 2 .
Values for the parameters of the proposed algorithm after the tuning phase.

Table 3 .
Parameter settings for q.

Table 4 .
Results of instances obtained by the proposed algorithm.

Table 5 .
Results of instances obtained by company's algorithm and the proposed algorithm.