利用遺傳算法求解函數極小值

思想

遺傳算法的根本思想就是達爾文的適者生存法則。

使用二進制編碼（也就是基因），對要進行優化的問題的某個屬性進行編碼。對於更適應環境的個體它有更大的概率（選擇）能夠將自己的基因遺傳給下一代（交叉）。

同時遺傳算法還允許個體的基因有一定的概率發生突變（突變），這樣可以豐富基因庫，使得可以跳出局部最優，找到全局最優。

步驟

以找二元函數\(f(x,y)\)最大值為例。

1. 確定群體數量n，並隨機生成n個個體，對於每個個體都有兩個長度為m的二進制串代表他們的DNA鏈，分別代表x、y.
2. 確定適應度函數，在這個問題中適應度函數就是\(f(x,y)\)，函數值越大代表其適應度越高.
3. 將每個個體的基因序列轉化為十進制值（二進制轉十進制），並按照比例映射到規定的x,y區間中，用適應度函數計算每個個體的適應度.
4. 選擇：通過輪盤賭的方式選出n對個體\(A_i,B_i,i \in [1,n]\)。
5. 交叉：對於每一對個體\(A_i,B_i\)，隨機選擇兩個[1,m]之間的數字\(p_1,p_2\)，作為兩個個體遺傳片段的分割點。例如新個體的x基因是由\(A_i\)的前\(p_1\)長的序列和\(B_i\)的后\(m-p_i\)長的序列遺傳得到的。
6. 變異：對於每個新產生的個體，隨機一個[0,1]的數字，若其小於變異概率，則隨機選中新個體中的一個基因中的一個單點，對其進行取反操作.
7. 用新產生的群體取代上一代的群體.
8. 重復3、4、5、6、7步操作T次，T為提前設置好的進化次數。

這樣最終群體代表的x、y就是函數\(f(x,y)\)的最大值對應坐標位置。

代碼

import numpy as np
import random

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

ITERATION_TIME = 100                # 迭代次數
POP_NUM = 100                       # 種群中個體的數量
GENE_LEN = 50                       # 基因的長度
DIMENSION = 2                       # 函數中自變量的個數
BASE = 10.0                         # 加在計算出的函數之中，防止出現函數值為負數的情況
MUTATION_RATE = 2                   # 變異概率
GENE_MAX_VALUE = 1 << GENE_LEN      # 基因能表示的最大值
RANGE = [                           # 函數自變量的范圍
    [-6.0, 3.0],
    [-3.0, 6.0]
]

def rand():
    return random.randint(0, 1333331)

class DRAW:
    def __init__(self):
        self.tot = 1

    def fitness(self, x):
        y = x[1]
        x = x[0]
        res = 3 * (1 - x) ** 2 * np.exp(-(x ** 2) - (y + 1) ** 2) - 10 * (x / 5 - x ** 3 - y ** 5) * np.exp(
            -x ** 2 - y ** 2) - 1 / 3 ** np.exp(-(x + 1) ** 2 - y ** 2)
        res += BASE
        return res

    def draw(self, points):
        if DIMENSION != 2:
            return
        fig = plt.figure()
        ax = Axes3D(fig)
        x = np.arange(RANGE[0][0], RANGE[0][1], 0.1)
        y = np.arange(RANGE[1][0], RANGE[1][1], 0.1)
        X, Y = np.meshgrid(x, y)
        Z = self.fitness([X, Y])
        ax = plt.subplot(111, projection='3d')
        ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow')
        xx, yy, zz = [], [], []
        for i in points:
            xx.append(i[0])
            yy.append(i[1])
            zz.append(self.fitness([i[0], i[1]]))
        ax.scatter(xx, yy, zz)
        ax.set_zlabel('Z')
        ax.set_ylabel('Y')
        ax.set_xlabel('X')
        plt.show()

class POPULATION(DRAW):
    def __init__(self):
        super().__init__()
        self.pop = self.getRandomPop()

    def getRandomPop(self):
        res = []
        for i in range(POP_NUM):
            t = []
            for j in range(DIMENSION):
                t.append([(rand() % 2) for k in range(GENE_LEN)])
            res.append(t)
        return res

    def BinaryToDecimal(self, b):
        res = 0
        for i in b:
            res = res * 2 + i
        return res

    def getDecimal(self):
        Decimal = []
        for unit in self.pop:
            t = []
            for i in range(DIMENSION):
                t.append(self.BinaryToDecimal(unit[i]))
            Decimal.append(t)
        return Decimal

    def decode(self, Decimal):
        DECODE = []
        for i in range(POP_NUM):
            t = []
            for j in range(DIMENSION):
                p = float(Decimal[i][j]) / GENE_MAX_VALUE
                p = (RANGE[j][1] - RANGE[j][0]) * p
                t.append(p + RANGE[j][0])
            DECODE.append(t)
        return DECODE

    def getFitness(self):
        Decimal = self.getDecimal()
        Decimal = self.decode(Decimal)
        fitness = []
        for i in Decimal:
            fitness.append(self.fitness(i))
        return fitness

    def getSelectRate(self):
        fitness = self.getFitness()
        sum = 0.0
        for i in fitness:
            sum = sum + i
        selectRate = []
        for i in fitness:
            selectRate.append(i / sum)
        return selectRate

    def select(self):
        selectRate = self.getSelectRate()
        p = np.random.choice([i for i in range(POP_NUM)], p = selectRate)
        return p

    def mutation(self, unit):
        if (rand() % 100 < MUTATION_RATE):
            r = rand() % DIMENSION
            p = rand() % GENE_LEN
            unit[r][p] = unit[r][p] ^ 1
        return unit

    def crossover(self, t1, t2):
        t1 = self.pop[t1]
        t2 = self.pop[t2]
        t = [t1, t2]
        res = []
        for i in range(DIMENSION):
            p = rand() % GENE_LEN
            tt = []
            r = rand() & 1
            for j in range(p):
                tt.append(t[r][i][j])
            for j in range(p, GENE_LEN):
                tt.append(t[r ^ 1][i][j])
            res.append(tt)
        return res


    def crossover_mutation(self):
        NEW_POP = []
        for j in range(POP_NUM):
            t1 = self.select()
            t2 = self.select()
            t = self.crossover(t1, t2)
            t = self.mutation(t)
            NEW_POP.append(t)
        self.pop = NEW_POP

    def print(self):
        Decimal = self.getDecimal()
        Decimal = self.decode(Decimal)
        res = [0.0 for i in range(DIMENSION)]
        for i in Decimal:
            for j in range(DIMENSION):
                res[j] = res[j] + i[j]
        for i in range(DIMENSION):
            res[i] = res[i] / POP_NUM
        for i in res:
            print(i, end=' ')
        print()


def main():
    pop = POPULATION()
    for i in range(ITERATION_TIME):
        pop.crossover_mutation()
        if i % 2 == 0:
            pop.draw(pop.decode(pop.getDecimal()))
            pop.print()
    pass


if __name__=='__main__':
    main()

運行結果

最一開始產生的隨機群體，他是隨機分布的。

我們目測也大概能看出極值點是在(0, 2)這附近。

最終結果，程序找到極值點在(-0.06622970029525414, 1.8151667597621541)這個位置，與我們目測值相符。