clear all
clc
close all
%% 初始化參數(shù)
d=5;                % 維度
options.lk=-32*ones(1,d);   % 下限
options.uk=32*ones(1,d);    % 上限
options.m=50; % 種群數(shù)量
options.MAXITER=500; % 最大迭代次數(shù)
options.n=length(options.uk);    % 問題的維度。
options.ObjFunction=@Ackley; % 調(diào)用目標(biāo)函數(shù)
options.Display_Flag=1; % 在迭代中顯示結(jié)果的標(biāo)志
options.run_parallel_index=0;
options.run=10;
 
if options.run_parallel_index
    stream = RandStream('mrg32k3a');
    parfor index=1:options.run
        set(stream,'Substream',index);
        RandStream.setGlobalStream(stream)
        [bestX, bestFitness, bestFitnessEvolution,nEval]=BH_v1(options);
        bestX_M(index,:)=bestX;
        Fbest_M(index)=bestFitness;
        fbest_evolution_M(index,:)=bestFitnessEvolution;
    end
else
    rng('default')
    for index=1:options.run
        [bestX, bestFitness, bestFitnessEvolution,nEval]=BH_v1(options);
        bestX_M(index,:)=bestX;
        Fbest_M(index)=bestFitness;
        fbest_evolution_M(index,:)=bestFitnessEvolution;
    end
end
 
%% 可視化
[a,b]=min(Fbest_M);
figure
plot(1:options.MAXITER,fbest_evolution_M(b,:))
xlabel('迭代次數(shù)')
ylabel('適應(yīng)度值')
 
fprintf(' MIN=%g  MEAN=%g  MEDIAN=%g MAX=%g  SD=%g \n',...
    min(Fbest_M),mean(Fbest_M),median(Fbest_M),max(Fbest_M),std(Fbest_M))
 

function [F, lb, ub, FGO] = Ackley(x)
%% 目標(biāo)函數(shù)值
if (nargin==0)
    F=[];
    d=2;                % 維度
    lb=-32*ones(1,d);   % 下限
    ub=32*ones(1,d);    % 上限
    FGO=0;              % 總體最優(yōu)值
else    
    n=size(x,2);
    for ix=1:size(x,1)
        x0=x(ix,:);
        F(ix) = -20*exp(-0.2*sqrt(1/n*sum(x0.^2)))-...
                    exp(1/n*sum(cos(2*pi*x0)))+20+exp(1);
    end
end

function [bestX, bestFitness, bestFitnessEvolution,nEval]=BH_v1(options)
%--------------------------------------------------------------------------
% 黑洞優(yōu)化算法
%--------------------------------------------------------------------------
% 在搜索空間中初始化具有隨機位置的恒星群。
% 開始。
%   對于每一顆星，評估目標(biāo)函數(shù)。
%   選擇具有最佳適應(yīng)值的最佳恒星作為黑洞。
%   根據(jù)公式更改每顆恒星的位置。
%   如果一顆恒星到達一個比黑洞成本更低的位置，交換它們的位置。
%   如果一顆恒星穿過黑洞的視界，在搜索空間中的隨機位置用一顆新恒星替換它。
%   如果滿足終止標(biāo)準(zhǔn)(最大迭代次數(shù)或足夠好的適應(yīng)度)，則退出循環(huán)。
% 結(jié)束循環(huán)
%--------------------------------------------------------------------------
 
ObjFunction=options.ObjFunction; % 目標(biāo)函數(shù)的命名
n=options.n;    % 該問題的維度
uk=options.uk;   % 在第k維上界。
lk=options.lk;  % 在第k的下界。
m=options.m; % m：樣本點數(shù)
MAXITER=options.MAXITER; % 最大迭代次數(shù)
nEval=0;
[x,xBH,iBH,ObjFunctionValue]=Initialize(options);
nEval=nEval+size(x,1);
for iteration =1:MAXITER
    %%  根據(jù)公式改變每顆星的位置。
    for i = 1 : m
        if i ~= iBH
            landa=rand;
            for k = 1 : n
                if landa<0.5
                    x(i,k)=x(i,k) + rand*(xBH(k)- x(i,k));
                else
                    x(i,k)=x(i,k) + rand*(xBH(k)- x(i,k));
                end
            end
        end
    end
    %% 如果恒星到達一個比黑洞成本低的位置，交換它們的位置
    ObjFunctionValue=feval(ObjFunction,x);
    nEval=nEval+size(x,1);
    %如果一顆恒星越過黑洞的事件視界，則在搜索空間中隨機位置用一顆新恒星代替
    R=ObjFunctionValue(iBH)/sum(ObjFunctionValue);
    %     R=exp(-n*ObjFunctionValue(iBH)/sum(ObjFunctionValue))
    for i = 1 : m
        Distance(i)=norm(xBH- x(i,:));
    end
    [x,ObjFunctionValue]=NewStarGeneration(x,Distance,R,options,iBH,ObjFunctionValue);
    [x]=bound(x,lk,uk);
    [xBH,iBH]=argmin(x,ObjFunctionValue,options);
    
    %--------------------------------------------------------------------------------
    bestFitnessEvolution(iteration)=ObjFunctionValue(iBH);
    %--------------------------------------------------------------------------------
    
    
    if options.Display_Flag==1
        fprintf('迭代次數(shù)N為 %g 最優(yōu)適應(yīng)度為 %g\n',iteration,ObjFunctionValue(iBH))
    end
    
end
bestX=xBH;
bestFitness=ObjFunctionValue(iBH);
end
 
function [x,xBH,iBH,ObjFunctionValue]=Initialize(options)
ObjFunction=options.ObjFunction; % the name of the objective function.
n=options.n;    % n: dimension of the problem.
uk=options.uk;  % up: upper bound in the kth dimension.
lk=options.lk;  % lp: lower bound in the kth dimension.
m=options.m;    % m: number of sample points
 
for i = 1 : m
    for k = 1 : n
        landa=rand;
        x(i,k) = lk(k) + landa*(uk(k) - lk(k));
    end
end
% x(end+1,:)=x0;
ObjFunctionValue=feval(ObjFunction,x);
[index1,index2]=sort(ObjFunctionValue);
x=x(index2(1:m),:);
xBH=x(1,:);
iBH=1;
ObjFunctionValue=ObjFunctionValue(index2(1:m));
end
 
function [xb,ib,xw,iw]=argmin(x,f,options)
[minf,ib]=min(f);
xb=x(ib,:);
[maxf,iw]=max(f);
xw=x(iw,:);
end
 
 
function [x,ObjFunctionValue]=NewStarGeneration(x,Distance,R,options,iBH,ObjFunctionValue)
ObjFunction=options.ObjFunction;
n=options.n;    
uk=options.uk;  
lk=options.lk;  
index=find(Distance<R);
for i=1:length(index)
    if index(i) ~= iBH
        for k = 1 : n
            x(i,k) = lk(k) + rand*(uk(k) - lk(k));
        end
        ObjFunctionValue(i)=feval(ObjFunction,x(i,:));
    end
end
end
function [x]=bound(x,l,u)
for j = 1:size(x,1)
    for k = 1:size(x,2)
        % 如果上限越界
        if x(j,k) > u(k),
            x(j,k) = u(k);
        end
        % 如果下限越界
        if x(j,k) < l(k),
            x(j,k) = l(k);
        end
    end
end
end