Analysis of the shortest relay queue policy in a cooperative random access network with collisions

The scope of this work is twofold: On the one hand, strongly motivated by emerging engineering issues in multiple access communication systems, we investigate the performance of a slotted-time relay-assisted cooperative random access wireless network with collisions and with join the shortest queue relay-routing protocol. For this model, we investigate the stability condition, and apply different methods to derive the joint equilibrium distribution of the queue lengths. On the other hand, using the cooperative communication system as a vehicle for illustration, we investigate and compare three different approaches for this type of multi-dimensional stochastic processes, namely the compensation approach, the power series algorithm (PSA), and the probability generating function (PGF) approach. We present an extensive numerical comparison of the compensation approach and PSA, and discuss which method performs better in terms of accuracy and computation time. We also provide details on how to compute the PGF in terms of a solution of a Riemann-Hilbert boundary value problem.

A typical relay-based cooperative wireless network operates as follows: There exists a network of a finite number of source users, a finite number of relay nodes and a common destination node. The source users transmit packets to the destination node with the cooperation of the relays. If a direct transmission of a user's packet to the destination fails, a cooperation strategy among sources and relays is employed to specify the relays that will store the blocked packet in their buffer. Relays are responsible for the transmission of the blocked packets to the destination, e.g., [46]. In a wireless network, transmission failures occur either due to packet collisions, or due to channel fading/noise and attenuation. In both cases the packet has to be re-transmitted at a later slot. The former case occurs when more than one node transmit simultaneously, while the latter one refers to the probabilistic nature of transmissions, see e.g., [46,50,51].
In this work, we consider a simple relay-based cooperative wireless network with a single source, two infinite capacity relay nodes, and a common destination with collisions under a load balancing relay scheme. We assume that due to deep fading and bad channel quality it is impossible for the source to communicate with the destination through a direct link, and thus, cooperation within the relays is imperative. Furthermore, we assume that the relays and the source user are sufficiently close such that the packets transmitted over the channel are always correctly received by the relays. The employed cooperation strategy among source and relays is queue-based, with as ultimate goal to minimise the total transmission time, i.e., the time that is needed to transmit a packet from the source to the destination. To this end, we consider join the shortest queue policy as it seems to be the most appropriate for such a wireless network [39,40,60,63].

Related work
For networks with cooperation the first point of interest concerns the characterisation of the stable throughput region, aka the stability region. For small, simple networks the stability region can be fully determined, see, e.g., [17,43,58], while for large, general networks, only bounds of these regions are known [12,38,49,55]. For a thorough overview of several techniques used for the derivation of the stability region the interested reader is referred to [24]. An alternative way to derive stability conditions is to use the concept of dominant systems, under which the network of interest is (stochastically) compared to a simpler one that is easier to analyse, see, e.g., [48] and the references therein. Another powerful tool to investigate necessary and sufficient stability conditions for work-conserving queueing networks is the use of fluid models, see [11,25,28,52,53].
Next to the characterisation of the stability region, the delay performance (i.e., the investigation of the joint relay queue length distribution) is another crucial performance measure in random access networks, which recently has regained attention due to the rapid development of real-time applications that require delay-based guarantees [27]. However, the impact of interacting queues causes severe mathematical difficulties, and there are very few studies that deal with the stability of the queueing delay. In [42], by performing an appropriate truncation of the infinite Markov chain, the authors approximated the steady-state probability vector, within any desired degree of accuracy, whereas in [56], diffusion approximations were applied. In [17,41,47] the PGF of the joint (relay) queue length distribution was obtained in terms of the solution of a boundary value problems. In [13,62] fluid models were used to investigate the delay analysis of random access networks. In [54] bounds for the queueing delay in a random access network with N > 2 nodes were also derived.
Note that, in the vast majority of the above mentioned literature, each user node in the network has its own dedicated traffic. Alternatively, the inclusion of cooperation under certain criteria gives rise to the introduction of load balancing techniques that forward the packets to specific relays with ultimate goal the optimisation of the overall network performance; e.g., [37].
Load balancing schemes are used to improve scalability and overall system throughput in dis-tributed systems. Such schemes improve the system's performance by dividing the work effectively among the participating nodes. Under certain conditions (identical servers, exponential service times, and service protocols that do not depend on the service requirement, such as the FCFS), the so-called join the shortest queue (JSQ) policy has several strong optimality properties: The JSQ policy minimises the overall mean delay among the class of load balancing policies that do not have any advance knowledge of the service requirements [19,40,60].
For small scale systems with two queues, the JSQ policy has been extensively studied. The JSQ policy was introduced in [30], in the context of two parallel (exponential) servers and a Poisson stream of arrivals. The first major steps towards its exact analysis were made in [22,33], using a uniformization approach that established that the PGF of the joint equilibrium distribution of the queue lengths is meromorphic (i.e., the equilibrium distribution can be written as a linear combination -finite or infinite -of product-form terms). For an extensive treatment of the JSQ system using a generating function approach, the interested reader is referred to [15,21]. An alternative approach that is not based on generating functions is the compensation approach, that directly solves the balance equations and obtains the equilibrium distribution in the form of a linear combination -finite or infinite -of product-form terms, see [3,4]. It is worth noting, that although there exist multiple analyses for the JSQ system, this line of work is very challenging, see [1] for a comparison of some of the approaches and an exposition of their strengths and limitations. Using the exact analysis for the case of two servers, Gupta et al. [29] obtained the first approximate analysis of the marginal equilibrium distribution of the queue lengths in a server farm with JSQ routing protocol, processor sharing service discipline and generally distributed service times.
A mainly numerical approach, that has been successfully used for among others JSQ systems, is the power-series algorithm (PSA), [8,10], which assumes that the equilibrium distribution can be written as a function of a system parameter (such as the load). In a different direction, the precedence relation (PR) method, see, e.g., [61], can be used to systematically construct bounds for any performance measure of the Markov chain under consideration. The strength of the PSA and the PR methods is that they are not restricted to two queues, but apply equally well to more than two queues. A weak point, however, is that the theoretical foundation of PSA is still incomplete and that PR produces no exact results, but bounds.
For large scale systems, the exact analysis is elusive and the vast majority of the literature focuses on asymptotic regimes, and heavy-traffic delay optimality in steady state. The most common asymptotic regimes include fluid limits, heavy traffic limits, mean field limits and the Halfin-Whitt (qualityefficiency-driven) regime. Under the assumption of a Poisson arrival stream with rate λ, N identical (exponential) servers, each with service rate µ, the above regimes can be viewed as different limits of the quantity λ − N µ, after appropriately rescaling by some function of N . See the seminal paper of Foschini and Salz [23] for the case N = 2, and [6] for the mean field limit analysis of a large underloaded parallel server model.

