2022-07-07

FRITによるデータ駆動型制御器チューニング

この記事ではFRITの手続きと「なぜそれで上手く行くのか?」を軽く軽く説明します。

1.はじめに
2.FRITとは?
3.おわりに
参考文献

1.はじめに

人間、10年・20年も生きていれば何度も制御器のパラメータを調整する機会があると思います。~~(要出典)~~

自分もNHK ロボコンで延々とPIDのパラメータチューニングをしていた記憶があります。
人力で調整するのって嫌じゃないですか? 私は嫌です。

この記事では自動でパラメータを調整する手法について説明します。

2.FRITとは?

FRIT(Fictitious Reference Iterative Tuning)とは,
実際のシステムへの入力と出力の1組の時系列データのみを用いて,
オフラインで制御器のパラメータ調整を行う手法.

閉ループ系を安定化する制御器のパラメータを一つ見つければ,
そのパラメータでの入出力データを用いて, 目標応答に近い応答をする制御器のパラメータを計算できる.

何度も制御対象に入力を加えてデータを回収したりする事無くパラメータを調整できる所が嬉しい.
閉ループ系での応答でパラメータ計算を行うので, 単体では不安定な制御対象にも適用できる.
非線形計画になる事と初期パラメータが必要な事が難点.

類似した手法としてVRFTが知られる.

問題設定

下のFig. 1のようなシステムを考える.

制御対象 $G$ は線形時不変で動特性は未知, 制御器 $C$ は $C(\rho, z) = \dfrac{b_\mu z^\mu + + ... + b_1 z + 1 }{a_\nu z^\nu + ... + a_1 z + a_0}$ として、
$\rho = \begin{bmatrix} a_0 & ... & a_\nu & b_1 & ... & b_\mu \end{bmatrix}$ を調整可能なパラメータとする.

このシステムの応答を目標応答 $y_d = T_d r$ に近づける事を目的とする.
( $T_d$ : 目標応答伝達関数)

FRITの具体的な手続き

初期パラメータを $\rho_{ini}$ とする.

1. 制御器 $C(\rho_{ini})$ によってFig. 1のシステムを制御した時の実際の入出力データを収集する
(入出力データをそれぞれ $u_{ini} := u(\rho_{ini}),\ \ y_{ini} := y(\phi_{ini})$ と置く. )

(※ 参照信号 $r$ から出力 $y$ までの閉ループ伝達関数は $T(\rho) = \dfrac{ GC(\rho) }{ 1+GC(\rho) }$ と書けるので,
対象への入力 $u$ 及び出力 $y$ についても $\rho$ をパラメータとした関数と見なせる.)

2. 評価関数 $J_F(\rho) := ||y_{ini} - T_d \tilde{r}(\rho)||^2$ を最小化するパラメータ $\rho$ を求める.
ただし疑似参照信号 $\tilde{r}(\rho)$ を $\tilde{r}(\rho) := C(\phi)^{-1}u_{ini} + y_{ini}$ と定義する.

以上の手続きがFRITのアルゴリズムである.

なぜこれで上手く行くのか(概略)

$y_{ini} (= T(\rho) \tilde{r}(\rho))$ に $T_d \tilde{r}(\rho)$ を近づける事で, $T_d$ と $T(\rho)$ を近づける感じのイメージ.

まず疑似参照信号 $\tilde{r}(\rho)$ については, $y_{ini} = Gu_{ini}$ より, 一般の $\rho$ について
$T(\rho)\tilde{r}(\rho) = \dfrac{GC(\rho)}{1+GC(\rho)} \tilde{r}(\rho) = \dfrac{Gu_{ini} }{1+GC(\rho) } + \dfrac{ GC(\rho) }{1+GC(\rho) }y_{ini} = \dfrac{ 1 + GC(\rho) }{1+GC(\rho) }y_{ini} = y_{ini}$ が成立.

① $J_F(\rho^*)=0$ が成立する場合
$J_F(\rho^*) = 0$ なる $\rho^*$ については, $y_{ini} = T_d \tilde{r}(\rho^*)$ が成立するので,
$y_{ini} = T(\rho^*) \tilde{r}(\rho^*) = T_d \tilde{r}(\rho^*)$ となる.

同じ信号 $\tilde{r}(\rho^*)$ に対して同じ出力 $y_{ini}$ を出力している事から,
システムの伝達特性も $\tilde{r}(\rho)$ の帯域で同じ特性を持つと言える.

以上から, $C(\rho^*)$ を実装して, 元の参照信号 $r$ を印加する事で $y(\rho^*) = \dfrac{ GC(\rho^*) }{ 1+GC(\rho^*) } \simeq T_d r$ となる.

② $J_F(\rho) = 0$ となる $\rho$ が存在しない場合
これは $T_d$ が $T(\rho) = \dfrac{GC(\rho)}{1+GC(\rho)}$ で表せる $T(\rho)$ の集合に含まれない場合に対応する.
(例えばPID制御器を用いると, その表現力ゆえにほとんどの場合には $T_d$ は上記の集合に含まれない. )

$\tilde{r}(\rho) = \dfrac{1}{T(\rho)}y_{ini}$ を用いて, $J_F(\rho) = \| y_{ini} - T_d \tilde{r}(\rho) \|^2 = \left\| \left(1 - \dfrac{T_d}{T(\rho)} \right) y_{ini} \right\|^2$ と書き換えると,

これは初期データの下での $T_d$ と $T(\rho)$ の相対的誤差の最小化と解釈できる.

また, Parsevalの定理を用いて変換すると, 離散時間フーリエ変換は $z = e^{i\omega}$ としたZ変換に対応する事から,
$J_F(\rho) = \dfrac{1}{2\pi} \int_{-\pi}^\pi \left\| \left( 1 - \dfrac{T_d(e^{i\omega})}{T(\rho, e^{i\omega})} \right) Y_{ini}(e^{i\omega}) \right\|^2 d\omega$ とできるので,
$Y_{ini}(e^{i\omega})$ が大きくなる周波数帯域での $T_d$ と $T(\rho)$ の相対的誤差の最小化と解釈できる.

この事から, パラメータの計算は初期データの周波数特性 $Y_{ini}(e^{i\omega})$ に大きく影響を受ける事が言える.

3.おわりに

対になる手法としてVRFTが知られています. (こちらは開ループ系に適用. 二次計画で解ける.)

Parsevalの定理周りがあまり詳しくないので説明が怪しいかもです.
(前提が抜けてたり, 間違いがあったりしたら教えて下さい)

(追記)
2023年のキャチロボ(ロボット競技)で実戦投入し,
円筒座標系アームの $r, \theta$ 制御用モーターの位置制御PIDゲインの調整に使用しました.

直感的な話ではありますが, 与えた目標応答に応じて良いパラメータが手軽に得られ,
(PIDなので)その後の調整も容易な非常に使い勝手の良い手法だと思います.

こちらの記事の実装が参考になります.
qiita.com

参考文献

www.jstage.jst.go.jp
(↑ここの説明を要約したものがこの記事に対応しています.)

2021-12-23

MCMCの画像のノイズ除去への応用

adventar.org
この記事はKyoto Universuty Advent Calendar 2021の23日目の記事です.

元々は逆強化学習について記事を書こうと考えていたのですが,
時間が無かったので別の内容について雑な記事を書きました. 申し訳ございません.
授業か何かでやった内容が面白かったので記事にまとめました.

はじめに
背景というか定義
MCMC
- Metropolis-Hastings法
- Gibbs sampling
MCMCを応用した画像のノイズ除去
- 原理
- 結果と実装
雑な感想と近況報告

はじめに

詳細釣り合いと極限分布の関係を雑に雑に述べて, 次にMCMCの代表的な手法を説明する.
最後にMCMCの応用として画像のノイズ除去を述べる.

