最適化アルゴリズムを実装していくぞ（くじらさんアルゴリズム）



# はじめに

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

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


# くじらさんアルゴリズム
## 概要

くじらさんアルゴリズム(The Whale Optimization Algorithm; WOA)はクジラの捕食行動に着目して作られたアルゴリズムです。
クジラはバブルネットという群行動によって獲物を捕食します。
バブルネットは図のように円を描きながら獲物を追い詰めて捕食する行動のようです。

![Bubble-net-whale-hunting-behavior-23.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/396150/e509a992-bb42-a8be-e62d-178bf3b70b45.png)


参考
・[くじらさんアルゴリズム](https://qiita.com/shinjikato/items/afd754a4dfca7efc4cd1)
・[Whale Optimization Algorithm](https://en.wikiversity.org/wiki/Whale_Optimization_Algorithm)


## アルゴリズム

獲物に近づく、獲物を探す、円を描くの３つの行動のどれかを選んで獲物を見つけます。


* アルゴリズムのフロー

![draw2-WOA.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/396150/5bc2d524-73b5-fdd4-1ea1-bd5e73d5813e.png)


* 用語の対応

|問題|くじらさんアルゴリズム|
|---|---|
|問題への入力|クジラの位置≒獲物の位置|
|評価値|クジラの位置における評価値|



* ハイパーパラメータに関して

|変数名|意味|備考|
|---|---|---|
|a_decrease|aの減少値|最初はランダムに移動し、後半は獲物に近づく。aはその移行速度|
|logarithmic_spiral|対数螺旋の係数|大きいほど回るときに大きく動く|




## 獲物に近づく

獲物に近づくと書いていますが実際に獲物がどこにいるかは分からないので、現状で最も評価が高いクジラを獲物と扱い、そこに近づきます。

まずは $\vec{A}$ と $\vec{C}$ を出します。

$$ \vec{A} = 2a\vec{r_1} - a$$

$$ \vec{C} = 2\vec{r_2}$$

$\vec{r_1}$ と $\vec{r_2}$ は座標と同じ次元の 0～1 の乱数です。
$a$ は初期値2から毎step減少していく変数です。($2 \geqq a \geqq 0$)

ここで $\vec{A}$ のノルム（$|\vec{A}|$）が1以下なら獲物に近づく行動をし、
違う場合は獲物を探す行動をします。

獲物に近づく場合の数式は以下です。

$$ \vec{D} = |\vec{C} \vec{X_{best}}(t) - \vec{X}(t)| $$
$$ \vec{X}(t+1) = \vec{X_{best}}(t) - \vec{A}\vec{D} $$

$\vec{X}$ が自分の座標を表し、$\vec{X_{best}}$ が最も評価の高いクジラです。
コードだと以下になります。

``` python
import numpy as np

a = 2  # 毎step減少する

r1 = np.random.rand(problem.size)  # 0-1の乱数を生成
r2 = np.random.rand(problem.size)  # 0-1の乱数を生成

A = (2.0 * np.multiply(a, r1)) - a
C = 2.0 * r2

i = 対象のクジラのindex

if np.linalg.norm(A) < 1:
    pos = np.asarray(whales[i].getArray())  # 対象クジラの座標
    best_pos = np.asarray(best_whale.getArray())  # 最良クジラの座標

    D = np.linalg.norm(np.multiply(C, best_pos) - pos)
    pos = best_pos - np.multiply(A, D)
```

## 獲物を探す

獲物に近づくで説明した $|A|$ が1以上の場合は獲物に近づきます。
近づき方は"獲物に近づく"場合と同じで、対象のクジラが最も評価の高いクジラからランダムなクジラに変わっただけです。

$$ \vec{D} = |\vec{C} \vec{X_p}(t) - \vec{X}(t)| $$
$$ \vec{X}(t+1) = \vec{X_p}(t) - \vec{A}\vec{D} $$

$\vec{X_p}$ はランダムなクジラの座標となります。


``` python
import numpy as np

i = 対象のクジラのindex

if np.linalg.norm(A) >= 1:
    pos = np.asarray(whales[i].getArray())  # 対象クジラの座標
    p_pos = np.asarray(whales[random.randint(0, len(whales)-1)].getArray())  # ランダムなクジラの座標

    D = np.linalg.norm(np.multiply(C, p_pos) - pos)
    pos = p_pos - np.multiply(A, D)
```


## 円を描く

円は獲物（評価値最大のクジラ）を基準に以下の式で描きます。

$$ \vec{D'} = |\vec{X_{best}}(t) - \vec{X}(t)| $$
$$ \vec{X}(t+1) = \vec{D'} e^{bL} \cos(2 \pi L) + \vec{X_{best}}(t) $$

$b$ は対数螺旋の係数(はlogarithmic_spiral)、$L$ は[-1, 1]のランダムな値です。


``` python
import numpy as np

i = 対象のクジラのindex

pos = np.asarray(whales[i].getArray())  # 対象クジラの座標
best_pos = np.asarray(best_whale.getArray())  # 最良クジラの座標

# 回る
D = np.linalg.norm(best_pos - pos)
L = np.random.uniform(-1, 1, problem.size)  # [-1,1]の乱数
_b = logarithmic_spiral
pos = np.multiply(np.multiply(D, np.exp(_b*r)), np.cos(2.0*np.pi*r)) + best_pos

```


## スパイラルのグラフ

円を描くがどういう移動をしているかを見てみました。

グラフのx軸が-1から1の乱数($L$）で、y軸が移動距離です。
数式だと以下の部分だけをグラフ化しています。

$$ spiral = \vec{D'} e^{bL} \cos(2 \pi L) $$

* ユークリッド距離による変化

$b$(logarithmic_spiral) を1に固定して、$\vec{D'}$のみを変化させた場合のグラフです。
$\vec{D'}$ は自分と獲物(最良クジラ)とのユークリッド距離を表しているので、獲物との距離による移動の変化を表しているといえます。

![plot_WOA_spiral1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/396150/a5304af9-2d4e-a2eb-eea4-652a2f1a9ad9.png)


離れているほど大きく移動する可能性が高い事がわかりますね。


* 対数螺旋係数による変化

今度は$b$(logarithmic_spiral)を変化させた場合です。
$D'$ は 1 に固定しています。

![plot_WOA_spiral2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/396150/bd9da0e3-2655-b8a5-fc58-21476f34cf9d.png)


$D'$ の場合と同じように $b$ が大きくなれば変化が大きいようです。



## コード全体

コード全体です。
コードを作成するにあたり[ここのコード](https://github.com/docwza/woa/blob/master/src/whale_optimization.py)を参考にしています。

``` python
import math
import random

import numpy as np

class WOA():
    def __init__(self, 
        whale_max,         # クジラの数
        a_decrease=0.001,  # 変数aの減少値
        logarithmic_spiral=1,  # 対数螺旋の係数
    ):
        self.whale_max = whale_max
        self.a_decrease = a_decrease
        self.logarithmic_spiral = logarithmic_spiral

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

        self.best_whale = None
        self.whales = []
        for _ in range(self.whale_max):
            o = problem.create()
            self.whales.append(o)

            if self.best_whale is None or self.best_whale.getScore() < o.getScore():
                self.best_whale = o.copy()
        self._a = 2


    def step(self):

        for whale in self.whales:
            pos = np.asarray(whale.getArray())

            if random.random() < 0.5:
                r1 = np.random.rand(self.problem.size)  # 0-1の乱数
                r2 = np.random.rand(self.problem.size)  # 0-1の乱数

                A = (2.0 * np.multiply(self._a, r1)) - self._a
                C = 2.0 * r2

                if np.linalg.norm(A) < 1:
                    # 獲物に近づく
                    new_pos = self.best_whale.getArray()
                else:
                    # 獲物を探す
                    new_pos = self.whales[random.randint(0, len(self.whales)-1)].getArray()
                new_pos = np.asarray(new_pos)

                D = np.linalg.norm(np.multiply(C, new_pos) - pos)
                pos = new_pos - np.multiply(A, D)

            else:
                # 回る
                best_pos = np.asarray(self.best_whale.getArray())

                D = np.linalg.norm(best_pos - pos)
                L = np.random.uniform(-1, 1, self.problem.size)  # [-1,1]の乱数
                
                _b = self.logarithmic_spiral
                pos = np.multiply(np.multiply(D, np.exp(_b*L)), np.cos(2.0*np.pi*L)) + best_pos

            whale.setArray(pos)
            if self.best_whale.getScore() < whale.getScore():
                self.best_whale = whale.copy()

        self._a -= self.a_decrease
        if self._a < 0:
            self._a = 0

```

# ハイパーパラメータ例

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


|問題|a_decrease|logarithmic_spiral|whale_max|
|---|---|---|---|
|EightQueen|0.020366541568838378|1.5723091675627572|45|
|function_Ackley|0.012846425355295324|1.817054209171554|49|
|function_Griewank|0.014055305155620233|1.5507180701447845|48|
|function_Michalewicz|0.019651164901908023|1.7671279692872293|46|
|function_Rastrigin|0.0173318428180831|0.7578390758302801|46|
|function_Schwefel|0.007491644624065206|1.6917050129680944|9|
|function_StyblinskiTang|0.024426401837372255|1.9474471085818506|50|
|function_XinSheYang|0.005722874915111857|0.4779624408790928|24|
|g2048|0.043666294831303784|1.09050483219953|37|
|LifeGame|0.0054058667884643585|0.06834119701477159|50|
|OneMax|0.014922301994610046|0.8650190784964947|27|
|TSP|0.0287871077043457|1.3855183848189636|45|



# 実際の動きの可視化

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

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

``` python
WOA(N, a_decrease=2/50, logarithmic_spiral=1)
```


## function_Ackley

* 1次元

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


* 2次元

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





## function_Rastrigin

* 1次元

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


* 2次元

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




## function_Schwefel

* 1次元

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


* 2次元

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




## function_StyblinskiTang

* 1次元

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


* 2次元

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



## function_XinSheYang

* 1次元

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


* 2次元

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




# あとがき

ぐるぐる渦巻きながらいい感じに集まっていきますね。
ただ、一度収束してしまうとそこから脱出することはないので a の値をどう制御するかが重要そうです。