WAITING TIME CHARACTERISTICS IN CYCLIC QUEUES

Sanne R. Smits; Ivo Adan; Ton G. de Kok

doi:10.1017/S0269964804183022

WAITING TIME CHARACTERISTICS IN CYCLIC QUEUES

Published online by Cambridge University Press: 01 July 2004

Sanne R. Smits ,

Ivo Adan and

Ton G. de Kok

Show author details

Sanne R. Smits: Affiliation:
Department of Technology Management, Eindhoven University of Technology, 5600 MB, Eindhoven, The Netherlands
Ivo Adan: Affiliation:
Department of Mathematics and Computer Science, Eindhoven University of Technology, 5600 MB, Eindhoven, The Netherlands, E-mail: iadan@win.tue.nl
Ton G. de Kok: Affiliation:
Department of Technology Management, Eindhoven University of Technology, 5600 MB, Eindhoven, The Netherlands

Article contents

Abstract
1. INTRODUCTION
2. MODEL DESCRIPTION
3. MOMENT-ITERATION METHOD
4. EXPONENTIAL INTERARRIVAL TIMES
5. NUMERICAL RESULTS
6. APPLICATION
APPENDIX
References

Rights & Permissions

Abstract

In this article, we study a single-server queue with FIFO service and cyclic interarrival and service times. An efficient approximative algorithm is developed for the first two moments of the waiting time. Numerical results are included to demonstrate that the algorithm yields accurate results. For the special case of exponential interarrival times, we present a simple exact analysis.

Information

Type: Research Article
Information: Probability in the Engineering and Informational Sciences , Volume 18 , Issue 3 , July 2004 , pp. 299 - 313

DOI: https://doi.org/10.1017/S0269964804183022 [Opens in a new window]
Copyright: © 2004 Cambridge University Press

1. INTRODUCTION

The present study concerns a multiclass queuing model with cyclic interarrival and service times. This model can be used, for example, when the inflow of customers depends on the day of the week or on the hour of the day. More specifically, this model arises in the modeling of a manufacturing system producing replenishment orders for stock locations that are controlled by periodic order-up-to policies. In that situation, the interarrival times are typically deterministic, and the order sizes are location dependent.

Cyclic queuing models have been studied by Morrice and Gajulapalli [3]. They considered a model with cyclic exponential interarrival and service times and derived bounds and exact results for the mean number in the system. Cohen [1] also studied the cyclic model and presented functional equations for the stationary waiting time distributions. These equations formulate a Hilbert boundary value problem, which can be solved if the Laplace–Stieltjes transforms (LSTs) of all the interarrival time distributions or all of the service time distributions are rational.

The central equation in this article is Lindley's equation for the waiting times. It is used to derive an iterative method to compute approximations for the first two moments of the waiting time. This method extends the one developed by de Kok [2]; it is applicable under generally distributed interarrival and service times and produces accurate results. Lindley's equation is also used to exactly determine the moments of the waiting time in the case of exponential (or Erlang) interarrival times.

The article is organized as follows. Section 2 provides a description of the model and introduces some notation. In Section 3, we present the moment-iteration method. In Section 4, we treat the special case of exponential (or Erlang) interarrival times and present an exact method for the computation of the moments of the waiting time. Numerical results are presented in Section 5, where we compare the approximations produced by the moment-iteration method with the exact results for the model with Erlang interarrival and service times and with simulation results for uniform and discrete interarrival time distributions. Finally, Section 6 is devoted to an application of the moment-iteration method.

2. MODEL DESCRIPTION

We consider a single-server queue providing FIFO service to N types of customer, numbered 1,…,N. The customers arrive in a cyclic pattern: first a type 1 customer, then one of type 2, then type 3 until type N, and then the cycle repeats. Define, for 1 ≤ i ≤ N and k ≥ 1, the following:

A_i,k is the time between the arrival of the kth type i customer and the previous arrival.

B_i,k is the service time of the kth type i customer.

W_i,k is the waiting time of the kth type i customer.

S_i,k is the sojourn time of the kth type i customer (= W_i,k + B_i,k).

For each i, both {A_i,k}_k≥1 and {B_i,k}_k≥1 are sequences of independent identically distributed (i.i.d.) random variables. The two sequences are also independent of each other and independent of the sequences for different i's. For stability, we assume that the traffic intensity ρ is less than one:

where the generic random variables A_i and B_i have the same distribution as A_i,k and B_i,k, respectively. It is readily verified that for k ≥ 1,

where (x)⁺ = max(x,0). These equations are the starting point for the iterative method presented in Section 3. This method approximates the first two moments of the stationary waiting times:

the limits exist by virtue of stability condition (1).

3. MOMENT-ITERATION METHOD

