Enhanced Filtering Algorithms for the Euclidean Traveling Salesperson Problem and its variants in Constraint Logic Programming

TPLP: Page 1–8.

© The Author(s), 2021. Published by Cambridge University Press 2021

1

doi:10.1017/xxxxx

Enhanced Filtering Algorithms for the Euclidean Traveling Salesperson Problem and its variants in Constraint Logic Programming ALESSANDRO BERTAGNON Department of Environmental and Prevention Sciences, University of Ferrara, Italy

MARCO GAVANELLI

arXiv:2608.10881v1 [cs.AI] 11 Aug 2026

Department of Engineering, University of Ferrara, Italy

Abstract The Traveling Salesperson Problem (TSP) is one of the best-known problems in computer science and arises in many engineering applications, such as smart vehicles and intelligent transportation systems. In the “Euclidean” case, each node is defined by its coordinates in the plane and distances are computed using the Euclidean metric. In the Constraint Programming (CP) literature, the Euclidean TSP is typically addressed by computing the full distance matrix and treating it as a general case; however this approach ignores the geometric information carried by the points’ coordinates. In this work, we propose new filtering algorithms, implemented in Constraint Logic Programming (CLP), that exploit such geometric information to achieve stronger constraint propagation than existing approaches. Moreover, we show how this methodology can be extended to other Euclidean variants of the TSP, including the Euclidean Generalized Traveling Salesperson Problem (EGTSP), which is relevant in practical routing and logistics applications. Experimental results demonstrate the computational advantages of the proposed approach. KEYWORDS: Euclidean Traveling Salesperson Problem, Constraint Logic Programming, Geometric Filtering Algorithms, Euclidean Generalized Traveling Salesperson Problem

1 Introduction Many engineering applications, including intelligent transportation systems, logistics, smart vehicles, and pickup-and-delivery services, require solving complex routing problems. These problems involve a vehicle or fleet of vehicles reaching a set of destinations while minimizing travel distance to save time, money, energy, and, in some cases, reduce pollution. Although each application has its own peculiarities, the problem formulation shares an omnipresent core, named in the literature the Traveling Salesperson Problem (TSP); the name comes from the problem’s most famous formulation: “A salesman has to visit a set of cities, each of which must be visited only once, and he wants to minimize the length of the tour”. More formally, given a weighted graph, the TSP requires to compute the minimum cost cycle that visits each vertex exactly once. Most routing problems are intractable, typically NP-hard. While heuristic approaches can quickly provide good solutions in some cases, systematic algorithms are preferred in

2

Bertagnon A. and Gavanelli M.

others. These algorithms can find the optimal solution, complete with a proof of optimality (in the eXplainable AI spirit) given enough time. Clearly, obtaining the provably optimal solution adds significant value. The basic TSP problem attracted the attention of researchers in different research areas, including operations research and artificial intelligence, with solution methods ranging from Integer Linear Programming, to Constraint Programming, to metaheuristics, just to name a few (Caseau and Koppstein 1993; Aguiar-Melgarejo et al. 2015; Ducomman et al. 2016; Boudreault and Quimper 2021). Beside the basic TSP problem, engineering applications involve a series of other requirements; the need to take into account such requirements spawned a plethora of extensions to the TSP, e.g., requiring each node to be visited in a time window (TSP with time windows), assuming that more than one agent (or vehicle) is available (Multiple TSP), or that not all nodes must be visited. For example, in the Generalized Traveling Salesperson Problem (GTSP) (Laporte et al. 1987; Laporte and Nobert 1983) (also known as set TSP) the nodes of the graph are partitioned into various sets, called clusters; an optimal solution of the GTSP is a minimum cost cycle that visits at least one node in each cluster. The GTSP models situations in which it is necessary to collect resources in an environment where several equivalent providers are available for each resource type, while minimizing travel time, or energy requirements from a battery. In the theoretical formulation of the TSP, the weights on the arcs can take any positive value without restrictions; however in most practical cases it is reasonable to assume that the distance function obeys the triangular inequality, where going directly from a node A to a node B cannot be more expensive than taking a detour to a third node C: such case is named metric TSP in the literature. Even more, in many industrial applications, the nodes are associated to points in the plane, they are connected through straight lines, and the distance is the length of the segment connecting the two nodes; this case is named the Euclidean Traveling Salesperson Problem (ETSP), and many instances in the TSPLIB (Reinelt 1991) set of benchmarks fall in this category. The techniques proposed in this paper are intended for complete Euclidean instances. They do not cover routing settings in which crossings are unavoidable, or in which replacing crossing connections with uncrossed ones would produce infeasible links, as may occur in real road networks. Such settings are outside the scope of the present work. From the computational complexity point of view, both the metric and the Euclidean TSP are NP-hard (Garey et al. 1976), even though, differently from the general TSP, the Euclidean TSP admits a Polynomial Time Approximation Scheme (PTAS) (Arora 1996; Mitchell 1999). Currently, the best solver for the TSP is the mathematical programming solver Concorde (Applegate et al. 2001). Concorde is based on Integer Linear Programming (ILP) techniques such as branch-and-cut and it was used to optimally solve instances with tens of thousands of nodes. However, most of the algorithms developed for solving the TSP, including Concorde, can only be used when addressing the basic TSP, whereas when there are other additional constraints, it becomes necessary to adopt more general solving techniques. A general technique developed in Artificial Intelligence to solve Constraint Optimization Problems (COPs) is constraint propagation (Mackworth 1977). The idea was so successful that new languages were embedding it as part of their operational semantics;

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

3

for example the Constraint Logic Programming (CLP) (Jaffar and Maher 1994) class of languages was defined to declaratively address constraint problems, and many languages of this class exploited constraint propagation to efficiently solve constraint problems. The CLP research area was later extended to include also imperative and object-oriented programming languages, generating the Constraint Programming (CP) research area. In CP, programmers can easily define constraint optimization problems, including side constraints, without being limited to linear or convex constraints or objective functions. This flexibility makes CP ideal for solving industrial optimization problems. Over the years, numerous solvers have been developed, and today, a variety of open-source and commercial solvers are available and widely used in both academic and commercial applications. Many works in the CP literature address the TSP (Caseau and Laburthe 1997; Kaya and Hooker 2006; Benchimol et al. 2012; Fages and Lorca 2012; Fages et al. 2016; Deudon et al. 2018); the usual way to tackle Euclidean TSPs is to compute the distance matrix and address the problem as a general TSP. On the other hand, this approach disregards all the additional information available in the ETSP: the coordinates of the points to be visited are known, and geometrical concepts (straight line segments, angles, etc.) can be defined in the Euclidean plane. An interesting research question is: “Is it possible to exploit the additional information intrinsically available in Euclidean instances in order to speed up the search?” The initial answer was presented in a conference publication Bertagnon and Gavanelli (2020), where a significant speedup was obtained in solving Euclidean instances in CLP, by exploiting efficient constraint propagation algorithms tailored for Euclidean TSP instances. This article significantly extends that work (Bertagnon and Gavanelli 2020), providing deeper insights, comprehensive analysis, and additional findings. As a first contribution, we provide a more precise description of the introduced algorithms and we extend the experimental campaign by testing the proposed algorithms on a large dataset of TSP instances, including publicly available and synthetic instances. Now, a further research question could be whether such enhanced pruning can also be applied beyond the TSP: as a second contribution, this article shows that the additional pruning based on geometric properties is not confined only to the case of the Euclidean TSP, but can also be applied to other variants of routing problems. For instance, we show how the proposed techniques can be applied to the Euclidean Generalized Traveling Salesperson Problem (EGTSP). The underlying no-crossing and convex-hull ordering properties are classical geometric results; our contribution lies in translating them into dedicated CP/CLP constraints and propagators, and in extending this approach to the EGTSP. We believe that this work could inspire researchers to extend the devised algorithms also to other routing problems, and also to inspire new propagation algorithms to exploit the additional information present in geometric instances deeply. The rest of the article is organized as follows. In Section 2 we introduce the necessary definitions and notations. In Section 3 we discuss related works. In Section 4 we provide a description of our techniques that can exploit the geometric information present in the ETSP to improve propagation. In Section 5 we show how geometric information can also be exploited in other variants of routing problems; we take EGTSP as an example. We show the computational results in Section 6 and, finally, in Section 7 we conclude.

4

Bertagnon A. and Gavanelli M. 2 Preliminaries and notation

A Constraint Satisfaction Problem (CSP) is a triple P = ⟨X , D, C⟩ where X is a set of n decision variables {x1 , x2 , . . . , xn }, D is a set of domains {Dom(x1 ), Dom(x2 ), . . . , Dom(xn )} and C is a set of constraints {c1 , c2 , . . . , cm }. Each domain Dom(xi ) is the set of all possible values that can be assigned to the variable xi . Each constraint ci consists of a pair ⟨Ri , Si ⟩ where Ri is a relation between the variables Si participating in the constraint. A solution for a CSP P is an instantiation that assigns a value to each variable in X and satisfies all the constraints C. While in a CSP all solutions are equally good, in many practical cases an objective function discriminates between solutions. A Constraint Optimization Problem (COP) is a pair ⟨P, f ⟩ where P = ⟨X , D, C⟩ is a CSP and f a function f : Dom(x1 ) × · · · × Dom(xn ) → R that associates a value to every solution of P. A solution s of P is an optimal solution of the COP ⟨P, f ⟩ if and only if there is no solution s′ of P such that f (s′ ) < f (s). In order to speed up the search of solutions, constraint programming systems alternate search with constraint propagation; each constraint is implemented by a collection of propagators. Each propagator’s objective is to observe the domains of the variables involved in the constraint and, as soon as a value is removed from the domain of a variable, it checks if some values from other domains have become unsupported and can be removed (pruned ). A propagator that is selected for execution is said to be awakened, or activated. The execution of a propagator can have one of three possible outcomes. In the first case, the propagator realizes that the constraint has no solution, e.g., because a domain becomes empty as a result of deleting some values. In this case a failure is returned and the search algorithm undoes the last choice(s) and attempts different assignments. In the second case, the propagator finds that the constraint is entailed i.e., it is satisfied whatever the values assumed by the variables among those remaining in the domains. In this case, the propagator can be deallocated. In the third case, neither of the previous two cases has occurred and so the propagator at the end of the execution is suspended waiting for an event to occur. The TSP is defined on a weighted graph G = (V, E, w), where V is a set of nodes (|V | = n), E is a set of edges, and w : E 7→ R+ . A path in G is a sequence pvs0 −vsk = vs0 es0 ,s1 vs1 . . . esk−1 ,sk vsk such that: (i) vs0 , vs1 , . . . , vsk ∈ V and are all distinct, and (ii) es0 ,s1 , es1 ,s2 , . . . , esk−1 ,sk ∈ E. Since a path is uniquely identified by the sequence of its nodes (or of its edges) in the proper order, we will often write paths as sequences of nodes to simplify the notation. Pk−1 The length of a path p is the sum of the weights of its edges: L(p) = i=0 w(esi ,si+1 ). A circuit c is a sequence obtained by appending esk ,s0 to a path pvs0 −vsk . A circuit that visits each vertex exactly once is often referred as Hamiltonian circuit. In the Euclidean case, let P = {P1 , . . . , Pn } be a set of points, where Pi = (xi , yi ). The graph associated with P is GP = (P, E P , wP ), where E P = {ei,j ≡ (Pi , Pj ) | Pi , Pj ∈ P, i ̸= j} and w(Pi , Pj ) = d(Pi , Pj ), where d is the Euclidean distance.

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

