はじめに

これまで、強化学習のアルゴリズムを使ってコントローラを設計するときに、コントローラを静的なネットワークで、実装してきました。すなわち、現時刻の観測値を入力し、操作量を出力するネットワークで、コントローラを実装してきました。

操作量を算出するために時系列データを利用できる場合でも、自己回帰構造を持つネットワーク(以下、RNN)を使うことを敬遠してきました。いまから振り返ってみると、以下のような背景があると考えます。

静的なネットワークで実装したコントローラで、議論に耐えられる結果を出してきた
静的なネットワークでは目標を実現できなくて、RNNを使うことで克服したというタスクを経験したことがない
強化学習のタスクで、RNNのハイパーパラメタを調整した経験がない

現代制御理論のある枠組みでは、コントローラを、 (1)観測出力に基づいてプラントの内部状態を復元するカルマンフィルタと、 (2)内部状態から操作出力を決める状態フィードバックゲインの２つから構成します[1]。カルマンフィルタと状態フィードバックゲインはそれぞれ独立に設計されます。信号の流れに注目すると、コントローラは、まず、カルマンフィルタで、観測出力の過去の時系列データを処理し、次に、状態フィードバックで、内部状態から操作出力を返します。

RNNで実装したコントローラも、観測出力の過去の時系列データを入力して、操作出力を返します。ただし、強化学習のアルゴリズムに従ってRNNを学習する場合には、観測出力の系列から操作出力までのプロセスを、一気通貫で学習します。すなわち、現代制御理論の枠組みと異なり、カルマンフィルタと状態フィードバックを分けて設計しません。

このことから、システムの内部状態の一部を直接観測できないタスクに取り組むことで、静的なネットワークで実装したコントローラでは達成できないけれども、 RNNで実装した場合には達成できる課題を見つけることができるのでは、という疑問を持ちました。

ここでは、格子状に結合した振動系を制御対象として、強化学習のアルゴリズムでコントローラを設計するタスクに取り組んだ結果を報告します。

タスクの定義

制御対象

結合振動系を制御対象とします。振動系の一部の状態は測定されているけれども、残りの振動子の状態は観測されておらず、制御に利用できないものとします。

結合振動子の状態方程式を以下の通り定義します。

\begin{align*}
\dot{\theta}_{0} &= \omega_{0},\\
\dot{\theta}_{1} &= \omega_{1},\\
&\cdots,\\
\dot{\theta}_{n-1} &= \omega_{n-1},\\
J_{0}\dot{\omega}_{0} &= P_{\rm source} - E_{0}E_{1} B_{01}\sin(\theta_{0} - \theta_{1}),\\
J_{1}\dot{\omega}_{1} &= - E_{1}E_{2} B_{12}\sin(\theta_{1} - \theta_{2}) + E_{0}E_{1} B_{01}\sin(\theta_{0} - \theta_{1}),\\
&\cdots\\
J_{n-1}\dot{\omega}_{n-1} &= E_{n-2}E_{n-1} B_{n-2~n-1}\sin(\theta_{n-2} - \theta_{n-1}) - P_{\rm load},
\end{align*}

$n$は振動子の個数です。 $\theta_{0},\theta_{1},\cdots$, $\omega_{0},\omega_{1},\cdots$, $J_{0},J_{1},\cdots$, $E_{0}, E_{1}, \cdots$ は、それぞれ、各振動子の位相角、角周波数、慣性モーメント、電源電圧です。また、$B_{01}, B_{12}, \cdots$は、それぞれ、振動子と振動子の間のサセプタンスです。 $P_{\rm source}, P_{\rm load}$は、それぞれ、電源から投入する電力潮流および負荷が消費する電力です。

電源の電力潮流は、負荷電力と操作量$u$の和とします。すなわち、コントローラは、電源と負荷の電力潮流の差分を操作します。

\begin{align*}
P_{\rm source} = P_{\rm load} + u.
\end{align*}

操作の目標は、負荷が要求する需要電力の変動に対して、位相角の角周波数の変動を抑えることです。