Equation (2) relates the waiting time of a customer to the sojourn time of the previous customer. From this equation we get the following expression for the nth moments of W_i,k:

Here, we concentrate on the first two moments (so n = 1,2). If the first two moments of the sojourn time of the previous customer are known and we fit a tractable distribution to these two moments, then the above expressions with the fitted distribution can be used to compute an approximation for the first two moments of the waiting (and sojourn) time of the present customer. For the two-moment fit, we may use a mixed Erlang or hyperexponential distribution, depending on whether the squared coefficient of variation is less or greater than one (see, e.g., Tijms [4]). More specifically, let E(S) and c_S² denote the mean and squared coefficient of variation of the sojourn time of the previous customer. If 1/k ≤ c_S² ≤ 1/(k − 1) for some k = 2,3,…, then the mean and squared coefficient of variation of the mixed Erlang distribution

with density

matches E(S) and c_S², provided the parameters p and μ are chosen as

If c_S² > 1, then the mean and squared coefficient of variation of the hyperexponential distribution

with density

matches E(S) and c_S², provided the parameters μ₁, μ₂, p₁, and p₂ are chosen as

This distribution is called the hyperexponential distribution with the gamma normalization because it has the property that its third moment as well matches that of the gamma distribution with mean E(S) and squared coefficient of variation c_S².

The above procedure is then repeated for the next customer and so on. The resulting iteration scheme is as follows.

Iteration scheme

1. Initially, set E [W_i,0] = E [W_i,0²] = 0 for i = 1,…,N and set i = k = 1.

2. Fit a tractable distribution to the first two moments of the sojourn time of the previous customer: Compute the first two moments and squared coefficient of variation of S_N,k−1 if i = 1 and of S_i−1,k if i > 1. Then, according to (4) and (5), the fitted distribution

is a mixture of two Erlang distributions with the same scale parameter if the squared coefficient of variation is less than one; otherwise it is a hyperexponential distribution with the gamma normalization.

3. Compute E [W_i,k] and E [W_i,k²] according to (3), with F_{S_N,k−1}(·) and F_{S_i−1,k}(·) replaced by the fitted distributions

, respectively.

4. If i < N, then set i = i + 1 and go to step 2; if i = N, compute the two sums

. If both are sufficiently small, then stop and use E [W_i,k] and E [W_i,k²] as the approximation for E [W_i] and E [W_i²] , respectively, for i = 1,…,N; otherwise set k = k + 1 and i = 1 and go to step 2.

4. EXPONENTIAL INTERARRIVAL TIMES

In this section, we consider the special case that the interarrival time A_i is exponentially distributed with mean 1/λ_i, i = 1,…,N. Let W_i and S_i be the waiting time and the sojourn time in steady state of a type i customer, with LSTs W_i(s) and S_i(s), respectively. Letting k → ∞ in (2), it follows that

where, by convention, a type 0 customer is the same as a type N customer. From these equations we get

for the transforms. By the memoryless property, the overshoot (A_i − S_i−1|A_i > S_i−1) is again exponential with parameter λ_i, so

Note that P(A_i > S_i−1) is the probability that a type i customer does not have to wait. Let us write

(Π_i is the probability of waiting). Further, using S_i(s) = W_i(s)B_i(s) where B_i(s) is the LST of B_i, the equations for the transforms W_i(s) can be written in the form

From these equations we can solve W_N(s), yielding

Of course, the other transforms W_i(s) are given by similar (symmetrical) expressions. To determine the unknown probabilities Π_i, we proceed as follows. First, we have to satisfy

For the denominator in (7), it holds that

and for the numerator,

Hence, from (8) we get

Further, since W_N(s) is well defined for Re(s) ≥ 0, it follows that whenever the denominator in (7) vanishes for some s with Re(s) ≥ 0, the numerator should also vanish. In the Appendix, we will prove that the denominator has exactly N zeros s with Re(s) ≥ 0, say s₀(= 0),s₁,…,s_N−1. We assume that these zeros are all distinct. Because the numerator of (7) must also vanish at s = s₁,…,s_N−1, we obtain the following equations:

Together with (10), this forms a set of N equations for N waiting probabilities Π₁,…,Π_N; it has a unique solution, because under the condition of stability (1), there is a unique stationary waiting time distribution and thus also a unique solution W_i(s). This completes the determination of the transforms W_i(s), as given by (7).

In the remainder of this section, we show how the moments of the waiting times can be determined. As a starting point, we take (6). Substituting the Taylor series

we obtain

Equating the coefficients of s^k on the left- and right-hand sides yields

and for k > 1,

Adding (12) over all i = 1,…,N gives the equation

This can be rewritten as (replace k by k + 1)

