A Benchmark Test Suite for the Electric Capacitated Vehicle Routing Problem

Severa1 logistic companies started utilizing electric vehicles (EVs) in their daily operations to reduce greenhouse gas pollution. However, the limited driving range of EVs may require visits to recharging stations during their operation. These potential visits have to be addressed, avoiding unnecessary long detours. We formulate the electric capacitated vehicle routing problem (E-CVRP), which incorporates the possibility of EVs visiting a recharging station while satisfying the delivery demands of customers. The energy consumption of the EVs is proportional to their cargo load which is an important constraint in real-world logistics applications. A new set of benchmark instances is proposed for the E-CVRP. As solution methods to these new benchmarks, we apply the ant colony optimization metaheuristic method and an exact method. Experimental results on the ECVRP demonstrate the high complexity of the problem and the efficiency of the applied metaheuristic solution method.


I. INTRODUCTION
In recent years, there is a growing interest of logistic companies in utilizing electric vehicles (EVs) for their daily operations due to the greenhouse effect. Therefore, a problem of routing a fleet of EVs has emerged, namely the electric vehicle routing problem (EVRP) [1]. The EVRP aims to find the best possible routes for a fleet of EVs so that a set of customers is visited once and only once by an EV, ensuring that the EVs will never run out of energy. Several variations of the EVRP have been introduced that mainly differ on their constraints such as the EVRP with time windows [2], the EVRP with pick up and delivery [3], the EVRP with multiple depots [4], the EVRP with service times [5] and the EVRP with backuals [6]. Refer to the comprehensive survey in [7] for further details.
The common constraint of all EVRP variations, is the consideration of visiting a recharging station due to the limited driving range of EVs. In this work, we formulate the electric capacitated vehicle routing problem (E-CVRP) as a mixedinteger linear program (MILP) where each EV has a limited cargo load to satisfy the delivery demands of the customers This work has been supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No. 739551 (KIOS CoE) and from the Government of the Republic of Cyprus through the Directorate General for European Programmes, Coordination and Development. and the energy consumption rate of the EV is proportional to that cargo load. Existing EVRP formulations assume that the energy consumption rate is fixed [2], [8], or it is estimated using longitudinal dynamics model [9]. In this work we emphasize on the effect the cargo load of the EV has on its energy consumption.
In addition, a new benchmark set of E-CVRP problem instances is proposed to represent the complexity of the problem. The allocation of the charging stations in the benchmark instances is optimized to ensure energy coverage for EVs from any customer locations. Experiments are conducted in order to provide an initial set of results to this new benchmark set. As E-CVRP extends the well-known capacitated vehicle routing problem (CVRP), the high complexity of the problem makes exact solution methods inefficient for solving realistically sized problem instances [10]. Hence, as an approximation solution method, we apply the ant colony optimization (ACO) metaheuristic [11], which is successfully used to solve world-scale instances of similar problems [12], [13]. The newly generated E-CVRP benchmark set consists of small-sized instances, that they can be solved (within a reasonable time) using the proposed MILP formulation with the Gurobi solver [14] in order to asses the performance of approximation solution methods (i.e., the deviation of the approximation solution from an upper bound solution that could be the optimal), and realistic-sized instances in order to study the efficiency of approximation solution methods.
The rest of the paper is organized as follows. Section II describes the E-CVRP and provides the MILP formulation of the problem. Section III describes the two aforementioned approaches for solving the E-CVRP. Section IV describes a method to generate E-CVRP problem instances from existing CVRP problem instances. Initial results obtained on the newly generated benchmark set are presented in Section V. Section VI gives the conclusions and future work of the paper.

