遗传算法求三元函数极值(python)-采用二进制编码

想看实数编码编码的
博客地址遗传算法求三元函数极值(python)-采用实数编码
*

遗传算法求三元函数极值(python)-采用二进制编码

本文的遗传算法采用二进制编码求三元函数极值
所求函数为

要想使用遗传算法，首要任务是进行编码

传统的 GA 中, DNA 我们能用一串二进制来表示, 比如:
DNA1 = [1, 1, 0, 1, 0, 0, 1]
DNA2 = [1, 0, 1, 1, 0, 1, 1]
这里，我们仍然使用二进制编码，但是如何与我们的问题对应起来呢？
为什么会要用二进制编码, 我们之后在下面的内容中详细说明这样编码的好处. 但是长成这样的 DNA 并不好使用. 如果要将它解码, 我们可以将二进制转换成十进制, 比如二进制的 11 就是十进制的 3. 这种转换的步骤在程序中很好执行. 但是有时候我们会需要精确到小数, 其实也很简单, 只要再将十进制的数浓缩一下就好. 比如我有 1111 这么长的 DNA, 我们产生的十进制数范围是 [0, 15], 而我需要的范围是 [-1, 1], 我们就将 [0, 15] 缩放到 [-1, 1] 这个范围就好.
我们知道二进制很容易转十进制，再区间压缩以下，这样一个DNA和一个解一一映射。

def translateDNA(pop):
    return pop.dot(2 ** np.arange(DNA_SIZE)[::-1]) / float(2**DNA_SIZE-1) * X_BOUND[1]

例如，1 0 1 0 1 0 0 1 0 0 => (4+32+128+256)/(1024-1)*(5-0)=3.3
上面针对的单变量的编码，求解单变量的函数极值
其完整代码如下，这里照搬了别人的代码，后面我会给出引用地址，
用它主要是参考，并修改为求解三元函数极值的代码
单变量求函数极值：

"""
Visualize Genetic Algorithm to find a maximum point in a function.
可视化遗传算法去找到一个函数的最高点
"""
import numpy as np
import matplotlib.pyplot as plt
 
DNA_SIZE = 10            # DNA length
POP_SIZE = 100           # population size，种群中个体数目
CROSS_RATE = 0.8         # mating probability (DNA crossover)，0.8的概率进行交叉配对
MUTATION_RATE = 0.003    # mutation probability，变异强度
N_GENERATIONS = 200      #迭代次数
X_BOUND = [0, 5]         # x upper and lower bounds，指定x的取值范围
 
 
def F(x):
    return np.sin(10*x)*x + np.cos(2*x)*x     # to find the maximum of this function
 
 
# find non-zero fitness for selection
#我们都需要一个评估好坏的方程, 这个方程通常被称为 fitness适应度.
#为了找到下面这个曲线当中的最高点. 那么这个 fitness 方程可以定义为高度, 越高的点, fitness 越高.
def get_fitness(pred):
    return pred + 1e-3 - np.min(pred)#因为如果直接返回pred可能是负值，而我们在计算概率的时候不能为负值。
    #要进行处理，np.min表示取最小，为最大的负数，可以使全部只变成正的；1e-3为了让float进行相除防止小数点后的数被省略
 
 
# convert binary DNA to decimal and normalize it to a range(0, 5)
#对基因的翻译，如这里函数，x轴是实数，这里解释了如何将遗传0、1序列翻译成实数。用十进制二进制转换
#pop (population)是一个储存二进制 DNA 的矩阵, 他的 shape 是这样 (pop_size, DNA_size)
#这里DNA_SIZE，X_BOUND是超参数
def translateDNA(pop):
    return pop.dot(2 ** np.arange(DNA_SIZE)[::-1]) / float(2**DNA_SIZE-1) * X_BOUND[1]
    #dot()函数是矩阵乘,而*则表示逐个元素相乘
	#np.arange()函数返回一个有终点和起点的固定步长的排列
	#pop.dot(2 ** np.arange(DNA_SIZE)[::-1])已经转换成十进制
	#但是需要归一化到0~5,如有1111这么长的DNA,要产生的十进制数范围是[0, 15], 而所需范围是 [-1, 1],就将[0,15]缩放到[-1,1]这个范围
	#a[::-1]相当于 a[-1:-len(a)-1:-1]，也就是从最后一个元素到第一个元素复制一遍。所以你看到一个倒序
	#np.arange(DNA_SIZE)[::-1]得到10,9,8,...,0
 
