Dynamic Shortest Paths Methods for the Time-Dependent TSP

: The time-dependent traveling salesman problem (TDTSP) asks for a shortest Hamiltonian tour in a directed graph where (asymmetric) arc-costs depend on the time the arc is entered. With trafﬁc data abundantly available, methods to optimize routes with respect to time-dependent travel times are widely desired. This holds in particular for the traveling salesman problem, which is a corner stone of logistic planning. In this paper, we devise column-generation-based IP methods to solve the TDTSP in full generality, both for arc-and path-based formulations. The algorithmic key is a time-dependent shortest path problem, which arises from the pricing problem of the column generation and is of independent interest—namely, to ﬁnd paths in a time-expanded graph that are acyclic in the underlying (non-expanded) graph. As this problem is computationally too costly, we price over the set of paths that contain no cycles of length k . In addition, we devise—tailored for the TDTSP—several families of valid inequalities, primal heuristics, a propagation method, and a branching rule. Combining these with the time-dependent shortest path pricing we provide—to our knowledge—the ﬁrst elaborate method to solve the TDTSP in general and with fully general time-dependence. We also provide for results on complexity and approximability of the TDTSP. In computational experiments on randomly generated instances, we are able to solve the large majority of small instances (20 nodes) to optimality, while closing about two thirds of the remaining gap of the large instances (40 nodes) after one hour of computation.


Introduction
The traveling salesman problem (TSP) is among the best-studied combinatorial optimization problems (see [1,2] for summaries), which is in part due to a wide area of applications in logistics, manufacturing, telecommunications, and more. Considerable effort has been put into theoretical analysis of the problem, and implementations of Branchand-Bound-based codes capable of computing optimal tours for instances consisting of several thousand nodes. Many heuristic solution techniques have been proposed as well, notably the Lin-Kernighan heuristic and its adaptations by Helsgaun [3], and numerous approaches based on genetic algorithms (see [4] for a comparative study). A prominent generalization of the TSP is the asymmetric traveling salesman problem (ATSP), which allows for different costs for the two directions in which a connection may be traversed.
There is, however, a drawback when it comes to real-world traffic problems, namely the fact that the TSP assumes that the cost required for utilizing an arc in the tour is constant, independent of the time at which the arc is traversed. This assumption does not generally hold in urban areas, where congestion, and, as a result, travel times fluctuate considerably. The problem of computing shortest paths in the presence of time-dependent travel times has been studied before, leading to significant algorithmic advances [5] rivaling those made for the shortest path problem itself.
As for the planning of time-dependent tours, several extensions of the ATSP have been examined in the literature (see below). The resulting formulations are, however, generally focused on the structure of the time-dependent cost functions rather than that of the underlying network itself. Conversely, we study the time-dependent TSP (TDTSP for short) focusing on paths through the network, leading to different Branch-Cut-and-Price formulations with time-dependent shortest path problems being solved during the column generation.
The TSP is very versatile in applications. For solving TSP, specialized IP methods have proven particularly fruitful. The importance of time-dependent cost functions is obvious at least for logistic applications. Thus, development of specialized IP methods for the TDTSP in full generality, with fully general time-dependence of the arc lengths, is of fundamental importance. To this end, we develop both arc-and path-based formulations for the TDTSP. Both formulations naturally lead to a high number of variables, naturally lending themselves to be solved by a column generation approach. It turns out that a pivotal step in this approach is the pricing problem, which yields a time-dependent shortest path problem of independent interest. This problem is to find and price paths in the time expansion of a graph G, which are acyclic when projected in the underlying graph G.

Preliminaries
The TDTSP is a generalization of the ATSP to the case of time-dependent cost functions. Specifically, we let D = (V, A) be a directed graph with n: = |V| ≥ 3 vertices, and θ max ∈ N a fixed time horizon for Θ: = {0, . . . , θ max }. The travel time for an arc a ∈ A is given by a function c a : Θ → N. For each sequence of arcs (a 1 , . . . , a k ) with a k = (u k , v k ) and v k = u k+1 , we can recursively define an arrival time θ arr (a 1 , . . . , a k ) := c u 1 ,v 1 (0), if k = 1, and θ arr (a 1 , . . . , a k−1 ) + c u k ,v k (θ arr (a 1 , . . . , a k−1 )) otherwise.
A feasible solution of the TDTSP is a tour beginning at a source vertex s ∈ V, visiting every other vertex exactly once, which then returns to s. The TDTSP asks for a feasible solution T = (a 1 , . . . , a n ) minimizing the arrival time θ arr (a 1 , . . . , a n ) back at s, which we will also denote by c(T). The specification of a source vertex is necessary for the TDTSP, as opposed to the ATSP, due to the time-dependency of the travel times c. Several special properties of travel time functions play an important role in time-dependent versions of combinatorial problems: The first-in-first-out (FIFO) property stipulates that the traveler who first enters an arc is also the first to leave it again. Formally, it must hold for each a ∈ A that θ + c a (θ) ≤ θ + c a (θ ) ∀ θ, θ ∈ Θ, θ ≤ θ .
Shortest paths with respect to time-dependent costs can be computed efficiently using a variant of Dijkstra's algorithm if the cost functions satisfy the FIFO property [6,7].
Secondly, several well-known results (e.g., [8,9]) state that the symmetric version of the TSP can be approximated in the case of metric cost coefficients, that is, cost coefficients satisfying the triangle inequality. The definition of the triangle inequality can be easily generalized to the time-dependent case-a set of travel time functions satisfies the timedependent triangle inequality if it holds for each u, v, w ∈ V with θ, θ + c uv (θ) ∈ Θ that θ + c uw (θ) ≤ θ + c uv (θ) + c vw (θ + c uv (θ)). (1)

