A heuristic for the Location Routing Problem with Stochastic
Demand and a maximum workday recourse
L.M. van der Heide
s2341026
May 22, 2018
Abstract
The Location Routing Problem with Stochastic Demand (LRPSD) concerns the problem where depot locations and vehicle routes need to be determined simultaneously, while the uncertain customer demands only become known upon arrival. This necessitates recourse actions where a vehicle is replenished reactively or preventively at the depot. We introduce a recourse cost for routes exceeding the workday time. This means that the amount of time a route takes, including travelled distance and time serving customers exceeds the time of a workday. We design an Adaptive Large Neighbourhood Search algorithm using especially defined operators to solve the LRPSD. Its quality is benchmarked against a lowerbound and instances from the literature for the LRPSD without the workday recourse. We find that the ALNS yields competitive results compared to the literature, even finding 8 new best known results. Moreover, when the workday recourse is incorporated, we find solutions that are more easily implemented. The performance of the workday recourse cost is also compared to a constraint on the maximum customer demands served by one vehicle. We find that the incorporation of the workday recourse cost yields solutions resulting in shorter distances travelled per vehicle and less depot capacity violations.
1 Introduction 2
1
Introduction
We consider the problem of a company that needs to send products to customers. It, therefore, needs to simultaneously determine a set of open depots and vehicle routes, or sequence in which customers are visited. This is, however, complicated by the uncertainty of customer demands that only become known at arrival of a vehicle. This means that routes are scheduled before the customers’ demands are certain. One vehicle is used for each route and all vehicles have the same maximum inventory. Since a customer’s demand only becomes known upon arrival of a vehicle, it is possible that its demand exceeds the remaining vehicle inventory. This is called a route failure, which causes a reactive replenishment. When a reactive replenishment occurs, the vehicle must return to the depot to replenish and revisit the same customer afterwards. After a visit to a customer, it is also possible to schedule a preventive replenishment, which means that the vehicle could replenish at the depot and then continue to the next customer in the route. This could prevent a route failure. The total routing costs are defined as the total expected distance travelled, including distances travelled due to reactive and preventive replenishments plus the fixed cost of using a vehicle. The routes scheduled before demand is known are denoted a priori routes. They differ from the actual routes since even though the sequence of customer visits is known, the vehicle may replenish reactively or preventively at the depot between customers.
Due to the preventive and reactive replenishments, the time needed to complete a route is unknown for each vehicle. However, the driver of such a vehicle should only be scheduled to work for a certain number of hours each day. Therefore, an additional overtime fee is paid when the route time, including service time at the customers, exceeds the workday duration. This balances the trade-off between of potentially long routes and overtime cost or many drivers on short routes. Furthermore, the depots may have insufficient capacity to fulfil the customers’ demands. When the depot capacity is exceeded, additional products can be provided from another source at an additional cost. The three costs that are generated due to the unknown customer demand are called recourse costs.
The combination of determining depot locations and corresponding routes with uncertain customer de-mand is called the Location Routing Problem with Stochastic Dede-mand (LRPSD) (Marinakis et al., 2016), a combination of the Location Routing Problem (LRP), in which deterministic demand is considered and the Vehicle Routing Problem with Stochastic Demand (VRPSD), for which only one depot is considered with unknown or stochastic customer demands. As the LRP and the VRPSD are NP-hard, this also holds for their combination, the LRPSD. The LRP and the VRPSD have both been researched extensively in the literature. For an overview of LRP literature we refer to Nagy and Salhi (2007) or more recently Ritzinger et al. (2016). Overviews of the VRPSD are provided in Oyola et al. (2016a) and in Braekers et al. (2016) and of its solution methods in Oyola et al. (2016b). In Gendreau et al. (2016) future research direction for the VRPSD are mentioned. The field of stochastic vehicle routing is still in the early stages and significant research on new models for stochastic vehicle routing is possible and necessary.
In this paper, we consider the LRPSD where customer demands only become known upon arrival in which depot locations and a priori vehicle routes need to be determined simultaneously. Due to the unknown customer demands, we include preventive and reactive recourse actions for route failures and incorporate a recourse cost for routes exceeding the workday time. Additionally, recourse actions for exceeding depot capacity are implemented. We find solutions consisting of opened depots and routes using a newly defined specialisation of the Adaptive Large Neighbourhood Search (ALNS) heuristic, which is benchmarked on special cases from the literature and optimal solutions for small instances.
1 Introduction 3
are considered is that of Huang (2015). It presents a multi-compartment pick-up and delivery version of the LRPSD with multiple products. Since it considers separate fleets for pick-up and delivery, the routing for each can be determined separately. Similar to our setting, they consider a cost per unit when the depot capacity is exceeded. In Klibi et al. (2010) the location transportation problem with stochastic demand is considered. Here, instead of scheduling routes, another company is hired to make the deliveries. The costs of these deliveries depend on the type, size and location of the demands and not on the distances between customers and the depot. Again only reactive replenishments are allowed. Currently, workday recourse actions as considered in our paper are not implemented in the literature. Our first assumption is that a vehicle can be replenished preventively as well as reactively. In the literature, it is often assumed that replenishments are only implemented when a route failure occurs (Huang, 2015), Klibi et al. (2010). This is a simplification from reality, where information on the current inventory of a vehicle is available at all times. When preventive replenishments are allowed without a restriction on route lengths, Yang et al. (2000) show that the best vehicle allocation is one vehicle per depot. The resulting solution is then impractical to implement, since drivers cannot work shifts that are too long. Hence, adding a maximum duration constraint yields a more practical and difficult problem. Our second assumption is an additional recourse cost when the time necessary to implement the realisation of the a priori route exceeds a workday. Another method to ensure routes are not too long is to assume that the probability of a route failure is smaller than a certain threshold. This is a much stronger version of our assumption that a cost is incurred if a route exceeds a certain length. This assumption is made by Huang (2015), who use it to determine the necessary number of routes. In practice, the probability of a route failure is only part of what influences the time a route takes. Therefore, the routes created using this assumption might be unnecessarily short. We, therefore, also include the time necessary for serving customers and the total distance travelled.
Our third assumption is that for demand exceeding the maximum depot capacity an additional price is paid. In current literature it is instead commonly assumed that the sum of the maximum inventories for vehicles departing from a depot is smaller than the depot capacity. Then the demand should never exceed the depot capacity, as long as the maximum vehicle inventories are not exceeded. Examples of papers with this assumption are Marinakis et al. (2016) and Marinakis (2015). However, in the problem considered in this paper this does not apply due to the unknown customer demands when the a priori routes are scheduled. This means that the vehicle capacities can be exceeded in the actual route due to the stochastic demand. Therefore, an assumption on the sum of vehicle capacities does not ensure that there is enough capacity at the depot. We add a recourse cost of exceeding depot capacity to incorporate this risk in the objective function.
Our fourth assumption is that customer demands only become known upon arrival at the customer. This has been implemented in both the VRPSD literature Braekers et al. (2016) and the LRPSD literature Huang (2015). In other stochastic LRP settings it is not the quantity of product demanded that is unknown, it is the occurrence of a demand that is unknown. This is also called the Location Routing Problem with Stochastic Customers (LRPSC), as in Albareda-Sambola et al. (2007). In the LRPSC a posteriori routes are determined after the occurrence of demand becomes known, whereas we only find these by implementing the a priori routes and implementing replenishments. Stochastic demand is also modelled using fuzzy demands in Nadizadeh and Nasab (2014) and Mehrjerdi and Nadizadeh (2013), among others. Finally, demand scenarios are implemented to model the unknown demand in Mete and Zabinsky (2010), among others. Both approaches simplify calculations, but they result in less precise costs. Another source of stochasticity in the LRP could be the travel time of routes. This is considered in, for instance, Fazel-Zarandi et al. (2013). Even though customer demand is then known the travel times are uncertain, which could occur due to traffic jams.
2 Problem Description 4
3 we translate the problem to a mathematical model. Then in Section 4 we describe the ALNS we use to solve this model. In Section 5 we compare the performance of our algorithm to the current literature. Finally, in Section 6 we give a conclusion and summarise our findings.
2
Problem Description
The LRPSD is defined on a set I = {1, . . . , m} of potential depot locations, each with a maximum capacity Qi and opening cost Oi for i ∈ I. Binary variables yi indicate an open depot if they equal one
and a closed depot if they equal zero. We consider the set J = {m + 1, . . . , m + n} of customers that need to be visited. These customers j are served by depot i ∈ I, which is determined using the variable zij
which is one if customer j is served by depot i and zero otherwise. The number of customers considered is n. We also consider the set of nodes V = I ∪ J , that is all customers and depots. Then one a priori route is scheduled for each vehicle to satisfy the customers’ demands. We assume that a route starts and ends at the same depot and that all customers are visited. There is a set K = {1, . . . , K} of identical vehicles with maximum inventory level L. This means that at most L units can be delivered before the vehicle needs to be replenished. The fixed cost for using a vehicle is F . F is calculated using the indicator variable ak which equals 1 if vehicle k ∈ K is scheduled and 0 otherwise. The route lengths can
be calculated using the distances between pairs of customers, and customers and depots dij, for i, j ∈ V.
We assume that a customer is served by one route only, so no split deliveries occur where multiple vehicles each deliver part of the demand. Then the variables xijk indicate the sequence of depots and customers
for all routes and are one if on route k ∈ K the vehicle travels from node i ∈ V to node j ∈ V.
The random demand of customer j is denoted by ξj and only becomes known at the arrival of a vehicle.
We assume that ξj are independently and identically distributed with a known discrete distribution.
Moreover, we assume that ξj ∈ L = {0, . . . , L} meaning a customer demand never exceeds the maximum
vehicle inventory. Realisations of this variable are denoted by Djω, where ω ∈ Ω is the scenario or possible
combination of demands for all of the customers. Ω is then the set of all possible scenarios. Then the inventory in the vehicle between customers j ∈ J and u ∈ J in scenario ω ∈ Ω for vehicle k ∈ K is denoted by Ijukω
Since demands are unknown when the routes and depots are determined, we need to schedule a priori routes. These are predetermined routes, scheduled before the customers’ demands become known. In the remainder of this paper we will use ‘routes’ to denote a priori routes. Since the realized sum of demands in a vehicle route may exceed the maximum vehicle inventory, the routes may fail when they are implemented. Therefore, we calculate the cost of a route as the expected route length, where reactive and preventive replenishments are implemented when route failures would occur. When the sum of demands realisations of the customers allocated to a depot exceeds its capacity, a depot failure occurs. Finally, a workday failure occurs when the sum of the distances and the number of customers visited by a vehicle multiplied by the customer service time exceeds the workday duration.
For each of these failures there are different recourse actions and corresponding costs. The recourse action for a depot failure is having an external party deliver additional products, which is modelled as an incurred cost per unit c as in Huang (2015). The cost of this recourse action is calculated using the variable qωi, which is the amount of customer demands exceeding the depot capacity of depot i ∈ I in
scenario ω ∈ Ω.
2 Problem Description 5
serving customers. We then define overtime as the amount of time this exceeds W , the workday length. The overtime fee is then w per unit of overtime. The variable vωk is then the amount of overtime for
each scenario ω ∈ Ω and route k ∈ K.
Finally, the recourse action for a route failure is to replenish the vehicle at the depot, as in Albareda-Sambola et al. (2007) and Marinakis et al. (2016). This reactive replenishment incurs recourse costs, which are the routing costs of going from the customer where this failure occurred to the depot and back to the customer to continue the route. We set the indicator variable rujω equal to 1 when a reactive
replenishment occurs between customers u and j ∈ J for demand scenario ω ∈ Ω and 0 otherwise. The additional route length due to a reactive replenishment between customer u ∈ J and customer j ∈ J in route k ∈ K for scenario ω ∈ Ω is then denoted by hujωk. This reactive replenishment is scheduled
due to a route failure at customer j, however, for modelling simplification the corresponding variables are defined between two customers. We also allow for preventive replenishments, meaning that after each visit to a customer the vehicle can replenish at the depot and continue the route from there. The indicator variable pujω is equal to 1 if a preventive replenishment is scheduled between customers u ∈ J
and j ∈ J in scenario ω ∈ Ω and is equal to 0 otherwise. The additional route length due to a preventive replenishment being scheduled between customer u ∈ J and customer j ∈ J in route k ∈ K for scenario ω ∈ Ω is denoted by gujωk. Finally, we define Rujω as the total quantity replenished due to either a
preventive or reactive replenishment between customers u ∈ J and j ∈ J in demand scenario ω ∈ Ω. Similar preventive and reactive replenishments are also included in Marinakis et al. (2016), Marinakis (2015) and Bianchi et al. (2006). Since preventive replenishments are only allowed when they are expected to be less costly than continuing the route with the current inventory, we need to calculate the remaining expected route costs. The function f (lj) calculates the cost of continuing the route from customer j with inventory l. Since in the remainder of this route preventive and reactive replenishments may occur we need to determine the replenishments between customers i and j in the remainder of the route, given the demand distribution of customer j. For this we define the indicator variable sljb, which equals one
if there is a preventive replenishment scheduled before customer j ∈ J who demands b ∈ L and the inventory before visiting the customer is l ∈ L.
The LRPSD then consists of three interacting parts: determining the set of open depots, allocating customers to the depots and determining routes to all customers, while minimizing total (expected) costs. The necessary sets, parameters, functions and variables are presented in Tables 1-2.
Tab. 1: Sets Sets Explanation
I All nodes corresponding to depots {1,. . . , m}
J All nodes corresponding to customers {m+1,. . . , m + n} V All nodes corresponding to either depots or customers, I ∪ J K All possible vehicles {1, . . . , K}
Ω All possible demand scenarios
L All possible vehicle inventory or demand levels {0, . . . , L} Tab. 2: Function Function Explanation
3 Mathematical model 6
Tab. 3: Parameters Parameter Explanation
m Number of depots n Number of customers
K Maximum number of vehicles
L Maximum inventory level of a vehicle Oi Opening costs of depot i ∈ I
dij distance from i to j for i, j ∈ V
F Fixed cost for using a vehicle
M A sufficiently big value: 2 maxi∈I,j∈Jdij
c Recourse cost per unit for exceeding depot capacity w Overtime fee per hour for exceeding workday length
Djω Realisation of demand for customer j ∈ J in scenario ω ∈ Ω
Qi Maximum capacity of depot i ∈ I
t Time necessary to serve a customer s Speed of vehicles in time per distance W Workday length
ξj Demand for customer j ∈ J , follows known discrete probability distribution function
Tab. 4: Variables Variable Explanation
yi Binary variable, equals 1 if depot i ∈ I is open, 0 otherwise
xijk Binary variable, equals 1 if route k ∈ K goes from customer or depot i to j for i, j ∈ V
ak Binary variable, equals 1 if route k ∈ K is scheduled, 0 otherwise
qωi The amount of depot capacity exceeded in scenario ω ∈ Ω for depot i ∈ I
vωk The amount of workday length exceeded for scenario ω ∈ Ω for route k ∈ K
gujωk The cost of the preventive replenishment between u ∈ V and j ∈ J in demand scenario
ω ∈ Ω for vehicle k ∈ K
hujωk The costs of the reactive replenishment between u ∈ V and j ∈ J in demand scenario
ω ∈ Ω for vehicle k ∈ K
zij Binary variable, equals 1 if customer j ∈ J is served by depot i ∈ I
Iijω The inventory carried from node i to j, i, j ∈ V in scenario ω ∈ Ω
Rujω The amount replenished due to a reactive or preventive replenishment between u ∈ J and
j ∈ J in scenario ω ∈ Ω
rujω Equals 1 if a route failure occurs between customers u ∈ J and j ∈ J in scenario ω ∈ Ω
pujω Equals 1 if a preventive replenishment is scheduled between customers u ∈ V and j ∈ J
in scenario ω ∈ Ω
sljb Equals 1 if a preventive replenishment is used before customer j ∈ J who demands b ∈ L
if the inventory is l ∈ L sufficient inventory
3
Mathematical model
3.1
Non-linear model
3 Mathematical model 7
constraints are introduced to specify the LRPSD.
The objective function of the non-linear mathematical model for the LRPSD is as follows: minX i∈I Oiyi+ X k∈K X u∈V X v∈V xuvkduv+ F X k∈K ak +X ω∈Ω cX i∈I qωi+ w X k∈K vωk+ X u∈J X j∈J X k∈K (gujωk+ hujωk) P (ω) (1)
In (1) the cost corresponding to a solution is calculated. First, the opening cost and the costs of the route are calculated. Then the depot failure costs are determined and the overtime fee is calculated for all possible demand combinations of the customers. The preventive and reactive replenishment costs are then added for each possible demand scenario. The total of these recourse costs is multiplied by the probability of this demand scenario occurring, P (ω). These costs are then minimized subject to the following constraints: X k∈K X i∈V xijk= 1 ∀j ∈ J (2) X k∈K X i∈V xjik= 1 ∀j ∈ J (3) X i∈I zij = 1 ∀j ∈ J (4) zij ≤ yi ∀i ∈ I, j ∈ J (5) ak ≥ xijk ∀k ∈ K, i, j ∈ V (6) X i,j∈S xijk≤ |S| − 1 ∀S ⊆ V, S 6= ∅, k ∈ K (7) X i∈I X j∈J xijk≤ 1 ∀k ∈ K (8) X u∈V
(xiuk+ xujk) ≤ 1 + zij ∀i ∈ I, j ∈ V, k ∈ K (9)
X j∈V xjuk= X j∈V xujk ∀k ∈ K, u ∈ V (10) qωi≥ X j∈J Djωzij− Qiyi ∀i ∈ I, ω ∈ Ω (11) vωk+ W + t ≥ X u∈V X v∈V (xuvk(t + sduv)) + s X i∈J X j∈J (gijωk+ hijωk) ∀k ∈ K, ω ∈ Ω (12) gujωk≥ X i∈I
zij(dui+ dij− duj) − M (1 − pujω) − M (1 − xujk) ∀u, j ∈ J , ω ∈ Ω (13)
hujωk≥
X
i∈I
zij(dji+ dij) − M (1 − rujω) − M (1 − xujk) ∀u, j ∈ J , ω ∈ Ω (14)
Iujω ≤ L X k∈K xujk ∀u ∈ V, j ∈ J , ω ∈ Ω (15) Iujω ≤ X v∈V Ivuω− Duω X k∈K
xujk+ Rujω ∀u, j ∈ J , ω ∈ Ω (16)
Iujω ≥ Djω
X
k∈K
3 Mathematical model 8
rujω+ pujω≤
X
k∈K
xujk ∀u, j ∈ J , ω ∈ Ω (18)
Rujω≤ (rujω+ pujω)L ∀u, j ∈ J , ω ∈ Ω (19)
pujω(f (L, j) + dui+ dij) < duj+ f (Iujω− rujωRujω, j) ∀u, j ∈ J , ω ∈ Ω (20)
f (l, j) =X k∈K X u∈J xjuk X b∈L h P (l ≥ b)(1 − sljb)(dju+ f (l − b, u)) + sljb( X i∈I zij(dji+ diu) + f (L, u))+ P (l < b) X i∈I zij(dji+ dij) + (1 − sljb)(dju+ f (L − b + l, j)) + sljb( X i∈I (dji+ diu) + f (L, u)) i +X k∈K X i∈I xjikdji ∀j ∈ J , l ∈ L (21) xuvk∈ {0, 1} ∀u, v ∈ V, k ∈ K (22) gujωk, hujωk∈ R+ u, j ∈ J , ω ∈ Ω, k ∈ K (23)
pujω, rujω∈ {0, 1} ∀u, j ∈ J , ω ∈ Ω (24)
yi, zij∈ {0, 1} ∀i ∈ I, j ∈ J (25)
Rujω, Ivjω ∈ L ∀v ∈ V, u, j ∈ J , k ∈ K, ω ∈ Ω (26)
sljb∈ {0, 1} ∀l, b ∈ L, j ∈ J (27)
qωi, vωk∈ R+0 ∀k ∈ K, ω ∈ Ω, i ∈ I. (28)
Constraints (2)– (10) are constraints that are commonly used in LP formulations for the LRP. Constraints (2) ensure that all customers are visited exactly once and Constraints (3) ensure that only one vehicle leaves each customer. Constraints (4) ensure all customers are allocated to a depot, whereas (5) ensure that customers can only be allocated to open depots. Moreover, Constraints (6) ensure that routes can only be used if the vehicle is scheduled. Constraints (7) are included to prevent subtours. Constraints (8) ensure that all edges from one customer to another can be used only once. Then constraints (9) are implemented to ensure that a customer can only be assigned to a depot if there is a route leaving from that depot that reaches the customer. Constraints (10) ensure the flow conservation, meaning the same route enters and leaves a node. The following constraints are newly designed to formulate the mathematical problem of the LRPSD as it is considered in this paper. Constraints (11) and (12) ensure proper calculation of, respectively, the depot failure recourse cost and the recourse cost of exceeding the work day. Constraints (11) ensure qωi is at least equal to the amount of demand that exceeds the depot
capacity or zero. Constraints (12) ensure that vωk is at least equal to the amount of time necessary for
the route of vehicle k, including the time spent at the customer, minus the work day length. Constraints (13) determine the costs of preventive replenishments. The costs are equal to returning to the depot and continuing the route from the next customer, if a preventive replenishment is scheduled. For its calculation we need M , which is a sufficiently big number to ensure the cost of a replenishment is only added when it is implemented. For this problem, we set it to 2 maxi∈I,j∈Jdij. Constraints (14) determine
4 Solution algorithm 9
responsive replenishments. Together with the constraint on the maximum of the inventory they ensure a correct replenishment quantity. Constraints (20) ensure that preventive replenishments can only occur when they seem profitable with all the information known when customer j is visited. This means that the costs of a preventive replenishment are less than of continuing the route. The cost of a preventive replenishment is the distance to the depot and continuing the route at the next customer with L inventory. The cost of continuing the route is the cost of going to the next customer with the current inventory without preventive replenishment. The current inventory between customers u and j is Iujω− rujωRujω,
since we want to compare with the inventory if there is no preventive replenishment. Finally constraint (21) defines the function that determines the expected costs of continuing the route from customer j onwards with inventory I, including possible (preventive) replenishments. First, we determine the next customer in the route, customer u. Then we define two cases, either the inventory is sufficient for customer j or a reactive replenishment will take place. If the inventory is sufficient, we still need to determine if we want to preventively replenish after visiting customer j. If the inventory l is insufficient, we know the inventory is replenished to L − b + l, since the remaining demand b was subtracted from L. Then again afterwards we need to decide on either continuing the route or replenishing preventively. Since the decision for a preventive replenishment is made when the demand of customer j is known, the variables indicating the expected preventive recourse sljb depend on the demand b.
Lastly, constraints (22)–(28) show the requirements on the used variables.
3.2
LRPSD Relaxation
To linearise the mathematical model for the LRPSD, we have to relax the assumption that customer demands become known at arrival. This assumption is relaxed so that all customer demands become known at the start of a route. Therefore, the a priori routes are scheduled before customer demands become known and hence route, depot and workday failures may still occur. However, since the reali-sations of demands become known before the reactive and preventive replenishments are scheduled, we can optimally determine when to replenish. We no longer need to specify rules on when a preventive replenishment can take place. Therefore, we can remove constraints (20) and (21) from the mathematical model. The resulting LRPSD relaxation yields a lower bound to the LRPSD.
4
Solution algorithm
To solve the mathematical problem that was defined in Section 3, we use an Adaptive Large Neighbour-hood Search (ALNS) heuristic. The classical ALNS heuristic is an iterative search and was first defined by Ropke and Pisinger (2006a) and further developed in Pisinger and Ropke (2007). It has been imple-mented for problems related to the LRPSD such as the VRP, the VRPSD and the LRP. Examples of successful implementations of the ALNS for the VRPSD are Lei et al. (2012), Luo et al. (2016) and others. Since the LRPSD consists of two parts: the LRP problem and the VRPSD we also consider papers on the LRP that use ALNS. This has been done by Hemmelmayr (2015), Hemmelmayr et al. (2012) Ko¸c (2016), Tunalıo˘glu et al. (2016) and Yang and Sun (2015). Since the ALNS has been implemented successfully on both halves of the considered problem, it appears a viable option for the LRPSD as well.
4 Solution algorithm 10
part of the partial solution. Then a repair operator repairs the partial solution by inserting the customers from the customer pool into the partial solution, which then becomes a complete solution.
The general outline of our ALNS is shown in Algorithm 1. Let x denote the matrix with entries equal to xijk, or x = (xijk) ∈ {0, 1}(m+n)×(m+n)×K. Similarly, let z denote the matrix with entries equal to
zij, or z = (zij) ∈ {0, 1}m×n. Then (x, z) yields a solution to the LRPSD, since the values of all other
variables can be determined using these. Then G(x, z) is the cost corresponding to the solution (x, z) of the LRPSD. An alternative calculation of G(x, z) is presented in Section 4.1.
Algorithm 1 Adaptive Large Neighbourhood Search outline
1: Construct initial solution (x, z)
2: Set G∗= G(x, z) and (x∗, z∗) = (x, z)
3: Set b = 0, br= 0, bT = 0, bτ = 0
4: while b < B, b − Bb < Bn, time < MT do 5: Determine q
6: Use w−i to determine o−i
7: Add q customers to customer pool using o−i
8: Use w+i to determine o+i 9: Implement o+i to form (x0, z0) 10: if G(x0, z0) < G∗ then 11: Set G∗= G(x0, z0) 12: Set (x∗, z∗) = (x0, z0) 13: Set Bb = b 14: bτ = 0 15: else 16: if (x’,z’) is accepted then 17: Set (x, z) = (x0, z0) 18: bτ = 0 19: Update s−i and s+i 20: if br= ρ then
21: Update wi− and w+i using s−i and s+i
22: Set br= 0 23: if bτ = τ then 24: Reset wi+, wi−, si+ and s−i 25: if bτ = 2τ then 26: Set (x, z) = (x∗, z∗) 27: if bT = BT then 28: Update T 29: Set bT = 0 30: Increment b, br, bT, bτ 31: end
4 Solution algorithm 11
lowest cost we found so far and we start the next iteration with this solution. If G(x0, z0) > G∗ we use a simulated annealing acceptance criterion dependent on a temperature T to determine if we accept this solution and use it in the next iteration. If it is not accepted, we start the next iteration from the original solution (x, z). The simulated annealing acceptance criterion is explained in Section 4.3. Next, the performance of the destroy and repair operators is used to update their scores, s−i and s+i respectively. These scores are used to update the operator weights every ρ iterations. If no new solutions are accepted for τ iterations, the scores and weights are reset to 0 and 1 respectively. If no new solutions are accepted for 2τ iterations, we continue with the solution yielding the lowest costs so far. Then, the temperature used for the simulated annealing acceptance criterion is updated every BT iterations.
We end the iteration by incrementing b, br, bT and bτ. Finally, the ALNS is finished when no solutions
yielding a lower cost have been found for Bn iterations, when the maximum number of iterations B is
reached or the time limit MT is exceeded.
In the next sections, the different aspects of the ALNS are explained in more detail.
4.1
Objective
To solve the model shown in Section 3 using the ALNS, we redefine the objective function. This enables us to calculate the same costs using recursion with a non-linear objective. The recursion lies in the routing cost calculation, as will be explained later on. Moreover, it ensures constraints (11) - (21) are no longer necessary.
For the formulation of the non-linear objective, we introduce some new notation. We use the same formulation as in Bianchi et al. (2006), with a slightly altered notation for clarification. We let C(x) denote the recourse function, which is the stochastic part of the objective function. Then, we define fk,u(l) as the expected length of route k from the node u onwards with an initial inventory level of l. The
calculation of fk,u(l) is shown below. We also define fk,u(l)c and fk,u(l)p as the expected route length
of route k from node u onwards with initial inventory level l when the a priori route is continued or a preventive replenishment is scheduled after node u, respectively. We let φ(j) denote the customer in the route following customer j and θ(k) denote the final customer in route k. Moreover, ψ(k) indicates the depot corresponding to vehicle k. An overview of the functions used in the objective function is presented in Table 4.1.
Function Explanation G(x, z) Objective function C(x) Recourse function
fkj(l) Function indicating the expected route length from customer or depot j in route k with
an inventory level l.
fk,u(l)c Function indicating the expected route length from customer or depot j in route k with
an inventory level l when no preventive recourse is scheduled.
fk,u(l)p Function indicating the expected route length from customer or depot j in route k with
an inventory level l when a preventive recourse is scheduled. φ(j) Function indicating the customer in the route following customer j ψ(k) Function indicating to which depot route k is connected
θ(k) Function indicating the final customer in route k
4 Solution algorithm 12
Similar to the original objective function, we calculate the depot opening costs and the vehicle costs. The vehicle costs are now calculated as the number of times a vehicle leaves a depot multiplied by the fixed cost of using a vehicle. Finally, C(x) is the recourse function which can be defined as follows:
C(x) = E " wX k∈K sfk,ψ(k)(L) + t X u∈V X j∈J xujk− W − t + # + E " cX i∈I X j∈J (ξjzij) − Qiyi + # +X k∈K fk,ψ(k)(L) (30)
The recourse function consists of three separate parts. The first part concerns the workday failure. We calculate the expected number of hours a route takes by multiplying the speed with the route length and adding the number of customers multiplied by t. If this total exceeds the number of hours in a workday W , a recourse cost per exceeded hour is added to the objective function. We do not need to add the distances travelled in the a priori route separately, as in the original formulation, since these are included in the function fk,i(l). Note that the expected value of the complete workday failure sum is calculated
and not of only the routes.
The second part is the cost of exceeding the depot capacity. Here we add up all the demands of customers allocated to the same depot and subtract the depot capacity. If the resulting value is positive, we multiply it with the costs of exceeding the depot capacity.
The third part concerns the routing costs. Since the recourse actions of a route failure are preventive or reactive replenishments at the depot, the routing costs include these additional distances. Let l be the remaining load after completion of service at customer j. Then the full calculation of fk,ψ(k)(L) is as
follows: fk,j(l) =min {fk,jc (l), f p k,j} (31) fk,jc (l) =dj,φ(j)+ P (ξφ(j)≤ l)fk,φ(j)(l − ξφ(j)) + P (ξφj> l) 2dφ(j),ψ(k)+ fk,φ(j)(l + L − ξφ(j)) (32) fk,jp =dj,ψ(k)+ dψ(k),φ(j)+ fk,φ(j)(L − ξφ(j)) (33) fk,θ(k)(l) =dθ(k),ψ(k), l ∈ L (34)
In equation (31) the preventive replenishing is introduced by choosing the least costly option: continuing the a priori route or returning to the depot to replenish preventively. Then the expected route length of continuing the route are presented in (32). It consists of the distance to the next customer in the route and the costs of either continuing the route from there, or returning to the depot for a reactive replenishment and continuing the route. In (33) the expected costs of preventive replenishments are calculated by adding the distance to the depot and from the depot to the next customer, after which the expected routing costs for all following customers are added. Finally, in (34) the cost of continuing the route after the final customer, meaning returning to the depot, are given. Using these equations, the expected routing costs of the remainder of the route can be calculated recursively. Starting from the final customer we find the routing cost, then using this we calculate the costs starting from the penultimate customer.
4.2
Initial solution
4 Solution algorithm 13
exceeds the sum of expected demands of all customers. Then we use one of the repair operators called greedy insertion, which is explained in Section 4.5, to design an initial solution. Similarly, in Ko¸c (2016) and Tunalıo˘glu et al. (2016) the initial solution is created using greedy insertion. In Hemmelmayr (2015) the total number of depots to be opened for the initial solution is a random number between the total demand divided by the average depot capacity and the total number of depots. This allows for more flexibility in the initial solution. However, by randomly opening depots we also achieve this, however, generally we will have fewer initially opened depots. They then design the initial solution by assigning customers to their nearest depot and construct routes using a savings algorithm. However, using greedy insertion is less time consuming, especially for sizeable experiments.
4.3
Acceptance criterion and stopping criterion
If a solution is accepted, or used as a basis for the next iteration, depends on a Simulated Annealing (SA) temperature acceptance criterion. This SA acceptance criterion was also implemented in Ropke and Pisinger (2006b), Pisinger and Ropke (2007) and others. The probability of accepting solution (x, z) is: e−(G(x,z)−G(x∗,z∗))/T, where T is the temperature used and (x∗, z∗) is the solution yielding the smallest
G(x, z) that has been found. The initial temperature is set to 2G(x1, z1), where (x1, z1) denotes the
initial solution. T is multiplied by α every BT iterations during the heuristic.
Finally, the ALNS is stopped after the maximum number of iterations B is reached, or when there has been no reduction in costs for Bn iterations.
4.4
Adaptive operator weights
Which destroy and repair operators are used is determined by two separate roulette wheel mechanisms, where the probability of choosing operator oiis based on its past performance through the operator weight
wi. A similar technique is used in, among others, Ropke and Pisinger (2006a). After each iteration, the
“score” si of the used operator i is increased by σ1 if the resulting solution is the best so far, σ2 if the
solution is accepted and has a lower cost than the previous accepted solution and by σ3 if the solution
is accepted but has a higher cost than the previous. Then, every segment, or ρ iterations, the weights of the operators are updated based on their scores and the scores are set to zero again. The weights are updating using the weight sensitivity parameter ω as follows: wit= (1 − ω)wit−1+ ωsi
Ni. Where w
t
i denotes
the updated weight, wt−1i denotes the old weight of operator i and Ni is the number of times operator
i was implemented in the past segment. We use the weight sensitivity parameter ω to determine how sensitive the operator weight is to the operator score of the past segment.
These operator weights are then used to find the probability of choosing an operator in the following way. Let O the set of repair (destroy) operators. Then the probability of choosing repair (destroy) operator oi is given by Pwi
i∈O
wi. Hence, the probability of choosing an operator with a high score is increased.
4 Solution algorithm 14
4.5
Destroy and repair operators
The ALNS searches a different solution neighbourhood every iteration by first destroying the current solution to the LRPSD using a destroy operator and then repairing it using a repair operator to find an improved solution. In the literature several operators appear. From the papers on VRPSD or LRP using ALNS we have selected common operators and implemented them. Moreover, we designed some destroy operators specifically for our problem. First, we explain the destroy operators in more detail.
4.5.1 Destroy operators
Destroy operators are used to destroy the current solution. This means that q customers are removed from the solution and added to the customer pool. If customer j is removed from solution (x, z), xuvk= 1,
for k ∈ K and u ∈ V such that xujk = 1 and for v ∈ V such that xjvk= 1. Then all other xujk and xjuk
are set to 0. Moreover, all zij = 0 for i ∈ I. This results in a partial solution where customer j is not
serviced. This customer pool is then used by the repair operator, which inserts the customers from the customer pool into the destroyed solution. The value of q is determined by randomly selecting a value from the Uniform distribution between the bounds {1, πn}. Where π is a parameter that determines the maximum percentage of customers added to the customer pool. We consider the following destroy operators.
Random removal starts by randomly selecting q customers and putting them in the customer pool. All n customers are equally likely to be selected. This operator is used in almost all papers using ALNS. Worst removal determines which q customers have the highest removal gain and adds those to the customer pool. We use a perturbation ϑ(j) which is random in [−10, 10] for each customer j to versify the search to prevent remaining in local optima. We define the removal gain is defined as the difference in total costs when the customer is included in the solution or not. Let G0(x, z, j) be defined as the cost of a solution where customer j has been removed from the solution defined by (x, z). Then the removal gain is defined as ∆G−j= G(x, z) − G0(x, z, j). We add the perturbation ϑ to this and then remove the
q customers with the highest ∆Gj+ ϑ(j) from the solution.
Related removal is started by determining a random seed customer j. Then j and the q − 1 most related customers are removed. Relatedness is defined based on nearness in distance, differences in expected demand and current location in the routes. The level of relatedness between customers u and j is called rujand calculated as follows: ruj = duj+ |λu− λj| + |eu− ej| + |Ku− Kj| + |ou− oj|. We define
λj as the mean demand of customer j and ej= {i : zij = 1} or the depot serving customer j. Moreover,
Kj is the route by which customer j is visited; in this route customer j is the oj-th to be visited. This
same definition has been used in the literature as well, for instance in Hemmelmayr et al. (2012). The Close depot operator is started by randomly determining which depot i to close (or keep closed) by setting the corresponding yito zero. Then, we add all customers that were served from this depot, i.e.
all j ∈ J for which zij = 1, to the customer pool. If none were served, none are added to the customer
pool. If all depots are closed afterwards, a random depot u is chosen and yu is set to 1.
The Depot swap operator starts, similar to current literature, with a random closed depot i and uses a roulette wheel selection of another depot to be opened. Let νvbe the probability of opening depot v ∈ I.
The probabilities of the roulette wheel depend on the distance to the closed depot. Hence, νv =Pdiv
u∈Idiu.
The customers j for which zij= 1 are added to the customer pool. At this point no customers are served
by the newly opened depot.
4 Solution algorithm 15
we add the q nearest customers to the customer pool. We define the nearness of customers j to depot i by the distance to the depot dij. Again, we have an open depot that does not serve any customers.
The Depot feasibility removal is an operator that has not been used in the literature before and was designed for the specific problem considered here. Firstly, for each depot i the probability of exceeding the depot capacity ζi is calculated as follows: ζi = P
P
j∈J
(ξjzij) > Qi
. Then the depots are ranked from highest to lowest ζi. We then start with the depot with the highest ζi and randomly determine
the number of customers between 0 and q to remove and randomly select this many customers to add to the customer pool. Then we move to the depot which is the next most likely to fail and use the same procedure. This is repeated until q customers have been removed or all open depots are considered. It is possible that less than q customers are added to the customer pool after this operator.
The Vehicle feasibility removal operator has also been designed especially for the current problem. It operates similarly to the previous one, however, it considers the probability of a route failure. Firstly, for each route k the probability of a route failure is determined as follows: κk = P
P u∈V P j∈J xujkξj > L . Then the routes are ranked from highest κk to lowest. Then, starting from the most likely to fail route,
we remove a random number of customer using the uniform distribution from 0 to q and add them to the customer pool. We continue doing so at the next most likely to fail route until q customers are added to the customer pool or there are no more routes.
The Workday feasibility removal operator has also been designed especially for the current problem. In addition to operating similarly to the previous operator, it also considers the probability of a workday failure. Firstly, for each route the probability of exceeding the workday is determined as follows: χk =
P sfk,ψ(k)(L) + t P u∈V P j∈J xujk− t > W !
. Then starting from the highest χkroute, we remove a random
number of customers using the uniform distribution from 0 to q and add them to the customer pool. We continue doing so at the next highest χk route until q customers are removed or there are no more routes.
4.5.2 Repair operators
Repair operators are used to insert the customers from the customer pool into the destroyed solution again. Let G(x, z, u, v, j) denote the cost of the solution where customer j is added in the solution represented by (x, z) between nodes u ∈ V and v ∈ J ∪ {u}. There are three different situations when customer j is inserted between nodes u and v. Firstly, j is inserted between two customers, or after a customer and before a depot: u ∈ J and v ∈ V. Then the solution changes as follows: zij = ziu, for
all i ∈ I. Moreover, for k ∈ K such that xuvk = 1 we set xujk and xjvk 1, while xuvk is now set to 0.
Secondly, j can be inserted after a depot and before a customer: u ∈ I and v ∈ J . Then, zij = ziv,
for all i ∈ I. Again, for k ∈ K such that xuvk = 1 we set xujk and xjvk to 1, while xuvk is now set to
0. Thirdly, customer j can be inserted between two depots v = u ∈ I. This means that customer j is visited by a new route, meaning zuj= 1 and zij = 0, for all i ∈ I/{u}. Moreover, xujk and xjuk are set
to 1 for k = arg min
k∈K
P
u∈U
P
v∈V
xuvk = 0. This means that the smallest value for k, for which no route yet
exists. We now consider the following repair operators that insert the customers from the customer pool into the solution.
Greedy insertion consists of randomizing the order of customers in the customer pool and inserting them in the least costly position. Let G(x, z, j)∗ = min
u∈V,v∈J ∪{u}
5 Experiments 16
and Pisinger (2006a), one of the first ALNS papers, customers are ordered based on their insertion costs. However, the proposed version of the operator is more recently implemented in the literature and is less time consuming. In Yang and Sun (2015) the insertion costs are defined as the distance to its predecessor and successor. However, since this would influence the probability of a route failure, we cannot use this same approximation. In Ko¸c (2016) and Ropke and Pisinger (2006a) a version of this operator is included, where a perturbation is added to the objective function. This would induce randomness and diversification. However, since we add customers in a random order, this is already incorporated. Regret insertion starts by determining for each customer in the customer pool the regret values ∆G(x, z, j) which determine the order in which customers are inserted in the least costly position. Then, again G(x, z, j)∗ = min
u∈V,v∈J ∪{u}G(x, z, u, v, j) denotes the cost of the solution where customer j is
in-serted in the least costly position in the route. Similarly, we define G(x, z, j)2 as the cost of the solution
where customer j is inserted in the second least costly position. Then the regretvalue of customer j is ∆G(x, z, j) = G(x, z, j)2− G(x, z, j)∗. Then we consider customers in the order from biggest to smallest
regret value and again insert them in the solution in the position where G(x, z, u, v, j) is minimized. This position might be different from the one considered in the order, since other customers are previously added. We now implement the same principle as in greedy insertion, only adding customers in this pre-specified order. In Ropke and Pisinger (2006a) only different routes are considered, while Hemmelmayr et al. (2012) use the same approach as we do and consider all locations. In the literature different levels of regret insertion are used. For instance the regret-3 operator reinserts every removed customer with the largest regret value first. This regretvalue is now determined by looking at the difference between the best and the 3’rd best insertion. Since determining the best insertion position is time consuming, we do not use this version of regret insertion.
Random insertion inserts the customers from the customer pool in a random location. For a randomly chosen customer j in the customer pool, a random depot i and two nodes u and v in route k for which xuvk = 1 are chosen. Then customer j in inserted into the solution between nodes u and v. Similarly,
all other customers in the customer pool are inserted randomly in the solution (x, z). This insertion method does not compare the costs corresponding to different locations and is widely implemented in the literature.
5
Experiments
5 Experiments 17
capacity as an upper bound on the customer demands, these instances are too time consuming to solve with our assumptions.
In order to gain meaningful results, the following parameters were tuned using instance Prins23: π, the upper bound on q the number of customers to add to the customer pool ,ρ the number of iterations between weight updates, τ when to reset the weights, the decrease in SA temperature α, the SA temper-ature is updated every BT iterations, the scores σ
1, σ2, σ3, the weight sensitivity ω. We tried different
settings and found that the most suitable combination for the LRPSD is: {π, ρ, τ, α, σ1, σ2, σ3, ω, BT} =
{0.4, 20, 100, 0.5, 2000, 20, 10, 5, 0.1}. We set R = 25000, MT equal to 20 hours and quit the heuristic
if no solution is accepted for more than Bn = 5000 iterations. All destroy and repair operators used
in the ALNS were also tested by excluding them one at a time and comparing the solution quality and speed. We found that is was beneficial to incorporate all operators. Moreover, the values of the operator weights were printed every ρ iterations. Generally, we found that the operator weights decreased as the number of iterations increased. This means that fewer solutions were accepted later on in the algorithm. The ALNS therefore, converges to an optimum over time. For the destroy operators, we found that generally the depot related operators had smaller weights than the other operators. The newly designed destroy operators, the vehicle feasible removal, day feasible removal and the depot feasible removal, all performed similarly to the standard destroy operators such as the worst removal and related removal. Moreover, we found that the random insertion operator had overall a lower weight than the greedy and regret insertion operators. However, at the end of the algorithm its weight did increase, whereas the other weights decreased. Therefore, we included it in our ALNS.
Since the calculation of the objective function for large numbers of customers and large vehicle capac-ities is time consuming, an approximation of the objective function is used during the ALNS for the larger instances. Normally, the customer demand has as upper bound the vehicle capacity. When the approximation is used, we do not use demands that have a probability of less than 0.01% in the recursive calculation of the expected route length. Moreover, the value for π is set to 0.2, which means fewer customers are added to the customer pool by destroy operators. We verify the objective approximation by comparing the solution costs of the ALNS with the objective approximation to the normal objective function. We find both the speed of iterations and the solution quality were improved for large instances. Therefore, this approximation is implemented for all Tuzun instances, and the Prins instances with 200 customers or 100 customers and maximum vehicle inventory of 150.
To validate our ALNS heuristic we compare its results with two different benchmarks. First, we compare it to a lowerbound found by relaxing the assumption that customer demands only become known upon arrival. Second, we compare the results of the ALNS for a special case of our problem to the best known results in the literature found by another heuristic presented in Marinakis et al. (2016). Then we calculate the value of the stochastic solution for the LRPSD and compare our new workday assumption to one in the literature.
5.1
Lowerbound
5 Experiments 18
cost for the depot failure and workday failure are not considered.
We cannot solve this model for all possible demand scenarios using the distribution functions of al customer demands due to the computational complexity. Therefore, we use Sample average approximation (SAA) based on a Monte Carlo approximation to solve the linear program to optimality using Cplex. The SAA has also been applied in, among others, Kleywegt et al. (2002). When at most one route is considered, the model contains 1048575 subtour elimination constraints which are added as lazy constraints to the solver. This means that these constraints are only added when they are violated by a solution found by the solver. Due to the computational complexity, we only consider 50 scenarios. This might result in a solution that is not optimal when the complete distributions of all demands are considered. Moreover, we need to specify the maximum number of routes used in the solution. We let the maximum number of routes be equal to the number of routes found when the ALNS was used to solve this problem. Since we consider the assumptions from Marinakis et al. (2016), the maximum number of routes for these instances is 1.
We compare the cost of the solution for the LRPSD found by the ALNS to the optimal solution of the relaxation found using the SAA to show the quality of our heuristic. However, we can only solve small problems to optimality. In Table 5 the costs for the Prins instances with 20 customers are shown. The epsilon that is reported is the percentage difference between the best integer found by the solution, which is reported in the table as lower bound, and the best bound found the Cplex. We find that the ALNS finds solutions comparable to the lower bound. The difference between the costs found by the ALNS and those found by the SAA can be explained by three things: the quality of the ALNS, the relaxation of the unknown demand assumption and the approximation of the distribution of demand by the SAA scenarios. If the ALNS is not well-tuned, the final solutions created are significantly costlier than an optimal solution. The relaxation of the unknown demand should ensure that the cost found using the SAA is a lowerbound to the LRPSD solved using the ALNS. However, since it is a lower bound, an optimal solution to the LRPSD will have a greater cost than an optimal solution to the relaxation. Finally, the SAA is only an approximation. If the scenarios included in the model are not representative of the actual demand distributions the SAA is not effective. For instance, if all demand realisations are lower than the average customer demand the cost found by the SAA is lower than that of an optimal solution. Conversely, if the scenarios are pessimistically chosen, meaning all demands are greater than the average customer demands, the cost found by the SAA is greater than that of an optimal solution. This complicates the comparison of the lowerbound and the ALNS cost. An additional complication is that the relative difference between the best integer objective and the objective of the best node remaining for the instance Prins1 is relatively big, since Cplex was stopped due to exceeding the memory capacity of 46Gb. This means that the cost of an optimal solution might be lower than the one reported here. Tab. 5: Objective function for the SAA and heuristic solution for the Prins et al. benchmark we created.
Instance Lowerbound Epsilon ALNS cost ∆% Prins1 40757.81 41.79% 41927 2.87% Prins2 31867.89 25.62% 30468 -4.30% Prins3 34803.54 36.81% 34100 -2.02% Prins4 30455.11 22.18% 29500 -3.14%
5.2
Literature benchmark
5 Experiments 19
vehicle capacities can be exceeded as is shown in Yang et al. (2000). The best known results (BKR) of the LRPSD with preventive and reactive replenishments in the literature are presented in Marinakis et al. (2016). The costs of our heuristic and the best known results for the Prins and Tuzun instances are presented in Table 6. We define ∆% = 100ALN S−BKR
BKR , where ALN S denotes the cost of the ALNS
and BKR denotes the cost of the BKR. The runs of the ALNS for Prins took 30 minutes for small instances and were terminated at 20 hours for the larger instances. The runtimes of the benchmark were not reported. The Tuzun dataset consists of instances with more customers and, therefore, the ALNS was terminated after 20 hours for every instance.
When the Prins dataset is considered, the ALNS found a solution that is significantly less costly than the BKR reported in Marinakis et al. (2016) for all instances. However, the reduction in cost is larger than would be expected and is possibly due to an unreported assumption made by Marinakis et al.. This assumption might be related to exceeding the depot capacity, since the Tuzun instance contains uncapacitated depots and for that dataset we do find similar results. Since the number of opened depots is not presented in Marinakis et al. (2016), we cannot confirm this suspicion. This unknown assumption makes it impossible to compare the solutions found by the Marinakis benchmark and our ALNS. In Section 8.1 of the Appendix, we show the solution quality of two different assumptions that could have been made by Marinakis et al.. However, these also result in costs significantly different from the BKR. Since the costs found by the relaxation in the previous section are similar to the ones found by the ALNS, we know that for the assumptions presented in Marinakis and included in our model the ALNS finds suitable solutions.
5 Experiments 20
Tab. 6: Comparison ALNS and BKR presented in Marinakis et al. (2016)
5 Experiments 21
5.3
LRPSD with exceeding workday recourse cost
We now present the costs of the newly defined LRPSD with recourse costs for a workday failure. Both travelling time and time at a customer are included in the workday time in this model. This makes it more flexible to different routing situations and more realistic than only including travelling time or the number of customers. For the calculation of our depot and workday exceeding recourse costs, we need to define some parameters. We set the recourse cost per unit for exceeding depot capacity c equal to 150. The time necessary to serve a customer t equal to 10, the time necessary to travel equal to the distance times s = 1, the workday length W equal to 15000 and finally, the overtime fee for exceeding workday length w is set to 10. These values were set to these levels after experimenting with various combinations. In Table 7 the cost of the LRPSD with workday recourse cost is shown for the Prins and Tuzun datasets. It can be seen that the ALNS performs quite well for smaller instances in Prins, which are sometimes even stopped before the maximum number of iterations is reached in a relatively small amount of time. However, when 200 customers are considered, the maximum time limit is often reached before 25000 iterations can be executed. It can also be seen that when the maximum vehicle inventory is 150, for all even numbered Prins instances, the cost are significantly smaller than when the maximum vehicle inventory is 70. However, the number of implemented iterations is smaller when the maximum vehicle inventory is 150. This is due to the fact that the vehicle capacities form an upper bound on the customer demands. When this upper bound is lower, fewer computations need to be made for the calculation of the expected route length, which results in a faster algorithm. Even when the objective approximation is used in the ALNS for the larger Prins instances, the time limit is reached before the maximum number of iterations. In Figures 1 and 2 the effect of greater maximum inventories on the route lengths is shown. In Figure 1 the maximum inventory is 70, while in Figure 2 it is 150. The depot locations are indicated with squares and the customer locations with circles. These figures show the effect of preventive and reactive replenishments and the workday exceeding recourse cost on the scheduled a priori route. When the maximum inventory is decreased, the probability of a route failure increases ceteris paribus, which has two separate effect. First, the cost of travelling to the depot for a replenishment is added. Second, this additional travelling time is included in the workday time.
For the Tuzun dataset, which contains relatively bigger instances, the time limit is reached for all in-stances. Here all maximum vehicle inventories are set to 150, which increases the computational com-plexity of the problem. In this dataset, all odd numbered instances contain 10 depots, whereas the even numbered instances contain 20. It can be seen that the even numbered instances are generally more time consuming, meaning fewer iterations were implemented before the time limit was reached. When more depots are included, there are more options that need to be considered in some of the operators, such as the depot swap removal operator.
5 Experiments 22 0 10 20 30 40 50 0 10 20 30 40 50 x coordinates y coordinates
Fig. 1: Route for the LRPSD for Prins7.
0 10 20 30 40 50 0 10 20 30 40 50 x coordinates y coordinates
Fig. 2: Route for the LRPSD for Prins8.
0 10 20 30 40 50 0 10 20 30 40 50 x coordinates y coordinates
Fig. 3: Route for Prins7 using the assumptions in Marinakis et al. (2016). 0 10 20 30 40 50 0 10 20 30 40 50 x coordinates y coordinates
5 Experiments 23
5 Experiments 24
5.4
Value stochastic solution
To show the importance of including the stochasticity of the LRPSD in the objective function, we calculate the value of the stochastic solution (VSS). The VSS indicates the benefit, or cost reduction, due to including stochastics in our model. To find the VSS we replace all stochastic variables by their means, meaning all customer demands are equal to the expected value with probability 1. Then we use the ALNS to solve this problem, for this solution the “actual” stochastic costs are then calculated. The difference between these costs and the costs found when the LRPSD including stochastic customer demands was solved is the VSS. In this situation preventive (and reactive) replenishments are still allowed, however, the route failures are known when the routes are scheduled. The resulting costs are presented in Table 10.
From Table 10 in Section 8.2 of the Appendix it can be seen that the VSS can take values up to 155080, which is 36% of the total costs for the LRPSD in that instance. On average the VSS for the instances considered is 39680.76. The VSS shows that it is imperative to include the stochastic elements of the LRPSD in the model even though this increases the computational complexity.
5.5
Maximum work day
This paper specifies a new recourse cost to incorporate the dilemma of either scheduling many vehicles or paying overtime due to routes exceeding the workday time into the objective function. This workday failure recourse cost combines the time necessary to visit customers with the time for routing, including replenishments.
In Huang (2015) a constraint on the probability of a route failure is used to incorporate the same dilemma. They assume Ph P j∈J ξj P u∈V xujk
5 Experiments 25 0 10 20 30 40 50 0 10 20 30 40 50 x coordinates y coordinates
Fig. 5: Route using exceeding workday recourse cost for Prins8. 0 10 20 30 40 50 0 10 20 30 40 50 x coordinates y coordinates
Fig. 6: Route using Huang (2015) constraint for Prins8.
Tab. 8: Huang constraint
Instance Huang Workday recourse
Distance Depots Routes Distance Depots Routes
6 Conclusion 26
6
Conclusion
We design an ALNS heuristic for the LRPSD with an additional recourse cost for routes that exceed the maximum workday time. The LRPSD is the problem where depot locations and a priori vehicle routes need to be determined simultaneously, while customer demands only become known upon arrival. We consider three different recourse costs: The cost of exceeding depot capacity, preventive and reactive replenishments at the depot for exceeding vehicle inventory and finally, the cost of exceeding the workday time. This final recourse cost incorporates the dilemma of short routes with many drivers or long routes where drivers need to work overtime into the LRPSD objective. In the workday we include the time spent at a customer as well as the time spent driving the route. To solve this problem we present a new ALNS which includes newly designed remove operators especially for the LRPSD. These new remove operators are based on the exceeding depot capacity, route failures and workday recourse costs.
Firstly, we validate our heuristic by calculating a lower bound on the objective value. We did this by relaxing the assumption that customer demands become known upon arrival. In this relaxation, the customer demands become known after the a priori routes are scheduled and before the vehicle departs from the depot. Using a Sample average approach we generated scenarios for customer demands and solved the resulting LP using Cplex for small instances. We found that the ALNS found solutions that are comparable to the lower bound.
Secondly, we use a benchmark from the literature to validate the ALNS. To compare our heuristic to the one presented in Marinakis et al. (2016), we remove the fixed vehicle costs and the costs of exceeding depot capacity and the workday failure cost. We find that our heuristic yields competitive results for the Tuzun dataset, finding 8 new best known results, and significantly lower results for the Prins dataset. This last difference might be due to an unreported assumption on customer demands exceeding the depot capacity.
We present the results of our ALNS for the LRPSD with exceeding workday recourse and find that for the Prins dataset the ALNS performs well. However, the ALNS is often terminated when the maximum time is exceeded for larger instances in the Tuzun dataset. In the Prins dataset the effect of a greater vehicle capacity is clear in both an increase in computation time and a reduction in cost due to a reduction in route failures. We compare solutions with the workday exceeding constraint to those for the Marinakis assumptions. It is clear that the solutions with the workday exceeding constraint are more sensitive to maximum vehicle capacities. This results in shorter routes and more opened depots than when the Marinakis assumptions are used.
We determine the value of the stochastic solution by comparing the stochastic cost of a solution to the cost of a solution found by assuming all customer demands equal their means. We find that including the stochastic demands in the model could lead to a reduction of up to 36%. Hence it is imperative that the stochastic demands are included in the model.
Finally, we determine the effect of the exceeding workday recourse cost. We compare the performance of our recourse cost to a constraint implemented in Huang (2015). This constraint also incorporates the problem of long routes into the objective value by ensuring that the probability of a route failure is at most 75%. We find that when the Huang constraint is used, the distance travelled and the number of routes scheduled are greater than when the workday recourse cost is used. When we compare the solutions of both approaches, it is clear that the constraint implemented by Huang emphasises short routes with few customers, whereas we consider both the distance travelled and the number of customers. By including both aspects are in our exceeding workday recourse, it is more flexible and more representative of the actual time spent during a workday.
7 Acknowledgements 27
each other’s customers. Currently, we assume all vehicles drive their own route, however, in practice if only one product is delivered vehicles could swap customers during their routes. This more dynamic situation has yet to be considered in research. Moreover, approximations of the costs of a route including the preventive and reactive routing costs should be researched. The exact cost calculations are very time consuming for large instances, which deteriorates the performance of the ALNS.
7
Acknowledgements
Computations were made on the supercomputer “briaree” from “Universit´e Laval”, managed by Calcul Qu´ebec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), the minist`ere de l’ ´Economie, de la science et de l’innovation du Qu´ebec (MESI) and the Fonds de recherche du Qu´ebec - Nature et technologies (FRQ-NT).
8
Appendix
8.1
Benchmark Prins dataset potential assumptions
When the Prins dataset is used for the LRP, fixed vehicle costs are added to the objective function. However, in the objective function presented in Marinakis et al. (2016), these are not present. It is possible that these were added in the calculation. From Table 9 we can see that even when the vehicle costs are included, the cost found by the ALNS is significantly smaller than the cost found by Marinakis et al..
Tab. 9: Potential assumptions by Marinakis et al. (2016)
Include vehicle cost Capacity assumption
Instance ALNS cost # Depots BKR ∆% ALNS cost # Depots BKR ∆%
8 Appendix 28
In Marinakis et al. (2016), depot failures are not considered, since they assume that the number of vehicles assigned to a depot multiplied by the maximum vehicle inventory is smaller than the depot capacity. However, since they do not consider a workday recourse cost, the maximum vehicle inventories can be exceeded. Therefore, the depot capacities can also be exceeded. We now consider a simplistic constraint that reduces depot failures. It states thatP
j∈Jλjzij ≤ Qi for all i ∈ I. This means that the
sum of mean demands for all customers assigned to depot i should not exceed the depot capacity. however, from Table 9, we can see that when this assumption is included, the resulting costs differ significantly from the cost found by Marinakis et al..
8.2
Value of Stochastic solution
In Table 10 the value of the stochastic solution is presented for the Prins dataset. Tab. 10: VSS for dataset Prins
8 Appendix 29
Tab. 11: Datasets Prins by Prins et al. (2004) and Tuzun by Tuzun and Burke (1999)
