in Number system versus time:
one simulation, exponential
service rate =60
average time in system 0.20435
average number in system 5.0243
Varying the service times
Compare exponential service times, constant
service times and uniform service times, all
with the same mean. Which appears to give
the smaller queue length?
Question 13
Consider one-server queues, all with exactly
the same expected interarrival time and
expected service times. Which system will
have the shortest expected queue length?
A: M/M/1
B: Poisson arrivals, uniformly distributed service
times
C: Poisson arrivals, constant service times
D: Constant interarrival times, constant service
times
E: constant interarrival times, exponential service
times
Possible Queue Disciplines:
Queue discipline:
FIFO=first in first out.
SIRO=service in random
order
LIFO=last in, first out
SPT= shortest
processing time first
PR=priority queue
(certain arrival types are
designated higher
priority and serviced
first- may bump lower
priority individuals in
service)
Queues. The A/S/c/N/K
notation
First symbol: A= arrival
time distribution is one
of
M=Poisson (Markov)
arrivals =exponential
interarrival times
D=constant
(Deterministic)
interarrival times
G=General interarrival
time distribution
Second symbol: S indicates
service time distribution
M=exponential service
times.
D=constant
G=general.
Third symbol:
c=number of servers
(=1,2,infinity)
Fourth: N=system capacity
(e.g. size of waiting room).
Fifth: K= calling population
size (may be infinite)
Two server queues: series and
parallel systems
Examples: The M/M/2 Queue: Java code by Lev Givon
(Columbia)
http://www.columbia.edu/~leg22/e6711/markovapplet.html
Series or Tandem Queue
Two queues in series.
Variables for tandem queue
(Ross, 6.3)
Counter variables:
NA=number of arrivals by time t
ND=number of departures by time t
Output Variables
A1 (n)=arrival time of customer n
A2(n)=arrival time of customer n it server 2
D(n)=departure time of customer n
Event List
tA=time of next arrival, t1=time of next service completion at server
1, t2 =time of next service completion at server 2.
SS variables
n1,n2 = number of customers at servers 1,2 (including in service)
Algorithm for tandem queue
If tA<min(t1,t2) the next
event is an arrival
Update clock t to tA
Increment NA
Generate new arrival time
tA
Increment n1
If server 1 started empty,
generate service time for
arrived customer
Collect output
If t1<min(tA,t2) the next
event is a departure from
server 1
Update clock t to t1
decrement n1
If n1=0, set t1= ,
otherwise generate service
time completion for next
customer at server 1
Collect output
Algorithm for tandem queue
(cont)
If t2<min(tA,t1) the next event is a
departure from server 2
Update clock t to t2
decrement n2
If n2=0, set t1= 1, otherwise generate service
time completion for next customer at server 2
Collect output
Parallel Servers
Example: a repair facility.
Two mechanics are hired to
maintain four robot paint
machines used in the
manufacture of
automobiles. The paint
machines each have
exponential times between
breakdowns. Average time
between breakdown=1 hour
for each machine.
For each server, time to
repair a machine is
exponential, mean=1/2
hour.
(a) What is the
equilibrium distribution?
(b) What is the long-run
average number of
inoperable machines?
Transition diagram for a
M/M/2/4/4 repair facility.
Nodes labeled with number
of broken machines
Verifying and validating the
simulation
Two Questions:
is the simulation
representing the model
correctly (model
verification)
does the model
accurately represent real
system? (model
validation)
The steady-state
probabilities can be
used to validate the
simulation. Compare
the empirical frequency
of the values 0,1,2,3,4
to their long-run
probabilities and the
average number in the
system to its expected
value.
A simulation of the repair facility
function [ET,NB,time_in_state]=repair_facility(N0,T)
t=0; Nt=N0;
%Nt=current number broken between 0 and 4
ET=[0]; NB=[N0];
%initialize event times ET and number broken NB vector.
lam=[4 3 2 1 .00000000001]; mu=[.00000000001 2 4 4 4];
% replaced lambda=0 by lambda very small for coding convenience
while t<T
TB=-log(rand)/lam(Nt+1); % time till next breakdown
TR=-log(rand)/mu(Nt+1); % time till next repair
t=t+min(TR,TB);
%update clock
ET=[ET t];
%update vector of event times
Nt=Nt+sign(TR-TB);
%number broken after time t.
NB=[NB Nt];
% update the number broken vector at time t.
end
d=diff(ET); time_in_state=[];
for n=length(d); i=0:4
time_in_state=[time_in_state sum(d(NB(1:n)==i))];
end
time_in_state=time_in_state/sum(time_in_state); % estimates probability of
states
Compare theoretical and
observed times in states
[ET,NB,TS]=repair_facility(0,30);
[ET,NB,TS]=repair_facility(0,5000);
TS =
0.1883
0.3714
0.2728
0.1313
0.0363
True steady-state values= [1 2 3/2 3/4 3/16]*16/87
= 0.1839 0.3678 0.2759 0.1379 0.0345
How are these steady state probabilities obtained?
Determining the steady state
distribution of the number in system
Let n be the probability of n broken machines
n = the rate of transition from n to n + 1
n = rate of transition from n to n-1.
Then rate leaving state n = rate entering state n.
n (n + n ) = n 1n 1 + n +1 n +1 , for n = 1,2,..
0 0 = 11
i
=1
i
Determining the steady state distribution of
the number in system (cont)
Solve
these
equations
for i .
Obtain
n
=
where
i
0
0 1 ... n 1
,n 1
1 2 ... n
i
=1
0
is determined
by
Steady state distribution of the number in
M/M/2/4/4
1
=
0
3
=
0
4
=
0
where
queue
4
2
4
2
4
2
0
3
4
3
4
2
=
0
43
24
2
4
21
44
is determined
by
i
= 1 so
i
16
[1
87
Expected
=
16
87
2
3/2
number
(2 + 3 + 9/4
+ 3/4)
3/4
inoperable
=
116
87
3/16]
=
= 1 . 33
Simulating the reliability of a
system: e.g. power grid.
Montreal Power Grid
Montreal receives power if
1,2,7 or
1,6,5,8 or
1,2,3,4,8 or
1,6,5,4,3,7 are all
operating. Suppose times
until failure of each of these
components under certain
circumstance (e.g. ice
storm) independent
System still operating
at time t if
min(T1 , T2 , T7 ) > t or
min(T1 , T6 , T5 , T8 ) > t or
min(T1 , T2 , T3 , T4 , T8 ) > t
min(T1 , T6 , T5 , T4 , T3 , T7 ) > t
where Ti = failure time of link i.
source
sink
SYSTEM FAILURETIME =
max(min(T1 , T2 , T7 ), min(T1 , T6 , T5 , T8 ), min(T1 , T2 , T3 , T4 , T8 ), min(T1 , T6 , T5 , T4 , T3 , T7 ))
Simulation of Montreal power grid
function L=lifepowr(T)
% input vector T of length 8,
lifetimes of components of montreal
% power grid. Outputs lifetime L of
the system.
L1=min(T([1 2 7]));
L2=min(T([1 5 6 8]));
L3=min(T([1 2 3 4 8]))
L4=min(T([1 3 4 5 6 7]));
L=max([L1,L2,L3,L4]);
Run Montreal power grid
simulation.
Assume exponential lifetimes, mean=20
years. (Matlab code)
L=[];
for i=1:10000
L=[L lifepowr(exprnd(20,1,8))];
end
hist(L,50)
(or try L=[L lifepowr(unifrnd(0,40,1,8))]; )
(scriptlifepowr)
Results:
Exponential lifetimes:
min(L) =0.0185
mean(L)= 9.1255
Uniform lifetimes:
min=0.018349 mean=
13.0282
System has Longer expected
lifetime
What if there is some dependence
among the exponential random
variables?
I used what is called a Normal copula model to produce the
dependence: i.e. I generated correlated () normal random
variables Xi and put
Zi = - ln( (Xi)) to get correlated exponential Zi with mean .
When =0.9 and = 20
L=[]; x=corexprnd(20,.9,10000,8);
for i=1:10000
L=[L lifepowr(x(i,:))];
end
hist(L,50)
This (copula) method for generating dependence is very common in
credit risk management
Results: now a 10% change
life<1.3
Although mean lifetime goes up, so
does the probability lifetime is <2
years
Up to around 0.13.
Number of simulations for
given accuracy
95 % CI for parameter is mean 2 S /
n
where S 2 is sample variance obtained from pilot simulation .
To ensure estimator within ,
(with confidence
set = 2 S /
around 95%)
n
2S
Solve for n . n
2
95% confidence interval for
mean system life
M=mean(L)
S=sqrt(var(L))
M+[-1 1]*2*S/sqrt(length(L))
to achieve accuracy to 2 decimals, need
delta=.01
sample size=(2*S/.01)^2
TOO LARGE-- we need VARIANCE
REDUCTION (better simulations)
Simulate the expected lifetime
of network below
Note flow is unidirectional on link 3
Simulating Lifetimes
Assume lifetimes are uniform[0,20]
L=[];
for i=1:10000
T=20*rand(1,5);
L=[L max([min(T([1 5])) min(T([2 3 5])) min(T([2 4]))])];
end
mean(L)
sqrt(var(L)/length(L))
hist(L,50)
Number of simulations for
given accuracy.
95 % CI for parameter is mean 2 S /
n
where S 2 is sample variance obtained from pilot simulation .
To ensure estimator within ,
(with confidence
set = 2 S /
around 95%)
n
2S
Solve for n . n
2
95% confidence interval for
mean system life
M=mean(L)
S=sqrt(var(L))
M+[-1 1]*2*S/sqrt(length(L))
to achieve accuracy to 2 decimals, need
delta=.01
sample size=(2*S/.01)^2
TOO LARGE-- we need VARIANCE
REDUCTION (better simulations)