Fundamental contribution.
In this work, we extend for the first time the framework of the compensation approach to a new class of random walks, viz., that violates one of the main assumptions of the compensation approach that of transitions being allowed only to neighbouring states. We show, using the system at hand as a vehicle for illustration, that the compensation approach can be extended to random walks in the quadrant that obey the following conditions: • 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.
• Forbidden steps: No transitions from interior states to the North, North-East, and East.
• Bounded transitions: Only transitions to a bounded region.
Note that the compensation approach was developed to satisfy a more restrictive setting than the bounded transitions, by allowing only transition to the nearest neighbours, see, e.g., [3,4,5]. This is a very promising result, since it seems to be possible to extend the compensation approach to random walks with large steps (extending the results obtained in [2]) and to queueing systems with (bounded) batch arrivals and/or (bounded) batch departures.

Paper structure
The paper is structured as follows: In Section 2, we describe the system under consideration in detail, and provide its stability condition. The computation of the equilibrium distribution of the joint queue lengths using the compensation approach is performed in Section 3. In Section 4, we apply PSA, while in Section 5, we provide a detailed analysis on how to derive the probability generating function of the equilibrium joint queue length distribution in terms of a solution of a Riemann-Hilbert boundary value problem. Section 6 is devoted to the comparison of the compensation approach and PSA method. Finally, we present conclusions and possible generalisations in Section 7.

The model
We consider a relay-assisted cooperative random access wireless network composed of a saturated source user, that transmits packets to a common destination node, under the cooperation of two relay nodes. The relays are equipped with infinite capacity buffers (queues), and they assist the user by transmitting the packets that failed to reach the destination. Packets have equal length (i.e., they consist of the same fixed number of bits [14]), and time is divided into slots corresponding to the transmission time of a packet. We consider an Early Arrival System (EAS), under which at the beginning of a slot packets arrive and they are routed to the relays according to the join the shortest relay queue (JSRQ) policy. On the other hand, departures are scheduled at the end of the slot.
Packet arrivals are assumed i.i.d. Bernoulli random variables from slot to slot, with the average number of arrivals being λ packets per slot. Upon the arrival of a packet, the source and the relays cooperate as follows: When the source transmits a packet, it is forwarded to the least loaded relay, i.e., to the relay with the smallest number of backlogged packets. Then, the relay node sends an acknowledgement to the source and takes over the responsibility of delivering the packet to the destination node by storing it in its queue. Such a protocol helps to keep a fair balance among the relays, as well as, it enhances the energy conservation of the relays (the relay node is usually a battery operated wireless device). In case the numbers of packets in the relay queues are equal, a packet is routed to relay r with probability π r , r = 1, 2. At the end of each slot, relay r (if it is non-empty) transmits a packet to the destination node with probability a r , r = 1, 2. If both relays transmit at the same slot, a collision occurs, and both packets have to be retransmitted in a later slot. If only one relay transmits, then the destination node successfully decodes it, sends an acknowledgement to the corresponding relay, and the packet exits the network. The acknowledgements are assumed to be error-free, instantaneous, and broadcasted to all relevant nodes. The nodes remove the successfully transmitted packets from their queues, while unsuccessful packets are retained.
In order to enhance the readability of the paper, and only for this reason, we focus hereon at the symmetric system, under which a r = a, r = 1, 2, i.e., both relays have identical transmission parameters, and π r = 1/2, r = 1, 2. Note that the analysis that follows can be directly generalised to the asymmetric system case, however this would render the notation more complicated, which would severely impact the readability of the paper.
Let Q r (n) be the number of stored packets at the buffer of relay r, r = 1, 2, at the beginning of the n-th slot, n ≥ 0. Then {Q(n), n ≥ 0} := {(Q 1 (n), Q 2 (n)), n ≥ 0} is a discrete time Markov chain with state space S = {(i, j) : i, j ≥ 0}. The corresponding probability transition diagram is depicted in Figure 1.
Remark 2.1. To better understand how the above probabilities are computed, we consider for example the case of p V 1,0 . This probability captures the transition from a state (i, j) ∈ V to a state (i + 1, j). In this case, at the beginning of a slot, due to the EAS system, with probability λ there is a new arrival. For the transition (i, j) → (i + 1, j) to occur it is necessary on top of the arrival that no departure occurs. The latter event happens with probabilityā 2 + a 2 . Thus, p V 1,0 = λ(ā 2 + a 2 ). The other transition probabilities are obtained in a similar fashion.

