Transcript Design , Layout, & Location of Facilities

Facilities Design
S.S. Heragu
Industrial Engineering Department
University of Louisville
Appendix:
Introduction to Queuing, Queuing
Network, and Simulation Modeling
A queuing system
Input Process
Queue
Service Process
CCCC
S
S
S
C: customers; S: servers
Departure
Process
Queuing Models are descriptive
models
•
•
•
•
What is the expected number of parts waiting
in a queue?
What is the expected time a part spends
waiting in a queue?
What is the probability that a machine will be
idle?
What is the probability of a queue being filled
to capacity?
Elements of a queuing system
-
Arrival process
Service process
Departure process
Queue discipline
Balking
Modeling of the Arrivals and Service Processes
•
•
•
Ti - ith interarrival time (or service time for the
ith customer)
Assume Ti to be an independent, continuous
random variable F.
Let the probability density function (pdf),
expected value and variance of F be f(t), E(F)
and V(F), respectively. Then,
c

0
c
P( F  c)   f (t )dt and P( F  c)   f (t )dt
Modeling of the Arrivals and Service Processes
•
•
Suppose that F follows an exponential
distribution with parameter λ.
The pdf f(t) for an exponential distribution is
a strictly decreasing function of t (t 0). So,
P0  F   t  Pt  F  t   t for any t and  t  0
•
It is not difficult to infer from this property
that F is likely to take on a very small value,
perhaps near zero.
Modeling of the Arrivals and Service Processes
•
•
•
•
Let F represent service time
Time to serve a customer is typically small,
but an occasional customer requires
extensive service
This service time distribution can be
approximated as an exponential distribution
Many, but not all systems exhibit such a
characteristic
Modeling of the Arrivals and Service Processes
•
PDF of an exponential distribution
e  t for t  0
f (t )  
0 for t  0
•
Cumulative probabilities for an exponential
distribution
c
c
0
0
P( F  c)   f (t )dt   et dt  1  ec
Modeling of the Arrivals and Service Processes
•
Verify the previous result by substituting x =
-λt. Then, dt = -(1/λ)dx. Therefore,
c
c
 e dt   e dx  e
 t
x
0
x c
0
 e
t c
0
 1  e c
0

P( F  c)   e dt  e
 t
 t 
c
 e  c
c

E ( F )   tf (t )dt
c
•

V ( F )   t  E (t ) f (t )dt
2
c
Mean can be shown to be 1/λ and Variance
can be shown to be 1/λ2 for exponential dist
Property 1 of the Exponential Distribution
• Memoryless
P(F  t   t F  t )  P(F   t ), provided t,  t  0
P( F  t   t F  t ) P( F  t   t )
P( F  t   t F  t ) 

P( F  t )
P( F  t )
•But
P( F  t )  et and P(F  t   t )  e (t  t )
•Hence,
P( F  t   t ) e  (t  t )
P( F  t   t F  t ) 
  t  P( F   t )
P( F  t )
e
Property 2 of the Exponential Distribution
• If interarrival time is exponentially distributed,
then the arrival process is Poisson
•Consider an arrival process {Nt = n}, where Nt =
n is the number of arrivals occur during time
interval t (t 0).
•Assume that N0 = 0 and the arrival process
satisfies the following three conditions
Property 2 of the Exponential Distribution
•Assume that N0 = 0 and the arrival process
satisfies the following three conditions.
-(1) The probability of an arrival occurring
between times t and  t depends only upon the
length and does not depend upon either the
number of arrivals occurring until time t or the
specific value of t. The probability is equal to
 t  o( t ), where o( t ) is a small quantity in
comparison to  t especially as tends to zero.
In other words,
o( t )
lim
0
 t 0  t
Property 2 of the Exponential Distribution
•Assume that N0 = 0 and the arrival process
satisfies the following three conditions.
- (2) The probability of > 1 arrival occurring
during a very small time  t interval = o( t )
- (3) The number of arrivals occurring in
nonoverlapping intervals are independent.
E.g., the number of arrivals between the first
twenty time units has no bearing on the
 t between the next twenty,
number of arrivals
thirty or forty time units.
Property 2 of the Exponential Distribution
•Assume P{Nt = n} is the probability of Nt = n
arrivals occurring during time interval t
•Probability of –ve arrivals is zero. Then,
P  Nt  t  n  P  Nt  n
lim
   P  Nt  n  P  Nt  n  1 for n  0,1, 2,...
 t 0
t
•LHS above is nothing but differentiation of P{Nt
= n} with respect to t. These infinite linear, first
order differential equations can be solved and
P{Nt = n} can be shown to be
 t 
n
e
 t
/ n ! for n=0,1,2,...
Property 2 of the Exponential Distribution
•Given the arrival process described by P(Nt = n),
assume F is the random variable for the time
between successive arrivals – the interarrival time
•To show that F follows an exponential
distribution, note that P{t T} = P{NT = 0}; i.e., the
probability of the next arrival exceeding time any
specific value T is equivalent to having no arrival
in time T.
P{t  T }   f (t )dt and P{N (t )  0}  e 
•Taking derivatives of both sides, we get:
f(t) = λe-λT, which is an exponential distribution

