Terminology

This document uses some DES-specific terminology, e.g., event, state, entity, process or attribute. Such standard terms can be easily found in any textbook about DES (refer to Banks (2005), for instance). There are, however, some -specific terms, and some elements that require further explanation to understand the package architecture.

Resource

A passive entity, as it is commonly understood in standard DES terminology. However, resources are conceived with queuing systems in mind, and therefore they comprise two internal self-managed parts:

Server

which, conceptually, represents the resource itself. It has a specified capacity and can be seized and released.

Queue

A priority queue of a certain size.

Manager

An active entity, i.e., a process, that has the ability to adjust properties of a resource (capacity and queue size) at run-time.

Source

A process responsible for creating new arrivals with a given interarrival time pattern and inserting them into the simulation model.

Arrival

A process capable of interacting with resources or other entities of the simulation model. It may have some attributes and prioritization values associated and, in general, a limited lifetime. Upon creation, every arrival is attached to a given trajectory.

Trajectory

An interlinkage of activities constituting a recipe for arrivals attached to it, i.e., an ordered set of actions that must be executed. The simulation model is ultimately represented by a set of trajectories.

Activity

The individual unit of action that allows arrivals to interact with resources and other entities, perform custom routines while spending time in the system, move back and forth through the trajectory dynamically, and much more.

Architecture

Extending an package (or any other piece of software written in any interpreted language) with compiled code poses an important trade-off between performance and flexibility: placing too much functionality into the compiled part produces gains in performance, but degrades modelling capabilities, and vice versa. The following lines are devoted to discuss how this trade-off is resolved in .

Figure sketches a UML (Unified Modelling Language) description of the architecture, which constitutes a process-based design, as in many modern DES frameworks. We draw the attention now to the classes (depicted in white).

The first main component is the class. It comprises the event loop and the event queue, which will be addressed in the next section. The provides methods for scheduling and unscheduling events. Moreover, it is responsible for managing simulation-wide entities (e.g., resources and arrival sources) and facilities (e.g., signaling between processes and batches) through diverse unordered maps:

• Maps of resources and processes (sources, arrivals and managers) by name.
• A map of pending events, which allows to unschedule a given process.
• Maps of signals subscribed by arrivals and handlers defined for different signals.
• Maps for forming batches of arrivals, named and unnamed.

This class also holds global attributes and monitoring information. Thus, monitoring counters, which are derived from the class, are centralized, and they register every change of state produced during the simulation time. There are five types of built-in changes of state that are recorded by calling ’s methods:

• An arrival is accepted into a resource (served or enqueued). The resource notifies about the new status of its internal counters.
• An arrival leaves a resource. The resource notifies the new status of its internal counters, and the arrival notifies start, end and activity times in that particular resource.
• A resource is modified during runtime (i.e., a change in the capacity or queue size). The resource notifies the new status of its internal counters.
• An arrival modifies an attribute, one of its own or a global one. The arrival notifies the new value.
• An arrival leaves its trajectory by exhausting the activities associated (considered as finished) or because of another reason (non-finished, e.g., it is rejected from a resource). The arrival notifies global start, end and activity times.

As mentioned in the previous section, there are two types of entities: passive ones () and active ones (processes , and ). Sources create new arrivals, and the latter are the main actors of the simulation model. Managers can be used for dynamically changing the properties of a resource (capacity and queue size). All processes share a method that is invoked by the event loop each time a new event is extracted from the event list.

There is a fourth kind of process not shown in Figure , called . It is a generic process that executes a given function once, and it is used by arrivals, resources, activities and the simulator itself to trigger dynamic actions or split up events. A is for instance used under the hood to trigger reneging or to broadcast signals after some delay.

The last main component, completely isolated from the , is the class. This abstract class represents a clonable object, chainable in a double-linked list to form trajectories. Most of the activities provided by derive from it. is another abstract class (not depicted in Figure ) which is derived from . Any activity supporting the definition of sub-trajectories must derive from this one instead, such as , or . All the activities must implement the virtual methods and .

Finally, it is worth mentioning the couple of blue circles depicted in Figure . They represent the points of presence of in the core, i.e., where the core interfaces back with to execute custom user-defined code.

