A metaheuristic for the Two-Echelon Vehicle Routing Problem
with Covering Options
D.L.J.U. Enthoven
Master’s Thesis Econometrics, Operations Research and Actuarial Studies Supervisor: Prof. dr. K.J. Roodbergen
A metaheuristic for the Two-Echelon Vehicle Routing
Problem with Covering Options
D.L.J.U. Enthoven
15th July 2018
Abstract
1
Introduction
Over the next decades, new challenges will emerge in parcel delivery logistics, due to the increasing world population living in cities, the growing popularity of e-commerce and the desire for faster delivery (Savelsbergh and Van Woensel, 2016). As a consequence, more delivery vehicles are required to satisfy all demand. The increased freight activity, resulting in more traffic congestion and air pollution, will negatively effect the quality of life in cities (Demir et al., 2015). Therefore, the challenge for logistical service providers is to introduce new distribution systems that preserve the liveability of cities, while simultaneously ensuring fast deliveries.
In order to reduce these negative externalities, we examine a delivery structure where parcels are not directly delivered to the customers from the depot, which is typically located far away from the customers. Instead, trucks with large capacities will transport the parcels from the depot to satellite locations and covering locations, intermediate facilities positioned between the depot and the customers. Satellite locations are small depots, from where the customers can be serviced by city freighters, small vehicles that are more environmentally friendly than trucks and better equipped to perform deliveries in densely populated areas (Van Der Helm, 2015). Covering locations are facilities, for example parcel lockers, where the parcels of nearby customers can be stored. Thereafter, customers can collect the parcels themselves at a convenient moment. The benefit of these locations for distributors is that these customers do not have to be reached directly (Deutsch and Golany, 2018; Arnold et al., 2018). A representation of this framework is depicted in Fig. 1.
Since each parcel traverses through one of the two intermediate facilities before reaching the customers, two levels of routing are created: the routes from the depot to the intermediate facilities (first echelon) and the routes from the satellite locations to the customers (second echelon). In the classic two-echelon vehicle routing problem (2E-VRP1; Perboli et al., 2011) only satellite locations
were available. One of the advantages of this framework is that different vehicles, trucks and city freighters, with preferred characteristics to tackle the echelon specific difficulties can be employed. Our setting extends the 2E-VRP by introducing covering locations as an additional delivery option. Since customers can collect the parcels themselves at the covering locations, the probability of failed deliveries is reduced and because less parcels have to pass through the satellite locations, the amount of city freighters required is decreased.
Depot Covering locations Satellite locations Customers 1st echelon routes 2nd echelon routes Cover radius
Figure 1: Example of a solution for the 2E-VRP-CO. The first echelon routes (trucks) and second echelon
routes (city freighters) are represented by solid and dotted lines, respectively. The satellite locations and covering locations which are available but currently not used are depicted in grey.
extra consolidation step, which is undesirable in a market where high speed delivery is important. Secondly, city freighters have a limited cargo capacity. Therefore they are less suited to replenish covering locations which typically serve multiple customers.
The goal of the 2E-VRP-CO problem is to determine (1) how to service all the customers, (2) which facilities of all available locations to employ and (3) how to construct the vehicle routes on the first and second echelon. The objective of the problem is to make these decisions in such a way that the routing and connection costs are minimized.
We propose a mathematical formulation for the 2E-VRP-CO. This formulation is used to solve small instances of the problem to optimality and is tested against alternative formulations. We show that the problem is a generalization of the 2E-VRP and the simultaneous facility location and vehicle routing problem (SFL-VRP) as introduced by Veenstra et al. (2018) without duration constraints. We develop an adaptive large neighborhood search (ALNS) metaheuristic to efficiently tackle larger instances. To assess the quality of the heuristic, computational experiments are performed on the well known instances for the 2E-VRP and the SFL-VRP. Finally, we provide a detailed discussion on the effect of different problem specific parameters and other managerial insights.
To the best of our knowledge, this configuration of routing with covering options has not been not studied before. A closely related problem is examined by Zhou et al. (2018). In their paper parcels from multiple depots have to be delivered to the customers through either home delivery or pickup locations, within a two echelon framework. Although the pickup lockers show similarities to the covering locations in our problem there are two main differentiating features. First, the covering locations we examine can only be visited by trucks, for the reasons mentioned before, while Zhou et al. (2018) solely allow the city freighters to access their pickup locations. Second, we allow for split deliveries on the first echelon, where Zhou et al. (2018) do not. Consequently their problem is not a generalization of the 2E-VRP. Zhou et al. (2018) propose a hybrid multi-population genetic algorithm to construct the vehicle routes.
the location routing problem (LRP, Zhou et al., 2016; Stenger et al., 2012), which generalizes the VRP by providing the option to select the locations of the depots. The 2E-VRP-CO differs from these three papers in that connection costs are considered and no fixed opening costs for covering locations are incurred. Furthermore, vehicles have to be routed on two levels instead of only one. Moreover, Veenstra et al. (2018) force all customers within the covering radius of an opened covering locations to be served by that location, while we still allow city freighters to satisfy these customers.
As the 2E-VRP-CO generalizes the 2E-VRP, we provide an overview of the most relevant 2E-VRP papers (Appendix A). For the 2E-VRP exact solution methods are proposed by Baldacci et al. (2013), Jepsen et al. (2013) and Santos et al. (2014), while Hemmelmayr et al. (2012) and Breunig et al. (2016) respectively introduce an ALNS and LNS (large neighborhood search) heuristic. Also extensions on the 2E-VRP that include more practical considerations or desired properties are investigated. Grangier et al. (2016) and Anderluh et al. (2017) expand the problem by considering time synchronization constraints at the satellite locations. They propose an ALNS and a greedy randomized adaptive search procedure (GRASP) with path relinking heuristic respectively. Two-echelon vehicle routing with environmental considerations is examined by Wang et al. (2017). They proposed a variable neighborhood search (VNS) heuristic to solve the problem.
The remainder of this paper is organized as follows. In Section 2, we define and develop the 2E-VRP-CO formally and we show how it generalizes the 2E-VRP and SFL-VRP problems. In Section 3 we describe our proposed ALNS heuristic. In Section 4 the results of the numerical experiments and managerial insights are presented. Conclusions are stated in Section 5.
2
Problem Statement and mathematical formulation
The 2E-VRP-CO is defined on a directed graph G = (V, A), where V is the set of all vertices. Let vertex V0= {0}denote the depot, VLthe set of covering locations, VS the set of satellite locations
and VC the set containing all customers, such that V = V0∪ VL∪ VS∪ VC. The set consisting of
the arcs between all vertices is expressed by A. The length of each arc (i, j) ∈ A, denoted by dij,
is nonnegative and satisfies the triangle inequality. Each customer, i ∈ VC, has a positive demand of q
i. As discussed, two different types of vehicles
are available to satisfy all requirements. On the first echelon at most m1trucks can be employed
to deliver the parcels from the depot to the intermediate facilities. All trucks have a capacity limit of Q1. The number of city freighter available to deliver the parcels from the satellite locations
to the customers is equal to m2. For both trucks and city freighters, the execution of multiple
trips by the same vehicle is not allowed and not all vehicles have to be used. All trucks start and end their route at the depot and each city freighter starts and ends its route at the same satellite location. Each city freighter can be assigned to one of the available satellite locations, provided that at most ¯mk city freighters can depart from each satellite location k ∈ VS. Each customer can
only be visited by one city freighter or serviced via one covering location, hence split deliveries on the second echelon are not allowed. On the other hand both types of intermediate facilities can be serviced by multiple trucks. The driving cost associated with traversing over an arc (i, j) in A is denoted by cij. This cost is proportional to the length and the specific type of vehicle passing over
arc (i, j), for all (i, j) ∈ A. In addition, covering location j can service customers located within a radius of rjfrom its position. The cost of connecting a customer to a covering location is expressed
The decision variables used to optimize the 2E-VRP-CO can be separated in three groups. First,
xij and w1ijare integer variables describing respectively the amount of trucks and the amount of
parcels on the first echelon arcs (i, j). In accordance, for the second echelon arcs these values are governed by yijkand wijk2 , where k specifies the satellite location from which the city freighter
departs. Since customers can only be serviced by at most one city freighter, yijkis a binary variable.
Second, binary variable zs
kidenotes if customer i is serviced through satellite location k. On the
other hand zc
kiindicates if customer i is serviced by covering location k. Finally, three auxiliary
variables are introduced to mitigate the formulation. Binary variable vkdenotes if covering location
kis used, while Ds k and D
t
krepresent the amount of parcels passing through each satellite location
and covering location k. An overview of the notation used in this paper is provided in Table 1.
The notation and flow based mixed integer linear program which we present are based on the formulation by Perboli et al. (2011) for the 2E-VRP. The covering location specific constraints and their notation are build upon the SLF-VRP model presented by Veenstra et al. (2018). For the 2E-VRP-CO we propose the following formulation:
X i∈V0∪VL∪VS wij1 − X i∈V0∪VL∪VS wji1 = Dsj ∀j ∈ VS (15) X i∈V0∪VL∪VS wij1 − X i∈V0∪VL∪VS wji1 = Dcj ∀j ∈ VL (16) X i∈V0∪VL∪VS wi01 − X i∈V0∪VL∪VS w10i= X i∈VC −qi (17) w1ij ≤ Q1xij ∀i, j ∈ V0∪ VL∪ VS (18) X i∈VC∪k Q2ijk− X i∈VC∪k Q2jik = zkjs dj ∀k ∈ VS, ∀j ∈ VC, (19) X i∈VC Q2ikk− X i∈VC Q2kik= −Djs ∀k ∈ VS (20)
w2ijk≤ Q2yijk ∀k ∈ VS, ∀i, j ∈ VS∪ VC (21)
ykij∈ {0, 1} ∀k ∈ VS, ∀i, j ∈ VS∪ VC (22) xij ∈ N0 ∀i, j ∈ V0∪ VL∪ VS (23) zkis ∈ {0, 1} ∀k ∈ VS, ∀i ∈ VC (24) zkic ∈ {0, 1} ∀k ∈ VL, ∀i ∈ VC (25) vk∈ {0, 1} ∀k ∈ VL (26) w1ij ≥ 0 ∀i, j ∈ V0∪ VL∪ VS (27) w2ijk≥ 0 ∀k ∈ VS, ∀i, j ∈ VS∪ VC (28)
The objective function (1) minimizes the sum of the driving cost of the first and second echelon routes, and the connection costs between customers and covering locations. Constraints (2) - (4) guarantee that the amount of trucks and city freighters used do not exceed the predefined limits. Constraints (5) - (8) are flow conservation constraints ensuring that the same amount of vehicles arrive and depart at each location. Constraints (9) and (10) govern that each covering location and satellite location can only service customers if it is visited by at least one truck. Constraint (11) assigns each customer to a satellite location or a covering location. Constraint (12) manages that only customers within the covering range can be serviced by a covering location. Constraints (10) - (12) govern the covering locations and are adaptations of constraints (1) - (6) by Veenstra et al. (2018)
Furthermore, constraints (13) - (21) supervise the amount of parcels traversing over all the arcs. Additionally, these constraints prevent the existence of subtours. Within constraints (13) and (14) the auxiliary variables denoting the requested demand of parcels at the satellite locations and covering locations are defined. Constraints (15) - (17) control the parcel flow on the first echelon arcs, similarly constraints (19) and (20) control the parcel flow on the second echelon. Inequalities (18) and (21) ensure that for both types of vehicles the capacity limits are not exceeded. Finally, constraints (22) - (28) define the nature and the range of the decision variables.
Table 1: Table of notation
Sets V = set of all nodes
V0 = the depot
VL = set of covering locations
VS = set of satellite locations
VC = set of customer locations
Parameters
m1 = The amount of trucks available
m2 = The amount of city freighters available
¯
mk = The maximum number of city freighters allowed to depart from
satellite location k
Q1 = The capacity of a truck
Q2 = The capacity of a city freighter
cij = The transportation cost from node i to node j
`ji = The connection cost between covering location j and customer i
α = the distance dependent connection cost factor
β = the fixed connection cost factor
qi = The demand of customer i
rj = The radius around covering location j in which customers can be
served.
dij = The distance between locations i and j.
Decision Variables
w1
ij = Flow passing through the first echelon arc (i, j)
wijk2 = Flow passing through the second echelon arc (i, j) and originating from satellite location k
xij = Number of trucks using the first echelon arc (i, j)
yijk =
1 if the second echelon arc (i, j) is used by a city freighter starting from satellite location k;
0 otherwise
zkis =
(
1 if customer i is served by satellite location k; 0 otherwise
zkic =
(
1 if customer i is served by covering location k; 0 otherwise
vk =
(
1 if covering location k is used; 0 otherwise
Dsj = The amount of parcels passing through each satellite location j
Dc
2.1
Valid Inequalities
To strengthen the mathematical formulation, we impose additional valid inequalities. We test cuts derived from the VRP and 2E-VRP literature and accept the cuts that provide improvements for the 2E-VRP-CO. In particular we apply two sets of cuts. Namely, flow-based cuts, which provide tighter bounds to existing constraints, and edge cuts, derived from subtour elimination constraints. We introduce flow based cuts, additional valid inequalities that provide tighter bounds for existing constraints. For example, constraint (21) can be strengthened by considering that each customer reduces the flow by an amount equal to its demand. Moreover, we state that a vehicle returning to its starting position should have no excess demand. Therefore, the following inequalities are valid: w2ijk≤ (K2− q i)yijk ∀k ∈ VS, ∀i, j ∈ VC (29) X i∈VL∪VS Q1i0= 0 (30) X i∈VC Q2ikk= 0 ∀j ∈ VS (31) w1ii= 0, xii= 0 ∀i ∈ V0∪ VL∪ VS (32) w2iik= 0, yii= 0 ∀i ∈ VS ∪ VC, ∀k ∈ VS (33)
Moreover, all flow decision variables going from node i to the same node i are fixed to zero by valid inequalities (32) and (33). Since this second set of valid inequalities provide tighter bounds without introducing additional constraints, we implement them without providing additional tests.
Furthermore, edge cuts introduce well-known subtour elimination constraints which are commonly used in the travelling salesman problem and the vehicle routing problem. They can be expressed as follows:
X
i,j∈VC 0
yijk ≤ |VC 0| − 1 ∀k ∈ VS, ∀VC 0⊂ VC, 2 ≤ |VC 0| ≤ |VC| − 2, (34)
where VC 0is a subset of the customers. The intuition behind these constraints is that subtours are
prevented if the amount of routes within a subset of customers is lower than the cardinality of that subset, that is the number of customers within the subset. The subsets with cardinality of one are excluded since constraints (32) and (33) already prevent that arcs from one location to itself are used.
Based on Perboli and Tadei (2010), we strengthen these inequalities by taking into account by which satellite location the customers are being serviced. The cardinality of the subset of customers VC 0can be expressed by a sum of variables of subsets from VC 0, such that the sum of variables z
kj
is equal to the amount of customers in VC 0minus one. Therefore, constraint (34) can be written
where VC 0is a subset of customers. The amount of total valid inequalities (34) and (35) increases
exponentially for higher levels of cardinality. Therefore we test for which level the best improve-ments can be found. Results from these tests can be found in Section 4.4.1. Similar inequalities could also be considered for the first echelon variables xij. However, these additional inequalities
are not included because they only provide marginal improvements, due to the limited number of satellite locations and covering locations relative to the amount of customers.
2.2
Generalizations of the 2E-VRP-CO
Proposition 1. The 2E-VRP-CO with VL= ∅is equivalent to the 2E-VRP.
Proof: Let an instance be given with VL = ∅. As no covering locations are available in this
instance all customers have to be reached by city freighters. Therefore, decision variables zc ki, vk
and Dc
jas well as constraints (11, 12, 14 and 16) become redundant. The objective of the problem
becomes to solely minimize the routing costs of the trucks and city freighters. The remaining model is equivalent to the 2E-VRP as presented by Perboli et al. (2011).
Proposition 2. The 2E-VRP-CO with VS = {0}, `
ji = 0 ∀j ∈ VL, ∀i ∈ VC and additional
con-straints xij ≤ 1 ∀j ∈ VL, ∀i ∈ V0∪ VS∪ VL (36) rjvj≤ dij+ X k∈VL zkic M ∀j ∈ VL, ∀i ∈ VC (37)
is equivalent to the SFL-VRP as introduced by Veenstra et al. (2018) without duration constraints. Proof: Let an instance of the specific 2E-VRP-CO be given as described above. In this instance
there is only one satellite location and it is located at the same position as the depot. The single satellite location is automatically covered by the depot, since the distance and therefore driving cost are reduced to zero. Therefore, the first echelon trucks now only have to service the covering locations. As mentioned in the introduction the covering locations can be regarded as lockers and are therefore similar to the medication lockers as applied in Veenstra et al. (2018). Constraint (36) replaces constraint (23), allowing only one truck to traverse on a first echelon arc. For this specific problem routing constraints (1) - (9) are identical to constraints (7) - (10) and (16) - (19) and cargo constraints (13) - (21) are equivalent to constraints (11) - (14) and (20) - (23) by Veenstra et al. (2018). Therefore the first and second echelon routes are identical to locker and patient routes, as presented in Veenstra et al. (2018).
Remaining is the difference between the lockers in Veenstra et al. (2018) and the covering locations as presented in this paper. One distinction is that Veenstra et al. (2018) force all customers within reach of an opened covering location to be serviced by that location. In the 2E-VRP-CO this restriction is not made. To account for this difference additional constraint (37) can be introduced, which is a combination of constraints (2) and (6) in Veenstra et al. (2018). Therefore, inequalities (37) and (10) - (12) are identical to constraints (1) - (6) by Veenstra et al. (2018). Another distinction is that Veenstra et al. (2018) consider fixed opening costs for the covering locations instead of connection costs. To account for this difference all connection costs could be set to zero,
`ji= 0, ∀j ∈ VL, ∀i ∈ VC. To take the fixed opening cost into account the cost of traversing on
an arc to an opened covering location j could be increased, such that cij ← Fj+ cij, where Fj
costs cannot be incurred multiple times, because each covering location can only be reached by a single truck due to constraint (36). Therefore the SFL-VRP, as presented by Veenstra et al. (2018), without duration constraints is as a special case of the 2E-VRP-CO.
3
Adaptive large neighborhood search heuristic
In this section we present our adaptive large neighborhood search (ALNS) heuristic to solve the 2E-VRP-CO. ALNS was first proposed by Ropke and Pisinger (2006) as an extension on the large neighborhood search (LNS) heuristic introduced by Shaw (1997) and the closely related ruin-and-recreate principle by Schrimpf et al. (2000). The aim of both ALNS and LNS is to improve solutions by iteratively applying destroy and repair operators such that a large range of potential solutions is examined. ALNS differs from the LNS heuristic in that the probability of selecting one of the operators is determined based on its past successes. Thereby, instance specific characteristics are indirectly taken into account. ALNS and LNS have gained a lot of popularity within routing problems as a result of its ability to explore different large neighborhoods around a current solution. Among these instances are many two echelon vehicle routing problems: the standard 2E-VRP (Hemmelmayr et al., 2012; Breunig et al., 2016), 2E-VRP with satellite synchronization (Grangier et al., 2016) and 2E-VRP with charging stations (Breunig et al., 2018). The structure of our algorithm, based on the framework by Hemmelmayr et al. (2012), is presented in Algorithm 1.
The procedure starts with creating an initial solution (line 1). Furthermore, the procedure initializes the operator probabilities π, specifying the probability that a specific operator is utilized, and the counters denoting the amount of iterations since the last local and global improvements were found,
ilocal, iglobalare initialized (line 2 - 3). Subsequently, we enter an iterative phase where destroy
and repair operators are used to change the solution (lines 4 - 27). First the current solution is destroyed (lines 6 - 11). A distinction is made between two types of destroy operators; DLare
the large impact operators, which change the configuration of satellite locations and covering locations that can be used, and DS are the operators that only effect a limited part of the solution.
A large destroy operator is only used after ωgraceiterations without improvements. The iterations
between the large destroy operators are referred to as the grace period. The grace period ensures that a new configuration of opened intermediate facilities can be examined thoroughly. A new solution is accepted based on a simulated annealing criterion or if a large destruction operator has just been employed (line 13 - 16). Each accepted solution is passed through a local search step (line 14). We always accept a new solution after large impact destroy operators for two reasons. First, since these operators change the solutions dramatically the repair operators are typically not capable to provide a good solution directly. Therefore, a local search is necessary to further improve the solution. Second, an exploration of each new configuration of opened satellite locations and covering locations is guaranteed. In lines (17 - 20) we check whether the new found solution is better than the best solution found during the current grace period. In that case the variable ilocalis reset to 0, thereby extending the current grace period. In lines (21 - 23) we test
whether the newly constructed solution is better than the best solution found within the entire process. If so, we store the solution and reset the variable iglobal. If there are ωglobal iterations
without finding a new global improvement we restart the iterative procedure from the current best solution found (lines 24 - 26). At the end of each iteration the operator probabilities π and the counter variables ilocaland iglobalare updated. The iterations continue until the stopping criterion
Algorithm 1 ALNS Heuristic 1: sbest, s∗, s ←InitialSolution 2: π ←InitializeOperatorProbabilities 3: ilocal, iglobal← 0 4: while StoppingCriterionNotReached do 5: s0← s
6: if ilocal< ωgracethen 7: s0 ← Destroy(s0, D S, π) 8: else 9: s0 ← Destroy(s0, D L, π) 10: ilocal← 0 11: end if 12: s0← Repair(s0, π) 13: if AcceptanceCriterion(s0, s) or i local = 0then 14: s0 ← LocalSearch(s0) 15: s ← s0 16: end if 17: if f (s0) < f (s∗)or i local= 0then 18: s∗← s0 19: ilocal← 0 20: end if 21: if f (s0) < f (sbest)then 22: sbest ← s0 23: end if
24: if iglobal> ωglobalthen
25: s ← sbest
26: end if
27: π ←UpdateOperatorProbabilities(π)
28: ilocal← ilocal+ 1, iglobal← iglobal+ 1 29: end while
30: Result sbest
The remainder of this section is structured as follows. Section 3.1 provides details on how the initial solution is created. In Sections 3.2 and 3.3 an overview of the destroy and repair operators is presented. Section 3.4 describes the acceptance mechanism. Section 3.5 specifies the local search procedures employed on the newly accepted solution. Finally, we explain in Section 3.6 the adaptive updating mechanism for the operator selection probabilities.
3.1
Initial solution
insertion position is available in any of the routes. In such a situation a new initial solution is created from scratch.
3.2
Destroy operators
Our proposed ALNS heuristic uses 13 different destroy operators. The first six operators are large impact operators. They change the set of satellite locations and covering locations available. Of these methods, the Satellite removal/Satellite opening/Satellite Swap operators are adaptations of the similar named operators by Hemmelmayr et al. (2012). Based on these operators, we introduce three new operators for the covering locations: Covering removal/Covering opening/Covering Swap. The remaining operators are small impact operators. Operators 7 - 11 are common operators within ALNS heuristics for vehicle routing problems (Ropke and Pisinger, 2006; Hemmelmayr et al., 2012). The Least Used Vehicle Removal operator is adapted from Grangier et al. (2016).
The small impact operators select q customers to be removed from the current solution, where for each iteration q is a random integer within range [q, ¯q]. Most of the removal operators sort the customers according to a specific criterion, such that the customers with low scores have a higher chance to be selected for removal. To prevent that an operator continuously selects the same customers for removal a randomness process is introduced, based on Franceschetti et al. (2017). Consider a ranking of |VC| customers according to some metric. The node selected to be removed is
the y-th customers, with y = dU [0, 1]ρ·|VC|e, where ρ is a measure of randomness. This randomness
procedure can be repeated until q customers are selected. High values of ρ lead to a higher probability of selecting nodes with smaller indices. Therefore, more randomness is introduced for low values of ρ. Note that for ρ = 1 the selection method is completely random.
We now provide a detailed description of the removal operators used in our ALNS heuristic: 1. Satellite Removal
The Satellite removal operator selects one satellite location at random and closes it. All customers serviced by this satellite location are removed from the current solution. Until a closed satellite location is opened at a later iteration, it cannot service any customers. This method allows to diversify our search over a framework consisting of only a subset of the satellite locations. To use only a fraction of all available satellite locations could result in better solutions if the decrease in first echelon routing costs is larger than the increase in second echelon routing costs. A safety mechanism is implemented such that always enough satellite locations are open to service all customers. This prevents situations where not enough satellite locations are available for a feasible solution. This operator deviates from the similar named operator by Hemmelmayr et al. (2012), since we do not open another satellite location immediately afterwards.
2. Satellite Opening
The Satellite opening operator selects a random satellite location from the set of closed satellite locations to open. To intensify the search around the newly opened location, the nearest q customers are removed from the current solution.
3. Satellite Swap
the q − qclosenearest customers to the newly opened location. We encourage swaps between
satellite locations which are positioned close together. Therefore, the probability to open a satellite location is inversely related to the distance with the newly closed satellite location. 4. Covering Removal
The Covering Removal operator is a new operator designed specifically for the 2E-VRP-CO. This operator changes the set of available covering locations and works similar to the Satellite Removal operator for the satellite locations. A covering location to be closed is selected at random from the set of opened covering locations. All customers currently serviced by the closed covering location are removed from the current solution.
5. Covering Opening
The Covering Opening operator is a new operator designed specifically for the 2E-VRP-CO. This operator increases the set of available covering locations, similar to the Satellite Opening operator for the satellite locations. From the set of closed covering locations one covering location is selected to be opened. To intensify the search around this location, all customers within the covering range of the opened covering location are removed from the solution. 6. Covering Swap
The Covering Swap operator is a new operator designed specifically for the 2E-VRP-CO. First an open covering location is closed and thereafter a different covering location, from the set of closed covering locations, is opened. All customers currently serviced by the closed covering location are removed from the solution. Additionally, customers in the covering range of the opened covering location are removed from the solution until q customers are removed in total. We encourage swaps between covering locations which are positioned close together. Therefore, the probability to open a covering location is inversely related to the distance with the newly closed covering location.
7. Random Removal
The Random Removal operator selects an amount of q customers at random and removes them from the solution.
8. Related Removal
The Related Removal operator removes customers who are positioned near each other. First a single customer is selected at random, the seed customer. Similar to the related removal operator proposed by Hemmelmayr et al. (2012) and Breunig et al. (2016) we remove the other q − 1 customers based on the distance to the seed customer. The nearest customers to the seed customer have the highest probability of being selected.
9. Worst Removal
The Worst Removal operator sorts the customers according to the removal gain. This is calculated as the difference in objective value between the current solution and the solution without the customer. Provided the ranking of customers based on their removal gain, q customers are selected by the randomization process described before.
10. Minimum Quantity Removal
The Minimum Quantity Removal operator removes customers with low demand. The ran-domization process is used to select the customers to remove. The idea behind this operator is that customers with low demand can be moved more easily to other routes in the solution, since they are less restricted by the capacity constraints.
11. Random Route Removal
city freighter from the current solution. If the selected route has more than q customers the route is only partially removed. If the selected route has less than q customers, additional routes are selected until q customers are removed.
12. Least Used Vehicle Removal
The Least Used Vehicle operator selects all customers from the route with the lowest capacity to be removed. As the Random Route Removal operator this method is reiterated until an amount of q customers are removed. The idea behind the operator is to minimize the amount of vehicles used. Since, routes which only service a few customers are on average less efficient. This method is based on the similar named method by Grangier et al. (2016) and the ‘Remove single node routes’ by Breunig et al. (2016).
13. Related Customer Removal with Delivery Method
The related Customer Removal with Delivery Method is introduced specifically for this problem. The operator is similar to the Related Customer Removal operator in that a seed customer is selected and other customers are chosen to be removed based on their distance to the seed customer. This operator differs from the Related Customer Removal operator in that it only selects nearby customers who are serviced by the same method as the seed customer. Therefore all customers removed are serviced either via covering locations or via satellite locations. The intuition behind this operator is that it intensifies the search regarding the current intermediate facilities used.
3.3
Repair operators
After the the destroy phase is executed, the removed customers are reinserted into the solution by repair operators. The repair operators only consider feasible insertion options where customers are serviced via one of the open intermediate facilities. The overall framework of the repair phase is similar to Breunig et al. (2016). First all customers are inserted into the solution. Thereafter, knowing the requested demand at the satellite location and covering locations, the first echelon routes are reconstructed. The ALNS heuristic we propose adopts three insertion operators. The
Greedy Insertion operator is a standard ALNS heuristic within the 2E-VRP literature. The Greedy Insertion Perturbed operator is adopted from Hemmelmayr et al. (2012). For the 2E-VRP-CO, we
introduce the Greedy Insertion Regret, which can be regarded as hybrid between the common
Greedy Insertion and K-Regret Insertion operators.
Preliminary experiments were performed on three other commonly used methods. The Random
Insertion operator, which positions customers in random order at random locations, was not able to
provide good solutions for the 2E-VRP-CO. The Basic Greedy and the K-Regret heuristic, introduced by Ropke and Pisinger (2006), extend the computation time significantly without meaningful improvements in solution quality. Therefore, these three operators are omitted from our ALNS heuristic.
We provide a detailed description of the three insertion operators used in our ALNS heuristic: 1. Greedy Insertion
The Greedy Insertion operator inserts all customers in random order at the position that minimizes the insertion cost.
2. Greedy Insertion Perturbed
introduction of noise in the selection step allows for additional solutions to be explored. Based on Hemmelmayr et al. (2012) we apply a perturbation factor on the interval [1 − τ, 1 + τ ] to the insertion costs.
3. Greedy Insertion Regret
The Greedy Insertion Regret is a new operator, which can be regarded as a hybrid between the Greedy Insertion operator and the well known K-Regret Insertion operator. This operator deviates from the Greedy Insertion operator in that the customers are not inserted in a random order. This operator prioritizes the customers who are most likely difficult to insert in the current routes, as does the K-Regret Insertion operator. If the best insertion location is no longer feasible for these customers the alternatives could be very expensive. For each customer the priority is proportional to the distance to the nearest customer which is already routed. The customers who are located furthest away have a higher probability to be routed first. The customers within range of a covering location are always assigned the lowest priority, since they can always be serviced via that specific covering location. This operator is different from the K-Regret Insertion operator in that the insertion costs of all customers to insert is not recalculated each time a customer has been inserted, thereby the computation time is reduced significantly.
After the second echelon routes are reconstructed, we recreate the first echelon routes. At first, we determine if the current first echelon routes are still feasible. If that is not the case we reconstruct the first echelon routes from the ground up. This is preferred over adapting the current first echelon routes, since the changes in the second echelon routes can have a substantial effect on the used intermediate facilities. One of the essential characteristic of the 2E-VRP-CO is that the intermediate facilities can be supplied by multiple trucks. Otherwise, no feasible solutions would be possible if the demanded quantity at an intermediate location would be larger than the capacity of a truck and no other intermediate facilities are available. We implement a preprocessing step used by Breunig et al. (2016) and Wang et al. (2017). For any intermediate location with a requested demand larger than the truck capacity, pseudo duplicate locations with a demand equal to the truck capacity are created, until the remaining demand is smaller than a full truckload. For all pseudo locations, direct back-and-forth trips from the depot are created. Routes for the remaining intermediate facilities are created using the same insertion operators as mentioned above. For each insertion operator the possibility exists that there is no feasible insertion location for one of the customers, due to the order in which the customers are reinserted into the solution. In this case the insertion phase is restarted with a new customer insertion order until a feasible solution is found. The likelihood of the occurrence of such a situation is dependent on the amount and capacity of the available city freighters.
3.4
Acceptance criteria
A new solution s0 is accepted based on three measures. Firstly, s0 is always accepted if it is an
improvement over the current solution, s. Secondly, s0is always accepted if a large destroy operator
is used in its creation, in that case the grace period counter igraceis equal to zero. Finally, we accept
some of the solutions with a worse objective value than the current solution based on a simulated annealing acceptance criterion (Kirkpatrick et al., 1983). The advantage of this method is local optima can be departed, by accepting solutions which do not result in improvements immediately. Therefore, a larger part of the solution space can be examined. The probability of acceptance is equal to e−(f (s0)−f (s))/T
the ALNS heuristic the temperature is initialized at T = Tstart. At each iteration the temperature
decreases by a factor κ ∈ [0, 1], denoted as the cooling rate.
3.5
Local search
A local search is performed on all new solutions accepted by the acceptance criteria, described in Section 3.4. First the second echelon routes are improved and afterwards the first echelon routes are revised. For both echelons the same local search operators are used. The local search on the second echelon routes is performed for each satellite location separately to prevent extreme first echelon route changes. Therefore, the requested demand at each satellite location and covering location remains unchanged and the first echelon routes are not impacted.
The local search consists of the following operators, which are performed sequentially:
1-0-Exchange, Intra-Swap, Intra-2-Opt, Inter-Swap and Inter-2-Opt. The 1-0-Exchange operator removes
each customer in a random order and reinserts the customer at the best location in one of the routes departing from the same satellite location. Afterwards two intra route operators are employed. The Intra-Swap operator changes the position of two customers within the same route. The
Intra-2-Opt removes two arcs and reconnects the remaining fragments of the routes optimally.
Two similar operators are also performed between 2 different routes originating from the same location, namely the Inter-Swap and Inter-2-Opt operators. The Inter-Swap operator changes the positions from two customers within different routes. The Inter-2-Opt is similar to the Intra-2-Opt operator, instead the arcs are removed from different routes. Each operator is executed once and we select the change that provides the largest improvement in objective value, without violating any of the constraints.
3.6
Updating operator probabilities
To select the destroy and repair operators, we use a roulette wheel selection principle with weights assigned to all operators. The probability of selecting an operator j is determined by weight of parameter j devided by the sum of the weights of all operators in the same type as operator j. We update the weights by using the same procedure as Pisinger and Ropke (2007). After a segment of 100 iterations the operators weights are updated based on a weighted average of their past value and the performance score of the operator over the segment. The performance is dependent on the amount of times an operator is selected and the score accumulated over the segment, nj and
ψjfor operator j. The score is accumulated by three methods: (1) if the new solution is better
than the best global solution found, the score increases by σ1, (2) if the new solution is better
than the current solution, the score increases by σ2and (3) if the new solution is worse than the
current solution but accepted by the simulated annealing acceptance criterion, the score increases by σ3. In all other cases the score remains unchanged. Hereby, we reward the operators that find
new improved solutions and operators that diversify the search to some extent. Both the destroy and repair operators obtain the same score increase, since we cannot distinguish the effect of the two operators on the new solution. Finally, the weights are updated by the procedure
θi← (1 − λ)θi+ λ
ψi
ni
, (38)
where θiis the weight of operator i and λ is the reaction factor, which controls how fast the weight
4
Computational experiments
This section presents the computational experiments conducted on well known benchmark instances for the 2E-VRP by Breunig et al. (2016), the instrances for the SFL-VRP by Veenstra et al. (2018) and on newly created 2E-VRP-CO instances. The mathematical model, described in Section 2, is solved using CPLEX 12.8 and the ALNS heuristic is programmed in R. The code was executed single threaded on one core. Section 4.1 provides descriptions of the sets of instances considered for the three different problems mentioned above. We describe the parameter setting of the ALNS heuristic in Section 4.2. In Section 4.3 we describe the performance of the ALNS heuristic on the benchmark instances for the 2E-VRP and the SFL-VRP. Section 4.4 reports the results of the heuristic on the 2E-VRP-CO instances in comparison to the branch-and-bound algorithm, presented in Section 2. A sensitivity analysis on the used destroy and repair operators is presented in Section 4.5. Finally, managerial insights are given in Section 4.6.
4.1
Description of the benchmark sets
In this section we describe the available sets of benchmark instances for the 2E-VRP and the SFL-VRP, which are special cases of the 2E-VRP-CO as described in Section 2.2. Furthermore we explain how new instances specifically for the 2E-VRP-CO are created from these sets. A schematic overview of the characteristics of the benchmark sets is given in Appendix B. It describes the number of customers (|C|), covering locations (|L|), satellite locations (|S|), trucks (|T |) and city freighters (|CF |) for the instances in the benchmark sets. As described in Section 2.2, the locker and patient routes for the SFL-VRP by Veenstra et al. (2018) can be considered as truck and city freighter routes respectively.
For the 2E-VRP we consider three well known sets of benchmark instances, namely ‘Set 2’, ‘Set 3’ and ‘Set 5’. The first two sets are adaptions proposed by Perboli et al. (2011) on some of the Christofides and Eilon instances for the VRP. Each instance contains 2 or 4 satellite locations, which are positioned at the coordinates of randomly selected customers. In both sets instances with 22, 31 and 50 customers are available. Unfortunately, two different versions of the instances of set 2 and 3 with 50 customers are in circulation, due to inconsistencies with respect to the satellite locations and nomenclature. Consequently, these different versions are referred to as subset ‘b’ and ‘c’ to identify the different versions in circulation. Optimal results for ‘Set 2’ and ‘Set 3’ are reported by Baldacci et al. (2013). The ‘Set 5’ instances are introduced by Hemmelmayr et al. (2012). This set contains instances with up to 200 customers, the largest instances available for the 2E-VRP. So far these instances have proven a difficult challenge and only 3 of the 18 best known solutions are proven optimal. In our analysis we use the 2E-VRP benchmark instances as provided by Breunig et al. (2016)2.
For the SFL-VRP 117 randomly generated instances are available. The amount of customers within these instances ranges from 30 to 150. However, we only used the first 72 instances for our comparison analysis, since for the other instances no best known solutions are reported. For each configuration of customers there are 9 instances with up to 50 covering locations, referred to as ‘lockers’ by Veenstra et al. (2018). The customers and lockers are positioned randomly on a square. The depot is placed at a random location near the centre of the square. The patient routes are penalized by a factor 10 relative to the locker routes. For further details concerning the exact configuration of these instances we refer to Veenstra et al. (2018).
Finally we present new sets of instances specifically created for the 2E-VRP-CO, since this problem has not been studied before. We construct the sets by expanding the three sets of 2E-VRP benchmark instances with additional covering locations. For each instance we introduce a number of covering locations equal to the amount of satellite locations. Each covering location is positioned at the same location as a randomly selected customer, with the restraint that the customer is not within the covering range of another covering location. This positioning method is selected because of four reasons: (1) it ensures that for each positive covering distance there is always at least one customer within the covering range, (2) it increases the likelihood that a covering location is placed near a group of clustered customers, (3) it governs that the covering locations are spread over the entire region and (4) it provides parallels to the placement method of satellite locations for the ‘Set 2’ and ‘Set 3’ instances for the 2E-VRP. In case there is no feasible location left to place a covering location, since all customers are already within range of another covering location, we place the remaining covering locations at a random position on the customer plane. The customer plane is the smallest square that contains all customers. Thus the customer plane covers all x-coordinates between xminand xmaxand y-coordinates between and yminand ymax. Where xminand xmaxare
the smallest and largest x-coordinates of all customers and likewise yminand ymax are the smallest
and largest y-coordinates of all customers. With ddiagwe denote the length of the diagonal line
between points [xmin, ymin]and [xmax, ymax]. The covering distance is identical for all covering
locations and is equal to 25% of ddiag. The connection cost, incurred if a customer i is serviced by
a covering location j, is `ji= αdji+ βddiag, where α is the distance dependent connection cost
factor and β a fixed connection cost factor. For all instances we set α = β = 0.25. Customers who are outside the covering range of a covering location cannot be serviced by that location. The corresponding connection costs are therefore set to infinity. No changes are made to the demand of each customer, the vehicle characteristics and the coordinates of the depot, satellite locations and customers. Therefore, these values are as reported for the corresponding 2E-VRP instances.
4.2
ALNS heuristic parameter tuning
This section provides details on the settings used for the parameters of the ALNS heuristic. The parameters are tuned on the 72 instances of ‘Set 2’, ‘Set 3’ and ‘Set 5’ of the 2E-VRP, the 72 random instances by Veenstra et al. (2018) and the 72 newly created instances for the 2E-VRP-CO. Special attention is paid to the effect of the parameters on the problems with 100 or more customers. We started the parameter tuning by considering the parameter settings described by Hemmelmayr et al. (2012). Subsequently, we iteratively tested a range of values for each parameter, only accepting improvements.
The final values of the heuristic parameters are as follows. At the start of the heuristic 75% of the satellite locations and 30% of the covering locations are opened at random. At each iteration a random amount of customers between q = max(0.2|C|, 5) and ¯q = min(0.4|C|, 40)
are removed from the current solution. The grace period, ωgrace, is equal to 250 iterations. The
current best solution after 10.000 iterations without a new global improvement. At that moment the operator selection probabilities are also reset. As a trade-off between solution quality and computation time the algorithm ends after 100.000 iterations.
4.3
Comperative Results on the 2E-VRP and the SFL-VRP instances
To validate the quality of our proposed heuristic we perform benchmark tests on the 2E-VRP and SFL-VRP instances. The results of the benchmark tests are provided in Tables 2 and 3. Our obtained results are compared to the best known solutions (BKS) for the instances within the literature. The best known solutions that are proven optimal are represented with an asterisk. For each instance the best and average solution found over 5 runs by our ALNS heuristic are reported, as well as the gap in percentages with the best known solution in the literature, the computation time, t(s), and the average time until the optimal solution is found, t∗(s). We note that computation speed may
not be competitive to the current heuristics proposed for these problems, since our heuristic is composed for a more generalized problem and is programmed in R.
Table 2 reports results of our ALNS heuristic on the 2E-VRP instances. Over all instances we find an average gap of 0.42% to the best known solution. For the instances of ‘Set 2’ and ‘Set 3’ our heuristic finds all the best known solutions. Over these sets our heuristic is on average 5 times slower then the heuristic proposed by Breunig et al. (2016). For the instances of ‘Set 5’ we increase the total amount of iterations to 250.000 and we perform restarts from the current best solution found after 15.000 iterations without finding a new improved solution. For the instances of ‘Set 5’ our heuristic finds the best known solution of one of the instances and we report an average gap of 1.7%. Over the instances of ‘Set 5’ our heuristic is on average 3 times slower then the heuristic proposed by Breunig et al. (2016).
To be able to solve the SFL-VRP we make a few adjustments to our ALNS heuristic. To take the duration constraints into account, the repair and local search phases are extended with additional validity checks. Any change to a solution which would result in a violation of the duration constraints is not allowed. Furthermore, multiple destroy operators are removed since they do not provide any benefits within the SFL-VRP setting. Considering that there is only one satellite location, positioned at the depot, the Satellite Removal, Satellite Opening and Satellite Swap operators are removed. The Minimum Quantity Removal operator is excluded, since no information on the amount of parcels required by each customers is available. Moreover, the Random Route
Removal and Least Used Vehicle Removal operators are discharged since we observed that for all but
Table 2: Results for the 2E-VRP benchmark instances of Breunig et al. (2016).
Inst. |C| |S| |T | |CF | BKS Our ALNS heuristic
Avg. 5 Best Gap BKS (%) t(s) t∗(s) Set 2a E-n22-k4-s6-17 21 2 3 4 417.07∗ 417.07 417.07 0.00 250 0 E-n22-k4-s8-14 21 2 3 4 384.96∗ 384.96 384.96 0.00 258 1 E-n22-k4-s9-19 21 2 3 4 470.60∗ 470.60 470.60 0.00 234 2 E-n22-k4-s10-14 21 2 3 4 371.50∗ 371.50 371.50 0.00 229 0 E-n22-k4-s11-12 21 2 3 4 427.22∗ 427.22 427.22 0.00 229 1 E-n22-k4-s12-16 21 2 3 4 392.78∗ 392.78 392.78 0.00 262 2 E-n33-k4-s1-9 32 2 3 4 730.16∗ 730.16 730.16 0.00 342 1 E-n33-k4-s2-13 32 2 3 4 714.63∗ 714.63 714.63 0.00 307 9 E-n33-k4-s3-17 32 2 3 4 707.48∗ 707.48 707.48 0.00 316 1 E-n33-k4-s4-5 32 2 3 4 778.74∗ 778.74 778.74 0.00 323 59 E-n33-k4-s7-25 32 2 3 4 756.85∗ 756.85 756.85 0.00 313 2 E-n33-k4-s14-22 32 2 3 4 779.05∗ 779.05 779.05 0.00 292 12 Set 2b E-n51-k5-s2-4-17-46 50 4 4 5 530.76∗ 530.76 530.76 0.00 351 162 E-n51-k5-s2-17 50 2 3 5 597.49∗ 597.49 597.49 0.00 323 29 E-n51-k5-s4-46 50 2 3 5 530.76∗ 530.76 530.76 0.00 346 27 E-n51-k5-s6-12 50 2 3 5 554.81∗ 554.81 554.81 0.00 310 78 E-n51-k5-s6-12-32-37 50 4 4 5 531.92∗ 531.92 531.92 0.00 314 55 E-n51-k5-s11-19 50 2 3 5 581.64∗ 582.94 581.64 0.00 331 12 E-n51-k5-s11-19-27-47 50 4 4 5 527.63∗ 527.63 527.63 0.00 314 93 E-n51-k5-s27-47 50 2 3 5 538.22∗ 538.22 538.22 0.00 360 8 E-n51-k5-s32-37 50 2 3 5 552.28∗ 552.28 552.28 0.00 330 3 Set 2c E-n51-k5-s2-4-17-46 50 2 3 5 601.39∗ 601.39 601.39 0.00 293 128 E-n51-k5-s2-17 50 2 3 5 601.39∗ 601.39 601.39 0.00 304 65 E-n51-k5-s4-46 50 2 3 5 702.33∗ 703.49 702.33 0.00 321 76 E-n51-k5-s6-12 50 2 3 5 567.42∗ 567.42 567.42 0.00 320 56 E-n51-k5-s6-12-32-37 50 4 4 5 567.42∗ 671.26 567.42 0.00 310 150 E-n51-k5-s11-19 50 2 3 5 617.42∗ 617.42 617.42 0.00 334 17 E-n51-k5-s11-19-27-47 50 4 4 5 530.76∗ 530.76 530.76 0.00 326 50 E-n51-k5-s27-47 50 2 3 5 530.76∗ 530.76 530.76 0.00 345 6 E-n51-k5-s32-37 50 2 3 5 752.59∗ 752.59 752.59 0.00 271 143 Set 3a E-n22-k4-s13-14 21 2 3 4 526.15∗ 526.15 526.15 0.00 255 2 E-n22-k4-s13-16 21 2 3 4 521.09∗ 521.09 521.09 0.00 244 2 E-n22-k4-s13-17 21 2 3 4 496.38∗ 496.38 496.38 0.00 237 1 E-n22-k4-s14-19 21 2 3 4 498.80∗ 498.80 498.80 0.00 227 5 E-n22-k4-s17-19 21 2 3 4 512.80∗ 512.80 512.80 0.00 219 2 E-n22-k4-s19-21 21 2 3 4 520.42∗ 520.42 520.42 0.00 198 3 E-n33-k4-s16-22 32 2 3 4 672.17∗ 672.30 672.17 0.00 286 86 E-n33-k4-s16-24 32 2 3 4 666.02∗ 666.02 666.02 0.00 289 22 E-n33-k4-s19-26 32 2 3 4 680.36∗ 680.36 680.36 0.00 306 75 E-n33-k4-s22-26 32 2 3 4 680.37∗ 680.37 680.37 0.00 301 76 E-n33-k4-s24-28 32 2 3 4 670.43∗ 670.43 670.43 0.00 290 20 E-n33-k4-s25-28 32 2 3 4 650.58∗ 650.58 650.58 0.00 290 6 Set 3b & 3c E-n51-k5-s12-18 50 2 3 5 690.59∗ 690.92 690.59 0.00 333 117 E-n51-k5-s12-41 50 2 3 5 683.05∗ 690.05 683.05 0.00 313 93 E-n51-k5-s12-43 50 2 3 5 710.41∗ 710.41 710.41 0.00 340 18 E-n51-k5-s13-19 50 2 3 5 560.73∗ 560.73 560.73 0.00 314 13 E-n51-k5-s13-42 50 2 3 5 564.45∗ 564.45 564.45 0.00 294 21
Table 2 – Continued from previous page
Inst. |C| |S| |T | |CF | BKS Our ALNS heuristic
Avg. 5 Best Gap BKS (%) t(s) t∗(s)
E-n51-k5-s13-44 50 2 3 5 564.45∗ 564.45 564.45 0.00 342 16 E-n51-k5-s39-41 50 2 3 5 728.54∗ 730.17 728.54 0.00 299 128 E-n51-k5-s40-41 50 2 3 5 723.75∗ 728.82 723.75 0.00 286 135 E-n51-k5-s40-42 50 2 3 5 746.31∗ 746.31 746.31 0.00 294 84 E-n51-k5-s40-43 50 2 3 5 752.15∗ 752.34 752.15 0.00 300 50 E-n51-k5-s41-42 50 2 3 5 771.56∗ 772.50 771.56 0.00 290 132 E-n51-k5-s41-44 50 2 3 5 802.91∗ 807.35 802.91 0.00 296 114 Set 5 100-5-1 100 5 5 32 1564.46∗ 1622.82 1611.94 3.03 1331 863 100-5-1b 100 5 5 15 1108.62 1129.26 1119.64 0.99 1402 1082 100-5-2 100 5 5 32 1016.32∗ 1027.07 1023.35 0.69 1473 1020 100-5-2b 100 5 5 15 782.25 792.01 782.25 0.00 1570 1130 100-5-3 100 5 5 30 1045.29∗ 1050.44 1048.46 0.30 1365 1067 100-5-3b 100 5 5 16 828.54 830.88 828.99 0.05 1638 1022 100-10-1 100 10 5 35 1124.93 1153.73 1136.89 1.06 1727 1515 100-10-1b 100 10 5 18 916.25 929.57 925.02 0.96 1924 1090 100-10-2 100 10 5 35 990.58 1022.56 1019.26 2.90 1518 1254 100-10-2b 100 10 5 18 768.61 783.45 781.36 1.66 1769 1271 100-10-3 100 10 5 32 1043.25 1077.68 1059.07 1.52 1606 1145 100-10-3b 100 10 5 17 850.92 871.45 865.31 1.69 1756 1067 200-10-1 100 10 5 62 1556.79 1624.00 1618.67 3.97 5751 4925 200-10-1b 100 10 5 30 1187.62 1244.06 1231.65 3.71 4282 3588 200-10-2 100 10 5 63 1365.74 1430.58 1411.22 3.33 6324 5528 200-10-2b 100 10 5 30 1002.85 1032.11 1021.97 1.91 4626 4398 200-10-3 100 10 5 63 1787.73 1827.98 1814.10 1.48 5694 5161 200-10-3b 100 10 5 30 1197.90 1230.16 1213.97 1.34 4031 3224 Avg. Set 2 578.27 581.81 578.27 0.00 305 41 Avg. Set 3 632.30 636.08 632.30 0.00 285 51 Avg. Set 5 1118.81 1148.88 1139.62 1.70 2766 2242 Avg. 734.46 743.73 739.66 0.42 914 595
Table 3: Results for the SFL-VRP benchmark instances of Veenstra et al. (2018).
Inst. |C| |L| Veenstra et al. (2018) BKS Our ALNS heuristic
Avg. 5 Best t(s) Avg. 5 Best Gap Veenstra (%) Gap BKS (%) t(s) t∗(s)
R1 30 10 3884.0 3884 10 3884∗ 3884.0 3884 0.00 0.00 390 18 R2 30 15 2897.0 2897 11 2897∗ 2897.0 2897 0.00 0.00 335 7 R3 30 20 2050.0 2050 11 2050∗ 2050.0 2050 0.00 0.00 339 7 R4 30 25 3496.0 3496 12 3496∗ 3496.0 3496 0.00 0.00 330 41 R5 30 30 3841.0 3841 15 3841∗ 3849.0 3841 0.00 0.00 397 72 R6 30 35 2767.0 2767 17 2767∗ 2767.0 2767 0.00 0.00 352 35 R7 30 40 2129.0 2129 16 2129∗ 2137.6 2129 0.00 0.00 363 75 R8 30 45 1012.4 1012 17 1012∗ 1024.0 1012 0.00 0.00 351 184 R9 30 50 455.0 455 18 455∗ 457.4 455 0.00 0.00 334 97 R10 40 10 4388.0 4388 20 4388∗ 4388.0 4388 0.00 0.00 520 16 R11 40 15 3874.0 3874 19 3874∗ 3874.0 3874 0.00 0.00 413 30 R12 40 20 3531.0 3531 17 3531∗ 3531.0 3531 0.00 0.00 383 56 R13 40 25 4495.0 4495 22 4495∗ 4495.0 4495 0.00 0.00 464 56
Table 3 – Continued from previous page
Inst. |C| |L| Veenstra et al. (2018) BKS Our ALNS heuristic
Avg. 5 Best t(s) Avg. 5 Best Gap Veenstra (%) Gap BKS (%) t(s) t∗(s)
R14 40 30 4399.0 4399 26 4399∗ 4399.0 4399 0.00 0.00 461 39 R15 40 35 3964.0 3964 27 3964∗ 3970.4 3964 0.00 0.00 481 124 R16 40 40 2789.0 2789 28 2789∗ 2801.4 2789 0.00 0.00 435 128 R17 40 45 2641.0 2641 25 2641∗ 2642.8 2641 0.00 0.00 369 169 R18 40 50 3390.0 3390 30 3390∗ 3390.2 3390 0.00 0.00 458 206 R19 50 10 5259.0 5259 31 5259∗ 5259.0 5259 0.00 0.00 524 18 R20 50 15 4768.0 4768 34 4768∗ 4768.0 4768 0.00 0.00 481 92 R21 50 20 5058.0 5058 32 5058∗ 5058.0 5058 0.00 0.00 552 50 R22 50 25 2351.0 2351 31 2351∗ 2360.8 2351 0.00 0.00 376 106 R23 50 30 4040.0 4040 38 4040∗ 4040.0 4040 0.00 0.00 512 62 R24 50 35 3959.0 3959 38 3959∗ 3959.0 3959 0.00 0.00 530 76 R25 50 40 3697.0 3697 36 3697∗ 3700.2 3697 0.00 0.00 511 106 R26 50 45 4120.0 4120 43 4120∗ 4120.0 4120 0.00 0.00 509 151 R27 50 50 2715.2 2715 38 2715∗ 2738.2 2721 0.22 0.22 453 211 R28 60 10 6314.5 6245 59 6170∗ 6170.0 6170 -1.20 0.00 860 8 R29 60 15 4289.0 4289 43 4289∗ 4289.0 4289 0.00 0.00 531 37 R30 60 20 5245.0 5245 54 5245∗ 5259.8 5245 0.00 0.00 586 55 R31 60 25 2540.0 2540 41 2540∗ 2540.0 2540 0.00 0.00 442 125 R32 60 30 3854.0 3854 49 3854∗ 3854.0 3854 0.00 0.00 517 88 R33 60 35 2960.0 2960 45 2960∗ 2975.6 2960 0.00 0.00 429 176 R34 60 40 4830.0 4830 72 4830∗ 5605.0 4830 0.00 0.00 591 90 R35 60 45 3562.0 3562 55 3562∗ 2569.6 3562 0.00 0.00 513 174 R36 60 50 2677.0 2677 48 2677∗ 2680.0 2677 0.00 0.00 387 155 R37 70 10 5605.0 5605 82 5605∗ 5605.0 5605 0.00 0.00 713 26 R38 70 15 5157.0 5157 73 5157∗ 5157.0 5157 0.00 0.00 517 25 R39 70 20 5447.0 5447 76 5447∗ 5447.0 5447 0.00 0.00 583 130 R40 70 25 4738.0 4738 118 4738∗ 4738.0 4738 0.00 0.00 538 111 R41 70 30 5083.0 5083 93 5083∗ 5083.0 5083 0.00 0.00 597 214 R42 70 35 1732.0 1732 104 1732∗ 1732.0 1732 0.00 0.00 408 96 R43 70 40 4053.0 4053 97 4053 4059.4 4053 0.00 0.00 518 236 R44 70 45 3365.0 3365 91 3365 3375.0 3365 0.00 0.00 540 265 R45 70 50 2287.0 2286 13 2286 2321.2 2286 0.00 0.00 491 251 R46 80 10 6244.0 6244 108 6244∗ 6244.0 6244 0.00 0.00 772 26 R47 80 15 5162.0 5162 99 5162 5162.0 5162 0.00 0.00 544 76 R48 80 20 5997.0 5997 119 5997 5997.0 5997 0.00 0.00 739 200 R49 80 25 4323.0 4323 93 4323 4323.0 4323 0.00 0.00 566 179 R50 80 30 4155.0 4155 104 4145∗ 4155.0 4155 0.00 0.24 518 84 R51 80 35 6181.0 6181 97 6181 6181.0 6181 0.00 0.00 729 386 R52 80 40 4065.0 4065 91 4065 4070.2 4065 0.00 0.00 630 276 R53 80 45 2424.0 2424 123 2424∗ 2440.0 2424 0.00 0.00 521 295 R54 80 50 3623.0 3623 108 3623 3700.2 3623 0.00 0.00 710 552 R55 90 10 7453.0 7453 173 7453 7502.2 7453 0.00 0.00 1005 383 R56 90 15 5919.0 5919 139 5919 5919.0 5919 0.00 0.00 714 83 R57 90 20 4668.0 4668 111 4668 4668.0 4668 0.00 0.00 494 65 R58 90 25 6003.0 6003 137 6003 6003.0 6003 0.00 0.00 711 148 R59 90 30 4638.0 4638 131 4638 4640.4 4638 0.00 0.00 617 305 R60 90 35 4867.0 4867 160 4867 4866.0 4862 -0.10 -0.10 701 156 R61 90 40 3681.0 3681 110 3681 3683.4 3681 0.00 0.00 575 206 R62 90 45 2553.0 2553 134 2553 2568.0 2553 0.00 0.00 510 368 R63 90 50 1586.0 1586 125 1586∗ 1608.2 1600 0.88 0.88 495 312 R64 100 10 7034.0 7034 219 7034 7050.0 7034 0.00 0.00 1004 203
