混合蛙跳算法

1.算法介绍

本算法参考：https://wenku.baidu.com/view/f91e94ec011ca300a6c390b4.html

2.算法步骤描述

2.1算法流程图

图片来自：https://wenku.baidu.com/view/3ddad81b4b7302768e9951e79b89680203d86b12.html

3.算法代码

%SLAFA.m文件
%
% F:青蛙种群数量(pupolation)
% M:青蛙族群(memeplexes)
% N:每个族群中青蛙的数量
% max_iteration:满足混合最大迭代次数
% NN:每个族群中最大进化次数
% q:构建子族群时选择的青蛙数量
% SMAX:子族群迭代时最大步长
% fobj:目标函数
% dim:目标函数维度
% la:目标函数求解区间最左端
% lb:目标函数求解区间最右端
%
%
function BestCost = MySLFA(F,M,N,max_iteration,NN,q,SMAX,fobj,dim,la,lb)
    %0.初始化。种群大小，族群数量以及族群中青蛙的数量。(F = M*N)
    %1.生成一个初始种群
    empty_individual.Position =[];%变量的随机取值
    empty_individual.Cost = [];%适应度函数值
    Pop = repmat(empty_individual,[F,1]);%提高运算速度
    %初始化种群
    for i = 1:F
        Pop(i).Position = unifrnd(la,lb,[1 dim]);%初始位置
        Pop(i).Cost = fobj(Pop(i).Position);%适应度函数
    end
    %2.对青蛙种群按照适应度函数进行排序
    Pop = SortPupolation(Pop);
    %此时记录下全局最优为PX
    PX = Pop(1);
    %定义一个变量记录求解得到的最优值
    BestCost = nan(max_iteration,1);%预先分配内存,提高运算速度
    for it = 1:max_iteration
        %3.对种群进行划分，分成M个族群，每个族群中青蛙数量为N
        I = reshape(1:F,M,N);%将种群分成M个组   
        %4.对每个族群实行局部最优化搜索
        for im = 1:M
            %对第im族群进行局部搜索
            mPop = Pop(I(im,:));
            %第im族群中进行最大进化次数为NN
            for in = 1:NN
                %4-3.按照概率权重选择q只青蛙组成子族群
                mPopNum = numel(mPop);%第im中群中青蛙的数量,即为N
                P = 2*(mPopNum+1-(1:mPopNum))/(mPopNum*(mPopNum+1));
                L = RandSelectFrogByProbility(P,q);
                subPop = mPop(L);
                %排序,记录最好值PB和最坏值PW
                [subPop,ind] = SortPupolation(subPop);
                L = L(ind);
                PB = subPop(1).Position;
                PW = subPop(end).Position;

                %定义标记符号
                ImprovementStep2 = false;%表示不需要进行第二步改良,即不需要执行4-5步骤
                Censorship = false;%表示已经产生更好的值,不需要执行4-6步骤
                %4-4.更新最坏值的位置,并检查位置的值是否在[la,lb]中
                %本步骤可以看做改良第一步
                NewFrog = subPop(end);
                NewFrog.Position = PW+CalculationDeltaS(PB,PW,SMAX);
                if IsInInteral(NewFrog.Position,la,lb)
                    %更新后的值在区间[la,lb]中,计算适应度函数值
                    NewFrog.Cost = fobj(NewFrog.Position);
                    if NewFrog.Cost < subPop(end).Cost 
                        %新产生的值比原来的值好，则用新青蛙替换原来的青蛙
                        subPop(end) = NewFrog;
                        %去4-7

                    else
                        %去4-5
                        ImprovementStep2 = true;
                    end
                else
                    %去4-5
                    ImprovementStep2 = true;
                end
                %4-5.如果4-4不能产生更好的结果,才执行本步骤
                %本步骤可以看做改良第二步
                if ImprovementStep2
                    NewFrog2 = subPop(end);
                    %使用全局最优重新计算步长增量
                    NewFrog2.Position = PW+CalculationDeltaS(PX.Position,PW,SMAX);
                    if IsInInteral(NewFrog2.Position,la,lb)
                        %更新后的值在区间[la,lb]中,计算适应度函数值
                        NewFrog2.Cost = fobj(NewFrog2.Position);
                        if NewFrog2.Cost < subPop(end).Cost 
                             %新产生的值比原来的值好，则用新青蛙替换原来的青蛙
                            subPop(end) = NewFrog2;
                            %去4-7

                        else
                            %去4-6
                            Censorship = true;
                        end
                    else
                        %去4-6
                        Censorship = true;
                    end
                end%end ImprovementStep2

                %4-6.上述方法不能产生更好的结果，则重新生成一只青蛙
                if Censorship
                    %在可行位置的地方重新生成青蛙的位置
                    [xx,yy]= FindALocation(mPop);%计算当前族群中的可行位置
                    subPop(end).Position = unifrnd(xx,yy);
                    subPop(end).Cost = fobj(subPop(end).Position);
                end%end Censorship
                %4-7.更新族群,并按照适应度函数排序
                mPop(L) = subPop;
                SortPupolation(mPop);%更新排序
            end%end in
            %更新种群的青蛙
             Pop(I(im,:)) = mPop;
        end%end im

        %5.对种群重新洗牌,再次划分
        Pop = SortPupolation(Pop);
        %更新全局最优解
        PX = Pop(1);
        %6.检查收敛性
        BestCost(it) = PX.Cost;
    end
