Alejandro Montoya
email: jmonto36@eafit.edu.co

Christelle Guéret
email: christelle.gueret@univ-angers.fr

Jorge E Mendoza
email: jorge.mendoza@univ-tours.fr

Juan G Villegas
email: juan.villegas@udea.edu.co

The electric vehicle routing problem with nonlinear charging function

Keywords: Vehicle routing problem, Electric vehicle routing problem with nonlinear charging function, Iterated local search (ILS), Matheuristic

Electric vehicle routing problems (E-VRPs) extend classical routing problems to consider the limited driving range of electric vehicles. In general, this limitation is overcome by introducing planned detours to battery charging stations. Most existing E-VRP models assume that the battery-charge level is a linear function of the charging time, but in reality the function is nonlinear. In this paper we extend current E-VRP models to consider nonlinear charging functions. We propose a hybrid metaheuristic that combines simple components from the literature and components specifically designed for this problem. To assess the importance of nonlinear charging functions, we present a computational study comparing our assumptions with those commonly made in the literature. Our results suggest that neglecting nonlinear charging may lead to infeasible or overly expensive solutions. Furthermore, to test our hybrid metaheuristic we propose a new 120-instance testbed. The results show that our method performs well on these instances.

Introduction

In the last few years several companies have started to use electric vehicles (EVs) in their operations. For example, La Poste operates at least 250 EVs and has signed orders for an additional 10,000 [START_REF] Kleindorfer | Fleet renewal with electric vehicles at La Poste[END_REF]; and the French electricity distribution company ENEDIS runs 2,000 EVs, accounting for 10% of their fleet in 20161 . Despite these developments, the large-scale adoption of EVs for service and distribution operations is still hampered by technical constraints such as battery charging times and limited battery capacity. For the most common EVs used in service operations, the minimum charging time is 0.5 h and the battery capacity is around 22 kWh. The latter leads to a nominal driving range of 142 km [START_REF] Pelletier | Battery electric vehicles for goods distribution: a survey of vehicle technology, market penetration, incentives and practices[END_REF]. In reality, the driving range could be significantly lower because the energy consumption increases with the slope of the road, the speed, and the use of peripherals [START_REF] De Cauwer | Energy consumption prediction for electric vehicles based on real-world data[END_REF]. For instance, [START_REF] Restrepo | Performance testing of electric vehicles on operating conditions in Bogotá DC, Colombia[END_REF] documented that the heating and air conditioning respectively reduce the driving range of an EV by about 30% and 8% per hour of use.

Automakers and battery manufacturers are investing significant amounts of capital and effort into the development of new technology to improve EV autonomy and charging time. For instance, General Motors (GM) reinvested USD 20 million into the GM Global Battery Systems Lab to help the company developing new battery technology for their vehicles [START_REF] Marcacci | GM triples size of its michigan electric vehicle battery laboratory (CT exclusive)[END_REF]. The results of these efforts, however, are transferred only slowly to commercially available EVs. In the meantime, companies using EVs in their daily operations need fleet management tools that can take into account limited driving ranges and slow charging times [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF]. To respond to this challenge, around 2012 the operations research community started to study a new family of vehicle routing problems (VRPs): the so-called electric VRPs (E-VRPs) [START_REF] Afroditi | Electric vehicle routing problem with industry constraints: Trends and insights for future research[END_REF][START_REF] Pelletier | 50th anniversary invited article-goods distribution with electric vehicles: Review and research perspectives[END_REF]). These problems consider the technical limitations of EVs. Because of the short driving range, E-VRP solutions frequently include routes with planned detours to charging stations (CSs). The need to detour usually arises in rural and semi-urban operations, where the distance covered by the routes on a single day is often higher than the driving range.

As has been the case for other optimization problems inspired by practical applications, research in E-VRPs started with primarily theoretical variants and is slowly moving toward problems that better capture reality. In general, E-VRP models make assumptions about the EV energy consumption, the charging infrastructure ownership, the capacity of the CSs, and the battery charging process. Most E-VRPs assume that energy consumption is directly and exclusively related to the traveled distance. However, as mentioned before, the consumption depends on a number of additional factors. To the best of our knowledge only [START_REF] Goeke | Routing a mixed fleet of electric and conventional vehicles[END_REF] and [START_REF] Lin | Electric vehicle routing problem[END_REF] use consumptions computed over actual road networks taking into account the EV parameters and their loads.

Similarly, most E-VRP models implicitly assume that the charging infrastructure is private. In this context, the decision-maker controls access to the CSs, so they are always available. However, in reality, mid-route charging is often performed at public stations and so the availability is uncertain. To our knowledge only [START_REF] Sweda | Adaptive routing and recharging policies for electric vehicles[END_REF] and [START_REF] Kullman | Electric vehicle routing with mid-route recharging and uncertain charging station availability[END_REF] deal with public infrastructure and consider uncertainty in CS availability.

CS capacity is another area in which current E-VRP models are still a step behind reality. All existing E-VRP research that we are aware of assumes that the CSs can simultaneously handle an unlimited number of EVs. In practice, each CS is usually equipped with only a few chargers. In some settings this assumption may be mild (e.g., a few geographically distant routes and private CSs). However, in most practical applications CS capacity plays a restrictive role.

Finally, in terms of the battery charging process, E-VRP models make assumptions about the charging policy and the charging function approximation. The former defines how much of the battery capacity can (or must) be restored when an EV visits a CS, and the latter models the relationship between battery charging time and battery level. In this paper, we focus on these assumptions.

In terms of the charging policies, the E-VRP literature can be classified into two groups: studies assuming full and partial charging policies. As the name suggests, in full charging policies, the battery capacity is fully restored every time an EV reaches a CS. Some studies in this group assume that the charging time is constant [START_REF] Conrad | The recharging vehicle routing problem[END_REF][START_REF] Erdogan | A green vehicle routing problem[END_REF][START_REF] Montoya | A multi-space sampling heuristic for the green vehicle routing problem[END_REF]. This is a plausible assumption in applications where the CSs replace a (partially) depleted battery with a fully charged one. Other researchers, including [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF], [START_REF] Goeke | Routing a mixed fleet of electric and conventional vehicles[END_REF], [START_REF] Schneider | An adaptive VNS algorithm for vehicle routing problems with intermediate stops[END_REF], [START_REF] Desaulniers | Exact algorithms for electric vehicle-routing problems with time windows[END_REF], [START_REF] Hiermann | The electric fleet size and mix vehicle routing problem with time windows and recharging stations[END_REF], [START_REF] Lin | Electric vehicle routing problem[END_REF], [START_REF] Szeto | Artificial bee colony approach to solving the electric vehicle routing problem[END_REF], consider full charging policies with a linear charging function approximation (i.e., the battery level is assumed to be a linear function of the charging time). In their models, the time spent at each CS depends on the battery level when the EV arrives and on the (constant) charging rate of the CS. In partial charging policies, the level of charge (and thus the time spent at each CS) is a decision variable. To the best of our knowledge, all existing E-VRP models with partial charging consider linear function approximations [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF][START_REF] Sassi | Iterated tabu search for the mix fleet vehicle routing problem with heterogenous electric vehicles[END_REF][START_REF] Bruglieri | A variable neighborhood search branching for the electric vehicle routing problem with time windows[END_REF][START_REF] Schiffer | The electric location routing problem with time windows and partial recharging[END_REF][START_REF] Desaulniers | Exact algorithms for electric vehicle-routing problems with time windows[END_REF][START_REF] Keskin | Partial recharge strategies for the electric vehicle routing problem with time windows[END_REF].

