#### Abstract

Two generalizations of the traveling salesman problem in which sites change their position in time are presented. The way the rank of different trajectory lengths changes in time is studied using the rank diversity. We analyze the statistical properties of rank distributions and rank dynamics and give evidence that the shortest and longest trajectories are more predictable and robust to change, that is, more stable.

#### 1. Introduction

Imagine that a certain product must be delivered, in the most efficient way, to stores in a region. When stores are fixed at given sites in space, finding the shortest path that links them all is the target of the* traveling salesman problem* (TSP) [1, 2], which here we refer to as static TSP. But what happens if the product must be delivered to moving targets such as peddlers? The TSP becomes time dependent. Another variation of the TSP, related to the peddlers’ scenario, occurs when one or more sites disappear 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 TSPs of this sort have several real-world applications, such as vehicle routing [3] and machine scheduling problems [4].

Many variations of the TSP have been analyzed in recent decades [1, 2]. Previous research related to the TSP has focused mainly on producing algorithms to find shortest paths but, to our knowledge, the properties of longer trajectories have not been discussed. In the present work we study the statistical properties of all trajectories in two generalizations of the TSP with time-dependent sites: the TSP* with moving sites* (bTSP), where sites can be interpreted as “boats” that move gradually in a region, and the TSP* with reallocation of sites* (rTSP), where sites move discontinuously across space (Figure 1(a)). In the peddlers’ example above, peddlers might not move during one day (so trajectories do not change) but on the next they may have a different position (possibly modifying the shortest path, such that a new optimization process is needed to find it). That would be equivalent to assume that trajectory changes in time-dependent TSPs occur at a much slower time scale than the traversal of the paths by salesmen. If we rank trajectories by their length, we can analyze how the properties of this ranking change in time with measures commonly used in the study of hierarchy dynamics in complex systems, such as the rank distribution and the rank diversity [5, 6].

**(a)**

**(b)**

Since the rank distribution was popularized in the 1930s by Zipf [7], it has been used to characterize complex systems of different nature [8]. We have recently discussed the cases of six Indo-European languages [5] and of six sports and games [6], where we found that does not follow Zipf’s law and varies slightly across cases. For all of these complex systems, we also studied how ranks change in time by means of the rank diversity . Explicitly, is the number of different elements that have rank within a given period of time . In the complex systems we have studied that, up until now, the rank diversity can be well approximated by a sigmoid curve.

#### 2. Results

##### 2.1. Static and Time-Dependent TSPs

Consider the static TSP with 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 be characterized by a nonunique 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, so the same trajectory could have been labeled “CDEAB.” Starting from his home city the salesman visits each site only once and returns home; that is, he follows a trajectory. There are different trajectories, and the usual problem is to find the shortest one within this set. In this work we go further and rank each trajectory according to its length, so the shortest one has the highest rank (), and the longest one has the lowest rank (). To study the rank distribution of the system, is taken as the inverse of trajectory length. The rank distribution of TSP trajectories is presented in Figure 1(b) for several random configurations corresponding to different values of . The results differ from Zipf’s law, but for all cases they have a similar shape, which resembles a beta distribution. Notice that any other choice of , like minus the length, or any decreasing function of the length, would change the shape of the curves presented in Figure 1(b). However, the rest of the results presented in this paper, except for Figure 7, would remain unchanged.

