The Hybrid Electric Vehicle—Traveling Salesman Problem with time windows

The Hybrid Electric Vehicle - Traveling Salesman Problem with Time Windows

Journal Pre-proof

The Hybrid Electric Vehicle - Traveling Salesman Problem with Time Windows Christian Doppstadt, Achim Koberstein, Daniele Vigo PII: DOI: Reference:

S0377-2217(19)31071-9 https://doi.org/10.1016/j.ejor.2019.12.031 EOR 16239

To appear in:

European Journal of Operational Research

Received date: Accepted date:

28 August 2019 21 December 2019

Please cite this article as: Christian Doppstadt, Achim Koberstein, Daniele Vigo, The Hybrid Electric Vehicle - Traveling Salesman Problem with Time Windows, European Journal of Operational Research (2019), doi: https://doi.org/10.1016/j.ejor.2019.12.031

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Highlights • We model customer deliveries with time windows by hybrid electric vehicles • We develop a new heuristic solution method for this optimization problem • We create a set of new benchmark instances for this optimization problem • With various numerical experiments we challenge the heuristics performance • Comparing our results with regular combustion vehicles indicate saving potential

1

The Hybrid Electric Vehicle - Traveling Salesman Problem with Time Windows Christian Doppstadt∗, Achim Koberstein European University Viadrina, Germany

Daniele Vigo University di Bologna, Italy

Abstract In this paper, we extend the Hybrid Electric Vehicle – Traveling Salesman Problem (HEV-TSP) that deploys hybrid electric vehicles for customer delivery tours, by considering that customers must be served within given time windows. This feature makes the problem very difficult to solve. We developed a Variable Neighborhood Search based heuristic solution method, which is able to handle hybrid electric vehicle problems with a realistic number of customers. We introduce a large set of benchmark instances, representing typical delivery areas for small package shipping companies. Exact solutions for instances with a small number of customers are calculated by formulating the problem as an integer linear program and solving the instances with the standard solver CPLEX. The proposed heuristic achieves optimal solutions on the small and good quality solutions on larger instances. Furthermore, the results show that the profitability of hybrid electric vehicles highly depends on the structure of the delivery area and the number of customers to serve. Therefore, our heuristic does not only serve to support decision makers in the daily tour planning, but also in the evaluation of the profitability of hybrid electric vehicles for a specific delivery area structure. Keywords: Traveling Salesman; Hybrid Electrical Vehicles; Time Windows; ∗ Corresponding

author

Preprint submitted to European Journal of Operational Research

December 27, 2019

Small Package Shipping

1. Introduction Reducing exhaust gases and sooty particles is one of the most important challenges for the future, especially, within urban areas, where the health of the inhabitants is endangered. One way to archive this, is the use of electric vehicles (EVs) instead of vehicles with internal combustion engines (ICEs). However, in the present, state of art EVs have some major disadvantages: their range is limited and additional infrastructure for recharging the battery is required. These limitations can be eliminated or reduced by using hybrid electric vehicles (HEVs), which combine combustion and electric engines. HEVs are widely used as passenger cars, but up to now they are relatively uncommon as delivery vehicles. With our approach we want to support decision makers in small package shipping (SPS) companies to evaluate the benefits of HEVs for their specific requirements, depending on the structure of a given delivery area. Therefore, we introduce a new optimization problem dealing with HEVs for delivery tours and combine it with time window constraints for the customer deliveries. We developed a heuristic solution method, which is able to calculate the costs and savings for a given practical data set, for this type of problem. Moreover, we solve small instances of this problem exactly, to evaluate the solution quality of our heuristic. With our approach, we extend the ’Hybrid Electric Vehicle – Traveling Salesman Problem’ (HEV-TSP) presented in Doppstadt et al. (2016) by including time windows for the customer deliveries. This means, that a customer can only be served after a given starting time and that the service has to be completed before a given end time is reached. We name the resulting problem ’Hybrid Electric Vehicle –Traveling Salesman Problem with Time Windows’ (HEV-TSPTW). The time windows amend an additional real world requirement to the base problem, as in today’s society customers prefer to schedule a delivery time which fits them best. The aim within this new setting is still to determine which arcs

3

to use in which mode of operation to minimize the total costs. To the best of our knowledge, this problem statement has not been studied before. Mancini (2017) is considering hybrid vehicles, however, instead of charging the battery while driving, the batteries may be recharged at recharging stations during the delivery tour. Recently, additional research on pure electric vehicles has been performed, see e.g., Schneider et al. (2014), Desaulniers et al. (2016). In our model, the operation modes of the HEVs have different costs and travel times for each arc within a given delivery network. Most HEVs distinguish four different modes: pure combustion, pure electric, charging the battery while driving in combustion mode, and a boost mode, which combines combustion and electric engine (ARADEX AG, 2019). We follow these four default modes in our methodology. For an extensive review on the technical aspects of electric, hybrid, and fuel cell vehicles we refer to Chan (2002, 2007). Four modes of operation extends the graph of a classical TSP to a multigraph with four arcs between all vertices. Adding time windows, excludes some or even many arcs for feasible tours, depending on the arcs used before. These problem modifications change the structure of the TSP in a way that requires new steps within the solution process, which we provide with our new heuristic. An example of such a multigraph-network with a depot (0) and three customers (1-3) is given in Figure 1. The time window for the depot defines two times: one for when the vehicle is allowed to leave the depot and the other for when the vehicle must have completed its tour. Therefore, the time window for the depot implies a maximum tour duration. Having more than one link between each pair of vertices dramatically increases the number of possible tours for the already NP-hard TSP. Additionally, the time windows narrow down the set of feasible solutions within these round tours. In this setting, it is very difficult to find good or even valid solutions. Therefore, we introduce a heuristic solution method to be able to solve problems with a practical number of customers. The heuristic uses a Variable Neighborhood Search (VNS) as a diversification mechanism to reach distant regions within the solution space. A general overview on VNS can be found in 4

[b1,'e1]' Depot' Customer'

[b0,'e0]'

[b3,'e3]' electric' charging' combus1on' boost'

[bx,'ex]' Time'Window'

[b2,'e2]' Figure 1: Delivery Network with Different Modes of Operation and Time Windows

(Hansen et al., 2010; Gendreau & Potvin, 2010). Soon after, da Silva & Urrutia (2010) proposed a refined two-phase VNS that outperformed the two previous approaches for large TSPTWs instances. Similar to our introduction of time windows, Savelsbergh (1985) extended the TSP to the Traveling Salesman Problem with Time Windows (TSPTW) with regular combustion vehicles. He describes an efficient local search concept to evaluate the feasibility of a solution considering time windows. A compressed annealing local search heuristic for the TSPTW is presented in Ohlmann & Thomas (2007), whereas, L´ opez-Ib´ an ˜ez & Blum (2010) use a beam-search together with a learning component inspired by ant colony optimization, outperforming all previous solution methods. The remainder of the paper is structured as follows. In Section 2 we present a formal description of the problem and define an integer linear program. The concept of our heuristic solution method is explained in Section 3. Section 4 deals with the numerical experiments we conducted, including a description on how the benchmark instances were generated. In Section 5 we conclude the paper and give an outlook on further research topics to extend and improve the present approach. 5

2. Problem Description We define the HEV-TSPTW on a directed network G = (V, A), with V = {0, 1, ..., n, n + 1} as set of vertices, including the starting (0) and ending depot (n + 1) and A as set of arcs. To model time and charging level of the battery, we use a flow formulation as presented in Ascheuer et al. (2001). Note that the ending depot is a copy of the starting depot, as required for the flow model. For every vertex in V , we have a given service time si , representing the time required to serve a customer. Service times for the starting and ending depot simulate the time the driver additionally needs at the depot, e.g., for loading parcels into the vehicle. Beside the service time, each node has a time window assigned, consisting of a start time bi and an end time ei . If the driver reaches a node before the start time, it is allowed to wait until the time window opens. The waiting time does not cause any additional costs. However, as this is increasing the total route duration, a long waiting time increases the probability to violate a time window at a later point or to exceed the maximum tour duration. The end time indicates when the service at the customer has to be finished and therefore, the time difference between end and start time has to be greater than or equal to the service time. This is in contrast to other time window based routing problems, where the vehicle only has to reach a customer by the end time, e.g., the classical Vehicle Routing Problem with Time Windows (VRPTW) Desaulniers et al. (2014). In our opinion we provide a more graspable view on real world requirements, as from a customer satisfaction point of view, it is more important that the service is finished within a given time window than started. However, as the service times are fixed and known in advance, both formulations are mathematically equivalent and our representation can easily be transformed into the other by subtracting the service time from the time window end time for each customer. The end time of the depot replaces a maximum tour duration and can be seen as a limit of the working hours for a driver. For the set of arcs A, all arcs are associated with the different costs c and travel times t for each mode of operation. These costs and travel times are

6

indicated as cc and tc for the pure combustion mode, ce and te for the electric mode, ccc and tcc for the charging mode, and cb and tb for the boost mode. To increase the readability of the model, we define a subset of A as actually usable arcs A0 by starting node i and ending node j, where i ∈ {0, 1, ..., n} and j ∈ {1, 2, ..., n + 1} and i 6= j. We set a maximum battery capacity lmax , a charging rate rc , and a discharging rate rd for the HEV. The charging of the battery is increased by the product of the charging rate and the travel time for a given arc, if it is driven in charging mode. Likewise, the charging of the battery is decreased by the product of the discharging rate and the travel time, if an arc is driven in electric or boost mode. For each arc we use a binary decision variable xij , indicating whether the arc is used or not. We have additional variables (xc , xe , xcc , and xb ) to indicate the mode of operation used on that specific arc. As commonly used in flow formulations, the time of arrival at node i is given by the continuous variable wi . Similarly, the battery charging level at customer i is indicated by li . The objective for the HEV-TSPTW is to minimize the total costs for the arcs used within the delivery tour:

min

X

cc b b c c e e ccc ij xij + cij xij + cij xij + cij xij

(i,j)∈A0

To ensure the feasibility of a solution, we impose the following eighteen constraints:

7

b c e 0 xij , xcc ij , xij , xij , xij ∈ {0, 1} ∀(i, j) ∈ A b c e xcc ij + xij + xij + xij = xij n X i=0

n+1 X

j=1 n X j=1 n X

∀(i, j) ∈ A0

(1) (2)

xij = 1 ∀j = 1, ..., n; i 6= j

(3)

xij = 1 ∀i = 1, ..., n; i 6= j

(4)

x0j = 1

(5)

xi(n+1) = 1

(6)

i=1

w0 = 0

(7) ∀(i, j) ∈ A0

(8)

xbij (max{wi , bi } + si + tbij − wj ) = 0 ∀(i, j) ∈ A0

(9)

xcij (max{wi , bi } + si + tcij − wj ) = 0 ∀(i, j) ∈ A0

(10)

xeij (max{wi , bi } + si + teij − wj ) = 0 ∀(i, j) ∈ A0

(11)

wj + sj ≤ ej

(12)

cc xcc ij (max{wi , bi } + si + tij − wj ) = 0

∀j = 1, ..., n + 1

l0 = 0

(13)

cc c 0 xcc ij (min{li + tij r , lmax } − lj ) = 0 ∀(i, j) ∈ A

(14)

xbij (li − tbij rd − lj ) = 0 ∀(i, j) ∈ A0

(15)

xcij (li − lj ) = 0 ∀(i, j) ∈ A0

(16)

xeij (li − teij rd − lj ) = 0 ∀(i, j) ∈ A0

(17)

0 ≤ li ≤ lmax

(18)

∀i = 0, ..., n + 1

The first constraint defines the decision variables for the use of an arc as binary. Along with constraint (2) it is ensured that an arc can only be used in one mode of operation. Every customer must only have one in- and one outgoing arc, which is guaranteed by constraints (3) and (4). For the depot this requirement is identical, but it needs a slightly different formulation because of 8

the duplication of the depot as starting and ending depot, which is given in (5) and (6). We define the starting time at the depot as zero (7). The arrival time w at customer j is defined by the arrival time w at the preceding customer i plus the service time s at i and the driving time t between i and j in the mode used. If the opening of the time window b at i is later than the arrival w at i, the arrival time w at j is calculated as the opening of the time window b at i plus service time at i and driving time between i and j instead, which implies a waiting time before serving i of opening of the time window b at i minus arrival w at i (bi − wi ). These time flow constraints are imposed for each mode of operation by (8)-(11). Due to readability, constraints (8)-(11) are not yet linear. However, they can be linearized by standard big M inequality formulations. It is required that the service at any node j is completed before the end of the time window (12). Note, we assume for all nodes j that the opening of its time window plus the service time is equal or smaller than the end of the time window (bj + sj ≤ ej ), otherwise, the problem would be infeasible. As we are only considering HEVs and not Plug-in Hybrid Electric Vehicles (PHEV), we set the initial charging of the battery to zero (13). Note that, considering a different initial charing status resulting from a PHEV or from the remaining charge of a previous drive can be easily implemented by modifying constraint (13). While driving from node i to j in the charging mode, the charging of the battery is increased by the product of driving time and charging rate (14). Moreover, if the charging of the battery would exceed the maximum battery capacity lmax , the charging of the battery is set to the maximum charging level, which is also included in constraints (14). For the boost mode, the charging is decreased by the product of discharging rate and driving time between nodes i and j (15). The same applies for the pure electric mode (17). While driving in combustion mode, the charging between two nodes remains equal (16). Similar to constraints (8)-(11), constraints (14)-(17) are not linearized to keep the model in a more compact form. Finally, constraints (18) ensure that the charging of the battery remains between 0 and the maximum capacity 9

[41,'61]' Depot' Customer'

[23,'43]'

[00,'80]'

electric' charging' combus1on' boost'

[bx,'ex]' Time'Window'

[72,'82]' Figure 2: Example of a Solution in Different Modes of Operation

during the whole tour. To point out the structure of the problem definition, we give an example of a possible solution for a HEV route in Figure 2. The first arc is used to charge the battery. The second arc is driven in pure electric mode, which gains savings compared to the other modes. As the charge is sufficient for the rest of the tour, the third arc is driven in combustion mode, which is cheaper than the charging mode. Finally, the last arc has to be driven in the boost mode to fulfill the time window with respect to the tour duration constraint. This is still feasible, as the battery has a remaining charge from charging it on the first arc.

3. Solution Approach In this section we describe how our heuristic solution method is designed. It is divided into two phases; the first one deals with the construction of an initial solution and the second one demonstrates how we subsequently improve the solution with a parallelized VNS, see (Davidovic et al., 2013). Figure 3 gives a first overview of the process, which is described in the following subsections.

10

Ini$aliza$on)

Best)Global)) Solu$on)

Parallel)VNS)

Shake)k8$mes) 28Opt)with)Mode) Change)

