On Outage Probability and Ergodic Rate of Downlink Multi-User Relay Systems with Combination of NOMA, SWIPT, and Beamforming

Abstract

A downlink of multi-user non-orthogonal multiple access (NOMA) relay systems is considered. To improve system performance, relay nodes are used to forward signals from the base station (BS) to the end users, and they are wirelessly powered by energy harvesting from the radio frequency transmitted from the BS. Moreover, beamforming is applied at the BS based on multiple antennas and relay nodes consist of one transmit antenna but several receive ones. The system performance is demonstrated through closed-form expressions of outage probability (OP) and ergodic rate (ER) over Rayleigh fading channels. The proposed system is investigated in two cases of perfect and imperfect successive interference cancellation (SIC), and the imperfect channel state information condition is also taken into consideration. The OP and ER are calculated in many scenarios and the optimal time fraction of energy harvesting corresponding to the minimum OP is discussed. The exact and approximate theoretical results are compared with the simulation result to confirm the proposed theoretical analysis method.

1. Introduction

The Internet of Things (IoT) has received an increasing attention from both industry and academia for recent years [1,2,3]. It is considered to be a crucial means for wireless connections in the era of the fourth industry. The IoT was also included in the advanced wireless standards such as the fourth generation (4G) and has been considering for the fifth generation (5G) of mobile communications [4,5]. Most IoT devices are battery operated or can be located in the deep ocean, outback, and military devices, therefore, battery recharging methods can be economically infeasible or unredeemable. As a result, the IoT system needs to be improved its performance such as lifetime, capacity, spectral and energy efficiencies to meet with requirements of future communication systems.

To support a large multi-user system such as the IoT, the non-orthogonal multiple access (NOMA) is significantly considered due to its high bandwidth efficiency [6,7,8]. These works investigated the performance of the multi-user NOMA system and solved problems through simulation. The coexistence of both OMA and NOMA methods for communication capability on a communication device in wireless networks is proposed in [9], i.e., according to requirement of Quality of Service (QoS), users decide to use OMA or NOMA method for information transmission. This scheme would be risky because of inaccurate recognition of QoS. On the other hand, the problems of mode selection, dynamic user association, and power optimization have been studied for in-band full duplex (IBFD) and NOMA operating networks [10].

There has been a research that takes multi-beam multiple input single output (MISO) NOMA into consideration [11]. However, the author assumed that each cluster has only two users named as the near user and far user. We are going to extend this work with more users in the same beam.

A combination of the NOMA and energy harvesting (EH) was considered in some papers [12,13,14,15,16]. In these papers, the authors proposed the simultaneous wireless information and power transfer (SWIPT) for the NOMA network in which the near user works as the EH relay to assist the far user for improving both network lifetime and spectral efficiency.

Moreover, there have been several papers researching on a combination of the radio frequency (RF) EH and relay communication to enhance the system performance [17,18,19,20,21]. For example, in [17], the dual-hop decode-and-forward (DF) relaying network employing a time switching-based relaying mechanism was analyzed and the proposed system was evaluated via the outage probability. In [18], the authors concentrated on the amplify-and-forward (AF) relaying system and investigated the fixed-gain and variable-gain relaying schemes to deal with interference. In [19], the multiple input multiple output (MIMO) full-duplex system was considered and the optimal throughput was analyzed. Since antennas are configured at relay nodes, the energy consumption of relay nodes depends on the number of antennas.

The small cell station can harvest energy transmitted by a macro-cell base station or by solar panels [22]. This work is close to our proposed system; however, it only investigated EH for the small cell station of a two-tier heterogeneous small cell network, and the authors did not take the NOMA aspect into consideration. On the other hand, to harvest energy from RF signal or solar panels, the base station uses different architectures, it is hard to combine in one device. The problem of EH efficiency from RF is shown in [23], in this work, the authors proposed a low complexity energy transfer and transmission adaptation algorithm in order to determine the optimal power station’s transmission power. However, this work considered an orthogonal frequency multiple access scheme.

A combining SWIPT-NOMA downlink relaying protocol is a potential approach to achieve high energy efficiency and performance, i.e., prolonging the lifetime and improving the throughput of wireless networks, especially for designing the future 5G and 6G networks. To the best of our knowledge, there is a lack of studies on a combination of the SWIPT and NOMA in relaying network systems.

Due to a big gap between the active sensitivity of energy harvester and that of decoder (i.e., −10 dBm for EH and −60 dBm for information receiving), the SWIPT-based network is only appropriate for short-distance transmission. To overcome this limitation, the authors of [24,25,26] proposed the systems in which individual EH and information receivers are used depending on the time switching (TS) and power slitting (PS) to supply power for wireless devices. In particular, [27] proposed three wireless power transfer policies. Namely, the cooperative power beacons (PBs) power transfer, the best PBs power transfer, and the nearest PBs power transfer. Since the SWIPT with PBs is fully controllable, it can provide various applications with high quality of service (QoS) requirements. The RF EH from PBs was proposed and investigated in [28,29]. However, [28,29], and also other previous works [25,26] only considered point-to-point communication systems and frequencies used to transmit PBs and information during EH period are different. In these papers, the PBs are transmitted continuously; however, this approach is not an optimal solution in the sense of saving the PBs power and requires a high practical cost as well.

It is clear that an application of the EH from the surrounding RF to supply power for wireless devices is an avocational research for the next generation networks, 5G and 6G. In our initial works [30,31], we have investigated the NOMA beamforming relay system and then concentrated on a combination of the SWIPT, NOMA and beamforming technique for the downlink of the multi-user system at the beginning step.

The contributions of this paper are summarized as follows

• We propose a downlink MIMO-NOMA with SWIPT relay system which aims to improve the energy and spectrum efficiencies of the system. To reduce the complexity and mitigate inter-user interference, users in the system are divided into several clusters and a suitable user clustering mode is found out. The relays are employed to mitigate the effect of fading and reduce the number of training sequences for estimating the CSI at the base station.
• We derive the closed-form expressions of the outage probability (OP) and ergodic rate (ER) for each user in any cluster. For the practical purpose, we investigate the proposed MIMO-NOMA with SWIPT relay system over Rayleigh channel. The time duration of EH is optimized in the sense of minimum OP. The analysis result provides insights to understand the performance of the MIMO-NOMA system with SWIPT through the mathematical expressions.
• We consider the system in the case of imperfect channel state information (CSI) caused by the downlink channel estimation error and both cases of perfect and imperfect successive interference cancellation (SIC). The result shows that the system performance is reduced significantly due to the imperfect CSI and imperfect SIC. All analysis results are compared with simulation results to confirm the correctness of the derived mathematical expressions.

The rest of the paper is organized as follows. Section 2 represents the proposed system model, illustrating CSI requirement, channel model, and energy harvesting technique. Section 3 concentrates on performance analysis, deriving the closed-form expressions of the OP and ER. Section 4 demonstrates numerical results and discussion. Finally, Section 5 concludes the paper.

For the convenience, in the following

Table 1. The notations.

Notation	Description
$\mathrm{Pr}$	Probability
$F_X\left(x\right)$	Cumulative distribution function (CDF)
$f_X\left(x\right)$	Probability density function (PDF)
$\mathcal{CN}(\mu ,{\sigma }^2)$	A circularly symmetric complex Gaussian RV x with mean $\mu$ and variance ${\sigma }^2$
$\mathbb{E}\left\{·\right\}$	The statistical expectation operator
$\mathsf{\Gamma }(·)$	Gamma function [32]
${\mathcal{K}}_n\left(·\right)$	The second kind of Bessel function, order n [32]
${\mathrm{E}}_n\left(z\right)$	Exponential integral function n [32]
$G_{pq}^{mn}\left(x{|}_{b_s}^{a_r}\right)$	Meijer’s G-Function [9.3] [32]
M	Number of relay nodes
N	Number of users in each cluster
$N_t,N_r$	Transmission antennas of BS and reception antennas of relay
${\rho }_{\mathcal{A}}$	Correlation coefficient
$\alpha$	Time switching ratio
$\eta$	Conversion efficiency
r	Data rate threshold

, we provide notations along with their descriptions used in this paper.

2. System Model

