遗传算法求极值(英文版,GA for Extremum)

最新推荐文章于 2022-06-25 11:04:28 发布

原创最新推荐文章于 2022-06-25 11:04:28 发布 · 833 阅读

5 ·

CC 4.0 BY-SA版权

文章标签：

#python #算法

研究生课程专栏收录该内容

13 篇文章

订阅专栏

本文介绍了遗传算法（GA），它由美国的John Holland在20世纪70年代提出，可用于解决优化问题，如今也应用于机器学习等领域。通过Python代码演示了用GA求解函数最小值，包括编码种群、解码、选择新种群、交叉与变异等步骤，最终能在短时间内找到高精度的全局最优解。

Introduction and thinking of GA

GA(Genetic Algorithm) was proposed by John Holland of America in 1970s. GA is simulated natural selection as Reproduction, Crossover, Mutation in DNA. When we calculate the complex questions of optimization, GA will be used successfully in short time.

Nowadays, GA can also be uesd in machine learning, signal processing, adaptive control.
\section{selecting the function, optimizing by GA}

The codes are consisted with main function and some subfunction.This is a simple example to demonstrate finding the minimum of function with Python. To simplify the assumption, we can use the simple function:
\begin{equation}
y=\sin x + 0.01 \exp ^x
\end{equation}
\subsection{length of chromosome}
Input the $p r e c i s i o n a n d d e f i n i t i o n d o m a i n$ (精度和定义域),then we can get the length of chromosomes. Here we use binary encoding.

# 根据解的精度确定染色体(chromosome)的长度
# 需要根据决策变量的上下边界来确定
def getEncodedLength(delta=0.0001, boundarylist=[]):
    # 每个变量的编码长度
    lengths = []
    for i in boundarylist:    # i在第一个循环取boundarylist的第一行1和2列，第二次取第二行1和2列
        lower = i[0]
        upper = i[1]
        # lamnda 代表匿名函数f(x)=0,50代表搜索的初始解
        res = fsolve(lambda x: ((upper - lower) * 1 / delta) - 2 ** x - 1, 50)
        length = int(np.floor(res[0]))
        lengths.append(length)
    return lengths
    pass

Code the population

Definite the population size is 10, so the number of chromosomes is 10. The function of these codes is inputing length and quantity, outputing initial chromosomes.

def getIntialPopulation(encodelength, populationSize):
    # 随机化初始种群为0
    chromosomes = np.zeros((populationSize, sum(encodelength)), dtype=np.uint8)
    for i in range(populationSize):
        chromosomes[i, :] = np.random.randint(0, 2, sum(encodelength))
    return chromosomes

Decoded

After decoding the genome generated by the above code, the function independent variable can be obtained:

def decodedChromosome(encodelength, chromosomes, boundarylist, delta=0.0001):
    populations = chromosomes.shape[0]
    variables = len(encodelength)
    decodedvalues = np.zeros((populations, variables))
    for k, chromosome in enumerate(chromosomes):
        chromosome = chromosome.tolist()
        start = 0
        for index, length in enumerate(encodelength):
            # 将一个染色体进行拆分，得到染色体片段
            power = length - 1
            demical = 0
            for i in range(start, length + start):
                demical += chromosome[i] * (2 ** power)
                power -= 1
            lower = boundarylist[index][0]
            upper = boundarylist[index][1]
            decodedvalue = lower + demical * (upper - lower) / (2 ** length - 1)
            decodedvalues[k, index] = decodedvalue
            start = length # 开始去下一段染色体的编码
    return decodedvalues

Finding the mathematic function

def getFitnessValue(func, chromosomesdecoded):
    population, nums = chromosomesdecoded.shape
    fitnessvalues = np.zeros((population, 1)) # 初始化种群的适应度值为0
    for i in range(population):  # 计算适应度值
        fitnessvalues[i, 0] = func(chromosomesdecoded[i, :])
    # 计算每个染色体被选择的概率
    probability = fitnessvalues / np.sum(fitnessvalues)
    # 得到每个染色体被选中的累积概率
    cum_probability = np.cumsum(probability)
    return fitnessvalues, cum_probability

Select new population

Input the code of chromosomes and cum_probability，random a number in 0 $∼\sim$ 1, compared with the cum_probability. Preserve, reproduct and outpute the chromosome if the cum>random.

def selectNewPopulation(chromosomes, cum_probability):
    m, n = chromosomes.shape
    newpopulation = np.zeros((m, n), dtype=np.uint8)
    # 随机产生M个概率值
    randoms = np.random.rand(m)
    for i, randoma in enumerate(randoms):
        logical = cum_probability >= randoma
        index = np.where(logical == 1)   # 寻找logical里面等于1的项，将其保留
        newpopulation[i, :] = chromosomes[index[0][0], :]
    return newpopulation
    pass

Crossover& Mutation