Remark 2.2.
Note that if both relay queues are non empty, then the successful transmission rate from each of them equalsāa. If one of them is empty then the other transmits with rate a. This demonstrates that the setting under consideration is a non-work-conservative setting. The only exception is the case 2āa = a ⇐⇒ a = 1/2 =ā (or equivalently the case 2āa =ā). Due to the non-work conservation setting, our model incorporates two features, that of the JSRQ and that of the "coupled processors" [20]. The combination of the JSRQ feature and the coupled processor feature considerably complicates the analysis.

Stability condition
Let (E k x , E k y ) denote the mean jump vector of {Q(n), n ≥ 0} in the angles k = H, V or in the rays k = H , V , D. Then, it is readily derived that Following the analysis in [35], we prove the following theorem.
Equivalently, the stability condition for the system can be written in terms of the system load Proof. The proof consists of two parts: in the first part, we show that the stated condition is sufficient and in the second part, we show that the stated condition is necessary.

Sufficiency.
To this purpose, we use Foster's criterion [21]: A Markov chain is ergodic provided that there exists a positive function f (x, y) on Z 2 + , a number > 0, and a finite set A ∈ Z 2 + such that with (θ x , θ y ) a random vector distributed as the one step jump of the Markov chain from state (x, y).
Assume that (5) holds and consider the function f (x, y) = x 2 + y 2 , which satisfies, and for some 0 > 0, and for all (x, y) : Assume now that x > 0, y = 0. Then since E H x < 0, for some 1 > 0, and for all x sufficiently large.
The case y > x is symmetric to x > y and further details are omitted. For x = y > 0, since , we can check (7) using again (8). Then, Foster's criterion applies and the chain is ergodic.
Necessity. It is sufficient to show that there exists a function f (x, y) on Z 2 + and a constant c > 0 such that Similarly, if x = 0, y > 0, Thus, f (x, y) satisfies (9) and the chain is non-ergodic.
Remark 2.4. Theorem 2.3 can be generalised to the asymmetric case following [35]. This analysis would reveal that the stability condition is

Equilibrium analysis: compensation approach
The compensation approach is developed by Adan et al. in a series of papers [3,4,5] and aims at a direct solution for the equilibrium joint queue length distribution, for a sub-class of two-dimensional random walks on the lattice of the first quadrant obeying the following conditions: • 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.
• Forbidden steps: No transitions from interior states to the North, North-East, and East. • Step size: Only transitions to neighbouring states.
It exploits the fact that the balance equations in the interior of the quarter plane are satisfied by a linear combination (finite or infinite) of product-form terms, the parameters of which satisfy a kernel equation, and that need to be chosen such that the balance 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.
As evident from Figure 1, the model at hand violates the first two conditions mentioned above. In order to apply the compensation approach, it is necessary to transform the state space (similarly as in the case of the classical JSQ model [3,4]). More concretely, we employ the following transformatioñ The corresponding probability transition diagram is depicted in Figure 2. Note that the above transformation has led to a random walk that only violates the nearest neighbour condition for the application of the compensation approach, but the new random walk has bounded transitions, see Figure 2. Despite violating the nearest neighbour condition, we show, in this paper, that the compensation approach can still be applied and lead to an equilibrium joint queue length distribution expressed in the form of a linear combination (infinite) of product-form terms. Let denote the equilibrium joint queue length distribution of the transformed random walk {Q(n), n ≥ 0}.

