Time-Dependent Traveling Salesman Problem with Multiple Time Windows

—The TSP, VRP and OP problems with time constraints have one common sub-problem – the task of ﬁnding the minimum route duration for a given order of customers. While much work has been done on routing and scheduling problems with time windows, to this date only few articles considered problems with multiple time windows. Moreover, since the assumption of constant travel time between two locations at all times is very unrealistic, problems with time-dependent travel were introduced and studied. Finally, it is also possible to imagine some situations, in which the service time changes during the day. Again, both issues have been investigated only in conjunction with single time windows. In this paper we propose a novel algorithm for computing minimum route duration in traveling salesman problem with multiple time windows and time-dependent travel and service time. The algorithm can be applied to wide range of problems in which a traveler has to visit a set of customers or locations within speciﬁed time windows taking into account the trafﬁc and variable service/visit time. Furthermore, we compare three metaheuristics for computing a one-day schedule for this problem, and show that it can be solved very efﬁciently.


I. INTRODUCTION
I N this paper we focus our work on time-dependent routing and scheduling problem with multiple time windows.The problem consist of an agent (tourist, sales representative, etc.) whose aim/duty is to visit a predefined set of customers/locations (e.g. points of interest).Each customer/location may define many time windows which indicates the availability during the day.Furthermore, we assume, that the travel time between the customers/locations changes due to the traffic.In this work we also assume different service/visit time in different time windows.
In this work we describe the theory and algorithms for computing one-day schedule in time-dependent traveling salesman problem with multiple time windows for application in many well known operational research problems such as vehicle routing problem (VRP) (see [2], [3], [8], [14]), orienteering problem (OP) (see [15], [16]), or generally traveling salesman problem (TSP) with complex time constraints.While much work has been done on mixed routing and scheduling problems with time windows, to this date only few articles considered problems with multiple time windows (cf.[2]).
Throughout this article, we will denote a sequence of customers as a route, while we use the term schedule to denote a route with fixed visit times.
The paper is organized as follows: in Section 2 we describe the problem, and discuss additional issues arising from multiple time windows, and time-dependent travel and service time.Section 3 describes preprocessing of time windows and presents the minimum route duration algorithm.In Section 4 we explain in details the algorithms used for computing the one-day schedule in time-dependent traveling salesman problem with multiple time windows.Section 5 shows the results of our numerical experiments.Finally, some concluding remarks are given in Section 6.

II. PROBLEM DESCRIPTION
We consider a Time-Dependent Traveling Salesman Problem with Multiple Time Windows (TDTSPMTW) with the following features: 1) each customer can define multiple time windows during which he is available and can be serviced; 2) the service time can be different in every time window of the customer; 3) the travel time depends on the traffic time zone, in which the transit actually occurs; 4) starting and ending depots are treated as customers so that they also have time windows.
The TDTSPMTW problem can be defined as follows.

A. Problem notation
Let I = {1, . . ., n} be the set of customers i ∈ I that are visited by the traveling salesman.Let W i , be the set of i-th customer time windows j ∈ W i , during which the visit can take place.The set of time windows of all the customers will be denoted by W = i∈I W i .Thus, [a j i , b j i ] will denote the j-th time window of i-th customer, where a j i is the beginning, and b j i is the end of the time window, and the service time of i-th customer in j-th time window will be denoted by s j i .Let Z, k ∈ Z, be the set of traffic time zones [p k , q k ], where p k is the beginning, and q k is the end of the traffic time zone, and t k i be the travel time from i-th customer to the next one in the sequence, in k-th time zone.Notice, that we deliberately define travel to the next customer instead of to the current one -this significantly simplifies the considerations.The visit at i-th customer will be denoted by [α i , β i ], where α i is the time of arrival at the customer, and β i is the departure time of the visit.

