Showing newest posts with label TSP. Show older posts
Showing newest posts with label TSP. Show older posts

Monday, February 2, 2009

Greedy randomized adaptive search procedure: Step-wise construction with local search

The Iterated Local Search algorithm introduced the combination of a biased global search algorithm followed by a local search to refine the results. This general problem solving pattern is also the basis of the Greedy Randomized Adaptive Search Procedure (GRASP).

The GRASP involves alternating a greedy adaptive solution construction procedure with a local search algorithm to refine the coarser search result. The approach was designed for combinatorial problems where the individual components of a candidate solution can be discriminated in the context of the goal (such as cost). This greedy step-wise construction involves the preparation of Restricted Candidate List (RCL) of components from which to select from, and the random selection of a component from the RCL to add to the growing solution. This process continues until a viable candidate solution is constructed and can be passed onto a local search procedure for evaluation.

The power of GRASP is in the preparation of the RCL. Typically the value of individual components to a final solution is used to restrict entry to the list (a threshold). Alternative approaches involve a fixed RCL size (of the top components), and an adaptive threshold or size over time. Management of the RCL requires careful attention, ensuring sufficient flexibility in the selection of components for a solution, whilst biasing the search towards a global optimum supplemented with a local search technique. It is interesting to note that the use of a strict RCL that selects the best component during each step of construction will result in a greedy search.

Greedy randomized adaptive search procedure
This guide provide a demonstration of the GRASP applied to a standard instance of the Traveling Salesman Problem. This classical combinatorial optimization provides a good example for the step-wise construction technique used by the GRASP.

The Berlin52TSP class provides an implementation of a small sized TSP with 52 cities. The problem data is stored as constants within the class. The COORDINATES constant provides a 2D array of city coordinates and the OPTIMAL_TOUR constant stores a representation of the optimal tour (permutation) of cities that used in the algorithms termination criterion. The initialize() constructor in turn calls the build_distance_matrix method to pre-calculate a matrix of the Euclidean distance between all combinations of cities. The city-to-city distance calculation represents the most expensive aspect of candidate solution evaluations which is eased via the pre-calculation requiring only a lookup of city ID's to determine the intervening distance. The evaluate(permutation) makes use of this pre-calculation via the dist(c1, c2) method, evaluating candidate solutions in the form a 51-city permutation where each city ID is specified between 1 and 52 inclusively.
class Berlin52TSP
OPTIMAL_TOUR = [1,49,32,45,19,41,8,9,10,43,33,51,11,52,14,13,47,26,
27,28,12,25,4,6,15,5,24,48,38,37,40,39,36,35,34,44,46,16,29,50,20,
23,30,2,7,42,21,17,3,18,31,22]

COORDINATES = [[565, 575],[25, 185],[345, 750],[945, 685],[845, 655],
[880, 660],[25, 230],[525, 1000],[580, 1175],[650, 1130],[1605, 620],
[1220, 580],[1465, 200],[1530, 5],[845, 680],[725, 370],[145, 665],
[415, 635],[510, 875], [560, 365],[300, 465],[520, 585],[480, 415],
[835, 625],[975, 580],[1215, 245],[1320, 315],[1250, 400],[660, 180],
[410, 250],[420, 555],[575, 665],[1150, 1160],[700, 580],[685, 595],
[685, 610],[770, 610],[795, 645],[720, 635],[760, 650],[475, 960],
[95, 260],[875, 920],[700, 500],[555, 815],[830, 485],[1170, 65],
[830, 610],[605, 625],[595, 360],[1340, 725],[1740, 245]]

attr_reader :num_cities

def initialize()
@num_cities = COORDINATES.length
@distance_matrix = Array.new(@num_cities) {Array.new(@num_cities)}
build_distance_matrix
@optimal_tour_length = evaluate(OPTIMAL_TOUR) # calculate
end

def build_distance_matrix
# symmetrical matrix along the diag
puts "pre-calculating the TSP distance matrix"
@distance_matrix.each_with_index do |row, i|
row.each_index do |j|
row[j] = euc_2d(COORDINATES[i], COORDINATES[j])
end
end
end

