Lab: Queuing Theory with R Simulation

Overview

In this lab, you will use R to simulate simple queuing systems and connect your results to basic queuing theory concepts: arrivals, service, waiting time, and server capacity.

We will work through three main parts:

1. Simulating Poisson arrivals (interarrival times).
2. Simulating an M/M/1 queue (single server).
3. Simulating an M/M/c queue (multiple servers).

You can run this lab in RStudio by knitting the R Markdown file.

---

Part 1: Poisson Arrivals

We assume customers arrive according to a Poisson process with rate \( \). The interarrival times follow an Exponential distribution with rate \( \).

lambda <- 4      # arrivals per hour
n <- 50          # number of arrivals to simulate

interarrival <- rexp(n, rate = lambda)
arrival_times <- cumsum(interarrival)

head(interarrival)

## [1] 0.18879546 0.29541069 0.03642668 0.03494882 0.10901716 0.72374213

head(arrival_times)

## [1] 0.1887955 0.4842062 0.5206328 0.5555817 0.6645988 1.3883409

plot(arrival_times, type = "h", lwd = 3,
     xlab = "Customer index", ylab = "Arrival time",
     main = "Simulated arrival times (Poisson process)")

Try changing lambda and re-running the chunk. How does the spacing of arrivals change?

---

Part 2: Simulating an M/M/1 Queue

We now add a single server with Exponential service times with rate \( \).

set.seed(123)

lambda <- 3   # arrival rate
mu <- 6       # service rate
n <- 1000      # number of customers

interarrival <- rexp(n, lambda)
service <- rexp(n, mu)
arrival <- cumsum(interarrival)

start_service <- numeric(n)
departure <- numeric(n)
wait <- numeric(n)

start_service[1] <- arrival[1]
departure[1] <- start_service[1] + service[1]

for (i in 2:n) {
  start_service[i] <- max(arrival[i], departure[i-1])
  wait[i] <- start_service[i] - arrival[i]
  departure[i] <- start_service[i] + service[i]
}

mean_waitall <- mean(wait)
mean_waitall

## [1] 0.1119213

# second half of the vector
mean_wait_2nd <- mean(wait[(n/2):n])
mean_wait_2nd

## [1] 0.1198462

plot(wait, type = "h", lwd = 3,
     main = "Waiting times in an M/M/1 queue",
     xlab = "Customer index", ylab = "Waiting time")

Compute the theoretical average waiting time in queue for an M/M/1 system:

\[ W_q = . \]

rho <- lambda / mu
Wq_theory <- lambda / (mu * (mu - lambda))
rho

## [1] 0.5

Wq_theory

## [1] 0.1666667

Compare mean_wait from simulation with Wq_theory. Are they close?

Try changing lambda and mu so that utilization rho gets closer to 1. What happens to the waiting times?

---

Part 3: Simulating an M/M/c Queue

For a multi-server system, we can use the queuecomputer package to simplify the logic.

# install.packages("queuecomputer") # if not installed
library(queuecomputer)

set.seed(123)
lambda <- 10    # arrival rate
mu <- 5         # service rate
c <- 3          # number of servers
n <- 20        # number of customers

arrivals <- cumsum(rexp(n, lambda))
service <- rexp(n, mu)

results <- queue_step(arrivals, service, servers = c)

summary(results)

## Total customers:
##  20
## Missed customers:
##  0
## Mean waiting time:
##  0.146
## Mean response time:
##  0.369
## Utilization factor:
##  0.80316720291063
## Mean queue length:
##  1.81
## Mean number of customers in system:
##  3.99

plot(results)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

Try modifying:

• c (number of servers),
• lambda (arrivals), and
• mu (service rate)

and examine how the summary statistics and plot change.

---

Reflection

In a few sentences, reflect on:

• How simulation helped you see queuing behavior (waiting, congestion) that is hard to understand from formulas alone.
• One real-world system (airline, hospital, store, website, etc.) where you think a simple queuing model could improve performance.

arrivals <- rpois(10000, lambda = 4)
hist(arrivals, col="skyblue", main="Poisson(4) Simulation")

dpois(1,4)

## [1] 0.07326256

dpois(3,4)

## [1] 0.1953668

dpois(4,4)

## [1] 0.1953668

exp distribution

dexp(0.5, rate = 4)

## [1] 0.5413411

pexp(0.5, rate = 4)

## [1] 0.8646647