遺傳算法的C語言實現(二)
上一次我們使用遺傳算法求解了一個較為復雜的多元非線性函數的極值問題,也基本了解了遺傳算法的實現基本步驟。這一次,我再以經典的TSP問題為例,更加深入地說明遺傳算法中選擇、交叉、變異等核心步驟的實現。而且這一次解決的是離散型問題,上一次解決的是連續型問題,剛好形成對照。
首先介紹一下TSP問題。TSP(traveling salesman problem,旅行商問題)是典型的NP完全問題,即其最壞情況下的時間復雜度隨著問題規模的增大按指數方式增長,到目前為止還沒有找到一個多項式時間的有效算法。TSP問題可以描述為:已知n個城市之間的相互距離,某一旅行商從某一個城市出發,訪問每個城市一次且僅一次,最後回到出發的城市,如何安排才能使其所走的路線最短。換言之,就是尋找一條遍歷n個城市的路徑,或者說搜索自然子集X={1,2,...,n}(X的元素表示對n個城市的編號)的一個排列P(X)={V1,V2,....,Vn},使得Td=∑d(Vi
再來說針對TSP問題使用遺傳算法的步驟。
(1)編碼問題:由於這是一個離散型的問題,我們采用整數編碼的方式,用1~n來表示n個城市,1~n的任意一個排列就構成了問題的一個解。可以知道,對於n個城市的TSP問題,一共有n!種不同的路線。
(2)種群初始化:對於N個個體的種群,隨機給出N個問題的解(相當於是染色體)作為初始種群。這裏具體采用的方法是:1,2,...,n作為第一個個體,然後2,3,..n分別與1交換位置得到n-1個解,從2開始,3,4,...,n分別與2交換位置得到n-2個解,依次類推。(如果這樣還不夠初始種群的數量,可以再考慮n,n-1,...,1這個序列,然後再按照相同的方法生成,等等)
(3)適應度函數:設一個解遍歷初始行走的總距離為D,則適應度fitness=1/D.即總距離越高,適應度越低,總距離越低(解越好),適應度越高。
(4) 選擇操作:個體被選中的概率與適應度成正比,適應度越高,個體被選中的概率越大。這裏仍然采用輪盤賭法。
(5) 交叉操作:交叉操作是遺傳算法最重要的操作,是產生新個體的主要來源,直接關系到算法的全局尋優能力,這裏采用部分映射交叉。比如對於n=10的情況,對於兩個路徑: 1 2 4 5 6 3 9 10 8 7
3 9 7 6 8 10 5 1 2 4
隨機產生兩個[1,10]之間的隨機數r1,r2,代表選擇交叉的位置,比如r1=2,r2=4,如上圖標紅的位置,將第一個個體r1到r2之間的基因(即城市序號)與第二個個體r1到r2之間的基因交換,交換之後變為:
1 9 7 6 6 3 9 10 8 7
3 2 4 5 8 10 5 1 2 4
黃色部分表示交叉過來的基因,這個時候會發現可能交叉過來的基因與原來其他位置上的基因有重復,容易直到,第一個個體重復基因的數目與第二個個體重復基因的數目是相同的(這裏都是3個),只需要把第一個個體重復的基因(用綠色標識)與第二個個體重復的基因做交換,即可以消除沖突。消除沖突之後的解如下:
1 9 7 6 5 3 2 10 8 4
3 2 4 5 8 10 6 1 9 7
(6)變異操作:變異操作采取對於一個染色體(即個體)隨機選取兩個基因進行交換的策略。比如對於染色體:
2 4 6 10 3 1 9 7 8 5
隨機選取了兩個位置p1=3,p2=8(標紅位置),交換這兩個位置的基因,得到新的染色體為:
2 4 7 10 3 1 9 6 8 5
(7) 進化逆轉操作:這個是標準的遺傳算法沒有的,是我們為了加速進化而加入的一個操作。這裏的進化是指逆轉操作具有單向性,即只有逆轉之後個體變得更優才會執行逆轉操作,否則逆轉無效。具體的方法是,隨機產生[1,10](這裏仍然以10個城市為例)之間的兩個隨機數r1和r2(其實也是允許相同的,只是r1,r2相同之後,逆轉自然無效,設置交叉變異都是無效的,但是這不會經常發生),然後將r1和r2之間的基因進行反向排序。比如對於染色體:
1 3 4 2 10 9 8 7 6 5
r1=3,r2=5,它們之間的基因(標紅位置)反向排列之後得到的染色體如下:
1 3 10 2 4 9 8 7 6 5
根據以上的步驟,我們就可以比較容易寫出用遺傳算法求解TSP問題的具體代碼了,這裏仍然使用C語言。先以規模比較小的城市為例,這裏取14個,城市之間的距離會直接在代碼中給出。代碼如下:
/* *遺傳算法(GA) 解決TSP 問題 *案例參考自《MATLAB 智能算法30個案例分析》 *本例以14個城市為例,14個城市的位置坐標如下(括號內第一個元素為X坐標,第二個為縱坐標):1:(16.47,96.10) 2:(16.47,94.44) 3:(20.09,92.54) *4:(22.39,93.37) 5:(25.23,97.24) 6:(22.00,96.05) 7:(20.47,97.02) 8:(17.20,96.29) 9:(16.30,97.38) 10:(14.05,98.12) 11:(16.53,97.38) *12:(21.52,95.59) 13:(19.41,97.13) 14:(20.09,92.55) *遺傳算法實現的步驟為:(1)編碼 (2) 種群初始化 (3) 構造適應度函數 (4) 選擇操作 (5) 交叉操作 (6) 變異操作 (7) 進化逆轉操作 * 具體實現的步驟這裏不詳細說,參考《MATLAB 智能算法30個案例分析》P38 - P40 * update in 16/12/4 * author:Lyrichu * email:[email protected] */ #include<stdio.h> #include<stdlib.h> #include<math.h> #include<time.h> #define maxgen 200 // 最大進化代數 #define sizepop 100 // 種群數目 #define pcross 0.6 // 交叉概率 #define pmutation 0.1 // 變異概率 #define lenchrom 14 // 染色體長度(這裏即為城市個數) double city_pos[lenchrom][2] = {{16.47,96.10},{16.47,94.44},{20.09,92.54},{22.39,93.37},{25.23,97.24},{22.00,96.05},{20.47,97.02}, {17.20,96.29},{16.30,97.38},{14.05,98.12},{16.53,97.38},{21.52,95.59},{19.41,97.13},{20.09,92.55}}; // 定義二維數組存放14個城市的X、Y坐標 int chrom[sizepop][lenchrom]; // 種群 int best_result[lenchrom]; // 最佳路線 double min_distance; // 最短路徑長度 // 函數聲明 void init(void); // 種群初始化函數 double distance(double *,double *); // 計算兩個城市之間的距離 double * min(double *); // 計算距離數組的最小值 double path_len(int *); // 計算某一個方案的路徑長度,適應度函數為路線長度的倒數 void Choice(int [sizepop][lenchrom]); // 選擇操作 void Cross(int [sizepop][lenchrom]); // 交叉操作 void Mutation(int [sizepop][lenchrom]); // 變異操作 void Reverse(int [sizepop][lenchrom]); // 逆轉操作 // 種群初始化 void init(void) { int num = 0; while(num < sizepop) { for(int i=0;i<sizepop;i++) for(int j=0;j<lenchrom;j++) chrom[i][j] = j+1; num++; for(int i=0;i<lenchrom-1;i++) { for(int j=i+1;j<lenchrom;j++) { int temp = chrom[num][i]; chrom[num][i] = chrom[num][j]; chrom[num][j] = temp; // 交換第num個個體的第i個元素和第j個元素 num++; if(num >= sizepop) break; } if(num >= sizepop) break; } // 如果經過上面的循環還是無法產生足夠的初始個體,則隨機再補充一部分 // 具體方式就是選擇兩個基因位置,然後交換 while(num < sizepop) { double r1 = ((double)rand())/(RAND_MAX+1.0); double r2 = ((double)rand())/(RAND_MAX+1.0); int p1 = (int)(lenchrom*r1); // 位置1 int p2 = (int)(lenchrom*r2); // 位置2 int temp = chrom[num][p1]; chrom[num][p1] = chrom[num][p2]; chrom[num][p2] = temp; // 交換基因位置 num++; } } } // 距離函數 double distance(double * city1,double * city2) { double x1 = *city1; double y1 = *(city1+1); double x2 = *(city2); double y2 = *(city2+1); double dis = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); return dis; } // min()函數 double * min(double * arr) { static double best_index[2]; double min_dis = *arr; double min_index = 0; for(int i=1;i<sizepop;i++) { double dis = *(arr+i); if(dis < min_dis) { min_dis = dis; min_index = i; } } best_index[0] = min_index; best_index[1] = min_dis; return best_index; } // 計算路徑長度 double path_len(int * arr) { double path = 0; // 初始化路徑長度 int index = *arr; // 定位到第一個數字(城市序號) for(int i=0;i<lenchrom-1;i++) { int index1 = *(arr+i); int index2 = *(arr+i+1); double dis = distance(city_pos[index1-1],city_pos[index2-1]); path += dis; } int last_index = *(arr+lenchrom-1); // 最後一個城市序號 int first_index = *arr; // 第一個城市序號 double last_dis = distance(city_pos[last_index-1],city_pos[first_index-1]); path = path + last_dis; return path; // 返回總的路徑長度 } // 選擇操作 void Choice(int chrom[sizepop][lenchrom]) { double pick; double choice_arr[sizepop][lenchrom]; double fit_pro[sizepop]; double sum = 0; double fit[sizepop]; // 適應度函數數組(距離的倒數) for(int j=0;j<sizepop;j++) { double path = path_len(chrom[j]); double fitness = 1/path; fit[j] = fitness; sum += fitness; } for(int j=0;j<sizepop;j++) { fit_pro[j] = fit[j]/sum; // 概率數組 } // 開始輪盤賭 for(int i=0;i<sizepop;i++) { pick = ((double)rand())/RAND_MAX; // 0到1之間的隨機數 for(int j=0;j<sizepop;j++) { pick = pick - fit_pro[j]; if(pick<=0) { for(int k=0;k<lenchrom;k++) choice_arr[i][k] = chrom[j][k]; // 選中一個個體 break; } } } for(int i=0;i<sizepop;i++) { for(int j=0;j<lenchrom;j++) chrom[i][j] = choice_arr[i][j]; } } //交叉操作 void Cross(int chrom[sizepop][lenchrom]) { double pick; double pick1,pick2; int choice1,choice2; int pos1,pos2; int temp; int conflict1[lenchrom]; // 沖突位置 int conflict2[lenchrom]; int num1,num2; int index1,index2; int move = 0; // 當前移動的位置 while(move<lenchrom-1) { pick = ((double)rand())/RAND_MAX; // 用於決定是否進行交叉操作 if(pick > pcross) { move += 2; continue; // 本次不進行交叉 } // 采用部分映射雜交 choice1 = move; // 用於選取雜交的兩個父代 choice2 = move+1; // 註意避免下標越界 pick1 = ((double)rand())/(RAND_MAX+1.0); pick2 = ((double)rand())/(RAND_MAX+1.0); pos1 = (int)(pick1*lenchrom); // 用於確定兩個雜交點的位置 pos2 = (int)(pick2*lenchrom); while(pos1 > lenchrom -2 || pos1 < 1) { pick1 = ((double)rand())/(RAND_MAX+1.0); pos1 = (int)(pick1*lenchrom); } while(pos2 > lenchrom -2 || pos2 < 1) { pick2 = ((double)rand())/(RAND_MAX+1.0); pos2 = (int)(pick2*lenchrom); } if(pos1 > pos2) { temp = pos1; pos1 = pos2; pos2 = temp; // 交換pos1和pos2的位置 } for(int j=pos1;j<=pos2;j++) { temp = chrom[choice1][j]; chrom[choice1][j] = chrom[choice2][j]; chrom[choice2][j] = temp; // 逐個交換順序 } num1 = 0; num2 = 0; if(pos1 > 0 && pos2 < lenchrom-1) { for(int j =0;j<=pos1-1;j++) { for(int k=pos1;k<=pos2;k++) { if(chrom[choice1][j] == chrom[choice1][k]) { conflict1[num1] = j; num1++; } if(chrom[choice2][j] == chrom[choice2][k]) { conflict2[num2] = j; num2++; } } } for(int j=pos2+1;j<lenchrom;j++) { for(int k=pos1;k<=pos2;k++) { if(chrom[choice1][j] == chrom[choice1][k]) { conflict1[num1] = j; num1++; } if(chrom[choice2][j] == chrom[choice2][k]) { conflict2[num2] = j; num2++; } } } } if((num1 == num2) && num1 > 0) { for(int j=0;j<num1;j++) { index1 = conflict1[j]; index2 = conflict2[j]; temp = chrom[choice1][index1]; // 交換沖突的位置 chrom[choice1][index1] = chrom[choice2][index2]; chrom[choice2][index2] = temp; } } move += 2; } } // 變異操作 // 變異策略采取隨機選取兩個點,將其對換位置 void Mutation(int chrom[sizepop][lenchrom]) { double pick,pick1,pick2; int pos1,pos2,temp; for(int i=0;i<sizepop;i++) { pick = ((double)rand())/RAND_MAX; // 用於判斷是否進行變異操作 if(pick > pmutation) continue; pick1 = ((double)rand())/(RAND_MAX+1.0); pick2 = ((double)rand())/(RAND_MAX+1.0); pos1 = (int)(pick1*lenchrom); // 選取進行變異的位置 pos2 = (int)(pick2*lenchrom); while(pos1 > lenchrom-1) { pick1 = ((double)rand())/(RAND_MAX+1.0); pos1 = (int)(pick1*lenchrom); } while(pos2 > lenchrom-1) { pick2 = ((double)rand())/(RAND_MAX+1.0); pos2 = (int)(pick2*lenchrom); } temp = chrom[i][pos1]; chrom[i][pos1] = chrom[i][pos2]; chrom[i][pos2] = temp; } } // 進化逆轉操作 void Reverse(int chrom[sizepop][lenchrom]) { double pick1,pick2; double dis,reverse_dis; int n; int flag,pos1,pos2,temp; int reverse_arr[lenchrom]; for(int i=0;i<sizepop;i++) { flag = 0; // 用於控制本次逆轉是否有效 while(flag == 0) { pick1 = ((double)rand())/(RAND_MAX+1.0); pick2 = ((double)rand())/(RAND_MAX+1.0); pos1 = (int)(pick1*lenchrom); // 選取進行逆轉操作的位置 pos2 = (int)(pick2*lenchrom); while(pos1 > lenchrom-1) { pick1 = ((double)rand())/(RAND_MAX+1.0); pos1 = (int)(pick1*lenchrom); } while(pos2 > lenchrom -1) { pick2 = ((double)rand())/(RAND_MAX+1.0); pos2 = (int)(pick2*lenchrom); } if(pos1 > pos2) { temp = pos1; pos1 = pos2; pos2 = temp; // 交換使得pos1 <= pos2 } if(pos1 < pos2) { for(int j=0;j<lenchrom;j++) reverse_arr[j] = chrom[i][j]; // 復制數組 n = 0; // 逆轉數目 for(int j=pos1;j<=pos2;j++) { reverse_arr[j] = chrom[i][pos2-n]; // 逆轉數組 n++; } reverse_dis = path_len(reverse_arr); // 逆轉之後的距離 dis = path_len(chrom[i]); // 原始距離 if(reverse_dis < dis) { for(int j=0;j<lenchrom;j++) chrom[i][j] = reverse_arr[j]; // 更新個體 } } flag = 1; } } } // 主函數 int main(void) { time_t start,finish; start = clock(); // 開始計時 srand((unsigned)time(NULL)); // 初始化隨機數種子 init(); // 初始化種群 int best_fit_index = 0; //最短路徑出現代數 double distance_arr[sizepop]; double dis; for(int j=0;j<sizepop;j++) { dis = path_len(chrom[j]); distance_arr[j] = dis; } double * best_index = min(distance_arr); // 計算最短路徑及序號 min_distance = *(best_index+1); // 最短路徑 int index = (int)(*best_index); // 最短路徑序號 for(int j=0;j<lenchrom;j++) best_result[j] = chrom[index][j]; // 最短路徑序列 // 開始進化 double * new_arr; double new_min_dis; int new_index; for(int i=0;i<maxgen;i++) { Choice(chrom); // 選擇 Cross(chrom); //交叉 Mutation(chrom); //變異 Reverse(chrom); // 逆轉操作 for(int j=0;j<sizepop;j++) distance_arr[j] = path_len(chrom[j]); // 距離數組 new_arr = min(distance_arr); new_min_dis = *(new_arr+1); //新的最短路徑 if(new_min_dis < min_distance) { min_distance = new_min_dis; // 更新最短路徑 new_index =(int)(*new_arr); for(int j=0;j<lenchrom;j++) best_result[j] = chrom[new_index][j]; // 更新最短路徑序列 best_fit_index = i+1; // 最短路徑代數 } } finish = clock(); // 計算結束 double duration = ((double)(finish-start))/CLOCKS_PER_SEC; // 計算耗時 printf("本程序使用遺傳算法求解規模為%d的TSP問題,種群數目為:%d,進化代數為:%d\n",lenchrom,sizepop,maxgen); printf("得到最短路徑為:%d-->%d-->%d-->%d-->%d-->%d-->%d-->%d-->%d-->%d-->%d-->%d-->%d-->%d\n",best_result[0],best_result[1],best_result[2], best_result[3],best_result[4],best_result[5],best_result[6],best_result[7],best_result[8],best_result[9],best_result[10],best_result[11], best_result[12],best_result[13]); printf("最短路徑長度為:%lf,得到最短路徑在第%d代.\n",min_distance,best_fit_index); printf("程序耗時:%lf秒.\n",duration); return 0; }
這裏取種群數目為100,最大進化代數為200,在ubuntu16.04下使用gcc編譯器得到結果如下:
經過多次求解發現,在到了一定代數之後最優解保持不變,而且多次求解得到的最優路徑相同,可以認為求得了問題的實際最優解。但是這裏城市個數只有14個,屬於規模較小的TSP問題,我決定再將問題規模變大來測試遺傳算法優化的性能。這裏選用規模為31的中國TSP問題,首先保持種群數目和進化代數不變,得到結果如下:
可以發現遺傳算法得到的最短路徑大概為16136.633514,而已知的中國TSP問題的最優解為15377,還是有一定的差距,並且算法有的時候收斂過早,陷入了局部最優解,並且穩定性較差,每次運行得到的結果波動較大。為了解決這個問題,我們決定將種群數目和進化代數增大,種群數目變為500,進化代數變為1000,重新運行程序得到的結果如下:
多次運行程序給出的最優解相近,波動很小,其中最好的一次為上圖中的15380.515324,與真實最優解15377十分接近,這說明在增大種群規模以及進化代數之後遺傳算法的優化能力又得到了進一步地提高,而且不易陷入局部最優解,總是能找到接近最優解的次優解,這也充分說明了遺傳算法對於求解TSP問題還是十分有效的,也說明了遺傳算法的普適性。
當然,遺傳算法也有其局限性,比如對於大規模的尋優問題容易早熟,陷入局部最優等等,如果有機會我後面還會再補充這方面的內容,以及遺傳算法在其他領域的應用。
遺傳算法的C語言實現(二)