def evaluate(permutation)
dist = 0
permutation.each_with_index do |c1, i|
c2 = (i==@num_cities-1) ? permutation[0] : permutation[i+1]
dist += dist(c1,c2)
end
return dist
end

# expects city numbers [1,52]
def dist(c1, c2)
@distance_matrix[c1-1][c2-1]
end

def euc_2d(c1, c2)
# As defined in TSPLIB'95 (EUC_2D)
Math::sqrt((c1[0] - c2[0])**2 + (c1[1] - c2[1])**2).round
end

def is_optimal?(scoring)
scoring == optimal_score
end

def optimal_score
@optimal_tour_length
end

# true if s1 is better score than s2
def is_better?(s1, s2)
s1 < s2 # minimizing
end
end
The Solution class provides a generic container for candidate solutions. Problem specific solution data (in this case permutations) are stored in the immutable @data instance variable, and a mutable solution costing is stored in the @score instance variable.
class Solution
attr_reader :data
attr_accessor :score

def initialize(data)
@data = data
@score = 0.0/0.0 # NaN
end

def to_s
"[#{@data.inspect}] (#{@score})"
end
end
The GRASP class provides an implementation of the Greedy Randomized Adaptive Search Procedure. The initialize(max_iterations) constructor for the algorithm sets up the maximum number of iterations for the algorithms termination criteria, as well as the algorithm specific @max_rcl_size that defines the maximum size of the Restricted Candidate List each step during construction, defaulted to 5 components. The search(problem) method provides the entry point to the search that executes the main loop until the termination criterion should_stop?(curr_it, problem) returns true. The main loop involves first a call to greedy_randomized_solution(current, problem) to construct a candidate solution, followed by a call to local_search_solution(candidate, problem) to refine the result.

The greedy_randomized_solution(solution, problem) provides a simplistic realization of the GRASP with a fixed size RCL. The construction starts with a random city and proceeds by first preparing a set of viable (unused) cities to which to connect. The candidates are ordered by their distance from the 'current' city, the list is trimmed to be no longer than the @max_rcl_size instance variable, and finally a random component from the list is selected and appended to the permutation.

The local search algorithm defined in the local_search_solution(solution, problem) method involves the application of the classical TSP local search called 2-opt (the two_opt_solution(solution) method) where two cities are selected and reconnected. This procedure is repeated 15 times as an approximation of locating the local optimum for the provided candidate solution.
class GRASP
attr_accessor :max_iterations, :max_rcl_size
attr_reader :best_solution

def initialize(max_iterations)
@max_iterations = max_iterations
@max_rcl_size = 5 # default
end

# execute the algorithm
def search(problem)
# initialize the search
@best_solution = nil
current = generate_initial_solution(problem)
evaluate_candidate_solution(current, problem)
curr_it = 0
begin
# greedy randomized solution
candidate = greedy_randomized_solution(current, problem)
evaluate_candidate_solution(candidate, problem) # eval
# local search
candidate = local_search_solution(candidate, problem)
# greedy acceptance
current = candidate if problem.is_better?(candidate.score, current.score)
curr_it += 1
end until should_stop?(curr_it, problem)
return @best_solution
end

def should_stop?(curr_it, problem)
(curr_it >= @max_iterations) or problem.is_optimal?(best_solution.score)
end

def generate_initial_solution(problem)
all = Array.new(problem.num_cities) {|i| (i+1)}
permutation = Array.new(all.length) {|i| all.delete_at(rand(all.length))}
return Solution.new(permutation)
end

def greedy_randomized_solution(solution, problem)
# construct a valid perm
perm = []
perm << rand(problem.num_cities) + 1 # random starting point
last = perm[0]
while perm.length < problem.num_cities
# create restricted candidate list for components
all = Array.new(problem.num_cities) {|i| (i+1)}
# ensure we can only add unused components
rcl = all - perm
# order by distance, asc (best have smallest distance)
rcl.sort! {|i,j| problem.dist(last, i) <=> problem.dist(last, j)}
# trim to fixed size
rcl.pop until rcl.length <= @max_rcl_size
# select random
last = rcl[rand(rcl.length)]
perm << last
end
return Solution.new(perm)
end

