System description

In this project we are modelling a queue network on an example of a shoe store. The customers are arriving randomly (i.e., exponential inter-arrival times). There are 5 workers at the store (servers), where 3 of them are helping customers with trying on shoes, and 2 of them are cashiers. Therefore, we are going to have 2 queueus with random exponential service times.

Upon customer’s arrival at the store, they examine the queue’s length, leaving if it’s too long (i.e., balking). If they enter the queue, they either wait until it’s their turn, or they leave if they have been waiting for too long (i.e. reneging). After the customer has seized the worker’s attention, he tries on the shoes, repeating the process if he is not satisfied which is calculated by a simple sample function. If the customer is done trying on the shoes, he releases the worker and either continues to the cashier or leaves based on a probability calculated by his dissatisfaction.

Both queues follow a FIFO/FCFS discipline and they are both finite in queue capacity. The servers differ by their service times and queue capacities. Due to the model’s exponential arrival and service rates, multiple servers and finite queue length’s, we are dealing with a M/M/c/k system.

library(simmer)
library(simmer.plot)
######################################################
# Simulation parameters
######################################################

simulate <- 8 * 60 # Simulation time - one day of work
set.seed(1234)
######################################################
# Parameters of the M/M/c/k queue
######################################################

smu <- 5 # Server service rate
omu <- 10 # Operator service rate
amu <- 5 # Arrival rate
sc <- 2 # Server capacity
mis <- 10 # Max in server
oc <- 3 # Operator capacity
mwo <- 5 # Max with operator
owt <- 10 # Max wait time before reneging from operator
swt <- 20 # Max wait before reneging from server
sqc <- mis - sc # Server queue capacity for balking
oqc <- mwo - oc # Operator queue capacity for balking
######################################################
# Time functions
######################################################

sst <- function() abs(rexp(1, 1/smu)) # Server service time
ost <- function() abs(rexp(1, 1/omu)) # Operator service time
iat <- function() abs(rexp(1, 1/amu)) # Inter-arrival time
omwt <- function() abs(rexp(1, 1/owt)) # Customer max wait time with operator before reneging
smwt <- function() abs(rexp(1, 1/swt)) # Customer max wait time with server before reneging
######################################################
# Utility functions
######################################################

cd <- function() sample(1:33, 1) # Customer dissatisfaction
op_or_se = function(pct) round(ifelse(runif(1, 0, 100) < pct, 1, 2)) # Branching decider (operator or server)

Simmer environment and server trajectory

######################################################
# Trajectories and simulation
######################################################

# Best practice to predefine the simmer environment
store <- simmer()

# Server (cashier) trajectory
server <- trajectory() %>%
  # Renegading from server in 'smwt' = function() abs(rexp(1, 1/swt)) where 'swt' is max wait time before reneging from cashier
  renege_in(smwt,
            # Exit point after reneging (smwt has passed)
            out = trajectory() %>%
              log_(function() {
                paste("Waited", now(store) - get_attribute(store, "start_time"), "minutes. Reneging...")
              })) %>%
  # Either seizeing the cashier or entering the queue if it's formed
  seize("server", 
        continue = FALSE,
        # Decides not to enter the queue because it's too long - Balking
        reject =
          trajectory() %>%
          log_("The server queue is too long. Balking..")) %>%
  # Customer has seized cashier's attention before 'swmt' has expired, therefore, we abort the renegation
  renege_abort() %>%
  # Finish in 'sst' - server service time
  timeout(sst) %>%
  # Release the cashier resource
  release("server")

Operator trajectory

# Operator (worker) trajectory
operator <- trajectory() %>%
  # Renegading from server in 'omwt' = function() abs(rexp(1, 1/owt)) where 'owt' is max wait time before reneging from worker
  renege_in(omwt,
            # Exit point after reneging (omwt has passed)
            out = trajectory() %>%
              log_(function() {
                paste("Waited", now(store) - get_attribute(store, "start_time"), "minutes. Reneging...")
              })) %>%
  # Either seizeing the worker or entering the queue if it's formed
  seize("operator",
        continue = FALSE,
        # Decides not to enter the queue because it's too long - Balking
        reject = 
          trajectory() %>%
          log_("The queue is too long. Balking...")) %>%
  # Customer has seized worker's attention before 'owmt' has expired, therefore, we abort the renegation
  renege_abort() %>%
  # Setting the attribute representing customer's satisfaction after trying on the shoes.
  # We set the mod to '+' if the process is repeated
  set_attribute("customer_dissatisfaction", cd, mod="+") %>%
  # Finish in 'ost' - operator service time
  timeout(ost) %>%
  # Go back 2 steps - set_attribute(...) - and check whether customer's dissatisfaction is less than 30
  # If it's less than 30, we go back 2 steps, repeating the process
  # Check function ignores the times parameter (not specified here, therefore, default to infinite - Inf)
  rollback(2, check = function() get_attribute(store, "customer_dissatisfaction") < 30) %>%
  # Customer will leave the store based on a probability calculated by his satisfaction divided by 100 (getting the percentile)
  # Ex. customer_dissatisfaction = 43, prob = 43/100 = .43
  leave(prob = function() get_attribute(store, "customer_dissatisfaction") / 100,
        # Necessary to set 'keep_seized' to false because unlike renege_* functions, leave() function's default behaviour is 'keep_seized = TRUE'
        keep_seized = FALSE) %>%
  # Release the worker resource
  release("operator") %>%
  # Join the server trajectory (cashier)
  # Necessary because of the branching options introduced in the customer trajectory
  join(server)

