ニューラルネットの重みを実数値遺伝的アルゴリズムで最適化してみる ～その2：倒立振子～

#はじめに
**ニューラルネット**の重みを**実数値遺伝的アルゴリズム**(以下、実数値GA)も用いて最適化するという手法(**NeuroEvolution**)を**強化学習**に適用してみました。

強化学習の対象は、古典的なベンチマーク問題である[倒立振子](https://ja.wikipedia.org/wiki/%E5%80%92%E7%AB%8B%E6%8C%AF%E5%AD%90)を選択しました。台車駆動型倒立振子といわれている方です。
はじめは自分でシミュレート環境を作成していましたが、OpenAI Gymにも倒立振子の環境があるようなのでそちらを用いました。

#OpenAI Gymとは
イーロン・マスクらが参加しているOpen AIが提供している強化学習アルゴリズムの検証プラットフォームです。
[https://gym.openai.com/](https://gym.openai.com/)

各ユーザのスコアは以前はAPIでOpenAIのサーバに直接送信して掲載していたようですが、現在(2019/7)ではgithubのwikiを編集して入力するようになっているようです。
[https://github.com/openai/gym/wiki/Leaderboard](https://github.com/openai/gym/wiki/Leaderboard)

##OpenAI Gymのインストール
Windows Anaconda環境で、「Anaconda Prompt」から以下のコマンドでインストールするだけで、倒立振子は動作することまで確認できました。アタリのゲームや3Dを使う場合は他のライブラリのインストールが必要かもしれません。

``` 
pip install gym
```

Google Colaboratoryでは、初期環境でインストール済みでした。

##倒立振子
OpenAI Gymの倒立振子のソースは以下にあります。
[https://github.com/openai/gym/blob/master/gym/envs/classic_control/cartpole.py](https://github.com/openai/gym/blob/master/gym/envs/classic_control/cartpole.py)

###レギュレーション
ソースのコメントとソースから読み取れる主なレギュレーションは以下の通りです。

- 状態(State)

| No       |  状態名     |  最小値   |最大値     |
|:---------|:------------|:---------|:---------|
| 0        | カートの位置(x) | -4.8     |4.8      |
| 1        | カートの速度(x') | -∞       |∞       |
| 2        | ポールの角度(θ) |-24deg    | 24deg   |
| 3        | ポールの速度(θ') | -∞       |  ∞       |

- 行動(Action)
    - 0 ： 左移動 (左方向に10Nの力を加える)
    - 1 ： 右移動 (右方向に10Nの力を加える)
- 状態初期値 ： 各状態で[-0.05～0.05]でランダムで。各状態の最小最大の5%ではないようです。なので、カートの位置では誤差のレベル、ポールの角度では最悪値で0.05[rad]→2.8[deg]になります。
- 終了条件
    - ポールの角度が-12～12deg外
    - カートの位置が-2.4～2.4外
    - 200step到達(1stepは0.02秒なので、4秒)
- 報酬(Reward)：1step毎に報酬+1 。なので、最大200の報酬が得られることになります。


#学習モデル
##ニューラルネットの構成
3層のニューラルネットを用いました。
- 入力層 ：4ニューロン
- 隠れ層 ：10ニューロン (活性化関数：シグモイド関数)
- 出力層 ：2ニューロン (活性化関数：Softmax関数)

![ニューラルネット_倒立振子.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/438288/6f2d6a7f-4b18-3130-2f7e-c76f9ae13a54.jpeg)

重み・バイアスの範囲：-5.0～5.0

##実装
``` nn.py
import numpy as np

# シグモイド 関数
def sigmoid(x):
    sigmoid_range = 34.538776394910684
    #オーバーフロー対策
    if x.any() <= -sigmoid_range:
        return 1e-15
    elif x.any() >= sigmoid_range:
        return 1.0 - 1e-15
    else:
        return 1.0 / (1.0 + np.exp(-x) )


# ソフトマックス 関数
def softmax(x):
    xmax = np.max(x) 
    exp_x = np.exp(x - xmax) 
    y = exp_x / np.sum(exp_x)
    return y

# Neural Network
class NN():

    #int    layer       : 層の数
    #int    input_size  : 入力層のニューロン数
    #int[]  hidden_size : 隠れ層のニューロン数
    #int    output_size : 出力層のニューロン数
    def __init__(self ,layer , input_size , hidden_size , output_size ):
        assert (layer >= 2 and input_size > 0 and output_size > 0) , 'ERROR:NN.init()'
        assert (layer-2 == len(hidden_size)) , 'ERROR:NN.init() hidden_size size'
        
        self.layer = layer
        self.input_unit  = input_size
        self.hidden_unit = hidden_size
        self.output_unit = output_size
            
        self.params = {}
        self.params['W1'] = np.zeros((input_size , hidden_size[0]) , dtype = 'float32')
        self.params['b1'] = np.zeros( hidden_size[0] , dtype = 'float32')
        self.params['W2'] = np.zeros((hidden_size[0] ,output_size), dtype = 'float32')
        self.params['b2'] = np.zeros( output_size , dtype = 'float32')
        
            
    #重み、バイアスの合計
    def get_weights_size(self) :
        w = 0
        for l in range(self.layer - 1):
            w_name = 'W' + str(l + 1)
            b_name = 'b' + str(l + 1)
            w += self.params[w_name].size
            w += self.params[b_name].size
            
        return w

    #重みの セット
    #numpy.ndarray weights : 重みとバイアス
    def set_weights(self , weights):
        w = 0
        
        for l in range(self.layer - 1):
            w_name = 'W' + str(l + 1)
            b_name = 'b' + str(l + 1)
            weight = weights[w: w+self.params[w_name].size]
            self.params[w_name] = weight.reshape(self.params[w_name].shape)
            w += self.params[w_name].size
            
            bias = weights[w: w+self.params[b_name].size]
            self.params[b_name] = bias.reshape(self.params[b_name].shape)
            w += self.params[b_name].size

        assert(weights.size == w) , \
            'Error:set_weights() weight size is invalid  weights = %d  nn = %d' %(weights.size , w)
        
    # 推論 

    def predict(self , x):
        W1 = self.params['W1']
        b1 = self.params['b1']
        W2 = self.params['W2']
        b2 = self.params['b2']
        
        i1 = np.dot(x,W1) + b1
        o1 = sigmoid(i1)
        i2 = np.dot(o1,W2) + b2
        y  = softmax(i2)
        
        return y
```

##実数値GAの構成
- 世代交代モデル：MGG
    - 世代交代数：2個体(エリート選択：1個体、ルーレット選択：1個体)
    - 致死個体(任意のパラメータが探索範囲を超えた個体) への対応：評価対象から外す
- 交叉法：UNDX (α=0.5、β=0.35)
- 学習パラメータ
    - 最大世代数：1000世代
    - 集団数：50個体
    - 生成子個体数: 100個体/世代
- 適合度：最小化問題にするために、OpenAI Gymの報酬をマイナスにしています。

##実装

``` nnga.py
import numpy as np
#OpenAIGym
import gym

from nn import NN

#重み の範囲
WEIGHT_MIN  = -5.0
WEIGHT_MAX  =  5.0

#評価値小 ＝ 優良個体
EVAL_TYPE_MIN = 0
#評価値大 ＝ 優良個体
EVAL_TYPE_MAX = 1

#層の数
LAYER_SIZE = 3
#入力層のニューロン数
INPUT_SIZE = 4
#出力層のニューロン数
OUTPUT_SIZE = 2

HIDDEN_SIZE = [10]

# 20step
SIMULATE_STEP = 200


#角度[rad]を-π～πに変換
def rad_conversion(rad):
    rad += np.pi;
    rad %= 2*np.pi;
    if rad < 0:
        rad += np.pi
    else:
        rad -= np.pi
    return rad

  
# Neural Network
class NNGA():
    #探索範囲取得
    def get_range() :
        ret = [WEIGHT_MIN , WEIGHT_MAX]
        return ret
    
    #評価タイプ取得
    def get_eval_type() :
        ret = EVAL_TYPE_MIN
        return ret
    
    def __init__(self ):
        self.nn = NN(LAYER_SIZE , INPUT_SIZE ,  HIDDEN_SIZE , OUTPUT_SIZE )
        
        print("Use Open AI Gym CartPole")
        self.env = gym.make('CartPole-v0')
        
    #重み、バイアスの合計
    def get_weights_size(self) :
        return self.nn.get_weights_size()
    
    #評価
    #weights : 重み
    def eval_NNGA(self ,weights ) :
        #重みの セット
        self.nn.set_weights(weights)
        
        #評価用シミュレート
        eval = self.eval_simulate()
        return eval
        
    #最終結果シミュレート
    #weights : 重み
    def finally_simulate(self , weights , filename, use_gym_render = False):
        #重みの セット
        self.nn.set_weights(weights)
        
        # シミュレート
        self.save_simulate( filename , use_gym_render)
    
    ##########################################################
    ### 倒立振子
    ##########################################################
    # 1stepシミュレート
    def simulate_1step(self , step ):
        inputs = np.reshape( self.state , (1 , INPUT_SIZE))
        
        #NNを使って推論 
        y = self.nn.predict(inputs)
        # 振り子アクション判定
        if(y[0][0] > y[0][1]):
            action = 0
        else:
            action = 1
            
        # 振り子シミュレート
        state , reward ,done ,dummy = self.env.step(action)
        state[2] =  rad_conversion(state[2])
        reward = -reward
        
        return state , reward

    #評価
    def eval_simulate(self ) :
        #振り子初期化
        self.state = self.env.reset()
        
        eval = 0
        for step in range(SIMULATE_STEP):
            self.state ,reward = self.simulate_1step(step )
            eval += reward
            
        return eval
    

    #シミュレート 動画作成
    def save_simulate(self ,  filename , use_gym_render) :
        #振り子初期化
        self.state = self.env.reset()
        
        eval = 0
        for step in range(SIMULATE_STEP):
            if (use_gym_render == True):
                self.env.render()
            
            self.state , reward=self.simulate_1step(step)
            eval += reward
            
            print('step={0}  state={1} reward={2}'.format(step,self.state,eval))
        
```

```mgg_test.py
import numpy as np

from nnga import NNGA
from Individual import Population
import crossover as cross
from mgg import MGG

# 交叉方法指定
crossType = cross.CrossoverType.UNDX
#致死個体処置方法
deadlyType = MGG.DeadlyIndividualType.Remove

#GAパラメータ
MAX_GENERATION  = 1000   #最大世代数
CROSS_TIMES     = 50    #交差回数
POPULATION_SIZE = 50    #個体サイズ

#MGG パラメータ
#エリート選択数
MGG_ELITE_SELECT_SIZE = 1

#ChoromIndex
CROME_SIZE = 1     #染色体(Chromosome) の 個数
VALUE_SIZE = 1     #評価値の数
cIndex = 0


if __name__ == "__main__":
    # NN
    nnga = NNGA()
    PARAMETER_SIZE = nnga.get_weights_size()
    
    # 評価関数に従った設定を取得
    dRange   = NNGA.get_range()
    evalType = NNGA.get_eval_type()
    
    mgg = MGG(dRange[0],dRange[1],dRange[0],dRange[1], crossType , 
              CROSS_TIMES , MGG_ELITE_SELECT_SIZE , deadlyType)

    #集団
    population = Population(POPULATION_SIZE , PARAMETER_SIZE , CROME_SIZE , VALUE_SIZE)
    
    # 交差方法初期化
    mgg.init_cross([PARAMETER_SIZE])

    #集団の初期化
    mgg.init_population(population)
    
    #初期集団 評価
    for ind in population.fIndividual:
        ind.fValue = nnga.eval_NNGA(ind.fChrom[cIndex])
    
    # 実行
    for gen  in  range(MAX_GENERATION+1):
        #MGGと交叉の実行
        mgg.make_children(population) 

        # 作成した子の評価
        for ind in mgg.fChildren.fIndividual:
            ind.fValue = nnga.eval_NNGA(ind.fChrom[cIndex])

        #MGGと交叉の実行
        mgg.select_children(population , evalType )

    #最終集団 評価
    for ind in population.fIndividual:
        ind.fValue = nnga.eval_NNGA(ind.fChrom[cIndex])
    #最終結果シミュレート
    bestIndex    = population.get_best_individual(evalType )
    nnga.finally_simulate (population.fIndividual[bestIndex].fChrom[cIndex] , 'sim' , True)
```

実数値GAの実装は[過去の記事](https://qiita.com/SwitchBlade/items/8f07b78f3922607a7ca3)と同じです。

#学習結果
##初期集団
初期集団の中の最良個体です。少しがんばりますが、どこかに行ってしまいます。
![NNGA_CartPoleGym_L4-[10]-2_g1000_UNDX_20190720_2244-0gen_Gym.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/438288/199c25aa-67eb-ec89-86b8-1eb9a8241ad9.gif)

##1000世代
うまくバランスが取れました！！
![NNGA_CartPoleGym_L4-[10]-2_g1000_UNDX_20190720_2244-F_Gym.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/438288/a0a8f70d-078a-45c3-ee06-80ec70cdccc9.gif)

開始早々20世代目には、最良個体は-200(最良値)になってます。
![NNGA_CartPoleGym_L4-[10]-2_g1000_UNDX_20190720_2244.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/438288/722de369-7dbd-8cad-6e8c-daf6bef77921.png)

致死個体(交差の結果、重み・バイアスが-5.0～5.0から外れた個体)は評価対象から外すとしているので、この間の評価回数は187回でした。遺伝的アルゴリズムはどうしても無駄な評価が増えるのですが、他の人のスコアと比べて、比較的ましな回数でした。


#まとめ
初期角度もほぼ垂直ですし、各stepでの外乱もないので、比較的簡単に解けてしまいました。
次はもうちょっとハードな環境でやってみたいと思います。。