python遺傳演算法(詳解)

學習程式碼來源於：遺傳演算法python

一.主要思想

遺傳演算法是根據達爾文的“適者生存，優勝劣汰”的思想來找到最優解的額，其特點是所找到的解是全域性最優解，相對於蟻群演算法可能出現的區域性最優解還是有優勢的。

二.主要名詞

個體（染色體）：一個染色體代表一個具體問題的一個解，一個染色體包含若干基因。

基因：一個基因代表具體問題解的一個決策變數。種群：多個個體（染色體）構成一個種群。即一個問題的多組解構成了解的種群。

我們的目的就是讓種群中”優勝劣汰“，最終只剩下一個最優解。接下來介紹最基本遺傳演算法，只用了選擇，交叉，變異三種遺傳運算元。

三.主要步驟

1）種群初始化。我們需要首先通過隨機生成的方式來創造一個種群，一般該種群的數量為100~500，這裡我們採用二進位制將一個染色體(解)編碼為基因型。隨後用進位制轉化，將二進位制的基因型轉化成十進位制的表現型。

2）適應度計算(種群評估)。這裡我們直接將目標函式值作為個體的適應度。

3）選擇(複製)操作。根據種群中個體的適應度大小，通過輪盤賭等方式將適應度高的個體從當前種群中選擇出來。其中輪盤賭即是與適應度成正比的概率來確定各個個體遺傳到下一代群體中的數量。

      具體步驟如下：

     (1)首先計算出所有個體的適應度總和Σfi。

     (2)其次計算出每個個體的相對適應度大小fi/Σfi，類似於softmax。

     (3)再產生一個0到1之間的隨機數，依據隨機數出現在上述哪個概率區域內來確定各個個體被選中的次數。

4）交叉(交配)運算。該步驟是遺傳演算法中產生新的個體的主要操作過程，它用一定的交配概率閾值(pc，一般是0.4到0.99)來控制是否採取單點交叉，多點交叉等方式生成新的交叉個體。

     具體步驟如下：

     (1)先對群體隨機配對。

     (2)再隨機設定交叉點的位置。

     (3)再互換配對染色體間的部分基因。 

5）變異運算。該步驟是產生新的個體的另一種操作。一般先隨機產生變異點，再根據變異概率閾值(pm,一般是0.0001到0.1)將變異點的原有基因取反。

6）終止判斷。如果滿足條件(迭代次數，一般是200~500)則終止演算法，否則返回step2。

四.程式碼實現

1.種群初始化     

(1)隨機生成若干個數的二進位制染色體。

# -*-coding:utf-8 -*-
#目標求解2*sin(x)+cos(x)最大值
import random
import math
import matplotlib.pyplot as plt
#初始化生成chromosome_length大小的population_size個個體的二進位制基因型種群
def species_origin(population_size,chromosome_length):
    population=[[]]
    #二維列表，包含染色體和基因
    for i in range(population_size):
        temporary=[]
        #染色體暫存器
        for j in range(chromosome_length):
            temporary.append(random.randint(0,1))
            #隨機產生一個染色體,由二進位制陣列成
        population.append(temporary)
            #將染色體新增到種群中
    return population[1:]
            # 將種群返回，種群是個二維陣列，個體和染色體兩維

(2)編碼 

將二進位制的染色體基因型編碼成十進位制的表現型。

#從二進位制到十進位制
 #input:種群,染色體長度
def translation(population,chromosome_length):
    temporary=[]
    for i in range(len(population)):
        total=0
        for j in range(chromosome_length):
            total+=population[i][j]*(math.pow(2,j))
            #從第一個基因開始，每位對2求冪，再求和
            # 如：0101 轉成十進位制為：1 * 20 + 0 * 21 + 1 * 22 + 0 * 23 = 1 + 0 + 4 + 0 = 5
        temporary.append(total)
        #一個染色體編碼完成，由一個二進位制數編碼為一個十進位制數
    return temporary
   # 返回種群中所有個體編碼完成後的十進位制數

2.適應度計算