which is valid for k ≥ 1. Using (11)–(13), all moments E [W_i^k] can now be computed recursively. First, note that the addition of (11) over all i gives an identity; therefore, we can omit one of these equations. Then, together with (13) for k = 1, we have a set of N equations, from which the first moments E [W_i], i = 1,…,N, can be computed. To find the second moments, we use (12) and (13) for k = 2, where we can again omit one equation in (12). Also, note that the lower (first) moments occurring in these equations are now known. The third moments follow from (12) and (13) for k = 3, and so on. It will be clear that we can repeat using (12) and (13) to successively compute all moments (if they exist).

Remark 4.1 (Erlang Interarrival Times): If the interarrival time A_i is Erlang-r(λ_i) distributed, then we can think of the interarrival time as consisting of r subinterarrivals. By associating with each subarrival an arrival of a customer, where the first r − 1 customers have zero service times and the last, the rth customer, has a service time B_i, we can analyze this case along the same lines as the exponential case.

Remark (Zeros): From (1) and (9) it follows that s₀ = 0 is a simple zero of the denominator in (7), but it might happen that some of the other zeros coincide. For example, if s₁ = s₂, then s = s₁ is also a double zero of the numerator in (7); thus, an additional equation for the probabilities Π_i can be obtained by requiring that the derivative of the numerator also vanishes at s = s₁.

5. NUMERICAL RESULTS

In this section, we validate the approximative moment-iteration method by simulation (and exact results, if available). Table 1 shows 30 different settings for the mean and squared coefficient of variation of the interarrival and service times for a model with 2 customer types. In each setting, the distributions of the interarrival and service times are matched to the mean and squared coefficient of variation according to the recipe described in Section 3. The results are listed in Table 2; the simulation results are based on 10 independent replicas of 6 × 10⁶ arrivals.

Different Settings for a Model with Two Customer Types

Results for Different Settings of a Model with Two Customer Types

Note that the waiting time characteristics of the two customer types can differ substantially when their interarrival time and service time characteristics are different. Table 2 indicates that the approximation method produces accurate estimates for the first two moments of the waiting times. This is confirmed by a much more extensive investigation, the results of which are now summarized.

We consider models with deterministic, uniform, exponential or (mixed) Erlang, and hyperexponential distributions for the interarrival and service times. The models are indicated with the well-known three-letter code A/B/1; it means that all customer types have the same type of interarrival time distribution, indicated by A, and the same type of service time distribution, indicated by B, but, of course, the parameters of the distributions depend on the customer type. For each case of interarrival and service time distributions, we randomly generate 10³ settings of the parameters, and each setting will be evaluated for 2, 5, and 25 customer types. If only mean and squared coefficients of variation are specified, we fit a distribution according to the recipe described in Section 3. Table 3 gives an overview of the parameter settings, where U(a,b) denotes the uniform distribution on (a,b), and U(S) denotes the (discrete) uniform distribution on the points in the set S.

Parameter Settings

Note that, according to Table 3, we will generate settings with different traffic intensities. We divide the settings into three categories: low load (0.4 < ρ < 0.6), medium load (0.6 ≤ ρ < 0.8), and high load (ρ ≥ 0.8). We compare the estimates produced by the moment-iteration method with the exact results of Section 4 in case of exponential or Erlang interarrival times and with simulation results otherwise. Again, the simulation results are based on 10 independent replicas of 6 × 10⁶ arrivals. In Tables 4 and 5, we display percentage errors; δ^{E [W]} and max δ^{E [W]} denote, respectively, the mean and maximum percentage error of the average of the waiting time over all customer types and all settings generated in a category, and δ^σ(W) and max δ^σ(W) denote, respectively, the average and maximum percentage error of the standard deviation of the waiting time over all customer types and all settings generated in a category.

Average Percentage Errors of the Moment-Iteration Method

Maximum Percentage Errors of the Moment-Iteration Method

We can conclude that the moment-iteration algorithm produces accurate results, especially for high loads; typically, the mean of the waiting times is more accurately estimated than the standard deviation. For low loads, the percentage error of the moments of the waiting time might be high; the absolute error, however, will be modest (in comparison with the service times). Further, it seems that the type of distribution of the interarrival times does not influence the quality of the estimates produced by the moment-iteration method; the traffic intensity and the squared coefficients of variation of the service times are more crucial.

6. APPLICATION

In this section, we apply the moment-iteration method to the manufacturing problem mentioned in Section 1. We consider a manufacturing system with N stockpoints, one for each item, and one production facility, which produces all the items in a FIFO order. The objective is to determine the sojourn time of the replenishment orders placed by the stockpoints. Figure 1 is a schematic representation of the model for N = 4.

Schematic representation of the model.

