旅行商问题2-OPT算法的并行与优化
GCC-6.2.0
OpenMPI/2.0.0
OpenMp 4.5 (2015-11)
介绍
废话不多说，查阅下面链接。
串行2-OPT的思路如下：
假如我们有{0, 1, 2, 3, 4, 5}这5个城市，初始路线为
0-1-2-3-4-5-0。我们通过两个指针i和k各指向一个城市进行遍历。每次遍历，我们会尝试对当前的**(i,k)对**做一次2optSwap（后面简称reverse），将ik之间的路线反转。
这个操作将会重复，直到我们无法找到更短的路径，我们便找到了局部最优解，也就是最终的结果。
- 例如我们遍历到i=1, k=3。我们对这个(i,k)对进行一次reverse，得到一个新路线：
0-3-2-1-4-5-0，如果新路线可以比之前的路线更短，我们就记下当前的i和k，以及reverse之后能缩短的距离。然后我们继续遍历，最后选择最短的那一条结果进行reverse。 - 这里原版的串行化2-opt会在找到第一个可以缩短路径的(i, k)之后马上reverse，然后将i和k设为初始值（i=1,k=2）从头开始遍历，并不会每次遍历所有城市然后取当前的最优解。如果使用原版，我们将无法对比并行优化的结果，因为每次的跳转完全随机且取决于具体的数据。
- 在计算新旧路线的距离差时有一个细节，即我们不需要计算路线的总长。以上面i=1,k=3为例，我们只需要计算distance(0, 1), distance(0, 3), distance(1, 4), distance(3, 4)。而其他路线是不变的。
2-opt-WikiPedia的伪代码，
//最直观的伪代码
procedure 2optSwap(route, i, k) {
1. 将 route[0] 到 route[i-1] 按顺序加入new_route
2. 将 route[i] 到 route[k] 倒序加入new_route
3. 将 route[k+1] 及之后的元素按顺序加入new_route
return new_route;
}
//稍微贴近实现的伪代码
repeat until no improvement is made {
start_again:
best_distance = calculateTotalDistance(existing_route)
for (i = 1; i <= number of nodes eligible to be swapped - 1; i++) {
for (k = i + 1; k <= number of nodes eligible to be swapped; k++) {
new_route = 2optSwap(existing_route, i, k)
new_distance = calculateTotalDistance(new_route)
if (new_distance < best_distance) {//如上所述，本文中的实现在这里稍有区别
existing_route = new_route
best_distance = new_distance
goto start_again
}
}
}
}
串行2-opt的C语言实现，参照伪代码并修改为每次遍历所有城市，选择这次遍历中找到的最短路径进行reverse：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <sys/time.h>
#include <stdint.h>
#define swap(a, b) {float tmp = a; a = b; b = tmp;}
void reverse(int len, int* route) {
for (int i=0; i<len/2; i++) {
swap(route[i], route[len-1-i]);
}
}
int get_distance(int x1,int y1, int x2, int y2){
return (int)sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
}
int TSPDistance(int size, int route[size+1], int x[size], int y[size]) {
int sum = 0;
for (int i=1;i<size+1;i++) {
sum += get_distance(x[route[i-1]], y[route[i-1]], x[route[i]], y[route[i]]);
}
return sum;
}
void twoOptTSP(int size, int route[size+1], int x[size], int y[size]) {
//initial route
for (int i=0;i<size;i++) {
route[i] = i;
}
route[size]=0;
int change;
int minchange;
int mini;
int mink;
int improvement = 0;//c
LOOP: do {//if negative changes, then it is still improving
minchange=0;
mini=0;
mink=0;
improvement = 0;
for (int i = 1; i < size-1; i++) {
for (int k = i + 1; k < size; k++) {
change = get_distance(x[route[i-1]], y[route[i-1]], x[route[k]], y[route[k]])
+ get_distance(x[route[i]], y[route[i]], x[route[k+1]], y[route[k+1]])
- get_distance(x[route[i-1]], y[route[i-1]], x[route[i]], y[route[i]])
- get_distance(x[route[k]], y[route[k]], x[route[k+1]], y[route[k+1]]);
if (change < minchange) {
improvement = 1;
minchange=change;
mini=i;
mink=k;
}
}
}
// printf("minchange %d i %d k %d\n", minchange, mini, mink);
reverse(mink-mini+1, &route[mini]);
} while(improvement);
int curr_distance = TSPDistance(size, route, x, y);
printf("curr_distance %d\n", curr_distance);
}
针对这个串行版本，我们将设计和实现2-OPT的分布式并行算法。算法用C语言实现，我们会使用openMP进行单机多核并行，使用openMPI进行多机并行。
实现思路
并行化这个算法，最直观的切入点便是：将指针i或k的遍历分配给不同的核。
仔细分析后，我们可以发现，因为k的初始值为i+1，所以随着i的增大，k的遍历数量会越来越小。
如果我们将i简单拆分给不同的机器，那么每个机器的工作量是不一样的。而工作量不同，必然导致有的机器先到达同步点，等待其他机器完成计算，造成资源浪费。因此我们要进行负载均衡。
i | 1 2 3 4 | 5 6 7 8 | 9 10 11 12 | 13 14 15 16 |
machine A B C D
本文中使用的方案是：像高斯求和等差数列一样，从头部和尾部同时分派i的任务。这样，每台机器需要计算的(i,k)对是相同的, 如果城市数无法整除机器数，最多只会造成一台机器的资源浪费。即：
i | 1 2 | 3 4 | 5 6 | 7 8 | 9 10 | 11 12 | 13 14 | 15 16 |
machine A B C D D C B A
如果机器是多核，我们可以进一步将k的遍历拆分给每一个核。
在遍历完i,k之后，我们得到了每个核中的最小值。我们需要将他们汇总在机器A，即master上进行比较，获得全局最小值。然后将reverse的结果分发回每个机器。
这里有一个细节，即我们仅发送最后要进行变更的i和k，而不发送master变换后的整个路径。这样可以最大程度减少通信时间。
以上，我们就完成了算法的并行化。接下来我们来对这个并行算法进行针对性优化。
优化
仔细观察后我们就能够发现，2-opt算法中存在大量互相不干扰的(i,k)对（下面简称兼容），遍历一次然后reverse其中最优的一个，和在同一次遍历中reverse多个兼容的ik对，得到的结果是一样的。下面是几个兼容的情况：
(i,k)范围不重叠
假设以下的情况：
我们有路径：0-1-2-3-4-5-6-7-8-9-0，我们在master汇总了三个(i, k)对：(1,2), (4,5), (7,8)。按照上面的算法，我们只会选择那个能得到最短的路径的ik进行reverse。
但实际上，我们换了其中一个ik对，假如是(1,2)，reverse后得到0-2-1-3-4-5-6-7-8-9-0。然后我们在下一次遍历后，我们仍然会得到另外两个ik对，即(4,5), (7,8)。而如果我们直接对上面的三个ik同时进行reverse，就节省了两次遍历。
因此，如果我们得到的多个ik之间的范围不重叠，就可以同时进行reverse。
- 这里有一个细节，如果我们获得的是
(1,2), (3,4)， 他们看似兼容，实际上需要计算的城市是不同的。 - reverse
(1,2)之前，路径为0-1-2-3-4-5-6-7-8-9-0，计算reverse(3,4)能缩短的距离，需要计算的是distance(2,3), distance(4,5), disntance(2,4), disntance(3,5) - reverse
(1,2)之后，路径变成了0-2-1-3-4-5-6-7-8-9-0，这时再计算(3,4)能缩短的距离，则需要计算distance(1,3), distance(4,5), distance(1,4), distance(3,5)
(i,k)有包含关系
假如我们还是有路径``0-1-2-3-4-5-6-7-8-9-0`，这次我们在master汇总了这样三个ik对：
(1,6), (3,4), (4,7)。仔细观察会发现，(3,4)是被包含在(1,6)里面的。但是(4,7)又与(1,6)不兼容。如果reverse(4,7)之后路线最短，那我们显然应该直接reverse(4,7)，然后进行下一次遍历。可以，并且reverse(1,6)之后路线是最短的，我们是不是可以同时reverse(1,6), (3,4)呢？
我们需要分情况讨论：
首先，我们reverse(1,6)之后的结果是：0-6-5-4-3-2-1-7-8-9-0。那么下一次遍历，我们需要比较reverse(3,4)和(4,7)谁更短。
如果(3,4)更短，我们就可以同时reverse(1,6), (3,4)，反之则不行。
实现方法
将上面的情况推到一般，我们实际需要做的就是对所有汇总到master的结果按照reverse后路线长度从短到长进行排序。求出从下标0开始的最长子序列，序列中的每一个元素都与在它之前的所有元素兼容。
比如说排序后结果为[(1,2), (5,6), (4,8), (10,11)]，我们从(5,6)开始遍历，发现(5,6)与(1,2)兼容，而(4,8)虽然与(1,2)兼容，但不与(5,6)兼容。所以我们就不用再检查(10,11)了，即使它与第一和第二个元素兼容，我们也不能保证在reverse第一和第二个元素后，不会出现比(10,11)更优且不兼容的解。
当然了，如果我们不在意最后的结果是否与串行化算法的结果一致，我们也可以直接reverse所有与第一个元素兼容的元素（因为烦人的死线，本文中仅实现了这种算法，当然这个得到的结果是不那么严谨的）。
伪代码
//parallel pseudocode
2opt(x, y, route):
do:
//phase 1
minchange[number of threads] = (0, 0, 0);//三个值分别为：(change in distance, i, k)
for ii = 0 to (size-1)/2 do in parallel:
processor p thread t does:
i = ii + 1 //负载均衡
for k = i + 1 to size-1:
change = get_change(route, i, k)
if (change < minchange[0])
minchange[t] = (change, i, k);
i = size-2-ii //负载均衡
for k = i + 1 to size-1:
change = get_change(route, i, k)
if (change < minchange[t][0])
minchange = (change, i, k);
//phase 2
gather minchange from all processors to processor_results of processor 0
processor 0 does:
add changes that can be performed together to send_results
broadcast send_results to all processors
//phase 3
all processor does locally:
for change in send_results:
reverse(route, change[1], change[2])
while there is improvement in distance
//寻找与最优解兼容的解
buildChangeMat(idx, minchange, mat):
If mat[idx].length == 0:
add minchange to mat[idx];
return;
Check idx:
If (not compatible):
return;
Else if (compatible):
add minchange to mat[idx];
return;
Else if (inside range of an element):
buildChangeMat(element, minchange, mat);
串行版，并行版和编译运行的shell代码，我都贴到文件里，有需要自取。
不太严谨的实验结果
| 城市数量 | 串行结果 （距离） | 并行结果（距离） |
|---|---|---|
| 1000 | 25287 | 25077 |
| 2000 | 35549 | 35373 |
| 3000 | 42951 | 43108 |
| 5000 | 55260 | 55240 |
Table 1 两台机器每个12核的运行结果
这里可以看到结果是不一样的。因为我们实现的是不严谨的版本。
| 城市数量 | 串行时间（s） | 并行时间（s） | Speedup | 并行效率 |
|---|---|---|---|---|
| 1000 | 16.40 | 0.50 | 32.80 | 1.37 |
| 2000 | 133.45 | 1.54 | 86.66 | 3.61 |
| 3000 | 473.46 | 4.42 | 107.12 | 4.46 |
| 5000 | 2259.01 | 51.02 | 44.28 | 1.85 |
Table 2 两台机器，每个12核的运行时间
| 城市数量 | 串行结果（距离） | 并行结果（距离） |
|---|---|---|
| 5000 | 55260 | 55235 |
Table 3 四台机器每个12核的运行结果
| 城市数量 | 串行时间（s） | 并行时间（s） | Speedup | 并行效率 |
|---|---|---|---|---|
| 5000 | 2249.84 | 20.81 | 108.11 | 2.25 |
Table 4 四台机器每个12核的运行时间
最后
总的来说，这个算法应该还有优化的空间。目前每个线程只会提供一个(i,k)对。我们可以把它变成两个或者更多，这样我们就有可能获得更多兼容的解。同时兼容的情况也可能不止上文中提到的两种。
我们可以赠机器只有n个局部最小值，我们可能错过了一些可以同时在一轮迭代中reverse的机会。如果能够储存更多的候选路径，我们可能可以获得更大的speedup。
不过，这个并行算法的缺憾也很直观：
- 一个是speedup取决于数据是否提供了在一轮迭代中进行多次变换的可能，所以speedup并不稳定。不过如果城市数量足够大，肯定会在早期的遍历中获得巨大的speedup。
- 另一个就是没有实现与串行化结果一致的版本。我们无法获得与串行算法相同的值，也就无法更准确地计算speedup。因为串行和并行的结果不同，所以我们无法确保双方的总计算量相同，有可能并行算法多进行了或者少进行了几轮迭代，才得出了最终的speedup。
- 不过总的来说，我们也可以看到这个并行算法取得的speedup有时完全超出预期，并行效率（=speedup/总核数）远大于1。
扫一扫
1548
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
