Introduction:

Simulation modeling solves real-world problems safely and efficiently. It provides an important method of analysis which is easily verified, communicated, and understood. Across industries and disciplines, simulation modeling provides valuable solutions by giving clear insights into complex systems.

Abstract

A discrete-event simulation (DES) models a system as a discrete sequence of events.The execution of an event at a given time will lead to a change in the system’s state. One can configure the combination of events in a system in such a way that it can simulate an actual process. The applications of DES-based simulation are broad. If used to simulate a process one can for example gain insights into the process’ risk, efficiency and effectiveness. Also, by simulation of an alternative configuration of a process, one can proactively estimate the benefits of changes to the process. This in turn allows one to get clear insights into the benefits of process redesign strategies (e.g., extra resources). simmer is a process-oriented and trajectory-based discrete-event simulation package for R. The architecture encloses a robust and fast simulation core written in C++ with automatic monitoring capabilities. It provides a rich and exible R API that revolves around the concept of trajectory, a common path in the simulation model for entities of the same type. It is Designed to be a generic framework, it leverages the power of Rcpp to boost performance and make DES modelling in R not only effective but also efficient. As a noteworthy characteristic, simmer exploits the concept of trajectory: a common path in the simulation model for entities of the same type. As a modelling framework it is flexible and simple to use, and leverages the chaining/piping workflow introduced by the ‘magrittr’ package. The development of the’simmer’ package started in the second half of 2014. The initial need for a DES framework for R came up in projects related to process optimization in healthcare facilities. Most of these cases involved patients following a clear trajectory through a care process. This background is not unimportant, as it led to the adoption and implementation of a trajectory concept at the very core of simmer’s DES engine. This strong focus on clearly defined trajectories is somewhat innovative and, more importantly, very intuitive. Over time, the simmer package has seen significant improvements and has been at the forefront of DES for R. Although it is the most generic DES framework, it is however not the only R package which delivers such functionality. For example, the ’SpaDES’ package (Chubaty and McIntire 2016) focuses on spatially explicit discrete models, and the queue computer package implements an efficient method for simulating queues with arbitrary arrival and service time.

People enter the registration desk, and there is a registrar for the purpose, and they entered the facility first. They pick up the test kit from the test staff and then get into the test bays. There are four test bays. They occupy one of the test bays. Take the swab and return to the test staff and return the swab and then go back to the same bay there and clean the squat, clean the bay and then leave. So, until they clean the bay, they actually occupy pretty much that test bay. Once this test staff receives the return swab, they perform the test and the result is sent to them. So that’s the process. These other process steps and these other resources which you struck test bays for and test off one.

Objective:

The main concern of the study is to perform the ‘Discrete Event simulation’ to model a system as a discrete sequence of events. The execution of an event at a given time will lead to a change in the system’s state. One can configure the combination of events in a system in such a way that it can simulate an actual process.The study allows us to understand how many people can go through the process and complete and come out of the end in a minute. There is a measure call ‘Process Cycle efficiency’, sometimes referred to as the “value added ratio” is a measurement of the amount of the value-added time is a process will figure out the time spent inside by each person spend in the actual activity and the remaining time is waiting time.

Methodology:

The discrete event simulation (DES) of covid-19 lateral test process has performed using R programming language for statistical computing and graphics supported by the R core team and the R Foundation for statistical computing. The simulation has performed using “simmer” package. The simmer package brings discrete-event simulation to R, it is a process-oriented and trajectory-based discrete-event simulation package for R. It is designed as a generic yet powerful process-oriented framework. The architecture encloses a robust and fast simulation core written in C++ with automatic monitoring capabilities. It leverages the power of ‘Rcpp’ (The ‘Rcpp’ package provides R functions as well as C++ classes which offer a seamless integration of R and C++. Many R data types and objects can be mapped back and forth to C++ equivalents which facilitates both writing of new code as well as easier integration of third-party libraries.) to boost performance and make DES modelling in R not only effective but also efficient.

R Simmer Packages & Functions

# R Simmer Packages & Functions

pkgs = c("simmer","simmer.plot", "simmer.bricks","data.table") 
check = pkgs %in% installed.packages()
if(!check[1]) install.packages("simmer")
if(!check[2]) install.packages("simmer.plot")
if(!check[3]) install.packages("simmer.bricks")
if(!check[4]) install.packages("data.table")
rm(pkgs, check)
library(simmer)
library(simmer.plot)

## Loading required package: ggplot2

## 
## Attaching package: 'simmer.plot'

## The following objects are masked from 'package:simmer':
## 
##     get_mon_arrivals, get_mon_attributes, get_mon_resources

library(simmer.bricks)
library(data.table)

