1. 基础模块
1.1 读取数据
标准的tsp格式可以用tsplib95库进行读取:
import numpy as np
import tsplib95
import matplotlib.pyplot as plt
problem = tsplib95.load('ALNS/examples/xqf131.tsp')
coords = np.array([(coord[0],coord[1]) for city, coord in problem.node_coords.items()])
solution = tsplib95.load('ALNS/examples/xqf131.opt.tour')
optimal = data.trace_tours(solution.tours)[0]
1.2 距离计算
使用numba进行加速,因此这里重新写了距离的计算方法:
@njit
def dist(a, b):
xd,yd =coords[a]-coords[b]
return round(math.sqrt(xd*xd+yd*yd))
@njit
def totaldist(tour):
s = 0
for i in np.arange(len(tour)):
s+=dist(tour[i],tour[(i+1)%len(tour)])
return s
如果内存足够大,甚至可以提前计算好距离矩阵:
from scipy.spatial.distance import cdist
distA=cdist(coords,coords,metric='euclidean')
distMatrix = np.round(distA).astype(np.int32)
@njit
def dist(a, b):
return distMatrix[a][b]
在Monalisa_100K中,距离矩阵保存下来是37.3G大小。
2. 启发式算法
2.1 贪心算法
从节点0开始,每次从未访问的点中,选出与当前节点最近的一个点。
@njit
def greedyConstruct(coords):
num = len(coords)
uv = np.arange(1,num) # unvisited
newTour = np.array([0]*(num))
for i in np.arange(1,num):
cur = newTour[i-1]
mindist = sdist(cur,uv[0])
ind = 0
for uvi,j in enumerate(uv):
dis = dist(cur,j)
if dis < mindist:
mindist = dis
ind = uvi
newTour[i] = uv[ind]
uv = np.delete(uv,ind)
return newTour
对Monalisa问题,用时20min10s(如果用内存缓存,则耗时10s),结果为6886142,画出来如下:
2.2 2-opt算法
2-opt本质是一个边交换启发式算法,非常简单。首先用myopic的方式生成一条路线,然后逐点使用2-opt算法进行优化。
首先是顺序遍历、myopic的方式:
@njit
def two_opt(tour):
num = len(tour)
for i in range(num - 1):
for j in range(i + 2, num):
change = dist(tour[i], tour[j]) + dist(tour[i+1], tour[(j+1)%num]) \
- dist(tour[i], tour[i+1]) - dist(tour[j], tour[(j+1)%num])
if change < -1:
return i,j
@njit
def simple_iter(iter_num,tour):
for i in range(iter_num):
i,j = two_opt(tour)
tour[i+1:j+1] = tour[i+1:j+1][::-1]
return tour
from tqdm import tqdm
for i in tqdm(range(1000)):
tour = random_iter(1000,tour)
print(totaldist(tour))
2.3 swap算法
swap算法是一个点交换算法,效率相对较低
@njit
def swap(tour):
num = len(tour)
for i in range(num - 3):
for j in range(i + 2, num-1):
change = dist(tour[i], tour[j+1]) + dist(tour[j+1], tour[i+2]) \
+ dist(tour[j], tour[(i+1)])+ dist(tour[i+1], tour[(j+2)%num]) \
- dist(tour[i], tour[i+1]) - dist(tour[i+1], tour[i+2]) \
- dist(tour[j], tour[(j+1)])- dist(tour[j+1], tour[(j+2)%num])
if change < -1:
return i,j
@njit
def simple_iter(iter_num,tour):
for i in range(iter_num):
i,j = swap(tour)
tour[i+1],tour[j+1] = tour[j+1],tour[i+1]
return tour
%time tour = simple_iter(100,tour)
2.4 insert算法
insert是一个点边交换的算法
@njit
def swap(tour):
num = len(tour)
for i in range(num - 3):
for j in range(i + 2, num-1):
change = dist(tour[i], tour[j+1]) + dist(tour[j+1], tour[i+2]) \
+ dist(tour[j], tour[(i+1)])+ dist(tour[i+1], tour[(j+2)%num]) \
- dist(tour[i], tour[i+1]) - dist(tour[i+1], tour[i+2]) \
- dist(tour[j], tour[(j+1)])- dist(tour[j+1], tour[(j+2)%num])
if change < -1:
return i,j
@njit
def simple_iter(iter_num,tour):
for i in range(iter_num):
i,j = swap(tour)
tour[i+1],tour[j+1] = tour[j+1],tour[i+1]
return tour
%time tour = simple_iter(100,tour)
2.5 post Process
相邻的k个点遍历所有可能性,找到最优解。首先是CPU版本,使用cache缓存加快计算速度
# post-process
import itertools
from functools import lru_cache as cache
k = 9
num = len(tour)
@cache(None)
def dist(a, b):
return round(np.linalg.norm(coords[a]-coords[b]))
def optimal_subpath(startindex,tour):
def cost(path):
s = 0
for i in np.arange(len(path)-1):
s+=dist(path[i],path[i+1])
return s+dist(path[0],tour[startindex-1])+dist(path[-1],tour[startindex+k])
nodes = [tour[i] for i in range(startindex,startindex+k)]
subpath = min(itertools.permutations(nodes, k), key=cost)
return cost(subpath)-cost(nodes),subpath
def getRes(s,tour):
for m in tqdm(range(s,num-k-1)):
r,p = optimal_subpath(m,tour)
if r<-1:
return m,p,r
return -1,-1,-1
s = 1
while s!=-1:
s,path,r = getRes(s,tour)
tour[s:s+k] = path
print(i,totaldist(tour))
然后是GPU版本
@cuda.jit
def postcudakernel(tour,coords,paths,z,dz):
def dis(a,b):
x1,y1 = coords[a]
x2,y2 = coords[b]
x = x1-x2
y = y1-y2
return math.sqrt(x*x+y*y)
k = cuda.grid(1)
num = len(z)
if k < num:
path = paths[0]
minindex = 0
firstdist = 0
for j in range(len(path)-1):
firstdist += dis(tour[(k+path[j])%num],tour[(k+path[j+1])%num])
mindist = firstdist
for i in range(1,len(paths)):
path = paths[i]
curdist = 0
for j in range(len(path)-1):
curdist += dis(tour[(k+path[j])%num],tour[(k+path[j+1])%num])
if curdist > mindist:
break
if curdist<mindist:
mindist = curdist
minindex = i
z[k] = minindex
dz[k] = firstdist-mindist
2.6 GA算法
GA可以当做是多条路的insert算法-把第二条路径的某些点insert到第一条路径里面。
首先,获取N个解,随机选取2条,比如:
tour1= [1,2,3,4,5,6,7,8]
tourt2=[5,6,7,8,4,3,2,1]
再随机选择两个点,比如3、6,则新的tour3 = [x,x,x,8,4,3,x,x]
。如果起点大于终点,那么把起点减去路径的长度即可。
剩下的补齐即可,变为tour3 = [1,2,5,8,4,3,6,7]
3. 启发式算法的加速思路
3.1 随机化
将顺序搜索变为随机搜索即可,下面是2-opt的随机化版本
from numpy.random import randint
@njit
def two_opt_random(tour,iternum):
num = len(tour)
for k in range(iternum):
i = randint(num - 2)
j = randint(i + 2, num)
change = dist(tour[i], tour[j]) + dist(tour[i+1], tour[(j+1)%num]) \
- dist(tour[i], tour[i+1]) - dist(tour[j], tour[(j+1)%num])
if change < -1:
return i,j
return 0,0
@njit
def random_iter(iter_num,tour):
for k in range(iter_num):
i,j = two_opt_random(tour,10000)
tour[i+1:j+1] = tour[i+1:j+1][::-1]
return tour
%time tour = random_iter(1000,tour)
3.2 并行化
改造成并行化之后,已经有GPU版本的雏形了。
from numpy.random import randint
from numba import prange,config
config.NUMBA_DEFAULT_NUM_THREADS=2000
iter2 = 1000000
@njit(parallel=True,fastmath = True)
def two_opt_kernel(tour,x,y):
num = len(tour)
z = np.zeros(iter2).astype(np.float32)
for k in prange(iter2):
i,j = x[k],y[k]
z[k] = dist(tour[i], tour[j]) + dist(tour[(i+1)%num], tour[(j+1)%num]) \
- dist(tour[i], tour[(i+1)%num]) - dist(tour[j], tour[(j+1)%num])
return z
@njit
def two_opt_random(tour,iternum):
num = len(tour)
for k in prange(iternum):
x = np.random.randint(0,num,iter2).astype(np.int32)
y = (np.random.randint(2,num,iter2).astype(np.int32) +x)%num
z = two_opt_kernel(tour,x,y)
ind = np.argmin(z)
if z[ind] < -1:
i,j =x[ind],y[ind]
if i < j:
tour[i+1:j+1] = tour[i+1:j+1][::-1]
else:
tour[j+1:i+1] = tour[j+1:i+1][::-1]
return tour
from tqdm import tqdm
for i in tqdm(range(1000)):
tour = random_iter(1000,tour)
print(totaldist(tour))
2-opt速度非常快,Monalisa_100K使用cpu,从贪心算法的初始解开始,经过10min便下降到了6034440。
3.3 上GPU
这里有几个注意点:1. cuda函数内不能新建数组,所有的数组需要在cpu上创建好,然后传入gpu。2. cuda里面的函数不能嵌套。3. numpy函数支持及其有限,比如numpy.linaq就不支持,很多函数需要自己重写。
这里主要使用2-opt的GPU版本
from tqdm import tqdm
num = 100000
@cuda.jit
def twooptkernel(tour,coords,x,y,z):
def dis(a,b):
x1,y1 = coords[a]
x2,y2 = coords[b]
return math.sqrt(math.pow(x1-x2,2)+math.pow(y1-y2,2))
k = cuda.grid(1)
if k < len(z):
i,j = x[k],y[k]
z[k] = dis(tour[i], tour[j]) + dis(tour[(i+1)%num], tour[(j+1)%num]) \
- dis(tour[i], tour[(i+1)%num]) - dis(tour[j], tour[(j+1)%num])
tour = greedyConstruct(num)
iter1,iter2 = 1000000,10000
for inum in tqdm(range(iter1)):
if inum%1000==0:
print(totaldist(tour))
x = np.random.randint(0,num,iter2).astype(np.int32)
y = (np.random.randint(2,num,iter2).astype(np.int32) +x)%num
z = np.zeros(iter2).astype(np.float32)
twooptkernel[iter2//16, 16](tour,coords,x,y,z)
cuda.synchronize()
ind = np.argmin(z)
if z[ind] < 0:
i,j =x[ind],y[ind]
if i < j:
tour[i+1:j+1] = tour[i+1:j+1][::-1]
else:
tour[j+1:i+1] = tour[j+1:i+1][::-1]
从贪心算法开始的图如下:用时20min20s,下降趋势图如下,明显看出下降趋势开始减慢:
这时候的路线图如下,总距离为6037645。
我们再来看下github上另一个算法给出的结果,用了11.8小时得到了5896819的结果(比我们10min的结果优3%,比我们20min的结果优2.39%)。