5

We use Pi Pj to denote the segment in the plane with the extremes Pi and Pj . Since in the Euclidean case every edge of a graph corresponds to a segment in the plane between the corresponding endpoints, we will often confuse the edge ei,j with the corresponding segment Pi Pj . Also, we will sometimes confuse the index i ∈ V with the corresponding ←−→ point Pi in the plane. We denote with Pi Pj the (infinite straight) line passing through points Pi and Pj , and with ∠Pi Pj Pk the counter-clockwise angle formed by the segments Pi Pj and Pj Pk with vertex in Pj from Pi to Pk . A crossing in a path p is defined as a common point Pq ∈ / P (i.e., internal to two edges) that is shared by two (or more) edges of p or a common point Pr ∈ P (i.e. an endpoint of some edge) that is shared by three (or more) edges of p. In the CP literature, three main representations have been devised for modelling the Hamiltonian circuit problem and the TSP: the permutation representation, the successor representation, and the set variable representation Benchimol et al. (2012). The last was also extended to the graph representation Dooms et al. (2005); Fages and Lorca (2012); Fages et al. (2016). In the permutation representation, n variables P osi are introduced P os = (P os1 , P os2 , . . . , P osn ) each with initial domain V ; variable P osi represents the i-th node that is visited. For example, if n = 5 and P os = (3, 5, 4, 2, 1) the corresponding tour will be (3, 5, 4, 2, 1, 3). The constraint model for the permutation representation includes an alldifferent(P os) constraint Régin (1994) on the P os array of variables, that ensures that each node is visited only once. The set variable representation is based on a set variable Puget (1992); Gervet (1993) that represents the edges that form the tour. In the successor representation, n variables Next i are defined Next = (Next 1 , Next 2 , . . . , Next n ); variable Next i represents the node that follows node i in the circuit, and its initial domain is {1, . . . , n} \ {i}. For example, if n = 5 and Next = (3, 5, 4, 2, 1) the corresponding tour will be (1, 3, 4, 2, 5, 1). The constraint model of the successor representation includes, as in the permutation representation, an alldifferent(Next) constraint on the Next array of variables, that in this case ensures that each node has exactly one incoming edge (degree constraints). The constraint model for the successor representation also includes a circuit(Next) Beldiceanu and Contejean (1994); Caseau and Laburthe (1997); Kaya and Hooker (2006) constraint (sometimes called nocycle) that avoids subtours, i.e., cycles of length less than n (subtour elimination constraints). In some cases, the constraint model includes, as redundant representation, also a set of Prev variables: Prev i represents the node that precedes node i in the circuit. Often, constraint models include redundant representations to obtain additional pruning. Both the alldifferent and the circuit constraints are already implemented in most constraint logic programming systems. Hereafter, unless otherwise specified, we refer always to the successor representation. 3 Related Work As already introduced, the Traveling Salesperson Problem is probably one of the most studied COP in the scientific literature. Many approaches proposed for its resolution have been published over time, even if we limit ourselves only to the Constraint Programming

6

Bertagnon A. and Gavanelli M.

literature. Usually, when it comes to exploit the objective function to prune the search space, CP formulations are not as effective as Integer Programming models, but they are more flexible in dealing with side constraints. To date in the CP literature, the most efficient solving methods are those that combine constraint propagation with ILP techniques. These hybrid methods yield the solving performance of CP models, under certain conditions, comparable to that obtained by ILP models. Various works propose to use relaxations of the TSP to prune suboptimal branches; the classical relaxations of the TSP are the assignment problem and the one-tree relaxation. Caseau and Laburthe (1997) propose a simple and effective rule for the circuit constraint (called nocycle). The rule consists in finding a path of mandatory edges of length at most n − 1 and removing the edge between the two endpoints of the path. They also propose to filter values based on the objective function, using the assignment-based and the spanning tree relaxations. Kaya and Hooker (2006) propose a new filtering rule based on separator graph, able to remove nonhamiltonian edges, for the circuit constraint. Francis and Stuckey (2014) consider various propagation algorithm for circuit and provide for each so-called explanations in the context of a lazy clause generation solver. Fages and Lorca (2012) show how properties of the reduced graphs (obtained by introducing a node for each Strongly Connected Components (SCC) of the original graphs) associated to Asymmetric TSPs can be used to improve the Minimum Spanning Tree (MST) relaxation. Pesant et al. (1998) address the TSP with Time Windows (TSPTW) and exploit the circuit constraint together with the MST relaxation. Focacci et al. (2002a;b) propose a filtering based on reduced costs for optimization constraints, and in particular use the assignment problem and the minimum spanning forest relaxation. The groundbreaking work in this area is that by Benchimol et al. (2012): it is the first in which a CP model was able to solve large TSP instances. In their work, they use a variety of techniques. They propose an implementation of the weighted circuit constraint that includes the Held and Karp (1970) scheme, iterates a Lagrangian relaxation to obtain a high-quality one-tree, and uses it to remove edges similarly to reduced costs filtering. It also identifies as mandatory those edges that, if removed, would increase the current lower bound over the quality of the incumbent solution. To find quickly a solution, they first run the Lin-Kernighan-Helsgaun algorithm (Lin and Kernighan 1973; Helsgaun 2000). In the experiments with asymmetric TSPs, they also use additive bounding (Fischetti and Toth 1992) to combine both the 1-tree and the assignment problem relaxations. In Fages et al. (2016), the authors further improve by casting the problem in CP(Graph) and by means of improved search heuristics (i.e., Last Conflict heuristic). More recently, Boudreault and Quimper (2021) extend the work of Benchimol et al. (2012) by proposing an enhanced CP-based Lagrangian relaxation approach. This method improves the filtering capabilities of the weighted circuit constraint by locally adjusting the Lagrangian multipliers to maximize the number of filtered values. They also develop two novel algorithms, SIMPLE and α-SETS, which provide stronger filtering and more efficient pruning of the search space. Isoart and Régin (2019) design a propagator based on the search of k-cutsets. The combination of this constraint with the Weighted Circuit Constraint (WCC) constraint has resulted in a significant reduction in the computation time.

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

7

We have already pointed out how from the computational complexity point of view, both the Metric and the Euclidean TSP are NP-hard (Garey et al. 1976); nonetheless, differently from the general TSP, the Euclidean TSP admits a PTAS (Arora 1996). A PTAS is an algorithm that, given an instance of an optimization problem and fixed a parameter ε > 0, produces in polynomial time a solution within a factor 1 + ε of the optimum. Note that the computation time can still be an exponential function of 1/ϵ. Notwithstanding these theoretically important results, solvers addressing the general TSP (e.g., Concorde) are faster in practice. Finally, concerning the GTSP, most approaches in the literature are based on ILP or metaheuristics. To the best of our knowledge, the use of CP for the GTSP has not previously been explored. Fischetti et al. studied the face structure of the polytope resulting from an ILP formulation of the problem (Fischetti et al. 1995) and proposed exact and heuristic separation algorithms for several classes of face-defining inequalities, including comb inequalities (Fischetti et al. 1997). These algorithms were used within a branch-and-cut framework. Laporte et al. (1987) addressed the problem using a subtour-elimination constraint relaxation algorithm. They also considered the asymmetric case, which they solved using a modified version of the branch-and-bound algorithm for the asymmetric TSP proposed by Carpeneto and Toth (1980). More recent works have focused on approximate methods. Pop et al. (2010) proposed a hybrid genetic algorithm for the GTSP and compared it with the random-key genetic algorithm of Snyder and Daskin (2006) and the memetic algorithm of Gutin and Karapetyan (2010). Their results showed competitive solution quality and computation time.

4 Avoiding crossings Among the geometric information that can be exploited to solve Euclidean TSPs, a wellknown result in the literature is that the optimal solution of a metric TSP (thus, also of a Euclidean TSP) in the plane cannot include two edges that cross each other (see Figure 1). We say that a set of points are all aligned if they lie on the same straight line. Theorem 4.1 (Flood (1956)). Let c∗ be an optimal tour of a metric TSP. Then, for each ei,j , ek,l ∈ c∗ such that {i, j, k, l} are all different and not all aligned, the segments Pi Pj ∩ Pk Pl = ∅. Proof. By contradiction, suppose that the optimal tour c∗ contains two segments ei,j , ek,l such that {i, j, k, l} are all different and Pi Pj ∩ Pk Pl = {Q}. Without loss of generality, let c∗ = Pi ei,j Pj pj−k Pk ek,l Pl pl−i Pi . Consider now the tour † c = Pi ei,k Pk pk−j Pj ej,l Pl pl−i Pi , where pk−j is the path pj−k reversed. The difference L(c∗ ) − L(c† ) = w(ei,j ) + w(ek,l ) − w(ei,k ) − w(ej,l ) = (d(Pi , Q) + d(Q, Pj ))+(d(Pk , Q)+d(Q, Pl ))−d(Pi , Pk )−d(Pj , Pl ) ≥ 0 applying the triangle inequality to the triangles (Pi , Pk , Q) and (Pj , Pl , Q). Considering that the points (Pi , Pk , Q) and (Pj , Pl , Q) are not aligned, the inequality becomes strict, so L(c† ) < L(c∗ ), that contradicts the fact that c∗ is optimal.

8

Bertagnon A. and Gavanelli M. Pk

Pi

b

b

pj k

pi l Q b b

Pj

Pl

Fig. 1: A self-crossing circuit. Theorem 4.1 suggests that instead of taking Pi Pj and Pk Pl , a shorter tour chooses the dotted edges Pi Pk and Pj Pl 4.1 Introducing nocrossing constraint From Theorem 4.1 follows that during the search for an optimal TSP it is possible to avoid those Hamiltonian paths that include crossing edges; this can be seen as a dominance rule (Jouglet and Carlier 2011). Based on this classical result, we propose the nocrossing constraint, which imposes that two segments in the TSP should not cross each other. In the successor representation, it is defined as follows: nocrossingi,j (Next i , Next j ) =
= (ni , nj ) ∈ Dom(Next i ) × Dom(Next j )| Pi Pni ∩ Pj Pnj ⊂ {Pi , Pj } 

(1)

Since the successor representation includes the alldifferent(Next) constraint, the variables N exti and N extj cannot assume the same value so the only type of intersection allowed by the nocrossing constraint, presented in Equation 1, is on one of the extremes i.e. PN exti ≡ Pj or PN extj ≡ Pi . The nocrossing constraint is a binary constraint, i.e., it involves exactly two decision variables, and it is also a redundant constraint as its inclusion in the model does not change the set of solutions. constraints, one for each pair For a TSP instance with n nodes, we introduce n(n−1) 2 of nodes in the graph. The nocrossing constraint can be implemented through a pair of propagators: one propagator removes values from the domain of Next i based on the values in the domain of Next j , while the other propagator propagates changes in the other direction. The nocrossing constraint propagator could be implemented naively, for example using the table constraint (Zhou 2009; Bessière and Régin 1997; Gent et al. 2007; Lecoutre 2008; Katsirelos and Walsh 2007) or the propia library (Le Provost and Wallace 1993). A table constraint consists of a table (usually a list of tuples) of values that the involved variables must, or must not, assume. However, these implementations would be inefficient, because the constraint wakes up most of the time without being able to propagate. In addition, with the table constraint one should initially compute large tables, containing, for all pairs of edges in the graph, if they cross or not. A naive propagator also tends to wake up every time a value is removed from the domain of Next i in order to remove inconsistent values from the domain of Next j . Propagating such constraint, with a naive propagator, would have the usual cost of arc con-

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

