基於粒子群演算法求解TSP問題（JAVA）__C語言

一、TSP問題

TSP問題（Travelling Salesman Problem）即旅行商問題，又譯為旅行推銷員問題、貨郎擔問題，是數學領域中著名問題之一。假設有一個旅行商人要拜訪n個城市，他必須選擇所要走的路徑，路徑的限制是每個城市只能拜訪一次，而且最後要回到原來出發的城市。路徑的選擇目標是要求得的路徑路程為所有路徑之中的最小值。

TSP問題是一個組合最佳化問題。該問題可以被證明具有NPC計算複雜性。TSP問題可以分為兩類，一類是對稱TSP問題（Symmetric TSP），另一類是非對稱問題（Asymmetric TSP）。所有的TSP問題都可以用一個圖（Graph）來描述：

V={c1, c2, …, ci, …, cn}，i = 1,2, …, n，是所有城市的集合.ci表示第i個城市，n為城市的數目；

E={(r, s): r,s∈ V}是所有城市之間串連的集合；

C = {crs: r,s∈ V}是所有城市之間串連的成本度量（一般為城市之間的距離）；

如果crs = csr, 那麼該TSP問題為對稱的，否則為非對稱的。

一個TSP問題可以表達為：

求解遍曆圖G = (V, E, C)，所有的節點一次並且回到起始節點，使得串連這些節點的路徑成本最低。

二、粒子群演算法

1、基本思想

粒子群演算法簡稱PSO，它的基本思想是類比鳥群的捕食行為。設想這樣一個情境：一群鳥在隨機搜尋食物。在這個地區裡只有一塊食物。所有的鳥都不知道食物在那裡。但是他們知道當前的位置離食物還有多遠。那麼找到食物的最優策略是什麼呢。最簡單有效就是搜尋目前離食物最近的鳥的周圍地區。

PSO從這種模型中得到啟示並用於解決最佳化問題。PSO中，每個最佳化問題的解都是搜尋空間中的一隻鳥。我們稱之為“粒子”。所有的粒子都有一個由被最佳化的函數決定的適應值(fitness value)，每個粒子還有一個速度決定他們飛翔的方向和距離。然後粒子們就追隨當前的最優粒子在解空間中搜尋。

PSO 初始化為一群隨機粒子(隨機解)。然後通過迭代找到最優解。在每一次迭代中，粒子通過跟蹤兩個"極值"來更新自己。第一個就是粒子本身所找到的最優解，這個解叫做個體極值pBest。另一個極值是整個種群目前找到的最優解，這個極值是全域極值gBest。另外也可以不用整個種群而只是用其中一部分作為粒子的鄰居，那麼在所有鄰置中的極值就是局部極值。

2、粒子公式

在找到這兩個最優值時，粒子根據如下的公式來更新自己的速度和新的位置：

v[i] = w * v[i] + c1 * rand() * (pbest[i] - present[i]) + c2 * rand() * (gbest - present[i])    

present[i] = present[i] + v[i] 

其中v[i]代表第i個粒子的速度，w代表慣性權值,c1和c2表示學習參數，rand()表示在0-1之間的隨機數,pbest[i]代表第i個粒子搜尋到的最優值,gbest代表整個叢集搜尋到的最優值,present[i]代表第i個粒子的當前位置。

3、個人見解

截止目前為止，粒子群最佳化演算法大體上可以分為兩大類，一類是最初的基本粒子群最佳化演算法，一類是改進後的廣義粒子群最佳化演算法，其實PSO最初的設計主要用於處理連續最佳化問題，如求函數極值，在複雜的組合最佳化問題上它的應用相當有限，後來經過眾多學者的改進才將其應用於求解TSP和單機調度之類的問題。廣義粒子群演算法模型和遺傳演算法相當類似，目前網上有關於粒子群演算法求解TSP的很多論文或代碼都是基於廣義粒子群演算法的，說簡單點就是進化思想，用交叉變異代替了基本粒子群演算法的迭代公式，當然他們也還是有粒子群最佳化的本質思想的，如與全域最優編碼交叉，與局部最優編碼交叉，變異等都是源自於基本粒子群演算法的迭代公式。

