Aurélien Froger

Jorge E Mendoza

Ola Jabali

Gilbert Laporte

Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions

published or not. The documents may come

Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions

Introduction

Electric vehicles (EVs) offer a viable alternative to internal combustion engine vehicles (ICEVs), and their market share is steadily increasing, partly as a result of several governmental and private policies. In parallel, the study of problems rising from using EVs by operational researchers has gained importance in recent years, and has given a rise to a new class of problems called electric vehicle routing problems (E-VRPs). These problems differ from the classical vehicle routing problem (VRP) in several ways since the EV technological constraints need to be accounted for. Indeed, the driving autonomy of EVs is limited by their battery capacity, which can be recovered at charging stations (CSs). These are much scarcer than conventional refueling stations for ICEVs, which means that EVs often may need to perform detours to reach a CS. The latter is especially true in the context of mid-haul or long-haul routing [START_REF] Montoya | Electric Vehicle Routing Problems: models and solution approaches[END_REF][START_REF] Schiffer | Logistics Networks with Intermediate Stops: Designing Innovative and Green Solutions[END_REF]. In addition, charging the battery of an EV takes much longer than refilling the tank of an ICEV. Therefore, E-VRPs involve additional decisions such as where and how long to charge.

The charging infrastructure, which includes the CSs used to recuperate energy of EVs, can be either public or private. The use of the former by private companies may result in uncertainty with respect to CS availability.

Therefore, the use of public infrastructure to charge EVs in a routing context has received little attention. We refer the reader to [START_REF] Sweda | Adaptive routing and recharging policies for electric vehicles[END_REF] and [START_REF] Kullman | Dynamic Electric Vehicle Routing with Mid-route Recharging and Uncertain Availability[END_REF] for further details. Since the uncertainty with respect to CS availability is neglected in the majority of the E-VRP literature, we conclude that most of the research on E-VRPs implicitly assumes that the charging infrastructure is privately owned by the EV operator.

We adopt this assumption in this paper as well.

Research on E-VRPs primarily started with the green vehicle routing problem (Green VRP) introduced by [START_REF] Erdogan | A green vehicle routing problem[END_REF]. In the Green VRP, a fleet of homogeneous alternative fuel vehicles must visit a set of customers from a single depot while traveling a minimum total distance, subject to a route duration constraint. The vehicles have a limited fuel tank capacity but can detour to CSs to extend their driving range.

The tank is fully recharged in constant time at each visit to a CS, an assumption that is appropriate for the vehicles operating on alternative fuels such as biofuels and liquid nitrogen. With respect to EVs, such an assumption makes sense if batteries are replaced by fully loaded batteries at swapping stations. To our knowledge, however, battery swapping does not seem to be on the path of massive adoption. Indeed, the company Better Place that offered battery swapping services declared bankruptcy in 20131 and Tesla's only battery swapping station (part of a pilot program) shut down in 2016 2 . Several variants of the Green VRP have been studied. Thus [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF] defined the electric routing problem with time windows and cargo capacity constraints (E-VRPTW). In the E-VRPTW, the batteries are fully recharged at CSs. However, the charging time is no longer assumed to be constant, it is linearly dependent on the state of charge (SoC) of the EV upon its arrival at the CS. Finally, no route duration limit is considered. In the remainder of this paper we use the term SoC to refer to the battery level in the EV.

The above problems assume that whenever a charging activity is performed, it is concluded only when the battery is fully charged. This assumption yields more tractable models, which is partly the reason behind its adoption by earlier versions of E-VRPs. Nonetheless, assuming that the battery is fully charged whenever a charging operation takes place is rather conservative. In view of the fact that the time needed to reach a fully charged battery may be substantial when charging operations occur en route, [START_REF] Felipe | A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges[END_REF] introduced a partial recharging policy. Compared to a full charging policy, a partial charging policy allows savings in both energy consumption and time. The authors assumed a linear charging process and considered CSs with multiple technologies. Each technology has a distinct charging rate, which reflects the variety of chargers one can encounter in practice (e.g., slow, medium, and fast). They also considered capacitated vehicles and route duration constraints. Furthermore, [START_REF] Desaulniers | Exact Algorithms for Electric Vehicle-Routing Problems with Time Windows[END_REF] tackled a version of the E-VRPTW in which the fleet size is fixed and partial charging is allowed.

In reality, battery charging functions of EVs are not linear (see [START_REF] Pelletier | Battery degradation and behaviour for electric vehicles: Review and numerical analyses of several models[END_REF] for further details). In order to model a more realistic charging process, [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] introduced the electric vehicle routing problem with nonlinear charging function (E-VRP-NL). In addition to considering a partial charging policy, the authors proposed a piecewise linear approximation of the charging process. They also showed that previously adopted linear charging functions yield infeasible routes with respect to route maximum-duration limit, when projected on the more realistic nonlinear charging functions. Moreover, even if the routes produced under linear charging functions are feasible, when subjected to nonlinear charging functions they tend to be of poor quality. The objective in the E-VRP-NL is to minimize the total time needed to serve all the customers, which is influenced by the charging decisions because it accounts for the time spent at CSs and for the nonlinear nature of the charging process. Using a temporal-based objective is all the more relevant since CSs may offer multiple technologies which impact the duration of the charging operations. To solve the E-VRP-NL, [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] developed a hybrid metaheuristic that combines an iterated local search heuristic with the solution of a set partitioning problem. Specifically, the first stage of the algorithm aims to build a high-quality and diverse pool of routes, while the second stage assembles solutions by selecting a subset of routes from the pool.

Research on the E-VRPs is now moving towards more accurate representations of EV characteristics which pose new modeling and algorithmic challenges. One of the most challenging issues when solving E-VRPs is making good charging decisions. In general, charging operations influence feasibility and in many cases they also increase costs, due to detours to CSs and the time spent at those CSs. Depending on the underlying E-VRP, deciding where to charge and how much to charge can be very complex. To better deal with these decisions, some authors have studied the fixed route vehicle charging problem (FRVCP). Given a fixed sequence of customers, the FRCVP consists in optimally inserting charging operations within it. As such, the FRVCP can be viewed as a subproblem of E-VRPs, which is pivotal to many heuristics that must compute the cost of a given route. [START_REF] Montoya | A multi-space sampling heuristic for the green vehicle routing problem[END_REF] solved a FRVCP with the full charging policy and constant charging time.

They formulated the problem as a constrained shortest path problem and solved it using the Pulse algorithm [START_REF] Lozano | On an exact method for the constrained shortest path problem[END_REF]. [START_REF] Roberti | The electric traveling salesman problem with time windows[END_REF] described a labeling algorithm to solve the FRVCP arising in the case where all the CSs have the same charging technology, the charging function is linear, and partial charging is allowed. [START_REF] Hiermann | The electric fleet size and mix vehicle routing problem with time windows and recharging stations[END_REF], [START_REF] Schiffer | An adaptive large neighborhood search for the location-routing problem with intra-route facilities[END_REF], [START_REF] Hiermann | Routing a mix of conventional, plug-in hybrid, and electric vehicles[END_REF] considered the FRVCP as a sub-problem of the E-VRP-TW with single charging technology and linear charging function. They tackled a special version of the problem where at most one CS can be inserted between any pair of non-CS nodes. Note that [START_REF] Hiermann | The electric fleet size and mix vehicle routing problem with time windows and recharging stations[END_REF] also assumed a full charging policy. Finally, [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] considered the FRVCP as a subproblem of the E-VRP-NL. They presented a mixed integer linear programming (MILP) model in which they assume that an EV can visit at most one CS between any pair of non-CS nodes and a heuristic.

The scientific contribution of this paper is twofold. First, we propose two new mixed integer linear programming (MILP) models for the E-VRP-NL, which is central to the E-VRP literature in that it realistically models the EV charging operations via piecewise linear approximations of nonlinear functions, while accounting for various charging technologies and allowing partial recharge. The first model results from replacing the node-based tracking of the time and SoC from the [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] model by an arc-based tracking of these quantities. We demonstrate that this new formulation provides an equal or better LP relaxation bound.

To avoid replicating the CS nodes, as was the case in both previously discussed models, we then define the problem on a multigraph, which allows us to derive a second model for the E-VRP-NL. We explain why this formulation better fits the problem from a modeling perspective, and we also demonstrate its superiority through extensive computational experiments. Building on the second new E-VRP-NL model, the second scientific contribution of this paper lies in the development of three new algorithms for the FRVCP: a direct approach that solves a MILP formulation of the problem, a heuristic algorithm, and an exact labeling algorithm. These three algorithms are compared through an extensive set of computational experiments.

The remainder of this paper is organized as follows. We formally describe the E-VRP-NL in Section 2. In Section 3, we present a new MILP model for the problem based on its "classical" definition on a simple digraph.

In Section 4, we introduce an alternative MILP model for the E-VRP-NL using the concept of paths visiting sequences of CSs. In Section 5, we describe a new MILP model, a heurisitic and an exact labeling algorithm for the FRVCP. We present computational results in Section 6, followed by conclusions in Section 7.

2 The electric vehicle routing problem with nonlinear charging function

The electric vehicle routing problem with nonlinear charging function, introduced by [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF], is defined as follows. Let I be the set of customers and let F be the set of CSs at which the vehicles can fully or partially recharge their battery. Each customer i I has a service time g i . The customers are served by an homogeneous fleet of EVs, which is unlimited. Each EV has a battery of capacity Q (expressed in kWh) and a maximum tour duration T max . At the start of the planning horizon, the EVs are located at a single depot (denoted by node 0), from which they leave fully charged. Traveling from a location i (the depot, a customer, or a CS) to another location j incurs a driving time t ij ¥ 0 and an energy consumption e ij ¥ 0. The triangular inequality holds for both the driving time and the energy consumption. Each CS i has a charging technology (e.g., slow, moderate, fast) associated with a piecewise linear concave charging function Φ i p∆q that maps, for an empty battery, the time ∆ spent charging at i to the SoC of the vehicle when it leaves i. If q is the SoC of an EV when arriving at i, then the SoC of the EV upon departing from i after charging during a time equal to ∆ is given by Φ i p∆ Φ ¡1

i pqqq. We denote by B i t0, . . . , b i u the ordered set of breakpoints of the piecewise linear approximation of the charging curve of CS i. Let c ik and a ik be the charging time and the SoC of breakpoint k B i of the CS i. Let ρ ik be the slope of the segment joining the points pc i,k¡1 , a i,k¡1 q and pc ik , a ik q (i.e., ρ ik pa ik ¡ a i,k¡1 q{pc ik ¡ c i,k¡1 q. Figure 1 shows the piecewise linear approximation for a CS i yielding a power of 22 kW charging a vehicle equipped with a 16 kWh battery.

Figure 1: Piecewise linear approximation for a CS yielding a power of 22 kW charging a 16 kWh battery.

Adapted from [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF].

Feasible solutions to the E-VRP-NL satisfy the following conditions: 1) each customer is served exactly once by a single vehicle, 2) each route starts and ends at the depot, 3) each route satisfies the maximum-duration limit T max , and 4) each route is energy-feasible (i.e., the battery level of an EV upon arriving at a location or departing from it lies between 0 and Q). The objective of the E-VRP-NL is to minimize the total time needed to serve all customers, including driving, service, and charging time.

3 Two CS replication-based models for the E-VRP-NL

The E-VRP-NL is defined on a digraph G pV, Aq, where V t0uI F I is the set of nodes and A tpi, jq : i, j V, i $ ju is the set of arcs connecting nodes of V . Node 0 represents the single depot. The set F I contains β i copies of each CS i F (i.e., |F I | °iF β i ). The value of β i corresponds to an upper bound on the number of times i can be visited. These copies are introduced for modeling convenience: each visit of an EV to a CS is modeled as a visit to a distinct copy of the CS, i.e., each node in F I can be restricted to be visited at most once. Therefore, the use of CS copies facilitates the tracking of the SoC and of the time. We denote by F I i F I the set containing the β i copies of CS i (i.e., |F I i | β i and F I iF F I i ). We assume that F I i is an ordered set and that its elements are numbered from 1 to β i . In the remainder of this paper, depending on the context, we refer to an element of F I or F I i as a CS copy. Furthermore, the graph G is not necessary complete, since some arcs cannot be traversed due to the SoC constraints. Therefore, we apply a preprocessing procedure to reduce its size (see A).

For the sake of completeness and clarity, we recall in §3.1 the model introduced in [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF].

In §3.2 we introduce an alternative formulation, which uses a different approach to model the time and SoC tracking.

First model for the E-VRP-NL

The formulation of [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] involves the following decision variables. The binary variable x ij is equal to 1 if and only if an EV travels on arc pi, jq A. The continuous variables τ j and y j track the time and SoC of the EV when it departs from node j V . The continuous variables q i and q i specify (according to the piecewise linear approximation of the charging function of i) the SoC of an EV when it enters and leaves CS copy i F I . The variables c i and c i are the scaled start and end times for charging an EV according to the charging function of CS copy i. The continuous variable ∆ i c i ¡ c i represents the duration of the charging operation performed at CS copy i. The binary variables w ik and w ik are equal to 1 if and only if the SoC is between a i,k¡1 and a ik , with k B i zt0u, when the EV enters and leaves CS copy i, respectively. Finally, the continuous variables λ ik and λ ik are the coefficients associated with the breakpoint pc ik , a ik q in the piecewise linear approximation, when the EV enters and leaves CS copy i, respectively. Specifically, the λ ik variables enable the expression of pc i , q i q as a convex combination of the breakpoints pc ik , a ik q, where k B i . Similarly, the λ ik variables enable the expression of pc i , q i q as a convex combination of the breakpoints pc ik , a ik q, where k B i . The [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] model for the E-VRP-NL, denoted as rF CS rep N ode s, is as follows:

rF CS rep N ode s minimize pi,jqA t ij x ij i F I ∆ i (1)
subject to pi,jqA

x ij 1 i I (2)
pi,jqA

x ij ¤ 1 i F I (3) pj,iqA x ji ¡ pi,jqA x ij 0 i V (4) e ij x ij ¡ p1 ¡ x ij qQ ¤ y i ¡ y j ¤ e ij x ij p1 ¡ x ij qQ pi, jq A : i V, j I (5) e ij x ij ¡ p1 ¡ x ij qQ ¤ y i ¡ q j ¤ e ij x ij p1 ¡ x ij qQ pi, jq A : i V, j F I (6) y i ¥ e i0 x i0 pi, 0q A (7)
y i q i i F I (8) y 0 Q (9) q i ¤ q i i F I (10) q i ķBi λ ik a ik i F I (11) c i ķBi λ ik c ik i F I (12) ķBi λ ik ķBizt0u w ik i F I (13) ķBizt0u w ik pi,jqA x ij i F I (14) λ i0 ¤ w i1 i F I (15) λ ik ¤ w ik w i,k 1 i F I , k B i zt0, b i u (16) λ ibi ¤ w ibi i F I (17) q i ķBi λ ik a ik i F I (18) c i ķBi λ ik c ik i F I (19) ķBi λ ik ķBizt0u w ik i F I (20) ķBizt0u w ik pi,jqA x ij i F I (21) λ i0 ¤ w i1 i F I (22) λ ik ¤ w ik w i,k 1 i F, I k B i zt0, b i u (23) λ ibi ¤ w ibi i F I (24) ∆ i c i ¡ c i i F I (25) τ i pt ij g j qx ij ¡ T max p1 ¡ x ij q ¤ τ j pi, jq A : i V, j I (26)
τ i ∆ j t ij x ij ¡ pS max T max qp1 ¡ x ij q ¤ τ j pi, jq A : i V, j F I (27) τ j t j0 ¤ T max j V (28) τ 0 0 (29) τ l ¡ ∆ l ¥ τ j ¡ ∆ j i F, j, l F I i , j l (30)
ph,jqA

x hj ¥ ph,lqA x hl i F, j, l F I i , j l (31)
x ij t0, 1u pi, jq A (32)

τ i ¥ 0, y i ¥ 0 i V (33) w ik t0, 1u, w ik t0, 1u i F I , k B i zt0u (34) λ ik ¥ 0, λ ik ¥ 0 i F I , k B i (35) q i ¥ 0, q i ¥ 0, ∆ i ¥ 0, c i ¥ 0, c i ¥ 0 i F I . ( 36
)
The objective function (1) minimizes the total driving and charging time. Constraints (2) ensure that each customer is served once, while constraints (3) state that every CS copy is visited at most once. Constraints (4) impose the flow conservation conditions. Constraints ( 5) and ( 6) track the SoC of an EV at each node.

Constraints (7) ensure that if an EV travels between a customer or CS copy and the depot, then it has sufficient energy to reach its destination. Constraints (8) reset the battery tracking to q i upon departure from CS copy i. Constraint (9) states that the SoC of an EV leaving the depot is Q. Constraints (10) couple the SoC of an

EV upon arrival at a CS with its SoC at departure. Constraints ( 11)-( 17) define the SoC and the charging time of an EV upon arrival at CS copy i, based on the piecewise linear approximation of its charging function.

Similarly, constraints ( 18)-( 24) define the SoC and charging time of an EV upon departure from CS copy i.

Constraints (25) define the time spent at any CS copy. Constraints ( 26) and ( 27) track the departure time from each node. The parameter S max max iF tc ibi u corresponds to the maximum charging time, i.e., the time to fully charge an empty battery at the slowest CS. Constraints ( 28) and ( 29) ensure that the EVs return to the depot no later than T max . Even if they are not necessary, constraints (30) and ( 31) avoid potential symmetries between the CS copies. These constraints ensure that the copies of CS i are visited in the order in which they appear in F i (i.e, a charging operation at j F I i must start after a charging operation at l F I i if j ¡ l). Finally, constraints (32)-( 36) define the domains of the decision variables.

In order to strengthen the previous formulation, we add the following valid inequalities:

y i ¥ min lF t0u te il u i V zt0u. ( 37
)
These inequalities state that the SoC of an EV upon departing from a node must be sufficient to reach the depot or the nearest CS.

Second model for the E-VRP-NL

