Queuing models in R, $\lambda$ Little

Question

It's noted that the number of folks in a stationary system will maintain an average equal to the rate of arrival multiplied by the mean of the service distribution.

The formula $L = \lambda w$ is valid for any queueing model in steady-state, where $L$ and $w$ are long-term steady-state average values respectively, and $\lambda$ denotes the arrival rate to the system.

We can add up the total service time in the system as follows:

$$\sum w_{j} = \sum i T_i$$

where $T_i$ represents the time units in which $i$ entities were in the system.

But things are always more interesting with sample sizing and simulating results, so say in $5000$ iterations we estimate a state of a system in 1-min intervals, between the first minute and the last arrival at a determined time.

So suppose we use a random interarrival rate of $\lambda = 2$ per minute and the service distribution is $N(8,1)$ minutes for the system.

How can I simulate this model in R, using rexp() and rnorm()? I also want to display in ts() to show in time plot.

Answer 1

Not directly answering your question of how to code it manually but for discrete simulation of queues in R I would strongly recommend the simmer package. The minimal code for your example would look like this (adapted from the tutorial).

library(simmer)
library(simmer.plot)
lambda <- 2

queue <- trajectory() %>%
  seize("server", amount=1) %>%
  timeout(function() {rnorm(1, 8, 1)}) %>%
  release("server", amount=1)

env <- simmer() %>%
  add_resource("server", capacity=1) %>%
  add_generator("arrival", queue, function() {rnorm(1, lambda)}) %>%
  run(until=100000)


resources <- get_mon_resources(env)
arrivals <- get_mon_arrivals(env)

plot(resources, metric = "usage", items = "server")
plot(arrivals, metric = "flow_time")

However, looking at the plot of flow times I am unsure if you conceived a stable queuing system - for mean = 1 in the service time it looks better:

Service Time Mean 8:

Service Time Mean 1:

Answer 2

So far this is what I have come up with

lambda <- 2
interarrivals <- rexp(5000,lambda)  ## (2 items per minute)

Provided the $\mu$ we expect that the interarrivals is about half a minute

mean(interarrivals) <- 0.516
service.times <- rnorm(5000, mean=8,sd=1)

where the service distribution is $N(8,1)$

arrival.times <- cumsum(interarrivals)
departure.times <- arrival.times+service.times

Where I am having issues is determining how many individuals remain in the system at time 1-min to the last.

I am having issues conceptualizing the model - this is what you are supposed to estimate (Little's law), which states that rate of arrival = 2 units per minute × mean of the service distribution = 8.