9

sistency (Mackworth 1977) propagation for a single constraint of O(d2 ) (if d is the size of the domains) in each activation of the constraint. From the definition of arc consistency and Equation 1, a value v ∈ Dom(Next j ) can be removed from Dom(Next j ) only if Pj Pv intersects all possible segments originating from Pi . A necessary condition for this is that all segments originating from Pi lie on the ←−→ same half-plane with respect to the line Pi Pj . Theorem 4.2. Let o, u ∈ Dom(Next i ) such that Po and Pu do not lie on the line l passing through Pi and Pj and are in different half-planes with respect to the line l. Then, any value k ∈ Dom(Next j ) such that Pk ̸∈ l is arc consistent with respect to the constraint nocrossingi,j (Next i , k). Proof. By contradiction, suppose there exists k ∈ Dom(Next j ) that is inconsistent with the constraint nocrossingi,j (Next i , k). Since it is inconsistent, the segment Pj Pk must cross all the segments Pi Pz such that z ∈ Dom(Next i ), so in particular it crosses both Pi Po and Pi Pu . But the intersection Io ≡ Pi Po ∩ Pj Pk lies on a different half-plane from the one that hosts Iu ≡ Pi Pu ∩ Pj Pk , so Pk lies in both half-planes, meaning that Pk ∈ l: absurd. Theorem 4.2 does not cover the case in which one of the points Po , Pu lies on the line ←−→ Pi Pj . The following proposition deals with the cases where three points are aligned. Proposition 4.3. Given a graph G whose nodes do not lie all on the same line; let a, b and c be three nodes of G such that Pc ∈ Pa Pb . Then, segment Pa Pb is not in the optimal TSP. Proof. The proof is similar to that of Theorem 4.1 as this proposition is no more than a special case. Let Pc be the point that lies on segment Pa Pb , by contradiction suppose that segment Pa Pb appears in the optimal TSP. The optimal TSP c∗ contains two points: one immediately preceding Pc , which we will denote by Pc′ , and one immediately following, denoted by Pc′′ . Let c∗ = Pa ea,b Pb pb−c′ Pc′ ec′ ,c Pc ec,c′′ Pc′′ pc′′ −a . Now, point Pc must be on segment Pc′ Pc′′ , if the points were not aligned it would be better to take segments Pc′ Pc′′ , Pa Pc and Pc Pb . Consider now the tour c† = Pa ea,c′ Pc′ pc′ −b Pb eb,c Pc ec,c′′ Pc′′ pc′′ −a . The difference L(c∗ ) − L(c† ) = w(ea,b ) + w(ec′ ,c ) − w(ea,c′ ) − w(eb,c ). At the beginning, we assumed that Pc is a point on segment Pa Pb , so we can rewrite w(ea,b ) as w(ea,c ) + w(ec,b ) obtaining L(c∗ ) − L(c† ) = w(ea,c ) + w(ec′ ,c ) − w(ea,c′ ) = d(Pa , Pc ) + d(Pc′ , Pc ) − d(Pa , Pc′ ) which is positive by the triangle inequality, contradicting the fact that segment Pa Pb can be part of c∗ . Given Proposition 4.3, it is possible to remove from the domain of each variable Next a all the values b such that the segment Pa Pb contains another node of the graph G. For the remainder of this article, we will assume that this pre-processing step has been performed before the search begins; this assumption simplifies the following discussion. Algorithm 1 sketches the algorithm of a propagator for the nocrossing constraint; it is awakened when the domain of variable Next i is reduced, and it performs propaga-

10

Bertagnon A. and Gavanelli M.

Algorithm 1 nocrossing propagator 1: function nocrossing propagator(i, Next i , j, Next j )

←−→ U nder ← select one element in Dom(Next i ) s.t. PU nder is under the line Pi Pj if there is no such element then nocrossing propagator phase 2(i, Next i , j, Next j ) ←−→ 5: else Over ← select one element in Dom(Next i ) s.t. POver is over the line Pi Pj 6: if there is no such element then 7: nocrossing propagator phase 2(i, Next i , j, Next j ) 8: else suspend waiting for either Over or U nder to be removed from Dom(Next i ) 9: function nocrossing propagator phase 2(i, Next i , j, Next j ) 10: α ← min{αx | x ∈ Dom(Next i )} 11: Let xiα be the value in Dom(Next i ) corresponding to α 12: β ← max{βx | x ∈ Dom(Next i )} 13: Let xβi be the value in Dom(Next i ) corresponding to β 14: for all y j ∈ Dom(Next j ) s.t. αyj < α do 15: if βyj > β then 16: remove y j from Dom(Next j ) 17: if |Dom(Next i )| > 1 ∧ |Dom(Next j )| > 1 then 18: suspend waiting for either xiα or xiβ to be removed from Dom(Next i ) 2: 3: 4:

tion to possibly reduce the domain of Next j ; to fully implement the constraint, another symmetric propagator would be imposed in the reverse direction (from Next j to Next i ). As long as the value j is contained in the domain of the variable Next i , according to Theorem 4.1, no propagation can take place because the point Pj lies on the same straight line as the point Pi ; for this reason, the propagator described in Algorithm 1 is imposed only when the value j is removed from Dom(Next i ) . In the first phase, lines 1-8 of Algorithm 1, the propagator is suspended and waits for ←−→ all elements in the domain of Next j to lie in the same half-plane with respect to Pi Pj . To do so, one element Over ∈ Dom(Next j ) and one U nder ∈ Dom(Next j ) that lie, ←−→ respectively, over and under the line Pi Pj are selected. If one of them does not exist, all possible segments originating from Pj lie on the same half-plane and the control passes to the next phase; otherwise, the propagator suspends waiting for one of the two Over and U nder elements to be removed from Dom(Next j ). This strategy mimics, in a sense, the idea of watched literals originally proposed in SAT solvers (Moskewicz et al. 2001) and currently used in several solvers, and significantly reduces the number of activations of the propagator. ←−→ Checking if a point Pa lies in the half-plane under or over the line Pb Pc amounts to check the sign of the (projection on the z axis) of the cross product (Pa − Pb ) × (Pc − Pb ). When all segments in the domain of the variable Next i lie on the same half-plane ←−→ with respect to the line Pi Pj , we should, for each element x ∈ Dom(N extj ), check if the segment Pj Px crosses all segments exiting from the point Pi . Naively performing the crossing check would require a number O(n2 ) of operations each time the propagator wakes up (typically an exponential number of times while exploring the search space). In order to reduce the number of crossing checks, it is possible to sort the elements Px ∈ Dom(N exti ) from smallest (denoted by the symbol α) to largest according to the ∠Pi Pj Px angle. Figure 2 shows that a segment Pj Pt (Pt ∈ Dom(Next j )) could intersect all segments originating from Pi only if the angle ∠Pi Pj Pt is smaller than α.

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

11

b

Pt b b

Pz b

↵t ↵ b

b

Pj

Pi

Fig. 2: An arrow from Px to Py means that y ∈ Dom(Next x ). Dashed lines are plotted to show the angles. In addition, if such a segment does not exist then there is no value to be removed from the Dom(Next j ) and we can suspend the propagator waiting for the value corresponding to α to be removed from Dom(Next i ). This turns out to be a necessary condition for the nocrossing propagator to remove values from the domain of Next j and is further detailed in Theorem 4.4. Theorem 4.4. For each k ∈ Dom(Next i ) ∪ Dom(Next j ), let αk = ∠Pi Pj Pk . If all the ←−→ elements q ∈ Dom(Next i ) lie on the same half-plane with respect to the line Pi Pj , then a necessary condition for a segment Pj Pt originating from Pj and reaching an element t ∈ Dom(Next j ) to cross all segments Pi Pq originating from Pi is that αt ≤ α, where α = min{αq | q ∈ Dom(Next i )}. Proof. Consider a coordinate system centered into Pj , with the abscissa pointing toward Pi and such that all the points in Dom(Next i ) have non-negative ordinate (see Figure 2). By contradiction, suppose that αt > α; we prove that there is a segment originating from Pi that does not intersect with Pj Pt . Let z ∈ Dom(Next i ) such that α = ∠Pi Pj Pz . In polar coordinates, the segment Pi Pz is seen from Pj with angles between 0 and α. All the points on the segment Pj Pt are seen under the angle αt . Since αt > α, there is no intersection between Pi Pz and Pj Pt . The condition of Theorem 4.4 is only necessary, a counterexample is shown in Figure 3, where segment Pj Pt does not cross all segments that leave the node Pi . The following theorem provides a sufficient condition for a segment to cross all segments originating from Pi . Theorem 4.5. For each k ∈ Dom(Next i ) ∪ Dom(Next j ), let βk = ∠Pk Pi Pj . Assume ←−→ all the elements q ∈ Dom(Next i ) lie on the same half-plane with respect to the line Pi Pj . Suppose there exists t ∈ Dom(Next j ) such that αt < α. Then a sufficient condition for segment Pj Pt to cross all segments Pi Pq such that q ∈ Dom(Next i ) is that βt > β, where β = max{βq | q ∈ Dom(Next i )}. Proof. By contradiction, suppose ∃k ∈ Dom(Next i ) such that Pi Pk does not intersect −−→ Pj Pt (Figure 4). Since αt < α ≤ αk , Pk and Pi lie on different sides of the ray Pj Pt .

12

Bertagnon A. and Gavanelli M. Pz b

Pk b

b

Pt b

↵ ↵t b

b

Pj

Pi

Fig. 3: The condition in Thm 4.4 is not sufficient: αt < α but Pj Pt does not cross all segments exiting from Pi . b

Pk b

Pt b

b

↵k

t

↵ ↵t

k

b

b

Pj

Pi

Fig. 4: An arrow from Px to Py means that y ∈ Dom(Next x ). Dashed lines are plotted to show the angles. In order not to have a crossing between Pi Pk and Pj Pt , both Pj and Pt must lie on −−→ the same half-plane with respect to Pi Pk . Since βj = 0, and all points in the domain of Next i are above the x axis, βj < βk , thus to be on the same half-plane, also βt < βk . But βt > β ≥ βk : contradiction. Thanks to Theorems 4.4 and 4.5 we can get a better complexity compared to the naive algorithm. Note that all angles can be pre-computed before starting the search, which avoids the computation of trigonometric functions during the search. It is also possible to pre-compute the elements in the (initial) domains of the two variables, sorted according to their angle α, and then store them in a linked list. In this way, the computation of the minimum in line 10 of Algorithm 1 is equivalent to finding the first element in the linked list that belongs to Dom(Next i ); assuming that the domain membership is checked in constant time, the minimum on line 10 can be found in O(d) amortized time on one branch of the search tree (if d is the number of elements in the domain). In the same way, the elements in Dom(Next i ) can be sorted by their angle β, and line 12 is also executed in O(d) amortized time on a branch of the search tree. The loop in lines 14

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

