模擬退火算法(SA)和迭代局部搜索(ILS)求解TSP的Java代碼分享

大家好呀！我們你們好久不見(jiàn)的。。。咳咳，初次見(jiàn)面的小編！

之前重新整理了ILS的代碼，有人留言問(wèn)能不能提供java版。

正好最近在學(xué)啟發(fā)式算法和java，為了造福人類小編打算提供模擬退火法和迭代局部搜索求解TSP的java版本，方便一些不喜歡C++的同鞋~~

代碼是基于我自己寫(xiě)的版本，但我是學(xué)習(xí)了公眾號(hào)推文之后寫(xiě)的，同時(shí)有參照原文代碼，因此雖然與C++代碼有些許區(qū)別，但總體是類似的。

本文不再贅述TSP或者SA，ILS了，有需要請(qǐng)點(diǎn)擊下方鏈接學(xué)習(xí)（一看就懂的那種哦！）：

干貨 | 用模擬退火(SA, Simulated Annealing)算法解決旅行商問(wèn)題

干貨|迭代局部搜索算法(Iterated local search)探幽（附C++代碼及注釋）

不多說(shuō)了，開(kāi)始看代碼吧~！

SA求解TSP的JAVA代碼

SA分為四個(gè)類：MainRun，Data，Path，SimulatedAnnealing。

?

MainRun是程序的入口。

package SA;
/** * 主函數(shù)運(yùn)行類 */
public class MainRun {
  public static void main (String args[])  {    long begintime = System.nanoTime();
    Data.PrintData();    SimulatedAnnealing.SA();    SimulatedAnnealing.PrintPath();
    long endtime = System.nanoTime();    double usedTime= (endtime - begintime)/(1e9);    System.out.println();    System.out.println("程序耗時(shí)："+usedTime+"s");  }}

Data是數(shù)據(jù)類，存放SA和TSP的數(shù)據(jù)。

?

package SA;
/** * 數(shù)據(jù)類，包括： * TSP城市坐標(biāo)，采用柏林52城 * SA系數(shù)。 */

public class Data {  public static final double  T0=50000.0, // 初始溫度  T_min=1e-8,  // 結(jié)束溫度  q=0.98,   // 退火系數(shù)  K=1;  //公式中的常數(shù)K   public static final int L=1000,  // 每個(gè)溫度時(shí)的迭代次數(shù)，即鏈長(zhǎng)  N=52;  // 城市數(shù)量
  public static double[][] city_pos= //柏林52城城市坐標(biāo)，最優(yōu)解7542  {  { 565,575 },{ 25,185 },{ 345,750 },{ 945,685 },{ 845,655 },  { 880,660 },{ 25,230 },{ 525,1000 },{ 580,1175 },{ 650,1130 },  { 1605,620 },{ 1220,580 },{ 1465,200 },{ 1530,5 },  { 845,680 },{ 725,370 },{ 145,665 },  { 415,635 },{ 510,875 },{ 560,365 },{ 300,465 },{ 520,585 },{ 480,415 },  { 835,625 },{ 975,580 },{ 1215,245 },{ 1320,315 },{ 1250,400 },{ 660,180 },  { 410,250 },{ 420,555 },{ 575,665 },{ 1150,1160 },{ 700,580 },{ 685,595 },  { 685,610 },{ 770,610 },{ 795,645 },{ 720,635 },{ 760,650 },{ 475,960 },  { 95,260 },{ 875,920 },{ 700,500 },{ 555,815 },{ 830,485 },{ 1170,65 },  { 830,610 },{ 605,625 },{ 595,360 },{ 1340,725 },{ 1740,245 }  };
  public static void PrintData(){    System.out.println("模擬退火算法，初始溫度T0="+T0+        ",降溫系數(shù)q="+q+",每個(gè)溫度迭代"+L+"次");  }
}

Path是路徑類，打包處理路徑的靜態(tài)方法。

?

