Greedy Randomized Adaptive Search 算法超详细解析，附代码实现TSP问题求解
Greedy Randomized Adaptive Search 算法超详细解析，附代码实现TSP问题求解
这两天刚好看到这个算法，然后就写一写吧。贪心随机自适应搜索虽然算是一个比较简单的启发式，但是效果也非常不错的。
01 概述
Greedy Randomized Adaptive Search，贪婪随机自适应搜索（GRAS），是组合优化问题中的多起点元启发式算法。
在算法的每次迭代中，主要由两个阶段组成：构造（construction）和局部搜索（ local search）。
构造（construction）阶段主要用于生成一个可行解，而后该初始可行解会被放进局部搜索进行邻域搜索，直到找到一个局部最优解为止。
02 整体框架
如上面所说，其实整一个算法的框架相对于其他算法来说还算比较简单明了，大家可以先看以下整体的伪代码：
(GRAS伪代码)
GRAS主要由两部分组成：
1) Greedy_Randomized_Construction：在贪心的基础上，加入一定的随机因素，构造初始可行解。
2) Local Search：对上面构造的初始可行解进行邻域搜索，直到找到一个局部最优解。
然后再多说两句：
Repair是什么鬼？有时候由于随机因素的加入，Greedy_Randomized_Construction阶段生成的解不一定都是可行解，所以为了保证下一步的Local Search能继续进行，加入repair算子，对解进行修复，保证其可行。
不是说自适应（Adaptive）吗？我怎么没看到Adaptive 的过程？ 别急，这个后面具体举例的时候会详细讲到。
03 举个例子说明
为了大家能更加深入理解该算法，我们举一个简单的例子来为大家详细讲解算法的流程。
假如有一个有限的集合E = {1, 2, 3, ……, n}。可行解的集合
，（我知道你们想问
是什么，由E的所有子集作为元素构成的集合。）目标函数
，求目标函数的最小值以及对应的解。
好了，相信大家都看懂上面的问题了（看不懂也别问我，摊手）。对于上述问题，我们来一步一个脚印用GRAS来求解之，来，跟紧小编的脚步……
强调了很多次，GRAS由两部分组成：Greedy_Randomized_Construction和Local Search，所以，在求解具体问题的时候，完成这两部分的设计，然后按照第二节所示的框架搭起来就可以。
3.1 Greedy_Randomized_Construction
这里还是老规矩，先上伪代码给大家看看，然后我们再进行讲解，毕竟对于算法来说，伪代码的作用不言而喻。
Greedy_Randomized_Construction
第1行，一开始解是一个空集。
第2行，初始化候选元素的集合，这里候选元素是指能放进Solution的元素（也就是目前Solution里面没有的），比如１，２，３……。
第3行，对候选集合的每个元素进行评估，计算将元素x放入Solution会导致目标函数f改变量delta_x。
第5行，根据delta_x对各个元素排序，选取部分较好的候选元素组成RCL表（贪心性体现在这里）。
第6行，随机在RCL中选取一个元素放进Solution。（算法的随机性）
第8、9行，更新候选元素集合，然后对每个元素进行重新评估计算delta值。（算法的自适应性体现在这里）
相信经过上面如此详细的介绍，大家都懂了吧！
3.2 Local Search
关于Local Search方面的内容，相信大家学习heuristic这么久了，就不用我多说什么了吧：
(Local Search)
简单看一下伪代码即可，主要是邻域算子的设计，然后就是在邻域里面进行搜索，找到一个局部最优解为止。
然后关于邻域搜索，有best-improving or first-improving strategy 两种策略，这个下次有时间出个专题给大家讲明白一些相关概念吧。
04 再论Greedy_Randomized_Construction
前面我们说了，Greedy_Randomized_Construction用于生成初始解，既然是Greedy_Randomized两个结合体，那么肯定就有一个权重分配的问题。
即，是Greedy成分多一点呢？还是Randomized成分多一点好呢？因此，为了控制这两个小老弟的权重，防止某个家伙在该过程中用力过猛导致解不那么好的情况，我们引入一个参数α：
其他部分就不再多说，可以看到，上面的α参数主要是控制RCL（什么是RCL回头去看！）的长度：
当α=0时，纯贪心，只能选取最优的候选元素。
当α=1时，纯随机，所有候选元素都可随机选。
05 代码实现
由于小编精力有限，就不从头写一遍了，从GitHub上找了一个感觉还不错的算法给大家，也是求解TSP问题的。
不过说实在的，python写算法的速度是很慢的，无论是速度还是算法架构等方面都不推荐大家用matlab或者python这种脚本性的语言写大型优化算法。
运行结果如下：
（Berlin52）
代码算例以及相关运行结果请移步留言区。
############################################################################
# Created by: Prof. Valdecy Pereira, D.Sc.
# UFF - Universidade Federal Fluminense (Brazil)
# email: valdecy.pereira@gmail.com
# Course: Metaheuristics
# Lesson: Local Search-GRASP
# Citation:
# PEREIRA, V. (2018). Project: Metaheuristic-Local_Search-GRASP, File: Python-MH-Local Search-GRASP.py, GitHub repository: <https://github.com/Valdecy/Metaheuristic-Local_Search-GRASP>
############################################################################
# Required Libraries
import pandas as pd
import random
import numpy as np
import copy
import os
from matplotlib import pyplot as plt
# Function: Tour Distance
def distance_calc(Xdata, city_tour):
distance = 0
for k in range(0, len(city_tour[0])-1):
m = k + 1
distance = distance + Xdata.iloc[city_tour[0][k]-1, city_tour[0][m]-1]
return distance
# Function: Euclidean Distance
def euclidean_distance(x, y):
distance = 0
for j in range(0, len(x)):
distance = (x.iloc[j] - y.iloc[j])**2 + distance
return distance**(1/2)
# Function: Initial Seed
def seed_function(Xdata):
seed = [[],float("inf")]
sequence = random.sample(list(range(1,Xdata.shape[0]+1)), Xdata.shape[0])
sequence.append(sequence[0])
seed[0] = sequence
seed[1] = distance_calc(Xdata, seed)
return seed
# Function: Build Distance Matrix
def buid_distance_matrix(coordinates):
Xdata = pd.DataFrame(np.zeros((coordinates.shape[0], coordinates.shape[0])))
for i in range(0, Xdata.shape[0]):
for j in range(0, Xdata.shape[1]):
if (i != j):
x = coordinates.iloc[i,:]
y = coordinates.iloc[j,:]
Xdata.iloc[i,j] = euclidean_distance(x, y)
return Xdata
# Function: Tour Plot
def plot_tour_distance_matrix (Xdata, city_tour):
m = Xdata.copy(deep = True)
for i in range(0, Xdata.shape[0]):
for j in range(0, Xdata.shape[1]):
m.iloc[i,j] = (1/2)*(Xdata.iloc[0,j]**2 + Xdata.iloc[i,0]**2 - Xdata.iloc[i,j]**2)
m = m.values
w, u = np.linalg.eig(np.matmul(m.T, m))
s = (np.diag(np.sort(w)[::-1]))**(1/2)
coordinates = np.matmul(u, s**(1/2))
coordinates = coordinates.real[:,0:2]
xy = pd.DataFrame(np.zeros((len(city_tour[0]), 2)))
for i in range(0, len(city_tour[0])):
if (i < len(city_tour[0])):
xy.iloc[i, 0] = coordinates[city_tour[0][i]-1, 0]
xy.iloc[i, 1] = coordinates[city_tour[0][i]-1, 1]
else:
xy.iloc[i, 0] = coordinates[city_tour[0][0]-1, 0]
xy.iloc[i, 1] = coordinates[city_tour[0][0]-1, 1]
plt.plot(xy.iloc[:,0], xy.iloc[:,1], marker = 's', alpha = 1, markersize = 7, color = 'black')
plt.plot(xy.iloc[0,0], xy.iloc[0,1], marker = 's', alpha = 1, markersize = 7, color = 'red')
plt.plot(xy.iloc[1,0], xy.iloc[1,1], marker = 's', alpha = 1, markersize = 7, color = 'orange')
return
# Function: Tour Plot
def plot_tour_coordinates (coordinates, city_tour):
coordinates = coordinates.values
xy = pd.DataFrame(np.zeros((len(city_tour[0]), 2)))
for i in range(0, len(city_tour[0])):
if (i < len(city_tour[0])):
xy.iloc[i, 0] = coordinates[city_tour[0][i]-1, 0]
xy.iloc[i, 1] = coordinates[city_tour[0][i]-1, 1]
else:
xy.iloc[i, 0] = coordinates[city_tour[0][0]-1, 0]
xy.iloc[i, 1] = coordinates[city_tour[0][0]-1, 1]
plt.plot(xy.iloc[:,0], xy.iloc[:,1], marker = 's', alpha = 1, markersize = 7, color = 'black')
plt.plot(xy.iloc[0,0], xy.iloc[0,1], marker = 's', alpha = 1, markersize = 7, color = 'red')
plt.plot(xy.iloc[1,0], xy.iloc[1,1], marker = 's', alpha = 1, markersize = 7, color = 'orange')
return
# Function: Rank Cities by Distance
def ranking(Xdata, city = 0):
rank = pd.DataFrame(np.zeros((Xdata.shape[0], 2)), columns = ['Distance', 'City'])
for i in range(0, rank.shape[0]):
rank.iloc[i,0] = Xdata.iloc[i,city]
rank.iloc[i,1] = i + 1
rank = rank.sort_values(by = 'Distance')
return rank
# Function: RCL
def restricted_candidate_list(Xdata, greediness_value = 0.5):
seed = [[],float("inf")]
sequence = []
sequence.append(random.sample(list(range(1,Xdata.shape[0]+1)), 1)[0])
for i in range(0, Xdata.shape[0]):
count = 1
rand = int.from_bytes(os.urandom(8), byteorder = "big") / ((1 << 64) - 1)
if (rand > greediness_value and len(sequence) < Xdata.shape[0]):
next_city = int(ranking(Xdata, city = sequence[-1] - 1).iloc[count,1])
while next_city in sequence:
count = np.clip(count+1,1,Xdata.shape[0]-1)
next_city = int(ranking(Xdata, city = sequence[-1] - 1).iloc[count,1])
sequence.append(next_city)
elif (rand <= greediness_value and len(sequence) < Xdata.shape[0]):
next_city = random.sample(list(range(1,Xdata.shape[0]+1)), 1)[0]
while next_city in sequence:
next_city = int(random.sample(list(range(1,Xdata.shape[0]+1)), 1)[0])
sequence.append(next_city)
sequence.append(sequence[0])
seed[0] = sequence
seed[1] = distance_calc(Xdata, seed)
return seed
# Function: 2_opt
def local_search_2_opt(Xdata, city_tour):
tour = copy.deepcopy(city_tour)
best_route = copy.deepcopy(tour)
seed = copy.deepcopy(tour)
for i in range(0, len(tour[0]) - 2):
for j in range(i+1, len(tour[0]) - 1):
best_route[0][i:j+1] = list(reversed(best_route[0][i:j+1]))
best_route[0][-1] = best_route[0][0]
best_route[1] = distance_calc(Xdata, best_route)
if (best_route[1] < tour[1]):
tour[1] = copy.deepcopy(best_route[1])
for n in range(0, len(tour[0])):
tour[0][n] = best_route[0][n]
best_route = copy.deepcopy(seed)
return tour
# Function: GRASP
def greedy_randomized_adaptive_search_procedure(Xdata, city_tour, iterations = 50, rcl = 25, greediness_value = 0.5):
count = 0
best_solution = copy.deepcopy(city_tour)
while (count < iterations):
rcl_list = []
for i in range(0, rcl):
rcl_list.append(restricted_candidate_list(Xdata, greediness_value = greediness_value))
candidate = int(random.sample(list(range(0,rcl)), 1)[0])
city_tour = local_search_2_opt(Xdata, city_tour = rcl_list[candidate])
while (city_tour[0] != rcl_list[candidate][0]):
rcl_list[candidate] = copy.deepcopy(city_tour)
city_tour = local_search_2_opt(Xdata, city_tour = rcl_list[candidate])
if (city_tour[1] < best_solution[1]):
best_solution = copy.deepcopy(city_tour)
count = count + 1
print("Iteration =", count, "-> Distance =", best_solution[1])
print("Best Solution =", best_solution)
return best_solution
######################## Part 1 - Usage ####################################
X = pd.read_csv('Python-MH-Local Search-GRASP-Dataset-01.txt', sep = '\t') #17 cities = 1922.33
seed = seed_function(X)
lsgrasp = greedy_randomized_adaptive_search_procedure(X, city_tour = seed, iterations = 5, rcl = 5, greediness_value = 0.5)
plot_tour_distance_matrix(X, lsgrasp) # Red Point = Initial city; Orange Point = Second City # The generated coordinates (2D projection) are aproximated, depending on the data, the optimum tour may present crosses.
Y = pd.read_csv('Python-MH-Local Search-GRASP-Dataset-02.txt', sep = '\t') # Berlin 52 = 7544.37
X = buid_distance_matrix(Y)
seed = seed_function(X)
lsgrasp = greedy_randomized_adaptive_search_procedure(X, city_tour = seed, iterations = 10, rcl = 15, greediness_value = 0.5)
plot_tour_coordinates (Y, lsgrasp) # Red Point = Initial city; Orange Point = Second City