1. 程式人生 > >Java版GA_TSP(我的第一個Java程序)

Java版GA_TSP(我的第一個Java程序)

結果 end figure 數列 fit 部分 遺傳 ret 平衡

  嗯哼,第一次寫博客,準確說是第一次通過文字的方式記錄自己的工作,閑話少敘,技術汪的博客就該直奔技術主題(關於排版問題,會在不斷寫博客的過程中慢慢學習,先將就著用吧,重在技術嘛~~~)。

  遺傳算法(Genetic Algorithm, GA),作為很多人接觸智能優化算法的第一個算法,互聯網的資料不可謂不多,但不是本文的重點,故在此不再過細展開,只簡單說下大概思想:根據現代生物學理論 “物競天擇,適者生存” 原理,不斷淘汰適應能力差的個體,模擬生物進化過程。大致步驟為:

  1. 生成一個初始種群(Initial Population);
  2. 通過計算種群中各個個體(Individual)的適應度(Fitness)大小來表示種群中各個個體的對環境的適應能力;
  3. 根據 “適者生存” 原則選擇(Select)部分個體繁衍(即Cross、Mutation操作)生成子代種群;
  4. 判斷是否滿足進化結束條件(即算法終止條件). 若滿足,結束進化,輸出結果;否則,重復執行步驟2~3.

  [註]本文用於求解TSP的遺傳算法並非經典遺傳算法,而是針對TSP特征寫的一個簡化版的遺傳算法,求解質量比較low,畢竟是第一個java程序,還是要預留寫改進空間麽不是~~~

  旅行商問題(Travel Salesman Problem, TSP)作為運籌學分支組合優化中的經典問題,一直備受關註。大白話的說法是:一個旅行商從某一城市出發,希望以最短的路程完成對N個城市的巡回訪問,並最終回到出發城市。專業一點的說法是:在有向圖中,從某一點開始,以最小代價補充不漏的遍歷所有點,並返回初始點。

  關於Java必須多說兩句。首先,Java語言是新學的,之前一直用的MATLAB,然而實驗室的工作站帶MATLAB還可以,自己的破筆記本實在帶不動。感覺每次用MATLAB運行代碼,我的小本都會經歷一場生死,從點擊“Run” 360小球噌的一下竄上95%那一瞬間,仿佛全世界就是剩下了小本呼哧呼哧的掙紮。

  本著人道主義精神,終於在某個深夜卸載了MATLAB,頓時感覺小本的重量都輕了許多(可能是幻覺吧)。MATLAB是卸載了,但畢業還要寫代碼啊,總不能不畢業吧,所以之後選用Python試了試,相比動輒十幾個G的MATLAB,幾十M的Python用起來確實不錯,在加上Python語言本身確實簡單,對代碼汪想法的快速程序表達非常友好,以及那句“人生苦短,我用Python”,使得我一度認為Python應該就是我的最終選擇了,但是……但是……但是用Python跑了GA&TSP後我就動搖了,運行速度著實堪憂,可能是自己技術太渣,即使使用Pypy也不能達到我的預期,而作為對算法實現效率有著苛刻要求的組合優化問題,為了把其他論文中的求解結果給PK下去,所以……所以……所以我又拋棄了Python。

  之所以沒有選用高效率的且大一就學過的C語言,則是因為實在搞不定C語言的內存管理,對於指針的使用更是隨心而動,致使最後Debug的時候就像二哈咬刺猬,都不知道從哪下口……

  嗯哼,然後就是Java了,選擇Java就是在使用了排除法後的選擇,兼顧了我筆記本的性能缺陷與運行效率。說了這麽多,終於到達戰場,幹貨出場:

技術分享圖片