In general, the charging functions are nonlinear, because the terminal voltage and current change during the charging process. This process is divided into two phases. In the first phase, the charging current is held constant, and thus the battery level increases linearly with time. The first charging phase continues until the battery's terminal voltage increases to a specific maximum value (see Figure 1). In the second phase, the current decreases exponentially and the terminal voltage is held constant to avoid battery damage. The battery level then increases concavely with time [START_REF] Pelletier | Goods distribution with electric vehicles: Battery degradation and behaviour modeling[END_REF]. Although the shape of the charging functions is known, devising analytical expressions to model them is complex because they depend on factors such as current, voltage, self-recovery, and temperature [START_REF] Wang | Estimation of state of charge of batteries for electric vehicles[END_REF]). The battery level is then described by differential equations. Since such equations are difficult to incorporate into E-VRP models, researchers rely on approximations of the actual charging functions. Bruglieri et al. (2014) use a linear approximation that considers only the linear segment of the charging function, i.e., between 0 and (around) 0.8Q, where Q represents the battery capacity. This approximation avoids dealing with the nonlinear segment of the charging function (i.e., from (around) 0.8Q to Q). He henceforth refer to this approximation as FS. [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF], [START_REF] Sassi | Vehicle routing problem with mixed fleet of conventional and heterogenous electric vehicles and time dependent charging costs[END_REF], [START_REF] Bruglieri | A variable neighborhood search branching for the electric vehicle routing problem with time windows[END_REF], [START_REF] Desaulniers | Exact algorithms for electric vehicle-routing problems with time windows[END_REF], [START_REF] Schiffer | The electric location routing problem with time windows and partial recharging[END_REF], and [START_REF] Keskin | Partial recharge strategies for the electric vehicle routing problem with time windows[END_REF] approximate the whole charging function using a linear expression. They do not explain how the approximation is calculated, but two options can be considered. In the first (L1) the charging rate of the function corresponds to the slope of its linear segment (see Figure 2b). This approximation is optimistic: it assumes that batteries charge to the level Q faster than they do in reality. In the second approximation (L2) the charging rate is the slope of the line connecting the first and last observations (see Figure 2c) of the charging curve. This approximation tends to be pessimistic: over a large portion of the curve, the charging rate is slower than in reality. In this paper, we study a new E-VRP that captures the nonlinear behavior of the charging process using a piecewise linear approximation. The main contributions of this research are fivefold. First, we introduce the electric vehicle routing problem with nonlinear charging functions (E-VRP-NL). Second, we propose a hybrid metaheuristic, which combines simple components from the literature and components specifically designed for this new problem. Third, we propose a set of realistic and publicly available instances. Fourth, we demonstrate through extensive computational experiments the importance of better approximating the actual battery charging function. Fifth, we analyze our solutions and provide some insight into the characteristics of good E-VRP-NL solutions.

The remainder of this paper is organized as follows. Section 2 formally introduces the E-VRP-NL. Section 3 describes our hybrid metaheuristic, and Section 4 presents the computational experiments. Finally, Section 5 concludes the paper and discusses future research.

Electric vehicle routing problem with nonlinear charging function

Problem description

Let I be the set of nodes representing the customers, F the set of CSs, and 0 a node representing the depot. Each customer i ∈ I has a service time p i . The E-VRP-NL is defined on a directed and complete graph G = (V, A), where V = {0} ∪ I ∪ F and F contains the set F and β copies of each CS (i.e., |F | = |F | × (1 + β)). The value of 1 + β corresponds to the number of times that each CS can be visited. Let A = {(i, j) : i, j ∈ V, i = j} be the set of arcs connecting vertices of V . Each arc (i, j) has two associated nonnegative values: a travel time t ij and an energy consumption e ij . The customers are served using an unlimited and homogeneous fleet of EVs. All the EVs have a battery of capacity Q (expressed in kWh) and a maximum tour duration T max . It is assumed that the EVs leave the depot with a fully charged battery, and that all the CSs can handle an unlimited number of EVs simultaneously. Feasible solutions to the E-VRP-NL satisfy the following conditions: each customer is visited exactly once; each route satisfies the maximum-duration limit; each route starts and ends at the depot; and the battery level when an EV arrives at and departs from any vertex is between 0 and Q.

Since the traveled distance is directly related to the energy consumption, most work on E-VRPs with a homogeneous fleet focuses on minimizing the total distance [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF][START_REF] Desaulniers | Exact algorithms for electric vehicle-routing problems with time windows[END_REF][START_REF] Hiermann | The electric fleet size and mix vehicle routing problem with time windows and recharging stations[END_REF][START_REF] Keskin | Partial recharge strategies for the electric vehicle routing problem with time windows[END_REF]. However, this objective function neglects the impact of charging operations. This may lead to solutions that charge the batteries more than needed or charge them even when their level is high. These decisions directly affect the battery's long-term degradation cost (which according to [START_REF] Becker | The vehicle relocation problem for the one-way electric vehicle sharing: An application to the Milan case[END_REF] can be three times the energy cost) and the charging fees at CSs [START_REF] Bansal | Charging of electric vehicles: Technology and policy implications[END_REF]. To better capture the impact of charging operations, in the E-VRP-NL we minimize the total travel and charging time. This objective function was studied by [START_REF] Zündorf | Electric vehicle routing with realistic recharging models[END_REF] and [START_REF] Liao | The electric vehicle touring problem[END_REF] for related routing problems.

Modeling of battery charging functions

Each CS i ∈ F has a charging mode (e.g., slow, moderate, fast) that is associated with a charging function g i (q i , ∆ i ). This function maps the charge level when the vehicle arrives at i (q i ) and the time spent charging at i (∆ i ) to the charge level when the vehicle leaves i. To avoid handling a two-dimensional function, we use the transformation proposed by [START_REF] Zündorf | Electric vehicle routing with realistic recharging models[END_REF]. Let ĝi (l) be the charging function when q i = 0 and the battery is charged for l time units; g i (q i , ∆ i ) is estimated as ĝi (∆ i + ĝ-1 (q i )). Note that ∆ i = ĝi -1 (o i ) -ĝi -1 (q i ), where o i is the charge level when the vehicle leaves i.

