Computationally Efficient Simulation of Queues: The R Package queuecomputer

Large networks of queueing systems model important real-world systems such as MapReduce clusters, web-servers, hospitals, call centers and airport passenger terminals. To model such systems accurately, we must infer queueing parameters from data. Unfortunately, for many queueing networks there is no clear way to proceed with parameter inference from data. Approximate Bayesian computation could offer a straightforward way to infer parameters for such networks if we could simulate data quickly enough. We present a computationally efficient method for simulating from a very general set of queueing networks with the R package queuecomputer. Remarkable speedups of more than 2 orders of magnitude are observed relative to the popular DES packages simmer and simpy. We replicate output from these packages to validate the package. The package is modular and integrates well with the popular R package dplyr. Complex queueing networks with tandem, parallel and fork/join topologies can easily be built with these two packages together. We show how to use this package with two examples: a call center and an airport terminal.


Introduction
The queues we encounter in our everyday experience, where customers wait in line to be served by a server, are a useful analogy for many other processes. We say analogy because the word customers could represent: MapReduce jobs (Lin, Zhang, Wierman, and Tan 2013); patients in a hospital (Takagi, Kanai, and Misue 2017); items in a manufacturing system (Dallery and Gershwin 1992); calls to a call center (Gans, Koole, and Mandelbaum 2003); shipping containers in a seaport (Kozan 1997) or even cognitive tasks (Cao 2013). Similarly, server could represent: a compute cluster; medical staff; machinery or a customer service representative at a call center. Queueing systems can also be networked together to form queueing networks. We can use queueing networks to build models of processes such as provision of internet services (Sutton and Jordan 2011), passenger facilitation at international airports (Wu and Mengersen 2013) and emergency evacuations (Van Woensel and Vandaele 2007). Clearly queueing systems and queueing networks are useful for understanding important real-world systems.
Performance measures for a given queueing system can often only be derived through simulation. Queues are usually simulated with discrete event simulation (DES; Rios Insua, Ruggeri, and Wiper 2012, p. 226). In DES changes in state are discontinuous. The state is acted upon by a countable list of events at certain times which cause the discontinuities. If the occurrence of an event is independent of everything except simulation time it is determined; otherwise, it is contingent (Nance 1981).
Popular DES software packages are available in many programming languages including: the R (R Core Team 2020) package simmer (Ucar, Smeets, and Azcorra 2018), the Python (Rossum et al. 2011) package simpy (Lünsdorf and Scherfke 2013) and the Java (Gosling 2000) package JMT (Bertoli, Casale, and Serazzi 2009). DES packages are often so expressive that they can be considered languages in their own right, indeed the programming language SIMULA (Dahl and Nygaard 1966) is a literal example of this.
The R package queuecomputer (Ebert 2020), which is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=queuecomputer, implements an algorithm that can easily be applied to a wide range of queueing systems and networks of queueing systems. It is vastly more computationally efficient than existing approaches to DES. We term this new computationally efficient algorithm queue departure computation (QDC). Computational efficiency is important because if we can simulate from queues quickly, then we can embed a queue simulation within an approximate Bayesian computation (ABC) algorithm (Sunnåker, Busetto, Numminen, Corander, Foll, and Dessimoz 2013) and estimate queue parameters for very complicated queueing models in a straightforward manner.
In Section 2 we review the literature on queueing theory and develop notation used throughout this paper. In Section 3 we present the QDC algorithm and compare it to DES. We demonstrate usage of the package in Section 4. Details of implementation and usage are discussed in Section 5. The package is validated in Section 6 by replicating results from DES packages simpy and simmer. We compare computed performance measures from the output of a queuecomputer simulation to theoretical results for M/M/2 queueing systems. We benchmark the package in Section 7 and compare computation time with simpy and simmer. Examples in Section 8 are used to demonstrate how the package can be used to simulate a call center and an international airport terminal.

