I'm user:15414362. Please merge my old account user:2931409 to it.
This post is for recovery of my SO account.
離散化を行うモデルの検討・実装
1 本記事で扱う「離散化」の定義
まず、語弊がないよう、本記事で扱う「離散化」を定義しておきます。まず、本記事では、数直線上のある領域を何分割かにし、点で区切ることを「離散化」と言うこととします。一般的に、数直線上の領域[X1, X2]をN分割するときは、等分割し、分割した領域のそれぞれの長さは (X2 - X1)/N になります。本記事では、このように離散化で領域をN分割する時、領域をN等分に分割するのではなく、分割後のN個の領域がそれぞれ異なる長さを持つように分割する方法について記載します。領域を異なる長さの領域に分割することができると、離散化が絡む問題の最適化を行う時に、離散化の方法を含め最適化できるようになるかもしれません。
2 モデルの概要
以下のような数直線を考えます。
図1: 数直線(1)
これを以下のように分割数N=3でd1,d2,d3の領域に分割することを考えます。
図2: 数直線(2)
すると、以下の式が成立します。
これは行列で表すと以下になります。
ここで、行列Tを以下のように定義します。この行列Tは下三角行列となっています。
また、x1,x2,x3をX, d1,d2,d3をDとして以下のように定義します。
すると先ほどのx1,x2,x3を表す式は以下のようになります。
この式の、d1,d2,d3が決まれば、 x1,x2,x3 も決まりますので、d1,d2,d3が学習すべきパラメータになると推測できます。
3 モデルの詳細
モデルの概要 では数直線上の[0,10]をざっくり分割することだけ考えていました。この章では、領域[s(start),e(end)]をN分割することを考えます。
モデルの概要 を参考に新たにモデルを以下のように定義します。
\(\lvert D \rvert\) は \(d1 + d2 + d3 + ... + dN\) です。 \(D\) をこれで割り算することで[0, 1]の範囲の値をとるように正規化しています。それを \(e -s\) でスケール、 \(s\) でシフトして、[s, e]の値をとるようにしています。\(T\) は 下三角行列です。このモデルを使用して、[s, e] をN分割することを考えます。
例えば、 区間[-5,5] を10分割する場合
この式の \(x10\) を計算すると 5 となります。これは [s, e] の e ですので、以下になります。
e を加えたので s も加えてみましょう。
これは、数直線で表すと以下になります。今回の例では \(s = -5\) 、 \(e = 5\) です。
図3: 数直線(3)
式をみるとパラメータ \(d10\) が分数の分子にはなく、分母に現れているのが興味深いです。 \(x9\) の式は右辺に \(\frac{d1 + d2 + d3 + d4 + d5 + d6 + d7 + d8 + d9}{d1 + d2 + d3 + d4 + d5 + d6 + d7 + d8 + d9 + d10}\) を含みます。この式は分母に \(d10\) があるおかげで 1 にならず, \(d10\) のための隙間を残しているように見えます。
結局のところ、左辺が s, e となっている式は学習と関係ありません。それを踏まえると、最終的に [s, e] を N分割する時に、学習に使用する式は以下になります。
試しに、 N=2 のときこの式がどうなるか見てみましょう。
これは、区間[s, e]をd1:d2に内分する数学の公式と同じになっています。そのため、問題ないことがわかります。(もしかすると本モデルは内分の公式を拡張したものになっているのかもしれません)
4 モデルまとめ
[s, e]をN分割する作成したモデルをもう一度以下に記載します。
これは、以下のように領域[s,e]をN分割します。
図4: パラメータと数直線の関係
5 使ってみる
このモデルを使用して、[-5, 5] の区間を N=10 で均等に分割することを考えます。 d1, d2, …, d10 はランダムに初期化し、均等に分割した場合との差分(MSE)を損失関数として定義し、それを最小化します。ソースコードは以下のようになります。
5.1 ソースコード
import numpy as np import torch N = 10 s, e = -5, 5 # モデル D = torch.rand((N,), requires_grad=True) T = torch.tensor(np.tri(N-1).astype(np.float32)) # 教師データ t = torch.tensor(np.linspace(s, e, N+1, dtype=np.float32)).view(-1, 1)[1:N] # 学習準備 optimizer = torch.optim.Adam([D]) lossfn = torch.nn.MSELoss() # 学習 for i in range(10000): optimizer.zero_grad() X = (e - s) * torch.matmul(T, D.view(-1, 1)[:N-1])/D.sum() + s if i==0: print("Before Training X:") print(X) loss = lossfn(X, t) loss.backward() optimizer.step() print("After Training X:") print(X)
6 考察
7 まとめ
離散化を行うモデルの検討・実装を行いました。
おやすみなさい。
gymモジュール Pendulum-v0 で、振り子を頂上で止まらせる
今回は gym モジュールの Pendulum-v0 の環境において、強化学習により振り子を頂上でとまらせます。
1 前回の記事
前回の記事は以下の通りです。本記事を読むにあたり参考になるかもしれません。
1.1 gymモジュール Pendulum-v0 で振り子を等速度で回転させる
1.2 gymモジュール Pendulum-v0 で振り子を揺らす
1.3 gymモジュール Pendulum-v0 で振り子をトルクを与えて静止させる
2 勉強会で発表したスライド
勉強会でも発表しています。資料は以下におきました。
3 今回行うこと
強化学習により振り子を頂上でとまらせます。
4 解決方法
4.1 報酬の定義
ここでは、報酬を \(\cos{\theta}\) のように定義します。頂上に近ければ高い報酬を、遠ければ低い報酬を与えて強化学習を行えば頂上で止まるようにしています。また、運動エネルギーが少ないほど高い報酬を与えるよう、\(\dot \theta^2\)で割り算しています。理論上は \(\cos{\theta}\) でも良いと思いますが、何度か試したところ、\(\cos{\theta}\)を３乗し、さらに運動エネルギーを考慮した報酬の式の方が頂上でとまることが多かったです。これが効率良い報酬の式だとは思いませんし、物理学的に頂上でとまるのに適した報酬の式があると思います。
| 報酬の定義 |
|---|
| \(r=\frac{\cos{\theta}^3}{\dot{\theta}^2}\) |
| 振り子の位置 | \(\theta\) | \(\cos{\theta}\) |
|---|---|---|
| 一番上 | \(0^\circ\) | 1 |
| 中間 | \(90^\circ\) | 0 |
| 一番下 | \(180^\circ\) | -1 |
4.2 状態遷移確率、方策の表の作成
状態遷移確率、方策の表の作成について説明します。
4.2.1 表の作成
- 理論
状態遷移確率を求めます。振り子の状態\(s\)に行動としてトルク τ を与えると状態が\(s'\)になります。今回は以下のような表形式で、状態\(s\)の時にトルク τ を与えた時の次の状態\(s'\)と\(r\)(報酬)を表すこととします。
表3: 状態、トルク、次の状態、報酬 \(s\) \(a\) \(r\) \(s'\) (\(\theta\), \(\dot \theta\)) τ \(\frac{\cos{\theta}^3}{\dot{\theta}^2}\) (\(\theta'\), \(\dot \theta'\)) パラメータのとりうる値の範囲は以下になります。
表4: パラメータのとりうる値の範囲 \(\theta\) \(\dot \theta\) τ \(r(=\frac{\cos{\theta}^3}{\dot{\theta}^2})\) \(\theta'\) \(\dot \theta'\) \([0^\circ,360^\circ)\) \([-8, 8]\) \([-2, 2]\) \([-\infty, \infty]\) \([0^\circ, 360^\circ°)\) \([-8, 8]\) \(\dot \theta\) や τ のパラメータの上限、下限値は、 gym モジュールの Pendulum-v0 の環境でデフォルトで設定されている値です。これらのパラメータのとりうる値の範囲の詳細については以下を参照してください。
報酬 \(r(=\frac{\cos{\theta}^3}{\dot{\theta}^2})\)については、\(\dot \theta\) が 0 を取りうるため、計算すると \(\infty\) になる可能性があります。実際のプログラムでは \(\infty\) にならないように微小な値を分母に加算しています。
上記をもとに\(s\) , \(a\) , \(r\) , \(s'\)の表(状態遷移確率)を作成します。パラメータの値は連続の値ですが、今回はこれを扱い易くするため離散化します。
振り子の運動は物理的な現象ですので初期状態が同じなら必ず同じ結果になります。(\(\theta\), \(\dot \theta\), τ ) で行動すれば確率100%で一つの状態(\(\theta'\), \(\dot \theta'\))に遷移します。しかし、今回は離散化するため、そうならない時もありますので注意が必要です。
- 実装
以下のコードにより、状態遷移確率、方策の表を作成します。
- コード
import numpy as np import pandas as pd ############################ # パラメータの表を作成する # ############################ from itertools import product class N(): theta = 2 theta_dot = 2 tau = 2 actions = tau ################################## # 表に使用するインデックスの定義 # ################################## # theta, tau, theta_dot の範囲を定義 bin_theta = np.linspace(0, 360, N.theta + 1) bin_theta_dot = np.linspace(-8, 8, N.theta_dot + 1) bin_tau = np.linspace(-2, 2, N.tau + 1) # 角度、速度、トルク を離散化(pd.IntervalIndex.from_breaks を使用) theta_index = pd.IntervalIndex.from_breaks(bin_theta) theta_dot_index = pd.IntervalIndex.from_breaks(bin_theta_dot) tau_index = pd.IntervalIndex.from_breaks(bin_tau) # インデックスの作成 index_0 = pd.MultiIndex.from_product( [theta_index, theta_dot_index, tau_index], names=["theta_cat", "theta_dot_cat", "tau_cat"]) index_1 = pd.MultiIndex.from_product([theta_index, theta_dot_index], names=["theta_2_cat", "theta_dot_2_cat"]) ############ # P の作成 # ############ # 「角度、速度、トルク」 - 「次の角度、次の速度」の表を作る P = pd.DataFrame(0, index=index_0, columns=index_1) ################## # 方策の表を作成 # ################## policy = pd.Series(1 / N.actions, index=index_0) ######## # 表示 # ######## print("状態遷移確率の表") print(P) print("方策の表") print(policy)
- コード
- 表のサンプル
コードを実行すると以下のような状態遷移確率の表(P) と 方策の表(policy) が得られます。実際に学習に使用する表は大きすぎるため、サイズを小さくして表示しています。
- 状態遷移確率の表(P)
表5: 状態遷移確率(P) \(\theta'\) (0.0, 180.0] (180.0, 360.0] \(\dot \theta'\) (-8.0, 0.0] (0.0, 8.0] (-8.0, 0.0] (0.0, 8.0] \(\theta\) \(\dot \theta\) τ (0.0, 180.0] (-8.0, 0.0] (-2.0, 0.0] 0 0 0 0 (0.0, 2.0] 0 0 0 0 (0.0, 8.0] (-2.0, 0.0] 0 0 0 0 (0.0, 2.0] 0 0 0 0 (180.0, 360.0] (-8.0, 0.0] (-2.0, 0.0] 0 0 0 0 (0.0, 2.0] 0 0 0 0 (0.0, 8.0] (-2.0, 0.0] 0 0 0 0 (0.0, 2.0] 0 0 0 0 - 方策の表(policy)
表6: 方策の表(policy) \(\theta\) \(\dot \theta\) τ Prob (0.0, 180.0] (-8.0, 0.0] (-2.0, 0.0] 0.5 (0.0, 2.0] 0.5 (0.0, 8.0] (-2.0, 0.0] 0.5 (0.0, 2.0] 0.5 (180.0, 360.0] (-8.0, 0.0] (-2.0, 0.0] 0.5 (0.0, 2.0] 0.5 (0.0, 8.0] (-2.0, 0.0] 0.5 (0.0, 2.0] 0.5 状態遷移確率の実際の値についてはシミュレータを動かすまでわからないのですが、ここでは、その値を求めるのは後回しにし、ひとまず 0 にしています。方策の表については、どの方策が最適化わからないので、全ての行動をまんべんなく行うよう(上の表で Prob の値を 0.5 にしています)にしています。このような方策の表の値を強化学習のアルゴリズムにより報酬が最大になるように調整します。
- 状態遷移確率の表(P)
4.2.2 状態遷移確率の実際の値を求める
- 考え方
(\(\theta\), \(\dot \theta\), τ ) のそれぞれの要素は実数です。今回は、これを扱い易くするため、離散化します。以下のような、ビンの数で分割して、それぞれのビンの中間の値を使用することとします。全部で25*50*25通りのパターンになります。一つ一つのパターンについて、シミュレータに(\(\theta\), \(\dot \theta\), τ )を与え、(\(\theta'\), \(\dot \theta'\))を得ます。このときの、(\(\theta\), \(\dot \theta\), τ ) と (\(\theta'\), \(\dot \theta'\)) の値を、 状態遷移確率Pの表にキーとして与え、そのバリューを +1 します。最後に、 P の各行を25(行動数) で割り、確率にすれば、実際のPが得られます。
表7: (\(\theta\), \(\dot \theta\), τ ) の分割数 分割数 \(\theta\) 25 \(\dot \theta\) 50 τ 25 - コード
以下のコードにより、状態遷移確率(P) の実際の値を求めています。
######################## # シミュレータを動かす # ######################## ############### # gymの初期化 # ############### # 環境の作成 env = gym.make('Pendulum-v0') # 環境の初期化(明示的に行う必要がある) env.reset() for theta_cat, theta_dot_cat, tau_cat in Psas.index: # 振り子の theta, speed を設定 # ビンの中間の値を取得するため Interval オブジェクトの mid プロパティを使用 env.env.state = np.deg2rad(theta_cat.mid), theta_dot_cat.mid obser, r, done, info = env.step([tau_cat.mid]) theta_2, theta_dot_2 = env.env.state # デフォルトで rad なので degree に変換 theta_2 = np.rad2deg(theta_2) theta_2 = f(theta_2) # カウントを +1 する Psas.loc[(theta_cat, theta_dot_cat, tau_cat), (theta_2, theta_dot_2)] += 1 # Psas 確率に変換 Psas = Psas / Psas.sum(axis=1).values.reshape(-1,1)
4.2.3 方策の初期化
以下のコードにより、方策(PI) を初期化しています。
# 方策の初期化 PIas = pd.DataFrame(1/N.actions, index=index_0, columns=["PI_tau_Prob"])
4.3 学習
学習について説明します。
4.3.1 理論
学習はベルマン方程式を解析的に解くことにより行います。詳細な説明は 参考文献 を参考にしてください。ざっくりと説明すると、以下の計算を行うことで状態価値Vと行動価値Qを求め、それを用いて方策を更新することで学習を行います。
(1)(2)(3)(4)は 参考文献 現場で使える！Python深層強化学習入門 強化学習と深層学習による探索と制御 P36, P32 より引用
4.3.2 コード
- 状態遷移確率(P) と方策の確率を合わせたものを求める
以下のコードにより、状態遷移確率(P) と方策の確率を合わせたものを求めています。状態遷移確率と方策の行列の各要素の積をとり、行動 τ について和をとっています。
# 状態遷移確率と方策の確率を合わせたものを求める P_PI_ss = Psas * PIas.values.reshape(-1, 1) P_PI_ss = P_PI_ss.sum(level=[0, 1])
- 状態遷移確率と方策から報酬の期待値を求める
以下のコードにより、状態遷移確率と方策から報酬の期待値を求めます。状態遷移確率と方策と報酬の行列の各要素の積をとり、行動 τ について和をとっています。
# R_PI_s を求める R_PI_s = Psas * PIas.values.reshape(-1, 1) * Rsas R_PI_s = R_PI_s.sum(axis=1).sum(level=[0, 1])
- ベルマン方程式を解き状態価値を求める
以下のコードにより、ベルマン方程式を解き状態価値を求めます。ベルマン方程式を解くのは以下の行列計算で完了します。計算は逆行列の計算を含みます。
# 状態価値 V を求める(ベルマン方程式を解析的に解く) size = len(R_PI_s) V_values = np.linalg.inv(np.eye(size) - gamma*P_PI_ss.values).dot(R_PI_s) V_s = pd.Series(V_values, index=R_PI_s.index)
- 行動価値を求める
以下のコードにより、行動価値を求めます。行動価値は即時報酬、状態価値を使用して計算できます。
# 行動価値 Q を求める Q_sa = Psas*(Rsas + gamma*V_s.values.reshape(1,-1)) Q_sa = Q_sa.sum(axis=1)
- 方策を更新する
以下のコードにより、方策を更新します。ここでは、行動価値が最大になる行動を確率1でとるように方策を更新します。ここまでの 状態遷移確率(P) と方策の確率を合わせたものを求める 〜 方策を更新する までが学習に必要な１ステップです。
# 方策を更新 # 状態のグループから行動価値が最大となるインデックスを求める # https://stackoverflow.com/questions/27842613/pandas-groupby-sort-within-groups largest_Qsa = Q_sa.groupby(level=[0,1], group_keys=False).nlargest(1) PIas.loc[:] = 0 for (theta_cat, theta_dot_cat), df in largest_Qsa.groupby(level=[0, 1]): PIas.loc[df.index] = 1
- 方策が更新されなくなるか、指定回数を超えたら方策の更新を打ち切る
方策が更新されなくなるか、指定回数を超えたら方策の更新を打ち切ります。方策を更新したのに、前の方策と同じになった場合は、その次も同じ、その次もと永遠と同じ方策になるため更新を打ち切ります。今回は、指定回数を超えても終了するようにプログラムしています。
# 方策が更新されなくなったら終了 if (pre_PIas.values==PIas.values).all(): break # 上限回数を超えても終了する cnt += 1 if cnt > 10: break
5 コード全体
作成したコードの全体がこちらです。
import numpy as np import pandas as pd import gym import time # 乱数のシードを設定 seed = 100 np.random.seed(seed) # 離散化に使用する分割数を定義 class N(): theta = 25 theta_dot = 50 tau = 25 actions = tau # 角度を [0-360) にする関数を定義 f = lambda x: x - np.floor(x/360)*360 ###################################### # 状態遷移確率の表の実際の値を求める # ###################################### ################## # ビンを作成する # ################## # theta, tau, theta_dot の範囲を定義 bin_theta = np.linspace(0, 360, N.theta + 1) bin_theta_dot = np.linspace(-8.01, 8.01, N.theta_dot + 1) bin_tau = np.linspace(-2.01, 2.01, N.tau + 1) # 角度、速度、トルク を離散化(pd.IntervalIndex.from_breaks を使用) theta_index = pd.IntervalIndex.from_breaks(bin_theta) theta_dot_index = pd.IntervalIndex.from_breaks(bin_theta_dot) tau_index = pd.IntervalIndex.from_breaks(bin_tau) # インデックスの作成 index_0 = pd.MultiIndex.from_product( [theta_index, theta_dot_index, tau_index], names=["theta_cat", "theta_dot_cat", "tau_cat"]) index_1 = pd.MultiIndex.from_product([theta_index, theta_dot_index], names=["theta_2_cat", "theta_dot_2_cat"]) print("Psasを作成") # 状態遷移確率用 Psas = pd.DataFrame(0, index=index_0, columns=index_1) print("Rsasを作成") # 報酬は theta が 0°に近いほどよいかつ、運動エネルギーが小さいほどよいようにした。 data = [] for theta_2_cat, theta_dot_2_cat in Psas.columns: data.append(np.cos(np.deg2rad(theta_2_cat.mid))**3/(theta_dot_2_cat.mid**2 + 0.00001)) arr = np.zeros((len(Psas),1)) + np.array(data) Rsas = pd.DataFrame(arr, index=index_0, columns=index_1) ######################## # シミュレータを動かす # ######################## ############### # gymの初期化 # ############### # 環境の作成 env = gym.make('Pendulum-v0') # 環境の初期化(明示的に行う必要がある) env.reset() for theta_cat, theta_dot_cat, tau_cat in Psas.index: # 振り子の theta, speed を設定 # ビンの中間の値を取得するため Interval オブジェクトの mid プロパティを使用 env.env.state = np.deg2rad(theta_cat.mid), theta_dot_cat.mid obser, r, done, info = env.step([tau_cat.mid]) theta_2, theta_dot_2 = env.env.state # デフォルトで rad なので degree に変換 theta_2 = np.rad2deg(theta_2) theta_2 = f(theta_2) # カウントを +1 する Psas.loc[(theta_cat, theta_dot_cat, tau_cat), (theta_2, theta_dot_2)] += 1 # Psas 確率に変換 Psas = Psas / Psas.sum(axis=1).values.reshape(-1,1) # 方策の初期化 PIas = pd.DataFrame(1/N.actions, index=index_0, columns=["PI_tau_Prob"]) cnt = 0 gamma = 0.99 while True: print(cnt) # 状態遷移確率と方策の確率を合わせたものを求める P_PI_ss = Psas * PIas.values.reshape(-1, 1) P_PI_ss = P_PI_ss.sum(level=[0, 1]) # R_PI_s を求める R_PI_s = Psas * PIas.values.reshape(-1, 1) * Rsas R_PI_s = R_PI_s.sum(axis=1).sum(level=[0, 1]) # pre_V_s = V_s.copy() pre_PIas = PIas.copy() # 状態価値 V を求める(ベルマン方程式を解析的に解く) size = len(R_PI_s) V_values = np.linalg.inv(np.eye(size) - gamma*P_PI_ss.values).dot(R_PI_s) V_s = pd.Series(V_values, index=R_PI_s.index) # 行動価値 Q を求める Q_sa = Psas*(Rsas + gamma*V_s.values.reshape(1,-1)) Q_sa = Q_sa.sum(axis=1) # 方策を更新 # 状態のグループから行動価値が最大となるインデックスを求める # https://stackoverflow.com/questions/27842613/pandas-groupby-sort-within-groups largest_Qsa = Q_sa.groupby(level=[0,1], group_keys=False).nlargest(1) PIas.loc[:] = 0 for (theta_cat, theta_dot_cat), df in largest_Qsa.groupby(level=[0, 1]): PIas.loc[df.index] = 1 # 状態価値の最小値を表示 print(V_s.min()) # 方策が更新されなくなったら終了 if (pre_PIas.values==PIas.values).all(): break # 上限回数を超えても終了する cnt += 1 if cnt > 10: break first = False ########################################## # 求めた方策に従い、シミュレータを動かす # ########################################## # 環境の作成 env = gym.make('Pendulum-v0') env.seed(1) for __ in range(10): # 環境の初期化(明示的に行う必要がある) env.reset() for _ in range(400): # 現在の状態を取得 theta, theta_dot = env.env.state # デフォルトで rad なので degree に変換 theta = np.rad2deg(theta) theta = f(theta) tau = PIas.loc[theta].loc[theta_dot].idxmax().values[0].mid env.step([tau]) env.render() # time.sleep(0.03) env.close()
6 出力結果
出力結果は以下になります。
図1. が コード全体 実行結果ですが、頂上で止まっていません。ここから乱数のシードやコードを変えて試行錯誤を繰り返し手動でチューニングを行い学習させたところ、それらしい動きになりました。
7 考察
振り子の動きを見ると、左右に振って徐々に上昇させていることがわかります。これは、一方向に τ の上限のトルクを加えても、頂上まで上昇することがあまりないため、上昇させるには左右に振るのが定石だからです。(頂上付近だったり、初速度が大きい場合は一方向のトルクでも頂上にたどり着きます)また、左右に振るだけでも、振り子が一番下にいる確率を下げれますので、これだけで、マイナスの報酬を少なくできます。少なくともこの左右に振る動きを自動的に獲得することができました。ただし、頂上で静止するのは難しいようです。振り子が垂直の時、重力と釣り合うため、他の状態に比べ、少しのトルクで大きく動きます。そこに離散化による誤差の影響が加わり学習が難しくなっているかもしれません。ただ、上記図2.のように、パラメータによっては、うまく学習するときもあるようです。離散化の時に\(\dot \theta\)を他の分割数を他のパラメータの分割数より多くしているのは、他の分割数をあげるよりより、\(\dot \theta\)の分割数をあげる方が良い動きが見られたからです。離散化の方法も学習に大きな影響を及ぼすことがわかりました。出力結果をみると、真の最適な方策に収束しているようには思いませんが、それに近い方策になっていると思います。簡単な報酬を定義しただけで、アルゴリズムが自動的に最適に近い方策を出してくれるので便利です。
8 まとめ
Pendulum-v0 について、報酬を独自に定義し、連続値を離散化し状態遷移確率などを表で表し、ベルマン方程式を解くことにより、方策の学習を行いました。学習済のパラメータを用いて、シミュレーションを10回繰り返したところ、10回全てで振り子を頂上で止まらせることができました。(図２のアニメーションを参照)
gymモジュール Pendulum-v0 で振り子をトルクを与えて静止させる
前回に続き、gym モジュールの Pendulum-v0 の理解を深めます。今回は、手で振り子を水平状態にし、その状態をキープするようにトルク与えて、振り子を静止させるシミュレーションを行います。今回も、強化学習は行いません。シミュレータの理解を深めます。
Table of Contents
1 前回の記事
前々回、前回の記事は以下の通りです。本記事を読むにあたり参考になるかもしれません。
1.1 gymモジュール Pendulum-v0 で振り子を等速度で回転させる
1.2 gymモジュール Pendulum-v0 で振り子を揺らす
2 今回行うこと
今回は、手で振り子を水平状態にし、その状態をキープするようにトルク与えて、振り子を静止させるシミュレーションを行います。この日本語を、プログラム的に表すと以下になります。
| No. | 日本語 | プログラム | 説明 |
|---|---|---|---|
| 1 | 手で振り子を水平状態にする | initial_theta = pi/2 | 水平状態に相当する角度は pi/2 のため、振り子の初期の角度を pi/2にします。手をそっと離すため、初期の速度は0にします。 |
| 2 | その状態をキープするようにトルク与えてシミュレーションする | env.step([-5]) | 手を離したあと、物理法則に従って振り子が動いて欲しいため env.step を使用します。引数の値はトルクを表し、引数を[-5]とすることで振り子を静止させます。([-5]にする理由は後述) |
表1. 今回行うこと
手を離したときに振り子の状態をキープするために重力が振り子に与える力と同じ力を与える必要があります。結論を言ってしまうとトルクとして-5を与えれば良いのですが、これをどう導き出したのか以下に記載します。
最初に、使用する物理的な変数を以下の表にまとめます。
| No. | 変数名 | 説明 |
|---|---|---|
| 1 | m | 質量 |
| 2 | g | 重力 |
| 3 | l | 棒の長さ |
表2. 物理的な変数の説明
まずは、水平状態で振り子にかかる力を求めます。振り子が水平の時、垂直方向にかかる力は重力しかないので、振り子にかかる力以下になります。
F = m * g
このFと同じ力を、重力と反対方向にかければ力が釣り合い、静止します。これは、振り子に重力と反対方向に g の加速度を与えてあげれば良いということです。
gym の Pendulum の環境では振り子に与える力を、支点からの距離と力の積であるトルクで指定する必要があるため、上記 F と釣りあうトルクを求めます。
トルクを求めるには、回転の中心である支点からの距離lに、その点にかかっている力 F を掛け算します。
T = F * l
回転体が糸と球でできているなら、この式で正しいですが、シミュレータをよく見みると、回転体は棒でできていることがわかります。質量の全てが棒の先端にあるわけではないため、棒の回転に必要なトルクは、そうである場合に比べて少なくてすむはずです。この場合、回転体の重心の位置(支点からl/2 離れた点)に質量mの球があるという問題に置き換えて考えます。以上より、今回の振り子(振り棒?)の重力を支えるために必要なトルクは以下になります。
T = F * l / 2
T = m * g * l / 2
最後に、実際にプログラムで指定すべきトルクの値を求めるため、gymモジュール の Pendulum で定義されている m, g, l の実際の値を調べます。調べたところそれぞれのデフォルトの値は以下になっていました。
| No. | 変数 | デフォルト値 |
|---|---|---|
| 1 | m | 1 |
| 2 | g | 10 |
| 3 | l | 1 |
表3. gymモジュール の Pendulum で定義されている m, g, l の実際の値
よって、今回与えるべきトルク T は以下になります。(重力の反対方向に力を与えるため、マイナスにしています)
T = -1 * 10 * 1 / 2
T = -5
注意点としては、Pendulum ではトルクの上限がデフォルトで 2 までとなっています。今回は、トルクとして 5 を与える必要があるので、この上限を引きあげる必要があります。上限の設定は max_torque 変数を使用します。
3 コード
今回書いたコードは以下になります。
import gym import time from math import pi env = gym.make('Pendulum-v0') # 強化学習の文脈における「環境」の作成 env.reset() # 環境の初期化(明示的に行う必要がある) env.env.max_torque = 5 # トルクの上限を引き上げる <---!前回との差分はこの行 # 初期状態 initial_theta = pi/2 initial_speed = 0 theta = initial_theta speed = initial_speed env.env.state = theta, speed # 振り子の theta, speed を設定 env.render('human') # 描画 time.sleep(0.033) # 待つ、 0.033 は 30fps にするため for i in range(300): # 30fps で 10秒間表示 env.step([-5]) # トルク -5 で振り子を次の位置に更新する <---!前回との差分はこの行 env.render('human') # 描画 time.sleep(0.033) # 待つ、 0.033 は 30fps にするため env.close()
4 コードの説明
env.step([x]) で、トルクx を振り子に与えて、シミュレーションを 1ステップ進めています。今回は env.step([-5])で、トルク-5を与えて、シミュレーションを1ステップ進めています。そして、 env.env.max_torque = 5 でトルクの上限を引き揚げています。前回との差分は以下の通りです。
| No. | 行数 | 前回 | 今回 | 説明 |
|---|---|---|---|---|
| 1 | 7 | - | env.env.max_torque = 5 | トルクの上限を5に引き揚げる |
| 2 | 18 | env.step([0]) | env.step([-5]) | トルク-5を与え、シミュレーションを１ステップ進める |
表4. 前回のコードとの差分
5 出力結果
出力結果は以下になります。
水平状態で重力と釣り合うようなトルクを与えているため、動かないことがわかります。
ただ、動かないと、本当にプログラムが正しく動作しているかわからないため、トルクをわずかに減らして T =-4.9 として実行して見ます。振り子がフラフラし、プログラムが動作していることがわかります。
6 まとめ
gymモジュールの Pendulum-v0 の環境を使用し、 手で振り子を水平状態にし、その状態をキープするようにトルク与えて、振り子を静止させるシミュレーションを行いました。
以上です。
gymモジュール Pendulum-v0 で振り子を揺らす
前回に続き、gym モジュールの Pendulum-v0 の理解を深めます。今回は、 手で振り子を水平状態にし、離すことに相当するシミュレーションを行います。今回も、強化学習は行いません。シミュレータの理解を深めます。
1 前回の記事
前回の記事はこちらです。gymとは何か、pendulumとは何かについてはこちらをご覧ください。
2 今回行うこと
今回は、 手で振り子を水平状態にし、離すことに相当するシミュレーションを行います。この日本語を、プログラム的に表すと以下になります。
| No | 日本語 | プログラム | 説明 |
|---|---|---|---|
| 1 | 手で振り子を水平状態にする | initial_theta = pi/2 | 水平状態に相当する角度は pi/2 のため、振り子の初期の角度を pi/2にします |
| 2 | 振り子から手を離す | initial_speed = 0 | 手をそっと離すため、速度が0のため、振り子の初速度を 0 にします |
| 3 | シミュレーションする | env.step([0]) | 手を離したあと、物理法則に従って振り子が動いて欲しいため env.step を使用します。引数の値はトルクを表し、今回は[0]とします。(トルクを与えると手を離すだけだけでなくなってしまいます |
表1. 今回行うこと
3 コード
今回書いたコードは以下になります。
import gym import time from math import pi env = gym.make('Pendulum-v0') # 強化学習の文脈における「環境」の作成 env.reset() # 環境の初期化(明示的に行う必要がある) # 初期状態 initial_theta = pi/2 initial_speed = 0 theta = initial_theta speed = initial_speed env.env.state = theta, speed # 振り子の theta, speed を設定 env.render('human') # 描画 time.sleep(0.033) # 待つ、 0.033 は 30fps にするため for i in range(300): # 30fps で 10秒間表示 env.step([0]) # トルク 0 で振り子を次の位置に更新する env.render('human') # 描画 time.sleep(0.033) # 待つ、 0.033 は 30fps にするため env.close()
4 コードの説明
env.step([x]) で、トルクx を振り子に与えて、シミュレーションを 1ステップ進めることができます。今回は env.step([0])で、トルク0、すなわちトルクを与えず、シミュレーションを1ステップ進めています。
5 出力結果
こちらからトルクは与えませんが、重力はデフォルトでかかっている状態のため、シミュレーション結果をみると現実世界の振り子のような動きになることがわかります。減衰する様子は見られないため、摩擦はないことがわかります。
6 まとめ
gymモジュールの Pendulum-v0 の環境を使用し、 振り子を水平状態にし、離すことに相当するシミュレーションを行いました。現実世界の振り子のような動きになることがわかりました。
以上です。
gymモジュール Pendulum-v0 で振り子を等速度で回転させる
gym モジュールの Pendulum-v0 とは以下のような振り子のシミュレータです。
今回は、このシュミレータを用いて、振り子を等速度で回転させるプログラムを書きます。振り子をアナログ時計の秒針のように一定の速度でグルグル回すことを目的とします。このプログラムを作ることで gymモジュール の API に触れ、gymの理解を深めます。振り子をグルグル回してAPI に触れるだけです。強化学習は使いません。もう一度言いますが、強化学習で振り子を等速度に制御する話ではありません。
Table of Contents
1 そもそも gym とは
gym とは公式サイトによると以下のように記載されています。
| 公式サイト(https://gym.openai.com/)の説明 | 翻訳 |
|---|---|
| Gym is a toolkit for developing and comparing reinforcement learning algorithms. It supports teaching agents everything from walking to playing games like Pong or Pinball. | Gym とは reinforcement learning(強化学習) アルゴリズムの開発や比較のためのツールキットです.歩行からPongやピンボールまでエージェントの教育をサポートします。 |
表1. gymとは
強化学習の理論の確認や練習に使えそうです。
2 Pendulum-v0 とは
Pendulum とは日本語では振り子を意味します。さらに v-0 はバージョン0ということです。gym の Pendulum-v0 はシンプルな振り子の環境を表しています。なんとなくv-1, v-2 なども存在しそうなことがわかります。
3 コード
今回書いたコードは以下になります。 １秒間に 45° だけ回転します。 それを8回繰り返して１周させています。
import gym import time from math import pi env = gym.make('Pendulum-v0') # 強化学習の文脈における「環境」の作成 env.reset() # 環境の初期化(明示的に行う必要がある) # 初期状態 theta = 0 env.env.state = theta, 0 # 振り子の theta, speed を設定 env.render('human') # 描画 time.sleep(1) # 1秒待つ for i in range(8): theta += pi/4 # theta env.env.state = theta, 0 # 振り子の theta, speed を設定 env.render('human') # 描画 time.sleep(1) # 1秒待つ env.close()
コード1. gymモジュール Pendulum-v0 で振り子を等速度で回転させる
4 コードの説明
env.env.state がポイントで、調べたところ、この変数が振り子の角度、振り子の速さの二つの情報を保持しているようでした。今回のプログラムでは、この変数に、振り子の角度を代入することで、振り子の角度を設定しています。最後に env.render() を実行し振り子のビジュアルな状態をディスプレイ上に書き出します。
5 出力結果
以下のように、振り子が１秒間に 45° だけ回転し、１週します。
6 まとめ
gymモジュールの Pendulum-v0 の環境を使用し、 振り子を等速度で回転させました。(厳密にいうと等速度に回転しているように見えるようにプログラムした)env.env.state で振り子の位置(角度)を指定できることがわかりました。gymモジュールの API をほんの少し理解することができました。
以上です、おやすみなさい。
混合分布が確率であることについて
混合分布が確率であることについて書きます。混合分布とは、確率質量関数または確率密度関数を線形結合させたものだと理解しています。本を読んでいたところ、結合係数の和が1のとき、混合分布も確率質量関数または確率密度関数として解釈することができるということでした。個人的にこれがなぜかわからなかったので備忘録としてまとめました。
Table of Contents
1 わからなかったこと
まず、わからなかったことについて書きます。「このように考えたらわかるようになった」ということについては次章に書きます。
話を簡単にするため、ここでは、離散的な事象についての確率質量関数に限定して話を進めます。例えば、四面体のサイコロと、いかさまサイコロの二つの確率質量関数の分布が以下だったとします。※サイコロがいかさまか否かは本題ではありません。
普通のサイコロ:
| X | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| P1(X) | 0.25 | 0.25 | 0.25 | 0.25 |
いかさまサイコロ:
| X | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| P2(X) | 0.1 | 0.2 | 0.3 | 0.4 |
このように、離散的な確率分布であれば、イメージしやすいです。まずは、二つの確率分布がイメージできました。これら二つの確率質量関数をそれぞれ P1(X), P2(X) とします。
※Xは確率変数
この二つの確率分布を用いて混合分布の関数 P3 を作成すると以下になります。
P3(X) = w1*P1(X) + w2*P2(X)
※w1, w2 はスカラで、結合係数
| X | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| P3(X) | w1*0.25 + w2*0.1 | w1*0.25 + w2*0.2 | w1*0.25 + w2*0.3 | w1*0.25 + w2*0.4 |
このとき
w1 + w2 = 1
であれば
P3(1) + P3(2) + P3(3) + P3(4)
が成り立ち、P3 は確率質量関数となる。
上記のがなぜそう言えるのかがわかりませんでした
2 このように考えたらわかるようになった
P3(1) + P3(2) + P3(3) + P3(4)
を計算すると確かに 1 になるので確率なのですが、次のように式を変形し
| X | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| P3(X) | w1*0.25 | w1*0.25 | w1*0.25 | w1*0.25 |
| + | + | + | + | |
| w2*0.1 | w2*0.2 | w2*0.3 | w2*0.4 |
このように分解して見ると理解しやすかったです。
そして、
w1*0.25 + w1*0.25 + w1*0.25 + w1*0.25
は
w1*(0.25 + 0.25 + 0.25 + 0.25)
すなわち
w1*1.0
すなわち
w1
結局のところ、「足して1になるベクトルの要素のそれぞれに w1 をかけて、合計を求める」ことは「1にw1をかける」ことと同じということです。
Python でこれを表現すると以下のようになるでしょうか、、、
import numpy as np from sympy import symbols, solve, Eq, simplify w1 = symbols("w1") a1, a2, a3, a4 = symbols("a1 a2 a3 a4") # [0.25, 0.25, 0.25, 0.25] を [a1, a2, a3, a4] としています。 vec = np.array([a1, a2, a3, a4]) s = sum(w1*vec) print(s.factor()) # a1 + a2 + a3 + a4 = 1 の条件を s に反映して s を求める left = a1 + a2 + a3 + a4 a1_new = solve(Eq(left, 1), a1)[0] # s に 条件を反映 print(simplify(s.subs({a1:a1_new})))
w1*(a1 + a2 + a3 + a4) w1
同様に
w2*0.1 + w2*0.2 + w2*0.3 + w2*0.4
は
w2(0.1 + 0.2 + 0.3 + 0.4)
すなわち
w2*1.0
すなわち
w2
です。
これらをまとめると
P3(1) + P3(2) + P3(3) + P3(4)
は
w1 + w2
と等しくなり
w1 + w2 = 1
であるならば
P3(1) + P3(2) + P3(3) + P3(4) = w1 + w2 = 1
となるため、P1(X) と P2(X) を結合係数 w1, w2(ただし w1 + w2 = 1 のとき) で線形結合させた関数 P3(X)は確率質量関数としてみなせる。
以上です。おやすみなさい。
