差分进化算法求解函数最值问题

声明：
1.本文源代码来自书目《智能优化算法及其MATLAB实例（第3版）》，目的在于为MATLAB初学者提供更简明的代码解析，方便读者了解算法及MATLAB编程基本原理。 2.文中代码每一行后都有相应注释，因此本文是一篇适合所有代码水平的学习者阅读的文章。如果觉得有帮助，麻烦点个赞哦！

下面3.1、3.2、3.3三个例子都是函数求最值问题，且均采用差分进化算法编程实现，基本代码框架相同，因此集中于一篇文章，便于读者对照阅读学习

例3.1

计算函数f(x)=x12+x22+···xn2(-20<=xi<=20)的最小值，其中个体x的维数n=10。这是一个简单的平方和函数，只有一个极小点在0处取得为0。
（在此之前例2.2中，同样的这样一个问题我们使用了遗传算法求解
戳它：https://blog.csdn.net/dhdjenrjj/article/details/115534284?spm=1001.2014.3001.5501）

基于基本差分进化算法的解题代码如下：

%%%%%差分进化算法求函数极值%%%%%%%%
%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%
clear all;%清除所有变量
close all;%清图
clc;%清屏
NP=50;%种群数量
D=10;%每个向量的维度
G=200;%最大进化代数
F0=0.4;%初始变异算子，由于此问题中我们采用自适应变异法（F值会发生变化）
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:NPOb(m)=func1(x(:,m));%调用代码段末尾的目标函数式进行计算
end
trace(1)=min(Ob);
%%%%%%%%%%差分进化循环%%%%%%%%%%%%
for gen=1:G%%%%%%%变异操作%%%%%%%%%%%%%%%%%%%%自适应变异算子%%%%%%%%%%%lamda=exp(1-G/(G+1-gen));F=F0*2^(lamda);%变异采取随机方式x=rand,差向量个数y=1%for m=1:NP%依次将序号值为1-NP的向量作为目标向量，分别生成相应变异向量r1=randi([1,NP],1,1);%从1-NP中随机选择一个序号作为r1while(r1==m)r1=randi([1,NP],1,1)%保证r1不等于目标向量序号值mendr2=randi([1,NP],1,1);%从1-NP中随机选择一个序号作为r1while(r2==m)|(r2==r1)r2=randi([1,NP],1,1)%保证r2不等于目标向量序号值m与r1endr3=randi([1,NP],1,1);%从1-NP中随机选择一个序号作为r1while(r3==m)|(r3==r1)|(r3==r2)r3=randi([1,NP],1,1)%保证m,r1,r2,r3互不相等endv(:,m)=x(:,r1)+F*(x(:,r2)-x(:,r3));%生成第gen代种群中序号为m的向量对应的变异向量end%%%%%%%%交叉操作%%%%%%%%%%%%%%r=randi([1,D],1,1);%从1-D生成随机整数for n=1:Dcr=rand(1);%生成【0,1】随机实数if(cr<=CR)|(n==r)%限定这个条件，能够保证至少有一位发生变异u(n,:)=v(n,:);%满足条件则将原种群每个向量相应位替换为变异数组中相应位置的数elseu(n,:)=x(n,:);endend%%%%%%%%%%边界条件的处理%%%%%%%%%%for n=1:Dfor m=1:NPif (u(n,m)<Xx)|(u(n,m)>Xs)%检查种群中每个向量每个维度的值u(n,m)=rand*(Xs-Xx)+Xx;%如果超出范围，则随机生成一个合法值替换越界的值endendend %%%%%%%%%选择操作（生成下一代）%%%%%%for m=1:NPOb1(m)=func1(u(:,m));%对当前一代的种群中每个向量计算目标函数值进行评价endfor m=1:NPif Ob1(m)<Ob(m)x(:,m)=u(:,m);%在上一代x(gen-1)与本代选择种群u(gen)中各向量中选择目标函数值较优的，作为本代生成的种群x(gen)endendfor m=1:NPOb(m)=func1(x(:,m));%计算新种群目标函数值end[minOb,index]=min(Ob);trace(gen+1)=min(Ob);%记录最优个体进化曲线if min(Ob(m))<yz%若最优个体目标函数值小于阈值，则退出差分进化循环breakend
end
disp(['最终最优个体为'])
x(:,index)
disp(['函数极小值为'])
minOb
%%%%%%%%%%%画图%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(trace);
xlabel('迭代次数')
ylabel('目标函数值')
title('DE目标函数曲线')
%%%%%%%%%%适应度函数%%%%%%%%%%%%%%%%%%%
function result=func1(x)
summ=sum(x.^2);
result=summ;
end