0
 t
Property 2 of the Exponential Distribution
•
•
•
•
Mean of the Poisson distribution given by E{Nt}
= λt is equal to the variance
Exponential is the only distribution to have
such a property
Much of the above discussion holds for the
service process also. To verify, the reader has
to replace the terms interarrival time and arrival
with service time and service completion,
respectively
Notice that the notation used for the arrival
process parameter is λ, whereas for the service
process it is μ
Kendall-Lee notation for Queuing models
•
A / B / C / D / E / F, where
- A denotes the nature of the arrival process
- B denotes the service time distribution
- C is the number of parallel servers (S 1)
- D is the queue discipline (FCFS, LCFS,
SIRO, GD, etc.)
- E is the max number of customers allowed
in the system (C 1)
- F denotes the size of the calling
population (finite or infinite)
Queuing models that can be solved exactly
•
•
•
•
•
The calling population is infinite
Some finite population models, like the
M/M/1/GD/K/K model - typically called the
machine repairman model
The queue capacity is infinite
For some special models, analytic solution is
possible for only FCFS and priority
disciplines
The probability distribution of the interarrival
and service times are exponential
Birth-death queuing model
•
•
The birth-death process is a continuous time
Markov chain or stochastic process in which
the state of the system at any continuous
time (rather than discrete points in time) is a
nonnegative integer
Consider a system at some state at time t, in
which we have Nt = n customers. Assume
that the pdf of the remaining time until next
arrival (birth) and pdf of the remaining time
until next service completion (death or
departure) are exponential with parameters
λn, μn, respectively, for n = 0, 1, 2, …
Birth-death queuing model
•
•
•
•
Only one event (either a birth or a death) can
occur in any period of time Δt. Birth and
death rates are independent of each other,
but state dependent
Because of the exponential assumption, birth
and death occur randomly
Figure above (rate diagram) shows an infinite
birth-death process in which birth and death
rates are state dependent
Analysis of most birth-death processes
computationally feasible only when the
system has reached steady state
Birth-death queuing model
•
•
•
•
•
Mean rate of entering any state n = Mean rate
of leaving that state
Balance equation or flow conservation
equation
Denote En(t), Ln(t) as the number of times the
process enters and leaves state n by time t
We have |En(t)-Ln(t)| < 1
Dividing the left- and right-hand sides by t
and letting t -> infinity, we get
E (t )
L (t )
. Hence,
En (t ) Ln (t )
lim n
 lim n
lim
t 
t

t
0
t 
t
t 
t
Birth-death queuing model
•
•
•
•
•
•
•
Pj is steady-state probability of being in state j
For state 0, μ1P1 = λ0P0. Therefore, P1 = (λ0/μ1)P0
For state 1, μ2P2 + λ0P0 = λ1P1 + μ1P1 Therefore,
- P2 = (λ1/μ2)P1 + (1/μ2)(μ1P1 - λ0P0) =
(λ1/μ2)(λ0/μ1)P0 + (λ0/μ2)(P0 - P0) = (λ1/μ2)(λ0/μ1)P0
For state 2, μ3P3 + λ1P1 = λ2P2 + μ2P2 Therefore,
- P3 = (λ2/μ3)P2 + (1/μ3)(μ2P2 - λ1P1) =
(λ2/μ3)(λ1/μ2)(λ0/μ1)P0
………
Let cn = (λn-1/μn)(λn-2/μn-1)…(λ0/μ1) for n = 1, 2, …
Then, Pn = cnP0, for n=1, 2, ...
Birth-death queuing model

P
n 0
n
 1.

1


Pn   Pn  1. So, 1   cn  P0  1 and P0 
.



n 1
 n1 
1   cn 
 n1 
L =
Lq =
Ls =
W =
W q=
Ws =
λ =
μ =

average number of customers in system
average number of customers in queue
average number of customers in service
average time a customers spends in system
average time a customer spends in queue
average time a customer spends in service
mean arrival rate
mean service rate
Birth-death queuing model

L   nPn
n 0
•

Lq   (n  S ) Pn
nS
Ls  L  Lq
Using Little’s Law, which states that L = λW, we
can calculate W, Wq and Ws
M/M/1/GD/inf/inf queuing model
•
•
•
•
•
•
•
•
•
Let λ and μ be the arrival and service rate
If λ > μ, the queue will grow infinitely; ρ = λ/μ is called the
traffic intensity or serve utilization
ρ must be less than 1
λj = λ; μj = μ, for i=1, 2, ... and μ0 = 0
P1 = (λ/μ)P0, P2 = (λ/μ)2P0, ..., Pn = (λ/μ)nP0. Substituting ρ =
(λ/μ) and assuming 0 < ρ < 1, we get P1 = ρP0, P2 = ρ2P0, ...,
Pn = ρnP0, …
P0+P 1+…+P n+… = 1. So, P0(1 + ρ + ρ 2 + … + ρn + …) = 1
To evaluate (1 + ρ + ρ 2 + … + ρn + …), let S = (1 + ρ + ρ 2 + …
+ ρn + …). Then ρS = (ρ + ρ 2 + … + ρn + …)
Clearly, ρS - S = 1. Hence, S = 1/(1-ρ). Thus, P0(1+ρ+ρ 2+
…+ρn +…) = P0/(1-ρ) = 1. Therefore, P0 = 1-ρ
Similarly, P1 = ρ(1-ρ), P2 = ρ2(1-ρ), …, Pn = ρn(1-ρ), …
M/M/1/GD/inf/inf queuing model



