#.rs.restartR()
#===============================================================================
# Modelled and coded by Ebert Brea
#===============================================================================
# This simulation model shows a queue M/M/1/k with dropped entities (customers);
# for what is taken into account in the statistical report, the entities
# that leave the queue when the queue reaches its maximum capacity k.
#===============================================================================
cat("\f")
rm(list=ls())
library(xtable)
library(latex2exp)
## Warning: package 'latex2exp' was built under R version 4.1.2
library(stargazer)
## Warning: package 'stargazer' was built under R version 4.1.2
##
## Please cite as:
## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
simulaMM1K <- function(lambda,mu,cap,ElapsedTime,NoSeed){
# ========================================================
# lambda: arrival rate
# mu: service rate
# ElapsedTime: period time simulation
# NoSeed: number of seed
# =======================================================
#lambda=10; mu=lambda/0.8; cap=5; ElapsedTime=400; NoSeed=20
set.seed(NoSeed)
library(simmer)
library(parallel)
env<-simmer()
BetweenArrival<- function(){rexp(1,lambda)}
ServiceTime<- function(){rexp(1,mu)}
CapacityQueue<-cap
path<-function(){ ifelse( get_queue_count(env, "Server") < CapacityQueue ,1, 2) }
#===============================================================================#%
Start<-trajectory("Start")%>%
set_global("Counter", 1, mod="+",init=-1) %>%
set_attribute("My_ID", function() get_global(env, "Counter") )%>%
set_attribute("QueueLength", function() get_queue_count(env, "Server") )%>%
set_attribute("Time_in",function(){now(env)})
Q1<-trajectory("Queue 1") %>%
seize("Server", amount=1) %>%
timeout(ServiceTime) %>%
release("Server", amount=1)
# %>%log_("from Q1")
Q2<-trajectory()%>%
set_global("Leave", 1, mod="+",init=0)
#%>% log_(function() {paste("Leaving Network, because the Queue 1 has", get_queue_count(env, "Server_1"), "customers")})
Network<-trajectory()%>%
join(Start)%>%
# log_(function() {paste("Queue 1 has", get_queue_count(env, "Server_1"), "customer(s) at my arrival")})%>%
branch(option = path,
continue=c(TRUE,TRUE),
join(Q1),
join(Q2))
env<-simmer()
env<-env %>%
add_resource("Server", capacity=1, queue_size = CapacityQueue) %>%
add_generator("Customer", Network, BetweenArrival, mon=2) %>%
run(ElapsedTime)
release_all(Network, "Server")
resources <- get_mon_resources(env)
arrivals <- get_mon_arrivals(env,per_resource=TRUE)
attribites <-get_mon_attributes(env)
Arrival.df<-data.frame(get_mon_arrivals(env))
Attribute.df<-data.frame(get_mon_attributes(env))
Ndata<-length(Attribute.df[,4])
tmp<-numeric(Ndata)
QL.Obs<-1
for(ii in 1:Ndata){
if(Attribute.df$key[ii]=="QueueLength"){tmp[QL.Obs]<-Attribute.df$value[ii];QL.Obs<-QL.Obs+1}
}
QL<-numeric(QL.Obs)
for(i in 1:QL.Obs){QL[i]<-tmp[i]}
DropPerformance <- c( sum(Attribute.df$key=="Leave"),
sum(Arrival.df$finished),
mean(Arrival.df$end_time - Arrival.df$start_time),
mean(QL))
# DropPerformance <- data.frame(attribites)
return(DropPerformance)
}
#==========================================================
# Multireplications
#==========================================================
NoReplications <- 30
L<-10 # Lambda
Rho<-0.8
#==============================================================================
Result <- matrix(data = NA, nrow = NoReplications, ncol = 4, byrow = TRUE)
for(r in 1:NoReplications){
Result[r,] <-simulaMM1K(lambda=L, mu=L/Rho, cap=4, ElapsedTime=400,NoSeed=r)
}
## Warning: package 'simmer' was built under R version 4.1.2
df.Result <- data.frame(Result)
names(df.Result) <- c("Leave","Finish","System.Time","Queue.Length")
summary(df.Result$Leave)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 289.0 342.2 361.0 366.0 392.8 445.0
summary(df.Result$System.Time)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1772 0.1838 0.1856 0.1857 0.1892 0.1933
summary(df.Result$Queue.Length)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.021 1.104 1.159 1.142 1.181 1.241
var(df.Result$Queue.Length)
## [1] 0.002930475
hist(df.Result$Leave,main="Dropped customers",
xlab="Number of customers",
col=blues9,labels=TRUE)
RatioDropped <-df.Result$Leave/df.Result$Finish
df.Result <- data.frame(RatioDropped,Result)
names(df.Result) <- c("Leave/Finish", "Leave","Finish","System.Time","Queue.Length")
summary(RatioDropped)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.07294 0.08540 0.09084 0.09112 0.09761 0.10792
t.test(df.Result$System.Time,mu = 0.186)
##
## One Sample t-test
##
## data: df.Result$System.Time
## t = -0.35767, df = 29, p-value = 0.7232
## alternative hypothesis: true mean is not equal to 0.186
## 95 percent confidence interval:
## 0.1842235 0.1872477
## sample estimates:
## mean of x
## 0.1857356
t.test(df.Result$Queue.Length, mu = 1.14)
##
## One Sample t-test
##
## data: df.Result$Queue.Length
## t = 0.2403, df = 29, p-value = 0.8118
## alternative hypothesis: true mean is not equal to 1.14
## 95 percent confidence interval:
## 1.122161 1.162589
## sample estimates:
## mean of x
## 1.142375
t.test(L*df.Result$System.Time, df.Result$Queue.Length,
paired = TRUE, var.equal = TRUE, conf.level = 0.95,mu=0.7)
##
## Paired t-test
##
## data: L * df.Result$System.Time and df.Result$Queue.Length
## t = 2.2425, df = 29, p-value = 0.03274
## alternative hypothesis: true difference in means is not equal to 0.7
## 95 percent confidence interval:
## 0.7013178 0.7286435
## sample estimates:
## mean of the differences
## 0.7149807
#========================================================
stargazer(df.Result, type = "text",
summary.stat = c("min", "p25", "median", "p75", "max", "mean", "sd"),
title="Performance of a queue M/M/1/k")
##
## Performance of a queue M/M/1/k
## ====================================================================
## Statistic Min Pctl(25) Median Pctl(75) Max Mean St. Dev.
## --------------------------------------------------------------------
## Leave/Finish 0.073 0.085 0.091 0.098 0.108 0.091 0.009
## Leave 289 342.2 361 392.8 445 366.033 39.146
## Finish 3,896 3,963.8 3,995 4,062.5 4,146 4,013.400 64.792
## System.Time 0.177 0.184 0.186 0.189 0.193 0.186 0.004
## Queue.Length 1.021 1.104 1.159 1.181 1.241 1.142 0.054
## --------------------------------------------------------------------
stargazer(df.Result, type = "latex",
summary.stat = c("min", "p25", "median", "p75", "max", "mean", "sd"),
title="Performance of a queue M/M/1/k")
##
## % Table created by stargazer v.5.2.3 by Marek Hlavac, Social Policy Institute. E-mail: marek.hlavac at gmail.com
## % Date and time: dom., nov. 27, 2022 - 11:32:17 a. m.
## \begin{table}[!htbp] \centering
## \caption{Performance of a queue M/M/1/k}
## \label{}
## \begin{tabular}{@{\extracolsep{5pt}}lccccccc}
## \\[-1.8ex]\hline
## \hline \\[-1.8ex]
## Statistic & \multicolumn{1}{c}{Min} & \multicolumn{1}{c}{Pctl(25)} & \multicolumn{1}{c}{Median} & \multicolumn{1}{c}{Pctl(75)} & \multicolumn{1}{c}{Max} & \multicolumn{1}{c}{Mean} & \multicolumn{1}{c}{St. Dev.} \\
## \hline \\[-1.8ex]
## Leave/Finish & 0.073 & 0.085 & 0.091 & 0.098 & 0.108 & 0.091 & 0.009 \\
## Leave & 289 & 342.2 & 361 & 392.8 & 445 & 366.033 & 39.146 \\
## Finish & 3,896 & 3,963.8 & 3,995 & 4,062.5 & 4,146 & 4,013.400 & 64.792 \\
## System.Time & 0.177 & 0.184 & 0.186 & 0.189 & 0.193 & 0.186 & 0.004 \\
## Queue.Length & 1.021 & 1.104 & 1.159 & 1.181 & 1.241 & 1.142 & 0.054 \\
## \hline \\[-1.8ex]
## \end{tabular}
## \end{table}
hist(RatioDropped,
include.lowest = TRUE,
right = TRUE,
breaks = nclass.Sturges,
main=paste("Ratio of dropped customers with","Rho =",Rho ),
xlab="Ratio",
col=rgb(0.5, 0.9, 0),
labels=TRUE,
prob=FALSE)
# =========================================================
# Experiments of dropped Customers
# =========================================================
NoReplications <- 10
mu <- seq(2,5,by=1)
lambda<-seq(2,5,by=1)
k<-4
SimTime<-400
Nr<-3
N1<-length(lambda)
N2<-length(mu)
#========================================================
# For writing a statistical report
#========================================================
LoDropped <-numeric(N1*N2)
LoAttended <-numeric(N1*N2)
LoRatioDropAtt <-numeric(N1*N2)
LoSystemTime <-numeric(N1*N2)
LoQueueLength <-numeric(N1*N2)
mDropped <-numeric(N1*N2)
mAttended <-numeric(N1*N2)
mRatioDropAtt <-numeric(N1*N2)
mSystemTime <-numeric(N1*N2)
mQueueLength <-numeric(N1*N2)
UpDropped <-numeric(N1*N2)
UpAttended <-numeric(N1*N2)
UpRatioDropAtt <-numeric(N1*N2)
UpSystemTime <-numeric(N1*N2)
UpQueueLength <-numeric(N1*N2)
Mu <- numeric(N1*N2)
Lambda <- numeric(N1*N2)
MultiResult = array(data=NA, dim=c(N1*N2,4,NoReplications))
# MultiResult is a 3D matrix, e.i., num [1:N1*N2, 1:3, 1:NoReplications]
n <- 0
for(i in 1:N1){
for(j in 1:N2){
n <- n+1
for(r in 1:NoReplications){
MultiResult[n,,r] <- simulaMM1K(lambda=lambda[i], mu=mu[j], cap=k, ElapsedTime=SimTime, NoSeed=r)
}
Lambda[n] <- lambda[i]
Mu[n] <- mu[j]
}
}
for(n in 1:(N1*N2)){
LoDropped[n] <-t.test(MultiResult[n,1,])$conf.int[1]
LoAttended[n] <-t.test(MultiResult[n,2,])$conf.int[1]
LoRatioDropAtt[n] <- t.test( (MultiResult[n,2,]-MultiResult[n,1,])/MultiResult[n,2,] )$conf.int[1]
LoSystemTime[n] <- t.test(MultiResult[n,3,])$conf.int[1]
LoQueueLength[n] <- t.test(MultiResult[n,4,])$conf.int[1]
mDropped[n] <-mean(MultiResult[n,1,])
mAttended[n] <-mean(MultiResult[n,2,])
mRatioDropAtt[n] <- mean( (MultiResult[n,2,]-MultiResult[n,1,])/MultiResult[n,2,] )
mSystemTime[n] <- mean(MultiResult[n,3,])
mQueueLength[n] <- mean(MultiResult[n,4,])
UpDropped[n] <-t.test(MultiResult[n,1,])$conf.int[2]
UpAttended[n] <-t.test(MultiResult[n,2,])$conf.int[2]
UpRatioDropAtt[n] <- t.test( (MultiResult[n,2,]-MultiResult[n,1,])/MultiResult[n,2,] )$conf.int[2]
UpSystemTime[n] <- t.test(MultiResult[n,3,])$conf.int[2]
UpQueueLength[n] <- t.test(MultiResult[n,4,])$conf.int[2]
}
df.ST <- data.frame(Lambda,Mu,mDropped,mAttended,mRatioDropAtt,mSystemTime,mQueueLength)
print(df.ST)
## Lambda Mu mDropped mAttended mRatioDropAtt mSystemTime mQueueLength
## 1 2 2 120.8 788.5 0.8472907 1.2478667 1.6025268
## 2 2 3 37.2 801.7 0.9538535 0.7199044 0.7901123
## 3 2 4 11.2 805.9 0.9863279 0.4522702 0.3930539
## 4 2 5 3.0 803.2 0.9962711 0.3168875 0.2218722
## 5 3 2 426.0 1197.1 0.6441162 1.1857804 2.5918564
## 6 3 3 185.2 1193.5 0.8450527 0.8150157 1.5997734
## 7 3 4 77.6 1195.3 0.9351679 0.5556477 0.9499252
## 8 3 5 37.2 1208.9 0.9693971 0.4015357 0.6226144
## 9 4 2 815.8 1604.1 0.4916643 1.0134043 3.0935281
## 10 4 3 474.6 1593.4 0.7023515 0.8145411 2.3323181
## 11 4 4 257.9 1593.0 0.8383335 0.6178971 1.6329986
## 12 4 5 137.6 1589.3 0.9135665 0.4625281 1.1153454
## 13 5 2 1187.6 1991.4 0.4038981 0.8676875 3.3395759
## 14 5 3 839.1 1996.9 0.5798331 0.7588603 2.8157925
## 15 5 4 528.8 1987.5 0.7339547 0.6215195 2.1829451
## 16 5 5 324.3 1987.1 0.8370655 0.4938210 1.6319491
print(xtable(df.ST), include.rownames = TRUE)
## % latex table generated in R 4.1.0 by xtable 1.8-4 package
## % Sun Nov 27 11:32:31 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrr}
## \hline
## & Lambda & Mu & mDropped & mAttended & mRatioDropAtt & mSystemTime & mQueueLength \\
## \hline
## 1 & 2.00 & 2.00 & 120.80 & 788.50 & 0.85 & 1.25 & 1.60 \\
## 2 & 2.00 & 3.00 & 37.20 & 801.70 & 0.95 & 0.72 & 0.79 \\
## 3 & 2.00 & 4.00 & 11.20 & 805.90 & 0.99 & 0.45 & 0.39 \\
## 4 & 2.00 & 5.00 & 3.00 & 803.20 & 1.00 & 0.32 & 0.22 \\
## 5 & 3.00 & 2.00 & 426.00 & 1197.10 & 0.64 & 1.19 & 2.59 \\
## 6 & 3.00 & 3.00 & 185.20 & 1193.50 & 0.85 & 0.82 & 1.60 \\
## 7 & 3.00 & 4.00 & 77.60 & 1195.30 & 0.94 & 0.56 & 0.95 \\
## 8 & 3.00 & 5.00 & 37.20 & 1208.90 & 0.97 & 0.40 & 0.62 \\
## 9 & 4.00 & 2.00 & 815.80 & 1604.10 & 0.49 & 1.01 & 3.09 \\
## 10 & 4.00 & 3.00 & 474.60 & 1593.40 & 0.70 & 0.81 & 2.33 \\
## 11 & 4.00 & 4.00 & 257.90 & 1593.00 & 0.84 & 0.62 & 1.63 \\
## 12 & 4.00 & 5.00 & 137.60 & 1589.30 & 0.91 & 0.46 & 1.12 \\
## 13 & 5.00 & 2.00 & 1187.60 & 1991.40 & 0.40 & 0.87 & 3.34 \\
## 14 & 5.00 & 3.00 & 839.10 & 1996.90 & 0.58 & 0.76 & 2.82 \\
## 15 & 5.00 & 4.00 & 528.80 & 1987.50 & 0.73 & 0.62 & 2.18 \\
## 16 & 5.00 & 5.00 & 324.30 & 1987.10 & 0.84 & 0.49 & 1.63 \\
## \hline
## \end{tabular}
## \end{table}
write.csv(x=df.ST, file="dfST.csv")
#========================================================
df.Report <- data.frame(Lambda, Mu,
LoDropped, mDropped, UpDropped,
LoAttended, mAttended, UpAttended,
LoRatioDropAtt, mRatioDropAtt, UpRatioDropAtt,
LoSystemTime, mSystemTime, UpSystemTime,
LoQueueLength, mQueueLength, UpQueueLength)
write.csv(x=df.Report, file="FullReport.csv")
print(df.Report)
## Lambda Mu LoDropped mDropped UpDropped LoAttended mAttended UpAttended
## 1 2 2 105.171639 120.8 136.428361 767.2490 788.5 809.7510
## 2 2 3 29.031138 37.2 45.368862 783.1495 801.7 820.2505
## 3 2 4 6.818703 11.2 15.581297 781.2840 805.9 830.5160
## 4 2 5 1.738229 3.0 4.261771 783.8786 803.2 822.5214
## 5 3 2 413.292483 426.0 438.707517 1183.9767 1197.1 1210.2233
## 6 3 3 163.521027 185.2 206.878973 1172.6433 1193.5 1214.3567
## 7 3 4 65.491881 77.6 89.708119 1174.6108 1195.3 1215.9892
## 8 3 5 28.570794 37.2 45.829206 1181.1852 1208.9 1236.6148
## 9 4 2 786.458156 815.8 845.141844 1574.5234 1604.1 1633.6766
## 10 4 3 452.279030 474.6 496.920970 1563.2130 1593.4 1623.5870
## 11 4 4 229.292779 257.9 286.507221 1564.8202 1593.0 1621.1798
## 12 4 5 123.678796 137.6 151.521204 1562.3833 1589.3 1616.2167
## 13 5 2 1145.893063 1187.6 1229.306937 1954.6007 1991.4 2028.1993
## 14 5 3 815.225101 839.1 862.974899 1967.5937 1996.9 2026.2063
## 15 5 4 502.490595 528.8 555.109405 1969.0257 1987.5 2005.9743
## 16 5 5 289.075098 324.3 359.524902 1957.8695 1987.1 2016.3305
## LoRatioDropAtt mRatioDropAtt UpRatioDropAtt LoSystemTime mSystemTime
## 1 0.8303676 0.8472907 0.8642138 1.2150016 1.2478667
## 2 0.9444996 0.9538535 0.9632073 0.6980260 0.7199044
## 3 0.9814883 0.9863279 0.9911675 0.4339600 0.4522702
## 4 0.9946989 0.9962711 0.9978433 0.3051808 0.3168875
## 5 0.6337887 0.6441162 0.6544437 1.1586513 1.1857804
## 6 0.8281062 0.8450527 0.8619992 0.7906492 0.8150157
## 7 0.9254248 0.9351679 0.9449110 0.5346896 0.5556477
## 8 0.9628191 0.9693971 0.9759751 0.3810710 0.4015357
## 9 0.4806388 0.4916643 0.5026899 0.9943041 1.0134043
## 10 0.6924412 0.7023515 0.7122619 0.7984517 0.8145411
## 11 0.8218587 0.8383335 0.8548083 0.5967063 0.6178971
## 12 0.9056682 0.9135665 0.9214648 0.4483050 0.4625281
## 13 0.3912238 0.4038981 0.4165723 0.8532835 0.8676875
## 14 0.5702885 0.5798331 0.5893777 0.7454953 0.7588603
## 15 0.7211832 0.7339547 0.7467262 0.6088358 0.6215195
## 16 0.8211173 0.8370655 0.8530137 0.4810203 0.4938210
## UpSystemTime LoQueueLength mQueueLength UpQueueLength
## 1 1.2807317 1.5104841 1.6025268 1.6945695
## 2 0.7417828 0.7180287 0.7901123 0.8621959
## 3 0.4705803 0.3385667 0.3930539 0.4475411
## 4 0.3285943 0.1959176 0.2218722 0.2478267
## 5 1.2129095 2.5494491 2.5918564 2.6342638
## 6 0.8393822 1.5154180 1.5997734 1.6841287
## 7 0.5766058 0.8946711 0.9499252 1.0051792
## 8 0.4220004 0.5675126 0.6226144 0.6777161
## 9 1.0325046 3.0673660 3.0935281 3.1196902
## 10 0.8306305 2.2702008 2.3323181 2.3944353
## 11 0.6390880 1.5417208 1.6329986 1.7242764
## 12 0.4767512 1.0607825 1.1153454 1.1699082
## 13 0.8820915 3.3032825 3.3395759 3.3758693
## 14 0.7722253 2.7709039 2.8157925 2.8606811
## 15 0.6342032 2.1235213 2.1829451 2.2423688
## 16 0.5066217 1.5502433 1.6319491 1.7136550
#=============================================================
# Fitting confidence level of the interval
#=============================================================
for(j in 1:4){
for(n in 1:(N1*N2)){
lo<-t.test(MultiResult[n,j,],conf.level = 0.95)$conf.int[1]
cl<-0.95
h<-1e-3
while(lo<0){
cl<-cl-h
lo<-t.test(MultiResult[n,j,],conf.level = cl)$conf.int[1]
}
if(cl<0.95)
{print("============================================");
print("Intervalo de confianza que debe ser revisado");
print("============================================");
print(paste("n =",n,"; ","j =",j));print(paste("Conf.Level =",cl));
print(t.test(MultiResult[n,j,],conf.level =cl) );
#LoDropped[n] <- t.test(MultiResult[n,j,],conf.level = cl)$conf.int[1];
#UpDropped[n] <- t.test(MultiResult[n,j,],conf.level = cl)$conf.int[2]
}
}
}
#==============================================================
n<-4
j<-1
tsdf<-NoReplications-1
v<-var(MultiResult[n,j,])
s<-sqrt(v)
library(graphics)
plot(function(x) dt((x-mDropped[n])/s, df = tsdf), -8, 10, ylim = c(0, 0.5),
main = "Non-central t - Density", yaxs = "i",ylab=TeX('$F_{T}(t)$'),xlab=TeX('$t$'))
dts<-function(x) dt((x-mDropped[n])/s, df = tsdf)
area <- abs(integrate(dts, LoDropped[n], 0)$value)
# New confident interval
alpha<-0.95-area
print(paste("P{ 0 < mu","<",signif( UpDropped[n],digits = 5) ,"} =",signif(0.95-area,digits = 3)))
## [1] "P{ 0 < mu < 4.2618 } = 0.624"