Customer trajectory

# Customer trajectory
customer <-
  trajectory() %>%
  log_("Customer enters") %>%
  # Utility attribute to calculate customer's time spent in the system
  set_attribute("start_time", function() {now(store)}) %>%
  # Branch represents a crosspath upon customer's arrival at the store
  # 'op_or_se' function represents a probability of choosing a certain path
  # We set it to 90, giving the customer 90% chance that he's in there because he wants to buy shoes,
  # proceeding to the operator (worker) trajectory, and 10% chance that he's in there just to get
  # some random information, proceeding to the server (cashier) trajectory
  branch(option = function() op_or_se(90),
         continue = TRUE,
         join(operator),
         join(server) %>%
           # If the customer choosed the server trajectory,
           # we again set the probability to 10% that he is going to
           # also buy shoes, even though he just came for random information
           # The other (most likely) path is just an exit point after getting random information
           branch(option = function() op_or_se(10),
                  continue = TRUE,
                  join(operator),
                  trajectory() %>%
                    log_("Left after information from server"))) %>%
  # Process finishes succesfully, and customer exit's the store
  log_("Completed")

Environment and simulation

# Environment and simulation
store %>%
  # Adding the server (cashier) resource
  add_resource("server", capacity = sc,queue_size = sqc) %>%
  # Adding the operator (worker) resource
  add_resource("operator", capacity = oc,queue_size = oqc) %>%
  # Generating random customer with 'iat' inter-arrival times
  # 'mon' is set to 2 to enable attribute tracking
  add_generator("customer", customer, iat, mon = 2) %>%
  # We run (simulate) the model for 'simulate' time
  run(simulate)

Customer trajectory plot - possible outcomes

Data analysis

######################################################
# Simulation data
######################################################

# Monitoring store attributes and ordering them by key
attributes_monitor <- store %>%
  get_mon_attributes() %>%
  .[order(.$key),]

# Monitoring store arrivals and transforming wait time by subtracting
# end, start and activitiy time
# Ordering by start times
customer_monitor <- store %>%
  get_mon_arrivals() %>%
  transform(wait = end_time - start_time - activity_time) %>%
  .[order(.$start_time),]

# Monitoring store resource and ordering them by time
resource_monitor <- store %>%
  get_mon_resources() %>%
  .[order(.$time),]

# Variable to be plotted representing store arrivals
customers <- store %>%
  get_mon_arrivals()

# Variable to be plotted representing store resources
servers <- store %>%
  get_mon_resources()

# Variable to be plotted representing store attributes
attributes <- store %>%
  get_mon_attributes()

mean_waiting_time <- mean(customer_monitor$wait)

queue_state <- head(resource_monitor$queue, -1)
server_state <- head(resource_monitor$server, -1)

time_state_lasted <- diff(resource_monitor$time)
time_at_end <- max(resource_monitor$time)

mean_active_customers <- sum(server_state * time_state_lasted) / time_at_end
mean_waiting_customers <- sum(queue_state * time_state_lasted) / time_at_end

# Extracting the average customer dissatisfaction from the attributes variable
keys <- unique(attributes$key)
vals <- sapply(keys, function(x) { sum(attributes[attributes$key==x,]$value) })
out <- data.frame(key=keys, val=vals)
avg_customer_dissatisfaction = out$val[[2]] / length(unique(attributes$name))

cat(" Average waiting = ", mean_waiting_customers, "\n",
    "Average active  = ", mean_active_customers, "\n",
    "Average customer dissatisfaction = ", avg_customer_dissatisfaction, "/ 100")

Plotting data