待ち行列M/M/1での平均系内客数を求めるプログラム

Question

待ち行列M/M/1モデルでの平均系内客数の求め方がわかりません。

Python : Jupyter notebookで待ち行列M/M/1モデルでの平均系内客数を求めるプログラムを書いているのですが、計算結果が理論値とかけ離れてしまいます。エラーメッセージは出なかったので根本的な計算のアルゴリズムが違うのだと思います。どこが間違っているのか、ご指導いただけたら幸いです。

該当のソースコード

import simpy
import numpy as np
import random

#客の到着率、サービス率は指数分布に従う

env = simpy.Environment()

time = [0] #客が出入りしたときの時刻
stay = 0 #客の滞在人数
stay_time = 0 #滞在人数×滞在時間

# 到着イベント
def arrive():
    while True:
        global stay, stay_time
        yield env.timeout(random.expovariate(0.8 )) # 平均到着率 
        env.process(queue())

# 待ち行列に並ぶ
def queue():
    global time, stay, stay_time
    time.append(env.now)
    stay_time += stay * (time[1] - time[0])
    stay += 1 #店に到着した
    
    
    while len(time) > 1:
        time = time[1:]
        yield env.timeout(random.expovariate(1.0)) # 平均サービス率
        time.append(env.now)
        stay_time += stay * (time[1] - time[0])
        stay -= 1 #店から出た
        time = time[1:]
        
# 実行
env.process(arrive())
env.run(until = 1000000)

print(stay_time / 1000000)

試したこと

平均到着率(単位時間に来る客の人数)をλ、平均サービス率(単位時間にサービスを提供する人数)をμ、ρ=λ/μとすると、平均系内客数の理論値はρ/(1-ρ)で求まります。例えばρ=0.8でプログラムを実行したところ、理論値４に対して実行結果が0.799と全く違う値となりました。他にも色々なρの値で試したのですが、ρとほぼ同じ値がそのまま実行結果として出てきました。

補足情報（FW/ツールのバージョンなど）

平均系内客数は以下の画像のやり方で求めようと試みました。この式をプログラムにしたつもりなのですが、一体どこが間違っているのでしょうか。皆様よろしくお願いします。

Answer 1

import simpy
import numpy as np
import random

# 到着イベント
def arrive():
    global time, stay, canserve
    while True:
        yield env.timeout(random.expovariate(0.8)) # 平均到着率 
        time.append(env.now)
        stay += 1
        if(canserve):
            env.process(queue())

# 待ち行列に並ぶ
def queue():
    global stay, stay_time, time, canserve
    canserve = False
    while len(time) > 0:
        yield env.timeout(random.expovariate(1.0)) # 平均サービス率
        stay_time += env.now - time[0]
        time = time[1:]
        stay -= 1 #店から出た
    canserve = True
    
time = [0] #客が出入りしたときの時刻
stay = 0 #客の滞在人数
stay_time = 0 #滞在人数×滞在時間
canserve = True
    
#客の到着率、サービス率は指数分布に従う
env = simpy.Environment()

# 実行
env.process(arrive())
env.run(until = 1000000)

#print(np.average(stay_time))
print(np.sum(stay_time)/1000000)

4.012869628764746

Comments

[tarouyamadaaaaa]

ありがとうございます！二つ質問があるのですがよろしいでしょうか？
1. 変数conserveはどのような役割を持つのでしょうか？
2. 画像の方法では時間の区間×人数を積算していますが、Penpen7様のプログラムでは時間の区間のみを積算しています。なぜこの方法でうまくいくのでしょうか？
重ね重ね申し訳ありませんが、よろしくお願いします

[Penpen7]

1. サンプルコード上では、
while len(time) > 0:
yield env.timeout(random.expovariate(1.0)) # 平均サービス率
stay_time += env.now - time[0]
time = time[1:]
stay -= 1 #店から出た
の処理を重複して実行しないようにしていたため、実行している間はFalseにし、実行していない時はTrueにすることで、同じ機能を実装しました。(can_service_beginという名前がいいかもしれません)
2. 平均系内客数の式の総和部分は上のグラフの積分値です。つまり面積値ということです。私がやったのは、全面積を求めるための分割方法を変えただけです。どのような分割方法かというと、ある人が入った時間と出た時間の差を横の長さ、縦の長さを1とした長方形で分けてやるものです。グラフに補助線を引いて考えてみてください。そうすると、人数は不要になり、単にある人が到着した時間と出た時間を引いて足すだけで面積を求めることができます。

[tarouyamadaaaaa]

ありがとうございます。助かりました

[Penpen7]

よかったです。この質問は解決済みとしてクローズしていただけますか？

[tarouyamadaaaaa]

承知しました