遺傳算法的C語言實現(一):以非線性函數求極值爲例

      之前搞數學建模的時候，研究過(其實也不算是研究，只是大概瞭解)一些人工智能算法，好比前面已經說過的粒子羣算法(PSO),還有著名的遺傳算法(GA),模擬退火算法(SA),蟻羣算法(ACA）等。當時懂得很是淺，只會copy別人的代碼(通常是MATLAB),改一改值和參數，東拼西湊就拿過來用了，根本沒有搞懂的其內在的原理究竟是什麼。這一段時間，我又從新翻了一下當時買的那本《MATLAB智能算法30個案例分析》，重讀一遍，發現這本書講的仍是很是不錯的，不只有現成算法的MATLAB實現，並且把每一種算法的實現原理都基本講清楚了，我再次看的時候，之前不懂的(也許之前太急功近利了)如今竟然基本所有弄懂了。真是讓我感到很是開心~~~因此我萌發了一個衝動，想用如今本身比較熟悉的C語言，把這本書上的案例從新實現一遍，完全弄懂每一個算法的原理。我打算不按期地更新幾篇來介紹其中幾種常見的智能算法和其應用(固然我也不止會看這本書，還會參考別的書籍和論文),儘可能把每種算法的原理給講明白，也算是對我接觸這些智能算法的一個總結吧。

     今天首先介紹遺傳算法(genetic algorithm,GA)。

     遺傳算法是一種進化算法，其基本原理是模仿天然界中的生物「物競天擇，適者生存」的進化法則，把問題參數編碼爲染色體，再利用迭代的方式進行選擇、交叉、變異等運算法則來交換種羣中染色體的信息，最終生成符合優化目標的染色體。

     在遺傳算法中，染色體對應的是數據或者數組，一般是由一維的串結構數據來表示，串上各個位置對應基因的的取值。基因組成的串就是染色體，或者稱爲基因型個體。必定數量的個體組成了種羣，種羣中個體數目的大小稱爲種羣大小，也稱爲種羣規模，而各個個體對環境的適應程度稱爲適應度。

　　標準遺傳算法的步驟以下：

　　(1) 編碼：遺傳算法在搜索解空間以前須要將解數據表示成遺傳空間的基因型串結構數據，這些串結構數據的不一樣組合構成了不一樣的染色體。

　　(2)初始化：即生成初始種羣。具體的作法是隨機生成N個初始的染色體(解空間的解),每個染色體其實就至關於一個個體，N個個體構成了一個初始種羣。遺傳算法以這N個個體做爲初始值開始進化。

　　(3) 適應度評估:適應度代表個體或者解的優劣性。不一樣的問題，定義適應度函數的方式也不一樣。好比若是是求一個函數的最大值，那麼通常某個解對應的函數值越大，那麼這個個體(解)的適應度也就越高。

     （4）選擇:選擇的目的是爲了從當前種羣中選出優良的個體，使它們有機會成爲父代繁殖下一代子孫。遺傳算法經過選擇過程體現這一思想，進行選擇的原則是適應性強的個體爲下一代貢獻一個或者多個後代的機率大。這體現了達爾文的適者生存原則。

　　(5)交叉：交叉操做是遺傳算法中最主要的遺傳操做。經過交叉能夠獲得新一代個體，新個體組合了其父代的個體特徵。交叉體現了信息交換的思想。

　　(6)變異:變異首先在羣體中隨機選擇一個個體，對於選中的個體以必定的機率（一般是比較小的機率，這與天然界一致，天然界的變異都是小几率事件）隨機改變染色體中某個基因的值。

有的時候除了選擇選擇、交叉、變異這三種操做以外，咱們還會針對具體的問題加入其它的操做（好比逆轉之類），可是選擇、交叉、變異是全部的遺傳算法都共同的擁有的遺傳操做。

本文以遺傳算法常見的一個應用領域——求解複雜的非線性函數極值問題爲例，來講明如何用代碼(這裏是用C語言，固然你能夠換作別的任何一種語言)來具體實現上述的遺傳算法操做。求解的函數來自《MATLAB智能算法30個案例分析》，即：

f = -5*sin(x1)*sin(x2)*sin(x3)*sin(x4)*sin(x5) - sin(5*x1)*sin(5*x2)*sin(5*x3)*sin(5*x4)*sin(5*x5) + 8

其中x1,x2,x3,x4,x5是0~0.9*PI之間的實數。該函數的最小值爲2,當x1,x2,x3,x4,x5都取PI/2時獲得。

使用C語言實現的遺傳算法求解代碼以下：

/* * 遺傳算法求解函數的最優值問題 * 參考自《MATLAB智能算法30個案例分析》 * 本例的尋優函數爲:f = -5*sin(x1)*sin(x2)*sin(x3)*sin(x4)*sin(x5) - sin(5*x1)*sin(5*x2)*sin(5*x3)*sin(5*x4)*sin(5*x5) + 8 * 其中x1,x2,x3,x4,x5是0~0.9*PI之間的實數。該函數的最小值爲2,當x1,x2,x3,x4,x5都取PI/2時獲得。 * update in 16/12/3 * author:Lyrichu * email:919987476@qq.com */ #include<stdio.h> #include<stdlib.h> #include<math.h> #include<time.h>
#define PI 3.1415926 //圓周率
#define sizepop 50 // 種羣數目
#define maxgen 500 // 進化代數
#define pcross 0.6 // 交叉機率
#define pmutation 0.1 // 變異機率
#define lenchrom 5 // 染色體長度
#define bound_down 0 // 變量下界，這裏由於都相同，因此就用單個值去代替了，若是每一個變量上下界不一樣，也許須要用數組定義
#define bound_up (0.9*3.1415926) // 上界
double chrom[sizepop][lenchrom]; // 種羣數組
double fitness[sizepop]; //種羣每一個個體的適應度
double fitness_prob[sizepop]; // 每一個個體被選中的機率(使用輪盤賭法肯定)
double bestfitness[maxgen]; // 每一代最優值
double gbest_pos[lenchrom]; // 取最優值的位置
double average_best[maxgen+1]; // 每一代平均最優值
double gbest; // 全部進化中的最優值
int gbest_index; // 取得最優值的迭代次數序號 // 目標函數
double fit_func(double * arr) { double x1 = *arr; double x2 = *(arr+1); double x3 = *(arr+2); double x4 = *(arr+3); double x5 = *(arr+4); double func_result = -5*sin(x1)*sin(x2)*sin(x3)*sin(x4)*sin(x5) - sin(5*x1)*sin(5*x2)*sin(5*x3)*sin(5*x4)*sin(5*x5) + 8; return func_result; } // 求和函數
double sum(double * fitness) { double sum_fit = 0; for(int i=0;i<sizepop;i++) sum_fit += *(fitness+i); return sum_fit; } // 求最小值函數
double * min(double * fitness) { double min_fit = *fitness; double min_index = 0; static double arr[2]; for(int i=1;i<sizepop;i++) { if(*(fitness+i) < min_fit) { min_fit = *(fitness+i); min_index = i; } } arr[0] = min_index; arr[1] = min_fit; return arr; } // 種羣初始化
void init_chrom() { for(int i=0;i<sizepop;i++) { for(int j=0;j<lenchrom;j++) { chrom[i][j] = bound_up*(((double)rand())/RAND_MAX); } fitness[i] = fit_func(chrom[i]); // 初始化適應度
 } } // 選擇操做
void Select(double chrom[sizepop][lenchrom]) { int index[sizepop]; for(int i=0;i<sizepop;i++) { double * arr = chrom[i]; fitness[i] = 1/(fit_func(arr)); // 本例是求最小值，適應度爲目標函數的倒數，即函數值越小，適應度越大
 } double sum_fitness = 0; for(int i=0;i<sizepop;i++) { sum_fitness += fitness[i]; // 適應度求和
 } for(int i=0;i<sizepop;i++) { fitness_prob[i] = fitness[i]/sum_fitness; } for(int i=0;i<sizepop;i++) { fitness[i] = 1/fitness[i]; // 恢復函數值
 } for(int i=0;i<sizepop;i++) // sizepop 次輪盤賭
 { double pick = ((double)rand())/RAND_MAX; while(pick < 0.0001) pick = ((double)rand())/RAND_MAX; for(int j=0;j<sizepop;j++) { pick = pick - fitness_prob[j]; if(pick<=0) { index[i] = j; break; } } } for(int i=0;i<sizepop;i++) { for(int j=0;j<lenchrom;j++) { chrom[i][j] = chrom[index[i]][j]; } fitness[i] = fitness[index[i]]; } } // 交叉操做
void Cross(double chrom[sizepop][lenchrom]) { for(int i=0;i<sizepop;i++) { // 隨機選擇兩個染色體進行交叉
        double pick1 = ((double)rand())/RAND_MAX; double pick2 = ((double)rand())/RAND_MAX; int choice1 = (int)(pick1*sizepop);// 第一個隨機選取的染色體序號
        int choice2 = (int)(pick2*sizepop);// 第二個染色體序號
        while(choice1 > sizepop-1) { pick1 = ((double)rand())/RAND_MAX; choice1 = (int)(pick1*sizepop); //防止選取位置過界
 } while(choice2 > sizepop-1) { pick2 = ((double)rand())/RAND_MAX; choice2 = (int)(pick2*sizepop); } double pick = ((double)rand())/RAND_MAX;// 用於判斷是否進行交叉操做
        if(pick>pcross) continue; int flag = 0; // 判斷交叉是否有效的標記
        while(flag == 0) { double pick = ((double)rand())/RAND_MAX; int pos = (int)(pick*lenchrom); while(pos > lenchrom-1) { double pick = ((double)rand())/RAND_MAX; int pos = (int)(pick*lenchrom); // 處理越界
 } // 交叉開始
            double r = ((double)rand())/RAND_MAX; double v1 = chrom[choice1][pos]; double v2 = chrom[choice2][pos]; chrom[choice1][pos] = r*v2 + (1-r)*v1; chrom[choice2][pos] = r*v1 + (1-r)*v2; // 交叉結束
            if(chrom[choice1][pos] >=bound_down && chrom[choice1][pos]<=bound_up && chrom[choice2][pos] >=bound_down && chrom[choice2][pos]<=bound_up) flag = 1; // 交叉有效，退出交叉，不然繼續下一次交叉
 } } } // 變異操做
void Mutation(double chrom[sizepop][lenchrom]) { for(int i=0;i<sizepop;i++) { double pick = ((double)rand())/RAND_MAX; // 隨機選擇一個染色體進行變異
        int choice = (int)(pick*sizepop); while(choice > sizepop-1) { pick = ((double)rand())/RAND_MAX; int choice = (int)(pick*sizepop);  // 處理下標越界
 } pick = ((double)rand())/RAND_MAX; // 用於決定該輪是否進行變異
        if(pick>pmutation) continue; pick = ((double)rand())/RAND_MAX; int pos = (int)(pick*lenchrom); while(pos > lenchrom-1) { pick = ((double)rand())/RAND_MAX; pos = (int)(pick*lenchrom); } double v = chrom[i][pos]; double v1 = v - bound_up; double v2 = bound_down - v; double r = ((double)rand())/RAND_MAX; double r1 = ((double)rand())/RAND_MAX; if(r >= 0.5) chrom[i][pos] = v - v1*r1*(1-((double)i)/maxgen)*(1-((double)i)/maxgen); else chrom[i][pos] = v + v2*r1*(1-((double)i)/maxgen)*(1-((double)i)/maxgen); // 注意這裏寫法是不會越界的，因此只用一次變異就能夠了
 } } // 主函數
int main(void) { time_t start,finish; start = clock(); // 程序開始計時
    srand((unsigned)time(NULL)); // 初始化隨機數種子
    init_chrom(); // 初始化種羣
    double * best_fit_index = min(fitness); int best_index =(int)(*best_fit_index); gbest = *(best_fit_index+1); // 最優值
    gbest_index = 0; average_best[0] = sum(fitness)/sizepop; //初始平均最優值
    for(int i=0;i<lenchrom;i++) gbest_pos[i] = chrom[best_index][i]; // 進化開始
    for(int i=0;i<maxgen;i++) { Select(chrom); // 選擇
        Cross(chrom); //交叉
        Mutation(chrom); //變異
        for(int j=0;j<sizepop;j++) { fitness[j] = fit_func(chrom[j]); } double sum_fit = sum(fitness); average_best[i+1] = sum_fit/sizepop; // 平均最優值
        double * arr = min(fitness); double new_best = *(arr+1); int new_best_index = (int)(*arr); // 新最優值序號
        if(new_best < gbest) { gbest = new_best; for(int j=0;j<lenchrom;j++) { gbest_pos[j] = chrom[new_best_index][j]; } gbest_index = i+1; } } finish = clock(); // 程序計算結束
    double duration = ((double)(finish-start))/CLOCKS_PER_SEC; printf("程序計算耗時:%lf秒\n.",duration); printf("遺傳算法進化了%d次，最優值爲:%lf,最優值在第%d代取得,此代的平均最優值爲%lf.\n",maxgen,gbest,gbest_index,average_best[gbest_index]); printf("取得最優值的地方爲(%lf,%lf,%lf,%lf,%lf).\n",gbest_pos[0],gbest_pos[1],gbest_pos[2],gbest_pos[3],gbest_pos[4]); return 0; }

說明：

（1）這裏是求函數的最小值，函數值越小的個體適應度越大，因此這裏採用倒數變換，即把適應度函數F(x)定義爲1/f(x)(f(x)爲目標函數值),這樣就保證了適應度函數F(x)越大，個體是越優的。(固然咱們已經保證了f（x）始終是一個正數)

(2) 選擇操做採用的是經典的輪盤賭法，即每一個個體被選爲父代的機率與其適應度是成正比的，適應度越高，被選中的機率越大，設pi爲第i個個體被選中的機率，則

pi=Fi/(∑Fi),Fi爲第i個個體的適應度函數。

(3)交叉操做：因爲個體採用實數編碼，因此交叉操做也採用實數交叉法，第k個染色體a_k和第l個染色體a_l在第j位的交叉操做方法爲:

a_kj=a_ij(1-b)+a_ijb;  a_lj=a_lj(1-b)+a_kjb,其中b是[0,1]之間的隨機數。k，l，j都是隨機產生的，即隨機選取兩個父代，再隨機選取一個基因，按照公式完成交叉操做。

(4)變異操做：第i個個體的第j個基因進行變異的操做方法是：

r>=0.5時，a_ij=a_ij+(a_max-a_ij)*f(g)

r<0.5時，a_ij=a_ij+(a_min-a_ij)*f(g)

其中a_max是基因a_ij的上界，a_min是a_ij的下界，f(g)=r2*(1-g/G_max）²，r2是一個[0,1]之間的隨機數，g是當前迭代次數，G_max是進化最大次數，r是[0,1]之間的隨機數。

能夠看出，當r>=0.5時，a_ij是在a_ij到a_max之間變化的；r<0.5時，a_ij是在a_min到a_ij之間變化的。這樣的函數其實能夠構造不止一種，來完成變異的操做。只要符合隨機性，以及保證基因在有效的變化範圍內便可。

我在ubuntu16.04下，使用gcc編譯器運行獲得的結果以下：

從結果上來看，種羣數目爲50，進化代數爲500時，獲得的最優解爲2.014102已經很是接近實際的最優解2了，固然，若是把種羣數目再變大好比500，進化代數變多好比1000代，那麼獲得的結果將會更精確，種羣數目爲500，進化代數爲1000時的最優解以下：

能夠看出，增長種羣數目以及增長進化代數以後，最優值從2.014102變爲2.000189,比原來的精度提升了兩位，可是相應地，程序的運行時間也顯著增長了(從0.08秒左右增長到2秒左右)。因此在具體求解複雜問題的時候(好比有的時候會碰到函數參數有幾十個之多),咱們就必須綜合考慮最優解的精度以及運算複雜度(運算時間)這兩個方面，權衡以後，取合適的參數進行問題的求解。

　　固然，這裏只是以求解多元非線性函數的極值問題爲例來講明遺傳算法的具體實現步驟，遺傳算法的實際運用很是普遍，不只僅能夠求解複雜函數極值，並且在運籌學等衆多領域有着很是普遍的應用，也常常會和其餘的智能算法或者優化算法結合來求解問題。這個後面我還會再較爲詳細地論述。