Crossover and mutation are not necessary. There are probabilities of them:
$P_c = 0.4 \sim 0.99\\ Mutation probability P_m = 0.0001 \sim 0.1$

The crossover probability is the probability of individual crossover. Mutation probability is the probability of mutation of each gene (encoded in binary, 1 to 0 and o to 1). Considered there are too many genes in a population, probability $×\times$ sample can be thought as the number of each crossover or mutation.

def crossover(population, Pc=0.8): # 新种群交叉
    # 根据交叉概率计算需要进行交叉的个体个数
    m, n = population.shape
    numbers = np.uint8(m * Pc)
    # 确保进行交叉的染色体个数是偶数个
    if numbers % 2 != 0:
        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, :] = population[i, :]
    # crossover
    while len(index) > 0:
        a = index.pop()
        b = index.pop()
        # 随机产生一个交叉点
        crossoverPoint = random.sample(range(1, n), 1)
        crossoverPoint = crossoverPoint[0]
        # one-single-point crossover
        updatepopulation[a, 0:crossoverPoint] = population[a, 0:crossoverPoint]
        updatepopulation[a, crossoverPoint:] = population[b, crossoverPoint:]
        updatepopulation[b, 0:crossoverPoint] = population[b, 0:crossoverPoint]
        updatepopulation[b, crossoverPoint:] = population[a, crossoverPoint:]
    return updatepopulation
    pass
    
    def mutation(population, Pm=0.005): # 染色体变异
        updatepopulation = np.copy(population)
        m, n = population.shape
        # 计算需要变异的基因个数
        gene_num = np.uint8(m * n * Pm)
        mutationGeneIndex = random.sample(range(0, m * n), gene_num)
        # 确定每个将要变异的基因在整个染色体中的基因座(即基因的具体位置)
        for gene in mutationGeneIndex:
            # 确定变异基因位于第几个染色体
            chromosomeIndex = gene // n
            # 确定变异基因位于当前染色体的第几个基因位
            geneIndex = gene % n
            if updatepopulation[chromosomeIndex, geneIndex] == 0:
                updatepopulation[chromosomeIndex, geneIndex] = 1
            else:
                updatepopulation[chromosomeIndex, geneIndex] = 0
        return updatepopulation
        pass

Fitness

def fitnessFunction():
    return lambda x: - np.sin(x[0]) - 0.01 * np.exp(x[1])
    pass

Main

def main(max_iter=5000):
    optimalSolutions = []
    optimalValues = []
    decisionVariables = [[-5, 5], [-5, 5]] # 决策变量的取值范围
    # 得到染色体编码长度
    lengthEncode = getEncodedLength(boundarylist=decisionVariables)
    for iteration in range(max_iter):
        # 得到初始种群编码
        chromosomesEncoded = getIntialPopulation(lengthEncode, 10)
        # 种群解码
        decoded = decodedChromosome(lengthEncode, chromosomesEncoded, decisionVariables)
        # 得到个体适应度值和个体的累积概率
        evalvalues, cum_proba = getFitnessValue(fitnessFunction(), decoded)
        # 选择新的种群
        newpopulations = selectNewPopulation(chromosomesEncoded, cum_proba)
        # 进行交叉操作
        crossoverpopulation = crossover(newpopulations)
        # 进行变异操作
        mutationpopulation = mutation(crossoverpopulation)
        # 将变异后的种群解码，得到每轮迭代最终的种群
        final_decoded = decodedChromosome(lengthEncode, mutationpopulation, decisionVariables)
        # 适应度评价
        fitnessvalues, cum_individual_proba = getFitnessValue(fitnessFunction(), final_decoded)
        # 搜索每次迭代的最优解，以及最优解对应的目标函数的取值
        optimalValues.append(np.max(list(fitnessvalues)))
        index = np.where(fitnessvalues == max(list(fitnessvalues)))
        optimalSolutions.append(final_decoded[index[0][0], :])
    # 搜索最优解
    optimalValue = np.max(optimalValues)
    optimalIndex = np.where(optimalValues == optimalValue)
    optimalSolution = optimalSolutions[optimalIndex[0][0]]
    return optimalSolution, optimalValue

Finaly:

solution, value = main()
print('最优解: x1, x2')
print(solution[0], solution[1])
#Measuring run time
elapsedtime = timeit.timeit(stmt=main, number=1)
print('Searching Time Elapsed:(S)', elapsedtime)

run, the solution can be got： $x_{1}=4.719386587319752, x_{2}=-4.731898985275044$ ，minf(x)： $- 0.9998874195924371 ， S e a r c h i n g T i m e E l a p s e d : (S) 3.7949634000000003$ 。
The actual solution of minf(x) is about $- 0.9999326$ ，as $x_{1}=4.71239$ ,$ x_{2}=-5$. Compared with GA, it can find the global optimal solution with high precision and short time (about 3.79s) after 5000 $i t e r a t i o n s$ (迭代)