Trajectory stability in the traveling salesman problem
Sergio Sánchez, Germinal Cocho, Jorge Flores, Carlos Gershenson, Gerardo Iñiguez, Carlos Pineda
TTrajectory stability in the traveling salesmanproblem
Sergio S ´anchez , Germinal Cocho , Jorge Flores , Carlos Gershenson , GerardoI ˜niguez , and Carlos Pineda Instituto de F´ısica, Universidad Nacional Aut ´onoma de M ´exico, 01000 M ´exico D.F., Mexico Instituto de Investigaciones en Matem ´aticas Aplicadas y en Sistemas, Universidad Nacional Aut ´onoma de M ´exico,01000 M ´exico D.F., Mexico Centro de Ciencias de la Complejidad, Universidad Nacional Aut ´onoma de M ´exico, 01000 M ´exico D.F., Mexico SENSEable City Lab, Massachusetts Institute of Technology, 02139 Cambridge, MA, USA ITMO University, 199034 St. Petersburg, Russian Federation Department of Computer Science, Aalto University School of Science, 00076 Aalto, Finland Faculty of Physics, University of Vienna, 1090 Wien, Austria * Corresponding author email: [email protected]
ABSTRACT
Two generalizations of the traveling salesman problem in which sites change their position in time are presented. The way therank of different trajectory lengths changes in time is studied using the rank diversity. We analyze the statistical properties ofrank distributions and rank dynamics and give evidence that the shortest and longest trajectories are more predictable androbust to change, that is, more stable.
Introduction
Imagine that a certain product must be delivered, in the most efficient way, to N stores in a region. When stores are fixed atgiven sites in space, finding the shortest path that links them all is the target of the traveling salesman problem (TSP) , whichhere we refer to as static TSP. But what happens if the product must be delivered to moving targets such as peddlers? TheTSP becomes time dependent. Another variation of the TSP, related to the peddlers’ scenario, occurs when one or more sitesdisappear and then reappear at other locations. Imagine, for example, a monkey looking for fruits that grow in a set of trees.What happens if one of the trees ceases producing fruits? The monkey, then, has to look for another tree. Time-dependent TSPsof this sort have several real-world applications, such as vehicle routing and machine scheduling problems .Many variations of the TSP have been analyzed in recent decades . Previous research related to the TSP has focusedmainly on producing algorithms to find shortest paths but, to our knowledge, the properties of longer trajectories have notbeen discussed. In the present work we study the statistical properties of all trajectories in two generalizations of the TSP withtime-dependent sites: the TSP with moving sites (bTSP), where sites can be interpreted as ‘boats’ that move gradually in aregion, and the TSP with reallocation of sites (rTSP), where sites move discontinuously across space [Figure 1(a)]. In thepeddlers’ example above, peddlers might not move during one day (so trajectories do not change) but on the next they may havea different position (possibly modifying the shortest path, such that a new optimization process is needed to find it). That wouldbe equivalent to assume that trajectory changes in time-dependent TSPs occur at a much slower time scale than the traversal ofthe paths by salesmen. If we rank trajectories by their length, we can analyze how the properties of this ranking change in timewith measures commonly used in the study of hierarchy dynamics in complex systems, such as the rank distribution f ( k ) andthe rank diversity d ( k ) .Since the rank distribution was popularized in the 1930s by Zipf , it has been used to characterize complex systems ofdifferent nature . We have recently discussed the cases of six indo-european languages and of six sports and games , where wefound that f ( k ) does not follow Zipf’s law and varies slightly across cases. For all of these complex systems, we also studiedhow ranks change in time by means of the rank diversity d ( k ) . Explicitly, d ( k ) is the number of different elements that haverank k within a given period of time T . In the complex systems we have studied up until now, the rank diversity can be wellapproximated by a sigmoid curve. a r X i v : . [ n li n . AO ] A ug a) static TSPTSP withmoving sites TSP withreallocation of sitesa bcd e one step a bcd ea b ce d time (b) k . . . . . . . f ( k ) N = 6 N = 7 N = 8 N = 9 N = 10 Figure 1.
Time-dependent traveling salesman problem with moving sites and with reallocation of sites. Rankdistribution for the TSP. (a) Diagram showing the two generalizations of the TSP considered here. By allowing nodes to shiftin time, the static TSP (center) becomes time dependent. In the TSP with moving sites (left), these are displaced in a smallinterval δ t . In the TSP with reallocation of sites (right), they are displaced to random positions after one time step. In theseTSP variations, nodes do not move during path traversals, that is, the optimization problem of finding paths is solved for a staticsituation. Shortest paths are shown with lines between sites. (b) Rank distribution according to the inverse path length f ( k ) , for N = N =
10 sites in particular but representative instances of the TSP.
Results
Static and time-dependent TSPs
Consider the static TSP with N sites. We shall label each site by a, b, c, and so on. A trajectory is a closed path on these sites,so each trajectory can characterized by a non-unique string of site labels. For example, the shortest trajectory in Figure 1(a)for the static TSP (top) is “aedcb”. Note that one could have started the string at any site, and even changed the direction, sothe same trajectory could have been labeled “cdeab”. Starting from his home city the salesman visits each site only once andreturns home, that is, he follows a trajectory. There are ( N − ) ! / k = k = k max = ( N − ) ! / f ( k ) is taken as the inverse of trajectory length. The rank distribution of TSP trajectories is presented in Figure 1(b) forseveral random configurations corresponding to different values of N . The results differ from Zipf’s law, but for all cases theyhave a similar shape, which resembles a beta distribution.We now study the problem of stability of trajectories in the TSP. Suppose the location of sites varies slowly over time, sothat the salesman can assume a static scenario for each travel, but the shortest trajectory, and in fact the rank of all trajectories,might change if the configuration of sites is modified enough. A toy model that captures this situation is the bTSP. Assume all N sites are allowed to move within a 1 × X and Y components ofthe velocity of each site are random with a uniform distribution between ±
1. The boats move with a constant velocity until theyreach a confining wall, where they bounce elastically [see Figure 1(a) and accompanying video in Supplementary Information].Let us consider a different time-dependent perturbation of the TSP, the rTSP. Here N sites are located in the unit squarewith a uniform random probability, and all trajectories are ranked as explained above. This is the “base” configuration. Thenone site is chosen at random and reallocated to a random position in the unit square, after which trajectory ranks are calculatedagain [Figure 1(a)]. After several iterations, we can explore this time-dependent process with the rank diversity. Even though acontinuous time dynamics does not exist as in the bTSP, we can still ask how many different trajectories occupy a given rank k ,so the rank diversity is well defined. Temporal evolution of trajectory ranks
We first explore how trajectory ranks change with time in both the bTSP and rTSP (Figure 2). For the bTSP, trajectories initiallyat the extremes (i.e. high or low k ) tend to remain in place, until at some point they detach and explore a wider range of ranks[Figure 2(a)]. That is, in the edges of the plot (high and low k ), rank as a function of k tends to be a horizontal straight line, a) t k (b) step number k Figure 2.
Time dependence of trajectory ranks.
Time dependence of trajectory ranks for the models depicted in Figure1(a), where particular instances of N = k axis as time goes by. We show results for the (a) bTSP, where trajectories spread slowly in time, and for the (b) TSP withreallocation of sites [rTSP], where there are drastic changes from one time step to another, as well as a tendency for keepinginitially shorter/longer trajectories short/long.whereas for middle k the corresponding lines vary more. For the rTSP a slightly different behavior is observed [Figure 2(b)].Here we note that time is an auxiliary variable with no relevant meaning, since we are simply selecting random positions torelocate randomly chosen sites. As these positions are taken as a series of random and uncorrelated values, reordering such asequence in time makes no statistical difference. Thus, it is reasonable to expect that trajectory ranks for the rTSP vary morethan for the bTSP. Despite this fact, we also observe that the trajectories near the extremes tend to keep the same values of k ,while intermediate trajectories do not exhibit such regularity, just like in the bTSP.In order to quantify the way in which trajectory ranks change in time, we shall use the rank diversity. Moreover, to developsome intuition on this measure and understand how it depends on the parameters chosen to calculate it, we present severalnumerical experiments for both models in Figure 3. Note that to calculate the rank diversity we have to choose an appropriatetime span T and a number m of time points at which observations are made. Alternatively, the time interval δ t betweenobservations can be chosen instead of m , where T = m δ t . In many real-world systems the time interval δ t is determined bydata availability; to calculate the diversity of words in English throughout the centuries, for example, one may use Google’sn-gram dataset, which implies a time interval δ t of one year (or an integer multiple of that). However, in the bTSP both thetotal time T and the time interval δ t can be chosen at will.Let us analyze the situations depicted in Figure 3. First, we study the bTSP with a fixed total time evolution T and varying m , thus effectively changing the value of δ t = T / m [Figure 3(a)]. As m increases, the time between observations decreases,and the boats move less, so the positions of the boats barely change in a single time interval δ t , and rank diversity diminishes.In fact, the diversity cannot be larger than the maximum number of trajectories divided by the number of time slots, k max / m . If,on the other hand, m is small enough, k max / m is larger than 1 and d ( k ) tends to saturate around 1. This can be summarized inthe formula d max = min (cid:26) , ( N − ) !2 m (cid:27) . (1)In the second scenario, we fix the number of observations and vary the total time T . As a limiting case, when T is much largerthan the time a boat needs to travel from one wall to another, the correlation between site configurations is lost, so for mostranks diversity is close to its maximum; only close to the shortest and longest trajectories we observe smaller diversities [Figure3(b)]. Since both N and m are fixed, in this case d max is the same for all T values. We also analyze the case of a fixed value of δ t . The value of m now changes, as in Figure 3(a), but simultaneously T varies, as in Figure 3(b), due to the relation δ t = T / m .As m increases, the effect on diversity seems similar to the one in which T is constant [Figure 3(c)]. Finally, we perform asimilar analysis for the rTSP [Figure 3(d)]. As m increases, the time between observations is not shortened. Each new positionof the selected site is uncorrelated with the previous one, regardless of the value of m . We thus expect and observe a similareffect as the one of Figure 3(c). In this case, the bound for different values of m does not change, since ( N − ) ! / m >
500 1000 1500 2000 25000.00.20.40.60.81.0 d ( k ) T = 2 • m = 200 • m = 5000 • m = 10 m = 10 • T = 2 • T = 10 • T = 100 δt = 5 × − • m = 200 • m = 5000 • m = 10 • m = 20 • m = 100 • m = 500 (a) bTSP (d)(c) bTSPbTSP rTSPrTSPbTSP(b)(e) (f) • n = 1 • n = 3 • n = 5 • n = 1 • n = 3 • n = 5 k Figure 3.
Parameter dependence of rank diversity.
We analyze how changing parameter values affects the rank diversityfor two time-dependent versions of the TSP. For panels (a)-(d) we consider N =
8. Starting with the bTSP and n = N , we: (a)vary the number of observations m while keeping the total time T =
2; (b) fix the number of time steps m = while varyingthe total time; and (c) fix the value of δ t = × − while changing the number of time steps (and thus the total time). For therTSP, we also (d) show how changing the number of time steps alters the diversity. Curves are averages over 100 realizations.Panels (e) and (f) respectively show the rank diversity for the bTSP when n sites are allowed to move, and for the rTSP when n sites are reallocated at each time step, both with N = m = T =
10, which leads to δ t = .
1. We choose this value of T so that (e) and (f) resemble each other. Results in (f) are averages over 1000 realizations.Horizontal lines correspond to d max , the maximum value of the diversity given in Eq. (1). In (b), (d), (e) and (f) all curves havethe same bound, whereas in (a) and (c) color is used to indicate the bounds of each curve. wo useful generalizations of the bTSP and the rTSP are now considered, since it will allow us to relate the behavior ofboth models. Assume that only some sites move in the bTSP. That is, instead of all site moving in the plane, n move and N − n are static. The cases analyzed before thus correspond to n = N . In a similar way, consider an rTSP where instead of reallocatinga single site, n sites are moved. The case n = N thus corresponds to a total reallocation of the system, while up to now wehave only discussed the value n =
1. Diversity for these generalizations seems to behave in a similar way as for the case n = n , as each of the n sites has a two dimensional space to move. For a giveninitial configuration, in the case n = N = k = , k = k = ρ k of the areascorresponding to random realizations [Figure 4(d-f)]. From these histograms we see that ρ ≈ ρ k max , but ρ is clearly differentfrom ρ k max / .It is also useful to study the expected value of the stability area for different ranks, (cid:104) A (cid:105) ρ k , as a function of the number ofsites (Figure 5). There is a different scaling for the extremal and intermediate rankings, so we expect that the difference inareas seen for the case N = ρ k is relatively large, the average value (cid:104) A (cid:105) ρ k can be determined with high accuracy (the statistical error for points in Figure5 are comparable to the size of the points). We can also see an even-odd effect for the longest trajectory associated with thechange in topology of a star-like configuration for an even or odd number of sites. This hints at a geometric understanding oftime-dependent TSPs which, however, is outside the scope of the present study. Connecting the bTSP and rTSP
Consider the rTSP with n =
1. For a particular rank k , each position in the unit square yields a trajectory. We can thus ‘paint’the unit square with colors corresponding to trajectories. We show examples with N = k = k = k =
6. Thisindicates already that d ( ) < d ( ) if a large ensemble is taken (so that errors due to finite sampling are small enough).For the bTSP with n =
1, the moving boat will cover the whole unit square uniformly for most choices of the velocities.The condition for having a uniform covering is that the vertical and horizontal speeds are incommensurable, which alwaysholds for this model. When this occurs, time averages yield the same value as space averages, i.e., the system is ergodic . Forsuch long times the whole unit square will be visited, i.e. the moving site will visit all colored areas. Therefore, the bTSP hasthe same rank diversity as the rTSP when the sampling m is the same (otherwise they are related by a constant factor).For a larger number of moving sites, n >
1, a similar reasoning holds. However, the configuration space is now 2 n dimensional, and thus visualizing it is more challenging. The condition for ergodicity still holds when all pairs of velocities areincommensurable. However, when the number of boats that move is increased, the diversity changes as in Figure 3(e) (wherethree different values of n are shown). When n = N − N , the diversity is formally maximum and we obtain no relevantinformation. In fact, when the time over which we calculate the diversity in the bTSP increases, d ( k ) tends to 1 [Figure 3(b)].The bTSP and rTSP are comparable since both explore a fraction of the 2 N -dimensional configuration space. The bTSPexplores a line of finite length in such a space, or a 2 N dimensional hypercube embedded in the configuration space if the wholeensemble of velocities is considered. The rTSP, on the other hand, explores a 2 n dimensional hyperplane embedded in the sameconfiguration space. Overall, both models behave in a similar fashion.Each realization of the TSP can be seen as a point in a 2 N -dimensional configuration space, where every pair of axis definesthe coordinates of each particle. The optimization problem is then different for each point in the configuration space. We haveanalyzed the stability of the solutions of the TSP under changes of the location of the point defining the configuration. Wehave further shown that the stability properties are similar for the two time-dependent generalizations of the TSP consideredhere. We have also stated under what conditions the behavior of both models is identical. We thus expect that these results areapplicable to other perturbations of the TSP that involve small variations in configuration space. Discussion
We have studied the statistical properties of rank distributions and rank dynamics of a novel variation of the traveling salesmanproblem where nodes shift their position in time. This allows us to explore the stability of trajectory ranking, which is related to . . . . . . . . . . . . . . . . . . .
00 0 . . . . . . . . . . . . . AA A ρ ( A ) ρ ( A ) ρ ( A ) (a) (b) (c)(d) (e) (f) k = 1 k = 180 k = 360 (cid:104) A (cid:105) ρ = 0 . (cid:104) A (cid:105) ρ = 0 . (cid:104) A (cid:105) ρ = 0 . Figure 4.
Stability of ranks under reallocation of sites. (a) Particular instance of the rTSP with N = n =
1. Blackpoints represent fixed sites, whereas the white one represents the site to be reallocated. The shortest trajectory is indicated as adashed red line. When the white point is reallocated in the green area A , the shortest trajectory remains unchanged, while if anode is reallocated in any part of the blue area, the shortest trajectory changes. (b) Same realization of the rTSP, but for thetrajectory with an intermediate rank ( k = k = A for N = k , as well as its mean value, (cid:104) A (cid:105) ρ k . - - - - - - l og (cid:104) A (cid:105) ρ k N • k = 1 • k = k max / • k = k max Figure 5.
Scaling in rank stability.
Scaling of the average value of the area for which a change in position will not changethe rank of a trajectory, for ranks k = k max / k max , and for different system sizes. a) k = 1 (b) k = 2 (c) k = 6 Figure 6.
Stability areas for different trajectories. (a) rTSP with N = n =
1, and three possible reallocations of thewhite point. The shortest trajectory, k =
1, is indicated as a solid, dashed and dotted line, corresponding to the threereallocations. Regions of varying color correspond to different trajectories for k = k = k =
6, respectively.the predictability of perturbation effects. Figure 5 shows the average probability of a trajectory maintaining its rank of 1000instances of the bTSP (with one site moving) as the number N of sites increases. We see that this probability, which reflects thestability of ranks to perturbations, decreases with N . However, the decrease is much lower for the shortest and longest pathsthan for the middle trajectories. This reflects the fact that the rank diversity is also lowest for the extreme paths. Moreover, thestability decreases much faster for the middle trajectories. Thus, we find that the shortest and longest trajectories are morepredictable and robust. This claim would benefit from analytic solutions for the numerical results presented here, but these arebeyond the scope of the present work.We also note that rank diversity curves are all symmetric, as shown in Figure 3. However, in previous studies we haveshown that the rank diversity d ( k ) has almost the same form for languages, sports and games. This is not symmetric and can beadjusted by a sigmoid in lognormal scale. Still, from other datasets we have also noticed symmetric rank diversity curves. Thedifference between symmetric and non-symmetric rank diversity curves seems to be related with the degree of ‘openness’ ofa system. In time-dependent TSPs the system is ‘closed’, as all trajectories are considered at all times. However, in sports,languages and other real-world phenomena, elements enter and exit the system. More precisely, elements enter and exit thesubset of the system available in datasets. If one plots rank diversity curves of ‘open’ systems in linear scale, they are verysimilar to the symmetric curves of closed systems .Changes in optimization problems pose challenges when change is faster than optimization, as solutions might be obsolete .These non-stationary problems are a common feature of complex systems . Our analysis suggests that the way in which statespaces of non-stationary problems change is not uniform. This implies that once an optimal solution is found, we can expect itto be more stable than non-optimal solutions. A more detailed exploration of the changes of state spaces in time will providefurther insight, and we trust that the methods presented here will contribute to this effort. Data Availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author onreasonable request.
References Lawler, E. L., Lenstra, J. K. & Rinnooy Kan, A. H. The traveling salesman problem (1985). Gutin, G. & Punnen, A.
The Traveling Salesman Problem and Its Variations . Combinatorial Optimization (Springer US,2007). URL https://books.google.com.mx/books?id=pfRSPwAACAAJ . Malandraki, C. & Daskin, M. S. Time dependent vehicle routing problems: Formulations, properties and heuristicalgorithms.
Transportation science , 185–200 (1992). Bigras, L.-P., Gamache, M. & Savard, G. The time-dependent traveling salesman problem and single ma-chine scheduling problems with sequence dependent setup times.
Discrete Optimization , 685 – 699 . DOIhttp://doi.org/10.1016/j.disopt.2008.04.001. Cocho, G., Flores, J., Gershenson, C., Pineda, C. & S´anchez, S. Rank diversity of languages: Generic behavior incomputational linguistics.
PLoS ONE , e0121898 (2015). URL http://dx.doi.org/10.1371%2Fjournal.pone.0121898 . DOI 10.1371/journal.pone.0121898. Morales, J. A. et al.
Generic temporal features of performance rankings in sports and games.
EPJ Data Science , 33 (2016).URL http://dx.doi.org/10.1140/epjds/s13688-016-0096-y . DOI 10.1140/epjds/s13688-016-0096-y. Zipf, G. K.
Selective Studies and the Principle of Relative Frequency in Language (Harvard University Press, Cambridge,MA, USA, 1932). Clauset, A., Shalizi, C. R. & Newman, M. E. Power-law distributions in empirical data.
SIAM Review , 661–703 (2009). Yaglom, A. M.
An introduction to the theory of stationary random functions. (Prentice-Hall, Inc., New Jersey., 1962).
Gershenson, C. Facing complexity: Prediction vs. adaptation. In Massip, A. & Bastardas, A. (eds.)
Complexity Perspectiveson Language, Communication and Society , 3–14 (Springer, Berlin Heidelberg, 2013). URL http://arxiv.org/abs/1112.3843 . Gershenson, C. The implications of interactions for science and philosophy.
Foundations of Science , 781–790 (2013).URL http://arxiv.org/abs/1105.2827 . DOI 10.1007/s10699-012-9305-8. Acknowledgements
We acknowledge support from UNAM-PAPIIT Grant No. IN111015.
Author contributions statement
J.F. and C.P. conceived the ideas behind this work. G.C. obtained the rank distribution for the static TSP. S.S. and C.P. carriedout the numerical simulations and statistical tests. G.I. provided the framework to contextualize the results. J.F., C.P., G.I., andC.G. wrote the manuscript. All authors discussed and reviewed the manuscript.
Additional information
Competing financial interests