Kei Minagawa's Blog

皆川圭(@keimina)のブログ、Pythonで試したことを書いていきます

gymモジュール Pendulum-v0 で、振り子を頂上で止まらせる

gymモジュール Pendulum-v0 で、振り子を頂上で止まらせる

今回は gym モジュールの Pendulum-v0 の環境において、強化学習により振り子を頂上でとまらせます。

1 前回の記事

前回の記事は以下の通りです。本記事を読むにあたり参考になるかもしれません。

1.1 gymモジュール Pendulum-v0 で振り子を等速度で回転させる

1.2 gymモジュール Pendulum-v0 で振り子を揺らす

1.3 gymモジュール Pendulum-v0 で振り子をトルクを与えて静止させる

2 勉強会で発表したスライド

勉強会でも発表しています。資料は以下におきました。

speakerdeck.com

3 今回行うこと

強化学習により振り子を頂上でとまらせます。

4 解決方法

今回は、強化学習の枠組みで、振り子が「頂上で止まる」ようにします。本章では、実装したコードの部分的な説明も行います。コード全体は コード全体 を参照してください。

4.1 報酬の定義

ここでは、報酬を \(\cos{\theta}\) のように定義します。頂上に近ければ高い報酬を、遠ければ低い報酬を与えて強化学習を行えば頂上で止まるようにしています。また、運動エネルギーが少ないほど高い報酬を与えるよう、\(\dot \theta^2\)で割り算しています。理論上は \(\cos{\theta}\) でも良いと思いますが、何度か試したところ、\(\cos{\theta}\)を3乗し、さらに運動エネルギーを考慮した報酬の式の方が頂上でとまることが多かったです。これが効率良い報酬の式だとは思いませんし、物理学的に頂上でとまるのに適した報酬の式があると思います。

表1: 報酬の定義
報酬の定義
\(r=\frac{\cos{\theta}^3}{\dot{\theta}^2}\)


表2: 振り子の位置と角度
振り子の位置 \(\theta\) \(\cos{\theta}\)
一番上 \(0^\circ\) 1
中間 \(90^\circ\) 0
一番下 \(180^\circ\) -1

4.2 状態遷移確率、方策の表の作成

状態遷移確率、方策の表の作成について説明します。

4.2.1 表の作成

  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 の環境でデフォルトで設定されている値です。これらのパラメータのとりうる値の範囲の詳細については以下を参照してください。

    github.com

    報酬 \(r(=\frac{\cos{\theta}^3}{\dot{\theta}^2})\)については、\(\dot \theta\) が 0 を取りうるため、計算すると \(\infty\) になる可能性があります。実際のプログラムでは \(\infty\) にならないように微小な値を分母に加算しています。

    上記をもとに\(s\) , \(a\) , \(r\) , \(s'\)の表(状態遷移確率)を作成します。パラメータの値は連続の値ですが、今回はこれを扱い易くするため離散化します。

    振り子の運動は物理的な現象ですので初期状態が同じなら必ず同じ結果になります。(\(\theta\), \(\dot \theta\), τ ) で行動すれば確率100%で一つの状態(\(\theta'\), \(\dot \theta'\))に遷移します。しかし、今回は離散化するため、そうならない時もありますので注意が必要です。

  2. 実装

    以下のコードにより、状態遷移確率、方策の表を作成します。

    • コード
    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)
    
  3. 表のサンプル

    コードを実行すると以下のような状態遷移確率の表(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 にしています)にしています。このような方策の表の値を強化学習アルゴリズムにより報酬が最大になるように調整します。

4.2.2 状態遷移確率の実際の値を求める

