嗯哼,第一次寫博客,准確說是第一次通過文字的方式記錄自己的工作,閑話少敘,技術汪的博客就該直奔技術主題(關於排版問題,會在不斷寫博客的過程中慢慢學習,先將就着用吧,重在技術嘛~~~)。
遺傳算法(Genetic Algorithm, GA),作為很多人接觸智能優化算法的第一個算法,互聯網上關於遺傳算法的資料不可謂不多,但由於其不是本文的重點,故在此不過細展開,只簡單說下大概思想:根據現代生物學理論 “物競天擇,適者生存” 原理,不斷淘汰適應能力差的個體,模擬生物進化過程。大致步驟為:
- 生成一個初始種群(Initial Population);
- 通過計算種群中各個個體(Individual)的適應度(Fitness)大小來表示種群中各個個體的對環境的適應能力;
- 根據 “適者生存” 原則選擇(Select)部分個體繁衍(即Cross、Mutation操作)生成子代種群;
- 判斷是否滿足進化結束條件(即算法終止條件). 若滿足,結束進化,輸出結果;否則,重復執行步驟2~3.
[注]本文用於求解TSP的遺傳算法並非經典遺傳算法,而是針對TSP特征寫的一個簡化版的遺傳算法,求解質量比較low,畢竟是第一個java程序,還是要預留些改進空間么不是~~~,下一篇文章將針對本文所用的遺傳算法進行改進,給個傳送門(Java版GA_TSP (2))。
旅行商問題(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類:用來計算各客戶點間的歐氏距離。原打算也放在GATSP類中,但考慮到有些問題中距離矩陣是給定的,而不是通過公式計算得到的,就單獨拿出來了。
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 1/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 = 1/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寫遺傳算法~~~
[完]