用于全局优化的分数阶Lèvy飞行蝙蝠算法源码复现

function [bestX,fMin] = FLFBA()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 设置默认参数
M = 1000;
pop = 10;
dim = 2;
gama = 0.9; %0.9
alpha = 0.99; %0.99
r0Max = 1; %1
r0Min = 0; %0
AMax = 2; %2
AMin = 0; %0.5
freqMax = 2; %2
freqMin = 0; %0.1
%%%%%%%%%%%%%%%%%%%%%%%
wMax = 0.9; %0.9
wMin = 0.2; %0.2
xsi_init = 0.6; %0.6
n = 2; %2
% 参数设置
lb= -5 * ones(1,dim );   % 最低边界
ub= 5 * ones(1,dim );    % 最高边界

alfa=0.632; 
beta=3/2;   % Lèvy飞行轨迹参数
sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*...
    beta*2^((beta-1)/2)))^(1/beta); % Lèvy飞行轨迹参数

Vhist = zeros(10,pop,dim);  % 速度的历史值保存列表
vLb = 0.6 * lb; %0.6
vUb = 0.6 * ub; %0.6
r = rand( pop, 1 ) .* 0.2 + 0; 
r0 = rand( pop, 1 ) .* ( r0Max - r0Min ) + r0Min;
A = rand( pop, 1 ) .* ( AMax - AMin ) + AMin;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initialization
for i = 1 : pop
    x( i, : ) = lb + (ub - lb) .* rand( 1, dim );
    v( i, : ) = rand( 1, dim );
    fit( i ) = Fobt(x(i,:)');
end
pFit = fit; % 个体最优适应度值
pX = x;     % 个体最优位置

[ fMin, bestIndex ] = min( fit );  % fMin表示全局最优
                   % bestX表示fMin对应的位置
bestX = x( bestIndex, : );
Vhist(1,:,:)= v;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 迭代开始
for iteration = 1 : M
    %%%%%%%%%%%%%%%%%%%%%% 迭代参数%%%%%%%%%%%%%%%%%%%%%%%%%%%
    freq = rand( pop, dim ) .* ( freqMax - freqMin ) + freqMin;
    w = (wMax - wMin) * ( M - iteration )/(1.0 * M) + wMin; %Inertia weight
    xsi1 = 1 + ((xsi_init-1)*((M-iteration)/M)^n);
    xsi2 = 1-xsi1;
    meanA = mean( A );
    meanP = mean(pX);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i = 1 : pop
        if rand < 0.5 %0.5
            q1 = randi([1 pop]);
            q2 = randi([1 pop]);
            X1 = pX(q1,:);
            X2 = pX(q2,:);
            newX=kill_bat(pX(i,:),X1,X2); % Eq(24-25)
            [pX(i,:),pFit(i)]=select_bat(pX(i,:),newX,pFit(i));
            %%%%%%%%%%%% Levy飞行轨迹 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
            u=randn(1,dim)*sigma;
            VV=randn(1,dim);
            step=u./abs(VV).^(1/beta);
            stepsize = 0.01*step.*randn(1,dim); %0.01
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            x( i, : ) = pX(i,:) + stepsize .* ... %pX(i,:)
                abs(meanP - x(i,:));    %Eq(26)           
        else   % 使用分数导数
            t1 = randi([1 10]);
            q = randi([1 pop]);
            while q== i
                q = randi([1 pop]);
            end
            %%%%%%%%%%%%%%% fractional derivative %%%%%%%%%%%%%%%%%%
            Vh = reshape(Vhist(:,i,:),10,dim); 
            vout = Frac_Diff_Der(alfa,Vh,t1); %Eq(27)
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            v(i,:) = w.* vout + freq(i,:) .* xsi1 .*(bestX-x(i,:))+...
                freq(i,:).* xsi2 .*(pX(i,:)-x(q,:)); %Eq(22) 
            v(i,:) = Bounds(v(i,:),vLb,vUb);
            x(i,:) = x(i,:) + v(i,:);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        end
        
        % 局部搜索
        if rand > r(i)
            %%%%%%%%%%%%%%%% 另外一种Levy飞行轨迹 %%%%%%%%%%%%%%%%%%%%
            u=randn(1,dim)*sigma;
            VV=randn(1,dim);
            step=u./abs(VV).^(1/beta);
            stepsize = 0.01*step;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            randnValueA = stepsize .*randn(1,dim).* ...
                (abs(A(i)-meanA)+realmin);
            x(i,:) = bestX .* (1+randnValueA); %Eq(32)
        end
        x(i,:) = Bounds(x(i,:),lb,ub);
        fit(i) = Fobt(x(i,:)');
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 更新个人最优适应度值和全局最优个体
    
    for i = 1 : pop
        if fit( i ) < pFit( i )
            pFit( i ) = fit( i );
            pX( i, : ) = x( i, : );
        end
        
        if( pFit( i ) < fMin && rand < A(i) )
            fMin = pFit( i );
            bestX = pX( i, : );
            A(i) = A(i) * alpha;
            r(i) = r0(i) * ( 1 - exp( -gama * iteration ) );
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    Vhist(2:10,:,:) = Vhist(1:9,:,:);
    Vhist(1,:,:) = v;
    
    fbest = fMin; %#ok<NASGU>
    FunMin(iteration) = fMin;
end
% End of the main program

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 以下函数与主程序相关
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Application of simple limits/bounds
function s = Bounds( s, Lb, Ub)
% 简单极限/界限的应用
temp = s;
I = temp < Lb;
temp(I) = Lb(I);

% 应用上界向量
J = temp > Ub;
temp(J) = Ub(J);
% 更新此移动
s = temp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X,fits]=select_bat(X,newX,fits)
% 评价所有个体

   fnew= Fobt(newX');
    if fnew<=fits
        fits =fnew;
        X=newX;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function newX=kill_bat(X,X1,X2)
 % 一小部分更差的点以一定概率被发现
 stepsize=rand*(X1-X2);
 newX=X+stepsize;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [vout] = Frac_Diff_Der(alfa,V,r)

switch(r)
    case 1
        vout = alfa*V(1,:);
    case 2
        vout = alfa*V(1,:)+(1/factorial(2))*V(2,:);
    case 3 
        vout = alfa*V(1,:)+(1/factorial(2))*V(2,:)+...
            (1/factorial(3))*alfa*(1-alfa)*V(3,:);
    case 4
        vout = alfa*V(1,:)+(1/factorial(2))*V(2,:)+...
            (1/factorial(3))*alfa*(1-alfa)*V(3,:)+ ...
            (1/factorial(4))*alfa*(1-alfa)*(2-alfa)*V(4,:);
    case 5
        vout = alfa*V(1,:)+(1/factorial(2))*V(2,:)+...
            (1/factorial(3))*alfa*(1-alfa)*V(3,:)+ ...
            (1/factorial(4))*alfa*(1-alfa)*(2-alfa)*V(4,:)+ ...
            (1/factorial(5))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*V(5,:);
    case 6
        vout = alfa*V(1,:)+(1/factorial(2))*V(2,:)+...
            (1/factorial(3))*alfa*(1-alfa)*V(3,:)+ ...
            (1/factorial(4))*alfa*(1-alfa)*(2-alfa)*V(4,:)+ ...
            (1/factorial(5))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*V(5,:)+...
            (1/factorial(6))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*V(6,:);
    case 7
       vout = alfa*V(1,:)+(1/factorial(2))*V(2,:)+...
            (1/factorial(3))*alfa*(1-alfa)*V(3,:)+ ...
            (1/factorial(4))*alfa*(1-alfa)*(2-alfa)*V(4,:)+ ...
            (1/factorial(5))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*V(5,:)+...
            (1/factorial(6))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*V(6,:)+ ...
            (1/factorial(7))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*(5-alfa)*V(7,:);
    case 8
        vout = alfa*V(1,:)+(1/factorial(2))*V(2,:)+...
            (1/factorial(3))*alfa*(1-alfa)*V(3,:)+ ...
            (1/factorial(4))*alfa*(1-alfa)*(2-alfa)*V(4,:)+ ...
            (1/factorial(5))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*V(5,:)+...
            (1/factorial(6))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*V(6,:)+ ...
            (1/factorial(7))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*(5-alfa)*V(7,:)+...
            (1/factorial(8))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*(5-alfa)*(6-alfa)*V(8,:);
    case 9 
        vout = alfa*V(1,:)+(1/factorial(2))*V(2,:)+...
            (1/factorial(3))*alfa*(1-alfa)*V(3,:)+ ...
            (1/factorial(4))*alfa*(1-alfa)*(2-alfa)*V(4,:)+ ...
            (1/factorial(5))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*V(5,:)+...
            (1/factorial(6))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*V(6,:)+ ...
            (1/factorial(7))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*(5-alfa)*V(7,:)+...
            (1/factorial(8))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*(5-alfa)*(6-alfa)*V(8,:)+...
            (1/factorial(9))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*(5-alfa)*(6-alfa)*(7-alfa)*V(9,:);
    case 10
        vout = alfa*V(1,:)+(1/factorial(2))*V(2,:)+...
            (1/factorial(3))*alfa*(1-alfa)*V(3,:)+ ...
            (1/factorial(4))*alfa*(1-alfa)*(2-alfa)*V(4,:)+ ...
            (1/factorial(5))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*V(5,:)+...
            (1/factorial(6))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*V(6,:)+ ...
            (1/factorial(7))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*(5-alfa)*V(7,:)+...
            (1/factorial(8))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*(5-alfa)*(6-alfa)*V(8,:)+...
            (1/factorial(9))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*(5-alfa)*(6-alfa)*(7-alfa)*V(9,:)+...
            (1/factorial(10))*alfa*(1-alfa)*(2-alfa)*(3-alfa)*...
            (4-alfa)*(5-alfa)*(6-alfa)*(7-alfa)*(8-alfa)*V(10,:);
    otherwise
                disp('Incorrect order, must be between 1 and 10');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ù
 function y = Fobt(x)
 npar = length(x);
 z = 0;
  for i =2:npar
      z = z+x(i)^2;
  end
  y = x(1)^2+1000000 * z;
参考文献：(2020) "Fractional Lèvy flight bat algorithm for global optimisation", Int. J. Bio-Inspired Computation, Vol. 15, No. 2, pp.100–112