7 min read

Discrete Event Simulation - Part 2

#prepare the libraries
library(tidyverse)
library(simmer)
library(simmer.plot)
library(intervals)

Simmer for R

R is my language of choice. Sorry, Python, I love you too, but life is too short to have two favorite languages.

In my short experiment with Simpy and Simmer, Simmer shines in one area: monitoring. In Simpy a list can be past to the simulation, and it is flexible to choose what to code.

With Simmer, three commands get_mon_arrivals(), get_mon_attributes(), get_mon_resources() can generate the story book of the simulation. The story is in data.frame format, and can be easily processed in R.

Let’s use the Machie Shop as an Example.

Instead of print out the events in user interface with log_(), we make the simulation invisible, and add an attribute to log the state (working, fail).

machine <- function(mttf, repair_time){
      machie <- trajectory() %>% 
      # machine_state: 1-working, 0-fail
      set_attribute("machine_state", 1) %>% 
        timeout(function() rexp(1, 1/mttf)) %>% 
        set_attribute("machine_state", -1, mod="+") %>%
        seize("repairman") %>% 
        timeout(repair_time) %>% 
        release("repairman") %>% 
        set_attribute("machine_state", 1, mod="+") %>% 
        rollback(6)
  } 
plot(machine(10,1))

Make it a function.

simmer_relia <- function(sim_time, mon_ = "attributes") {
  env <- simmer() %>% 
  # use resource to track system state
  # all items in series will require same "repair" resource. 
  # set capacity to inf, assuming no queue for repair
  # if repair resource busy -> system down
  # for n-oo-k items, if resource server > k-n ->system down 
  add_resource("repairman", capacity = Inf) %>% 
  add_generator("machine_a", machine(100, 10), at(0), mon = 2) %>% 
  add_generator("machine_b", machine(100, 5), at(0), mon = 2) %>% 
  run(sim_time) %>% 
  invisible
  if (mon_=="resources"){
    sim_log <- get_mon_resources(env)
  } else if (mon_=="attributes") {
    sim_log <- get_mon_attributes(env)
  }
  sim_log
}

Loop for replications

N_RUN <-  10
SIM_END = 500
ptm <- proc.time()
machine_state_log <- 1:N_RUN %>% 
  map_df(function(x) cbind(simmer_relia(SIM_END), tibble(i_run = x)))
simmer_time <- proc.time() - ptm
tail(machine_state_log)
##         time       name           key value replication i_run
## 197 359.7661 machine_a0 machine_state     0           1    10
## 198 369.7661 machine_a0 machine_state     1           1    10
## 199 399.4862 machine_b0 machine_state     0           1    10
## 200 404.4862 machine_b0 machine_state     1           1    10
## 201 413.9976 machine_b0 machine_state     0           1    10
## 202 418.9976 machine_b0 machine_state     1           1    10

Test

simu_attributes <- simmer_relia(500)
print(simu_attributes)
##          time       name           key value replication
## 1    0.000000 machine_a0 machine_state     1           1
## 2    0.000000 machine_b0 machine_state     1           1
## 3    6.565448 machine_a0 machine_state     0           1
## 4   16.565448 machine_a0 machine_state     1           1
## 5   53.187506 machine_a0 machine_state     0           1
## 6   63.187506 machine_a0 machine_state     1           1
## 7   87.176643 machine_b0 machine_state     0           1
## 8   92.176643 machine_b0 machine_state     1           1
## 9  129.372387 machine_a0 machine_state     0           1
## 10 139.372387 machine_a0 machine_state     1           1
## 11 165.108210 machine_b0 machine_state     0           1
## 12 170.108210 machine_b0 machine_state     1           1
## 13 213.701007 machine_a0 machine_state     0           1
## 14 223.701007 machine_a0 machine_state     1           1
## 15 255.782351 machine_a0 machine_state     0           1
## 16 265.782351 machine_a0 machine_state     1           1
## 17 266.917983 machine_b0 machine_state     0           1
## 18 271.917983 machine_b0 machine_state     1           1
## 19 447.143880 machine_a0 machine_state     0           1
## 20 457.143880 machine_a0 machine_state     1           1
simu_resources <- simmer_relia(500, mon = "resources")
plot(simu_resources, metric = "usage", "repairman", items = "server", steps = TRUE)

Monitoring: resources vs attributes

