はじめに

SimPyというPythonの離散事象シミュレーション用のパッケージを見つけて試してみたら気に入ったので自分用の備忘録も兼ねて使い方をまとめていく．

SimPyのドキュメント
SimPyのソース
第1回の記事（はじめの一歩）
第2回の記事（リソースを理解しよう：Container編）
第3回の記事（リソースを理解しよう：StoreとResource編）
第4回の記事（アニメーションを追加しよう）

今回は最終回として，これまでの回で何度かとりあげた在庫管理の例を題材にして，簡単なテーブル型のQ学習で行動政策を獲得していく在庫管理者エージェントを導入してみよう．

サンプルコードと概説

さっそくサンプルコードをみていこう．第2回や第4回にとりあげたのとほぼ同じ在庫管理モデルのシミュレーションである．Containerリソースを継承したCustomContainerクラスで対象システムをモデル化している．

今回新たに追加したのはコストの側面である．在庫管理者が発注を出すたびに発注コスト（ORDER_COST）がかかり，品物を保管しておくために1個1期あたり所定の在庫保管コスト(STOCK_COST)がかかる．さらに，顧客を待たせた場合は1個1期あたり所定のペナルティコスト（PENALTY）が必要になる．

在庫管理者はこれらのコスト（の割引現在価値）がなるべく小さくなるように，発注のタイミングと量を決められるように行動政策を学習したい．

配送業者のプロセス（deliverer()）と顧客のプロセス（customer()）は第2回や第4回のものとほぼ同じであり，recorder()プロセスは描画のためのデータを記録するためのもの，main()関数の後半は描画のためのコードなので，それらの詳細は省略する．

以下では，在庫管理者エージェント（Agentクラス）に主に注目する．

import numpy as np
import matplotlib.pyplot as plt
import simpy
import math

STATES = 25  # state in ((,-1], [0, 1], [2, 3], ..., [46,))
ACTIONS = 10  # action in (0, 3, 6, 9, 12, 15, 18, 21, 24, 25)
ORDER_COST = 1
STOCK_COST = 0.025
PENALTY = 0.1
DISCOUNT = 0.9
LEARNING = 0.2

class CustomContainer(simpy.Container):
    def __init__(self, env, capacity=float('inf'), init=0):
        self.env = env
        self.orders = []  # list of back orders
        super(CustomContainer, self).__init__(env, capacity, init)

    @property
    def shortage(self):  # total amount requested by waiting customers
        num = 0
        for customer in self.get_queue:
            num += customer.amount
        return num

    def get_state(self):  # return encoded state number
        total = self.level +sum(self.orders) -self.shortage
        return max(0, min(math.floor(total /2) +1, STATES -1))

    def get_reward(self, action):  # negative reward or cost
        reward = self.level *STOCK_COST +self.shortage *PENALTY
        if action > 0:
            reward += ORDER_COST
        return reward

class Agent:
    def __init__(self, env):
        self.env = env
        self.epsilon = 0.1
        self.recent_rewards = [0] *100
        self.Q = np.random.rand(STATES *ACTIONS).reshape(STATES, ACTIONS) *10
        self.Q[0][0] = math.inf
        for a in range(1, ACTIONS):
            self.Q[STATES -1][a] = math.inf

    @property
    def average_reward(self):
        return sum(self.recent_rewards) /len(self.recent_rewards)

    def e_greedy(self, state):
        if state == 0:  # you should order at least some amount
            if self.epsilon <= np.random.rand():
                return self.Q[state].argmin()  # greedy
            else:
                return np.random.randint(1, ACTIONS)  # random
        elif state == STATES -1:  # you cannot order
            return 0
        elif self.epsilon <= np.random.rand():
            return self.Q[state].argmin()  # greedy
        else:
            return np.random.randint(0, ACTIONS)  # random

    def move(self):
        state_to = self.env.model.get_state()
        while True:
            state_from = state_to
            action = self.e_greedy(state_from)
            if action > 0:  # when you order
                self.env.model.orders.append(action *3)
                self.env.process(deliverer(self.env))  # activate deliverer
                if not self.env.record.triggered:
                    self.env.record.succeed()  # record log
            yield self.env.timeout(1)  # periodic inventory check
            state_to = self.env.model.get_state()
            reward = self.env.model.get_reward(action)
            self.recent_rewards.append(reward)
            self.recent_rewards.pop(0)
            self.update_Q(state_from, action, state_to, reward)

    def update_Q(self, state_from, action, state_to, reward):
        Q_to_min = self.Q[state_to].min()
        self.Q[state_from][action] += LEARNING *(reward +DISCOUNT *Q_to_min -self.Q[state_from][action])

    def update_epsilon(self):
        self.epsilon *= 0.9

