chapter_tsp.tex

\chapter{Shortest Hamiltonian cycle (Travelling salesman problem)}

The travelling salesman problem (TSP) asks the following question: "Given a list of cities and the distances between each pair of cities, what is the shortest possible route that visits each city exactly once and returns to the origin city?"

We will assume that there are roads (edges) between all cities (complete graph) and that the distances are Euclidean distances (Euclidean TSP).

We will write an approximation algorithm for TSP, which will be based on the minimum spanning tree (MST) algorithm.

\section{Approximation}

Implement the following 2-approximation algorithm (that means that the length of our solution will be better than 2 times the length of an optimal solution).
\begin{enumerate}
    \item Find minimal spanning tree of our graph (use built-in Sage function \verb`min_spanning_tree`).
    \item Run DFS on this tree.
    \item Take vertices in the order of increasing (DFS) start time.
\end{enumerate}

\medskip
\begin{sageCell}
def TSP_approximation(G):
    """
    Returns Hamiltonian cycle (travelling salesman circuit) using a 2-approximation algorithm.
    """
    mst = G.min_spanning_tree(by_weight=True)
    T = Graph(mst) # graph (tree) from edges

    # DFS with times
    r = T.vertices(sort=False)[0]
    _, start, _ = DFS_with_times(T, r)
    sort_start = sorted(list(start.items()), key=lambda p: p[1])
    cycle = []
    length = 0
    sort_start.append(sort_start[0])
    for i in range(1, len(sort_start)):
        u = sort_start[i - 1][0]
        v = sort_start[i][0]
        cycle.append((u, v))
        length += G.edge_label(u, v)
    return cycle, length
\end{sageCell}

\subsubsection*{Example}

Create the test example, a complete graph with 10 vertices with given vertex coordinates.

\medskip
\begin{sageCell}
def distance(a, b):
    """
    Return Euclidean distance between a = (ax, ay) and b = (bx, by)
    """
    ax, ay = a
    bx, by = b
    return math.sqrt((bx - ax)**2 + (by - ay)**2)

def set_euclidean_distances(G):
    """
    Set Euclidean distances as edge weights (labels)
    """
    pos = G.get_pos()
    for (u, v) in G.edges(sort=False, labels=False):
        G.set_edge_label(u, v, distance(pos[u], pos[v]))

H = graphs.CompleteGraph(10)
H.set_pos({0: [8, 1], 1: [0, 8], 5: [1, 0], 2: [5, 3], 3: [1.5, 7], 4: [2, 4],
  6: [6, 2], 7: [3, 1], 8: [2, 2], 9: [3, 3]})
set_euclidean_distances(H)
\end{sageCell}

\begin{sageCell}
    cycle, length = TSP_approximation(H)
\end{sageCell}
\begin{sageCell}
    cycle, length
\end{sageCell}
\begin{outCell}
   ([(0, 6),
     (6, 2),
     (2, 9),
     (9, 4),
     (4, 3),
     (3, 1),
     (1, 8),
     (8, 5),
     (5, 7),
     (7, 0)], 27.70534328046342)
\end{outCell}

Compare with the optimal solution, computed using the built-in Sage function\\ \verb`traveling_salesman_problem`.

\medskip
\begin{sageCell}
    def cycle_length(cycle, G):
        length = 0
        for u, v in cycle:
            length += G.edge_label(u, v)
        return length

    opt_cycle = H.traveling_salesman_problem(use_edge_labels=True)
    opt_length = cycle_length(opt_cycle.edges(sort=False, labels=False),
      opt_cycle)
    opt_length
\end{sageCell}
\begin{outCell}
    27.540829118717557
\end{outCell}

\section{Iterative improvement}

\subsection{2-changes on intersecting segments}

Improve a solution iteratively using 2-changes on \emph{intersecting} segments.

A \emph{2-change} is a transformation of a cycle by removing two non-consecutive edges and adding two other edges such that the resulting graph is still a cycle.

\medskip
\begin{sageCell}
def iterate_2_changes_intersecting(cycle, G, n=1000):
    """
    Iterate by eliminating intersections by 2-changes. Make at most n iterations
    """
    for k in range(n):
        inter = find_intersection(cycle, G)
        if inter == None:
            return cycle
        i, j = inter
        cycle = perform_2_change(cycle, i, j)
    return cycle