#这里进行优胜劣汰的选择步骤
#适者生存的 select() 很简单, 我们只要按照适应程度 fitness 来选 pop 中的 parent 就好. fitness 越大, 越有可能被选到.
def select(pop, fitness):    # nature selection wrt pop's fitness
    idx = np.random.choice(np.arange(POP_SIZE), size=POP_SIZE, replace=True,p=fitness/fitness.sum())
	#这里概率不能为负，所以pred要进行非负处理
	#replace表示抽样后是否放回，这里为True表示有放回，则可能会出现相同的索引值
    # p 就是选它的比例，按比例来选择适应度高的,也会保留一些适应度低的，因为也可能后面产生更好的变异
    #np.random.choice表示从序列中取值  np.arange()函数返回一个有终点和起点的固定步长的排列
    return pop[idx]
 
#繁衍，交叉父母的基因
def crossover(parent, pop):     # mating process (genes crossover)
    if np.random.rand() < CROSS_RATE: #这里是0.8的概率父亲会选择一个母亲进行交叉配对
        i_ = np.random.randint(0, POP_SIZE, size=1)                           #select another individual from pop选择母亲索引一个
        cross_points = np.random.randint(0, 2, size=DNA_SIZE).astype(np.bool) #得到一行[01001100]也是0、1为了选择哪些点进行交叉;然后进行布尔化
        parent[cross_points] = pop[i_, cross_points]
		#布尔型数组可以用于数组索引，布尔型数组长度必须跟被索引的轴长度一致
		#生成布尔数组可以组合应用多个布尔条件,使用&(和),|(或)之类的布尔算数运算符，python的关键字and和or在布尔型数组中无效
		#parent[cross_points]即parent列表中取出cross_points为True地方的值&&&&&！！！！
		#【母亲是pop的i_索引行DNA，选出母亲对应在cross_points为TRUE的地方的值】赋给【父亲DNA对应在cross_points选出为TRUE的地方的值】。
    return parent
 
#繁衍，有变异的基因会出现
#将某些 DNA 中的 0 变成 1, 1 变成 0
def mutate(child):
    for point in range(DNA_SIZE):
        if np.random.rand() < MUTATION_RATE:
            child[point] = 1 if child[point] == 0 else 0
    return child
 
#产生离散均匀分布的整数，若high不为None时，取[low,high)之间随机整数，否则取值[0,low)之间随机整数。
pop = np.random.randint(2, size=(POP_SIZE, DNA_SIZE))   # initialize the pop DNA
#pop = np.=random.randint(0,2,(1,DNA_SIZE).repeat(POP_SIZE,axis=0))这里是生成了一样的DNA，后面也可以随着变异变成不一样的
 
#[[01001100],
# [10111100],
# ...]
plt.ion()       # something about plotting开启图像交互模式
 
x = np.linspace(*X_BOUND, 200)
#linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)
#X_BOUND = [0, 5],要产生200个样本点
#返回固定间隔的数据。他将返回num个等间距的样本，在区间[start,stop]中。其中，区间的结束端点可以被排除在外(用endpoint标识)
plt.plot(x, F(x))
 
for _ in range(N_GENERATIONS):
    F_values = F(translateDNA(pop))    # compute function value by extracting DNA传入到F函数
 
    # something about plotting
    if 'sca' in globals(): sca.remove()
    sca = plt.scatter(translateDNA(pop), F_values, s=200, lw=0, c='red', alpha=0.5); plt.pause(0.05)#plt.pause表示显示秒数
 
    # GA part (evolution)
    fitness = get_fitness(F_values) #计算适应度
    print("Most fitted DNA: ", pop[np.argmax(fitness), :])
    pop = select(pop, fitness)#这里选出了另外一种population
    pop_copy = pop.copy()# 备个份
    for parent in pop: #这里parent为遍历pop，一次为其中一行，而这里的pop是从原pop中按适应度概率有放回的选出了POP_SIZE行
        child = crossover(parent, pop_copy)#繁衍
        child = mutate(child) #进行变异
        parent[:] = child       # parent is replaced by its child# 宝宝变大人
 
plt.ioff(); plt.show()
 
#在使用matplotlib的过程中，不能像matlab一样同时开几个窗口进行比较，可以采用交互模式，但是放在脚本里运行一闪而过，图像并不停留
#python可视化库matplotlib有两种显示模式：阻塞（block）模式&交互（interactive）模式
#在交互模式下：plt.plot(x)或plt.imshow(x)是直接出图像，不需要plt.show()
#如果在脚本中使用ion()命令开启了交互模式，没有使用ioff()关闭的话，则图像不会常留。防止这种情况，需要在plt.show()之前加上ioff()命令。
#在阻塞模式下：打开一个窗口以后必须关掉才能打开下一个新的窗口。这种情况下，默认是不能像Matlab一样同时开很多窗口进行对比的
#plt.plot(x)或plt.imshow(x)是直接出图像，需要plt.show()后才能显示图像

关于三元函数代码的修改：
我的想法是这样的，不知道对不对，不对的话，请大佬留言指出，首先要改的是种群pop产生的矩阵，我的代码是这样的：

 pop = np.random.randint(0,2, size=(POP_SIZE, DNA_SIZE*3)) #matrix (POP_SIZE, DNA_SIZE*3)

第一个参数表示产生0-2范围的随机整数，不包括2

为了说明方便，这里先假设POP_SIZE=10,DNA_SIZE=4,则DNA_SIZE*3=12，因为有三个变量x,y,z,每个变量占一个DNA_SIZE,则DNA总长度是12位，
即[010011011010],
则上面产生的种群矩阵为
[[1 0 0 0 1 0 0 1 0 0 0 0]
[0 1 1 1 1 1 1 1 1 0 1 1]
[0 0 1 0 1 1 1 1 0 1 1 0]
[1 1 0 0 1 1 1 0 1 1 0 0]
[0 0 1 1 0 1 1 1 1 0 0 1]
[1 0 0 1 0 1 1 0 1 0 0 0]
[0 0 1 0 1 1 0 1 0 0 1 0]
[1 1 1 1 0 1 1 0 0 1 1 0]
[0 1 0 0 0 0 1 0 1 1 1 1]
[1 1 1 1 1 1 1 0 1 1 0 0]]

共10行，一行表示一个二进制编码表示的DNA，
矩阵的行数为种群数目
在上面的种群矩阵中，一行表示一个个体，而每个个体是有x,y,z三个基因片段所组成的，在这里
取前DNA_SIZE=4列表示x
所以x_pop为
[[1 0 0 0]
[0 1 1 1]
[0 0 1 0]
[1 1 0 0]
[0 0 1 1]
[1 0 0 1]
[0 0 1 0]
[1 1 1 1]
[0 1 0 0]
[1 1 1 1]],
然后依次类推，中间4列表示y,最后4列表示z
初始种群的个体即使都一样也无所谓，因为后面会通过交叉，配对等会变得不一样。
在translateDNA这块修改如下，
代码如下：

def translateDNA(pop): #pop表示种群矩阵，一行表示一个二进制编码表示的DNA，矩阵的行数为种群数目
    x_pop = pop[:,0:DNA_SIZE]#取前DNA_SIZE个列表示x
    y_pop = pop[:,DNA_SIZE:DNA_SIZE*2]#取中间DNA_SIZE个列表示y
    z_pop = pop[:,DNA_SIZE*2:]#取后DNA_SIZE个列表示z
    #print(x_pop)
    #print(y_pop)
    #print(z_pop)
    #print("\n")

    '''pop:(POP_SIZE,DNA_SIZE)*(DNA_SIZE,1) --> (POP_SIZE,1)'''#二进制--->十进制
    x = x_pop.dot(2**np.arange(DNA_SIZE)[::-1])/float(2**DNA_SIZE-1)*(X_BOUND[1]-X_BOUND[0])+X_BOUND[0]
    y = y_pop.dot(2**np.arange(DNA_SIZE)[::-1])/float(2**DNA_SIZE-1)*(Y_BOUND[1]-Y_BOUND[0])+Y_BOUND[0]
    z = z_pop.dot(2**np.arange(DNA_SIZE)[::-1])/float(2**DNA_SIZE-1)*(Z_BOUND[1]-Z_BOUND[0])+Z_BOUND[0]
    #print(x,z)
    return x,y,z

有两种求解三元函数极值的代码，
第一种是按照我自己的想法，在select操作中，关于概率的修改。
下面是我的求解三元函数极值的完整代码：

# x1*x2-x1*x2+x3
import numpy as np
import random


DNA_SIZE = 24
POP_SIZE = 100
CROSSOVER_RATE = 0.8
MUTATION_RATE = 0.015
N_GENERATIONS = 100
X_BOUND = [3.0,5.0]#x1
Y_BOUND = [2.1,6.7]#x2
Z_BOUND = [1.2,9.6]#x3


def F(x, y,z):
   val=x*x-x*y+z
   #print(val)
   '''
   for index in range(len(val)):
      if val[index] <0.2:
         val[index]=0.2
         '''
   return val



def get_fitness(pop): 
    x,y ,z= translateDNA(pop)
    pred = F(x, y,z)
    return pred