def deliverer(env):
    yield env.timeout(3)  # delivery lead time = 3
    env.model.put(env.model.orders.pop(0))

def customer(env):
    while True:
        time_to = np.random.exponential(1)
        yield env.timeout(time_to)
        how_many = np.random.randint(1, 8)  # mean = 4
        env.model.get(how_many)
        if not env.record.triggered:
            env.record.succeed()  # record log

def recorder(env):  # process for recording log for visualization
    _t = []
    env.record = env.event()
    while True:
        yield env.record
        _t.append(env.now)
        env.y11.append(env.model.level)
        env.y12.append(sum(env.model.orders))
        env.y13.append(env.model.shortage)
        env.t.append(env.now)
        env.z.append(env.agent.average_reward)
        if env.now > 200:
            t_min = env.now -200
            env.y11 = [
                env.y11[i] for i in range(len(_t)) if _t[i] > t_min
                ]
            env.y12 = [
                env.y12[i] for i in range(len(_t)) if _t[i] > t_min
                ]
            env.y13 = [
                env.y13[i] for i in range(len(_t)) if _t[i] > t_min
                ]
            _t = [_t[i] for i in range(len(_t)) if _t[i] > t_min]
            env.x = [_t[i] -max(_t) +200 for i in range(len(_t))]
        else:
            env.x = _t
        env.record = env.event()

def main():
    env = simpy.Environment()
    env.model = CustomContainer(env)
    env.agent = Agent(env)
    env.process(recorder(env))
    env.process(customer(env))
    env.process(env.agent.move())
# ---------- code for visualization ----------
    env.x = []
    env.y11 = []
    env.y12 = []
    env.y13 = []
    env.t = []
    env.z = []

    fig = plt.figure(1, figsize=(12, 8))

    ax1 = fig.add_subplot(221)
    ax1.set_xlabel('time')
    ax1.set_ylabel('cost')
    ax1.set_xlim(0, 50000)
    ax1.set_ylim(0, 2)
    line1, = ax1.plot(env.t, env.z, label='average cost')
    ax1.legend()
    ax1.grid()

    ax2 = fig.add_subplot(222)
    ax2.set_xlabel('time')
    ax2.set_ylabel('number')
    ax2.set_xlim(0, 200)
    ax2.set_ylim(0, 60)
    line21, = ax2.plot(env.x, env.y11, label='at hand')
    line22, = ax2.plot(env.x, env.y12, label='ordered')
    line23, = ax2.plot(env.x, env.y13, label='shortage')
    ax2.legend()
    ax2.grid()

    ax3 = fig.add_subplot(223)

    for t in range(1, 1000):
        env.run(until=t*50)  # stepwise execution
        if t % 50 == 0:
            env.agent.update_epsilon()
            print('epsilon = {}'.format(env.agent.epsilon))
        line1.set_data(env.t, env.z)
        line21.set_data(env.x, env.y11)
        line22.set_data(env.x, env.y12)
        line23.set_data(env.x, env.y13)
        heatmap = ax3.imshow(env.agent.Q, vmin=2, vmax=10, cmap='jet', aspect=0.25)
        bar = plt.colorbar(heatmap, ax=ax3)
        plt.pause(0.1)
        bar.remove()
    plt.show()
# ---------- ---------- ---------- ----------

if __name__ == "__main__":
    main()

Agentクラスは，在庫管理者の行動政策を標準的なテーブル型のモデルフリーQ学習で獲得していくエージェントになっている．状態は，現在の在庫量（level）とバックオーダ（oeders）の総量の和から，待たされている顧客の要求量の総和（shortage）を減算したスカラー値を離散化したもの，行動は，簡単のため，3の倍数に限定した発注量である．

Q値を乱数で初期化し，ε-greedyで行動しながら，標準的な方法でQ値を更新していっていることがわかる．この記事は強化学習（Q学習）の解説を意図したものではないので，学習の詳細にはこれ以上踏み込まないことにする．

ここで，注目してほしいことは，このε-greedyに従った在庫管理者エージェントの行動政策が，SimPyのプロセスとしてコード化されており，それが，main()の中で

env.process(env.agent.move())

とするだけで，簡単にシミュレーションに統合できていることである．

まとめ

ここまで5回に渡ってPythonの離散事象シミュレーション用のパッケージSiｍPyについてまとめてきた．もし最終回まで目を通してくださった方がいたとしたら有り難いと思う．何かの参考になれば幸いです．

この記事のプログラムのソース（GitHub）