Python和Matlab如何實現(xiàn)蝙蝠算法

這篇文章主要介紹“Python和Matlab如何實現(xiàn)蝙蝠算法”的相關(guān)知識，小編通過實際案例向大家展示操作過程，操作方法簡單快捷，實用性強，希望這篇“Python和Matlab如何實現(xiàn)蝙蝠算法”文章能幫助大家解決問題。

1 前言

蝙蝠算法是2010年楊教授基于群體智能提出的啟發(fā)式搜索算法，是一種搜索全局最優(yōu)解的有效方法。該算法基于迭代優(yōu)化，初始化為一組隨機解，然后迭代搜尋最優(yōu)解，且在最優(yōu)解周圍通過隨機飛行產(chǎn)生局部新解，加強局部搜索速度。該算法具有實現(xiàn)簡單、參數(shù)少等特點。

該算法主要用于目標(biāo)函數(shù)尋優(yōu)，基于蝙蝠種群利用產(chǎn)生的聲波搜索獵物和控制飛行方向的特征來實現(xiàn)函數(shù)的尋優(yōu)。以一只蝙蝠作為基本單元，且每只蝙蝠都有一個適應(yīng)值來對函數(shù)解空間進行優(yōu)化。每只蝙蝠可以調(diào)整自身發(fā)射聲波的響度、頻率等對空間進行搜索，使整個種群的活動逐步由無序變?yōu)橛行颉?strong>但蝙蝠算法在尋優(yōu)末段容易陷入局部的極值，本文引入變速度權(quán)重因子修正系數(shù)(和粒子群參數(shù)調(diào)整的自適應(yīng)過程相似)，盡可能避免局部極值的困境，從而達到全局尋優(yōu)的效果。

2 蝙蝠算法原理細(xì)講

首先在自變量范圍內(nèi)產(chǎn)生蝙蝠的隨機位置；

然后每個蝙蝠向周圍某處移動，飛行靠聲波反饋來更變方向，又向周圍移動是隨機的，蝙蝠下一時刻移動到某一位置有一定的出現(xiàn)頻率，所以本文在運動公式中加了聲波和頻率兩個因素。每一時刻的位移看作是一次飛行，每次飛行的距離很短（距離長短反映出搜素精度）。

每只蝙蝠初始位置是隨機的，在初始位置時刻其中一只蝙蝠對應(yīng)的函數(shù)最值，算法開始會記錄該位置，然后所有蝙蝠逐漸向該位置靠近，飛行方向大致向當(dāng)前最值方向（即該方向的位置該蝙蝠的出現(xiàn)頻率更高），但是在飛行方向也是隨機飛行的，相當(dāng)于逐步搜索過去。蝙蝠飛行過程中不會因為當(dāng)前飛行出現(xiàn)的最大值位置而改變種群的飛行趨勢，這樣可以盡量避免算法陷入局部極值。

算法每飛行一次就記錄一次種群所處位置的最值，最后找出記錄中的最值。

算法有兩個參數(shù)可以影響最終的結(jié)果：種群數(shù)量和飛行次數(shù)。其中種群數(shù)量的影響是最大的。 

3 詳細(xì)步驟

3.1 初始化相關(guān)參數(shù)

蝙蝠的位置為Xi，飛行速度Vi，聲音響度為Ai，頻率yi范圍，設(shè)有目標(biāo)函數(shù)為 ：

3.2 更改脈沖頻率產(chǎn)生的解并更變蝙蝠位置與飛行速度

蝙蝠i在t-1時的位置和飛行速度分別表示為和，群體當(dāng)前找到的最優(yōu)位置為。接著根據(jù)自身發(fā)出不同的音響搜尋獵物，通過接受反饋信息來調(diào)整位置xi和飛行速度v(i)。其飛行的速度更變公式如下：

 w(t)其中為時刻變速慣性權(quán)重因子，作用是使蝙蝠的前期搜索對后期搜索提供參照，wmax為w(t)的最大值、wmin為w(t)的最小值；，一般取2，Tmax為最大迭代次數(shù)；為當(dāng)前位置最優(yōu)解；y(i)為頻率滿足正態(tài)均勻分布的一個隨機數(shù)，β是一個隨機變量，且。開始運行時，蝙蝠在隨機進行頻率分配。

為控制蝙蝠所處位置在自變量范圍內(nèi)，本文針對該情況設(shè)置了邊界規(guī)則：如果下次運動的位置超出了自變量范圍，那么下次飛行的位置為投影在的邊界上的位置。