Shake)k8$mes) 28Opt)with)Mode) Change)

Hill)Climbing) Insert)

Hill)Climbing) Insert)

Hill)Cimbing) χ8Mode)Change))

Hill)Cimbing) χ8Mode)Change))

Best)Local) Solu$on)

Best)Local) Solu$on)

Figure 3: Overview of Heuristic Solution Method

3.1. Initialization Phase Our initialization phase consists of three steps, each dealing with a different criterion relevant for a solution in our problem setting. Subsequently, we describe these steps in detail. Creating feasible solutions becomes quite difficult, taking time windows into account, especially, when they are tight and most or all customers are associated with such a tight time window. A solution could also be infeasible due to an invalid charging level of the battery. This means, that power is consumed while driving an arc in electrical or boost mode, but the power has not been generated driving one or more arcs in charging mode before. Finally, the tour could also exceed the total tour duration. As these factors highly narrow down the solution space, we allow infeasible solutions during the whole solution process and measure the possible types of violation, which allows us to compare even infeasible tours in terms of their quality and guide the heuristic towards feasible solutions. As time windows are considered to be the most critical constraint, we try to address this problem first during the construction of our initial solution. Most critical for the feasibility of a time window, is the end time of a customer. In contrast to the start time, where waiting is allowed, the end time is a fixed

11

[00,)80]) Depot)

[32,)52])

[53,)63]) combus'on)

[31,)71])

[00,)80])

[bx,)ex]) Time)Window)

Customer)

Figure 4: Initialization Phase A: Order Nodes by End Time Window

limit for each customer. This leads to the first step within our initialization method. We order the customers by their time window end time in ascending order, not considering the distances, costs or any other metrics. Customers without a time window are inserted between the customer with the latest time window end time (max value for ei with i ∈ A0 ) and before returning to the depot (n+1). If two customers have the same end time, their order is indeterminate and we consider them as equal. If one order is superior to the other due to other factors, it will most likely be detected during the second step within the initialization phase. However, our heuristic will deliver a deterministic order for the vertices at this point. The customers ordered in this way, are now linked in combustion mode, which is most likely to generate a feasible tour, as no constraint concerning the battery is violated. Figure 4 illustrates this step. Some time windows are not critical for a route. This means that moving the customers associated with these time windows does not effect the relevant tour parameters making a tour infeasible (e.g., later time windows). However, moving these customers might improve other tour metrics, especially the total costs. To detect this kind of savings, we introduce the second step. During a hill-climbing process, we use Insertion-moves to modify the position of single customers in the tour. This means that firstly, three arcs in a tour are removed; the two arcs in and out of the customer that is moved and the arc between the two other customers, where the moved customer is inserted. Afterwards, the route is reconnected. The parameters for accepting a new solution are that the current solution is not worsened in terms of time window violation and maximum tour duration violation, but the total costs are reduced. As the

12

[00,$80]$ Depot$

[31,$71]$

[32,$52]$ combus@on$

[53,$63]$

[00,$80]$

[bx,$ex]$ Time$Window$

Customer$

Figure 5: Initialization Phase B: Hill-Climbing Customer Insertion

costs are related to the distances driven, the distances driven are implicitly reduced, too. Charging violation does not have to be considered at this point, as all arcs are still driven in combustion mode. We repeat this procedure until no further improvement is found. In Figure 5, we give an example for the result of a single insertion move based on the solution of Figure 4, moving node 1 between the depot and node 2, assuming that the total costs for this solution are reduced. Up to this point, we left the battery charging out of consideration. During the initialization we want to cover two aspects regarding the charging. First, we try to make a solution feasible, which might require that the distance between two customers must be driven in boost mode due to a time window constraint. This requires that arcs must be driven in charging mode before. Second, the initial intention of HEVs is to generate tours that cause less costs than tours driven with conventional vehicles. Therefore, we are exchanging the mode of arcs to reduce the total costs. Our process is able to cover both aspects at once. However, to give higher priority to the first one, a feasible solution is always evaluated better than an infeasible solution, even if the costs of the infeasible solution are lower. To exchange the operation mode of arcs, we use the χ-Mode Change (χ-MC) presented in Doppstadt et al. (2016). The idea is to remove the driving mode of χ arcs at a time and test all other possible mode combinations in terms of feasibility and objective value. The χ-Mode Change (χ-MC) is the basis for the third and final step during the initialization. Similar to the customer insertion procedure, we also run the χ-MC as a hill-climbing process. The χ-MC tries to exchange all combinations of χ-arcs

13

[00,$80]$ Depot$ Customer$

[31,$71]$

[32,$52]$ electric$ charging$ combus6on$

[53,$63]$

[00,$80]$

[bx,$ex]$ Time$Window$

Figure 6: Initialization Phase C: Hill-Climbing χ-Mode Change

simultaneously, within the given solution. If no further improvement is found the initialization terminates and returns the best solution found. An example of a 3-MC mode is given in Figure 6, where the first two arcs from the solution in Figure 5 are changed to charging mode and the third arc is changed to electric mode. We determine the actual number for χ during the numerical studies. Note, that a solution can still be infeasible at this point. Although, in this case, the infeasibility will most likely be caused by a tour duration violation, it might also be impossible to cover the distance between two subsequent customers within the available time. The following improvement phase considers this fact and has some means to move towards feasible solutions beside reducing the total solution costs. These means are described in the next subsection. 3.2. Improvement Phase The improvement phase of our heuristic is based on the VNS, as introduced by Mladenovi´c & Hansen (1997). Attempts to simply adjust the Iterated Tabu Search (ITS) introduced in Doppstadt et al. (2016) for the HEV-TSP to the additional requirements of the HEV-TSPTW failed and we were not able, to achieve decent solution quality. Therefore, we decided to use a VNS variant instead, which has proven to provide very good results for the TSPTW (e.g., da Silva & Urrutia (2010)). To make use of the design of recent Central Processor Units (CPUs) in modern computer systems, having several cores, we designed a parallelization process around the basic VNS (P-VNS). This process creates a number – provided as parameter – of subprocesses, which run as independent threads and also takes

14

care of synchronizing and evaluating the results once all subprocesses are completed. These subprocesses can be seen as the actual VNS part of the heuristic, which is shown as simplified pseudo-code in Algorithm 1 and described in detail below. Algorithm 1 Variable Neighborhood Search 1: best local Solution ← input Solution, I ← input I, K ← input K 2:

i ← 1, k ← 1

3:

for i ≤ I do

4: 5:

for k ≤ K do

local Solution0

←

random

←

best

Solution

←

best

Solution

Solution

neighborhood

k

of

hill-climbing

Insert

on

hill-climbing

χ-MC

on

best local Solution 6:

local Solution00 local Solution0

7:

local Solution000 local Solution00

8:

if chooseSolutionAorB(local Solution000 , best local Solution) then best local Solution ← local Solution000

9:

k←1

10: 11:

else k ←k+1

12: 13:

end if

14:

end for

15:

i←i+1

16:

end for

17:

return best local Solution Beside the input solution, which has been generated during the initialization

phase, the VNS subprocess requires two additional parameters. Let I denote the number of iterations and K the number of neighborhoods to search (line 1). The control variables i and k are initialized with 1 (line 2). Subsequently, for each iteration and each neighborhood three steps are repeated: 15

[11,$51]$

[23,$43]$

[00,$80]$

[25,$45]$

[11,$51]$

[23,$43]$

[25,$45]$

[72,$82]$

[74,$84]$

[76,$86]$

[00,$80]$

[72,$82]$

[74,$84]$

Depot$ Customer$

[76,$86]$

electric$ charging$ combus=on$ boost$

[bx,$ex]$ Time$Window$

Figure 7: Improvement Phase A: Random 2-Opt with Mode Change

1. Draw a random solution out of the neighborhood k starting from the best local Solution and assign it to local Solution0 (Line 5). As this step is intended as a diversification mechanism, the neighborhood structure is based on a simple, but for the HEV-TSPTW very destructive, operator. We apply a random 2-Opt with Mode Change (Doppstadt et al., 2016) to the best local Solution. This means that we randomly delete two arcs from the solution and reconnect the partial tours with two new arcs, each chosen randomly. This leads to two possible new round tours, depending on the direction to start; one of these solutions is again randomly selected. This random 2-Opt with Mode Change is now repeated k-times. A 2-Opt move destroys a solution for a problem with time windows particularly heavily, as the tour is partially reversed. Figure 7 shows an example of applying the random 2-Opt with Mode Change to a solution with six customers twice. The example illustrates how destructive this operation actually can be, as only one out of seven arcs remains untouched compared to the original solution. 2. Improve local Solution0 by applying a hill-climbing Insert-Move to it and store the result as local Solution00 (Line 6). This step is comparable to 16

the hill-climbing Insert-Move during the initialization phase, apart from the fact that the charging of the resulting tour has to be considered and checked during the improvement phase, as the tour is not driven in pure combustion mode anymore. 3. Apply a hill-climbing χ-MC on local Solution00 and save the result as local Solution000 (Line 7). Thereby, we increase the likelihood to use arcs in a cheaper mode of operation and utilize the advantages of the HEV. This step is also comparable to the χ-MC during the initialization. After completing these steps, we compare the result of local Solution000 with the best local Solution (Line 8), assign local Solution000 to best local Solution and reset k = 1 if an improvement was found (Line 9/10). If no improvement was found, k is increased by 1 (Line 12). Note that during this iterative process i is only increased if no new best local Solution has been found in any neighborhood k ∈ K (Line 15). The algorithm terminates when the maximum number of iterations I is reached and returns the best local Solution found. The subroutine call chooseSolutionAorB used in (Line 8) compares two solutions by applying several metrics to determine the better solution. As this function is highly relevant for the whole process, we also describe it detailed in Algorithm 2. A feasible solution is always better than an infeasible one (Line 2). If both solutions are feasible (Line 4), the solution with the lower total costs is considered as the better one (Line 5). If the costs are actually equal (Line 7), we compare the final charging of the battery (Line 8), as a battery with a higher charging level increases the potential for further improvement. For two infeasible solutions (line 12), the one with the lower amount of total time window violation (Line 13), as sum of all violations of the tour, is regarded as the better one. If the time window violations are also equal (Line 15), we compare the final charging of the battery again (Line 16). The result of each VNS subprocess is finally evaluated within the P-VNS component, introduced in Figure 3, using the same chooseSolutionAorB func-

