David L Cortés-Murcia
email: david.cortesmurcia@utt.fr

Caroline Prodhon
email: caroline.prodhon@utt.fr

H Murat Afsar
email: murat.afsar@utt.fr

The electric vehicle routing problem with time windows, partial recharges and satellite customers

Keywords: electric vehicles, vehicle routing, iterated local search

In this paper, a new variant of the electric vehicle routing problem is presented. Since recharging time is commonly considered as idle time, the aim is to take advantage of it by allowing customer visits by an alternative mode of transport while the electric vehicle is at a recharging station. This is particularly pertinent in a city logistics context. A mathematical model is proposed, as well as an Iterated Local Search metaheuristic. The Iterated Local Search is reinforced by adding Variable Neighborhood Descent and set partitioning. The proposed method is tested on public instances for E-VRPTWPR. Finally, it is shown that allowing satellite customers enables us to take advantage of recharging times and to reduce time spent at recharging stations.

Introduction

In present-day logistics, reducing negative impacts on society and on the environment is a difficult challenge for companies. According to the European Commission (2016), about 25% of Europe's greenhouse gas (GHG) emissions are generated by transport activities. In addition, the main difficulty lies in the fact that most of the fleet used by companies, is composed of conventional internal combustion vehicles (ICVs). The ICVs contribute to air pollution, because of the emissions of carbon dioxide (CO 2 ), nitrogen oxides (NO x ), elemental carbon and organic carbon; and contribute to noise pollution as well. As reported by [START_REF] Nesterova | State of the art of the electric freight vehicles implementation in city logistics -Update 2015[END_REF], since the World Health Organization stated that poor air quality is a serious health risk, one of the major short term concerns for local authorities is to improve local air quality. Subsequently, the Electric Vehicles (EVs) get more and more attention from companies and governments because they can play an important role in achieving this goal. For a number of years now, there have been a growing number of projects and national initiatives by European countries aiming to facilitate the use of EVs. Likewise, there are private initiatives by companies like UPS, La Poste and TNT Express who have included EVs in their delivery operations [START_REF] Nesterova | State of the art of the electric freight vehicles implementation in city logistics[END_REF].