# User Defined Function To Log Both Resource & Activity Name (Using Attributes) 
VisitStep <- function(activity, duration = 0, resource = "0", n = 0) {
  t1 = trajectory() %>% set_attribute(activity, 0)     # enter the queue
  t2 = trajectory() %>% seize(resource, n)             # seize the resource
  t3 = trajectory() %>% set_attribute(activity, 1)     # begin the activity
  t4 = trajectory() %>% timeout(duration)              # activity duration
  t5 = trajectory() %>% release(resource, n)           # release the resource
  t6 = trajectory() %>% set_attribute(activity, 2)     # leave the activity
  if(resource == "0") return(join(t1,t3,t4,t6)) 
  return(join(t1,t2,t3,t4,t5,t6)) 
}

# Number of Arrivals Complete, Throughput, PCE, Histogram of Process Time, Flow Time
SimSummary <- function(x) {
  id = unique(x[Activity == "Exit"]$ID)
  x = x[ID %in% id]
  y = x[Activity == "Exit"]
  TH = length(id)/(max(y$Finish)-min(y$Finish))
  x[, ProcTime := Finish - Begin]
  x[, WaitTime := Begin - Arrival]
  x[, FlowTime := Finish - Arrival]
  x = x[, .(ProcTime = sum(ProcTime), 
            WaitTime = sum(WaitTime),
            FlowTime = sum(FlowTime)),
        by = ID]
  cat("\n",
      paste("Number of Arrivals Processed :", length(id)), "\n",
      paste("Avg. Throughput              :", round(TH,2)), "\n",
      paste("Process Cycle Efficiency     :", round(sum(x$ProcTime)/sum(x$FlowTime),2)), "\n","\n")
  
  p = data.frame(time = x$ProcTime); p$name = "ProcTime"
  f = data.frame(time = x$FlowTime); f$name = "FlowTime"
  ggplot(rbind(p,f), aes(time, fill = name)) + geom_density(alpha = 0.4)
}

# Step Plot for Arrivals & Departures
EventLogPlot <- function(x) {
  xentr = x[Activity == "Enter"]
  xexit = x[Activity == "Exit"]
  n = trunc(max(xexit$Finish))+1
  dt = data.frame(time = 1:n)
  dt$arrivals = 0; dt$departures = 0
  for(i in 0:n) {
    dt$arrivals[i] = nrow(xentr[Arrival <= i])
    dt$departures[i] = nrow(xexit[Finish <= i])
  }
  ggplot(dt, aes(time)) +
    geom_step(aes(y = arrivals), colour="blue") +
    geom_step(aes(y = departures), colour="red")    
}

# Dumbbell Chart for Arrivals & Departures
EventLogPlot2 <- function(x) {
  id = unique(x[Activity == "Exit"]$ID)
  x = x[ID %in% id]
  x = x[, .(Arr = min(Arrival), 
            Dep = max(Finish)),
        by = ID]
  FlowTime = x$Dep - x$Arr
  n = length(unique(x$ID))
  xmax = (trunc(max(x$Dep)/10)+1)*10
  plot(x$Dep, col = "white",
       xlab = "Time",
       ylab = "Arrivals",
       yaxt = "n",
       xlim = c(0,xmax),
       ylim = c(0,n+1))
  for(i in 1:n) {
    j = (n+1-i)
    lines(c(x$Arr[i],x$Dep[i]),c(j,j), col = "grey50")
    points(x$Arr[i], j, pch = 20, cex = 0.8, col = "blue")
    points(x$Dep[i], j, pch = 15, cex = 0.6, col = "maroon")
  }
  return(FlowTime)
}

# Dumbbell Plot Given Environment Via Arrivals Table
ArrivalsPlot <- function(env) {
  x = data.table(get_mon_arrivals(env, per_resource = TRUE, ongoing = FALSE))
  x = x[order(start_time)]
  x = x[, .(Arr = min(start_time), 
            Dep = max(end_time)),
        by = name]
  n = length(unique(x$name))
  xmax = (trunc(max(x$Dep)/10)+1)*10
  plot(x$Dep, col = "white",
       xlab = "Time",
       ylab = "Arrivals",
       yaxt = "n",
       xlim = c(0,xmax),
       ylim = c(0,n+1))
  for(i in 1:n) {
    j = (n+1-i)
    lines(c(x$Arr[i],x$Dep[i]),c(j,j), col = "grey50")
    points(x$Arr[i], j, pch = 20, cex = 0.8, col = "blue")
    points(x$Dep[i], j, pch = 15, cex = 0.6, col = "maroon")
  }
}

Sys.setenv(TZ = 'GMT')

cat("\014")

we will create statistical variables which are required for simulation & we will also create the parameters which necessary for the simulation process.These are the different parameters we define ‘Arrival’,‘RegTime’, ‘SwabTime’, ‘CleanTime’, ‘TestTime’, ‘SimTime’. SimTime referred duration of simulation which is 240 minutes.

Parameters Required for Simulation

