【機器學習】遺傳算法(Genetic Algorithm)的Python實現(xiàn)

??本文章用Python實現(xiàn)了基本的優(yōu)化遺傳算法并用類進行了封裝

一、遺傳算法概述

??遺傳算法(Genetic Algorithm)是模擬達爾文生物進化論的自然選擇和遺傳學機理的生物進化過程的計算模型,是一種通過模擬自然進化過程搜索最優(yōu)解的方法。
??遺傳算法的降解可以看有史以來最容易理解的遺傳算法,用動畫展現(xiàn)的原理。
??遺傳算法的一些基本實現(xiàn)可以看莫煩進化算法,用python打好了大致架構并應用了有趣的題目

二、遺傳算法實現(xiàn)

1.首先引入需要的包

#jupyter matplotlib內置繪圖的魔術命令
%matplotlib inline
#matplotlib繪圖、3D相關包
from matplotlib import pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
#用于動態(tài)顯示圖像 
from IPython import display

#運算、操作處理相關包
import numpy as np
import pandas as pd

2.打好大致框架

class GA():

    def __init__(self, nums, bound, func, DNA_SIZE=None, cross_rate=0.8, mutation=0.003):
        pass
#將編碼后的DNA翻譯回來(解碼)
    def translateDNA(self):
        pass
#得到適應度
    def get_fitness(self, non_negative=False):
        pass
#自然選擇
    def select(self):
        pass
#染色體交叉
    def crossover(self):
        pass
#基因變異
    def mutate(self):
        pass
#進化
    def evolution(self):
        pass
#重置
    def reset(self):
        pass
#打印當前狀態(tài)日志
    def log(self):
        pass
#一維變量作圖
    def plot_in_jupyter_1d(self, iter_time=200):
        pass

2.init分析及實現(xiàn)

??作為構造方法,應該確定構造一個對象需要什么,并驗證得到的數(shù)據(jù)
??我們希望得到初始的數(shù)據(jù)nums,每個變量的邊界bound, 運算需要的函數(shù)func,DNA的大?。ɑ驂A基對的個數(shù))DNA_SIZE,染色體交叉的概率cross_rate,基因變異的概率mutation
??其中得到的nums為二維N * M列表,N進化種群的總數(shù),MDNA個數(shù)(變量的多少)。
??boundM * 2的二維列表,形如:[(0, 5), (-11.2, +134), ...]。
??func為方法(function)對象,既可以用def定義出的,也可以是lambda表達式定義的匿名函數(shù)。
??DNA_SIZE可以指定大小,為None時會自動指派。
??cross_ratemutation也有各自的默認值

編碼及解碼方式
將固定的二進制整數(shù)位數(shù)映射到范圍內

encoding\_num = (num-var\_min)\times{\frac{2^{DNA\_SIZE}}{var\_len}}\\ decoding\_num = \frac{encoding\_num}{\frac{2^{DNA\_SIZE}}{var\_len}} + var\_min其中,encoding\_numdecoding\_num分別代表編碼和解碼數(shù),num為需要編碼的數(shù),var\_min為當前的范圍下限,DNA\_SIZE為DNA長度,既堿基對個數(shù),var\_len為當前變量取值范圍的長度,既上限-下限。

# input:
#     nums: m * n  n is nums_of x, y, z, ...,and m is population's quantity
#     bound:n * 2  [(min, nax), (min, max), (min, max),...]
#     DNA_SIZE is binary bit size, None is auto
def __init__(self, nums, bound, func, DNA_SIZE=None, cross_rate=0.8, mutation=0.003):
    nums = np.array(nums)
    bound = np.array(bound)
    self.bound = bound
    if nums.shape[1] != bound.shape[0]:
        raise Exception(f'范圍的數(shù)量與變量的數(shù)量不一致, 您有{nums.shape[1]}個變量,卻有{bound.shape[0]}個范圍')

    for var in nums:
        for index, var_curr in enumerate(var):
            if var_curr < bound[index][0] or var_curr > bound[index][1]:
                raise Exception(f'{var_curr}不在取值范圍內')

    for min_bound, max_bound in bound:
        if max_bound < min_bound:
            raise Exception(f'抱歉,({min_bound}, {max_bound})不是合格的取值范圍')

    #所有變量的最小值和最大值
    #var_len為所有變量的取值范圍大小
    #bit為每個變量按整數(shù)編碼最小的二進制位數(shù)
    min_nums, max_nums = np.array(list(zip(*bound)))
    self.var_len = var_len = max_nums-min_nums
    bits = np.ceil(np.log2(var_len+1))

    if DNA_SIZE == None:
        DNA_SIZE = int(np.max(bits))
    self.DNA_SIZE = DNA_SIZE

    #POP_SIZE為進化的種群數(shù)
    self.POP_SIZE = len(nums)
    POP = np.zeros((*nums.shape, DNA_SIZE))
    for i in range(nums.shape[0]):
        for j in range(nums.shape[1]):
    #編碼方式:
            num = int(round((nums[i,j] - bound[j][0]) * ((2**DNA_SIZE) / var_len[j])))
    #用python自帶的格式化轉化為前面空0的二進制字符串,然后拆分成列表
            POP[i,j] = [int(k) for k in ('{0:0'+str(DNA_SIZE)+'b}').format(num)]
    self.POP = POP
    #用于后面重置(reset)
    self.copy_POP = POP.copy()
    self.cross_rate = cross_rate
    self.mutation = mutation
    self.func = func

