Steady-state analysis of shortest expected delay routing

We consider a queueing system consisting of two non-identical 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 non-identical 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 for numerical calculations and it also provides insight in the asymptotic behavior of the equilibrium probabilities as one of the state coordinates tends to infinity.


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 q 1 and q 2 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 q 1 + 1 and in queue 2 is (q 2 + 1)/s.The SED routing policy assigns an arriving customer to queue 1 if q 1 + 1 < (q 2 + 1)/s and to queue 2 if q 1 + 1 > (q 2 + 1)/s.In case the expected delays in both queues are equal, i.e. q 1 + 1 = (q 2 + 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 e.g.[12,Theorem 1], ρ := λ/(1 + s) < 1. (1.1) q 2 q 1 + 1 = (q 2 + 1)/s Transition rate diagram of the Markov process on the state space (q 1 , q 2 ) with s = 3 2 .Only transitions corresponding to arrivals that reach or cross the dashed line are shown.
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 (q 1 , q 2 ) with the transition rate diagram as in Figure 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(q 1 + 1) = q 2 + 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, there are many situations in which we have limited knowledge of the service times of the waiting customers.Such systems include the 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,10,17,18,21,24,26,27,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 [15,29].For the two non-identical server setting, Hajek [23] 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 [19] 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 models is evident from the analysis of the JSQ policy for two identical parallel exponential servers [22].The first major steps towards its analysis were made in [17,26], 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 readers is referred to [11,16].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,5,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 e.g.[1,3,4,5,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 Section 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 Section 3. Some numerical results are presented in Section 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 Section 6.

Variable Expression
Interpretation q 1 Number of customers (groups of size 1) in queue 1 q 2 Number of customers in queue 2 q 2 /s Number of groups of size s in queue 2 m m = min(q 1 , q 2 /s ) Minimum number of groups in queue 1 and 2 n n = q 2 /s − q 1 Difference between number of groups in queue 2 and 1 r r = mod(q 2 , s) Number of customers in queue 2 that are not in a group

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(q 1 + 1) = q 1 + 1, see Figure 1.In this section we transform the two-dimensional state space (q 1 , q 2 ) 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 2.1.
In queue 2 we count the number of groups of size s and denote it as j, i.e. j = q 2 /s , and we denote the number of remaining customers as r, i.e. r = mod(q 2 , 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 js + 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 q 1 + 1 < j + (r + 1)/s the arriving customer joins queue 1 and if q 1 + 1 > j + (r + 1)/s the arriving customer joins queue 2. In case the expected delays in both queues are equal, i.e. q 1 + 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 m = min(q 1 , j) and the difference between queue 2 and queue 1, i.e. n = j − q 1 .Using this notation, the SED system is formulated as a three-dimensional Markov process with state space {(m, n, r) | m ∈ N 0 , 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 with x T 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.
The transition rates are given by (m, n, r) Three-dimensional (m, n, r) transition rate diagram.Note that the third dimension is perpendicular and gives rise to s-layers.corresponding to arrivals, and corresponding to service completions.Figure 2(a) displays the transition rate diagram for the three-dimensional state space.The transition rates are described by the matrices A x,y in the positive quadrant and B x,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 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 The equilibrium equations can be written in matrix-vector form.We partition the state space as illustrated in Figure 2(b).For the interior of the positive and negative quadrant we have the following inner equations 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 The vertical boundary equations are Finally, for the three remaining boundary states near the origin, we have (A 0,0 + I)p(0, 1) + A 0,−1 p(0, 2) + A −1,1 p(1, 0) + A 0,1 p(0, 0) = 0, (2.9) (B 0,0 + sM (0,0) )p(0, −1) + B 0,1 p(0, −2) + B −1,−1 p(1, 0) + B 0,−1 p(0, 0) = 0, (2.10) (B 0,0 + I + sM (0,0) )p(0, 0) + A 0,−1 p(0, 1) + B 0,1 p(0, −1) = 0. (2.11) Remark 2.1.(Rational s) A system with a rational service rate s = s 2 s 1 can also be analyzed.In that case, one needs to consider a system with two servers and service rates s 1 and s 2 .Similar to our analysis at the start of Section 2, one denotes the number of groups of size s 1 in queue 1 as i and the number of groups of size s 2 in queue 2 as j.Then, let r n ∈ {0, 1, . . ., s n −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 (r 1 , r 2 ) ∈ {0, 1, . . ., s 1 − 1} × {0, 1, . . ., s 2 − 1}.This state space description leads to a transition rate diagram that has a similar structure as the one seen in Figure 2(a).In this sense, a system with a rational service rate can be analyzed using the approach described in this paper.

Evolution of the compensation approach and our contribution
In this section we use the 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 Section 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 (q 1 , q 2 ) ∈ N 2 0 , where q i is the number of customers at queue i, including a customer possibly in service.By defining m = min(q 1 , q 2 ) and n = q 2 − q 1 , 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 Figure 3(a).
For the symmetric JSQ model in [4], the initial solution is of the form η 0 α m 0 β n 0 , where η 0 is a coefficient and α 0 and β 0 are the compensation parameters that satisfy the kernel equation Figure 4: For a system with two identical servers and JSQ routing, the compensation approach generates in each compensation step a single compensation term.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 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 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: where C is the normalization constant.Figure 4 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 non-identical with speeds 1 and s.The symmetry argument used earlier does not hold anymore and one needs to consider the complete half-plane, see Figure 3(b) for the transition rate diagram.Note that the half-plane consists of two 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 Figure 5.
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 walk in the positive and negative quadrant are mirror images, which permits to perform the analysis only on the positive quadrant, see Figure 3(c) for the transition rate diagram.In [2] the authors extend the compensation approach to a three-dimensional setting.Due to the threedimensional state space, each compensation term takes the form α m β n i + (α, β), 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 Figure 6.
Our contribution.The model at hand is defined on an s-layered half plane, thus requires Figure 5: For the asymmetric JSQ model, the compensation approach generates different compensation terms for the positive (straight arrow) and the negative quadrant (snaked arrow).

• • • initial solution vertical compensation horizontal compensation
Figure 6: For the symmetric JSQ model with Erlang-s arrivals, the compensation approach generates s compensation terms in a HCS and just one compensation term in a VCS. 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 Figure 7. Additionally, the product-form solutions take the form α m β n i + (α, β) or α m β −n i − (α, β) 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, 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 has 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 5.11, the horizontal compensation step is described in Lemma 5.13, and the vertical compensation step is described in Lemma 5.12, see Figure 7.

Numerical results
Expression (3.3) is amenable for numerical calculations after applying truncation.For m ≥ 0, n ≥ 1, (c) ρ = 0.9 Figure 8: Heat plot of the equilibrium distribution mass for the SED system with s = 3, q = 0.4, and varying ρ, determined according to (4.1) with L = 16.The heat plot shows where the probability mass is located (darker colour means more mass).The dashed lines are q 1 + 1 = (q 2 + 1)/s and q 1 = q 2 /s.
where the empty sum −1 l=0 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 q 1 + 1 = (q 2 + 1)/s and the line where the expected waiting time is equal q 1 = q 2 /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 equilibrium distribution for this model and varying ρ is given in Figure 8 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 Figure 9 shows that away from the origin, the convergence of the series is very fast, but the convergence is also quite fast for states close to the origin.Note that the distance of a state (m, n) to the origin is directly related to the rate of convergence of the series expression (3.3).In particular, it seems to be a function of m + |n|: faster convergence further away from the origin.This property is formally proven in Section 5.7 and can be exploited for numerical computations.
Remark 4.1.(No curse of dimensionality) Take a triangular set of states where M is some non-negative integer.Technically, M needs to be strictly larger than some non-negative integer N , but we do not go into the details here; the lower bound N is described in Section 5.7.For states outside T M , 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) / ∈ T M is relatively small.For states (m, n) ∈ T M , 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) / ∈ T M are known.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) / ∈ T M , see Figure 9.So it is evident that the compensation approach does not suffer from the curse of dimensionality.

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 Section 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 Section 5.3.The vertical compensation step is outlined in Section 5.4.Section 5.5 describes the horizontal compensation procedure.We formalize the resulting solution in terms of a sequence of product-forms in Section 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.

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 productform 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, Section 4].
The related model is constructed as follows.We start from the state space of the original model in Figure 2(a) and bend the vertical axis as shown in Figure 10.Note that for this modified model, m ∈ Z.We add an additional state (−1, 0, s − 1) with a transition rate Â0,1 = se 0 to state (−1, 1) and a transition rate B0,−1 = (1 + s)ρqe s−1 to state (−1, −1), where e i 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, Â0,−1 = se T 0 .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, B0,1 = 1 T , 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.
In Lemma 5.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 5.2.For ρ < 1, the equilibrium distribution of the modified model is of the form (5.6) Proof.Let k denote the total number of customers in the system, i.e.
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 pk denote the probability of having k customers in the system.From the balance principle between states k and k + 1 we obtain pk+1 = ρp k , k > s. (5.8) (5.9) Combining the last two results immediately yields α = ρ 1+s .
Combining Lemmas 5.1 and 5.2 yields the following result for the original model.
with q(n) = (q(n, 0), q(n, 1), . . ., q(n, s − 1)) T non-zero and ∞ n=−∞ ρ −(1+s)|n| q(n, r) < ∞, r = 0, 1, . . ., s − 1. (5.11) The solution obtained in Proposition 5.3 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 5.4, 5.6 and 5.9, 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.
Lemma 5.6.The equation det(D + (α, β)) = 0 assumes the form and can be rewritten to where u i is the i-th root of unity of u s = 1 and β 1/s is the principal root.
Theorem 5.8.(de Smit [14]) Let A(z) = (a i,j (z)) and B(z) = (b i,j (z)) be complex n × n matrices, where B(z) is diagonal.The elements a i,j (z) and b i,j (z), 0 ≤ i, j ≤ n − 1 are meromorphic functions in a simply connected region S in which T is the set of all poles of these functions.C is a rectifiable closed Jordan curve in S \ T .Let x B and x A+B be the number of zeros inside C of det(B(z)) and det(A(z) + B(z)), respectively, and y B and y A+B the number of poles inside C of det(B(z)) and det(A(z) + B(z)) (zeros and poles of higher order are counted according to this order).

Compensation on the vertical boundary
The initial solution (5.43) does not satisfy the positive and negative vertical boundary.We consider a term ηα m β n i + (α, β) 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 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 s i=1 ν i α m i β n i + (α i , β) to the original term, such that the linear combination . 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 5.7 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 non-zero.
Let us now consider a term from the negative quadrant ηα m β −n i − (α, β) 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 case of the initial solution (5.43), this would be the solution on the negative quadrant.As before, the compensation approach dictates to add a compensation term ν s+1 α m s+1 β −n i − (α s+1 , β) to the original term, such that the linear combination where 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 5.12.(Vertical compensation) (i) Consider the product form ηα m β n i + (α, β) 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 + (α, Then there exists a coefficient ν 1 such that 2) and (2.7).The coefficient ν 1 satisfies (5.52) (ii) Consider a product form ηα m β n i − (α, β) 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 non-zero vector satisfying (5.17).Then there exists a coefficient ν s+1 such that 3) and (2.8).The coefficient ν s+1 is given by . (5.54) Proof.(i) Use (5.12) to establish that (5.55) Then using Remark 5.7 we can establish that there exists a single ν i that is non-zero.We label the single non-zero coefficient as ν 1 and find it from (5.48) with the help of (5.55).

Compensation on the horizontal boundary
The solution obtained after compensation on the vertical boundary, as outlined in Lemma 5.12, 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 quadrant.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.
(ii) The proof is identical to the proof of (i).

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.
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: (I) The initial solution is determined from Lemma 5.11.
(V) For a vertical compensation step with fixed β l,i : (i) If the index i is not a multiple of s + 1, then the compensation terms are determined according to Lemma 5.12(i).
(ii) Otherwise, the compensation terms are determined according to Lemma 5.12(ii).
(H) For a horizontal compensation step with fixed α l,i : (i) If the index i is not a multiple of s + 1, then the compensation terms are determined according to Lemma 5.13(i).
(ii) Otherwise, the compensation terms are determined according to Lemma 5.13(ii).

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∈N 0 , i=1,2,...,(s+1) l and {β l,i } l∈N 0 , 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.
Proof.(i) Follows immediately from Lemmas 5.6 and 5.9.(ii) For a fixed α, let β 1 , β 2 , . . ., β s be the roots of (5.24)  and let C = {α ∈ C | |α| ≤ α 0,1 = ρ 1+s }.Using Lemmas 5.6 and 5.9 we have t(α) < 1 and in particular, since ρ ∈ (0, 1), we have that C is a closed proper subset of the unit disc and thus c 1 := max α∈C t(α) < 1.One can reason along the same lines and obtain a bound c 2 for a fixed β.The exponential decrease in terms of the level l follows from the fact that each β that is generated from an α satisfies |β| < |α|c 1 and each α that is generated from a β satisfies |α| < |β|c 2 .Thus, we define c := c 1 c 2 .
(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 5.16.
(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 B. (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 5.15, 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.Proof.(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 it's descendants as R
Remark 5.19.The index N is the minimal non-negative 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) 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.10.Then L N is a set of states for which the series in (5.65) converge absolutely.The restricted stochastic process on the set L N is an irreducible Markov process, whose associated equilibrium equations are identical to the equations of the original unrestricted process on the set L N , expect for the equilibrium equation of state (N, 0, s − 1).Hence, the process restricted to L N is ergodic so that the series in (5.65) can be normalized to produce the equilibrium distribution of the restricted process on L N .Since the set N 0 × Z × {0, 1, . . ., s − 1} \ L N 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 4.1.
Remark 5.21.(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 Section 4 and in particular Figure 9. Denote a triangular set of states as T (i) Determine the minimal non-negative 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 .
(ii) Select integers M and K such that N < M < K.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.

Conclusion
We have studied a queueing system with two non-identical 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 non-identical 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, 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 threedimensional 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 threedimensional Markov processes.These conditions will be comparable to the three conditions in Section 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 multi-server systems with a SED routing protocol.These approximations can be derived along the same lines as in [21].An approximate performance analysis for a system with two servers can be found in [28].
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 [27].Since the compensation approach determines the complete expansion of each equilibrium probability, one can approximate rare events with arbitrary precision.
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),A 0,1 + αA −1,1 h − α As it turns out, when we let α ↓ 0 in the system of linear equations (B.1)-(B.3),we obtain a degenerate system of equations.By properly scaling the system, one obtains a non-degenerate system of equations.
In this appendix we adopt the following notation to indicate the limiting value of a ratio x lim := lim This concludes the proof of part (c).

Figure 1 :
Figure1: Transition rate diagram of the Markov process on the state space (q 1 , q 2 ) with s =3  2 .Only transitions corresponding to arrivals that reach or cross the dashed line are shown.

Figure 2 :
Figure 2: Transition rate diagram and state space partitioning.
Identical servers and Erlang arrivals with s phases (s-layered transition rate diagram)[2]

Figure 3 :
Figure 3: Simplified transition rate diagrams on the state space (m, n) for JSQ systems with two servers.

Figure 7 :
Figure7: For the SED system, the compensation approach generates different compensation terms for the positive (straight arrow) or the negative quadrant (snaked arrow).The number of terms generated in the positive and negative quadrant are different.

Figure 9 :
Figure 9: The number of compensation steps L for each state (m, n) such that the resulting p L (m, n) is accurate according to (4.2).The dashed line is m + |n| = 3. Parameters are s = 4, ρ = 0.8, and q = 0.4.

Lemma 5 . 1 .
For ρ < 1, the equilibrium distribution of the modified model exists and is of the form

( i )
For every α with |α| ∈ (0, 1), equation (5.25) has exactly one root β i inside the open circle of radius |α| for each i.Furthermore, all s roots β i are distinct.(ii) For every β with |β| ∈ (0, 1), equation (5.25) has exactly one root α i inside the open circle of radius |β| for each i.Furthermore, all s roots α i are distinct.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: (a) For every fixed α with |α| ∈ (0, 1), equation (5.24) has exactly s roots β inside the open circle of radius |α|.(b) These s roots β are distinct.(c) For every fixed α with |α| ∈ (0, 1), equation (5.25) has at least one root β i inside the open circle of radius |α| for each i.Combining these three steps proves that for every fixed α with |α| ∈ (0, 1), equation (5.25) has at exactly one root β i inside the open circle of radius |α| for each i and the β i 's are distinct.(i)(a) Equation (5.24) is a polynomial of degree 2s in β and we will show that exactly s roots are inside the open circle of radius |α|.Other possible roots inside the open unit circle appear not to be useful, since these will produce a divergent solution, as follows from (5.11).Divide both sides of (5.24) by α 2s and set z = β/α to obtain f (z) s − αs s z s+1 = 0 (5.26) with f (z) = (1 + s)(ρ + 1)z − (1 + s)ρz 2 − 1.Then f (z) has the two roots

Figure 11 :
Figure 11: Indexing of the compensation terms for a subtree.A dashed rectangle indicates according to which lemma the compensation terms are generated.A straight arrow indicates that the roots α or β are obtained through (5.24) and a snaked arrow through (5.35).

Theorem 5 . 20 .
For all states (m, n, r), m ∈ N 0 , 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 non-negative 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, Section 11], but we include it here for completeness.Define L N = {(m, n) | m ∈ N 0 , n ∈ Z, m + |n| > N } ∪ {(N, 0, s − 1)} and note the similarity with the set of states defined in the proof of Lemma 5.1 and the associated transition rate diagram in Figure

Figure 12
serves as a visual aid.

m
+ |n| = M determine p(m, n) using the equilibrium equations m + |n| = K determine p(m, n) using the compensation approach

Figure 12 :
Figure 12: Regions of the efficient numerical scheme.

Table 1 :
Interpretation of the state variables of the Markov process.