set.seed(1234)
Arrival = function() {return(rexp(1, 1/2.5))}          # Inter Arrival Time - Exponential Distribution
RegTime = function() {return(runif(1, 0.5, 1.5))}      # Registration Time - Uniform Distribution
SwabTime = function() {return(rnorm(1, 2, 1/3))}       # Swab Taking Time - Normal Distribution
CleanTime = function() {return(rnorm(1, 2, 0.33))}     # Bay Cleaning Time - Normal Distribution
TestTime = function() {return(rnorm(1, 1.5, 0.25))}    # Test Process Time - Normal Distribution
SimTime = 240

A. Define Trajectory

The process introduces the trajectory i.e., people entering to register then pick up the test kit & then take the swab. They are entering to the test staff to return the swab. Go back to the bay, clean the bay and leave. Now, after this process ends, there is a parallel path here. There’s a split in the process where the return swab is processed and the performance of the test is done and the result is sent. In this case we’re going through this step register, pick up test, take swab, return swab and clean bay and leave. Our concern doesn’t lie to branch off the process. The discrete event simulation has performed so that there is a single simple trajectory which ends with clean bay. Then people can leave. That’s what I’m going to present now.

We also need to define the resource is the register on the test staff and the base. and we also define the arrivals that arrive inside the process.

 set.seed(1234) 
 # A Define Process Trajectory
  t = trajectory() %>%
  # Registration Step
  seize("Registrar", 1) %>%  
  timeout(RegTime) %>%    
  release("Registrar",1) %>%
  
  # Pickup Test Kit
  seize("Test Staff", 1) %>%   
  timeout(0.5) %>%
  release("Test Staff", 1) %>% 
  
  # Take Swab, Return Swab & Clean Bay Steps
  seize("Bay", 1) %>%  
  timeout(SwabTime) %>%
  timeout(0.5) %>%       
  timeout(CleanTime) %>%
  release("Bay",1)

1) What is Trajectory:

Here we have a trajectory defined with the t. So, t is the variable into which all the trajectories defined every step within the processes defined separately as we written above. And they’re all connected with this pipe operator.so if you want to change the process later on for something else, you could simply use this template. But you need to ensure that this syntax followed properly.

2) Registration Step:

The Registration step is essential. When we go to the registration step, a registrar is seized. That is, they engaged and a particular some amount of time is elapsed & once the registration process is done, the registrar is released. That means they are available for the next person to come and engaged. So, we have registrar seized and then time out registration time. So, this Regtime is actually defined previously. So, it gets a statistical variable and stores in the Regtime time value, and every time you call it, it gives a different value. But statistically, according to the distribution that we have already defined previously. & then we release the register. So, this is the registration step.

3) Pickup Test Kit:

The next is Pickup test kit. If you look at the process, the test staff is engaged for this step, seize is the engaging command. So, you seize the test stuff your timeout a half a minute that we have fixed it. timeout written over here is not statistical variable we fixed it here. Then we release the test staff.

4) Take Swab, Return Swab & Clean Bay Steps:

Next, we take the swab, returned the swab and clean the bay. So, when we take the swab, obviously, before taking the swab, you have to engage the bay. Because the swabs taken in the bay when you go there, you engage the way you seize the bay out of the four. You seize one of the bay, perform this task & then don’t release it without releasing it go back to the test staff. Ideally, we engage the test staff, finished the written swamp step and release the test staff. Go back to the bay, clean the bay and release the bay and leave.

##Print & Plot Trajectory
  print(t)

## trajectory: anonymous, 11 activities
## { Activity: Seize        | resource: Registrar, amount: 1 }
## { Activity: Timeout      | delay: function() }
## { Activity: Release      | resource: Registrar, amount: 1 }
## { Activity: Seize        | resource: Test Staff, amount: 1 }
## { Activity: Timeout      | delay: 0.5 }
## { Activity: Release      | resource: Test Staff, amount: 1 }
## { Activity: Seize        | resource: Bay, amount: 1 }
## { Activity: Timeout      | delay: function() }
## { Activity: Timeout      | delay: 0.5 }
## { Activity: Timeout      | delay: function() }
## { Activity: Release      | resource: Bay, amount: 1 }

  plot(t, verbose = TRUE)

B. Define Resources, Arrivals:

The next step is to define resources & arrivals. So before defining resources an arrival we need to create simulation environment.

1) Create Simulation Environment:

A simulation environment is defined as a programming environment of a computer, that is dedicated to systems simulation and that takes care for a flexible an intelligent interfacing between a user and the system to be experimentally studied.

# B. Define Resources & Arrivals
set.seed(1234)
# Create Simulation Environment
env <- simmer()

# Define Resources
env %>% add_resource("Registrar", 1)

## simmer environment: anonymous | now: 0 | next: 
## { Monitor: in memory }
## { Resource: Registrar | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }

env %>% add_resource("Bay", 4)

## simmer environment: anonymous | now: 0 | next: 
## { Monitor: in memory }
## { Resource: Registrar | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Resource: Bay | monitored: TRUE | server status: 0(4) | queue status: 0(Inf) }

