rush.R

library(ggplot2)
library(parallel)

options(mc.cores = detectCores() - 1)

RNGkind("L'Ecuyer-CMRG")
set.seed(666)

# Function that calculates (by simulation) the overflow probability of the
# system
prob <- function(
  n = 80,    # Number of employees likely to clock out during the rush period
  p = 300,   # Duration of rush period (seconds)
  d = 5,     # Time required for clocking out (seconds)
  k = 16,    # Number of time clocks
  s = 10^5   # Number of simulations
) {

  # The system can be overloaded only if the number of employees is stricly
  # greater than the number of devices
  if (n <= k) stop("The system cannot be overloaded")

  # Each employee arrives randomly in the rush period according to a uniform
  # distribution. The clocking out times are stored in a matrix.
  t <- matrix(runif(s * n) * p, nrow = s, ncol = n)

  # Sort each line of the clocking out time matrix (sorting by time for each
  # simulation)
  t <- t(apply(t, 1, sort))

  # For each employee sorted in the order of arrival, we calculate the time
  # difference with the employee with k more ranks in the order of arrival. If
  # one of these time differences is less than the time needed to clocking out,
  # then the system is overflowed.
  o <- sapply(1:(n - k), function(j) t[, j + k] - t[, j] <= d)

  # Proportion of overflowding
  mean(apply(o, 1, any))
}

# Creation of a table containing the probabilities according to different
# numbers of devices. The other parameters are fixed. Default values are
# defined in the function.
prob_tbl <- do.call(rbind, mclapply(1:16, function(k) {
  data.frame(number_of_devices = k, probability = prob(k = k))
}))

# Figure
prob_fig <- ggplot(prob_tbl, aes(x = number_of_devices, y = probability)) +
  geom_point() +
  geom_line() +
  labs(x = "Number of devices", y = "Probability of overflowding")

# Export figure in working directory
png("rush.png")
print(prob_fig)
dev.off()