Pythonを用いた離散事象シミュレーション入門 #2

# はじめに

前回は，離散事象シミュレーションのメカニズムと，その最も基本的な実装方法について理解した．具体的には，イベントカレンダ（事象を生起タイミングの順に並べたリスト）を保持しておき，その先頭から順に事象を1つずつ取り出して，そのタイプに応じた処理を実行していくという流れであった．サンプルコードでは，モデルの`run()`メソッドの中にwhileループがあり，そのループを周回するごとに，イベントカレンダから先頭の事象を取り出し，そのタイプをif文で確認し，タイプに応じた処理に誘導する，というステップを実行していた．この実装方法は，離散事象シミュレーションのメカニズムをそのまま素直にコード化する，最もナイーブなコーディングスタイルであるともいえる（以降，このスタイルをイベント志向と呼ぶことがある）．

実は，このイベント志向のコーディングスタイルには（少なくとも）2つの問題がある．1つ目は，事象のタイプ数や対象システムの構成要素数などに応じて，`run()`メソッドが肥大化していくという点である．2つ目は，対象システムの状態を変化させる処理が，その変化を引き起こす主体や変化を引き受ける客体ではなく，変化を引き起こすシグナル，すなわち事象ごとにまとめられるため，（主体や客体の視点から）互いに関連する一連の変化が，細かな部分に分割され，複数の異なる箇所に記述されることになるという点である．これらはいずれも，コード全体の見通しを阻害するものであり，その影響は，対象システムの規模が大きくなるにつれてどんどん深刻になっていく．

これらの問題を解消するための方法の1つとして，プロセス志向のコーディングスタイルを採用することが挙げられる．PythonのSimPyモジュールを導入すると，簡単にこのプロセス志向のスタイルで開発を進めることができるようになる．今回は，このモジュールを使った離散事象シミュレーションのコーディングの基礎を身につけよう．

# SimPyの導入

## SimPyの概要とインストール

最初にすべきことはSimpyモジュールのインストールである．Google Colaboratotyでは，下のように，pipを使って簡単にインストールすることができる（なお，自分のローカルな環境にインストールする際には行頭の!は不要）．

```console
! pip install simpy
```

SimPyの主な構成要素は，core.pyの中にあるEnvironmentクラス，events.pyの中にあるEventクラスとそのサブクラス群，resourcesディレクトリの中にあるリソース関連のクラス群，の3つだと考えればよいだろう．また，これらに加えてジェネレータとして実装されるプロセス関数・メソッドが重要な役割を果たす．リソース関連のクラス群については次回に取り上げることにして，今回はそれ以外の3つに基づいて，プロセス志向コーディングの基本をおさえよう．

## シミュレーション環境と事象

Environmentクラスは，シミュレーション時間の管理やイベントカレンダの操作など，離散事象シミュレーションのための最も基本的な機能を提供してくれる．したがって，SimPyを用いてシミュレーションモデルを開発する際には，シミュレーション環境（=Environmentクラスのインスタンス）を必ず1つ生成することになる．一方，Eventクラスは，事象を表現するためのクラスであり，後で見るようにいくつかのサブクラスも用意されている．

ある事象（=Eventクラスのインスタンス）をイベントカレンダに登録することをtriggerするという．通常のEventクラスの事象の場合，それはtriggerされたのと同時刻に生起することになる．一方，triggerしてから所定の時間遅れ後に生起するような事象を利用したいことも多い．その場合は，Eventクラスのサブクラスである，Timeoutクラスを利用すればよい．

SimPyでは，ある事象が生起した際に実行される処理をその事象のコールバックと呼び，各事象`e`にコールバックの集合`e.callbacks`を付与することで，事象`e`の生起に伴って発生する処理がまとめて実行されるようになっている．簡単な例を見てみよう．

```python
import random
import simpy

class Skelton:
    def __init__(self, env):
        self.env = env  # pointer to the SimPy environment
        self.count = 0  # an example state variable

    def update(self, e):
        self.count += 1  # increment the event counter

    def print_state(self):
        print('{} th event occurs at {}'.format(self.count, round(self.env.now)))

    def run(self, horizon):
        while True:
            e = simpy.Timeout(self.env, random.expovariate(1))  # create an Timeout instance
            e.callbacks.append(self.update)  # register update() method in e's callbacks
            if self.env.now > horizon:  # if horizen is passed
                break  # stop simulation
            else:
                self.print_state()
                self.env.step()  # process the next event

env = simpy.Environment()
model = Skelton(env)
model.run(200)
```

