やってみた!

やってみた!

試したことを中心に、書評や興味のあること、思ったこととか

PyBullet-HumanoidFlagrunHarderBulletEnv-v0(2)

2020/1/25改正
 学習継続時に早くalphaが収束するようalpha、log_alpha、alpha_optimizerを保存するように変更しました。gpu有、無しの両環境で保存データを共有できるようモデル読み込み時にmap_location=deviceを追加しました。

2020/1/23改正
 BATCH_SIZEを128から論文と同じ256に変更したら学習が飛躍的に改善しました。あわせてSIGMA_DECAY を 3000から10に、num_episodesを9000に変更しました。

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

 DDPGではHumonoidFlagrunをうまく学習させられなかったので、SAC(Soft Actor-Critic)で学習させてみました。以下まとめ。

  •  SAC動かすのにかなり手間取ったので、まだHumanoidFlagrunHarderを学習させることができていません(できないかも・・・汗)。本記事は、やっと学習するようになったHumanoidBulletEnv-v0を例に紹介します。
  • 学習に8時間かけて、やっと、すぐには転ばず、足を出す動作を覚えました。さらに時間をかければもっと安定するようになるでしょう。学習しているのを確認するだけなら2,3時間だけでもいいと思いますが、なにせ時間がかかります。

1.SAC(Soft Actor-Critic) 

 SAC(Soft Actor-Critic)は2018年に発表された比較的新しい方式で、DDPGと同様連続値を扱うことができる強化学習法のひとつです。以下の論文は初期のSACの計算手順の改良版について書かれたもので、今回はこれをベースにしました。本論文にはSACはDDPGや他の手法と比較して学習の進みが早いといった特長があることが示されています。

 本論文はGoogleカルフォルニアバークレー校の共同作業のようです。

 こちらも参考にさせてもらいました。なおSpinning Upの内容は改良前のSAC計算方法なので若干異なりますが、SACの式の解説等、基本的考え方は参考になります。

Soft Actor-Critic — Spinning Up documentation

GitHub - openai/spinningup: An educational resource to help anyone learn deep reinforcement learning.

GitHub - ku2482/soft-actor-critic.pytorch: A PyTorch Implementation of Soft Actor-Critic.

 次の章でとにかくコードを動かしてみますが、その前にSACのイメージだけでも。

  • 方策policyの出力(行動)を確率的なものとします。policyの出力は中心値に対してgauss分布でばらつくものとし、中心値とばらつきの大きさエントロピーニューラルネットワークで推定します。
  • 価値関数は、報酬に加え方策のエントロピー(ばらつきの大きさ)の大きさを加えます。これにより、より大きなエントロピーを持つ行動を学習するようになります。不安定な極所的最適解が学習から排除され、学習が安定して進むものと思われます。

     例えば水上の飛び石の上を渡っているとすると、DDPGでは単純に乱数で選んで学習するので小さな石を選び、うまくいったらそれを学習してしまうのに対し、SACはある幅をもった行動の価値を高いと考えるので、常に大きな石を選ぶように学習します。

    f:id:akifukka:20200121192727j:plain
     この積み重ねでSACはより安定した行動を学習し、極所最適解にはまったりせず学習が進むということだと思われます。

2.とにかく動かしてみる

 SACの細かい話は次回以降にして、とにかく動かしてみましょう。なお、ある程度動きが改善されてるなと感じるくらいまで学習させるには8~10時間くらいの覚悟が必要です。

(1) Colaboratoryを立ち上げる

 ランタイム-ランタイムのタイプを変更でGPUにします。

(2)ライブラリをインストールします。

 MountainCarContinuousも動かせるようBox2Dもインストールしています。

!apt-get update
!apt-get -qq -y install xvfb freeglut3-dev ffmpeg
!pip3 -q install pyglet
!pip3 -q install pyopengl
!pip3 -q install pyvirtualdisplay
!apt-get install x11-utils
!pip3 install pybullet

!pip3 install 'gym[Box2D]'
!pip3 install  'gym[classic_control]'

(3)Colaboratoryのランタイムを再起動します。

 インストールしたライブラリを有効にするためランタイムを再起動します。

 ランタイム-ランタイムの再起動

(4)Googleドライブをマウントします。

 学習したモデルを保存するためGoogleドライブをマウントします。
 実行するとリンクが表示されるのでリンクをクリックして、アカウントを選択して表示されたコードをCtrl-Cでコピーして、Enter your authorization code:の四角のエリアにCtrl-Pでペーストして、Entキーを押します。

from google.colab import drive
drive.mount('/content/drive')

(5) ネットワーク、オプティマイザ、リプレイメモリ等を生成し初期化します。

#初期化コード
import gym
import pybullet_envs
import math
import random
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from collections import namedtuple
#from itertools import count

import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
import torchvision.transforms as T

from pyvirtualdisplay import Display
from scipy.interpolate import interp1d

display = Display(visible=0, size=(1024, 768))
display.start()
import os
#pyvirtualdisplayの仕様変更で次の文はエラーになるのでコメントアウトします(2021/6/5修正)
#os.environ["DISPLAY"] = ":" + str(display.display) + "." + str(display.screen)

env = gym.make('HumanoidBulletEnv-v0')
#env = gym.make('HumanoidFlagrunBulletEnv-v0')
#env = gym.make('HumanoidFlagrunHarderBulletEnv-v0')
#env = gym.make('MountainCarContinuous-v0')

# if gpu is to be used
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

#リプレイメモリー
Transition = namedtuple('Transition',
                        ('state', 'action', 'next_state', 'reward','done'))

class ReplayMemory(object):
    def __init__(self, capacity):
        self.capacity = capacity
        self.memory = []
        self.position = 0

    def push(self, *args):
        """Saves a transition."""
        if len(self.memory) < self.capacity:
            self.memory.append(None)
        self.memory[self.position] = Transition(*args)
        self.position = (self.position + 1) % self.capacity

    def sample(self, batch_size):
        return random.sample(self.memory, batch_size)

    def __len__(self):
        return len(self.memory)
   
    def reset(self):
        self.memory = []
        self.position = 0   

#Actor(Policy) Network squashed Gaussian Policy
LOG_STD_MAX = 2
LOG_STD_MIN = -20
EPS = 1e-8

class Policy_net(nn.Module):
    def __init__(self, NumObs, NumAct,NumHidden):
        super(Policy_net, self).__init__()
        self.ln1 = nn.Linear(NumObs, NumHidden)
        self.ln2 = nn.Linear(NumHidden, NumHidden)
        self.ln3 = nn.Linear(NumHidden, NumHidden)

        #mu 中心値算出用
        self.muln = nn.Linear(NumHidden, NumAct)

        #log_std エントロピー×ー1算出用
        self.log_stdln = nn.Linear(NumHidden, NumAct)

    def forward(self, x):
        x = F.relu(self.ln1(x))
        x = F.relu(self.ln2(x))
#層数が多いと収束が遅いことがあるので層数減らした。収束が早くなったかは未確認。
#        x = F.relu(self.ln3(x))

        mu = self.muln(x)
        log_std = torch.tanh(self.log_stdln(x))
        log_std = LOG_STD_MIN + (LOG_STD_MAX - LOG_STD_MIN) * (log_std + 1) * 0.5

        log_std_exp = log_std.exp()
        #πの計算(action) ノイズはガウス分布
        noise = torch.randn_like(mu)
        pi = mu + noise * log_std_exp
        #Gaussian Policy類推度の計算(Log-Likelihood)
        pre_sum = -0.5 * (((pi-mu)/(log_std_exp + EPS))**2 + 2 * log_std + np.log(2 * np.pi))
        logp_pi = torch.sum(pre_sum,1,keepdim=True)

        #squash tanhで範囲をー1~1につぶす
        mu, pi, logp_pi = apply_squashing_func(mu, pi, logp_pi)

        return mu,pi,logp_pi

#傾きを保持したままクリッピングだけを行う
def clip_but_pass_gradient(x, low=-1., high=1.):
    clip_high = (x > high).float()
    clip_low = (x < low).float()
    return x + ((high - x)*clip_high + (low - x)*clip_low).detach()

def apply_squashing_func(mu, pi, logp_pi):
    mu = mu.tanh()
    pi = pi.tanh()

    logp_pi = logp_pi - torch.sum((clip_but_pass_gradient(1 - pi**2,low = 0, high = 1) + 1e-6).log(),1,keepdim=True)
    return mu, pi, logp_pi

#Critic(Q) Network Q値の計算
class Q_net(nn.Module):
    def __init__(self, NumInput, NumHidden):
        super(Q_net, self).__init__()
        self.ln1 = nn.Linear(NumInput, NumHidden)
        self.ln2 = nn.Linear(NumHidden, NumHidden)
        self.ln3 = nn.Linear(NumHidden, NumHidden)
        self.ln4 = nn.Linear(NumHidden, 1)
        self.ln4.weight.data.uniform_(-3e-3, 3e-3)

    def forward(self, input):
        x = F.leaky_relu(self.ln1(input))
        x = F.leaky_relu(self.ln2(x))
#層数が多いと収束が遅いことがあるので層数減らした。収束が早くなったかは未確認。
#        x = F.leaky_relu(self.ln3(x))
        y = self.ln4(x)
        return y

def optimize_model():
    global alpha

    if len(memory) < BATCH_SIZE:
        return
    transitions = memory.sample(BATCH_SIZE)
    batch = Transition(*zip(*transitions))

    #batchを一括して処理するためtensorのリストを2次元テンソルに変換
    state_batch = torch.stack(batch.state,dim=0)
    action_batch = torch.stack(batch.action,dim=0)
    reward_batch = torch.stack(batch.reward,dim=0)
    next_state_batch = torch.stack(batch.next_state,dim=0)
    done_batch = torch.stack(batch.done,dim=0)

    #Qターゲット
    next_mus,next_actions,next_logp_pis = policy_net(next_state_batch)
    next_q1s = q1_target_net(torch.cat([next_state_batch, next_actions], 1))
    next_q2s = q2_target_net(torch.cat([next_state_batch, next_actions], 1))
    next_qs = torch.min(next_q1s, next_q2s)

    q_targets = reward_batch + GAMMA * (1.0-done_batch) * (next_qs - alpha * next_logp_pis)

    #Q1ネット
    #loss関数
#    q1_loss = 0.5 * F.mse_loss(q1_net(torch.cat([state_batch, action_batch],1)),q_targets) 中身がわかるような記述に変更
    q1_loss = 0.5 * torch.mean((q1_net(torch.cat([state_batch, action_batch],1)) - q_targets)**2)
    #ネットの学習
    q1_optimizer.zero_grad()
    #誤差逆伝搬
    q1_loss.backward(retain_graph=True)
    #重み更新
    q1_optimizer.step()

    #Q2ネット
    #loss関数
#    q2_loss = 0.5 * F.mse_loss(q2_net(torch.cat([state_batch, action_batch],1)),q_targets) 中身がわかるような記述に変更
    q2_loss = 0.5 * torch.mean((q2_net(torch.cat([state_batch, action_batch],1)) - q_targets)**2)
    #ネットの学習
    q2_optimizer.zero_grad()
    #誤差逆伝搬
    q2_loss.backward(retain_graph=True)
    #重み更新
    q2_optimizer.step()

    #policyネットのloss関数
    mus,actions,logp_pis = policy_net(state_batch)
    q1s = q1_net(torch.cat([state_batch,actions],1))
    q2s = q2_net(torch.cat([state_batch,actions],1))
    qs = torch.min(q1s, q2s)
    p_loss = torch.mean(alpha * logp_pis - qs)

    #policyネットの学習
    p_optimizer.zero_grad()
    #誤差逆伝搬
    p_loss.backward(retain_graph=True)
    #重み更新
    p_optimizer.step()

    #エラー"one of the variables needed for gradient computation has been modified by an inplace operation"が出るようになった
    #policy更新された後にlogp_pisを更新していないため。再計算追加(2021/6/6)
    mus,actions,logp_pis = policy_net(state_batch)
    alpha_loss = -torch.mean(log_alpha*(target_entropy + logp_pis))

    #log_alpha更新
    alpha_optimizer.zero_grad()
    #誤差逆伝搬
    alpha_loss.backward(retain_graph=True)

    #log_alpha更新
    alpha_optimizer.step()
    #alpha更新
    alpha = log_alpha.exp()

    tau = 0.005
    #q1,2 targetネットのソフトアップデート
    #学習の収束を安定させるためゆっくり学習するtarget netを作りloss関数の計算に使う。
    #学習後のネット重み×tau分を反映させ、ゆっくり追従させる。
    for target_param, local_param in zip(q1_target_net.parameters(), q1_net.parameters()):
      target_param.data.copy_(tau*local_param.data + (1.0-tau)*target_param.data)

    for target_param, local_param in zip(q2_target_net.parameters(), q2_net.parameters()):
      target_param.data.copy_(tau*local_param.data + (1.0-tau)*target_param.data)