On the one hand, the usage of EVs leads to some advantages [START_REF] Schiffer | The electric location routing problem with time windows and partial recharging[END_REF]. It represents cleaner transport because they have no local GHG emissions and produce minimal noise in comparison to ICVs. Moreover, a zero emission balance can be obtained if renewable energy sources are used to generate the electricity. In addition, the operational costs and the maintenance costs for EVs have been reported to be lower than those of ICVs (Davis and Figliozzi, 2013;[START_REF] Pelletier | 50th anniversary invited articlegoods distribution with electric vehicles: review and research perspectives[END_REF]. On the other hand, it is essential to address some disadvantages. For example, acquisition cost are higher for EVs than for ICVs. Also, limited driving range and long charging time affect operational decisions. In addition, the insufficiency of charging infrastructure is considered an obstacle. Despite that, due to cutting-edge technology, governmental subsidies, supportive policies, and regulations for GHG Emissions, EVs are now competitive and are an alternative to ICVs (Davis and Figliozzi, 2013;[START_REF] Nesterova | State of the art of the electric freight vehicles implementation in city logistics -Update 2015[END_REF].

To deal with the effects on the operational decisions, studies on Vehicle Routing Problems (VRP) consider the increase in EVs usage (Afroditi et al., 2014). First, Erdogan and Miller-Hooks (2012) proposed a Green Vehicle Routing Problem (G-VRP). In this paper, the fleet is composed of Alternative Fuel Vehicles (AFV). Later, the Electric Vehicle Routing Problem with Time Windows (E-VRPTW) was presented by [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF]. In this case, fleet is composed exclusively of EVs and customers have time windows. Subsequently, the E-VRP is studied by different authors who consider: mixed fleets with electric and non-electric vehicles, linear and non-linear recharging models, linear and non-linear energy consumption models and decisions related to location of recharging stations (RSs). A brief discussion of the related literature is presented in Section 2. In general, all these studies intend to provide decision tools, and to encourage the transition to green energies, especially to EVs.

Even so, the intra-route charging time is often considered as a downtime because charging times could be so long that they affect the total operational time in comparison with the ICVs routing. Thus, the intention of our paper is to propose a novel variant of the E-VRP which aims to take advantage of the intra-route charging time. In this variant, we allow customer visits by an alternative means of transport while the EV is at a RS. In other words, for each visit to a RS, it is possible to visit one customer while the vehicle is in the charging process. These customers are called satellite customers. For simplicity, in this document we call the alternative mode as by walking but it could represent any type of alternative vehicle (bikes, drones, segways, etc.). In addition, partial recharges, time windows, vehicle capacity, battery capacity and customer demand are considered.

In general, RSs are managed as parking lots because of the long recharging times. Thus, the operative cost of visiting a RS is given by the amount of energy recharged but also by the time a parking lot is occupied. There are even car parks with RS where the charging fee is according to time. In the proposed problem, the parking time can be longer than the recharging time because of satellite customers. In the proposed model, with respect to recharge, it is assumed that the state of charge (SoC) linearly increases in accordance with a recharge rate and the time spent at RS. Therefore, the objective function is the minimization of the total time spent at the RSs. This objective implies the reduction of the total distance given that recharging time depends on energy requirement and that, in return, depends on the distance travelled. Additionally, it is considered that going by an alternative mode is commonly slower than going with the EV. So, this objective function avoids undesirable scenarios (for practitioners) where, because of reducing total distance, walking times increases and the route duration is affected.

To sum up, the proposed variant allows the intra-route charging time to become an advantage and to reduce the total distance performed by the EVs. This variant is more pertinent for small-package shipping or no capacity service industries in a city logistics context. The variant is called the Electric Vehicle Routing Problem with Time Windows and Satellite Customers (E-VRPTWsc).

The paper is structured as follows: in Section 2, a review of related literature is discussed. In Section 3, the Electric Vehicle Routing Problem with Time Windows and Satellite Customers (E-VRPTWsc) and the corresponding mathematical model are presented. An Hybrid-ILS meta-heuristic is proposed as a solution method and is described in Section 4. The test instances proposed by [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF] are adapted and computational results are shown in Section 5. Also, a comparison with the best known solutions (BKS) of the E-VRPTW with Partial Recharging (E-VRPTWPR) and the benefit of developing this variant is highlighted. Finally, the conclusions of the paper and an overview of future research is given in Section 6.

Literature review

A short state of the art related to our problem is presented. Interested readers are referred to [START_REF] Pelletier | 50th anniversary invited articlegoods distribution with electric vehicles: review and research perspectives[END_REF] and [START_REF] Schiffer | Vehicle routing and location routing with intermediate stops: A review[END_REF]. Those surveys provide an overview of transportation with EV. In this section we first review the papers on vehicle routing with alternative fuel. Second, we refer to E-VRP and its variants.

Conrad and Figliozzi (2011) extend the CVRP-TW by including vehicles with limited range. The vehicles are allowed to recharge at customer locations to extend the driving range. Two objective functions are proposed: the minimization of the number of routes, and the minimization of total cost related to the travel distance, service time and vehicle recharging. Two recharging policies are proposed: full charge and partial charge equivalent to 80% of full battery. As solution method the authors propose a heuristic based on an iterative construction and improvement algorithm for the CVRP-TW.

Soon afterward, one of the first papers where recharging points were modeled as dedicated points is in the Green Vehicle Routing Problem (G-VRP) considered by Erdogan and Miller-Hooks (2012). These points have to be visited to extend vehicle's driving range. As mentioned before, in G-VRP the fleet is composed of AFVs and the objective is to minimize the total distance. The authors use full refueling policy. They introduce two constructive heuristics, the Modified Clarke and Wright Savings heuristic and the Density-Based Clustering Algorithm. Other works address the green VRP: Felipe et al. ( 2014) who present several heuristics which are used within a non deterministic Simulated Annealing framework; and [START_REF] Montoya | A multi-space sampling heuristic for the green vehicle routing problem[END_REF] who design a competitive two-phase heuristic.

To the best of our knowledge, the first reference to the electric vehicle routing problems can be found in [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF]. The authors introduce the E-VRPTW. It is devoted to managing EVs, visits to recharging stations and customer time windows. They propose a hierarchical objective function which aims to minimize the number of vehicles and the total traveled distance. The recharging time is introduced as a linear function and full charge is mandatory. The energy consumption is proportional to the distance travelled. They present a hybrid heuristic that combines a tabu search heuristic with a variable neighborhood search algorithm. Furthermore, the test instances based on [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF], are developed. Later, distance, minimizing number of vehicles, number of charging stations sited, and minimizing total cost as well. A mathematical model is presented and a set of small instances is solved using commercial solvers. Later, [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF] propose the Location Routing Problem with Intraroute Facilities which is a generalization of the ELRP-TWPR. Large instances of the ELRP-TWPR are solved by using an ALNS which is enhanced by local search and dynamic programing components.

Recently, a problem focused on reducing the inefficiency of RS in a LRP context can be found in Schiffer et al. (2018). In the problem it is considered that at the intermediate stops either energy or freight or both can be replenished. As solution method they apply an ALNS. New instances are proposed and the performance of the ALNS is tested by solving instances of the ELRP-TWPR.

Hiermann et al. ( 2019) introduce an E-VRP combining conventional, plug-in hybrid, and electric vehicles. They use a metaheuristic approach based on a genetic algorithm hybridized with an integer programming. They prove the competitiveness of their algorithm solving E-VRPTW and E-VRPTWPR instances and some new best known solutions are provided.

All the previous related works, including our proposed problem, assume that the EV batteries are charged in a conductive way. So, it means that the vehicle is plugged into a electric power source and it must wait until the battery is sufficiently charged. However, there is another research field which works with battery swapping as a charging system. [START_REF] Yang | Battery swap station location-routing problem with capacitated electric vehicles[END_REF] presents an electric vehicle battery swap station location routing problem (BSS-EV-LRP). The problem includes location decisions of the battery swapping stations from a set of candidate places. Later, Hof et al. ( 2017) uses an Adaptive Variable Neighborhood Search algorithm to solve the BSS-EV-LRP. They improve the results presented by [START_REF] Yang | Battery swap station location-routing problem with capacitated electric vehicles[END_REF]. Another problem working with battery swapping stations is the Dial-a-Ride Problem with Electric Vehicles and battery swapping stations (DARP-EV) presented by Masmoudi et al. (2018). It consists of designing vehicle routes and schedules to serve a set of pre-specified transport requests during a certain planning horizon using EVs. A variant of DARP is presented by Bongiovanni et al. (2019). They introduce the electric Autonomous Dial-a-Ride Problem which considers the use of electric autonomous vehicles (e-ADARP). Even though in this variant EVs are charged in a conductive way, autonomous vehicles can operate on a non-stop schedule. Thus, the service is not limited by the driver's shift and charging times do not represent a driver idle time.

The use of battery swapping stations significantly reduces the recharging time, and it is considered as an option for reducing the idle time of the conductive recharge. However, battery swapping requires vehicles using automated battery stations and it presents technical barriers that complicate its implementation. Some technical barriers include huge infrastructure cost for the swapping stations, space and cost associated with large stock of batteries, standardization of vehicle batteries, and the risk of battery damage due to excessive swapping [START_REF] Pelletier | 50th anniversary invited articlegoods distribution with electric vehicles: review and research perspectives[END_REF].

As far as we know, none of the E-VRPs previously presented cover all aspects of our problem. Most of the existing E-VRPs do not propose strategies to reduce the impact of the long recharging times associated with the intraroute stops. Thus, we propose a routing problem which takes E-VRPTW constraints, partial recharges as well as the possibility of serving a customer during the recharging operation.

Problem statement

The Electric Vehicle Routing Problem with Time Windows and Satellite Customers (E-VRPTWsc) is defined on a complete directed graph G = (V , A) with a set of vertices V = {V ∪ F ∪ {0, N + 1}} and a set of arcs given by A = {(i, j)|i, j ∈ V , i = j}. Let V = {1, ..., N } be the set of customers, F = {0, ..., M } be the set of recharging stations and F be the set including dummy vertices that represent the multiple visits to vertices on F . Vertices 0 and N + 1 denote the depot as departing and arriving node respectively. There is a fixed fleet size of P homogeneous EVs with a cargo load capacity of Q and a battery capacity B. Each customer i has an associated demand q i , a service time s i and a hard time window [e i , l i ]. Energy consumption of an EV traveling through an arc (i, j) ∈ A is determined by the arc distance d ij and the consumption rate γ. The E-VRPTWsc seeks to find at most P routes departing from the depot visiting all customers and coming back to the depot. The routes include the stops at the RS and the customers visited by walking when needed. For each RS visit, at most one customer can be visited by walking. The objective function aims to minimize the total time spent at the RS. This time is composed of: the time to set up the fleet at the beginning of the day and the time spent at RSs during the operation. First, it is assumed that EVs are discharged at beginning of the day. Thus, set up time depends on the initial SoC required for each vehicle to perform the route and the inverse recharging rate r. The positive variable IC i establishes the initial SoC that is required by the vehicle that departs from the depot and arrives at node i ∈ V .

Second, the time spent at the RS i during the intraroute charging is represented by the positive variable T RT i , and it depends on the recharging time and the walking time. The recharging time t i is assumed linear and it depends both on the inverse recharging rate r and the difference between the required energy and SoC on arrival at i ∈ F . The walking time t ij is computed when node j ∈ V is assigned as satellite customer departing from RS i ∈ F . In this situation, t ij takes into account the walking speed ratio α, s j and [e j , l j ]. Finally, T RT i is computed as max{t i , j∈V t ij } for each node i ∈ F .

Thus, the problem can be formulated as a mixed-integer linear program (MILP). For every arc (i, j) ∈ A the boolean decision variable x ij is equal to 1 if arc (i, j) is traversed, 0 otherwise. Moreover, for the set of arcs {(i, j)|i ∈ F , j ∈ V } the boolean decision variable z ij is defined. It is equal to 1 if the customer j is visited by walking from recharging station i, 0 otherwise. Variable u i defines the remaining cargo, τ i defines the arrival time and y i defines the SoC on arrival at vertex i ∈ V . Variable w i is the amount of energy recharged at recharging station i ∈ F . Finally, the recharging time t i is computed as rw i . The walking time t ij is computed as τ j +s j +αt ji -τ i . Table 1 summarizes sets, variables, and parameters of the model.

Using this notation, the E-VRPTWsc can be formulated as the following integer program:

min i∈V rIC i + i∈F T RT i (1)
Subject to:

j∈V ,i =j x ij + j∈F z ji = 1 ∀i ∈ V (2) j∈V ,i =j x ij ≤ 1 ∀i ∈ F (3) j∈V ,i =j x ij - j∈V ,i =j x ji = 0 ∀i ∈ V (4) i∈V x 0i ≤ P (5) u i -q i x ij + Q(1 -x ij ) ≥ u j ∀i ∈ V , ∀j ∈ V , i = j (6) u i -(q k z ik ) + Q(2 -x ij -z ik ) ≥ u j ∀i ∈ F , ∀j ∈ V , ∀k ∈ V, i = j, j = k (7) u i ≤ Q ∀i ∈ V (8) rw i ≤ T RT i ∀i ∈ F (9) τ j + (s j + αt ji )z ij -l N +1 (1 -z ij ) -τ i ≤ T RT i ∀i ∈ F , ∀j ∈ V (10) y j + cd 0j -B(1 -x 0j ) ≤ IC j ∀j ∈ V (11) τ i + (s i + t ij )x ij -l N +1 (1 -x ij ) ≤ τ j ∀i ∈ V , ∀j ∈ V , i = j (12) τ i + αt ij z ij -l N +1 (1 -z ij ) ≤ τ j ∀i ∈ F , ∀j ∈ V ( 13
)
τ i + rw i + t ij x ij -(l N +1 + rB)(1 -x ij ) ≤ τ j ∀i ∈ F , ∀j ∈ V , i = j (14) τ k + (αt ik + s k )z ik + t ij x ij -(l N +1 + rB)(2 -x ij -z ik ) ≤ τ j ∀i ∈ F , ∀j ∈ V , ∀k ∈ V, i = j, j = k (15) e i ≤ τ i ≤ l i ∀i ∈ V ( 16
)
y i -γd ij x ij + B(1 -x ij ) ≥ y j ∀i ∈ V, ∀j ∈ V , i = j ( 17
)
y i + w i -γd ij x ij + B(1 -x ij ) ≥ y j ∀i ∈ F , ∀j ∈ V , i = j (18) 0 ≤ y i ≤ B ∀i ∈ V (19) 0 ≤ y i + w i ≤ B ∀i ∈ F (20) i∈V z ji ≤ 1 ∀i ∈ F (21) i∈V ,i =j x ij ≥ i∈V z ji ∀j ∈ F (22) x ij ∈ {0, 1} ∀i, j ∈ V , i = j (23) z ij ∈ {0, 1} ∀i ∈ F , j ∈ V (24) u i , τ i , IC i , y i ≥ 0 ∀i ∈ V (25) T RT i , w i ≥ 0 ∀i ∈ F (26)
The objective function (1) minimizes the total recharging time at depot and at recharging stations. Constraints (2) state that all customers have to be visited once while recharging stations and dummy vertices cannot be visited more than once, constraints (3). Flow conservation constraints (4) guarantee for each vertex that the number of incoming arcs is equal to the number of outgoing arcs. Fleet size is considered by constraints (5). Load flow and fulfillment of demand are represented by constraints ( 6) for all the vertices and by constraints (7) for those vertices visited after visiting a recharging station, where demand of a customer visited by walking affects the remaining load. Vehicle capacity is restricted by constraints (8). Total time spent at recharging stations is governed by t i represented by constraints (9), and by t ij defined by constraints (10). The SoC required for each route at depot is represented in constraints (11). For the arrival times at each vertex three situations have to be considered: constraints (12) link arrival time for a vertex j in a vehicle route, constraints (13) link arrival time for a customer j visited by walking while constraints ( 14) and ( 15) link arrival time for a vertex j visited after a recharging station i. Time windows are respected by constraints (16).

Battery level at a vertex following a customer visit is set by constraints (17) while constraints (18) set battery level at vertex after a recharging station visit. Battery capacity is guaranteed by constraints ( 19) and (20). Constraints (21) establish the maximum number of satellite customers that could be attended from each visit at a recharging station. Relation between variables x and z is represented in constraints ( 22). Finally, ( 23)-( 26) are the definition domain of variables.

Solution method

The E-VRPTWsc is an extension of the VRP and belongs to NP-hard problems. Our aim is to solve large instances in reasonable times, and it is well known that the use of metaheuristics is an efficient way to handle these type of problems. Accordingly, our solution method is based on an Iterated Local Search (ILS) metaheuristic framework.

Algorithm 1 Hybrid ILS

1: procedure Hybrid ILS(MaxILSIter,NbNoImp)

2: S ← 0; 3: best cost ← ∞; 4: i ← 0, j ← 0; 5: Ω ← ∅; 6:
S ← RandomizeParallelInsertion(); Ω ← Ω ∪ S;

11: i ← i + 1, j ← j + 1; 12:
if f (S) ≤ best cost then To enhance this basic framework, our ILS is hybridized with a Variable Search Descending algorithm (VND) as local search component and with a set partitioning model as postoptimization phase. VND is a variant of the variable neighborhood search metaheuristic systematically changing the neighborhood each time no improvement is found by the current local search [START_REF] Mladenović | Variable neighborhood search[END_REF]. The set partitioning component uses an MILP to select a subset of routes (columns in the set partitioning model) among a set of feasible routes, in such a way that each customer is visited once. In the proposed algorithm, these feasible routes are generated and stored during the ILS phase. Adding a set partitioning component to an ILS structure has been used to solve other VRP variants demonstrating good performance: [START_REF] Subramanian | A hybrid algorithm for the heterogeneous fleet vehicle routing problem[END_REF] for the Heterogeneous Fleet Vehicle Routing Problem, [START_REF] Villegas | A matheuristic for the truck and trailer routing problem[END_REF] for the truck and trailer routing problem and [START_REF] Morais | Iterated local search heuristics for the vehicle routing problem with cross-docking[END_REF] for the the Vehicle Routing Problem with Cross-Docking.

Algorithm 1 presents the pseudo-code of the Hybrid ILS proposed. First, an initial solution S is generated using a randomized parallel insertion heuristic, as described in Subsection 4.3. In this step, time window and battery limit violations are allowed. To manage infeasible solutions, the objective function is penalized (see Subsection 4.1). Second, the perturbation component consists of removing a set of η customers and reinserting them at new positions. It is explained in detail in Subsection 4.4. Third, the local search is carried out by the VND component detailed in Subsection 4.5. In short, the VND is composed of classical operators adapted to manage RSs and satellite customers (see Subsection 4.2). Once the perturbation and the VND procedures are applied, a solution S is obtained. S is composed of a set of routes which are added to Ω, the set of "best" routes. The ILS component (perturbation + local search) is executed until M axILSIter iterations or N bN oImp iterations without improvement are reached. Finally, the set partitioning model is executed using Ω as the set of columns. This is explained in Subsection 4.6.

Generalized cost function

Since the problem is highly constrained (time windows, fleet size, battery and capacity), allowing the solution search to visit infeasible regions may be promising to facilitate the exploration, and this option is exploited in our approach.

To handle infeasible routes, the objective function of a solution S considers the total time spent at RSs, time window penalty, and battery penalty. This strategy is successfully implemented by [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF] for the location-routing problem with intraroute facilities. They test their method with various problems including the electric locationrouting problem with time windows and partial recharges. Their generalized cost function and penalty functions are taken and adapted, to be able to handle the satellite customers. For the Hybrid ILS, the generalized cost function is defined as:

f gen (S) = f (S) + λ(T W (S) + BT (S)) (27)
Where f (S) denotes the original objective function (Eq. 1). T W (S) represents time window violations and BT (S) battery violations. Those are scaled by the penalty factor λ. The magnitude of λ allows the local search to stay for shorter or longer periods in the unfeasible region. Thus, a large λ guides the search directly to feasible regions, while with a small λ to reach a feasible solution could take longer. This penalty approach is proven to be efficient (see, e.g., Goeke and Schneider (2015) for the E-VRPTW, [START_REF] Vidal | A hybrid genetic algorithm with adaptive diversity management for a large class of vehicle routing problems with time-windows[END_REF] for a large class of vehicle routing problems with time windows, and [START_REF] Nagata | A penalty-based edge assembly memetic algorithm for the vehicle routing problem with time windows[END_REF] for the vehicle routing problem with time windows).

Since a solution S is defined by a set of P routes (S = {π 1 , π 2 , ..., π P }), the total time spent at RSs f (S), the total time window penalty T W (S) and total battery penalty BT (S) are calculated by summing total time at RSs per route or route penalties for the routes of solution S, as follows:

f (S) = P j=1 f (π j ) (28) T W (S) = P j=1
T W (π j ) ( 29)

BT (S) = P j=1
BT (π j ) (30)

Move evaluation

Evaluating moves implies: computing the change in the objective function; verifying vehicle capacity feasibility, and computing time window and battery penalties. It is well known that the vehicle capacity can be checked in O(1). As explained by Goeke and Schneider (2015), feasibility in battery capacity can be also verified in O(1) if energy consumption is independent of cargo load. Likewise, [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF] shows that time window and battery penalties could be computed in O(1) for the EVRP-TWPR when a node is inserted between two partial routes. But if route segments with two or more nodes are inserted between two partial routes, penalties have to be extended for the route segment first.

In the proposed variant, walking time affects the time spent at an RS. It means that route's minimum duration depends on travel time, recharging time, walking times and time windows. The complexity of handling the time variation affects the evaluation of the objective function and the penalties. As mentioned above, in the proposed method time window and battery penalties evaluations are adapted from those proposed by [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF]. To evaluate time window and battery penalties they propose a corridor-based penalty approach. The time traveling concept (Time Warp) is used and it allows violations to be computed without overpenalizing route segments. This methodology permits inter-route moves to be evaluated in constant time, while for intra-route moves it requires updating information on the route segment affected by the movement. The latter can be achieved by having an efficient move evaluation. Therefore, forward functions, backward functions, and concatenation operators are needed and they are presented below.

Forward functions

For evaluating the proposed functions all resource dependencies are modeled in time units. The notation used by [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF] is preserved and extended to the new variables. Thus, time needed to recharge the energy consumed on an arc (i, j) is denoted as In E-VRP problems with partial recharges, the RSs usually do not have short time windows (even any time window) or fixed service time. However, in the proposed model if an RS has assigned a customer i ∈ V to be visited by walking, it is necessary to consider customer time window and service time. For this reason, for each RS visit j ∈ F a service time (walking time) s j , an earliest arrival time e j , and a latest arrival time l j are computed before each local search iteration. This is done using the following equations ( 31) -( 33).

h ij = r * d ij
s j = αt ji + s i + αt ij if a customer i ∈ V is visited "by walking" 0 otherwise (31)
e j = max{e i -αt ji , e j } if a customer i ∈ V is visited "by walking" e j otherwise (32)

l j = min{l i -αt ji , l i } if a customer i ∈ V is visited "by walking" l i otherwise (33)
In the following, the equations to extend each variable a x with x ∈ {rw, cost, min, max, sl, rt, add} from vertex i to vertex j are presented. A brief description of how variables move forward on a route from vertex i to vertex j is given, and the main differences from the original functions is explained. However, for a detailed description of each variable and the different penalty cases we refer the reader to [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF].

To compute the objective function, two new variables are added. First, a rw i keeps track of the cumulated extra recharging time spent at previous RSs due to walking at leaving vertex i. So, if a rw i > 0, it means that the EV spent more time at the RS because of walking activity than the time required to recharge the energy consumed up to node i. It is obtained by adding the walking time s j when a RS j is visited, and subtracting the energy consumed on the arc (i, j).

a rw j = min{s j + max{a rw i -h ij , 0}, H} if j ∈ F max{a rw i -h ij , 0} otherwise (34)
Second, a cost i indicates the cumulated time spent at RSs when leaving vertex i. It increases when visiting a customer by walking or because of the energy consumed on the arc (i, j) when a rw i = 0 (when no extra energy is recharged in previous RSs).

a cost j = a cost i + max{h ij -a rw i , 0} + s j if j ∈ F a cost i + max{h ij -a rw i , 0} otherwise (35)
The earliest arrival time a min j for vertex j and the slack a sl ij are adapted. In the proposed version, those variables include e i and s i for the cases when an extension is made from a RS where a customer is visited by walking.

a min j = max{e j , max{e i , ãmin i } + s i + t ij } + a add ij if i ∈ F max{e j , ãmin i + s i + t ij } + a add ij otherwise ( 36
)
a sl ij = max{e j -(max{e i , ãmin i } + s i + t ij ), 0} if i ∈ F max{e j -(ã min i + s i + t ij ), 0} otherwise (37)
In the same way, variables a max j , a rt j , and a add ij are modified by including e i and s i for the cases when an extension is made from a RS where a satellite customer is visited. However, because of the changes on variable a min j , a new slack must be considered. It is defined by the expression max{e i -ãmin i , 0}. Following the strategy observed with the walking time s i and a sl ij , this new slack time is included and considered as a recharging time.

a max j =      max{e j , max{e i , ãmin i } + max{a rt i -s i -max{e i -ãmin i , 0}, 0} + s i + t ij } if i ∈ F max{e j , ãmax i + s i + t ij , 0} otherwise (38)
a rt j = min{H, max{a rt i -s i -a sl ij -max{e i -ãmin i , 0}, 0} + h ij } if i ∈ F min{H, max{a rt i -min{a sl ij , ãmax i -ãmin i }, 0} + h ij } otherwise ( 39
)
a add ij = max{max{a rt i -s i -a sl ij -max{e i -ãmin i , 0}, 0} + h ij -H, 0} if i ∈ F max{max{a rt i -min{a sl ij , ãmax i -ãmin i }, 0} + h ij -H, 0} otherwise (40)
Time window penalties and battery penalties can be calculated following the original functions. Likewise, the shifts between a min i and a max i to ãmin i and ãmin i remain the same. The only thing to keep in mind is when i ∈ F or v ∈ F (when a RS is visited) l i and l v must be replaced with l i and l v respectively in equations ( 41)-( 44).

--→ T W (π i ) = v∈π i max{min{a max v , a min v } -l v , 0} (41)
-→ BT (π i ) = v∈π i max{a min v -a max v , 0} (42)
ãmin i = min{a min i , a max i , l i } (43) ãmax i = min{l i , min{a min i , a max i , l i } + max{a max i -a min i , 0}} (44)
Initial values for these variables to estimate forward functions on a route are a cost 0 =a rw 0 = a min 0 =a max 0 =ã min 0 =ã max 0 =a rt j =0. With these adjusted functions, the corridor based-approach presented by [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF] is able to address the forward functions including partial and time-dependent charges, and satellite customers. The backward functions are adapted as well, following the same strategy used for the forward functions. The original notation is preserved, so backward functions are referred to as b x with x ∈ {cost, rw, min, max, sl, rt, add} and as bmax i and bmin i as well. Finally, these functions are listed in Appendix A.

Concatenation operators

The formulae for time-efficient route concatenation are also adjusted and presented below. A new formula to estimate the objective function is presented as well.

For a route π e = 0, ..., x, v, y, ..., n + 1 , which is constructed by inserting vertex v between two partial routes π 1 = 0, ...x and π 2 = y, ..., n + 1 , the objective function of route π e is computed using Equation ( 45). In that case it is necessary to compute the forward and backward functions for node v.

f (π e ) = a cost v + b cost v -min{a rw v , b rw v } (45)
The time window penalty for route π e is calculated in constant time using Equation ( 46). It sums the penalty of each partial path. Additional penalties are added if the time window [e v , l v ] or the overall route duration is violated. Note that in the case where v ∈ F , e v and l v must be replaced with e v and l v .

T W (π e ) = --→ T W (π 1 ) + ← -- T W (π 2 ) + max{0, a min v -l v -max{0, a min v -a max v }}+ max{0, max{e v , min{a min v , a max v , l v }} -b min v -max{0, b max v -b min v }} (46)
The battery penalty for route π e is calculated in constant time using Equation (47). It includes the penalty of each partial path. Additional penalties are added if the new route contains a route segment where the energy required exceeds the battery capacity. If an RS is located at vertex v, additional recharging time or additional walking time is included. However, additional time window penalties due to time at RS are already included in Equation ( 46) and ignored within ( 47) and ( 48

BT (π e ) = -→ BT (π 1 ) + ←- BT (π 2 ) + max{0, a min v -a max v } + max{0, b max v -b min v } + C (47) C =                          max{0, max{0, a rt v -max{0, e v -(a min v -max{0, a min v -a max v })}} + max{0, b rt v -s v -max{0, b min v + max{0, b max v -b min v } -l v }} -H- min{max{0, a rt v -max{0, e v -a min v + max{0, a min v -a max v }}}, max{0, min{l v , max{b max v , b min v }} -max{e v , min{a min v , a max v }}}}} if v ∈ F max{0, a rt v + b rt v -H -min{H, max{0, b min v -b max v } + max{0, a max v -a min v }, max{0, min{l v , max{b min v , b max v }} -min{a min v , a max v }}}} otherwise ( 48
)
For a route π e = 0, ..., x, y, ..., n + 1 , constructed by concatenating two partial routes π 1 = 0, ...x and π 2 = y, ..., n + 1 , the objective function of route π e is computed using Equation ( 45). In that case it is necessary to compute the forward functions for node y. Following the strategy applied above, the time window penalty can be determined by ( 49), while the battery penalty is determined with ( 50) and (51).

T W (π e ) = --→ T W (π 1 ) + ← -- T W (π 2 ) + max{0, a min y -l y -max{0, a min y -a max y }}+ max{0, max{e y , min{a min y , a max y , l y }} -bmin y } (49) BT (π e ) = -→ BT (π 1 ) + ←- BT (π 2 ) + max{0, a min y -a max y } + D (50) D =                          max{0, max{0, a rt y -max{0, e y -(a min y -max{0, a min y -a max y })}} + max{0, b rt y -s v -(b min v -l v )} -H- min{max{0, a rt v -max{0, e v -a min v + max{0, a min v -a max v }}}, max{0, min{l v , b min v } -max{e v , min{a min v , a max v }}}}} if v ∈ F max{0, a rt v + b rt v -H -min{H, max{0, b min v -b max v } + max{0, a max v -a min v }, max{0, min{l v , max{b min v , b max v }} -min{a min v , a max v }}}} otherwise ( 51
)
With these expressions, the objective function, the time window penalties, and the energy violations can be derived in constant time for every inter-route neighborhood structure considered in this document. It must be remembered that for the intra-route movements the evaluation is done in linear time.

Initial solution

The Hybrid ILS starts by generating an initial solution using a greedy randomized parallel insertion heuristic. First, P routes are opened. Second, the set of customers is sorted in decreasing order based on their latest arrival time l i . Then, a random customer from the first φ customers is selected. Using the generalized cost function as criteria (see Subsection 4.1), the node is inserted following a best insertion strategy. During the procedure the vehicle capacity constraint is always respected. This procedure is repeated until all customers are assigned to a route. Afterwards, some battery and time window infeasible routes can remain. Finally, the local search is applied (see Subsection 4.5) to repair the initial solution. Note that, at the end, the initial solution may still contain some time window and battery limit violations which will be managed through the penalties during the ILS phase.

Perturbation

The perturbation consists of removing η random customers, then grouping them into a set of removed nodes. Later, one node is selected randomly and it is inserted in a arbitrary position. The new position must be different from the original location, and the capacity constraint must be guaranteed. It is done one by one for all the removed customers. If a customer cannot be inserted in any route because of capacity constraint, the process is restarted.

Variable Neighbourhood Descent -VND

The VND is implemented as local search strategy following a best improvement strategy. It is composed of a neighborhood set based on classic operators (Relocate, Swap and 2-OPT*). Additionally, operators Insert RS and Remove RS for inserting and removing visits to recharging stations are included. These operators are used successfully by [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF] and Goeke and Schneider (2015) showing high performance. The operators included are described below.

• Relocate: Removes one node from its current route and reinserts it elsewhere.

• Swap: Exchanges two non consecutive nodes.

• 2-OPT*: Removes two edges from two different routes and reconnects them by two new edges.

• InsertRS: Inserts a visit to a RS. There are two types of insertion. First a simple visit to a recharging station is evaluated. Second, an insertion where a customer visit is replaced by a RS visit, and the customer is assigned to be visited by walking.

• Remove RS: Removes a visit to a RS. If the RS has assigned a customer to be visited by walking, the evaluation considers that the customer is inserted in the position where the RS was.

Relocate and Swap are applied for intra-and inter-route moves, 2-OPT* is designated for inter-route moves. Swap is determined to handle only customers. Relocate handle customers and for intra-route moves it deals with RSs. And as it is described above, 2-OPT* is only applied for inter-route moves.

Additionally, new specific operators for assigning or removing customers serviced by walking are included. These operators still rely on the Relocate and Swap strategies.

• BecomeSatellite: Assign a customer to be visited by walking from an already visited RS.

• BecomeByVeh: Insert a customer that was assigned to be visited by walking into a the vehicle route.

• SwapVehSat: Exchange between a customer that is visited by walking and a customer that is visited in the vehicle route.

The above-mentioned movements are applied for intra-and inter-route changes. In the VND, first the inter-routes moves are evaluated, because they are where evaluations are performed in constant time. Those are followed by InsertRS and the intra-route versions of the moves. At the end RemoveRS looks for ways to clean unneeded RS stops from the solution.

Set partitioning component

The set partitioning model is solved once ILS finishes. Let Ω be the set of routes stored during the ILS phase. Let ρ k be the cost of the route r k ∈ Ω and ξ ik = 1 if route r k visits node i and 0 otherwise. Thus, the E-VRPTWsc can be described with the following set partitioning model: min

r k ∈Ω ρ k θ k (52) r k ∈Ω θ k ≤ P (53) r k ∈Ω ξ ik θ k = 1 ∀i ∈ V (54) θ k ∈ {0, 1} ∀r k ∈ Ω (55)
The binary decision variable θ k is equal to 1 if route r k is selected in the solution, and equal to 0 if not. Constraint (53) limits the number of vehicles that are used. Finally, constraints (54) guarantee that customers must be served once.

Computational experiments

A set of instances for E-VRPTWsc can be easily obtained from the E-VRPTW instances as the E-VRPTW is a special case of the E-VRPTWsc where no customer is assigned to be visited by walking and full charges are performed.

The set of instances proposed by [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF] for the E-VRPTW are created based on the benchmark instances proposed by [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF] for the VRPTW. The instances are classified into three groups: clustered customer distribution (c), random customer distribution (r) and a mixture between random/clustered customer distribution (rc).

Since the proposed variant includes an alternative transportation mode, it is meant to be suitable in a city logistics framework. Thus, the c instances are interesting to analyze because they represent geographical scenarios where distances, times and speed allow us to take advantage of the recharging times. A brief analysis is provided in Section 5.4 and Section 5.5.

For large-size instances, three sets of instances are created by fixing α parameter to 1.5, 2 and 4. It represents scenarios with three different types of walking speeds, in other words, three types of alternative modes. For small size instances parameter α is fixed at 2. NB α = 2 means that the time of traveling through an arc by walking takes twice the time that is required by the EV to perform the same distance. The number of vehicles is established using the values of the updated BKS for the E-VRPTWPR. To the best of our knowledge, these BKS are found in the results presented by Keskin and C ¸atay ( 2016), Schiffer and[START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF]Hiermann et al. (2019).

Section 5.1 presents the parameter tuning. The influence of the initial solution quality is analyzed in Section 5.2. Results for the large size E-VRPTWPR (no walking allowed) are shown in Section 5.3 to prove the competitiveness of the Hybrid-ILS. In Section 5.4 and 5.5 results for small and large size E-VRPTWsc instances are presented. A comparison with the BKS of the E-VRPTWPR is also presented, to evaluate the pertinence of the variant.

The hybrid-ILS is implemented in C++. The method is executed on a machine with Intel Core i5-5300U processor with 2.30GHz speed and 8 GB RAM.

Parameter setting

The hybrid-ILS proposed has four parameters to be tuned: φ, η, λ ,and M axILSIter (see Section 4). The influence of modifying φ, η and λ parameters is evaluated by running a subset of 14 instances five times with different parameter combinations. The size of the subset of candidate customers to be inserted φ is used to build the initial solution, and it is expressed as a percentage of the total number of customers. η represents the number of customers to be relocated in the perturbation, and it is expressed as a percentage of the total number of customers. Penalties are scaled by λ. Since penalties are expressed in time units, λ is fixed as a percentage of the T max (time horizon). The different combinations of η and λ parameters tested are presented in Table 3. Also presented are average objective function GAP percentage of the best of 5 (∆f best ), and average objective function GAP percentage of the average of 5 (∆f ). In the test M axILSIter is fixed at 5000 iterations and φ in 0.05|V |. Average computational time t is presented in minutes. Average number of iterations required to reach the best solution of each run is presented as κ and average time to reach the best solution of each run is noted by t best . GAPs are measured with respect to the best solution found for each instance during the entire test.

The results show that smaller GAPs can be obtained by setting λ between 1% and 5%, and fixing η between 15% and 20%. Because the quality of solutions is not very different, and based on the execution time parameter, the configuration λ = 0.01T max and η = max{2, 0.15|V |} is selected. η parameter is presented in this way because in small instances with 5 customers, 15% of total customers represents just one customer and the perturbation is not large enough. Once these parameters are fixed, the influence of the number of iterations M axILSIter is analyzed by plotting the behavior of the ILS component on 5 runs of the algorithm for four different instances.

Figure 1 shows the results of the test. Based on that information, the parameter M axILSIter is set at 4000 iterations and N bN oImp at 2500 iterations without improvement.

Initial solution impact

The influence of the initial solution is analyzed in two ways. First, the influence of the randomness on the construction of the initial solution is tested versus a purely greedy construction. Table 4 shows the different φ values tested, and ∆f best and ∆f values introduced above. T W min show the average of the smallest penalty value found in the initial solution of the 5 runs. T W min the average penalty value in the initial solution. Second, by looking at Figure 1 it is possible to see how the algorithm performs by using initial solutions with different qualities. The graphs do not show a clear link between the quality of the initial solution and the quality of the final solution.

According to the results, φ parameter is fixed on 0.05|V |. This is because, on average, configurations with small φ values produce initial solutions with fewer penalties and to reach better solutions. Likewise, results show that having a randomness in the initial solution allow us to reach better solutions. Table 4: Results using different parameter settings for φ on a subset of the E-VRPTWsc instances.

Results for the EVRP-TWPR

To validate the performance of the proposed Hybrid ILS approach, this method is used to solve EVRP-TWPR instances. Thus, in the VND component every evaluation involving nodes visited by walking are dropped.

The best solutions, obtained after executing the Hybrid ILS 10 times for each instance, are reported in Table 5. Those results are compared with the best known solutions (BKS) found by state-of-the-art algorithms: the ALNS proposed in Keskin and C ¸atay (2016) (KC ¸), the ALNS presented in Schiffer and Walther (2018) (S&W) and the hybrid genetic algorithm presented in Hiermann et al. (2019) (HGA). As none of the methods finds the BKS for all the instances, an individual comparison with each method is presented as well. The results shown in the table for other authors also correspond to the best solution found after 10 runs.

The results in Table 5 show that our algorithm is able to reach solutions for the E-VRPTWPR with an average gap of -0.34% with regard to BKS. Moreover, it is able to BKS denotes the best known solution. For each solution method, the number of vehicles m, the total distance ϕ D and the percentage gap ∆ D of the best solution obtained in 10 runs of our algorithm are presented. For our solution method average run-times are presented in column t (in minutes).

BKS KC ¸1 S&W 2 HGA 3 Hybrid ILS 4 Inst. m ϕ a D m ϕ b D m ϕ c D m ϕ d D m ϕ d D ∆ a D ∆ b D ∆ c D ∆ d D t c101
1 Intel Xeon E5 processor with 3.30 GHz speed and 32 GB RAM.

2 Intel Core i7 processor with 3.60 GHz speed and 16 GB RAM.

3 Intel Xeon 2643 processor with 3.30 GHz speed and 4 GB RAM.

4 Intel Core i5-5300U processor with 2.30 GHz speed and 8 GB RAM.

-Solution uses different number of vehicles.

Table 5: Results for the large-size EVRP-TWPR instances allowing satellite customers can reduce the total charging time by 9.66% and the vehicle distance by 10.66%, on average, compared to the BKS for the E-VRPTWPR.

Comparing the reductions obtained between the c sets of 5, 10 and 15 customers, a slightly decrease of the distance and recharging time improvements can be observed. It is obvious that the potential reduction in distance and time is correlated with the number of visits to RSs and number of customers. The larger the ratio between number of visits to RSs and number of customers, the bigger the reduction in vehicle distance and total energy consumed can be achieved, with the proposed variant.

A second numerical test is done with the small size instances. The performance of the hybrid ILS and the commercial software CPLEX on small instances of the E-VRPTWsc is compared. The mathematical model presented in Section 3 is used in this test. In the model the set F determines the number of RS copies and, implicitly, the number of times that a RS can be visited. On one hand, in accordance with the problem definition, one RS can be visited several times. It can be represented by a large number of RS copies. On the other hand, the larger the F set is the more complex the graph is. For solving the instances with CPLEX an iterative procedure is followed. First the model is solved with no RS copies. If the problem is solved to optimality, a new iteration adding one copy of each RS is done. The process is repeated until the optimal value of two consecutive iterations is equal, or until the current iteration cannot be solved to optimality with a time limit of 7200 seconds. The comparison is presented in Table 7. In column ϕ R proven optimal solutions are underlined and best solutions are in bold. Bold values in columns ∆ D and ∆ R show that Hybrid ILS finds a feasible solution better than feasible solution provided by CPLEX (for non optimal solutions).

RS, in less than 2 seconds.

Results on the large E-VRPTWsc instances

The proposed Hybrid ILS is used to solve large-size instances. The results of the computational test are shown in Table 8. The results are presented in three sets according to the α value. Results demonstrate that allowing visits by walking could reduce the total time at RSs and the total distance of the EVs in the three proposed scenarios compared with the BKS for the E-VRPTWPR (note that the BKS are updated according to the results obtained in Section 5.3). Doing an analysis by subset of instances, the faster the alternative mode is (smaller values of α), the larger is the reduction obtained in total distance and total time at RS. Strictly speaking, when α is 1.5 it is possible to reduce by 3.53% total time spent at RS and by 3.75% the total distance of EVs. For α equal to 2 total time spent at RS is reduced 2.33% and the total distance 2.56%. And a reduction of 0.79% on the total time spent at RS and of 0.82% on the total distance is achieved when α is 4. According to this, it is clear that a vehicle that moves faster allows you to get better results, as would be expected.

In Appendix B one table for each set of instances is presented. In each table, number of visits to RSs, number of satellite customers, recharging times and walking times are presented for the solutions presented in Table 8. Likewise, column ϕ R -ϕ CT shows the amount of extra time spent at the RS due to walking. If we look at instance c107 in Table B.9, we can find the largest value of additional time due to visiting a satellite customer. Comparing the 77.03 units of time with the total recharge time, this represents 2.3% of the recharging time. When the difference is equal to zero, it means that walking time is a hidden time that does not affect the duration of the route. Therefore, it is verified that the variant is taking advantage of the dead times caused by the recharging activity.

Looking for the type of geographical distribution, the c instances are where greater reductions in time at RSs and distance can be obtained. It follows the behavior observed in the small instances and shows the benefit of allowing satellite customers when using EVs in clustered scenarios or in a city logistics context.

On average, the proposed method obtains high quality solutions in 3.30 minutes. The results show that when α is 1.5 the algorithm takes longer than when α is 4. This happens because with better walking speeds the number of potential customers that can be visited by walking increases. Finally, the Hybrid-ILS takes longer to solve E-VRPTWsc instances compared to the results reported for the EVRP-TWPR. This is caused by the walking evaluations.

Conclusions and perspectives

In this paper the Electric Vehicle Routing Problem with Time Windows and Satellite Customers is introduced as a novel variation of the E-VRPTW. The objective function is the minimization of the time spent at recharging stations. This also allows a reduction in the total vehicle distance. A mathematical model is presented and a solution method based on an ILS structure is designed. A generalized cost function, forward functions, backward functions and concatenation operators are proposed to evaluate most of the local search movements in constant time.

Computational experiments show that, allowing the visits to the customer by an alternative mode enable us to: 1) reduce total recharging time, which means a reduction in total vehicle distance; and 2) avoid part of the idle times during the intra-route recharging activities. The efficiency of the algorithm is recognized by the tests on a special case, E-VRPTWPR where new BKS solutions are found and the average computational time is less than the time reported for State-of-the-art algorithms.

  procedure ILS is an effective method for several vehicle routing problems. As Lourenço et al. (2003) detail, it generates a sequence of local optima by alternating local search and perturbation.

  , and time needed to recharge the full battery capacity is represented by H = r * B. The variables required to estimate the objective function, the time window penalty and the battery penalty are presented in

  ) to avoid overpenalization. Because of the changes on the variables, the expressions max{0, e vl v }, which represent new slack times (for v ∈ F ), are included as recharging times.