package SA;import static SA.Data.*;import static java.lang.Math.*;
/** * 路徑類，打包處理路徑的靜態(tài)方法： * 計(jì)算兩點(diǎn)間距離 distance * 計(jì)算路徑長(zhǎng)度 path_len * 產(chǎn)生新解（新路徑）create_new */
public class Path {  public static double distance(double[] p1,double[] p2)  //計(jì)算兩點(diǎn)間距離 {    double dis=0;    dis=sqrt(pow(p1[0]-p2[0],2)+pow(p1[1]-p2[1],2));    return dis;  }
  public static double path_len(int[] city_list) //計(jì)算路徑長(zhǎng)度 {    double path=0;    for (int i=0;i<(N-1);i++)      path+=distance(city_pos[city_list[i]],city_pos[city_list[i+1]]);    path+=distance(city_pos[city_list[0]],city_pos[city_list[N-1]]) ;    return path;  }
  public static void create_new(int[] city_list)   //產(chǎn)生新解{     int i=(int) (random()*N);     int j=(int) (random()*N);    while(j==i)      j=(int) (random()*N);
    int temp;     temp=city_list[i];     city_list[i]=city_list[j];     city_list[j]=temp;???}}

SimulatedAnnealing開(kāi)始模擬退火。

package SA;import static SA.Data.*;   import static SA.Path.*;import static java.lang.Math.*;import java.util.*;
/** * 模擬退火過(guò)程 */
public class SimulatedAnnealing {
  private static int[] bestpath=new int[N];  private static int count=0;  private static double f2;
  public static void SA() //模擬退火{    for(int i=0;i      bestpath[i]=i;
     int[] newpath=Arrays.copyOf(bestpath,bestpath.length);   // 新解
     double f1;    f2=path_len(bestpath);     double T=T0;
     while(T>T_min)     {       for (int i=0;i      {         create_new(newpath);         f1=path_len(newpath);         double de=f1-f2;         if (de<0)         {           System.arraycopy(newpath, 0, bestpath, 0, bestpath.length);          f2=f1;           }        else          {           double r = random();           if(exp(de/(K*T))           {             System.arraycopy(newpath, 0, bestpath, 0,bestpath.length);             f2=f1;             }             else               System.arraycopy(bestpath, 0, newpath, 0, bestpath.length);         }       }       T*=q;         count++;     }  }
  public static void PrintPath() //打印最優(yōu)解{    System.out.println("共降溫"+count+"次，得到的TSP最優(yōu)距離為"+f2+"路徑為");     for(int j=0;j     {       if (j==0)         System.out.print(bestpath[j]);       else         System.out.print("--->"+bestpath[j]);     }  }}

ILS求解TSP的JAVA代碼

ILS分為MainRun，City，Delta，Perturbation，LocalSearch，Solution。

?

主函數(shù)運(yùn)行類，包括了ILS總方法。

?

package ILSTSP;import static ILSTSP.City.*; import static ILSTSP.Solution.*; import static ILSTSP.Perturbation.*;import static ILSTSP.LocalSearch.*;
/** * 主函數(shù)運(yùn)行類，包括了ILS總方法以及計(jì)時(shí)功能。 */
public class MainRun {  public static void main(String args[])  {    int max_iterations = 600;    int max_no_improve = 50;    long begintime = System.nanoTime();
    Solution best_solution=new Solution();    iterated_local_search(best_solution, max_iterations, max_no_improve);
    System.out.println();    System.out.println("搜索完成！最優(yōu)路線總長(zhǎng)度 = "+best_solution.cost);    System.out.println("最優(yōu)訪問(wèn)城市序列如下：" );      for (int i = 0; i < CITY_SIZE;i++)        System.out.print(best_solution.permutation[i]+"\t");
      long endtime = System.nanoTime();    double usedTime= (endtime - begintime)/(1e9);    System.out.println();    System.out.println("程序耗時(shí)："+usedTime+"s");  }
  static void iterated_local_search(Solution best_solution, int max_iterations, int max_no_improve)  {      Solution current_solution = new Solution();
      //獲得初始隨機(jī)解      random_permutation(best_solution.permutation);
      best_solution.cost = cost_total(best_solution.permutation);      local_search(best_solution, max_no_improve); //初始搜索
      for (int i = 0; i < max_iterations; i++)      {          perturbation(best_solution,current_solution); //擾動(dòng)+判斷是否接受新解          local_search(current_solution, max_no_improve);//繼續(xù)局部搜索
          //找到更優(yōu)解          if (current_solution.cost < best_solution.cost)          {              for (int j = 0; j < CITY_SIZE; j++)              {                  best_solution.permutation[j] = current_solution.permutation[j];              }              best_solution.cost = current_solution.cost;          }
          System.out.println("迭代搜索 "+i+ " 次\t" +          "最優(yōu)解 = "+ best_solution.cost+          "當(dāng)前解 = "+ current_solution.cost               );      }    }}

City類存放城市數(shù)據(jù)，柏林52城坐標(biāo)。

?

package ILSTSP;
import static java.lang.Math.*;
/** * TSP數(shù)據(jù)類， * 存放柏林52城坐標(biāo)， * 兩點(diǎn)間距離計(jì)算distance_2city。 */
public class City {
  public static int CITY_SIZE=52;  // 城市數(shù)量
  public static int[][] cities=    {    { 565,575 },{ 25,185 },{ 345,750 },{ 945,685 },{ 845,655 },    { 880,660 },{ 25,230 },{ 525,1000 },{ 580,1175 },{ 650,1130 },    { 1605,620 },{ 1220,580 },{ 1465,200 },{ 1530,5 },    { 845,680 },{ 725,370 },{ 145,665 },    { 415,635 },{ 510,875 },{ 560,365 },{ 300,465 },{ 520,585 },{ 480,415 },    { 835,625 },{ 975,580 },{ 1215,245 },{ 1320,315 },{ 1250,400 },{ 660,180 },    { 410,250 },{ 420,555 },{ 575,665 },{ 1150,1160 },{ 700,580 },{ 685,595 },    { 685,610 },{ 770,610 },{ 795,645 },{ 720,635 },{ 760,650 },{ 475,960 },    { 95,260 },{ 875,920 },{ 700,500 },{ 555,815 },{ 830,485 },{ 1170,65 },    { 830,610 },{ 605,625 },{ 595,360 },{ 1340,725 },{ 1740,245 }    };
  public static double distance_2city(int[] c1,int[] c2)  //計(jì)算兩點(diǎn)間距離 {    double dis=0;    dis=sqrt(pow(c1[0]-c2[0],2)+pow(c1[1]-c2[1],2));    return dis;  }}

Solution類，存放解，即路徑。

?

package ILSTSP;
import static java.lang.Math.*; import static ILSTSP.City.*;
public class Solution {  public int[] permutation=new int[CITY_SIZE]; //城市排列    public double cost;  
  //獲取隨機(jī)城市排列    public static void random_permutation(int[] cities_permutation){        int n=CITY_SIZE;      int[] temp=new int[CITY_SIZE];      for(int i=0;i         temp[i]=i;       for(int i=0;i      {        int r=(int) random()*n;        cities_permutation[i]=temp[r];        temp[r]=temp[n-1];        n--;      }      cities_permutation[CITY_SIZE-1]=temp[0];    }
    public static double cost_total(int[] cities_permutation){        double total_distance = 0;        for (int i = 0; i < CITY_SIZE-1; i++)            total_distance += distance_2city(cities[cities_permutation[i]], cities[cities_permutation[i+1]]);        total_distance += distance_2city(cities[cities_permutation[0]], cities[cities_permutation[CITY_SIZE-1]]);        //最后一個(gè)城市和第一個(gè)城市計(jì)算距離        return total_distance;    }}

擾動(dòng)類。

?

package ILSTSP;import static ILSTSP.Solution.*; import static ILSTSP.City.*;import static java.lang.Math.*; 
/** * 擾動(dòng)類，跳出局部。 */
public class Perturbation {  //擾動(dòng)  public static void perturbation(Solution best_solution, Solution current_solution){      double_bridge_move(best_solution.permutation, current_solution.permutation);      current_solution.cost = cost_total(current_solution.permutation);  }
  //將城市序列分成4塊，然后按塊重新打亂順序。  //用于擾動(dòng)函數(shù)  private static void double_bridge_move(int[] cities_permutation, int[]  new_cities_permutation){    int[] pos=new int[5];    pos[0]= 0;     pos[1]= (int) random()*(CITY_SIZE/3)+1;      pos[2]= (int) (random()*(CITY_SIZE/3)+CITY_SIZE/3);      pos[3]= (int) (random()*(CITY_SIZE/3)+(CITY_SIZE/3)*2);    pos[4]= CITY_SIZE;
    int n=4;    int[] random_order=new int[4],        temp=new int[4];    for(int i=0;i<4;i++)       temp[i]=i;     for(int i=0;i<3;i++)    {      int r=(int) (random()*n);      random_order[i]=temp[r];      temp[r]=temp[n-1];      n--;    }    random_order[3]=temp[0];
    int deadprotect=0;    do    {      int i=0;        for (int j=pos[random_order[0]];j        {          new_cities_permutation[i]=cities_permutation[j];          i++;        }
        for (int j=pos[random_order[1]];j        {          new_cities_permutation[i]=cities_permutation[j];          i++;        }
        for (int j=pos[random_order[2]];j         {          new_cities_permutation[i]=cities_permutation[j];          i++;        }
        for (int j=pos[random_order[3]];j        {          new_cities_permutation[i]=cities_permutation[j];          i++;        }
        deadprotect++;        if (deadprotect==5) break;
      }while(AcceptanceCriterion(cities_permutation, new_cities_permutation));  }
  //判斷接受準(zhǔn)則  private static boolean AcceptanceCriterion(int[] cities_permutation, int[] new_cities_permutation){    int AcceptLimite=500;    double c1=cost_total(cities_permutation);    double c2=cost_total(new_cities_permutation)-AcceptLimite;    if (c1      return false;    else      return true;  }}

局部搜索類

?

package ILSTSP;
import static ILSTSP.City.*; import static ILSTSP.Delta.*;
/** * 局部搜索類， * 采用反轉(zhuǎn)i到j(luò)之間城市序列的鄰域動(dòng)作。 */
public class LocalSearch {  //本地局部搜索，邊界條件 max_no_improve  //best_solution最優(yōu)解  //current_solution當(dāng)前解  public static void local_search(Solution best_solution, int max_no_improve){      int count = 0;      int i, k;      double inital_cost = best_solution.cost; //初始花費(fèi)      double now_cost = 0;
      Solution current_solution = new Solution(); //為了防止爆棧……直接new了，你懂的
    for (i = 0; i < CITY_SIZE - 1; i++)      {         for (k = i + 1; k < CITY_SIZE; k++)      {        Delta[i][k] = calc_delta(i, k, best_solution.permutation);          }      }
      do      {          //枚舉排列          for (i = 0; i < CITY_SIZE - 1; i++)          {              for (k = i + 1; k < CITY_SIZE; k++)              {                  //鄰域動(dòng)作                  two_opt_swap(best_solution.permutation, current_solution.permutation, i, k);          now_cost = inital_cost + Delta[i][k];                  current_solution.cost = now_cost;                  if (current_solution.cost < best_solution.cost)                  {                      count = 0; //better cost found, so reset                      for (int j = 0; j < CITY_SIZE; j++)                      {                          best_solution.permutation[j] = current_solution.permutation[j];                      }                      best_solution.cost = current_solution.cost;                      inital_cost = best_solution.cost;                      Update(i, k, best_solution.permutation);                  }              }          }          count++;      } while (count <= max_no_improve);  }
  //鄰域動(dòng)作 反轉(zhuǎn)index_i <-> index_j 間的元素  private static void two_opt_swap(int[] cities_permutation, int[] new_cities_permutation, int index_i, int index_j){      for (int i = 0; i < CITY_SIZE; i++)      {          new_cities_permutation[i] = cities_permutation[i];      }
      swap_element(new_cities_permutation, index_i, index_j);  }
  //顛倒數(shù)組中下標(biāo)begin到end的元素位置  private static void swap_element(int[] p, int begin, int end){      int temp;      while (begin < end)      {          temp = p[begin];          p[begin] = p[end];          p[end] = temp;          begin++;          end--;      }  }}

Delta類，在原推文中未講解，這里略微講解一下。

?

Delta實(shí)際上是對(duì)局部搜索的過(guò)程進(jìn)行去重優(yōu)化。

Delta[i][j]中存放的數(shù)字表示反轉(zhuǎn)i，j間路徑后對(duì)總距離的改變量。（我們之前沒(méi)有定義計(jì)算總距離的方法，因?yàn)檫@次改為了由Delta計(jì)算總距離）

?

對(duì)calc_delta部分，每次翻轉(zhuǎn)以后沒(méi)必要再次重新計(jì)算Delta值，只需要在翻轉(zhuǎn)的頭尾做個(gè)小小處理。

?

比如：

有城市序列?? 1-2-3-4-5 總距離 = d_12 + d_23 + d_34 + d_45 + d_51 = A

翻轉(zhuǎn)后的序列 1-4-3-2-5 總距離 = d_14 + d_43 + d_32 + d_25 + d_51 = B

由于 d_ij 與 d_ji是一樣的，所以B也可以表示成 B = A – d_12 – d_45 + d_14 + d_25。

?

對(duì)Update部分，每次翻轉(zhuǎn)以后不需要依次更新所有Delta值，有一部分是可以忽略的。

?

比如：

對(duì)于城市序列1-2-3-4-5-6-7-8-9-10，如果對(duì)3-5應(yīng)用了鄰域操作2-opt （即翻轉(zhuǎn)）， 事實(shí)上對(duì)于Delta 1-2、6-10是不需要重復(fù)計(jì)算的。

package ILSTSP;
import static ILSTSP.City.*;
/** * Delta類， * 對(duì)局部搜索的過(guò)程進(jìn)行去重優(yōu)化。 * Delta[i][j]數(shù)組中存放的數(shù)字表示反轉(zhuǎn)i，j間路徑后對(duì)總距離的改變量。 */
public class Delta {  public static double[][] Delta=new double[CITY_SIZE][CITY_SIZE];
  public static double calc_delta(int i, int k,  int[] tmp)  {    /*                  以下計(jì)算說(shuō)明：                  對(duì)于每個(gè)方案，翻轉(zhuǎn)以后沒(méi)必要再次重新計(jì)算總距離                  只需要在翻轉(zhuǎn)的頭尾做個(gè)小小處理
                  比如：                  有城市序列   1-2-3-4-5 總距離 = d12 + d23 + d34 + d45 + d51 = A                  翻轉(zhuǎn)后的序列 1-4-3-2-5 總距離 = d14 + d43 + d32 + d25 + d51 = B                  由于 dij 與 dji是一樣的，所以B也可以表示成 B = A - d12 - d45 + d14 + d25                  下面的優(yōu)化就是基于這種原理      */      double delta=0;      if((i==0)&&(k==CITY_SIZE-1))         delta=0;      else    {      int i2=(i-1+CITY_SIZE)%CITY_SIZE;          int k2=(k+1)%CITY_SIZE;           delta = 0                  - distance_2city(cities[tmp[i2]], cities[tmp[i]])                  + distance_2city(cities[tmp[i2]], cities[tmp[k]])                  - distance_2city(cities[tmp[k]], cities[tmp[k2]])                  + distance_2city(cities[tmp[i]], cities[tmp[k2]]);    }      return delta;  }
  public static void Update(int i, int k,  int[] tmp)  {    /*    去重處理，對(duì)于Delta數(shù)組來(lái)說(shuō)，對(duì)于城市序列1-2-3-4-5-6-7-8-9-10，如果對(duì)3-5應(yīng)用了鄰域操作2-opt ， 事實(shí)上對(duì)于    7-10之間的翻轉(zhuǎn)是不需要重復(fù)計(jì)算的。所以用Delta提前預(yù)處理一下。
    當(dāng)然由于這里的計(jì)算本身是O（1） 的，事實(shí)上并沒(méi)有帶來(lái)時(shí)間復(fù)雜度的減少（更新操作反而增加了復(fù)雜度）     如果delta計(jì)算 是O（n）的，這種去重操作效果是明顯的。      */    if (i!=0 && k != CITY_SIZE - 1){      i --; k ++;      for (int j = i; j <= k; j ++){        for (int l = j + 1; l < CITY_SIZE; l ++){          Delta[j][l] = calc_delta(j, l, tmp);        }      }
      for (int j = 0; j < k; j ++){        for (int l = i; l <= k; l ++){          if (j >= l) continue;          Delta[j][l] = calc_delta(j, l, tmp);        }      }    }// 如果不是邊界，更新(i-1, k + 1)之間的     else{      for (i = 0; i < CITY_SIZE - 1; i++)          {            for (k = i + 1; k < CITY_SIZE; k++)        {          Delta[i][k] = calc_delta(i, k, tmp);            }          }      }// 邊界要特殊更新 ??}}

欲下載本文相關(guān)的代碼及算例，請(qǐng)移步留言區(qū)。

?贊 賞?

---The End---