几何平均优化器(Geometric Mean Optimizer, GMO)算法简介及matlab代码

最新推荐文章于 2024-04-08 08:54:05 发布

Lzh_023016

最新推荐文章于 2024-04-08 08:54:05 发布

阅读量954

点赞数

分类专栏： Matlab 文章标签：算法 matlab 开发语言

本文链接：https://blog.csdn.net/Lzh_023016/article/details/132291524

版权

Matlab 专栏收录该内容

20 篇文章 3 订阅

订阅专栏

几何平均优化器(Geometric Mean Optimizer, GMO)

几何平均优化器是一种新的元启发式算法，它模拟了数学中几何平均算子的独特性质。该算子可以同时评估搜索空间中搜索个体的适应度和多样性。在GMO中，在求解优化问题时，将某个体的对立面的标度目标值的几何平均值分配给该个体作为其权重，代表其在搜索过程中指导其他个体的总资格。
在这里插入图片描述

GMO算法的数学模型

GMO是一种元启发式算法，利用多个搜索个体在相互社交互动时的行为。必须在任何优化算法中描述这些个体应该如何相互作用以产生最佳结果。
$X_{i} =\left (x _{i1}, x _{i2}, \cdots ,x _{iD}\right ) ，V_{i} =\left (v _{i1}, v _{i2}, \cdots ,v _{iD}\right )$
上式分别表示种群个体的位置及速度。
为了解决优化问题，GMO可以通过单一指标同时评价搜索个体在搜索空间中的适应度和多样性，称为双适应度指标 (dual-fitness index, DFI) 。为了计算DFI，参照模糊逻辑的概念，几个变量之间的相似性可以通过将归属于这些变量的模糊隶属函数(MF)值相乘来表示。这种乘法被称为基于积的Larsen隐函数，广泛应用于模糊系统。n个隶属度(MD)的几何平均值可以表示为：
$\sqrt[n]{\mu _{1} \times\mu _{2} \times\cdots \times \mu _{n} }$
为了定义MFs，利用下面的公式：
$f\left ( x \right ) =\frac{1}{1+e^{a\left ( x-c \right ) } } ;a< 0$

将x代入个人迄今最佳搜索个体计算出的目标函数值，可根据式计算出该个体的模糊值MF：
$MF_{j}^{t} =\frac{1}{1+exp\left [ -\frac{4}{\sigma ^{t}\sqrt{e} } \times \left ( Z_{best,j}^{t} -\mu ^{t} \right ) \right ] };j=1,2,\cdots ,N$
然后，为搜索个体计算一个新的指标，我们将其命名为双适应度指数(dual-fitness index, DFI)：
$DFI_{i}^{t} =MF_{1}^{t} \times \cdots \times MF_{i-1}^{t}\times MF_{i+1}^{t}\times \cdots \times MF_{N}^{t} = {\textstyle \prod_{j=1,j\neq i}^{N}MF_{j}^{t} }$

为了使所有个体最佳到目前为止的个体都有助于为每个个体生成唯一的全局引导个体，定义了所有相反的个体最佳到目前为止的个体的加权平均值，其权重为分配给它们的DFI值。
$Y_{i}^{t} =\frac{\sum_{j=1,j\neq i}^{N}DFI_{j}^{t}\times X_{j}^{best} }{\sum_{j=1}^{N}DFI_{j}^{t}+\varepsilon }$
因此，当在寻找引导个体时解决精英主义问题时，可转换为：
$Y_{i}^{t} =\frac{\sum_{j\in Nbest,j\neq i}^{N}DFI_{j}^{t}\times X_{j}^{best} }{\sum_{j\in Nbest}^{N}DFI_{j}^{t}+\varepsilon }$
为了使Y具有更强的随机性，更好地保护Y的多样性，这些个体在GMO中发生突变。突变方案选择高斯突变。将这种形式的突变施加到个体上。
$Y_{i,mut}^{t} =Y_{i}^{t} +\omega \times randn\times \left ( Std_{max}^{t} -Std^{t} \right )$
最后，根据下式推导出各搜索个体的更新方程：
$V_{i}^{t+1} =\omega \times V_{i}^{t} +\varphi \times \left ( Y_{i,mut}^{t} -X_{i}^{t} \right )$
$\varphi =1+\left ( 2\times randn-1 \right ) \times \omega$
$X_{i}^{t+1} =X_{i}^{t} +V_{i}^{t+1}$
$\omega =1-\frac{t}{t_{max} }$