17

Algorithm 2 Comparison of Solutions 1: function chooseSolutionAorB(Solution A, Solution B) 2: 3: 4: 5:

if A is feasible and B is infeasible then return true else if A is feasible and B is feasible then if Costs of A < Costs of B then return true

6: 7:

else if Costs of A = Costs of B then if Final Charging of A > Final Charging of B then

8:

return true

9:

end if

10: 11: 12: 13:

end if else if A is infeasible and B is infeasible then if Time Window Violation of A < Time Window Violation of B then return true

14: 15:

else if Time Window Violation of A = Time Window Violation of B then if Final Charging of A > Final Charging of B then

16:

return true

17:

end if

18: 19:

end if

20:

end if

21:

return false

22:

end function

18

tion to determine the best global Solution, which is then returned as the final result. 4. Numerical Studies To evaluate our heuristic, we provide numerical results of different analyses in the following five subsections. First, we give detailed information on how the benchmark instances were created. Subsection 2 shows the exact results for the small instances with up to 10 customers derived with CPLEX. Afterwards, we provide different results produced with our heuristic, covering parameter tuning, different parameter setups and their impact on objective value and runtime, related to instance characteristics. In Subsection 4 we evaluate the performance of the heuristic by comparing its results to CPLEX and another heuristic. Finally, we generate pure combustion mode and ’arc splitting’ solutions to interpret the impact of HEVs on potential savings. 4.1. Benchmark Instances To create a set of test instances for the HEVTSPTW, we took the existing set of benchmark instances without time windows (Doppstadt et al., 2016), which are based on real world delivery tours. These 36 base instances differ in the number of customers (8, 10, 20, and 50) and the location of the depot. In the first subset, the depot is adjoining the delivery area; in the second subset the depot is 28km and in the third subset 57km apart from the delivery area. The delivery area itself has in all cases a coverage of approximately 25km2 . For each customer and depot location there are three different instances, varying service times and the exact customer locations within the delivery area. Therefore, the costs and travel times of these instances differ. For additional information on the original instances, we refer to the article mentioned above. We extended this existing set of instances by adding time windows for the depot – replacing the maximum tour duration of 8 hours – and for the customers. For each of the 36 existing instances, we added six different time window densities, resulting in a total of 216 instances. An overview of the base instance 19

TW 0 TW 1 TW 2 TW 3 TW 4 TW 5 Instance Cust. TWs Dist. Cust. TWs Dist. Cust. TWs Dist. Cust. TWs Dist. Cust. TWs Dist. Cust. TWs Dist. HEVTSPTW 1 08 1 8 0 0 8 2 0 8 6 0 8 2+ 0 8 4+ 0 8 8+ 0 HEVTSPTW 1 08 2 8 0 0 8 2 0 8 6 0 8 2+ 0 8 4+ 0 8 8+ 0 HEVTSPTW 1 08 3 8 0 0 8 2 0 8 6 0 8 2+ 0 8 4+ 0 8 8+ 0 HEVTSPTW 2 08 1 8 0 28 8 2 28 8 6 28 8 2+ 28 8 4+ 28 8 8+ 28 HEVTSPTW 2 08 2 8 0 28 8 2 28 8 6 28 8 2+ 28 8 4+ 28 8 8+ 28 HEVTSPTW 2 08 3 8 0 28 8 2 28 8 6 28 8 2+ 28 8 4+ 28 8 8+ 28 HEVTSPTW 3 08 1 8 0 57 8 2 57 8 6 57 8 2+ 57 8 4+ 57 8 8+ 57 HEVTSPTW 3 08 2 8 0 57 8 2 57 8 6 57 8 2+ 57 8 4+ 57 8 8+ 57 HEVTSPTW 3 08 3 8 0 57 8 2 57 8 6 57 8 2+ 57 8 4+ 57 8 8+ 57 HEVTSPTW 1 10 1 10 0 0 10 3 0 10 8 0 10 2+ 0 10 5+ 0 10 10+ 0 HEVTSPTW 1 10 2 10 0 0 10 3 0 10 8 0 10 2+ 0 10 5+ 0 10 10+ 0 HEVTSPTW 1 10 3 10 0 0 10 3 0 10 8 0 10 2+ 0 10 5+ 0 10 10+ 0 HEVTSPTW 2 10 1 10 0 28 10 3 28 10 8 28 10 2+ 28 10 5+ 28 10 10+ 28 HEVTSPTW 2 10 2 10 0 28 10 3 28 10 8 28 10 2+ 28 10 5+ 28 10 10+ 28 HEVTSPTW 2 10 3 10 0 28 10 3 28 10 8 28 10 2+ 28 10 5+ 28 10 10+ 28 HEVTSPTW 3 10 1 10 0 57 10 3 57 10 8 57 10 2+ 57 10 5+ 57 10 10+ 57 HEVTSPTW 3 10 2 10 0 57 10 3 57 10 8 57 10 2+ 57 10 5+ 57 10 10+ 57 HEVTSPTW 3 10 3 10 0 57 10 3 57 10 8 57 10 2+ 57 10 5+ 57 10 10+ 57 HEVTSPTW 1 20 1 20 0 0 20 8 0 20 18 0 20 2+ 0 20 10+ 0 20 20+ 0 HEVTSPTW 1 20 2 20 0 0 20 8 0 20 18 0 20 2+ 0 20 10+ 0 20 20+ 0 HEVTSPTW 1 20 3 20 0 0 20 8 0 20 18 0 20 2+ 0 20 10+ 0 20 20+ 0 HEVTSPTW 2 20 1 20 0 28 20 8 28 20 18 28 20 2+ 28 20 10+ 28 20 20+ 28 HEVTSPTW 2 20 2 20 0 28 20 8 28 20 18 28 20 2+ 28 20 10+ 28 20 20+ 28 HEVTSPTW 2 20 3 20 0 28 20 8 28 20 18 28 20 2+ 28 20 10+ 28 20 20+ 28 HEVTSPTW 3 20 1 20 0 57 20 8 57 20 18 57 20 2+ 57 20 10+ 57 20 20+ 57 HEVTSPTW 3 20 2 20 0 57 20 8 57 20 18 57 20 2+ 57 20 10+ 57 20 20+ 57 HEVTSPTW 3 20 3 20 0 57 20 8 57 20 18 57 20 2+ 57 20 10+ 57 20 20+ 57 HEVTSPTW 1 50 1 50 0 0 50 23 0 50 48 0 50 2+ 0 50 25+ 0 50 50+ 0 HEVTSPTW 1 50 2 50 0 0 50 23 0 50 48 0 50 2+ 0 50 25+ 0 50 50+ 0 HEVTSPTW 1 50 3 50 0 0 50 23 0 50 48 0 50 2+ 0 50 25+ 0 50 50+ 0 HEVTSPTW 2 50 1 50 0 28 50 23 28 50 48 28 50 2+ 28 50 25+ 28 50 50+ 28 HEVTSPTW 2 50 2 50 0 28 50 23 28 50 48 28 50 2+ 28 50 25+ 28 50 50+ 28 HEVTSPTW 2 50 3 50 0 28 50 23 28 50 48 28 50 2+ 28 50 25+ 28 50 50+ 28 HEVTSPTW 3 50 1 50 0 57 50 23 57 50 48 57 50 2+ 57 50 25+ 57 50 50+ 57 HEVTSPTW 3 50 2 50 0 57 50 23 57 50 48 57 50 2+ 57 50 25+ 57 50 50+ 57 HEVTSPTW 3 50 3 50 0 57 50 23 57 50 48 57 50 2+ 57 50 25+ 57 50 50+ 57

Table 1: Overview of the 216 Benchmark Instances’ Characteristics

data is given in Table 1. For each base instance (HEVTSPTW 1 08 1 - HEVTSPTW 3 50 3) and time window density (TW 0 - TW 5) combination, the table provides the number of customers (Cust.), the number of time windows (TWs) and the distance from the depot to the delivery area (Dist.). The number of TWs for each time window density are calculated as follows: For TW 0 |T W s| = 0, for TW 1 |T W s| = 0.5 ∗ |Cust.| − 2, for TW 2 |T W s| = |Cust.| − 2, for TW 3 |T W s| = 2, for TW 4 |T W s| = 0.5 ∗ |Cust.| and for TW 5 |T W s| = |Cust.|. The + for the TWs in TW 3 - TW 5 indicates that these instances contain a special pair of customers with very late time windows. This combination will, with a high probability, cause waiting time and might require the boost mode between these two customers or on the way back to the depot. For a single base instance the customers and their TWs are always identical, meaning that e.g., HEVTSPTW 1 08 1 TW 5 contains the same two special TWs as

20

HEVTSPTW 1 08 1 TW 3. Additionally, the regular TWs are also identical for each customer per instance (HEVTSPTW 1 08 1 TW 2 contains the TWs from HEVTSPTW 1 08 1 TW 1). Therefore, the TWs from TW 2 and TW 3 sum up to the TWs of TW 5 and TW 5 is the superset of TWs for any other time window density setting for the same base instance. The depot is neither contained in |Cust.| nor in |T W s| and the starting time at the depot is always set to 0 and the end time to 480 minutes after start, which is equal to a maximum duration of 8 hours. As in Savelsbergh (1985), we randomly generate time windows for a subset of customers. In our procedure, first we generate the time windows for the TW 5 instance and then randomly remove customers to create classes with fewer time windows. We set the time window width to 30 minutes for all instances and ensure for the two late time windows that the distance between the associated customers and the depot can be driven in boost mode. As the time windows are rounded to full integer numbers, it is potentially still possible to use other modes of operation. The remaining time windows are generated with a uniform distribution over an interval that can be reached from the depot in a feasible solution. As this procedure does not guarantee that the instance has a feasible solution, we run our ’Initialization Phase’ with the last two arcs fixed to boost mode and check the feasibility. If the instance is infeasible, the time window generation process is restarted. All 216 instances are available as dataset for downloading on Mendeley Data (Doppstadt et al., 2019). 4.2. Results of Exact Solution Approach All numerical results for the exact and heuristic computations were generated on the same workstation, running with an Intel Core i7-4790 3.6 GHz and 24 GB Ram on Windows 10 having only system processes running in parallel to our tests. CPLEX version is 12.6. To solve the problem with CPLEX we limited the maximum calculation time to 100 hours per instance. With this limitation, we were able to find exact 21

Instance Objective Dist. Cust. TWs. Var. 0 8 0 1 1830.69 0 8 2 1 1849.37 0 8 6 1 1849.37 0 8 2+ 1 1849.37 0 8 4+ 1 1849.37 0 8 8+ 1 1849.37 0 8 0 2 1553.15 0 8 2 2 1553.15 0 8 6 2 1553.15 0 8 2+ 2 1553.15 0 8 4+ 2 1553.15 0 8 8+ 2 1553.15 0 8 0 3 1435.74 0 8 2 3 1436.09 0 8 6 3 1460.31 0 8 2+ 3 1436.09 0 8 4+ 3 1436.09 0 8 8+ 3 1460.31

Time 188 200 23 7 7 1 586 567 139 40 23 7 1938 2971 147 104 55 1

Instance Objective Dist. Cust. TWs. Var. 28 8 0 1 7189.85 28 8 2 1 7189.85 28 8 6 1 7189.85 28 8 2+ 1 7202.07 28 8 4+ 1 7202.07 28 8 8+ 1 7202.07 28 8 0 2 7140.95 28 8 2 2 7174.85 28 8 6 2 7174.85 28 8 2+ 2 7174.85 28 8 4+ 2 7174.85 28 8 8+ 2 7174.85 28 8 0 3 7292.25 28 8 2 3 7292.25 28 8 6 3 7292.25 28 8 2+ 3 7292.25 28 8 4+ 3 7292.25 28 8 8+ 3 7292.25