The function ĝi (l) is concave (Bruglieri et al. 2014[START_REF] Hõimoja | An ultrafast ev charging station demonstrator[END_REF][START_REF] Pelletier | Goods distribution with electric vehicles: Battery degradation and behaviour modeling[END_REF] with an asymptote at Q. Similarly to [START_REF] Zündorf | Electric vehicle routing with realistic recharging models[END_REF], we argue that ĝi (l) can be accurately approximated using piecewise linear functions. We support our claim using the data provided by [START_REF] Uhrig | E-mobility in car parks-guidelines for charging infrastructure expansion planning and operation based on stochastic simulations[END_REF]. We fit piecewise linear functions to their data and obtain approximations with an average relative absolute error of 0.90%, 1.24%, and 1.90% for CSs of 11, 22, and 44 kW, respectively. Figure 3 shows the piecewise linear approximation for a CS i of 22 kW charging a vehicle equipped with a battery of 16 kWh. In the plot, c ik and a ik represent the charging time and the charge level for the breakpoint k ∈ B of the CS i ∈ F , where B = {0, .., b} is the set of breakpoints of the piecewise linear approximation.

𝑐 𝑖1 𝑐 𝑖2 𝑐 𝑖3 𝑐 𝑖0 𝑎 𝑖0 𝑎 𝑖1 𝑎 𝑖2 𝑎 𝑖3 = 𝑄

Battery level

Time

Figure 3: Real data vs. piecewise linear approximation for a CS with 22 kW charging a battery of 16 kWh.

Illustrative example

Figure 4 presents a numerical example illustrating the E-VRP-NL. The figure depicts a solution to an instance with 7 customers and 3 CSs. The CSs have different technologies (slow and fast), and each technology has a specific piecewise-linear charging function. The charging function maps the battery levels q i and o i to the charging times s i and d i to estimate the time spent at the CS i ∈ F (∆ i ). In this example, Route 1 does not visit any CS, because its total energy consumption is less than the battery capacity. Route 2 visits CS 8: the EV arrives at the CS with a battery level q 8 = 1.0, and it charges the battery to a level o 8 = 6.0. To estimate the time spent at the CS, we use the piecewise-linear charging function: the charging times associated with q 8 and o 8 are s i = 0.8 and d i = 6.0, so the time spent at CS 8 is ∆ 8 = 6.0 -0.8 = 5.2. The duration of Route 2 is the sum of the travel time (13.0), the charging time (5.2), and the service time (1.0), i.e., 19.2, which is less than T max. The cost of this route is 18.2 (travel time + charging time). Finally, Route 3 visits CSs 10 and 9, and it spends ∆ 10 = 7.2 and ∆ 9 = 1.6 time units charging in these CSs, respectively.

Mixed-integer linear programming formulation

To help the reader understand the E-VRP-NL, we now provide a mixed integer linear programming (MILP) formulation of the problem. The MILP uses the following decision variables: variable x ij is equal to 1 if an EV travels from vertex i to j, and 0 otherwise. Variables τ j and y j track the time and charge level when the EV departs from vertex j ∈ V . Variables q i and o i specify the charge levels when an EV arrives at and departs from CS i ∈ F , and s i and d i are the associated charging times. Variable ∆ i = d i -s i represents the time spent at CS i ∈ F . Variables z ik and w ik are equal to 1 if the charge level is between a i,k-1 and a ik , with k ∈ B \ {0}, when the EV arrives at and departs from CS i ∈ F respectively. Finally, variables α ik and λ ik are the coefficients of the breakpoint k ∈ B in the piecewise linear approximation, when the EV arrives at and departs from CS i ∈ F respectively. The MILP formulation follows:

min i,j∈V t ij x ij + i∈F ∆ i (1) subject to j∈V,i =j x ij = 1, ∀i ∈ I (2) j∈V,i =j x ij ≤ 1, ∀i ∈ F (3) j∈V,i =j x ji - j∈V,i =j x ij = 0, ∀i ∈ V (4) e ij x ij -(1 -x ij )Q ≤ y i -y j ≤ e ij x ij + (1 -x ij )Q, ∀i ∈ V, ∀j ∈ I (5) e ij x ij -(1 -x ij )Q ≤ y i -q j ≤ e ij x ij + (1 -x ij )Q, ∀i ∈ V, ∀j ∈ F (6) y i ≥ e i0 x i0 , ∀i ∈ V (7)
y i = o i , ∀i ∈ F (8) y 0 = Q (9) q i ≤ o i , ∀i ∈ F ( 10
)
q i = k∈B α ik a ik , ∀i ∈ F ( 11
)
s i = k∈B α ik c ik , ∀i ∈ F (12) k∈B α ik = k∈B z ik , ∀i ∈ F (13) k∈B z ik = j∈V x ij , ∀i ∈ F (14) α ik ≤ z ik + z i,k+1 , ∀i ∈ F , ∀k ∈ B \ {b} (15) α ib ≤ z ib , ∀i ∈ F (16)
o i = k∈B λ ik a ik , ∀i ∈ F (17)
d i = k∈B λ ik c ik , ∀i ∈ F (18) k∈B λ ik = k∈B w ik , ∀i ∈ F (19) k∈B w ik = j∈V x ij , ∀i ∈ F (20) λ ik ≤ w ik + w i,k+1 , ∀i ∈ F, ∀k ∈ B \ {b} (21) λ ib ≤ w ib , ∀i ∈ F (22) ∆ i = d i -s i , ∀i ∈ F (23) τ i + (t ij + p j )x ij -T max (1 -x ij ) ≤ τ j , ∀i ∈ V, ∀j ∈ I (24)
τ i + ∆ j + t ij x ij -(S max + T max )(1 -x ij ) ≤ τ j , ∀i ∈ V, ∀j ∈ F (25) τ j + t j0 ≤ T max , ∀j ∈ V (26) τ 0 ≤ T max (27) x ij = 0, ∀i, j ∈ F : m ij = 1 (28) τ i ≥ τ j , ∀i, j ∈ F : m ij = 1, j ≤ i (29) τ j ≤ T max i∈V x ij , ∀j ∈ F (30) i∈V x ih ≥ j∈V x jf , ∀h, f ∈ F : m hf = 1, h ≤ f (31) x ij ∈ {0, 1}, ∀i, j ∈ V (32) τ i ≥ 0, y i ≥ 0 ∀i ∈ V (33) z ik ∈ {0, 1}, w ik ∈ {0, 1}, α ik ≥ 0, λ ik ≥ 0, ∀i ∈ F , ∀k ∈ B ( 34
)
q i ≥ 0, o i ≥ 0, s i ≥ 0, d i ≥ 0, ∆ i ≥ 0, ∀i ∈ F (35)
The objective function (1) seeks to minimize the total time (travel times plus charging times). Constraints (2) ensure that each customer is visited once. Constraints (3) ensure that each CS copy is visited at most once. Constraints (4) impose the flow conservation. Constraints ( 5) and ( 6) track the battery charge level at each vertex. Constraints (7) ensure that if the EV travels between a vertex and the depot, it has sufficient energy to reach its destination. Constraints (8) reset the battery tracking to o i upon departure from CS i ∈ F . Constraint (9) ensures that the battery charge level is Q at the depot. Constraints (10) couple the charge levels when an EV arrives at and departs from any CS. Constraints ( 11)-( 16) define the charge level (and its corresponding charging time) when an EV arrives at CS i ∈ F (based on the piecewise linear approximation of the charging function). Similarly, constraints ( 17)-( 22) define the charge level (and its corresponding charging time) when an EV departs from CS i ∈ F . Constraints (23) define the time spent at any CS. Constraints ( 24) and ( 25) track the departure time at each vertex, where S max = max i∈F {c ib }. Constraints ( 26) and ( 27) ensure that the EVs return to the depot no later than T max . Constraints ( 28) and ( 31) help to avoid the symmetry generated by the copies of the CSs. Parameter m ij is equal to 1 if i and j ∈ F represent the same CS. Finally, constraints (32)-( 35) define the domain of the decision variables.

Solving the E-VRP-NL

Lenstra & Kan (1981) demonstrated that the classical VRP, commonly known as the CVRP, is NP-hard. Since the CVRP is a special case of our E-VRP-NL, the latter is also NPhard. We therefore propose a metaheuristic approach. Like many metaheuristics for other VRPs, our approach explores new solutions by building new routes or applying changes (moves) to existing ones. In this process, the algorithm makes sequencing and charging decisions. The former fix the order in which the route visits its assigned customers, while the latter determine where and how much to charge the EV serving the route. Sequencing and charging decisions can be made either simultaneously or in two phases (sequencing first, charging second). We use the latter option. To make charging decisions, we solve what we call the fixed route vehicle charging problem, or simply FRVCP. In a nutshell, the problem consists in i) inserting charging stations into a fixed sequence of customers and ii) deciding how much to charge at each inserted station. Since the FRVCP plays a key role in our approach, we discuss it in the next subsection before introducing our metaheuristic.

The fixed-route vehicle-charging problem

The FRVCP is a variant of the well-known fixed-route vehicle-refueling problem (FRVRP). The FRVRP seeks the minimum-cost refueling policy (which fuel stations to visit and the refueling quantity at each station) for a given origin-destination route [START_REF] Suzuki | A variable-reduction technique for the fixed-route vehicle-refueling problem[END_REF]. Most of the research into the FRVRP and its variants applies only to internal combustion vehicles (which have negligible refueling times), but a few extensions to EVs have been reported.

Most of these extensions assume full charging policies [START_REF] Montoya | A multi-space sampling heuristic for the green vehicle routing problem[END_REF][START_REF] Hiermann | The electric fleet size and mix vehicle routing problem with time windows and recharging stations[END_REF][START_REF] Liao | The electric vehicle touring problem[END_REF]). To our knowledge, only [START_REF] Sweda | Optimal recharging policies for electric vehicles[END_REF] assume a partial charging policy. Their problem differs from ours in three fundamental ways: i) they do not take into account the charging times (because their objective is to minimize the energy and degradation costs), ii) they do not deal with maximum route duration constraints, and iii) their CSs are already included in the fixed route and no detours are to be planned.

In the FRVCP the objective is to find the charging decisions (where and how much to charge) that minimize the sum of the charging times and detour times while satisfying the following conditions: the level of the battery when the EV arrives at any vertex is nonnegative; the charge in the battery does not exceed its capacity; and the route satisfies the maximum-duration limit. Since the FRVRP is NP-hard [START_REF] Suzuki | A variable-reduction technique for the fixed-route vehicle-refueling problem[END_REF]) and the FRVCP generalizes the FRVRP, we can conclude that the FRVCP is also NP-hard.

Let Π = {π(0), π(1), ..., π(i), ..., π(j), ...π(n r -1), π(n r )} be the fixed route, where π(0) and π(n r ) represent the depot and π(1),

• • • , π(i), • • • , π(j), • • • , π(n r -1) the customers.
The fixed route has a total time t, which is the sum of the travel times plus the service times. Note that by definition i) Π does not visit CSs, and ii) it is energy infeasible (i.e., requires more than one full battery to complete)3 . The feasibility of Π may be restored by inserting visits to CSs. As mentioned in Section 3, each CS j ∈ F has a piecewise-linear charging function defined by a set of breakpoints B. Each segment of the piecewise linear function is defined between breakpoints k -1 and k ∈ B, has a slope ρ jk (representing a charging rate), and is bounded by the battery levels a jk-1 and a jk . The values e π(i-1)π(i) and t π(i-1)π(i) represent the energy consumption and the travel time between vertices π(i-1) and π(i) ∈ Π. Similarly, e π(i-1)j and t π(i-1)j represent the energy consumption and the travel time between vertex π(i -1) ∈ Π and CS j ∈ F , and e jπ(i) and t jπ(i) represent the energy consumption and the travel time between CS j ∈ F and vertex π(i) ∈ Π. Figure 5 presents an example illustrating the FRVCP using the fixed route corresponding to Route 3 in the example introduced in Section 2. Figure 5a shows the fixed route with 3 customers. Figure 5b shows all the possible CS insertions into the fixed route (the arcs in bold correspond to the FRVCP solution). Figure 5c shows the energy-feasible route resulting from solving the FRVCP. We formulate the FRVCP as an MILP using the following decision variables: variable

ε π(i)j is equal to 1 if the EV charges at CS j ∈ F before visiting vertex π(i) ∈ Π. Variable φ π(i) tracks the battery level. If ε π(i)j = 0, φ π(i) is the battery level when the EV arrives at vertex π(i). If ε π(i)j = 1, φ π(i)
is the battery level when the EV arrives at CS j ∈ F immediately before visiting vertex π(i). Variable θ π(i)jk is equal to 1 if the EV charges on the segment defined by breakpoints k -1 and k ∈ B at CS j ∈ F before visiting vertex π(i) ∈ Π. Finally, variables δ π(i)jk and µ π(i)jk are (respectively) the amount of energy charged and the battery level when the charging process finishes on the segment between breakpoints k -1 and k ∈ B at CS j ∈ F before the visit to vertex π(i) ∈ Π. The MILP formulation of the FRVCP follows:

min π(i)∈Π\{π(0)} j∈F k∈B\{0} δ π(i)jk ρ jk + π(i)∈Π\{π(0)} j∈F ε π(i)j (t π(i-1)j + t jπ(i) -t π(i-1)π(i) ) (36) subject to φ π(1) = Q - j∈F ε π(1)j e π(0)j -e π(0)π(1)   1 - j∈F ε π(1)j   (37) φ π(i) = φ π(i-1) + j∈F k∈B\{0} δ π(i-1)jk - j∈F ε π(i-1)j e jπ(i-1) - j∈F ε π(i)j e π(i-1)j -e π(i-1)π(i)   1 - j∈F ε π(i)j   ∀π(i) ∈ Π \ {π(0), π(1), π(n r )} (38) φ π(nr) = φ π(nr-1) + j∈F k∈B\{0} δ π(nr-1)jk + j∈F k∈B\{0} δ π(nr)jk - j∈F ε π(nr-1)j e jπ(nr-1) - j∈F ε π(nr)j (e π(nr-1)j + e jπ(nr) ) -e π(nr-1)π(nr)   1 - j∈F ε π(nr)j   (39) φ π(nr-1) + j∈F k∈B\{0}
δ π(nr-1)jk -j∈F e jπ(nr-1) ε π(nr-1)j -j∈F e π(nr-1)j ε π(nr)j ≥ 0 (40)

µ π(i)j1 = φ π(i) + δ π(i)j1 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (41) µ π(nr)j1 = φ π(nr-1) + l∈F k∈B\{0} δ π(nr-1)lk - l∈F e lπ(nr-1) ε π(nr-1)l -ε π(nr)j e π(nr-1)j + δ π(nr)j1 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (42) µ π(i)jk = µ π(i)j,k-1 + δ π(i)jk ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0, 1} (43) µ π(i)jk ≥ a jk-1 θ π(i)jk ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0, 1} (44)
µ π(i)jk ≤ a jk θ π(i)jk + (1 -θ π(i)jk )Q ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (45) j∈F ε π(i)j ≤ 1, ∀π(i) ∈ Π \ {π(0)} (46) θ π(i)jk ≤ ε π(i)j ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (47) δ π(i)jk ≤ θ π(i)jk Q ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (48) t + π(i)∈Π\{π(0)} j∈F k∈B\{0} δ π(i)jk ρ jk + π(i)∈Π j∈F ε π(i)j (t π(i-1)j + t jπ(i) ) -t π(i-1)π(i) ) ≤ T max (49) φ π(i) ≥ 0, ∀π i ∈ Π \ {π(0)} (50) ε π(i)j ∈ {0, 1}, ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F (51) θ π(i)jk ∈ {0, 1} ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (52) δ π(i)jk ≥ 0 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (53) µ π(i)jk ≥ 0 ∀π(i) ∈ Π \ {π(0)}, ∀j ∈ F, ∀k ∈ B \ {0} (54)
The objective function (36) seeks to minimize the total route time (including charging and detour times). Constraints (37)-( 40) define the battery level when the EV arrives at vertex 41)-( 43) define the battery level when the EV finishes charging at CS j ∈ F in the segment between breakpoints k -1 and k ∈ B before visiting vertex π(i) ∈ Π. Constraints ( 44)-( 45) ensure that if the EV charges on a given segment, the battery level is between the values of its corresponding break points (a j,k-1 and a jk ). Constraints ( 46) state that only one CS is visited between any two vertices of the fixed route. Constraints (47) ensure that the EV uses only segments of the visited CSs. Likewise, constraints (48) ensure that the EV charges only at the selected segments of the visited CSs. Constraint (49) represents the duration constraint of the route. Finally, constraints ( 50)-( 54) define the domain of the decision variables.