#action,observationの要素数取得
n_actions = env.action_space.shape[0]
n_observations = env.observation_space.shape[0]
target_entropy = -torch.Tensor([n_actions]).to(device)
print('num actions,observations,target_entropy',n_actions,n_observations,target_entropy)

#ネットワーク
#policyネットとそのtargetネット
policy_net = Policy_net(n_observations, n_actions,256).to(device)
#Qネット,Vネットとそのtargetネット
q1_net = Q_net(n_observations + n_actions, 256).to(device)
q2_net = Q_net(n_observations + n_actions, 256).to(device)
q1_target_net = Q_net(n_observations + n_actions, 256).to(device)
q2_target_net = Q_net(n_observations + n_actions, 256).to(device)
#Qターゲットネット初期化(コピー)
q1_target_net.load_state_dict(q1_net.state_dict())
q2_target_net.load_state_dict(q1_net.state_dict())

#学習用optimizer生成
p_optimizer = optim.Adam(policy_net.parameters(),lr=3e-4)
q1_optimizer = optim.Adam(q1_net.parameters(),lr=3e-4)
q2_optimizer = optim.Adam(q2_net.parameters(),lr=3e-4)
log_alpha = torch.zeros(1,requires_grad=True,device=device)
alpha = log_alpha.exp()
alpha_optimizer = optim.Adam([log_alpha], lr=3e-4)

#学習用リプレイメモリ生成
memory = ReplayMemory(1000000)
memory.reset()

#総数カウント用
total_episode = 0
total_step = 0

 (6)学習

 学習します。1000エピソードごとにネットワークをGoogleドライブに保存します。ファイル名はmodel_Humanoid_SACxxxxで、xxxxはエピソード数です。

 num_epsode=20000を書き換えて学習episode数を変更します。最初は1000~5000くらいにして、様子を見た方がいいかも知れません。この学習コードだけ繰り返して実行すると、前の状態から学習を継続してくれます。

 また、後で説明するモデル読み込んだ後、学習コードを実行するとそこから継続して学習できます。

 学習を途中で止めても、再実行で継続できます(ごくたまに失敗しますが)。

 (5)の初期化コードを実行すると直前までの学習モデルや、ファイルから読み込んだモデルが初期化されてしまうので注意が必要です。

 なおColaboratoryは90分セッションが切れているとリセットされてしまうということで、次の記事の通りChromeのアドオンAuto Refreshを使います。

Google Colaboratoryの90分セッション切れ対策【自動接続】 - Qiita

#初期化コードを再実行せず、このコードだけを繰り返し実行すると直前の学習状態、
#もしくはファイルから読み込んだモデルを続けて学習できる。

#BATCH_SIZE = 128
BATCH_SIZE = 256
#Qネット学習時の報酬の割引率
GAMMA = 0.99
#alphaは自動計算のためコメントアウト
#alpha = 0.05

#episode数に対するノイズの減少係数
#OUNoise オリジナルのSACでは使わないが試行錯誤でいれてみた。
#ちゃんと確認できていないがMountainCarでは効果あるかも。
SIGMA_START = 1.0 #最初
SIGMA_END = 0.0 #最後
SIGMA_DECAY = 10 #このepisode回数で約30%まで減衰
SIGMA_DECAY_END = SIGMA_DECAY * 1.2 #このepisode以降はOUNoiseを使用しない
theta = 0.08

#50 フレーム/sec グラフ描画用
FPS = 50
#学習エピソード数。num_epsode=20000を変更する。
num_episodes = 9000
display_on = 1 # num_episodeの最後の回の動画データ保存
save_model_on = 1 # 0 モデル自動セーブしない 1:モデル自動セーブする
t_max = 1400 # 1エピソードでの最大ステップ数。ステップ数がこの値以上になったら強制的にenvをリセットして次のエピソードに移行

noise =np.array([random.uniform(-0.5,0.5) for i in range(n_actions)],dtype = np.float)

for i_episode in range(num_episodes):
    #whileループ内初期化
    observation = env.reset()
    done = False
    reward_total = 0.0
    logp_pi_total = 0.0
    t = 0
    #グラフ用データ保存領域初期化
    frames = []
    observations = []
    actions = []
    ts = []
    for i in range(n_observations):
      observations.append([])
    for i in range(n_actions):
      actions.append([])

    sigma = SIGMA_END + (SIGMA_START - SIGMA_END) * math.exp(-1. * total_episode / SIGMA_DECAY)

    total_episode += 1

    while not done and t < t_max:

        with torch.no_grad():
          mu,action,logp_pi = policy_net(torch.from_numpy(observation).float().reshape(1,-1).to(device))

          mu = mu.cpu().data.numpy().reshape(-1)
          action = action.cpu().data.numpy().reshape(-1)
          logp_pi = logp_pi.cpu().data.numpy().reshape(-1)[0]
          logp_pi_total += logp_pi

          #OUNoise 最初だけ。total_episode > SIGMA_DECAY_ENDではOUNoiseを加算しない。
          noise = noise - theta * noise + sigma * np.array([random.uniform(-1.0,1.0) for i in range(len(noise))])
          if total_episode < SIGMA_DECAY_END :
            action = mu + noise

          #episodeの最後は純粋にネットワークの性能を取得するためノイズ無し(中心値mu)-------------------
          if (i_episode == num_episodes -1 ):
            action = mu
          action = np.clip(action, -1, 1)

        #物理モデル1ステップ---------------------------
        next_observation, reward, done, info = env.step(action)
        reward_total = reward_total + reward

        #学習用にメモリに保存--------------------------
        tensor_observation = torch.tensor(observation,device=device,dtype=torch.float)
        tensor_action = torch.tensor(action,device=device,dtype=torch.float)
        tensor_next_observation = torch.tensor(next_observation,device=device,dtype=torch.float)
        tensor_reward = torch.tensor([reward],device=device,dtype=torch.float)
        tensor_done =  torch.tensor([done],device=device,dtype=torch.float)
        memory.push(tensor_observation, tensor_action, tensor_next_observation, tensor_reward,tensor_done)

        #observation(state)更新------------------------
        observation = next_observation

        #学習 現バージョンは毎ステップごとに学習。
        #episode終了毎にまとめて学習させることもできるが計算が遅くなったような気がしたので。
#        if ((t == t_max -1) or done) and (i_episode != num_episodes -1):
#          for j in range(t):
        optimize_model()

        #データ保存------------------------------------
        if (i_episode == num_episodes -1) and (display_on == 1) :
#          print(i_episode,observation,mu,action,_logp_pi)
          #動画
          frames.append(env.render(mode = 'rgb_array'))
          #グラフ 現状未使用
          for i in range(n_observations):
            observations[i].append(observation[i])
          for i in range(n_actions):
            actions[i].append(action[i])
          ts.append(t/FPS)

        #時間を進める----------------------------------
        t += 1
        total_step += 1
        # end while loop ------------------------------

    #10episode毎に、最後のepisodeの報酬、エントロピーを表示
    if (total_episode % 10 == 0) or (i_episode == num_episodes -1 ):
      print('episode,step,reward,entropy=',total_episode,total_step,reward_total,-logp_pi_total/t)
    
    #1000episodeごとにモデルをGoogleドライブに保存
    if (total_episode % 1000 == 0) and (save_model_on == 1):
      torch.save({
        'epoch':total_episode,
        'step':total_step,
        'log_alpha':log_alpha,
        'alpha':alpha,
        'policy':policy_net.state_dict(),
        'q1':q1_net.state_dict(),
        'q2':q2_net.state_dict(),
        'q1_target':q1_target_net.state_dict(),
        'q2_target':q2_target_net.state_dict(),
        'p_optimizer':p_optimizer.state_dict(),
        'q1_optimizer':q1_optimizer.state_dict(),
        'q2_optimizer':q2_optimizer.state_dict(),
        'alpha_optimizer':alpha_optimizer.state_dict(),
      },'/content/drive/My Drive/model_Humanoid_SAC' + str(total_episode) )
# end for loop --------------------------------------

実行すると10episodeごとに次のように表示されます。

episode,step,reward,entropy= 10 188 -63.53772342780692 0.4764575527773963
episode,step,reward,entropy= 20 374 -53.021724941738086 10.449289242426554
episode,step,reward,entropy= 30 581 -67.98715715128392 8.465339357202703
 

 ここで、各数値の意味は次の通りです。

episode:総エピソード数、step:総ステップ数、entropy:エントロピー

 最初はすぐに転倒してしまうため、1episodeのstepが少ないので、stepはあまり増えません(上の例では、10episodeで200stepしか進まないので平均すると1episode=20step程度で転倒しています)。

 なお、humanoidでは、エントロピーが-17になるようalpha(エントロピー調整係数(Qを求める際に、報酬にエントロピーを加算する時にエントロピーに乗ずる係数)、 inverse of reward scaleとも呼ばれる)を自動調整しているので、学習が進むとエントロピーもー17近辺に近づいていきます。

 (7)画像を表示

 画像を生成します。またgif動画ファイルを作ってGoogleドライブに保存します。

import matplotlib.pyplot as plt
import matplotlib.animation
import numpy as np
from IPython.display import HTML
from matplotlib.animation import PillowWriter

fig = plt.figure(figsize=(5, 4))
plt.axis('off')
plt.subplots_adjust(left=0, right=1, bottom=0, top=1)

images = []
for f in frames:
  image = plt.imshow(f)
  images.append([image])
ani = matplotlib.animation.ArtistAnimation(fig, images, interval=30, repeat_delay=1)
ani.save("/content/drive/My Drive/images.gif", writer="pillow")

HTML(ani.to_jshtml())

  (8)モデルを読み込み

 Googleドライブに自動保存されたモデルを読み込みます。そこから学習を継続することも可能です。ただしリプレイメモリの内容は保存されていないので、過去の試行例は引き継がれないので、純粋な意味では条件は変わります。

#Googleドライブからモデルを読み込む。ファイル名は都度変更すること。
checkpoint = torch.load('/content/drive/My Drive/model_Humanoid_SAC16000', map_location=device)
total_episode = checkpoint['epoch']
total_step = checkpoint['step']
log_alpha = checkpoint['log_alpha']
alpha = checkpoint['alpha']
policy_net.load_state_dict(checkpoint['policy'])
q1_net.load_state_dict(checkpoint['q1'])
q2_net.load_state_dict(checkpoint['q2'])
q1_target_net.load_state_dict(checkpoint['q1_target'])
q2_target_net.load_state_dict(checkpoint['q2_target'])
p_optimizer.load_state_dict(checkpoint['p_optimizer'])
q1_optimizer.load_state_dict(checkpoint['q1_optimizer'])
q2_optimizer.load_state_dict(checkpoint['q2_optimizer'])
alpha_optimizer.load_state_dict(checkpoint['alpha_optimizer'])

3.結果です。

  約連続98時間学習した結果です。論文では100万stepの学習で報酬5000だとか・・・。当方、137万stepで報酬291・・・。
 BATCH_SIZEを128から256に変更したところ大きく改善しました。9000 episodesで約10時間かかります。ちなみに最大ステップ数1400を変更するとrewardも増えます。

