matlab新手快速上手3（差分进化算法）

本文用经典差分进化框架模板，对matlab新手友好，快速上手看懂matlab代码，快速应用实践，源代码在文末给出。

差分进化算法定义：

差分进化算法（Differential Evolution，简称DE算法）是一种优化算法，由Storn和Price于1997年提出。这个算法主要用于解决连续优化问题。

DE算法的基本思想是模拟生物进化中的个体竞争和适者生存的原理。它通过维护一个种群，每个个体都代表问题的一个解，并通过不断的迭代和交叉变异来寻找最优解。

DE算法的核心操作包括选择、交叉和变异：

选择：通过比较个体的适应度来选择优秀的个体进入下一代。
变异：随机选择种群中的个体，并基于差分向量进行变异操作，产生新的个体。
交叉：将变异后的个体与原始个体进行交叉操作，产生新的解。

DE算法具有较好的全局搜索能力和较快的收敛速度，在许多优化问题中表现良好。它已经被广泛应用于工程优化、信号处理、机器学习等领域。

开始编程：

参数与子函数定义：

function DE
%---------------------------------- 参数设置 -------------------------------
NP=100;                            %种群规模
D=10;                              %变量个数
MinX=-30;                          %范围下限
MaxX=30;                           %范围上限
alpha=0.6;                         %缩放因子
beta=0.6;                          %缩放因子
CR=0.8;                            %交叉概率
Error=0.001;                         %限定精度
Max_N=1000;                        %限定代数
flagc=[0,Max_N];                   %收敛标志
%---------------------- 粒子位置初始化 ----------------------
for i=1:1:DX(:,i)=MinX+(MaxX-MinX)*rand(NP,1);selfX(:,i)=X(:,i);
end
%----------------------- 计算函数值 -------------------------
F=fun(X);
selfF=F;
%---------------------------- 子函数1：目标函数 ----------------------------
function F=fun(X)
for i=1:1:size(X,1)for j=1:1:size(X,2)x(j)=X(i,j);endfor j=1:1:size(X,2)-1temp(j)=100*(x(j+1)-x(j)^2)^2+(x(j)-1)^2;endF(i)=sum(temp);
end
%--------------------- 求最优解 ---------------------------
[Bestf,Indexf]=sort(F,2);        %对NP个函数值排序Bestfi=Bestf(1);                 %第1个函数值最小
Bestp=Indexf(1);                 %最优粒子序号

粒子位置初始化：

此函数生成一个X矩阵，为NP行，D列的矩阵，矩阵元素为0~1的随机值，每一行代表一个个体，列的元素代表每个自变量，selfX矩阵与它相等。

子函数定义：

根据fun(X)得知，子函数为经典的罗森布鲁克函数： $f(x) = \sum_{i=1}^{D-1} \left[ 100 \times (x_{i+1} - x_i^2)^2 + (x_i - 1)^2 \right]$ ，此函数在x = 1处时，函数最小为0。差分进化算法就是寻找此函数找到最小值的点，最终输出的X中最优个体中的元素应该都接近于1。

附：罗森布鲁克函数（百度百科）

程序主体：

%--------------------- 程序主循环开始 ----------------------
for gen=1:1:Max_Ntime(gen)=gen;%----------------- 变异操作 ---------------------------for i=1:1:NPflag1=ceil(rand*NP);while(flag1==i)flag1=ceil(rand*NP);endflag2=ceil(rand*NP);while((flag2==i)|(flag2==flag1))flag2=ceil(rand*NP);endX(i,:)=X(i,:)+alpha*(X(Bestp,:)-X(i,:))+beta*(X(flag1,:)-X(flag2,:));end%----------------- 交叉操作 ---------------------------for i=1:1:NPtemp=rand(1,D);for j=1:1:Dif temp(j)>CRX(i,j)=selfX(i,j);end                endend%----------------- 计算函数值 ---------------------------F=fun(X);%----------------- 选择操作 ---------------------------for i=1:1:NPif F(i)>=selfF(i)F(i)=selfF(i);X(i,:)=selfX(i,:);            endend%----------------- 迭代更新 ---------------------------for i=1:1:NPselfF(i)=F(i); selfX(i,:)=X(i,:);            end%----------------- 求最优解 ---------------------------[Bestf,Indexf]=sort(F,2);        %对NP个函数值排序Bestfi=Bestf(1);                 %第1个函数值最小Bestp=Indexf(1);                 %最优粒子序号%----------------------------- 记录结果 --------------------------------result(gen)=Bestfi;if mod(gen,10)==0disp(sprintf('代数：%d -------- 结果：%f',gen,Bestfi));plot(time,result,'r');axis([1,Max_N,0,100]);xlabel('迭代步数');ylabel('优化结果');drawnow;pause(0.1);endif Bestfi<Error break;end
end
disp(' ');
disp(sprintf('迭代步数：%d -------- 优化结果：%f',gen,Bestfi));
disp(X);