Time 395 283 5 382 69 1 481 355 8 507 63 1 86 101 10 88 35 1

Instance Objective Dist. Cust. TWs. Var. 57 8 0 1 12706.60 57 8 2 1 12706.60 57 8 6 1 12724.80 57 8 2+ 1 12706.60 57 8 4+ 1 12706.60 57 8 8+ 1 12724.80 57 8 0 2 12687.90 57 8 2 2 12713.00 57 8 6 2 12713.00 57 8 2+ 2 12713.00 57 8 4+ 2 12713.00 57 8 8+ 2 12713.00 57 8 0 3 12708.60 57 8 2 3 12708.60 57 8 6 3 12738.00 57 8 2+ 3 12738.00 57 8 4+ 3 12738.00 57 8 8+ 3 12738.00

Time 541 376 22 752 103 1 596 849 26 807 206 1 741 783 35 1283 247 7

Table 2: CPLEX Results of 8 Customer Instances

solutions for all instances with 8 customers and for almost all instances with 10 customers. Note that due to the required node duplication of the depot, the instances with 8 customers have actually 10 nodes for the model formulation and 12 nodes for the instances with 10 customers, respectively. In Table 2 we provide the CPLEX results for the instances with 8 customers. Beside the instance characteristics (distance to depot, number of customers, time window class and customer locations variant) and the objective value found, which is the optimum for all 8 customer instances, we also provide the required solution time rounded up to full seconds. The results for the 10 customer instances are given in Table 3. The instances where CPLEX was not able to find the optimal solution are marked with an asterisk (*). For these five instances the objective value is the lower bound found until the time limit was reached. The results indicate that, in general, the more restricted the solution space is by time windows, the easier the instances are to solve for CPLEX. Especially, the customer pair with the late time windows is decreasing the required calculation time for most of the instances (TW 3 - TW 5 compared to TW 0 - TW 2). However, this statement is not true for all of the instances, e.g., solving HEVTSPTW 3 10 1 TW 1 took significantly longer than HEVTSPTW 3 10 1 TW 0, 22

Instance Instance Instance Objective Time Objective Time Objective Time Dist. Cust. TWs. Var. Dist. Cust. TWs. Var. Dist. Cust. TWs. Var. 0 10 0 1 1798.94 37464 28 10 0 1 7308.15 252926 57 10 0 1 12747.20 288003 0 10 3 1 1798.94 49111 28 10 3 1 7314.32 81083 57 10 3 1 12759.20 322796 0 10 8 1 1801.46 1468 28 10 8 1 7329.18 1039 57 10 8 1 12759.20 687 0 10 2+ 1 1798.94 508 28 10 2+ 1 6944.72 * 360000 57 10 2+ 1 12759.20 251960 0 10 5+ 1 1798.94 382 28 10 5+ 1 7314.32 6559 57 10 5+ 1 12759.20 8192 0 10 10+ 1 1801.46 11 28 10 10+ 1 7329.18 67 57 10 10+ 1 12759.20 12 0 10 0 2 1203.25 * 360002 28 10 0 2 7362.93 44797 57 10 0 2 12489.00 * 360002 0 10 3 2 1601.73 82567 28 10 3 2 7362.93 20729 57 10 3 2 12738.40 321980 0 10 8 2 1601.73 508 28 10 8 2 7362.93 1076 57 10 8 2 12738.40 626 0 10 2+ 2 1601.73 3772 28 10 2+ 2 7362.93 59937 57 10 2+ 2 12445.20 * 360003 0 10 5+ 2 1601.73 317 28 10 5+ 2 7362.93 6720 57 10 5+ 2 12738.40 26316 0 10 10+ 2 1601.73 1 28 10 10+ 2 7362.93 141 57 10 10+ 2 12738.40 33 0 10 0 3 1478.90 345616 28 10 0 3 7290.82 11869 57 10 0 3 12935.50 217161 0 10 3 3 1531.45 312300 28 10 3 3 7290.82 16341 57 10 3 3 12935.50 343006 0 10 8 3 1531.45 143 28 10 8 3 7290.82 953 57 10 8 3 12935.50 4844 0 10 2+ 3 1531.45 29565 28 10 2+ 3 7290.82 15772 57 10 2+ 3 12648.60 * 360000 0 10 5+ 3 1531.45 19094 28 10 5+ 3 7290.82 1403 57 10 5+ 3 12935.50 14576 0 10 10+ 3 1531.45 48 28 10 10+ 3 7290.82 5 57 10 10+ 3 12935.50 80

Table 3: CPLEX Results of 10 Customer Instances

which has 2 instead of 0 TWs. At this point we would also like to point out the massive increase of overall runtime comparing the 8 and 10 customer instances. Solving the whole set of 8 customer instances took approximately 4.8 hours. In contrast, solving the set of 10 customer instances took almost 60 days, even without finding the optimal solution for all instances. The results already illustrate the reason why a heuristic is required to solve instances within a reasonable calculation time and to provide a feasible solution at all. 4.3. Results of Heuristic Solution Approach In this subsection, we apply different heuristic solution approaches. As a first step, we perform tests to tune the parameters of our heuristic. Based on the parameter tuning, we run two different setups for all instances with our HEV-TSPTW heuristic. Parameter Tuning The goal for the heuristic is to achieve a solution quality close to the exact solutions, while requiring only a fraction of the calculation time. To approach this goal, we need to define several parameters for the heuristic setup. The first

23

Instance Dist. Cust. TWs. Var. 0 8 6 3 28 8 8+ 2 57 8 2 1 0 10 0 3 28 10 2+ 2 57 10 3 2 0 20 20+ 3 28 20 0 3 57 20 10+ 2 0 50 50+ 2 28 50 0 3 57 50 2+ 2

3-MC Objective 1513.43 7196.07 12815.27 1725.81 7454.36 12885.34 1611.32 8137.58 13393.22 2499.34 9117.61 14260.57

Time 0.0 0.0 0.0 0.2 0.1 0.1 2.9 2.7 3.0 136.6 108.3 105.8

Objective 1492.94 7196.07 12750.25 1729.79 7375.8 12812.55 1610.3 8141.16 13319.64 2497.82 9094.77 14261.58

4-MC Gap 3-MC -1.35% 0.00% -0.51% 0.23% -1.05% -0.56% -0.06% 0.04% -0.55% -0.06% -0.25% 0.01%

Time 2.4 1.9 3.0 17.2 9.2 8.9 329.0 375.1 672.5 30578.3 55704.2 45328.1

Objective 1489.51 7176.51 12747.8 1723.32 7378.2 12814.16 1610.3 8141.95 13322.65 2497.82 9095.41 14262.24

5-MC Gap 3-MC -1.58% -0.27% -0.53% -0.14% -1.02% -0.55% -0.06% 0.05% -0.53% -0.06% -0.24% 0.01%

Time 4.9 5.0 4.2 68.5 30.1 25.2 1953.7 1579.9 3718.1 190276.0 624888.8 1051873.5

Table 4: Comparison of different χ-Mode Changes for Initializations Phase

parameter to set is χ for the χ-Mode operator during the initialization phase. To avoid an over-tuning of the heuristic, we limit the set for the tuning to 12 out of the 216 instances (5.56 %). We ran the initialization process with χ = 3, χ = 4 and χ = 5. The results are given in Table 4. For each instance, we provide the objective value, the runtime in seconds and for the 4-MC and 5-MC we also provide the percentage objective value gap to the 3-MC. All generated solutions for the 12 instances are valid. In terms of calculation time, the results show that the time from 3-MC to 4-MC increases up to a factor of over 500 and another factor of up to over 20 from 4-MC to 5-MC. In contrast, the maximum gained solution improvement for 4-MC and 5-MC compared to 3-MC is below 1.6 %. This led us to the decision to chose the 3-MC not only during the initialization, but also during the improvement phase. To define the parameter setup for the improvement phase itself, we run the same 12 instances as before with a total of 48 different parameter combinations, defined by four different parameters. The first parameter is the number of threads (Th.) to be executed in parallel. As the CPU of our system has four physical cores, we selected four and eight as potential values for this parameter. Next, we had to define the number of iterations (Itr.) running within each thread. As the overall number of iterations is based on threads multiplied by iterations per threads, we used 25, 50 and 100 as values here, to allow a

24

Param. Setup Th. Itr. K. Im. 4 25 4 BI 8 100 16 FI Param. Setup Th. Itr. K. Im. 4 25 4 BI 8 100 16 FI Param. Setup Th. Itr. K. Im. 4 25 4 BI 8 100 16 FI

Dist. 0 Cust. 8 TWs. 6 Var. 3 Dist. 28 Cust. 8 TWs. 8+ Var. 2 Avg. Agg. Std. Avg. Agg. Std. Obj. Time Dev. Obj. Time Dev. 1460.31 0.6 0.00 7174.85 0.4 0.00 1460.31 9.4 0.00 7174.85 10.6 0.00

Dist. 57 Cust. 8 TWs. 2 Var. 1 Avg. Agg. Std. Obj. Time Dev. 12708.00 0.5 2.79 12706.60 11.6 0.00

Dist. 0 Cust. 10 TWs. 0 Var. 3 Avg. Agg. Std. Obj. Time Dev. 1489.41 1.2 21.02 1478.90 22.4 0.00

Dist. 28 Cust. 10 TWs. 2+ Var. 2 Dist. 57 Cust. 10 TWs. 3 Var. 2 Dist. 0 Cust. 20 TWs. 20+ Var. 3 Dist. 28 Cust. 20 TWs. 0 Var. 3 Avg. Agg. Std. Avg. Agg. Std. Avg. Agg. Std. Avg. Agg. Std. Obj. Time Dev. Obj. Time Dev. Obj. Time Dev. Obj. Time Dev. 7362.93 1.0 0.00 12738.38 0.9 0.00 1608.37 17.6 0.00 7717.78 16.7 3.16 7362.93 24.2 0.00 12738.38 28.4 0.00 1608.37 530.1 0.00 7712.54 385.0 1.70 Dist. 57 Cust. 20 TWs. 10+ Var. 2 Dist. 0 Cust. 50 TWs. 50+ Var. 2 Dist. 28 Cust. 50 TWs. 0 Var. 3 Dist. 57 Cust. 50 TWs. 2+ Var. 2 Avg. Agg. Std. Avg. Agg. Std. Avg. Agg. Std. Avg. Agg. Std. Obj. Time Dev. Obj. Time Dev. Obj. Time Dev. Obj. Time Dev. 13287.61 17.9 0.00 2497.78 820.1 0.00 8749.51 761.4 50.57 14288.27 749.0 93.94 13287.61 796.7 0.00 2497.78 37670.2 0.00 8505.83 12939.3 39.53 14041.25 18379.0 33.65

Table 5: Comparison of different Parameter Setups for Improvement Phase

comparison of the two different thread settings with an equal number of total iterations. E.g., 8 threads with 25 iterations have the same total number of iterations as 4 threads with 50 iterations. As third parameter we use the VNS specific K-value (K), which defines how many neighborhoods to search within each iteration. Basically, for our P-VNS the higher the K-value the more the previous solution is destructed before going back into the improvement phase at the end of each iteration step. As last parameter, we define an improvement strategy (Im.), distinguishing between first improvement (FI) and best improvement (BI). This parameter affects all improvement processes within the heuristic and determines if after finding a first improved solution the current improvement iteration is terminated (FI) or if the whole neighborhood is explored and the solution with the largest/best improvement is finally selected (BI). Due to the randomness factor within the P-VNS, we ran the test five times for each instance with every setup of the 48 combinations. Two setups showed an outstanding behavior compared to the other ones. The results for these parameter combinations and aggregated run information are given in Table 5, the full results can be found in the supplementary material. We provide the average objective value, the aggregated runtime in seconds and the standard deviation of the objective value over the five runs. Unsurprisingly, the overall results show that with an increase of the instance size, the number of iterations

25