背景というか定義

概要しか書けていません. 申し訳無いです...

・離散時間マルコフ連鎖
定義 : 状態空間 $\mathbb{S}$ 上の確率過程 $\{X_n:n \in \mathbb{Z}\}$ が $\forall n \in \mathbb{N}$ と状態列 $j_0,...,j_n \in \mathbb{S}$ に対して,
$P(X_{n+1}=j | X_n=i, X_{n-1}=j_{n-1},...,X_0 = j_0) = P(X_{n+1}=j|X_n = i)$ を満たす時,
$\{X_n\}$ を状態空間 $\mathbb{S}$ 上の離散時間マルコフ連鎖と呼ぶ.

※ マルコフ性 : 次の状態 $X_{t+1}$ が現在の状態 $X_t$ にのみ依存する事

・推移確率行列
定義 : 行列 $\mathbf{P}=(P_{i,j})_{i,j \in \mathcal{S}} \geq \mathbf{O},\ \ P_{i,j}=P(X_{n+1}=j|X_n = i)$ をマルコフ連鎖 $\{X_n\}$ の推移確率行列と呼ぶ.
要は $i \rightarrow j$ に遷移する確率を要素 $P_{i,j}$ とする行列の事.

推移確率行列 $\mathbf{P}$ については以下が成立する事が知られている.
1. $\mathbf{P}$ の各行の和は $1$ となる. ( $\mathbf{Pe}=\mathbf{e}$ )
2. $\forall n, m \in \mathbb{Z}_+$ に対して, $\mathbf{P}^{n+m} = \mathbf{P}^n\mathbf{P}^m$ が成立

2つめはチャップマン=コルモゴロフ方程式と呼ばれる. これはマルコフ性から証明できる.

・既約(連結クラス)
定義 : 状態部分集合 $\mathbb{C \subset S}$ が既約とは, $\forall (i,j)\in \mathbb{C}^2$ に対して, $i$ から $j$ に到達可能という事.
雑にまとめるとどの状態間も行き来できるという事.

・再帰的
定義 : $\tau_j$ を状態 $j$ への初到達時刻と置く. $P(\tau_i < \infty | X_0 = i) = \sum_{n=1}^{\infty} P(\tau_i = n|X_0 = i)=1$ が成立する時,
状態 $i$ は再帰的であると言う. 任意の状態が再帰的な時, マルコフ連鎖は再帰的と言う.

※ 再帰的かつ $\mathbb{E}[\tau_i|X_0 = i$ < \infty]も成立する時, 状態 $i$ は正再帰的であると言う.
※ 規約なマルコフ連鎖に置いては, 「ある状態 $i$ が(正)再帰的 $\Rightarrow$ 全ての状態が(正)再帰的」が成立.
(互いの状態の到達可能性から示せる.)

・周期性
定義 : 状態 $i$ に対して, $\{ n \in \mathbb{N}\ |\ P(X_n = i|X_0 = i) = \mathbf{P}^{(n)}_{i,i} > 0\}$ の最大公約数 $d_i$ を状態 $i$ の周期と呼ぶ.

※ 周期 $d_i$ が $1$ に等しい時, 非周期的であると言う.
※ 互いに到達可能な状態は周期が同じになる.

・定常分布
定義 : $\mathbf{\pi P } = \mathbf{\pi}$ を満たす確率ベクトル $\mathbf{\pi} = (\pi_i)_{i \in \mathbb{S}}$ の事.

※ 常に存在するとは限らない.
※ 規約かつ再帰的なマルコフ連鎖は定常分布を持つ.
※ 定常分布と極限分布は異なる(定常分布に到達するとは限らない)

・極限分布
定義 : 任意の初期分布 $\mathbf{\pi}^(0)$ に対して, $\lim_{n \rightarrow \infty}\mathbf{\pi}^{(n)}$ が同じ確率ベクトル $\mathbf{p} = (p_s)_{s \in \mathbb{S}}$ に収束する時,
$\mathbf{p}$ を極限分布と呼ぶ. 推移確率行列がエルゴート的(規約かつ非周期的で再帰的な事)な時には極限分布と定常分布が一致する. ( $\lim_{n\rightarrow \infty}P^{(n)}_{i,j} = \pi_j,\ \ \ \forall (i,j) \in \mathbb{S}$ )

・詳細釣り合い
定義 : 推移確率行列 $\mathbf{P}$ に対して, $\mu_i P_{i,j} = \mu_j P_{j,i}$ なる方程式の事.
これが成立する時, $\sum_{i \in \mathbb{S}} \mu_i P_{i,j} = \sum_{i \in \mathbb{S}} \mu_j P_{j,i} = \mu_j$ が成立するので
$\mathbf{\mu}$ は $\mathbf{\mu P} = \mathbf{\mu}$ を満たす. (これは定常測度である事に注意)

これは, 詳細釣り合いを満たす様にマルコフ連鎖を構成してやれば目的の分布を定常分布に持つ様にできる事を意味している. MCMCはこれを利用している.

MCMC

求めたい分布を定常分布に持つマルコフ連鎖を用いて, 分布からのサンプルを生成する手法のこと.
十分に分布が収束するまでの期間(バーンイン)のサンプルデータは破棄する.

Metropolis-Hastings法

===================================
提案分布を $h$ とし, 近似したい分布の密度関数を $g$ とする. エルゴード性を仮定.
0. $\mathbf{X}^{(0)} := \mathbf{x}$ : 初期値
1. $\mathbf{V}^{(t)}|\mathbf{x}^{(t)} \sim h(\mathbf{v}^{(t)}|\mathbf{x}^{(t)})$
2. $a^{(t)}=\min\left\{1, \frac{h(\mathbf{X}^{(t)}|\mathbf{V}^{(t)})g(\mathbf{V}^{(t)})}{h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)})g(\mathbf{X}^{(t)})} \right\}$
3. $\left\{\begin{array}{lll} P(\mathbf{X}^{(t+1)} = \mathbf{V}^{(t)}) = a^{(t)} \\ P(\mathbf{X}^{(t+1)} = \mathbf{x}^{(t)}) = 1 - a^{(t)} \end{array}\right.$ として, $\mathbf{V}^{(t)}$ を移動先の点として受け入れるか確率的に決定する. -> 1に戻る.
===================================

詳細釣り合いの成立について,
1. $\mathbf{x}^{(t)} = \mathbf{x}^{(t+1)}$ の時, 自明
2. $\mathbf{x}^{(t)} \neq \mathbf{x}^{(t+1)}$ の時, $p(\mathbf{V}^{(t)}|\mathbf{X}^{(t)}) = \min\left\{1, \frac{h(\mathbf{X}^{(t)}|\mathbf{V}^{(t)})g(\mathbf{V}^{(t)})}{h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)})g(\mathbf{X}^{(t)})} \right\}h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)})$ と書ける.
なので,
$\begin{equation}\begin{split}p(\mathbf{V}^{(t)}|\mathbf{X}^{(t)})g(\mathbf{X}^{(t)}) &= \min\left\{1, \frac{h(\mathbf{X}^{(t)}|\mathbf{V}^{(t)})g(\mathbf{V}^{(t)})}{h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)})g(\mathbf{X}^{(t)})} \right\}h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)}) g(\mathbf{X}^{(t)}) \\ &= \min\{h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)}) g(\mathbf{X}^{(t)}),\ h(\mathbf{X}^{(t)}|\mathbf{V}^{(t)})g(\mathbf{V}^{(t)})\} \\ &= p(\mathbf{X}^{(t)}|\mathbf{V}^{(t)})g(\mathbf{V}^{(t)}) \end{split}\end{equation}$
が成立する.
以上より詳細釣り合いが示された. また, この時の分布は $\mathbf{X}^{(t)}\sim g(\mathbf{x})$ に従う.