Related Work
Several generalizations of the ATSP have been considered in the literature, such as the TSP with time windows [10,11], or the class of vehicle routing problems (VRPs) [12]. An interesting novel approach to the TSP and its variant is related to path-based formulations-the flow conservation constraints of the TSP ensure that every solution corresponds to a set of cycles in the underlying graph, making it possible to reformulate the problem in terms of path variables and solving it using a Branch-Cut-and-Price framework [13]. The same holds for the VRP [14] and other TSP-related problems [15]. By itself, this reformulation is not stronger than the original formulation due to Dantzig et al. [16]. It is, however, possible to restrict the set of path variables in order to exclude paths which are not tours. Of course, any such modification alters the pricing problem, generally having an adverse effect on its computational tractability, requiring a balance between the quality of the relaxation and the computational effort required to solve the pricing problem. Promising approaches include the generation of k-cycle-free [17], as well as so-called ng paths [18].
An early extension of the TSP to incorporate time-dependency, due to Picard and Queyranne [19], generalizes the travel time of an arc (i, j) ∈ A to be sequence-dependent, i.e., a function c ij (k) (k = 1, . . . , n). This sequence-dependent variant of the TSP, which we refer to as the SDTSP, has since been studied both theoretically [15] and experimentally [19,20]. Notably, there is a correspondence between this sequence-dependent generalization of the TSP and the Minimum Latency Problem (MLP), which asks for a tour minimizing the sum of waiting times of customers. The correspondence has inspired some mixed-integer programming (MIP) formulations [21] for the MLP. The SDTSP is also closely related to identical machine scheduling, in particular P|| ∑ w j T j (see [22]).
Some attempts have been made [23,24] to solve the TDTSP itself using Branch-and-Cut algorithms, with the focus of exploiting a special structure of the time-dependent cost functions. Specifically, the travel-time model introduced in [25] assumes that travel times are determined by a piecewise constant travel speed function, leading to a three-index formulation where the total number of variables depends on the complexity of the travel speed functions in terms of the number of their break points. While the approach seems successful for highly structured travel speed functions, the computational tractability apparently degrades for more irregular instances (see, for instance, the computational results in [23]). Another approach, going by the name of Dynamic Discretization Discovery, has been proposed [26,27] specifically for the TDTSP with time windows, which relies on iteratively refining a discretization of the time horizon according to optimal solutions of coarser discretizations determined during previous stages. The approach essentially substitutes one large integer program with several well-chosen smaller ones, exploiting the combination of time-dependent travel times and time windows. Conversely, we pursue an approach that does not rely on any special condition for the travel times.
A preliminary version of this paper [28] explored the (ultimately unsuccessful) usage of machine-learning techniques in order to solve the TDTSP.

Structure of This Paper
We begin by examining the complexity of the TDTSP in Section 2, establishing that the problem is hard to approximate even if an instance satisfies both the FIFO property and the time-dependent triangle inequality. On the other hand, we establish that an ATSP approximation algorithm can be used to approximate the TDTSP under some more restrictive circumstances.
We proceed in Section 3 by introducing several formulations for the TDTSP which do not need any specific structure of the underlying travel time functions. The formulation is based on a time expansion of the original graph, resulting in a potentially large number of variables and constraints, necessitating some form of column generation. We describe the specific pricing problem, which corresponds to computing shortest time-dependent paths, that is, shortest paths in a time-expanded network. We augment the approach by computing time-dependent k-cycle-free paths, using a dual stabilization technique in order to decrease the number of required pricing iterations. We add several valid inequalities, a propagation method, a custom branching rule, and primal heuristics in order to improve the solution process.
We demonstrate the effectiveness of our approach in Section 4 based on a computational experiment on several instances of differing sizes, providing a conclusion in Section 5.

Complexity
As a generalization of the TSP, the TDTSP is N P-hard itself. Recall that while there exists no α-approximation for any α > 1 for the general TSP, the metric TSP can be approximated [29] (p. 557). This does not hold for the metric TDTSP: Theorem 1. There is no α-approximation algorithm for any α > 1 for the TDTSP unless P = N P, even if the FIFO property and the time-dependent triangle inequality are satisfied.
Proof. Suppose there exists an α-approximation algorithm ALG for some α ≥ 1. We show that algorithm ALG could be used to solve the HAMILTONIAN CYCLE problem on an undirected graph G = (V, E). To this end, let D = (V, A) be the bidirected complete graph with costs Note that G is Hamiltonian if, and only if D contains a tour with costs of exactly n with respect to the cost function d. Let M := αn + 1, s ∈ V be an arbitrary source vertex, and θ max := n · M be a time horizon. For θ ∈ Θ and a ∈ A, let be a set of time-dependent cost functions. These cost functions satisfy both the FIFO property and the time-dependent triangle inequality (1). We distinguish two cases with respect to the TDTSP instance (D, c, s, θ max ): -G has a Hamilton cycle, and therefore, a tour T = (a 1 , . . . , a n ) with costs of n exists in D: In this case, the values d(a i ) for i = 1, . . . , n are all equal to one, implying that the optimal solution to the TDTSP instance has an arrival time of n. ALG yields a tour T apx with an arrival time of, at most, αn < M. Consequently, the first case of (2) is never attained, implying that each arc has a travel time of one, and ALG correctly identifies a Hamiltonian cycle in G. -G does not possess a Hamiltonian cycle, and every tour T = (a 1 , . . . , a n = (u, s)) has a cost of at least n + 1 with respect to d: As a result, there must be one arc a k ∈ T having a travel time of at least two. If k < n, then T arrives at u at time θ ≥ n, and the travel time of a n , and consequently the arrival time of T is at least M. If, on the other hand, k = n, then T arrives at u at time n − 1 and arc a n has a travel time of more than one, which is only possible if d uv = 2, yielding a travel time of M for a n . In any case, the arrival time of T is at least M > αn. Based on this distinction, G has a Hamilton cycle if, and only if ALG produces a tour with an arrival time of n.