also needs to be increased to find the overall best solutions and therefore, more calculation time is required. The two selected parameter setups reflect that our focus is on both solution quality and runtime. With 4 Th., 25 Itr., 4 K and BI the P-VNS had the shortest runtime for seven out of the 12 instances, while finding the best solution. Subsequently, we refer to this setup as ’runtime optimized’. On the other hand, 8 Th., 100 Itr., 18 K and FI was the only setup to find the best solution for all of the 12 instances. We refer to this as ’quality optimized’ setup. Impact of the Number of Customers In the next step, we ran the heuristic with the runtime optimized and the quality optimized parameter setup on all 216 instances five times and calculated overall best and average solutions per instance as well as the aggregated runtime for the five runs. Based on these results, we clustered the instance set by number of customers, to evaluate the impact on the solution quality. Figure 8 shows two boxplots with the percentage gap between best and average solution on the y-axis and the number of customers on the x-axis. Part (a) are the results for the runtime optimized setup whereas Part (b) are the results for the quality optimized setup. The figure illustrates that overall solution quality is already very stable for the runtime optimized setup but unsurprisingly becomes even more stable for the quality optimized runs, except for a small number of outliers. However, in both cases the deviation between best and average solutions becomes noticeably higher for the 50 customer instances than for the other groups even if it is below 1 % for the majority of the instances. In Figure 9 we provide the aggregated runtime for the five runs in seconds on the y-axis and the same clustering as before. In contrast to Figure 8 where the scales are identical, the scale for the runtime on the y-axis is 50 times larger for the quality optimized setup in Part (b) compared to Part (a). It becomes obvious, that the required runtime for the 8, 10 and 20 customer instances is insignificant compared to the instances with 50 customers. The total runtime for the runtime optimized results was around 12.8 hours, whereas, the overall 26

5

4

4 % Gap Best / Avg. Solution

% Gap Best / Avg. Solution

5

3

2

1

3

2

1

0

8

10 20 Number of Customers

0

50

8

(a) Runtime Optimized Results

10 20 Number of Customers

50

(b) Quality Optimized Results

Figure 8: Customer Clustered Gap of Best and Avg. Objective Value 1200

60000 50000 Aggregated Runtime (sec.)

Aggregated Runtime (sec.)

1000 800 600 400 200

40000 30000 20000 10000

0

0

8

10 20 Number of Customers

50

8

(a) Runtime Optimized Results

10 20 Number of Customers

50

(b) Quality Optimized Results

Figure 9: Customer Clustered Runtime

runtime for quality optimized tests reached 429.8 hours, which is almost 18 days. Impact of Depot to Delivery Area Distance Beside the impact of the number of customers on quality and runtime, we were also interested in the influence of the distance between delivery area and the

27

depot. For this analysis, we clustered the results by the three distance groups considering the subset of the 50 customer instances to avoid any side effects caused by this parameter. Figure 10 provides the two corresponding boxplots

5

5

4

4 % Gap Best / Avg. Solution

% Gap Best / Avg. Solution

for the gap between best and average solution. On first sight, it appears that

3

2

1

0

3

2

1

0

28 Distance to Depot (km)

0

57

(a) Runtime Optimized Results

0

28 Distance to Depot (km)

57

(b) Quality Optimized Results

Figure 10: 50 Customer Depot Distance Clustered Gap of Best and Avg. Objective Value

the instances where the depot is adjoining the delivery areas are more difficult to solve, as the deviation is higher for both parameter setups. Nevertheless, there is also another explanation. The overall objective values are higher for the instances where the depot is farther away from the delivery area as their solutions contain two significantly longer arcs (into and out of the delivery area). Therefore, a comparable absolute deviation in the objective value leads to a higher percentage gap for the instances where the distance is small. Figure 11 is also based on the data of the 50 customers only and provides the runtime on the y-axis. In our opinion, the distance to the depot does not have a significant impact on the runtime. The only noticeable finding is that the range of the runtime for the runtime optimized results on the last group (57 km) is smaller than for any other group. However, we could not find a reasonable explanation for this behavior and it might just be a coincident.

28

1200

60000 50000 Aggregated Runtime (sec.)

Aggregated Runtime (sec.)

1000

800

600

40000 30000 20000

400 10000

0

28 Distance to Depot (km)

57

0

(a) Runtime Optimized Results

28 Distance to Depot (km)

57

(b) Quality Optimized Results

Figure 11: 50 Customer Depot Distance Runtime

Impact of Time Window Groups In this section, we investigate the impact of the time window groups on quality and runtime. To this end, we cluster the results by the six categories of time windows. Limiting the scope to 50 customers instances becomes essential in this case to ensure that the groups have the exact number of time windows. Figure 12 provides the two boxplots for the solution quality analysis. The results indicate that the strict time window combination in group 2 and 5 makes the instances easier to solve. This is consistent with the results generated with CPLEX, where these instances required significantly lower runtimes. As for the heuristic the runtime is basically dependent on the parameter setup, the results do not show a lower runtime but a reduced deviation. Figure 13 provides the runtime results for the time window clustering. Similarly to Figure 11, we do not see an impact of the types of time windows on the runtimes. The full results for the heuristic runs can again be found in the supplementary material.

29

5

4

4 % Gap Best / Avg. Solution

% Gap Best / Avg. Solution

5

3

2

1

3

2

1

0

0

1

2 3 Time Window Group

4

0

5

(a) Runtime Optimized Results

0

1

2 3 Time Window Group

4

5

(b) Quality Optimized Results

Figure 12: 50 Customer Time Window Clustered Gap of Best and Avg. Objective Value 1200

60000 50000 Aggregated Runtime (sec.)

Aggregated Runtime (sec.)

1000

800

600

40000 30000 20000

400 10000

0

1

2 3 Time Window Group

4

5

0

(a) Runtime Optimized Results

1

2 3 Time Window Group

4

5

(b) Quality Optimized Results

Figure 13: 50 Customer Time Window Clustered Runtime

4.4. Performance Evaluation of P-VNS As it is impossible to evaluate the performance of a sole heuristic in terms of solution quality or runtime, we use two types of comparisons in this subsection to provide such a performance evaluation. First, we compare the results for the 8 and 10 customer instances generated by CPLEX with the P-VNS results. 30

Runtime Instance Optimized Dist. Cust. TWs. Var. Gap Factor 0 8 0 1 0.00% 0.0 0 8 2 1 0.00% 0.0 0 8 6 1 0.00% 0.0 0 8 2+ 1 0.00% 0.1 0 8 4+ 1 0.00% 0.1 0 8 8+ 1 0.00% 0.4 0 8 0 2 0.00% 0.0 0 8 2 2 0.00% 0.0 0 8 6 2 0.00% 0.0 0 8 2+ 2 0.00% 0.0 0 8 4+ 2 0.00% 0.0 0 8 8+ 2 0.00% 0.1 0 8 0 3 0.00% 0.0 0 8 2 3 0.68% 0.0 0 8 6 3 0.00% 0.0 0 8 2+ 3 0.00% 0.0 0 8 4+ 3 0.00% 0.0 0 8 8+ 3 0.00% 0.3 28 8 0 1 0.00% 0.0 28 8 2 1 0.00% 0.0 28 8 6 1 0.00% 0.1 28 8 2+ 1 0.00% 0.0 28 8 4+ 1 0.00% 0.0 28 8 8+ 1 0.00% 0.7 28 8 0 2 0.00% 0.0 28 8 2 2 0.00% 0.0 28 8 6 2 0.00% 0.0 28 8 2+ 2 0.00% 0.0 28 8 4+ 2 0.00% 0.0 28 8 8+ 2 0.00% 0.7 28 8 0 3 0.00% 0.0 28 8 2 3 0.00% 0.0 28 8 6 3 0.00% 0.0 28 8 2+ 3 0.00% 0.0 28 8 4+ 3 0.00% 0.0 28 8 8+ 3 0.00% 0.5 57 8 0 1 0.00% 0.0 57 8 2 1 0.00% 0.0 57 8 6 1 0.00% 0.0 57 8 2+ 1 0.00% 0.0 57 8 4+ 1 0.00% 0.0 57 8 8+ 1 0.00% 0.4 57 8 0 2 0.00% 0.0 57 8 2 2 0.00% 0.0 57 8 6 2 0.00% 0.0 57 8 2+ 2 0.00% 0.0 57 8 4+ 2 0.00% 0.0 57 8 8+ 2 0.00% 0.4 57 8 0 3 0.00% 0.0 57 8 2 3 0.00% 0.0 57 8 6 3 0.00% 0.0 57 8 2+ 3 0.18% 0.0 57 8 4+ 3 0.18% 0.0 57 8 8+ 3 0.00% 0.1

Quality Optimized Gap Factor 0.00% 0.0 0.00% 0.0 0.00% 0.4 0.00% 1.5 0.00% 1.5 0.00% 9.6 0.00% 0.0 0.00% 0.0 0.00% 0.1 0.00% 0.2 0.00% 0.4 0.00% 1.3 0.00% 0.0 0.00% 0.0 0.00% 0.1 0.00% 0.1 0.00% 0.2 0.00% 7.3 0.00% 0.0 0.00% 0.0 0.00% 2.2 0.00% 0.0 0.00% 0.1 0.00% 19.0 0.00% 0.0 0.00% 0.0 0.00% 1.3 0.00% 0.0 0.00% 0.1 0.00% 20.5 0.00% 0.1 0.00% 0.1 0.00% 0.9 0.00% 0.1 0.00% 0.3 0.00% 9.3 0.00% 0.0 0.00% 0.0 0.00% 0.5 0.00% 0.0 0.00% 0.1 0.00% 11.8 0.00% 0.0 0.00% 0.0 0.00% 0.4 0.00% 0.0 0.00% 0.1 0.00% 11.1 0.00% 0.0 0.00% 0.0 0.00% 0.3 0.00% 0.0 0.00% 0.0 0.00% 1.8

Runtime Quality Instance Optimized Optimized Dist. Cust. TWs. Var. Gap Factor Gap Factor 0 10 0 1 0.00% 0.0 0.00% 0.0 0 10 3 1 0.00% 0.0 0.00% 0.0 0 10 8 1 0.00% 0.0 0.00% 0.0 0 10 2+ 1 0.00% 0.0 0.14% 0.0 0 10 5+ 1 0.00% 0.0 0.14% 0.1 0 10 10+ 1 0.00% 0.1 0.00% 2.1 0 10 0 2 24.71% * 0.0 24.71% * 0.0 0 10 3 2 0.00% 0.0 0.02% 0.0 0 10 8 2 0.00% 0.0 0.00% 0.1 0 10 2+ 2 0.00% 0.0 0.02% 0.0 0 10 5+ 2 0.00% 0.0 0.02% 0.1 0 10 10+ 2 0.00% 1.1 0.00% 26.9 0 10 0 3 0.00% 0.0 0.00% 0.0 0 10 3 3 0.00% 0.0 0.00% 0.0 0 10 8 3 0.00% 0.0 0.00% 0.2 0 10 2+ 3 0.00% 0.0 0.00% 0.0 0 10 5+ 3 0.00% 0.0 0.00% 0.0 0 10 10+ 3 0.00% 0.0 0.00% 0.5 28 10 0 1 0.00% 0.0 0.00% 0.0 28 10 3 1 0.00% 0.0 0.00% 0.0 28 10 8 1 0.00% 0.0 0.00% 0.0 28 10 2+ 1 5.05% * 0.0 5.05% * 0.0 28 10 5+ 1 0.00% 0.0 0.00% 0.0 28 10 10+ 1 0.00% 0.0 0.00% 0.4 28 10 0 2 0.00% 0.0 0.00% 0.0 28 10 3 2 0.00% 0.0 0.00% 0.0 28 10 8 2 0.00% 0.0 0.00% 0.0 28 10 2+ 2 0.00% 0.0 0.00% 0.0 28 10 5+ 2 0.00% 0.0 0.00% 0.0 28 10 10+ 2 0.00% 0.0 0.00% 0.2 28 10 0 3 0.00% 0.0 0.00% 0.0 28 10 3 3 0.00% 0.0 0.00% 0.0 28 10 8 3 0.00% 0.0 0.00% 0.0 28 10 2+ 3 0.00% 0.0 0.00% 0.0 28 10 5+ 3 0.00% 0.0 0.00% 0.0 28 10 10+ 3 0.00% 0.2 0.00% 4.2 57 10 0 1 0.00% 0.0 0.00% 0.0 57 10 3 1 0.00% 0.0 0.00% 0.0 57 10 8 1 0.00% 0.0 0.00% 0.0 57 10 2+ 1 0.00% 0.0 0.00% 0.0 57 10 5+ 1 0.00% 0.0 0.00% 0.0 57 10 10+ 1 0.00% 0.1 0.00% 2.4 57 10 0 2 1.86% * 0.0 1.86% * 0.0 57 10 3 2 0.00% 0.0 0.00% 0.0 57 10 8 2 0.00% 0.0 0.00% 0.1 57 10 2+ 2 2.30% * 0.0 2.30% * 0.0 57 10 5+ 2 0.00% 0.0 0.00% 0.0 57 10 10+ 2 0.00% 0.0 0.00% 1.1 57 10 0 3 0.00% 0.0 0.00% 0.0 57 10 3 3 0.00% 0.0 0.00% 0.0 57 10 8 3 0.00% 0.0 0.00% 0.0 57 10 2+ 3 2.22% * 0.0 2.22% * 0.0 57 10 5+ 3 0.00% 0.0 0.00% 0.0 57 10 10+ 3 0.00% 0.0 0.00% 0.4

