This tutorial is adapted from a tutorial for the Python 2 package
‘SimPy’, here.
Users familiar with SimPy may find this tutorial helpful for
transitioning to simmer.
In this tutorial we model a simple bank with customers arriving at random. We develop the model step-by-step, starting out simply, and producing a running program at each stage.
A simulation should always be developed to answer a specific question; in these models we investigate how changing the number of bank servers or tellers might affect the waiting time for customers.
We first model a single customer who arrives at the bank for a visit, looks around at the decor for a time and then leaves. There is no queueing. First we will assume his arrival time and the time he spends in the bank are fixed.
The arrival time is fixed at 5, and the time spent in the bank is fixed at 10. We interpret ‘5’ and ‘10’ as ‘5 minutes’ and ‘10 minutes’. The simulation runs for a maximum of 100 minutes, or until all the customers that are generated complete their trajectories.
Note where these constants appear in the code below.
library(simmer)
customer <-
trajectory("Customer's path") %>%
log_("Here I am") %>%
timeout(10) %>%
log_("I must leave")
bank <-
simmer("bank") %>%
add_generator("Customer", customer, at(5))
bank %>% run(until = 100)
#> 5: Customer0: Here I am
#> 15: Customer0: I must leave
#> simmer environment: bank | now: 15 | next:
#> { Monitor: in memory }
#> { Source: Customer | monitored: 1 | n_generated: 1 }
bank %>% get_mon_arrivals()
#> name start_time end_time activity_time finished replication
#> 1 Customer0 5 15 10 TRUE 1The short trace printed out by the get_mon_arrivals
function shows the result. The program finishes at simulation time 15
because there are no further events to be executed. At the end of the
visit, the customer has no more actions and no other objects or
customers are active.
Now we extend the model to allow our customer to arrive at a random simulated time though we will keep the time in the bank at 10, as before.
The change occurs in the arguments to the add_generator
function. The function rexp draws from an exponential
distribution with the given parameter, which in this case is
1/5. See ?rexp for more details. We also seed
the random number generator with 10211 so that the same sequence of
random numbers will be drawn every time the script is run.
library(simmer)
set.seed(10212)
customer <-
trajectory("Customer's path") %>%
log_("Here I am") %>%
timeout(10) %>%
log_("I must leave")
bank <-
simmer("bank") %>%
add_generator("Customer", customer, at(rexp(1, 1/5)))
bank %>% run(until = 100)
#> 7.8393: Customer0: Here I am
#> 17.8393: Customer0: I must leave
#> simmer environment: bank | now: 17.8393046225526 | next:
#> { Monitor: in memory }
#> { Source: Customer | monitored: 1 | n_generated: 1 }
bank %>% get_mon_arrivals()
#> name start_time end_time activity_time finished replication
#> 1 Customer0 7.839305 17.8393 10 TRUE 1The trace shows that the customer now arrives at time 7.839305. Changing the seed value would change that time.
Our simulation does little so far. To consider a simulation with several customers we return to the simple deterministic model and add more customers.
The program is almost as easy as the first example (A customer
arriving at a fixed time). The main change is in the
add_generator function, where we generate three customers
by defining three start times. We also increase the maximum simulation
time to 400 in the run function. Observe that we need only
one definition of the customer trajectory, and generate several
customers who follow that trajectory. These customers act quite
independently in this model.
Each customer also stays in the bank for a different, but non-random,
amount of time. This is an unusual case, so it requires an unusual bit
of R code. The timeout function accepts either a constant
waiting time, or a function that is called once per customer and returns
a single value (e.g. rexp(1, 1/10)). So we need to create a
function that returns a different, but specific, value each time it is
called. We define a function called loop to do this. The
loop function returns another function. In the example, we
store that function under the name x. When called, the
function x returns one of (in the example) 7, 10, and 20,
in sequence, wrapping around when it is called a fourth time.
# Function to specify a series of waiting times, that loop around
loop <- function(...) {
time_diffs <- c(...)
i <- 0
function() {
if (i < length(time_diffs)) {
i <<- i+1
} else {
i <<- 1
}
return(time_diffs[i])
}
}
x <- loop(10, 7, 20)
x(); x(); x(); x(); x()
#> [1] 10
#> [1] 7
#> [1] 20
#> [1] 10
#> [1] 7The technical term for the loop function is a ‘closure’.
How closures work is beyond the scope of this vignette; if you wish to
learn more, then Advanced
R by Hadley Wickam has a good explanation.
When we use loop in the timeout function,
we don’t need to assign its output to a name; it can be an ‘anonymous’
function. It will be called whenever a customer is about to wait and
needs to know how long to wait for. Because only three customers are
generated, and their first step is the timeout step, they are assigned
7, 10, and 20 in order.
Note that this code differs from the SimPy in the order that
customers are defined. Here, the first customer to be defined is the one
who arrives at 2 and waits for 7. That is because the arguments to
at() must be in ascending order.
library(simmer)
# Function to specify a series of waiting times in a loop
loop <- function(...) {
time_diffs <- c(...)
i <- 0
function() {
if (i < length(time_diffs)) {
i <<- i+1
} else {
i <<- 1
}
return(time_diffs[i])
}
}
customer <-
trajectory("Customer's path") %>%
log_("Here I am") %>%
timeout(loop(7, 10, 20)) %>%
log_("I must leave")
bank <-
simmer("bank") %>%
add_generator("Customer", customer, at(2, 5, 12))
bank %>% run(until = 400)
#> 2: Customer0: Here I am
#> 5: Customer1: Here I am
#> 9: Customer0: I must leave
#> 12: Customer2: Here I am
#> 15: Customer1: I must leave
#> 32: Customer2: I must leave
#> simmer environment: bank | now: 32 | next:
#> { Monitor: in memory }
#> { Source: Customer | monitored: 1 | n_generated: 3 }
bank %>% get_mon_arrivals()
#> name start_time end_time activity_time finished replication
#> 1 Customer0 2 9 7 TRUE 1
#> 2 Customer1 5 15 10 TRUE 1
#> 3 Customer2 12 32 20 TRUE 1Alternatively, we can create three different customer trajectories for the three waiting times. This is best done by creating an initial template, and then modifying the waiting time for each copy.
Note that the order in which the customers are defined does not matter this time, and we can also name each customer.
library(simmer)
# Create a template trajectory
customer <-
trajectory("Customer's path") %>%
log_("Here I am") %>%
timeout(1) %>% # The timeout of 1 is a placeholder to be overwritten later
log_("I must leave")
# Create three copies of the template
Klaus <- customer
Tony <- customer
Evelyn <- customer
# Modify the timeout of each copy
Klaus[2] <- timeout(trajectory(), 10)
Tony[2] <- timeout(trajectory(), 7)
Evelyn[2] <- timeout(trajectory(), 20)
# Check that the modifications worked
Klaus
#> trajectory: Customer's path, 3 activities
#> { Activity: Log | message: Here I am, level: 0 }
#> { Activity: Timeout | delay: 10 }
#> { Activity: Log | message: I must lea..., level: 0 }
Tony
#> trajectory: Customer's path, 3 activities
#> { Activity: Log | message: Here I am, level: 0 }
#> { Activity: Timeout | delay: 7 }
#> { Activity: Log | message: I must lea..., level: 0 }
Evelyn
#> trajectory: Customer's path, 3 activities
#> { Activity: Log | message: Here I am, level: 0 }
#> { Activity: Timeout | delay: 20 }
#> { Activity: Log | message: I must lea..., level: 0 }
bank <-
simmer("bank") %>%
add_generator("Klaus", Klaus, at(5)) %>%
add_generator("Tony", Tony, at(2)) %>%
add_generator("Evelyn", Evelyn, at(12))
bank %>% run(until = 400)
#> 2: Tony0: Here I am
#> 5: Klaus0: Here I am
#> 9: Tony0: I must leave
#> 12: Evelyn0: Here I am
#> 15: Klaus0: I must leave
#> 32: Evelyn0: I must leave
#> simmer environment: bank | now: 32 | next:
#> { Monitor: in memory }
#> { Source: Klaus | monitored: 1 | n_generated: 1 }
#> { Source: Tony | monitored: 1 | n_generated: 1 }
#> { Source: Evelyn | monitored: 1 | n_generated: 1 }
bank %>% get_mon_arrivals()
#> name start_time end_time activity_time finished replication
#> 1 Tony0 2 9 7 TRUE 1
#> 2 Klaus0 5 15 10 TRUE 1
#> 3 Evelyn0 12 32 20 TRUE 1Again, the simulations finish before the 400 specified in the
run function.
Another change will allow us to have more customers. To make things clearer we do not use random numbers in this model.
The change is in the add_generator function, where we
use a convenience function from_to to create a sequence of
start times for five customers, starting at time 0, with an interarrival
time of 10 between each customer. One idiosyncracy of the syntax is that
no arrival is created on the to time, so we give it as 41,
one unit after the last arrival to be generated. Another is that the
interarrival time must be specified as a function, hence we define a
constant function function() {10}
library(simmer)
customer <-
trajectory("Customer's path") %>%
log_("Here I am") %>%
timeout(12) %>%
log_("I must leave")
bank <-
simmer("bank") %>%
add_generator("Customer", customer, from_to(0, 41, function() {10}))
bank %>% run(until = 400)
#> 0: Customer0: Here I am
#> 10: Customer1: Here I am
#> 12: Customer0: I must leave
#> 20: Customer2: Here I am
#> 22: Customer1: I must leave
#> 30: Customer3: Here I am
#> 32: Customer2: I must leave
#> 40: Customer4: Here I am
#> 42: Customer3: I must leave
#> 52: Customer4: I must leave
#> simmer environment: bank | now: 52 | next:
#> { Monitor: in memory }
#> { Source: Customer | monitored: 1 | n_generated: 5 }
bank %>% get_mon_arrivals()
#> name start_time end_time activity_time finished replication
#> 1 Customer0 0 12 12 TRUE 1
#> 2 Customer1 10 22 12 TRUE 1
#> 3 Customer2 20 32 12 TRUE 1
#> 4 Customer3 30 42 12 TRUE 1
#> 5 Customer4 40 52 12 TRUE 1We now extend this model to allow arrivals at random. In simulation this is usually interpreted as meaning that the times between customer arrivals are distributed as exponential random variates. There is little change in our program. The only difference between this and the previous example of a single customer generated at a random time is that this example generates several customers at different random times.
The change occurs in the arguments to the add_generator
function. The function rexp draws from an exponential
distribution with the given parameter, which in this case is
1/10. See ?rexp for more details. We also seed
the random number generator with 1289 so that the same sequence of
random numbers will be drawn every time the script is run. The 0 is the
time of the first customer, then four random interarrival times are
drawn, and a final -1 stops the generator.
The reason why we cannot use the from_to function here
is that we want to control the number of arrivals that are generated,
rather than the end-time of arrival generation.
library(simmer)
set.seed(1289)
customer <-
trajectory("Customer's path") %>%
log_("Here I am") %>%
timeout(12) %>%
log_("I must leave")
bank <-
simmer("bank") %>%
add_generator("Customer", customer, function() {c(0, rexp(4, 1/10), -1)})
bank %>% run(until = 400)
#> 0: Customer0: Here I am
#> 12: Customer0: I must leave
#> 14.3114: Customer1: Here I am
#> 26.3114: Customer1: I must leave
#> 26.5588: Customer2: Here I am
#> 35.8506: Customer3: Here I am
#> 38.2706: Customer4: Here I am
#> 38.5588: Customer2: I must leave
#> 47.8506: Customer3: I must leave
#> 50.2706: Customer4: I must leave
#> simmer environment: bank | now: 50.2705868901582 | next:
#> { Monitor: in memory }
#> { Source: Customer | monitored: 1 | n_generated: 5 }
bank %>% get_mon_arrivals()
#> name start_time end_time activity_time finished replication
#> 1 Customer0 0.00000 12.00000 12 TRUE 1
#> 2 Customer1 14.31136 26.31136 12 TRUE 1
#> 3 Customer2 26.55882 38.55882 12 TRUE 1
#> 4 Customer3 35.85064 47.85064 12 TRUE 1
#> 5 Customer4 38.27059 50.27059 12 TRUE 1So far, the model has been more like an art gallery, the customers entering, looking around, and leaving. Now they are going to require service from the bank clerk. We extend the model to include a service counter that will be modelled as a ‘resource’. The actions of a Resource are simple: a customer requests a unit of the resource (a clerk). If one is free, then the customer gets service (and the unit is no longer available to other customers). If there is no free clerk, then the customer joins the queue (managed by the resource object) until it is the customer’s turn to be served. As each customer completes service and releases the unit, the clerk can start serving the next in line.
The service counter is created with the add_resource
function. Default arguments specify that it can serve one customer at a
time, and has infinite queueing capacity.
The seize function causes the customer to join the queue
at the counter. If the queue is empty and the counter is available (not
serving any customers), then the customer claims the counter for itself
and moves onto the timeout step. Otherwise the customer
must wait until the counter becomes available. Behaviour of the customer
while in the queue is controlled by the arguments of the
seize function, rather than by any other functions. Once
the timeout step is complete, the release
function causes the customer to make the counter available to other
customers in the queue.
Since the activity trace does not produce the waiting time by
default, this is calculated and appended using the
transform function.
library(simmer)
set.seed(1234)
bank <- simmer()
customer <-
trajectory("Customer's path") %>%
log_("Here I am") %>%
seize("counter") %>%
log_(function() {paste("Waited: ", now(bank) - get_start_time(bank))}) %>%
timeout(12) %>%
release("counter") %>%
log_("Finished")
bank <-
simmer("bank") %>%
add_resource("counter") %>%
add_generator("Customer", customer, function() {c(0, rexp(4, 1/10), -1)})
bank %>% run(until = 400)
#> 0: Customer0: Here I am
#> 0: Customer0: Waited: 0
#> 12: Customer0: Finished
#> 25.0176: Customer1: Here I am
#> 25.0176: Customer1: Waited: 0
#> 27.4852: Customer2: Here I am
#> 27.551: Customer3: Here I am
#> 37.0176: Customer1: Finished
#> 37.0176: Customer2: Waited: 9.53241116646677
#> 44.9785: Customer4: Here I am
#> 49.0176: Customer2: Finished
#> 49.0176: Customer3: Waited: 21.466591599012
#> 61.0176: Customer3: Finished
#> 61.0176: Customer4: Waited: 16.0391307005949
#> 73.0176: Customer4: Finished
#> simmer environment: bank | now: 73.0175860496223 | next:
#> { Monitor: in memory }
#> { Resource: counter | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
#> { Source: Customer | monitored: 1 | n_generated: 5 }
bank %>%
get_mon_arrivals() %>%
transform(waiting_time = end_time - start_time - activity_time)
#> name start_time end_time activity_time finished replication waiting_time
#> 1 Customer0 0.00000 12.00000 12 TRUE 1 0.000000
#> 2 Customer1 25.01759 37.01759 12 TRUE 1 0.000000
#> 3 Customer2 27.48517 49.01759 12 TRUE 1 9.532411
#> 4 Customer3 27.55099 61.01759 12 TRUE 1 21.466592
#> 5 Customer4 44.97846 73.01759 12 TRUE 1 16.039131Examining the trace we see that the first two customers get instant service but the others have to wait. We still only have five customers, so we cannot draw general conclusions.
This is a simple change to the model in that we retain the single service counter but make the customer service time a random variable. As is traditional in the study of simple queues we first assume an exponential service time.
Note that the argument to timeout must be a function,
otherwise it would apply a constant timeout to every customer.
library(simmer)
set.seed(1269)
customer <-
trajectory("Customer's path") %>%
log_("Here I am") %>%
seize("counter") %>%
log_(function() {paste("Waited: ", now(bank) - get_start_time(bank))}) %>%
# timeout(rexp(1, 1/12)) would generate a single random time and use it for
# every arrival, whereas the following line generates a random time for each
# arrival
timeout(function() {rexp(1, 1/12)}) %>%
release("counter") %>%
log_("Finished")
bank <-
simmer("bank") %>%
add_resource("counter") %>%
add_generator("Customer", customer, function() {c(0, rexp(4, 1/10), -1)})
bank %>% run(until = 400)
#> 0: Customer0: Here I am
#> 0: Customer0: Waited: 0
#> 3.68365: Customer1: Here I am
#> 3.89529: Customer2: Here I am
#> 4.91521: Customer0: Finished
#> 4.91521: Customer1: Waited: 1.23156707156024
#> 10.961: Customer3: Here I am
#> 12.4143: Customer4: Here I am
#> 21.0552: Customer1: Finished
#> 21.0552: Customer2: Waited: 17.1598927142276
#> 65.3116: Customer2: Finished
#> 65.3116: Customer3: Waited: 54.3505873236467
#> 80.1134: Customer3: Finished
#> 80.1134: Customer4: Waited: 67.6990166719428
#> 93.7286: Customer4: Finished
#> simmer environment: bank | now: 93.7285504163997 | next:
#> { Monitor: in memory }
#> { Resource: counter | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
#> { Source: Customer | monitored: 1 | n_generated: 5 }
bank %>%
get_mon_arrivals() %>%
transform(waiting_time = end_time - start_time - activity_time)
#> name start_time end_time activity_time finished replication
#> 1 Customer0 0.000000 4.915215 4.915215 TRUE 1
#> 2 Customer1 3.683648 21.055183 16.139968 TRUE 1
#> 3 Customer2 3.895290 65.311551 44.256368 TRUE 1
#> 4 Customer3 10.960963 80.113355 14.801804 TRUE 1
#> 5 Customer4 12.414338 93.728550 13.615196 TRUE 1
#> waiting_time
#> 1 0.000000
#> 2 1.231567
#> 3 17.159893
#> 4 54.350587
#> 5 67.699017This model with random arrivals and exponential service times is an example of an M/M/1 queue and could rather easily be solved analytically to calculate the steady-state mean waiting time and other operating characteristics. (But not so easily solved for its transient behavior.)
When we introduce several counters we must decide on a queue discipline. Are customers going to make one queue or are they going to form separate queues in front of each counter? Then there are complications - will they be allowed to switch lines (jockey)? We first consider a single queue with several counters and later consider separate isolated queues. We will not look at jockeying.
Here we model a bank whose customers arrive randomly and are to be served at a group of counters, taking a random time for service, where we assume that waiting customers form a single first-in first-out queue.
The only difference between this model and the single-server model is
in the add_resource function, where we have increased the
capacity to two so that it can serve two customers at once.
library(simmer)
set.seed(1269)
customer <-
trajectory("Customer's path") %>%
log_("Here I am") %>%
seize("counter") %>%
log_(function() {paste("Waited: ", now(bank) - get_start_time(bank))}) %>%
timeout(function() {rexp(1, 1/12)}) %>%
release("counter") %>%
log_("Finished")
bank <-
simmer("bank") %>%
add_resource("counter", 2) %>% # Here is the change
add_generator("Customer", customer, function() {c(0, rexp(4, 1/10), -1)})
bank %>% run(until = 400)
#> 0: Customer0: Here I am
#> 0: Customer0: Waited: 0
#> 3.68365: Customer1: Here I am
#> 3.68365: Customer1: Waited: 0
#> 3.89529: Customer2: Here I am
#> 4.91521: Customer0: Finished
#> 4.91521: Customer2: Waited: 1.01992474790691
#> 10.961: Customer3: Here I am
#> 12.4143: Customer4: Here I am
#> 19.8236: Customer1: Finished
#> 19.8236: Customer3: Waited: 8.86265226572255
#> 34.6254: Customer3: Finished
#> 34.6254: Customer4: Waited: 22.2110816140186
#> 48.2406: Customer4: Finished
#> 49.1716: Customer2: Finished
#> simmer environment: bank | now: 49.1715828181142 | next:
#> { Monitor: in memory }
#> { Resource: counter | monitored: TRUE | server status: 0(2) | queue status: 0(Inf) }
#> { Source: Customer | monitored: 1 | n_generated: 5 }
bank %>%
get_mon_arrivals() %>%
transform(waiting_time = end_time - start_time - activity_time)
#> name start_time end_time activity_time finished replication
#> 1 Customer0 0.000000 4.915215 4.915215 TRUE 1
#> 2 Customer1 3.683648 19.823616 16.139968 TRUE 1
#> 3 Customer3 10.960963 34.625419 14.801804 TRUE 1
#> 4 Customer4 12.414338 48.240615 13.615196 TRUE 1
#> 5 Customer2 3.895290 49.171583 44.256368 TRUE 1
#> waiting_time
#> 1 0.000000
#> 2 0.000000
#> 3 8.862652
#> 4 22.211082
#> 5 1.019925The waiting times in this model are much shorter than those for the single service counter. For example, the waiting time for Customer3 has been reduced from nearly 17 minutes to less than 9. Again we have too few customers processed to draw general conclusions.
Each counter is now assumed to have its own queue. The programming is more complicated because the customer has to decide which queue to join. The obvious technique is to make each counter a separate resource.
In practice, a customer might join the shortest queue. We implement
this behaviour by first selecting the shortest queue, using the
select function. Then we use seize_selected to
enter the chosen queue, and later release_selected.
The rest of the program is the same as before.
library(simmer)
set.seed(1014)
customer <-
trajectory("Customer's path") %>%
log_("Here I am") %>%
select(c("counter1", "counter2"), policy = "shortest-queue") %>%
seize_selected() %>%
log_(function() {paste("Waited: ", now(bank) - get_start_time(bank))}) %>%
timeout(function() {rexp(1, 1/12)}) %>%
release_selected() %>%
log_("Finished")
bank <-
simmer("bank") %>%
add_resource("counter1", 1) %>%
add_resource("counter2", 1) %>%
add_generator("Customer", customer, function() {c(0, rexp(4, 1/10), -1)})
bank %>% run(until = 400)
#> 0: Customer0: Here I am
#> 0: Customer0: Waited: 0
#> 23.7144: Customer1: Here I am
#> 23.7144: Customer1: Waited: 0
#> 30.4011: Customer2: Here I am
#> 32.4163: Customer3: Here I am
#> 33.941: Customer1: Finished
#> 33.941: Customer3: Waited: 1.52471543798104
#> 39.5302: Customer3: Finished
#> 48.8559: Customer4: Here I am
#> 48.8559: Customer4: Waited: 0
#> 55.3965: Customer4: Finished
#> 62.1037: Customer0: Finished
#> 62.1037: Customer2: Waited: 31.7026170465295
#> 62.3885: Customer2: Finished
#> simmer environment: bank | now: 62.3885266177193 | next:
#> { Monitor: in memory }
#> { Resource: counter1 | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
#> { Resource: counter2 | monitored: TRUE | server status: 0(1) | queue status: 0(Inf) }
#> { Source: Customer | monitored: 1 | n_generated: 5 }
bank %>%
get_mon_arrivals() %>%
transform(service_start_time = end_time - activity_time) %>%
.[order(.$start_time),]
#> name start_time end_time activity_time finished replication
#> 4 Customer0 0.00000 62.10372 62.1037152 TRUE 1
#> 1 Customer1 23.71444 33.94103 10.2265939 TRUE 1
#> 5 Customer2 30.40110 62.38853 0.2848114 TRUE 1
#> 2 Customer3 32.41632 39.53020 5.5891677 TRUE 1
#> 3 Customer4 48.85593 55.39645 6.5405163 TRUE 1
#> service_start_time
#> 4 0.00000
#> 1 23.71444
#> 5 62.10372
#> 2 33.94103
#> 3 48.85593
bank %>%
get_mon_resources() %>%
.[order(.$time),]
#> resource time server queue capacity queue_size system limit replication
#> 1 counter1 0.00000 1 0 1 Inf 1 Inf 1
#> 2 counter2 23.71444 1 0 1 Inf 1 Inf 1
#> 3 counter1 30.40110 1 1 1 Inf 2 Inf 1
#> 4 counter2 32.41632 1 1 1 Inf 2 Inf 1
#> 5 counter2 33.94103 1 0 1 Inf 1 Inf 1
#> 6 counter2 39.53020 0 0 1 Inf 0 Inf 1
#> 7 counter2 48.85593 1 0 1 Inf 1 Inf 1
#> 8 counter2 55.39645 0 0 1 Inf 0 Inf 1
#> 9 counter1 62.10372 1 0 1 Inf 1 Inf 1
#> 10 counter1 62.38853 0 0 1 Inf 0 Inf 1The results show that the customers chose the counter with the smallest number. Unlucky Customer2 who joins the wrong queue has to wait until Customer0 finishes at time 62.10372, and is the last to leave. There are, however, too few arrivals in these runs, limited as they are to five customers, to draw any general conclusions about the relative efficiencies of the two systems.
We now demonstrate how to calculate average waiting times for our
customers. In the original SimPy version of this tutorial, this involved
using ‘Monitors’. In simmer, data is returned by the
get_mon_* family of functions, as has already been
demonstrated. Here, we simply summarise the data frame returned by the
get_mon_arrivals function, using standard R functions.
We also increase the number of customers to 50 (find the number ‘49’ in the code).
library(simmer)
set.seed(100005)
customer <-
trajectory("Customer's path") %>%
seize("counter") %>%
timeout(function() {rexp(1, 1/12)}) %>%
release("counter")
bank <-
simmer("bank") %>%
add_resource("counter", 2) %>%
add_generator("Customer", customer, function() {c(0, rexp(49, 1/10), -1)})
bank %>% run(until = 1000)
#> simmer environment: bank | now: 650.615510598544 | next:
#> { Monitor: in memory }
#> { Resource: counter | monitored: TRUE | server status: 0(2) | queue status: 0(Inf) }
#> { Source: Customer | monitored: 1 | n_generated: 50 }
result <-
bank %>%
get_mon_arrivals() %>%
transform(waiting_time = end_time - start_time - activity_time)The average waiting time for 50 customers in this 2-counter system is more reliable (i.e., less subject to random simulation effects) than the times we measured before but it is still not sufficiently reliable for real-world decisions. We should also replicate the runs using different random number seeds. The result of this run is:
To get a number of independent measurements we must replicate the
runs using different random number seeds. Each replication must be
independent of previous ones, so the environment (bank) must be
redefined for each run, so that the random interarrival times in the
add_generator function are generated from scratch.
We take the chunks of code that build the environment (bank) and run
the simulation, and wrap them in the mclapply function from
the ‘parallel’ package. This function runs each simulation in parallel,
using the available cores in the computer’s processor. Because we use
seeds for reproducability, we pass these to the function that runs the
simulation (function(the_seed)).
library(simmer)
library(parallel)
customer <-
trajectory("Customer's path") %>%
seize("counter") %>%
timeout(function() {rexp(1, 1/12)}) %>%
release("counter")
mclapply(c(393943, 100005, 777999555, 319999772), function(the_seed) {
set.seed(the_seed)
bank <-
simmer("bank") %>%
add_resource("counter", 2) %>%
add_generator("Customer", customer, function() {c(0, rexp(49, 1/10), -1)})
bank %>% run(until = 400)
result <-
bank %>%
get_mon_arrivals() %>%
transform(waiting_time = end_time - start_time - activity_time)
paste("Average wait for ", sum(result$finished), " completions was ",
mean(result$waiting_time), "minutes.")
}) %>% unlist()
#> [1] "Average wait for 49 completions was 3.4958501806687 minutes."
#> [2] "Average wait for 29 completions was 0.584941574937737 minutes."
#> [3] "Average wait for 32 completions was 1.00343842138177 minutes."
#> [4] "Average wait for 42 completions was 7.41231393636088 minutes."The results show some variation. Remember, though, that the system is still only operating for 50 customers, so the system may not be in steady-state.