Transcript Matrix Analytic Methods

Matrix Analytic Methods Some Real Life Applications
David Lucantoni
DLT Consulting, L.L.C
www.DLTconsulting.com
Queueing Theory and Network Applications, Korea University, 2011
Objectives of this Talk
 Review Matrix Analytic Methods pioneered by Marcel Neuts
and others
 Present real life examples where Matrix Analytic Methods
have been used successfully in industry
 Congestion due to Internet Dial Up
 Application to Common Channel Signaling (SS7)
 Modeling WiFi Hotspot Traffic
 Multiplexing bursty traffic in an ATM or IP network
 Call Acceptance for a VoIP System
 Transient Results and applications
Queueing Theory and Network Applications, Korea University, 2011
Objectives of this Talk (Cont’d)
 No new mathematical analyses
 A detailed description of any of these models would take the
whole talk
 Will only give general description of model parameters
 Detailed definitions and analyses are widely available in the literature
 Will give some matrix results along with the corresponding scalar
(exponential or Poisson) results
Queueing Theory and Network Applications, Korea University, 2011
Analysis of Matrix Analytic Models
 Models are analyzed by purely probabilistic arguments (e.g.,
Markov Renewal Theory)
 Explicit expressions lead to very stable numerical algorithms for
computation of performance metrics
 Matrix analytic methods avoid
 the complex analysis used in previous approaches
 the necessity of numerically computing the roots of an analytic
function
 Numerical instability when roots are close to each other or multiple
Queueing Theory and Network Applications, Korea University, 2011
Phase-Type (PH) Distributions
 Will only consider continuous-time distributions here
 Phase-type distributions contain many well-known and
popular distributions
 Exponential distribution, (m=1)
 Erlang distribution
 Hyper-exponential distribution
Queueing Theory and Network Applications, Korea University, 2011
Phase-Type Distributions (Cont’d)
 Defined as the time till absorption in an (m+1)-state Markov
process with 1 absorbing state
 Let T be the m×m matrix of infinitesimal rates between
transient states, α, the vector of probabilities of starting in a
particular transient state and T0 be the rates into the
absorbing state
 Let X be a phase-type random variable and F(x)=P(X≤x)
 Then
F(x) = 1- e- l x (exponential) F(x) = 1- a eTxT 0 (Phase - type)
Queueing Theory and Network Applications, Korea University, 2011
Assuring Emergency Services Access (e.g.,
911) in the presence of long holding times
due to Internet dial-up calls
 The Problem






No dial-tone condition - “Customers hear thin air.”
VERY SERIOUS – e.g., 90/day life threatening emergency calls .
Priorities non-implementable
Preliminary suspicion on maintenance activities.
Could it be congestion ? Red flag – internet usage !
Model solved using a finite quasi-birth-and-death process (QBD)
 [Ramaswami, et al., 2005]
Queueing Theory and Network Applications, Korea University, 2011
Cable Telephony
TDM/T-1
ADM
HDT
GR-303
Backbone
Transport OC-48+
SS7
DACS
GR-303
ADM
DWDM
Proprietary
Protocol
CLASS 5
SWITCH
CMTS
DWDM
Master Headend or
Backbone
PSTN
Primary Hub
(100K - 200K HHP) Each
Network
Interface
Unit
EMS
Systems
DOCSIS
Protocol
DHCP/DNS
TFTP
Servers
Cable
Modem
Fiber Node
Transport
Central
Office
Zone 4
Zone 5
 16 Per Secondary
Hub
Zone 1
Zone 2
Zone 3
Queueing Theory and Network Applications, Korea University, 2011
Effect of Long Holding Times (Cont’d)
 Preliminary Analysis
 problem in evening hours 7 PM – 11 PM
 no discernible spatial pattern
 Push-backs
 Engineering followed “standard” procedures
 Only 5-8% of calls are internet dial-ups
 Field measurements “do not support” congestion hypothesis