Table 6: Comparison of CPLEX with P-VNS Results

Second, we use the results from Doppstadt et al. (2016) to compare the results for the instances without time windows to provide at least a partial comparison between different types of heuristics. Comparison of CPLEX and P-VNS For the comparison between the CPLEX results and P-VNS we use both, runtime optimized and quality optimized setup P-VNS results. In Table 6 we provide these values in an aggregated format. The ’Gap’-columns are the per-

31

centage gaps between the CPLEX and the two P-VNS objective values. All instances unsolved by CPLEX are marked with an asterisk. In these cases the gap is measured as the difference between the lower bound found by CPLEX after 100 hours and the P-VNS solution. ’Factor’ is the time factor comparing the CPLEX and the aggregated P-VNS runtimes for the five runs. For three of the 8 customer instances the runtime optimized P-VNS did not find the optimal solution, for five of the 10 customer instances the quality optimized setup had a small gap compared to the optimal CPLEX solution. Overall, the PVNS was able to find all solutions proven to be optimal by CPLEX. For the instances not solved by CPLEX a final statement regarding the solution quality is not possible. However, during all numerical studies with any used setup and configuration, we never found any better solution with our heuristic. In terms of runtimes the results provide some interesting findings. In general, the runtimes for the heuristic, no matter if runtime or quality optimized setup, are insignificant compared to the CPLEX runtimes. For some of the instances with the highest time window density (TW 5) this changes slightly. The runtime optimized P-VNS is still faster in all but one case (Factor 1.1), whereas, the quality optimized setup requires up to almost 27 times as long as CPLEX. The reason for this is basically that CPLEX seems to have an advantage whenever an instance is highly restricted and the runtime massively drops in these cases while the runtimes for the P-VNS are relatively constant. Note that the CPLEX solutions have some small precision error in the objective value. This is the case when the objective values have five integer digits and affects the second position after decimal point. This means that for the same solution CPLEX and the P-VNS have a small deviation in the objective value. Manual checks have shown that the objective value generated with the P-VNS is the exact one. Evaluation of Traveling Salesman Instances without Time Windows In this subsection we compare partial results of the P-VNS with ITS results for HEV-TSP to get an additional evaluation of the P-VNS solution quality. 32

Runtime Instance Optimized Dist. Cust. TWs. Var. Gap Factor 0 8 0 1 0.00% 0.0 0 8 0 2 0.00% 0.0 0 8 0 3 0.00% 0.0 28 8 0 1 0.00% 0.0 28 8 0 2 0.00% 0.0 28 8 0 3 0.00% 0.0 57 8 0 1 0.00% 0.0 57 8 0 2 0.00% 0.0 57 8 0 3 0.00% 0.0 0 10 0 1 0.00% 0.0 0 10 0 2 0.00% 0.0 0 10 0 3 0.00% 0.0 28 10 0 1 0.00% 0.0 28 10 0 2 0.00% 0.0 28 10 0 3 0.00% 0.0 57 10 0 1 0.00% 0.0 57 10 0 2 0.00% 0.0 57 10 0 3 0.00% 0.0

Quality Optimized Gap Factor 0.00% 0.4 0.00% 0.4 0.00% 0.4 0.00% 0.4 0.00% 0.4 0.00% 0.4 0.00% 0.5 0.00% 0.4 0.00% 0.5 0.00% 0.4 0.00% 0.4 0.00% 0.4 0.00% 0.3 0.00% 0.4 0.00% 0.3 0.00% 0.4 0.00% 0.4 0.00% 0.4

Instance Dist. Cust. TWs. Var. 0 20 0 1 0 20 0 2 0 20 0 3 28 20 0 1 28 20 0 2 28 20 0 3 57 20 0 1 57 20 0 2 57 20 0 3 0 50 0 1 0 50 0 2 0 50 0 3 28 50 0 1 28 50 0 2 28 50 0 3 57 50 0 1 57 50 0 2 57 50 0 3

Runtime Quality Optimized Optimized Gap Factor Gap Factor 0.00% 0.0 0.00% 1.1 0.00% 0.1 0.00% 1.4 0.00% 0.1 0.01% 1.0 0.00% 0.0 0.00% 0.8 0.02% 0.0 0.00% 0.9 0.11% 0.0 0.06% 0.8 0.06% 0.0 -0.05% 1.0 0.02% 0.0 0.00% 1.1 0.00% 0.0 0.00% 1.2 9.16% 0.1 0.17% 2.3 8.66% 0.1 3.23% 2.2 3.31% 0.1 -0.11% 2.3 0.70% 0.1 0.28% 0.7 2.90% 0.1 1.69% 1.4 3.47% 0.1 0.73% 1.1 1.01% 0.1 -0.61% 1.4 1.29% 0.1 0.64% 1.5 0.89% 0.1 0.44% 1.7

Table 7: Comparison of Iterated Tabu Search with P-VNS for Instances without Time Windows

As the results from the ITS paper are only available for instances without time windows, the scope of this experiment is limited to these instances (TW 0). Both tests were executed on the same hardware, so the runtimes are also comparable. Again, we use the runtime and the quality optimized setup. The results of this comparison are pictured in Table 7 in the same compact format as for the previous test. During these tests, we were able to find three new best solutions for the HEV-TSP with the P-VNS; one for a 20 customer instance and two for 50 customer instances. The solutions for ITS and P-VNS are equal for the 8 and 10 customer instances with significantly lower runtimes for the runtime optimized setup and are still between 50 and 70 % lower runtimes for the quality optimized setup. For the 20 customer instances the objective values of the PVNS are akin to the ITS, whereas, for the 50 customer instances, except for the new best solutions, the ITS performs better than the P-VNS, especially, compared to the runtime optimized setup. For the 20 and 50 customer instances the runtimes for the runtime optimized runs are still very low, for the quality optimized runs the results in terms of required runtime are inconclusive and the time factor range is between 0.7 and 2.3. Considering the new best solutions and the relatively lower runtime, the performance of the P-VNS seems to be reasonably good, particularly, as it has been tailored for the HEV-TSPTW and not the HEV-TSP. 33

4.5. Evaluation of Hybrid Electric Vehicles and Intra Arc Mode Change Finally, we are looking into more operational aspects of the HEV-TSPTW. In this subsection, we initially generate solutions in pure combustion mode. Next, we use a different method to partially drive the route in electric mode, by splitting an arc and changing the mode of operation at this point. In the last experiment, we combine the P-VNS with this arc splitting method. As the quality gap between the ’runtime optimized’ and the ’quality optimized’ setup is relatively low, we decided to proceed with the ’runtime optimized’ setup for the consecutive tests. Pure Combustion Vehicles Results Beside the actual results of HEV-TSPTW, one of the most interesting aspects of HEVs is the potential savings compared to pure combustion vehicles. To be able to provide such a comparison, we conduct an additional numerical experiment by limiting our P-VNS heuristic to pure combustion mode. This means that we actually solve our instances as TSPTW. This allows us to compare hybrid and combustion vehicle solutions and estimate the cost savings based on the instances. To achieve this, the P-VNS was modified in the following sense. Looking back at Figure 3, the Shake move is performed without Mode Change and the Hill Climbing χ-Mode Change is simply skipped. Therefore, the mode of an arc based on the pure combustion solution from the initialization phase is never changed. At this point we admit that there are potentially better heuristics available for the TSPTW as shown in the introduction, however, as this is a concurrent experiment and we just want to point out potential saving, this approach is providing sufficiently good solution quality. To point out the most interesting finding of this test, we clustered the data again by distance between depot and delivery area. The boxplot in Figure 14 shows the percentage savings of the P-VNS (runtime optimized) solution compared to the pure combustion vehicle tests. Results indicate potential savings of up to 13 %; especially, the closer the delivery area is to the depot, the higher are the savings. This can be explained by long arcs into the delivery area, which are always driven in com34

% Combustion Solution / Runtime Optimized P-VNS

12 10 8 6 4 2 0 0

28 Distance to Depot (km)

57

Figure 14: 50 Customer Combustion Vehicles Results

bustion or charging mode and therefore, contribute significantly to the costs, but cannot generate reasonable savings. This means that the potential savings of the HEV mostly depend on the relation between the distance to the depot and within the delivery area. In contrast to the delivery area, the time windows seem to have no significant impact on the savings. The experiment took less than 1.8 hours for all instances and five runs, strengthening the earlier statement on the complexity increase due to the different modes of operation. Full results can be found in the supplementary material. Combustion Solution plus Intra Arc Mode Change Results A limitation of our approach is that the mode of operation of the vehicle can not be changed intra-arc. (Doppstadt et al., 2016) introduced a method to overcome this limitation and drive an arc in charging and electric mode in such a way, that the change between the two modes is done exactly at the point so that the battery charging is equal at the beginning and the end of this arc. Unfortunately, this has an impact on the required driving time for the arc, which might result in an infeasible solution in the case of the HEV-TSPTW due to a time window constraint violation. Therefore, we need to validate the feasibility of the Intra Arc Mode Change (IAMC) after each move and only

35

apply it, when the resulting solution is still valid. In this test setup, we start with a pure combustion solution as in the previous test, and apply the IAMC to all arcs as a post-processing step afterwards. The results are given in Table 8. Overall runtime for this test was similar to the previous one, which means that the additional effort for the IAMC is negligibly small. In accordance with the findings from the pure combustion test, the instances where the delivery area is far away from the depot are able to generate the highest savings, as the long arcs can now be driven in mixed mode. For several of the instances with the delivery area close to the depot, the regular P-VNS is superior to the IAMC applied to the pure combustion solution. This means, that the P-VNS is capable of finding mode combinations that generate overall higher savings. P-VNS plus Intra Arc Mode Change Results As a final test with our heuristic, we combined the P-VNS with the runtime optimized parameters and the IAMC to check if the initial solutions still have potential for additional improvement and to evaluate if the additional ’effort’ in terms of calculation time for the P-VNS is superior to the combustion plus IAMC. The setup with five runs is the same as for the previous tests and results are given in Table 9. The test took slightly more than 12.6 hours, so the time is relatively close to the time for the pure runtime optimized P-VNS setup; in fact for some of the instances the runtime was even slightly faster (factor below 1.0). Although, the combustion plus IAMC was able to provide significant savings compared to the P-VNS, the combination of P-VNS plus IAMC post processing outperforms any other test we conducted and should be the preferred setup for real world problems. 5. Conclusion and Outlook In this paper we present a new optimization problem considering delivery tours with time windows to serve customers with hybrid electric vehicles. We introduce a large set of 216 benchmark instances, representing real world delivery tours from a depot location. Small instances with 8 and 10 customers 36