Decay rate
From the exact analysis of the system, see Section 3.2 and the result therein, it turns out that, although the model at hand combines both the JSRQ feature and the coupled processor feature, the dominating feature is that of the JSRQ. This is already noticeable in the decay result of the following proposition, indicating a behaviour similar to that of the classical JSQ, cf. [3].
with c, d 0 , δ constants that do not depend on k.
Proposition 3.1 can be intuitively understood by comparing the transformed model with the corresponding Geo/Geo/2 model with one queue, Bernoulli arrivals with success probability λ, two identical servers, each with Bernoulli service with success probabilityāa. Let Q 1 and Q 2 denote the equilibrium queue lengths in the original random walk, and let Q denote the equilibrium queue length in the Geo/Geo/2 model with one queue. It is then expected that P(Q 1 + Q 2 = k) and P(Q = k) have the same decay rate, since both systems will work at full capacity whenever the total number of customers grows large. Moreover, since the JSRQ protocol constantly aims at balancing the lengths of the two queues over time, it is expected that, for large values of k, For the standard Geo/Geo/2 model with one queue, it is well known that Combining (19) and (20) leads to the following conjectured behaviour of the tail probability of the minimum queue length for some positive constant C. Hence, the decay rate of the tail probabilities for min{Q 1 , Q 2 } is conjectured to be equal to the square of the decay rate of the tail probabilities of Q. Proposition 3.1 states this conjecture, for the case when the difference of the queue sizes is fixed.
Furthermore, one can determine the decay rate of the marginal distribution of min(Q 1 , Q 2 ). The latter does not follow immediately from Proposition 3.1, because the summation over the difference in queue sizes, which can be unbounded, requires a formal justification. In this paper, we derive an exact expression for π k,l , which, among other things renders the following result.
with c, d 0 , δ constants that do not depend on k.
The proofs of Propositions 3.1 and 3.2, along with several other asymptotic and exact results, are given in Section 3.2.

The compensation procedure
In this section, we obtain the equilibrium joint queue length distribution using the compensation approach. This approach yields an explicit expression by directly exploiting the balance equations, without any transforms.
The compensation approach as has been briefly discussed in the introduction attempts to solve the balance equations by a linear combination of product-form terms. This is achieved by first characterising a sufficiently rich basis of product-form solutions satisfying the balance equations in the interior of the state space. Subsequently this basis is used to construct a linear combination that also satisfies the equations for the boundary states. Note that the basis contains uncountably many elements. Therefore a procedure is needed to select the appropriate elements. This procedure is based on a compensation argument (which explains the name of the method): after introducing the first term, countably many terms may subsequently be added so as to alternatingly compensate for the error on one of the two boundaries. The main steps of the analysis are briefly outlined below.
The solutions to Equation (23) form the basis. In particular the points on the curve (23) and inside the region 0 < |γ|, |δ| < 1 characterise a continuum of product-forms satisfying the inner equations.
Step 2 Construct a linear combination of elements in this rich basis, which is a formal solution to the balance equations. Here the word formal is used to indicate that (at this stage) we do not treat the convergence of the solution. This aspect is treated in Step 3. The formal solution is constructed as follows: (a) The construction of a linear combination starts with a suitable initial term γ 0 that satisfies the interior of the state space and also the balance equations (11)- (13). In Lemma B.1, in Appendix B, we prove that γ 0 = ρ 2 . Then, from Equation (23) we obtain the unique δ 0 , with |δ 0 | < |γ 0 |, such that satisfies Equations (10)- (13), where the double turnstile symbol ( ) is used to signify that π k,l semantically entails the form d 0 γ k 0 δ l 0 . The uniqueness of the δ-root is proven in Lemma B.2. The constant d 0 can be set equal to one, this will be corrected in Step 4 with the computation of the normalisation constant.
(b) The starting tuple (γ 0 , δ 0 ) violates Equation (14) on the vertical boundary. To compensate for this error, we add a new product-form term coming from the basis, such that the sum of the two terms satisfies the balance equations in all states on the vertical boundary. In particular, it is easy to show that the new tuple is (γ 1 , δ 0 ). The new γ-root is uniquely determined from Equation (23), with |δ 1 | < |γ 0 |. The uniqueness of the γ-root is proven in Lemma B.2. Then, satisfies (14) if with d 0 known constant from the previous step.
(c) Adding the new term violates the balance equations (11)-(13), hence we compensate for this error by adding a product-form solution (γ 1 , δ 1 ) satisfying (23), and (11)- (13), such that The three unknowns e 0,1 , e 1,1 and d 1 can be computed from the following system of three equations with (d) We continue in this manner until we construct the entire formal series and π 0,0 , π 0,1 , and π 0,2 are obtained from Equations (15)- (17). Note that Equation (26) can be equivalently written as follows Step 3 Prove that the formal solution (26) and (27) converges. This is split up into two parts: i) we first show in Proposition B.3 that the sequences {γ i } i∈N and {δ i } i∈N converge to zero exponentially fast, and ii) we show in Proposition B.6 that the formal solution converges absolutely in all states. The Propositions and their proofs are in Appendix B.
Step 4 Determine the normalisation constant.
Performing the steps described above for the compensation approach leads to the following main result for the equilibrium distribution.
with c denoting the normalisation constant. The sequences {γ i } i∈N , {δ i } i∈N , {c i } i∈N , {d i } i∈N , and {e i } i∈N are obtained recursively based on the analysis of Steps 1-4 above.
Clearly, from this result we can derive similar expressions for other performance characteristics such as the mean queue lengths, the correlation between the queue lengths, the mean waiting time, etc.