π(i) ∈ Π if ε π(i)j = 0; or at CS j ∈ F before visiting vertex π(i) ∈ Π, if ε π(i)j = 1. Constraints (
Our metaheuristic for the E-VRP-NL solves the FRVCP at various steps. It uses two different approaches: a commercial solver running on the MILP formulation introduced above (boosted by tailored preprocessing strategies) and a greedy heuristic. For the sake of brevity, these approaches are not discussed in the main body of the paper; full details can be found in AppendixA.

Hybrid metaheuristic

To solve the E-VRP-NL we developed a hybrid metaheuristic combining an iterated local search (ILS) and a heuristic concentration (HC). The former is a metaheuristic that starts by generating an initial solution (with a constructive heuristic). This solution is then improved by a local search procedure. At each iteration of the ILS, the best current solution is perturbed, and a new ILS iteration starts from the perturbed solution. More details of the ILS can be found in [START_REF] Lourenço | Iterated local search: Framework and applications[END_REF]. The HC is an approach that tries to build a global optimum using parts of the local optima found during a heuristic search procedure [START_REF] Rosing | Heuristic concentration: Two stage solution construction[END_REF].

Our ILS+HC starts from an initial solution generated using a sequence-first split-second approach. The latter uses a nearest-neighbor heuristic [START_REF] Rosenkrantz | An analysis of several heuristics for the traveling salesman problem[END_REF] to build a TSP tour visiting all the customers and a splitting procedure to find an E-VRP-NL solution. Then, at each iteration of the ILS we improve the current solution using a variable neighborhood descent (VND; [START_REF] Mladenović | Variable neighborhood search[END_REF] with three local search operators: relocate, 2-Opt, and global charging improvement (GCI). At the end of each ILS iteration, we update the best solution and add the routes of the local optimum to a pool of routes Ω ⊂ R, where R is the set of all feasible routes. To diversify the search, we concatenate the routes of the local optimum to build a new TSP tour, and then perturb the new TSP tour. We start a new ILS iteration by splitting the perturbed TSP tour. After K iterations the ILS ends, and we carry out the HC. In this phase, we solve a set partitioning problem over the set of routes Ω to obtain an E-VRP-NL solution. In the remainder of this section, we describe the main components of our method.

Split

To extract a feasible solution from a TSP tour, we use an adaptation of the splitting procedure introduced by [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]. The splitting procedure builds a directed acyclic graph G * = (V * , A * ) composed of the ordered vertex set V * = (v 0 , v 1 , ..., v i , ..., v n ) and the arc set A * . Vertex v 0 = 0 is an auxiliary vertex, and each vertex v i represents the customer in the ith position of the TSP tour. Arc (v i , v i+nr ) ∈ A * represents a feasible route r v i ,v i+nr with a travel time t rv i ,v i+nr , starting and ending at the depot and visiting customers in the sequence v i+1 to v i+nr .

Note that since the TSP tour includes only customers, route r v i ,v i+nr may be energyinfeasible; in that case, we try to repair it by solving an FRVCP. If inserting CSs into r v i ,v i+nr increases the duration of the route beyond T max , we do not include the arc associated with the route in A * . Finally, to obtain a feasible E-VRP-NL solution, the splitting procedure finds the set of arcs (i.e., routes) along the shortest path connecting 0 and v n in G * .

Figure 6 illustrates the tour splitting procedure using the example from Section 2. Figure 6a shows the TSP tour. Figure 6b shows the auxiliary graph G * , where the arcs in bold correspond to the shortest path. Figure 6c shows the solution found by the splitting procedure.

Variable neighborhood descent

To improve the solution generated by the splitting procedure we use a VND based on three local search operators. The first two operators, namely, relocate and 2-Opt4 , focus on the sequencing decisions. In other words, these two operators alter only the sequence of customers and do not insert, remove, or change the position of CSs. To update the charging times after a relocate or 2-Opt move we use the rule proposed by [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF]: when visiting a CS, charge the strict minimum amount of energy needed to continue to the next CS (or the depot if there is no CS downstream). If reaching the next CS (or the depot) is impossible, even with a fully charged battery, the move is deemed infeasible. Similarly, if after updating the charging times the resulting route is infeasible in terms of the maximumduration limit, the move is discarded. It is worth noting that the [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF] rule is optimal when all the CSs are homogeneous; however, this is not the case in our E-VRP-NL.

As its name suggests, the third operator, GCI, focuses on the charging decisions. GCI is applied to every route visiting at least one CS. First, it removes from the route all CS visits. If the resulting route is energy-feasible, it stops. If the route is energy-infeasible, it solves an FRVCP to optimize the charging decisions for that route.

Perturb

To diversify the search, we concatenate the routes of the current best solution to build a TSP tour. Then, we perturb the resulting tour with a randomized double bridge operator [START_REF] Lourenço | Iterated local search: Framework and applications[END_REF]) and apply the split procedure to obtain a new E-VRP-NL solution. The randomized double bridge operator cuts four arcs and introduces four new ones. Figure 7 illustrates the steps of the perturbation procedure.

Heuristic concentration

Finally, the HC component solves a set partitioning formulation over the pool of routes Ω: min R⊆Ω r∈R t r : ∪ r∈R = V ; r i ∩ r j = 0 ∀r i , r j ∈ R . The objective is to select the best subset of routes from Ω to build the set of routes R (i.e., the final solution) guaranteeing that each customer will be visited by exactly one route.

Computational experiments

In this section, we present three computational studies. The first study assesses the benefits of better approximating the battery charging function. The second study evaluates the performance of our ILS+HC. The third study analyzes the characteristics of the best known solutions (BKSs) found by our ILS+HC. The goal of this analysis is to provide researchers with insight that may be useful in the design of new solution methods for the E-VRP-NL.

Test instances for the E-VRP-NL

We generated a new 120-instance testbed built using real data for EV configuration and battery charging functions. To ensure feasibility, we opted to generate our instances instead of adapting existing datasets from the literature. To build the instances we first generated 30 sets of customer locations with {10, 20, 40, 80, 160, 320} customers. For each instance size, we generated 5 sets of customer locations. We located the customers in a geographic space of 120 x 120 km using either a random uniform distribution, a random clustered distribution, or a mixture of both. For each of the 30 sets of locations we chose the customer location strategy using a uniform probability distribution. Our main motivation for choosing the 120 x 120 km area was to build instances representing a semi-urban operation.

For each of the 30 sets of locations we built 4 instances varying the level of charging infrastructure availability and the strategy used to locate the CSs. We considered two levels of charging infrastructure availability: low and high. To favor feasibility, for each combination of number of customers and infrastructure availability level we handpicked the number of CSs as a proportion of the number of customers. We located the CSs either randomly or using a simple p-median heuristic. Our p-median heuristic starts from a set of randomly generated CS locations and iteratively moves those locations trying to minimize the total distance between the CSs and the customers. We included three types of CSs: slow, moderate, and fast. For each CS we randomly selected the type using a uniform probability distribution.

The EVs in our instances are Peugeot iOns. This EV has a consumption rate of 0.125 kWh/km, and a battery of 16 kWh. As mentioned in Section 1, the energy consumption on an arc (e ij ) depends on various factors. For simplicity we followed the classical approach in the literature and assumed that this consumption is simply the EV's consumption rate multiplied by the arc's length. To generate the charging functions we fit piecewise linear functions to the real charging data provided by [START_REF] Uhrig | E-mobility in car parks-guidelines for charging infrastructure expansion planning and operation based on stochastic simulations[END_REF]. Figure 8 depicts our piecewise linear approximations. Finally, we set the maximum route duration for every instance to 10 h. Our 120 instances are publicly available at www.vrp-rep.org5 .

Benefits of better approximating the charging function

To assess the value of a charging function approximation that captures the nonlinear behavior of the process, we conducted an experiment comparing our approximation with those commonly used in the literature. The experiment consists in solving a subset of instances with four charging function approximations (i.e., FS, L1, L2, and our piecewise linear -hereafter called PL) and comparing the solutions in terms of objective function and feasibility. Since PL generalizes FS, L1, and L2, any method for the E-VRP-NL can be adapted to work with the other three approximations by a manipulation of the input data. To avoid the bias introduced by the solution method, we compare only optimal solutions delivered by the MILP introduced in Section 2 (running on Gurobi 5.6). This choice restricted the size of the instances used in the experiment to 10 customers and 3 CSs. We believe, however, that our conclusions hold for any instance size. All the experiments were run on a computing cluster with 2.33 GHz Inter Xeon E5410 processors with 16 GB of RAM running under Linux Rocks 6.1.1.

As mentioned in Section 3, the MILP formulation uses β copies of the CSs to model multiple visits to the same CS. Although several authors have used this strategy [START_REF] Conrad | The recharging vehicle routing problem[END_REF][START_REF] Erdogan | A green vehicle routing problem[END_REF][START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF][START_REF] Sassi | Vehicle routing problem with mixed fleet of conventional and heterogenous electric vehicles and time dependent charging costs[END_REF][START_REF] Goeke | Routing a mixed fleet of electric and conventional vehicles[END_REF][START_REF] Hiermann | The electric fleet size and mix vehicle routing problem with time windows and recharging stations[END_REF], they do not explain how the value of β is set. It is worth noting that β plays an important role in the definition of the solution space, and therefore it restricts the optimal solution of the model. For instance, an optimal solution found with β = 3 may not be optimal for β = 4. In practice, there is no restriction on the number of times that a CS can be visited, but large values of β result in models that are computationally intractable. To overcome this difficulty, we designed an iterative procedure to solve the MILP formulation for increasing values of β. Starting with β = 0, at each iteration our procedure (i) tries to solve the MILP formulation to optimality with a time limit of 100 h, and (ii) sets β = β + 1. The procedure stops when the time limit is reached or an iteration ends with a solution s β satisfying f (s β ) = f (s * β-1 ), where f (•) denotes the objective function and * an optimal solution.

Table 1 presents the results. Since the PL approximation is the closest to reality, the results obtained using the other approximations are compared with reference to the results of the PL. For each charging function approximation, we give the objective function value (of ), the percentage gap between of and the PL solution (G), the number of routes in the solution (r), and the value of β. Since in practice the charging time is controlled by the nonlinear charging function, we evaluate the charging decisions of the L1 and L2 solutions a posteriori using the PL approximation. The last rows of Table 1 summarize the results. We present, for each approximation, the average and maximum percentage gap, the number of solutions employing more EVs than in the PL solution, and the number of infeasible solutions.

In the FS approximation EVs can charge their batteries to only around 80% of the actual capacity. Artificially constraining the capacity may force the EVs to detour to CSs more often than necessary when traveling to distant customers. Because the maximum route duration is limited, the time spent detouring and recharging the battery reduces the number of customers that can be visited. Consequently, more routes may be needed to service the same number of customers. Our results confirm this intuition: in 3 of the 20 instances the FS approximation increases the number of routes. Furthermore, in practice some distant customers may not be included in routes unless the EVs can fully use their battery capacity. In our experiments, 9 instances become infeasible with the FS approximation. In conclusion, although the FS approximation simplifies the problem (avoiding the nonlinear segment of the charging function) it may lead to solutions that are infeasible or have larger fleets and are (on average) 2.70% more expensive.

As mentioned before, L1 assumes that batteries charge faster than they do in reality (Figure 2b). As a consequence, routes based on L1 may in practice need more time to reach the planned charge levels. The extra time may make a route infeasible if there is little slack in the duration constraint. Indeed, the post-hoc evaluation shows that for 14 instances, the L1 solutions are infeasible in practice. On the other hand, L2 assumes that batteries charge slower than in reality (Figure 2c). Overestimating the charging times does not lead to feasibility issues, but the resulting routes may be overly conservative. For instance, in our experiments L2 leads to solutions that are (on average) 1.45% more expensive, and it increases the number of routes in two instances. We implemented our ILS in Java (jre V.1.8.0) and used Gurobi Optimizer (version 5.6.0) to solve the FRVCP and the set partitioning problem in the HC component. We set a time limit of 800 s in Gurobi to control the running time of the HC phase. All the experiments were run on a computing cluster with 2.33 GHz Inter Xeon E5410 processors with 16 GB of RAM running under Linux Rocks 6.1.1. The results of our ILS+HC are computed over 10 runs. Each replication of the experiments was run on a single processor.

Parameter settings

Three main parameters define the configuration of our ILS+HC: the number of iterations (K) and the methods used to solve the FRVCP in i) the split procedure (Split(•)) and ii) the GCI neighborhood (GCI(•)). We conducted a computational study to fine-tune these parameters. Table 2 sumarizes the values tested for each parameter.

To avoid overfitting, we conducted the parameter tuning on a training instance set. The latter is made up of 24 additional instances of 6 different sizes (in terms of number of customers): 10, 20, 40, 80, 160, 320. For each size we generated 4 instances. Figure 9 presents the results of our parameter tuning. The X coordinate represents K, the Y coordinate represents the CPU time, the circle diameter represents the average gap6 with respect to the BKS, and the color represents the combination of methods used to solve the FRVCP.

As expected, with higher values of K the algorithm achieves better results: between the best configuration with K = 40 and that with K = 120 there is a difference of more than 1.1% in the average gap. To simplify the discussion, we will focus on the results obtained with K = 120, but the conclusions are valid for any K. Not surprisingly, with {Split(S), GCI(S)} our ILS delivers the best results but also consumes the most CPU time. The opposite is true with {Split(H), GCI(H)}: the algorithm is fast, but it delivers poor solutions. This result highlights the importance of making optimal charging decisions when solving the E-VRP-NL (and E-VRPs in general). With {Split(H), GCI(S)} the method obtains better results than with {Split(S), GCI(H)}. This result was expected and is consistent with the notion that an aggressive local search is more important than an excellent initial solution generator. After analyzing our results, we decided that {K = 80, Split(H), GCI(S)} is the configuration that offers the best trade-off between solution accuracy and computational performance. The experiments reported in the remainder of the paper used this configuration.