#save args對象保留參數(shù):
#        bound                取值范圍
#        var_len              取值范圍大小
#        POP_SIZE             種群大小
#        POP                  編碼后的種群[[[1,0,1,...],[1,1,0,...],...]]
#                             一維元素是各個種群,二維元素是各個DNA[1,0,1,0],三維元素是堿基對1/0
#        copy_POP             復制的種群,用于重置
#        cross_rate           染色體交換概率
#        mutation             基因突變概率
#        func                 適應度函數(shù)

3.translateDNA分析及實現(xiàn)

??作為翻譯函數(shù),我們希望它能把被編碼過的種群POP翻譯(解碼)回來
??方法:先把0101列表轉化為10進制數(shù),再通過解碼公式解碼

def translateDNA(self):
    W_vector = np.array([2**i for i in range(self.DNA_SIZE)]).reshape((self.DNA_SIZE, 1))[::-1]
    binary_vector = self.POP.dot(W_vector).reshape(self.POP.shape[0:2])
    for i in range(binary_vector.shape[0]):
        for j in range(binary_vector.shape[1]):
            binary_vector[i, j] /= ((2**self.DNA_SIZE)/self.var_len[j])
            binary_vector[i, j] += self.bound[j][0]
    return binary_vector

??首先W_vector是帶權向量,形如[..., 256, 128, 64, 32, 16, 8, 4, 2, 1],類似二進制轉十進制的過程,將二進制對應位數(shù)累加即可。
??可以簡化為矩陣乘法,首先,POP的大小(shape)為(POP_SIZE, VAR_QUANTITY, DNA_SIZE),假設有100個種群,每個種群有3個DNA(變量數(shù)量),每個DNA有20個堿基對,那么對應的W_vector是
\left[\begin{matrix} 2^{19} \\ 2^{18} \\ 2^{17} \\ …… \\ 2^{2} \\ 2^{1} \\ 2^{0} \end{matrix}\right]
形狀為(DNA_SIZE, 1),與POP相乘得到的矩陣為(POP_SIZE, VAR_QUANTITY, 1),最后一維只有一個數(shù)字顯然不是我們想要的,用reshape((POP_SIZE, VAR_QUANTITY))可以將最后一維去掉,我們也回到了最開始(輸入)的維度。
??最終,返回翻譯回來的矩陣

3. get_fitness分析及實現(xiàn)

??我們希望通過get_fitness得到當前的適應值,有可能懷有疑問,直接將translateDNA的結果交給適應函數(shù)func函數(shù)不就是適應值嗎?的確這樣很簡單,但是有時候我們需要對適應值進行進一步處理,因此,我們封裝一層函數(shù),易于在有需求時修改代碼。

def get_fitness(self, non_negative=False):
    result = self.func(*np.array(list(zip(*self.translateDNA()))))
    if non_negative:
        min_fit = np.min(result, axis=0)
        result -= min_fit 
    return result

??我們在后面看到一個需求,就是有時候我們需要非負的適應值,因此我們加了一個帶默認值參數(shù)non_negative,假如需要非負適應值時,只要傳入True就可以了
??首先將含有一堆[x0, x1, ..]變量的列表轉置,變?yōu)楹衃POP0_x0, POP1_x1, ..]這種形式,將其解包后傳給func,得到的result是一維行向量。
??如果要求非負,只要令每個適應值減去最小值即可。結果返回適應值向量
??注意:本實例的適應度函數(shù)func不需要判斷(if),如果適應度函數(shù)中有判斷(分段),此處需要遍歷數(shù)組而不能直接傳入func進行矩陣運算。