In this paper, a downlink of the multi-user (MUs) relaying communication network is investigated. In this system, the BS communicates with users via a support of relay (R) nodes as illustrated in Figure 1. There are M relay nodes that are deployed in the cell and have constrained energy. To improve the energy efficiency, the relay node harvests energy from surrounding RF. The BS node and each relay node have $N_t$ and $N_r$ antennas, respectively, satisfying the following condition $(N_t\ge (M-1))$. While users have a single antenna because of the limited size, the BS can be considered to be having massive antennas. The block diagram of the first hop from the BS node to the relay node is shown in Figure 2, in which the signal, which is transmitted to all users in the mth cluster, is denoted as ${\mathbf{x}}_{\mathrm{S},\mathit{m}}$. The signal ${\mathbf{x}}_{\mathrm{S},\mathit{m}}$ is multiplied by $N_t\times 1$ beamforming vector ${\mathbf{w}}_m$ before being transmitted over channel, and the beamforming vector ${\mathbf{w}}_m$ is assumed to be normalized $\parallel {\mathbf{w}}_m{\parallel }^2=1$.

Users that are close to a relay node, are grouped into a cluster with the help of the GPS algorithm (The spatial direction of users can be found via various methods/technologies such as GPS or user location tracking algorithms.). We assume that there are N users in every cluster and the NOMA scheme is applied to improve the spectrum efficiency, i.e., signals are superposed in the power domain at the BS. Moreover, the NOMA techniques and DF protocol are performed at the relay node and users to remove the inter-user interference. In the next section, the proposed multi-relay NOMA scheme is described in detail.

The operation of the MIMO-NOMA system with SWIPT and beamforming in DF protocol can be summarized as follows. A complete signal transmission procedure from the BS to the ${\mathrm{D}}_n$ under assistance of SWIPT relays in the proposed MIMO-NOMA system takes three time slots. In the first time slot, the relay harvests energy from broadcasts of BS over downlink channel. The duration time of this slot is denoted by $\alpha T$. In the second time slot, the BS sends the signal to relays after performing superposition coding based on the characteristic of each user. In the third time slot, $R_m$ re-encodes the received signals and then assigns the power for these signals before forwarding them to the users.

2.1. CSI Condition and Channel Model

For transmitting beamforming, the BS requires the CSI. In the proposed model, the BS receives the CSI by up-link channel estimation. Herein, at the beginning of a working section, relay nodes send pilot symbols to the BS. Since the BS and relay node are imperfectly synchronized, all pilot signals transmitted from the relay node to the BS are non-orthogonal and imperfect. As a result, the perfect CSI obtained at the BS is absolutely challenging owing to both causes as channel estimation error and feedback delay. Assuming that the minimum mean squared error (MMSE) estimation is used, the difference between the estimated channel matrix ${\widehat{\mathbf{h}}}_m$ and the actual channel matrix ${\mathbf{h}}_m$ can be illustrated as the following equation [33]

$$\begin{array}{c}\hfill {\mathbf{h}}_m={\widehat{\mathbf{h}}}_m+\mathbf{e},\end{array}$$

(1)

here $\mathbf{e}$ denotes the channel estimation error vector with independent and identically distributed (i.i.d.) zero mean and unit variance complex Gaussian distributed entries.

To allocate power for users, the relay node also require CSI. The difference between the estimated channel coefficient ${\widehat{g}}_n$ and actual channel coefficient $g_n$ can be also demonstrated as

$$\begin{array}{c}\hfill g_n={\widehat{g}}_n+e,\end{array}$$

(2)

here e denotes an error of the channel coefficient estimation between the $R_m$ and the ${\mathrm{D}}_{\mathit{n}}$.

In addition, ${\widehat{\mathsf{\Omega }}}_{\mathcal{A}}={\mathsf{\Omega }}_{\mathcal{A}}-{\sigma }_{{\epsilon }_{\mathcal{A}}}^2$ with $\mathcal{A}\in \{{\mathrm{SR}}_{\mathit{m}},{\mathrm{RD}}_{\mathit{n}}\}$ is the channel gain norm of the ${\widehat{\mathbf{h}}}_m$. Let $0\le {\rho }_{\mathcal{A}}\le 1$ is the correlation between coefficients of channel estimation error and coefficients of the perfect channel. Let ${\sigma }_{{\epsilon }_{\mathcal{A}}}^2={\rho }_{\mathcal{A}}{\mathsf{\Omega }}_{\mathcal{A}}$ and ${\widehat{\mathsf{\Omega }}}_{\mathcal{A}}=(1-{\rho }_{\mathcal{A}}){\mathsf{\Omega }}_{\mathcal{A}}$ denote the variance norm and the variable of channel estimation error, respectively.

In this paper, the zero-force beamforming (ZFBF) is employed at the BS to deal with a trade-off between the complexity of implementation and performance of the system. The weight ${\mathbf{w}}_m$ for the mth cluster is designed to cancel the interference from the other clusters.

Every antenna of the BS sends the superposition code, which consists of N signals of the mth cluster and can be described by $x_{\mathrm{m}}$.

$$\begin{array}{c}\hfill x_{\mathrm{m}}=\sum _{n=1}^N\sqrt{a_n{P}_{\mathrm{S}}}{x}_{m,n},\end{array}$$

(3)

here $a_n$ illustrates the power allocation coefficient of the $(m,n)$th user with ${\sum }_{n=1}^N{a}_n=1$. At the output of antenna, this signal is built with the beam vector, ${\mathbf{w}}_m$, that depends on the ZFBF. As a result, the signal corresponding to the mth cluster is given by the following equation.

$$\begin{array}{c}\hfill {\mathbf{x}}_{\mathrm{S},\mathrm{m}}={\mathbf{w}}_m{x}_{\mathrm{S},\mathrm{m}}={\mathbf{w}}_m\sum _{n=1}^N\sqrt{a_n{P}_{\mathrm{S}}}{x}_{m,n}.\end{array}$$

(4)

The $N_t\times 1$ signal vector ${\mathbf{x}}_{\mathrm{S}}={\mathbf{x}}_{\mathrm{S},1}+\cdots +{\mathbf{x}}_{\mathrm{S},\mathrm{M}}$ is transmitted to all relay nodes over a single cell. Hence, the received signal vector at the relay node can be calculated as

$$\begin{array}{cc}\hfill {\mathbf{y}}_R& ={\mathbf{h}}_m{\mathbf{x}}_S+{\mathbf{n}}_m\hfill \\ \hfill & =\underset{m\mathrm{th} \mathrm{beam}}{\underbrace{{\mathbf{h}}_m{\mathbf{w}}_m\sum _{n=1}^N\sqrt{a_n{P}_{\mathrm{S}}}{x}_{\mathrm{m},\mathrm{n}}}}+\underset{\mathrm{other} \mathrm{beam}}{\underbrace{{\mathbf{h}}_m\sum _{k\ne m}{\mathbf{x}}_{s,k}}}+{\mathbf{n}}_m,\hfill \end{array}$$

(5)

where ${\mathbf{n}}_m=[n_1,\cdots ,n_N]\in {\mathbb{C}}^{1\times N}$ with $n_m\sim \mathcal{CN}(0,{\sigma }_m^2)$ denotes the additive white Gaussian noise (AWGN) at the mth relay (The channel gain formulation, the interference model, the signaling overhead can be found in [34]).

According to the NOMA technique, the power domain is used for multiple access of the MUs at the same time and frequency, i.e., different power levels are allocated to different users. Consequently, ${\mathbf{w}}_m$ can be considered to be a projection of ${\mathbf{h}}_m$ in a null space of interfered channels according to the mth cluster. To maximize the channel gain and mitigate the inter-cluster interference, the ${\mathbf{w}}_m$ can be defined as

$$\begin{array}{c}\hfill {\mathbf{w}}_m=\frac{{\mathsf{\Pi }}_m{\mathbf{h}}_m}{\parallel {\mathsf{\Pi }}_m{\mathbf{h}}_m\parallel },\end{array}$$

(6)

here ${\mathsf{\Pi }}_m={\mathbf{I}}_N-{\mathbf{H}}_m{\left({\mathbf{H}}_m^H{\mathbf{H}}_m\right)}^{-1}{\mathbf{H}}_m^H$, ${\mathbf{H}}_m$ denotes an extended channel matrix that does not consist of ${\mathbf{h}}_m$ and its construction is represented as

$$\begin{array}{c}\hfill {\mathbf{H}}_m={[{\mathbf{h}}_1,{\mathbf{h}}_2,\cdots ,{\mathbf{h}}_{m-1},{\mathbf{h}}_{m+1},\cdots ,{\mathbf{h}}_M]}^T.\end{array}$$