13

- 16 searches the linked list and stops as soon as an element is found with an angle greater than or equal to α; assuming that the removal of a domain element takes place in constant time, and since the comparison on line 15 takes a constant time, the entire loop has O(d) complexity for each activation of the propagator of phase 2 (to be compared with the O(d2 ) of the naive propagator). Since the propagator is activated when at least one element from Dom(Next i ) is removed, this propagator is woken up at most O(d) times in a branch of the search tree, which gives phase 2 a general complexity of O(d2 ) amortized time in one branch of the search tree. 4.2 Convex hull reasoning and clockwise constraint A useful consequence of Theorem 4.1 is given in the following Corollary 4.6 and is based on the concept of convex hull. The convex hull H(P) of a set of points P in a Euclidean space is the minimum convex set containing all the points. In the plane it corresponds to a convex polygon, and it is completely defined by its vertices. Corollary 4.6 (Deineko et al. (1994)). “Assuming that not all cities lie on one line, an optimal tour has the property that the cities on the boundary of the convex hull of the cities are visited in their cyclic order”. We denote with ∂H(P) = ⟨H0 , H1 , . . . , H|∂H(P)|−1 ⟩ the sequence of vertexes on the boundary of the convex hull in clockwise order. We also define a successor function S∂H(P) (Hi ) = H(i+1) mod |∂H(P)| to denote the successor of the point Hi in the sequence ∂H(P). Whenever it is clear from the context, we omit the sequence and indicate the successor of Hi with S(Hi ) instead of S∂H(P) (Hi ). To compute ∂H(P), we used the widely known Andrew’s monotone chain algorithm (Andrew 1979), with complexity O(n log n). The algorithm calculates, with complexity O(n), the upper and lower hull, which are then combined to form ∂H(P). However, this requires that the points are first sorted with respect to their x-coordinates, and also with respect to their y-coordinates in case of a tie, hence the complexity already mentioned. To simplify the exposition, the following pruning is presented in the context of symmetry breaking constraints, although similar reasoning might be performed also while not breaking symmetries. In the successor representation, the same TSP can be represented by two symmetrical solutions that differ just for the order (clockwise or counter-clockwise) in which the nodes are visited. One way to break this symmetry is to fix one direction (e.g., clockwise); in such a case, the convex hull reasoning is an efficient way to impose the clockwise order. Based on this classical cyclic-order property, we devised three ways to exploit the information about the hull for propagation. The simplest is to impose that the successor of a convex hull vertex cannot be another vertex member of ∂H(P) except for the one that immediately follows it, see Equation 2. ∀i ∈ [0, |∂H(P)| − 1], Dom(Next Hi ) ∩ ∂H(P) ⊆ {S(Hi )}

(2)

Equation 2 represents a set of unary constraints, which are in practice equivalent to reducing the initial domain of the variables, so no overhead is produced during the search (see Figure 5).

14

Bertagnon A. and Gavanelli M. Hi+1 b

Hi b

b b b b

b b

b b b b

b

Fig. 5: The successor of a convex hull vertex cannot be another vertex on the boundary of the convex hull except for the one that immediately follows it. b

b b b b b

b b

b b b b

H i ⌘ Ph

b

Fig. 6: In order to visit nodes in a clockwise order, the angle between the incoming edge and the outgoing edge of a convex hull vertex cannot be positive (it must be between −π and 0). The second way is reasoning on the angle formed by the incoming and the outgoing arcs in hull vertexes: in order to visit nodes in a clockwise order, the angle between the incoming edge and the outgoing edge of Hi ≡ Ph cannot be positive (it must be between −π and 0) or, stated otherwise, it must correspond to a right turn (see Figure 6). The easiest implementation consists of waiting until Next h becomes instantiated with a value; when the arc outgoing from Ph is fixed, the value h is removed from the domain of all other variables Next i such that the angle ∠Pi Ph PNext h would correspond to a left turn. If the constraint model also contains the Prev variables, the angle formed by the outgoing edge and any reference direction must be smaller than the angle that the incoming edge forms with the same reference direction. This produces a propagator in the same spirit of the classical less-than propagator: simply compute the minimum angle in Dom(Next h ) and remove from Dom(Prev h ) the elements associated with a smaller (or

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

15

Hi b

b

Hi+1

b

Pj

b

b b

b b

b b b b

b

Fig. 7: Each path originating from a convex hull vertex cannot reach any convex hull vertex except for the one immediately following it. Implementation is inspired by the circuit constraint (Caseau and Laburthe 1997) but performs more powerful pruning. equal) angle (Algorithm 2). A symmetric propagator takes care of the opposite direction (from Prev h to Next h ). Again, note that all angles are pre-computed before search, and the search for the minimum takes O(d) amortized time over one branch of the search tree. Algorithm 2 clockwise angle propagator Require: Ph to be on the boundary of the convex hull 1: function clockwise angle propagator(h, Next, Prev ) 2: m ← min{αv | v ∈ Dom(Next h )} 3: Dom(Prev h ) ← Dom(Prev h ) \ {i | αi ≤ m} 4: if Next h and Prev h are not instantiated with a value then 5: suspend until m is removed from Dom(Next h )

The third way results from stating that any path starting from a convex hull vertex cannot reach any other convex hull vertex except the one that directly follows it. Put it more precisely, each vertex in a path originating from a point Hi cannot reach any vertex in ∂H(P) except for S(Hi ) (see Figure 7). The propagator is imposed for each pair (Hi , S(Hi )). The implementation of this propagator is inspired by the circuit constraint (Caseau and Laburthe 1997) but performs more powerful pruning thanks to the convex hull reasoning. If a partial path has been defined starting from Hi up to a node j, and such path does not contain vertex S(Hi ), then the variable Next j cannot take any value in ∂H(P) except for S(Hi ) (Algorithm 3). If the partial path reaches the next vertex of the convex hull, the constraint is entailed. In the propagator for the circuit constraint (Caseau and Laburthe 1997), instead, from the domain of the variable Next j only the initial value (in our case, Hi ) of the path is removed. For each vertex Hi ∈ ∂H(P), the propagator is initially invoked with Start = Hi and End = S(Hi ). In the recursive calls, End remains unchanged, while Start advances along the currently instantiated path.

16

Bertagnon A. and Gavanelli M.

Algorithm 3 hull path propagator Require: PEnd to be on the boundary of the convex hull. 1: function hull path propagator(Start, End, ∂H(P), Next) 2: if Start == End then 3: return true 4: if Next Start is instantiated with a value then 5: hull path propagator(Next Start , End, ∂H(P), Next) 6: else 7: remove (∂H(P) \ {End}) from Dom(Next Start ) 8: suspend waiting for Next Start to become instantiated with a value

4.2.1 Extension of the convex hull reasoning Now that we have propagators that exploit the knowledge about the convex hull, we wish to extend their applicability to the points in the interior of the hull as well. One consequence of the absence of crossings is that the optimal TSP is a simple polygon, which is a closed polygonal chain of line segments that do not cross each other, and it divides the plane into exactly two areas: an internal and an external area. Now imagine cutting the optimal TSP with two vertical lines (parallel to the y axis): the stripe between the two lines will contain alternate internal and external areas. The borders (i.e., the parts of the circuit inside the stripe) will be visited alternately from left to right (or clockwise) and from right to left (or counter-clockwise). To exploit this informal intuition for pruning, we provide the following theorem: Theorem 4.7. Suppose that (e.g., during search) a partial path ps−e has been defined, starting in node s and ending in node e. Consider the polygon Q delimited by such path and by the segment Ps Pe , and suppose that such polygon is a simple polygon, i.e., no two edges intersect. Suppose that the partial path ps−e touches its vertexes in clockwise order. Let F ⊂ V be the set of nodes whose corresponding points lie in the interior of polygon I Q, I = F ∪ {s, e} and let ∂H(I) = ⟨H0I ≡ e, H1I , . . . , Hk−1 , Hki ≡ s⟩ be the sequence of vertexes of its convex hull, in counter-clockwise order. Suppose that the convex hull H(I) does not intersect the path ps−e , except for the endpoints Ps and Pe . Then any non self-crossing tour containing the path ps−e reaches the vertexes in ∂H(I) in the order H0I , . . . HkI . Proof. By contradiction, suppose that a non self-crossing circuit c∗ ⊇ ps−e reaches the vertexes of H I in an order different from their sequence order. W.l.o.g., suppose that I I after vertex HjI , the next vertex of H I reached by the tour c∗ is Hj+2 ; i.e., vertex Hj+1 I is not reached in the order of H I . Let pHjI −Hj+2 ⊂ c∗ be the path connecting HjI and I Hj+2 . Consider the polygon R (dotted, in Fig. 8) delimited by:

(i) the path ps−e ; (ii) the perimeter of the convex hull ∂H(I) from Pe to HjI ; I (iii) the path pHjI −Hj+2 ;

I (iv) the perimeter of the convex hull ∂H(I) from Hj+2 to Ps .

17

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

I Clearly, the point Hj+1 lies in the interior of R. In order to reach it without self∗ I . On the other hand, c∗ cannot cross crossings, c cannot pass through ps−e nor pHjI −Hj+2 I ∂H(I) between Pe and HjI (and from Hj+2 to Ps ) because, by definition of convex hull, there are no vertexes to be visited that are interior to polygon Q and that do not belong I to ∂H(I). So, c∗ cannot reach Hj+1 .

ps e b

b b

I Hj+1

I Hj+2

b

b

HjI b b b

b

b b

H(I) b

b

b

b b

b

Ps

Q

Pe b b

R

Fig. 8: From Theorem 4.7: the path ps−e is the current assignment. Polygon Q is delimited by ps−e and the (dashed) segment Pe Ps . I contains the points strictly inside Q plus Ps and Pe ; H(I) (grey in the picture) is its convex hull. Polygon R is delimited by ps−e and the dotted segments. If we find an internal hull H(I) (see Figure 8), we can then apply the previous propagators (Algorithms 2 and 3) also to ∂H(I), with the obvious care that if ps−e reaches points in clockwise order, then ∂H(I) will be reached in counter-clockwise order (and vice versa). We maintain all partial paths during the search, and we compute a convex hull for each of these paths. In this way, we are orthogonal to heuristics (we can use any search heuristic without invalidating our propagation). In the implementation, we applied the pruning on internal hulls only when the polygon Q is convex. Checking the convexity amounts to check that each turn in the path is on the same side (a right turn, on a clockwise path), it makes it easier to find when a point is inside the polygon and also allows us to avoid checking that H(I) does not intersect ps−e . The time complexity of our implementation of the extended convex hull is O(n2 ) to compute the points inside the polygon, then we use Andrew’s monotone chain (O(n log n)) (Andrew 1979) to find the hull. We currently recompute ∂H(I) from scratch after each decision.

18

Bertagnon A. and Gavanelli M. 5 Geometric reasoning for the Euclidean Generalized TSP