GMO算法的流程图

在这里插入图片描述

GMO算法的Matlab代码

（1）调用程序

% geometric mean optimizer(GMO)
% Reference:GMO: geometric mean optimizer for solving engineering problems.
% Code by Luzhenhui Yangtze University

clc
clear
tic
run=1;                                             % GWO算法独立运行次数，一般为30次。
np=50;                                             % 种群个体数目
maxit=1000;                                        % 算法的迭代次数
Function_name='F1';                                % 优化目标函数 F1~F23
epsilon=0;                                         % 其值根据目标函数的类型确定
[lb,ub,nx,fobj]=Objective_Function(Function_name); % 获取目标函数的相关信息
varmax=ub*ones(1,nx);                              % 变量取值范围上界（矩阵：1*nx）
varmin=lb*ones(1,nx);                              % 变量取值范围下界（矩阵：1*nx）
limvel=0.1;                                        % 搜索空间中的最大距离与形成最大速度的比率 
velmax=limvel*(varmax(1,1:nx)-varmin(1,1:nx));     % 速度的上界
velmin=-velmax;                                    % 速度的下界
z_iter_main=zeros(run,maxit);
z_final_main=zeros(run);
pos_final_main=zeros(run,nx);
x1=zeros(maxit);
y1=zeros(maxit);
% 反复运行GWO，“run”次
for nrun=1:run
    [z_iter,z_final,pos_final]=GMO(np,nx,maxit,varmax,varmin,velmax,velmin,epsilon,fobj);
     z_iter_main(nrun,1:maxit)=z_iter(1:maxit);
     z_final_main(nrun)=z_final;
     pos_final_main(nrun,1:nx)=pos_final(1:nx);
end
% 显示完整的结果
disp(['The final statistical results calculated when implementing the GMO  for ',num2str(run),' times are as follows:']);
disp(['The average of the final objective function values calculated over ',num2str(run),' times = ',num2str(mean(z_final_main(1:run)))]);
disp(['The median of the final objective function values calculated over ',num2str(run),' times = ',num2str(median(z_final_main(1:run)))]);
disp(['The best of the final objective function values calculated over ',num2str(run),' times = ',num2str(min(z_final_main(1:run)))]);
disp(['The standard deviation of the final objective function values calculated over ',num2str(run),' times = ',num2str(std(z_final_main(1:run)))]);
% 绘制GWO的迭代收敛曲线
for i=1:maxit
    x1(i)=i;sum1=0;
    for j=1:run
        sum1=sum1+z_iter_main(j,i);
    end
    y1(i)=sum1/run;
end
semilogy(x1,y1,'-r')
xlabel('Iteration');
ylabel('Average best-so-far');
legend('GMO');
hold on
time_gmo = toc;
disp(['Elapsed time of running the GMO for ',num2str(run),' times = ',num2str(time_gmo),' seconds']);

（2）GMO算法主程序

% geometric mean optimizer(GMO)
% Reference:GMO: geometric mean optimizer for solving engineering problems.
% Code by Luzhenhui Yangtze University

% 输出
% z_iter：每一次迭代的最优适应度值，其可以用于迭代收敛图的绘制
% z_final：全局最优适应度值
% pos_final：全局最优适应度值所对应的最优个体
% 输入
% np：算法种群个体数目
% nx：变量的维度
% maxit：算法的最大迭代次数
% varmax：变量的取值范围上界(1*nx)
% velmax：定义速率的上限
% velmin：定义速率的下限
% epsilon：对于优化目标函数的类型，进行不同的赋值
% fobj：优化的目标函数
function [z_iter,z_final,pos_final] = GMO(np,nx,maxit,varmax,varmin,velmax,velmin,epsilon,fobj)
pp_pbest=zeros(np,nx);
pp_kbest=zeros(np,nx);
pp_guide=zeros(np,nx);
mutant=zeros(1,nx);
stdev2=zeros(nx);
index=zeros(np);
fit=zeros(maxit,np);
DFI=zeros(np);
optimal_pos=zeros(1,nx);
z_pbest=zeros(np);
z_kbest=zeros(np);
z_optimal=inf*ones(maxit);
pos_final=zeros(nx);
z_iter=zeros(maxit);
kbest_max=np;
kbest_min=2;
% Starting the optimization process
it=1;
% 初始化种群个体及“速率”
[pp,pv]=Initialization(np,nx,varmax,varmin,velmax,velmin);
% 计算种群个体的适应度值并确定最优值及其对应的最优个体
for j=1:np
    z=fobj(pp(j,1:nx));
    z_pbest(j)=z;
    pp_pbest(j,1:nx)=pp(j,1:nx);
