A hybrid ant colony optimization algorithm for a multi-objective vehicle routing problem with flexible time windows

Accepted Manuscript

A hybrid ant colony optimization algorithm for a multi-objective vehicle routing problem with flexible time windows Huizhen Zhang, Qinwan Zhang, Liang Ma, Ziying Zhang, Yun Liu PII: DOI: Reference:

S0020-0255(16)32120-X https://doi.org/10.1016/j.ins.2019.03.070 INS 14410

To appear in:

Information Sciences

Received date: Revised date: Accepted date:

20 December 2016 23 March 2019 26 March 2019

Please cite this article as: Huizhen Zhang, Qinwan Zhang, Liang Ma, Ziying Zhang, Yun Liu, A hybrid ant colony optimization algorithm for a multi-objective vehicle routing problem with flexible time windows, Information Sciences (2019), doi: https://doi.org/10.1016/j.ins.2019.03.070

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

A hybrid ant colony optimization algorithm for a multi-objective vehicle routing problem with flexible time windows Qinwan Zhang†

Liang Ma‡

March 22, 2019

Abstract

Ziying Zhang§

Yun Liu¶k

CR IP T

Huizhen Zhang∗

M

AN US

Abstract: In this paper, we present a multi-objective vehicle routing problem with flexible time windows (MOVRPFlexTW). In this problem, a fleet of vehicles can service a set of customers earlier and later than the required time with a given tolerance. This flexibility enables a logistics company to save distribution costs at the expense of customer satisfaction. This paper uses a direct interpretation of the vehicle routing problem with flexible time windows (VRPFlexTW) as a multi-objective problem, where the total distribution costs (including travel costs and fixed vehicle costs) are minimized and the overall customer satisfaction is maximized. We propose a solution strategy based on ant colony optimization and three mutation operators, which incorporates the concept of Pareto optimality for multi-objective optimization. The performance of the proposed approach was evaluated using the well-known benchmark Solomon’s problems. Our experimental results show that the suggested approach is effective, because it provides solutions that are comparative to the best known results in the literatures. Finally, we not only highlight the advantages of MOVRPFlexTW when compared with the single objective VRPFlexTW, but also illustrate its applicability and validation by analyzing a real case.

Introduction

PT

1

ED

Key words: Vehicle routing problem, Flexible time windows, Ant colony optimization algorithm, Multi-Objective Optimization, Mutation operator.

AC

CE

The vehicle routing problem (VRP) is a well-known NP-hard combinatorial problem arising in transportation, logistics, scheduling and supply chain management. The basic VRP determines an optimal set of routes for a set of vehicles to deliver goods and/or services from a supply depot to a set of geographically dispersed customers in a cost-effective manner. As the VRP appears in real life, it may have several classes of additional constraints, such as limits on the vehicle capacity [46, 48], time windows for serving customers [3, 4, 6, 5], route lengths, or the number of hours worked by a driver or a distribution clerk. For a recent complete review on the classification of different VRP variants, see [29, 9]. As with basic VRP, most VRP variants are known to be NP-hard. The capacitated vehicle routing problem with (hard) time windows (VRPTW) is an important variant of the VRP, and has been drawing extensive attention from researchers. It aims to find the optimal set ∗

Corresponding author: [email protected], School of Management,University of Shanghai for Science and Technology, Mail box 459, 516 Jungong Road, Shanghai 200093,China † School of Management,University of Shanghai for Science and Technology ‡ School of Management,University of Shanghai for Science and Technology § School of Materials Engineering, Shanghai University of Engineering Science ¶ School of Management,University of Shanghai for Science and Technology

1

ACCEPTED MANUSCRIPT

of routes for a set of identical vehicles with restricted capacity so that all customers, whose demands are known, are serviced exactly once by a single vehicle within each time window defined by the earliest and latest times. Moreover, each vehicle should start and end its route at the given depot. The vehicle is permitted to arrive at the customer point before the earliest time, and wait at no cost until service becomes allowable, but it is not permitted to arrive at the customer point after the latest time.

CR IP T

The VRPTW requires the time window constraints to be strictly satisfied. However, in some real-world situations( e.g. newspaper delivery, non-emergency parcel delivery, furniture delivery), the time window constraints can be violated to a certain extent. These violations of time windows are considered in the vehicle routing problem with soft time windows (VRPSTW) and the vehicle routing problem with flexible time windows (VRPFlexTW). The VRPSTW assumes that some or all customer time windows can be violated [2], and customers are available to be serviced at any point of the planning horizon [43]. The VRPFlexTW assumes that some or all customer time windows can be violated to a certain extent [43]. The VRPFlexTW allows the vehicles to serve the customers before and after the earliest and latest times with a given tolerance. Therefore, the VRPSTW can be considered as a special case of the VRPFlexTW where the violation of time windows is unbounded, i.e., infinite flexible bounds [43].

PT

ED

M

AN US

The violation of time window constraints may lead to reducing the total distance traveled and the number of vehicles used for the service at the expense of customer satisfaction [43]. Thus, when minimizing the routing costs and maximizing the customer satisfaction (or minimizing the customer inconvenience) affected by the time window violations are defined, these two different objectives frequently conflict. Customer inconvenience for being visited too late and (possibly) too early with respect to the desired time window is typically quantified and the objective functions in VRPSTW and VRPFlexTW are generally modeled as a weighted combination of routing costs and a penalty corresponding to the customer inconvenience [41]. However, how to measure the customer inconvenience and quantify the relative weights of routing and customer inconvenience are still open research questions [41]. For these reasons, adopting a multi-objective point of view can be advantageous and the optimization process needs to provide a range of solutions that represent the trade-offs between the objectives [21], rather than a single solution. Therefore, this article studies the multi-objective vehicle routing problem with flexible time windows (MOVRPFlexTW). This problem is based on the VRPFlexTW proposed in [43] which aggregates the travel costs, fixed costs for using the vehicles to delivery service, and penalty costs incurred for early and late service into a single objective function. MOVRPFlexTW aims to minimize the routing costs, which consist of travel costs and the fixed costs of the vehicles used for the service, while also maximizing the customer satisfaction.

AC

CE

To find some trade-off solutions between the routing costs and customer satisfaction, in this paper we propose a new Pareto-based multi-objective approach that combines the ant colony optimization (ACO) algorithm and mutation operators to solve MOVRPFlexTW. The motivation of using hybrid ACO (HACO) in our work is multifold. Firstly, the ACO algorithm is a swarm based meta-heuristic approach for solving combinatorial optimization problems which has produced good solutions within a reasonable computation time (e.g., multidimensional knapsack problem [39], traveling salesman problem [13] and job shop scheduling problem), especially many ACO-based approaches have been proposed to solve different variants of the VRP in the literatures (e.g., vehicle routing problem with simultaneous pickup and delivery [28], heterogeneous vehicle routing problem with mixed backhaul [49] and multi-compartment vehicle routing problem [38]). Secondly, swarm optimization algorithms maintain a pool of solutions during the search process. Therefore, swarm optimization algorithms are particularly suited for finding many trade-off solutions in a single run [31]. Thirdly, a metaheuristic is a general algorithmic framework for finding solutions to optimization problems within which local heuristics are guided and adapted to effectively explore the solution space [25]. Finally, different meta-heuristics have already been applied to other variants of VRP and provided good results [37, 46, 32, 7, 26]. All these factors have strongly motivated us to employ hybrid ACO algorithm to

2

ACCEPTED MANUSCRIPT

solve MOVRPFlexTW for the first time.

CR IP T

The main contributions of this paper are as follows. Firstly, we have introduced and modeled the MOVRPFlexTW. Secondly, to produce high-quality solutions, we have developed an ACO metaheuristic based algorithm to solve the MOVRPFlexTW problem. We have not only hybridized the ACO meta-heuristic with mutation operators which exploit the neighborhood search space, but also incorporated the concept of Pareto optimality for multi-objective optimization into this hybrid algorithm. This hybridization is an effective way to find trade-off solutions for the MOVRPFlexTW. Finally, we have conducted extensive experiments using publicly available benchmark instances to assess the operational gains of using multiple objectives and flexible time window. We have also analyzed the performance of our algorithm experimentally. As reported later, our approach can effectively generate trade-off solutions to MOVRPFlexTW.

2

AN US

The remainder of this paper is organized as follows. Section 2 contains a brief overview of the relevant literatures. Section 3 formally describes MOVRPFlexTW. Section 4 presents our Pareto-based hybrid ant colony optimization algorithm for solving MOVRPFlexTW. Section 5 describes our computational experiments that investigated the performance of the proposed algorithm and the advantageous of MOVRPFlexTW, respectively. Section 6 applies the presented MOVRPFlexTW on a real case study. Section 7 presents our concluding remarks.

Literature review

2.1

M

The aim of this section is to provide a brief but comprehensive overview of previous work related to the vehicle routing problem with time windows using multi-objective optimization techniques.

The vehicle routing problem with time windows

AC

CE

PT

ED

There have been extensive investigations into VRPTW because of its high computational complexity and usefulness in real life, especially with respect to the extremely important economic impact on all logistic systems. Some exact methods have been proposed for the VRPTW, including dynamic programming, column generation, and Lagrangian relaxation based methods [18]. However, exact methods often perform poorly when solving intermediate and large VRPTW instances. Thus, heuristic and meta-heuristic methodologies, which seek approximate solutions in polynomial time instead of exact methods with intolerable computing cost, are very suitable for solving larger VRPTW instances. Many different heuristics and meta-heuristics for solving the VRPTW have been proposed, such as genetic algorithm [22], simulated annealing [5], ACO, particle swarm optimization, and tabu search. Along with the basic features of each solution method, experimental results for benchmark test problems have also been presented and analyzed in the literatures. Br¨aysy and Gendreau surveyed route construction, local search algorithms, and meta-heuristics in [10] and [11]. Interested readers can refer to the review papers [18, 3, 11] for details of mathematical formulations, exact methods and meta-heuristic algorithms for the VRPTW. VRPSTW and VRPFlexTW are relaxations of VRPTW (in which time windows are treated as hard constraints). The time windows are unboundedly relaxed in VRPSTW and relaxed according to a fixed tolerance in VRPFlexTW. That is, the desired time windows can be violated at the cost of customer satisfaction in VRPSTW and VRPFlexTW. In particular, customers are assumed to be available for the service at any moment of the planning horizon in VRPSTW. VRPFlexTW is distinct from VRPSTW in that the former considers a bounded time windows violation. Vehicles in VRPFlexTW are allowed to provide both early and late services with a limited penalty, and to wait at no cost in the case of early 3

ACCEPTED MANUSCRIPT

arrival (until the extended earliest time) without any waiting time limitations [43]. When compared with VRPTW, VRPSTW and VRPFlexTW have the largest and larger feasible solution spaces [43], respectively.

CR IP T

Compared with the VRPTW, there have been significantly fewer investigations into VRPSTW and VRPFlexTW. The majority of literatures on VRPSTW and VRPFlexTW consider a linear penalty function for the violations of time windows. Chiang and Russell [14] and Balakrishnan [2] developed different tabu searches for the VRPSTW problem. Taillard et al. [44] also proposed a tabu search for the VRPSTW, in which a linear penalization is only applied to late visits. Calvete et al.[12] applied goals programming to VRPSTW with a linear penalization for early and late visits. Ibaraki et al. [24] measured time window violations as a nonconvex, piecewise linear, and time dependent function. Zare-Reisabadi et al. [50] developed an ant colony algorithm with local searches for the VRPSTW with a linear penalization for early and late visits. Bhusiri et al. [8] solved VRPSTW using column generation, which penalized both early and late arrivals, but the arrival time was bounded in an outer time window. Salani et al. [41] solved VRPSTW using a branch-and-cut-and-price technique. Koskosidis et al. [30] proposed a heuristic algorithm to decompose VRPSTW into an assignment component and a series of routing and scheduling components.

Multi-objective VRPTW

ED

2.2

M

AN US

Figliozzi [19] proposed an iterative route construction and improvement algorithm for the vehicle routing problem, which considered bounded early and late visits. Although this problem is called VRPSTW in [19], in fact it is a VRPFlexTW. Tas¸ et al. [43] developed a solution procedure that constructs feasible vehicle routes for VRPFlexTW via a tabu search algorithm. Qureshi et al. [35, 36] developed a column generation based exact algorithm for the vehicle routing and scheduling problem with semi soft time windows (VRPSSTW). The VRPSSTW presented in [35, 36] only considered an upper bound on the late time window violations, so this problem can be viewed as a special case of the VRPFlexTW.

AC

CE

PT

Jozefowiez et al. [27] classified the different objectives of variants of the VRP according to the component of the problem with which they are associated, i.e. the route, the resources, and the node/arc activity. The most common objective related to the route is to minimize the cost of the generated solutions with respect to distance traveled and/or time required. Other authors aimed to reduce the imbalance among the workloads of the vehicles, where workloads can be expressed as the number of customers visited, the quantity of goods delivered [6], the time required, or the route length [6]. Most studies dealing with objectives related to node/arc activity replaced the time windows with an objective that minimizes either the number of violated time windows constraints [19] or the total wait time of the customers and/or drivers due to earliness or lateness, or that maximizes customer satisfaction. The main objectives related to resources can be roughly classified into vehicle-related objectives and goods-related objectives. The majority of vehicle-related objectives focus on minimizing the number of vehicles or maximizing vehicle cost-effectiveness in terms of time or capacity. Goods-related objectives typically consider the nature of the goods. For instance, if the merchandise is perishable, we should minimize the deterioration. Many techniques have recently been proposed for solving multi-objective VRPTW. These strategies can be divided into three general categories: hierarchical techniques, aggregating approaches, and Pareto methods [6, 5]. Hierarchical techniques, which establish a hierarchy among the objectives to be optimized at the same time, intrinsically consider one objective as more important than others [5]. Aggregating approaches include all of objectives in a single objective function using mathematical transformations. As with hierarchical techniques, aggregating approaches are also unable to find all the Pareto optimal solutions. Additionally, the main disadvantage of aggregating approaches is that 4

ACCEPTED MANUSCRIPT

