论文拾萃|Solution-based tabu search求解Max-Minsum DP（附代码及详细注释）

void Swap_Tabu_Search(int Sol[], int NSol[], double P[], double* F)
{
 //register表示使用CPU内部寄存器，可以加快变量读写速度
 register int swap, NM;
 register double p, f_c, p_min;
 register double tabu_best_fc, best_fc; //禁忌状态的最优解及当前最优解
 int non_improve = 0; // the stop conditon of TS
 int best_x[50] = { 0 }, best_y[50] = { 0 }, x = 0, y = 0;
 int tabu_best_x[50] = { 0 }, tabu_best_y[50] = { 0 };

 NM = N - M;
 int Hx1 = 0, Hx2 = 0, Hx3 = 0;
 int Hx11 = 0, Hx22 = 0, Hx33 = 0;
 double current_time, starting_time;
 //公式hash funcions
 for (int i = 0;i < M;i++)
 {
  Hx1 += W1[Sol[i]];
  Hx2 += W2[Sol[i]];
  Hx3 += W3[Sol[i]];
 }

 Build_Potential(Sol, NSol, P);
 //拷贝，将当前解作为最优解
 f_best = f;
 for (int i = 0;i < M;i++) {
  SS_BEST.S[i] = Sol[i];
  //cout << SS_BEST.S[i] << " ";
 }
 //cout << endl;
 for (int i = 0;i < NM;i++) {
  SS_BEST.NS[i] = NSol[i];
  //cout << SS_BEST.NS[i] << " ";
 }
 //cout << endl;
 for (int i = 0;i < N;i++) {
  SS_BEST.P[i] = P[i];
  //cout << SS_BEST.P[i] << " ";
 }
 //cout << endl;

 //记录时间
 starting_time = clock();
 int iter = 0;
 current_time = (double)(1.0 * (clock() - starting_time) / CLOCKS_PER_SEC);

 //the number of tabu neighbors and non-tabu neighbors
 int num_tabu_best = 0, num_best = 0;

 while (current_time < time_limit)
 {
  tabu_best_fc = -9999999.0;
  best_fc = -9999999.0;
  num_tabu_best = 0;
  num_best = 0;
  //排序以得高质量子集X、Y，否则全搞m*(n-m)太大了
  for (int i = 0;i < M;i++)
   dataM[i] = P[Sol[i]];
  Quick_Sort(dataM, M, Sol);

  for (int i = 0;i < NM;i++)
   dataNM[i] = P[NSol[i]];
  Quick_Sort(dataNM, NM, NSol);

  //a) Evaluating the neighborhood
  for (int x = 0.3 * M;x >= 0;x--)
  {
   for (int y = 0.7 * NM;y < NM;y++)
   {
    p_min = 9999999.0;

    //其他P[Sol[i]]需将与Sol[x]的距离转换为与Sol[y]的距离
    for (int i = 0;i < x;i++)
    {
     p = P[Sol[i]] - D[Sol[i]][Sol[x]] + D[Sol[i]][NSol[y]];
     p_min = p < p_min ? p : p_min;
    }
    //自身P[Sol[i]]的改变
    p = P[NSol[y]] - D[NSol[y]][Sol[x]];
    p_min = p < p_min ? p : p_min;
    //其他P[Sol[i]]需将与Sol[x]的距离转换为与Sol[y]的距离
    for (int i = x + 1;i < M;i++)
    {
     p = P[Sol[i]] - D[Sol[i]][Sol[x]] + D[Sol[i]][NSol[y]];
     p_min = p < p_min ? p : p_min;
    }

    f_c = p_min;

    Hx11 = Hx1 + (W1[NSol[y]] - W1[Sol[x]]);
    Hx22 = Hx2 + (W2[NSol[y]] - W2[Sol[x]]);
    Hx33 = Hx3 + (W3[NSol[y]] - W3[Sol[x]]);

    //禁忌状态
    if (H1[Hx11 % L] == 1 && H2[Hx22 % L] == 1 && H3[Hx33 % L] == 1)
    {
     if (f_c > tabu_best_fc) //为打破禁忌做准备
     {
      tabu_best_x[0] = x;
      tabu_best_y[0] = y;
      tabu_best_fc = f_c;
      num_tabu_best = 1;
     }
     else if (f_c == tabu_best_fc && num_tabu_best < 50)
     {
      tabu_best_x[num_tabu_best] = x;
      tabu_best_y[num_tabu_best] = y;
      num_tabu_best++;
     }
    }
    else
    {
     if (f_c > best_fc)
     {
      best_x[0] = x;
      best_y[0] = y;
      best_fc = f_c;
      num_best = 1;
     }
     else if (f_c == best_fc && num_best < 50)
     {
      best_x[num_best] = x;
      best_y[num_best] = y;
      num_best++;
     }
    }
   }
  }
  //b. Moves
  if ((num_tabu_best > 0 && tabu_best_fc > best_fc && tabu_best_fc > f_best + 1.0e-7) || num_best == 0) // aspiration criterion
  {
   //无路可走，打破禁忌
   f = tabu_best_fc;
   int select = rand() % num_tabu_best;
   int u = tabu_best_x[select];
   int v = tabu_best_y[select];
   //u,v交换
   Swap_Update_Potential(P, Sol[u], NSol[v]);

   Hx1 = Hx1 + W1[NSol[v]] - W1[Sol[u]];
   Hx2 = Hx2 + W2[NSol[v]] - W2[Sol[u]];
   Hx3 = Hx3 + W3[NSol[v]] - W3[Sol[u]];

   H1[Hx1 % L] = 1;
   H2[Hx2 % L] = 1;
   H3[Hx3 % L] = 1;

   swap = Sol[u];
   Sol[u] = NSol[v];
   NSol[v] = swap;
  }
  // non tabu
  else
  {
   f = best_fc;
   int select = rand() % num_best;
   int u = best_x[select];
   int v = best_y[select];
   Swap_Update_Potential(P, Sol[u], NSol[v]);

   Hx1 = Hx1 + W1[NSol[v]] - W1[Sol[u]];
   Hx2 = Hx2 + W2[NSol[v]] - W2[Sol[u]];
   Hx3 = Hx3 + W3[NSol[v]] - W3[Sol[u]];
   H1[Hx1 % L] = 1;
   H2[Hx2 % L] = 1;
   H3[Hx3 % L] = 1;

   swap = Sol[u];
   Sol[u] = NSol[v];
   NSol[v] = swap;
  }
  //c. print
  iter++;
  current_time = (double)(1.0 * (clock() - starting_time) / CLOCKS_PER_SEC);
  if (f >= f_best + 1.0e-7) //浮点数无法精确表达
  {
   if (f > f_best + 1.0e-7)
   {
    f_best = f;
    for (int i = 0;i < M;i++)
     SS_BEST.S[i] = Sol[i];
    for (int i = 0;i < NM;i++)
     SS_BEST.NS[i] = NSol[i];
    for (int i = 0;i < N;i++)
     SS_BEST.P[i] = P[i];
    SS_BEST.f = f;
    non_improve = 0;
    time_to_target = current_time;
   }
   else if (f == f_best)
    non_improve++;
  }
  else non_improve++;
 }
 for (int i = 0;i < M;i++)  Sol[i] = SS_BEST.S[i];
 for (int i = 0;i < NM;i++) NSol[i] = SS_BEST.NS[i];
 for (int i = 0;i < N;i++)  P[i] = SS_BEST.P[i];
 *F = f_best;
}