II. THE ELECTRIC CAPACITATED VEHICLE ROUTING PROBLEM (E-CVRP) A. Problem Description
The E-CVRP is a challenging N P-hard combinatorial optimization problem as it is an extension of the conventional 978-1-7281-6929-3/20/$31.00 ©2020 IEEE CVRP problem incorporating additional hard constraints [10]. The problem can be described as follows: given a fleet of EVs with limited battery charge level and limited cargo load capacity, we need to find the best possible route for each EV, starting and ending to the central depot, satisfying the delivery demands of a set of customers.
The E-CVRP can be defined on a complete weighted graph G = (N, A), where N = {0} ∪ I ∪ F is a set of nodes and A = {(i, j) | i, j ∈ N, i = j} is a set of arcs connecting these nodes. A non-negative value d ij is associated with each arc which represents the euclidean distance between nodes i and j. Node 0 denotes the central depot. The set I ⊂ N denotes the set of customers, where each customer i ∈ I is assigned a positive value δ i indicating the customer's delivery demand 1 . The set F ⊂ N denotes the set of β i node copies of each charging station i ∈ F (i.e., |F | = i∈F β i ), which are used to permit multiple visits to each charging station i ∈ F (if required) [15]. The upper bound on the number of node copies for each charging station is equal to β i = 2|I|, because in the worst case an EV for each customer is needed and a visit to a charging station before and after serving it is required [16].
Each EV 2 has a cargo load, u i (0 ≤ u i ≤ C), where C is the maximal cargo load of an EV and u i is the remaining cargo load of an EV on arrival at node i. In addition, each EV has a battery charge level y i (0 ≤ y i ≤ Q), where Q is the maximal battery charge level of an EV and y i is the remaining battery charge level of an EV on arrival at node i. Each traveled arc consumes the amount h i d ij of the remaining battery charge level, where h i denotes the energy consumption rate of an EV traversing that arc. The energy consumption rate is proportional to the current cargo load of an EV (i.e., the heavier the EV the more energy will be consumed) defined as follows: where r is a constant value of the energy consumption rate, u i is the remaining cargo load of an EV on arrival at node i, and C is the maximal cargo load of an EV. There are other minor factors that may affect the energy consumption rate [5], but in this work we focus on the weight of the EV which is one of the most important factors [17]. The objective function of the E-CVRP is to find a set of EV routes that minimize the total distance traveled where: • Every customer is visited once and only once by exactly one EV. • All EVs start and end at the depot. • All EVs always start fully charged from the depot, e.g., y 0 = Q. • All EVs always start with fully loaded cargo from the depot, e.g., u 0 = C. • For every EV route the total delivery demand of customers does not exceed the EV's maximal cargo load C. 1 For the central depot and charging stations the demand is set to zero, i.e., δ 0 = 0 ∧ δ i = 0, ∀i / ∈ I. 2 A homogeneous fleet of EVs is considered.
• The battery charge level on each traveled arc is never negative. • EVs always leave a charging station fully charged. • The charging stations can be visited multiple times by any EV.

B. MILP Problem Formulation
To solve the E-CVRP, a binary decision variable is defined, i.e., x ij , ∀i, j ∈ N , denoting whether arc (i, j) is traversed by an EV (i.e., x ij = 1) or not (i.e., x ij = 0). The mathematical model of E-CVRP is formulated as: where Eq. (2a) defines the E-CVRP objective function, Eq. (2b) enforces the connectivity of customer visits, Eq. (2c) handles the connectivity of visits to recharging stations, Eq. (2d) establishes the flow conservation by guaranteeing that at each vertex, the number of incoming arcs is equal to the number of outgoing arcs, Eq. (2e), Eq. (2f) and Eq. (2g) guarantee delivery demand fulfillment at all customers by assuring a non-negative cargo load upon arrival at any vertex, Eq. (2h) ensures that EVs start with full cargo load from the depot, Eq. (2i), Eq. (2j) and Eq. (2k) ensure that the battery charge level never falls below 0, Eq. (2l) ensures that EVs start fully charge from the depot, and Eq. (2m) is the binary decision variable described previously.

III. SOLVING THE E-CVRP
In order to provide an initial set of results (e.g., optimal values, upper bounds, or best known values) for each problem instance of the new benchmark set we use the MILP formulation in Eqs. (2) and the ACO metaheuristic, described in Algorithm 1, as solution methods for the E-CVRP.

A. Using the MILP Formulation
Although the formulation defined in Eqs. (2) fully describes the E-CVRP problem, it cannot be handled by standard mathematical optimization solvers, such as the Gurobi solver [14]. This is because that formulation does not adhere to any MILP standard form, because of the constraints in Eq. (2i) and Eq. (2j) as they involve a product of a continuous and a binary variable, that is, u i x ij . Recall that the current cargo load u i is used to calculate the variable energy consumption rate h i in Eq (1). To cope with this issue, the continuous variable m ij is introduced to replaced the product, m ij = u i x ij , so that it will be equivalently transformed to a set of MILP inequalities according to the big "M" notation [28] as follows: with M ≥ C.
To summarize, the problem presented in Eqs.
(2) can be reformulated into a MILP problem by replacing the product of variables (i.e., u i x ij ) in Eq. (2i) and Eq. (2j) with m ij and adding Eqs. (3). As a result, we get a transformed mathematical program that involves only linear constraints.

