Math 365 Queues

Math 365 Queues The objective of this lab is to further study queueing as a type of birth and death process using a combination of simulations and theoretical analysis. 1 Simulation of M/M/1 Queue In general for birth and death queueing processes, we let λn denote the arrival rate when there are currently n customers, and µn denote the departure rate, where both arrivals and departures are treated as exponential processes. That is, if there are currently n customers, the time until the next event is an exponential random variable with rate λn + µn ; the probability that the next event n n and the probability that the next event is a departure is λnµ+µ . is an arrival is λnλ+µ n n The M/M/1 queue has constant rates λn = λ for n ≥ 0 and µn = µ for n ≥ 1 (µ0 = 0). We’ll start by running some simulations of an M/M/1 queue to see some examples of what such processes look like for the positive recurrent case (λ < µ) in which we expect to repeatedly return to an empty queue. Feel free to vary the values of λ and µ in the simulation. If it seems to be running too long, halt it and try a different pair of values (make sure λ < µ so you don’t get exploding queues). First initialize some variables: # M/M/1 queue simulator lambda <- 1 # arrival rate mu <- 4 # service rate duration <- 10000 # total duration of the simulation t <- 0 # current time in the simulation queue <- 0 # start with empty queue s <- 0 # running sum for computing average queue length We assume the queue is initially empty, so generate a first arrival to the queue with rate λ: # first arrival to start process T1 <- rexp(1,rate=lambda) currentqueue <- 1 eventsTime <- T1 t <- T1 nEvents <- 1 # total number of events that have occurred The next event can be either an arrival or departure. The time T1 until that next event is an exponential random variable with rate λ + µ. To determine whether that event is an arrival or a λ departure, generate a random number p between 0 and 1; if p < λ+µ , then the event is an arrival, otherwise it’s a departure. while (t<duration) { nEvents <- nEvents+1 if(currentqueue>0) { # nonempty queue T1 <- rexp(1,rate=lambda+mu) # time until next event # is event an arrival or departure? p <- runif(1,0,1) queue[nEvents] <- currentqueue # how many have been in queue before this new event currentqueue <- ifelse(p<lambda/(lambda+mu), currentqueue+1, # arrival currentqueue-1) # departure } else { # empty queue T1 <- rexp(1,rate=lambda) queue[nEvents] <- currentqueue currentqueue <- 1 } t <- t+T1 # time of next arrival eventsTime[nEvents] <- T1 # inter-event time s <- s+T1*queue[nEvents] # spent T1 at nth queue length } Note that the R code above also has to check at each step whether the queue is empty, in which case only an arrival is possible. Once the simulation is complete, you can plot the queue length over time: plot (cumsum(eventsTime),queue,type="s ", xlab="Time ", ylab=" Queue length", main=paste("M/M/1 Simulation, mean queue length=",s/t)) 2 Analysis of M/M/1 Queue What do we expect the average queue length to be? What proportion of the time do we expect the queue to be empty over the long run? What is the average time customers will spend in the queue? We can do some straightforward analysis to answer these questions. Begin by calculating the limiting probability distribution π(n) for n ≥ 0. Exercise 1 For the M/M/1 queue, π(n) satisfies the equilibrium conditions (λ + µ)π(n) = λπ(n − 1) + µπ(n + 1) for n ≥ 1 and λπ(0) = µπ(1). First solve for π(1), then for π(2), and then P recognize the pattern that is occurring to find the general solution in terms of π(0). Enforce ∞ n=0 π(n) = 1 to find π(0) (use the geometric series formula to simplify). State the resulting formula for π(n), which gives the long-run proportion of time that the queue length equals n (equivalently, the probability at any given time t that the queue length is n). Exercise 2 Run a simulation and check that the observed proportion of time in which the queue length is 0 agrees well with the theoretical value for π(0) found in Exercise 1. An easy way to compute this value in R is the following: sum(eventsTime[which(queue==0)])/t Exercise 3 Compare the observed proportion of time in which the queue length is 1 with the theoretical value for π(1). Exercise 4 Let L be the expected value of the queue length (average number of customers): L= ∞ X nπ(n). n=0 Find a simple formula in your solution for π(n) from Exercise 1 and simplifying P for Lnby substituting r using the formula ∞ nr = (which can be obtained by taking a derivative of the geometric n=1 (1−r)2 series formula). Exercise 5 Compare the observed average queue length (s/t in R) with the theoretical value for L. Exercise 6 Finally, the average time W each customer spends in line is the average queue length divided by the rate of arrivals: W = L/λ. This relation is known as “Little’s law” State a formula for W using the result of Exercise 4.