B. Traffic time zones
Before we can proceed with explaining the minimum route duration algorithm, the problem of traffic time zones has to be accommodated.Since the customers can already have multiple time windows, we can take advantage of this property and create additional, "virtual" time windows so that the travel time in each window is well-defined.For every time window we have to check whether it lies in one traffic time zone, or maybe spans across multiple zones.In the latter case, the original time window has to be divided into smaller, overlapping windows.We shall explain this on the following example.Let [a, b] be a time window with service time s, that spreads over two time zones: [p 1 , q 1 ] with travel time t 1 and [p 2 , q 2 ] with travel time t 2 , such that p 1 < a < q 1 = p 2 < b < q 2 .Than, we need to divide the original window into the following ones: [a, q 1 ] with service time s and travel time t 1 and [p 2 − s, b] with service time s and travel time t 2 .Notice, that the beginning of the second window is brought forward by the service time, because we consider travel time to the next customer in sequence, and the transit starts as soon as the current customer has been serviced.Thus, let V = i∈I V i be the set of virtual time windows of customers i ∈ I. Finally, we can define a function γ : V → Z that maps the virtual time windows into the time zones.Since each virtual time window falls within one traffic time zone, we know that the function is well-defined.

C. The normalized formulation
Having the travel time unambiguously defined for every time window we can normalize the model similarly to [8], [16].Thus, we merge the service time and the travel time into one parameter d j i = s j i + t γ(j) i , denoting the visit duration.At the same time, we postpone the ending of every (virtual) time window by the travel time associated with it, i.e. windows ]. Notice, that from this point on in the article the departure time will have new meaning, i.e., the moment the salesman reaches the next customer in the sequence.
The master problem of TDTSPMTW, i.e. the problem of finding optimal sequence of customers to be visited during one day, can be defined as follows: The problem of finding the minimum route duration for a given sequence of customers (i.e. the subproblem of TDT-SPMTW) can be formulated as follows.

E. Subproblem formulation
Let us define the main decision variables and explain their meaning.For the time windows selection we define: 1 if the visit at customer i takes place within the time window j 0 else (2) Using these notations, we write a mathematical model for the TDTSPMTW subproblem.The objective is to minimize the route duration (calculated as a difference between departure time of last customer and arrival time at first customer): j∈Vi Constraints (4) handle the time windows in a classical way (the departure time must be within a time window) with the noticeable addition of the upper index, since we have multiple time windows.Constraints (5) ensure that exactly one time window per customer is chosen.The auxiliary constraints (6) compute the start of every visit (arrival time α).Finally, constraints (7) forbid to start the next visit before the current customer departure time, so that the visits do not overlap.The y variables are binary, and departure times β are non-negative real variables (8).

III. MINIMUM ROUTE DURATION ALGORITHM
The minimum route duration algorithm that we have developed requires a feasible solutions to start with.Hence, the preprocessing of the virtual time windows is needed.

A. Preprocessing
In order for the algorithm to work, unnecessary time windows from the bottom-right (see Algorithm 1) and top-left (see Algorithm 2) corners have to be removed (if you imagine the first window of the first customer in the bottom-left corner of a diagram, and the last window of the last customer in the topright corner).Similar approach has already been proposed in [16], but in our formulation the visit duration may be different in subsequent time windows of a customer, hence we can dismiss only the ones that lead to unfeasible schedule.
Algorithm 1 TDTSPMTW: top-down preprocessing return false {schedule not feasible} 6: for j ∈ V i do 8: end for 10: α ← max j∈Vi {b j i − d j i } 11: end for 12: return true Algorithm 2 TDTSPMTW: bottom-up preprocessing return false {schedule not feasible} 6: for j ∈ V i do 8: end for The top-down preprocessing algorithm (Algorithm 1) begins with calculating the latest possible arriving at the last customer (denoted by α).Than, starting from the second to last customer, virtual time windows are removed, for which the visit starting at the beginning of the window would exceed α (line 3).If after this step there is no more time windows, the schedule (sequence of customers) is not feasible, and the algorithm terminates.Otherwise, the remaining windows are tightened so that the ending of each window does not exceed α (line 8).Finally, the α value is recalculated for the current customer (line 10) and the procedure starts over with previous customer.
The bottom-up preprocessing algorithm (Algorithm 2) works almost identically as top-down.The differences are as follows: • the earliest possible departure time from the customer (denoted by β) is taken into consideration (lines: 1, 10); • algorithm iterates from the second customer to the last one; • windows are removed if the beginning of visit that starts as late as possible overlaps the departure time β (line 3); • the remaining windows are tightened so that they do not begin earlier than β (line 8).The minimum route duration algorithm that we propose in this paper consists of iteratively reviewing schedules of which the one with the shortest duration is chosen.The procedure of constructing each schedule is divided into two phases.

