【Matlab】智能优化算法_平衡优化器算法EO
- 1.背景介绍
- 2.数学模型
- 2.1 初始化和功能评估
- 2.2 平衡池和候选者（Ceq）
- 2.3 指数项（F）
- 2.3 生成率（G）
- 3.文件结构
- 4.伪代码
- 5.详细代码及注释
- 5.1 EO.m
- 5.2 Get_Functions_details.m
- 5.3 initialization.m
- 5.4 main.m
- 6.运行结果
- 7.参考文献
1.背景介绍
EO方法的灵感来源于控制体积上的简单混合动态质量平衡，其中使用质量平衡方程来描述控制体积中非反应成分的浓度，作为其各种源和汇机制的函数。质量平衡方程为进入、离开和在控制体积中产生的质量守恒提供了基础物理。表示一般质量平衡方程的一阶常微分方程被描述为：
C是控制体积内的浓度（V），V（dC/dt）是控制体积内的质量变化率。是控制体积内的质量变化率，Q是进出控制体积的体积流量，Ceq代表平衡状态下的浓度，其中控制体积内没有生成，G是控制体积内的质量生成率。当V（dC/dt） 达到零时，就达到了一个稳定的平衡状态。通过对公式（1）的重新排列，可以求出dC/dt与Q/V的函数关系；其中Q/V代表停留时间的倒数，这里称为λ，或周转率（即λ=Q/V）。随后，公式（1）也可以重新排列，以解决控制体积（C）中的浓度与时间（t）的关系：
方程（3）显示了方程（2）随时间的积分：
这导致：
在等式（4）中，F的计算如下：
其中t0和C0是初始开始时间和浓度，取决于积分间隔。等式（4）可用于估计具有已知周转率的控制体积中的浓度，或者用于使用具有已知生成率和其他条件的简单线性回归来计算平均周转率。
在EO中，粒子类似于溶液，浓度类似于粒子在PSO算法中的位置。如等式（4）所示，有三个项表示粒子的更新规则，每个粒子通过三个单独的项更新其浓度。第一项是平衡浓度，定义为从一个称为平衡池的池中随机选择的迄今为止最好的溶液之一。第二项与粒子和平衡状态之间的浓度差有关，这起到了直接搜索机制的作用。这个术语鼓励粒子全局搜索域，充当探索者。第三个术语与生成率有关，生成率主要扮演开发者或解决方案细化者的角色，尤其是小步骤，尽管有时也会扮演探索者的角色。每个术语及其影响。
2.数学模型
2.1 初始化和功能评估
与大多数元启发式算法类似，EO使用初始种群来启动优化过程。初始浓度是基于搜索空间中具有均匀随机初始化的粒子数量和维度构建的，如下所示：
C initial i是第i个粒子的初始浓度矢量，Cmin和Cmax表示维度的最小值和最大值，Rand i是[0，1]区间内的随机矢量，n是作为总体的粒子数。评估粒子的适应度函数。
2.2 平衡池和候选者（Ceq）
平衡状态是算法的最终收敛状态，它被期望为全局最优。在优化过程开始时，没有关于平衡状态的知识，并且仅确定平衡候选者以提供粒子的搜索模式。基于不同类型案例问题下的不同实验，这些候选粒子是在整个优化过程中识别出的迄今为止最好的四个粒子加上另一个粒子，其浓度是上述四个粒子的算术平均值。这四个候选者有助于EO具有更好的勘探能力，而平均值有助于开发。候选的数量是任意的，并且基于优化问题的类型。可以使用其他数量的候选者（例如3或5）。例如，GWO使用三种迄今为止最好的候选者（阿尔法狼、贝塔狼和伽马狼）来更新其他狼的位置。
每次迭代中的每个粒子都通过在以相同概率选择的候选者中随机选择来更新其浓度。例如，在第一次迭代中，第一个粒子基于Ceq（1）更新其所有浓度；然后，在第二次迭代中，它可以基于Ceq（ave）更新其浓度。在优化过程结束之前，每个粒子都将经历更新过程，所有候选解决方案都将为每个粒子接收大致相同数量的更新。
2.3 指数项（F）
对主要浓度更新规则有贡献的下一项是指数项（F）。该术语的准确定义将有助于EO在勘探和开采之间取得合理的平衡。由于实际控制体积中的周转率可以随时间变化，因此假设λ是[0，1]区间内的随机向量。
其中，时间t被定义为迭代（Iter）的函数，因此随着迭代次数的增加而减少：
其中，Iter和Max_Iter分别表示当前迭代次数和最大迭代次数，a2是用于管理利用能力的常数值。为了通过降低搜索速度来保证收敛性，同时提高算法的探索和利用能力，本研究还考虑：
其中a1是控制勘探能力的常数值。a1越高，勘探能力越好，因此开采性能越低。同样，a2越高，开采能力越好，勘探能力越低。第三个组成部分，符号（r−0.5），对勘探和开发方向的影响。r是介于0和1之间的随机向量。对于本文随后解决的所有问题，a1和a2分别等于2和1。这些常数是通过测试函数子集的经验测试来选择的。然而，这些参数可以根据需要针对其他问题进行调整。
等式（11）显示了等式（8）的修订版本，将等式（10）替换为等式（（8）：
2.3 生成率（G）
生成率是所提出的算法中最重要的术语之一，通过改进开发阶段来提供精确的解决方案。例如，将生成率描述为一阶指数衰减过程的一个多用途模型定义为：
其中G0是初始值，k表示衰减常数。为了具有更可控和系统的搜索模式，并限制随机变量的数量，本研究假设k=λ，并使用先前导出的指数项。因此，生成速率方程的最终集合如下：
其中r1和r2是[0，1]中的随机数，并且GCP向量是通过重复由等式（15）得到的相同值来构造的。在这个方程中，GCP被定义为生成率控制参数，它包括生成项对更新过程的贡献的可能性。指定有多少粒子使用生成项来更新其状态的这种贡献的概率由另一个称为生成概率（GP）的项决定。这种贡献的机制由等式（14）和（15）决定。方程（15）发生在每个粒子的水平上。例如，如果GCP为零，则G等于零，并且在没有生成率项的情况下更新该特定粒子的所有维度。在GP=0.5的情况下，勘探和开发之间实现了良好的平衡。最后，EO的更新规则如下：
其中F在等式（11）中定义，V被认为是单位。
等式（16）中的第一项是平衡浓度，其中第二项和第三项表示浓度的变化。第二个术语负责全局搜索空间以找到最佳点。这个术语对探索有更大的贡献，从而利用了浓度的大变化（即平衡和样品颗粒之间的直接差异）。当它找到一个点时，第三项有助于使解决方案更加准确。因此，这个术语对开采有更大的贡献，并受益于由生成率术语控制的浓度的小变化（等式（13））。根据粒子浓度和平衡候选者以及周转率（λ）等参数，第二项和第三项可能具有相同或相反的符号。相同的符号使变异变大，这有助于更好地搜索整个域，相反的符号使变化变小，有助于局部搜索
尽管第二项试图找到离均衡候选者相对较远的解，而第三项试图细化更接近候选者的解，但这种情况并不总是发生。第三项分母中的小周转率（例如≤0.05）增加了其变化，也有助于在某些维度上进行探索。图1展示了这些术语如何促进勘探和开发的一维版本。C1−Ceq代表方程中的第二项。（16）而Ceq−λC1代表第三项（G是G0的函数）。生成速率项（等式（13）-（15））控制这些变化。因为λ随每个维度的变化而变化，所以这种大的变化只发生在λ值较小的维度上。值得一提的是，这一功能类似于进化算法中的突变算子，极大地帮助EO开发解决方案。
图2显示了在所提出的算法中，样本粒子上所有平衡候选者的协作以及它们如何一个接一个地影响浓度更新的概念草图。由于平衡候选者的拓扑位置在初始迭代中是不同的，并且指数项会产生大的随机数，因此这种逐步更新的过程有助于粒子在搜索中覆盖整个域。相反的情况发生在最后的迭代中，当候选者通过类似的配置围绕最优点时。在这些时候，指数项会产生较小的随机数，这有助于通过提供较小的步长来细化解决方案。这个概念也可以作为超空间扩展到更高的维度，由此浓度将随着粒子在n维空间中的运动而更新。
3.文件结构
EO.m % 平衡优化器算法
Get_Functions_details.m % 基准的全部信息和实现
initialization.m % 初始化
main.m % 主函数
4.伪代码
5.详细代码及注释
5.1 EO.m
function [Convergence_curve,Ave,Sd]=EO(Particles_no,Max_iter,lb,ub,dim,fobj,Run_no)
for irun=1:Run_no
Ceq1=zeros(1,dim); Ceq1_fit=inf;
Ceq2=zeros(1,dim); Ceq2_fit=inf;
Ceq3=zeros(1,dim); Ceq3_fit=inf;
Ceq4=zeros(1,dim); Ceq4_fit=inf;
C=initialization(Particles_no,dim,ub,lb);
Iter=0; V=1;
a1=2;
a2=1;
GP=0.5;
while Iter<Max_iter
for i=1:size(C,1)
Flag4ub=C(i,:)>ub;
Flag4lb=C(i,:)<lb;
C(i,:)=(C(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
fitness(i)=fobj(C(i,:));
if fitness(i)<Ceq1_fit
Ceq1_fit=fitness(i); Ceq1=C(i,:);
elseif fitness(i)>Ceq1_fit && fitness(i)<Ceq2_fit
Ceq2_fit=fitness(i); Ceq2=C(i,:);
elseif fitness(i)>Ceq1_fit && fitness(i)>Ceq2_fit && fitness(i)<Ceq3_fit
Ceq3_fit=fitness(i); Ceq3=C(i,:);
elseif fitness(i)>Ceq1_fit && fitness(i)>Ceq2_fit && fitness(i)>Ceq3_fit && fitness(i)<Ceq4_fit
Ceq4_fit=fitness(i); Ceq4=C(i,:);
end
end
%---------------- Memory saving-------------------
if Iter==0
fit_old=fitness; C_old=C;
end
for i=1:Particles_no
if fit_old(i)<fitness(i)
fitness(i)=fit_old(i); C(i,:)=C_old(i,:);
end
end
C_old=C; fit_old=fitness;
%-------------------------------------------------
Ceq_ave=(Ceq1+Ceq2+Ceq3+Ceq4)/4; % averaged candidate
C_pool=[Ceq1; Ceq2; Ceq3; Ceq4; Ceq_ave]; % Equilibrium pool
t=(1-Iter/Max_iter)^(a2*Iter/Max_iter); % Eq (9)
for i=1:Particles_no
lambda=rand(1,dim); % lambda in Eq(11)
r=rand(1,dim); % r in Eq(11)
Ceq=C_pool(randi(size(C_pool,1)),:); % random selection of one candidate from the pool
F=a1*sign(r-0.5).*(exp(-lambda.*t)-1); % Eq(11)
r1=rand(); r2=rand(); % r1 and r2 in Eq(15)
GCP=0.5*r1*ones(1,dim)*(r2>=GP); % Eq(15)
G0=GCP.*(Ceq-lambda.*C(i,:)); % Eq(14)
G=G0.*F; % Eq(13)
C(i,:)=Ceq+(C(i,:)-Ceq).*F+(G./lambda*V).*(1-F); % Eq(16)
end
Iter=Iter+1;
Convergence_curve(Iter)=Ceq1_fit;
Ceqfit_run(irun)=Ceq1_fit;
end
display(['Run no : ', num2str(irun)]);
display(['The best solution obtained by EO is : ', num2str(Ceq1,10)]);
display(['The best optimal value of the objective funciton found by EO is : ', num2str(Ceq1_fit,10)]);
disp(sprintf('--------------------------------------'));
end
Ave=mean(Ceqfit_run);
Sd=std(Ceqfit_run);
5.2 Get_Functions_details.m
function [lb,ub,dim,fobj] = Get_Functions_details(F)
switch F
case 'F1'
fobj = @F1;
lb=-100;
ub=100;
dim=30;
case 'F2'
fobj = @F2;
lb=-100;
ub=100;
dim=30;
case 'F3'
fobj = @F3;
lb=-100;
ub=100;
dim=30;
case 'F4'
fobj = @F4;
lb=-100;
ub=100;
dim=30;
case 'F5'
fobj = @F5;
lb=-30;
ub=30;
dim=30;
case 'F6'
fobj = @F6;
lb=-100;
ub=100;
dim=30;
case 'F7'
fobj = @F7;
lb=-1.28;
ub=1.28;
dim=30;
case 'F8'
fobj = @F8;
lb=-500;
ub=500;
dim=30;
case 'F9'
fobj = @F9;
lb=-5.12;
ub=5.12;
dim=30;
case 'F10'
fobj = @F10;
lb=-32;
ub=32;
dim=30;
case 'F11'
fobj = @F11;
lb=-600;
ub=600;
dim=30;
case 'F12'
fobj = @F12;
lb=-50;
ub=50;
dim=30;
case 'F13'
fobj = @F13;
lb=-50;
ub=50;
dim=30;
case 'F14'
fobj = @F14;
lb=-65.536;
ub=65.536;
dim=2;
case 'F15'
fobj = @F15;
lb=-5;
ub=5;
dim=4;
case 'F16'
fobj = @F16;
lb=-5;
ub=5;
dim=2;
case 'F17'
fobj = @F17;
lb=[-5,0];
ub=[10,15];
dim=2;
case 'F18'
fobj = @F18;
lb=-2;
ub=2;
dim=2;
case 'F19'
fobj = @F19;
lb=0;
ub=1;
dim=3;
case 'F20'
fobj = @F20;
lb=0;
ub=1;
dim=6;
case 'F21'
fobj = @F21;
lb=0;
ub=10;
dim=4;
case 'F22'
fobj = @F22;
lb=0;
ub=10;
dim=4;
case 'F23'
fobj = @F23;
lb=0;
ub=10;
dim=4;
end
end
% F1
function o = F1(x)
o=sum(x.^2);
end
% F2
function o = F2(x)
o=sum(abs(x))+prod(abs(x));
end
% F3
function o = F3(x)
dim=size(x,2);
o=0;
for i=1:dim
o=o+sum(x(1:i))^2;
end
end
% F4
function o = F4(x)
o=max(abs(x));
end
% F5
function o = F5(x)
dim=size(x,2);
o=sum(100*(x(2:dim)-(x(1:dim-1).^2)).^2+(x(1:dim-1)-1).^2);
end
% F6
function o = F6(x)
o=sum(abs((x+.5)).^2);
end
% F7
function o = F7(x)
dim=size(x,2);
o=sum([1:dim].*(x.^4))+rand;
end
% F8
function o = F8(x)
o=sum(-x.*sin(sqrt(abs(x))));
end
% F9
function o = F9(x)
dim=size(x,2);
o=sum(x.^2-10*cos(2*pi.*x))+10*dim;
end
% F10
function o = F10(x)
dim=size(x,2);
o=-20*exp(-.2*sqrt(sum(x.^2)/dim))-exp(sum(cos(2*pi.*x))/dim)+20+exp(1);
end
% F11
function o = F11(x)
dim=size(x,2);
o=sum(x.^2)/4000-prod(cos(x./sqrt([1:dim])))+1;
end
% F12
function o = F12(x)
dim=size(x,2);
o=(pi/dim)*(10*((sin(pi*(1+(x(1)+1)/4)))^2)+sum((((x(1:dim-1)+1)./4).^2).*...
(1+10.*((sin(pi.*(1+(x(2:dim)+1)./4)))).^2))+((x(dim)+1)/4)^2)+sum(Ufun(x,10,100,4));
end
% F13
function o = F13(x)
dim=size(x,2);
o=.1*((sin(3*pi*x(1)))^2+sum((x(1:dim-1)-1).^2.*(1+(sin(3.*pi.*x(2:dim))).^2))+...
((x(dim)-1)^2)*(1+(sin(2*pi*x(dim)))^2))+sum(Ufun(x,5,100,4));
end
% F14
function o = F14(x)
aS=[-32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32;,...
-32 -32 -32 -32 -32 -16 -16 -16 -16 -16 0 0 0 0 0 16 16 16 16 16 32 32 32 32 32];
for j=1:25
bS(j)=sum((x'-aS(:,j)).^6);
end
o=(1/500+sum(1./([1:25]+bS))).^(-1);
end
% F15
function o = F15(x)
aK=[.1957 .1947 .1735 .16 .0844 .0627 .0456 .0342 .0323 .0235 .0246];
bK=[.25 .5 1 2 4 6 8 10 12 14 16];bK=1./bK;
o=sum((aK-((x(1).*(bK.^2+x(2).*bK))./(bK.^2+x(3).*bK+x(4)))).^2);
end
% F16
function o = F16(x)
o=4*(x(1)^2)-2.1*(x(1)^4)+(x(1)^6)/3+x(1)*x(2)-4*(x(2)^2)+4*(x(2)^4);
end
% F17
function o = F17(x)
o=(x(2)-(x(1)^2)*5.1/(4*(pi^2))+5/pi*x(1)-6)^2+10*(1-1/(8*pi))*cos(x(1))+10;
end
% F18
function o = F18(x)
o=(1+(x(1)+x(2)+1)^2*(19-14*x(1)+3*(x(1)^2)-14*x(2)+6*x(1)*x(2)+3*x(2)^2))*...
(30+(2*x(1)-3*x(2))^2*(18-32*x(1)+12*(x(1)^2)+48*x(2)-36*x(1)*x(2)+27*(x(2)^2)));
end
% F19
function o = F19(x)
aH=[3 10 30;.1 10 35;3 10 30;.1 10 35];cH=[1 1.2 3 3.2];
pH=[.3689 .117 .2673;.4699 .4387 .747;.1091 .8732 .5547;.03815 .5743 .8828];
o=0;
for i=1:4
o=o-cH(i)*exp(-(sum(aH(i,:).*((x-pH(i,:)).^2))));
end
end
% F20
function o = F20(x)
aH=[10 3 17 3.5 1.7 8;.05 10 17 .1 8 14;3 3.5 1.7 10 17 8;17 8 .05 10 .1 14];
cH=[1 1.2 3 3.2];
pH=[.1312 .1696 .5569 .0124 .8283 .5886;.2329 .4135 .8307 .3736 .1004 .9991;...
.2348 .1415 .3522 .2883 .3047 .6650;.4047 .8828 .8732 .5743 .1091 .0381];
o=0;
for i=1:4
o=o-cH(i)*exp(-(sum(aH(i,:).*((x-pH(i,:)).^2))));
end
end
% F21
function o = F21(x)
aSH=[4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7;2 9 2 9;5 5 3 3;8 1 8 1;6 2 6 2;7 3.6 7 3.6];
cSH=[.1 .2 .2 .4 .4 .6 .3 .7 .5 .5];
o=0;
for i=1:5
o=o-((x-aSH(i,:))*(x-aSH(i,:))'+cSH(i))^(-1);
end
end
% F22
function o = F22(x)
aSH=[4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7;2 9 2 9;5 5 3 3;8 1 8 1;6 2 6 2;7 3.6 7 3.6];
cSH=[.1 .2 .2 .4 .4 .6 .3 .7 .5 .5];
o=0;
for i=1:7
o=o-((x-aSH(i,:))*(x-aSH(i,:))'+cSH(i))^(-1);
end
end
% F23
function o = F23(x)
aSH=[4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7;2 9 2 9;5 5 3 3;8 1 8 1;6 2 6 2;7 3.6 7 3.6];
cSH=[.1 .2 .2 .4 .4 .6 .3 .7 .5 .5];
o=0;
for i=1:10
o=o-((x-aSH(i,:))*(x-aSH(i,:))'+cSH(i))^(-1);
end
end
function o=Ufun(x,a,k,m)
o=k.*((x-a).^m).*(x>a)+k.*((-x-a).^m).*(x<(-a));
end
5.3 initialization.m
function [Cin,domain]=initialization(SearchAgents_no,dim,ub,lb)
Boundary_no= size(ub,2); % numnber of boundaries
% If the boundaries of all variables are equal and user enter a signle
% number for both ub and lb
if Boundary_no==1
Cin=rand(SearchAgents_no,dim).*(ub-lb)+lb;
domain=ones(1,dim)*(ub-lb);
end
% If each variable has a different lb and ub
if Boundary_no>1
for i=1:dim
ub_i=ub(i);
lb_i=lb(i);
Cin(:,i)=rand(SearchAgents_no,1).*(ub_i-lb_i)+lb_i;
end
domain=ones(1,dim).*(ub-lb);
end
5.4 main.m
clear all
clc
tic;
Run_no=30; % Number of independent runs
Particles_no=30; % Number of particles
Max_iteration=500; % Maximum number of iterations
Function_name='F1';
[lb,ub,dim,fobj]=Get_Functions_details(Function_name);
[Convergence_curve,Ave,Sd]=EO(Particles_no,Max_iteration,lb,ub,dim,fobj,Run_no);
display(['The average objective function is : ', num2str(Ave,7)]);
display(['The standard deviation is : ', num2str(Sd,7)]);
toc;
6.运行结果
7.参考文献
[1]Faramarzi A,Heidarinejad M,Stephens B, et al. Equilibrium optimizer: A novel optimization algorithm[J]. Knowledge-Based Systems,2020,191©.