表の作成 の 表5で、状態遷移確率の表を作成しましたが、値を全て 0 にしていました。ここでは、実際にシミュレータを動かして実際の値を求めます。

  1. 考え方

    (\(\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
  2. コード

    以下のコードにより、状態遷移確率(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を求め、それを用いて方策を更新することで学習を行います。

\begin{eqnarray} \bigl[\mathbf{P}^\pi\bigr]_{ss'} &=& \sum_a\pi(a|s)p(s'|s,a)\\ \bigl[\mathbf{R}^\pi\bigr]_{s} &=& \sum_a\pi(a|s)\sum_{s'}p(s'|s,a)r(s,a,s')\\ \mathbf{v} &=& \left(\mathbf{1} - \gamma\mathbf{P}^\pi\right)^{-1}\mathbf{R}^\pi\\ q_{\pi}(s,a) &=& \sum_{s'}p(s'|s,a)\biggl[r(s,a,s') + \gamma v_{\pi}(s')\biggr] \end{eqnarray}

(1)(2)(3)(4)は 参考文献 現場で使える!Python深層強化学習入門 強化学習と深層学習による探索と制御 P36, P32 より引用

4.3.2 コード

  1. 状態遷移確率(P) と方策の確率を合わせたものを求める

    以下のコードにより、状態遷移確率(P) と方策の確率を合わせたものを求めています。状態遷移確率と方策の行列の各要素の積をとり、行動 τ について和をとっています。

    # 状態遷移確率と方策の確率を合わせたものを求める
    P_PI_ss = Psas * PIas.values.reshape(-1, 1)
    P_PI_ss = P_PI_ss.sum(level=[0, 1])
    
  2. 状態遷移確率と方策から報酬の期待値を求める

    以下のコードにより、状態遷移確率と方策から報酬の期待値を求めます。状態遷移確率と方策と報酬の行列の各要素の積をとり、行動 τ について和をとっています。

    # 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])
    
  3. ベルマン方程式を解き状態価値を求める

    以下のコードにより、ベルマン方程式を解き状態価値を求めます。ベルマン方程式を解くのは以下の行列計算で完了します。計算は逆行列の計算を含みます。

    # 状態価値 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)
    
  4. 行動価値を求める

    以下のコードにより、行動価値を求めます。行動価値は即時報酬、状態価値を使用して計算できます。

    # 行動価値 Q を求める
    Q_sa = Psas*(Rsas + gamma*V_s.values.reshape(1,-1))
    Q_sa = Q_sa.sum(axis=1)
    
  5. 方策を更新する

    以下のコードにより、方策を更新します。ここでは、行動価値が最大になる行動を確率1でとるように方策を更新します。ここまでの 状態遷移確率(P) と方策の確率を合わせたものを求める方策を更新する までが学習に必要な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
    
  6. 方策が更新されなくなるか、指定回数を超えたら方策の更新を打ち切る

    方策が更新されなくなるか、指定回数を超えたら方策の更新を打ち切ります。方策を更新したのに、前の方策と同じになった場合は、その次も同じ、その次もと永遠と同じ方策になるため更新を打ち切ります。今回は、指定回数を超えても終了するようにプログラムしています。

    # 方策が更新されなくなったら終了
    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 出力結果

出力結果は以下になります。

https://raw.githubusercontent.com/keimina/gif/master/gymモジュール_Pendulum-v0_方策反復法_get_real_P_2.gif

図1. コード全体の実行結果

図1. が コード全体 実行結果ですが、頂上で止まっていません。ここから乱数のシードやコードを変えて試行錯誤を繰り返し手動でチューニングを行い学習させたところ、それらしい動きになりました。

https://github.com/keimina/gif/blob/master/gymモジュール_Pendulum-v0_方策反復法_get_real_P_good_2.gif?raw=true

図2. 手動チューニングを行い学習させた結果

7 考察

振り子の動きを見ると、左右に振って徐々に上昇させていることがわかります。これは、一方向に τ の上限のトルクを加えても、頂上まで上昇することがあまりないため、上昇させるには左右に振るのが定石だからです。(頂上付近だったり、初速度が大きい場合は一方向のトルクでも頂上にたどり着きます)また、左右に振るだけでも、振り子が一番下にいる確率を下げれますので、これだけで、マイナスの報酬を少なくできます。少なくともこの左右に振る動きを自動的に獲得することができました。ただし、頂上で静止するのは難しいようです。振り子が垂直の時、重力と釣り合うため、他の状態に比べ、少しのトルクで大きく動きます。そこに離散化による誤差の影響が加わり学習が難しくなっているかもしれません。ただ、上記図2.のように、パラメータによっては、うまく学習するときもあるようです。離散化の時に\(\dot \theta\)を他の分割数を他のパラメータの分割数より多くしているのは、他の分割数をあげるよりより、\(\dot \theta\)の分割数をあげる方が良い動きが見られたからです。離散化の方法も学習に大きな影響を及ぼすことがわかりました。出力結果をみると、真の最適な方策に収束しているようには思いませんが、それに近い方策になっていると思います。簡単な報酬を定義しただけで、アルゴリズムが自動的に最適に近い方策を出してくれるので便利です。

8 まとめ

Pendulum-v0 について、報酬を独自に定義し、連続値を離散化し状態遷移確率などを表で表し、ベルマン方程式を解くことにより、方策の学習を行いました。学習済のパラメータを用いて、シミュレーションを10回繰り返したところ、10回全てで振り子を頂上で止まらせることができました。(図2のアニメーションを参照)

9 参考文献

9.1 現場で使える!Python深層強化学習入門 強化学習と深層学習による探索と制御

9.2 これからの強化学習

www.morikita.co.jp

以上です。おやすみなさい。

gymモジュール Pendulum-v0 で振り子をトルクを与えて静止させる

前回に続き、gym モジュールの Pendulum-v0 の理解を深めます。今回は、手で振り子を水平状態にし、その状態をキープするようにトルク与えて、振り子を静止させるシミュレーションを行います。今回も、強化学習は行いません。シミュレータの理解を深めます。

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を与え、シミュレーションを1ステップ進める

表4. 前回のコードとの差分

5 出力結果

出力結果は以下になります。

f:id:keimina:20200308170625p:plain

図1. 出力結果1

水平状態で重力と釣り合うようなトルクを与えているため、動かないことがわかります。

ただ、動かないと、本当にプログラムが正しく動作しているかわからないため、トルクをわずかに減らして T =-4.9 として実行して見ます。振り子がフラフラし、プログラムが動作していることがわかります。

f:id:keimina:20200308160812g:plain

図2. 出力結果2

6 まとめ

gymモジュールの Pendulum-v0 の環境を使用し、 手で振り子を水平状態にし、その状態をキープするようにトルク与えて、振り子を静止させるシミュレーションを行いました。

以上です。

gymモジュール Pendulum-v0 で振り子を揺らす

前回に続き、gym モジュールの Pendulum-v0 の理解を深めます。今回は、 手で振り子を水平状態にし、離すことに相当するシミュレーションを行います。今回も、強化学習は行いません。シミュレータの理解を深めます。

1 前回の記事

前回の記事はこちらです。gymとは何か、pendulumとは何かについてはこちらをご覧ください。

keimina.hatenablog.jp

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 出力結果

こちらからトルクは与えませんが、重力はデフォルトでかかっている状態のため、シミュレーション結果をみると現実世界の振り子のような動きになることがわかります。減衰する様子は見られないため、摩擦はないことがわかります。

f:id:keimina:20200307194138g:plain

図1. 出力結果

6 まとめ

gymモジュールの Pendulum-v0 の環境を使用し、 振り子を水平状態にし、離すことに相当するシミュレーションを行いました。現実世界の振り子のような動きになることがわかりました。

以上です。

gymモジュール Pendulum-v0 で振り子を等速度で回転させる

gym モジュールの Pendulum-v0 とは以下のような振り子のシミュレータです。

f:id:keimina:20200306000655p:plain

図1. Pendulum-v0

今回は、このシュミレータを用いて、振り子を等速度で回転させるプログラムを書きます。振り子をアナログ時計の秒針のように一定の速度でグルグル回すことを目的とします。このプログラムを作ることで gymモジュール の API に触れ、gymの理解を深めます。振り子をグルグル回してAPI に触れるだけです。強化学習は使いません。もう一度言いますが、強化学習で振り子を等速度に制御する話ではありません。

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 コード

今回書いたコードは以下になります。 1秒間に 45° だけ回転します。 それを8回繰り返して1周させています。

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 出力結果

以下のように、振り子が1秒間に 45° だけ回転し、1週します。

f:id:keimina:20200306001238g:plain

図2. 出力結果

6 まとめ

gymモジュールの Pendulum-v0 の環境を使用し、 振り子を等速度で回転させました。(厳密にいうと等速度に回転しているように見えるようにプログラムした)env.env.state で振り子の位置(角度)を指定できることがわかりました。gymモジュールの API をほんの少し理解することができました。

以上です、おやすみなさい。

混合分布が確率であることについて

混合分布が確率であることについて書きます。混合分布とは、確率質量関数または確率密度関数を線形結合させたものだと理解しています。本を読んでいたところ、結合係数の和が1のとき、混合分布も確率質量関数または確率密度関数として解釈することができるということでした。個人的にこれがなぜかわからなかったので備忘録としてまとめました。

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)は確率質量関数としてみなせる。

以上です。おやすみなさい。

if文を自動生成する方法


if文を自動生成する方法について記載します。


1 したいこと



以下のような表があったとき、入力 (a, b) を f(a, b) に変換する関数を自動で作成したいと思います。


a b f(a,b)
0 0 1
0 1 0
1 0 1
1 1 0


表1.


この関数を手動でプログラムすると以下のようになるかと思います。


コード:


def f(a, b):
if a==0 and b==0:
ret = 1
if a==0 and b==1:
ret = 0
if a==1 and b==0:
ret = 0
if a==1 and b==1:
ret = 1
return ret


上記の関数の出力が表1.の出力と同じになるか確認します。


コード:


from itertools import product
for x in product([0, 1], [0, 1]):
line = map(str, x + (f(*x),))
print("|".join(line))
'''
コメント:
itertools.product はドット積のパターンを生成するのに使用しています。
for文を入れ子にしなくて良くなるので便利です。'''


出力:


0|0|1
0|1|0
1|0|0
1|1|1


したいことをまとめると、、、


上記の表1.のような、入力と出力の関係を定義した表が与えられたとき、上記「def f(…):」の中でプログラムしたようなif文を自動生成生成したいということです。


2 解決策 (1)



解決策 (1) では if文を自動生成するプログラムを自作します。以下のように実装しました。


コード:


import pandas as pd

# 入力となる表を DataFrame として作成する
df = pd.DataFrame(
[[0, 0, 1],
[0, 1, 0],
[1, 0, 0],
[1, 1, 1]],
columns=['a', 'b', 'ret']
)

# 評価を入力とし用いて、関数fを定義するための文字列を作成する
arg_t = ",".join(df.columns[:-1])
if_t = ""
for l in df.values:
s = " and ".join("%s==%d"%(col, i) for col, i in zip(df.columns[:-1], l[:-1]))
if_t += " if %s:\n"%s
if_t += " ret = %d\n"%l[-1]


f_t = "def f({}):\n"\
"{}"\
" return ret".format(arg_t, if_t)

# 上で定義した文字列をインタプリンタで評価し、関数fを上書きする
exec(f_t)

# 関数fが期待通りの動作になっているか確認する
from itertools import product
for x in product([0, 1], [0, 1]):
line = map(str, x + (f(*x),))
print("|".join(line))


出力:


0|0|1
0|1|0
1|0|0
1|1|1


if文をテキストとして書いてそれを exec関数で評価して、関数を定義しています。その関数を使用して表と同じ出力を出力しています。


※exec関数は、文字列をインタプリンタの入力として評価します。予期しないバグの温床になるので使用しないようにしましょう。また、任意のコードを実行可能なので実行するのはセキュリティ上の危険が伴います。ソースコードをよく理解した上で使用しましょう。

念のため生成した関数を表示してみます.


コード:


print('====================作成した関数の表示====================')
print(f_t)


出力:


====================作成した関数の表示====================
def f(a,b):
if a==0 and b==0:
ret = 1
if a==0 and b==1:
ret = 0
if a==1 and b==0:
ret = 0
if a==1 and b==1:
ret = 1
return ret


問題なく自動生成されていることがわかります。


3 解決策 (2)



解決策としては 解決策 (1) でも問題ないのですが、もう少しスマートな方法はないでしょうか。調べたところ sympy.SOPform 関数このような問題を解決するのにうってつけな関数でした。sympy.SOPform 関数を使用すると、表1.のような表を与えると、それに紐づく論理式を組み立ててくれます。さらに、冗長な論理式であった場合は冗長でない論理式にしてくれます。


コード:


import pandas as pd
from sympy import symbols
from sympy.logic import SOPform

# 入力となる表を DataFrame として作成する
df = pd.DataFrame(
[[0, 0, 1],
[0, 1, 0],
[1, 0, 0],
[1, 1, 1]],
columns=['a', 'b', 'ret']
)

ss = symbols(" ".join(df.columns[:-1]))
idx = df.iloc[:, -1]==1

if any(idx):
variables = ss
minterms = df[idx].values[:,:-1].tolist()
g = SOPform(variables, minterms)# 出力が 1 になるパターンだけを第二引数に設定する


# 評価を入力とし用いて、関数fを定義するための文字列を作成する
arg_t = ",".join(df.columns[:-1])
f_t = "def f({0}):\n"\
" args_bak = {0}\n"\
" {0} = symbols('{0}')\n"\
" expr = {1}\n"\
" return int(bool(expr.subs(dict(zip([{0}], args_bak)))))".format(arg_t, str(g))

exec(f_t)

# 関数fが期待通りの動作になっているか確認する
for line in df.values.tolist():
line = line[:-1] + [f(*line[:-1])]
print("|".join(str(i) for i in line))


出力:


0|0|1
0|1|0
1|0|0
1|1|1


作成した関数も見てみます。


コード:


print(f_t)


出力:


def f(a,b):
args_bak = a,b
a,b = symbols('a,b')
expr = (a & b) | (~a & ~b)
return int(bool(expr.subs(dict(zip([a,b], args_bak)))))


「(a & b) | (~a & ~b)」が sympy.SOPform 関数で作成した論理式です。自動で論理式を作成してくれて大変便利です。


次は、もっと大きい表で試してみます。まず、先ほどより大きめの表を作成します。


コード:


import numpy as np
from itertools import product

np.random.seed(1)
ret = np.random.choice([0, 1], 16)
columns = ['a', 'b', 'c', 'd']
df = pd.DataFrame(list(product([0,1], [0,1], [0,1], [0,1])), columns=columns)
df['ret'] = ret
print(df)


出力:


a b c d ret
0 0 0 0 0 1
1 0 0 0 1 1
2 0 0 1 0 0
3 0 0 1 1 0
4 0 1 0 0 1
5 0 1 0 1 1
6 0 1 1 0 1
7 0 1 1 1 1
8 1 0 0 0 1
9 1 0 0 1 0
10 1 0 1 0 0
11 1 0 1 1 1
12 1 1 0 0 0
13 1 1 0 1 1
14 1 1 1 0 1
15 1 1 1 1 0


この表に対して、 前の手順と同じことをします。


コード:


from sympy import symbols
from sympy.logic import SOPform

ss = symbols(" ".join(df.columns[:-1]))
idx = df.iloc[:, -1]==1

if any(idx):
variables = ss
minterms = df[idx].values[:,:-1].tolist()
g = SOPform(variables, minterms)

# 評価を入力とし用いて、関数fを定義するための文字列を作成する
arg_t = ",".join(df.columns[:-1])
f_t = "def f({0}):\n"\
" args_bak = {0}\n"\
" {0} = symbols('{0}')\n"\
" expr = {1}\n"\
" return int(bool(expr.subs(dict(zip([{0}], args_bak)))))".format(arg_t, str(g))

exec(f_t)

print('========関数fが期待通りの動作になっているか確認する========')
# 関数fが期待通りの動作になっているか確認する
lst = []
for line in df.values.tolist():
line = line[:-1] + [f(*line[:-1])]
print("|".join(str(i) for i in line))
lst.append(line)
print('===================df と lst を比較する====================')
print('df と lst の', '比較OK' if df.values.tolist() == lst else '比較NG')
print('====================作成した関数の表示====================')
print(f_t)


出力:


========関数fが期待通りの動作になっているか確認する========
0|0|0|0|1
0|0|0|1|1
0|0|1|0|0
0|0|1|1|0
0|1|0|0|1
0|1|0|1|1
0|1|1|0|1
0|1|1|1|1
1|0|0|0|1
1|0|0|1|0
1|0|1|0|0
1|0|1|1|1
1|1|0|0|0
1|1|0|1|1
1|1|1|0|1
1|1|1|1|0
===================df と lst を比較する====================
df と lst の 比較OK
====================作成した関数の表示====================
def f(a,b,c,d):
args_bak = a,b,c,d
a,b,c,d = symbols('a,b,c,d')
expr = (b & ~a) | (~a & ~c) | (b & c & ~d) | (b & d & ~c) | (a & c & d & ~b) | (~b & ~c & ~d)
return int(bool(expr.subs(dict(zip([a,b,c,d], args_bak)))))


作成した関数の出力と、表で定義した関数の出力が同じであることを確認し、問題なく、動作していることがわかります。


4 sympy.SOPform 関数について理解を深める



sympy.SOPform 関数を使用して自動で求まった論理式ですが、以下の sympy.logic.boolalg.to_dnf 関数を通しても同じなります。


コード:


from sympy.logic.boolalg import to_dnf
print(to_dnf(expr, simplify=True))


出力:


(b & ~a) | (~a & ~c) | (b & c & ~d) | (b & d & ~c) | (a & c & d & ~b) | (~b & ~c & ~d)


to_dnf関数とはなんでしょうか?関数のドキュメントをみてみます。


sympy.logic.boolalg.to_dnf(expr, simplify=False)

Convert a propositional logical sentence s to disjunctive normal
form. That is, of the form ((A & ~B & …) | (B & C & …) | …) If
simplify is True, the expr is evaluated to its simplest DNF form using
the Quine-McCluskey algorithm.


引数 simplify=True のとき Quine-McCluskey アルゴリズムで最も単純なDNF形式に変換されるとあります。参考文献No4. によると Quine–McCluskey アルゴリズムはブール関数を最小化するために使用されるアルゴリズムとのことです。また、参考文献No2. No3. から DNF とは論理式を表現する時に使える演算子を制限した論理式のフォーマットのようなものだと解釈しました。このような、分野をブール関数の簡単化(参考文献No6.)というようです。


まとめると、解決策 (2) では sympy.SOPform 関数を使用して、真理値表からブール関数を自動的に構築するだけでなく、そのブール関数の簡単化を行ったということになります。


5 参考文献



参考文献は以下の通りです。

No. タイトル URL 内容
1 Logic https://docs.sympy.org/latest/modules/logic.html sympy.SOPform関数について
2 連言標準形 https://ja.wikipedia.org/wiki/連言標準形 CNFについて
3 選言標準形 https://ja.wikipedia.org/wiki/選言標準形 DNFについて
4 Quine–McCluskey algorithm https://en.wikipedia.org/wiki/Quine–McCluskey_algorithm Quine–McCluskey algorithmについて
5 クワイン・マクラスキー法 https://ja.wikipedia.org/wiki/クワイン・マクラスキー法 クワイン・マクラスキー法について
6 ブール関数 https://ja.wikipedia.org/wiki/ブール関数 ブール関数について、簡単化について



以上です。

sympyで「... = ...」の式を「... = 0」の式にする方法

sympyで「… = …」の式を「… = 0」の式にする方法を記載します。

1 行いたいこと

まず、行いたいことを明確にします。例えば、以下の式があったとします。

a = b

上記の式の両辺から b を引く(b を左辺に移項する)と以下のようになります。

a - b = 0

今回行いたいことは、上記のように両辺から右辺の式を引く、または左辺の式を引くことにより「… = 0」の式を作成することです。

2 解決策

以下が解決策です。まず、 a = b を表す等式を作成する方法は以下になります。

from sympy import Symbol, Eq
a = Symbol('a')
b = Symbol('b')
left = a
right = b
e = Eq(left, right)
print(e)
Eq(a, b)

上記は、a と b を表す Symbol オブジェクトを left, right 変数に格納しています。a = b を表す等式を変数 e に格納するということを行っています。等式の作成には sympy.Eq 関数を使用しています。

両辺から右辺の式を引き算すれば「… = 0」の式になります。これを、sympy で書くと以下になります。

# 両辺から右辺の式 right を引く
e2 = Eq(left - right, right - right)
# 上記プログラムは以下と同じ
e2 = Eq(left - right, 0)
print(e2)
Eq(a - b, 0)

sympyで「… = …」の式を「… = 0」の式にする方法の説明は以上です。物足りないので練習問題を追加しておきます。

3 練習問題

練習として以下を解くことにします。

3.1 【問 1】

問:

以下の等式を移項して「... = 0」の等式にするにはどうすれば良いでしょうか?

  a = (b - 1.08614074 * a)**2

答え:

# Symbol の定義
a = Symbol('a')
b = Symbol('b')
# 右辺、左辺の定義
left = a
right = (b - 1.08614074 * a)**2
# (変形前)等式の定義
e = Eq(left, right)
print("変形前: ", e)
# (変形後)等式の定義
e2 = Eq(left - right, 0)
print("変形後: ", e2)
変形前:  Eq(a, (-1.08614074*a + b)**2)
変形後:  Eq(a - (-1.08614074*a + b)**2, 0)

3.2 【問 2】

問:

【問 1】の式を b について解く(b = f(a) の形に変形する)にはどうす
れば良いでしょうか?

※式を解くには sympy.solve 関数を使用します。第一引数に等式を指定します。
【問 1】のソースコードの e, e2 どちらを使用しても問題ありません)第二引
数にどのSymbolについて解くかを指定します。戻り値は要素が解のリストです。

答え:

from sympy import solve
# sympy.solve 関数を使用して b について解く
answer = solve(e, b)
# 上記プログラムは以下と同じ
answer = solve(e2, b)
print(answer)
[-sqrt(a) + 1.08614074*a, sqrt(a) + 1.08614074*a]

以上です。