def local_search_solution(solution, problem)
# greedy iterated 2-opt
15.times do
candidate = two_opt_solution(solution)
evaluate_candidate_solution(candidate, problem)
if problem.is_better?(candidate.score, solution.score)
solution = candidate
end
end
return solution
end

def two_opt_solution(solution)
perm = Array.new(solution.data) # copy
# select a sub-sequence
c1, c2 = rand(perm.length), rand(perm.length)
c2 = rand(perm.length) while c1 == c2
# ensure c1 is low and c2 is high
c1, c2 = c2, c1 if c2 < c1
# reverse sub-sequence
perm[c1...c2] = perm[c1...c2].reverse
return Solution.new(perm)
end

def evaluate_candidate_solution(solution, problem)
solution.score = problem.evaluate(solution.data)
# keep track of the best solution found
if @best_solution.nil? or
problem.is_better?(solution.score, @best_solution.score)
@best_solution = solution
puts "> new best: #{solution.score}"
end
end
end
The algorithm can be executed by first creating an instance of the problem, the algorithm, and passing the problem to the algorithms search method. The algorithm will displaying its progress during the search and return the best result found. The global random number generated is seeded with a fixed number to ensure a consistent sequence of random numbers are used which can prove useful during testing.
srand(1) # set the random number seed to 1
algorithm = GRASP.new(500) # limit to 500 iterations
problem = Berlin52TSP.new # create a problem
best = algorithm.search(problem) # execute the search
puts "Best Solution: #{best}" # display the best solution
The demonstrated GRASP provides much opportunity for extension. One may modify the implementation of the TSP problem instance to provide a class that can load problems form definition files provided by the TSPLIB. The size of the RCL list may be tuned to vary the greediness of solution construction. Further, the preparation of the RCL may be modified to make use of the current working solution passed into the function, perhaps using existing selected components to bias the selection of those components during construction. Finally, adaptive RCL management procedures may be adopted as well as alternative local search procedures for the TSP.

Source Code
Download: GRASP.rb

Further Reading

Papers
Sites

Have you found a bug? Know of a great reference? Got an idea for a guide? Please send an email to jasonb@InspiredAlgorithms.com.

Need Help? Want to discuss this post? Please visit the Inspired Algorithms User Group.

Monday, January 19, 2009

Tabu Search: Using explicit memory

A problem with local search techniques is that after locating a local optima, they will continue to re-sample the same optima and neighboring solutions in the search space over and over again until a stopping condition is triggered. Tabu Search is an algorithm to address this problem by adding memory to a search algorithm, preventing or penalizing the algorithm for re-sampling (cycling) the same solutions and in turn encouraging further exploration of the search space.

The core search algorithm embedded within a tabu search generates candidate solutions in the neighborhood of a current working solution. A tabu list is maintained that holds a set of recently visited neighboring solutions (or components thereof) that may not be visited again (they are taboo). Moves are held in the tabu list for a defined tabu tenure that is typically fixed, proportional to the problem size, or generated randomly within a pre-defined range.

The tabu list is referred to as the short term memory of a search algorithm. A searches intermediate term memory refer to the factors that direct the process to promising areas (expected payoff) of the search space inferred from past experience, referred to as intensification. Examples of such intermediate memory processes include the inclusion of common components from recently visited good quality solutions, adapting the tabu tenures when locally optimal solutions are located, and restarting the search (clearing the tabu list) with a locally optimal solution. Long term memory processes are used to counter the intensification of intermediate memory by encouraging exploration in the search space, referred to as diversification. Examples include the use of penalties per component that are increased with the frequency of inclusion, and changing the acceptance criteria to allow non-improving moves.

Finally, so-called aspiration criteria may be used to permit neighbourhood moves on the tabu list (overcome the tabu). A standard aspiration criteria is to accept a tabu move if is better than a the current working position or better than the best solution found by the search to date.