BATCH_SIZEが128:episode:16352、step:1378063 reward:291
BATCH_SIZEが256:episode:9002、step:1626911 reward:563

f:id:akifukka:20200125075100g:plain
つづく
 SACについて少し解説を書いた後、HumanoidFlagrunに挑戦する予定です。まだ全く手が付いていないので時間かかりそうですが・・。

 なにせ今回、最初はSACの初期の論文ベースで実装して、いきなりFlagrunをやってみたのですが全然ダメで。バグなのか課題が難しすぎるのか判断つかず2週間ほど迷走して、やっと現状までたどりついたところなので・・・。では!

 強化学習 カテゴリーの記事一覧 - やってみた! 

PyBullet-HumanoidFlagrunHarderBulletEnv-v0(1)

 次は3Dの物理シミュレータを使ってみます。以前はOpen AI Gymで使える3D物理環境は有料のMuJoCo用だけでしたが、今では無料で使えるPyBullet用環境(env)もあるということなので、こちらを使ってみます。

 PyBulletはErwin Cumansさんらが開発したオープンソースの3D物理シミュレーションツールです。ロボット業界で広く使われている物理モデル記述法であるURDFに対応しており応用範囲は広そうです。Erwin Cumansさんは今はGoogleで働いているそうな。

 もともとMuJoCo用に用意されたOpen AI Gymのenvが、PyBullet用に書き直されGithub上のPyBulletの以下のフォルダに収納されています。

bullet3/examples/pybullet/gym/pybullet_envs/

 手始め?に、あらかじめ用意されているenv、HumanoidFlagrunHarderBulletEnv-v0を試してみます。実は、DDPGで1週間ほど試行錯誤しながら学習させてみたのですが、いつまでたっても歩けるようにならず・・・。ということで今回は学習済みサンプルを使っての紹介です。

1.HumanoidFlagrunHarderBulletEnv-v0

 ソースコードはbullet3/examples/pybullet/gym/pybullet_envs/gym_locomotion_envs.py

です。人型のロボットが赤い球状のFlag(旗)を走って追いかける動作を学習するための環境です。

(1)動かしてみる

 PyBulletの以下フォルダに学習済みactorを使ったデモプログラムが用意されています。このままだとColaboratoryで絵を描けないので改造して動かします。

bullet3/examples/pybullet/gym/pybullet_envs/examples/enjoy_TF_HumanoidFlagrunHarderBulletEnv_v1_2017jul.py

①Colaboratoryノートブックの作成

 ColaboratoryでPython3の新しいノートブックを作成、ランタイム-ランタイムのタイプを変更でハードウェアアクセラレータにGPUをセットします。

 ノートブックに適当な名前をつけます。

②ライブラリをインストールします

 今までと違い、最後にpybulletをインストールする行が追加されています。

!apt-get update
!apt-get -qq -y install xvfb freeglut3-dev ffmpeg
!pip3 -q install pyglet
!pip3 -q install pyopengl
!pip3 -q install pyvirtualdisplay
!apt-get install x11-utils
!pip3 install pybullet

③Runtimeを再起動

 インストールしたライブラリを有効にするためランタイムを再起動します。

 ランタイム-ランタイムの再起動

④学習済みactorのクラス定義ファイルをGithubからコピー

!wget https://raw.githubusercontent.com/bulletphysics/bullet3/master/examples/pybullet/gym/pybullet_envs/examples/enjoy_TF_HumanoidFlagrunHarderBulletEnv_v1_2017jul.py

⑤環境を動かします

from enjoy_TF_HumanoidFlagrunHarderBulletEnv_v1_2017jul import SmallReactivePolicy
import gym
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.animation
from PIL import Image
from IPython.display import HTML
from pyvirtualdisplay import Display

display = Display(visible=0, size=(1024, 768))
display.start()

import os
#pyvirtualdisplayの仕様変更で次の文はエラーになるのでコメントアウトします(2021/6/5修正)
#os.environ["DISPLAY"] = ":" + str(display.display) + "." + str(display.screen)

env = gym.make('HumanoidFlagrunHarderBulletEnv-v0')

pi = SmallReactivePolicy(env.observation_space, env.action_space)

frames = []
t = 0
done = False
reward_total = 0.0

observation = env.reset()

while not done and t < 1400:

  action = pi.act(observation)

  #物理モデル1ステップ---------------------------
  observation, reward, done, info = env.step(action)
  reward_total = reward_total + reward

  #動画
  frames.append(env.render(mode = 'rgb_array'))

  #時間を進める----------------------------------
  t += 1
  # end while loop ------------------------------

print('reward=',reward_total)

fig = plt.figure(figsize=(6, 6))
plt.axis('off')

images = []
for f in frames:
  image = plt.imshow(f)
  images.append([image])
ani = matplotlib.animation.ArtistAnimation(fig, images, interval=30, repeat_delay=1)

HTML(ani.to_jshtml())

f:id:akifukka:20200125081754g:plain
(2)報酬について

 ソースコード、bullet3/examples/pybullet/gym/pybullet_envs/gym_locomotion_envs.pyを読むと報酬は次の総和で与えられるようです。詳しくはWalkerBaseBulletEnvのdef step内のself.rewardsから内容を追いかけることができます。

  • self._alive
    倒れている:-1
    立っている:1~2の間でbodyのz座標による。通常歩行している高さなら2。
    robot_locomotors.pyのclass HumanoidFlagrunHarder、def alive_bonus参照
  • progress(進捗)
    potentialの差(進捗分)。potentialは次の値で与えられる。ロボットが立ち上がるのは姿勢の報酬があるため。
    姿勢      :立っている200、寝ている100
    目標との距離  :距離/step時間 potentialの進捗が目標接近速度になる
    class HumanoidFlagrunHarder、def calc_potential参照。
  • electricity_cost
    関節の速度とトルクから計算したDCモータに対するコスト(負の値になる)
    gym_locomotion_envs.pyのWalkerBaseBulletEnvのdef step参照。
  • joints_at_limit_cost
    関節が限界まで曲がっている、伸びきっている時の減点。
    gym_locomotion_envs.pyのWalkerBaseBulletEnvのdef step参照。
  • feet_collision_cost
    たぶん未使用で0のまま

 

今回のまとめ

 3Dの物理シミュレータPyBulletを使ったOpen AI Gym用env、HumanoidFlagrunHarderBulletEnv-v0を学習済みサンプルを使い、Google Colaboratoryで動かしてみました。

 また、DDPGで学習に挑戦してみましたが、いつまでたっても歩けるようになりませんでした(うまくいかなかったので記事にはしていません)。

 ということで、単純なDDPGではなくて、その改良型

 Twin Delayed DDPG — Spinning Up documentation

 Soft Actor-Critic — Spinning Up documentation

といったあたりに挑戦する必要がありそうです・・・。時間かかりそうですね。DDPGでもふらふら歩けるようになると思っていましたが甘かった。では!

つづく

強化学習 カテゴリーの記事一覧 - やってみた! 

Open AI Gym Box2D BipedalWalkerをColaboratoryで動かしてみる(7)

 今回はおまけということで、DDPGに教師を追加してみました。記事の最初の方で作成したPD制御を教師として、DDPGの経験処理中に行動をアシストすると、学習に何か効果があるかを試しました。

 結果、最初は教師の影響を受けて大股で歩こうとしていたものの、最終的にはちょこちょこ走りになってしまいました。学習回数も教師無しより余計にかかってしまいました。これは教師自体がへたくそ(すぐつまづく)だったこともあると思いますが。。

 その他、Colaboratoryが動画保存中に異常終了することが多かったので、学習だけ実行してネットワークをGoogleドライブに保存し、後でネットワークを読み込んで動画を作るように等、少し改良しました。

 BipedalWalkerは本記事で最後なので、これまでの記事の繰り返しになりますが、やり方を一通り最初からまとめておきます。

1. Colaboratoryを立ち上げる

 ランタイム-ランタイムのタイプを変更でGPUにします。

2.ライブラリをインストールします。

!apt-get update
!pip3 install box2d-py
!pip3 install 'gym[Box2D]'
!pip3 install  'gym[classic_control]'
!apt-get -qq -y install xvfb freeglut3-dev ffmpeg
!pip3 -q install pyglet
!pip3 -q install pyopengl
!pip3 -q install pyvirtualdisplay
!apt-get install x11-utils

 3.Colaboratoryのランタイムを再起動します。

 インストールしたライブラリを有効にするためランタイムを再起動します。

 ランタイム-ランタイムの再起動

4.Googleドライブをマウントします。

 学習したモデルを保存するためGoogleドライブをマウントします。
 実行するとリンクが表示されるのでリンクをクリックして、アカウントを選択して表示されたコードをCtrl-Cでコピーして、Enter your authorization code:の四角のエリアにCtrl-Pでペーストして、Entキーを押します。

from google.colab import drive
drive.mount('/content/drive')

 5.学習前の処理

 ネットワーク、オプティマイザ、リプレイメモリ等を生成し初期化します。ちなみに、BatchNorm処理は試したのですが結局使っていません(コメントアウトしています)。

import gym
import math
import random
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from collections import namedtuple
from itertools import count
from PIL import Image

import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
import torchvision.transforms as T

from pyvirtualdisplay import Display
from scipy.interpolate import interp1d

display = Display(visible=0, size=(1024, 768))
display.start()
import os
#pyvirtualdisplayの仕様変更で次の文はエラーになるのでコメントアウトします(2021/6/5修正)
#os.environ["DISPLAY"] = ":" + str(display.display) + "." + str(display.screen)

#Bipedalwalkerのv2は無くなったのでv3に変更(2021/6/5)
env = gym.make('BipedalWalker-v3')
#env = gym.make('MountainCarContinuous-v0')

# if gpu is to be used
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

#リプレイメモリー
Transition = namedtuple('Transition',
                        ('state', 'action', 'next_state', 'reward','done'))

class ReplayMemory(object):
    def __init__(self, capacity):
        self.capacity = capacity
        self.memory = []
        self.position = 0

    def push(self, *args):
        """Saves a transition."""
        if len(self.memory) < self.capacity:
            self.memory.append(None)
        self.memory[self.position] = Transition(*args)
        self.position = (self.position + 1) % self.capacity

    def sample(self, batch_size):
        return random.sample(self.memory, batch_size)

    def __len__(self):
        return len(self.memory)
   
    def reset(self):
        self.memory = []
        self.position = 0   

#Actor(Polycy) Network
class Policy_net(nn.Module):
    def __init__(self, NumObs, NumAct,NumHidden):
        super(Policy_net, self).__init__()
        self.ln1 = nn.Linear(NumObs, NumHidden)
        self.ln2 = nn.Linear(NumHidden, NumHidden)
        self.ln3 = nn.Linear(NumHidden, NumHidden)
        self.ln4 = nn.Linear(NumHidden, NumAct)
        self.ln4.weight.data.uniform_(-3e-3, 3e-3)

    def forward(self, x):
        x = F.relu(self.ln1(x))
        x = F.relu(self.ln2(x))
        x = F.relu(self.ln3(x))
        y = torch.tanh(self.ln4(x))
        return y

#Critic(Q) Network Q値の計算
class Q_net(nn.Module):
    def __init__(self, NumObs, NumAct,NumHidden1,NumHidden2):
        super(Q_net, self).__init__()
        self.ln1 = nn.Linear(NumObs, NumHidden1)
        self.ln2 = nn.Linear(NumHidden1+NumAct, NumHidden2)
        self.ln3 = nn.Linear(NumHidden2, NumHidden2)
        self.ln4 = nn.Linear(NumHidden2, 1)
        self.ln4.weight.data.uniform_(-3e-3, 3e-3)
        self.bn1 = nn.BatchNorm1d(num_features=NumHidden2)

    def forward(self, state, action):
        x = F.leaky_relu(self.ln1(state))
        x = torch.cat((x, action), dim=1)
#        x = F.leaky_relu(self.bn1(self.ln2(x)))
        x = F.leaky_relu(self.ln2(x))
        x = F.leaky_relu(self.ln3(x))
        y = self.ln4(x)
        return y

