Hatice Calik
email: hatice.calik@kuleuven.be

Ammar Oulamara
email: oulamara@loria.fr

Caroline Prodhon
email: caroline.prodhon@utt.fr

Said Salhi
email: s.salhi@kent.ac.uk

The electric location-routing problem with heterogeneous fleet: Formulation and Benders decomposition approach

Keywords: Environmental Logistics, Electric Vehicle Routing, Recharging Station Location, Integer Programming, Benders Decomposition

In this paper, we focus on a problem that requires location of recharging stations and routing of electric vehicles in a goods distribution system. The goods are disseminated from a depot and distributed to the customers via electric vehicles with limited capacity. Differently from the classical vehicle routing problem, the vehicles have battery restrictions that need to be recharged at some stations if a trip is longer than their range. The problem reduces to finding the optimal locations of the recharging stations and their number to minimize the total cost, which includes the routing cost, the recharging cost, and the fixed costs of opening stations and operating vehicles. We propose a novel mathematical formulation and an efficient Benders decomposition algorithm embedded into a two-phase general framework to solve this environmental logistics problem. Phase I solves a restricted problem to provide an upper bound for the original problem which is later solved in Phase II. Between the two phases, an intermediate processing procedure is introduced to reduce the computations of the Phase II problem. This is achieved by a combination of the Phase I upper bound and several lower bounds obtained via exploiting the underlying network structure. Our approach solves the problem in a general setting with nonidentical stations and vehicles by allowing multiple visits to the stations and partial recharging. The computational study provides both managerial and methodological insights.

Introduction

The transport sector is responsible, to a large extent, for energy consumption and greenhouse gas emissions. According to the European Environment Agency (2018), the energy consumption of road transport increases by 32% from 1990 to 2016 in the EEA-33. To tackle environmental and energy challenges, several countries are considering the prospect of carbon neutrality over the next 30 years, with the objective of discouraging the sale of vehicles emitting greenhouse gases. The implementation of such a strategy has already begun with the introduction of low-emission zones (LEZ), where vehicles with higher emissions either cannot enter the area or have to pay a high penalty. For instance, the traffic pollution charge in London LEZ is £100 per day for larger vans and minibuses and rises to double this amount for lorries, buses, and coaches. Vehicles with alternative fuels, such as electric vehicles (EVs) and hydrogen vehicles, provide credible solutions for achieving the carbon neutrality target.

Unlike the hydrogen vehicle, which is currently at the experimental stage, and consequently having an exorbitant cost, the EV has reached an industrial maturity that makes it competitive compared to the combustion vehicle. However, as indicated by [START_REF] Davis | A methodology to evaluate the competitiveness of electric delivery trucks[END_REF] and [START_REF] Sassi | Electric vehicle scheduling and optimal charging problem: complexity, exact and heuristic approaches[END_REF], EVs are still facing weaknesses related to their availability, purchase price, and battery management. From a logistics point of view, there are still weaknesses that are worth pointing at. These include (i) The limited choice of light duty EVs offered by the car industry. These vehicles are mainly needed in the last-mile logistics.

(ii) The limited EV driving range. For instance, for light duty EVs, the range is between 120km and 180km. Note that the range can depend on topology of the road as well as weather and driving conditions.

(iii) The long charging time. The time to fully charge a vehicle can take up to 8 hours depending on the capacity of the battery pack and chargers' level.

(iv) The lack of availability of charging infrastructures in existing road networks. Although all these weaknesses are manageable in practice, the cost of EV presents a barrier to their extensive use. An opportunity to reduce the vehicle's price is focusing on the development of those markets that are ready to adopt such a green-based strategy. Such markets allow a largescale production of EVs which can consequently lead to the reduction of vehicle costs. Last-mile logistics provide this opportunity to speed up the market penetration of EVs. In such markets, an EV has the advantage of meeting the requirement of low-emission zones that are mainly located in city centers. Here, the distances covered in last-mile logistics are either within its range or it requires one charging session along the route only. Furthermore, even though the acquisition cost for EVs is usually higher than the combustion engine vehicles, this difference can be offset at the operational cost of EV usage. This is because a high utilization of EVs favors their TCO (Total Cost Ownership) since their operating costs (maintenance, tax, fuel, and depreciation) are low compared to those of their counterparts.

In this paper, we consider a goods distribution system that utilizes EVs. This is a system where the operating companies have access to their own recharging stations (private) or subscribe to a contract that warrants access (without queuing restrictions) to certain recharging stations which have to be selected. Similar business models are considered by [START_REF] Yang | Battery swap station location-routing problem with capacitated electric vehicles[END_REF] for battery swap stations and by Schiffer and Walther (2017a,b) for recharging stations. In these types of business models, the operators need to decide on both the location and the routing aspects. As location and routing decisions are interdependent, they need to be handled simultaneously to operate an overall system in the most profitable way [START_REF] Salhi | The effect of ignoring routes when locating depots[END_REF]. It may be argued that it is difficult to integrate operational decisions such as routing into strategic decisions like locating facilities. Though this is a critical issue, studies dealing with this dilemma showed that an intelligent way of incorporating the results of the integration can be very useful. For instance, [START_REF] Salhi | Consistency and robustness in location-routing[END_REF] conduct a robustness analysis leading to a conclusion that integrated models constantly provide higher quality solutions and they are as reliable as 'locate first -route second' methods.

In our study, we consider the vehicles to depart from a single depot. We also assume there is a sufficient number of charging stations and electrical grid capacity. This is to ensure that all vehicles are fully charged before their departure from the depot. However, we may need to recharge them during their trips if the total energy consumption to visit certain customers is larger than the battery capacity. Once a station is opened, it might be visited multiple times by any vehicle. As we allow partial recharging, the vehicles do not need to be fully recharged. Besides, we do not impose any restrictions on the types of stations or vehicles. In other words, we allow the use of heterogeneous vehicles and stations that might have different location-dependent costs.

The problem is to decide on the number and location of stations, the number of vehicles needed, the amount of recharging needed for each vehicle, and the route(s) for visiting all the customers. The objective is to minimize the total cost which includes the variable cost of routing and recharging as well as the fixed costs of opening stations and operating vehicles. In this study, we develop an exact method based on a new formulation which utilizes disaggregated commodity flows to express sub-tour elimination and capacity restrictions [START_REF] Yaman | Formulations and valid inequalities for the heterogeneous vehicle routing problem[END_REF][START_REF] Baldacci | Routing a heterogeneous fleet of vehicles, in: The vehicle routing problem: latest advances and new challenges[END_REF][START_REF] Salhi | Adaptation of some vehicle fleet mix heuristics[END_REF][START_REF] Salhi | Incorporating vehicle routing into the vehicle fleet composition problem[END_REF]. There are several applications in the literature where flow based formulations with capacity constraints are successfully solved using a Benders decomposition approach. These include the hub location-routing problem studied by [START_REF] De Camargo | A new formulation and an exact approach for the many-to-many hub location-routing problem[END_REF] and network design problems by [START_REF] Fortz | An improved Benders decomposition applied to a multi-layer network design problem[END_REF], [START_REF] Botton | Benders decomposition for the hop-constrained survivable network design problem[END_REF][START_REF] Calik | A Benders decomposition based framework for solving cable trench problems[END_REF]. See also other relevant Benders decomposition applications for the location of EV recharging stations in car sharing systems [START_REF] Fortz | A Benders decomposition method for locating stations in a one-way electric car sharing system under demand uncertainty[END_REF], under probabilistic travel range [START_REF] Lee | Benders-and-price approach for electric vehicle charging station location problem under probabilistic travel range[END_REF], and with plug-in hybrid EVs [START_REF] Arslan | A benders decomposition approach for the charging station location problem with plug-in hybrid electric vehicles[END_REF]; in the survey by [START_REF] Costa | A survey on Benders decomposition applied to fixed-charge network design problems[END_REF] for fixed-charge network design problems; and in the book by [START_REF] Birge | Introduction to stochastic programming[END_REF] for stochastic programming problems. This motivates us to apply a Benders decomposition algorithm leading to very successful results for solving the small size instances which are shown to be challenging by the preceding study of Schiffer and Walther (2017b). To the best of our knowledge, the heterogeneous fleet feature which makes the problem relatively much more challenging to handle by exact or heuristic methods has not been considered within the location routing framework of combined recharging station location and EV routing problems. The mathematical models introduced by Schiffer and Walther (2017) for the homogeneous fleet variants cannot be utilized or easily adapted to solve the heterogeneous fleet variants. Given that the models for homogeneous variants have difficulty in solving even small instances with 5 customers, there is a clear need for more efficient exact methods which can solve relatively larger instances (e.g. 10-15 customers) and enable performance evaluations of heuristic methods to be developed in the future.