Queueing theory
Queueing theory is the study of queueing systems and originated from the work of Agner Krarup Erlang in 1909 to plan infrastructure requirements for the Danish telephone system (Thomopoulos 2012, p. 2). An overview of the notation and definitions used in this paper is given in Table 1 This vector represents the time at which each server 1 : K will next be free. We consider this vector to be the state of the system. A queueing system is defined as follows. Each customer i = 1, 2, . . . has an arrival time a i (or equivalently an inter-arrival time δ i = a i − a i−1 , a 0 = 0) and an amount of time they require with a server, called the service time s i . Typically a server can serve only one customer at a time. A server which is currently serving another customer is said to be unavailable, a server without a customer is available. If all servers are unavailable when a customer arrives then customers must wait in the queue until a server is available. Detailed introductions to queueing systems can be found in standard texts such as Bhat (2015).
The characteristics of a queueing system are expressed with the notation of Kendall (1953). This notation has since been extended to six characteristics: • f δ , inter-arrival distribution; • f s , service distribution; • K, number of servers ∈ N; • C, capacity of system ∈ N; • n, customer population ∈ N; and • R, service discipline.
Choices for inter-arrival and service distributions are denoted by "M" for exponential and independently distributed, "GI" for general and independently distributed and "G" for general without the independence assumption. The capacity of the system C refers to the maximum number of customers within the system at any one time 1 . Customers are within the system if they are being served or waiting in the queue. The customer population n is the total number of customers including those outside of the system (yet to arrive or already departed). The service discipline R defines how customers in the queue are allocated to available servers. The most common service discipline is first come first serve (FCFS). To specify a queueing system, these characteristics are placed in the order given above and separated by a forward slash "/". The simplest queueing system is exponential in distribution for both the inter-arrival δ i iid ∼ exp(λ) ∀i ∈ 1 : n and service processes s i iid ∼ exp(µ) ∀i ∈ 1 : n, where λ and µ are exponential rate parameters. Additionally, K is set to 1, C and n are infinite, and R is FCFS. It is denoted by M/M/1/∞/∞/FCFS, which is shortened to M/M/1.
Parameter inference for this system was considered first by Clarke (1957), estimators were derived from the likelihood function. This likelihood is later used by Muddapur (1972) to derive the joint posterior distribution. Bayesian inference for queueing systems is summarized in detail by Rios Insua et al. (2012).
Managers and planners are less interested in parameter inference and more interested in performance measures such as: N (t), the number of customers in system at time t;B, the average number of busy servers; ρ, the resource utilization; andw, the average waiting time for customers. If λ < Kµ the queueing system will eventually reach equilibrium and distributions of performance measures become independent of time.
In the case of a M/M/K system, these performance measures have equilibrium distributions which can be derived analytically. These derivations are found in standard queueing theory textbooks (Lipsky 2008;Thomopoulos 2012). For instance, the limit probability of N customers in the system P(N ) is where ρ, the resource utilization, is defined as λ Kµ . For an M/M/K system this is equal to the expected number of busy servers divided by the total number of servers E(B) K (Cassandras and Lafortune 2009, p. 451). The expected number of customers in the system is (Bhat 2015) and the expected waiting time is If the parameters of f δ and f s are uncertain, then we must turn to predictive distributions for estimates of performance measures, which are computed analytically for M/M/K queues (Equations 2 and 3). Predictive distributions of performance measures using Bayesian posterior distributions are derived by Armero (1994); Armero and Bayarri (1999). Jackson (1957) was one of the first to consider networks of queueing systems. In a Jackson network, there is a set of J queueing systems. After a customer is served by queueing system j, they arrive at another queueing system with fixed probability p j,k . Customers leave the system with probability 1 − J k=1 p j,k . Other examples of queueing networks include the tandem (Glynn and Whitt 1991), parallel (Hunt and Foote 1995) and the fork/join (Kim and Agrawala 1989) topologies.
In a tandem queueing network, customers traverse an ordered series of queues before departing the system. Real examples of such systems include airport terminals, internet services and manufacturing systems. In a parallel network, the set of customers is partitioned into different queueing systems, with their own two vectors a and s. In a fork/join network each task (another term for customer) is split into a number of subtasks which are to be completed by distinct parallel servers. The task is only complete once all its associated subtasks arrive at the join point. One example of this situation is that of arriving passengers and their bags at an airport.
Most models of queueing systems assume time-invariant inter-arrival and service processes. In practice, many real-world queues have inter-arrival processes which are strongly timedependent, such as: call centers (Weinberg, Brown, and Stroud 2007;Brown, Gans, Mandelbaum, Sakov, Shen, Zeltyn, and Zhao 2005), airport runways (Koopman 1972) and hospitals (Brahimi and Worthington 1991). In the case of the M/M/1 queue, we can adapt the notation to M (t)/M (t)/1 to represent exponential processes where parameters λ(t) and µ(t) change with time. Such queueing systems are referred to as dynamic queueing systems.
In general, analytic solutions do not exist for dynamic queueing systems (Malone 1995;Worthington 2009). Green, Kolesar, and Svoronos (1991) showed that using stationary queueing systems to model dynamic queueing systems leads to serious error even if deviation from stationarity is slight. The problem is compounded once we consider queueing networks. Understanding long-term and transient behavior of such queues can only be achieved with approximation methods or simulation. We now detail the QDC algorithm, a computationally efficient method for simulating queueing systems.

