Implementing a queue class
Who said that queues are boring ? Sure implementing a queue class in R is rather interesting.
My intention is to define a queue class which eases the analysis for scenarios where empirical data may also be used for inter arrival and service time, hence not limited to well known specific probability distributions.
As a result, all queue types classification as M/M, G/M, etc. is handled by the generation process of the arrival time and service time, which shall be input parameters for the queue performance computation.
Based on that, to create an object of class queue shall be enough to specify the length of the internal queue buffer, which is used to store incoming customers waiting to be serviced. FIFO is the default scheduling policy of the queue class I am going to implement. So far, I am not introducing further ones just for simplicity.
As output data I like to have the following:
- inter-arrival time calendar as input
- service time calendar as input
- a data frame where the events sequence is outlined
- a data frame where per each customer the arrival and departure timeline is reported together with the experienced sojourn time and queue wait
- plots of the density of two above
- plot of the events timeline (arrival/departure mix)
- plot of the the queue length time evolution
- plot of the sojourn time evolution
- plot of the queue wait time evolution
- summary of various statistics information about the queue performance
I will take advantage of the S4 class system.
Specifically the queue class shall have aggregated the instance of another class, the queueperf class, which define and implements the queue performance. The queue performance stores the events sequence dataframe and the timeline one which are computed by the performance queue method.
Overall I am going to define and implement:
class queue:
- constructor input parameter: queue buffer length
- methods: performance, plot, summary
- slots: length, arrivaltime, servicetime, performance
class queueperf:
- constructor input parameters: events and timeline dataframes as a list
- methods: plot, summary
- slots: performance
All that leads to the following code:
compute.idletime <- function(event.df) {
len <- length(event.df$queue)
zero.q <- which(event.df$queue.length[-len] == 0)
delta.t <- sapply(zero.q, function(x) { event.df$time[x+1]-event.df$time[x] })
sum(delta.t)
}
compute.totalservicetime <- function(event.df) {
len <- length(event.df$queue)
event.df$time[len]
}
queueperf <- setClass(
"queueperf",
slots = c(
performance = "list"
),
prototype=list(
performance = list(events=data.frame(), timeline=data.frame())
),
validity=function(object)
{
return(TRUE)
}
)
setMethod(f="print",
signature="queueperf",
definition=function(x, ...)
{
}
)
## Creating a generic function for 'print' from package 'base' in the global environment
## [1] "print"
setMethod(f="summary",
signature="queueperf",
definition=function(object, ...)
{
idle <- compute.idletime(object@performance$events)
out1 <- paste("\nidle time:", sep=" ", idle)
servicetime.all <- compute.totalservicetime(object@performance$events)
out2 <- paste("busy time:", sep=" ", servicetime.all-idle)
out3 <- paste("total service time:", sep= " ", servicetime.all)
out <- paste(out1, out2, out3, sep="\n")
cat(out)
cat("\nqueue wait summary:\n")
print(summary(object@performance$timeline$queue.wait))
cat("\nqueue length summary:\n")
print(summary(object@performance$events$queue.length))
cat("\nsojourn time summary:\n")
print(summary(object@performance$timeline$sojourn.time))
}
)
## Creating a generic function for 'summary' from package 'base' in the global environment
## [1] "summary"
setMethod(f="plot",
signature="queueperf",
definition=function(x, y, ...)
{
events <- x@performance$events
old.par <- par("mfrow")
par(mfrow=c(2,2))
lbl <- c("A" = "arrival", "D"= "departure")
f <- factor(events$type)
lev <- levels(f)
plot(events$time, type='n', col=f,
xlab="event", ylab="time",
main ="Events timeline")
points(events$time, cex=0.5, col=f)
legend("topleft", fill = f, cex = 0.7,
legend=lbl[lev], bty="n")
plot(events[,"time"], events[,"queue.length"],
type='h', xlab="time", ylab="queue length",
main="Queue length time evolution")
timeline = x@performance$timeline
plot(timeline[,"sojourn.time"],
type='h', xlab="arrival event", ylab="sojourn time",
main="Sojourn time evolution")
plot(timeline[,"queue.wait"],
type='h', xlab="arrival event", ylab="queue wait",
main="Queue wait time evolution")
par(mfrow=old.par)
}
)
## Creating a generic function for 'plot' from package 'graphics' in the global environment
## [1] "plot"
queue <- setClass(
"queue",
slots = c(
length = "numeric",
performance = "queueperf",
arrival.time = "numeric",
service.time = "numeric"
),
prototype=list(
length = 1
),
validity=function(object)
{
if (object@length > 0) { return(TRUE) }
return(FALSE)
}
)
setGeneric(name="performance",
def=function(object, arrival.time, service.time = NULL)
{
standardGeneric("performance")
}
)
## [1] "performance"
setMethod(f="performance",
signature="queue",
definition=function(object, arrival.time, service.time)
{
arrival.cum <- cumsum(arrival.time)
length <- length(arrival.cum)
if (is.null(service.time) == FALSE && length != length(service.time)) {
stop("Arrival time and service time vectors length mismatch")
}
departure.cum <- vector(mode="numeric", length=length)
queue.length <- vector(mode="numeric", length=length)
object@service.time <- service.time
object@arrival.time <- arrival.time
# first arrival, service and departure
departure.cum[1] <- arrival.time[1] + service.time[1]
queue.length[1] <- 1
# computation of the cumulative departure time
# if the arrival time preceeds the current in service customer departure
# the i-th customer needs to wait, otherwise it can be serviced right away
for (i in 2:length) {
if (arrival.cum[i] < departure.cum[i-1]) {
departure.cum[i] <- departure.cum[i-1] + service.time[i]
} else {
departure.cum[i] <- arrival.cum[i] + service.time[i]
}
}
arrival <- cbind(time=arrival.cum, type=1)
departure <- cbind(time=departure.cum, type=-1)
status <- rbind(arrival, departure)
order <- order(status[,"time"])
events <- status[order,]
queue.length <- cumsum(events[,2])
df <- data.frame(events, queue.length)
df$type[df$type==1] <- "A"
df$type[df$type==-1] <- "D"
df2 <- data.frame(cbind(arrival.cum, departure.cum,
sojourn.time = departure.cum - arrival.cum))
wait <- as.numeric(format(df2[,"sojourn.time"] - service.time, scientific=FALSE))
wait2 <- round(wait, 4)
df2 <- data.frame(cbind(df2, queue.wait=wait2))
object@performance <- queueperf(performance=list(events=df, timeline=df2))
return (object)
}
)
## [1] "performance"
setMethod(f="summary",
signature="queue",
definition=function(object, ...)
{
lambda <- 1/mean(object@arrival.time)
mu <- 1/mean(object@service.time)
out1 <- paste("estimated lambda:", sep=" ", lambda)
out2 <- paste("estimated mu:", sep= " ", mu)
out3 <- paste("estimated rho:", sep= " ", lambda/mu)
out <- paste(out1, out2, out3, sep="\n")
cat(out)
summary(object@performance)
}
)
## [1] "summary"
setMethod(f="plot",
signature="queue",
definition=function(x, y, ...)
{
old.mfrow = par("mfrow")
par(mfrow=c(1,2))
plot(density(x@arrival.time), main = "Inter-arrival time density", ylab="density")
plot(density(x@service.time), main = "Service time density", ylab="density")
readline("Hit <Return> to see next plot: ")
plot(x@performance)
par(mfrow=old.mfrow)
}
)
## [1] "plot"
Now, I am going to:
- allocate a queue object
- compute the inter arrival time
- compute the service time
- compute the queue performance
- show plot and summary of the results
All that for a M/M/1 queue class instance.
a <- queue(length=100)
suppressPackageStartupMessages(library(distr))
set.seed(1023)
arrival.ia <- r(Exp(rate=2))(100)
arrival.ia <- round(arrival.ia, 4)
str(arrival.ia)
## num [1:100] 0.781 0.397 0.982 1.164 0.315 ...
service.time <- r(Exp(rate=2.5))(100)
service.time <- round(service.time, 4)
str(service.time)
## num [1:100] 0.307 0.806 0.056 0.459 0.477 ...
qp <- performance(a, arrival.ia, service.time)
str(qp)
## Formal class 'queue' [package ".GlobalEnv"] with 4 slots
## ..@ length : num 100
## ..@ performance :Formal class 'queueperf' [package ".GlobalEnv"] with 1 slot
## .. .. ..@ performance:List of 2
## .. .. .. ..$ events :'data.frame': 200 obs. of 3 variables:
## .. .. .. .. ..$ time : num [1:200] 0.781 1.088 1.178 1.984 2.16 ...
## .. .. .. .. ..$ type : chr [1:200] "A" "D" "A" "D" ...
## .. .. .. .. ..$ queue.length: num [1:200] 1 0 1 0 1 0 1 2 1 2 ...
## .. .. .. ..$ timeline:'data.frame': 100 obs. of 4 variables:
## .. .. .. .. ..$ arrival.cum : num [1:100] 0.781 1.178 2.16 3.324 3.64 ...
## .. .. .. .. ..$ departure.cum: num [1:100] 1.09 1.98 2.22 3.78 4.26 ...
## .. .. .. .. ..$ sojourn.time : num [1:100] 0.307 0.806 0.056 0.459 0.62 ...
## .. .. .. .. ..$ queue.wait : num [1:100] 0 0 0 0 0.143 ...
## ..@ arrival.time: num [1:100] 0.781 0.397 0.982 1.164 0.315 ...
## ..@ service.time: num [1:100] 0.307 0.806 0.056 0.459 0.477 ...
qp@performance
## An object of class "queueperf"
## Slot "performance":
## $events
## time type queue.length
## 1 0.7808 A 1
## 2 1.0883 D 0
## 3 1.1777 A 1
## 4 1.9837 D 0
## 5 2.1599 A 1
## 6 2.2159 D 0
## 7 3.3243 A 1
## 8 3.6397 A 2
## 9 3.7831 D 1
## 10 3.8857 A 2
## 11 4.2597 D 1
## 12 4.6172 D 0
## 13 4.9729 A 1
## 14 5.0802 D 0
## 15 5.7123 A 1
## 16 5.7368 D 0
## 17 6.1063 A 1
## 18 7.3714 A 2
## 19 7.3889 A 3
## 20 7.4819 A 4
## 21 7.6019 D 3
## 22 7.8885 A 4
## 23 8.2363 D 3
## 24 8.7880 A 4
## 25 8.8562 D 3
## 26 8.8928 D 2
## 27 9.2165 A 3
## 28 9.4101 A 4
## 29 9.7184 D 3
## 30 9.7249 A 4
## 31 9.8230 D 3
## 32 9.9547 A 4
## 33 9.9996 A 5
## 34 10.1840 A 6
## 35 10.4480 A 7
## 36 10.5103 D 6
## 37 10.5442 A 7
## 38 10.7654 D 6
## 39 11.2256 A 7
## 40 11.3103 D 6
## 41 11.5488 D 5
## 42 11.5666 A 6
## 43 11.6372 D 5
## 44 11.6661 D 4
## 45 11.7388 D 3
## 46 11.8696 D 2
## 47 11.9438 D 1
## 48 12.1686 D 0
## 49 12.5632 A 1
## 50 12.7728 A 2
## 51 13.3042 A 3
## 52 13.3106 A 4
## 53 13.4005 A 5
## 54 13.6664 D 4
## 55 13.8718 A 5
## 56 13.9162 D 4
## 57 13.9688 A 5
## 58 13.9727 D 4
## 59 13.9994 D 3
## 60 14.1985 D 2
## 61 14.2226 D 1
## 62 14.8827 D 0
## 63 15.3715 A 1
## 64 15.7387 A 2
## 65 16.0412 A 3
## 66 16.3486 D 2
## 67 16.8209 D 1
## 68 16.8370 A 2
## 69 17.0100 A 3
## 70 17.0401 A 4
## 71 17.3258 A 5
## 72 17.5671 A 6
## 73 17.7596 D 5
## 74 17.7649 D 4
## 75 17.8770 D 3
## 76 17.9215 D 2
## 77 18.0418 D 1
## 78 18.7658 A 2
## 79 19.2570 D 1
## 80 19.4772 A 2
## 81 19.5083 D 1
## 82 19.5622 D 0
## 83 19.7107 A 1
## 84 19.8526 D 0
## 85 20.7625 A 1
## 86 20.9414 D 0
## 87 21.4008 A 1
## 88 21.7578 A 2
## 89 22.0574 D 1
## 90 22.3271 A 2
## 91 22.6150 D 1
## 92 23.0105 D 0
## 93 23.0258 A 1
## 94 23.1921 D 0
## 95 23.4783 A 1
## 96 23.9126 D 0
## 97 24.1362 A 1
## 98 24.3772 A 2
## 99 24.6027 D 1
## 100 25.0094 A 2
## 101 25.4663 D 1
## 102 25.4971 A 2
## 103 26.0343 D 1
## 104 26.1485 D 0
## 105 26.6532 A 1
## 106 26.6989 D 0
## 107 26.7140 A 1
## 108 26.7567 A 2
## 109 27.3624 D 1
## 110 27.4381 D 0
## 111 27.6567 A 1
## 112 28.1931 A 2
## 113 28.2861 D 1
## 114 28.3317 A 2
## 115 28.4973 D 1
## 116 28.7545 D 0
## 117 29.2003 A 1
## 118 29.3864 A 2
## 119 29.7037 A 3
## 120 29.7298 D 2
## 121 29.9772 D 1
## 122 29.9979 A 2
## 123 30.3465 A 3
## 124 30.4511 A 4
## 125 30.6643 D 3
## 126 31.0524 D 2
## 127 31.3780 D 1
## 128 31.4398 D 0
## 129 31.5667 A 1
## 130 31.8654 D 0
## 131 32.0470 A 1
## 132 32.3607 A 2
## 133 32.7069 A 3
## 134 32.9096 A 4
## 135 33.4253 D 3
## 136 33.6367 A 4
## 137 34.0520 D 3
## 138 35.1546 A 4
## 139 35.8088 D 3
## 140 35.9437 A 4
## 141 36.1147 D 3
## 142 36.7328 A 4
## 143 37.3675 D 3
## 144 37.7979 A 4
## 145 37.9639 D 3
## 146 38.0984 A 4
## 147 38.1659 D 3
## 148 38.4045 D 2
## 149 38.5867 D 1
## 150 39.0014 D 0
## 151 39.3553 A 1
## 152 39.3563 A 2
## 153 39.9060 A 3
## 154 40.2063 D 2
## 155 40.4903 D 1
## 156 40.7980 A 2
## 157 41.6902 D 1
## 158 41.7235 A 2
## 159 42.2867 A 3
## 160 42.4280 A 4
## 161 42.4394 A 5
## 162 42.5217 A 6
## 163 42.7379 A 7
## 164 43.0888 D 6
## 165 43.3991 A 7
## 166 43.4902 D 6
## 167 43.4927 D 5
## 168 43.5402 D 4
## 169 43.6019 A 5
## 170 44.3074 D 4
## 171 44.3345 D 3
## 172 44.3561 A 4
## 173 44.6330 A 5
## 174 44.7682 A 6
## 175 45.0346 D 5
## 176 45.0658 A 6
## 177 45.0681 D 5
## 178 45.1558 D 4
## 179 45.2077 A 5
## 180 45.4672 D 4
## 181 45.9107 A 5
## 182 46.0041 A 6
## 183 46.4021 D 5
## 184 46.7747 A 6
## 185 46.9354 D 5
## 186 47.3696 D 4
## 187 47.6280 D 3
## 188 47.7294 A 4
## 189 47.9620 D 3
## 190 48.3034 D 2
## 191 48.3473 D 1
## 192 48.4076 A 2
## 193 48.9265 D 1
## 194 49.0667 A 2
## 195 49.2309 D 1
## 196 49.2851 D 0
## 197 49.6158 A 1
## 198 49.6603 A 2
## 199 49.7648 D 1
## 200 49.9844 D 0
##
## $timeline
## arrival.cum departure.cum sojourn.time queue.wait
## 1 0.7808 1.0883 0.3075 0.0000
## 2 1.1777 1.9837 0.8060 0.0000
## 3 2.1599 2.2159 0.0560 0.0000
## 4 3.3243 3.7831 0.4588 0.0000
## 5 3.6397 4.2597 0.6200 0.1434
## 6 3.8857 4.6172 0.7315 0.3740
## 7 4.9729 5.0802 0.1073 0.0000
## 8 5.7123 5.7368 0.0245 0.0000
## 9 6.1063 7.6019 1.4956 0.0000
## 10 7.3714 8.2363 0.8649 0.2305
## 11 7.3889 8.8562 1.4673 0.8474
## 12 7.4819 8.8928 1.4109 1.3743
## 13 7.8885 9.7184 1.8299 1.0043
## 14 8.7880 9.8230 1.0350 0.9304
## 15 9.2165 10.5103 1.2938 0.6065
## 16 9.4101 10.7654 1.3553 1.1002
## 17 9.7249 11.3103 1.5854 1.0405
## 18 9.9547 11.5488 1.5941 1.3556
## 19 9.9996 11.6372 1.6376 1.5492
## 20 10.1840 11.6661 1.4821 1.4532
## 21 10.4480 11.7388 1.2908 1.2181
## 22 10.5442 11.8696 1.3254 1.1946
## 23 11.2256 11.9438 0.7182 0.6440
## 24 11.5666 12.1686 0.6020 0.3772
## 25 12.5632 13.6664 1.1032 0.0000
## 26 12.7728 13.9162 1.1434 0.8936
## 27 13.3042 13.9727 0.6685 0.6120
## 28 13.3106 13.9994 0.6888 0.6621
## 29 13.4005 14.1985 0.7980 0.5989
## 30 13.8718 14.2226 0.3508 0.3267
## 31 13.9688 14.8827 0.9139 0.2538
## 32 15.3715 16.3486 0.9771 0.0000
## 33 15.7387 16.8209 1.0822 0.6099
## 34 16.0412 17.7596 1.7184 0.7797
## 35 16.8370 17.7649 0.9279 0.9226
## 36 17.0100 17.8770 0.8670 0.7549
## 37 17.0401 17.9215 0.8814 0.8369
## 38 17.3258 18.0418 0.7160 0.5957
## 39 17.5671 19.2570 1.6899 0.4747
## 40 18.7658 19.5083 0.7425 0.4912
## 41 19.4772 19.5622 0.0850 0.0311
## 42 19.7107 19.8526 0.1419 0.0000
## 43 20.7625 20.9414 0.1789 0.0000
## 44 21.4008 22.0574 0.6566 0.0000
## 45 21.7578 22.6150 0.8572 0.2996
## 46 22.3271 23.0105 0.6834 0.2879
## 47 23.0258 23.1921 0.1663 0.0000
## 48 23.4783 23.9126 0.4343 0.0000
## 49 24.1362 24.6027 0.4665 0.0000
## 50 24.3772 25.4663 1.0891 0.2255
## 51 25.0094 26.0343 1.0249 0.4569
## 52 25.4971 26.1485 0.6514 0.5372
## 53 26.6532 26.6989 0.0457 0.0000
## 54 26.7140 27.3624 0.6484 0.0000
## 55 26.7567 27.4381 0.6814 0.6057
## 56 27.6567 28.2861 0.6294 0.0000
## 57 28.1931 28.4973 0.3042 0.0930
## 58 28.3317 28.7545 0.4228 0.1656
## 59 29.2003 29.7298 0.5295 0.0000
## 60 29.3864 29.9772 0.5908 0.3434
## 61 29.7037 30.6643 0.9606 0.2735
## 62 29.9979 31.0524 1.0545 0.6664
## 63 30.3465 31.3780 1.0315 0.7059
## 64 30.4511 31.4398 0.9887 0.9269
## 65 31.5667 31.8654 0.2987 0.0000
## 66 32.0470 33.4253 1.3783 0.0000
## 67 32.3607 34.0520 1.6913 1.0646
## 68 32.7069 35.8088 3.1019 1.3451
## 69 32.9096 36.1147 3.2051 2.8992
## 70 33.6367 37.3675 3.7308 2.4780
## 71 35.1546 37.9639 2.8093 2.2129
## 72 35.9437 38.1659 2.2222 2.0202
## 73 36.7328 38.4045 1.6717 1.4331
## 74 37.7979 38.5867 0.7888 0.6066
## 75 38.0984 39.0014 0.9030 0.4883
## 76 39.3553 40.2063 0.8510 0.0000
## 77 39.3563 40.4903 1.1340 0.8500
## 78 39.9060 41.6902 1.7842 0.5843
## 79 40.7980 43.0888 2.2908 0.8922
## 80 41.7235 43.4902 1.7667 1.3653
## 81 42.2867 43.4927 1.2060 1.2035
## 82 42.4280 43.5402 1.1122 1.0647
## 83 42.4394 44.3074 1.8680 1.1008
## 84 42.5217 44.3345 1.8128 1.7857
## 85 42.7379 45.0346 2.2967 1.5966
## 86 43.3991 45.0681 1.6690 1.6355
## 87 43.6019 45.1558 1.5539 1.4662
## 88 44.3561 45.4672 1.1111 0.7997
## 89 44.6330 46.4021 1.7691 0.8342
## 90 44.7682 46.9354 2.1672 1.6339
## 91 45.0658 47.3696 2.3038 1.8696
## 92 45.2077 47.6280 2.4203 2.1619
## 93 45.9107 47.9620 2.0513 1.7173
## 94 46.0041 48.3034 2.2993 1.9579
## 95 46.7747 48.3473 1.5726 1.5287
## 96 47.7294 48.9265 1.1971 0.6179
## 97 48.4076 49.2309 0.8233 0.5189
## 98 49.0667 49.2851 0.2184 0.1642
## 99 49.6158 49.7648 0.1490 0.0000
## 100 49.6603 49.9844 0.3241 0.1045
The queue performance provides both events and timeline dataframes. Based on those dataframes, herein I show the resulting plots and summary of our queue instance.
plot(qp)
## Hit <Return> to see next plot:
summary(qp)
## estimated lambda: 2.01368094836318
## estimated mu: 2.41728649921317
## estimated rho: 0.833033630485519
## idle time: 7.83489999999999
## busy time: 42.1495
## total service time: 49.9844
## queue wait summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.1016 0.6066 0.7188 1.1000 2.8990
##
## queue length summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 1.00 2.00 2.62 4.00 7.00
##
## sojourn time summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0245 0.6507 1.0070 1.1320 1.5880 3.7310
The queue plots comprise the probability density functions for arrivals and service time and the queue performance related ones, which are additional four plots:
- the events timeline plot shows the sequence of arrival and departure events during the time
- the queue length time evolution shows per each time tick the corresponding queue buffer length
- the sojourn time evolution outlines per each customer arrival event the sojourn time within the queue system (queue wait + service time)
- the queue wait time shows per each customer the time spent inside the buffer queue
The summary provides with additional details about the queue performance in terms of load factor, idle time, busy time and summaries of queue wait, queue length and sojourn time.
Further, if I need to run the same analysis for a G/G/1 queue, I would have the following.
arrival.ia <- r(Norm(mean=2, sd=0.2))(100)
arrival.ia <- round(arrival.ia, 4)
str(arrival.ia)
## num [1:100] 1.72 2.38 1.84 2.14 1.84 ...
service.time <- r(Norm(mean=1.7, sd=0.8))(100)
service.time <- round(service.time, 4)
str(service.time)
## num [1:100] 0.93 1.813 1.568 0.707 1.028 ...
qp <- performance(a, arrival.ia, service.time)
plot(qp)
## Hit <Return> to see next plot:
summary(qp)
## estimated lambda: 0.504312629450748
## estimated mu: 0.545381464338324
## estimated rho: 0.924697046795673
## idle time: 13.8475000000001
## busy time: 185.0783
## total service time: 198.9258
## queue wait summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.07067 0.91740 0.92240 1.52600 3.18200
##
## queue length summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 1.00 1.00 1.33 2.00 3.00
##
## sojourn time summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6361 1.9090 2.8750 2.7560 3.5070 5.0790
Hence, you may see how much flexible is our queue class definition with respect to the arrival and service time definitions and specifications.
No comments:
Post a Comment
Note: Only a member of this blog may post a comment.