A multi-start heuristic approach for the split-delivery vehicle routing problem with minimum delivery amounts

A multi-start heuristic approach for the split-delivery vehicle routing problem with minimum delivery amounts

Transportation Research Part E 88 (2016) 11–31 Contents lists available at ScienceDirect Transportation Research Part E journal homepage: www.elsevi...

1MB Sizes 10 Downloads 178 Views

Transportation Research Part E 88 (2016) 11–31

Contents lists available at ScienceDirect

Transportation Research Part E journal homepage: www.elsevier.com/locate/tre

A multi-start heuristic approach for the split-delivery vehicle routing problem with minimum delivery amounts Anthony Fu-Wha Han ⇑, Yu-Ching Chu 1 Department of Transportation and Logistics Management, National Chiao Tung University, 1001 Ta Hsueh Rd., Hsinchu 300, Taiwan

a r t i c l e

i n f o

Article history: Received 22 September 2015 Received in revised form 28 December 2015 Accepted 29 January 2016

Keywords: Vehicle routing Split delivery Minimal delivery amount Variable neighborhood descent Node-ejection chains Heuristics

a b s t r a c t We propose a new multi-start solution approach for the split-delivery vehicle routing problem with minimum delivery amounts (SDVRP-MDA). Initial solutions are generated by both node-insertion and route-addition procedures with a single parameter to control the restart. These solutions are then improved by a variable neighborhood descent metaheuristic with a novel search operator inspired by node-ejection chains. We test the proposed approach with 32 benchmark instances for four different minimum delivery fractions. Using the proposed algorithm, out of 128 cases tested, we find 81 best known solutions and 34 new best solutions; overall, we find 43 new best solutions. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction This paper addresses the split-delivery vehicle routing problem with minimum delivery amounts (SDVRP-MDA), in which the demand of a customer can be split and delivered to the customer by several vehicles on different routes, but the customers each require a minimum delivery amount (MDA) every time they are visited. The SDVRP-MDA, to which attention was first drawn by Gulczynski et al. (2010), is a restrained version of the split-delivery vehicle routing problem (SDVRP) introduced earlier by Dror and Trudeau (1989). The SDVRP and the SDVRP-MDA are both variants of the conventional vehicle routing problem (VRP). By allowing split deliveries without MDA requirements, the distance traveled and the number of vehicles employed can be potentially reduced by up to 50% (Archetti et al., 2006a, 2008a). Split deliveries may thus result in significant cost savings for in-house logistics operations. However, when shipped to outside customers, split deliveries incur extra work and create distractions for those customers who are visited more than once. These distractions are undesirable to customers, unless the amount or value of each delivery reaches a certain level of satisfaction for the customer. Gulczynski et al. (2010) thus proposed the SDVRP-MDA with an additional MDA constraint to the SDVRP. In practice, the MDA can be imposed as a monetary value by the distributor. For example, according to Gulczynski et al. (2010), Pizza Hut in the USA delivered orders placed by phone or online if the order was at least $8.00 (approximately). Similarly, the biggest chain store 7-Eleven in Taiwan (7-ELEVEN Taiwan, 2015) and Domino’s Pizza in Singapore (Domino’s Pizza Singapore, 2015) will do the same if the order is at least NT$300 (approximately $9.10) and SG$15 (approximately $10.73), respectively. Moreover, Everhealth Pharmaceuticals (Everhealth Pharmaceuticals, 2015) in Taiwan, and Tesco Groceries in ⇑ Corresponding author. Tel.: +886 35731680. 1

E-mail addresses: [email protected], [email protected] (A.F.-W. Han), [email protected] (Y.-C. Chu). Tel.: +886 35731682.

http://dx.doi.org/10.1016/j.tre.2016.01.014 1366-5545/Ó 2016 Elsevier Ltd. All rights reserved.

12

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

the UK (Tesco PLC, 2015) will deliver orders which are at least NT$2,000 (approximately $60.62) and £40 (approximately $61.30), respectively. Gulczynski et al. (2010) also proposed a minimum delivery fraction p (0 6 p 6 1) to determine the MDA for each customer i, i.e., MDAi = dp  di e where di is the demand of customer i. In this study, we set p in the domain of [0, 0.5], i.e., 0 6 p 6 0:5. This is because, when p > 0:5, it is impossible to have two deliveries each with a fraction of more than half, and the problem thus reduces to a conventional VRP. Accordingly, the VRP and the SDVRP can be regarded as special cases of the SDVRP-MDA with p > 0:5 and p = 0, respectively. Fig. 1 illustrates an example of a simple network with three customers, where the edge labels are distances, the node labels in parentheses are demands (or delivery amounts), and the vehicle capacity is 100. For p > 0:5, the problem is equivalent to a VRP. In its optimal solution (Fig. 1a), there are three routes without split deliveries, and the total distance traveled is 48. For p = 0, the problem is equivalent to a SDVRP; and its optimal solution (Fig. 1b) shows that customer 2 (denoted by two nodes in gray) receives two deliveries from two separate routes. In this SDVRP solution, the total distance traveled (i.e., 38), is less than that in the VRP solution (i.e., 48). In Fig. 1(c), the optimal solution for the SDVRP-MDA (with p = 0.2) contains two routes, and the total distance traveled (i.e., 43) is longer than that in the SDVRP solution (i.e., 38). Due to p = 0.2, each customer must have at least 20% of its demand delivered on a single route. As shown in Fig. 1(c), customer 3 receives two deliveries, and the smaller amount of these two is 30, which meets the MDA requirement (i.e., 30 > MDA3 = 19 = 95  0.2). On the other hand, the solution in Fig. 1(b) is infeasible for p = 0.2, because the delivery amount on one of the two routes for customer 2 is only 5, which is less than MDA2 (i.e., 7 = 35  0.2). Similar to the SDVRP that was proved to be NP-hard (Dror and Trudeau, 1990), the SDVRP-MDA is also NP-hard. Despite the abundant literature on the SDVRP, the existing literature on the SDVRP-MDA is limited. To the best of our knowledge, the only method available for solving the SDVRP-MDA in the literature is the hybrid heuristic method proposed by Gulczynski et al. (2010). In this hybrid method, conventional VRP solutions were generated first using a modified savings algorithm (Yellow, 1970). These initial VRP solutions were then improved by an endpoint mixed integer program with minimum delivery amounts (EMIP-MDA) to provide the first solutions for the SDVRP-MDA. These solutions were further improved by an enhanced record-to-record travel algorithm (ERTR) developed by Groër et al. (2011). Overall, the method was called EMIP-MDA + ERTR. Gulczynski et al. (2010) also proposed a new set of test problems for the SDVRP-MDA by modifying the 21-instance set introduced by Chen et al. (2007) for the SDVRP. These instances were applied to test EMIP-MDA + ERTR for four different p values (i.e., p = 0.1, 0.2, 0.3, and 0.4) with good results. Recently, Xiong et al. (2013) proved that, in general, the maximal saving in travel cost for the SDVRP-MDA is the same as that for the SDVRP (i.e., as much as 50%). The only exception is for a special case when p = 0.5. In such a case, the potential saving is as much as one third. In this paper, we propose a new multi-start two-phased variable neighborhood descent heuristic approach for solving the SDVRP-MDA. This multi-start approach is implemented with a parameter a (0 6 a 6 1). For each a in a restart, initial solutions are first generated by a split-delivery route construction algorithm, which we call SRC. These initial solutions are then improved by a variable neighborhood descent (VND) procedure. Overall, the proposed approach is referred to as SRC + VND, which incorporates the following unique features for the SDVRP-MDA: (1) two different yet complementary constructive procedures to enhance the diversification of the search, (2) an adaptive procedure to construct MDA-restrained splitdelivery routes using a mix of node insertion and route addition steps, and (3) a novel neighborhood search operator inspired by node-ejection chains to enhance the intensification of the search. We have tested our solution algorithms with all available benchmark problems of the SDVRP-MDA. Out of the 32 instance problems each tested with four different p values (i.e., p = 0.1, 0.2, 0.3, and 0.4), we have obtained 81 existing best solutions and 34 new ones from the proposed SRC + VND algorithm. Moreover, the proposed algorithm relies on only one parameter and requires minimal tuning efforts.

Fig. 1. Comparison of the VRP, SDVRP and SDVRP-MDA. Nodes in gray denote split deliveries.

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

13

The remainder of the paper is organized as follows. In Section 2, we review the literature on the SDVRP. In Section 3, we present a mathematical formulation for the SDVRP-MDA with a proof of the integrality property for its optimal solution. The proposed SRC + VND solution algorithm with a novel neighborhood search operator are described in Section 4. The computational results are presented in Section 5. Finally, Section 6 concludes the paper.

2. Literature review on the SDVRP The earliest heuristics for solving the SDVRP were proposed by Dror and Trudeau (1989, 1990) which incorporated k-split interchanges and route additions. Since then, many different approaches have been developed to solve the SDVRP. Among these approaches, tabu search seems to be mostly adopted. Ho and Haugland (2004) presented the first study which applied tabu search to solve the problem with time windows. Archetti et al. (2006b) implemented a three-phased tabu search algorithm called SPLITABU and improved most of the instances considered by Dror and Trudeau (1989). Aleman and Hill (2010) presented a tabu search with a vocabulary building approach called TSVBA. This TSVBA is a population-based search approach that uses a set of solutions to determine attractive attributes with which to construct new solutions. Recently, Berbotto et al. (2014) proposed a tabu search heuristic called randomized granular tabu search (RGTS), in which a granular neighborhood is defined by the most promising arcs in a quality solution. Both the TSVBA and the RGTS demonstrated their competitiveness with other approaches in the SDVRP literature. Metaheuristics other than tabu search have also been applied to solve the SDVRP. Using evolutionary methods, Mota et al. (2007) and Campos et al. (2008) took a scatter search approach to solve the SDVRP. Boudia et al. (2007) presented an effective memetic algorithm with population management techniques called MAjPM and improved most of the benchmark instances they tested. Derigs et al. (2010) applied five metaheuristic methods, namely simulated annealing, threshold accepting, record-to-record travel, attribute based hill climber (ABHC), and attribute based local beam search for solving the SDVRP, and found that the ABHC method showed the best performance. Aleman et al. (2010) proposed a threephased heuristic called ICA + VND. Both the ABHC and the ICA + VND heuristics improved the best known solution (BKS) for several benchmark instances in their tests. Hybrid heuristics, which integrate metaheuristics with integer programs to facilitate the search, have also contributed to solving the SDVRP. Chen et al. (2007) developed a hybrid approach called EMIP + VRTR, which combined an endpoint mixed integer program (EMIP) and the variable-length neighbor list record-to-record travel algorithm (VRTR) proposed earlier by Li et al. (2005). Archetti et al. (2008b) presented an optimization-based heuristic using an integer program to explore promising parts of the search space identified by a tabu search heuristic. Note that the EMIP + VRTR approach actually provides the fundamental framework for the later work on EMIP-MDA + ERTR by Gulczynski et al. (2010) in solving the SDVRP-MDA. Some exact-solution approaches were developed for determining and improving the bounds of the optimal solution for the SDVRP. Belenguer et al. (2000) performed a polyhedral analysis, explored several valid inequalities, and implemented a cutting-plane algorithm. Jin et al. (2007) presented a two-stage exact algorithm, while Jin et al. (2008) improved the bounds generated by Belenguer et al. (2000) with a column generation algorithm. Recently, Archetti et al. (2014) presented two branch-and-cut algorithms that provide lower bounds to the optimum. These algorithms were able to solve most of the instances with up to 50 customers and two instances with 75 and 100 customers. lP m N In a SDVRP, the minimal possible number of vehicles, K min , can be determined as K min ¼ i¼1 di =Q . When split deliveries are not restricted, it can be easily verified that a feasible solution always exists for the SDVRP with K min vehicles. Accordingly, Archetti et al. (2011) defined the problem with two versions: SDVRP-LF and SDVRP-UF. The former is for a fleet size limited to K min , and the latter is for an unlimited fleet size. Unless specified otherwise, the SDVRP in the literature usually refers to the unlimited case (i.e., SDVRP-UF). More recently, Silva et al. (2015) presented a multi-start iterative local search based heuristic with an innovative perturbation mechanism called SplitILS to solve the SDVRP-LF and SDVRP-UF. For both cases, it was found that the SplitILS outperformed almost all of the existing methods reported in the literature and improved the BKS for 243 out of 324 instances tested. Note that in the SDVRP-MDA, only the case with an unlimited fleet is considered since the limited-fleet case is not valid. More about the modeling, solution approaches, variants, and applications of the SDVRP can be found in a complete survey by Archetti and Speranza (2012).

3. Mathematical formulation and integrality property of the SDVRP-MDA The SDVRP-MDA can be defined on a complete graph G = (V; A), where V ¼ f0g [ N is the set of nodes and A is the set of arcs, i.e., A ¼ fði; jÞji; j 2 V; i–jg. The node 0 represents the depot and N = {1, 2, . . . , n} represents the set of n customers. The distance of an arc ði; jÞ 2 A is cij . As mentioned earlier, p denotes the minimum delivery fraction (0 6 p 6 0:5), and the minimal delivery amount for a customer i 2 N is MDAi=dp  di e. Assume that there is a homogeneous fleet of m vehicles each with a capacity of Q. Every vehicle must start and end its route at the depot. For each customer i, the demand must be satisfied, and the amount delivered each time must exceed MDAi. Let xvij be equal to 1 if vehicle v travels directly from i to j, and

14

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

equal to 0 otherwise, and qiv be the amount of the demand of i delivered by vehicle v. The problem is to find a set of routes (in terms of xvij ) and related delivery amounts qiv that gives a minimum total travel distance. Based on the formulation of the SDVRP proposed by Archetti et al. (2006b), we present a mixed integer programing formulation (P) for the SDVRP-MDA as follows:

min

ðPÞ

n X n X m X cij xvij i¼0 j¼0

ð1Þ

v ¼1

subject to: n X m X xvij P 1; i¼0

v ¼1

for j ¼ 0; . . . ; n;

n n X X xvik  xvkj ¼ 0; i¼0

ð2Þ

for k ¼ 0; . . . ; n;

v ¼ 1; . . . ; m;

ð3Þ

j¼0

XX xvij 6 jSj  1;

for

v ¼ 1; . . . ; m;

S # N;

ð4Þ

i2S j2S

qiv 6 di

n X

xvij ;

v ¼ 1; . . . ; m;

for i ¼ 1; . . . ; n;

ð5Þ

j¼0 m X qiv ¼ di ;

v ¼1

n X qiv 6 Q;

for i ¼ 1; . . . ; n;

for

ð6Þ

v ¼ 1; . . . ; m;

ð7Þ

i¼1

qiv P MDAi ; xvij 2 f1; 0g;

v ¼ 1; . . . ; m;

for i ¼ 1; . . . ; n; for i ¼ 0; . . . ; n;

j ¼ 0; . . . ; n;

v ¼ 1; . . . ; m:

ð8Þ ð9Þ

Constraints (2)–(4) are the standard constraints for the VRP; constraints (2) ensure that each customer is visited at least once; (3) are the flow conservation constraints, and (4) are the subtours elimination constraints. Constraints (5)–(7) specify how the demands of the customers are allocated among the vehicles. More precisely, constraints (5) impose the restriction that customer i can be served by vehicle v only if v passes through i. Constraints (6) ensure that the entire demand of each customer is satisfied, while constraints (7) require that the amount delivered by each vehicle does not exceed the capacity. Constraints (8) define the MDA restrictions for the amounts delivered to each customer by a single vehicle. We also present the integrality property of the optimal solution for the SDVRP-MDA as follows: Theorem 1. If P is feasible, and Q ; di and MDAi are all integers, then there always exists an optimal integer solution to P. Proof. Given an optimal solution of P, let xvij  and qiv be the corresponding optimal solutions. Let B denote the set of pairs P (i; v ) with i 2 N, v 2 f1; . . . ; mg, such that nj¼0 xvij  P 1. For each customer i, define SðiÞ ¼ fv jði; v Þ 2 Bg, i.e., the set of vehicles

that pass through i. Similarly, define for each vehicle v the set Dðv Þ ¼ fijði; v Þ 2 Bg, i.e., the set of customers to whom deliveries are made by v. Given the solution values of xvij  , constraints (2), (3), (4) and (9) are satisfied, and constraints (5) can be

covered by constraints (6) and (8) with the (i; v ) relation embedded in B. The remaining constraints (6)–(8), which are related only to qiv , can be written as follows:

X v 2SðiÞ

qiv ¼ di ;

for i ¼ 1; . . . ; n;

qiv 6 Q;

for

X i2Dðv Þ

qiv P MDAi ;

v ¼ 1; . . . ; m;

for ði; v Þ 2 B:

0

ð6 Þ

0

ð7 Þ

0

ð8 Þ

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

15

We set yv i ¼ qiv and use it to rewrite the above three constraints as follows:

X

v 2SðiÞ

yv i ¼ di ;

for i ¼ 1; . . . ; n;

yv i 6 Q;

for

X

i2Dðv Þ

yv i P MDAi ;

ð10Þ

v ¼ 1; . . . ; m;

ð11Þ

for ði; v Þ 2 B:

ð12Þ

Constraints (10)–(12) now represent a set of constraints for a minimum cost flow problem, which has a bipartite network structure with supply-only nodes, i.e., the v nodes in (11), demand-only nodes, i.e., the i nodes in (10), and the flow variables P yv i in (12). This minimum cost flow problem is feasible, i.e., mQ P ni¼1 di , because it is based on an optimal solution of P. Therefore, when the supply (i.e., Q ), the demand (i.e., di ) and the lower bound of the flow variables (i.e., MDAi) are all integers, the integrality property of network flow solutions indicates that it always has an integer optimal flow solution (Ahuja et al., 1993). In other words, for a given optimal solution of xvij  for P, we can use constraints (10)–(12) to define a corresponding minimum-cost network subproblem for assigning the flows of yv i (i.e., qiv ). In addition, for this subproblem, there always exists a set of integer solutions qiv for P as long as Q ; di and MDAi are all integers. h 4. Solution algorithm of SRC + VND Our proposed solution algorithm, which is called SRC + VND, can be summarized by the pseudocode presented in Table 1. This algorithm has a single parameter a (0 6 a 6 1), which is controlled by Da for the multi-start implementation of SRC + VND. Specifically, the algorithm first starts with a = 0 (line 3), and then restarts itself each time with a different a, i.e., a þ Da (line 16) until a > 1. This a has an important role in our route construction procedures, and will be defined in Section 4.1. For a specific a, the solution process has a two-phased structure. In the first phase, the SRC constructive algorithm is used to generate initial solutions. These initial solutions are then improved by a VND procedure in the second phase. We include two procedures of SRC (i.e., SRC0 and SRC1) to construct diversified initial solutions. Accordingly, two initial solutions (s0 and s1 ) and two improved solutions (s00 and s01 ) are generated (lines 5–8). The better one of s00 and s01 defines the final solution s obtained for the given a in this run (lines 9–12), and the incumbent solution is updated accordingly (lines 13–15). The algorithm then proceeds to restart with a different a (line 16), and repeats itself until a > 1. In the end, the best solution obtained from all a tested provides the final solution sbest . For simplicity, the terms customer and node will be used interchangeably hereafter in this paper, as will the terms distance and cost. A customer whose demand is not fully served during the solution process will be referred to as an unfulfilled customer. 4.1. Multi-start initial solution construction with SRC As mentioned earlier, our proposed SRC includes two constructive procedures, namely, SRC0 and SRC1. Both of these procedures have a common structure and can be described by the pseudocode presented in Table 2. In line 1 of this table, a declarative parameter type is used to define a specific constructive procedure, i.e., type = 0 for SRC0, and type = 1 for SRC1. The entire SRC constructive algorithm starts with the initialization of the set of routes R, the list of unfulfilled customers L, the unsatisfied demand ui , and the set of forbidden routes F i for an unfulfilled customer i (lines 2–4). The initial route Table 1 The pseudocode of SRC + VND.

16

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31 Table 2 The pseudocode of SRC.

r1 is generated by connecting the depot with the farthest node f from the depot, i.e., r 1 ¼ ð0  f  0Þ, and removing f from the unfulfilled list L (lines 5–7). This initial route r 1 will be treated differently for SRC0 and SRC1. For SRC0, r1 is added directly to R and becomes available for further route construction work (lines 8–9). On the other hand, r1 is kept intact until the end of the algorithm when all customers are fully served in SRC1 (lines 23–24). The intention of keeping r1 intact in SRC1 is twofold: (1) to improve diversification with a significant different pattern of initial solutions generated from SRC0, and (2) to enhance intensification by providing more space for improvement in the VND than SRC0. After the construction of r 1 , the SRC continues to a node selection process (lines 10–20). The key to this process is the weighted combination cost (line 18) defined as follows:

W i ¼ aIi  ð1  aÞc0i ;

for i 2 L

ð13Þ

where a represents the weighting parameter (0 6 a 6 1) of W i , and Ii denotes the cheapest insertion cost of node i, i.e., Ii ¼ minr2R Iir , where Iir is the cost of inserting node i in the best position of an existing route r. This weighted combination cost is computed for all unfulfilled nodes. The node that yields the lowest cost is selected to proceed with an adaptive procedure (line 22) to construct SDVRP-MDA routes using a mix of node insertion and route addition steps. This adaptive route construction procedure is referred to as NIRA (for Node Insertion and Route Addition), and will be introduced shortly. The entire SRC constructive algorithm repeats itself until all customers are fully served, i.e., L = / (and the first route r 1 is settled for SRC1). Note that during the node selection process, some routes, i.e., those forbidden routes in F i (see line 14), must be excluded from R to prevent the constructive algorithm from diverging. This will be explained in the following description of NIRA. The NIRA procedure is an adaptive route-construction subroutine for an unfulfilled node i which attempts insertion into a route r. Table 3 presents the pseudocode of NIRA. When the (cheapest) insertion cost Iir is smaller than 2c0i , a node-insertion process is implemented (lines 3–18). Otherwise, a new route ð0  i  0Þ is created and added to the set of existing routes R (lines 19–22). For the node-insertion process, if the vehicle on route r has sufficient spare capacity to accommodate ui , i.e., the unfulfilled demand of node i, this process is completed in a straightforward manner (lines 3–4). However, if this is not the case, an overloading occurs, and we must remove at least Or , i.e., the overloaded amount of the vehicle on route r, from a certain node v on route r. This node v is determined as follows (line 7):

v ¼ arg

minfc0p jqpr P Or g; p2r

ð14Þ

Specifically, among those p nodes whose individual removal can release the overloading of r, i.e., p 2 r and qpr P Or , the one nearest to the depot is chosen as the removing node v. The nearest-to-depot node is selected for removal because when it is processed later in the succeeding iterations of NIRA, it tends to incur a lower cost than the other nodes in r due to its proximity to the depot. Once the removing node v is selected, we need to decide how much to remove from qv r in order to release route r from overloading, and to ensure that subsequent deliveries to v meet the MDA requirements. The amount qv r can be removed

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

17

Table 3 The pseudocode of NIRA.

either completely or partly from route r. The latter (the partial removal case) is preferred to the former (the complete removal case) because the more of qv r that remains in r, the more efficiently the routes are constructed during this process. Furthermore, for the partial removal case, qv r can be considered to consist of two parts: one is the amount that remains in r, which must exceed MDAv , and the other one is the amount that is to be removed from route r. The latter must satisfy the minimal shift amount (MSAv) defined as follows:

MSAv ¼ maxðOr ; MDAv Þ;

ð15Þ

where Or is the least amount required to release r from overloading, and MDAv is required for later deliveries to the unfulfilled node v. Only when qv r is sufficient to cover both the MDAv and MSAv, i.e., (qv r  MSAv) P MDAv, the MSAv part of qv r will be removed from r (lines 10–12). Otherwise, qv r is completely removed from r (lines 14–17). At the end of both cases, node v will be returned to the unfulfilled list L for later processing. At the same time, r must be marked as a forbidden route for v (line 18). Otherwise, the same (v ; r) pair tends to be selected again for a later step as the input to the NIRA procedure for v. Once this happens, it is very likely that v will be rejected from r once again. This will result in cyclic iterations which cannot converge. An obvious case, which will directly lead to such a divergence, occurs when the node v selected (line 7) is the same as the input node i (line 1) in a NIRA execution.

4.2. Solution improvement using VND with a new operator of SNEC Our solution algorithm employs a VND approach for solution improvement. The VND is the deterministic version of the variable neighborhood search (VNS) metaheuristic proposed by Mladenovic´ and Hansen (1997). The basic idea of VND (or VNS, in general) is to successively explore a set of predefined neighborhood structures to provide a better solution. Given Table 4 The pseudocode of VND.

18

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

N k (k = 1, . . . , kmax) a finite set of pre-selected neighborhood structures, N k ðxÞ the set of solutions in the kth neighborhood of x, and s an initial solution, our VND can be described by the pseudocode in Table 4. For each of the neighborhood structures, the algorithm keeps searching for an improved solution s0 until no further improvement can be made, and then proceeds to search with the next neighborhood structure (lines 4–9). The process repeats itself until an improved solution cannot be found with all the kmax neighborhood structures (lines 2–10). In the end, the final solution is a local optimum with respect to all neighborhoods N k , and thus the chances of finding a global optimum are higher than by using a single neighborhood structure. We employ five neighborhood structures in our VND, which include four conventional VRP local search operators and an innovative operator inspired by node-ejection chains. The first two neighborhood structures (N 1 and N 2 ) are for intra-route improvement:  2-opt (N 1 ): This arc-exchange operator was introduced by Lin and Kernighan (1973), in which two nonadjacent arcs are deleted and another two are added to generate an improved new route.  Or-opt (N 2 ): This node-shift operator was introduced by Or (1976), in which one, two or three consecutive nodes are removed and inserted in another position of the route for improvement. The next two neighborhood structures (N 3 and N 4 ) are for inter-route improvement:  2-opt⁄ (N 3 ): This arc-exchange operator was first introduced by Potvin and Rousseau (1995) for solving the VRP with time windows, and thus it does not accept inverted segments in the modified routes. We have relaxed this restriction because our problem does not have time window constraints. Fig. 2 illustrates the 2-opt⁄ implemented in our VND, where two arcs (i, i + 1) and (j, j + 1) are removed (Fig. 2a), and reconnected to generate two new routes: one accepts non-inverted segments (Fig. 2b) and the other one accepts inverted segments (Fig. 2c). The squares and circles in Fig. 2 represent the depots and customer nodes, respectively.  k-interchanges (N 4 ): The original node-exchange k-interchanges scheme (Osman, 1993) considers the exchanging of up to k customers between two routes. To limit the number of possibilities, we have considered k = 2. According to Cordeau and Laporte (2005), these exchanges are better described as swap operations of (k1, k2) (with k2 6 k1 6 k), in that each represents an operation where k1 customers are shifted from route 1 to route 2, and k2 customers are removed from route 2 to route 1. In our implementation, we have only adopted the swap operations (1, 2) and (2, 2). This is because the swap operations (1, 0) and (1, 1) can be covered in the fifth neighborhood structure as described in what follows, and the swap (2, 0) might be achieved by two separate (1, 0) operations. Some additional work is necessary for applying N 1 through N 4 to the SDVRP-MDA. First, the MDA requirement must be met for all moves involved in these neighborhood operations. Moreover, when a better neighborhood solution is found in N 3 (i.e., 2-opt⁄) or N 4 (i.e., k-interchanges), the modified routes must be checked to see if they contain customers visited twice with two deliveries on a single route. If this is the case, a repair step must be called to combine those two deliveries on the route. The fifth neighborhood structure (N 5 ) is formed by a completely new inter-route node exchange operator, i.e., the splitdelivery node ejection chains (SNEC), which is inspired by node-ejection chain ideas. The first ejection chain procedure was proposed by Glover (1992) for the traveling salesman problem. Rego (2001) adopted the ejection chain idea and modified

i

j+1

j

i+1

r1 i

j+1

r2

r1

(b)

r2

i j

i+1

j+1

r2

r1

(a) j

Fig. 2. Illustration of the 2-opt⁄ implemented.

(c)

i+1

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

19

it to the local search operator node-ejection chains for the VRP. These node-ejection chains, in general, are a series of segments (called triplets by the author), each consisting of three consecutive nodes in a route. In its operation, all the middle nodes in the triplets are selected to proceed with a multi-node exchange or insertion process to find a better neighborhood solution. Unlike the node-ejection chains operator, which proceeds on a pre-determined structure of triplets, the new SNEC operator is a highly adaptive process. Specifically, once a node i is ejected from route r 1 , how the nodes and routes are selected for subsequent processing, and how much the delivery amounts are to be transferred between and among routes depend on the chain reactions triggered by the ejected node i. Fig. 3 illustrates all the possible moves involved during the adaptive process of the SNEC operator. The operator begins with ejecting a node i from route r1 to a different route r2 . This includes two possible initial states in which node i is either served by only r1 , as shown in Fig. 3(a1), or served by r1 and r2 with two deliveries, as shown in Fig. 3(a2). The operator then inserts node i into the cheapest location of route r 2 . If r 2 is not overloaded, then node i is settled with r2 as depicted in Fig. 3(b). However, if r 2 is overloaded, the operator proceeds to transfer some or the entire delivery amount from a different node j (j–i) on r2 to another route r3 (r 3 –r 2 ) in order to maintain the loading feasibility of r2 . Furthermore, if the vehicle on r 3 has sufficient spare capacity, then node j (with its delivery amount qjr2 ) is completely shifted from r 2 to r 3 , as shown in Fig. 3(c1). Otherwise, qjr2 can only be partly transferred to r3 , and node j will be served by two routes, i.e., r 2 and r 3 , as shown in Fig. 3(c2). During this process, SNEC checks for all possible i, r 1 , r2 and, if necessary, all j and r 3 to find a best neighborhood solution. More details of SNEC are described with a pseudocode in Appendix A.

Fig. 3. Illustration of the SNEC operator.

Table 5 Special-case scenarios of SNEC and the corresponding operators. Special-case scenario

1 2 3 4 5 a b c d e

Corresponding operators

Start

End

Remark

(a1) (a1) (a1) (a2) (a2)

(b) (c1) (c2) (c2) (c2)

– r3 r3 r3 r3

Ho and Haugland (2004). Aleman et al. (2010). Derigs et al. (2010). Berbotto et al. (2014). Silva et al. (2015).

¼ r1 ¼ r1 ¼ r1 – r1

relocatea, shiftb, Relocated, Shift(1, 0)e exchangea, swapb, EXCHANGEc, Exchanged, Swap (1, 1)e shift⁄b, Swap (1, 1)⁄e relocate splita, ExchSplitd DelSplitNewd

20

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

Note that in SNEC, the first route must always be different from the second route (i.e., r 1 –r 2 ), and the second and third routes must also be different (i.e., r 2 –r 3 ). Nevertheless, the third route r3 , which appears to be different from r1 as shown in Fig. 3(c1) and (c2), can be the same as r1 . The design of adaptive moves in SNEC makes it a rather comprehensive and general operator compared to many existing operators that have been applied for solving the SDVRP. In fact, these operators can be regarded as special cases of SNEC. Such operators include relocate, exchange, and relocate split by Ho and Haugland (2004); shift, swap, and shift⁄ by Aleman et al. (2010); EXCHANGE by Derigs et al. (2010); Relocate, Exchange, ExchSplit, and DelSplitNew by Berbotto et al. (2014); and Shift (1, 0), Swap(1, 1), and Swap (1, 1)* by Silva et al. (2015). Table 5 summarizes five special-case scenarios of SNEC, each of which corresponds to some existing operators. For example, the fourth special-case scenario, in which SNEC strictly operates from Fig. 3(a2) and terminates at Fig. 3(c2) with r 3 ¼ r1 , corresponds to the search operations of both relocate split (Ho and Haugland, 2004) and ExchSplit (Berbotto et al., 2014). Accordingly, SNEC tends to explore a larger neighborhood than those corresponding local search operators listed in Table 5. We believe that SNEC plays an important role in the intensification of the search for our VND.

5. Results and analysis Benchmark instances are limited for the SDVRP-MDA. At this time, only two sets of problem instances are available for the SDVRP-MDA. The first set (Set G) was proposed by Gulczynski et al. (2010), which is based on the 21-instance problem set (Set C) introduced earlier by Chen et al. (2007) for the SDVRP. For Set G (and Set C), all instances were designed with symmetrically distributed customers. Due to this symmetric design, Gulczynski et al. (2010) obtained near-to-optimal solutions composed of some pre-designed route patterns by inspecting the instance problems visually. In order to keep the routes in these visually estimated solutions unchanged, the authors devised a systematic way to modify the original demands (in Set C) for the corresponding instances in Set G according to different minimum delivery fractions. Except for these modified demands, all other attributes of the instances in Set G, which include the number and location (i.e., node coordinates) of customers, and the vehicle capacity (i.e., 100), were kept the same as those in Set C. For each instance, the number of customers, which ranges from 8 to 288, is also reflected in the name of the instance. For example, instance MD10_64 represents the 10th problem in Set G which has 64 customers. The second set (Set B) of problem instances was first introduced by Belenguer et al. (2000) for the SDVRP. This problem set was adopted by Gulczynski et al. (2010) for the test of the SDVRP-MDA. In the instances of this set, the coordinates of nodes were set to correspond to those of the instances eil51, eil76, and eil101 from TSPLIB, and the vehicle capacity Q was set to 160. The demands were generated from a uniform distribution in the interval [aQ, bQ] with six scenarios defined by the following range bounds (a, b): (0.01, 0.1), (0.1, 0.3), (0.1, 0.5), (0.1, 0.9), (0.3, 0.7), and (0.7, 0.9). These six demand scenarios were denoted by D1 thorough D6, respectively, and these notations were included as a part of the name for each instance. For example, instance S76D3 represents a problem with 75 customers and a demand range bound of D3, i.e., (0.1, 0.5). Note that there were originally 14 instances in Set B, but only 11 were tested by Gulczynski et al. (2010) for the SDVRP-MDA. The full data of the instances in Set B were made available by Belenguer et al. (2000) from http://www.uv.es/belengue/ SDVRP/sdvrplib.tar.gz. As for the instances of Set G, Gulczynski et al. (2010) provided a computer program (generator) to generate the associated attribute data of the instances. However, for different computer executions, this generator may yield real number values with different decimal places in the node coordinates for an instance. This tends to yield different cost data for the same instance due to rounding errors. In order to avoid this problem, we have adopted the specific data of those coordinates with two decimal places from Gulczynski (2010). We have also made this data of Set G available from https://db. tt/ifRgKVuc, for future benchmark tests for the SDVRP-MDA. We have coded the proposed solution algorithms in C# and have implemented them on a personal computer with a 3.4 GHz i7-3770 processor, 16 GB of RAM and the Windows 7 operating system. For parameter tuning and VND setting, we first set Da ¼ 0:01 and test all 120 permutations of the five neighborhood structures, i.e., N 1 , N 2 , N 3 , N 4 and N 5 , with all instance problems in Set G and Set B. The best permutation found is N 1 , N 4 , N 2 , N 3 and N 5 (i.e., 2-opt, k-interchanges, Or-opt, 2-opt⁄ and SNEC). We then use this VND permutation to test for a better setting of Da. The results are shown in Fig. 4. In general, the smaller Da is, the better our algorithm performs yet the more computation time that is required. Considering both accuracy and speed, we choose Da ¼ 0:008, for which the average gap (i.e., the average percentage deviation from the BKS) is 0.20% and the average computation time is 37.2 s. With Da ¼ 0:008, we test again for all VND permutations, and find that the best VND permutation is still N 1 , N 4 , N 2 , N 3 and N 5 . We thus adopt these settings for SRC + VND, and report our test results in the following sections.

5.1. Test results for the case of SDVRP (p = 0) The SDVRP is a special case of the SDVRP-MDA for p = 0. Therefore, a good algorithm for the SDVRP-MDA should also work for the SDVRP. As reviewed in Section 2, an abundance of research has been devoted to the SDVRP during the last two and a half decades. Moreover, quite a few outstanding studies with breakthrough results were published most recently, e.g., Berbotto et al. (2014), Archetti et al. (2014), and Silva et al. (2015). All these works, which include heuristic and

21

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

exact-solution approaches as well as the limited fleet (LF) and unlimited fleet (UF) cases, are compiled to update the BKS for our SDVRP test. Tables 6 and 7 present the results of our SDVRP test for Set C and Set B, respectively. The proposed SRC + VND algorithm is compared with EMIP-MDA + ERTR and four other very competitive algorithms for the SDVRP, namely, TSVBA (Aleman and Hill, 2010), RGTS (Berbotto et al., 2014), SplitILS-LF, and SplitILS-UF (Silva et al., 2015). The SRC + VND is found to compare favorably with most of these best available algorithms with an average gap of 0.13% for Set C (21 instances), and such a gap of 1.07% for Set B (11 instances). Although it is second to SplitILS-LF and SplitILS-UF, the proposed SRC + VND outperforms EMIP-MDA, RGTS and TSVBA. When compared with EMIP-MDA + ERTR alone, SRC + VND outperforms 22 instances (i.e., 13 for Set C, and 9 for Set B) and equals 7 for all the 32 instances tested. Furthermore, for each of the eight instances with more than 100 customers, SRC + VND provides a better solution than EMIP-MDA + ERTR. As to the computation speed, SRC + VND executes much faster than EMIP-MDA + ERTR. For Set C, the average computation time is 75.7 s for SRC + VND, while that for EMIP-MDA + ERTR is 1004.3 s. Similarly, for Set B, the average computation times for SRC + VND and EMIP-MDA + ERTR are 21.3 and 487.9 s, respectively. In Table 6, we find slightly different objective (total cost) values in the best solutions for three instances in Set C due to rounding errors in their cost data. These three instances are SD12_80, SD16_144, and SD17_160. 5.2. The SDVRP-MDA test results for Set G Table 8 presents the results of our SDVRP-MDA test for Set G with p = 0.1, 0.2, and 0.3. The BKS of this test are the visually estimated solutions from Gulczynski et al. (2010). In the table, initial solution results obtained from SRC are also included for comparison due to their near-to-BKS performance. The results show that SRC and SRC + VND both perform better than EMIPMDA + ERTR. For each group of 21 instances tested with a different p value, we find 11 BKS from SRC, and all 21 BKS from SRC + VND; and in each group, the average gap for SRC (and SRC + VND) is no more than 0.32%. On the other hand, for EMIPMDA + ERTR in each of the three test groups, no more than 8 BKS are obtained and the average gap is larger than 0.93%. In Table 8, we observe that seven instances yield slightly different objective values from the same best solutions due to rounding errors in their cost data. These seven instances are MD10_64, MD15_144, MD17_160, MD18_160, MD19_192, MD20_240, and MD21_288. The results of our SDVRP-MDA test for Set G with p = 0.4 are presented in Table 9. We put these results in a separate table because the benchmark solutions for this test group (p = 0.4) are not entirely the same as those (visually estimated solutions) for the other three groups (p = 0.1, 0.2 and 0.3). Specifically, when tested with p = 0.4, EMIP-MDA + ERTR resulted in better solutions than those visually estimated solutions for two instances MD6_32 and MD16_144 (Gulczynski et al., 2010). The BKS of these two instances are thus different for Tables 8 and 9. Our constructive algorithm SRC no longer performs better than EMIP-MDA + ERTR in the test group for p = 0.4. As shown in Table 9, SRC performs with an average gap of 3.33%, which is larger than the gap of 1.00% for EMIP-MDA + ERTR. On the other hand, the proposed SRC + VND algorithm outperforms EMIP-MDA + ERTR significantly with an average gap of 0.06%. For the 21 instances tested with p = 0.4, we find 18 BKS from SRC + VND. In addition, we obtain from SRC + VND a new best solution for MD21_288, which has the biggest problem size (with 288 customers) in our test. This improved BKS is underlined in Table 9. We also observe that, for p = 0.4, six instances yield slightly different objective values from the same best solutions due to rounding errors in their cost data. As shown in Table 9, these six instances are MD14_120, MD15_144, MD16_144, MD17_160, MD18_160, and MD20_240.

50

0.00 Time time(s) (seconds)

gap (%) (%) Dev.

42.13s

40 -0.08%

32.06s

30

-0.10

26.62s

gap (%)

time (seconds)

-0.05

37.23s

20.62s

20

-0.15 -0.15% 13.48s

-0.17%

-0.19% -0.18%

10

-0.20

-0.20%

0

-0.25 0.0400 (1/25)

0.0200 (1/50)

0.0133 (1/75)

0.0100 (1/100)

0.0080 (1/125)

Δα Fig. 4. Results of the parameter test of Da.

0.0067 (1/150)

22

Table 6 Results of the SDVRP test for Set C. Set G instances

Average Avg. gap (%)

EMIP-MDA + ERTR

TSVBAa

RGTSb

SplitILS-LFc

SplitILS-UFc

SRC + VND

Cost

Timed (s)

Cost

Timee (s)

Costf

Timef,g (s)

Costh

Timeh,i (s)

Costh

Timeh,i (s)

Cost

Timej (s)

228.28k 708.28k 430.58k 631.05k 1390.57k 831.24k 3640.00k 5068.28k 2044.20k 2684.88k 13280.00k 7213.61m 10110.58 10715.53m 15093.85m 3381.26m 26493.56 14197.80m 19989.95m 39635.51m 11271.06n

228.28 708.28 430.58 631.05 1390.57 831.24 3640.00 5092.36 2044.20 2704.69 13358.30 7256.77 10141.80 10780.00 15216.30 3382.16 26640.70 14357.80 20348.20 39902.80 11436.70

3.10 424.10 20.80 454.70 655.50 656.00 669.90 678.30 662.30 700.90 709.00 729.00 729.20 1718.10 1664.60 1654.80 1785.30 1723.00 1787.60 1839.80 1825.10

228.28 708.28 430.58 631.05 1390.57 831.24 3640.00 5068.28 2071.03 2747.83 13280.00 7213.62 10110.58 10802.87 15153.45 3446.43 26493.56 14323.04 20157.10 39722.86 11458.76

0.08 0.40 0.47 1.87 1.02 3.72 0.80 1.44 10.71 10.82 4.70 18.26 8.87 79.11 50.01 237.64 87.53 213.32 190.03 387.67 2532.80

228.28 708.28 430.58 633.98 1401.72 861.12 3640.00 5068.28 2058.03 2709.12 13280.00 7213.62 10129.52 10783.00 15158.30 3481.21 26512.51 14293.49 20154.32 39703.32 11369.31

0 0 0 0 3 2 3 2 1 6 16 17 63 39 106 54 122 56 299 571 345

228.28 708.28 430.58 631.05 1390.57 831.24 3640.00 5068.28 2044.43 2684.88 13280.00 7216.34 10110.58 10722.73 15102.85 3395.16 26499.23 14202.85 20000.54 39648.42 11357.62

0.05 0.58 0.59 2.16 5.90 5.62 13.74 24.07 35.86 81.76 136.43 179.19 168.07 432.26 658.54 580.27 484.43 676.77 1261.95 1518.12 4326.99

228.28 708.28 430.58 631.05 1390.57 831.24 3640.00 5068.28 2044.20 2684.88 13280.00 7216.60 10110.58 10723.79 15105.90 3394.48 26499.32 14205.07 20007.52 39647.61 11365.37

0.05 0.63 0.62 2.26 6.07 5.81 14.12 24.93 38.78 101.10 152.42 210.71 189.45 479.85 731.98 930.72 577.29 834.60 1524.67 1563.38 5034.56

228.28 708.28 430.58 631.05 1390.57 831.24 3640.00 5068.28 2055.54 2693.19 13280.00 7213.62 10110.58 10770.15 15138.63 3381.29 26493.58 14273.46 20070.07 39690.02 11278.96

0.08 0.09 0.06 0.20 0.62 0.48 1.31 2.28 1.68 4.66 9.45 10.41 19.22 38.16 64.19 45.08 100.40 108.08 194.08 405.10 584.32

9001.91

9072.51 0.47

1004.34

9043.31 0.50

182.92

9038.95 0.63

9009.23 0.07

504.45

9010.17 0.08

591.62

9017.97 0.13

75.71

Boldface denotes the best solutions. a Aleman and Hill (2010). b Berbotto et al. (2014). c Silva et al. (2015). d P4, 3.0 GHz. e P4, 2.8 GHz. f The average of 4 executions. g 4 GB, 2.1 GHz. h The average of 20 executions. i i7, 2.93 GHz. j i7-3770, 3.4 GHz. k The BKS was proved to be optimal by Archetti et al. (2014). l The difference in the BKS of these instances is due to rounding errors in their cost data. m Silva et al. (2015). n Derigs et al. (2010).

81.19

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

SD1_8 SD2_16 SD3_16 SD4_24 SD5_32 SD6_32 SD7_40 SD8_48 SD9_48 SD10_64 SD11_80 SD12_80l SD13_96 SD14_120 SD15_144 SD16_144l SD17_160l SD18_160 SD19_192 SD20_240 SD21_288

BKS

23

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31 Table 7 Results of the SDVRP test for Set B. Set B instances

BKS

S51D2 S51D3 S51D4 S51D5 S51D6 S76D2 S76D3 S76D4 S101D2 S101D3 S101D5 Average Avg. gap (%)

EMIP-MDA + ERTR

TSVBAa

RGTSb

SplitILS-LFc

SplitILS-UFc

SRC + VND

Cost

Timed (s)

Cost

Timee (s)

Costf

Timef,g (s)

Costh

Timeh,i (s)

Costh

Timeh,i (s)

Cost

Timej (s)

708.42k,l 947.97m 1560.88m 1333.67m 2169.10k 1087.40k 1427.81k 2079.76k 1378.43k 1874.81k 2791.22k

717.35 969.18 1580.79 1356.37 2186.29 1105.19 1442.61 2104.87 1397.38 1921.67 2852.01

60.8 56.6 662.6 660.0 677.8 31.4 164.7 980.0 402.0 961.1 709.8

718.69 969.78 1628.20 1362.19 2236.16 1128.15 1472.92 2180.13 1409.03 1947.62 2910.71

31.66 18.75 19.77 15.39 14.38 60.44 51.13 53.56 219.52 132.19 131.16

723.32 970.89 1614.90 1385.31 2215.41 1113.43 1461.20 2103.05 1417.40 1907.92 2898.50

1 4 14 3 2 10 15 14 21 19 14

709.54 949.96 1563.25 1333.85 2174.71 1089.69 1429.01 2080.76 1386.45 1881.26 2795.73

9.98 14.15 59.96 32.41 83.79 74.51 88.72 173.55 129.94 277.62 696.64

709.49 950.12 1563.29 1333.67 2177.78 1089.45 1429.26 2081.16 1386.03 1880.62 2795.36

11.20 15.74 56.28 36.69 62.55 69.36 96.50 188.38 151.66 317.29 572.13

711.42 958.41 1577.56 1339.19 2190.72 1103.37 1444.46 2091.92 1396.11 1915.25 2821.41

3.85 4.34 6.49 5.94 9.16 15.13 18.42 26.15 46.52 45.86 52.12

1578.13

1603.06 1.56

487.89

1633.05 3.22

68.00

1619.21 2.57

10.64

1581.29 0.20

149.21

1581.48 0.21

143.43

1595.44 1.07

21.27

Boldface denotes the best solutions. a Aleman and Hill (2010). b Berbotto et al. (2014). c Silva et al. (2015). d P4, 3.0 GHz. e P4, 2.8 GHz. f The average of 4 executions. g 4 GB, 2.1 GHz. h The average of 20 executions. i i7, 2.93 GHz. j i7-3770, 3.4 GHz. k Silva et al. (2015). l The BKS was proved to be optimal by Archetti et al. (2014). m Archetti et al. (2014).

Regarding the computation speed, the proposed SRC + VND algorithm executes much faster than EMIP-MDA + ERTR. Overall, for Set G, the average computation time across four test groups (shown in Tables 8 and 9) is 47.3 s for SRC + VND compared to approximately 933.0 s for EMIP-MDA + ERTR.

5.3. The SDVRP-MDA test results for Set B Table 10 summarizes the results of our SDVRP-MDA test for Set B. Eleven instances are tested in each of the four groups with different p values, resulting in 44 different benchmark cases tested for Set B. For each of these cases, the BKS is the solution obtained from EMIP-MDA + ERTR (Gulczynski et al., 2010). It is important to note that, for a Set-B instance tested with a specific p, the best solution may be obtained from a different test with a larger p. For example, when an instance is tested with p = 0.1, the best solution for this instance may come from its test with p = 0.2, 0.3, or 0.4. This is because solutions from a more restricted problem are always feasible for the original problem. We find that, in the test group with p = 0.1, there are four instances (i.e., S76D3_75, S76D4_75, S101D3_100, and S101D5_100) whose BKS can be obtained from EMIP-MDA + ERTR tested with p = 0.2. Similarly, in the test group with p = 0.3, there are two instances (i.e., S76D3_75 and S101D2_100) whose BKS can be obtained from EMIP-MDA + ERTR tested with p = 0.4. Accordingly, we have adjusted the BKS for the above-mentioned six instances in Table 10. Let SRC + VNDall denote the implementation of our algorithm with all 120 permutations of the five neighborhood search operators in the VND procedure. The results of SRC + VNDall are also presented in Table 10 for comparison. These results are included because SRC + VNDall yields significantly better solutions than SRC + VND in the test for Set B. However, this is not the case for Set G, where SRC + VNDall performs approximately the same as SRC + VND, with an average gap of 0.014%. The results of SRC + VNDall are thus not shown in Tables 8 and 9. Table 10 shows that SRC + VND and SRC + VNDall both demonstrate a dominant performance over EMIP-MDA + ERTR. Specifically, out of the total of 44 cases tested for Set B, we find 33 improved BKS from SRC + VND, and 42 from SRC + VNDall. The only two exceptions, for which SRC + VNDall cannot find a new best solution, are S51D6_50 with p = 0.3 and S51D4_50 with p = 0.4. As mentioned earlier, the best solution of an instance in Set B can be obtained from a corresponding test for a larger p. We find one such case from our test results. In the test group with p = 0.2, the best solution for S51D5_50 can be obtained from its test for p = 0.3. This new best solution is given in the rightmost column of Table 10 under the label of ‘‘Best solution adjusted”.

24

Table 8 Results of the SDVRP-MDA test for Set G: p = 0.1, 0.2 and 0.3. Set G instance

BKS (Visually estimated solutions)

p = 0.1

p = 0.2

EMIP-MDA + ERTRa

SRC

SRC + VND

p = 0.3

228.28 720.00 430.58 631.05 1402.40 831.24 3588.28 5040.00 2044.20 2684.89 13200.00 7150.58 10042.40 10711.07 15004.22 3631.30 26362.36 14200.92 19964.87 39484.21 11645.47

0.09 0.11 0.11 0.22 0.47 0.47 0.86 1.37 1.47 3.29 6.08 5.99 10.90 23.45 43.09 40.09 52.03 70.28 126.77 219.84 488.91

228.28 720.00 430.58 631.05 1402.40 831.24 3588.28 5060.00 2063.50 2704.89 13280.00 7182.93 10130.60 10733.10 15116.40 3865.24 26519.50 14559.20 20300.40 40102.30 12438.60

228.28 720.00 430.58 640.01 1402.40 831.24 3588.28 5040.00 2080.07 2706.51 13200.00 7150.58 10042.40 10720.04 15049.01 3661.89 26362.36 14271.20 19998.45 39528.99 11691.69

52.18

9137.55 1.17 7

9016.38 0.30 11

9016.38 0.30 11

30.29

8999.92 0.00 21

228.28 720.00 430.58 640.01 1402.40 831.24 3588.28 5040.00 2080.07 2706.51 13200.00 7150.58 10042.40 10720.04 15078.66 3661.89 26362.36 14271.20 20018.45 39558.64 11691.69

0.06 0.11 0.11 0.22 0.44 0.42 0.70 1.12 1.22 2.67 4.76 4.90 8.28 15.90 26.94 28.45 36.07 37.10 64.09 125.41 266.68

228.28 720.00 430.58 631.05 1402.40 831.24 3588.28 5040.00 2044.20 2684.89 13200.00 7150.58 10042.40 10711.07 15004.22 3631.30 26362.36 14200.92 19964.87 39484.21 11645.47

0.06 0.11 0.16 0.23 0.48 0.48 0.86 1.40 1.65 3.48 6.04 6.02 11.51 21.92 43.57 36.97 50.64 57.86 108.30 212.05 428.44

48.26

9158.46 1.42 8

9020.16 0.32 11

Boldface denotes the best solutions. a According to Gulczynski et al. (2010), the computer times of the tests with different p0 s are approximately the same. b P4, 3.0 GHz. c i7-3770, 3.4 GHz. d The difference in the BKS of these instances is due to rounding errors in their cost data.

SRC + VND

30.72

8999.92 0.00 21

SRC + VND

29.79

8999.92 0.00 21

47.25

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

0.08 0.11 0.08 0.20 0.41 0.39 0.69 1.08 1.17 2.56 4.56 4.63 7.94 15.24 26.05 27.55 35.02 36.21 63.93 129.03 279.21

933.03

228.28 720.00 430.58 631.05 1402.40 831.24 3588.28 5040.00 2063.50 2710.64 13334.10 7170.58 10112.40 10836.30 15172.10 3962.67 26646.50 14420.20 20355.70 40018.30 12652.90

228.28 720.00 430.58 631.05 1402.40 831.24 3588.28 5040.00 2044.20 2684.89 13200.00 7150.58 10042.40 10711.07 15004.22 3631.30 26362.36 14200.92 19964.87 39484.21 11645.47

228.28 720.00 430.58 640.01 1402.40 831.24 3588.28 5040.00 2080.07 2706.51 13200.00 7150.58 10042.40 10720.04 15049.01 3661.89 26362.36 14271.20 19998.45 39528.99 11691.69

9129.06 0.93 6

0.09 0.12 0.12 0.23 0.47 0.47 0.83 1.42 1.45 3.31 6.04 6.35 10.83 23.40 41.50 38.06 52.59 65.89 118.44 211.35 430.53

0.08 0.06 0.11 0.20 0.42 0.39 0.69 1.09 1.19 2.70 4.66 4.79 8.21 16.02 27.30 29.09 36.43 37.21 65.10 128.45 280.93

Timec (s)

8999.92

Timec (s)

Timec (s)

Cost

Average Avg. gap (%) BKS found

Cost

Cost

Timec (s)

1.06 291.81 5.52 227.62 644.03 644.95 645.52 646.91 645.16 648.84 662.44 659.81 667.02 1644.98 1640.47 1625.83 1654.72 1639.72 1649.31 1691.84 1656.14

Timec (s)

Timec (s)

Cost

228.28 720.00 430.58 631.05 1402.43 831.24 3588.28 5060.00 2074.12 2691.69 13220.00 7182.93 10111.80 10845.90 15180.70 3755.70 26628.40 14477.80 20432.20 40202.50 12014.60

Cost

Cost

Timeb (s)

228.28 720.00 430.58 631.05 1402.40 831.24 3588.28 5040.00 2044.20 2684.88 13200.00 7150.58 10042.40 10711.06 15004.20 3631.30 26362.46 14200.90 19964.86 39484.23 11645.50

SRC

SRC

Cost MD1_8 MD2_16 MD3_16 MD4_24 MD5_32 MD6_32 MD7_40 MD8_48 MD9_48 MD10_64d MD11_80 MD12_80 MD13_96 MD14_120 MD15_144d MD16_144 MD17_160d MD18_160d MD19_192d MD20_240d MD21_288d

EMIPMDA + ERTR Cost

EMIPMDA + ERTR Cost

25

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31 Table 9 Results of the SDVRP-MDA test for Set G: p = 0.4. Set G instance

BKS

Estimated solutions

p = 0.4 EMIP-MDA + ERTRa

MD1_8 MD2_16 MD3_16 MD4_24 MD5_32 MD6_32 MD7_40 MD8_48 MD9_48 MD10_64 MD11_80 MD12_80 MD13_96 MD14_120d MD15_144d MD16_144d MD17_160d MD18_160d MD19_192 MD20_240d MD21_288 Average Avg. gap (%) BKS improved

SRC

SRC + VND

Cost

Timeb (s)

Cost

Timec (s)

Cost

Timec (s)

228.28 720.00 430.58 631.05 1402.40 831.24 3588.28 5040.00 2044.20 2684.88 13200.00 7150.58 10042.40 10711.07 15004.22 3443.69 26362.35 14200.92 20032.12 39484.21 11457.88

0.06 0.08 0.09 0.17 0.45 0.39 0.90 1.58 1.25 2.95 6.91 6.88 12.60 22.48 39.42 30.42 56.88 52.03 97.58 199.45 340.04

228.28 720.00 430.58 631.05 1402.40 830.26 3588.28 5040.00 2044.20 2684.88 13200.00 7150.58 10042.40 10711.06 15004.20 3445.50 26362.46 14200.90 19964.86 39484.23 11645.50

228.28 720.00 430.58 631.05 1402.40 831.24 3588.28 5040.00 2044.20 2684.88 13200.00 7150.58 10042.40 10711.06 15004.20 3631.30 26362.46 14200.90 19964.86 39484.23 11645.50

228.28 720.00 430.58 631.05 1414.75 830.26 3588.28 5040.00 2059.03 2708.80 13240.00 7260.01 10233.50 10865.10 15202.80 3445.50 26904.70 14447.60 20608.90 40551.40 11909.10

1.06 291.81 5.52 227.62 644.03 644.95 645.52 646.91 645.16 648.84 662.44 659.81 667.02 1644.98 1640.47 1625.83 1654.72 1639.72 1649.31 1691.84 1656.14

228.28 753.14 430.58 640.01 1402.40 831.24 3810.26 5363.84 2080.07 2684.89 13796.98 7695.17 10592.22 11057.13 15796.02 3661.89 27456.64 14724.31 20826.25 41220.15 11701.90

0.06 0.08 0.09 0.17 0.39 0.33 0.69 1.12 1.09 2.42 5.12 4.99 8.69 16.26 27.77 23.84 38.55 37.67 65.13 144.55 278.59

8991.03

8999.92 0.26 –

9158.08 1.00

933.03

9369.21 3.33 0

31.31



8985.25 0.06 1

41.55

Boldface denotes the best solutions, and the improved BKS is underlined. a The computer times of the tests with different p0 s were approximately the same. b P4, 3.0 GHz. c i7-3770, 3.4 GHz. d The difference in the BKS of these instances is due to rounding errors in their cost data.

The computation times associated with EMIP-MDA + ERTR for Set B are not available in the literature; thus, we cannot conduct a formal comparison for time efficiency. However, based on the results from Set G, it is reasonable to expect that SRC + VND would still execute faster than EMIP-MDA + ERTR for Set B. Regarding the computation time, we also find that the proposed SRC + VND behaves differently from EMIP-MDA + ERTR. As described in Gulczynski et al. (2010), the computation times associated with different p values are approximately the same. However, our results indicate that the average computation times for SRC + VND (and SRC + VNDall) exhibit a downward trend as p increases. Specifically, in Tables 8 and 9 (Set G), the average computation times for SRC + VND are 52.2, 48.3, 47.3, and 41.6 s for p = 0.1, 0.2, 0.3, and 0.4, respectively. Similarly, in Table 10 (Set B), these times are 21.7, 18.9, 16.6, and 14.7 s for p = 0.1, 0.2, 0.3, and 0.4, respectively. This behavior suggests that SRC + VND (and SRC + VNDall) may work more efficiently for problems with more restricted minimal delivery fractions, i.e., larger p values. For reference, we present the detailed results of the new best solution that we have found for S76D4_75 (p = 0.2) in Appendix B. On https://db.tt/40JBGP7w, we offer the detailed results of the solutions that we have found for all the 128 cases tested for 32 benchmark instances (in Sets G and B) with four p values. As mentioned earlier, of the 14 instances in Set B introduced by Belenguer et al. (2000) for the SDVRP, only 11 instances were tested for the SDVRP-MDA by Gulczynski et al. (2010). The three untested instances are S51D1_50, S76D1_75, and S101D1_100. For completeness, we have also tested these three instances for p = 0.1, 0.2, 0.3, and 0.4. The results of this additional test are given in Appendix C as a benchmark for future research. 5.4. Results from different implementations with SRC0 and SRC1 We have used two procedures (i.e., SRC0 and SRC1) in the SRC constructive algorithm. For every instance tested, the better solution obtained from SRC0 and SRC1 defines the initial solution for SRC. Similarly, the better solution of SRC0 + VND and SRC1 + VND defines the final solution for SRC + VND. In this section, we analyze the effects of SRC0 and SRC1 on different implementations with VND or VNDall. The results of these different implementations are summarized in Table 11. We first examine the initial solution performance related to SRC0, SRC1, and SRC. While the two constructive procedures SRC0 and SRC1 consume approximately the same computation time, SRC0 demonstrates a dominant performance over SRC1. For all instances tested with four p values, the average gap for SRC0 is 3.70%, which is much lower than that for SRC1, (i.e., 5.97%). The SRC0 performs particularly well for Set G with an average gap of 1.22%, suggesting that SRC0 tends to work more effectively for problems with symmetrically distributed nodes.

26

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

Table 10 Results of the SDVRP-MDA test for Set B: p = 0.1, 0.2, 0.3 and 0.4. p

0.1

0.2

0.3

0.4

Set B instance

BKS

EMIP-MDA + ERTRa

SRC + VND

Cost

Gap (%)

Cost

Gap (%)

Timeb (s)

Cost

Gap (%)

Timeb (s)

SRC + VNDall

Best solution adjusted

S51D2_50 S51D3_50 S51D4_50 S51D5_50 S51D6_50 S76D2_75 S76D3_75 S76D4_75 S101D2_100 S101D3_100 S101D5_100

717.35 969.99 1588.91 1373.98 2225.51 1106.86 1446.48c 2118.86c 1398.13 1929.96c 2862.14c

717.35 969.99 1588.91 1373.98 2225.51 1106.86 1457.40 2123.16 1398.13 1930.86 2862.34

0.00 0.00 0.00 0.00 0.00 0.00 0.75 0.20 0.00 0.05 0.01

713.30 953.17 1569.89 1344.68 2197.64 1100.90 1448.39 2095.95 1407.35 1900.60 2851.33

0.57 1.73 1.20 2.13 1.25 0.54 0.13 1.08 0.66 1.52 0.38

3.99 4.84 6.68 5.49 7.96 15.85 20.78 25.96 48.20 46.29 53.01

709.77 949.78 1569.89 1335.55 2197.04 1095.53 1440.56 2095.95 1395.42 1899.66 2839.49

1.06 2.08 1.20 2.80 1.28 1.02 0.41 1.08 0.19 1.57 0.79

439.65 561.03 775.82 666.33 916.53 1782.50 2362.41 3040.63 5600.83 5404.45 6325.75

Average BKS improved

1612.56 –

1614.05 –

0.09

1598.47 9

0.87

21.73

1593.51 11

1.23

2534.18

S51D2_50 S51D3_50 S51D4_50 S51D5_50 S51D6_50 S76D2_75 S76D3_75 S76D4_75 S101D2_100 S101D3_100 S101D5_100

717.35 969.99 1593.69 1377.99 2285.37 1116.64 1446.48 2118.86 1398.87 1929.96 2862.14

717.35 969.99 1593.69 1377.99 2285.37 1116.64 1446.48 2118.86 1398.87 1929.96 2862.14

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

712.59 959.49 1589.32 1359.22 2240.14 1116.07 1450.04 2112.62 1401.57 1917.56 2866.29

0.66 1.08 0.27 1.36 1.98 0.05 0.25 0.29 0.19 0.64 0.14

3.79 4.26 6.19 4.66 5.07 13.53 17.07 22.89 44.32 43.31 42.96

709.77 955.45 1578.15 1351.12 2240.14 1104.27 1445.63 2103.78 1396.78 1904.09 2861.88

1.06 1.50 0.98 1.95 1.98 1.11 0.06 0.71 0.15 1.34 0.01

424.37 511.86 709.12 553.29 619.05 1621.75 2095.96 2837.47 5544.72 5319.77 5379.26

Average BKS improved

1619.76 –

1619.76 –

0.00

1611.36 8

0.52

18.91

1604.64 11

0.99

2328.78

S51D2_50 S51D3_50 S51D4_50 S51D5_50 S51D6_50 S76D2_75 S76D3_75 S76D4_75 S101D2_100 S101D3_100 S101D5_100

717.35 969.99 1597.89 1383.71 2301.51 1116.64 1453.17d 2151.49 1398.88d 1939.96 2901.76

717.35 969.99 1597.89 1383.71 2301.51 1116.64 1453.25 2151.49 1401.85 1939.96 2901.76

0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.21 0.00 0.00

716.17 966.68 1597.10 1356.64 2307.14 1118.41 1446.81 2123.06 1407.19 1921.27 2887.90

0.17 0.34 0.05 1.96 0.24 0.16 0.44 1.32 0.59 0.96 0.48

3.09 4.32 4.43 4.02 2.84 14.20 15.23 17.88 40.76 37.44 38.45

711.56 959.16 1595.40 1349.55 2307.14 1108.20 1446.81 2123.06 1396.99 1918.41 2873.36

0.81 1.12 0.16 2.47 0.24 0.76 0.44 1.32 0.14 1.11 0.98

370.07 508.91 528.17 488.59 357.24 1676.14 1912.85 2255.93 5108.96 4678.31 4691.72

Average BKS improved

1630.21 –

1630.49 –

0.02

1622.58 8

0.43

16.61

1617.24 10

0.82

2052.44

S51D2_50 S51D3_50 S51D4_50 S51D5_50 S51D6_50 S76D2_75 S76D3_75 S76D4_75 S101D2_100 S101D3_100 S101D5_100

717.35 978.41 1612.30 1389.32 2402.35 1116.64 1453.17 2167.27 1398.88 1942.94 2904.61

717.35 978.41 1612.30 1389.32 2402.35 1116.64 1453.17 2167.27 1398.88 1942.94 2904.61

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

717.31 962.21 1619.87 1373.82 2395.33 1113.78 1453.43 2143.90 1403.70 1935.08 2898.98

0.01 1.66 0.47 1.12 0.29 0.26 0.02 1.08 0.34 0.40 0.19

3.10 3.70 3.40 3.32 1.37 13.51 14.21 16.50 39.42 31.09 32.15

713.60 959.16 1613.06 1360.01 2395.33 1112.58 1446.88 2129.97 1397.75 1922.10 2890.76

0.52 1.97 0.05 2.11 0.29 0.36 0.43 1.72 0.08 1.07 0.48

367.89 445.00 405.11 399.73 174.85 1656.92 1712.80 2035.99 4862.50 3938.04 3950.67

Average BKS improved

1643.93 –

1643.93 –

0.00

1637.95 8

0.38

14.71

1631.02 10

0.82

1813.59

Boldface denotes the best solutions, and each improved BKS is underlined. a The computation time is not given. b i7-3770, 3.4 GHz. c From EMIP-MDA + ERTR for p = 0.2. d From EMIP-MDA + ERTR for p = 0.4. e From SRC + VNDall for p = 0.3.

1349.55e

27

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31 Table 11 Results of initial and final solutions from different implementations. Set G

Set B

Gap (%) Initial solutions SRC0 SRC1 SRC

Time (s)

All instances

Gap (%)

Time (s)

Gap (%)

Time (s)

1.22 2.99 1.07

15.17 15.36 30.53

8.43 11.67 8.28

1.04 1.03 2.07

3.70 5.97 3.55

10.31 10.43 20.75

Final solutions SRC0 + VND SRC1 + VND SRC + VND

0.01 0.03 0.01

23.16 24.15 47.31

0.24 0.40 0.55

8.51 9.47 17.99

0.09 0.12 0.20

18.12 19.10 37.23

SRC0 + VNDall SRC1 + VNDall SRC + VNDall

0.01 0.02 0.01

2802.02 2915.81 5717.83

0.69 0.90 0.96

1040.51 1141.74 2182.25

0.24 0.30 0.34

2196.50 2305.97 4502.47

10 SRC0 SRC1 8

SRC Phase I: Initial Solutions

SRC0+VND

5.97 %

SRC1+VND

Average gap (%)

6

SRC+VND

3.70 %

4

3.55 % 6.09 %

2

3.79 % Phase II: Final Solutions

-0.09 %

-0.20 %

0 -0.12 % -2 0

10

20

30

40

Average CPU time (seconds) Fig. 5. Competitiveness and complementarity of SRC0 + VND and SRC1 + VND.

Although SRC1 performs poorly compared to SRC0, and contributes very little in Phase I for the initial solution construction, this situation, however, changes dramatically in Phase II. Specifically, the final solution obtained from SRC1 + VND becomes quite competitive compared with that from SRC0 + VND. In fact, for all the instances tested with four different p values, SRC1 + VND even outperforms SRC0 + VND, and the average gap for SRC1 + VND (i.e., 0.12%) is even lower than that for SRC0 + VND (i.e., 0.09%). This come-from-behind performance of SRC1 + VND is illustrated in Fig. 5. The improvement in the average gap from the initial to the final solutions for SRC1 + VND is 6.09% compared to that of 3.79% for SRC0 + VND. Another interesting finding is the complementarity of SRC0 and SRC1 in the VND implementation for producing final solutions. By the solution complementarity of two procedures, we mean that one procedure is likely to produce good solutions, while the other is not. Accordingly, when these two procedures are combined for use, they tend to produce much better solutions than when they are used alone. In our case, when SRC0 + VND and SRC1 + VND are implemented alone, the average gaps (for all cases tested) are 0.09% and 0.12%, respectively. When both of them are used in SRC + VND, they yield better results with an average gap of 0.20%. Although the absolute value of the difference in these two gaps appears to be insignificant, the relative improvement is actually remarkable. Specifically, the relative improvement from 0.12% to 0.20% is 66.7%, which is significant. The results of different implementations with VNDall are also given in Table 11. These results show that implementations with VNDall yield better solutions, but require much more execution time than VND. For all instances tested with four p values, the average gap is 0.34% and the average time is 4502.5 s for SRC + VNDall, compared to that gap of 0.20% and time of 37.2 s for SRC + VND. For the final solutions of SRC1 + VNDall, the come-from-behind behavior (over SRC0 + VNDall) associated with SRC1 still prevails, but the solution complementarity of SRC0 + VNDall and SRC1 + VNDall becomes less significant than that of SRC0 + VND and SRC1 + VND.

28

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

Overall, the results in Table 11 indicate that the proposed SRC + VND implementation seems to be the most timeeffective, which yields solutions with an average gap of 0.20% in 37.2 s. If accuracy is of greater concern than time, then SRC1 + VNDall may provide a good alternative, which is capable of generating nearly the same accuracy of results as SRC + VNDall with an average gap of 0.30% in about half the time of SRC + VNDall, i.e., 2306 s. On the other hand, if time is of greater concern, then SRC1 + VND tends to be a better choice, which can provide good results with an average gap of 0.12% in approximately only 19 s. 6. Conclusions Split deliveries, while beneficial to in-house logistic operations, usually incur extra work and more distraction to outside customers. Good solution methods for the SDVRP-MDA tend to provide a way to resolve this dilemma. Nevertheless, thus far, the research devoted to the SDVRP-MDA has been quite limited and much less than that conducted on the SDVRP.

Table A The pseudocode of SNEC.

29

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

In this paper, we propose a multi-start variable neighborhood search approach for solving the SDVRP-MDA. The proposed solution algorithms are tested with all 32 available benchmark instances for four different minimum delivery fractions (i.e., p = 0.1, 0.2, 0.3, and 0.4). Out of these 128 cases tested, we have found 81 and improved 34 best known solutions with the proposed SRC + VND algorithm in an average of about 37.2 s. The proposed algorithm is also tested for the SDVRP with competitive performance compared to several of the best algorithms available in the literature. The results of this study suggest that a relatively simple multi-start heuristic with good search operators can be as effective as sophisticated metaheuristics. For future work, we expect that the proposed constructive algorithm SRC and operator SNEC can be applied with other metaheuristics to develop more practical solution tools for the SDVRP-MDA and other related problems.

Table B The new best solution of S76D4_75 with p = 0.2. Benchmark set: B Instance: S76D4_75 BKS cost: 2118.86 Our best solution: 2103.78 Personal computer: 3.4 GHz i7-3770 processor, 16 GB of RAM, Windows7

a

Route no.

Customer (delivery amount) serviced in sequence

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

48(7)a, 47(18)a, 36(21)a, 60(98), 37(16)a, 0 33(14)a, 64(121), 22(25), 0 25(39)a, 55(55), 18(26), 50(35), 0 20(37), 70(87), 37(36)a, 0 36(21)a, 71(139), 0 14(29), 59(101), 11(30)a, 0 47(68)a, 69(89), 0 23(35), 56(73), 1(52)a, 0 58(27)a, 10(33), 31(100), 0 6(17)a, 43(126), 33(17)a, 0 11(80)a, 66(72), 38(8)a, 0 28(37)a, 61(121), 0 40(24)a, 25(85)a, 9(51), 0 3(139), 0 57(42)a, 15(118), 0 2(26)a, 62(60), 73(74), 0 58(10)a, 38(17)a, 65(105), 11(28)a, 0 13(65)a, 57(93)a, 0 51(18)a, 49(141), 0

Minimum delivery fraction: p = 0.2 Algorithm: SRC + VNDall Computational time: 2255.93 s

Route no.

Customer (delivery amount) serviced in sequence

20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

30(5)a, 21(143), 28(10)a, 0 13(37)a, 54(41), 19(82), 0 48(21)a, 5(36), 29(84), 4(18)a, 0 2(20)a, 74(103), 30(19)a, 4(13)a, 0 8(17)a, 53(143), 0 67(143), 0 12(17)a, 39(143), 0 6(22)a, 63(131), 0 51(14)a, 16(25), 24(44), 44(20), 32(47), 40(9)a, 0 58(12)a, 72(143), 12(5)a, 0 68(138), 0 8(24)a, 35(133), 0 34(24), 46(31), 52(49), 27(50), 0 4(34)a, 45(126), 0 7(143), 0 26(24), 17(135), 0 75(124), 0 41(96), 42(27), 1(34)a, 0

This denotes a split delivery.

Table C Additional test results of three instances for Set B. p

0.1

0.2

0.3

0.4

Set B instance

SRC + VND

SRC + VNDall

Best solution adjusted

Cost

Timea (s)

Cost

Timea (s)

S51D1_50 S76D1_75 S101D1_100 Average

462.45 611.36 734.75 602.85

1.47 12.68 35.12 16.42

459.50 602.75 732.06 598.10

175.80 1548.46 4243.33 1989.20

S51D1_50 S76D1_75 S101D1_100 Average

462.45 605.82 735.86 601.38

1.44 13.14 33.63 16.07

459.50 600.19 730.87 596.85

179.38 1533.65 4089.32 1934.12

S51D1_50 S76D1_75 S101D1_100 Average

462.45 603.36 738.44 601.41

1.45 12.65 31.59 15.23

459.50 602.75 730.11 597.45

186.76 1521.21 3757.31 1821.76

S51D1_50 S76D1_75 S101D1_100 Average

459.50 611.02 738.42 602.98

1.59 11.61 29.75 14.32

459.50 600.19 731.27 596.99

194.15 1416.64 3699.46 1770.08

Boldface denotes the best solutions. a i7-3770, 3.4 GHz. b From SRC + VNDall for p = 0.2. c From SRC + VNDall for p = 0.3. d From SRC + VNDall for p = 0.4.

600.19b,d 730.11c

730.11c

600.19d

30

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

Appendix A Table A presents the pseudocode of SNEC. In Section 4.2, SNEC is described with nodes i and j, and routes r 1 , r2 , r 3 . For a given r1 and a given i (lines 2–3), the procedure checks for all possible r 2 (in terms of w1 as in line 5) and, if necessary, all j (in terms of v as in line 14) and r3 (in terms of w2 as in line 16) to improve the solution. The results of all these checks are classified into four categories with different Shifttype values, corresponding to those scenarios depicted in Fig. 3. Specifically, for Shifttype = 0 (Fig. 3a1 and a2), no improvement can be found by ejecting i from r1 to w1 (lines 6, 7 and 12); for Shifttype = 1 (Fig. 3b), the total cost is improved with a feasible move of ejecting i from r 1 to w1 (lines 7–11); for Shifttype = 2 (Fig. 3c1), the total cost is improved by ejecting i from r1 to w1 , and shifting a node v entirely (with all its qv w1 ) from w1 to w2 (lines 14–23); and (4) for Shifttype = 3 (Fig. 3c2), the total cost is improved by ejecting i from r1 to w1 , and shifting a node v partly (with a part of its qv w1 ) from w1 to w2 (lines 14–19, and lines 24–30). For each check, the incumbent solution is updated with r 1 , r 2 , and r3 accordingly (lines 31–46). Finally, the best neighborhood solution of SNEC is returned (line 47), when a complete iteration of all loops with r 1 , i, w1 , v and w2 has been checked. The most complicated part of the SNEC procedure is to verify the case of Shifttype = 3, i.e., the case for shifting part of qv w1 from w1 to w2 . In this case, we must consider both the maximum amount which can be shifted out of w1 , i.e., min(SQ w2 , qv w1  MDAv ), and the minimum amount that must be transferred to w2 , i.e., max(Ow1 , MDAv  qv w2 ), where SQ w2 and Ow1 are the spare capacity of w2 , and the overloaded amount of w1 , respectively. Only when the former is no less than the latter (line 25), the procedure proceeds to validate Shifttype = 3 (lines 26–30). Appendix B Table B presents the detailed results of the improved BKS that we have found for S76D4_75 with p = 0.2. In this solution, the maximum number of nodes that are visited more than once on a single route is four. These four nodes are 36, 37, 47 and 48 on route 1. The maximal number of split deliveries for a single node is three, which occurs for nodes 11 and 58, i.e., both nodes are visited three times on three different routes. Appendix C Table C summarizes our test results of those three Set-B instances that were not tested by Gulczynski et al. (2010). Each of these instances is tested with four p values, i.e., p = 0.1, 0.2, 0.3 and 0.4. References 7-ELEVEN Taiwan, 2015. President Chain Store Inc., Takeout service order . Ahuja, R.K., Magnanti, T.L., Orlin, J.B., 1993. Network Flows: Theory, Algorithms, and Applications. Prentice-Hall, Inc. Aleman, R.E., Hill, R.R., 2010. A tabu search with vocabulary building approach for the vehicle routing problem with split demands. Int. J. Metaheurist. 1 (1), 55–80. Aleman, R.E., Zhang, X.H., Hill, R.R., 2010. An adaptive memory algorithm for the split delivery vehicle routing problem. J. Heurist. 16 (3), 441–473. Archetti, C., Speranza, M.G., 2012. Vehicle routing problems with split deliveries. Int. Trans. Oper. Res. 19 (1–2), 3–22. Archetti, C., Savelsbergh, M.W.P., Speranza, M.G., 2006a. Worst-case analysis for split delivery vehicle routing problems. Transport. Sci. 40 (2), 226–234. Archetti, C., Speranza, M.G., Hertz, A., 2006b. A tabu search algorithm for the split delivery vehicle routing problem. Transport. Sci. 40 (1), 64–73. Archetti, C., Savelsbergh, M.W.P., Speranza, M.G., 2008a. To split or not to split: that is the question. Transport. Res. Part E 44 (1), 114–123. Archetti, C., Speranza, M.G., Savelsbergh, M.W.P., 2008b. An optimization-based heuristic for the split delivery vehicle routing problem. Transport. Sci. 42 (1), 22–31. Archetti, C., Bianchessi, N., Speranza, M.G., 2014. Branch-and-cut algorithms for the split delivery vehicle routing problem. Eur. J. Oper. Res. 238 (3), 685– 698. Belenguer, J.M., Martinez, M.C., Mota, E., 2000. A lower bound for the split delivery vehicle routing problem. Oper. Res. 48 (5), 801–810. Berbotto, L., García, S., Nogales, F., 2014. A randomized granular tabu search heuristic for the split delivery vehicle routing problem. Ann. Oper. Res. 222 (1), 153–173. Boudia, M., Prins, C., Reghioui, M., 2007. An effective memetic algorithm with population management for the split delivery vehicle routing problem. In: Bartz-Beielstein, T., Blesa Aguilera, M., Blum, C., Naujoks, B., Roli, A., Rudolph, G., Sampels, M. (Eds.), Hybrid Metaheuristics. Springer, Berlin, Heidelberg, pp. 16–30. Campos, V., Corberán, A., Mota, E., 2008. A scatter search algorithm for the split delivery vehicle routing problem. In: Fink, A., Rothlauf, F. (Eds.), Advances in Computational Intelligence in Transport, Logistics, and Supply Chain Management. Springer, Berlin, Heidelberg, pp. 137–152. Chen, S., Golden, B., Wasil, E., 2007. The split delivery vehicle routing problem: applications, algorithms, test problems, and computational results. Networks 49 (4), 318–329. Cordeau, J.-F., Laporte, G., 2005. Tabu search heuristics for the vehicle routing problem. In: Sharda, R., Voß, S., Rego, C., Alidaee, B. (Eds.), Metaheuristic Optimization via Memory and Evolution. Springer, US, pp. 145–163. Derigs, U., Li, B., Vogel, U., 2010. Local search-based metaheuristics for the split delivery vehicle routing problem. J. Oper. Res. Soc. 61 (9), 1356–1364. Domino’s Pizza Singapore Pte. Ltd., 2015. Terms of use . Dror, M., Trudeau, P., 1989. Savings by split delivery routing. Transport. Sci. 23 (2), 141–145. Dror, M., Trudeau, P., 1990. Split delivery routing. Nav. Res. Logist. 37 (3), 383–402. Everhealth Pharmaceuticals Co. Ltd., 2015. News . Glover, F., 1992. New ejection chain and alternating path methods for traveling salesman problems. In: Balci, O., Sharda, R., Zenios, S.A. (Eds.), Computer Science and Operations Research: New Developments in Their Interfaces. Pergamon, Oxford, pp. 491–507. Groër, C., Golden, B., Wasil, E., 2011. A parallel algorithm for the vehicle routing problem. INFORMS J. Comput. 23 (2), 315–330. Gulczynski, D., 2010. Integer Programming-Based Heuristics for Vehicle Routing Problems. University of Maryland. Gulczynski, D., Golden, B., Wasil, E., 2010. The split delivery vehicle routing problem with minimum delivery amounts. Transport. Res. Part E 46 (5), 612– 626.

A.F.-W. Han, Y.-C. Chu / Transportation Research Part E 88 (2016) 11–31

31

Ho, S.C., Haugland, D., 2004. A tabu search heuristic for the vehicle routing problem with time windows and split deliveries. Comput. Oper. Res. 31 (12), 1947–1964. Jin, M., Liu, K., Bowden, R.O., 2007. A two-stage algorithm with valid inequalities for the split delivery vehicle routing problem. Int. J. Prod. Econ. 105 (1), 228–242. Jin, M., Liu, K., Eksioglu, B., 2008. A column generation approach for the split delivery vehicle routing problem. Oper. Res. Lett. 36 (2), 265–270. Li, F., Golden, B., Wasil, E., 2005. Very large-scale vehicle routing: new test problems, algorithms, and results. Comput. Oper. Res. 32 (5), 1165–1179. Lin, S., Kernighan, B.W., 1973. An effective heuristic algorithm for the traveling-salesman problem. Oper. Res. 21 (2), 498–516. Mladenovic´, N., Hansen, P., 1997. Variable neighborhood search. Comput. Oper. Res. 24 (11), 1097–1100. Mota, E., Campos, V., Corberán, Á., 2007. A new metaheuristic for the vehicle routing problem with split demands. In: Cotta, C., Hemert, J. (Eds.), Evolutionary Computation in Combinatorial Optimization. Springer, Berlin, Heidelberg, pp. 121–129. Or, I., 1976. Traveling Salesman-Type Combinatorial Problems and Their Relation to the Logistics of Regional Blood Banking. Northwestern University. Osman, I., 1993. Metastrategy simulated annealing and tabu search algorithms for the vehicle routing problem. Ann. Oper. Res. 41 (4), 421–451. Potvin, J.-Y., Rousseau, J.-M., 1995. An exchange heuristic for routing problems with time windows. J. Oper. Res. Soc. 46 (12), 1433–1446. Rego, C., 2001. Node-ejection chains for the vehicle routing problem: sequential and parallel algorithms. Parallel Comput. 27 (3), 201–222. Silva, M.M., Subramanian, A., Ochi, L.S., 2015. An iterated local search heuristic for the split delivery vehicle routing problem. Comput. Oper. Res. 53 (0), 234– 249. Tesco PLC., 2015. Delivery Options . Xiong, Y., Gulczynski, D., Kleitman, D., Golden, B., Wasil, E., 2013. A worst-case analysis for the split delivery vehicle routing problem with minimum delivery amounts. Optimiz. Lett. 7 (7), 1597–1609. Yellow, P.C., 1970. A computational modification to the savings method of vehicle scheduling. Oper. Res. Quart. 21 (2), 281–283.