Fixed number of servers
QDC can be considered as a multiserver extension to an algorithm presented by Lindley (1952). For a single server queueing system, the departure time of the ith customer is: since the customer either waits for a server or the server waits for a customer. The algorithm (not the paper) was, surprisingly, not extended to multiserver systems until Krivulin (1994). However with each new customer i the algorithm must search a growing i + 1 length vector. This algorithm, therefore, scales poorly, with computational complexity O(n 2 ), where n is the number of customers. Kin and Chan (2010) adapted the original algorithm of Kiefer and Wolfowitz (1955) to an O(nK) algorithm for multiserver tandem queues with blocking, that is G/G/K/C queueing systems where C is the maximum capacity number of customers in the queueing systems.
QDC can also be viewed as a computationally efficient solution to the set of equations presented in Sutton and Jordan (2011, p. 259) for FCFS queueing systems. There is a single queue served by a fixed number of K servers. The ith customer observes a set of times b i = {b ik |k ∈ 1 : K} which represents the times when each server will next be available. The customer i selects the earliest available server p i = argmin(b i ) from b i . The departure time for the ith customer is, therefore, d i = max(a i , b p i ) + s i , since the server must wait for the customer or the customer must wait for the server. The QDC algorithm for a fixed number of servers (Algorithm 1) pre-sorts the arrival times. Rather than assigning a b i for each customer i to form the matrix b ∈ M n×K , QDC considers b as a continually updated K length vector representing the state of the system. This algorithm is simple and computationally efficient. At each iteration of the loop, we need only search b, a K length vector for the minimum element in code line 8. In the language of DES, we consider b as the system state and a as the event list, which are all determined events. This differs from conventional DES approaches to modeling queueing systems where the queue length is the system state, and both a and d constitute the event list, where the events of a are determined and the events of d are continually updated and therefore contingent.
Algorithm 1 can simulate any queue of the form G(t)/G(t)/K/∞/n/FCFS where K and n can be made arbitrarily large. Furthermore, the inter-arrival and service distributions can be of completely general form and even have a dependency structure between them. Since the arrival and service times are supplied by the user rather than sampled in-situ, the algorithm "decouples" statistical sampling from queue computation. This frees the user to simulate queues of arbitrarily complex f δ,s , where K is fixed.
Sort (a, s) in terms of a (ascending).

3:
Create vector p ∈ N n .
for i ∈ 1 : n do 8: Put (a, d, p) back to original (input) ordering of a.

