In Section 2.1 we modeled time as progressing in discrete chunks. However, we can also model queueing systems in continuous time, so that jobs can arrive at any moment in time and have arbitrary service times. In this section, we develop a set of recursions to construct the waiting times of jobs served in their order of arrival.¹¹ That is, according to the FIFO scheduling rule.

Let us imagine that a machine²² A single server. starts working on a one-hour job at 9 am, and there is no job in the queue. When the next job arrives before \(10\) am, this second job has to wait in the queue until the first job finishes at 10 am. If, however, this second job arrives after 10 am, it finds the server free, so its service can start immediately upon arrival.

More generally, suppose that we are given an ordered sequence of arrival times³³ Moments at which the system state is inspected or changes are often called epochs. \(\{A_{k}\}\) and a set of iid service times \(\{S_{k}\}\), so that job \(k\) arrives at time \(A_{k}\) and needs an amount \(S_{k}\) of service. The departure times \(\{D_{k}\}\) follow from the recursion⁴⁴ Note that we assume that job \(k\) `reveals' its required service at the moment it arrives.

\begin{align*} \As_{k}&= \max\{A_k, D_{k-1}\}, & D_{k} &= \As_{k} +S_k, & D_{0} &=0, \tag{3.1.1} \end{align*}

where \(\As_{k}\) is the time the service of job \(k\) starts⁵⁵ That is, when this job moves from the queue to the server. , and \(D_{0}=0\) is meant to start the recursions. Why is this true? The service of job \(k\) cannot start before it arrives at \(A_k\) or before job \(k-1\) leaves the server at \(D_{k-1}\). Thus, the first moment the service of job \(k\) can start is \(\max\{A_k, D_{k-1}\}\). When started, it takes \(S_{k}\) to get served.

Once \(A_{k}\) and \(D_{k}\) are known, we compute the sojourn time, i.e., the total time in the system, and the waiting time, i.e., the time in queue⁶⁶ But not at the server. , as

\begin{align*} \J_k &= D_{k} - A_k, & \W_k &= \As_k - A_k = \J_{k} - S_{k}. \end{align*}

While it is often simple to measure arrival times in practical situations, it is harder to generate \(\{A_{k}\}\) directly with a simulator. For this we let a random number generator generate iid rvs \(\{X_{k}\}\) that represent the inter-arrival times between jobs. The arrival times can then be computed recursively from Eq. (1.2.2) as \(A_k = A_{k-1} + X_k\), \(A_{0} = 0\). And conversely, from Eq. (1.3.1), if \(\{A_{k}\}\) is known, \(X_{k} = A_{k}-A_{k-1}\). Interestingly, if \(\{X_{k}\}\) and \(\{S_{k}\}\) are given, we need not first compute \(\{A_{k}\}\) and \(\{D_{k}\}\) using Eq. (3.1.1) to find the sojourn times and waiting times. In fact, suppose that job \(k-1\) has to wait a time \(\W_{k-1}\) in queue and spends \(S_{k-1}\) in service, then its sojourn time \(\J_{k-1} =\W_{k-1}+ S_{k-1}\). If \(X_k\) time units elapse between the arrivals of job \(k-1\) and \(k\), then job \(k\) must wait in queue⁷⁷ As we start counting jobs for \(k=1, 2, \ldots\), \(\W_{0}\) has no physical meaning.

\begin{align*} \W_{k} &= [\W_{k-1} + S_{k-1}- X_k]^+, & \W_{0} &=0, & \J_{k} &= \W_{k} + S_{k}. \tag{3.1.2} \end{align*}

However, if only \(\J_k\) is needed, use \(\J_{k} = [\J_{k-1} - X_k]^+ + S_k\), \(\J_{0} =0\).

Later we will need the system state at arbitrary times, not only at arrivals. The number of arrivals \(A(t)\) during \([0, t]\) can be computed from \(\{A_{k}\}\) as⁸⁸ Read the RHS as the maximum of all \(k\) such that \(A_{k} \leq t\).

\begin{equation*} A(t):= \sum_{k=1}^{\infty}\1{A_k\leq t} = \max\{k: A_k \leq t\}. \tag{3.1.3} \end{equation*}

The function \(t\to A(t)\) is right-continuous.⁹⁹ Define \(f(x-) = \lim_{y\uparrow x} f(y)\) and \(f(x+) := \lim_{y\downarrow x} f(y)\). A function \(f\) is right-continuous if \(f(x) = f(x+)\) for all \(x\), and left-continuous if \(f(x) = f(x-)\) for all \(x\).

Conversely, we can retrieve the arrival times \(\{A_{k}\}\) from the process \(\{A(t)\}\). For example, if \(A(s) = k-1\) and \(A(t) = k\), then the arrival time \(A_k\) of the \(k\)th job must lie somewhere in \((s,t]\). Specifically,

\begin{equation*} A_k = \min\{t: A(t) \geq k\}, \quad A(0) = 0. \end{equation*}

And if \(\{A_{k}\}\) is known, then \(X_k = A_k - A_{k-1}\). Thus, the arrival times, the inter-arrival times, and the number of arrivals as a function of time can be derived from one another.

Similarly, for the number of jobs that started service, \(\As(t)\), and the departures \(D(t)\) we have

\begin{align*} \As(t) &= \max\{k : \As_k \leq t\}, & \As_{k} &=\min\{t: \As(t) \geq k\}, \\ D(t) &= \max\{k : D_k \leq t\}, & D_{k} &=\min\{t: D(t) \geq k\}. \end{align*}

The work in the system \(W(t)\) is the time needed to clear the remaining service of all jobs present at time \(t\). To construct \(\{W(t)\}\) graphically, we draw lines that start at points \((A_k, \J_k)\) and go down with slope -1 until the lines hit the \(x\)-axis.¹⁰¹⁰ The slope is \(-1\) because a single server processes one hour of work per hour. If the system is empty, the work in the system remains at zero until an arrival occurs. To express \(W(t)\) as a function of time, note that \(A(t)\) is the last arrival at \(t\), thus, \(D_{A(t)}\) is the time this job leaves. Therefore,

\begin{equation*} W(t)= [D_{A(t)} - t]^+. \tag{3.1.4} \end{equation*}

The work in the system is equal to the so-called virtual waiting time. For example, if \(W(t) = 5\) hours¹¹¹¹ It takes 5 hours to complete all jobs in the system at time \(t\). , then a job arriving at \(t\) would have a waiting time of \(5\) hours. Fig. 1 shows how the above concepts relate to each other.

Once we have the arrival and departure processes, it is easy to compute the number of jobs in the system at time \(t\), as illustrated in Fig. 2,

\begin{equation*} \L(t) = A(t) - D(t) + \L(0), \tag{3.1.5} \end{equation*}

where \(\L(0)\) is the number of jobs in the system at time \(t=0\); typically \(\L(0)=0\). Recalling that a job can be either waiting or in service, we distinguish \(\L(t)\) (in the system), \(\QQ(t)\) (in queue), and \(\Ls(t)\) (in service):

\begin{align*} \L(t) &= \QQ(t) + \Ls(t), & \Ls(t) &= \As(t) - D(t) = \L(t) - \QQ(t). \end{align*}

Note that statistics obtained from a sample path \(t \to \QQ(t)\) can be quite different from the samples \(k\to \QQ(A_k-)\) as observed by arriving jobs, and likewise for \(L(t)\) and \(L(A_{k}-)\).¹²¹² We need to be careful about the left and right limits at jump epochs.

The following indicator function

\begin{equation*} I_k(s) := \1{A_{k}\leq s < D_{k}}, \tag{3.1.6} \end{equation*}

represents the presence of job \(k\) at time \(s\).¹³¹³ \(I_{k}(s)=1\) if job \(k\) is present in the system at time \(s\) and zero otherwise. Clearly, \(\int_{0}^{t} I_{k}(s) \d s\) corresponds to the total time job \(k\) spends in the system up to time \(t\) while summing over all jobs present at time \(s\) gives \(L(s)\). Therefore,

\begin{align*} \J_k &= \int_0^\infty I_k(s)\d s, & \L(t) &= \sum_{k=1}^\infty I_k(t). \tag{3.1.7} \end{align*}

With this we can relate sojourn times and the number of jobs in the system. When \(n\) is the index of the last departure of a simulation, we have on the one hand,¹⁴¹⁴ At time \(D_n\) all \(n\) jobs have departed.

\begin{align*} \frac{1}{D_{n}} \sum_{k=1}^n \int_{0}^{D_{n}}I_{k}(s) \d s = \frac{1}{D_{n}} \sum_{k=1}^n \J_{k}, \end{align*}

while on the other hand,

\begin{align*} \frac{1}{D_{n}}\int_{0}^{D_{n}} \sum_{k=1}^n I_{k}(s) \d s = \frac{1}{D_{n}}\int_{0}^{D_{n}} \L(s) \d s = \bar L. \end{align*}

As the expressions at the LHSs are equal, the sum over all sojourn times divided by the time of the last departure equals the time-average number of jobs in the system:

\begin{align*} \bar L \approx \frac{1}{D_{n}} \sum_{k=1}^n \J_{k}. \tag{3.1.8} \end{align*}

We finish this section with constructing the recursions for the waiting times in a multi-server queue. Jobs line up in queue \(i\), \(i=1, \ldots, c\) that is served by a server working with rate \(\mu(i)\). When job \(k\) arrives, it sees a vector \(w_k = (w_{k}(1), \ldots, w_{k}(c))\) of waiting times at the queues. Suppose the job selects the line with the shortest waiting time.¹⁵¹⁵ This is not necessarily the same as the shortest queue. Let \(e(i)\) denote the \(i\)th unit vector¹⁶¹⁶ A vector with a 1 at place \(i\) and zeros elsewhere. , and \(\mathbf{1} = \sum_{i=1}^{c} e_{i}\). Then, in analogy with Eq. (3.1.2), the waiting time \(w_{k+1}\) updates as¹⁷¹⁷ \([\cdot]^+\) applies element-wise to the vector.

\begin{align*} s_{k} &= \argmin_{i \in \{1, \ldots, c\}}\{w_{k}(i)\},& w_{k+1} &= [w_{k} + S_{k} e(s_{k})/\mu(s_{k}) - X_{k}\mathbf{1}]^+, \tag{3.1.9} \end{align*}

where \(s_k\) is the queue selected by job \(k\).

While the above recursions are charming in their simplicity, it is hard to generalize these recursions to more realistic queueing situations involving control. In later sections we discuss how to use discrete-event simulation to analyze such harder cases.