Numerical implementation of the compensation approach
In this section, we discuss how to numerically implement the equilibrium distribution analysis based on the compensation approach. The equilibrium distribution π k,l , k ≥ 0, l ≥ 0, is written as a linear combination of product-form terms, cf. (28), (29). The first step of the numerical implementation is to consider the first few terms of the series expression for π k,l . More concretely, let where here N ca denotes the truncation level of the series expression obtained using the compensation approach. From the analysis of Proposition B.6, it follows that the inclusion of more terms of the series expression improves to a desired accuracy the computation of the equilibrium distribution, i.e., the bigger the value of N ca the better the approximation of π k,l . The constant N ca is determined such that Nca k,l=0 with ε the precision error. This is included as the stopping criterium for the algorithm, i.e., we start with N ca = 1 and as long as Equation (32) is not satisfied, we increase the value of N ca by one.
In Algorithm 1, we provide all the necessary information for the numerical implementation of the compensation approach.

Applicability of the compensation approach in case of bounded transitions
In Section 3, we mentioned that the model at hand violates the nearest neighbour condition for the applicability of the compensation approach. Nonetheless, the analysis performed demonstrated that the compensation approach can be generalised to cover a larger class of random walks permitting transitions not only to the nearest neighbours, but to a bounded region of neighbouring states. From our analysis, it becomes clear that for random walks with the structure of the system at hand and bounded quasi-birth-death transitions along the rays R L = {(k, l) : 2k + l = L}, L ≥ L 0 , the analysis performed in Section 3.2 carries out to a tee. We have confirmed this in the system depicted in Figure 3.

Equilibrium analysis: Power series algorithm implementation
In this section, we show how the power series algorithm (PSA) can be used to analyse the relay-based cooperative communication network with collisions and with JSRQ protocol. PSA is an algorithmic procedure which is often used to numerically obtain the performance measures of multi-dimensional queueing models, which fit in the class of quasi birth-and-death processes. The intrinsic idea behind PSA is the transformation of the non-recursively solvable set of balance equations into a recursively solvable set of equations by adding one dimension into the state space. This is achieved by expressing the equilibrium distribution as power series in some variable based on the model parameters, this allows the calculation of the equilibrium joint queue length distribution. PSA was first applied by Beneš [7] and thereafter by Hooghiemstra et al. [32], and it was further developed by Blanc and co-authors, see, e.g., [9,10,34]. Although this procedure lacks in theoretical foundation, it has been very successfully applied to several systems with multiple queues. One of the objectives of this work is to provide an extensive numerical comparison between the compensation approach and PSA. For this reason, we show how PSA is applied to the two relay model at hand. Figure 3: The one step transition probabilities diagram for a few representative states Remark 4.1. Note that PSA is a powerful numerically oriented procedure, that can be applied to the asymmetric case and it can be generalised to any (finite) number of relays.

Computation scheme
For the analysis that follows, we use the transformed model, see Section 3, and the balance equations (10)- (17). For the application of the PSA procedure, we observe the following property. ρ n+k+l β(n, k, l), k, l ≥ 0, under the following assumptions: A1. 0 ≤ ρ < 1 (stability condition); A2. π k,l can be analytically continued as a function of ρ into a domain which includes the disk |ρ − 1 2 | ≤ 1 2 .
As illustrated in [9,59] and the references therein, this property is valid for among others any quasi birth-anddeath system, the JSQ system, the coupled processor system, covering also the model under consideration.
In what follows, we apply PSA to the model at hand and derive a recursive, computational scheme for the equilibrium joint queue length distribution. For more details on PSA, the interested reader is referred to [9,18,32,59], where one can find the detailed implementation of PSA for a plethora of systems. Following the steps of PSA in [10], we obtain and solve a recursive set of equations for the coefficients β(n, k, l) in (33). From this, the equilibrium distribution is computed, as well as any performance measures derived from it.
In Appendix C, we present all details of the method. For illustration purposes, we have chosen a = 1 2 , so as to simplify the notation and enhance the readability of the method. The computation scheme of the PSA is summarised in the next paragraphs.
We first substitute the power-series expansion (33) into the balance equations (10)-(17) and the normalisation equation. This leads to a polynomial expression in ρ equal to zero, which we use to equate the corresponding powers of ρ, and obtain a recursion in the coefficients β(n, k, l), n ≥ 0, (k, l) ∈Ŝ. This results in a computational scheme through which we can compute the coefficients β(n, k, l). The computation of the coefficients β(n, k, l) is equivalent, in terms of mathematical complexity, to solving a system of equations, cf. Equations (61)-(68), for a finite large value of n ≥ 0.
Having computed the coefficients β(n, k, l), we use them to approximate, to any degree of required precision, the equilibrium distribution, using (33), and thus, we can compute many performance measures by writing them as a power series in ρ. From a theoretical perspective, PSA can be used to compute the performance measures in a symbolic fashion, however, in practice, this can be achieved only for coefficients of very small order. This is mainly due to the fact that it is required to solve symbolically a recursive scheme, which is increasingly hard as the order of ρ increases. The numerical computation of the performance measures, for given values of the system parameters, is possible in principle up to an arbitrary precision, by applying the recursion until the desired precision is achieved. However, note that the given values of the system parameters need to be chosen so as to be within the radius of convergence of (33). In order to expand the range of values (within the stability region), for which PSA can produce accurate numerical results, Keane et al. [32] propose the application of a bilinear conformal transformation of the real interval [0, 1) onto itself with G ≥ 0. Using the above conformal mapping, the power series (33) is written as a power series over the parameter θ π k,l = ∞ n=0 θ n+k+l u(n, k, l), k, l ≥ 0.
The derivation of the recurrence relations for the coefficients u(n, k, l) of the new power series (35) is similar to the process for β(n, k, l), see Appendix C.

