No CrossRef data available.
We use cookies to distinguish you from other users and to provide you with a better experience on our websites. Close this message to accept cookies or find out how to manage your cookie settings.
We use cookies to distinguish you from other users and to provide you with a better experience on our websites. Close this message to accept cookies or find out how to manage your cookie settings.
Published online by Cambridge University Press: 01 July 2004
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.
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.
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:
Ai,k is the time between the arrival of the kth type i customer and the previous arrival.
Bi,k is the service time of the kth type i customer.
Wi,k is the waiting time of the kth type i customer.
Si,k is the sojourn time of the kth type i customer (= Wi,k + Bi,k).
For each i, both {Ai,k}k≥1 and {Bi,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 Ai and Bi have the same distribution as Ai,k and Bi,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).
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 Wi,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 cS2 denote the mean and squared coefficient of variation of the sojourn time of the previous customer. If 1/k ≤ cS2 ≤ 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 cS2, provided the parameters p and μ are chosen as
If cS2 > 1, then the mean and squared coefficient of variation of the hyperexponential distribution
with density
matches E(S) and cS2, provided the parameters μ1, μ2, p1, and p2 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 cS2.
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 [Wi,0] = E [Wi,02] = 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 SN,k−1 if i = 1 and of Si−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 [Wi,k] and E [Wi,k2] according to (3), with FSN,k−1(·) and FSi−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 [Wi,k] and E [Wi,k2] as the approximation for E [Wi] and E [Wi2] , respectively, for i = 1,…,N; otherwise set k = k + 1 and i = 1 and go to step 2.
In this section, we consider the special case that the interarrival time Ai is exponentially distributed with mean 1/λi, i = 1,…,N. Let Wi and Si be the waiting time and the sojourn time in steady state of a type i customer, with LSTs Wi(s) and Si(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 (Ai − Si−1|Ai > Si−1) is again exponential with parameter λi, so
Note that P(Ai > Si−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 Si(s) = Wi(s)Bi(s) where Bi(s) is the LST of Bi, the equations for the transforms Wi(s) can be written in the form
From these equations we can solve WN(s), yielding
Of course, the other transforms Wi(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 WN(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 s0(= 0),s1,…,sN−1. We assume that these zeros are all distinct. Because the numerator of (7) must also vanish at s = s1,…,sN−1, we obtain the following equations:
Together with (10), this forms a set of N equations for N waiting probabilities Π1,…,Π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 Wi(s). This completes the determination of the transforms Wi(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 sk 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 [Wik] 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 [Wi], 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 Ai 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 Bi, we can analyze this case along the same lines as the exponential case.
Remark (Zeros): From (1) and (9) it follows that s0 = 0 is a simple zero of the denominator in (7), but it might happen that some of the other zeros coincide. For example, if s1 = s2, then s = s1 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 = s1.
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 × 106 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 103 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 × 106 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.
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 Ri 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 (Ri = R). The processing times of replenishment orders from stockpoint i are assumed to be i.i.d. random variables with mean E [Bi] and standard deviation σ(Bi). The first time the inventory position is inspected is denoted by ti0 (i = 1,…,N). The values of Ri, E [Bi], σ(Bi), and ti0 (i = 1,…,N) are listed in Table 6.
Input Parameters
Now, we are interested in the sojourn time or production lead time Si of replenishment orders of stockpoint i; obviously, the production lead time Si is required for setting an appropriate stock level.
The time Ai between the arrival of an item i replenishment order and item i − 1 replenishment order follows from ti0 and Ri (i = 1,…,N). Because Ai is deterministic, we are dealing with a cyclic D/G/1 queue. Given Ai and the first two moments of Bi, we can determine the first two moments of Wi and Si (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 cB1 = 0.37, cB2 = 0.32, cB3 = 0.20 and cB4 = 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.
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 Bi(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 maxi λi and radius δ + maxi λ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 C0. There are no other zeros in the right half plane since for all s with Re(s) ≥ 0 outside C0 we have
To complete the proof, we must remove the initial assumption that for some ε > 0, the transforms Bi(s) are analytic for all s with Re(s) > −ε. To this end, first consider, instead of Bi, the truncated service times min(Bi,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 maxi λi and radius maxi λi; this may be useful for the numerical calculation of the zeros.
Different Settings for a Model with Two Customer Types
Results for Different Settings of a Model with Two Customer Types
Parameter Settings
Average Percentage Errors of the Moment-Iteration Method
Maximum Percentage Errors of the Moment-Iteration Method
Schematic representation of the model.
Input Parameters
Numerical Results
You have
Access