個體適應度與其對應的個體表現型x的目標函式值相關聯，x越接近於目標函式的最優點，其適應度越大，從而其存活的概率越大。反之適應度越小，存活概率越小。這就引出一個問題關於適應度函式的選擇，本例中，函式值總取非負值(刪去了。。。)，以函式最大值為優化目標，故直接將目標函式作為適應度函式。這裡我們直接將目標函式2*sin(x)+cos作為個體適應度。如果，你想優化的是多元函式的話，需要將個體中基因型的每個變數提取出來，分別帶入目標函式。比如說：我們想求x1+lnx2的最大值。基因編碼為4位編碼，其中前兩位是x1，後兩位是x2。那麼我們在求適應度的時候，需要將這兩個值分別帶入f(x1)=x,f(x2)=lnx。再對f(x1)和f(x2)求和得到此個體的適應度，具體的在此篇中有詳解。在本篇中，雖然染色體長度為10,但是實際上只有一個變數。

# 目標函式相當於環境 對染色體進行篩選，這裡是2*sin(x)+cos(x)
def function(population,chromosome_length,max_value):
    temporary=[]
    function1=[]
    temporary=translation(population,chromosome_length)
    # 暫存種群中的所有的染色體(十進位制)
    for i in range(len(temporary)):
        x=temporary[i]*max_value/(math.pow(2,chromosome_length)-1)
        #一個基因代表一個決策變數，其演算法是先轉化成十進位制，然後再除以2的基因個數次方減1(固定值)。
        function1.append(2*math.sin(x)+math.cos(x))
        #這裡將2*sin(x)+cos(x)作為目標函式，也是適應度函式
    return function1

關於這裡的x的生成方式，據資料說只是一個編碼規則。於是，我對其進行了幾種實驗。

如下，這種是刪了max_value的輸出結果：

這種少了除號後面的結果：

這種是隻保留temporay[i]的結果：

從以上結果我們可以看出，這幾種結果對最後的優化結果都沒有什麼影響。但是對其優化速度和平滑性都有影響。這也正是適應度函式選擇的意義，選的好的話就可以加快優化效果。

補充：2018.6.16

由下面的選擇操作可知，適應度是選擇操作的主要參考依據，適應度函式(Fitness Function)的選取直接影響到遺傳演算法的收斂速度以及能否找到最優解。因而適應度函式的選擇問題在遺傳演算法中是一項很值得研究的課題。一般情況下，關於適應度與目標函式的選擇有以下這兩種關係:

也就是說，我們要在每一輪迭代(進化)時，要將所有個體的適應度函式值都要遍歷，然後得到最/大小值。此時每一輪的適應度函式都是在變化的。之前，我們程式中的個體適應度函式是不變化的。實際上要得到更準確的結果，除了以上基本變化，還有以下的適應度變化方法：

我們對個體的適應度調整的目的有兩個：

一是維持個體之間的合理差距，加速競爭。

二是避免個體之間的差距過大，限制競爭。

此外，還有適應度共享等調節方式，此處不再贅述。來自於，選擇和適應度函式。

另外，在這段程式中，採用一種巧妙的方式限定了表現型的大小。比如，我們這裡限制最大值為30,基因長度為6,則30/（2**6-1）=0.47,我們這裡取一個表現型的最大值-2**6=64，則64*0.47=30.47,大致(這個詞..有點不好意思說出來)符合我們的限制範圍。

3.選擇操作

(1).只保留非負值的適應度/函式值(不小於0)

def fitness(function1):
    fitness1=[]
    min_fitness=mf=0
    for i in range(len(function1)):
        if(function1[i]+mf>0):
            temporary=mf+function1[i]
        else:
            temporary=0.0
        # 如果適應度小於0,則定為0
        fitness1.append(temporary)
        #將適應度新增到列表中
    return fitness1

(2).首先計算出所有個體的適應度總和Σfi

#計算適應度和
def sum(fitness1):
    total=0
    for i in range(len(fitness1)):
        total+=fitness1[i]
    return total

#計算適應度斐波納挈列表，這裡是為了求出累積的適應度
def cumsum(fitness1):
    for i in range(len(fitness1)-2,-1,-1):
        # range(start,stop,[step])
        # 倒計數
        total=0
        j=0
        while(j<=i):
            total+=fitness1[j]
            j+=1
        #這裡是為了將適應度劃分成區間
        fitness1[i]=total
        fitness1[len(fitness1)-1]=1

(3).再產生一個0到1之間的隨機數，依據隨機數出現在上述哪個概率區域內來確定各個個體被選中的次數。

