Abstract
We consider two coupled queues with a generalized processor sharing service discipline. The second queue has a much smaller Poisson arrival rate than the first queue, while the customer service times are of comparable magnitude. The processor sharing server devotes most of its resources to the first queue, except when it is empty. The fraction of resources devoted to the second queue is small, of the same order as the ratio of the arrival rates. We assume that the primary queue is heavily loaded and that the secondary queue is critically loaded. If we let the small arrival rate to the secondary queue be , where , then in this asymptotic limit the number of customers in the first queue will be large, of order , while that in the second queue will be somewhat smaller, of order . We obtain a two-dimensional diffusion approximation for this model and explicitly solve for the joint steady state probability distribution of the numbers of customers in the two queues. This work complements that in (Morrison, 2010), which the second queue was assumed to be heavily or lightly loaded, leading to mean queue lengths that were or , respectively.
1. Introduction
The study of two coupled queues is a fundamental problem in queueing theory and applied probability. Classic examples include the shortest queue problem [1–3], the longer queue problem [4], fork-join models [5–7] and two coupled queues with generalized processor sharing [8–10], which is the subject of the present investigation. Computing the joint probability distribution for these models typically leads to functional equations that may sometimes be recast as boundary value problems [11], such as Dirichlet and Riemann-Hilbert problems.
Generalized processor sharing (GPS) models have become quite popular in recent years, as they provide scheduling algorithms that yield both service differentiation among different customer classes and also gains from statistical multiplexing. Some recent investigations and applications of such models appear in [12–15], where they are used, for example, for flow control in integrated service networks.
We consider here two parallel queues with respective Poisson arrival rates and , where . Thus, the arrivals to the second queue are much less frequent than those to the first, and we immediately scale the second arrival rate by , thus introducing . The service times are assumed to be exponentially distributed in both queues, with respective means and 1. Thus, we are taking the unit of time as the mean service time in the second queue. There is a single processor sharing server that works at unit rate and devotes of its capacity to the first queue and the remaining to the second queue, provided both queues are nonempty. If one queue is empty, the processor devotes all of its capacity to the other queue. The total load is given by , and we assume that the system is in heavy traffic so that this quantity will be close to 1. Hence, we define from the relation and assume that , so that the system is stable. This also means that the first queue is heavily loaded. The second queue has traffic intensity , and we say it is underloaded if , overloaded if and critically loaded if (more precisely ). The underloaded and overloaded cases were analyzed in [16]. We denote by (resp., ) the number of customers in the first (resp., second) queue and the joint steady state probability distribution will be denoted by .
For the underloaded case most of the mass occurs on the scale and , so there will tend to be only a few customers in the second queue. Asymptotically, has a product form behavior, with an exponential distribution in and a geometric distribution in (see [16]). This analysis was recently extended to an arbitrary number of parallel queues by Morrison and Borst [17], as long as one queue is heavily loaded and all the others are underloaded (with similar assumptions about arrival rates and processor-sharing factors as above). For the overloaded case most of the probability mass occurs for both , and in [16] a diffusion limit of the form is obtained. Here may be characterized as the solution to a parabolic PDE, in the variables and . Here, we will analyze the critically loaded case, which also leads to a diffusion limit with now , where will satisfy an elliptic PDE. Thus the critically loaded case has and leads to a somewhat more difficult problem than either the underloaded or overloaded cases.
We will obtain the PDE for as a limiting case of the difference equation(s) satisfied by , and explicitly solve the PDE by transform methods. We will obtain detailed results for the marginal distributions and , as well as the mean queue lengths. We will also obtain other approximations to that are valid on scales where and/or . In particular, we shall show that is on the scale and .
Previous work on this model includes Fayolle and Iasnogorodski [8] (see also [10]) and the more recent study of Guillemin and Pinchon [9]. There the authors consider the double generating function and obtain a functional equation for the boundary values and . This is ultimately converted to a Dirichlet problem, which is solved to yield the boundary values of in terms of elliptic integrals. One of the authors (J. A. Morrison) has verified that by analyzing the results in [8] and [9], in the asymptotic limit we consider and with the scaling and , and then inverting the double transform, we also obtain our main approximation . However, this involves a lengthy calculation that takes far more work than our direct approach, which consists of deriving a limiting PDE and solving it, along with appropriate boundary conditions. Also, this direct approach should work for other models of this type, including ones with finite capacities of customers, and with 3 or more coupled queues.
Other recent work on diffusion approximations for generalized processor sharing models includes Ramanan and Reiman [18] (see also references therein). This work, however, is more concerned with theoretical aspects of the diffusion approximations, such as convergence of the discrete problem to a certain diffusion process. Here, our focus is on obtaining the explicit solution to the limiting diffusion equation that arises from the balance equations. It is highly likely that this equation can be interpreted as the Kolmogorov forward equation for some appropriate diffusion process, but we do not consider such “process level” aspects here, as our approach is largely analytical.
Our approach has the merit that it can be used to compute correction terms to the diffusion approximations, and we do this in some cases here. Also, we treat scales other than the basic diffusion scale, where, for example, some of the variables remain discrete, which we then relate to the diffusion scale by asymptotic matching. This type of analysis is needed, for example, to accurately compute boundary probabilities.
From a mathematical viewpoint, the diffusion approximation we obtain (i.e., ) is somewhat nonstandard in that the density vanishes as and the approximation breaks down for . Also, the corner behavior of the problem is much different than what is typical. In [19] we analyzed a more general version of this model in another heavy traffic limit, assuming that the arrival rates and processor-sharing factors were of comparable magnitude. That analysis led to an elliptic PDE that was more complicated than the one obtained here, but probably more representative of typical diffusion approximations to two coupled queues, such as those considered in [20–22].
The present scaling limit leads to a separable, elliptic PDE in the variables and . Since the boundary conditions are somewhat simpler than those for the diffusion model in [19], we are able to obtain a more explicit solution to this equation, using classical transform theory [23]. We then evaluate the solution in various limiting cases, to obtain even simpler results that yield more insight into model behavior.
Yet another analysis of the model considered here is done in [24], but there it was assumed that both the arrival and service rates of the secondary customers are small, while the server devotes comparable resources to each queue.
The paper is organized as follows. In Section 2, we state the problem more precisely and give the balance equations satisfied by . In Section 3, we summarize all of our main results. The derivations are given in Section 4 for the scale and in Section 5 for the other ranges of .
Throughout the paper we will use the notation to mean , to mean , and to mean that is bounded for sufficiently close to .
2. Formulation
We consider two parallel infinite capacity queues for different traffic classes. The jobs arrive as Poisson processes with rates for the primary class, and for the secondary class, where . Hence, the secondary jobs arrive much less frequently than the primary ones. Moreover, it is assumed that the primary and secondary jobs have exponentially distributed service requirements with mean service times and 1, respectively, where . The server works at unit rate, and if neither queue is empty devotes fractions and of its effort to the primary and secondary queues, respectively, where . The corresponding service rates are and . If one queue is empty the server works at unit rate on the other queue, so that this model is work conserving. It is assumed that the primary queue is heavily loaded, with Moreover, we assume that the secondary queue is critically loaded, with where may have either sign. Thus, the asymptotic limit we consider has with and fixed, and then and vary with so that (2.1) and (2.2) hold.
Since the system is stable. Let denote the stationary probability that there are jobs in the primary queue and jobs in the secondary queue. Then, the balance equations satisfied by are The normalization condition is
The mean number of jobs in the primary and secondary queues are From Little’s Law, the corresponding mean waiting times are and . From [8], the conservation law of Kleinrock implies, using (2.1) and (2.2), that
3. Summary of Results
We consider and , with We then have
Proposition 3.1.
This gives the limiting density of , which applies for fixed values of and . We next evaluate this density in various limiting cases, such as and , to gain more insight into its structure, and to verify consistency with results in [16].
Corollary 3.2. If , and , then
Corollary 3.3. If and , then where the complementary error function is given by [25]
Corollary 3.4. If and , then
Corollary 3.5. If and , then
Remark 3.6. This matches with [16, Result 7], with and .
Corollary 3.7. If and , then
Remark 3.8. This matches with [16, Result 4], with and .
In Propositions 3.9 and 3.16 below, we give the limiting marginals of and , which apply for and fixed, respectively. Then, we simplify these marginal densities in various limiting cases.
Proposition 3.9. The scaled lowest order asymptotic approximation to the stationary distribution of the number of jobs in the secondary queue is where
Corollary 3.10. If and , then which is consistent with (3.3).
Corollary 3.11. If and , then
Corollary 3.12. If , then which is consistent with (3.6).
Corollary 3.13. If and , then which is consistent with (3.7).
Corollary 3.14. If and , then which is consistent with (3.8).
Remark 3.15. This matches with [16, Result 3], with and .
Proposition 3.16. The scaled asymptotic approximation to the stationary distribution of the number of jobs in the primary queue is where
Corollary 3.17. If , then
Remark 3.18. This matches with [16, Result 7] for , and with [16, Result 2] for since, for ,
We next give expansions for the mean queue lengths, for (with fixed).
Proposition 3.19. The lowest order asymptotic approximation to the stationary mean secondary queue length is where
Remark 3.20. We have verified from (3.17) that so that from (3.16), which is consistent with (2.9).
Corollary 3.21. If , then These match with [16, Result 1 and Corollary 10], respectively.
We next give some asymptotic results for that apply for fixed and and/or . We also give the “corner” behavior as .
Proposition 3.22.
(i) with and ,
(ii) with for , or for ,
(iii),
(iv),
(v) with ,
(vi) with ,
Remark 3.23. The results show that for the density is asymptotically of product form in a sector, and distinctly nonproduct form in the complimentary sector. For the product form behavior is absent. Item (v) shows that has an integrable singularity near the corner, while item (vi) shows that as . The second integral in (3.32) may be expressed in terms of a modified Bessel function, using the identity
The approximation is only valid for and . For other ranges of , other expansions must be constructed, and these we summarize below.
Proposition 3.24. , ,
,
,
Remark 3.25. We comment that for and so that the diffusion approximation still applies. For and item (i) still applies, and then , which is independent of and can be used to estimate the boundary probabilities , which are for . Note that in item (ii) and for in item (iii), is , which is larger than the order of magnitude of the diffusion approximation on the scale. But, the total mass in the range in (ii) is , while that in ranges (i) and (iii) is .
4. Analysis of the Main Diffusion Approximation
If we let and use (3.1), (2.1), and (2.2), in (2.3) and (2.4), we obtain to lowest order and the boundary condition We will discuss the second boundary condition, along , after (4.16).
We let It follows that
To solve (4.4) and (4.5) we apply a transform in the variable. Using the theory of distributions and Green’s functions for ordinary differential equations (see [23, p. 294, exercise 4.24]), we have the following transform pair: Here, the constant appears only when and the term corresponds to a single discrete eigenvalue in the spectral theory. By multiplying (4.7) by and integrating from and , we find that . Applying (4.6) with and , we let then integration by parts and the boundary condition (4.5) leads to Hence, from (4.4), so that Applying the inversion formula (4.7) Let Then, from (4.1) and (4.2), But, from the normalization condition (2.7), (3.1) and the Euler-Maclaurin summation formula [26], , so that Hence, from (4.3) and (4.13), Also, if we use (3.1), (2.1), and (2.2), in (2.5), we obtain the lowest order boundary condition , for . We conclude that and are proportional to a delta function at and hence, from (4.8) and (4.16), that . Proposition 3.1 follows from (4.3), (4.11), (4.12), and (4.15).
We may rewrite the integral in (3.2) as We then deform the contour of integration to one around a cut in the -plane from to , and let For there is a contribution from the pole at , and we obtain For and , the main contribution to the integral comes from , and Corollary 3.2 follows.
If and , then the integral in (4.19) is approximated by which leads [27] to Corollary 3.3. Corollaries 3.4 and 3.5 follow from the asymptotic approximation [25] , and the limiting value . If and , then from (4.19), we obtain the approximation and hence Corollary 3.7.
Proposition 3.9 follows from Proposition 3.1, (3.1) and the Euler-Maclaurin summation formula [26]. Next, from (4.19) For and , the main contribution to the integral comes from , and Corollary 3.10 follows. If and , then which leads [27] to Corollary 3.11. Corollaries 3.12, 3.13, and 3.14 follow directly from Corollaries 3.4, 3.5, and 3.7, respectively.
We now consider the first order correction term in (3.1). If we use (3.1), (2.1), and (2.2), in (2.3) and (2.4), we obtain for and , and the boundary condition We let It follows, from (4.15) and (4.24)–(4.26), that Hence, from (4.1) at , where
From Proposition 3.1, It follows that From (3.1), (4.13), (4.26), (4.28), and the Euler-Maclaurin summation formula [26], Hence, Proposition 3.16 follows in an elementary manner from (4.31) and (4.33).
If we let in (4.22), we obtain where is given by (3.21). Hence, The evaluation of the integral in (4.35) is routine, but tedious, and Proposition 3.19 follows from (3.1), since
We next establish the various asymptotic formulas in Proposition 3.22. We first note that the integrand in Proposition 3.1 is an even function of , and if is viewed as complex it has a simple pole at and branch points at Then, we can represent, for any , as the contour integral Here, is a horizontal contour in the -plane, on which The condition in (4.39) insures that if the pole at lies below the contour .
If we can shift the contour to the real -axis, and then (4.38) becomes the same as (3.2). If the pole must be taken into account in making this shift, and the residue from this pole yields the exponential terms in the left hand side of (3.2).
To evaluate (4.38) for we employ the saddle point method. There is a saddle point where so that The saddle is on the imaginary axis and the directions of steepest descent are . Then shifting the contour into another horizontal contour through leads to (3.24). Such a shift is always permissible if , but if we need the saddle to lie above the pole, that is, , and this occurs precisely when . We thus obtain the condition in item (i) of Proposition 3.22. If and the pole dominates the saddle point contribution, and we obtain (3.23).
For and a different analysis is needed, as now so the saddle approaches the real axis, where the integrand in (4.38) has as simple zero. In this case (which applies only if ), we shift back to the real axis and expand the integrand for . Using we thus obtain where . Evaluating the integral(s) in (4.43) leads to (3.25). If we consider the opposite limit, where and , then the saddle point approaches the upper branch point at . But by deforming to an integral about the branch cut we can show that the final result coincides with the expansion of (3.25) for , which is given by (3.28).
Next, we consider the corner behavior of as . Now, the main contribution to the integral will come from where is large. From (3.2) for , we then obtain which yields (3.30), and this can be shown to remain valid for .
Finally, we fix and let . Simply, setting in (3.2) leads to a divergent integral. However, an integration by parts leads to, for , where and Expanding (4.46) for and noting that, by contour integration (if ), we write the integrand as By using (4.47), identifying the terms in (4.48), and explicitly performing the differentiation with respect to , we ultimately obtain (3.31) and (3.32). This completes the (sketched) derivation of Proposition 3.22.
5. Analysis of Boundary and Corner Regions
We analyze cases where and/or . While these carry mass that is asymptotically small, they must be considered to insure that is properly normalized to higher orders in , and to compute higher order approximations to the moments. Also, to determine we used the boundary condition , and analysis of cases where and are small will allow us to examine this condition more carefully.
First we observe from (3.31) that vanishes linearly as , which indicates a nonuniformity in the asymptotics. We first consider the scale with and set Here, is a constant that will be determined by asymptotic matching. From (2.3), in terms of the variables and , we have Thus, the leading term must satisfy, since , and we write as which is a linear function of . The expansion for must satisfy the boundary condition (2.5), which is To leading order, (5.4) implies that so that (5.3) becomes . But, as cannot match to , since vanishes as . This indicates that we must analyze another scale, which has and (with fixed).
We thus set and By rewriting (2.3) on the scale, we obtain The limiting form of (5.6) as is , so we write On the scale the boundary condition (2.5), using also and , becomes To leading order as , (5.8) implies that Combining (5.7) and (5.9), we thus write the approximation on the scale as We then determine and by asymptotically matching (5.10), as , to as . Noting that , we find, in view of (3.31), that This establishes (3.34) in Proposition 3.24. Expression (5.10) remains valid for as then and we obtain , which is independent of and consistent with our previous analysis for .
We consider some limiting cases of the scale result in (3.34). As from (3.28) we obtain and then (3.32) yields, after some computations, It follows from (3.34) that The above results are consistent with those of Morrison in [16] for the case , where it was shown that Hence, for our approximation agrees with (5.14).
We will next consider the scale and matching this to the (or ) scale(s) will require the behavior of (3.34) as . For the second integral in the definition of in (3.32) dominates, and we obtain so that Also, will be less singular than , and, since , for the approximation becomes of order on the scale.
On the scale we will show that is asymptotically larger than for , which we just inferred to be . We thus set On the scale, the limiting form of (2.3), with (5.17), becomes while the boundary condition (2.5) yields Hence, is linear in , and writing we then obtain from (2.5), which we sum to get Next, we consider the boundary condition (2.4) and the corner condition (2.6). Since , we find that and also satisfies (5.22). Hence, we have and . The corner condition (2.6), with (5.17), yields to leading order and at the next order It follows from (5.24) that . Also, (2.5) with becomes so that (5.19) holds to leading order even if . To summarize, we have obtained It remains to determine and the constant . By letting and and matching to the scale, we conclude that as While these relations are consistent with (5.21), they do not determine for .
We thus examine another scale, that has and . By matching to (5.26) as , we expect that will be on this scale, so we set Then, from (2.3), we obtain to leading order while the boundary condition (2.4) leads to
Thus (5.29) corresponds to driftless diffusion in the variable and a driftless random walk in the discrete variable. We may replace (5.30) by the “artificial” boundary condition , by extending (5.29) to hold also at .
To analyze (5.29) consider the discrete Green’s function problem where is a “spectral” parameter. This is easily solved to yield By integrating (5.32) over a loop in the complex -plane that encircles the branch cut , we obtain where we evaluated the contour integral explicitly. Here, is the Kronecker delta symbol. From (5.33), we infer that for any sequence which leads to the transform pair Returning to (5.29) we set and applying the transform in (5.35) to (5.29) and (5.30) leads to the ODE Using and requiring also to decay as , we thus have and hence We determine the function by asymptotic matching.
First we consider the matching of as to as . For the major contribution to the integral in (5.39) will come from near , and the Laplace method yields Setting in (3.2) and letting , or, equivalently, using (3.30), we conclude that . Now, we match the and scales.
For as should agree with as , which is possible only if . Thus, in (5.39) must be orthogonal to for all and hence . The constant must equal , since .
Having determined we proceed to let in (5.39). If so that in (5.26). If Thus, by matching as to as we conclude that Then, we easily obtain Now, in (5.20) are determined, as is in (5.26), and we have established (3.36) in Proposition 3.24. Also, using in (5.39) leads to (3.35) in Proposition 3.24.
Thus, we have used a singular perturbation analysis to approximate on scales other than the scale. While these carry mass that is , the size of may actually be larger than , as is evident from (3.35).
Acknowledgment
The work of C. Knessl. was partly supported by NSF Grant no. DMS 05-03745 and by NSA Grant no. H 98230-08-1-0102.