Remark 1 (Dynamic Programming).
It is well-known [30] that the (asymmetric) TSP can be solved by using a dynamic programming approach. Let C(S, v) be the smallest cost of an (s, v)-path consisting of the vertices S ⊆ V with s, v ∈ S, and s = v. Then, C(S, v) satisfies the following recursive relationship: The cost of an optimal tour is then given by min v =s C(V, v) + c(v, s), and can be determined by computing the values C(S, v) in increasing order of |S|, yielding an O(2 n · n 2 ) algorithm for the ATSP.
In order to adapt the approach to the TDTSP, we can define C(S, v) as the minimum arrival time of any (s, v)-path departing from s at θ = 0, consisting of the vertices S ⊆ V with s, v ∈ S and s = v. If the TDTSP instance in question satisfies the FIFO property, C(S, v) satisfies a similar relationship: These relations yield an algorithm for the TDTSP with the same running time as its counterpart for the ATSP. Note that the FIFO property is key here, since it ensures that only the shortest path for fixed (S, v) needs to be considered for subsequent computations.

Approximation for Special Cases
Letč : A → N be the static underestimator of the time-dependent cost function c, defined asč If the time-dependent cost functions are of low variance, the underestimator yields TDTSP-approximations based on an underlying ATSP: Then, any α-approximation of the ATSP yields an (αλ)-approximation of the TDTSP.
Proof. Let T opt , T apx be the optimal and α-approximate tour with respect to the costsč and T T opt be the optimal TDTSP tour. We have that Since the ATSP is inapproximable in general, further assumptions, such as a metric lower boundč uv , are necessary to obtain any approximation results. What is more, we show in Appendix A that while it is possible to generalize the well-known one-tree relaxation [31] from the ATSP to the TDTSP, even the computation of the corresponding lower bound is an N P-hard problem.

A Branch-and-Price Algorithm
Our aim in the next section is to provide a state-of-the-art mixed integer programming (MIP)-based algorithm to solve the TDTSP in full generality. We begin by formally defining time-expanded graphs and establishing both an arc-based and a path-based formulation, since it is not immediately clear which approach will work better. Both formulations consist of large numbers of variables, necessitating column generation approaches. Since feasible tours correspond to acyclic paths through the time-expanded graph, the key to solving the TDTSP efficiently is: How can we find a shortest path through the time-expanded graph which is acyclic on the original graph?
Optimizing over the set of acyclic paths, which we denote by P * , corresponds to solving the TDTSP itself. The set P of all paths through the time-expanded graph is much easier to handle, leading to a very fast pricing algorithm at the cost of a substantially weaker relaxation. To balance the computational effort of the pricing step and the strength of the relaxation, we generate k-cycle-free paths, that is, paths containing no cycles of size ≤ k, where larger values of k increase the computational effort while improving the relaxation.
After discussing the subject of column generation, we strengthen both formulations based on valid inequalities. To this end, we review several classes of valid inequalities for the SDTSP and other related problems, adapting the class and their respective separation algorithms to the TDTSP. Lastly, we make several improvements relating to well-known MIP techniques, such as adding a branching rule and several primal heuristics.
We denote an arc (u θ , v θ ) by (u, v, θ) and assume from now on that c uv (θ) > 0 for all (u, v) ∈ A, θ ∈ T (v). It follows that D T is acyclic.
Example 1. Figure 1 shows a directed graph with travel times for each arc and its time expansion. Any tour on D can be embedded into D T as an (s 0 , s θ )-path.

An Arc-Based Formulation
We formulate the TDTSP based on binary variables for each arc of D T similar to the SDTSP formulation in [19]. The formulation has two kinds of constraints: A number of covering constraints ensures that each vertex in V has exactly one outgoing arc in D T , while flow conservation constraints guarantee that any feasible solution consists of a single (s 0 , s θ )-path: (TDTSP-A) Remark 2. By virtue of the flow conservation constraints in (TDTSP-A), any solution of the program or its LP-relaxation can be decomposed into a set of paths leading from vertex s 0 to vertices s θ for some values θ ⊆ T (s). As a result, there exist several equivalent linear objectives, for example the arrival time objective, given by

Relation to the Static ATSP
There is a strong relationship between the TDTSP and the ATSP. Let y : A → R ≥0 be the compound flow traversing an arc (u, v) ∈ A with respect to a feasible solution (x a ) a∈A T of (TDTSP-A), given as The covering constraints and the flow conservation yield the well-known two matching equations, The covering constraints and the integrality of the variables x imply that the values y are binary for any feasible solution x. However, a correct static ATSP formulation still requires subtour elimination constraints (SECs) of the form Since D T is acyclic, any solution of (TDTSP-A) is guaranteed to satisfy the additional SECs. Equivalently, flow augmentation techniques [32] for strengthening ATSP formulations are redundant for (TDTSP-A). Conversely, SECs are not necessarily satisfied by fractional solutions. Thus, (TDTSP-A) can be strengthened by separating SECs with respect to the underlying static ATSP (see Section 3.4).
Valid inequalities for the ATSP can be included in the TDTSP by first computing the compound flow y * of a solution x * of the LP-relaxation of (TDTSP-A), executing some separation algorithm yielding a valid inequality in the compound variables y, and finally expressing this inequality in terms of the original problem variables x.
Any feasible TDTSP solution is feasible for the underlying ATSP, and generic ATSP solutions do not necessarily produce feasible solutions of the TDTSP. Specifically, no tour T = (a 1 , . . . , a n ) with θ arr (T) > θ max can be embedded into D T . The complete description of the TDTSP in terms of compound variables y can be obtained by adding forbidden path constraints of the form ∑ a∈P y a ≤ k − 1 ∀ P = (a 1 , . . . , a k ) : θ arr (P) > θ max .
As a result, facet-defining ATSP inequalities, while valid, are not necessarily facetdefining for the TDTSP.

