差分进化算法_差分进化算法入门及实例应用

这个就当是自己学习的整理总结吧。

一、差分进化算法理论

差分进化算法是基于群体智能理论的优化算法，通过群体内个体间的合作与竞争而产生的智能优化搜索算法。

1.1 差分进化特点

（1）结构简单，容易使用。主要通过差分变异算子来进行遗传操作。

（2）性能优越。具有较好的可靠性，鲁棒性，高效性。

（3）自适应性。差分变异算子可以是一个常数，也可以是一个具有变异步长和搜索方向的自适应能力。

（4）具有内在并行性，可协同搜索。在同一要求下，差分进化算法拥有更快的收敛速度。

1.2 差分进化算法操作

（1）初始化

差分进化算法利用NP个维数为D的实数值参数向量，作为每一代的种群。一般初始化种群符合均匀概率分布。

（2）变异

基本变异向量由下式产生

要求r1,r2,r3，i互不相同，所以NP必须大于4，F是变异算子通常在0-2之间，它控制偏差向量的放大作用。

通俗的理解就是把两个向量的差值乘上变异算子加给第三个向量作为新的变异向量。

（3）交叉

为了增加干扰向量的多样性，引入交叉操作。

randb(j)表示产生[0,1]之间随机数第j个估计值，rnbr表示一个随机选择的序列。CR是交叉算子。

通俗理解就是如果随机产生的randb(j)小于CR或者j=r,那么就将变异后的种群放入选择群体中，如果不是就就将原来的种群放入选择种群中。

（4）选择

为了决定选择种群中的向量是否能成为下一代的一员，试验向量与当前的目标向量进行比较，如果出现目标函数被最小化，那么具有最小目标函数的向量将在下一代出现。试验向量只与一个个体进行比较，而不是所有个体。

（5）边界条件处理

一种方法是将超过边界的向量使用可行域中随机产生的参数向量代替。

例如：u(n,m)=rand*(Xs-Xx)+Xx;

另一种方法是，进行边界吸收处理，直接放栗子好理解

 if u(n,m) 
    
      Xs u(n,m)=Xs; end

二、改进的差分进化算法

（1）自适应差分进化算法

主要是使用了自适应算子。在基本的差分进化算法中，变异算子经常取常数，比较难准确确定，变异率太大，全局最优解低，变异率小，群体多样性下降，易出现‘早熟’的现象。我们可以设计这样的变异算子：

这样开始的时候变异算子为2F0，在初期可以保持多样性，防止早熟。随着进展，变异算子降低最后变为F0，避免最优解遭到破坏。

还可以设计一个随机范围的交叉算子：

这样交叉算子的均值就在0.75，保持了群体多样性。

（2）离散差分进化算法

采用的是浮点数编码，用到了floor这个函数向下取整即可。

三、差分进化算法流程

四、参数说明

种群数量NP：规模越大多样性越好，但是增加计算难度，一般在5D~10D之间，而且必须大于4。

变异算子F：

,决定偏差向量的放大比例。F=0.5是一个比较好的初始选择，如果过早收敛可以增大F或者NP

交叉算子CR：

,CR越大，发生交叉的可能性就越大，一般CR=0.1，较大的CR会加速收敛。

终止条件：一般当目标函数值小于阈值时程序终止，一般取

。

五、实例练习

（1）计算函数

的最小值，其中个体x的维数n=10，这是一个简单的平方和函数，只有一个极小点x=（0,0,…,0），理论上最小值分f(0,0,…,0)=0。

解：使用的是自适应进化算法，所以变异算子和变异算子不是常数

clear all; close all; clc; NP=50; %种群数量 D=10; %变量的维数 G=200; %最大进化代数 F0=0.4; %初始变异算子 CR=0.1; %交叉算子 Xs=20; %上限 Xx=-20; %下限 yz=10^-6; %阈值 %赋初值 x=zeros(D,NP); %初始种群 v=zeros(D,NP); %变异种群 u=zeros(D,NP); %选择种群 x=rand(D,NP)*(Xs-Xx)+Xx; %赋初值 %计算目标函数 for m=1:NP Ob(m)=func1(x(:,m)); end trace(1)=min(Ob); %差分进化循环 for gen=1:G %变异操作 %自适应变异算子 lamda=exp(1-G/(G+1-gen)); F=F0*2^lamda; %r1,r2,r3,m互不相同 for m=1:NP r1=randi([1,NP],1,1); while (r1==m) r1=randi([1,NP],1,1); end r2=randi([1,NP],1,1); while (r2==m)||(r1==r2) r2=randi([1,NP],1,1); end r3=randi([1,NP],1,1); while (r3==m)||(r2==r3)||(r1==r3) r3=randi([1,NP],1,1); end v(:,m)=x(:,r1)+F*(x(:,r2)-x(:,r3)); end %交叉操作 r=randi([1,D],1,1); for n=1:D cr=rand(1); if (cr<=CR)||(n==r) u(n,:)=v(n,:); else u(n,:)=x(n,:); end end %边界条件处理 for n=1:D for m=1:NP if (u(n,m) 
    
      Xs) u(n,m)=rand*(Xs-Xx)+Xx; end end end %选择操作 for m=1:NP Ob1(m)=func1(u(:,m)); end for m=1:NP if Ob1(m)

