今天整理之前寫的程式碼，發現在做數模期間寫的用Python實現的遺傳演算法，感覺還是挺有意思的，就拿出來分享一下。

首先遺傳演算法是一種優化演算法，通過模擬基因的優勝劣汰，進行計算（具體的演算法思路什麼的就不贅述了）。大致過程分為初始化編碼、個體評價、選擇，交叉，變異。

以目標式子 y = 10 * sin(5x) + 7 * cos(4x)為例，計算其最大值

首先是初始化，包括具體要計算的式子、種群數量、染色體長度、交配概率、變異概率等。並且要對基因序列進行初始化

1. pop_size = 500# 種群數量
2. max_value = 10# 基因中允許出現的最大值
3. chrom_length = 10
   # 染色體長度
4. pc = 0.6# 交配概率
5. pm = 0.01# 變異概率
6. results = [[]]      # 儲存每一代的最優解，N個二元組
7. fit_value = []      # 個體適應度
8. fit_mean = []       # 平均適應度
9. pop = geneEncoding(pop_size, chrom_length)  

其中genEncodeing是自定義的一個簡單隨機生成序列的函式，具體實現如下

1. def geneEncoding(pop_size, chrom_length):  
2.     pop = [[]]  
3.     for i in range(pop_size):  
4.         temp = []  
5.         for j in range(chrom_length):  
6.             temp.append(random.randint(0, 1))  
7.         pop.append(temp)  
8.     return pop[1:]  

編碼完成之後就是要進行個體評價，個體評價主要是計算各個編碼出來的list的值以及對應帶入目標式子的值。其實編碼出來的就是一堆2進位制list。這些2進位制list每個都代表了一個數。其值的計算方式為轉換為10進位制，然後除以2的序列長度次方減一，也就是全一list的十進位制減一。根據這個規則就能計算出所有list的值和帶入要計算式子中的值，程式碼如下

1. # 0.0 coding:utf-8 0.0
2. # 解碼並計算值
3. import math  
4. def decodechrom(pop, chrom_length):  
5.     temp = []  
6.     for i in range(len(pop)):  
7.         t = 0
8.         for j in range(chrom_length):  
9.             t += pop[i][j] * (math.pow(2, j))  
10.         temp.append(t)  
11.     return temp  
12. def calobjValue(pop, chrom_length, max_value):  
13.     temp1 = []  
14.     obj_value = []  
15.     temp1 = decodechrom(pop, chrom_length)  
16.     for i in range(len(temp1)):  
17.         x = temp1[i] * max_value / (math.pow(2, chrom_length) - 1)  
18.         obj_value.append(10 * math.sin(5 * x) + 7 * math.cos(4 * x))  
19.     return obj_value  

有了具體的值和對應的基因序列，然後進行一次淘汰，目的是淘汰掉一些不可能的壞值。這裡由於是計算最大值，於是就淘汰負值就好了

1. # 0.0 coding:utf-8 0.0
2. # 淘汰（去除負值）
3. def calfitValue(obj_value):  
4.     fit_value = []  
5.     c_min = 0
6.     for i in range(len(obj_value)):  
7.         if(obj_value[i] + c_min > 0):  
8.             temp = c_min + obj_value[i]  
9.         else:  
10.             temp = 0.0
11.         fit_value.append(temp)  
12.     return fit_value  

然後就是進行選擇，這是整個遺傳演算法最核心的部分。選擇實際上模擬生物遺傳進化的優勝劣汰，讓優秀的個體儘可能存活，讓差的個體儘可能的淘汰。個體的好壞是取決於個體適應度。個體適應度越高，越容易被留下，個體適應度越低越容易被淘汰。具體的程式碼如下

1. # 0.0 coding:utf-8 0.0
2. # 選擇
3. import random  
4. def sum(fit_value):  
5.     total = 0
6.     for i in range(len(fit_value)):  
7.         total += fit_value[i]  
8.     return total  
9. def cumsum(fit_value):  
10.     for i in range(len(fit_value)-2, -1, -1):  
11.         t = 0
12.         j = 0
13.         while(j <= i):  
14.             t += fit_value[j]  
15.             j += 1
16.         fit_value[i] = t  
17.         fit_value[len(fit_value)-1] = 1
18. def selection(pop, fit_value):  
19.     newfit_value = []  
20.     # 適應度總和
21.     total_fit = sum(fit_value)  
22.     for i in range(len(fit_value)):  
23.         newfit_value.append(fit_value[i] / total_fit)  
24.     # 計算累計概率
25.     cumsum(newfit_value)  
26.     ms = []  
27.     pop_len = len(pop)  
28.     for i in range(pop_len):  
29.         ms.append(random.random())  
30.     ms.sort()  
31.     fitin = 0
32.     newin = 0
33.     newpop = pop  
34.     # 轉輪盤選擇法
35.     while newin < pop_len:  
36.         if(ms[newin] < newfit_value[fitin]):  
37.             newpop[newin] = pop[fitin]  
38.             newin = newin + 1
39.         else:  
40.             fitin = fitin + 1
41.     pop = newpop  

以上程式碼主要進行了3個操作，首先是計算個體適應度總和，然後在計算各自的累積適應度。這兩步都好理解，主要是第三步，轉輪盤選擇法。這一步首先是生成基因總數個0-1的小數，然後分別和各個基因的累積個體適應度進行比較。如果累積個體適應度大於隨機數則進行保留，否則就淘汰。這一塊的核心思想在於：一個基因的個體適應度越高，他所佔據的累計適應度空隙就越大，也就是說他越容易被保留下來。

選擇完後就是進行交配和變異，這個兩個步驟很好理解。就是對基因序列進行改變，只不過改變的方式不一樣

交配：

1. # 0.0 coding:utf-8 0.0
2. # 交配
3. import random  
4. def crossover(pop, pc):  
5.     pop_len = len(pop)  
6.     for i in range(pop_len - 1):  
7.         if(random.random() < pc):  
8.             cpoint = random.randint(0,len(pop[0]))  
9.             temp1 = []  
10.             temp2 = []  
11.             temp1.extend(pop[i][0:cpoint])  
12.             temp1.extend(pop[i+1][cpoint:len(pop[i])])  
13.             temp2.extend(pop[i+1][0:cpoint])  
14.             temp2.extend(pop[i][cpoint:len(pop[i])])  
15.             pop[i] = temp1  
16.             pop[i+1] = temp2  

變異：

1. # 0.0 coding:utf-8 0.0
2. # 基因突變
3. import random  
4. def mutation(pop, pm):  
5.     px = len(pop)  
6.     py = len(pop[0])  
7.     for i in range(px):  
8.         if(random.random() < pm):  
9.             mpoint = random.randint(0, py-1)  
10.             if(pop[i][mpoint] == 1):  
11.                 pop[i][mpoint] = 0
12.             else:  
13.                 pop[i][mpoint] = 1

整個遺傳演算法的實現完成了，總的呼叫入口程式碼如下

1. # 0.0 coding:utf-8 0.0
2. import matplotlib.pyplot as plt  
3. import math  
4. from calobjValue import calobjValue  
5. from calfitValue import calfitValue  
6. from selection import selection  
7. from crossover import crossover  
8. from mutation import mutation  
9. from best import best  
10. from geneEncoding import geneEncoding  
11. print'y = 10 * math.sin(5 * x) + 7 * math.cos(4 * x)'
12. # 計算2進位制序列代表的數值
13. def b2d(b, max_value, chrom_length):  
14.     t = 0
15.     for j in range(len(b)):  
16.         t += b[j] * (math.pow(2, j))  
17.     t = t * max_value / (math.pow(2, chrom_length) - 1)  
18.     return t  
19. pop_size = 500# 種群數量
20. max_value = 10# 基因中允許出現的最大值
21. chrom_length = 10# 染色體長度
22. pc = 0.6# 交配概率
23. pm = 0.01# 變異概率
24. results = [[]]      # 儲存每一代的最優解，N個二元組
25. fit_value = []      # 個體適應度
26. fit_mean = []       # 平均適應度
27. # pop = [[0, 1, 0, 1, 0, 1, 0, 1, 0, 1] for i in range(pop_size)]
28. pop = geneEncoding(pop_size, chrom_length)  
29. for i in range(pop_size):  
30.     obj_value = calobjValue(pop, chrom_length, max_value)        # 個體評價
31.     fit_value = calfitValue(obj_value)      # 淘汰
32.     best_individual, best_fit = best(pop, fit_value)        # 第一個儲存最優的解, 第二個儲存最優基因
33.     results.append([best_fit, b2d(best_individual, max_value, chrom_length)])  
34.     selection(pop, fit_value)       # 新種群複製
35.     crossover(pop, pc)      # 交配
36.     mutation(pop, pm)       # 變異
37. results = results[1:]  
38. results.sort()  
39. X = []  
40. Y = []  
41. for i in range(500):  
42.     X.append(i)  
43.     t = results[i][0]  
44.     Y.append(t)  
45. plt.plot(X, Y)  
46. plt.show()  

最後呼叫了一下matplotlib包，把500代最優解的變化趨勢表現出來。

完整程式碼可以在github 檢視