library(simmer)
library(simmer.plot)
library(parallel)
set.seed(1234)

The M/M/1 system

In Kendall’s notation, an M/M/1 system has exponential arrivals (M/M/1), a single server (M/M/1) with exponential service time (M/M/1) and an inifinite queue (implicit M/M/1/\(\infty\)). For instance, people arriving at an ATM at rate \(\lambda\), waiting their turn in the street and withdrawing money at rate \(\mu\).

Let us recall the basic parameters/performance measures of this system:

\[\begin{aligned} \rho &= \frac{\lambda}{\mu} &&\equiv \mbox{Server} \hspace{0.5ex} \mbox{utilization,}\hspace{0.5ex} \mbox{traffic}\hspace{0.5ex} \mbox{intensity} \\ L &= \frac{\rho}{1-\rho} &&\equiv \mbox{Average} \hspace{0.5ex}\mbox{number} \hspace{0.5ex}\mbox{of} \hspace{0.5ex}\mbox{customers} \hspace{0.5ex}\mbox{in} \hspace{0.5ex}\mbox{the} \hspace{0.5ex}\mbox{system} \hspace{0.5ex}\mbox{(queue}\hspace{0.5ex} \mbox{+} \hspace{0.5ex} \mbox{server)} \\ W &= \frac{L}{\lambda} &&\equiv \mbox{Average} \hspace{0.5ex}\mbox{time} \hspace{0.5ex}\mbox{in} \hspace{0.5ex}\mbox{the}\hspace{0.5ex} \mbox{system}\hspace{0.5ex} \mbox{(queue}\hspace{0.5ex} \mbox{+ server)} \hspace{0.5ex}\mbox{[Little's} \hspace{0.5ex}\mbox{law]} \\ \end{aligned}\]

whenever \(\rho < 1\), so the system is in steady state. If \(\rho < 1\) is not true, it means that the system is unstable: there are more arrivals than the server is capable of handling, and the queue will grow indefinitely.

The simulation of an M/M/1 system is quite simple using simmer:

lambda <- 2
mu <- 4
rho <- lambda/mu # = 2/4

mm1.trajectory <- trajectory() %>%
  seize("resource", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("resource", amount=1)

mm1.env <- simmer() %>%
  add_resource("resource", capacity=1, queue_size=Inf) %>%
  add_generator("arrival", mm1.trajectory, function() rexp(1, lambda)) %>% 
  run(until=10000)

COMMANDS in package “simmer.plot”

plot(x, what = c(“resources”, “arrivals”, “attributes”), metric = NULL, …)

Arguments

x - a single simmer environment or a list of environments representing several replications.

what - type of plot, one of c(“resources”, “arrivals”, “attributes”).

metric - specific metric for each type of plot.

       what = "resources" one of c("usage", "utilization").

       what = "arrivals" one of c("activity_time", "waiting_time", "flow_time").

… further arguments for each kind of plot.

       what   = "resources" all metrics names the name of the resource(s) 
       (a single string or a character vector) to show.
       
       metric = "usage" items the components of the resource to be plotted,
       one or more of c("system", "queue", "server").
       
       steps if TRUE, shows the instantaneous usage instead of the cumulative average.

The simmer.plot package provides convenience plotting functions to quickly visualise the usage of a resource over time, for instance. Below, we can see how the simulation converges to the theoretical average number of customers in the system.

# Theoretical value
mm1.L <- rho/(1-rho)

# Evolution of the average number of customers in the system
plot(mm1.env, "resources", "usage", "resource", items="system") +
  geom_hline(yintercept=mm1.L)

Also, we can see how the simulation converges to the theoretical average number of customers in the queue.

# Theoretical value
mm1.Lq <- (rho**2)/(1-rho)


# Evolution of the average number of customers in the queue
plot(mm1.env, "resources", "usage", "resource", items="queue") +
  geom_hline(yintercept=mm1.Lq)

Also, it is possible to visualise, for instance, the instantaneous usage of individual elements by playing with the parameters items and steps.

plot(mm1.env, "resources", "usage", "resource", items=c("queue", "server"), steps=TRUE) +
  xlim(0, 20) + ylim(0, 4)
#> Warning: Removed 80954 rows containing missing values (geom_path).

#> Warning: Removed 80954 rows containing missing values (geom_path).

Experimentally, we obtain the time spent by each customer in the system and we compare the average with the theoretical expression.

mm1.arrivals <- get_mon_arrivals(mm1.env)
mm1.t_system <- mm1.arrivals$end_time - mm1.arrivals$start_time
lambda_eff<- lambda
mm1.W <- mm1.L / lambda_eff # Little's Law
mm1.W ; mean(mm1.t_system)
#> [1] 0.5
#> [1] 0.5072202

It seems that it matches the theoretical value pretty well. Let’s take a closer look, by looking at the 95% confidence intervals.

Note (see Wikipedia) In statistics, a confidence interval (CI) is a type of interval estimate (of a population parameter) that is computed from the observed data. The confidence level is the frequency (i.e., the proportion) of confidence intervals that contain the true value of their corresponding parameter. In other words, if confidence intervals are constructed using a given confidence level in an infinite number of independent experiments, the proportion of those intervals that contain the true value of the parameter will match the confidence level. Confidence intervals consist of a range of values (interval) that act as good estimates of the unknown population parameter. However, the interval computed from a particular sample does not necessarily include the true value of the parameter. Since the observed data are random samples from the true population, the confidence interval obtained from the data is also random.


The replications needed for the construction of the CI can be done with standard R tools:

envs <- mclapply(1:50, function(i) {
  simmer() %>%
    add_resource("resource", capacity=1, queue_size=Inf) %>%
    add_generator("arrival", mm1.trajectory, function() rexp(1, lambda)) %>%
    run(2000) %>%
    wrap()
}, mc.set.seed=TRUE)

t_system <- get_mon_arrivals(envs) %>%
  dplyr::mutate(t_system = end_time - start_time) %>%
  dplyr::group_by(replication) %>%
  dplyr::summarise(mean = mean(t_system))

n_system <-get_mon_resources(envs)%>%
  dplyr::group_by(replication) %>%
  dplyr::summarise(mean = sum(head(system, -1) * diff(time)) / tail(time, 1))

t.test(t_system$mean)
#> 
#>  One Sample t-test
#> 
#> data:  t_system$mean
#> t = 150.79, df = 49, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#>  0.4966873 0.5101047
#> sample estimates:
#> mean of x 
#>  0.503396
t.test(n_system$mean)
#> 
#>  One Sample t-test
#> 
#> data:  n_system$mean
#> t = 135.01, df = 49, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#>  0.997687 1.027836
#> sample estimates:
#> mean of x 
#>  1.012761

Finally, the inverse of the mean difference between arrivals is the effective rate, which matches (approx.) the real lambda because there are no rejections.

lambda; 1/mean(diff(subset(mm1.arrivals, finished==TRUE)$start_time))
#> [1] 2
#> [1] 2.027843

The M/M/c system

In Kendall’s notation, an M/M/c system has exponential arrivals (M/M/c), c servers (M/M/c) with exponential service time (M/M/c) and an inifinite queue (implicit M/M/c/\(\infty\)). For instance, people arriving at rate \(\lambda\) at a bank with c tellers, waiting their turn and being served at rate \(\mu\).

Below are the basic parameters/performance measures of this system:

\[\begin{aligned} \mbox{Denote by} \hspace{2ex} \rho &= \frac{\lambda}{\mu} &&\equiv \mbox{the traffic intensity} \\ \frac{\rho}{c} &= \frac{\lambda}{c \mu} &&\equiv \mbox{Server utilization}\\ \pi_0 &= \left(\sum_{j=0}^{c-1} \frac{\rho^j}{j!} + \frac{\rho^c}{c!}\frac{1}{(1 - \frac{\rho}{c})}\right)^{-1} &&\equiv \mbox{probability system (queue + servers) is empty} \\ L_q &= \frac{\rho^c}{c!}\pi_0 \frac{\frac{\rho}{c}}{(1 - \frac{\rho}{c})^2} &&\equiv \mbox{Average number of customers in the queue} \\ L_s &= \rho &&\equiv \mbox{Average number of customers in service} \\ W_q &= \frac{L_q}{\lambda} &&\equiv \mbox{Average time in the system (queue + servers) [Little's law]} \\ \end{aligned}\]

whenever \(\frac{\rho}{c} < 1\), so the system is in steady state. If \(\frac{\rho}{c} < 1\) is not true, it means that the system is unstable: there are more arrivals than the servers are capable of handling, and the queue will grow indefinitely.

The simulation of an M/M/c system is simple using simmer:

lambda <- 5
mu <- 4
c<-2
rho <- lambda/mu
rho_c <- lambda/(c*mu) # = 5/8

mmc.trajectory <- trajectory() %>%
  seize("resource", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("resource", amount=1)

mmc.env <- simmer() %>%
  add_resource("resource", capacity=2, queue_size=Inf) %>%
  add_generator("arrival", mmc.trajectory, function() rexp(1, lambda)) %>% 
  run(until=2000)

Our simmer.plot package provides convenience plotting functions to quickly visualise the usage of a resource over time, for instance. Down below, we can see how the simulation converges to the theoretical average number of customers in the system.

# Theoretical values
ind <- 0:(c-1)

pi_0 <-1/(sum(rho**ind/factorial(ind)) + (rho**c/factorial(c))*(1/(1-rho/c)))
mmc.L <- rho +  ((rho**c/factorial(c))*pi_0)*((rho/c)/((1-rho/c)**2))
mmc.Lq <- mmc.L - rho

# Evolution of the average number of customers in the system
plot(mmc.env, "resources", "usage", "resource", items="system") +
  geom_hline(yintercept=mmc.L) + ylim(0, 3)

# Evolution of the average number of customers in the queue
plot(mmc.env, "resources", "usage", "resource", items="queue") + 
  geom_hline(yintercept=mmc.Lq) + ylim(0, 1.5)

It is possible also to visualise, for instance, the instantaneous usage of individual elements by playing with the parameters items and steps.

plot(mmc.env, "resources", "usage", "resource", items=c("queue", "server"), steps=TRUE) +
  xlim(0, 20) + ylim(0, 5)
#> Warning: Removed 39794 rows containing missing values (geom_path).

#> Warning: Removed 39794 rows containing missing values (geom_path).

Experimentally, we obtain the time spent by each customer in the system and we compare the average with the theoretical expression.

mmc.arrivals <- get_mon_arrivals(mmc.env)
mmc.t_system <- mmc.arrivals$end_time - mmc.arrivals$start_time
lambda_eff<-lambda
mmc.W <- mmc.L / lambda_eff
mmc.W ; mean(mmc.t_system)
#> [1] 0.4102564
#> [1] 0.407997

It seems that it matches the theoretical value pretty well. Again, looking at the confidence intervals. Replication can be done with standard R tools:

envs <- mclapply(1:50, function(i) {
  simmer() %>%
    add_resource("resource", capacity=c, queue_size=Inf) %>%
    add_generator("arrival", mmc.trajectory, function() rexp(1, lambda)) %>%
    run(500) %>%
    wrap()
}, mc.set.seed=TRUE)

t_system <- get_mon_arrivals(envs) %>%
  dplyr::mutate(t_system = end_time - start_time) %>%
  dplyr::group_by(replication) %>%
  dplyr::summarise(mean = mean(t_system))

n_system <-get_mon_resources(envs)%>%
  dplyr::group_by(replication) %>%
  dplyr::summarise(mean = sum(head(system, -1) * diff(time)) / tail(time, 1))

t.test(t_system$mean)
#> 
#>  One Sample t-test
#> 
#> data:  t_system$mean
#> t = 102.76, df = 49, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#>  0.4016513 0.4176745
#> sample estimates:
#> mean of x 
#> 0.4096629
t.test(n_system$mean)
#> 
#>  One Sample t-test
#> 
#> data:  n_system$mean
#> t = 88.139, df = 49, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#>  2.004372 2.097904
#> sample estimates:
#> mean of x 
#>  2.051138

Finally, the inverse of the mean difference between arrivals is the effective rate, which matches (approx.) the real lambda because there are no rejections.

lambda; 1/mean(diff(subset(mmc.arrivals, finished==TRUE)$start_time))
#> [1] 5
#> [1] 5.03123

M/M/c/k systems

An M/M/c/k system keeps exponential arrivals and service times, but has more than one server in general and a finite queue, which often is more realistic. For instance, a router may have several processor to handle packets, and the in/out queues are necessarily finite. We consider the case when \(\lambda \ne c\mu\).

Below are the basic parameters/performance measures of this system:

\[\begin{aligned} &\rho = \frac{\lambda}{\mu} &&\equiv \mbox{the traffic intesity} \\ %\frac{\rho}{c} &= \frac{\lambda}{c \mu} &&\equiv \mbox{Server utilization}\\ \pi_0 &= \left(\sum_{j=0}^{c-1} \frac{\rho^j}{j!} + \frac{\rho^c}{c!}\frac{1-(\frac{\rho}{c})^{K-c+1}}{(1 - \frac{\rho}{c})}\right)^{-1} &&\equiv \mbox{probability system (queue + servers) is empty} \\ L_q &= \pi_0 \frac{\rho^{c+1}}{c! c}\frac{1}{(1 - \frac{\rho}{c})^2}\left[1-\left(\frac{\rho}{c}\right)^{K-c+1} - (K-c+1) \left(\frac{\rho}{c}\right)^{K-c} \left(1 - \frac{\rho}{c}\right) \right] &&\equiv \mbox{Average number of customers in the queue} \\ L &= L_q +c - \pi_0 \sum_{n=0}^{c-1} \frac{(c-n) \rho^n}{n!} &&\equiv \mbox{Average number of customers in the system} \\ W &= \frac{L}{\lambda_{eff}} &&\equiv \mbox{Average time in the system [Little's law]} \\ \end{aligned}\]

This is the simulation of an M/M/2/3 system (2 server, 1 position in queue).

lam <-4
mu <- 1 

K<-3
c<-2

S<-2000

mm23.trajectory <- trajectory() %>%
  seize("resource", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("resource", amount=1)

mm23.env <- simmer() %>%
  add_resource("resource", capacity = 2, queue_size = 1) %>%
  add_generator("arrival", mm23.trajectory, function() rexp(1, lam)) %>%
  run(until=S)


# Theoretical values

rho<-lam/mu

ind <- 0:(c-1)

pi_0 <-1/(sum(rho**ind/factorial(ind)) + (rho**c/factorial(c))*
            (1 - (rho/c)**(K-c+1))/(1-rho/c))

mm23.Lq <- pi_0*((rho**(c+1))/(factorial(c)*c*(1-(rho/c))**2))*
  (1 - (rho/c)**(K-c+1)- (K-c+1)*(rho/c)**(K-c)*(1-(rho/c)))

mm23.L <- mm23.Lq + c - pi_0* sum((c-ind)*(rho**ind)/factorial(ind))

We can see how the simulation converges to the theoretical average number of customers in the system.

# Evolution of the average number of customers in the system
plot(mm23.env, "resources", "usage", "resource", items= "system") +
  geom_hline(yintercept=mm23.L)

We can see how the simulation converges to the theoretical average number of customers in the queue.


# Evolution of the average number of customers in the queue
plot(mm23.env, "resources", "usage", "resource", "queue") + 
  geom_hline(yintercept=mm23.Lq)

It is possible also to visualise, for instance, the instantaneous usage of individual elements by playing with the parameters items and steps.

plot(mm23.env, "resources", "usage", "resource", items=c("queue", "server"), steps=TRUE) +
  xlim(0, 20) + ylim(0, 3)
#> Warning: Removed 14146 rows containing missing values (geom_path).

#> Warning: Removed 14146 rows containing missing values (geom_path).

Experimentally, we obtain the time spent by each customer in the system and we compare the average with the theoretical expression.

mm23.arrivals <- get_mon_arrivals(mm23.env)
mm23.t_system <- subset(mm23.arrivals, finished ==TRUE)

lambda_eff<-1/mean(diff(subset(mm23.arrivals, finished==TRUE)$start_time))

mm23.t_system$delay<- mm23.t_system$end_time -mm23.t_system$start_time

mm23.W <- mm23.L/lambda_eff
mm23.W 
#> [1] 1.313265
mean(mm23.t_system$delay)
#> [1] 1.31924

It seems that it matches the theoretical value pretty well. Again, looking at the confidence intervals.

envs <- mclapply(1:100, function(i) {
  simmer() %>%
    add_resource("resource", capacity=2, queue_size=1) %>%
    add_generator("arrival", mm23.trajectory, function() rexp(1, lam)) %>%
    run(2000) %>%
    wrap()
}, mc.set.seed=TRUE)

t_system <- get_mon_arrivals(envs) %>%
  dplyr::mutate(wait = end_time - start_time) %>%
  dplyr::group_by(replication) %>%
  dplyr::summarise(mean = mean(wait[finished]))

n_system <-get_mon_resources(envs)%>%
  dplyr::group_by(replication) %>%
  dplyr::summarise(mean = sum(head(system, -1) * diff(time)) / tail(time, 1))

t.test(t_system$mean)
#> 
#>  One Sample t-test
#> 
#> data:  t_system$mean
#> t = 567.67, df = 99, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#>  1.302206 1.311341
#> sample estimates:
#> mean of x 
#>  1.306774
t.test(n_system$mean)
#> 
#>  One Sample t-test
#> 
#> data:  n_system$mean
#> t = 1305.8, df = 99, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#>  2.337708 2.344823
#> sample estimates:
#> mean of x 
#>  2.341265

Finally, the inverse of the mean difference between arrivals is the effective rate, which is different from the real lambda because there are rejections in the system.

lam; 1/mean(diff(subset(mm23.arrivals, finished==TRUE)$start_time))
#> [1] 4
#> [1] 1.785494

In this case, there are rejections when the queue is full.

mm23.arrivals <- get_mon_arrivals(mm23.env)
mm23.arrivals %>%
  dplyr::summarise(rejection_rate = sum(!finished)/length(finished))
#>   rejection_rate
#> 1      0.5469782

M/M/c/c systems (BCC (blocked customers cleared) system )

An M/M/c/c system keeps exponential arrivals and service times, has more than one server in general and no queue. We consider the case when \(\lambda \ne c\mu\).

Below are the basic parameters/performance measures of this system: same as for M/M/c/K, but set K = c

\[\begin{aligned} &\rho = \frac{\lambda}{\mu} &&\equiv \mbox{the traffic intesity} \\ %\frac{\rho}{c} &= \frac{\lambda}{c \mu} &&\equiv \mbox{Server utilization}\\ \pi_0 &= \left(\sum_{j=0}^{c-1} \frac{\rho^j}{j!} + \frac{\rho^c}{c!}\frac{1-(\frac{\rho}{c})^{K-c+1}}{(1 - \frac{\rho}{c})}\right)^{-1} &&\equiv \mbox{probability system (queue + servers) is empty} \\ L_q &= \pi_0 \frac{\rho^{c+1}}{c! c}\frac{1}{(1 - \frac{\rho}{c})^2}\left[1-\left(\frac{\rho}{c}\right)^{K-c+1} - (K-c+1) \left(\frac{\rho}{c}\right)^{K-c} \left(1 - \frac{\rho}{c}\right) \right] &&\equiv \mbox{Average number of customers in the queue} \\ L &= L_q +c - \pi_0 \sum_{n=0}^{c-1} \frac{(c-n) \rho^n}{n!} &&\equiv \mbox{Average number of customers in the system} \\ W &= \frac{L}{\lambda_{eff}} &&\equiv \mbox{Average time in the system [Little's law]} \\ \end{aligned}\]

This is the simulation of an M/M/3/3 system (3 servers, no queue).

lam <-4
mu <- 1 

K<-3
c<-3

S<-2000

mm33.trajectory <- trajectory() %>%
  seize("resource", amount=1) %>%
  timeout(function() rexp(1, mu)) %>%
  release("resource", amount=1)

mm33.env <- simmer() %>%
  add_resource("resource", capacity = 3, queue_size = 0) %>%
  add_generator("arrival", mm33.trajectory, function() rexp(1, lam)) %>%
  run(until=S)


# Theoretical values

rho<-lam/mu

ind <- 0:(c-1)

pi_0 <-1/(sum(rho**ind/factorial(ind)) + (rho**c/factorial(c))*
            (1 - (rho/c)**(K-c+1))/(1-rho/c))

mm33.Lq <- pi_0*((rho**(c+1))/(factorial(c)*c*(1-(rho/c))**2))*
  (1 - (rho/c)**(K-c+1)-(K-c+1)*(rho/c)**(K-c)*(1-(rho/c)))

mm33.L <- mm33.Lq + c - pi_0* sum((c-ind)*(rho**ind)/factorial(ind))

We can see how the simulation converges to the theoretical average number of customers in the system.

# Evolution of the average number of customers in the system
plot(mm33.env, "resources", "usage", "resource", items= "system") +  
  geom_hline(yintercept=mm33.L)

The length of the queue in this queueing model is always zero, so no need for pictorial illustration.

It is possible also to visualise, for instance, the instantaneous usage of individual elements by playing with the parameters items and steps.

plot(mm33.env, "resources", "usage", "resource", items=c("queue", "server"), steps=TRUE) +
  xlim(0, 20) + ylim(0, 3)
#> Warning: Removed 17188 rows containing missing values (geom_path).

#> Warning: Removed 17188 rows containing missing values (geom_path).

Experimentally, we obtain the time spent by each customer in the system and we compare the average with the theoretical expression.

mm33.arrivals <- get_mon_arrivals(mm33.env)
mm33.t_system <- subset(mm33.arrivals, finished ==TRUE)

lambda_eff<-1/mean(diff(subset(mm33.arrivals, finished==TRUE)$start_time))

mm33.t_system$delay<- mm33.t_system$end_time -mm33.t_system$start_time

mm33.W <- mm33.L/lambda_eff
mm33.W 
#> [1] 1.008864
mean(mm33.t_system$delay)
#> [1] 1.020125

It seems that it matches the theoretical value pretty well. Again, looking at the confidence intervals:

envs <- mclapply(1:100, function(i) {
  simmer() %>%
    add_resource("resource", capacity=3, queue_size=0) %>%
    add_generator("arrival", mm33.trajectory, function() rexp(1, lam)) %>%
    run(2000) %>%
    wrap()
}, mc.set.seed=TRUE)

t_system <- get_mon_arrivals(envs) %>%
  dplyr::mutate(wait = end_time - start_time) %>%
  dplyr::group_by(replication) %>%
  dplyr::summarise(mean = mean(wait[finished]))

n_system <-get_mon_resources(envs)%>%
  dplyr::group_by(replication) %>%
  dplyr::summarise(mean = sum(head(system, -1) * diff(time)) / tail(time, 1))

t.test(t_system$mean)
#> 
#>  One Sample t-test
#> 
#> data:  t_system$mean
#> t = 608.53, df = 99, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#>  0.9969781 1.0035010
#> sample estimates:
#> mean of x 
#>   1.00024
t.test(n_system$mean)
#> 
#>  One Sample t-test
#> 
#> data:  n_system$mean
#> t = 1277.8, df = 99, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#>  2.194040 2.200865
#> sample estimates:
#> mean of x 
#>  2.197452

Finally, the inverse of the mean difference between arrivals is the effective rate, which is different from the real lambda, because there are rejections in the system.

lam; 1/mean(diff(subset(mm33.arrivals, finished==TRUE)$start_time))
#> [1] 4
#> [1] 2.177877

In this case, there are rejections when all servers are busy.

mm33.arrivals <- get_mon_arrivals(mm33.env)
mm33.arrivals %>%
  dplyr::summarise(rejection_rate = sum(!finished)/length(finished))
#>   rejection_rate
#> 1      0.4546025