Steady-state analysis of shortest expected delay routing

Queueing Systems, Aug 2016

We consider a queueing system consisting of two nonidentical exponential servers, where each server has its own dedicated queue and serves the customers in that queue FCFS. Customers arrive according to a Poisson process and join the queue promising the shortest expected delay, which is a natural and near-optimal policy for systems with nonidentical servers. This system can be modeled as an inhomogeneous random walk in the quadrant. By stretching the boundaries of the compensation approach we prove that the equilibrium distribution of this random walk can be expressed as a series of product forms that can be determined recursively. The resulting series expression is directly amenable to numerical calculations and it also provides insight into the asymptotic behavior of the equilibrium probabilities as one of the state coordinates tends to infinity.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs11134-016-9497-7.pdf

Steady-state analysis of shortest expected delay routing

Queueing Syst Steady-state analysis of shortest expected delay routing Jori Selen 0 1 2 Ivo Adan 0 1 2 Stella Kapodistria 0 1 2 Johan van Leeuwaarden 0 1 2 B Jori Selen 0 1 2 Stella Kapodistria 0 1 2 0 Department of Industrial Engineering & Innovation Sciences, Eindhoven University of Technology , Eindhoven , Netherlands 1 Department of Mechanical Engineering, Eindhoven University of Technology , Eindhoven , Netherlands 2 Department of Mathematics and Computer Science, Eindhoven University of Technology , Eindhoven , Netherlands We consider a queueing system consisting of two nonidentical exponential servers, where each server has its own dedicated queue and serves the customers in that queue FCFS. Customers arrive according to a Poisson process and join the queue promising the shortest expected delay, which is a natural and near-optimal policy for systems with nonidentical servers. This system can be modeled as an inhomogeneous random walk in the quadrant. By stretching the boundaries of the compensation approach we prove that the equilibrium distribution of this random walk can be expressed as a series of product forms that can be determined recursively. The resulting series expression is directly amenable to numerical calculations and it also provides insight into the asymptotic behavior of the equilibrium probabilities as one of the state coordinates tends to infinity. - Mathematics Subject Classification 60K25 · 68M20 · 60J27 1 Introduction In this paper we analyze the performance of a system with two servers under the shortest expected delay (SED) routing policy. This routing policy assigns an arriving customer to the queue that has the shortest expected delay (sojourn time), where delay refers to the waiting time plus the service time. This policy arises naturally in various application areas, and poses considerable mathematical challenges. In particular, we focus on a queueing system with two servers, where each server has its own queue with unlimited buffer capacity. All service times are independent and exponentially distributed, but the two servers have different service rates, i.e., respectively, 1 and s. In both queues customers are served in FCFS-order. Customers arrive according to a Poisson process with rate λ and upon arrival join one of the two queues according to the following mechanism: Let q1 and q2 be the number of customers in queue 1 and 2, respectively, including a possible customer in service. For an arriving customer, the expected delay in queue 1 is q1 + 1 and in queue 2 is (q2 + 1)/s. The SED routing policy assigns an arriving customer to queue 1 if q1 + 1 < (q2 + 1)/s and to queue 2 if q1 + 1 > (q2 + 1)/s. In case the expected delays in both queues are equal, i.e., q1 + 1 = (q2 + 1)/s, the arriving customer joins queue 1 with probability q and queue 2 with probability 1 − q. Once a customer has joined a queue, no switching or jockeying is allowed. The service times are assumed independent of the arrival process and the customer decisions. This system is stable if and only if, see, for example, [13, Theorem 1], ρ := λ/(1 + s) < 1. (1.1) We refer to this specific queueing system as the SED system. The SED system can be modeled as a two-dimensional Markov process on the states (q1, q2) with the transition rate diagram as in Fig. 1. From the transition rate diagram, it is evident that this state description leads to an inhomogeneous random walk in the quadrant, making an exact analysis extremely difficult. The inhomogeneous behavior occurs along the line s(q1 + 1) = q2 + 1 and it can be expected that the solution structure of the stationary probabilities above and below this line will be different. Moreover, for s = 1, this line divides the state space in unequal proportions, further increasing the complexity of an exact analysis. Under the assumption of identical service rates, SED routing becomes join the shortest queue (JSQ) routing which is known to minimize the mean expected delay [29]. If the service rates of the two servers are different, then JSQ is not optimal. As Whitt [29] points out, if the system if fully observable and we are a priori knowledgeable of the service times of waiting customers, for example, by looking at the shopping baskets in a supermarket, then the natural choice is to join the queue promising the smallest sojourn time in the system, instead of the shorter of the queues. However, this type of information is too detailed and might not always be available. In particular, q1 + 1 = (q2 + 1)/s q2 λ λ(1−q) λ λ λ λ λq λ λ λ λ q1 there are many situations in which we have limited knowledge of the service times of the waiting customers. Such systems include teller waiting lines, production facilities, communication networks, etc. If the only available information is that of the number of waiting customers at each queue, then it is quite natural, although not necessarily optimal, to choose the shortest queue routing [4,5,9,11,16,17,20,23,25,26,30,31]. However, if on top of the number of waiting customers we can estimate the expected service times at the two queues, then the natural choice is to route the customers according to SED, rather than JSQ. SED routing seems to be a natural choice, but in practice this choice is not always optimal, since it does not minimize the mean stationary delay [14,29]. For the two nonidentical server setting, Hajek [22] solves a Markov decision problem and proves that the optimal routing policy is of the threshold type. However, SED routing exhibits relatively good performance at both ends of the range of system utilizations, but performs slightly worse than other policies at medium utilizations; see [7]. Furthermore, Foschini [18] has shown that SED routing is asymptotically optimal and results in complete resource pooling in the heavy traffic limit. State-dependent routing policies such as JSQ and SED are typically difficult to analyze. For the stationary behavior of queueing systems with SED routing very little is known. The difficulty in analyzing this type of model is evident from the analysis of the JSQ policy for two identical parallel exponential servers [21]. The first major steps towards its analysis were made in [16,25], using a uniformization approach that established that the generating function of the equilibrium distribution is meromorphic. Thus, a partial fraction expansion of the generating function of the joint stationary queue length distribution would in principle lead to a representation of the equilibrium distribution as an infinite linear combination of product forms. For an extensive treatment of the JSQ system using a generating function approach, the interested reader is referred to [12,15]. An alternative approach that is not based on generating functions is the compensation approach [4]. This approach directly solves the equilibrium equations and leads to an explicit solution. The essence of the approach is to first characterize the product forms satisfying the equilibrium equations for states in the inner region of the two-dimensional state space and then to use the product forms in this set to construct a linear combination of product forms which also satisfies the boundary conditions. The construction is based on a compensation argument: after introducing the starting product form, new product forms are added so as to alternately compensate for the errors on the two boundaries. The compensation approach has been developed in a series of papers [1,4–6] and aims at a direct solution for a class of two-dimensional random walks on the lattice of the first quadrant that obey the following conditions: (i) Step size: only transitions to neighboring states. (ii) Forbidden steps: no transitions from interior states to the North, North-East, and East. (iii) Homogeneity: the same transitions occur according to the same rates for all interior points, and similarly for all points on the horizontal boundary, and for all points on the vertical boundary. The approach exploits the fact that the balance equations in the interior of the quarter plane are satisfied by linear (finite or infinite) combinations of product forms that need to be chosen such that the equilibrium equations on the boundaries are satisfied as well. As it turns out, this can be done by alternatingly compensating for the errors on the two boundaries, which eventually leads to an infinite series of product forms. The SED queueing system in this paper is more complicated than the JSQ and other classical queueing systems, see, for example, [1,3–6], since the two-dimensional random walk that describes the SED system exhibits inhomogeneous behavior in the interior of the quadrant. In this paper, we show that the compensation approach can nevertheless be further developed to overcome the obstacles caused by the inhomogeneous behavior of the random walk. This leads to a solution for the stationary distribution in the form of a tree of product forms. The only other work in this direction is [3], which considers SED routing for two identical parallel servers with Poisson arrivals and Erlang distributed service times. The crucial difference with our setting is that we do not focus on generalizing service times, but instead consider servers with different service rates. The remainder of the paper is organized as follows. In Sect. 2 we introduce the model in detail and describe the equilibrium equations. We discuss the compensation approach and its methodological extensions together with our contribution in Sect. 3. Some numerical results are presented in Sect. 4. Section 5 applies the compensation approach to determine the equilibrium distribution of the SED system as a series of product-form solutions. Finally, we present some conclusions in Sect. 6. 2 Equilibrium equations The Markov process associated with the SED system has an inhomogeneous behavior in the interior of the quadrant, specifically along the line s(q1 + 1) = q1 + 1; see Fig. 1. In this section we transform the two-dimensional state space (q1, q2) to a half-plane with a finite third dimension. For this state description, we show that the theoretical framework of the compensation approach can be extended and in this way we determine the equilibrium distribution of the SED system. We will henceforth assume that s is a positive integer number. The service rate s could also be chosen to be rational and the analysis would be similar, but notationally more difficult. We further elaborate on this point in Remark 1. In queue 2 we count the number of groups of size s and denote it by j , i.e., j = q2/s , and we denote the number of remaining customers by r , i.e., r = mod(q2, s). Clearly, a single group in queue 2 requires the same expected amount of work as a single customer in queue 1. The total number of customers in queue 2 is thus j s + r and for an arriving customer the expected delay in queue 2 is j + (r + 1)/s. In terms of these variables, SED routing works as follows: if q1 + 1 < j + (r + 1)/s, the arriving customer joins queue 1, and if q1 + 1 > j + (r + 1)/s, the arriving customer joins queue 2. If the expected delays in both queues are equal, i.e., q1 + 1 = j + (r + 1)/s, the arriving customer joins queue 1 with probability q and queue 2 with probability 1 − q. For convenience we introduce the length of the shortest queue as m = min(q1, j ) and the difference between queue 2 and queue 1, i.e., n = j − q1. Using this notation, the SED system is formulated as a three-dimensional Markov process with state space {(m, n, r ) | m ∈ N0, n ∈ Z, r = 0, 1, . . . , s − 1}. Under the stability condition (1.1) the equilibrium distribution exists. Let p(m, n, r ) denote the equilibrium probability of being in state (m, n, r ) and let p(m, n) = ( p(m, n, 0), p(m, n, 1), . . . , p(m, n, s − 1))T (2.1) with xT the transpose of a vector x. Throughout the paper, we use bold lowercase letters or numbers for vectors and uppercase Latin letters for matrices. For convenience, we have listed all state variables and their interpretation in Table 1. Variable q1 q2 q2/s m n r The transition rates are given by ⎧ (m + 1, n − 1, r ), m ≥ 0, n > 0, r = 0, 1, . . . , s − 1, (m, n, r ) −(−1+−−s)→ρ ⎪⎨ (m, n, r + 1), m ≥ 0, n ≤ 0, r = 0, 1, . . . , s − 2, ⎪⎩ (m + 1, n + 1, 0), m ≥ 0, n < 0, r = s − 1, (1+s)ρ(1−q) (m, 0, s − 1) −−−−−−→ (m, 1, 0), m ≥ 0, (1+s)ρq (m, 0, s − 1) −−−−−−→ (m, −1, s − 1), m ≥ 0, corresponding to arrivals, and 1 (m, n, r ) −→ (m − 1, n + 1, r ), m > 0, n ≥ 0, r = 0, 1, . . . , s − 1, (m, n + 1, r ), m ≥ 0, n < 0, r = 0, 1, . . . , s − 1, ⎧ (m, n, r − 1), s ⎪⎨ (m, n − 1, s − 1), (m, n, r ) −→ m ≥ 0, n ∈ Z, r = 1, 2, . . . , s − 1, m ≥ 0, n > 0, r = 0, ⎪⎩ (m − 1, n − 1, s − 1), m > 0, n ≤ 0, r = 0, corresponding to service completions. Figure 2a displays the transition rate diagram for the three-dimensional state space. The transition rates are described by the matrices Ax,y in the positive quadrant and Bx,y in the negative quadrant, where the pair (x , y) indicates the step size in the (m, n)-direction. Let I be the s × s identity matrix, M (x,y) be an s × s binary matrix with element (x , y) equal to one and zeros elsewhere, and L be an s × s subdiagonal matrix with elements (x , x − 1), x = 1, 2, . . . , s − 1 equal to one and zeros elsewhere. For consistency with the indexing of the vector p(m, n), indexing of a matrix starts at 0. The transition rate matrices take the form A1,−1 = (1 + s)ρ I, A−1,1 = B0,1 = I, A0,−1 = B−1,−1 = s M (s−1,0), A0,0 = −(1 + s)(ρ + 1)I + s L T , A0,1 = (1 + s)ρ(1 − q)M (0,s−1), B1,1 = (1 + s)ρ M (0,s−1), B0,−1 = (1 + s)ρq M (s−1,s−1), B0,0 = A0,0 + (1 + s)ρ L . The equilibrium equations can be written in matrix-vector form. We partition the state space as illustrated in Fig. 2b. For the interior of the positive and negative quadrant we have the following inner equations: A0,0p(m, n) + A1,−1p(m − 1, n + 1) B0,0p(m, n) + B1,1p(m − 1, n − 1) + A0,−1p(m, n + 1) + A−1,1p(m + 1, n − 1) = 0, m ≥ 1, n ≥ 2, + B0,1p(m, n − 1) + B−1,−1p(m + 1, n + 1) = 0, m ≥ 1, n ≤ −2. (2.2) (2.3) The equilibrium equations corresponding to the states on the horizontal axis, or directly adjacent to the horizontal axis, are referred to as the horizontal boundary equations and are given by n A−1,1 A0,0 A0,0 + I Finally, for the three remaining boundary states near the origin, we have ( A0,0 + I )p(0, 1) + A0,−1p(0, 2) + A−1,1p(1, 0) + A0,1p(0, 0) = 0, ( B0,0 + s M (0,0))p(0, −1) + B0,1p(0, −2) m ≥ 1, n = 1 m ≥ 1, n = 0 +B−1,−1p(1, 0) + B0,−1p(0, 0) = 0, (B0,0 + I + s M (0,0))p(0, 0) + A0,−1p(0, 1) + B0,1p(0, −1) = 0, corresponding to the states (0, 1), (0, −1) and (0, 0), respectively. s2 can also be Remark 1 (Rational s) A system with a rational service rate s = s1 analyzed. In that case, one needs to consider a system with two servers and service rates s1 and s2. Similar to our analysis at the start of Sect. 2, one denotes the number of groups of size s1 in queue 1 as i and the number of groups of size s2 in queue 2 as j . Then, let rn ∈ {0, 1, . . . , sn − 1} denote the number of remaining customers in queue n = 1, 2. Based on the aforementioned construction, a single group in either queue 1 or 2 requires the same expected amount of work. Lastly, set m = min(i, j ), n = j − i and the third finite dimension is a lexicographical ordering of the states (r1, r2) ∈ {0, 1, . . . , s1 − 1} × {0, 1, . . . , s2 − 1}. This state space description leads to a transition rate diagram that has a similar structure as the one seen in Fig. 2a. In this sense, a system with a rational service rate can be analyzed using the approach described in this paper. (2.10) (2.11) 3 Evolution of the compensation approach In this section we use the following abbreviations: vertical boundary (VB); vertical compensation step (VCS); horizontal boundary (HB); and horizontal compensation step (HCS). The compensation approach is used for the direct determination of the equilibrium distribution of Markov processes that satisfy the three conditions mentioned in Sect. 1. The key idea is a compensation procedure: the equilibrium distribution can be represented as a series of product-form solutions, which is generated term by term starting from an initial solution, such that each term compensates for the error introduced by its preceding term on one of the boundaries of the state space. In this section, we motivate why the SED system requires a fundamental extension of the compensation approach. We do so by first describing the evolution of the compensation approach through a series of models and present for each model the corresponding methodological contribution. Finally, we describe the extension required for the SED system. The compensation approach was pioneered by Adan et al. [4], for a queueing system with two identical exponential servers, both with rate 1, and JSQ routing. Such a queueing system can be modeled as a Markov process with states (q1, q2) ∈ N02, where qi is the number of customers at queue i , including a customer possibly in service. By defining m = min(q1, q2) and n = q2 − q1, one transforms the state space from an inhomogeneous random walk in the quadrant to a random walk in the half-plane that is homogeneous in each quadrant. Since the two quadrants are mirror images of each other, it is not needed to determine the equilibrium probabilities in both quadrants; it suffices to do so in the positive quadrant. The transition rate diagram of the Markov process is shown in Fig. 3a. For the symmetric JSQ model in [4], the initial solution is of the form η0α0m βn, 0 where η0 is a coefficient and α0 and β0 are the compensation parameters that satisfy the kernel equation n n n 1 1 2ρ m m 1 s s 1 2ρ 2ρ I 2sρL I 2sρM(0,s−1) which is obtained by substituting αm βn in the equilibrium equations of the interior of the positive quadrant and dividing by common powers. This initial solution satisfies the equilibrium equations in the interior and on the HB (there is only one such solution). In order to compensate for the error on the VB, one adds the compensation term ν0α1m β0n such that η0α0m β0n + ν0α1m β0n satisfies the equilibrium equations in the interior and on the VB. The compensation parameter α1 with α1 < β0 is generated from (3.1) for a fixed β = β0. The coefficient ν0 satisfies a linear equation and is a function of η0, α0, β0, and α1. The resulting solution violates the equilibrium equations on the HB. Hence, one adds another compensation term η1α1m β1n such that ν0α1m β0n + η1α1m β1n satisfies the equilibrium equations in the interior and on the HB, where β1 and η1 are determined in a similar way as for the VCS. Repeating the compensation steps leads to a series expression for the equilibrium probabilities that satisfies all equilibrium equations: ∞ l=0 ∞ l=0 where C is the normalization constant. Figure 4a displays the way in which the compensation parameters are generated. The first extension is presented in [5] for the asymmetric JSQ model, i.e., the servers are now assumed to be nonidentical with speeds 1 and s. The symmetry argument used earlier does not hold anymore and one needs to consider the complete half-plane; see Fig. 3b for the transition rate diagram. Note that the half-plane consists of two α0,1 α0,1 vertical compensation α1,1 α1,2 α1,1 . . . α1,s vertical compensation (c) horizontal compensation quadrants with different transition rates that are coupled on the horizontal axis. The approach in this case is an extension of the approach introduced in [4]. In a VCS, one compensates solutions that satisfy the positive inner equations on the positive VB as well as solutions that satisfy the negative inner equations on the negative VB. Two kernel equations (one for each quadrant) are used to generate the αs, and the coefficients satisfy different linear equations. For a HCS, each product-form solution that satisfies the positive inner equations generates a single β for the positive quadrant and a single β for the negative quadrant. Accordingly, a product-form solution that satisfies the negative inner equations is compensated on the HB. Thus, in this case the generation of compensation parameters has a binary tree structure; see Fig. 4b. A further extension of the compensation approach is presented in [2] for a model with two identical servers, Erlang-s arrivals and JSQ routing. The state description is enhanced by adding a finite third dimension that keeps track of the number of completed arrival phases. The random walks in the positive and negative quadrant are mirror images, which permits us to perform the analysis only on the positive quadrant; see Fig. 3c for the transition rate diagram. In [2] the authors extend the compensation approach to a three-dimensional setting. Due to the three-dimensional state space, each compensation term takes the form αm βni+(α, β), where i+(α, β) is a vector of coefficients of dimension s (equal to the number of arrival phases). In each HCS, s different parameters β are generated instead of just one. A graphical representation of the generation of compensation terms, which has an s-fold tree structure, is depicted in Fig. 4c. 3.1 Our contribution The model at hand is defined on an s-layered half-plane, thus requiring that we further extend the compensation approach. Similarly to [5], we need to account for the two quadrants by considering two kernel functions (one for each quadrant). Furthermore, in accordance with [2], in every HCS, a total of s different parameters β are generated for the positive quadrant and a single β for the negative quadrant. This leads to a (s +1)fold tree structure for the compensation parameters as depicted in Fig. 5. Additionally, the product-form solutions take the form αm βni+(α, β) or αm β−ni−(α, β) depending on whether they are defined in the positive or negative quadrant, respectively. The resulting solution for the equilibrium distribution is, for m ≥ 0, n ≥ 1, For m ≥ 0, n ≤ −1, where C is the normalization constant and d(i ) := (i − 1)(s + 1). The first subscript l is the level at which a parameter resides, starting at level l = 0 (the initial solution). Within a level, the parameters are differentiated by using an additional index i . A horizontal compensation step and the initial solution have coefficients η and a vertical compensation step has coefficients ν. Additionally, a vector h is generated in each horizontal compensation step. The initial solution is described in Lemma 6, the horizontal compensation step is described in Lemma 8, and the vertical compensation step is described in Lemma 7; see Fig. 5. 4 Numerical results Expression (3.3) is amenable for numerical calculations after applying truncation. For m ≥ 0, n ≥ 1, pL (m, n) = C where the empty sum l−=10 is 0. Here, L = 0 indicates only the initial solution and for instance L = 3 indicates an initial solution, a vertical, horizontal, and another vertical compensation. Naturally, as L increases, the approximation becomes more accurate. We perform several numerical experiments that verify that under SED routing the joint queue length process concentrates between the line where the expected delays in both queues are equal, q1 + 1 = (q2 + 1)/s, and the line where the expected waiting time is equal, q1 = q2/s, using the equilibrium distribution. To this end, we consider a system with service rates 1 and s = 3, q = 0.4 and set L = 16 in (4.1). The For m ≥ 0, n = 0, l=0 i=1 5 q1 10 equilibrium distribution for this model and varying ρ is given in Fig. 6 in the form of a heat plot, supporting our claim. Next, to demonstrate the rate of convergence of the series in (3.3), we derive the number of compensation steps L for which the equilibrium probabilities p(m, n) are considered sufficiently accurate: As a measure of accuracy we compute for each state (m, n) the minimum number of compensation steps L such that Remark 2 (No curse of dimensionality) Take a triangular set of states TM = {(m, n) | m ∈ N0, n ∈ Z, m + |n| ≤ M }, where M is some nonnegative integer. Technically, M needs to be strictly larger than some nonnegative integer N , but we do not go into the details here; the lower bound N is described in Sect. 5.7. For states outside TM , we can use (4.1) to compute p(m, n). Since the number of compensation steps L required to achieve accurate results according to (4.2) decreases with m + |n|, the number of compensation steps L for each state (m, n) ∈/ TM is relatively small. For states (m, n) ∈ TM , a linear system of equilibrium equations needs to be solved to determine the equilibrium probabilities p(m, n), where one uses that p(m, n), (m, n) ∈/ TM , are known. 1 m As an example, for s = 4, ρ = 0.8, and q = 0.4, the choice M = 3 implies that L = 1, which gives accurate results for p(m, n), (m, n) ∈/ TM ; see Fig. 7. So it is evident that the compensation approach does not suffer from the curse of dimensionality. 5 Applying the compensation approach 5.1 Outline In the following subsections we describe the main steps in constructing the equilibrium distribution of the SED system. First, we develop some preliminary results in Sect. 5.2, showing that the inner equations have a product-form solution and we determine one of the two parameters of the product-form solution explicitly. Using these preliminary results, we determine the unique initial solution in Sect. 5.3. The vertical compensation step is outlined in Sect. 5.4. Section 5.5 describes the horizontal compensation procedure. We formalize the resulting solution in terms of a sequence of product forms in Sect. 5.6. This sequence of compensation terms grows as a (s + 1)-fold tree and the problem is in showing that the sequence converges. Section 5.7 is devoted to the issue of convergence. A−1,1 A0,0 Aˆ0,−1 Bˆ0,1 A−1,1 A0,0 A0,−1 A1,−1 A−1,1 Fig. 8 s-layered transition rate diagram of the modified model. Note that the third dimension, i.e., r , is perpendicular to the page. The nonfilled dot, graphed at position (−1, 0), corresponds to the single state (−1, 0, s − 1) 5.2 Preliminary results We conjecture that the inner equations have a product-form solution. To this end, we examine a related model that has the same behavior in the interior as the original model. We then show that the equilibrium distribution of this related model can be expressed as a product-form solution. Moreover, modeling this process as a quasi-birth–and– death (QBD) queue, we obtain a closed form expression for one of the parameters of the product-form solution. This procedure is closely related to the one in [2, Sect. 4]. The related model is constructed as follows. We start from the state space of the original model in Fig. 2a and bend the vertical axis as shown in Fig. 8. Note that for this modified model, m ∈ Z. We add an additional state (−1, 0, s − 1) with a transition rate Aˆ 0,1 = se0 to state (−1, 1) and a transition rate Bˆ0,−1 = (1 + s)ρqes−1 to state (−1, −1), where ei is a column vector of zeros of length s with a one at position i . The transitions from the states on the diagonal m + n = 0, n ≥ 1 are kept consistent with the transitions from a state in the positive interior of the SED model, where the downward transitions are redirected to the state (−1, 0, s − 1), specifically, Aˆ 0,−1 = se0T . Similarly, the transitions from the states on the diagonal m − n = 0, n ≤ −1 are kept consistent with the ones in the negative interior of the SED model, where the upward transitions are redirected to the state (−1, 0, s − 1), specifically, Bˆ0,1 = 1T , where 1 is a column vector of ones of size s. The characteristic feature of the modified model is that its equilibrium equations for m + |n| = 0, |n| ≥ 2 are exactly the same as the ones in the interior and the equilibrium equations for n ∈ {−1, 0, 1} are exactly the same as the ones on the horizontal boundary of the original model. In this sense, the modified model has no “vertical boundary” equations. We next present two lemmas. The first lemma states that a product-form solution exists for the modified model, while in the second lemma we identify the geometric term of the product-form expression. Let pˆ(m, n) = ( pˆ(m, n, 0), pˆ(m, n, 1), . . . , pˆ(m, n, s − 1))T denote the equilibrium distribution of the modified model. Lemma 1 For ρ < 1, the equilibrium distribution of the modified model exists and is of the form pˆ(m, n) = αm qˆ (n), m + |n| ≥ 0, n ∈ Z, with α ∈ (0, 1) and qˆ (n) = (qˆ (n, 0), qˆ (n, 1), . . . , qˆ (n, s − 1))T such that ∞ Proof First notice that the modified model is stable whenever the original model is stable and that pˆ(m, n) satisfies the inner and horizontal boundary equations (2.2)–(2.6) for all m +|n| = 0, except for state (0, 0). Observe that the modified model restricted to the area {(m, n) | m ∈ N0, n ∈ Z, m + |n| ≥ n0} ∪ {(n0 − 1, 0, s − 1)}, n0 = 1, 2, . . . embarked by two lines parallel to the diagonal axes yields exactly the same process. Hence, we can conclude that pˆ (m + 1, n) = αpˆ (m, n), m + |n| ≥ 0, n ∈ Z, (5.1) (5.2) (5.3) (5.4) (5.5) and therefore Finally, we observe that pˆ(m, n, r ) = αm qˆ (n, r ), m + |n| ≥ 0, n ∈ Z, r = 0, 1, . . . , s − 1. ∞ n=−∞ α−|n|qˆ (n, r ) = ∞ n=−∞ pˆ(−|n|, n, r ) < 1, which concludes the proof. In Lemma 1 we have shown that the equilibrium distribution of the modified model has a product form which is unique up to a positive multiplicative constant. In the next lemma we determine the unique α in (5.1). More concretely, the unique α is equal to ρ1+s . Lemma 2 For ρ < 1, the equilibrium distribution of the modified model is of the form pˆ (m, n) = ρ(1+s)m qˆ (n), m + |n| ≥ 0, n ∈ Z. (5.6) k = (1 + s)m + sn + r, n ≥ 0, (1 + s)m − n + r, n < 0. Then, the number of customers in the system forms a Markov process and for all states k > s, the transitions are given by (i) from state k to state k + 1 with rate (1 + s)ρ; and (ii) from state k + 1 to state k with rate 1 + s. Let pˆk denote the probability of having k customers in the system. From the balance principle between states k and k + 1 we obtain pˆk+1 = ρ pˆk , k > s. Furthermore, pˆk+s+1 = = = α = α pˆk . (m,n,r) (1+s)m+sn+r=k+s+1 (m,n,r) (1+r)(m−1)+sn+r=k (l,n,r) (1+s)l+sn+r=k Proof Let k denote the total number of customers in the system, i.e., Combining the last two results immediately yields α = ρ1+s . Combining Lemmas 1 and 2 yields the following result for the original model. Proposition 1 For ρ < 1, the inner and horizontal boundary equations (2.2)–(2.6) have a unique solution of the form p(m, n) = ρ(1+s)m q(n), m ≥ 0, n ∈ Z, with q(n) = (q(n, 0), q(n, 1), . . . , q(n, s − 1))T nonzero and pˆ(m, n, r ) + αm qˆ (n, h) + (m,n,r) (1+s)m−n+r=k+s+1 (m,n,r) (1+s)(m−1)−n+r=k (l,n,r) (1+s)l−n+r=k αl qˆ (n, h) + α αl qˆ (n, r ) pˆ(m, n, r ) αm qˆ (n, r ) (5.7) (5.8) (5.9) (5.10) (5.11) ∞ n=−∞ ρ−(1+s)|n|q(n, r ) < ∞, r = 0, 1, . . . , s − 1. The solution obtained in Proposition 1 satisfies the inner and horizontal boundary equations. However, we still need to specify the form of the vector q(n). It will become apparent, from Lemmas 3, 4 and 5, that the form of the vector q(n) is entirely different in the positive quadrant, on the horizontal axis, and the negative quadrant. Correctly identifying q(n) will result in the initial solution that satisfies the equilibrium equations of the interior and the horizontal axis. In the following lemmas we describe the form of a solution satisfying the inner equations. D+(α, β)i+(α, β) = 0, with D+(α, β) = αβ A0,0 + αβ2 A0,−1 + α2 A−1,1 + β2 A1,−1 and the eigenvector i+(α, β) = (i+(α, β, 0), i+(α, β, 1), . . . , i+(α, β, s − 1))T satisfies, for r = 0, 1, . . . , s − 1, i+(α, β, r ) i+(α, β, 0) = . (ii) The product form p(m, n) = αm β−n i−(α, β), m ≥ 0, n ≤ −1 is a solution of the inner equations of the negative quadrant (2.3) if D−(α, β)i−(α, β) = 0, with D−(α, β) = αβ B0,0 + αβ2 B0,1 + α2 B−1,−1 + β2 B1,1 and the eigenvector i−(α, β) = (i−(α, β, 0), i−(α, β, 1), . . . , i−(α, β, s − 1))T satisfies, for r = 0, 1, . . . , s − 1, i−(α, β, r ) i−(α, β, 0) = Ψ (α, β, ψ−(β))ψ+(β)r − Ψ (α, β, ψ+(β))ψ−(β)r Ψ (α, β, ψ−(β)) − Ψ (α, β, ψ+(β)) , with ψ±(β) = and (1 + s)(ρ + 1) − β ± (β − (1 + s)(ρ + 1))2 − 4s(1 + s)ρ 2s β Ψ (α, β, ψ ) = β − (1 + s)(ρ + 1) + sψ + α (1 + s)ρψ s−1 = (1 + s)ρψ −1 β ψ s − 1 . α Proof (i) Inserting the product form p(m, n) = αm βni+(α, β) into (2.2) and dividing by αm−1βn−1 results in (5.12). For det(D+(α, β)) = 0, the rank of D+(α, β) is s − 1 and thus we have one free variable, which allows us to express i+(α, β, r ) in terms of i+(α, β, 0). The eigenvector i+(α, β) follows from the system of linear equations (5.12), which reads, for r = 0, 1, . . . , s − 2, α2 +β2(1+s)ρ −αβ(1+s)(ρ +1) i+(α, β, r )+αβsi+(α, β, r +1) = 0, (5.20) and gives (5.14). (5.12) (5.13) (5.14) (5.15) (5.16) (5.17) (5.18) (5.19) We construct a solution i−(α, β, r ) = (c+ψ+(β)r + c−ψ−(β)r )i−(α, β, 0) and thus require c+ + c− = 1. The two roots ψ±(β) follow from substituting i−(α, β, r ) = ψr i−(α, β, 0) in (5.22) and dividing by ψr−1i−(α, β, 0), yielding sψ 2 + (β − (1 + s)(ρ + 1))ψ + (1 + s)ρ = 0. The constants c+ and c− follow from (5.21) and the simplification follows from the fact that i−(α, β, r ) = ψr i−(α, β, 0) satisfies (5.23). Remark 3 (Normalized eigenvectors) Without loss of generality we henceforth set the first elements i+(α, β, 0) = i−(α, β, 0) = 1 for all α and β. Since normalization of the equilibrium distribution follows afterwards, one can arbitrarily select the value of the first element of the eigenvectors i+(α, β) and i−(α, β). We wish to determine the αs and βs with 0 < |α|, |β| < 1, for which (5.12) and (5.15) have nonzero solutions i+(α, β) and i−(α, β), respectively. Equivalently, we determine αs and βs that satisfy det(D+(α, β)) = 0 or det(D−(α, β)) = 0. We first investigate the location and the number of zeros of det(D+(α, β)) = 0. Lemma 4 The equation det(D+(α, β)) = 0 assumes the form αβ(1 + s)(ρ + 1) − β2(1 + s)ρ − α2 s − β(αβs)s = 0 and can be rewritten as where ui is the i -th root of unity of us = 1 and β1/s is the principal root. Proof Dividing (5.24) by (αβs)s and taking the s-th root reduces (5.24) to (5.25). (i) The proof consists of three main steps. Furthermore, by using (5.15) for the summands one obtains (5.44). Equation (5.45) follows from substituting (5.43) into (2.6) and using (5.44). One can arbitrarily choose ηs+1 due to the normalization that follows at the end. 5.4 Compensation on the vertical boundary The initial solution (5.43) does not satisfy the positive and negative vertical boundary. We consider a term ηαm βni+(α, β) that satisfies the positive inner balance equations and stems from a solution that satisfies both the inner and horizontal boundary equations. For now, we refer to this term as the original term. In the case of the initial solution (5.43), this would be one of the s product-form solutions of the positive quadrant. The idea behind the compensation approach is to add compensation terms is=1 νi αim βni+(αi , β) to the original term, such that the linear combination ηαm βni+(α, β) + νi αim βni+(αi , β), m ≥ 0, n ≥ 1, (5.47) s i=1 s i=1 satisfies both (2.7) and (2.2). Substituting (5.47) in (2.7) gives ηβn−1V+(α, β)i+(α, β) + νi βn−1V+(αi , β)i+(αi , β) = 0, n ≥ 2, (5.48) where V+(α, β) = β( A0,0 + I ) + β2 A0,−1 + α A−1,1. Note that we indeed require that the original term and the compensation terms share the same β and therefore we obtain the s roots α1, α2, . . . , αs with |αi | < |β| from (5.24). Clearly, (5.47) now satisfies the positive inner equations (2.2). This leaves the s coefficients νi to satisfy the s equations (5.48). However, from Remark 4 we deduce that there exists a specific i for which i+(α, β) = i+(αi , β) and thus there is only one coefficient that is required to be nonzero. Let us now consider a term from the negative quadrant ηαm β−ni−(α, β) that satisfies the negative inner balance equations and stems from the same solution that satisfies both the inner and horizontal boundary equations. Let us now refer to this term as the original term. In the case of the initial solution (5.43), this would be the solution on the negative quadrant. As before, the compensation approach dictates adding a m compensation term νs+1αs+1β−ni−(αs+1, β) to the original term, such that the linear combination ηαm β−ni−(α, β) + νs+1αs+1β−ni−(αs+1, β), m ≥ 0, n ≤ −1, m (5.49) satisfies both (2.8) and (2.3). Substituting (5.49) in (2.8) gives ηβ−n+1V−(α, β)i−(α, β) + νs+1β−n+1V−(αs+1, β)i−(αs+1, β) = 0, n ≤ −2, (5.50) where V−(α, β) = β(B0,0 +s M (0,0))+β2 B0,1+α B−1,−1. Note that we indeed require that the original term and the compensation term share the same β and therefore we obtain the root αs+1 with |αs+1| < |β| from (5.35). This ensures that the linear combination (5.49) satisfies the negative inner equations (2.3). Even though there is only one coefficient νs+1 for s equations, it turns out that this provides enough freedom to satisfy (5.50). Lemma 7 (Vertical compensation) (i) Consider the product form ηαm βni+(α, β) that satisfies the positive inner equations (2.2) that stems from a solution that satisfies the inner and horizontal boundary equations. For this α and fixed β, let α1 be the root that satisfies i+(α, β) = i+(α1, β) with |α1| < |β|. Then there exists a coefficient ν1 such that p(m, n) = ηαm βni+(α, β) + ν1α1m βni+(α1, β), m ≥ 0, n ≥ 1, satisfies (2.2) and (2.7). The coefficient ν1 satisfies ν1 = −η 1 − βα (1 + s)ρ 1 − αβ1 (1 + s)ρ . (ii) Consider a product form ηαm βni−(α, β) that satisfies the negative inner equations (2.3) that stems from a solution that satisfies the inner and horizontal boundary equations. For this β, let αs+1 be the root of (5.35) with |αs+1| < |β| and let i−(αs+1, β) be the corresponding nonzero vector satisfying (5.17). Then there exists a coefficient νs+1 such that p(m, n) = ηαm β−ni−(α, β) + νs+1αs+1β−ni−(αs+1, β), m ≥ 0, n ≤ −1, m satisfies (2.3) and (2.8). The coefficient νs+1 is given by νs+1 = −η s − βα (1 + s)ρi−(α, β, s − 1) β (1 + s)ρi−(αs+1, β, s − 1) s − αs+1 . Proof (i) Use (5.12) to establish that V+(α, β)i+(α, β) = β2 β β I − α A1,−1 i+(α, β) = β 1 − α (1 + s)ρ i+(α, β). (5.55) Then using Remark 4 we can establish that there exists a single νi that is nonzero. We label the single nonzero coefficient as ν1 and find it from (5.48) with the help of (5.55). (5.51) (5.52) (5.53) (5.54) (ii) Use (5.15) to establish that β V−(α, β)i−(α, β) = β s M (0,0) − α (1 + s)ρ M (0,s−1) i−(α, β). (5.56) Thus, (5.50) reduces to a single equation and νs+1 follows directly. 5.5 Compensation on the horizontal boundary The solution obtained after compensation on the vertical boundary, as outlined in Lemma 7, does not satisfy the horizontal boundary equations. So, we need to compensate for the error on the horizontal boundary by adding new terms. The compensation procedure on the horizontal boundary has a few differences from the one described in the previous section. We outline these differences by informally treating the compensation of a product-form term of the positive quadrant. The difference in the compensation procedure is due to the fact that adding compensation terms only for the positive quadrant does not make the solution satisfy the horizontal boundary equations. Intuitively, this originates from the fact that the horizontal boundary, i.e., n = 0, is connected to both the positive and negative quadrants. Thus, for the horizontal compensation step, we need to add product-form terms for the complete positive half-plane. It turns out that these product-form terms are nearly identical to the initial solution, which indeed satisfied both the inner and horizontal boundary equations. The same procedure holds for a product-form term of the negative quadrant. The formal compensation on the horizontal boundary is outlined in the next lemma. Lemma 8 (Horizontal compensation) (i) Consider the product form ναm βni+(α, β) that satisfies the positive inner equations (2.2) that stems from a solution that satisfies the inner and vertical boundary equations. For this α, let β1, β2, . . . , βs be the roots of (5.24) with |βi | < |α|, i = 1, 2, . . . , s and let i+(α, βi ) be the corresponding nonzero vectors satisfying (5.14). Symmetrically, for this α, let βs+1 be the root of (5.35) with |βs+1| < |α| and let i−(α, βs+1) be the corresponding nonzero vector satisfying (5.17). Then there exists a nonzero vector h and coefficients η1, η2, . . . , ηs+1 such that ⎧ ναm βni+(α, β) p(m, n) = ⎪⎪⎪⎪⎨⎪⎪ + s ηi αm βini+(α, βi ), m ≥ 0, n ≥ 1, i=1 ⎪⎪⎪ αm h, m ≥ 0, n = 0, ⎪ ⎪⎩⎪ ηs+1αm βs−+n1i−(α, βs+1), m ≥ 0, n ≤ −1, satisfies (2.2)–(2.6). (5.57) The vector h and the coefficients η1, η2, . . . , ηs+1 satisfy αsh(0) + (1 + s)ρqh(s − 1) − ηs+1αs = 0. ⎪ [b]ναm β−ni−(α, β) ⎪⎪⎩ + ηs+1αm βs−+n1i−(α, βs+1) , m ≥ 0, n ≤ −1, The vector h and the coefficients η1, η2, . . . , ηs+1 satisfy i=1 +ηs+1βs+1 B1,1 + α B0,1 i−(α, βs+1) = −νβ B1,1 + α B0,1 i−(α, β), (5.63) αsh(0) + (1 + s)ρqh(s − 1) − ηs+1αs = ναs. Proof (i) Equations (5.58) and (5.60) follow from exactly the same reasoning as outlined in the proof of Lemma 6. Equation (5.60) follows from substituting (5.57) in the negative horizontal boundary equations (2.5) and using (5.15). Note that this indeed reduces to a single equation. (ii) The proof is identical to the proof of (i). (ii) For m ≥ 0, (iii) For m ≥ 0, n ≤ −1, p(m, 0) ∝ 5.6 Constructing the equilibrium distribution The compensation approach ultimately leads to the following expressions for the equilibrium distribution of the SED system, where the symbol ∝ indicates proportionality. (i) For m ≥ 0, n ≥ 1, p(m, n) ∝ νl+1,i(s+1)αl+1,i(s+1)βl−,in(s+1)i−(αl+1,i(s+1), βl,i(s+1)). m We briefly describe the indexing of the compensation terms, which grows as a tree. We indicate the level at which a parameter resides with l, starting at level l = 0. Within a level, we differentiate between parameters by using an additional index i . The procedure for generating terms is as follows: Recall the definition d(i ) := (i − 1)(s + 1) and note that d(i ) + s + 1 = i (s + 1). Informally, the indexing of subsequent terms is visualized in Fig. 9. Note that the first term in the compensation approach is α0,1 = ρ1+s , c.f. Lemma 6. 5.7 Absolute convergence In this subsection we prove that the series for the equilibrium probabilities, as formulated in (5.65), are absolutely convergent. Before we are able to do that, we need some preliminary results. Specifically, we establish that the sequences {αl,i }l∈N0, i=1,2,...,(s+1)l and {βl,i }l∈N0, i=1,2,...,(s+1)l+1 decrease exponentially fast and uniformly in the levels of the parameter tree. Furthermore, we investigate the asymptotic behavior of (ratios of) the parameters αl,i , βl,i , the coefficients, and the elements of the (eigen)vectors. Corollary 1 Let α¯l := maxi=1,2,...,(s+1)l |αl,i |, β¯l := maxi=1,2,...,(s+1)l+1 |βl,i |. Then, (i) ρ1+s = |α0,1| > β¯0 > α¯ 1 > β¯1 > α¯ 2 > β¯2 > · · · (ii) There exists c ∈ (0, 1), such that 0 < α¯l , β¯l < cl . Proof (i) Follows immediately from Lemmas 4 and 5. (ii) For a fixed α, let β1, β2, . . . , βs be the roots of (5.24) with |βi | < |α|, i = 1, 2, . . . , s and let βs+1 be the root of (5.35) with |βs+1| < |α|. We define t (α) := i=1,m2,a..x.,s+1 |βi /α| (5.66) and let C = {α ∈ C | |α| ≤ α0,1 = ρ1+s }. Using Lemmas 4 and 5 we have t (α) < 1 and in particular, since ρ ∈ (0, 1), we have that C is a closed proper subset of the unit disk and thus c1 := maxα∈C t (α) < 1. One can reason along the same lines and obtain a bound c2 for a fixed β. The exponential decrease in terms of the level l follows from the fact that each β that is generated from an α satisfies |β| < |α|c1 and each α that is generated from a β satisfies |α| < |β|c2. Thus, we define c := c1c2. In Corollary 1 we established that both sequences {αl,i }l∈N0, i=1,2,...,(s+1)l and {βl,i }l∈N0, i=1,2,...,(s+1)l+1 tend to zero as l → ∞. Thus, letting α ↓ 0 or β ↓ 0 is equivalent to letting l → ∞. In what follows, whenever the specific dependence on the level l is not needed, we opt to use the simpler notation of Sects. 5.2–5.5. The next lemma presents the limiting behavior of the compensation parameters and their associated eigenvectors. Lemma 9 (i) For a fixed α, let β1, β2, . . . , βs be the roots of (5.24) with |βi | < |α|, i = 1, 2, . . . , s and let βs+1 be the root of (5.35) with |βs+1| < |α|. Then, as α ↓ 0, (a) the ratio βi /α → v−, i = 1, 2, . . . , s with v− < 1, where v− is the root of v2(1 + s)ρ − v(1 + s)(ρ + 1) + 1 = 0. (5.67) (b) the eigenvector i+(α, βi ) → e0, i = 1, 2, . . . , s. Lemma 8(i) or (ii) Lemma 7(ii) (c) the ratio βs+1/α → w− with w− < 1, where w− is the root of w2((1 + s)ρ)s − wss (ψ+(0)s + ψ−(0)s ) + ss = 0, (5.68) where ψ±(0) is defined in (5.18). (d) for r = 0, 1, . . . , s − 1, the elements of the eigenvector i−(α, βs+1, r ) → ψ+(0)r . (ii) For a fixed β, let α1, α2, . . . , αs be the roots of (5.24) with |αi | < |β|, i = 1, 2, . . . , s and let αs+1 be the root of (5.35) with |αs+1| < |β|. Then, as β ↓ 0, (a) the ratio αi /β → 1/v+, i = 1, 2 . . . , s with v+ > 1, where v+ is the root of (5.67). (b) the eigenvector i+(αi , β) → e0, i = 1, 2 . . . , s. (c) the ratio αs+1/β → 1/w+ with w+ > 1, where w+ is the root of (5.68). (d) for r = 0, 1, . . . , s − 1, the elements of the eigenvector i−(αs+1, β, r ) → ψ−(0)r . Proof (i)(a) Set v = βi /α and let α ↓ 0 in (5.26) to establish (5.67). The roots of (5.67) are defined in (5.27) and satisfy 0 < v− < 1 and v+ > 1. (i)(b) We have that i+(α, β, 0) = 1 for all α and for α ↓ 0 the other elements in (5.14) go to ((1 + s)(ρ + 1) − v−(1 + s)ρ − 1/v−)/s, which is equal to zero by (i)(a). (i)(c) Set w = βs+1/α, divide (5.35) by α2 and let α ↓ 0 to establish (5.68). We note that ψ±(·) is a continuous function, so that indeed for α ↓ 0, ψ±(αw) → ψ±(0). Furthermore, (5.68) is a quadratic equation in w with roots |w−| < 1 and |w+| > 1, where the inequalities follow from Lemma 5, and also satisfy w−, w+ ∈ R+. (i)(d) Again, set w = βs+1/α. Recall the definition of the eigenvector in (5.17). We aim at showing that for α ↓ 0, Ψ (α, βs+1, ψ−(αw)) → 0 and Ψ (α, βs+1, ψ+(αw)) tends to some nonzero constant. As α ↓ 0, we have that Ψ (α, βs+1, ψ ) → (1 + s)ρψ −1 w−ψ s − 1 . (5.69) Note that ψ−(0) ∈ (0, 1) and ψ+(0) > 1. Since w− < 1, the limiting value in (5.69) can only equal zero in the case that w− = 1/ψ+(0)s . Thus, we have established that Ψ (α, βs+1, ψ−(αw)) tends to a nonzero constant as α ↓ 0. We next verify that w− = 1/ψ+(0)s is a solution to (5.68), which proves that Ψ (α, βs+1, ψ+(αw)) → 0 as α ↓ 0. Substituting w = 1/ψ+(0)s in the left-hand side of (5.68) yields ψ+(10)2s ((1 + s)ρ)s − ψ+1(0)s ss (ψ+(0)s + ψ−(0)s ) + ss , 1 = ψ+(0)2s ((1 + s)ρ)s − ss (ψ+(0)ψ−(0))s = 0, = ψ+(10)2s ((1 + s)ρ)s − ss ψ+(0)s (ψ+(0)s + ψ−(0)s ) + ss ψ+(0)2s , (5.70) where we used ψ+(0)ψ−(0) = (1 + s)ρ/s, as found in (5.23). (ii) The proof is identical to the proof of (i). Finally, we describe the limiting behavior of the coefficients. In the following lemma we introduce the variable γ associated with a product-form solution that satisfies the positive inner equations, and the variable θ associated with a product-form solution that satisfies the negative inner equations. Lemma 10 (i) Consider the setting of Lemma 7(i). Then, as β ↓ 0, ν1 1 − v−(1 + s)ρ η → − 1 − v+(1 + s)ρ =: γν . (ii) Consider the setting of Lemma 7(ii). Then, as β ↓ 0, lim α↓0 ν |a j (r )| =: γη, (b) with (5.71) (5.72) (5.73) (5.74) (5.75) (5.76) (5.77) (5.78) (5.79) (5.80) where a j = (a j (0), a j (1), . . . , a j (s − 1))T is the solution to the following linear system of equations for a specific j , with b j a column vector of size s with unknowns, v−W a j + Lb j = −v+w j − γηs+1 w−ψ+(0)s−1e0, − W a j + (1 + s)ρ(1 − q)M (0,s−1)b j = w j w j := (iv) Consider the setting of Lemma 8(ii). Then, as α ↓ 0, (a) h ν → 0. v−W c + Ld = − w+ψ−(0)s−1 + θηs+1 w−ψ+(0)s−1 e0, − W c + (1 + s)ρ(1 − q)M (0,s−1)d = 0. Proof (i) Divide both sides of (5.52) by η and let α ↓ 0. (ii) The proof is identical to the proof of (i). (iii) Due to the length of this proof, the proof has been relegated to Appendix 1. (iv) The proof is identical to the proof of (iii). We have established the limiting behavior of ratios of compensation parameters, coefficients, eigenvectors, and the vector on the horizontal axis. As we have seen in Lemma 9, the limiting values of the eigenvectors i+(α, β) and i−(α, β) are finite. So, we can bound the absolute value of both of them by a constant not depending on α or β. In doing so, one does not need to take into account the eigenvectors i+(α, β) and i−(α, β) to establish absolute convergence. Based on this observation and the asymptotic results derived above, we now show that the series appearing in (5.65) are absolutely convergent. Theorem 2 There exists a positive integer N such that (i) the series converges for all m ≥ 0, |n| ≥ 1 with m + |n| > N . (ii) the series ∞ (s+1)l+1 m αl,i |hl,i | converges for all m ≥ N . (5.81) (5.82) (5.83) (5.84) (5.85) (5.86) (5.87) (iv) the series m+|n|>N p(m, n, r ), r = 0, 1, . . . , s − 1 converges absolutely, where p(m, n, r ) is given in (5.65). Proof For this proof, we introduce a partitioning of the elements in a matrix of size s + 1, namely I = {( j, k) | j = 1, 2, . . . , s, k = 1, 2, . . . , s}, H = {( j, k) | j = s + 1, k = 1, 2, . . . , s}, V = {( j, k) | j = 1, 2, . . . , s, k = s + 1}, and E = {( j, k) | j = s + 1, k = s + 1}. (i) We can view the series as an infinite tree with s + 1 roots and each term has s + 1 children. We define the ratios of a term with each of its descendants as Rl(,1i), j→k (m, n) := m n |ηl+1,d(d(i)+ j)+k | αl+1,d(i)+ j β| | l+1,d(d(i)+ j)+k m n |ηl,d(i)+ j | αl,i βl|,d|(i)+ j . (5.89) n n m+|n| We multiply the ratio Rl(,1i), j→k (m, n) by ννll,,dd((ii))++ jj ααll||++n||11,,dd((ii))++ jj ααll||,,nii|| ββllm,,dd+((ii|n))++| jj , yielding (5.88) Rl(,1i), j→k (m, n) = ηl+1,d(d(i)+ j)+k νl,d(i)+ j αlm++1,|nd|(i)+ j νl,d(i)+ j ηl,d(i)+ j βlm,d+(i|n)+| j βlm,d+(i|n)+| j βl|+|1,d(d(i)+ j)+k αl|,i| n n × αlm,i+|n| αl|+n|1,d(i)+ j βl|,nd|(i)+ j . (5.90) As l → ∞ we obtain by Lemmas 9 and 10 that liml→∞ Rl(,1i), j→k (m, n) ≤ R(j1→) k (m, n), where the inequality is element-wise, and R(j1→) k (m, n) = ⎧ |γη||γν ||v−/v+|m+|n|, ( j, k) ∈ I, ⎪⎪⎪⎨ |θη||θνs+1 ||v−||n||w−/w+|m |1/w+||n|, ( j, k) ∈ H, ⎪ |γηs+1 ||γν ||v−/v+|m |1/v+||n||w−||n|, ( j, k) ∈ V, ⎪⎪⎩ |θηs+1 ||θνs+1 ||w−/w+|m+|n|, ( j, k) ∈ E . The convergence of the series is determined by the spectral radius of the corresponding matrix of ratios R(1)(m, n) := (R(j1→) k (m, n)) j,k=1,2,...,s+1. If the spectral radius of the matrix R(1)(m, n) is less than one, the series converges. For more details, see, for example, [5, Sect. 9] or [2, proof of Theorem 5.2]. Observe that |v−|, |1/v+|, |w−|, |1/w+| < 1, but the values of the limiting ratios |γν |, |θνs+1 |, |γη|, |θη|, |γηs+1 |, |θηs+1 | can be larger than one. The spectral radius of the matrix R(1)(m, n) depends on m + |n| and thus we can find an integer N with m + n > N for which the spectral radius is less than one, ensuring the series converges for m + |n| > N . (5.91) (ii) Using a similar analysis as in (i), we define the ratio as Rl(,2i), j→k (m, n) := m n |νl+2,d(d(i)+ j)+k | αl+2,d(d(i)+ j)+k β| | l+1,d(d(i)+ j)+k . m n |νl+1,d(i)+ j | αl+1,d(i)+ j βl|,d|(i)+ j For l → ∞ we have that liml→∞ Rl(,2i), j→k (m, n) ≤ R(j2→) k (m, n) and ⎧ |γη||γν ||v−/v+|m+|n|, ( j, k) ∈ I, ⎪⎪⎨⎪ |θη||γν ||v−||n||v−/v+|m |1/w+||n|, ( j, k) ∈ H, ⎪ |γηs+1 ||θνs+1 ||1/v+||n||w−||n||w−/w+|m , ( j, k) ∈ V, ⎪⎪⎩ |θηs+1 ||θνs+1 ||w−/w+|m+|n|, ( j, k) ∈ E . Lemmas 10(iii)(a) and (iv)(a) show that |hl,i |/|νl,i | → 0 and thus we can bound it from above and need not consider it when proving convergence of the series. Proving convergence of The spectral radius of the matrix R(2)(m, n) := (R(j2→) k (m, n)) j,k=1,2,...,s+1 depends on m + |n| and thus we can find an integer N with m + n > N for which the spectral radius is less than one. (iii) We can rewrite (5.87) as m m αl,i |hl,i | = α0,1 |h0,1| + R(j2→) k (m, n) = establishes convergence of (5.87). Exploiting the similarity with (5.86), we define the ratio R(3) l,i, j→k (m) := Rl(,2i), j→k (m, 0). Hence, the limiting ratios can be bound from above: liml→∞ Rl(,3i), j→k (m) ≤ R(j3→) k (m) = R(j2→) k (m, 0). The spectral radius of the matrix R(3)(m) := (R(j3→) k (m)) j,k=1,2,...,s+1 depends on m and thus we can find an integer N with m ≥ N for which the spectral radius is less than one. (iv) This follows straightforwardly, along the same lines as in [5, Sect. 9]. Remark 6 We note that the series in Theorem 2(i) (without the absolute values) corresponds to the sum of the first series in (5.65a) and the first series in (5.65c). The series in Theorem 2(ii) (without the absolute values) corresponds to the sum of the second series in (5.65a) and the second series in (5.65c). So, proving convergence of the series in Theorem 2 establishes convergence of the series in (5.65). Remark 7 The index N is the minimal nonnegative integer for which the spectral radii of the matrices R(1)(m, n) and R(2)(m, n) are both less than one for m + |n| > N and the spectral radius of R(3)(m) is less than one for m ≥ N . So, for all states (m, n) (5.92) (5.93) (5.94) (5.95) Table 2 The index N for fixed q = 0.4 and varying s and ρ s ρ N 2 with m + |n| > N and (N , 0) the series in (5.65) converges. In general, the index N is small. In Table 2 we list the index N for fixed q, while varying the values of s and ρ. Note that the area of convergence might be bigger than the one based on N , since N is based on R(1)(m, n), R(2)(m, n), and R(3)(m) which are upper bounds on the rate of convergence. We are now in the position to formulate the main result of the paper. Theorem 3 For all states (m, n, r ), m ∈ N0, n ∈ Z, r = 0, 1, . . . , s − 1 and m + |n| > N including m = N , n = 0, the equilibrium probabilities p(m, n) are proportional, up to a multiplicative constant C , to the expressions in (5.65), see (3.3), where C is the normalization constant of the equilibrium distribution. The remaining p(m, n), m + |n| ≤ N are determined by the finite system of equilibrium equations for the states (m, n) with m + |n| ≤ N , where N is the minimal nonnegative integer for which the spectral radii of the matrices R(1)(m, n) and R(2)(m, n) are both less than one for m + |n| > N and the spectral radius of R(3)(m) is less than one for m ≥ N . Proof This proof is similar to the proof in [2, proof of Theorem 5.3] and [5, Sect. 11], but we include it here for completeness. Define LN = {(m, n) | m ∈ N0, n ∈ Z, m + |n| > N } ∪ {(N , 0, s − 1)} and note the similarity with the set of states defined in the proof of Lemma 1 and the associated transition rate diagram in Fig. 8. Then LN is a set of states for which the series in (5.65) converge absolutely. The restricted stochastic process on the set LN is an irreducible Markov process, whose associated equilibrium equations are identical to the equations of the original unrestricted process on the set LN , except for the equilibrium equation of state (N , 0, s − 1). Hence, the process restricted to LN is ergodic so that the series in (5.65) can be normalized to produce the equilibrium distribution of the restricted process on LN . Since the set N0 × Z × {0, 1, . . . , s − 1} \ LN is finite, it follows that the original process is ergodic and relating appropriately the equilibrium probabilities of the unrestricted and restricted process completes the proof. The following remark is in line with Remark 2. Remark 8 (An efficient numerical scheme) The following scheme exploits the fact that the rate of convergence of the series in (5.65) increases as m + |n| increases. So, further away from the origin, fewer compensation steps L are needed to achieve the same accuracy according to (4.1) and (4.2). This property was seen in Sect. 4 and in particular Fig. 7. Denote a triangular set of states as Tx := {(m, n) | m ∈ N0, n ∈ Z, m + |n| ≤ x }. Figure 10 serves as a visual aid. (i) Determine the minimal nonnegative integer N for which the spectral radii of the matrices R(1)(m, n) and R(2)(m, n) are both less than one for m + |n| > N and the spectral radius of R(3)(m) is less than one for m ≥ N . determine p(m, n) using the equilibrium equations m determine p(m, n) using the compensation approach m + |n| = M n m + |n| = K (ii) Select integers M and K such that N < M < K . (iii) Determine p(m, n), (m, n) ∈ TK \ TM according to (4.1) with C = 1 and L the minimal integer such that (4.2) holds. (iv) Determine p(m, n), (m, n) ∈ TM from the equilibrium equations for the states (m, n) ∈ TM . (v) Normalize the equilibrium distribution by dividing each equilibrium probability by the sum (m,n)∈TK p(m, n)1. The integer L in step (iii) depends on M . As M increases, L decreases or stays constant, but the size of the system of equilibrium equations in step (iv) increases. This tradeoff is clearly in favor of selecting a larger M : decreasing L decreases the number of computation steps exponentially, whereas the size of the system of equilibrium equations increases polynomially with M . As K increases, the equilibrium probabilities become more accurate. K can be chosen arbitrarily large; it has little to no impact on the performance. 6 Conclusion We have studied a queueing system with two nonidentical servers with service rates 1 and s ∈ N, respectively, Poisson arrivals and the SED routing policy. This policy assigns an arriving customer to the queue with the smallest expected delay, i.e., waiting time plus service time. The SED routing policy is a natural and simple routing policy that balances the load for the two nonidentical servers. Although not always optimal, SED performs well at both ends of system utilization range. Moreover, SED routing is asymptotically optimal in the heavy traffic regime. The SED system can be modeled as an inhomogeneous random walk in the quadrant. By appropriately transforming the state space, we mapped the two-dimensional state space into a half-plane with a finite third dimension. The random walks on each quadrant are different, yet homogeneous inside each quadrant. Extending the compensation approach to this three-dimensional setting, we showed in this paper that the equilibrium distribution of the joint queue length can be represented as a series of product-form solutions. These product-form solutions are generated iteratively to compensate for the error introduced by its preceding product-form term. The analysis presented in this paper proves that the compensation approach can be applied in the context of a three-dimensional state space. We believe that a similar analysis can be used to investigate general conditions for applicability of the compensation approach for three-dimensional Markov processes. These conditions will be comparable to the three conditions in Sect. 1. Furthermore, the compensation approach can possibly be extended to the case where all three dimensions are infinite, paving the way for performance analysis of higher-dimensional Markov processes. The insights gained for the SED system with two servers, specifically the series expressions for the equilibrium probabilities, can be used to develop approximations of the performance of heterogeneous multiserver systems with a SED routing protocol. These approximations can be derived along the same lines as in [20]. An approximate performance analysis for a system with two servers can be found in [27]. Another interesting direction for future research is to study rare events or tail probabilities in the SED system, in a similar way as done for JSQ systems in [26]. Since the compensation approach determines the complete expansion of each equilibrium probability, one can approximate rare events with arbitrary precision. Acknowledgments The authors would like to thank A. J. E. M. Janssen for proving step (c) of Lemma 4(i). This work was supported by an NWO free competition grant, an ERC starting grant, and the NWO Gravitation Project NETWORKS. Proof of step (c) of Lemma 4(i) The following argument was communicated to us by A. J. E. M. Janssen. Our goal is to show that for every α with |α| ∈ (0, 1) the equation has at least one root βi with |βi | < |α| for each i . To that end, we consider the more general formulation E σ = f (z) = zt (v+ − z)(z − v−), with z = β/α, E = (1+ss)ρ > 0, t = s +s1 ∈ (1, 2] and v± defined in (5.27) with the properties 0 < v− < 1, v+ > 1, v+ + v− = . One retrieves the original form (7.1) by setting σ = ui α1/s . The proof consists of three main steps. (a) We note that f (z) is analytic near z = v−. So, by the Lagrange inversion theorem, we have, in a neighborhood of σ = f (v−) = 0, z = g(σ ) = v− + = v− + ∞ σ n dn−1 z − v− n=1 n! dzn−1 zEt (v+ − z)(z − v−) ∞ (σ/E )n dn−1 znt dzn−1 (v+ − z)n z=v− . n Note that g(·) is analytic in a neighborhood of σ = f (v−) = 0. Let g1(z) = znt and g2(z) = (v+ − z)−n so that we can write where dn−1 dzn−1 g1(z)g2(z) = dn−k−1 dk dzn−k−1 g1(z) dzk g2(z), dn−k−1 dzn−k−1 g1(z) = n−k−2 i=0 (nt − i ) znt−(n−k−1), ddzkk g2(z) = (n(+n −k −1)1!)! (v+ − z)−(n+k). (7.2) < 1, (7.10) Fig. 11 The function σ = f (z) for z > 0 with ρ = 0.8 and s = 2 znt dzn−1 (v+ − z)n z=v− = n−1 (n + k − 1)! k=0 k!(n − 1 − k)! n−k−2 i=0 (nt − i ) vnt−(n−k−1) − (v+ − v−)n+k > 0. (7.8) Hence, the power series g(σ ) has positive power series coefficients. (b) We now extend the convergence range of g(σ ) from a neighborhood of 0 to the entire unit disk. This is achieved by using the Pringsheim theorem. We are allowed to do so because the coefficients of the power series of g(σ ) are positive. According to the Pringsheim theorem we can limit attention to the positive real line. With reference to Fig. 11, we let σmax = max{ f (z) | z ≥ v−}. Note that for z ∈ [v−, g(σmax)) the function f (z) is strictly increasing and is analytic by (7.2), and thus its inverse, g(σ ), is analytic for σ ∈ [ f (v−), f (g(σmax))) = [0, σmax). Clearly, σmax > 1, hence by the Pringsheim theorem g(σ ) is analytic for |σ | < σmax and thus inside the unit disk. (c) It is important to see where g(σmax) lies. We find z = g(σmax) from For t ∈ (1, 2) we take the positive solution of the quadratic equation in (7.9) and obtain after some rewriting d 0 = dz 1 zt (v+ − z)(z − v−) 1 = zt+1 v+v−t − (t − 1)(v+ + v−)z − (2 − t )z2 1 t 1 + ρ = zt+1 (1 + s)ρ − (t − 1) ρ z − (2 − t )z2 . z = g(σmax) = 1 + ρ + 2 (1 + ρ)2 + 4ρ(s − 1) which also holds for t = 2. So, for a σ with |σ | < σmax, we have from the positivity of the power series coefficients that |g(σ )| ≤ g(|σ |) ≤ g(σmax) < 1. By taking σi = ui α1/s , i = 1, 2, . . . , s with |σi | < 1 we see that the power series g(σi ) converges and |z| = |g(σi )| < 1. Finally, we can conclude that for every α with |α| ∈ (0, 1), (7.1) has at least one root βi with |βi | < |α| for each i . Proof of Lemma 10(iii) In this appendix we prove the convergence of the ratio of coefficients used in the horizontal compensation step for a solution that satisfies the positive inner equations. Specifically, we consider the setting of Lemma 7(i). We wish to establish the limiting values, or an upper bound on the absolute value of the following ratios, as α ↓ 0: (a) h/ν, (b) ηs+1/ν, (c) ηi /ν, i = 1, 2, . . . , s. We reiterate below the equations that determine these ratios, which can be found in the main body of the paper in (5.58)–(5.60), s i=1 A0,1 + α A−1,1 h − α ηi βi A1,−1 + α A0,−1 i+(α, βi ) + ηs+1βs+1 B1,1 + α B0,1 i−(α, βs+1) = −νβ A1,−1 + α A0,−1 i+(α, β), (8.2) − ηs+1αs + αsh(0) + (1 + s)ρqh(s − 1) = 0. As it turns out, when we let α ↓ 0 in the system of linear equations (8.1)–(8.3), we obtain a degenerate system of equations. By properly scaling the system, one obtains a nondegenerate system of equations. In this appendix we adopt the following notation to indicate the limiting value of a ratio: x lim x := lαi↓m0 ν . (a) Divide (8.3) by ν and rearrange to obtain h(s − 1) ν 1 = α q (1 + s)ρ ηs+1 h(0) ν − ν . (8.1) (8.3) (8.4) (8.5) By letting α ↓ 0 in (8.5) we get hlim(s − 1) = 0. Next, divide (8.2) by α and ν, which yields A1,−1 + α A0,−1 i+(α, β). (8.6) We let α ↓ 0 and use that A1,−1 = (1 + s)ρ I and B1,1 = (1 + s)ρ M (0,s−1). This gives B0,0hlim + v−(1 + s)ρ ηilime0 + ηsli+m1w−(1 + s)ρψ+(0)s−1e0 = −v+(1 + s)ρe0. (8.7) The matrix B0,0 is a tri-diagonal matrix and together with hlim(s − 1) = 0 we find that hlim = 0. (b) Together with hlim = 0, equation (8.7) reduces to the single equation ηilim + ηsli+m1w−ψ+(0)s−1 = −v+. Now, divide (8.1) by α and ν, yielding h h A0,1 αν + A−1,1 ν − i=1 ηi i+(α, βi ) = i+(α, β). Note that A0,1h/(αν) = (1 + s)ρ(1 − q)h(s − 1)/(αν)e0, so that, together with (8.5) and A−1,1 = I , (8.9) reduces to s 1 −q q ηsν+1 − h(ν0) i=1 ηi i+(α, βi ) = i+(α, β). ηilim + s 1 −q q ηsli+m1 = 1. One can solve the linear system of equations (8.8) and (8.11) for the two unknowns s i=1 ηilim and ηsli+m1, where the solution of ηsli+m1 is given in Lemma 10(iii)(b). (c) We wish to establish a scaling of each equation in (8.1)–(8.3) so that, in the limit for α ↓ 0, we obtain a nondegenerate system of equations. From (5.14) and (5.25) we know that i+(α, βi , r ) = ui βir/s , r = 0, 1, . . . , s − 1 and from Remark 4 it is clear α := 1 α−1/s α−2/s · · · α−(s−1)/s T . To be able to do that, we need some new notation. Let us introduce the column vectors η := η1 η2 · · · ηs T , and h˜ := h(0) h(1) α1/s α2/s · · · h(sα−1) T , and the matrices W (α) := i+(α, β1) i+(α, β2) · · · i+(α, βs ) , W˜ (α) := βα1 i+(α, β1) βα2 i+(α, β2) · · · βαs i+(α, βs ) . Using the introduced vectors and matrices we can write ηi η ν i+(α, βi ) = W (α) , and ν ηi βi η ν α i+(α, βi ) = W˜ (α) ν . i=1 We introduce the symbol ◦ for element-wise (Hadamard) multiplication, i.e., for column vectors x and y, we have that z = x ◦ y means that element z(r ) = x (r )y(r ). Below we list some useful element-wise multiplications that we will use: 1 β1 1/s u1 α 1 β2 1/s u2 α ∃ j : α ◦ i+(α, β) = 1 βα 1/s u j · · · βα (s−1)/s usj−1 T , α ◦ h = α1/s h˜ . Moreover, for α ↓ 0, we determine that α ◦ W (α) → W , α ◦ W˜ (α) → v−W , and α ◦ i+(α, β) → w j , where W and w j are given in (5.78)–(5.79). (8.17) (8.18) (8.19) We are now in a position to multiply both (8.6) and (8.9) element-wise by α. For the element-wise multiplication of (8.6) we get h α ◦ B0,0 ν + α ◦ + ηsν+1 βsα+1 α ◦ B1,1 + α B0,1 i−(α, βs+1) = − α α ◦ A1,−1 + α A0,−1 i+(α, β). β lim ηi α↓0 ν |ηljim(r )|. This concludes the proof of part (c). 1. Adan, I., Boxma, O., Resing, J.: Queueing models with multiple waiting lines. Queueing Syst. 37(1–3), 65–98 (2001) 2. Adan, I., Kapodistria, S., van Leeuwaarden, J.: Erlang arrivals joining the shorter queue. Queueing Syst. 74(2–3), 273–302 (2013) (8.20) (8.21) Lh˜ lim + v−W ηlim = −v+w j − γηs+1 w−ψ+(0)s−1e0, (1 + s)ρ(1 − q)M (0,s−1)h˜ lim − V ηlim = w j . We have constructed 2s equations, (8.23)–(8.24), for the 2s unknowns h˜ lim and ηlim. The solution to this system of equations depends on j . So, let h˜ lim and ηlim be the j j solution to (8.23)–(8.24) for a specific j . Then, 3. Adan, I., Wessels, J.: Shortest expected delay routing for Erlang servers. Queueing Syst. 23(1–4), 77–105 (1996) 4. Adan, I., Wessels, J., Zijm, W.: Analysis of the symmetric shortest queue problem. Stoch. Model. 6(4), 691–713 (1990) 5. Adan, I., Wessels, J., Zijm, W.: Analysis of the asymmetric shortest queue problem. Queueing Syst. 8(1), 1–58 (1991) 6. Adan, I., Wessels, J., Zijm, W.: A compensation approach for two-dimensional Markov processes. Adv. Appl. Probab. 25, 783–817 (1993) 7. Banawan, S., Zeidat, N.: A comparative study of load sharing in heterogeneous multicomputer systems. In: Proceedings of the 25th Annual Symposium on Simulation, pp. 22–31. IEEE Computer Society Press, Los Alamitos (1992) 8. Bateman, P., Diamond, H.: On the oscillation theorems of Pringsheim and Landau. In: Bambah, R.P., Dumir, V.C., Hans-Gill, R.J. (eds.) Number Theory, pp. 43–54. Springer, Basel (2000) 9. Blanc, J.: Bad luck when joining the shortest queue. Eur. J. Oper. Res. 195(1), 167–173 (2009) 10. de Bruijn, N.: Asymptotic Methods in Analysis, vol. 4. Courier Corporation, New York (1970) 11. Cohen, J.: Analysis of the asymmetrical shortest two-server queueing model. Int. J. Stoch. Anal. 11(2), 115–162 (1998) 12. Cohen, J., Boxma, O.: Boundary Value Problems in Queueing System Analysis. Elsevier, Amsterdam (2000) 13. Coombs-Reyes, J.: Customer allocation policies in a two server network: stability and exact asymptotics. Ph.D. thesis, Georgia Institute of Technology (2003) 14. Ephremides, A., Varaiya, P., Walrand, J.: A simple dynamic routing problem. IEEE Trans. Autom. Control 25(4), 690–693 (1980) 15. Fayolle, G., Iasnogorodski, R., Malyshev, V.: Random Walks in the Quarter-Plane: Algebraic Methods, Boundary Value Problems and Applications, vol. 40. Springer Science & Business Media, New York (1999) 16. Flatto, L., McKean, H.: Two queues in parallel. Commun. Pure Appl. Math. 30(2), 255–263 (1977) 17. Foley, R., McDonald, D.: Join the shortest queue: stability and exact asymptotics. Ann. Appl. Probab. 11, 569–607 (2001) 18. Foschini, G.: On heavy traffic diffusion analysis and dynamic routing in packet switched networks. In: Chadey, K., Reiser, M. (eds.) Computer Performance, pp. 499–513. NorthHolland, Amsterdam (1977) 19. Gould, H.: The Girard-Waring power sum formulas for symmetric functions, and Fibonacci sequences. Fibonacci Q 37, 135–140 (1999) 20. Gupta, V., Harchol-Balter, M., Sigman, K., Whitt, W.: Analysis of join-the-shortest-queue routing for web server farms. Perform. Eval. 64(9), 1062–1081 (2007) 21. Haight, F.: Two queues in parallel. Biometrika 45(3–4), 401–410 (1958) 22. Hajek, B.: Optimal control of two interacting service stations. IEEE Trans. Autom. Control 29(6), 491–499 (1984) 23. Halfin, S.: The shortest queue problem. J. Appl. Probab. 22, 865–878 (1985) 24. Hille, E.: Analytic Function Theory, vol. 1. Ginn, Boston (1959) 25. Kingman, J.: Two similar queues in parallel. Ann. Math. Stat. 32, 1314–1323 (1961) 26. Li, H., Miyazawa, M., Zhao, Y.: Geometric decay in a QBD process with countable background states with applications to a join-the-shortest-queue model. Stoch. Model. 23(3), 413–438 (2007) 27. Selen, J., Adan, I., Kapodistria, S.: Approximate performance analysis of generalized join the shortest queue routing. In: W. Knottenbelt, K. Wolter, A. Busic, M. Gribaudo, P. Reinecke (eds.) VALUETOOLS, 9th EAI International Conference on Performance Evaluation Methodologies and Tools (2016) 28. de Smit, J.: The queue G I /M/s with customers of different types or the queue G I /Hm /s. Adv. Appl. Probab. 15, 392–419 (1983) 29. Whitt, W.: Deciding which queue to join: some counterexamples. Oper. Res. 34(1), 55–62 (1986) 30. Yao, H., Knessl, C.: On the infinite server shortest queue problem: symmetric case. Stoch. Model. 21(1), 101–132 (2005) 31. Yao, H., Knessl, C.: On the shortest queue version of the Erlang loss model. Stud. Appl. Math. 120(2), 129–212 (2008)


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs11134-016-9497-7.pdf

Jori Selen, Ivo Adan, Stella Kapodistria, Johan van Leeuwaarden. Steady-state analysis of shortest expected delay routing, Queueing Systems, 2016, 309-354, DOI: 10.1007/s11134-016-9497-7