env %>% add_resource("Test Staff", 1)

## simmer environment: anonymous | now: 0 | next: 
## { Monitor: in memory }
## { Resource: Registrar | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Resource: Bay | monitored: TRUE | server status: 0(4) | queue status: 0(Inf) }
## { Resource: Test Staff | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }

# Generate Arrivals
env %>% add_generator("Student", t, Arrival)

## simmer environment: anonymous | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: Registrar | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Resource: Bay | monitored: TRUE | server status: 0(4) | queue status: 0(Inf) }
## { Resource: Test Staff | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Source: Student | monitored: 1 | n_generated: 0 }

Here we’ve defined the arrivals. Here We’ve used the pipe operator to join the add_generator() command into environment. Then students or customers are coming into the environment, that’s why it’s written in text “Student”. We need to define their arriving in to the trajectory that we’ve defined previously i.e., t. Here we need to provide inter arrival time of that and therefore this Arrival is the variable that we have defined previously.

C. Run Simulation:

# C. Run Simulation
set.seed(1234)
env %>% run(until = SimTime)

## simmer environment: anonymous | now: 240 | next: 240.035296530097
## { Monitor: in memory }
## { Resource: Registrar | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Resource: Bay | monitored: TRUE | server status: 4(4) | queue status: 0(Inf) }
## { Resource: Test Staff | monitored: TRUE | server status: 1(1) | queue status: 0(Inf) }
## { Source: Student | monitored: 1 | n_generated: 109 }

D. Store Results:

Now our job is to store the results. So, we have used the functions get_mon_resources() for resources and get_mon_arrivals() for arrivals. These two functions are available within the R-simmer. You just need to provide the environment (env) into get_mon_resources().

# D. Store Results
set.seed(1234)
resources = get_mon_resources(env)
arrivals = get_mon_arrivals(env, per_resource = TRUE)

# Display Arrivals Table
head(arrivals)

##       name start_time  end_time activity_time   resource replication
## 1 Student0   6.254397  7.615312     1.3609154  Registrar           1
## 2 Student0   7.615312  8.115312     0.5000000 Test Staff           1
## 3 Student1   6.871294  8.347862     0.7325505  Registrar           1
## 4 Student1   8.347862  8.847862     0.5000000 Test Staff           1
## 5 Student2   7.572847  9.541454     1.1935913  Registrar           1
## 6 Student2   9.541454 10.041454     0.5000000 Test Staff           1

# Display path of single individual
subset(arrivals, name == "Student8")

##        name start_time end_time activity_time   resource replication
## 22 Student8   26.23661 27.30139     0.9849912  Registrar           1
## 23 Student8   27.30139 27.80139     0.5000000 Test Staff           1
## 31 Student8   27.80139 31.76554     3.9641504        Bay           1

arrivals = get_mon_arrivals(env, per_resource = FALSE)
hist(arrivals$end_time - arrivals$start_time)

Note:

if, the difference between end_time & start_time is greater than activity time (end_time - start_time ≥ activity_time) it refers; that particular student is waiting for that resource before the activity gets done.

E. Plot Charts & Interpret Results:

Now we are going to interpret our results of simulation by plotting charts. So, first thing we’re going to do is to plot dumbbell chart of arrivals. This gives you the intuitive view of results of the simulation. ### I. Plot Dumbbell Chart of Arrivals (Custom Function)

####N.B: This arrivals plot is a custom function. It means that this function isn’t provided R-simmer library. The function is written in R Simmer Functions.When we use the source fucntion then all the functions written in that R script will be loaded.

# Plot Dumbbell Chart of Arrivals (Custom Function)
set.seed(1234)
ArrivalsPlot(env)

Here the plot is all about arrival VS time where 0 to 240 is the simulation time. The blue is the starting point. More precisely, the blue denotes the arrival time of different students that arrive in various time. And red dots are the exits. The line between blue dots & red dots is path. This dumbbell chart gives the idea how simulation is going on. There is a random statistical nature of the arrivals because we use statistical variables, not uniformed variables. they come like this randomly according to that average value according to the exponential distribution and they get processed and they go off here. So this is a statistical way in which we can show the results.

II. Plot Utilization Chart

# Plot Utilization Chart
set.seed(1234)
plot(resources, metric = "utilization", c("Registrar", "Test Staff", "Bay") )

So, these are the resources that are being employed in this simulation. It means that the total time available during the 240 minutes like the registrar is busy 40% of the time test stuff 20% and Bay are busy 40% of the time. In case of Bay, it’s not 40% of time here. We need to multiplied four times. So, we have multiplied 4 times with 240. All the bays are busy. And this is the level of business of the resources. If the resources of the fully busy like they’re literally Working end to end without any gap off breathing time in between. Then the resource utilization will go 100%. When resources will go up all the utilization then queues will form. It means that incoming students will wait for the resource to become available and therefore queues will form as the queues form waiting time will increase because of which the total flow time will also be more. so that’s what this resource utilization is.