Queueing Theory and Network Applications, Korea University, 2011
Phase Type Model
(using the EM-algorithm) ; fit uses ~4M data points
Phase Type Fit
Probability Observed
0.1
5.4
0.2
12
0.3
21
0.4
32.4
0.5
48
0.6
72.6
0.7
120.6
0.8
232.8
0.9
601.8
0.95
1237.8
0.99
3952.2
0.995
6074.4
Fitted
5.8
12.7
21.2
32.2
47.7
72.2
120.2
235.9
597.6
1268.0
4035.2
7168.4
Queueing Theory and Network Applications, Korea University, 2011
The Length Biasing Effect
 Remaining holding time distribution: [1-F(x)] / μ
RESIDUAL DISTRIBUTION
h(x)=[1-F(x)]/μ
Probability Observed Residual
0.1
5.4
40.8
0.2
12
116.1
0.3
21
238.0
0.4
32.4
423.7
0.5
48
703.0
0.6
72.6
1124.6
0.7
120.6
1809.7
0.8
232.8
3285.3
0.9
601.8
7028.2
0.95
1237.8 11073.8
0.99
3952.2 20489.2
Median
Mean
48
297
703.0
2367.0
Queueing Theory and Network Applications, Korea University, 2011
Empirical Validation
 Empirical data on residuals (~11 K observations) vs. Model
Queueing Theory and Network Applications, Korea University, 2011
Impact of Long Holding Times
 Insensitivity relates only to blocking over an infinite horizon.
 Heavy tails
 PERSISTENCE of congestion.
So, short term performance is quite different.
“Standard” procedures do not work !
 Bad periods compensated by long very good periods. Good
situation for control !!! (See the paper and patents for
possible controls.)
Queueing Theory and Network Applications, Korea University, 2011
Several Other Real Life
Applications
 Excessive Link Oscillations in the Common Channel
Signaling (SS7) Network [Ramaswami and Wang,
1993]
 Solved using Phase-type distributions
 Modeling and characterization of large-scale Wi-Fi
traffic in public hot-spots [Ghosh, et al., 2011]
 Solved using an M/G/∞ where G is ln PH distribution
Queueing Theory and Network Applications, Korea University, 2011
The Batch Markovian Arrival Process
(BMAP)
 Natural generalization of the Poisson process
 Contains many well-known arrival processes (in all cases,
correlated batch sizes are easily handled)
 Poisson process (m=1)
 PH renewal process
 Markov modulated Poisson process (MMPP)
 Superposition of BMAPs is again a BMAP
 Useful for modeling arrivals to an internal node in a network
Queueing Theory and Network Applications, Korea University, 2011
The Batch Markovian Arrival Process (Con’t)
 Consider an m-state Markov process and assume that at each
transition of this process, there is a probability of a batch arrival
that depends both on the state before and after the transition
 Let Dn, n≥1, denote the infinitesimal rate of a batch arrival of size
n, keeping track of the underlying states before and after the
transition
 The generating function D(z) is given by
D(z) = å n=0 Dn z n
¥
For Poisson : D(z) = - l + l z
Queueing Theory and Network Applications, Korea University, 2011
BMAP/G/1 Queue
 The BMAP/G/1 queue is extremely useful for modeling the performance of
broadband packet networks
 Traffic is bursty
 Service times are NOT phase-type
 Asynchronous Transfer Mode (ATM)
 Fixed length cells result in deterministic service times
 IP networks
 Finite packet sizes  finite mixture of deterministic service times
 BMAP models very bursty traffic and is closed under superposition, i.e.,
modeling the output of nodes entering another node
Queueing Theory and Network Applications, Korea University, 2011
Busy Period of the BMAP/G/1 Queue
 Let h(s) be the Laplace-Stieltjes transform of the
service time distribution
 Let G(z) be the probability generating function of the
number of customers served during a busy period
(keeping track of the arrival phase). Then we have
G(z) = zh(l - lG(z))
(M / G /1)
G(z) = zh(-D[G(z)]) (BMAP / G /1)
Queueing Theory and Network Applications, Korea University, 2011
Busy Period of the BMAP/G/1 Queue
(Cont’d)
 Define G to be G(1) and g to be the stationary