Tabu Search for the Traveling Salesman Problem
This guide applies the tabu search to a classical combinatorial optimization problem called the Traveling Salesman Problem (TSP). The general class of problem involves finding an order to visit a set of cities to minimize the overall distance traveled.

The Berlin52TSP class provides a definition of a TSP instance called Berlin 52 (with 52 nodes or cities). The class is large given that it has the city coordinates hard coded in the COORDINATES class constant, as well as the optimal permutation in the OPTIMAL_TOUR constant. The evaluate(permutation) expects a valid permutation of 52 non-repeating cities to visit as an array of integers (between 1 and 52 inclusively) and is evaluated as the rounded Euclidean distance of the cycle.
class Berlin52TSP
OPTIMAL_TOUR = [1,49,32,45,19,41,8,9,10,43,33,51,11,52,14,13,47,26,
27,28,12,25,4,6,15,5,24,48,38,37,40,39,36,35,34,44,46,16,29,50,20,
23,30,2,7,42,21,17,3,18,31,22]

COORDINATES = [[565, 575],[25, 185],[345, 750],[945, 685],[845, 655],
[880, 660],[25, 230],[525, 1000],[580, 1175],[650, 1130],[1605, 620],
[1220, 580],[1465, 200],[1530, 5],[845, 680],[725, 370],[145, 665],
[415, 635],[510, 875], [560, 365],[300, 465],[520, 585],[480, 415],
[835, 625],[975, 580],[1215, 245],[1320, 315],[1250, 400],[660, 180],
[410, 250],[420, 555],[575, 665],[1150, 1160],[700, 580],[685, 595],
[685, 610],[770, 610],[795, 645],[720, 635],[760, 650],[475, 960],
[95, 260],[875, 920],[700, 500],[555, 815],[830, 485],[1170, 65],
[830, 610],[605, 625],[595, 360],[1340, 725],[1740, 245]]

attr_reader :num_cities

def initialize()
@num_cities = COORDINATES.length
@optimal_tour_length = evaluate(OPTIMAL_TOUR) # calculate
end

def evaluate(permutation)
dist = 0
permutation.each_with_index do |c1, i|
c2 = (i==@num_cities-1) ? permutation[0] : permutation[i+1]
dist += euc_2d(COORDINATES[c1-1], COORDINATES[c2-1])
end
return dist
end

def euc_2d(c1, c2)
# As defined in TSPLIB'95 (EUC_2D)
Math::sqrt((c1[0] - c2[0])**2 + (c1[1] - c2[1])**2).round
end

def is_optimal?(scoring)
scoring == optimal_score
end

def optimal_score
@optimal_tour_length
end

# true if s1 is better score than s2
def is_better?(s1, s2)
s1 < s2 # minimizing
end
end
The Solution class provides a container for candidate solutions, with an immutable @data instance variable for problem specific solution representation and a mutable @score instance variable to hold its costing. The overridden to_s method provides a pretty string representation of the solution for display on the command line. Important to this example is the overridden ==(other) method that is used to test equality between solution objects. Here, the method is specialized assuming the @data holds an array, and compares two Solution objects based on the equality of the permutations they hold. Permutation comparison is not dependent on the order the permutation is recorded, so equality is checked in both directly and with one of the permutations reversed. This equality method is used to assess whether a solution has already been added to the tabu list (later on). Finally, a mutable @count instance variable is maintained to use in the frequency a tabu solution is desired (later).
class Solution
attr_reader :data
attr_accessor :score, :count

def initialize(data)
@data = data
@score = 0.0/0.0 # NaN
@count = 0
end

def ==(other)
# permutation is direction independant
@data == other.data or @data.reverse == other.data
end

def to_s
"[#{@data.inspect}] (#{@score})"
end
end
The TabuSearchAlgorithm class offers an implementation of a iterative 2-opt search algorithm with a tabu list and aspiration criteria for a given TSP instance. The initialize(max_iterations) constructor records the maximum number of iterations the algorithm will execute and sets additional tabu search specific configuration. The @tabu_tenure instance variable defines the the size of the tabu list under the assumption that one candidate is added to the list each iteration. The @aspiration_frequency instance variable defines the maximum number of requests that may be made for a tabu solution on the tabu list before it may be reused (prematurely liberated from the tabu list).

