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:

Simulating Poisson arrivals (interarrival times).
Simulating an M/M/1 queue (single server).
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

# ---- FULL PIPELINE IN ONE CHUNK ----

# Load packages
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

# Remove any mistaken object named df (prevents mutate error)
if (exists("df")) rm(df)

## Warning in rm(df): object 'df' not found

# Read your CSV (edit the file path!)
df <- read.csv("/Users/chrisbergersen/Library/CloudStorage/OneDrive-Personal/ANLC-320 Applied Stats/sales.csv", stringsAsFactors = FALSE)

# Add Hours, Date_only, Weekday columns
df <- df %>%
  mutate(
    Hours = as.integer(substr(Time, 1,2)),
    Date_only = mdy(substr(Date, 1, 10)),
    Weekday = wday(Date_only, label = TRUE, abbr = TRUE)
  )

# Preview result
head(df)

##    Invoice.ID Branch      City Customer.type Gender           Product.line
## 1 750-67-8428      A    Yangon        Member Female      Health and beauty
## 2 226-31-3081      C Naypyitaw        Normal Female Electronic accessories
## 3 631-41-3108      A    Yangon        Normal   Male     Home and lifestyle
## 4 123-19-1176      A    Yangon        Member   Male      Health and beauty
## 5 373-73-7910      A    Yangon        Normal   Male      Sports and travel
## 6 699-14-3026      C Naypyitaw        Normal   Male Electronic accessories
##   Unit.price Quantity  Tax.5.    Total      Date  Time     Payment   cogs
## 1      74.69        7 26.1415 548.9715  1/5/2019 13:08     Ewallet 522.83
## 2      15.28        5  3.8200  80.2200  3/8/2019 10:29        Cash  76.40
## 3      46.33        7 16.2155 340.5255  3/3/2019 13:23 Credit card 324.31
## 4      58.22        8 23.2880 489.0480 1/27/2019 20:33     Ewallet 465.76
## 5      86.31        7 30.2085 634.3785  2/8/2019 10:37     Ewallet 604.17
## 6      85.39        7 29.8865 627.6165 3/25/2019 18:30     Ewallet 597.73
##   gross.margin.percentage gross.income Rating Hours  Date_only Weekday
## 1                4.761905      26.1415    9.1    13 2019-01-05     Sat
## 2                4.761905       3.8200    9.6    10 2019-03-08     Fri
## 3                4.761905      16.2155    7.4    13 2019-03-03     Sun
## 4                4.761905      23.2880    8.4    20 2019-01-27     Sun
## 5                4.761905      30.2085    5.3    10 2019-02-08     Fri
## 6                4.761905      29.8865    4.1    18 2019-03-25     Mon

library(dplyr)
library(ggplot2)

# Summarize weekly totals

weekly_sales <- df %>%
group_by(Weekday) %>%
summarise(Total_Sales = sum(gross.income, na.rm = TRUE))

# Horizontal bar chart with labels

ggplot(weekly_sales, aes(x = reorder(Weekday, Total_Sales),
y = Total_Sales,
fill = Total_Sales)) +
geom_col() +
geom_text(aes(label = round(Total_Sales, 0)),
hjust = -0.1, size = 3) +  # labels to the right of bars
coord_flip() +
scale_fill_gradientn(colors = c("red", "orange", "yellow", "green", "blue", "purple")) +
labs(title = "Weekly Total Sales Breakdown",
x = "Week",
y = "Total Sales") +
theme_minimal(base_size = 14) +
theme(legend.position = "none")

library(ggplot2)
library(dplyr)

# Summarize sales by weekday

weekday_sales <- df %>%
group_by(Weekday) %>%
summarise(Total_Sales = sum(gross.income, na.rm = TRUE))

# Horizontal bar chart with labels

ggplot(weekday_sales, aes(x = reorder(Weekday, Total_Sales),
y = Total_Sales,
fill = Total_Sales)) +
geom_col() +
geom_text(aes(label = round(Total_Sales, 0)),
hjust = -0.1, size = 3) +  # labels to the right of bars
coord_flip() +
scale_fill_gradientn(colors = c("pink","orange","yellow","green","blue","purple")) +
labs(title = "Total Sales by Weekday",
x = "Weekday",
y = "Total Sales") +
theme_minimal(base_size = 14) +
theme(legend.position = "none")

library(ggplot2)
library(dplyr)

# Summarize sales by hour

hourly_sales <- df %>%
group_by(Hours) %>%
summarise(Total_Sales = sum(gross.income, na.rm = TRUE))

# Horizontal bar chart with labels

ggplot(hourly_sales, aes(x = factor(Hours),
y = Total_Sales,
fill = Total_Sales)) +
geom_col() +
geom_text(aes(label = round(Total_Sales, 0)),
vjust = -0.3, size = 3) +  # labels above bars
scale_fill_gradientn(colors = c("red","orange","yellow","green","blue","purple")) +
labs(title = "Total Sales by Hour",
x = "Hour of Day",
y = "Total Sales") +
theme_minimal(base_size = 14) +
theme(legend.position = "none")

library(ggplot2)
library(dplyr)

# Summarize sales by weekday and hour

heatmap_sales <- df %>%
group_by(Weekday, Hours) %>%
summarise(Total_Sales = sum(gross.income, na.rm = TRUE), .groups = "drop")

# Ensure weekdays are ordered Mon → Sun

heatmap_sales$Weekday <- factor(heatmap_sales$Weekday,
levels = c("Mon","Tue","Wed","Thu","Fri","Sat","Sun"))

# Heatmap with labels

ggplot(heatmap_sales, aes(x = factor(Hours), y = Weekday, fill = Total_Sales)) +
geom_tile(color = "white") +
geom_text(aes(label = round(Total_Sales, 0)), size = 3) +  # labels in tiles
scale_fill_gradientn(colors = c("yellow","orange","red","purple","blue")) +
labs(title = "Sales Heatmap by Weekday and Hour",
x = "Hour of Day",
y = "Weekday") +
theme_minimal(base_size = 14) +
theme(legend.title = element_text(size = 12),
legend.text = element_text(size = 10))

library(ggplot2)
library(dplyr)

# Summarize sales by weekday and hour

grouped_sales <- df %>%
group_by(Weekday, Hours) %>%
summarise(Total_Sales = sum(gross.income, na.rm = TRUE), .groups = "drop")

# Ensure weekdays are ordered Mon → Sun

grouped_sales$Weekday <- factor(grouped_sales$Weekday,
levels = c("Mon","Tue","Wed","Thu","Fri","Sat","Sun"))

# Grouped bar chart with labels

ggplot(grouped_sales, aes(x = Weekday, y = Total_Sales, fill = factor(Hours))) +
geom_col(position = position_dodge(width = 0.9)) +
geom_text(aes(label = round(Total_Sales, 0)),
position = position_dodge(width = 0.9),
vjust = -0.3, size = 3) +  # labels above bars
scale_fill_manual(values = rainbow(length(unique(grouped_sales$Hours)))) +
labs(title = "Total Sales by Weekday and Hour",
x = "Weekday",
y = "Total Sales",
fill = "Hour of Day") +
theme_minimal(base_size = 14)