L =  n Pn   n  (1   )  (1   )  n n
n
•
n=0
n 0
n=0
The summation can be determined by substituting

S  =  n n    2 2  3 3  ...
n 0
•
•
•
Then S    S  =    2  ...   n  ...
RHS of the above = ρS and S was shown to be = 1/(1-ρ)
Hence, we get S    S  =  and S  =  . Therefore,
2
1   
1   
L
•
(1   ) 

( /  )




2
(1   )
(1   ) (1   /  ) (    )
Lq can also be determined in this way:



n 1
n 1
n 1
Lq   (n  1) Pn   nPn   Pn  L  (1  P0 )  L  1  1    L  
Lq 

(   )
 


      2
2


(   ) 
 (   )
 (   )

 1
 
2

Ls  L  Lq 



 
(   )  (   )
 (   )  (   )  
M/M/1/GD/inf/inf queuing model
•
Using Little's formula, we can then determine, W, Wq and Ws
as
1

W 
Wq 
 (   )
 (   )
L
Ws 
1

Example 1
Find the following values:
• 1) The percentage of the time that the press is idle.
• 2) The average number of parts in the queuing system.
• 3) The average queue length.
• 4) The throughput time of the system.
• 5) The amount of time spent in the queue.
Example 1 Solution
10
•
•
•
•
•
•
•
•
0
10
1
10
2
10
3
M/M/1 queuing system
12
12
12
12
λ = 10 per hour and μ = 12 per hour
Utilization factor is ρ = λ/μ = 0.833
The probability that the press is idle is the probability that
there are no jobs being served or waiting. P0 = 1-ρ = 0.167 or
the press is idle 16.7% of the time
The average number of parts in the queuing system is L =
ρ/(1-ρ) = 5 parts. This is also the average Work-In-Process
(WIP) inventory
The average queue length is Lq = ρ2/(1-ρ) = 4.167 parts
The throughput time (or average time spent by a part in the
system) is W = 1/(μ(1-ρ)) = 0.5 hours
The amount of time spent in the queue is Wq = ρ/(μ(1-ρ)) =
0.4167 hours = 25 minutes
. . .
Property 3 of the Exponential Distribution
•The minimum of independent exponential
random variables is also exponential
•Let F1, F2, …, Fn be independent exponential
random variables with parameters μ1, μ2, … μn
•Let Fmin = min {F1, F2, …, Fn}. Then for any t > 0,
•P(Fmin > t) = P{F1>t}P{F2>t}…P{Fn>t} =
e
1 t
e
2 t
 e n t = e[ 1 +  2 ++  n ]t .
Property 3 of the Exponential Distribution
•Suppose n types of customers, with the ith type of
customer having an exponential interarrival time
distribution with parameter λi, arrive at a queuing system
•Assume that an arrival has just taken place
•Then, from the no-memory property, it follows that the time
remaining until the next arrival is also exponential
•Using property 3, it is easy to see that the interarrival time
for the entire queuing system (which is the minimum
amongst all interarrival times) has an exponential
distribution with parameter S
 i .
i 1
M/M/S/GD/inf/inf queuing model
•
•
•
•
•
•
•
Let λ and μ be the arrival and service rate
If n < S, all customers are in service
If n > S, all S servers are busy and n-S customers are waiting
From property 3, service rate for the entire queuing system
is
.
 utilization
The
factor for an M/M/S system is ρ = λ/Sμ
Clearly, μn = μ{min{n,S}}
Hence, μn = nμ for n = 1, 2, …, S and μn = Sμ for n = S+1, S+2,
…, Also, λn = λ for n = 0, 1, 2, …, Hence, Pn = P0cn, for n = 1,
2, … S. We get
S
i 1
i
n
 
  P  
Pn  P0  
   n   for n  1, 2,..., S
1 2 n  n !   
M/M/S/GD/inf/inf queuing model
S
 
    P0      
Pn  P0  



     
1 2 S  S  S   S !     S  
P0 
1

c
n 0
•
n
•
1

S
c  c
n 0
n
nS
n

n
 P0    
=
 for n  S  1, S  2,..., 
nS  
 S !S    
1
( /  ) n (  /  ) S


n!
S!
n 0
S 1
  



nS  S  

nS
The second summation in the denominator is of the form
1+x+x2+x3+… where x = λ/(Sμ). This summation is equal to
1/(1-x) or 1/(1-(λ/Sμ)). Hence,
P0 
•

nS
1
( /  ) (  /  ) S


n!
S!
n 0
S 1
n


1
 1   /( S  ) 


