De Cuypere, Eline

De Turck, Koen

Fiems, Dieter

English

Eline De Cuypere

Koen De Turck

Dieter Fiems

Crossref

Operations Research Letters

A Maclaurin-series expansion approach to multiple paired queues

Ghent University Academic Bibliography

A Maclaurin-series expansion approach to multiple paired queuesEline De Cuypere, Koen De Turck, Dieter Fiems∗Ghent University, Department of Telecommunications and Information ProcessingAbstractMotivated by kitting processes in assembly systems, we consider a Markovian queueing system with K paired finite-capacity buffers. Pairing means that departures from the buffers are synchronised and that service is interrupted if anyof the buffers is empty. To cope with the inherent state-space explosion problem, we propose an approximate numericalalgorithm which calculates the first L coefficients of the Maclaurin series expansion of the steady-state probability vectorin O(KLM) operations, M being the size of the state space.1. IntroductionWe consider a system of K queues, each queue havingfinite capacity. Let Ci denote the capacity of the ith queue.Moreover, for each of the queues, customers arrive in ac-cordance with an independent Poisson process, let λi > 0denote the arrival rate in queue i. Departures from thedifferent queues are paired which means that there are si-multaneous departures from all queues with rate µ as longas all queues are non-empty. If one of the queues is empty,there are no departures.The queueing system at hand is motivated by kittingprocesses in assembly systems. A kitting process collectsthe necessary parts for a given end product in a con-tainer prior to assembly. While conceptually simple, kit-ting comes with many advantages. Kitting clearly miti-gates storage space requirements at the assembly stationsince no part inventories need to be kept there. Moreover,parts are placed in proper positions in the container suchthat assembly time reductions can be realised [3, 9]. Akitting process is obviously related to a paired queueingsystem: the inventories of the different parts that go intothe kit correspond to the different buffers, the kitting timecorresponds to the service time and kitting is blocked ifone or more parts are missing.Paired queueing systems have been studied by variousauthors. Harrison [4] studies stability of paired queueingunder very general assumptions: K ≥ 2 infinite-capacitybuffers, generally distributed interarrival times at the dif-ferent buffers and generally distributed service times. Heshows that it is necessary to impose a restriction on thesize of the buffer to ensure stability of the queueing sys-tem. In particular, the distribution of the vector of wait-ing times (in the different queues) of the components of a∗Corresponding author (Dieter.Fiems@UGent.be), Ghent Univer-sity, Department of Telecommunications and Information ProcessingSt-Pietersnieuwstraat 41, 9000 Gent, Belgium.paired customer is shown to be defective. The inherent in-stability was also demonstrated in [8] where the excess —the difference between the queue sizes — is studied in thetwo-queue case. Assuming finite capacity buffers, Hoppand Simon developed a model for a two-buffer kitting pro-cess with exponentially distributed processing times forkits and Poisson arrivals [5]. The exponential service timesand Poisson arrival assumptions were later relaxed in [12]and [2], respectively. For paired queueing systems withmore than two finite buffers, the size of the state-spaceof the associated Markov chain grows quickly, even forthe case of Poisson arrivals and exponential service times.Hence, most authors focus on approximations; a recent ac-count on approximations of multi-buffer paired queueingsystems can be found in [10]. Also the present letter in-vestigates approximations for multi-buffer paired queueingsystems. In particular, we propose a numerical evaluationmethod for Markovian paired queueing systems which re-lies on a Maclaurin-series expansion of the steady-stateprobability vector. For an overview on the technique ofseries expansions in stochastic systems, which is knownunder the names light traffic analysis or stochastic pertur-bation, we refer the reader to the surveys in [1, 7]. Finally,we note that the paired queueing system somewhat re-sembles a fork-join queueing system; see e.g. [6] and thereferences therein. However, in fork-join queueing systemsboth arrivals and departures in the different buffers aresynchronised, which leads to entirely different dynamics.2. Maclaurin-series expansionAs arrivals in the different queues are modelled by Pois-son processes and the service time distribution is expo-nential, the state of the system is described by a vectori ∈ C whose kth element corresponds to the queue size ofthe kth buffer. Here C = C1 × . . . × CK denotes the statespace of this continuous-time Markov chain (CTMC), withCk = {0, 1, . . . , Ck} being the set of possible levels of queuePreprint submitted to Elsevier January 28, 2014k. Let pi(i) be the steady-state probability of state i, i ∈ C.These steady-state probabilities satisfy the following set ofbalance equations,pi(i1, i2, . . . , iK)(µK∏`=11{i`>0} +K∑`=11{i`<C`}λ`)=pi(i1 + 1, i2 + 1, . . . , iK + 1)µK∏`=11{i`<C`}+K∑`=1pi(i1, . . . , i`−1, i` − 1, i`+1, . . . , iK)λ`1{i`>0} , (1)for all i = (i1, i2, . . . , iK) ∈ C and where 1{x} is the indica-tor function which equals one if x is true and equals zerootherwise. While the former system of equations is easilysolved if there are only a few queues with low capacity,the size of the state space explodes for even a moderatenumber of queues and reasonable queue capacities and adirect solution is computationally infeasible.To mitigate this state space explosion problem, we relyon a Maclaurin series expansion in µ. It is shown in the ap-pendix that pi(i) is analytic in µ = 0 and therefore admitsthe representation,pi(i) =∞∑n=0pin(i)µn ,for 0 ≤ µ < µ0 and for i ∈ C. Here µ0 is a non-negativevalue for which a lower bound is provided in the appendix.Substituting the former expression in the balance equa-tions yields,∞∑n=0pin(i1, i2, . . . , iK)µn×(µK∏`=11{i`>0} +K∑`=11{i`<C`}λ`)=∞∑n=0pin(i1 + 1, i2 + 1, . . . , iK + 1)µn+1K∏`=11{i`<C`}+∞∑n=0K∑`=1pin(i1, . . . , i`−1, i`−1, i`+1, . . . , iK)λ`µn1{i`>0} .For i ∈ C∗ = C\{c}, with c .= [C1, C2, . . . , CK ], comparingthe terms in µ0 on both sides of the former equation yields,pi0(i1, i2, . . . , iK) = 0 ,whereas comparing the terms in µn for n > 0 gives,pin(i1, i2, . . . , iK) =1∑K`=1 1{i`<C`}λ`×(1{n>0}pin−1(i1 + 1, i2 + 1, . . . , iK + 1)K∏`=11{i`<C`}+K∑`=1pin(i1, . . . , i`−1, i` − 1, i`+1, . . . , iK)λ`1{i`>0}− 1{n>0}pin−1(i1, i2, . . . , iK)K∏`=11{i`>0}). (2)For i = c, such a comparison does not yield an expressionfor pin(i). To determine the remaining unknown, we invokethe normalisation condition:∑i∈Cpi0(i) = 1 ,∑i∈Cpin(i) = 0 .Solving for pin(c) then yields,pi0(c) = 1 , pin(c) = −∑i∈C∗pin(i)Once the series expansion of the steady state distribu-tion has been obtained, the expansions of various perfor-mance measures directly follow. Let X ∼ pi, then for aperformance measure J = E[f(X)] we have,J =∑i∈Cf(i)pi(i) =∑i∈Cf(i)∞∑n=0pin(i)µn=∞∑n=0∑i∈Cf(i)pin(i)µn =∞∑n=0Jnµn , (3)for 0 ≤ µ < µ0 with,Jn =∑i∈Cf(i)pin(i) .The interchange of the summations is justified by the finite-ness of C and the convergence of ∑n pin(i)µn for all i ∈C. As such, any term Jn in the expansion of a perfor-mance measure J can be calculated from the correspond-ing vector pin of the expansion of the steady-state vec-tor. Performance measures of interest include amongstothers the `th order moment of the queue content of thekth queue (f(i) = i`k), the blocking probability (f(i) = 1−∏Kj=1 1{ij>0}) and the throughput (f(i) = µ∏Kj=1 1{ij>0}).Computational complexity. From (2), calculation of pin(i)takes at most K + 2 additions and one division (assumingthe rate sums are known). Hence, the computational com-plexity of calculating pin is O(KM), with M = |C| the sizeof the state space. Having the same complexity for everyadditional term in the expansion, calculating the first Lcoefficients then has complexity O(KML).As the size of the state space is very large, limited mem-ory consumption is equally important. To limit memoryconsumption to the size of storing only one steady-statevector one can proceed as follows. Assuming one is mainlyinterested in the expansion of a number of performancemeasures, note that once the nth term of the expansion ofthe steady state vector is determined, the corresponding2terms in the expansions of various performance measurescan be determined as well; see (3). Hence, there is no needto keep track of previous terms of the expansion of steady-state probabilities unless they are required for further cal-culations of coefficients of steady state probabilities. From(2) one sees that pin(i) is expressed in terms of pin−1(j),with j larger then i (lexicographically). This means thatthe coefficients of the vector pin−1 can be overwritten pro-gressively during the calculation of pin and memory foronly one vector of size M is needed.3. Numerical resultsTo illustrate our series expansion approach, we nowassess its accuracy by means of some numerical examples.First, consider a system with K = 5 paired queues, eachqueue having capacity C = 10. Moreover, the arrival in-tensity at each queue is equal to λ = 1. Hence, the pairedqueueing system is symmetric and performance measuresare equal for all queues. Figures 1(a) and 1(b) depictthe mean queue content and the blocking probability ina queue versus the service rate µ, respectively. For bothfigures, series expansions of various orders are depicted asindicated (L = 1, 2, 5 for figure 1(a) and L = 10, 11, 12for figure 1(b)), as well as simulation results which allowfor assessing the accuracy of the series expansions. As ex-pected, the mean queue content decreases and the blockingprobability increases as the service rate µ increases. More-over, for µ = 0, the queues are completely filled as thereis no service. From figure 1(a), it is observed that the ap-proximation method at hand is accurate for low orders ofthe expansion (L = 5) whereas more terms are needed toaccurately determine the blocking probability (L = 12);see figure 1(b). As the computation time of the series ex-pansion is linear in the number of terms in the expansion,accurately assessing the blocking probability takes morethan twice the computation time of assessing the meanqueue content.Figure 2(a) depicts the mean of the queue content ofthe first and second queue out of 5 paired queues, whereasfigure 2(b) depicts the corresponding variances. For bothfigures, the expansion of order N = 20 is compared withsimulation results. The capacity equals 10 for all queues,and the arrival intensity in all but the first queue equalsλi = 1, i = 2, . . . , 5. The arrival rate in the first queue islowered to λ1 = 0.8. In comparison with the symmetricpaired queueing system of figure 1(a), the mean queue con-tent increases for the second queue. This does not come asa surprise. Decreasing the arrival rate in the first queue im-plies that this queue is empty more often, thereby blockingservice in the other queues. Finally, note that the varianceincreases for increasing µ, µ = 0 corresponds to the casethat the queue content deterministically equals the queuecapacity for all queues, hence the variance is zero.Appendix: Convergence of the power seriesWe now justify the series expansion. The basic ideasin this section date back to the seminal work of Schweitzer[11]. The series expansion is validated by explicitly con-structing such an expansion. We first introduce some ad-ditional notation and the basic notion of the deviation ma-trix of a CTMC.Let pi(µ) denote the steady state solution [pi(i)]i∈C ofthe balance equations (1). We have made the dependenceof pi(µ) on µ explicit for ease of notation. The balanceequations can then be written in matrix notation as fol-lows,pi(µ)Q(µ) = pi(µ)(Q0 + µQ1) = 0 , (4)where Q(µ) is the |C| × |C| generator matrix of the CTMCand where Q0 and Q1 are known matrices that do not de-pend on µ. In view of the system assumptions it is readilyseen that Q(0) = Q0 only has one recurrent state, i.e. c(the full state) is recurrent and all the others are tran-sient. Therefore, the stationary vector pi(0) exists, withstate pi(0)(c) = 1 and pi(0)(i) = 0 for i ∈ C∗.Let D0 be the deviation matrix of the CTMC withgenerator matrix Q0,D0 =∫ ∞0(P0(t)−Π0)dt . (5)Here the family {P0(t) = exp(Q0t), t ≥ 0} is the Markovsemigroup of the CTMC, and Π0 = limt→∞ P0(t) = 1′pi(0),1′ being a column vector of ones. As the state-space C isfinite, the deviation matrix is well defined. Moreover, thedeviation matrix satisfies D01′ = 0 — the row sums arezero — and,D0Q0 = Q0D0 = Π0 − I . (6)Theorem 1. The solution pi(µ) of the CTMC adheres tothe following power series expansion,pi(µ) =∞∑k=0(pi(0)(Q1D0)k)µk , (7)for 0 ≤ µ < µ0, µ−10 being the spectral radius of Q1D0.Moreover, µ0 is bounded from below by µ∗0 and µ∗1,µ∗0 =(2∫ ∞0(1−K∏k=1F (t;Ck, λk))dt)−1≥(2K∑k=1Ckλk)−1= µ∗1 ,with F being the Erlang distribution,F (t;Ck, λk) = 1−Ck−1∑n=01n!e−λkt(λkt)n .30.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9µ7.07.58.08.59.09.510.0Mean number of itemsN = 1N = 2N = 5simulation(a)0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7µ0.000.010.020.030.04Blocking probabilityN = 10N = 11N = 12simulation(b)Figure 1: Mean queue content and blocking probability for a symmetric paired queueing system.0.0 0.1 0.2 0.3 0.4 0.5 0.6µ7.07.58.08.59.09.510.0Mean number of itemsbuffer 1buffer 2simulation buffer 1simulation buffer 2(a)0.0 0.1 0.2 0.3 0.4 0.5 0.6µ0246810Variance of the number of itemsbuffer 1buffer 2simulation buffer 1simulation buffer 2(b)Figure 2: Mean (a) and variance (b) of the queue content of an asymmetric paired queueing system.4Proof. Multiplying (4) by D0 and invoking (6) yields,pi(µ)(Q0 + µQ1)D0 = pi(µ)(Π0 − I) + pi(µ)µQ1D0 = 0 .Moreover, we have pi(µ)Π0 = pi(µ)1′pi(0) = pi(0), such that,pi(µ)(I − µQ1D0) = pi(0) .The spectral radius of µQ1D0 is µ/µ0. Hence for µ < µ0,(I−µQ1D0) is invertible and the Neumann series convergesto the inverse,∞∑k=0(µQ1D0)k = (I − µQ1D0)−1 .Combining the previous expressions immediately yieldsthe series expansion (7).As all elements but the last column of Π0 are zero,only the last column of D0 may contain negative values;see (5). Moreover, the row sums of D0 are zero, hence thelast column is equal in absolute value to the sum of theother columns. The entries in the last column of D0 havethe following interpretation,[D0]ic = −∫ ∞0(1− [P0(t)]ic)dt = −E[Ti] ,where Ti is a random variable denoting the time it takes toreach the full state c from state i (assuming no departures).This interpretation shows that γ.= E[T0] ≥ E[Ti] for alli ∈ C where 0 denotes the empty state.The time to fill up the ith queue is Erlang distributedwith Ci stages and rate λi and the time to fill up all queuesis the maximum of K Erlang distributed random variables.Therefore, the cumulative distribution of T0 is the productofK Erlang distributions and γ is calculated by integratingthis distribution,γ =∫ ∞0(1−K∏k=1F (t;Ck, λk))dt .Moreover, the maximum of K non-negative random vari-ables is bounded from above by the sum of these randomvariables, which yields the following crude upper boundfor γ,γ ≤K∑k=1Ckλk, (8)the kth term in the sum on the right-hand side correspond-ing to the mean time to fill up the kth queue.As the row sums of Q1 are zero (Q(µ) is a generatormatrix for every µ), we have Q1Π0 = 0. Moreover, for anyinduced matrix norm, we have ‖Q1D0‖ ≥ µ0. Therefore,we find,µ−10 ≤ ‖Q1D0‖ = ‖Q1(D0 + γΠ0)‖ ≤ ‖Q1‖‖D0 + γΠ0‖ .Particularly using the maximum absolute row sum norm,we have ‖Q1‖ = 2; [Q1]ii = −1 if all queues are non-empty in state i and 0 if this not the case such that the2 3 4 5 6 7 8 9 10C0.00.10.20.30.40.50.60.70.80.9Convergence domainµ0µ ∗1µ ∗0Figure 3: spectral radius µ0 and lower bounds µ∗0 and µ∗1.corresponding row sums equal 2 and 0 respectively. In viewof the definition of γ, one easily verifies that the matrixD0 + γΠ0 has no negative entries. Recalling that D0 haszero row sums, this shows that all row sums of D0 + γΠ0equal γ: ‖D0 + γΠ0‖ = γ and,1µ0≤ 2γ .= 1µ∗0,which proves the lower bound µ∗0 for µ0. The lower boundµ∗1 follows from µ−10 ≤ 2γ and the crude bound (8) forγ.To illustrate Theorem 1, figure 3 depicts µ0, the spec-tral radius of Q1D0 and the lower bounds µ∗0 and µ∗1 fora system with K = 3 paired queues, each queue having avarying capacity from 2 to 10. As the figure shows, thebounds are much smaller than the convergence radius. Itshould be noted that both bounds are easy-to-derive butalso rather loose bounds on the convergence radius. Thebounds above can be made tighter by (1) not relying onthe submultiplicative property of the matrix norm; (2) amatrix norm which is more adapted to this model. Boththese approaches quickly lead to lengthy calculations andwe consider them to be outside of the scope of the paper.AcknowledgementsThe second author is a post-doctoral fellow with theResearch Foundation, Flanders. This research has beenfunded by the Interuniversity Attraction Poles Programmeinitiated by the Belgian Science Policy Office.References[1] B. B laszczyszyn, T. Rolski, V. Schmidt. Advances in Queue-ing: Theory, Methods and Open Problems, chapter Light-trafficapproximations in queues and related stochastic models. CRCPress, Boca Raton, Florida, 1995.5[2] E. De Cuypere, D. Fiems. Performance evaluation of a kittingprocess. In: Proceedings of the 18th International Conferenceon Analytical and Stochastic Modelling Techniques and Appli-cations (ASMTA 2011), pp. 175–188, Venice, June 2011.[3] B. Johansson, M. Johansson. High automated kitting systemfor small parts: a case study from the Volvo Uddevalla plant. In:Proceedings of the 23rd International Symposium on Automo-tive Technology and Automation, pp. 75–82, Vienna, Austria,1990.[4] J. Harrison. Assembly-like queues. Journal Of Applied Proba-bility 10:354–367, 1973.[5] W.J. Hopp, J.T. Simon. Bounds and heuristics for assembly-likequeues. Queueing Systems 4:137–156, 1989.[6] S.S. Ko, R.F. Serfozo. Response times in M/M/s fork-join net-works. Advances in Applied Probability 36(3):854–871, 2004.[7] I. Kovalenko. Rare events in queueing theory. A survey. Queue-ing systems 16(1):1–49, 1994.[8] G. Latouche. Queues with paired customers. Journal of AppliedProbability 18(3)684–696, 1981.[9] R. Ramakrishnan, A. Krishnamurthy. Analytical approxima-tions for kitting systems with multiple inputs. Asia-PacificJournal of Operations Research 25(2):187–216, 2008.[10] R. Ramakrishnan, A. Krishnamurthy. Performance evaluationof a synchronization station with multiple inputs and popula-tion constraints. Computers & Operations Research 39:560–570,2012.[11] P.J. Schweitzer. Perturbation Theory and Finite MarkovChains, Journal of Applied Probability 5(2):401–413, 1968.[12] M. Takahashi, H. Osawa, T. Fujisawa. On a synchronizationqueue with two finite buffers. Queueing Systems 36:107–23,2000.6

https://biblio.ugent.be/publication/4324157/file/8521901.pdf

A Maclaurin-series expansion approach to multiple paired queues

Abstract

Similar works

Full text

Available Versions

Crossref

Ghent University Academic Bibliography