三、粒子群最佳化演算法求解TSP問題 關於基本粒子群最佳化演算法的使用，可參考最近這篇文章， 自話粒子群演算法 ，既然是求解TSP，使用基本的粒子群演算法迭代公式肯定是不行的，在這裡我也不想寫一個與遺傳演算法差不多的粒子群演算法，確實沒必要，只要看瞭解遺傳演算法，分分鐘就能寫出來，在這裡我想使用一種，比較特別的迭代方式，先來看看它的迭代公式：

在這裡我不多囉嗦了，想要詳細瞭解，可能你需要看一下這篇論文：http://download.csdn.net/detail/wangqiuyun/6373499 有點古老是吧，沒關係，要的是它的思想。

我們使用TSP問題依然來自於tsplib上的att48，這是一個對稱TSP問題，城市規模為48，其最優值為10628.其距離計算方法下圖所示：

好，下面是具體代碼：

package noah;import java.io.BufferedReader;import java.io.FileInputStream;import java.io.IOException;import java.io.InputStreamReader;import java.util.ArrayList;import java.util.Random;public class PSO {private int bestNum;private float w;private int MAX_GEN;// 迭代次數private int scale;// 種群規模private int cityNum; // 城市數量，編碼長度private int t;// 當前代數private int[][] distance; // 距離矩陣private int[][] oPopulation;// 粒子群private ArrayList<ArrayList<SO>> listV;// 每科粒子的初始交換序列private int[][] Pd;// 一顆粒子曆代中出現最好的解，private int[] vPd;// 解的評價值private int[] Pgd;// 整個粒子群經曆過的的最好的解，每個粒子都能記住自己搜尋到的最好解private int vPgd;// 最好的解的評價值private int bestT;// 最佳出現代數private int[] fitness;// 種群適應度，表示種群中各個個體的適應度private Random random;public PSO() {}/** * constructor of GA *  * @param n *            城市數量 * @param g *            運行代數 * @param w *            權重 **/public PSO(int n, int g, int s, float w) {this.cityNum = n;this.MAX_GEN = g;this.scale = s;this.w = w;}// 給編譯器一條指令，告訴它對被批註的代碼元素內部的某些警告保持靜默@SuppressWarnings("resource")/** * 初始化PSO演算法類 * @param filename 資料檔案名，該檔案儲存體所有城市節點座標資料 * @throws IOException */private void init(String filename) throws IOException {// 讀取資料int[] x;int[] y;String strbuff;BufferedReader data = new BufferedReader(new InputStreamReader(new FileInputStream(filename)));distance = new int[cityNum][cityNum];x = new int[cityNum];y = new int[cityNum];for (int i = 0; i < cityNum; i++) {// 讀取一行資料，資料格式1 6734 1453strbuff = data.readLine();// 字元分割String[] strcol = strbuff.split(" ");x[i] = Integer.valueOf(strcol[1]);// x座標y[i] = Integer.valueOf(strcol[2]);// y座標}// 計算距離矩陣// ，針對具體問題，距離計算方法也不一樣，此處用的是att48作為案例，它有48個城市，距離計算方法為偽歐氏距離，最優值為10628for (int i = 0; i < cityNum - 1; i++) {distance[i][i] = 0; // 對角線為0for (int j = i + 1; j < cityNum; j++) {double rij = Math.sqrt(((x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j])* (y[i] - y[j])) / 10.0);// 四捨五入，取整int tij = (int) Math.round(rij);if (tij < rij) {distance[i][j] = tij + 1;distance[j][i] = distance[i][j];} else {distance[i][j] = tij;distance[j][i] = distance[i][j];}}}distance[cityNum - 1][cityNum - 1] = 0;oPopulation = new int[scale][cityNum];fitness = new int[scale];Pd = new int[scale][cityNum];vPd = new int[scale];/* * for(int i=0;i<scale;i++) { vPd[i]=Integer.MAX_VALUE; } */Pgd = new int[cityNum];vPgd = Integer.MAX_VALUE;// nPopulation = new int[scale][cityNum];bestT = 0;t = 0;random = new Random(System.currentTimeMillis());/* * for(int i=0;i<cityNum;i++) { for(int j=0;j<cityNum;j++) { * System.out.print(distance[i][j]+","); } System.out.println(); } */}// 初始化種群，多種隨機產生辦法void initGroup() {int i, j, k;for (k = 0; k < scale; k++)// 種群數{oPopulation[k][0] = random.nextInt(65535) % cityNum;for (i = 1; i < cityNum;)// 粒子個數{oPopulation[k][i] = random.nextInt(65535) % cityNum;for (j = 0; j < i; j++) {if (oPopulation[k][i] == oPopulation[k][j]) {break;}}if (j == i) {i++;}}}/* * for(i=0;i<scale;i++) { for(j=0;j<cityNum;j++) { * System.out.print(oldPopulation[i][j]+","); } System.out.println(); } */}void initListV() {int ra;int raA;int raB;listV = new ArrayList<ArrayList<SO>>();for (int i = 0; i < scale; i++) {ArrayList<SO> list = new ArrayList<SO>();ra = random.nextInt(65535) % cityNum;for (int j = 0; j < ra; j++) {raA = random.nextInt(65535) % cityNum;raB = random.nextInt(65535) % cityNum;while (raA == raB) {raB = random.nextInt(65535) % cityNum;}// raA與raB不一樣SO s = new SO(raA, raB);list.add(s);}listV.add(list);}}public int evaluate(int[] chr) {// 0123int len = 0;// 編碼，起始城市,城市1,城市2...城市nfor (int i = 1; i < cityNum; i++) {len += distance[chr[i - 1]][chr[i]];}// 城市n,起始城市len += distance[chr[cityNum - 1]][chr[0]];return len;}// 求一個基本交換序列作用於編碼arr後的編碼public void add(int[] arr, ArrayList<SO> list) {int temp = -1;SO s;for (int i = 0; i < list.size(); i++) {s = list.get(i);temp = arr[s.getX()];arr[s.getX()] = arr[s.getY()];arr[s.getY()] = temp;}}// 求兩個編碼的基本交換序列，如A-B=SSpublic ArrayList<SO> minus(int[] a, int[] b) {int[] temp = b.clone();/* * int[] temp=new int[L]; for(int i=0;i<L;i++) { temp[i]=b[i]; } */int index;// 交換子SO s;// 交換序列ArrayList<SO> list = new ArrayList<SO>();for (int i = 0; i < cityNum; i++) {if (a[i] != temp[i]) {// 在temp中找出與a[i]相同數值的下標indexindex = findNum(temp, a[i]);// 在temp中交換下標i與下標index的值changeIndex(temp, i, index);// 記住交換子s = new SO(i, index);// 儲存交換子list.add(s);}}return list;}// 在arr數組中尋找num，返回num的下標public int findNum(int[] arr, int num) {int index = -1;for (int i = 0; i < cityNum; i++) {if (arr[i] == num) {index = i;break;}}return index;}// 將數組arr下標index1與下標index2的值交換public void changeIndex(int[] arr, int index1, int index2) {int temp = arr[index1];arr[index1] = arr[index2];arr[index2] = temp;}// 二維數組拷貝public void copyarray(int[][] from, int[][] to) {for (int i = 0; i < scale; i++) {for (int j = 0; j < cityNum; j++) {to[i][j] = from[i][j];}}}// 一維數組拷貝public void copyarrayNum(int[] from, int[] to) {for (int i = 0; i < cityNum; i++) {to[i] = from[i];}}public void evolution() {int i, j, k;int len = 0;float ra = 0f;ArrayList<SO> Vi;// 迭代一次for (t = 0; t < MAX_GEN; t++) {// 對於每顆粒子for (i = 0; i < scale; i++) {if(i==bestNum) continue;ArrayList<SO> Vii = new ArrayList<SO>();//System.out.println("------------------------------");// 更新速度// Vii=wVi+ra(Pid-Xid)+rb(Pgd-Xid)Vi = listV.get(i);// wVi+表示擷取Vi中size*w取整個交換序列len = (int) (Vi.size() * w);//越界判斷//if(len>cityNum) len=cityNum;//System.out.println("w:"+w+" len:"+len+" Vi.size():"+Vi.size());for (j = 0; j < len; j++) {Vii.add(Vi.get(j));}// Pid-XidArrayList<SO> a = minus(Pd[i], oPopulation[i]);ra = random.nextFloat();// ra(Pid-Xid)+len = (int) (a.size() * ra);//越界判斷//if(len>cityNum) len=cityNum;//System.out.println("ra:"+ra+" len:"+len+" a.size():"+a.size());for (j = 0; j < len; j++) {Vii.add(a.get(j));}// Pid-XidArrayList<SO> b = minus(Pgd, oPopulation[i]);ra = random.nextFloat();// ra(Pid-Xid)+len = (int) (b.size() * ra);//越界判斷//if(len>cityNum) len=cityNum;//System.out.println("ra:"+ra+" len:"+len+" b.size():"+b.size());for (j = 0; j < len; j++) {SO tt= b.get(j);Vii.add(tt);}//System.out.println("------------------------------Vii.size():"+Vii.size());// 儲存新ViilistV.add(i, Vii);// 更新位置// Xid’=Xid+Vidadd(oPopulation[i], Vii);}// 計算新粒子群適應度，Fitness[max],選出最好的解for (k = 0; k < scale; k++) {fitness[k] = evaluate(oPopulation[k]);if (vPd[k] > fitness[k]) {vPd[k] = fitness[k];copyarrayNum(oPopulation[k], Pd[k]);bestNum=k;}if (vPgd > vPd[k]) {System.out.println("最佳長度"+vPgd+" 代數："+bestT);bestT = t;vPgd = vPd[k];copyarrayNum(Pd[k], Pgd);}}}}public void solve() {int i;int k;initGroup();initListV();// 每顆粒子記住自己最好的解copyarray(oPopulation, Pd);// 計算初始化種群適應度，Fitness[max],選出最好的解for (k = 0; k < scale; k++) {fitness[k] = evaluate(oPopulation[k]);vPd[k] = fitness[k];if (vPgd > vPd[k]) {vPgd = vPd[k];copyarrayNum(Pd[k], Pgd);bestNum=k;}}// 列印System.out.println("初始粒子群...");for (k = 0; k < scale; k++) {for (i = 0; i < cityNum; i++) {System.out.print(oPopulation[k][i] + ",");}System.out.println();System.out.println("----" + fitness[k]);/*ArrayList<SO> li = listV.get(k);int l = li.size();for (i = 0; i < l; i++) {li.get(i).print();}System.out.println("----");*/}// 進化evolution();// 列印System.out.println("最後粒子群...");for (k = 0; k < scale; k++) {for (i = 0; i < cityNum; i++) {System.out.print(oPopulation[k][i] + ",");}System.out.println();System.out.println("----" + fitness[k]);/*ArrayList<SO> li = listV.get(k);int l = li.size();for (i = 0; i < l; i++) {li.get(i).print();}System.out.println("----");*/}System.out.println("最佳長度出現代數：");System.out.println(bestT);System.out.println("最佳長度");System.out.println(vPgd);System.out.println("最佳路徑：");for (i = 0; i < cityNum; i++) {System.out.print(Pgd[i] + ",");}}/** * @param args * @throws IOException */public static void main(String[] args) throws IOException {System.out.println("Start....");PSO pso = new PSO(48, 5000, 30, 0.5f);pso.init("c://data.txt");pso.solve();}}

運行結果截圖：

四、總結 對於這個實驗結果，其實我個人是不能接受的，但是目前我只能把問題歸結於那篇論文提出的迭代公式上，不過那論文年代久遠，它上面就試了14個點，這裡是48個，差距肯定是有的，或許這裡再進行一些調整，沒準效果會好些，但是目前我本人時間有限，還來不及深入實驗，只能讓感興趣的人來進行深入剝解。