Substituting ρS=(λ/μ), we get
Thus,
P0 
1
S 1
(  S )n (  S ) S


n!
S!
n 0
 P0 (  S )n
for n  1, 2,..., S

n!
P0  
n
 P0 (  S ) for n  S  1, S  2,...
 S ! S n  S
 1 
 1  


M/M/S/GD/inf/inf queuing model
•
We now find Lq


nS
i 0
Lq   (n  S ) Pn   iPS i
•
•
•
•
•
•
(  S )S  i P0 (  S )S P0  i (  S ) S P0 
 i

 i  S !(1   )2
S!
S ! i 0
i 0

The last summation above denoted as S’ in the M/M/1 model was
previously shown to be equal to ρ/(1-ρ)2
Wq = Lq/λ, W=Wq+1/μ, L=λW = λ(Wq+1/μ)= Lq+λ/μ
Results for the basic M/M/1 and M/M/S models are summarized in
the Table A.1
It should be noted that ρ=λ/μ for the single server model, but
ρ=λ/(μS) for the multiple server case. In both models, Pn=cnP0.
Formulae for cn are summarized below
Single server model: cn=ρn, for n=1,2,…, P0 = 1-ρ
Multiple server model:
 1   n
   for n  1, 2,..., S
 n!   
cn  
n
 1 
 S ! S n  S    for n  S  1, S  2,..., 
 

P0 
1
S 1
( /  ) n (  /  ) S


n
!
S!
n 0


1
 1   /( S  ) 


Shortcut formula - M/M/S/GD/inf/inf
model
•
•
Note that the P0 calculation which itself is required to calculate Lq
and Wq is somewhat tedious
Sakasegawa (1977) came up with an approximate formula to
calculate Wq from which Lq can be obtained using Little’s law. The
Wq formula is shown next

 / S 
Wq 
S 1   / S 
2 ( S 1) 1
Property 4 of the Exponential Distribution
•If we have multiple arrivals in a queuing system
and each interarrival time has an exponential
distribution, the interarrival time for the system as
a whole is also exponential with parameter equal
to the sum of the parameters for each interarrival
•Reverse also true
Example 2
In an effort to reduce WIP inventory, the press in Example 1 is
modified to have two servers. Find the following:
• The average queue length.
• The average number of parts in the queuing system (WIP).
• The expected amount of time spent in the queue.
• The throughput of the system.
Example 2 Solution
10
0
1
12
•
•
•
•
•
•
•
•
10
10
2
24
10
. . .
3
24
24
M/M/2 model: Rate diagram shown above
λ = 10 per hour and μ = 12 per hour. The 2 servers result in a
utilization factor of ρ = λ/Sμ = 0.417
The new idle time of the press is P0 = 1/[(ρS)S/S!(1-ρ) +
Σ((ρS)n/n!)] = 7/17 = 0.412, or 41.2% of the time
The expected queue length in the two server system is Lq =
[ρ(ρS)SP0]/[S!(1-ρ)2] = 7/40 = 0.175 parts
The expected WIP is L = Lq + λ/μ = 121/120 = 1.0083 parts
The expected amount of time spent in the queue is Wq =
[(ρS)SP0]/[S!Sμ(1-ρ)2] = 7/400 = 0.0175 hours = 1.05 minutes
The throughput time of the system is W = Wq + 1/μ = 0.10083
hours = 6.05 minutes
Veriy results using Sakasegawa (1977) formula
Example 3
•
•
•
•
A workstation in a company receives parts from two sources - a subcontractor
and a machining center within the company. The arrival rate of each is Poisson
with parameter λ/2. The company has purchased another identical workstation
in order to increase throughput. Each has a service rate of μ. The company has
the following three ways of designing the system.
- Have all the parts from the subcontractor visit the first workstation, and all
the parts from the internal machining center visit the recently purchased
workstation for processing
- Because the parts coming from the subcontractor are identical to those
from the internal machining center, and the two workstations can be
combined into one so as to have a service rate of 2μ, the company can
have parts from both sources join a single queue and visit the combined
workstation
- Combine the two arrivals into a single queue but keep the workstations
separate
Determine the effective arrival and service rates of the above systems
Determine which system would result in the most waiting time for the parts.
Which system would result in the least waiting time?
If reducing WIP becomes a more important criterion that waiting time in the
system, show that system (c) is preferred to system (b) above. Assume that
WIP is measured as the mean number of parts in the queue
Example 3 Solution
•
System (a) is an M/M/1 queue with parameters λ/2 and μ. System (b)
is an M/M/1 queue with parameters λ and 2μ. System (c) is an M/M/2
queue with parameters λ and μ. The three systems are depicted
below.
Example 3 Solution
•
•
•
•
•
The arrival rate is λ/2 + λ/2 = λ for (b) and (c). The effective arrival and service
rates are λ and 2μ for all three systems.
For an M/M/1 system with parameters λ and μ, the waiting time in a queue (W)
and length of the queue (Lq) are: W =1/(μ(1-ρ)) = 1/(μ-λ); Lq=λ2/μ(μ-λ).
For an M/M/2 system, P0 =1/[(2ρ)2/(2!(1-ρ)) + Σ1n=0 ((2ρ)n/n!)] =1/[(2ρ2)/(2!(1-ρ) +
(2ρ)0/0! + (2ρ)1/1!] =[2(1-ρ)]/[2+4ρ-2ρ-4ρ2+4ρ2)] = (1-ρ)/(1+ρ); Lq=[(2ρ)2P0ρ]/[2!(1ρ)2] = [(2ρ)2(1-ρ)ρ]/[2!(1+ρ)(1-ρ)2] = 2ρ3/(1-ρ2) = λ3/[μ(4μ2-λ2)]; Wq=Lq/λ = λ2/[μ(4μ2λ2)]; W=Wq + 1/μ = λ2/[μ(4μ2-λ2)]+ 1/μ = 4μ/(4μ2-λ2)
For system (a), therefore, W(a) = 1/(μ-λ/2) = 2/(2μ-λ). Lq(a) = λ2/4μ(μ-λ/2) = λ2/2μ(2μλ). For system (b), W(b) = 1/(2μ-λ). Lq(b) = λ2/2μ(2μ-λ). For system (c), W(c) =
4μ/(4μ2-λ2). Lq(c) = λ3/[μ(4μ2-λ2)]. To prevent an infinite queue, λ must be less
than 2μ (i.e. ρ<1). Therefore, W(c) = 4μ/(4μ2-λ2) = [4μ(2μ+λ)][1/(2μ-λ)] > 1/(2μ-λ) =
W(b). Thus, W(a)>W(c) >W(b). Hence, system (a) would result in the most waiting
time and system (b) the least.
The WIP in the queue of system (c) is given by Lq(c) = [λ2/(2μ(2μ-λ))][2λ/(2μ+λ)] <
[λ2]/[2μ(2μ-λ)] = Lq(b). When WIP in the queue is the criterion system (c) is
preferred to system (b), but when manufacturing lead time is considered, (b) is
better than (c).
Finite Queue Size Models
•
•
Single server model:
  n for n = 1,2,...,C