end
for i=1:nx
    stdev2(1,i)=std(pp_pbest(1:np,i));
end
max_stdev2=max(stdev2(1,1:nx));
ave=mean(z_pbest(1:np));        % 计算目前个体的平均值
stdev=std(z_pbest(1:np));       % 计算目前个体的方差
kbest=kbest_max-(kbest_max-kbest_min)*(it/maxit);
n_best=round(kbest);
% 公式(3)
for j=1:np
    index(j)=j;
end
for j=1:np
    prod=1;
    for jj=1:np
        if jj~=j
            prod=prod*(1/(1+exp((-4)/(stdev*sqrt(exp(1)))*(z_pbest(jj)-ave)))); 
        end
    end
    fit(it,j)=prod;
end
for j=1:np-1
    for jj=j+1:np
        if fit(it,jj)>fit(it,j)
            c1=fit(it,j);
            fit(it,j)=fit(it,jj);
            fit(it,jj)=c1;
            c2=index(j);
            index(j)=index(jj);
            index(jj)=c2;
        end
    end
end
% 精英解
sum1=0;
for j=1:n_best
    z_kbest(j)=z_pbest(index(j));
    pp_kbest(j,1:nx)=pp_pbest(index(j),1:nx);
    sum1=sum1+fit(it,j);
end 
% 公式（5）
for j=1:np
    pp_guide(j,1:nx)=zeros(1,nx);
    for jj=1:n_best
        if index(jj)~=j
            DFI(jj)=fit(it,jj)/(sum1+epsilon);
            pp_guide(j,1:nx)=pp_guide(j,1:nx)+DFI(jj).*pp_kbest(jj,1:nx);
        end
    end
end
% 确定最优值及其对应的最优个体
for j=1:np
    if z_pbest(j)<z_optimal(it)
        z_optimal(it)=z_pbest(j);
        optimal_pos(it,1:nx)=pp_pbest(j,1:nx);
    end
end
z_iter(it)=z_optimal(it); % 保存当前迭代的全局最优适应度值
% 算法主要迭代步骤
while it<maxit
    it=it+1;        
    w=1-(it/maxit); % 公式(9)
    for j=1:np 
        mutant(1,1:nx)=pp_guide(j,1:nx)+w*randn(1,nx).*(max_stdev2-stdev2(1,1:nx));   % 公式(6)
        pv(j,1:nx)=w*pv(j,1:nx)+(ones(1,nx)+(2*rand(1,nx)-ones(1,nx))*w).*(mutant(1,1:nx)-pp(j,1:nx));   % 公式(7)
        % 速度的边界处理（确保在规定的取值范围内）
        flag4lbv=pv(j,:)<velmin(1,:);
        flag4ubv=pv(j,:)>velmax(1,:);
        pv(j,:)=(pv(j,:)).*(~(flag4lbv+flag4ubv))+velmin.*flag4lbv+velmax.*flag4ubv;
        
        pp(j,:)=pp(j,:)+pv(j,:);    % 公式(8)
        
        % 种群个体的边界处理
        flag4lbp=pp(j,:)<varmin(1,:);
        flag4ubp=pp(j,:)>varmax(1,:);
        pp(j,:)=(pp(j,:)).*(~(flag4lbp+flag4ubp))+varmin.*flag4lbp+varmax.*flag4ubp; 
        pv(j,:)=(pv(j,:)).*(ones(1,nx)-2*(flag4lbp+flag4ubp)); 
        % 计算当前种群个体的适应度值，确定最优值及其相对应的最优个体
        z=fobj(pp(j,1:nx));
        if z<z_pbest(j)
            z_pbest(j)=z;
            pp_pbest(j,1:nx)=pp(j,1:nx);
        end
    end
    for i=1:nx
        stdev2(1,i)=std(pp_pbest(1:np,i));
    end
    max_stdev2=max(stdev2(1,1:nx));
    % 计算均值及方差
    ave=mean(z_pbest(1:np));
    stdev=std(z_pbest(1:np));
    % 确定精英迄今为止最佳解决方案的数量
    kbest=kbest_max-(kbest_max-kbest_min)*(it/maxit);
    n_best=round(kbest);
    % 公式(3)
    for j=1:np
        index(j) = j;
    end
    for j=1:np
        prod=1;
        for jj=1:np
            if jj~=j
                prod=prod*(1/(1+exp((-4)/(stdev*sqrt(exp(1)))*(z_pbest(jj)-ave))));
            end
        end
        fit(it,j)=prod;
    end
    for j=1:np-1
        for jj=j+1:np
            if fit(it,jj)>fit(it,j)
                c1=fit(it,j);
                fit(it,j)=fit(it,jj);
                fit(it,jj)=c1;
                c2=index(j);
                index(j)=index(jj);
                index(jj)=c2;
            end
        end
    end
    sum1=0;
    for j=1:n_best
        z_kbest(j)=z_pbest(index(j));
        pp_kbest(j,1:nx)=pp_pbest(index(j),1:nx);
        sum1=sum1+fit(it,j);
    end
    for j=1:np
        pp_guide(j,1:nx)=zeros(1,nx);
        for jj=1:n_best
            if index(jj)~=j
                DFI(jj)=fit(it,jj)/(sum1+epsilon);
                pp_guide(j,1:nx)=pp_guide(j,1:nx)+DFI(jj).*pp_kbest(jj,1:nx); % Eq.(5)
            end
        end
    end
    % 确定目前的最优适应度值及其对应的个体
    for j=1:np
        if z_pbest(j)<z_optimal(it)
            z_optimal(it)=z_pbest(j);
            optimal_pos(it,1:nx)=pp_pbest(j,1:nx);
        end
    end
    % 保存当前迭代的最优适应度值
    z_iter(it)=z_optimal(it);
