Geometric reasoning on the euclidean traveling salesperson problem in answer set programming 1

Abstract

The Traveling Salesperson Problem (TSP) is one of the best-known problems in computer science. Many instances and real world applications fall into the Euclidean TSP special case, in which each node is identified by its coordinates on the plane and the Euclidean distance is used as cost function. It is worth noting that in the Euclidean TSP more information is available than in the general case; in a previous publication, the use of geometric information has been exploited to speedup TSP solving for Constraint Logic Programming (CLP) solvers. In this work, we study the applicability of geometric reasoning to the Euclidean TSP in the context of an ASP computation. We compare experimentally a classical ASP approach to the TSP and the effect of the reasoning based on geometric properties. We also compare the speedup of the additional filtering based on geometric information on an ASP solver and a CLP on Finite Domain (CLP(FD)) solver.

Keywords

Answer set programming euclidean traveling salesperson problem experimental comparison of ASP and CLP(FD)

1 Introduction

The Traveling Salesperson Problem (TSP) is a classical COP that, given a graph with non-negative weights on the edges, involves finding the minimal cost cycle that visits each node exactly once. Many instances and real world applications fall into a special case called Euclidean TSP, as also evidenced by the many Euclidean instances present in the famous TSPLIB [32] benchmark set. In the Euclidean TSP each node is identified by its coordinates in the Euclidean plane, and the weight is the Euclidean distance.

The usual way to tackle Euclidean TSPs is to compute the distance matrix and address the problem as a general TSP. However, in the Euclidean TSP more information is available than in the general case: the coordinates of the nodes are available and geometric concepts, like angles and straight lines, can be defined in the Euclidean plane. In a previous publication [5], we showed that the use of geometric information can be exploited to speedup TSP solving in CLP solvers. The work started from the observation that, due to the triangular inequality, in a Euclidean TSPs two edges cannot cross each other, since eliminating the crossing would result in another cycle of shorter length. This observation was then extended and exploited in the development of efficient algorithms that reduce the search space.

Beside CLP, ASP is another logic programming paradigm that due to the expressiveness of the language and the availability of efficient solving systems [10 , 34] is used in advanced applications.

In this paper we study the applicability of geometric reasoning to the Euclidean TSP in the context of ASP. We propose an encoding of Euclidean TSP as an ASP program with reasoning based on geometric properties, aiming to eliminate crossings also by leveraging properties associated with the convex hull of the point set. We present and empirically evaluate different ASP encodings for computing the convex hull of sets of points.

To assess the performance and applicability of our proposed techniques we compare how solving time is affected when different types of geometric information are exploited in the encoding. We show that the obtained speedup is often between 1 and 3 orders of magnitude. The experiments involve both randomly generated instances and real-life instances such as those from the aforementioned TSPLIB [32]. On real-life instances, we compare the solution quality after running both methods for the same time, and obtained that the tour length obtained by applying geometric reasoning was 75% with respect to that obtained by the basic ASP approach.

Finally, we compare the speedup of the additional filtering based on geometric information on an ASP solver and a CLP on Finite Domain (CLP(FD)) solver; in both cases the speedup is computed with respect to the usual approach of computing the distance matrix and then solving the TSP with the usual techniques (both in ASP and in CLP(FD)).

2 Preliminaries

An ASP program Π consists of a finite set of clauses interpreted in the Stable Model semantics [22]. Each clause in an implication H ← B. The full ASP syntax is reported in [9]. The head H of the clause can be an atom of the form a (t₁, …, t_n) (where t₁, … t_n are terms), a choice atom {a (t₁, …, t_n)} or an aggregate atom. The body B is a set of literals of the form either a or not a where again a can be an atom or an aggregate atom. An aggregate atom can be A {t₁, …, t_m : c₁, …, c_m} ∘ n where A can be #sum, #min, #max or #count, and ∘ can be a relational operator.

Clauses with empty body are called facts, while clauses with empty head are called ICs, and their head is intended to be false.

An interpretation I of a program Π is a subset (not necessarily proper) of the set of atoms occurring in Π. The atoms in I are considered true in the interpretation, while all remaining atoms are false.

The reductΠ^I of a program Π with respect to the interpretation I is a set of clauses defined as follows. For each rule r ∈ Π containing a negative literal not a such that a ∈I, r does not belong to Π^I. For all other rules r ∈ Π, all negative literals are removed. Clearly, the reduct Π^I is a positive program. An interpretation I is a stable model of program Π if it is a minimal model of Π^I.

Since each node in a graph of a Euclidean TSP is associated with a point in a plane, in the rest of the paper we will often use the terms “point” and “node” indistinctly.

3 Filtering techniques for the Euclidean TSP in ASP

As stated in previous sections, Euclidean TSP instances include the coordinates on the plane of the points; such coordinates in ASP can be represented by facts point(I,X,Y) where I is a numerical identifier of the point and X,Y are its coordinates in the Euclidean plane. Instances also include cost(A,B,C) facts where C is the Euclidean distance between points A and B. A cycle on the graph, so a solution, will be represented by a set of cycle(A,B) atoms.

The basic Euclidean TSP encoding in ASP is presented in Listing 1. This encoding includes the degree constraint that ensures that, in the solution, the degree of each node is equal to two. Line 1 ensures that each point has exactly one outgoing edge while line 2 ensures that each point has exactly one incoming edge. Lines 4–6 state that starting from the point with smallest I all other points can be reached through the cycle; this ensures that all points are visited. Line 7 is the objective function; we want to minimize the cost of the cycle.

{ cycle(A,B): cost(A,B,_) }=1:- point(A,_,_).

{ cycle(A,B): cost(A,B,_) }=1:- point(B,_,_).

