#=====================================================
# This is a simulation model of a three-queue network,
# which has been coded under R by Ebert Brea.
# Emails: ebertbrea@gmail.com, ebert.brea@ucv.ve
#=====================================================
rm(list=ls())
library("simmer")
library(parallel)
library(simmer.plot)
## 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 <- 100 # Number of replications (Samples)
ElapsedTime <-1000 # Elapsed time of simulation
set.seed(1234)
###################################################
# Parameters of the Queue M/M/C/S
###################################################
lambda <- 2 # Arrival rate
mu1 <- 4 # Service rate of 1
QC1 <- Inf # Capacity of the queue 1
NServ1 <- 1 # Number of service 1
mu2 <- 4 # Service rate of 2
QC2 <- Inf # Capacity of the queue 2
NServ2 <- 1 # Number of service 2
mu3 <- 4 # Service rate of 3
QC3 <- Inf # Capacity of the queue 3
NServ3 <- 1 # Number of service 3
############################################
# Time Functions
############################################
BetweenArrival<- function() rexp(1,lambda)
ServiceTime1 <- function() rexp(1,mu1)
ServiceTime2 <- function() rexp(1,mu2)
ServiceTime3 <- function() rexp(1,mu3)
############################################
Q2.traj <- trajectory() %>%
seize("Q2.resource", amount=1) %>%
timeout(ServiceTime2) %>%
release("Q2.resource", amount=1)
Q3.traj <- trajectory() %>%
seize("Q3.resource", amount=1) %>%
timeout(ServiceTime3) %>%
release("Q3.resource", amount=1)
Q1.traj <- trajectory() %>%
seize("Q1.resource", amount=1) %>%
timeout(ServiceTime1) %>%
release("Q1.resource", amount=1) %>%
branch(option = function() {round(ifelse(runif(1,0,100)<50,1,2))},
continue=c(FALSE, TRUE),
join(Q2.traj),
join(Q3.traj)%>%rollback(amount=7, times = 1,check = NULL)
)
plot(Q1.traj)
############################################
# For running "NoReplic" replications
############################################
envs <- mclapply(1:NoReplic, function(i) {
simmer("QueueNetwork") %>%
add_resource("Q1.resource",capacity=NServ1,queue_size=QC1) %>%
add_resource("Q2.resource",capacity=NServ2,queue_size=QC2) %>%
add_resource("Q3.resource",capacity=NServ3,queue_size=QC3) %>%
add_generator("arrival", Q1.traj, BetweenArrival) %>%
run(ElapsedTime) %>%
wrap()
})
arrivals <- get_mon_arrivals(envs)
plot(arrivals, metric = "flow_time")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
plot(arrivals, metric = "activity_time")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
plot(arrivals, metric = "waiting_time")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
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 = 116.54, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 1.478295 1.529507
## sample estimates:
## mean of x
## 1.503901
summary(QN.data[["mean"]])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.196 1.421 1.505 1.504 1.563 1.891
hist(QN.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(QN.data[["mean"]],
main= "Empirical cumulative distribution",
ylab="CDF",xlab="Mean",adj=0.5,
verticals=TRUE,col.01line=2,cex=0,pch=".")