the weights are difficult to determine precisely [6, 5]. Pareto methods apply the notation of Pareto dominance to evaluate the quality of a solution or to compare solutions [27], so they can overcome the drawbacks of two other methods and are becoming increasingly popular. Pareto methods are mainly used in evolutionary and swarm algorithms. Evolutionary and swarm algorithms create a pool of solutions (called a population) during the search process. Therefore, these algorithms are particularly suited for obtaining a set of trade-off solutions for multi-objective problems [31]. Many pioneers have proposed Pareto-based hybrid methods for multi-objective VRPTW by combining evolutionary/swarm algorithms and local searches, heuristics, and/or exact methods [6, 5, 22].

3

AN US

CR IP T

There are few studies on multi-objective VRPFlexTW. Therefore, a major research effort is necessary in this field, including hybrid methods that take advantage of different heuristic evolutionary/swarm algorithms, exact methods, and local searches to improve the quality of the solutions and/or to reduce the computing time. In this paper, we aim to solve the multi-objective vehicle routing problem with flexible time windows using a Pareto-based hybrid ACO algorithm, which minimizes the total routing costs and maximizes the average customer satisfaction. This formulation allows us to obtain solutions with different customer satisfaction levels and costs incurred by the traveled distance and the number of vehicles used, so that the decision maker can decide which solution is the most suitable according to certain criteria.

A multi-objective optimization model for the VRPFlexTW

ED

M

A logistics company will typically make appointments with customers using order information or by telephone before distributing the goods or service. However, this often leads to distribution vehicles arriving earlier or later than the appointment time because of, for example, an unreasonable distribution scheduling or traffic jams. In real life, customers have a certain tolerance to services arriving earlier and later than the arranged time. This does not lead to significant losses to the logistics company in the short term, but an early or postponed service can decrease customer satisfaction and have a negative effect on the income and reputation of the company in the long term. Therefore, a logistics company should aim to improve customer satisfaction while minimizing the total traveled distance and number of vehicles.

PT

In this section, we consider a multi-objective VRPFlexTW formulation that aims to simultaneously maximize customer satisfaction while the vehicle routes are subject to capacity and time windows constraints and minimize the costs related to the traveled distance and number of vehicles.

AC

CE

The VRPTW can be modeled as a directed complete graph G(V, E) [18], with a vertex set V = {0, 1, · · · , n} and an edge set E = {hi, ji|i, j ∈ V and i 6= j}. Here, the central depot is denoted as 0, and customers are represented as i(i ∈ V and i 6= 0). M = {1, 2, · · · , m} is a set of vehicles. Each arc in the network describes the connection between two vertices. Each vehicle with limited capacity bk (k ∈ M ) and maximum route time rk starts from and returns back to the depot after serving a number of customers. Each customer i ∈ V (i 6= 0) has a variable demand di and must be visited only once by exactly one of the vehicles. The sum of all the demands serviced by vehicle k must be less than or equal to its capacity bk . Similarly, the total travel time of vehicle k cannot exceed its allowed route time rk . Additionally, the satisfaction µi (zi ) (zi is the time when vehicle arrives at customer i) is equal to 1 if customer i(i ∈ V ) is serviced within the appointed time interval [ei , li ]. A customer i cannot be served later than the extended latest time Li . A vehicle arriving at customer i earlier than the extended earliest time Ei incurs waiting time wi . Although the waiting time is permitted at no cost, customer satisfaction µi (zi ) decreases to 0 when wi > 0. (ei − Ei ) and (Li − li ) are customer i’s tolerances for being served earlier and later than the appointed time, respectively. Figure 1 shows the

5

ACCEPTED MANUSCRIPT

𝜇

1 complete satisfaction

𝑙𝑖

𝑒𝑖

𝐿𝑖

𝑧

CR IP T

𝐸𝑖

Figure 1: The customer satisfaction function.

customer satisfaction function. The satisfaction of customers serviced within the intervals [Ei , ei ] and [li , Li ] is assumed to vary linearly with the arrival time zi in this paper, and are calculated using zi < Ei , Ei ≤ zi < ei , ei ≤ zi ≤ li , li < zi ≤ Li , zi > Li .

AN US

  0,    z −E   eii −Eii , 1, µi (zi ) =  Li −zi    Li −li ,   0,

(1)

M

For the depot’s time window we assume that E0 and e0 are equal to 0, and l0 is equal to the extended latest time, L0 . This assumption implies that all vehicles start from the depot at time 0 and should return back to it before L0 . For simplicity, we set L0 = rk for any k ∈ M , that is, the maximum route time rk is also the depot’s time window for vehicle k.

ED

Before presenting the mathematical formulation, we define the following notations: hk is the transportation cost per unit distance of vehicle k,

PT

fk is the fixed cost incurred for using vehicle k, cij is the distance between vertex i and vertex j,

CE

si is the service time at vertex i,

wi is the waiting time at vertex i,

AC

tij is the time required for traveling from vertex i to vertex j. Decision variables xijk