Figure 1 :

 1 Figure 1: (Color online) Behavior of ILS component on four E-VRPTWsc instances

Table 1 :

 1 Variable and parameter definitions of the E-VRPTWsc model.

	Sets
	V	Set of customer vertices
	F	Set of dummy vertices that represents the visits to RS on F
	V	Set of nodes, recharging visits and depot nodes V = {V ∪ F ∪ {0, N + 1}}
	Parameters
	0, N + 1 Depot nodes
	P	Number of electric vehicles
	Q	Load capacity
	B	Battery capacity
	q i	Demand of customer i
	s i	Service time of customer i
	e i	Earliest start of service at node i
	l i	Latest start of service at node i
	d ij	Distance between vertices i and j
	t ij	Travel time between vertices i and j
	γ	Consumption rate
	r	Inverse recharging rate
	α	Walking speed ratio
	Decision Variables
	t i	Recharging time spend at RS i
	t ij	Walking time at RS i caused by visiting customer j by walking
	IC i	Initial SoC that is required by the vehicle that departs from the depot and arrives at node i
	T RT i	Time spent at the RS i during the intraroute charging
	τ i	Arrival time at node i
	u i	Remaining cargo on arrival at node i
	w i	Amount of energy recharged at RS i
	x ij	Boolean variable indicating if arc (i, j) is traversed
	y i	SoC on arrival at node i
	z ij	Boolean variable indicating if the customer j is visited by walking from RS i