The search(problem) provides the entry point for the search for a given TSP instance. A random starting point is used to initialize the search provided by the generate_initial_solution(problem) method. The algorithm loop generally involves generating neighbouring candidate solutions one at a time, evaluating a solution, and using it as the current working position greedily if it has an improved problem specific scoring. Neighboring candidate solutions are generated by calling the two_opt_solution(solution) method for a given current working point. This standard local search for the TSP was chosen for its speed, effectiveness, and simplicity and involves selecting two break points in a permutation and reversing the sub-sequence between the points.

The generation of neighboring candidate solutions loops until a candidate solution is generated that is not already on the tabu list. Importantly, this loop increments the number of times each tabu solution is requested, liberating a specific candidate solution if it's request frequency exceeds @aspiration_frequency, providing a basic aspiration criteria. The tabu list is maintained at the bottom of the main algorithm loop, first adding the current candidate solution regardless of whether it is accepted or not, and then whittling down the tabu_list to @tabu_tenure when required (treating the tabu_list array like a first-in-last-out FILO queue).
class TabuSearchAlgorithm
attr_accessor :max_iterations
attr_reader :best_solution

def initialize(max_iterations)
@max_iterations = max_iterations
@tabu_tenure = 100 # maximum number of solutions on the list
@aspiration_frequency = 5 # max tries to visit a taboo before allowing
end

# execute a search on the provided problem
def search(problem)
tabu_list = [] # prep the tabu list
# random starting point
current = generate_initial_solution(problem)
tabu_list << current # make candidate taboo
evaluate_candidate(current, problem)
curr_it = 0
begin
# generate new neighbour that is not tabo
selection_ok = true
begin
selection_ok = true # optimism
candidate = two_opt_solution(current) # generate
# check tabu
if tabu_list.include?(candidate)
other = tabu_list.find {|t| t==candidate}
if (other.count+=1) >= @aspiration_frequency
puts "> we have an aspiration!!!"
tabu_list.delete_if {|t| t==candidate} # ensure no duplicates
else
puts "> rejected tabu"
selection_ok = false
end
end
end until selection_ok
evaluate_candidate(candidate, problem) # evaluate
# manage short term memory
tabu_list << candidate # make candidate taboo
tabu_list.delete_at(0) while tabu_list.length > @tabu_tenure
# greedy acceptance
current = candidate if problem.is_better?(candidate.score, current.score)
curr_it += 1
end until should_stop?(curr_it, problem)
return @best_solution
end

def should_stop?(curr_it, problem)
(curr_it >= @max_iterations) or problem.is_optimal?(best_solution.score)
end

def generate_initial_solution(problem)
all = Array.new(problem.num_cities) {|i| (i+1)}
permutation = Array.new(all.length) {|i| all.delete_at(rand(all.length))}
return Solution.new(permutation)
end

def two_opt_solution(solution)
perm = Array.new(solution.data) # copy
# select a sub-sequence
c1, c2 = rand(perm.length), rand(perm.length)
c2 = rand(perm.length) while c1 == c2
# ensure c1 is low and c2 is high
c1, c2 = c2, c1 if c2 < c1
# reverse sub-sequence
perm[c1...c2] = perm[c1...c2].reverse
return Solution.new(perm)
end

def evaluate_candidate(solution, problem)
solution.score = problem.evaluate(solution.data)
# keep track of the best solution found
if @best_solution.nil? or
problem.is_better?(solution.score, @best_solution.score)
@best_solution = solution
puts " > new best: #{solution.score}"
end
end
end
The algorithm is demonstrated by first initializing the global random number generator to 1 to ensure the same sequence of random numbers on each execution (useful for testing). The problem is created and algorithm is initialized with a maximum of 2000 iterations. The search is executed by passing the problem instance to the algorithms search function, and after the execution of the run the best solution is displayed.
srand(1) # set the random number seed to 1
algorithm = TabuSearchAlgorithm.new(2000) # limit to 2000 iterations
problem = Berlin52TSP.new # create a problem
best = algorithm.search(problem) # execute the search
puts "Best Solution: #{best}" # display the best solution
During the execution of the algorithm, the scoring of the best solution is displayed each time a new best is located. The loop of the main algorithm displays a message each time a tabu solutions aspiration criteria is triggered or when a generated candidate solution is denied based on its existence in the tabu list. You will observe that as the search gets closer to the global optimum, the number of candidates denied increases and a number of liberation events occur as the aspiration criteria are triggered.