We now introduce a modeling alternative for the E-VRP-NL. Formulation rF CS rep N ode s tracks the time and the SoC of the EVs with node-indexed variables. We therefore refer to constraints ( 5)-( 9), ( 26)-( 29), ( 31

t ij x ij i F I ∆ i subject to
(2) ¡ (4), (10) ¡ (25), ( 31), ( 32), (34) ¡ (36)

y 0j Qx 0j p0, jq A (38)
pi,jqA

y ij ¡ pi,jqA e ij x ij pj,lqA y jl j I (39)
pi,jqA

y ij ¡ pi,jqA e ij x ij q j j F I (40) pj,lqA y jl q j j F I (41) y ij ¤ ¢ Q ¡ min lF t0u te li u x ij pi, jq A (42) pi,jqA pτ ij pt ij g j q x ij q pj,lqA τ jl j I (43) pi,jqA pτ ij t ij x ij q ∆ j pj,lqA τ jl j F I (44) τ ij ¤ pT max ¡ t ij ¡ g j ¡ t j0 q x ij pi, jq A : j I (45) τ ij ¤ T max ¡ t ij ¡ ∆ min j ¡ t j0 ¨xij pi, jq A : j F I (46) pj,hqA τ jh ¡ ∆ j ¥ pl,hqA τ lh ¡ ∆ l i F, j, l F I i : j l (47) τ ij ¥ 0, y ij ¥ 0 pi, jq A. (48)
Constraints (38) state that the EVs leave the depot fully charged. Constraints (39) track the SoC of an EV at each customer. Similarly, constraints (40) track the SoC of an EV when it arrives at a CS copy, and constraints (41) track the SoC of the EV when it leaves a CS copy. Constraints (42) couple the variables y ij and x ij . Constraints ( 43) and ( 44) track the departure time from each customer and CS copies, respectively.

Constraints ( 45) and ( 46) couple the τ ij and x ij variables. Specifically, if an EV traverses an arc pi, jq, then its departure time must guarantee that it returns to the depot without exceeding the tour maximum-duration limit. The parameter ∆ min j represents the minimum duration of a charging operation at j F I , and is equal to the time needed to charge the energy quantity min l,l I V ztju:pl,jqApj,l I qA te lj e jl I ¡ e ll Iu if the EV has an empty battery. It is a lower bound on the time spent charging to recover the energy consumed to make the detour to j. Constraints (47) break symmetries created by the introduction of CS copies. Constraints (48) define the domains of the newly introduced decision variables.

In order to strengthen the previous formulation, we add the following valid inequalities:

y ij ¥ ¢ e ij min lF t0u te jl u x ij pi, jq A : i $ 0. ( 49
)
Constraints ( 49) state that if an EV traverses arc pi, jq, then its SoC when leaving i must be sufficient to traverse the arc and then reach the closest CS or depot.

According to the experiments conducted by [START_REF] Ascheuer | Solving the Asymmetric Travelling Salesman Problem with time windows by branch-and-cut[END_REF] on the asymmetric traveling salesman problem with time windows, MILP solvers tend to perform better on arc-based tracking constraints, similar to (39)-( 48), than on those based on node-based tracking constraints when there are few time windows. Therefore, there is ground to believe that reformulating the time and SoC tracking in rF CS rep N ode s using this new model may yield a better computational performance. Proposition 3.1 supports this claim for the two previously presented formulations.

Proposition 3.1. Let P arc denote the polyhedron associated with the linear relaxation of rF CS rep Arc s, and let Let P node denote the polyhedron associated with the linear relaxation of rF CS rep N ode s. Then P arc P node .

Proof. See B.

Valid inequalities for the first two models

We introduce a new set of valid inequalities to strengthen the previous rF CS rep N ode s and rF CS rep Arc s formulations.

Let ec min l be the minimum amount of energy that an EV must charge at CS copy l. This corresponds to the charge needed to recover the minimum energy consumed to make a detour to l. Specifically, ec min l min i,jV ztlu:pi,lqApl,jqA te il e lj ¡e ij u. Assuming an empty battery upon arrival at the CS, we can therefore state the following valid constraints:

∆ i ¥ Φ i pec min i q pi,jqA x ij i F I . ( 50
)
4 Path-based model without CS copies for the E-VRP-NL

A major drawback of the two previous formulations defined on a digraph is the need to replicate CS nodes.

Indeed, to ensure that no optimal solutions are cut off, the number of copies may have to be very large, which affects the solution time. We can find examples where this value is equal to 4|I| (see C). However, working with an order of |I| copies for each CS is impractical. Despite this limitation, the vast majority of MILP models for E-VRPs listed in the literature are based on the CS-replication paradigm. To our knowledge, only [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] proposed a procedure to set the number of CS copies needed to be used in their model. They used the same number β of CS copies for every CS. Starting with β 0, they solved the MILP formulation to optimality with a time limit of T (set to 100 hours in their experiments); β was then incremented and the resulting MILP was solved. If β ¥ 1 and either the time limit was reached or the optimal solution obtained with β l had the same value as the one obtained with β l ¡1, they stopped the procedure and set β l ¡1.

This procedure is not optimal when the time limit is reached. Moreover, even when an instance is solved to optimality, this procedure is not optimal (see D for examples).

We propose an alternative model that does not use CS copies. This new model is based on enumerating the sequence of CS nodes that can be visited by an EV between each pair of customer or depot nodes. Given two nodes (customer or depot) i and j, we call a CS path (CSP) a simple path starting from i visiting a sequence of CSs (possibly empty) and ending at j.

CSPs were introduced by [START_REF] Andelmin | Electric vehicle routing with realistic recharging models[END_REF] as "e-paths" in the context of the E-VRPTW with partial recharging. They were successfully used by [START_REF] Bartolini | An exact algorithm for the green vehicle routing problem[END_REF] (as "refuel paths") and constituted the main building block of an exact algorithm for the Green VRP. They were also applied by Roberti and Wen (2016) (as "recharging paths") to the electric traveling salesman problem with time windows. [START_REF] Koç | The green vehicle routing problem: A heuristic based exact solution approach[END_REF] and [START_REF] Leggieri | A practical solution approach for the green vehicle routing problem[END_REF] employed a similar concept. In the two latter studies, however, only one CS can be visited between any pair of nodes. This assumption may cut off optimal solutions to the Green VRP (see [START_REF] Bartolini | An exact algorithm for the green vehicle routing problem[END_REF] for an example). Because in the Green VRP the vehicles are fully charged in constant time upon visiting a CS, the SoC of an EV after traversing a refuel path and the time spent traveling on a path can then be computed a priori. Allowing partial charging combined with the non-linearity of the charging functions and/or the consideration of multiple technologies (any one of the two latter is a sufficient condition), like in the E-VRP-NL, necessitates deciding how much to charge at each CS and thus tracking the time and SoC throughout the CSPs.

We develop the concept of CSP for the E-VRP-NL. We redefine the problem on a directed multigraph r G p r V , r Aq, where r V t0u I, and r

A is the set of arcs connecting the nodes of r V . More specifically, an arc in à represents a CSP p, starting in oppq Ṽ and ending in dppq Ṽ . Let n p denote the number of CSs in p and let L p t0, 1, ¤ ¤ ¤ , n p ¡ 1u be the ordered set of CS positions in p. We denote by µ p plq the CS at position l L p . We note that L p r indicates that p does not visit any CS between oppq and dppq. Let also e p and t p be the energy consumption and the driving time associated with CSP p P . Given two nodes i, j r

V , we define P ij as the set of CSPs connecting i to j, and we define P i,j r V ,i$j P ij as the set of all CSPs.

Third model for the E-VRP-NL

Our path-based formulation of the E-VRP-NL involves the following decisions variables. The binary variable

x p is equal to 1 if and only if an EV travels on CSP p P . The continuous variables τ p and y p track the time and SoC of an EV when it departs from node oppq to node dppq on CSP p. If no vehicle travels on the CSP, then both variables take the value of 0. For a path p, the continuous variables q pl and q pl specify the SoC of an EV when it enters and leaves µ p plq, respectively. We model the piecewise linear approximation of the charging function following the same principle as in the first two models. The variables c pl and c pl are the scaled start and end times for charging an EV according to the charging function of CS µ p plq. The continuous variable ∆ pl represents the duration of the charging operation performed at µ p plq. The binary variables w plk and w plk are equal to 1 if and only if the SoC is between a µpplq,k¡1 and a µpplq,k , with k B i zt0u, when the EV enters and leaves CS µ p plq, respectively. The continuous variables λ plk and λ plk are the coefficients associated with the breakpoint pc µpplq,k , a µpplq,k q in the piecewise linear approximation, when the EV enters and leaves CS µ p plq, respectively. The path-based formulation of the E-VRP-NL, denoted as rF P ath Arc s, is as follows:

rF P ath Arc s min pP ¤ ¥ t p x p ļ Lp ∆ pl (51) subject to j r V ,i$j p Pij x p 1 i I (52) j r V ,i$j p Pji x p ¡ j r V ,i$j p Pij x p 0 i r V (53) j r V ,j$i p Pji ¤ ¥ y p ¡ e p x p ļ Lp pq pl ¡ q pl q j r V ,j$i p Pij y p i I (54)
y p ¡ e oppq,µpp0q x p q p0 p P : |L p | $ 0 (55) q p,l¡1 ¡ e µppl¡1q,µpplq x p q pl p P, l L p zt0u

(56) c pl , c pl , q pl , q pl , ∆ pl ¥ 0 p P, l L p (79)

y p Qx p i r V zt0u, p P 0i (57)
y p ¡ e p x p ļ Lp pq pl ¡ q pl q ¥ 0 i r V zt0u, p P i0 (58)
y p ¤ Qx p p P (59) j r V zt0u,j$i p Pji τ p j r V ,j$i p Pji ¤ ¥ t p x p ļ Lp ∆ pl g i j r V ,j$i p Pij τ p i I (60)
τ p ļ Lp ∆ pl ¤ T max ¡ t p ¡ g dppq ¡ t
w plk t0, 1u, w plk t0, 1u p P, l L p , k B i zt0u

(80)

λ plk ¥ 0, λ plk ¥ 0 p P, l L p , k B i . (81)
The

Dominated paths

Theoretically the set r A contains a very large number of arcs that grows quickly with the number of CSs and the number of customers. However, many of these arcs will not be part of an optimal solution. Hence, this section is devoted to identifying CSPs that can be eliminated using dominance rules. Defining these rules is a complex task since we need to account for the partial charging policy, the various charging technologies, and the piecewise linearity of the charging functions.

We first introduce the concept of recharging path (RP). We call a RP a CSP for which the SoC of the EV upon departing from each CS -except the last visited one -is known. Indeed, the SoC required at the end of a RP is unknown, because the destination of the EV after traveling the RP is itself unknown. We denote by pp, φq a RP associated with CSP p where φ is a vector specifying the SoC of the EV upon departing from each CS visited in position l L p ztn p ¡ 1u. Specifically, φ l is the SoC of the EV upon departing from CS µ p plq, with l L p ztn p ¡ 1u. We refer to φ as the SoC target vector.

Let us now assume that the SoC of the EV when departing from the origin of each CSP is known. We denote this SoC by q and refer to it as the initial SoC. A SoC target vector is said to be feasible for a CSP with respect to q if it leads to an energy-feasible path, that is a path in which the SoC is non-negative and less than Q when the vehicle enters or leaves any of its nodes. Let CS be a CS other than the last on a CSP (i.e., a CS visited in position l L p ztn p ¡ 1u when traveling a CSP p). Since the number of SoC targets at the CSs of a RP is infinite, there is an unbounded number of RPs associated with a given CSP that contains two or more CSs. However, [START_REF] Zündorf | Electric vehicle routing with realistic recharging models[END_REF] showed that for a given sequence of CS visits between two non-CS nodes, we only need to consider a limited number of SoC targets for the CSs. This is due to the piecewise linearity of the charging functions. Specifically, the number of SoC targets that should be considered at a CS is restricted by the number of breakpoints of its charging function. The possibilities are the following: 1) One can leave a CS with just enough energy to arrive to the next CS with an empty battery. 2) One can also stop charging at a CS at any breakpoint of the charging function if the resulting SoC is larger than the SoC of the EV upon arriving at the CS. Therefore, we apply these results to RPs, and by doing so we limit the number of RPs associated with every CSP.

We now define a dominance rule between two RPs. To this end, we exploit the concept of SoC functions introduced by [START_REF] Zündorf | Electric vehicle routing with realistic recharging models[END_REF] to solve an electric vehicle shortest path problem, which consists in finding a minimum duration route satisfying battery constraints and allowing the EV to recharge en route. A SoC function SoC q pp,φq , for a RP pp, φq and a given initial SoC q, maps the duration of the path t to the final SoC of the EV at destination dppq.

Since we assume that all charging functions are piecewise linear, a SoC function can be defined by an ordered set of supporting points tpt 1 , y 1 q, ¤ ¤ ¤ , pt k , y k qu sorted in non-decreasing order of duration. These points are derived from the breakpoints of the charging function of the last CS in a RP according to p, φ and q.

The value t 1 represents the minimal time required by an EV to travel on the RP according to the battery constraints (no undercharging or overcharging is allowed and the SoC is sufficient to travel on each arc) and q.

A SoC function represented by tpt 1 , y 1 q, ¤ ¤ ¤ , pt n , y n qu is defined as follows:

SoC q pp,φq ptq 6 9 9 9 8 9 9 9 7 Figure 2 depicts a SoC function. In this example, the first supporting point of the curve is (1, 1.5). This means that the minimum duration to travel on the RP is equal to 1, and in this case, the SoC at the destination of the EV is equal to 1.5. If the SoC needed at dppq is 3, the duration of the corresponding path is equal to 2 (the duration is increased due to a larger amount of energy charged at a CS previously visited by the RP). A dominance rule between RPs having the same origin and destination can then be defined as follows:

¡V if 0 ¤ t t 1 pt ¡ t k¡1 qpy k ¡ y k¡1 q t k y k¡1 if t k¡1 ¤ t t k , k t2, 3, ..., nu y n if t ¥ t n .
Definition 4.1. Let pp 1 , φ 1 q and pp 2 , φ 2 q be two RPs with the same origin o opp 1 q opp 2 q and the same destination d dpp 1 q dpp 2 q, and let q be the initial SoC of the EV at o. RP pp 1 , φ 1 q dominates RP pp 2 , φ 2 q with respect to q if SoC q pp1,φ1q ptq ¥ SoC q pp2,φ2q ptq for every t ¥ 0 and there exists t I ¥ 0 such that SoC q pp1,φ1q pt I q ¡ SoC q pp2,φ2q pt I q. We denote this dominance by pp 1 , φ 1 q ¡ q pp 2 , φ 2 q.

This definition states that a RP dominates another only if for every possible duration the EV reaches the destination with higher or equal SoC. The dominance of RP pp 1 , φ 1 q over RP pp 2 , φ 2 q can be established by comparing the value of SoC q pp1,φ1q and SoC q pp2,φ2q for all abscissas of the supporting points of both functions.

For convenience, we also provide an equivalent definition to Definition 4.1 based on the comparison of the inverse of the SoC function. Let Θ q pp,φq and Θ q pp,φq be the maximum and minimum SoC that can be achieved at the destination of RP pp, φq if the initial SoC is equal to q. The inverse of function SoC q pp,φq , denoted pSoC q pp,φq q ¡1 , exists only in the interval rΘ q pp,φq , Θ q pp,φq s. We define its extension, denoted InvSoC q pp,φq , over the interval r0, Qs as follows:

InvSoC q pp,φq pqq 6 9 9 8 9 9 7 pSoC q pp,φq q ¡1 pΘ q pp,φq q if 0 ¤ q Θ q pp,φq pSoC q pp,φq q ¡1 pqq if Θ q pp,φq ¤ q ¤ Θ q pp,φq

V if Θ q pp,φq q ¤ Q . ( 83
)
Note that since it is not possible to achieve a SoC 0 ¤ q Θ q pp,φq , the first case assumes that the time required to reach such a q is equivalent to the time needed to reach Θ q pp,φq . An equivalent definition to Definition of 4.1 then reads:

Definition 4.2. Let pp 1 , φ 1 q and pp 2 , φ 2 q be two RPs with the same origin o opp 1 q opp 2 q and the same destination d dpp 1 q dpp 2 q, and let q be the initial SoC of the EV. RP pp 1 , φ 1 q dominates RP pp 2 , φ 2 q with respect to q if InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq for every q and there exists qI such that InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. This is also denoted by pp 1 , φ 1 q ¡ q pp 2 , φ 2 q.

This definition states that a RP dominates another only if for every possible SoC at the destination it has a shorter or equal duration. As in Definition 4.1, the dominance of RP pp 1 , φ 1 q over RP pp 2 , φ 2 q can be established by comparing the value of InvSoC q pp1,φ1q and InvSoC q pp2,φ2q for all abscissas of the supporting points of both functions.

Definition 4.3. CSP p 1 is said to be non-dominated if there exists at least one RP pp 1 , φ 1 q and an initial SoC q, such that for every RP pp 2 , φ 2 q with the same origin o opp 1 q opp 2 q and the same destination d dpp 1 q dpp 2 q we have pp 2 , φ 2 q £ q pp 1 , φ 1 q (i.e., pp 2 , φ 2 q does not dominate pp 1 , φ 1 q for an initial SoC q). The initial SoC q can be restricted for a RP pp, φq to take its value in the interval rq M IN pp,φq , q M AX pp,φq s. Let pp,φq plq denote the SoC upon arrival at CS µ p plq for RP pp, φq (with l L p zt0u). This quantity does not depend on the initial SoC q. The values of q M IN pp,φq and q M AX pp,φq can then be defined as follows: The value q M IN pp,φq corresponds to the minimum SoC needed at the origin to make RP pp, φq energy-feasible. The first element in the definition of q M AX pp,φq corresponds to the maximum SoC achievable at the origin of the RP.

q M IN
The second element considers that if the initial SoC allows the EV to reach a CS with a SoC larger than or equal to that achievable by visiting previous CSs, then the RP that goes directly from the origin to that CS dominates the current RP (because the driving time and energy consumption satisfy the triangular inequality).

It is impractical to compute the non-dominated RPs for every possible initial SoC using Definition 4.3, since there is an infinite number of possible initial SoCs (except at the depot). To resolve this issue, we show in Lemma 4.1 how given a SoC function and an initial SoC q, we are able to compute the SoC functions for other initial SoCs q I for a given RP. Based on this result, we demonstrate in Proposition 4.1 how the dominance relationship between two different RPs with the same origin and destination can be efficiently verified.

Lemma 4.1. Let pp, φq be a RP and let q rq M IN pp,φq , q M AX pp,φq s. If the initial SoC at oppq is equal to q I with q ¤ q I ¤ q M AX pp,φq , then the inverse of the SoC function InvSoC q I pp,φq can be computed as follows:

• If n p $ 0 (i.e., p contains at least one CS) d q r0, Qs, InvSoC q I pp,φq pqq 5

InvSoC q pp,φq pqq ¡ λ p qq I if q ¥ q I ¡ e p InvSoC q
pp,φq pq I ¡ e p q ¡ λ p qq I if q q I ¡ e p (84)

where λ p qq I Φ ¡1 µpp0q pq I ¡ e oppq,µpp0q q ¡ Φ ¡1 µpp0q pq ¡ e oppq,µpp0q q corresponds to the difference in time spent charging q I ¡ q at the first CS of the path.

• If n p 0 (i.e., p does not contain any CS) d q r0, Qs, InvSoC q I pp,φq pqq 5

t p if q I ¡ e p ¥ q V if q I ¡ e p q (85) Proof. See B.
For each pair of values q 1 and q 2 such that q M IN pp1,φ1q ¤ q 1 ¤ q M AX pp1,φ1q and q M IN pp2,φ2q ¤ q 2 ¤ q M AX pp2,φ2q , we define two functions, denoted by T pp1,φ1q pp2,φ2q and T pp1,φ1q pp2,φ2q . Intuitively, these functions represent the maximum and minimum differences in time for the two RPs to reach the same SoC at the destination. They are defined as follows:

T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q 6 8 7

max qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tInvSoC q1 pp1,φ1q pqq ¡ InvSoC q2 pp2,φ2q pqqu n p1 ¡ 0 n p2 ¡ 0 max qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tt p1 ¡ InvSoC q2 pp2,φ2q pqqu n p1 0 n p2 ¡ 0 T pp1,φ1q
pp2,φ2q pq, q 1 , q 2 q 6 8 7

min qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tInvSoC q1 pp1,φ1q pqq ¡ InvSoC q2 pp2,φ2q pqqu n p1 ¡ 0 n p2 ¡ 0 min qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tInvSoC q1 pp1,φ1q pqq ¡ t p2 u n p1 ¡ 0 n p2 0
where Ω pp1,φ1q pp2,φ2q pqq is an interval equal to rmaxt0, q ¡ minte p1 , e p2 uu, mintΘ q1 pp1,φ1q , Θ q2 pp2,φ2q us. We note that for a CSP with n p 0, since there are no charging decisions involved, there is a single RP. This RP cannot be dominated for every possible value of the initial SoC by other RPs, which are associated with other CSPs, because of the triangular inequality. We therefore never verify if a RP p with n p 0 is dominated. However, such a RP may dominate other RPs. In this case, the initial SoC of RP p with n p 0 must be such that the SoC at the destination is larger than the maximum SoC attainable at the destination by the other RPs. We also note that given two RPs pp 1 , φ 1 q and pp 2 , φ 2 q with the same origin and the same destination, such that q M IN pp1,φ1q ¤ q M IN pp2,φ2q , the dominance for these two RPs must be verified only for q ¥ q M IN pp2,φ2q . Finally, we note that for a RP pp, φq such that n p ¡ 0, reaching the maximum SoC at the destination is established when the EV is fully charged at the last CS of p. We note that this is independent of the value of the initial SoC q.

