mog project: 2011.11

Python: M/M/c 待ち行列問題を解く

待ち行列理論の M/M/c モデル(＝M/M/c/∞モデル、M/M/1モデルの一般化)の問題を Python で解く。
コードを簡単にするため、ビルトインの float をオーバーライドして0除算を可能にしている。

入力：トランザクション数 {λ∈0以上の浮動小数点数}　（＝平均到着間隔 Ta の逆数）
　　　　平均サービス時間 {Ts∈0以上の浮動小数点数}　（＝平均サービス率 μ の逆数）
　　　　サービス窓口数 {c∈1以上100以下の自然数}

出力：arr(arrival rate) ＝トランザクション数＝ λ
　　　dep(departure rate) ＝平均サービス率＝ μ
　　　c ＝サービス窓口数
　　　util ＝サービス利用率＝ ρ ＝ λ / cμ
　　　TmQue ＝平均待ち行列時間
　　　TmSys ＝平均応答時間（ターンアラウンドタイム）
　　　NmQue ＝待ち行列の長さの平均
　　　NmSys ＝待ち系の長さの平均　　　

・mmc_queue.py

# Copyright (c) 2011 Mog Project. All rights reserved.
# -*- coding: utf-8 -*-
"""
In the mathematical theory of random processes, the M/M/c queue is
a multi-server queue model. It is a generalisation of the M/M/1 queue.

    Arrivals are a Poisson process
    Service time is exponentially distributed
    There are c servers
    The length of queue in which arriving users wait before being served is
        infinite
    The population of users (i.e. the pool of users), or requests, available
        to join the system is infinite

This type of system arises in many situations, including lines at a bank,
    phone queueing systems, the application of computer resources, etc...

http://en.wikipedia.org/wiki/M/M/c_queue
"""
import math
from decimal import Decimal

INF = float('inf')
class Float(float):
  def __div__(self, other):
    if other:
      return float.__div__(self, other)
    elif self == 0.0 or math.isnan(self):
      return float('nan')
    return float('inf') if self > 0.0 else float('-inf')

class InputError(Exception): pass
class MMcQueue(object):
  def __init__(self, arrival_rate, service_time, num_servers=1):
    """
    @param arrival_rate: 0 and over
    @param service_time: 0 and over
    @param num_servers: int between 1 and 100, inclusive
    """
    self.c = num_servers  # Number of servers.
    if not (isinstance(num_servers, int) and 1 <= num_servers <= 100):
      raise InputError('num_servers must be int between 1 and 100, inclusive.')
    if arrival_rate < 0 or service_time < 0:
      raise InputError('arrival_time and service_time must be 0 and over.')

    self.arrival_rate = Float(arrival_rate)
    service_time = Float(service_time)
    self.service_rate = Float(1.0) / service_time
    self.util = self.arrival_rate / self.service_rate / self.c

    if self.util == 0.0:
      self.tw = self.lq = self.lw = 0.0
      self.tq = service_time
      return
    if self.util >= 1.0:
      self.tw = self.tq = self.lq = self.lw = INF
      return

    crc = pow(self.c * self.util, self.c)
    factc = math.factorial(self.c)
    x = 0.0
    for i in range(self.c):
      x += pow(self.c * self.util, i) / math.factorial(i)
    x += crc / (factc * (1.0 - self.util))
    x = crc / x / (factc * (1.0 - self.util))

    # Expected number of customers waiting in queue – not in service
    self.lw = x * self.util / (1.0 - self.util)
    # Average waiting time in queue
    self.tw = Float(self.lw) / self.arrival_rate
    # Average time in system (queued + serviced)
    self.tq = self.tw + service_time
    # Expected number of customers in system
    self.lq = self.arrival_rate * self.tq
    return

  def __str__(self):
    return ('arr=%.3f, dep=%.3f, c=%d, util=%.3f, ' +
        'TmQue=%.3f, TmSys=%.3f, NmQue=%.3f, NmSys=%.3f') \
        % (self.arrival_rate, self.service_rate, self.c, self.util,
        self.tw, self.tq, self.lw, self.lq)

def main():
  # boundary values
  print MMcQueue(0, 0)
  print MMcQueue(1, 0)
  print MMcQueue(0, 1)
  print MMcQueue(1, 1)
  print MMcQueue(INF, 0)
  print MMcQueue(INF, 1)
  print MMcQueue(0, INF)
  print MMcQueue(1, INF)
  print MMcQueue(INF, INF)
  print MMcQueue(1, 99, 100)

  # error cases
#  print MMcQueue(-1, 1)
#  print MMcQueue(1, -1)
#  print MMcQueue(1, 1, 0)
#  print MMcQueue(1, 1, 101)

  # loop sample
  i = Decimal(0.0)
  while i <= 2.0:
    print MMcQueue(i, 10, 15)
    i += Decimal(0.1)
  for i in range(1, 10):
    print MMcQueue(2, 2, i)

if __name__ == '__main__':
  main()

・出力結果