Conditional case
Suppose that the number of servers that customers can use changes throughout the day. This reflects realistic situations where more servers are rostered on for busier times of the day. We say that for a certain time t, the customers have a choice of K(t) open servers from K. This means that there are K(t) servers rostered-on for time t. We define the term closed as the opposite of open.
We represent the number of open servers throughout the day as a step function. Time is on the positive real number line and is partitioned by L knot locations The number of open servers in each epoch is represented by a L + 1 length vector y = (y 1 , . . . , y L+1 ) ∈ N L+1 0 . If we assume that none of the service times s span the length of more than one epoch ( then we need to consider a change in state over at most 1 knot location. This step function is determined input by the user. Like the arrival and service times (a, s), it is changeable by the user before the simulation but not during the simulation.
We close server k by writing an ∞ symbol to b k ensuring that no customer can use that server. If the server needs to be open again at time t, we write t to b k allowing customers to use that server. Since x now corresponds to changes in b, it is part of the event list along with a. The entire event list is still determined and need not be updated mid-simulation.
This algorithm can simulate queues of form Sort (a, s) in terms of a (ascending). 3: Create vector p ∈ N n .

11:
for i ∈ 1 : n do 12: for k ∈ 1 : K do 13: where Condition 4 does not hold, we must use the less computationally efficient but more general unconditional algorithm below.

Unconditional case
If Condition 4 does not hold or if, otherwise, we wish to control exactly which servers are open at what time then we must use a less computationally efficient algorithm (Algorithm 4). Each server k has its own partition of L k knot locations x k = (x k,1 , . . . , x k,L k ) ∈ R L k + and each y k = (y k,1 , . . . , y k,L k+1 ) is an alternating sequence of 0 and 1s of length L k+1 indicating whether the server is open or closed respectively for the associated epoch. The vector c is used slightly differently to how it is used in Sutton and Jordan (2011). We use it to represent the time at which each server is next available for the current customer i, given the current system state b. It is the output of the next_fun function.

This algorithm can simulate queueing systems of form G(t)/G(t)/K(t)/∞/n/FCFS, where K(t) refers to the number of open servers changing with time.
In addition, we can specify which particular servers are available when, not just how many and we are not bound by Condition 4. Once again we note that b can be considered as the system state and the event list is formed by a and x 1 , . . . , x K , where K in this context is the total number of distinct servers. This function can be called with the queue_step function in queuecomputer by supplying a 'server.list' object to the servers argument. For the rest of this paper we focus on Algorithms 1 and 2 for their relative conceptual simplicity and computational efficiency.

Discussion
With the algorithms so far presented, we can simulate from a very general set of queueing systems G(t)/G(t)/K(t)/∞/n/FCFS in a computationally efficient manner. In contrast to the algorithm of Kin and Chan (2010), the state vector b is written over in each iteration. The memory usage for QDC, therefore, scales with O(n) rather than O(nK).
Tandem queueing networks can be simulated by using the output of one queueing system as the input to the next queueing system. We demonstrate this idea with the international airport terminal example in Section 8.2. Fork/join queueing networks are addressed in the next section where we explain the implementation details of queuecomputer with regards to the QDC algorithm.

Usage
The purpose of the package queuecomputer is to compute, deterministically, the output of a queueing system given the arrival and service times for all customers. The most important function is queue_step. The first argument to queue_step is a vector of arrival times, the second argument is a vector of service times and the third argument specifies the servers available.
R> library("queuecomputer") R> arrivals <-cumsum(rexp (100) The output of a queue_step function is a 'queue_list' object. We built a summary method for objects of class 'queue_list', which we now demonstrate. If the last element of y is zero, it is possible that some customers will never be served, this is the "Missed customers" output. The performance measures that follow are the mean waiting timew, the mean response timer = d − a, the observed utilization factorB/K, the mean queue length and the mean number of customers in the system respectively. The utilization factorB/K takes into account the changing number of open servers K(t) where Algorithm 2 is used. We now explain the implementation details of the package.

Implementation
The for loops within Algorithms 1 and 2 are written in C++ with the Armadillo library (Sanderson and Curtin 2016). The C++ for loops are called using the R packages Rcpp (Eddelbuettel and François 2011) and RcppArmadillo (Eddelbuettel and Sanderson 2014). We use R to provide wrapper functions for the C++ code.
The queue_step calls the more primitive queue function which is a wrapper for S3 methods which implement Algorithms 1, 2 or 4 depending on the class of the object supplied to the server argument of queue_step. If server inherits from numeric, then queue runs Algorithm 1, if it is a 'server.stepfun' object then queue runs Algorithm 2, and if it is a 'server.list' object then queue runs Algorithm 4. The queue function computes departure times d and server allocations p and the queue_step function adds additional output such as waiting times and queue lengths which are used in summary and plot methods.
To simulate fork/join networks, the queuecomputer function wait_step provides a simple wrapper to the base function pmax.int, this function computes the maximum of each row for a set of two equal length numeric vectors. The vectors represent the departure times for each subjob and the departure time for the entire job is the maximum of each subjob.
In simmer and simpy users supply generator functions for simulating δ and service times s, the user enters the set of input parameters θ I for these generator functions and starts the simulation. The inter-arrival time is re-sampled after each arrival and the service time is sampled when the server begins with a new customer. This makes it difficult to model queues where distributions for inter-arrival times do not make sense: like the immigration counter for an airport, where multiple flights generate customers; or when arrival times and service times are not independent. In queuecomputer sampling is "decoupled" from computation, the user samples a and service times s using any method. The outputs d and p are then computed deterministically.
We now demonstrate the validity of queuecomputer's output by replicating results from the DES packages simmer and simpy. We then replicate equilibrium analytic results of performance measures for the M/M/2 queue.

Comparison with simmer and simpy
To demonstrate the validity of the algorithm we consider a M/M/2/∞/1000/FCFS queue. If QDC is valid for any M/M/K queueing system, then it is valid for any G(t)/G(t)/K queueing system. This is because any non-zero (a, s) could conceivably come from two exponential distributions, even if the probability of the particular realization is vanishingly small. We replicate exact departure times computed with the simmer and simpy packages using queuecomputer. First, we generate a and s to be used as input to all three packages. R> set.seed(1) R> n_customers <-10^4 R> lambda_a <-1/1 R> lambda_s <-1/0.9 R> interarrivals <-rexp(n_customers, lambda_a) R> arrivals <-cumsum(interarrivals) R> service <-rexp(n_customers, lambda_s) We now input these objects into the three scripts using queuecomputer, simmer, or simpy. First, we run the queuecomputer script. The queuecomputer_output object is sorted in ascending order so that the departure times can be compared to the DES packages.
R> queuecomputer_output <-queue_step(arrivals = arrivals, + service = service, servers = 2) R> head(sort(depart(queuecomputer_output))) [1] 1.340151 2.288112 2.639976 2.796572 3.249794 5.714967 The DES packages simmer and simpy are not built to allow users to input (a, s) directly. Rather, the user supplies parameters for f δ and f s so that inter-arrival and service times can be sampled at each step when needed. To allow simmer and simpy to accept pre-sampled input (a, s) we use generator functions instead of rexp(rate) or random.expovariate(rate) calls in R and Python respectively, details of this work can be found in the supplementary material. We create an interface to simmer so that it can be called in the same way as queuecomputer.
R> simmer_output <-simmer_step(arrivals = arrivals, + service = service, servers = 2) R> head(simmer_output) [1] 1.340151 2.288112 2.639976 2.796572 3.249794 5.714967 The same departure times are observed. Similarly in Python we create an interface to simpy so that it can be called in a similar way to queuecomputer.

])
A check of all three sorted vectors of d from each package revealed that all were equal to within 5 significant figures for every d i , i = 1 : 1000.

