最適化アルゴリズムを実装していくぞ（差分進化）





# はじめに

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

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



# 差分進化

差分進化（Differential evolution：DE）は、生物の進化をもとに考案された遺伝的アルゴリズムと似ている手法です。

大きく、Mutation(突然変異)、Crossover(交叉)、Selection(生存選択)の3つのフェーズがあります。
遺伝的アルゴリズムみたいにフェーズ毎にアルゴリズムがいろいろあるようですが、本記事では一番簡単なアルゴリズムで実装しています。


参考
・[差分進化法でハイパーパラメータチューニング](https://qiita.com/n-suzuki/items/b8d4ccc4b6936120567e)
・[差分進化(Wikipedia)](https://ja.wikipedia.org/wiki/%E5%B7%AE%E5%88%86%E9%80%B2%E5%8C%96)


* アルゴリズムのフロー

![draw2-DE.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/396150/882111cf-0140-3c77-f23f-a08102778409.png)



* 用語の対応

|問題|差分進化|
|---|---|
|入力値の配列|エージェントの位置|
|入力値|エージェント|
|評価値|エージェントの評価値|


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

|変数名|意味|所感|
|---|---|---|
|crossover_rate|交叉する確率|変異ベクトルになる確率|
|scaling|変異ベクトルのスケール率(長さ)|大きいほど移動範囲が大きい|


## Mutation(突然変異)

ランダムに選んだ3人のエージェントから変異ベクトルを作成します。

``` python
import random

i = 自分を表すindex

# i番目を除いた3エージェントをランダムに選択
r1, r2, r3 = random.sample([ j for j in range(len(agents)) if j != i ], 3)
pos1 = agents[r1].getArray()
pos2 = agents[r2].getArray()
pos3 = agents[r3].getArray()
```

変異ベクトルは以下の式で出します。

$$ \vec{v} = \vec{x_1} + F (\vec{x_2} - \vec{x_3}) $$

$F$ は変異差分の適用率(スケール因子:scaling factor)を表す0～2の実数です。
pythonコードだと以下になります。

``` python
# ベクトル計算しやすいようにnumpy化
import numpy as np
pos1 = np.asarray(pos1)
pos2 = np.asarray(pos2)
pos3 = np.asarray(pos3)

m_pos = pos1 + scaling * (pos2 - pos3)
```



## Crossover(交叉)

一様交叉(binary crossover)で交叉します。
ある確率(crossover_rate)で成分を入れ替える交叉となります。

また交叉は、変異ベクトル側の要素を1成分だけ必ず取り入れます。

``` python
import random

# 変異ベクトルで交叉させる(一様交叉)
pos = agent.getArray()
ri = random.randint(0, len(pos))  # 1成分は必ず変異ベクトル
for j in range(len(pos)):
    if  ri == j or random.random() < self.crossover_rate:
        pos[j] = m_pos[j]
    else:
        pass  # 更新しない

# 新しい位置のエージェントを作成
new_agent = problem.create(pos)
```


## 生存選択

交叉でできた新しいエージェントと今のエージェントを比べ、更新されていれば置き換えます。

``` python
# 優れている個体なら置き換える
if agents[i].getScore() < new_agent.getScore():
    agents[i] = new_agent
```



## コード全体

コード全体です。
クラス化しているので上記コードと少し違います。

``` python
import math
import random

import numpy as np

class DE():
    def __init__(self,
        agent_max,           # エージェント数
        crossover_rate=0.5,  # 交叉率
        scaling=0.5,         # 差分の適用率
    ):
        self.agent_max = agent_max
        self.crossover_rate = crossover_rate
        self.scaling = scaling

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

        # 初期位置の生成
        self.agents = []
        for _ in range(self.agent_max):
            self.agents.append(problem.create())


    def step(self):

        for i, agent in enumerate(self.agents):

            # iを含まない3個体をランダムに選択
            r1, r2, r3 = random.sample([ j for j in range(len(self.agents)) if j != i ], 3)
            pos1 = self.agents[r1].getArray()
            pos2 = self.agents[r2].getArray()
            pos3 = self.agents[r3].getArray()

            # 3個体から変異ベクトルをだす
            pos1 = np.asarray(pos1)
            pos2 = np.asarray(pos2)
            pos3 = np.asarray(pos3)
            m_pos = pos1 + self.scaling * (pos2 - pos3)

            # 変異ベクトルで交叉させる(一様交叉)
            pos = agent.getArray()
            ri = random.randint(0, len(pos))  # 1成分は必ず変異ベクトル
            for j in range(len(pos)):
                if  ri == j or random.random() < self.crossover_rate:
                    pos[j] = m_pos[j]
                else:
                    pass  # 更新しない

            # 優れている個体なら置き換える
            new_agent = self.problem.create(pos)
            self.count += 1
            if agent.getScore() < new_agent.getScore():
                self.agents[i] = new_agent


```




# ハイパーパラメータ例

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


|問題|agent_max|crossover_rate|scaling|
|---|---|---|---|
|EightQueen|5|0.008405098138137779|1.7482804860765253|
|function_Ackley|36|0.4076390525351224|0.2908895854800526|
|function_Griewank|14|0.27752386128521395|0.4629100940098222|
|function_Michalewicz|12|0.1532879607238835|0.0742830755371933|
|function_Rastrigin|28|0.33513859646880306|0.0754225020709786|
|function_Schwefel|13|0.00032331965923372563|0.13153649005308807|
|function_StyblinskiTang|39|0.21247741932099348|0.08732185323441227|
|function_XinSheYang|33|0.0955103914325307|0.008270969294347359|
|LifeGame|39|0.6612227467897149|1.136453380180552|
|OneMax|4|0.1190487045395953|1.1581036102901494|
|TSP|23|0.41212989299137665|0.014644735558753091|


# 実際の動きの可視化

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

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

``` python
DE(N, crossover_rate=0, scaling=0.4)
```



## function_Ackley

* 1次元

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


* 2次元

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





## function_Rastrigin

* 1次元

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


* 2次元

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




## function_Schwefel

* 1次元

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


* 2次元

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




## function_StyblinskiTang

* 1次元

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


* 2次元

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



## function_XinSheYang

* 1次元

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


* 2次元

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



# あとがき

いい感じに収束していきますね。
ただランダム移動がないので局所解に陥りやすい気はします。