cn = 
 0 for n  C
 ( S )n /n!

n
nS
c n =  [( S ) ] / [ S! S ]
0

Multiple server model:
C 1
eff    Pn   (1  PC )
for n  1,2,..., S
for n  S  1, S  2,...,C
for n  C
n 0
Model
M/M/1/GD/C/inf
M/M/S/GD/C/inf
Po
(1-ρ)/(1-ρC+1)
1/[ΣSn=1(ρS)n/n! + ((ρS)S/S!)(ΣCn=S+1 ρn-S]
L
ρ/(1-ρ) - (C+1)ρC+1/(1-ρC+1)
Lq+λ(1-PC)/μ
Lq
L+P0-1
(ρS)SP0ρ[1-ρC-S-(C-S)ρC-S(1-ρ)]/[S!(1-ρ)2]
W
L/(λ(1-PC))
L/(λ(1-PC))
Wq
Lq/(λ(1-PC))
Lq/(λ(1-PC))
Finite Source Models
•
•
Single server model:
 N ! n
for n  0,1, 2,..., N

cn   ( N  n)!
0
for n  N

Multiple server model:
Model

N ! n
 ( N  n)!n !  n for n  0,1, 2,..., S


N ! n
cn  
for n  S  1, S  2,..., N
nS n
 ( N  n)! S ! S 
0
for n  N


M/M/1/GD/N/N
M/M/S/GD/N/N
P0
1/[ΣNn=0(λnN!/(μn(N-n)!)]
1/[ΣS-1n=0(λnN!)/(μnn!(N-n)!)+ΣNn=S(λnN!)/(μnS!Sn-S(N-n)!)]
L
N-μ(1-P0)/λ
Lq+ΣS-1n=0nPn + S(1-ΣS-1n=0Pn)
Lq
N-(λ+μ)(1-P0)/λ
ΣNn=S(n-S)Pn
W
L/(λ(N-L))
L/(λ(N-L))
Wq
Lq/(λ(N-L))
Lq/(λ(N-L))
Nonexponential Models
•
•
L
M/G/1:
M/Ek/1:
Lq
2
2
2
2 2
2
2
2
2
 1 +k 
λσ + ρ
λ /(k μ ) + λ / μ
λ +λ k
=
=
=
=

2( 1  ρ)
2( 1  λ / )
2k (μ  λ)  2k 
Wq =
•
G/G/m:
 2 2   2
2(1   )
 1+ k 
=

  2k 
Lq
 ca2 + cs2 

W q  
2







 (    ) 




.
 (    ) 
2


λ


 μ(μ  ) 
Queuing Networks
•
•
•
Open
Closed
Semi-Open
Open Networks
•
•
•
•
Jackson network has the following properties.
The arrival process from the external node to node i is Poisson with
mean rate of γi.
The service time at node i follows an exponential distribution with
parameter μi.
pi0 is the probability that a part will exit the network after completion
of processing at node i and pij is the probability that a part will visit
node j after completion of processing at node i. pi0 and pij are
assumed to be known and independent of the state of the system.
λ
Machine 1
µ1
Machine 2
µ2
Two-machine System
λ
Two-node, Open Jackson Network
n1-1, n2+1
n1+1, n2
The flow balance (rate in equals rate out) equations
for any {n1, n2} pair greater than 0
(λ+μ1+μ2)P(n1,n2) = λP(n1-1,n2) + μ2P(n1,n2+1) + n1-1, n2
μ1P(n1+1,n2-1)
For the remaining states (0,0), (n1,0), (0,n2), the
following flow equations are also easy to derive.
λP(0,0) = μ2P(0,1)
(λ+μ1)P(n1,0) = μ2P(n1,1) + λP(n1-1,0)
(λ+μ2)P(0,n2) = μ1P(1,n2-1) + μ2P(0,n2+1)
Also, Σ Σ P(n1,n2) = 1
n1, n2-1
n1, n2
n1, n2+1
n1+1, n2-1
Two-node, Open Jackson Network
n1-1, n2+1
n1+1, n2
n1, n2
P(n1,n2) = ρ1n1(1-ρ1)ρ2n2(1-ρ2) = Pn1Pn2
P(0,0) = (1-ρ1)(1-ρ2) = P0P0
n1, n2-1
n1-1, n2
n1, n2+1
n1+1, n2-1
General Open Jackson Network
•
•
•
•
•
•
•
For each node i, i=1, 2, …, m, let:
γi
mean rate of the Poisson arrival process from the external
world to node i
1/μi
mean (exponentially distributed) service time at node i
Si
number of identical servers at node i
pi0
probability that a part will exit the network after completion of
processing at node i
pij
probability that a part will visit node j after completion of
processing at node i.
Because the arrival at a node is made up of external and internal
arrivals - the latter occurring as a result of service completions at
other internal nodes - the arrival rate at node i, i=1, 2, …, m, is:
m
i   i   p ji  j
j 1
Example 4
•
Consider a system consisting of a pre-heat furnace, a forge press, and
an inspection station. An ingot arrives at the furnace to be heated and
then sent to the press. After pressing, the ingot is transferred to
inspection. Ingots are delivered to the furnace at a rate of 12 per hour.
The furnace heats the ingots at a rate of 15 per hour. The press and
inspection each have capacities of 18 and 21 ingots per hour. All
interarrivals are assumed to be exponential. The three machine network
is shown below.
λ = 12
•
•
Furnace
Press
Inspect
µ = 15
µ = 18
µ = 21
Determine the expected number of jobs in the furnace (F), press (P), and
inspection (I) as well as the WIP and processing time of the system.
Suppose the capacity of the queue in front of the press is 5 ingots, is
the current design sufficient?
Solution to Example 4
•
•
•
•
•
•
•
•
•
•
•
•
•
•
λ = 12
Furnace
Press
Inspect
µ = 15
µ = 18
µ = 21
The arrival rate for both P and I is λ = 12 ingots per hour
λF = 12 per hour μF = 15 per hour
ρF = 0.80
λP = 12 per hour μP = 18 per hour
ρP = 0.67
λI = 12 per hour μI = 21 per hour
ρI = 0.57
Using the equations in Table A.1, we find the operating characteristics
as follows.
WF = 1/(μF(1-ρF)) = 1/(15(1 - 0.8)) = 0.33 hours = 20 minutes
LF = ρF/(1-ρF) = 0.8/(1 - 0.8) = 4 ingots
WP = 1/(μP(1-ρP)) = 1/(18(1 - 0.67)) = 0.167 hours = 10 minutes
LP = ρP/(1-ρP) = 0.67/(1 - 0.67) = 2 ingots
WI = 1/(μI(1-ρI)) = 1/(21(1 - 0.57)) = 0.11 hours = 6.67 minutes
LI = ρI/(1-ρI) = 0.57/(1 - 0.57) = 1.33 ingots
The system WIP and processing time is the sum of the individual
values.
L = LF + LP + LI = 4 + 2 + 1.33 = 7.33 ingots
W = WF + WP + WI = 20 + 10 + 6.67 = 36.67 minutes
Solution to Example 4
•
•
•
To find out if a design with a maximum of 5 ingots is sufficient, we must
find the probability that more than 5 ingots are in the queue. Recall that
Pn, the probability of n jobs being at a machine in steady state is given
by the equation, Pn=(1-ρ)ρn, n 1. Therefore,
Pn>5 = 1-Pn<5 = 1-(1-0.67)(1.00+0.67+0.672+0.673+0.674+0.6755) = 0.09046.
This design will result in the queue in front of the press being at full
capacity 9.05% of the time. When the queue is full there will be a
blocking problem in the queuing system and the steady-state results in
part 1 will not hold because they were derived under the assumption
that each machine has an infinite queue in front of it.
λ = 12
Furnace
Press
Inspect
µ = 15
µ = 18
µ = 21
Example 5
•
Consider the open network queuing system in Example 4. Suppose that
the rejection rate at inspection is 10%. Rejected parts are sent back to the
furnace and are reworked through the entire system. This new network is
shown in Figure A.11. Determine the expected number of jobs in the
furnace (F), press (P), and inspection (I) as well as the WIP and
processing time of the system.
λR
λ = 12
Furnace
Press
Inspect
µ = 15
µ = 18
µ = 21
λA
Solution to Example 5
Step 1: Calculate the arrival rate at each center.
• λF=γF+λR, where λF = effective arrival rate to the furnace
• γF = external arrival rate
• λR = internal arrival rate (rejected parts) = P(R)λF
• We also have,
• P(R) = the probability of rejection;
• P(A) = the probability of acceptance = 1-P(R)
• Therefore,
• λF= γF/(1-P(R))= γF/P(A)
• The arrival and service rates and utilization at each center are:
• λF=13.33 per hour
μF=15 per hour
ρF=0.889
• λP=13.33 per hour
μP=18 per hour
ρP=0.741
• λI=13.33 per hour μI=21 per hour
ρI=0.635
• Utilization factors are now calculated using the effective arrival rates
Solution to Example 5
Step 2: Analyze each machine in the system independently.
• WF=1/(μF(1-ρF)) = 1/(15(1 - 0.889)) = 0.6 hours = 36 minutes
• LF=ρF/(1-ρF) = 0.889/(1 - 0.889) = 8 ingots
• WP=1/(μP(1-ρP)) = 1/(18(1 - 0.741)) = 0.214 hours = 12.9 minutes
• LP=ρP/(1-ρP) = 0.741/(1 - 0.741) = 2.86 ingots
• WI=1/(μI(1-ρI)) = 1/(21(1 - 0.635)) = 0.124 hours = 7.5 minutes
• LI=ρI/(1-ρI) = 0.635/(1 - 0.635) = 1.74 ingots
Step 3: Combine the results from each center to analyze the performance of
the entire system.
• vF=λF/γF = 13.33/12 = 1.111 vP=vI = 1.111
• L=LF+LP+LI = 8 + 2.86 + 1.74 = 12.6 ingots
• W=vFWF+vPWP +vIWI = 1.111(36 + 12.9 + 7.5) = 62.7 minutes
• The "rework" loop increases the WIP and cycle time by 71.9%, and
70.9%, respectively
Example 6
•
The figure below shows a four-machine open network queuing system
consisting of a furnace (F), a press (P), a rolling mill (M), and an
inspection area (I). Each machine has a single server and the rejection
rate is 5%. The part routing matrix is shown below. Given that the external
arrival rate is λ=12, and μF=15, μP=12, μM=18, and μI=21, determine the WIP
and processing time of the system
Press
Furnace
Inspect
Mill
From/To
F
P
M
I
Exit
F
-
0.3
0.7
-
-
P
-
-
-
1.0
-
M
-
0.4
-
0.6
-
I
0.05
--
-
-
0.95
Solution to Example 6
Step 1: Calculate the effective arrival rate.
• λF=λ+0.05λI; λP=0.3λF+0.4λM; λM=0.7λF; λI=1.0λP+0.6λM=λF
• From Example 5, we know that λF = λ/P(A) = 12/0.95 = 12.63. Solving these
equations yields the effective arrival rates presented next. The service rates
and utilizations at each machine center are also shown.
• λF=12.63
μF=15
ρF=0.842
• λP=7.33
μP=12
ρP=0.611
• λM=8.84
μM =18
ρM=0.491
• λI=12.63
μI=21
ρI=0.602
Step 2: Analyze each machine in the system independently. All the machines
have only one server, therefore, we will use Table A.1 for the operating
characteristics of an M/M/1 queue.
• WF=1/(μF(1-ρF)) = 1/(15(1 - 0.842)) = 0.422 hours = 25.33 minutes
• LF=ρF/(1-ρF) = 0.842/(1 - 0.842) = 5.33 ingots
• WP=1/(μP(1-ρP)) = 1/(15(1 - 0.842)) = 0.422 hours = 25.33 minutes
• LP=ρP/(1-ρP) = 0.611/(1 - 0.611) = 1.57 ingots
• WM=1/(μM(1-ρM)) = 1/(18(1 - 0.491)) = 0.109 hours = 6.6 minutes
• LM=ρM/(1-ρM) = 0.491(1 - 0.491) = 0.97 ingots
• WI=1/(μI(1-ρI)) = 1/(21(1 - 0.602)) = 0.119 hours = 7.2 minutes
• LI=ρI/(1-ρI) = 0.602(1 - 0.602) = 1.51 ingots
Solve two key equations simultaneously
•
•
•
•
Parametric Decomposition (P-D) Method
- Two parameters – mean and scv
- Can be applied for any general, but known distribution
Decompose network into stochastically independent queues
- Capture network interaction effects
- Solve two sets of equations simultaneously
Solve each node separately
- Apply G/G/m model
Combine results to get network-specific performance measures
Key Network Operations
•
Departures
a
d  a
d
cd2  1  (1   2 )(ca2  1) 
2
departure
•
Merging of Arrivals
i
 a   i
a
i

 
c  w  i
i   k
 k
2
a
m
(cs2  1)

 2
ci  1  w


superposition
•
Splitting of Departures
i
a
i  pi a
ci2  pi c 2  1  pi
split
n
•
a 2   1, 2 d1
 jk 
oi
b
i
i
jl1 jkl
Y
i1 l 2
n oi
b Y
i
kl
Batching
i1 l1
2
2
ca2
 max0,1,2 1 1,2cd1

i
jkl
Key Network Operations
•
•
p
n
i 1
j 1
ˆk   0i k   ˆ j p jk  jk
Mean arrival rate into each node:
m
2
2
SCV of inter-arrival times into each node: cak  ak   cai bik , k  1,2,...,m
i 1

2
1  p jk   p jk  2j x j  


p

c

1

p




0
k
0
k
0
k
jk
jk

where


j 1


ak  1  wk  n

 p max   1, 0

jk
jk
 

j 0
n
bjk  wk pjk p jk  jk (1   2j )
uk 
1
  p 
n
j 0
2
jk
, p 'jk 
ij
and
j
x j  1
p
 jk 
oi
 ˆ
i 1 l 1
p oi
i
jk
 ijk,lY jli Yki,l 1
 ˆ
i 1 l 1
i
jk
Y jli Yki,l 1
1
max csj2 ,0.2  1

mj 
Closed Queuing Network
Queue 1
Queue
Robot
Machine 1
Queue 1
Machine 1
.
.
.
.
.
.
Queue k
Machine k
Queue 2
Machine 2
Queue k + 1 Machine k + 1
Queue
Robot
.
.
.
.
.
.
Queue m
Machine m
Mean Value Analysis (MVA) Algorithm
Step 1: The first step is to estimate an initial value of Lij
Step 2: Determine the throughput time Wij using the equation
Wij 
 N j  1  Lij



i j  N j  Si i j
1

 L
    ir r
 r  j  Si i

, i  1,2,...,m; j  1,2...,r

Step 3: Determine the production rate Xj for each machine j
Xr 
Nj
m
v W
i 1
ij
ij
Step 4: Determine the queue length Lij using the equation Lij=Xj(vijWij)
Step 5: Compare the Lij value calculated in step 4 with the previous value. If
the new value is within a desired range of the previous value, stop. If not,
go to step 2.
Example 7
•
•
•
•
Consider a machine shop which manufactures two parts. Part 1 visits
machines 1 and 3, and part 2 visits machines 1 and 2. The service rates for
machines 1, 2, and 3 are 10 parts per hour, 15 parts per hour, and 12 parts
per hour, respectively. The service rates are not affected by the part being
serviced. JIT/Kanban considerations limit the WIP of parts 1 and 2 to 5 and
4 units respectively.
1) Determine the amount of time each part spends waiting and being served
at each machine.
2) Determine the number of each part waiting and being served at each
machine.
3) Based on the solution to parts 1 and 2, what conclusions can be drawn
about the system design?
Solution to Example 7
•
•
•
•
Step 1: Determine initial values of Lij. The number of each part type in the
system is N1=5, and N2=4. Each part type visits only 2 of the 3 machines.
Therefore, L11=L31 =N1/2=2.5 and L12=L22=N2/2=2. Because part type 1 does not
visit machine 2, and part type 2 does not visit machine 3, L21=L32=0
Step 2: Determine the throughput time Wij
- W11 = 1/10 + [4/5][2.5/10] + 2/10 = 0.50 = 30 minutes
- W21 = 0
- W31 = 1/12 + [4/5][2.5/12] = 0.25 = 15 minutes
- W12 = 1/10 + [3/4][2/10] + 2.5/10 = 0.50 = 30 minutes
- W22 = 1/12 + [3/4][2/15] = 0.25 = 0.167 = 10 minutes
- W32 = 0
Step 3: Determine the production rate Xr. In this example, vij =1 for all i and r. - We get X1 = 5/[30+15] = 1/9 = 0.111 parts per minute
- X2 = 4/[30+10] = 1/10 = 0.10 parts per minute
Step 4: Determine the queue length Lij using Lij=Xr(vijWij) for all i and j.
- L11 = (0.111)(30) = 3.33 parts; L21 = 0; L31 = (0.111)(15) = 1.67 parts
- L12 = (0.10)(30) = 3 parts; L22 = (0.10)(10) = 1 part; L32 = 0
- Notice that L11+L31=5, and L12+L22=4
Solution to Example 7
•
•
Step 5: Compare the value of Lij calculated in step 4 with the previous value.
The values of Lij warrant another iteration. Therefore, we return to step 2 and
repeat until the new calculated values of Lij are within the recommended
range (5%) of the previous values. The solution below is found after four
additional iterations.
- L11 = 4.28 parts
W11 = 47.8 minutes
- L21 = 0
W21 = 0
- L31 = 0.72 parts
W31 = 8.1 minutes
- L12 = 3.60 parts
W12 = 47.5 minutes
- L22 = 0.40 parts
W22 = 5.3 minutes
- L32 = 0
W32 = 0
The current design is creating a build up in front of machine 1. Approximately
88% [(4.28+3.60)/(5+4)] of the WIP is either waiting or being served at machine
1.