Performance of the hybrid metaheuristic

Since the E-VRP-NL is a new problem, there are no results or algorithms to benchmark against. To get an idea of the quality of the solutions delivered by our ILS+HC, we compared its results with the results obtained by Gurobi using the MILP formulation introduced in Section 2. Gurobi reported integer solutions for 45 of the 120 instances; 27 of those solutions are proven optima. The instances with proven optima are the twenty 10-customer instances and seven 20-customer instances. Table 3 summarizes the results of our ILS+HC on the instances with proven optima using five metrics: the number of optimal solutions found, the average and maximum gap with respect to the optimal solution, the average best gap 7 , and the average computing time. The results suggest that our ILS+HC is able to deliver high- 7 The best gap is the gap between the best solution found over 10 runs and the optimal solution. quality solutions for the E-VRP-NL. It matched the 27 optimal solutions, and the average and maximum gaps show that our method is stable.

Table 4 summarizes the results of our ILS+HC on the 18 instances in which Gurobi found an integer solution but did not prove optimality. We use the metrics introduced above, replacing the number of optimal solutions by the number of BKSs found and computing the average and maximum gaps with respect to the BKSs. On these instances, our ILS+HC found all the BKSs (18/18) and reported an average gap of 1.09.

Since we do not have benchmark solutions for the remaining 93 instances, we do not discuss in this section the results of our ILS+HC on those instances. To allow future comparisons with our method, we report in AppendixB detailed results for each of the 120 instances. These results are also available at http://www.vrp-rep.org. 17.17

Characteristics of good E-VRP-NL solutions

We analyze in this section the characteristics of the BKSs found in our experiments. We aim to provide the reader with insight that may be useful for the design of new solution methods for the E-VRP-NL.

In total, our BKSs are made up of 1,426 routes. Our first analysis concerns the fraction of those routes that exploit mid-route charging. Figure 10 presents the percentage of routes with and without visits to CSs, grouped by instance size. The data shows that on average 71.47% of the routes in the BKSs visit at least one CS. This percentage is roughly the same for each instance size. This figure provides two insights. First, mid-route charging is a key element of good E-VRP-NL solutions (probably because it gives algorithms flexibility in assigning the customers to the routes). Second, since charging decisions affect most of the routes making up a good solution, they play a critical role in its quality.

Number of customers

Routes with CSs Routes without CSs The second analysis concerns the number of mid-route charges per route. Figures 11a and11b present histograms of the number of mid-route charges per route and the maximum number of mid-route charges per route in a solution. Figure 11a shows that for those routes performing mid-route charging, 58.58% charge once, 40.00% twice, 1.43% three times. Although a large proportion of the routes perform a single mid-route charge, 85.83% of the solutions contain at least one route with more than one mid-route charge (Figure 11b). These figures suggest that models and methods for the E-VRP-NL can benefit from relaxing the upper bound of one visit to a CS per route that is sometimes used in the literature. The third analysis concerns the energy recovered through mid-route charges. Figure 12 presents a histogram of the average battery level (in % of the total battery capacity) after a mid-route charge. The numbers show that over 90% of the mid-route charges are partial charges (i.e., they do not fully charge the battery). This figure indicates the importance of embedding components capable of making partial charging decisions into E-VRP-NL solution methods8 . A second interesting observation from the data displayed in Figure 12 is the percentage (around 12%) of mid-route charges that restore the battery to above 80% of its capacity. As mentioned in Section 2, a common assumption in the E-VRP literature is that the battery can be charged only in the linear segment of the charging curve (which ends at roughly 80% of the capacity). Our data suggests that good E-VRP-NL solutions often include routes with mid-route charges that take the battery level into the nonlinear part of the charging function.

Conclusion and future work

In this paper, we have introduced a new E-VRP that captures the nonlinear charging behavior of the charging process using a piecewise linear approximation: the electric vehicle routing problem with nonlinear charging function (E-VRP-NL). To solve the problem we proposed an ILS enhanced with HC. At the heart of our method is a neighborhood scheme that solves a new variant of the FRVCP. This problem consists in optimizing the charging decisions (where and how much to charge) of a route serving a fixed sequence of customers. We conducted an experiment comparing our charging function approximation to those in the literature. Our results show that neglecting the nonlinear nature of the charging function may lead to infeasible or overly expensive solutions. To assess the performance of our ILS+HC, we ran it on a new set of instances. We compared its solutions with the solutions obtained by Gurobi using the MILP formulation of the E-VRP-NL. The results suggest that our ILS+HC is able to deliver high-quality solutions for the E-VRP-NL. Finally, we analyzed the solutions delivered by our method to provide fellow researchers with insight into the characteristics of good E-VRP-NL solutions. Our analysis concluded that good solutions tend to use multiple mid-route charges, exploit partial recharges, and employ the nonlinear segment of the battery charging function.

Interesting research directions include designing alternative E-VRP-NL methods (both exact and heuristic) to give a better basis for comparison with our results. Another interesting possibility would be to develop approaches for the FRVCP that offer a different trade-off between accuracy and efficiency than those of the two approaches proposed in this paper. Last but not least, it would be interesting to extend the E-VRP-NL to consider capacitated CSs. To adapt our ILS+HC to this new variant two components would need to be re-designed: the GCI neighborhood and the HC. The former would need to solve the FRCVP simultaneously for all routes ensuring that CS capacity constraints are satisfied. The latter would need to guarantee that CS capacities are satisfied in the set-partitioning formulation. In ongoing research, we are working on these extensions.

Acknowledgement

The authors would like to thank the Universidad EAFIT scientific computing center (APOLO) for its support for the computational experiments. This research was partly funded in Colombia by Universidad EAFIT, Programa de Movilidad Doctoral hacia Francia (Colfuturo -Emb. de Francia -ASCUN -Colciencias -Min. de Educación), Programa Enlaza Mundos (Alcaldía de Medellín), and Universidad de Antoquia through project UdeA-2016-10705; and in France by Agence Nationale de la Recherche through project e-VRO (ANR-15-CE22-0005-01).

AppendixA. Solution approaches for the fixed-route vehicle-charging problem

As mentioned above, we propose two approaches for the FRVCP. The first solves an MILP formulation using a commercial solver, and the second solves the problem using a greedy heuristic adapted from the literature. This appendix describes these approaches.

AppendixA.1. MILP formulation running on commercial solver

The model introduced in Section 2.4 may seem too complex to be efficiently solved using out-of-the-box software. However, in practice the number of customers per route and the number of available CSs tend to be low, so the resulting MILP formulations are within the scope of commercial solvers. For instance, on a problem with 8 CSs, the MILP formulation for a fixed route serving 10 customers has 539 continuous variables, 352 integer variables, and 1,263 constraints9 . We optimally solved this problem using Gurobi Optimizer (version 5.6.0) in 0.06 s. We therefore decided to embed into our ILS+HC a component that uses a commercial solver for the FRVCP.

To further reduce the size of the MILP formulation and improve the solver's performance, we propose four preprocessing strategies that eliminate infeasible CS insertions. Our strategies rely on the following two premises: (i) the energy consumption and the travel time between vertices satisfy the triangular inequality; and (ii) since the piecewise-linear charging function is concave (i.e., ρ j,k-1 ≥ ρ jk ), the first segment has the fastest charging rate.

We propose four strategies. The first three filter CS insertions that are infeasible independently of how the customers are sequenced in the routes. These strategies are applied only once before running the ILS+HC. The fourth strategy filters CS insertions that are infeasible for a specific fixed route.

Strategy 1: This strategy estimates the minimum time τ needed to visit CS j ∈ F between two vertices i and h ∈ I ∪ {0}. This time is defined as the sum of a lower bound on the travel time (u) and a lower bound on the charging time (v) of any route serving customers i and h. Note that in any route serving customer i, t 0i is a lower bound on the route's duration from the start to customer i. Similarly, note that in any route serving customer h, t h0 is a lower bound on the time to complete a route from vertex h. From these two observations we can compute the minimum duration of a route visiting CS j between vertices i and h as u = t 0i + t ij + t jh + t h0 + p i + p h . To compute the minimum charging time v, we need to compute the minimum amount of energy (ec) that an EV coming from vertex i and traveling to vertex h must charge at CS j. This amount is the charge needed to recover the energy consumed to make the detour to j, i.e., ec = e ij + e jh -e ih . Because the battery level when the EV arrives at i in any E-VRP-NL solution is unknown a priori, we consider that the battery is charged at j using the fastest charging rate (ρ 0j ). Then v = ec ρ 0j . It is clear that if τ = u + v > T max then any route visiting CS j between customers i and h is infeasible. We therefore forbid this insertion in our MILP.

