以前在學校的時候，寫了一個博客介紹遺傳算法，並通過Matlab實現了該算法。但是很多讀者反饋說代碼運行不起來。
我因為現在沒有Matlab工具了，而且Matlab還是挺貴的，所以還是用Python實現以下遺傳算法，供大家學習。
寫這個博客主要是因為我在學校學習遺傳算法的時候，沒學會，只是大概知道，但是完全用不起來。但是這個算法還是很有用的，所以還是想着記錄一下，讓在課堂上沒學會的同學，可以在這里學會，並且可以着手。
1、遺傳算法介紹
遺傳算法，模擬達爾文進化論的自然選擇和遺傳學機理的生物進化構成的計算模型，一種不斷選擇優良個體的算法。談到遺傳，想想自然界動物遺傳是怎么來的，自然主要過程包括染色體的選擇，交叉，變異（不明白這個的可以去看看生物學），這些操作后，保證了以后的個體基本上是最優的，那么以后再繼續這樣下去，就可以一直最優了。
2、解決的問題
先說說自己要解決的問題吧，遺傳算法很有名，自然能解決的問題很多了，在原理上不變的情況下，只要改變模型的應用環境和形式，基本上都可以。但是遺傳算法主要還是解決優化類問題，尤其是那種不能直接解出來的很復雜的問題，而實際情況通常也是這樣的。
本部分主要為了了解遺傳算法的應用，選擇一個復雜的二維函數來進行遺傳算法優化，函數顯示為y=10*sin(5*x)+7*abs(x-5)+10,這個函數圖像為：
Python實現代碼：
import numpy as np import matplotlib.pyplot as plt x=np.arange(0,10,0.001) y=10*np.sin(5*x) + 7*np.abs(x-5) + 10 plt.plot(x,y) plt.show()
怎么樣，還是有一點復雜的吧，當然你還可以任意假設和編寫，只要符合就可以。那么現在問你，要你一下求出最大值你能求出來嗎？（這個貌似可以，很容易寫出來，求個導數然后，----如果再復雜一點估計就不行了）這類問題如果用遺傳算法或者其他優化方法就很簡單了，為什么呢？說白了，其實就是計算機太笨了，同時計算速度又超快，舉個例子吧，我把x等分成100萬份，再一下子都帶值進去算，求出對應的100萬個y的值，再比較他們的大小，找到最大值不就可以了嗎，很笨吧，人算是不可能的，但是計算機可以。而遺傳算法也是很笨的一個個搜索，只不過加了一點什么東西，就是人為的給它算的方向和策略，讓它有目的的算，這也就是算法了。
3、如何開始
不明白遺傳算法的會問怎么開始呢？恩，其實每個算法都有自己的開始方式，遺傳算法也是，首先是選擇個體了。我們知道一個種群中可能只有一個個體嗎？不可能吧，肯定很多才對，這樣相互結合的機會才多，產生的后代才會多種多樣，才會有更好的優良基因，有利於種群的發展。那么算法也是如此，當然個體多少是個問題，一般來說20-100之間我覺得差不多了。那么個體究竟是什么呢？在我們這個問題中自然就是x值了。其他情況下，個體就是所求問題的變量，這里我們假設個體數選100個，也就是開始選100個不同的x值，不明白的話就假設是100個猴子吧。好了，現在有了100個猴子組成的一個種群，那么這個種群應該怎么發展才能越來越好？說到這，我們想想，如何定義這個越來越好呢？這個應該有一個評價指標吧。在我們的這個問題中，好像是對應的y值越大越好是吧。我們甚至可以給他們排個名來決定哪些好哪些不好。我們把這個叫做對於個體的適應度，這應該算是算法的后半部分才對。
4、編碼
好，什么是編碼？其實編碼就是把自變量（x）換一下形式而已，在這個形式下，更容易操作其他過程（比如交叉，變異什么的）而已。舉個例子吧，加入我們取x=1, 2, 3，我們可以把x編碼成x=a, b, c,我就說123對應就是abc，為什么要這樣做呢？比如問題里面你能夠獲取的都是abc的組合之類的，那么這樣編碼以后，你就可以再返回去成123來操作了。一般的編碼都是些什么呢？二進制編碼，自然數編碼，矩陣編碼。。等等，很多，不詳細寫了。而用的最多的可以說是二進制編碼，感覺這和人體DNA，基因的排列很相似。想想DNA怎么排的？不就是在兩條長鏈上一對一排的嗎？那么什么是二進制編碼？很簡單，就是1,0,1,0對應的來組合排列而已。比如：1100100010, 0011001001等等，這些都是位數長度為10的二進制編碼。再想想1在計算機的二進制形式是什么？如果八位來表示的話，是不是就是0000 0001；8是不是就是0000 1000；以此類推，那么我們這里也是這樣，把對應的x值換算成這種編碼形式，我們這里可以看到x的范圍是0-10吧，如果按照計算機這樣的方式是不是到0000 1010這里就完事了？想想這樣多短，前面四位都沒有用上多浪費呀，那么要想都用上怎么辦呢？也很簡單，我們把0000 0001不認為是1不就可以了嗎？因為1111 1111是255，那么如果說每一份為1/255的話，那么0000 0001不就是1/255了嗎？這個時候1怎樣表示了？不就是1111 1111嗎？好了我們把范圍擴大一些吧，每一份不是1/255,而是1/255*10,那么這個時候最大值是多少？是不是10，恩，這樣x編碼的范圍不就在0-10之間了嗎？這里就又有問題了，想想這樣做的話，x的最小精度是多少？就是1/255*10.雖然很小，但是在0-1/255*10之間的數你能不能取到？無論如何都取不到吧。那么就又來一個問題，怎樣去擴大這個精度呢？如果要保持0-10不變的話，只能增加位數了，把9位編碼，10位，20位，100位，哇，夠大了吧，變成100個0,1組合，很恐怖吧，事實上，究竟是多少要視情況而定，一般20位左右感覺就可以了，雖然說越大越好，但是太大了消耗內存，速度慢了，不值得。本題中，我們設置它為一個變量，先暫時取為10來實驗。好了，如果還不明白為什么要編碼看下面的吧。知道了交叉與變異，你就知道了。
5、關於交叉與變異
先說變異，什么是變異？簡單，基因發生突變就叫變異，有了編碼的概念，那就在編碼基礎上來說變異。首先就講最簡單的變異，就是個體的變異。現在以10位長的編碼來說，比如把x=3編碼一下，隨便假設為11000 10010吧，好了，在變異操作時，假設第5位變異了（說一下變異就是一位或者多位1或0變成0或1，也只能0,1之間變，沒辦法啊），那么這個時候變成什么了？是不是11001 10010再反編碼回去成x是多少呢？那肯定不是3了，變了呀，是多少是肯定可以反算回去的，這里懶得算了，就假設為3.213吧，發沒發現，這樣一來，x是不是變了？既然變了就好啊，帶到原函數（適應度函數）里面比較這兩個x值對應的哪個y值大一些，如果后面變異后的大一些是不是就是說產生了好的變異啊，就可以在下一次個體選擇的時候選擇它了。那么想想很多x來一起變異會怎么樣呢？肯定會生成很多很多的解吧，反復這么做會怎么樣呢？只要每次都保留最優解的話，我來個循環100萬次，也總能找到解吧，當然這么多次得花多久，也不太合適，這還只是一個點位在進行變異，如果每次我讓多個點位變異呢？哇，又不可思議了，變化更大了吧。當然，變異不止如此，更多的去看專業的論文吧。知道了變異是干什么的，剩下的都好說了，好了，這還是變異，想想自然界遺傳中除了變異還有什么？交叉吧，那么交叉又是什么？
學過生物的都知道，動物交配時，部分染色體干什么了？是不是交叉了？就是把相應部分的基因交換了，你的給我，我的給你，很有愛吧。再以編碼為例，比如現在隨便從100個x值中選取兩個吧，假設正好選中了x=3和4，對應的編碼假設是11001 10101 和00101 01011，那么怎么交叉呢？我們知道每次交叉的染色體通常是一塊一塊的，恩，這里在算法設計上也來一塊一塊的吧。比如說就把位置在2,3,4號的編碼給整體交叉了，那么x=3對應的位置是100吧，x=4對應的位置是010吧，好，交換以后x=3對應的位置就變成了010，x=4對應的位置就變成100，加回去就變成什么了？x=3是不是就是10101 10101，x=4是不是就是01001 01011了。而現在，把他們再反編碼回去還是x=3和x=4嗎？顯然又不是了吧（當然也有小概率是一樣的吧，很小）。那是什么？不想算，還是假設吧，假設為3.234和4.358吧，好了新的個體是不是又來了？恩，同理，帶到適應度函數里面去吧，再取優秀個體，完事。同樣，有些專門研究這種算法的開發出來各種各樣的交叉方式，什么一個個體的前3個與后一個個體的后三個交叉，中間幾位來交叉等等，總之就是生產新個體。而這樣做的目的在哪呢？無非是三個字，隨機性，充分保證生產新個體具有隨機性，你說你的x=3變異后為3.2,3.2距離3那么近，在一些存在局部最優解的問題上就永遠跳不出局部最優解，相反，你的x=1一下子變異成了x=5，哇，好大的變化，一下從這頭到了那頭，這對於算法的廣闊搜索能力來說是非常好的。
交叉前：
交叉后：
講完了這部分，現在知道了為什么要編碼了吧？如果你不編碼，你說你想要你的x=3怎么去變異？怎么去交叉？當然也不是沒有方法，比如你生成一個小的隨機數加到x=3上，但是你想想這兩種方法哪一個更具有隨機性、普遍性？顯然的。而更多的時候交叉與變異是在一起操作的，先交叉，再變異是普遍遺傳算法的操作步驟。
6、關於選擇的問題
說完了上面的部分，再說說選擇吧，選擇是什么？就是優勝劣汰。好的留下來，差的走人，在自然界中直接gg了是吧。不停地選擇使的種群一直朝着較好的方向行走。
對應到本問題來說，遺傳算法的選擇是什么樣子？在前面說到，每次交叉或者變異是不是產生了新的個體？如果這些個體都保留下來的話，種群是要爆炸的，第一次循環可能有100個x，第二次循環就200個個體，再來那么10萬次，哇哦，多少了，好多。顯然不可能吧。而且在算法里面，我們還規定的是每次循環都必須保證都是100個個體。那么必須在200個個體中剔除100個吧。好了，問題來了，如何剔除呢？有人說很簡單，排名吧，取前100號x不就可以了嗎？排名這個東西真的好嗎？我就不信，憑什么差一點的不能選上，搞不好在下一次變異中一下子沖到了第一呢？這個問題在選擇上也有一些對應的規則，最通用的就是輪盤賭法，簡單來說就是一種概率選擇法（當然還有許多其他的方法，感興趣的自己搜相關的文獻吧，我也沒用過）。什么是輪盤賭法呢？就是把對應所有y值（適應度函數值）加起來，再用各自的y值去除以這個sum值，這樣是不是誰的概率大誰的概率小就很清楚了？然后再隨機生成一個0-1的概率值p，誰在p的范圍里面是不是就選擇誰，比如說x=3時在100個x中y的值最大，那么選擇它的概率是不是就最大，比如說是0.1（0.1小嗎？不小了好吧，想想其他的會是什么，都比0.1小，那么從概率上講，選100次的話，是不是就有10次選擇x=3，其他的都不足10次是吧，那么在下一次100個種群個體中就有10個x=3了，再來一回可能就有20個x=3的個體了。再就是30個，最后就只剩下100個x=3的個體了，它自己在哪里交叉變異已經沒有意義了，如果到了這個時候就意味着這個算法可以結束了）。再詳細點，如下圖所示吧：現在要在下面三個大類中選取100個x個體，輪盤賭轉100次以后，是不是個體數落在s3中的個體多一些，選擇的原理就是這樣，再不明白直接后面的程序吧。
7、還差點什么呢？
至此，感覺也差不多了吧，選擇完后再重復上述步驟交叉，變異等等，那么什么時候是個頭了？很簡單，辦法就是迭代次數，迭代10次看一下結果，20次，看一下結果，30次，40次，100次，當次數達到一定程度以后，優秀的個體越來越多，大都集中在最優解附近，即使變異或者交叉了也是在這個最優解附近，沒有影響的。在下一次選擇就又變回來了，那么至此就真的結束了。比如說先來結果吧，該問題按我的思路做完后，迭代100次變成什么樣子了？看圖：
代碼：
import numpy as np import random from scipy.optimize import fsolve import matplotlib.pyplot as plt import heapq import time # 求染色體長度 def getEncodeLength(decisionvariables, delta): # 將每個變量的編碼長度放入數組 # 這個要看有幾個自變量，假如自變量的個數是3，那么返回的lengths的元素就有3個 lengths = [] for decisionvar in decisionvariables: uper = decisionvar[1] low = decisionvar[0] # res返回一個數組 res = fsolve(lambda x:((uper-low)/delta-2**x+1), 30) # ceil()向上取整 length = int(np.ceil(res[0])) lengths.append(length) return lengths def getinitialPopulation(length, populationSize): # 根據DNA長度和種群大小初始化種群 chromsomes = np.zeros((populationSize, length), dtype=np.int) for popusize in range(populationSize): # np.random.randint產生[0,2)之間的隨機整數，第三個數表示隨機數的數量 chromsomes[popusize, :] = np.random.randint(0, 2, length) return chromsomes # 染色體解碼得到表現形的解 def getDecode(population,encodeLength,decisionVariables, delta): # 把染色體計算成數字 populationsize = population.shape[0] length = len(encodeLength) decodeVariables = np.zeros((populationsize, length), dtype=np.float) for i,populationchild in enumerate(population): start = 0 for j, lengthchild in enumerate(encodeLength): power = lengthchild - 1 decimal = 0 for k in range(start, start+lengthchild): decimal += populationchild[k] * (2**power) power = power - 1 start = lengthchild lower = decisionVariables[j][0] uper = decisionVariables[j][1] decodevalue = lower+decimal*(uper-lower)/(2**lengthchild-1) decodeVariables[i][j] = decodevalue return decodeVariables # 得到每個個體的適應度值及累計概率 def getFitnessValue(func, decode): popusize, decisionvar = decode.shape fitnessValue = np.zeros((popusize, 1)) for popunum in range(popusize): fitnessValue[popunum][0] = func(decode[popunum][0]) probability = fitnessValue/np.sum(fitnessValue) cum_probability = np.cumsum(probability) return fitnessValue, cum_probability def selectNewPopulation(decodepopu,cum_probability): m,n = decodepopu.shape newPopulation = np.zeros((m,n)) for i in range(m): randomnum = np.random.random() for j in range(m): if(randomnum < cum_probability[j]): newPopulation[i] = decodepopu[j] break return newPopulation def crossNewPopulation(newpopu, prob): m, n = newpopu.shape numbers = np.uint8(m*prob) if numbers %2 != 0: numbers = numbers + 1 updatepopulation = np.zeros((m,n),dtype=np.uint8) index = random.sample(range(m),numbers) for i in range(m): if not index.__contains__(i): updatepopulation[i] = newpopu[i] j = 0 while j < numbers: crosspoint = np.random.randint(0,n,1) crosspoint = crosspoint[0] updatepopulation[index[j]][0:crosspoint] = newpopu[index[j]][0:crosspoint] updatepopulation[index[j]][crosspoint:] = newpopu[index[j+1]][crosspoint:] updatepopulation[index[j+1]][0:crosspoint] = newpopu[index[j+1]][0:crosspoint] updatepopulation[index[j+1]][crosspoint:] = newpopu[index[j]][crosspoint:] j = j + 2 return updatepopulation def mutation(crosspopulation,mutaprob): mutationpopu = np.copy(crosspopulation) m,n = crosspopulation.shape mutationnums = np.uint8(m*n*mutaprob) mutationindex = random.sample(range(m*n),mutationnums) for geneindex in mutationindex: row = np.uint8(np.floor(geneindex/n)) colume = geneindex % n if mutationpopu[row][colume] == 0: mutationpopu[row][colume] =1 else: mutationpopu[row][colume] = 0 return mutationpopu def findMaxPopulation(population, maxevaluation, maxSize): maxevalue = maxevaluation.flatten() maxevaluelist = maxevalue.tolist() maxIndex = map(maxevaluelist.index, heapq.nlargest(100, maxevaluelist)) index = list(maxIndex) colume = population.shape[1] maxPopulation = np.zeros((maxSize, colume)) i = 0 for ind in index: maxPopulation[i] = population[ind] i = i + 1 return maxPopulation def fitnessFunction(): return lambda x: 10*np.sin(5*x) + 7*np.abs(x-5) + 10 def main(): optimalvalue = [] optimalvariables = [] decisionVariables = [[0.0, 10.0]] delta = 0.0001 EncodeLength = getEncodeLength(decisionVariables, delta) print(EncodeLength) initialPopuSize = 100 population = getinitialPopulation(sum(EncodeLength), initialPopuSize) print(population) # 最大進化代數 maxgeneration = 30 # 交叉概率 prob = 0.8 # 變異概率 mutationprob = 0.1 # 新生成的種群數量 maxPopuSize = 100 for generation in range(maxgeneration): decode = getDecode(population, EncodeLength, decisionVariables, delta) evaluation, cum_proba = getFitnessValue(fitnessFunction(), decode) newpopulations = selectNewPopulation(population, cum_proba) crossPopulation = crossNewPopulation(newpopulations, prob) mutationpopulation = mutation(crossPopulation, mutationprob) totalpopalation = np.vstack((population, mutationpopulation)) final_decode = getDecode(totalpopalation, EncodeLength, decisionVariables,delta) final_evaluation, final_cumprob = getFitnessValue(fitnessFunction(), final_decode) population = findMaxPopulation(totalpopalation, final_evaluation, maxPopuSize) optimalvalue.append(np.max(final_evaluation)) index = np.where(final_evaluation == max(final_evaluation)) optimalvariables.append(list(final_decode[index[0][0]])) decodeTemp = getDecode(population, EncodeLength, decisionVariables, delta) evaluationTemp, cum_probaTemp = getFitnessValue(fitnessFunction(), decodeTemp) plt.plot(decodeTemp,evaluationTemp,"*") x = np.arange(0.0, 10, 0.001) y = 10 * np.sin(5 * x) + 7 * np.abs(x - 5) + 10 plt.plot(x, y) plt.show() time.sleep(1) x = [i for i in range(maxgeneration)] y = [optimalvalue[i] for i in range(maxgeneration)] plt.plot(x,y) plt.show() optimalval = np.max(optimalvalue) index = np.where(optimalvalue == max(optimalvalue)) optimavar = optimalvariables[index[0][0]] return optimalval, optimavar if __name__ == '__main__': optval, optvar = main() print("x1:",optvar[0]) print("maxValue:",optval)