Hence for simplicity we denote Θ q pp,φq by Θ pp,φq . Based on these observations, the dominance test between the two RPs (where q M IN pp1,φ1q ¤ q M IN pp2,φ2q ) needs to be carried for q rQ pp1,φ1q pp2,φ2q , Q pp1,φ1q pp2,φ2q s (when the first element is less than or equal to the second element) where Q pp1,φ1q pp2,φ2q and Q pp1,φ1q pp2,φ2q are defined as follows:

Q pp1,φ1q pp2,φ2q 6 9 9 8 9 9 7

q M IN pp2,φ2q if n p1 ¥ 1 n p2 ¥ 1 maxtq M IN pp2,φ2q , Θ pp2,φ2q e opp1q,dpp1q u if n p1 0 n p2 ¥ 1 maxtq M IN pp2,φ2q , Θ pp1,φ1q e opp2q,dpp2q u if n p1 ¥ 1 n p2 0 Q pp1,φ1q pp2,φ2q mintq M AX pp1,φ1q , q M AX pp2,φ2q u Observe that if Q pp1,φ1q
pp2,φ2q ¡ Q pp1,φ1q pp2,φ2q , then we cannot derive any dominance relationship. Moreover, if n p1 ¡ 0 and Θ pp1,φ1q Θ pp2,φ2q , then the EV can have a higher SoC after traversing pp 2 , φ 2 q compared to the one obtained after traversing pp 1 , φ 1 q, which means that there does not exist q such that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q. Similarly, if n p2 ¡ 0 and Θ pp2,φ2q Θ pp1,φ1q , then there does not exist q such that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q. The following result allows us to compute intervals for the initial SoC for which there may exist a dominance between RPs.

Proposition 4.1. Let pp 1 , φ 1 q and pp 2 , φ 2 q be two RPs with the same origin o opp 1 q opp 2 q and the same destination d dpp 1 q dpp 2 q such that q M IN pp1,φ1q ¤ q M IN pp2,φ2q . For k t1, 2u, let q k be a parameter equal to Q pp1,φ1q pp2,φ2q

if n p k 0 and q M IN pp k ,φ k q otherwise. If Q pp1,φ1q pp2,φ2q ¤ Q pp1,φ1q
pp2,φ2q , the following statements hold:

• Case 1: n p1 0 or Θ pp1,φ1q ¥ Θ pp2,φ2q :

if for an initial SoC q rQ pp1,φ1q pp2,φ2q , minte p1 , e p2 us, we have T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¤ ∆ p1 pqq ¡ ∆ p2 pqq and T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ∆ p1 pqq ¡ ∆ p2 pqq, then we can state that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q if for an initial SoC q rminte p1 , e p2 u, Q pp1,φ1q pp2,φ2q ss, we have maxtT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1qu ¤ ∆ p1 pqq¡∆ p2 pqq and mintT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1qu ∆ p1 pqq¡∆ p2 pqq, then we can state that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q

• Case 2: n p2 0 or Θ pp2,φ2q ¥ Θ pp1,φ1q

if for an initial SoC q rQ pp1,φ1q pp2,φ2q , minte p1 , e p2 us, we have T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¥ ∆ p1 pqq ¡ ∆ p2 pqq and T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we can state that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q.

if for an initial SoC q rminte p1 , e p2 u, Q pp1,φ1q pp2,φ2q ss, we have mintT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2qu ¥ ∆ p1 pqq¡∆ p2 pqq and maxtT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2qu ¡ ∆ p1 pqq¡∆ p2 pqq, then we can state that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q.

where ∆ p pqq, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , αq, and U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , αq are defined as follows:

∆ p pqq 5 Φ ¡1 µpp0q pq ¡ e oppq,µpp0q q if n p ¥ 1 0 if n p 0 U pp1,φ1q
pp2,φ2q pq, q 1 , q 2 , αq 6 9 9 8 9 9 7

InvSoC q1 pp1,φ1q pmaxt0, q ¡ minte p1 , e p2 uuq ¡ InvSoC q2 pp2,φ2q pmaxt0, q ¡ e p2 uq if n p1 ¥ 1 t p1 ¡ InvSoC q2 pp2,φ2q pmaxt0, q ¡ e p2 uq if n p1 0 α 1 t p1 ¡ InvSoC q2 pp2,φ2q pq ¡ minte p1 , e p2 uq if n p1 0 α 2 U pp1,φ1q
pp2,φ2q pq, q 1 , q 2 , αq 6 9 9 8 9 9 7

InvSoC q1 pp1,φ1q pmaxt0, q ¡ e p1 uq ¡ InvSoC q2 pp2,φ2q pmaxt0, q ¡ minte p1 , e p2 uuq if n p2 ¥ 1 InvSoC q1 pp1,φ1q pq ¡ minte p1 , e p2 uq ¡ t p2 if n p2 0 α 1 InvSoC q1 pp1,φ1q pmaxt0, q ¡ e p1 uq ¡ t p2 if n p2 0 α 2 Proof. See B.
Using this proposition, we design an algorithm that enumerates all non-dominated CSPs between any pair of nodes i, j r V belonging to the multigraph. We first build a directed graph Ĝij pN ij , A ij q, where N ij ti, ju F . For each CS node l F , we create the arcs pi, lq and pl, jq. For each pair of CSs l, l I , we also create the arc pl, l I q. We also consider the arc pi, jq. We apply the preprocessing rules described in A on Ĝij . We then solve on this graph a multi-objective constrained shortest path problem from i to j where 1) the two objectives are the minimization of the total duration of the path (driving + charging time) and the maximization of the final SoC; 2) there is a constrained resource: the energy (i.e., the SoC of the EV must always be in the interval r0, Qs); and 3) the SoC upon leaving i and the SoC needed upon reaching j are unknown. To this end, we use a label-correcting algorithm where each label represents a RP. A SoC function is associated with each RP pp, φq assuming that the initial SoC is equal to q M IN pp,φq . Non-dominated CSPs are then associated with Pareto optimal RPs. We outline the pseudo-code of our label-correcting algorithm in Algorithm 1.

We initialize the algorithm with one label at node j representing a RP associated with arc (i,j). We also initialize the algorithm with one label for each l F such that pi, lq A ij . This label represents a partial RP associated with partial CSP (i,l). At each iteration of the algorithm, we select the label (among the unprocessed labels) with the minimum duration and we extend it along the arcs of Ĝij . Selecting this specific label is intended to reduce the number of labels the algorithm generates. If we extend a label along an arc pl, l I q, then the new label consists of the translation of the SoC function along the vector (t ll I, ¡e ll I). Note that some supporting points of the new SoC function (after the extension along arc pl, l I q) can have a negative Algorithm 1: ComputeCSPs(i,j) input : Ĝij output: a set containing all the non-dominated CSPs between i and j 1 forall l N ij ztiu do L 1 plq Ð r (L 1 plq stores the unprocessed labels associated with node l) ; 2 forall l N ij ztiu do L 2 plq Ð r (L 2 plq stores the processed labels associated with node l) ; 3 Add to L 2 pjq the RP associated with arc pi, jq 4 R Ð r (R stores nodes from N ij ztiu with unprocessed labels) 5 forall l F such that pi, lq A ij do 6 Add to L 1 plq the label corresponding to the RP starting at i and going directly to l

7 end 8 R Ð tl F : pi, lq A ij u 9 while R $ r do l Ð minElement(R) (gets from
R the node which has associated the unprocessed label with the minimum total time) pp l , φ l q Ð deleteMin(L 1 plq) (removes from L 1 plq the label with the minimum total time)

if L 1 plq r then R Ð Rztlu ; Update(L 2 plq, pp l , φ l q)
forall pl, l I q Ĝij do if IsExtendable(pp l , φ l q,l I ) then (checks if RP pp, φq can be extended to node l I ) pp l I, φ l Iq ÐExtendppp l , φ l q, l I q (extends RP pp l , φ l q to node l I ) K ÐCreateLabelsppp l I, φ l Iqq (if l I is a CS and p l I contains more than a CS, then it creates the relevant labels according to the supporting points of the SoC function associated with RP pp l I, φ l Iq, otherwise K only contains pp l I, φ l Iq)

forall pp, φq K do Update(L 1 pl I q, pp, φq) (inserts RP pp, φq in L 1 pl I q if it is not dominated according to Proposition 4.1, and verifies if previously inserted labels become dominated)

end if L 1 pl I q has changed and l I R then R Ð R tl I u ; end end end P Ð r
forall pp j , φ j q L 2 pjq do if p j P then P Ð P tp j u ; end return P y-coordinate, which implies violating a resource constraint. Therefore, we do not consider these points. If all the supporting points of the new SoC function have a negative y-coordinate, then traversing the arc is impossible. If l I is a CS, we need to set the SoC target value at the previous CS (if there is any in the current path). We note that setting this SoC target value requires knowing the successive CS (if it exists). Stopping charging at the previous CS encountered in the path and starting charging at the current CS is only relevant at a supporting point of the current SoC function [START_REF] Zündorf | Electric vehicle routing with realistic recharging models[END_REF]. Hence we create one new label for each supporting point of the current SoC function in order to explore the possibility of switching over to the new CS at that point: Note that we can restrict the creation of new labels to the case where the slope of the charging function at the new CS is larger than that of the current SoC function. For each of these labels, the new SoC function is derived from the charging function associated with the new CS. Moreover, whenever we obtain a new label at a node, we check whether the corresponding partial RP is non-dominated before accepting it and we verify if previously generated labels associated with the same node become dominated. Specifically, using the result of Proposition 4.1 we can identify intervals for the SoC at the departure from i for which the newly generated label is dominated by or dominates a processed or an unprocessed label associated with the same node. If a RP pp, φq becomes dominated for every SoC q rq M IN pp,φq , q M AX pp,φq s, then it can be safely removed from the set of labels. The algorithm terminates when there are no more labels to set at the destination node.

To speed up the algorithm, we do not extend a label pp, φq associated with a node l to a node l I if at least one of the following conditions is met: 1) the sum of the minimum duration associated with the SoC function at the current label (i.e., the x-coordinate of the first supporting point), the driving time t ll I, and the quantity t 0i t l I j t j0 is larger than the sum of the duration limit T max and the quantity Φ ¡1 µpp0q pQq; 2) If l I is a CS node, then while scanning the nodes of the route backwards we encounter the same node l I without having scanned before a CS with a faster technology (note that a CS can be visited more than once in a path). The latter condition is due to the positivity of the energy consumption and driving time.

The fixed route vehicle charging problem

Building on the insights obtained by the previous sections, we now focus on the fixed route vehicle charging problem (FRVCP) which is a subproblem of the E-VRP-NL. We recall that FRVCPs are fundamental subproblems of E-VRP heuristics. Let Π be a route that does not visit any CS. The objective of the FRVCP is to determine the charging operations (where to insert the CSs and how much to charge), in order to minimize the total route duration while satisfying the following conditions: the customers in the resulting route are visited in the initial order, and the resulting route is energy-feasible and satisfies the maximum-duration limit. We restrict our attention to routes that cannot be traveled by an EV without the need to recharge the battery. We therefore need to insert CSs in the route to make it energy-feasible. Figure 3 shows an example of the FRVCP with a route serving four customers and the possibility of visiting two CSs (as many times as necessary). In Here we propose three alternative algorithms to solve the FRVCP. Let Π p πp0q, πp1q, ¤ ¤ ¤ , πpiq, ¤ ¤ ¤ , πpnpΠq ¡ 1q, πpnpΠqq q be a fixed route, where πp0q and πpnpΠqq represent the depot and the remaining πpiq are the visited customers. We first describe in §5.1 a preprocessing procedure to reduce the size of the search space. We then present in §5.2 a MILP formulation for the problem. Contrary to the formulation by [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF], our model accounts for visits to multiple CSs between any pair of non-CS nodes. Finally, we present in §5.3 an exact labeling algorithm for this formulation. [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] also introduced a greedy heuristic for the FRVCP. Contrary to their formulation, their heuristic allows visiting a number of CSs between two successive nodes of Π. Building on their work, we develop in §5.4 a new FRVCP heuristic.

Preprocessing procedure

Our preprocessing procedure works as follows. We first compute, according to Π and T max , the maximum time that is possible to spend detouring and charging at CSs. Using this value, we identify a priori the set of feasible CS insertions and the set of feasible consecutive CS visits between every pair of consecutive nodes of Π. We also remove additional infeasible connections between the nodes of Π and the CSs. More specifically, we perform the following preprocessing steps: 1) Let ∆tpΠq T max ¡ °npΠq¡1 i0 t πpiq,πpi 1q be the maximum time that is possible to spend detouring and charging at CSs; 2) Let πpiq be the i th node in Π and j a CS. For i t0, . . . , npΠq ¡ 1u, if t πpiq,j t j,πpi 1q ¡ t πpiq,πpi 1q ¡ ∆tpΠq, then we forbid the insertion of CS j between nodes πpiq and πpi 1q; 3) Let j, j I be two different CSs. If t πpiq,j t j,j I t j I ,πpi 1q ¡ t πpiq,πpi 1q ¡ ∆tpΠq, then we forbid visiting CS j I after j between nodes πpiq and πpi 1q; 4) Let i, j be two nodes that belong to Π or to F . If min lF t0u te li u e ij min lF t0u te jl u ¡ Q, then we forbid visiting node j directly after node i (since it violates the battery constraints).

A path-based model for the FRVCP

Building on the path-based formulation proposed in Section 4, we now propose a MILP formulation for the FRVCP. Let P pΠq be the set of CSPs between two consecutive nodes of route Π. To generate these paths, we use the labeling algorithm described in Algorithm 1. The construction of the underlying graph takes into account the results of the preprocessing procedure. We use the same decisions variables as in formulation rF P ath Arc s but we define them only for paths belonging to P pΠq. The binary variable x p is equal to 1 if and only if an EV travels on CSP p P pΠq. For a path p, continuous variable y p track the time at which the EV leaves oppq. The variables c pl and c pl are the scaled start and end times for charging an EV according to the charging function of CS µ p plq. Continuous variables q pl and q pl specify (according to the piecewise linear approximation of the charging function) the SoC of an EV when it enters and leaves µ p plq, i.e., the CS at position l L p . The continuous variable ∆ pl represents the duration of the charging operation performed at µ p plq. The binary variables w plk and w plk are equal to 1 if and only if the SoC is between a µpplq,k¡1 and a µpplq,k , with k B i zt0u, when the EV enters and leaves CS µ p plq, respectively. The continuous variables λ plk and λ plk are the coefficients associated with the breakpoint pc µpplq,k , a µpplq,k q in the piecewise linear approximation, when the EV enters and leaves CS µ p plq, respectively. Let e p and t p be the energy consumption and the driving time associated with path p P . A formulation of the FRVCP for a route Π, denoted as rFR P ath pΠqs, is as follows:

rFR P ath pΠqs minimize pP pΠq ¤ ¥ t p x p ļ Lp ∆ pl (86) subject to pP πpiq,πpi 1q x p 1, i t0, ¤ ¤ ¤ , npΠq ¡ 1u (87) pP πpi¡1q,πpiq ¤ ¥ y p ¡ e p x p ļ Lp pq pl ¡ q pl q Pπpiq,πpi 1q y p , i t1, ¤ ¤ ¤ , npΠq ¡ 1u (88) y p ¡ e oppq,µpp0q x p q p0 , p P pΠq : |L p | $ 0 (89) q p,l¡1 ¡ e µppl¡1q,µpplq x p q pl , p P pΠq, l L p zt0u (90) pP πpnpΠq¡1q,πpnpΠqq ¤ ¥ y p ¡ e p x p ¡ ļ Lp pq pl ¡ q pl q ¥ 0, (
λ plk ¤ w plk w pl,k 1 p P pΠq, l L p , k B µpplq zt0, b µpplq u (107) λ plb µp plq ¤ w plb µp plq p P pΠq, l L p (108) ∆ pl c pl ¡ c pl p P pΠq, l L p (109)
x p t0, 1u, p P pΠq (110)

y p ¥ 0 p P pΠq (111)
c pl , c pl , q pl , q pl , ∆ pl ¥ 0 p P pΠq, l L p (112)

w plk t0, 1u, w plk t0, 1u p P, l L p , k B i zt0u

(113)

λ plk ¥ 0, λ plk ¥ 0 p P, l L p , k B i . (114)
The objective function (86) minimizes the total driving and charging time in the route. Constraints (87) ensure the selection of a CSP between each pair of successive non-CS nodes. The remaining constraints are similar to those introduced for model rF P ath Arc s.

An exact labeling algorithm for the FRVCP

We now present an exact labeling algorithm for the FRVCP. We apply it on a directed graph GpΠq that contains the nodes of the initial route Π and the CS nodes. To enforce the order of the visits of the nodes that belong to Π, we replicate the CS nodes in such a way that one copy of every CS is associated with one arc of the initial route Π. Specifically, for every value of i between 0 and npΠq ¡ 1, we create a copy lpiq of CS l F . We then create the arc pπpiq, πpi 1qq, the arcs pπpiq, lpiqq and plpiq, πpi 1qq for each CS l F , and the arcs plpiq, l I piqq for each pair of CS l, l I F . It is worth mentioning that we only add these arcs to the graph if they still exist after running the preprocessing routines introduced in §5.1. Figure 4 depicts the graph considered for the example of Figure 3.

The FRVCP can then be defined as a constrained shortest path problem from πp0q to πpnpΠqq on GpΠq.

The objective is to minimize the path duration while the SoC acts as a resource constraint (i.e., the SoC must always be in the interval r0, Qs). To solve this problem, we propose the label-setting algorithm described in Algorithm 3, which shares many similarities with Algorithm 1. Algorithm 3 uses RPs as labels and exploits the concept of SoC functions presented in §4.2. Relaxing an arc and setting a label at a CS follows the same procedure as the one described in in §4.2. The algorithm terminates as soon as a label for node npΠq is set.

The quantity of energy charged at the last CS of this RP is computed such that the EV reaches the depot with an empty battery.

To improve the efficiency of the algorithm, we apply two different strategies. We make use of lower bounds.

Since we know the sequence of nodes to visit, we compute beforehand at each node j GpΠq lower bounds using Dijkstra's algorithm. We first compute a lower bound Ý Ñ t j on the time and a lower bound Ý Ñ e j on the energy required to finish the route by running a backward search from πpnpΠqq. If Ý Ñ e j is larger than Q, we can increase the value of Ý Ñ t j by assuming that the difference Q ¡ Ý Ñ e j will be charged using the fastest charging rate.

Specifically, we set the value of Ý Ñ t j to Ý Ñ t j maxt0, Ý Ñ e j ¡ Qu{ρ ¦ where ρ ¦ corresponds to the steepest slope for a segment of the piecewise linear approximation of the charging functions. We also compute a lower bound Ð Ý e j on the energy consumption from the departing depot to node j by running a forward search from πp0q. Then, we do not extend a label associated with a node i to a node j if one of the following four conditions is met: 1) the maximum SoC an EV can have upon departure from j is strictly less than the energy needed to reach to the nearest CS according to the energy consumption, or is less than the SoC it could have attained by traveling on the route to j without visiting any CS. Specifically, the difference between the maximum of the SoC function and the energy consumed e ij is strictly smaller than maxt min

lF 0 te jl u, Q ¡ Ð Ý e j u;
2) The sum of the minimum duration associated with the SoC function at the current label (i.e., the x-coordinate of the first supporting point), the driving time t ij , and Ý Ñ t j is larger than the duration limit T max ; 3) If j is a CS node, then while scanning the nodes of the route backward we encounter the same node j without having scanned before the depot, a customer, or a CS with a faster technology; 4) If j is a CS node, the value of the inverse of the SoC function evaluated in Ý Ñ e j is less than the sum of the minimum time associated with the current label and the time increase due to the detour to j.

Finally, to avoid enumerating all the RPs from i to j, we use the dominance criterion introduced in Definition 4.1 (for a SoC upon departing to the depot equal to Q) to discard dominated RPs throughout the algorithm.

A heuristic algorithm for the FRVCP