III. Plot Usage Chart (Average)

# Plot Usage Chart (Average)
set.seed(1234)
plot(resources, metric = "usage", "Registrar")

This is the resource usage. Here we have three different colored line. It says in use how much percentage is in use in a system, server and queue. So, server meaning the resource is serving. And while the resource is serving, that much (shown in figure) is queuing in front of the resource. So, in total, what is queuing & what is being served is the total on the blue line. The blue line is the sum of these red line and the green line. And the red line is the queue.

IV. Plot Usage Chart (Instantaneous)

# Plot Usage Chart (Instantaneous) 
plot(resources, metric = "usage", "Registrar", items = "server", steps = TRUE)

This is the actual usage. So, the intuitive meaning of usage is either zero or one. This is the step wise depiction of that usage. This is instantaneous.

Advanced Trajectory - Parallel Paths

F. Model Parallel Paths Using Clone Synchronize Function:

Here we will look at advanced trajectory for parallel paths. Before dive into execute advanced trajectory, first we should understand why do we need advanced trajectory & parallel paths. Because you see, the student comes from register. Pick up test, take swab returns, swab. At this point, the student returns to the bay and leaves while the students that sample that he has returned will go through a parallel path at off perform test and then to the result and therefore the process needs to split there. So, we want to model now the trajectory properly. So, we will model the trajectory into three segments. The first segment from here register pickup test takes for and returns for these first four steps with the model as t1. t2 will be Clean Bay and then goes to leave and t3 will we perform tests. And we will join up all the three trajectories t1,t2,t3 into a single trajectory using couple of commands which actually helps us to split the incoming student into two and go into two different parts. That’s the intent off this particular task.

First, we reset the environment otherwise the command will not work.

# Reset Simulation 
set.seed(1234)
reset(env)

## simmer environment: anonymous | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: Registrar | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Resource: Bay | monitored: TRUE | server status: 0(4) | queue status: 0(Inf) }
## { Resource: Test Staff | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Source: Student | monitored: 1 | n_generated: 0 }

Now our task is to build up the trajectories. Here we’ve used visit command instead of using Seize, Timeout, Release.The approach is like seize the resource, wait for some time (which is timeout) , then release the resource. All these three steps are encapsulated in a single command visit.

Advanced Trajectory - Parallel Paths

set.seed(1234)
t1 = trajectory() %>%
  
  # Registration Step
  visit("Registrar", RegTime, 1) %>%
  
  # Pickup Test Kit
  visit("Test Staff", 0.5, 1) %>%
  
  # Take Swab (Note - Bay Not Released)
  seize("Bay", 1) %>%  
  timeout(SwabTime) %>%
  
  # Return Swab
  visit("Test Staff", 0.5, 1) 

# Here the path splits into parallel routes

t2 = trajectory() %>%
  
  # Clean Bay
  timeout(CleanTime) %>%
  
  # Release Bay (seized in t1)
  release("Bay",1) 

t3 = trajectory() %>%
  # Perform Test
  visit("Test Staff", TestTime, 1) 


# Join Trajectories
t = trajectory() %>%
  
  # path before the split
  join(t1) %>%
  
  # split into parallel paths t2, t3
  clone(n = 2, t2, t3) %>%
  
  # join parallel paths 
  synchronize(wait = TRUE) 

print(t)

## trajectory: anonymous, 18 activities
## { Activity: Seize        | resource: Registrar, amount: 1 }
## { Activity: Timeout      | delay: function() }
## { Activity: Release      | resource: Registrar, amount: 1 }
## { Activity: Seize        | resource: Test Staff, amount: 1 }
## { Activity: Timeout      | delay: 0.5 }
## { Activity: Release      | resource: Test Staff, amount: 1 }
## { Activity: Seize        | resource: Bay, amount: 1 }
## { Activity: Timeout      | delay: function() }
## { Activity: Seize        | resource: Test Staff, amount: 1 }
## { Activity: Timeout      | delay: 0.5 }
## { Activity: Release      | resource: Test Staff, amount: 1 }
## { Activity: Clone        | n: 2 }
##   Fork 1, continue,  trajectory: anonymous, 2 activities
##   { Activity: Timeout      | delay: function() }
##   { Activity: Release      | resource: Bay, amount: 1 }
##   Fork 2, continue,  trajectory: anonymous, 3 activities
##   { Activity: Seize        | resource: Test Staff, amount: 1 }
##   { Activity: Timeout      | delay: function() }
##   { Activity: Release      | resource: Test Staff, amount: 1 }
## { Activity: Synchronize  | wait: 1 }

plot(t, verbose = TRUE)

Here t is joined up trajectories by defining empty trajectory using pipe operator & to which we add this t1 using join command. Once t1 is done we can’t join t2 & then t3. Because t2 & t3 is a parallel path. So, using clone command we split t2 & t3 into two parallel paths. Once they have split, they have gone through the respective paths then they have to join back. synchronize(wait = TRUE) is the command to join back the parallel paths.

