動的計画法(2) M/M/1の計算

 数式で漸化式で表現されれば、再帰関数を利用して動的計画法で計算できるということでした。今回は待ち行列のM/M/1の定常分布を動的計画法で求めていきます。

(1) 待ち行列 M/M/1の表現

式の表現は、[1]やhttps://www.win.tue.nl/~iadan/que/h4.pdfを見て欲しいのですが、リンクの式(2), (3), (4)がM/M/1の定常分布を出すための平衡方程式です。これをDPで解いていくのですが、その前にこの平衡方程式を解けば、定常分布の理論解が得られます。理論解は式(6)です。要するにここでは、理論解で得られた値とDPで得られた結果を比較して精度を確認していきます。

(2) 理論解の導出

	#M/M/1の定常分布を返す
	def getStationary(lmd, mu, n):
	pi = 1 - (lmd / mu)
	pi *= (lmd / mu ) ** n
	return pi
	#M/M/1での理論解
	lmd = 1.5
	mu = 2.0
	n = 10
	pi = getStationary(lmd, mu, n)
	print(pi)

view raw MM1理論値 hosted with ❤ by GitHub

この条件で、系内人数が10人となる確率は0.014078378677368164となります。

(2) DPでの算出

次に漸化式を利用して、再帰関数を作り、DPで同じ確率を算出します。漸化式を整理すると、３項間漸化式になり、pi(0)とpi(1)を初期値で設定します。

これをプログラムにします。

	#DPでM/M/1を解いてみる
	def getDPStationary(lmd, mu, n):
	if(n == 1):
	return (lmd / mu) * (1 - (lmd / mu))
	elif (n == 0):
	return (1 - (lmd / mu))
	else:
	return (lmd + mu)/mu * getDPStationary(lmd, mu, n-1) - (lmd / mu) * getDPStationary(lmd, mu, n-2)
	pi = getDPStationary(lmd, mu, n)
	print(pi)

view raw MM1_DP hosted with ❤ by GitHub

それなりの精度で出てきます。0.014078378677368164です。

このように漸化式であらわされているものは、再帰関数を用いて動的計画法として解いていくと便利です。

参考

[1] Bolch, G., Greiner, S., De Meer, H., & Trivedi, K. S. (2006). Queueing networks and Markov chains: modeling and performance evaluation with computer science applications. John Wiley & Sons.   (P105 Birth Death(λ_i、μ_i) P214でMM1(λ、μ)を与えている)