probability vector of G
 It is easy to show that g is also the stationary
probability vector of the phase of the arrival process
during idle periods
 For M/G/1, the probability that the system is empty is 1ρ; for BMAP/G/1, the probability that the system is
empty and that the arrival process is in phase j is the jth
element of (1-ρ)g
Queueing Theory and Network Applications, Korea University, 2011
Virtual Waiting Time of the
BMAP/G/1 Queue
 Let W(s) be the Laplace-Stieltjes transform of the
virtual waiting time distribution (keeping track of the
arrival phase)
 Let g be the stationary probability vector of the
stochastic matrix G
 Then we have
s(1- r )
W (s) =
s - l + l h(s)
M / G /1
W (s) = s(1- r )g[sI + D(h(s))]-1
BMAP / G /1
Queueing Theory and Network Applications, Korea University, 2011
Stationary Queue Length at Departures
 Let Xn be the stationary probability that there are n
customers in the queue at departures
 Let X(z) =
n
,
X
z
å 0 n then
¥
(1- r )(z -1)A(z)
X(z) =
z - A(z)
M / G /1
X(z) = (1- r )l -1gD(z)A(z)[zI - A(z)]-1
BMAP / G /1
Queueing Theory and Network Applications, Korea University, 2011
Computing these Distributions
 We compute these distributions by numerically
inverting the transforms [Choudhury, Lucantoni and
Whitt, 1994]
 These algorithms are for inverting multi-dimensional
transforms which we use later for computing transient
distributions
Queueing Theory and Network Applications, Korea University, 2011
Known Results for the BMAP/G/1
Queue
 Stationary distributions [Lucantoni, 1993]
 Queue length distribution at arrivals, departures, arbitrary time
 Waiting time distributions
 Transient distributions (given the appropriate initial conditions)
[Lucantoni, Choudhury, and Whitt,1994], [Lucantoni, 1998]
 P(system empty at time t)
 Workload at time t
 Queue length at time t
 Delay of nth arrival
 Queue length at nth departure
 P(nth departure occurs less than or equal to time x)
Queueing Theory and Network Applications, Korea University, 2011
Modeling the Input to an Interior
Network Node
 Consider the output from several nodes which all go to
the same node
 Ideal candidate for the MAP/G/1 queue
Queueing Theory and Network Applications, Korea University, 2011
Multiplexing Bursty Traffic in an ATM
Network
 ATM requirement is that the probability of
blocking must be less than 10-9

Tails of distributions are the most
important performance measure
 Superposition of 64 on-off sources each
modeled as a 2-state MMPP
 Exact Analysis
 Far from Poisson
 Exponential tail does not become relevant
until blocking is 10-40
Queueing Theory and Network Applications, Korea University, 2011
Advantage of MAP/G/1 Approach over MatrixGeometric Approach When service time has
Many Phases
 We computed the tail of the waiting time distribution in an
MMPP64/E1024/1 queue by an by a distribution that is close to an
E1024 and matches the first three moments of the deterministic
distribution
 The matrix-geometric approach would require using matrices of
order 65,536
 The MAP/G/1 approach looks at the process embedded at
departures; since we don’t have to keep track of the phase of
service, the matrices are only of order 64
Queueing Theory and Network Applications, Korea University, 2011
Advantage of MAP/G/1 (Cont’d)
 G satisfies
¥
G=ò e
0
D[G]x
dH(x)
 When H is E1024, the integral on the right hand side is
computed by inverting a 64×64 matrix 10 times
 Much easier than computing the matrix exponential
Queueing Theory and Network Applications, Korea University, 2011
IP Data Traffic is Bursty
 Sample Measurements
4500
from a previous
assignment
 Highly bursty
 Mean = 720 pkts/sec
4000
3500
Packets
 Sampled every second
Packet Counts
5000
3000
2500
2000
1500
1000
500
0
0
1
2
Hours
Queueing Theory and Network Applications, Korea University, 2011
3
4
Simulated Poisson Traffic
Packet Counts - Poisson
800
780
 Poisson traffic looks