Gibbs sampling

===================================
$\mathbf{x} = [(\mathbf{x})_1,...,\ (\mathbf{x})_m$ ], $(\mathbf{x})_{-i} := [(\mathbf{x})_1,...,(\mathbf{x})_{i-1},(\mathbf{x})_{i+1},...,(\mathbf{x})_m$ ]と置く. エルゴード性を過程
1. 適当に $i$ を選択する. (ランダムとか順番とか)
2. $(\mathbf{X}^{(t+1)})_i | \mathbf{x}^{(t)}\ \sim\ \ g((\mathbf{x})_i|(\mathbf{x}^{(t)})_{-i}), \ \ \ (\mathbf{x}^{(t+1)})_{-i} = (\mathbf{x}^{(t)})_{-i}$ として定める. --> 1.に戻る
(条件付き分布 $g((\mathbf{x})_i|(\mathbf{x})_{-i}) := \frac{g(\mathbf{x})}{\int g(\mathbf{x})d(\mathbf{x})_i}$ である. )
===================================

※ これはMH法の特殊なケースに対応.
$h_i([(\mathbf{x})_1,...,(\mathbf{x})_{i-1},(\mathbf{v})_i,(\mathbf{x})_{i+1},...,(\mathbf{x})_m]|\mathbf{x}) = g((\mathbf{v})_i | (\mathbf{x})_{-i}) = \frac{g(\mathbf{v})}{\int g(\mathbf{x})d(\mathbf{x})_i}$
これに対して, $\frac{h_i(\mathbf{x|v})g(\mathbf{v})}{h_i(\mathbf{v|x})g(\mathbf{x})} = \frac{g((\mathbf{x})_i|(\mathbf{v})_{-i})g(\mathbf{v})}{g((\mathbf{v})_i|(\mathbf{x})_{-i})g(\mathbf{x})} = \frac{\frac{g(\mathbf{x})}{\int g(\mathbf{x})d(\mathbf{x})_i}g(\mathbf{v})}{\frac{g(\mathbf{v})}{\int g(\mathbf{x})d(\mathbf{x})_i}g(\mathbf{x})} = 1$ より,
選択確率が1のMH法に対応している事が言える.

※ 実際の応用では1の $i$ の選択は番号順にやる事が多いらしい.

MCMCを応用した画像のノイズ除去

膨張・収縮を使った方が実用的レベルでは速いと思います.

原理

ノイズが乗る前の画像(未知)を $H$ , ノイズが乗った後の画像(既知)を $V$ と置く.
尤度関数 $P(H|V)$ に従って画像を生成.

同時確率は $P(H, V) = \frac{1}{Z}\exp \left( \alpha \sum_{i} \sum_j v_{ij} h_{ij} + \beta \sum_{i',j' \in N(i,j)} h_{ij}h_{i'j'} \right)$ として定義する.
ただし, $Z$ は正規化項, $N(i,j)$ は $(i,j)$ の近傍. $\alpha, \beta$ はハイパーパラメータ.

画像データの配列の要素は $\pm 1$ が入るので, この様に定義すると, 周囲の点と同じ値を持つほど確率が上がり,
持たないほど下がる様に同時確率を構成できる.

結果と実装

ギブスサンプリングを用いて尤度関数を近似する.
結果は以下の通り.

↑ノイズを載せた画像(よく見るアレ)

↑ノイズ除去後の画像

実装は以下.

import numpy as np
from tqdm import tqdm
import matplotlib.pyplot as plt

def load_image(filename):
	tmp = plt.imread(filename)
	img_gray = 0.299*tmp[:,:,2] + 0.587*tmp[:,:,1] + 0.114*tmp[:,:,0]  # grayscale
	img = np.where(img_gray>0.5, 1, -1)
	fixed_img = np.zeros([img.shape[0]+2, img.shape[1]+2])
	fixed_img[1:-1, 1:-1] = img
	return fixed_img

# V:ノイズ有りの画像, H:元の画像
# 点(i,j)でのP(H|V) := 1/Z exp(v[i,j]*h[i,j] + \sum_{k,l \in N(i,j)}h[i,j]h[k,l] )  ただしN(i,j)は(i,j)の近傍
# H[i,j]をh_ijと書く事にすると, P(h_ij=1 |V,H_N(i,j))が求まる.

def denoise(noised_pict, total_step, burnin_step, alpha, beta):
	ans = np.random.choice([1, -1], size=noised_pict.shape)
	x_max, y_max = ans.shape[0], ans.shape[1]
	for _ in tqdm(range(burnin_step+total_step)):
		for i in range(1, x_max-1):
			for j in range(1, y_max-1):
				tmp = beta*np.sum([ans[i-1+l][j-1+m] for l in range(3) for m in range(3)]) + alpha*noised_pict[i][j]
				prob_h = np.exp(tmp)/(np.exp(tmp)+np.exp(-tmp))  # H[i,j]=1となる条件付き確率
				ans[i][j] = 2*(np.random.random_sample() < prob_h) -1.0
	return ans

if __name__ == "__main__":
	file_name = "./img_noised.png"
	tmp = load_image(file_name)
	plt.imshow(tmp, cmap="gray")
	plt.show()
	ans = denoise(tmp, 1000, 100, 1.0, 0.5)
	plt.imshow(ans, cmap="gray")
	plt.show()

雑な感想と近況報告

詳細釣り合いと定常分布の関係を書いた記事があまり多くなかった気がするので書きました.

2021年はずっとロボコンに追い回されていました. 夏休みを返して...(恨み)

2021-08-10

Soft Actor-Critic(SAC)でdonkeycarのシミュを走らせてみる

強化学習アルゴリズム

東大の松尾先生の研究室主催のセミナー(RLSP2021)の復習に書いた記事です.
1回生で参加し, 知識不足ではありましたがとても勉強になりました. ありがとうございました.

f:id:Emile11235:20210810191909p:plain

記号の定義
SACとは?
environmentについて
実装
感想など
出典,参考など

記号の定義

$\mathcal{S, A}$ : 状態, 行動
$p: \mathcal{S \times A} \rightarrow \mathbb{R}$ : 遷移確率関数
$\pi : \mathcal{S \rightarrow A}$ : 方策
$r : \mathcal{S \times A} \rightarrow \mathbb{R}$ : 報酬関数

SACとは?

概要

エントロピー項 $H(X)$ を追加した目的関数 $J(\theta) := \mathbb{E}[\sum_{t} r(s_t, a_t) + \alpha H(\pi_\phi(\cdot | s_t))]$ を用いる手法.
状態空間, 行動空間が連続な場合にも適用可能.

※大抵は目的関数は $J(\theta) := \mathbb{E}[\sum_{t} r(s_t, a_t)]$
※エントロピー項の定義は $H(X) := - \sum_{x \in U} P_X(x) \log P_X(x)$

Criricは $Q$ 関数( $Q_\theta : \mathcal{S \times A} \rightarrow \mathbb{R}$ )を近似し, Actorは $\pi_\phi$ ( $Q$ 関数のsoftmax方策)を近似する.
ActorはGaussian Policyを用い, 分布の平均値と分散を出力する事で行動の分布を生成する.

学習とそれぞれの損失関数

・Critic( $Q_\theta$ )の損失関数

Criticはエントロピー項を考慮した $Q$ 関数を近似する.
具体的にはベルマン作用素 $B^\pi Q_\theta(s_t,a_t) := r(s_t, a_t) + \gamma \mathbb{E}_{s_{t+1}\ \ \sim p}[ V(s_{t+1}) ]$ を用いて更新する.
※状態価値関数は $V(s_t) := \mathbb{E}_{a_t \ \sim \pi}[ Q(s_t, a_t) - \alpha \log \pi(a_t | s_t) ]$ と定義される