4. select分析及實現(xiàn)

??select意為自然選擇,適者生存,不適者被淘汰,因此適應度低的種群會有更大概率消失,當然這也意味著,就算是有再高的適應度,也是有幾率被淘汰的(除非所有種群已經(jīng)達到統(tǒng)一),我們希望這個函數(shù)能實現(xiàn)一次對內部成員POP的自然選擇

def select(self):
    fitness = self.get_fitness(non_negative=True)
    self.POP = self.POP[np.random.choice(np.arange(self.POP.shape[0]), size=self.POP.shape[0], replace=True, p=fitness/np.sum(fitness))]

??首先,我們通過get_fitness獲得適應度向量。
??然后是一大串程序,我們從choice開始分析。numpy.random.choice是numpy中的隨機抽取函數(shù),可以對序列中的元素進行隨機抽取,第一個參數(shù)為待抽取序列,size參數(shù)為抽取的數(shù)量,replace參數(shù)為是否允許重復,p一個列表,列表中的每個元素為第一個參數(shù)序列中每個參數(shù)被選中的概率
??因此我們并非是將原有的元素按概率進行剔除,而是類似新抽取一個列表,適應度大的種群將有更大的幾率存活并在生態(tài)圈(POP)中占有更大的比重
??在程序中我們不是直接抽取POP,而是對POP進行編號(np.arange),對編號序列進行抽取,生態(tài)圈大小維持原來的大小,允許重復,否則自然選擇不會起作用,而每個種群的概率是:\frac{fit_{current}}{\sum_{i=1}^{n}{fit_i}}將所有種群組成的列表傳給p,得到的是優(yōu)勝者的編號,類似[12, 3, 45, 23, 12, 12,...]這樣。
??將返回的矩陣傳給原來的POP,POP為numpy.array,給item方法(等同于numpy.array對象取[ ]操作符)傳遞list<int>返回結果為列表中每個元素標號的元素,如上面的列表編號,將返回[a[12], a[3], a[45], a[23], a[12], a[12], ...]
??可能有人注意到了,第一行的get_fitness給non_negative傳遞了True,這是因為運算概率時如果有負值會出現(xiàn)負概率,而choice就會拋出:ValueError: probabilities are not non-gegative,意為概率不能為負,因此通過傳入?yún)?shù)令get_fitness返回負值

5. crossover分析及實現(xiàn)

??作為染色體交叉函數(shù),我們希望這個函數(shù)能對生態(tài)圈進行一次交叉操作。
??首先對每個種群檢測是否交換DNA序列。
??如果交換,則在生態(tài)圈中隨機找到一個與其交叉的種群。
??隨機確定交叉的位數(shù),然后進行交叉。

def crossover(self):
    for people in self.POP:
        if np.random.rand() < self.cross_rate:
            i_ = np.random.randint(0, self.POP.shape[0], size=1)
            cross_points = np.random.randint(0, 2, size=(len(self.var_len) ,self.DNA_SIZE)).astype(np.bool) 
            people[cross_points] = self.POP[i_, cross_points]

??其中i_就是提供交叉DNA種群的編號,cross_point是交叉點,最后一行是類似掩膜(mask)的操作,將交叉點的堿基對交給被交叉的種群。

6. mutation分析及實現(xiàn)

??mutation是變異函數(shù),對所有生態(tài)圈種群,找到所有DNA,檢測所有堿基,按概率變異,因此對三維列表POP進行遍歷。

def mutate(self):
    for people in self.POP:
        for var in people:
            for point in range(self.DNA_SIZE):
                if np.random.rand() < self.mutation:
                    var[point] = 1 if var[point] == 0 else 1

7. evolution分析及實現(xiàn)

??進化函數(shù)所要做的是對生態(tài)圈進行一次整體進化,解析出來,只需要進行自然選擇、染色體交叉、基因變異各一次就好,在前面都抽象好的當下,實現(xiàn)十分容易。

def evolution(self):
    self.select()
    self.crossover()
    self.mutate()

8. log分析及實現(xiàn)

