运筹系列61：TSP问题GPU算法初探

IE06

已于 2023-02-17 08:28:48 修改

阅读量1.3k

点赞数 6

分类专栏：运筹学文章标签： python

于 2021-03-29 00:38:33 首次发布

本文链接：https://blog.csdn.net/kittyzc/article/details/115289910

版权

运筹学专栏收录该内容

93 篇文章 360 订阅

订阅专栏

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%)。

IE06

关注

6
点赞
踩
5

收藏

觉得还不错? 一键收藏
2
评论
运筹系列61：TSP问题GPU算法初探

1. 基础模块1.1 读取数据标准的tsp格式可以用tsplib95库进行读取：import numpy as npimport tsplib95import matplotlib.pyplot as pltdata = tsplib95.load('ALNS/examples/xqf131.tsp')coords = np.array([(coord[0],coord[1]) for city, coord in data.node_coords.items()])solution = t
复制链接

扫一扫

专栏目录