基因演算法解决函数极值问题－python版
算法 ，基因演算法，python，numpy
1.求二元函数在以下区间的近似最大值
f(x , y) = 21.5 + sin(4 * PI * x) + sin(20 * PI * y)
(3.0<=x<=12.1,4.1<=y<=5.8)
2.直接上代码
import numpy as np
import random
from math import ceil, log, pow, pi, sin
class GeneticAlogorithm:
#约束变量
a1, b1, a2, b2=3.0,12.1,4.1,5.8
#存储相关值
precison=10000
output=[0.0,0.0,0.0]
def init(self,M,T,Pc,Pm):
self.M,self.T,self.Pc,self.Pm=M,T,Pc,Pm
#码位，二进制编码长度
self.xLen=int(ceil(log((self.b1-self.a1)*self.precison,2)))
self.yLen=int(ceil(log((self.b2-self.a2)*self.precison,2)))
self.Len=self.xLen+self.yLen
#存储种群基因
self.population=np.random.random(0, 2, size=(self.M, self.Len))
self.value=np.zeros(self.M)
def decodegroup(self):
x=y=0.0
tmp_value=[]
for i in range(self.M):
tmp=0.0
for j in range(0, self.xLen):
tmp=tmp+self.population[i][j] * pow(2, j)
x = tmp * (self.b1- self.a1) / (pow(2, self.xLen) - 1) + self.a1
#
tmp=0.0
for j in range(self.xLen, self.Len):
tmp=tmp+self.population[i][j] * pow(2, j-self.xLen)
y = tmp * (self.b2- self.a2) / (pow(2, self.yLen) - 1) + self.a2
cal_value=21.5 + x * sin(4 * pi * x ) + y * sin(20 * pi * y)
tmp_value.append(cal_value)
if cal_value > self.output[2]:
self.output=[x, y, cal_value]
###
self.value=np.array(tmp_value)
def makecopy(self):
#计算所有的值的和
value_sum = np.sum(self.value)
#计算适应度
self.value = np.add.accumlate(self.value)
value_pi = self.value / value_sum
value_pi = np.add.acculate(value_pi)
#轮盘赌法
choice=np.zeros(self.M)
for i in range(self.M):
for j in range(self.M):
if random.random() < value_pi[j]:
choice[i]=j
break
self.population=self.population[choice]
def doswitch(self):
'''
单点交叉算子 交叉概率 Pc'''
tmp = np.random.permutation(self.M)
self.population[tmp]
for i in range(0, self.M, 2):
index=int(random.randint(1,self.M-1) * self.Pc)
self.population[i][-index:],self.population[i+1][-index:]=self.population[i+1][-index:],self.population[i][-index:]
def domutation(self):
'''
平均每代会有cnt个基因发生突变'''
count=int(ceil(self.Pm * self.M * self.Len))
for i in range(count):
index=random.randint(0, self.M * self.Len-1)#[a,b]
m,n=index/self.Len, index%self.Len
self.population[m][n]= 1 if self.population[m][n]==1 else 0
def run(self):
self.init(1000, 1700, 0.1345, 0.0001234)
self.decodegroup()
for i in range(self.T):
self.makecopy
self.doswitch()
self.domutation()
self.decodegroup()
print self.output
if __name__=='__main__':
hao=GeneticAlogorithm()
hao.run()3.代码解释
3.1导入相关库
import numpy as np
import random
from math import ceil, log, pow, pi, sin3.2相关约束变量
#约束变量
a1, b1, a2, b2=3.0,12.1,4.1,5.8
#保留精度的位数，最后其实没用，还是默认系统的位数，如果真的保留五位数，最后最大值没有区分度
precison=10000
#存储近似最大值，以及近似最大值时的x和y
output=[0.0,0.0,0.0]3.3相关初始化
def init(self,M,T,Pc,Pm):
self.M,self.T,self.Pc,self.Pm=M,T,Pc,Pm
#码位，二进制编码长度
self.xLen=int(ceil(log((self.b1-self.a1)*self.precison,2)))
self.yLen=int(ceil(log((self.b2-self.a2)*self.precison,2)))
self.Len=self.xLen+self.yLen
#存储种群基因
self.population=np.random.random(0, 2, size=(self.M, self.Len))
self.value=np.zeros(self.M)- self.M：种群初始大小（数量）
- self.T：种群演化的代数
- self.Pc：种群基因交差概率
- self.Pm：种群基因突变概率
- self.xLen：代表X的基因长度
- self.yLen：代表Y的基因长度
- self.Len：种群一个个体的基因长度，码长
- self.population：存储一个种群的全部基因，就是01序列
- self.value：存储一个种群的表现，怎么说，放在函数就是函数值，至于为什么要把一代群体所有的值都保存下来，直接保留一个最大值不可以吗？不可以，需要根据这一代种群所有的value来计算每个个体的适应度
至于xLen，yLen：
2^ m <= (b-a) * presion <= 2^(m+1)解释一下这个函数：
np.random.random(0, 2, size=(self.M, self.Len))
该函数返回一个self.M ＊ self.Len的二维数组，由随机的0，1填充
type(walks)
type ‘numpy.ndarray’
3.4解码并更新近似最大值
def decodegroup(self):
x=y=0.0
tmp_value=[]
for i in range(self.M):
tmp=0.0
for j in range(0, self.xLen):
tmp=tmp+self.population[i][j] * pow(2, j)
x = tmp * (self.b1- self.a1) / (pow(2, self.xLen) - 1) + self.a1
#
tmp=0.0
for j in range(self.xLen, self.Len):
tmp=tmp+self.population[i][j] * pow(2, j-self.xLen)
y = tmp * (self.b2- self.a2) / (pow(2, self.yLen) - 1) + self.a2
cal_value=21.5 + x * sin(4 * pi * x ) + y * sin(20 * pi * y)
tmp_value.append(cal_value)
if cal_value > self.output[2]:
self.output=[x, y, cal_value]
###
self.value=np.array(tmp_value)二进制转十进制这就不用说了，说说如何转换成区间的值：
x = tmp * (self.b1- self.a1) / (pow(2, self.xLen) - 1) + self.a1
基础值＋占的份数＊每一份的占的值＝在这个区间的实际值
3.5选择运算
def makecopy(self):
#计算所有的值的和
value_sum = np.sum(self.value)
#计算适应度
self.value = np.add.accumlate(self.value)
value_pi = self.value / value_sum
value_pi = np.add.acculate(value_pi)
#轮盘赌法
choice=np.zeros(self.M)
for i in range(self.M):
for j in range(self.M):
if random.random() < value_pi[j]:
choice[i]=j
break
self.population=self.population[choice]
如何选择见上面图片。采自这里
3.6交差运算
def doswitch(self):
'''
单点交叉算子 交叉概率 Pc'''
tmp = np.random.permutation(self.M)
self.population[tmp]
for i in range(0, self.M, 2):
index=int(random.randint(1,self.M-1) * self.Pc)
self.population[i][-index:],self.population[i+1][-index:]=self.population[i+1][-index:],self.population[i][-index:]- 不清楚np.random.permutation(self.M)的意义
>>> tmp = np.random.permutation(15)
>>> tmp
array([ 0, 10, 13, 12, 5, 2, 6, 9, 3, 8, 4, 14, 1, 11, 7])- 二维数组self.population根据tmp切片，打比方说，就是现在第2列装的是原来的第10列
tmp = np.random.permutation(self.M)
self.population[tmp]3.7变异运算
def domutation(self):
'''
平均每代会有cnt个基因发生突变'''
count=int(ceil(self.Pm * self.M * self.Len))
for i in range(count):
index=random.randint(0, self.M * self.Len-1)#[a,b]
m,n=index/self.Len, index%self.Len
self.population[m][n]= 1 if self.population[m][n]==1 else 0- 根据种群的总基因数和种群突变概率算出平均每一代种群的突变基因个数count
- 随机count次，让某一个体的某一位基因突变，0->1，1->0
3.8程序运行流程图
def run(self):
self.init(1000, 1700, 0.1345, 0.0001234)
self.decodegroup()
for i in range(self.T):
self.makecopy
self.doswitch()
self.domutation()
self.decodegroup()
print self.output注：程序效率不高，并没有很接近实际最大值，不过先这样吧，，也没想明白道理，毕竟python也刚入门，之所以用python写，一是：为了练习python的numpy库，二是：这个库以及python的切片索引功能来解释这个基因演算法条理清楚，代码简洁。之后会有C#版，效率更高，敬请期待。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