最终，程序运行结果如下：

例3.2

求函数f(x,y)=3cos(xy)+x+y的最小值，其中x的取值范围为[-4,4]。这是一个有多个局部极值的函数。

基于基本差分进化算法的MATLAB程序如下：

%%%%%%%%%%%差分进化算法求二元函数最值%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;%清除所有变量
close all;%清图
clc;%清屏
NP=20;%种群数量
D=2;%变量的维数为2，因为目标函数自变量x,y是两个
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:NPOb(m)=func2(x(:,m));
end
trace(1)=min(Ob);
%%%%%%%%差分进化循环%%%%%%%%%%%%%%%%%%%%%%%%%
for gen=1:G%%%%%%%%%%%变异操作%%%%%%%%%%%%%%%%%%%%%%变异采取随机方式x=rand,差向量个数y=1%for m=1:NP%依次将序号值为1-NP的向量作为目标向量，分别生成相应变异向量r1=randi([1,NP],1,1);%从1-NP中随机选择一个序号作为r1while(r1==m)r1=randi([1,NP],1,1)%保证r1不等于目标向量序号值mendr2=randi([1,NP],1,1);%从1-NP中随机选择一个序号作为r1while(r2==m)|(r2==r1)r2=randi([1,NP],1,1)%保证r2不等于目标向量序号值m与r1endr3=randi([1,NP],1,1);%从1-NP中随机选择一个序号作为r1while(r3==m)|(r3==r1)|(r3==r2)r3=randi([1,NP],1,1)%保证m,r1,r2,r3互不相等endv(:,m)=x(:,r1)+F*(x(:,r2)-x(:,r3));%生成第gen代种群中序号为m的向量对应的变异向量end%%%%%%%%交叉操作%%%%%%%%%%%%%%r=randi([1,D],1,1);%从1-D生成随机整数for n=1:Dcr=rand(1);%生成【0,1】随机实数if(cr<=CR)|(n==r)%限定这个条件，能够保证至少有一位发生变异u(n,:)=v(n,:);%满足条件则将原种群每个向量相应位替换为变异数组中相应位置的数elseu(n,:)=x(n,:);endend%%%%%%利用边界吸收进行边界条件的处理%%%%%%%%%%for n=1:Dfor m=1:NPif u(n,m)<Xxu(n,m)=Xx;endif u(n,m)>Xsu(n,m)=Xs;%若超出所给范围，则使得该值改为与之最相近的边界值endendend%%%%%%%%%选择操作%%%%%%%%%%%%%%%%%%%%%for m=1:NPOb1(m)=func2(u(:,m));%计算选择种群的各向量适应度endfor m=1:NPif Ob1(m)<Ob(m)x(:,m)=u(:,m);%选择两矩阵中编号为m的位置处，目标函数值更优的作为下一代endendfor m=1:NPOb(m)=func2(x(:,m));%更新下一代适应度值endtrace(gen+1)=min(Ob);
end
[SortOb,Index]=sort(Ob);%按目标函数值升序排列
best=x(:,Index);
X=best(:,1)%最小值点
Y=min(Ob)%最小值
%%%%%%%%%%%%%画图%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(trace)
xlabel('迭代次数')
ylabel('目标函数值')
title('DE目标函数曲线')
%%%%%%%%%%适应度函数%%%%%%%%%%%%%%%%%%%%%%%%%
function value=func2(x)
value=3*cos(x(1)*x(2))+x(1)+x(2);%该函数的输入参数为x是一个二维向量，x(1)即x,x(2)即y
end