Our methodological contributions are twofold:

-to propose a new mixed integer programming formulation for this strategic electric locationrouting problem and -to develop a Benders decomposition algorithm embedded in a novel two-phase framework to solve the problem to proven optimality by making use of several lower and upper bounds.

The rest of the paper is organized as follows: Section 2 gives an informative review on the related works. In Section 3, we provide the notation used throughout the paper and present our mathematical formulation. In Section 4, we propose our Benders decomposition algorithm followed by Sections 5 and 6 describing the implementation and the intermediate reduction process, respectively. In Section 7, we provide the setting and present the results of our computational study. We conclude in Section 8 with a summary of our findings and a highlight of some future research directions.

Related work

Location of recharging stations can be seen as a facility location problem. The purpose is then to decide on the optimal number and locations of facilities given the position of customers to serve. In this vein, [START_REF] He | Incorporating institutional and spatial factors in the selection of the optimal locations of public electric vehicle charging facilities: A case study of beijing, china[END_REF] present a case study in Beijing, China. Their objectives are to incorporate the local constraints of supply and demand on public EV charging stations into facility location models, and to compare the optimal locations from three different location models: the set covering model, the maximal covering location model, and the p-median model. [START_REF] Liu | Locating multiple types of charging facilities for battery electric vehicles[END_REF] address the optimal location of multiple types of charging facilities, including dynamic wireless charging facilities and different levels of plug-in charging stations. Their tri-level program first treats the model as a black-box optimization, which is then solved by an efficient approximation model.

However, as raised in [START_REF] Salhi | The effect of ignoring routes when locating depots[END_REF], facility location and routing decisions are interdependent and should be tackled simultaneously. In the general case where both vehicles and depots are capacitated, the problem is known as the capacitated location routing problem (CLRP). The aim here is to i) define which depots must be opened, ii) assign each serviced node (customer) to one and only one depot and, iii) route the vehicle to serve the nodes, in such a way that the sum of the depot cost and the total routing cost is minimized. Many papers appeared in the subject and more particularly during the last decade, as shown in surveys by [START_REF] Nagy | Location-routing: Issues, models and methods[END_REF]; [START_REF] Prodhon | A survey of recent research on location-routing problems[END_REF], and [START_REF] Schneider | A survey of the standard location-routing problem[END_REF]. To solve this NP-hard problem to proven optimality, branch-and-cut [START_REF] Belenguer | A branch-and-cut method for the capacitated location-routing problem[END_REF] and set partitioning based [START_REF] Akca | A branch-and-price algorithm for combined location and routing problems under capacity restrictions[END_REF][START_REF] Contardo | A branch-and-cut-and-price algorithm for the capacitated location-routing problem[END_REF] algorithms are available in the literature. Additionally, several new efficient metaheuristics are proposed. These include a cooperative Lagrangean relaxation-granular tabu search heuristic by [START_REF] Prins | Solving the capacitated location-routing problem by a cooperative Lagrangean relaxation-granular tabu search heuristic[END_REF], an adaptive large-neighborhood search (ALNS) by [START_REF] Hemmelmayr | An adaptive large neighborhood search heuristic for two-echelon vehicle routing problems arising in city logistics[END_REF], and a three-phase matheuristic by [START_REF] Contardo | A grasp+ilp-based metaheuristic for the capacitated location-routing problem[END_REF]. Other studies cover a multiple ant colony optimization algorithm [START_REF] Ting | A multiple ant colony optimization algorithm for the capacitated location routing problem[END_REF]) and a two-phase hybrid heuristic [START_REF] Escobar | A two-phase hybrid heuristic algorithm for the capacitated location-routing problem[END_REF]. Very recently, a tree-based search algorithm by [START_REF] Schneider | Large composite neighborhoods for the capacitated locationrouting problem[END_REF] and a Genetic Algorithm by [START_REF] Lopes | A simple and effective evolutionary algorithm for the capacitated location-routing problem[END_REF] are proposed.

The integration of the location of recharging stations with the routing decision, also called electric location-routing problem (ELRP), is relatively recent though it can lead to massive environmental benefits. To the best of our knowledge, the first study of simultaneous vehicle routing and charging station location for commercial EVs is presented in a conference paper in 2012 by [START_REF] Worley | Simultaneous vehicle routing and charging station siting for commercial electric vehicles[END_REF]. Then, [START_REF] Yang | Battery swap station location-routing problem with capacitated electric vehicles[END_REF] introduce the interesting battery swap station location-routing problem, where the charge is completely fulfilled at each stop. The authors develop two heuristic approaches. The problem is revisited by [START_REF] Hof | Solving the battery swap station location-routing problem with capacitated electric vehicles using an avns algorithm for vehicle-routing problems with intermediate stops[END_REF] who adapt an interesting and powerful adaptive variable neighborhood search (AVNS) heuristic originally dedicated to the vehicle routing problem (VRP) with intermediate depots. Recently, [START_REF] Zhang | A novel location-routing problem in electric vehicle transportation with stochastic demands[END_REF] introduce a battery swap station location-routing problem with stochastic demand and solve this problem by developing a hybrid algorithm combining binary particle swarm optimization and variable neighborhood search. Another relevant study by [START_REF] Yıldız | A branch and price approach for routing and refueling station location model[END_REF] introduce a branch and price algorithm for routing and refueling station location problem.

The first paper dealing with partial recharge may come from [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF], and is dedicated to a Green Vehicle Routing Problem (G-VRP). In G-VRP the fleet is composed of Alternative Fuel Vehicles (AFV) where, in addition to the routing of each EV, the amount of energy recharged and the technology used must also be determined. However, the location aspect is not considered. Constructive and improving heuristics are embedded in a Simulated Annealing framework. The partial recharging policies are then reused showing that they may considerably improve the routing decisions as noted by [START_REF] Keskin | Partial recharge strategies for the electric vehicle routing problem with time windows[END_REF]. Thus, Schiffer and Walther (2017b) extend the problem by including the location of charging stations which leads to the electric location routing problem with time windows and partial recharging (ELRP-TWPR). The authors focus on a problem with a single type of vehicle and multiple visits to the stations. They propose a mathematical formulation based on Miller-Tucker-Zemlin type constraints, supported by several preprocessing steps to eliminate the arcs that violate time windows, capacity, and battery restriction constraints. The Location Routing Problem with Intraroute Facilities which is a generalization of the ELRP-TWPR is explored by Schiffer and Walther (2017a) where large instances are solved using an ALNS which is enhanced by local search and dynamic programming components. A lower bounding procedure integrated to this ALNS algorithm by [START_REF] Schiffer | Designing sustainable mid-haul logistics networks with intra-route multi-resource facilities[END_REF] provides improved results for solving the ELRP-TWPR.

Our problem can be considered as a generalization of the electric vehicle routing problem (EVRP) with location decisions or an electric location-routing problem (ELRP) with a heterogeneous fleet, multi-type stations, multi-visit, and partial recharging. The EVRP literature is not extensively discussed here except those papers considering the location decisions. However, we refer the reader to [START_REF] Pelletier | 50th anniversary invited article-goods distribution with electric vehicles: review and research perspectives[END_REF] for an overview on the EVRP studies. In the next section we provide the notation and a mathematical formulation of the problem.