The goal of this section is to design an effective and efficient heuristic algorithm, which could be periodically called inside a local search algorithm. Such an algorithm would prevent the search from ignoring regions of the search space that are feasible with respect to route duration constraints. The algorithm borrows some ideas from the FRVCP heuristic of [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF]. We describe the general framework of our heuristic in Algorithm 2. We refer the reader to E.2 for the algorithmic details on the different procedures of the heuristic.

We note that this heuristic inserts CSs into Π.

Throughout the algorithm, we compute the SoC Y pΠ, iq of the EV upon arrival at node i Π assuming the following heuristic (H) charging policy: at each CS, the EV is charged in such a way that its departure SoC is equal to the minimum between the battery capacity Q and the energy required to reach the next depot or CS in Π. We then define the energy deficit at node i as maxt0, ¡Y pΠ, iqu.

The algorithm starts with procedure locateCS1(¤) that scans the nodes in the route (see Algorithm 4). If the arc connecting two successive nodes does not exist (this rarely happens but it requires special attention), then the algorithm tries to insert a CS between these nodes by calling procedure insertCS(¤) (see Algorithm 5). Inserting a CS between two nodes follows the same scheme throughout the algorithm. We evaluate the insertion of each CS by calling function getLBTime(¤) which computes a lower bound on the duration to travel the resulting route. The lower bound value is the sum of two components: the driving time and a lower bound on the time to charge the energy deficit observed upon arrival at CSs and at the return to the depot. This latter bound is computed by dividing the energy deficit by the fastest charging rate of the instance. The energy deficit, computed using function getEnergyDeficit(¤), is equal to °iΠpF t0uq maxt0, ¡Y pΠ, iqu. We retain the insertion that leads to the minimum lower bound value. We forbid any insertion if after executing it the energy deficit computed at the inserted CS and at the next CS that is part of the route, if any, or the depot is larger than 0. Note that it may happen that the algorithm selects an insertion leading to an energy infeasible route rather than an insertion leading to an energy-feasible route. Throughout the execution of the algorithm, we therefore store the energy-feasible route with the minimum duration, assuming the H charging policy, returned by function getTime(¤).

Whenever there exists an arc between each pair of successive nodes, we try to insert additional CSs by calling procedure locateCS2(¤) (see Algorithm 6). Here, we iteratively test the insertion of a CS between two successive nodes in the current route and we proceed to the insertion yielding the minimum lower bound on the duration to travel the route. We try to insert additional CSs until we cannot improve the duration of the current route. Specifically, this procedure terminates when the lower bound on the duration of the current route is larger than or equal to the duration of the best energy-feasible route.

We then consider the best energy-feasible route computed during the previous steps of the algorithm.

Calling procedure defineChargingAmount(¤) (see Algorithm 7), we check whether it is possible to transfer some charge between the CSs of the route. Indeed, it may be profitable in some cases to charge at a CS more than the energy required to reach the next CS. To this end, we scan the CSs starting by the end of the route, and we try to transfer some energy from a CS to the CS that precedes it. A transfer of energy reduces the duration of the route if the slope of the charging function associated with the departure SoC of the EV at the previous CS is larger than the slope of the charging function associated with the arrival SoC of the EV at the current CS. Finally, we check whether the maximum-duration limit is exceeded by the optimized route.

Algorithm 2: ProcedureFRVCP(Π 0 ) input : A route Π 0 tπp0q, πp1q, ¤ ¤ ¤ , πpiq, ¤ ¤ ¤ , πpjq, ¤ ¤ ¤ , πpnpΠq ¡ 1q, πpnpΠqqu output: A couple (f, Π) where the route Π is an energy-feasible route whose duration does not exceed the duration limit

Tmax if boolean f is equal to true, or the initial route Π 0 if boolean f is equal to false 1 Apply preprocessing rules on the multigraph induced by Π and 6 Computational experiments