このベルマン作用素はエントロピーの最大化を含む様に修正されているが,
$r'(s_t, a_t) := r(s_t, a_t) - \alpha \log \pi (a_t | s_t)$ を報酬関数として用いると解釈すれば,
通常のベルマン期待作用素と同じ形になり, 同等の性質が言える.
(targetの方策に基づいた $Q$ 関数に収束. 方策の評価をすると解釈できる.)

以上から, Criticは損失関数として
$J_Q(\theta) := \mathbb{E}_{(s_t,a_t)\ \sim D} [ \frac{1}{2} (Q_\theta (s_t, a_t) - (r(s_t, a_t) + \gamma \mathbb{E}_{s_{t+1}\ \ \sim p} [ V_{\bar{\theta}} (s_{t+1}) ] )\ )^2 ]$
を用いて学習する. (ベルマン作用素を適用した $Q$ 関数に最も近くなるように最小二乗法をかけてる)

※ $\bar{\theta} := \min_{\theta^{target}_i}\ \ [ V_{\theta^{target}_i}\ \ (s_{t+1}) ] = \min_{\theta^{target}_i}\ \ [ Q_{\theta^{target}_i}\ \ (s_t,a_t) - \alpha \log \pi_\phi (a_t | s_t) ]$ であり,
double networkとTarget Q-networkと呼ばれる手法を用いている.
(不確かな時は楽観的にというヒューリスティックと急激な変化を抑える意図がある)

・Actor( $\pi$ )の損失関数

$\pi$ は $Q$ 関数のsoftmax方策を近似する.
方策は $\pi_{new} := \mathrm{argmin}_{\pi' \in \Pi} D_{KL}(\pi'(\cdot | s_t) || \frac{\exp(Q_\theta \ \ (s_t, \cdot) )}{Z_\theta (s_t)})$ として更新する.
(方策の改善と解釈できる.)