Notation and Problem Formulation (PF)

Consider a given network with a set of customer locations and a set of potential charging station locations, from which we are required to select a subset of stations. Each customer should be served by a vehicle originating from the depot and each vehicle can perform a single trip. The vehicles have a battery restriction and they have to visit one or more among the selected charging stations before the battery is depleted if a trip longer than their range is to be traversed. In addition, we have a fleet of heterogeneous vehicles with a limited number of vehicles of each type. We first provide the notation and a scheme for allowing multiple visits which is then followed by the new formulation.

Notation

In this section, we list the parameters and the decision variables as follows:

Parameters: G = (N, A): the given network.

A: the arc set. Decision variables:

y j = 1 if station j ∈ J is open, 0 otherwise x k ij = 1 if arc (i, j) is traversed by vehicle k ∈ K, 0 otherwise z k j is the amount of energy recharged at station j ∈ J for vehicle k ∈ K b k ij is the battery level of vehicle k ∈ K at node i ∈ N before leaving for node j ∈ N expressed in kWh l k ij is the cumulative load of vehicle k ∈ K at node i ∈ N before leaving for node j ∈ N .
In the remaining of this paper, we assume I ⊂ J but all the methods can be easily adapted to the case where I \ J = ∅ by simply defining y j and z k j variables for all j ∈ I ∪ J, k ∈ K and setting y j = z k j = 0, ∀j ∈ I \ J.

A novel mechanism that caters for multiple visits

In order to allow multiple visits to a station, we perform the following interesting and powerful modification on our input network. For each potential station, we create as many dummy copies as the number of customers (Steps 2-3 in Algorithm 1). If this potential station is also a demand node, the demand of the first copy is identical to the demand of the potential station whereas the demand of the remaining copies is set to zero. Similarly, the fixed cost of the first copy is identical to the fixed cost of the station whereas it is set to zero for the remaining copies (Step 4). The arc set of the modified network includes all the arcs of the original network. Additionally, all distinct node pairs in the modified network are connected via a direct arc except the pairs which are the copies of the same station (Step 5). Finally, we define an arc set for each vehicle which contains all arcs in the modified network whose energy consumption is not greater than the range of the vehicle and the total demand of its endpoints is not greater than the freight capacity of the vehicle (Step 6).

Algorithm 1: Network modification

Step 1: Let |I| be the number of demand nodes.

Step 2: Create |I| copies of each station.

Step 3: Form set J A j = {j 1 , j 2 , . . . , j |I| } for each j ∈ J and J A = j∈J J A j .

Step 4: For each j ∈ J, set

f j 1 = f j ; d j 1 = d j and f j i = d j i = 0, i = 2, . . . , |I| where j 1 , . . . , j |I| ∈ J A j . Step 5: Let N E = J A ∪ {0} and A E = A ∪ {(i, j) : i, j ∈ N E ; i = j; ¬(i, j ∈ J A l for some l ∈ J)}.
Step 6: Define

A k = {(i, j) ∈ A E : e ij ≤ β k ; d i + d j ≤ Q k } for k ∈ K.

Mathematical formulation PF

The problem can be formulated with commodity flows as follows:

min k∈K (i,j)∈A k c ij x k ij + k∈K j∈J A r k z k j + j∈J A f j y j + k∈K (0,i)∈A k v k x k 0i (1) s.t. y i ≤ y j , i ∈ J A j : i = j (2) i∈J A x k 0i ≤ 1, k ∈ K (3) k∈K (j,i)∈A k x k ji = 1, i ∈ J A : d i > 0 (4) (j,i)∈A k x k ji ≤ y i , i ∈ J A : d i = 0, k ∈ K (5) (i,j)∈A k x k ij - (j,i)∈A k x k ji = 0, i ∈ J A , k ∈ K (6) (i,j)∈A k (l k ij -d i x k ij ) = (j,i)∈A k l k ji , i ∈ J A , ∀k ∈ K (7) j∈J A l k 0j = 0, k ∈ K (8) l k ij ≤ Q k x k ij , k ∈ K, (i, j) ∈ A k (9) (i,j)∈A k e ij x k ij - j∈J A z k j ≤ β k , k ∈ K (10) (i,j)∈A k b k ij = (j,i)∈A k (b k ji -e ji x k ji ) + z k i , i ∈ I ∪ J A , k ∈ K (11) z k j ≤ β k y j , j ∈ J A , k ∈ K (12) z k j ≤ β k (i,j)∈A k x k ij , j ∈ J A , k ∈ K (13) b k 0j = β k x k 0j , k ∈ K, (0, j) ∈ A k (14) b k ij ≤ β k x k ij , k ∈ K, (i, j) ∈ A k (15) b k ij ≥ e ij x k ij , k ∈ K, (i, j) ∈ A k (16) y j ∈ {0, 1}, j ∈ J A (17) x k ij ∈ {0, 1}, k ∈ K, (i, j) ∈ A k (18) l k ij ≥ 0, k ∈ K, (i, j) ∈ A k (19) b k ij ≥ 0, k ∈ K, (i, j) ∈ A k (20) z k j ≥ 0, j ∈ J A , k ∈ K. (21)
The objective function (1) minimizes the total sum of routing costs, charging costs, fixed costs of opening stations, and fixed cost of using vehicles. If a zero-demand copy of a station is opened, Constraints (2) force the original copy of this node to be opened and therefore, ensure that the costs of the stations are counted in the objective function. By Constraints (3), we restrict the number of trips by each vehicle to at most one. Constraints (3)-( 9) together ensure that each client is served by a unique vehicle that starts its trip at the depot and the capacities of vehicles are respected. Constraints ( 5) ensure that a zero-demand copy of any station is visited only if that station is open. We ensure the elimination of sub-tours for each vehicle trip via the load (flow) preservation constraints ( 7)-( 9). Note that Constraints (9) also ensure that the vehicle freight capacities are always respected. Battery restriction on the vehicles are imposed by Constraints ( 10) and ( 11). Constraints ( 12) and ( 13) avoid recharging of a vehicle at a node that has no station and that is not visited by that vehicle, respectively.

We initialize the battery level for each vehicle to 100% by Constraints ( 14). For each arcvehicle pair, Constraints (15) restrict the amount of battery level with full battery level if the arc is traversed by the vehicle and set it to zero otherwise; Constraints (16) make sure that the battery level is larger than the energy consumption on the arc that will be traversed by the vehicle. Finally, Constraints ( 17)-( 21) represent the binary and non-negativity restrictions on the decision variables.

Benders Decomposition Algorithm (BDA)

Our mathematical formulation can be solved by using a Benders decomposition [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF]) framework that we briefly described here. The details of our algorithm are presented next. The classical Benders decomposition method aims to solve a mixed integer program (MIP) with a group of integer variables and a group of continuous variables by decomposing the MIP into a master problem (MP) with all integer variables and a series of subproblems with continuous variables. For each feasible solution of MP, a subproblem (SP) is constructed by fixing the values of all the integer variables in the MIP to the value obtained from the master problem. Each extreme ray and extreme point of the dual of this SP provides a so called feasibility and an optimality cut, respectively, for the MP. Since the full enumeration of the extreme points and extreme rays is impractical, the cutting plane procedures are usually employed for the generation and the addition of these cuts.

The classical Benders decomposition method could suffer from slow convergence especially if the subproblem is large in size. On the other hand, the method would perform relatively efficiently if the subproblem can be decomposed further into smaller and easy-to-solve subproblems as in multi-commodity, multi-period, or multi-scenario problems [START_REF] Birge | Introduction to stochastic programming[END_REF]. Motivated by this fact, we aim to further decompose our problem into |K| smaller problems, each one corresponding to a single vehicle trip. For this purpose, we decide to keep y, x, l variables in the master problem and z, b variables in the subproblems. The separation method for the optimality cuts plays a crucial role in efficient Benders implementations. To speed up our implementation, we adopt the recently developed high-performing strategy of C ¸alık and Fortz (2019) and modify our formulation accordingly to obtain only feasibility cuts from the dual subproblems. In order to achieve this, we introduce an additional non-negative decision variable w k , ∀k ∈ K and make a slight modification to our model to ensure that w k takes value j∈J A z k j , ∀k ∈ K. The modified formulation (PF2), as given below, is defined by Constraints ( 2)-( 21) and ( 22)-( 24):

