MM1-GM1.R

#=====================================================
# Simulation model of a sequence of two queues.
# Email: ebertbrea@gmail.com
#=====================================================

rm(list=ls())
library("simmer")

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

library(parallel)
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.1.2

library(simmer.plot)

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

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

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

library(knitr)

## Warning: package 'knitr' was built under R version 4.1.3

env <- simmer()
###################################################################
#                Parameters of the Simulation
####################################################################
OneRepl=FALSE     # TRUE: for running just one simulation replication;
                  # FALSE: it will run NoReplic simulation replication 
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|mu6|mu7|
muRate<-c( 4 , 4 , 4 , 4 , 4,  4,  4)  # Service rate of kth service
#========| Q1| Q2| Q3| Q4| Q5| Q6| Q7|
CapQ <- c(Inf,Inf,Inf,Inf,Inf, Inf, Inf) # Capacity of each kth queue
#========| S1| S2| S3| S4| S5| S6| S7|
NServ <-c(  1,  1,  1,  1,  1,  1,  1) # Number of servers at each kth queue

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


Nq <- integer(7)# c(  0,  0,  0,  0,  0, 0, 0)
N <- integer(7)

#########################################################
#           Time Functions
#########################################################
BetweenArrival<- function() rexp(1,lambda)
ServiceTime<- function(k){rexp(1,muRate[k])}

############################################

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

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

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(5)) %>%
  release(Resource[5], amount=1)

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

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


# A Sequence of two Queues, namely: Q1 & Q2
Network<-trajectory()%>%
  join(Q1)%>%
  join(Q2)


plot(Network)

if(OneRepl){
  ######################################################### 
  #     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_generator("arrival", Network, BetweenArrival) %>%
    run(ElapsedTime)
  
  resources <- get_mon_resources(Network.env)
  plot(resources, metric = "utilization")
  plot(resources, metric = "usage")


QNr.data <-
  get_mon_arrivals(Network.env) %>%
  dplyr::summarise(mean = mean(end_time - start_time))
}


#########################################################
#     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_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 = 41.645, df = 9, p-value = 1.323e-11
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.9291452 1.0358852
## sample estimates:
## mean of x 
## 0.9825152

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

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8371  0.9451  0.9839  0.9825  1.0231  1.0873

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=".")