We develop a restricted dynamical programming heuristic for a complicated traveling salesman problem: a)cities are grouped into clusters, resp. Generalized TSP; b)precedence constraints are imposed on the order of visiting the clusters, resp. Precedence Constrained TSP; c)the costs of moving to the next cluster and doing the required job inside one are aggregated in a minimax manner, resp. Bottleneck TSP; d)all the costs may depend on the sequence of previously visited clusters, resp. Sequence-Dependent TSP or Time Dependent TSP. Such multiplicity of constraints complicates the use of mixed integer-linear programming, while dynamic programming (DP) bene ts from them; the latter may be supplemented with a branch-and-bound strategy, which necessitates a \DP-compliant" heuristic. The proposed heuristic always yields a feasible solution, which is not always the case with heuristics, and its precision may be tuned until it becomes the exact DP.
Our object is a particularly complicated version of the well-known Traveling Salesman Problem (TSP), which combines several its generalizations that are usually treated separately: Bottleneck TSP [23, Ch. 15], Generalized TSP [23, Ch. 12], Precedence Constrained TSP [40, 18, 29, 8], Sequence-Dependent3 TSP [3] (as a generalization of the more well-known Time-Dependent TSP [20]). Nevertheless, their combination is neither a purely scholastic e ort nor art for art's 3 Note that the Sequence Dependent designation is mostly applied to scheduling problems, see [4], and has a di erent meaning: the cost of present action only depends on the cost of the previous one, not on the entire previous sequence. To the best of the author's knowledge, the only term for dependence on the whole (previous) sequence was State Dependent, introduced in [35] for a scheduling problem concerning printed circuit board assembly; nevertheless, we prefer to follow to the `Sequence Dependent' designation of [3]. sake. Its possible applications include printed circuit boards design (a less general version, without precedence constraints, was considered in [25]) and optimization of the cycle time for industrial robots (for a survey of robotic task sequencing, refer to [2]). See the relevant reviews of TSP and its variations in [34, 37, 31, 23, 30]. For a most recent treatment of Precedence Constrained TSP (TSP-PC), also known as Sequential Ordering Problem (SOP), see [19].
A Terminological note We call the variation of TSP we consider in this paper a Sequence-Dependent Precedence Constrained Bottleneck Generalized TSP (SDBGTSP-PC); however, since the sequence dependence does not currently play an important part in our model, we will mostly omit that designation and refer to our problem as BGTSP-PC, bearing the possible sequence dependence in mind. The designation TSP-PC for Precedence Constrained TSP is borrowed from [8]; we nd it appealing, since it poses no risk of confusion with the very di erent Prize Collecting TSP, which is also abbreviated PCTSP [23, Ch. 14], makes an explicit reference to the \ordinary" TSP (in contrast with SOP), and does not speci cally mention asymmetricity (our method is symmetricity-agnostic), in contrast with PCATS [6]. Another important issue is that we consider an open problem of the TSP family, i.e., return to origin is not mandated; to the best of our knowledge, the open TSP was rst posed in [15], with a reference to a 1965 report by N.Deo and S.L.Hakimi, as \Shortest Hamiltonian Chain Problem".
The rest of the paper is as follows: in Sect. 1 we describe general notation and de nitions; Sect. 2 describes the problem statement, and Sect. 3 completes the description with de nitions of dynamic programming subproblems and the Bellman Equation. Section 4 describes the exact dynamic programming for our problem and Sect. 5 describes our experience of shared memory parallelization of this algorithm. Section 6 discusses the proposed heuristic, reports and compares the results of experiments with the parallel implementation of the exact algorithm and the proposed heuristic. Sections 1{4 follow the most recent paper on exact solution of the problem [13], 1
We employ the standard set-theoretic notation (quanti ers, propositional connectives, etc.); , denotes equality by de nition. Each set, all elements of which are sets themselves, is called a family. For every two objects a and b, denote by fa; bg the (unique) set that contains a, b, and nothing else. In the case a = b, this yields a singleton fag = fbg. We employ the standard Kuratovskii ordered pair de nition: for two arbitrary objects u and v, their ordered pair (OP) is de ned by (u; v) , ffug; fu; vgg; its rst element is u and the second one is v. For an OP z, pr1(z) denotes its rst element and pr2(z) denotes its second element; these are uniquely de ned by the condition z = (pr1(z); pr2(z)); in case z 2 A B, where A and B are sets, we have pr1(z) 2 A and pr2(z) 2 B. We employ the usual canonical representation of ordered triplet [16, x 1.3]: for three objects a, b, and c, we assume (a; b; c) , ((a; b); c). A similar convention is used for Cartesian product of three sets: for arbitrary sets A, B, and C, we have A B C , (A B) C [16, x 1.3]; it obviously means that (x; y) 2 A B C 8x 2 A B 8y 2 C. In connection with this, let us also recall the convention concerning the notation for values of function of three variables: for sets A, B, C, and D, function h : A B C ! D and elements 2 A B and 2 C, in accordance with the above-mentioned representation of A B C, it is valid to consider the element h( ; ) 2 D to be de ned.
As usual, [0; 1[ , f 2 Rj0 g (R is the real line). For each nonempty set S, denote by R+[S] the set of all (nonnegative) functions from S to [0; 1[. As usual, N , f1; 2; : : :g. Assume N0 , f0g [ N and p; q , fi 2 N0j(p i) ^ (i q)g 8p 2 N0 8q 2 N0. Note that the latter de nition yields the empty set if p > q.
For a nonempty nite set K, let jKj 2 N be the power of the set K; then, let (bi)[K] denote the set of all bijections of the \interval" 1; jKj onto K. In particular, for a xed N 2 N, let P , (bi)[1; N ] be the set of all permutations of the \interval" 1; N ; for each 2 P, there exists a permutation 1 2 P such that ( 1(k)) = 1( (k)) = k 8k 2 1; N . Denote by P(H) (P0(H)) the family of all (all nonempty) subsets of set H; let Fin(H) be the family of all nite sets from P0(H). 2
Here and below, x a nonempty set X, where everything happens, a point x0 2 X, which is called the base, a natural number N , N 2, which is the main dimension parameter, sets
M1 2 Fin X; : : : ; MN 2 Fin X; referred to as megalopolises, and relations
M1 2 P0(M1
M1); : : : ; MN 2 P0(MN
MN ):
For j 2 1; N , OPs z 2 Mj describe the possible ways of conducting interior
jobs inside the megalopolis Mj : pr1(z) determines the entry point and pr2(z)
determines the exit point. The scheme of movements is as follows:
x0
!
!
!
pr1 z(
M(jin) , fpr1(z) : z 2 Mj g 8j 2 1; N ; M(jout) , fpr2(z) : z 2 Mj g 8j 2 1; N :
X , Xin ,
x0
x0 [ x0 [ [ N [ Mi(in) i=1
N
[ Mi ;
i=1
N
[ Mi(out) ;
i=1
[ f?g ;
clearly, Xin X Xout X. Like (
The author is aware of the three means to de ne precedence constraints: { Partial order [46]; { Directed acyclic graph [18]; { Nondescript binary relation (a set of OPs) [11].
Clearly, these approaches produce the same results (all may be reduced to a
kind of binary relation, note also the result of [1] on path information) and their
4 the terms \city" and \point" are used interchangeably
(
Recall that P = (bi)[1; N ] is the set of all (complete) routes; it is a nonempty set of cardinality jP j = N !. In terms of address pairs, the set of feasible routes is expressed as follows (see [11, Pt. 2]):
A ,
n
for an address pair z 2 K and \times" t1 2 1; N and t2 2 1; N . In addition to
the route, we also choose the track, or trajectory, which is determined (
At last, we can proceed to the de nition of the quality criterion for our problem. To account for the in uence of the set of pending tasks on the cost function (recall that our problem is sequence dependent ), we will need the symbol N , P0(1; N ) (we call an element of N a task set ). Now, let us introduce the following N + 1 cost functions c 2 R+(Xout
Xin
N); c1 2 R+(Xin
N); : : : ; cN 2 R+(Xin
(
N
max
t20;N 1
c(pr2(zt); pr1(zt+1); f (s) : s 2 t + 1; N g)+
+ c (t+1)(zt+1; f (s) : s 2 t + 1; N g) :
(
The interior jobs setting introduced in [10] provides a means for universal expression of what we are to do inside a megalopolis, thereby unifying the Generalized TSP [23, Ch. 13] and Clustered TSP [14] approaches: Generalized or International TSP. Each cluster is to be visited exactly once, i.e., only one city is to be visited per cluster. To adapt our statement to these requirements, we set 8i 2 1; N Mi , (b; b) : b 2 Mi , i.e., we mandate exit at the point of entry for every megalopolis, and set the rudimentary zero interior job costs: 8i 2 1; N 8b 2 Mi8K 2 N ci(b; b; K) = 0.
Clustered TSP. For each cluster, all of its cities must be visited contiguously before proceeding to the next cluster, i.e., we have an open TSP (Shortest Hamiltonian Chain Problem [15]) inside each cluster thus we can never exit a cluster at the city we used to enter it, hence 8i 2 1; N Mi , (a; b) : (a; b 2 Mi) ^ (b 6= a) , and the costs of interior jobs are set to the costs of the respective open TSP tours starting at the respective entry points.
The problem statement for BGTSP-PC is as follows: C( ) (zi)i=0 =
N
max t20;N 1 c(pr2(zt); pr1(zt+1); f (s) : s 2 t + 1; N g)+
+ c (t+1)(zt+1; f (s) : s 2 t + 1; N g) ; C( ) (zi)i=0 ! min;
N 2 A; (zi)iN=0 2 Z( ): (BGTSP-PC) Denote the value (extremum) of the problem by V ,
V , min min 2A (zi)iN=02Z( )
C( ) (zi)i=0 2 [0; 1[;
N
(
A dynamic programming (DP henceforth) solution consists of embedding a problem into a family of similar problems and obtaining a relation between the extrema and solutions of less di cult problems of the family and the full problem; this relation is expressed through a Bellman (recurrence) equation (see the general description of the method in [38, Ch. 9]). Note that the method below is an example of backwards DP akin to [7], whereas the forward DP [24, x1.2] seems to be more common as applied to the problems of the TSP family.
Precedence constraints pose a challenge when we attempt to reduce the full problem to a series of subproblems since they are formulated with respect to a complete set of megalopolises (in non-generalized TSP, cities) 1; N , and it is not immediately clear how to `extend' precedence constraints to its subsets. The idea is to regard the set of megalopolises, which is to be visited in some subproblem, as a `pre x' (the forward approach of [24, 33, 46, 8]) of some feasible route, or a `su x' (the backwards approach of [7, 11]) thereof. Naturally, such pre x or su x has to be ordered consistently with the precedence constraints on a full route. Since we chose the backwards approach, in the following de nitions we default to su xes, with a possible passing reference to pre xes.
To properly describe the `feasible su xes', we need several new de nitions. First of all, a subset K of the complete task set 1; N is considered feasible if 8z 2 K (pr1(z) 2 K) ) (pr2(z) 2 K). For pre xes, the implication is reversed, see [24, x1.3]. Any task set encountered below is to be assumed feasible unless explicitly stated otherwise, i.e., let us rede ne N to be not the set of all task sets but the set of all feasible task sets.
Infeasible task sets do not in uence the solution of the problem, and it is this lack of in uence that makes it possible to apply DP to precedence constrained problems of the TSP family thanks to a reduction of state space. For a quanti cation of this lack of in uence and the respective reduction of state space, refer to [46, 47, 43].
For an arbitrary task set K 2 N, consider the set [K] of the address pairs that are \saturated" in K: We call the mapping I a crossing-out rule. For a task set Ke , it speci es the possible \entry points", from which it is possible to start a walk through Ke . From the partially ordered set perspective, the mapping I produces a speci c maximum antichain, retaining, for each chain present in Ke , only its minimum element. In a forward procedure, the respective mapping would retain the maximum element, see [46].
We can now de ne the set of partial routes through a task set K 2 N, which
we described above as `feasible su xes':
(I
bi)[K] ,
This set is non-empty for feasible K 2 N (see the proof in [11, Pr. 2.2.2,2.2.3]).
Note that feasible complete routes (
The complete problem is to visit the set of megalopolises 1; N , starting from a separate point, the base x0. Subproblems resemble the full problem inasmuch as they involve visiting a feasible subset K 2 N, starting from some x 2 Xout. However, there has to be an additional restriction on the base point for subproblems: the movement from the base to an element of K we decide to visit rst must be feasible with respect to precedence constraints. Once more, the mapping I is used to formally express this feasibility: let x 2 Xout belong to Mi(out) for some i 2 1; N n K. Clearly, the movement from x to the element of K we will visit rst is to be considered feasible if and only if the megalopolis Mi may occupy the rst place in a partial route over fig [ K; thus, x 2 Mi(out) is said to be a feasible base for the task set K if and only if K [ fig is a feasible set and i 2 I K [ fig .
We also need to de ne the set of partial tracks that agree with a given partial
route. For K 2 N, a feasible base x 2 Xout, and 2 (I bi)[K], let
Z(x; K; ) ,
(zi)jiK=0j 2 ZK
z0 = (f?g; x) & zt 2 M (t) 8t 2 1; jKj
: (
( ) (zi)jiK=0j , CK
max
t20;jKj 1
c pr2(zt); pr1(zt+1);
(s) : s 2 t + 1; jKj
+
+ c (t+1) zt+1;
(s) : s 2 t + 1; jKj
:
(
min min 2(I bi)[K] (zi)jiK=0j2Z(x;K; )
( ) (zi)jiK=0j 2 [0; 1[:
CK
(
Let us supplement the de nition of v(x; K) with the trivial case K = ?; to nd out which x 2 Xout we can use, we may actually apply the de nition of feasible base to the empty set. This yields the points x belonging to all megalopolises that are not senders, i 2 1; N j 8z 2 K i 6= pr1(z) , and thus quali ed to terminate a feasible (full) route. We need this trivial case to prime the recurrence procedure for of V ; set
v(x; ?) , 0 8x 2 Xout: This initial condition serves the needs of the open TSP-like problems (see [15]). To accommodate for the more oft-cited closed TSP setting, where it is mandated to return to the starting point after walking through 1; N , it is only necessary to replace zero costs with the costs of going from x to the starting point. It is equally easy to introduce a more general terminal cost function reminiscent of that in the optimal control theory; it is explored in [11].
At long last, we are prepared to state the Bellman equation: for K 2 N and a feasible base x 2 Xout, v(x; K) = + cj z; K ; v pr2(z); K n fjg o: (BF) For the proof, refer to [13]; it is not too complex yet rather laborious, in no small part due to the generalized clustered character of the problem. 4
In this section we describe the exact DP solution of (BGTSP-PC), on which the
heuristic is based. Recall that earlier we have rede ned the set N to be the set
of all feasible task sets. Yet, we nd it more appealing to use a di erent symbol
in this section: let G , K 2 P0 1; N 8z 2 K (pr1(z) 2 K) ) (pr2(z) 2 K)
be the set of all feasible task sets; note the exclusion of the empty task set. It is
regarded as \feasible", but it is more straightforward to treat it separately. We
now proceed to describe the procedure of generating and traversing the set of all
feasible states, where to traverse a state means to obtain its value (
First of all, partition the feasible task sets according to their cardinality, n Gs ,
K 2 G
o s = jKj 8s 2 1; N : Note the boundary elements of this partition: the last one, GN = 1; N is the singleton re ecting the complete task set. The rst element G1 is the set of all \nonsenders", which are eligible to terminate a feasible route; clearly, no \sender" megalopolis is eligible to do so. The remaining elements of the partition are de ned in a recurrent way,
Gs 1 = nK n ftg : K 2 Gs; t 2 I(K)
o
8s 2 2; N :
(
Let us now describe the procedure that constructs the states based on feasible task sets. Consider the feasible expansion of a feasible set,
J (K) , j 2 1; N n K
fjg [ K 2 Gs+1 2 P0(1; N n K);
(
M[K] , [ M(jout);
j2J (K) [
K2Gs and construct the appropriate states,
D[K] , (x; K) : x 2 M[K] : Collecting D[K] for all feasible K of a certain cardinality, we obtain the layers of state space,
Ds ,
Ds[K] 2 P0(Xout
Gs) 8s 2 1; N 1; let us supplement this de nition with the two boundary cases D0 , (x; f?g) : x 2 M(jout); j 2 G1 and DN , (x0; 1; N ) .
From the feasible task sets, the state space layers inherit a connection between
the elements of neighbouring cardinality:
y; K n fkg 2 Ds 1 8s 2 1; N 8K 2 Gs 8k 2 I(K) 8y 2 M(kout):
(
Technically, we consider the restrictions of v( ; ) (
vs(x; K) , v(x; K) 8(x; K) 2 Ds: Then, we say that for all s 2 1; N , the function vs 2 R+[Ds] is obtained from the function vs 1 2 R+[Ds 1] through \strati ed" (BF), vs(x; K) = 5 In an ordinary, non-generalized TSP, there is no need for a separate set M[K], or rather, this set would match the feasible expansion of K.
Note that feasible sets and, therefore, feasible states, are generated
top-tobottom (
After the value V = vN x0; 1; N of the complete problem is found through backwards traversal of (BFs), we need to obtain the actual solution of the problem, the optimal route and track. Starting with vN , at each step from vs to vs 1, append js 2 I(K) and zs 2 Mjs that yield the corresponding extremum in (BFs) to the end of optimal route and track. There may be multiple optimal routes and tracks, more so for a bottleneck problem, and it is possible to obtain all of them through this procedure; however, we were only concerned with nding some optimal solution. i : 5
Parallel implementation of exact DP The algorithm speci ed in the previous section was implemented in C++ with the work sharing done through the OpenMP shared memory multiprocessing API. A similar parallel implementation was reported in [21] for the additive GTSP-PC without sequence dependence; a similar non-OpenMP based on C# threads was reported in [27] for the same problem. A di erent parallelization strategy for GTSP-PC was reported in [12]. The divide-and-conquer approach of [33] may be applied to GTSP-PC in a way that will not invalidate precedence constraints [44]; the author is not aware of applications of the divide-and-conquer strategy to DP for precedence-constrained problems.
Let us rst recall the general structure of the algorithm:
1. Prime the feasible state generation with the layer D0 of states with empty
task sets. Generate the feasible states D1; : : : ; DN through procedure (
Recall that, for each k 2 1; N , to generate the next state space layer Dk
with the aid of bottom-up procedure (
The parallel work sharing through OpenMP happens in the treatment of the current set of feasible task sets Gi. Each of the feasible task sets can be processed independently, and since the states corresponding to di erent task sets never coincide, no race conditions occur in generation and pricing of the feasible states set Di. On the contrary, it is possible that some K0 2 Gi+1 could be generated from both K1 2 Gi and K2 2 Gi, K1 6= K2, a trivial example of which is the nal GN ; see the paragraph below for an exception to this rule. To avoid race conditions, an omp critical directive had to be placed around the operation of writing to Gi+1 to prevent another thread from writing there at the same time. Better algorithms for generating the feasible task sets would eliminate this race condition altogether since they never generate a new feasible task set more than once, for more details on those algorithms, see the list [9, Apx. 2.2] and the corresponding papers. The one we used was related to the algorithm from [33] and hence at most quadratic in the total number of feasible task sets.
Speaking pedantically, the precedence constraints K may in fact mandate a total order on 1; N . In that case the route is actually xed from the start and only an optimal track is to be determined; a Generalized TSP under such conditions is known as Serdyukov TSP, and there is also Ordered Cluster TSP [23, p. 26], with the usual di erence between Generalized and Clustered TSP as regards the cities to be visited in each megalopolis. In this case, it is futile to do a parallel processing of elements of Gi since each of them is a singleton. Still, the corresponding state space layers Di do not degenerate into singletons and it is possible to share the work by pricing the newly generated elements of each Di in parallel.
Data structures. Our principal data structure was nested hash table implemented with the aid of std::unordered map (hence C++11). Its type can be speci ed as std::vector <std::unordered map <uint32 t, std::unordered map <uint16 t, float> > >, where uintXX t stands for XX-bit unsigned integer type. Semantically, to access the cost of walking through K 2 G [ f?g starting from x 2 M[K], i.e., v(x; K), we had to call the subscript operator three times: [|K|][K][x]. The task sets were coded in the standard bit-masking way; our implementation provided for up to 31 megalopolises; zero was reserved for the base. The cities were numbered such that the numbers for all cities in each megalopolis were consecutive. A nested structure let us avoid the unnecessary repetition of 32-bit task set labels at every state that corresponds to the same task set, of which there is a non-negligible number in generalized problems; for non-generalized problems, the at structure as speci ed in, for example, [36], may be justi ed. An additional bene t of nested structure was the ability to generate the next (in the sense of task set cardinality) feasible states layer by only going through the corresponding set of feasible task sets, without resorting to examining the whole list of present feasible states, which is considerably larger, more so in a generalized problem such as ours.
We adopted hash table as a base data structure because of its amortized constant-time search and the lack of a need to remove elements once installed; since both keys we used were unsigned integer numbers, we did not have to invent our own hash function. Our rst experiments were with std::map, which o ered logarithmic search time, but the switch to std::unordered map made the computation time decrease by a factor of 4 (we admit to only comparing their performance on a single test problem), at which point we settled with the hash table. On the ip side, the use of a hash table did little to help decrease the memory footprint, this perennial problem of dynamical programming for TSPs; it also became rather di cult to predict that footprint. We did not do rigorous memory pro ling, but we can still say that 4GB of RAM were enough for the solution of the 27{25{25 problem, and 30{25{25 (see the problem descriptions in Sect. 6) took more than 15GB and less than 40GB. The use of hash table also mandated the implementation of OpenMP tasks work-sharing construct, which only became available in the third version of the shared memory API, with a single task being the processing of a single feasible task set, i.e, the generation and pricing of its corresponding states and the generation of \successor" feasible task sets. The pseudocode of OpenMP tasks implementation is listed below, where the blocks are de ned by indentation instead of braces as would be t C++.
F#oprraegamcah osm2p 1p;aNrallel default (shared)
1 #pragma omp single nowait
#pragma omp task untied firstprivate(K) Generate J (K)
#pragma omp critical (expand) foo[|K|+1].add(K [ j )
Compute vs(x; K) #pragma omp critical (costwrite)
foo[jKj][K][x]:=vs(x; K) 5.1
Experiment For problem descriptions, refer to Subsect. 6.1. Here and below, the absolute computation time and the scaling factor (serial run time divided by the run time considered) are expressed as \MM:SS || Ratio". The times speci ed include input-output operations. The following results were reported in [42]: srl par-4 par-8 27{10{25{NO 04:58 ||1 01:25 ||3.506 00:44 ||6.773 27{20{25{NO 38:45 ||1 10:34 ||3.667 05:18 ||7.311 27{25{25{NO 73:32 ||1 20:27 ||3.596 10:13 ||7.197
They were all obtained on the Uran supercomputer (for details, see http: //parallel.uran.ru/node/6 [in Russian]) at IMM UrB RAS. The Uran supercomputer ran 64-bit Scientific Linux 6.4; the compiler used was GCC 4.4.7, optimization level -O2. Evidently, the speed-up was near-linear, as much as it could scale for a shared-memory implementation. Eight was the maximum number of cores per node that was available at the time.
Computation times from a more recent implementation of the algorithm, which featured a streamlined code, a probable source of improvement through aiding optimizing compiler rather than trying to actually optimize the code, are listed below: srl par-2 par-3 par-4 27{10{25{NO 02:57 ||1 01:27 ||2.034 01:00 ||2.95 00:46 || 3.848 27{20{25{NO 22:34 ||1 11:28 ||1.968 08:08 ||2.775 05:53 || 3.836 27{25{25{NO 42:56 ||1 22:23 ||1.918 14:48 ||2.901 11:23 || 3.772 These computations were carried out on the author's PC (Intel Core i3-3450) in 64-bit Windows 7 environment. The compiler used was Intel C++ 15 Update 1 since Microsoft Visual C++ would not support the required OpenMP 3.0 API; the compilation parameters were copied over from the default \Release" con guration, with the obvious addition of a ag that allows OpenMP code generation. The superlinear speedup as seen on the two-core calculation for 27{10{25{NO is most probably caused by OpenMP tasks overhead.
As evident from the two tables above, the proposed shared-memory parallel implementation scheme that works through OpenMP tasks provides a reasonable speedup for the problems considered, with the main barrier being the number of cores on a single node that could work through OpenMP. It seems to be possible to make a similar message-passing implementation, however, the latter seems to be sensible only for the problems of greater magnitude than those speci ed above, either through a greater number of megalopolises or cities therein, or less strict precedence constraints, as it would certainly impose a greater overhead. However, a rudimentary message-passing implementation might be used to divide the work between two nodes with the aim of using only half as much memory on each node as required for a single-node (shared memory) implementation, see [33, 44]. 6
Restricted DP heuristic. Experiment The heuristic we implement was proposed for Time-Dependent TSP in [36], and then reviewed and implemented as `a framework for solving realistic Vehicle Routing Problems' [22]. It consists of retaining only H, H 2 N, best, as measured by their values, states at each state space layer; we call H the depth parameter or simply depth. Naturally, if H surpasses the cardinality of the most populous state space layer, this heuristic turns into exact dynamic programming. On the other hand, xing H = 1 yields the well-known nearest neighbour heuristic.
In contrast with the exact algorithm, the base data structure chosen for the heuristic was std::map. Our intention was to implement a heuristic that is reasonably fast and has a minimum memory footprint. We did not calculate the Pareto optimum, but we gured that the linear search time associated with the use of a common array for states will not meet our de nition of \fast" and thereby took the compromise data structure. We used (std::deque) for keeping a \sorted list" of states when searching for the H best at each layer, but linear search time was not an issue here since removal only happened from the end (i.e., the worst state would go when a better one was obtained), and the addition of a state mandated re-sorting the container (done through std::sort) anyway. Another sacri ce made to diminish the memory footprint was the lack of bu ering when generating and pricing the new states, i.e., they were examined one at a time.
The heuristic is based on a restricted recurrence relation reminiscent of the exact Bellman function in (BFs). All \restricted" items in the expressions below and elsewhere are denoted by placing a tilde above the respective notation used for the \exact" items. + cj z; K ; vs 1 pr2(z); K n fjg e
o; (BgF s) vs(x; K) = e min min maxnc x; pr1(z); K j2eI(K) z2 Mej v0 x^; f?g = 0; x^; f?g 2 D0; e
The main di erence between (BFs) and (BgF s) lies in a restriction of I( ) and Mj : whereas in the exact procedure, for each layer D1; : : : ; DN , the existence of vs 1 pr2(z); Knfjg is guaranteed for all j 2 I(K) and z 2 Mj , the minimization over which yields vs(x; K), it is not always the case for a restricted procedure. It is entirely possible that some state l; K nflg 2 Ds 1 did not make it into the H best that formed Des 1 Ds 1; its (heuristic) value ve l; K nflg was not retained and could not be used in the computation of ves(x; K). Hence the need for eI(K) I(K) and Mej Mj that retain only the elements j 2 I(K) and (zin; zout) 2 Mj for which the states zout; K n fjg were among the H best that formed Des 1. Since in the restricted case the domain over which the minimization is conducted becomes smaller, we obviously have vs(x; K) ves(x; K) 8s 2 1; N for all (x; K) for which ves was calculated, thus, (BgF s) provides an upper bound Ve for the value V . The outline of the algorithm is as follows: 1. Prime the algorithm with Ge0 = G0, De0 = D0, and cfG1 = G1. The depth requirement is not imposed on De0 because all the states have the value of zero. 2. For each l 2 1; N 1 { For each K 2 Gel
Generate its feasible expansion J (K) For each j 2 J (K)
For every xout 2 M(jout) (a) Calculate v(xout; K) e (b) If jDelj = H and 9(y; K0) 2 Del : v(y; K0) > ve(xout; K),
e remove (y; K0) from Del and add (x; K) to Del. { For each (x; K) 2 Del
Let i 2 J (K) be such that x 2 Mi(out). Add K [ fig to Gel+1. 3. Calculate v(x0; 1; N ) = Ve .
e
4. Recover a route and track that conform to Ve (when substituted into (
The recovery procedure that obtains a route and track yielding at most Ve
when substituted into quality criterion (
For an upper estimate of complexity of the algorithm, let us x some more constants. Recall that we have N megalopolises plus the base. Let each megalopolis have at most m cities. Let b denote the constant time required to calculate the summary cost of exterior movement c x; pr1(z) and interior job ci pr1(z); pr2(z) ; the latter assumption is correct if all possible interior jobs are either simple or pre-calculated. Assume the time to compute the maximum of two values is also included in b.
Let us now consider, how long does it take to compute the heuristic value for a state (x; K). The generation of a feasible state K takes at most N 2 operations (see [9, Apx. 2.2]); then, we have to actually compute v(x; K), to which end we e need to consider all i 2 eI(K), of which there are never more than N , and at most m2 pairs (zin; zout) which form the set Mi; for each such pair, it takes b to calculate maxnc(x; zin) + ci(zin; zout); vjKj 1 zout; K n fig o: Thus, to calculate (BgF s) for a state (x; K), it takes at most N 2 N m2 b = N 3m2b.
There are N + 1 layers, each having at most H states, with the exception
of D0, which is not constrained by H, but still has at most N states. The last
layer DN always has a single state x0; 1; N . Thus, the computation of Ve takes
at most
(N + (N
2)H + 1)N 3m2b = O(N 4Hm2b):
(
The search for a conforming solution necessitates the examination of, in the
worst case, all (N + (N 2)H + 1) states. The examination consists of checking if
it was indeed the given state (x; K) that led to v y; K [ fjg ; thus, it is the same
as actually calculating v y; K [ fjg without the need to generate it, hence the
cost (N +(N 2)H +1)N m2b, which does not change the O-value in (
The cost of exterior movement was speci ed as Euclidean distance in X; for the sequence dependent-case, it was multiplied by a trivially sequence-dependent coe cient a(jKj) = 1 + N jKj . The cost of interior jobs was the Manhattan
N norm jj jj of movement from the \entry point" into the megalopolis to the \exit point" through its center (for two plane vectors x = (x1; x2); y = (y1; y2), the Manhattan norm is jjx yjj = jx1 y1j + jx2 y2j); this type of interior jobs is closer to the classical Generalized TSP than to Clustered TSP.
Megalopolises were modeled as equal radius disks, and the cities were placed on the circumference with equal angular distances between them (which, obviously, depended on the number of cities in the megalopolis); megalopolises were distributed randomly. Below, dimensions of the problems are encoded in the form \X{Y{Z{W", where X is the number of megalopolises, Y is the number of cities per megalopolis, Z is the number of precedence constraints, and W is either \SD" for the sequence dependence as speci ed above or \NO" for problems without sequence dependence. In the two 30{25{25{* problems, the geometry is the same (cities have the same coordinates). All data sets are available from the author on request.
Since it is the exact form of precedence constraints and not just their number
that in uences the complexity of solving the appropriate problem via dynamic
programming [46, 47, 43], we list them below, in our preferred address pairs form:
(
It is not reasonable to directly relate these dimension parameters to the top results for TSP-PC [19] since the combination of generalized nature of the problem and precedence constraints on megalopolises, in absence of a requirement to visit of all cities of a megalopolis, precludes a direct transformation into a generic (bottleneck) TSP-PC. Still, the number of cities taken into account is rather considerable for a highly constrained problem, n1 = 30 25 = 750 cities and n2 = 27 25 = 675 cities, respectively. Computation times are speci ed in the HH:MM:SS or MM:SS format.
All algorithms for all problems were encoded in C++11; exact algorithms were parallelized with the aid of the shared memory multiprocessing API OpenMP 3.0. Since a heuristic is meant to be simple to compute, we did not make a parallel implementation, yet, it is possible make such an implementation with the same means. The exact programs were run on the Uran supercomputer (for details, see http://parallel.uran.ru/node/6 [in Russian]) at IMM UrB RAS, and the heuristic was run on the author's PC (Intel Core-i4-3450, 16 GB RAM). The Uran supercomputer ran 64-bit Scientific Linux 6.4; the compiler used was GCC 4.4.7, optimization level -O2. The author's PC ran 64-bit Windows 7, the compiler used was Microsoft Visual C++ 2013 with default optimization options (\Release" con guration). 6.2
27{25{25{NO An exact solution was reported at the conference [42]. The value of the problem was V = 341:962. It took the Uran supercomputer 01:13:32 on a single core, 00:20:27 on 4 cores, and 00:10:13 on 8 cores to arrive at this conclusion; each core ran a single thread. The computation times relate as 7.197:2.002:1. The single-core to four-core relation is 3.596:1. 6.3
30{25{25{NO An exact solution was reported at the conference [42]. The value of the problem was V = 316:68. It took the Uran supercomputer 01:42:48 to arrive at this conclusion on 8 cores. 6.4
30{25{25{SD An exact solution of this problem was reported in [13]; the value of the problem was V = 376:63, and the computation took the Uran supercomputer 01:46:34 on 12 cores. 9 3 0 37 .1 0 0 5 86 .13 0 : 2 1 3 1
3 0 1 0 52 .9 0 4 19 .23 0 : 1 0 4 1 7 6 4 0 0 :002 .182 .46 1 0 1 3 2 1 :002 .192 .49 3 1 0 8 2 .4 1 0 7 2 . 1 0 7 2 . 6 3 e l 3 6 e l a 0 :02 .7 b b T 0 0 02 .76 a 8 5 7 8 3 8 7 3 8 3 03 .0 0 0 66 .05
: 1 0 9 3 1 8 :003 .072 .39 1 0 1 3 e
V H im Ve =
T Ve
02 .0 0 0 66 .56
: 1 0 9 2 7 5 6.5
Conclusion For all of the problems considered, the proposed heuristic found near-optimal solutions in a reasonable amount of time. The memory footprint of the heuristic was quite small as compared to that of the exact procedure: at most 100MB were necessary for 20000-deep solution of 30{25{25{NO, which is quite low as far as DP is concerned. Two problems were solved to optimality (30{25{25{NO and 30{25{25{SD), and one (27{25{25{NO) was solved to within 13% of optimal. The run time did not exceed 30:00.
Our intention for implementing a DP-compliant heuristic was the possible use of the latter in a Morin{Marsten branch-and-bound strategy for dynamic programming [39] to overcome the memory limitations and possibly improve the computation times. The results of experiments with the heuristic can be considered proof-of-concept: for small depth parameters (up to 250), the computation times stayed reasonably small while the result exhibited a marked improvement over the greedy algorithm (the H = 1 column) and the larger depth parameters yielded near-optimal results. Thus, a large depth parameter may yield a decent upper bound in a reasonable time. To nally implement a branch-and-bound solution, we still need a lower bound for the problem. We are not aware of lower bound algorithms that speci cally target precedence-constrained BTSP or generalized BTSP, therefore, the search has to start at general-purpose lower bound algorithms for plain BTSP; for a most recent treatment of such, refer to [32].
It is also interesting to note how the heuristic tends to stick to a locally best solution as the depth increases beyond a certain number (H = 150; H = 250 in 27{25{25{NO, H = 150 in 30{25{25{NO and 30{25{25{SD).
Acknowledgement. This work was supported by the Russian Foundation for Basic Research (project no. 13-08-00643). The author would also like to express his gratitude to the anonymous reviewers who pointed out the important shortcomings of this paper and made it possible to rectify them in the nal version.