(P F 2) min k∈K i∈N E j∈N E :i =j c ij x k ij + k∈K r k w k + j∈J A f j y j + k∈K i∈J A v k x k 0i (22) s.t. w k = j∈J A z k j , ∀k ∈ K (23) w k ≥ 0, ∀k ∈ K (24) (2) -(21).
When solving PF2 in a Benders fashion, we employ a branch-and-cut framework which keeps y, x, l, w variables in the master problem (MP) and z, b variables in the subproblems. 2) -( 9), ( 17) -( 19), ( 24)

(M P ) min (22) s.t. (
w k ≥ i∈N E j∈N E :j =i e ij x k ij -β k , ∀k ∈ K. ( 25
)
Let (y, x, l, w) be the vector of variable values in the solution obtained from the master problem. One can easily observe that if w k = 0, then, no recharging is needed for the corresponding vehicle trip and (y, x, l, w) is feasible for PF2. On the other hand, if w k > 0, we construct and solve the dual of the subproblem SP k (y, x, l, w) for every k ∈ K. Note that when w k > 0, an optimal solution to the original problem should satisfy Equation (32), which is helpful in the following mathematical manipulations leading to an efficient implementation.

SP k (y, x, l, w) min 0 (26) s.t. j∈J A z k j = w k , (27)
z k j ≤ β k y j , ∀j ∈ J A (28) z k j ≤ β k i∈N E :i =j x k ij , ∀j ∈ J A (29) j∈N E :j =i b k ij = j∈N E :j =i (b k ji -e ji x k ji ) + z k i , ∀i ∈ J A (30) b k 0j = β k x k 0j , (31)
b k j0 = e jo x k j0 , ( 32
) b k ij ≤ β k x k ij , i, j ∈ N E : i = j (33) b k ij ≥ e ij x k ij , i, j ∈ N E : i = j (34) b k ij ≥ 0, i, j ∈ N E : i = j (35) z k j ≥ 0, ∀j ∈ J (36)
Note that Constraints (10) of PF is ensured by Constraints ( 25) and ( 27) as 28)-( 36) is non-empty and w k > 0, then

j∈J A z k j = w k ≥ i∈N E j∈N E :j =i e ij x k ij -β k , ∀k ∈ K. Moreover,
j∈J A z k j = i∈N E j∈N E :j =i e ij x k ij -β k . Proof. z k i = j∈N E :j =i b k ij - j∈N E :j =i (b k ji -e ji x k ji ), ∀i ∈ J A by (30). Moreover, i∈J A β k x k 0i = β k since w k > 0. i∈J A z k i = i∈J A j∈N E :j =i b k ij - i∈J A j∈N E :j =i (b k ji -e ji x k ji ) = i∈J A b k i0 - i∈J A b k 0i + i∈J A j∈J A :i =j (b k ij -b k ji ) + i∈J A j∈N E :j =i e ji x k ji = i∈J A e k i0 x k i0 - i∈J A β k x k 0i + i∈J A j∈N E :j =i e ji x k ji = i∈N E j∈N E :j =i e ji x k ji -β k
After elimination of equality constraints and a few mathematical manipulations on the remaining subproblem, we obtain the following SP k in canonical maximization form for each k ∈ K: max 0 (37)

s.t. - j∈J z k j ≤ -w k , ( 38
)
z k j ≤ β k y j , ∀j ∈ J A (39) z k j ≤ β k i∈J A :i =j x k ij , ∀j ∈ J A (40) z k j + i∈J A :i =j b k ij - i∈J A :i =j b k ji ≤ e k j0 x k j0 -β k x k 0j + i∈N E :i =j e ij x k ij , ∀j ∈ J A (41) -z k j - i∈J A :i =j b k ij + i∈J A :i =j b k ji ≤ β k x k 0j -e j0 x k j0 - i∈N E :i =j e ij x k ij , ∀j ∈ J A (42) b k ij ≤ β k x k ij , ∀i, j ∈ J A : i = j (43) -b k ij ≤ -e ij x k ij , ∀i, j ∈ J A : i = j (44) b k ij ≥ 0, ∀i, j ∈ J A : i = j (45) z k j ≥ 0, ∀j ∈ J A (46)
Let α, δ j , π j , γ j ,ρ j , φ ij , ij be the dual variables associated with constraints (38)-( 44), respectively. Then, we can write the equivalent dual problem D k (y, x, l, w) for each k ∈ K as follows:

D k (y, x, l, w) min -w k α + j∈J A β k y j δ j + i∈N E j∈J A :i =j β k x k ij π j + j∈J A e j0 x k j0 γ j - j∈J A β k x k 0j γ j + i∈N E j∈J A :i =j e ij x k ij γ j + j∈J A β k x k 0j ρ j - j∈J A e j0 x k j0 ρ j - i∈N E j∈J A :i =j e ij x k ij ρ j + i∈J A j∈J A :i =j β k x k ij φ ij - i∈J A j∈J A :i =j e ij x k ij ij (47) s.t. -α + δ j + π j + γ j -ρ j ≥ 0, ∀j ∈ J A (48)
-

γ i + γ j + ρ i -ρ j + φ ij -ij ≥ 0, ∀i, j ∈ J A : i = j (49) α ≥ 0, ( 50
)
δ j , γ j , π j , ρ j ≥ 0, ∀j ∈ J A (51) φ ij , ij ≥ 0, ∀i, j ∈ J A : i = j (52)
In order to avoid solving the same dual problem twice (once for detecting unboundedness and once for obtaining a feasibility cut), we solve a bounded dual problem instead. This will imply unboundedness of D k (y, x, l, w) if the optimal value is negative (see Lemma 4.2). To do so, we bound variables α, γ j , ρ j , ∀j ∈ J A , and ij , ∀i, j ∈ J A : i = j by 1 from above. Let us refer to this bounded dual problem as D B k (y, x, l, w). If the optimal value of D B k (y, x, l, w) is negative valued, we add the feasibility cut (53) to MP to cut the current solution (y, x, l, w).

-αw

+ j∈J A β k δ j y j + i∈N E j∈J A :i =j β k π j x k ij + j∈J A e j0 γ j x k j0 - j∈J A β k γ j x k 0j + i∈N E j∈J A :i =j e ij γ j x k ij + j∈J A β k ρ j x k 0j - j∈J A e j0 ρ j x k j0 - i∈N E j∈J A :i =j e ij ρ j x k ij + i∈J A j∈J A :i =j β k φ ij x k ij - i∈J A j∈J A :i =j e ij ij x k ij ≥ 0 ( 53
)
Lemma 4.2. If the optimal value of D B k (y, x, l, w) is negative, then D k (y, x, l, w) is unbounded.

Proof. Let ψ be an optimal solution to D B k (y, x, l, w) and let g(ψ) < 0 be the value of this solution. For any positive constant υ, υψ is a feasible solution for D k (y, x, l, w). Then, for an arbitrarily large υ, g(υψ) = υg(ψ) < 0 will be an arbitrarily small solution value for D k (y, x, l, w) which implies that D k (y, x, l, w) is unbounded.

Implementation Details -General Framework