(7)

The condition to transmit the signal from the BS to the mth relay node is that there is a dimension greater than zero in the null space of ${\mathbf{H}}_m$. As a result, we have

$$\begin{array}{c}\hfill {\mathbf{h}}_m^H{\mathbf{w}}_j=0,\forall m\ne j.\end{array}$$

(8)

In other words, the condition in (8) is used to cancel completely interference during of the first hop of inter-relay. In addition, owning to the assistance the multiple relays, the size of matrix for designing beamforming is not so large and the computational complexity is bearable.

In this paper, a beam is assumed to be designed perfectly. Therefore, the output signal of the SIC structure at the mth relay node can be illustrated as

$$\begin{array}{cc}\hfill {\mathbf{y}}_{R,SIC}& ={\mathbf{h}}_m{\mathbf{w}}_m\sum _{n=1}^N\sqrt{a_n{P}_{\mathrm{S}}}{x}_{m,n}+{\mathbf{n}}_m\hfill \\ \hfill & =\underset{\mathrm{desired} \mathrm{signal} \mathrm{of} (m,n)\mathrm{th} \mathrm{user}}{\underbrace{{\mathbf{h}}_m{\mathbf{w}}_m\sqrt{a_n{P}_{\mathrm{S}}}{x}_{m,n}}}+\underset{\mathrm{interference} \mathrm{of} \mathrm{other} \mathrm{users}}{\underbrace{{\mathbf{h}}_m{\mathbf{w}}_m\sum _{i=n+1}^N\sqrt{a_i{P}_{\mathrm{S}}}{x}_{m,i}}}+\underset{\mathrm{interference} \mathrm{of} \mathrm{imperfect} \mathrm{SIC}}{\underbrace{{\mathbf{h}}_m{\mathbf{w}}_m\sum _{k=1}^{n-1}\sqrt{{\xi }_1{a}_k{P}_{\mathrm{S}}}{x}_{m,k}}}+{\mathbf{n}}_m.\hfill \end{array}$$

(9)

The term ${\mathbf{h}}_m{\mathbf{w}}_m{\displaystyle \sum _{k=1}^{n-1}}\sqrt{{\xi }_1{a}_k{P}_{\mathrm{S}}}{x}_{m,k}$ is equal to zero if the SIC is perfect.

2.2. Energy Harvesting

In this work, the time switching (TS) is investigated owning to its higher performance and easier implementation than the power splitting when these two protocols are applied for the multi-user MIMO networks with the EH technique as demonstrated in [35]. The TS architecture [36] as in Figure 3 is supposed for the operation of the relay node (The proposed analysis approach can also be applied to the power spitting EH model [37]). In this figure, T is a block time to transmit a certain block of information from the source node to the destination node. The communication process is split into three time slots, the first slot, $\alpha T$ ($0\le \alpha \le 1$), for energy harvesting at the relay node, the second one, $(1-\alpha )T/2$, for the S-to-R information transmission and the third one, $(1-\alpha )T/2$, for the R-to-D information transmission. In case the relay node does not have an energy buffer to store the harvested energy, the harvest-use (HU) architecture [38] can be used.

Assuming that the EH process is only performed from the received signals in $\alpha T$ time interval of each period. As a result, the harvested energy can be calculated as ([36], Equation (2)).

$$\begin{array}{c}\hfill E_h=\alpha T\eta P_{\mathrm{S}}{\left|{\mathbf{h}}_m{\mathbf{w}}_m\right|}^2,\end{array}$$

(10)

where $\eta$, $0<\eta <1$, illustrates the energy conversion efficiency and depends on the harvester quality and $P_{\mathrm{S}}$ denotes the transmission power of the source.

According to the HU architecture, the relay node uses total harvested energy in the harvesting phase to forward the received data to users. Hence, the transmission power of the relay node is given by

$$\begin{array}{cc}\hfill P_{\mathrm{R}}& =\frac{2\alpha \eta P_{\mathrm{S}}}{(1-\alpha )}|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2=\varphi P_{\mathrm{S}}{\left|{\mathbf{h}}_m{\mathbf{w}}_m\right|}^2,\hfill \end{array}$$

(11)

here $\varphi =2\alpha \eta /(1-\alpha )$.

The large-scale fading coefficient of the $(m,n)$th user is denoted by $\sqrt{d_n^{-\beta }}$, where $d_n$ is the distance between the relay node and the $(m,n)$th user and $\beta$ is the path loss factor. Moreover, ${\tilde{g}}_n$ represents the small-scale fading coefficient of the $(m,n)$th user, leading to the fact that $g_n={\tilde{g}}_n\sqrt{d_n^{-\beta }}$ and $g_n\sim \mathcal{C}\mathcal{N}(\mu ,{\mathsf{\Omega }}_{{\mathrm{RD}}_n})$, where ${\mathsf{\Omega }}_{{\mathrm{RD}}_n}=\mathbb{E}\left\{\right|g_n{|}^2\}$. The AWGN noise vector at the $(m,n)$th user is $w_{D_n}\sim \mathcal{C}\mathcal{N}(0,N_0)$ with $N_0$ is the noise variance at the receiver. All channel gains have the i.i.d Rayleigh distribution. Without any loss of generality, we assume that $d_1>d_2,\dots ,>d_N$. Hence, the channel gains are sorted according to the ascending order $|g_1{|}^2<,\dots ,<{\left|g_N\right|}^2$.

During the second time slot, the relay node re-encodes and forwards the messages to the user. Since the transmission power of the relay node, which depends on the quality of the energy harvester, is always smaller than that of the source. Therefore, it is necessary to re-allocate power to be fair performance of users. The output signal of the SIC architecture of the $(m,n)$th user is given as

$$\begin{array}{cc}\hfill y_{D_n}=& \underset{\mathrm{desired} \mathrm{signal} \mathrm{of} (m,n)\mathrm{th} \mathrm{user}}{\underbrace{g_n\sqrt{b_n{P}_{\mathrm{R}}}{x}_{m,n}}}+\underset{\mathrm{interference} \mathrm{of} \mathrm{other} \mathrm{users}}{\underbrace{g_n{\sum }_{i=n+1}^N\sqrt{b_i{P}_{\mathrm{R}}}{x}_{m,i}}}+\underset{\mathrm{interference} \mathrm{of} \mathrm{imperfect} \mathrm{SIC}}{\underbrace{{\sum }_{k=1}^{n-1}\sqrt{{\xi }_2{b}_k{P}_{\mathrm{R}}}{g}_n{x}_{m,k}}}+w_{D_n}.\hfill \end{array}$$

(12)

At the relay node, the SIC removes the interference signal of the other users having stronger power, therefore, in the case of perfect SIC, the signal to interference plus noise ratio (SINR) of the $(m,n)$th user at the relay node is expressed by ${\gamma }_{R,n}$, and from (9), it can be calculated as

$$\begin{array}{c}\hfill {\gamma }_{R,n}=\frac{P_{\mathrm{S}}{a}_n|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2}{{\sum }_{i=n+1}^N{P}_{\mathrm{S}}{a}_i|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2+{\sigma }_{R,n}^2}.\end{array}$$

(13)

In the case of imperfect SIC, due to channel estimation error and remaining interference, the term ${\mathbf{h}}_m{\mathbf{w}}_m{\displaystyle \sum _{k=1}^{n-1}}\sqrt{{\xi }_1{a}_k{P}_{\mathrm{S}}}{x}_{m,k}$ in (9) is not equals to zero and it depends on the quality of the SIC structure. As a result, the instantaneous SINR in this case is given by

$$\begin{array}{c}\hfill {\gamma }_{R,n}^{ipSIC}=\frac{P_{\mathrm{S}}{a}_n|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2}{{\displaystyle \sum _{i=n+1}^N}{P}_{\mathrm{S}}{a}_i|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2+{\left|{\mathbf{h}}_m{\mathbf{w}}_m\right|}^2{\xi }_1{\displaystyle \sum _{k=1}^{n-1}}{a}_k{P}_{\mathrm{S}}+{\sigma }_{\mathrm{R}}^2},\end{array}$$

(14)

where ${\xi }_1$ represents the impact level of residual interference at the relay.

From (12), the SINR at the $(m,n)$th user in the case of perfect SIC is calculated as follows. Note that, in case the power for $(m,i)$th user is larger than that for $(m,n)$th user, i.e., $n\ge i$ or $b_n\le b_i$, then the SINR of $x_{m,i}$ at the $(m,n)$th user can also be expressed as