B. Phase 1: constructing feasible sub-schedule
The first phase consist of computing a feasible sub-schedule that starts from a given customer, and ends with the last one (see Algorithm 3).This procedure requires an index of for j ∈ V i do 5: end if end for 12: end for 13: end if 14: return τ * a customer and an initial departure time.Repeatedly, among the virtual time windows that have enough room to fit the visit after the given departure time (b j i − d j i ≥ β), the earliest visit departure time (calculated as: max{a j i , β} + d j i ) is searched for.

C. Phase 2: constructing dominant sub-schedule
This phase is based on the notion of the dominant solutions (see [8], [14], [16]).To put it simply, a schedule with starting time α 1 dominates schedule with starting time α 2 , if α 1 > α 2 and at the same time ending time β 1 = β 2 (cf.[8]).In our procedure, instead of fixed visit departure time, we use the arrival time (obviously, α = β − d).Nevertheless, the principle is the same -we repeatedly search for the latest possible starting time of a visit (calculated as: min{b j i , α}−d j i ) among the virtual time windows, that are suitable (a j i + d j i ≤ α), i.e. windows in which the currently considered visit does not overlap the initial one.

D. The main algorithm
The algorithm enumerates schedules, constructing them in a particular fashion.For each virtual time window, first, a feasible sub-schedule is constructed with the initial departure time set to the end of the current time window (line 5).Secondly, a dominant sub-schedule is constructed with the initial arrival time set to the latest possible start of the visit in the considered window (line 6).The route duration is than computed as the difference between the departure time from last customer and arrival time at the first customer (line 7).The best solution found during the process is stored (lines 8-10) and returned (line 14).For the algorithm overview see Algorithm 5.
Although the main procedure of the algorithm looks selfexplanatory, the reason it finds the optimal solution is not trivial.Savelsbergh [14] has introduced the concept of forward time slack to postpone the beginning of service at a given customer.It can be proven, that the optimal schedule can be postponed until one of the visits ends with the time window.Hence, by iteratively reviewing schedules one by one so that every possible time window with visit at the end of it is taken into consideration, an optimal schedule is found by our algorithm.

E. Computational complexity
The preprocessing procedures have both O(|V|) time complexity.The sub-schedule construction procedures have together O(|V|) time complexity (they consider disjoint sets of time windows).The main procedure has O(|V|) time complexity (every time window is taken into consideration).Hence, the total time complexity of the algorithm is O(|V| 2 ).

IV. TDTSPMTW MASTER PROBLEM ALGORITHMS
For solving the TDTSPMTW master problem, i.e. computing the one-day schedule, we have chosen three different

A. Simulated annealing
Simulated annealing was first introduced by Kirkpatrick [10], while Černý [1] pointed out the analogy between the annealing process of solids and solving combinatorial problems.Researchers have been studying the application of the SA algorithm in various fields of optimization.Koulamas [11] presented a survey of operational research problems in which the heuristic was applied.The effectives of the algorithm was also inspected in particular by Hurkała and Hurkała [5], [6], and also Hurkała and Śliwiński [7].
The optimization process of the simulated annealing algorithm can be described in the following steps.Before the algorithm can start, an initial solution is required.Then, repeatedly, a candidate solution is randomly chosen from the neighborhood of the current solution.If the candidate solution is the same or better than the current one, it is accepted and replaces the current solution.A worse solution than the current one still has a chance to be accepted with, so called, acceptance probability.This probability is a function of difference between objective values of both solutions and depends on a control parameter taken from the thermodynamics, called temperature.The temperature is decreased after a number of iterations, and the process continues as described above.The optimization is stopped either after a maximum number of iterations or when a minimum temperature is reached.The best solution found during the annealing process is considered final.For the algorithm overview see Algorithm 6.

