A memetic algorithm with optimal recombination for the asymmetric travelling salesman problem

Abstract

We propose a new memetic algorithm with optimal recombination for the asymmetric travelling salesman problem (ATSP). The optimal recombination problem (ORP) is solved in a crossover operator based on a new exact algorithm that solves the ATSP on cubic digraphs. A new mutation operator makes random jumps in 3-opt or 4-opt neighborhoods. The initial population is constructed by means of greedy constructive heuristics. The 3-opt local search is used to improve the initial and the final populations. A computational experiment on the TSPLIB instances shows that the proposed algorithm yields results competitive to those of other well-known algorithms for ATSP and confirms that the ORP may be used successfully in memetic algorithms.

Appendix

The appendix contains an outline of the algorithm of Eppstein from [5], a proof of Lemma 1 (which uses some ideas from the algorithm of D. Eppstein), and the details about the restarting rule from [6].

1.1 Algorithm of D. Eppstein

For the sake of completeness, here we present the recursive algorithm of D. Eppstein proposed in [5] for listing all Hamiltonian cycles in a graph of a degree three.

Let $G=(V,E)$ be a simple graph with maximum degree three, F be the set of forced edges and ${U:=E\backslash F}$. The algorithm of D. Eppstein enumerates all Hamiltonian cycles of G, containing all forced edges, as shown in Algorithm 1. At each stage we firstly reduce the size of the input graph without branching (Step 1). Secondly, we choose an edge to branch on, and consider two subproblems (Step 2), the one in which the edge is forced to be in the solution and the one in which it is excluded. These two subproblems are solved recursively. Steps of the algorithm either return or reduce the input graph to one or two smaller cubic graphs.

Algorithm 2. Listing all Hamiltonian Cycles Given the Forced Edges

Input: Simple cubic graph G with a set of forced edges F.

Output: All Hamiltonian cycles of the graph G, containing all forced edges.

• 1. Repeat the following steps until one of the steps returns or none of them applies:

  • (a) If G contains a vertex with degree zero or one, or if F contains three edges meeting at a vertex, backtrack.

  • (b) If G contains only three vertices and F consists of a Hamiltonian cycle, output the cycle and backtrack.

  • (c) If G contains a vertex with degree two, add its incident edges to F.

  • (d) If G contains a triangle xyz, and the non-triangle edge incident to x belongs to F, add edge yz to F.

  • (e) If F contains exactly two edges $\{x,y\}$ and $\{y,z\}$ meeting at some vertex y, remove from G that vertex and any other edge incident to it, and replace the two edges by single edge $\{x,z\}$. If this contraction would lead to two parallel edges in G, remove the other edge from G.

• 2. If F is nonempty, let $\{x,y\}$ be any edge in F and $\{y,z\}$ be an adjacent edge in $G\backslash F$. Otherwise, if F is empty, let $\{y,z\}$ be any edge in G. Call the algorithm recursively on the two graphs $G,\ F \cup \{y,z\}$ and $G \backslash \{y,z\},\ F$.

Theorem 2 is a combination of results proved in [5].

Theorem 2

[5] The number of Hamiltonian cycles in $G=(V,F\cup U)$ is $O(2^{|U|/4})$ in the worst case and these cycles can be listed be means of Algorithm 2 in $O(2^{|U|/4})$ time and linear space. Each step of Algorithm 2 can be implemented in constant time per step.

1.2 Proof of Lemma 1 from section 2.1

Consider the forced ATSP on a directed cubic graph ${\mathcal G}=({{\mathcal {V}}},{{\mathcal {A}}})$, where ${{\mathcal {V}}}$ is the set of vertices and ${{\mathcal {A}}}$ is the set of arcs. Each vertex $v\in {{\mathcal {V}}}$ has a degree at most three. Every arc $(v_i,v_j)\in {{\mathcal {A}}}$ is associated with a length (weight) $c(v_i,v_j)\ge 0$. Let ${{\mathcal {F}}}$ be the set of forced arcs in the graph ${{\mathcal {G}}}$. We propose an exact algorithm for finding a minimum length Hamiltonian circuit in ${{\mathcal {G}}}$.

Initially the graph ${{\mathcal {G}}}$ is modified by means of the following preprocessing procedure. For each arc $(u,v)\in {{\mathcal {G}}}$, we check whether there is an arc (v, u) in the graph ${{\mathcal {G}}}$. If such an arc (v, u) is not present in ${{\mathcal {A}}}$, then we add it to the graph and set ${c(v,u):=+\infty }$. If such an arc (v, u) is in ${{\mathcal {A}}}$, then we identify the following four cases:

• 1) if $(u,v)\in {{\mathcal {F}}}$ and $(v,u)\in {{\mathcal {F}}}$, then we report that the forced ATSP on digraph ${{\mathcal {G}}}$ is infeasible;

• 2) if $(u,v)\in {{\mathcal {F}}}$, but $(v,u)\notin {{\mathcal {F}}}$, then ${c(v,u):=+\infty }$;

• 3) if $(v,u)\in {{\mathcal {F}}}$, but $(u,v)\notin {{\mathcal {F}}}$, then ${c(u,v):=+\infty }$;

• 4) if $(u,v)\notin {{\mathcal {F}}}$ and $(v,u)\notin {{\mathcal {F}}}$, then do nothing.

If an arc (u, v) is forced in ${{\mathcal {G}}}$ after the preprocessing, then the length of (v, u) is infinity. A solution to the modified problem may be obtained through enumeration of all feasible solutions to the TSP with forced edges on a supplementary graph ${\bar{G}} =({{\mathcal {V}}}, {\bar{E}})$, where a pair of vertices u, v is connected iff these vertices were connected by a pair of arcs in the digraph ${{\mathcal {G}}}$. An edge $\{u, v\}\in {\bar{E}}$ is assumed to be forced if $(u, v)\in {{\mathcal {F}}}$ or $(v, u)\in {{\mathcal {F}}}$. A set of forced edges in ${\bar{G}}$ will be denoted by ${\bar{F}}$. Given a set of forced arcs ${{\mathcal {F}}},$ a minimum length Hamiltonian circuit in ${{\mathcal {G}}}$ may be found using Algorithm 3, which is an adaptation of the algorithm proposed by Eppstein [5] (see Algorithm 2).

Algorithm 3. Finding a Minimum Length Hamiltonian Circuit Given the Forced Arcs

Input: Directed cubic graph ${{\mathcal {G}}}$ and the corresponding undirected graph ${\bar{G}}$ with a set of forced edges ${\bar{F}}$.

Output: A minimum length Hamiltonian circuit of the graph ${{\mathcal {G}}}$, containing all forced arcs, if such circuit exists.

• 1. Let $c_{\mathrm {bst}}$ denote the best found value of objective function. Initially set ${c_{bst}:=+\infty }$.

• 2. Repeat the following steps until one of the steps returns or none of them applies:

  • (a) If ${\bar{G}}$ contains a vertex with degree zero or one, or if ${\bar{F}}$ contains three edges meeting at a vertex, backtrack.

  • (b) If graph ${\bar{G}}$ contains only three vertices and ${\bar{F}}$ consists of a Hamiltonian cycle, compute thelengths of the circulations within this cycle in both directions using arc lengths from graph ${{\mathcal {G}}}$. If necessary, update the value of $c_{\mathrm {bst}}$ and the circulation corresponding to this value, and backtrack. Deleted vertices of the original input graph ${{\mathcal {G}}}$ will be added to the best found circuit at backtracking stages.

  • (c) If ${\bar{G}}$ contains a vertex with degree two, add its incident edges to ${\bar{F}}$.

  • (d) If ${\bar{G}}$ contains a triangle xyz, and the non-triangle edge incident to x belongs to ${\bar{F}}$, add edge $\{y,z\}$ to ${\bar{F}}$.

  • (e) If ${\bar{F}}$ contains exactly two edges $\{x,y\},$$\{y,z\}$ meeting at some vertex y, add edge $\{x,z\}$ to ${\bar{G}}$, and add arcs (x, z),  (z, x) to ${{\mathcal {G}}}$. If this leads to two parallel edges in ${\bar{G}}$ (or parallel arcs in ${{\mathcal {G}}}$), remove the other edge from ${\bar{G}}$ (or the other arc from ${{\mathcal {G}}}$). The lengths of arcs (x, z) and (z, x) are computed as follows: $c(x,z):=c(x,y)+c(y,z),$$c(z,x):=c(z,y)+c(y,x)$. Remove vertex y from ${\bar{G}}$ and ${{\mathcal {G}}}$.

• 3. If ${\bar{F}}$ is nonempty, let $\{x,y\}$ be any edge in ${\bar{F}}$ and let $\{y,z\}$ be an adjacent edge in ${\bar{G}}\backslash {\bar{F}}$. Otherwise, if ${\bar{F}}$ is empty, let $\{y,z\}$ be any edge in ${\bar{G}}$. Call the algorithm recursively on the following two input data sets: ${{\mathcal {G}}}$, ${\bar{G}},\ {\bar{F}} \cup \{y,z\}$ and ${{\mathcal {G}}}\backslash \{(y,z)\cup (z,y)\}$, ${\bar{G}}\backslash \{y,z\},\ {\bar{F}}$.

• 4. If ${c_{\mathrm {bst}}<+\infty }$, then return the Hamiltonian circuit corresponding to $c_{\mathrm {bst}}$. Otherwise no solution exists.

Algorithm 3 is recursive and at each call it either returns or reduces the input graph to two smaller cubic graphs by choosing arcs to branch on. In contrast to the algorithm from [5], our algorithm operates with a directed graph. It enumerates all Hamiltonian circuits of this graph and finds a minimum length one. Due to the structure of the initial graph ${{\mathcal {G}}}$, at most one circuit of finite length on three vertices is identified every time at Step 2b. The circuit corresponds to a circulation in the initial graph ${{\mathcal {G}}}$, containing all forced arcs.

Whenever a backtracking occurs in Algorithm 3, all vertices, edges and arcs deleted at the current stage are restored and these vertices are added to the best found circuit if it is needed. This leads to the same number of additional steps as has been performed at the current stage. So, each step of Algorithm 3 has a time complexity O(1), and the total number of steps is $O(2^{ (|{\bar{E}}|-|{\bar{F}}|)/4 })$ as for the algorithm from [5] (see Theorem 2). The time complexity of the preprocessing is just $O(|{\mathcal V}|)$, therefore, Lemma 1 holds.

1.3 Restarting rule

Schnabel census method was developed in biometrics for statistical estimation of the size of animal populations. According to this method, one takes repeated samples of size $n_0$ from a population and counts the number of distinct animals seen. Often it is assumed that the probability of catching any particular animal is the same. The sampled animals are marked, unless they were marked previously, and returned back into the population. Then a statistical estimate for the total number $\nu$ of individuals in population is computed on the basis of the total number of distinct animals marked in all the samples. In this paper, we apply the Schnabel census method to estimate the number of values that a discrete random variable may take with non-zero probability.

Let a parameter r define the length of the historical period considered for statistical analysis in the $\mathcal MA$ restart rule. Given a value of r, we assume that during the r latest iterations, all new offspring in the $\mathcal MA$ obeyed the same distribution and they may be treated as the sampled animals in the Schnabel census method. Then we apply the Schnabel census method in order to estimate the number $\nu$ of different solutions that may be visited with a positive probability, assuming that the current distribution of offspring remains unchanged.

We also assume that in the latest r iterations of a $\mathcal MA$, the observed sample consists of r independent offspring solutions. Let us define the random variable K as the number of distinct solutions in this sample. We make another simplifying assumption that all solutions, that may be generated with the current distribution, have equal probabilities. Then, for any fixed $\nu ,$ the random variable K has the following distribution:

$$\begin{aligned} \Pr \{K=k\} = \frac{\nu !}{(\nu -k)!}\frac{S(r,k)}{\nu ^r} , \end{aligned}$$

where $S(r,k)=\frac{1}{k!} \sum _{s=0}^k (-1)^k {k \atopwithdelims ()s} (k-s)^r$ is the Stirling number of the second kind. This distribution is also known as the Arfwedson distribution. The maximum likelihood estimate ${\hat{\nu }}^{\mathrm{ML}}$ for the unknown $\nu$ is

$$\begin{aligned} {\hat{\nu }}(r,k)=\text{ argmax } \left\{ \frac{\nu !}{(\nu -k)! \nu ^r}\right\} , \end{aligned}$$

(1)

where k is the number of different solutions actually generated on the latest r iterations. The value ${\hat{\nu }}^\mathrm{ML}={\hat{\nu }}(r,k)$ may be found from (1) by the standard one-dimensional optimization methods. The technical details as well as further literature references may be found e.g. in [23].

In this paper, we use the restart rule proposed in [6], which restarts an algorithm as soon as the estimate ${\hat{\nu }}^{\mathrm{ML}}$ becomes equal to k. The rationale behind this rule is that once the equality ${\hat{\nu }}^{\mathrm{ML}}=k$ is satisfied, most likely there are no more non-visited solutions in the area where the GA population spent the latest r iterations. In such a situation, it is more appropriate to restart the GA rather than to wait till the population distribution will significantly change by the evolutionary mechanisms.

The value of parameter r is chosen adaptively as follows: Whenever the best found solution is improved, we set r to be the population size N. If the best incumbent was not improved during the latest 2r iterations, then the value of r is doubled. Here we reset r to the population size when the best found solution is improved, assuming that if the best incumbent has been improved, this means that the population has reached a new unexplored area and the length of the historic period for analysis should be reduced. To reduce the CPU cost, the termination condition is checked only when the value of r is updated.