arr=0.000, dep=inf, c=1, util=0.000, TmQue=0.000, TmSys=0.000, NmQue=0.000, NmSys=0.000
arr=1.000, dep=inf, c=1, util=0.000, TmQue=0.000, TmSys=0.000, NmQue=0.000, NmSys=0.000
arr=0.000, dep=1.000, c=1, util=0.000, TmQue=0.000, TmSys=1.000, NmQue=0.000, NmSys=0.000
arr=1.000, dep=1.000, c=1, util=1.000, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=inf, dep=inf, c=1, util=nan, TmQue=nan, TmSys=nan, NmQue=nan, NmSys=nan
arr=inf, dep=1.000, c=1, util=inf, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=0.000, dep=0.000, c=1, util=nan, TmQue=nan, TmSys=nan, NmQue=nan, NmSys=nan
arr=1.000, dep=0.000, c=1, util=inf, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=inf, dep=0.000, c=1, util=inf, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=1.000, dep=0.010, c=100, util=0.990, TmQue=87.394, TmSys=186.394, NmQue=87.394, NmSys=186.394
arr=0.000, dep=0.100, c=15, util=0.000, TmQue=0.000, TmSys=10.000, NmQue=0.000, NmSys=0.000
arr=0.100, dep=0.100, c=15, util=0.067, TmQue=0.000, TmSys=10.000, NmQue=0.000, NmSys=1.000
arr=0.200, dep=0.100, c=15, util=0.133, TmQue=0.000, TmSys=10.000, NmQue=0.000, NmSys=2.000
arr=0.300, dep=0.100, c=15, util=0.200, TmQue=0.000, TmSys=10.000, NmQue=0.000, NmSys=3.000
arr=0.400, dep=0.100, c=15, util=0.267, TmQue=0.000, TmSys=10.000, NmQue=0.000, NmSys=4.000
arr=0.500, dep=0.100, c=15, util=0.333, TmQue=0.000, TmSys=10.000, NmQue=0.000, NmSys=5.000
arr=0.600, dep=0.100, c=15, util=0.400, TmQue=0.002, TmSys=10.002, NmQue=0.001, NmSys=6.001
arr=0.700, dep=0.100, c=15, util=0.467, TmQue=0.008, TmSys=10.008, NmQue=0.005, NmSys=7.005
arr=0.800, dep=0.100, c=15, util=0.533, TmQue=0.028, TmSys=10.028, NmQue=0.022, NmSys=8.022
arr=0.900, dep=0.100, c=15, util=0.600, TmQue=0.080, TmSys=10.080, NmQue=0.072, NmSys=9.072
arr=1.000, dep=0.100, c=15, util=0.667, TmQue=0.204, TmSys=10.204, NmQue=0.204, NmSys=10.204
arr=1.100, dep=0.100, c=15, util=0.733, TmQue=0.474, TmSys=10.474, NmQue=0.522, NmSys=11.522
arr=1.200, dep=0.100, c=15, util=0.800, TmQue=1.064, TmSys=11.064, NmQue=1.277, NmSys=13.277
arr=1.300, dep=0.100, c=15, util=0.867, TmQue=2.478, TmSys=12.478, NmQue=3.222, NmSys=16.222
arr=1.400, dep=0.100, c=15, util=0.933, TmQue=7.223, TmSys=17.223, NmQue=10.112, NmSys=24.112
arr=1.500, dep=0.100, c=15, util=1.000, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=1.600, dep=0.100, c=15, util=1.067, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=1.700, dep=0.100, c=15, util=1.133, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=1.800, dep=0.100, c=15, util=1.200, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=1.900, dep=0.100, c=15, util=1.267, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=2.000, dep=0.500, c=1, util=4.000, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=2.000, dep=0.500, c=2, util=2.000, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=2.000, dep=0.500, c=3, util=1.333, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=2.000, dep=0.500, c=4, util=1.000, TmQue=inf, TmSys=inf, NmQue=inf, NmSys=inf
arr=2.000, dep=0.500, c=5, util=0.800, TmQue=1.108, TmSys=3.108, NmQue=2.216, NmSys=6.216
arr=2.000, dep=0.500, c=6, util=0.667, TmQue=0.285, TmSys=2.285, NmQue=0.570, NmSys=4.570
arr=2.000, dep=0.500, c=7, util=0.571, TmQue=0.090, TmSys=2.090, NmQue=0.180, NmSys=4.180
arr=2.000, dep=0.500, c=8, util=0.500, TmQue=0.030, TmSys=2.030, NmQue=0.059, NmSys=4.059
arr=2.000, dep=0.500, c=9, util=0.444, TmQue=0.010, TmSys=2.010, NmQue=0.019, NmSys=4.019

参考：
http://en.wikipedia.org/wiki/M/M/c_queue
http://tomari.org/main/java/machim.html

mog project

11.30.2011

Disk failures in the real world

11.28.2011

One-liner IsPrime & GCD in C++

11.20.2011

Compare regions in the same file with Sakura Editor

11.13.2011

HttpWatch on Firefox

11.12.2011

M/M/c (M/M/1) queue with Python

11.10.2011

Useful techniques for AWK

11.03.2011

How to break out of all loops in Python