The propagation mechanisms developed in Section 4 are not limited to the Euclidean TSP, but can also be extended to other routing problems and variants of TSP that use the Euclidean distance as a metric. As an example, consider the Generalized Traveling Salesperson Problem (GTSP); in the GTSP, also known as set TSP, the set V of nodes of the graph is partitioned into C1 , . . . , Cm subsets where each subset of nodes Ci is called cluster and Cj ∩ Ck = ∅ for all j, k ∈ {1, . . . , m}. We use C(Pi ) to denote the cluster associated with point Pi . The objective is to compute the minimum cost cycle that visits each cluster at least once. The GTSP has a wide number of applications (Noon and Bean 1991; Laporte et al. 1996) including sequencing of computer files, routing of welfare customers through governmental agencies, design of ring networks, flexible manufacturing scheduling, airport selection and routing for courier planes, and postal routing. The Euclidean Generalized Traveling Salesperson Problem (EGTSP), similarly to the ETSP, uses Euclidean distance as its metric and the coordinates of the points on the plane are known. As in the ETSP, we employed the successor representation to model the EGTSP: the constraint model (Equation 3) involves n variables Next i . Since in the EGTSP not all nodes are visited, we take the convention proposed in Francis and Stuckey (2014) that unvisited nodes have themselves as successors (i.e. Next i = i); clearly this requires a modified version of the circuit constraint Francis and Stuckey (2014) that accepts sub-tours of length 1 (Equation 3b). The model includes the alldifferent constraint (Equation 3a) to avoid two nodes from having the same successor, and a constraint imposing that at least one node is selected for each cluster (Equation 3c). alldifferent(Next)

(3a)

set circuit(Next) X |Next i − i| ≥ 1 ∀k ∈ {1, . . . , m}

(3b) (3c)

i∈Ck

Having introduced the constraint model, we now investigate how the geometric propagation developed for the ETSP can be adapted to the EGTSP. In EGTSP we cannot directly apply the pruning introduced in Section 4.2 because the vertexes on the boundary of the convex hull may not be visited. Obviously in case such points are visited then they must be visited in the correct order, as already shown for the ETSP. In the following we will show how, under certain conditions, such pruning can also be extended to other “internal” points. To simplify the following discussion, Definition 5.1 introduces the concept of set of neighbours of a point on the boundary of the convex hull. Definition 5.1. Given Hi ∈ ∂H(P) a point Pj ∈ C(Hi ) is called a neighbour of Hi if and only if the segment Hi Pj does not intersect H(P \ C(Hi )). The set N (Hi ) is the set of all points Pj ∈ C(Hi ) that are neighbours of Hi (by definition Hi ∈ N (Hi )). In the example in Figure 9, points have a color representing the cluster they belong to. Among the points in the green cluster, only Pk is not a neighbor of Hi , because the segment Hi Pk intersects the convex hull H(P \ C(Hi )) of the non-green points.

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

19

∂H(P) !

Pk !

!

H(P \ C(Hi )) !

!

Hi

S(Hi )

Fig. 9: From Definition 5.1 Pk is not a neighbor of Hi because the segment Hi Pk intersects the convex hull H(P \ C(Hi )). Now we can extend the reasoning of the convex hull to the GTSP. Theorem 5.2. Let Pg ≡ Hi be a point in ∂H(P), and let Pb ≡ S(Hi ) be its successor in the sequence ∂H(P). Let Cg and Cb be the clusters containing Pg and Pb , respectively, with Cg ̸= Cb . Let ΦPg = H(N (Pg ) ∪ N (Pb )) be the convex hull of the neighbours of Pg and Pb , and ′ let ΦPg = H(P \ (C(Pg ) ∪ C(Pb ))) be the convex hull of all points belonging to clusters other than Cg and Cb . ′ If ΦPg ∩ ΦPg = ∅, then no point in N (Pb ) can have a point in N (Pg ) as its successor in any optimal EGTSP. Proof. Let points Pg′ and Pb′ be neighbour of point Pg and neighbour of point Pb respectively, and let c∗ be an optimal EGTSP that includes both the point Pg′ and the point Pb′ . Firstly, we want to prove that both points Pg′ and Pb′ are in ∂H(c∗ ). By definition of neighbour of a point on the boundary of the convex hull, no point in the segment Pg Pg′ is in H(P \ C(Pg )), so every point Pj of the segment is on the border ∂H((P \ C(Pg )) ∪ {Pj }). In an optimal EGTSP there is only one point for each cluster so if we have selected Pg′ there cannot be any other point belonging to cluster Cg , it follows that c∗ ⊂ H((P \ C(Pg )) ∪ {Pg′ }); since Pg′ ∈ ∂H((P \ C(Pg )) ∪ {Pg′ }), also Pg′ ∈ ∂H(c∗ ). A similar argument holds for the point Pb′ , so we proved that they are both points in ∂H(c∗ ). Suppose now, by contradiction, that Pb′ has Pg′ as its successor in c∗ (see Figure 10a). Since Pg′ is in the border ∂H(c∗ ), the angle formed by the segment Pb′ Pg′ and the segment from Pg′ to its successor must correspond to a right turn (it must be between −π and 0). It follows that the successor of Pg′ in c∗ , denoted by Pk , must be to the right ←−−→ of the straight line Pb′ Pg′ for c∗ to be optimal. According to the position of the point Pk , it is necessary to distinguish 3 different cases.

20

Bertagnon A. and Gavanelli M.

Algorithm 4 hull set propagator 1: function find neighbours(∂H(P), Next) 2: for each Hi ∈ ∂H(P) do 3: N (Hi ) ← ∅ 4: for each Pj ∈ C(Hi ) do 5: if Hi Pj does not intersect H(P \ C(Hi )) then 6: N (Hi ) ← N (Hi ) ∪ Pj

for each Hi ∈ ∂H(P) do hull set propagator(∂H(P), Next, N (Hi ), N (S(Hi ))) 9: function hull set propagator(∂H(P), Next, N (Hi ), N (S(Hi ))) 10: ΦHi ← H(N (Hi ) ∪ N (S(Hi ))) ′ 11: ΦHi ← H(P \ (C(Hi ) ∪ C(S(Hi )))) ′ 12: if ΦHi ∩ ΦHi = ∅ then 13: for each Pj ∈ N (S(Hi )) do 14: remove N (Hi ) from Dom(Next Pj ) 7: 8:

15: 16:

else ′ suspend waiting for ∂H(ΦHi ) to be modified ′

′

In the first case the point Pk lies in the interior of polygon Pb Pb Pg Pg and therefore ′ Pk ∈ ΦPg but simultaneously Pk ∈ ΦPg because C(Pk ) ̸= C(Pg ) and C(Pk ) ̸= C(Pb ) ′

contradicting the condition ΦPg ∩ ΦPg = ∅ (see Figure 10b). In the second case both segments Pk Pg′ and Pk Pb′ intersect the segment Pg Pb so it is not true that Pb is the successor of Pg in ∂H(P) (see Figure 10c). If the first two cases do not occur then either Pk Pb′ intersects Pg Pg′ or Pk Pg′ intersects Pb Pb′ . In the case where Pk Pb′ intersects Pg Pg′ it means that H(P \ C(Pg )) intersects Pg Pg′ so Pg′ by definition cannot be a neighbour of Pg (see Figure 10d). For the same reason Pb′ is not neighbour of Pb if Pk Pg′ intersects Pb Pb′ .

The implementation of a propagator exploiting Theorem 5.2 is given in Algorithm 4. In the proof of the Theorem 5.2 we also proved that neighbours of a point in ∂H(P), if selected, will be on the boundary of the convex hull of the optimal solution so it is possible to impose the clockwise constraint (see Section 4.2) on them as well. In the experimental evaluation, we applied only the reasoning on the angle formed by the incoming and the outgoing arcs in hull vertexes, using the easiest implementation, which consists of waiting until variables become instantiated with a value.

6 Experimental Evaluation To assess the effectiveness of the proposed algorithms, we devised experiments for both Euclidean TSPs and Euclidean Generalized TSPs. All algorithms are implemented in the ECLi PSe CLP language Schimpf and Shen (2012). All tests were run on ECLi PSe v. 7.0, build #48, with a time limit of 1800s on Intel® Xeon® E5-2630 v3 CPUs running at 2.4GHz, using only one core and with 2GB of reserved memory. The ECLi PSe Constraint Programming System is distributed as open-source software1 . The implementation uses 1 https://eclipseclp.org/

21

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

! !

Pg ! !

Pg ! !

!

!

Φi

!

Φi

!

Φ Pg

! !

Φ Pg

Pk !

!

Pg

Pg Pb !

Pb !

Pb

Pb (a)

(b)

! !

Pg ! !

Pk

Pg ! ! !

!

Φi

Φi ! !

! !

Pg

Pg Pb !

Pb ! Pk

Pb

Pb (c)

(d)

Fig. 10: From Theorem 5.2 Pb′ cannot have Pg′ as successor in any optimal EGTSP if ′ ΦPg ∩ ΦPg = ∅. the suspension and event-based wake-up facilities provided by ECLi PSe . These facilities are used to react to domain modifications or variable instantiation and to reschedule the propagators only when relevant changes occur. However, the filtering rules are not specific to ECLi PSe . They can be ported to other CLP or CP systems that support user-defined propagators, inspection and modification of finite domains, and event-based activation on domain changes.

22

Bertagnon A. and Gavanelli M. 6.1 Euclidean TSP

Euclidean TSP constraint model is based on the successor representation described in the preliminaries (with the circuit constraint and alldifferent Puget (1998) for improved pruning) with both Next and Prev variables, that are linked through the inverse constraint. In this section we compare different constraint models extending the one presented in the preliminaries by adding additional constraints. In the following we denote by GEO the constraint model that adds the nocrossing constraint, the removal of aligned points according to Proposition 4.3, and the clockwise constraint that implements the propagation described in Section 4.2. In order to show that the pruning we provide is not subsumed by that of state-of-theart techniques, we implemented in ECLi PSe , in the successor representation, also the Held and Karp bound with pruning based on Reduced and Marginal Costs, as proposed by (Benchimol et al. 2012) (shown with HK-RMC in the following). The purpose of this comparison is to isolate the incremental effect of the proposed geometric propagation within a common CLP implementation. As search strategies, we use the MAX-REGRET (Caseau and Laburthe 1997) and the state-of-the-art LC FIRST MAX COST (Fages et al. 2016), based on Last Conflict (Lecoutre et al. 2009). As (Benchimol et al. 2012), we also experimented injecting the upper bound given by the Lin-Kernighan-Helsgaun (LKH) (v. 2.0.7) algorithm. We experimented both with structured and with random (unstructured) instances. 6.1.1 Structured instances Structured instances model the position of real locations or sets of points from real-world problems. We consider instances taken from the TSPLIB (Reinelt 1991), the Concorde website2 and the CITIES dataset3 up to 100 nodes. These sources provide various types of instances, amongst which Euclidean ones, represented as sets of points in the plane, and geographical ones, represented as sets of points (with latitude and longitude) on the surface of the Earth. In order to have a more complete set of instances, we decided to consider also geographical instances in which the cities to be visited lie on a limited part of the geoid, so that the geographical distance can be approximated with the Euclidean distance. Table 1 reports the results for the 19 structured instances for which at least one configuration reached the optimal solution within the 1800-second time limit. Six of the 25 tested instances were omitted because all compared configurations reached the time limit. Each instance listed in the table is represented by its name, followed by the solving time (in seconds) and the number of search nodes required to find the optimal solution and prove its optimality. Aside from simple instances, the constraint models that contain geometric filtering are the ones that optimally solve the instances in the shortest time. These results show the positive interaction between the geometric filtering and that carried out by HK-RMC alone. By analysing the instances that were not solved to optimality, we found that our additional pruning is more effective during the proof of optimality, rather than on finding good initial solutions. 2 http://www.math.uwaterloo.ca/tsp/world/countries.html 3 https://people.sc.fsu.edu/