760
bursty
Packets
740
720
700
 However, note the scale
of the y-axis
680
660
640
0
5
10
Minutes
Queueing Theory and Network Applications, Korea University, 2011
15
Poisson Process Compared to
Data
Packet Counts - Poisson
3000
3000
2500
2500
2000
1500
2000
1500
1000
1000
500
500
0
0
5
10
Minutes
Packet Counts - Measurements
3500
Packets
Packets
3500
15
0
0
5
10
Minutes
Queueing Theory and Network Applications, Korea University, 2011
15
Fit an D-MMPP to Data
 [Heyman, Lucantoni,
2003]
 Choose λ1 so that the
peak is at the 95th
percentile of a Poisson
distribution, i.e.,
l1 + 2 l1 = peak
Then
l1 =
(
)
1+ peak -1
2
Queueing Theory and Network Applications, Korea University, 2011
Fit an MMPP to Data (Cont’d)
 Lower bound of data covered by λ1 is
l1 - 2 l1
 Let this be the upper bound covered by λ2, i.e.,
l2 + 2 l2 = l1 - 2 l1
 Then
l2 =
(
l1 - 2
)
2
 Continue…
Queueing Theory and Network Applications, Korea University, 2011
Fit an MMPP to Data (Cont’d)
 Compute the transition probabilities by
number _ of _ transitions _ from _ i _ to _ j
pij =
number _ of _ transitions _ out _ of _ i
 For performance computations convert to a continuous
time MMPP (see paper)
Queueing Theory and Network Applications, Korea University, 2011
Fit an MMPP to Data
 Q-Q plot: plots quantiles
 Should be linear if both
come from same
distribution
 Very good fit
Queueing Theory and Network Applications, Korea University, 2011
Simulated MMPP Traffic
Packet Counts - MMPP
3000
3000
2500
2500
2000
1500
2000
1500
1000
1000
500
500
0
0
5
10
Minutes
Packet Counts - Measurements
3500
Packets
Packets
3500
15
0
0
5
10
Minutes
Queueing Theory and Network Applications, Korea University, 2011
15
Another Example
Simulated MMPP Packet Counts
500
500
400
400
300
300
200
200
100
100
0
0
5
10
Minutes
Measurement Packet Counts
600
Packets
Packets
600
15
0
0
5
10
Minutes
Queueing Theory and Network Applications, Korea University, 2011
15
Call Acceptance for a VoIP System
 A startup company was offering VoIP service to clients by buying a fixed
bandwidth virtual circuit (“fat pipe”) through the Internet
 Wanted to know how big the fat pipe should be to handle a certain
number of voice calls and still meet the Quality of Service (QoS) e.g., tail
of the delay distribution
 Measurements were taken for IP voice calls for different codecs, e.g.,
G.723, G.729, etc., and different measurement intervals
 We developed a program where the user inputs the file containing the
measurements, the measurement interval and the target delay
Queueing Theory and Network Applications, Korea University, 2011
Call Acceptance for a VoIP System (Cont’d)
 We fit an MMPP to the measurement data and by solving
the delay using an MMPP/D/1 queue
 If the computed QoS is greater or less than the targeted
QoS, we adjust the service rate and re-compute; continue
adjusting the rate using a binary search until the computed
QoS is within a given tolerance
 The program then outputs the size of the required fat pipe
and the required buffer size
Queueing Theory and Network Applications, Korea University, 2011
Notes on Transient Results
 Generalized every result derived by Takács for the
transient M/G/1 queue [Takács, 1962]
 All results are direct matrix analogues of the scalar
results
 Using only probabilistic arguments, resulted in
expressions that can be computed numerically without
computing the roots of an analytic equation
Queueing Theory and Network Applications, Korea University, 2011
Sample Transient Result (emptiness function)
 Consider a MAP/G/1 queue (similar results hold for the
BMAP/G/1 queue)
 Let Px0(t) be an m×m matrix where the (i,j) entry is the