We now study the problem of stability of trajectories in the TSP. Suppose the location of sites varies slowly over time, so that the salesman can assume a static scenario for each trajectory, 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 sites are allowed to move within a square as if they were boats. In the TSP with moving sites, the and components of the velocity of each site are random with a uniform distribution among . The boats move with a constant velocity until they reach a confining wall, where they bounce elastically [see Figure 1(a) and accompanying video in https://www.youtube.com/watch?v=murdnwzRHu4].

Let us consider a different time-dependent perturbation of the TSP, the rTSP. Here sites are located in the unit square with a uniform random probability, and all trajectories are ranked as explained above. This is the “base” configuration. Then one site is chosen at random and reallocated to a random position in the unit square, after which trajectory ranks are calculated again (Figure 1(a)). After several iterations, during which the same site is randomly reallocated, we can explore this time-dependent process with the rank diversity. Even though a continuous time dynamics does not exist as in the bTSP, we can still ask how many different trajectories occupy a given rank , so the rank diversity is well defined.

##### 2.2. 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 initially at the extremes (i.e., high or low ) 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 ), rank as a function of tends to be a horizontal straight line, whereas for middle 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 to relocate randomly chosen sites. As these positions are taken as a series of random and uncorrelated values, reordering such a sequence in time makes no statistical difference. Thus, it is reasonable to expect that trajectory ranks for the rTSP vary more than for the bTSP. Despite this fact, we also observe that the trajectories near the extremes tend to keep the same values of , while intermediate trajectories do not exhibit such regularity, just like in the bTSP.

**(a)**

**(b)**

In order to quantify the way in which trajectory ranks change in time, we shall use the rank diversity. Moreover, to develop some intuition on this measure and understand how it depends on the parameters chosen to calculate it, we present several numerical experiments for both models in Figure 3. Note that to calculate the rank diversity we have to choose an appropriate time span and a number of time points at which observations are made. Alternatively, the time interval between observations can be chosen instead of , where . In many real-world systems the time interval is determined by data availability; to calculate the diversity of words in English throughout the centuries, for example, one may use Google’s n-gram dataset, which implies a time interval of one year (or an integer multiple of that). However, in the bTSP both the total time and the time interval can be chosen at will.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Let us analyze the situations depicted in Figure 3. First, we study the bTSP with a fixed total time evolution and varying , thus effectively changing the value of (Figure 3(a)). As increases, the time between observations decreases, and the boats move less, so the positions of the boats barely change in a single time interval , and rank diversity diminishes. In fact, the diversity cannot be larger than the maximum number of trajectories divided by the number of time slots, . If, on the other hand, is small enough, is larger than 1 and tends to saturate around 1. This can be summarized in the formulaIn the second scenario, we fix the number of observations and vary the total time . As a limiting case, when is much larger than the time a boat needs to travel from one wall to another, the correlation between site configurations is lost, so for most ranks diversity is close to its maximum; only close to the shortest and longest trajectories we observe smaller diversities (Figure 3(b)). Since both and are fixed, in this case is the same for all values. We also analyze the case of a fixed value of . The value of now changes, as in Figure 3(a), but simultaneously varies, as in Figure 3(b), due to the relation . As increases, the effect on diversity seems similar to the one in which is constant (Figure 3(c)). Finally, we perform a similar analysis for the rTSP (Figure 3(d)). As increases, the time between observations is not shortened. Each new position of the selected site is uncorrelated with the previous one, regardless of the value of . We thus expect and observe a similar effect as the one of Figure 3(c). In this case, the bound for different values of does not change, since for all cases shown.

Two useful generalizations of the bTSP and the rTSP are now considered, since this will allow us to relate the behavior of both models. Assume that only some sites move in the bTSP. That is, instead of all site moving in the plane, move and are static. The cases analyzed before thus correspond to . In a similar way, consider an rTSP where instead of reallocating a single site, sites are moved (the same sites through each iteration). The case thus corresponds to a total reallocation of the system, while up to now we have only discussed the value . Diversity for these generalizations seems to behave in a similar way as for the case , as seen in Figures 3(e) and 3(f). In fact, for the bTSP one may even consider a single moving site and the previous conclusions still hold. We have obtained similar results for other variations of the bTSP, such as periodic boundary conditions and different ways of choosing the initial conditions and velocities.

Let us analyze geometrically a particular, but random, instance of the rTSP with the goal of gaining some insight. The configuration space of such problem has dimension , as each of the sites has a two-dimensional space to move. For a given initial configuration, in the case we can plot the locus of positions to which a given point may travel, so the trajectory associated with a given rank remains unchanged. In Figures 4(a)–4(c) we take as an example, so the shortest, middle, and longest trajectories correspond to , and , respectively. The example is representative in the sense that, for the extremal trajectories, the aforementioned locus has a big area and the middle one has a very small area, but particular values change broadly from realization to realization. We can also obtain a histogram of the probability density of the areas corresponding to random realizations (Figures 4(d)–4(f)). From these histograms we see that , but is clearly different from .

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

It is also useful to study the expected value of the stability area for different ranks, , as a function of the number of sites (Figure 5). There is a different scaling for the extremal and intermediate rankings, so we expect that the difference in areas seen for the case is exponentially larger for bigger systems. Although the standard deviation of the distribution is relatively large, the average value can be determined with high accuracy (the statistical error for points in Figure 5 are comparable to the size of the points). We can also see an even-odd effect for the longest trajectory associated with the change in topology of a star-like configuration for an even or odd number of sites. This hints at a geometric understanding of time-dependent TSPs which, however, is outside the scope of the present study.

##### 2.3. Connecting the bTSP and rTSP

Consider the rTSP with . For a particular rank , 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 for the shortest trajectory in Figure 6(a) and for and in Figures 6(b) and 6(c), respectively. Diversity is the number of different colors in the picture divided by the number of observations. Notice that there is a qualitative difference between the cases and . This indicates already that if a large ensemble is taken (so that errors due to finite sampling are small enough).

**(a)**

**(b)**

**(c)**

For the bTSP with , 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 always holds for this model. However, notice that the system will not go through all possible configurations of the TSP, as boats are fixed. Still, when the aforementioned condition occurs, time averages yield the same value as space averages; that is, the system is* ergodic* [9] in an appropriately chosen manifold. For such long times the whole unit square will be visited; that is, the moving site will visit all colored areas. Therefore, the bTSP has the same rank diversity as the rTSP when the sampling is the same (otherwise they are related by a constant factor). For example, if an instance of the bTSP is chosen as in Figure 6, with one white point indicating the position of the moving site, for , will saturate to , equal to for the rTSP if the random sampling includes points in all different areas. For , and for the intermediate ranking .

For a larger number of moving sites, , a similar reasoning holds. However, the configuration space is now -dimensional; thus visualization is more challenging. The condition for ergodicity still holds when all pairs of velocities are incommensurable. However, when the number of boats that move is increased, the diversity changes as in Figure 3(e) (where three different values of are shown). When or , the diversity is formally maximum and we obtain no relevant information. In fact, when the time over which we calculate the diversity in the bTSP increases, tends to 1 (Figure 3(b)).

The bTSP and rTSP are comparable since both explore a fraction of the -dimensional configuration space. The bTSP explores a line of finite length in such a space, or a -dimensional hypercube embedded in the configuration space if the whole ensemble of velocities is considered. The rTSP, on the other hand, explores a -dimensional hyperplane embedded in the same configuration space. Overall, both models behave in a similar fashion.

Each realization of the TSP can be seen as a point in a -dimensional configuration space, where every pair of axis defines the coordinates of each particle. The optimization problem is then different for each point in the configuration space. We have analyzed the stability of the solutions of the TSP under changes of the location of the point defining the configuration. We have further shown that the stability properties are similar for the two time-dependent generalizations of the TSP considered here. We have also stated under what conditions the behavior of both models is identical. We thus expect that these results are applicable to other perturbations of the TSP.

#### 3. Discussion

We have studied the statistical properties of rank distributions and rank dynamics of a novel variation of the traveling salesman problem where nodes shift their position in time. This allows us to explore the stability of trajectory ranking, which is related to the predictability of perturbation effects. Figure 5 shows the average probability of a trajectory maintaining its rank of 1000 instances of the bTSP (with one site moving) as the number of sites increases. We see that this probability, which reflects the stability of ranks to perturbations, decreases with . However, the decrease is much lower for the shortest and longest paths than for the middle trajectories. This reflects the fact that the rank diversity is also lowest for the extreme paths. Moreover, the stability decreases much faster for the middle trajectories. Thus, we find that the shortest and longest trajectories are more predictable and robust.

Further light on the problem can be shed by studying the probability density of a random path of a random set of points. This is equivalent to studying the probability density of the trajectory length and is closely related to the derivative of as a function of (see Figure 1(b)). In Figure 7 we see a low density of trajectories for the shortest and longest lengths and an approximately, but clearly not exact, symmetrical shape. Under the assumption that perturbations alter all trajectories in a statistically equivalent way by changing their length by a fixed value, trajectories in the lower and higher ranks would change their rank less than the ones in intermediate ranks. With some effort, an analytical analysis in this direction might shed some light on the problem, but this is beyond 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 have shown that the rank diversity has almost the same form for languages, sports, and games. This is not symmetric and can be adjusted by a sigmoid in lognormal scale. Still, from other datasets we have also noticed symmetric rank diversity curves. The difference between symmetric and nonsymmetric rank diversity curves seems to be related to the degree of “openness” of a 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 the subset of the system available in datasets. If one plots rank diversity curves of “open” systems in linear scale, they are very similar to the symmetric curves of closed systems [5, 6].

Changes in optimization problems pose challenges when change is faster than optimization, as solutions might be obsolete [10]. These nonstationary problems are a common feature of complex systems [11]. Our analysis suggests that the way in which state spaces of nonstationary problems change is not uniform. This implies that once an optimal solution is found we can expect it to be more stable than nonoptimal solutions. A more detailed exploration of the changes of state spaces in time will provide further insight, and we trust that the methods presented here will contribute to this effort.

#### Disclosure

Earlier versions of this work were presented at the Complex Systems Conference in Cancun, Mexico, in September 2017 and at the 6th International Conference on Complex Networks and Their Applications in Lyon, France, in November 2017. An abstract was included in a booklet for both conferences, but a full version of this article has not been published.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The authors acknowledge support from UNAM-PAPIIT Grants nos. IN111015 and IG100518, CONACyT Grant no. 285754, and the Fundación Marcos Moshinsky.