Table 2 .

 2

	a cost i a rw i a min i a max i a sl ij a rt i a add ij ãmin i ∧ ãmax i	Cumulated recharging time at leaving vertex i Cumulated extra recharging time due to walking at leaving vertex i Earliest allowed arrival time at a vertex i Arrival time at vertex i if as much energy as possible was recharged at preceding RS Slack between vertex i and vertex j due to time windows when traveling from i to j. Inverse residual battery capacity at vertex i Additional energy that has to be charged at the preceding RS to travel arc (i, j). If a violation occurs a min i and a max i are shifted to ãmin i and ãmax i . It prevents repeated
		penalization.

Table 2 :

 2 Forward functions variables

Table 3 :

 3 Results using different parameter settings for η and λ on a subset of the E-VRPTWsc instances.

  number of vehicles, total distance performed by the vehicles ϕ D , total time spent at RS ϕ R , gap percentage ∆ D for total distance, gap percentage ∆ R for total recharging time. Average run-times are presented in column t (in minutes).

	EVRP-TWPR BKS			EVRP-TWsc α = 1.5					EVRP-TWsc α = 2					EVRP-TWsc α = 4
	Inst. m	ϕD	ϕR	m	ϕD	∆D	ϕR	∆R	t	m	ϕD	∆D	ϕR	∆R	t	m	ϕD	∆D	ϕR	∆R	t
	c101 12 1043.38 3537.06	12 992.71	-4.86 3365.30 -4.86 2.61	12 994.53	-4.68 3371.49 -4.68 2.00	12 1030.64 -1.22 3493.88 -1.22	1.72
	c102 11 1017.70 3450.00	11 913.03 -10.29 3103.82 -10.03 2.38	11 934.65	-8.16 3174.41 -7.99 1.85	11 978.12 -3.89 3344.67 -3.05	1.90
	c103 10 971.19 3292.33	10 847.33 -12.75 2880.61 -12.51 3.13	10 851.85 -12.29 2936.13 -10.82 2.90	10 933.29 -3.90 3186.07 -3.23	3.10
	c104 10 884.38 2998.04	10 802.48	-9.26 2722.25 -9.20 3.75	10 820.17	-7.26 2780.38 -7.26 3.39	10 856.49 -3.15 2912.83 -2.84	3.13
	c105 11 1015.79 3443.52	11 922.40	-9.19 3175.66 -7.78 3.12	11 978.55	-3.67 3317.31 -3.67 2.81	11 1007.48 -0.82 3415.35 -0.82	2.02
	c106 11 1009.33 3421.61	11 904.69 -10.37 3132.49 -8.45 3.07	11 936.53	-7.21 3237.32 -5.39 3.20	11 990.76 -1.84 3358.68 -1.84	2.12
	c107 10 1046.50 3547.63	10 895.37 -14.44 3112.34 -12.27 3.13	10 931.43 -11.00 3228.79 -8.99 3.13	10 1008.47 -3.63 3418.71 -3.63	2.64
	c108 10 1022.48 3466.20	10 875.34 -14.39 2998.60 -13.49 3.27	10 884.49 -13.50 3048.56 -12.05 3.01	10 966.62 -5.46 3276.84 -5.46	2.73
	c109 10 940.38 3187.90	10 821.30 -12.66 2819.95 -11.54 3.99	10 837.99 -10.89 2912.25 -8.65 3.70	10 906.86 -3.57 3074.24 -3.57	3.45
	c201 4 629.95 1436.28	4	608.51	-3.40 1387.39 -3.40 2.39	4	610.38	-3.11 1392.90 -3.02 1.67	4	619.44 -1.67 1412.32 -1.67	1.52
	c202 4 629.95 1436.28	4	606.66	-3.70 1383.17 -3.70 2.68	4	610.63	-3.07 1392.23 -3.07 1.98	4	619.00 -1.74 1411.32 -1.74	1.89
	c203 4 629.95 1436.28	4	601.88	-4.46 1372.28 -4.46 3.53	4	606.59	-3.71 1383.02 -3.71 3.31	4	615.47 -2.30 1403.27 -2.30	2.44
	c204 4 628.91 1440.21	3	581.98	-7.46 1332.74 -7.46 3.47	4	599.74	-4.64 1373.40 -4.64 3.94	4	614.16 -2.35 1406.43 -2.35	4.21
	c205 4 629.95 1442.58	4	601.60	-4.50 1377.65 -4.50 2.61	4	603.90	-4.13 1382.93 -4.13 2.55	4	616.17 -2.19 1411.03 -2.19	2.12
	c206 4 629.95 1442.58	4	597.25	-5.19 1367.70 -5.19 3.09	4	599.62	-4.81 1373.13 -4.81 3.06	4	611.54 -2.92 1403.91 -2.68	4.03
	c207 4 629.95 1442.58	4	596.95	-5.24 1367.01 -5.24 3.63	4	598.90	-4.93 1371.48 -4.93 3.33	4	614.18 -2.50 1406.80 -2.48	4.29
	c208 4 629.95 1442.58	4	589.38	-6.44 1349.69 -6.44 3.15	4	593.58	-5.77 1361.21 -5.64 3.40	4	611.93 -2.86 1401.36 -2.86	3.03
	r101 18 1606.98 771.35	18 1579.05 -1.74	761.75	-1.24 3.41	18 1600.74 -0.39	768.36	-0.39 3.17	18 1605.27 -0.11 770.53 -0.11	2.86
	r102 15 1461.23 701.39	15 1400.79 -4.14	673.85	-3.93 3.65	15 1472.88	0.80	709.28	1.13	3.39	15 1457.14 -0.28 699.63 -0.25	3.04
	r103 13 1212.37 581.94	13 1166.23 -3.81	563.98	-3.09 3.81	13 1197.93 -1.19	575.01	-1.19 3.56	13 1210.53 -0.15 581.06 -0.15	3.29
	r104 11 1051.41 473.14	11 1014.91 -3.47	459.17	-2.95 3.91	11 1028.61 -2.17	463.80	-1.97 3.73	11 1050.22 -0.11 472.60 -0.11	3.75
	r105 14 1347.80 646.94	14 1304.98 -3.18	626.39	-3.18 3.66	14 1311.36 -2.70	631.14	-2.44 3.56	14 1341.09 -0.50 643.72 -0.50	3.25
	r106 13 1256.19 602.97	13 1204.79 -4.09	579.67	-3.87 3.83	13 1224.60 -2.51	589.10	-2.30 3.61	13 1254.83 -0.11 602.32 -0.11	3.43
	r107 12 1108.47 498.81	11 1066.91 -3.75	480.31	-3.71 3.37	11 1088.56 -1.80	493.82	-1.00 3.25	11 1101.31 -0.65 495.59 -0.65	3.04
	r108 11 1020.52 479.64	10 968.39	-5.11	457.57	-4.60 3.65	10 995.50	-2.45	470.33	-1.94 3.69	11 1019.03 -0.15 478.95 -0.15	3.62
	r109 12 1185.77 545.46	12 1131.62 -4.57	523.76	-3.98 4.08	12 1147.24 -3.25	530.34	-2.77 3.68	12 1171.93 -1.17 539.09 -1.17	3.54
	r110 11 1070.99 481.95	11 1031.28 -3.71	465.82	-3.35 4.28	11 1040.83 -2.82	469.68	-2.55 4.09	11 1067.26 -0.35 480.27 -0.35	3.95
	r111 11 1072.46 493.33	11 1034.68 -3.52	475.95	-3.52 3.58	11 1046.19 -2.45	481.89	-2.32 4.13	11 1065.85 -0.62 490.29 -0.62	3.96
	r112 11 1001.79 460.82	10 969.09	-3.26	446.15	-3.18 4.46	10 978.38	-2.34	450.06	-2.34 3.73	11 998.42 -0.34 459.27 -0.34	3.19
	r201 3 1255.81 200.93	3 1249.36 -0.51	199.90	-0.51 3.77	3 1250.29 -0.44	200.05	-0.44 3.54	3 1251.29 -0.36 200.21 -0.36	3.12
	r202 3 1051.46 136.69	3 1040.71 -1.02	135.29	-1.02 2.94	3 1043.66 -0.74	135.68	-0.74 2.84	3 1047.45 -0.38 136.17 -0.38	2.52
	r203 3 895.54 143.29	3	885.79	-1.09	141.73	-1.09 3.56	3	884.99	-1.18	141.60	-1.18 3.49	3	893.19 -0.26 142.91 -0.26	2.83
	r204 2 780.91	93.71	2	767.87	-1.67	92.18	-1.63 5.06	2	768.80	-1.55	92.26	-1.55 4.94	2	781.22	0.04	93.75	0.04	4.98
	r205 3 987.22 148.08	3	973.71	-1.37	146.06	-1.37 3.47	3	978.27	-0.91	146.74	-0.91 2.61	3	986.39 -0.08 147.96 -0.08	2.71
	r206 3 922.70 156.86	3	913.82	-0.96	155.35	-0.96 3.59	3	930.09	0.80	158.11	0.80	3.29	3	929.47	0.73	158.01	0.73	3.14
	r207 2 843.20	92.75	2	839.48	-0.44	92.34	-0.44 4.50	2	843.00	-0.02	92.73	-0.02 4.57	2	844.49	0.15	92.89	0.15	3.97
	r208 2 736.12 103.06	2	726.83	-1.26	101.76	-1.26 4.55	2	728.57	-1.03	102.00	-1.03 4.14	2	733.36 -0.37 102.67 -0.37	3.31
	r209 3 863.36 138.14	3	858.36	-0.58	137.39	-0.54 3.20	3	863.63	0.03	138.18	0.03	3.15	3	865.07	0.20	138.41	0.20	3.03
	r210 3 843.36 134.94	3	832.39	-1.30	133.18	-1.30 3.28	3	837.98	-0.64	134.08	-0.64 2.92	3	842.81 -0.07 134.85 -0.07	3.00
	r211 2 826.88	90.96	2	824.23	-0.32	90.67	-0.32 4.91	2	822.03	-0.59	90.42	-0.59 4.65	2	838.67	1.43	92.25	1.43	3.97
	rc101 15 1648.99 626.62	15 1646.23 -0.17	625.57	-0.17 3.34	15 1646.24 -0.17	625.69	-0.15 2.93	15 1658.89 0.60	630.38	0.60	2.64
	rc102 14 1510.16 573.86	14 1496.72 -0.89	568.75	-0.89 3.13	14 1507.28 -0.19	572.77	-0.19 2.67	14 1509.09 -0.07 573.46 -0.07	2.46
	rc103 12 1346.83 511.80	12 1321.30 -1.90	502.10	-1.90 3.24	12 1348.50	0.12	512.55	0.15	2.91	12 1345.21 -0.12 511.18 -0.12	2.74
	rc104 11 1175.06 446.52	11 1150.18 -2.12	439.90	-1.48 3.00	11 1173.11 -0.17	445.78	-0.17 3.51	11 1173.46 -0.14 445.92 -0.14	2.70
	rc105 14 1446.30 549.60	13 1424.16 -1.53	541.18	-1.53 3.02	14 1440.89 -0.37	547.54	-0.37 2.69	14 1446.30 0.00	549.60	0.00	2.66
	rc106 13 1383.14 526.66	13 1377.88 -0.38	523.60	-0.58 3.15	13 1382.21 -0.07	525.24	-0.27 2.86	13 1387.60 0.32	527.29	0.12	2.68
	rc107 12 1244.83 473.03	12 1237.82 -0.56	470.37	-0.56 3.09	12 1243.08 -0.14	472.37	-0.14 2.89	12 1240.71 -0.33 471.47 -0.33	2.83
	rc108 11 1154.14 438.57	11 1150.76 -0.29	437.29	-0.29 3.15	11 1154.11	0.00	438.56	0.00	2.58	11 1158.81 0.40	440.35	0.40	2.68
	rc201 4 1443.07 202.03	4 1431.83 -0.78	200.46	-0.78 2.92	4 1440.82 -0.16	201.72	-0.16 2.69	4 1441.75 -0.09 201.85 -0.09	2.38
	rc202 3 1403.32 154.37	3 1392.35 -0.78	153.16	-0.78 3.91	3 1411.69	0.60	155.29	0.60	3.80	3 1405.73 0.17	154.63	0.17	3.63
	rc203 3 1060.32 148.44	3 1062.11	0.17	148.70	0.17	3.61	3 1051.91 -0.79	147.27	-0.79 3.60	3 1067.60 0.69	149.46	0.69	3.74
	rc204 3 884.75 168.10	3	872.88	-1.34	165.85	-1.34 3.36	3	874.10	-1.20	166.08	-1.20 3.39	3	890.17	0.61	169.13	0.61	2.57
	rc205 3 1249.56 187.43	3 1222.27 -2.18	183.34	-2.18 4.82	3 1246.91 -0.21	187.04	-0.21 4.36	3 1244.23 -0.43 186.64 -0.43	4.49
	rc206 3 1187.40 154.36	3 1173.62 -1.16	152.57	-1.16 3.73	3 1199.77	1.04	155.97	1.04	3.75	3 1198.47 0.93	155.80	0.93	3.48
	rc207 3 985.67 137.99	3	971.08	-1.48	135.95	-1.48 3.74	3	987.01	0.14	138.87	0.63	3.85	3	988.50	0.29	138.39	0.29	3.56
	rc208 3 833.12 149.96	3	825.31	-0.94	148.56	-0.94 4.32	3	836.23	0.37	150.52	0.37	3.71	3	838.20	0.61	150.88	0.61	3.27
	Average					-3.75		-3.53 3.52			-2.56		-2.33 3.31			-0.82		-0.79	3.08