\end{sageCell}
(See the code for \verb`find_intersection` and \verb`perform_2_change` at the end of this chapter.)

\subsubsection*{Example}

\begin{sageCell}
    cycle = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 7), (7, 8), (8, 9),
      (9, 0)]
    new_cycle = iterate_2_changes_intersecting(cycle, H, 10)
    (cycle_length(new_cycle, H), cycle_length(cycle, H))
\end{sageCell}
\begin{outCell}
    (30.367268577721603, 46.94180782091779)
\end{outCell}

\subsection{2-changes on random edges}

Improve a solution iteratively using 2-changes on \emph{random}
non-adjacent cycle edges.

\medskip
\begin{sageCell}
def iterate_2_changes(cycle, G, n):
    min_length = cycle_length(cycle, G)
    min_cycle = cycle
    for k in range(n):
        i = randint(0, len(min_cycle) - 1)
        add = randint(2, len(min_cycle) - 2)
        j = (i + add) % len(min_cycle)
        new_cycle = perform_2_change(min_cycle, i, j)
        new_length = cycle_length(new_cycle, G)
        if new_length < min_length:
            min_length = new_length
            min_cycle = new_cycle
    return min_cycle
\end{sageCell}

\subsubsection*{Example}

\begin{sageCell}
    new_cycle = iterate_2_changes(cycle, H, 50000)
    (cycle_length(new_cycle, H), cycle_length(cycle, H))
\end{sageCell}
\begin{outCell}
    (27.669330960262805, 46.94180782091779)
\end{outCell}
\begin{sageCell}
    (cycle_length(new_cycle, H), opt_length)
\end{sageCell}
\begin{outCell}
    (27.669330960262805, 27.540829118717557)
\end{outCell}


\subsection{Code of auxiliary functions}

\begin{sageCell}
def perform_2_change(cycle, i, j):
    """
    Perform a 2-change on a (hamiltonian) cycle for edges with
    (non-consecutive) indices i and j. Cycle is a list of edges
    """
    if i > j:
        i, j = j, i
    e1 = cycle[i]
    e2 = cycle[j]
    v1, u1 = e1
    v2, u2 = e2
    result = []
    revert = False
    for k in range(i):
        result.append(cycle[k])
    result.append((v1, v2))
    for k in reversed(range(i + 1, j)):
        result.append(tuple(reversed(cycle[k])))
    result.append((u1, u2))
    for k in range(j + 1, len(cycle)):
        result.append(cycle[k])
    return result
\end{sageCell}

Intersection of two segments.
\begin{sageCell}
def find_intersection(cycle, G):
    """
    Find indices of two non-consecutive cycle edges which intersect and None if there are none
    """
    pos = G.get_pos()
    for i in range(len(cycle)):
        ei = cycle[i]
        l = len(cycle) if i > 0 else len(cycle) - 1
        for j in range(i + 2, l):
            ej = cycle[j]
            if do_intersect(pos[ei[0]], pos[ei[1]], pos[ej[0]], pos[ej[1]]):
                return (i, j)
    return None

def on_segment(p, q, r):
    if ((q[0] <= max(p[0], r[0])) and (q[0] >= min(p[0], r[0])) and
        (q[1] <= max(p[1], r[1])) and (q[1] >= min(p[1], r[1]))):
        return True
    return False

def orientation(p, q, r):
    val = (float(q[1] - p[1]) * (r[0] - q[0])) - (float(q[0] - p[0]) * (r[1] - q[1]))
    if (val > 0):
        return 1
    elif (val < 0):
        return 2
    else:
        return 0

# Check if two segments intersect
# https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/
def do_intersect(p1, q1, p2, q2):
    o1 = orientation(p1, q1, p2)
    o2 = orientation(p1, q1, q2)
    o3 = orientation(p2, q2, p1)
    o4 = orientation(p2, q2, q1)

    if ((o1 != o2) and (o3 != o4)):
        return True
    if ((o1 == 0) and on_segment(p1, p2, q1)):
        return True
    if ((o2 == 0) and on_segment(p1, q2, q1)):
        return True
    if ((o3 == 0) and on_segment(p2, p1, q2)):
        return True
    if ((o4 == 0) and on_segment(p2, q1, q2)):
        return True
    return False
\end{sageCell}