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

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

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

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

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

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

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

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

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

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

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

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

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

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

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 Rφ 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

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

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

New Analytic Result Theorem 2 (Gunther 2016†) The missing terms correct the denominator of (7) as follows: Rφ 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

Example 1 (M/M/8) From Table of polynomials for m = 8 Rφ 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

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

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

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 Rφ 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

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

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) Rφ 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

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

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

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