B. Using the ACO Metaheuristic
The MAX -MIN Ant System (MMAS) [19] variant is used which is one of the most-studied ACO variants and has already proven its good performance on relevant EVRPs [5], [20]. In this section we describe the application of MMAS to the E-CVRP, including the proposed method to construct feasible solutions.
1) Initialization: A colony of ω ants is initially positioned at the central depot, i.e., component 0. All the solution components of the problem are associated with a pheromone trail value which is uniformly initialized at the start of the execution as follows: where τ ij is the pheromone trail value associated with arc (i, j) connecting solution components i and j, and τ 0 is the initial pheromone trail value. A good value for τ 0 was found to be 1/ρC nn , where ρ is the evaporation rate (more details are provided later on) and C nn is the solution quality of the solution generated by the nearest-neighbor heuristic [21].
2) Solution Construction: Each ant k represents a complete E-CVRP solution T k (i.e., the routes of all EVs) and makes selections biased by the existing pheromone trails and some heuristic information associated with the solution components of the problem, until all the components (i.e., customers) are selected.
The probability distribution with which ant k selects the next customer component j from solution component i is defined as follows: where τ ij and η ij = 1/d ij are, respectively, the existing pheromone trail and the heuristic information available a priori between components i and j. α and β are the two parameters which determine the relative influence of τ ij and η ij , respectively. N k i is the set of unselected customers components for the k-th ant adjacent to component i. Note that the depot and charging station components are not included in the N k i set (the selection of these components is described later in Section III-B3).
The selected customer component j must satisfy the EV's capacity constraints defined in Eqs. (2e) and (2g), and the energy constraints defined in Eqs. (2i), (2j), and (2k). It is worth mentioning that the construction method proposed in this section can be utilized in other constructive or local search heuristics to construct feasible solutions when addressing similar EVRPs. The following section describes in more detail the handling of the aforementioned constraints.
3) Constraints Handling: To construct a feasible solution the next customer component j to be selected, using Eq. (5), must satisfy all the following criteria: 1) The delivery demand of the customer component j, i.e., δ j , must not cause the cargo load to fall below zero, that is: 2) The energy consumption to travel from solution component i to customer component j must not cause the energy level, i.e., y i , to fall below zero, that is: 3) The customer component j must have at least one charging station or the central depot within its energy range [5], that is: When the next customer violates the delivery demand criteria in Eq. (6), then the depot component is added to the E-CVRP solution to close the route of the EV (i.e., denotes the return of the EV to the central depot). In case this return to the depot violates the energy criteria in Eq. (7), then the EV visits a recharging station prior to its return to the depot. When the next customer violates any of the energy criteria in Eq. (7) or the range criteria in Eq. (8), then the closest energy recharging station s from the current solution component i to the (potentially) next customer j is selected as follows: and is added to the E-CVRP solution (i.e., denotes the visit to an energy recharging station). Then, the next customer to for k = 1 to ω do 8: T k ← InsertDepotComponent() 10: while ant k has not selected all customers do 11: i ← GetLastInsertedComponentFrom(T k ) 12: if j satisfies all the criteria then 15: T k ← InsertCustomerComponent(j) 16: else if j violates the capacity constraint in Eq.
if C k is better than C ib then if arc (i, j) ∈ T best then 40: end if 49: end while 50: OUTPUT: T bs be visited will be most probably j if all three criteria defined in Eq. (6), Eq. (7) and Eq. (8) are now satisfied, otherwise another customer satisfying all the aforementioned criteria will be selected. The possibility of not having enough energy to travel to any energy recharging station or back to the central depot is eliminated because the customer component must always satisfy the range criteria in Eq. (8).
4) Pheromone Update: In the MMAS variant [19], [21], [22], the pheromone trails are updated by first decreasing the pheromone trails on all arcs (using pheromone evaporation), and then increasing the pheromone trails on the arcs of the solution constructed by the best ant (using pheromone deposit). The pheromone evaporation is applied as follows: where ρ (ρ ∈ (0, 1]) is the evaporation rate. After pheromone evaporation, the best ant deposits pheromone on the arcs of its solution components as follows: where ∆τ best ij = 1/C best is the amount of pheromone that the best ant deposits and C best is the quality of the best solution T best . The "best" ant that is allowed to deposit pheromone may be either the best-so-far ant 3 , in which C best = C bs , or the iteration-best ant, in which case C best = C ib . These two ants are allowed to deposit pheromone in an alternate way. More precisely, the iteration-best ant deposits pheromone at each iteration and the best-so-far ant deposits pheromone occasionally (i.e., every 25 iterations in this work, more details are provided in [22]). Pheromone evaporation enables the algorithm to gradually decrease the pheromone trails on the arcs that belong to low-quality solutions so that the pheromone trails on the arcs that belong to high-quality solutions consist of higher values (e.g., they will be increased regularly because of pheromone deposit).
Since only the best ant is allowed to deposit pheromone the search is focused on a specific area of the search space. Hence, there is consequent danger of getting trapped in search stagnation 4 from the early stages of the optimization process. To address this issue, the lower and upper limits of the pheromone trail values are explicitly imposed, preventing in this way the excessive growth of pheromone trails on the arcs of the best ant, as follows: where τ min and τ max are, respectively, the minimum and maximum pheromone trails. The τ max value is bounded Algorithm 2 Convert CVRP to E-CVRP 1: INPUT: Classical CVRP benchmark instance 2: Read initial benchmark problem of size n 3: Define maximal battery charge level (i.e., Q) 4: Add charging stations 5: OUTPUT: Electric CVRP benchmark instance by 1/ρC bs , where C bs is the solution quality of the bestso-far ant 5 . The τ min value is set to τ min = τ max (1 − n √ 0.05)/((avg − 1) n √ 0.05) where avg is the average number of different choices available to an ant at each solution construction step and n is the size of the problem instance. Note that the τ max and τ min values are updated whenever whenever a new best-so-far ant is discovered.
In addition, all the pheromone trails are uniformly reinitialized to the τ max value whenever the stagnation behavior occurs or when no improved solution is found for a given number of iterations. The stagnation behavior is detected using the λ-branching factor that measures the distribution of the pheromone trail values [23]. The λ-branching factor is given by the number of arcs incident to node i satisfying the following condition: τ ij ≥ τ i min + λ(τ i max − τ i min ), where τ i max and τ i min are, respectively, the maximum and minimum pheromone values on arcs incident to node i, and λ ∈ [0, 1] is a constant. The average λ-branching factor from all nodes' λ-branching factors indicates whether the search has entered stagnation behavior or not.