$${\gamma }_{D,i}=\{\begin{array}{}\mathrm{(15a)}& \frac{P_{\mathrm{R}}{b}_i{\left|g_n\right|}^2}{{\sum }_{k=i+1}^N{b}_k{P}_{\mathrm{R}}|g_n{|}^2+{\sigma }_{D,n}^2},\hfill & \mathrm{if} i < N,\hfill \mathrm{(15b)}& \frac{P_{\mathrm{R}}{b}_i{\left|g_n\right|}^2}{{\sigma }_{D,n}^2},\hfill & \mathrm{if} i = N.\hfill \end{array}$$

When $x_{m,i}$ is decoded successfully, i.e., ${\gamma }_{D,i}\ge {\gamma }_{{\mathrm{th}}_{\mathit{i}}}$, the $x_{m,i}$ is removed from the $(m,n)$th user’s SIC structure, here ${\gamma }_{{\mathrm{th}}_{\mathit{i}}}$ is the threshold of SINR of the user. The SIC performs continuously at the $(m,n)$th user until its own signal is decoded successfully.

In the case of imperfect SIC, the SINR of $x_{m,i}$ at the $(m,n)$th user with $n\ge i$ is given as

$$\begin{array}{c}\hfill {\gamma }_{D,i}^{ipSIC}=\frac{P_{\mathrm{R}}{b}_i{\left|g_n\right|}^2}{{\sum }_{k=i+1}^N{b}_k{P}_{\mathrm{R}}|g_n{|}^2+{\sum }_{j=1}^{i-1}{\xi }_2{b}_j{P}_{\mathrm{R}}{\left|g_n\right|}^2+{\sigma }_{D,n}^2}.\end{array}$$

(16)

where ${\xi }_2$ represents the impact level of residual interference at the $(m,n)$th user.

Moreover, with the DF protocol, the end-to-end SINR is defined as the minimum value of the two hops, $\mathrm{BS}\to \mathrm{R}$ and $\mathrm{R}\to {\mathrm{D}}_{\mathit{n}}$, consequently, we have

$$\begin{array}{c}\hfill {\gamma }_{e2e}=min({\gamma }_{\mathrm{R},\mathit{n}},{\gamma }_{\mathrm{D},\mathit{n}}).\end{array}$$

(17)

3. Performance Analysis

3.1. Outage Probability

The outage probability (OP) is defined as the probability of an event in which the transmission rate of the system is lower than the minimum required data rate. In this section, the OP of the considered $(m,n)$th user is derived in two cases of the perfect and imperfect SICs.

Let $r_1$ and $r_n$ (bit/s/Hz) denote the minimum required data rate of two links, $\mathrm{BS}\to \mathrm{R}$ and $\mathrm{R}\to {\mathrm{D}}_{\mathit{n}}$, respectively. To simplify calculations, set $r_1=r_n=r$. Therefore, the OP of the $(m,n)$th user can be defined as

$$\begin{array}{cc}\hfill {\mathrm{OP}}_{{\mathrm{D}}_{\mathit{n}}}=& Pr\left[\frac{1-\alpha }{2}{log}_2\left(1+{\gamma }_{e2e}\right)\right]\hfill \\ \hfill =& Pr\left({\gamma }_{e2e}<2^{\frac{2r}{1-\alpha }}-1\right)\hfill \end{array}$$

(18)

By substituting ${\gamma }_{e2e}$ from (17) into (18), the OP of $(m,n)$th user can be calculated as

$$\begin{array}{cc}\hfill {\mathrm{OP}}_{{\mathrm{D}}_{\mathit{n}}}& =Pr\left[min({\gamma }_{{\mathrm{R}}_{\mathit{n}}},{\gamma }_{{\mathrm{D}}_{\mathit{n}}})<{\gamma }_{\mathrm{th}}\right]\hfill \\ \hfill & =1-Pr\left({\gamma }_{{\mathrm{R}}_{\mathit{n}}}\ge {\gamma }_{\mathrm{th}},{\gamma }_{{\mathrm{D}}_{\mathit{n}}}\ge {\gamma }_{\mathrm{th}}\right),\hfill \end{array}$$

(19)

here ${\gamma }_{\mathrm{th}}=2^{\frac{2r}{1-\alpha }}-1$.

Based on the equations of SNR given in (13) and (15a), the OP in case of perfect SIC can be calculated as in (20).

$$\begin{array}{c}\hfill {\mathrm{OP}}_{{\mathrm{D}}_{\mathit{n}}}=1-Pr\left(\frac{P_{\mathrm{S}}{a}_n|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2}{{\sum }_{i=n+1}^N{P}_{\mathrm{S}}{a}_i|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2+{\sigma }_{\mathrm{R}}^2}\ge {\gamma }_{\mathrm{th}},\frac{P_{\mathrm{R}}{b}_n{\left|g_n\right|}^2}{{\sum }_{i=n+1}^N{b}_i{P}_{\mathrm{R}}{\left|g_n\right|}^2+{\sigma }_{\mathrm{D},n}^2}\ge {\gamma }_{\mathrm{th}}\right),\end{array}$$

(20)

In the case of imperfect SIC, from (14) and (16), the outage probability is derived as follows.

$$\begin{array}{c}\hfill {\mathrm{OP}}_{D_{\mathit{n}}}^{ipSIC}=1-Pr\left(\begin{array}{c}\frac{P_{\mathrm{S}}{a}_n|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2}{{\sum }_{i=n+1}^N{P}_{\mathrm{S}}{a}_i|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2+{\left|{\mathbf{h}}_m{\mathbf{w}}_m\right|}^2{\xi }_1{\sum }_{k=1}^{n-1}{a}_k{P}_{\mathrm{S}}+{\sigma }_{\mathrm{R}}^2}\ge {\gamma }_{\mathrm{th}},\\ \frac{P_{\mathrm{R}}{b}_n{\left|g_n\right|}^2}{{\sum }_{i=n+1}^N{b}_i{P}_{\mathrm{R}}|g_n{|}^2+{\sum }_{k=1}^{n-1}{\xi }_2{b}_k{P}_{\mathrm{R}}{\left|g_n\right|}^2+{\sigma }_{\mathrm{D},n}^2}\ge {\gamma }_{\mathrm{th}}\end{array}\right).\end{array}$$

(21)

Set $X=|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2$ and $Y=|g_n{|}^2$ for notation convenience. Since the normalized ${\mathbf{w}}_m$ is designed independently with ${\mathbf{h}}_m$, the $|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2$ is Chi-square distributed with $2K$ degree of freedom (If a variable is random and has Chi-square distribution with K degree of freedom, it becomes sum of Rayleigh distribution ([39], p. 16)), and $K=N_t(N_r-(M+1))$[33], and $|g_n{|}^2$ is Chi-square distributed with two degree of freedom. Hence, the CDF and PDF of X are given as

$$\begin{array}{c}\hfill F_X\left(x\right)=1-exp\left(-\frac{x}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right)\sum _{k=0}^{K-1}\frac{1}{k!}{\left(\frac{x}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right)}^k,\end{array}$$

(22)

$$\begin{array}{c}\hfill f_X\left(x\right)=\frac{x^{K-1}}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}exp\left(-\frac{x}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right).\end{array}$$

(23)

Based on the ordered statistics, the PDF of channel gain $|g_n{|}^2$ can be written as ([40], Equation (6.58)).

$$\begin{array}{cc}\hfill f_{|g_n{|}^2}\left(y\right)=& \frac{N!}{(N-n)!(n-1)!}\sum _{j=0}^{N-n}{(-1)}^j\left(\genfrac{}{}{0pt}{}{N-n}{j}\right)f_{|g_j{|}^2}\left(y\right){\left[F_{|g_j{|}^2}\left(y\right)\right]}^{n+j-1},\hfill \end{array}$$

(24)

where

$$\begin{array}{c}\hfill f_{|g_j{|}^2}\left(y\right)=\frac{1}{(1-\rho ){\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}exp\left(-\frac{y}{(1-\rho ){\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}\right),\end{array}$$

(25)

$$\begin{array}{c}\hfill F_{|g_j{|}^2}\left(y\right)=1-exp\left(-\frac{y}{(1-\rho ){\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}\right).\end{array}$$

(26)

The corresponding CDF of $|g_n{|}^2$ is given by

$$\begin{array}{cc}\hfill F_{|g_n{|}^2}\left(y\right)=& \frac{N!}{(N-n)!(n-1)!}\sum _{j=0}^{N-n}\frac{{(-1)}^j}{n+j}\left(\genfrac{}{}{0pt}{}{N-n}{j}\right){\left[1-exp\left(-\frac{y}{(1-\rho ){\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}\right)\right]}^{n+j},\hfill \end{array}$$

(27)

Using the binomial expansion for (27) we have

$$\begin{array}{cc}\hfill F_{|g_n{|}^2}\left(y\right)=& \frac{N!}{(N-n)!(n-1)!}\sum _{j=0}^{N-n}\frac{{(-1)}^j}{n+j}\left(\genfrac{}{}{0pt}{}{N-n}{j}\right)\sum _{\ell =0}^{n+j}{(-1)}^{\ell }\left(\genfrac{}{}{0pt}{}{n+j}{\ell }\right)exp\left(-\frac{\ell y}{(1-\rho ){\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}\right).\hfill \end{array}$$

(28)

As a result, the OP in two cases of the perfect and imperfect SICs is evaluated through the CDF and PDF of both S-R and R-D hops.

The Equation (20) can be represented as

$$\begin{array}{c}\hfill {\mathrm{OP}}_{{\mathrm{D}}_{\mathit{n}}}\stackrel{\mathsf{\Delta }}{=}1-Pr\left(X\ge \frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{R}}^2}{P_{\mathrm{S}}(a_n-{\gamma }_{\mathrm{th}}\tilde{a})},XY\ge \frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{D},n}^2}{P_{\mathrm{S}}\varphi (b_n-{\gamma }_{\mathrm{th}}\tilde{b})}\right),\end{array}$$

(29)

where $\varphi =\frac{2\alpha \eta }{1-\alpha }$, $\tilde{a}={\sum }_{i=n+1}^N{a}_i$ and $\tilde{b}={\sum }_{i=n+1}^N{b}_i$. The $\stackrel{\mathsf{\Delta }}{=}$ means that conditions, the $a_n>{\gamma }_{\mathrm{th}}\tilde{a}$ and $b_n>{\gamma }_{\mathrm{th}}\tilde{b}$, need to be satisfied. If not, $a_n\le {\gamma }_{\mathrm{th}}\tilde{a}$ or $b_n\le {\gamma }_{\mathrm{th}}\tilde{b}$, the outage always happens because $X,Y\in (0,\infty )$ is considered only. Hence, the power allocated for the ${\mathrm{D}}_n$ is more than the power of the others.

Based on the theory of joint probability of two random variables [41], the Equation (29) can be represented as

$$\begin{array}{cc}\hfill {\mathrm{OP}}_{{\mathrm{D}}_{\mathit{n}}}& =1-{\int }_u^{\infty }\left[1-F_Y\left(\frac{v}{x}\right)\right]f_X\left(x\right)dx\hfill \\ \hfill & =1-\underset{{\mathcal{I}}_1}{\underbrace{{\int }_u^{\infty }{f}_X\left(x\right)dx}}+\underset{{\mathcal{I}}_2}{\underbrace{{\int }_u^{\infty }{F}_Y\left(\frac{v}{x}\right)f_X\left(x\right)dx}},\hfill \end{array}$$

(30)

where $u=\frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{R}}^2}{P_{\mathrm{S}}(a_n-{\gamma }_{\mathrm{th}}\tilde{a})}$ and $v=\frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{D},n}^2}{P_{\mathrm{S}}\varphi (b_n-{\gamma }_{\mathrm{th}}\tilde{b})}$.

Substituting (23) into (30), we obtain ${\mathcal{I}}_1$ as

$$\begin{array}{c}\hfill {\mathcal{I}}_1=\frac{1}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}{\int }_u^{\infty }{x}^{K-1}exp\left(-\frac{x}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right)dx.\end{array}$$

(31)

Thank to the help of ([32], Equation (3.351.2)), we have

$$\begin{array}{c}\hfill {\mathcal{I}}_1=\frac{1}{\mathsf{\Gamma }\left(K\right)}\mathsf{\Gamma }\left(K,\frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{R}}^2}{P_{\mathrm{S}}(a_n-{\gamma }_{\mathrm{th}}\tilde{a})(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right).\end{array}$$

(32)

Substituting (23) and (28) into (30), we obtain ${\mathcal{I}}_2$ as in (33).

$$\begin{array}{cc}\hfill {\mathcal{I}}_2=& \frac{N!}{(N-n)!(n-1)!}\sum _{j=0}^{N-n}\frac{{(-1)}^j}{n+j}\left(\frac{N-n}{j}\right)\sum _{\ell =0}^{n+j}{(-1)}^{\ell }\left(\genfrac{}{}{0pt}{}{n+j}{\ell }\right)\frac{1}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}\hfill \\ \hfill & \times {\int }_u^{\infty }{x}^{K-1}exp\left(-\frac{\ell v}{x(1-\rho ){\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}\right)exp\left(-\frac{x}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right)dx.\hfill \end{array}$$

(33)

To calculate exactly ${\mathcal{I}}_2$, we can use expansion in series for exponential function as in [42], i.e., $exp\left(-\frac{a}{x}\right)={\sum }_{t=0}^{\infty }\frac{{(-1)}^t}{t!}{\left(\frac{a}{x}\right)}^t$. Thus, ${\mathcal{I}}_2$ is rewritten by

$$\begin{array}{cc}\hfill {\mathcal{I}}_2=& \frac{N!}{(N-n)!(n-1)!}\sum _{j=0}^{N-n}\frac{{(-1)}^j}{n+j}\left(\genfrac{}{}{0pt}{}{N-n}{j}\right)\sum _{\ell =0}^{n+j}{(-1)}^{\ell }\left(\genfrac{}{}{0pt}{}{n+j}{\ell }\right)\frac{1}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}\hfill \\ \hfill & \times \sum _{t=0}^{\infty }\frac{{(-1)}^t}{t!}{\left(\frac{\ell v}{(1-\rho ){\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}\right)}^t\underset{{\mathsf{\Delta }}_1}{\underbrace{{\int }_u^{\infty }{x}^{q}exp\left(-\frac{x}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right)dx}},\mathrm{where}q=K-1-t.\hfill \end{array}$$

(34)

To calculate ${\mathsf{\Delta }}_1$ in (34), we consider two cases of power index, i.e., $q\ge 0$ and $q<0$. Using ([32], Equation (3.351.2)) for the case of $q\ge 0$ and exponential integral function for the case of $q<0$. The exponential integral function have been defined by Schloemilch as in ([43], [5.1.4]).

$$\begin{array}{c}\hfill E_n\left(z\right)={\int }_1^{\infty }\frac{e^{-zt}}{t^n}dt.\end{array}$$

(35)

Finally, we have ${\mathsf{\Delta }}_1$ as in (36) and (37).

$$\begin{array}{cc}\hfill & {\mathsf{\Delta }}_1={\left(\frac{1}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right)}^{-q-1}\mathsf{\Gamma }\left(q+1,\frac{u}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right),\mathrm{if} q\ge 0.\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & {\mathsf{\Delta }}_1={\left(\frac{1}{u}\right)}^{q-1}{E}_q\left(\frac{u}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right),\mathrm{if} q<0\hfill \end{array}$$

(37)

In the case of the transmit power is high enough, i.e., $u=\frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{R}}^2}{P_{\mathrm{S}}(a_n-{\gamma }_{\mathrm{th}}\tilde{a})}\to 0$, we have ${\mathcal{I}}_1=1$. Hence, the closed-form of ${\mathrm{OP}}_{{\mathrm{D}}_{\mathit{n}}}$ can be simplified as follows.

