1次遅れ系を制御するコントローラを強化学習で設計した振り返り


# はじめに

機械学習のアルゴリズムを評価するときに、パラメーターの調整やネットワークの設計といった試行錯誤にかける時間が、最近、不足しているな、と感じています。もっともらしい結果、期待した結果、議論に耐えられる結果が、たまたま得られると、もうそれで満足してしまって、結論を出して、報告し、アルゴリズムの評価を完了する、ということを繰り返してきました。これでは、試行錯誤のプロセスを効率よく回せるようにならないな、と反省しています。

議論に耐えられるような「キラキラ」した結果を報告しないといけない、「キラキラ」に至るまでの「ドロドロ」した試行錯誤の過程は、報告する必要がない、と思っています。そう考えているから、試行錯誤に対する振り返りが、おざなりのままになったのだろうな、と考えています。

ここでは、アルゴリズムを評価する途中の「ドロドロ」した試行錯誤の過程を、なるべく丁寧にすくい上げて白日の下にさらし、振り返りと知見の共有の機会にしたいと思います。

# ゴール

強化学習のアルゴリズムを応用して、ある制御対象を制御要件に合うようにコントローラを設計するタスクに取り組んだ活動の記録を報告します。タスクに取り組むプロセスの課題を挙げること、また、強化学習のアルゴリズムを応用する中で獲得した知見を報告します。

# タスク