この例では，前回のスケルトンモデルとほぼ同じ機能をSimPyのEnvironmentクラスとTimeoutクラスを利用して再現している．前回自作していたEventクラスやCalendarクラスは（それらに相当する機能をSimPyが提供してくれるので）不要である．末尾の3行を見てほしい．シミュレーション環境（`=env`）を生成した後，それを引数として，対象システムのモデル（`=model`）を生成している．そして，そのモデルの`run()`メソッドを，`horizon=200`で実行している．

Skeltonクラスの中身を確認すると，`run()`メソッドにはwhileループがあり，その各周回で，Timeoutクラスの事象（`=e`）を生成して，そのコールバックの集合（`=e.callbacks`）に，`update()`メソッドを登録していることがわかる．なお，`update()`メソッドは，`count`をインクリメントするだけのダミーである．また，事象のコールバックは，その事象を唯一の引数とする関数（厳密には，呼び出し可能オブジェクト）の形式でなければならない．

Timeoutクラスの事象を生成する際の第1引数は，対応するシミュレーション環境`env`，第2引数は，時間遅れの長さ（上の例では指数分布に従う乱数で与えられている）である．なお，Timeout事象は，それを生成した際に自動的にtriggerされる（が，通常のEventクラスの事象は，後述するように，明示的にtriggerする必要がある）．

シミュレーション時間は，シミュレーション環境`env`の変数`now`で管理されている（上の`run()`メソッドの中からは`self.env.now`で参照できるようになっている）．この値が，引数として渡された`horizon`以上であれば，whileループを抜け，シミュレーションを終了する．そうでなければ，シミュレーション環境の`step()`メソッドを呼んでいるが，これはイベントカレンダから先頭の事象`e`を1つ取り出して生起させる（つまり，`e.callbacks`に含まれているコールバックを順に実行していく）という処理に対応している．

## プロセス関数・メソッド

上の例のSkeltonクラスは，機能の一部をシミュレーション環境に任せているため，前回と比べるとかなりシンプルになっている．ただしそれだけでは，基本的・共通的な機能はSimPyが面倒を見てくれるので，自分でコーディングしなければならない部分が少なくなるということにしか過ぎない．実は，SimPyを導入することの本質的なメリットは，むしろその先にこそあるといえる．

この本質的なメリットをもたらすものが，プロセス関数・メソッドである．これによって，SimPyでは，プロセス志向でコーディングしていくことが可能になるのである．次に，その基本的な仕組みを例を用いて説明していこう．下の例を見てほしい．

```python
class Skelton2:
    def __init__(self, env):
        self.env = env  # pointer to the SimPy environment
        self.count = 0  # an example state variable

    def print_state(self):
        print('{} th event occurs at {}'.format(self.count, round(self.env.now)))

    def process_method(self):  # an example process method
        while True:
            self.print_state()
            yield simpy.Timeout(self.env, random.expovariate(1))
            self.count += 1  # corresponding to Skelton's update()

def process_func(env):  # an example process function
    while True:
        env.model.print_state()
        yield simpy.Timeout(env, random.expovariate(1))
        env.model.count += 1  # corresponding to Skelton's update()

env = simpy.Environment()
env.model = Skelton2(env)
# simpy.Process(env, process_func(env))  # when using process function
simpy.Process(env, env.model.process_method())  # when using process method
env.run(until=200)
```

これは，上でみた例をプロセス関数・メソッドを用いて書き直したものである．Skeltonクラスにあった`run()`メソッド（と`update()`メソッド）がなくなり，Skelton2クラスには，`process_method()`というメソッドが新たに登場していることに気がついたと思う．これがプロセスメソッドである．なお，このプロセスメソッドは利用せず，代わりに，同じ機能を果たすプロセス関数（上の例では，`process_func()`関数）を用いてもよい（この例では双方とも用意されているが，実際にはどちらか一方だけでよい）．

`process_method()`や`process_func()`の中にyield文があることからわかるように，これらはPythonのジェネレータになっている．通常の関数やメソッドがreturnで結果を返して終了するのに対して，ジェネレータはyieldで結果を返すとそこで一時停止するだけで，終了はしない．そして，後で再開命令のシグナルを受け取ると，yield文の先から処理を再開する．