m

Table 8 :

 8 Results for large-size instances for the E-VRPTWsc

  number of vehicles, total distance performed by the vehicles ϕD, total time spent at RS ϕR, number of visits to RS nRS, number of satellite customers St., charging time ϕCT , and ϕW T total walking time are presented. Average run-times are presented in column t (in minutes). Finally, column ϕR -ϕCT shows the amount of extra time spent at the RS due to walking TableB.10: Detailed results for the large-size E-VRPTWsc instances with α = 2

	Inst.	m	ϕ D	ϕ R	t	nRS St.	ϕ CT	ϕ W T	ϕ R -ϕ CT
	c101	12 994.53 3371.49 2.00	15	13 3371.45 1930.15	0.04
	c102	11 934.65 3174.41 1.85	12	11 3168.48 1591.99	5.93
	c103	10 851.85 2936.13 2.90	13	12 2887.77 1817.09	48.36
	c104	10 820.17 2780.38 3.39	14	13 2780.38 1940.96	0.00
	c105	11 978.55 3317.31 2.81	18	13 3317.27 1795.43	0.04
	c106	11 936.53 3237.32 3.20	16	13 3174.82 2025.13	62.50
	c107	10 931.43 3228.79 3.13	14	13 3157.55 2045.75	71.24
	c108	10 884.49 3048.56 3.01	10	8 2998.41 1415.12	50.15
	c109	10 837.99 2912.25 3.70	14	12 2840.79 1858.28	71.46
	c201	4 610.38 1392.90 1.67	7	7 1391.66 881.30	1.24
	c202	4 610.63 1392.23 1.98	7	7 1392.23 907.55	0.00
	c203	4 606.59 1383.02 3.31	9	8 1383.03 1226.58	0.00
	c204	4 599.74 1373.40 3.94	9	9 1373.40 1260.51	0.00
	c205	4 603.90 1382.93 2.55	8	8 1382.93 1209.07	0.00
	c206	4 599.62 1373.13 3.06	8	8 1373.13 1249.95	0.00
	c207	4 598.90 1371.48 3.33	8	8 1371.48 1174.30	0.00
	c208	4 593.58 1361.21 3.40	8	8 1359.30 1228.58	1.91
	r101	18 1600.74 768.36 3.17	22	4	768.36	79.57	0.00
	r102	15 1472.88 709.28 3.39	22	8	706.98 173.86	2.30
	r103	13 1197.93 575.01 3.56	18	6	575.01 128.45	0.00
	r104	11 1028.61 463.80 3.73	14	9	462.87 200.14	0.93
	r105	14 1311.36 631.14 3.56	20	12 629.45 254.50	1.68
	r106	13 1224.60 589.10 3.61	17	9	587.81 203.34	1.29
	r107	11 1088.56 493.82 3.25	14	10 489.85 221.40	3.97
	r108	10 995.50 470.33 3.69	15	9	467.88 196.74	2.45
	r109	12 1147.24 530.34 3.68	13	8	527.73 169.17	2.60
	r110	11 1040.83 469.68 4.09	13	8	468.37 172.59	1.30
	r111	11 1046.19 481.89 4.13	14	7	481.25 134.22	0.65
	r112	10 978.38 450.06 3.73	13	7	450.06 152.11	0.00
	r201	3 1250.29 200.05 3.54	7	4	200.05	80.52	0.00
	r202	3 1043.66 135.68 2.84	5	5	135.68	98.19	0.00
	r203	3 884.99 141.60 3.49	7	4	141.60	91.48	0.00
	r204	2 768.80	92.26 4.94	3	3	92.26	71.43	0.00
	r205	3 978.27 146.74 2.61	4	4	146.74	80.93	0.00
	r206	3 930.09 158.11 3.29	6	4	158.11	90.62	0.00
	r207	2 843.00	92.73 4.57	3	2	92.73	36.94	0.00
	r208	2 728.57 102.00 4.14	4	4	102.00	69.20	0.00
	r209	3 863.63 138.18 3.15	4	4	138.18	82.41	0.00
	r210	3 837.98 134.08 2.92	4	4	134.08	80.04	0.00
	r211	2 822.03	90.42 4.65	3	2	90.42	38.31	0.00
	rc101	15 1646.24 625.69 2.93	21	5	625.57 104.82	0.11
	rc102	14 1507.28 572.77 2.67	19	2	572.77	38.42	0.00
	rc103	12 1348.50 512.55 2.91	16	4	512.43	97.24	0.12
	rc104	11 1173.11 445.78 3.51	13	7	445.78 153.27	0.00
	rc105	14 1440.89 547.54 2.69	14	2	547.54	40.00	0.00
	rc106	13 1382.21 525.24 2.86	16	2	525.24	46.42	0.00
	rc107	12 1243.08 472.37 2.89	12	3	472.37	62.00	0.00
	rc108	11 1154.11 438.56 2.58	12	2	438.56	40.00	0.00
	rc201	4 1440.82 201.72 2.69	7	3	201.71	71.46	0.00
	rc202	3 1411.69 155.29 3.80	4	1	155.29	26.49	0.00
	rc203	3 1051.91 147.27 3.60	6	2	147.27	53.46	0.00
	rc204	3 874.10 166.08 3.39	7	4	166.08 102.31	0.00
	rc205	3 1246.91 187.04 4.36	9	4	187.04	84.97	0.00
	rc206	3 1199.77 155.97 3.75	7	2	155.97	44.97	0.00
	rc207	3 987.01 138.87 3.85	5	4	138.18	86.82	0.68
	rc208	3 836.23 150.52 3.71	6	5	150.52 115.82	0.00
	Average			3.31			932.39 529.15	5.91