Algorithm 6 Simulated annealing
Require: Initial solution s 1 1: s * ← s 1 2: for i = 1 to N do 3: for t = 1 to N const do 4: if δ ≤ 0 or e −δ/kτ > random(0, 1) then τ ← τ * α 14: end for 15: return s * The main building block of the simulated annealing is the temperature decrease (also known as the cooling process), which consists of decreasing the temperature by a reduce factor.The parameters associated with this mechanism are as follows: 1) Initial temperature.
2) Function of temperature decrease in consecutive iterations.
3) The number of iterations at each temperature (Metropolis equilibrium).4) Minimum temperature at which the algorithm terminates or alternatively the maximum number of iterations as the stopping criterion.Let τ be the temperature and α be the reduce factor.Then the annealing scheme can be represented as the following recursive function: where i is the number of current iteration in which the cooling schedule takes place.Second building block of SA that has to be customized for a particular problem is the acceptance probability function, which determines whether to accept or reject candidate solution that is worse than the current one.The most widely used function is the following: where δ = E(s 2 ) − E(s 1 ) is the difference between the objective value (denoted by E) of the candidate (s 2 ) and the current solution (s 1 ), and k is the Boltzmann constant found by: where δ 0 is an estimated difference between objective values of two solutions, p 0 is the initial value of the acceptance probability and τ 0 is the initial temperature.Notice that we use 74 POSITION PAPERS OF THE FEDCSIS.Ł ÓD Ź, 2015 decimal logarithm rather than natural, which is most widely seen in the literature.

