{'fraction_non_alphanumeric': 0.03366457328820053, 'fraction_numerical': 0.01117798796216681, 'mean_word_length': 4.109442776109443, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability Analysis of TDD Networks Revisited: A trade-off between Complexity and Precision

In this paper, we revisit the stability region of a cellular time division duplex (TDD) network. We characterize the queuing stability region of a network model that consists of two types of communications: (i) users communicating with the base station and (ii) users communicating with each other by passing through the base station. When a communication passes through the base station (BS) then a packet cannot be delivered to the destination UE until it is first received by the BS queue from the source UE. Due to the relaying functionality at the BS level, a coupling is created between the queues of the source users and the BS queues. In addition, contrarily to the majority of the existing works where an ON/OFF model of transmission is considered, we assume a link adaptation model (i.e. multiple rate model) where the bit rate of a link depends on its radio conditions. The coupling between the queues as well as the multiple rate model are the main challenges that highly increase the complexity of the stability region characterization. In this paper, we propose a simple approach that permits to overcome these challenges and to provide a full characterization of the exact stability region as a convex polytope with a finite number of vertices. An approximated model is proposed for reducing the computational complexity of the exact stability region. For the multi-user scenario, a trade-off is established between the complexity and the preciseness of the approximated stability region compared to the exact one. Furthermore, numerical results are presented to corroborate our claims.


I. INTRODUCTION
Time division duplex (TDD) systems use a single frequency band for both uplink (UL) and downlink (DL) traffic which offers the flexibility to adjust the UL and DL channels based on their respective traffic demand. TDD provides dynamic UL and DL bandwidth allocation which allows the network to combine spectrum bands and achieve greater spectral efficiency when customers need it most. With the growth of various new applications with high traffic and data rate demand, new techniques have been investigated to fulfill the requirements of future cellular networks. The severity of this situation will increase with fifth generation (5G) cellular networks that demand the support of higher data rates and lower latency. By enabling a seamlessly adaptation of the spectrum bands and an efficient UL/DL load asymmetry, TDD can be used for improving the capacity of dense area with high mobile data demand. Besides, motivated by the use of Massive MIMO, TDD is mainly adopted in 5G systems.
The physical layer study of cellular networks provides the information-theoretic capacity of these scenarios while assuming that the queues are saturated (not empty). In this case, the rate region is defined as the set of the achievable bit rates of the system of saturated queues. However, integrating the network level (e.g. [1] and [2]) to the study has demonstrated important gains in terms of throughput and delay. Under bursty traffic arrivals, the stability region becomes a relevant measure of the queues' bit rates. The stability region is defined as the set of arrival rates that can be supported by the network under the stability constraint of the queues (i.e. as defined in section II, a queue is considered stable if its mean arrival rate is lower than its mean departure rate). The stability region has several definitions as one can see in [3]. The difference between capacity region and stability region has been the subject of interests of several studies (e.g. [4] and [5]).
Queuing stability region is the performance metric used in this work. The consideration of bursty traffic and queuing analysis is motivated by the following scenario. In TDD cellular networks, any pair of users communicate with each other by passing through the base station (BS) where the UL and DL parts of this communication compete for the same spectrum. Therefore, the BS plays the role of a relay that receives packets from the first device (source) on the Uplink (UL) stacks them in a buffer and then transmits them to the second device (destination) on the Downlink (DL) using an opportunistic scheduling algorithm. In a TDD system, the UL and DL transmissions are performed on different timeslots, hence a packet cannot be transmitted on the downlink (BS-to-destination) until it is first received by the BS on the uplink (source-to-BS). During some time-slots, even when the channel state of the BS-destination link (DL) is favorable, it may happen that the buffer at the BS is empty (i.e no need to schedule this DL). This coupling cannot be captured by a simple performance analysis at the physical layer. Hence, a traffic pattern must be included in the analysis in order to provide a more realistic evaluation of the network capacity. The particularity of this study is the simultaneous consideration of UL and DL communications which implies a relaying operation at the BS level and a coupling between the system of queues (i.e. contrarily to the majority of the existing works where either UL or DL communications are studied).
Due to the relaying aspects at the BS level, it is interesting to examine the existing works regarding cooperative and relaying networks. The approach of cooperative relaying and multi-hop wireless networks has been studied at: the physical layer (e.g. [6] and [7]) and the traffic layer (e.g. [8] and [1]). However, these works consider only simple scenarios (i.e. three-node network or one relay and one destination scenario). Indeed, the majority of the works in the area of relaying networks consider an ALOHA random access system: (i) [1] and [9] use stochastic dominance technique, (ii) [10] considers energyharvesting capabilities and computes the relaying parameter that maximizes the stable throughput rate of the source and (iii) [11] proposes two queuing strategies for cooperative networks. Furthermore, wireless multiple-access system with probabilistic channel receptions was studied in [12] for a network of N sources communicating with one destination; the impact of a protocol-level cooperation is investigated in such wireless network. Authors in [8] have introduced a new cognitive multiple-access protocol in a relay assisted network. Both [8] and [12] take into account the time division multiple access (TDMA) as a scheduling policy. Furthermore, the majority of these studies consider simple networks (i.e. consisting of one source to destination communication aided by a relay node) due to the difficulty that one can face while characterizing the stability region for multi-user cooperative networks.
Both the bursty traffic and the relaying role of the BS lead to a coupling between the queues in the system such that the service rate of each queue will depend on the state of the other queues. The stability region characterization of the system of interacting queues has been a challenging problem and has received the attention of researchers (i.e. especially for the case of multiple-access channel networks). For the slotted ALOHA system, several approximated models were proposed (e.g. [13], [14] and [15]) and some bounds for the ergodicity region have been obtained (e.g. [16] and [17]). Several tools were proposed in the literature to overcome the challenges of the system of interacting queues. The stochastic dominance technique was introduced in [18]. It consists of elaborating a simple way for capturing the interaction between the queues by studying simple auxiliary systems of queues that dominate the system of interest. This technique is especially applied to simple scenarios in order to: (i) obtain the stability region (e.g. [19], [20] and [21]) and (ii) characterize the stability bounds for the system of interacting queues (e.g. [22]). The stability of a coupled system of queues where the service rate of each queue depends on the number of customers in all the other queues was extended in [23]. Their stability approach consists of deriving marginal drift criteria for multi-class birth and death processes. This analysis was adopted in [24] for proposing a coupled processors model that was applied for the study of underlay device-to-device communications in [25]. Furthermore, to deal with the coupling between the queues, authors in [26] and [27] propose a decomposition model and an iterative approach in Stochastic Petri Nets.
The majority of the previous works consider ALOHA random access (i.e. without a centralized scheduling) contrariwise to our scheme where a TDD cellular system with a centralized TDMA scheduling is investigated. Preceding efforts were limited to simple scenarios (i.e. three-node scenario or single destination scenario) with an elementary ON/OFF model of transmission. However, in this work we analytically characterize the stability region of TDD cellular systems (i.e. both UL and DL communications) with a TDMA scheduling and where both multi-user and multiple rate model are taken into account.
In this paper, we consider the case of a discrete-time slotted TDD system with 2K + U users (UE), each of which has a buffer of infinite capacity. U users communicate with the BS (denoted by UE2BS communications in the sequel) and K pairs of users communicate with each other by passing through the BS (denoted by UE2UE communications in the sequel). Under a TDMA policy, we characterize the stability region of the network. As we mentioned before, what makes the problem challenging is the BS functionality as a relay between the UL and the DL of the UE2UE communications (i.e. a packet cannot be transmitted on the DL if it is not received by the BS on the UL). Hence, the performance of the queues will depend on the state of the queues (empty or not) at the BS level.
Most of the works concerning the analysis of stability region in cellular network are focused on the strict consideration of either uplink or downlink communications. However, in our scenario we consider UE2UE communications that take place on both uplink and downlink in such a way that the BS plays the role of a relay between the source UE and the destination UE. We revisit the queuing analysis of such TDD scenario in order to characterize its stability region as a convex polytope with a limited number of vertices and to derive the existing threshold between the complexity and the precision of this characterization. The key particularities of this work are summarized as follows: • In the literature, the stability region analysis of interacting queues with relaying functionality was mainly analyzed either: (i) based on stochastic dominance technique for finding the stability region (mainly for aloha systems with three-user scenario and ON/OFF transmission model) or at least describing the necessary and sufficient conditions for queuing stability or (ii) via numerical and simulation results. However, in this paper we approach the problem in a different way and we transform the multidimensional Markov Chain that models the network to multiple one dimensional (1D) Markov Chain models (as detailed in the sequel). After evaluating the different 1D Markov Chains, we characterize the exact stability region as a convex polytope with a limited number of vertices. This exact stability region turns to have a high complexity. Therefore, we propose approximated models characterized by having an explicit analytic form and a low complexity while a high precision is guaranteed (i.e compared to the exact stability region). We start by considering the simple scenario that consists of one UE2UE and one UE2BS communication. In theorem 2, we characterize the exact stability region of this simple scenario which turns to have a non explicit complex form. Hence, we propose a closed-form tight upper bound of the exact stability region (see theorem 4). Moreover, for the general scenario of multiple UE2UE and UE2BS communications, we evaluate the exact stability region and we discuss its computational complexity (see theorem 6). Thus, we propose two techniques to respectively reduce the complexity of the exact stability region (see theorems 7 and 10). We deduce the trade-off between the complexity and the precision of the stability region computation. The remainder of this paper is organized as follows. Section II describes the system model. The stability region analysis for a three-UEs simple case is provided in Section III in order to have a clear perception of the advantages of overlay D2D compared to cellular communications. The stability region of the multi-UE general case is presented in Section IV as a simple convex polytope. Numerical results are presented in Section V. Section VI concludes the paper whereas the proofs are provided in the appendices. 1

II. SYSTEM MODEL
In this study, the considered scenario consists of two sets of communications: (i) set of UE2BS communications (where users are transmitting packets to the BS e.g. to contact users in other cells or to have internet connections, etc.) and (ii) set of UE2UE communications (between pair of UEs). UE2UE communications are performed by passing through the BS.
We consider a single cell scenario with K UE2UE communications and U UE2BS communications. In other terms, we suppose K + U communications and 2K + U users in the cell. We denote by UE i,s and UE i,d the pair source and destination users corresponding to the i th UE2UE communication (for all 1 ≤ i ≤ K) and by UE j,u the user corresponding to the j th uplink UE2BS communication (for all 1 ≤ j ≤ U ). Let us describe the cellular scenario illustrated in Figure 1.
In this scenario, if two devices want to communicate with each other, they must exchange their packets through the BS. The communication between UE i,s and UE i,d is performed through the BS (for all 1 ≤ i ≤ K) such that the BS transmits packets to the destination user UE i,d only if it receives them from the source user UE i,s . Therefore, for each UE2UE communication corresponds a buffer at the base station that could be empty during some time slots. Hence, a coupling exists between the queues such that the service rate of the users' queues depend on the state being empty or not of the BS which makes the queuing stability analysis challenging. In this scenario the following links exist: link i,s : UE i,s -BS; link i,d : BS -UE i,d and link j,u : UE j,u -BS (with 1 ≤ i ≤ K and 1 ≤ j ≤ U ). The considered network consists of two type of communications UE2UE communications and UE2BS communications. In practice, these two types of communications coexist with downlink communications (i.e. between the BS and the users) that aim to deliver downlink traffic to mobile users (e.g. internet connections or receive calls from users in other cells etc.). The extension of this work to the scenario where this downlink traffic is taken into account is straightforward.

A. Priority Policies
We call priority policy the sorting of the communications' priorities according to which the users are chosen for transmission. In other terms, among the users that are able to transmit (which means have some packets to transmit and have the required radio conditions to do it), the UE that is chosen to transmit is the one that has the highest priority according to the considered priority policy. Hence, a user is scheduled only when all the more prioritized users are not able to transmit. We denote by Ω Γ the set of all the possible priority policies. Γ ∈ Ω Γ denotes a priority policy according to which the users are chosen for transmission. One can see that for N communications the number of possible priority policies is given by the number of existing permutations: N !.
Note that any other scheduling policy is nothing but a convex combination of these priority policies. For this reason, our work is based on studying these priority policies that characterize the corner points of the stability region. Any other scheduling corresponds to an interior point of the characterized stability region.
These priority policies allows us to avoid the multidimensional Markov Chain modeling of the interacting queues. Thus, for a given priority policy, each queue can be modeled by a one dimensional (1D) Markov chain. Our approach transforms the multidimensional Markov Chain, that captures the dependency between the queues, to a 1D Markov Chain model for a given priority policy. However, the modeling that we propose remains challenging due to the coupling between the queues. As one can see later, an additional analysis is required for capturing the interaction between the queues.

B. System of queues
We consider a system of queues to describe the studied scenario. UE j,u (for all 1 ≤ j ≤ U ) communicate with the BS through an uplink cellular communication and the queue of user UE j,u is represented by Q j,u . The communication between UE i,s and UE i,d (for all 1 ≤ i ≤ K) is represented by two consecutive queues: the uplink queue Q i,s of UE i,s and the download queue Q i,BS (see figure (2)). The BS does not transmit to UE i,d unless it has received at least one packet from UE i,s . This coupling between the queues induces that the service rates of all the queues Q i,s and Q i,u depend on the state (empty / not empty) of each Q i,BS which makes the queuing stability analysis challenging. The users' queues are assumed saturated. The traffic arriving to the queues Q i,s and Q j,u is time varying, i.i.d. over time and with rate respectively equal to λ i,s and λ j,u for 1 ≤ i ≤ K and 1 ≤ j ≤ U . The traffic arriving to Q i,BS is nothing but the departure from Q i,s .
The traffic departure from the users' queues is also time varying and depends on the queues' states (empty or not), the scheduling allocation decision and the time varying channel conditions. For a given priority policy Γ, the average service rates of the queues Q i,s , Q i,BS and Q j,u are respectively denoted by µ i,s (Γ) and µ i,d (Γ) and µ j,u (Γ). The vector that describes the service rate of the users' queues for a given policy Γ is the following: The vector that describe the arrival rates of the users' queues is given by:

C. Set of bit rates
We consider an adaptive modulation scheme such that the transmission rates is improved by exploiting the channel state at the transmitter. The SNR values are divided into a finite set of intervals [S 1 , ..., S M ] where the j th interval S j = [γ j , γ j−1 ] (called hereinafter j th state) is characterized by two SNR thresholds γ j−1 and γ j such that γ j−1 > γ j . Here, the adaptive modulation consists on considering a finite set [r 1 , ..., r M ] of bit rates as a mapping of the SNR intervals. Therefore, if a link has a SNR within the state S j than the bit rate if this transmission is r j . It is worth mentioning that to transmit at a rate r j in the downlink or the uplink, the SNR states and thresholds may be different. In order to deal with that we use γ U L j and γ DL j to describe the SNR thresholds for respectively the uplinks and downlinks.

D. Probabilities for channel quality
The channel between any two nodes in the network is modeled as a Rayleigh fading channel that remains constant during one time slot and changes independently from one time slot to another based on a complex Gaussian distribution with zero mean and unit variance. The received SNR for a link i is given by: where h i is the fading coefficient, d i is the distance between source and destination, P i is the transmission power, β is the path loss exponent and N 0 is the noise. We denote by p n i the probability that the SNR of the link i is within n th interval (i.e. in state S n ): Given a complex Gaussian distribution of the channel h i then: Note that P [SN R i ≥ +∞] = 0 and P [SN R i ≥ −∞] = 1. Let us consider the following notation

E. Stability region definition
Based on the definition in [3], the i th queue is stable if its length Q i (t) satisfies: It follows from Loyne's theorem [28] that if the arrival and service process of the i th queue are strictly jointly stationary then this queue is stable if λ i < µ i where λ i and µ i denote the mean arrival and service rate of the queue Q i . Hence for a given scheduling policy ξ, the stability region R (ξ) is given as the closure set of arrival rate vectors for which all the queues are stables under the scheduling policy ξ.
where ≺ is the component-wise inequality, λ is the vector of the queues' average arrival rates and µ (ξ) is the vector of the queues' average service rates under the scheduling policy ξ.
The stability region R is defined as the union of the stability region for all the feasible scheduling policies (denoted by Ω ξ ).

F. Our approach for computing the stability region
We based our approach on the fact that all the existing scheduling policies can be written as a convex combination of the priority policies. This gives the idea to study the stability region based on the defined priority policies. Thus, the stability region R is given by the union of the R (Γ) over all the feasible priority policies (denoted by Ω Γ ). Note that any region R (ξ) for any other scheduling ξ ∈ Ω Γ is a subset of the region R (Γ) due to the fact that this scheduling policy ξ can be written as a convex combination of the set of the priority policies Ω Γ .
Hence, for characterizing the stability region we have to increase as much as possible the service rates of the users' queues. For the i th UE2UE communications, this means increasing the service rate of the UL side (UE i,s -BS) which implies the increase of the arrival rate at the DL side (BS-UE i,d ) and risks the loss of the stability of the BS queue Q i,BS . Recall that the service rates of both UL and DL sides are coupled through the scheduling. Instability of Q i,BS means that many packets will not be delivered to UE i,d and the network becomes unstable (it means that average delay is not finite).
Here, the challenge is to characterize the stability region by finding the priority policies that maximize the service rates of the UL side ( hence achieve the border points of the stability region) while guaranteeing the Q BS stability. This coupling between the stability region and Q BS is very challenging in general in queuing theory (due to scheduling) and especially in the context of time varying wireless channels (i.e. fast fading) where the allocated bit rate changes according to the time varying channel state.
If we are able to find the policies that achieve the corner points of the stability region and we denote these policies by Ω * Γ then this set of policies is sufficient to describe the stability region of the system. Therefore, for each priority policy Γ ∈ Ω * Γ , we find the region that capture the coupling between the queues and guarantee their stability. Hence, based on the definition (1) of the operator co, the stability region can be simply characterized by these corner points and given by: {x i }, we define co (S) as it follows: This work is divided into two sections. In the first section, we consider a simple scenario that consists of three users and contains one UE2UE communication and one UE2BS communication. For this three-UEs scenario, we consider a set of three bit rates [r 1 , r 2 , r 3 ] such that the transmission rate of a link takes a value within this set depending on its channel state. We start by providing the exact stability region of this scenario. This stability region does not have an explicit form and turns be computationally complex (especially due to the consideration of a set of three bit-rates and not a ON/OFF model). In order to reduce complexity, we find an upper bound of this stability region and prove analytically that it is a very close approximation of the stability region. This upper bound that has a simple and explicit form reduces the complexity of the stability region computation. This simple scenario brings several interesting insights on the stability region of these cellular scenarios as well as the complexity of such study. This motivates us to elaborate a simple, explicit and good approximation of this stability region.
In the second section, we consider the general scenario with K UE2UE communications and U UE2BS communications. For this general scenario, a set of two bit rates [r 1 , r 2 ] is considered due to the computation complexity of considering a set of three bit rates for this general case. Similarly to the simple case, we provide the exact stability region of the general scenario and we show the complexity of such computation. Hence, we propose a tight and explicit approximation that decreases the complexity of the stability region. The main novelty of this work is to capture, for such cellular scenarios, the trade-off between the complexity and the precision of the stability region calculation.

A. 3-UEs scenario description
The simple scenario (SS) consists of considering a single cell wireless network with three users in the cell UE 1,s , UE 1,d , UE 1,u respectively denoted by UE s , UE d , UE u . We assume that user UE s wants to communicate with UE d , while UE u is transmitting packets to the BS. The communication between UE s and UE d is performed through the BS. Hence, we consider the following three links: link s UE s -BS, link d : BS -UE d , link u : UE u -BS. The 3-UEs scenario is illustrated in Figure 3. The BS does not transmit to UE d unless it has received at least one packet from UE s . This coupling between the queues induces that the service rates of Q s and Q u depend on the state (empty / not empty) of Q BS which makes the queuing stability analysis challenging. The traffic arriving to Q BS is nothing but the departure from Q s . The traffic departure is also time varying and depends on the scheduling allocation decision and the time varying channel conditions. The departure average rates from queues Q s and Q u are respectively denoted by µ s and µ d .
For this simple scenario, a set of three bit rates is considered [r 1 , r 2 , r 3 ] with r 1 = kr 2 and r 3 = 0 (with k ∈ N + ). For this set of bit rate corresponds a set of SNR intervals [S 1 , S 2 , S 3 ]. As we mentioned before in equation (1), we use the notation p n i to describe the probability that the SNR of the link i is within the n th SNR interval (i.e. in state S n ). We suppose a complex Gaussian distribution of the channel h i ; hence these probabilities can be easly derived from equation (2). For the downlink BS-UE d , these probabilities are the following: For both uplinks UE s -BS and UE u -BS, these probabilities are the following: with i = s for the UE s -BS link and i = u for the UE u -BS link.
In the aforementioned expressions, d s is the distance between U E s and the BS, d u is the distance between U E u and the BS, d d is the distance between U E d and the BS. In the sequel, we refer to the three-UEs scenarios with 3 bit rates (r 1 , r 2 and r 3 ) by the simple scenario SS .

B. Organization
For the 3-UEs scenario study, the organization is as follows: section III-C provides a theoretical analysis of the exact stability region (SR) of the system. We note that deriving this SR is computationally complex. Therefore, we propose in section III-D an explicit and simple expression of an upper bound (UB) for this SR. Hence, in section III-E we compare the SR to its UB by expressing the maximum difference between both of them. Thus, this UB turns to be a simple and very close approximation of the SR.

C. Exact stability region for the 3-UEs scenario
In this section, we characterize the exact stability region for the 3-UEs cellular scenario. Assuming a user scheduling such that only one communication is possible in each timeslot then for a given timeslot, only one of the following communications is scheduled: UE u -BS or UE s -UE d . Recall that Γ is the priority policy according to which the links are sorted. We know that in order to characterize the stability region it is sufficient to consider its corner points that correspond to the "extreme policies" where the priority is always given to the same communication or when the priority is always given to the communication that has the better channel state. In this paper we denote by Ω ss Γ the set of the priority policies corresponding to the corner points of the stability region. Here, we only indicate that this set contains 6 policies {Γ 1 , Γ 2 , Γ 3 , Γ 4 , Γ 5 , Γ 6 } and that finding the values for this set of policies leads to the characterization of the stability region. We exhaustively describe these policies in the appendix A. However, the problem is not simply solved by finding the subset Ω ss Γ . Actually, for each priority policy Γ ∈ Ω ss Γ , the challenge remains in capturing the coupling between the queues due to the relaying functionality of the BS between the UL and DL traffic of the UE2UE communications.
Indeed, the priority policies' approach transforms the multidimensional Markov Chain, that models the interacting queues, to a one dimensional (1D) Markov chain for each priority policy. In order to characterize the stability region, all the priority policies should be considered or at least the priority policies that achieve the corner points of the stability region. Therefore, the queues can be modeled by a 1D Markov Chain per priority policy. Nonetheless, the proposed modeling (1D Markov Chain), for each priority policy, remains challenging since the service rates of the queues still depend on the state of the BS queue due to the coupling between the UL and DL traffic of the UE2UE communications. Thus, the following additional analysis is required for capturing the dependency between the queues.
For a given priority policy Γ, the service rate of the queues depend not only on the SNR states of the links but also on the state (empty or not) of the BS. For this reason, we consider a queuing theory approach in order to capture this coupling between the service rate of the queues (Q s and Q u ) and the state of Q BS . In section V, we expose examples that show the impact of queuing approach on the performance analysis of the network.
Q BS might be empty at some timeslot, therefore when UE s -UE d communication is scheduled then the choice between UL (UE s -BS) or DL (BS-UE d ) depends not only on the SNR states of these two links but also on the state (empty or not) of Q BS . In order to take that into account, we introduce a new parameter α called fraction vector. It is used to divide the time where the UE 1 -to-UE 3 communication is scheduled between its UL and DL such that the stability of the queue Q BS remains satisfied.
The UE s to UE d communication, modeled as a chain of queues Q s -Q BS , might be in 4 different SNR states: (S 1 , S 1 ) , (S 1 , S 2 ) , (S 2 , S 1 ) , (S 2 , S 2 ) where the SNR state S i is defined in section II. For each couple of SNR states for the links (UE s -BS, BS-UE d ), we define a parameter α i (with 0 ≤ α i ≤ 1) that corresponds to the fraction of time that the resources are allocated to the UL (UE s to BS) whereas 1 − α i corresponds to the fraction of time that the resources are allocated to the DL (BS to UE d ). It gives that α = (α 1 , α 2 , α 3 , α 4 ). The fraction vector α is considered only for these 4 couples of SNR states where both UL and DL of the UE s -to-UE d communication have a SNR state different than S 3 thus a non zero rate. Only for these cases, the concept of fraction parameter makes a sense for dividing resources between both UL and DL that are able to transmit. However, for the other combinations of couple SNR states (where at least one SNR state is equal to S 3 ), this concept does not hold because only the link with a positive rate is able to transmit.
For different priority policies of UE u -BS and UE s -UE d communications, different fraction of resources will be allocated to the UL and DL of the chain Q s -Q BS which corresponds to different values of α ∈ [0, 1] 4 . For each priority policy, we find the optimal fraction vector α * = (α * 1 , α * 2 , α * 3 , α * 4 ) that achieves one corner point of the stability region. Finding α * allows us to avoid the need to vary α ∈ [0, 1] 4 for each priority policy in order to obtain the corresponding corner point. To make expression simpler, we use the following notation for a given priority policy Γ: • U (Γ) and V (Γ) as the probabilities that UE s transmits respectively at rate r 1 and r 2 when Q BS is empty.
• W (Γ) and X (Γ) as the probabilities that UE u transmits respectively at rate r 1 and r 2 when Q BS is empty.
• Y (Γ) and Z (Γ) the probabilities that UE u transmits respectively at rate r 1 and r 2 when Q BS is not empty.
Theorem 2. The stability region for the 3-UEs cellular scenario is given by (λ s , λ u ) ∈ R ss c s.t.: where Ω ss Γ is the set of the priority policies for the simple scenario that achieve the corner points of the stability region (with |Ω ss Γ | = 6). The queues' service rates µ s (α, Γ) and µ u (α, Γ) are respectively given by (3) and (4) for all the priority policies Γ ∈ Ω ss Γ .
With Π 0 BS (Γ) the probability that the queue Q BS is empty under the priority policy Γ and the six values of the parameters U, V, W, X, Y, Z corresponding to the six border priority policies Γ ∈ Ω ss Γ are given in table I. We can resume the procedure of finding the stability region by the algorithm below. The complexity of the following computation as well as the non existence of an explicit form of the exact stability region comes from the consideration of a set of three bit rates which makes complicated the Markov Chain model of the queue Q BS . Thus, finding the probability that this queue is empty Π 0 should be deduced from the solution of a system of linear equation. Due to that, we do not have an explicit form of the service rates of the queues. Hence, we are not able to find the optimum fraction vector α * for each priority policy that achieves the corner point corresponding to this policy and avoids us to consider all the values α ∈ [0, 1] 4 .
Algorithm 1 Procedure of the exact stability region for all α ∈ [0, 1] 4 do 3: Find the probability Π 0 by solving the system of linear equations 4: Deduce the service rates of the queues 5: end for 6: end for The computational complexity of the above form of the stability region comes from two factors: (i) the need of varying the fraction vector α within [0, 1] 4 in order to find the stability region and (ii) the calculation of the probability Π 0 BS that the queue Q BS is empty for each value of α. In fact, considering a set of three bit rates makes the Markov Chain modeling of the BS queue Q BS complicated as shown in the proof of the above theorem. It yields to a complicated computation of the probability Π 0 BS which is deduced by solving a linear system of equation for each value of the fraction vector α. On the aim of reducing the complexity and having an explicit form, we propose a simple analytic expression of an approximation of this stability region.

D. Approximated stability region
As we mentioned before, finding the exact stability region needs a lot of computation (for finding Π 0 BS ) and an exhaustive variation of α ∈ [0, 1] 4 . These are the main two causes of the computational complexity of the exact stability region. Actually, the exact stability region has not an explicit analytic form due to the fact that it is deduced from the numerical solution of a system of equations (for finding Π 0 BS ) for each value of α ∈ [0, 1] 4 . Therefore, it is important to propose an epsilon-close upper bound ( -UB) that has an explicit form and that is simple to compute. An epsilon-close upper bound presents an upper bound of the exact stability region that has epsilon as the maximum distance between the exact region and its upper bound. Therefore, we start by simplifying the Markov Chain modeling of the queue Q BS in order to derive an approximation of the stability region that has a close-form analytic expression and that is simple to compute. After computing the small gap between the exact and approximated stability region, the latter one turns to be a close upper bound of the exact stability region.
Considering three different bit rates generates a complicated Markov Chain for modeling the queue Q BS (Upper side of Fig. 4). As one can see in the proof of the theorem 2, finding the probability that Q BS has no packet to transmit (Π 0 BS which is mandatory for finding the stability region) is challenging for this Markov Chain model and has not an explicit form. The reader is invited to refer to the proof of theorem 2 for more details on this Markov Chain.
Therefore, we propose a simple approximated Markov Chain for modeling the queue Q BS (Lower side of Fig. 4). This approximated model is a simple birth and death Markov Chain where passing from state x i to x i+1 (receiving a packet) corresponds to the average probability of receiving a packet at both rates r 1 and r 2 . Moreover, passing from state x i to x i−1 (transmitting a packet) corresponds to the average probability of transmitting a packet at both rates r 1 and r 2 . We recall that the multiple rate model is still taken into consideration at the transition probabilities level of the approximated Markov Chain. This approximated model has mainly two advantages: (i) an explicit form of the probability of being at state 0 and (ii) an explicit form of the fraction optimal fraction vector α that achieves the corner point of the stability region. Hence, compared to the exact stability region, this approximation permits to avoid the two main reasons of complexity as it follows: computation of Π 0 BS by applying an explicit formula and not by solving a system of equations and the lack of the need of varying α within [0, 1] 4 due to the explicit expression of the optimal fraction vector that achieves the corner point of this approximated stability region. The details of the explicit expression of the approximated model are given in the proofs of lemma 3 and theorem 4. In the sequel we characterize the approximated stability region by proceeding in two steps: (i) we derive in lemma (3) the condition that α should satisfy for each priority policy then (ii) we provide in theorem (4) the corner points that correspond to the priority policies Γ ∈ Ω ss Γ and that describe this approximated stability region.
We distinguish the results of the approximation from that of the exact stability region by adding the following indication˜over all the used notation. Recall that the set of bit rates is given by [r 1 = kr 2 , r 2 , r 3 = 0]. Lemma 3. For a given priority policy Γ, the fraction vector (α 1 ,α 2 ,α 3 ,α 4 ) should verify: Theorem 4. The approximated stability region for the 3-UEs cellular scenario is the set of (λ s , λ u ) ∈R ss c such that:R where the queues' service ratesμ s (Γ) andμ u (Γ) are respectively given by (6) and (7), S α is a limited subset that is simply computed by (8). Γ ∈ Ω ss Γ where Ω ss Γ is the set of priority policies that achieve the border of the stability region (with |Ω ss Γ | = 6). The values of the parameters U (Γ), V (Γ), W (Γ), X (Γ), Y (Γ) and Z (Γ) for the six priority policies Γ ∈ Ω ss Γ are given in table I. Recall that Proof. See Appendix-C.

E. Comparison real and approximation
The exact stability region has not an explicit form because of the complicated Markov Chain of the queue Q BS that has not an explicit expression of 0 Bs (but a solution of a system of equations). Moreover, this complexity is accentuated by the need of varying α within all the interval [0, 1] 4 in order to elaborate the exact stability region. However, the approximated stability region has an explicit expression of both parameters: (i) probability that the BS queue is empty and (ii) the optimal fraction vector that achieves the corner point of this region. Hence, the importance of this approximated model lies on the proposition of a simple analytic form of the stability region. Furthermore, we verify that the approximated stability region is an epsilon-close upper bound of the exact stability region. To do so, we demonstrate analytically that the relative error between these two regions is positive and bounded.
Theorem 5. For the 3-Users cellular scenario, the approximated stability regionR ss c is a close upper bound of the exact stability region R ss c with a maximum relative error * R . Therefore, R ss c is bounded as it follows: and the following parameters depend on the priority policy Γ: •α * the optimum fraction vector computed analytically in S α (given by (8)) and that achieves the corner point of the approximated stability region. • a 01 and a 02 are the arrival probabilities at Q BS when this queue is empty with a respective rates r 1 and r 2 (see equations (13) and (14)). • a 11 (α * ) and a 12 (α * ) are the arrival probabilities at Q BS when this queue is not empty with a respective rates r 1 and r 2 (see equations (15) and (16)). • b 11 (α * ) and b 12 (α * ) are the departure probabilities from Q BS at rate r 1 and r 2 respectively. They depend on the fraction vector α (see equations (11) and (12)).
Proof. See Appendix-D.
Hence, for each priority policy we find (Γ) which correspond to the deviation between both the exact corresponding corner point and its approximation. We verified that (Γ) ≥ 0 for all Γ ∈ Ω ss Γ , henceR ss c presents an upper bound of the exact stability region R ss c . As one can see in the numerical section, the epsilon difference between the exact stability region and its upper bound is small. Hence, we highly reduce the complexity by finding an explicit and close upper bound of the exact stability region.

IV. MULTI-UE SCENARIO
For the general scenario, we consider K UE2UE communications between pair of UEs (UE i,s and UE i,d ) and U UE2BS communications between BS and UE i,u . That means that in total 2K + U users are considered in the cell. For the multi-user case, the study is applied for a set of two bit rates [r 1 , r 2 ] with r 2 = 0. Unfortunately, considering a set of 3 bit rates as in section III is complex and hard to compute and to simplify. Thus, we study the case of two bit rates. Even for this case, the characterization of the exact stability region, given by theorem 6, remains computationally complex for the multi-user case. Therefore, we limit the complexity of the exact stability region in theorem 7). However, the complexity remains important. Hence, we propose in theorem 10 a simple approximation of the exact stability region that is characterized by the following: (i) highly reducing the complexity of the exact stability region and (ii) being an -close approximation of the exact stability region (it means that the maximum distance between the approximated and exact stability regions is equal to a small number ). The importance of this result lies on shifting a very complex problem to a simple approximated one with a low complexity and a high preciseness. A trade-off exists between the precision of the stability region and the complexity of characterizing this region.
In order to to simplify the presentation of the results in the sequel, we use the following notation of the transmission probabilities at rate r 1 : We denote by Γ the priority policy under which the users are sorted. The scheduling policy consists on choosing a UE if and only if all the more prioritized UEs are not able to transmit. We denote by Ω Γ the set of all the possible sorting of the users. the number of existing priority policies consists of the number of the possible permutation of K + U communications. Hence, |Ω Γ | = (K + U )!. We define the following two sets for each communication i for all 1 ≤ i ≤ K + U : The vector that describe the arrival rates of the users' queues is given by: Here, we derive the stability region of the multi-UE cellular scenario where the devices that want to communicate with each other must exchange their packets through the BS. We consider the following notation of the service rate: • µ i,s (Γ): service rate of the uplink of the i th UE2UE communication • µ j,u (Γ): service rate of the downlink of the j th UE2BS communication The traffic departure is time varying and depends on the scheduling allocation decision and the time varying channel conditions. The departure average rates from all the users queues for the cellular scenario is denoted by: We assume a user scheduling such that only one communication is possible in each timeslot. Q i,BS might be empty at some timeslot, therefore when the UE i,s -UE i,d communication is scheduled then the choice between uplink (UE i,s to BS) or downlink (BS to UE i,d ) depends not only on the SNR states of these two links but also on the state (empty or not) of the corresponding queue at BS Q i,BS . In order to take that into account, we introduce in the analysis a new parameter α i ∈ [0, 1] for each UE2UE communication (1 ≤ i ≤ K) that describes the fraction of time that the resources are respectively allocated to the uplink (UE i,s to BS). In other terms, 1 − α i is the fraction of time that the resources are allocated to the downlink (BS to UE i,d ). The fraction vector is mainly introduced for guaranteeing the stability region of the queues at the BS level. We denote by α the fraction vector that considers all the UE2UE communications. α = [α 1 , α 2 , ..., α K ]