Strategy 2: This strategy computes a lower bound on the remaining energy that an EV must have on arrival at customer i to be able to visit CS j next. Note that in terms of energy remaining, the best way to reach vertex i is to visit it immediately after fully charging at CS c(i) = arg min l∈F ∪{0} e li . If Q -e c(i)i < e ij any route visiting j after i is energy-infeasible and we can safely forbid this insertion in our MILP.

Strategy 3: Note that o = Q -e ji is a lower bound on the energy remaining when an EV arrives at customer i immediately after charging at CS j. Note also that o must be enough to at least close the route (reach the depot) or reach the closest CS in terms of energy consumption. If o < e ic(i) , where c(i) = arg min l∈F ∪{0} e il , any route visiting j before i is energy-infeasible and we can safely forbid this insertion.

Strategy 4: This strategy estimates a lower bound on the new duration of a given fixed route if CS j ∈ F is inserted between vertices π(i) and π(i + 1) ∈ Π, i = n r . This bound, t , is defined as the sum of the new travel time (u) and a lower bound on the charging time (v). It is easy to see that u = t + t π(i)j + t jπ(i+1) -t π(i)π(i+1) . Similarly to Strategy 1, to estimate v, we consider that the battery is charged at j using the fastest charging rate. Therefore, v = ec ρ 0j , where ec = e π(i)j + e jπ(i+1) -e π(i)π(i+1) is the charge needed to recover the energy consumed in the detour to j. If t = (u + v) > T max , inserting j between vertices π(i) and π(i + 1) leads to an infeasible route, so we can safely forbid this insertion.

Applying our preprocessing strategies to the MILP formulation for the 10-customer route of the example above, we reduce the model to 71 continuous variables, 40 integer variables, and 165 constraints. We solved the model in Gurobi in 0.02 s.

AppendixA.2. Greedy heuristic

Existing metaheuristics for E-VRPs use various approaches to make charging decisions. One popular approach is the recharge relocation operator proposed by [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF] for the green vehicle routing problem with multiple technologies and partial recharges (GVRP-MTPR). This approach considers the insertion of only one CS per route. Starting from an energy-feasible fixed-route the procedure first deletes the current CS. Then, it tries to improve the charging decisions by inserting each CS into each arc of the fixed-route. To decide how much energy to charge at the inserted CS, the algorithm applies a simple rule: charge the minimum amount of energy needed to reach the depot (i.e., to complete the route). We propose here a heuristic to solve the FRVCP based on this approach.

Our heuristic has two phases: location of CSs and setting of the charge. In the first phase, it iteratively inserts CSs into the arcs of the fixed-route Π to ensure feasibility in terms of energy. In the second phase, it improves the charging decisions by adjusting the energy charged at each visited CS. Algorithm 1 describes our heuristic. The approach uses four important procedures: trackBattery(•), sumNegative(•), totalTime(•), and copyAndInsert(•). Procedure trackBattery(•) computes the battery level Y i at each vertex i ∈ Π, assuming that the EV fully charges its battery at each visited CS. Note that Y i may take negative values. Procedure sumNegative(•) computes the sum of the battery levels with negative values (i.e., s = i∈Π min{0, Y i }). Procedure totalTime(•) computes the total time t of the route (assuming a full charging policy). Finally, procedure copyAndInsert(•) takes as input a fixed route, a CS, and a position in the route; it returns a copy of the fixed route with the CS inserted at the given position.

The heuristic starts with the location phase (lines 2-27). After computing s for the current fixed-route Π, it enters the outer loop (lines 7-27). In each pass through the inner loop (lines 8-25), the heuristic i) evaluates the insertion of a CS into each arc of Π assuming that the EV fully charges its battery, and ii) selects the insertion that maximizes s (lines 14-18). If s = 0 (i.e., the route is energy-feasible), it selects the insertion that minimizes t (lines 19-23). Then, it performs the selected insertion (line 26). If the route Π is still energy-infeasible (i.e., s < 0), the heuristic starts again at line 8 and tries to insert additional CSs until feasibility in terms of energy is reached.

In the charge-setting phase (line 28), the heuristic invokes procedure ruleMinEnergy(Π) to set the energy charged at each CS according to the [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF] rule. Finally, it evaluates if the route satisfies the maximum-duration constraint (lines 29-33). It returns a boolean variable indicating whether or not the fixed-route is feasible (f ) and the route Π with the newly inserted CSs. if s = 0 and t < t then B.5 gives the results reported by our ILS+HC on the E-VRP-PNL instances. We compare our results with the results obtained by Gurobi running the MILP introduced in Section 2. For each instance, we give the problem name 10 and the BKS reported by Gurobi or ILS+HC.

For the results obtained with Gurobi, we report the best solution (Best) and the gap with respect to the BKS (G). For the results obtained with the ILS+HC, we report the best solution, the average solution (Avg.), and the average computing time (CPU) in seconds over ten runs. For the ILS+HC, we provide the gap of the average solution and best solution with reference to the BKS. The last rows of the table summarize the average and maximum gaps with respect to the BKSs, the number of times each method found the BKS, and the average and maximum running times. Values in bold indicate that a method found the BKS.

Figure 1 :

 1 Figure 1: Typical charging curve, where i, u, and SoC represent the current, terminal voltage, and state of charge respectively. The SoC is equivalent to the battery level. (Source: Hõimoja et al. 2012).

Figure 2 :

 2 Figure 2: Linear approximations in the literature vs. real data provided by Uhrig et al. (2015) 2 .

ChargingFigure 4 :

 4 Figure 4: Example of a feasible E-VRP-NL solution.

Fixed

  Figure 5: Example of a fixed-route vehicle-charging problem.

Figure 6 :

 6 Figure 6: Splitting a TSP tour into an E-VRP-NL solution.

Figure 7 :

 7 Figure 7: Perturbation procedure example using double bridge operator.

Figure 9 :

 9 Figure 9: Results of the parameter tuning: testing 20 possible configurations.

Figure 10 :

 10 Figure 10: Percentage of the routes with/without visits to CSs by instance size.

  of the maximum number of visited CSs in the routes of each solution.

Figure 11 :

 11 Figure 11: Analysis of the number of visits to CSs.

Figure 12 :

 12 Figure 12: Histogram of the average battery level (in % of total battery capacity) after a mid-route charge.

Table 1 :

 1 Comparison of our charging function approximation (piecewise linear approximation) with charging function approximations from the literature

	Instance	PL	FS			L1	L2
						Solution	Evaluation	Solution	Evaluation
		of of of r r r β β β	of of of	r r r β β β	of of of G(%) G(%) G(%) r r r β β β of of of	r r r β β β	of of of G(%) G(%) G(%) r r r β β β
	tc0c10s2cf1	19.75 3 2	NFS NFS NFS NFS 19.61 3 2 NFE NFE NFE NFE 20.50 3 2 20.22 2.38 3 2
	tc0c10s2ct1	12.30 2 0	12.61 2.52 3	0 12.22 2 0 12.42 0.98 2	0 12.46 2 0 12.30 0.00 2 0
	tc0c10s3cf1	19.75 3 2	NFS NFS NFS NFS 19.61 3 2 NFE NFE NFE NFE 20.50 3 2 20.22 2.38 3 2
	tc0c10s3ct1	10.80 2 0 10.80 0.00 2	0 10.79 2 0 11.03 2.13 2	0 10.97 2 0 10.80 0.00 2 0
	tc1c10s2cf2	9.03 3 0	9.03 0.00 3	0	9.03 3 0 9.12 1.00 3	0	9.14 3 0 9.03 0.00 3 0
	tc1c10s2cf3	16.37 3 2	NFS NFS NFS NFS 15.99 3 1 NFE NFE NFE NFE 16.89 3 2 16.37 0.00 3 2
	tc1c10s2cf4	16.10 3 2	NFS NFS NFS NFS 15.66 3 2 NFE NFE NFE NFE 16.43 3 2 16.23 0.81 3 2
	tc1c10s2ct2	10.75 3 1 10.75 0.00 3	1 10.75 3 0 10.76 0.09 3	0 10.94 3 0 10.78 0.28 3 0
	tc1c10s2ct3	13.17 2 2	15.98 21.34 3	2 13.06 2 2 NFE NFE NFE NFE 13.60 2 2 13.17 0.00 2 2
	tc1c10s2ct4	13.83 2 1	NFS NFS NFS NFS 13.34 2 1 NFE NFE NFE NFE 14.17 2 1 14.17 2.46 2 1
	tc1c10s3cf2	9.03 3 0	9.03 0.00 3	0	9.03 3 0 9.12 1.00 3	0	9.14 3 0 9.03 0.00 3 0
	tc1c10s3cf3	16.37 3 1	NFS NFS NFS NFS 15.99 3 1 NFE NFE NFE NFE 16.89 3 2 16.37 0.00 3 2
	tc1c10s3cf4	14.90 3 1	NFS NFS NFS NFS 14.56 2 1 NFE NFE NFE NFE 15.18 3 0 15.18 1.88 3 0
	tc1c10s3ct2	9.20 3 0	9.20 0.00 3	0	9.19 3 0 NFE NFE NFE NFE 10.80 3 0 10.57 14.89 3 0
	tc1c10s3ct3	13.02 2 0	13.07 0.38 2	1 12.98 2 0 13.16 1.08 2	0 13.60 2 0 13.02 0.00 2 0
	tc1c10s3ct4	13.21 2 0	13.58 2.80 3	1 12.92 2 1 NFE NFE NFE NFE 13.71 2 0 13.21 0.00 2 0
	tc2c10s2cf0	21.77 3 3	NFS NFS NFS NFS 14.53 2 2 NFE NFE NFE NFE 22.78 4 4 22.15 1.75 4 4
	tc2c10s2ct0	12.45 3 2 12.45 0.00 3	2 12.44 3 3 NFE NFE NFE NFE 12.93 3 2 12.45 0.00 3 2
	tc2c10s3cf0	21.77 3 2	NFS NFS NFS NFS 14.53 2 2 NFE NFE NFE NFE 23.02 4 3 22.20 1.98 4 3
	tc2c10s3ct0	11.51 3 0 11.51 0.00 3	0 11.50 3 0 NFE NFE NFE NFE 11.92 3 0 11.54 0.26 3 0
	Avg. Difference (%)	2.70			1.04	1.45
	Max. Difference (%)	21.34			2.13	14.89
	Solutions with larger fleet		3			0	2
	Infeasible solutions	9			14		0
	NFS: Infeasible solution, NFE: Infeasible evaluation
	G(%) = (of -ofP L)/ofP L × 100
	4.3. Results of our ILS+HC on E-VRP-NL instances
	4.3.1. Experimental environment