end
z_final=z_iter(maxit);
pos_final(1:nx)=optimal_pos(maxit,1:nx);
end

（3）算法初始化

% geometric mean optimizer(GMO)
% Reference:GMO: geometric mean optimizer for solving engineering problems.
% Code by Luzhenhui Yangtze University

% 初始化种群个体及速度
function [pp,pv] = Initialization(np,nx,varmax,varmin,velmax,velmin)
pp=zeros(np,nx);
pv=zeros(np,nx);

for j=1:np
    pp(j,1:nx)=(varmax-varmin).*rand(1,nx)+varmin;
    pv(j,1:nx)=(velmax-velmin).*rand(1,nx)+velmin;
end
end

（4）测试函数

% geometric mean optimizer(GMO)
% Reference:GMO: geometric mean optimizer for solving engineering problems.
% Code by Luzhenhui Yangtze University

function [lb,ub,dim,fobj] = Objective_Function(F)

switch F
    case 'F1'
        fobj = @F1;
        lb=-100;
        ub=100;
        dim=30;
        
    case 'F2'
        fobj = @F2;
        lb=-10;
        ub=10;
        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=100;
        
    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)
so = 0;
x = x - so*ones(1,size(x,2));
o=sum(x.^2);
end

% F2

function o = F2(x)
so = 0;
x = x - so*ones(1,size(x,2));
o=sum(abs(x))+prod(abs(x));
end

% F3

function o = F3(x)
so = 0;
x = x - so*ones(1,size(x,2));    
dim = size(x,2);
o=0;
for i=1:dim
    o=o+sum(x(1:i))^2;
end
end

% F4

function o = F4(x)
so = 0;
x = x - so*ones(1,size(x,2));
o=max(abs(x));
end

% F5

function o = F5(x)
so = 0;
x = x - so*ones(1,size(x,2));
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)
so = 0;
x = x - so*ones(1,size(x,2));
o=sum(floor((x+.5)).^2);
end

% F7

function o = F7(x)
so = 0;
x = x - so*ones(1,size(x,2));
dim=size(x,2);
o=sum([1:dim].*(x.^4))+rand;
end

% F8

function o = F8(x)
so = 0;
x = x - so*ones(1,size(x,2));
o=sum(-x.*sin(sqrt(abs(x))));
end

% F9

function o = F9(x)
so = 0;
x = x - so*ones(1,size(x,2));
dim=size(x,2);
o=sum(x.^2-10*cos(2*pi.*x))+10*dim;
end

% F10

function o = F10(x)
so = 0;
x = x - so*ones(1,size(x,2));
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)
so = 0;
dim=size(x,2);
x = x - so*ones(1,dim);
o=sum(x.^2)/4000-prod(cos(x./sqrt([1:dim])))+1;
end

% F12

function o = F12(x)
so = 0;
x = x - so*ones(1,size(x,2));
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)
so = 0;
x = x - so*ones(1,size(x,2));
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