MM1k.R

#=====================================================
# Simulation model of an M/M/n/c/Inf/FIFO queue.
# 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 <- 100     # 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|
muRate<-c( 4 , 4 )  # Service rate of kth service
#========| Q1| Q2|
CapQ <- c(Inf,Inf) # Capacity of each kth queue
#========| S1| S2|
NServ <-c(  1,  1) # Number of servers at each kth queue

Resource <-c("Q1Res","Q2Res")
Queue<-c("Queue1","Queue2")


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

#########################################################
#           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)


# A Queue, namely: Q1 
Network<-trajectory()%>%
  join(Q1)


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_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_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"]],mu=0.5,alternative = c("two.sided"))

## 
##  One Sample t-test
## 
## data:  QN.data[["mean"]]
## t = -0.81485, df = 99, p-value = 0.4171
## alternative hypothesis: true mean is not equal to 0.5
## 95 percent confidence interval:
##  0.4847334 0.5063780
## sample estimates:
## mean of x 
## 0.4955557

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

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3426  0.4643  0.4930  0.4956  0.5206  0.6740

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