A Path-Based Formulation
Any feasible solution of the TDTSP problem (TDTSP-A) can be decomposed into (s 0 , s θ )-paths for some values θ ⊆ T (s). These paths correspond to cycles containing vertex s in D. We denote by P the set of all (s 0 , s θ )-paths not containing any vertex s θ with θ = θ . For each path P ∈ P and each vertex v ∈ V, we let χ v,P be the set of vertices in V T that have an outgoing arc contained in P, and α v,P its cardinality, that is, The problem can be reformulated in terms of path variables: (TDTSP-P) Any solution of this formulation consists of a single variable x P set to 1 and all others set to 0, in which case P must be contained in P * , that is, corresponding to a tour. Any fractional solution consists of, at most, n different paths in P which need not be tours in D.
The resulting system is small in terms of the number of constraints at the expense of the number of variables, necessitating a column generation [33] approach.

Column Generation
Let (λ v ) v∈V be the dual variables corresponding to the covering constraints in (TDTSP-P). The reduced cost of a variable x P in the corresponding LP-relaxation is given as which can be rewritten as a function A T → R via The pricing problem therefore corresponds to a shortest path problem in D T , which can be solved in linear time since D T is acyclic.

Lagrangean Pricing & Dual Stabilization
It is fairly straightforward to derive a Lagrangean relaxation of the LP-relaxation of (TDTSP-P), which is solved during the pricing problem. The constraint ∑ P∈P α s,P · x P = 1 is equivalent to ∑ P∈P x P = 1, since every path in P leaves s 0 exactly once. By penalizing all covering constraints in the objective and removing all but this constraint, we obtain the Lagrangean relaxation If there exists a path with negative reduced costs with respect to λ, an optimal solution of L P (λ) is obtained by setting x * P = 1, where P ∈ P has the minimum reduced cost. If there is no such path, we set x * ≡ 0. In any case, we obtain a lower bound on the LPrelaxation of (TDTSP-P) during the pricing. It is easy to adapt this approach to (TDTSP-A), enabling us to use path-based pricing in this case as well.
We use the Lagrangean relaxation in order to perform a dual stabilization based on the weighted Dantzig-Wolfe decomposition method introduced in [22]. The weighted Dantzig-Wolfe decomposition method judges the quality of dual values according to the value of their Lagrangean relaxation. By maintaining a center of stability and only tentatively moving the center towards the current dual values, the technique can significantly decrease the time required to solve the LP-relaxations of the TDTSP formulations.

Pricing Acyclic Paths
Many paths which are generated in the pricing do not share much resemblance with tours in the underlying graph D: On the one hand, certain paths only contain few vertices and lead almost immediately back to s θ . We will address this problem later using the propagation of lower bounds. On the other hand, paths frequently contain cycles with respect to D, that is, they contain two different copies A path of at least k arcs contains a k-cycle if any of its subpaths is a k-cycle. A path not containing a j-cycle for j ≤ k is called k-cycle-free. Naturally, a k-cycle-free path is also j-cycle-free for all j < k, and it is a tour if, and only if it is (n − 1)-cycle-free.
In order to strengthen the formulations, we restrict the set of variables in the pathbased formulation (TDTSP-P) to the subset P k ⊆ P of k-cycle-free paths, noting that P =: P 1 ⊆ P 2 . . . ⊆ P n−1 = P * .
Similarly, the sets P k yield increasingly tighter Lagrangean relaxations, i.e., satisfying The restriction to the set P k requires us to adapt the pricing procedure to k-cycle-free paths. The k-cycle-free shortest problem has previously been examined [17], leading to a two-cycle-free labeling scheme with linear running time, as well as an algorithm computing k-cycle-free paths with a running time of O((k!) 2 ) in the parameter k. In order to improve the overall performance of the column generation, we use the weighted Dantzig-Wolfe decomposition described above based on the values L P k (λ).
Adapting the original formulation (TDTSP-A) is not as straightforward, since we have no control regarding the path decomposition producing the values x uv,θ . As a compromise, we propose to generate new columns according to the arcs of k-cycle-free paths. Since we cannot guarantee that the resulting flow can be decomposed into k-cycle-free paths, we cannot ensure that L P k (λ) ≤ LP (TDTSP-A). Nonetheless, it still holds that L P k (λ) ≤ (TDTSP-A), which is sufficient to ensure that, if no more k-cycle-free paths can be found, the restricted LP is a relaxation of (TDTSP-A). We are, however, unable to apply the weighted Dantzig-Wolfe decomposition to the arc-based formulation.

Valid Inequalities
We consider several classes of valid inequalities, some of which are well-known ATSP inequalities, whereas others are either adaptations of SDTSP inequalities or newly derived ones. We express the inequalities in terms of the arc variables x from (TDTSP-A), understanding that converting them to (TDTSP-P) is straightforward. The separation is performed according to a fractional solution x * with compound values y * given by

ATSP Inequalities
Apart from the subtour elimination constraints, probably the best-known family of facet-defining inequalities for the ATSP goes by the name of D + k -inequalities [2,34]. These inequalities are derived based on the compound variables y on a complete directed graph D = (V, A). To simplify notation, for sets S, T ⊆ V we let The D + k -inequality for a sequence (v 1 , . . . , v k ) of 2 ≤ k < n distinct vertices is given by The separation of D + k -inequalities involves the enumeration of possible sequences in a Branch-and-Bound-like fashion. Nonetheless, the separation works well in practice, since many of the possible sequences can be pruned.

