Showing newest posts with label traveling salesman problem. Show older posts
Showing newest posts with label traveling salesman problem. 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.