def optimize_model():
    if len(memory) < BATCH_SIZE:
        return
    transitions = memory.sample(BATCH_SIZE)
    batch = Transition(*zip(*transitions))

    #batchを一括して処理するためtensorのリストを2次元テンソルに変換
    state_batch = torch.stack(batch.state,dim=0)
    action_batch = torch.stack(batch.action,dim=0)
    reward_batch = torch.stack(batch.reward,dim=0)
    next_state_batch = torch.stack(batch.next_state,dim=0)
    done_batch = torch.stack(batch.done,dim=0)

    #Qネットのloss関数計算
    next_actions = policy_target_net(next_state_batch)
    next_Q_values =q_target_net(next_state_batch ,next_actions)
    expected_Q_values = (next_Q_values * GAMMA)*(1.0-done_batch) + reward_batch

    Q_values =  q_net(state_batch ,action_batch)

    #Qネットのloss関数
    q_loss = F.mse_loss(Q_values,expected_Q_values)

    #Qネットの学習
    q_optimizer.zero_grad()
    #誤差逆伝搬
    q_loss.backward()
    #重み更新
    q_optimizer.step()

    #policyネットのloss関数
    actions = policy_net(state_batch)
    p_loss = -q_net(state_batch,actions).mean()

    #policyネットの学習
    p_optimizer.zero_grad()
    #誤差逆伝搬
    p_loss.backward()
    #重み更新
    p_optimizer.step()

    tau = 0.001
    #targetネットのソフトアップデート
    #学習の収束を安定させるためゆっくり学習するtarget netを作りloss関数の計算に使う。
    #学習後のネット重み×tau分を反映させ、ゆっくり追従させる。
    for target_param, local_param in zip(policy_target_net.parameters(), policy_net.parameters()):
      target_param.data.copy_(tau*local_param.data + (1.0-tau)*target_param.data)

    for target_param, local_param in zip(q_target_net.parameters(), q_net.parameters()):
      target_param.data.copy_(tau*local_param.data + (1.0-tau)*target_param.data)

#action,observationの要素数取得
n_actions = env.action_space.shape[0]
n_observations = env.observation_space.shape[0]
print('num actions,observations',n_actions,n_observations )

#ネットワーク
#policyネットとそのtargetネット
policy_net = Policy_net(n_observations, n_actions,128).to(device)
policy_target_net = Policy_net(n_observations, n_actions,128).to(device)
#Qネットとそのtargetネット
q_net = Q_net(n_observations, n_actions, 64, 128).to(device)
q_target_net = Q_net(n_observations, n_actions, 64, 128).to(device)
#学習用optimizer生成
p_optimizer = optim.Adam(policy_net.parameters(),lr=1e-3)
q_optimizer = optim.Adam(q_net.parameters(),lr=3e-3, weight_decay=0.0001)
#学習用リプレイメモリ生成
memory = ReplayMemory(1000000)
memory.reset()
#ノイズの大きさ計算用、最初は大きくして学習が進んだら小さくするためのカウンタ
steps_done = 0

6.学習

 学習します。1000エピソードごとにネットワークをGoogleドライブに保存します。ファイル名はmodelxxxxtで、xxxxはエピソード数です。

 教師の影響を徐々に減衰するようにしてあって、減衰の仕方はSIGMA_T_START、SIGMA_T_END、SIGMA_T_DECAYで調整します。行動actionは次の式で計算します。

action += sigma_t * (action1 - action)

 sigma_tがそのEPSODEにおける減衰率、action1は教師値、actionはpolicyネットワークの予想する最適行動です。教師とpolicyネットワークの予想値との差を、sigma_tの割合だけ、actionを教師側に近い値に修正しています。教師が手を取って介助しているイメージです。

 ノイズはそのままだと教師の影響が薄れてしまうため若干小さめにしました。

#角度制御ゲイン
 #P(比例)項
pgain = np.array([4,4,4,4], dtype = 'float')
 #D(微分)項
dgain = np.array([0.2,0.2,0.2,0.2], dtype = 'float')

#脚の指令角スケジュール
xsch = np.array([0,1,2,3,4,5,6,7,8,9,10],dtype = 'float')
yhip1  = np.array([ 0.021, -0.39, -0.58, -0.54, 0.45,    0.18, 0.88,  0.78, 0.60,  0.20, 0.021],dtype = 'float')
yknee1 = np.array([  0.17,  0.85,  0.95,  0.92, -0.38, -0.05,  -0.34,  -0.090,  0.074, 0.043,  0.17],dtype = 'float')
yhip2  = np.array([ 0.18,  0.88,  0.78,  0.60,  0.20, 0.021, -0.39, -0.58, -0.54, 0.45,    0.18],dtype = 'float')
yknee2 = np.array([ -0.05,-0.34,-0.090,  0.074,  0.043 ,  0.17,  0.85, 0.95,  0.92, -0.38, 0.05],dtype = 'float')
numx = xsch.shape[0]

hip1demmand = interp1d(xsch, yhip1)
knee1demmand = interp1d(xsch, yknee1)
hip2demmand = interp1d(xsch, yhip2)
knee2demmand = interp1d(xsch, yknee2)

legangle = np.array([0,0,0,0], dtype = 'float')
leganglespeed = np.array([0,0,0,0], dtype = 'float')
legangledemand = np.array([0,0,0,0], dtype = 'float')

speed = 4.2

#BATCH_SIZE = 128
BATCH_SIZE = 128
#Qネット学習時の報酬の割引率
GAMMA = 0.98
#steps_doneに対するノイズの減少係数
SIGMA_START = 0.5 #最初
SIGMA_END = 0.3 #最後
SIGMA_DECAY = 800 #この回数で約30%まで減衰

#steps_doneに対する教師の減少係数
SIGMA_T_START = 0.4 #最初
SIGMA_T_END = 0.02 #最後
SIGMA_T_DECAY = 3000 #この回数で約30%まで減衰

#50 フレーム/sec グラフ描画用
FPS = 50
#学習数。num_epsode=200を変更する。
num_episodes = 7000

for i_episode in range(num_episodes):
    #whileループ内初期化
    x = 0.0
    observation = env.reset()
    done = False
    reward_total = 0.0
    t = 0
    #グラフ用データ保存領域初期化
    frames = []
    observations = []
    actions = []
    ts = []
    for i in range(n_observations):
      observations.append([])
    for i in range(n_actions):
      actions.append([])

    #NOISE
    noise =np.array([random.uniform(-0.5,0.5) for i in range(n_actions)],dtype = np.float)
    theta = 0.08
    sigma = SIGMA_END + (SIGMA_START - SIGMA_END) * math.exp(-1. * steps_done / SIGMA_DECAY)
    #TEACHER
    sigma_t = SIGMA_T_END + (SIGMA_T_START - SIGMA_T_END) * math.exp(-1. * steps_done / SIGMA_T_DECAY)
    steps_done += 1
    display_on = 0
    
    while not done and t < 700:
        #指令値生成------------------------------------
        x += speed/FPS
        if x>10.0 :x-=10
        legangledemand[0] = hip1demmand(x)
        legangledemand[1] = knee1demmand(x)
        legangledemand[2] = hip2demmand(x)
        legangledemand[3] = knee2demmand(x)

        #脚の角度 PD制御
        action1 = pgain * (legangledemand - legangle) - dgain * leganglespeed
        action1 = np.clip(action1, -0.6, 0.6)

        #----------------------------------------------
        with torch.no_grad():
          action = policy_net(torch.from_numpy(observation).float().to(device))
          action = action.cpu().data.numpy()

        if (i_episode != num_episodes -1 ):
          action += sigma_t * (action1 - action)

        #最後は純粋にネットワークのデータを取得するためノイズ無し-------------------
        if (i_episode != num_episodes -1 ):
          #OUNoise
          noise = noise - theta * noise + sigma * np.array([random.uniform(-1.0,1.0) for i in range(len(noise))])
          action += noise

        action = np.clip(action, -1, 1)

        #物理モデル1ステップ---------------------------
        next_observation, reward, done, info = env.step(action)
        reward_total = reward_total + reward

        #脚の角度,角速度,胴体角取り出し----------------
        legangle[0] = observation[4] #HIP1
        leganglespeed[0] = observation[5]
        legangle[1] = observation[6] #KNEE1
        leganglespeed[1] = observation[7]
        legangle[2] = observation[9] #HIP2
        leganglespeed[2] = observation[10]
        legangle[3] = observation[11] #KNEE2
        leganglespeed[3] = observation[11]

        #学習用にメモリに保存--------------------------
        tensor_observation = torch.tensor(observation,device=device,dtype=torch.float)
        tensor_action = torch.tensor(action,device=device,dtype=torch.float)
        tensor_next_observation = torch.tensor(next_observation,device=device,dtype=torch.float)
        tensor_reward = torch.tensor([reward],device=device,dtype=torch.float)
        tensor_done =  torch.tensor([done],device=device,dtype=torch.float)
        memory.push(tensor_observation, tensor_action, tensor_next_observation, tensor_reward,tensor_done)

        #observation(state)更新------------------------
        observation = next_observation

        #学習してpolicy_net,target_neを更新
        optimize_model()
        #データ保存------------------------------------
        if (i_episode == num_episodes -1) and (display_on == 1) :
          #動画
          frames.append(env.render(mode = 'rgb_array'))
          #グラフ
          for i in range(n_observations):
            observations[i].append(observation[i])
          for i in range(n_actions):
            actions[i].append(action[i])
          ts.append(t/FPS)

        #時間を進める----------------------------------
        t += 1
        # end while loop ------------------------------
    if (steps_done % 10 == 0) or (i_episode == num_episodes -1 ):
#      print('episode,reward=',i_episode,reward_total)
      print('episode,reward=',steps_done,reward_total)
    if (steps_done % 1000 == 0):
      torch.save({
        'epoch':steps_done,
        'policy':policy_net.state_dict(),
        'policy_target':policy_target_net.state_dict(),
        'p_optimizer':p_optimizer.state_dict(),
        'q':q_net.state_dict(),
        'q_target':q_target_net.state_dict(),
        'q_optimizer':q_optimizer.state_dict(),
      },'/content/drive/My Drive/model' + str(steps_done) + 't' )
  # end for loop --------------------------------------

7.ネットワーク読み込み

  読み込みたいネットワークにあわせて赤字は都度変更してください。

checkpoint = torch.load('/content/drive/My Drive/model4000t')
steps_done = checkpoint['epoch']
policy_net.load_state_dict(checkpoint['policy'])
policy_target_net.load_state_dict(checkpoint['policy_target'])
p_optimizer.load_state_dict(checkpoint['p_optimizer'])
q_net.load_state_dict(checkpoint['q'])
q_target_net.load_state_dict(checkpoint['q_target'])
q_optimizer.load_state_dict(checkpoint['q_optimizer'])

8.学習したネットワークを動かして動画を保存

 6.学習のソースを一部変更して動かします。

①num_episodesを1にする。

 #学習数。num_epsode=200を変更する。
 num_episodes = 1

②画像保存スイッチを入れる。

 display_on = 1

 ③実行します。

9.画像を表示

import matplotlib.pyplot as plt
import matplotlib.animation
import numpy as np
from IPython.display import HTML

fig2 = plt.figure(figsize=(6, 6))
plt.axis('off')

images = []
for f in frames:
  image = plt.imshow(f)
  images.append([image])
ani = matplotlib.animation.ArtistAnimation(fig2, images, interval=30, repeat_delay=1)

HTML(ani.to_jshtml())