3.3 搜尋局部最優(yōu)解

3.4 通過蝙蝠多次飛行產(chǎn)生多個新解，進行全局搜索，若得到的新解

那么接受該解。

3.5 排列所有蝙蝠的位置，并找出當(dāng)前最優(yōu)值及對應(yīng)的位置

3.6 設(shè)當(dāng)前最優(yōu)解為，然后使所有蝙蝠繼續(xù)向下一時刻運動，并返回步驟2重新計算。

3.7 時刻結(jié)束，輸出:最優(yōu)解

4 Python實現(xiàn)

4.1 代碼

#=========導(dǎo)入相關(guān)庫===============
import numpy as np
from numpy.random import random as rand
 
#========參數(shù)設(shè)置==============
# objfun:目標(biāo)函數(shù) 
# N_pop: 種群規(guī)模，通常為10到40
# N_gen: 迭代數(shù)
# A: 響度（恒定或降低） 
# r: 脈沖率（恒定或減?。?nbsp;
# 此頻率范圍決定范圍
# 如有必要，應(yīng)更改這些值 
# Qmin: 頻率最小值
# Qmax: 頻率最大值
# d: 維度
# lower: 下界
# upper: 上界
 
def bat_algorithm(objfun, N_pop=20, N_gen=1000, A=0.5, r=0.5,
    Qmin=0, Qmax=2, d=10, lower=-2, upper=2):
 
    N_iter = 0 # Total number of function evaluations
 
    #=====速度上下限================
    Lower_bound = lower * np.ones((1,d))
    Upper_bound = upper * np.ones((1,d))
 
    Q = np.zeros((N_pop, 1)) # 頻率
    v = np.zeros((N_pop, d)) # 速度
    S = np.zeros((N_pop, d))
 
    #=====初始化種群、初始解=======
    # Sol = np.random.uniform(Lower_bound, Upper_bound, (N_pop, d))
    # Fitness = objfun(Sol)
    Sol = np.zeros((N_pop, d))
    Fitness = np.zeros((N_pop, 1))
    for i in range(N_pop):
        Sol[i] = np.random.uniform(Lower_bound, Upper_bound, (1, d))
        Fitness[i] = objfun(Sol[i])
 
    #====找出初始最優(yōu)解===========
    fmin = min(Fitness)
    Index = list(Fitness).index(fmin)
    best = Sol[Index]
 
    #======開始迭代=======
    for t in range(N_gen):
 
        #====對所有蝙蝠/解決方案進行循環(huán) ======
        for i in range(N_pop):
            # Q[i] = Qmin + (Qmin - Qmax) * np.random.rand
            Q[i] = np.random.uniform(Qmin, Qmax)
            v[i] = v[i] + (Sol[i] - best) * Q[i]
            S[i] = Sol[i] + v[i]
 
            #===應(yīng)用簡單的界限/限制====
            Sol[i] = simplebounds(Sol[i], Lower_bound, Upper_bound)
            # Pulse rate
            if rand() > r:
                # The factor 0.001 limits the step sizes of random walks
                S[i] = best + 0.001*np.random.randn(1, d)
 
            #====評估新的解決方案 ===========
            # print(i)
            Fnew = objfun(S[i])
            #====如果解決方案有所改進，或者聲音不太大，請更新====
            if (Fnew <= Fitness[i]) and (rand() < A):
                Sol[i] = S[i]
                Fitness[i] = Fnew
 
            #====更新當(dāng)前的最佳解決方案======
            if Fnew <= fmin:
                best = S[i]
                fmin = Fnew
 
        N_iter = N_iter + N_pop
 
    print('Number of evaluations: ', N_iter)
    print("Best = ", best, '\n fmin = ', fmin)
 
    return best
 
 
def simplebounds(s, Lower_bound, Upper_bound):
 
    Index = s > Lower_bound
    s = Index * s + ~Index * Lower_bound
    Index = s < Upper_bound
    s = Index * s + ~Index * Upper_bound
 
    return s
 
 
#====目標(biāo)函數(shù)=============
def test_function(u):
    a = u ** 2
    return a.sum(axis=0)
 
 
if __name__ == '__main__':
    # print(bat_algorithm(test_function))
    bat_algorithm(test_function)

4.2 結(jié)果

5 Matlab實現(xiàn)

5.1 代碼