A. Converting CVRP to E-CVRP
In this section, we provide the algorithmic steps followed to convert a CVRP problem instance to an E-CVRP problem instance with the basic procedures outlined in Algorithm 2. More specifically, at the first step we import the classical CVRP problem instance where each customer node i ∈ I represents a point in the Cartesian plane (e.g., lines 1-2). In the sequel, the maximum battery charge level Q is defined as follows: whered is the maximum euclidean distance between the depot and all other customer nodes, i.e.,d = arg max i∈I {d 0i }, and γ is a user defined constant. The goal of stations placement (Algorithm 2, line 4) is to find the minimum number of charging stations together with their positions so that all customer nodes are covered by at least one charging station. This procedure ensures that problem instances with remote customer nodes stay feasible. In particular, it is assumed that each charging station covers an area of a circle with a specific radius R = φQ, where φ is a user selected constant. Increasing the value of φ results in a larger coverage area of charging stations. Each customer i ∈ I can be covered from at least one charging station placed 5 C bs is initialized with the solution quality of the estimated optimal tour C nn . within the area of the circle that has radius R with customer i being the center of the circle. All potential positioning points are in the convex-hull area of all customer points that can be computed by the Quickhull method [18]. The set P denotes any potential discrete point in the above convex area in which a charging station can be potentially placed (excluding all the points of customer nodes). Furthermore, the set A i , ∀i ∈ I, determines all the potential charging stations points in set P that cover each customer i ∈ I (i.e., the area of the circle).
To optimally solve the described charging station placement problem we introduce a MILP formulation that aims to select the minimum number (together with their position) to cover each customer node. Let the variable ψ j ∈ {0, 1}, j ∈ P , denote whether the potential point j ∈ P is selected to place a charging station (ψ j = 1) or not (ψ j = 0). The related mathematical formulation is as follows: where Eq. (14a) is the objective function as described previously, Eq. (14b) is responsible to ensure that each customer i ∈ I is covered from at least one charging station and Eq. (14c) defines the binary decision variable. Finally, in line 5 of Algorithm 2 the defined maximum battery charge level value and the positioning of the charging stations are used to generate the E-CVRP benchmark instance.