10.まとめ

 教師をつけて試しにDDPGで学習をしてみました。結果は次の通りです。

  • 学習したpolicyネットワークは最初は教師の影響を受けて、教師に近い動き(大股歩き)をしようとしました。
  • 教師の影響が強いと、もともと教師自体がすぐにつまづくためか、学習したpolicyもすぐにつまづいてしまいました(10000回も学習すると教師よりは多少ましでしたが・・)
  • SIGMA_T_ENDを0.2→0.02に変更して教師の影響を小さくすると、徐々にちょこちょこ走りになっていきました。
  • 結局、教師無しと比較して、学習回数が減ったとか、大股で歩けるようになったとかいう具体的なメリットは得られませんでした

 今回はあまりいい結果ではありませんでしたが、DDPG学習に教師の影響を受けさせられることを確かめられたので、どこかで使えることもあるかもしれません。ということで前向きに考えることにします!(大分試すのに時間をとられましたが・・・)。

おしまい

 BipedalWalkerの記事はこれでおしまいです!。今回の記事は長く続きましたが、最後までお付き合い頂きありがとうございました。

 強化学習 カテゴリーの記事一覧 - やってみた!

Open AI Gym Box2D BipedalWalkerをColaboratoryで動かしてみる(6)

 前回はDDPGをざくっと解説してみました。

 今回はDDPGでBipedalWalkerを学習させてみます。

1.BipedalWalker-v2の報酬について

 Open AI Gymのgithubサイトにあるソースリスト

 https://github.com/openai/gym

 gym-envs-box2d-bipedal_walker.pyの関数stepを見ると報酬は次のようになっています。

  • ゲームオーバ(体が地面に接地、もしくは領域外まで後ろに下がる)
    reward = -100
  • 上記以外
    shaping = 130*pos[0]/SCALE-5.0*abs(state[0])
    reward = shaping - prev_shaping - 0.00035 * MOTORS_TORQUE * np.clip(np.abs(a), 0, 1)

    進んだ距離:前に進みゴール位置で報酬=300
    姿勢   :水平で報酬0、傾いている分だけ減点
    モータ  :モータトルクを使った分だけ報酬を減点

2.DDPGを試してみる

 動かしてみる(4)で紹介したソースリストのコメントアウトを変更するとBipedalWalkerを動かせます。

(1)環境を変更 

変更前
#env = gym.make('BipedalWalker-v2')
env = gym.make('MountainCarContinuous-v0')

変更後
env = gym.make('BipedalWalker-v2')
#env = gym.make('MountainCarContinuous-v0')

(2)エピソード回数、獲得報酬表示

変更前

        # end while loop ------------------------------
    print('episode,reward=',i_episode,reward_total)

変更後

        # end while loop ------------------------------
    if (steps_done % 10 == 0) or (i_episode == num_episodes -1 ):
#      print('episode,reward=',i_episode,reward_total)
      print('episode,reward=',steps_done,reward_total)

・学習回数が多いので毎回表示ではなく10回に1回表示に変更。
・学習の様子を見て追加で学習させた時の総エピソード数を知るため、実行のたびに0からカウントされるi_epsodeのかわりに、リセットされないsteps_doneを表示。

(3)パラメータ変更

#学習用リプレイメモリ生成
memory = ReplayMemory(1000000)


#BATCH_SIZE = 128
BATCH_SIZE = 128
#Qネット学習時の報酬の割引率
GAMMA = 0.98
#steps_doneに対するノイズの減少係数
SIGMA_START = 1.0 #最初
SIGMA_END = 0.3 #最後
SIGMA_DECAY = 800 #この回数で約30%まで減衰

#50 フレーム/sec グラフ描画用
FPS = 50
#学習数。num_epsode=200を変更する。
num_episodes = 500

(4)最大ステップ数削減
 計算時間がかかるので最大ステップ数を1400→700に変更しました。

while not done and t < 700:
#指令値生成------------------------------------

(5)結果
 学習を何回か繰り返しました。うまくいったケースでは学習数約1000回でぎこちなくスキップして進んでいきました。ですが、さらに学習を進め1500回では、すぐに転んでしまいました。

1000回

f:id:akifukka:20191226153720g:plain

episode,reward= 1000 35.9231961772401

 

2000回
反対側の脚も使うようになってきてちょっと期待できそうです。

f:id:akifukka:20191226160713g:plain

episode,reward= 2000 -51.45643730989602

 

2500回
ちょこちょこ走り

f:id:akifukka:20191226163953g:plain

episode,reward= 2500 208.7962870610689

 3000回、3500回
ちょこちょこ走りのままですが、すぐに転んでしまいました。

4000回
ちょこちょこ結構走りましたが、2500回の総報酬には届かず。

6000回

f:id:akifukka:20191226214324g:plain

episode,reward= 6000 256.1380908660748

7000回
 while文のループを増やしたらゴールまで行けました。6000回の時にみられたおっとっとの動きもなく、安定したちょこちょこ走りでした。

f:id:akifukka:20191226215201p:plain
episode,reward= 7001 306.013638338901

(6)一応ソースリストを貼り付けておきます。

 ライブラリのインストールはOpen AI Gym Box2D BipedalWalkerをColaboratoryで動かしてみる(4)を参照してください。

学習前の処理

import gym
import math
import random
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from collections import namedtuple
from itertools import count
from PIL import Image

import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
import torchvision.transforms as T

from pyvirtualdisplay import Display
from scipy.interpolate import interp1d

display = Display(visible=0, size=(1024, 768))
display.start()
import os
#pyvirtualdisplayの仕様変更で次の文はエラーになるのでコメントアウトします(2021/6/5修正)
#os.environ["DISPLAY"] = ":" + str(display.display) + "." + str(display.screen)

#Bipedalwalkerのv2は無くなったのでv3に変更(2021/6/5)
env = gym.make('BipedalWalker-v3')
#env = gym.make('MountainCarContinuous-v0')

# if gpu is to be used
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

#リプレイメモリー
Transition = namedtuple('Transition',
                        ('state', 'action', 'next_state', 'reward','done'))

class ReplayMemory(object):
    def __init__(self, capacity):
        self.capacity = capacity
        self.memory = []
        self.position = 0

    def push(self, *args):
        """Saves a transition."""
        if len(self.memory) < self.capacity:
            self.memory.append(None)
        self.memory[self.position] = Transition(*args)
        self.position = (self.position + 1) % self.capacity

    def sample(self, batch_size):
        return random.sample(self.memory, batch_size)

    def __len__(self):
        return len(self.memory)
   
    def reset(self):
        self.memory = []
        self.position = 0   

#Actor(Polycy) Network
class Policy_net(nn.Module):
    def __init__(self, NumObs, NumAct,NumHidden):
        super(Policy_net, self).__init__()
        self.ln1 = nn.Linear(NumObs, NumHidden)
        self.ln2 = nn.Linear(NumHidden, NumHidden)
        self.ln3 = nn.Linear(NumHidden, NumHidden)
        self.ln4 = nn.Linear(NumHidden, NumAct)
        self.ln4.weight.data.uniform_(-3e-3, 3e-3)

    def forward(self, x):
        x = F.relu(self.ln1(x))
        x = F.relu(self.ln2(x))
        x = F.relu(self.ln3(x))
        y = torch.tanh(self.ln4(x))
        return y

#Critic(Q) Network Q値の計算
class Q_net(nn.Module):
    def __init__(self, NumObs, NumAct,NumHidden1,NumHidden2):
        super(Q_net, self).__init__()
        self.ln1 = nn.Linear(NumObs, NumHidden1)
        self.ln2 = nn.Linear(NumHidden1+NumAct, NumHidden2)
        self.ln3 = nn.Linear(NumHidden2, NumHidden2)
        self.ln4 = nn.Linear(NumHidden2, 1)
        self.ln4.weight.data.uniform_(-3e-3, 3e-3)
        self.bn1 = nn.BatchNorm1d(num_features=NumHidden2)

    def forward(self, state, action):
        x = F.leaky_relu(self.ln1(state))
        x = torch.cat((x, action), dim=1)
#        x = F.leaky_relu(self.bn1(self.ln2(x)))
        x = F.leaky_relu(self.ln2(x))
        x = F.leaky_relu(self.ln3(x))
        y = self.ln4(x)
        return y

def optimize_model():
    if len(memory) < BATCH_SIZE:
        return
    transitions = memory.sample(BATCH_SIZE)
    batch = Transition(*zip(*transitions))

    #batchを一括して処理するためtensorのリストを2次元テンソルに変換
    state_batch = torch.stack(batch.state,dim=0)
    action_batch = torch.stack(batch.action,dim=0)
    reward_batch = torch.stack(batch.reward,dim=0)
    next_state_batch = torch.stack(batch.next_state,dim=0)
    done_batch = torch.stack(batch.done,dim=0)

    #Qネットのloss関数計算
    next_actions = policy_target_net(next_state_batch)
    next_Q_values =q_target_net(next_state_batch ,next_actions)
    expected_Q_values = (next_Q_values * GAMMA)*(1.0-done_batch) + reward_batch

    Q_values =  q_net(state_batch ,action_batch)

    #Qネットのloss関数
    q_loss = F.mse_loss(Q_values,expected_Q_values)

    #Qネットの学習
    q_optimizer.zero_grad()
    #誤差逆伝搬
    q_loss.backward()
    #重み更新
    q_optimizer.step()

    #policyネットのloss関数
    actions = policy_net(state_batch)
    p_loss = -q_net(state_batch,actions).mean()

    #policyネットの学習
    p_optimizer.zero_grad()
    #誤差逆伝搬
    p_loss.backward()
    #重み更新
    p_optimizer.step()

    tau = 0.001
    #targetネットのソフトアップデート
    #学習の収束を安定させるためゆっくり学習するtarget netを作りloss関数の計算に使う。
    #学習後のネット重み×tau分を反映させ、ゆっくり追従させる。
    for target_param, local_param in zip(policy_target_net.parameters(), policy_net.parameters()):
      target_param.data.copy_(tau*local_param.data + (1.0-tau)*target_param.data)

    for target_param, local_param in zip(q_target_net.parameters(), q_net.parameters()):
      target_param.data.copy_(tau*local_param.data + (1.0-tau)*target_param.data)

#action,observationの要素数取得
n_actions = env.action_space.shape[0]
n_observations = env.observation_space.shape[0]
print('num actions,observations',n_actions,n_observations )

#ネットワーク
#policyネットとそのtargetネット
policy_net = Policy_net(n_observations, n_actions,128).to(device)
policy_target_net = Policy_net(n_observations, n_actions,128).to(device)
#Qネットとそのtargetネット
q_net = Q_net(n_observations, n_actions, 64, 128).to(device)
q_target_net = Q_net(n_observations, n_actions, 64, 128).to(device)
#学習用optimizer生成
p_optimizer = optim.Adam(policy_net.parameters(),lr=1e-3)
q_optimizer = optim.Adam(q_net.parameters(),lr=3e-3, weight_decay=0.0001)
#学習用リプレイメモリ生成
memory = ReplayMemory(1000000)
memory.reset()
#ノイズの大きさ計算用、最初は大きくして学習が進んだら小さくするためのカウンタ
steps_done = 0

学習する。エピソード1000になっているので、何度か繰り返すか、値を変更して実行してください。

#BATCH_SIZE = 128
BATCH_SIZE = 128
#Qネット学習時の報酬の割引率
GAMMA = 0.98
#steps_doneに対するノイズの減少係数
SIGMA_START = 1.0 #最初
SIGMA_END = 0.3 #最後
SIGMA_DECAY = 800 #この回数で約30%まで減衰

#50 フレーム/sec グラフ描画用
FPS = 50
#学習数。num_epsode=200を変更する。
num_episodes = 1000

for i_episode in range(num_episodes):
    #whileループ内初期化
    observation = env.reset()
    done = False
    reward_total = 0.0
    t = 0
    #グラフ用データ保存領域初期化
    frames = []
    observations = []
    actions = []
    ts = []
    for i in range(n_observations):
      observations.append([])
    for i in range(n_actions):
      actions.append([])

    noise =np.array([random.uniform(-0.5,0.5) for i in range(n_actions)],dtype = np.float)
    theta = 0.08
    sigma = SIGMA_END + (SIGMA_START - SIGMA_END) * math.exp(-1. * steps_done / SIGMA_DECAY)
    steps_done += 1

    while not done and t < 700:
        #指令値生成------------------------------------
        sample = random.random()

        with torch.no_grad():
          action = policy_net(torch.from_numpy(observation).float().to(device))
          action = action.cpu().data.numpy()

        #最後は純粋にネットワークのデータを取得するためノイズ無し-------------------
        if (i_episode != num_episodes -1 ):
          #OUNoise
          noise = noise - theta * noise + sigma * np.array([random.uniform(-1.0,1.0) for i in range(len(noise))])
          action += noise

        action = np.clip(action, -1, 1)

        #物理モデル1ステップ---------------------------
        next_observation, reward, done, info = env.step(action)
        reward_total = reward_total + reward

        #学習用にメモリに保存--------------------------
        tensor_observation = torch.tensor(observation,device=device,dtype=torch.float)
        tensor_action = torch.tensor(action,device=device,dtype=torch.float)
        tensor_next_observation = torch.tensor(next_observation,device=device,dtype=torch.float)
        tensor_reward = torch.tensor([reward],device=device,dtype=torch.float)
        tensor_done =  torch.tensor([done],device=device,dtype=torch.float)
        memory.push(tensor_observation, tensor_action, tensor_next_observation, tensor_reward,tensor_done)

        #observation(state)更新------------------------
        observation = next_observation

        #学習してpolicy_net,target_neを更新
        optimize_model()
        #データ保存------------------------------------
        if i_episode == num_episodes -1 :
          #動画
          frames.append(env.render(mode = 'rgb_array'))
          #グラフ
          for i in range(n_observations):
            observations[i].append(observation[i])
          for i in range(n_actions):
            actions[i].append(action[i])
          ts.append(t/FPS)

        #時間を進める----------------------------------
        t += 1
        # end while loop ------------------------------
    if (steps_done % 10 == 0) or (i_episode == num_episodes -1 ):
#      print('episode,reward=',i_episode,reward_total)
      print('episode,reward=',steps_done,reward_total)
  # end for loop --------------------------------------

 動画を表示する

import matplotlib.pyplot as plt
import matplotlib.animation
import numpy as np
from IPython.display import HTML

fig2 = plt.figure(figsize=(6, 6))
plt.axis('off')

images = []
for f in frames:
  image = plt.imshow(f)
  images.append([image])
ani = matplotlib.animation.ArtistAnimation(fig2, images, interval=30, repeat_delay=1)

HTML(ani.to_jshtml())

 BipedalWalker、とりあえず(思ったよりも)何とか動きました。

つづく
 強化学習 カテゴリーの記事一覧 - やってみた!

Open AI Gym Box2D BipedalWalkerをColaboratoryで動かしてみる(5)

  前回はDDPG(Deep Deterministic Policy Gradient)でMountainCarContinuous挑戦し、無事学習して山登りに成功しました。(BipedalWalkerは手強いので後回しです・・・)

 今回は中身について、ざっくりですが解説してみます。

 1.DDPG(Deep Deterministic Policy Gradient)

 f:id:akifukka:20191223201154j:plain

 DQNでは、ある環境の時に行動a1,a2・・・を取った時の価値Q1,Q2・・(未来に得られる報酬の合計)を推定するQネットワークを学習します。ここで環境やQxは連続値を取ります。学習が完了したら、Q1,Q2・・・を比較して最も大きな価値の行動axを選択します。ネットワークの構造上、行動a1,a2・・・は離散的になります。点数を限りなく増やせば連続値を模擬できますが、あまり現実的ではありません。

 これに対しDDPGでは行動aもQ値を計算する入力にします。ここでの行動aは連続値です。ただ、このままだと学習後にどの行動aを選ぶかといった時に、aの値を変化させて価値Qが最大になるところを探し出す必要があります。

 じゃあということで、最初から価値Qが最大になる行動aを推定するpolicyネットワークを追加しちゃえばというのがDDPGです。

 policyネットワークはQネットワークの値が最大になる様、Qネットワークを損失関数に使って学習させます。

2.Qネットワーク(Q(s,a))の学習

 'Bellman equation'という考え方を使います。

f:id:akifukka:20191224201708j:plain

 価値関数Qを、これからもらえる報酬の総和で定義します。Qの引数は1項の図を見てわかるようにobservationの状態sと行動aです。

 ただし、そのまま報酬を未来永劫に渡って足し続けるとQが無限大に発散してしまうので、報酬に割引率を乗じて意味のある数字に収束するようにします(①)。すると、Qは次のステップQ(st+1,at+1)と割引率γを使って書くことができます(②)。

 次に左辺=0となるように変形します(③)。で、2乗して損失関数Lを定義して(④)、損失関数が最小になるように学習すればQを学習できることになります。

 ちなみに時間tでゴールした場合は、それ以降は報酬をもらえない、すなわちQt+1=0になるはずなので、その模擬も含めて損失関数は最終的に次のように書けます。

f:id:akifukka:20191223221032j:plain

 赤字はtarget関数です。これは価値関数Qをそのまま使うと、学習で自身の更新の影響を受け安定して収束しなくなるのを防ぐため、価値関数Qより少し遅れてゆっくり変化するようにしたものです。

 割引率は0.999等、1に近い値を使います。1に近ければ近いほど遠い将来の報酬の影響を受け、値が小さければ比較的短い時間範囲の報酬までしか影響を受けません。

 Qネットワークは上記の損失関数(loss関数)が最小になるように学習させます。
環境st、その時に取った行動at、その結果の環境st+1は実際に1ステップ動かした時の値(経験値)を使えます。1ステップの経験値からat+1の情報は得られないので、後で説明する行動を予測するpolicyネットワークを使って計算します。この時使用するpolicyネットワークも実際のpolicyネットワークに遅れて学習するtarget関数です。

 ソースリストのQネットワークの学習はdef optimize_model()中で行っています。該当する箇所を次に示します。

#Qネットのloss関数計算
next_actions = policy_target_net(next_state_batch)
next_Q_values =q_target_net(next_state_batch ,next_actions)
expected_Q_values = (next_Q_values * GAMMA)*(1.0-done_batch) + reward_batch

Q_values =  q_net(state_batch ,action_batch)

#Qネットのloss関数
q_loss = F.mse_loss(Q_values,expected_Q_values)

#Qネットの学習
q_optimizer.zero_grad()
#誤差逆伝搬
q_loss.backward()
#重み更新
q_optimizer.step()

 F.mse_loss()は差の2乗平均を算出する関数で、q_optimizer.zero_grad()~q_optimizer.step()で学習します。

 なお、複数の経験(ステップ)を一括して学習するため、環境s、行動a等は複数の経験分をstack(テンソルを複数まとめる)しています。

3. policyネットワーク(policy(s))の学習

 policyネットワークはQが最大となる行動を推測するよう学習させます。pytorchは損失関数が最小になるよう学習させるため、Q値に-1を乗じた損失関数を定義、値が最小になるよう学習させます。

  該当する箇所を次に示します。policyネットワーク(Target)はQネットの学習、経験の取得に使用するため、Qネットと交互に学習させ、並行して学習させます。

#policyネットのloss関数
actions = policy_net(state_batch)
p_loss = -q_net(state_batch,actions).mean()

#policyネットの学習
p_optimizer.zero_grad()
#誤差逆伝搬
p_loss.backward()
#重み更新
p_optimizer.step()

4.Target関数のソフトアップデート

 Target関数は前述の通り実際のQネット、policyネットにゆっくり追従するネットワークで、それぞれQネット、policyネットを学習した後、次のように更新します。

tau = 0.001
#targetネットのソフトアップデート
#学習の収束を安定させるためゆっくり学習するtarget netを作りloss関数の計算に使う。
#学習後のネット重み×tau分を反映させ、ゆっくり追従させる。
for target_param, local_param in zip(policy_target_net.parameters(), policy_net.parameters()):
  target_param.data.copy_(tau*local_param.data + (1.0-tau)*target_param.data)

for target_param, local_param in zip(q_target_net.parameters(), q_net.parameters()):
  target_param.data.copy_(tau*local_param.data + (1.0-tau)*target_param.data)

 ネットワークの重みをそのままコピーするのではなく、tau(この場合は1/1000)だけ学習後の値を使い、残りの割合分は現状の値をそのまま使うことで、学習の効果がゆっくり浸透するようにしています(1次遅れと同じイメージ)。

5.リプレイメモリ

 直前の経験だけを使って学習すると、直近のpolicyが出す行動に偏りQネットの学習もその領域に偏ってしまうため、過去の経験をメモリにためておき満遍なく学習するようにしています。リプレイメモリは有限の大きさにしておき、あまりに古くて使えない経験は新しい経験に上書きされます。

6.ノイズ

 経験を取得するのにpolicyの出力する行動を使うと、行動が偏ってしまい、学習もpolicyが出力する行動のごく狭い範囲に留まってしまいます。そのため極所最適解に陥り、他に良い行動があっても発見できまないまま終わってしまいます。

 DDPGでは、より良い行動を探すため、policyが出力する行動にノイズを混ぜて幅広く経験させる様、工夫されています。

7.ネットワーク定義

 今回はQネット、policyネット共に4層、隠れ要素128のネットワークにしてみました。層数、隠れ要素数ともに特にこの値にした理由があるわけではないので、変更して試してみるのもいいと思います。

8.まとめ

 DDPGのアルゴリズムの概要を纏めました。学習を効率良く行うには、報酬の与え方(今回はOpen AI Gymなので変更しない)、割引率、ノイズの大きさ、ネットワーク構造といったあたりが鍵になると思われます。

 次回は再びBipedalWalkerを動かします。

 強化学習 カテゴリーの記事一覧 - やってみた!

Open AI Gym Box2D BipedalWalkerをColaboratoryで動かしてみる(4)

改正2019.12.26

 ソースリスト中でsteps_done +=1の位置をwhileループ(各ステップ計算ループ)からepsodeのforループに移動しました(バグ)。このバグのためSIGMA_DECAYがほとんど効かず、すぐにノイズが小さくなっていました。あわせてSIGMA_DECAYの設定値も150→30に見直しました。

-----------------------------------

 前回は歩行パターンのスケジュールを作り、PD制御で脚の角度を制御し何とか2,3歩歩くことができました。

 ここからはpytorchで深層強化学習を実装していきます。が、BipedalWalkerはあまりに手強いので、手始めにMountainCarContinuousで試してみます。

 主に次のサイトを参考にさせてもらいました。

  実装はDDPG(Deep Deterministic Policy Gradient)です。元の論文はこちら。

 また、こちらにOPEN AIによる解説があります。

Deep Deterministic Policy Gradient — Spinning Up documentation

1.DDPG 

 DDPGは強化学習の手法のひとつDPG(Deterministic Policy Gradient)を発展させ、DeepLerningを適用したものです。元となったDPGの論文はこちら。

Deterministic Policy Gradient Algorithms,Silver et al., 2014

 DDPGの特徴は次の通りです。

  • action(制御の指令値)に連続値(discreteではなく)を扱うことができますDQN(Deep Q Network)の制御指令値はDiscrete値なので応用範囲が限定されますが、DDPGは連続値なので一般的な制御に応用することができます。
  • actor-critic アルゴリズムを採用します。actor(policy)は行動方策を推定し、critic(Q関数)は行動の価値を推定するものです。これらをディープネットワークで定義して同時に学習していきます。

2.動かしてみます

 DDPGの詳細は後回しにして、まずはどんなものか動かしてみます。

(1)Corabolatoryの準備

 Colaboratoryを立ち上げます。「ランタイム」ー「ランタイムのタイプを変更」でハードウェアアクセラレータ「GPU」を選択します。

(2)必要なライブラリのインストール

必要なライブラリをインストールします。

!apt-get update
!pip3 install box2d-py
!pip3 install 'gym[Box2D]'
!pip3 install  'gym[classic_control]'
!apt-get -qq -y install xvfb freeglut3-dev ffmpeg
!pip3 -q install pyglet
!pip3 -q install pyopengl
!pip3 -q install pyvirtualdisplay
!apt-get install x11-utils

(3)ランタイムを再起動します。

 インルトールした内容を反映させるため、一旦ランタイムを再起動します。
「ランタイム」-「ランタイムの再起動」

(4)学習前の処理(リプレイメモリ、ネットワークの定義、初期化)

 次のコードを実行します。

import gym
import math
import random
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from collections import namedtuple
from itertools import count
from PIL import Image

import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
import torchvision.transforms as T

from pyvirtualdisplay import Display
from scipy.interpolate import interp1d

display = Display(visible=0, size=(1024, 768))
display.start()
import os
#pyvirtualdisplayの仕様変更で次の文はエラーになるのでコメントアウトします(2021/6/5変更)
#os.environ["DISPLAY"] = ":" + str(display.display) + "." + str(display.screen)

#Bipedalwalkerのv2は無くなったのでv3に変更(2021/6/5)
#env = gym.make('BipedalWalker-v3')
env = gym.make('MountainCarContinuous-v0')

# if gpu is to be used
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

#リプレイメモリー
Transition = namedtuple('Transition',
                        ('state', 'action', 'next_state', 'reward','done'))

class ReplayMemory(object):
    def __init__(self, capacity):
        self.capacity = capacity
        self.memory = []
        self.position = 0

    def push(self, *args):
        """Saves a transition."""
        if len(self.memory) < self.capacity:
            self.memory.append(None)
        self.memory[self.position] = Transition(*args)
        self.position = (self.position + 1) % self.capacity

    def sample(self, batch_size):
        return random.sample(self.memory, batch_size)

    def __len__(self):
        return len(self.memory)
   
    def reset(self):
        self.memory = []
        self.position = 0   

#Actor(Polycy) Network
class Policy_net(nn.Module):
    def __init__(self, NumObs, NumAct,NumHidden):
        super(Policy_net, self).__init__()
        self.ln1 = nn.Linear(NumObs, NumHidden)
        self.ln2 = nn.Linear(NumHidden, NumHidden)
        self.ln3 = nn.Linear(NumHidden, NumHidden)
        self.ln4 = nn.Linear(NumHidden, NumAct)
        self.ln4.weight.data.uniform_(-3e-3, 3e-3)

    def forward(self, x):
        x = F.relu(self.ln1(x))
        x = F.relu(self.ln2(x))
        x = F.relu(self.ln3(x))
        y = torch.tanh(self.ln4(x))
        return y

#Critic(Q) Network Q値の計算
class Q_net(nn.Module):
    def __init__(self, NumObs, NumAct,NumHidden1,NumHidden2):
        super(Q_net, self).__init__()
        self.ln1 = nn.Linear(NumObs, NumHidden1)
        self.ln2 = nn.Linear(NumHidden1+NumAct, NumHidden2)
        self.ln3 = nn.Linear(NumHidden2, NumHidden2)
        self.ln4 = nn.Linear(NumHidden2, 1)
        self.ln4.weight.data.uniform_(-3e-3, 3e-3)

    def forward(self, state, action):
        x = F.leaky_relu(self.ln1(state))
        x = torch.cat((x, action), dim=1)
        x = F.leaky_relu(self.ln2(x))
        x = F.leaky_relu(self.ln3(x))
        y = self.ln4(x)
        return y

def optimize_model():
    if len(memory) < BATCH_SIZE:
        return
    transitions = memory.sample(BATCH_SIZE)
    batch = Transition(*zip(*transitions))

    #batchを一括して処理するためtensorのリストを2次元テンソルに変換
    state_batch = torch.stack(batch.state,dim=0)
    action_batch = torch.stack(batch.action,dim=0)
    reward_batch = torch.stack(batch.reward,dim=0)
    next_state_batch = torch.stack(batch.next_state,dim=0)
    done_batch = torch.stack(batch.done,dim=0)

    #Qネットのloss関数計算
    next_actions = policy_target_net(next_state_batch)
    next_Q_values =q_target_net(next_state_batch ,next_actions)
    expected_Q_values = (next_Q_values * GAMMA)*(1.0-done_batch) + reward_batch

    Q_values =  q_net(state_batch ,action_batch)

    #Qネットのloss関数
    q_loss = F.mse_loss(Q_values,expected_Q_values)

    #Qネットの学習
    q_optimizer.zero_grad()
    #誤差逆伝搬
    q_loss.backward()
    #重み更新
    q_optimizer.step()

    #policyネットのloss関数
    actions = policy_net(state_batch)
    p_loss = -q_net(state_batch,actions).mean()

    #policyネットの学習
    p_optimizer.zero_grad()
    #誤差逆伝搬
    p_loss.backward()
    #重み更新
    p_optimizer.step()

    tau = 0.001
    #targetネットのソフトアップデート
    #学習の収束を安定させるためゆっくり学習するtarget netを作りloss関数の計算に使う。
    #学習後のネット重み×tau分を反映させ、ゆっくり追従させる。
    for target_param, local_param in zip(policy_target_net.parameters(), policy_net.parameters()):
      target_param.data.copy_(tau*local_param.data + (1.0-tau)*target_param.data)

    for target_param, local_param in zip(q_target_net.parameters(), q_net.parameters()):
      target_param.data.copy_(tau*local_param.data + (1.0-tau)*target_param.data)

#action,observationの要素数取得
n_actions = env.action_space.shape[0]
n_observations = env.observation_space.shape[0]
print('num actions,observations',n_actions,n_observations )

#ネットワーク
#policyネットとそのtargetネット
policy_net = Policy_net(n_observations, n_actions,128).to(device)
policy_target_net = Policy_net(n_observations, n_actions,128).to(device)
#Qネットとそのtargetネット
q_net = Q_net(n_observations, n_actions, 64, 128).to(device)
q_target_net = Q_net(n_observations, n_actions, 64, 128).to(device)
#学習用optimizer生成
p_optimizer = optim.Adam(policy_net.parameters(),lr=1e-3)
q_optimizer = optim.Adam(q_net.parameters(),lr=3e-3, weight_decay=0.0001)
#学習用リプレイメモリ生成
memory = ReplayMemory(10000)
memory.reset()
#ノイズの大きさ計算用、最初は大きくして学習が進んだら小さくするためのカウンタ
steps_done = 0

(5)学習します

 所要時間は10分くらいでしょうか。表示されるrewardは報酬の合計値で、山登り成功の場合は80~95点くらいになります。学習中は制御指令値にノイズを加算しているので完全な実力値とは言えません。epsodeの最後の回はノイズ無しなので、ほぼ実力値になります。

 学習epsodeが200回の場合、0から始めて199回目のrewardが実力相当の値です。学習数は必要に応じて変更してください。

#BATCH_SIZE = 128
BATCH_SIZE = 128
#Qネット学習時の報酬の割引率
GAMMA = 0.999
#step_doneに対するノイズの減少係数
SIGMA_START = 1.0 #最初
SIGMA_END = 0.3 #最後
SIGMA_DECAY = 30 #この回数で約30%まで減衰

#50 フレーム/sec グラフ描画用
FPS = 50
#学習数。num_epsode=200を変更する。
num_episodes = 200

for i_episode in range(num_episodes):
    #whileループ内初期化
    observation = env.reset()
    done = False
    reward_total = 0.0
    t = 0
    #グラフ用データ保存領域初期化
    frames = []
    observations = []
    actions = []
    ts = []
    for i in range(n_observations):
      observations.append([])
    for i in range(n_actions):
      actions.append([])

    noise =np.array([random.uniform(-0.5,0.5) for i in range(n_actions)],dtype = np.float)
    theta = 0.08
    sigma = SIGMA_END + (SIGMA_START - SIGMA_END) * math.exp(-1. * steps_done / SIGMA_DECAY)
    steps_done += 1
while not done and t < 1400: #指令値生成------------------------------------ sample = random.random() with torch.no_grad(): action = policy_net(torch.from_numpy(observation).float().to(device)) action = action.cpu().data.numpy() #最後は純粋にネットワークのデータを取得するためノイズ無し------------------- if (i_episode != num_episodes -1 ): #OUNoise noise = noise - theta * noise + sigma * np.array([random.uniform(-1.0,1.0) for i in range(len(noise))]) action += noise action = np.clip(action, -1, 1) #物理モデル1ステップ--------------------------- next_observation, reward, done, info = env.step(action) reward_total = reward_total + reward #学習用にメモリに保存-------------------------- tensor_observation = torch.tensor(observation,device=device,dtype=torch.float) tensor_action = torch.tensor(action,device=device,dtype=torch.float) tensor_next_observation = torch.tensor(next_observation,device=device,dtype=torch.float) tensor_reward = torch.tensor([reward],device=device,dtype=torch.float) tensor_done = torch.tensor([done],device=device,dtype=torch.float) memory.push(tensor_observation, tensor_action, tensor_next_observation, tensor_reward,tensor_done) #observation(state)更新------------------------ observation = next_observation #学習してpolicy_net,target_neを更新 optimize_model() #データ保存------------------------------------ if i_episode == num_episodes -1 : #動画 frames.append(env.render(mode = 'rgb_array')) #グラフ for i in range(n_observations): observations[i].append(observation[i]) for i in range(n_actions): actions[i].append(action[i]) ts.append(t/FPS) #時間を進める---------------------------------- t += 1 # end while loop ------------------------------ print('episode,reward=',i_episode,reward_total) # end for loop --------------------------------------

 (6)学習結果を見てみる

次を実行させると動画を見ることができます。

import matplotlib.pyplot as plt
import matplotlib.animation
import numpy as np
from IPython.display import HTML

fig2 = plt.figure(figsize=(6, 6))
plt.axis('off')

images = []
for f in frames:
  image = plt.imshow(f)
  images.append([image])
ani = matplotlib.animation.ArtistAnimation(fig2, images, interval=30, repeat_delay=1)

HTML(ani.to_jshtml()) 

f:id:akifukka:20191222113715j:plain


(7)ネットワークを可視化してみる

 デバッグのためQネット(critic)、policyネット(actor)を可視化してみました。少し謎のところもありますが・・・。

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

X, Y = np.mgrid[-1:1.1:0.1, -1:1.1:0.1]

#Q map
Z = np.zeros((21,21,21),dtype = float)
Za = np.zeros((21,21),dtype = float)
position = -1.0
for i in range(21):#position loop
  speed = -1.0
  for j in range(21):#speed loop
    a = -1.0
    a_qmax = -1.0
    qmax = -10000.0
    for k in range(21):#a loop
      q = q_net( torch.tensor([[position,speed]]).float().to(device), torch.tensor([[a]]).float().to(device)).to("cpu").detach().numpy()[0][0]
      Z[i][j][k] = q
      if (qmax < q) :
        qmax = q 
        a_qmax = a
      a = a + 0.1
    Za[j][i] = a_qmax
    speed = speed + 0.1
  position = position + 0.1

#Qmap
fig = plt.figure(figsize=(18, 6), facecolor="w")
ax = fig.add_subplot(121, projection="3d")
surf = ax.plot_surface(X, Y, Z[10])
ax.set_xlabel("speed",fontsize = 14)
ax.set_ylabel("a",fontsize = 14)
ax.set_zlabel("Q",fontsize = 14)
ax.set_title("position 0.0")

ax = fig.add_subplot(122, projection="3d")
surf = ax.plot_surface(X, Y, Za)
ax.set_ylabel("position",fontsize = 14)
ax.set_xlabel("speed",fontsize = 14)
ax.set_zlabel("a",fontsize = 14)
ax.set_title("calc from Qmap")

plt.show()

#policy map
Zpolicy = np.zeros((21,21),dtype = float)
speed = -1.0
for i in range(21):
  position = -1.0
  for j in range(21):
    Zpolicy[i][j] = policy_net(torch.tensor([[position,speed]]).float().to(device)).to("cpu").detach().numpy()[0][0]
    position = position + 0.1
  speed = speed + 0.1

fig = plt.figure(figsize=(9, 6), facecolor="w")
ax = fig.add_subplot(111, projection="3d")
ax.set_xlabel("speed",fontsize = 14)
ax.set_ylabel("position",fontsize = 14)
ax.set_zlabel("a",fontsize = 14)
surf = ax.plot_surface(X, Y, Zpolicy,color="green")
ax.set_title("policy")
plt.show()

 

f:id:akifukka:20191222114605j:plain

 左上の図はposition=0.0(図の真ん中なので谷より若干右より?)における価値関数Qネットの計算値です。

 Q = Q_net(observation, action)

 observation=[position,speed]

 a=action (+で右に押す、-で左に押す)

 aが変わっても以外にQの値は大きくは変わりませんでした。

 右上の図は行動aを変化させながらQを計算して、最もQが大きくなるaを描画したものです。

 左下の図は行動aを決めるpolicyネットの計算値です。この図から次のことがわかります。

・速度が-の時(台車が左に動いている)は-に押す(左に押す)

・速度が+の時(台車が右に動いている)は+に押す(右に押す)

 右上のQネットから算出した行動予測の図とpolicyネットを計算した左下の図はほぼ一致すると思っていたのですが、どうも一致しないようです。学習の初期ではほぼ一致してますが、学習が進むにつれ違いが出てきます。

 今のところ、この現象は理解できていません・・・。が、学習自体は上手くいっているようなので・・。

ーーーーーーーーーーーーーーーーーーーーーーーー

つづく!

 次回はDDPGの概要を理解できる様、ざっくり解説をしてみます。BipedalWalkerはその後に挑戦ということで(なんせ先に中身を理解しておかないと、太刀打ちできなさそうなので)。

 強化学習 カテゴリーの記事一覧 - やってみた!

Open AI Gym Box2D BipedalWalkerをColaboratoryで動かしてみる(3)

 前回の続きです。歩行のスケジュールをチューニングした結果、2,3歩歩けるようになりました。これ以上は胴体の傾きをつかって制御するなどもう一工夫必要そうです。 

  スケジュール制御での歩行は一旦このくらいにしておきます。とりあえず現状の歩行を以下に報告します。

(2)歩行スケジュールのチューニング

 チューニング結果の歩行スケジュールを次に示します。

xsch = np.array([0,1,2,3,4,5,6,7,8,9,10],dtype = 'float')
yhip1  = np.array([ 0.021, -0.39, -0.58, -0.54, 0.45,    0.18, 0.88,  0.78, 0.60,  0.20, 0.021],dtype = 'float')
yknee1 = np.array([  0.17,  0.85,  0.95,  0.92, -0.38, -0.05,  -0.34,  -0.090,  0.074, 0.043,  0.17],dtype = 'float')
yhip2  = np.array([ 0.18,  0.88,  0.78,  0.60,  0.20, 0.021, -0.39, -0.58, -0.54, 0.45,    0.18],dtype = 'float')
yknee2 = np.array([ -0.05,-0.34,-0.090,  0.074,  0.043 ,  0.17,  0.85, 0.95,  0.92, -0.38, 0.05],dtype = 'float')

 前回のソースリストで上記スケジュールを差替え、次の赤文字のspeedを5→4.2に変更します。

for i in range(1):
    ・
    ・
    speed = 5

 歩き方はこんな感じ。スキップのような歩き方になりました。実行するたび若干変わりますが、たぶん地面のでこぼこが毎回ちょっとずつ違うんじゃないかと思います。

f:id:akifukka:20191212220343j:plain

 さすがに完全手動でチューニングするのは面倒なので、乱数でスケジュールをちょっとずつ変更してはBipedalWalkerの報酬(reward)が多ければ採用、小さければ不採用として、ある程度自動でチューニングさせました。ただし、報酬が安定しない(地面のでこぼこの違い)ので、最終的にどの数値を採用するかはセンス?適当に選びました・・・。

 つぎはぎで、あまり参考になりませんが乱数でチューニングするのに使ったソースリストを残しておきます。

import torch
import numpy as np
import gym
import random
from pyvirtualdisplay import Display
from scipy.interpolate import interp1d

display = Display(visible=0, size=(1024, 768))
display.start()
import os
#pyvirtualdisplayの仕様変更で次の文はエラーになるのでコメントアウトします(2021/6/5修正)
#os.environ["DISPLAY"] = ":" + str(display.display) + "." + str(display.screen)

#Bipedalwalkerのv2は無くなったのでv3に変更(2021/6/5)
#env = gym.make('BipedalWalker-v3')

#角度制御ゲイン
 #P(比例)項
pgain = np.array([4,4,4,4], dtype = 'float')
 #D(微分)項
dgain = np.array([0.2,0.2,0.2,0.2], dtype = 'float')

#脚の指令角スケジュール
xsch = np.array([0,1,2,3,4,5,6,7,8,9,10],dtype = 'float')
yhip1  = np.array([ 0.021, -0.39, -0.58, -0.54, 0.45,    0.18, 0.88,  0.78, 0.60,  0.20, 0.021],dtype = 'float')
yknee1 = np.array([  0.17,  0.85,  0.95,  0.92, -0.38, -0.05,  -0.34,  -0.090,  0.074, 0.043,  0.17],dtype = 'float')
yhip2  = np.array([ 0.18,  0.88,  0.78,  0.60,  0.20, 0.021, -0.39, -0.58, -0.54, 0.45,    0.18],dtype = 'float')
yknee2 = np.array([ -0.05,-0.34,-0.090,  0.074,  0.043 ,  0.17,  0.85, 0.95,  0.92, -0.38, 0.05],dtype = 'float')
numx = xsch.shape[0]

hip1demmand = interp1d(xsch, yhip1)
knee1demmand = interp1d(xsch, yknee1)
hip2demmand = interp1d(xsch, yhip2)
knee2demmand = interp1d(xsch, yknee2)

graphschys =[]
graphschys.append([hip1demmand(0)])
graphschys.append([knee1demmand(0)])
graphschys.append([hip2demmand(0)])
graphschys.append([knee2demmand(0)])
graphschxs =[0.0]

for graphschx in np.arange(0.1,10.0,0.1):
  graphschxs.append(graphschx)
  graphschys[0].append(hip1demmand(graphschx))
  graphschys[1].append(knee1demmand(graphschx))
  graphschys[2].append(hip2demmand(graphschx))
  graphschys[3].append(knee2demmand(graphschx))

#モデルタイムステップ(1秒間に50フレーム = 0.02s)
FPS = 50.0

speed = 4.2
bestspeed = 4.2
bestreward = -9999.0

prevyhip1 = yhip1.copy()
prevyknee1 = yknee1.copy()
prevyhip2 = yhip2.copy()
prevyknee2 = yknee2.copy()

for i in range(30):
    obs = env.reset()
    done = False
    frames = []
    observations = []
    actions = []
    ts = []
    legangle = np.array([0,0,0,0], dtype = 'float')
    leganglespeed = np.array([0,0,0,0], dtype = 'float')
    legangledemand = np.array([0,0,0,0], dtype = 'float')

    x = 0.0
    t = 0
    leg1gcount = 0
    leg2gcount = 0
    leg1touchground =0
    leg2touchground =0
    reward = 0.0

    while not done and t < 1000:
        #動画用画面保存-------------------------------
        frames.append(env.render(mode = 'rgb_array'))

        #指令値生成-----------------------------------
        x += speed/FPS
        if x>10.0 :x-=10

        #接地でスケジュール強制的に進める
#        if x>0.5 and x<2.5 and leg2touchground == 1:x = 4.5 - x
#        if x>5.5 and x<7.5 and leg1touchground == 1:x = 9.5 - (x - 5)

        legangledemand[0] = hip1demmand(x)
        legangledemand[1] = knee1demmand(x)
        legangledemand[2] = hip2demmand(x)
        legangledemand[3] = knee2demmand(x)

        #角度指令値上書き。スケジュールで動かす時はコメントアウトすること
#        legangledemand = np.array([1,-1,-1,-1], dtype = 'float')

        #脚の角度 PD制御
        action = pgain * (legangledemand - legangle) - dgain * leganglespeed

        #2足歩行物理モデル1ステップ(0.02s)-------------
        observation, r, done, info = env.step(action)
        reward = reward + r

        #脚の角度,角速度,胴体角取り出し----------------
        legangle[0] = observation[4] #HIP1
        leganglespeed[0] = observation[5]
        legangle[1] = observation[6] #KNEE1
        leganglespeed[1] = observation[7]
        legangle[2] = observation[9] #HIP2
        leganglespeed[2] = observation[10]
        legangle[3] = observation[11] #KNEE2
        leganglespeed[3] = observation[11]

        #leg1 接地判定 チャタリング防止
        if observation[8] == 1:
          if leg1gcount <3 :leg1gcount += 1
        else:
          leg1gcount = 0
        
        if leg1gcount == 3:
          leg1touchground = 1
        else:
          leg1touchground = 0

        #leg2 接地判定 チャタリング防止
        if observation[13] == 1:
          if leg2gcount <3 :leg2gcount += 1
        else:
          leg2gcount = 0;

        if leg2gcount == 3:
          leg2touchground = 1
        else:
          leg2touchground = 0

        #グラフ用データ保存----------------------------
        if t == 0 :
          for i in range(24):
            observations.append([observation[i]])
          for i in range(4):
            actions.append([action[i]])
        else :
          for i in range(24):
            observations[i].append(observation[i])
          for i in range(4):
            actions[i].append(action[i])
        ts.append(t/FPS)

        #時間を進める----------------------------------
        t += 1
  
    if reward > bestreward:
#      bestspeed = speed
      bestreward = reward.copy()
      print(yhip1)
      print(yknee1)
      print(yhip2)
      print(yknee2)
      print(i,bestreward)
    else:
      #元に戻す
      yhip1 = prevyhip1.copy()
      yknee1 = prevyknee1.copy()
      yhip2 = prevyhip2.copy()
      yknee2 = prevyknee2.copy()

    #一時退避
    prevyhip1 = yhip1.copy()
    prevyknee1 = yknee1.copy()
    prevyhip2 = yhip2.copy()
    prevyknee2 = yknee2.copy()

    prev1 = random.randint(0,1)
    prev21 = random.randint(0,9)
    #両足同じスケジュールにする
    prev22 = prev21 + 5
    if prev22 >9:
      prev22 = prev22 - 10 

    if prev1 == 0:
      yhip1[prev21] += random.uniform(-0.1,0.1)
      yhip2[prev22] = yhip1[prev21]
      #端点は同じ値を維持
      yhip1[10] = yhip1[0]
      yhip2[10] = yhip2[0]
    elif prev1 == 1:
      yknee1[prev21] += random.uniform(-0.1,0.1)
      yknee2[prev22] = yknee1[prev21]
      #端点は同じ値を維持
      yknee1[10] = yknee1[0]
      yknee2[10] = yknee2[0]

#    speed = bestspeed + random.uniform(-1, 1)
    hip1demmand = interp1d(xsch, yhip1)
    knee1demmand = interp1d(xsch, yknee1)
    hip2demmand = interp1d(xsch, yhip2)
    knee2demmand = interp1d(xsch, yknee2)

#end for loop -----------------------------------------

#元に戻す
yhip1 = prevyhip1
knee1 = prevyknee1
yhip2 = prevyhip2
knee2 = prevyknee2

 ちなみにspeedも乱数で変更して報酬が多くなるよう選びました。

 pythonは慣れないこともあり変数への '=' が代入では無く、新たに作られるオブジェクトへの参照だと知らなくて少し手間取りました。最初、スケジュールを退避したはずが、退避できずにおかしな動きをしていました。copy()を追加、新らたなオブジェクトを作ることで正しく動くようになりましたが、今後も変数の使い方は注意しないと。

 a+=1、a-=1といった演算も新たにオブジェクトを作ってくれないということで、pythonに慣れるまで気を付けないといけなさそうです。(a=a+1の場合は右辺で新たなオブジェクトを生成してくれるとのこと)。

 さて、スケジュール制御はここまでとして、次回からニューラルネットワークを使っていく予定です。今日は短いけどおしまい。

つづく

 強化学習 カテゴリーの記事一覧 - やってみた!