end
%计算当前族群中的可行位置
function[x,y] = FindALocation(mPop)
    nPop = numel(mPop);
    x = mPop.Position;
    y = mPop.Position;
    for i = 2:nPop
        x = min(x, mPop(i).Position);
        y = max(y, mPop(i).Position);
    end
end
%判断向量x的值是否在区间[la,lb]中
function flag = IsInInteral(x,la,lb)
    flag = all(x >= la) && all(x <= lb);
end
%计算子族群每次更新时的步长增量
function s = CalculationDeltaS(PB,PW,SMAX)
    %PB是最好的位置
    %PW是最坏的位置
    %SMAX为最大步长
    LL = numel(PB);
    siga = rand();%生成一个(0-1)之间的随机数
    s = zeros(size(PB));%接收步长结果
    for i = 1:LL
        t = PB(i)-PW(i);
        if t>=0
            %min{int[rand(PB(i) ? PW(i))],SMAX}
            res = min(floor(siga*t),SMAX);
            %res = min(siga*t,SMAX);
            s(i) = res;
        else
            %max{int[rand(PB(i) ? PW(i))],?SMAX}
            res = max(ceil(siga*t),-SMAX);
            %res = max(siga*t,-SMAX);
            s(i) = res;
        end
    end
end
%根据概率权重选择q只青蛙组成子族群
function L = RandSelectFrogByProbility(P,q)
    L = zeros(q,1);
    for i = 1:q
        L(i) = randsample(numel(P),1, true, P);
    end
end
%种群按照适应度函数从小到大排序
function[pop,index] = SortPupolation(pop)
    cost = [pop.Cost];
    [~,index] = sort(cost);
    pop = pop(index);   
end

%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=-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=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;
        ub=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

%main.m文件
clear;clc;
% F:青蛙种群数量(pupolation)
% M:青蛙族群(memeplexes)
% N:每个族群中青蛙的数量
% max_iteration:满足混合最大迭代次数
% NN:每个族群中最大进化次数
% q:构建子族群时选择的青蛙数量
% SMAX:子族群迭代时最大步长
% fobj:目标函数
% dim:目标函数维度
% la:目标函数求解区间最左端
% lb:目标函数求解区间最右端

%F1
[lb,ub,dim,fobj] = Get_Functions_details('F3');
BestCost = MySLFA(100,5,20,1000,3,6,5,fobj,dim,lb,ub);
semilogy(BestCost, 'LineWidth', 2);
xlabel('Iteration');
ylabel('Best Cost');
grid on;

注意：本算法是根据论文(Muzaffar Eusuff , Kevin Lansey & Fayzul Pasha)Shuffled frog-leaping algorithm: a memetic metaheuristic for discrete optimization中提出的步骤进行复现。欢迎大家提出问题，若有问题欢迎联系:ronggong20@qq.com。