Two Coupled Queues with Vastly Different Arrival Rates : Critical Loading Case

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 O ε , where 0 ≤ ε 1, then in this asymptotic limit the number of customers in the first queue will be large, of order O ε−1 , while that in the second queue will be somewhat smaller, of order O ε−1/2 . 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 O ε−1 or O 1 , respectively.


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 Advances in Operations Research different customer classes and also gains from statistical multiplexing.Some recent investigations and applications of such models appear in [12][13][14][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 0 < ε 1.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 1/μ 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 1 − εκ 1 − O ε 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 λ/μ εσ 1 − εω and assume that ω > 0, 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 σ/κ < 1, overloaded if σ/κ > 1 and critically loaded if σ/κ ≈ 1 more precisely σ/κ 1 O √ ε .The underloaded and overloaded cases were analyzed in 16 .We denote by N 1 resp., N 2 the number of customers in the first resp., second queue and the joint steady state probability distribution will be denoted by p m, n Prob N 1 m, N 2 n .For the underloaded case most of the mass occurs on the scale m O ε −1 and n O 1 , so there will tend to be only a few customers in the second queue.Asymptotically, p m, n has a product form behavior, with an exponential distribution in εN 1 and a geometric distribution in N 2 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 m, n O ε −1 , and in 16 a diffusion limit of the form p m, n ∼ ε 2 P * ζ, τ ε 2 P * εm, εn is obtained.Here P * 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 p m, n ∼ ε 3/2 φ 0 ζ, w ε 3/2 φ 0 εm, √ εn , where φ 0 will satisfy an elliptic PDE.Thus the critically loaded case has N 2 O ε −1/2 and leads to a somewhat more difficult problem than either the underloaded or overloaded cases.
We will obtain the PDE for φ 0 ζ, w as a limiting case of the difference equation s satisfied by p m, n , and explicitly solve the PDE by transform methods.We will obtain detailed results for the marginal distributions Prob N 1 m ∞ n 0 p m, n and Prob N 2 n ∞ m 0 p m, n , as well as the mean queue lengths.We will also obtain other approximations to p m, n that are valid on scales where m o ε −1 and/or n o ε −1/2 .In particular, we shall show that p m, n is O ε on the scale m O ε −1/2 and n O 1 .
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 F x, y m,n x m y n p m, n and obtain a functional equation for the boundary values F x, 0 and F 0, y .This is ultimately converted to a Dirichlet problem, which is solved to yield the boundary values of F 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 x 1 − O ε and y 1 − O √ ε , and then inverting the double transform, we also obtain our main approximation p m, n ∼ ε 3/2 φ 0 ζ, w .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., φ 0 ζ, w is somewhat nonstandard in that the density vanishes as ζ → 0 and the approximation breaks down for ζ O √ ε .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 εm and √ εn.
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 p m, n .In Section 3, we summarize all of our main results.The derivations are given in Section 4 for the scale m O ε −1 , n O ε −1/2 and in Section 5 for the other ranges of m, n.
Throughout the paper we will use the notation f x ∼ g x to mean lim x → x 0 f x / g x 1, f x o g x to mean lim x → x 0 f x /g x 0, and f x O g x to mean that |f x /g x | is bounded for x sufficiently close to x 0 .

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 0 < ε 1.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 1/μ and 1, respectively, where μ O 1 .The server works at unit rate, and if neither queue is empty devotes fractions 1 − εκ and εκ of its effort to the primary and secondary queues, respectively, where κ O 1 .The corresponding service rates are 1 − εκ μ 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 ε → 0 with ω, μ, κ and δ fixed, and then σ and λ vary with ε so that 2.1 and 2.2 hold.
Since ω > 0 the system is stable.Let p m, n denote the stationary probability that there are m jobs in the primary queue and n jobs in the secondary queue.Then, the balance equations satisfied by p m, n are λ εσ p 0, 0 μp 1, 0 p 0, 1 .

2.6
The normalization condition is The mean number of jobs in the primary and secondary queues are Advances in Operations Research 5 From Little's Law, the corresponding mean waiting times are E N 1 /λ and E N 2 / εσ .From 8 , the conservation law of Kleinrock implies, using 2.1 and 2.2 , that 2.9

Summary of Results
We consider

3.1
We then have.

3.2
This gives the limiting density of N 1 , N 2 , which applies for fixed values of ζ and w.We next evaluate this density in various limiting cases, such as w → ∞ and δ → ±∞, to gain more insight into its structure, and to verify consistency with results in 16 .
where the complementary error function is given by [25] Erfc

3.8
Remark 3.8.This matches with 16, Result 4 , with τ w In Propositions 3.9 and 3.16 below, we give the limiting marginals of N 2 and N 1 , which apply for w 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 3.17 Corollary 3.17.
Proposition 3.19.The lowest order asymptotic approximation to the stationary mean secondary queue length is
We next give some asymptotic results for φ 0 ζ, w that apply for δ fixed and ζ and/or w → ∞.We also give the "corner" behavior as ζ, w → 0, 0 . dβ.

3.32
Remark 3.23.The results show that for δ < 0 the density φ 0 is asymptotically of product form in a sector, and distinctly nonproduct form in the complimentary sector.For δ > 0 the product form behavior is absent.Item v shows that φ 0 has an integrable singularity near the corner, while item vi shows that φ 0 O ζ as ζ → 0. The second integral in 3.32 may be expressed in terms of a modified Bessel function, using the identity dZ 2UK 1 U .

3.33
The approximation p m, n ≈ ε 3/2 φ 0 ζ, w is only valid for m O ε −1 and n O ε −1/2 .For other ranges of m, n , other expansions must be constructed, and these we summarize below.

3.36
Remark 3.25.We comment that for n O 1 and ζ > 0 m O ε −1 p m, n ∼ ε 3/2 φ 0 ζ, 0 so that the diffusion approximation still applies.For m O 1 and w > 0 item i still applies, and then p m, n ∼ ε 2 P − w , which is independent of m and can be used to estimate the boundary probabilities p 0, n , which are O ε 2 for n O ε −1/2 .Note that in item ii and for n 0 in item iii , p m, n is O ε , which is larger than the order of magnitude O ε 3/2 of the diffusion approximation on the ζ, w scale.But, the total mass in the range in ii is O √ ε , while that in ranges i and iii is O ε .

Analysis of the Main Diffusion Approximation
and the boundary condition We will discuss the second boundary condition, along ζ 0, after 4.16 .

Advances in Operations Research
We let To solve 4.4 and 4.5 we apply a transform in the w 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 C appears only when A < 0 and the term Ae Ax corresponds to a single discrete eigenvalue in the spectral theory.By multiplying 4.7 by e Ax and integrating from x 0 and x ∞, we find that C 2 ∞ 0 F x e Ax dx.Applying 4.6 with B β and x w, we let then integration by parts and the boundary condition 4.5 leads to Hence, from 4.4 , Also, if we use 3.1 , 2.1 , and 2.2 , in 2.5 , we obtain the lowest order boundary condition ∂φ 0 /∂w 0, w 0, for w > 0. We conclude that φ 0 0, w and Φ 0 0, w are proportional to a delta function at w 0 and hence, from 4.8 and 4.16 , that Ω 0 0, β ωβ.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 i/2κ δ 2 κμω 2 to i∞, and let

Advances in Operations Research
For δ < 0 there is a contribution from the pole at β −iδ/ 2κ , and we obtain for ζ > 0 and w > 0, and the boundary condition We let It follows, from 4.15 and 4.24 -4.26 , that

Advances in Operations Research
It follows that  where c is given by 3.21 .Hence,

4.35
The evaluation of the integral in 4.35 is routine, but tedious, and Proposition 3.19 follows from 3.1 , since

4.36
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 β −iδ/ 2κ and branch points at Then, we can represent, for any δ, φ 0 as the contour integral

4.38
Here, C is a horizontal contour in the β-plane, on which The condition in 4.39 insures that if δ < 0 the pole at β i|δ|/ 2κ lies below the contour C.
If δ > 0 we can shift the contour to the real β-axis, and then 4.38 becomes the same as 3.2 .If δ < 0 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 ζ, w → ∞ we employ the saddle point method.There is a saddle point where

4.41
The saddle is on the imaginary axis and the directions of steepest descent are arg β − ib s 0, π.Then shifting the contour C into another horizontal contour through ib s leads to 3.24 .Such a shift is always permissible if δ > 0, but if δ < 0 we need the saddle to lie above the pole, that is, b s > |δ|/ 2κ , and this occurs precisely when w/ζ > |δ|/ μω .We thus obtain the condition in item i of Proposition 3.22.If δ < 0 and w/ζ < |δ|/ μω the pole dominates the saddle point contribution, and we obtain 3.23 .
For ζ → ∞ and w O 1 a different analysis is needed, as now b s → 0 so the saddle approaches the real axis, where the integrand in 4.38 has as simple zero.In this case which applies only if δ > 0 , we shift C back to the real axis and expand the integrand for β → 0. Using

Advances in Operations Research
we thus obtain where Δ ω 2 δ 2 / κμ .Evaluating the integral s in 4.43 leads to 3.25 .If we consider the opposite limit, where ζ O 1 and w → ∞, then the saddle point approaches the upper branch point at ib .But by deforming C to an integral about the branch cut we can show that the final result coincides with the expansion of 3.25 for w/ζ 1, which is given by 3.28 .Next, we consider the corner behavior of φ 0 as ζ, w → 0. Now, the main contribution to the integral will come from where |β| is large.From 3.2 for δ > 0, we then obtain 4.44 which yields 3.30 , and this can be shown to remain valid for δ < 0. Finally, we fix w and let ζ → 0. Simply, setting ζ 0 in 3.2 leads to a divergent integral.However, an integration by parts leads to, for δ > 0, where Δ ω 2 δ 2 /κμ and

4.46
Expanding 4.46 for ζ → 0 and noting that, by contour integration if δ > 0 , we write the integrand as

Analysis of Boundary and Corner Regions
We analyze cases where m o ε −1 and/or n o ε −1/2 .While these carry mass that is asymptotically small, they must be considered to insure that p m, n is properly normalized to higher orders in ε, and to compute higher order approximations to the moments.Also, to determine φ 0 ζ, w we used the boundary condition φ 0 0, w ωδ w , and analysis of cases where ζ and w are small will allow us to examine this condition more carefully.
First we observe from 3.31 that φ 0 ζ, w vanishes linearly as ζ → 0, which indicates a nonuniformity in the asymptotics.We first consider the scale m O 1 with w > 0 and set p m, n ε ν 1 P m, w; ε ∼ ε ν 1 P m, w .

5.1
Here, ν 1 is a constant that will be determined by asymptotic matching.From 2.

5.6
The limiting form of 5.6 as ε → 0 is P SS 0, so we write P S, w SP w P − w .

5.8
To leading order as ε → 0, 5.8 implies that P w 0, w μP S 0, w 0. 5.9 Combining 5.7 and 5.9 , we thus write the approximation on the S, w scale as We then determine ν 2 and P w by asymptotically matching 5. We consider some limiting cases of the S, w scale result in 3.34 .As δ → ∞ from 3.28 we obtain b ∼ δ/ 2κ and then 3.32 yields, after some computations,

5.12
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 p m, n ∼ ε 2 P − w agrees with 5.14 .
We will next consider the scale m, n O 1 and matching this to the S, w or m, w scale s will require the behavior of 3.34 as w → 0. For w → 0 the second integral in the definition of a w in 3.32 dominates, and we obtain a w ∼ 2 w 2 , w −→ 0 , 5.15 so that Also, P − w will be less singular O w −1 than P w , and, since w n/ √ ε, for w → 0 the approximation ε 2 √ εmP w P − w becomes of order O ε 3/2 on the m, n scale.On the scale m, n O 1 we will show that p m, 0 is asymptotically larger than p m, n for n 1, which we just inferred to be O ε 3/2 .We thus set p m, n

5.17
On the m, n scale, the limiting form of 2.3 , with 5.17 Hence, Q m, n is linear in m, and writing we then obtain from 2.5 ,

5.26
It remains to determine Q n and the constant γ.By letting m S/ √ ε and n → ∞ and matching to the S, w scale, we conclude that as n → ∞ While these relations are consistent with 5.21 , they do not determine Q n for n O 1 .
We thus examine another scale, that has m S/ √ ε O ε −1/2 and n O 1 .By matching to 5.26 as m → ∞, we expect that p m, n will be O ε on this scale, so we set p m, n ∼ εQ S, n .

5.28
Then, from 2.3 , we obtain to leading order where ν is a "spectral" parameter.This is easily solved to yield

5.32
By integrating 5.32 over a loop in the complex ν-plane that encircles the branch cut Re ν ∈ 0, 4 , we obtain  For S → ∞ the major contribution to the integral in 5.39 will come from near Ω 0, and the Laplace method yields Q S, n ∼ 2 π μ κ f 0 S .

5.40
Setting w 0 in 3.2 and letting ζ → 0, or, equivalently, using 3.30 , we conclude that f 0 ω.Now, we match the S, n and m, n scales.

5.43
Now, Q , Q − in 5.20 are determined, as is γ in 5.26 , and we have established 3.36 in Proposition 3.24.Also, using f Ω ω cos Ω/2 in 5.39 leads to 3.35 in Proposition 3.24.Thus, we have used a singular perturbation analysis to approximate p m, n on scales other than the ζ, w scale.While these carry mass that is o 1 , the size of p m, n may actually be larger than O ε 3/2 , as is evident from 3.35 .

2 u
cosh u exp − w/2κ δ 2 κμω 2 cosh u − c cosh u c cosh 2 u − c 2 du, 4.34 47 , identifying the O ζ 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.
, becomes 2μQ m, n μQ m − 1, n μQ m 1, n , 5.18 Advances in Operations Research while the boundary condition 2.5 yields

Thus 5 .
.29 while the boundary condition 2.4 leads to μQ SS S, 0 κ Q S, 1 − Q S, 29 corresponds to driftless diffusion in the S variable and a driftless random walk in the discrete n variable.We may replace 5.30 by the "artificial" boundary condition Q S, −1 Q S, 0 , by extending 5.29 to hold also at n 0. To analyze 5.29 consider the discrete Green's function problem G n 1 G n − 1 ν − 2 G n −δ n, n 0 ,
To leading order, 5.4 implies that P 0, w P 1, w so that 5.3 becomes P m, w P − w .But, as m → ∞, ε ν 1 P − w cannot match to ε 3/2 φ 0 ζ, w , since φ 0 vanishes as ζ mε → 0. This indicates that we must analyze another scale, which has m This establishes 3.34 in Proposition 3.24.Expression 5.10 remains valid for m O 1 as then S√ εm → 0 and we obtain p m, n ∼ ε 2 μ ∞ w P u du, which is independent of m and consistent with our previous analysis for m O 1 .