There is plenty of room for extension on this example. One may change the configuration of the tabu search, increasing or decreasing both the tabu tenure as well as the aspiration frequency. Also experimentation of the search algorithm with and without a tabu list may be interesting to see if and when the short term memory of this tabu search example are useful. Alternative representations may be used for the tabu list such as the selection of cities to swap in the local search. Alternatively, the local search may be adjusted to construct tours in a step-wise (per-city) manner and maintain tabu information per-edge in an adjacency matrix. This matrix would then discourage or prevent the use of specific edges during solution construction. Finally, additional memory processes may be introduced such as search-restarts with the best found solution for intermediate memory, and the solution acceptance criteria may be made less greedy to allow non-improving moves to be taken.

Source Code
Download: TabuSearch.rb

Further Reading


Papers:
Sites:

Have you found a bug? Know of a great reference? Got an idea for a guide? Please send an email to jasonb@InspiredAlgorithms.com.

Need Help? Want to discuss this post? Please visit the Inspired Algorithms User Group.

Tuesday, January 13, 2009

Iterated Local Search: Cycling global and local search

A Hill Climbing algorithm such as the one demonstrated in Localized Random Search is an example of a Local Search Algorithm. A local search involves moving through a search space from one neighboring candidate solution to the next, such as from vertex to vertex along the edges of a graph. The Iterated Local Search algorithm involves the repeated application of a local search algorithm applied to the candidate solutions found by a broader search process that involves a biased random walk through the search space.

The algorithm works by first selecting a starting point for the search either randomly or via a domain specific construction heuristic. The starting position is then optimized to an approximation of the local optimum using a local search strategy. The algorithm loop involves three steps: a perturbation of the current position, the application of the local search to the perturbation, and an acceptance decision of whether or not the locally optimizing candidate solution should replace the current working position for the search.

A good starting point is important for shorter algorithm runs, whereas this importance degrades as the length of the run is increased because the algorithm is given more opportunity to improve. The Iterative Local Search algorithm resembles random restart of a local search technique, although out-performs this technique given the dependence between each new starting position of the search. The perturbation of the current working position provides long jumps in the search space that are further refined by the local search. If the jumps are too large, performance degrades as the behaviors resembles random restart. If the jumps are too small in the search space, the broader search technique resembles the local search algorithm.

The goal is to find a balance between the perturbation and local search mechanisms within the algorithm such that the local search technique cannot easily undo the longer jumps made by the perturbations. There is a relationship between the perturbations and the acceptance criteria for replacing the current working position. Always accepting improved points may result in a greedy algorithm that only accepts small steps and cannot explore distant regions of the search space. Longer jumps made by the perturbation mechanism may require a more flexible (less greedy) acceptance criteria for the search to progress.

Iterated Local Search of the Traveling Salesman Problem
The Iterated Local Search algorithm is suited to problem domains where the good solutions in the search space are clustered. The approach is more commonly applied to combinatorial optimization problems like the Traveling Salesman Problem (TSP) where this clustering of good solutions is manifest as the sharing of common edges. This guide demonstrates the application of the Iterated Local Search algorithm applied to a standardized instance of the TSP from the TSP library (TSPLIB).