また, $\mathbb{E}_{s_t \ \sim D}[ D_{KL}(\pi'(\cdot | s_t) || \frac{\exp(Q_\theta \ \ (s_t, \cdot) )}{Z_\theta (s_t)}) ] = \mathbb{E}_{a_t\ \ \sim \pi}[ \alpha \log \pi_\phi(a_t | s_t) - Q_\theta (s_t, a_t) ] + const.$
より, Actorは損失関数として $J_\pi(\phi) = \mathbb{E}_{a_t\ \ \sim \pi}[ \alpha \log \pi_\phi(a_t | s_t) - Q_{\theta'} (s_t, a_t) ]$ を用いて学習する.

※ $Z_\theta$ はsoftmax関数の正規化項, $\theta' := \min_{\theta_i}\{ Q_{\theta_i}(s_t, a_t) \}$ .

ただし, $a_t \sim N(\mu_t, \sigma_t)$ として行動を生成すると微分が出来なくなるので,
行動を $\epsilon_t \sim N(0, I)$ を用いて, $a_t := \mu_t + \epsilon_t \circ \sigma_t$ として行動を生成する事で, 微分できる様にする.
(reparameterization trickと呼ばれる手法)

全体の流れ

SACではExperience ReplayとTarget Q-Networkを用いる事に注意.
ActorとCriticの関係は以下の様になる.
f:id:Emile11235:20210810021217p:plain

全体の流れとしては,
1. 現在の方策から行動を生成し, $(s_t, a_t, s_{t+1}, r_t)$ をバッファ $B$ に収集.
2. $\{s_t, a_t, s_{t+1}, a_t\}$ をバッファ $B$ からサンプリング.
3. Q関数を損失関数 $J_Q(\theta)$ を用いて, $\theta \leftarrow \theta - \lambda_Q \nabla_{\theta_i} J_Q(\theta_i) \ \ \ (i=1,\ 2)$ として更新
4. 方策を損失関数 $J_\pi(\phi)$ を用いて, $\phi \leftarrow \phi - \lambda_\pi \nabla_\phi J_\pi(\phi)$ として更新
5. 何回かQ関数と方策を更新したら,
target Q-functionも $\theta^{target}_i \leftarrow \tau \theta^{target}_i + (1-\tau)\theta^{target}_i \ \ \ \ (i = 1,2)$ として更新

environmentについて

gym-donkeycarを利用した.
github.com

ドキュメントに書いてある仕様と実際の仕様が異なっていた気がします.
(どこか忘れてしまいました.)

実装

この実装では, envから得られる状態にVAEで前処理をした物を3frame分集めて1つのstateにした.
また, actionも3回連続で入力した.

(VAEはRLSP2021で同じチームだったikedaさん(Ryunosuke Ikeda (@ImR0305) | Twitter)が実装してくださいました. ありがとうございます.)

github.com

f:id:Emile11235:20210810033059g:plain
↑走行時の映像

※VAEはEncoderとDecoderを持っておいて入力をEncodeした物をDecodeした物と入力の誤差が小さく成る様にEncoder, Decoderを学習させる事で入力をEncodeする事で次元削減できるという手法.

f:id:Emile11235:20210815214554p:plain
↑criticのlossのグラフ

f:id:Emile11235:20210815214558p:plain
↑actorのlossのグラフ

f:id:Emile11235:20210815214602p:plain
↑報酬のグラフ

何かおかしいです...
(何が原因か分かる方がいらっしゃれば教えてください...!)

感想など

RLSP2021はとても楽しかったし勉強になった.

出典,参考など

1. 元論文
arxiv.org
2. 松尾研のスプリングセミナーの授業スライド

2020-12-21

A2Cの理論と実装

強化学習

この記事はKyoto University Advent Calendar 2020の記事の12/21日分として書かれた物です.
adventar.org
先日のd0ra1998先輩の記事には改めて京都の四季の美しさに気付かされました.
先輩方の面白い記事が並んでいる中でこんな駄文を投下するのは気が引けますが....

A2Cってなに?
どういう仕組み?
- Actor-Criticとは?
- 損失関数の構造
実装
結果
嵌った所
参考にした物
ポエム

A2Cってなに?

A3Cよりあとに提案された,Advantage, Actor-Criticの2つを用いる手法です.(Aが2つなのでA2C.)
experience replayを用いないので, RNN を使えるといったメリットもあります.

学習の概要は,

一つのネットワークに対して複数のAgentを準備しておいて,
(それぞれのAgentに対して)行動選択 -> 新しいState, Reward を観測 -> (規定回数まで)行動選択に戻る
-> 蓄積した行動,State,Rewardを用いて学習

を繰り返すという流れです.

Agentが共通というのがA3Cとの違いです.(とても実装が楽)
それ以外の学習自体はA3Cと全く同じ.

どういう仕組み?

Actor-Criticとは?

行動選択を行うActor関数と評価を行うCritic関数の2つを用いて学習する手法です.

今回はモデルとしてニューラルネットワークを用いて,
Actorは $\pi(a|s;\theta )$ を近似する様に,Criticは $V^{\pi}(s; \theta_v)$ を近似する様に学習する様にします.

一般には $\theta ,\ \theta_v$ は別々に準備するのですが,
今回は画像を入力に用いるので,途中までは同じ層を用いて最後に分岐させるようにします.

損失関数の構造

$Value\ Loss := \frac{1}{2} \sum \{ (\sum_{i=0}^{n-1} \gamma^i r_{t+i} + \gamma^n V(s_{t+n}; \theta_v) ) - V(s_t; \theta_v) \}^2$ を
$V(s; \theta_v )$ の教師信号に用います.

これはforward viewに基づいた $n$ -step rewardを使用していて,
1つの報酬がそれより前の状態にも伝わるので,報酬に時系列的なずれのあるような環境で有効です.

( $n$ -stepのTD誤差と同じです.)

$Policy\ Loss := - \sum \{ A^{\pi}(s) \log (\pi (a|s) ) \} + \beta H(\pi)$ をActorの教師信号に用います.

$A^{\pi}(s) := (\sum_{i=0}^{n-1} \gamma^i r_{t+i} + \gamma^n V(s_{t+n}; \theta_v) ) - V(s_t; \theta_v)$ はAdvantageと呼ばれていて,

方策勾配定理より,
$f(\theta) = \sum_{s \in \mathcal{S}}p_{s_0}(s)V^{\pi_{\theta}}(s)$ とした時の方策勾配が,一般の関数 $b: \mathcal{S} \rightarrow \mathbb{R}$ について,
$\ \ \ \nabla_{\theta }f(\theta) = \mathbb{E}^{\pi}[ \nabla_{\theta} \log{\pi(a|s; \theta)(Q^{\pi}(s,a) - b(s))}$ ]となる事を利用しています.

この $A^{\pi}$ は, $Q^{\pi}(s,a) - V^{\pi}(s)$ を近似していて,行動のその状態での相対的良さを意味しています.
$\ \ V^\pi(s)$ を引く事で方策勾配の分散を抑えられて学習が安定するというメリットもあります.

後ろの項は,エントロピー項と呼ばれ,
$H(\pi) := - \sum \pi(a|s) \log \pi(a|s)$ と定義されます.

この項を $Policy\ Loss$ に加える事で,局所最適解に学習初期段階で収束する事を防ぐ事が期待できます.
$\ \ \beta$ はハイパーパラメータで,この項の寄与の具合を調整します.

今回は共通のニューラルネットを用いて学習するので,
上記の項をまとめて,

$Loss := \alpha Value\ Loss + Policy\ Loss \\ = \frac{ \alpha }{2} \sum \{ (\sum_{i=0}^{n-1} \gamma^i r_{t+i} + \gamma^n V(s_{t+n}; \theta_v) ) - V(s_t; \theta_v) \}^2 - \sum \{ A^{\pi}(s) \log (\pi (a|s) ) \} + \beta H(\pi)$

$\ \$ を全体の損失関数として使用します.

実装

Open AIのBreakOutを環境に使用しました.
詳しくは下のコードを読んでください.(main2.py, model.py, env.py を使用しています.)
github.com

結果

直近50回のスコア平均が310を安定して超える様になりました.

嵌った所

1.Advantage LossのAdvantage項に勾配を流さない

2.Value LossのAdvantage項の $V(s_{t+k})$ 項に勾配を流さない

3.Adamで学習する時には,各LayerをHe normalで初期化すると良い

4.rewardのクリッピング

5.reset後の開始時にランダムな長さだけ,stateを動かした方が学習が偏らない

~~6.あちこちで実装をミスっていた~~

参考にした物

実装が綺麗で分かりやすい記事です. (本当にありがとうございました)
horomary.hatenablog.com

損失関数等がまとまってて見やすいです.
medium.com

A3Cの論文,それなりにお気持ちfullでスッキリします.
https://arxiv.org/pdf/1602.01783.pdf

お勉強につかったので.
www.kspub.co.jp

ポエム

ハイパーパラメータや細かい実装の工夫に学習が左右されるのが厄介ですね.

ここまで読んでくださってありがとうございました.

2020-12-13

DDPGでPendulum-v0を解く

強化学習

もし間違っている所やおかしいなと思う所があれば指摘してください

DDPGとは
実装
参考文献

DDPGとは

概要

DDPGは,Actor-Criticの構造を取り,
学習には soft-target や Experiments Replay といった手法を用いる手法で,連続値での制御を行います.

$Critic$ では $Q(s,a | \theta^Q)$ を学習し, $Actor$ では, $\pi_d (s| \theta^\pi)$ を学習します,

soft-target は fixed-target に含まれる手法であり,

今回であれば $Actor,\ Critic$ 両方のネットワークにそれぞれ target ネットワークを準備しておきます.

詳しくは後で述べますが, Actor は確率的方策ではなく決定論的方策を学習していて,
行動方策には出力値にノイズを載せた物を使います.

この特徴から, off-policy な手法であると言えます.

※Experiments Replayは,行動方策と学習方策の一致する on-policy な Agentの学習には使えません.
(以前の履歴は以前の方策に基づいているので.)

今回の実装では,ノイズには平均 0 のガウシアンノイズを用いました.

決定論的な方策勾配定理

証明は~~面倒くさい~~長いのでこの記事では説明しません.

結果だけを述べると,

方策勾配は, $\mathbb{E}_{s_t \sim \rho^\pi} [\nabla_{\theta^\pi} Q(s,a | \theta^Q)| s = s_t, a = \pi_d (s_t | \theta^\pi)]$ として近似されます.

なので, $\theta^{\pi}$ の学習は, $\theta^\pi_{t+1} = \theta^\pi_t + \alpha \mathbb{E}[ \nabla_{\theta^\pi} Q(s_t, \pi_d(s_t | \theta^\pi)) ]$ によって行われます.
( $\alpha$ : 学習率)

詳細は下記の論文を見てください.

arxiv.org

損失関数

$Critic$ の損失関数には
$L_{critic} = \sum_{t} \{ r_t + \gamma Q (s_{t+1}, \pi_d(s_{t+1}| \theta^{\pi_{target}}) | \theta^{Q_{target} }) - Q(s_t, a_t | \theta^Q) \}^2 )$ を用います.

$Actor$ の損失関数には,
$L_{actor} = - \sum_{t} Q(s_t, \pi_d(s_t | \theta^\pi) | \theta^Q )$ を用います.

これを $\theta^\pi$ について微分してやれば,先程の方策勾配が求まる事に注意してください.

実装時には,
1. $target$ ネットワークには,勾配を流さない様にする
2.方策勾配の計算時に求まった勾配で $Critic$ を更新しない.

に気をつけてください.

soft-target

fixed-targetでは一定間隔で fixed-target と本体のネットワークを完全に同期していましたが,

この手法では,

$\theta^{Q_{target}} \longleftarrow \tau \theta^{Q} + (1 - \tau) \theta^{Q_{target}}$

$\theta^{\pi_{target}} \longleftarrow \tau \theta^{\pi} + (1 - \tau) \theta^{\pi_{target}}$

によって緩やかにネットワークを同期するようにします.
(同期は,ネットワークの更新と同じ位高頻度で行います.)

これによって学習が安定化されるそうです.

実装

Pendulum-v0の問題設定

状態は, $( \cos(\theta),\ \sin(\theta),\ \dot{\theta} )$ の組で与えられ,

行動は, $\ a\ (-2 \leq a \leq 2)$ の連続値で与えます.

また,報酬は, $r = - (\theta^2 + 0.1 \dot{\theta}^2 + 0.001 a^2)$ で計算されます.

これまでの記事で扱った問題との最大の違いは行動を連続値で扱わなくてはならない事です.

ソースコード

github.com

main.pyとmodel.py に書かれています.

結果

学習回数は200回程度です.

参考文献

arxiv.org

↓方策勾配の導出 (分かりやすい)
sykwer.hatenablog.jp

↓実装が分かりやすい
horomary.hatenablog.com

2020-08-20

ベルマン作用素の基本と性質

強化学習アルゴリズム

この記事上での定義
- ベルマン方程式再び
ベルマン作用素とその性質
- 2つのベルマン作用素
  - 定義
- 性質
まとめ

間違っていたら教えてください

この記事上での定義

・作用素 : 関数から関数への写像
・ $\mathbb{R}^{S} : S \longrightarrow \mathbb{R}$
・ $S$ : 状態集合
・ $A$ : 行動の集合
・ $g(s,a): S \times A \rightarrow \mathbb{R}$ : 状態 $s$ で行動 $a$ を取った時に手に入る報酬を与える関数.
・ $p_t(s'|s,a) : S \times A \times S \rightarrow [0,1]$ : 状態 $s$ で行動 $a$ を取った時に,状態 $s'$ に遷移する確率を与える関数.
・ $\pi(s,a) : S \times A \rightarrow [0,1]$ : 方策の事.状態 $s$ で行動 $a$ を取る確率を与える関数.

ベルマン方程式再び

定義だけ書きます.
ベルマン方程式 : $V^{\pi}(s) =\displaystyle{\sum_{a \in A} \pi(s,a)(g(s,a)+\gamma \sum_{s' \in S}p_t(s'|s,a)V^{\pi}(s'))}$
ベルマン最適方程式 : $V^*(s) = \displaystyle{\max_{a \in A}\{g(s,a) + \gamma \sum_{s' \in S}p_t(s'|s,a)V^*(s')\}}$

ベルマン作用素とその性質

定義 $\rightarrow$ 性質と証明 $\rightarrow$ の順番に書きます.

2つのベルマン作用素

MDP $M=\{S,A,p_{s_0},p_t,g\}$ に対してのお話である事に注意してください.

定義

ベルマン期待作用素 $B_{\pi}:\mathbb{R}^S \rightarrow \mathbb{R}^S$ は,
$(B_{\pi}V)(s) := \displaystyle{\sum_{a \in A}\pi(a|s)(g(s,a)+\gamma \sum_{s' \in S}p_t(s'|s,a)V(s')),\ \forall s \in S}$ と,

ベルマン最適作用素 $B_* : \mathbb{R}^S \rightarrow \mathbb{R}^S$ は,
$(B_*V)(s) := \displaystyle{ \max_{a \in A}\{ g(s,a) + \gamma \sum_{s' \in S}p_t(s'|s,a)V(s') \},\ \forall s \in S }$ と定義します.

ベルマン方程式たちとそっくりな見た目です.
この様に定義する嬉しい理由はこれらの作用素の持つ性質に有ります(まとめに書いてあります).

性質

1.単調性

$v(s) \leq v'(s),\ \ \forall s \in S$ がならば,
・ $(B^n_* v)(s) \leq (B^n_* v')(s),\ \ \forall s \in S,\forall n \in \mathbb{N}$ と,
・ $(B^n_{\pi} v)(s) \leq (B^n_{\pi} v')(s),\ \ \forall s \in S,\forall n \in \mathbb{N}$ が成立.

証明(上の式)

数学的帰納法を用いて成立を示す.
成立を示す式を①とする.
(1) $n = 0$ の時, $v(s) \leq v'(s),\ \ \forall s \in S$ は明らかに成立するので,①の成立が示された.

(2) $n = k (\neq 0), k \in \mathbb{N}_0$ の時の成立を仮定し, $n = k+1$ の時の成立を示す.

$p_t(s'|s,a) > 0$ より, $p_t(s'|s,a)(B^k_*v)(s') \leq p_t(s'|s,a)(B^k_*v')(s'),\ \ \forall s \in S$ の成立が言える.
ゆえに, $\gamma \geq 0$ である事を踏まえて,
$\displaystyle{ g(s,a) + \gamma \sum_{s' \in S}p_t(s'|s,a)(B^k_*v)(s') \leq g(s,a) + \gamma \sum_{s' \in S}p_t(s'|s,a)(B^k_*v')(s'),\ \ \forall s \in S }$
より,
$\displaystyle{ \max_{a \in A}\{g(s,a) + \gamma \sum_{s' \in S}p_t(s'|s,a)(B^k_*v)(s') \} \leq \max_{a \in A} \{g(s,a) + \gamma \sum_{s' \in S}p_t(s'|s,a)(B^k_*v')(s') \} }$ ,
つまり $(B^{k+1}_*v)(s) \leq (B^{k+1}_*v')(s),\ \ \forall s \in S$ が成立する.
以上より, $n=k+1$ でも成立.

以上から,数学的帰納法より,①の成立が示された.
下の式も全く同様に証明できる.

2.ベルマン作用素はバラせる.

関数 $v$ に定数 $b$ を加えた関数を $(v+b):\mathbb{R}^{S}$ と定義する.
任意の $s \in S,n \in \mathbb{N}_0$ に対して,
・ $(B^n_*(v+b))(s) = (B^n_*v)(s) + \gamma^n b$ ,
・ $(B^n_{\pi}(v+b))(s) = (B^n_{\pi}v)(s) + \gamma^n b$ が成立.

証明(上の式)

$(v+b)(s) = v(s) + b,\ \ \forall s \in S$ が成立する事に留意して,数学的帰納法を用いて示す.
成立を示したい式を①と置く.

(1) $n = 0$ の時,
$(v+b)(s) = v(s) + b$ より,明らかに①は成立する.

(2) $n = k$ の時の成立を仮定して, $n = k+1$ でも成立する事を示す.
$(B^k_*(v+b))(s) = (B^k_*v)(s) + \gamma^k b,\ \ \forall s \in S$ が成立するから,これの両辺に $B_*$ を適用して,
$(B^{k+1}_* (v+b))(s) = B_*(B^{k}_*v + \gamma^k b)(s) = (B^{k+1}_* v)(s) + \gamma^{k+1}b,\ \ \forall s \in S$ が成立する.
ゆえに, $n = k+1$ でも成立する.

以上より,数学的帰納法から,①の成立が示された.
下の式も全く同様に証明できる.

3.収束性

任意の有界の状態関数 $s : \mathbb{R}^S$ に対して,
(1) $V^*(s) = \lim_{k \rightarrow \infty}(B^k_{*}v)(s),\ \forall s \in S$
(2)非定常な方策系列 $\Pi':=\{\pi_0,\pi_1,...,\pi_{k-1}\} \in \Pi_k^M$ のベルマン期待作用素 $B^k_{\Pi' }$ について,
$\Pi := \{\Pi',\pi_k,\pi_{k+1},...\} \in \Pi_k^M$ に対して, $V^{\Pi}(s) = \displaystyle{ \lim_{k \rightarrow \infty}(B^k_{\Pi}v)(s),\ \ \forall s \in S}$ が成立.

証明:

(1)関数 $v$ は $V^*$ の有界関数なので,
$|V^*(s)-(s)| \leq b,\ \ \forall s \in S$ なる定数 $b \in \mathbb{R}$ が存在する.
よって, $V^*(s) - b \leq v(s) \leq V^*(s) + b$ が成立.

この両辺に $B_*$ を $k$ 回適用して,単調性より,
$B_*^k(V^*-b)(s) \leq (B_*^kv)(s) \leq B_*^k(V^*+b)(s),\ \ \forall s \in S$ ,さらに変形して,
$(B_{*}^{k}V^{*})(s)-\gamma ^{k}b \leq (B_{*}^{k}v)(s) \leq (B_{*}^{k}V^{*})(s)+\gamma ^{k}b$ が成立.

これに対して, $k \rightarrow \infty$ とすれば,はさみうちの原理より, $\displaystyle{ \lim_{k \rightarrow \infty}(B^k_*v)(s) = V^*(s) },\ \ \forall s \in S$ が成立.
以上から,示された.

(2)
$V^{\Pi}(s) = \mathbb{E}[ \sum_{t=0}^{k-1}\gamma^tg(S_t,A_t)|S_0 = s ] + \mathbb{E} [ \sum_{t=k}^{\infty}\gamma^t g(S_t,A_t)|S_0 = s ],\ \ \forall s \in S$ より,
$g(s,a)$ は有界なので, $|g(s,a)| \leq R_{\rm{max}},\ \forall s \in S,\forall a \in A$ なる $R_{\rm{max}}$ が存在するから,
$\displaystyle{ -\sum_{t=k}^{\infty}\gamma^t R_{\rm{max}} \leq V^{\Pi}(s) - \mathbb{E}[\sum_{t=0}^{k-1}\gamma^t g(S_t,A_t)|S_0 = s ] \leq \sum_{t=k}^{\infty}\gamma^t R_{\rm{max}},\ \ \forall s\in S}$ が成立.
ゆえに, $\displaystyle{ \frac{-\gamma^k R_{\rm{max}}}{1-\gamma} \leq V^{\Pi}(s) - \mathbb{E}[\sum_{t=0}^{k-1}\gamma^t g(S_t,A_t)|S_0 = s ] \leq \frac{\gamma^k R_{\rm{max}}}{1-\gamma},\ \ \forall s \in S }$ が成立.

また, $(B_{\Pi}^kv)(s)$ は式の意味から,
$\displaystyle{ (B_{\Pi}^kv)(s) = \mathbb{E}[\sum_{t=0}^{k-1}\gamma^t g(S_t,A_t) + \gamma^k v(S_k)|S_0 = s ] = \mathbb{E}[\sum_{t=0}^{k-1}\gamma^t g(S_t,A_t)|S_0 = s ] + \gamma^k \mathbb{E}[v(S_k)|S_0 = s] }$
$v$ は有界なので, $|v(s)| \leq V_{max},\ \ \forall s \in S$ なる定数 $V_{max} \in \mathbb{R}$ が存在.
ゆえに, $\displaystyle{ -\gamma^k V_{max} \leq (B_{\Pi}^kv)(s) - \mathbb{E}[\sum_{t=0}^{k-1}\gamma^t g(S_t,A_t)|S_0 = s ] \leq \gamma^k V_{max},\ \ \forall s \in S }$ が成立.

以上の2つを合わせて,
$\displaystyle{ - \gamma^k(V_{max} + \frac{R_{max}}{1-\gamma} ) \leq (B_{\Pi}^kv)(s) \leq \gamma^k (V_{max} + \frac{R_{max}}{1- \gamma } ) }$ が成立するから, $k \rightarrow \infty$ とすれば,
はさみうちの原理より(2)の成立が示された.

4.一意性

(1)ベルマン最適方程式の解になる関数 $v:\mathbb{R}^S$ は,
$(B_*v)(s) = v(s),\ \forall s \in S$ を満たし,この様になる $v$ は $V^*$ のみである.
(2)定常な方策 $\pi$ についてのベルマン期待方程式の解になる関数 $v : \mathbb{R}^S$ は,
$(B^{\pi}v)(s) = v(s),\ \forall s \in S$ を満たし,この様になる $v$ は $V^{\pi}$ のみである.

証明
(1)
$\exists s\in S,V_{*}(s)\neq v'(s)$ かつ $(B_{*}v')(s)=v'(s)$ を満たす $v':S\rightarrow \mathbb{R}$ の存在を仮定する.
この時, $B^*v'=v'$ より, $\displaystyle{ V^*(s) = \lim_{k \rightarrow \infty}(B_*^kv')(s) = v'(s),\ \forall s \in S}$ が成立.

これは仮定に矛盾する.ゆえに,示された.

(2):(1)と全く同じ様に示せる.
ちなみに, $\pi$ が定常方策でなければならないのは,
定常では無い時,同じ状態 $s$ であっても時間ステップによって $p_t(s'|s,a)$ が異なるからである.

5.縮小性

任意の有界関数 $v:S\rightarrow \mathbb{R},v':S\rightarrow \mathbb{R}$ と $k\in \mathbb{N_0},0 \leq \gamma < 1$ に対して,
(1)ベルマン最適作用素 $B_{*}$ について,
$\max_{s\in S}|(B_*^kv)(s)-(B_*^kv')(s)|\leq \gamma^{k}\max_{s\in S}|v(s)-v'(s)|$
が成立.
(2)任意の $\pi \in \Pi$ のベルマン最適作用素 $B_{\pi}$ について
$\max_{s\in S}|(B^k_{\pi}v)(s)-(B^k_{\pi}v')(s)| \leq \gamma^{k}\max_{s\in S}|v(s)-v'(s)|$
が成立.

証明
(1) $\alpha := \min_{s \in S}\{v(s) - v'(s)\},\ \beta := \max_{s \in S}\{v(s)-v'(s)\}$ と置く.
式の意味から, $v'(s) + \alpha \leq v(s) \leq v'(s) + \beta$ が成立する.
この両辺に $B_*$ を繰り返し適用して
$(B_*^kv')(s) + \gamma^k \alpha \leq (B_*^kv)(s) \leq (B_*^kv')(s) + \gamma^k \beta$ が成立.

ここで, $\epsilon := \max\{ |\alpha|,|\beta|\}$ とすると,
$(B_*^kv')(s) - \gamma^k \epsilon \leq (B_*^kv)(s) \leq (B_*^kv')(s) + \gamma^k \epsilon$ が成立.
ゆえに, $-\gamma^k \epsilon \leq (B_*^kv)(s) - (B_*^kv')(s) \leq \gamma^k \epsilon$ が成立する.

以上より, $| (B_*^kv)(s) - (B_*^kv')(s) | \leq \gamma^k \epsilon = \gamma^k \max_{s\in S}|v(s)-v'(s)|$ が成立する事が示された.

(2)(1)とまっっったく同じ.

ちなみに,このような写像( $B^*$ とか)の事を縮小写像と呼ぶ.
$L_{\infty}$ ノルムで考えれば $\max$ が付いている理由が分かる.

まとめ

・上に書いた5つの性質から,有界な関数 $v: \mathbb{R}^S$ にベルマン最適作用素やベルマン作用素を適用する事で,
$v$ が $V^*$ や $V^{\pi}$ に収束する事が分かる.

・これらの性質が考察の基礎になる.(例えば,近似ベルマン作用素が本当に収束するかを考える時とか)

参考文献　
www.kspub.co.jp
↑とてもわかり易く説明されています.

2020-06-20

MountainCar-v0をfixed target Q-networkを用いて解く

強化学習

もし間違っている所やおかしいなと思う所があれば指摘してください

MountainCar-v0の問題設定
fixed target Q-networkについて
- 定義再び
- どういう手法か
その他の工夫
- Experience Replay
実装
結果

MountainCar-v0の問題設定

目的:台車を右上の旗がある所に到達させる.

$Action$ は離散値で $0(push left),1(no push),2(push right)$ の3種類.

$State$ は連続値で $position([-1.2,0.6$ )]と $velocity([-0.07,0.07$ )]の2種類.

ここを見るとわかりやすい.
github.com

fixed target Q-networkについて

定義再び

報酬関数 $g: S \times A \longrightarrow \mathbb{R}$ .
状態価値関数 $V_n^{\pi}(s): S \longrightarrow \mathbb{R}$ は $V^{\pi}(s):=\sum_{t=0}^{\infty}\gamma^t \mathbb{E}^{\pi}[g(S_t,A_t)|S_0=s ]$
行動価値関数 $Q_n^{\pi}(s):S \times A \longrightarrow \mathbb{R}$ は $Q^{\pi}(s,a) := \sum_{t=0}^{\infty} \gamma^t \mathbb{E}^{\pi}[g(S_t,A_t)|S_0=s,A_0=a ]$
と表される.

上の定義から,
$Q_n^{\pi}(s,a) = \mathbb{E}^{\pi}[g(S_t,A_t) + \gamma V_n^{\pi}(S_{t+1}) |S_t=s,A_t=a ]$ と変形できるので,

状態価値関数 $V_n^{\pi}(s)$ の更新にベルマン最適作用素 $B_*$ を用いて
$V_{n+1}^{*}(s) = V_{n+1}^{\pi}(s) = B_*V_{n}^{*}(s) = \max_{a \in A}\{ \mathbb{E}^{\pi}[g(S_t,A_t) + \gamma V_n^*(S_{t+1}) |S_t=s,A_t=a] \}$
と更新する様にすると, $V_{n+1}^{\pi}(s) = \max_{a\in A}Q_n^{\pi}(s,a))$ が成立する.

定義から, $Q_{n+1}^{\pi}(s,a) = \mathbb{E}^{\pi}[g(S_t,A_t) + \gamma \max_{a' \in A}Q_n^{\pi}(S_{t+1},a') | S_t=s,A_t=a]$
と更新できる.

どういう手法か

上のQ値の更新式, $Q_{n+1}^{\pi}(s,a) = \mathbb{E}^{\pi}[g(S_t,A_t) + \gamma \max_{a' \in A}Q_n^{\pi}(S_{t+1},a') | S_t=s,A_t=a]$ に対して,
少し工夫をする.

具体的に言うと,Q値を算出するのに用いるニューラルネットワークを2種類準備する.
そのうちの一方(train network)は更新が発生する度に更新をするのに対して,もう一方(target network)は基本的には更新しない.
ただし,target networkは一定周期(今回の実装では,１つのEpisodeが終了するタイミング)でtarget networkとtrain networkを同期させる.

train networkのQ値を $\overset{\sim}{Q}$ ,target networkのQ値を $\underset{\sim}{Q}$ と表す様にすると,
$\overset{\sim}{Q}$ の学習の教師信号には, $g(s,a)+\gamma \max_{a' \in A}\underset{\sim}Q(s',a')$ を用いる.

この様にする事で, Q値が不安定になることを防げ, 学習が早くなる.

その他の工夫

Experience Replay

学習に用いるデータをランダムにシャッフルすることで,時系列順のステップに生じる強い相関を緩和する手法.
確率的勾配降下法と同じお気持ち.

実装

import gym
import numpy as np
from gym import wrappers  # gymの画像保存
import tensorflow as tf
from tensorflow import keras
from collections import deque
from collections import namedtuple
from keras.models import Sequential
from keras.layers import Dense
from keras.optimizers import Adam
import random

train_repeat_times = 2000
Exp = namedtuple("Experience", ["s", "a", "r", "n_s", "d"])
# state, act ,reward , next state, done
#  observe, rew, done, info = self.env.step(action)

class Agent:
    actions = [0, 1, 2]
    gamma = 0.99
    eps = 1.0
    eps_decay = 0.05
    eps_min = 0.01
    iteration_times = 201  # game is 201steps.
    buffer_get_size = 32  # bufferの大きさが32になったらtrain.
    exp = deque(maxlen=20000)
    env = gym.make("MountainCar-v0")
    state = env.reset()
    episode_times = 0  # 現在のepisode数

    def __init__(self):
        self.train_network = self.create_network()  # 毎回更新する.(fixed target Q-network)
        self.target_network = self.create_network()  # 1episodeごとに更新する.
        print(self.target_network.inputs)
        self.target_network.set_weights(self.train_network.get_weights())
        # self.env = wrappers.Monitor(self.env, "/home/emile/Videos/", video_callable=(lambda ep: ep % 50 == 0))

    def create_network(self):
        hide_size = 16
        random.seed(0)
        model = Sequential()
        print(self.env.observation_space)
        model.add(Dense(hide_size, activation='relu', input_shape=self.env.observation_space.shape))
        model.add(Dense(hide_size*2, activation='relu'))
        model.add(Dense(self.env.action_space.n, activation='linear'))
        model.compile(loss='mse', optimizer=Adam(lr=0.001))
        model.summary()
        return model

    def eps_change(self):
        self.eps = max(self.eps-self.eps_decay, self.eps_min)

    def get_best_action(self, state):
        return np.argmax(self.predict2(state)[0])

    def predict2(self, states):  # 評価値を算出.
        stat = np.vstack([[states]])
        return self.train_network.predict(stat)

    def get_policy_action(self, state):  # epsilon-greedy
        tmp = random.uniform(0.0, 1.0)
        if tmp > self.eps:
            return self.get_best_action(state)
        else:  # random
            return np.random.randint(0, self.env.action_space.n)

    def update2(self):
        if len(self.exp) < self.buffer_get_size:
            return
        experiences = random.sample(self.exp, self.buffer_get_size)
        #self.exp.clear() <- ここを外したら上手く動き始めた...
        states = np.vstack([e.s for e in experiences])
        state_n = np.vstack([e.n_s for e in experiences])
        est = self.train_network.predict(states)
        fut = self.target_network.predict(state_n)  # Q_{train}(s,a)更新用
        for i, e in enumerate(experiences):
            rew = e.r
            if not e.d:
                rew += max(fut[i]) * self.gamma
            est[i][e.a] = rew
        self.train_network.fit(states, est, epochs=1, verbose=0)

    def play_step(self):
        action = self.get_policy_action(self.state)
        new_state, rew, done, info = self.env.step(action)  # 新しい状態を取得
        if new_state[0] >= 0.4:
            rew += 0.5
        elif new_state[0] >= 0.5:
            rew += 10.0
        e = Exp([self.state], action, rew, [new_state], done)
        self.exp.append(e)
        self.state = new_state
        return done, rew

    def play_episode(self, flag):
        self.state = self.env.reset()  # envを初期化
        ac_score = 0.0

        for tim in range(self.iteration_times):
            if flag:  # Trueなら映像を表示.
                self.env.render()
            done, reward = self.play_step()
            ac_score += reward
            self.update2()
            if done:
                break

        self.target_network.set_weights(self.train_network.get_weights())  # 両方のnetworkの同期.
        self.episode_times += 1
        if tim < 198:
            print("Success!.")
        print("eps is {},episode_times is {},length is {}times, score is {}".format(self.eps, self.episode_times, tim, ac_score))
        self.eps_change()

    def learn2(self):
        for i in range(train_repeat_times):
            flag = False
            if (i % 20) == 0:
                flag = True
            self.play_episode(flag)
        print("end.")


player = Agent()
player.learn2()

結果

結果には多少差があります

1.はじめに

2.FRITとは?

問題設定

FRITの具体的な手続き

なぜこれで上手く行くのか(概略)

3.おわりに

参考文献

はじめに

背景というか定義

Metropolis-Hastings法

Gibbs sampling

MCMCを応用した画像のノイズ除去

原理

結果と実装

雑な感想と近況報告

記号の定義

SACとは?

概要

学習とそれぞれの損失関数

・Critic()の損失関数

・Actor()の損失関数

全体の流れ

environmentについて

実装

感想など

出典,参考など

A2Cってなに?

どういう仕組み?

Actor-Criticとは?

損失関数の構造

実装

結果

嵌った所

参考にした物

ポエム

DDPGとは

概要

決定論的な方策勾配定理

損失関数

soft-target

実装

Pendulum-v0の問題設定

結果

参考文献

この記事上での定義

ベルマン方程式再び

ベルマン作用素とその性質

2つのベルマン作用素

定義

性質

1.単調性

2.ベルマン作用素はバラせる.

3.収束性

4.一意性

5.縮小性

まとめ

MountainCar-v0の問題設定

fixed target Q-networkについて

定義再び

どういう手法か

その他の工夫

Experience Replay

実装

結果

・Critic( $Q_\theta$ )の損失関数

・Actor( $\pi$ )の損失関数