コントローラは、電源の位相角と角周波数$\theta_{0}, \omega_{0}$を観測できるものの、電送路上の振動子の状態$\theta_{1},\cdots,\theta_{n-1}, \omega_{1},\cdots,\omega_{n-1}$を観測できない、と仮定します。各時点における観測出力の値もしくは過去の系列のデータを使って、コントローラは操作を決定します。

慣性モーメント、サセプタンス、電源電圧を、結合子全体にわたって、すべて1とします。

\begin{align*}
&J_{0} = J_{1} = \cdots = J_{n-1} = 1,\\
&E_{0} = E_{1} = \cdots = E_{n-1} = 1,\\
&B_{01} = B_{12} = \cdots = B_{n-2~n-1}= 1.\\
\end{align*}

評価関数

制御性能は、角周波数の変動の大きさをシミュレーションする期間にわたって平均したもので評価します。

\begin{align*} \tag{1}
\frac{1}{T}\int_{0}^{T}\dfrac{1}{n}\sum_{i=0}^{n-1}|\omega_{i}(t)|  {\rm d}t.
\end{align*}

ベンチマーク

最適レギュレータで制御したときの制御性能を、ベンチマークとします。

振動子の間の位相差が十分小さいと仮定して、正弦関数を１次式で近似し、さらに、すべての振動子の状態を観測できると仮定することで、最適レギュレータの設計法を応用し、状態フィードバックゲインを算出しました。

コントローラ

コントローラのネットワークとコスト関数を説明します。

ネットワークの構成について

コントローラを、 Soft Actor-Critic [2],[3]で実装しました。

静的なネットワークの場合には、 actor network および critic network をそれぞれ、以下の通り実装しました。


actor_net = actor_distribution_network.ActorDistributionNetwork(
	  observation_spec,
	  action_spec,
	  fc_layer_params=(),
	  continuous_projection_net=(
		  tanh_normal_projection_network.TanhNormalProjectionNetwork))

critic_net = critic_network.CriticNetwork(
		(observation_spec, action_spec),
		observation_fc_layer_params=None,
		action_fc_layer_params=None,
		joint_fc_layer_params=(32,),
		kernel_initializer='glorot_uniform',
		last_kernel_initializer='glorot_uniform')

RNNの場合には、以下の通りとしました。


actor_net = ActorDistributionRnnNetwork(
	  observation_spec,
	  action_spec,
	  lstm_size = (16,),
	  input_fc_layer_params = (),
	  output_fc_layer_params = (),
	  continuous_projection_net=(
		  tanh_normal_projection_network.TanhNormalProjectionNetwork))

critic_net = CriticRnnNetwork(
		(observation_spec, action_spec),
		lstm_size = (16,),
		observation_fc_layer_params=None,
		action_fc_layer_params=None,
		joint_fc_layer_params=(32,),
		kernel_initializer='glorot_uniform',
		last_kernel_initializer='glorot_uniform')

コスト関数の設定

２種類のコスト関数を実装しました。

１つ目のコスト関数は、すべての振動子にわたる角周波数の絶対値の平均値に、操作量の絶対値を加えたものです。負荷が変動したときに角周波数の動揺を抑えたいため、角周波数の大きさをコスト関数に加えました。また、コントローラの学習を正則化するために、操作量の絶対値もコスト関数に加えました。

\begin{align*}
c_{1}(t) = \dfrac{1}{n}\sum_{i=0}^{n-1} |\omega_{i}(t)| + |u(t)|
\end{align*}

２つ目のコスト関数は、１つ目のコスト関数に、振動子全体にわたる偏差の平均値（以下、位相のオフセット）の絶対値を加えたものです。それぞれの振動子の偏差は、角周波数を積分したものです。したがって、負荷変動後の角周波数の動揺を抑制するために、角周波数の積分もコスト関数に加えていることになります。

\begin{align*}
c_{2}(t) = \dfrac{1}{n}\sum_{i=0}^{n-1} |\omega_{i}(t)| + |u(t)| + \left|\dfrac{1}{n}\sum_{i=0}^{n-1} \theta_{i}(t)\right|
\end{align*}