With simmer, by adding user-defined attributes, we can log the simualtion flexibly. In the context of the system reliability, since only one “repariman” resource is requested when ANY of the machines fail, the “repairman” resource can be used directly for simple system reliability metrics. If the it is series system, whenever the repairman is busy, the system is down. Or it is n-out-of-k, then >= n simultaneous “repairs” implies system down.

The attribute logs the state of all the “machines”, it is more flexible to use the attributes. We are interest in the event sequence of the system down time, the Intervals pacakge, which provides a collections of operation on intervals over the real number line (R), is helpful the determine the system state base on the component state.

In the example above, we have the event sequence of the machine_a and machine_b, so for a parellel system, the intersection of downtime interval is the system down time. and for a series system, union.

Convert event sequences to state intervals

machine_state_interval <- function(machine_state_log){
  machine_state_log %>%
    group_by(i_run, name) %>% 
    arrange(time) %>% 
    # end_time is lead(start_time) (remove leading NA), last one is sim_end
    mutate(end_time = c(head(lead(time), -1), SIM_END)) %>% 
    ungroup() %>% 
    dplyr::select(name, value, start_time = time, end_time, i_run) %>% 
    arrange(i_run, name, start_time)
}

state_interval <- machine_state_interval(machine_state_log)

machine_down_interval <- function(state_interval, i, machine_name){
  fail_log <- state_interval %>% 
    filter(i_run == i, name == machine_name) %>% 
    filter(value == 0)
  Intervals(matrix(c(fail_log$start_time,fail_log$end_time), ncol = 2))
} 

parallel -> intersection of the interval

parallel_down_log <- tibble(start_time = as.numeric(), end_time = as.numeric(), i_run = as.integer())

for (i in 1:N_RUN){
  parallel_down_interval <- interval_intersection(
    machine_down_interval(state_interval, i, "machine_a0"), 
    machine_down_interval(state_interval, i, "machine_b0") )
  if(length(parallel_down_interval) > 0){
    print(i)
    print(parallel_down_interval)
    parallel_down_log <- rbind(parallel_down_log, 
                            cbind(as.matrix(parallel_down_interval), tibble(i_run = i)))
  }
}
## [1] 2
## Object of class Intervals
## 1 interval over R:
## [308.263006831024, 313.263006831024]
## [1] 3
## Object of class Intervals
## 1 interval over R:
## [464.573652200062, 469.573652200062]
## [1] 6
## Object of class Intervals
## 2 intervals over R:
## [122.675988327856, 124.459376529469]
## [354.835283285587, 358.009080867873]
## [1] 9
## Object of class Intervals
## 1 interval over R:
## [101.486453557692, 106.486453557692]
## [1] 10
## Object of class Intervals
## 1 interval over R:
## [68.2142139412463, 71.3290806114674]
names(parallel_down_log) <- c("start", "end", "i_run")
print(parallel_down_log)
##       start       end i_run
## 1 308.26301 313.26301     2
## 2 464.57365 469.57365     3
## 3 122.67599 124.45938     6
## 4 354.83528 358.00908     6
## 5 101.48645 106.48645     9
## 6  68.21421  71.32908    10

Series -> union of the interval

series_down_log <- tibble(start_time = as.numeric(), end_time = as.numeric(), i_run = as.integer())

for (i in 1:N_RUN){
  series_down_interval <- interval_union(
    machine_down_interval(state_interval, i, "machine_a0"), 
    machine_down_interval(state_interval, i, "machine_b0") )
  if(length(series_down_interval) > 0){
    series_down_log <- rbind(series_down_log, 
                             cbind(as.matrix(series_down_interval), tibble(i_run = i)))
  }
}
names(series_down_log) <- c("start", "end", "i_run")
head(series_down_log,20)
##        start       end i_run
## 1   56.64423  66.64423     1
## 2  109.02726 114.02726     1
## 3  155.79494 165.79494     1
## 4  211.19745 216.19745     1
## 5  237.11120 242.11120     1
## 6  255.50268 260.50268     1
## 7  282.21540 287.21540     1
## 8  292.21647 302.21647     1
## 9  388.29990 393.29990     1
## 10 418.07750 428.07750     1
## 11 433.47568 438.47568     1
## 12 461.82842 471.82842     1
## 13 171.84651 176.84651     2
## 14 197.45198 202.45198     2
## 15 274.54367 284.54367     2
## 16 308.23356 318.23356     2
## 17 327.31578 332.31578     2
## 18 372.34737 377.34737     2
## 19 484.74549 489.74549     2
## 20  14.21094  19.21094     3