A. Organization
In the sequel, we start by characterizing the exact stability region of the multi-UE scenario in section IV-B. This region is given by considering all the feasible priority policies Γ ∈ Ω Γ and all the values of α ∈ [0, 1] K are considered. Since computing this region is highly complex. Thus in section IV-C, we start by reducing the complexity by limiting the number of fraction vector α that should be considered for characterizing the stability region. Furthermore, an epsilon-close approximation with a simple explicit form and a very low complexity is proposed. A trade-off is elaborated between the complexity and the precision of the stability region computation.

B. Exact stability region
The theorem below gives the exact stability region of the multi-UE scenario. The complexity of the exact stability region is discussed directly below the theorem 6.
Theorem 6. For the multi-UE cellular scenario, the exact stability region is given by the set of λ ∈ R c such that: . Ω Γ is the set of all the possible priority policies (with |Ω Γ | = (K + U )!). The elements of µ (Γ) which are µ i,s (Γ) and µ j,u (Γ) (for 1 ≤ i ≤ K, 1 ≤ j ≤ U ) are respectively given by (9) and (10).
Proof. See Appendix-E.
The exact stability region is computed by following the instructions below: * Compute the service rates of the queues * Deduce a point of the exact stability region • Exact stability region as the convex hull of all these points Based on the steps above for the exact stability region computation, we deduced the high complexity of this approach. Indeed, the complexity of the exact stability region computation comes from two factors: the high number of priority policies (depending on the number of communications) as well as the large number of the fraction vector α values that should be considered for each priority policy. Actually, the number of the considered computing the exact stability region pushes us to consider all the the possible priority policies Γ ∈ Ω Γ . Having K + U communications means that there exists (K + U )! possible policies to sort them. Further, for each priority Γ, we should vary α ∈ [0, α * ]. Suppose that for each α i Thus, for each priority policy we consider L K values of the fraction vector α. It is clear that bigger is L higher is the precision of the stability region (in numerical section L ≥ 10 3 ). we deduce that even by considering only a set of two bit rates [r 1 , r 2 = 0], the complexity of the stability region computation remains high: L K (K + U )!.
Therefore, we start reducing the complexity of the exact stability region by verifying that only the border points of α i ∈ {0, α * i } (for 1 ≤ i ≤ K) are sufficient for characterizing the exact stability region. Furthermore, we propose an approximation of the exact stability region that highly decreases the complexity with a tight loss of precision. This study show a trade-off between the precision of the stability region computation and its complexity.