最终程序运行效果如下：

例3.3

用离散差分进化算法求函数f(x,y)=-((x²+y-1)²+(x+y²-7)²)/200+10的最大值，其中x和y都是取值为-100到100之间的整数。

MATLAB程序如下：

%%%%%%%%%%%离散差分进化算法求二元函数最值%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;%清除所有变量
close all;%清图
clc;%清屏
NP=20;%种群数量
D=2;%变量的维数为2，因为目标函数自变量x,y是两个
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:NPOb(m)=func3(x(:,m));
end
trace(1)=max(Ob);
%%%%%%%%差分进化循环%%%%%%%%%%%%%%%%%%%%%%%%%
for gen=1:G%%%%%%%%%%%变异操作%%%%%%%%%%%%%%%%%%%%%%变异采取随机方式x=rand,差向量个数y=1%for m=1:NP%依次将序号值为1-NP的向量作为目标向量，分别生成相应变异向量r1=randi([1,NP],1,1);%从1-NP中随机选择一个序号作为r1while(r1==m)r1=randi([1,NP],1,1)%保证r1不等于目标向量序号值mendr2=randi([1,NP],1,1);%从1-NP中随机选择一个序号作为r1while(r2==m)|(r2==r1)r2=randi([1,NP],1,1)%保证r2不等于目标向量序号值m与r1endr3=randi([1,NP],1,1);%从1-NP中随机选择一个序号作为r1while(r3==m)|(r3==r1)|(r3==r2)r3=randi([1,NP],1,1)%保证m,r1,r2,r3互不相等endv(:,m)=floor(x(:,r1)+F*(x(:,r2)-x(:,r3)));%离散DE算法要向下取整生成第gen代种群中序号为m的向量对应的变异向量end%%%%%%%%交叉操作%%%%%%%%%%%%%%r=randi([1,D],1,1);%从1-D生成随机整数for n=1:Dcr=rand(1);%生成【0,1】随机实数if(cr<=CR)|(n==r)%限定这个条件，能够保证至少有一位发生变异u(n,:)=v(n,:);%满足条件则将原种群每个向量相应位替换为变异数组中相应位置的数elseu(n,:)=x(n,:);endend%%%%%%利用边界吸收进行边界条件的处理%%%%%%%%%%for n=1:Dfor m=1:NPif u(n,m)<Xxu(n,m)=Xx;endif u(n,m)>Xsu(n,m)=Xs;%若超出所给范围，则使得该值改为与之最相近的边界值endendend%%%%%%%%%选择操作%%%%%%%%%%%%%%%%%%%%%for m=1:NPOb1(m)=func3(u(:,m));%计算选择种群的各向量适应度endfor m=1:NPif Ob1(m)>Ob(m)x(:,m)=u(:,m);%选择两矩阵中编号为m的位置处，目标函数值更优的作为下一代endendfor m=1:NPOb(m)=func3(x(:,m));%更新下一代适应度值endtrace(gen+1)=max(Ob);
end
[SortOb,Index]=sort(Ob);%按目标函数值升序排列,Index是一个按照适应度升序排列的序号数组
sort_x=x(:,Index);
X=sort_x(:,end)%最大值点
Y=max(Ob)%最大值
%%%%%%%%%%%%%画图%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(trace);
xlabel('迭代次数')
ylabel('目标函数值')
title('DE目标函数曲线')
%%%%%%%%%%适应度函数%%%%%%%%%%%%%%%%%%%%%%%%%
function value=func3(x)
value=-((x(1).^2+x(2)-1).^2+(x(1)+x(2).^2-7).^2)/200+10;%该函数的输入参数为x是一个二维向量，x(1)即x,x(2)即y
end

最终程序运行结果如下：