遺傳算法示意圖

  下面上代碼:

  Data類:用來放客戶點的坐標,也可以用I/O流從外部文件中獲取。

 1 package GA;
 2 
 3 public class Data {
 4     public static double[][] XY(){
 5         double [][] xy = new double [][] {
 6             {    1304    ,    2312    }    ,
 7             {    3639    ,    1315    }    ,
 8             {    4177    ,    2244    }    ,
 9             {    3712    ,    1399    }    ,
10             {    3488    ,    1535    }    ,
11             {    3326    ,    1556    }    ,
12             {    3238    ,    1229    }    ,
13             {    4196    ,    1004    }    ,
14             {    4312    ,    790    }    ,
15             {    4386    ,    570    }    ,
16             {    3007    ,    1970    }    ,
17             {    2562    ,    1756    }    ,
18             {    2788    ,    1491    }    ,
19             {    2381    ,    1676    }    ,
20             {    1332    ,    695    }    ,
21             {    3715    ,    1678    }    ,
22             {    3918    ,    2179    }    ,
23             {    4061    ,    2370    }    ,
24             {    3780    ,    2212    }    ,
25             {    3676    ,    2578    }    ,
26             {    4029    ,    2838    }    ,
27             {    4263    ,    2931    }    ,
28             {    3429    ,    1908    }    ,
29             {    3507    ,    2367    }    ,
30             {    3394    ,    2643    }    ,
31             {    3439    ,    3201    }    ,
32             {    2935    ,    3240    }    ,
33             {    3140    ,    3550    }    ,
34             {    2545    ,    2357    }    ,
35             {    2778    ,    2826    }    ,
36             {    2370    ,    2975    }    
37         };
38         return xy;
39     }
40 }

  DistanceMatrix類:用來計算各客戶點間的歐氏距離。原打算也放在GSTSP類中,但考慮到有些問題中距離矩陣是給定的,而不是通過公式計算得到的,就單獨拿出來了。

 1 package GA;
 2 
 3 public class DistanceMatrix {
 4     public static double[][] DistMatrix(double[][] xy){
 5         //計算距離矩陣
 6         int N = xy.length;
 7         double[][] DM = new double[N][N];
 8         for (int i = 0; i < N; i++) {
 9             for (int j = 0; j < N; j++) {
10                 DM[i][j] = Math.hypot(xy[i][0] - xy[j][0],xy[i][1] - xy[j][1]);
11             }
12         }
13         return DM;
14     }
15 }

  

  Pop類:用來放置對Pop的操作方法,包括計算適應度的fit()函數(貌似這個命名不規範,但是我懶呀,我就不改,你能怎樣???)和尋找並記錄最有個體的Max()函數(寫這個函數以及初始化種群的時候非常懷戀MATLAB,MATLAB中自帶的的max()能夠同時獲取最大值的索引,初始化種群就更簡單了,一個randperm()解決所有問題。)

 1 package GA;
 2 
 3 public class Pop {
 4     
 5     public static double fit(int[] Ind) {
 6         //計算適應度函數
 7         double[][] xy = Data.XY();
 8         double[][] DM = DistanceMatrix.DistMatrix(xy);
 9         double dist = 0.0;
10         int s = Ind.length;
11         
12         for (int i0 = 0; i0 < s - 1; i0++) {
13             dist += DM[Ind[i0] - 1][Ind[i0 + 1] - 1];
14         }
15         dist += DM[Ind[Ind.length - 1] - 1][Ind[0] - 1];
16         return dist;
17     }
18     
19     public static double[] Max(double[] Fit) {
20         //尋找最優個體及相應的位置
21         double[] max = new double[2];
22         double maxFit = 0.0;
23         double maxIndex = -1;
24         for (int i = 0; i < Fit.length; i++) {
25             if (Fit[i] > maxFit) {
26                 maxFit = Fit[i];
27                 maxIndex = (double)i;
28             }
29         }
30         max[0] = maxFit;
31         max[1] = maxIndex;
32         return max; 
33     }
34 }

  Sharking類:存放擾動算子,想著以後其他智能算法中也能用得著,就一次性造好輪子,以後就能直接用了,寫這幾個擾動算子的時候也想誇誇MATLAB強大的輪子庫, 如fliplr() 就是很好的輪子,還有對數組的各種切割、拼接神奇操作,讓我只能一邊痛苦的造著輪子,一邊默念“效率第一,Java無敵”……

 1 package GA;
 2 
 3 import java.util.Random;
 4 
 5 public class Sharking {
 6     
 7     public static int[] Swap(int[] S) {
 8         //交換操作
 9         
10         Random rand = new Random();
11         int I = rand.nextInt(S.length);
12         int J = rand.nextInt(S.length);
13         
14         int tmp = S[I];
15         S[I] = S[J];
16         S[J] = tmp;
17         
18         return S;
19     }
20     
21     public static int[] Flip(int[] S) {
22         //翻轉操作
23         
24         int[] S0 = new int [S.length];
25         
26         Random rand = new Random();
27         int tmpI = rand.nextInt(S.length);
28         int tmpJ = tmpI;
29         while(tmpI==tmpJ) {
30             tmpJ = rand.nextInt(S.length);
31         }
32         int I = Math.min(tmpI, tmpJ);
33         int J = Math.max(tmpI, tmpJ);
34     
35         for (int i = 0; i < S0.length;i++) {
36             if (i >= I && i <= J) {
37                 S0[i] = S[I+J-i];
38             }else {
39                 S0[i] = S[i];
40             }
41         }
42         return S0;
43     }
44     
45     public static int[] Insert(int[] S) {
46         //插入操作
47         
48         int[] S0 = new int [S.length];
49         
50         Random rand = new Random();
51         int tmpI = rand.nextInt(S.length);
52         int tmpJ = tmpI;
53         while(tmpI==tmpJ) {
54             tmpJ = rand.nextInt(S.length);
55         }
56         int I = Math.min(tmpI, tmpJ);
57         int J = Math.max(tmpI, tmpJ);
58         
59         for (int i = 0; i < S0.length;i++) {
60             if (i >= I && i < J) {
61                 S0[i] = S[i+1];
62             }else if(i==J){
63                 S0[i] = S[I];
64             }else{
65                 S0[i] = S[i];
66             }
67         }
68         return S0;
69     }
72 }

  GATSP類:本來準備將算法單獨弄一個類,後來……後來……我懶呀、我懶呀…..., 所以現在的這個類又臭又長,既包括算法,也包括其他輪子沒有幹的所有事情,比如說結果輸出。

  1 package GA;
  2 
  3 import java.util.Random;
  4 
  5 public class GATSP {
  6     
  7     public static double[] Cusume(double[] A) {
  8         //適應於輪盤賭的累加器
  9         double[] cus = new double[A.length+1];
 10         
 11         cus[0]  = 0.0;
 12         for (int i = 0; i < A.length; i++) {
 13             cus[i+1] = cus[i] + A[i];
 14         }
 15         return cus;
 16     }
 17     
 18     public static double Sum(double[] Arr) {
 19         //求和
 20         double sum = 0.0;
 21         for (int i = 0;i <Arr.length;i++) {
 22             sum += Arr[i];
 23         }
 24         return sum;
 25     }
 26     
 27     
 28     
 29     public static void main(String[] args) {
 30         
 31         long startTime=System.currentTimeMillis();
 32         
 33         //參數列表
 34         //31城市TSP最優解15377.711
 35         int MaxGen = 500;
 36         int PopSize =200;
 37         double[][] xy = Data.XY();
 38         int N = xy.length;
 39         
 40         double[][] DM = DistanceMatrix.DistMatrix(xy);
 41         int[][] Pop = new int[PopSize][N];
 42         double[] Trace = new double[MaxGen];
 43         Pop nowPop = new Pop();
 44         double bs = 1e10;
 45         int[] BS = new int[N];
 46         
 47         
 48         //生成初始種群
 49         for (int p = 0; p < PopSize; p++) {
 50             for (int j = 0; j < N;j++) {
 51                 Pop[p][j] = j + 1;
 52             }
 53             //隨機生成初始個體
 54             for (int k = 0; k < N;k++) {
 55                 Random rand = new Random();
 56                 int r = rand.nextInt(N);
 57                 int tmp;
 58                 tmp = Pop[p][r];
 59                 Pop[p][r] = Pop[p][k];
 60                 Pop[p][k] = tmp;
 61             }
 62         }
 63         
 64         //進入叠代
 65         for (int gen = 0; gen < MaxGen;gen++) {
 66             // 計算種群適應度
 67             double[] Fit = new double[PopSize];
 68             int[] indiv = new int[N];
 69             
 70             for (int p = 0; p < PopSize;p++) {
 71                 //取一個個體
 72                 for (int j = 0; j < N;j++) {
 73                     indiv[j] = Pop[p][j];
 74                 }
 75                 Fit[p] = nowPop.fit(indiv);
 76             }
 77             
 78             //更新最優個體以及最優個體的適應度
 79             double[] SortAfterFit = new double[PopSize];//拷貝一份適應度數組
 80             for (int i = 0; i < PopSize;i++) {
 81                 SortAfterFit[i] = Fit[i];
 82             }
 83             double[] BestS = nowPop.Max(Fit);
 84             double tmpbs = BestS[0];            //當前代最優解(最優個體的適應度)
 85             if (tmpbs < bs) {
 86                 bs = tmpbs;
 87                 int BestIndex = (int)BestS[1];
 88                 for (int i = 0; i < N; i++) {
 89                     BS[i] =Pop[BestIndex][i];            //最優個體
 90                 }
 91             }
 92             Trace[gen] = bs;
 93             
 94             //輪盤賭選擇
 95             double[] cusFit0 = Cusume(Fit);//數組長度變為N+1,補充了首個元素0.0
 96             double sumFit = Sum(Fit);
 97                 //歸一化
 98             double[] cusFit = new double[cusFit0.length];
 99             for (int i = 0; i < cusFit.length; i++) {
100                 cusFit[i] =cusFit0[i] / sumFit;
101             }
102             
103             int[][] newPop = new int[PopSize][N];
104             for (int q = 0;q < N; q++) {
105                  newPop[0][q] = BS[q];
106              }
107             for (int p = 1; p < PopSize; p++) {
108                 double rand = Math.random();
109                 for (int r = 0; r < PopSize; r++) {
110                     if (rand > cusFit[r] && rand <= cusFit[r+1]) {
111                         for (int q = 0;q < N; q++) {
112                             newPop[p][q] = Pop[r][q];
113                         }
114                     }
115                 }
116             }
117             
118             //擾動操作
119             for (int p = 0; p < PopSize; p++) {
120                 double R = Math.random();
121                 
122                 int[] S = new int[N];
123                 for (int i = 0; i < N; i++) {
124                     S[i] = newPop[p][i];
125                 }
126                 
127                 int[] S0 = new int[N];
128                 if (R < 0.33) {
129                     S0 = Sharking.Swap(S);
130                 }else if (R > 0.67) {
131                     S0 = Sharking.Insert(S);
132                 }else {
133                     S0 = Sharking.Flip(S);
134                 }
135                 
136                 for (int i = 0; i < N; i++) {
137                     newPop[p][i] = S0[i];
138                 }
139             }
140             
141             //更新種群
142             for (int p = 0; p < PopSize; p++) {
143                 for (int i = 0; i < N; i++) {
144                     Pop[p][i] = newPop[p][i];
145                 }
146             }
147         }//結果叠代
148         
149         long endTime=System.currentTimeMillis();
150         //結果輸出
151         System.out.println("經過"+MaxGen+"次叠代,最短路徑長度為:"+bs);
152         System.out.println("程序用時 "+(double)(endTime - startTime)/1000+"秒.");
153         double bs0 = 15377.711;
154         System.out.println("與最優解的誤差為"+(bs-bs0)/bs0*100+"%.");
155         for (int i = 0; i < MaxGen; i++) {
156             System.out.println(i+1+"  "+Trace[i]);
157         }
158         
159         System.out.println("相應的最短路徑為");
160         for (int i = 0; i < N; i++) {
161             System.out.print(BS[i]+"->");
162         }
163         System.out.print(BS[0]);
164     }
165 }

  本來劇情發展到這裏,應該是求解結果展示,然後就可以領盒飯了,無奈結果太渣,不忍心拿出來,所以劇情稍微轉下,多說幾句關於MATLAB與Java:

  本來應該還有一個PlotFigure類,然而Java學的不精,暫時還不會用,說到這裏,又懷戀起MATLAB來,其強大的可視化結果能力,確實很有吸引力,想想之前只要用MATLAB寫上幾行代碼就可以在每次跑完程序有一堆酷酷的圖出來,比如說算法叠代收斂圖、最優路徑圖什麽的,現在只有“28->23->29->19->……”這種結果,心裏實在不是很平衡,所以,我決定……下次一定要把Java的繪圖學會。

Ps:還是截圖紀念下第一次Java寫遺傳算法~~~

技術分享圖片

  [完]

Java版GA_TSP(我的第一個Java程序)