Numerical implementation of PSA
In order to compare, the two approaches, namely the compensation approach and PSA, we sketch below the algorithmic implementation of the latter. Similarly to the compensation approach, we truncate the series expression (35) to the N psa -th term.

Equilibrium analysis: The probability generating function approach
In the following, we show how one can apply the probability generating function (PGF) approach to analyse the model at hand. The analysis through the PGF approach leads to a functional equation, whose solution usually presents formidable difficulties. However, for the two-dimensional case, techniques have been developed to reduce the problem of the solution of the functional equation to standard problems from the theory of boundary value problems, see, e.g., [15,21]. Even in these cases, the analysis is lengthy and complicated, involving sophisticated complex analysis, Riemann surfaces, and the determination of some conformal mapping. In most cases, this requires numerical analysis, which makes the formal solutions to the boundary value problems less insightful. This drawback is overcome by the use of the compensation approach and the PSA method developed in the previous sections, revealing their superiority with respect to the PGF approach.
Our aim hereon is to provide the basic steps of the analysis. To this purpose, we consider the original process {Q(n), n ≥ 0}, cf. Section 2. Similarly to the analysis of the classical JSQ model, we need to account for the different transition patterns above and below the diagonal of the first quadrant, which can be seen as "the third boundary", see, e.g., [15,Chapter III], [21,Chapter 10], and [36]. Let denote the equilibrium joint queue length distribution of the original random walk {Q(n), n ≥ 0}. From the balance equations, we obtain after tedious, but straightforward algebra with The main steps of the approach are summarised below and follow considerably the lines in [15].

Comparison of methods
In this section, we compare the compensation approach and PSA method on the basis of their performance for the system at hand. The comparison is performed in terms of numerical accuracy (absolute difference in the obtained results by both methods) and time complexity. Algorithms 1 and 2 are designed so as to compute any performance measure, which depends on the equilibrium distribution, given a desired precision ε.

Comparison of numerical accuracy
With the compensation approach, we can compute an explicit expression of the equilibrium distribution as a (infinite) linear combination of product-form terms. We have proven, that the corresponding truncated linear combinations provide asymptotic expansions which improve as we include more terms. In this perspective, we can relatively control the error of the approach. However, for PSA, to the best of our knowledge, there exists no error bound due to the lack of theoretical support for this approach, see [9,10]. The error produced by the PSA implementation can be controlled somewhat by including more terms of the series. On top of that, it is sometimes unclear when this method diverges [34], but the radius of convergence of the power series can be extended using a transformation, cf. (34).
In Table 1, we depict the total expected sojourn time of a packet measured in time slots and the correlation coefficient between the queue lengths of the two relays, for various values of the load ρ. The total expected sojourn time of a packet is computed as , and as expected, as ρ increases, so do the values of E[S]. The correlation coefficient between the queue relays is computed as . For the computations performed, we choose ρ = 0.1, 0.4, 0.7, 0.9, 0.95, so as to cover lighter and heavier traffic results. Furthermore, we choose G = 1 for the PSA implementation.
As expected, as ρ increases the correlation coefficient tends to one, and as the values of Table 1 indicate, it almost behaves like a linear function of ρ. Furthermore, it is evident, from the numerical results, that both approaches produce similar outcomes (as long as the value for ρ is away from the stability region boundary ρ = 1) differing by approximately as much as the range of the precision, cf. column | CA − PSA |. However, as ρ approaches one, PSA becomes highly unstable, while the compensation approach is producing accurate results in the entire stability region. The numerical instability of PSA can be explained observing that as ρ → 1, θ → 1, cf. Equation (34), indicating that the power series expression is approaching the boundary of the region of convergence. This can be overcome, to a certain degree, by further increasing the value G ≥ 1.
Note that the compensation approach only works in the case of two relays, while PSA can be extended to any finite number of relays. As such, we conclude that the compensation approach is more suitable when working with two relays, while PSA is generalisable to a system with more than two relays.  Table 1: Total expected sojourn time of a packet and the correlation between the queue relays.