このように，プロセス関数・メソッドは，Eventクラスのインスタンスをyieldする形で定義されたジェネレータであり，SimPyでは，これをプロセス志向コーディングのためのトリックとして利用している．具体的には，プロセス関数・メソッドがある事象`e`をyieldすると，`e.callbacks`に，そのプロセス関数・メソッドの再開命令が自動的に追加されるようになっているのである．

プロセス関数・メソッドは，yieldした事象が生起すると再開されるので，その再開後の部分にその事象によって引き起こされる状態変化（この例では，`count`のインクリメント）を直接記述しておけばよいことになる．したがって，この例では，`update()`メソッドをコールバックの集合に登録することは不要になっている．この例のように，1つのTimeout事象と単純な状態変化（`count`のインクリメント）だけではメリットは実感しにくいかもしれないが，複数の事象の影響を受けながら，複雑に状態変化が進んでいくようなプロセスも，これによって直感的に記述できるようになる．

なお，作成したプロセス関数・メソッドがシミュレーション内で実行されるようにするためには，それをシミュレーション環境に登録しておかなければならない．これを行っているのが，下から2行目（やコメントアウトされている3行目）である．具体的には，Processクラスのインスタンスを作成していることがわかる．この際に，該当のプロセスをスタートさせるシグナルを発する事象（Initialize事象）を生成し，triggerするという処理が，裏で自動的に実行されている．

また，一番下の行にある，シミュレーション環境の`run()`メソッドは，`step()`メソッドを繰り返すラッパーである．`run(until=時刻)`あるいは`run(until=事象)`として，ある時刻まであるいはある事象が生起するまでシミュレーションを進めることができる．この例では，シミュレーション時間が200になるまでシミュレーションを進めている．

## 複数プロセスの相互作用

複数のプロセス関数・メソッドを定義して，同じシミュレーションの中で互いに関連付けながら実行していくことができる．ここではその例をみておこう．下に簡単な例を示す．

```python
class Skelton3(Skelton):
    def __init__(self, env):
        super().__init__(env)

    def main_process(self):
        while True:
            self.print_state()
            yield self.env.timeout(random.expovariate(1))  # shortcut for simpy.Timeout()
            self.count += 1
            if self.count %3 == 0:
                self.env.signal4A.succeed()  # signal for resuming sub process A

    def sub_process_A(self):
        self.env.signal4A = self.env.event()  # create the first signal
        while True:
            yield self.env.signal4A
            print('> sub process A is resumed at {}'.format(round(self.env.now)))
            self.env.signal4A = self.env.event()  # create the next signal
            if self.count %5 == 0:
                self.env.process(self.sub_process_B())  # register sub process B

    def sub_process_B(self):
        print('>> sub process B is started at {}'.format(round(self.env.now)))
        yield self.env.timeout(10)  # shortcut for simpy.Timeout()
        print('>> sub process B is finished at {}'.format(round(self.env.now)))

env = simpy.Environment()
env.model = Skelton3(env)
env.process(env.model.main_process())  # shortcut for simpy.Process()
env.process(env.model.sub_process_A())  # shortcut for simpy.Process()
env.run(until=200)
```

Skelton3クラスの中に，`main_process()`，`sub_pricess_A()`，`sub_process_B()`という3つのプロセスメソッドが定義されている．これらのうち，`main_process()`メソッドは，末尾の2行を除くと，Skelton2クラスの`process_method()`メソッドとほぼ同じである．なお，シミュレーション環境の`timeout()`メソッドは，`simpy.Timeout()`へのショートカットであり，引数が1つで済むためよく用いられる．

追加されている末尾の2行では，`count`の値が3で割り切れる際に，ある処理を実行していることがわかる．ここに，シミュレーション環境の`signal4A`は，`sub_process_A()`メソッドの1行目（および5行目）で生成されているEventクラスのインスタンス，すなわち事象である．そして，事象の`succeed()`メソッドは，それをtriggerするという処理を実行するものである．したがって，この箇所は，`count`が3で割り切れるたびに`signal4A`をtriggerする，という機能を果たしていることになる．

次に，`sub_process_A()`メソッドの方を見てほしい．3行目でこの事象をyieldしていることから，このメソッドはこの箇所で一時停止することになる．そして，`main_process()`メソッドの方で`signal4A`がtriggerされ，シミュレーション環境がこの事象を生起させると，`sub_process_A()`メソッドが再開される，という流れになっている．この流れは，複数のプロセス関数・メソッドを関連付けるための典型的な方法の1つである．

