An improved Ant Colony System for the Sequential Ordering Problem
Rafał Skinderowicza
a University

of Silesia, Intitute of Computer Science,
Będzińska 39, 41-205 Sosnowiec, Poland

arXiv:1705.01076v1 [] 2 May 2017

©2017. This manuscript version is made available under the CC-BY-NC-ND 4.0 license
http://creativecommons.org/licenses/by-nc-nd/4.0/

Abstract
It is not rare that the performance of one metaheuristic algorithm can be improved by incorporating
ideas taken from another. In this article we present how Simulated Annealing (SA) can be used
to improve the efficiency of the Ant Colony System (ACS) and Enhanced ACS when solving the
Sequential Ordering Problem (SOP). Moreover, we show how the very same ideas can be applied
to improve the convergence of a dedicated local search, i.e. the SOP-3-exchange algorithm. A
statistical analysis of the proposed algorithms both in terms of finding suitable parameter values and
the quality of the generated solutions is presented based on a series of computational experiments
conducted on SOP instances from the well-known TSPLIB and SOPLIB2006 repositories. The
proposed ACS-SA and EACS-SA algorithms often generate solutions of better quality than the
ACS and EACS, respectively. Moreover, the EACS-SA algorithm combined with the proposed
SOP-3-exchange-SA local search was able to find 10 new best solutions for the SOP instances
from the SOPLIB2006 repository, thus improving the state-of-the-art results as known from the
literature. Overall, the best known or improved solutions were found in 41 out of 48 cases.
Keywords: Ant Colony System, Simulated Annealing, Sequential Ordering Problem,
combinatorial optimization

1. Introduction
In recent years, a large number of metaheuristic optimization algorithms (MOAs) has been proposed, and some of these were created based on inspiration drawn from natural phenomena [24].
Examples of these are the Ant Colony System algorithm that was inspired by the foraging behavior of certain species of ants and Simulated Annealing (SA) with some ideas taken from metallurgy [37, 11]. Metaheuristics are often applied to find solutions of an acceptable quality to difficult
combinatorial optimization problems, particularly NP-complete ones. A good example is the Sequential Ordering Problem (SOP) which consists in finding a minimum weight Hamiltonian path
on a directed graph with weights that is subject to precedence constraints among the nodes. Although less time-consuming than the exact approaches, MOAs differ in their efficiency, which can
sometimes be improved by combining ideas taken from other MOAs.

Email address: rafal.skinderowicz@us.edu.pl (Rafał Skinderowicz)

Accepted Manuscript. DOI: 10.1016/j.cor.2017.04.012

1.1. Contributions
The main aim of the paper is to show that Simulated Annealing could be used to improve
the convergence speed of ACS and Enhanced ACS (EACS) algorithms. The proposed solution is
easy to implement and does not increase the algorithms’ asymptotic complexity. Moreover, we
developed a modified version of the SOP-3-exchange local search (LS) heuristic as proposed by
Gambardella et al. [21] for the SOP. The modification, again, includes ideas taken from the SA to
allow the algorithm to escape local optima. A thorough experimental evaluation on a number of
SOP instances from well-known datasets confirms the efficiency of the proposed algorithms. In fact,
in several cases we obtained results of a better quality than those from state-of-the-art methods in
the literature [23, 26].
The paper is organized as follows: Section 2 focuses on recent ideas of improving the efficiency
of the ACS (and related algorithms), some of which include the SA. We also briefly discuss recent
work on solving the SOP. In Sec. 3 we describe our approach to improve the convergence speed of
the ACS by using the SA. Section 4 presents a similar approach, but to improve the local minima
escape ability of the SOP-3-exchange local search algorithm which is paired with the ACS and the
EACS when solving the SOP. Section 5 presents the results of computational experiments conducted
on two sets of SOP instances. The last section contains the conclusions and some ideas for further
work.
2. Related work
Multiple modifications to the ACO family of algorithms have been proposed in the literature.
Most of them refer to pheromone update rules and parameter tuning. Hassani et al. [15] proposed a
modified global pheromone update rule for the ACS in which not only the global best ant but also
all ants with inferior solutions may update the pheromone with probability calculated according to
the acceptance rule of the SA. The initial temperature was set arbitrarily to 100 and an exponential
cooling schedule was applied. The limited computational experiments showed that in most cases
the algorithm achieved better results than the Ant System. Bouhafs et al. [7] proposed a two-phase
approach based on the SA and ACS to solve the Capacitated Location-Routing Problem. The
SA was used to find facility locations while the ACS was used to solve the corresponding location
routing problem. In most cases The algorithm was able to improve the best-known solutions to the
problem.
In most of the ACO and SA combinations the latter plays the role of a local search used to
improve the quality of the solutions generated by the ants. Behnamian et al. [5] proposed a hybrid of
the ACO, SA and Variable Neighborhood Search algorithms for solving parallel-machine scheduling
problems. The SA was used to guide the dedicated LS. A successful combination of the ACS and
the SA was proposed by Ayob and Jaradat [4] for solving course timetabling problems. The SA
was used along with the Tabu Search to improve solutions generated by the ACS. The results
for the proposed algorithm were of better quality than those of the ACS alone or of the MAXMIN Ant System. The SA was again used as the LS for the ACS by Wassila and Boukra [43].
The approach was slightly faster but comparable in terms of solutions quality, than other natureinspired metaheuristics for the intrusion detection problem. Similarly, the SA played the role of an
LS improving the results generated by the ants in the ACS solving the Vehicle Routing Problem
with Time Windows [9]. Chen and Chien [10] proposed a complex hybrid of four metaheuristics,
namely of the Genetic Algorithm, SA, ACS, and the Particle Swarm Optimization, for solving the
TSP. The SA played the role of a mutation operator in the GA part of the hybrid. In a paper
2

by Xi et al. [44] the solutions generated by the ant system were later improved by the SA when
solving the 3D/2D fixed-outline floor planning problem. McKendall and Shang [32] used the SA
as the LS method in one of their Hybrid Ant System algorithms for solving the dynamic facility
layout problem. The resulting algorithm was able to improve some of the best known results for the
problem. Similarly, the solutions generated by the ACO were a starting point for the SA solving
the problem of managing energy resources considering intensive use of electric vehicles [42]. The
combined approach produced solutions of quality better than of the SA or ACO alone.
2.1. Sequential Ordering Problem
The Sequential Ordering Problem is a generalization of the Asymmetric TSP (ATSP). The goal
is to find the shortest Hamiltonian path from a starting city (source node) to a destination city
(final node) by going through each of the remaining cities (nodes) exactly once. Moreover, some
cities have to be visited before others. Due to precedence constraints, the problem is sometimes
referred to as the Precedence Constrained Traveling Salesman problem (PCATS). The SOP can be
viewed as a scheduling problem in which many jobs have to be scheduled on a single machine. The
processing times for the jobs are given along with the setup times between pairs of jobs. Also, some
jobs have to be completed before others. The goal is to minimize the total makespan [17]. Other
real-world problems that can be modeled as an instance of the SOP include the Single Vehicle
Routing Problem with pick-up and delivery constraints or the routing of a stacker crane in an
automatic storage system [2].
An instance of the SOP can be described using a graph G = (V, E) where V is a set of nodes
containing the starting us and the final vf nodes, and E = {(u, v)|u, v ∈ V, u 6= v} is a set of
weighted directed edges. Additionally, a precedence graph H = (V, R) is given, where R defines
the precedence constraints, i.e. an edge (u, v) ∈ R if node u has to precede node v in every feasible
solution. By definition, the starting node us precedes every other node, i.e. (us , v) ∈ R ∀v ∈
V \ {us }, and the final node uf has to be preceded by all other nodes, i.e. (u, vf ) ∈ R ∀u ∈ V \ {vf }.
The precedence graph, R, has to be acyclic for feasible solutions to exist.
On the one hand, the precedence constraints make solving the SOP more difficult than the
ATSP because the solution construction algorithms have to check for the precedence constraints.
On the other hand, the precedence constraints may limit the number of feasible solutions, which
can be beneficial for the exact methods, such as the branch-and-cut [29]. All of the SOP instances
from the SOPLIB2006 repository with the largest relative number of precedence constraints (60%)
were solved to optimality, whereas those in which the constraints concerned 15% percent of all edges
remained unsolved [26].
The SOP was introduced along with a mathematical model and exact algorithms by Escudero
et al. [16]. Several exact approaches for solving the SOP have since been proposed. Escudero et
al. applied Lagrangian Relaxation to solve the SOP, which the authors called the Relax-and-Cut
algorithm [18]. Hernádvölgyi et al. proposed a branch-and-bound algorithm with the lower bounds
obtained from homomorphic abstractions of the original states space [29]. The authors solved to
optimality several instances (with 40-50 vertices) from the TSPLIB repository.
Gouveia et al. proposed a cutting plane algorithm with the SOP formulations involving additional exponential-sized sets of inequalities [25]. The authors were able to improve the best known
lower bounds for many SOP instances from the TSPLIB repository and to solve to optimality instance p43.4 (the calculations took 15282 sec.). Later, Guveia and Ruthmair solved to optimality
several SOP instances from the SOPLIB2006 repository by using the branch-and-cut algorithms
combined with several preprocessing methods, heuristics, and separation routines [26]. They used a
3

single core of an Intel Xeon E5540 or E5649 both with a 2.53GHz clock. For most of the instances
the optima were found under an hour, but for 12 instances no optima were found with the time
limit set to 24 hours.
The exact methods are time consuming, particularly if the size of the problem reaches a few
hundred nodes, hence much of the research has been focused on heuristic algorithms for the SOP.
Guerriero and Mancini proposed a parallel roll-out heuristic in which several threads simultaneously visit different portions of the solution space and periodically exchange information about the
solutions found [27]. The algorithm was able to match the best-known solutions for most of the
SOP instances from the TSPLIB repository, although its main drawback was a high computational
cost.
Gambardella et al. proposed a combination of the Ant Colony System and a novel LS procedure
called the SOP-3-exchange [21]. The resulting algorithm, denoted as HAS-SOP, allowed to improve
many best-known results for many SOP instances from the TSPLIB repository. Monetamanni et
al. added to the HAS-SOP a Heuristic Manipulation Technique which creates and adds artificial
precedence constraints to the original problem [35]. The method led to better results, particularly
for large SOP instances.
A discrete Particle Swarm Optimization hybridized with the SOP-3-exchange heuristic was
proposed by Anghinolfi et al. [1]. The algorithm was able to improve many of the best results
presented in [21, 35]. Later, Gamabrdella et al., basing their findings on an analysis of the drawbacks
of the HAS-SOP algorithm, proposed an improved ACS version called the Enhanced Ant Colony
System (EACS) [22]. The two main changes were proposed. First, the construction phase of the
EACS used information about the best solution found so far. Second, the LS was run only if the
current solution was within 20% of the best found solution. The EACS was able to further improve
some of the best results obtained by Anghionlfi et al. [1] and to date remains one of the most
efficient methods for solving the SOP.
3. Improving ACS Convergence with Simulated Annealing
Ant colony optimization (ACO) is probably the best-known algorithm that was inspired by the
foraging behavior of ants in nature. It is a population-based meta-heuristic algorithm that is often
used to solve difficult combinatorial and continuous optimization problems. In general, it does not
guarantee that the optimal solution will be found, but solutions that are found are often of good
enough quality (for practical use) [14].
In the ACO, a number of artificial agents (ants) construct iteratively complete solutions to an
optimization problem. An ant starts with an empty solution and, in subsequent steps, extends it
with components selected from the set of all available components. Each component has an associated pheromone trail and a heuristic value. The higher the product of the pheromone concentration
(value) and the heuristic value is, the higher the probability that it will be selected by the ant.
In nature, ants communicate indirectly with one another by depositing small amounts of chemical
substances called pheromones, e.g., an ant that has found a food source marks the path to the nest
with small amounts of pheromone. The pheromone trail attracts other ants and leads them to the
food source. The more ants that repeat the process, the higher the concentration of the pheromone
trail becomes, hence the process becomes autocatalytic. The pheromone evaporates with time,
so the pheromone concentration does not increase indefinitely. The ACO algorithms use artificial
pheromone trails, with the pheromone concentration represented as real numbers. The set of all