m

improve the best known solution in 25 instances. Comparing our method with the ALNS presented by Keskin and C ¸atay (2016) the average gap is -0.99%. Likewise, comparing our method with the results obtained by [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF] the average gap is 0.13%. Finally, the average gap of our method versus the hybrid genetic algorithm, presented in Hiermann et al. (2019), is -1.58%.

It is difficult to perform an exact comparison with respect to computational time because of the clear difference between the characteristics of the workstations used by the other authors (see footnotes Table 5). Also note that upper bounds on the number of vehicles are required in our solution method and the numbers of vehicles reported in the BKS were used. Keskin and C ¸atay (2016) algorithm report an average run time of 16.77 min solving the full set of large instances. [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF] ALNS show an average time of 4.31 minutes. HGA from Hiermann et al. (2019) reports an average computational time of 9.96 minutes. Finally, the average time of our Hybrid-ILS is 2.56 minutes. We can see that Hybrid-ILS has high performance compared to the state of the art methods and it is considerably faster with a computer that seems to have less computational power.

Results on the small E-VRPTWsc instances

In this subsection the benefit of visiting customers by walking is analyzed, by comparing the results of the EVRP-TWPR with the results of E-VRPTWsc. EVRP-TWPR best known solutions are taken from the literature, and results for E-VRPTWsc are obtained by executing the hybrid ILS 10 times. The results on small-size instances are shown in Table 6. EVRP The number of vehicles m, number of visited RS st., the total distance performed by the vehicles ϕD, the total time spent at RS ϕR, the percentage gap ∆D for the total distance and the percentage gap ∆R for the total recharging time are presented. Average run-times are presented in column t (in seconds).