コード全体の下から2行目，3行目を見ると，`main_process()`メソッド，`sub_process_A()`メソッドは共に，シミュレーション開始前にシミュレーション環境に登録されていることがわかる．なお，シミュレーション環境の`process()`メソッドは，`simpy.Process()`へのショートカットであり，こちらも，引数が1つで済むためよく用いられる．

したがって，シミュレーションが始まるとこれらのプロセスは自動的に開始され，上で定めた相互作用に則って進んでいくことになる（具体的には，まず`main_process()`メソッドがスタートし，yieldまで進み，一時停止した後，`sub_process_A()`メソッドがスタートし，yieldまで進み，一時停止する．その後は，Timeout事象が生起すると，`main_process()`メソッドが再開され，その中で`signal4A`が生起すると（その次に`main_process()`メソッドが一時停止した後に）`sub_process_A()`メソッドが再開される，という具合である）．

次に，`sub_process_B()`メソッドの方を見てみよう．こちらは，whileループをもたない，単発のプロセスになっていることがわかる．このプロセスの実行はどのように制御されているのだろうか．実は，`sub_process_A()`メソッドの中にその謎が隠れている．末尾の2行を見てほしい．`count`が5で割り切れる際に，`sub_process_B()`メソッドをシミュレーション環境に登録していることがわかる．これを受けて，このプロセスが自動的に実行されるようになるわけである．このように，シミュレーション環境への新たなプロセスの登録は，シミュレーション開始前だけではなく，開始後の任意の時点でも行うことができる．この流れもまた，複数のプロセス関数・メソッドを関連付けるための典型的な方法の1つである．

## 少し発展的な話題（への入口）

### 事象のvalueとok

事象`e`は`value`という変数をもっている．`e.value`のデフォルト値は`None`であるが，それに（`None`以外の）値をセットして，プロセス関数・メソッドに渡すことができる．そのためには，事象`e`をトリガーする際に，

```python
e.succeed(valueにセットしたい値)
```

とする（Timeout事象の場合は，インスタンス生成時に，キーワード引数として，「`value=valueにセットしたい値`」のように指定する）．そして，プロセス関数・メソッド側で，yieldの箇所に，

```python
v = yied e
```

と書けば，`v`に`e.value`の値が入るという仕組みである．

さらに，事象`e`は`ok`という変数ももっている．事象`e`をtriggerするときに`succeed()`メソッドを利用すると，自動的に，`e.ok=True`となる．これは，`succeed()`メソッドの名称からもわかるように，その事象が成功裏に生起したことを表す．

実は，事象`e`をtriggerするには他にも，`e.fail(exception)`や`e.trigger(event)`というメソッドを使うこともできる．前者では，`e.ok=False`となり，その事象の生起が何らかの意味で失敗したことを示唆する．このメソッドを使った場合は，`e.value`に`exception`に指定された例外が入り，事象`e`が処理される際にその例外が発生する（ので，待ち受けしているプロセス関数・メソッドなどで例外処理を行う）．また，後者では，事象`e`の`ok`や`value`の値は，引数として渡された別の事象`event`と同じにセットされる．

### 事象待受けのあれこれ

プロセス関数・メソッドで複数の事象の論理結合を待受けすることができる．その際は，and結合には`&`，or結合には`|`をそれぞれ利用する．例えば，3つの事象`e1`，`e2`，`e3`があったとして，

```python
values = yield (e1 | e2) & e3
```

のようにできるということである．このとき，`values`は，各事象の`value`のOrderedDictになる（もちろん，各事象の`value`の値が不要であれは，「`values=`」は書かなくてよい）．

逆に，同じ事象を複数のプロセス関数・メソッドで待受けしてもよい．この場合は，その事象のコールバックの集合に（自動的に）再開命令が追加された順に，それらのプロセスが再開されていくことになる．

### Processクラスについて

プロセス関数・メソッドを登録する際に，Processクラスのインスタンスを作成していた．これを，

```python
p = simpy.Process(env, process_func())
```

などのように，後で参照できるようにしておくと便利なことがある．

実は，ProcessクラスはEventクラスを継承しているので，これも事象の一種であるとみなせる．すなわち，上の`p`を事象として扱うことができる（returnした際にtriggerされたとみなされ，戻り値があればそれが`value`の値となる）．