set.seed(1234)
# Create Environment, Add Resources, Generator (Arrivals)
env <- simmer() %>%
  add_resource("Registrar", 1)  %>%
  add_resource("Bay", 4)   %>%
  add_resource("Test Staff", 1) %>%
  add_generator("Student", t, Arrival)

# Run Simulation
env %>% run(until = SimTime)

## simmer environment: anonymous | now: 240 | next: 240.497373316456
## { Monitor: in memory }
## { Resource: Registrar | monitored: TRUE | server status: 1(1) | queue status: 0(Inf) }
## { Resource: Bay | monitored: TRUE | server status: 4(4) | queue status: 24(Inf) }
## { Resource: Test Staff | monitored: TRUE | server status: 1(1) | queue status: 6(Inf) }
## { Source: Student | monitored: 1 | n_generated: 118 }

Plot Dumbbell Chart of Arrivals:

set.seed(1234)
#Plot Dumbbell Chart of Arrivals (Custom Function)
ArrivalsPlot(env)

This the result of this simulation. Here you can see the blue lines are the arrivals & as you progress in time the length of stay within the system is actually increasing because the rate at which the arrivals are coming are faster than the rate at which departures are leaving.

Store Results:

# Store Results
resources = get_mon_resources(env)
arrivals = get_mon_arrivals(env, per_resource = TRUE)
arrivals2 = get_mon_arrivals(env)
subset(arrivals, name == "Student0")

##        name start_time  end_time activity_time   resource replication
## 1  Student0   6.254397  7.615312      1.360915  Registrar           1
## 2  Student0   7.615312  8.115312      0.500000 Test Staff           1
## 7  Student0  10.258353 10.758353      0.500000 Test Staff           1
## 8  Student0  10.758353 12.102739      1.344386 Test Staff           1
## 10 Student0   8.115312 13.082877      4.967565        Bay           1

Plot Flow Time:

# Plot Flow Time
plot(arrivals2, metric="flow_time")

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Here in this dataset start time and end time for every student. There is no resource, it is not by resource and therefore we can just subtract start time from the end time and plot the Flow time.

G. Capturing Process Step Data & Advanced Visualization:

head(arrivals)

##       name start_time  end_time activity_time   resource replication
## 1 Student0   6.254397  7.615312     1.3609154  Registrar           1
## 2 Student0   7.615312  8.115312     0.5000000 Test Staff           1
## 3 Student1   6.871294  8.347862     0.7325505  Registrar           1
## 4 Student1   8.347862  8.847862     0.5000000 Test Staff           1
## 5 Student2   7.572847  9.541454     1.1935913  Registrar           1
## 6 Student2   9.541454 10.041454     0.5000000 Test Staff           1

Here name of the resources, start time, end time is given. But if you go back to the process map, you’ll see that test staff has three different activities. So, there is no way for us to figure out what activity is going on at any point in time with this sort of table. And therefore, we want to capture the activity names. So that we will go for different approach.

As we have performed simulation work previously, that’s why we need to reset the simulation again.

Reset Simulation: reset(env)

Here we have used visitstep() function to execute this task. This visitStep function is a user defined function. When we load the source file it will be loaded itself.

Syntax:

VisitStep(Activity_Name, Duration, Resource, Res_Qty)

visit step function can be used like this, you provide the activity name first, the duration of the activity, the resource required to engage in that activity and the quantity of the resource (like whether do you need one resource two resource etc.).

VisitStep(Activity_Name, Duration)

You can also use this partially activity name and duration without mentioning the resource. When you don’t mention the resource, no resource will be used but the activity name will be logged and the duration which will be mentioned within the parenthesis which is VisitStep(Activity_Name, Duration) will be spent. That’s how it works.

Log Activity Names - “VisitStep” Function:

# Log Activity Names - "VisitStep" Function
# *********************************************************
set.seed(1234)
# Reset Simulation 
reset(env)

## simmer environment: anonymous | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: Registrar | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Resource: Bay | monitored: TRUE | server status: 0(4) | queue status: 0(Inf) }
## { Resource: Test Staff | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Source: Student | monitored: 1 | n_generated: 0 }

# Student Trajectory - Alternative Approach (Easier)
# "VisitStep" is a user defined function 
# VisitStep(Activity_Name, Duration, Resource, Res_Qty)
# VisitStep(Activity_Name, Duration)

t1 = join(
  VisitStep("Enter", 0),           # Always use "Enter" step
  VisitStep("Register", RegTime, "Registrar", 1),
  VisitStep("Pick TestKit", 0.5, "Test Staff", 1),
  trajectory () %>% seize("Bay", 1),
  VisitStep("Take Swab", SwabTime),
  VisitStep("Return Swab", 0.5, "Test Staff", 1)
)