Instance Instance Instance Instance Gap Factor Gap Factor Gap Factor Gap Factor Dist. Cust. TWs. Var. Dist. Cust. TWs. Var. Dist. Cust. TWs. Var. Dist. Cust. TWs. Var. 0 8 0 1 -1.11% 0.2 0 10 0 1 1.02% 0.1 0 20 0 1 1.93% 0.1 0 50 0 1 -4.42% 0.0 0 8 2 1 -2.14% 0.2 0 10 3 1 1.02% 0.2 0 20 8 1 1.93% 0.2 0 50 23 1 7.02% 0.2 0 8 6 1 -2.14% 0.2 0 10 8 1 0.88% 0.2 0 20 18 1 1.93% 0.2 0 50 48 1 6.90% 0.2 0 8 2+ 1 -2.14% 0.2 0 10 2+ 1 1.02% 0.2 0 20 2+ 1 1.93% 0.2 0 50 2+ 1 2.53% 0.1 0 8 4+ 1 -2.14% 0.4 0 10 5+ 1 1.02% 0.2 0 20 10+ 1 1.93% 0.3 0 50 25+ 1 7.02% 0.1 0 8 8+ 1 -2.14% 0.3 0 10 10+ 1 0.88% 0.3 0 20 20+ 1 1.93% 0.2 0 50 50+ 1 6.90% 0.2 0 8 0 2 -2.68% 0.1 0 10 0 2 -1.97% 0.1 0 20 0 2 2.16% 0.1 0 50 0 2 -2.53% 0.0 0 8 2 2 -2.68% 0.2 0 10 3 2 -2.16% 0.1 0 20 8 2 2.19% 0.2 0 50 23 2 6.13% 0.1 0 8 6 2 -2.68% 0.2 0 10 8 2 -2.16% 0.2 0 20 18 2 2.19% 0.2 0 50 48 2 6.23% 0.1 0 8 2+ 2 -2.68% 0.2 0 10 2+ 2 -2.16% 0.2 0 20 2+ 2 2.19% 0.2 0 50 2+ 2 -1.55% 0.1 0 8 4+ 2 -2.68% 0.2 0 10 5+ 2 -2.16% 0.2 0 20 10+ 2 2.19% 0.2 0 50 25+ 2 6.15% 0.1 0 8 8+ 2 -2.68% 0.2 0 10 10+ 2 -2.16% 0.2 0 20 20+ 2 2.19% 0.2 0 50 50+ 2 6.23% 0.2 0 8 0 3 -1.64% 0.1 0 10 0 3 1.38% 0.1 0 20 0 3 3.46% 0.1 0 50 0 3 2.61% 0.0 0 8 2 3 -2.36% 0.2 0 10 3 3 -2.12% 0.1 0 20 8 3 3.48% 0.1 0 50 23 3 5.81% 0.1 0 8 6 3 -3.38% 0.2 0 10 8 3 -2.12% 0.2 0 20 18 3 3.36% 0.2 0 50 48 3 5.19% 0.1 0 8 2+ 3 -1.67% 0.2 0 10 2+ 3 -2.12% 0.2 0 20 2+ 3 3.48% 0.1 0 50 2+ 3 0.61% 0.1 0 8 4+ 3 -1.67% 0.2 0 10 5+ 3 -2.12% 0.2 0 20 10+ 3 3.48% 0.2 0 50 25+ 3 4.89% 0.1 0 8 8+ 3 -3.38% 0.3 0 10 10+ 3 -2.12% 0.2 0 20 20+ 3 3.36% 0.2 0 50 50+ 3 5.19% 0.2 28 8 0 1 -6.45% 0.2 28 10 0 1 -6.18% 0.2 28 20 0 1 -5.32% 0.1 28 50 0 1 -4.40% 0.1 28 8 2 1 -6.45% 0.3 28 10 3 1 -6.27% 0.2 28 20 8 1 -5.32% 0.2 28 50 23 1 -3.12% 0.3 28 8 6 1 -6.45% 0.2 28 10 8 1 -6.48% 0.2 28 20 18 1 -5.32% 0.2 28 50 48 1 -2.78% 0.2 28 8 2+ 1 -6.64% 0.2 28 10 2+ 1 -6.27% 0.3 28 20 2+ 1 -5.32% 0.2 28 50 2+ 1 -4.98% 0.1 28 8 4+ 1 -6.64% 0.2 28 10 5+ 1 -6.27% 0.3 28 20 10+ 1 -5.33% 0.2 28 50 25+ 1 -2.75% 0.5 28 8 8+ 1 -6.64% 0.3 28 10 10+ 1 -6.48% 0.3 28 20 20+ 1 -5.32% 0.3 28 50 50+ 1 -1.05% 0.2 28 8 0 2 -6.85% 0.2 28 10 0 2 -6.52% 0.2 28 20 0 2 -5.54% 0.1 28 50 0 2 -6.83% 0.0 28 8 2 2 -7.36% 0.2 28 10 3 2 -6.52% 0.3 28 20 8 2 -5.52% 0.2 28 50 23 2 -3.58% 0.2 28 8 6 2 -7.36% 0.2 28 10 8 2 -6.52% 0.3 28 20 18 2 -5.52% 0.2 28 50 48 2 -3.58% 0.2 28 8 2+ 2 -7.36% 0.3 28 10 2+ 2 -6.52% 0.2 28 20 2+ 2 -5.52% 0.2 28 50 2+ 2 -7.70% 0.1 28 8 4+ 2 -7.36% 0.3 28 10 5+ 2 -6.52% 0.3 28 20 10+ 2 -5.52% 0.2 28 50 25+ 2 -3.58% 0.2 28 8 8+ 2 -7.36% 0.3 28 10 10+ 2 -6.52% 0.3 28 20 20+ 2 -5.52% 0.3 28 50 50+ 2 -3.56% 0.2 28 8 0 3 -6.98% 0.2 28 10 0 3 -6.05% 0.1 28 20 0 3 -5.42% 0.1 28 50 0 3 -7.29% 0.1 28 8 2 3 -6.98% 0.3 28 10 3 3 -6.05% 0.2 28 20 8 3 -5.39% 0.2 28 50 23 3 -3.60% 0.3 28 8 6 3 -6.98% 0.3 28 10 8 3 -6.05% 0.2 28 20 18 3 -5.31% 0.2 28 50 48 3 -3.37% 0.2 28 8 2+ 3 -6.98% 0.3 28 10 2+ 3 -6.05% 0.3 28 20 2+ 3 -5.31% 0.2 28 50 2+ 3 -5.37% 0.1 28 8 4+ 3 -6.98% 0.2 28 10 5+ 3 -6.05% 0.3 28 20 10+ 3 -5.31% 0.2 28 50 25+ 3 -3.49% 0.4 28 8 8+ 3 -6.98% 0.3 28 10 10+ 3 -6.05% 0.3 28 20 20+ 3 -5.31% 0.2 28 50 50+ 3 -3.37% 0.2 57 8 0 1 -6.78% 0.2 57 10 0 1 -7.09% 0.2 57 20 0 1 -6.38% 0.1 57 50 0 1 -6.81% 0.0 57 8 2 1 -6.78% 0.2 57 10 3 1 -7.19% 0.2 57 20 8 1 -6.38% 0.2 57 50 23 1 -5.50% 0.1 57 8 6 1 -6.93% 0.3 57 10 8 1 -7.19% 0.2 57 20 18 1 -6.38% 0.2 57 50 48 1 -4.83% 0.2 57 8 2+ 1 -6.78% 0.2 57 10 2+ 1 -7.19% 0.2 57 20 2+ 1 -6.27% 0.2 57 50 2+ 1 -6.93% 0.1 57 8 4+ 1 -6.78% 0.2 57 10 5+ 1 -7.19% 0.3 57 20 10+ 1 -6.38% 0.2 57 50 25+ 1 -4.83% 0.2 57 8 8+ 1 -6.93% 0.2 57 10 10+ 1 -7.19% 0.3 57 20 20+ 1 -6.38% 0.2 57 50 50+ 1 -4.83% 0.2 57 8 0 2 -7.08% 0.2 57 10 0 2 -6.91% 0.2 57 20 0 2 -6.39% 0.1 57 50 0 2 -6.51% 0.0 57 8 2 2 -7.29% 0.2 57 10 3 2 -7.02% 0.2 57 20 8 2 -6.37% 0.1 57 50 23 2 -5.34% 0.1 57 8 6 2 -7.29% 0.3 57 10 8 2 -7.02% 0.3 57 20 18 2 -6.37% 0.2 57 50 48 2 -5.38% 0.2 57 8 2+ 2 -7.29% 0.3 57 10 2+ 2 -7.02% 0.2 57 20 2+ 2 -6.37% 0.2 57 50 2+ 2 -7.05% 0.2 57 8 4+ 2 -7.29% 0.3 57 10 5+ 2 -7.02% 0.3 57 20 10+ 2 -6.37% 0.2 57 50 25+ 2 -5.38% 0.2 57 8 8+ 2 -7.29% 0.3 57 10 10+ 2 -7.02% 0.3 57 20 20+ 2 -6.37% 0.3 57 50 50+ 2 -5.38% 0.2 57 8 0 3 -7.25% 0.2 57 10 0 3 -7.06% 0.1 57 20 0 3 -6.31% 0.1 57 50 0 3 -6.10% 0.0 57 8 2 3 -7.25% 0.2 57 10 3 3 -7.06% 0.2 57 20 8 3 -6.31% 0.2 57 50 23 3 -5.68% 0.1 57 8 6 3 -7.50% 0.3 57 10 8 3 -7.06% 0.2 57 20 18 3 -6.31% 0.2 57 50 48 3 -5.19% 0.2 57 8 2+ 3 -7.69% 0.3 57 10 2+ 3 -7.06% 0.2 57 20 2+ 3 -6.31% 0.2 57 50 2+ 3 -6.10% 0.1 57 8 4+ 3 -7.69% 0.4 57 10 5+ 3 -7.06% 0.3 57 20 10+ 3 -6.31% 0.2 57 50 25+ 3 -5.16% 0.2 57 8 8+ 3 -7.50% 0.4 57 10 10+ 3 -7.06% 0.3 57 20 20+ 3 -6.31% 0.3 57 50 50+ 3 -5.20% 0.2

Table 8: Combustion plus Intra Arc Mode Change Results

37

