Matlab实现黑洞优化算法的示例代码

clear all
clc
close all
%% 初始化参数
d=5;                % 维度
options.lk=-32*ones(1,d);   % 下限
options.uk=32*ones(1,d);    % 上限
options.m=50; % 种群数量
options.MAXITER=500; % 最大迭代次数
options.n=length(options.uk);    % 问题的维度。
options.ObjFunction=@Ackley; % 调用目标函数
options.Display_Flag=1; % 在迭代中显示结果的标志
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('迭代次数')
ylabel('适应度值')
 
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)
%% 目标函数值
if (nargin==0)
    F=[];
    d=2;                % 维度
    lb=-32*ones(1,d);   % 下限
    ub=32*ones(1,d);    % 上限
    FGO=0;              % 总体最优值
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)
%--------------------------------------------------------------------------
% 黑洞优化算法
%--------------------------------------------------------------------------
% 在搜索空间中初始化具有随机位置的恒星群。
% 开始。
%   对于每一颗星，评估目标函数。
%   选择具有最佳适应值的最佳恒星作为黑洞。
%   根据公式更改每颗恒星的位置。
%   如果一颗恒星到达一个比黑洞成本更低的位置，交换它们的位置。
%   如果一颗恒星穿过黑洞的视界，在搜索空间中的随机位置用一颗新恒星替换它。
%   如果满足终止标准(最大迭代次数或足够好的适应度)，则退出循环。
% 结束循环
%--------------------------------------------------------------------------
 
ObjFunction=options.ObjFunction; % 目标函数的命名
n=options.n;    % 该问题的维度
uk=options.uk;   % 在第k维上界。
lk=options.lk;  % 在第k的下界。
m=options.m; % m：样本点数
MAXITER=options.MAXITER; % 最大迭代次数
nEval=0;
[x,xBH,iBH,ObjFunctionValue]=Initialize(options);
nEval=nEval+size(x,1);
for iteration =1:MAXITER
    %%  根据公式改变每颗星的位置。
    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('迭代次数N为 %g 最优适应度为 %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