C. Precision versus complexity
The main challenge is to reduce the complexity of the stability region computation while guaranteeing its precision. As we have discussed before, the complexity of the exact stability region is mainly caused by the high number of the following two parameters that should considered for the characterization of this region: (i) priority policies and (ii) fraction vector α values. These two factors of complexity are respectively studied in this section. In the first part, we reduce the complexity by limiting the values of the fraction vector α that should be considered for each priority policy. We prove that a limited set of values of the fraction vector α is sufficient for characterizing the exact stability region. However, the complexity remains high due to the large number of the priority policies that should be considered (K + U )!. Hence, in the second part we propose an epsilon-close approximation of the exact stability region of the symmetric case (see definition (8). This approximation highly reduces the complexity by limiting the number of the considered priority policies while guaranteeing a high precision.
Theorem 7. For the multi-UE cellular scenario, the exact stability region is the set of λ ∈ R c such that R c can be simplified as it follows: where the elements of µ (Γ) which are µ i,s (Γ) and µ j,u (Γ) (for 1 ≤ i ≤ K, 1 ≤ j ≤ U ) are respectively given by (9) and (10). Where the limited set S α of 2 K elements reduce the complexity.
Proof. See Appendix-F. In theorem 7, the complexity of the exact stability region is decreased due to the fact that for each priority policy Γ ∈ Ω Γ , the border values of each α i are sufficient for characterizing the corner point of the exact stability region that corresponds to Γ. Indeed, for a given priority policy Γ ∈ Ω Γ , the service rates vector µ (α, Γ) for any α ∈ [0, α * ] can be written as a convex combination of the service rates vectors µ (α, Γ) with α ∈ S α . Actually, the L values of each α i that we should consider in theorem 6 for each priority policy Γ is reduced to 2 values in theorem 7 (0 and α * i ). The proof of the theorem 7 shows that varying α within the finite set S α is sufficient for the characterization of the exact stability region. Hence, the computational complexity is reduced from L K (K + U )! to 2 K (K + U )!.
However, the complexity of the exact stability region remains high because of the large number of the priority policies that should be considered ((K + U )!). Hence, for the symmetric case (defined in (8)), we reduce further the complexity by limiting the number of the considered priority policies for the characterization of the stability region. An -approximation, such that the maximum distance between the approximated and exact stability region is limited to , is proposed. The main importance of this approximation is that it highly reduces the complexity (in terms of number of the considered priority policies and number of the fraction vector values) while guaranteeing the precision of the result.
Definition 8. The symmetric case consists of considering all the UEs of the cell at the same distance d from the BS. Hence, this symmetric case is defined by the following values: p s = p i,s = p j,u and Definition 9.R is called an -approximation of a stability region R (with 0 ≤ ≤ 1) iff the following is verified:R ⊆ R ⊆R + The symmetric case is defined as the case where all the UEs are at the same distance d from the base station. For this case, we prove in theorem 10 that we can additionally reduce the computation complexity by limiting the number of policies Γ that should be considered in order characterize the stability region. Thus, we avoid the computation complexity that comes from the consideration of all the policies Γ ∈ Ω Γ (i.e. for a network of K + U communications, it exists|Ω Γ | = (K + U )! priority policies overall). Therefore, we propose an -approximation of the exact stability region that highly reduces the complexity while guaranteeing a high precision.
Theorem 10. For the symmetric case of the multi-UE cellular scenario with a given couple {p s , p d }, an -approximation of the stability region (for a given ) is the set of λ ∈R c such that:     Hence, the exact stability region can be bounded as it follows: where Ω K0 Γ is the set of the feasible priority policies of the subsets of K 0 communications among all the K + U communications (with |Ω K0 Γ | = (K+U )! (K+U −K0)! ), S K0 α is the set of the values α i ∈ {0, α * i } where i corresponds to the UE2UE communications within these subsets of K 0 elements. The elements of µ (Γ) which are µ i,s (Γ) and µ j,u (Γ) (for 1 ≤ i ≤ K, 1 ≤ j ≤ U ) are respectively given by (9) and (10). Note that the value of K 0 is limited to (K + U ) which corresponds to the total number of communications.
Proof. See Appendix-G.
The proposed approximation in theorem 10 is defined by considering all the following cases: (i) the priority policies of all the subset of K 0 elements within the K + U communications, (ii) for each priority policies, the UE2UE communications within these K 0 elements have a fraction element that varies within its two border values (0 and α * i ). The value of K 0 depends on the precision of the approximation; it increases when decreases. Even for a tight value of , K 0 stays too much lower than K + U which reduces the computational complexity to: 2 K0 (K+U )! (K+U −K0)! . We can see that K 0 is constant for a given triplet {p s , p d , } and does not depend on the number of communications. Therefore, the importance of this approximation is to reduce the exponential complexity O 2 K for characterizing the exact stability region to a polynomial complexity O K K0 .

V. NUMERICAL RESULTS
Simulations are performed in order to validate our theoretical results and numerically validate the trade-off between the precision and the complexity of the stability region computation. We consider a cell of radius R = 500 m. We consider a 10 MHz LTE-like TDD system. When a user is scheduled, it transmits over 50 subcarriers called Resource Blocks (RB). The considered bit rates per RB are {r 1 = 400, r 2 = 200, r 3 = 0} kbps/RB. The unit of the shown results is the arrival rate per RB [kbps/RB] from which the total bit rate can be deduced by multiplying by the number of allocated RB. The total transmission power (i.e. over all RBs) of the BS is P DL = 40 W and that of the mobiles is P U L = 0.25 W (see [29]). The noise power is equal to −195 dB/Hz for the downlink DL and −199 dB/Hz for the uplink U L. . We consider the pathloss model specified in [29]. The SNR thresholds are respectively γ U L 1 = 9.5 dB, γ U L 2 = 2.5 dB, γ DL 1 = 7.5 dB and γ DL 2 = 1.5 dB. These values are practical values based on a throughput-SNR mapping for a 10 MHz E-UTRA TDD network (see [30]).

A. Three-UEs scenario
We start by presenting the results of three-UEs network in order to validate the theoretical results corresponding to this scenario. For clarity of the presentation, we consider that the three UEs UE s , UE u and UE d are at the same distance d from the BS. This assumption is just for simplifying the illustration of the results however the theoretical results can be applied for any distribution of the users in the cell.
1) Approximated stability region: We illustrate the approximated stability region of the three-UEs scenario for different distances d between the users and the BS. This study is done for two purposes: (i) validating that the computed set of α in equation (8) as well as the considered set of priority policies Ω ss Γ are sufficient for characterizing the approximated stability region and no need to consider all the values α ∈ [0, 1] and all the priority policies Γ ∈ Ω Γ and (ii) elaborating the effect of the distance d on the performance of the users.
In Fig. 6, for different distances d between UEs and BS d = {100, 350, 500}, we plot in Red the stability region for the approximated three-UEs scenario obtained by exhaustive search (all Γ ∈ Ω Γ and all α ∈ [0, 1] 4 ) and we compare it to the one plotted in Green and obtained from Theorem 4. We find that both curves coincides which verifies that the latter theorem presents the stability region for the approximated three-UEs scenario. The importance of this theorem is the fact of providing a simple and close form expression of this stability region which is a close approximation of the computational complex exact stability region (as proved in theorem 5and numerically showed in the next subsection V-A2). Fig. 6 gives also insights on the effect of the distance d on the three-UEs stability region and shows how the performance of the system degrades when the distance d increases.
2) Exact vs. approximated stability region : In this section we numerically validate thatR ss c (given by (4)) is a close approximation of the stability region for the three-UEs scenario R ss c (given by (2). We present in Fig.7 that the stability region and its approximation coincide for different distances d between UEs and BS. Simulations considering different users positions verify that the approximated stability region is a good approximation for the exact stability for the three-UEs scenario.
3) Queuing impact: We start by showing the advantages of the queuing theory approach in our analysis. For this purpose, we consider a random positioning of 3 UEs in the cell and we apply the two performance evaluation approaches on the cellular scenario: (i) taking into account the queuing aspects and the coupling between the queues (ii) without considering any coupling between the queues. For the first approach, we compute the stability region as the set of the arrival rate vectors to the sources that can be stably supported by the network considering all the possible policies, hence we consider the coupling between the queues (w/ coupling). For the second approach, we assume that the BS queues have a full buffer and we compute the rate region that describes the achievable data rates depending on the the channel states of  the links and without any bursty traffic neither coupling between the queues (w/o coupling). Comparing these two results in Fig. 8 verifies that introducing the traffic pattern and the queues' coupling have an effect on the performance of the cellular scenario in terms of stability region. Considering that the queues have a full buffer, the DL (BS to UE 3 ) has always packets to transmit and it will always be considered by the scheduler. This causes less scheduling for the UL thus a lower bit rate over the UL. However, when coupling between the queues is considered, the DL can be scheduled only if at least one packet is received by the BS from UE 1 . Due to that, when no packet are buffered at the DL, this latter will not be scheduled and the UL will be able to transmit more packets and will improve its bit rate. This explains the gain that the queuing approach offers to the performance evaluation of the scenario.