The Berlin52TSP class contains a the city coordinate data and optimal tour data as COORDINATES and OPTIMAL_TOUR constants respectively, taken from the TSPLIB definition files for the Berlin52 problem instance. The initialize() constructor sets the accessible @num_cities instance variable and calculates the @optimal_tour_length used by the is_optimal?(scoring) method for assessing scores to see if an optimal tour has been discovered. The evaluate(permutation)method calculates the distance of a given tour permutation which is an array of non-repeating integers that specify cities between 1 and 52 inclusively. Evaluation involves calculating the Euclidean distance for each of the city pairs using the problem definition coordinates. The euc_2d(c1, c2) method matches the Euclidean calculation specified by the TSPLIB documentation, most notably involving the rounding of calculated distances to integers for this symmetrical TSP instance.
class Berlin52TSP
OPTIMAL_TOUR = [1,49,32,45,19,41,8,9,10,43,33,51,11,52,14,13,47,26,
27,28,12,25,4,6,15,5,24,48,38,37,40,39,36,35,34,44,46,16,29,50,20,
23,30,2,7,42,21,17,3,18,31,22]

COORDINATES = [[565, 575],[25, 185],[345, 750],[945, 685],[845, 655],
[880, 660],[25, 230],[525, 1000],[580, 1175],[650, 1130],[1605, 620],
[1220, 580],[1465, 200],[1530, 5],[845, 680],[725, 370],[145, 665],
[415, 635],[510, 875], [560, 365],[300, 465],[520, 585],[480, 415],
[835, 625],[975, 580],[1215, 245],[1320, 315],[1250, 400],[660, 180],
[410, 250],[420, 555],[575, 665],[1150, 1160],[700, 580],[685, 595],
[685, 610],[770, 610],[795, 645],[720, 635],[760, 650],[475, 960],
[95, 260],[875, 920],[700, 500],[555, 815],[830, 485],[1170, 65],
[830, 610],[605, 625],[595, 360],[1340, 725],[1740, 245]]

attr_reader :num_cities

def initialize()
@num_cities = COORDINATES.length
@optimal_tour_length = evaluate(OPTIMAL_TOUR) # calculate
end

def evaluate(permutation)
dist = 0
permutation.each_with_index do |c1, i|
c2 = (i==@num_cities-1) ? permutation[0] : permutation[i+1]
dist += euc_2d(COORDINATES[c1-1], COORDINATES[c2-1])
end
return dist
end

def euc_2d(c1, c2)
# As defined in TSPLIB'95 (EUC_2D)
Math::sqrt((c1[0] - c2[0])**2 + (c1[1] - c2[1])**2).round
end

def is_optimal?(scoring)
scoring == optimal_score
end

def optimal_score
@optimal_tour_length
end

# true if s1 is better score than s2
def is_better?(s1, s2)
s1 < s2 # minimizing
end
end
The Solution class is a classical container for candidate solutions, holding a TSP permutation in the immutable @data instance variable and managing a mutable permutation scoring in the @score instance variable.
class Solution
attr_reader :data
attr_accessor :score

def initialize(data)
@data = data
@score = 0.0/0.0 # NaN
end

def to_s
"[#{@data.inspect}] (#{@score})"
end
end
The IteratedLocalAlgorithm class is quite large. The search(problem) method is generalized and executes the iterated local search algorithm on the provided problem definition. It starts by calling generate_initial_solution(problem) to prepare a starting point for the search and refining it with a call to local_search_solution(current, problem). The algorithms main loop is a text book implementation involving repeated rounds of generating a perturbed versions of the current working solution, refining it with a local search procedure, and in this case using a greedy (candidate must be better) acceptance criteria.

The meat of the algorithm are specific to TSP permutations. The generate_initial_solution(problem) generates a random permutation as the starting point of the search. The perturb_solution(solution) method generates a 4-opt variation of a given permutation, also known as a double-bridge move. Basically the permutation is split into 4 pieces and reordered. This function may be called repeatedly if more drastic perturbations are desired, perhaps on large problem instances. The local_search_solution(solution, problem) method is an iterative application of the 2-opt procedure, that evaluates generated solutions and greedily maintains the best candidate solution discovered. The 2-opt procedure involves breaking the tour permutation into two parts (two break points) and reconnecting the tour back together by switching the end points (reversing the sequence between the break points). This is a classical and fast local search procedure for TSP and here it is repeated 30 times.
class IteratedLocalAlgorithm
attr_accessor :max_iterations
attr_reader :best_solution