Replicate theoretical results for M/M/3
We use a M/M/3/∞/5 × 10 6 /FCFS simulation in queuecomputer to replicate theoretical equilibrium results for key performance indicators for a M/M/2/∞/∞/FCFS queueing system. We set λ to 1 and set µ to 2.

Theoretical results
We first note that the traffic intensity is ρ of 2/3 = 0.6, which should correspond to the average number of busy servers. The probability of N customers in the system is given by Equation 1. We perform this computation up to N = 20 and display the results in Figure 1. The expected waiting time E(w) is 0.4 and the expected number of customers in the system E(N ) is 2.8.

Simulation results
The inputs a and s must first be generated.
R> MM3 <-queue_step(arrivals = arrivals, service = service, servers = K) R> summary(MM3) We see that the observed time average number of busy servers is 0.6661402 which is close to 0.6 the value for ρ. We can see that the observed mean waiting time is close to the expected mean waiting time. The expected number of customers in the system, from the distribution P(N ) is close to the observed number of customers in the system. The entire distribution of P(N ) is replicated in Figure 1.

Method
The compare the computational efficiency of each package we compute the departure times from a M/M/2/∞/n/FCFS queueing system, with λ = 1 and µ = 1.1. To understand how n affects computation time we repeat the experiment 100 times for n = 10 2 , 10 3 , 10 5 and 10 6 . We also repeat the experiment at n = 10 7 for queuecomputer. We compare the median time taken for each combination of package and n.
The simulation was conducted on a system with Intel Core i7-6700 CPU @ 3.

