最適化アルゴリズムを実装していくぞ（コウモリアルゴリズム）



# はじめに

最適化アルゴリズムの実装シリーズです。
まずは[概要](https://qiita.com/pocokhc/items/07b698cc426cadb3a64e)を見てください。

コードは[github](https://github.com/pocokhc/metaheuristics)にあります。


# コウモリアルゴリズム

コウモリアルゴリズム(Bat Algorithm)はコウモリが超音波により獲物や障害物の位置を特定する反響定位行動に着目したアルゴリズムです。

コウモリは、速度、周波数、パルス率、音量の4つの要素を踏まえ行動を決めます。

・[進化計算アルゴリズム入門 生物の行動科学から導く最適解](https://www.amazon.co.jp/%E9%80%B2%E5%8C%96%E8%A8%88%E7%AE%97%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0%E5%85%A5%E9%96%80-%E7%94%9F%E7%89%A9%E3%81%AE%E8%A1%8C%E5%8B%95%E7%A7%91%E5%AD%A6%E3%81%8B%E3%82%89%E5%B0%8E%E3%81%8F%E6%9C%80%E9%81%A9%E8%A7%A3-%E5%A4%A7%E8%B0%B7%E7%B4%80%E5%AD%90-ebook/dp/B07DQFVK1H)



* アルゴリズムのフロー


![draw-Bat.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/396150/e258f932-e055-2cd2-77d1-f4f850270bce.png)

![draw-Bat2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/396150/ad9b515a-ab12-5d11-f661-c23379acb677.png)


* 用語の対応

|問題|コウモリアルゴリズム|
|---|---|
|入力値の配列|コウモリの位置|
|入力値|コウモリ|
|評価値|コウモリの位置における評価値|


* コード内で使う変数の意味

|変数名|意味|所感|
|---|---|---|
|problem|任意の問題のクラス（問題編を参照）||
|bat_max|コウモリの数||
|frequency_min|周波数の最低値|最良コウモリへの移動割合の下限値|
|frequency_max|周波数の最大値|最良コウモリへの移動割合の上限値|
|good_bat_rate|良いコウモリの割合||
|volume_init|音量の初期値|音量が大きいほど乱数性が強くなり広域を探索|
|volume_update_rate|音量の減少率|音量は時間が進むほど小さくなる|
|pulse_convergence_value|パルス率の収束値|パルス率は良いコウモリの近くへ移動するかの確率、収束値が1だと最終的に移動しなくなる|
|pulse_convergence_speed|パルス率の収束の速さ|0.5だと10回、0.2だと20回ぐらいで収束する|

## 最良コウモリに向かって近づく

最良コウモリへの移動は、最良コウモリの位置と自分の位置の差分と周波数を元に速度を更新し、その速度を使って移動します。

$$
\vec{v_i}(t+1) = \vec{v_i}(t) + f_{req} \bigl( \vec{g}(t) - \vec{x_i}(t) \bigr)
$$

$$
\vec{x_i}(t+1) = \vec{x_i}(t) + \vec{v_i}(t+1)
$$


$\vec{v_i}$ が速度、$\vec{x_i}$が位置、$\vec{g}$ が最良コウモリの位置になります。


また、周波数は下限($f_{min}$)～上限($f_{max}$)の間における乱数となります。

$$
f_{req} = f_{min} + (f_{max} - f_{min}) \times rand[0,1]
$$


``` python
import random
import numpy as np

# 周波数
frequency = frequency_min + (frequency_max - frequency_min) * random.random()

pos_best = np.asarray(best_bat.getArray())
pos = np.asarray(bat["bat"].getArray())

# 速度を計算
bat["v"] = bat["v"] + frequency * (pos_best - pos)

# 位置を更新し、新しい位置のコウモリを作成
new_bat1 = problem.create(pos + bat["v"])
```


## 良いコウモリの近傍に移動する

良いコウモリへの移動ですが、まずパルス率と乱数を比較し、乱数値以上なら移動をしません。

移動する場合は、まず良いコウモリを選びます。
良いコウモリは全コウモリの内、上位数%のコウモリの中からランダムに選択します。

そして、良いコウモリの近くに移動します。

$$
\vec{x_i}(t+1) = \vec{b}(t) + \bar{A} \vec{\epsilon}
$$

ここで $\bar{A}$ は全コウモリの平均音量、$\vec{\epsilon}$ は-1～1の乱数ベクトルになります。


``` python
import random
import numpy as np

new_bat2 = None

# 乱数で良いコウモリの近くに移動
if bat["pulse"] < random.random():

    # ソートする
    bats.sort(key=lambda x: x["bat"].getScore())

    # いいコウモリをランダムに選ぶ
    r = random.randint(1, good_bat_num)
    pos_good = np.asarray(bats[-r]["bat"].getArray())
    
    # 全コウモリの平均音量
    volume = np.average([x["volume"] for x in self.bats])
    
    rn = np.random.uniform(-1, 1, len(pos_good))  # -1,1の乱数
    new_pos2 = pos_good + volume * rn
    new_bat2 = problem.create(new_pos2)
```


## 更新

パルス率($r_i$)の更新式は以下です。

$$
r_i(t+1) = R_0(1 - \epsilon^{-\gamma t} )
$$

$R_0$ はパルス率の収束値(定数)、$\gamma$ はパルス率の収束の速さ(定数)です。

また、音量($A_i$)の更新式は以下です。

$$
A_i(t+1) = \alpha A_i(t)
$$

$\alpha$ は音量を減少させる割合(定数)です。

更新はアルゴリズムのフロー図通りに更新していきます。

``` python
import random

volume_update_rate = 0.9
pulse_convergence_value = 0.5
pulse_convergence_speed = 0.9

def _calcPulse(count):
    return pulse_convergence_value * (1-math.exp(-pulse_convergence_speed*count))

# ランダムに新しい位置を生成
new_bat3 = problem.create()

score1 = new_bat1.getScore()
score2 = None if new_bat2 is None else new_bat2.getScore()
score3 = new_bat3.getScore()

if (score2 is None or score1 > score2) and score1 > score3 :
    # 音量の乱数
    if random.random() < bat["volume"]:
        if score2 is None or score2 <= score3:
            # ランダムな位置で更新
            bat["bat"] = new_bat3
        else:
            # 良いコウモリ付近の位置で更新
            bat["bat"] = new_bat3

        # パルス率の更新
        bat["pulse"] = _calcPulse(step_count)

        # 音量の更新
        bat["volume"] = volume_update_rate * bat["volume"]

    else:
        # 最良コウモリへの位置で更新
        bat["bat"] = new_bat1
else:
    # 最良コウモリへの位置で更新
    bat["bat"] = new_bat1

# ステップカウント+1
step_count += 1
```

### パルス率について

パルス率ですが、収束スピードによって以下のように変化します。(pulse_convergence_value=1の場合)
0.2でも２０回ぐらいで収束と、結構収束するのが早いので小さい値がいいかもしれません。

![bat_pulse_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/396150/f14c9853-0429-4cec-e0ff-ef503f3228b8.png)



## コード全体


``` python
import math
import random

import numpy as np

class Bat():
    def __init__(self,
        bat_max,
        frequency_min=0,
        frequency_max=1,

        good_bat_rate=0.2,
        volume_init=1.0,
        volume_update_rate=0.9,
        pulse_convergence_value=0.5,
        pulse_convergence_speed=0.9,
    ):
        self.bat_max = bat_max
        self.frequency_min = frequency_min
        self.frequency_max = frequency_max

        self.good_bat_num = int(bat_max * good_bat_rate + 0.5)
        if self.good_bat_num > bat_max:
            self.good_bat_num = bat_max
        if self.good_bat_num < 1:
            self.good_bat_num = 1
        self.volume_init = volume_init
        self.volume_update_rate = volume_update_rate
        self.pulse_convergence_value = pulse_convergence_value
        self.pulse_convergence_speed = pulse_convergence_speed


    def init(self, problem):
        self.problem = problem
        self.count = 0

        self.step_count = 0
        self.best_bat = None
        self.bats = []
        for _ in range(self.bat_max):
            o = problem.create()

            d = {
                "bat": o,
                "pulse": self._calcPulse(0),
                "volume": self.volume_init,
                "v": np.zeros(problem.size)
            }
            self.bats.append(d)

            if self.best_bat is None or self.best_bat.getScore() < o.getScore():
                self.best_bat = o

    def _calcPulse(self, count):
        return self.pulse_convergence_value*(1-math.exp(-self.pulse_convergence_speed*count))
        
    def step(self):
        self.bats.sort(key=lambda x: x["bat"].getScore())

        for bat in self.bats:

            # 周波数
            frequency = self.frequency_min + (self.frequency_max - self.frequency_min) * random.random()
            
            # 最良のコウモリに近づく
            pos_best = self.best_bat.getArray()
            pos = bat["bat"].getArray()
            bat["v"] = bat["v"] + frequency * (pos_best - pos)
            new_bat1 = self.problem.create(pos + bat["v"])

            # 音量判定
            if random.random() >= bat["volume"]:
                # 新しい位置で更新
                bat["bat"] = new_bat1
            
                # 最良コウモリチェック
                if self.best_bat.getScore() < bat["bat"].getScore():
                    self.best_bat = bat["bat"]
                continue

            # 乱数で良いコウモリの近くに移動
            new_bat2 = None
            if bat["pulse"] < random.random():
            
                # 全コウモリの平均音量
                volume = np.average([x["volume"] for x in self.bats])
                
                # いいコウモリを選ぶ
                r = random.randint(1, self.good_bat_num)
                pos_good = self.bats[-r]["bat"].getArray()
                
                rn = np.random.uniform(-1, 1, len(pos_good))  # -1,1の乱数
                new_bat2 = self.problem.create(pos_good + volume * rn)
            
            # ランダムに生成
            new_bat3 = self.problem.create()
                
            # 新しい位置の比較
            score1 = new_bat1.getScore()
            score2 = None if new_bat2 is None else new_bat2.getScore()
            score3 = new_bat3.getScore()
            if (score2 is None or score1 > score2) and score1 > score3 :
                if score2 is None or score2 <= score3:
                    bat["bat"] = new_bat3
                else:
                    bat["bat"] = new_bat2
                
                # パルス率の更新
                bat["pulse"] = self._calcPulse(self.step_count)

                # 音量の更新
                bat["volume"] = self.volume_update_rate * bat["volume"]
            else:
                bat["bat"] = new_bat1

            # 最良コウモリチェック
            if self.best_bat.getScore() < bat["bat"].getScore():
                self.best_bat = bat["bat"]

        self.step_count += 1


```

# ハイパーパラメータ例

各問題に対して optuna でハイパーパラメータを最適化した結果です。
最適化の1回の試行は、探索時間を2秒間として結果を出しています。
これを100回実行し、最適なハイパーパラメータを optuna に探してもらいました。
※frequency_minは0で固定しています

|問題|bat_max|frequency_max|frequency_min|good_bat_rate|pulse_convergence_speed|pulse_convergence_value|volume_init|volume_update_rate|
|---|---|---|---|---|---|---|---|---|
|EightQueen|21|0.16056311709148877|0.0|0.18895324827898902|0.7743073337572011|0.015299647553258189|0.37862199762341764|0.5219340167573981|
|function_Ackley|6|0.5807409179883148|0.0|0.8431315762963505|0.6168807797723693|0.41730375451852453|0.41343060460565206|0.15049946819845028|
|function_Griewank|8|0.9452402185165402|0.0|0.7201234515965615|0.7260546767194784|0.130926854366477|0.013897440965911972|0.018889941226347226|
|function_Michalewicz|49|0.4179217659403093|0.0|0.6361175345153247|0.5018316263776195|0.42459903007426236|0.4968382725357193|0.1065841508592828|
|function_Rastrigin|4|0.6140740614258697|0.0|0.08321822782531468|0.37791620723418984|0.5942334579901081|0.2272362640804505|0.6090017290739149|
|function_Schwefel|5|0.9974629382587714|0.0|0.48057178703432885|0.017850539098788976|0.8624933090356682|0.008382588455720222|0.9676085790533966|
|function_StyblinskiTang|9|0.9919883606063762|0.0|0.3096861786862085|0.6695497098867679|0.8832509763124465|1.3802416508441109|0.24307207412221196|
|LifeGame|17|0.17520950919099773|0.0|0.5512117484751708|0.3239752160390367|0.5877290654832432|0.6901442886406529|0.7052284742486757|
|OneMax|47|0.4118670851332011|0.0|0.4420165018547974|0.1217770472809899|0.8862495030963856|1.0185236105175342|0.5188797824478796|
|TSP|8|0.44386885244368085|0.0|0.19840772635097428|0.24715859136402973|0.3995313116893353|0.40899114232584255|0.00012730831532623|


# 実際の動きの可視化

1次元は6個体、2次元は20個体で50step実行した結果です。
赤い丸がそのstepでの最高スコアを持っている個体となります。

パラメータは以下で実行しました。

``` python
Bat(N, frequency_min=0, frequency_max=0.05, good_bat_rate=0.1, volume_init=0.5, pulse_convergence_value=0.9, pulse_convergence_speed=0.1)  
```



## function_Ackley

* 1次元

![function_Ackley_Bat_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Ackley_Bat_2.gif?raw=true)


* 2次元

![function_Ackley_Bat_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Ackley_Bat_3.gif?raw=true)



## function_Rastrigin

* 1次元

![ffunction_Rastrigin_Bat_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Rastrigin_Bat_2.gif?raw=true)


* 2次元

![function_Rastrigin_Bat_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Rastrigin_Bat_3.gif?raw=true)




## function_Schwefel

* 1次元

![function_Schwefel_Bat_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Schwefel_Bat_2.gif?raw=true)


* 2次元

![function_Schwefel_Bat_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_Schwefel_Bat_3.gif?raw=true)




## function_StyblinskiTang

* 1次元

![function_StyblinskiTang_Bat_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_StyblinskiTang_Bat_2.gif?raw=true)


* 2次元

![function_StyblinskiTang_Bat_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_StyblinskiTang_Bat_3.gif?raw=true)



## function_XinSheYang

* 1次元

![function_XinSheYang_Bat_2.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_XinSheYang_Bat_2.gif?raw=true)


* 2次元

![function_XinSheYang_Bat_3.gif](https://github.com/pocokhc/metaheuristics/blob/main/imgs/function_XinSheYang_Bat_3.gif?raw=true)





# あとがき

近傍に移動した後に速度が残ったままになるので、それで広域を探索している感じがします。
また、速度が基本残ったままになるので後半ほど移動がすごいおおざっぱになるような…。
ハイパーパラメータの調整が結構シビアかもしれません。