The stockpoints are controlled by periodic order-up-to policies. The periodic order-up-to policy operates as follows. Every R_i units, the inventory position of stockpoint i is inspected and a replenishment order is placed at the production facility to raise the inventory position up to the order-up-to level. The inventory position is defined as the physical inventory level plus the stock on order minus the backorders. We assume that the review period of each item is identical (R_i = R). The processing times of replenishment orders from stockpoint i are assumed to be i.i.d. random variables with mean E [B_i] and standard deviation σ(B_i). The first time the inventory position is inspected is denoted by t_i⁰ (i = 1,…,N). The values of R_i, E [B_i], σ(B_i), and t_i⁰ (i = 1,…,N) are listed in Table 6.

Input Parameters

Now, we are interested in the sojourn time or production lead time S_i of replenishment orders of stockpoint i; obviously, the production lead time S_i is required for setting an appropriate stock level.

The time A_i between the arrival of an item i replenishment order and item i − 1 replenishment order follows from t_i⁰ and R_i (i = 1,…,N). Because A_i is deterministic, we are dealing with a cyclic D/G/1 queue. Given A_i and the first two moments of B_i, we can determine the first two moments of W_i and S_i (i = 1,…,N) by using the moment-iteration method presented in Section 3. The results are presented in Table 7.

Numerical Results

Note that the traffic intensity of the cyclic D/G/1 queue describing the manufacturing system is 0.91. Hence, Table 4 indicates that in this case (with N = 4), we can expect errors in the mean waiting time close to 1% and in its standard deviation close to 4%; this seems sufficiently accurate for practical purposes. From Table 6, it follows that c_B₁ = 0.37, c_B₂ = 0.32, c_B₃ = 0.20 and c_B₄ = 0.20. In Table 7, we see that replenishment orders from stockpoint 1 clearly benefit from the smaller variation in processing times of replenishment orders from stockpoints 3 and 4.

APPENDIX

In this appendix, we will prove that the denominator of (7) has N zeros s with Re(s) ≥ 0. To do so, we will use Rouche's theorem, which reads as follows. Let f (s) and g(s) be analytic functions inside and on a smooth contour C and suppose that |g(s)| < | f (s)| on C. Then, f (s) and f (s) + g(s) have the same number of zeros inside C. Of course, we may also replace f (s) + g(s) by g(s) − f (s) in this formulation; actually, that is the form we will use.

We first assume that there is some ε > 0 such that the transforms B_i(s),i = 1,…,N, are analytic for all s with Re(s) > −ε. This assumption holds, for example, for service time distributions with a finite support or an exponential tail. Now, take

So g(s) − f (s) is the denominator of (7). As a contour, we take the circle C_δ with center max_i λ_i and radius δ + max_i λ_i with 0 < δ < ε. It is easily verified that for all s on C_δ,

Note that

Hence, the stability condition (1) states that g′(0) > f′(0). Thus, for sufficiently small δ > 0, it holds that g(−δ) < f (−δ), which implies together with (A.1) that |g(s)| < | f (s)| for all s on C_δ. Rouche's theorem now guarantees that f (s) and g(s) − f (s) have the same number of zeros inside C_δ. Since f (s) has N zeros inside C_δ, the same holds for g(s) − f (s). Letting δ tend to zero, it follows that g(s) − f (s) has exactly N zeros inside or on the circle C₀. There are no other zeros in the right half plane since for all s with Re(s) ≥ 0 outside C₀ we have

To complete the proof, we must remove the initial assumption that for some ε > 0, the transforms B_i(s) are analytic for all s with Re(s) > −ε. To this end, first consider, instead of B_i, the truncated service times min(B_i,K), where K > 0 is some constant. For these truncated service times, the claim for the zeros of the denominator (7) holds; by letting K tend to infinity, the claim also follows for the original service time distributions.

Remark: In fact, we not only showed that the existence of N zeros s with Re(s) ≥ 0 but also that they are located inside or on the circle with center max_i λ_i and radius max_i λ_i; this may be useful for the numerical calculation of the zeros.

References

REFERENCES

Cohen, J.W. (1996). Periodic Pollaczek waiting time processes. In C.C. Heyde, Y.V. Prohorov, R. Pyke, & S.T. Rachev (eds.), Athens Conference on Applied Probability and Time Series Analysis: Applied Probability in Honor of J.M. Gani, pp. 361–377. New York: Springer-Verlag.

de Kok, A.G. (1989). A Moment-iteration method for approximating the waiting-time characteristics of the GI/G/1 queue. Probability in the Engineering and Informational Sciences 3: 273–287.Google Scholar

Morrice, D.J. & Gajulapalli, R.S. (1994). A single server queue with cyclically indexed arrival and services rates. Queueing Systems 15: 165–198.Google Scholar

Tijms, H.C. (1994). Stochastic models: An algorithmic approach. Chichester: Wiley.