Incompatibilites
Incompatibilites between binary variables, collected in the so-called incompatibility graph, have proven to be highly useful in order to generate strong inequalities for mixed-integer programs in general (see [35]). For the ATSP, two arcs (u, v) = (u , v ) are incompatible if they share one or both of their endpoints, that is, The two common classes of inequalities derived from incompatibility graphs are clique and odd-cycle inequalities. Cliques in the incompatibility graph of the ATSP correspond to the sets δ + (v) and δ − (v). For the (TDTSP-A), clique inequalities are already implied by the constraints.
Odd-cycle inequalities for the ATSP are known as odd closed alternating trail inequalities [36], odd CATS for short. An odd CAT corresponds to a sequence a 1 , . . . , a 2k+1 of distinct arcs in A such that arcs a i and a i+1 share either head or tail. The odd CAT inequality corresponding to these arcs is given by Odd CAT inequalities for the ATSP can be separated heuristically by computing shortest paths in an auxiliary bipartite graph. For the TDTSP, these cuts correspond to odd cycles of cliques rather than odd cycles of simple incompatibilities.

Odd Path-Free Inequalities
Let S ⊆ V \ {s} be a set of vertices of the underlying graph, V T (S) := {u θ ∈ V T | u ∈ S} the corresponding vertices in V T , and A T (S) the arcs of the subgraph of D T induced by V T (S). The SEC corresponding to S implies that the set A T (S) can contain, at most, |S| − 1 arcs.
We obtain another class of inequalities by restricting ourselves to subsets with odd cardinality, i.e., |S| = 2k + 1. A subset A ⊆ A T is path-free with respect to S if it contains no nontrivial paths in D T , a nontrivial path consisting of more than one arc. If a set is path-free, than any arc in the intersection of A with a tour through D T connects two distinct vertices. As a result, the intersection can contain, at most, k arcs, yielding the following odd path-free inequality (see Figure 2): We separate these inequalities heuristically by first determining promising subsets S, and then computing the optimal set A for each promising subset S. In order to determine a promising subset S, we note that the odd path-free inequality is strongest for small values of k. We therefore restrict ourselves to the case of k = 1, determining several promising sets based on the amount of induced flow y * (S). To avoid very similar inequalities, we only consider the best promising 3-set containing each vertex v ∈ D \ {s}. For each candidate set S, we compute an optimal set A using an integer program.

Lifted Subtour Elimination Inequalities
Subtour elimination constraints can be used to strengthen formulation (TDTSP-A). They can be separated in polynomial time by solving a series of flow problems on the underlying graph D. To further strengthen SECs for the SDTSP, it was observed in [15] that any tour has to leave S sufficiently early in order to be able to reach the vertices in V \ S and return to s within the time horizon. We proceed to adapt the approach to the TDTSP. Letθ be such that θ ≥ max{θ | there exists a tour T leaving S for the first time at θ}.
Then, the following lifted subtour elimination constraint (LSEC) is valid for (TDTSP-A): The value ofθ is maximized for a tour which first serves S, then V \ S and returns to s immediately afterwards. Computing the optimum value ofθ is intractable in practice, as it would involve the solution of a TDTSP on S itself.
Instead of maximizing the time spent in S, we therefore minimize the time spent in V \ S. Since we do not know the optimum value ofθ, we do not know the values of the cost function c. Hence, we relax the problem by replacing the time-dependent cost functions c uv (θ) by its static underestimatorč uv . Minimizing the length of a static tour leaving s, traversing V \ S, and finally returning to s is easily formulated as an ATSP, yielding an optimal value ofθ. The latest departure timeθ can then be determined as θ max −θ.
This estimation of the optimalθ can be rather rough if the cost functions c differ significantly from their underestimatorsč. As a result, the computational effort required to solve the ATSP does not seem merited. We choose to instead bound the static tour cost from below by computing an arborescence of minimum weight.

Cycle Inequalities
While it is possible to generate variables for the (TDTSP-A) according to k-cyclefree paths, solutions of the LP relaxation generally do not possess a k-cycle-free path decomposition. We adapt the cycle inequalities introduced in [15] for the SDTSP to the TDTSP, cutting off solutions based on the elimination of k-cycles. Consider a path P = (a = (u, u 1 , θ), a 1 = (u 1 , v 1 , θ 1 ), . . . , a k = (u k , v k , θ k )) in D T satisfying the following properties: 1.
The arcs a 1 , . . . , a k form a k-cycle C not containing s.
Let T be a tour entering the cycle C via arc a. Since T is k-cycle-free, it can contain, at most, j ≤ k − 1 of the arcs in C, leaving the cycle via an arc (u j , w, θ j ) with w = u j+1 . It also holds that w cannot coincide with u for < j or u (otherwise T would contain a cycle of length ≤ j). Based on these observations, the cycle inequality associated with P, defined as is valid for (TDTSP-A). We separate these inequalities using a simple heuristic: Based on x * we compute a decomposition into paths P 1 , . . . , P k ⊆ P. For each path P i , we determine whether it contains a k-cycle, adding the corresponding inequality if necessary.
Note that cycle inequalities were reexamined in [37], leading to the stronger class of Time-Dependent Cycle Inequalities (TDCIs) for the SDTSP. The techniques required to derive these inequalities are, however, specific to the SDTSP and do not generalize to the TDTSP.

Unitary AFCs
Unitary admissible flow constraints (AFCs) were introduced in [15] for the SDTSP. In terms of the TDTSP, they can be derived as follows: Summing up the flow conservation constraints of (TDTSP-A) over a set S ⊆ V T not containing any vertices s θ yields the equation x(δ − (S)) = x(δ + (S)). Based on an incoming arc (u, v, θ) ∈ δ − (S), it can be relaxed to x uv,θ ≤ x(δ + (S)). This constraint by itself is obviously redundant. However, any tour T that enters S via (u, v, θ) has to leave S using some arc (u , v , θ ) ∈ δ + (S) such that v = u, v, yielding the unitary AFC inequality To separate these inequalities, we first observe that stronger unitary AFC inequalities correspond to smaller cuts, that is, a set S with (u, v, θ) ∈ δ − (S ) and δ + (S ) ⊆ δ + (S) produces a stronger inequality than S. As a result, any vertex in S except for v θ+c uv (θ) can be removed, provided that it does not contain any arc entering from another vertex in S. We can therefore assume that S only contains vertices reachable from v θ+c uv (θ) . Similarly, we can also assume that S contains no copies of u and v, other than v θ+c uv (θ) itself. Hence, the separation of unitary AFC inequalities can be conducted by solving a min-cut problem for each arc (u, v, θ) with x * uv,θ > 0, with capacities based on the values of x * .