$$\begin{array}{cc}\hfill {\mathrm{OP}}_{{\mathrm{D}}_{\mathit{n}}}\approx & \frac{N!}{(N-n)!(n-1)!}\sum _{j=0}^{N-n}\frac{{(-1)}^j}{n+j}\left(\genfrac{}{}{0pt}{}{N-n}{j}\right)\sum _{\ell =0}^{n+j}{(-1)}^{\ell }\left(\genfrac{}{}{0pt}{}{n+j}{\ell }\right)\frac{1}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}\hfill \\ \hfill & \times {\int }_0^{\infty }{x}^{K-1}exp\left(-\frac{\ell v}{x(1-\rho ){\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}-\frac{x}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right)dx.\hfill \end{array}$$

(38)

Based on ([32], Equation (3.471.9)), an approximation of ${\mathcal{I}}_2$ is evaluated as

$$\begin{array}{cc}\hfill {\mathrm{OP}}_{{\mathrm{D}}_{\mathit{n}}}\approx & \frac{N!}{(N-n)!(n-1)!}\sum _{j=0}^{N-n}\frac{{(-1)}^j}{n+j}\left(\genfrac{}{}{0pt}{}{N-n}{j}\right)\sum _{\ell =0}^{n+j}{(-1)}^{\ell }\left(\genfrac{}{}{0pt}{}{n+j}{\ell }\right)\frac{1}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}\hfill \\ \hfill & \times 2{\left(\frac{\ell v{\mathsf{\Omega }}_{\mathrm{SR}}}{{\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}\right)}^{\frac{K}{2}}{\mathcal{K}}_K\left(\sqrt{\frac{4\ell v}{{(1-\rho )}^2{\mathsf{\Omega }}_{\mathrm{SR}}{\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}}\right),\hfill \end{array}$$

(39)

where ${\mathcal{K}}_K$ expresses the modified Bessel function of the second kind of order K.

On the other hand, in the case of imperfect SIC, (21) can be rewritten as

$$\begin{array}{c}\hfill {\mathrm{OP}}_{{\mathrm{D}}_{\mathit{n}}}^{\mathrm{ipSIC}}\stackrel{\mathsf{\Delta }}{=}1-Pr\left(\begin{array}{c}X\ge \frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{R}}^2}{P_{\mathrm{S}}[a_n-{\gamma }_{\mathrm{th}}({\beta }_1+{\beta }_2)]},\\ XY\ge \frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{D},n}^2}{P_{\mathrm{S}}\varphi [b_n-{\gamma }_{\mathrm{th}}({\psi }_1+{\psi }_2)]}\end{array}\right),\end{array}$$

(40)

where ${\beta }_1={\sum }_{i=n+1}^N{a}_i$, ${\beta }_2={\xi }_1{\sum }_{k=1}^{n-1}{a}_k$ and ${\psi }_1={\sum }_{i=n+1}^N{b}_i$, ${\psi }_2={\xi }_2{\sum }_{k=1}^{n-1}{b}_k$.

Similarly to the case of perfect SIC, based on the theory of joint probability of two random variable function, we have the OP under imperfect SIC as

$$\begin{array}{cc}\hfill {\mathrm{OP}}_{{\mathrm{D}}_{\mathit{n}}}^{\mathrm{ipSIC}}& =1-{\int }_{\mu }^{\infty }\left[1-F_Y\left(\frac{\lambda }{x}\right)\right]f_X\left(x\right)dx\hfill \\ \hfill & =1-\underset{{\mathcal{J}}_1}{\underbrace{{\int }_{\mu }^{\infty }{f}_X\left(x\right)dx}}+\underset{{\mathcal{J}}_2}{\underbrace{{\int }_{\mu }^{\infty }{F}_Y\left(\frac{\lambda }{x}\right)f_X\left(x\right)dx}},\hfill \end{array}$$

(41)

where $\mu =\frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{R}}^2}{P_{\mathrm{S}}[a_n-{\gamma }_{\mathrm{th}}({\beta }_1+{\beta }_2)]}$ and $\lambda =\frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{D},n}^2}{P_{\mathrm{S}}\varphi [b_n-{\gamma }_{\mathrm{th}}({\psi }_1+{\psi }_2)]}$.

Similar to the case of perfect SIC and based on CDF and PDF given in (22) and (23), the result of imperfect SIC is described as follows.

$$\begin{array}{c}\hfill {\mathcal{J}}_1=\frac{1}{\mathsf{\Gamma }\left(K\right)}\mathsf{\Gamma }\left(K,\frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{R}}^2}{P_{\mathrm{S}}[a_n-{\gamma }_{\mathrm{th}}({\beta }_1+{\beta }_2)](1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right).\end{array}$$

(42)

Besides, ${\mathcal{J}}_2$ can be obtained as follows.

$$\begin{array}{cc}\hfill {\mathcal{J}}_2=& \frac{N!}{(N-n)!(n-1)!}\sum _{j=0}^{N-n}\frac{{(-1)}^j}{n+j}\left(\genfrac{}{}{0pt}{}{N-n}{j}\right)\sum _{\ell =0}^{n+j}{(-1)}^{\ell }\left(\genfrac{}{}{0pt}{}{n+j}{\ell }\right)\frac{1}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}\hfill \\ \hfill & \times \sum _{t=0}^{\infty }\frac{{(-1)}^t}{t!}{\left(\frac{\ell \lambda }{(1-\rho ){\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}\right)}^t\underset{{\mathsf{\Delta }}_2}{\underbrace{{\int }_{\mu }^{\infty }{x}^{q}exp\left(-\frac{x}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right)dx}},\mathrm{where}q=K-1-t.\hfill \end{array}$$

(43)

By replacing the inner parameters for the corresponding variables, we have ${\mathsf{\Delta }}_2$ as in (44) and (45).

$$\begin{array}{cc}\hfill & {\mathsf{\Delta }}_2={\left(\frac{1}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}}\right)}^{-q-1}\mathsf{\Gamma }\left(q+1,\frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{R}}^2}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}{P}_{\mathrm{S}}[a_n-{\gamma }_{\mathrm{th}}({\beta }_1+{\beta }_2)]}\right),\mathrm{if} q\ge 0.\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill & {\mathsf{\Delta }}_2={\left(\frac{P_{\mathrm{S}}[a_n-{\gamma }_{\mathrm{th}}({\beta }_1+{\beta }_2)]}{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{R}}^2}\right)}^{q-1}{E}_q\left(\frac{{\gamma }_{\mathrm{th}}{\sigma }_{\mathrm{R}}^2}{(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}{P}_{\mathrm{S}}[a_n-{\gamma }_{\mathrm{th}}({\beta }_1+{\beta }_2)]}\right),\mathrm{if} q<0\hfill \end{array}$$

(45)

3.2. Ergodic Rate

Before investigating the ergodic rate of the system, we examine the result that had been shown in [44], i.e., the transmission power of the relay is often smaller than the transmission power of the BS because the relay node is supplied by harvesting energy. In Figure 4, the probabilities of ${\gamma }_{\mathrm{SR}}$ and ${\gamma }_{\mathrm{RD}}$ are compared. It is clear that ${\gamma }_{\mathrm{RD}}$ is always less than ${\gamma }_{\mathrm{SR}}$. Therefore, the probability that ${\gamma }_{{\mathrm{RD}}_{\mathit{n}}}<{\gamma }_{\mathrm{th}}$ is always higher than the probability that ${\gamma }_{\mathrm{SR}}<{\gamma }_{\mathrm{th}}$. This ensures that $min({\gamma }_{\mathrm{SR}},{\gamma }_{{\mathrm{RD}}_{\mathit{n}}})={\gamma }_{{\mathrm{RD}}_{\mathit{n}}}$ and it can be considered to be a crucial factor that limits the performance of such systems.

Based on the above analysis and Claude Shannon theory, the instantaneous data rate of $(m,n)$th user is given by

$$\begin{array}{cc}\hfill R_{m,n}=& \frac{1-\alpha }{2}{log}_2\left(1+\frac{P_{\mathrm{R}}{b}_n{\left|g_n\right|}^2}{\tilde{b}{P}_{\mathrm{R}}|g_n{|}^2+{\psi }_2{P}_{\mathrm{R}}{\left|g_n\right|}^2+{\sigma }_{\mathrm{D},n}^2}\right)\hfill \\ \hfill =& \frac{1-\alpha }{2}{log}_2\left(\frac{P_{\mathrm{R}}{\left|g_n\right|}^2(\tilde{b}+{\psi }_2+b_n)+{\sigma }_{\mathrm{D},n}^2}{\tilde{b}{P}_{\mathrm{R}}|g_n{|}^2+{\psi }_2{P}_{\mathrm{R}}{\left|g_n\right|}^2+{\sigma }_{\mathrm{D},n}^2}\right).\hfill \end{array}$$

(46)

We can rewrite (46) as

$$\begin{array}{cc}\hfill R_{m,n}=& \frac{1-\alpha }{2}{log}_2\left(1+\varphi \mathsf{\Psi }|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2{\left|g_n\right|}^2(\tilde{b}+{\psi }_2+b_n)\right)-\frac{1-\alpha }{2}{log}_2\left(1+\varphi \mathsf{\Psi }|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2{\left|g_n\right|}^2(\tilde{b}+{\psi }_2)\right),\hfill \end{array}$$

(47)

where $\mathsf{\Psi }=\frac{P_{\mathrm{S}}}{{\sigma }_{{\mathrm{D}}_{\mathit{n}}}^2}$ is the received signal to noise ratio.

To simplify, we let ${\mathsf{\Gamma }}_1=\varphi \mathsf{\Psi }|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2{\left|g_n\right|}^2(\tilde{b}+{\psi }_2+b_n)$ and ${\mathsf{\Gamma }}_2=\varphi \mathsf{\Psi }|{\mathbf{h}}_m{\mathbf{w}}_m{|}^2{\left|g_n\right|}^2(\tilde{b}+{\psi }_2)$. Thus, we have CDFs of ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ given as

$$\begin{array}{cc}\hfill F_{{\mathsf{\Gamma }}_1}\left({\gamma }_1\right)=1-& \sum _{n=1}^N\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{2{(-1)}^{n-1}}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}{\left(\frac{n{\gamma }_1{\mathsf{\Omega }}_{\mathrm{SR}}}{{\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}\varphi \mathsf{\Psi }\mathcal{B}}\right)}^{\frac{K}{2}}{\mathcal{K}}_K\left(\sqrt{\frac{4n{\gamma }_1}{{(1-\rho )}^2{\mathsf{\Omega }}_{\mathrm{SR}}\varphi \mathsf{\Psi }\mathcal{B}{\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}}\right),\hfill \end{array}$$

(48)

and

$$\begin{array}{cc}\hfill F_{{\mathsf{\Gamma }}_2}\left({\gamma }_2\right)=1-& \sum _{n=1}^N\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{2{(-1)}^{n-1}}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}{\left(\frac{n{\gamma }_2{\mathsf{\Omega }}_{\mathrm{SR}}}{{\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}\varphi \mathsf{\Psi }\mathcal{C}}\right)}^{\frac{K}{2}}{\mathcal{K}}_K\left(\sqrt{\frac{4n{\gamma }_2}{{(1-\rho )}^2{\mathsf{\Omega }}_{\mathrm{SR}}\varphi \mathsf{\Psi }\mathcal{C}{\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}}\right),\hfill \end{array}$$

(49)

where $\mathcal{B}=(\tilde{b}+b_n+{\psi }_2)$, $\mathcal{C}=(\tilde{b}+{\psi }_2)$. $F_{{\mathsf{\Gamma }}_1}\left(\gamma \right)$ and $F_{{\mathsf{\Gamma }}_2}\left(\gamma \right)$ are solved in Appendix A.

Based on (47), the ergodic rate can be calculated as

$$\begin{array}{cc}\hfill \mathbb{E}\left[R_{m,n}\right]=& \underset{{\mathrm{C}}_1}{\underbrace{\frac{1-\alpha }{2}\mathbb{E}\left[{log}_2\left(1+{\mathsf{\Gamma }}_1\right)\right]}}-\underset{{\mathrm{C}}_2}{\underbrace{\frac{1-\alpha }{2}\mathbb{E}\left[{log}_2\left(1+{\mathsf{\Gamma }}_2\right)\right]}}.\hfill \end{array}$$

(50)

Expression (50) follows the strictly monotonically increasing property of the logarithm function for non-negative real numbers.

$$\begin{array}{cc}\hfill {\mathrm{C}}_1& =\frac{1-\alpha }{2ln2}\underset{0}{\overset{\infty }{\int }}{log}_2(1+{\mathsf{\Gamma }}_1)f_{{\mathsf{\Gamma }}_1}\left({\gamma }_1\right)d{\gamma }_1.\hfill \end{array}$$

(51)

Based on the properties of the expected value of a random variable, using the integration-by-parts method, we have

$$\begin{array}{cc}\hfill {\mathrm{C}}_1& =\frac{1-\alpha }{2ln2}\underset{0}{\overset{b}{\int }}{log}_2(1+{\mathsf{\Gamma }}_1)d[F_{{\mathsf{\Gamma }}_1}\left({\gamma }_1\right)-1]\hfill \\ \hfill & =\frac{1-\alpha }{2ln2}\left({log}_2(1+{\mathsf{\Gamma }}_1)[F_{{\mathsf{\Gamma }}_1}\left(b\right)-1]\right)-\frac{1-\alpha }{2ln2}\underset{0}{\overset{b}{\int }}\frac{F_{{\mathsf{\Gamma }}_1}\left({\gamma }_1\right)-1}{1+{\gamma }_1}d{\gamma }_1.\hfill \end{array}$$

(52)

$F_{{\mathsf{\Gamma }}_1}\left({\gamma }_1\right)$ is CDF of ${\mathsf{\Gamma }}_1$, thus when b is a finite value, it is seen that ${log}_2(1+{\mathsf{\Gamma }}_1)[F_{{\mathsf{\Gamma }}_1}\left(b\right)-1]=0$. Thus, ${\mathrm{C}}_1$ is given as

$$\begin{array}{cc}\hfill {\mathrm{C}}_1& =\frac{1-\alpha }{2ln2}\underset{0}{\overset{\infty }{\int }}\frac{1-F_{{\mathsf{\Gamma }}_1}\left({\gamma }_1\right)}{1+{\gamma }_1}d{\gamma }_1.\hfill \end{array}$$

(53)

Similar to the working step in (51), we can rewrite the second part as the following equation.

$$\begin{array}{cc}\hfill {\mathrm{C}}_2& =\frac{1-\alpha }{2ln2}\underset{0}{\overset{\infty }{\int }}{log}_2(1+{\mathsf{\Gamma }}_2)f_{{\mathsf{\Gamma }}_2}\left({\gamma }_2\right)d{\gamma }_2\hfill \\ \hfill & =\frac{1-\alpha }{2ln2}\underset{0}{\overset{\infty }{\int }}\frac{1-F_{{\mathsf{\Gamma }}_2}\left({\gamma }_1\right)}{1+{\gamma }_2}d{\gamma }_2,\hfill \end{array}$$

(54)

where the $F_{{\mathsf{\Gamma }}_1}\left(\gamma \right)$ and $F_{{\mathsf{\Gamma }}_2}\left(\gamma \right)$ are given in (48) and (49).

With the help of ([32], Equation (7.811.5)) and ([32], Equation (9.343)), after some manipulations, we have ${\mathrm{C}}_1$ and ${\mathrm{C}}_2$ as

$$\begin{array}{cc}\hfill {\mathrm{C}}_1=& \frac{1-\alpha }{2ln2}\sum _{n=1}^N\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{{(-1)}^{n-1}}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}{\left(\frac{n{\mathsf{\Omega }}_{\mathrm{SR}}}{{\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}\varphi \mathsf{\Psi }\mathcal{B}}\right)}^{\frac{K}{2}}{\mathrm{G}}_{1,3}^{3,1}\left({\left.\frac{n}{{(1-\rho )}^2{\mathsf{\Omega }}_{\mathrm{SR}}\varphi \mathsf{\Psi }\mathcal{B}{\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}\right|}_{-\frac{K}{2},-\frac{K}{2},\frac{K}{2}}^{-\frac{K}{2}}\right),\hfill \end{array}$$

(55)

and

$$\begin{array}{cc}\hfill {\mathrm{C}}_2=& \frac{1-\alpha }{2ln2}\sum _{n=1}^N\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{{(-1)}^{n-1}}{\mathsf{\Gamma }\left(K\right){\left[(1-\rho ){\mathsf{\Omega }}_{\mathrm{SR}}\right]}^K}{\left(\frac{n{\mathsf{\Omega }}_{\mathrm{SR}}}{{\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}\varphi \mathsf{\Psi }\mathcal{C}}\right)}^{\frac{K}{2}}{\mathrm{G}}_{1,3}^{3,1}\left({\left.\frac{n}{{(1-\rho )}^2{\mathsf{\Omega }}_{\mathrm{SR}}\varphi \mathsf{\Psi }\mathcal{C}{\mathsf{\Omega }}_{{\mathrm{RD}}_{\mathit{j}}}}\right|}_{-\frac{K}{2},-\frac{K}{2},\frac{K}{2}}^{-\frac{K}{2}}\right).\hfill \end{array}$$

(56)

The proof is shown in Appendix A.

4. Numerical Results

In this section, the performance of the proposed SWIPT-NOMA system in terms of the outage probability and ergodic rate is numerically analyzed. The simulation parameters are set as follows. The Monte Carlo simulations is run with $2\times {10}^{14}$ trials. The variance of noise component at all receiving nodes is assumed to be constant, ${\sigma }^2=1$. The average channel gain $\mathbb{E}\left\{\right|{\mathbf{h}}_m{|}^2\}={\mathsf{\Omega }}_{{\mathrm{SR}}_m}=\mathbb{E}\left\{\right|g_1{|}^2\}={\mathsf{\Omega }}_{{\mathrm{R}}_m{\mathrm{D}}_1}=1$, $\mathbb{E}\left\{\right|g_2{|}^2\}={\mathsf{\Omega }}_{{\mathrm{R}}_m{\mathrm{D}}_2}=2$ and $\mathbb{E}\left\{\right|g_3{|}^2\}={\mathsf{\Omega }}_{{\mathrm{R}}_m{\mathrm{D}}_3}=4$.

The threshold data rates of ${\mathrm{D}}_n$: $r_1=r_2=r_3=1[\mathrm{b}/\mathrm{s}/\mathrm{Hz}]$. The energy conversion efficiency coefficient of $\mathrm{R}$: $\eta =0.85$. The number of relay nodes: $M=3$ and each cluster has three users severed instantaneously by each relay node. The power allocation coefficient of the nth user: $a_n=(N-n+1)/\mu$, where $\mu$ is chosen such that ${\sum }_{n=1}^N\sqrt{a_n{P}_S}=1$. To simplify the system design and settings, we assume that the power allocation coefficients at the BS and the relay nodes are the same.

Firstly, users are grouped randomly into clusters according to Poisson distribution in two cases of the $\lambda$ parameter values, $\lambda$ can be considered to be an expected value of the number of such events, in which an event is that a user appears in a cluster. Only three users with the best channel conditions are supported by the relay node at the same time, and the OPs of every user in both cases of $\lambda$ are illustrated in Figure 5. As can be seen from this figure, if the transmission power is high enough, the OPs evaluated by exact analysis and approximation are the same, demonstrating that the proposed methods to calculate the OP are significantly reasonable and actually accurate. Moreover, while the transmission power is assigned for signal of user 1 with the highest value and user 3 with the lowest value, the OP of user 3 is the lowest, in other words, user 3 has the best quality. The reason is that the total interference affected to user 1 consists of signals of both user 2 and user 3, and its own noise, while the total interference influenced on user 3 contains its own noise only owing to the perfect SIC. Furthermore, the system performance is improved if the $\lambda$ parameter increases. The reason is explained as when the $\lambda$ rises, the number of users appears in a cluster goes up, resulting in the selection of three users with the best channel conditions is more efficient.

In Figure 6, the proposed system is investigated in such scenario that users are fixed and randomly grouped according to Poison rule with the average value of the $\lambda$ parameter. It is clear that when the number of transmitting antennas of the BS increases, the system performance is enhanced. There are two reasons, i.e., the diversity gain of the link between BS-R is improved and the harvested energy at the relay node rises or transmission power of the relay node rises and then the messages are sent from the relay node more dependably. In addition, the system performance in case of random user distribution may be better than that in case of fixed user distribution. The reason is that if users are fixed, the channel gains are likely to be fixed and then the selection of three users with best channel condition is less effective than that in case of random distribution. Finally, the results in exact analysis and approximation are absolutely the same at high enough SNR region, once again, confirming the accuracy of the proposed methods.

The effect of channel estimation error on the OP of the system is illustrated in Figure 7. The curves show that increasing $\rho$, i.e., increasing estimation error, results in a reduction of the OP performance. This result demonstrates that the channel estimation error influences significantly on the system performance. In fact, the instantaneous CSI not only affects to coding gain but also causes leakage beam also. Finally, the Monte Carlo simulation guarantees the correctness of the exact and approximation analytical results.

In Figure 8, the OP of users versus SNR in two cases of perfect and imperfect SICs are represented. In this figure, the OP of imperfect SIC is higher than that of perfect SIC and the OP of user 1 are the same in two cases. In fact, user 1 does not employ the SIC structure, it decodes its own signal by considering the other users’ signals as interference. However, the residual power of user 1 impacts to the decoding process of the other users. Different from the case of perfect SIC, the quality of user 3 is least of all among users in the scenario of imperfect SIC. It is because the residual power, that plays as interference, at user 3, it is the largest. The provided simulation results confirm the accuracy of the proposed analytical method, and the approximation is the same with simulation in the medium and high SINR regions.

Figure 9 demonstrates the impact of the time fraction of EH, $0<\alpha <1$, on the OP of users. In fact, the optimal $\alpha$ is a sophisticated function of the channel and system parameters. The smaller $\alpha$ results in the smaller energy for relay nodes, however, the larger $\alpha$ brings the smaller transmission power for communication of the BS. Therefore, if the $\alpha$ increases from 0 to 1, the OP reduces to an optimal value and then increases again. As a result, the optimal $\alpha$ corresponding to the smallest OP is determined by this plot. In the figure, the minimum values of OP for every user are different although the setting parameters are the same. Moreover, the optimal time duration for EH of user 3 is the largest. The reason is that the power allocation for user 3 is the smallest, consequently, it needs longer time duration for EH to ensure the transmission power of the relay node to forward signals.

Figure 10 demonstrates the ER of users versus SINR in the case of perfect SIC. In this figure, the ER can be considered to be a function of the mean SINR, based on the derived results in (55) and (56). Figure 10 shows that the ERs of users 1 and 2 increase negligibly in the low SINR region and remain stable in the high SINR region; however, that of user 3 increases exponentially. There is a trade-off between the complexity and the performance, i.e., user 1 directly detected its own signal, whereas the users 2 and 3 have the first-order and second-order of the SIC, respectively.

The sum ER of the system is demonstrated in Figure 11 for two cases, i.e., perfect and imperfect SICs. Firstly, the simulation results are the same with analytical results. Secondly, the sum ER in the case of perfect SIC is higher than that in the case of imperfect SIC, especially the gap of them increases when the SINR increases. The reason can be explained as if SINR increases by rising the transmission power, the interference due to imperfect SIC also increases, consequently, the sum ER with imperfect SIC goes up less significantly than that with perfect SIC. Finally, the ERs of user 1 are the same in both cases because when user 1 decodes the signal, the SIC scheme is not used.

Figure 12 demonstrates the sum ER of the system versus SNR with the following settings. The number of users is 12, there are three typical user clustering modes, and the imperfect SIC coefficients are ${\xi }_1={\xi }_2=0.06$. As shown in Figure 12, the first mode has the best performance due to the lowest cross user interference. In contrast, the third mode achieves the worst performance because when the number of users in the cluster increases, the cross user interference rises. It is worth noting that the number of users and the number of clusters determine the number of beams and the order of SIC. As a result, we can dynamically select the user clustering mode according to system parameters and channel conditions for balancing the overall performance and the computational complexity. In addition, the proposed multiple relay aided MIMO-NOMA scheme and a MIMO-NOMA scheme without relay are compared in Figure 12. It is obvious that the proposed scheme has much better performance than the scheme without relay. In fact, the scheme without relay has lower channel gains than the proposed scheme, i.e., ${\mathsf{\Omega }}_{\mathrm{BS}-{\mathrm{D}}_1}=0.5,{\mathsf{\Omega }}_{\mathrm{BS}-{\mathrm{D}}_2}=1$ and ${\mathsf{\Omega }}_{\mathrm{BS}-{\mathrm{D}}_3}=2$, because of the longer distance from point to point (BS, R or D). Moreover, in the proposed scheme, the inter-relay interference can be effectively canceled.

5. Conclusions

In this paper, a combination design of NOMA, SWIPT, and beamforming is proposed for the downlink of multi-user relaying systems. To analyze and evaluate the system performance, the outage probability and ergodic capacity are derived based on analytics and statistics. The proposed system is investigated obviously in two cases of the perfect and imperfect SIC structures, and the effect of the imperfect CSI is also taken into consideration. Moreover, an optimal time fraction of energy harvesting is found out to minimize the outage probability and the reasonable user clustering mode is discussed. The analytical and simulation results demonstrate the accuracy of the proposed analysis method. In the future, the authors will investigate the system with a full-duplex model to improve the spectrum efficiency and apply the power splitting scheme to the SWIPT with the finite block-length.