Table 6: Results for the Small Instances for the E-VRPTWsc

Results demonstrate important reductions in terms of total charging time and total distance by EVs when the customers are clustered. However, the random and randomclustered instances might represent geographical scenarios where the large distances between customers and RS are not suitable for going by walking, as it can be expected. This reinforces the idea that our variant is focused on a city logistics framework. According to The number of vehicles m, the total distance performed by the vehicles ϕD, the total time spent at RS ϕR, the percentage gap ∆D for the total distance and the percentage gap ∆R for the total recharging time are presented. Average run-times are presented in column t and CPLEX execution time in column t (in seconds).

OM * Out of Memory.

Table 7: Results for small-size instances for the E-VRPTWsc

The results show that the Hybrid-ILS is able to find all the optimal solutions obtained by CPLEX. Likewise, the hybrid ILS always find a solution better than the upper bound provided by CPLEX after being executed for 2 Hours. On average it provides solutions with a GAP of -2.31% in the total distance and a GAP of -1.79% in the total time spent at

In the current version of the problem, only one customer can be visited during the recharging of the EV. One of the research perspectives is to explore the possibility of visiting more than one customer which will bring the concept of secondary routes. Likewise, considering different charging fees for RS and depot are part of our research perspectives. , and bmin i are derived from the forward penalty functions. The new functions and the adapted functions to evaluate a route backward along arc (i, j) from vertex j to vertex i are stated below:

For the following functions (A.8 -A.11) if i ∈ F , e i must be replaced with e i . bmin