Additional Techniques
The addition of cutting planes already significantly strengthens the TDTSP formulations. There are, however, several other techniques which can be used to speed up the computation of the optimal tour in a Branch-and-Bound framework:

Propagation
During the solution process, we have a dual bound θ based on the value of the LP-relaxation of the current Branch-and-Bound node and a primal bound θ based on the best-known integral solution (we let θ := −∞ and θ := +∞ if the bounds are not available). The LP-relaxation frequently contains (s 0 , s θ )-paths containing only a few arcs each, leading back to s at θ < θ. Similarly, there are long paths which cannot be part of an optimal solution. Formally, variable x vw,θ can be fixed to zero if θ > θ or (θ + c vw (θ) < θ and w = s).
We propagate any improvement in the bounds θ and θ by fixing the appropriate variables. We also include the propagation into the pricing problem.

Compound Branching
Branching on variables x uv,θ likely yields a highly unbalanced Branch-and-Bound tree, since fixing a variable x uv,θ to one is a strong restriction, while fixing it to zero is a very weak one. We instead add the binary compound flow variables y and coupling constraints (3) to the formulations and assign branching priorities in order to force the solver to branch on compound flow variables whenever possible, leading to a more balanced tree. Since any binary solution y * corresponds to a tour, it is never necessary to branch on the variables x. Since the number of compound variables is generally much lower than the number of arcs in A T , the addition of the compound variables is unlikely to significantly affect computational performance.

Primal Heuristics
We use an incremental construction heuristic to obtain a path P starting at s 0 . At each turn, we append an arc leading to vertices whose counterparts in D are still unexplored by P, finishing when the path forms a tour. The selection of the arcs is based on the current fractional solution (x * , y * ). Arcs that are fixed to zero are always disregarded. If there are multiple possible arcs to be added, we score each arc (u, v, θ), and then choose one with probability proportional to the score. We use three different scoring functions: - The inverse of the travel time c uv (θ) (breaking ties arbitrarily); - The variable value of x * uv,θ using travel times to break ties or - The compound value y * uv using the same tie-breaking rule. Since this construction is computationally inexpensive, we can increase the chance of finding an improved tour by using all scoring functions and performing several iterations with different random seeds.

Computational Experiments
We proceed to report on the empirical findings on a set of artificially generated TDTSP instances. All experiments were carried out based on an implementation in the C++ programming language that compiled with full optimization. We used version 7.0.1 of the SCIP [38] optimization suite and version 9.0.0 of GUROBI [39] as an underlying LP solver. We ran the experiments on an AMD Epyc 7742 processor clocked at up to 3.40 GHz, and imposed a time limit of 3600 s for all computations. Each individual computation was carried out by a single-threaded process.

Instances
We generate several problem instances, each given by a complete directed graph and cost functions associated with its arcs. We embed the vertices of the graph into {0, . . . , 100} 2 and introduce (symmetric) costs c a using rounded-down Euclidean distances between the points of the embedding. Based on the costs c a , we generate a piecewise linear function f a (θ): N → Z with M ∈ N breaking points θ 1 < θ 2 < . . . < θ M within Θ:

1.
We let f a (0): = 0 and fix the slope at zero to +1.