4

pheromone trails is usually called a pheromone memory and plays a crucial role in the performance
of the ACO family of algorithms [12, 14].
For the TSP (and related problems) the problem is usually modeled by using a graph G(V, E).
An artificial ant constructs its solution starting from a randomly selected node. In subsequent steps
it moves from the current node to one of the unvisited neighbor nodes by using the corresponding
edge. The pheromone trails τuv are deposited on the edges, (u, v) ∈ E, of graph G and together
with a priori knowledge about the problem, reflected in the heuristic values associated with each
edge ηuv they guide the construction process.
The Ant Colony System is an improved version of the Ant System by Dorigo et al. [13]. In the
ACS the ant located at node i selects a next node j according to a pseudo-random proportional
rule [14]:
(
arg maxl∈Jki [τil ] · [ηil ]β , if q ≤ q0
j=
(1)
J,
if q > q0 ,
where ηil is a cost associated with an edge (i, l), τil is the value of the pheromone trail on edge
(i, l), Jki is a set of available (candidate) nodes of ant k, and q0 is a parameter, 0 ≤ q ≤ 1. J is a
node (city) selected according to the probability distribution defined by:
P (J|i) = P

[τiJ ] · [ηiJ ]β
.
β
l∈J i [τil ] · [ηil ]

(2)

k

The choice defined by Eq. 1 depends on the value of parameter q0 . If the randomly drawn
number q is lower than the parameter q0 , then the choice is greedy and the ant selects the node
to which an edge leads with the maximum product of pheromone trail τij and heuristic ηij values.
Otherwise q ≥ q0 and the choice is random with the probability distribution given by Eq. 2. The
first case is often referred to as exploitation of the knowledge gathered by the ants (in the pheromone
memory). Usually, a value of q0 close to 1 (often 0.9 and above) leads to good quality results in
a shorter period of time as compared to the base ACO algorithm [14]. Some authors even use a
higher value calculated as follows: q0 = 1 − |Vs | , where parameter s is the number of nodes that
should be selected randomly with the probability defined by Eq. 2 [21].
During construction of the solutions the ants in the ACS update the values of the pheromone
trails on the traversed edges. Each ant, after making a move from node u to node v, applies a local
pheromone update rule that decreases the amount of pheromone on edge (a, b) according to:
τab ← (1 − ψ) · τab + ψ · τ0 ,

(3)

where ψ is a parameter regulating evaporation of the pheromone over time and τ0 is the initial
pheromone level. The rationale behind formula (3) is that it lowers the probability of selecting the
same nodes by subsequent ants, hence it increases variety in the constructed solutions.
The global pheromone update performed after the ants have completed construction of their
solutions is more important. The update rule results in the increase in pheromone levels on trails
corresponding to the best solution found so far (Sbest ) and its value by Lbest . For each (u, v) ∈ Sbest ,
the pheromone changes according to the formula:
τuv ← (1 − ρ) · τuv + ρ · ∆τuv ,

(4)

where ∆τab = L−1
best and ρ ∈ (0, 1) is a parameter regulating the strength of the pheromone increase.
The global pheromone update ensures that edges belonging to the current best solution have higher
5

probabilities of being selected in the algorithm’s subsequent iterations. The global best solution
is used during the global pheromone update because it leads to slightly better solutions than the
iteration best solution [14].
In order to further shorten the computation time of Eq. 1, the so-called candidate set is used
which consists of the nearest neighbors of the current node. The size of the candidate set is usually
in the range of 10 to 25 [12, 14]. For comparison, the size n of the problem is often two or three
orders of magnitude larger, hence the use of candidate sets further limits exploration of the solution
search space. The candidate sets are a greedy heuristic based on the observation that good quality
solutions are comprised mainly of short edges. If all of the candidate set elements are already a part
of the constructed solution the ant selects one of the remaining (unvisited) nodes. The candidate
set is usually computed at the beginning and does not change. Randall and Montgomery [38]
investigated the idea of dynamic candidate set updates for the TSP and the Quadratic Assignment
Problem (QAP). The dynamic versions resulted in solutions of better quality but also significantly
increased the computation time of the whole algorithm.
3.1. Enhanced ACS
The Enhanced ACS algorithm proposed by Gambardella et al. is an efficient metaheuristic for
the SOP [22]. It differs from the ACS in two ways. The first is a modified solution construction
phase which is much more focused on the best solution found so far. Instead of direct application
of Eq. 1 an ant selects the node which follows the current node in the best solution so far (if
the random number q is lower than the parameter q0 ). Only if the node is already a part of the
constructed solution does the ant consider other nodes, i.e. it selects the edge with the maximum
product of the pheromone and heuristic values. If q ≥ q0 , the selection process from the ACS is
used. Parameter q0 usually has a value of 0.9 or higher, hence this modification significantly speeds
up the construction process although it limits the exploration capability of the EACS, and without
a strong LS, it achieves results of lower quality than the ACS [23].
The second modification is strong integration of the solution construction phase with the LS.
The LS is run only if the cost of the current solution is within 20% of the best solution found so
far. Also, the LS is initialized so that only elements of the current solution which are out of order
with respect to the best solution are placed on the so-called don’t push stack. The elements of the
stack are the starting points for the LS. This increases the emphasis on areas of the solution search
space that are potentially unexplored.
A slightly modified version of the EACS was proposed by Ezzat [19]. The main difference
concerns the choice of the next node in the solution construction process. At first it tries to select
the node v which follows the current node u in the best solution so far. If v is already a part
of the constructed solution, it selects the next node with the probability defined by Eq. 3. This
change favors exploration and makes the algorithm less exploitative than the EACS but still more
exploitative than the ACS. Later Ezzat et al. adapted the EigenAnt algorithm to solve the SOP [20].
The computational experiments showed that the proposed algorithm had performance comparable
to the EACS.
3.2. Simulated Annealing
Simulated Annealing is one of the most well-known general metaheuristic optimization methods.
It was inspired by the Monte Carlo method of sampling the states of a (physical) thermodynamic
system. In the SA, a solution to the optimized problem is equivalent to a state of the thermodynamic
system, and its quality corresponds to the system’s current energy [39]. The SA works as follows:
6

starting from an initial solution X0 , a sequence of solutions (Xi ), X ∈ S is generated, where S is a
set of all feasible solutions. Given a current solution Xi , a candidate solution Yi is generated and
its cost C(Yi ) is calculated. The next solution Xi+1 is selected according to:


Yi , if C(Yi ) < C(Xi ),
Xi+1 = Yi , with probability pi if C(Yi ) ≥ C(Xi ),
(5)


Xi , otherwise.
Probability pi is defined as pi = exp (−(C(Yi ) − C(Xi ))/Ti ), where Ti > 0 is called a temperature.
The physical analogy on which the SA is based requires that the system be kept close to a thermal
equilibrium as the temperature is lowered.
The most often used cooling schedule is the exponential schedule of the form: Ti+1 = λTi , where
λ is a parameter. In fact, the exponential cooling schedule usually lowers the temperature too fast
for the system to reach a near-equilibrium state and does not guarantee convergence to the global
optimum. Nevertheless, it is useful in practice because it is easy to implement and often leads to
good quality solutions if the computation time is limited [8]. More advanced cooling schedules have
been proposed; the two well-known ones are the adaptive cooling schedule by Huang et al. [39] and
the efficient cooling schedule by Lam [33].
3.3. Combining ACS with Simulated Annealing
The ACS generally offers a better convergence speed than the Ant System or ACO [12]. This
stems, among others, from the more exploitative solution construction process and the global
pheromone update rule that places emphasis on the best solution found so far. This usually speeds
up the process of finding good quality solutions but also makes escaping local minima very difficult.
Simulated Annealing, on the other hand, offers a simple solution to escape the local minima. We
propose how to combine the ACS and SA to enhance the ACS search process while maintaining
its exploitation oriented nature. The proposed algorithm, ACS with the SA (ACS-SA in short),
can be summarized as follows. The ACS search process is guided (in part) by the pheromone
trail values. At the end of each iteration the global pheromone update rule increases the values of
the pheromone trails corresponding to the components (edges) of the current best solution (global
best). In the proposed ACS-SA algorithm the global update rule uses instead an active solution
which may not necessarily be the best solution found so far. At the end of every iteration each of
the solutions generated by the ants is compared with the active solution. If the new solution is of
better quality, it replaces the current active solution. Otherwise, the new solution may still replace
the active solution but with a probability defined by the Metropolis criterion known from the SA.
While the ACS is always focused on the neighborhood of the best solution found so far and can
become trapped in a local optimum for a long time, the proposed ACS-SA has a greater chance of
escaping the local optima by shifting focus to the solution with a higher cost.
Figure 1 presents a pseudocode of the proposed ACS-SA algorithm. The major part of the algorithm does not differ from the ACS, i.e. the only differences are related to temperature initialization
(line 1), the cooling schedule (line 26) and the active solution selection process (line 19). Inclusion
of the SA into the ACS results in a more exploratory search process, but it may also lead to a
prolonged examination of areas of the solution space that contain solutions of a poor quality. This
is prevented by allowing the current global best solution to be selected as the active solution with
a probability of 0.1 (line 20). This heuristic might not be necessary if a more advanced cooling
schedule is used. The present work is intended to be proof of the concept that the SA may be used
7

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

T ← SA_compute_initial_temperature(γ)
// Calculate T0
for i ← 1 to #iterations do
for j ← 1 to #ants do
routeAnt(j) [1] ← U {1, #nodes}
// Start from randomly selected nodes
end
for k = 2 to #nodes do
// Build complete solutions
for j ← 1 to #ants do
routeAnt(j) [k] ← select_next_node(routeAnt(j) )
local_pheromone_update(routeAnt(j) [k − 1], routeAnt(j) [k])
end
end
for j ← 1 to #ants do
// Local update on the closing edges
local_pheromone_update(routeAnt(j) [#nodes], routeAnt(j) [1])
end
local_best ← select_best(routeAnt(1) , routeAnt(2) , . . . , routeAnt(#ants) )
if Cost(local_best) < Cost(global_best) then
global_best ← local_best
end
SA_select_solution(active_solution, [routeAnt(1) , routeAnt(2) , . . . , routeAnt(#ants) ])
if U (0, 1) < 0.1 then
// Allow a greedy ACS update with a small probability
global_pheromone_update(global_best)
else
global_pheromone_update(active_solution)
end
end
T ←T ·λ
// Lower annealing temperature
Figure 1: Ant Colony System combined with Simulated Annealing (ACS-SA).

8

1
2
3
4
5
6
7
8
9
10
11

Procedure SA_select_solution(active_solution,solutions)
for i ← 1 to #ants do
if Cost(solution[i]) < Cost(active_solution) then
active_solution ← solution[i]
else
∆C ← Cost(solution[i]) − Cost(active_solution)
if U (0, 1) < e−∆C/T then
active_solution ← solution[i]
end
end
end
Figure 2: Simulated Annealing-related procedure to update the active_solution.

to improve the convergence speed of the ACS, hence the geometric cooling schedule was adapted
for its simplicity. In future work a more advanced schedule, e.g. an adaptive cooling schedule by
Lam [33], could be applied.
Figure 2 presents the active solution selection procedure. The process iterates over a set of
solutions built by the ants. If the cost of an ant’s solution is lower than the cost of the active
solution, it replaces the active solution (lines 3–5 in Fig. 2). Otherwise, the ant’s solution (of
a worse quality) may replace the active solution with a probability calculated according to the
Metropolis criterion from the SA (lines 6–7). As the temperature is lowered, the probability of
accepting a worse solution goes down to 0 and the process becomes equivalent to that of the ACS.
The initial temperature T0 plays an important role in the SA. In our work we applied the idea of
an adaptive temperature calculation which was proposed in [3]. The calculation requires a sample
of randomly generated solutions whose values (costs) are used to calculate the initial temperature
according to:
∆C + 3σ∆C
T0 =
,
(6)
ln (1/γ)
where ∆C is the mean of absolute differences between the costs of consecutive pairs of solutions from
the sample, σ∆C is the sample standard deviation and γ is a parameter denoting the probability of
accepting a worse solution, i.e. with a higher cost. The idea behind Eq. 6 is based on the central
limit theorem which states that the mean of a large sample of independent random variables is
approximately normally distributed, hence, almost all (approx. 99.7%) absolute differences between
the quality of randomly generated solutions fall in the range of (∆C − 3σ∆C , ∆C + 3σ∆C ). Knowing
the approximation of the highest difference in quality between a pair of solutions allows to calculate
the initial temperature so that the probability of accepting a worse solution is γ.
Although the temperature initialization requires additional computations, it does not increase
the asymptotic complexity of the ACS algorithm. In our experiments a sample of 1000 random
solutions was used due to a negligible additional cost; however, a much smaller number could also
be acceptable.
3.4. Combining Enhanced ACS with the SA
As described in Sec. 3.1, the EACS differs only slightly from the ACS. The differences are minor
and concern the solution construction process and the LS application, hence it is straightforward
to apply exactly the same ideas to incorporate the SA into the EACS as in the proposed ACS-SA
algorithm. Due to its more (i.e. relative to the ACS) exploitative behavior, the EACS is even
9

more susceptible to getting trapped in the local minima, hence it should also benefit from the SA
component. The resulting algorithm will henceforth be denoted as the EACS-SA.
4. Efficient Local Search for the SOP
Even though the ACS, MMAS and related algorithms perform competitively to other nature inspired metaheuristics their convergence can still be improved with a problem-specific local
search [14]. When combined with the LS, the ACS is responsible for finding a candidate solution,
while the aim of the LS is to improve it by performing small changes leading to a neighboring
solution of a better quality. In this section we start with a description of the state-of-the-art LS
heuristic for the SOP and later propose a modified version which incorporates the SA component.
4.1. SOP-3-exchange
Gambardella et al. [21] proposed an efficient LS heuristic for the SOP called the SOP-3-exchange.
It adapts the 3-opt heuristic known from the TSP to the SOP without an increase in algorithm
time complexity. The SOP-3-exchange belongs to the family of edge-exchange procedures, in which
a new solution is generated by replacing k existing edges with another set of k edges for which the
cost of the solution is lower. This operation is usually called k-exchange, and the value of k can be
fixed (typically 2 or 3) or can vary as in the Lin-Kernighan heuristic [28]. Starting from the initial
solution and applying the k-exchange iteratively until no further improving exchange exists leads
to a k-optimal solution. This process requires in the worst-case scenario, O(nk ) time.
During a k-exchange procedure k existing edges are removed producing k disjoined paths which
are then reconnected with k new edges. In some cases, reconnection of the paths requires that
some of them be reversed, e.g. in the case of a 2-opt move and a closed path < 0, . . . , i − 1, i, i +
1, . . . , h − 1, h, h + 1, . . . , n − 1 >; there are two possible ways to reconnect the subpaths after
removal of the (i, i + 1) and (h, h + 1) edges, namely < . . . , h + 1, i + 1, . . . , h − 1, h, i, i − 1, . . . > and
< . . . , i − 1, i, h, h − 1, . . . , i + 1, h + 1, . . . >; both require a reversal of the subpath. The reversal,
however, is problematic for the SOP because the distances between the nodes are asymmetric, hence
the length of the path after the reversal should be recalculated what requires O(n) time. Because
the cost of a k-opt move should be calculated in a constant time an efficient implementation of the
k-opt heuristic for the SOP should be restricted only to path-preserving exchanges [21].
The smallest k that allows a path-preserving exchange is k = 3, denoted as the path-preserving3-exchange shown in Fig. 3. By removing the (h, h + 1), (i, i + 1) and (j, j + 1) edges and adding
(h, i + 1), (j, h + 1) and (i, j + 1) edges the two neighboring subpaths are swapped, thus preserving
the relative order of the elements. After performing the path-preserving-3-exchange one would still
need to verify if the precedence constraints for the two subpaths are preserved. This requires O(n2 )
time in the general case but can be avoided as in the method proposed by Gambardella et al. [21].
There are two necessary procedures to reduce the computation time. The first requires keeping the
lexicographic order while searching for the path-preserving-3-exchange. The second, is the use of a
labeling method.
Figure 3 shows how the route changes when applying the path-preserving-3-exchange. The tree
indices h, i, and j (h < i < j) define two sub-paths in the route: left (h + 1, i) and right (i + 1, j).
The subpaths are swapped as a result of performing the exchange, i.e. the right path comes before
the left path. This can only happen if there are no precedence constraints between the considered
node and the nodes in the left path.

10

The path-preserving-3-exchange as proposed by Gambardella et al. [21] consists of forward and
backward searches for feasible path-preserving-3-exchanges. The forward search involves incrementing j iteratively, thus increasing the length of the right path by one. This requires that only the
precedence constraints be checked between the elements of the left path and the node considered
for inclusion into the right path. Eventually, a precedence constraint is hit and the procedure is
repeated with the left path being extended with a single element (by incrementing i) and the right
path set to a single element, i.e. (j), j = i + 1. After all of the possibilities are exhausted h is
incremented and the process repeats for all possible i and j values (i < j < n). This leads to O(n3 )
possible pairs of subpaths each requiring O(n) constraints verification, hence a total complexity of
O(n4 ).
The cost of constraints verification can be reduced to O(1) due to the labeling procedure. The
procedure works as follows. Each time the left subpath is extended with a new node u (during
the sop-3-exchange), mark(v) is set to count for every node v for which there exists a precedence
constraint between u and v. The count is a variable initially set to 0 and incremented each time
the left path grows, i.e. h is incremented. Thanks to the procedure, each time the right path is
extended with a node x one needs only to check the value of mark(x). If the value equals the count,
then the node at index j (in the right path) has to be visited after the nodes in the left path, hence
the two paths cannot be swapped. This reduces the complexity of the whole search for a feasible
path-preserving-3-exchange to O(n3 ), which is asymptotically equal to the complexity of a 3-opt
heuristic used to solve the TSP.
The forward search for the path-preserving-3-exchange considers only exchanges defined by
indices i, j and h such that 0 < h < i < j < n, where n is the number of nodes. The backward
search is analogous to the forward search but the left and right paths are expanded in the direction
of decreasing indices, i.e. the left path "moves" from the end of the sequence to the beginning.
Summarizing, the time complexity of finding a single profitable path-preserving-3-exchange using the described procedure is O(n3 ). This is still expensive as the procedure is applied (to a
single solution) in a loop until no further improving move is found and it has to be repeated for
the subsequent solutions. Gambardella et al. [21] proposed two additional changes to reduce the
algorithm’s computation time. The first is to limit the search to only a subset of all potential moves.
By default the SOP-3-exchange for each index h considers all valid i and j indices. Assuming most
of the changes will involve relatively short paths, the i values can be restricted to h + 1, h + 2, h + 3
for the forward procedure and h − 1, h − 2, h − 3 for the backward procedure. This version was
named OR-exchange [21]. The second change involves the use of two additional heuristics, i.e. don’t
look bits and don’t push stack. Don’t look bits is a data structure that was proposed by Bentley [6]
which works as follows. A bit is associated with each node of the solution. At the beginning all
the bits are turned off and are turned on when the SOP-3-exchange starts looking for a profitable
exchange originating from the node. If the don’t look bit is turned on, the corresponding node is
ignored by the subsequent SOP-3-exchange searches until the node is involved in a profitable pathpreserving-3-exchange. Then all six pivot nodes (h, h + 1, i, i + 1, j, j + 1) are turned off. Use of the
don’t look bits aims to focus the search on the changing parts of the solution. The purpose behind
the use of the don’t push stack is similar – it contains nodes h to be selected as a starting point of a
path-preserving-3-exchange. At the beginning the stack is initialized with all of the nodes. During
the search a node, h, is removed from the stack and if the feasible move originating from node h is
found the six nodes involved in the exchange are pushed onto the stack (if they do not belong to it
already). An additional benefit of using the don’t push stack is that the linear order in which the
nodes are considered during the search for a profitable path-preserving-3-exchange is broken.

11

left
0

h

h+1

right
i

i+1

j

j+1

n

Figure 3: Path-preserving-3-exchange for the SOP

1
2
3
4
5
6
7

Procedure: SOP_3_exchange
Input
: route[0..n]
Input
: best_route[0..n]
S ← init_dont_push_stack(route, best_route)
while S 6= ∅ do
h ← pop_stack(S)
if forward_SOP_3_exchange(h,route) 6= True then
backward_SOP_3_exchange(h,route)
end
end

// A solution to improve
// The current best solution

Figure 4: Pseudocode of the SOP-3-exchange procedure.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

Function: forward_SOP_3_exchange
Input
: h ∈ {0, 1, . . . , n − 2}
// An index of the SOP-3-exchange start element
Input
: route[0..n]
// A SOP solution being improved
Output : True if an acceptable change was found, False otherwise
succeded ← False
label[0..n] ← ∅
for i ← h+1 to n-1 do
foreach node u preceded by route[i] do
// Labelling procedure
label[u] = h
// Node u comes after route[i]
end
best_delta ← 0
best_j ← ∅
j ← i+1
while label[j] 6= h do
delta ← calculate_cost_change(h,i,j)
if is_move_accepted(delta,best_delta) then
best_j ← j
best_delta ← delta
end
j ← j+1
end
if best_j 6= ∅ then
perform_exchange(h, i, best_j)
// Exchange paths (h + 1, i) and (i + 1, best_j)
succeded ← True
break
end
end

Figure 5: Forward part of the SOP-3-exchange algorithm that searches for a feasible move in the forward direction
starting at a given node h. The search in the backward direction is analogous with indices i and j being decreased
instead of increasing.

12

Figure 4 presents a pseudocode of the SOP-3-exchange algorithm. It starts with an initialization
of the don’t push stack and repeatedly searches for a feasible move. First, it searches in the forward
direction (line 4) and if it fails, the backward search is applied (line 5). The pseudocode of the search
in the forward direction is shown in Fig. 5 (the search in the backward direction is analogous). It
starts with a given index h that denotes the starting point of a possible path-preserving-3-exchange
and searches for the remaining two points, denoted by indices i and j. The labeling procedure is
applied incrementally (lines 4–6). The function is_move_accepted in line 12 simply checks if the
proposed decrease in the solution value is better than the current best, but it can be replaced by a
more advanced criterion as will be shown later.
4.2. Improving SOP-3-exchange Efficiency with SA
The SOP-3-exchange LS is efficient in improving solutions generated by the ants; however, the
improvement process is greedy and only better (downhill) moves are accepted. It makes it possible
to reach a local optimum quickly; however, it also makes it unable to find any better solution that
would require making at least one uphill move. Similarly to our idea of incorporating the SA into
the ACS and EACS algorithms, we propose to include the SA decision process into the SOP-3exchange in order to make it more explorative. The proposed modification is easy to implement as
it only requires to modify the greedy condition as to whether to accept a given subpath exchange in
the forward search for a path-preserving-3-exchange (line 12 in Fig. 5) (analogously in the backward
search). The pseudocode of the proposed modification is shown in Fig. 6. The decision whether
to accept the proposed move (subpath exchange) is made based on the change (decrease) in the
solution value and the value of the best move found so far. If the proposed move is better than
the current best, it is always accepted. Otherwise, if it results in the same decrease of the solution
length then it is accepted with a probability of 10% (lines 4–5 in Fig. 6). It allows to accept moves
which do not change the solution value but which result in a different relative order of the solution
nodes. Finally, if the proposed move is worse than the best move found so far, it is accepted with
the probability calculated using the Metropolis criterion, as in the SA.
Similarly to the ACS-SA, there are two parameters related to the SA component of the proposed
SOP-3-exchange-SA algorithm, namely λLS and γLS . The former is used in the geometric cooling
schedule to lower the temperature TLS , while γLS is related to an initial probability of accepting
a worse move. There is, however, a slight difference in the temperature initialization relative to
the ACS-SA. In the ACS-SA the initial temperature is calculated based on a sample of differences
in the quality (length) of the randomly generated solutions, just before the main computations.
In the SOP-3-exchange-SA the sample comes from the values of the differences in the solution
quality (delta values) resulting from the subsequent path-preserving-3-exchanges considered during
the initial runs of the SOP-3-exchange-SA (lines 11–16 in Fig. 6). In other words, there is no
dedicated temperature initialization phase in the SOP-3-exchange-SA and the sample of delta values
is collected on the run in order not to slow down the whole algorithm. After the sample of 105
(a value found experimentally) is collected, the initial value of temperature TLS is calculated, and
in subsequent invocations of the SOP-3-exchange-SA the temperature is reset to this initial value
without recalculating.
5. Computational Experiments
A series of computational experiments was conducted in order to evaluate the performance
of the proposed algorithms. In the first part of the experiments we focused on the efficiency of
13

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

Function: is_move_accepted_SA
Input: change
// A decrease in a solution length if the move is applied
Input: best_change
// The value of the best move found so far
Output: True if the move value should be accepted, False
∆ ← change − best_change
// ∆ is used to compare the proposed move with the current best
accept ← False
if ∆ > 0 then
accept ← True
// Always accept a better move
else if ∆ = 0 and U (0, 1) < 0.1 then
// Move does not change solution value
accept ← True
// but changes relative order of the nodes
else
// Worse move, apply the Simulated Annealing criterion
if temperature TLS was initialized then
set accept to True with prob. e−∆/TLS
TLS ← λLS TLS
// Lower the temperature
else temperature was not yet initialized
D ← D∪∆
// Collect a sample of ∆ values
if |D| ≥ 105 then
// Calc. the initial temp. if the collected sample is large enough
TLS ← calc_initial_temperature(D)
end
end

Figure 6: Simulated Annealing-related criterion used to accept a solution change (move) in the SOP-3-exchange-SA
local search algorithm.

the ACS-SA and EACS-SA used alone, i.e. without the problem-specific LS. In the second part
the focus was placed on the efficiency of the algorithms coupled with the SOP-3-exchange and
SOP-3-exchange-SA LS heuristics.
The ACS and EACS require that a number of parameters be set. Based on preliminary computations and suggestions from the literature we used the following settings in our experiments:
number of ants, m = 10; β = 0.5; ψ = 0.01 and ρ = 0.1, and local and global pheromone evaporation ratio, respectively; q0 = n−20
n , where n is the size of the problem. The computations were
repeated 30 times for each configuration of the parameter values and the problem instance. The
computations were conducted on a machine equipped with a Xeon E5-2680v3 12 core CPU clocked
at 2.5GHz, although a single core was used per run. All algorithms were implemented in C++ and
compiled with the GNU compiler with the -Ofast switch1 .
5.1. ACS-SA Parameter Tuning
The first part of the experiments was focused on the behavior of the ACS-SA algorithm depending on the values of the SA-related parameters. The proposed ACS-SA algorithm uses a simple
exponential cooling schedule Tk = T0 · λk , where λ < 1 is the cooling factor and T0 is the initial temperature. Although the exponential cooling schedule does not guarantee convergence to a
global optimum, it has the advantage of being easy to implement and often performs well in practice [31]. Preliminary computations showed that the most important factor for the performance
of the ACS-SA was the λ parameter which directly influences the speed of the SA convergence.
The best performance was observed for λ ≥ 0.999, for which the probability of accepting worse
quality solutions and, hence, escaping local minima remained high for a relatively long time. It is
not without significance that the algorithm was run for 105 iterations, and for shorter/longer runs
1 The

source code is available at https://github.com/RSkinderowicz/AntColonySystemSA

14

Table 1: Table containing the p-values of the post-hoc pairwise comparison between the results of the ACS-SA with
various (λ, γ) values (shown in the second row) according to the non-parametric, two-sided multiple comparison
procedure by Mack and Skillings at α = 0.05 [30]. The test itself corrects for the Type I family-wise error. The +/symbol after a value denotes that the configuration in a row was significantly better/worse than the configuration in
a column.

A
B
C
D
E
F
G
H
I

A
(0.999, 0.1)

B
(0.999, 0.5)

C
(0.999, 0.9)

D
(0.9995, 0.1)

E
(0.9995, 0.5)

F
(0.9995, 0.9)

G
(0.9999, 0.1)

H
(0.9999, 0.5)

I
(0.9999, 0.9)

–
0.9970
1.0
0.3863
0.9740
0.9595
0.0004+
0.4397
0.0205-

0.9970
–
0.9765
0.8856
1.0
1.0
0.0093+
0.9155
0.0010-

1.0
0.9765
–
0.2235
0.9019
0.8678
0.0001+
0.2639
0.0509

0.3863
0.8856
0.2235
–
0.9703
0.9817
0.4191.0
1.0
< 0.0001−

0.9740
1.0
0.9019
0.9703
–
1.0
0.0265+
0.9813
0.0002-

0.9595
1.0
0.8678
0.9817
1.0
–
0.0346+
0.9891.0
0.0002-

0.00040.00930.00010.4191.0
0.02650.0346–
0.3668
< 0.0001−

0.4397
0.9155
0.2639
1.0
0.9813
0.9891.0
0.3668
–
< 0.0001−

0.0205+
0.0010+
0.0509
< 0.0001+
0.0002+
0.0002+
< 0.0001+
< 0.0001+
–

a smaller/higher λ value could prove better. In fact, the λ value could be calculated based on the
total number of iterations if used in practice [31]. A number of "promising" values was selected
for a more thorough investigation, namely λ ∈ {0.999, 0.9995, 0.9999}. The initial temperature T0
was calculated for each problem instance during the initialization phase so that the probability of
accepting a worse solution (an uphill move) at the beginning was approx. equal to the specified
probability γ (a parameter, independent of a problem instance). The mean difference between the
successive solutions was estimated based on a sample of randomly generated solutions. In our experiments we considered γ ∈ {0.1, 0.5, 0.9} leading to a total of 9 combinations of λ and γ. The
algorithm was run for a total of 14 SOP instances from the TSPLIB repository, namely: ft53.1,
ft53.4, ft53.3, ft53.2, ft70.4, ft70.3, ft70.2, ft70.1, prob.100, kro124p.4, kro124p.3, kro124p.2 and
kro124p.1.
We used statistical tests to verify if the results for the various values of parameters differed
significantly. The proposed experimental design can be viewed as a two-way (two-factor) layout in
which the main factor is the combination of λ and γ values, while the second factor (also called a
blocking factor) is the problem instance (13 instances in our case) [30]. More specifically, the design
can be described as a randomized block design with an equal number of replications per treatmentblock combination. A suitable non-parametric (distribution-free) statistical test was proposed by
Mack and Skillings and is an equivalent of a parametric two-way ANOVA [34]. The null hypothesis,
H0 , which is of our interest here is that of no differences in the medians (of the solution quality)
for algorithms with various λ and γ values considered here (a total of 9 combinations). Rejecting
the null hypothesis would mean that the different values of the λ and γ parameters lead various
performance of the ACS-SA. The test requires that the Mack-Skillings statistic be computed (MS )
which is then compared with a critical value msα at the α level of significance (α = 0.05 in our
case) [30]. The null hypothesis H0 is rejected if MS ≥ msα . In our case MS ≈ 72.68 while the
critical value ms0.05 ≈ 15.23, hence H0 was rejected, providing rather strong evidence that the
values of λ and γ have a significant impact on the quality of the results generated by the ACS-SA.
This is an expected result because the value of λ should have a strong effect on the search trajectory
of the SA.
After the rejection of H0 , we can apply a post-hoc test to compare the individual pairs of algorithm results obtained for respective pairs of λ and γ values. A suitable asymptotically distributionfree, two-sided, multiple comparison procedure using within-block ranks was proposed by Mack and

15

22
20
18
16
14
12
10
8
6
4
500

ACS
ACS-SA*
ACS-SA**
ACS-SA***

20500

40500

60500

80500

100500

Mean relative error [%]

Mean relative error [%]

Skillings [30, 34]. Table 1 contains the final p-values of the pairwise comparison. As can be observed, in most cases there were no significant differences between the results of the ACS-SA with
the various λ and γ values. The only exception was configuration λ = 0.9999 and γ = 0.1, for which
the results were significantly better 6 out of 8 times. On the other hand, configuration λ = 0.9999
and γ = 0.9 was worse 7 out of 8 times. This shows that the SA component of the ACS-SA has the
strongest influence if the temperature is decreased slowly. It is important to properly adjust the
initial probability γ of accepting a worse quality solution and, hence, the initial temperature T0 . If
the probability is high, the algorithm easily accepts inferior solutions, particularly at the beginning
of the computations, and drifts away from the good quality solutions. It is worth emphasizing that
these observations are valid for the computation budget (time) used in the experiments; greatly
increasing the computation time could show even better convergence for higher γ values. Figure 7
shows the convergence plots for the ACS-SA with various λ and γ levels: for λ = 0.999 the temperature drops relatively quickly and convergence of the ACS-SA resembles that of the ACS. For
λ = 0.9999 the temperature drops more slowly and the algorithm has a greater chance of escaping
the local minima for a longer period of time. By increasing the initial temperature (as for γ = 0.9)
we can extend the initial "free wandering" phase at the expense of slower convergence.

120500

Iteration

20
18
16
14
12
10
8
6
4
2
500

ACS
ACS-SA*
ACS-SA**
ACS-SA***

10500

20500

30500

40500

50500

60500

Iteration

(a) ACS-SA convergence for the kro124p.2 instance.

(b) ACS-SA convergence for the ft53.4 instance.

Figure 7: Convergence of the ACS and ACS-SA for the two SOP instances. The ACS-SA was run with different
(λ, γ) values, namely: ACS-SA* – (0.9999, 0.1), ACS-SA** – (0.9999, 0.9) and ACS-SA*** – (0.999, 0.1). The first
500 iterations were skipped for clarity.

5.2. ACS-SA and EACS-SA Performance
The first experiment showed that the SA component indeed had a significant impact on ACSSA search convergence. In the subsequent experiment we focused on a comparison between the
ACS-SA relative to the ACS. We also considered the EACS and the EACS combined with the SA
(EACS-SA). Both the ACS-SA and the EACS-SA were run with λ = 0.9999 and γ = 0.1, chosen
based on the previous experiment. To make the comparison fair, all of the algorithms were run
with a time limit of 60 seconds of CPU time. Although the limit was relatively low it was sufficient
to detect differences in the performance of the algorithms. A total of 20 SOP instances of varying
size were selected from the TSPLIB repository.
The boxplots of the mean solution error are shown in Fig. 8 and Fig. 9. The differences between
the quality of the solutions generated by the algorithms are clearly visible. For the smaller instances,
performance of the ACS and EACS was relatively similar and, in most cases, worse than that of the

16

20
ACS
ACS-SA
EACS
EACS-SA

Mean relative error [%]

18
16
14
12
10
8
6
4
2

p.
2
kr

o1

24
o1
kr

24

p.
1

4
0.
ft7

3
0.
ft7

2
ft7

0.

1
0.
ft7

4
3.
ft5

3
3.
ft5

2
3.
ft5

ft5

3.

1

0

Problem

Figure 8: Boxplot of the mean solution error for the ACS, ACS-SA, EACS and EACS-SA algorithms for the (smaller)
SOP instances from the TSPLIB repository.
45
ACS
ACS-SA
EACS
EACS-SA

Mean relative error [%]

40
35
30
25
20
15
10
5

78
a
g3
rb

58
a
g3
rb

41
a
g3
rb

23
a
g3
rb

53
a
g2
rb

50
a
g1
rb

09
a
rb

g1

0
10
b.
pr
o

4p
.4
o1
2
kr

kr

o1
2

4p
.3

0

Problem

Figure 9: Boxplot of the mean solution error for the ACS, ACS-SA, EACS and EACS-SA algorithms for the (larger)
SOP instances from the TSPLIB repository.

ACS-SA and EACS-SA, respectively. The differences became more distinct for larger instances (up
to 380 nodes), for which the EACS outperformed the ACS in most cases. The ACS-SA generally beat
the ACS but even better performance was achieved by the EACS-SA version, particularly for the
largest instances. The results were compared statistically to make the comparison more complete.
For each problem instance, the Kruskal-Wallis non-parametric one-way analysis of variance test (an
extension of the Mann-Whitney U test) with α = 0.05 was applied to check the hypothesis that
the results of the four algorithms came from the same distribution. The hypothesis was rejected in
19 out of 20 cases meaning that the results of the algorithms differed significantly. In such cases
a post-hoc test was applied to compare all pairs of results. For this purpose the Bonferroni-Dunn
test was employed with a family-wise Type I error correction (αF W = 0.05) [40]. The results are
summarized in Tab. 2. For each pair of algorithms, only the final verdict is shown with a letter
indicating the algorithm that achieved significantly better results than the others. The ACS-SA
outperformed the ACS in 12 cases, while never generating worse results. The SA component is
17

Table 2: Summary of results obtained by the ACS, ACS-SA, EACS and EACS-SA algorithms on a set of 20 SOP
instances from the TSPLIB repository. The left-most part of the table contains the mean solution lengths along with
the standard deviations shown in the braces. The last 6 columns contain a summary of the statistical comparison
between the respective pairs of algorithms according to the two-sided, non-parametric Bonferroni-Dunn test with
a family-wise Type I error correction (αF W = 0.05). The capital letter indicates the algorithm which obtained
significantly better results than the others. Hyphens denote that there were no significant differences between the
results of the respective algorithms.
Problem

ACS (A)

ACS-SA (B)

EACS (C)

EACS-SA (D)

ft53.1
ft53.2
ft53.3
ft53.4
ft70.1
ft70.2
ft70.3
ft70.4
kro124p.1
kro124p.2
kro124p.3
kro124p.4
prob.100
rbg109a
rbg150a
rbg253a
rbg341a
rbg323a
rbg358a
rbg378a

7857 (161.8)
8713 (159.5)
11506 (578.9)
14744 (201.8)
40437 (458.0)
42263 (454.5)
44674 (667.0)
56098 (325.5)
42166 (757.8)
44548 (1113.2)
53915 (1832.6)
80373 (1120.1)
1485 (87.8)
1111 (9.9)
1872 (11.8)
3156 (15.1)
3103 (60.8)
3590 (29.9)
3284 (81.1)
3483 (54.0)

7673 (18.5)
8522 (107.0)
11417 (506.4)
14704 (81.7)
40054 (223.6)
41629 (409.0)
44388 (333.3)
56146 (127.0)
41313 (572.8)
43049 (764.4)
51411 (321.2)
79708 (922.3)
1489 (75.2)
1107 (11.4)
1855 (9.7)
3123 (13.6)
2989 (40.8)
3517 (21.0)
3128 (34.9)
3347 (38.9)

7818 (162.9)
8647 (256.8)
11271 (605.5)
14779 (128.2)
40692 (588.6)
42396 (664.4)
44589 (570.0)
55593 (564.1)
42324 (988.9)
44270 (1318.3)
53313 (1678.5)
81204 (1064.1)
1505 (76.6)
1093 (9.7)
1817 (11.6)
3101 (19.8)
2990 (37.2)
3540 (27.8)
3087 (52.0)
3442 (59.8)

7702 (46.1)
8348 (156.6)
11418 (544.0)
14639 (101.1)
40150 (345.5)
41710 (355.0)
43946 (436.6)
55305 (362.9)
41768 (731.4)
43529 (809.9)
51351 (963.8)
78973 (746.0)
1438 (66.8)
1067 (7.9)
1788 (10.5)
3041 (9.0)
2821 (31.5)
3393 (31.3)
2883 (36.3)
3184 (37.5)

A vs B

A vs C

A vs D

B vs C

B vs D

B
B
B
B
B
B
B
B
B
B
B
B

C
C
C
C
C
C
C
-

D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D

B
B
B
C
B
B
B
B
C
C
B

D
D
D
D
D
D
D
D
D
D
D
D

C vs D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D
D

actually more beneficial than the exploitation-oriented heuristics introduced in the EACS which
generated significantly better quality results only in 7 cases as compared to the ACS. The EACSSA outperformed both the EACS and the ACS in 16 out of 20 cases. The greatest difference could
be observed particularly for the larger SOP instances.
400000
ACS
ACS-SA
EACS
EACS-SA

350000

Iterations

300000
250000
200000
150000
100000
50000

a
78

a
g3
rb

58

a
g3
rb

41

a
g3
rb

23

a
g3
rb

53

a
g2
rb

50

a
g1
rb

09

0

g1
rb

10
b.

p.
4

pr
o

24

24
kr
o1

kr
o1

p.
3

0

Problem

Figure 10: Boxplot of the mean number of iterations made by the ACS, ACS-SA, EACS and EACS-SA algorithms
for SOP instances selected from the TSPLIB repository.

The Simulated Annealing component in both the ACS-SA and the EACS-SA does not increase

18

the asymptotic time complexity of the algorithms. Only the initial temperature calculation requires
a number of random solutions to be constructed, while the main ACS loop is little affected by the
Metropolis rule and the cooling schedule computations. Figure 10 shows the mean number of
iterations performed by each of the considered algorithms within a time limit of 60 sec. As can be
observed, the number of iterations depends mostly on the size of the problem instance, while the
differences between the algorithms are relatively small. The EACS and EACS-SA are faster than
the other two algorithms due to the less expensive solution construction process which builds a new
solution by reusing significant parts of a solution from the previous iteration.
5.3. SOP-3-exchange-SA Parameter Tuning
Similarly to the ACS-SA and EACS-SA the SOP-3-exchange-SA algorithm has two more,
SA-related, parameters, namely λLS and γLS . Based on preliminary computations, several values were preselected for further investigation, namely λLS ∈ {0.8, 0.9, 0.95, 0.99} and γLS ∈
{0.1, 0.5, 0.9}. All 12 combinations of the parameters values were considered. For each combination the EACS algorithm with the SOP-3-exchange-SA was run on a set of 14 instances of sizes
from 400 to 700 selected from the SOPLIB2006 repository, namely: R.400.100.15, R.400.100.30,
R.400.1000.15, R.400.1000.30, R.500.100.15, R.500.1000.1, R.500.1000.15, R.500.1000.30,
R.600.100.15, R.600.100.30, R.600.1000.15, R.700.100.30, R.700.1000.1, R.700.1000.15.
The non-parametric Mack-Skillings test was used to verify if there were any significant differences
between the results for the different λLS and γLS values, similarly to Sec. 5.1. The null hypothesis
H0 stating that there were no differences between the medians of the solutions’ quality produced
for the different parameter values was rejected if M S ≥ msα , where M S is the Mack-Skillings
statistic and msα is the critical value at specified level of significance α. In our case, M S ≈ 1840.75
and ms0.05 ≈ 19.66, hence H0 was rejected, providing strong evidence for the significant differences
between the quality of results of the EACS-SA with SOP-3-exchange-SA obtained for the different
λLS and γLS values.
A post-hoc multiple comparison test by Mack and Skillings [30] was applied to find out for
which values of the parameters the results were of better quality. Table 3 contains the computed
p-values, where a value at the intersection of the i-th row and j-th column denotes the p-value of the
comparison between the results obtained for values of λLS and γLS corresponding to the i-th row and
j-th column, respectively. An analysis of the test results revealed that for λLS = 0.99 and γLS = 0.1
the results were significantly better than for any other combination of values. Simultaneously, the
worst configuration, in terms of solution quality, was λLS = 0.8 and γLS = 0.1, hence the γLS
parameter is of lower significance than λLS , which directly influences the speed of the temperature
decrease in the SA component of the SOP-3-exchange-SA. Generally, the best results were obtained
for λLS equal to 0.95 and 0.99.
5.4. Comparison of algorithms
The last part of the experiments concerned the performance of ACS, ACS-SA, EACS and EACSSA combined with the LS algorithms, i.e. SOP-3-exchange and SOP-3-exchange-SA. This gives a
total of 8 algorithm combinations. To make the comparison fair, the algorithms were run with
the same time limit of 120 seconds and the same values of parameters (where appropriate). The
algorithms were run on SOP instances (48 in total) from the SOPLIB2006 repository [35]. It is
worth noting, that the EACS with the SOP-3-exchange is the current state-of-the-art metaheuristic
for the SOP [23].

19

Table 3: Table containing the p-values of the post-hoc pairwise comparison between the results of the EACS-SA
with the SOP-3-exchange-SA local search algorithm with various (λLS , γLS ) values according to the non-parametric,
two-sided multiple comparison procedure by Mack and Skillings at α = 0.05 [30]. The +/- symbol after a value
denotes that the results for the configuration in a row were significantly better/worse than those obtained for the
configuration in a column.

A
B
C
D
E
F
G
H
I
J
K
L

A
(0.8, 0.1)

B
(0.8, 0.5)

C
(0.8, 0.9)

D
(0.9, 0.1)

E
(0.9, 0.5)

F
(0.9, 0.9)

G
(0.95,0.1)

H
(0.95,0.5)

I
(0.95, 0.9)

J
(0.99, 0.1)

K
(0.99, 0.5)

L
(0.99, 0.9)

–
< 0.0001+
< 0.0001+
0.0131+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+

< 0.0001–
< 0.0001+
0.8878
0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+

< 0.0001< 0.0001–
< 0.00011.0000
< 0.0001+
1.0000
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+

0.01310.8878
< 0.0001+
–
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+

< 0.00010.00011.0000
< 0.0001–
< 0.0001+
0.9990
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+

< 0.0001< 0.0001< 0.0001< 0.0001< 0.0001–
< 0.00010.0309+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+

< 0.0001< 0.00011.0000
< 0.00010.9990
< 0.0001+
–
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+
< 0.0001+

< 0.0001< 0.0001< 0.0001< 0.0001< 0.00010.0309< 0.0001–
0.0077+
< 0.0001+
0.0014+
0.3709

< 0.0001< 0.0001< 0.0001< 0.0001< 0.0001< 0.0001< 0.00010.0077–
0.0073+
1.0000
0.9712

< 0.0001< 0.0001< 0.0001< 0.0001< 0.0001< 0.0001< 0.0001< 0.00010.0073–
0.0331< 0.0001-

< 0.0001< 0.0001< 0.0001< 0.0001< 0.0001< 0.0001< 0.00010.00141.0000
0.0331+
–
0.8270

< 0.0001< 0.0001< 0.0001< 0.0001< 0.0001< 0.0001< 0.00010.3709
0.9712
< 0.0001+
0.8270
–

Table 4: Table containing the p-values of the post-hoc pairwise comparison between the results of the ACS, ACS-SA,
EACS and EACS-SA algorithms according to the non-parametric, two-sided multiple comparison procedure by Mack
and Skillings at α = 0.05 [30]. The +/- symbol after a value denotes that the results for the algorithm in a row were
significantly better/worse than those obtained for the algorithm in a column. The LS1 and LS2 subscripts denote
the local search method used, i.e. the SOP-3-exchange and SOP-3-exchange-SA, respectively.
Algorithm
ACSLS1 (A)
ACSLS2 (B)
ACS-SALS1 (C)
ACS-SALS2 (D)
EACSLS1 (E)
EACSLS2 (F)
EACS-SALS1 (G)
EACS-SALS2 (H)

A

B

C

D

E

F

G

<0.0001+
<0.0001<0.0001<0.0001+
<0.0001+
<0.0001+
<0.0001+

<0.0001<0.0001<0.0001<0.0001+
<0.0001+
<0.0001+
<0.0001+

<0.0001+
<0.0001+
1.0000
<0.0001+
<0.0001+
<0.0001+
<0.0001+

<0.0001+
<0.0001+
1.0000
<0.0001+
<0.0001+
<0.0001+
<0.0001+

<0.0001<0.0001<0.0001<0.0001<0.0001+
0.9999
<0.0001+

<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001-

<0.0001<0.0001<0.0001<0.00010.9999
<0.0001+
<0.0001+

20

H
<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001+
<0.0001-

A quick analysis of the obtained results showed noticeable differences in the efficiency of the investigated algorithms. The experiment design allows to check for statistically significant differences
between the algorithms by using the non-parametric Mack-Skillings test for a two-factor layout.
The first factor is the algorithm that is applied while the second (blocking) factor is the SOP instance that is solved. The null hypothesis H0 of our interest is that there are no differences in
the quality of the solutions generated by the algorithms. The rejection of H0 would mean that
the algorithms differ in the quality of generated solutions. The critical value for the test at level
of significance equal to 0.05 is ms0.05 ≈ 14.03 and the Mack-Skillings statistic is M S ≈ 8843.11,
meaning that M S > ms0.05 , hence H0 was rejected.
Rejection of the null hypothesis allows us to apply a post-hoc test (also proposed by Mack
and Skillings [30]) to perform a pairwise comparison of the algorithms. The resulting p-values are
shown in Tab. 4. As can be observed, all values are either close to 0 or close to 1, meaning that
the differences are either sharp or nonexistent, respectively. Not surprisingly, all EACS variants obtained significantly better results than the ACS-based algorithms. The most efficient algorithm was
the EACS with the SOP-3-exchange-SA LS, which obtained results that were significantly better
than any of the other remaining algorithms. The second best was the EACS-SA with the SOP3-exchange-SA LS. Out of the four ACS variants the ACS with SOP-3-exchange-SA performed
better than the other three, thus confirming the efficiency of the proposed SOP-3-exchange-SA
LS. The worst performing were the ACS-SA with the SOP-3-exchange and the ACS-SA with the
SOP-3-exchange-SA. The poor performance of the ACS-SA variants can be explained by the weakened emphasis on the exploitation which admittedly increases the probability of escaping from
local optima but also slows the overall convergence of the algorithm, which is clearly visible if the
computational budget is modest, as in the experiment conducted here (120 sec.).
Even though some of the algorithms can be seen as generally more efficient than others, this is
not true in every case, as can be observed in Tab. 5, in which the sample mean and sample standard
deviation values are presented for the EACS and EACS-SA algorithms. The two most efficient, in
terms of solution quality, were the EACS with the SOP-3-exchange and the EACS-SA with the
SOP-3-exchange LS. While the former achieved lower mean values for more problem instances,
the latter performed particularly well for instances of a size up to 500. For the largest instances
(R.600.* and R.700.* ), the EACS with the SOP-3-exchange obtained the lowest mean values in 14
out of 16 cases. This suggests that the EACS-SA algorithm did not have enough time to converge
within the specified time limit.
Similar observations can be made from the analysis of the best solutions found by the algorithms
presented in Tab. 6. The table also contains the values of the best-known solutions; some of which
were obtained by Gouveia and Ruthmair by using an exact method (branch-and-cut) [26] and by
Papapanagiotou et al. [36], while the rest by metaheuristics, including the EACS with the SOP3-exchange [23]. For the 18 SOP instances, all four algorithms were able to find the best-known
solution at least once per 30 runs. For the 10 SOP instances new best solutions were found by
the proposed algorithms. The EACS with the SOP-3-exchange-SA found the new best solutions
for 6 instances, i.e. R.300.1000.15, R.500.100.15, R.500.100.15, R.600.1000.15, R.700.100.15, and
R.700.1000.15. The EACS-SA with the SOP-3-exchange-SA found the new best solutions for 4
instances, namely: R.300.100.15, R.400.100.15, R.400.1000.15 and R.500.1000.15. Overall, the
best known or improved solutions were obtained by at least one of the algorithms in 37 out of
48 cases (77%). All algorithms struggled most with instances in which the number of precedence
constraints was smallest, i.e. 1% (instances R.*.*.1 ) which suggests that there is still some room
for improvement of the LS algorithms.

21

Table 5: Sample mean and standard deviation values for the EACS and EACS-SA algorithms with the SOP-3exchange and SOP-3-exchange-SA LS heuristics obtained for the SOP instances from the SOPLIB2006 repository [35].
The smallest means in a row are written in bold.

Instance

EACS
+SOP-3-exchange
Avg.
Std. dev.

EACS
+SOP-3-exchange-SA
Avg.
Std. dev.

EACS-SA
+SOP-3-exchange
Avg.
Std. dev.

EACS-SA
+SOP-3-exchange-SA
Avg.
Std. dev.

R.200.100.1
R.200.100.15
R.200.100.30
R.200.100.60
R.200.1000.1
R.200.1000.15
R.200.1000.30
R.200.1000.60

74.5
1935.5
4216.0
71749.0
1457.2
21766.2
41196.0
71556.0

3.8
31.0
0.0
0.0
17.9
298.1
0.0
0.0

69.4
1837.1
4216.0
71749.0
1459.0
20771.4
41196.0
71556.0

2.4
25.0
0.0
0.0
21.5
201.0
0.0
0.0

67.6
1916.9
4216.0
71749.0
1433.9
21577.6
41196.0
71556.0

1.9
27.8
0.0
0.0
9.8
399.5
0.0
0.0

66.6
1828.6
4216.0
71749.0
1467.2
21104.0
41196.0
71556.0

2.2
56.5
0.0
0.0
12.2
881.0
0.0
0.0

R.300.100.1
R.300.100.15
R.300.100.30
R.300.100.60
R.300.1000.1
R.300.1000.15
R.300.1000.30
R.300.1000.60

45.2
3328.1
6124.2
9726.0
1424.0
31556.9
54223.3
109482.9

3.9
54.8
10.9
0.0
29.4
647.4
72.7
36.3

36.0
3177.4
6120.0
9726.0
1432.2
29713.0
54147.0
109471.0

2.7
18.7
0.0
0.0
33.3
441.0
0.0
0.0

39.7
3289.9
6122.3
9726.0
1436.7
31196.1
54176.5
109488.3

3.2
31.0
5.3
0.0
21.6
637.0
30.7
41.3

31.3
3209.9
6120.0
9726.0
1464.2
30002.0
54147.0
109471.0

2.4
97.4
0.0
0.0
24.3
1319.5
0.0
0.0

R.400.100.1
R.400.100.15
R.400.100.30
R.400.100.60
R.400.1000.1
R.400.1000.15
R.400.1000.30
R.400.1000.60

35.9
4270.9
8167.6
15228.0
1504.7
42436.7
85320.8
140816.0

4.8
67.9
9.9
0.0
29.9
632.7
117.5
0.0

23.6
3969.8
8165.0
15228.0
1528.6
39854.6
85157.5
140816.0

4.1
29.3
0.0
0.0
24.8
317.6
62.9
0.0

33.1
4184.8
8166.0
15228.0
1521.9
41258.3
85188.6
140816.0

4.6
52.5
0.2
0.0
22.6
516.9
77.1
0.0

19.0
3930.5
8165.0
15228.0
1561.3
39293.0
85128.6
140816.0

2.2
28.2
0.0
0.0
31.9
204.3
3.3
0.0

R.500.100.1
R.500.100.15
R.500.100.30
R.500.100.60
R.500.1000.1
R.500.1000.15
R.500.1000.30
R.500.1000.60

25.4
5801.2
9709.9
18255.3
1524.8
54644.0
99465.1
178247.6

3.3
91.3
20.4
5.2
38.8
760.2
209.0
108.6

11.6
5417.1
9668.3
18240.0
1556.9
50503.3
99038.0
178212.0

2.6
54.5
5.4
0.0
34.3
360.3
42.4
0.0

25.3
5717.8
9696.0
18253.4
1542.5
53719.0
99166.6
178268.9

4.4
86.4
13.2
7.7
34.1
632.9
78.9
130.1

9.9
5449.8
9740.5
18240.0
1591.3
50310.4
99022.0
178212.0

2.6
254.7
87.2
0.0
38.3
420.0
160.9
0.0

R.600.100.1
R.600.100.15
R.600.100.30
R.600.100.60
R.600.1000.1
R.600.1000.15
R.600.1000.30
R.600.1000.60

21.3
6280.1
12504.9
23293.0
1582.7
61492.4
127460.6
214608.0

4.0
114.5
15.7
0.0
36.0
964.8
384.8
0.0

4.9
5621.5
12467.9
23293.0
1679.3
56638.3
126798.6
214608.0

2.1
40.9
5.2
0.0
40.6
389.9
3.5
0.0

22.6
6224.2
12499.7
23293.0
1657.7
62859.3
128011.6
214641.1

2.9
100.5
11.8
0.0
104.0
3426.5
305.2
59.0

14.3
5797.5
12569.7
23293.0
1749.3
61581.3
129319.9
214608.0

8.8
511.5
34.9
0.0
47.2
7760.0
505.6
0.0

R.700.100.1
R.700.100.15
R.700.100.30
R.700.100.60
R.700.1000.1
R.700.1000.15
R.700.1000.30
R.700.1000.60

16.6
7835.3
14585.7
24108.7
1539.8
72271.5
135922.5
245602.3

3.5
110.4
28.2
7.3
45.9
934.4
441.5
49.5

1.7
7195.4
14510.3
24102.0
1644.1
66738.6
134495.7
245589.0

0.7
65.3
0.8
0.0
36.0
479.7
40.0
0.0

26.4
8964.0
14693.2
24131.3
1718.5
79815.4
136468.1
245848.2

6.0
775.0
36.5
13.6
129.8
7934.4
347.6
154.1

30.3
7865.6
14792.3
24102.0
1930.4
83061.5
138703.9
245589.0

7.5
854.7
37.7
0.0
83.8
9842.3
477.9
0.0

22

Table 6: Best solution values obtained by the EACS and EACS-SA algorithms with the SOP-3-exchange and SOP3-exchange-SA LS heuristics for the SOP instances from the SOPLIB2006 repository [35]. The smallest values in a
row were written in bold. Values improving upon the best results known from the literature were underlined.
Best known
[23, 26, 36]

EACS
+SOP-3-exchange

EACS
+SOP-3-exchange-SA

EACS-SA
+SOP-3-exchange

R.200.100.1
R.200.100.15
R.200.100.30
R.200.100.60
R.200.1000.1
R.200.1000.15
R.200.1000.30
R.200.1000.60

61
1792
4216
71749
1404
20481
41196
71556

67
1890
4216
71749
1426
21113
41196
71556

64
1796
4216
71749
1423
20481
41196
71556

63
1868
4216
71749
1416
20946
41196
71556

64
1792
4216
71749
1437
20481
41196
71556

R.300.100.1
R.300.100.15
R.300.100.30
R.300.100.60
R.300.1000.1
R.300.1000.15
R.300.1000.30
R.300.1000.60

26
3161
6120
9726
1294
29183
54147
109471

39
3251
6120
9726
1363
30295
54147
109471

31
3154
6120
9726
1382
29068
54147
109471

32
3235
6120
9726
1389
30099
54147
109471

28
3152
6120
9726
1413
29111
54147
109471

R.400.100.1
R.400.100.15
R.400.100.30
R.400.100.60
R.400.1000.1
R.400.1000.15
R.400.1000.30
R.400.1000.60

13
3906
8165
15228
1343
43268
85128
140816

28
4139
8165
15228
1459
41067
85128
140816

17
3918
8165
15228
1488
39148
85128
140816

26
4088
8165
15228
1491
40457
85128
140816

14
3883
8165
15228
1496
38963
85128
140816

R.500.100.1
R.500.100.15
R.500.100.30
R.500.100.60
R.500.1000.1
R.500.1000.15
R.500.1000.30
R.500.1000.60

4
5361
9665
18240
1316
50725
98987
178212

17
5615
9687
18240
1464
53082
99072
178212

8
5305
9665
18240
1487
49907
98987
178212

17
5575
9675
18240
1476
52630
99002
178212

6
5315
9665
18240
1519
49719
98987
178212

R.600.100.1
R.600.100.15
R.600.100.30
R.600.100.60
R.600.1000.1
R.600.1000.15
R.600.1000.30
R.600.1000.60

1
5684
12465
23293
1337
57237
126798
214608

15
6087
12481
23293
1500
59660
126798
214608

2
5548
12465
23293
1611
55725
126798
214608

16
6008
12484
23293
1518
60194
127531
214608

4
5568
12496
23293
1644
56283
128214
214608

R.700.100.1
R.700.100.15
R.700.100.30
R.700.100.60
R.700.1000.1
R.700.1000.15
R.700.1000.30
R.700.1000.60

1
7311
14510
24102
1231
66837
134474
245589

10
7629
14534
24102
1462
70219
134712
245589

1
7109
14510
24102
1583
66008
134474
245589

18
8011
14620
24114
1588
71484
135841
245589

16
7191
14714
24102
1713
68009
137784
245589

Instance

23

EACS-SA
+SOP-3-exchange-SA

The results as presented above confirm that the proposed incorporation of the SA into the main
algorithm (EACS) and into the local search (SOP-3-exchange) is able to improve the quality of
the generated solutions to the SOP. In order to further clarify the differences between the existing
approach, i.e. the EACS with the SOP-3-exchange LS, and the proposed EACS-SA with the SOP3-exchange-SA, both were run on SOP instances from the SOPLIB2006 repository, however the
computation time was increased to 600 seconds. This is a five-fold increase vs the time limit used
in the experiments presented above. By giving the algorithms more time, we lower the risk of one
algorithm dominating an other because of the limited time. The results are presented in Tab. 7. In
most cases the results of the EACS-SA with the SOP-3-exchange-SA were of a better quality than
those obtained for the EACS with the SOP-3-exchange, although the relative differences between
the algorithms depended on the SOP instance that was being solved. The results were checked
for a statistically significant differences using the non-parametric Wilcoxon rank sum test at a
significance level of α = 0.05 (the respective p-values are reported in the table). In 33 out of 48
(69%) cases (instances) the solutions generated by the EACS-SA with the SOP-3-exchange-SA were
of a significantly better quality than those generated by the EACS with the SOP-3-exchange. In
4 cases (8%) the results of the former algorithm were significantly worse and in 11 cases (23%) no
significant differences were observed.
Taking into account the best solutions generated during 30 executions of the algorithms for
each of the SOP instances considered, the proposed algorithm reached the best-known results in 31
cases, and in 10 cases new best solutions were found. Because of the increased computation time
limit, in 7 out of those 10 cases the results were improvement over those presented in Tab. 6. To
summarize, the best-known or improved results were generated for 41 out of the 48 (85%) SOP
instances considered here. The EACS with the SOP-3-exchange found the best known results in 18
cases; however, no new best solutions were found in any case.
All of the SOP instances for which the EACS-SA with the SOP-3-exchange-SA generated significantly worse results than the EACS with the SOP-3-exchange are of the form R.*.1000.1 what
suggests either an overall inferior convergence of the former algorithm for this kind of instances,
or an insufficient time limit to match the convergence of the latter algorithm. In fact, the second
possibility seems to be true because for the smallest of the R.*.1000.1 instances, i.e. R.200.1000.1,
the EACS-SA with the SOP-3-exchange-SA generated significantly better results, and for the second smallest instance, i.e. R.300.1000.1, there were no significant differences between the results
of the two algorithms. To confirm our assumption, both algorithms were run for the instances:
R.300.1000.1, R.400.1000.1, R.500.1000.1, R.600.1000.1, and R.700.1000.1 but with the time
limit increased to 1200 seconds (doubled) per run. The results are presented in Tab. 8. As can
be seen, the advantage of the EACS with the SOP-3-exchange over the EACS-SA with the SOP3-exchange-SA disappeared, and both algorithms generated results of a similar quality. Statistical
comparison based on the non-parametric Wilcoxon rank sum test showed no significant differences
for the R.400.1000.1, R.500.1000.1 and R.600.1000.1 instances. Surprisingly, the increased time
limit allowed the EACS-SA to obtain significantly better results for the two remaining instances,
i.e. R.300.1000.1, and R.700.1000.1, although the advantage was small relative to the optimum
values.
Considering all the results, the efficiency of the proposed algorithm (in terms of the quality of
solutions) was statistically significantly better than the original approach for approx. 73% of the
SOP instances, while never being worse. However, a sufficient computation time is necessary to
reach this level of performance. In most cases 600 seconds was enough, whereas for a few instances
the limit of 1200 seconds was necessary. In practical applications, the algorithm could be sped up

24

by using parallel computations.
5.5. Speed comparison

Mean iterations

300000
EACSLS1
EACSLS2
EACS-SALS1
EACS-SALS2

250000
200000
150000
100000
50000
0
200

300

400

500

600

700

Problem size

Figure 11: Average number of iterations performed by the EACS and EACS-SA algorithms vs the size of the SOP
instance. The LS1 and LS2 subscripts denote the local search method used, i.e. the SOP-3-exchange and SOP-3exchange-SA, respectively.

The algorithms differ not only in the quality of the generated solutions but also in the relative
speed. The SA component does not affect the asymptotic time complexity of the ACS and EACS but
it may influence the solution search "trajectory", thus possibly impacting the runtime, particularly
if a local search is used. The SOP-3-exchange tries to improve a solution by searching only for the
improving changes (moves) and its time complexity depends on the relative order of nodes in the
solution. If the solution changes slightly from iteration to iteration, the runtime shortens because of
the focusing only on the changed parts of the solution. In contrast, the SOP-3-exchange-SA, due to
the SA component, may also accept a number of worse (up-hill) moves, hence the overall runtime
should increase. Figure 11 shows a bar plot of the average number of iterations performed for
the EACS and EACS-SA with both LS variants vs the size of the SOP instance. As expected, the
algorithms with the SOP-3-exchange-SA were slower than the algorithms with the SOP-3-exchange.
The fastest algorithm was the EACS with the SOP-3-exchange, beating the EACS-SA with the same
LS. Interestingly, the slowest algorithm was the EACS with the SOP-3-exchange-SA; it was even
slower than the EACS-SA with the same LS. This is probably due to the fact that the EACS can
relatively easily get trapped in a "deep" local minimum from which an escape is difficult even if the
SOP-3-exchange-SA accepts a number of up-hill moves. On the other hand, the EACS-SA focuses
on a larger number of different solutions during the search, some of which are less time-consuming
to improve by the LS. Finally, the larger the size of the instance, the lower the number of iterations
performed by the algorithms.
6. Summary
The Ant Colony System and particularly its enhanced version (EACS) are competitive metaheuristics whose efficiency has been shown in a number of cases [14, 21, 23]. Nevertheless, we have
shown that the search process of the ACS and EACS can still be improved with ideas taken from
Simulated Annealing. Specifically, instead of increasing the pheromone values based on the current
best solution, the proposed ACS-SA and EACS-SA algorithms increase the pheromone values based
on the current active solution that is chosen from among all the solutions constructed by the ants.
The active solution may not necessarily be the current best solution as it is selected probabilistically
25

Table 7: Results of the EACS with the SOP-3-exchange (I) and EACS-SA with the SOP-3-exchange-SA (II) algorithms for SOPLIB2006 instances with the time limit set to 600 sec. Verdict denotes the algorithm for which
the obtained results were of a significantly better quality than the results of the other algorithm according to the
non-parametric Wilcoxon rank-sum test at a level of significance α = 0.05. Cases for which there was no significant
difference are marked with a "–".

Instance

Best known
[23, 26, 36]

EACS + SOP-3-exchange (I)
Avg.

Std. dev.

EACS-SA + SOP-3-exchange-SA (II)

Best

Avg.

Std. dev.

p-value

Verdict
I vs II

Best

R.200.100.1
R.200.100.15
R.200.100.30
R.200.100.60
R.200.1000.1
R.200.1000.15
R.200.1000.30
R.200.1000.60

61
1792
4216
71749
1404
20481
41196
71556

71.8
1915.7
4216
71749
1445.4
21586.1
41196
71556

2.8
33.6
0
0
15.3
373.6
0
0

67
1849
4216
71749
1408
20837
41196
71556

63
1821.6
4216
71749
1423.2
20639.2
41196
71556

0.9
47.3
0
0
6.8
200.2
0
0

62
1792
4216
71749
1407
20481
41196
71556

<0.0001
<0.0001
1
1
<0.0001
<0.0001
1
1

II
II
II
II
-

R.300.100.1
R.300.100.15
R.300.100.30
R.300.100.60
R.300.1000.1
R.300.1000.15
R.300.1000.30
R.300.1000.60

26
3161
6120
9726
1294
29183
54147
109471

43.8
3297.7
6120.9
9726
1398.3
31132.5
54179.3
109490.8

5
29.3
1
0
26.9
542.8
50.8
45.1

37
3238
6120
9726
1356
30013
54147
109471

27.8
3174.2
6120
9726
1394.1
29309.7
54147
109471

1.7
47.7
0
0
19.9
193.5
0
0

26
3152
6120
9726
1355
29026
54147
109471

<0.0001
<0.0001
<0.0001
1
0.5691
<0.0001
<0.0001
0.0214

II
II
II
II
II
II

R.400.100.1
R.400.100.15
R.400.100.30
R.400.100.60
R.400.1000.1
R.400.1000.15
R.400.1000.30
R.400.1000.60

13
3906
8165
15228
1343
43268
85128
140816

33.2
4184.1
8165.9
15228
1471.3
41821.2
85253.2
140816

4.2
50.4
0.3
0
24.7
504.7
105.6
0

25
4084
8165
15228
1426
40869
85146
140816

14.6
3902.6
8165
15228
1497
39207.1
85128
140816

1.9
12.6
0
0
25
148
0
0

13
3883
8165
15228
1454
38963
85128
140816

<0.0001
<0.0001
<0.0001
1
0.0002
<0.0001
<0.0001
1

II
II
II
I
II
II
-

R.500.100.1
R.500.100.15
R.500.100.30
R.500.100.60
R.500.1000.1
R.500.1000.15
R.500.1000.30
R.500.1000.60

4
5361
9665
18240
1316
50725
98987
178212

24
5739
9696.1
18255.9
1465.6
53818.4
99240.5
178259.5

4
76.8
9.1
4.3
30.9
726.2
117.9
123.1

13
5578
9687
18240
1416
52549
99018
178212

4.9
5348.7
9667.6
18240
1490.7
49815.3
98987
178212

0.9
27
1
0
31.5
197.8
0
0

4
5284
9665
18240
1438
49504
98987
178212

<0.0001
<0.0001
<0.0001
<0.0001
0.0043
<0.0001
<0.0001
0.0418

II
II
II
II
I
II
II
II

R.600.100.1
R.600.100.15
R.600.100.30
R.600.100.60
R.600.1000.1
R.600.1000.15
R.600.1000.30
R.600.1000.60

1
5684
12465
23293
1337
57237
126798
214608

17.4
6143.7
12484.6
23293
1515.7
60100.8
127106.1
214608

3
114.4
13.7
0
27.7
785.7
339.8
0

12
5919
12468
23293
1446
58426
126798
214608

1.7
5544
12465
23293
1540.2
55891.4
126798
214608

1
30.1
0
0
25.6
446.3
0
0

1
5493
12465
23293
1492
55213
126798
214608

<0.0001
<0.0001
<0.0001
1
0.001
<0.0001
<0.0001
1

II
II
II
I
II
II
-

R.700.100.1
R.700.100.15
R.700.100.30
R.700.100.60
R.700.1000.1
R.700.1000.15
R.700.1000.30
R.700.1000.60

1
7311
14510
24102
1231
66837
134474
245589

10.9
7707
14563.9
24103.7
1463
70778.6
135480
245589

2.2
90.4
24.7
4.3
34.6
784.6
339.9
0

7
7560
14516
24102
1382
69410
134942
245589

1
7125.3
14510
24102
1484.9
65935.9
134474
245589

0
52.3
0
0
29.7
436.9
0
0

1
7021
14510
24102
1423
65305
134474
245589

<0.0001
<0.0001
<0.0001
0.0215
0.028
<0.0001
<0.0001
1

II
II
II
II
I
II
II
-

26

Table 8: Results of the EACS with the SOP-3-exchange (I) and EACS-SA with the SOP-3-exchange-SA (II) algorithms for the selected SOPLIB2006 instances with the time limit set to 1200 sec. The meaning of the columns is as
before.

Instance
R.300.1000.1
R.400.1000.1
R.500.1000.1
R.600.1000.1
R.700.1000.1

Best known
[23, 26, 36]
1294
1343
1316
1337
1231

EACS + SOP-3-exchange (I)
Avg.
1397.7
1456.0
1463.1
1504.2
1445.6

Std. dev.
27.9
26.1
28.8
30.2
26.3

Best

EACS-SA + SOP-3-exchange-SA (II)
Avg.

1360
1414
1416
1455
1398

1361.3
1462.1
1459.1
1508.8
1430.2

Std. dev.
13.7
24.9
20.7
26.0
24.1

p-value

Verdict
I vs II

<0.0001
0.3254
0.7006
0.3477
0.0281

II
–
–
–
II

Best
1339
1421
1411
1447
1381

by using the Metropolis criterion from the SA. This change weakens the exploitative focus of the
ACS and EACS, thus increasing the chance of escaping local optima. The computational experiments on a set of SOP instances from the TSPLIB repository and subsequent statistical analyses
have shown that in most cases the resulting ACS-SA and EACS-SA algorithms perform significantly
better than the original algorithm.
An efficient local search heuristic is necessary for state-of-the-art performance in solving the
SOP. Based on the same SA inspirations, we proposed an enhanced version of the state-of-the-art
SOP-3-exchange heuristic by Gambardella [21]. The resulting SOP-3-exchange-SA algorithm is
more resilient to getting trapped in local minima, at the expense of increased computation time.
The computational experiments conducted on a set of 48 SOP instances sized from 200 to 700
showed that the proposed EACS and EACS-SA with the SOP-3-exchange and SOP-3-exchange-SA
local searches are in many cases able to find solutions of better quality than the original EACS
with the SOP-3-exchange (a current state-of-the-art metaheuristic for the SOP), within the same
computation time limit. In fact, new, best solutions were obtained for 10 instances. In total, the
best known or improved solutions were obtained at least once for a total of 85% of the SOP instances
considered here.
Although the proposed modifications are easy to implement and improve the performance of the
original algorithms, they have some minor drawbacks. First, they increase the computation time
relative to the original algorithms. Second, they require to set the values of the new parameters
related to the SA cooling schedule (λ and γ). Also, relatively poor performance for SOP instances
with a small number (1%) of precedence constraints shows that there is still room for improvement,
both in the ACS-SA, EACS-SA and local search methods.
In the future, a more advanced cooling schedule could be used to improve the convergence
of the SA component of the proposed algorithms. A good candidate seems to be the adaptive
cooling schedule that was proposed by Lam [33], although it requires a complex parameter setting
and a method of controlling how much the subsequent solutions differ from one another. An
interesting idea could also be to activate the SA component only if search process stagnation is
detected. Because the proposed fusion between the ACS and SA is problem-agnostic one could try
to apply it to solve other difficult combinatorial optimization problems. The performance of the
proposed algorithms in terms of computation time could also be improved with the help of parallel
computations, as the ACS is susceptible to parallelization even with modern GPUs [41].
Acknowledgments: This research was supported in part by PL-Grid Infrastructure.

27

References
[1] Davide Anghinolfi, Roberto Montemanni, Massimo Paolucci, and Luca Maria Gambardella. A hybrid particle
swarm optimization approach for the sequential ordering problem. Computers & OR, 38(7):1076–1085, 2011.
[2] Norbert Ascheuer. Hamiltonian path problems in the on-line optimization of flexible manufacturing systems.
PhD thesis, Technische Universität Berlin, 1995.
[3] Mir M. Atiqullah. An efficient simple cooling schedule for simulated annealing. In Antonio Laganà, Marina L.
Gavrilova, Vipin Kumar, Youngsong Mun, Chih Jeng Kenneth Tan, and Osvaldo Gervasi, editors, Computational Science and Its Applications - ICCSA 2004, International Conference, Assisi, Italy, May 14-17, 2004,
Proceedings, Part III, volume 3045 of Lecture Notes in Computer Science, pages 396–404. Springer, 2004.
[4] Masri Binti Ayob and Ghaith M. Jaradat. Hybrid ant colony systems for course timetabling problems. In
Proceedings of the 2nd Conference on Data Mining and Optimization, DMO 2009, Universiti Kebangsaan
Malaysia, 27-28 October 2009, pages 120–126. IEEE, 2009.
[5] J. Behnamian, M. Zandieh, and S.M.T. Fatemi Ghomi. Parallel-machine scheduling problems with sequencedependent setup times using an aco, {SA} and {VNS} hybrid algorithm. Expert Systems with Applications,
36(6):9637 – 9644, 2009.
[6] Jon Louis Bentley. Fast algorithms for geometric traveling salesman problems. INFORMS Journal on Computing, 4(4):387–411, 1992.
[7] Lyamine Bouhafs, Amir Hajjam, and Abder Koukam. A combination of simulated annealing and ant colony
system for the capacitated location-routing problem. In Bogdan Gabrys, Robert J. Howlett, and Lakhmi C. Jain,
editors, Knowledge-Based Intelligent Information and Engineering Systems, 10th International Conference,
KES 2006, Bournemouth, UK, October 9-11, 2006, Proceedings, Part I, volume 4251 of Lecture Notes in
Computer Science, pages 409–416. Springer, 2006.
[8] OLIVIER CATONI. Rough large deviation estimates for simulated annealing: Application to exponential
schedules. The Annals of Probability, 20(3):1109–1146, 1992.
[9] Chia-Ho Chen and Ching-Jung Ting. A hybrid ant colony system for vehicle routing problem with time windows.
Journal of the Eastern Asia Society for Transportation Studies, 6:2822–2836, 2005.
[10] Shyi-Ming Chen and Chih-Yao Chien. Solving the traveling salesman problem based on the genetic simulated
annealing ant colony system with particle swarm optimization techniques. Expert Syst. Appl., 38(12):14439–
14450, 2011.
[11] Zbigniew J. Czech, Wojciech Mikanik, and Rafal Skinderowicz. Implementing a parallel simulated annealing
algorithm. In Roman Wyrzykowski, Jack Dongarra, Konrad Karczewski, and Jerzy Wasniewski, editors, Parallel
Processing and Applied Mathematics, 8th International Conference, PPAM 2009, Wroclaw, Poland, September
13-16, 2009. Revised Selected Papers, Part I, volume 6067 of Lecture Notes in Computer Science, pages 146–155.
Springer, 2009.
[12] Marco Dorigo and Luca Maria Gambardella. Ant colony system: a cooperative learning approach to the traveling
salesman problem. IEEE Trans. Evolutionary Computation, 1(1):53–66, 1997.
[13] Marco Dorigo, Vittorio Maniezzo, and Alberto Colorni. Ant system: optimization by a colony of cooperating
agents. IEEE Transactions on Systems, Man, and Cybernetics, Part B, 26(1):29–41, 1996.
[14] Marco Dorigo and Thomas Stützle. Ant colony optimization. MIT Press, 2004.
[15] Amir Hajjam El Hassani, Abder Koukam, and Lyamine Bouhafs. A hybrid ant colony system approach for the
capacitated vehicle routing problem and the capacitated vehicle routing problem with time windows. INTECH
Open Access Publisher, 2008.
[16] LF Escudero. An inexact algorithm for the sequential ordering problem. European Journal of Operational
Research, 37(2):236–249, 1988.

28

[17] LF Escudero, Monique Guignard, and Kavindra Malik. A lagrangian relax-and-cut approach for the sequential
ordering problem with precedence relationships. Annals of Operations Research, 50(1):219–237, 1994.
[18] L.F. Escudero, Monique Guignard, and Kavindra Malik. A lagrangian relax-and-cut approach for the sequential
ordering problem with precedence relationships. Annals of Operations Research, 50(1):219–237, 1994.
[19] Ahmed Ezzat and Ashraf M. Abdelbar. A less exploitative variation of the enhanced ant colony system applied
to SOP. In Proceedings of the IEEE Congress on Evolutionary Computation, CEC 2013, Cancun, Mexico,
June 20-23, 2013, pages 1917–1924. IEEE, 2013.
[20] Ahmed Ezzat, Ashraf M. Abdelbar, and Donald C. Wunsch II. A bare-bones ant colony optimization algorithm
that performs competitively on the sequential ordering problem. Memetic Computing, 6(1):19–29, 2014.
[21] Luca Maria Gambardella and Marco Dorigo. An ant colony system hybridized with a new local search for the
sequential ordering problem. INFORMS Journal on Computing, 12(3):237–255, 2000.
[22] Luca Maria Gambardella, Roberto Montemanni, and Dennis Weyland. An enhanced ant colony system for the
sequential ordering problem. In Diethard Klatte, Hans-Jakob Lüthi, and Karl Schmedders, editors, Operations
Research Proceedings 2011, Selected Papers of the International Conference on Operations Research (OR 2011),
August 30 - September 2, 2011, Zurich, Switzerland, pages 355–360. Springer, 2011.
[23] Luca Maria Gambardella, Roberto Montemanni, and Dennis Weyland. Coupling ant colony systems with strong
local searches. European Journal of Operational Research, 220(3):831–843, 2012.
[24] Fred W Glover and Gary A Kochenberger. Handbook of metaheuristics, volume 57. Springer Science & Business
Media, 2006.
[25] Luis Gouveia and Pierre Pesneau. On extended formulations for the precedence constrained asymmetric traveling
salesman problem. Networks, 48(2):77–89, 2006.
[26] Luis Gouveia and Mario Ruthmair. Load-dependent and precedence-based models for pickup and delivery
problems. Computers & OR, 63:56–71, 2015.
[27] Francesca Guerriero and M. Mancini. A cooperative parallel rollout algorithm for the sequential ordering
problem. Parallel Computing, 29(5):663–677, 2003.
[28] Keld Helsgaun. An effective implementation of K-opt moves for the Lin-Kernighan TSP heuristic. PhD thesis,
Roskilde University. Department of Computer Science, 2006.
[29] István T Hernádvölgyi. Solving the sequential ordering problem with automatically generated lower bounds. In
Operations Research Proceedings 2003, pages 355–362. Springer, 2004.
[30] Myles Hollander, Douglas A Wolfe, and Eric Chicken. Nonparametric statistical methods. John Wiley & Sons,
2013.
[31] Lester Ingber. Simulated annealing: Practice versus theory. Mathematical and computer modelling, 18(11):29–
57, 1993.
[32] Alan R. McKendall Jr. and Jin Shang. Hybrid ant systems for the dynamic facility layout problem. Computers
& OR, 33:790–803, 2006.
[33] Jimmy Lam. An efficient simulated annealing schedule. Yale University, New Haven, CT, 1988.
[34] Gregory A Mack and John H Skillings. A friedman-type rank test for main effects in a two-factor anova. Journal
of the American Statistical Association, 75(372):947–951, 1980.
[35] Roberto Montemanni, D. H. Smith, and Luca Maria Gambardella. A heuristic manipulation technique for the
sequential ordering problem. Computers & OR, 35(12):3931–3944, 2008.
[36] Vassilis Papapanagiotou, J Jamal, Roberto Montemanni, Ghassan Shobaki, and Luca Maria Gambardella. A
comparison of two exact algorithms for the sequential ordering problem. In Systems, Process and Control
(ICSPC), 2015 IEEE Conference on, pages 73–78. IEEE, 2015.

29

[37] Martín Pedemonte, Sergio Nesmachnow, and Héctor Cancela. A survey on parallel ant colony optimization.
Appl. Soft Comput., 11(8):5181–5197, 2011.
[38] Marcus Randall and James Montgomery. Candidate set strategies for ant colony optimisation. In Ant Algorithms, Third International Workshop, ANTS 2002, Brussels, Belgium, September 12-14, 2002, Proceedings,
pages 243–249, 2002.
[39] F Romeo, Vincentelli Ak Sangiovanni, and Md Huang. An efficient general cooling schedule for simulated
annealing. PROCEEDING OF IEEE INTERNATIONAL CONFERENCE ON COMPUTER AIDED DESIGN,
1986.
[40] David J Sheskin. Handbook of parametric and nonparametric statistical procedures. crc Press, 2003.
[41] Rafal Skinderowicz. The gpu-based parallel ant colony system. J. Parallel Distrib. Comput., 98:48–60, 2016.
[42] Tiago Sousa, Zita Vale, Joao Paulo Carvalho, Tiago Pinto, and Hugo Morais. A hybrid simulated annealing
approach to handle energy resource management considering an intensive use of electric vehicles. Energy, 67:81
– 96, 2014.
[43] Guendouzi Wassila and Boukra Abdelmadjid. A hybrid intrusion detection approach using ant colony system
and simulated annealing (acs-sa). In Proceedings of the International Conference on Intelligent Information
Processing, Security and Advanced Communication, IPAC ’15, pages 55:1–55:5, New York, NY, USA, 2015.
ACM.
[44] Qi Xu, Song Chen, and Bin Li. Combining the ant system algorithm and simulated annealing for 3d/2d
fixed-outline floorplanning. Appl. Soft Comput., 40:150–160, 2016.

30