~jburkardt/datasets/cities/cities.html

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

23

Table 1: Comparing filtering algorithms on structured instances with time limit 1800s. For each instance we report total solving time and number of explored nodes to reach the optimal solution and prove its optimality.

instances uk12 burma14 ulysses16 ulysses22 wg22 bayg29 wi29 dj38 dantzig42 att48 uscap50 eil51 berlin52 kn57 wg59 st70 eil76 rat99 rd100

LC FIRST MAX COST HK-RMC GEO+HK-RMC CPU CPU time nodes time nodes 0.04 0.06 0.07 0.12 0.14 0.45 0.23 0.42 2.16 3.21 0.64 T/O 0.86 8.53 1.13 1411.63 72.56 T/O 25.81

12 14 16 22 22 30 29 38 48 66 50 52 180 59 3934 227 117

0.05 0.07 0.09 0.15 0.18 0.44 0.32 0.65 1.60 2.38 1.22 523.92 1.67 5.67 1.92 327.79 33.20 T/O 21.79

12 14 16 22 22 35 29 38 43 55 50 1384 52 120 59 985 156 107

MAX-REGRET HK-RMC GEO+HK-RMC CPU CPU time nodes time node 0.04 0.05 0.07 0.12 0.13 0.41 0.23 0.42 4.83 8.02 0.64 313.53 0.84 8.24 1.35 T/O 160.09 T/O 52.99

1 1 1 1 1 3 1 1 28 36 1 916 1 38 1 230 80

0.04 0.07 0.08 0.16 0.17 0.43 0.32 0.65 2.11 3.10 1.22 84.43 1.52 5.14 1.8 T/O 31.8 436.97 21.94

1 1 1 1 1 2 1 1 3 10 1 195 1 18 1 26 574 22

6.1.2 Random instances Random or unstructured instances consist of points with random coordinates on the plane. In the literature (Johnson and McGeoch 2007; Fischetti and Toth 1992), typically two types of instances are considered: uniform and clustered. Instances of the uniform class are obtained by placing points in the plane according to a uniform distribution, while clustered instances consist of clusters of points, whose centres are uniformly distributed in the plane. Random instances have been generated through the R-package netgen (Bossek 2015) and its functions: generateRandomNetwork, generateClusteredNetwork. The function generateRandomNetwork generates a random graph by placing points uniformly distributed in a 106 ×106 square. The function generateClusteredNetwork initially generates cluster centers by Latin Hypercube Sampling (LHS) to ensure that the clusters are placed separately from each other. Then it distributes points to the clusters according to a normal distribution having as mean vector the centers of the clusters and as variance the distance to the center of the nearest neighbor cluster. We randomly generated instances from 20 to 100 nodes in steps of 2, in both classes. For each size in the uniform class we generated 32 instances, for a total of 1312 instances. For the clustered class we tested n.cluster = 3, 5, 8, 11 as number of clusters; for each size and number of clusters we generated 16 instances, for a total of 2624 instances. For the

24

Bertagnon A. and Gavanelli M.

uniform class, 276 of the 1312 generated instances (21.0%) were excluded from the plots and aggregate statistics because all compared configurations reached the 1800-second time limit. For the clustered class, this occurred for 960 of the 2624 generated instances (36.6%). In Figure 11, we present log-log scatter plots and boxplots comparing the GEO+HKRMC model, which includes our filtering algorithm, to the HK-RMC model when solving uniform instances with the LC FIRST MAX COST search strategy. Figures 11a and 11b compare the solving time, while Figures 11c and 11d compare the search nodes. Figure 11a shows that the solving time is significantly reduced when the GEO+HK-RMC model is used compared to HK-RMC. The reduction in solving time is evident across all node ranges, with the most substantial improvements observed for larger instances. On average, there is a 70% reduction in solving time when our filtering is used. Figure 11b provides a boxplot representation of the solving time across different instances size; note that, even on a logarithmic scale, the gap between the median solving times increases with the problem size, reaching an order of magnitude for instances with 80 or more nodes, despite the saturating effect of the timeout encountered by the competitor method for larger instance sizes. Figure 11c and 11d show the number of search nodes required to reach a solution. Similar to the solving time, the GEO+HK-RMC model consistently requires fewer search nodes than the HK-RMC model, with an average reduction of 59%. Notably, this reduction holds even though, in timeout cases, the node count includes only the nodes explored up to the timeout threshold, potentially underestimating the effort required by the HK-RMC model and thus disadvantaging the GEO+HK-RMC model. In Figure 12, we present similar log-log scatter plots and boxplots comparing the GEO+HK-RMC and HK-RMC constraint models when using the MAX-REGRET search strategy. As with LC FIRST MAX COST, the GEO+HK-RMC model achieves substantial improvements in both solving time and number of search nodes. Specifically, it reduces solving time by an average of 70% and the number of search nodes by 71%, with nearly all instances lying below the diagonal in the scatter plots. Notably, several small instances are solved by the GEO+HK-RMC model using just one search node. Similar issues with timeouts apply here as well. When an instance times out, we only count the nodes explored up to the time limit. This likely underestimates how many nodes the HK-RMC model would actually need, making the comparison unfairly harder for our model. Table 2 presents the number of instances that reached the timeout limit of 1800 seconds during the experiments for various configurations of nodes and constraint models. The results indicate that for small instances (20 to 38 nodes), both algorithms performed similarly well with almost no instances reaching the timeout. As the dimension increased, the GEO+HK-RMC algorithm consistently outperformed HK-RMC, particularly with the LC FIRST MAX COST search strategy. For example, in the 90-100 node range, GEO+HK-RMC with LC FIRST MAX COST had 119 timeouts compared to 144 of HK-RMC.

Enhanced Filtering Algorithms for the ETSP and its variants in CLP # nodes

103

103

20-28

25

HK-RMC GEO+HK-RMC

40-48

102

102

50-58

CPU Time [s]

GEO+HK-RMC

30-38

60-68 70-78 80-88

101

90-100

100

101

100

10−1

10−1

10−1

100

101

102

103

20-28 30-38 40-48 50-58 60-68 70-78 80-88 90-100 # nodes

HK-RMC

(a) CPU Time

(b) CPU Time

# nodes

HK-RMC

20-28

GEO+HK-RMC

30-38 40-48

103

50-58

Search Nodes

GEO+HK-RMC

103

60-68 70-78 80-88 90-100

102

102

101

101 101

102 HK-RMC

(c) Search Nodes

103

20-28 30-38 40-48 50-58 60-68 70-78 80-88 90-100 # nodes

(d) Search Nodes

Fig. 11: Comparing solving time 11a, 11b and search nodes 11c, 11d of HK-RMC and GEO+HK-RMC models on randomly-generated Euclidean TSP instances belonging to the uniform class when LC FIRST MAX COST search strategy is used. The solving time limit was set to 1,800 seconds. Each data point in these figures corresponds to an instance. We removed instances where both models ran into timeout.

26

Bertagnon A. and Gavanelli M. # nodes

103

103

20-28

HK-RMC GEO+HK-RMC

40-48

102

102

50-58

CPU Time [s]

GEO+HK-RMC

30-38

60-68 70-78 80-88

101

90-100

100

101

100

10−1

10−1

10−1

100

101

102

103

20-28 30-38 40-48 50-58 60-68 70-78 80-88 90-100 # nodes

HK-RMC

(a) CPU Time 104

(b) CPU Time 104

# nodes 20-28

HK-RMC GEO+HK-RMC

30-38 40-48

103

50-58

Search Nodes

GEO+HK-RMC

103

60-68 70-78

102

102

80-88 90-100

101

101

100

100 100

101

102 HK-RMC

(c) Search Nodes

103

104

20-28 30-38 40-48 50-58 60-68 70-78 80-88 90-100 # nodes

(d) Search Nodes

Fig. 12: Comparing solving time 12a, 12b and search nodes 12c, 12d of HK-RMC and GEO+HK-RMC models on randomly-generated Euclidean TSP instances belonging to the uniform class when MAX-REGRET search strategy is used. The solving time limit was set to 1,800 seconds. Each data point in these figures corresponds to an instance. We removed instances where both models ran into timeout.

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

27

Table 2: Number of uniform instances that were not solved within the timeout during experiments grouped by size (shown as nodes range). The experiments compared the HK-RMC and GEO+HK-RMC constraint models, with the MAX-REGRET and LC FIRST MAX COST search strategies.

# nodes

HK-RMC

20-28 30-38 40-48 50-58 60-68 70-78 80-88 90-100

0 1 4 34 57 90 111 156

MAX-REGRET GEO+HK-RMC 0 0 1 13 29 68 78 131

LC FIRST MAX COST HK-RMC GEO+HK-RMC 0 5 12 35 51 84 105 144

0 1 3 26 34 65 68 119

Experiments on clustered instances, as shown in Figure 13, confirm that the use of geometric filtering significantly reduces the average solving time by 70% when compared to the HK-RMC model alone. Figure 13a presents a log-log scatter plot comparing the solving time of the GEO+HKRMC model with the HK-RMC model across different node ranges and cluster counts when LC FIRST MAX COST search strategy is used. Each point represents an instance, and the results show a clear reduction in solving time for the GEO+HK-RMC model. Figure 13b shows a similar scatter plot but focuses on the number of search nodes. The GEO+HK-RMC model consistently requires fewer search nodes to reach a solution, indicating more efficient search performance. Figures 13c and 13d replicate the analysis using the MAX-REGRET search strategy. Figure 13c shows the solving time comparison, while Figure 13d presents the number of search nodes. The patterns observed are similar to those with the LC FIRST MAX COST strategy, further validating the robustness of the GEO+HK-RMC model. Specifically, the GEO+HK-RMC model provided substantial reductions in the number of explored nodes: a 43% reduction with the LC FIRST MAX COST strategy and a 75% reduction with the MAX-REGRET search strategy. From Figure 13, it can also be observed that the number of clusters does not significantly affect the solving performance of the two constraint models. This suggests that the improvements provided by the geometric filtering are robust across varying cluster configurations. Table 3 presents the number of clustered instances that reached the timeout limit of 1800 seconds during the experiments, categorized by different configurations of clusters, dimensions, and constraint model. The trend concerning the instance size mirrors that observed in the uniform instances. For smaller instances (up to 68 nodes), the MAX-REGRET strategy generally performs better. However, for larger instances, the LC FIRST MAX COST strategy shows superior performance. As for the constraint models, the GEO+HK-RMC model consistently outperforms the HK-RMC model across different configurations. When analyzing the number of clusters, a clear pattern emerges: instances with a higher number of clusters are generally easier to solve, as reflected by the reduced

28

Bertagnon A. and Gavanelli M. # nodes 20-28

GEO+HK-RMC

30-38 10

2

40-48 50-58

104

n.cluster  3  5  8  11

30-38 40-48

60-68 70-78 101

80-88 90-100

10

# nodes 20-28