最外层循环：

大循环表示迭代Max_N次，time存储每次循环的迭代次数。

变异操作：

for循环遍历NP次，表示对每个个体都进行变异，flag1表示一个[1，NP]之间的随机整数，while控制这个整数不能与当前个体序号相等，flag1与flag2也不能相同，也就是防止选择相同的两个个体导致变异不明显。后续的X(i,:)的操作就是让当前个体加上0.6倍的与种群最优个体的差值向量（这个就是差分变异），再加上0.6倍的两个随机个体的差值向量。这样就实现了一个个体的变异操作。

交叉操作：

交叉操作很明显，就是根据交叉概率，将当前的个体的某些参数变为没变异之前的参数

选择操作：

选择操作就是根据适应度进行选择，将适应度低的个体变为上一代的个体。

迭代更新：

将这一代的个体更新到self中，为下一次迭代做准备，也就是在下一代更新后，将下一代与这个self做对比。

后续就是排序，绘图，输出，等操作，就是将个体按照适应度排序，将最优个体进行输出

源代码：

%---------------------------------- 程序说明 -------------------------------%                          该程序实现了基本差分变异算法%---------------------------------- 程序正文 -------------------------------
function DE
%---------------------------------- 参数设置 -------------------------------
NP=100;                            %种群规模
D=10;                              %变量个数
MinX=-30;                          %范围下限
MaxX=30;                           %范围上限
alpha=0.6;                         %缩放因子
beta=0.6;                          %缩放因子
CR=0.8;                            %交叉概率
Error=0.001;                         %限定精度
Max_N=1000;                        %限定代数
flagc=[0,Max_N];                   %收敛标志
%---------------------- 粒子位置初始化 ----------------------
for i=1:1:DX(:,i)=MinX+(MaxX-MinX)*rand(NP,1);selfX(:,i)=X(:,i);
end
%X=MinX+(MaxX-MinX)*rand(NP,D);
%selfX = X;
%----------------------- 计算函数值 -------------------------
F=fun(X);
selfF=F;
%--------------------- 求最优解 ---------------------------
[Bestf,Indexf]=sort(F,2);        %对NP个函数值排序Bestfi=Bestf(1);                 %第1个函数值最小
Bestp=Indexf(1);                 %最优粒子序号
%--------------------- 程序主循环开始 ----------------------
for gen=1:1:Max_Ntime(gen)=gen;%----------------- 变异操作 ---------------------------for i=1:1:NPflag1=ceil(rand*NP);while(flag1==i)flag1=ceil(rand*NP);endflag2=ceil(rand*NP);while((flag2==i)|(flag2==flag1))flag2=ceil(rand*NP);endX(i,:)=X(i,:)+alpha*(X(Bestp,:)-X(i,:))+beta*(X(flag1,:)-X(flag2,:));end%----------------- 交叉操作 ---------------------------for i=1:1:NPtemp=rand(1,D);for j=1:1:Dif temp(j)>CRX(i,j)=selfX(i,j);end                endend%----------------- 计算函数值 ---------------------------F=fun(X);%----------------- 选择操作 ---------------------------for i=1:1:NPif F(i)>=selfF(i)F(i)=selfF(i);X(i,:)=selfX(i,:);            endend%----------------- 迭代更新 ---------------------------for i=1:1:NPselfF(i)=F(i); selfX(i,:)=X(i,:);            end%----------------- 求最优解 ---------------------------[Bestf,Indexf]=sort(F,2);        %对NP个函数值排序Bestfi=Bestf(1);                 %第1个函数值最小Bestp=Indexf(1);                 %最优粒子序号%----------------------------- 记录结果 --------------------------------result(gen)=Bestfi;if mod(gen,10)==0disp(sprintf('代数：%d -------- 结果：%f',gen,Bestfi));plot(time,result,'r');axis([1,Max_N,0,100]);xlabel('迭代步数');ylabel('优化结果');drawnow;pause(0.1);endif Bestfi<Error break;end
end
disp(' ');
disp(sprintf('迭代步数：%d -------- 优化结果：%f',gen,Bestfi));
disp(X);
%---------------------------- 子函数1：目标函数 ----------------------------
function F=fun(X)
for i=1:1:size(X,1)for j=1:1:size(X,2)x(j)=X(i,j);endfor j=1:1:size(X,2)-1temp(j)=100*(x(j+1)-x(j)^2)^2+(x(j)-1)^2;endF(i)=sum(temp);
end