The general framework of our algorithm mainly consists of two phases. In Phase I, we solve the problem with at most one visit to each station. This is done by including only one copy of each station in BDA models MP and D k (), k ∈ K. We refer to the BDA solving this restricted problem as BDA 1 . In the second phase, we focus on the general problem that allows multiple visits to stations. Between the two phases, we perform an intermediate reduction procedure (See Section 6) to decrease the size of the problem in Phase II. The aim is to cut as much as possible without eliminating any potential solution that is better than the one in Phase I. We provide a brief summary of the general framework in Section 5.2. For clarity of presentation, we use the notation 'BDA' throughout the paper to refer to both the algorithm and the formulation.

Through our preliminary experiments, we observed that our algorithm has a better convergence behavior if we introduce a high quality initial feasible solution to our master problem. In order to achieve this, we first perform a 'Step 0' process where we solve our BDA formulation via a CPLEX option that allows stopping after finding the first integer feasible solution. We also introduce a partial warm start solution to CPLEX by opening all potential stations. In our experiments, CPLEX usually finds a solution with all stations opened. We then improve this solution by closing some of the stations. This removal process is a greedy approach based on checking the energy consumption between three consecutive stations and then closing the intermediate one if the battery level is sufficient to go from the first one to the third one. Finally, we introduce the set of open stations of this improved solution as a partial warm start solution for our Phase I problem and solve BDA 1 with the valid inequalities given next in Section 5.1.

Valid Inequalities for Phase I

Let N V min be a lower bound on the number of vehicles needed for any feasible solution. We can obtain such a lower bound by solving a bin packing problem (BPP) as follows. Define s k = 1 if vehicle k is used, 0 otherwise and a ik = 1 if the request of customer i is provided by vehicle k, otherwise. Constraints (55) assign each customer to a vehicle while Constraints (56) ensure that these assignments respect the capacities of vehicles.

(BP P ) N V min = min k∈K s k (54) s.t. k∈K a ik = 1, ∀i ∈ I (55) i∈I d i a ik ≤ Q k s k , ∀k ∈ K (56) s k ∈ {0, 1}, ∀k ∈ K (57) a ik ∈ {0, 1}, ∀i ∈ N, k ∈ K. ( 58
)
We can detect the infeasibility due to insufficient freight capacity by solving BPP. Our preliminary experiments revealed that introducing Constraint (59), which enforces using at least N V min vehicles, usually reduces the solving time. This observation has led us to include this constraint in our computations for every model of Phase I and Phase II.

k∈K j∈J A x k 0j ≥ N V min (59)
When we solve BDA to optimality with at most one visit to each station (BDA 1 ), we include the following sets of valid inequalities to our master problem:

k∈K j:(i,j)∈A k x k ij ≤ 1, i ∈ J A : d i > 0 (60) j:(i,j)∈A k x k ij - j:(0,j)∈A k x k 0j ≤ 0, i ∈ N E , ∀k ∈ K (61) i∈J A x k i0 ≤ 1, ∀k ∈ K (62) j:(j,i)∈A k x k ij ≤ y i , i ∈ J A : d i = 0, ∀k ∈ K (63) y i ≤ k∈K j∈N E x ij , ∀i ∈ J A : d i = 0 (64) y j ≤ k∈K i∈N E x ij , ∀j ∈ J A : d j = 0 (65) w k ≤ j∈J A β k y j ∀k ∈ K. (66)
Constraints ( 60) restrict the number of arcs entering a demand node to at most one. Constraints (61) ensure that an arc is visited by a vehicle only if that vehicle leaves the depot. Constraints (62) make sure that each vehicle enters the depot at most once. Constraints (63) forbid leaving a zero-demand copy of a station if it is not open while Constraints ( 64) and ( 65) forbid opening zero-demand station copies if they are not visited by any vehicle. Constraints (66) limit the total recharging for each vehicle by full battery charging times the number of open stations.

Even though most of these constraints are implied by the original constraints, their inclusion improves the time performance of our algorithm considerably.

Let (y 1 , x 1 , l 1 , w 1 ) be the solution with value Z 1 that we obtain from Phase I. After the intermediate process which will be explained in Section 6, we proceed to Phase II to solve a reduced problem via BDA with the valid inequalities of Section 5.2 below. We introduce y 1 as a partial warm start solution to the Phase II problem.

Valid Inequalities for Phase II

When we apply BDA for the last time with all possible copies of potential stations, in addition to the valid inequalities ( 59),( 61)-( 66), we also introduce the following set of valid inequalities to break the symmetry between the copies of stations:

i:(i,j)∈A k x k ij ≤ i:((j-1),i)∈A k
x k (j-1)i , ∀k ∈ K, j is the m th copy of some j 1 :

d j 1 > 0, m ≥ 3 (67) i:(i,j)∈A k x k ij ≤ i:((j-1),i)∈A k
x k (j-1)i , ∀k ∈ K, j is the m th copy of some j 1 :

d j 1 = 0, m ≥ 2. ( 68
)
Constraints ( 67) and ( 68) make sure that an additional copy of any station is visited by a vehicle only if the preceding copy is visited by the same vehicle. Exceptionally, the second copy (the first non-original copy), might be visited by a vehicle not serving the original copy if it is a demand node.

General framework

Below we give a brief summary of the general framework of our algorithm:

Step 0: Solve BDA 1 to obtain a feasible solution (y 0 , x 0 , l 0 , w 0 ) (not necessarily optimal).

Close the redundant stations of (y 0 , x 0 , l 0 , w 0 ) in a greedy manner and obtain (y, x, l, w).

Step 1: Phase I: Solve BDA 1 ∪( 59)-( 66) with partial warm start y to obtain the optimal solution (y 1 , x 1 , l 1 , w 1 ).

Step 2: Apply the intermediate process (see Section 6) to reduce the size of BDA∪( 59)-(68).

Step 3: Phase II: Solve the reduced BDA∪( 59)-( 68) with partial warm start y 1 to obtain the optimal solution.

Intermediate Reduction Process

Creating multiple copies of stations leads to a large-size formulation and excessive solving times. We develop a two phase method that solves our Benders formulation initially for a single copy of each station. Based on the value Z 1 of the solution obtained at this stage, we apply an intermediate processing procedure that checks the availability of a solution with multiple copies of stations that has a smaller objective value than Z 1 . This is an iterative procedure that proceeds by increasing the number of copies considered, say m, one by one and applies lower bound checking steps. The aim of this procedure is to check whether there exists a solution of BDA with exactly m copies for some station j whose cost is lower than Z 1 . Lemma 6.1. Let Z LB (m,j,k) be a lower bound on the cost when exactly m copies of station j is visited by vehicle k. If Z LB (m,j,k) ≥ Z 1 , ∀j, k, then, there exists no solution with m copies of any station whose value is less than Z 1 .

Proof. Any feasible solution to a minimization problem provides an upper bound. Therefore, the value of any feasible solution as described in Lemma 6.1 has to be greater than or equal to Z LB (m,j,k) ≥ Z 1 . Lemma 6.2. If there exists some lower bound Z LB (m,j,k) such that Z LB (m ,j,k) ≥ Z LB (m,j,k) , ∀m ≥ m and if Lemma 6.1 holds for such Z LB (m,j,k) of m for all j, k, then, there exist no solution with more than or equal m copies of any station whose value is less than Z

1 . Proof. Z LB (m ,j,k) ≥ Z LB (m,j,k) ≥ Z 1 , ∀j, k, m
such that m ≥ m by Lemma 6.1. Then, there exists no solution with m copies of any station whose value is less than Z 1 .

Below we give the details on how we obtain a lower bound that satisfies Lemma 6.2. Let us consider a potential station j. If we use exactly m copies of this station, it means that we visit at least m -1 different customers with some vehicle k. This leads to a partial network structure as 0 . . . j . . . i 1 . . . j . . . i 2 . . . . . . i m-1 . . . j . . . 0.

Let E m and R m be the amount of energy consumption and the amount of recharging needed, respectively, when we visit station j exactly m times by some vehicle k 1 . Now, we consider two cases:

Case 1: all customers are visited by k 1 .

Case 2: some customers are visited by other vehicle(s). For any of Case 1 or Case 2, the following observation holds:

Observation:

(i) E m > E base m = (m -1)β k 1 + e j0 and R m > R base m = (m -2)β k 1 + e j0 if m is even. (ii) E m > E base m = (m -1)β k 1 and R m > R base m = (m -2)β k 1 if m is odd.
In Figure 1, we illustrate this observation for m = 2 and m = 3. In this figure, if E 2 ≤ β k 1 + e j0 , we would not need to visit j twice. Similarly, if E 3 ≤ 2β k 1 , it would be redundant to visit j three times. In Case 1, we calculate the minimal possible energy consumption on such a partial network. More explicitly, we define E C m = e 0j + e ji 1 + e i 1 j + e ji 2 + e i 2 j + . . . + e ji m-1 + e j0 where i 1 , . . . , i m-1 are m -1 closest customers to j (e.g. see Figure 2). For this particular case, we can further obtain a lower bound on the total energy consumption by constructing a 1-tree obtained via a minimum spanning tree (MST) which spans the union set of all customers and m copies of j and that is connected to the depot node with two minimal arcs. Let the energy consumption on this 1-tree be

E T m and E m = max{E C m , E T m , E base m }.
We obtain a lower bound C on the total routing cost similarly. Then, R = max{E m -β k 1 , 0, R m } gives us the amount of recharging needed for this partial network and Z LB = R × r k 1 + C + f j + v k 1 gives us a lower bound on the cost of routing all customers by vehicle k via visiting j m times or more.

When we look at Case 2, we investigate all possible vehicle combinations that need to be considered by iteratively increasing the number of additional vehicles. If we find a combination with k vehicles whose lower bound is less than Z 1 , we do not check the combinations with more than k vehicles. Let us assume that in addition to k 1 , we use K * = {k 2 , . . . , k l }. This means that we are visiting a different customer by each additional vehicle. Therefore, we add 0 . . . i k h . . . 0 as a connected component to our partial network for each vehicle k h ∈ K * where i k h is the closest customer to depot which is not served by preceding vehicles. Let E(k h ) and R(k h ) denote the total energy consumption and the amount of recharging needed, respectively, for the connected component for vehicle k h ∈ K * . In a similar fashion to that of Case 1, we calculate the total energy consumption

E C m (K * ) = E m + k h ∈K * E(k h ), the amount of recharging needed R C m (K * ) = R + k h ∈K * R(k h )
, the total routing cost and hence a lower bound Z LB1 K * on the total cost with the corresponding vehicle combination.

In order to obtain another lower bound Z LB2 K * from the 1-tree constructed with the total energy consumption E T m (K * ), this time, we use β = k∈K * ∪k 1 β k as the total battery available in our calculation for the amount of recharging needed, that is, R m (K * ) = max{max{E C m (K * ), E T m (K * )}β, 0} and r = min k∈K * ∪k 1 r k as the unit recharging cost. Our bound Z LB K * for the corresponding combination is defined as 3,4, and 5 for sample 1-tree constructions for visiting two copies of j with 1, 2, and 3 vehicles, respectively. When k vehicles are used, additional k -1 copies of the depot are created and added to the set of nodes to find an MST. The 1-tree is then constructed by connecting this MST to the remaining (original) copy of the depot via two minimal arcs. Note that there do not exist any arcs between the copies of the same stations or the depot. The 1-trees constructed in this manner provide a lower bound for the shortest-length Hamiltonian cycles as in Figure 6, which also provide a lower bound on the lengths of corresponding VRP tours.

Z LB K * = max{Z LB1 K * , Z LB2 K * }. See Figures
Figure 3: Construction of the 1-tree for obtaining LB2 when visiting m = 2 copies of j with one vehicle.

If min{Z LB , Z LB K * } ≥ Z m-1 for every (k 1 , K * ) combination, then, the value of any solution visiting j no less than m times will be no better than Z m-1 . So, in further iterations, we do not need more than m -1 copies of j. If this holds for all stations, we can terminate the iterative checking procedure and solve our algorithm BDA with at most m -1 copies for each station.

Additional speed-up mechanism:

When solving BDA for this last time, we further apply variable fixing by using the information we obtained from this iterative procedure. More explicitly, if it is decided that we do not need more than l copies at a given station j, we fix all y values to zero for all those copies of j. Similarly, if we decided that visiting more than l copies of station j with vehicle k 1 is not optimal, then we set all x variables of those copies to zero for vehicle k 1 .

Figure 4: Construction of the 1-tree for obtaining LB2 when visiting m = 2 copies of j with two vehicles.

Figure 5: Construction of the 1-tree for obtaining LB2 when visiting m = 2 copies of j with three vehicles.

Figure 6: Illustrative Hamiltonian cycles for calculating the 1-tree lower bound LB2 on the cost of visiting m = 2 copies of j with 1,2, or 3 vehicles 7 Computational Study

Experimental setting and data instances

In order to test our methods, we generated problem instances based on the data set provided by [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF] for EVRP. This data set has 36 different instances with 5, 10, and 15 customers (12 instances for each customer size). The number of potential stations vary between 7 and 22. From these instances, we retrieved the demand and network information (node coordinates). In the remainder of the paper, we will refer to these instances as 'networks' for clarity reasons. The vehicle freight capacities are equal to 200 in the original data. We introduced additional levels of capacities (80, 100), especially, to test relatively smaller instances. Similarly, for the battery capacities, we conducted tests for low, medium, and high capacities (10, 16, 22 kWh) to avoid extremely loose values on the tests of small problems [START_REF] Sassi | Electric vehicle scheduling and optimal charging problem: complexity, exact and heuristic approaches[END_REF].

In our experiments, we use IBM ILOG CPLEX 12.8 in a Java environment. We run our tests on a server with Intel(R) Xeon(R) CPU E5-2640 v3 at 2.60 GHz processor and 16 cores. For each experiment, we set a memory limit of 16 GB and a time limit of 3600 seconds for instances with |I| = 5 and 10800 seconds for the others. BDA and PF run using a single thread.

We assume that the system will be equipped with fast charging facilities. Note that as there are no time windows or maximum travelling restrictions, introducing slow and fast charging facilities in each potential station would lead to optimal solutions with only cheaper type of charging facilities. As we do not have data on location-dependent fixed station costs for the instance networks, we utilize a single type of stations in our experiments. The lifetime of a charging facility is estimated to be 3 years and it is 5 years for the vehicles. When we calculate the fixed costs of opening stations and purchasing/leasing vehicles, we divide their costs by the number of days within their lifetime. We approximately obtain f j = 8 e as the fixed cost of opening stations and v k =16, 26, and 36 e as the fixed cost of low, medium, and high capacity vehicles, respectively. Let l ij be the distance between nodes i, j ∈ N ; then, we set c ij = l ij × 0.03 (cents/km), r k = r = 0.07 (cents/km) for k ∈ K, e ij = l ij × 135 (Wh/km). In order to better tackle the precision issues of CPLEX, we multiplied all the cost values by 100. We conduct experiments on instances with heterogeneous fleets of 2, 3, and 4 vehicles of three different types, see Table 1. In order to observe the value of using heterogeneous fleets compared to homogeneous fleets, we also conduct experiments with homogeneous fleets of 1, 2, 3, and 4 vehicles for each vehicle type. For each network, we test the problem with heterogeneous fleets shown in Table 2 as well as the homogeneous fleets generated. Certain network-fleet combinations are infeasible and we exclude them from our computational analysis.

-2 2 1,2 K4V1-1-1-2 4 1,1,1,2 K2V1-3 2 1,3 K4V1-1-1-3 4 1,1,1,3 K2V2-3 2 2,3 K4V1-1-2-2 4 1,1,2,2 K3V1-1-2 3 1,1,2 K4V1-2-2-2 4 1,2,2,2 K3V1-1-3 3 1,1,3 K4V1-1-3-3 4 1,1,3,3 K3V1-2-2 3 1,2,2 K4V1-1-2-3 4 1,1,2,3 K3V1-2-3 3 1,2,3 K4V1-2-2-3 4 1,2,2,3 K3V1-3-3 3 1,3,3 K4V1-2-3-3 4 1,2,3,3 K3V2-2-3 3 2,3,3