F 2 pf, Π, Π best q Ð locateCS1(Π 0 ) 3 if f true getLBTime(Π)¤ Tmax then 4 Π best Ð locateCS2(Π,Π best ) 5 if Π best $ null then 6 (f, Π best ) Ð reviseChargingAmount(Π best
We used Gurobi 7.5.0 (through its Java API) to solve the MILP models. All tests were performed on a single thread with 12 GB and on a cluster of 27 computers, each having 12 cores and two Intel(R) Xeon X5675 3.07

GHz processors. In all the tables, the CPU time is rounded to the nearest integer.

The first aim of our computational experiments is to assess and compare the performance of the CS replication-based and the path-based models, which are presented in Sections 3 and 4, for solving small size instances to the E-VRP-NL. These results are presented in §6.1. The second aim of our computational experiments is to assess the quality of the algorithms developed in Section 5 for the FRVCP. These results are presented in §6.2.

Results for E-VRP-NL formulations

For the E-VRP-NL we considered the twenty 10-customer and the twenty 20-customer instances of the Montoya et al. ( 2017) testbed (publicly available at http://vrp-rep.org). In these instances, both the energy consumption and travel times are obtained by multiplying the arc length by constant factors (consumption per unit distance and speed). It is worth noting that none of our results relies on the assumption that the energy consumption linearly depends on the travel time (this is only a characteristic of the instances). We ran the MILP models with a three-hour time limit.

CS replication-based models

We first tested the CS replication-based models rF CS rep N ode s and rF CS rep Arc s presented in Section 3. The number of copies of each CS i F was set to an integer value β ¥ 1 (i.e., β i β, di F ).

Table 1 reports for each model and each value of β the number of instances proven infeasible by the solver within the time limit (#Inf), the number of instances reaching the time limit without finding any feasible solution (#Unk), the number of instances with a feasible solution not proven to be optimal (#Opt), the number of instances with a feasible solution proven to be optimal (#Opt), the average CPU time in seconds (Time (s)) for the instances in (#Opt), and the average gap (Gap) for the instances in (#Opt). We compute the gap as pz ¡ z LB q{z, where z is the objective function value of the best integer solution returned by the solver, and z LB is the best lower bound retrieved by the solver running the corresponding model. The detailed results for all the tested instances are reported in F.

We first observe that for small values of β the rF CS rep Arc s model cannot always prove the infeasibility of some instances. This may be due to the relatively large number of variables it contains (around twice as many as the rF CS rep N ode s model for the tested instances). Second, the results confirm that the value of β influences the feasibility of the instances. Since we do not know any procedure for fixing β, while guaranteeing optimality, we can only compare models for identical values of β. The best results are obtained with arc-based tracking constraints for which the number of optimal solutions is larger and the solution time is significantly reduced, compared to node-based tracking. Moreover, when the time limit is reached, the best solution returned by the arc-based models has a better quality. Finally, our results illustrate the difficulty of optimally solving even small instances with the CS replication-based models.

Path-based model

We also tested the path-based model rF P ath Arc s presented in Section 4. We first assessed the impact of the dominance rule described in §4.2. Applying this rule reduces the average number of paths by a factor of two for the 10-customer instances (on average 647 paths were generated without the dominance rule versus 270 with it), and by a factor of five for the 20-customer instances (on average 6512 paths were generated without the dominance rule versus 1314 with it). This leads to more tractable instances: the average number of paths between each couple of non-CS nodes is approximately 2.4 and 3.2 for the 10-customer and 20-customer instances, respectively. Our results show that the MILP solver can optimally solve all the 10-customer instances within an average time of seven minutes. The solver fails to solve most of the instances with 20 customers due to the large size of the models. Doubling the number of customers from 10 to 20 multiplies the number of paths by five (from 270 to 1314, on average).

Figure 5 compares the CS replication-based and path-based models. Specifically, the figure shows the number of instances optimally solved according to the solution time for the different models. We conclude that the MILP solver performs better on the path-based model than on the CS replication-based models. Moreover, we recall that using the path-based model ensures that optimal solutions are not cut off.

Results for FRVCP algorithms

Since the FRVCP problem is typically solved within local search algorithms, we wanted to compare the speed of the algorithms we developed and the quality of the routes they provide. We tested the three algorithms described in Section 5 to solve the FRVCP: our formulation rFR P ath pΠqs with the MILP solver, our labeling algorithm, and our heuristic (referred to as MILP PATH, LABEL, and HEURISTIC). As a basis of comparison, Moreover, in order to assess the relevance of considering the insertion of multiple CSs between each pair of non-CS nodes of the initial route, we also solved the FRVCP-1. To this end, we adapted our path-based MILP formulation by only considering paths with at most one CS. We also adapted the graph on which we applied our labeling algorithm. These modifications yield a new MILP model and a new labeling algorithm (referred to as MILP PATH 1 and LABEL 1) for this special version of the FRVCP. We also implemented the FRVCP formulation of [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] (referred to as MILP M2017). No computational time limit was imposed on these tests.

We first tested the algorithms on the pool of routes used by [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] in the second phase of their metaheuristic. This pool consists of all routes encountered at each local optima solution found by an iterated local search procedure. The pool contains a total of 29,443 routes having between 3 and 15 nodes, which may visit between 2 and 24 CSs. For each algorithm, in Table 3 we report the number of routes for which the duration is below the duration limit T max (#Feas), the number of routes solved to optimality (#Opt), the average CPU time in milliseconds (Time (ms)), the average gap (Gap), and the maximum observed gap over all routes (Gap max).

Considering the detailed results presented in F.2 and F.3, we conclude that for all algorithms the number of CSs and the number of nodes have a relatively small impact on the solution time and quality. We note that the two exact algorithms for the FRVCP, i.e., MILP PATH and LABEL, were able to solve all 29,443 routes to optimality within an average runtime of less than 100 milliseconds. LABEL is more than 90 times faster than MILP PATH. HEURISTIC produces competitive results with a small average gap in a very short time. It also outperforms HEURISTIC M2017 in terms of gap. However, both heuristics may fail to find existing feasible solutions. Such failures are due to the fact that feasibility often coincides with optimality. Since our heuristic found many more optimal solutions, compared to the FRVCP heuristic of [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF], it found more feasible solutions. In particular, HEURISTIC failed to produce a feasible solution in 2.1% of the routes, whereas the HEURISTIC M2017 failed to produce a feasible solution in 13.7% of the routes. Although LABEL outperforms HEURISTIC, the latter is easier to implement than the former which can make it attractive for some applications.

We observe that MILP PATH 1 solves FRVCP-1 much faster than MILP M2017. However, LABEL 1 substantially outperforms both of these formulations. Comparing FRVCP with FRVCP-1, we notice that restricting the insertion of at most one CS between any pair of non-CS nodes does not impact the quality of the solutions in most cases. Specifically, for 28,085 out of the 29,443 routes, the optimal objective function value considering either the FRVCP or the FRVCP-1 are the same. However, for the remaining 1358 routes, the charging decisions obtained by solving the FRVCP-1 are not optimal and in the worst case they can be quite far from the optimal ones. Specifically, let z ¦ and z ¦ 1 correspond to the average optimal objective function value obtained when solving the FRVCP and the FRVCP-1 on these 1358 routes. We observe that pz ¦ 1 ¡z ¦ q{z ¦ 1 is approximately 0.01%. However, in the worst case (i.e., for a particular route), the deviation between the two optimal objective function values is 7.68%. We also tested our labeling algorithm of §5.3 on the 1426 routes making up the best solutions reported by [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] for 120 instances with 10, 20, 40, 80, 160, or 320 customers (20 instances for each instance size). We recall that the charging decisions for these routes were found using the FRVCP formulation of [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF], which assumes at most one CS visit between two successive non-CS nodes. Therefore, we want to assess if the results can be improved by solving the more general version of the FRVCP, which allows multiple consecutive visits to CSs. We report in Table 4 the name of the instances for which a better solution has been identified, the number of routes (#R ¦ ) for which the charging decisions have been revised

over the number of routes (#R) in the initial solution. We show the maximum gap improvement between a route that has been revised and the initial route (Max Gap R ¦ ). We also report the average gap with respect to the previous best known solution (Gap). We found 23 new best E-VRP-NL solutions. On these instances, taking optimal charging decisions improved the solutions of Montoya et al. ( 2017) by approximately 0.09%.

This shows that charging decisions are critical in the E-VRP-NL. Therefore, there is ground to believe that the use of the labeling algorithm inside a metaheuristic algorithm deserves further study. The results presented in this section open an interesting research avenue. The only way to obtain an exact evaluation of a move in local search algorithms for the E-VRP-NL is to re-optimize the charging decisions (i.e., where and how much to charge) for the route(s) involved in the move. Specifically, one must decide the charging decisions such that each route has the lowest possible cost. This is a daunting task because it implies solving one (inter route moves) or more (intra route moves) FRVCPs. As [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] pointed out, exactly evaluating the moves using MILPs is intractable (our findings only confirm their observation). For this reason, most existing approaches use fast-to-compute proxy move evaluations. As a matter of fact, only the most sophisticated approaches use MILPs or exact algorithms at certain points (e.g., every time they reach a local optimum) to improve the charging decisions. A notable exception is the multi-start local search heuristic proposed by [START_REF] Andelmin | A multi-start local search heuristic for the green vehicle routing problem based on a multigraph reformulation[END_REF] for the Green VRP. This heuristic works on a multigraph (where each arc represents a refuel path) avoiding the need to explicitly consider charging stations. Our results suggest that because of its good computational performance, our labeling algorithm may be used to exactly evaluate a large number of moves (if not all) in local search algorithms without penalizing the overall performance of the method.

Conclusions

We Most of the available E-VRP models, as well as the arc-based formulation, are based on pre-setting the number of CS copies. We have demonstrated that generating a limited number of copies compromises the quality of the solution by possibly eliminating optimal solutions or by yielding an infeasible model. With more copies, the models become larger and more difficult to solve, even for small-size instances. To overcome these drawbacks, we have proposed a third formulation for the E-VRP-NL based on the concept of CSPs between non-CS nodes.

Using dominance rules to discard unpromising paths, we have shown through our computational experiments that this third formulation systematically yields better results than the first two formulations. All 10-customer instances could be solved to optimality within a much shorter time. As expected, solving our models using a commercial solver allows us to solve only small-size instances. However, most E-VRPs have been handled

by node-based formulations, and for these cases the arc-based model provides a highly relevant alternative formulation. Furthermore, the introduction of the path-based model provides a new promising avenue for the development of exact and heuristic algorithms for the E-VRP-NL.

Building on the formulations we have developed for the E-VRP-NL, we have also introduced heuristic and exact algorithms to solve the fixed route vehicle charging problem (FRVCP). Given an initial energy infeasible route, this problem consists in finding the charging decisions that minimize the total duration of the route. We have proposed a new MILP formulation, an exact labeling algorithm, and a heuristic algorithm. Computational experiments have shown that taking optimal or near-optimal charging decisions is crucial in producing high quality E-VRP-NL solutions. Indeed, we have improved 23 out of 120 best known E-VRP-NL solutions by allowing the visit of multiple CSs between two non-CS nodes. Since the labeling algorithm is able to optimally solve the FRVCP within a very short time, it can be easily embedded in future E-VRP-NL heuristics.

A Graph preprocessing

Some arcs from the arc set A can be removed without cutting off the optimal solution. Specifically, we introduce the following preprocessing steps to reduce the size of the arc set:

1. Considering i, j V , if min lF t0u te li u e ij min lF t0u te il u ¡ Q then we infer that the maximum possible SoC at i the EV is not sufficient for traversing arc pi, jq and reaching the nearest CS or the depot. Therefore, such an arc is removed from A.

2. Let t ij be the duration of a shortest path from i to j in G (after applying preprocessing rule 1) with respect to arc duration t ij . Let ec min ij be the minimum charging amount necessary to build a route traversing arc pi, jq. Then ec min ij maxt0, pe 0i e ij e j0 q ¡ Qu. We also define t charge ij as a lower bound on the charging time necessary to perform a route traversing the arc pi, jq. Namely,

t
D G G G F G G G E D G G G F G G G E
The value of ρ ¦ corresponds to the steepest slope for a segment of the piecewise linear approximation of the charging functions (i.e., ρ ¦ max lF tρ l1 u). ¡ T max , we remove the arc pi, jq from set A, since traversing the arc pi, jq leads to a violation of the route maximum-duration limit.

B Proofs

Proof of Proposition 3.1. Assume that the arc-based tracking constraints are satisfied by setting every variable x ij , y ij , q j , q j , τ ij , and ∆ j to x ij , y ij , q j , q j , τ ij , and ∆ j , respectively. We now prove that the node-based tracking constraints are satisfied if: 1) the variables x ij are set to values x ij for every pi, jq A, 2) the variables q j , q j , ∆ j are set to values q j , q j , ∆ j for every j F , and 3) the variables y i , τ i are set °pi,lqA y il and °pi,lqA τ il for every i V zt0u. We also set the value of τ 0 to 0 and the value of y 0 to Q.

• Let pi, jq A such that j I. Then,

-if i V zt0u, y i ¡ y j pi,lqA y il ¡ pj,lqA y jl y i ¡ y j pi,lqA y il ¡ pl,jqA
py lj ¡ e lj x lj q from constraints (39)

y i ¡ y j e ij x ij pi,lqA:l$j y il ¡ pl,jqA:l$i py lj ¡ e lj x lj q y i ¡ y j ¤ e ij x ij pi,lqA:l$j y il y i ¡ y j ¤ e ij x ij
pi,lqA:l$j pQ ¡ min l I Ft0u te l I i uqx il from constraints ( 42)

y i ¡ y j ¤ e ij x ij Q pi,lqA:l$j x il y i ¡ y j ¤ e ij x ij Qp1 ¡ x ij q from constraints (2)
and

y i ¡ y j pi,lqA y il ¡ pj,lqA y jl y i ¡ y j pi,lqA y il ¡ pl,jqA
py lj ¡ e lj x lj q from constraints (39)

y i ¡ y j e ij x ij ¡ pl,jqA:l$i y lj pi,lqA:l$j y il pl,jqA:l$i e lj x lj y i ¡ y j ¥ e ij x ij ¡ pl,jqA:l$i y lj y i ¡ y j ¥ e ij x ij ¡ pl,jqA:l$i pQ ¡ min l I Ft0u
te l I l uqx lj from constraints ( 42)

y i ¡ y j ¥ e ij x ij ¡ Q pl,jqA:l$i x lj y i ¡ y j ¥ e ij x ij ¡ Qp1 ¡ x ij q.
from constraints (2)

-if i 0, y 0 ¡ y j Q ¡ pj,lqA y jl y 0 ¡ y j Q ¡ pl,jqA
py lj ¡ e lj x lj q from constraints (39)

y 0 ¡ y j Q e 0j x 0j ¡ y 0j ¡ pl,jqA:l$0 py lj ¡ e lj x lj q y 0 ¡ y j ¤ e 0j x 0j Q ¡ Qx 0j
from constraints (38)

y 0 ¡ y j ¤ e 0j x 0j Qp1 ¡ x 0j q from constraints (2)
and

y 0 ¡ y j Q ¡ pj,lqA y jl y 0 ¡ y j Q ¡ pl,jqA
py lj ¡ e lj x lj q from constraints (39)

y 0 ¡ y j Q e 0j x 0j ¡ pl,jqA y lj pl,jqA:l$0 e lj x lj y 0 ¡ y j ¥ Q e 0j x 0j ¡ Q pl,jqA
x lj from constraints ( 42)

y 0 ¡ y j ¥ Q e 0j x 0j ¡ Q from constraints (2)
y 0 ¡ y j ¥ e 0j x 0j y 0 ¡ y j ¥ e 0j x 0j ¡ Qp1 ¡ x 0j q.
Constraints ( 5) are therefore satisfied.

• Let pi, jq A such that j F I . Then,

-if i V zt0u, y i ¡ q j pi,lqA y il ¡ pl,jqA
py lj ¡ e lj x lj q from constraints ( 40)

y i ¡ q j e ij x ij
pi,lqA:l$j y il ¡ pl,jqA:l$i py lj ¡ e lj x lj q

y i ¡ q j ¤ e ij x ij
pi,lqA:l$j y il

y i ¡ q j ¤ e ij x ij Q pi,lqA:l$j
x il from constraints (42)

y i ¡ q j ¤ e ij x ij Qp1 ¡ x ij q from constraints (3)
and

y i ¡ q j pi,lqA y il ¡ pl,jqA
py lj ¡ e lj x lj q from constraints (40)

y i ¡ q j e ij x ij ¡ pl,jqA:l$i y lj pi,lqA:l$j y il pl,jqA:l$i e lj x lj y i ¡ q j ¥ e ij x ij ¡ pl,jqA:l$i y lj y i ¡ q j ¥ e ij x ij ¡ Q pl,jqA:l$i
x lj from constraints (42)

y i ¡ q j ¥ e ij x ij ¡ Qp1 ¡ x ij q.
from constraints (3)

-if i 0, y 0 ¡ q j Q ¡ pl,jqA
py lj ¡ e lj x lj q from constraints (40)

y 0 ¡ q j Q e 0j x 0j ¡ y 0j ¡ pl,jqA:l$0 py lj ¡ e lj x lj q y 0 ¡ q j ¤ e 0j x 0j Q ¡ Qx 0j
from constraints (38)

y 0 ¡ q j ¤ e 0j x 0j Qp1 ¡ x 0j q from constraints (3)
and

y 0 ¡ q j Q ¡ pl,jqA
py lj ¡ e lj x lj q from constraints (40)

y 0 ¡ q j Q e 0j x 0j ¡ pl,jqA y lj pl,jqA:l$0 e lj x lj y 0 ¡ q j ¥ Q e 0j x 0j ¡ Q pl,jqA
x lj from constraints ( 42)

y 0 ¡ q j ¥ Q e 0j x 0j ¡ Q from constraints (3) y 0 ¡ q j ¥ e 0j x 0j y 0 ¡ q j ¥ e 0j x 0j ¡ Qp1 ¡ x 0j q.
Constraints ( 6) are therefore satisfied.

• Let i F I , then y i pi,lqA

y il y i q i .
from constraints ( 41)

Constraints ( 8) are therefore satisfied.

• Let pi, 0q A, then y i pi,lqA

y il y i ¥ y i0 y i ¥ e i0 x i0 .
from constraints (49) Constraints ( 7) are therefore satisfied.

• Let pi, jq A such that j I. Then,

-if i V zt0u, τ j ¡ τ i pj,lqA τ jl ¡ pi,lqA τ il τ j ¡ τ i pl,jqA
pτ lj pt lj g j q x lj q ¡ pi,lqA:

τ il from constraints (43) τ j ¡ τ i pt ij g j q x ij
pl,jqA:l$i pτ lj pt lj g j q x lj q ¡ pi,lqA:l$j

τ il τ j ¡ τ i ¥ pt ij g j q x ij ¡ pi,lqA:l$j τ il τ j ¡ τ i ¥ pt ij g j q x ij ¡ T max pi,lqA:l$j
x il from constraints (45)

τ j ¡ τ i ¥ pt ij g j q x ij ¡ T max p1 ¡ x ij q.
from constraints ( 2)

-if i 0, τ j ¡ τ 0 pj,lqA τ jl τ j ¡ τ 0 pl,jqA
pτ lj pt lj g j q x lj q from constraints (43)

τ j ¡ τ 0 pt 0j g j q x 0j τ 0j
pl,jqA:l$0 pτ lj pt lj g j q x lj q τ j ¡ τ 0 ¥ pt 0j g j q x 0j τ j ¡ τ 0 ¥ pt 0j g j q x 0j ¡ T max p1 ¡ x 0j q.

from constraints (2)

Constraints ( 26) are therefore satisfied.

• Let pi, jq A such that j F I . Then,

-if i V zt0u, τ j ¡ τ i pj,lqA τ jl ¡ pi,lqA τ il τ j ¡ τ i pl,jqA pτ lj t lj x lj q ∆ j ¡ pi,lqA τ il from constraints (44) τ j ¡ τ i t ij x ij ∆ j pl,jqA:l$i pτ lj t lj x lj q ¡ pi,lqA:l$j τ il τ j ¡ τ i ¥ t ij x ij ∆ j ¡ pi,lqA:l$j τ il τ j ¡ τ i ¥ t ij x ij ∆ j ¡ T max
pi,lqA:l$j

x il from constraints (45)

τ j ¡ τ i ¥ t ij x ij ∆ j ¡ T max p1 ¡ x ij q
from constraints (2)

τ j ¡ τ i ¥ t ij x ij ∆ j ¡ pT max S max qp1 ¡ x ij q.
if i 0,

τ j ¡ τ 0 pj,lqA τ jl τ j ¡ τ 0 pl,jqA
pτ lj t lj x lj q ∆ j from constraints (44)

τ j ¡ τ 0 t 0j x 0j ∆ j τ 0j pl,jqA:l$0 pτ lj t lj x lj q τ j ¡ τ 0 ¥ t ij x 0j ∆ j τ j ¡ τ 0 ¥ t ij x 0j ∆ j ¡ pT max S max qp1 ¡ x 0j q.
Constraints ( 27) are therefore satisfied.

• Let j V zt0u. Then, if j I, τ j t j0 pj,lqA

τ jl t j0 τ j t j0 ¤ pj,lqA pT max ¡ t jl ¡ p l ¡ t l0 q x jl t j0 τ j t j0 ¤ pj,lqA pT max ¡ t jl ¡ t l0 q x jl t j0
from the triangular inequality

τ j t j0 ¤ pj,lqA pT max ¡ t j0 q x jl t j0 τ j t j0 ¤ T max ¡ t j0 t j0
from constraints (2)

τ j t j0 ¤ T max and if j F I , τ j t j0 pj,lqA τ jl t j0 τ j t j0 ¤ pj,lqA T max ¡ t jl ¡ t l0 ¡ ∆ min j ¨xjl t j0 τ j t j0 ¤ pj,lqA pT max ¡ t j0 q x jl t j0
from the triangular inequality

τ j t j0 ¤ T max ¡ t j0 t j0
from constraints (2)

τ j t j0 ¤ T max .
Constraints ( 28) are therefore satisfied. By definition of the values τj , constraints ( 29) and ( 30) are also satisfied. The non-negativity constraints (33) are satisfied. Hence all the node-based tracking constraints are satisfied.

Proof of Lemma 4.1. Let us now consider a RP pp, φq containing at least one CS. If the SoC of the EV upon departing from the origin is equal to q I rather than q ¤ q I , then the EV arrives at the first CS µ p p0q with a SoC equal to q I ¡ e oppq,µpp0q rather than q ¡ e oppq,µpp0q . The inverse of the SoC function InvSoC q pp,φq includes the time charging the EV from q ¡ e oppq,µpp0q up to φp0q, which is larger than or equal to q I ¡ e oppq,µpp0q . Since the EV arrives with a larger SoC at µ p p0q when starting with a SoC equal to q I , the time spent by the EV to reach a given SoC is reduced by the time spent at µ p p0q charging the EV from q ¡ e oppq,µpp0q to q I ¡ e oppq,µpp0q . This time is equal to λ p qq I Φ ¡1 µpp0q pq I ¡ e oppq,µpp0q q ¡ Φ ¡1 µpp0q pq ¡ e oppq,µpp0q q. The inverse of the SoC function InvSoC q I pp,φq can therefore be built from the function InvSoC q pp,φq by subtracting λ p qq I from the time needed to reach a given SoC. For a SoC that is smaller than the minimum SoC that can be reached at the destination, the time is consistently set (according to the definition of the inverse of the SoC function) to the time needed to reach that minimum. Now consider a RP pp, φq that does not contain any CS. The inverse of the SoC function is equal to t p if the SoC q is less than or equal to the SoC of the EV upon arriving at the destination dppq. Otherwise, the EV can never reached this SoC and therefore the value of the SoC function is set to V.

Proof of Proposition 4.1.

Case 1: If n p1 0 or Θ pp1,φ1q ¥ Θ pp2,φ2q , we have to consider two cases: n p1 ¡ 0 and n p1 0 (note that n p2 ¡ 0 since we never check the dominance of a RP that does not contain any CS).

• Case 1.1: n p1 ¡ 0 (q 1 e opp1q,µp 1 p0q and q 2 e opp2q,µp 2 p0q ) First, observe that if q rQ pp1,φ1q pp2,φ2q , Q pp1,φ1q pp2,φ2q s is such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¤ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¤ ∆ p1 pqq ¡ ∆ p2 pqq ñ max qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tInvSoC q1 pp1,φ1q pqq ¡ InvSoC q2 pp2,φ2q pqqu ¤ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pqq ¡ InvSoC q2 pp2,φ2q pqq ¤ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pqq ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¤ InvSoC q2 pp2,φ2q pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq Ω pp1,φ1q pp2,φ2q pqq, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq from lemma 4.1 ñ dq rmaxt0, q ¡ minte p1 , e p2 uu, Θ q2 pp2,φ2q s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq

(115)
Second, observe that if q rQ pp1,φ1q pp2,φ2q , Q pp1,φ1q pp2,φ2q s is such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ∆ p1 pqq ¡ ∆ p2 pqq ñ min qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tInvSoC q1 pp1,φ1q pqq ¡ InvSoC q2 pp2,φ2q pqqu Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pq I q ¡ InvSoC q2 pp2,φ2q pq I q Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I Ω pp1,φ1q pp2,φ2q pqq, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q from lemma 4.1 ñ hq I rmaxt0, q ¡ minte p1 , e p2 uu, Θ q2 pp2,φ2q s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q (116)

Let q rQ pp1,φ1q pp2,φ2q , minte p1 , e p2 us be such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¤ ∆ p1 pqq ¡ ∆ p2 pqq and T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ∆ p1 pqq ¡ ∆ p2 pqq. Equation ( 115) is equivalent to dq r0, Θ pp2,φ2q s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq and Equation ( 116) is equivalent to hq I r0, Θ pp2,φ2q s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. This means that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q. Let q rminte p1 , e p2 u, Q pp1,φ1q pp2,φ2q s be such that maxtT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1qu ¤ ∆ p1 pqq ¡ ∆ p2 pqq and mintT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1qu ∆ p1 pqq¡∆ p2 pqq. We consider the two following subcases:

1. If q is such that minte p1 , e p2 u q ¤ maxte p1 , e p2 u, then -If minte p1 , e p2 u e p1 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ¤ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ InvSoC q2 pp2,φ2q p0q ¤ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¤ InvSoC q2 pp2,φ2q p0q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p1 s, InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¤ InvSoC q2 pp2,φ2q pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq from lemma 4.1

Since Equation ( 115) holds, we therefore have InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq for every q r0, Θ pp2,φ2q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p1 q Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q InvSoC q2 pp2,φ2q pq ¡ e p1 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p1 s, InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p1 s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q from lemma 4.1 (117)

Since either Equation ( 116) or Equation ( 117) holds, there exists qI r0, Θ pp2,φ2q s such that InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q. -If minte p1 , e p2 u e p2 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ¤ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p2 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¤ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p2 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¤ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p2 s, InvSoC q1 pp1,φ1q pqq ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¤ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq from lemma 4.1

Since Equation ( 115) holds, we therefore have InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq for every q r0, Θ pp2,φ2q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q
pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q p0q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q p0q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, InvSoC q1 pp1,φ1q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q from lemma 4.1 (118)

Since either Equation ( 116) or Equation ( 118) holds, there exists qI r0, Θ pp2,φ2q s such that InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q. 2. If q is such that maxte p1 , e p2 u q ¤ Q pp1,φ1q pp2,φ2q , then -If minte p1 , e p2 u e p1 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ¤ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¤ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¤ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq rq ¡ e p2 , q ¡ e p1 s, InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¤ InvSoC q2 pp2,φ2q pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq rq ¡ e p2 , q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq from lemma 4.1 ñ dq r0, q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq

Since Equation (115) holds, we therefore have InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq for every q r0, Θ pp2,φ2q s.

Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p1 q Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q InvSoC q2 pp2,φ2q pq ¡ e p1 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p1 s, InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p1 s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q from lemma 4.1

Since either Equation (116) or Equation (119) holds, there exists qI r0, Θ pp2,φ2q s such that InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q. -If minte p1 , e p2 u e p2 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ¤ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p2 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¤ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p2 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¤ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq rq ¡ e p1 , q ¡ e p2 s, InvSoC q1 pp1,φ1q pqq ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¤ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq rq ¡ e p1 , q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq from lemma 4.1 ñ dq r0, q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq

Since Equation ( 115) holds, we therefore have InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq for every q r0, Θ pp2,φ2q s.

Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q
pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, InvSoC q1 pp1,φ1q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q from lemma 4.1

Since either Equation (116) or Equation (120) holds, there exists qI r0, Θ pp2,φ2q s such that InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q.

• Case 1.2: n p1 0 (q 2 e opp2q,µp 2 p0q ) First, observe that if q rQ pp1,φ1q pp2,φ2q , Q pp1,φ1q pp2,φ2q s is such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¤ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¤ ∆ p1 pqq ¡ ∆ p2 pqq ñ max qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tt p1 ¡ InvSoC q2 pp2,φ2q pqqu ¤ ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq Ω pp1,φ1q
pp2,φ2q pqq, t p1 ¡ InvSoC q2 pp2,φ2q pqq ¤ ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq Ω pp1,φ1q pp2,φ2q pqq, t p1 ¤ InvSoC q2 pp2,φ2q pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq Ω pp1,φ1q pp2,φ2q pqq, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq from lemma 4.1 ñ dq rmaxt0, q ¡ minte p1 , e p2 uu, Θ q2 pp2,φ2q s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq

Second, observe that if q rQ pp1,φ1q pp2,φ2q , Q pp1,φ1q pp2,φ2q s is such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ∆ p1 pqq ¡ ∆ p2 pqq ñ max qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tt p1 ¡ InvSoC q2 pp2,φ2q pqqu ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I Ω pp1,φ1q pp2,φ2q pqq, t p1 ¡ InvSoC q2 pp2,φ2q pq I q ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I Ω pp1,φ1q pp2,φ2q pqq, t p1 InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I Ω pp1,φ1q pp2,φ2q pqq, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q from lemma 4.1 ñ hq I rmaxt0, q ¡ minte p1 , e p2 uu, Θ q2 pp2,φ2q s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q (122)

Let q rQ pp1,φ1q pp2,φ2q , minte p1 , e p2 us be such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¤ ∆ p1 pqq ¡ ∆ p2 pqq and T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ∆ p1 pqq ¡ ∆ p2 pqq. Equation ( 121) is equivalent to dq r0, Θ pp2,φ2q s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq and Equation ( 122) is equivalent to hq I r0, Θ pp2,φ2q s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. This means that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q. Let q rminte p1 , e p2 u, Q pp1,φ1q pp2,φ2q s be such that maxtT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1qu ¤ ∆ p1 pqq ¡ ∆ p2 pqq and mintT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1qu ∆ p1 pqq¡∆ p2 pqq. We consider the two following subcases:

1. If q is such that minte p1 , e p2 u ¤ q ¤ maxte p1 , e p2 u, then -If minte p1 , e p2 u e p1 , then

U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ¤ ∆ p1 pqq ¡ ∆ p2 pqq ñ t p1 ¡ InvSoC q2 pp2,φ2q p0q ¤ ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ t p1 ¤ InvSoC q2
pp2,φ2q p0q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p1 s, t p1 ¤ InvSoC q2 pp2,φ2q pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq from lemma 4.1

Since Equation ( 121) holds, we therefore have InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq for every q r0, Θ pp2,φ2q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq ñ t p1 ¡ InvSoC q2 pp2,φ2q pq ¡ e p1 q ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ t p1 InvSoC q2 pp2,φ2q pq ¡ e p1 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p1 s, t p1 InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p1 s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q from lemma 4.1

Since either Equation (122) or Equation ( 123) holds, there exists qI r0, Θ pp2,φ2q s such that InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q. -If minte p1 , e p2 u e p2 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ¤ ∆ p1 pqq ¡ ∆ p2 pqq ñ t p1 ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¤ ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ t p1 ¤ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq from lemma 4.1

Since Equation ( 121) holds, we therefore have InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq for every q r0, Θ pp2,φ2q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq ñ t p1 ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ t p1 InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, t p1 InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q from lemma 4.1

Since either Equation (122) or Equation ( 124) holds, there exists qI r0, Θ pp2,φ2q s such that InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q. 2. If q is such that maxte p1 , e p2 u ¤ q ¤ Q pp1,φ1q pp2,φ2q , then

-If minte p1 , e p2 u e p1 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ¤ ∆ p1 pqq ¡ ∆ p2 pqq ñ t p1 ¡ InvSoC q2
pp2,φ2q pq ¡ e p2 q ¤ ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ t p1 ¤ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq rq ¡ e p2 , q ¡ e p1 s, t p1 ¤ InvSoC q2 pp2,φ2q pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq rq ¡ e p2 , q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq from lemma 4.1 ñ dq r0, q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq

Since Equation ( 121) holds, we therefore have InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq for every q r0, Θ pp2,φ2q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq ñ t p1 ¡ InvSoC q2 pp2,φ2q pq ¡ e p1 q ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ t p1 InvSoC q2 pp2,φ2q pq ¡ e p1 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p1 s, t p1 InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p1 s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q from lemma 4.1

Since either Equation (122) or Equation (125) holds, there exists qI r0, Θ pp2,φ2q s such that InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q. -If minte p1 , e p2 u e p2 , then

U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ¤ ∆ p1 pqq ¡ ∆ p2 pqq ñ t p1 ¡ InvSoC q2
pp2,φ2q pq ¡ e p2 q ¤ ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ t p1 ¤ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq from lemma 4.1

Since Equation ( 121) holds, we therefore have InvSoC q pp1,φ1q pqq ¤ InvSoC q pp2,φ2q pqq for every q r0, Θ pp2,φ2q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 1q ∆ p1 pqq ¡ ∆ p2 pqq ñ t p1 ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ t p1 InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, t p1 InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q from lemma 4.1

Since either Equation (122) or Equation ( 126) holds, there exists qI r0, Θ pp2,φ2q s such that InvSoC q pp1,φ1q pq I q InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 1 , φ 1 q ¡ q pp 2 , φ 2 q.

Case 2: If n p2 0 or Θ pp2,φ2q ¥ Θ pp1,φ1q , we have to consider two cases: n p2 ¡ 0 and n p2 0 (note that n p1 ¡ 0 since we never check the dominance of a RP that does not contain any CS).

• Case 2.1: n p2 ¡ 0 (q 1 e opp1q,µp 1 p0q and q 2 e opp2q,µp 2 p0q ) First, observe that if q rQ pp1,φ1q pp2,φ2q , Q pp1,φ1q pp2,φ2q s is such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¥ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ô min qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tInvSoC q1 pp1,φ1q pqq ¡ InvSoC q2 pp2,φ2q pqqu ¥ ∆ p1 pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pqq ¡ InvSoC q2 pp2,φ2q pqq ¥ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pqq ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ InvSoC q2 pp2,φ2q pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq Ω pp1,φ1q pp2,φ2q pqq, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq from lemma 4.1 ñ dq rmaxt0, q ¡ minte p1 , e p2 uu, Θ q1 pp1,φ1q s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq

Second, observe that if q rQ pp1,φ1q pp2,φ2q , Q pp1,φ1q pp2,φ2q s is such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¡ ∆ p1 pqq ¡ ∆ p2 pqq ô max qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tInvSoC q1 pp1,φ1q pqq ¡ InvSoC q2 pp2,φ2q pqqu ¡ ∆ p1 pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pq I q ¡ InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I Ω pp1,φ1q pp2,φ2q pqq, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q from lemma 4.1 ñ hq I rmaxt0, q ¡ minte p1 , e p2 uu, Θ q1 pp1,φ1q s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q (128)

Let q rQ pp1,φ1q pp2,φ2q , minte p1 , e p2 us be such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¥ ∆ p1 pqq ¡ ∆ p2 pqq and T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¡ ∆ p1 pqq ¡ ∆ p2 pqq. Equation ( 127) is equivalent to dq r0, Θ pp1,φ1q s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq and Equation ( 128) is equivalent to hq I r0, Θ pp1,φ1q s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q. This means that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q. Let q rminte p1 , e p2 u, Q pp1,φ1q pp2,φ2q s be such that mintT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2qu ¥ ∆ p1 pqq¡∆ p2 pqq and maxtT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2qu ¡ ∆ p1 pqq ¡ ∆ p2 pqq. We consider the two following subcases:

1. If q is such that minte p1 , e p2 u ¤ q ¤ maxte p1 , e p2 u, then -If minte p1 , e p2 u e p1 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p1 q ¥ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ InvSoC q2 pp2,φ2q pq ¡ e p1 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p1 s, InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ InvSoC q2 pp2,φ2q pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq from lemma 4.1

Since Equation ( 127) holds, we therefore have InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq for every q r0, Θ pp1,φ1q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q
pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ InvSoC q2 pp2,φ2q p0q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ InvSoC q2 pp2,φ2q p0q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p1 s, InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p1 s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q from lemma 4.1

Since either Equation ( 128) or Equation ( 129) holds, there exists qI r0, Θ pp1,φ1q s such that InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q. -If minte p1 , e p2 u e p2 , then pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq for every q r0, Θ pp1,φ1q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we have: U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q p0q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q p0q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, InvSoC q1 pp1,φ1q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q from lemma 4.1 (130)

Since either Equation (128) or Equation (130) holds, there exists qI r0, Θ pp1,φ1q s such that InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q. 2. If q is such that maxte p1 , e p2 u ¤ q ¤ Q pp1,φ1q pp2,φ2q , then -If minte p1 , e p2 u e p1 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p1 q ¥ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ InvSoC q2 pp2,φ2q pq ¡ e p1 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq rq ¡ e p2 , q ¡ e p1 s, InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ InvSoC q2 pp2,φ2q pqq ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq rq ¡ e p2 , q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq from lemma 4.1 ñ dq r0, q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq

Since Equation ( 127) holds, we therefore have InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq for every q r0, Θ pp1,φ1q s.

Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q
pp2,φ2q pq, q 1 , q 2 , 2q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I rq ¡ e p2 , q ¡ e p1 s, InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ InvSoC q2 pp2,φ2q pq I q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I rq ¡ e p2 , q ¡ e p1 s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q from lemma 4.1 ñ hq I r0, q ¡ e p1 s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q (131)

Since either Equation ( 128) or Equation ( 131) holds, there exists qI r0, Θ pp1,φ1q s such that InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q. -If minte p1 , e p2 u e p2 , then pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq for every q r0, Θ pp1,φ1q s.

Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q
pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p2 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p2 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, InvSoC q1 pp1,φ1q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ hq I r0, q ¡ e p2 s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q from lemma 4.1 (132)

Since either Equation (128) or Equation (132) holds, there exists qI r0, Θ pp1,φ1q s such that InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q.

• Case 2.2: n p2 0 (q 1 e opp1q,µp 1 p0q ) First, observe that if q rQ pp1,φ1q pp2,φ2q , Q pp1,φ1q pp2,φ2q s is such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¥ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

T pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ô min qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tInvSoC q1 pp1,φ1q pqq ¡ t p2 u ¥ ∆ p1 pqq ñ dq Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pqq ¡ t p2 ¥ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ñ dq Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pqq ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ t p2 ñ dq Ω pp1,φ1q pp2,φ2q pqq, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq from lemma 4.1 ñ dq rmaxt0, q ¡ minte p1 , e p2 uu, Θ q1 pp1,φ1q s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq

Second, observe that if q rQ pp1,φ1q pp2,φ2q , Q pp1,φ1q pp2,φ2q s is such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

T pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq ô max qΩ pp 1 ,φ 1 q pp 2 ,φ 2 q pqq tInvSoC q1 pp1,φ1q pqq ¡ t p2 u ¡ ∆ p1 pqq ñ hq I Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pq I q ¡ t p2 ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ñ hq I Ω pp1,φ1q pp2,φ2q pqq, InvSoC q1 pp1,φ1q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ t p2 ñ hq I Ω pp1,φ1q pp2,φ2q pqq, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q from lemma 4.1 ñ hq I rmaxt0, q ¡ minte p1 , e p2 uu, Θ q1 pp1,φ1q s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q (134)

Let q rQ pp1,φ1q pp2,φ2q , minte p1 , e p2 us be such that T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¥ ∆ p1 pqq ¡ ∆ p2 pqq and T pp1,φ1q pp2,φ2q pq, q 1 , q 2 q ¡ ∆ p1 pqq ¡ ∆ p2 pqq. Equation ( 133) is equivalent to dq r0, Θ pp1,φ1q s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq and Equation ( 134) is equivalent to hq I r0, Θ pp1,φ1q s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q. This means that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q. Let q rminte p1 , e p2 u, Q pp1,φ1q pp2,φ2q s be such that mintT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2qu ¥ ∆ p1 pqq¡∆ p2 pqq and maxtT pp1,φ1q pp2,φ2q pq, q 1 , q 2 q, U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2qu ¡ ∆ p1 pqq ¡ ∆ p2 pqq. We consider the two following subcases:

1. If q is such that minte p1 , e p2 u ¤ q ¤ maxte p1 , e p2 u, then -If minte p1 , e p2 u e p1 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ t p2 ¥ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ t p2 ñ dq r0, q ¡ e p1 s, InvSoC q1 pp1,φ1q pqq ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ t p2 ñ dq r0, q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq from lemma 4.1

Since Equation ( 133) holds, we therefore have InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq for every q r0, Θ pp1,φ1q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we have: U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ t p2 ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ t p2 ñ hq I r0, q ¡ e p1 s, InvSoC q1 pp1,φ1q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ t p2 ñ hq I r0, q ¡ e p1 s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q from lemma 4.1 (135)

Since either Equation (134) or Equation ( 135) holds, there exists qI r0, Θ pp1,φ1q s such that InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q. -If minte p1 , e p2 u e p2 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q p0q ¡ t p2 ¥ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ñ InvSoC q1 pp1,φ1q p0q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ t p2 ñ dq r0, q ¡ e p2 s, InvSoC q1 pp1,φ1q pqq ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ t p2 ñ dq r0, q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq from lemma 4.1

Since Equation ( 133) holds, we therefore have InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq for every q r0, Θ pp1,φ1q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q
pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p2 q ¡ t p2 ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p2 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ t p2 ñ hq I r0, q ¡ e p2 s, InvSoC q1 pp1,φ1q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ t p2 ñ hq I r0, q ¡ e p2 s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q from lemma 4.1 (136)

Since either Equation (134) or Equation (136) holds, there exists qI r0, Θ pp1,φ1q s such that InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q. 2. If q is such that maxte p1 , e p2 u ¤ q ¤ Q pp1,φ1q pp2,φ2q , then -If minte p1 , e p2 u e p1 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ t p2 ¥ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ t p2 ñ dq r0, q ¡ e p1 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq from lemma 4.1

Since Equation ( 133) holds, we therefore have InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq for every q r0, Θ pp1,φ1q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we have: U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ t p2 ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ t p2 ñ hq I r0, q ¡ e p1 s, InvSoC q pp1,φ1q pqq I ¡ InvSoC q pp2,φ2q pq I q from lemma 4.1

(137)

Since either Equation ( 134) or Equation ( 137) holds, there exists qI r0, Θ pp1,φ1q s such that InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q. -If minte p1 , e p2 u e p2 , then U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ t p2 ¥ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q qq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ t p2 ñ dq rq ¡ e p1 , q ¡ e p2 s, InvSoC q1 pp1,φ1q pqq ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ t p2 ñ dq rq ¡ e p1 , q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq from lemma 4.1 ñ dq r0, q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq

Since Equation ( 133) holds, we therefore have InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq for every q r0, Θ pp1,φ1q s. Observe that if q is such that U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq, then we have:

U pp1,φ1q
pp2,φ2q pq, q 1 , q 2 , 2q ¡ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p2 q ¡ t p2 ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q qq ñ InvSoC q1 pp1,φ1q pq ¡ e p2 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ t p2 ñ hq I r0, q ¡ e p2 s, InvSoC q1 pp1,φ1q pq I q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ t p2 ñ hq I r0, q ¡ e p2 s, InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q from lemma 4.1

Since either Equation (134) or Equation (138) holds, there exists qI r0, Θ pp1,φ1q s such that InvSoC q pp1,φ1q pq I q ¡ InvSoC q pp2,φ2q pq I q. Therefore, we can conclude that pp 2 , φ 2 q ¡ q pp 1 , φ 1 q.

C Example requiring a large number of CS copies Figure 6 depicts an example in which a CS is visited 4|I| times in the optimal solution. In this example, we have one depot (denoted by 0), one customer (denoted by 1), and two CSs (denoted by 2 and 3). The service time at the customer is equal to 0.5, the battery capacity of the EVs is equal to 10, and T max is set to 40.

The optimal solution to the E-VRP-NL contains one route (0, 2, 3, 2, 1, 2, 3, 2, 0) and the total time is equal to 37.7. In this solution, the SoC of the EV upon arriving for the first time at CS 2 is equal to 1. To visit customer 1, the EV must necessarily leave CS 2 with a fully charged battery. Fully recharging the battery from 1 to full capacity at CS 2 takes 11.2 time units. However, by detouring to CS 3, the EV can obtain a fully charged battery in less time. Indeed, if the EV charges one energy unit at CS 2, 10 energy units at CS 3, and finally two energy units at CS 2, then the EV is fully charged in 11 time units (taking the round trip time to CS 3 into account). Moreover, after visiting customer 1, the EV arrives at CS 2 with an empty battery. To return to the depot, the EV must necessarily leave CS 2 with a SoC equal to 90%. Charging nine energy units at CS 2 takes in that case 10.4 time units. However, if the EV charges two energy units at CS 2, 10 energy units at CS 3, and finally one energy unit at CS 2, then the target SoC for the EV at CS 2 is obtained in 10.2 time units. In this example with only one customer, CS 2 is visited four times.

E.2 Algorithmic details on the heuristic algorithm to solve the FRVCP

Ð Π 0 t Ð V f Ð false while true do for i 0 ¤ ¤ ¤ npΠq ¡ 1 do (f', Π I , Π best ) Ð insertCS(i,Π,Π best ) if f' true then f Ð true t I Ð getLBTime(Π I ) if t I t then t Ð t I Π Ð Π I end end end if f true then Π Ð Π else return Π best end end return null
Algorithm 7: reviseChargingAmount(Π 0 ) input : A route Π 0 visiting m CS(s) numbered in the order of visits from 1 to m (let f piq be the i-th CS visited) output: (f, Π) where Π is an energy-feasible route derived from Π 0 after optimizing the charging decisions, and f is equal to true if the maximum-duration limit is not exceeded, false otherwise. F Detailed computational results

Π Ð Π 0 ; /*

F.1 E-VRP-NL

We represent each instance using the code tcγ 1 cγ 2 sγ 3 cγ 4 #, where γ 1 is the method used to place the customers (i.e., 0) randomization, 1) combination of randomization and clustering, 2) clustering), γ 2 is the number of customers, γ 3 is the number of CSs, γ 4 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 , , 3 , 4). The symbol "Inf" means that the instance was proven infeasible, whereas the symbol "-" means that no feasible solution was found by the solver but the instance was not proven infeasible. The CPU time reported in column "Time (s)" in the following tables is in seconds and rounded to the nearest integer.

F.1.1 CS replication-based models

F.2 FRVCP

The average CPU time reported in column "Time (ms)" in the following tables is in milliseconds and rounded up to the nearest tenth.

F.3 FRVCP 1

The average CPU time reported in column "Time (ms)" in the following tables is in milliseconds and rounded up to the nearest tenth.

  ), (33), and (37) as node-based tracking constraints. Here we will model the time and SoC tracking by introducing arcbased variables. More precisely, we replace variables τ j and y j with the continuous variables τ ij and y ij , which respectively represent the time and SoC of an EV when it departs from node i V to travel on arc pi, jq A. If no vehicle travels on arc pi, jq, then both variables are 0. We model the time and SoC tracking through the following arc-based tracking constraints which replace the node-based tracking constraints in rF CS rep N ode s. The second model for the E-VRP-NL, denoted as rF CS rep Arc s, is as follows:

Figure 2 :

 2 Figure 2: Example of a SoC function associated with the supporting points (1, 1.5), (2, 3), and (5, 4).

  the FRVCP we allow an EV to visit several CSs between two successive non-CS nodes. Many studies focus on a more restricted version of the FRVCP in which an EV can visit at most one CS between two successive non-CS nodes. We denote this restricted version by FRVCP-1.

Figure 3 :

 3 Figure 3: Example of a fixed route vehicle charging problem with route (0,2,7,5,3,0).

Figure 4 :

 4 Figure 4: Graph considered when applying the labeling algorithm to solve the fixed route vehicle charging problem with route (0,2,7,5,3,0).

Figure 5 :

 5 Figure 5: Performance charts of the different models on small size instances

  have compared three MILP formulations for the electric vehicle routing problem with nonlinear charging function (E-VRP-NL), two of which are new. The first formulation was proposed by Montoya et al. (2017), where tracking the time and the SoC of each route was achieved through the use of node-based variables. We have first proposed a new formulation that uses arc-based tracking variables. Computational experiments have shown that this alternative tracking strategy drastically improves the results. This can be explained by the fact that the arc-based formulation had a much tighter LP relaxation gap, compared to the node-based formulation.

  plk ¤ w plk w pl,k 1 p P, l L p , k B µpplq zt0, b µpplq u

	∆ pl c pl ¡ c pl p P, l L p	(76)
	x p t0, 1u p P		(77)
	τ p ¥ 0, y p ¥ 0 p P		(78)
						dppq,0 ¨xp p P	(61)
	q pl c pl	ķB µp plq ķB	ķB µp plq ķB µp plq λ plk	λ plk a µpplqk λ plk c µpplqk ķB µp plq zt0u	p P, l L p p P, l L p w plk p P, l L p	(62) (63) (64)
		µp plq zt0u	w plk x p p P, l L p	(65)
	λ pl0 ¤ w pl1 p P, l L p	(66)
	λ plk ¤ w plk w pl,k 1 p P, l L p , k B µpplq zt0, b µpplq u	(67)
	λ plb µp plq ¤ w plb µp plq p P, l L p q pl λ plk a µpplqk p P, l L p ķB µp plq c pl λ plk c µpplqk p P, l L p ķB µp plq ķB ķB µp plq λ plk w plk p P, l L p µp plq zt0u ķB	(68) (69) (70) (71)
		µp plq zt0u	w plk x p p P, l L p	(72)
	λ i0 ¤ w pl1 p P, l L p	(73)
							(74)
	λ plb µp plq ¤ w plb µp plq p P, l L p	(75)

λ

  λ plk ¤ w plk w pl,k 1 p P pΠq, l L p , k B µpplq zt0, b µpplq u (100)λ plb µp plq ¤ w plb µp plq p P pΠq, l L p

	c pl	λ plk c µpplqk ķB µp plq zt0u w plk x p p P pΠq, l L p p P pΠq, l L p ķB µp plq λ plk w plk p P pΠq, l L p µp plq zt0u ķB µp plq ķB	(96) (97) (98)
	λ pl0 ¤ w pl1 p P pΠq, l L p	(99)
							(101)
	q pl c pl	ķB µp plq ķB	ķB µp plq ķB µp plq λ plk	λ plk a µpplqk λ plk c µpplqk ķB µp plq zt0u	p P pΠq, l L p p P pΠq, l L p w plk p P pΠq, l L p	(102) (103) (104)
		µp plq zt0u
							91)
	y p Qx p , p P πp0q,πp1q	(92)
	y p ¤ Qx p , p P pΠq ¤		(93)
	pP pΠq q pl	¥ t p x p ķB µp plq	ļ Lp λ plk a µpplqk ∆ pl	npΠq¡1 i1 p P pΠq, l L p g πpiq ¤ T max	(94) (95)

w plk x p p P pΠq, l L p (105) λ i0 ¤ w pl1 p P pΠq, l L p (106)

Table 1 :

 1 Proposition 3.1 demonstrates that the LP relaxation value of rF CS rep Arc s can never be less than that of rF CS rep N ode s. In our experiments, we computed the average relative gap between the linear relaxation value and the best lower bound provided by the linear relaxation of any of the formulations. The models with arc-based tracking variables yields tighter bounds than those with node-based tracking variables. Specifically, let z LB Arc and z LB N ode correspond to the average LP relaxation of rF CS rep Detailed computational results for the CS replication-based models on the 10-customer and 20-

	customer instances
	|I| β Model 10 1 rF CS rep Arc rF CS rep N ode s s 2 rF CS rep Arc s rF CS rep N ode s 3 rF CS rep Arc s rF CS rep N ode s 4 rF CS rep Arc s rF CS rep N ode s 20 1 rF CS rep Arc s rF CS rep N ode s 2 rF CS rep Arc s rF CS rep N ode s 3 rF CS rep Arc s rF CS rep N ode s 4 rF CS rep Arc s rF CS rep N ode s	#Inf #Unk #Opt #Opt Time (s) 10 0 0 10 5 10 0 0 10 252 5 0 0 15 433 5 0 4 11 699 31.8% Gap ---0 0 0 20 821 -0 0 8 12 1289 32.6% 0 0 4 16 215 14.4% 0 0 11 9 1423 34.0% 8 1 5 6 632 16.4% 3 6 11 0 1274 24.8% 2 4 8 6 1323 8.8% 1 5 13 1 5663 24.0% 2 3 10 5 2233 9.8% 0 5 14 1 8706 28.9% 0 5 10 5 1200 10.9% 0 5 15 0 -27.0%

Arc

s and rF CS rep N ode s. We observe that pz LB Arc ¡z LB N ode q{z LB Arc is approximately 57%. This seems to explain why the solver yields a better performance running on arc-based models.

Table 2 :

 2 Table 2 reports the same measures as in Table1for the path-based model. See F for the detailed results for each instance. Computational results for the path-based model on the 10-customer and 20-customer instances

	|I| #Inf #Unk #Opt #Opt Time (s) Gap
	10	0	0	0	20	229	-
	20	0	0	15	5	489 14.0%

Table 3 :

 3 Computational results for the different algorithms on 29,443 routes

	FRVCP

Table 4 :

 4

	Instance	#R ¦ /#R Max Gap R ¦ Gap	Instance	#R ¦ /#R Max Gap R ¦ Gap
	tc0c40s8cf0	2/7	-4.48% -0.765%	tc1c160s24cf3	1/17	-1.03%	-0.066%
	tc2c40s5cf2	1/6	-0.20% -0.024%	tc1c320s24cf2	2/36	-1.00%	-0.045%
	tc2c40s8cf2	1/6	-0.20% -0.025%	tc1c320s24cf3	1/30	-0.46%	-0.011%
	tc0c80s8cf1	1/10	-0.03% -0.003%	tc1c320s38cf2	1/33	-0.34%	-0.008%
	tc1c80s12cf2	1/8	-0.24% -0.038%	tc1c320s38ct3	1/30	-0.31%	-0.008%
	tc2c80s8cf4	1/10	-1.86% -0.235%	tc2c320s24cf0	1/38	-0.01% -0.0004%
	tc2c80s8ct3	1/9	-0.08% -0.009%	tc2c320s24ct4	1/32	-2.57%	-0.096%
	tc0c160s16cf4	3/18	-0.43% -0.074%	tc2c320s38cf4	1/32	-0.19%	-0.006%
	tc0c160s16ct4	3/18	-0.31% -0.052%	tc2c320s38ct0	2/41	-0.31%	-0.006%
	tc0c160s24cf4	3/18	-0.44% -0.076%	tc2c320s38ct1	1/28	-0.02%	-0.001%
	tc0c160s24ct4	4/18	-3.37% -0.197%	tc2c320s38ct4	1/32	-0.09%	-0.004%
	tc1c160s16cf3	1/17	-5.64% -0.347%

E-VRP-NL solutions of

[START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF]

improved by applying the labeling algorithm on each route

  A route Π 0 tπp0q, πp1q, ¤ ¤ ¤ , πpiq, ¤ ¤ ¤ , πpnpΠq ¡ 1q, πpnpΠqqu output: (f, Π, Π best ) where Π is the best route obtained if if boolean f is equal to true, or the initial route Π 0 otherwise. A route Π 0 tπp0q, πp1q, ¤ ¤ ¤ , πpiq, ¤ ¤ ¤ , πpn ¡ 1q, πpnqu A position i where to insert a CS of set FThe current best energy-feasible route Π best output: (f, Π, Π best ). Π is the route yielding the minimum lower bound on the duration to travel it if f is equal to true, or the initial route Π 0 otherwise. Π best represents the current best energy-feasible route and may have been updated pΠ I , jq ¥ 0 Y pΠ I , nextpΠ I , jqq ¥ 0 Algorithm 6: locateCS2(Π 0 , Π best ) input : A route Π 0 tπp0q, πp1q, ¤ ¤ ¤ , πpiq, ¤ ¤ ¤ , πpjq, ¤ ¤ ¤ , πpn ¡ 1q, πpnquThe current best energy-feasible route Π best output: Π best where Π best represents the energy-feasible route with the minimum duration under the H charging rule policy Π

	Algorithm 4: locateCS1(Π 0 )
		Π best represents the current best energy-feasible route with the minimum duration (under the H charging policy)
		that may have been computed during the procedure
	Π Ð Π 0
	Π best Ð null i Ð 0, n Ð npΠq Algorithm 5: insertCS(i, Π 0 , Π best )
	while i n do
	if the arc connecting πpiq and πpi 1q was eliminated during the preprocessing procedure then
		pf, Π, Π best q Ð insertCS(i,Π,Π best )
		if f false then
		return (false, Π 0 , Π best ) by the procedure
	Π Ð Π 0	end n Ð n 1
	else
		i Ð i 1
	end
	end
	return (true, Π, Π best )
		©
		then
		if tI t then
		t Ð tI Π Ð Π I
		end
		end
		end
	end
	end
	return (f,Π,Π best )

input : input :

s Ð getEnergyDeficit(Π 0 ) t best Ð mintTmax, getTime(Π best )u t Ð V, f Ð false for j F do

if the arcs connecting πpiq to j and j to πpi 1q were not eliminated during the preprocess then

Π I Ð copyAndInsert(Π 0 , j, i) s I Ð getEnergyDeficit(Π I ) tI Ð getLBTime(Π I ) ; /*

Let nextpΠ, jq be the first CS or depot node following node j in route Π */ if tI t best s I s ¡ Y f Ð true end if s I 0 then t I Ð getTime(Π I ) if t I t best then t best Ð t I Π best Ð Π I

  Let SoC 1 piq, SoC 2 piq, and chpiq be the arrival and departure SoC and the amount of energy charged at the Let b 1 piq and b 2 piq be the segment associated with the SoC SoC 1 piq and SoC 2 piq for the i-thCS visited */ i Ð m ¡ 1, t Ð getTime(Π best ) while i ¥ 1 do while ρ f piq,b 2 piq ¡ ρ f pi 1q,b 1 pi 1q chpi 1q ¡ 0 doTransfer from the pi 1q-th to the i ¡ th CS of route Π the amount of energy mintchpi 1q, a f piq,b 2 piq ¡ SoC 2 piqu Update accordingly SoC 1 pi 1q, SoC 2 piq, chpi 1q, b 2 piq, b 1 pi 1q and the duration t

	i-th CS visited	*/
	; /* end
	i Ð i ¡ 1
	end
	if t ¤ Tmax then
	return (true, Π)
	end
	return (false, Π 0 )

Table 7 :

 7 Detailed computational results on the 10-customer instances for the CS replication-based models

	(β 1)
	Instance	rF CS rep Arc Obj Bound LP Bound Time (s) Obj Bound LP Bound Time (s) s rF CS rep N ode s
	tc0c10s2cf1	Inf	Inf	Inf	0	Inf	Inf	6.02
	tc0c10s2ct1 12.30	12.30	9.05	1 12.30	12.30	6.33
	tc0c10s3cf1	Inf	Inf	12.06	0	Inf	Inf	6.02
	tc0c10s3ct1 10.80	10.80	7.72	1 10.80	10.80	5.13
	tc1c10s2cf2	9.03	9.03	6.72	1 9.03	9.03	5.56
	tc1c10s2cf3	Inf	Inf	Inf	0	Inf	Inf	3.42
	tc1c10s2cf4	Inf	Inf	Inf	0	Inf	Inf	6.47
	tc1c10s2ct2 10.76	10.76	6.73	6 10.76	10.76	4.58
	tc1c10s2ct3	Inf	Inf	9.83	0	Inf	Inf	3.18
	tc1c10s2ct4	Inf	Inf	Inf	0	Inf	Inf	6.23
	tc1c10s3cf2	9.03	9.03	6.72	1 9.03	9.03	5.56
	tc1c10s3cf3	Inf	Inf	Inf	0	Inf	Inf	3.42
	tc1c10s3cf4 14.95	14.95	10.58	1 14.95	14.95	5.80
	tc1c10s3ct2 9.20	9.20	6.36	11 9.20	9.20	4.46
	tc1c10s3ct3 13.02	13.02	7.14	4 13.02	13.02	3.16
	tc1c10s3ct4 13.21	13.21	8.85	4 13.21	13.21	6.09
	tc2c10s2cf0	Inf	Inf	9.48	0	Inf	Inf	1.04
	tc2c10s2ct0	Inf	Inf	4.67	13	Inf	Inf	0.89
	tc2c10s3cf0	Inf	Inf	9.48	0	Inf	Inf	1.04
	tc2c10s3ct0 11.51	11.51	4.21	63 11.51	11.51	0.89

Table 8 :

 8 Detailed computational results on the 10-customer instances for the CS replication-based models

	(β 2)
	Instance	rF CS rep Arc Obj Bound LP Bound Time (s) Obj Bound LP Bound Time (s) s rF CS rep N ode s
	tc0c10s2cf1	Inf	Inf	11.74	2	Inf	Inf	6.02
	tc0c10s2ct1 12.30	12.30	8.94	19 12.30	12.30	6.33
	tc0c10s3cf1 20.50	20.50	11.55	42 20.50	20.50	6.02
	tc0c10s3ct1 10.80	10.80	7.69	3 10.80	10.80	5.13
	tc1c10s2cf2	9.03	9.03	6.72	3 9.03	9.03	5.56
	tc1c10s2cf3	Inf	Inf	10.79	72	Inf	Inf	3.42
	tc1c10s2cf4 16.14	16.14	12.15	5 16.14	16.14	6.47
	tc1c10s2ct2 10.75	10.75	6.73	34 10.75	9.29	4.58
	tc1c10s2ct3	Inf	Inf	8.52	38	Inf	Inf	3.18
	tc1c10s2ct4 13.83	13.83	10.68	5 13.83	13.83	6.23
	tc1c10s3cf2	9.03	9.03	6.72	3 9.03	9.03	5.56
	tc1c10s3cf3 16.94	16.94	10.71	81 16.94	16.94	3.42
	tc1c10s3cf4 14.90	14.90	10.57	15 14.90	14.90	5.80
	tc1c10s3ct2 9.20	9.20	6.36	57	9.2	8.21	4.46
	tc1c10s3ct3 13.02	13.02	7.07	36 13.08	6.80	3.16
	tc1c10s3ct4 13.21	13.21	8.75	22 13.21	13.21	6.09
	tc2c10s2cf0	Inf	Inf	8.94	255	Inf	Inf	1.04
	tc2c10s2ct0 13.84	13.84	4.67	1165 13.84	13.84	0.89
	tc2c10s3cf0	Inf	Inf	8.94	6156	Inf	Inf	1.04
	tc2c10s3ct0 11.51	11.51	4.2	654 11.51	5.17	0.89

Table 9 :

 9 Detailed computational results on the 10-customer instances for the CS replication-based models

	(β 3)
	Instance	rF CS rep Arc Obj Bound LP Bound Time (s) Obj Bound LP Bound Time (s) s rF CS rep N ode s
	tc0c10s2cf1 19.75	19.75	11.74	21 19.75	19.75	6.02
	tc0c10s2ct1 12.30	12.30	8.94	21 12.30	12.30	6.33
	tc0c10s3cf1 19.75	19.75	11.55	1685 19.75	19.75	6.02
	tc0c10s3ct1 10.80	10.80	7.69	9 10.80	10.80	5.13
	tc1c10s2cf2	9.03	9.03	6.72	4 9.03	9.03	5.56
	tc1c10s2cf3 16.37	16.37	10.79	39 16.37	16.37	3.42
	tc1c10s2cf4 16.10	16.10	12.14	10 16.10	16.10	6.47
	tc1c10s2ct2 10.75	10.75	6.73	108 10.75	8.42	4.58
	tc1c10s2ct3 13.17	13.17	8.48	11 13.17	13.17	3.18
	tc1c10s2ct4 13.83	13.83	10.68	17 13.83	13.83	6.23
	tc1c10s3cf2	9.03	9.03	6.72	6 9.03	9.03	5.56
	tc1c10s3cf3 16.37	16.37	10.71	529 16.37	12.76	3.42
	tc1c10s3cf4 14.90	14.90	10.57	37 14.90	14.90	5.80
	tc1c10s3ct2 9.20	9.20	6.36	80 9.20	7.82	4.46
	tc1c10s3ct3 13.02	13.02	7.07	504 13.02	4.90	3.16
	tc1c10s3ct4 13.21	13.21	8.75	34 13.21	13.21	6.09
	tc2c10s2cf0 22.12	22.12	8.94	3647 22.12	19.56	1.04
	tc2c10s2ct0 12.45	12.45	4.67	334 12.45	7.14	0.89
	tc2c10s3cf0 22.12	22.12	8.94	8604 22.12	16.95	1.04
	tc2c10s3ct0 11.51	11.51	4.20	712 11.51	4.34	0.89

Table 10 :

 10 Detailed computational results on the 10-customer instances for the CS replication-based models

	(β 4)
	Instance	rF CS rep Arc Obj Bound LP Bound Time (s) Obj Bound LP Bound Time (s) s rF CS rep N ode s
	tc0c10s2cf1 19.75	19.75	11.74	96 19.75	19.75	6.02
	tc0c10s2ct1 12.30	12.30	8.94	81 12.30	12.30	6.33
	tc0c10s3cf1 19.75	17.68	11.55	10800 19.75	14.47	6.02
	tc0c10s3ct1 10.80	10.80	7.69	37 10.80	10.80	5.13
	tc1c10s2cf2	9.03	9.03	6.72	12 9.03	9.03	5.56
	tc1c10s2cf3 16.37	16.37	10.79	82 16.37	13.86	3.42
	tc1c10s2cf4 16.10	16.10	12.14	26 16.10	16.10	6.47
	tc1c10s2ct2 10.75	10.75	6.73	122 10.75	8.33	4.58
	tc1c10s2ct3 13.17	13.17	8.48	41 13.17	13.17	3.18
	tc1c10s2ct4 13.83	13.83	10.68	17 13.83	13.83	6.23
	tc1c10s3cf2	9.03	9.03	6.72	10 9.03	9.03	5.56
	tc1c10s3cf3 16.37	16.37	10.71	1663 16.37	11.95	3.42
	tc1c10s3cf4 14.90	14.90	10.57	76 14.90	14.90	5.80
	tc1c10s3ct2 9.20	9.20	6.36	137 9.20	7.66	4.46
	tc1c10s3ct3 13.02	13.02	7.07	207 13.08	5.28	3.16
	tc1c10s3ct4 13.21	13.21	8.75	45 13.21	11.56	6.09
	tc2c10s2cf0 21.77	19.62	8.94	10800 21.77	12.70	1.04
	tc2c10s2ct0 12.45	12.45	4.67	789 12.45	6.32	0.89
	tc2c10s3cf0 21.77	16.86	8.94	10800 21.83	11.98	1.04
	tc2c10s3ct0 11.51	9.83	4.20	10800 11.56	4.90	0.89

Table 11 :

 11 Detailed computational results on the 20-customer instances for the CS replication-based models

	(β 1)
	Instance	rF CS rep Arc Obj Bound LP Bound Time (s) Obj Bound LP Bound Time (s) s rF CS rep N ode s
	tc0c20s3cf2	Inf	Inf	Inf	0	Inf	Inf	6.89
	tc0c20s3ct2	Inf	Inf	11.65	460	-	-	6.96
	tc0c20s4cf2	Inf	Inf	Inf	0	Inf	Inf	6.86
	tc0c20s4ct2	Inf	Inf	10.89	10800	-	-	6.62
	tc1c20s3cf1 20.75	17.76	14.00	10800 20.78	14.68	8.02
	tc1c20s3cf3 16.80	16.10	9.81	10800 16.84	11.55	5.20
	tc1c20s3cf4	Inf	Inf	13.21	18	Inf	Inf	9.52
	tc1c20s3ct1 24.42	19.23	14.85	10800 24.42	16.20	8.02
	tc1c20s3ct3 12.60	12.60	9.35	189 12.67	10.12	5.45
	tc1c20s3ct4 16.24	16.24	13.62	33 16.24	14.55	10.18
	tc1c20s4cf1 19.29	15.82	11.85	10800 19.68	12.97	7.63
	tc1c20s4cf3 16.80	16.80	9.81	6474 16.91	10.02	5.19
	tc1c20s4cf4 19.09	19.09	13.18	226 19.09	18.15	9.52
	tc1c20s4ct1 23.19	17.59	12.74	10800 23.19	14.54	7.65
	tc1c20s4ct3 14.43	14.43	9.69	553 14.43	10.54	5.23
	tc1c20s4ct4 17.14	17.14	13.00	51 17.14	16.41	10.08
	tc2c20s3cf0	Inf	Inf	11.32	1	-	-	2.03
	tc2c20s3ct0	Inf	Inf	11.56	12	-	-	1.95
	tc2c20s4cf0	Inf	Inf	11.13	813	-	-	2.03
	tc2c20s4ct0	Inf	Inf	11.43	18	-	-	1.94

Table 12 :

 12 Detailed computational results on the 20-customer instances for the CS replication-based models

	(β 2)
	Instance	rF CS rep Arc Obj Bound LP Bound Time (s) Obj Bound LP Bound Time (s) s rF CS rep N ode s
	tc0c20s3cf2	Inf	Inf	14.89	4	-	-	6.89
	tc0c20s3ct2 17.08	15.34	11.64	10800 17.31	12.45	6.96
	tc0c20s4cf2	Inf	Inf	14.54	6	Inf	Inf	6.86
	tc0c20s4ct2 16.99	14.92	10.87	10800 17.21	11.97	6.62
	tc1c20s3cf1 17.53	17.53	13.91	8909 17.63	14.25	8.02
	tc1c20s3cf3 16.44	14.21	9.72	10800 16.53	10.07	5.19
	tc1c20s3cf4 17.00	17.00	13.21	87 17.00	17.00	9.41
	tc1c20s3ct1 19.49	19.00	14.74	10800 19.71	15.73	8.02
	tc1c20s3ct3 12.60	12.60	9.24	1086 12.60	9.29	5.45
	tc1c20s3ct4 16.21	16.21	13.35	94 16.21	14.44	10.17
	tc1c20s4cf1 16.40	15.42	11.8	10800 16.48	12.92	7.63
	tc1c20s4cf3 16.44	14.30	9.72	10800 16.94	10.04	5.19
	tc1c20s4cf4 17.00	17.00	13.17	174 17.01	15.18	9.41
	tc1c20s4ct1 18.36	17.25	12.6	10800 18.98	13.53	7.65
	tc1c20s4ct3 14.43	13.42	9.53	10800 14.43	10.43	5.23
	tc1c20s4ct4 17.00	17.00	12.78	222 17.00	15.73	10.02
	tc2c20s3cf0	-	-	11.18	10800	-	-	2.03
	tc2c20s3ct0	-	-	11.04	10800	-	-	1.95
	tc2c20s4cf0	-	-	11.02	10800	-	-	2.03
	tc2c20s4ct0	-	-	10.97	10800	-	-	1.94

Table 13 :

 13 Detailed computational results on the 20-customer instances for the CS replication-based models (β 3)

	Instance	rF CS rep Arc Obj Bound LP Bound Time (s) Obj Bound LP Bound Time (s) s rF CS rep N ode s
	tc0c20s3cf2	Inf	Inf	14.89	9074	-	-	6.89
	tc0c20s3ct2 17.25	15.06	11.63	10800 17.33	12.32	6.96
	tc0c20s4cf2	Inf	Inf	14.54	3619	-	-	6.86
	tc0c20s4ct2 16.99	14.76	10.87	10800 17.04	10.25	6.62
	tc1c20s3cf1 17.49	16.98	13.91	10800 17.49	14.31	8.02
	tc1c20s3cf3 16.44	14.04	9.72	10800 16.48	9.94	5.19
	tc1c20s3cf4 17.00	17.00	13.21	110 17.00	17.00	9.41
	tc1c20s3ct1 19.38	18.78	14.73	10800 20.06	15.36	8.02
	tc1c20s3ct3 12.60	12.60	9.24	1748 12.60	9.24	5.45
	tc1c20s3ct4 16.21	16.21	13.24	239 16.21	14.16	10.17
	tc1c20s4cf1 16.38	15.62	11.8	10800 16.40	12.96	7.63
	tc1c20s4cf3 16.44	13.83	9.72	10800 16.48	9.67	5.19
	tc1c20s4cf4 17.00	17.00	13.17	215 17.87	14.83	9.41
	tc1c20s4ct1 18.02	17.17	12.59	10800 18.45	13.64	7.65
	tc1c20s4ct3 14.43	13.34	9.52	10800 14.43	10.41	5.23
	tc1c20s4ct4 17.00	17.00	12.78	626 18.36	14.41	10.02
	tc2c20s3cf0	-	-	11.18	10800	-	-	2.03
	tc2c20s3ct0 25.79	20.90	10.91	10800 25.84	10.17	1.95
	tc2c20s4cf0	-	-	11.02	10800	-	-	2.03
	tc2c20s4ct0	-	-	10.87	10800	-	-	1.94

Table 14 :

 14 Detailed computational results on the 20-customer instances for the CS replication-based models (β 4)

	Instance	rF CS rep Arc Obj Bound LP Bound Time (s) Obj Bound LP Bound Time (s) s rF CS rep N ode s
	tc0c20s3cf2	-	-	14.89	10800	-	-	6.89
	tc0c20s3ct2 17.08	14.96	11.63	10800 17.31	12.19	6.96
	tc0c20s4cf2	-	-	14.54	10800	-	-	6.86
	tc0c20s4ct2 16.99	14.63	10.87	10800 17.72	10.16	6.62
	tc1c20s3cf1 17.49	16.94	13.91	10800 17.49	14.25	8.02
	tc1c20s3cf3 16.44	13.99	9.72	10800 16.95	9.36	5.19
	tc1c20s3cf4 17.00	17.00	13.21	126 17.00	16.82	9.41
	tc1c20s3ct1 18.94	18.42	14.73	10800 19.68	14.93	8.02
	tc1c20s3ct3 12.60	12.60	9.24	4177 12.67	9.29	5.45
	tc1c20s3ct4 16.21	16.21	13.14	357 16.23	13.70	10.17
	tc1c20s4cf1 16.38	15.17	11.8	10800 16.38	12.91	7.63
	tc1c20s4cf3 16.80	13.41	9.72	10800 16.89	10.30	5.19
	tc1c20s4cf4 17.00	17.00	13.17	463 17.00	14.99	9.41
	tc1c20s4ct1 17.81	16.60	12.59	10800 17.87	13.73	7.65
	tc1c20s4ct3 14.43	13.28	9.51	10800 14.50	10.40	5.23
	tc1c20s4ct4 17.00	17.00	12.78	876 17.17	14.31	10.02
	tc2c20s3cf0	-	-	11.18	10800	-	-	2.03
	tc2c20s3ct0 25.79	20.85	10.91	10800 26.24	10.13	1.95
	tc2c20s4cf0	-	-	11.02	10800	-	-	2.03
	tc2c20s4ct0	-	-	10.87	10800	-	-	1.94
	F.1.2 Path-based model

Table 15 :

 15 Detailed computational results on the 10-customer instances for the path-based model

	Instance	#CSPs	rF P ath Arc s Obj Bound LP Bound Time (s)
	tc0c10s2cf1	197 19.75	19.75	12.26	11
	tc0c10s2ct1	245 12.30	12.30	9.36	10
	tc0c10s3cf1	197 19.75	19.75	12.26	11
	tc0c10s3ct1	358 10.80	10.80	8.25	7
	tc1c10s2cf2	192 9.03	9.03	6.81	4
	tc1c10s2cf3	204 16.37	16.37	11.26	56
	tc1c10s2cf4	253 16.10	16.10	12.74	10
	tc1c10s2ct2	283 10.75	10.75	6.99	118
	tc1c10s2ct3	249 13.17	13.17	9.01	4
	tc1c10s2ct4	259 13.83	13.83	11.15	4
	tc1c10s3cf2	192 9.03	9.03	6.81	4
	tc1c10s3cf3	218 16.37	16.37	11.25	48
	tc1c10s3cf4	336 14.90	14.90	11.17	13
	tc1c10s3ct2	448 9.20	9.20	6.61	45
	tc1c10s3ct3	350 13.02	13.02	7.33	27
	tc1c10s3ct4	317 13.21	13.21	8.99	10
	tc2c10s2cf0	238 21.77	21.77	10.96	348
	tc2c10s2ct0	265 12.45	12.45	4.94	794
	tc2c10s3cf0	238 21.77	21.77	10.96	350
	tc2c10s3ct0	366 11.51	11.51	4.41	2716

Table 16 :

 16 Detailed computational results on the 20-customer instances for the path-based model

	Instance	#CSPs	rF P ath Arc s Obj Bound LP Bound Time (s)
	tc0c20s3cf2	939 27.49	22.9	15.34	10800
	tc0c20s3ct2	1284 17.08	14.93	11.66	10800
	tc0c20s4cf2	1036 27.47	22.97	15.34	10800
	tc0c20s4ct2	1350 16.99	14.98	11.40	10800
	tc1c20s3cf1	1189 17.49	17.13	14.08	10800
	tc1c20s3cf3	1491 16.48	13.28	9.92	10800
	tc1c20s3cf4	1257 17.00	17.00	14.95	111
	tc1c20s3ct1	1255 18.94	18.55	14.98	10800
	tc1c20s3ct3	1373 12.60	12.60	9.39	1529
	tc1c20s3ct4	1104 16.21	16.21	13.91	163
	tc1c20s4cf1	1304 16.38	16.05	12.58	10800
	tc1c20s4cf3	1607 16.84	13.19	9.92	10800
	tc1c20s4cf4	1322 17.00	17.00	14.69	134
	tc1c20s4ct1	1516 18.02	16.67	13.72	10800
	tc1c20s4ct3	1645 14.43	13.60	9.69	10800
	tc1c20s4ct4	1330 17.00	17.00	14.54	506
	tc2c20s3cf0	1437 24.68	18.31	11.95	10800
	tc2c20s3ct0	1318 25.80	21.29	11.74	10800
	tc2c20s4cf0	1353 24.67	17.78	11.81	10800
	tc2c20s4ct0	1195 26.02	20.45	11.46	10800

Table 17 :

 17 Computational results for the different algorithms on 29,443 routes according to the number of CSs

	|F| Algorithm	#Feas #Opt	Gap	Time (ms)	|F| Algorithm	#Feas #Opt	Gap	Time (ms)
	2	MILP PATH	440	440 0.000%	14.2	8	MILP PATH	5316 0.000%	178.4
		LABEL	440	440 0.000%	0.8		LABEL	5316 0.000%	0.7
		HEURISTIC	440	435 0.006%	0.3		HEURISTIC	5225 0.174%	0.3
		HEURISTIC M2017	418	339 0.599%	0.2		HEURISTIC M2017	4412 1.453%	0.2
	3	MILP PATH	1569 1569 0.000%	27.7	12 MILP PATH	3218 0.000%	76.6
		LABEL	1569 1569 0.000%	0.7		LABEL	3218 0.000%	1.0
		HEURISTIC	1559 1475 0.055%	0.3		HEURISTIC	3202 0.199%	0.3
		HEURISTIC M2017	1498 1229 0.692%	0.2		HEURISTIC M2017	2996 0.959%	0.2
	4	MILP PATH	1079 1079 0.000%	27.1	16 MILP PATH	7861 0.000%	69.0
		LABEL	1079 1079 0.000%	0.7		LABEL	7861 0.000%	0.9
		HEURISTIC	1078 1047 0.025%	0.3		HEURISTIC	7613 0.236%	0.2
		HEURISTIC M2017	1046	876 0.807%	0.1		HEURISTIC M2017	6424 0.690%	0.2
	5	MILP PATH	1982 1982 0.000%	24.0	24 MILP PATH	7978 0.000%	123.1
		LABEL	1982 1982 0.000%	0.6		LABEL	7978 0.000%	1.5
		HEURISTIC	1948 1763 0.154%	0.2		HEURISTIC	7750 0.228%	0.3
		HEURISTIC M2017	1737 1422 0.903%	0.1		HEURISTIC M2017	6867 0.534%	0.3

Table 18 :

 18 Computational results for the different algorithms on 29,443 routes according to the number of nodes in the initial route (including the depot at the beginning and at the end)

	|npΠq| Algorithm	#Feas #Opt	Gap	Time (ms)	|npΠq| Algorithm	#Feas #Opt	Gap	Time (ms)
	3	MILP PATH	58	58 0.000%	19.7	10	MILP PATH	5828 5828 0.000%	99.9
		LABEL	58	58 0.000%	0.5		LABEL	5828 5828 0.000%	1.0
		HEURISTIC	58	51 0.285%	0.1		HEURISTIC	5631 4427 0.193%	0.2
		HEURISTIC M2017	58	49 0.926%	0.1		HEURISTIC M2017	4843 3505 0.541%	0.2
	4	MILP PATH	298	298 0.000%	23.4	11	MILP PATH	4991 4991 0.000%	166.6
		LABEL	298	298 0.000%	0.6		LABEL	4991 4991 0.000%	1.0
		HEURISTIC	298	278 0.090%	0.2		HEURISTIC	4905 4250 0.105%	0.2
		HEURISTIC M2017	292	246 1.303%	0.1		HEURISTIC M2017	4640 3658 0.284%	0.2
	5	MILP PATH	710	710 0.000%	29.4	12	MILP PATH	3942 3942 0.000%	77.5
		LABEL	710	710 0.000%	0.7		LABEL	3942 3942 0.000%	1.1
		HEURISTIC	706	613 0.250%	0.2		HEURISTIC	3907 3579 0.040%	0.3
		HEURISTIC M2017	678	512 2.028%	0.1		HEURISTIC M2017	3813 3322 0.139%	0.2
	6	MILP PATH	1227 1227 0.000%	43.6	13	MILP PATH	1523 1523 0.000%	77.5
		LABEL	1227 1227 0.000%	0.9		LABEL	1523 1523 0.000%	1.2
		HEURISTIC	1193 1010 0.298%	0.3		HEURISTIC	1515 1441 0.023%	0.2
		HEURISTIC M2017	1158	825 2.431%	0.2		HEURISTIC M2017	1501 1425 0.020%	0.2
	7	MILP PATH	2084 2084 0.000%	59.5	14	MILP PATH	318	318 0.000%	55.2
		LABEL	2084 2084 0.000%	0.9		LABEL	318	318 0.000%	1.3
		HEURISTIC	2051 1648 0.291%	0.3		HEURISTIC	317	306 0.012%	0.2
		HEURISTIC M2017	1825 1159 1.911%	0.2		HEURISTIC M2017	318	301 0.012%	0.2
	8	MILP PATH	3411 3411 0.000%	76.8	15	MILP PATH	23	23 0.000%	16.8
		LABEL	3411 3411 0.000%	0.9		LABEL	23	23 0.000%	1.9
		HEURISTIC	3318 2543 0.309%	0.3		HEURISTIC	23	22 0.007%	0.2
		HEURISTIC M2017	2571 1659 1.753%	0.2		HEURISTIC M2017	23	23 0.000%	0.2
	9	MILP PATH	5030 5030 0.000%	104.0
		LABEL	5030 5030 0.000%	0.9
		HEURISTIC	4893 3624 0.314%	0.3
		HEURISTIC M2017	3678 2420 1.077%	0.2

Table 19 :

 19 Computational results for the different algorithms on 29,443 routes according to the number of CSs

	|F| Algorithm	#Feas #Opt	Gap	Time (ms)	|F| Algorithm	#Feas #Opt	Gap	Time (ms)
	2	MILP PATH 1	440	440 0.000%	9.9	8	MILP PATH 1	5316 5316 0.000%	23.1
		LABEL 1	440	440 0.000%	0.3		LABEL 1	5316 5316 0.000%	0.3
		MILP M2017	440	440 0.000%	13.1		MILP M2017	5316 5316 0.000%	39.2
	3	MILP PATH 1	1569 1569 0.000%	14.4	12 MILP PATH 1	3218 3218 0.000%	30.2
		LABEL 1	1569 1569 0.000%	0.4		LABEL 1	3218 3218 0.000%	0.5
		MILP M2017	1569 1569 0.000%	19.7		MILP M2017	3218 3218 0.000%	57.6
	4	MILP PATH 1	1079 1079 0.000%	14.8	16 MILP PATH 1	7861 7861 0.000%	31.7
		LABEL 1	1079 1079 0.000%	0.3		LABEL 1	7861 7861 0.000%	0.6
		MILP M2017	1079 1079 0.000%	21.4		MILP M2017	7861 7861 0.000%	58.4
	5	MILP PATH 1	1982 1982 0.000%	14.0	24 MILP PATH 1	7978 7978 0.000%	43.5
		LABEL 1	1982 1982 0.000%	0.3		LABEL 1	7978 7978 0.000%	1.0
		MILP M2017	1982 1982 0.000%	23.3		MILP M2017	7978 7978 0.000%	82.8

Table 20 :

 20 Computational results for the different algorithms on 29,443 routes according to the number of nodes in the initial route (including the depot at the beginning and at the end)

	|npΠq| Algorithm	#Feas #Opt	Gap	Time (ms)	|npΠq| Algorithm	#Feas #Opt	Gap	Time (ms)
	3	MILP PATH 1	58	58 0.000%	8.1	10	MILP PATH 1	5828 5828 0.000%	33.4
		LABEL 1	58	58 0.000%	0.2		LABEL 1	5828 5828 0.000%	0.6
		MILP M2017	58	58 0.000%	14.1		MILP M2017	5828 5828 0.000%	64.9
	4	MILP PATH 1	298	298 0.000%	9.0	11	MILP PATH 1	4991 4991 0.000%	36.3
		LABEL 1	298	298 0.000%	0.2		LABEL 1	4991 4991 0.000%	0.7
		MILP M2017	298	298 0.000%	18.7		MILP M2017	4991 4991 0.000%	53.7
	5	MILP PATH 1	710	710 0.000%	11.1	12	MILP PATH 1	3942 3942 0.000%	36.7
		LABEL 1	710	710 0.000%	0.2		LABEL 1	3942 3942 0.000%	0.8
		MILP M2017	710	710 0.000%	28.6		MILP M2017	3942 3942 0.000%	38.1
	6	MILP PATH 1	1227 1227 0.000%	14.1	13	MILP PATH 1	1523 1523 0.000%	33.4
		LABEL 1	1227 1227 0.000%	0.3		LABEL 1	1523 1523 0.000%	0.9
		MILP M2017	1227 1227 0.000%	41.3		MILP M2017	1523 1523 0.000%	22.2
	7	MILP PATH 1	2084 2084 0.000%	19.6	14	MILP PATH 1	318	318 0.000%	30.8
		LABEL 1	2084 2084 0.000%	0.4		LABEL 1	318	318 0.000%	1.0
		MILP M2017	2084 2084 0.000%	58.7		MILP M2017	318	318 0.000%	17.1
	8	MILP PATH 1	3411 3411 0.000%	24.7	15	MILP PATH 1	23	23 0.000%	13.8
		LABEL 1	3411 3411 0.000%	0.4		LABEL 1	23	23 0.000%	0.8
		MILP M2017	3411 3411 0.000%	67.6		MILP M2017	23	23 0.000%	5.5
	9	MILP PATH 1	5030 5030 0.000%	30.2
		LABEL 1	5030 5030 0.000%	0.5
		MILP M2017	5030 5030 0.000%	70.3
	F.4 E-VRP-NL

Table 21 :

 21 New BKS obtained for the E-VRP-NL by applying the labeling algorithm on each route of the solutions reported by[START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF]

	Instance	Old BKS New BKS
	tc0c40s8cf0	31.282	31.045
	tc2c40s8cf2	27.147	27.141
	tc2c40s5cf2	27.543	27.536
	tc0c80s8cf1	45.227	45.225
	tc1c80s12cf2	29.544	29.532
	tc2c80s8cf4	49.287	49.171
	tc2c80s8ct3	32.315	32.312
	tc0c160s16cf4	82.924	82.863
	tc0c160s16ct4	82.365	82.323
	tc0c160s24ct4	80.956	80.796
	tc0c160s24cf4	81.442	81.380
	tc1c160s16cf3	71.757	71.509
	tc1c160s24cf3	68.558	68.512
	tc1c320s24cf2	152.131	152.063
	tc1c320s24cf3	117.475	117.462
	tc1c320s38cf2	141.631	141.620
	tc1c320s38ct3	116.075	116.065
	tc2c320s24cf0	182.454	182.453
	tc2c320s24ct4	121.936	121.820
	tc2c320s38cf4	122.325	122.318
	tc2c320s38ct0	190.974	190.963
	tc2c320s38ct1	94.534	94.533
	tc2c320s38ct4	121.662	121.657

https://en.wikipedia.org/wiki/Better_Place\_(company) -Retrieved

2018-04-10 2 http://www.teslarati.com/tesla-shuts-down-battery-swap-program-for-superchargers/ -Retrieved 2018-04-10

We also remove all the arcs between CS copies of the same CS. More specifically, for every i F and j, l F I i , we remove the arc pj, lq from set A.

U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q p0q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¥ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q p0q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p2 s, InvSoC q1 pp1,φ1q pqq ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq r0, q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq from lemma 4.1Since Equation (127) holds, we therefore have InvSoC q

U pp1,φ1q pp2,φ2q pq, q 1 , q 2 , 2q ¥ ∆ p1 pqq ¡ ∆ p2 pqq ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¥ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ InvSoC q1 pp1,φ1q pq ¡ e p1 q ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq rq ¡ e p1 , q ¡ e p2 s, InvSoC q1 pp1,φ1q pqq ¡ Φ ¡1 µp 1 p0q pq ¡ e o,µp 1 p0q q ¥ InvSoC q2 pp2,φ2q pq ¡ e p2 q ¡ Φ ¡1 µp 2 p0q pq ¡ e o,µp 2 p0q q ñ dq rq ¡ e p1 , q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqq from lemma 4.1 ñ dq r0, q ¡ e p2 s, InvSoC q pp1,φ1q pqq ¥ InvSoC q pp2,φ2q pqqSince Equation (127) holds, we therefore have InvSoC q

 [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF]to set the number β of CS copies is not necessarily optimal. We show this for both the case with nonlinear (example 1) and the case with linear (example 2) charging function approximations.

In the first example (denoted example 1), we have one depot (denoted by 0), two customers (denoted by 1 and 2), and one CS (denoted by 3). The service time at the customers is equal to 0.5, the battery capacity of the EVs is equal to 10 and T max is set to 20. Tables 5a and5b show the driving time and energy consumption between the nodes. Figure 5c depicts the charging function at the CS. For β 0 (i.e., the CS cannot be visited) and β 1 (i.e., the CS can be visited at most once), the optimal solutions contains two routes p0, 1, 0q and p0, 2, 0q and the total time is equal to 17. According to the [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] procedure, β will then be fixed to 0. However, for β 2 (i.e., the CS can be visited at most twice), the optimal solution contains one route p0, 2, 3, 1, 3, 0q and the total time is equal to 16.8. In this solution, the EV arrives at CS 3 with an empty battery. We charge the battery up to 50% and 30% during the first and second visits, which takes 0.5 and 0.3 time units, respectively. Now consider the case where the charging function is linear. We introduce another example (denoted as example 2). We have one depot (denoted by 0), four customers (denoted by 1, 2, 3, and 4), and one CS (denoted by 5). The service time at the customers is equal to 0.5, the battery capacity of the EVs is equal to 10 and T max is set to 10. Tables 6a and6b show the driving times and energy consumption between the nodes.

Figure 6c shows the charging function at the CS 5.

Table 6: Data of example 2. For β 0 and β 1, the optimal solutions contains two routes p0, 1, 4, 0q and p0, 2, 3, 0q and the total time is equal to 10. According to the [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] procedure, β will then be fixed to 0. However, for β 2, the optimal solution contains one route p0, 1, 5, 2, 3, 5, 4, 0q and the total time is equal to 9.65. At the first visit to the CS, the arrival SoC is 35% and we charge the battery up to 100%, which takes 0.325 time units. At the second visit to the CS, the arrival SoC of the EV is 0 and we charge the battery up to 65%, which takes 0.325 time units.

Note that if we assume a full charging policy and a constant charging time δ ¥ 0 in Example 2, the solutions remain the same. The solution for β 0 or 1 has an objective value equal to 10. The solution for β 2 is equal to 9 2δ. Therefore a better solution is found as soon as δ 0.5.

In conclusion, the procedure of [START_REF] Montoya | The electric vehicle routing problem with nonlinear charging function[END_REF] is not guaranteed to be optimal for constant, linear, and nonlinear charging functions.

E Algorithmic details

E.1 Algorithmic details on the labeling algorithm for the FRVCP Algorithm 3: LabelingAlgo(Π) input : a route Π that is not energy-feasible output: a route Π I that is energy-feasible if found, the initial route otherwise forall l Π tFu do L 1 plq Ð r (L 1 plq stores the unprocessed labels associated with node l) ; forall l Π tFu do L 2 plq Ð r (L 2 plq stores the processed labels associated with node l) ; Add to L 1 pΠp0qq the label corresponding to a RP containing only the departing depot Πp0q R Ð tΠp0qu (R stores nodes with unprocessed labels) while R $ r do l Ð minElement(R) (removes from R the node associated with a RP pp, φq with the minimum total time) pp l , φ l q Ð deleteMin(L 1 plq) (removes from L 1 plq the label with the minimum total time) if l ΠpnpΠqq then if p l contains a CS then Set the SoC target value at the last CS of p l in such a way the EV arrives with an empty battery at l end Let Π I be the route corresponding to label pp l , φ l q return Π I end if L 1 plq r then R Ð Rztlu ; Update(L 2 plq, pp l , φ l q) forall pl, l I q Ĝij do if IsExtendable(pp l , φ l q,l I ) then (checks if RP pp, φq can be extended to node l I ) pp l I, φ l Iq ÐExtendppp l , φ l q, l I q (extends RP pp l , φ l q to node l I ) K ÐCreateLabelsppp l I, φ l Iqq (if l I is a CS and p l I contains more than a CS, then it creates relevant labels according to the supporting points of the SoC function associated with RP pp l I, φ l Iq, otherwise K only contains pp l I, φ l Iq)

forall pp, φq K do Update(L 1 pl I q, pp, φq) (inserts RP pp, φq in L 1 pl I q if it is not dominated according to the dominance rule introduced in Definition 4.1, and verifies if previously inserted labels become dominated) end if L 1 pl I q has changed and l I R then R Ð R tl I u ; end end end return Π