GEO+HK-RMC

103

10

50-58

3

60-68 70-78 80-88 90-100

102

0

10−1 10−1

101

100

101

102

103

101

102

HK-RMC

# nodes 20-28

GEO+HK-RMC

30-38 2

40-48 50-58

104

n.cluster  3  5  8  11

30-38 40-48

103

70-78 80-88 90-100

50-58

n.cluster  3  5  8  11

60-68 70-78

102

80-88 90-100

101

100

10−1 10−1

# nodes 20-28

60-68 101

104

(b) Search Nodes

GEO+HK-RMC

103

103 HK-RMC

(a) CPU Time [s]

10

n.cluster  3  5  8  11

100

100

101 HK-RMC

(c) CPU Time [s]

102

103

100

101

102

103

104

HK-RMC

(d) Search Nodes

Fig. 13: Comparing solving time and search nodes of HK-RMC and GEO+HK-RMC models on randomly-generated Euclidean TSP instances belonging to the clustered class varying the number of clusters (shown by the different colors). Figures 13a and 13b show the results obtained using the LC FIRST MAX COST search strategy while Figures 13c and 13d show the results obtained using the MAX-REGRET search strategy. The solving time limit was set to 1,800 seconds. We removed instances where both models ran into timeout. number of timeouts. This suggests that increasing the number of clusters simplifies the problem structure, making instances more tractable within the given time limit.

29

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

Table 3: Number of clustered instances that reached the timeout limit of 1800 seconds during the experiments grouped by different number of clusters and size (shown as nodes range). The experiments compared the HK-RMC and GEO+HK-RMC algorithms, with the MAX-REGRET and LC FIRST MAX COST search strategies.

MAX-REGRET HK-RMC GEO+HK-RMC

LC FIRST MAX COST HK-RMC GEO+HK-RMC

# clusters

# nodes

3

20-28 30-38 40-48 50-58 60-68 70-78 80-88 90-100

0 9 28 55 67 74 78 96

0 4 23 50 63 72 76 95

0 11 30 52 64 67 76 91

0 8 28 51 63 67 75 90

5

20-28 30-38 40-48 50-58 60-68 70-78 80-88 90-100

0 0 1 16 27 44 61 85

0 0 0 12 25 38 50 78

0 0 2 15 25 35 53 80

0 0 1 13 25 31 45 77

8

20-28 30-38 40-48 50-58 60-68 70-78 80-88 90-100

0 0 0 2 11 14 36 38

0 0 0 0 5 10 29 34

0 0 0 2 10 15 30 36

0 0 0 0 5 8 25 32

11

20-28 30-38 40-48 50-58 60-68 70-78 80-88 90-100

0 0 0 2 4 18 20 45

0 0 0 0 2 9 11 40

0 0 0 2 10 15 30 36

0 0 0 0 5 8 25 32

6.2 Euclidean Generalized TSP The instances for the GTSP in the literature (Laporte and Nobert 1983; Laporte et al. 1987; Fischetti et al. 1997) are normally created from existing structured (see Section 6.1.1) or random (see Section 6.1.2) instances. Starting from random TSP instances, we experimented with three approaches to partition points into clusters.

30

Bertagnon A. and Gavanelli M.

Once a number of clusters c is fixed, the uniform approach (inspired by (Laporte and Nobert 1983; Laporte et al. 1987)) distributes with uniform probability the points among the clusters. The clustering approach was proposed by Fischetti et al. (1997). For an instance with n points, the number of clusters is set to c = ⌈n/5⌉. Then c points, that will represent the centres of the clusters, are selected to be as far away as possible from each other. The remaining points are then associated with the nearest center. The grid approach was also proposed by Fischetti et al. (1997). Consider the minimum rectangle with sides parallel to the coordinate axes enclosing all points in the instance. The rectangle is subdivided in an L × L grid, and each grid cell containing at least one point is a cluster. The value L is determined so that each cluster contains an average number µ of nodes, let CLUSTER(H) be the number of clusters corresponding to an H × H grid, L is defined as the smallest integer such that CLUSTER(L) ≥ n/µ. As in the case of the ETSP, we compare two constraint models based on the successor representation. The first model, denoted by ECLP, includes the basic model for the EGTSP described in Equation 3. The second model, called GEO, adds the geometric pruning based on the convex hull introduced in Section 5. For all experiments we used the max-regret Caseau and Laburthe (1997) search strategy. We randomly generated instances from 20 to 50 nodes, using the three previously introduced approaches. For instances generated using the grid approach we set the parameter µ = 5. For each size and each approach we generated 16 instances, for a total of 1488 instances. Of these 1488 instances, 424 (28.5%) were excluded from the plots and aggregate statistics because all compared configurations reached the 1800-second time limit. Figure 14 compares the experimental results of the two models using a cactus plot, showing the number of solved instances as a function of the CPU time. The GEO model consistently outperforms the reference ECLP model by solving more instances within the same time limit. In addition to solving a higher number of instances, the GEO model also achieves a consistent reduction in both computation time and number of explored nodes across all instance types. If we consider instances generated with the clustering approach the average solving time with the GEO model is reduced by 76% and the average number of nodes explored is reduced by 77%. For grid type instances the results are similar as the reduction of the average solving time results 67% together with a reduction of the average number of explored nodes by 66%. For uniform type instances the average solving time is reduced by 19% and the average number of search nodes by 17%. These slightly lower results can be attributed to the distribution of points that makes it difficult to have sets of neighbours of significant size in the instances. Since the filtering algorithm relies on sets of neighbours their absence weakens the pruning. However, the results overall show that our approach is definitely useful in improving solving performance. Overall, the main practical benefit of the proposed geometric propagators is a reduction of the search space, rather than merely a lower constant cost per search node. On the random ETSP instances, the number of explored nodes is reduced by 43–75%, depending on the instance class and search strategy. This reduction is accompanied by lower solving times and fewer timeouts, especially for the larger tested instances. A similar behaviour

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

31

2,000 ECLP uniform GEO uniform ECLP grid GEO grid ECLP clustering

1,500

CPU time [s]

GEO clustering

1,000

500

0 200

220

240

260

280

300

320

340

360

380

400

# solved instances

Fig. 14: Comparing number of solved instances by ECLP and GEO models on Euclidean GTSP instances generated using different approaches to partition points into clusters. The solving time limit was set to 1,800 seconds. We removed instances where both models ran into timeout. is observed for the EGTSP, where the geometric model explores fewer nodes and solves more instances within the time limit. The complexity improvement discussed for the nocrossing propagator concerns its implementation relative to a naive propagator; we do not claim an improvement in the asymptotic complexity of solving the TSP.

7 Conclusions In this article, we showed the effectiveness of the use of geometric information in routing problems. While the most common approach in Constraint Programming and related research areas is to address the Euclidean Traveling Salesperson Problem with the same constraint model as the general TSP, we showed that the use of the additional information naturally available in many TSP instances can provide further pruning and speed up the search. In fact, in the ETSP the coordinates of the points in the plane are known, and while the common approach is to use such knowledge only to compute the distance matrix and disregarding it in the constraint model, we employed it to strengthen the constraint model with additional redundant constraints that effectively reduce the search space. This article is a revised and extended version of a conference publication Bertagnon and Gavanelli (2020), which, to the best of our knowledge, was the first work to exploit the classical no-crossing and convex-hull ordering properties through dedicated propagation mechanisms in a successor-based CP/CLP model for the Euclidean TSP.

32

Bertagnon A. and Gavanelli M.

However, while Bertagnon and Gavanelli (2020) showed the effectiveness of the idea only in the TSP problem (without any variation that can occur in real life problems), in this article we showed how the use of geometric information can be exploited also in the EGTSP, achieving a speedup of up to 4.2× in solving time and a reduction by a similar factor in the number of explored nodes. We provided a comprehensive and detailed exposition of the techniques we developed; a thorough and systematic experimentation campaign showed the effectiveness of the proposed algorithms. As no CP models for the EGTSP were found in the literature, we proposed a first constraint model; again, the addition of geometric reasoning let us reduce the average solving time to 1/3 on clustered instances and to 4/5 on uniform instances. Our results are not yet comparable with those obtained by state-of-the-art TSP solvers, such as Concorde, but on the other hand Concorde cannot solve GTSP instances. Also note that the implementation of the proposed algorithms is based on a declarative approach and on the open-source Constraint Logic Programming solver ECLi PSe . Indeed, the aim of this research was to show that exploiting geometric information can result in stronger pruning rather than proposing a very fast and efficient algorithm which would probably have required an implementation in an imperative language. We believe that the same techniques can nevertheless also be exploited in solvers implemented in imperative languages. One limitation of the current work is that it is not universally applicable to all TSP variants, such as the TSP with Time Windows, since its optimal solutions may contain crossings; but crossings avoidance could also be interesting per se in some applications. As future work, we plan to investigate extensions of the proposed techniques to other routing models in which the geometric structure and the feasible connections preserve the assumptions required by the proposed pruning.

Acknowledgements Alessandro Bertagnon and Marco Gavanelli are members of the Gruppo Nazionale Calcolo Scientifico-Istituto Nazionale di Alta Matematica (GNCS-INdAM).

Competing interests The authors declare none.

References Aguiar-Melgarejo, P., Laborie, P., and Solnon, C. A time-dependent no-overlap constraint: Application to urban delivery problems. In Michel, L., editor, Integration of AI and OR Techniques in Constraint Programming - 12th International Conference, CPAIOR 2015, Barcelona, Spain, May 18-22, 2015, Proceedings 2015, volume 9075 of Lecture Notes in Computer Science, pp. 1–17. Springer. Andrew, A. 1979. Another efficient algorithm for convex hulls in two dimensions. Information Processing Letters,, 9, 5, 216–219. Applegate, D., Bixby, R. E., Chvátal, V., and Cook, W. J. TSP cuts which do not conform to the template paradigm. In Jünger, M. and Naddef, D., editors, Computational

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

33

