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]

以上です。

ヘビサイド関数を使って矩形関数を作る方法(Sympy)

0 or 1 をとるシンプルな関数として、ヘビサイド関数があります。本記事は、 Python のライブラリである Sympy の Heaviside 関数を使って、矩形関数を実装したときの記事です。ヘビサイド関数を用いて矩形関数を作ることができました。

1 参考文献

以下のサイトが参考になりました。

oeis.org

2ビサイド関数とは

ビサイド関数とは以下のようなグラフを出力する関数です。

f:id:keimina:20200111151158p:plain

図1. ヘビサイド関数

3 矩形関数とは

矩形関数とは以下のようなグラフを出力する関数です。

f:id:keimina:20200111151257p:plain

図2. 矩形関数

4 矩形関数をヘビサイド関数で表すには

矩形関数をヘビサイド関数で表すには、簡単に言ってしまうと、上に凸で正の値をとる2次関数の値をヘビサイド関数に与えると、矩形関数になります。以下にこれを実装したコードを記載します。

5 コード

矩形関数をヘビサイド関数で実装したコードは以下になります。

# ライブラリ読み込み(数値計算、グラフ出力用)
import numpy as np
import matplotlib.pyplot as plt

# ライブラリ(数式モデル構築用)
from sympy import Heaviside

# ヘビサイド関数は sympy.Heaviside を使用する
# numpy のベクトルに関数を適用できるように np.vectorize を使用する
f = lambda x: Heaviside(-(x - 2)*(x + 2), 0)
f = np.vectorize(f)

# 代入する x の値を定義する
xdata = np.linspace(-6, 6, 100)

# 代入する y の値を求める
ydata = f(xdata)

# 求めた xdata, ydata をプロットする
plt.plot(xdata, ydata)
plt.show(False)

6 解説

ビサイド関数は sympy.Heaviside を使用します。 sympy.Heaviside ではx=0 の時の戻り値を第2引数で設定することができます。上記コードでは第2引数に 0 を指定していることに注意します。

7 出力

以下のグラフが出力されます。-(x - 2) * (x + 2) が正になる範囲、すなわち -2 < x < 2 の範囲で 1 に、それ以外の範囲では 0 になることがわかります。

f:id:keimina:20200111151257p:plain

図3. 矩形関数

8 まとめ

sympy ライブラリの関数 Heaviside を用いて、 矩形関数を実装しました。上に凸で正の値をとる2次関数の値をヘビサイド関数に与えると、矩形関数になることがわかりました。

以上です。