Lecture 6
Stefanka Chukova
6. Introduction to queues (queueing systems)
Kendall-Lee notations for queues
In 1951 Kendall and Lee introduced a very convenient notation for the description of the queueing systems. It consists of six labels, put in a specific order \(1/2/3/4/5/6\):
arrivals
service
number of servers
service discipline
maximum number of customers in the system
size of the customer population
Moreover, the following default notations were accepted:
M - interarrival (or service) times are iid exponential
D - interarrival (or service) times are deterministic
G - interarrival (service) times are iid and governed by some general distribution
Often, only the first three labels are present. Then the remaining parameters assume their default values, which are a general service discipline, infinite number of customers in the system, infinite customer population.
Example:
M/M/3/FCFS/7 means, iid exponential interarrival arrival times, iid exponential service times, 3 parallel servers, first-come first-served discipline, maximum 7 customers in the system, which includes the customers in service and waiting in the line. This system has limited queue length with 4 waiting positions.
What is M/M/1?
Performance measures
Average number of jobs in the system \(L\)
Average number of jobs in the waiting line \(L_q\)
Average number of jobs in service \(L_s\)
Average time spent in the system \(W\)
Average time spent in the waiting line \(W_q\)
Average time spent in service \(W_s\)
M/M/1 queueing system
Back to M/M/1 queueing system.
- a single-server,
the interarrival times \(T\) are iid exponential with parameter \(\lambda\),
the service times \(S\) are iid exponential with parameter \(\mu\).
Let us consider a particular M/M/1 system with \[\lambda = 3|\mbox{hr}, \hspace{1cm} \mbox{and} \hspace{1cm} \mu = 5|\mbox{hr}.\] We need to check that \(\rho = \frac{\lambda}{\mu} <1\) to make sure that the system will operate properly (i.e., the number of arrivals in the system will not increase to infinity). In technical terms, this is the condition that will assure that the system will reach its steady state performance. \(\rho\) is called the traffic intensity.
Also, from the Queueing Theory is known that the expected number of jobs (or arrivals) \(L\) in a properly operating M/M/1 system is
\[ L = \frac{\rho}{1 - \rho}.\]
M/M/1 - the Simulation Model (with \(\lambda = 3\), \(\mu = 5\))
the state of the system at time \(t\) is \(N(t)\) - the number of jobs in the system
- Possible states \[ \{0,1, 2, \ldots, n, \ldots\}\]
Events - an arrival or a completion of service The event generation method will be described later
The state transition formula is (departure is equivalent to service completion)
- There are two methods for advancing the simulation clock
Fixed small time increment \(\Delta t\) (not considered in SCIE201)
Advance time by a small fixed amount equal to \(\Delta t\)
Update the system: determine what events have occurred during the elapsed time interval \(\Delta t\), what is the resulting state of the system.
Next-event increment
Advance time to the next event occurring in the system
Update the system by determining the new state and by randomly generating the time till the next event
We will be discussing the next-event increment scenario for M/M/1 system. We proceed as follows:
for the event arrival - generate a uniform random number \(r_A \in (0,1)\) and convert it to an \(exp(\lambda)\) observation.
for the event departure - generate a uniform random number \(r_D \in (0,1)\) and convert it to an \(exp(\mu)\) observation.
So, next question is: given a uniform random number \(r\), how to generate observations from an exponential distribution with parameter \(\lambda\)? It can be shown that:
is an observation from \(X \sim exp(\lambda)\).
For convenience of the next computations, we need to convert the arrival rate \(\lambda = 3\) per hour to rate per minute, i.e., \(\lambda = \frac{3}{60} = \frac{1}{20}\) per minute. Similarly, the service rate \(\mu = \frac{5}{60} = \frac{1}{12}\) per minute.
So, we need some uniform random numbers \(r_A\) and \(r_D\) to generate observations from the interarraval time \(T\) and the service time \(S\) of the M/M/1 system. We will use the following sequence of uniform random numbers:
\[ r_A =\{ 0.096, 0.569, 0.764, 0.492, 0.950, 0.610, 0.145\}\] and \[r_D = \{0.665, 0.842, 0.224, 0.552, 0.590, 0.041 \}\]
Next, we will generate a sequence of the interarrival times using the sequences of the uniform random numbers given above:
\(r_{A_1} = 0.096\) then \(T_1 = -\frac{1}{\frac{1}{20}}\,\, log (1-0.096) = 2.019\).
\(r_{A_2} = 0.569\) then \(T_2 = -\frac{1}{\frac{1}{20}}\,\, log (1-0.569) = 16.833\).
\(r_{A_3} = 0.764\) then \(T_3 = -\frac{1}{\frac{1}{20}}\,\, log (1-0.764) = 28.878\).
and so on …
Next we will generate a sequence of the departure times:
\(r_{D_1} = 0.665\) then \(S_1 = -\frac{1}{\frac{1}{20}}\,\, log (1-0.665) = 13.123\).
\(r_{D_2}= 0.842\) then \(S_2 = -\frac{1}{\frac{1}{20}}\,\, log (1-0.842) = 22.142\).
and so on …
Use the uniform random numbers given in the Table below to simulate this M/M/1 system (with \(\lambda = 3\), \(\mu = 5\)).
| t, min | N(t) | \(r_A\) | Next IAT | \(r_D\) | Next ST | Next A | Next D | Next event |
|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 0.096 | 2.019 | - | - | 2.019 | - | A |
| 2.019 | 1 | 0.569 | 16.833 | 0.665 | 13.123 | 18.852 | 15.142 | D |
| 15.142 | 0 | - | - | - | - | 18.852 | - | A |
| 18.852 | 1 | 0.764 | 28.878 | 0.842 | 22.142 | 47.730 | 40.994 | D |
| 40.994 | 0 | - | - | - | - | 47.730 | - | A |
| 47.730 | 1 |
A possible R code to simulate M/M/1 queue - next-event approach.
# Discrete queue simulation - next-event increment
mm1next<-function(T){ # T the time to run the simulation
### mm1 simulation using next-event method
### return list of: t - times of system state change and q - number of jobs in system
rateARR<-3.0 # arrival rate
rateSER<-5.0 # service rate
clock <- 0.0 # timer (simulation clock)
n<-0 # number of jobs in the system
q<-c(0.0) # list of number of jobs in the system
t<-c(clock) # related times
nextA <- clock+rexp(1,rate=rateARR)
nextD <- nextA + rexp(1,rate=rateSER)
while (clock <= T){
clock <- min(nextA,nextD)
if (clock==nextA){ # the clock is at arrival time
n<-n+1 # increase the number of jobs by 1
cat(clock,n,' Arrival \n')
nextA <- clock+rexp(1,rateARR) # next arrival time
if (n == 1) {nextD=clock+rexp(1,rateSER)} # if only one job next departure time
}
else{ # the clock is at departure time
if (n>0){
n<-n-1
cat(clock,n,' Depart\n')
if (n >=1) nextD<-clock+rexp(1,rateSER) # there are jobs waiting for service
# and next service starts immediately
else nextD <- 99999 # n = 0, no waiting jobs -
# not sure when the next job will be served and depart -
# assign something very large to nextD
}
else cat("Error: ",clock,n,nextA,nextD,'\n') # n == 0 the clock can not be at departure time
}
q<-append(q,n)
t<-append(t,clock)
}
return(list(t=t,q=q))
}
set.seed(121)
reslt <- mm1next(3.0)## 0.4300122 1 Arrival
## 0.4639046 2 Arrival
## 0.7054606 1 Depart
## 0.7498767 2 Arrival
## 0.8478262 1 Depart
## 0.8924036 2 Arrival
## 0.9706018 3 Arrival
## 1.183617 4 Arrival
## 1.282534 3 Depart
## 1.285379 2 Depart
## 1.31546 1 Depart
## 1.353962 2 Arrival
## 1.5295 1 Depart
## 1.643032 2 Arrival
## 1.696698 1 Depart
## 1.779501 0 Depart
## 2.122024 1 Arrival
## 2.184859 0 Depart
## 3.057353 1 Arrivaldf<-data.frame(reslt)
print(df)## t q
## 1 0.0000000 0
## 2 0.4300122 1
## 3 0.4639046 2
## 4 0.7054606 1
## 5 0.7498767 2
## 6 0.8478262 1
## 7 0.8924036 2
## 8 0.9706018 3
## 9 1.1836172 4
## 10 1.2825336 3
## 11 1.2853794 2
## 12 1.3154596 1
## 13 1.3539616 2
## 14 1.5294995 1
## 15 1.6430315 2
## 16 1.6966979 1
## 17 1.7795010 0
## 18 2.1220239 1
## 19 2.1848592 0
## 20 3.0573526 1plot(reslt$t,reslt$q, type="s", ylab="Number in the system", xlab = "Time")
Computing the average number of jobs in the system - the idea of time average
The command mean() is useful for estimating the average of a set of isolated observations such as waiting times or times in the system for different arrivals/jobs. These are independently measured values, which exist only when waiting finishes or an arrival leaves the system.
But many variables are not of this form: for example the number of customers in a bank. These “continuous values” have the characteristic that they exist at all times – though their value may be zero. They are “continuous” in time though they may be discrete in value (such as the number of customers).
One way to find the average level of a variable such as the average length of a queue or the average number in the system (which is the same as the average state of the system, which is the same as the average level of the system) is to observe the value at regular intervals and take the mean of those observations. This will certainly give an unbiased estimate of the mean number of jobs in the system. On the other hand we know that these values can only change at event times (as illustrated in the graph above). But since we have all the times and corresponding state of the system, we can reconstruct the complete time-graph. Thus we can calculate the average state of the system precisely.
This is the idea of time average, so we need the time-weighted average of the graph above. This is the total area under the graph divided by the total time.
# Discrete queue simulation - next-event increment
mm1next<-function(T){ # T the time to run the simulation
### mm1 simulation using next-event method
### return list of: t - times of system state change and q - number of jobs in system
rateARR<-3.0 # arrival rate
rateSER<-5.0 # service rate
clock <- 0.0 # timer (simulation clock)
n<-0 # number of jobs in the system
q<-c(0.0) # list of number of jobs in the system
t<-c(clock) # related times
nextA <- clock+rexp(1,rate=rateARR)
nextD <- nextA + rexp(1,rate=rateSER)
while (clock <= T){
clock <- min(nextA,nextD)
if (clock==nextA){ # the clock is at arrival time
n<-n+1 # increase the number of jobs by 1
#cat(clock,n,' Arrival \n')
nextA <- clock+rexp(1,rateARR) # next arrival time
if (n == 1) {nextD=clock+rexp(1,rateSER)} # if only one job next departure time
}
else{ # the clock is at departure time
if (n>0){
n<-n-1
#cat(clock,n,' Depart\n')
if (n >=1) nextD<-clock+rexp(1,rateSER) # there are jobs waiting for service
# and next service starts immediately
else nextD <- 99999 # n = 0, no waiting jobs -
# not sure when the next job will be served and depart -
# assign something very large to nextD
}
else cat("Error: ",clock,n,nextA,nextD,'\n') # n == 0 the clock can not be at departure time
}
q<-append(q,n)
t<-append(t,clock)
}
return(list(t=t,q=q))
}
set.seed(121)
# computing the average number of jobs in the system - the so called time Average.
df<-data.frame(mm1next(5000.0))
df1<-data.frame(t1=diff(df$t), q1=head(df$q, -1))
df1$new <- df1$t1*df1$q1
L <- sum(df1$new)/tail(df$t, 1)
print(L)## [1] 1.494495