Combinatorial Optimization, Optimal or Provably Near-Optimal Solutions 2001, volume 2241 of Lecture Notes in Computer Science, pp. 261–304. Springer. Arora, S. Polynomial time approximation schemes for Euclidean TSP and other geometric problems. In Proceedings of 37th Conference on Foundations of Computer Science 1996, pp. 2–11. Beldiceanu, N. and Contejean, E. 1994. Introducing global constraints in CHIP. Math. Comput. Model.,, 20, 12, 97–123. Benchimol, P., van Hoeve, W. J., Régin, J., Rousseau, L., and Rueher, M. 2012. Improved filtering for weighted circuit constraints. Constraints,, 17, 3, 205–233. Bertagnon, A. and Gavanelli, M. Improved filtering for the euclidean traveling salesperson problem in CLP(FD). In The Thirty-Fourth AAAI Conference on Artificial Intelligence, AAAI 2020, The Thirty-Second Innovative Applications of Artificial Intelligence Conference, IAAI 2020, The Tenth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2020, New York, NY, USA, February 7-12, 2020 2020, pp. 1412–1419. AAAI Press. Bessière, C. and Régin, J. Arc consistency for general constraint networks: Preliminary results. In Proceedings of the Fifteenth International Joint Conference on Artificial Intelligence, IJCAI 97, Nagoya, Japan, August 23-29, 1997, 2 Volumes 1997, pp. 398–404. Morgan Kaufmann. Bossek, J. 2015. netgen: Network generator for combinatorial graph problems. R package version,, 1. Boudreault, R. and Quimper, C.-G. Improved cp-based lagrangian relaxation approach with an application to the tsp. In Zhou, Z.-H., editor, Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21 2021, pp. 1374–1380. International Joint Conferences on Artificial Intelligence Organization. Main Track. Carpeneto, G. and Toth, P. 1980. Some new branching and bounding criteria for the asymmetric travelling salesman problem. Management Science,, 26, 7, 736–743. Caseau, Y. and Koppstein, P. A rule-based approach to a time-constrained traveling salesman problem. In 2nd Int. Symp. on AI and Mathematics 1992 1993. Bellcore Technical Memorandum. Caseau, Y. and Laburthe, F. Solving small TSPs with constraints. In Naish, L., editor, Logic Programming, Proc. of the Fourteenth International Conference on Logic Programming, Leuven, Belgium, July 8-11, 1997 1997, pp. 316–330. MIT Press. Deineko, V. G., van Dal, R., and Rote, G. 1994. The convex-hull-and-line traveling salesman problem: A solvable case. Inf. Process. Lett.,, 51, 3, 141–148. Deudon, M., Cournut, P., Lacoste, A., Adulyasak, Y., and Rousseau, L. Learning heuristics for the TSP by policy gradient. In van Hoeve, W. J., editor, Integration of Constraint Programming, Artificial Intelligence, and Operations Research - 15th International Conference, CPAIOR 2018, Delft, The Netherlands, June 26-29, 2018, Proceedings 2018, volume 10848 of LNCS, pp. 170–181. Springer. Dooms, G., Deville, Y., and Dupont, P. CP(Graph): Introducing a graph computation domain in constraint programming. In van Beek, P., editor, Principles and Practice of Constraint Programming - CP 2005, 11th International Conference, CP 2005, Sitges, Spain, October 1-5, 2005, Proceedings 2005, volume 3709 of Lecture Notes in Computer Science, pp. 211–225. Springer. Ducomman, S., Cambazard, H., and Penz, B. Alternative filtering for the weighted circuit constraint: Comparing lower bounds for the TSP and solving TSPTW. In Schuurmans, D. and Wellman, M. P., editors, Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence, February 12-17, 2016, Phoenix, Arizona, USA 2016, pp. 3390–3396. AAAI Press. Fages, J. and Lorca, X. 2012. Improving the asymmetric TSP by considering graph structure. CoRR,, abs/1206.3437.

34

Bertagnon A. and Gavanelli M.

Fages, J., Lorca, X., and Rousseau, L. 2016. The salesman and the tree: the importance of search in CP. Constraints,, 21, 2, 145–162. Fischetti, M., González, J. J. S., and Toth, P. 1995. The symmetric generalized traveling salesman polytope. Networks,, 26, 2, 113–123. Fischetti, M., González, J. J. S., and Toth, P. 1997. A branch-and-cut algorithm for the symmetric generalized traveling salesman problem. Oper. Res.,, 45, 3, 378–394. Fischetti, M. and Toth, P. 1992. An additive bounding procedure for the asymmetric travelling salesman problem. Math. Program.,, 53, 173–197. Flood, M. M. 1956. The traveling-salesman problem. Operations Research,, 4. Focacci, F., Lodi, A., and Milano, M. 2002a. Embedding relaxations in global constraints for solving TSP and TSPTW. Ann. Math. Artif. Intell.,, 34a, 4, 291–311. Focacci, F., Lodi, A., and Milano, M. 2002b. A hybrid exact algorithm for the TSPTW. INFORMS Journal on Computing,, 14b, 4, 403–417. Francis, K. G. and Stuckey, P. J. 2014. Explaining circuit propagation. Constraints,, 19, 1, 1–29. Garey, M. R., Graham, R. L., and Johnson, D. S. Some NP-complete geometric problems. In Proceedings of the Eighth Annual ACM Symposium on Theory of Computing 1976, STOC ’76, pp. 10–22, New York, NY, USA. ACM. Gent, I. P., Jefferson, C., Miguel, I., and Nightingale, P. Data structures for generalised arc consistency for extensional constraints. In Proceedings of the Twenty-Second AAAI Conference on Artificial Intelligence, July 22-26, 2007, Vancouver, British Columbia, Canada 2007, pp. 191–197. AAAI Press. Gervet, C. New structures of symbolic constraint objects: sets and graphs. In Third Workshop on Constraint Logic Programming (WCLP’2003) 1993. Gutin, G. Z. and Karapetyan, D. 2010. A memetic algorithm for the generalized traveling salesman problem. Nat. Comput.,, 9, 1, 47–60. Held, M. and Karp, R. M. 1970. The traveling-salesman problem and minimum spanning trees. Operations Research,, 18, 1138–1162. Helsgaun, K. 2000. An effective implementation of the Lin-Kernighan traveling salesman heuristic. European Journal of Operational Research,, 126, 1, 106–130. Isoart, N. and Régin, J.-C. Integration of structural constraints into TSP models. In International Conference on Principles and Practice of Constraint Programming 2019, pp. 284–299. Springer. Jaffar, J. and Maher, M. J. 1994. Constraint logic programming: A survey. J. Log. Program.,, 19/20, 503–581. Johnson, D. S. and McGeoch, L. A. Experimental analysis of heuristics for the stsp. In The traveling salesman problem and its variations 2007, pp. 369–443. Springer. Jouglet, A. and Carlier, J. 2011. Dominance rules in combinatorial optimization problems. European Journal of Operational Research,, 212, 3, 433–444. Katsirelos, G. and Walsh, T. A compression algorithm for large arity extensional constraints. In Bessiere, C., editor, Principles and Practice of Constraint Programming - CP 2007, 13th International Conference, CP 2007, Providence, RI, USA, September 23-27, 2007, Proceedings 2007, volume 4741 of Lecture Notes in Computer Science, pp. 379–393. Springer. Kaya, L. G. and Hooker, J. N. A filter for the circuit constraint. In Benhamou, F., editor, Principles and Practice of Constraint Programming - CP 2006 2006, volume 4204 of Lecture Notes in Computer Science, pp. 706–710. Springer. Laporte, G., Asef-Vaziri, A., and Sriskandarajah, C. 1996. Some applications of the generalized travelling salesman problem. The Journal of the Operational Research Society,, 47, 12, 1461–1467.

Enhanced Filtering Algorithms for the ETSP and its variants in CLP

35

Laporte, G., Mercure, H., and Nobert, Y. 1987. Generalized travelling salesman problem through n sets of nodes: the asymmetrical case. Discrete Applied Mathematics,, 18, 2, 185–197. Laporte, G. and Nobert, Y. 1983. Generalized travelling salesman problem through n sets of nodes: An integer programming approach. INFOR: Information Systems and Operational Research,, 21, 1, 61–75. Le Provost, T. and Wallace, M. 1993. Generalized constraint propagation over the CLP scheme. J. Log. Program.,, 16, 3, 319–359. Lecoutre, C. Optimization of simple tabular reduction for table constraints. In Stuckey, P. J., editor, Principles and Practice of Constraint Programming, 14th International Conference, CP 2008, Sydney, Australia, September 14-18, 2008. Proceedings 2008, volume 5202 of Lecture Notes in Computer Science, pp. 128–143. Springer. Lecoutre, C., Saı̈s, L., Tabary, S., and Vidal, V. 2009. Reasoning from last conflict(s) in constraint programming. Artif. Intell.,, 173, 18, 1592–1614. Lin, S. and Kernighan, B. W. 1973. An effective heuristic algorithm for the traveling-salesman problem. Operations Research,, 21, 2, 498–516. Mackworth, A. K. 1977. Consistency in networks of relations. Artif. Intell.,, 8, 1, 99–118. Mitchell, J. S. B. 1999. Guillotine subdivisions approximate polygonal subdivisions: A simple polynomial-time approximation scheme for geometric TSP, k-MST, and related problems. SIAM J. Comput.,, 28, 4, 1298–1309. Moskewicz, M. W., Madigan, C. F., Zhao, Y., Zhang, L., and Malik, S. Chaff: Engineering an efficient SAT solver. In Proceedings of the 38th Annual Design Automation Conference 2001, DAC ’01, pp. 530–535, New York, NY, USA. ACM. Noon, C. E. and Bean, J. C. 1991. A lagrangian based approach for the asymmetric generalized traveling salesman problem. Operations Research,, 39, 4, 623–632. Pesant, G., Gendreau, M., Potvin, J., and Rousseau, J. 1998. An exact constraint logic programming algorithm for the traveling salesman problem with time windows. Transportation Science,, 32, 1, 12–29. Pop, P. C., Matei, O., and Sabo, C. A new approach for solving the generalized traveling salesman problem. In Blesa, M. J., Blum, C., Raidl, G. R., Roli, A., and Sampels, M., editors, Hybrid Metaheuristics - 7th International Workshop, HM 2010, Vienna, Austria, October 1-2, 2010. Proceedings 2010, volume 6373 of Lecture Notes in Computer Science, pp. 62–72. Springer. Puget, J. A fast algorithm for the bound consistency of alldiff constraints. In Mostow, J. and Rich, C., editors, Proceedings of the Fifteenth National Conference on Artificial Intelligence and Tenth Innovative Applications of Artificial Intelligence Conference, AAAI 98, IAAI 98, July 26-30, 1998, Madison, Wisconsin, USA 1998, pp. 359–366. AAAI Press / The MIT Press. Puget, J.-F. Pecos: a high level constraint programming language. In Proceedings of SPICIS 1992, volume 92. Régin, J. A filtering algorithm for constraints of difference in CSPs. In Hayes-Roth, B. and Korf, R. E., editors, Proceedings of the 12th National Conference on Artificial Intelligence, Seattle, WA, USA, July 31 - August 4, 1994, Volume 1. 1994, pp. 362–367. AAAI Press / The MIT Press. Reinelt, G. 1991. TSPLIB - A traveling salesman problem library. INFORMS Journal on Computing,, 3, 4, 376–384. Schimpf, J. and Shen, K. 2012. Ecli pse - from LP to CLP. TPLP,, 12, 1-2, 127–156. Snyder, L. V. and Daskin, M. S. 2006. A random-key genetic algorithm for the generalized traveling salesman problem. Eur. J. Oper. Res.,, 174, 1, 38–53. Zhou, N.-F. Encoding table constraints in CLP(FD) based on pair-wise AC. In Hill, P. M. and Warren, D. S., editors, Logic Programming, 25th International Conference, ICLP 2009,

36

Bertagnon A. and Gavanelli M.

Pasadena, CA, USA, July 14-17, 2009. Proceedings 2009, volume 5649 of Lecture Notes in Computer Science, pp. 402–416. Springer.

Related documents

• A characterization of efficiently compilable constraint languages#927083
• Knowledge by constraint#740428
• Solving Fuzzy Satisfiability via Mixed-Integer Non-Linear Programming#785147
• Solving Fuzzy Satisfiability via Mixed-Integer Non-Linear Programming#783557
• Introduction to QUDO, Tensor QUDO and HOBO formulations: Qudits, Equivalences, Knapsack Problem, Traveling Salesman Problem and Combinatorial Games#189316