In summary, the core is responsible for all the heavy tasks, i.e., managing the event loop, the event list, generic resources and processes, collecting all the statistics, and so on. And still, it provides enough flexibility to the user for modelling the interarrival times from and execute any custom user-defined code through the activities.

The event queue

The event queue is the most fundamental part of any DES software. It is responsible for maintaining a list of events to be executed by the event loop in an ordered fashion by time of occurrence. This last requirement establishes the need for a data structure with a low access, search, insertion and deletion complexity. A binary tree is a well-known data structure that satisfies these properties, and it is commonly used for this purpose. Unfortunately, binary trees, or equivalent structures, cannot be efficiently implemented without pointers, and this is the main reason why pure is very inefficient for DES.

In , the event queue is defined as a multiset, a kind of associative container implemented as a balanced tree internally. Apart from the efficiency, it was selected to support event unscheduling through iterators. Each event holds a pointer to a process, which will be retrieved and run in the event loop. Events are inserted in the event queue ordered by 1) time of occurrence and 2) priority. This secondary order criterion is devoted to solve a common issue for DES software called the simultaneity problem.

The trajectory object

A trajectory can be defined as a recipe and consists of an ordered set of activities. The idea behind this concept is very similar to the idea behind for data manipulation (Wickham, François, and Müller 2019). To borrow the words of H. Wickham, “by constraining your options, it simplifies how you can think about” discrete-event modelling. Activities are verbs that correspond to common functional DES blocks.

The method instantiates the object, and activities can be appended using the operator:

traj0 <- trajectory() %>%
  log_("Entering the trajectory") %>%
  timeout(10) %>%
  log_("Leaving the trajectory")

The trajectory above illustrates the two most basic activities available: displaying a message () and spending some time in the system (). An arrival attached to this trajectory will execute the activities in the given order, i.e., it will display “Entering the trajectory”, then it will spend 10 units of (simulated) time, and finally it will display “Leaving the trajectory”.

The example uses fixed parameters: a string and a numeric value respectively. However, at least the main parameter for all activities (this is specified in the documentation) can also be what we will call a dynamical parameter, i.e., a function. This thus, although not quite useful yet, is also valid:

traj1 <- trajectory() %>%
  log_(function() "Entering the trajectory") %>%
  timeout(function() 10) %>%
  log_(function() "Leaving the trajectory")

Also, trajectories can be split apart, joined together and modified:

traj2 <- join(traj0[c(1, 3)], traj0[2])
traj2[1] <- traj2[3]
traj2

## trajectory: anonymous, 3 activities
## { Activity: Timeout      | delay: 10 }
## { Activity: Log          | message: Leaving th..., level: 0 }
## { Activity: Timeout      | delay: 10 }

There are many activities available. We will briefly review them by categorizing them into different topics.

env <- simmer()

traj <- trajectory() %>%
  set_attribute("weight", 80) %>%
  set_attribute("weight", 1, mod="+") %>%
  log_(function() paste0("My weight is ", get_attribute(env, "weight")))

patient <- trajectory() %>%
  log_("arriving...") %>%
  seize(
    "doctor", 1, continue = c(TRUE, FALSE),
    post.seize = trajectory("accepted patient") %>%
      log_("doctor seized"),
    reject = trajectory("rejected patient") %>%
      log_("rejected!") %>%
      seize("nurse", 1) %>%
      log_("nurse seized") %>%
      timeout(2) %>%
      release("nurse", 1) %>%
      log_("nurse released")) %>%
  timeout(5) %>%
  release("doctor", 1) %>%
  log_("doctor released")

env <- simmer() %>%
  add_resource("doctor", capacity = 1, queue_size = 0) %>%
  add_resource("nurse", capacity = 10, queue_size = 0) %>%
  add_generator("patient", patient, at(0, 1)) %>%
  run()

## 0: patient0: arriving...
## 0: patient0: doctor seized
## 1: patient1: arriving...
## 1: patient1: rejected!
## 1: patient1: nurse seized
## 3: patient1: nurse released
## 5: patient0: doctor released

traj <- trajectory() %>%
  deactivate("dummy") %>%
  timeout(1) %>%
  activate("dummy")

simmer() %>%
  add_generator("dummy", traj, function() 1) %>%
  run(10) %>%
  get_mon_arrivals()

##     name start_time end_time activity_time finished replication
## 1 dummy0          1        2             1     TRUE           1
## 2 dummy1          3        4             1     TRUE           1
## 3 dummy2          5        6             1     TRUE           1
## 4 dummy3          7        8             1     TRUE           1

env <- simmer()

traj <- trajectory() %>%
  branch(
    option = function() round(now(env)), continue = c(FALSE, TRUE),
    trajectory() %>% log_("branch 1"),
    trajectory() %>% log_("branch 2")) %>%
  clone(
    n = 2,
    trajectory() %>% log_("clone 0"),
    trajectory() %>% log_("clone 1")) %>%
  synchronize(wait = TRUE) %>%
  log_("out")

env %>%
  add_generator("dummy", traj, at(1, 2)) %>%
  run() %>% invisible

## 1: dummy0: branch 1
## 2: dummy1: branch 2
## 2: dummy1: clone 0
## 2: dummy1: clone 1
## 2: dummy1: out

hello <- trajectory() %>%
  log_("Hello!") %>%
  timeout(1) %>%
  rollback(amount = 2, times = 2)

simmer() %>%
  add_generator("hello_sayer", hello, at(0)) %>%
  run() %>% invisible

## 0: hello_sayer0: Hello!
## 1: hello_sayer0: Hello!
## 2: hello_sayer0: Hello!

roller <- trajectory() %>%
  batch(10, timeout = 5, permanent = FALSE) %>%
  seize("rollercoaster", 1) %>%
  timeout(5) %>%
  release("rollercoaster", 1) %>%
  separate()

t_blocked <- trajectory() %>%
  trap(
    "you shall pass",
    handler = trajectory() %>%
      log_("got a signal!")) %>%
  log_("waiting...") %>%
  wait() %>%
  log_("continuing!")

t_signal <- trajectory() %>%
  log_("you shall pass") %>%
  send("you shall pass")

simmer() %>%
  add_generator("blocked", t_blocked, at(0)) %>%
  add_generator("signaler", t_signal, at(5)) %>%
  run() %>% invisible

## 0: blocked0: waiting...
## 5: signaler0: you shall pass
## 5: blocked0: got a signal!
## 5: blocked0: continuing!

bank <- trajectory() %>%
  log_("Here I am") %>%
  renege_in(
    5, out = trajectory() %>%
      log_("Lost my patience. Reneging...")) %>%
  seize("clerk", 1) %>%
  renege_abort() %>%
  log_("I'm being attended") %>%
  timeout(10) %>%
  release("clerk", 1) %>%
  log_("Finished")

simmer() %>%
  add_resource("clerk", 1) %>%
  add_generator("customer", bank, at(0, 1)) %>%
  run() %>% invisible

## 0: customer0: Here I am
## 0: customer0: I'm being attended
## 1: customer1: Here I am
## 6: customer1: Lost my patience. Reneging...
## 10: customer0: Finished

The simulation environment

The simulation environment manages resources and sources, and controls the simulation execution. The method instantiates the object, after which resources and sources can be appended using the operator:

env <- simmer() %>%
  add_resource("res_name", 1) %>%
  add_generator("arrival", traj0, function() 25) %>%
  print()

## simmer environment: anonymous | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: res_name | monitored: TRUE | server status: 0(1) | queue... }
## { Source: arrival | monitored: 1 | n_generated: 0 }

Then, the simulation can be executed, or , a stop time:

env %>% run(until=30)

## 25: arrival0: Entering the trajectory

## simmer environment: anonymous | now: 30 | next: 35
## { Monitor: in memory }
## { Resource: res_name | monitored: TRUE | server status: 0(1) | queue... }
## { Source: arrival | monitored: 1 | n_generated: 2 }

There are a number of methods for extracting information, such as the simulation time (), future scheduled events (), and getters for obtaining resources’ and sources’ parameters (capacity, queue size, server count and queue count; number of arrivals generated so far). There are also several setters available for resources and sources (capacity, queue size; trajectory, distribution).

A object can be and re-run. However, there is a special method, , intended to extract all the information from the object encapsulated into a environment and to deallocate that object. Thus, most of the getters work also when applied to wrapped environments, but such an object cannot be reset or re-run anymore.

Monitoring and data retrieval

There are three methods for obtaining monitored data (if any) about arrivals, resources and attributes. They can be applied to a single simulation environment or to a list of environments, and the returning object is always a data frame, even if no data was found. Each processed simulation environment is treated as a different replication, and a numeric column named is added to every returned data frame with environment indexes as values.

get_mon_arrivals()

Returns timing information per arrival: of the arrival, , , (time not spent in resource queues) and a flag, , that indicates whether the arrival exhausted its activities (or was rejected). By default, this information is referred to the arrivals’ entire lifetime, but it may be obtained on a per-resource basis by specifying .

get_mon_resources()

Returns state changes in resources: name, instant of the event that triggered the state change, count, count, , , count ( + ) and system ( + ).

get_mon_attributes()

Returns state changes in attributes: of the attribute, instant of the event that triggered the state change, name of that identifies the attribute and .

The following sections aim to provide some basic modelling examples. The topics addressed are queuing systems, replication, parallelization and some best practices. We invite the reader to learn about a broader selection of activities and modelling techniques available in the package vignettes, which cover the use of attributes, loops, batching, branching, shared events, reneging and advanced uses of resources among others.

Queuing systems

The concept of trajectory developed in emerges as a natural way to simulate a wide range of problems related to Continuous-Time Markov Chains (CTMC), and more specifically to the so-called birth-death processes and queuing systems. Indeed, not only provides very flexible resources (with or without queue), branches, delays and arrival sources, but they are bundled in a very comprehensive framework of verbs that can be chained with the pipe operator. Let us explore the expressiveness of a trajectory using a traditional queuing example: the M/M/1. The package vignettes include further examples on M/M/c/k systems, queueing networks and CTMC models.

In Kendall’s notation (Kendall 1953), an M/M/1 system has exponential arrivals (M/M/1), a single server (M/M/1) with exponential service time (M/M/1) and an infinite queue (implicit M/M/1/∞). For instance, people arriving at an ATM at rate λ, waiting their turn in the street and withdrawing money at rate μ. These are the basic parameters of the system, whenever ρ < 1:

If ρ ≥ 1, it means that the system is unstable: there are more arrivals than the server is capable of handling and the queue will grow indefinitely. The simulation of an M/M/1 system is quite simple using :

set.seed(1234)
lambda <- 2
mu <- 4
rho <- lambda/mu

mm1.traj <- trajectory() %>%
  seize("mm1.resource", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("mm1.resource", amount=1)

mm1.env <- simmer() %>%
  add_resource("mm1.resource", capacity=1, queue_size=Inf) %>%
  add_generator("arrival", mm1.traj, function() rexp(1, lambda)) %>%
  run(until=2000)

After the parameter setup, the first code block defines the trajectory: each arrival will seize the resource, wait some exponential random time (service time) and release the resource. The second code block instantiates the simulation environment, creates the resource, attaches an exponential generator to the trajectory and runs the simulation for 2000 units of time. Note that trajectories can be defined independently of the simulation environment, but it is recommended to instantiate the latter in the first place, so that trajectories are able to extract information from it (e.g., the simulation time).

As a next step, we could extract the monitoring information and perform some analyses. The extension package (Ucar and Smeets 2019b) provides convenience plotting methods to, for instance, quickly visualize the usage of a resource over time. Figure gives a glimpse of this simulation using this package. In particular, it shows that the average number of customers in the system converges to the theoretical value given by Equation .

Detail of the resource usage.

Replication and parallelization

Typically, running a certain simulation only once is useless. In general, we will be interested in replicating the model execution many times, maybe with different initial conditions, and then perform some statistical analysis over the output. This can be easily achieved using standard tools, e.g., or similar functions.

Additionally, we can leverage the parallelized version of , , provided by the package ( Core Team 2019), to speed up this process. Unfortunately, parallelization has the shortcoming that we lose the underlying objects when each thread finishes. To avoid losing the monitored data, the method can be used to extract and wrap these data into a pure object before the object is garbage-collected.

The following example uses and to perform 100 replicas of the M/M/1 simulation from the previous section (note that the trajectory is not redefined):

library("parallel")
set.seed(1234)

mm1.envs <- mclapply(1:100, function(i) {
  simmer() %>%
    add_resource("mm1.resource", capacity=1, queue_size=Inf) %>%
    add_generator("arrival", mm1.traj, function() rexp(100, lambda)) %>%
    run(until=1000/lambda) %>% wrap()
}, mc.set.seed=FALSE)

With all these replicas, we could, for instance, perform a t-test over N, the average number of customers in the system:

mm1.data <- get_mon_arrivals(mm1.envs) %>%
  aggregate(end_time - start_time ~ replication, data=., mean)
t.test(mm1.data[[2]])

## 
##  One Sample t-test
## 
## data:  mm1.data[[2]]
## t = 94.883, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.4925143 0.5135535
## sample estimates:
## mean of x 
## 0.5030339

Best practices

DES modelling can be done in an event-by-event basis, but this approach is fairly tedious and mostly unpractical. Instead, modern process-oriented approaches commonly relate to the identification of resources and processes in a given problem, and the interactions between them. The package internally follows this paradigm and exposes generic resources and processes (arrivals, in terminology), so that the user can implement all the interactions as action sequences (trajectories).

There are usually multiple valid ways of mapping the identified resources and processes into the elements exposed by the API. For example, let us suppose that we would like to model an alarm clock beeping every second. In this case, the beep may be identified as a process, so that we have different beeps (multiple arrivals) entering a trajectory once per second:

beep <- trajectory() %>%
  log_("beeeep!")

env <- simmer() %>%
  add_generator("beep", beep, function() 1) %>%
  run(2.5)

## 1: beep0: beeeep!
## 2: beep1: beeeep!

But instead, identifying the alarm clock as the process is equally valid, and then we have a single alarm (single arrival) producing all the beeps in a loop:

alarm <- trajectory() %>%
  timeout(1) %>%
  log_("beeeep!") %>%
  rollback(2)

env <- simmer() %>%
  add_generator("alarm", alarm, at(0)) %>%
  run(2.5)

## 1: alarm0: beeeep!
## 2: alarm0: beeeep!

These are two common design patterns in for which the outcome is the same, although there are subtle differences that depend on the problem being considered and the monitoring requirements. Furthermore, as a model becomes more complex and detailed, the resulting mapping and syntax may become more artificious. These issues are shared in different ways by other frameworks as well, such as , and arise due to their generic nature.

Furthermore, the piping mechanism used in the API may invite the user to produce large monolithic trajectories. However, it should be noted that it is usually better to break them down into small manageable pieces. For instance, the following example parameterizes the access to a , where refers to arbitrary service times, and n servers are seized. Then, it is used to instantiate the trajectory shown in the former M/M/1 example:

xgn <- function(resource, G, n)
  trajectory() %>%
    seize(resource, n) %>%
    timeout(G) %>%
    release(resource, n)

(mm1.traj <- xgn("mm1.resource", function() rexp(1, mu), 1))

## trajectory: anonymous, 3 activities
## { Activity: Seize        | resource: mm1.resource, amount: 1 }
## { Activity: Timeout      | delay: function() }
## { Activity: Release      | resource: mm1.resource, amount: 1 }

Standard tools ( and the like) may also be used to generate large lists of trajectories with some variations. These small pieces can be concatenated together into longer trajectories using , but at the same time, they allow for multiple points of attachment of arrivals.

During a simulation, trajectories can interact with the simulation environment in order to extract or modify parameters of interest such as the current simulation time, attributes, status of resources (get the number of arrivals in a resource, get or set resources’ capacity or queue size), or status of sources (get the number of generated arrivals, set sources’ attached trajectory or distribution). The only requirement is that the simulation object must be defined in the same environment (or a parent one) before the simulation is started. Effectively, it is enough to detach the method from the instantiation (), namely, they should not be called in the same pipe. But, for the sake of consistency, it is a good coding practice to instantiate the simulation object always in the first place as follows:

set.seed(1234)
env <- simmer()

traj <- trajectory() %>%
  log_(function() paste0("Current simulation time: ", now(env)))

env <- env %>%
  add_generator("dummy", traj, at(rexp(1, 1))) %>%
  run()

## 2.50176: dummy0: Current simulation time: 2.50175860496223

Comparison with similar frameworks

A significant effort has been put into the design of in order to make it performant enough to run general and relatively large simulation models in a reasonable amount of time. In this regard, a relevant comparison can be made against other general-purpose DES frameworks such as and . To this effect, we retake the M/M/1 example from Section , which can be bundled into the following test:

test_mm1_simmer <- function(n, m, mon=FALSE) {
  mm1 <- trajectory() %>%
    seize("server", 1) %>%
    timeout(function() rexp(1, 1.1)) %>%
    release("server", 1)
  env <- simmer() %>%
    add_resource("server", 1, mon=mon) %>%
    add_generator("customer", mm1, function() rexp(m, 1), mon=mon) %>%
    run(until=n)
}

With the selected arrival rate, λ = 1, this test simulates an average of arrivals entering a nearly saturated system (ρ = 1/1.1). Given that generators are able to create arrivals in batches (i.e., more than one arrival for each function call) for improved performance, the parameter controls the size of the batch. Finally, the flag enables or disables monitoring.

Let us build now the equivalent model using , with base for random number generation. We prepare the benchmark from using the package (Gil Bellosta 2015) as follows:

rPython::python.exec("
import simpy, random, time

def test_mm1(n):
  def exp_source(env, lambd, server, mu):
      while True:
          dt = random.expovariate(lambd)
          yield env.timeout(dt)
          env.process(customer(env, server, mu))
  def customer(env, server, mu):
      with server.request() as req:
          yield req
          dt = random.expovariate(mu)
          yield env.timeout(dt)
  env = simpy.Environment()
  server = simpy.Resource(env, capacity=1)
  env.process(exp_source(env, 1, server, 1.1))
  env.run(until=n)

def benchmark(n, times):
  results = []
  for i in range(0, times):
    start = time.time()
    test_mm1(n)
    results.append(time.time() - start)
  return results
")

Equivalently, this can be done for and using the package (Gong 2017).⁴ Once more, controls the number of arrivals simulated on average:

rjulia::julia_init()
rjulia::julia_void_eval("
using SimJulia, Distributions

function test_mm1(n::Float64)
  function exp_source(env::Environment, lambd::Float64,
                      server::Resource, mu::Float64)
    while true
      dt = rand(Exponential(1/lambd))
      yield(Timeout(env, dt))
      Process(env, customer, server, mu)
    end
  end
  function customer(env::Environment, server::Resource, mu::Float64)
    yield(Request(server))
    dt = rand(Exponential(1/mu))
    yield(Timeout(env, dt))
    yield(Release(server))
  end
  env = Environment()
  server = Resource(env, 1)
  Process(env, exp_source, 1.0, server, 1.1)
  run(env, n)
end

function benchmark(n::Float64, times::Int)
  results = Float64[]
  test_mm1(n)
  for i = 1:times
    push!(results, @elapsed test_mm1(n))
  end
  return(results)
end
")

It can be noted that in both cases there is no monitoring involved, because either nor provide automatic monitoring as does. Furthermore, the resulting code for is more concise and expressive than the equivalent ones for and , which are very similar.

We obtain the reference benchmark with and 20 replicas for both packages as follows:

n <- 1e4L
times <- 20
ref <- data.frame(
  SimPy = rPython::python.call("benchmark", n, times),
  SimJulia = rjulia::j2r(paste0("benchmark(", n, ".0, ", times, ")")))

As a matter of fact, we also tested a small DES skeleton in pure provided in (Matloff 2011, 7.8.3 Extended Example: Discrete–Event Simulation in). This code was formalized into an package called , available on GitHub⁵ since 2014. The original code implemented the event queue as an ordered vector which was updated by performing a binary search. Thus, the execution time of this version was two orders of magnitude slower than the other frameworks. The most recent version on GitHub (as of 2017) takes another clever approach though: it supposes that the event vector will be short and approximately ordered; therefore, the event vector is not sorted anymore, and the next event is found using a simple linear search. These assumptions hold for many cases, and particularly for this M/M/1 scenario. As a result, the performance of this model is only  ∼ 2.2 times slower than . Still, it is clear that pure cannot compete with other languages in discrete-event simulation, and is not considered in our comparisons hereafter.

Finally, we set a benchmark for using , again with and 20 replicas for each test. Figure shows the output of this benchmark. is tested both in monitored and in non-monitored mode. The results show that the performance of is equivalent to and . The non-monitored shows a slightly better performance than these frameworks, while the monitored shows a slightly worse performance.

At this point, it is worth highlighting ’s ability to generate arrivals in batches (hence parameter m). To better understand the impact of batched arrival generation, the benchmark was repeated over a range of m values (1, …, 100). The results of the batched arrival generation runs are shown in Figure . This plot depicts the average execution time of the model with (red) and without (blue) monitoring as a function of the generator batch size m. The black dashed line sets the average execution time of the model to serve as a reference.

The performance with corresponds to what has been shown in Figure . But as m increases, performance quickly improves and becomes  ∼ 1.6 to 1.9 times faster than . Surprisingly, there is no additional gain with batches greater than 40-50 arrivals at a time, but there is no penalty either with bigger batches. Therefore, it is always recommended to generate arrivals in big batches whenever possible.

The simulation core provided by is quite fast, as we have demonstrated, but performance is adversely affected by numerous calls to . The practice of calling from is generally strongly discouraged due to the overhead involved. However, in the case of , it not only makes sense, but is even fundamental in order to provide the user with enough flexibility to build all kinds of simulation models. Nevertheless, this cost must be known, and taken into account whenever a higher performance is needed.

To explore the cost of calling from , let us define the following test:

test_simmer <- function(n, delay) {
  test <- trajectory() %>%
    timeout(delay)
  env <- simmer() %>%
    add_generator("test", test, at(1:n)) %>%
    run(Inf)
  arrivals <- get_mon_arrivals(env)
}

This toy example performs a very simple simulation in which arrivals are attached (in one shot, thanks to the convenience function ) to a trajectory at t = 1, 2, ..., n. The trajectory consists of a single activity: a timeout with some configurable that may be a fixed value or a function call. Finally, after the simulation, the monitored data is extracted from the simulation core to . Effectively, this is equivalent to generating a data frame of rows (see the example output in Table ).

As a matter of comparison, the following function produces the very same data using base :

test_R_for <- function(n) {
  name <- character(n)
  start_time <- numeric(n)
  end_time <- numeric(n)
  activity_time <- logical(n)
  finished <- numeric(n)
  for (i in 1:n) {
    name[i] <- paste0("test", i-1)
    start_time[i] <- i
    end_time[i] <- i+1
    activity_time[i] <- 1
    finished[i] <- TRUE
  }
  arrivals <- data.frame(name = name, start_time = start_time,
    end_time = end_time, activity_time = activity_time,
    finished = finished, replication = 1)
}

Note that we are using a loop to mimic the behavior of ’s internals, of how monitoring is made, but we concede the advantage of pre-allocated vectors to . A second base implementation, which builds upon the lapply() function, is implemented as the function:

test_R_lapply <- function(n) {
  as.data.frame(do.call(rbind, lapply(1:n, function(i) {
    list(name = paste0("test", i - 1), start_time = i,
      end_time = i + 1, activity_time = 1, finished = TRUE,
      replication = 1)
  })))
}

The , and functions all produce exactly the same data in a similar manner (cfr. Table ). Now, we want to compare how a delay consisting of a function call instead of a fixed value impacts the performance of , and we use and as yardsticks.

To this end, the package (Mersmann 2018) is used. The benchmark was executed with and 20 replicas for each test. Table shows a summary of the resulting timings. As we can see, is  ∼ 4.4 times faster than -based base and  ∼ 3.6 times faster than -based base on average when we set a fixed delay. On the other hand, if we replace it for a function call, the execution becomes  ∼ 6.5 times slower, or  ∼ 1.5 times slower than -based base . It is indeed a quite good result if we take into account the fact that base pre-allocates memory, and that is doing a lot more internally. But still, these results highlight the overheads involved and encourage the use of fixed values instead of function calls whenever possible.