マルチエージェント深層強化学習を勉強しよう　第一回-MADDPGの解説及び実装

#この記事を書く理由及び今後の目標について
マルチエージェントシステムについて研究している航空宇宙工学科の学生ですが。最近はマルチエージェント強化学習(MARL)について興味がわいたので、そのアルゴリズム手法を何回かに分けてまとめようと思いました。入門者の参考及び自分用のメモでもあるので、できるだけわかりやすく伝えていきたいと思います。
#概要
本日はActor-Critic手法として有名なDDPG(Deep Deterministic Policy Gradient)を拡張した手法である[MADDPG(Multi-Agent Deep Deterministic Policy Gradient)](https://arxiv.org/pdf/1706.02275.pdf)について紹介したいと思います。この記事はアルゴリズムの簡単な解説及びPytorchを用いる実装を示すが、具体的な理論については省略させていただきます。Actor-CriticやDDPGについてわからない人は以下の関連記事から読むのをお勧めします。
####関連記事及び参考Github
1.[[論文解説] Deterministic Policy Gradient Algorithms](https://qiita.com/ku2482/items/7fd5fb65666328209ace)
2.[Sarsa、Actor-Critic法の解説および倒立振子問題を解く](https://qiita.com/KokiSakano/items/0ab412ae9f69865f26d2)
3.[openaiにあるmaddpgの元コード](https://github.com/openai/maddpg/)(Tensowflowで実装されていました）
4.[訓練環境:Multi-Agent Particle Environments (MPE)](https://github.com/openai/multiagent-particle-envs)
#なぜマルチエージェント
強化学習は単一エージェントに対して著しく発展してきたが、現実世界では、エージェント常に複数存在します。より効率的に仕事をさせるためにはエージェントたちを協調的(たまには敵対)に制御する必要があります。(例えば工場内のUAVとUGVの協調作業など)MARLは近年では盛んに研究されています、単一エージェントと比べて、その行動や学習はとても難しいです。論文の中でも言及していますが、エージェントの数が増えれば増えるほど、学習する際の勾配は正しい方向に進む確率は指数関数的に下がります。

#アルゴリズム
##全体像
![fig_maddpg_paper.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/507786/11c5141f-ff15-66dc-b061-583807a15708.png)
論文にある図をそのまま引用させていただきますが。以上の図はアルゴリズムの全体像を示しています。この図は何の意味をしているかと言いますと、まずN個のエージェントはそれぞれ方策$\boldsymbol{\pi}= \\{ \boldsymbol{\pi_{1}},...,\boldsymbol{\pi_{N}}\\}$とします。当たり前ですが状態値$\boldsymbol{o_{i}}$をそれぞれの方策に渡すとそれぞれのエージェントがとるべき行動$\boldsymbol{a_{i}}$が吐き出されます。これはそれぞれのエージェント(actor)は**Decentralized**(分散)な方策をとることを意味しています。次にすべてのエージェントを管理できるN個の価値関数(Q関数）があります。この価値関数はすべてのエージェントから状態値$\boldsymbol{x=\\{o_{1},...o_{N}\\}}$と行動$\boldsymbol{a=\\{a_{1},...a_{N}\\}}$を受け取ってその価値を返します。言い換えると**Centralized**なcritic(批評家)であることを意味しています。つまり、N人のcriticがすべてのactor(演者)について評価し、actorたちはその評価に応じてそれぞれのとる方策を更新します。
##少しだけ数式を見よう
以上の議論から、まずactorたち最大化するべき目的関数はcriticからの評価であることは明らかです。それぞれの確率的な方策を構成するパラメータを$\boldsymbol{\theta_{i}}$とすると、
$$
J(\boldsymbol{\theta_{i}})=\mathbb{E_{s \sim p^{\mu},a_{i} \sim \pi_{i}}}[Q_{i}(\boldsymbol{x},a_{1},...a_{N})]
$$
目的関数は以上と書きます。両辺をパラメータである$\boldsymbol{\theta_{i}}$で微分すると、確率的な方策を更新する勾配は以下になります。
$$
\nabla_{\boldsymbol{\theta_{i}}} J(\boldsymbol{\theta_{i}})=\mathbb{E_{s \sim p^{\mu},a_{i} \sim \pi_{i} }}[\nabla_{\boldsymbol{\theta_{i}}}\log{\boldsymbol{\pi_{i}}(a_{i}|o_{i})}Q_{i}(\boldsymbol{x},a_{1},...a_{N})]
$$
しかしすでに知られているように、**MADDPG**は名前の通り**決定的方策をとるoff-Policy**の手法であるため、以上の式は以下のように書き換えられます。(off-Policyというのは簡単に言うと取るべき行動は現在の方策(Policy)に関係なく(Off)常に最善を尽くす手法のことです。)
方策の勾配式を決定的方策に書き換えると
$$
\nabla_{\boldsymbol{\theta_{i}}} J(\boldsymbol{\mu_{i}})=\mathbb{E_{x,a_{i} \sim D }}[\nabla_{\boldsymbol{\theta_{i}}}{\boldsymbol{\mu_{i}}(a_{i}|o_{i})}\nabla_{\boldsymbol{a_{i}}}Q_{i}^{\mu}(\boldsymbol{x},a_{1},...,a_{i}^{'},...,a_{N})|_{a _{i} ^{'}=\mu _{i}  (o)}]
$$
決定的な方策っていうのは状態値$o$をある関数$\mu(\cdot)$に渡すと決められた行動が出力されます、連続値空間でよく用いられます。それに対して確率的な方策っていうのは状態値$o$をある関数$\pi(\cdot)$に渡すとそれぞれの行動を行う離散的な確率が出力されます。またこの式に注意するべきところはエージェント$i$を評価する価値関数$Q _{i}(\cdot)$に渡す行動は$i$の行動$a _{i} ^{'}$のみ一ステップ分更新します。

ではcriticたちの評価どうやって正しく評価しているかどうかを見るでしょう？ここは慣例の**TD誤差**(説明は割愛します)の出番です。価値関数を更新するロス関数は以下に示します。
$$
\mathcal{L}(\theta _{i}) = \mathbb{E} _{x,a,r,x^{'}}[(Q _{i}^{\mu}(\boldsymbol{x},a _{1},...a _{N})-y^2)]
$$
ここで,$x,a,r,x^{'}$はあらかじめ経験を貯蓄していたバッファー$\mathcal{D}$からサンプルした値です。また、学習のターゲットである$y$は以下の式で求めます。
$$
y = r _{i} + \gamma Q _{i}^{\boldsymbol{\mu^{'}}}(\boldsymbol{x^{'}},a _{1}^{'},...,a _{i},...,a _{N}^{'})| _{a'=\mu'(o)}
$$
ここで$\boldsymbol{\mu^{'}}(\cdot)$と$\boldsymbol{Q^{'}}(\cdot)$ はターゲット関数になっていて、価値関数$\boldsymbol{ Q }(\cdot)$はこいつらに準じてsoft-updateしていく感じです。

##updateの手順
1.経験を貯めたReplay-Buffer$\mathcal{D}$からbatch-size分の$(x',a,r,x)$をサンプルします
2.criticのパラメータ$\theta_{c}$をロス関数の勾配で更新します
$$
\theta_{c} \leftarrow \theta_{c} + \alpha  \nabla\mathcal{L}(\theta _{i})
$$
3.actorのパラメーター$\theta _{a}$を目標関数の勾配で更新します

$$
\theta_{a} \leftarrow \theta _{a} + \alpha \nabla\mathcal{J}(\theta _{i})
$$

4.ターゲットネットワークのパラメータをsoft-updateします。(soft-updateというのは元のネットワークパラメータに準じてじわじわ更新していく感じです)

```math
\theta_{c} ^{'} \leftarrow \tau\theta_{c}  + (1-\tau)\theta_{c}^{'} \\
\theta_{a} ^{'} \leftarrow \tau\theta_{a}  + (1-\tau)\theta_{a}^{'} \\
```

#実装
数式見てもパッと来ないですよね！いざ実装してみましょう！(Pytorch初心者なので間違いがあったら指摘してください）

##ライブラリー
ここでmake_envというのは[訓練環境:Multi-Agent Particle Environments (MPE)](https://github.com/openai/multiagent-particle-envs)からインポートしたものです、各自ダウンロードしといてください！ほかはPytorchでよく使うものをインポートしています。

```python:lib.py
from make_env import make_env
import numpy as np
import copy
from collections import deque
import gym
import random
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torch.nn.utils import clip_grad_norm_
import matplotlib
import matplotlib.animation as animation
import matplotlib.pyplot as plt

#GPU環境の確認
device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')
print(device) #cuda:0
```
##Replay-Buffer
学習経験を貯めるBufferの実装、特に変わった部分なく、好みでdequeをつかっています。

```python:buffer.py
class ReplayBuffer:
    def __init__(self,memory_size=int(1e+6)):
        self.memory = deque([],maxlen=memory_size)
        self.is_gpu = torch.cuda.is_available
        return
    def cache(self,state,next_state,action,reward,done):
        if self.is_gpu:
            state = torch.tensor(state,dtype=torch.float).cuda()
            next_state = torch.tensor(next_state,dtype=torch.float).cuda()
            action = torch.tensor(action,dtype=torch.float).cuda()
            reward = torch.tensor(reward).cuda()
            done = torch.tensor([done]).cuda()
        else:
            state = torch.tensor(state,dtype=torch.float)
            next_state = torch.tensor(next_state,dtype=torch.float)
            action = torch.tensor(action,dtype=torch.float)
            reward = torch.tensor(reward)
            done = torch.tensor([done])
        self.memory.append((state,next_state,action,reward,done))
    def sample(self,batch_size=64):
        batch = random.sample(self.memory,batch_size)
        state,next_state,action,reward,done = map(torch.stack,zip(*batch))
        return state,next_state,action,reward.squeeze(),done.squeeze()
```

##モデル
モデルは隠れサイズが64の全結合層ネットワーク(MLP)によって構成されています。

```python:model.py
#Actorのネットワーク
class PolicyNetwork(nn.Module):
    def __init__(self,num_state,num_action,hidden_size=64):
        super(PolicyNetwork,self).__init__()
        self.fc1 = nn.Linear(num_state,hidden_size)
        self.fc2 = nn.Linear(hidden_size,hidden_size)
        self.fc3 = nn.Linear(hidden_size,num_action)
    
    def forward(self,x,index):
        x = x[:,index]
        h = F.relu(self.fc1(x))
        h = F.relu(self.fc2(h))
        action = self.fc3(h)
        return action

#Criticのネットワーク
class QNetwork(nn.Module):
    def __init__(self,num_state,num_action,agent_num,hidden_size=64,init_w=3e-3):
        super(QNetwork,self).__init__()
        input_size1 = num_state * agent_num
        input_size2 = hidden_size + num_action * agent_num
        self.fc1 = nn.Linear(input_size1,hidden_size)
        self.fc2 = nn.Linear(input_size2,hidden_size)
        self.fc3 = nn.Linear(hidden_size,1)
        self.fc3.weight.data.uniform_(-init_w, init_w)
        self.fc3.bias.data.uniform_(-init_w, init_w)
    
    def forward(self,states,actions):
        states = states.view(states.size()[0],-1)
        actions = actions.view(actions.size()[0],-1)
        h = F.relu(self.fc1(states))
        x = torch.cat([h,actions],1)
        h = F.relu(self.fc2(x))
        q = self.fc3(h)
        return q
```

##ノイズ
学習の時に行動にノイズをかかって学習しやすくするので、ノイズは**Ornstein–Uhlenbeck process**を採用しています。この実装は[DDPGでPendulum-v0を学習](https://qiita.com/sekishoku/items/f929d7737e804b9aec7a)から拝借させていただきました。

```python:noise.py
class OrnsteinUhlenbeckProcess:
    def __init__(self, theta=0.15, mu=0.0, sigma=0.2, dt=1e-2, x0=None, size=1, sigma_min=None, n_steps_annealing=1000):
        self.theta = theta
        self.mu = mu
        self.sigma = sigma
        self.dt = dt
        self.x0 = x0
        self.size = size
        self.num_steps = 0

        self.x_prev = self.x0 if self.x0 is not None else np.zeros(self.size)

        if sigma_min is not None:
            self.m = -float(sigma - sigma_min) / float(n_steps_annealing)
            self.c = sigma
            self.sigma_min = sigma_min
        else:
            self.m = 0
            self.c = sigma
            self.sigma_min = sigma

    def current_sigma(self):
        sigma = max(self.sigma_min, self.m * float(self.num_steps) + self.c)
        return sigma

    def sample(self):
        x = self.x_prev + self.theta * (self.mu - self.x_prev) * self.dt + self.current_sigma() * np.sqrt(self.dt) * np.random.normal(size=self.size)
        self.x_prev = x
        self.num_steps += 1
        return x
```
##エージェント
各種パラメータの設定は基本的に論文通りにしています、学習しやすくするために勾配のノルムをclipしました。

```python:maddpg.py
class MaddpgAgents:
    def __init__(self,observation_space,action_space,num_agent,gamma=0.95,lr=0.01,batch_size=1024,memory_size=int(1e6),tau=0.01,grad_norm_clipping = 0.5):
        self.num_state = observation_space 
        self.num_action = action_space
        self.n = num_agent
        self.gamma = gamma
        self.actor_group = [PolicyNetwork(self.num_state,self.num_action).to(device) for _ in range(self.n)]
        self.target_actor_group = copy.deepcopy(self.actor_group)
        self.actor_optimizer_group = [optim.Adam(self.actor_group[i].parameters(),lr=0.001
                                                 ) for i in range(self.n)]
        self.critic_group = [QNetwork(self.num_state,self.num_action,self.n).to(device) for _ in range(self.n)]
        self.target_critic_group = copy.deepcopy(self.critic_group)
        self.critic_optimizer_group = [optim.Adam(self.critic_group[i].parameters(),lr=lr) for i in range(self.n)]
        self.buffer = ReplayBuffer(memory_size=memory_size)
        self.loss_fn = torch.nn.MSELoss()
        self.batch_size = batch_size
        self.is_gpu = torch.cuda.is_available
        self.noise = OrnsteinUhlenbeckProcess(size=self.num_action)       
        self.grad_norm_clipping = grad_norm_clipping
        self.tau = tau
    
    @torch.no_grad()
    def td_targeti(self,reward,state,next_state,done,agent_index):
        next_actions = []
        for i in range(self.n):
            actionsi = torch.tanh(self.target_actor_group[i](state,i))
            actionsi = actionsi[:,np.newaxis,:]
            next_actions.append(actionsi)
        next_actions = torch.cat(next_actions,dim=1)
        next_q = self.target_critic_group[agent_index](next_state,next_actions)
        if self.n != 1:
            reward = reward[:,agent_index]
            done = done[:,agent_index]
        reward = reward[:,np.newaxis]
        done = done[:,np.newaxis]
        done = torch.tensor(done,dtype=torch.int)
        td_targeti = reward + self.gamma * next_q*(1.-done.data)
        return td_targeti.float()

    def update(self):
        for i in range(self.n):
            state,next_state,action,reward,done = self.buffer.sample(self.batch_size)
            td_targeti = self.td_targeti(reward,state,next_state,done,i)　
            current_q = self.critic_group[i](state,action)

            #criticの更新
            critic_loss = self.loss_fn(current_q,td_targeti)
            self.critic_optimizer_group[i].zero_grad()
            critic_loss.backward()
            clip_grad_norm_(self.critic_group[i].parameters(),max_norm=self.grad_norm_clipping)
            self.critic_optimizer_group[i].step()

            #actorの更新
            ac = action.clone()
            ac_up = self.actor_group[i](state,i)
            ac[:,i,:] = torch.tanh(ac_up)
            pr = -self.critic_group[i](state,ac).mean()
            pg = (ac[:,i,:].pow(2)).mean()
            actor_loss = pr + pg*1e-3
            self.actor_optimizer_group[i].zero_grad()
            clip_grad_norm_(self.actor_group[i].parameters(),max_norm=self.grad_norm_clipping)
            actor_loss.backward()
            self.actor_optimizer_group[i].step()

        #soft-update
        for i in range(self.n):
            for target_param, local_param in zip(self.target_actor_group[i].parameters(), self.actor_group[i].parameters()):
                target_param.data.copy_(self.tau * local_param.data + (1.0 - self.tau) * target_param.data)
            for target_param, local_param in zip(self.target_critic_group[i].parameters(), self.critic_group[i].parameters()):
                target_param.data.copy_(self.tau * local_param.data + (1.0 - self.tau) * target_param.data)
 
    #行動を取得する、学習の時はノイズかかるが実際はgreedy=Trueで貪欲探索します       
    def get_action(self,state,greedy=False):
        state = torch.tensor(state,dtype=torch.float).cuda()
        state = state[np.newaxis,:,:]
        actions = []
        for i in range(self.n):
            action = torch.tanh(self.actor_group[i](state,index=i))
            if not greedy:
                action += torch.tensor(self.noise.sample(),dtype=torch.float).cuda()
            actions.append(action)
        actions = torch.cat(actions,dim=0)
        return np.clip(actions.detach().cpu().numpy(),-1.0,1.0)

```
##学習と検証
ここは特に変わった部分がなく、強化学習によくある流れですね。

```python:train_and_test.py
#各種設定
num_episode = 20000  #学習エピソード数（論文では25000になっています）
memory_size = 100000  #replay bufferの大きさ
initial_memory_size = 100000  #最初貯める数
#ログ用の設定
episode_rewards = []
num_average_epidodes = 100

#今回採用した環境はsimple_spreadですが、ほかにもいろいろあります
env = make_env('simple_spread')
max_steps = 25  # エピソードの最大ステップ数
agent = MaddpgAgents(18, 5, num_agent=env.n,memory_size=memory_size)

#最初にreplay bufferにノイズのかかった行動をしたときのデータを入れる
state = env.reset()
for step in range(initial_memory_size):
    if step % max_steps == 0:
        state = env.reset()
    actions = agent.get_action(state)
    next_state, reward, done, _ = env.step(actions)
    agent.buffer.cache(state,next_state,actions,reward,done)
    state = next_state
print('%d Data collected' % (initial_memory_size))

for episode in range(num_episode):
    state = env.reset() 
    episode_reward = 0
    for t in range(max_steps):
        actions = agent.get_action(state)
        next_state, reward, done, _ = env.step(actions)
        episode_reward += sum(reward)
        agent.buffer.cache(state,next_state,actions,reward,done)
        state = next_state
        if all(done):
            break
    if episode % 4 == 0:
        agent.update()
    episode_rewards.append(episode_reward)
    if episode % 20 == 0:
        print("Episode %d finished | Episode reward %f" % (episode, episode_reward))

#累積報酬の移動平均を表示
moving_average = np.convolve(episode_rewards, np.ones(num_average_epidodes)/num_average_epidodes, mode='valid')
plt.plot(np.arange(len(moving_average)),moving_average)
plt.title('MADDPG: average rewards in %d episodes' % num_average_epidodes)
plt.xlabel('episode')
plt.ylabel('rewards')
plt.show()

env.close()

#結果を環境で試す１
state = env.reset()
episode_reward = 0
frames = []
env.render()
screen = env.render(mode='rgb_array')
frames.append(screen[0])
for t in range(max_steps):
    actions = agent.get_action(state,greedy=True)
    next_state, reward, done, _ = env.step(actions)
    episode_reward += sum(reward)
    agent.buffer.cache(state,next_state,actions,reward,done)
    env.render()
    screen = env.render(mode='rgb_array')
    frames.append(screen[0])
    state = next_state
print(episode_reward)
```
##学習結果
![maddpg_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/507786/9f737ae4-96bf-87dd-c424-fe683c65057a.png)
以上のコードを実行し、20000episodeで学習した結果はこの図になります。報酬は$-500$ぐらいで収束していますがまだまだハイパーパラメータ調整することで上昇する余地があると思います。

##全体のコード
https://colab.research.google.com/drive/1VHONrpvf_Xq79T8CrOb27IwAysmvhKiO?usp=sharing
Googl colabに全体のコードをつくりました、そのうちgithubでMARLプロジェクトを立ち上げようと思いますので興味ある人実行してみてください。(動画出力うまくいきませんでしたので、動画見ようと思うひとはローカルで実行したほうがいいですね！理由わかるひといたら教えてください...)

学習済みのイメージ
https://www.youtube.com/watch?v=8hBbDHOwVkY

#まとめ
今回はMADDPGについて紹介をしました、これを改善するためにまだまだいろんな手法がありますので今後逐次紹介していきたいと思います。記事書くのに慣れてないので大目にみていただけると助かります、ご指摘なども承っております。