#3.選擇種群中個體適應度最大的個體
def selection(population,fitness1):
    new_fitness=[]
    #單個公式暫存器
    total_fitness=sum(fitness1)
    #將所有的適應度求和
    for i in range(len(fitness1)):
        new_fitness.append(fitness1[i]/total_fitness)
    #將所有個體的適應度概率化,類似於softmax
    cumsum(new_fitness)
    #將所有個體的適應度劃分成區間
    ms=[]
    #存活的種群
    population_length=pop_len=len(population)
    #求出種群長度
    #根據隨機數確定哪幾個能存活

    for i in range(pop_len):
        ms.append(random.random())
    # 產生種群個數的隨機值
    ms.sort()
    # 存活的種群排序
    fitin=0
    newin=0
    new_population=new_pop=population

    #輪盤賭方式
    while newin<pop_len:
        if(ms[newin]<new_fitness[fitin]):
            new_pop[newin]=pop[fitin]
            newin+=1
        else:
            fitin+=1
    population=new_pop

這裡需要詳細解釋一下輪盤賭演算法,其實就是典型的幾何概率問題.這裡根據cumsum函式將個體的適應度[0.1,0.2,0.2,0.4.....]都劃分成區間[0.1,0.3,0.5,0.9...],如下圖:

然後,random產生的ms[]排序後為[0.08,0.12,0.15,0.4.....],就有如下的程式:

我們可以看出新的種群是如下分佈的:

這意味著某個體適應度越大,新的種群中該個體數量越多。

4.交叉

def crossover(population,pc):
#pc是概率閾值，選擇單點交叉還是多點交叉，生成新的交叉個體，這裡沒用
    pop_len=len(population)

    for i in range(pop_len-1):
        cpoint=random.randint(0,len(population[0]))
        #在種群個數內隨機生成單點交叉點
        temporary1=[]
        temporary2=[]

        temporary1.extend(pop[i][0:cpoint])
        temporary1.extend(pop[i+1][cpoint:len(population[i])])
        #將tmporary1作為暫存器，暫時存放第i個染色體中的前0到cpoint個基因，
        #然後再把第i+1個染色體中的後cpoint到第i個染色體中的基因個數，補充到temporary2後面

        temporary2.extend(pop[i+1][0:cpoint])
        temporary2.extend(pop[i][cpoint:len(pop[i])])
        # 將tmporary2作為暫存器，暫時存放第i+1個染色體中的前0到cpoint個基因，
        # 然後再把第i個染色體中的後cpoint到第i個染色體中的基因個數，補充到temporary2後面
        pop[i]=temporary1
        pop[i+1]=temporary2
        # 第i個染色體和第i+1個染色體基因重組/交叉完成

5.變異

#step4：突變
def mutation(population,pm):
    # pm是概率閾值
    px=len(population)
    # 求出種群中所有種群/個體的個數
    py=len(population[0])
    # 染色體/個體中基因的個數
    for i in range(px):
        if(random.random()<pm):
        #如果小於閾值就變異
            mpoint=random.randint(0,py-1)
            # 生成0到py-1的隨機數
            if(population[i][mpoint]==1):
            #將mpoint個基因進行單點隨機變異，變為0或者1
                population[i][mpoint]=0
            else:
                population[i][mpoint]=1

6.其他

# 將每一個染色體都轉化成十進位制 max_value為基因最大值，為了後面畫圖用
def b2d(b,max_value,chromosome_length):
    total=0
    for i in range(len(b)):
        total=total+b[i]*math.pow(2,i)
    #從第一位開始，每一位對2求冪，然後求和，得到十進位制數？
    total=total*max_value/(math.pow(2,chromosome_length)-1)
    return total

#尋找最好的適應度和個體
def best(population,fitness1):

    px=len(population)
    bestindividual=[]
    bestfitness=fitness1[0]

    for i in range(1,px):
   # 迴圈找出最大的適應度，適應度最大的也就是最好的個體
        if(fitness1[i]>bestfitness):

            bestfitness=fitness1[i]
            bestindividual=population[i]

    return [bestindividual,bestfitness]

7.主程式

population_size=500
max_value=10
# 基因中允許出現的最大值
chromosome_length=10
pc=0.6
pm=0.01

results=[[]]
fitness1=[]
fitmean=[]

population=pop=species_origin(population_size,chromosome_length)
#生成一個初始的種群