数値実験

パラメタ

以下の表に示す３つのケース、それぞれで、コントローラを学習しました。

case	コスト関数	コントローラのネットワーク
case1	$c_{1}$	静的なネットワーク
case2	$c_{2}$	静的なネットワーク
case3	$c_{2}$	RNN

コスト関数の設定がコントローラの学習に与える影響を調べるために、 case1,case2ともにコントローラを静的なネットワークで実装し、 case2 のみ、コスト関数に位相のオフセットを加えました。

コントローラが過去の入力系列を利用できることによる効果を調べるために、コスト関数を共通の設定として、 case2 では静的なネットワークで、 case3 ではRNNで、コントローラを実装しました。

手順

以下の手順でコントローラを学習しました。

actor network および critic network をそれぞれ初期化する
replay buffer を初期化する
結合振動子を初期状態にする。
コントローラと結合振動子をつないで閉ループシミュレーションを開始する。
閉ループシミュレーションの間に、一定の時間間隔$2\pi/8$で操作量を更新する。
閉ループシミュレーションの間に、操作量を更新するごとに、操作量、結合振動子の観測出力、コストをreplay buffer に保存する。
閉ループシミュレーションが終了したら、actor network, critic network のパラメタを強化学習のアルゴリズムで更新する。
手順4に戻る。

手順3から6までの工程を１回のエピソードとします。数値実験では、エピソードを16回繰り返し、それぞれのエピソードが完了するごとに、コントローラのパラメタを更新しました。コントローラの初期条件や負荷変動を乱数で決めているため、学習の結果にばらつきが出ます。安定した結果を得るために、同一のパラメタで、5つのコントローラを、それぞれ独立に学習しました。

結果と考察

case1

図1(a)は、学習した5つのコントローラそれぞれの学習曲線を表しています。横軸はエピソードの回数を表していて、縦軸は、結合振動系全体にわたる角周波数の大きさを表します。定義は、式(1)を参照してください。赤の実線は、それぞれのエピソードにおける、5つのコントローラそれぞれの制御性能の平均値を表します。破線は、ベンチマークの制御性能を表します。

図1(b)は、16エピソード学習を完了した（５つあるコントローラのある１つの）コントローラで振動系を閉ループ制御したときのシミュレーション結果です。 10個の振動子のうち、電源($\theta_{0},\omega_{0}$)、中心($\theta_{5},\omega_{5}$)、負荷($\theta_{9},\omega_{9}$)、それぞれの周波数と位相の時間応答を表します。

特徴は以下の通りです。

図1(a)から、コントローラの学習はおおむね安定していて、8エピソード以降、制御性能が収束しています。
ベンチマークの制御性能に比べて、強化学習で設計したコントローラの制御性能には改善の余地があります。
図1(b)から、位相角が振動を伴いながら増加あるいは減少しています。これは、角周波数に定常偏差が残っていることを示唆します。

図1(a) 学習した5つのコントローラそれぞれの学習曲線

図1(b) コントローラで振動系を閉ループ制御したときのシミュレーション結果

case2

図2(a),(b)は、それぞれ、学習した5つのコントローラそれぞれの学習曲線と学習したコントローラで振動系を閉ループシミュレーションしたときの応答を表します。図の見方は、case1の説明をご覧ください。

特徴は次の通りです。

図2(a)から、コントローラの学習は、case1と同様に安定していて、16エピソード経過した後には、既に制御性能が収束しています。
ベンチマークの制御性能に比べると改善の余地がありますが、case1に比べると大分改善されています。
図2(b)から、位相角が発散せずに一定の範囲で振動しています。位相角は周波数を積分したものですから、この結果は周波数の定常偏差が0であることを示唆します。

図2(a) 学習した5つのコントローラそれぞれの学習曲線

図2(b) コントローラで振動系を閉ループ制御したときのシミュレーション結果

考察: コスト関数の設計について

case1, case2 ともに静的なネットワークでコントローラを実装しました。すなわち、電源の周波数と位相を入力として、電源から供給する電力潮流を操作しました。したがって、case1,case2　ともに、同一の種類のデータをコントローラに入力しています。

