HW5: First JaamSim Models:

library(tidyverse)
library(kableExtra)
library(latex2exp)

4.10.2:

Create a model like 4-1 except use an arrival rate, \(\lambda\), of 120 entities per hour and service rate, \(\mu\), of 190 entities per hour. Model duration is 100 hours. Calculate exact values for the steady state time entities spend in the system and the expected number of entities processed in 100 hours.

Answer:

Since \(\lambda = 120\) entities per hour or \(\lambda = 120/60 = 2\) would be an exponential distribution with mean \(\frac{1}{\lambda}=0.5\) and \(\mu=190\) served per hour which would be an exponential distribution with mean \(\frac{1}{\mu}=\frac{6}{19}\).

To provide exact values, I’ll reuse my_MMC_func from Discussion Week 4.

my_MMC_func <- function(t_inter, t_svc, c, prob = 1, nm){
  lambda <- (1/t_inter)*prob                         # inter-arrival times
  mu <- 1/t_svc                                      # svc rate
  rho <- lambda/(c*mu)                               # utilization per enti
  i <- (0:(c-1))                                     # the index to sum over
  
  # prob system is empty i.e. p(0)
  p_empty <- 1/(((c*rho)^c)/(factorial(c)*(1-rho)) + sum(((c*rho)^i)/factorial(i)))
  
  Lq <- (rho*((c*rho)^c)*p_empty)/(factorial(c)*(1-rho)^2)#steady-st avg enti in q
  Wq <- Lq/lambda                                         # avg time in queue
  L <- Lq + rho*c                                         # steady-state avg of enti
  W <- Wq + t_svc                                         # steady-st avg time in sy
  Lq <- lambda*Wq                                         # steady-st avg enti in q
  
  a_df <- data.frame(key = c("rho", "L", "Lq", "W", "Wq"),
                     metric = c("Utilization $\\rho$",
                                "Number in System $L$",
                                "Number in Queue $L_q$",
                                "Time in System $W$",
                                "Time in Queue $W_q$"),
                     value = c(rho, L, Lq, W, Wq), 
                     stringsAsFactors = F) %>% 
    select(key, metric, value)
  
  return(a_df)
}

my_mm1 <-my_MMC_func(1/2, 6/19, 1, 1, "M/M/1")

results412 <- my_mm1 %>% select(metric, value)

As shown in the below table, the exact steady state value of time spent in-system is \(W=0.541\) minutes. To get the number of entities processed after 100 hours, I first calculate the number entities per hour and then multiplied by 100 hours: \(\frac{1}{2} \times 60 \times 100 =\) 3,000 entities per hour.

Results
metric value
Utilization \(\rho\) 0.632
Number in System \(L\) 1.714
Number in Queue \(L_q\) 1.083
Time in System \(W\) 0.857
Time in Queue \(W_q\) 0.541

4.10.3:

Run the simulation described in problem 4.10.2 for more than 100 replications, view the histogram \(W\).

Answer:

Below, please see a small screen cap of my JaamSim simulation.

JaamSim Model for 4.10.3

JaamSim Model for 4.10.3

Because I am using JaamSim to complete this assignment, I had to create a custom report via edits to the UniteTypeList and RunOutputList keywords of my simulation in order to caputre the data needed (see comments in code chunk below for example). In the function below, read_sim_rep, I read in the .dat file of the JaamSim report, clean up the column header info, and complete the Little’s Law calculations. Further, I store the results of each run in data table that I can then use to plot.

# helper function to get the student t half-tails. 
my_h <- function(x, ci=.95){
  n <- length(x)
  z_t <- qt( 1 - (1-ci)/2, df = n-1)
  z_t * sd(x)/sqrt(n)
}
# Below is a sample of the RunOutputList used to in this sim:  
#{ [Simulation].RunNumber } { [Simulation].SimTime } 
#{ [ServQueue].AverageQueueTime } { [ServQueue].QueueLengthAverage }
#{ [Serv].Utilisation } { [Serv].WorkingTime } { [Serv].NumberProcessed }
#{ [ArrivalDist].CalculatedMean } { [ServDist].CalculatedMean } 
#{ [ArrivalDist].SampleMean } { [ServDist].SampleMean }

