A Comparison of Algorithms for Finding an Efficient Theme Park Tour

The problem of efficiently touring a theme park so as tominimize the amount of time spent in queues is an instance of the Traveling Salesman Problem with Time-Dependent Service Times (TSP-TS). In this paper, we present a mixed-integer linear programming formulation of the TSP-TS and describe a branch-and-cut algorithm based on this model. In addition, we develop a lower bound for the TSP-TS and describe two metaheuristic approaches for obtaining good quality solutions: a genetic algorithm and a tabu search algorithm. Using test instances motivated by actual theme park data, we conduct a computational study to compare the effectiveness of our algorithms.


Introduction
Theme parks like Walt Disney World in Orlando, Florida, attract millions of tourists each year.While these parks provide great entertainment, a common complaint is the amount of time spent waiting in line for the various attractions.Whole industries have arisen to help tourists maximize their entertainment value by offering advice on how to optimally tour parks to minimize such waiting (e.g., [1,2]).Here, we look at the problem of finding a shortest route through an amusement park that takes into account the wait time of the attractions.This is an instance of the Traveling Salesman Problem with Time-Dependent Service Times (TSP-TS) [3].
The TSP-TS is a variation of the Traveling Salesman Problem (TSP) and was introduced by Tas et al. [3].Given a set of locations and distances between every pair of locations, the classical TSP seeks to find the shortest possible route that visits each location and returns to the starting location.For the TSP-TS, the duration of time spent at each location is defined as a function of the arrival time to that location.The objective is to minimize the total route duration, which consists of the sum of the total travel time and the total service time.In our setting, the service times represent attraction wait times and ride times.
Time-dependency in the TSP literature is typically addressed in terms of travel times.In particular, the Time-Dependent Traveling Salesman Problem (TDTSP) seeks to find the shortest tour through the locations when the time to travel depends not only on the distance but also on the time of day the route is traversed.Because the TDTSP incorporates a more realistic touring cost structure, it has been used to model several other applications including scheduling single-machine jobs with time-dependent setup costs [4,5], creating timetables for university exams to minimize backto-back exams [6], production planning for car assembly lines [7], satisfying product demands at minimum travel and purchasing costs [8], and vehicle routing with varying travel times such as within regions of congestion [9][10][11][12][13][14][15].Vander Wiel and Sahinidis [16,17] note that the TDTSP is NPhard, but little has been published in terms of heuristics, especially heuristics that incorporate a time-dependency.Multiple authors [4,5,12,[16][17][18] have proposed mixedinteger linear programs (MILP) for producing solutions to the TDTSP, but these exact solutions have generally been for (4) We create a new metric, , to guide tour construction and augmentation and we develop both a genetic algorithm and tabu search algorithm that can be used to find good quality solutions to the TSP-TS in an efficient manner.
(5) Motivated by actual theme park data, we introduce a collection of test problems that utilize different types of wait time distributions.
(6) Finally, we perform a detailed computational study to compare the effectiveness of our algorithms.
The remainder of this paper is organized as follows.In Section 2, we formally describe our version of the TSP-TS, propose a mixed-integer linear programming (MILP) formulation, and describe a branch-and-cut algorithm.In Section 3, we show that our version of the TSP-TS cannot be modeled as a TDTSP and describe a method for computing a lower bound on the optimal solution.Section 4 describes two different metaheuristic approaches for solving the TSP-TS using a genetic algorithm and a tabu search algorithm.We present our test data and computational results in Sections 5 and 6, respectively.Section 7 highlights our conclusion.

Problem Description and Formulation
In this section, we begin by introducing the notation that will be used throughout the paper.We then present a MILP model for our version of the TSP-TS and describe a branch-and-cut algorithm based on this formulation.

Notation.
Let  = (, ) be a connected digraph with node set  = {1, 2, . . ., } and arc set  = {(, ) | ,  ∈ ,  ̸ = }.Associated with each arc is a travel time from node  to node ,   , while associated with each node  is a service time function   (), where  corresponds to the arrival time to node .For our application, the nodes represent theme park attractions,   represents the time to walk from attraction  to attraction , and   () represents the sum of the wait time and ride time when arriving at attraction  at time .The TSP-TS aims to minimize the total tour duration, beginning and ending at node 1 (possibly representing the entrance to the park), including the total walking time and the total service time.
Tas et al. [3] considered both linear and quadratic service time functions   ().However, such functions are not realistic in a theme park scenario.In practice, service (wait) time data at each attraction is collected at  + 1 discrete times  0 ,  1 , . . .,   .In this paper, we consider two types of   () functions based on data collected at these times.The first is to define   () as a step function where   () is set to be the wait time recorded at attraction  at time  −1 for all  ∈ [ −1 ,   ).The second is to compute   () using linear interpolation of the wait times recorded at  −1 and   for all  ∈ [ −1 ,   ).While wait times computed using linear interpolation are more realistic, they are difficult to model using mixed-integer programming.However, such functions can easily be handled using metaheuristics.

Mathematical Model.
To formulate the TSP-TS as an MILP, we will assume that   () is a step function over the  intervals [ 0 ,  1 ), [ 1 ,  2 ), . . ., [ −1 ,   ].We refer to [ −1 ,   ) as time interval .Our model is formulated on an expanded graph where each arc (, ) from node  to node  is replaced by  parallel arcs from  to , one for each time interval.In this new network we define    to represent the service time (sum of the wait time and ride time) when arriving at node  during time period .
For notational convenience, we introduce the set  = {2, 3, . . ., } to represent the set of possible successor nodes and define  = {1, . . ., } to be the set of all time intervals.The decision variables for our formulation are as follows.
Using the models presented in [12,28] as a starting point, we formulate our version of the TSP-TS as follows. minimize subject to binary ∀,  ∈ ,  ̸ = ,  ∈  (10) Recall that we assume that our tour begins and ends at node 1.
The objective function (2) minimizes the arrival time back to node 1 as only one variable  1 can be positive.Constraints (3) and (4) ensure that each node is visited exactly once.Constraint (5) computes the arrival time to node  after leaving node 1, while constraints (6) compute the arrival time to node  visited after node  ∈ .Note that constraints (5) and ( 6) act as the subtour elimination constraints.The temporal constraints (7) and (8) ensure that the correct time interval  is chosen when visiting node  after node  (note that [ −1 ,   ) = [ −1 ,   − 1] since we assume   is an integer for all ).Finally, constraints (9) ensure that   > 0 only if one of    = 1 for any .To strengthen the continuous relaxation of the formulation, we can include additional restrictions designed to tighten the model as suggested in [12].Toward this end, we define the following sets: We can then add in the following new restrictions, where  is a lower bound on the optimal solution.
Constraints ( 12) and ( 13) tighten the formulation as follows.Suppose    = 1 so that we travel from node  to node  and arrive during time period .Then constraints (12) ensure that if you then visit node  immediately after node , then you must use a time interval that occurs after your arrival to .Suppose    = 1 so that you travel from node  to node  and arrive during time period .Then constraints (13) ensure that when you traveled from node  to node  you must use a time period that finishes service at  before you leave node  to travel to node .Finally, constraint (14) allows us to utilize a lower bound on the optimal solution, such as that computed in Section 3.

Branch-and-Cut Algorithm.
In this section we describe a branch-and-cut algorithm that we developed based on the MILP described in the previous section.Our algorithm begins with a preprocessing phase that first computes a lower bound  to use in constraint (14) using the method described in Section 3. We then utilize the metaheuristics of Section 4 to obtain an incumbent solution to seed the branch-and-bound procedure.The solution obtained from the metaheuristics is also used to determine an upper bound on the number of time intervals  that are needed within the MILP formulation.
Recall that between each pair of nodes there are  parallel temporal arcs, each of which is represented by a variable    .It is therefore advantageous to keep the number of temporal arcs to a minimum so as to reduce the size of the MILP formulation, which we can accomplish using the best tour found by the metaheuristics.
After the preprocessing phase, we submit the MILP (2), (3), ( 4), ( 5), ( 6), ( 7), ( 8), ( 9), (10), (12), (13), and ( 14) and the incumbent solution to the mixed-integer programming solver CPLEX, which was used to implement the branchand-cut algorithm.We utilize CPLEX's callback functionality to add cuts during the enumeration.Specifically, we incorporate additional subtour elimination constraints.Recall that constraints ( 5) and ( 6) act as the subtour elimination constraints for our formulation, which are similar to the Miller-Tucker-Zemlin constraints for the TSP [30].Note that these constraints are weaker than the typical, yet exponential in number, subtour elimination constraints for the TSP that were introduced by Dantiz, Fulkerson, and Johnson [31].We can extend these tighter subtour elimination constraints to our formulation as follows, where  is the node set of a subtour and || is the cardinality of : At each node of the branch-and-bound tree we identify violated subtour elimination constraints (15) of the LP relaxation using a separation technique based on a min-cut algorithm of [32] and subsequently add these constraints to the formulation.Even though these additional cuts are not necessary for the elimination of the subtours, the inequalities strengthen the linear programming relaxation of the problem.Note that the inequalities (15) added at any node of the enumeration tree are valid for all the other nodes because these inequalities are valid for the entire formulation (3), (4), ( 5), ( 6), ( 7), ( 8), ( 9), ( 10), ( 12), (13), and (14).Thus, at a given node, all the inequalities generated so far were incorporated into the formulation.

Computing a Lower Bound
We begin this section by showing that our version of the TSP-TS cannot be modeled as a TDTSP, and therefore we are unable to use the techniques that have been developed for computing lower bounds for the TDTSP.We then describe a method for computing a lower bound for the TSP-TS.First, we need to introduce some additional notation.We define each tour of the nodes as a permutation  = [1,  2 , . . .,   ], where   ∈ , 2 ≤  ≤ , and associated with tour  is a tour cost, () =    1 , representing the time it takes to traverse the tour as described by permutation  and to return to  1 = 1, the entrance.To break down this tour cost, recall we defined   −1   to be the time that the tour arrives at the th attraction, and let us define    to be the time that the tour departs the th attraction.Note we define  1 = 0. We can then recursively define the arrival and departure times for each of the attractions as follows.The arrival time at the th attraction for  > 1 is The departure time of the th attraction for  > 1 is Thus, Following the approach in [19] it seems reasonable to try to incorporate service time variability into travel times and use many of the results in the TDTSP literature to establish bounds for the TSP-TS.Much of this literature proposes to model the TDTSP by establishing a nonnegative traversal rate V  () for arc (, ) during times  in time period [ −1 ,   ).This rate is fixed during each time period but is allowed to change to V  (  ) when the time period changes from [ −1 ,   ) to [  ,  +1 ), where   ∈ [  ,  +1 ), even if the arc traversal is not complete.This modeling technique offers the benefits of establishing lower bounds for the TDTSP.However, as we show here, the TSP-TS cannot be modeled by using this technique.
Consider an instance of the TSP-TS with the following conditions during time periods [ −1 ,   ) and [  ,  +1 ): (1) It is possible to arrive at node  −1 at some time  0 , where  0 <   and arrive at node   at some time (2) The slope of the service time function for some node  −1 is -1 during time period [  ,  +1 ) so that, for instance, In this instance of the TSP-TS, if some tour traverses arc ( −1 ,   ) starting at time   , then it will arrive at node   at some time   =   +   −1 (  ) +   −1   .However, if some tour traverses arc ( −1 ,   ) beginning at  +1 , then it must arrive at node   at time Note in departing node  −1 at time   or at time  +1 we arrive at node   at time   in both cases.
Let us attempt to model this TSP-TS using a nonnegative traversal rate V  −1   () common in modeling the TDTSP.Since we arrive at node   at the same time,   , regardless of whether we begin traversing arc ( −1 ,   ) at time   or time  +1 , we cover the same distance, and we have This implies that V  −1   (  )[ +1 −   ] = 0, and thus V  −1   (  ) = 0. Hence, a tour in our TSP-TS cannot cover any ground along arc ( −1 ,   ) during the interval [  ,  +1 ).This means that if some tour in the TSP-TS finishes traversing arc ( −1 ,   ) at time , then  ≤   or  >  +1 .This contradicts the first condition of our TSP-TS instance, which implies that a tour can complete edge ( −1 ,   ) at time between   and  +1 .Thus, this instance of a TSP-TS cannot be modeled by converting service times to edge traversal rates as is common in modeling the TDTSP.Further, this scenario is common to amusement parks, especially among attractions that have batch servicing.
Cordeau et al. [19] use variable arc traversal speeds to find a lower bound for the TDTSP.They create this lower bound by assigning, in each time interval, each arc its maximum speed found over all time intervals.This creates an asymmetric Traveling Salesman Problem (ATSP) as all arcs now have a fixed traversal time.Solutions to this ATSP provide a lower bound to the TDTSP.However, as shown in [3], solutions to the TSP in the case of variable service times are not necessarily a lower bound for the TSP-TS.
For a preliminary lower bound on the TSP-TS, one may assign to each node  a wait time equal to the minimum wait time achieved at that node throughout the day,   = min    ().If we incorporate this time into every arc leaving node , the problem becomes an instance of the asymmetric TSP where arc (, ) has a weight of   +   .Let  0 () denote the time for a tour  to finish in this ATSP.Notice that, for any tour , () ≥  0 (), since every arc in the TSP-TS costs at least as much as the corresponding arc in the ATSP.Thus, for a tour  *  that is optimal to the ATSP problem,  0 () ≥  0 ( *  ), and as a result () ≥  0 ( *  ), for all tours .Thus  0 ( *  ) is a lower bound for the TSP-TS.
However, in amusement parks, wait times for different rides often increase together, and it is unlikely to create a tour with every ride achieving a minimum wait time.We can find an improved lower bound by incorporating the extra wait cost incurred by not visiting each attraction at its minimum wait time.Define    =    (  −1   ) −    to be the wait time above    for each node   in a tour .Then  has some total extra cost accumulated over all nodes, ∑  =1    .Consider a tour  * that is optimal to an instance of a TSP-TS.We would like to find a tighter lower bound of ( * ) than  0 ( *  ).To that end, let  = { * 1 ,  * 2 , . . .,  *  } be the set of nodes corresponding the first  attractions visited by the optimal tour  * .Let Π() be the set of all permutations of the nodes in , and let   = [  1 ,   2 , . . .,    ] be the permutation of  that achieves min ∈Π() ∑  =1    .Each permutation  ∈ Π() has a time that it requires to visit each of its  nodes.We call this time ().Let  * ∈ Π() be the permutation of  that achieves min ∈Π() ().That is, among all tours that begin with some permutation of the nodes in ,  * finishes the tour of the nodes in  fastest, while   minimizes the extra wait incurred in any tour of nodes in .For each node V ∈  −  we can redefine  V , call it   V (( * )), to be the minimum value of  V () over all times ( * ) ≤  ≤ (), where  is any solution to the TSP-TS.Thus node V ∈  −  incurs an extra wait of at

least 𝜇 󸀠
V (( * )) −  V because it must be visited after time ( * ), as shown in Figure 1.
To build a new lower bound on ( * ), we start with the observation that for the particular  used by  * , the extra wait incurred by the first  nodes of  * is bounded below by the minimum possible extra wait over all permutations of the nodes in : Additionally, the extra wait time for all nodes visited in  * after the th node must be at least as large as the difference of the updated minimum wait after ( * ) and the original minimum wait time: Starting with the idea of adding extra wait times incurred to the existing lower bound  0 ( *  ) and incorporating inequalities ( 21) and ( 22), we observe the following: Define the lower bound of the extra wait time of the nodes in  as Since we do not know which  nodes of  the optimal tour  * starts with, we exhaustively consider each of the (   ) possible sets  that  * may start with and examine all ! possible orderings of each  to find   ,  * , and , the set of  nodes that minimizes ().This would then form a lower bound   that improves upon the previous lower bound,  0 ( *  ).Continuing from inequality (23), we have As  * must start with some  nodes, it must incur extra wait based on some  that was considered, and so must have extra wait at least as large as the minimum, ().This extra wait was completely ignored in the ATSP.In using this approach to form the lower bound, we typically use 5 ≤  ≤ 7.

Metaheuristics
In this section we describe two different metaheuristics that can be used to solve our version of the TSP-TS.In particular, we develop both a genetic algorithm and a tabu search algorithm.
Several authors have proposed metaheuristic approaches to solving instances of the TDTSP.Testa et al. [26] implement eight different genetic operators on 50 randomly generated instances of the TDTSP and found certain combinations of genetic operators were effective at producing high-quality solutions.They show that the crossover operators edge recombination and cycle crossover along with high levels of mutation produce particularly good solutions on these instances and better solutions than can be generated via dynamic programming in a fraction of the time (see Section 4.1.4for crossover descriptions).Li et al. [25] use a chained Lin-Kernighan algorithm with a double-bridge kick mutator within a genetic algorithm framework to effectively produce solutions to problems with realistic traffic assumptions (including time periods and areas with traffic jam travel restrictions).
Before we discuss our metaheuristic approaches, we define a new metric for comparing possible adjustments to a tour that takes into account differences from the average wait times of attractions.Certain tour-construction heuristic approaches will attempt to build a tour making choices using local information at each step.Two nodes  and  with equal service times   () =   () can seem equally attractive to visit in time .However, a twenty minute wait for an attraction that has an average wait time of an hour at time  is more of a bargain than an attraction with an average wait of fifteen minutes.We let   be the average wait for attraction .Then for each time , a node  has a deviation from its average Nodes that have little variation from average can be visited during any time period without much bonus or penalty.However, nodes that show large deviation from average can potentially provide huge tour time savings and also carry enormous time penalties.For example, Figure 2 shows an attraction whose wait time is broken down into 15 minute time periods.The average wait time, approximately 40 minutes, is shown with the dashed line.We note that   ( node can be scheduled earlier or later in the tour, there are time-saving benefits (when   < 0), and there are penalties for scheduling it during the middle of the day (when   > 0).We also take into account an average travel time between nodes.Let   =   −   , where   is the average of all  − 1 travel times associated with edges starting at node  over all .Using Equations ( 16), (17), and ( 26) we define a metric that incorporates these averages of both wait times and travel times: To the best of our knowledge, this metric has not appeared in the literature.
4.1.Genetic Algorithms for the TSP-TS.Genetic algorithms can be customized to fit a problem and are normally defined by choices in the following components [33]: population representation and initialization, fitness evaluation, reproduction selection, and choice of genetic and replacement operators.
We define each candidate solution in the population to be represented by a permutation  of the  locations, which we refer to as a tour, and having fitness (), the tour duration given in (18).We keep 40 candidate tours at every generation.At each iteration of a genetic algorithm, we choose an operator via dynamic operator selection as described in [26].With this selection process, operators that have demonstrated success have a higher likelihood to be chosen to be used, where success is measured by how many direct children, grandchildren, etc. were inserted into the top half of the population as a result of this operator.Either one or two tours from the population are chosen at random, depending on the number of tours the operator takes as input.The operator produces an offspring tour.As in [26], we employ a ( + 1) reproduction approach.That is, if  is the candidate solution with the highest (worst) fitness value, and  is the new offspring tour, we compare the tour costs of  and  and keep the tour with the lower (better) fitness value.

Population Initialization.
In generating the initial population, we randomly generate a portion of the 40 tours.We combine these random tours with tours generated by more sophisticated algorithms.We tested many of the most commonly used construction heuristics, such as nearest neighbor and a dynamic programming algorithm proposed by Malandraki and Dial [21].We also tested a variant of nearest neighbor introduced here, called -Nearest Neighbor (NN).This heuristic is the same as nearest neighbor except at each iteration of the tour construction, if we are located at node  at time period , we choose a node  that minimizes the value of   () in (27).

Local Search Operators. Many genetic algorithms employ a local search heuristic to improve population tours.
It has been shown in [34] to be an effective way to improve solution quality.A common local search heuristic for the TSP is the -opt exchange heuristic.These heuristics remove  arcs of a current solution and replace them with  different arcs to reproduce a tour.This procedure is repeated until no improvements can be found.The 2-opt heuristic has been shown to be particularly effective for improving solutions for the TSP. Figure 3 shows a 2-opt improvement for a tour of eight attractions.
The 2-opt in Figure 3 reverses arcs (, ), (, ), and (, ) meaning that these attractions are visited at different times of day than previously, which can significantly change the tour length.In a time-independent case, this would not cause the same effect.
The Lin-Kernighan algorithm [35] is a variable -opt procedure that dynamically determines how many exchanges will be considered and has been shown to be very effective at producing high-quality solutions for both the symmetric TSP and the TDTSP.We considered three versions of the Lin-Kernighan (LK) heuristic differentiated by whether it took the first improvement, exhaustively searched to find the best 2opt improvement, or used the  metric.We found that the version that uses the first improving 2-opt, we refer to this as LKTD-1.

Mutators. Mutators take one tour as input and produce
a new tour by augmenting the selected tour in some fashion.We considered several mutators and our tests show the following to be effective for this problem.
(i) Uniform Order-based Mutation (UOM): This operator, described in [36], works on a single parent tour by performing a simple swap of two randomly selected nodes in that tour.
(ii) -Mutate (M): This mutation operator, like UOM, performs a single swap of nodes in a parent tour by finding the node after the arc that has the greatest penalty in terms of the  metric from Equation (27) and swaps it with a random other node.(iii) Double-Bridge Kick Operator (BK): As described in Li et al. [25], this operator, which cannot be achieved with 2-opts alone and which is useful in escaping local minima, randomly removes four arcs from the tour and relinks the nodes in a different manner, as shown in Figure 4.

Crossover Operators.
Crossover operators take two parent tours from the population and try to strategically combine them to produce an offspring.The crossover operators we found most effective were the following.
(i) Cycle Crossover (CC): This operator, described in [37], produces an offspring from two parents by ensuring any node appearing in position  in the offspring tour must appear in position  in at least one of the parent tours.
(ii) Edge Recombination (ER): First proposed for the TSP in [38], ER produces a single offspring tour from two parent tours by building a list of edges present in both parents and transferring these edges to the offspring.The offspring is then completed by adding in the remaining nodes offering preference to nodes with smaller numbers of outgoing edges.As noted in [26], the motivation behind ER is that in the TSP connections between locations contribute more to the tour length than the position of the locations in the tour.
(iii) Node Recombination (NR) via Path Relinking: This operator follows the strategy of path relinking in [39].A path is constructed between two parent tours by searching the neighborhood of the tours.Given two parent tours,  and , we search solutions that share nodes in the same position in the tour as both parents.To describe the process, Figure 5 shows two parent tours,  = [, , , , , , , ] and  = [, , , , , , , ], which share common nodes in positions 2, 3, and 7.In positions where  and  do not have common elements, these elements are changed in succession as follows.If  differs from  in position ,   is moved to its position in , say position , the element in   is moved to position  while all other nodes are held constant.We thus create a tour   with more nodes in common with  than  had in common with .For example, in Figure 5,  differs from  in position 1.The element in position 1 in tour , node , is moved to its position in , position 8, and the element in position 8 in , node , is moved to position 1 in  to create a new tour.This change is performed sequentially for each 1 ≤  ≤  where the tours have noncommon elements and each of these resulting tours is evaluated for its fitness.The tour with the best fitness after these switches is chosen as the new parent tour , tour [, , , , , , , ] in Figure 5, and the process is repeated until any switch results in tour .The best tour found in the intermediate steps is chosen as the offspring of  and .
(iv) Sort Crossover (SX): This operator is essentially NN where the selection set is restricted to the two parent tours.The SX operator takes two tours  and  and constructs the offspring tour  so that at each iteration ,   is chosen to be the node with the minimum  metric (equation ( 27)) out of the next node in  or in  that has not already appeared in the constructed tour .

Preliminary Computational Study.
To test the effectiveness of different local search operators, mutators, and crossover operators when applied to instances of the TSP-TS, we completed an extensive computational study on randomly generated instances of the TSP-TS (see Section 5 for details of our test problems).For the sake of brevity, we omit the specific details of our tests here.The most effective combination of initialization, local search operators, mutators, and crossover mutators is described as Algorithm A in Table 1.To provide a basis of comparison we also describe the genetic algorithm proposed in [26], identified as Algorithm B in Table 1.Note that all of our tests regarding genetic algorithms in Section 6 will focus on these two algorithms.

Tabu Search Algorithm.
In this section we describe the tabu search algorithm used to solve TSP-TS.Tabu search, originally introduced by Glover [40] and later formalized in [41][42][43], is a metaheuristic algorithm that is known to be quite effective for hard combinatorial optimization problems.Tabu search begins with an initial solution and uses a local search procedure to move from the current solution to the best one in its neighborhood, even if the move leads to a worsening of the objective function value.To prevent becoming stuck at a local optima plateau, attributes of solutions that have recently been visited are declared "tabu" for a certain number of iterations.

Search Neighborhoods.
In our tabu search, which is similar to [44], we use the nearest neighbor heuristic to generate an initial solution.Our search neighborhoods are defined by performing all possible swap and shift moves that are not tabu.Given a tour  = [ 1 ,  2 , . . .,   ], where   represents the attraction in position , the swap move chooses two attractions   and   and then exchanges the attractions so that each attraction is located in the position previously occupied by the other one.The shift move chooses two attractions  and  with  <  and a number  with  ≤ , which defines the number of attractions to move.It then relocates attractions  −  through  to , shifting any attractions between  and  to the left.In our algorithm, we consider all possible shift moves for every pair of positions  and  and for every value  ≤ .
The two neighborhoods defined by the swap and shift moves are determined separately and the best allowable moves from each neighborhood are compared to choose the one to be performed.During the search process, we use the standard aspiration criteria which allows a tabu move if it leads to a tour with a better solution than the best solution found thus far.Our tabu search stops after a specified amount of time has passed, which we will discuss further in Section 6.

Tabu Status and Tenures.
For our tabu search, we used a fixed tabu tenure of 7.5 ln() as suggested in [44].When a swap occurs, exchanging attractions  and , swapping any attractions located at the positions  and  is declared tabu.That is, swapping the index positions becomes tabu rather than swapping the attractions themselves.When a shift occurs, shifting  and  with  < , moving attractions  and  to their previous positions is made tabu.This is similar to how [44] handled making swaps tabu, but we found it useful when performing shifts.A shift move of  and  is considered tabu if moving attraction  to the position of  is tabu.Note that the swap and shift tabu lists are independent of each other.

Diversification Strategies.
To help drive the search into new regions of the solution space we use both light and strong diversification techniques.A light diversification is applied when the current tour's fitness value has not changed by over 0.5% in the past ⌊/4⌋ iterations or when the current tour's fitness value has not changed at all in the past five iterations.When the light diversification is applied, the exhaustive swap and shift neighborhoods are determined using an alternate fitness function that adds a small penalty to the true objective value.Specifically, for every tour   , the frequency penalty (  ) is equal to the sum of the iterations that each node  has been at index  in the current tour, divided by the total number of iterations.The new fitness value for the tour is then equal to F(  ) = (  ) + 0.1(  )(  )/√.This new fitness function urges the exhaustive searches to find tours outside of the local optima plateau.This type of function was first suggested by [45] and later used by [44].
A strong diversification is utilized when the current tour's fitness value has not changed by over 0.5% in the past ⌊/8⌋ iterations or when the best fitness value has not changed for over 9 iterations.When strong diversification is applied, a subset of ⌊/5⌋ nodes from the current tour are randomly permuted to yield a new current tour.To select the subset of nodes to be scrambled, we keep track of the number of iterations each node has been at the same index and select the ⌊/5⌋ nodes that have been at their current index the longest.This diversification can have a significant impact in terms of sending the algorithm into an unexplored area of the search space.

Test Instances
Within a theme park, the wait time data varies from attraction to attraction and from hour to hour.Using data supplied by touringplans.com[29], Figure 6 shows estimated wait times for 29 attractions at Walt Disney World's Magic Kingdom for the date October 20, 2014.As evidenced from the figure, some rides have high wait times throughout the day while others have wait times that peak early, such as thrill rides, peak midday, such as kid attractions, or exhibit a sawtooth pattern, such as bulk-entry shows.To develop test problems that mimic the types of wait time distributions seen at Walt Disney World, we generate a wait time   () at each attraction  that resembles either a uniform, sawtooth, or peaking distribution.We also generate a travel time   () between attractions  and .Because we focus on the scenario that the wait times,   (), dwarf the walking times,   (), we generated the walking time functions independent of .To generate   , each node was randomly assigned a position on a 1000 by 1000 grid, and the Euclidean distance was calculated between each pair of nodes.This result was divided by 100, to give a result generally between 0 and 15, but largely clustered around the 5-9 range.This aligned well with the travel times seen in real-world scenarios.
We used the actual theme park data to motivate the collection of distributions from which we generate wait times for networks of various sizes.For each distribution, we have a collection of parameters that we sample from the indicated range to give each node its own wait time distribution.There are five types of distributions we consider to generate a distribution on the interval  ∈ [0,   ].In each case, we compute a scaled distribution f() incorporating the average of the distribution over the interval [0,   ], , and a randomly chosen scaling factor, : (1) Uniform distribution: Let () = 1, where the parameter  is randomly chosen from the interval [5,15].

Journal of Applied Mathematics
(3) Bimodal distribution: Let where the following parameters are chosen randomly from the given intervals: (4) Beta Distribution: Let where  is chosen randomly from the set {2, 3, 4, 5} and  is chosen randomly from the interval [15,70].
The Beta and Erlang distributions can produce asymmetric distributions, so we also introduce a randomly chosen reflection option where we use (  − ) instead of ().Figure 7 shows examples of each of these distributions.As mentioned before, the distributions in Figure 7 attempt to mimic wait times that might typically occur in a theme park.The Beta distribution with  = 5 and  = 40 might describe an attraction that has a larger wait time at the beginning of the day diminishes throughout the day.Whereas the bimodal distribution represents an attraction that gets crowded during the early morning hours and during the afternoon but has lower wait times during the lunch hours.The sawtooth pattern might represent an attraction that takes in visitors at fixed points in time, like a show or movie attraction that runs on the hour and half-hour.
Using this collection of distributions, we developed ten different types of networks.These networks are described in

Computational Results
We begin by examining our results for the three metaheuristics: genetic algorithms A and B described in Section 4.1, and the tabu search algorithm described in Section 4.2.The algorithms were implemented in Java and executed on a Dell Precision T5610 Workstation equipped with dual 2.6 GHz processors and 32 GB of RAM running 64-bit Windows 7.
To allow for comparisons between the three algorithms, we utilized a fixed time limit of five minutes for each algorithm, as opposed to using a fixed number of generations for the genetic algorithms or a fixed number of iterations for the tabu search.
In Table 3 we report the minimum objective values (tour fitness () from equation ( 18)) found using the three different algorithms, along with the lower bounds computed using the method described in Section 3, for the ten networks of size 30.Within the table, the first column identifies the network type (as described in Table 2), while the second column gives the lower bound computed.The next three columns give the minimum objective values found by genetic algorithm A, A number of observations can be made from Table 3. First, note that when the service times were computed using linear interpolation, the tabu search algorithm was the most effective, achieving the minimum tour length found among the three algorithms for eight of the ten networks.However, when the service times were computed using a step function, genetic algorithm A was superior, achieving the minimum tour length found for nine out of the ten networks, while the tabu search algorithm performed poorly.Second, note that for Network 2, which consists entirely of wait times that are constant throughout the day, all three algorithms were able to achieve a tour length equal to the lower bound of 388, which must therefore be the optimal solution.This is not surprising in that when all the wait times are constant, the TSP-TS essentially reduces to a standard TSP, which is not difficult to solve for 30 nodes.Third, the quality of the lower bounds obtained varied significantly depending on the network.For example, the lower bounds for Networks 1, 4, 5, and 7 were quite poor in that the gap between the bound and the best tour found is quite large.We believe this is due to the abundance of attractions with Beta and Erlang distributions and the nature in which we formulate the lower bound.Recall, the lower bound tests permutations of 5 ≤  ≤ 7 attractions to find a permutation with low tour costs and combine that with minimum wait times within the remaining time window from the attractions left to tour.In most theme parks, this works well in achieving a good lower bound because most attractions achieve a minimum wait time during the first few time periods and all attractions see an increase in wait times as time period progress.In networks with left and right Erlang and Beta distributions, there are attractions achieving their minimum wait times during the first few time periods and others achieving their minimum in later time periods.Hence, this tends to create lower bound tours utilizing attractions with minimum wait times achieved in the first few time periods as part of the set  but still have many attractions achieving minimum wait times later, thus not increasing the lower bound substantially.
In Table 4 we report the same information as in Table 3, but for the size 50 networks.Genetic algorithm A and the tabu search algorithm performed similarly for both methods of computing the service times, achieving the minimum objective value found for about half of the network instances.Genetic algorithm B appears to be the weakest of the three algorithms, only achieving the minimum objective value found in three of the network instances for both methods of computing the service times.Note that for Network 2, which consists entirely of wait times that are constant gaps.Again, we believe this is due to the presence of one-sided Beta and Erlang distributions in large quantities.

Concluding Remarks
On crowded days at theme parks, the wait time of an attraction depends heavily on the time of day that the attraction is visited.Given that frustration over waiting in lines can ruin vacations, there is a need for effective solutions to the TSP-TS that incorporates the time-dependency of the problem.In this paper, we have extended the TSP-TS literature to include these real-world service time functions and showed how a lower bound on an optimal solution can be calculated.We investigated a variety of approaches to generate optimal and near-optimal solutions including a new metric used in construction approaches that gives preference for using attractions with high potential time savings, a branch-andcut algorithm for our proposed MILP, and metaheuristic approaches.Our genetic algorithm and tabu search produced promising and efficient results, but our computational study highlighted the need for more sophisticated methods to yield tighter lower bounds to aid the branch-and-cut algorithm.These results are encouraging, but there are further constraints that could be included in the model to make it more realistic.As it stands, our tours give patrons no time to eat or take breaks during the day, two highly desirable characteristics of theme park tours.In particular, meals could be incorporated as additional attractions where the wait times could be interpreted as the time it takes to complete a meal at a given establishment.Breaks could be implemented in a similar fashion.Each of these could have user-defined preferences for the times of day that these occur.
Further, theme parks have innovated around the wait time issue.Most parks now offer a limited number of queue priorities, which are passes that allow for shorter wait times where patrons are given a time window to arrive to receive priority entrance to the attraction.Sometimes these passes are unlimited, but patrons are charged extra for their use.In other instances, the number is limited, but the opportunity is included with the price of admission.With the introduction of these priority passes, an interesting question arises for future study.Mainly, how can patrons optimally choose a limited number of queue priorities to further reduce the time it takes to tour a theme park?

Figure 1 :
Figure 1: Extra wait time incurred by calculating a new minimum wait time   V (( * )) for the feasible window of time that node V will be visited, [( * ), ()].

Figure 5 :
Figure 5: Path relinking shown on a tour of eight attractions.Underlined nodes signify agreement between the two original parent tours.Boxed nodes signify nodes that are in the correct spot due to the swap performed during that step.

Table 1 :
[26]tic algorithms used for computational study.Algorithm A is our best combination of operators and Algorithm B is the best found in the literature[26].

Table 2 ,
which provides the percentage of attractions that are assigned to each type of distribution for each network type.

Table 2 :
Percentage of testing networks comprised by each type of distribution.

Table 3 :
Results for size 30 networks.
() are computed using linear interpolation.The last three columns give the same information for when the service times   () are computed using a step function.The boldfaced entries indicate the minimum value achieved over all cases.

Table 6 :
MIP gaps of minimum objective values found.