一方、case2では、コスト関数に、周波数の絶対値の平均に加えて、位相のオフセットの絶対値を加えました。周波数の変動を抑制するという目標であれば、 case1のように、コスト関数に周波数の絶対値が含まれていれば十分で、 case2のように、周波数を積分した位相のオフセットをコストに加える必要がないように思います。

実際には、位相のオフセットをコストに加えることで、はっきりとした効果が現れました。これは、次のように解釈できます。 (1) 外乱（負荷変動）がステップ状に変化するため、 (2)　周波数の偏差は、高周波数成分よりも、低周波成分に偏る。 (3) 位相のオフセットをコスト関数に加えることで、コストの周波数成分のうち、低周波成分に重みづけられた。 (4) コストに占める割合が相対的に大きい低周波数成分を重視してコントローラを設計することで、コストを効率よく抑えるコントローラを設計した。この解釈が正しければ、（低周波数成分に偏るといった）偏差の特徴を理解することが、コスト関数を設計するために必要である、と言えます。

以上から、コスト関数を、周波数の絶対値の平均といった制御要件から直感的に導いても、不十分である可能性を示しました。また、コストが低周波数成分に偏るといったドメイン知識を使ってコスト関数を設計すると、効率よくコントローラを学習できる可能性を示しました。特に、コントローラのネットワークを設計するときだけでなく、コスト関数を設計するときにも、偏差の低周波数成分をコストに含めるといった動特性に配慮が必要である可能性を示しました。

case3

学習した5つのコントローラそれぞれの学習曲線を図3(a)に、学習したコントローラで振動系を閉ループシミュレーションしたときの応答を図3(b)に示します。図の凡例については、case1を参照ください。

特徴は次の通りです。

図3(a)から、コントローラの学習は進んでいるものの、case1 や case2 に比べると、学習が不安定です。
また、5つのコントローラの間の制御性能のばらつきが、case1, case2 に比べて大きいです。
制御性能の平均は、case2 とおおむね同程度です。
図3(b)から、位相角は発散していません。これは、周波数の定常偏差を平均的に0にしたことを示唆します。

図3(a) 学習した5つのコントローラそれぞれの学習曲線

図3(b) コントローラで振動系を閉ループ制御したときのシミュレーション結果

考察: 制御対象の選び方について

RNNでコントローラを実装したcase3の制御性能が、静的なネットワークでコントローラを実装したcase2と比べて、よりベンチマークの制御性能に近づくことを期待していました。電源の角周波数と位相の過去の系列から、結合振動系の他の振動子の状態に関する情報を RNNで抽出することで、静的なネットワークにできない制御を実現できると期待したからです。実際には、case2 に比べて、case3で制御性能の改善が見られませんでした。

case2 の制御性能がベンチマークの制御性能に既に近いため、case3 でさらに改善できる余地が限られています。今回検討した結合振動系は、静的なネットワークに比べたRNNの利点を検討するのに適した例題ではなかった、と考えます。

また、RNNで実装したコントローラの学習が不安定になりました。 RNNの利点を検討する以前に、RNNで実装したコントローラを安定して学習する工夫が不十分だな、と考えています。

まとめ

静的なネットワークよりもRNNでコントローラを実装する方が適切なタスクを見つけるために、格子状に結合した振動系を制御対象として、強化学習のアルゴリズムでコントローラを学習するタスクに取り組みました。

結果として、コスト関数を、周波数の絶対値の平均といった制御要件から直感的に導いても、不十分である可能性を示しました。また、コストが低周波数成分に偏るといったドメイン知識を使ってコスト関数を設計すると、効率よくコントローラを学習できる可能性を示しました。特に、コントローラのネットワークを設計するときだけでなく、コスト関数を設計するときにも、偏差の低周波数成分をコストに含めるといった動特性に配慮が必要である可能性を示しました。今回検討した結合振動系は、静的なネットワークに比べたRNNの利点を検討するのに適した制御対象ではないことを示しました。