B. List-based threshold accepting
List-based threshold accepting algorithm (LBTA) introduced by Lee [12], [13] is an extent of the threshold accepting meta-heuristic, which belongs to the randomized search class of algorithms.The search trajectory crosses the solution space by moving from one solution to a random neighbor of that solution, and so on.Unlike the greedy local search methods which consist of choosing a better solution from the neighborhood of the current solution until such can be found (hill climbing), the threshold accepting allows choosing a worse candidate solution based on a threshold value.In the general concept of the threshold accepting algorithm it is assumed that a set of decreasing threshold values is given before the computation or an initial threshold value and a decrease schedule is specified.The rate at which the values decrease controls the trade-off between diversification (associated with large threshold values) and intensification (small threshold values) of the search.It is immensely difficult to predict how the algorithm will behave when a certain decrease rate is applied for a given problem without running the actual computation.It is also very common that the algorithm with the same parameters works better for some problem instances and significantly worse for others.These reflections led to the list-based threshold accepting branch of threshold accepting meta-heuristic.
In the list-based threshold accepting approach, instead of a predefined set of values, a list is dynamically created during a presolve phase of the algorithm.The list, which in a way contains knowledge about the search space of the underlying problem, is then used to solve it.
The first phase of the algorithm consists of gathering information about the search space of the problem that is to be solved.From an initial solution a neighbor solution is created using a move function (perturbation operator) chosen at random from a predefined set of functions.If the candidate solution is better than the current one, it is accepted and becomes the current solution.Otherwise, a threshold value is calculated as a relative change between the two solutions: and added to the list, where C(s i ) is the objective function value of the solution s i ∈ S, and S is a set of all feasible solutions.For this formula to work, it is silently assumed that C : S → R + ∪ {0}.This procedure is repeated until the specified size of the list is reached.For the algorithm overview see Algorithm 7.
The second phase of the algorithm is the main optimization routine, in which a solution to the problem is found.The algorithm itself is very similar to that of the previous phase.We start from an initial solution, create new solution from the neighborhood of current one using one of the move function, and compare both solutions.If the candidate solution is better, it becomes the current one.Otherwise a relative Algorithm 7 Creating the list of threshold values Require: Initial solution s 1 , list size S, set of move operators m ∈ M 1: i ← 0 2: while i < N do end if 12: end while 13: return list change is calculated.To this point algorithms in both phases are identical.The difference in the optimization procedure is that we compare the threshold value with the largest value from the list.If the new threshold value is larger, then the new solution is discarded.Otherwise, the new threshold value replaces the value from the list, and the candidate solution is accepted to next iteration.The best solution found during the optimization process is considered final.
The list-based threshold accepting algorithm also incorporates early termination mechanism: after a (specified) number of candidate solutions is subsequently discarded, the optimization is stopped, and the best solution found so far is returned.The optimization procedure of the list-based threshold accepting algorithm is shown in Algorithm 8.
The original LBTA algorithm does not have a solution space independent stopping criterion.If the number of subsequently discarded worse solutions is set too high, the algorithm will run for an unacceptable long time (it has been observed during preliminary tests).Hence, we propose to use additionally a global counter of iterations so that when a limit is reached, the algorithm terminates gracefully.
In the first phase of the list-based threshold accepting algorithm the list is populated with values of relative change between two solutions ∆ ≥ 0. After careful consideration, we believe that including zeros in the list is a misconception.In the actual optimization procedure, i.e. the second phase, the threshold value is computed only if the new solution is worse than the current one, which means that the calculated relative change will always have a positive value (∆ new > 0).The new threshold value is compared with the largest value from the list (T hmax ).Thus, we can distinguish three cases: 1) T hmax = 0: since thresholds are non-negative from definition, in this case the list contains all zero elements and it will not change throughout the whole procedure (T hmax is constant 2) T hmax > 0 and ∆ new < T hmax : the largest (positive) threshold value from the list T hmax is replaced by a smaller (positive) threshold value ∆ new .The number of zero elements in the list remains the same throughout the whole procedure and therefore is completely irrelevant to the optimization process.The effective list size is equal to the number of positive elements.3) T hmax > 0 and ∆ new ≥ T hmax : the new solution is discarded and the list remains unchanged.
The main idea behind the list is to control the diversification and intensification of the search process.In the early stage of the search the algorithm should allow to cover as much solution space as possible, which means that the thresholds in the list are expected to be large enough to make that happen.In the middle stage, the algorithm should slowly stop fostering the diversification and begin to foster the intensification of the search.In the end stage, the intensification should be the strongest, i.e. the list is supposed to contain smaller and smaller threshold values, which induces discarding of worse solution candidates.As a consequence, the algorithm is converging to a local or possibly even a global optimum.

C. Variable neighborhood search
Variable neighborhood search (VNS) is a metaheuristic algorithm proposed by Mledanović [9].This global optimization method is based on an idea of systematically changing the neighborhood in the descent to local minima and in the escape from the valleys which contain them.It has already been successfully used in different Vehicle Routing Problems.
The VNS algorithm consists of two building blocks: variable neighborhood descent (VND) and reduced variable neighborhood search (RVNS).
The optimization process of VND can be explained as follows.First, an initial solution is required.Within a given neighborhood a candidate solution is repeatedly generated and it replaces the current one if it is better.After a specified number of iterations (neighborhood size) the neighborhood is changed to the first one if a better than initial solution has been found.Furthermore, the best solution found so far replaces the initial solution.Otherwise, if the search resulted in no better solutions that the initial one, the current neighborhood is changed to the next one.Either way, the whole operation is ). Comparing a positive threshold value ∆ new against zero yields in discarding the candidate solution.The conclusions are as follows: a) it does not matter how many zeros are in the list, Algorithm 8 LBTA optimization procedure Require: Initial solution s 1 , thresholds list L, set of move operators m ∈ M 1: i ← 0 2: s * ← s 1 3: while i ≤ N do