probability that the system is empty at time t with the arrival
process is in phase j given the at time 0 the work in the
system is x and the arrival process is in phase i
 Let px0(s) be the Laplace-Stieltjes transform of Px0(t)
Queueing Theory and Network Applications, Korea University, 2011
Emptiness Function (Cont’d)
 Then
-(s+ l - lG(s)) x
e
px 0 (s) =
s + l - lG(s)
M / G /1
px 0 (s) = e -(sI -D[G(s)])x (sI - D[G(s)])-1
MAP / G /1
Queueing Theory and Network Applications, Korea University, 2011
Transient Performance Measures
 We compute the virtual waiting time (work in system) at
time t for a given queue length at time 0 by numerically
inverting the 2-dimensional transform, W(s1,s2)
Queueing Theory and Network Applications, Korea University, 2011
Sample Transient Results Virtual Delay Distribution
 Superposition of
four identical twostate MMPP’s
 Service times are
distributed as E16
 Utilization = 0.7
 queue length = 0 at
time 0
 queue length = 32
at time 0
Queueing Theory and Network Applications, Korea University, 2011
Sample Transient Results Unstable System
 Utilization (ρ) = 2.0
 Queue length = 0 at
time 0
Queueing Theory and Network Applications, Korea University, 2011
Potential Use of Transient Results
 Call acceptance algorithms
 Current call acceptance algorithms are based on stationary
probabilities
 Given the number of calls in the system, the time to reach a
stationary distribution can be MUCH longer than the interarrival and inter-departure times
 Call acceptance tables could be computed a-priori based on
decision time intervals, the current state of the system, the
average and peak bandwidth and burstiness parameter of the
arriving call
Queueing Theory and Network Applications, Korea University, 2011
References
1. V. Ramaswami, D. Poole, S. Ahn, S. Byers, and A. Kaplan, "Assuring emergency
services access: .providing dial tone in the presence of long holding time internet dialup calls," Interfaces, vol. 35, pp. 411-22, 2005
2. V. Ramaswami and J. L. Wang, "Analysis of the Link Error Monitoring Protocols in the
Common Channel Signaling Network," IEEE/ACM Transactions on Networking, vol.
1, pp. 34-47, 1993
3. Ghosh, R. Jana, V. Ramaswami, J. Rowland, and N. K. Shankaranarayanan, "Modeling
and characterization of large-scale Wi-Fi traffic in public hot-spots," in IEEE
INFOCOM, Shanghai, China, 2011
4. D. M. Lucantoni, "The BMAP/G/1 queue: A tutorial," in Models and techniques for
Performance Evaluation of Computer and Communications Systems, L. D. a. R.
Nelson, Ed., ed: Springer Verlag, 1993, pp. 330-58
5. D. M. Lucantoni, G. L. Choudhury, and W. Whitt, "The transient BMAP/G/1 queue,"
Stoch. Mod., vol. 10, pp. 145-82, 1994
6. D. M. Lucantoni, "Further transient analysis of the BMAP/G/1 queue," Stoch. Mod.,
vol. 14, pp. 461-78, 1998
Queueing Theory and Network Applications, Korea University, 2011
References (Cont’d)
7.
Heymen and Lucantoni, “Modeling Multiple IP Traffic Streams with Rate Limits,”
IEEE/ACM Trans. On Networking, 11, No. 6, 2003
8.
Choudhury, Lucantoni and Whitt, “Multidimensional transform inversion with
applications to the transient M/G/1 queue," Ann. Appl. Prob ., 4, No. 3, 719-740,
1994
7.
L. Takács, Introduction to the Theory of Queues. New York: Oxford University
Press, 1962
Queueing Theory and Network Applications, Korea University, 2011
Where to Find Additional References
 Some Lucantoni papers can be downloaded from
www.DLTconsulting.com
 Follow the links to my Publications page (or just Google
my name)
 More references are listed in those papers
Queueing Theory and Network Applications, Korea University, 2011