Instance Instance Instance Instance Gap Factor Gap Factor Gap Factor Gap Factor Dist. Cust. TWs. Var. Dist. Cust. TWs. Var. Dist. Cust. TWs. Var. Dist. Cust. TWs. Var. 0 8 0 1 -1.60% 1.0 0 10 0 1 0.00% 1.0 0 20 0 1 -0.23% 1.0 0 50 0 1 -7.94% 1.0 0 8 2 1 -2.28% 1.1 0 10 3 1 0.00% 0.9 0 20 8 1 -0.23% 1.0 0 50 23 1 -0.26% 1.0 0 8 6 1 -2.31% 1.1 0 10 8 1 -0.11% 1.0 0 20 18 1 -0.23% 1.0 0 50 48 1 -0.33% 1.0 0 8 2+ 1 -2.28% 1.0 0 10 2+ 1 0.00% 1.0 0 20 2+ 1 -0.19% 1.0 0 50 2+ 1 3.13% 0.9 0 8 4+ 1 -2.38% 1.2 0 10 5+ 1 0.00% 1.1 0 20 10+ 1 -0.23% 1.0 0 50 25+ 1 -0.26% 1.0 0 8 8+ 1 -2.31% 1.2 0 10 10+ 1 -0.11% 1.1 0 20 20+ 1 -0.23% 0.9 0 50 50+ 1 -0.33% 1.0 0 8 0 2 -3.09% 1.0 0 10 0 2 -2.62% 0.9 0 20 0 2 -0.25% 0.9 0 50 0 2 -0.87% 0.9 0 8 2 2 -3.25% 1.0 0 10 3 2 -2.27% 1.0 0 20 8 2 -0.26% 1.0 0 50 23 2 -0.18% 1.0 0 8 6 2 -3.25% 1.0 0 10 8 2 -2.27% 1.0 0 20 18 2 -0.26% 0.9 0 50 48 2 -0.40% 1.0 0 8 2+ 2 -3.25% 1.0 0 10 2+ 2 -2.27% 1.1 0 20 2+ 2 -0.11% 1.0 0 50 2+ 2 -4.85% 1.0 0 8 4+ 2 -3.25% 1.0 0 10 5+ 2 -2.27% 1.1 0 20 10+ 2 -0.25% 1.0 0 50 25+ 2 -1.01% 1.1 0 8 8+ 2 -3.25% 1.2 0 10 10+ 2 -2.27% 0.9 0 20 20+ 2 -0.26% 1.1 0 50 50+ 2 -0.40% 1.1 0 8 0 3 -2.38% 1.0 0 10 0 3 0.00% 1.0 0 20 0 3 -0.59% 0.9 0 50 0 3 -0.22% 0.9 0 8 2 3 -2.58% 1.0 0 10 3 3 -2.89% 1.0 0 20 8 3 -0.59% 0.9 0 50 23 3 0.10% 1.1 0 8 6 3 -3.38% 1.0 0 10 8 3 -2.89% 1.0 0 20 18 3 -0.42% 1.0 0 50 48 3 0.00% 1.1 0 8 2+ 3 -1.67% 1.0 0 10 2+ 3 -2.89% 1.0 0 20 2+ 3 -0.29% 1.0 0 50 2+ 3 -2.70% 1.1 0 8 4+ 3 -1.67% 1.1 0 10 5+ 3 -2.89% 1.0 0 20 10+ 3 -0.59% 1.0 0 50 25+ 3 0.39% 1.0 0 8 8+ 3 -3.38% 1.1 0 10 10+ 3 -2.89% 1.0 0 20 20+ 3 -0.42% 1.0 0 50 50+ 3 0.05% 1.1 28 8 0 1 -6.44% 1.0 28 10 0 1 -6.28% 1.1 28 20 0 1 -5.82% 1.0 28 50 0 1 -3.36% 1.0 28 8 2 1 -6.42% 1.0 28 10 3 1 -6.37% 1.1 28 20 8 1 -5.78% 0.9 28 50 23 1 -5.04% 0.9 28 8 6 1 -6.42% 0.8 28 10 8 1 -6.59% 1.0 28 20 18 1 -5.84% 1.0 28 50 48 1 -5.12% 1.1 28 8 2+ 1 -6.64% 0.8 28 10 2+ 1 -6.37% 1.1 28 20 2+ 1 -5.82% 0.9 28 50 2+ 1 -6.45% 1.0 28 8 4+ 1 -6.60% 1.0 28 10 5+ 1 -6.37% 1.1 28 20 10+ 1 -5.82% 1.0 28 50 25+ 1 -4.12% 1.3 28 8 8+ 1 -6.60% 1.0 28 10 10+ 1 -6.59% 1.1 28 20 20+ 1 -5.84% 1.0 28 50 50+ 1 -2.64% 0.9 28 8 0 2 -6.95% 1.0 28 10 0 2 -6.70% 1.1 28 20 0 2 -6.07% 0.9 28 50 0 2 -3.96% 0.8 28 8 2 2 -7.36% 1.2 28 10 3 2 -6.58% 1.1 28 20 8 2 -6.04% 1.0 28 50 23 2 -5.33% 0.8 28 8 6 2 -7.36% 1.1 28 10 8 2 -6.58% 1.0 28 20 18 2 -6.06% 1.1 28 50 48 2 -5.39% 1.1 28 8 2+ 2 -7.36% 1.3 28 10 2+ 2 -6.58% 1.0 28 20 2+ 2 -5.94% 0.9 28 50 2+ 2 -4.70% 1.1 28 8 4+ 2 -7.36% 1.1 28 10 5+ 2 -6.58% 1.0 28 20 10+ 2 -6.05% 0.9 28 50 25+ 2 -5.38% 0.9 28 8 8+ 2 -7.36% 1.1 28 10 10+ 2 -6.58% 0.9 28 20 20+ 2 -6.05% 1.1 28 50 50+ 2 -5.36% 1.0 28 8 0 3 -6.98% 1.1 28 10 0 3 -6.27% 1.1 28 20 0 3 -5.64% 0.9 28 50 0 3 -5.36% 1.0 28 8 2 3 -6.98% 1.0 28 10 3 3 -6.27% 1.0 28 20 8 3 -5.86% 1.0 28 50 23 3 -5.34% 0.9 28 8 6 3 -6.98% 1.1 28 10 8 3 -6.27% 1.1 28 20 18 3 -5.78% 0.9 28 50 48 3 -5.19% 1.1 28 8 2+ 3 -6.97% 1.0 28 10 2+ 3 -6.27% 1.0 28 20 2+ 3 -5.72% 0.9 28 50 2+ 3 -4.69% 0.9 28 8 4+ 3 -6.97% 0.8 28 10 5+ 3 -6.27% 1.1 28 20 10+ 3 -5.78% 0.9 28 50 25+ 3 -5.23% 1.1 28 8 8+ 3 -6.97% 0.9 28 10 10+ 3 -6.27% 1.1 28 20 20+ 3 -5.79% 0.9 28 50 50+ 3 -5.14% 1.1 57 8 0 1 -6.78% 1.0 57 10 0 1 -7.15% 0.9 57 20 0 1 -6.72% 1.0 57 50 0 1 -4.94% 0.9 57 8 2 1 -6.78% 1.0 57 10 3 1 -7.20% 0.9 57 20 8 1 -6.75% 1.0 57 50 23 1 -6.42% 1.1 57 8 6 1 -6.93% 1.0 57 10 8 1 -7.20% 1.0 57 20 18 1 -6.79% 1.0 57 50 48 1 -6.28% 1.0 57 8 2+ 1 -6.78% 0.9 57 10 2+ 1 -7.20% 0.9 57 20 2+ 1 -6.65% 0.9 57 50 2+ 1 -6.35% 0.9 57 8 4+ 1 -6.78% 1.0 57 10 5+ 1 -7.20% 1.0 57 20 10+ 1 -6.79% 0.9 57 50 25+ 1 -6.25% 0.9 57 8 8+ 1 -6.93% 0.9 57 10 10+ 1 -7.20% 1.0 57 20 20+ 1 -6.79% 1.0 57 50 50+ 1 -6.28% 1.1 57 8 0 2 -7.12% 1.0 57 10 0 2 -6.99% 1.1 57 20 0 2 -6.89% 0.9 57 50 0 2 -5.47% 1.0 57 8 2 2 -7.33% 1.1 57 10 3 2 -7.07% 1.0 57 20 8 2 -6.87% 0.9 57 50 23 2 -6.29% 1.0 57 8 6 2 -7.33% 1.0 57 10 8 2 -7.07% 1.0 57 20 18 2 -6.87% 1.0 57 50 48 2 -6.27% 1.0 57 8 2+ 2 -7.33% 1.2 57 10 2+ 2 -7.07% 1.1 57 20 2+ 2 -6.87% 0.9 57 50 2+ 2 -5.98% 1.1 57 8 4+ 2 -7.33% 1.1 57 10 5+ 2 -7.07% 1.0 57 20 10+ 2 -6.87% 0.9 57 50 25+ 2 -6.27% 1.0 57 8 8+ 2 -7.33% 1.0 57 10 10+ 2 -7.07% 1.0 57 20 20+ 2 -6.87% 1.0 57 50 50+ 2 -6.27% 1.0 57 8 0 3 -7.34% 1.0 57 10 0 3 -7.13% 1.0 57 20 0 3 -6.79% 0.9 57 50 0 3 -6.44% 0.8 57 8 2 3 -7.34% 1.0 57 10 3 3 -7.13% 1.0 57 20 8 3 -6.80% 1.0 57 50 23 3 -5.82% 0.9 57 8 6 3 -7.50% 1.1 57 10 8 3 -7.13% 1.0 57 20 18 3 -6.80% 1.0 57 50 48 3 -6.41% 1.1 57 8 2+ 3 -7.69% 1.3 57 10 2+ 3 -7.13% 1.1 57 20 2+ 3 -6.79% 1.0 57 50 2+ 3 -5.91% 0.9 57 8 4+ 3 -7.69% 1.0 57 10 5+ 3 -7.13% 1.1 57 20 10+ 3 -6.79% 0.9 57 50 25+ 3 -6.36% 1.0 57 8 8+ 3 -7.50% 1.2 57 10 10+ 3 -7.13% 1.0 57 20 20+ 3 -6.80% 1.1 57 50 50+ 3 -6.41% 0.9

Table 9: Runtime Optimized plus Intra Arc Mode Change Results

38

are solved exactly with CPLEX, based on a mixed-integer programming formulation. In addition to the small instances, we also solved 20 and 50 customer instances with a Variable Neighborhood Search heuristic. The results show that the heuristic is able to provide exact solutions for the small instances, in general, with a fraction of the required calculation time compared to CPLEX. The numerical studies indicate that also for the larger instances the solution quality is reasonably good, although, optimality cannot be proven. Additionally, our heuristic was able to find new best solutions for an existing problem, strengthening the impression of the sound solution quality. Solving the Traveling Salesman Problem with Time Windows with regular combustion vehicles and comparing this with our Hybrid Electric Vehicle results, we are able to estimate the saving such a Hybrid Vehicle is able to generate compared to a regular one. Applying an Intra Arc Mode Change on our Hybrid Vehicle solutions, it is possible to even increase these potential savings. Considering the insignificant calculation time for this additional step, this setup is the recommendation to solve the problem for a daily tour planning on an operational base. Although we see still potential to improve the solution quality of our heuristic, the consequential next step for the optimization of Hybrid Electric Vehicles would be the consideration of more than one vehicle, extending the Traveling Salesman Problem to a Vehicle Routing Problem. However, as the problem is already very difficult to solve for a single vehicle, a deeper analysis on the overall feasibility of such an approach is required.

References References ARADEX AG (2019). Retrofit electric drive kit for diesel delivery vehicles. Online. Ascheuer, N., Fischetti, M., & Gr¨ otschel, M. (2001). Solving the asymmetric

39

travelling salesman problem with time windows by branch-and-cut. Mathematical Programming, 90 , 475–506. Chan, C. C. (2002). The state of the art of electric and hybrid vehicles. Proceedings of the IEEE , 90 , 247–275. Chan, C. C. (2007). The state of the art of electric, hybrid, and fuel cell vehicles. Proceedings of the IEEE , 95 , 704–718. Davidovic, T., Crainic, T. G., & Guertin, F. (2013). Parallelization strategies for variable neighborhood search. In CIRRELT-2013-47 . Desaulniers, G., Errico, F., Irnich, S., & Schneider, M. (2016). Exact Algorithms for Electric Vehicle-Routing Problems with Time Windows. Operations Research, 64 , 1388–1405. Desaulniers, G., Madsen, O. B., & Ropke, S. (2014). Chapter 5: The vehicle routing problem with time windows. In Vehicle Routing chapter 5. (pp. 119– 159). Society for Industrial and Applied Mathematics. Doppstadt, C., Koberstein, A., & Vigo, D. (2016). The hybrid electric vehicle traveling salesman problem. European Journal of Operational Research, 253 , 825 – 842. Doppstadt, C., Koberstein, A., & Vigo, D. (2019). The hybrid electric vehicle - traveling salesman problem with time windows (HEV-TSPTW) instances. Mendeley Dataset. Gendreau, M., & Potvin, J.-Y. (2010). Handbook of Metaheuristics. (2nd ed.). Springer Publishing Company, Incorporated. Hansen, P., Mladenovi´c, N., & Moreno P´erez, J. A. (2010). Variable neighbourhood search: methods and applications. Annals of Operations Research, 175 , 367–407.

40

L´opez-Ib´ an ˜ez, M., & Blum, C. (2010). Beam-ACO for the travelling salesman problem with time windows. Computers & Operations Research, 37 , 1570 – 1583. Mancini, S. (2017). The hybrid vehicle routing problem. Transportation Research Part C: Emerging Technologies, 78 , 1–12. Mladenovi´c, N., & Hansen, P. (1997). Variable neighborhood search. Computers & Operations Research, 24 , 1097 – 1100. Ohlmann, J., & Thomas, B. (2007). A compressed annealing approach to the traveling salesman problem with time windows. Informs Journal on Computing - INFORMS , . Savelsbergh, M. W. P. (1985). Local search in routing problems with time windows. Annals of Operations Research, 4 , 285–305. Schneider, M., Stenger, A., & Goeke, D. (2014). The electric vehicle-routing problem with time windows and recharging stations. Transportation Science, 48 . da Silva, R. F., & Urrutia, S. (2010). A general vns heuristic for the traveling salesman problem with time windows. Discrete Optimization, 7 , 203 – 211.

41