Matlab實現黑洞優化算法的示例代碼

前言

應用的領域很多。

1.概述

1.1黑洞算法

根據黑洞現象原理首次提出BH 算法,它在傳統PSO基礎上引入瞭新的機制,有效地提高瞭收斂速度並防止瞭陷入局部極值的情況發生;但是該方法卻沒有提及如何確定黑洞邊界和如何處理吸收星體的問題. Hatamlou BH算法進行瞭完善,讓其更加接近於黑洞的自然現狀,使其具有黑洞現象的一般特質.該方法不但沿襲瞭原BH算法的一切特性,同時還解決瞭原BH 算法的黑洞邊界問題和吸收多餘星點的問題.

1.2黑洞搜索優化算法

BH 算法的思想較為完整地描述瞭自然界中黑洞現象的一般特性;以BH 算法具有的這些特性為基礎,BH算法的尋優搜索原理比較簡單,易於實現.BH算法主要是模擬實際黑洞現象,在一定搜索空間內隨機佈置一定數量的星體,並通過統計學手段確定和評估搜索空間內每一個星體的適應度函數,選出一個具有最佳適應度值的星體作為黑洞.這個黑洞邊界(以R為半徑)被視為當前全局最優解所在區域,而黑洞本身被視為當前全局最優解;BH 算法中的黑洞具有與自然界黑洞同樣的強吸引能力,在搜索域內所有其他星體都將向其靠攏且無法逃逸.星體被黑洞吸引並向黑洞靠攏的公式如下:

式中:xi(t)和xi+1(t+1)分別代表第i個星體在t和 t+1次迭代時的位置;rand表示0,1]之間的隨機數;Am是搜索空間內的黑洞位置;N是星體的數量.

在更替迭代過程中黑洞將進入其邊界以內的星體一一吸收.每當一個星體被吸收的同時,在搜索空間內隨機產生一個新的星體﹐星體的總數保持不變.黑洞邊界范圍計算公式如下:

式中:fBH是黑洞的適應度函數值;fi是第i個星體的適應度函數值.

在被黑洞吸引的過程中,如果有一個星體i的適應度函數值fi優於黑洞的適應度函數值fBH則說明該黑洞所處位置並非最優,需要黑洞和星體互換位置.然後算法將以新生成的黑洞為中心繼續運行下去,其他星體將通過位置更替公式(1)重新湧向新的黑洞位置,並被其吸收.

1.3黑洞搜索算法的實現過程

根據上述分析,得出黑洞算法的實現步驟:

步驟1初始化參數,在搜索空間內隨機選取一定數量星體;

步驟2計算每個星體的適應度函數值;

步驟3選擇一個具有最佳適應度值的星體作為黑洞;

步驟4根據式(1)更改每個星體的位置;

步驟5如果星體最優適應度函數值優於黑洞的適應度函數值,則交換二者位置;

步驟6 如果一個星體在黑洞的邊界范圍內,那麼這個星體將被吸收﹐黑洞的邊界半徑由式(2)計算,同時,在搜索空間內隨機產生一個新的星體;

步驟7當系統達到最大迭代次數,或者出現一個最好的適應度值時,程序停止運行,否則返回步驟4.

2.Matlab代碼實現

2.1主函數

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))
 

2.2目標函數

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

2.3黑洞優化算法 

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

3.結果展現

以上就是Matlab實現黑洞優化算法的示例代碼的詳細內容,更多關於Matlab黑洞優化算法的資料請關註WalkonNet其它相關文章!

推薦閱讀: