Enhancing the performance of downlink NOMA relaying networks by RF energy harvesting and data buffering at relay

Recently, non-orthogonal multiple access (NOMA) and energy harvesting (EH) from radio frequency (RF) have been considered as promising candidates for next-generation mobile communications. In this paper, we investigate a novel solution to enhance the reliability and the supply stability of a downlink NOMA relaying networks, in which we integrate two techniques: (i) simultaneous wireless information and power transfer, i.e. the relay node can harvest the energy from source signals and use this energy to help forward information from source node to two user nodes; and (ii) data buffer aid at relay node, i.e. the data packets received from the source can be stored in a buffer and then be re-transmitted to the destination nodes only when the channel condition is good. The performance of the proposed system is analyzed rigorously to derive the system outage probability (OP), throughput and the average packet delay. Furthermore, a power allocation optimization problem to minimize the OP is formulated and solution to this problem is also provided. In addition, the optimal transmission rates to maximize the throughput of each user are presented in numerical results. Monte Carlo simulations are conducted to verify the analytical results, which confirms that with data buffer at relay, the overall outage probability (OOP) can be reduced significantly.


Introduction
The non-orthogonal multiple-access (NOMA) is considered as a promising multiuser communication technique for the fifth-generation (5G) mobile network since it can achieve superior spectral efficiency, especially for applications in large-scale Internet-of-things (IoT) networks [1]. Unlike the orthogonal-multiple-access (OMA), NOMA users can share the same radio resources, including time and bandwidth. The key idea for this advantage is to employ the power domain or codes domain, where different users are distinguished by their power levels or different code [2]. Recent researcher have shown that NOMA can be applied to not only point-to-point but also relay networks [3]. While the application of NOMA to the pointto-point networks were well investigated, there is still increasing need for the case of cooperative relaying networks [4][5][6]. The authors of [4,5] proposed a dual-hop cooperative relaying scheme using NOMA, where two source nodes communicate with each other simultaneously via a common relay on the same frequency band. Especially, in [5], the author considered replacing the static relay by unmanned aerial vehicles (UAV). In this scheme, after receiving symbols transmitted in parallel by both sources with different power levels, the relay or UAV forwards the superposition coded composite signal using NOMA to two destinations. However, in this work power control for uplink multiple access was not considered. In 2020, Hoang et al. [6] proposed a multi-beamforming scheme for multiple clusters in downlink NOMA system, where the relay nodes are the cluster heads.
In addition, harvesting energy from the ambient environment has become a promising solution for energyconstrained electronic devices, which are conventionally supported by limited power sources such as battery [7][8][9][10]. In some special applications, charging the battery is too expensive or even impossible, e.g. sensor network works under toxic environment and body area network. Moreover, some natural energy sources such as solar and wind, and radio frequency (RF) can be also utilized as effective sources for energy harvesting (EH). Compared with other kinds of energy sources, the RF EH [11], also known as wireless energy transfer, has some advantages. Since RF EH is an active energy supply method, it can provide more reliable energy flow to guarantee the quality of service. On the other hand, buffer-aided (BA) relays have been shown to improve wireless conditions through increased diversity, thus resulting in reduced outages and better throughput [12].
Therefore, utilizing RF-EH technique together with NOMA scheme helps prolong the lifetime and improves the spectral utilization efficiency of the energy-constrained multi-user wireless relaying networks. The NOMA systems combining with RF EH were investigated in [13][14][15][16][17][18]. A simultaneous wireless information and power transfer (SWIPT) of NOMA networks were considered in [13], where base station simultaneously communicates with the users, namely, relay user and far user, by using finite block length data burst. The optimization problem to obtain the optimal power allocation factors to NOMA user fairness performance are proposed in [14]. The authors in [14] stated that the sum rates can be maximized when equal individual power allocation factors. The authors of [15] investigated a NOMA system in which the NOMA user near the source can act as the EH relays to assist the far NOMA users in forwarding the information. In these works, the authors considered the users reception signals in two different time slots, however, the service was simultaneously provided to users in the NOMA system. The authors in [16] considered a NOMA scheme in cognitive radio networks, where the relay node with SWIPT assists the source to transmit the superimposed composite signal to two receivers. They derived the closed-form of outage probabilities (OPs) for 2 receivers and minimized the system energy efficiency (EE). In [17], minimizing energy consumption is the ultimate goal in a SWIPT-NOMA downlink system with multiple near users and one far user, where the authors applied a golden section Dinkelbach algorithm to solve the optimization problem. However, there was no outage analysis in this work. A pair of best near and best far users were selected in [18] for improving the OP in a SWIPT-NOMA system, where the near user served as a relay for the far user. The common point of these above works is that they do not consider using a data buffer at the relay. The sum-throughput maximization problem by beamforming design for relaying NOMA networks with SWIPT was investigated in [19]. In this system, the users are grouped into multiple clusters based on cluster-specific beamforming. The results of this abovementioned works shown that the performance of NOMA outperform OMA scheme.
An energy BA EH relay was applied in cooperative communication system to improve sytem performance in [20], but in this work the authors didn't consider data BA. On the other hand, the work in [21] proposed a hybrid NOMA/OMA system with BA relay selection. As a result, two BA opportunistic relay selection algorithms were proposed. The aim of that work is to improve the outage performance and sum-rate of the system. The work in [22] studied the conventional cooperative NOMA system with BA relaying. Herein the authors considered an adaptive transmission scheme in which different working modes are employed in different time slots. The authors in [23] proposed a priority-based max-link relay selection for data BA decode-and-forward (DF) cooperative networks. In this work, the authors derived analytical expressions of outage probability (OP) and bit error rate (BER) to evaluate the system performance. In addition to NOMA downlink system investigation, the hybrid NOMA/OMA uplink system with the help of a BA relay was also considered in [24].
So far, most of previous works related to cooperative communication protocols and NOMA technique in literature have proposed the superposition signal coding at the source. Meanwhile, the relay node, which has a fixed power, only decodes and forwards signals to destinations. However, due to the random nature of wireless channel, the amount of energy harvested at the relay is usually very small and variable. Hence, the power allocation at the relay can reduce the feedback energy and guarantee the performance fairness for all users. To the best of our knowledge, the analysis of the effectiveness of using energy and/or data buffer at the relay is still moderated. In 2019, our research group has initially introduced the outage analysis of a relay-assisted NOMA networks with two users, in which SWIPT and data buffer aid could be applied at the relay [25]. Another work on SWIPT-NOMA systems that exploits the energy and data buffer at relay was proposed in [26], where the authors introduced a protocol of adaptive allocating the transmit powers to NOMA users and the power splitting factor of EH to minimize the total energy consumption. However, that paper assumes the ideal data buffer, i.e. the size of the buffer is sufficiently large to prevent the overflow, and the derivation of the OOP expression has not been carried out in that work. A summary of previous works on SWIPT-NOMA systems is shown in Table 1, in which their advantages and limitations are also highlighted. Motivated by the above facts, in this paper, we extend our research work by analyzing a different cooperative DF relaying scheme where a source transmits information packets to the relay, while relay broadcasts the modulation superposition signals to two users. Here, R employs time-switching-based EH prior to the communication with destinations. Based on the channel gain from the relay node to the destination node, the relay performs one of three power allocation schemes, i.e. fixed, adaptive, and optimal power allocation, for two users. Furthermore, NOMA technology is investigated for the cooperative transmission in term of two scenarios, i.e. with and without buffer aid at the relay node.
The main contributions of this paper are summarized as below: -We proposed a novel downlink NOMA relay system applying SWIPT, where designated relay node is either equipped with buffer or not. To harvest energy, the relay uses time switching protocol (TS). The system performance is investigated in both cases, i.e. perfect and imperfect SIC. In addition, the outage performance is considered for both linear and non-linear EH model. The performance of the system is improved, and additionally the spectral utilization efficiency and the lifetime of wireless networks will be enhanced. -The optimal power allocation at the relay node to minimize OP is also considered in this paper. Since the harvested energy is very small, the reallocation of the harvested energy after converting from the source RF signals is important for saving cost. Moreover, the optimal transmission rate to maximize the throughput of each user is also provided. -Markov chain model and state-transition matrix is used to describe the random process at the BA relay. It is shown that with BA relay the diversity gain is improved significantly.
-The system performance is demonstrated by OP and the throughput of the system over Rayleigh fading channel. We derive the closed-form expression to evaluate the rate of symbols and outage performance of SWIPT-NOMA system. The performance comparison between NOMA and OMA system for the case with/without BA relay is also provided. -The analytical results are validated by simulation. From the derived closed-form expressions, practical networks can be investigated in the in future.
Our proposed system can find applications in surveillance sensor networks for disaster detection or Internet of Things (IoT) networks, where installing fixed power lines or frequent replacing of batteries for a large number of nodes is inconvenient, for example, the sensor network that works in a toxic environment. Furthermore, the proposed system can be applied to sensor nodes in wireless body area networks for healthcare and other medical applications. Another application area that our work can be applied is RFID, which has caught intensive research and has been widely used for identification, tracking, and inventory management. Besides, the advantages of the proposed solutions in terms of low energy cost, reducing greenhouse effect, and prolonging lifetime are very useful for future mobile networks. The remaining of the paper is organized as follows. Sect. 2 presents the SWIPT-NOMA system model and channel model. The analysis of OP with and without buffer aid at the relay node is given in Sects. 4 and 4.2, respectively. The derivation of the system throughput is given in Sect. 4.3. The average packet delay is demonstrated in Sect. 4.5. Numerical results, which verify our analysis, are presented in Sect. 5. Finally, the conclusion is given in Sect. 6. For the convenience, we provide in Table 2 the notations along with their descriptions used in this paper.

System model
The NOMA downlink relay network that is investigated in this paper is shown in Fig. 1. In this network, a source ðSÞ wants to send its messages to two destinations ðD 1 Þ (the near user) and ðD 2 Þ (the far user) simultaneously with the help of a relay, ðRÞ, which is capable of harvesting energy from the source. It is assumed that S, D 1 , and D 2 have stable power supply, while R has no extra energy supply for forwarding task and needs to harvest energy from S. In addition, the relay node is assumed to have an information buffer to store the received messages [28]. We assume that the direct links between the source and both destinations are not available due to far distance or obstacles. All nodes are equipped with a single antenna and operate in a halfduplex mode. 1 The channels between two arbitrary nodes are subject to block and flat Rayleigh fading, which has only the smallscale component. This means that the channel coefficients are constant during each data transmission block, but vary from one block to another. We assume that the relay has perfect channel state information (CSI) of the S ! R link and R ! D 1 ; D 2 links at the beginning of each transmission block by applying some short pilot signals. Based on this set of information, R can decide whether it is ready to operate in transmitting or receiving mode [29].
As shown in Fig. 1, the complex channel fading coefficient of the link between S and R is denoted by The additive white Gaussian noise (AWGN) at R, D 1 and D 2 is respectively denoted by w A $ CN 0; N 0 ð Þ, where A 2 R; D 1 ; D 2 f g . Without loss of generality, we assume that the users' channel gains are sorted in the descending order as follows:

Energy harvesting and information transfer protocol
The time-switching relaying (TSR) protocol for EH [30] is applied in this work. Please refer to [30] for more detailed explanation. 2 Note that, the work in [30] only considered EH for the orthogonal multiple access system, i.e., the conventional relay system, and without data buffer aid. According to this protocol (shown in Fig. 2), each transmission block of duration T is divided into three time slots. In the first time slot of duration aT, the relay harvests the energy from the source. After that, the communication from the source (S) to the destinations (D 1 and D 2 ) consists of two time slots with equal length ð1 À aÞT=2, in which S can send the superposition signal to R and then R can forward 3 it to D 1 ; D 2 , respectively [31]. Here, 0 a 1 denotes the EH time-switching factor in one transmission block. Therefore, the amount of harvested energy at the relay for the case of linear model in the ith time slot is given by [32,33] where 0\g\1 denotes the energy conversion efficiency, which depends on the harvesting electric circuitry. We assume that all energy harvested is consumed by the relay for forwarding information signals to all users D i . 4 The processing power for the transmitting/receiving circuitry at the relay is generally negligible compared to the power used for information transfer and perhaps venial. So, from (1), the transmit power of the relay is given as Remark 1 In practice, the energy harvester will output a constant due to its circuit design. Such a non-linear model can represent the joint effect of non-linear phenomena caused by hardware constraints, including circuit sensitivity limitation and current leakage. The main reason for nonlinear EH model can be explained from the relationship between the input RF power and the output direct current of energy harvester. This nonlinear function is caused by nonlinear semiconductor devices such as diodes and transistors in the energy harvester circuit. Denote P th as the saturation power threshold of harvester, thus the transmit power of the relay node for the proposed nonlinear model using in this manuscript is given by

Data buffer at relay
In order to process the data the relay is equipped with a buffer for storing the decoded messages (a.k.a. packet) received from the source. For TS scheme the relay also has an energy storage device to store the harvested energy. 5 The communication procedure using data buffer at the relay can be described as follows. For each transmission block (with duration of T), the relay or source node is allowed to transmit data (in their assigned time slot as illustrated in Fig. 2) depending on the status of the relay buffers and the availability of the corresponding channel, so that it can provide a successful transmission or reception of one packet. 6 For example, if the S ! R link is in outage or the data buffer is full at one transmission block, then S should not send its message in the second time slot of that block.
If S is allowed to transmit, it will generate a transmit packet of size 2r 0 T bits, which is intended for two destination nodes D 1 and D 2 , to send to the relay. Here, r 0 is the desired transmission rate of the system. Each packet contains two segments, the first of which is used for the transmit symbols to D 1 and the other is used for the transmit symbols to D 2 . The relay buffer has L ! 2 storage units, each of which can store 2r 0 Tð1 À aÞ bits. The relay node decodes the received packet and stores it into the storage device. Each storage device is split into two parts of the same length, which are used to store the information symbols intended for D 1 and D 2 , respectively.
Assume that the source always has information symbols to transmit. At the relay, the number of storage units that have stored information symbols is denoted as l, where 0 l L. The availability of the links to transmit data depends on the value of l. In particular, for l ¼ 0, only the source-to-relay link is available because the relay has no information symbol to transmit. For 0\l\L, both the source-to-relay and relay-to-destinations links are available. For l ¼ L, only the relay-to-destinations links are available because the relay has no free storage units to store the information symbols received from the source.

Signal model
In each time slot, if the source (S) is selected to transmit data, it combines two messages x 1 and x 2 into a transmission packet by the superposed modulation in the different power levels. Then, the received signal at the relay is given by where P S denotes the transmit power of the source node, w R is AWGN at R, h 1 is the channel gain between source and relay and is modeled by Rayleigh distribution. So, jh 1 j 2 is exponentially distributed.
The signal-to-noise ratio (SNR) of the source-to-relay link is given by [28] In a specific time slot, if the relay (R) is selected, it transmits a modulated superposition information symbol x 2 stored in the buffer through multicast communication, where x 1 and x 2 denote the information symbols intended for D 1 and D 2 , respectively. a i ; i ¼ 1; 2 is the power allocation coefficient for D i , respectively. At the end of a time slot, the received signals at the destinations are given by Without loss of generality, it is assumed that the channel power gains are already sorted in descending order, i.e. jg 1 j 2 [ jg 2 j 2 . According to NOMA principle, the relay allocates more power for D 2 in order to balance the fairness of the system performance and to achieve the boundary of the capacity region for a pair of users [2]. Due to the broadcast nature of the wireless environment, the instantaneous signal-to-interference-and-noise ratio (SINR) of the R ! D 2 link is given by where the information symbol x 1 is treated as the interference at D 2 , i.e., SIC is not applied at D 2 . 7 At D 1 , the goal is to decode information symbol x 1 . By applying SIC technology [34]., D 1 can remove the detected information symbol x 2 from the received signals before decoding x 1 . From (5), the instantaneous SNR and SINR of the R ! D 1 link in the case of perfect SIC are expressed respectively as For the case of imperfect SIC, i.e., the residual interference power of x 2 has an impact on x 1 , the SINR at D 1 is given as where 0 j 1 the level of residual interference power of x 2 after SIC.  6 One packet consist of the decoded messages for x 1 and x 2 . 7 The order of SIC decoding merely depends on the order of SNR or equivalently, the magnitude of the channel gains, at D 1 and D 2 .

Problem definitions
To help the readers understand the multiple goals of this work, in this section we summarize the technical problems that we are going to tackle in this paper.

Overall outage probability analysis
The overall outage probability (OOP) of the system is defined as the probability that neither the source-to-relay link nor the relay-to-both-destination links is unavailable for transmission to achieve the target predefined transmission rate. For simplicity, we assume that the target transmission rates from the source to the relay and the relay to the destinations are the same and equal to r 0 . Hence, the outage event occurs when the instantaneous end-to-end capacity 1Àa 2 log 2 ð1 þ c e2e Þ is less than r 0 . The factor of 1Àa 2 appears due to the two consecutive time slots for communication between the source and the destination. It is easy to see that OOP is equal to the probability that the output SNR, c e2e , falls below a certain threshold, In the Sect. 4, we will derive the closedform expressions of the OOP for both cases: with or without the data buffer-aid at the relay.

Throughput analysis
The other goal of this work is to derive the throughput of the NOMA system with and without BA relay. The throughput is derived from the previous expression of OP. When the codeword is long enough (i.e., its duration is larger than the channel coherence time), the ergodic capacity should be considered as a performance factor. In contrast, when the codeword is short (i.e., its duration is less than the channel coherence time), the throughput should be considered instead of the ergodic capacity [35,36].
Assuming that the symbols x 1 and x 2 are transmitted at rates of r 1 and r 2 , respectively, over wireless channels with ideal channel coding, then the throughput for transmission of x ' , where ' 2 f1; 2g can be determined as where c x ' e2e is the end-to-end SINR at D ' . Then the total throughput of the considered NOMA without buffer aided relay system is We are going to derive the throughput expressions for both the non-data-buffer and with-data-buffer cases in Sect. 4.

Optimization of power allocation
To obtain the minimum OOP and minimum OP D 1 in (18) and (20) we formulate the OP minimization problems P 1 and P 2 as and P 2 : min where a 1 is the power allocation coeffiecient for D 1 (the power allocation coeffiecient for D 2 is a 2 ¼ 1 À a 1 ). The condition a 1 \ 1 1þc th indicates that the OP is less than one.

Optimal data transmission
To determine the optimal data transmission rate that maximizes the throughput of each user D i in (34), we formulate the following optimization problem: ð14Þ Due to the complexity of the throughput functions, in this paper we will use numerical method to find the optimal data rate. The optimal solutions are given for both perfect SIC and imperfect SIC cases in Sect. 5.

Average packet delay of the buffer-aided relay system
One of the important performance metric for the data buffer is the average packet delay. In this paper, we also aim to derive the closed-form expression of the total average packet delay in our considered system. To determine this quantity, we consider a period of time that is long enough, and denote the throughputs of S and R by s S and s R , respectively. Let's also define Q S and Q R as the average queuing length (in number of time slots) of S and R, respectively. Based on [37, Eq. (7)], Q S is derived by how the packets are delivered by S thoroughly. Then the average packet delay at the source is defined as and the average packet delay at the relay is expressed as Finally, the total average packet delay of the system is given by 4 Performance analysis

Outage probability without buffer aided relay
In this section, we investigate the system performance of the SWIPT NOMA downlink relaying in terms of the overall OP, the OP for each destination and the throughput of the system.

Overall outage probability
The following theorem provides the exact closed-form expression of the OOP and the approximation of the OP of the SWIPT-NOMA downlink relaying system.

Theorem 1
The OOP of the system when the relay knows both g 1 and g 2 is given by where Ei(x) denotes the exponential integral function [27], , and the con- At high SNR regime, OOP can be approximated by where K 1 ðÁÞ is the first-order modified Bessel function of the second kind.
Proof Please refer to ''Appendix A''. When the condition a 1 \ 1 1þc th holds, to ensure that the overall outage does not occur, we need to allocate more power to D 2 . With the assumption that SIC for information symbol x 2 at the D 1 is perfect, if the relay knows the channel responses g 1 and g 2 , it can adjust the power allocation coefficient a 1 to balance the OP of the relay-to-destinations links. h It should be noted that if jg 1 j 2 [ jg 2 j 2 , then ð1 À a 1 ÞP R jg 1 j 2 . This remark is very important for analyzing the OP expressions in the next part.

Outage probability at the destinations
In this section, we derive a closed-form expression of the OP at each destination. When one destination is in outage, the other still can detect its desired information symbol.
The system may switch to the conventional OMA system. However, in this case, the system performance will be degraded because the power allocation at the relay have been fixed for D 1 and D 2 . Theorem 2 provides closed-form expressions of OPs at D 1 and D 2 , respectively.

Theorem 2
The OP at D 1 and D 2 are given respectively as where Q max ¼max n 1 a 1 /P S ; 1Àa À1;i 2 f1;2g, and the condition a 1 1 1þn i , i 2 f1;2g holds. Here, r 1 and r 2 are the target transmission rates at D 1 and D 2 , respectively.
For the case of imperfect SIC at D 1 , the OP is given by where D max ¼ max n 1 /P S ða 1 À n 1 jð1 À a 2 ÞÞ ; Proof To obtain the OP expression of D 1 and D 2 , we first analyze the instantaneous SINR and SNR at D 1 and D 2 . For details of the proof, please refer to ''Appendix B''. h According to Theorems 1 and 2, the power allocation coefficients for NOMA system must be selected properly to guarantee the fairness between D 1 and D 2 . On another aspect, the impact of the transmission rate and imperfect SIC coefficient to the system OOP may be significant, as will be shown later in numerical results. Therefore, appropriate selection of the EH time-switching factor and the tranmission rate is necessary.

Optimal power allocation to minimize the outage probability
In this section, we study the power allocation problem to minimize the OP of the proposed SWIPT-NOMA system.

Theorem 3
The optimal power allocation coefficient a Ã 1 to minimize the OP in the EH-NOMA system is given by Proof The problems P 1 and P 2 defined in (12) and (13) are equivalent to the problems of maximizing respectively. In addition, we assume that c th ¼ n 1 , which means that the data rate of D 1 is equal to the system data rate.
; s:t:0 a 1 1; a 1 \ 1 1 þ c th : ð24Þ ; s:t:0 a 1 1; We consider two cases of the objective functions of P 1 and P 2 under the following conditions: Then, the above problems can be solved similarly as those in [38, 4.1]. After some mathematical manipulations, we obtain the optimal power allocation coefficient a Ã 1 as shown in Theorem (3). It should be noticed that the power allocation for D 1 is considered, subject to the condition that the decoding of data for D 2 is successful. 8 h

Analysis of the outage probability with buffer aided relay
In this section, we investigate the OP of the proposed system, where the buffer aid is employed at the relay node. For convenience, we assume that the source node always has data to transmit. We also consider the number of transmitted symbols as the number of transmitted packets. The general operation of the system with buffer aid relay is depicted as follow. The system chooses a node to transmit (source or relay) in a given time slot. To perform this, the information of the outage states of the S ! R link and R ! D i link is required. Therefore, the system uses one bit for feedback information from the destination to the relay. This information helps R known if the link R ! D i is in outage or not. One bit, which is fed back from the relay to source, is used to determine if the source should be in transmit or silent mode. In transmit mode, the source transmits the packets to the relay, then the relay decodes the packets and stores the decoded packets in its buffer. After that, the relay transmits the packets to the destinations node. If the source is selected to transmit but the S ! R link is in outage, the source should remains silent, and the outage occurs. Similarly, if the relay is selected to transmit but the R ! D i links is in outage, the relay remains silent. Unlike the case of without buffer-aid, in this case the outage event is defined as the probability that the relay does not receive and transmit, i.e. it remains silent. To describe the state transition of the BA relay, we denote the outage events of S ! R link and R ! D i link by O SR and O RD , respectively. The probabilities that the links are not in outage are 8 Moreover, if the relay knows the channel gains g 1 and g 2 and the total power factor a 1 þ a 2 is equal to one, we can allocate power for D 1 and D 2 by adapting to channel gains, i.e. a 1 ¼ jg 2 j 2 =ðjg 1 j 2 þ jg 2 j 2 Þ, to ensure the fairness in outage performance. O SR ¼ 1 À O SR and O RD ¼ 1 À O RD , respectively. In addition, R decision scheme is described as in Table 3. In Table 3, 'l' and 'L' respectively represents the number of packets stored in the buffer and the buffer size at the relay node. The outage and non-outage links are indicated by '0' and '1', respectively.
To calculate the OP of the system, from the Table 3, we build a Markov chain. We start at the initial state l ¼ 0 (i.e. when the buffer is empty). If the S À R link is in outage, which means that the source does not transmit, then the buffer will be empty. In other words, the buffer state moves from l ¼ 0 to l ¼ 0 with probability of O SR (Case B in Table 3).
When the S À R link is not in outage, we consider two cases. The first case is when the R À D link is in outage (Case D). Consequently, the relay receives the signal, making the buffer state moves from l ¼ 0 to l ¼ 1 with probability of ð1 À O SR ÞO RD . The second case is when the link RD is not in outage (Case G). The relay receives the signal, making the buffer state moves from l ¼ 0 to l ¼ 1 with probability of ð1 À O SR Þð1 À O RD Þ. Combining these two cases, the buffer state moves from l ¼ 0 to l ¼ 1 with probability of 1 À O SR . Similarly, we can obtain the probability of moving to the other states. Finally, we get the Markov chain showing the state transitions as depicted in Fig. 3. The whole operation of data buffer can be summarized in the Algorithm 1.
From Table 3 and the presented Markov chain, the OP of the considered system is calculated as where Prfl ¼ 0g and Prfl ¼ Lg are the probabilities of the events that the buffer is empty and full, respectively. To derive the OP of the system in (26), we define a state transition matrix A with size of ðL þ 1Þ Â ðL þ 1Þ for the Markov chain as follows, where A ij denotes the element of the ith row and jth column of the matrix A, which refers to the probability of moving from state i at time t to state j at time t þ 1.
For the case of L ¼ 5, matrix A is expressed as We should note that matrix A is not symmetric because the states are not symmetric and the transition probabilities are not the same. Then, the stationary distribution p of the Markov chain can be expressed as where I is an identity matrix, B is an ðL þ 1Þ Â ðL þ 1Þ matrix with all elements equal to 1, b ¼ ð 1 1 . . . 1 Þ T and denote as Theorem 4 With the BA relaying, the OP of the proposed system becomes Proof To determine the state transition matrix A, we need to derive O SR and O RD . We assume that the minimum data transmission rate from S ! R is r 0 , then the OP of S ! R link is defined as According to SIC principle, if D 1 is able to remove x 2 from its received signal, the OP of the link from the relay to the destination nodes is given by After doing some algebra, we get where A ¼

Throughput analysis
In this section, we derive the throughput of the NOMA system with and without BA relay. According to the definitive formula of the throughput [(34) and (11)], the total throughput of the considered NOMA without buffer aided relay system is where c (34) indicates that very high code rate r ' of the signal in each S À D ' communication link results in lower probability of successful decoding and degraded system throughput. In contrast, low r ' results in high probability of successful decoding, but the system throughput may not be as high as expected because of the low transmission rate r ' . The overall throughput of the considered NOMA system without buffer aided relay can be rewritten as

Remark 2 Equation
Remark 3 When using buffer at the relay to store data packets, the throughput of the relay is considered as a key Wireless Networks (2022) 28:1857-1877 1867 performance factor. The throughput of the BA relay is defined as transmission rate r 0 multiplied with the probability that the relay transmits successfully. If R transmits successful, the buffer state changes from state s i at time t to state s iÀ1 at time t þ 1. Thus, the throughput of the BA relaying NOMA can be found as

Optimal data transmission rate
As can be seen in (34) and (36), the system throughput is a function of the data transmission rate. It is necessary to determine the optimal rate that maximizes the throughput of each user. From (10), when the transmit power is given, the throughput of D 1 is given as where n 1 ¼ 2 r 1 À 1.
Since the channel coefficients are exponential random variables, the throughput of D 1 in the case of perfect SIC can be rewritten as and the throughput of D 1 in the case of imperfect SIC is Similarly, we have the throughput of D 2 as where n 2 ¼ 2 r 2 À 1. From (38), (39) and (40), we can use one-dimensional search techniques to find the transmission rate that maximizes the throughput of the D 1 and D 2 , respectively. The results will be shown in Sect. 5.

Average packet delay of the buffer-aided relay system
In this section, the average packet delay of the system, including the delay at the source and at the relay, is analyzed. When we consider a period of time that is long enough, we can assume that the sum of the packets transmitted by the relay is equal to the sum of the packets received by it. Therefore, both throughputs of S and R, denoted by s S and s R respectively, can be calculated as Let's define Q S and Q R as the average queuing length (in number of time slots) of S and R, respectively. Based on [37, Eq. (7)], Q S is derived by how the packets are delivered by S thoroughly. The source transmitted at most one packet in each time slot, thus Q S is Then the average packet delay at the source is determined as The average queuing length of R can be given as Thus, the average packet delay at the relay is expressed as Finally, the total average packet delay of the system is given by For example, from Table 3, R can be chosen in the Case F and Case G as SNR ! 1, i.e. OP ! 0. We assume that the system starts at state s 1 , (i.e., the buffer is empty, Fig. 3 The diagram of the Markov chain of buffer states at the relay node l ¼ 0). Based on Case G, the relay receives two packets in two time slots, so the buffer moves to state s 3 , (i.e., l ¼ 2). After that, in Case F, the relay transmits two packets to D 1 and D 2 in one time slot, the buffer moves to state s 2 (i.e., l ¼ 1). Then, based on Case G, the relay receives a packet and the buffer moves again to state s 3 , (i.e., l ¼ 2). Then, the buffer repeats the this cycle. Thus, we have the state transition probabilities is Then by using (45), we obtain D R ¼ 3 1ÀOP . When SNR ! 1, i.e., OP ! 0, we have D R ¼ 3. So, the total average packet delay of the NOMA system is D ¼ 4. The mathematical result matches with numerical result.

Numerical results
In this section, detailed numerical results are provided to illustrate the impact of power allocation on the performance of SWIPT-NOMA system with data BA relay, in terms of the OP and the throughput. For comparison, we also provide the performance of the SWIPT-OMA system with the same parameters. Configurations and parameters of the system are explained as follows. In this system, we only consider small-scale Rayleigh fading and assume that the terminals are stationary, where D 1 is closer to the relay nodes than D 2 . Hence, we need to allocate more power to D 2 than D 1 to ensure the user fairness. We perform 10 Â 2 14 independent trials for each Monte-Carlo simulation.
In this paper we consider three typical power allocation schemes as follows. Firstly, the optimal power allocation coefficient is derived in Theorem (3). Second, we used an adaptive power allocation scheme for D 1 and D 2 , in which the allocation factors are adjusted corresponding to channel gains, i.e. a 1 ¼ Efjg 2 j 2 g=ðEfjg 1 j 2 g þ Efjg 2 j 2 gÞ. Thirdly, the fixed power allocation scheme is considered, where the allocation factor for D 1 is fixed at a 1 ¼ 0:3 and the one for D 2 is 1 À a 1 .
The EH time-switching factor is a ¼ 0:3 and the energy conversion efficiency is g ¼ 0:8. The system data rate r 1 ¼ 1bps=Hz while r 2 ¼ r 0 ¼ 0:5 bps/Hz. The saturation power threshold of harvester varies as P th 2 f15; 20; 25g dB. The buffer size is set to L ¼ 25. In all following figures, 'Ana' represents the analytical results, while 'Sim' denotes the simulation results. On the other hand, 'Opt PA', 'Adapt PA', and 'Fix PA' represents the optimal power allocation, adapive power allocation, and fixed power allocation, respectively. The obtained numerical results show that the optimal power allocation can increase the system performance and the NOMA scheme significantly improves the spectrum utilization. Figure 4 illustrates the OOP in terms of the average SNR in dB for two cases, i.e., with and without data buffer aid at the relay node. As observed from Fig. 4, the OOP of the system (including the outage events of both D 1 and D 2 ) in the case of optimal power allocation outperforms the case of fixed power allocation. From Fig 4, we can see that the benefit of the optimal power allocation in terms of OOP compared to fixed and adaptive power allocation schemes is not significant. This is because we have selected fixed power allocation factor a 1 ¼ 0:3, which is approximated to the optimal power allocation factor a Ã 1 (when a ¼ 0:3 and r ¼ 0:5 returns a Ã 1 ¼ 0:2709). This is also approximated to the adaptive power allocation factor. In order to achieve the fairness of the overall system performance, the transmitter needs to allocate power according to the channel gains of R ! D 1 and R ! D 2 . Furthermore, the approximation results calculated from (19) are very close to the exact results obtained from (18), especially in the high SNR regime. Therefore, we can use (19) to easily calculate the OOP of the system. We can also see in Fig. 4 that the diversity order of the system with data buffer aid is equal to 2. Meanwhile, for the case of without BA relaying, the diversity order of the system is equal to 1. Hence, employing data buffer at the relay leads to the redution of OP, but it trades off with the packet delay. We can also see that the analytical results agree well with the simulation results. Figure 5 plots the OOP and the OP of D 1 and D 2 , respectively. The optimal power allocation coefficient as presented in Theorem 3 is used for both the cases with and without data buffer-aid at the relay node. We can see that the outage performance of D 1 is better than D 2 . This is because the distance from the relay to D 2 is longer than the one from the relay to D 1 [39] The OOP is calculated as the probability of the event that both D 1 and D 2 cannot decode their symbols successfully. The simulation and analytical results are in exellent match, which validates the correctness of the closed-form expression of (18), (20), and (21) as in the proposed analysis. From Fig. 5, we can observe that the joint OP of D 1 and D 2 is less than the individual OP of each D 1 and D 2 , which is reasonable. Figure 6 presents the OP D1 versus SNR in dB for different the power allocation scheme. The power allocation coefficient is fixed at a 1 ¼ 0:3 when investigating the outage performance. Moreover, we also conduct the optimal power allocation for the relay as described in Theorem 3. For the adaptive power allocation scheme, we set a 1 ¼ Efjg 2 j 2 g=ðEfjg 1 j 2 g þ Efjg 2 j 2 gÞ and a 2 ¼ 1 À a 1 . Again, we can see that analytical results are in excellent agreement with simulation results. From Fig. 6 we can see that the NOMA system has better outage performance than the OMA system. Different from Fig. 4, the gap of the curves plotted in Fig. 6 is more significant for the cases with fixed or adaptive power allocation and the case of optimal power allocation. The reason is that the probability of the event that both D 1 and D 2 are in outage is less than the probability of the each event that D 1 and D 2 is in outage. Figure 7 plots the OP D2 versus SNR in dB for different power allocation schemes and OP of D 2 in the case of nonlinear EH. From this figure we see that, in the cases of fixed and adaptive power allocation scheme, the outage performance of D 2 in NOMA system is worse than the OMA system. In contrast, the NOMA system with the optimal power allocation has better outage performance than the OMA system. For all schemes, the NOMA system provides better spectral efficiency because two users are served simultaneously. The gap between the curves plotted in Fig. 7 is also insignificant for the two cases with fixed and adaptive power allocation scheme. The reason is that we choose the power allocation coefficient for x 1 and x 2 based on the average channel gains from R to D 1 and D 2 , which leads to the fact that OP of D 1 and OP of D 2 are the same. In addition, we also see that the analytical results are in excellent agreement with the simulation results. On the other hand, in Fig. 7 we demonstrate the OP of D 2 for nonlinear EH model obtained by simulation. Here, we see that the OP for nonlinear model is saturated by a thresholds P th . As we can observe, the different P th results in different OP values. In particular, increasing the saturation threshold, which means more power is harvested at the relay, makes the OP of D 2 decreasing. In addition, for smaller P th , the saturation point appears at smaller SNR.
The effect of imperfect SIC to the system OP is illustrated in Fig. 8. In this case, we use the optimal power allocation scheme and the exact analysis of (22) for reference. The curves in Fig. 8 show that decreasing j, i.e. decreasing residual interference level, results in an improvement in outage performance.In particular, j ¼ 0, i.e., no SIC error, is corresponding to teh best OP of D 1 . This can be explained clearly as follows. Because x 2 is decoded first, then x 1 is decoded in the existence of residual interference by x 2 , the performance of D 1 is definitely affected by the imperfect SIC. In summay, imperfect SIC has no effect on the demodulation of the signal x 2 , but do has an impact on the demodulation of x 1 . In another aspect, Fig. 8 shows the comparison of OP at D 1 in the cases of non-linear and linear EH. It is indicated by Fig. 8 that, for linear EH the OP continuously decreases as the average    Fig. 6 OP of D 1 versus SNR, performance comparison of fixed, adaptive and optimal the power allocation SNR increases, while for non linear EH the OP reaches a saturated state as high values of average SNR. Figure 9 depicts the effect of the power allocation coefficient on the OP. It should be noted that we only set the coefficient for D 1 while the coefficient for D 2 is derived from the condition a 2 ¼ 1 À a 1 . As shown in Fig. 9, different data rate, r 2 , leads to different minimum values of OP. In this figure, we can see that when the transmission rate r 2 is reduced, the power allocation coefficient for D 2 decreases to get better system performance and guarantee the fairness in outage performance of D 1 and D 2 . Figure 10 demonstrates the throughput of NOMA system with and without buffer-aid at relay. As we can see, the system with buffer aid outperforms the ones without buffer aid. The reason is that when the buffer is used, the transmission is carried out with more reliable channel condition. The benefits of the BA relay is remarkable, especially at low SNR values. On the other hand, the throughput of NOMA system is also higher than the throughput of OMA system. It is because in OMA scheme only one packet is transmitted by relay in one time slot, while in NOMA scheme two packets can be transmitted by relay in one time slot. Figure 11 illustrates the throughput of NOMA systems with and without buffer aid at relay in imperfect SIC condition versus SNR in dB. The special case that j ¼ 0 is corresponding to the perfect SIC operation. On the other hand, the throughput of OMA system with and without buffer aid at relay are also plots to compare with NOMA systems. The results show that for both cases with and without buffer aid, the throughput of NOMA system is always higher than the throughput of OMA system. Moreover, the performance gap between the systems with and without buffer aid at high SNR regime is insignificant.
It is due to the trade-off between the probability of successful decoding and the transmission rate. Figure 12 illustrates the effect of the transmission rate on the throughput of each user in the considered SWIPT-NOMA relay system. We can see that the throughput of each user increases to the maximum point, then decreases to zero. This can be explained as follows. At the small data rate, the system is likely to decode the message suscessfully. Therefore, increasing data rate leads to higher throughput. But at some certain point, increasing data rate makes the decoding harder to succeed. As a result, the OP increases quickly and makes the throughput decrease. Another observation is that for both perfect and imperfect SIC cases, the throughput of D 1 is always higher than that of D 2 . It is due to the closer distance from D 1 to R than from D 2 to R, which enhances the probability of successful decoding the signal at D 1 . In addition, the optimal data rates for two users are also different. Specially, for the throughput of D 2 is maximized at the transmission rate of about 1 bps/Hz, while the optimal data rates for D 1 in imperfect and perfect SIC cases are 1.5 bps/Hz and 2 bps/ Hz, respectively. Because the path loss of the R ! D 2 link is larger than that of R ! D 1 link, the optimal transmission rate for D 2 should be smaller than the one for D 1 . In Fig. 13, we plot the average packet delay in our proposed system as a function of the average SNR in dB to verify the result in Sect. 4.5, where r 1 ¼ r 2 ¼ 0:5 bits/slot. As observed from Fig. 13, at high SNR regime, i.e. SNR ! 10 dB, we get the total average packet delay equal to 4 packets. This result confirms that the mathematical derivation in Sect. 4.5 is accurate. On the other hand, the average packet delay at the relay is larger than the average packet delay at the source. It is because when R is selected to transmit data, it must know both g 1 and g 2 , while in receiving mode, it only needs to know h. The delay in this scenario is acceptable when considering the corresponding outage performance and throughput.

Conclusion
In this paper, we proposed a NOMA cooperative relaying network with and without data buffer-aid at relay. In addition, the relay node harvests the energy from the source using the time-switching protocol. We focus on deriving the OP and throughput of the system over Rayleigh fading channels. Moreover, we proposed a power allocation scheme at the relay node, which aims to reduce the feedback cost. Specifically, power allocation only depends on the targeted point on the capacity region of the scheduled users, i.e. c th . Numerical results of the OPs and throughput show that the proposed NOMA downlink relaying system significantly outperforms the equivalent OMA system, especially in the case of imperfect SIC, where the outage performance is decreasing remarkably.
The data buffer-aid employed at the relay helps improve the outage performance of the system. However, the outage performance of the proposed system trades off with the permissible packet delay. The impact of fixed, adaptive and optimal power allocation scheme on the performance of EH-NOMA downlink relaying network was also investigated to confirm that, when optimal power allocation schem is applied, the OP of the system is improved significantly. Furthermore, to take into account the practical characteristics of EH circuits, we provide the simulation results for non-linear EH model. All closed-form With the derivation of closed-form expressions of OP and throughput as well as the optimal power allocation and transmission rate for different variations of the systems, this work provides a comprehensive performance analysis of SWIPT-NOMA systems with BA relay, which is distinguished from other works on the similar topic.
Our proposed relaying network can achieve two goals: (i) the energy efficiency is improved by the harvested energy from ambient RF environment. (ii) the spectrum utilizing efficiency in the proposed system is superior to that of the OMA system. In this model, all nodes are equipped with a single antenna. However, it can be developed for multiple-antenna systems in near-future work.

Appendix A
The goal of this appendix is to derive the OOP of the SWIPT-NOMA system over Rayleigh fading channels. The overall OP of the system can be expressed as From (4), we obtain the closed-form expression of the first term of (47) as To obtain the closed-form expression of the second term of (47), we rewrite the second term of (47) as From (49) and the condition c x 2 !x 1 D 1 [ c x 2 D 2 , we have: By substituting (4), (7), and (8) into (50) and denoting jh 1 j 2 ¼ X, jg 1 j 2 ¼ Y, and jg 2 j 2 ¼ Z, we obtain: As seen from (51), the outage always occurs if a 1 ! 1 1þc th .
Thus, allocating more power to the D 2 is required so that 1 À a 1 ð1 þ c th Þ [ 0 always holds. The condition a 1 \ 1 1þc th is used throughout this paper. For simplicity, we can rewrite (51) as where W min ¼ min c th a 1 /P S ; Based on the properties of conditional probability [40] and the assumption that the channel gains have exponential distributions, we have: where F Y ðyÞ ¼ 1 À e Ày X 2 and f X ðxÞ ¼ 1 X 1 e Àx X 1 are the CDF of X and the PDF of Y, respectively.
Substituting PDF of X into (53) yields By using the Taylor's series expansions of the exponential function, doing some algebras on (54), and using [27, 3.351.4], we obtain the second term of (47) as Fig. 13 The average packets delay versus average SNR; performance comparison for relay decision, priority to transmit and priority to receive ð55Þ We obtain the OP expression of the system by combining (48) and (55). When the transmission power is high, we have c th ( P S . Then, (54) can be approximated as

Appendix B
Due to the imperfect detection at the relay node, it may forward wrongly decoded signals to D 1 and D 2 and cannot apply SIC technique on symbol x 2 at D 1 . Hence, similarly to [41], for any modulation scheme, the dual-hop of the links S ! R ! D 1 or S ! R ! D 2 can be modeled as an equivalent one-hop channel whose output SINR X i , i 2 f1; 2g, at high SNR regime can be tightly approximated. Let X 1 and X 2 denote the SINRs obtained at D 1 and D 2 , respectively [4].
To find the OP of D 1 , from (57), we have the OP expression of D 1 as where Q max ¼ max n 1 a 1 /P S ; n 1 /P S ð1Àa 1 ð1þn 1 ÞÞ n o .
By using the conditional probability [40], we can rewrite (59) as Since the CDF and PDF of X and Y are exponential distribution functions, we have: By using the Taylor series expansions of the exponential function, after some manipulations of (61) using [27, 3.351.4], we obtain the expression of the OP of D 1 as presented in (20). For the case of imperfect SIC, we replace c x 1 D 1 from (9) into (59), then OP are derived similar to (59) to (61). After some manipulations, we have D max ¼ max n 1 /P S ða 1 Àn 1 jð1Àa 2 ÞÞ ; n 1 /P S ð1Àa 1 ð1þn 1 ÞÞ n o ; and this term is used instead of Q max in (61). Next, we calculate the OP expression at D 2 . With the given SINR at D 2 and the notation X 2 ¼ minðc R ; c x 2 D 2 Þ, we can write: Substituting (4) and (6) into (62) yields Then, by applying similar calculations as in Appendix A we can obtain the OP of D 2 as OP D 2 ¼ 1 À Z 1 n 2 P S 1 À F Z n 2 x/P S a 2 À a 1 n 2 ð Þ ! f X x ð Þdx À Á \c th À Á ¼ Pr a 1 P R jg 1 j 2 N 0 \c th ; ð1 À a 1 ÞP R jg 2 j 2 a 1 P R jg 2 j 2 þ 1 \c th ! ¼ Pr XY \ c th a 1 /P S ; XZ\ c th /P S ð1 À a 1 ð1 þ c th ÞÞ : Based on the definition of the conditional probability, we can write: where A ¼ c th a 1 /P S , B ¼ c th /P S ð1Àa 1 ð1þc th ÞÞ . After some manipulations, we get (33), which completes the proof of Theorem 4. h Funding There is no funding available for this work.
Data availability There is no other data and material rather than the one presented in this paper.
Code availability Not applicable.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Ethical approval This article does not contain any studies with human participants or animals performed by any of the authors.