Results and discussion
The median computation time for each package and for varying numbers of customers from 100 to 10 6 customers (up to 10 7 customers for queuecomputer) is shown in Figure 2. We observe phenomenal speedups for queuecomputer compared to both packages: compared to simpy speedups of 35 (at 100 customers) to 1000 (at 10 6 customers) are observed, and for simmer speedups of 50 (at 100 customers) to 300 (at 10 6 customers) are observed. The speedup is lower for smaller n since queuecomputer approaches a minimum computation time.
Simulating 10 million customers takes less than 1 second for queuecomputer. We see no reason why queues of different arrival and service distributions should not have similar speedups. This is because, as mentioned earlier, any non-negative (a, s) could come from two exponential distributions.
Clearly, QDC and its implementation queuecomputer are a more computationally efficient way to simulate queueing systems of the form G(t)/G(t)/K/∞/M/FCFS than conventional DES algorithms implemented by simpy and simmer.

Call center
We demonstrate queuecomputer by simulating a call center. The arrival time for each customer is the time that they called, and the service time is how long it takes for their problem to be resolved once they reach an available customer service representative. Let us assume that the customers arrive by a homogeneous Poisson process over the course of the day.
R> service <-rexp(20, 0.5) R> head(service) [1] 2.6669670 1.2434810 0.4197332 0.6188957 2.2118725 1.5483755 We put the arrival and service times into the queue_step function to compute the departure times. Here we have set the number of customer service representatives to two. The servers argument is used for this input.
R> queue_obj <-queue_step(arrivals, service, servers = 2, + labels = customers) R> head(queue_obj$departures_df) We can see that Johatam arrives first but leaves after Beatriz. This is possible because there are two servers. Johatam's service took so long that the next two customers were served by the other server. It is easy to see how the departure times were computed in this simple example. Johatam and Beatriz were the first customers for each server so we can compute their departure time by just adding their service times to their arrival times. Devante, however, had to wait for an available server, since he arrived after the first two customers arrived but before the first two customers departed. He must wait until one of these customers departs before he can be served. We add the departure time of the first customer of server 2 (Beatriz) to his service time to compute his departure time.

R> firstcustomers[2] + service[3]
[1] 3.600039 So the first two customers had no waiting time, but Devante had to wait for an available server. We can compute the waiting times for all three customers in this manner: R> depart(queue_obj)[1:3] -arrivals[1:3] -service[1:3] [1] 0.000000 0.000000 1.097774 The depart function is a convenience function for retrieving the departure times from a 'queue_list' object. The queue_step function returns a 'queue_list' object. There is a summary method for this class within the queuecomputer package, which can be accessed by calling summary(departures). The plot method in queuecomputer for 'queue_list' objects uses the plotting package gg-plot2 (Wickham 2009) to return a list of plots. We produce four plots: a histogram of the arrival and departure times ( Figure 3); a plot of the queue length and number of customers in the system over time ( Figure 4); a plot of the waiting and service times for each customer ( Figure 5); and a plot of the empirical cumulative distribution function for arrival and departure times ( Figure 6). These plots correspond to selections 2, 5 and 6 in the which argument, a similar API to the plot method for 'lm' objects in the stats package (R Core Team 2020).

