QueueArray2.R

#=====================================================
# Simulation model of a five-queue network, wherein has
# been included a bifurcation of path and a rollback. 
# 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

env <- simmer()
###################################################
#           Parameters of the Simulation
###################################################
NoReplic <- 10    # Number of replications (Samples)
ElapsedTime <-400 # Elapsed time of simulation 
set.seed(1234)
#########################################################
#       Parameters of the Queues
#########################################################
lambda <- 2                    # Arrival rate to the network
#========|mu1|mu2|mu3|mu4|mu5|
muRate<-c( 4 , 4 , 4 , 4 , 4)  # Service rate of kth service
#========| Q1| Q2| Q3| Q4| Q5|
CapQ <- c(Inf,Inf,Inf,Inf,Inf) # Capacity of each kth queue
#========| S1| S2| S3| S4| S5|
NServ <-c(  1,  1,  1,  1,  1) # Number of servers at each kth queue

Resource <-c("Q1Res","Q2Res","Q3Res","Q4Res","Q5Res")
Queue<-c("Queue1","Queue2","Queue3","Queue4","Queue5")

#########################################################
#           Time Functions
#########################################################
BetweenArrival<- function() rexp(1,lambda)
ServiceTime<- function(k){rexp(1,muRate[k])}
############################################
Q3<-trajectory() %>%
  seize(Resource[3], amount=1) %>%
  timeout(function() ServiceTime(3)) %>%
  release(Resource[3], amount=1)

Q4<-trajectory() %>%
  seize(Resource[4], amount=1) %>%
  timeout(function() ServiceTime(4)) %>%
  release(Resource[4], amount=1)

Q5<-trajectory() %>%
  seize(Resource[5], amount=1) %>%
  timeout(function() ServiceTime(i)) %>%
  release(Resource[5], amount=1)

Network<-trajectory()

for(i in 1:2){
  Network<-Network%>%
  seize(Resource[i], amount=1) %>%
  timeout(function() ServiceTime(i)) %>%
  release(Resource[i], amount=1)}

  Network<-Network%>%
  branch(option = function() {round(ifelse(runif(1,0,100)<50,1,2))},
         continue=c(TRUE, TRUE), 
         join(Q3),
         join(Q4)%>%
         rollback(amount=7, times = 1,check = NULL)%>%join(Q5)
         )

plot(Network)

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

#########################################################
#     For running just one simulation sample
#########################################################
Network.env<-simmer("QueueNetwork") %>%
  add_resource(Resource[1],capacity=NServ[1],queue_size=CapQ[1]) %>%
  add_resource(Resource[2],capacity=NServ[2],queue_size=CapQ[2]) %>%
  add_resource(Resource[3],capacity=NServ[3],queue_size=CapQ[3]) %>%
  add_resource(Resource[4],capacity=NServ[4],queue_size=CapQ[4]) %>%
  add_resource(Resource[5],capacity=NServ[5],queue_size=CapQ[5]) %>%
  add_generator("arrival", Network, BetweenArrival) %>%
  run(ElapsedTime)
  
resources <- get_mon_resources(Network.env)
plot(resources, metric = "utilization")

plot(resources, metric = "usage")

#########################################################
#     For running "NoReplic" simulation replications
#########################################################

envs <- mclapply(1:NoReplic, function(i) {
  simmer("QueueNetwork") %>%
    add_resource(Resource[1],capacity=NServ[1],queue_size=CapQ[1]) %>%
    add_resource(Resource[2],capacity=NServ[2],queue_size=CapQ[2]) %>%
    add_resource(Resource[3],capacity=NServ[3],queue_size=CapQ[3]) %>%
    add_resource(Resource[4],capacity=NServ[4],queue_size=CapQ[4]) %>%
    add_resource(Resource[5],capacity=NServ[5],queue_size=CapQ[5]) %>%
    add_generator("arrival", Network, BetweenArrival) %>%
    run(ElapsedTime) %>%
    wrap()
})

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

plot(resources, metric = "usage")

QN.data <-
  get_mon_arrivals(envs) %>%
  dplyr::group_by(replication) %>%
  dplyr::summarise(mean = mean(end_time - start_time))

t.test(QN.data[["mean"]])

## 
##  One Sample t-test
## 
## data:  QN.data[["mean"]]
## t = 20.872, df = 9, p-value = 6.231e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.945974 2.419079
## sample estimates:
## mean of x 
##  2.182527

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

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.821   2.013   2.121   2.183   2.271   2.959

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

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