Queueing Theory

A.I. McLeod, University of Western Ontario
Copyright ©, A.I. McLeod, 1999.

The object of queueing theory is to model and quantify the degree of congestion in a system such as customers waiting at a bank for service. These waiting line type systems have numerous applications in a wide variety of fields such as telephone and computer engineering, health care planning and management science.

Standard M/M/s Queue

In this queueing system, inter-arrival times and service times of each of s are independently exponentially distributed with mean parameters (λ^(-1)) and (μ^(-1)) respectively. It is assumed that the queue discipline is FIFO (first-in, first-out) and that there is no limit on the size of the queue. A customer enters the system and joins the queue if there is one and when his turn comes, the customer is served by the first available server. When this service is completed, the customer leaves the system. One measure of the activity in the system is the traffic intensity, ρ  =  λ/(μ s). Provided that ρ < 1, the system will eventually settle down to a steady-state no matter what initial condition is specified. It can be shown that the steady-state probability distribution, (Subscript[P, n])  =  Pr {N(t) = n}  are given by

\(TraditionalForm\`P\_0\  = \ 
        1/\([∑\+\(n = 0\)\%\(s - 1\)\((λ/μ)\)\^n\/\(n!\)\)\)+((((λ/μ))^s)/(s!)) 
          (1/(1 - ρ)) ]

(Subscript[P, n]) = {({{((((λ/μ))^n)/(n!)) (Subscript[P, 0]),  
                n ≤ s, }, {((((λ/μ))^n)/(s! (s^(n - s)))) 
                  (Subscript[P, 0]) ,  n ≥ s.}})

The steady state probabilities can be interpreted as the fraction of time the system spends in state n. In some applications such as telephony engineering the probability of queueing  (Subscript[P, q]) = 
        1 - ((Subscript[P, 0]) +  ...   + (Subscript[P, s - 1]))\(TraditionalForm\`\(\_\ \)\)is an important system performance characteristic. 
There are four other common performance measures: (Subscript[L, q]), (Subscript[W, q]), L, W. First,  the average number of customers waiting in the queue,

(Subscript[L, q]) = (∑_(n = s)^∞)(n - s) 
              (Subscript[P, n])=\(TraditionalForm\`\(\((λ/μ)\)\^s\/\(s!\)\) \(ρ\/\((1 - ρ)\)\^2\) P\_0\).

Second, the average time spend waiting in the queue until service starts is derived from Little's identity, (Subscript[W, q]) = (Subscript[L, q])/λ. Third, the average time spent in the system, W = (Subscript[W, q]) + (λ^(-1)). Fourth, the average number of customers in the system,

L = (∑_(n = 0)^∞) n (Subscript[P, n]) =  
          (Subscript[L, q]) + (λ/μ).

Here is a function to compute these six performance measures: ρ, (Subscript[P, q]),  (Subscript[L, q]),  (Subscript[W, q]), L, W.  As a byproduct of this function, another function P[ k ] is defined. P[ k ] is our steady-state probability (Subscript[P, k]).

\(MMs[λ_, μ_, \ s_]\  := \ 
      Module[{L, W, Lq, \ Wq, \ Pq, \ ρ, \ 
          r = λ/μ}, ρ = λ/\((s\ μ)\); 
		If[ρ ≥ 1, \ Print["Error: ρ≥1"]; 
          Return[{}]]; 
		P0 = 
          1/\((r^s/\((\(s!\)\ \ \((1 - ρ)\))\) + 
                Apply[Plus, \ 
                  Table[r^n/\(n!\), {n, 0, s - 1}]])\); 
		Clear[
          P]; 
		P[
            k_ /; \((k ≥ 0\  && \ 
                  k ≤ s)\)]\  := \ \((λ/μ)\)^
                k/\((\(k!\))\)\ P0; 
		P[
            k_ /; k > s]\  := \ \((λ/μ)\)^
                k/\((\(s!\)\ s^\((k - s)\))\)\ P0; 
		Pq = \ 
          1 - Apply[Plus, \ Table[P[k], {k, 0, s - 1}]]; 
		Lq = 
          r^s\ ρ\ P0/\((\(s!\)\ \((1 - ρ)\)^2)\); 
		Wq = 
          Lq/λ; \ W = Wq + 1/μ; \ 
        L = W\ λ; 
		{"ρ" -> ρ, "Pq" -> 
            Pq, \ "L" -> L, \ "W" -> W, \*""(Subscript[L, q])"" -> 
            Lq, \ \*""(Subscript[W, q])"" -> Wq\ }
		]\)

MMs[2, 3, 1]

MMs[2, 3, 2]

Plot[Wq /. MMs[x, 1, 1], {x, 0.1, 0.9},  
        PlotPoints -> 100];

Simulation Analysis: Standard  Case

Let N(t) be the number in the system at time t. So N(t) is a step function with increments when new arrivals occur and decrements when service completions occur. The total observed time spent in the system over the interval  [ 0, t]  is

T = (∫_0^t) N(
              s) ds.

So the average time spent in the system over the interval  [0,  t]  is W  =  T/A(t), where A(t) is the number of arrivals over the interval [0, t]. Assuming equilibrium of the system, we can apply Little's identity to deduce the average number in the system over the [0, t], L  =  W λ, where λ is the average arrival rate. In other words the mean interarrival time is (λ^(-1)). Let (μ^(-1)) be the average service rate then the average time a customer spends waiting until they receive service is (Subscript[W, q]) =  W - (μ^(-1)) and the mean number waiting in the queue is (Subscript[L, q]) = (Subscript[W, q]) λ. By computing the fraction of time spent in each state over the simulation we can obtain estimates of the steady state probabilities, (Subscript[P, k]).
As discussed in Lewis & Orav (1989, §3.2, §3.3) it is also possible to obtain the empirical distributions for the waiting times and queue lengths. However for many applications the average values are sufficient. By doing a very large number of simulations we avoid two problems. First, after a large number of iterations, the sample means have converged so we don't have to worry about trying to estimate their standard error or precision. We can check if this is the case simply by repeating our simulation again and seeing if there is any noticeable difference. Second, we bypass the issue of transient behaviour. Often simulation studies delete the first 100 iterations to remove the effect of startup. However by doing a large number of iterations, this problem also disappears. Again, if there is concern we can repeat the simulation trying various initial configurations.
With simulation analysis we can of course solving queueing problems that are not otherwise tractable. First we look at the standard case and then extend to an intractable case where the arrivals are exponential but autocorrelated. In the standard case we make the assumption that the queue length is unlimited, the queue discipline is FIFO, there are s servers, the inter-arrival times are independent with rate λ and the service times are exponential and independent with rate μ.

Needs["Statistics`ContinuousDistributions`"];

ArrivalTime[] := Random[ExponentialDistribution[2]];
ServiceTime[]  :=  
        Random[ExponentialDistribution[3]];

\(\(InitializeSystem[servers_:  1]\  := 
        Module[{}, 
	 (*\ initialization\ of\ queueing\ system\ *) 
	 (*\ 
            Nt\  = \ 
              number\ in\ the\ system\ when\ an\ event\ occurs\ *) 
	 (*\ 
            TIME\  = \ 
              the\ time\ since\ the\ system\ initialization\ *) 
	 (*\ 
            At\  = \ 
              number\ of\ arrivals\ to\ the\ system\ so\ far\ *) 
	 (*\ 
            i\  = \ counter, \ 
            number\ of\ events\ *) 
	 (*\ Λ\  = \ 
              empirical\ mean\ of\ inter - 
                arrival\ times\ *) 
	 (*\ Σ\  = \ 
              empirical\ mean\ of\ service\ duration\ times\ *) 
	 (*\ 
            f\  = \ number\ of\ arrivals\ who\ receive\ immediate\ service\ *) 
	 (*\ s\  = \ number\ of\ servers\ *) 
	 (*\ 
            X\  = \ time\ of\ next\ arrival\ *) 
	 (*\ 
            Y\  = \ service\ times\ remaining\ for\ those\ receiving\ service\ *) 	
	 (*\ n[x]\  = \ total\ time\ system\ is\ in\ state\ x, \ 
            x\  = \ 0, \ 1, \ 2, \  ... *) 
		Clear[n]; 
          Nt\  = \ \(TIME\  = \ \(T = \(At = \(i = \(Λ = \(Σ = \(f\_\(\(\ \)\(=\)\)\) \(n[x_] = 0\)\)\)\)\)\)\); 
          s = servers; X = ArrivalTime[]; \ 
          Y = {};
		]\ ;\)\ \ \ \)
\(\(NextEvent[] := \ 
		Which[\ 
          Nt == 0, \  (*\ if\ the\ system\ is\ empty, \ 
            next\ event\ is\ an\ arrival\ *) \ 
		\ T = T + Nt*X; \ 
          n[Nt] += X; \ TIME = TIME + X; At = At + 1; \ Nt = Nt + 1; 
          X = ArrivalTime[]; Λ = \((Λ\ i + 
                  X)\)/\(++i\); 
	\ \ \ Y = {Ym = 
                ServiceTime[]}; Σ = \((Σ \((i - 
                        1)\) + Ym)\)/i; \(\(++f\);\), 
		 (*\ 
            if\ next\ event\ is\ an\ arrival\ *) 
		X ≤ \((Ym = 
                Min[Y])\), \ T = T + Nt*X; \ n[Nt] += X; TIME = TIME + X; \ 
          Y -= X; 
	\(At++\); \ \(Nt++\); 
          X = ArrivalTime[]; Λ = \((Λ\ i + 
                  X)\)/\(++i\); 
			 (*\ if\ Nt ≤ s, \ 
            service\ starts\ immediately\ since\ at\ this\ point\ Nt\ has\ already\ been\ incremented*) 
		If[Nt ≤ s, \ \(++k\); 
            AppendTo[Y, \ 
              Ym\  = \ 
                ServiceTime[]]; \ Σ = \((Σ \((i - 
                          1)\) + Ym)\)/i], 
		 (*\ 
            next\ event\ must\ be\ a\ service\ completion\ *) 
			X > \((Ym = Min[Y])\), 
		\ T = T + Nt*Ym; n[Nt] += Ym; 
          TIME = TIME + Ym; \(Nt--\); X = X - Ym; \(i++\); 
          Y = Drop[Y, First[Position[Y, Ym]]]; Y -= Ym; 
          If[Nt ≥ s, 
            AppendTo[Y, \ 
              Ym\  = \ 
                ServiceTime[]]; \ Σ = \((Σ\ \((i - 1)\) + Ym)\)/i];
						];\)		\)
\(W[] := \ T/At; Wq[] := W[] - Σ; \ 
    L[] := W[]/Λ; \ Lq[] := Wq[]/Λ\ ; 
    Pr[k_] := n[k]/TIME;	\)
\(\(QueueSummary[] := 
        Module[{}, 
		Print["W  = ", \ 
            W[], \ \*""                              (Subscript[W, q]) = "", \ 
            Wq[]]; 
		Print["L  = ", \ 
            L[], \ \*""                              (Subscript[L, q]) = "", \ 
            Lq[]]; 
		Print["λ = ", \ 
            1/Λ, \ \*""                              (λ^(-1)) = "", \ Λ]; 
		Print["μ = ", \ 1/Σ, \ \*""                              (μ^(-1)) = "", \ Σ]; 
		Print["Probability of Queueing = ", \ 1.0 - k/At]; 
		Print["Total Simulated Time = ", \ 
            TIME]; 
		Print["Total Number of Iterations = ", \ 
            i]
	];\)	\)

InitializeSystem[1];

Do[NextEvent[], {20000}]

QueueSummary[]

Table[{k, Pr[k]}, {k, 0, 10}] // TableForm

InitializeSystem[2];

Do[NextEvent[], {20000}]

QueueSummary[]

The reader can verify that these averages do not change if we repeat the computation again with new random numbers.

Table[{k, Pr[k]}, {k, 0, 10}] // TableForm

Simulation Analysis: Autocorrelated Arrivals

We will use a model for autocorrelated exponentially distributed random variables was proposed by Tavares (1980).  For any  \(TraditionalForm\`ϕ, \ \ \(\(|\)\(\ \)\(\(ϕ | \(\(<\)\(1\)\)\)\(,\)\)\)\)  let the sequence of random variables (Subscript[X, t]) be defined by

(Subscript[X, t])={({{\(X\_0, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ t = 0, \)}, {min((Subscript[X, t - 1]), (Subscript[Y, t]))/ϕ,       t > 0, }})

where (Subscript[X, 0]) is exponentially distributed with rate parameter λ and (Subscript[Y, t]), t = 
        1, 2,  ...  are a sequence of independent and identically distributed exponential random variables with parameter λ (1 - ϕ)/ϕ. Then it can be shown (Tavares, 1980) that (Subscript[X, t]) is exponentially distributed with rate parameter λ for all t ≥ 0 and the sequence defines a first-order autoregressive process with autocorrelation function given by ρ(k)  =  (ϕ^k),  
      k ≥ 0,  where ρ(
          k) = γ(k)/γ(0), γ(k)  =  Cov {(Subscript[X, t]),  (Subscript[X, t + k])),  k ≥ 0.  When ϕ is large, the congestion in the system is tremendously increased. This makes sense a priori since the effect of strong positive autocorrelation causes the arrivals to occur in clusters. Similarly we would expect negative autocorrelation to decrease the congestion in the system.

Clear[ArrivalTime,  ServiceTime];
ServiceTime[] := Random[ExponentialDistribution[3]]; ϕ = 0.9;  
    Xt = Random[ExponentialDistribution[2]]; 
    ArrivalTime[]  := Xt =  
        Min[Xt, Yt = 
              Random[ExponentialDistribution[
                  2 (1 - ϕ)/ϕ]]]/ϕ;

InitializeSystem[1];

Do[NextEvent[], {20000}]

QueueSummary[]

InitializeSystem[2];

Do[NextEvent[], {20000}]

QueueSummary[]

Once again, it may be shown that these averages don't change appreciably if we repeat the computation again. We can compare the steady state distributions for the standard case and the autocorrelated case using a barchart.

Needs["Graphics`Graphics`"];

MMs[2, 3, 2];

BarChart[Table[P[k], {k, 0, 6}],  Table[Pr[k], {k, 0, 6}],  
        PlotRange -> All];

Recall P[k] and Pr[k] are steady-state probabilities from the functions MMs and QueueSummary. This plot shows that the fraction of time the system is idle is slightly larger in the autocorrelated case. Note that the probability of queueing does not equal 1 - (Subscript[P, 0]) - (Subscript[P, 1]) as in the independent case. To see that this is reasonably consider a system in operation for 24 hours (1440 minutes) with s = 1 servers and with bulk scheduled arrivals once per day. Each bulk is comprised of 100 customers who require 1 minute for service. Then the probability of queueing is 99/100 even though (Subscript[P, 0]) = 1340/1440. 
As a check, let's calculate L, the average number in the system, directly from the simulation probabilities,

Sum[k*Pr[k], {k, 0, 100}]