B. Multi-UE scenario
For the multi-UE scenario, without loss of generality and for clarity reasons, we illustrate the results for the symmetric case where all the users are at the same distance d from the BS.

1) Exact vs approximated stability region:
For the multi-UE scenario, we consider a cell containing 50 UE2UE communications with a uniform random drop in a cell of radius R = 500m such that all the UEs are at distance d = 350m. The performance metric that we use to compare between the stability regions of the exact symmetric case (expressed in theorem 7) and the approximated symmetric case (expressed in theorem 10) is the average service rate per user. It is equal to the sum of the service rates of all the users divided by the number of users. For d = 350m and = 0.1 we find K 0 = 3. Fig. 9 illustrates that the variation of this performance metric as function of the number K of UE2UE communications in the network. It shows that the -approximated stability region is a very tight approximation of the stability region. We deduce that the complexity is reduced from 2 50 × 50! 3 × 10 79 for the exact stability region (from theorem 7) to 2 3 × 50! 47! = 9 × 10 5 for the approximated stability region (from theorem 10). This illustrates how the exponential complexity O 2 50 is reduced to a polynomial one O 50 3 .
2) Trade-off complexity vs precision: For the symmetric case, where all the users are considered at the same distance d from the BS, we show in theorem 10 that the value of K 0 characterizes the number of priority policies that is sufficient to be studied in order to describe the stability region, hence the complexity of the approximated stability region. For three different ∈ {10 −3 , 10 −2 , 10 −1 }, we show in Fig.10 the variation of K 0 as function of the distance d between the users and the BS. We can see that for = 10 −2 , the maximum value of K 0 (achieved at the edge of the cell) is small and equal to 12. Thus, the complexity of the approximated stability region can be deduced by applying the following: In this paper, we have carried out a queuing analysis to study a TDD network involving two types of communications: 1-UE2UE communications passing through the BS and 2-UE2BS communications. TDMA user scheduling and link adaptation model are the main assumptions in this work. First, we have shown the interest of our queuing theory approach by verifying that it provides a more realistic performance evaluation of the network compared to a performance analysis strictly based on physical layer that ignores the dynamic effect of the traffic pattern. This is due to the coupling between the UEs' queues and the BS queues. Second, we have given the analytic expressions of the exact stability region for the general scenario of different UE2UE and UE2BS communications. Then, we have studied the complexity of such a computation. Finally, a highly accurate approximation is proposed such that simple analytic expressions of the stability regions are provided as convex polytopes with limited number of vertices (computationally feasible). The variation of these stability regions as a function of the queues and channels states is investigated. We examine the trade-off between the precision of the proposed approximations and their computational complexity.