提示一下size 函数的用法，size(x,1)返回行数，size(x,2)返回列数。

（2）求函数

的最小值，其中x的范围是[-4,4]，y的取值范围是[-4,4]，这是一个有多个局部极值的函数，函数图形如图所示。

解：

clear all; close all; clc; NP=20; %种群数量 D=2; %变量的维数 G=100; %最大进化代数 F=0.5; %变异算子 CR=0.1; %交叉算子 Xs=4; %上限 Xx=-4; %下限 %赋初值 x=zeros(D,NP); %初始种群 v=zeros(D,NP); %变异种群 u=zeros(D,NP); %选择种群 x=rand(D,NP)*(Xs-Xx)+Xx; %赋初值 %计算目标函数 for m=1:NP Ob(m)=func2(x(:,m)); end trace(1)=min(Ob); %差分进化循环 for gen=1:G %变异操作 %r1,r2,r3,m各不相同 for m=1:NP r1=randi([1,NP],1,1); while (r1==m) r1=randi([1,NP],1,1); end r2=randi([1,NP],1,1); while (r1==r2)||(r2==m) r2=randi([1,NP],1,1); end r3=randi([1,NP],1,1); while (r1==r3)||(r2==r3)||(r3==m) r3=randi([1,NP],1,1); end v(:,m)=x(:,r1)+F*(x(:,r2)-x(:,r3)); end %交叉操作 r=randi([1,D],1,1); for n=1:D cr=rand(1); if(cr 
    
      Xs u(n,m)=Xs; end end end %选择操作 for m=1:NP Ob1(m)=func2(u(:,m)); end for m=1:NP if Ob1(m)

结果：当x=-4,y=-3.9478时取到最小值-10.9734

（3）用离散差分算法求函数

的最大值，其中x取-100~100之间的整数，y取-100~100之间的整数，其函数图像如图所示:

解：

clear all; close all; clc; NP=20; %种群数量 D=2; %变量的维数 G=100; %最大进化代数 F=0.5; %变异算子 CR=0.1; %交叉算子 Xs=100; %上限 Xx=-100; %下限 %赋初值 x=zeros(D,NP); %初始种群 v=zeros(D,NP); %变异种群 u=zeros(D,NP); %选择种群 x=randi([Xx,Xs],D,NP); %赋初值 %计算目标函数 for m=1:NP Ob(m)=func3(x(:,m)); end trace(1)=max(Ob); %差分进化循环 for gen=1:G for m=1:NP r1=randi([1,NP],1,1); while (r1==m) r1=randi([1,NP],1,1); end r2=randi([1,NP],1,1); while (r1==r2)||(r2==m) r2=randi([1,NP],1,1); end r3=randi([1,NP],1,1); while (r3==m)||(r3==r2)||(r3==r1) r3=randi([1,NP],1,1); end v(:,m)=floor(x(:,r1)+F*(x(:,r2)-x(:,r3))); %floor向下取整 end %交叉操作 r=randi([1,D],1,1); for n=1:D cr=rand(1); if (cr<=CR)||(n==r) u(n,:)=v(n,:); else u(n,:)=x(n,:); end end %边界吸收 for n=1:D for m=1:NP if u(n,m) 
    
      Xs u(n,m)=Xs; end end end %选择操作 for m=1:NP Ob1(m)=func3(u(:,m)); end for m=1:NP if Ob1(m)>Ob(m) x(:,m)=u(:,m); end end for m=1:NP Ob(m)=func3(x(:,m)); end trace(gen+1)=max(Ob); end [SortOb,Index]=sort(Ob); X=x(:,Index); XBest=x(:,end); %最优变量 Y=max(Ob); %最优 %画图 figure plot(trace) xlabel('迭代次数') ylabel('目标函数值') title('DE目标函数曲线') %%适应度函数 function y=func3(x) y=-((x(1).^2+x(2)-1).^2+(x(1)+x(2).^2-7).^2)/200+10; end