轉 | 模擬退火算法（SA）和迭代局部搜索（ILS）求解TSP的Java代碼分享

以下文章來源於數據魔術師 ，作者周航

前言

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

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

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

代碼是基於我自己寫的版本，但我是學習了公眾號推文之后寫的，同時有參照原文代碼，因此雖然與C++代碼有些許區別，但總體是類似的。

本文不再贅述TSP或者SA，ILS了，有需要請點擊下方鏈接學習（一看就懂的那種哦！）：

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

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

不多說了，開始看代碼吧~！

代碼下載

欲下載本文相關的代碼及算例，請關注公眾號【程序猿聲】，后台回復【SAILSJAVA】不包括【】即可

SA求解TSP的JAVA代碼

SA分為四個類：MainRun，Data，Path，SimulatedAnnealing。

MainRun是程序的入口。

package SA;

/**
 * 主函數運行類
 */

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("程序耗時："+usedTime+"s");
  }
}

Data是數據類，存放SA和TSP的數據。

package SA;

/**
 * 數據類，包括：
 * TSP城市坐標，采用柏林52城
 * SA系數。
 */


public class Data {
  public static final double  T0=50000.0, // 初始溫度
  T_min=1e-8,  // 結束溫度
  q=0.98,   // 退火系數
  K=1;  //公式中的常數K 
  public static final int L=1000,  // 每個溫度時的迭代次數，即鏈長
  N=52;  // 城市數量

  public static double[][] city_pos= //柏林52城城市坐標，最優解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+
        ",降溫系數q="+q+",每個溫度迭代"+L+"次");
  }

}

Path是路徑類，打包處理路徑的靜態方法。

package SA;
import static SA.Data.*;
import static java.lang.Math.*;

/**
 * 路徑類，打包處理路徑的靜態方法：
 * 計算兩點間距離 distance
 * 計算路徑長度 path_len
 * 產生新解（新路徑）create_new
 */

public class Path {
  public static double distance(double[] p1,double[] p2)  //計算兩點間距離 
{
    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) //計算路徑長度 
{
    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)   //產生新解
{
     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開始模擬退火。

package SA;
import static SA.Data.*;   
import static SA.Path.*;
import static java.lang.Math.*;
import java.util.*;

/**
 * 模擬退火過程
 */

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<N;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<L;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))<r)
           {
             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() //打印最優解
{
    System.out.println("共降溫"+count+"次，得到的TSP最優距離為"+f2+"路徑為");
     for(int j=0;j<N;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。

主函數運行類，包括了ILS總方法。

package ILSTSP;
import static ILSTSP.City.*; 
import static ILSTSP.Solution.*; 
import static ILSTSP.Perturbation.*;
import static ILSTSP.LocalSearch.*;

/**
 * 主函數運行類，包括了ILS總方法以及計時功能。
 */

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("搜索完成！最優路線總長度 = "+best_solution.cost);
    System.out.println("最優訪問城市序列如下：" );
      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("程序耗時："+usedTime+"s");
  }

  static void iterated_local_search(Solution best_solution, int max_iterations, int max_no_improve)
  {
      Solution current_solution = new Solution();

      //獲得初始隨機解
      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); //擾動+判斷是否接受新解
          local_search(current_solution, max_no_improve);//繼續局部搜索

          //找到更優解
          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" +
          "最優解 = "+ best_solution.cost+
          "當前解 = "+ current_solution.cost 
              );
      }  
  }
}

City類存放城市數據，柏林52城坐標。

package ILSTSP;

import static java.lang.Math.*;

/**
 * TSP數據類，
 * 存放柏林52城坐標，
 * 兩點間距離計算distance_2city。
 */

public class City {

  public static int CITY_SIZE=52;  // 城市數量

  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)  //計算兩點間距離 
{
    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;  

  //獲取隨機城市排列
    public static void random_permutation(int[] cities_permutation)
{
        int n=CITY_SIZE;
      int[] temp=new int[CITY_SIZE];
      for(int i=0;i<CITY_SIZE;i++)
         temp[i]=i; 
      for(int i=0;i<CITY_SIZE-1;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]]);
        //最后一個城市和第一個城市計算距離
        return total_distance;
    }
}

擾動類。