このタスクでは、ある1次遅れ系を制御対象として、その出力を一定値に制御するコントローラを、Soft Actor-Critic(以下、SAC[1])により学習しました。最適レギュレータによる制御性能をベンチマークとしました。学習の目標は、このベンチマークに到達すること、としました。なお、SACは、強化学習のアルゴリズムの１つで、[tf_agents の実装](https://www.tensorflow.org/agents/api_docs/python/tf_agents/agents/SacAgent) を使いました。

# 制御対象

制御対象の状態方程式とコスト関数について説明します。

## 状態遷移

１次遅れ系の状態方程式を、以下の差分方程式で実装しました。

```math
\begin{align}
x(0) &= 0,\\
x(t+1) &= (1-\alpha)x(t) + \alpha u(t) + \beta w(t),~ t = 0 \cdots T-1,\\
y(t) &= x(t),~ t = 0 \cdots T,\\
\beta &= \sqrt{1 - (1-\alpha)^2}.
\end{align}
```
なお、x,u,y,wは、それぞれ、状態変数、操作変数、観測出力および外乱変数です。外乱は平均0、分散１の白色性正規雑音からサンプルしました。αはパラメータです。

この系は、時刻0から時刻Tまで、それぞれの時刻で操作と外乱の影響を受けながら、状態遷移します。この一連の時間発展を１つのエピソードとします。

## コスト関数

コスト関数を下記の通り設定しました。

```math
w\cdot\frac{1}{T}\sum_{t=1}^{T}y(t)^2 + (1-w)\cdot\frac{1}{T}\sum_{t=0}^{T-1}u(t)^2, {\rm where }~ w = 0.9.
```

# ベンチマーク

制御対象の１次遅れ系を、最適レギュレータにより制御したときのコスト関数をベンチマークとしました。なお、最適レギュレータは、参考文献[2]を参照しながら、設計しました。また、代数 Riccati 方程式は、[scipy.linalg.solve_discrete_are](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_discrete_are.html) を使って求解しました。


# コントローラ

コントローラを、[SAC agents](https://www.tensorflow.org/agents/api_docs/python/tf_agents/agents/SacAgent) により実装します。SAC agents を構成するネットワークのうち、actor network とcritic network は、それぞれ、操作量を決めるのアルゴリズムおよびコントローラの評価関数を近似する関数を実装します。actor network と critic network の設定を説明します。

## actor network

入力および出力を、それぞれ、設定値に対する観測出力の偏差および操作出力とします。偏差に一定のゲインをかけて、操作出力とします。このゲインをSac Agents により学習します。すなわち、ネットワークの入出力を線形のネットワークで実装し、ネットワークのパラメタを強化学習によって求めます。

## critic network

critic network は、３層のニューラルネットワークにより実装します。入力は、設定値に対する観測出力の偏差と操作出力の２変数です。すなわち、state-action value function を３層のニューラルネットワークにより近似します。

# 数値実験

数値実験の概要と結果を説明します。

## 手順

まず、数値実験の手順を説明します。

ある１つのコントローラを初期化してから、学習を終了するまでの手順は、以下の通りです。

1. コントローラのパラメタを初期化する。
2. replay buffer を初期化する。
3. コントローラを制御対象に閉ループ接続して、初期時刻から終了時刻まで、１エピソード分だけ、シミュレーションする。
4. シミュレーションする過程で、操作変数、観測出力、コストのログを、replay buffer に蓄積する。
5. エピソードが終了した後に、replay buffer に蓄積したデータで、コントローラのパラメタを更新する。
6. 工程 3.に戻る。
7. 工程3から6までを、16回繰り返す。

工程4から6は、１回のエピソードを完了するごとに、コントローラを学習する処理です。したがって、16回分のエピソードを実行しながら、コントローラを更新しました。

## パラメタ

工程5では、batch learning の処理を繰り返し呼び出すことで、コントローラのパラメタを更新します。batch learning を呼び出す回数（以下、反復回数）を、ハイパーパラメタ`num_train_iteration`とします。反復回数を8とした場合、32とした場合および128とした場合、それぞれで学習しました。反復回数が学習に与える影響を考察するためです。また、３つの反復回数、それぞれついて、5つのコントローラを独立に学習しました。したがって、合計15個のコントローラを学習しました。

## 結果

### 反復回数を32とした場合

まず、反復回数を32とした場合、すなわち、エピソードが終了するごとに、32回 batch learning のアルゴリズムを呼び出して、コントローラのパラメタを更新した場合の学習の結果を説明します。コントローラを５つ、それぞれ、独立に学習しました。コントローラそれぞれの学習曲線は、図1の通りになりました。横軸がエピソードの回数を表していて、縦軸がコストを表します。赤色で強調した曲線は、それぞれのエピソードで、5つのコントローラのコストを平均した値を示します。青色で強調した破線は、ベンチマークのコストです。以下に所見をまとめます。

* 11エピソードあるいは12エピソードまでに、いずれのコントローラも、概ね一定のコストに収束しました。
  * また、収束したコストは、ベンチマークのコストにおおむね一致しました。
  * 一定のコストに収束するまでの学習曲線のふるまいは、コントローラごとに異なりました。
  * しかし、5つのコントローラについてコストを平均すると、コストが単調に改善されました。
* 13エピソード以降、コントローラのコストのばらつきが大きくなりました。
  * また、コストの平均値は、ベンチマークのコストから徐々に乖離しました。

パフォーマンスが劣化した要因をまだ調べていませんが、次のように推測しています。(1)コントローラのパラメタが収束するに従って、それぞれのエピソードで、同じようなパターンの操作変数および観測出力のログがreplay buffer に蓄積し、actor network および critic network の学習に使うデータのばらつきが減少した。そして、(2)ばらつきが少ないデータで学習したことで、ネットワークがデータにoverfit し、コントローラのパフォーマンスにばらつきが発生した、と考えています。

13エピソード以降さらに学習を進めて、コストのばらつきが拡大と縮小を繰り返すようならば、この仮説を裏付けることができると考えます。

![fig_itr=32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2495166/35ad0a5e-93a6-4336-1ee5-cd186e908db4.png)

図1. 反復回数を32とした場合の学習曲線


### 反復回数を8とした場合

次に、反復回数を8とした場合の学習の結果を説明します。やはり、５つのコントローラをそれぞれ独立して学習しました。学習曲線は、図2の通りになりました。所見は以下の通りです。
* 16エピソード経過しても、コストが収束するまで、学習が進みませんでした。
* 5つのコントローラのコストを平均した値は、エピソードを重ねるに従って、概ね単調に改善されました。

エピソード１回あたり、32回反復した時に比べて、8回しか反復していないために、学習の速度が遅くなった、と解釈します。

![fig_itr=8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2495166/751d9932-b168-db72-0620-4da5cb8bc26c.png)

図2. 反復回数を8とした場合の学習曲線

### 反復回数を128とした場合

最後に、エピソードごとに128回繰り返しbatch learning のアルゴリズムを呼び出してコントローラを学習した結果を紹介します。コントローラを５つ、それぞれ独立に学習したときの学習曲線は、図3の通りになりました。所見は以下の通りです。

* ５つのコントローラのうち、３つのコントローラは、６エピソードまでに、ベンチマークのコストに収束しました。
  * ３つのコントローラのうち１つは、9エピソード以降、ハンチングを繰り返しながら、パフォーマンスが劣化しました。
* 残り２つのコントローラは、ハンチングを繰り返しながらコストを改善する傾向を示しました。

反復回数を8および32とした場合の学習曲線に比べて、不安定な挙動を示しました。要因を次のように推測しています。エピソード１回あたりのステップ数は100でした。すなわち、128という反復回数は、エピソード１回あたりのステップ数よりも大きな値でした。コントローラは、replay buffer から繰り返しデータを取り出して学習します。したがって、重複したデータを使ってコントローラを学習し、コントローラがデータに overfit し、パフォーマンスが不安定になったと解釈します。

![fig_itr=128.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2495166/c340161e-f585-53b0-9f6d-6aaddbc22146.png)

図3. 反復回数を128とした場合の学習曲線

### 結果と考察

#### 収束判定の必要性

ハイパーパラメタによって学習の経過が大きく異なるため、学習の収束を判定する基準を設定することが必要だと考えます。実際、反復回数を32とした場合には、あるエピソード以降、学習を進めると、かえってパフォーマンスが劣化します。また、反復回数を8とした場合には、コントローラを十分に学習するためには、元々設定した回数である16回以上、エピソードを繰り返す必要があります。さらに、今回の数値実験のようにベンチマークを利用できるとは限らないため、一般には、パフォーマンスの目標を学習する以前に設定できません。以上から収束を判定する条件がないと、学習のアルゴリズムやハイパーパラメタによるパフォーマンスの影響を分析するのが難しくなる、と考えます。

#### 学習曲線を目視で観察することと判断のバイアス

今回の数値実験の結果を振り返ってみると、反復回数は32が妥当だろうと考えます。学習曲線を観察して、コントローラのパフォーマンスがばらついていないか、単調に収束しているか、などを目視で確認すると、違和感なく判断できるかと思います。

個人的には、ハイパーパラメタを選ぶために学習曲線を目視で確認することが、アルゴリズムを評価する試行錯誤を非効率にする１つの要因である、と考えます。目視で判断することで、以下の例に示すように、判断にバイアスが入り込みやすくなる、と考えるためです。

* (1) まず、confirmation bias について考えます。confirmation bias は、仮説を裏付けるデータを重視して判断するバイアスです。学習を改善すると期待したハイパーパラメタの効果を検証する状況を想定します。目視で学習曲線を確認することで、効果があると期待したケースを良い結果であると判断しやすくなります。

* (2)次に、primacy bias について考えます。初めに観察した結果を重視して判断するバイアスです。反復回数=128の学習曲線を振り返ります。5つのうち、2つのコントローラは素早く収束し、パフォーマンスも安定していました。この2つのコントローラのいずれかを、一番初めに観測すると、反復回数=128が最も適切であると判断しやすくなります。

* (3)最後に、framing bias について考えます。一部の傾向を拡大して観察することで、全体の傾向を観察した場合と異なる判断をするバイアスです。反復回数=8では、学習が不十分であるために、学習曲線の前半部分だけを観察しました。学習の前半では、一般に、パフォーマンスが安定しないため、学習曲線にばらつきがあるように観察されます。そのため、反復回数=8の学習曲線は、反復回数=32の場合に比べて、必然的にばらつきます。このことは、どちらの反復回数を採用するべきか、今回の結果だけからでは、本来、判断できないことを示唆します。それにも関わらず、反復回数=8の学習曲線は、反復回数=32の場合に比べて、不安定であるという印象を与えるために、反復回数=32の方が妥当だと判断してしまいます。

# まとめ

機械学習のアルゴリズムを評価するプロセスを見直すために、強化学習のアルゴリズムによりコントローラを設計することで、ある一次遅れ系のプロセスを制御するタスクに取り組みました。結果として、(1)収束を判定する基準を検討する必要がある、また、(2)学習曲線を目視で判断することで判断のバイアスにつながり易い、という気づきを得ました。

皆さんは、どうやってバイアスを避けているのでしょうか。

# 参考文献

* [1] Haarnoja, T., Zhou, A., Abbeel, P., & Levine, S. (2018). Soft Actor-Critic: Off-Policy Maximum Entropy Deep Reinforcement Learning with a Stochastic Actor. ICML.
* [2] 萩原(1999) 『ディジタル制御入門』、コロナ社

# ソースコード

制御対象の一次遅れ系を実装したコードは以下の通りです。

```python
from copy import copy
from scipy.linalg import solve_discrete_are
import numpy as np


class Simulator(object):

    def __init__(self, T = 10, K = 1.0, discount = 0.9, maxIteration = 100, dv = 1.0, randomState = np.random.RandomState(), loggers = [list(),], costWeightOnState = 0.9):
        
        self._state = None
        self._u_prev = None
        self._episode_ended = False
        self._time = 0
        self.T = T
        self.K = K
        self.discount = discount
        self.maxIteration = maxIteration
        self.dv = dv
        self.beta = np.sqrt(1 - (1-1/T)**2)
        self.randomState = randomState
        self.loggers = copy(loggers)
        self.costWeightOnState = costWeightOnState

    def getActionBoundary(self):
        return {"minimum": (-1.,), "maximum": (1.,)}

    def getObservationBoundary(self):
        return {"minimum": (-1.,), "maximum": (1.,)}

    def getActionShape(self):
        return (1,)

    def getObservationShape(self):
        return (1,)
    
    def getObservation(self):
        boundary = self.getObservationBoundary()
        y = np.max((np.min((boundary["maximum"][0], self._state)), boundary["minimum"][0]))
        return np.array((y,), np.float32) # (1,)

    def reset(self):
        
        for logger in self.loggers:
            logger.clear()
        
        self._time = 0
        self._u_prev = 0.
        self._state = 0. 
        self._episode_ended = False

        for logger in self.loggers:
            logger.append((self._time, self._state, self._u_prev, self._episode_ended))
        
        return self.getObservation()
    
    def step(self, action):
        # action: (1,)
        
        if self._episode_ended == True:
            self.reset()
        
        cost = lambda x,u: np.sqrt(self.costWeightOnState * x**2 + (1-self.costWeightOnState) * u**2)

        u = action[0] # (,)
        w = self.randomState.randn() # (,)
        self._state = (1-1/self.T) * self._state + self.K/self.T * u + self.dv * self.beta * w
        self._u_prev = u
        self._time += 1
        
        reward = -cost(self._state, u)
        
        if self._time == self.maxIteration:
            discount = None
            self._episode_ended = True            
        else:                        
            discount = self.discount
            self._episode_ended = False        

        for logger in self.loggers:
            logger.append((self._time, self._state, self._u_prev, self._episode_ended))
        
        return self.getObservation(), reward, discount, self._episode_ended        
    
    def seed(self, seed):
        self.randomState.seed(seed)
    
    def getLoggers(self):
        return copy(self.loggers)

    def getGainLQR(self):
        
        A = np.array([[(1-1/self.T),],]) # (1,1)
        B = np.array([[self.K/self.T,],]) # (1,1)
        Q = np.eye(1) * np.sqrt(self.costWeightOnState) # (1,1)
        R = np.eye(1) * np.sqrt(1-self.costWeightOnState) # (1,1)
        
        P = solve_discrete_are(A,B,Q,R) # (1,1)
        K = - np.linalg.inv(R) @ B.T @ np.linalg.inv(A.T) @ (P-Q) # (1,1)
                
        return K
```