def translateDNA(pop): #pop表示种群矩阵，一行表示一个二进制编码表示的DNA，矩阵的行数为种群数目
    x_pop = pop[:,0:DNA_SIZE]#取前DNA_SIZE个列表示x
    y_pop = pop[:,DNA_SIZE:DNA_SIZE*2] ##取中间DNA_SIZE个列表示y
    z_pop = pop[:,DNA_SIZE*2:]取后DNA_SIZE个列表示z
    #print(x_pop)
    #print(y_pop)
    #print(z_pop)
    #print("\n")

    '''pop:(POP_SIZE,DNA_SIZE)*(DNA_SIZE,1) --> (POP_SIZE,1)'''#二进制--->十进制
    x = x_pop.dot(2**np.arange(DNA_SIZE)[::-1])/float(2**DNA_SIZE-1)*(X_BOUND[1]-X_BOUND[0])+X_BOUND[0]
    y = y_pop.dot(2**np.arange(DNA_SIZE)[::-1])/float(2**DNA_SIZE-1)*(Y_BOUND[1]-Y_BOUND[0])+Y_BOUND[0]
    z = z_pop.dot(2**np.arange(DNA_SIZE)[::-1])/float(2**DNA_SIZE-1)*(Z_BOUND[1]-Z_BOUND[0])+Z_BOUND[0]
    #print(x,z)
    return x,y,z


def mutation(child, MUTATION_RATE=0.003):
	if np.random.rand() < MUTATION_RATE: 				#以MUTATION_RATE的概率进行变异
		mutate_point = np.random.randint(0, DNA_SIZE*3)	#随机产生一个实数，代表要变异基因的位置
		child[mutate_point] = child[mutate_point]^1 	#将变异点的二进制为反转
		
def crossover_and_mutation(pop, CROSSOVER_RATE = 0.8):
	new_pop = []
	for father in pop:		#遍历种群中的每一个个体，将该个体作为父亲
		child = father		#孩子先得到父亲的全部基因（这里我把一串二进制串的那些0，1称为基因）
		if np.random.rand() < CROSSOVER_RATE:			#产生子代时不是必然发生交叉，而是以一定的概率发生交叉
			mother = pop[np.random.randint(POP_SIZE)]	#再种群中选择另一个个体，并将该个体作为母亲
			cross_points = np.random.randint(low=0, high=DNA_SIZE*3)	#随机产生交叉的点
			child[cross_points:] = mother[cross_points:]		#孩子得到位于交叉点后的母亲的基因
		mutation(child)	#每个后代有一定的机率发生变异
		new_pop.append(child)

	return new_pop



def select(pop, fitness):    # nature selection wrt pop's fitness
    #idx = np.random.choice(np.arange(POP_SIZE), size=POP_SIZE, replace=True,
    fitnew = fitness
    for index in range(len(fitnew)):
        if fitnew[index] <0:
           fitnew[index]=0
      
    p=(fitnew)/(fitnew.sum())

    idx = np.random.choice(np.arange(POP_SIZE), size=POP_SIZE, replace=True,
                           p=p)
    #print(idx)
    return pop[idx]

def print_info(pop):
	fitness = get_fitness(pop)
	max_fitness_index = np.argmax(fitness)
	print("max_fitness:", fitness[max_fitness_index])
	x,y,z = translateDNA(pop)
	print("最优的基因型：", pop[max_fitness_index])
	print("(x, y, z):", (x[max_fitness_index], y[max_fitness_index],z[max_fitness_index]))


if __name__ == "__main__":
   pop = np.random.randint(0,2, size=(POP_SIZE, DNA_SIZE*3)) #matrix (POP_SIZE, DNA_SIZE)
   #print(pop)
   for _ in range(N_GENERATIONS):#迭代N代
      x,y,z = translateDNA(pop)
      pop = np.array(crossover_and_mutation(pop, CROSSOVER_RATE))
      fitness = get_fitness(pop)
      pop = select(pop, fitness) #选择生成新的种群

	
   print_info(pop)

第二种方法是按照网上的意思：关于函数值：因为如果直接返回pred可能是负值，而我们在计算概率的时候不能为负值。
# 要进行处理，np.min表示取最小，为最大的负数，可以使全部只变成正的；1e-3为了让float进行相除防止小数点后的数被省略
其完整代码如下：

# x1*x2-x1*x2+x3
import numpy as np
import random

DNA_SIZE = 24
POP_SIZE = 100
CROSSOVER_RATE = 0.8
MUTATION_RATE = 0.015
N_GENERATIONS = 100
X_BOUND = [3.0, 5.0]  # x1
Y_BOUND = [2.1, 6.7]  # x2
Z_BOUND = [1.2, 9.6]  # x3


def F(x, y, z):
    val = x * x - x * y + z
    return val