clear
wmax=0.9;%慣性權(quán)重最大值
wmin=0.4;%慣性權(quán)重最小值
n=10000; % 群體大小
A=rand(1,n); % 聲音響度 (不變或者減小)
%% 頻率范圍
Qmin=0; % 最低頻率
Qmax=1; % 最高頻率
d=2;% 搜索變量的維數(shù)（即頻率和速度）
%% 初始矩陣
Q=zeros(n,1); % 頻率矩陣初始化
v=zeros(n,d); % 速度矩陣初始化，初始化意義就是產(chǎn)生一個初始矩陣
%% x自變量范圍
u=-3;
o=12.1;
% y自變量范圍
p=4.1;
l=5.8;
%% 初始化群體/解
for i=1:n
    Sol(i,1)=-3+(12.1+3)*rand(1,1);%x自變量范圍【-3，12.1】
    Sol(i,2)=4.1+(5.8-4.1)*rand(1,1);%y自變量【4.1，5.8范圍】
    %將隨機生成的兩個自變量帶入函數(shù)式
    Fitness(i)=Fun(Sol(i,:));%函數(shù)值
end
%% 尋找當(dāng)前最優(yōu)解
[fmax,I]=max(Fitness);
best=Sol(I,:);
T=100;%飛行次數(shù)
%% 開始飛行
for t=1:T
    for i=1:n,
        Q(i)=Qmin+(Qmin-Qmax)*rand;%rand均勻分布的隨機數(shù)
        %v(i,:)=v(i,:)+(Sol(i,:)-best)*Q(i);（原速度）
        w=(wmax-wmin)*exp(-2*(t/T)^2)+wmin;%慣性權(quán)重因子
        v(i,:)=w*v(i,:)+(Sol(i,:)-best)*A(i)*Q(i);%更改后的速度
        S(i,:)=Sol(i,:)+v(i,:);%位置移動
       %% 邊界問題，如果下次飛行超出自變量范圍外了，那么下次飛行的位置為投影在的邊界上的位置
        %x軸
        if S(i,1)>o
            S(i,1)=o;
        end
        if S(i,1)<u
            S(i,1)=u;
        end
        %y軸
        if S(i,2)>l
            S(i,2)=l;
        end
        if S(i,2)<p
            S(i,2)=p;
        end
        %% 評估該次飛行后產(chǎn)生的新解
        Fnew(i)=Fun(S(i,:));
    end
    [Fmax,Z]=max(Fnew);%找出該次飛行后產(chǎn)生的最大值
    C(t,:)=S(Z,:);
    FFnew(t)=Fmax;
end
[Ffmax,N]=max(FFnew);%找出整個飛行過程中的最大值
M=C(N,:)
Ffmax

%目標(biāo)函數(shù)
function z=Fun(u)
   z=21.5+u(1)*sin(4*pi*u(1))+u(2)*sin(20*pi*u(2));

5.2 結(jié)果 

5.3 展望

如果是其他函數(shù)怎么辦呢？

函數(shù)z=21.5+u(1)*sin(4*pi*u(1))+u(2)*sin(20*pi*u(2))中的兩個自變量對應(yīng)程序中的是Sol(i,1)=-3+(12.1+3)*rand(1,1);%x自變量范圍【-3，12.1】和Sol(i,2)=4.1+(5.8-4.1)*rand(1,1);%y自變量【4.1，5.8范圍】，兩自變量產(chǎn)生的是列矩陣，而程序中Sol(i,:) 提取的是矩陣中的行，所以如果對該函數(shù)增減自變量，或者想求其他含有多個自變量的函數(shù)，程序中只用修改以下程序部分，其他參數(shù)也可以自行更改（注：自變量的范圍和個數(shù)增加了，那么種群個數(shù)和飛行次數(shù)務(wù)必要增加）：

Sol(i,1)=-3+(12.1+3)*rand(1,1);%x自變量范圍【-3，12.1】

Sol(i,2)=4.1+(5.8-4.1)*rand(1,1);%y自變量【4.1，5.8范圍】

和

% x自變量范圍

u=-3;

o=12.1;

% y自變量范圍

p=4.1;

l=5.8;

關(guān)于“Python和Matlab如何實現(xiàn)蝙蝠算法”的內(nèi)容就介紹到這里了，感謝大家的閱讀。如果想了解更多行業(yè)相關(guān)的知識，可以關(guān)注億速云行業(yè)資訊頻道，小編每天都會為大家更新不同的知識點。