Analysis

We perform two types of analysis: (1) performance analysis of the model and the algorithm in Section 7.2.1 and (2) managerial insights for selecting the fleet types in Section 7.2.2. For simplicity, we provide average results for each network type in Table 3 but the detailed results can be found in the supplementary document. For illustration, we also display the results in Figures 7891011.

In these figures and tables, N S and N V show the number of stations opened and the number of vehicles used, respectively. The value 'g' represents the gap provided by CPLEX at the end of the time limit (0.00 if the problem is solved to proven optimality) for the corresponding model solved. Similarly, 't(s)' represents the total time spent in seconds for the corresponding model or algorithm if it includes additional processes. The value of the best solution obtained from a model is given under 'Obj'. The selected vehicle types are indicated with their Q k values. For example, [80,100] indicates that the corresponding solution selects one vehicle of type 1 and one vehicle of type 2.

Methodological observations and insights

Among 204 instances of |I| = 5 (|N | ∈ {8, 9}), our formulation PF can solve all but ten instances to optimality within one hour. The average time spent by PF on these instances is 523.93 seconds, which is much larger compared to that of BDA. On the other hand, we observe that our algorithm BDA can solve all the instances with |I| = 5 to optimality requiring a few seconds only. In fact, BDA is relatively faster than PF in every single instance.

For larger problems (|I| = 10, 15), optimality could not be guaranteed with PF within the time limit. In fact, PF could not even find a feasible solution for many instances. Based on these observations, for larger problems, we present the results for BDA only.

For the networks with 10 customers and up to 14 potential stations (|N | ∈ {13, 15}), BDA solves all Phase I instances to optimality in 188.53 seconds on average. The algorithm reaches either the time or memory limit before proven optimality in Phase II for 6 out of 200 instances in this category. These are instances either with fleet type K4V1-1-1-2 or K4V1-1-2-2, which are challenging also for larger instances with 15 customers, as observed in Figure 7. Figure 7 compares the average solving times (Phase I, Phase II, and Phase I+Phase II) and gaps across fleet types. In this figure, Phase I gaps and times are better indicators for challenging instances as Phase II gaps and solving times are not available due to memory limit for some instances. BDA is able to solve 65% of the instances with |I| = 15 (|N | ∈ {18, 19, 20, 21, 22, 23}) to proven optimality within the A common pattern we observe with the instances reaching the memory limit is that (1) their fleets consist of four vehicles with either one type 3 vehicle and three type 1 vehicles or only type 1 and type 2 vehicles and (2) Phase I solution uses all four vehicles. The intermediate processing procedure cannot reduce the Phase II problem sufficiently due to high level of degeneracy and low-cost and low-capacity vehicles in the fleet.

Overall, the algorithm is able to provide very high quality solutions within three hours. It solves 88% of all instances to optimality. The average gap among all instances is 0.03 only.

For several instances, Phase I and Phase II costs are not identical. We interpret this as follows:

• If Phase II problem does not reach proven optimality within the time limit and Phase I cost is less than Phase II cost, this indicates the benefit of solving a restricted version when the original problem is too difficult to solve within the available time and computational resources.

• If neither Phase I nor Phase II problems reach proven optimality within the time limit and Phase I cost is less than Phase II cost, either a solution using multiple copies of a station is found in Phase II, which is not feasible for Phase I problem, or Phase I solution helps Phase II problem in finding a better-quality solution faster.

Below are some further observations and future research directions for developing algorithms with improved performance.

• The algorithm needs more time to reach proven optimality as the size of the fleet increases. This is often because the dual bound is too weak and the majority of the time is spent for closing the gap. The dual bounds can be strengthened by using good valid inequalities. However, valid inequalities might also make the formulation heavier and more difficult to find feasible solutions. Therefore, it is in general more efficient to decompose the problem into smaller problems with smaller fleet configurations and solve them iteratively by updating the fleet configurations at each step. Obviously, there will be a trade-off between the number of small fleet configurations and the size of each configuration as in most decomposition methods.

• The focus in this paper was on heterogeneous-fleet problems with a limited number of vehicles.

For solving instances with homogeneous fleets, the efficiency of the algorithm can be improved by updating the intermediate process. A similar improvement-procedure update would also be helpful in solving the instances with an unlimited number of vehicles of each type. It is, indeed, interesting to note that small vehicles result in higher costs, an observation which is not that obvious. This is because they require potentially more charging and more stations to be opened which incur very high costs.

Managerial insights

Figure 9 shows the results of an analysis from another perspective where we calculate the averages over all fleet types for each network. In this figure, we can clearly observe that the number of stations needed changes a lot depending on the network type. This is then reflected in the cost. We also observe that the average cost for a smaller-network instance, for example, rc108c5, can be much higher than the cost for a larger-network instance, for example, c106c15. Moreover, we compare the average cost, number of opened stations and vehicles used with homogeneous and heterogeneous fleets of at least two vehicles in Figure 10. We observe that the average cost is much higher when using homogeneous fleets compared to heterogeneous ones, it is indeed twice as much for instances with |I| = 5. The average number of stations opened with homogeneous fleets is also larger for each network group. And in general, a similar conclusion can be made for the average number of vehicles used.

In Figure 11, we also show the average cost per homogeneous fleet, including the fleets with a single vehicle. Similar to previous observations, the average cost is lower when the fleet contains larger vehicles.

Below we provide some further key observations:

• For the majority of the instances where |I| = 5 and all types of vehicles are present in the heterogeneous fleet, medium-size vehicles (type 2) are not selected in the optimal solutions. This is not the case with larger networks.

• When the fleets are homogeneous, no instance uses four vehicles of type 3 and only four instances, two from networks rc202c15 and rc204c15 each, use three type 3 vehicles.

• When solving the instances with fleets of at least two vehicles, we observe that the optimal solutions serving all the demand via a single vehicle only uses the largest vehicle type (type 3). Although there exist several instances where it is feasible to serve all demand with a single vehicle of type 1 or type 2, such solutions are suboptimal and solutions with lower cost can be obtained using multiple vehicles.

Conclusion and Future Research Directions

In this paper, we introduce an electric location-routing problem with heterogeneous fleet and partial recharging. We initially propose a new mixed integer programming formulation for this problem. This is a formulation with three-index binary routing variables where the sub-tour elimination is enhanced via a group of load (flow) preservation constraints. We further utilize additional nonnegative continuous variables to satisfy battery restrictions and energy-related constraints. We test our formulation on small problem instances from the literature. Although the formulation is able to solve instances with 5 customers and up to 8 potential stations to optimality, we observe that its performance is limited when it comes to solving larger problems.

As we aim to solve this problem to optimality, we further develop a two-phase algorithm based on the Benders decomposition of our formulation. The first phase solves a restricted version of the problem that allows at most one visit to each station. By using the information obtained, the second phase problem, which is the generalized problem allowing multiple visits to any station, is reduced in size, making it relatively easier compared to the case with no a priori processing. This enhancement step allows us to solve 88% of all the instances with up to 15 customers and 22 potential stations to optimality. The average optimality gap over all other instances is negligible, just 0.03. In summary, our approach obtains very high quality solutions within the time limit.

We observe through our experimental study that the problem is usually harder to solve when the vehicle capacities are smaller. We also found that using small vehicles results in higher costs.

Though the main focus of this study is to present an exact method with proven optimality, this approach can be easily combined with additional procedures leading to powerful matheuristics to obtain near optimal solutions for larger instances, see [START_REF] Salhi | Heuristic search: The emerging science of problem solving[END_REF]. This problem can also be tackled by powerful metaheuristics whose performance can be evaluated using lower bounds obtained from the proposed method.