for i in range(population_size):#注意這裡是迭代500次
    function1=function(population,chromosome_length,max_value)
    fitness1=fitness(function1)
    best_individual,best_fitness=best(population,fitness1)
    results.append([best_fitness,b2d(best_individual,max_value,chromosome_length)])
     #將最好的個體和最好的適應度儲存，並將最好的個體轉成十進位制
    selection(population,fitness1)#選擇
    crossover(population,pc)#交配
    mutation(population,pm)#變異

results=results[1:]
results.sort()
X=[]
Y=[]
for i in range(500):#500輪的結果
    X.append(i)
    Y.append(results[i][0])
plt.plot(X,Y)
plt.show()

8.迭代結果

我們可以看到最後經過500輪的迭代後，其結果已經非常接近2.2了。

9.完整程式碼

# -*-coding:utf-8 -*-
#目標求解2*sin(x)+cos(x)最大值
import random
import math
import matplotlib.pyplot as plt
class GA(object):
#初始化種群 生成chromosome_length大小的population_size個個體的種群

    def __init__(self,population_size,chromosome_length,max_value,pc,pm):

        self.population_size=population_size
        self.choromosome_length=chromosome_length
        # self.population=[[]]
        self.max_value=max_value
        self.pc=pc
        self.pm=pm
        # self.fitness_value=[]



    def species_origin(self):
        population=[[]]
        for i in range(self.population_size):

            temporary=[]
        #染色體暫存器
            for j in range(self.choromosome_length):

                temporary.append(random.randint(0,1))
            #隨機產生一個染色體,由二進位制陣列成

            population.append(temporary)
            #將染色體新增到種群中
        return population[1:]
            # 將種群返回，種群是個二維陣列，個體和染色體兩維

    #從二進位制到十進位制
    #編碼  input:種群,染色體長度 編碼過程就是將多元函式轉化成一元函式的過程
    def translation(self,population):

        temporary=[]
        for i in range(len(population)):
            total=0
            for j in range(self.choromosome_length):
                total+=population[i][j]*(math.pow(2,j))
            #從第一個基因開始，每位對2求冪，再求和
            # 如：0101 轉成十進位制為：1 * 20 + 0 * 21 + 1 * 22 + 0 * 23 = 1 + 0 + 4 + 0 = 5
            temporary.append(total)
        #一個染色體編碼完成，由一個二進位制數編碼為一個十進位制數
        return temporary
   # 返回種群中所有個體編碼完成後的十進位制數



#from protein to function,according to its functoin value

#a protein realize its function according its structure
# 目標函式相當於環境 對染色體進行篩選，這裡是2*sin(x)+math.cos(x)
    def function(self,population):
        temporary=[]
        function1=[]
        temporary=self.translation(population)
        for i in range(len(temporary)):
            x=temporary[i]*self.max_value/(math.pow(2,self.choromosome_length)-10)
            function1.append(2*math.sin(x)+math.cos(x))

         #這裡將sin(x)作為目標函式
        return function1

#定義適應度
    def fitness(self,function1):

        fitness_value=[]

        num=len(function1)

        for i in range(num):

            if(function1[i]>0):
                temporary=function1[i]
            else:
                temporary=0.0
        # 如果適應度小於0,則定為0

            fitness_value.append(temporary)
        #將適應度新增到列表中

        return fitness_value

#計算適應度和

    def sum(self,fitness_value):
        total=0

        for i in range(len(fitness_value)):
            total+=fitness_value[i]
        return total

#計算適應度斐伯納且列表
    def cumsum(self,fitness1):
        for i in range(len(fitness1)-2,-1,-1):
        # range(start,stop,[step])
        # 倒計數
            total=0
            j=0

            while(j<=i):
                 total+=fitness1[j]
                 j+=1

            fitness1[i]=total
            fitness1[len(fitness1)-1]=1


#3.選擇種群中個體適應度最大的個體
    def selection(self,population,fitness_value):
        new_fitness=[]
    #單個公式暫存器
        total_fitness=self.sum(fitness_value)
    #將所有的適應度求和
        for i in range(len(fitness_value)):
            new_fitness.append(fitness_value[i]/total_fitness)
    #將所有個體的適應度正則化
        self.cumsum(new_fitness)
    #
        ms=[]
    #存活的種群
        population_length=pop_len=len(population)
    #求出種群長度
    #根據隨機數確定哪幾個能存活

        for i in range(pop_len):
            ms.append(random.random())
    # 產生種群個數的隨機值
    # ms.sort()
    # 存活的種群排序
        fitin=0
        newin=0
        new_population=new_pop=population

    #輪盤賭方式
        while newin<pop_len:
              if(ms[newin]<new_fitness[fitin]):
                 new_pop[newin]=population[fitin]
                 newin+=1
              else:
                 fitin+=1
        population=new_pop

#4.交叉操作
    def crossover(self,population):
#pc是概率閾值，選擇單點交叉還是多點交叉，生成新的交叉個體，這裡沒用
        pop_len=len(population)

        for i in range(pop_len-1):

            if(random.random()<self.pc):

               cpoint=random.randint(0,len(population[0]))
           #在種群個數內隨機生成單點交叉點
               temporary1=[]
               temporary2=[]

               temporary1.extend(population[i][0:cpoint])
               temporary1.extend(population[i+1][cpoint:len(population[i])])
           #將tmporary1作為暫存器，暫時存放第i個染色體中的前0到cpoint個基因，
           #然後再把第i+1個染色體中的後cpoint到第i個染色體中的基因個數，補充到temporary2後面

               temporary2.extend(population[i+1][0:cpoint])
               temporary2.extend(population[i][cpoint:len(population[i])])
        # 將tmporary2作為暫存器，暫時存放第i+1個染色體中的前0到cpoint個基因，
        # 然後再把第i個染色體中的後cpoint到第i個染色體中的基因個數，補充到temporary2後面
               population[i]=temporary1
               population[i+1]=temporary2
        # 第i個染色體和第i+1個染色體基因重組/交叉完成
    def mutation(self,population):
     # pm是概率閾值
         px=len(population)
    # 求出種群中所有種群/個體的個數
         py=len(population[0])
    # 染色體/個體基因的個數
         for i in range(px):
             if(random.random()<self.pm):
                mpoint=random.randint(0,py-1)
            #
                if(population[i][mpoint]==1):
               #將mpoint個基因進行單點隨機變異，變為0或者1
                   population[i][mpoint]=0
                else:
                   population[i][mpoint]=1

#transform the binary to decimalism
# 將每一個染色體都轉化成十進位制 max_value,再篩去過大的值
    def b2d(self,best_individual):
        total=0
        b=len(best_individual)
        for i in range(b):
            total=total+best_individual[i]*math.pow(2,i)

        total=total*self.max_value/(math.pow(2,self.choromosome_length)-1)
        return total

#尋找最好的適應度和個體

    def best(self,population,fitness_value):

        px=len(population)
        bestindividual=[]
        bestfitness=fitness_value[0]
        # print(fitness_value)

        for i in range(1,px):
   # 迴圈找出最大的適應度，適應度最大的也就是最好的個體
            if(fitness_value[i]>bestfitness):

               bestfitness=fitness_value[i]
               bestindividual=population[i]

        return [bestindividual,bestfitness]


    def plot(self, results):
        X = []
        Y = []

        for i in range(500):
            X.append(i)
            Y.append(results[i][0])

        plt.plot(X, Y)
        plt.show()

    def main(self):

        results = [[]]
        fitness_value = []
        fitmean = []

        population = pop = self.species_origin()

        for i in range(500):
            function_value = self.function(population)
            # print('fit funtion_value:',function_value)
            fitness_value = self.fitness(function_value)
            # print('fitness_value:',fitness_value)

            best_individual, best_fitness = self.best(population,fitness_value)
            results.append([best_fitness, self.b2d(best_individual)])
        # 將最好的個體和最好的適應度儲存，並將最好的個體轉成十進位制,適應度函式
            self.selection(population,fitness_value)
            self.crossover(population)
            self.mutation(population)
        results = results[1:]
        results.sort()
        self.plot(results)

if __name__ == '__main__':


   population_size=400
   max_value=10
   chromosome_length=20
   pc=0.6
   pm=0.01
   ga=GA(population_size,chromosome_length,max_value,pc,pm)
   ga.main()

10.參考文獻：

1.西安電子科學技術大學, 尚榮華. 選擇和適應度函式[EB/OL].

https://wenku.baidu.com/view/d4fd8ab20129bd64783e0912a216147917117e11.html.

2.百度百科. 適應度函式[EB/OL]. https://baike.baidu.com/item/適應度函式/20593164?fr=aladdin.

希望有志同道合的小夥伴關注我的公眾平臺，歡迎您的批評指正，共同交流進步。