差分进化算法原理
差分进化算法是一种随机的启发式搜索算法,简单易用,有较强的鲁棒性和全局搜索能力。
差分进化算法是一种自组织最小化方法,利用种群中随机选择的不同向量来干扰一个现有向量,种群中的每个向量均要受到干扰,其中种群干扰向量可独立进行的,由此说明,其进行是并行的。
差分进化算法一共有五个步骤:
- 初始化
- 变异
- 交叉
- 选择
- 边界条件处理
通过上述流程可获得一个收敛性非常好的结果,引导搜索结果向全局最优解逼近。
差分进化算法的流程详解:
初始化
差分进化算法利用NP个维数为D的实数值参数向量,将它们作为每一代的种群,个体表示为:
x
i
,
G
(
i
=
1
,
2
,
3
,
⋯
,
N
P
)
x_{i,G}(i=1,2,3,\cdots,NP)
xi,G(i=1,2,3,⋯,NP)
NP为种群规模,
i
i
i表示在种群中的序列,G表示进化代数,在此过程中种群规模始终保持不变。
通常 初始化种群是在给定约束条件内随机生成,一般假定所有随机生成种群符合均匀概率分布。如果可以预先得到问题的初步解,则可通过对初步解加入正态分布随机偏差来产生,可提高重建效率。
初始化式子为
x
j
i
,
0
=
r
a
n
d
[
0
,
1
]
×
(
x
j
U
−
x
j
L
)
+
x
j
L
x_{ji,0}=rand[0,1]\times (x_{j}^{U}-x_{j}^{L})+x_{j}^{L}
xji,0=rand[0,1]×(xjU−xjL)+xjL
x
j
L
x_{j}^{L}
xjL为下界限,
x
j
U
x_{j}^{U}
xjU为上界限。
代码如下:
#给定初始值
NP=50 #种群规模
D=2 #向量的维数
G=100 #最大进化代数
F=0.5 #变异算子
CR=0.1 #交叉算子
Xs=4 #上限
Xx=-4 #下限
################种群初始化##########
x=np.zeros((D,NP)) #初始化种群
v=np.zeros((D,NP)) #初始化变异种群
u=np.zeros((D,NP)) #初始化选择种群
x=np.random.rand(D,NP)*(Xs-Xx)+Xx #产生Xx~Xs间的随机数
变异
在基本差分进化算法中目标向量
x
i
,
G
x_{i,G}
xi,G由下式产生:
V
i
,
G
+
1
=
x
r
1
,
G
+
F
(
x
r
2
,
G
−
x
r
3
,
G
)
V_{i,G+1}=x_{r1,G}+F(x_{r2,G}-x_{r3,G})
Vi,G+1=xr1,G+F(xr2,G−xr3,G)
式子中
r
1
,
r
2
,
r
3
,
i
r1,r2,r3,i
r1,r2,r3,i 为
[
1
,
N
P
]
[1,NP]
[1,NP]上四个不同的序号。
变异算子
F
F
F为
[
0
,
2
]
[0,2]
[0,2]一个实常数因数,作用是起到控制偏差变量的放大作用
代码如下:
#############变异操作##########
for m in range(0,NP):
r1=random.randint(0,NP-1)
while (r1==m):
r1 = random.randint(0, NP-1)
r2 = random.randint(0, NP-1)
while (r2 == m) and (r2 == r1):
r2 = random.randint(0, NP-1)
r3 = random.randint(0, NP-1)
while (r3 == m) and (r3==r2) and (r3==r1):
r3 = random.randint(0,NP-1)
#差分
v[:,m]=x[:,r1]+F*(x[:,r2]-x[:,r3])
交叉
交叉为了增加干扰参数向量
U
i
,
G
+
1
U_{i,G+1}
Ui,G+1的多样性,交叉算子起到的作用类似一个阀值的作用。数学式子不便于看,我们直接看这部分的python代码
###########交叉操作############
r=np.random.randint(1,3)
for n in range(0,D):
cr=np.random.rand()
if (cr<=CR and n==r):
u[n,:]=v[n,:]
else:
u[n,:]=x[n,:]
cr即为交叉算子。
选择
为决定实验向量是否会成为下一代中的成员,按照贪婪策略将实验向量
U
i
,
G
+
1
U_{i,G+1}
Ui,G+1与目标向量
x
i
,
G
x_{i,G}
xi,G进行比较。注意比较时只和一个进行比较,并不是和种群中所有的个体进行比较。
Ob1=[]
for m in range(0,NP):
Ob1.append(fun_fitness(u[:,m]))
for m in range(0,NP):
if Ob1[m]<Ob[m]:
x[:,m]=u[:,m]
for m in range(0,NP):
Ob[m]=fun_fitness(x[:,m])
trace.append(min(Ob))
sort_Ob=sorted(Ob)
sort_Ob_index=get_index(Ob,sort_Ob)
x=x[:,sort_Ob_index]
X=x[:,0]
Y=min(Ob)
边界条件处理
边界条件处理有多种,比如可以将不符合的新个体用随机产生的参数向量代替,又比如可以进行边界吸收处理,直接将超过的设置为临界值。
此处展示代码使用方式为边界吸收处理。
#######################边界条件处理##############
for n in range(0,D):
for m in range(0,NP):
if u[n,m]<Xx:
u[n,m]=Xx
if u[n,m]>Xs:
u[n,m]=Xs
实例
求函数
f
(
x
,
y
)
=
3
c
o
s
(
x
y
)
+
x
+
y
f(x,y)=3cos(xy)+x+y
f(x,y)=3cos(xy)+x+y的最小值。
函数图像为:
求解python代码附上:
import numpy as np
from matplotlib import pyplot as plt
import random
from mpl_toolkits.mplot3d import Axes3D
#定义坐标轴
fig = plt.figure()
ax1 = plt.axes(projection='3d')
xx = np.arange(-4,4,0.02)
yy = np.arange(-4,4,0.02)
X, Y = np.meshgrid(xx, yy)
Z = 3*np.cos(X*Y)+X+Y
#作图
ax1.plot_surface(X,Y,Z,cmap='rainbow')
#ax3.contour(X,Y,Z, zdim='z',offset=-2,cmap='rainbow) #等高线图,要设置offset,为Z的最小值
plt.xlabel("X-axis")
plt.ylabel("Y-axis")
plt.title("3D image")
plt.show()
################求解过程############
#给定初始值
NP=50 #种群规模
D=2 #向量的维数
G=100 #最大进化代数
F=0.5 #变异算子
CR=0.1 #交叉算子
Xs=4 #上限
Xx=-4 #下限
################种群初始化##########
x=np.zeros((D,NP)) #初始化种群
v=np.zeros((D,NP)) #初始化变异种群
u=np.zeros((D,NP)) #初始化选择种群
x=np.random.rand(D,NP)*(Xs-Xx)+Xx #产生Xx~Xs间的随机数
#适应度函数
def fun_fitness(x):
value=3*np.cos(x[0]*x[1])+x[0]+x[1]
return value
def get_index(Ob,sort_Ob):
sort_Ob_index=[]
for i in range(0,len(Ob)):
qt=Ob.index(sort_Ob[i])
sort_Ob_index.append(qt)
return sort_Ob_index
#计算初始化种群对应的适应度函数值
Ob=[]
for m in range(0,NP):
Ob.append(fun_fitness(x[:,[m]]))
trace=[]
trace.append(min(Ob))
#######################差分进化主体循环#############
for gen in range(0,G):
#############变异操作##########
for m in range(0,NP):
r1=random.randint(0,NP-1)
while (r1==m):
r1 = random.randint(0, NP-1)
r2 = random.randint(0, NP-1)
while (r2 == m) and (r2 == r1):
r2 = random.randint(0, NP-1)
r3 = random.randint(0, NP-1)
while (r3 == m) and (r3==r2) and (r3==r1):
r3 = random.randint(0,NP-1)
#差分
v[:,m]=x[:,r1]+F*(x[:,r2]-x[:,r3])
###########交叉操作############
r=np.random.randint(1,3)
for n in range(0,D):
cr=np.random.rand()
if (cr<=CR and n==r):
u[n,:]=v[n,:]
else:
u[n,:]=x[n,:]
#######################边界条件处理##############
for n in range(0,D):
for m in range(0,NP):
if u[n,m]<Xx:
u[n,m]=Xx
if u[n,m]>Xs:
u[n,m]=Xs
###########选择操作##########
Ob1=[]
for m in range(0,NP):
Ob1.append(fun_fitness(u[:,m]))
for m in range(0,NP):
if Ob1[m]<Ob[m]:
x[:,m]=u[:,m]
for m in range(0,NP):
Ob[m]=fun_fitness(x[:,m])
trace.append(min(Ob))
sort_Ob=sorted(Ob)
sort_Ob_index=get_index(Ob,sort_Ob)
x=x[:,sort_Ob_index]
X=x[:,0]
Y=min(Ob)
#################画图##############
fig2=plt.figure(2)
ax2 = plt.axes()
plt.plot(trace)
plt.grid()
plt.xlabel("Iterations")
plt.ylabel("objective function")
plt.title("DE objective function curve")
plt.show()
求解结果:
**注意:**差分进化算法具有随机性,每次求解结果并不完全相同,展示的仅为作者本人本次运行的截图。