Comparison of time complexity
From the numerical implementation, it is notable that both approaches provide accurate results to a desired precision. In both approaches, the time complexity of the corresponding algorithm depends on ε through the determination of T ca or T psa , cf. Step 4 or Step 3, respectively, which we can control. Equating ε for both algorithms, we can compare the methods in terms of their algorithmic time complexity, so as to characterise the performance speed of the two approaches. We describe time complexity in the Big-O notation as a function of the input size, i.e., O(f (·)) is measured as the maximum number of elementary steps needed to perform the algorithm, provided that each step is executed in constant (or equal) time. This way, the time required for the algorithm is described independently of the numerical implementation.
We discuss the time complexity of the two approaches separately: For the compensation approach, see Algorithm 1, for a series truncation level N ca and state-space truncation T ca , the algorithm converges with order O N ca (T ca ) 2.376 . This order is obtained as follows: (i) In Steps 2 and 5, we need to compute recursively the 2(N ca + 1) roots of Equation (23). This can be done by using the Bisection method or the False position (aka regula falsi) method. With both methods, we are able to choose the bisection intervals within the interval (0, γ 0 ), which results in time complexity O(log(γ 0 )) for the computation of a single root. In order to compute all the roots, the Bisection method needs to be repeated at least 2(N ca + 1) times. Thus, the time complexity for the calculation of all the roots is of order O(N ca log(γ 0 )).
(ii) In Step 7, we need to compute recursively the coefficients c i , e i,o , e i,1 and d i , i = 0, 1, . . . , N ca . To this purpose, we need to solve a system of equations. The system is solved implementing the Coppersmith-Winograd algorithm [16], which has a complexity O(N 2.376 ca ).
(iv) In Step 9, we need to solve a system of equations. Using again the Coppersmith-Winograd algorithm. This step has a complexity of O((T ca /2) 2.376 ), and needs to be performed at least N ca times.
Thus, the time complexity of the entire step is of order O N ca T ca /2 2.376 .
(v) In the construction of Algorithm 1, we have considered a part of the state-space region in which we directly use Equation (31), cf.
Step 7, and a part of the state-space region in which we solve a linear system of the equations cf.
Step 8. This order is obtained as follows: (i) In Step 5, for 0 ≤ n ≤ N psa and 0 ≤ k, l ≤ T psa , we compute u(n, k, l). This step determines the dominant time complexity in PSA. These coefficients are obtained solving a system of equations (69)-(76), in Appendix C, by implementing for example the Coppersmith-Winograd algorithm [16]. Note that this system of equations reveals that (i) the series truncation level should be chosen such that N psa ≤ T psa , and (ii) the state-space truncation should contain all states 0 ≤ k + l ≤ T psa . All in all, this yields that the complexity is O (N psa + 1) (T psa + 1) 2 /4 2.376 .
In the above discussion, we have demonstarted that the compensation algorithm has better big-O time complexity than PSA.

Comparison of JSRQ to other routing protocols
The JSRQ routing protocol, balances the load between the two relays, but due to this balancing, it also seems to increase collisions. For this reason, in this section, we compare the JSRQ routing protocol to a single server system. To make the two systems comparable, we consider the arrival and departure probabilities of the single server system to be the same as in the JSRQ system, i.e., λ and a. Furthermore, in order to identify the comparable region for both models we use the corresponding stability conditions. The single server system is stable if λ < a, which implies a ∈ (λ, 1). The JSQR system is stable for λ < 2a(1 − a). The stability region of the JSQR system can be equivalently written as The comparison between the protocols is discussed on the basis of the following points: (i) stability region comparison: From the above discussion, it becomes evident that the length of the stability region of the single server system is 1 − λ and of the JSRQ system it is √ 1 − 2λ. For λ < 1/2, we observe that the region of the single server system contains the region of the JSQR system.
(ii) Total expected queue length comparison: For the single server system, the total expected queue length goes to infinity as a ↓ λ and to zero as a ↑ 1. For the JSQR system, the total expected queue length goes to infinity as a → a + or a → a − , defined in Equation (47). Furthermore, for λ < 1/2, it is easy to show that a − < λ < 1/2 < a + . Moreover, the two systems have identical total expected queue lengths if a = 1/2. Note that for a = 1/2 the load of the single server system, λā/λa, is equal to the load of the JSRQ system, λ(ā 2 + a 2 )/2λāa. All in all, the above analysis reveals that for a < 1/2 the JSRQ outperforms the single server system, for a = 1/2 the systems perform identically, while for a > 1/2 the single server system performs better than the JSRQ system. See Figure 4, for a depiction of the above general remarks. (iii) Consideration of collision: Our model incorporates the phenomenon of collisions, whereas the single server model does not consider any collisions between the packets. Hence in the single server model the probability of collisions between the packets is zero, whereas in the considered model when both relays transmit packets in the same slot, the packets need to be retransmitted in a later slot. Communication networks in practice adopt the technique of retransmission in case of collisions to avoid losses, hence the considered model represents reality more closely than a single server model.