For VisitStep(“Take Swab”, SwabTime) we have not mentioned any resource because I mentioned the resource in the previous step using the seize function. And we’ve seized one Bay which is a resource. Because, without using this( trajectory () %>% seize(“Bay”, 1),) approach If I put the bay here( VisitStep(“Take Swab”, SwabTime),)within this command then at the end of this command it will release the bay and we don’t want it to release the bay. as in the previous situation we occupy the bay at the time of the take swabs. But we don’t release the Bay. Come back to the test stuff. We return the swab and go back to the bay where we were and clean the bay and then leave. Therefore, the bay is not released. And for that reason, we basically used this partial command as in the second approach, and put in additional seize bay here.

t2 = join(
  VisitStep("Clean Bay", CleanTime),
  trajectory () %>% release("Bay", 1)
)

In t2 is basically clean the bay. Here also we have not provided Bay as the resource. We only spend the time to clean the bay and release the bay

t3 = VisitStep("Perform Test", TestTime, "Test Staff", 1)


# Join Trajectories
t = trajectory() %>%
  
  # path before the split
  join(t1) %>%
  
  # split into parallel paths t2, t3
  clone(n = 2, t2, t3) %>%
  
  # join parallel paths 
  synchronize(wait = TRUE) %>%
 
  join(VisitStep("Exit", 0))         # Always use "Exit" step

Whatever task performed here is already shown previously. The additional thing which is performed join(VisitStep(“Exit”, 0))

# ** Use "mon = 2" in add_generator function to log Activity Names **
set.seed(1234)
# Create Environment, Add Resources, Generator (Arrivals)
env <- simmer() %>%
  add_resource("Registrar", 1)  %>%
  add_resource("Bay", 4)   %>%
  add_resource("Test Staff", 1) %>%
  add_generator("Student", t, Arrival, mon = 2) # mon = 2 to log attributes

All of this is exactly what we have done previously. Except that we have created the trajectory in a different way and the reason is we want to create what is known as an event log. We want to log the activity names. So, we have created event log with the activity names.

set.seed(1234)
# Run Simulation
env %>% run(until = SimTime)

## simmer environment: anonymous | now: 240 | next: 240.3083143571
## { Monitor: in memory }
## { Resource: Registrar | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
## { Resource: Bay | monitored: TRUE | server status: 4(4) | queue status: 10(Inf) }
## { Resource: Test Staff | monitored: TRUE | server status: 1(1) | queue status: 4(Inf) }
## { Source: Student | monitored: 2 | n_generated: 100 }

# Plot Dumbbell Chart of Arrivals (Custom Function)
ArrivalsPlot(env)

# Store Results
resources = get_mon_resources(env)
arrivals = get_mon_arrivals(env, per_resource = TRUE)
arrivals2 = get_mon_arrivals(env, per_resource = FALSE)

Create Event Log With Activity Names:

eventlog = dcast(data.table(get_mon_attributes(env)), 
                 name + key ~ value, value.var = "time",
                 fun.aggregate = sum)
names(eventlog) = c("ID","Activity", "Arrival","Begin","Finish")
eventlog = eventlog[order(Arrival)]
print(eventlog)

##             ID     Activity    Arrival      Begin   Finish
##   1:  Student0        Enter   6.254397   6.254397 6.254397
##   2:  Student0     Register   6.254397   6.254397 7.615312
##   3:  Student1        Enter   6.871294   6.871294 6.871294
##   4:  Student1     Register   6.871294   7.615312 8.347862
##   5:  Student2        Enter   7.572847   7.572847 7.572847
##  ---                                                      
## 727: Student84 Perform Test 238.308314   0.000000 0.000000
## 728: Student85    Clean Bay 238.808314 238.808314 0.000000
## 729: Student85 Perform Test 238.808314   0.000000 0.000000
## 730: Student86    Clean Bay 239.308314 239.308314 0.000000
## 731: Student86 Perform Test 239.308314   0.000000 0.000000

So, we have the student ID. Students are coming in and the activity (entire register, pick up test) arrival, (Arrival time), begin (begin time) and Finish (finish time).The begin time i.e., you arrive, but you don’t begin the step, meaning you’re waiting in front of that process.

SIMULATION SUMMARY:

SimSummary(eventlog)

## 
##  Number of Arrivals Processed : 84 
##  Avg. Throughput              : 0.37 
##  Process Cycle Efficiency     : 0.38 
##

Here 84 people are processed. This is Avg. Throughput: 0.37 it means how many people can go through the process and complete and come out of the end in a minute.

Process Cycle Efficiency: 0.38 it means that 38% all the time spent inside by each student spend in the actual activity and the remaining time is waiting time.

So, for example, if Process Cycle Efficiency: 30%, it means that 0.3% is the real activity time and the remaining 70% is of waiting time.

It also produces a chart our Flow Time and Process Time. You can see the total processing time is a lot less compared to the total flow time, which is due to the processing time and waiting.

Step Plot of Arrivals & Departures:

EventLogPlot(eventlog)

We can also visualize it using this event log plot. and there is this blue line is the rate at which the arrivals coming in and the red line is the rate at which the arrival saw going out. If they stay parallel more or less closely, then the process is alright. But if they diverge, then we have a concern, especially the red line is below the blue line. Then it means that we are diverging and so that’s a concern.

Dumbbell Plot of Arrivals & Departures:

We’re going to plot this eventlogplot2() which is nothing just dumbbell plot. It outputs the flow time.

EventLogPlot2(eventlog)

##  [1]  6.828480  7.490039  8.312235  6.064584  8.706775  8.414767 10.978788
##  [8] 11.479597 10.905037 12.499068 10.520380  6.079888  5.936559  6.942756
## [15] 10.550767 12.294309 11.965431 14.319155 15.761152 13.514162 15.065991
## [22] 14.208768 17.928792 16.780212 18.083878 17.591043 21.472383 19.212524
## [29] 18.698697 20.117190 21.816977 21.955645 20.048484 20.256561 20.868571
## [36] 19.095199 20.077576 24.276871 27.459900 28.441094 21.613291 29.141425
## [43] 29.557627 25.244982 27.051990 29.149391 25.707625 27.026022 26.518318
## [50] 28.976521 25.984659 26.656100 26.843217 26.320816 25.979610 26.304081
## [57] 27.393285 32.840580 33.615532 34.725363 33.510939 37.417060 21.200504
## [64] 22.900074 21.983748 19.393120 18.923129 19.580881 22.031330 23.578413
## [71] 22.040893 25.082350 22.697200 25.208101 25.476748 27.271438 30.507511
## [78] 29.398591 30.168234 29.893776 32.394951 33.511583 34.758529 35.968368

Plot Flow Time:

FlowTime = EventLogPlot2(eventlog) #Here we have stored the above output

plot (FlowTime, type = "l", col = "blue")

Here index is the number of students that are coming in i.e., 84 people. The flow time of those students who are spending this much time which is shown in the plot.

Reference:

Abu-Taieh, E., & El Sheikh, A. R. (2007). Commercial Simulation Packages: A Comparative Study. International Journal of Simulation, 8(2), 66–76.
Dearie, J., & Warfield, T. (1976, July 12-14). The development and use of a simulation model of an outpatient clinic. Proceedings of the 1976 Summer computer Simulation Conference, Simulation Council, Washington, DC, (pp. 554-558).
Dexter, F., Macario, A., Traub, R., Hopwood, M., & Lubarsky, D. (1999). An Operating Room Scheduling Strategy to Maximize the Use of Operating Room Block Time: Computer Simulation of Patient Scheduling and Survey of Patients’ Preferences for Surgical Waiting Time. Anesthesia and Analgesia, 89, 7–20. doi:10.1097/00000539-199907000-00003
https://www.nhs.uk/conditions/coronavirus-covid-19/testing/how-to-do-a-test-at-home-or-at-a-test-site/how-to-do-a-rapid-lateral-flow-test/
https://jamanetwork.com/journals/jamanetworkopen/fullarticle/2767513
https://youtu.be/21WQB0E-6-M
Fone, D., Hollinghurst, S., & Temple, M. (2003). Systematic review of the use and value of computer simulation modeling in population health and health care delivery. Journal of Public Health Medicine, 25(4), 325–335. doi:10.1093/pubmed/fdg075
Discrete Event Simulation using R: Oncology Hospital Example

https://www.linkedin.com/pulse/discrete-event-simulation-using-r-oncology-hospital-example-kamal

Simulation of Covid-19 Lateral Test Process Using R simmer

Abhirup Moitra

2022-07-29

Introduction:

Abstract

Objective:

Methodology:

R Simmer Packages & Functions

Parameters Required for Simulation

A. Define Trajectory

1) What is Trajectory:

2) Registration Step:

3) Pickup Test Kit:

4) Take Swab, Return Swab & Clean Bay Steps:

B. Define Resources, Arrivals:

1) Create Simulation Environment:

C. Run Simulation:

D. Store Results:

Note:

E. Plot Charts & Interpret Results:

II. Plot Utilization Chart

III. Plot Usage Chart (Average)

IV. Plot Usage Chart (Instantaneous)

Advanced Trajectory - Parallel Paths

F. Model Parallel Paths Using Clone Synchronize Function:

Advanced Trajectory - Parallel Paths

Plot Dumbbell Chart of Arrivals:

Store Results:

Plot Flow Time:

G. Capturing Process Step Data & Advanced Visualization:

Syntax:

VisitStep(Activity_Name, Duration, Resource, Res_Qty)

VisitStep(Activity_Name, Duration)

Log Activity Names - “VisitStep” Function:

Create Event Log With Activity Names:

SIMULATION SUMMARY:

Step Plot of Arrivals & Departures:

Dumbbell Plot of Arrivals & Departures:

Plot Flow Time:

Reference: