4

Is there a way to use a (customized) routing engine together with the simmer package for discrete-event simulation? (or an alternative package)

Context: I'm running dicrete-event simulations (DES) with R. Till now all my simulations are built without using one of the R packages designed for DES. Since my code gets bigger and bigger (and performance worse) I'm thinking about switching to one of the R packages designed for DES.

For some portions of my code I see how I could switch it to simmer. But till now I couldn't figure out how to use a routing-logic together with resource dispatching.

Example: The following minimal example shows what kind of functionality I need (and couldn't figure out how to build with simmer).

Generate some data, events (jobs) and resources

set.seed(1)

events <- data.frame(
  id = 1:3L,
  t = sort(trunc(rexp(3) * 100)),
  position = runif(3),
  resource = NA,
  worktime = NA
)

resources <- data.frame(
  id = 1:2L,
  position = c(0.2, 0.8),
  t_free = 0
)

Simplified version of a routing logic: calculate the route based on position of event and resources. (For the example just points on a 1-D space between 0 and 1, in the real example a customized version of OSRM algorithm together with historical data..)

waytime <- function(events, resources, i) {
  trunc(abs(events$position[i] - resources$position[resources$id == events$resource[i]]) * 100)
}

Two versions of the simulation. sim just takes the first available resource with no thinking about the waytime. sim_nearest calculates waytimes for all free resources and dispatches to the closest one. sim_nearest is what I want in my real examples and don't know how to build using simmer.

sim <- function(events, resources) {
  for (i in 1:nrow(events)) {
    # Default dispatching: Use the first free vehicle
    events$resource[i] <- resources$id[resources$t_free <= events$t[i]][1]
    # Simulate event
    events$worktime[i] <- waytime(events, resources, i)
    resources$t_free[events$resource[i]] <- events$t[i] + events$worktime[i]
  }
  return(list(events = events, resources = resources))
}

sim_use_nearest <- function(events, resources) {
  for (i in 1:nrow(events)) {
    # Dispatching by position: Use the nearest free resource
    ids_free <- resources$id[resources$t_free <= events$t[i]]
    events$resource[i] <- resources$id[which.min(abs(resources$position[ids_free] - events$position[i]))]
    # Simulate event
    events$worktime[i] <- waytime(events, resources, i)
    resources$t_free[events$resource[i]] <- events$t[i] + events$worktime[i]
  }
  return(list(events = events, resources = resources))
}

Simulate the two alternatives:

res <- sim(events, resources)
res_use_nearest <- sim_use_nearest(events, resources)

See the differences:

res$events
# id   t  position resource worktime
#  1  14 0.9082078        1       70
#  2  75 0.2016819        2       59
#  3 118 0.8983897        1       69
res$resources
# id position t_free
#  1      0.2    187
#  2      0.8    134
res_use_nearest$events
# id   t  position resource worktime
#  1  14 0.9082078        2       10
#  2  75 0.2016819        1        0
#  3 118 0.8983897        2        9
res_use_nearest$resources
# id position t_free
#  1      0.2     75
#  2      0.8    127

Is it possible to generate the same results with simmer (or another R DES package)?

symbolrush
  • 7,123
  • 1
  • 39
  • 67

3 Answers3

3

Following you'll find a possible solution for your minimal example with the simmer package.

First we chose the alternative to simulate which is later used in set_attribute:

sim_first_available <- T
sim_use_nearest <- F

Generate the events and resources data as before. 

set.seed(1)

events <- data.frame(
  id = 1:3L,
  t = sort(trunc(rexp(3) * 100)),
  position = runif(3),
  resource = NA,
  worktime = NA
)

resources <- data.frame(
  id = 1:2L,
  position = c(0.2, 0.8),
  t_free = 0
)

Start simmer with the trajectory sim.

library(simmer)

sim <- trajectory() %>%

Then set t_free as a global attribute. At the first arrival (t = 14) you can use t_free from the resource data to initialize. At later arrivals use get_global to get the current t_free of the specific resource.

  set_global(paste0("t_free_res_", resources$id), function() {
    if (now(env) == 14) {return(resources$t_free) # Initialize parameters when first event arrives
    } else {
      get_global(env, paste0("t_free_res_", resources$id))
    }}) %>%

Now define the attributes for this event:

Based on the current simulation time choose the event_position from the data frame events.

  set_attribute(c("event_position","my_resource", "timeout"), function() {
    t <- now(env)
    event_position <- events$position[events$t == t]

my_resource is chosen acc. to the alternative you want to simulate. 

    t_free <- get_global(env, paste0("t_free_res_", resources$id))
    if (sim_first_available & !sim_use_nearest) {
      my_resource <- resources$id[t_free <= now(env)][1]
    } else if (!sim_first_available & sim_use_nearest){
      ids_free <- resources$id[t_free <= now(env)]
      my_resource <- resources$id[which.min(abs(resources$position[ids_free] - event_position))]
    }

Based on the resource_pos calculate the timeout for that resource and return the attributes:

resource_pos <- resources$position[resources$id == my_resource]
        timeout <- trunc(abs(event_position - resource_pos)*100)

        return(c(event_position, my_resource, timeout))
      }) %>%

Select the defined resource and seize it:

  select(resources = function() paste0("res_", get_attribute(env, "my_resource"))) %>%
  seize_selected(amount = 1) %>%

Now overwrite t_free of that resource by adding the timeout to the current simulation time.

  set_global(function() {
    paste0("t_free_res_", get_attribute(env, "my_resource"))
  }, function() {
    return(now(env) + get_attribute(env, "timeout"))
  }) %>%

Set the calculated timeout to the resource and release it again.

  timeout(function() get_attribute(env, "timeout")) %>% 
  release_selected(amount = 1)

Finally generate the events for the trajectory sim at the defined time intervals in events, add the resources and run the simulation.

env <- simmer()  %>%
  add_generator("event_", sim, at(events$t), mon = 2) %>%
  add_resource("res_1", capacity = 1) %>%
  add_resource("res_2", capacity = 1)

env %>% run()

print(get_mon_attributes(env))
print(get_mon_arrivals(env))
print(get_mon_resources(env))

Hope this helps.

Samy
  • 31
  • 2
3

Samy's approach is fine, but I would take a slightly different one (note that this isn't tested, because I didn't write the necessary routing_logic function):

library(simmer)

env <- simmer()

t <- trajectory() %>%
  seize("available_resources") %>%
  set_attribute(c("res_id", "delay"), routing_logic) %>%
  select(function() paste0("res_", get_attribute(env, "res_id"))) %>%
  seize_selected() %>%
  timeout_from_attribute("delay") %>%
  release_selected() %>%
  release("available_resources")

Note that "available_resources" (which must be a resource with a capacity equal to the number of resources you have) is like a token. Once seized, it means that there's some resource available. Otherwise, events just sit there and wait.

routing_logic() must be a function that selects a "res_id" based on some policy (e.g., first available or nearest), computes the delay and returns both values, which are stored as attributes. In that function, you may use get_capacity() to know the status of each resource without having to set t_free. You may also retrieve the position attribute for that event, which will be set automatically as follows:

set.seed(1)

events <- data.frame(
  t = sort(trunc(rexp(3) * 100)),
  position = runif(3)
)

resources <- data.frame(
  id = 1:2L,
  position = c(0.2, 0.8)
)

env %>% 
  add_dataframe("event_", t, events, mon=2, col_time="t", time="absolute") %>%
  add_resource("available_resources", capacity=nrow(resources))

for (id in resources$id) env %>%
  add_resource(paste0("res_", id), capacity=1, queue_size=0)

As you can see, I have directly connected the events data frame to the trajectory (you don't need the resource and worktime anymore; the former will be stored as the res_id attribute, and the latter will be automatically monitored by simmer and retrieved with get_mon_arrivals()). We specify that t is the time column, and the other one, position will be added to each event as an attribute, as I said before.

With this setup, you just need to redefine routing_logic() to achieve different policies and different results.

Iñaki Úcar
  • 935
  • 1
  • 5
  • 11
2

Iñaki's approch is very useful, as it uses features of the newest simmer version. Out of interest I completed his example with a routing logic and - as expected - the results are the same. Thanks for your Input Iñaki.

library(simmer)

env <- simmer()

t <- trajectory() %>%
  seize("available_resources") %>%
  set_attribute(c("res_id", "delay"), function() {
    # find available resources
    capacities <- numeric(nrow(resources))
    for (i in 1:length(capacities)) {
      capacities[i] <- get_server_count(env, paste0("res_", resources$id[i]))
    }
    available <- ifelse(capacities == 0, T, F)
    index_available <- which(available)
    # calculate the delay for available resources
    event_position <- get_attribute(env, "position")
    delay <- trunc(abs(event_position - resources$position[available])*100)
    # take the nearest available resource. 
    index <- index_available[which.min(delay)]
    return(c(index,min(delay)))
  }) %>%
  select(function() paste0("res_", get_attribute(env, "res_id"))) %>%
  seize_selected() %>%
  timeout_from_attribute("delay") %>%
  release_selected() %>%
  release("available_resources")
# --------------------------------------------------------------------
set.seed(1)

events <- data.frame(
  t = sort(trunc(rexp(3) * 100)),
  position = runif(3)
)

resources <- data.frame(
  id = 1:2L,
  position = c(0.2, 0.8)
)

env %>% 
  add_dataframe("event_", t, events, mon=2, col_time="t", time="absolute") %>%
  add_resource("available_resources", capacity=nrow(resources))
for (id in resources$id) env %>%
  add_resource(paste0("res_", id), capacity=1, queue_size=0)

env %>% run()
# --------------------------------------------------------------------
library(simmer.plot)
print(plot(get_mon_resources(env), metric = "usage", c("available_resources", "res_1", "res_2"), items = "server", steps = TRUE))
MichiSmith
  • 317
  • 1
  • 7