R> summary(queue_obj)
R> plot(queue_obj, which = c(2, 4, 5, 6)) Notice that in Figure 5, if we draw a horizontal line anywhere on the plot it will never pass through more than one green bar or more than one blue bar. This must be the case otherwise a server would be serving more than one customer at a time.

International airport terminal
The package integrates naturally with the popular data manipulation R package dplyr (Wickham, François, Henry, and Müller 2020). We demonstrate how to integrate queuecomputer and dplyr with a more complex airport terminal example than before (Figure 7). Passengers from a set of 120 flights disembark at the arrivals concourse and proceed through immigration using either the "smart gate" or the "manual gate" route, we therefore have two queues in parallel. The route taken (smart gate or manual gate) by each passenger is predetermined, but the server used by the passenger within these separate queueing systems is not.  Step function plot of customers in queue and system    : Diagram of larger airport scenario, there are 120 flights in total and two multiserver queueing systems operate in parallel. Passengers are pre-assigned to travel through either the "manual" or "smart gate" route through immigration. The passengers and bags are "forked" when each aircraft arrives and then "joined" at the baggage hall.
Their bags are unloaded from the flights and proceed to the baggage hall with a delay, the division of passengers and bags is a fork/join network. The bags and passengers are forked at the arrival concourse and joined at the baggage hall. After immigration, the passengers proceed on to the baggage hall where they pick up their bags.
We have a synthetic dataset of passengers ID from 120 flights FlightNo, with an average of 103.8 passengers per flight for a total of 20,758 passengers. The dataset includes (for each passenger ID): their flight number FlightNo, the arrival time of that flight arrival, the route taken (smart/manual gate) by that passenger route_imm, the arrival times to immigration after they walk through the terminal arrival_imm and the service time needed by the passenger at their immigration queueing system service_imm. Immigration processing is split into two routes with the route_imm variable. The "smart gate" route has 5 servers, whereas the "manual" route has 10 servers before time 600, 12 servers between time 600 and time 780 and 8 servers from time 780 onwards. We store this information in a new data frame called server_df.
The dataset is then processed as if it has been split in two. The column departures_imm represents the times at which passengers depart immigration after having been served either through the manual counter or smart gate. The column departures_bc represents the times that customers leave with their bags from the baggage hall. Waiting times can be summarized with the summarise function from dplyr, here we compute summaries of waiting times for each FlightNo and immigration route route_imm and a summary of waiting times only by route_imm. R> Passenger_df %>% group_by(FlightNo, route_imm) %>% + summarise( + waiting_imm = mean(departures_imm -service_imm -arrive_imm), + waiting_bc = mean(departures_bc -departures_imm)) Boxplot of waiting times Figure 9: Box plot of waiting times for each stage of passenger processing within the international airport terminal.

Conclusion
The R package queuecomputer implements QDC. It can be used to simulate any queueing systems or tandem network of queueing systems of general form G(t)/G(t)/K(t)/∞/n/FCFS. Fast algorithms for multiserver queueing systems have been proposed in the past (Krivulin 1994;Sutton and Jordan 2010;Kin and Chan 2010). These algorithms have generated little notice, even in the cases where their computational efficiency is demonstrated (Kin and Chan 2010). QDC is conceptually simpler, more efficient memory-wise and modular.
We validated QDC with analytic results and by replicating output generated by existing DES packages simpy and simmer. We observe speedups of up to 3 orders of magnitude. The speed of the package will allow queue simulations to be embedded within ABC algorithms, which will be addressed in future work. Unlike existing DES packages, sampling and departure time computation are clearly "decoupled" and therefore allow the user to simulate queueing systems with arrival and service time distributions of arbitrary complexity. The package integrates well with the data manipulation package dplyr and these two packages together allow the user to quickly and easily simulate queueing networks with parallel, tandem and fork/join topologies.