B. Description of E-CVRP Problem Instances
Four different sets of CVRP benchmark instances are used to generate the new E-CVRP benchmark instances in the experiments 6 , including the popular instances of Christofides and Eilon [24], Christofides et al. [25], and Fisher [26], and the recent instances of Uchua et al. [27]. The details for converting a CVRP instance to an E-CVRP instance are given in Section IV. The converted E-CVRP benchmark problem instances are described in Tables I, II, III, IV, where the columns (starting from the left) denote the name 7 of the problem instance, the number of customers (#customers), the number of charging stations (#stations), the minimum number of EV routes a solution can have 8 (#routes), the maximum cargo load of an EV (C), and the maximum battery charge level of an EV (Q). Note that all E-CVRP benchmark instances consist of a single central depot. To generate reasonable (i.e., enforcing EVs to visit charging stations) E-CVRP instances the γ and φ parameters are set to γ = 2 and φ = 0.125. Increasing the   value of γ and decreasing the value of φ results in a higher number of charging stations.

A. Experimental Setup
In order to provide an initial set of results to the new benchmark, experiments have been conducted with the MMAS approach described in Section III-B and the MILP approach defined in Section III-A. The experiments were performed under a Linux System with an Intel Core i7-3930K 3.20GHz processor with 12MB cache and 16GB RAM.
The MILP approach was solved using the Gurobi solver [14] and a time limit of approximately 3 weeks was imposed. The MMAS approach performs 25000n evaluations 9 for 10 independent executions (with different random seeds), where n is the problem size. The total distance traveled of the solution with the best quality together with the CPU time are recorded. The colony size of MMAS was set to ω = n ants. The evaporation rate was set to ρ = 0.02 and the two parameters of the decision rule were set to typical values as follows: α = 1 and β = 2.

B. Analysis of the Results
The experimental results for small-scale problem instances of MMAS and MILP approaches are given in Table V, whereas the experimental results for large-scale problem instance are given in Table VI. Note that the MILP approach was not able to provide a solution for the large-scale problem instances within the pre-determined time limit, and, thus, only the results of the MMAS approach are reported in Table VI. "U B" is the best upper bound within the time limit found by the MILP approach (in some cases the given upper bound is also the optimal). "best", "worst", "mean" and "stdev" are, respectively, the minimum, maximum, average and standard deviation values found by MMAS from its 10 independent runs. Furthermore, "t avg " is the CPU time required (averaged over 10 runs) for MMAS to find the solution with the best quality. The "U B" and "best" values also include the number of routes that the solution consists of within the brackets.
From Table V it can be observed that the MILP approach can find better solutions than the MMAS 10 approach on most of the smallest problem instances. Nevertheless, the deviation of the solution quality of the MMAS approach from the upper bound solution found by the MILP approach on E-n29-k4-s7, E-n30-k3-s7 and E-n35-k3-s5 problem instances (which is also the optimal solution for these cases) is very small. However, the MMAS approach requires a few seconds whereas the MILP approach requires weeks. As the problem size increases it can be observed that the MMAS approach can find better solutions than the MILP approach (e.g., on the E-n60-k5-s9 problem instance) because the pre-determined time appears to be short. This observation was expected considering the high complexity of the E-CVRP.
From Table VI it can be observed that the MMAS approach was able to find a solution for all E-CVRP problem instances, requiring a few minutes for most problem instances (with less than 300 nodes) and couple of hours for the remaining problem instances. This observation demonstrates the efficiency of the MMAS approach against the MILP    instances, e.g., the E-n37-k4-s4 problem instance in Ta-ble V. In addition, it verifies that ACO algorithms are effective solvers for these type of problems. Finally, the median solution (from the solutions of the 10 runs) generated by the MMAS approach is plotted in Fig. 1 for the E-n29-k4-s7 problem instance and the optimal solution generated by the MILP approach is plotted in Fig. 2 for the same problem instance. From Fig.1 it can be observed that the solution generated by the MMAS approach consists of four routes, in which all of them contain a single visit to a recharging station. On the contrary, from Fig. 2 the optimal solution consists of four routes in which one route contains two visits to recharging stations and one visit for each of the remaining three routes. Note that the solutions generated by both approaches can be operated by four EVs simultaneously or by a single EV sequentially. Although the median solution generated by the MMAS approach requires less visits to recharging stations than the optimal solution, its solution quality is inferior. This observation shows the challenge when addressing the E-CVRP because fewer detours to recharging stations does not necessarily lead to better solution quality.

VI. CONCLUSIONS
In this work, we present a new variation of the vehicle routing problem for determining minimum distance routes for EVs. The E-CVRP considers limited battery charge level and limited cargo load capacity for the EVs. The energy consumption of the EVs is based on the cargo load of the EVs, that is, the heavier the EV the more energy will be consumed. The problem is formulated as a MILP and it is solved using an exact solution method and an approximation solution method. The experimental results on a set of newly generated benchmark E-CVRP instances represent the complexity of the problem because the applied approaches experience challenges either in terms of efficiency or solution quality.
For future work, it will be interesting to design heuristic methods that are problem specific to this E-CVRP in order to further improve the solution quality. In fact, there is a lot of room for improvements on these benchmark problem instances considering the initial experimental results presented in this work.