$30 off During Our Annual Pro Sale. View Details »

Morphing M/M/m: A New View of an Old Queue

Morphing M/M/m: A New View of an Old Queue

In this centennial year of Erlang's 1917 paper that solved the M/M/m queue, we derive the residence time without resorting to the (by now) conventional machinery of probability theory and Markov chains. Addendum added June 27, 2018.

Dr. Neil Gunther

July 18, 2017
Tweet

More Decks by Dr. Neil Gunther

Other Decks in Research

Transcript

  1. Morphing M/M/m
    A new view of an old queue
    Dr. Neil Gunther — @DrQz
    Performance Dynamics
    July 3, 2018
    IFORS 2017 Conference, Quebec City, Canada
    Tuesday, July 18, 2017
    Session: TB-22, Room: 2104B
    Simulation, Stochastic programming, and Modeling Track
    SM
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 1 / 24

    View Slide

  2. Abstract
    TB-22 IFORS 2017 - Quebec City
    ⌅ TB-22
    Tuesday, 10:30-12:00 - 2104B
    Simulation, stochastic programming and
    modeling
    Stream: Simulation, stochastic programming and model-
    ing (contributed)
    Contributed session
    Chair: Benjamin Legros
    Chair: Dinesh Sharma
    1 - Mathematical analysis of machine repair problems with
    common cause failure, hot spares and multiple repair-
    men
    Dinesh Sharma
    We study the machine repairable system comprising M operating ma-
    chines, H spares and more than one repairman where "the partial server
    vacation" is applied on some of the repairmen. In this system, the first
    repairman never takes vacation and always available for servicing of
    failed machines while other repairmen goes to random length vacation
    whenever the number of failed machines are less than N, N +1 respec-
    tively. Machines may breakdown individually or due to common cause
    according to Poisson process. Vacation time and service time of repair-
    men follows the exponential distribution. Recursive approach is used
    to obtain the steady state probabilities. A cost model is developed to
    determine the optimum value of failed machine maintaining the sys-
    tem availability and other performance measures. Sensitivity analysis
    is investigated for optimal conditions and also analyzes the reliability
    characteristics of the system.
    2 - Unintended consequences of optimizing a queue disci-
    pline for a service level defined by a percentile of the
    waiting time
    Benjamin Legros
    4 - Single-period newsvendor problem under random end-
    of-season demand
    Subrata Mitra
    Newsvendor problems, which have attracted the attention of re-
    searchers since 1950’s, have wide applications in various indus-
    tries. There have been many extensions to the standard single-period
    newsvendor problem. In this paper, we consider the single-period,
    single-item and single-stage newsvendor problem under random end-
    of-season demand, and develop a model to determine the optimal order
    quantity and expected profit. We prove that the optimal order quantity
    and expected profit thus obtained are lower than their respective values
    obtained from the standard newsvendor formulation. We also provide
    numerical examples and perform sensitivity analyses to compute the
    extent of deviations of the ’true’ optimal solutions from the newsven-
    dor solutions. We observe that the deviations are most sensitive to the
    ratio of the means of the demand distributions. The deviations are also
    found sensitive to the contribution margin, salvage price, coe cients
    of variation of the demand distributions and correlation between sea-
    sonal and end-of-season demands. We provide broad guidelines for
    managers as to when the model developed in this paper should be used
    and when the standard newsvendor formulation would su ce to deter-
    mine the order quantity. Finally, we present the concluding remarks
    and directions for future research.
    ⌅ TB-23
    Tuesday, 10:30-12:00 - 2105
    MADM principles 2
    Stream: Multiple criteria decision analysis
    Invited session
    Chair: Jung-Ho Lu
    1 - A hybrid multiple attributes decision-making model for
    of the waiting time. This may create an incentive for managers to mod-
    ify the traditional first-come-first-served discipline of service. For this
    purpose, we consider the analysis of the M/M/s queue under the queue-
    ing discipline which minimizes a given percentile of the waiting time.
    We prove that a strict non-preemptive priority should be given to the
    oldest customer who has waited less than the acceptable waiting time.
    We derive closed-form expressions of the performance measures un-
    der this discipline, and evaluate the unintended consequences that this
    discipline may have on service levels and on sta ng decisions. In par-
    ticular, we show that although this discipline may reduce sta ng costs,
    it leads to excessive wait for non-prioritized customers.
    3 - Morphing M/M/m: A new view of an old queue
    Neil Gunther
    2017 is the centenary of A.K. Erlang’s paper on waiting times in an
    M/D/m queue. M/M/m queues are used to model call centers, multi-
    cores & the Internet. Unfortunately, those who should be using M/M/m
    models often don’t know applied probability theory. Our remedy de-
    fines a morphing approximation to M/M/m that’s accurate within 10%
    for typical applications+. The morphing residence-time formula is
    both simpler and more intuitive than the exact solution involving the
    Erlang-C function. We have also developed an animation of this mor-
    phing process. An outstanding challenge, however, has been to eluci-
    date the nature of the corrections that transform the approximate mor-
    phing solution to the exact Erlang solution. In this presentation, we
    show: 1) the morphing solutions correspond to the m-roots of unity
    in the complex z-plane; 2) the exact solutions can be expressed as
    a rational function with poles; 3) these poles lie inside the unit disk
    and converge around the Szego curve with increasing m-servers; 4)
    the correction factor for the morphing model is defined by the deflated
    polynomial; 5) the pattern of poles in the z-plane provides a conve-
    nient visualization of how the morphing solutions di↵er from the exact
    solutions.
    2
    88
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 2 / 24

    View Slide

  3. c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 3 / 24

    View Slide

  4. This can be a problem
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 3 / 24

    View Slide

  5. This can be a problem
    for the people who should be using it most
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 3 / 24

    View Slide

  6. This can be a problem
    for the people who should be using it most
    outside the OR community
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 3 / 24

    View Slide

  7. Motivation
    IT computer operators, administrators, developers, performance
    engineers, need to use queueing theory (OR)
    Many don’t have a background in applied probability theory (OR)
    Modeling often considered “unrealistic” by mgmt (OR)
    A significant issue for OR in general
    Simplest M/M/1 queue can be derived algebraically
    So can closed M/M/1/N/N (from Little’s law)
    Avoids probabilities, BD dynamics, Markov chains, etc.
    IT people like that
    How to derive M/M/m algebraically ?
    cf. Erlang’s approach of 1917
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 4 / 24

    View Slide

  8. Single Server M/M/1
    The “simplest” queue
    !
    S
    Mean residence time:
    R = S + W
    = S + QS
    = S + (λR)S
    = S + R(λS)
    = S + R ρ
    R =
    S
    1 − ρ
    (1)
    No probability theory used
    Mean waiting time
    W =
    ρS
    1 − ρ
    cf. A.K. Erlang 1909
    W =
    ρS
    2(1 − ρ)
    for “simpler” M/D/1 queue
    (Not 1st textbook example)
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 5 / 24

    View Slide

  9. Twin Server M/M/2
    Adding another server produces a less intuitive result
    !
    The Erlang solution is
    R2 =
    S
    1 − ρ2
    (2)
    Since ρ < 1, ρ2 << 1, so (1 − ρ2) > (1 − ρ)
    Hence, R2 < R1
    How to derive this algebraically?
    Think: ρ2 represents m = 2 server “dependency” of some kind
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 6 / 24

    View Slide

  10. Triple Server M/M/3
    Adding yet another server =⇒ triple dependency?
    !
    Guess m = 3 mean residence time is:
    R3 =
    S
    1 − ρ3
    Wrong! Correct expression is
    R3 = S +
    S
    3(1 − ρ)
    9ρ3
    2 + ρ(4 + 3ρ)
    (3)
    which is inscrutable and does not simplify
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 7 / 24

    View Slide

  11. Here’s What Happened
    Twin server case (2) was deceptively simple because of a coincidence
    Exact expression for R2
    — corresponding to (3) for R3
    — is:
    R2 = S +
    S
    2(1 − ρ)
    2ρ2
    1 + ρ
    (4)
    which just happens to simplify to (2)
    R2 =
    S
    1 − ρ2
    Last factor in (4) comes from the Erlang C function which is completely
    unintuitive and gets more complicated as m gets bigger.
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 8 / 24

    View Slide

  12. Erlang’s M/M/m Solution (1917)
    Exact equation for Rm
    is
    Rm = S +
    S
    m(1 − ρ)
    C(m, a) (5)
    where C(m, a) was published 1 by Agner Erlang exactly 100 years ago
    C(m, a) =
    am
    m!(1 − ρ)
    m−1
    k=0
    ak
    k!
    +
    am
    m!(1 − ρ)
    (6)
    a = mρ is the (telephonic) traffic intensity measured in Erlangs
    Arnold Allen: “Solving (6) is an unnatural act”
    Can code an algorithm but all intuition is lost
    1Tore O. Engset published later
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 9 / 24

    View Slide

  13. Miraculous Morphing
    Rm
    is captured intuitively in the following diagram and animation
    S/m
    !
    ?
    !
    m
    S
    A miracle happens
    Morphing queues
    As λ increases, M/M/m appears to “morph” from m parallel M/M/1
    queues into a single M/M/1 queue that is m-times faster
    Already know correct R1
    for M/M/1 from eqn. (1)
    Just need to formalize the “small miracle” in the middle
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 10 / 24

    View Slide

  14. M/M/2 Morphing Function
    Trial morphing factor that modifies R1
    for parallel M/M/1s (λ → 1
    2λ)
    R2 =
    S
    1 − 1
    2λS
    φ2(ρ)
    1 φ2(ρ) → 1 as 1
    2 λ → 0 (light traffic) for 2 parallel M/M/1 queues
    2 φ2(ρ) → 1
    2
    as 1
    2 λ → 1/S (heavy traffic) for 1 fast M/M/1 queue
    Suggests the functional form
    φ2(ρ) =
    1
    1 + ρ
    with ρ = λS/2
    R2 =
    S
    1 − ρ
    1
    1 + ρ
    =
    S
    1 − ρ2
    Correct! Pure algebra from diagram (no prob thy was used)
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 11 / 24

    View Slide

  15. Generalized Morphing Model
    Theorem 1 (Gunther 2004†)
    M/M/m morphing factor involves a truncated geometric series
    φm(ρ) =
    1
    1 + ρ + ρ2 + . . . + ρm−1
    φm(ρ) produces the approximate mean residence time for M/M/m

    m
    S
    1 − ρm
    (7)
    Intuitive derivation without probability theory
    Eqn. (7) is not identical to the Erlang solution
    Morphing Rφ
    m
    underestimates exact RE
    m

    Close enough for Guerrilla style operational estimates
    † See Analyzing Computer System Performance with Perl::PDQ, Springer (2004)
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 12 / 24

    View Slide

  16. Correcting the Morphing Approximation
    Morphing model and Erlang Rm
    identical for m = 1, 2
    Approximate for m ≥ 3
    Error generally < 10% 2
    Can also develop by playing The Multiserver Numbers Game (2011)
    Erlang Morphing
    m = {1, 2, 4, 8, 16, 128}
    0.0 0.2 0.4 0.6 0.8 1.0
    ρ
    1
    2
    3
    4
    R/S
    0
    0.05
    0.10
    0.15
    Can we find analytic corrections? Rφ
    m → RE
    m
    2
    TB-22 audience: Is there a bound on the morphing error? (See Addendum)
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 13 / 24

    View Slide

  17. Pieces of a Dream
    !
    "!
    (1#")!
    m
    !
    m
    !
    a = !S
    C(m, ρ) =
    B(m, ρ)
    1 − [1 − B(m, ρ)] ρ
    m
    S
    S/m
    !
    (1"#)!
    m
    S
    Many diagrams and simulations later it became clear that a visual/intuitive
    approach to finding morphing model corrections was doomed:
    C(m, ρ) affects Rm
    through waiting times Wm
    only
    C(m, ρ) arises from the Poisson dsn not a geometric series
    Correction terms in ρ are as complicated as C(m, ρ) itself!
    As ρ → 0, M/M/m behaves like an infinite server, not parallel queues
    Erlang’s 1917 paper contains no intuitive insights
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 14 / 24

    View Slide

  18. New Analytic Result
    Theorem 2 (Gunther 2016†)
    The missing terms correct the denominator of (7) as follows:

    m
    =
    S
    1 −
    cm
    Pm−1(ρ)
    ρm
    where Pm−1(ρ) is the deflated polynomial associated with
    Pm
    (ρ) = cm
    ρm + . . . + c3ρ3
    + c2ρ2
    + c1ρ + c0
    m Pm(ρ)
    1 −ρ + 1
    2 −ρ
    2
    + 1
    3 3ρ
    3
    + ρ
    2
    − 2ρ − 2
    4 8ρ
    4
    + 4ρ
    3
    − 3ρ
    2
    − 6ρ − 3
    5 125ρ
    5
    + 75ρ
    4
    − 20ρ
    3
    − 84ρ
    2
    − 72ρ − 24
    6 −54ρ
    6
    − 36ρ
    5
    + 30ρ
    3
    + 35ρ
    2
    + 20ρ + 5
    7 16807ρ
    7
    + 12005ρ
    6
    + 2058ρ
    5
    − 7350ρ
    4
    − 10920ρ
    3
    − 8280ρ
    2
    − 3600ρ − 720
    8 16384ρ
    8
    + 12288ρ
    7
    + 3584ρ
    6
    − 5376ρ
    5
    − 10080ρ
    4
    − 9240ρ
    3
    − 5355ρ
    2
    − 1890ρ − 315
    † See Erlang Redux Resolved! (July 28, 2016)
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 15 / 24

    View Slide

  19. Example 1 (M/M/8)
    From Table of polynomials for m = 8

    8 =
    S
    1 −
    16 384
    P7(ρ)
    ρ8
    Compute normalized R/S using Mathematica and its ErlangC function
    �������� ����{� = �� ρ = ����}�


    � - ρ�


    � - ��� �� ���
    -���-���� ρ-���� ρ�-���� ρ�-�� ��� ρ�-���� ρ�+���� ρ�+�� ��� ρ�
     ρ�

    � +
    �������[�� � ρ]
    � (� - ρ)
    
    �������� {�������� �������� �������}
    Morphing model Rφ
    8 = 1.11125 underestimates RE
    8
    (as expected)
    Corrected morphing model Rφ
    8
    agrees with Erlang RE
    8 = 1.17849
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 16 / 24

    View Slide

  20. Visualization in Complex Plane
    Analytic continuation into C: ρ → z = ρx
    + iρy
    Rm
    (z) poles determined by location of Pm
    (z) zeros
    Morphing approximation (m = 8)
    1 all zeros lie on unit disk
    2 symmetric roots of unity
    Erlang solution (m = 8)
    1 7 zeros inside unit disk
    2 asymmetric left & right half-planes
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 17 / 24

    View Slide

  21. Convergence of Erlang Zeros
    m = 256
    -1 -0.5 0.5 1
    Re()
    -

    Im()
    Szego curve
    Figure 1: RE
    m
    zeros lie inside the unit disk but outside the Szeg˝
    o bound
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 18 / 24

    View Slide

  22. Connection between Erlang and Szeg˝
    o
    1 In C: ρ → x + iy ≡ z, morphing model poles in the normalized residence
    time Rm
    = Rm
    /S are determined by a simple polynomial of degree m

    m
    (z) =
    1
    1 − zm
    (8)
    viz., the roots of unity zm = 1. Derivation requires a geometric series
    1 + z + z2 + · · · + zm truncated at m servers (Theorem 1).
    2 From Theorem 2, complex Erlang residence time can be written as
    RE
    m
    (z) =
    Pm−1(z)
    Pm
    (z)
    (9)
    Eq. (9) is a rational function associated with the Poisson distribution
    whose (real) CDF is
    F(k, m) = e−a
    m
    k=0
    ak
    k!
    with a = mρ. The summation is a truncated exponential series.
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 19 / 24

    View Slide

  23. Successive truncations of an (otherwise infinite) exponential series at
    n = 1, 2, . . . , can be regarded as polynomials of degree n
    However, the zeros of those successive polynomials move outside the
    unit disk |z| ≤ 1 in Figure 1 as n → ∞
    Szeg˝
    o (1926) tamed this behavior by scaling the roots: zn
    → ζn
    = zn
    /n
    He proved ζn
    converge, as n → ∞, to a complex function:
    |ζn
    e1−ζn | = 1 (10)
    Eq. (10) is the Szeg˝
    o bound: teardrop shaped (red) curve in Fig. 1
    Solve (10) numerically using the Lambert W function
    Remarkably, the zeros of (9) are self-scaling and therefore also
    converge on the Szeg˝
    o curve (no explicit rescaling required)
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 20 / 24

    View Slide

  24. Conclusion
    RE
    1
    has a purely algebraic derivation (no prob thy required)
    But RE
    m
    is inscrutible to derive in the same way
    Derived the morphing approximation Rφ
    m
    (Theorem 1)

    m
    can now be corrected to produce RE
    m
    (Theorem 2) 3
    Missing terms arise from a complicated polynomial Pm
    (ρ)
    Not what I had in mind when I started (analytically or visually)
    Complex analysis is not familiar to most IT people either
    But it’s correct and reveals how Rφ
    m
    is related algebraically to RE
    m
    Convenient visualization by extending ρ from R → C
    A new perspective for future work?
    Fairly certain Erlang didn’t look at M/M/m this way
    3
    TB-22 audience suggests publication in top-tier journal (along with proofs)
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 21 / 24

    View Slide

  25. Addendum (June 08, 2018)
    The morphing approximation (7)
    Rmorph
    =
    S
    1 − ρm
    (11)
    to the exact Erlang solution (5)
    Rexact
    = S + C(m, ρ)
    S
    m(1 − ρ)
    (12)
    has a relative error defined by
    ∆Rm,ρ
    =
    Rexact
    − Rmorph
    Rexact
    (13)
    A member of the IFORS audience asked if there was an upper bound on ∆R
    We can now answer in the affirmative
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 22 / 24

    View Slide

  26. The upper bound on ∆R (dashed curve) is given by
    ∆Rm,ρ
    <
    ln m1/4
    1 + ln m
    (14)
    10 1000 105 107
    m
    0.05
    0.10
    0.15
    0.20
    0.25
    ΔRm,ρ
    © 2018 Performance Dynamics
    ρ = 0.900000
    ρ = 0.990000
    ρ = 0.999000
    ρ = 0.999900
    ρ = 0.999990
    ρ = 0.999999
    Error bound
    All maxima in ∆R lie below (14) which asymptotically approaches 25%
    Prior to this result it could have been that ∆R reached 100%
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 23 / 24

    View Slide

  27. Questions?
    Performance Dynamics Company
    Castro Valley, California
    www.perfdynamics.com
    perfdynamics.blogspot.com
    facebook.com/PerformanceDynamics
    twitter.com/DrQz
    [email protected]
    +1-510-537-5758
    c 2018 Performance Dynamics Morphing M/M/m July 3, 2018 24 / 24

    View Slide