# Each RunOutputList needs a list of units UnitTypeList: 
#DimensionlessUnit TimeUnit TimeUnit DimensionlessUnit 
#DimensionlessUnit TimeUnit DimensionlessUnit TimeUnit 
#TimeUnit TimeUnit TimeUnit
# function reads custom report from JaamSim
read_sim_rep <- function(file_name){
  x <- read.delim(file_name, header = F, stringsAsFactors = F, skip = 2) %>% 
    rename(run = V1, sim_t = V2, Wq_t = V3, Lq = V4, rho = V5, svr_t_tot = V6,
           num_svd_tot = V7, lam_t = V8, mu_t = V9, slam = V10, smu = V11) %>%
    mutate(W_t = Wq_t + smu,
           lam_t = 1/lam_t,
           slam = 1/slam,
           smu = 1/smu,
           L = rho / (1 - rho), 
           n = length(Wq_t)) %>%
    select(run, lam_t, mu_t, Wq_t, W_t, Lq, L, rho, slam, smu, n) 
  
  stat_table <- data.frame(key = c("rho", "L", "Lq", "W", "Wq"), 
                           metric = c("Utilization $\\rho$", 
                  "Number in System $L$", 
                  "Number in Queue $L_q$", 
                  "Time in System $W$", 
                  "Time in Queue $W_q$"), stringsAsFactors = F)
    
  x_stats <- x %>% ungroup() %>% 
    select(Wq_t, W_t, Lq, L, rho, slam, smu, n) %>% 
    summarise_all(funs(mean, my_h)) %>%
    t() %>% as.data.frame() %>% 
    rownames_to_column() %>% arrange(rowname) %>% 
    mutate(key = str_extract(rowname, "^[^_]+(?=_)"), 
           measure = str_extract(rowname, "[^_]*$")) %>% 
    ungroup() %>% select(V1, key, measure) %>% 
    group_by(key, measure) %>% 
    spread(measure, V1)
  
  x_stat_table <- stat_table %>% 
    left_join(x_stats, by = "key") %>% 
    select(key, metric, mean, h)
  
  return(list("x_stat_table" = x_stat_table, "x" = x))
}

sim41 <- read_sim_rep("model43.dat")

results41 <- sim41$x_stat_table %>% select(metric, mean, h)

First, I look and confirm that my function is properly calculating all of the metrics properly:

100 Replications of 100 hour runs
metric mean h
Utilization \(\rho\) 0.631 0.002
Number in System \(L\) 1.715 0.012
Number in Queue \(L_q\) 1.082 0.017
Time in System \(W\) 0.856 0.008
Time in Queue \(W_q\) 0.541 0.008

When you think about it, a “SMORE” plot is just a fancy kinda violin plot. Below, I’ve included violin plots of \(L, L_q, W_q,\) and \(W\).

plot_sim41 <- sim41$x %>% select(L, Lq, Wq_t, W_t) %>% gather()

my_order <- c("L", "Lq", "Wq_t", "W_t")
plot_sim41$key <-  ordered(plot_sim41$key, my_order)

my_boxplot <- plot_sim41

  #geom_boxplot(aes(x= key, y=value))
my_boxplot %>% 
  ggplot(aes(y = value, x = key, group = key),  alpha = .2) +
  facet_grid(~key, scales = "free") +
  geom_violin(alpha = .8) +
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        legend.position="none")

But using R allows me to make histograms, too:

plot_sim41 %>% 
  ggplot(aes(value),  alpha = .2) +
  facet_grid(~key, scales = "free") +
  geom_histogram(alpha = .8) +
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        legend.position="none")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

4.10.9:

Answer:

I suspect the equivalent in JaamSim would be to edit the .cfg file directly. While that’s possible, I’m not sure it’s providing the same lesson. In place of the exercise, I’ve included the .cfg file for my Sim in the commented coded chunk below (all commented out as this is a file that the JaamSim.exe procedure reads to produce the model).