Table 2 :

 2 Values of the parameters evaluated in the experiment

	Parameter	Values
	K	{40, 60, 80, 100, 120}
	Split(•)	{S, H}
	GCI (•)	{S, H}
	S stands for commercial solver (see AppendixA.1)
	H stands for greedy heuristic (see AppendixA.2)

Table 3 :

 3 Summary results of our ILS+HC on instances with proven optima

	Metric	ILS + HC
	Number of optimal solutions	27/27
	Avg. Gap (%)	0.29
	Max. Gap (%)	1.87
	Avg. Best Gap (%)	0.00
	Avg. Time (s)	6.31

Table 4 :

 4 Summary results of our ILS+HC on instances without proven optima

	Metric	ILS + HC
	Number of BKSs	18/18
	Avg. Gap (%)	1.09
	Max. Gap (%)	2.95
	Avg. Best Gap (%)	0.00
	Avg. Time (s)

Table B

 B

		.5 -continued from previous page Table B.5 -continued from previous page
			Gurobi Gurobi		ILS+HC ILS+HC
	Instance Instance	BKS Best G (%) Best G (%) Avg. G (%) t (s) BKS Best G (%) Best G (%) Avg. G (%) t (s)
	tc0c20s3ct2 tc1c20s3ct3 Avg. Gap tc1c20s3ct4 Max. Gap tc2c20s4cf0 Found solution tc1c20s4cf1 Best tc0c20s4cf2 Avg. Time tc1c20s4cf3 Max. Time	17.08 17.11 12.65 12.68 16.21 16.21 24.67 25.36 45 16.39 16.40 27 27.48 -16.56 16.80	0.18 17.08 0.24 12.65 2.61 0.00 16.21 41.36 2.80 24.67 120 0.06 16.39 120 -27.48 1.45 16.56	0.00 17.13 0.00 12.72 0.00 0.00 16.25 0.00 0.00 24.69 0.00 16.40 0.00 27.61 0.00 16.80	0.29 0.55 1.32 0.25 4.44 0.08 0.06 0.47 1.45	8.49 8.86 5.16 14.63 13.47 849.55 12.81 8.69 7636.50
	tc1c20s4cf4	17.00 17.00	0.00 17.00	0.00 17.00	0.00	4.17
	tc2c20s4ct0	26.02	-	-26.02	0.00 26.02	0.00	15.25
	tc1c20s4ct1	18.25 18.25	0.00 18.25	0.00 18.32	0.38	16.14
	tc0c20s4ct2	16.99 17.21	1.29 16.99	0.00 17.10	0.65	9.33
	tc1c20s4ct3	14.43 14.43	0.00 14.43	0.00	14.5	0.49	7.99
	tc1c20s4ct4	17.00 17.00	0.00 17.00	0.00 17.00	0.00	6.08
	tc0c40s5cf0	32.67	-	-32.67	0.00 33.25	1.78	23.85
	tc1c40s5cf1	65.16	-	-65.16	0.00 66.03	1.34	44.01
	tc2c40s5cf2	27.54 38.93 41.36 27.54	0.00 27.67	0.47	31.64
	tc2c40s5cf3	19.74 21.04	6.59 19.74	0.00 20.18	2.23	16.85
	tc0c40s5cf4	30.77 36.47 18.52 30.77	0.00 31.49	2.34	33.33
	tc0c40s5ct0	28.72	-	-28.72	0.00 29.35	2.19	24.50
	tc1c40s5ct1	52.68	-	-52.68	0.00 53.36	1.29	58.52
	tc2c40s5ct2	26.91	-	-26.91	0.00 27.02	0.41	22.85
	tc2c40s5ct3	23.54	-	-23.54	0.00 23.77	0.98	26.48
	tc0c40s5ct4	28.63	-	-28.63	0.00 28.72	0.31	32.55
	tc0c40s8cf0	31.28	-	-31.28	0.00 32.02	2.37	33.59
	tc1c40s8cf1	40.75	-	-40.75	0.00 42.33	3.88	69.99
	tc2c40s8cf2	27.15 29.19	7.51 27.15	0.00 27.31	0.59	28.92
	tc2c40s8cf3	19.66 22.01 11.95 19.66	0.00 20.24	2.95	19.46
	tc0c40s8cf4	29.32	-	-29.32	0.00 29.86	1.84	43.05
	tc0c40s8ct0	26.35 30.29 14.95 26.35	0.00 26.89	2.05	28.54
	tc1c40s8ct1	40.56	-	-40.56	0.00 41.19	1.55	70.50
	tc2c40s8ct2	26.33	-	-26.33	0.00 26.71	1.44	25.64
	tc2c40s8ct3	22.71 23.51	3.52 22.71	0.00 23.23	2.29	25.25
	tc0c40s8ct4	29.20	-	-29.20	0.00 29.27	0.24	47.46
	tc0c80s8cf0	39.43	-	-39.43	0.00 39.86	1.09	56.41
	tc0c80s8cf1	45.23	-	-45.23	0.00 45.73	1.11 121.27
	tc1c80s8cf2	30.81	-	-30.81	0.00 31.83	3.31	50.99
	tc2c80s8cf3	32.44	-	-32.44	0.00 32.60	0.49	64.05
	tc2c80s8cf4	49.29	-	-49.29	0.00 49.69	0.81	99.84
	tc0c80s8ct0	41.90	-	-41.90	0.00 42.76	2.05	54.35
	tc0c80s8ct1	45.27	-	-45.27	0.00 45.85	1.28 129.66
	tc1c80s8ct2	31.74	-	-31.74	0.00 32.36	1.95	59.73
	tc2c80s8ct3	32.31	-	-32.31	0.00 32.55	0.74	65.15
	tc2c80s8ct4	44.83	-	-44.83	0.00 46.61	3.97 111.24
	tc0c80s12cf0	34.64	-	-34.64	0.00 35.59	2.74	57.24
	tc0c80s12cf1	42.90	-	-42.90	0.00 44.07	2.73	74.58
	tc1c80s12cf2	29.54	-	-29.54	0.00 30.73	4.03	61.34
	tc2c80s12cf3	31.97	-	-31.97	0.00 32.70	2.28	75.64
						Continued on next page

http://www.avere-france.org/Site/Article/?article_id=5644. Last accessed 11/16/2016.

[START_REF] Uhrig | E-mobility in car parks-guidelines for charging infrastructure expansion planning and operation based on stochastic simulations[END_REF] conducted experiments to estimate the charging time for different charge levels with two types of EVs and three types of CSs.

The optimal solution to an FRVCP solved over an energy-feasible fixed route is trivial: the original fixed route itself.

In our implementation we use intra-route and inter-route versions with best-improvement selection.

The instances will be made available after the completion of the reviewing process.

G(%) = (of -of BKS )/of BKS × 100

Note that until 2016 this was the exception rather than the rule in the E-VRP literature.

We randomly picked the route from the best solution found for a randomly selected instance.

10 tcαcβsµc #, where α α α is the method used to place the customers (i.e., 0: randomization, 1: clustering, 2: mixture of both), β β β is the number of customers, µ µ µ is the number of the CSs, is 't' if we use a p-median heuristic to locate the CSs and 'f' otherwise, and # # # is the number of the instance for each combination of parameters (i.e., # = 0, 1, 2, 3, 4).