reached(A):-

A = #min{ I: point(I,_,_) }.

reached(B):- cycle(A,B), reached(B).

:- point(B,_,_), not reached(B).

#minimize{ C,A,B: cycle(A,B), cost(A,B,C) }.

Listing 1: Euclidean TSP encoding in ASP

3.1 Avoiding intersections

The following is a well-known result in the literature.

Theorem 1. (Flood [16]). The optimal solution of a metric TSP in the plane cannot include two edges that cross each other.

Intuitively, in Fig. 1 instead of taking segments $\bar{P_{i} P_{j}}$ and $\bar{P_{k} P_{l}}$ , a shorter tour chooses the gray edges $\bar{P_{i} P_{k}}$ and $\bar{P_{l} P_{j}}$ .

Fig. 1

A self-crossing circuit. (Picture created with ASPECT [7]).

The Euclidean TSP is a special case of the metric TSP since the Euclidean distance satisfies the triangular inequality; from Theorem 1, one can avoid, during the search for an optimal Euclidean TSP, those paths that include crossing edges.

Avoiding crossings is particularly simple and declarative in ASP; in Listing 2, line 1 imposes that a pair of segments in the TSP should not cross each other (except, possibly, in one endpoint). The cross(A,B,C,D) predicate declaratively states that segment $\bar{AB}$ and segment $\bar{CD}$ cross each other. To verify the intersection of two segments given their endpoints, we implemented the Faster Line Segment Intersection algorithm proposed by F. Antonio [1]. The idea of this algorithm is to represent each point as a 2D vector (A = (x_a, y_a)). Now a generic point P on $\bar{AB}$ can be represented parametrically by a linear combination of A and B as follows: $P = α A + (1 - α) B$ where α is in the interval [0, 1] when representing a point on the segment $\bar{AB}$ . Finding the crossing point Q of two segments $\bar{AB}$ and $\bar{CD}$ reduces to solving the following system of equations: ${\begin{matrix} Q & = A + r (B - A) \\ Q & = C + s (D - C) \end{matrix}$ whose solution is $\begin{matrix} r = \frac{{Num}_{r}}{Den} \\ s = \frac{{Num}_{s}}{Den} \end{matrix}$ where $\begin{matrix} {Num}_{r} = (y_{a} - y_{c}) (x_{d} - x_{c}) - (x_{a} - x_{c}) (y_{d} - y_{c}) \\ Den = (x_{b} - x_{a}) (y_{d} - y_{c}) - (y_{b} - y_{a}) (x_{d} - x_{c}) \\ {Num}_{s} = (y_{a} - y_{c}) (x_{b} - x_{a}) - (x_{a} - x_{c}) (y_{b} - y_{a}) \end{matrix}$ and where the coordinates of a point P are named (x_p, y_p). However, we are not interested in explicitly calculating the intersection point, so we do not need r and s explicitly, it is sufficient to verify that r and s are in the range ] 0, 1 [ (it is easy to verify that the extreme values 0 and 1 correspond to either parallel segments or intersection in one of the extremes of one segment).

Clearly, in order to have r > 0 and s > 0, the numerators and denominator must all have the same sign (lines 10-12 in Listing 2); in order to have r < 1 and s < 1, the absolute value of the nominator must be less than the absolute value of the denominator (lines 13-14).

:- cycle(A,B), cycle(C,D), cross(A,B,C,D).

cross(A,B,C,D):-

point(A,Ax,Ay),

point(B,Bx,By),

point(C,Cx,Cy),

point(D,Dx,Dy),

Den =

(Bx-Ax)*(Dy-Cy)-(By-Ay)*(Dx-Cx),

NumR =

(Ay-Cy)*(Dx-Cx)-(Ax-Cx)*(Dy-Cy),

NumS =

(Ay-Cy)*(Bx-Ax)-(Ax-Cx)*(By-Ay),

Den != 0,

NumR * Den >0,

NumS * Den >0,

|NumR| < |Den|,

|NumS| < |Den|.

Listing 2: Constraint for avoiding crossing edges

It is interesting to notice that currently available intelligent grounders (such as gringo [20]) are able to compile the cross/4 predicate in Listing 2 into a series of facts, and the integrity constraint in line 1 is grounded into a series of denials that forbid the pairs of edges that intersect. In their turn, these denials help the solver reducing the search space without adding new choices to the code.

3.2 Reasoning on the convex hull

A useful consequence of Theorem 1 is based on the concept of convex hull and is given in Corollary 1. The convex hull 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. We denote with $H (P) = 〈 H_{0}, H_{1}, \dots, H_{| H (P) | - 1} 〉$ the sequence of nodes on the boundary of the convex hull.

Corollary 1. (Deineko, van Dal, and Rote [13]). Assuming that not all nodes of the graph lie on the same line, nodes on the boundary of the convex hull are always visited in their cyclic order in an optimal TSP.

Corollary 1 can be exploited in several ways to reduce the search space; the main idea is to define an order of visit (either clockwise or counterclockwise) of the tour, and impose that in all points belonging to the border of the convex hull such order is preserved. In a sense this type of pruning can be thought as a symmetry breaking constraint, and in the experimental evaluation it will be compared with another symmetry breaking technique that imposes one order of visit of the tour.

The points on the boundary of the convex hull will be identified by means of the hull(X) predicate. The counterclockwise visiting order of the hull will be encoded through convex_hull_next(A,B) facts, which state that point B follows point A on the boundary of the convex hull.

In [5] we devised three ways to exploit the information about the hull for propagation. In the following we report them together with the corresponding ASP encoding. The first constraint, implemented in Listing 3, 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.

:- not convex_hull_next(X,Y), cycle(X,Y), hull(X), hull(Y).

Listing 3: 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.

This integrity constraint imposes that edges connecting far away nodes in the border of the convex hull cannot be taken, and are immediately propagated by the unit propagation of the solver with no need of search.

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 counterclockwise order, the angle between the incoming edge and the outgoing edge of a node H_i on the boundary of the convex hull cannot be negative (it must be between 0 and π) or, stated otherwise, it must correspond to a left turn.

:- hull(X), cycle(From,X), cycle(X,To), turn_right(From,X,To).

turn_right(From,P,To):- point(From,X,Y), point(P,X0,Y0), point(To,X1,Y1), (Y-Y0)*(X1-X0)-(X-X0)*(Y1-Y0)<0.

Listing 4: In order to visit nodes in a counterclockwise order, the angle between the incoming edge and the outgoing edge of a convex hull vertex cannot be negative.

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 H_i on the boundary of the convex hull cannot reach any vertex in $H (P)$ except for $H_{(i + 1 \mod | H (P) |)}$ .

path_hull(A,B):- hull(A), cycle(A,B), not hull(B).

path_hull(A,C):- hull(A), cycle(B,C), not hull(C), path_hull(A,B).

:- not convex_hull_next(A,C), path_hull(A,B), cycle(B,C), hull(C).

Listing 5: Each path originating from a convex hull vertex cannot reach any convex hull vertex except for the one immediately following it.

Predicate path_hull(A,B) is true if from a vertex A on the border of the convex hull there is a path connecting it to an internal vertex B (i.e., a vertex that is not on the boundary of the convex hull) passing only through internal vertexes. The integrity constraint in line 6 imposes that such path cannot terminate in any vertex on the border of the convex hull, except the immediate successor of vertex A.

3.2.1 Computing the convex hull

Several algorithms exist to compute the border of the convex hull of a set of vertices; however we preferred to compute it with a declarative program directly in ASP. We can declaratively state that a vertex is on the border of the convex hull if it is not internal to any triangle having as vertexes any three nodes of the graph (line 1 in Listing 6). In order to check if a point P is internal to some triangle (predicate inside in Listing 6), it is enough to check if it stands always on the same side with respect to the line segments on the sides of the triangle; on its turn, this can be easily defined leveraging on the previously defined predicate turn_right.

This approach computes the set of vertices on the border of the hull in O (n⁴), when n is the number of vertexes; this is much higher than other available algorithms (e.g., Andrew’s monotone chain has O (n log n) complexity), but it is worth noting that all the computation can be performed by the grounder, and indeed gringo provides a ground set of facts so that during the search for an answer set it is not necessary to compute $H (P)$ . In future work we will explore other algorithms to compute the convex hull.

Once the points belonging to the boundary of the hull have been found, it is necessary to find what is their cyclic order. Given two points A and B on the boundary of the convex hull, point B is successor to point A, chosen the counterclockwise direction, if there are no other points to the “right” of segment $\bar{AB}$ . (see line 3).

Note that all the computation in Listing 6, concerning the calculation of the convex hull, is performed by the grounder therefore it does not affect solving performance.

hull(P):- not inside(P), point(P,_,_).

inside(P):- point(P1,_,_), point(P2,_,_), point(P3,_,_), turn_right(P,P1,P2), turn_right(P,P2,P3), turn_right(P,P3,P1).

convex_hull_next(A,B):- hull(A), not turn_right(A,B,_), hull(B), A != B.

Listing 6: Calculation of convex hull in ASP

3.2.2 Jarvis march

The program presented in Section 3.2.1 that computes the points in the border of the convex hull is highly declarative and stratified, and it is completely solved by efficient grounders such as gringo. On the other hand, it has clearly O (n⁴) complexity (if n is the number of points), and there exist more efficient algorithms.

We implemented in ASP a version inspired by the Jarvis march [26]. In a nutshell, the Jarvis march can be described as follows. First, a point on the border of the hull is found (for example, the point having smallest x coordinate, breaking possible ties by taking the point with smallest y coordinate). Afterwards, the following inductive algorithm is executed: starting from the last point A added to the set $H (P)$ , the point B that corresponds to a leftmost turn from A, i.e., such that there is no point C that is left with respect to the segment $\bar{AB}$ , also belongs to the border of the hull. In case of ties (i.e., in case there are three aligned points in the border of the hull), we take the point in the border that is closest to A. Such algorithm has complexity O (hn), where h is the number of points in the border of the hull.

The ASP implementation inspired by the Jarvis march is shown in Listing 7. The point with smallest x is declared a point in the border of the hull; the following ones are found (predicate convex_closure_next) by the inductive definition described earlier. The complexity is not the same as the Jarvis March implemented in an imperative language, but it is improved with respect to the implementation in Section 3.2.1: the fact that in lines 4 and 5 the point A is limited to be on the border of the hull (hull (A)) reduces complexity to O (n²h). On the other hand, the program in Listing 7 is not stratified, and gringo does not completely solve it.

hull(Leftmost):- point(Leftmost,Xl,Yl), not exists_point_smaller_x(Leftmost).

hull(B):- convex_closure_next(A,B).

convex_closure_next(A,B):- hull(A), not exists_external_point(A,B), point(B,_,_), A!=B.

exists_external_point(A,B):- hull(A), turn_right(C,A,B).

exists_external_point(A,B):- hull(A), aligned(A,C,B).

aligned(A,B,C):- point(A,Xa,Ya),

point(B,Xb,Yb), point(C,Xc,Yc),

(Ya-Yb)*(Xa-Xc)=(Ya-Yc)*(Xa-Xb),

Xa <Xb, Xb <Xc.

aligned(A,B,C):- point(A,Xa,Ya),

point(B,Xb,Yb), point(C,Xc,Yc),

(Ya-Yb)*(Xa-Xc)=(Ya-Yc)*(Xa-Xb),

Xa >Xb, Xb >Xc.

exists_point_smaller_x(A):- point(A,Xa,_), point(B,Xb,_), Xb <Xa.

exists_point_smaller_x(A):- point(A,Xa,Ya), point(B,Xb,Yb), Xb < = Xa, Yb <Ya.

Listing 7: Jarvis March in ASP

For this reason, we also experimented with a modification of the program in listing 7, which is stratified; the only difference is that in clauses 4 and 5 the atom hull(A) is removed. The complexity becomes O (n³) but gringo is now able to completely solve it and compute the convex hull.

3.3 Exploiting acyclicity checker

The clingo solver includes an efficient acyclicity checker [8]: edges identified with the #edge keyword must form an acyclic graph. Such acyclicity can be exploited in TSP solving as follows.

1 { cycle(A,B): cost(A,B,_) } 1:- point(A,_,_).

1 { cycle(A,B): cost(A,B,_) } 1:- point(B,_,_).

#edge(A,B): cycle(A,B), A != 1, B != 1.

#minimize{ C,A,B: cycle(A,B),cost(A,B,C) }.

Listing 8: Euclidean TSP encoding in ASP with acyclicity constraints

The rule in line 3 creates an atom #edge(A, B) for each edge of the TSP except for those insisting on one node (namely, node 1). Since clingo automatically imposes that the #edge predicate forms an acyclic graph, this line imposes that the set of edges in the TSP (excluding one node) form an acyclic graph, or, stated otherwise, that there are no sub-tours in the TSP.

4 Experimental evaluation

To assess the effectiveness of the proposed algorithms, we compared experimentally a classical ASP encoding for the Euclidean TSP with encodings using the reasoning based on geometric properties presented in the previous section. We also devised experiments to compare the speedup of the additional filtering based on geometric information on an ASP solver and a CLP(FD) solver. As CLP(FD) solver, we chose the ic library of ECLⁱPS^e v.7.0 build #54 [33] and as ASP solver we chose clingo 5.6.2 [21].

To generate realistic Euclidean TSP instances, we used the generator of the DIMACS challenge [12], that provides instances in two classes: uniform and clustered. We randomly generated instances from 10 to 50 nodes, in both classes. For each size and class, we generated 16 instances.

Uniform random generated instances consist of integer coordinate points uniformly distributed in a square of 10³ side. Randomly generated instances consist of clusters of points, whose centers are uniformly distributed in a square of 10³ side. Each point is then randomly associated with a cluster center and two normally distributed variables, each of which is then multiplied by 10³/# nodes, rounded, and added to the corresponding integer coordinate of the chosen center to obtain the coordinates of the point.

4.1 Comparing convex hull algorithms

In Section 3.2, we introduced three approaches for computing the convex hull in ASP. In Fig. 3, we compare the average grounding time of the three convex hull algorithms as the instance size varies. As depicted in the figure, the most efficient algorithm is labeled as JARVIS2 which is the stratified version of the JARVIS algorithm (see Listing 7). Computing the hull using the triangle method (see Listing 6), called INSIDE, achieves intermediate performance. It is noteworthy that for all methods, the grounding time remains below one second, even for the largest instances tested. The size of the ground program is almost the same for the two stratified methods (JARVIS2 and INSIDE), but it is larger for the JARVIS algorithm; for example, in the TSPLIB [32] instance att48, the size of the ground program is 41.147MiB using INSIDE, 50.833MiB for JARVIS and 41.148MiB for JARVIS2.

A second comparison made on solving times of randomly generated Euclidean TSP instances shows how the three algorithms are actually comparable; refer to the scatter plot in Fig. 2 for more details. In light of the results just presented, for the remainder of the experimental evaluation, we will consider the JARVIS2 model as the preferred hull computation method.

Fig. 2

Comparison of solving time using different methods for convex hull calculation.

Fig. 3

Comparison of grounding time of different methods for convex hull calculation.

4.2 Results on random generated Euclidean TSP instances

For simplicity, we summarized in Table 1 the structure of the various ASP programs tested in the experimental evaluation. The ASP_BASE program is the reference encoding for TSP in ASP (Listing 1), ASP_HULL and ASP_NOCROSS implement only convex hull-based pruning (Section 3.2) and crossing removal-based pruning (Section 3.1), respectively. Finally, ASP_GEOMETRIC implements both geometric reasoning. We also tested Euclidean TSP encodings based on acyclicity constraint called ASP_ACY_BASE and ASP_ACY_GEOMETRIC.

Table 1
ASP programs tested in the experimental evaluation. In all experiments we used the stratified variant of Listing 7

Name Configuration

ASP_BASE Listing 1 + symmetry breaking*

ASP_NOCROSS Listing 1 + 2

ASP_HULL Listing 1 + 3 + 4 + 5 + 7*

ASP_GEOMETRIC ASP_NOCROSS + ASP_HULL

ASP_ACYBASE Listing 8 + symmetry breaking

ASP_ACYGEOMETRIC Listing 8 + 2 + 3 + 4 + 5 + 7*

Name	Configuration
ASP_BASE	Listing 1 + symmetry breaking
ASP_NOCROSS	Listing 1 + 2
ASP_HULL	Listing 1 + 3 + 4 + 5 + 7*
ASP_GEOMETRIC	ASP_NOCROSS + ASP_HULL
ASP_ACYBASE	Listing 8 + symmetry breaking
ASP_ACYGEOMETRIC	Listing 8 + 2 + 3 + 4 + 5 + 7*

Since convex hull-based constraints also behave as symmetry breaking by imposing the direction of cycle, we added the following constraint

:- cycle(A,1), cycle(1,B), A>B.

to the ASP_BASE program to remove symmetries and have a fair comparison.

All tests were run with a time limit of 1800s on Intel^® Xeon^® Platinum 8358 running at 2.6GHz, using only one core and with 8GB of reserved memory.

Figure 4 compares ASP encodings without acyclicity constraints. In Fig. 4a, each point represents one instance solved with ASP_BASE encoding (whose running time is in the x axis) and with one of the encodings exploiting geometric reasoning (whose running time is reported in the y axis). Since all points stand below the y = x straight line, the geometric properties were able to speedup the running time in all instances. Note the log-log scale: speedups from 1 to 3 orders of magnitude were often obtained. Taken individually, both ASP_NOCROSS and ASP_HULL encodings showed a lower solving time when compared with ASP_BASE encoding. However, the best encoding turns out to be ASP_GEOMETRIC containing all the geometric reasoning presented in this paper.

Fig. 4

Comparison of Euclidean TSP ASP encodings using different geometric reasoning.

The graph in Fig. 4b shows how the average solving time varies with respect to the size of Euclidean TSP instances. Each point is calculated as the geometric mean of the solving times for the 32 instances of that size present in the dataset (16 uniform and 16 clustered). In particular, note how the use of the ASP_GEOMETRIC encoding allows solving instances containing up to 30% more nodes.

Figure 5 shows the effect of the acyclicity checker; the use of the acyclicity checker provided very small improvement, and mostly in conjunction with the use of geometric constraints.

Fig. 5

Comparison of ASP encoding using acyclicity constraint.

The comparison of grounding times for different encodings is depicted in Fig. 6. Specifically, the grounding time for the ASP_GEOMETRIC encoding increases quadratically as the number of cross/4 predicates grows quadratically. However, the time required for grounding never exceeds 10 seconds in the larger instances thus remaining negligible compared to the solving time.

Fig. 6

Comparison of grounding time of different Euclidean TSP ASP encodings.

We also compared the grounding size of the three constraints proposed in Section 3.2 to exploit the information about the convex hull. Figure 7 illustrates how the grounding size changes with varying problem sizes. For each number of nodes, we calculated the average grounding size obtained. The lines on the graph correspond to the visit-in-order constraint (Listing 3), the incoming angle constraint (Listing 4), and the path constraint (Listing 5), as well as the combination of all three constraints used together (see ASP_HULL). To compute the convex hull we used the stratified variant of Listing 7. As illustrated in the graph, the visitin- order constraint constraint is the cheapest, as it excludes only O (h²) edges (where h is the number of nodes on the border of the hull). Constraint incoming angle constraint excludes O (hn²) edges, while path constraint is clearly the most expensive, in terms of grounding size. Overall, the grounding size does not exceed 2.5 MiB with all constraints applied simultaneously.

Fig. 7

Comparison of average grounding size of the ASP encoding for the three constraints proposed to exploit the information about the convex hull.

4.3 Results on TSPLIB instances

We compared ASP_BASE and ASP_GEOMETRIC approaches also on instances taken from the TSPLIB [32], the Concorde website 1 and the CITIES dataset 2 up to 100 nodes. These sources provide various types of instances, including Euclidean instances, depicted as point sets in a two dimensional plane, and geographical instances, depicted as point sets (with latitude and longitude coordinates) on Earth’s surface. We included all Euclidean instances and additionally incorporated geographical instances where the cities to be visited are situated within a confined region of the Earth’s surface, allowing for the approximation of geographical distances using Euclidean distances.

All tests were run with a time limit of 86400s on Intel^® Xeon^® Platinum 8358 running at 2.6GHz, using only one core and with 8GB of reserved memory.

Figure 8 shows, for each tested instance, the quality of the best solution found within the timeout (the lower the better). For each method, the y-axis reports the ratio of objective function of the best found solution with respect to the optimal solution (obtained from the literature), so a value of 1 means that the ASP program was able to obtain the real optimal solution. The program using geometric filtering consistently achieves a superior solution quality, with only one exception. On average, the length of the tour obtained with the geometric reasoning was 75% that obtained with the basic version. This is a notable improvement for real life applications: for a company, saving 25% of the travel time and 25% of the fuel is a significant improvement (not to count the saving in polluting emissions). To better show how the solution quality improves varying the allotted running time, Figure 9 illustrates the average anytime behavior of the two approaches across all 23 tested instances. It is noteworthy that the ASP_GEOMETRIC model ensures better solutions compared to ASP_BASE within the timeout. Moreover, the quality of solutions improves significantly even during the initial stages of the solving procedure.

Fig. 8

Comparison of Euclidean TSP ASP encodings on TSPLIB instances. The chart illustrates the ratio between the best solution found within the timeout and the known optimal solution (lower is better).

Fig. 9

Performance profile of Euclidean TSP ASP encodings on TSPLIB instances. A higher curve suggests less efficient performance, while a lower curve reflects quicker convergence towards the optimal solution.

4.4 Comparing ASP and CLP(FD)

In CLP(FD), we adopted the successor representation. In the successor representation, n variables Next_i are defined Next = (Next₁, Next₂, …, 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 basic Euclidean TSP encoding includes an alldifferent(Next) constraint on the Next array of variables, which ensures that each node has exactly one incoming edge (degree constraints), and a circuit(Next) [2 , 27] constraint (sometimes called nocycle) that avoids subtours, i.e., cycles of length less than n. This encoding also includes, as redundant representation, a set of Prev variables: Prev_i represents the node that precedes node i in the circuit.

The CLP_GEOMETRIC program includes the basic CLP(FD) encoding for Euclidean TSP together with constraints for crossing removal and convex hull-based pruning, which in this program also acts as a symmetry breaking constraint. CLP_BASE program, on the other hand, besides the basic encoding, includes only the Next₁ < Prev₁ symmetry breaking constraint.

A comparison of how filtering based on geometric information performs differently on an ASP solver and a CLP(FD) solver is presented in Figure 10. We decided to compare the speedup obtained through the exploitation of geometric reasoning. In Figure 10a, where each point represents an instance of Euclidean TSP, the speedup obtained in the CLP(FD) solver (x-axis) is compared with the speedup obtained in the ASP solver (y-axis). On average, the use of geometric reasoning shows, on a single instance, a greater speedup in performance on the ASP solver than on the CLP(FD) solver.

Fig. 10

Comparison of geometric filtering on an ASP solver (clingo) and a CLP(FD) solver (ECLⁱPS^e).

Considering solving times, those recorded on the CLP(FD) solver still remain lower than those shown by the ASP solver. The cactus plot in Figure 10b shows in the x-axis the number of instances that could be solved (with proof of optimality) within the runtime plotted in the y-axis. It shows that CLP_BASE without geometric constraints even succeeds in solving 191 more instances than the encoding ASP_GEOMETRIC that uses geometric constraints instead. If we further compare CLP_GEOMETRIC using geometric reasoning with ASP_GEOMETRIC the gap becomes even more marked as the higher number of instances solved within the time limit increases from 191 to 321.

One difference between the CLP(FD) and the ASP encodings stands in the propagation they perform. In CLP(FD), we implemented a binary nocrossing constraint; each instance of the constraint involves two points and ensures that the edges exiting the two nodes do not cross each other. The propagators perform an arc-consistency propagation [5], that in the context of the nocrossing constraint means that if an outgoing edge AB from a node A crosses all possible edges exiting from node C, then the value B is removed from the domain of A (or, the edge AB cannot be taken).

In the ASP encoding, the crossing avoidance is delegated to the integrity constraint in line 1 of listing 2. Such a constraint ensures that if an edge AB is taken, then all edges CD that cross AB cannot be taken as edges of the TSP; such constraint is only activated when one of the involved edge atoms is true, i.e., when the solver tries to fix an edge as belonging to the TSP. Instead, unit propagation does not impose that if all edges exiting from a node C cross the edge AB, then edge AB cannot be taken. Of course a modern solver might perform more clever reasoning, involving both the integrity constraint 1 of listing 2 and clause 1 in listing 1, and conflict-directed clause learning might help in subsequent search steps, but unit propagation alone cannot detect that the edge AB should be removed as soon as it crosses all edges exiting from C.

5 Discussion

Considering the comparison between the ASP and the CLP(FD) approaches, two aspects must be noticed. Indeed, the CLP(FD) approach is faster, and able to solve more instances. On the other hand, the implementation of the propagators took several development time, was based on some non-trivial theorems that were key to reduce the computational complexity of the propagators [5], and takes about 1500 lines of code. The ASP approach, instead, is completely listed in this article, and is based on declarative considerations, without the need to develop complex propagation algorithms.

For example, in order to remove crossings in ASP, one only needs to define what a crossing is and add an integrity constraint, as seen in Listing 1. A similar approach could also be employed in CLP(FD), and indeed it would be a rapid prototyping approach to quickly implement the idea; on the other hand, in order to make it efficient a significant amount of development time would be required.

So while CLP may offer superior speed, it also entails greater complexity when compared to ASP. While we do not claim that our experience with this application is universally representative, we suggest that the improved readability of ASP encodings could be a reason that sheds light on a pertinent aspect of the evolving landscape of industrial applications, wherein ASP is steadily gaining traction owing to its balance of efficiency and simplicity.

In conclusions, we believe that the geometric reasoning on ASP provided a significant speedup and can help extending the applicability of ASP.

6 Related work

The encoding of the TSP in ASP is a classic one and can be fond, e.g., in [15], which presents also a variant not following the “guess and check” methodology; it does not provide computational results to compare the approaches.

The TSP was also addressed in extensions of ASP, such as CASP. Lierler [28] compares experimentally several CASP systems on the travelling salesman problem.

CASP was used to solve a Delivery Problem, in which robots have to move items inside a warehouse [30]. The extensive work describes in detail the problem, that can be considered as a problem of multi-agent path finding, and it includes a routing component within the warehouse, which is related to the TSP and in which possible routes are identified within a graph, and it also includes a scheduling component, as the exact timing in which each robot is in each location influences the synchronization of activities. In particular, the robots should not collide, and this is achieved by declaring conflict zones, that can be occupied by at most one robot at the same time. The graphs are very large, so the authors decide to reduce the search space in various ways, possibly excluding the optimal solution but allowing them to obtain good solutions quickly. The graphs they consider are planar, meaning that there are no crossings between edges. They successfully employ clingo[dl] [25].

In future work, we plan to implement our geometric reasoning also in CASP and study the speedup provided also in CASP implementations.

Several works address the TSP or some of its variants in CP. One of the most successful formulation is the so-called successor representation, in which a variable is associated to each node of the graph, and its value represents the following node in the TSP. Such formulation requires the alldifferent constraint [31] and also a constraint named circuit or nocycle, which imposes to avoid subtours (i.e., shorter cycles that do not involve all nodes). Several propagation algorithms have been proposed for the circuit constraint [2 , 27]. Other works address extensions of the TSP, such as the TSP with time windows. Pesant et al. [29] address the TSPTW and exploit the circuit constraint together with the minimum spanning tree relaxation. Focacci et al. [17, 18] propose reduced costs filtering to optimization constraints, and in particular use the assignment problem and the minimum spanning forest relaxation. Benchimol et al. [4] propose a Lagrangian relaxation scheme (based on the Held and Karp scheme) to perform improved reduced cost filtering in the TSP. Isoart et al. [24] introduce further pruning based on k-cutsets. Deudon et al. [14] trained a neural network to predict the most promising variable to branch on during search, based on the coordinates of the points of an Euclidean TSP.

7 Conclusion

In this paper we applied reasoning based on geometric constraints to the Euclidean TSP (that was developed in a previous publication [5] in the context of Constraint Logic Programming), a widely used case of the famous Traveling Salesperson Problem, within an ASP-based solving approach. The classical encoding of the TSP in ASP was significantly improved thanks to the geometric reasoning: without sacrificing the simplicity, declarativeness and succinctness of the the approach, speedups between 1 and 3 orders of magnitude were easily obtained. Experiments on real-life instances, such as those in the TSPLIB, have also demonstrated the effectiveness of the proposed approach even when obtaining the optimal solution is not required or feasible.

In future work, we plan to test additional ASP approaches for the Euclidean TSP, particularly by introducing difference logic and using the clingo[dl] CASP solver [25]. Another interesting research direction is the use of machine learning to select only those nocrossing constraints that are most effective; this approach was already studied in the CLP(FD) context [3] and provided a further speedup. Moreover, the use of Constraint ASP approaches could provide further improvements and also the geometric reasoning outlined in this paper could be adapted and applied to other routing problems on the Euclidean plane, such as the Euclidean Vehicle Routing Problem or the Euclidean Generalized TSP.

Finally, we plan to study how geometric constraints presented in this paper can be extended into non-Euclidean spaces such as TSP geographic instances where coordinates are expressed in terms of latitude and longitude.

Conflict of interest

The authors declare none.

Footnotes

Acknowledgements

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

Research funded by the Italian Ministry of University and Research through PNRR - M4C2 - Investimento 1.3 (Decreto Direttoriale MUR n. 341 del 15/03/2022), Partenariato Esteso PE00000013 - “FAIR - Future Artificial Intelligence Research” - Spoke 8 “Pervasive AI”, funded by the European Union under the NextGeneration EU programme.

References

Antonio

, Faster line segment intersection. In Kirk

, editor, Graphics Gems III, (1992), pp. 199–202, San Diego, CA, USA. Academic Press Professional, Inc. ISBN 0-12-409671-9.

Beldiceanu

, Contejean

, Introducing global constraints in CHIP, Math. Comput. Model.20(12) (1994), 97–123, December. ISSN 0895-7177.

Bellodi

, Bertagnon

, Gavanelli

, Zese

, Improving the efficiency of Euclidean TSP solving in constraint programming by predicting effective nocrossing constraints. In Baldoni

Matteo

, Bandini

Stefania

, editors, AIxIA 2020 - Advances in Artificial Intelligence -XIXth International Conference of the Italian Association forArtificial Intelligence, Virtual Event, November 25-27, 2020, Revised Selected Papers, volume 12414 of Lecture Notes in Computer Science, (2020), pp. 318–334. Springer. doi: 10.1007/978-3-030-77091-4_20.

Benchimol

, van Hoeve

W.J.

, Régin

J.-C.

, Rousseau

L.-M.

, Rueher

, Improved filtering for weighted circuit constraints, Constraints17(3) (2012), 205–233.

Bertagnon

, Gavanelli

, 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, pp. 1412–1419. AAAI Press, 2020. doi: 10.1609/aaai.v34i02.5498.

Bertagnon

, Gavanelli

, Geometric reasoning on the traveling salesperson problem: comparing answer set programming and constraint logic programming approaches. In Arias

, Batsakis

, Faber

, Gupta

, Pacenza

, Papadakis

, Robaldo

, Rückschloß

, Salazar

, Saribatur

Z.G.

, Tachmazidis

, Weitkämper

, Wyner

A.Z.

, editors, Proceedings of the International Conference on Logic Programming 2023 Workshops co-located with the 39th International Conference on Logic Programming (ICLP 2023), London, United Kingdom, July 9th and 10th, 2023, volume 3437 of CEUR Workshop Proceedings. CEUR-WS.org, 2023.

Bertagnon

, Gavanelli

, Zanotti

, ASPECT: Answer Set rePresentation as vEctor graphiCs in laTex. In Dovier

Agostino

, Formisano

Andrea

editors, Proceedings of the 38th Italian Conference on Computational Logic, Udine, Italy, June 21-23, 2023, volume 3428 of CEUR Workshop Proceedings. CEUR-WS.org, 2023. URL https://ceur-ws.org/Vol-3428/paper3.pdf.

Bomanson

, Gebser

, Janhunen

, Kaufmann

, Schaub.

, Answer set programming modulo acyclicity, Fundamenta Informaticae147 (2016), 63–91.

Calimeri

, Faber

, Gebser

, Ianni

, Kaminski

, Krennwallner

, Leone

, Maratea

, Ricca

, Schaub

, ASP-Core-2 input language format, Theory Pract. Log. Program.20(2) 2020. doi:10.1017/S1471068419000450

10.

Calimeri

, Manna

, Mastria

, Morelli

M.C.

, Perri

, Zangari

, I-DLVsr: A stream reasoning system based on I-DLV, Theory Pract. Log. Program.21(5) (2021), 610–628. doi: 10.1017/S147106842100034X.

11.

Caseau

, Laburthe

, Solving small TSPs with constraints. In Naish

Lee

, editor, Logic Programming, Proceedings of the Fourteenth International Conference on Logic Programming, Leuven, Belgium, July 8-11, 1997, pp. 316–330. MIT Press, 1997. ISBN 0-262-64035-X.

12.

Cirasella

, Johnson

D.S.

, McGeoch

L.A.

, Zhang

, The asymmetric traveling salesman problem: Algorithms, instance generators, and tests. In Buchsbaum

A.L.

, Snoeyink

, editors, Algorithm Engineering and Experimentation, Third International Workshop, ALENEX 2001, Washington, DC, USA, January 5-6, 2001, Revised Papers, volume 2153 of Lecture Notes in Computer Science, pp. 32–59. Springer, 2001. doi: 10.1007/3-540-44808-X_3

13.

Deineko

V.G.

, van Dal

, Rote

, , The convex-hull-and-line traveling salesman problem: A solvable case, Inf. Process. Lett.51(3) (1994), 141–148.

14.

Deudon

, Cournut

, Lacoste

, Adulyasak

, Rousseau

L.-M.

, Learning heuristics for the TSP by policy gradient. In Jan van Hoeve

, editor, Integration of Constraint Programming, Artificial Intelligence, and Operations Research - 15th International Conference, CPAIOR 2018, Delft, The Netherlands, June 26-29, 2018, Proceedings, volume 10848 of Lecture Notes in Computer Science, pp. 170–181. Springer, 2018. ISBN 978-3-319-93030-5.

15.

Faber

, Ricca

, Leone

, Answer set programming. In Wah

, editor, Wiley Encyclopedia of Computer Science and Engineering. Wiley, Hoboken, 2008.

16.

Flood

M.M.

, The traveling-salesman problem, Operations Research4, 1956.

17.

Focacci

, Lodi

, Milano

, Embedding relaxations in global constraints for solving TSP and TSPTW, Ann. Math. Artif. Intell.34(4) (2002), 291–311.

18.

Focacci

, Lodi

, Milano

, A hybrid exact algorithm for the TSPTW, INFORMS Journal on Computing14(4) (2002), 403–417.

19.

Francis

K.G.

, Stuckey

P.J.

, Explaining circuit propagation, Constraints19(1) (2014), 1–29. doi: 10.1007/s10601-013-9148-0.

20.

Gebser

, Harrison

, Kaminski

, Lifschitz

, Schaub

, Abstract gringo, Theory Pract. Log. Program.15(4-5) (2015), 449–463. doi: 10.1017/S1471068415000150

21.

Gebser

, Kaminski

, Kaufmann

, Schaub

, Multi-shot ASP solving with clingo, Theory Pract. Log. Program.19(1) (2019), 27–82. doi: 10.1017/S1471068418000054

22.

Gelfond

, Lifschitz

, The stable model semantics for logic programming. In Kowalski

Robert A.

, Bowen

Kenneth A.

, editors, ICLP, (1988), pp. 1070–1080. MIT Press.

23.

Giunchiglia

, Lierler

, Maratea

, Answer set programming based on propositional satisfiability, J. Autom. Reason.36(4) (2006), 345–377. doi: 10.1007/s10817-006-9033-2

24.

Isoart

, 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.

25.

Janhunen

, Kaminski

, Ostrowski

, Schellhorn

, Wanko

, Schaub

, Clingo goes linear constraints over reals and integers, Theory Pract. Log. Program.17(5-6) (2017), 872–888. doi: 10.1017/S1471068417000242

26.

Jarvis

R.A.

, On the identification of the convex hull of a finite set of points in the plane, Information Processing Letters2(1) (1973), 18–21. ISSN 0020-0190. doi: https://doi.org/10.1016/0020-0190(73)90020-3.

27.

Genç Kaya

, Hooker

J.N.

, A filter for the circuit constraint. In Frédéric Benhamou, editor, Principles and Practice of Constraint Programming - CP 2006, 12th International Conference, CP 2006, Nantes, France, September 25-29, 2006, Proceedings, volume 4204 of Lecture Notes in Computer Science, pp. 706–710. Springer, 2006. ISBN 3-540-46267-8.

28.

Lierler

, Constraint answer set programming: Integrational and translational (or smt-based) approaches, Theory and Practice of Logic Programming23(1) (2023), 195225. doi: 10.1017/S1471068421000478.

29.

Pesant

, Gendreau

, Potvin

J.-Y.

, Rousseau

J.-M.

, An exact constraint logic programming algorithm for the traveling salesman problem with time windows, Transportation Science32(1) (1998), 12–29.

30.

Rajaratnam

, Schaub

, Wanko

, Chen

, Liu

, Cao Son

, Solving an industrial-scale warehouse delivery problem with answer set programming modulo difference constraints, Algorithms16(4) (2023). ISSN 1999-4893. doi: 10.3390/a16040216.

31.

Régin

J.-C.

, A filtering algorithm for constraints of difference in CSPs. In Hayes-Roth

Barbara

, Korf

Richard E.

, editors, Proceedings of the 12th National Conference on Artificial Intelligence, Seattle, WA, USA, July 31 - August 4, 1994, Volume 1., pp. 362–367. AAAI Press / The MIT Press, 1994. ISBN 0-262-61102-3.

32.

Reinelt

, TSPLIB - A traveling salesman problem library, INFORMS Journal on Computing3(4) (1991), 376–384.

33.

Schimpf

, Shen

, Eclipse - from LP to CLP, Theory and Practice of Logic Programming12(1-2) (2012), 127–156.

34.

Syrjänen

, Niemelä

, The smodels system. In Eiter

, Faber

, Truszczynski

, editors, Logic Programming and Nonmonotonic Reasoning, 6th International Conference, LPNMR 2001, Vienna, Austria, September 17-19, 2001, Proceedings, volume 2173 of Lecture Notes in Computer Science, (2001), pp. 434–438. Springer. doi: 10.1007/3-540-45402-0_38