# 
# RecordEdits
# 
# Define ColladaModel { Axis  Grid100x100  model-model  person-model  register-model }
# Define DisplayEntity { XY-Grid  XYZ-Axis }
# Define EntityConveyor { GenToServ  ServToSink }
# Define EntityGenerator { Gen }
# Define EntitySink { Sink }
# Define ExponentialDistribution { ArrivalDist  ServDist }
# Define OverlayClock { Clock }
# Define OverlayText { Title }
# Define Queue { ServQueue }
# Define Server { Server1 }
# Define SimEntity { Proto }
# Define View { View1  View2 }
# 
# ArrivalDist UnitType { TimeUnit }
# ServDist UnitType { TimeUnit }
# 
# Simulation UnitTypeList { DimensionlessUnit  TimeUnit  TimeUnit  DimensionlessUnit  DimensionlessUnit  TimeUnit  DimensionlessUnit  TimeUnit  TimeUnit  TimeUnit  TimeUnit }
# 
# Simulation Description { 'Simulation run control inputs' }
# Simulation RunDuration { 6000  min }
# Simulation InitializationDuration {  }
# Simulation PauseCondition {  }
# Simulation GlobalSubstreamSeed { [Simulation].RunNumber }
# Simulation PrintReport { TRUE }
# Simulation RunOutputList { {  [Simulation].RunNumber  }  {  [Simulation].SimTime  }  {  [ServQueue].AverageQueueTime  }  {  [ServQueue].QueueLengthAverage  }  {  [**DeletedEntity**].Utilisation  }  {  [**DeletedEntity**].WorkingTime  }  {  [**DeletedEntity**].NumberProcessed  }  {  [ArrivalDist].CalculatedMean  }  {  [ServDist].CalculatedMean  }  {  [ArrivalDist].SampleMean  }  {  [ServDist].SampleMean  } }
# Simulation TickLength {  }
# Simulation RunIndexDefinitionList { 100 }
# Simulation StartingRunNumber { 1 }
# Simulation EndingRunNumber { 100 }
# Simulation DisplayedUnits { min }
# Simulation SnapToGrid { TRUE }
# Simulation RealTime { FALSE }
# Simulation RealTimeFactor { 200 }
# Simulation PauseTime {  }
# Simulation ShowModelBuilder { TRUE }
# Simulation ShowObjectSelector { TRUE }
# Simulation ShowInputEditor { TRUE }
# Simulation ShowOutputViewer { TRUE }
# Simulation ShowPropertyViewer { FALSE }
# Simulation ShowLogViewer { FALSE }
# 
# ArrivalDist RandomSeed { 1 }
# ArrivalDist MaxValue {  }
# ArrivalDist Mean { .5  min }
# ArrivalDist Position { -2.500000  -0.600000  0.000000  m }
# 
# Axis ColladaFile { <res>/shapes/axis_text.dae }
# 
# Clock Description { 'Simulation date and time (no leap years or leap seconds)' }
# Clock TextHeight { 10 }
# Clock StartingYear { 2014 }
# Clock DateFormat { 'yyyy-MMM-dd HH:mm:ss.SSS' }
# Clock ScreenPosition { 15  15 }
# Clock AlignBottom { TRUE }
# Clock FontColour { gray20 }
# Clock FontStyle { ITALIC }
# 
# Gen NextComponent { GenToServ }
# Gen FirstArrivalTime {  }
# Gen InterArrivalTime { ArrivalDist }
# Gen PrototypeEntity { Proto }
# Gen Position { -2.500000  0.500000  0.000000  m }
# 
# GenToServ NextComponent { ServQueue }
# GenToServ Position { -2.000000  0.400000  0.000000  m }
# GenToServ Points { {  -2.000  0.400  0.000  m  }  {  -1.600  0.100  0.000  m  } }
# 
# Grid100x100 ColladaFile { <res>/shapes/grid100x100.dae }
# 
# Proto Position { -2.500000  1.600000  0.000000  m }
# Proto Alignment { 0.0  0.0  -0.5 }
# Proto DisplayModel { person-model }
# 
# ServDist RandomSeed { 2 }
# ServDist MaxValue {  }
# ServDist Mean { .315789  min }
# ServDist Position { -0.000000  2.900000  0.000000  m }
# 
# ServQueue Position { -1.300000  0.100000  0.000000  m }
# 
# ServToSink NextComponent { Sink }
# ServToSink Position { 0.400000  1.400000  0.000000  m }
# ServToSink Points { {  0.400  1.400  0.000  m  }  {  3.000  0.300  0.000  m  } }
# 
# Server1 NextComponent { ServToSink }
# Server1 WaitQueue { ServQueue }
# Server1 ServiceTime { ServDist }
# Server1 Position { 0.000000  1.100000  0.000000  m }
# Server1 DisplayModel { model-model }
# 
# Sink Position { 3.400000  0.400000  0.000000  m }
# 
# Title Description { 'Title for the simulation model' }
# Title TextHeight { 18 }
# Title Format { Model4.10.3 }
# Title ScreenPosition { 15  15 }
# Title FontColour { 150  23  46 }
# Title FontStyle { BOLD }
# 
# View1 Description { 'Default view window' }
# View1 ViewCenter { -4.496528  3.114731  -0.451659  m }
# View1 ViewPosition { 2.506600  -1.736179  1.105573  m }
# View1 WindowSize { 940  621 }
# View1 ShowWindow { TRUE }
# View1 Lock2D { FALSE }
# View1 SkyboxImage { <res>/images/sky_map_2048x1024.jpg }
# 
# View2 ViewCenter { -2.841356  2.649494  -0.865653  m }
# View2 ViewPosition { -2.841356  2.649494  7.794602  m }
# View2 WindowPosition { 270  89 }
# View2 ShowWindow { FALSE }
# 
# XY-Grid Description { 'Grid for the X-Y plane (100 m x 100 m)' }
# XY-Grid Size { 100  100  m }
# XY-Grid DisplayModel { Grid100x100 }
# XY-Grid Show { TRUE }
# XY-Grid Movable { FALSE }
# 
# XYZ-Axis Description { 'Unit vectors' }
# XYZ-Axis Alignment { -0.4393409  -0.4410096  -0.4394292 }
# XYZ-Axis Size { 1.125000  1.1568242  1.1266404  m }
# XYZ-Axis DisplayModel { Axis }
# XYZ-Axis Show { FALSE }
# XYZ-Axis Movable { FALSE }
# 
# model-model ColladaFile { model.dae }
# 
# person-model ColladaFile { person.dae }
# 
# register-model ColladaFile { register.dae }

4.10.12:

Answer:

Below, please find a screen capture of my model working in 3D with fancy graphics mapped over them in the simulation:

A screen cap of my 3D animation in process

A screen cap of my 3D animation in process

The results of this simulation can be found under the answer to 4.10.3, above.