2.
The slope alternates between +1 and −1 with break points at θ i for i = 1, . . . , M. Based on f a , the time-dependent cost function c a : {0, . . . , θ max } → N is given as where the parameter λ > 1 controls which multiple of c a (θ)/c a can be attained (see Figure 3 for By construction, the function f satisfies the FIFO property independently of the choice of the time steps θ. We then use (time-dependent) shortest-path distances with respect to c(·) instead of c directly, in order to ensure that the (time-dependent) triangle inequality is satisfied as well.
We generate 50 small instances consisting of 20 vertices each, as well as 20 large instances consisting of 40 vertices based on different random seeds. For each instance, we compute an optimal tour T 0 with respect to the time-independent costs. The tour T 0 serves both as an initial feasible solution and as a means to determine a suitable time horizon θ max to construct a time-expanded graph D T .
There is a significant increase in the number of arcs during the time-expansion-while the original number of 20 vertices of the small instances leads to 380 arcs in the underlying graph, the number of arcs in the time-expanded graphs D T ranges from about 80,000 to over 130,000. For the large instances, the number of arcs in D T lies between 800,000 and 1,000,000.

Remark 3.
By means of Theorem 2, the choice of objective enabled us to compute a lower bound on the optimal objective values of our TDTSP instances. These bounds are oftentimes better than the objective values of the LP-relaxations of the root nodes. However, in order to evaluate our formulations, we do not use them during our computations.

Comparisons between Formulations
We begin by comparing the performance of the different formulations, namely, the arc-based formulation (TDTSP-A), both with and without path-based pricing, and its path-based counterpart (TDTSP-P). The increased size of the instances (see above) has immediate consequences regarding the computational performance (see Table 1 for the solution statistics for the large instances): Formulation (TDTSP-A) spends an average of about 13 min solving the root relaxations of the problem instances. Enabling column generation decreases this time to about 55 s, whereas the path-based formulation (TDTSP-P) averages at 17 s. While the difference in running times is less pronounced for the small instances, it is clear that column generation is necessary in order to solve larger TDTSP instances.
Despite the fact that the path-based formulation solves the initial relaxation the fastest, it explores only seven nodes on average during its execution, compared to 29 nodes explored by the arc-based formulation. Consequently, the average remaining gaps are 0.57 and 0.55, respectively. The gap here is defined as (p − d)/d, with p being the objective function value of the best-known feasible solution, and d the best-known lower bound.

Performance of Different Pricers
We proceed to evaluate the pricing methods introduced in Section 3.3. As a first step, we study the difference between pricing individual arcs and entire paths with respect to the arc-based formulation (TDTSP-A). While both approaches yield the same gaps at the root node, path-based pricing is substantially faster than pricing arcs, while yielding relaxations with substantially fewer variables.
We then examine the effect of pricing k-cycle-free paths, based on the algorithms introduced in [17], the experimental results being depicted in Table 2. Clearly, increasing the value of k increases the computational effort. On the other hand, larger values of k yield tighter relaxations. To obtain a realistic picture of the running times, we solved the LPrelaxation of (TDTSP-A) and measured the required running time, thereby accounting for deviations in the size of the programs arising from the use of different pricing techniques. After solving the relaxations, we compute the relative gaps between the different dual bounds and the primal bound obtained from the original path-based pricing method, thereby ignoring any distortions due to improved primal bounds.
We find that, independently of the formulation, pricing two-cycle-free paths incurs a negligible computational overhead of less than a minute, while closing almost 30% of the relative gap. Increasing k beyond two immediately results in a significant increase to the computational effort, requiring almost 40 min for large instances for the arc-based formulation. What is more, the effect of pricing k-cycle-free paths on the relative gap becomes less pronounced for larger values of k.
Regarding the formulations, the pricing procedure becomes significantly more difficult for the path-based formulation (TDTSP-P). While the effect is mitigated when we add dual stabilization, the difference remains significant. A likely explanation is that, when solving the arc-based formulation, the LP solver can compose the arcs of previously added paths in order to obtain new paths, which may not have been explicitly added before. As a result, the solver is quickly able to compute near-optimal dual values, thereby guiding the column generation. What is more, the difference between the formulations with respect to the relative gaps is very slight and does not merit the computational effort required.
Based on these observations we proposed to use the arc-based formulation while pricing two-cycle-free paths in order to achieve a reasonable trade-off between computational overhead and decrease in relative gap. We evaluate the separators introduced in Section 3.4 using the same approach as above, comparing the running times and relative gaps for all separators (see Table 3). In terms of the relative gap, the subtour inequalities, as well as their lifted counterparts (see Section 3.4.4) perform best, closing about 25% of the gap of the original formulation. The computational overhead is relatively lower compared to pricing three-or four-cyclefree paths. We evaluated miscellaneous combinations of different separators, finding the combination of lifted subtour and D + k inequalities to be most efficient. In order to judge the effectiveness of the LP-based construction heuristics introduced in Section 3.5.3, we solve the initial relaxations of instances, then apply the heuristics in order to improve the primal bound. Once more, we record both the execution time and the gap relative to the original dual bound. The results, shown in Table 4, demonstrate the effectiveness of these fairly simple heuristics-while the computational effort is minor, with increases in running time well below 10 s, up to 30% of the gap is closed solely based on improved primal bounds. Based on the previous experiments, we augment the arc-based formulation (TDTSP-A) with two-cycle-free path-based pricing, the combination of lifted subtour and D + k inequalities, and all of the primal heuristics. In addition, we use the propagation method introduced in Section 3.5.1 to fix binary variables whenever possible.
Based on the combined improvements, we are able to solve 46 of the 50 small instances to optimality. The remaining gap of the four unsolved instances averages 4.5%.
While we are unable to solve any of the large instances to optimality, we reduce the original remaining gap from 55% to 19% (see Tables 1 and 5). In order to compare our algorithm with a state-of-the-art solver, we used GUROBI as a black-box solver on the same instances. As a black-box solver, GUROBI uses neither a column-generation approach, nor the specialized heuristics and valid inequalities we introduced above, consequently averaging at a significantly larger gap of 44% on the large instances.

Conclusions and Future Work
In this paper, we have discussed several theoretical and empirical properties of the TDTSP. Since the TDTSP is a generalization of the ATSP, many of the complexity-specific theoretical results, such as N P-hardness and inapproximability, carry over to the TDTSP.
Unfortunately, several positive results regarding the ATSP are not retained in the TDTSP. Specifically, the TDTSP remains inapproximable even if a generalized triangle inequality is satisfied. Furthermore, even simple relaxations, such as time-dependent trees, cannot be used to determine combinatorial lower bounds on the TDTSP.
From a practitioner's point of view, the increase in problem size poses significant problems when trying to solve even moderate-sized instances of the TDTSP. The authors of [15] conclude that there are challenging instances of the SDTSP with less than 100 vertices. While the results date back some years, the increase in computational complexity is apparent even in the case of the SDTSP.
To be able to tackle the TDTSP, a sophisticated pricing routine is necessary. It is apparent that generating columns according to time-dependent paths is superior to generating columns for individual arcs. More recent advances regarding specialized shortest path variants facilitate the solution process by achieving improved dual bounds. Further improvements in this area will likely benefit the TDTSP as well.
The connection between TDTSP and ATSP yields a variety of feasible classes of inequalities which help to significantly improve dual bounds. Unfortunately, the generalizations of SDTSP-type inequalities do not perform equally well in comparison. Objective value propagation and primal heuristics decrease the gap even further. The primal heuristics profit from the connection to the ATSP, significantly outperforming the heuristics built into the solver itself.
We have focused on developing an algorithm capable of solving instances with fairly arbitrary travel times. It may be worth investigating the performance on instances with more regular travel time functions by conducting a computational study. Similarly, the existing formulations can be extended to incorporate time windows, suggesting further comparative experiments.

Data Availability Statement:
The implementation of the algorithm as well as the generator for all random instances is available at https://github.com/chrhansk/time-dependent-tsp.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Time-Dependent Spanning Trees
Relaxations play an important role in integer programming in general and the TSP in particular. They provide lower bounds which can be used to obtain quality guarantees for solutions. The prevalent relaxation of combinatorial problems formulated as integer programs is given by their fractional relaxations. In several cases however, it is possible to derive purely combinatorial relaxations. In case of the symmetric traveling salesman problem, a popular combinatorial relaxation is given by one-trees [31]. A one-tree with respect to a graph G = (V, E) with V = {1, . . . , n} consists of a spanning tree together with an additional edge joining vertex s: = 1 to another vertex v ∈ V. Since every tour is a one-tree, the one-tree of minimum cost provides a lower bound on the cost of a tour. The computation of such a one-tree in the static case involves the computation of a minimum spanning tree (MST), which can be carried out efficiently. The approach generalizes to the ATSP, in which case the one-tree corresponds to an arborescence with root s together with a single arc entering s from another vertex v.
An important question to ask is whether these approaches carry over to the timedependent case. We begin by generalizing the definition of spanning trees to the timedependent case. To this end, we assume throughout this subsection that D is bidirected, corresponding to an undirected graph G D , and that the time-dependent cost function c is symmetric, i.e., c uv (θ) = c vu (θ) holds for all (u, v) ∈ A, θ ∈ Θ, making the cost function c correspond to the undirected edges of G D . We consider a spanning tree T = (V, F) of G D . For each vertex v ∈ V, there exists a unique (s, v)-path P v in T, enabling us to define an arrival time θ arr T (u) induced by T via P u for each u ∈ V. Based on these arrival times we define the following: Definition A1. A time-dependent minimum spanning tree (TDMST) is a spanning tree T of G D minimizing c(T), which is defined as c(T) := ∑ {u,v}∈F c uv (θ arr T (u)).
The corresponding optimization problem asks for the value of a TDMST for an instance (G D , c, s). As was the case before, any tour becomes a spanning tree once we remove its last arc, the arrival time of the resulting path at its target coinciding with the cost of the tree. Consequently, the optimal value of the corresponding TDTSP instance is bounded from below by the cost of any minimum spanning tree. Unfortunately, approximating the TDMST is an N P-hard problem itself: Theorem A1. There is no α-approximation algorithm for any α > 1 for the TDMST problem unless P = N P.
Proof. Consider an instance of the 3SAT problem (see [29] (p. 391)), given in terms of a set X: = {x 1 , . . . , x n } of Boolean variables and a set Z: = {Z 1 , . . . , Z m } of clauses, each clause consisting of exactly three literals in X∪ X. We construct a suitable instance of the TDMST problem using a number of components. First we define a component A i for each variable x i ∈ X. The component is shown in Figure A1a, the edges being annotated with their (constant) travel times. We define a component B i for each clause Z i ∈ Z. To this end, let x i1 , x i2 , x i3 denote the variables whose literals appear in Z i and M ≥ 2n + 6m + 1. The edges in the component have the following travel times: 1.
The edges {w i1 , w i2 } and {w i2 , w i3 } have a constant travel time of 1.
The same holds true for the edges {v i2 , w i2 } and {v i3 , w i3 }.

3.
The travel time of the edge connecting x ik and v ik depends on whether x k or x k appears in the clause Z j . In the former case the travel time is constantly 1, whereas in the latter it is given by c x ik v ik (θ) := 0 if θ ≥ 2, and 2 otherwise.
The instance including the components, depicted in Figure A2, has a TDMST of cost less than M if and only if the 3SAT instance is satisfiable: Figure A2. The TDMST construction used to prove Theorem A1.
Consider a satisfying truth assignment of the 3SAT instance. For each variable x j set to true we choose the path (s, s j , x j , t j ) in component A j . If the variable is set to false we choose the path (s, s j , t j , x j ) to be part of the spanning tree. Thus, the arrival time of the tree at x j is 1 if x j is set to true and 2 otherwise. For each clause Z i we add the edges between the vertices corresponding to its variables and their respective v i counterparts. Since the clause is satisfied, the arrival time at one of v i1 , v i2 , or v i3 is at most 2, enabling us to add the edges {v ik , w ik }, {w i1 , w i2 }, and {w i2 , w i3 } at a cost of 1 each. The total cost of the resulting tree is at most 2n + 6m < M.
Conversely, consider a tree T with costs less than M. We first consider the (s, s j )-paths in T for all j = 1, . . . , n. If this path does not contain {s, s j } itself, then it must start by t traversing a variable gadget A j for j = j continued by a clause gadget B i , entering via v ik and leaving via v il . As was established above, the earliest arrival time at v ik is 2, implying that the time of arrival at w il is at least 4, resulting in a travel time of at least M for the traversal of edge {w il , v il }, thereby contradicting the assumption of c(T) < M. For the same reason, T must contain either the path (s, s j , x j ) or (s, s j , t j , x j ). Based on the optimality of T, the tree must contain (s, s j , x j , t j ) in the former case. We can therefore identify the choice of paths with a variable assignment as we did above.
To see that the assignment is feasible, we note that for each clause Z i one of the edges {v ik , w ik }, say {v i1 , w i1 }, is in T and directed away from v i1 . Since c(T) < M, the travel time of this edge must be less than M, implying that the arrival time of T at v i1 is at most 2. If x i1 appears in Z j , then the arrival time at x i1 is 1, and the variable is set to true satisfying Z i . If x i1 appears in Z j , then it follows that the arrival time at x i1 is 2, the variable is set to false and clause Z i is satisfied as well. Now assume the existence of an α-approximation for the TDMST problem. We construct the instance introduced above for M := α · (2n + 6m) and apply the approximation. If the resulting tree has costs less than M, the 3SAT instance is satisfiable. Otherwise, the optimal TDMST has costs at least M/α ≥ 2n + 6m, and the instance is not satisfiable.