论文拾萃|Solution-based tabu search求解Dynamic BDP
论文拾萃|Solution-based tabu search求解Dynamic BDP
“众所周知，既然是在春天，就不要去做秋天的事。”额，不对，拿错剧本了，众所周知管院男女比例令人羡慕，现如今这个班级内部消化问题有待商榷，本文中提到的二部图或对单身狗们有所启发。。。
Part 1 Dynamic bipartite drawing Problem
1.1 图的应用
世间的一切对象都可化为节点；世间一切关系都可化为节点间的一条线；从而组成了如梦幻泡影的图。将来的环球必定是图的世界。 图的应用领域极其广泛，例如 software engineering , VLSI circuit design , or information visualization 皆少不了图的身影。
1.2 研究背景及意义
增量图在图问题中处于一个很重要的地位。
在现实生活的从属关系网络中，个人和集体可以描画为点，而边则代表着个人与集体之间的关系。这些人际网络变化很快，每当有新的成员或集体出现，我们就希望新布局既美观又保持动态稳定性。
动态（增量）二部图问题目前研究较少，而论文作者提出的这种基于解的禁忌搜索（Solution-based Tabu Search）算法从未在图问题中应用过，这对于采用基于解的禁忌搜索研究排列问题而有指导意义。
1.3 解决问题
首先我们科普下二部图，也称二分图，如果我们能将一个图的节点集合分割成两个独立的子集A和B，并使图中的每一条边的两个节点一个来自A集合，一个来自B集合，我们就将这个图称为二分图。
在这个图中，一条边的两端节点，必须来自不同的集合。换句话说，每个节点的邻居，必须是与它来自不同集合的。集合一共有两个。二分图要求一条边的两个端点必须是属于两个不同的集合的：即邻居必须异色如果出现邻居之间同色，则false。
二部图问题（ Bipartite Drawing Problem，简称BDP）是 NP-hard 组合优化问题，其目标在于减少一个二部图中交点的个数，这个交点指的就是边的交点，例如将上图右下角二分图中2、4相连就会得到一个交点。而本问题中，我们考虑二部图问题的动态变化，这就成为了一个变种问题，称为Dynamic BDP（简称DBDP），即在一个给定的二部图中增加顶点及与之相关的边来得到一个新的二部图，值得强调的是，我们要保证新图和原图在一定程度上的相似性，也就是在新的图中，旧的点顺序不变。因而作者提出了一种结合adaptive memory的禁忌搜索方法即基于迭代解的禁忌搜索（ISB-TS）来高效率的搜索解空间，并将其 explict memory 与现有基于 attributive memory 的最佳方法进行了比较。
1）BDP 又是个啥？ 首先，阐述经典的二部图问题（BDP）：二部图被定义为 G=(V1,V2,E)，就是顶点集 V 可分割为两个互不相交的子集，并且图中每条边（E）依附的两个顶点都分属于这两个互不相交的子集 V1、V2，两个子集内的顶点不相邻。BDP 的一个解表示为 D = (π_1 , π_2 )，左点集中的顶点 u 和右点集中的顶点 v 所处位置即可用π_1(u)和π_2(v)来表达，若 u 在图中的位置先于 v，则π_1(u)<π_1(v)，反之 π_1(u)>π_1(v) 。
2）DBDP 的介绍 现在，我们介绍动态二部图问题（DBDP），通过增加两个点集AV_k及他们相关的边集AE_k ( k = 1 , 2 )得到一个增量图。表达方式同理，具体见图。
为方便表述，后文中将原图中的顶点称为原点，新加入的点称为增点。
正如上图所示，黑点为原点，须保持相对位置不变，白点为增点，可以任意移动。不难看出，从左图到右图的过程就改变了五个点的位置，你是不是很容易就接受了这种变化，正是因为布局改动很小，正如看电影一般，背景未变的情况下你便能更专注于人物对话。 图从左到右的这种改变让边交点数显著减少，DBDP 的目标便是找到一种“制图技巧”使边的交点数最少同时保持原点的相对位置。
在这个问题中，边和点都是已知的，也即对于增点而言，我们已知他们和原点之间的两两配对关系，点已经划分为两个集合，我们需要做的就是决定点的顺序，这个顺序就是决策变量，目标函数就是边的交点数。
Part 2 禁忌搜索算法(TS)再回顾
对于禁忌搜索各部分内容，公众号之前的推文中已经介绍过多次了。还不了解的同学可以先移步干货 | 到底是什么算法，能让人们如此绝望？了解算法细节。
Part 3 具体算法介绍
作者提出的ISB-TS算法采用了迭代局部搜索方法的基本方案。采用基于解的禁忌搜索来加强对解空间的搜索，并且采用了一种自适应扰动机制，使得每当得到局部最优，就会移动到一个新的搜索区域，避免陷入局部最优。而后者也是ISB-TS算法的侧重所在，在某个搜索阶段使解适当程度的多样化。其大体架构如下：
大体过程就是应用基于解的禁忌搜索方法来获得局部最优解S(第6行)。然后，更新到目前为止找到的最优解，并且更新连续的未改进局部最优的参数ξ(第7-12行)。在每次局部优化之后，ISB-TS根据搜索的状态，通过采用不同的扰动来尝试从一个局部最优移动到另一个局部最优。如果找到更好的局部最优值S，则ISB-TS切换到较弱的扰动，否则ISB-TS通过增加扰动强度L来更强烈地扰动S (第15-16行)。在访问了一定数量的局部最优解而没有改进迄今找到的最优解后，应用明显更强的扰动来最终将搜索推向搜索空间中的新的、更远的区域(第17-20行)。当确定扰动强度时，使用扰动算子在搜索空间的不同区域中生成新的解。基于解的禁忌搜索和自适应扰动机制迭代交替，直到满足停止标准(即最大计算时间T_{max})，并最终返回最佳解S_{BEST}作为最终结果(第24行)。
在这篇文章里，文献的方法即利用 Solution-based tabu search 来求解Dynamic bipartite drawing Problem，该算法依靠hash vectors来有效地确定候选解的禁忌状态，并通过一个受约束的交换邻域来保证算法的高计算效率。
3.1 初始解生成
ISB-TS 采用了贪婪随机自适应搜索(GRASP)来构造初始解。
对GRASP不熟悉的朋友可以通过推文Greedy Randomized Adaptive Search 算法超详细解析，附代码实现TSP问题求解了解一下。
利用贪婪随机自适应搜索，我们从一个包含所有原点的部分解出发（或者说未完成的解，原谅我表述不是很准确），然后不断迭代插入增点最后完成完整解的构造。其伪代码如下：
在每次迭代中，集合CL(即candidate list)记录AV中剩余顶点的所有对及其对应的可行位置，值得注意的是，在CL中生成顶点的受限候选列表RCL(第6行)。
3.2 禁忌搜索阶段
ISB-TS 算法采用约束邻域结构和 Fast neighborhood evaluation strategy 以及专门的基于解的禁忌策略得到局部最优。
关于去重优化以及哈希函数的应用这里就不再详细介绍，详见之前推文论文拾萃|Solution-based tabu search 求解 Max-Minsum DP（附代码及详细注释）
3.2.1 邻域结构
我们将一个顶点从它当前位置到新的位置的距离定义为距离门槛 Ф，超过这个门槛的移动将会被禁止，从而约束邻域动作。提出的约束邻域结构由两种邻域动作构成，即 CN(S)=N_1(S)∪N_2(S)。其中邻域动作N_1表示将一个增点插入到满足距离限制的另一个顶点的前面或后面。对应邻居解可表示为
邻域动作N_2表示将一个原点插入到其他任意一个增点的前面或后面同时保证原点的次序满足距离限制。
其中Insert(v,u)代表将顶点v从在当前解 S 中的位置插入到顶点u所在位置前后，若原先u在前，则插入到 u前面，反之插入到 u 后面，特别的，如果v是增点，邻域动作移动到的位置就有n_1-1(n_2-1)种可能，【友情提示：n_1、n_2为左右两点集顶点数目】，当然，我们只保留可行解，并且由于Ф的存在，时间复杂度就会被限制在 O(n^2)。
3.2.2 Fast neighborhood evaluation strategy
由于大部分顶点相对位置都未改变，所以去重优化便可派上用场，大幅度减少时间复杂度。我们给定同一点集的两个顶点u,v，NC_{uv} 表示与u点相连的边和与 v点相连的边的交点数，【强调：暗含限制-点 u在点 v 前面】，则我们可用矩阵 M_k(u,v) 记录边交点数目。
现在举个例子便于读者理解
顶点 3 相关边 (3,C) 和 (3,D) 以及顶点 5 相关边 (5,B) 的交点数为 2，则 M_1(3,5)=2，而这里我们就来解答之前强调 π_1(u)<π_1(v) 的原因，因为 M_1(5,3) 就代表着将顶点 5 插入到顶点 3 前面后相关边的交点数，不难验证为 0，具体如下图所示，读者不妨一一验证。
基于此，我们计算目标函数就会变得非常轻松，我也不多费口舌。
3.2.3 基于解的禁忌搜索
关于哈希函数的应用在另一篇推文中已有介绍，这里只放哈希值的计算公式
其中参数 k 用来定义哈希函数，参数 λ 为哈希向量的长度。值得一提的是，哈希值的更新可以在 O(1) 时间内完成，如下图所示。
禁忌搜索部分伪代码如下：
3.3 自适应扰动机制
自适应扰动机制是为了逃离局部最优以找到更优的解而提出，并且根据当前所处的阶段不断变换扰动相关参数来加强搜索的机制。它由 destruction 和 reconstruction 两部分构成，即从解 S 移除一些增点和迭代的重新插入当前解。伪代码如下：
Part 4 核心代码分享
以下为求解多维背包问题的核心代码（同样是two-stage tabu search algorithm）：
void Swap_tabu_search(Solution &S, double X)
{
int i, j, k, k0, x, y, v, u, iter ;
int K;
int count;
int assign, Ncount, Ncount1, Ncount2;
double pho, eta;
int num, num1;
long int sum_penalty;
long int delt, tabu_delt, non_delt;
int select, swap;
int non_improve = 0 ; //the stop condition of TS
long int tabu_best_fc, best_fc, f_c ;
int num_tabu_best, num_best; //the number of tabu neighbors and non-tabu neighbors
Neighbor best[ 50 ];
Neighbor tabu_best[ 50 ];
int Hx1, Hx2, Hx3;
int Hx11, Hx22, Hx33;
double current_time, starting_time;
eta = 0.3;
pho = 0.3;
Build_Information();
Hx1 = 0;
Hx2 = 0;
Hx3 = 0;
for(i=0;i<N;i++)
{
Hx1 += W1[i]*(S.X[i]);
Hx2 += W2[i]*(S.X[i]);
Hx3 += W3[i]*(S.X[i]);
}
f = S.f;
f_best = S.f;
S_BEST.IN = S.IN;
S_BEST.ON = S.ON;
S_BEST.f = S.f;
for(i=0;i<N;i++) S_BEST.S[i] = S.S[i];
for(i=0;i<N;i++) S_BEST.NS[i] = S.NS[i];
for(i=0;i<N;i++) S_BEST.X[i] = S.X[i];
for(i=0;i<M+Q;i++) S_BEST.IW[i] = S.IW[i];
iter = 0;
non_improve = 0;
starting_time = clock();
current_time = (double) (1.0*(clock()-starting_time)/CLOCKS_PER_SEC);
while( current_time < X )
{
tabu_best_fc = -99999999 ;
best_fc = -99999999 ;
num_tabu_best = 0 ;
num_best = 0 ;
num = 0;
num1 = 0;
//2) Evaluating the neighborhood N3
for(x = 0; x < S.IN; x++)
{
for(y = 0; y < S.ON; y++)
{
f_c = f + (P[S.NS[y]]-P[S.S[x]]) ;
if(f_c <= f_best) continue;
sum_penalty = 0;
for(i=0;i<M;i++) if(S.IW[i]+(R[i][S.NS[y]]-R[i][S.S[x]]) > B[i]) sum_penalty += (S.IW[i] + R[i][S.NS[y]] - R[i][S.S[x]] - B[i]);
for(i=M;i<M+Q;i++) if(S.IW[i]+(R[i][S.NS[y]]-R[i][S.S[x]]) < B[i]) sum_penalty += (B[i] - S.IW[i] - R[i][S.NS[y]] + R[i][S.S[x]]);
delt = -1.0*lanbda*sum_penalty;
Hx11 = Hx1 + (W1[S.NS[y]]- W1[S.S[x]]) ;
Hx22 = Hx2 + (W2[S.NS[y]]- W2[S.S[x]]) ;
Hx33 = Hx3 + (W3[S.NS[y]]- W3[S.S[x]]) ;
if( H1[Hx11%L]==1 && H2[Hx22%L]==1 && H3[Hx33%L]==1 ) // if it is tabu
{
num ++;
if( f_c + delt > tabu_best_fc )
{
tabu_best[ 0 ].x1 = x ;
tabu_best[ 0 ].y1 = y ;
tabu_best[ 0 ].type = 2 ;
tabu_best[ 0 ].f = f_c ;
tabu_best_fc = f_c + delt ;
tabu_delt = delt;
num_tabu_best = 1 ;
}
else if( f_c + delt == tabu_best_fc && num_tabu_best < 50 )
{
tabu_best[ num_tabu_best ].x1 = x;
tabu_best[ num_tabu_best ].y1 = y;
tabu_best[ num_tabu_best ].type = 2;
tabu_best[ num_tabu_best ].f = f_c;
num_tabu_best++ ;
}
}
else
{
num1++;
if( f_c + delt > best_fc )
{
best[ 0 ].x1 = x;
best[ 0 ].y1 = y;
best[ 0 ].type = 2;
best[ 0 ].f = f_c;
best_fc = f_c + delt;
non_delt = delt;
num_best = 1 ;
}
else if( f_c + delt == best_fc && num_best < 50 )
{
best[ num_best ].x1 = x ;
best[ num_best ].y1 = y ;
best[ num_best ].type = 2 ;
best[ num_best ].f = f_c;
num_best++ ;
}
}
}
}
//4) moves
if( num_best == 0 ) // aspiration criterion
{
select = rand() % num_tabu_best ;
f = tabu_best[select].f ;
u = tabu_best[ select ].x1;
v = tabu_best[ select ].y1;
Hx1 += (W1[S.NS[v]]-W1[S.S[u]]);
Hx2 += (W2[S.NS[v]]-W2[S.S[u]]);
Hx3 += (W3[S.NS[v]]-W3[S.S[u]]);
H1[Hx1%L] = 1;
H2[Hx2%L] = 1;
H3[Hx3%L] = 1;
for(i=0;i<M+Q;i++)
{
S.IW[i] += (R[i][S.NS[v]] - R[i][S.S[u]]);
}
S.X[S.S[u]] = 1 - S.X[S.S[u]];
S.X[S.NS[v]] = 1 - S.X[S.NS[v]];
swap = S.S[u] ;
S.S[u] = S.NS[v];
S.NS[v] = swap ;
S.f = f;
}
else // non tabu
{
select = rand() % num_best ;
f = best[ select ].f ;
u = best[ select ].x1;
v = best[ select ].y1;
Hx1 += (W1[S.NS[v]]-W1[S.S[u]]);
Hx2 += (W2[S.NS[v]]-W2[S.S[u]]);
Hx3 += (W3[S.NS[v]]-W3[S.S[u]]);
H1[Hx1%L] = 1;
H2[Hx2%L] = 1;
H3[Hx3%L] = 1;
for(i=0;i<M+Q;i++)
{
S.IW[i] += (R[i][S.NS[v]] - R[i][S.S[u]]);
}
S.X[S.S[u]] = 1 - S.X[S.S[u]];
S.X[S.NS[v]] = 1 - S.X[S.NS[v]];
swap = S.S[u] ;
S.S[u] = S.NS[v];
S.NS[v] = swap ;
S.f = f;
}
//5. print
iter ++ ;
current_time = (double) (1.0*(clock()-starting_time)/CLOCKS_PER_SEC);
if( f >= f_best )
{
sum_penalty = 0;
for(i=0;i<M;i++) if(S.IW[i] > B[i]) sum_penalty += (S.IW[i] - B[i]);
for(i=M;i<M+Q;i++) if(S.IW[i] < B[i]) sum_penalty += (B[i] - S.IW[i]);
//printf("sum =%d\n",sum_penalty);
if( f > f_best && sum_penalty == 0 )
{
f_best = f;
for( i = 0 ; i < N ; i ++ )
S_BEST.X[ i ] = S.X[ i ] ;
for( i = 0 ; i < N ; i ++ )
S_BEST.S[ i ] = S.S[ i ] ;
for( i = 0 ; i < N ; i ++ )
S_BEST.NS[ i ] = S.NS[ i ] ;
for( i = 0 ; i < M+Q ; i ++ )
S_BEST.IW[ i ] = S.IW[ i ] ;
S_BEST.IN = S.IN;
S_BEST.ON = S.ON;
S_BEST.f = f;
// printf("\n swap_move : %8d %d %d %d %lf %lf", iter, f, f_best, non_improve, pho, current_time);
non_improve = 0 ;
time_one_run = current_time;
}
else if ( f == f_best ) non_improve ++ ;
}
else non_improve ++ ;
}
for(i=0;i<N;i++) S.X[i] = S_BEST.X[i];
for(i=0;i<N;i++) S.S[i] = S_BEST.S[i];
for(i=0;i<N;i++) S.NS[i] = S_BEST.NS[i];
for(i=0;i<M+Q;i++) S.IW[i] = S_BEST.IW[i];
S.f = S_BEST.f;
S.IN = S_BEST.IN;
S.ON = S_BEST.ON;
} // With a fixed value of K.
欲下载本文相关代码，请移步留言区
参考文献：
B. Peng, D. Liu, Z. Lü, R. Martí and J. Ding "Adaptive memory programming for the dynamic bipartite drawing problem"Information Sciences 2020 Vol. 517 Pages 183-197
-The End-
文案&排版：朱正雄（华中科技大学管理学院本科一年级）
指导老师：秦虎老师 （华中科技大学管理学院）
审稿人：周航（华中科技大学管理学院本科二年级）
如对文中内容有疑问，欢迎交流。PS：部分资料来自网络。 如有需求，可以联系： 秦虎老师（华中科技大学管理学院：professor.qin@qq.com） 朱正雄（华中科技大学管理学院本科一年级：2627889552@qq.com）
周航（华中科技大学管理学院本科二年级：zh20010728@126.com）