def get_fitness(pop):
    x, y, z = translateDNA(pop)
    pred = F(x, y, z)
    return pred + 1e-3 - np.min(pred)  # 因为如果直接返回pred可能是负值，而我们在计算概率的时候不能为负值。
    # 要进行处理，np.min表示取最小，为最大的负数，可以使全部只变成正的；1e-3为了让float进行相除防止小数点后的数被省略


def translateDNA(pop):  # pop表示种群矩阵，一行表示一个二进制编码表示的DNA，矩阵的行数为种群数目
    x_pop = pop[:, 0:DNA_SIZE]  # 取前DNA_SIZE个列表示x
    y_pop = pop[:, DNA_SIZE:DNA_SIZE * 2]  # 取中间DNA_SIZE个列表示y
    z_pop = pop[:, DNA_SIZE * 2:]  # 取后DNA_SIZE个列表示z


    '''pop:(POP_SIZE,DNA_SIZE)*(DNA_SIZE,1) --> (POP_SIZE,1)'''  # 二进制--->十进制
    x = x_pop.dot(2 ** np.arange(DNA_SIZE)[::-1]) / float(2 ** DNA_SIZE - 1) * (X_BOUND[1] - X_BOUND[0]) + X_BOUND[0]
    y = y_pop.dot(2 ** np.arange(DNA_SIZE)[::-1]) / float(2 ** DNA_SIZE - 1) * (Y_BOUND[1] - Y_BOUND[0]) + Y_BOUND[0]
    z = z_pop.dot(2 ** np.arange(DNA_SIZE)[::-1]) / float(2 ** DNA_SIZE - 1) * (Z_BOUND[1] - Z_BOUND[0]) + Z_BOUND[0]
    # print(x,z)
    return x, y, z


def mutation(child, MUTATION_RATE=0.003):
    if np.random.rand() < MUTATION_RATE:  # 以MUTATION_RATE的概率进行变异
        mutate_point = np.random.randint(0, DNA_SIZE * 3)  # 随机产生一个实数，代表要变异基因的位置
        child[mutate_point] = child[mutate_point] ^ 1  # 将变异点的二进制为反转


def crossover_and_mutation(pop, CROSSOVER_RATE=0.8):
    new_pop = []
    for father in pop:  # 遍历种群中的每一个个体，将该个体作为父亲
        child = father  # 孩子先得到父亲的全部基因（这里我把一串二进制串的那些0，1称为基因）
        if np.random.rand() < CROSSOVER_RATE:  # 产生子代时不是必然发生交叉，而是以一定的概率发生交叉
            mother = pop[np.random.randint(POP_SIZE)]  # 再种群中选择另一个个体，并将该个体作为母亲
            cross_points = np.random.randint(low=0, high=DNA_SIZE * 3)  # 随机产生交叉的点
            child[cross_points:] = mother[cross_points:]  # 孩子得到位于交叉点后的母亲的基因
        mutation(child)  # 每个后代有一定的机率发生变异
        new_pop.append(child)

    return new_pop


def select(pop, fitness):  # nature selection wrt pop's fitness
    idx = np.random.choice(np.arange(POP_SIZE), size=POP_SIZE, replace=True,
                           p=fitness/fitness.sum())
    # print(idx)
    return pop[idx]


def print_info(pop):
    fitness = get_fitness(pop)
    max_fitness_index = np.argmax(fitness)
    # print("max_fitness:", fitness[max_fitness_index])
    x, y, z = translateDNA(pop)
    print("(x, y, z):", (x[max_fitness_index], y[max_fitness_index], z[max_fitness_index]))
    x=x[max_fitness_index]
    y=y[max_fitness_index]
    z=z[max_fitness_index]
    print("函数极值：", F(x,y,z))


if __name__ == "__main__":
    pop = np.random.randint(0, 2, size=(POP_SIZE, DNA_SIZE * 3))  # matrix (POP_SIZE, DNA_SIZE)
    # print(pop)
    for _ in range(N_GENERATIONS):  # 迭代N代
        x, y, z = translateDNA(pop)
        pop = np.array(crossover_and_mutation(pop, CROSSOVER_RATE))
        fitness = get_fitness(pop)
        pop = select(pop, fitness)  # 选择生成新的种群

    print_info(pop)

参考文献：
https://blog.csdn.net/yyyxxxsss/article/details/82464736
https://morvanzhou.github.io/tutorials/machine-learning/evolutionary-algorithm/2-01-genetic-algorithm/
https://www.cnblogs.com/lfri/p/12240355.html
https://blog.csdn.net/zxxxxxxy/article/details/105287743?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.nonecase&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.nonecase