In this section, we introduce Sakasegawa's formula, which gives an approximation for the expected waiting time for a \(G/G/c\) queue. We discuss the insights it offers and show how to estimate the squared coefficient of variation of the inter-arrival time from data.

Sakasegawa's formula provides a reasonable approximation¹¹ There is no known exact expression for the \(G/G/c\) queue. for \(\E \W\) of the \(G/G/c\) queue, and has the form

where the utilization²² Here is a subtle point. When there are multiple machines, the utilization of each machine is not necessarily equal to \(\rho\). For instance, if there is a preference to choose the `left-most' machine whenever it is free, then the utilization of this machine is larger than \(\rho\). Only when the machines have the same speed and routing is such that each machine receives a fraction \(\lambda/c\) of jobs, all machines in one station have the same utilization. of the station is \(\rho = \lambda \E S/c \in[0, 1)\), and \(C_a^2 = \V X/(\E X )^2\) and \(C_s^2=\V S/(\E S)^2\) are the scv of the inter-arrival time and service time, respectively.

Note that Sakasegawa's formula for the \(G/G/c\) queue reduces to the PK formula for \(\E{\W}\) for the \(M/G/1\) queue; just fill in \(c=1\) and set \(C_a^2=1\). Part of the intuition behind Sakasegawa's formula is based on the following idea. To generalize the formula for \(\E \W\) from the \(M/M/1\) queue to that of the \(M/G/1\) queue, we replace \((1+1)/2\) by \((1+C_{s}^{2})/2\); the result remains exact. Now, to go from the \(M/G/1\) queue to the \(G/G/1\) queue, replace \((1+C_{s}^{2})/2\) by \((C_a^2+C_s^2)/2\); this is no longer exact, but the approximation remains good for most practical purposes. For the details behind the power with \(c\), we refer to Sakasegawa's paper.

It is crucial to memorize the insights this formula offers. First, we see that \(\E{\W} \sim (1-\rho)^{-1}\). Consequently, when \(\rho\) is large, the waiting time is (very) long. And, not only is \(\E\W\) large, it is also extremely sensitive to the actual value of \(\rho\). Clearly, such situations must be avoided and therefore, when trying to improve a queueing system, the first focus should be on reducing \(\rho\).

Second, \(\rho = \lambda \E S/c\). Thus, when \(\rho\) must be made smaller, we have just three options.³³ And not more! We can reduce the arrival rate \(\lambda\) of jobs, for instance, by blocking demand or sending it elsewhere, such as to another machine.⁴⁴ And if customers have a choice, they will take their own measures, simply by going elsewhere. We can make \(\E S\) smaller by replacing a server with a faster one or by shifting some of the processing steps of a job to other servers, thereby making job sizes smaller. Finally, we can add servers, i.e., making \(c\) larger.⁵⁵ This is precisely what you see in supermarkets. If technically possible, adapting \(c\) to \(\lambda\) when \(\E S\) cannot be easily changed is a very effective mechanism for controlling waiting times.

Third, \(\E{\W} \sim C_a^2\) and \(\E{\W} \sim C_s^2\), which implies that when the job inter-arrival or service times are very variable, \(\E\W\) is large. Thus, after ensuring that \(\rho\) is sufficiently small, it becomes important to concentrate on reducing \(C_a^2\) and \(C_s^2\).⁶⁶ It works the other way too: systems with low variability can deal with higher load.

Finally, \(\E{\W} \sim \E S/c\). This means that from the perspective of a job in queue, average job service times are \(c\) times as short as `their own' service times.

Clearly, Sakasegawa's equation requires an estimate⁷⁷ If you are not interested in statistics, you can skip this part. of \(C_a^2\) and \(C_s^2\). Now it is not always easy in practice to determine the actual service-time distribution, since service times are often estimated by a planner rather than measured. Similarly, the actual arrival moments of jobs are often not registered, only the date, or perhaps the hour, that a customer arrived. Hence, it is often not possible to directly estimate \(C_a^2\) and \(C_s^2\) from the available information.

However, when the number of arrivals per period \(\{a_n\}\) has been logged for some time, we can use \(\{a_n\}\) to estimate \(C_a^2\) as⁸⁸ This can of course also be used to estimate \(C_s^2\).

where

We derive this result in steps.⁹⁹ The arguments are based on a book by Cox on renewal theory. If you don't like to extend your knowledge of statistics, you can skip this derivation.

First, recall some results from earlier sections. Let the inter-arrival time \(X\) have mean \(1/\lambda\) and variance \(\sigma^2\), so that

Let \(A_k\) be the arrival time of the \(k\)th arrival and \(A(t)\) the number of arrivals up to time \(t\). Consider the following useful relation between \(A(t)\) and \(A_k\), see Ex 1.2.7,

Since the inter-arrival times have finite mean and second moment by assumption, we use the central limit theorem to see that

where \(N(0,1)\) is a standard normal rv with distribution \(\Phi(\cdot)\). Similarly, for an \(\alpha\) yet to be determined,¹⁰¹⁰ This is a common trick: suppose that there exists a constant to take of the scaling, and then try an extra relation to `ferret out' the scaling constant.

Thus, \(\E{A(t)} = \lambda t\) and \(\V{A(t)}=\alpha^2 t\).

Using that \(\P{N(0,1) \leq y} = \P{N(0,1) > -y}\), we get

We can use this relation between the distributions of \(A(t)\) and \(A_k\) to see that \(\P{A_k >t_k } = \P{A(t_k) <k}\) where we define for ease

With this we get,

Since \((A(t_k) - \lambda t_k)/ \alpha \sqrt{t_k} \to N(0,1)\) as \(t_k \to \infty\), the above implies that

as \(t_k \to \infty\). Using the above definition of \(t_k\), the LHS of this equation can be written as

Since \(t_k \to \infty\) is implied by (and implies) \(k\to\infty\), we therefore want that \(\alpha\) is such that

as \(k\to\infty\). This is precisely the case when

Finally, for \(t\) large, or, by the same token, \(k\) large,

where the last equation follows from the above definition of \(C_a^2\).

a = 3.0
b = 9.0
EX = (b + a) / 2.0  # expected inter-arrival time
labda = 1.0 / EX  # per minute
VA = (b - a) * (b - a) / 12.0
CA2 = VA / (EX * EX)
print(f"{EX=}, {labda=:.4f}, {CA2=:.4f}")

ES = 5.0
VS = 4
CS2 = VS / (ES * ES)
rho = labda * ES
print(f"{CS2=:.4f}, {rho=:.4f}")

W = (CA2 + CS2) / 2.0 * rho / (1.0 - rho) * ES
print(f"{W=:.4f}")

EX=6.0, labda=0.1667, CA2=0.0833
CS2=0.1600, rho=0.8333
W=3.0417

W = (1.0 + CS2) / 2.0 * rho / (1.0 - rho) * ES
print(f"{W=:.4f}")

W=14.5000