A. Proof of theorem 2
Here, we find the exact stability region for the three-UEs scenario. For this aim, we proceed as it follows: step 1 models the queue Q BS as a Markov chain and expresses the arrival and service probabilities of this chain, step 2 computes the probability that the queue Q BS is empty, step 3 obtains the service rate of both queues Q s and Q u , step 4 combines the results of the previous steps to provide the exact stability region of the three-UEs scenario, step 5 verifies that having a set of arrival rates within this stability region is equivalent to the stability of the system of queues. The main challenge is to solve the complicated Markov Chain that models the queue Q BS in order to find the probability that this queue is empty and deduce the service rate of the queues in the systems.
The coupling between the queues leads to a multidimensional Markov Chain modeling of the system of queues. Since this model complicates the study of the stability region, we approach this problem in a different way. We use the priority policies definition (in section II) to transform the multidimensional model to a 1D Markov Chain (one dimensional). Hence, for a given priority policy, each queue can be modeled by a 1D Markov Chain. In order to characterize the stability region, all the priority policies should be considered or at least the priority policies that achieve the corner points of the stability region.
We recall the three-UEs scenario contains one UE2UE communication and one UE2BS communication. The UE2UE communication is modeled by the cascade of Q s and Q BS and the UE2BS communication is modeled by the queue Q u . We consider the set of three rates {r 1 , r 2 , r 3 } with r 1 = kr 2 and r 3 = 0. So, if P is the number of packets transmitted at r 2 then kP is the number of packets transmitted at r 1 . Assuming that the two sources queues Q s and Q u are saturated, we want to characterize the stability region.
Step 1: Markov chain model of Q BS Cellular communications are modeled as coupled processor sharing queues where the service rates of Q s (equivalent to the arrival rate of Q BS ) as well as the service rate of Q u depend on the state (empty or not) of Q BS . Let us study the Markov chain of the queue Q BS in order to deduce the probability of being empty which means the probability of having at least one packet.
Q BS corresponds to the queue at the BS side that can be modeled as a Markov chain with a transition probability from a state x i to a state x i+1 is equal to the probability of receiving P packets per slot which is equivalent to the transmission probability of the queue Q s at rate r 2 due to the fact that the packets arriving to Q BS correspond to the packets transmitted by UE s . Thus, the transition probability from a state x i to a state x i+k is equal to the probability of receiving kP packets per slot which is equivalent to the transmission probability of the queue Q s at rate r 1 . Moreover, the transition probability from a state x i to a state x i−1 is equal to the transmission probability at the downlink BS-UE d at rate r 2 and the transition probability from a state x i to a state x i−k is equal to the transmission probability at the downlink BS-UE d at rate r 1 . As we mentioned before, the scheduling decision depend on the state (empty or not) of Q BS which means that the behavior of Q s and Q u is coupled to the state of Q BS and the arrival probability when Q BS is empty (state 0) is not equal to the arrival probability when Q BS is not empty (state = 0). In figure 11, we present the discrete time Markov Chain with infinite states that describes the evolution of the queue Q BS . We start by studying the queue Q BS in terms of probability of transmitting and probability of receiving packets in order to deduce the probability that the queue Q BS is empty.
Step 1.a: Service probability of Q BS For Q BS , two service probabilities exist the first one b 11 corresponds to the probability of transmitting at rate r 1 and the second one b 12 corresponds to the probability of transmitting at rate r 2 .
Parameters U, V, W , X , Y , Z , M and N from table III are used to make expressions simpler. These parameters depend on the considered priority policy Γ and their values are deduced from that of the transmission probabilities (such that those in table II for Γ 3 ). The values of these parameters are specified in table I for all Γ ∈ Ω ss Γ . Hence, the service probabilities b 11 and b 12 of the queue Q BS can be written as it follows: For each priority policy Γ ∈ Ω ss Γ , we find the transmission probabilities of all the links (s, d and u) then we compute the parameters U, V, W , X , Y , Z , M and N using table III.
Step 1.b: Arrival probability of Q BS Due to the coupling between the queues Q s and Q u with the state (empty or not) of Q BS , we should differ between (i) a 01 and a 02 that denote the arrival probabilities receptively at rate r 1 and r 2 when Q BS is empty and (ii) a 11 and a 12 that denote the arrival probabilities receptively at rate r 1 and r 2 when Q BS is not empty. For a given priority policy Γ, these probabilities are given by: p i u qs(S1, Si, S1, Γ |QBS(t) > 0) = α1U Step 2: Probability that Q BS empty Let us consider a 0 = a 01 + a 02 , a 1 = a 11 + a 12 and b 1 = b 11 + b 12 . We solve the balance equations in order to find the stationary distribution Π of the Markov chain of Q BS and more specifically the probability that Q BS is empty Π 0 . Recall that the stationary distribution of a Markov Chain exists if and only if the stability condition is verified. Hence, finding the stationary distribution based on the balance equations is a necessary and sufficient condition for guaranteeing the stability of Q BS . As shown in figure (11), the balance equations at the different states i of the Markov chain are the following: For i > k : We look for a solution of the balance equations of the form Π n = x n . Then we construct a linear combination of these solutions (which also satisfying the boundary equations (17), (18), (19)) that verifies the normalization equation n Π n = 1.
We start by substituting Π n by x n in (21) and then dividing by x n−k . This yields the following polynomial equation: The first root x = 1 but this one is not a useful, since we must be able to normalize the solution of the equilibrium equations. Stability conditions requires that this polynomial has at least one root x with |x| < 1. Let us say we found R roots such that |x| < 1: . Then the stationary distribution is given by the following linear combination: Where Π n for 0 ≤ n ≤ k as well as the coefficients c i for 0 ≤ i ≤ R are computed by solving the following system of equations: (i) first R + K balance equations and (ii) the normalization equation: For each choice of the coefficients c k the linear combination satisfies (21). These coefficients add some freedom that can be used to also satisfy the initial balance equations for 0 ≤ i ≤ k as well as the normalization equation. Therefore, we have R+k+1 unknowns to determine in order to find the stationary distribution of the queue Q BS , these unknowns are: the R coefficients c i and the k + 1 probabilities Π i for 0 ≤ i ≤ k. These unknowns are found by solving the following linear system AX = B with X T = [Π 0 , Π 1 , Π 2 , ..., Π k−1 , Π k−1 , c 1 , c 2 , c 3 ..., c R−1 , c R ] (Note that X T is the transpose of the vector X). By substituting of (23) into the first R + k − 1 balance equation as well as the normalization equation, we deduce for k > 2 the matrix A and B of the linear system as it follows. (The special case of k = 2 is treated afterward. The result above holds for any k ∈ N + . However, we explicit the result for the case (r 1 = 2r 2 ⇒ k = 2) where k − 1 = 1 hence the case 2 ≤ i ≤ k − 1 should not be taken into the account for this especial case. Thus, the set of balance equations is the following: For i = 1: For i = 2: For i ≥ 3: By analogy to the general case of r 1 = kr 2 , we look for a solution of the balance equations of the form Π BS (n) = x n and then we construct a linear combination of these solutions also satisfying the boundary equations (24), (25), (26) as well as the normalization equation n Π BS (n) = 1. Substituting of Π BS (n) by x n into (27) and then dividing by x n−2 yield the following polynomial equation: Recall that the root x = 1 is not a useful one, since we must be able to normalize the solution of the equilibrium equations. Stability conditions requires that this polynomial has at least one root x with |x| < 1 say x 1 , x 2 , x 3 . We now consider the stationary distribution given by the following linear combination: and Π n for = n = 0, 1, 2 are given by the boundary equations. For each choice of the coefficients c k the linear combination satisfies (27). These coefficients add some freedom that can be used to also satisfy the equations (24), (25), (26)  It can be shown that the unknowns Π 0 , Π 1 and Π 2 , c 1 , c 2 , c 3 is found by solving the following linear problem AX = B such that: Solving the system of linear equations given by AX = B gives the values of Π i for 0 ≤ i ≤ k − 1 and of the coefficients c j for 1 ≤ j ≤ R. Therefore, we find the probability that the queue BS is empty which is Π 0 .
Step 3: Service rates of Q s and Q u We follow the procedure below, based on queuing theory analysis of the network capacity, in order to derive the stability region study of the cellular scenario. For the 3-UEs scenario, the stability region is characterized by computing the two service rates µ s and µ u of the UE2UE and UE2BS communications.
Step 3.a: Service rate of Q s Let us start with the service rate of the UE2UE communications µ s . If Q BS is empty then the service rate of Q s is denoted by µ 0 s , otherwise it is denoted by µ 1 s . Hence, the average service rate of Q s for a given priority policy Γ is computed as it follows: with µ 0 s (Γ) and µ 1 s (α, Γ) given by: In order to simplify expressions, we consider the following notation U and V that depend on the considered priority policy Γ and their values for the six priority policies Γ ∈ Ω ss Γ are specified in table I. Hence µ 0 s and µ 1 s can be written as it follows: For each priority policy Γ ∈ Ω ss Γ , we compute the probability that Q BS is empty Π 0 (α, Γ) from step 2 then we substitute the above expressions of µ 0 s (Γ) and µ 1 s (α, Γ) into the equation (30), we obtain (3) as the service rate µ s (α, Γ) of the queue Q s .
Step 3.b: Service rate of Q u After computing µ s , let us compute the second element of the stability region which is the service rate µ u of the UE2UE communication. If Q BS is empty then the service rate of Q u is denoted by µ 0 u otherwise it is denoted by µ 1 u . By analogy to the Q s , the service rate of Q u is computed as it follows: with µ 0 u (Γ) and µ 1 u (α, Γ) given by: Parameters W, X, Y and Z from table III are used to make expressions simpler. These parameters depend on the considered priority policy Γ and their values for the six priority policies Γ ∈ Ω ss Γ are specified in table I. Hence, µ 0 u and µ 1 u of the queue Q u can be written as it follows: For each priority policy Γ ∈ Ω ss Γ , we compute the probability that Q BS is empty Π 0 (α, Γ) from step 2 then we substitute the above expressions of µ 0 u (Γ) and µ 1 u (α, Γ) into the equation (33), we obtain (4) as the service rate µ u (α, Γ) of the queue Q u .
Step 4: Characterization of the stability region Combining the results of the previous steps, we deduce the exact stability region of the three-UEs scenario. Supposing that the arrival and service processes of Q s and Q u are strictly stationary and ergodic then their stability which is determined using Loyne's criterion is given by the condition that the average arrival rate is smaller than the average service rate.
The following procedure is pursue for capturing the stability region of the scenario. We start by considering a priority policies that correspond to the corner points of the stability region (Γ ∈ Ω ss Γ with |Ω ss Γ | = 6 ). Then, for this priority policy, we vary α ∈ [0, 1]. For each value of α, we find the probability that Q BS is empty in order to deduce the service rates of the queues in the system. This procedure is applied for all the priority policies Γ ∈ Ω ss Γ .
Therefore, the stability region for the 3-UEs cellular scenario is characterized by the set of mean arrival rates λ s and λ u inR ss c such that : where µ s (α, Γ) and µ u (α, Γ) are respectively given by (3) and (4).
Step 5: Arrival rates within the stability region is equivalent to the stability of the system of queues We prove that if the mean arrival rates λ s and λ u are in the region R ss c is equivalent to the stability of the system of queues. To do so we prove that being having λ ∈ R ss c gives the stability of the queues and vice versa. We suppose the following notation: • a i (t): arrival process at queue Q i for i = s, d, u • λ s and λ u the mean arrival rates at respectively the UE2UE and UE2BS communications. By the law of large numbers, we have with probability 1: The second moments of the arrival processes E a i (t) 2 are assumed to be finite.
where b i (t) represents the amount of service offered to channel i on slot t and is defined by a function b i (S (t) , Q (t) , Γ (t)): Further, by the law of large numbers, we have with probability 1: Step 5.a: λ ∈ R ss c ⇒ Stability of the queues for different scheduling policies Γ such that λ i ≤ µ * i . We have by the law of large numbers that: which gives that queue i is stable for i = s, d, u. We deduce that if λ ∈ R ss c then the system of queues is stable.
Step 5.b: Stability of the queues ⇒ λ ∈ R ss c If the queues are stable then each queue i ∈ {s, d, u} has a stationary distribution Π i . The mean service rate is given by: Knowing that each queue i ∈ {s, d, u} is stable then is a combination of the limit of the stability region R ss c hence λ ∈ R ss c .

B. Proof of lemma 3
Here, we propose an approximated model for the three-UEs scenario and we define the condition that the fraction vector α should verify in order to satisfy the stability of the approximated BS queuẽ Q BS . For this aim, we proceed as it follows: step 1 proposes the approximated Markov chain model of Q BS and gives the expression of the probability that this queue is empty as function of the arrival and service probabilities, step 2 provides the expressions of the arrival and service probabilities of this chain, step 3 deduces the stability condition that the corresponding fraction vectorα should verify. The main challenge is to find an approximated Markov Chain that satisfies the following criteria: simple model, with an explicit form of the stationary distribution, that reduces complexity and finally an approximation that provides close results to the exact stability region.
Step 1: Approximated system model Cellular communications are modeled as coupled processor sharing queues where the service rates of Q s (equivalent to the arrival rate of Q BS ) as well as the service rate of Q u depend on the state (empty or not) of Q BS . We have already studied the stability region of the cellular scenario by consider the real Markov chain representing the queue Q BS . This real Markov chain consider the transition probabilities of both rates r 1 and r 2 . We note that finding such stability region is done by varying the vector α ∈ [0, 1] 4 and solving for each α the linear system of equations (see subsection A) in order to find the probability that Q BS is empty. Let us present an approximation of the stability region for which we provide an explicit formula of the points (µ s , µ u ) as well as the exact fraction vectorα * for which the limit of the stability region is achieved. In other terms, this explicit formula avoid us the consideration of all the α ∈ [0, 1] 4 as well as the resolving of the system of equations.
To do so, we consider the average number of packets received and transmitted per time slot at the BS side in order to find the probability of receiving and transmitting a packet per time slot at the BS level. Thus, the goal is to have an approximated queueQ BS that can be modeled by a simple birth and death Markov chain. We remark first that the average number of packets received by the BS is given by P a 02 + kP a 01 when the queue is empty and this average is equal to P a 12 + kP a 11 otherwise, we note second that the average number of transmitted packets by the BS is equal to P b 12 + kP b 11 . In other terms, with probability a 12 + ka 11 the BS is receiving P packets per time slot and with probability b 12 + kb 11 the BS is transmitting P packets per time slot.
The approximation is based on the simplification of the Markov Chain model that helps us to reduce the complexity of the problem. Therefore, the approximated Markov Chain still take into account the multiple rate model at the level of the probabilities of transition. Indeed, the multiple rate model is integrated in such a way that the probability of transmitting P packets at rate r 1 is given by the sum of the two following terms: (i) k times the probability of having a transmission at a rate r 1 and (ii) the probability of having a transmission at a rate r 2 .
Step 1.a: Approximated Markov ChainQ BS As we mentioned before, the scheduling decision depend on the state (empty or not) of Q BS which means that the behavior of Q s and Q u is coupled to the state of Q BS . Thus, the arrival probability when Q BS is empty (state 0 of Markov chain) is not equal to the arrival probability when Q BS is not empty (state = 0 of Markov chain). The approximated queueQ BS corresponds to the queue at the BS side for which the transition probability from a state x i to a state x i+1 is equivalent to the probability of transmitting P packets per time slot by Q s which is equal to the probability of receiving P packets per time slot byQ BS which means a probability of a 02 + ka 01 whenQ BS is empty and a 12 + ka 11 otherwise. Moreover, the transition probability from a state x i to a state x i−1 is equivalent to the probability of transmitting P packets per time slot at the downlink BS-UE d which means a probability of b 12 + kb 11 . In figure 12, we present the discrete time Markov Chain with infinite states which describes the evolution of the approximated queueQ BS . We see that these probabilities depend on the priority policy Γ, but for clarity reasons of the expressions we may use the notation b BS , a 0 BS and a 1 BS instead of b BS (Γ), a 0 BS (Γ) and a 1 BS (α, Γ). For the approximated Markov chain we consider the fraction vectorα = [α 1 ,α 2 ,α 3 ,α 4 ].
Step 1.b: Probability thatQ BS is empty The probability thatQ BS is empty is deduced from the stationary distributionΠ of the approximated Markov Chain ofQ BS (figure 12) which is computed by applying the balance equations as follows: Knowing that the stationary distribution verifies ∞ n=1Π n = 1; hence the probability thatQ BS is empty for a given priority policy Γ is given by: Step 2: Service and arrival probabilities ofQ BS In order to find the stability conditions, we need to find the expressions of transitions probabilities of the approximated Markov chain ofQ BS .
Step 2.a: Service probability ofQ BS Using parameters U, V and N from table III and based on equations (11) and (12), the service probability b BS can be written as it follows: Step 2.b: Arrival probability ofQ BS We should differ between a 0 BS : arrival probability whenQ BS is empty and a 1 BS : arrival probability whenQ BS is not empty. Using parameters U, V and N from table III and based on equations (13), (14), (15) and (16), the probabilities of arrival a 0 BS and a 1 BS of the approximated queueQ BS respectively when it is empty or not empty can be written as it follows: The values of the parameter above U , V and N are given in table I for all the considered priority policies Γ that correspond to the corner point of the stability region (Γ ∈ Ω ss Γ ). Since the arrival probability a BS at the queueQ BS is given by: Step 3: Stability condition ofQ BS Supposing that the arrival and service processes ofQ BS are strictly stationary and ergodic then the stability ofQ BS can be determined using Loyne's criterion which states that the queue is stable if and only if the average arrival rate is strictly less than the average service rate. Then the stability condition that the fraction vectorα = (α 1 ,α 2 ,α 3 ,α 4 ) should verify in order to satisfy the stability of the queuẽ Q BS is given by: Since a 1 BS and b BS are given by (39) and (37) then the queueQ BS is stable if the fraction vectorα verifies (5).

C. Proof of theorem 4
Here, we study the approximated model for the three-UEs scenario and we derive its corresponding stability region. For this aim, we proceed similarly to the proof in appendix A. However, the first two steps for the approximated model (the modeling of the BS queue by an approximated Markov Chain as well as the computation of the probability that this queue is empty) are already proved in appendix B. Hence, based on these results, we pursue the procedure by the following: step 1 computes the service rate of both queues Q s and Q u , step 2 provides the optimal fraction vector that depends on the priority policy Γ ∈ Ω ss Γ and that achieves the corner point corresponding to each priority policy Γ, step 3 combines the previous results and provides the simple and explicit form of the approximated stability region and step 4, similarly to the step 5 of the appendix A, verifies that having a set of arrival rates within this stability region is equivalent to the stability of the approximated system of queues. Here, there are two main challenges: (i) finding the approximated Markov Chain model of the queue Q BS that will help us to achieve an explicit and simple form of the stability region and (ii) the computation of the optimal fraction vector that will reduce the the complexity of the problem.
Step 1: Service rate of Q s and Q u The optimization problem above is a linear-fractional programming (LFP) where the objective function is the ratio of two linear functions with a set of linear constraints. This LFP has a bounded and non-empty feasible region and can be transformed into an equivalent linear problem (LP) which the solution (solved by any LP solution method such that simplex algorithm) yields that of the original LFP problem.
The solution of the optimization problem above belongs to the scenarios where at least one of the two constraints (C 01 ,C 02 ) is reached which corresponds to different cases: (i) constraint C 02 is achieved ((α 1 ,α 2 ,α 3 ,α 4 ) ∈ {0, 1} 4 ) while the stability condition C 01 is verified without equality and (ii) the stability condition C 01 is achieved with (α 1 ,α 2 ,α 3 ,α 4 ) ∈ [0, 1] 4 . These cases will be studied in order to come up with the set, called S α , of potential solutions of the optimization problem (45) . Several cases will be detailed where from each one we can deduce an element of the set S α .
Case 1: For the first case, we start by supposing that at the optimal solution the constraint C 01 is verified but not achieved while the second constraint C 02 is reached. The optimum in this case is the combination ofα =α 1 ,α 2 ,α 3 ,α 4 ∈ {0, 1} 4 that maximizesμ s (α, Γ) while verifying the stability condition C 01 .
Case 2: For the second case, we suppose that the stability condition C 01 is achieved at the optimum. Hence, the optimization problem is converted to the problem (46) with F 1 = M and By analogy to (45), we solve the problem (46) by considering that at least one of the linear constraints is reached which corresponds to the two following cases: (i) the first one assumes that C 12 constraint is achieved and (ii) the other one considers that C 11 constraint is achieved. We study both cases then the optimum corresponds to the case that maximizesμ s (α, Γ).
We start by supposing that at the optimal solution the constraint C 11 is not achieved then the optimum in this case is the combination ofα 2 ,α 3 ∈ {0, 1} 2 that maximizesμ s (α, Γ) and for which it existsα 1 andα 4 ∈ [0, 1] that verify (46). max α 2 ,α 3 1 2 r 1 p 1 s U + r 2 p 2 s V 1 + (k−1)α 2 p 1 Case 3: For this case, we suppose that the constraint C 11 is active at the optimum by achieving either the lower bound F 1 (forα 1 =α 4 = 0) or the upper bound F 2 (for α 1 = α 4 = 1). For each bound F 1 and F 2 respectively,α 2 andα 3 are chosen in a way to verify respectively the following equations: For this case, we suppose that the constraint C 11 is active at the optimum then the optimization problem can be written as in (47).
The LFP problem (47) depends only fromα 2 where the objective function is linearly depending onα 2 with a positive coefficient. Hence, the objective function is maximized by increasing as much as possiblẽ α 2 such that the constraints (47) are verified.
Combining the results of the three cases above, for each priority policy Γ, we find the set S α to which corresponds the optimum α * that maximizesμ s (α, Γ) and achieves the corresponding corner point of the stability region. Hence, we hugely reduce the complexity by finding the explicit value of α * that gives the corner points that are sufficient to characterize the approximated stability region.
Step 3: Characterization of the approximated stability region Combining the results of the previous steps, we deduce the approximation of the stability region for the three-UEs scenario. Supposing that the arrival and service processes of Q s and Q u are strictly stationary and ergodic then the stability ofQ BS can be determined using Loyne's criterion which states that the queue is stable if and only if the average arrival rate is strictly less than the average service rate.
We proceed as it follows for characterizing the approximated stability region. We start by considering a priority policy that corresponds to a corner point of the region (Γ ∈ Ω ss Γ with |Ω ss Γ | = 6). Then, for the considered Γ, we find the corresponding optimum fraction vector α * within the finite set (8), . We deduce the service rate of the system of queues. Finally, these steps are applied over all the priority policies Γ ∈ Ω ss Γ . Thus, the approximated stability regionR ss c for the 3-UEs scenario is characterized by the set of the mean arrival rates (λ s , λ u ) inR ss c such that: whereμ s (α * , Γ) andμ u (α * , Γ) are respectively given by (6) and (7).
Step 4: Arrival rates within the stability region is equivalent to the stability of the system of queues Proving that having a set of mean arrival rates (λ s , λ u ) inR ss c is equivalent to the stability of the system queue is identical to the demonstration done in step 5 of the appendix A. For brevity, we avoid repetition and remove this proof.

D. Proof of theorem 5
Here, we compare between the real and approximated stability region for the three-UEs scenario. For this aim, we proceed similarly to the proof in appendix A. For a given priority policy Γ ∈ Ω ss Γ , we proceed as it follows : in step 1, we compare between the two Markov Chain model by computing the absolute error between the two values Π 0 andΠ 0 that respectively correspond to the probabilities that Q BS andQ BS are empty and step 2 provides the absolute error between the service rates µ s andμ s in order to veridy that the approximation is an upper bound of the stability region. In step 3 we deduce the distance between these two regions based on the maximum relative error between their six corner points that correspond to the set of the border priority policies Ω ss Γ . Finally, in step 4 we find the upper and the lower bounds of the stability region as function of its approximation.
Step 1: Comparison between the Markov chains of Q BS andQ BS We start by presenting the used notation for facilitating the readability of the expressions: • Let us keep in mind that the comparison between the two Markov chain models (exact and approximation) is done for each priority policy Γ ∈ Ω ss Γ , hence to improve the clarity of the expressions we use the following: a 01 = a 01 (α, Γ), a 02 = a 02 (α, Γ), a 11 = a 11 (α, Γ), a 12 = a 12 (α, Γ), b 11 = b 11 (α, Γ), b 12 = b 12 (α, Γ) • a 0 = a 01 + a 02 , a 1 = a 11 + a 12 , b 1 = b 11 + b 12 • Π n = Π n (α, Γ) probability that Q BS is at state n of the Markov chain which means that the queue contains n P-packets •Π n =Π n (α, Γ) probability thatQ BS is at state n of the approximated Markov chain which means that the approximated queueQ BS contains n P-packets.
The main challenges of the exact stability region computation lies on the following to facts: (i) this region has no closed form expression since the final solution requires solving a linear system and (ii) the complexity of the exact stability region is high. Therefore, the advantage of the proposed approximation is to have a close upper bound that has an simple and explicit solution.