目录
前言
之前写过多目标优化的情况,算是一个最基础的版本吧,只考虑了两个变量作为输入,绝对谈不上多个输入,今天这篇在之前的基础上,将输入量增加至6个,同时对交叉、变异函数进行封装,提升代码的规范性和可读性。
基础版多目标优化代码入门:
阅读本文之前默认大家已经具备了基本遗传算法的知识了,否则不建议阅读。遗传算法的基础知识我在之前的博客中给出了说明,这里不再赘述了。
代码正文
封装函数模块,简单说明一下吧。funcl是优化函数,见下方:
length:计算变量上下限的差值,这一步是为了二进制编码转换成实数写的。
cross:染色体交叉功能。
variation:染色体变异功能
translate:解码,转换成实数值。
def funclsin(x1,x2,x3,x4,x5,x6):
z=-25*(x1-2)**2-(x2-2)**2-(x3-1)**2-(x4-4)**2-(x5-1)**2
return z
def funclcos(x1,x2,x3,x4,x5,x6):
z=x1**2+x2**2+x3**2+x4**2+x5**2+x6**2
return z
def length(lis):
result=max(lis)-min(lis)
return result
def cross(lislocal,a,b):
temp=lislocal[a+1][b]
lislocal[a+1][b]=lislocal[a][b]
lislocal[a][b]=temp
return lislocal
def variation(lislocal,a,b):
lislocal[a][b]=(not lislocal[a][b])+0 #取反操作,python和matlab不一样
return lislocal
def translate(xlist,Ulist,ll):
m=0
x=0
for k in range(ll):
m=Ulist[k]*(2**(k-1))+m
x=xlist[0]+m*length(xlist)/(2**(ll-1))
return x
多目标优化过程,思路跟之前是一样的,不过是增加了输入量,所以代码看上去更多更复杂了,这也是为啥要强调代码封装的重要性。
import numpy as np
import matplotlib.pyplot as plt
import math as mt
#%%
pi=mt.pi
NP=50 #初始化种群数
L=15 #二进制位串长度
Pc=0.9 #交叉率
Pm=0.3 #变异率
G=20 #最大遗传代数
x1=[0,10]
x2=[0,10]
x6=[0,10]
x3=[1,5]
x5=[1,5]
x4=[0,6]
f11=np.random.randint(0,high=2,size=(NP,L)) #生成随机初始种群
f22=np.random.randint(0,high=2,size=(NP,L)) #生成随机初始种群,这里把xy分开生成
f33=np.random.randint(0,high=2,size=(NP,L))
f44=np.random.randint(0,high=2,size=(NP,L))
f55=np.random.randint(0,high=2,size=(NP,L))
f66=np.random.randint(0,high=2,size=(NP,L))
x11=np.zeros(100).tolist() #记录种群生成的x1~x6的十进制数值
x22=np.zeros(100).tolist()
x33=np.zeros(100).tolist()
x44=np.zeros(100).tolist()
x55=np.zeros(100).tolist()
x66=np.zeros(100).tolist()
tempx1=np.zeros(100).tolist()
tempx2=np.zeros(100).tolist()
tempx3=np.zeros(100).tolist()
tempx4=np.zeros(100).tolist()
tempx5=np.zeros(100).tolist()
tempx6=np.zeros(100).tolist()
f1trace=[] #记录第一个目标函数
f2trace=[] #记录第二个目标函数
diedazhongzhi=-1
#%%遗传算法循环
for i in range(G):
print(i)
f1=[]
f2=[]
nf=f11 #出新的种群,在其基础上进行交叉、变异
nf1=f22
nf2=f33
nf3=f44
nf4=f55
nf5=f66
for M in range(0,NP,2):
p=np.random.rand() #交叉
if p<Pc:
q=np.random.randint(0,high=2,size=(1,L))[0].tolist()
for j in range(L):
if q[j]==1:
nf=cross(nf,M,j)
nf1=cross(nf1,M,j)
nf2=cross(nf2,M,j)
nf3=cross(nf3,M,j)
nf4=cross(nf4,M,j)
nf5=cross(nf5,M,j)
j=1
while j<=(NP*Pm):
h=np.random.randint(1,high=NP) #变异
for k in range(round(L*Pm)):
g=np.random.randint(1,high=L)
nf=variation(nf,h,g)
nf1=variation(nf1,h,g)
nf2=variation(nf2,h,g)
nf3=variation(nf3,h,g)
nf4=variation(nf4,h,g)
nf5=variation(nf5,h,g)
j+=1
#交叉变异结束之后,新一代nf加上前代f共2*NP个个体进行筛选
newf=np.vstack((f11,nf)) #垂直方向累加两个f
newf1=np.vstack((f22,nf1))
newf2=np.vstack((f33,nf2))
newf3=np.vstack((f44,nf3))
newf4=np.vstack((f55,nf4))
newf5=np.vstack((f66,nf5))
for j in range(newf.shape[0]):
U=newf[j]
U1=newf1[j]
U2=newf2[j]
U3=newf3[j]
U4=newf4[j]
U5=newf5[j]
x11[j]=translate(x1,U,L) #将二进制解码为定义域范围内十进制
x22[j]=translate(x2,U1,L)
x33[j]=translate(x3,U2,L)
x44[j]=translate(x4,U3,L)
x55[j]=translate(x5,U4,L)
x66[j]=translate(x6,U5,L)
f1.append(funclsin(x11[j],x22[j],x33[j],x44[j],x55[j],x66[j]))
f2.append(funclcos(x11[j],x22[j],x33[j],x44[j],x55[j],x66[j]))
fs=[]
for j in range(len(f1)):
fs.append(f1[j])
fs.append(f2[j])
fs=np.array(fs)
fs=fs.reshape(len(f1),2) #把求解得到的适应度函数转换为多目标散点的形式
fs=fs.tolist()
ps=np.zeros((len(f1),len(f1))) #i,j=1,i支配j,我们希望点的值越小越好
for k in range(0,len(f1)):
for j in range(0,len(f1)):
if fs[k][0]<fs[j][0] and fs[k][1]<fs[j][1]:
ps[k][j]=1
elif fs[k][0]==fs[j][0] and fs[k][1]<fs[j][1]:
pass
elif fs[k][0]<fs[j][0] and fs[k][1]==fs[j][1]:
ps[k][j]=1
#这一步是这样的,ps是一个比较矩阵,如果100个点互相比(包括自己)
#如果被置1了,就说明被支配了,全部比完之后,纵向全为0,则说明没有点支配自己
jishu=[]
for k in range(len(ps)):
aa=0
for j in range(len(ps)):
if ps[j][k]==1:
aa+=1
jishu.append(aa)
if max(jishu)==0 and diedazhongzhi==-1:
print('迭代结束于{}代'.format(i))
diedazhongzhi=i+6
elif i==diedazhongzhi:
break
sortjishu=np.sort(jishu) #升序排序
index= np.argsort(jishu) #升序对应索引
aaa=list(set(sortjishu))
front=[[] for k in range(len(set(sortjishu)))]
for k in range(len(set(sortjishu))):
for j in range(len(index)):
if sortjishu[j]==aaa[k]:
front[k].append(index[j]) #记录对应的索引位置
#只是依靠分层的话,很可能最后会优化到一个点上
#所以需要加入拥挤距离来刻画点的分布情况
#f11到f66是记录考虑拥挤距离之后的染色体,在程序的末端部分添加。
f11=[]
f22=[]
f33=[]
f44=[]
f55=[]
f66=[]
for j in range(len(front)):
if len(f11)==NP:
break
tempf11=[]
tempf22=[]
tempf33=[]
tempf44=[]
tempf55=[]
tempf66=[]
for k in range(len(front[j])):
tempf11.append(newf[front[j][k]]) #把第J层的染色体添加进入tempf中
tempf22.append(newf1[front[j][k]])
tempf33.append(newf2[front[j][k]])
tempf44.append(newf3[front[j][k]])
tempf55.append(newf4[front[j][k]])
tempf66.append(newf5[front[j][k]])
tempf1=[]
tempf2=[]
for M in range(len(tempf11)): #计算第j层染色体的函数值,存在temf1.2中
U=newf[M]
U1=newf1[M]
U2=newf2[M]
U3=newf3[M]
U4=newf4[M]
U5=newf5[M]
tempx1[M]=translate(x1,U,L)
tempx2[M]=translate(x2,U,L)
tempx3[M]=translate(x3,U,L)
tempx4[M]=translate(x4,U,L)
tempx5[M]=translate(x5,U,L)
tempx6[M]=translate(x6,U,L)
tempf1.append(funclsin(tempx1[M],tempx2[M],tempx3[M],tempx4[M],tempx5[M],tempx6[M]))
tempf2.append(funclcos(tempx1[M],tempx2[M],tempx3[M],tempx4[M],tempx5[M],tempx6[M]))
#对f1计算其拥挤距离,此时与f2无关。
tempf1index=np.argsort(tempf1) #tempf1升序排序
distance=np.zeros(len(tempf1))
#将两端位置值置为正负无穷
#导致收敛到两端上,我现在只取两个极值点,避免重复
distance[tempf1index[0]]=float('inf')
distance[tempf1index[-1]]=float('inf')
#计算每个点在f1下的拥挤距离
for k in range(len(tempf1index)):
if abs(distance[tempf1index[k]])<10: #把极值点排除在外,极值点如果不排除呢?
distance[tempf1index[k]]=(tempf1[tempf1index[k+1]]-tempf1[tempf1index[k-1]])/(max(tempf1)-min(tempf1))
else:
continue
#对f2计算拥挤距离,与f1无关
tempf1index1=np.argsort(tempf2) #tempf2升序排序
distance1=np.zeros(len(tempf1))
#将两端位置值置为正负无穷
distance1[tempf1index1[0]]=float('inf')
distance1[tempf1index1[-1]]=float('inf')
for k in range(len(tempf1index1)):
if abs(distance1[tempf1index1[k]])<10: #把极值点排除在外
distance1[tempf1index1[k]]=(tempf2[tempf1index1[k+1]]-tempf2[tempf1index1[k-1]])/(max(tempf2)-min(tempf2))
else:
continue
sumdis=[] #当收敛到一个点时,拥挤距离就失去意义了
for k in range(len(distance)):
sundistance=distance[k]+distance1[k]
sumdis.append(sundistance)
disindex=np.argsort(sumdis)[::-1] #sumdis降序排序(拥挤距离越大越好)
for k in range(len(sumdis)):
f11.append(tempf11[disindex[k]])
f22.append(tempf22[disindex[k]])
f33.append(tempf33[disindex[k]])
f44.append(tempf44[disindex[k]])
f55.append(tempf55[disindex[k]])
f66.append(tempf66[disindex[k]])
if len(f11)==NP:
break
else:
continue
f1=[]
f2=[]
for j in range(len(f11)):
U=f11[j]
U1=f22[j]
U2=f33[j]
U3=f44[j]
U4=f55[j]
U5=f66[j]
x11[j]=translate(x1,U,L) #将二进制解码为定义域范围内十进制
x22[j]=translate(x2,U1,L)
x33[j]=translate(x3,U2,L)
x44[j]=translate(x4,U3,L)
x55[j]=translate(x5,U4,L)
x66[j]=translate(x6,U5,L)
f1.append(funclsin(x11[j],x22[j],x33[j],x44[j],x55[j],x66[j]))
f2.append(funclcos(x11[j],x22[j],x33[j],x44[j],x55[j],x66[j]))
plt.scatter(f1,f2)
plt.savefig(r'C:\Users\sc35131\Desktop\Genetic Algorithm\pic\\'+str(i),bbox_inches = 'tight',pad_inches = 0,dpi =100)
plt.close()
结果展示
给出第一代和第20代的对比图,我们设计的目标是希望优化函数越小越好。
第一代
第20代
优化效果还是比较明显的。
写在最后
这个脚本功能已经可以满足大多数比较简单的工程实际优化问题了,但是,可提升的空间还是很大的,具体如下:
- 代码的封装,更进一步,完全可以用class来写同样的功能,那时候代码就会简洁许多了,这也是我的小小目标吧,之后一定会把相同的功能再次复现一遍。(可能要等上一段时间了)
- 约束条件,我们现在的遗传算法就是给出了输入量的区间范围,然后不断地取交叉变异,选出符合我们函数目标的个体,但是,在实际问题中,其实是有很多的约束条件的。比如同样对于我们举例的这个函数,它的输入量之间就可以加入一些约束条件。
那如何在遗传算法中加入这些约束条件,就是后面我们需要讨论和解决的问题了。
假期结束了,继续加油。