また，事象`p`がtriggerされる前にその`interrupt()`メソッドを呼ぶことで，対応するプロセス関数・メソッドを中断（異常終了）させることができる．これによって，そのプロセス関数・メソッドがyieldで待っている事象`e`のコールバックの集合から，対応する再開命令が削除される．また，このプロセス関数・メソッドに例外`simpy.exceptions.Interrupt(cause)`が投げ込まれるので，それを受け取って処理することによって異常終了時の挙動を指定することができる．この`interrupt()`メソッドは事象e自体には影響を与えない（ので，例外処理後にその事象`e`を再度待受けしてもよい）．

# 簡単な在庫管理の例

もう少し具体的なイメージを掴んでもらうために，最後に，前回も取り上げた簡単な在庫管理の例を示しておこう．

```python
class Model:
    def __init__(self, env, op, oq, lt, init):
        self.env = env
        self.op = op  # ordering point
        self.oq = oq  # order quantity
        self.lt = lt  # replenishment lead time
        self.at_hand = init  # how many items you have at hand
        self.loss = 0  # opportunity loss
        self.orders = []  # list of back orders

    @property
    def total(self):
        return sum(self.orders) +self.at_hand

    def print_state(self):
        print('[{}] current level: {}, back order: {}, lost sales: {} '.format(round(self.env.now), self.at_hand, self.orders, self.loss))
        self.env.log.extend()

    def seller(self):
        while True:
            yield self.env.timeout(random.expovariate(1))
            if self.at_hand > 0:
                self.at_hand -= 1  # sell an item to the customer
                self.env.stocktake.succeed()  # activate the stocktaker
            else:
                self.loss += 1  # sorry we are out of stock
            self.print_state()  # state after dealing with each customer

    def stocktaker(self):
        self.env.stocktake = self.env.event()  # create the first signal
        while True:
            yield self.env.stocktake
            if self.total <= self.op:
                self.orders.append(self.oq)
                self.env.process(self.deliverer())  # activate deliverer
            self.env.stocktake = self.env.event()  # create the next signal

    def deliverer(self):
        self.print_state()  # state after an order is placed
        yield self.env.timeout(self.lt)
        if len(self.orders) > 0:
            self.at_hand += self.orders.pop(0)
        self.print_state()  # state after an order is fulfilled
```

前回のModelクラスと比較すると，`run()`メソッド（とその他のいくつかのメソッド）が削除され，3つのプロセスメソッドが新たに定義されていることがわかる．これらのプロセスメソッドは，それぞれ，ランダムに到着する顧客に対応する販売担当者，店頭在庫量を確認して必要に応じて発注を行う在庫管理者，発注を受けて商品を配送する配送担当者の働きにぞれぞれ対応している．これらの働きが混在して記述されていた前回の`run()`メソッドと比較すると，コードの見通しが良くなったと感じられるのではないだろうか．この効果は，対象システムの規模に伴って大きくなっていく．

SImPyを導入したのに応じて，Logクラスにも少し変更を加えておこう．

```python
import matplotlib.pyplot as plt

class Log:
    def __init__(self, env):
        self.env = env
        self.time = []
        self.at_hand = []
        self.loss = []
        self.total = []
        self.extend()

    def extend(self):
        self.time.append(self.env.now)
        self.at_hand.append(self.env.model.at_hand)
        self.loss.append(self.env.model.loss)
        self.total.append(self.env.model.total)

    def plot_log(self):
        plt.plot(self.time, self.at_hand, drawstyle = "steps-post")
        plt.xlabel("time (minute)")
        plt.ylabel("number of items")
        plt.show()
```

このシミュレーションモデルを実行してみるには，下のようにすればよい．

```python
env = simpy.Environment()
env.model = Model(env, 10, 20, 10, 20)  # op, oq, lt, init
env.log = Log(env)
env.process(env.model.seller())
env.process(env.model.stocktaker())
env.run(until=200)
env.log.plot_log()
```

# 演習課題

前回の演習課題で作成した，飲食店のランチタイムの様子を表現したシミュレーションモデルを，SIｍPyを使って，プロセス志向のコードに書き換えてみよう．

# まとめ

今回は，SimPyを導入し，それを用いてプロセス志向でシミュレーションモデルを構築していく方法の基礎を紹介した．次回は，リソース関連のクラス群とその使い方についてみていこう．

# リンク

* [SimPyモジュール](https://simpy.readthedocs.io/en/latest/)
* [SimPyのコード](https://gitlab.com/team-simpy/simpy)
* [この記事のGoogle Colab Notebook](https://github.com/j54854/myColab/blob/main/Rinko_SimPy_2.ipynb)
* [第1回](https://qiita.com/j54854/items/2d98aeec4974638c24ba)