package ILSTSP;
import static ILSTSP.Solution.*; 
import static ILSTSP.City.*;
import static java.lang.Math.*; 

/**
 * 擾動類，跳出局部。
 */

public class Perturbation {
  //擾動
  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塊，然后按塊重新打亂順序。
  //用於擾動函數
  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<pos[random_order[0]+1];j++)
        {
          new_cities_permutation[i]=cities_permutation[j];
          i++;
        }

        for (int j=pos[random_order[1]];j<pos[random_order[1]+1];j++)
        {
          new_cities_permutation[i]=cities_permutation[j];
          i++;
        }

        for (int j=pos[random_order[2]];j<pos[random_order[2]+1];j++)
         {
          new_cities_permutation[i]=cities_permutation[j];
          i++;
        }

        for (int j=pos[random_order[3]];j<pos[random_order[3]+1];j++)
        {
          new_cities_permutation[i]=cities_permutation[j];
          i++;
        }

        deadprotect++;
        if (deadprotect==5) break;

      }while(AcceptanceCriterion(cities_permutation, new_cities_permutation));
  }

  //判斷接受准則
  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<c2)
      return false;
    else
      return true;
  }
}

局部搜索類

package ILSTSP;

import static ILSTSP.City.*; 
import static ILSTSP.Delta.*;

/**
 * 局部搜索類，
 * 采用反轉i到j之間城市序列的鄰域動作。
 */

public class LocalSearch {
  //本地局部搜索，邊界條件 max_no_improve
  //best_solution最優解
  //current_solution當前解
  public static void local_search(Solution best_solution, int max_no_improve)
{
      int count = 0;
      int i, k;
      double inital_cost = best_solution.cost; //初始花費
      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++)
              {
                  //鄰域動作
                  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);
  }

  //鄰域動作 反轉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);
  }

  //顛倒數組中下標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實際上是對局部搜索的過程進行去重優化。

Delta[i][j]中存放的數字表示反轉i，j間路徑后對總距離的改變量。（我們之前沒有定義計算總距離的方法，因為這次改為了由Delta計算總距離）

對calc_delta部分，每次翻轉以后沒必要再次重新計算Delta值，只需要在翻轉的頭尾做個小小處理。

比如：

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

翻轉后的序列 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。

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

比如：

對於城市序列1-2-3-4-5-6-7-8-9-10，如果對3-5應用了鄰域操作2-opt （即翻轉）， 事實上對於Delta 1-2、6-10是不需要重復計算的。

package ILSTSP;

import static ILSTSP.City.*;

/**
 * Delta類，
 * 對局部搜索的過程進行去重優化。
 * Delta[i][j]數組中存放的數字表示反轉i，j間路徑后對總距離的改變量。
 */

public class Delta {
  public static double[][] Delta=new double[CITY_SIZE][CITY_SIZE];

  public static double calc_delta(int i, int k,  int[] tmp)
  {
    /*
                  以下計算說明：
                  對於每個方案，翻轉以后沒必要再次重新計算總距離
                  只需要在翻轉的頭尾做個小小處理

                  比如：
                  有城市序列   1-2-3-4-5 總距離 = d12 + d23 + d34 + d45 + d51 = A
                  翻轉后的序列 1-4-3-2-5 總距離 = d14 + d43 + d32 + d25 + d51 = B
                  由於 dij 與 dji是一樣的，所以B也可以表示成 B = A - d12 - d45 + d14 + d25
                  下面的優化就是基於這種原理
      */
      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)
  {
    /*
    去重處理，對於Delta數組來說，對於城市序列1-2-3-4-5-6-7-8-9-10，如果對3-5應用了鄰域操作2-opt ， 事實上對於
    7-10之間的翻轉是不需要重復計算的。所以用Delta提前預處理一下。

    當然由於這里的計算本身是O（1） 的，事實上並沒有帶來時間復雜度的減少（更新操作反而增加了復雜度） 
    如果delta計算 是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);
            }
          }  
    }// 邊界要特殊更新 
  }
}

這次的代碼是由剛學java的小白小編寫的，翻譯的過程中可能有一些處理不是很好，包括對類的分類等等。歡迎各位在留言區指正交流，小編虛心接受！

好了，這次的java代碼分享就到這里了。兩段代碼，容量不小。同鞋們可以分兩次看。
那么我們下次再見~