( 1, if vehicle k travels from vertex i to vertex j = 0, otherwise

yik

( 1, = 0,

if vertex i is served by vehicle k otherwise

Since the VRP was proposed by Dantzig and Ramser as a generalization of the traveling salesman problem, cost has mostly been associated with the traveled distance, but there are several other types 6

ACCEPTED MANUSCRIPT

of cost involved, such as waiting time, service costs(early and late arrivals)[42], the number of vehicles and delivery time [27]. For simplicity, in this paper we only consider the main routing costs resulted by the traveled distance and the number of vehicles. The traveled distance affects time and fuel costs, and the number of vehicles affects vehicle and labor costs [22]. The objective of the VRPFlexTW is to serve all the customers such that the following objectives are met and the following constraints are satisfied. Objectives:

CR IP T

• minimize the total routing costs; and • maximize the average customer satisfaction. Constraints:

• for every vehicle route, the total demand cannot exceed the capacity of the vehicle;

AN US

• the total travel time of each vehicle cannot exceed its maximum route time; • each customer is served exactly once by exactly one vehicle;

• each vehicle route starts from the depot and ends at the depot; and • the vehicle serves every customer within the required time frame.

M

Given the above defined objectives, parameters and decision variables, the problem can be formulated as follows: n

i=1

min

m X k=1

hk

n X n X

cij xijk +

i=0 j=0

m X

fk

(2) n X

s.t.

CE

i=1 m X

di yik ≤ bk

∀k ∈ {1, 2, · · · , m}

(4)

yik = 1

∀i ∈ {1, 2, · · · , n}

(5)

∀k ∈ {1, 2, · · · , m}

(6)

AC

k=1

n X j=1

n X

i=0 n X

(3)

x0jk

j=1

k=1

PT

n X

ED

max

1X µi (zi ) n

x0jk −

n X

xi0k = 0

i=1

xijk = yjk

∀k ∈ {1, 2, · · · , m}, ∀j ∈ {1, 2, · · · , n}

(7)

xijk = yik

∀k ∈ {1, 2, · · · , m}, ∀i ∈ {1, 2, · · · , n}

(8)

j=0

n X n X i=0 j=0

xijk (tij + si + wi ) ≤ rk

∀k ∈ {1, 2, · · · , m}

(9) (10)

w0 = s 0 = 0

7

ACCEPTED MANUSCRIPT

m X n X

xijk (zi + wi + si + tij ) = zj

k=1 i=0

Ei ≤ zi + wi ≤ Li

wi = max{0, Ei − zi }

xijk ∈ {0, 1}

yik ∈ {0, 1}

∀j ∈ {1, 2, · · · , n}

(11)

∀i ∈ {1, 2, · · · , n}

(12)

∀i ∈ {1, 2, · · · , n}

∀i, j ∈ {1, 2, · · · , n}, ∀k ∈ {1, 2, · · · , m} ∀i ∈ {1, 2, · · · , n}, ∀k ∈ {1, 2, · · · , m}

zi ≥ 0

∀i ∈ {1, 2, · · · , n}

(13) (14) (15) (16)

4

AN US

CR IP T

In this model, Objective (2) is to maximize the customer satisfaction. Objective (3) is to minimize the total routing costs, which consist of travel costs and fixed vehicle costs. Constraint (4) guarantees that the vehicle capacity is not exceeded; Constraint (5) ensures that each customer is served by exactly one vehicle; and Constraint (6) ensures that each route starts and ends at the depot. Constraints (7)-(8) guarantee that each customer is served exactly once. Constraint (9) ensures that the maximum route time is not exceeded; Constraint (10) defines the waiting and service time at the depot; Constraint (11) represents the relationship between the arrival time at a vertex and the departure time from its predecessor; Constraint (12) ensures that customers are served within the required time; and Constraint (13) defines the waiting time.

Multi-objective hybrid ant colony optimization algorithm

PT

ED

M

Hybridization enhances strengths and compensates weaknesses of two or more methods, with the aim of generating better solutions by combining the key elements of competing methodologies [6]. Most hybridized methods in literatures were designed by embedding one solution method into another one as a local search technique(e.g., hill-climbing method is incorporated into genetic algorithm to search for better routing solutions [22], the saving algorithm and λ-interchange mechanism are embedded into ACO to obtain improved solution method to VRPTW [15], simulated annealing is combined with evolutionary algorithm as an acceptance criterion to generate high quality non-dominated solutions for multi-objective VRPTW [6]). Some previous studies have shown the advantages of embedding local search techniques into evolutionary approaches for solving VRP [1].

AC

CE

ACO algorithms imitate the behavior of a real ant colony when searching for food. During the search, each ant marks the trail it is moving along by laying a substance called pheromone. Ants on a short path will be quicker to return to the nest, so more pheromone will be deposited on the shorter paths. The amount of pheromone in a path lets other ants know if it is promising. An important benefit of this algorithm is parallelism: several solutions can be built simultaneously, they exchange information during the procedure, and use information from previous iterations. With the aim of obtaining high quality, non-dominated solutions for the MOVRPFlexTW presented in Section 3, we propose a new approach to improve the efficiency of a Pareto-based multi-objective ACO algorithm using mutation operators. We can take advantage of the ACO algorithm and mutation operators to obtain a search procedure with a well-balanced exploration-exploitation behavior using hybridization. This proposed HACO starts with P ants at the depot. Each ant constructs a sloution to the problem according to the mechanism of solution construction. Mutation operators will then be used to the randomly selected solutions in order to exploit their neighborhood search space, and K new solutions will be obtained. An acceptance criterion will then be used to select P prospective solutions from the P + K solutions based on their fitness computed by evaluation function. The next step will be to use a Pareto-optimization method to select non-dominated solutions from the P accepted solutions. Whereafter, the pheromone trails will be updated to explore the search space of P 8

ACCEPTED MANUSCRIPT

ants and produce better subsequent solutions. The HACO will then reiterate through this process until a predefined number of iterations has been reached. The specific features and general process of the proposed HACO are described below.

4.1

Solution construction

CR IP T

When using HACO to solve the MOVRPFlexTW, each ant starts at the depot and constructs a route by incrementally selecting customers under some constraints. When the next selection in the current route is impossible, the ant returns to the depot and constructs another route by assigning the unrouted customers. This process is repeated until all customers are routed, and a feasible solution can be attained. At each construction step, each ant located at customer i uses a probabilistic formula to select the next customer j. That is, ( arg max{(τiu )α (ηiu )β (ψiu )γ } for u ∈ / tabul , if q ≤ q0 , j= (17) S, otherwise,

AN US

where S is obtained by roulette wheel rule with probability distribution   P (τij )α (ηij )β (ψij )γ , if j ∈ / tabul , (τiu )α (ηiu )β (ψiu )γ u∈tabu / l pij = 0, otherwise.

(18)

In (17) and (18), τiu is equal to the amount of pheromone on the path between the current customer i and possible customer u; ηiu denotes the visibility on edge (i, u), which is defined as the inverse of the distance between customer i and customer u, ηiu = c1iu ; ψiu = max(1,l1 −z l ) denotes the urgency u

iu

ED

M

l is the time when ant l arrives at customer location u from customer i; and of edge (i, u), where ziu α, β, and γ are parameters that determine the relative influences of the pheromone trail, the visibility value, and the urgency value, respectively. Customers already visited by ant l are stored in the set tabul and are not considered for selection. q is a random variable uniformly distributed in [0, 1] and q0 , 0 ≤ q0 ≤ 1, is a parameter that controls the exploitation against the exploration of ants during the searching process.

PT

The solutions can be constructed by using the following procedure. Algorithm 1 ( solution construction)

CE

1. Information:

• Input: 0 ≤ q0 ≤ 1, the number of ants P , the demand of each customer di for all i ∈ {1, 2, · · · , n}, the limited capacity bk and maximum route time rk for all k ∈ {1, 2, · · · , m}.

AC

• Output: P initial solutions.

2. Initialization: • For l = 1, 2, · · · , P , start ant l at the depot, and initialize the set of visited customers tabul = ∅ and the set of candidate customers Candl = {1, 2, · · · , n}. • Set k = 1 and l = 1.

3. For the given vehicle k (k ∈ {1, 2, · · · , m}), set its available capacity bbk = bk and route time rrk = rk .

4. Update the candidate list: According to tabul , bbk , rrk , the vehicle capacity constraints (4) and route time constraints (9), update Candl for ant l. 5. Criterion of returning to the depot: If Candl = ∅, then set k = k + 1 and go to Step 3(the ant l returns to the depot and constructs another route); otherwise, go to Step 6.

9

ACCEPTED MANUSCRIPT 0 1

2

3

0

4

5

6

7

0

8

9

10

0

9

10

11

12

13

14

0

05 11

12

13

14

3 1

0 1

314 0

2

2

4

5

6

7

0

8

4

0

3 0

1

4

14

7

11 10

8

0 12

13 1

9

2

4Figure 5 2: 6 A7typical 8 9vehicle 10 routing 11 problem. 12 13

2

3

0

4

5

6

7

14

8

AN US

10

0 1

6

7

11

3

5

CR IP T

12

13

6

2

0

8

9

10

0

11

12

13

14

0

3 from figure 2. (a) Sequence codes of the3solution 9

1

2

3

1

1

4

6 27

5

2

8

9

10

11

12

13

14

4 4 from figure 3(a). 5 14(b) Sequence codes of the customers

5

M

14

0

6

Figure 3: Sequence codes.3

8. 9.

4.2

ED

PT

7.

CE

6.

0 6 12 1 13 2 Select the13 next customer: If12Candl 6= ∅,11then select the next customer to be visited j(j ∈ / tabul and 7 4 j ∈ Candl ) according to Equations (17) and (18). 7 5 11 Set tabul = tabul ∪ {j}, rrk = rrk − tij 10 − sj − wj (i is the vertex visited by vehicle k before j) and 8 bbk = bbk − dj , and go to Step 8. 1410 08 6 If tabul = {1, 2, · · · , n}, then set l = l + 1 and go to Step 9; otherwise, go to Step 4. 12otherwise, set k9 = 1 and go to Step 3. Stopping criterion: If13 l > P , then stop; 9

11

1 2Operation 3 4 5 6 7 8 9 10 11 Mutation 10 1

2

3

4

5

6

7

8

9

10

7

12 11

13

14 138 14

12

AC

In this paper, we introduce the mutation operation into the ACO algorithm to diversify the population and prevent convergence to a local optimum. The aim of the mutation operation is to alter some 9 customers and produce a new solution. There are 3three mutation operators in this algorithm: the 0 1 2 3 0 9 8 7 0 6 5 4 0 10 12 operators 0 13 14randomly 0 interchange operator, shifting operator, and inverse operator. These 11 mutation 1 2 conduct customer exchanges without considering the constraints on vehicle capacity and route time. Let q1 , q2 , and q3 be such that 0 ≤ q1 , q2 , q3 ≤ 1 and q2 < q3 . Then, two random numbers qmu1 and 4 qmu2 are uniformly generated in [0, 1]. The solution constructed by ant l is not selected for mutation 1 2 3 4 5 6 7 8 9 10 11 125 13 14 if qmu1 > q1 , otherwise the interchange operator is applied if qmu2 ≤ q2 , the shifting operator is 14inverse operator is applied if qmu2 > q3 . applied if q2 < qmu2 ≤ q3 , or 0 6 To describe the mutation operators, we present an example solution in figure 3(a), which is a sequence 0 1 2 312 0 in 9 figure 8 7 2.0 The 6 mutation 5 4 operators 0 10 were 11 implemented 12 0 13 on 14 the0 codes transformed 13 from the solution 11

7

10 10 8

9

7

112

4 ACCEPTED 10 MANUSCRIPT 0 1

2

3

0 4 14

5

6

7

0

8

9

8

10

5 11

0

12

0 93 113 2

3

4

12 1 6 7

5

92 10

8

12

13

4

7

9 12

5

10

6

8

0

14

7

4

5

8

0 14 2 3 13

14

6

11

11 Exchanging 10

14

13

6 11

12

13

1

0

1

2

14 2 3 4 0 7 10 9 5 0 3 4 5 6 7 8 9 10

CR IP T

9 Inserting 0 7 Figure 4: Sequence acquired 11by using interchange operator. 10 6 8 0 11 11 12 8 13 14

Shifting 9 1

2 3 4 5 6 2 3 9 4 5

7 6

13

1

0

8 9 10 11 12 13 14 7 8 10 11 12 13 14

AN US

1

12

Figure 5: Sequence acquired by using shift operator.

3

sequence codes shown 14 5were 0 1 in 2figure 3 3(b), 0 9 which 0 2obtained 6 7 by8 eliminating 10 0 0s 11 (the12depot) 13 from 14 the0 solution shown in figure 3(a).

M

4 e.g., 1 and 14, 5 and 7, 6 and 9, Interchange operator. Randomly select several pairs of customers, 5 8 and 10. By exchanging these pairs of customers we obtain a new sequence. The result is shown in 14 figure 4. 0

6

ED

Shift operator. Randomly select two points, e.g., 4 and 9. The last customer between these two points (including these two points) is moved12 to the first position, and the other customers are separately 13 moved to the position next to them. The new sequence obtained using the shifting operator is shown 7 11 in figure 5.

PT

Inverse operator. Randomly select two points, e.g., 10 4 and 9. A new sequence is obtained by inversing 8 6. the subsequence between these two points. The result is shown in figure

AC

CE

Three new solutions can be obtained by inserting 0s into the new sequence acquired using the muta9 route satisfies the constraints on vehicle tion operation. In the new solutions shown in figure 7, each capacity and route time.

1

2

3

4

5

6

7

8

9

10

11

12

13

14

5

4

10

11

12

13

14

Inversing

1

2

3

9

8

7

6

Figure 6: Sequence acquired by using inverse operator.

0 1

2

3

0

9

8

7

0

11

6

5

4

0

10

11

12

0

13

14

0

9

1

2

14

3

4

2 3

5

6

7

8

9

10

11

12

13

14

ACCEPTED MANUSCRIPT Shifting

4

7

9

5

10

6

8

11

12

1 2 3 4 5 6 7 8 9 10 11 12 13 1 2 3 9 Inversing 4 5 6 7 8 10 11 12 0 14 2 3 4 0 7 9 5 0 10 6 8

13

1

14 13 14 0 11 12

13

1

0

13

14

0

(a) New solution obtained by using interchange operator.

1 2 0 1

3 2

9 3

8 0

7 6 5 4 10 9 4 5 0 6

11 7

12 8

13 10

14 0 11

12

(b) New solution obtained by using shift operator.

Inserting 0

0 1

2

3

0

9

8

7

0

6

5

4

0

10

11

12

0

13

14

0

Figure 7: New solutions.

4.3

Acceptance criterion

CR IP T

(c) New solution obtained by using inverse operator.

Ratiol =

AN US

Suppose that K offspring solutions are obtained by applying the mutation operators, then the number of solutions is P + K (where P is the number of ants and K ≤ P ). The P + K solutions are ranked according to their fitness Ratiol (1 ≤ l ≤ P + K) which is evaluated using Satisl , Costl

for 1 ≤ l ≤ P + K

(19)

M

where Satisl and Costl are the objective values related to solution l, which are calculated using objectives (2) and (3), respectively. Then, the first P solutions are used to update the pheromone trails described in Subsection 4.5.

Pareto-optimization method

ED

4.4

CE

PT

For a problem with more than one objective function, given any two solutions (S1 and S2 ) there are two possibilities: one dominates the other or they are indifferent. We say that solution S1 dominates (or is preferred to) S2 (S1  S2 , where ≺ denotes worse and  denotes better) if and only if S1 is better than S2 in at least one objective, and not worse in the others. Two solutions are called indifferent if neither dominates the other.

AC

We consider the two objectives of the MOVRPFlexTW presented in Section 3 to be equally important, so we can extend the above concepts to find the set of non-dominated solutions from the P accepted solutions.

4.5

Pheromone Updating

Updating the pheromone trails is a key element to adaptive learning technique of ACO and to improve the subsequent solutions. In this paper, after each iteration the pheromone trails of traveled edge (i, j) are updated using P X new old τij = (1 − ρ) × τij + ∆τijl (20) l=1

12

ACCEPTED MANUSCRIPT

where τijnew is the pheromone on edge (i, j) after updating; τijold is the pheromone on edge (i.j) before updating; ρ is a parameter that controls the speed of pheromone evaporation; P is the number of ants; and ∆τijl is the amount of pheromone that ant l deposits on edge (i, j). ∆τijl is defined as ∆τijl =

(

Q×Satisl Costl ,

if edge (i, j) belongs to T l

0,

otherwise

(21)

CR IP T

where Q is a constant related to the amount of pheromone deposited by the ants; and T l is the tour built by ant l. Equation (21) implies that edges belonging to better tours receive more pheromone, so they are more likely to be chosen by ants in subsequent iterations. Too many or too few pheromone trails may lead to premature convergence. Therefore, we set τmin and τmax (the minimum and maximum pheromone trails) for all pheromone trails τij to satisfy τmin ≤ τ ≤ τmax , following [15].

4.6

Hybrid ant colony optimization algorithm for MOVRPFlexTW

AN US

Based on the procedures presented in Subsection 4.1-4.5, the steps of the hybrid ant colony optimization (HACO) algorithm for solving MOVRPFlexTW are shown as follows. Algorithm 2 (HACO) 1. Information:

• Input: the parameters in the formulation (2)-(16),such as m, n, hk and so on.

M

• Output: a set of non-dominated solutions for MOVRPFlexTW (2)-(16). 2. Initialization:

ED

• Initialize the parameters, such as the number of ants P , the maximum number of iterations maxiter, the pheromone trails τij (0 ≤ i, j ≤ n), the visibility ηij (0 ≤ i, j ≤ n), 0 ≤ ρ, q0 , q1 , q2 , q3 ≤ 1(q2 < q3 ) and so on.

PT

• Set the counter to nc = 1, the set of global non-dominated solutions Gnon−dom = ∅ and the set of local non-dominated solutions Lnon−dom = ∅, respectively. 3. Generate solution: If nc ≤ maxiter, produce P solutions by using the solution construction method shown as Algorithm 1; otherwise, the algorithm terminates and outputs Gnon−dom .

CE

4. Mutation:

(a) Set l = 1.

(b) Uniformly generate two random numbers qmu1 and qmu2 in [0, 1].

AC

(c) If qmu1 ≤ q1 , the mutation operation is performed on the solution generated by ant l   interchange operator if qmu2 ≤ q2 ; shifting operator if q2 < qmu2 ≤ q3 ;   inverse operator if qmu2 > q3 ; otherwise, go to Step 4.(d).

(d) Set l = l + 1 and repeat Step 4.(b) and Step 4.(c) until l > P . 5. Acceptance criterion: Arrange the P + K solutions (P solutions obtained by P ants and K new solutions obtained using the mutation operation) in the order of descending fitness computed using equation (19), and accept the first P solutions. 6. Compute non-dominated solutions:

13

ACCEPTED MANUSCRIPT

(a) Compute the routing costs and average customer satisfaction related to each solution. (b) Find a set of non-dominated solutions from the P accepted solutions, and use it to update Lnon−dom . (c) Based on the set Gnon−dom ∪ Lnon−dom , generate a new set of non-dominated solutions and use it to update Gnon−dom . 7. Update pheromone trails: According to equations (20)-(21), update the pheromone trails τij (0 ≤ i, j ≤ n), and τij is replaced by τmax if τij > τmax , or by τmin if τij < τmin . 8. Set nc = nc + 1 and go to Step 3.

Computational analysis

CR IP T

5

5.1

AN US

To assess the advantages of MOVRPFlexTW and the performance of the proposed HACO, we applied them to a set of Solomon’s VRPTW benchmark examples. The time windows of these instances are revised to adapt them to the MOVRPFlexTW formulation presented in Section 3. The algorithm was implemented in Matlab R2010a and executed on a computer with an Intel Core(TM) i7-4790 CPU 3.60GHz, with 16.0GB of RAM.

Description of experimental problems

ED

M

Solomon’s problems include 56 data sets of size n = 100 that are available from Solomon’s web site 1 . These instances are categorized into C1, C2, R1, R2, RC1 and RC2. Customers of the instances in sets C1 and C2 are clustered either geographically or according to time windows. Customers of the instances in sets R1 and R2 are randomly located, whereas those in sets RC1 and RC2 have mixed characteristics of random locations and clusters. Moreover, the instances vary in vehicle capacity, fleet size, and maximum route times. Specially, the instances in C1, R1 and RC1 have narrow depot time windows so only a few customers can be served by the same vehicle. In contrast, instances in C2, R2 and RC2 have wider time windows so many customers can be serviced by the same vehicle. Like reference [43], we set the cost coefficients (hk , fk ) to (2.0, 400).

AC

CE

PT

The customers’ details in Solomon’s instances are given in sequence of customer index, X-coordinates, Y-coordinates, demand for load di ( i ∈ N ), ready time ei ( i ∈ N ), due time li ( i ∈ N ), and service time si ( i ∈ N ). All the customers in Solomon’s instances do not have any tolerance to unpunctual service, that is, all customers must be serviced between the ready time and the due time. To obtain instances that are relevant to the MOVRPFlexTW model in (2)-(15), we generated a flexible time window [Ei , Li ] for each customer i ∈ N with respect to the length of the original time window, where Ei = ei − wbli (li − ei ) and Li = li + wbui (li − ei ) (wbli and wbui are the fractions used to set the maximum allowed violations).

5.2 Sensitivity analysis and parameter settings A crucial feature of ACO is that the quality of the solution can be controlled by the parameters. In the field of evolutionary computing, two traditional approaches are usually employed to choose parameter values [21], following the scheme offered in [17]: parameter tuning, where (good) parameter values are established before the run of a given evolutionary algorithm. In this case, parameter values are fixed in the initialization stage and do not change while the evolutionary algorithm is running; parameter control, where (good) parameter values are established during the run of a given evolutionary 1

http://w.cba.neu.edu/ msolomon/problems.htm

14

ACCEPTED MANUSCRIPT

algorithm. In this case, parameter values are given an initial value when starting the evolutionary algorithm and they undergo changes while the evolutionary algorithm is running. In this paper the first approach is used. In order to illustrate how those parameters affect the overall performance of the HACO, the analyses between the approximate optimal solution and the crucial parameters, including α, β, γ, etc., are carried out on three randomly selected problems C101, R101, and RC101. When investigating the relationship between one of the parameters and the approximate optimal solution, the other parameters are kept unchanged. This investigation is executed by using the proposed HACO to solve the MOVRPFlexTW where wbli = wbui = 0.

CR IP T

(1) Parameter α Table 1. The approximate optimal objective values obtained by using different α. α

Approximate optimal objective values

1 2 3 4 5 6

C101 5657.88 5712.94 5711.35 5712.94 5729.30 5778.75

RC101 9273.84 9276.26 9352.64 9386.85 9400.46 9397.22

AN US

R101 10901.60 10924.61 11092.94 11220.32 11154.76 11228.54

M

If α = 0, the most urgent and closest customer is likely to be selected. In this case, HACO is same as the classical greedy random search algorithm. Table 1 presents the approximate optimal objective values obtained by using different parameters α. According to the results shown in Table 1, we can know that the larger α is, the more easily HACO gets a local optimal solution. During the experiments with HACO, it has been detected that the relatively good values of the parameter α are within the interval [1,3]. As shown in Table 1, for the problems C101, R101, and RC101, the best α obtained in our experiments was found to be equal to 1. (2) Parameter β β

Approximate optimal objective values

1 2 3 4 5 6

C101 5826.77 5705.90 5657.88 5657.88 5705.90 5775.90

R101 11854.14 11271.08 10901.60 11275.20 11093.09 11477.52

RC101 9863.11 9626.78 9359.63 9273.84 9274.47 9634.22

CE

PT

ED

Table 2. The approximate optimal objective values obtained by using different β.

AC

If β = 0, the next serviced customer is selected without using the heuristic factors η. In this case, HACO tends to converge to local optimums and does not perform well. Table 2 presents the approximate optimal objective values obtained by using different parameters β. For problem C101, the optimal objective value 5657.88 is obtained when β equals to 3 and 4. For problems R101 and RC101, their optimal objective values are obtained when β equals to 3 and 4, respectively. According to the numbers shown in Table 2, we can know that the relatively good values of the parameter β are within the interval [3,5]. (3) Parameter γ

If γ = 0, the next serviced customer is selected without using the heuristic factors ψ. In this case, HACO also tends to converge to local optimums. Table 3 presents the approximate optimal objective values obtained by using different parameters γ. For problem C101, the optimal objective value 15

ACCEPTED MANUSCRIPT Table 3. The approximate optimal objective values obtained by using different γ. γ

Approximate optimal objective values

1 2 3 4 5 6

C101 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88

R101 10901.60 11324.38 11668.15 11708.59 11834.21 11873.86

RC101 9299.68 9273.84 9384.50 9587.96 9600.13 9646.13

CR IP T

5657.88 is always obtained whenever the value of parameter γ changes in the prespecified interval [1,6]. For problem R101, the approximate optimal objective value gets worse with the increase of γ. In other words, the larger γ is, the more easily HACO gets a local optimal solution for problem R101. For problem RC101, the approximate optimal objective values are relatively good when parameter γ is within the interval [1,3], especially the optimal objective value 9273.84 is obtained when γ equals to 2. According to the numbers shown in Table 3, we can conclude that the relatively good values of the parameter γ are within the interval [1,3]. (4) Parameter ρ

C101 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88

R101 11160.75 11091.53 11110.01 11187.18 11009.76 11004.93 10901.60 10993.05 11203.36

RC101 9640.95 9540.96 9568.33 9488.56 9273.84 9366.63 9273.84 9336.70 9702.68

ED

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Approximate optimal objective values

M

ρ

AN US

Table 4. The approximate optimal objective values obtained by using different ρ.

AC

CE

PT

A reasonable amount of pheromone on the ants’ trails is dependent on the suitable value of parameter ρ. A smaller ρ increases the difficulty of obtaining the optimal solution, because the pheromone laid on better trails cannot distinguish the better routes. However, a larger ρ means that HACO is more likely to converge to a local optimum, because the amount of pheromone on the currently worst trails is dramatically less. Table 4 presents the approximate optimal objective values obtained by using different ρ. For problem C101, the optimal objective value 5657.88 is always obtained whenever the value of parameter ρ changes in the prespecified interval [0.1,0.9]. For problem R101, the approximate optimal objective values are relatively good when parameter ρ is within the interval [0.7,0.8]. For problem RC101, the optimal objective value 9273.84 is obtained when ρ equals to 0.5 and 0.7. From the numbers presented in Table 4, we can conclude that the HACO performs relatively well when ρ is in the interval [0.5, 0.8]. (5) Parameter q0 The edges can be classified into three families: i) edges belonging to the last best tour (BE, best edges), ii) edges that do not belong to the last best tour, but which did in one of the two preceding iterations (TE, testable edges), and iii) the remaining edges, that is, the ones that have never belonged to a best tour or have not in the last two iterations (UE, uninteresting edges) [16]. If q0 is close to 1, the HACO favors the BE and TE edges. Consequently, ants easily converge to a local optimum. If q0 is close to 0, HACO cannot effectively use the information from previous iterations. Consequently, the ants tend to routes with worse objective values. Table 5 presents the approximate optimal objective values 16

ACCEPTED MANUSCRIPT Table 5. The approximate optimal objective values obtained by using different q0 . q0

Approximate optimal objective values

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

C101 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88 5705.90 5705.90

RC101 9669.02 9661.53 9472.84 9450.39 9430.35 9273.84 9298.96 9281.10 9537.71

CR IP T

R101 11348.45 11405.97 11410.33 11043.63 11040.30 10901.60 10941.43 10975.51 11153.58

(6) Mutation rate q1

AN US

obtained by using different q0 . For problem C101, the optimal objective value 5657.88 is always obtained whenever the value of parameter q0 changes in interval [0.1,0.7]. For problems R101 and RC101, the approximate optimal objective values are relatively good when parameter q0 is within the interval [0.6,0.8], especially their optimal objective values are obtained when q0 equals to 0.6. From the numbers presented in Table 5, we can conclude that the HACO performs relatively well when q0 is in the interval [0.6, 0.8].

Table 6. The approximate optimal objective values obtained by using different q1 . q1

Approximate optimal objective values

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

C101 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88 5657.88

ED

M

R101 11405.97 11405.97 11410.33 11043.63 10975.51 10901.60 11040.30 11410.33 11405.97

RC101 9472.84 9661.53 9472.84 9450.39 9281.10 9273.84 9298.96 9430.35 9537.71

AC

CE

PT

One of the crucial things in the mutation operation is determining the proper mutation rate. Very small mutation rate might not prevent falling back into the visited local optima; however, large mutation rate might be very similar to a crude random multistart. Table 6 presents the approximate optimal objective values obtained by using different mutation rate q1 . For problem C101, the optimal objective value 5657.88 is always obtained whenever the value of parameter ρ changes in the prespecified interval [0.1,0.9]. For problems R101 and RC101, the approximate optimal objective values are relatively good when parameter q1 is within the interval [0.5,0.6], especially their optimal objective values are obtained when q1 equals to 0.6. From the numbers presented in Table 6, we can conclude that the HACO performs better when q1 is in the interval [0.5, 0.6].

Based on the above discussion and previous research, we implemented HACO with population size = 35, number of generation = 100, repetition = 10, mutation rate q1 = 0.6, α = 1, β = 4, γ = 2, q0 = 0.6, ρ = 0.7 and repetitions of interchanges = substring length for inversion = substring length for shifting = n/10, where n is the number of customers.

5.3

Analysis of results

We analyzed the experimental results from three different perspectives. Firstly, we compared the results obtained by HACO with previous results. Secondly, we investigated whether similar or better 17

ACCEPTED MANUSCRIPT

Table 7. Comparison of the original VRPTW solutions with the solutions obtained by implementing MOVRPFlexTW with the suggestive HACO for the category of type 1 (C1, R1 and RC1), where wbli = wbui = 0 (i ∈ {1, 2, · · · , n}). MOVRPFlexTW+HACO Best known

Best results #Veh.

Dis.

C101 C102 C103 C104 C105 C106 C107 C108 C109

10 10 10 10 10 10 10 10 10

828.94 828.94 828.06 824.78 828.94 828.94 828.94 828.94 828.94

10 10 10 10 10 10 10 10 10

R101 R102 R103 R104 R105 R106 R107 R108 R109 R110 R111 R112

19 17 14 10 14 12 11 10 12 11 10 10

1650.80 1434.00 1237.05 974.24 1377.11 1252.03 1100.52 960.26 1169.85 1112.21 1096.72 976.99

RC101 RC102 RC103 RC104 RC105 RC106 RC107 RC108

15 13 11 10 16 11 11 10

Avg.

11.62

Average in 10 runs Cost

#Veh.

Dis.

828.94 828.94 828.06 824.78 828.94 828.94 828.94 828.94 828.94

Gapdis (%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

5657.88 5657.88 5656.12 5649.56 5657.88 5657.88 5657.88 5657.88 5657.88

10 10 10 10 10 10 10 10 10

828.94 831.42 835.46 838.53 828.94 828.94 833.78 838.15 828.94

Gapdis (%) 0.00 0.30 0.89 1.67 0.00 0.00 0.58 1.11 0.00

Time (Sec.) 351 368 385 360 954 916 1267 1223 1159

19 17 14 10 14 12 11 10 12 11 10 10

1669.34 1443.18 1238.56 974.24 1382.50 1252.03 1108.80 960.88 1222.70 1114.56 1096.73 976.99

1.12 0.64 0.12 0.00 0.39 0.00 0.75 0.06 4.52 0.21 0.00 0.00

10938.68 9686.36 8077.12 5948.48 8365.00 7304.06 6617.60 5921.76 7245.40 6629.12 6193.46 5953.98

19 17 14 10 14 12 11 10 12 11 10 10

1675.45 1447.20 1239.40 974.24 1385.67 1258.61 1110.54 961.23 1225.60 1116.79 1096.73 987.17

1.49 0.92 0.19 0.00 0.62 0.53 0.91 0.10 4.77 0.41 0.00 1.04

1277 1265 1222 1225 1205 1187 1179 1187 1164 1135 1119 1073

1636.92 1470.26 1261.67 1135.48 1590.25 1427.13 1230.48 1142.66

15 13 11 10 16 11 11 10

1645.67 1476.50 1262.01 1135.48 1589.40 1424.73 1230.48 1142.66

1127.31

11.62

1131.13

CR IP T

Dis.

0.53 0.42 0.03 0.00 -0.05 -0.17 0.00 0.00

9291.34 8153.00 6924.02 6270.96 9578.80 7249.46 6860.96 6285.32

15 13 11 10 16 11 11 10

1648.76 1479.34 1263.25 1135.50 1592.45 1427.82 1234.73 1148.96

0.72 0.62 0.13 0.00 0.14 0.05 0.35 0.55

1235 1182 1185 1120 1114 1086 1138 1127

0.30

6910.54

11.62

1134.57

0.62

1048.36

M

#Veh.

AN US

Instance

ED

results were obtained when multiple objectives were optimized simultaneously, rather than when only one objective was considered. Finally, we investigated the effects of implementing different time windows.

PT

5.3.1 Comparison with previous results

AC

CE

To assess the performance of HACO in terms of applicability and solution quality, we ran some experiments using Solomon’s original VRPTW instances. Tables 7 and 8 present the best known solutions reported in the literatures and the solutions obtained using the proposed method for MOVRPFlexTW with wbli = wbui = 0 (i ∈ {1, 2, · · · , n}). According to customer satisfaction definition presented in Section 3, every customer’s satisfaction is equal to 1 if wbli = wbui = 0 (i ∈ {1, 2, · · · , n}). In this case, MOVRPFlexTW becomes a single objective problem. In Table 7 and Table 8, the column labeled Best known gives the best known published solutions from [21, 22]; MOVRPFlexTW+HACO gives the best and average results obtained using MOVRPFlexTW with HACO over 10 runs; #Veh. and Dis. give the number of vehicles and total distance, respectively; Cost gives the routing costs; Gapdis (%) gives the difference between the distance obtained by HACO and the best known distance, by HACO−best known distance that is Gapdis = distance obtained × 100%, which can relatively measure best known distance the efficiency of the proposed HACO; Time(Sec.) gives the average computational time (in seconds) of the HACO. Table 7 and Table 8 also give the average values of every column over all type 1 instances and type 2 instances, respectively.

Bold numbers indicate that the solutions obtained by HACO are very similar to or better than the best known solutions from the literatures, considering either the number of vehicles or traveled distance. The main result is that solutions obtained using the proposed HACO method were very similar to or 18

ACCEPTED MANUSCRIPT

Table 8. Comparison of the original VRPTW solutions with the solutions obtained by implementing MOVRPFlexTW with the suggestive HACO for the category of type 2 (C2, R2 and RC2), where wbli = wbui = 0 (i ∈ {1, 2, · · · , n}). MOVRPFlexTW+HACO Best known

Best results Dis.

#Veh.

Dis.

C201 C202 C203 C204 C205 C206 C207 C208

3 3 3 3 3 3 3 3

591.56 591.56 591.17 590.60 588.16 588.49 588.29 588.32

3 3 3 3 3 3 3 3

R201 R202 R203 R204 R205 R206 R207 R208 R209 R210 R211

5 4 4 3 3 3 3 2 3 3 4

1206.42 1091.21 935.04 789.72 994.42 833.00 814.78 726.82 855.00 954.12 761.10

RC201 RC202 RC203 RC204 RC205 RC206 RC207 RC208

6 4 4 3 4 3 4 4

Avg.

3.44

Average in 10 runs Cost

#Veh.

Dis.

591.56 591.56 591.17 598.72 588.16 588.90 592.46 588.32

Gapdis (%) 0.00 0.00 0.00 1.37 0.00 0.07 0.71 0.00

2383.12 2383.12 2382.34 2397.44 2376.32 2377.80 2384.92 2376.64

3 3 3 3 3 3 3 3

591.56 594.28 602.13 615.68 596.77 593.92 595.66 590.12

4 4 4 3 3 3 3 2 3 3 4

1252.37 1091.21 937.56 794.56 994.43 906.14 825.48 726.50 855.00 939.37 761.50

3.81 0.00 0.27 0.61 0.00 8.78 1.31 -0.04 0.00 -1.55 0.05

4104.74 3782.42 3475.12 2789.12 3188.86 3012.28 2850.96 2253.00 2910.00 3078.74 3123.00

1134.91 1181.99 1026.61 798.46 1300.25 1153.93 1040.67 785.93

4 4 4 3 4 3 4 4

1423.70 1263.53 1048.00 799.52 1410.30 1146.32 1140.67 828.71

25.45 6.90 2.08 0.13 8.46 -0.66 9.61 5.44

4447.40 4127.06 3696.00 2799.04 4420.60 3492.64 3881.34 3257.42

855.65

3.33

884.29

2.70

3101.91

Gapdis (%) 0.00 0.46 1.85 4.25 1.46 0.92 1.25 0.31

Time (Sec.) 1260 1356 1401 1283 1359 1342 1377 1484

CR IP T

#Veh.

4 4 4 3 3 3 3 2 3 3 4

1257.98 1091.21 938.43 796.67 994.43 906.64 827.63 727.85 855.00 943.18 762.72

4.27 0.00 0.36 0.88 0.00 8.84 1.58 0.14 0.00 1.15 0.21

1791 1886 1818 1917 1850 1905 1925 1918 1897 1935 1964

4 4 4 3 4 3 4 4

1428.53 1267.95 1053.71 803.76 1413.49 1153.68 1143.56 835.03

25.87 7.27 2.64 0.66 8.71 0.02 9.89 6.25

1837 1813 1823 1826 1837 1810 1810 1809

3.33

888.21

3.31

1712.19

AN US

Instance

AC

CE

PT

ED

M

better than the best known solutions in 38 out of the 56 tested instances. When compared with the best known solutions, HACO achieves the minimum number of vehicles and a small average difference with respect to the distance (0.30% for the category of type 1 in Table 7 and 2.70% for the category of type 2 in Table 8). As can be seen, HACO leads to new best results with the smaller routing distance for instances RC105 and RC106 in Table 7, and R208, R210 and RC206 in Table 8. For the instances C101-109, HACO can always find the best known solution. For RC201, although the distance produced by HACO is worse than the best known distance, there are four vehicles, which is better than the best known result. For the other 42 instances, the Gapdis varies between 0.00% and 9.61%. Additionally, the average computational time of HACO for every problem in sets C1, R1, and RC1 is less than 1300 seconds. The problems in sets C2, R2 and RC2 required slightly more computational time because of their wider time windows, which allow a more flexible arrangement in the routing construction process. Besides, a vehicle with longer route also takes up more computational time during the feasibility evaluation processes [45]. In general, the results obtained by the proposed HACO method are good in terms of solution quality and computational time, when compared with the best published results. To study the consistency and reliability of the results obtained by HACO, 10 different but repeated simulations with randomly generated initial populations have been performed for the Solomon’s 56 instances. The average values of the various simulation results are also represented in Table 7 and Table 8. As shown in Table 7 and Table 8, the average number of vehicles are equal to the best number for each instance and except instance RC201, the average values of Gapdis obtained from HACO for the 10 different but repeated simulation runs are found to vary within 0.00-9.89% from the best values. Therefore, we can conclude that the results obtained by HACO are rather consistent. It is also observed that the category of type 1 (C101-C109, R101-R112 and RC101-RC108) gives smaller average values of Gapdis as compared with the category of type 2 (C201-C208, R201-R211 and RC201-RC208), since the number of customers per route is smaller for the category of type 1 [45]. 19

ACCEPTED MANUSCRIPT

Table 9. Comparisons between the HACO and state-of-the-art algorithms from literatures for the category of type 1 (C1, R1 and RC1). Instance

HMOEA

MOGA

MACS-IH

HSFLA

ILNSA

LGA

#Veh. 10 10 10 10 10 10 10 10 10

Dis. 828.94 828.94 828.06 825.65 828.94 828.94 828.94 828.94 828.94

#Veh. 10 10 10 10 10 10 10 10 10

Dis. 828.94 828.94 828.06 824.78 828.94 828.94 828.94 828.94 828.94

#Veh. 10 10 10 10 10 10 10 10 10

Dis. 828.94 828.94 828.06 824.78 828.94 828.94 828.94 828.94 828.94

#Veh. 10 10 10 10 10 10 10 10 10

Dis. 828.94 828.94 839.37 838.98 828.94 842.10 828.94 832.74 828.94

R101 R102 R103 R104 R105 R106 R107 R108 R109 R110 R111 R112

18 18 14 10 15 13 11 10 12 12 12 10

1613.59 1454.68 1235.68 974.24 1375.23 1260.2 1085.75 954.03 1157.74 1104.56 1057.8 974.73

19 17 14 10 14 13 11 10 12 11 11 10

1690.28 1513.74 1237.05 1020.87 1415.13 1254.22 1100.52 975.34 1169.85 1112.21 1084.76 976.99

19 17 13 9 14 12 10 9 11 10 10 9

1645.79 1486.12 1292.68 1007.24 1377.11 1251.98 1104.66 960.88 1194.73 1118.59 1096.72 982.14

19 17 13 9 14 12 10 9 11 10 10 9

1650.80 1486.12 1292.67 1007.31 1377.11 1252.03 1104.66 960.88 1194.73 1118.84 1096.73 982.14

18 16 12 10 14 12 11 10 12 11 11 10

1612.29 1473.41 1279.37 1025.47 1377.95 1276.48 1153.86 990.82 1179.73 1113.10 1155.39 981.46

20 18 15 11 16 13 12 10 13 12 12 10

1646.90 1474.28 1222.68 989.53 1382.78 1250.11 1083.42 952.44 1160.69 1080.69 1057.64 965.00

19 17 14 10 14 12 11 10 12 11 10 10

1669.34 1443.18 1238.56 974.24 1382.50 1252.03 1108.80 960.88 1222.70 1114.56 1096.73 976.99

RC101 RC102 RC103 RC104 RC105 RC106 RC107 RC108

16 13 11 10 14 13 11 11

1641.65 1470.26 1267.86 1145.49 1589.91 1371.69 1222.16 1133.9

15 14 12 10 14 12 12 10

1636.92 1488.36 1306.42 1140.45 1616.56 1454.61 1254.26 1141.34

14 12 11 10 13 11 11 10

1696.94 1554.75 1261.67 1135.48 1629.44 1424.73 1230.48 1139.82

14 12 11 10 13 11 11 10

1696.95 1554.75 1261.67 1135.48 1629.44 1424.73 1230.48 1139.82

15 13 11 11 13 12 11 11

1671.54 1447.14 1313.79 1163.54 1502.48 1406.25 1278.96 1172.83

16 15 12 10 16 14 12 11

1660.55 1494.92 1276.05 1151.63 1556.21 1402.25 1212.83 1133.25

15 13 11 10 16 11 11 10

1645.67 1476.50 1262.01 1135.48 1589.40 1424.73 1230.48 1142.66

AN US

341 33046.17 202492.34

325 33047.37 196094.74

325 33052.76 196105.52

#Veh. 10 10 10 10 10 10 10 10 10

Dis. 828.94 828.94 828.06 824.78 828.94 828.94 828.94 828.94 828.94

334 33073.75 199747.50

358 32609.27 208418.54

337 32802.86 200405.72

M

344 32549.78 202699.56

Dis. 828.94 828.94 828.06 824.78 828.94 828.94 828.94 828.94 828.94

CR IP T

C101 C102 C103 C104 C105 C106 C107 C108 C109

Dis. 828.93 828.19 828.06 825.54 828.9 828.17 829.34 832.28 829.22

Tol. Tol. Cost

#Veh. 10 10 10 10 10 10 10 10 10

HACO

#Veh. 10 10 10 10 10 10 10 10 10

Dis. 591.58 591.56 593.25 595.55 588.16 588.49 588.88 588.03

CE

C201 C202 C203 C204 C205 C206 C207 C208

HMOEA #Veh. 3 3 3 3 3 3 3 3

MOGA

MACS-IH

HSFLA

ILNSA

LGA

HACO

#Veh. 3 3 3 3 3 3 3 3

Dis. 591.56 591.56 596.55 590.6 588.88 588.49 588.29 588.32

#Veh. 3 3 3 3 3 3 3 3

Dis. 591.56 591.56 591.17 590.6 588.88 588.49 588.29 588.32

#Veh. 3 3 3 3 3 3 3 3

Dis. 591.56 591.56 591.17 590.60 588.88 588.49 588.29 588.32

#Veh. 3 3 3 3 3 3 3 3

Dis. 591.56 591.56 591.17 591.17 588.88 588.49 588.29 591.39

#Veh. 3 3 3 3 3 3 3 3

Dis. 591.56 591.56 591.17 590.60 588.88 588.49 588.29 588.32

#Veh. 3 3 3 3 3 3 3 3

Dis. 591.56 591.56 591.17 598.72 588.16 588.90 592.46 588.32

PT

Instance

ED

Table 10. Comparisons between the HACO and state-of-the-art algorithms from literatures for the category of type 2 (C2, R2 and RC2).

5 4 4 3 3 3 3 2 3 5 4

1206.42 1091.21 935.04 789.72 1094.65 940.12 852.62 790.6 974.88 982.31 811.59

4 4 3 3 3 3 3 2 3 3 3

1268.44 1112.59 989.11 760.82 1084.34 919.73 825.07 773.13 971.7 985.38 833.76

4 3 3 2 3 3 2 2 3 3 2

1252.37 1191.7 939.54 825.52 994.42 906.14 890.61 726.75 909.16 939.34 892.71

4 3 3 2 3 3 2 2 3 3 2

1252.37 1191.70 939.50 825.52 994.43 906.14 890.61 726.82 909.16 939.37 885.71

5 4 3 3 3 3 3 3 3 3 3

1285.90 1195.10 980.72 801.68 1051.60 941.17 843.77 744.71 938.29 966.50 855.76

9 8 6 4 6 6 4 3 5 5 5

1156.30 1042.30 877.29 736.52 960.35 894.19 800.79 706.86 860.63 948.82 762.23

4 4 4 3 3 3 3 2 3 3 4

1252.37 1091.21 937.56 794.56 994.43 906.14 825.48 726.50 855.00 939.37 761.50

RC201 RC202 RC203 RC204 RC205 RC206 RC207 RC208

6 5 4 3 5 4 4 3

1134.91 1130.53 1026.61 879.82 1295.46 1139.55 1040.67 898.49

4 4 3 3 4 4 3 3

1423.73 1183.88 1131.78 806.44 1352.39 1269.64 1140.23 881.2

4 3 3 3 4 3 3 3

1406.91 1365.65 1049.62 798.41 1297.19 1146.32 1061.14 828.14

4 3 3 3 4 3 3 3

1406.94 1365.64 1049.62 798.46 1297.65 1146.32 1061.14 828.14

5 4 3 3 5 4 3 3

1399.20 1382.40 1097.60 828.55 1269.40 1123.00 1090.20 859.13

10 8 6 4 8 7 7 5

1281.63 1103.47 942.00 796.12 1168.89 1060.52 970.97 782.70

4 4 4 3 4 3 4 4

1423.70 1263.53 1048.00 799.52 1410.30 1146.32 1140.67 828.71

AC

R201 R202 R203 R204 R205 R206 R207 R208 R209 R210 R211

Tol. Tol. Cost

97 23740.70 86281.40

86 24437.61 83275.22

80 24140.51 80281.02

80 24134.11 80268.22

20

90 24377.19 84754.38

140 22571.45 101142.90

90 23875.72 83751.44

ACCEPTED MANUSCRIPT

CR IP T

We also compare the results obtained by the proposed HACO with the results obtained by 6 stateof-the-art algorithms from the literatures, which are hybrid multi-objective evolutionary algorithm (HMOEA) [45], multi-objective genetic algorithm (MOGA) [34], multiple ant colony system algorithm hybridized with insertion heuristics (MACS-IH) [4], hybrid shuffled frog leaping algorithm (HSFLA) [32], improved large neighborhood search algorithm (ILNSA) [23] and localized genetic algorithm (LGA) [47]. These 6 algorithms that are used in the comparisons are heuristic algorithms with very effective strategies for avoiding local optimum, like insertion heuristics, mutation procedure, neighborhood search. In addition to these 6 algorithms, there is a large number of researchers that solved the vehicle routing problem with time windows with different heuristic algorithms [33]. However, if we present analytically all of the results obtained by the heuristic algorithms from the literatures, a huge table (or tables) will be produced and the most important outcomes will not be sufficiently presented.

AN US

Table 9 and Table 10 give the comparisons between the HACO and the comparing algorithms. #Veh. and Dis. have the same meanings as shown in Tables 7 and 8. Bold numbers indicate that the solutions obtained by HACO are very similar to or better than the results obtained by the comparing algorithms. In Tables 9 and 10, the second row from the bottom gives the total accumulated number of vehicles and the total traveled distance for each one of the algorithms. The total routing costs for each one of the algorithms are also given in the last rows of Tables 9 and 10.

PT

ED

M

As shown in Tables 9 and 10, the proposed HACO produces excellent routing results with 27 instances (out of the Solomon’s 56 instances) achieving a smallest number of vehicles and lowest traveled distance as compared with the best solutions obtained by other 6 algorithms. From Table 9, we can also see that HACO performs as well as MACS-IH, HSFLA and LGA on C1 (C101-C109) instances. These 4 algorithms succeed in finding the best known solutions on C1 problems which have all customers located in geographical location and tight time windows. The total costs obtained by the proposed HACO for the Somolon’s 56 instances are better than that produced by HMOEA, MOGA, ILNSA and LGA. The number of vehicles and the total traveled distance are two objectives of the VRPTW. In [4] and [32], the results have the following hierarchical objectives: (1) minimizing the number of vehicles and (2)minimizing the total distance, which result in lower number of vehicles and smaller routing costs obtained by MACS-IH and HSFLA. As shown before, the proposed HACO does not use hierarchical computational procedure to solve the MOVRPFlexTW with wbli = wbui = 0 (i ∈ {1, 2, · · · , n}). Thus, we can say that the proposed HACO still appears relatively competitive.

AC

CE

To further investigate the performance of the algorithms, statistical tests are also conducted as experimental tests in the current study. All computations related to the statistical tests were performed with the statistical software SPSS. Given that nonparametric tests do not require explicit conditions to be conducted, the nonparametric Friedman test is employed to determine whether the differences observed between the proposed HACO and those comparing algorithms are statistically significant or not. A wider description of nonparametric tests and their usage for analysis of algorithms’ behavior are presented in [20]. The Gapcost (Equation (22)) of each instance is considered for the statistical tests. Gapcost =

cost obtained by the given algorithm − best cost × 100% best known cost

(22)

where the cost means the routing costs. Best cost is calculated on the basis of the best known solution presented in Table 7 and Table 8, and the cost obtained by the given algorithm is calculated on the basis of the solution presented in Table 9 and Table 10. Significance level α = 0.05 has been assumed. Table 11 presents average ranks calculated for each algorithm, the Friedman χ2F statistic, p-value, and decision about acceptance of zero hypothesis for 21

ACCEPTED MANUSCRIPT

Table 11. Results of the Friedman test. Average rank 4.71 4.35 2.37 2.42 4.63 5.62 3.91 131.18 0.00 No

CR IP T

Algorithm HMOEA MOGA MACS-IH HSFLA ILNSA LGA HACO χ2F p-value Accept ?

Table 12. Results of pairwise comparison using the Wilcoxon signed-ranks nonparametric test. HACO vs HMOEA MOGA MACS-IH HSFLA ILNSA LGA

p-value 0.034 0.317 1.000 1.000 0.068 1.000

Accept ? No Yes Yes Yes Yes Yes

Set of instances R2

C2

HMOEA MOGA MACS-IH HSFLA ILNSA LGA

0.225 0.686 0.273 0.273 0.500 0.273

Yes Yes Yes Yes Yes Yes

RC1

R1

HMOEA MOGA MACS-IH HSFLA ILNSA LGA

0.047 0.041 0.013 0.013 0.049 0.005

No No No No No No

RC2

HACO vs HMOEA MOGA MACS-IH HSFLA ILNSA LGA

p-value 0.013 0.041 0.063 0.063 0.016 0.003

Accept ? No No Yes Yes No No

HMOEA MOGA MACS-IH HSFLA ILNSA LGA

0.033 0.040 0.043 0.043 0.026 0.025

No No No No No No

HMOEA MOGA MACS-IH HSFLA ILNSA LGA

0.047 0.031 0.018 0.018 0.067 0.012

No No No No Yes No

M

AN US

Set of instances C1

ED

the whole set of instances. Because of the fact that the value of χ2F is equal to 131.18 (critical value of χ2F with 6 degrees of freedom is equal to 12.69), zero hypothesis has to be rejected at assumed level of significance. Moreover, the p-value (0.00 < 0.05) strongly suggests the existence of significant difference among the considered algorithms.

CE

PT

In order to detect differences between particular algorithms, the Wilcoxon signed-ranks nonparametric test [40] is employed to compare two algorithms for each instance (pairwise comparison). Results of the Wilcoxon tests are presented in Table 12, where the p-values are presented for each set of instances and for each pair of algorithms. Similar as previously, the last column indicates whether zero hypothesis should be accepted or rejected for each pair of compared algorithms within each set of instances.

AC

As results in Table 12 show, pairwise comparison of the HACO with other 6 comparing algorithms does not always confirm significant differences between it and other ones. Whereas significant differences are observed for each compared pair of algorithms for R1 and RC1 set of instances. It means that HACO outperforms other 6 comparing algorithms while solving instances belonging to sets R1 and RC1. For R2 and RC2 set of instances, differences can be also observed for 4 and 5 compared pair of algorithms, respectively. In the case of C1 set of instances, difference are observed between HACO and HMOEA. An interesting observation from the analysis is that for instances belonging to C2 set, zero hypotheses have to be accepted where HACO with other 6 heuristic algorithms are compared. It means that HACO and other 6 comparing algorithms statistically perform similarly while solving instances belonging to this set.

22

ACCEPTED MANUSCRIPT

5.3.2

Solutions obtained by applying HACO to solve MOVRPFlexTW

Table 13. Comparison of the single objective VRPFlexTW solutions with the solutions obtained by implementing MOVRPFlexTW with the suggestive HACO where wbli = wbui = 0.05 (i ∈ {1, 2, · · · , n}). MOVRPFlexTW+HACO Ins.

Single obj. VRPFlexTW

Dis.

Satis.

Cost

#Veh.

Dis.

828.94 838.04 864.64 848.98 828.94 828.94 828.94 828.94 828.94

10 10 10 10 10 10 10 10 10

828.94 828.94 838.56 856.04 828.94 828.94 828.94 828.94 828.94

1.000 1.000 0.944 0.968 1.000 1.000 1.000 1.000 1.000

5657.88 5657.88 5677.12 5712.08 5657.88 5657.88 5657.88 5657.88 5657.88

10.00 10.00 10.30 10.50 10.00 10.00 10.00 10.00 10.00

828.94 828.94 849.17 862.30 828.94 828.94 828.94 835.90 828.94

Time (Sec.) 428 497 508 485 1029 1352 1469 1584 1720

R101

19

1669.34

570 726

17.00

1530.74

1484

R103

14

1232.08

1539

14.30

1244.84

1297

R104 R105

10 15

994.97 1370.79

1644 1395

10.70 15.20

1007.63 1390.75

1384 1492

R106

13

1240.71

1413

12.40

1267.35

1373

R107

11

1091.73

1079

11.30

1135.46

1436

R108

10

952.87

786

10.50

966.34

1275

R109

13

1163.73

1069

R110

11

1073.11

1456

R111

12

1077.80

460

R112

10

967.98

673

10947.16 10973.08 9839.34 9893.52 8077.12 8106.50 5988.48 8365.00 8393.86 7304.06 7373.82 6617.60 6795.00 5921.76 5947.34 7245.40 7337.04 6635.30 6779.38 6608.74 6655.66 5959.98 5989.34

1598

1509.70

0.914 0.935 0.886 0.908 0.973 0.984 0.993 0.971 0.986 0.938 0.957 0.964 0.981 0.967 0.981 0.884 0.902 0.943 0.965 0.918 0.924 0.936 0.948

1684.63

17

1673.58 1686.54 1519.67 1546.76 1238.56 1253.25 994.24 1382.50 1396.93 1252.03 1286.91 1108.80 1197.50 960.88 973.67 1222.70 1268.52 1117.65 1189.69 1104.37 1127.83 979.99 994.67

19.00

R102

19 19 17 17 14 14 10 14 14 12 12 11 11 10 10 12 12 11 11 11 11 10 10

RC101

15

1629.99

1190 1078

1506.30

1368

11

1265.89

11.50

1283.50

1365

RC104

10

1152.08

10.60

1149.56

1307

RC105

15

1538.94

16.30

1603.27

1325

RC106

12

RC107

11

RC108

10

9291.34 9568.80 8159.04 8297.26 6931.78 6979.30 6277.88 6314.44 9578.80 9623.00 7254.26 7356.34 6860.96 7133.02 6279.64 6306.94

13.40

RC103

0.951 0.962 0.963 0.982 0.974 0.985 0.968 0.980 0.982 0.991 0.915 0.926 0.904 0.933 0.941 0.965

1496

1493.60

1645.67 1784.40 1479.52 1548.63 1265.89 1289.65 1138.94 1157.22 1589.40 1611.50 1427.13 1478.17 1230.48 1366.51 1139.82 1153.47

1706.40

13

15 15 13 13 11 11 10 10 16 16 11 11 11 11 10 10

15.70

RC102

12.00

1214.18

0.958

7228.37

11.79

1428 1321

1341.22

777

1194.74

775

1113.39

1989

CE

Ave.

775

1124.14

958.62

AN US

10 10 10 10 10 10 10 10 10

CR IP T

#Veh.

M

C101 C102 C103 C104 C105 C106 C107 C108 C109

Average in 10 runs

Time (Sec.) 87 492 629 1468 637 709 593 539 503

ED

Dis.

PT

#Veh.

Non-dominated solutions

12.70

1248.53

1304

11.60

1174.50

1322

12.10

1136.70

1317

10.00

993.20

1266

11.70

1463.28

1315

11.50

1318.49

1285

10.80

1144.53

1308

12.04

1154.38

1254.74

AC

In this subsection, we evaluate the benefits gained by MOVRPFlexTW. Tables 13-15 provide the solutions for the single objective VRPFlexTW (Single obj. VRPFlexTW) reported in the literature [43], as well as the non-dominated solutions and average results produced by the proposed HACO method for MOVRPFlexTW over 10 runs. These tables present the total distance (Dis.), number of vehicles (#Veh.), total routing cost (Cost), and average customer satisfactions (Satis.) for: (i) MOVRPFlexTW with wbli = wbui = 0.05; (ii) MOVRPFlexTW with wbli = wbui = 0.15; and (iii) MOVRPFlexTW with wbli = wbui = 0.25. In Tables 13-15, the solutions produced using MOVRPFlexTW and HACO are written in boldface if they are very similar to or better than the solutions of the single objective VRPFlexTW presented in [43]. Note that the single objective VRPFlexTW was solved using the tabu search algorithm, and its objective was to minimize the total cost, which consists of traveling costs, fixed vehicle costs, and penalty costs incurred for time window violations [43]. In this paper, time window violations are reflected by customer satisfaction, which is maximized as a separate objective. 23

ACCEPTED MANUSCRIPT

Problem types C2, R2, and RC2 are omitted here because these instances require more computational time if their wider time windows are further relaxed.

CR IP T

In this proposed approach, a decision maker’s preferences among the objectives (customer satisfaction and routing costs) are taken into consideration to choose among the solutions. Table 13 gives the results for MOVRPFlexTW with wbli = wbui = 0.05. In Table 13, a single Pareto solution was obtained for problems C101-C109 and R104. For C101, C102, C103, C105, C106, C107, C108, C109, and R104, their solutions obtained using MOVRPFlexTW were the same as or better than the solutions obtained using the single objective VRPFlexTW. For C104, although the distance obtained using MOVRPFlexTW was slightly worse than the distance obtained by using single objective VRPFlexTW, the difference was only 0.83%. For the other problems, different non-dominated solutions were obtained using MOVRPFlexTW. Compared with the number of vehicles obtained using the single objective VRPFlexTW, MOVRPFlexTW implemented using HACO performed better for R105, R106, R109, R111, and RC106.

AN US

Tables 14 and 15 present the solutions obtained for the MOVRPFlexTW with wbli = wbui = 0.15 for all i ∈ {1, 2, · · · , n} and with wbli = wbui = 0.25 for all i ∈ {1, 2, · · · , n}, respectively. Compared with the solutions obtained using the single objective VRPFlexTW, the non-dominated solutions of C102, C103 and RC103 in Table 14 and the non-dominated solutions of C102, C103, R106, R107, R109, R110, RC102, RC103, RC104, RC105 and RC108 in Table 15 are much better in terms of distance and the number of vehicles, but the MOVRPFlexTW cannot always provide similar or better solutions for all instances.

Effects of implementing different time windows

PT

5.3.3

ED

M

As shown in Tables 13-15, although a part of the non-dominated solutions obtained using MOVRPFlexTW are not better than the solutions obtained using the single objective VRPFlexTW, we should note that these alternate Pareto solutions are clearly comparable and a decision maker can decide which solution is more preferable based on his/her preferences. Moreover, the proposed MOVRPFlexTW eliminates the need to experiment with weights for the time window violations (as required in a weighted-sum approach), which is useful because poorly chosen weights can result in unsatisfactory solutions. Furthermore, the solutions and computational times in Tables 13-15 confirm that the proposed HACO method can provide good solutions in a reasonable computational time.

CE

From the numbers shown in Tables 13-15, it is easy to find that the average computational time is increased with the enlarged time windows. For example, the average computational time of solving the instance C101 where wbli = wbui = 0.05 is 428 seconds over 10 runs, while the average computational time of solving it where wbli = wbui = 0.15 and wbli = wbui = 0.25 are 678 seconds and 1017 seconds, respectively. This is because the wider relaxed time windows require more computational time.

AC

Comparing the average results in Table 13 and 14 with those in Table 15, it is clear that larger flexibility fractions reduce the total distance using fewer vehicles where the customer satisfaction is lower. We can similarly compare the average results in Table 13 with those presented in Table 14. This indicates that the decision maker can increase customer satisfaction at the expense of routing costs using narrow flexible time windows; on the contrary, the decision maker can decrease the routing costs at the expense of customer satisfaction using wide flexible time windows.

24

ACCEPTED MANUSCRIPT

Table 14. Comparison of the single objective VRPFlexTW solutions with the solutions obtained by implementing MOVRPFlexTW with the suggestive HACO where wbli = wbui = 0.15 (i ∈ {1, 2, · · · , n}). MOVRPFlexTW+HACO Non-dominated solutions

C101 C102

10 10

828.94 1007.31

Time (Sec.) 315 712

C103

10

972.41

617

C104 C105 C106 C107 C108 C109

10 10 10 10 10 10

855.29 828.94 828.94 828.94 828.94 828.94

1575 681 672 573 570 636

R101

18

1623.69

1196

R102

16

1425.14

1236

R103

14

1210.14

2112

R104

10

980.65

1523

R105

14

1343.93

1938

R106

12

1186.79

1435

R107

11

1050.86

1726

9

961.70

1390

R109

12

1128.15

1938

R110

11

1033.29

1723

R111

11

1068.23

1017

R112

10

965.95

648

RC101

15

1595.05

1512

RC102

13

1428.71

2127

RC103

12

1279.81

RC104

10

1123.85

RC105

13

1457.73

RC106

11

RC107

11

RC108

10

Ave.

11.48

742

917

PT

1128

1298.27

995

1183.25

1895

1094.27

1959

1112.00

1224.41

CE

Dis.

Satis.

Cost

#Veh.

Dis.

10 10 10 10 10 10 10 10 10 10 10

828.94 857.64 878.04 864.64 937.68 848.98 828.94 828.94 828.94 828.94 828.94

1.000 0.894 0.916 0.883 0.907 0.952 1.000 1.000 1.000 1.000 1.000

5657.88 5715.28 5756.08 5729.28 5875.36 5697.96 5657.88 5657.88 5657.88 5657.88 5657.88

10.0 10.4

828.94 864.37

Time (Sec.) 678 649

10.0

907.48

1008

10.2 10.3 10.0 10.3 10.6 10.3

867.54 833.64 828.94 835.46 843.74 835.84

895 1343 1799 1775 1795 1933

18 18 18 16 16 14 14 10 10 14 14 12 12 11 11 9 9 12 12 11 11 11 11 10 10

1654.87 1661.59 1670.31 1445.23 1478.56 1227.35 1234.32 980.65 987.60 1343.93 1356.87 1189.80 1207.99 1075.80 1123.30 961.70 987.70 1130.64 1167.86 1038.67 1089.60 1064.25 1087.43 963.20 981.90

0.836 0.874 0.891 0.783 0.805 0.942 0.961 0.954 0.967 0.918 0.937 0.881 0.903 0.928 0.943 0.868 0.908 0.847 0.872 0.854 0.881 0.848 0.869 0.886 0.905

10509.74 10523.18 10540.62 9290.46 9357.12 8054.70 8068.64 5961.30 5975.20 8287.86 8313.74 7179.60 7215.98 6551.60 6646.60 5523.40 5575.40 7061.28 7135.72 6477.34 6579.20 6528.50 6574.86 5926.40 5963.80

18.6

1663.20

1944

16.6

1463.84

1866

14.3

1229.84

1755

10.5

983.56

1865

14.8

1352.74

1934

12.0

1196.58

1637

11.0

1104.67

1698

15 15 13 13 11 11 10 10 13 13 11 11 11 11 11 10 10

1595.05 1632.08 1436.89 1467.28 1154.60 1189.50 1083.90 1134.60 1484.60 1509.80 1298.27 1328.94 1209.64 1223.57 1237.60 1104.52 1124.60

0.928 0.946 0.903 0.937 0.845 0.897 0.937 0.957 0.728 0.849 0.774 0.795 0.853 0.868 0.886 0.897 0.918

9190.10 9264.16 8073.78 8134.56 6709.20 6779.00 6167.80 6269.20 8169.20 8219.60 6996.54 7057.88 6819.28 6847.14 6875.20 6209.04 6249.20

11.75

1163.91

0.901

7029.70

ED

R108

Average in 10 runs

#Veh.

CR IP T

Dis.

AN US

#Veh.

10.1

965.40

1549

12.3

1159.85

1798

11.3

1065.74

1965

11.8

1083.59

1985

10.0

970.13

1436

15.2

1608.59

1988

13.1

1449.78

1976

11.2

1176.43

1798

10.0

1095.64

1838

14.2

1498.46

1848

11.0

1308.68

1766

11.3

1230.58

1895

10.3

1116.95

1696

11.78

1116.21

1658.96

A real case application

AC

6

Single obj. VRPFlexTW

M

Instance

In this section, a real case application has been developed for Vanguard logistics which is a leading company in the vegetable and food distribution sector in China, with the purpose of validating the multi-objective VRPFlexTW model. The problem consists of 16 delivery points (16 supermarkets around Beilin district, Xi’an) served directly from a depot (supply base of vegetable and food). Figure 8 illustrates the geographical position of 16 delivery points. Owing that the depot (longitude 108.611118 and latitude 34.117327) is located in Hu county and far away from the 16 delivery points, it is not marked in Figure 8. The distance of transportation between each pairs of points are listed in Table 16 (0 is the index of depot). The fleet of vehicles consists of several same rigid trucks with capacity of 25

ACCEPTED MANUSCRIPT

Table 15. Comparison of the single objective VRPFlexTW solutions with the solutions obtained by implementing MOVRPFlexTW with the suggestive HACO where wbli = wbui = 0.25 (i ∈ {1, 2, · · · , n}). MOVRPFlexTW+HACO Single obj. VRPFlexTW

Non-dominated solutions

10 10

828.94 1017.06

C103

10

1022.88

634

C104

10

933.47

1741

C105 C106 C107 C108 C109

10 10 10 10 10

828.94 828.94 828.94 828.94 828.94

532 597 496 689 511

R101

17

1597.96

1117

R102

15

1405.93

664

R103

13

1225.99

1013

R104

10

1022.51

1013

R105

13

1294.36

1490

R106

12

1207.99

914

R107

11

1062.88

902

R108

9

946.54

1524

R109

12

1110.40

1301

R110

11

1061.10

559

R111

10

1053.24

887

R112

10

971.28

1416

RC101 RC102

14 13

1521.56 1395.25

1205 1835

RC103

10

1221.53

RC104

10

1133.50

RC105

13

1425.89

RC106

11

RC107

10

RC108

10

Ave.

11.17

1053

764

PT

1173 588

1157.88

1252

1127.99

740

1109.34

948.55

Dis.

Satis.

Cost

#Veh.

Dis.

10 10 10 10 10 10 10 10 10 10 10 10 10

828.94 903.65 946.85 918.12 974.20 1013.52 933.47 946.36 828.94 828.94 828.94 828.94 828.94

1.000 0.814 0.825 0.894 0.903 0.913 0.903 0.912 1.000 1.000 1.000 1.000 1.000

5657.88 5807.30 5893.70 5836.24 5948.39 6027.04 5866.94 5892.72 5657.88 5657.88 5657.88 5657.88 5657.88

10.0 10.0

828.94 927.65

Time (Sec.) 1017 977

10.0

946.82

1469

10.0

940.55

1136

10.3 10.0 10.2 10.0 10.0

830.24 828.94 830.76 831.25 828.94

1795 2007 2104 2201 2594

17 17 15 15 13 13 10 10 13 13 13 12 12 11 11 9 9 11 11 10 10 10 10 10 10

1568.93 1614.65 1405.93 1438.93 1211.09 1268.76 1008.54 1039.46 1304.50 1367.90 1441.20 1147.60 1164.90 1027.85 1048.39 947.20 968.70 1228.20 1285.60 1001.40 1108.30 1053.24 1063.52 976.80 981.90

0.693 0.775 0.672 0.775 0.628 0.635 0.906 0.915 0.746 0.771 0.793 0.813 0.865 0.875 0.894 0.846 0.873 0.658 0.673 0.516 0.579 0.558 0.584 0.894 0.905

9937.86 10029.30 8811.86 8877.86 7622.18 7737.52 6017.08 6078.92 7809.00 7935.80 8082.40 7095.20 7129.80 6455.70 6496.78 5494.40 5537.40 6856.40 6971.20 6002.80 6216.60 6106.48 6127.04 5953.60 5963.80

17.8

1593.28

2568

15.6

1427.35

2365

13.7

1254.65

2290

10.0

1019.93

2229

13.4

1430.86

2498

12.0

1158.43

2009

11.0

1039.82

2014

9.0

950.48

2104

11.6

1277.90

2339

10.5

1106.85

2438

13 12 13 10 10 10 10 11 13 11 11 10 10 10 10

1788.40 1593.70 1394.60 1154.60 1175.20 1061.69 1124.20 1219.70 1359.60 1280.12 1298.27 1123.00 1173.47 1070.01 1071.60

0.509 0.458 0.716 0.764 0.782 0.903 0.946 0.485 0.697 0.615 0.774 0.597 0.625 0.784 0.796

11.11

1135.31

0.782

10.3

1060.25

2655

10.0

979.68

2186

8776.80 7987.40 7989.20 6309.20 6350.40 6123.38 6248.40 6839.40 7919.20 6960.24 6996.54 6246.00 6346.94 6140.02 6143.20

13.6 12.8

1793.82 1428.90

2646 2785

10.2

1170.98

2227

10.0

1087.32

2453

12.5

1298.30

2457

11.0

1294.57

2114

10.0

1162.49

2228

10.0

1070.30

2105

6715.90

11.22

1117.25

2138.21

AC

CE

1280.12

Average in 10 runs

#Veh.

CR IP T

C101 C102

Time (Sec.) 377 521

AN US

Dis.

ED

#Veh.

M

Ins.

80 tons and average speed of 60 kilometers per hour to deliver the supermarkets’ demands. Table 17 shows the depot and supermarkets’ details in sequence of supermarket index, extended earliest time to begin the service Ei , ready time ei , due time li , extended latest time to begin the service Li , daily demands and service time. As shown in Table 17, the time windows of almost all supermarkets are between 10:00 and 12:00 which are not rush hours, and therefore the traffic jam is not considered in this study. The maximum route travel time for each vehicle is 3 hours. The transportation cost and the fixed cost of using each vehicle are 3 Yuan per kilometer and 600 Yuan, respectively. Because the unpunctual services have negative effect on the sales of corresponding supermarket, all supermarkets usually require to be served between the ready time and the due time. However, owing to the smaller flow of customers during the traditional festival holidays such as Dragon Boat Festival, 26

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 8: Geographical positions of the 16 delivery points (supermarkets). Table 16. Distance of transportation between each pair of points (km). 2 8.000 1.650 0.000

3 13.000 4.750 4.218 0.000

4 7.000 2.042 1.773 5.507 0.000

5 16.000 0.936 5.068 3.230 6.495 0.000

6 10.000 3.111 4.859 6.842 4.133 7.873 0.000

7 14.000 7.112 5.520 8.129 4.952 8.934 9.843 0.000

8 11.000 4.956 3.364 5.973 3.388 6.778 7.687 2.382 0.000

M

1 7.000 0.000

ED

Index 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

9 15.000 6.524 6.012 1.855 7.439 3.321 7.598 6.402 5.001 0.000

10 15.000 6.545 6.033 1.796 7.460 3.192 7.619 6.423 5.570 0.267 0.000

11 5.000 3.339 4.383 7.070 3.125 8.101 3.474 4.837 7.137 9.064 8.796 0.000

12 14.000 7.105 6.534 4.436 6.287 7.639 9.536 5.281 2.907 4.232 3.976 9.920 0.000

13 60.000 5.144 6.730 6.295 6.775 6.862 3.431 8.765 10.000 8.289 8.021 6.388 9.887 0.000

14 15.000 8.540 8.008 3.975 9.396 3.828 10.000 8.418 7.226 2.546 2.847 11.000 6.600 9.629 0.000

15 10.000 8.590 9.730 8.410 9.816 6.026 6.447 11.000 13.000 8.653 9.224 9.954 11.000 5.135 6.106 0.000

16 80.000 6.745 7.885 7.498 7.978 5.824 4.609 9.968 11.000 9.025 9.224 8.116 11.000 2.987 8.067 2.013

CE

PT

Mid-Autumn Festival and Spring Festival, the time window constraints of supermarkets can be violated to a certain extent. The extended earliest time Ei and latest time Li of each supermarkets are presented in Table 17. Additionally, during the traditional festival holidays, the company needs to pay large extra allowance to the distribution clerks who still deliver demands. The decision maker of this company needs to design the vehicle routs with two different objectives: (1)minimizing the total routing costs, and (2) maximizing the total customer satisfaction.

AC

As stated before, the MOVRPFlexTW with wbli = wbui = 0 (i ∈ {1, 2, · · · , n}) is equivalent to the single objective vehicle routing problem with hard time window. To strictly satisfy the time window constraints of all delivery supermarkets, we firstly solve this problem by implementing MOVRPFlexTW with the proposed HACO where wbli = wbui = 0 (i ∈ {0, 1, 2, · · · , 16}), and the results are presented in Table 18.

From Table 18, we can know that the company usually needs to arrange 5 vehicles to deliver the demands within the appointed time windows [ei , li ] (i ∈ {0, 1, 2, · · · , 16}). Because of the narrow time windows, all the vehicles are not full-load when they dispatch from the depot. All the vehicles can return back to the depot before 12:00. According to the total distance 203.908km shown in the last row in Table 18 and the number of vehicles, we can obtain the daily routing cost which is around 3611.724 Yuan. 27

ACCEPTED MANUSCRIPT Table 17. Time windows and demands of each point. Ei

ei

li

Li

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

10:00 10:25 10:30 10:40 10:15 10:10 10:15 10:05 10:15 10:05 10:40 10:00 10:20 10:30 10:20 10:05 10:40

10:00 10:30 10:50 10:45 10:40 10:15 10:20 10:10 10:20 10:10 11:00 10:10 10:40 10:55 10:30 11:25 11:00

12:00 11:10 11:10 11:05 10:55 10:50 10:50 10:30 10:45 10:30 11:20 10:30 11:00 11:15 10:55 11:40 11:30

12:30 11:20 11:20 11:30 11:00 11:40 11:30 10:40 10:55 10:55 11:30 11:40 11:05 11:30 11:45 12:00 11:50

Demands (ton) 0 30 16 18 9 23 17 12 26 14 19 15 29 31 10 21 18

Service time (minute) 0 10 15 15 10 15 12 10 8 10 15 13 11 10 8 10 12

CR IP T

Point

Table 18. Result obtained by implementing MOVRPFlexTW to solve the proposed case where wbli = wbui = 0 (i ∈ {0, 1, 2, · · · , 16}). Route

1 2 3 4 5 Total

0-13-16-15-0 0-5-14-3-10-0 0-9-8-12-0 0-7-1-4-2-0 0-11-6-0

Distance (km) 75.000 40.599 36.908 32.927 18.474 203.908

Delivery time (hour) 1.790 1.560 1.090 1.295 0.724 6.459

Time of returning back to the depot 11:47 11:33 11:05 11:17 10:43

Sum of demands (ton) 70 70 69 67 32 308.000

AN US

Vehicle

PT

ED

M

For delivering the demands during the festival holidays, the MOVRPFlexTW provides a single Pareto solution which is presented in Table 19. The total customer satisfactions of each route are listed in the last column of Table 19. Owing to the violation of the appointed time windows [ei , li ] (i ∈ {0, 1, 2, · · · , 16}), the number of vehicles and the total distance presented in Table 19 (4 and 174.875km, respectively) are smaller than the corresponding results shown in Table 18, which result in smaller daily routing cost around 2924.625 Yuan. Furthermore, the total customer satisfaction is decreased from 16.00 to 10.65. Additionally, because of the decreased number of vehicles, the total delivery time is relatively increased and the time of returning back to the depot for each vehicle is also later in Table 19. After extending the time windows from [ei , li ] to [Ei , Li ], a part of vehicles, which is not full-load in the case of narrow time windows, can load more demands and serve more points. For example, the first vehicle in Table 18 returns to the depot after serving supermarket 15. After extending the time windows, supermarket 14 is added to this route (the first vehicle in Table 19) before returning to the depot.

CE

The above observations are very similar to the ones obtained in subsection 5.3.3. During the festival holidays, the decision maker can arrange a smaller number of vehicles to deliver the demands at the expense of supermarkets’ satisfaction using the flexible time windows.

AC

Table 19. Result obtained by implementing MOVRPFlexTW to solve the proposed case where the flexible time windows [Ei , Li ](i ∈ {0, 1, 2, · · · , 16}) are used. Vehicle

Route

1 2 3 4 Total

0-13-16-15-14-0 0-4-7-8-12-0 0-9-10-3-5-0 0-2-1-6-11-0

Distance (km) 86.106 31.241 36.293 21.235 174.875

Delivery time (hours) 2.100 1.300 1.767 1.550 6.717

Time of returning back to the depot 12:06 11:18 11:46 11:33

28

Sum of demands (ton) 80 76 74 78 308

Satisfaction 3.040 3.000 2.500 2.110 10.65

ACCEPTED MANUSCRIPT

7

Conclusion

In this paper, we considered a multi-objective VRPFlexTW, which is an extension of VRPTW, in which the customers are allowed to be serviced outside their required time interval with a certain tolerance. This problem also considers the customer satisfaction related to the route, which is represented as a convex function dependent on the vehicle arrival time. Compared with VRPTW, VRPFlexTW is more flexible and applicable to real life problems. This paper considered VRPFlexTW as a multi-objective problem, with the aim of minimizing the overall routing costs and maximizing the overall customer satisfaction.

AN US

CR IP T

The basic framework of the proposed approach for solving MOVRPFlexTW was inspired by the classical ACO meta-heuristic. However, this approach was extended using three mutation operators that enhance the local search ability and allow for random global explorations. These mutation operators prevent the ACO from getting trapped in local optimal solution. Meanwhile, using the concept of Pareto optimization, we can obtain better non-dominated solutions for MOVRPFlexTW. The proposed HACO was applied to solve the benchmark Solomon’s instances. Our results imply that the proposed HACO is quite sufficient as compared with the best published results. Because the MOVRPFlexTW uses more flexible time windows, more savings can be made for the logistic company at the expense of customer satisfaction. This trade-off must be carefully balanced by the decision makers with respect to the preferences of logistic company and the concerns of their customers.

PT

ED

M

A number of future research directions and ideas have arisen from this work. On the one hand, several new MOVRPFlexTW models can be obtained by reformulating objective functions, e.g., tripleobjective VRPFlexTW in which the total vehicles used for service and overall total traveled distance are minimized, and the overall customer satisfaction is maximized; and quadruple-objective VRPFlexTW in which the total vehicles used for service, overall total traveled distance and waiting time are minimized, and the overall customer satisfaction is maximized. On the other hand, due to the relatively high efficiency and generality of HACO, it can be easily adapted to solve some classic variants of VRP. For example, VRPTW (vehicle routing problem with hard time window), by not relaxing the time windows; HVRP ( heterogeneous fleet vehicle routing problem ), by separately checking the capacity constraint for each vehicle; and VRPSTW ( vehicle routing problem with soft time window), by limitlessly relaxing the customers’ time windows in the planning horizon. Therefore, we hope that the proposed HACO with required modification will be implemented to successfully solve these problems.

CE

Acknowledgement

AC

The authors would like to thank the anonymous referees and the editor for their constructive feedback. This work is supported by the National Natural Science Foundation of China (grant. 71401106), the Shanghai Natural Science Foundation (grant. 14ZR1418700), the Humanities and Social Sciences Foundation from the Ministry of Education of China (grants. 16YJA630037 and 19YJAZH064).

References [1] E. Alba and B. Dorronsoro. Computing nine new best-so-far solutions for capacitated vrp with a cellular genetic algorithm. Information Processing Letters, 98(6):225–230, 2006. [2] N. Balakrishnan. Simple heuristics for the vehicle routeing problem with soft time windows. Journal of the Operational Research Society, 44(3):279–287, 1993.

29

ACCEPTED MANUSCRIPT

[3] R. Baldacci, A. Mingozzi, and R. Roberti. Recent exact algorithms for solving the vehicle routing problem under capacity and time window constraints. European Journal of Operational Research, 218(1):1–6, 2012. [4] S. R. Balseiro, I. Loiseau, and J. Ramonet. An ant colony algorithm hybridized with insertion heuristics for the timedependent vehicle routing problem with time windows. Computers & Operations Research, 38(6):954–966, 2011.

CR IP T

[5] R. Ba˜ nos, J. Ortega, C. Gil, A. Fern´andez, and F. de Toro. A simulated annealing-based parallel multi-objective approach to vehicle routing problems with time windows. Expert Systems with Applications, 40(5):1696–1707, 2013. [6] R. Ba˜ nos, J. Ortega, C. Gil, A. L. M´ arquez, and F. de Toro. A hybrid meta-heuristic for multiobjective vehicle routing problems with time windows. Computers & Industrial Engineering, 65(2):286–296, 2013.

AN US

[7] A. K. Beheshti and S. R. Hejazi. A novel hybrid column generation-metaheuristic approach for the vehicle routing problem with general soft time window. Information Sciences, 316:598–615, 2015. [8] N. Bhusiri, A. G. Qureshi, and E. Taniguchi. The trade-off between fixed vehicle costs and timedependent arrival penalties in a routing problem. Transportation Research Part E: Logistics and Transportation Review, 62:1–22, 2014. [9] K. Braekers, K. Ramaekers, and I. V. Nieuwenhuyse. The vehicle routing problem: State of the art classification and review. Computers & Industrial Engineering, 99:300–313, 2016.

M

[10] O. Br¨aysy and M. Gendreau. Vehicle routing problem with time windows, part 1: Route construction and local search algorithms. Transportation Science, 39(1):104–118, 2005.

ED

[11] O. Br¨aysy and M. Gendreau. Vehicle routing problem with time windows, part 2: Metaheuristics. Transportation Science, 39(1):119–139, 2005.

PT

[12] H. I. Calvete, C. Gal´e, M. Oliveros, and B. S´anchez-Valverde. A goal programming approach to vehicle routing problems with soft time windows. European Journal of Operational Research, 177(3):1720–1733, 2007.

CE

[13] L. Chen, H.-Y. Sun, and S. Wang. A parallel ant colony algorithm on massively parallel processors and its convergence analysis for the travelling salesman problem. Information Sciences, 199:31–42, 2012.

AC

[14] W.-C. Chiang and R. A. Russell. A metaheuristic for the vehicle-routeing problem with soft time windows. Journal of the Operational Research Society, 55(12):1298–1310, 2004. [15] Q. Ding, X. Hu, L. Sun, and Y. Wang. An improved ant colony optimization and its application to vehicle routing problem with time windows. Neurocomputing, 98:101–107, 2012.

[16] M. Dorigo and L. M. Gambardella. Ant colony system: A cooperative approach to the traveling salesman problem. IEEE Transactions on Evolutionary Computation, 1(1):53–66, 1997. [17] A. E. Eiben and S. K. Smit. Parameter tuning for configuring and analyzing evolutionary algorithms. Swarm and Evolutionary Computation, 1(1):19–31, 2011. [18] N. A. El-Sherbeny. Vehicle routing with time windows: An overview of exact, heuristic and metaheuristic methods. Journal of King Saud University (Science), 22(3):123–131, 2010. 30

ACCEPTED MANUSCRIPT

[19] M. A. Figliozzi. An iterative route construction and improvement algorithm for the vehicle routing problem with soft time windows. Transportation Research Part C: Emerging Technologies, 18(5):668–679, 2010. [20] S. Garcia, D. Molina, M. Lozano, and F. Herrera. A study on the use of non-parametric tests for analyzing the evolutionary algorithms’s behaviour: a case study on the cec’2005 special session on real parameter optimization. Journal of Heuristics, 15:617–644, 2009.

CR IP T

[21] S. F. Ghannadpour, S. Noori, and R. Tavakkoli-Moghaddam. A multi-objective dynamic vehicle routing problem with fuzzy time windows: Model, solution and application. Applied soft computing, 14:504–527, 2014. [22] K. Ghoseiri and S. F. Ghannadpour. Multi-objective vehicle routing problem with time windows using goal programming and genetic algorithm. Applied Soft Computing, 10(4):1096–1107, 2010. [23] L. Hong. An improved lns algorithm for real-time vehicle routing problem with time windows. Computers & Operations Research, 39(2):151–163, 2012.

AN US

[24] T. Ibaraki, S. Imahori, M. Kubo, T. Masuda, T. Uno, and M. Yagiura. Effective local search algorithms for routing and scheduling problems with general time-window constraints. Transportation Science, 39(2):206–232, 2005. [25] S. Iqbal, M. Kaykobad, and M. S. Rahman. Solving the multi-objective vehicle routing problem with soft time windows with the help of bees. Swarm and Evolutionary Computation, 24:50–64, 2015.

M

[26] O. Jabali, R. Leus, T. Van Woensel, and T. de Kok. Self-imposed time windows in vehicle routing problems. OR Spectrum, 37(2):331–352, 2015.

ED

[27] N. Jozefowiez, F. Semet, and E. Talbi. Multi-objective vehicle routing problems. European Journal of Operational Research, 189(2):293–309, 2008.

PT

[28] C. B. Kalayci and C. Kaya. An ant colony system empowered variable neighborhood search algorithm for the vehicle routing problem with simultaneous pickup and delivery. Expert Systems with Applications, 66:163–175, 2016.

CE

[29] C. Koc¸, T. Bektas¸, O. Jabali, and G. Laporte. Thirty years of heterogeneous vehicle routing. European Journal of Operational Research, 249(1):1–21, 2016.

AC

[30] Y. A. Koskosidis, W. B. Powell, and M. M. Solomon. An optimization-based heuristic for vehicle routing and scheduling with soft time window constraints. Transportation Science, 26(2):69–85, 1992. [31] A. A. Kovacs, S. N. Parragh, and R. F. Hartl. The multi-objective generalized consistent vehicle routing problem. European Journal of Operational Research, 247(2):441–458, 2015. [32] J. Luo, X. Li, M.-R. Chen, and H. Liu. A novel hybrid shuffled frog leaping algorithm for vehicle routing problem with time windows. Information Sciences, 316:266–292, 2015. [33] Y. Marinakis, M. Marinaki, and A. Migdalas. A multi-adaptive particle swarm optimization for the vehicle routing problem with time windows. Information Sciences, 481:311–329, 2019. [34] B. Ombuki, B. J. Ross, and F. Hanshar. Multi-objective genetic algorithm for vehicle routing problem with time windows. Applied Intelligence, 24(1):17–30, 2006. 31

ACCEPTED MANUSCRIPT

[35] A. G. Qureshi, E. Taniguchi, and T. Yamada. An exact solution approach for vehicle routing and scheduling problems with soft time windows. Transportation Research Part E: Logistics and Transportation Review, 45(6):960–977, 2009. [36] A. G. Qureshi, E. Taniguchi, and T. Yamada. Exact solution for the vehicle routing problem with semi soft time windows and its application. Procedia - Social and Behavioral Sciences, 2(3):5931–5943, 2010. [37] A. Rajasekhar, N. Lynn, S. Das, and P. N. Suganthan. Computing with the collective intelligence of honey bees-a survey. Swarm and Evolutionary Computation, 32:25–48, 2017.

CR IP T

[38] M. Reed, A. Yiannakou, and R. Evering. An ant colony algorithm for the multi-compartment vehicle routing problem. Applied Soft Computing, 15:169–176, 2014. [39] Z.-G. Ren, Z.-R. Feng, and A.-M. Zhang. Fusing ant colony optimization with lagrangian relaxation for the multiple-choice multidimensional knapsack problem. Information Sciences, 182(1):15–29, 2012.

AN US

[40] S. Rostami, D. O’Reilly, A. Shenfield, and N. Bowring. A novel preference articulation operator for the evolutionary multi-objective optimisation of classifiers in concealed weapons detection. Information Sciences, 295:494–520, 2015. [41] M. Salani, M. Battarra, and L. M. Gambardella. Operations Research Proceedings 2014, chapter Exact Algorithms for the Vehicle Routing Problem with Soft Time Windows, pages 481–486. Springer International Publishing, 2016.

M

[42] D. Tas¸, N. Dellaert, T. Van Woensel, and T. de Kok. Vehicle routing problem with stochastic travel times including soft time windows and service costs. Computers & Operations Research, 40(1):214–224, 2013.

ED

[43] D. Tas¸, O. Jabali, and T. Van Woensel. A vehicle routing problem with flexible time windows. Computers & Operations Research, 52(Part A):39–54, 2014. [44] E. Taillard, P. Badeau, M. Gendreau, F. Guertin, and J.Y. Potvin. A tabu search heuristic for the vehicle routing problem with soft time windows. Transportation Science, 31(2):170–186, 1997.

PT

[45] K. C. Tan, Y. H. Chew, and L. H. Lee. A hybrid multiobjective evolutionary algorithm for solving vehicle routing problem with time windows. Computational Optimization and Applications, 34(1):115–151, 2006.

CE

[46] E. Teymourian, V. Kayvanfar, GH. M. Komaki, and M. Zandieh. Enhanced intelligent water drops and cuckoo search algorithms for solving the capacitated vehicle routing problem. Information Sciences, 334-335:354–378, 2016.

AC

[47] Z. Ursani, D. Essam, D. Cornforth, and R. Stocker. Localized genetic algorithm for vehicle routing problem with time windows. Applied Soft Computing, 11(8):5375–5390, 2011. [48] T. Vidal. Technical note: Split algorithm in o(n) for the capacitated vehicle routing problem. Computers & Operations Research, 69:40–47, 2016. [49] W. Wu, Y. Tian, and T. Jin. A label based ant colony algorithm for heterogeneous vehicle routing with mixed backhaul. Applied Soft Computing, 47:224–234, 2016. [50] E. Zare-Reisabadi and S. H. Mirmohammadi. Site dependent vehicle routing problem with soft time window: Modeling and solution approach. Computers & Industrial Engineering, 90:177– 185, 2015.

32