みなさんは、どういったタスクで、自己回帰構造を持ったネットワークを使って、制御しているのでしょうか。

参考文献

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

ソースコード

制御対象の結合振動系を実装したソースコードを共有します。


import numpy as np

class Oscillator(object):
    '''
    classdocs
    '''


    def __init__(self, n_oscillator, dt, simulation_period):
        '''
        Constructor
        '''
        
        self._n_oscillator = n_oscillator
        self._b = 1.0
        self._dt = dt
        self._cnt = None
        self._t = None
        self._x = None        
        self._nsteps = int(simulation_period/dt)
        
    def reset(self):        
        # (2*n_oscillator,)
        
        self._cnt = 0
        self._t = 0.
        self._x = np.zeros(2*self._n_oscillator) # (2*n_oscillator,)
        
        self.Theta = np.zeros((self._nsteps+1, self._n_oscillator)) # (*+1, n)
        self.Omega = np.zeros((self._nsteps+1, self._n_oscillator)) # (*+1, n)
        self.Psource = np.zeros(self._nsteps) # (*, 1)
        self.Pload = np.zeros(self._nsteps) # (*, 1)        
        self.Time =  np.zeros(self._nsteps+1) # (*+1,)
        
        self.Time[self._cnt] = self._t
        self.Theta[self._cnt,:] = self._x[:self._n_oscillator]
        self.Omega[self._cnt,:] = self._x[self._n_oscillator:]
        
    def step(self, p_source, p_load):
        
        self.Psource[self._cnt] = p_source
        self.Pload[self._cnt] = p_load
        
        k0 = self.f(self._x, p_source, p_load)
        k1 = self.f(self._x + k0 * self._dt/2, p_source, p_load)
        k2 = self.f(self._x + k1 * self._dt/2, p_source, p_load)
        k3 = self.f(self._x + k2 * self._dt, p_source, p_load)

        self._cnt += 1        
        self._t += self._dt    
        self._x += (k0+2*k1+2*k2+k3)/6 * self._dt # (2*n_oscillator,)

        self.Time[self._cnt] = self._t
        self.Theta[self._cnt,:] = self._x[:self._n_oscillator]
        self.Omega[self._cnt,:] = self._x[self._n_oscillator:]
        
        done = self._cnt == self._nsteps
        
        return done
        
    def f(self, x, p_source, p_load):        
        # x: (n_oscillator*2), x[:n_oscillator] ... theta, x[n_oscillator:] ... omega, 
    
        n_oscillator = self._n_oscillator
        b = self._b
        
        theta = x[:n_oscillator]
        omega = x[n_oscillator:]
        
        dxdt = np.zeros(2*n_oscillator)
        
        # d(theta)/dt = omega
        dxdt[:n_oscillator] = omega # (n_oscillator,)
    
        # d(omega)/dt = ...
        tmp = b * np.sin(theta[:-1] - theta[1:])
        dxdt[n_oscillator] += p_source
        dxdt[n_oscillator:-1] += -tmp
        dxdt[-1] += -p_load
        dxdt[n_oscillator+1:] += tmp
        
        return dxdt
        
    def get_state_eq_expression(self):
    
        # x = [theta, omega]
        # dxdt = [df/dtheta, df/domega]
        
        n_oscillator = self._n_oscillator
        b = self._b
        
        A = np.zeros(shape=(2*n_oscillator, 2*n_oscillator)) # (2n, 2n)
        B = np.zeros(shape=(2*n_oscillator,1)) # (2n,1)
        
        A[:n_oscillator,n_oscillator:] = np.eye(n_oscillator,n_oscillator)
        A[n_oscillator:,:n_oscillator] = -b * (2 * np.eye(n_oscillator,n_oscillator) - np.eye(n_oscillator,n_oscillator,1) - np.eye(n_oscillator,n_oscillator,-1))
        A[n_oscillator,0] = -b
        A[2*n_oscillator-1,n_oscillator-1] = -b    
        B[n_oscillator] = 1
        
        return A, B # (2*n_oscillator, 2*n_oscillator), (2*n_oscillator, 1)