The current problem can be extended to cater for several deterministic and stochastic variants that are worth exploring. These include the consideration of time windows, multiple depots and/or additional location decisions for the selection of depots, as well as periodicity or uncertainty in the customer demand. Moreover, the model and the algorithm we propose in this paper can be easily modified to solve the problem variants where vehicle-dependent energy consumption and routing costs are considered.

N

  = I ∪ J ∪ {0}: the set of all nodes I: the set of customer locations J: the set of potential locations for charging stations 0: the depot node K: the set of vehicles d i > 0: the demand of client i ∈ I c ij : the routing cost of traversing arc (i, j) ∈ A e ij : the energy consumption on arc (i, j) ∈ A expressed in kWh f j : the fixed cost of opening a charging station at node j ∈ J r k : the unit cost of recharging of vehicle k ∈ K v k : the fixed cost of operating vehicle k ∈ K Q k : the load capacity of vehicle k ∈ K β k : the battery capacity of vehicle k ∈ K expressed in kWh.

Figure 1 :

 1 Figure 1: Illustration of minimal energy consumption and recharging need for visiting m copies of j with vehicle k for m = 2, 3.

Figure 2 :

 2 Figure 2: Calculation of minimal energy consumption E C m on the partial network for visiting m = 3 copies of j with vehicle k 1 .

Figure 7 :

 7 Figure 7: Average BDA solving time (Phase I, Phase II, and Phase I+Phase II) and gaps per heterogeneous fleet.

Figure 8

 8 Figure 8 compares the average cost, number of stations opened and the number of vehicles used for each given heterogeneous fleet. As expected, the average cost in general decreases as the fleet size increases. In this Figure, most significant drops in cost occur when the fleet has a type-3 vehicle and another vehicle of type 2 or type 3. These are also the instances where the optimal or best known solutions open fewer number of stations. The highest number of stations are opened with fleets K2V1-2, K3V1-1-2, and K4V1-1-1-2 which have relatively smaller total freight capacity and total range. Thus, solutions to these instances also use more vehicles on average.

Figure 8 :

 8 Figure 8: Average cost, number of stations opened and number of vehicles used for each heterogeneous fleet type.

Figure 9 :

 9 Figure 9: Average cost, number of stations opened and number of vehicles used for each network type (heterogeneous fleets).

Figure 10 :

 10 Figure 10: A comparison of homogeneous and heterogeneous fleet: average cost, number of opened stations and vehicles used for each network size.

Figure 11 :

 11 Figure 11: Average cost for each homogeneous fleet.

  we can replace equality (27) with inequality

	j∈J A	z k j ≥
	w k due to Lemma 4.1.
	Lemma 4.1. If (

Table 1 :

 1 Freight capacity, battery capacity, and cost values for the three vehicle types considered.

	Vehicle type 1	2	3
	Q k	80 100 200
	β k	10 16	22
	v k	16 26	36

Table 2 :

 2 Heterogeneous fleets tested for each network of Schneider et al. (2014). ID Fleet size Vehicle types available Fleet ID Fleet size Vehicle types available K2V1

	Fleet

Table 3 :

 3 Average solution, gap, and solving time values for BDA and PF.

				BDA Phase I	BDA Phase II	PF
	Net. ID |I| |J|	Obj	g	t(s)	Obj	g	t(s) N S N V	Obj	g	t(s)
	c101C5		7 7127.48 0.00	1.32 7127.48 0.00	1.51 0.41 2.12 7127.48 0.00 446.07
	c103C5		7 4259.00 0.00	1.14 4259.00 0.00	1.42 0.00 1.35 4259.00 0.00 170.53
	c206C5		8 8553.95 0.00	2.55 8553.95 0.00	2.56 0.29 2.35 8553.96 0.00 819.52
	c208C5		7 5713.57 0.00	0.83 5713.57 0.00	1.03 0.18 1.35 5713.57 0.00	47.77
	r104C5		7 4214.04 0.00	1.09 4214.04 0.00	1.20 0.00 1.35 4214.04 0.00 182.98
	r105C5 r202C5	5	7 4233.90 0.00 7 3948.66 0.00	0.87 4233.90 0.00 0.92 3948.66 0.00	1.01 0.00 1.35 4233.90 0.00 1.03 0.00 1.59 3948.66 0.00	72.56 27.56
	r203C5		8 6675.67 0.00	2.31 6675.67 0.00	2.28 0.06 2.12 6675.67 0.00 494.06
	rc105C5		8 7421.29 0.00	2.59 7421.29 0.00	2.80 0.06 2.12 7421.28 0.02 1077.36
	rc108C5		8 20738.84 0.00	1.97 20738.84 0.00	2.68 0.65 2.47 20738.84 0.06 815.56
	rc204C5		8 5219.75 0.00	2.61 5219.75 0.00	2.88 0.00 2.00 5219.82 0.05 1548.00
	rc208C5		7 5860.37 0.00	0.62 5860.37 0.00	0.77 0.18 1.35 5860.37 0.06 585.22
	Average				0.00	1.57		0.00	1.76 0.15 1.79	0.02 523.93
	c101c10		14 18721.65 0.00	346.81 17880.34 0.00	259.39 0.81 2.88
	c104c10		13 20193.80 0.00	259.55 20146.89 0.00	235.76 0.82 2.82
	c202c10		14 8208.00 0.00	164.15 8208.00 0.00	215.44 0.56 2.38
	c205c10		12 14160.13 0.00	182.23 14160.13 0.00	374.76 0.59 2.47
	r102c10		13 6936.08 0.00	22.18 6936.08 0.00	20.87 0.06 2.12
	r103c10 r201c10	10	12 4247.15 0.00 13 5895.34 0.00	7.56 4247.15 0.00 237.09 5895.34 0.00	9.23 0.00 1.35 215.06 0.00 2.13
	r203c10		14 9780.68 0.00	151.44 9780.68 0.00	189.97 0.24 2.53
	rc102c10		13 27920.23 0.00	28.69 26254.23 0.00	32.92 1.31 3.00
	rc108c10		13 36706.26 0.00	320.01 36133.71 0.00	735.61 1.18 2.94
	rc201c10		13 14250.35 0.00	368.28 14250.35 0.00	352.18 0.47 2.59
	rc205c10		13 32217.30 0.00	175.72 31550.95 0.00	840.48 1.35 2.82
	Average				0.00 188.53		0.00 286.73 0.62 2.50
	c103c15		19 16455.88 0.07 5622.99 15546.37 0.06 4679.68 0.81 2.88
	c106c15		17 8316.83 0.01 1606.37 8316.83 0.01 1698.20 0.06 2.29
	c202c15		19 33255.84 0.09 4389.71 30922.17 0.06 4295.05 1.31 3.00
	c208c15		18 16909.28 0.08 4809.75 16396.69 0.06 4625.89 0.81 2.94
	r102c15		22 9356.76 0.03 3650.76 9356.76 0.04 4458.27 0.38 2.56
	r105c15 r202c15	15	20 9618.99 0.07 5267.99 9618.99 0.07 5524.17 0.25 2.63 20 17991.78 0.11 4650.14 15944.78 0.09 5091.58 1.20 2.87
	r209c15		19 10064.97 0.02 1237.37 9981.06 0.00 1050.26 0.31 2.56
	rc103c15		19 30891.86 0.17 6556.05 30072.29 0.22 7016.98 1.38 3.00
	rc108c15		19 44373.32 0.16 4815.57 41028.85 0.28 5865.23 1.60 3.00
	rc202c15		19 33912.00 0.16 6017.79 34087.56 0.18 7061.24 1.44 3.25
	rc204c15		21 28632.75 0.22 8130.34 19652.83 0.17 9054.50 1.46 3.31
	Average				0.10 4658.65		0.10 4833.38 0.90 2.85

time limits. The average gap is 0.10.

Acknowledgment: This work was supported by the French National Research Agency [grant number ANR-15-CE22-0017]. We thank the authors of Schneider et al. (2014) and Schiffer and Walther (2017b) for providing the data instances.