智能优化算法之差分进化算法

差分进化算法原理

差分进化算法是一种随机的启发式搜索算法，简单易用，有较强的鲁棒性和全局搜索能力。
差分进化算法是一种自组织最小化方法，利用种群中随机选择的不同向量来干扰一个现有向量，种群中的每个向量均要受到干扰，其中种群干扰向量可独立进行的，由此说明，其进行是并行的。
差分进化算法一共有五个步骤：

初始化
变异
交叉
选择
边界条件处理

通过上述流程可获得一个收敛性非常好的结果，引导搜索结果向全局最优解逼近。

差分进化算法的流程详解：

初始化
差分进化算法利用NP个维数为D的实数值参数向量，将它们作为每一代的种群，个体表示为： xi,G(i=1,2,3,⋯,NP)x_{i,G}(i=1,2,3,\cdots,NP)xi,G(i=1,2,3,⋯,NP)
NP为种群规模，iii表示在种群中的序列，G表示进化代数，在此过程中种群规模始终保持不变。
通常初始化种群是在给定约束条件内随机生成，一般假定所有随机生成种群符合均匀概率分布。如果可以预先得到问题的初步解，则可通过对初步解加入正态分布随机偏差来产生，可提高重建效率。
初始化式子为
xji,0=rand[0,1]×(xjU−xjL)+xjLx_{ji,0}=rand[0,1]\times (x_{j}^{U}-x_{j}^{L})+x_{j}^{L}xji,0=rand[0,1]×(xjU−xjL)+xjL
xjLx_{j}^{L}xjL为下界限，xjUx_{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间的随机数

变异
在基本差分进化算法中目标向量xi,Gx_{i,G}xi,G由下式产生：
Vi,G+1=xr1,G+F(xr2,G−xr3,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)
式子中r1,r2,r3,ir1,r2,r3,ir1,r2,r3,i 为 [1,NP][1,NP][1,NP]上四个不同的序号。
变异算子FFF为[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])

交叉
交叉为了增加干扰参数向量Ui,G+1U_{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即为交叉算子。

选择
为决定实验向量是否会成为下一代中的成员，按照贪婪策略将实验向量Ui,G+1U_{i,G+1}Ui,G+1与目标向量xi,Gx_{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]=Xxif u[n,m]>Xs:u[n,m]=Xs

实例

求函数f(x,y)=3cos(xy)+x+yf(x,y)=3cos(xy)+x+yf(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]=Xxif 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()

求解结果：
**注意：**差分进化算法具有随机性，每次求解结果并不完全相同，展示的仅为作者本人本次运行的截图。