Conclusions and possible extensions
In this work, we focused on the application oriented problem of characterising the queueing delay experienced in a slotted-time relay-assisted cooperative random access wireless network with collisions and join the shortest relay queue (JSRQ) routing protocol. Note that due to the collisions, there is strong interdependence among the queues of the relays, resulting in different service probabilities, when both relays are busy, than the service probability, when one of the relays is empty. Thus, the system at hand incorporates two features: the JSRQ feature and the coupled processor feature. For this system, we investigate the stability condition and apply three different methods for the computation of the equilibrium joint queue (relay) length distribution, namely the compensation approach , the power series algorithm (PSA), and the probability generating function approach. A detailed comparison of the compensation approach and PSA is presented. More importantly, we have applied the compensation approach to a random walk on the positive quadrant with transitions to a bounded region of neighbours, extending the framework of the compensation approach to a wider class of random walks (than the nearest neighbour for which the approach was originally developed).
In a future work, we plan to generalise this work in several directions. A challenging task is related to the equilibrium analysis of a cooperative network with a queue-aware transmission protocol under which each relay configures its transmission parameters based on the status of the other. Such a protocol serves towards self-aware and intelligent networks. Moreover, we also plan to characterise the delay using a multi-packet reception model instead of the collision channel model. Under such a scheme, we can have successful transmissions even if multiple nodes (relays or source(s)) transmit simultaneously, which will definitely improve the throughput performance of the network. An interesting challenging task is the investigation of the queueing delay at a random access network with an arbitrary number of relay nodes under the JSRQ routing policy. Finally, it will be interesting to investigate service policies by taking into account the state of the network, for the ultimate goal of improving the system performance.

Appendices A An alternative way for the derivation of the stability condition
Denote by A r (n) the number of arrivals at relay node r, at the beginning of slot n (i.e., A(n) = A 1 (n) + A 2 (n), with E[A(n)] = λ), and by S r (n) the number of departures from relay node r, r = 1, 2, at the end of slot n, n ≥ 0. Then, the queue evolution is as follows Equivalently, we can define a non-negative random variable U r (n), which denotes the unused service in a time slot, and rewrite the above equation as where U r (n) = max(0, S r (n) − A r (n) − Q r (n)).
Note that due to the symmetry of the model E(S 1 ) = E(S 2 ) := E(S). Note also that aā ≤ E(S r ) ≤ a, r = 1, 2. Thus,
Proof. We consider a modified model, which is closely related to the {Q(n), n ≥ 0} model and it has the same asymptotic behaviour. This modified model corresponds to a random walk on a slightly different grid, namely In the interior and on the horizontal boundary, the modified model has the same transition rates as the {Q(n), n ≥ 0} model. A characteristic feature of the modified model is that its balance equations for k + l = 0 are exactly the same as the ones in the interior (in this sense the modified model has no "vertical" boundary equations) and both models have the same stability region. Therefore, the balance equations for the modified model are Equations (10)- (13) for all k + l ≥ 0, k ∈ Z with only the equation for state (0, 0) being different due to the incoming rates from the states with k + l = 0, k ∈ Z.
To determine the γ we consider levels of the form L = {(k, l) : 2k + l = L} and letπ L = 2k+l=Lπ k,l . The balance equations between the levels are given by λ(ā 2 + a 2 )π L = 2λāaπ L+1 , L ≥ 1, Substituting (51) into (50) yields So far, we have shown that the equilibrium distribution of the modified model has a product form solution which is unique up to a positive multiplicative constant. Returning to the {Q(n), n ≥ 0} model, we can immediately assume that the solution of the balance equations (10)-(13) is identical to the expression for the modified model as given in (48). Furthermore, the above analysis implies that this product form is unique, since the equilibrium distribution of the modified model is unique.
= λ(a 2 +ā 2 ) +āa Then the lemma follows by applying Rouché's theorem to f (z) and g(z) on the unit circle.
Then, by applying the above iteratively yields It remains to prove (a) and (b) stated above.
It follows from Proposition B.3 that the sequences of {γ i } i∈N and {δ i } i∈N tend to zero as i → ∞. The limiting behaviour of the sequences is presented in Lemma B.4 and the limiting behaviour of the coefficients is treated in Lemma B.5.
Taking now the limit as γ → 0 (which also implies that δ → 0) and δ/γ → w yields (53). By applying Rouché's theorem in Equation (53), it makes it immediately evident that the resulting equation has one root inside and one root outside the unit circle. This completes the proof of assertion (i). The proof of assertion (ii) follows similarly.