??在進化的過程中,我們希望得到種群當前狀態(tài),因此需要一個日志函數(shù),在這里,我實現(xiàn)的是返回包含當前所有DNA值以及種群的適應度的Pandas.DataFrame

def log(self):
    return pd.DataFrame(np.hstack((self.translateDNA(), self.get_fitness().reshape((len(self.POP),1)))), 
                        columns=[f'x{i}' for i in range(len(self.var_len))]+['F'])
GA.log()

8. reset分析及實現(xiàn)

??將內部備份的copy_POP復制回來即可

    def reset(self):
        self.POP = self.copy_POP.copy()

9.plot_in_jupyter_1d分析及實現(xiàn)

??在jupyter-notebook中顯示多次迭代的動態(tài)圖形

    def plot_in_jupyter_1d(self, iter_time=200):
        is_ipython = 'inline' in matplotlib.get_backend()
        if is_ipython:
            from IPython import display

        plt.ion()
        for _ in range(iter_time):
            plt.cla()
            x = np.linspace(*self.bound[0], self.var_len[0]*50)
            plt.plot(x, self.func(x))
            x = self.translateDNA().reshape(self.POP_SIZE)
            plt.scatter(x, self.func(x), s=200, lw=0, c='red', alpha=0.5)
            if is_ipython:
                display.clear_output(wait=True)
                display.display(plt.gcf())

            self.evolution()

??iter_time為迭代次數(shù)

三、實驗

1.單變量實驗

??實驗目標(適應度)函數(shù)為:F(X) = X \times\sin{10·X} + X \times \cos{2·X}變量范圍為[0, 5],堿基對個數(shù)為10,生態(tài)圈初始隨機。

func = lambda x:np.sin(10*x)*x + np.cos(2*x)*x
ga = GA([[np.random.rand()*5] for _ in range(100)], [(0,5)], DNA_SIZE=10, func=func)
ga.plot_in_jupyter_1d()

可以看到點逐漸上升
開始
結果

2.多變量實驗

??多變量作圖沒有封裝,但類的抽象度很高,很容易使用,實驗所用的目標函數(shù)為:X·cos{2\pi Y+Y·\sin{2\pi·X}}生態(tài)圈初始隨機,范圍分別為[-2, 10]、[-2, 10],堿基對個數(shù)為20,交叉概率0.7,變異概率為0.01。

nums = list(zip(np.arange(-2, 2, 0.2), np.arange(-2, 2, 0.2)))
bound = [(-2, 2), (-2, 2)]
func = lambda x, y: x*np.cos(2*np.pi*y)+y*np.sin(2*np.pi*x)
DNA_SIZE = 20
cross_rate = 0.7
mutation = 0.01

ga = GA(nums=nums, bound=bound, func=func, DNA_SIZE=DNA_SIZE, cross_rate=cross_rate, mutation=mutation)
is_ipython = 'inline' in matplotlib.get_backend()
if is_ipython:
    from IPython import display

plt.ion() #打開交互式
for _ in range(200):
    plt.cla() #清空畫布
    X = Y = np.arange(-2, 2, 0.2)
    X, Y=np.meshgrid(X,Y)
    Z = func(X, Y)
    
    fig1=plt.figure()
    ax=Axes3D(fig1)
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.cm.coolwarm)#做曲面
    ax.set_xlabel('x0 label', color='r')
    ax.set_ylabel('x1 label', color='g')
    ax.set_zlabel('F label', color='b')#給三個坐標軸注明
    
    ax.scatter(*list(zip(*ga.translateDNA())), ga.get_fitness(), s=100, lw=0, c='red', alpha=0.5)
    
    if is_ipython:
        display.clear_output(wait=True)
        display.display(plt.gcf())

    ga.evolution()
開始
中間
結束

用打印日志的方式找到最大值,同理,也可以進行排序顯示:

ga.log().max()

#output:
#x0   -1.766613
#x1    1.765591
#F     3.265491
#dtype: float64

注意:每次運行結果可能有少量差異

四、后續(xù)維護

??基礎的架構已經(jīng)打完,可以在原有程序上進行改造。例如,將DNA_SIZE改成不等長方式,給每個重要變量做出get、set方法,豐富日志功能、信息,支持高緯度圖顯示等等

最后編輯于
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成,瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明:文章內容(如有圖片或視頻亦包括在內)由作者上傳并發(fā)布,文章內容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容