QueueMM1WithReplications.R

#=================================================
# This is a simulation model of a queue M/M/C/S, 
# which has been coded by Ebert Brea 
# Emails: ebertbrea@gmail.com, ebert.brea@ucv.ve
#=================================================
rm(list=ls())
library("simmer")

## Warning: package 'simmer' was built under R version 3.5.1

library(parallel)
library(simmer.plot)

## Warning: package 'simmer.plot' was built under R version 3.5.1

## Loading required package: ggplot2

## 
## Attaching package: 'simmer.plot'

## The following objects are masked from 'package:simmer':
## 
##     get_mon_arrivals, get_mon_attributes, get_mon_resources

###################################################
#           Parameters of the Simulation
###################################################
NoReplic <- 10    # Number of replications (Samples)
ElapsedTime <-200 # Elapsed time of simulation 
set.seed(1234)
###################################################
#       Parameters of the Queue M/M/C/S
###################################################
lambda <- 2 # Arrival rate
mu <- 4     # Service rate
QC <- 1     # Capacity of the queue
QS <- Inf   # Size of the queue

############################################
#           Time Functions
############################################
BetweenArrival<- function() rexp(1,lambda)
ServiceTime <- function() rexp(1,mu)
############################################

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

mm1.env <- simmer() %>%
add_resource("mm1.resource",capacity=QC,queue_size=QS) %>%
add_generator("arrival", mm1.traj, BetweenArrival)
  
##############################################
#          Replications  
##############################################

mm1.envs <- mclapply(1:NoReplic, function(i) {
    simmer("Queue_M/M/C/S") %>%
    add_resource("mm1.resource", capacity=QC, queue_size=QS) %>%
    add_generator("arrival", mm1.traj, BetweenArrival) %>%
    run(ElapsedTime) %>%
      wrap()
  })  
  
arrivals <- get_mon_arrivals(mm1.envs)
plot(arrivals, metric = "flow_time")

## Warning: package 'bindrcpp' was built under R version 3.5.1

## `geom_smooth()` using method = 'gam'

plot(arrivals, metric = "activity_time")

## `geom_smooth()` using method = 'gam'

plot(arrivals, metric = "waiting_time")

## `geom_smooth()` using method = 'gam'

resources <- get_mon_resources(mm1.envs)
plot(resources, metric = "utilization")

plot(resources, metric = "usage")

mm1.data <-
  get_mon_arrivals(mm1.envs) %>%
  dplyr::group_by(replication) %>%
  dplyr::summarise(mean = mean(end_time - start_time))
t.test(mm1.data[["mean"]])

## 
##  One Sample t-test
## 
## data:  mm1.data[["mean"]]
## t = 21.724, df = 9, p-value = 4.374e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.4348889 0.5359877
## sample estimates:
## mean of x 
## 0.4854383

summary(mm1.data[["mean"]])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4147  0.4411  0.4585  0.4854  0.5218  0.6487

hist(mm1.data[["mean"]],
     include.lowest = TRUE, right = TRUE,
     breaks="Sturges", col = rgb(0, 1, 0),
     main = paste("Histogram of Mean"), 
     xlab = "Mean",
     density = NULL)

plot.ecdf(mm1.data[["mean"]],
          main= "Empirical cumulative distribution",
          ylab="CDF",xlab="Mean",adj=0.5,
          verticals=TRUE,col.01line=2,cex=0,pch=".")

#Plot of resources. Two metrics available:
#  1) the usage of a resource over the simulation time frame.
#  2) the utilization of specified resources in the simulation.
#Plot of arrivals. Three metrics available:
#  1) activity_time.
#  2) waiting_time.
#  3) flow_time.
#Plot of attributes.