def initialize(max_iterations)
@max_iterations = max_iterations
end

# execute a iterated local search on the provided problem
def search(problem)
# random starting point
current = generate_initial_solution(problem)
@best_solution = current
evaluate_candidate_solution(current, problem)
# local search
local_best = local_search_solution(current, problem)
current = local_best if problem.is_better?(local_best.score, current.score)
curr_it = 0
begin
# generate perturbation
pert_solution = perturb_solution(current)
evaluate_candidate_solution(pert_solution, problem)
# local search
local_best = local_search_solution(pert_solution, problem)
# greedy acceptance
current = local_best if problem.is_better?(local_best.score, current.score)
curr_it += 1
end until should_stop?(curr_it, problem)
return @best_solution
end

def should_stop?(curr_it, problem)
(curr_it >= @max_iterations) or problem.is_optimal?(best_solution.score)
end

def generate_initial_solution(problem)
all = Array.new(problem.num_cities) {|i| (i+1)}
permutation = Array.new(all.length) {|i| all.delete_at(rand(all.length))}
return Solution.new(permutation)
end

def perturb_solution(solution)
data = solution.data
length = data.length
# double-bridge move (4-opt), break into 4 parts (a,b,c,d)
pos1 = 1 + rand(length / 4)
pos2 = pos1 + 1 + rand(length / 4)
pos3 = pos2 + 1 + rand(length / 4)
# put it back together (a,d,c,b)
perm = data[0...pos1] + data[pos3...length] +
data[pos2...pos3] + data[pos1...pos2]
return Solution.new(perm)
end

def local_search_solution(solution, problem)
# greedy iterated 2-opt
30.times do
candidate = two_opt_solution(solution)
evaluate_candidate_solution(candidate, problem)
if problem.is_better?(candidate.score, solution.score)
solution = candidate
end
end
return solution
end

def two_opt_solution(solution)
perm = Array.new(solution.data) # copy
# select a sub-sequence
c1, c2 = rand(perm.length), rand(perm.length)
c2 = rand(perm.length) while c1 == c2
# ensure c1 is low and c2 is high
c1, c2 = c2, c1 if c2 < c1
# reverse sub-sequence
perm[c1...c2] = perm[c1...c2].reverse
return Solution.new(perm)
end

def evaluate_candidate_solution(solution, problem)
solution.score = problem.evaluate(solution.data)
# keep track of the best solution found
if problem.is_better?(solution.score, @best_solution.score)
@best_solution = solution
puts " > new best: #{solution.score}"
end
end
end
The demonstration of the algorithm involves first seeding the global random number generator to one, and constructing both an instance of the problem and the algorithm. The algorithm is executed on the problem instance and the algorithm stop condition is triggered if the optimal solution is discovered or a maximum of 2000 iterations of the main algorithm loop are completed.
srand(1) # set the random number seed to 1
algorithm = IteratedLocalAlgorithm.new(1000) # limit to 1000 iterations
problem = Berlin52TSP.new # create a problem
best = algorithm.search(problem) # execute the search
puts "Best Solution: #{best}" # display the best solution
There is plenty of room for extension of this interesting guide. The problem definition may be made generic such that it loads problem instance descriptions and optimal tours from TSPLIB files. The problem class may be made computationally more efficient by pre-calculating the city distance matrix on construction and looking up these distance evaluations in the evaluate method. The use of a randomly generated starting point for the search suggests that the algorithm may need to be executed for a long time to approximate the global optima (longer than 2000 iterations). The algorithm may be adjusted to use a deterministically generated nearest-neighbor solution as the starting point of the search that expected to be a much better quality solution. Finally, more advanced local search procedures exist for the TSP (such as Lin-Kernighan), which may be integrated into the algorithm.

Source Code
Download: IteratedLocalSearch.rb

Further Reading

Papers
Sites

Have you found a bug? Know of a great reference? Got an idea for a guide? Please send an email to jasonb@InspiredAlgorithms.com

Need Help? Want to discuss this post? Please visit the Inspired Algorithms User Group.