On the Performance of Energy Harvesting Non-Orthogonal Multiple Access Relaying System with Imperfect Channel State Information over Rayleigh Fading Channels

Abstract

In this paper, we propose a non-orthogonal multiple access (NOMA) relaying system, where a source node communicates simultaneously with multiple users via the assistance of the best amplify-and-forward (AF) relay. The best relay is selected among N relays which are capable of harvesting the energy from radio frequency (RF) signals. We analyze the performance of the proposed NOMA relaying system in the conditions of imperfect channel state information (CSI) and Rayleigh fading by deriving the exact expressions of the outage probability (OP) and the approximate expression of the ergodic capacities of each user and the whole system. We also determine the optimal energy harvesting duration which minimizes the OP. Numerical results show that, for the same parameter settings, the performance of the proposed NOMA relaying system, especially the ergodic capacity of the whole system, outperforms that of the orthogonal-multiple-access (OMA) relaying system. Monte-Carlo simulations are used to validate the correctness of the analytical results.

1. Introduction

Nowadays, the Internet of Things (IoT) has received increasing attention from both industry and academia. It is considered an important mean for wireless connections in the fourth industrial revolution. IoT is also being used in the fourth generation (4G) mobile communications and will be applied to the fifth generation (5G). In order to support a large multiuser system such as IoT, the non-orthogonal multiple access (NOMA) is a very potential technique due to its high bandwidth efficiency [1,2]. Moreover, compared with conventional orthogonal multiple access (OMA) systems, such as time division multiple access (TDMA), code division multiple access (CDMA), orthogonal frequency division multiple access (OFDMA), the NOMA systems offer better fairness among users, even for users with weak channel conditions such as the cell-edge users. The fundamentals of the NOMA system can be found in Reference [3] while a study of NOMA system in cellular communication with machine-to-machine in IoT is given in Reference [4].

Recently, the power supply for terminal devices in wireless networks has become an important matter and has attracted much interest from researchers. Besides using the optimal power allocation for the fifth generation (5G) and sixth generation (6G) networks (6G network will start to enter the market by 2026 [5]) to reduce the power consumption [6], another promising method to improve the lifetime of communication devices is to generate electric power from some external energy sources such as solar, wind, and radio frequency (RF) signal to charge the batteries. Unfortunately, natural energy sources are not suitable for small-size mobile devices and in some cases they cannot be used in the healthcare monitoring networks and the sensor networks with real-time requirements. In contrast, the RF energy is often available due to its increased power density and availability, and is independent on environmental conditions, including weather, climate, and temperature. As the result, the RF energy harvesting (EH), also called simultaneous wireless information and power transfer (SWIPT), has been widely used compared with other kinds of energy harvesting techniques [7,8,9]. SWIPT has been applied not only in the point-to-point systems but also in relaying systems because deploying relays can improve the amount of harvested energy and the coverage area of wireless networks. The authors of References [10,11] investigated information and energy receiver architecture for SWIPT networks. Reference [11] especially considered a non-linear energy harvesting model which described the practical system well.

To prolong network lifetime and improve the spectral utilizing efficiency, NOMA is combined with SWIPT [12]. In Reference [13], the authors investigated the tradeoff among the energy efficiency, fairness, harvested energy, and system sum rate of NOMA systems in power domain. Investigation of an integrated wireless communication system including NOMA, full-duplex relaying, and energy harvesting techniques was conducted in Reference [14]. The authors of References [15,16] studied the system performance of cooperative NOMA systems and derived the expressions of outage probability in the conditions of perfect successive interference cancellation (SIC) and perfect channel state information (CSI). In Reference [17], the near users which are close to the base station will harvest the RF energy and forward signals to far users. The analysis results showed that if the time switching ratio in NOMA system with SWIPT is appropriately chosen, the diversity gain will not be impaired. The authors of Reference [18] proposed a NOMA system where source node communicates with two users via the assistance of the best relay with the RF energy harvesting capability. The exact expressions of the outage probability and throughput were used as the criteria to evaluate the system performance. The effects of power allocation and time switching ratio on the performance of multi-user NOMA system were investigated in Reference [19]. Specifically, the authors derived the outage probability expression and determined the optimal power allocation coefficient for two NOMA power allocation policies, namely NOMA with fixed power allocation (F-NOMA) and cognitive radio inspired NOMA (CR-NOMA). It was shown that when a reasonable power allocation coefficient is selected, higher system performance can be achieved in comparison with the conventional multi-user system.

We observe that all previous works only mentioned the case of perfect CSI and used only one relay to forward signals to multiple users. Moreover, although the partial relay selection has been widely studied in conventional wireless systems, it has not been analyzed in NOMA systems. Another observation is that the NOMA systems perform superimposing signals in power domain, thus they always require CSI to allocate power for all users. However, due to variation in the communication quality of wireless environment, the imperfect CSI may happen [20,21]. Perfect CSI exists if and only if the amount of feedback CSI from users to the base station is large and the length of the pilot sequences which are used to estimate channel is very long. Unfortunately, these conditions rarely happen in practice. Therefore, investigation of the impact of imperfect CSI on the relay selection and power allocation is vitally important to the design of practical NOMA systems.

Motivated by the above issues, in this paper we propose a downlink NOMA relaying system with partial relaying selection. In this system, source node transmits superposition modulated signals to multiple users via the assistance of the best relay. The best relay is chosen from a set of relays which are capable of harvesting RF energy and grouped by their locations. Based on the feedback CSI from all users, the source node performs power allocation and chooses the best communication link. The main contributions of this paper can be summarized as follows:

• We overcome the limitation of current multiple access techniques and the energy demand of wireless networks by proposing the downlink NOMA relaying system where the best relay is selected from a set of multiple RF energy harvesting relays.
• We study the system performance in terms of the outage probability and the ergodic capacity of each user and the whole system in the condition of imperfect CSI and Rayleigh fading. The imperfection of the CSI is modeled by the correlation coefficient and its impact on the system performance is investigated by using both analysis and simulation approaches. We also compare the outage performance and the ergodic capacity of the proposed NOMA relaying system with those of OMA relaying system.
• We determine the optimal time switching ratio to balance between the energy harvesting and the signal processing so that the outage probability can be minimized. All analysis results are validated by simulation results.

The rest of this paper is organized as follows. Section 2 describes the proposed downlink NOMA relaying system with partial relaying selection and time switching (TS) protocol. The analysis of the outage probability and ergodic capacity of the proposed system are presented in Section 3 and Section 4, respectively. Section 5 shows numerical results to evaluate the system performance. Finally, the conclusions are given in Section 6.

For the sake of clarity, the frequently used mathematical notations together with their descriptions are summarized in

Table 1. The mathematical notations used in this paper.

Notation	Description
$F_U\left(u\right)$	Cumulative distribution function (CDF)
$f_U\left(u\right)$	Probability density function (PDF)
$\mathcal{CN}(\mu ,{\sigma }^2)$	Circularly symmetric complex Gaussian distribution X with mean $\mu$ and variance ${\sigma }^2$
${\gamma }_{\mathrm{th}}$	Predefined outage threshold
$\mathbb{E}\left\{·\right\}$	Expectation operator
${\mathcal{K}}_n(·)$	Second order Bessel function n [22]
$I_0(·)$	Modified zero order Bessel function of first kind [22]
$\alpha$	Time switching ratio
$\eta$	Energy conversion efficiency
$\rho$	Channel correlation coefficient
T	Transmission period

.

2. System Model

Figure 1 illustrates the proposed downlink NOMA relaying system. In this system, source node $\mathrm{S}$ transmits the signals which are coded and superposed in power domain to multiple users ${\mathrm{D}}_m$, $m\in \{1,\cdots ,M\}$, via the assistance of the best relay which is selected from a set of relays ${\mathrm{R}}_n$, $n\in \{1,\cdots ,N\}$. The direct link $\mathrm{S}$-${\mathrm{D}}_m$ is assumed not available because the distance between $\mathrm{S}$ and ${\mathrm{D}}_m$ is larger than the coverage area of $\mathrm{S}$ or due to deep shadow fading.

We consider that all nodes are equipped with single antenna and operate in half-duplex mode. All channels between $\mathrm{S}$ and ${\mathrm{R}}_n$ and between ${\mathrm{R}}_n$ and ${\mathrm{D}}_m$ are influenced by block fading, that is, the symbol rate is larger than channel varying rate so that it can be considered as constant over each symbol duration. The communication links from $\mathrm{S}$ to ${\mathrm{R}}_n$ and from ${\mathrm{R}}_n$ to ${\mathrm{D}}_m$ are respectively modeled as complex Gaussian distributions with zero mean and variances ${\mathsf{\Omega }}_{1,n}$ and ${\mathsf{\Omega }}_{{\mathrm{R}}_n{\mathrm{D}}_m}$, that is, $h_{1,n}\sim \mathcal{CN}(0,{\mathsf{\Omega }}_{1,n})$ and $g_m\sim \mathcal{CN}(0,{\mathsf{\Omega }}_{{\mathrm{R}}_n{\mathrm{D}}_m})$. The Additive White Gaussian Noise (AWGN) at the relays and users are $w_{{\mathrm{R}}_n}\sim \mathcal{CN}(0,{\sigma }_{{\mathrm{R}}_n}^2)$ and $w_{{\mathrm{D}}_m}\sim \mathcal{CN}(0,{\sigma }_{{\mathrm{D}}_m}^2)$, respectively. Because of the time varying characteristic of wireless channel, its coherent time may be altered when the feedback delay is larger than the transmission block period of a symbol. Thus, the received CSI is always outdated at the transmitter, which often happens in practice [23,24].

We denote ${\rho }_i$, $0\le {\rho }_i\le 1$, $i\in \left\{1,2\right\}$, as the correlation coefficients between the past channel $h_{1,n}$ and the current channel ${\widehat{h}}_{1,n}$, similarly for ${\widehat{g}}_m$ and $g_m$. These coefficients can be considered as the measurements of the fluctuation rate of wireless channels and are related solely to the time delay. Based on the Markov chain, the relationship between ${\widehat{h}}_{1,n}$ and $h_{1,n}$ and between ${\widehat{g}}_m$ and $g_m$ can be presented as [25]

$$\begin{array}{c}\hfill {\widehat{h}}_{1,n}={\rho }_1{h}_{1,n}+\sqrt{1-{\rho }_1^2}{\epsilon }_{1,n},\end{array}$$

(1)

$$\begin{array}{c}\hfill {\widehat{g}}_m={\rho }_2{g}_m+\sqrt{1-{\rho }_2^2}{\epsilon }_m,\end{array}$$

(2)

where ${\epsilon }_{1,n}$ and ${\epsilon }_m$ are the circular symmetric complex Gaussian random variables which can be modeled as ${\epsilon }_{1,n}\sim \mathcal{C}\mathcal{N}(0,{\sigma }^2)$ and ${\epsilon }_m\sim \mathcal{C}\mathcal{N}(0,{\sigma }^2)$, respectively.

In this paper, a partial relay selection (PRS) scheme [26] is used to select the best relay from a set of relays. According to the PRS scheme, $\mathrm{S}$ continuously monitors the gain of $\mathrm{S}$–${\mathrm{R}}_n$ channels by using the feedback signal and selects the communication link that gives the largest instantaneous channel gain, that is,

$$\begin{array}{c}\hfill {\gamma }_b=arg\underset{n=1,2\cdots ,N}{max}{\gamma }_{1,n},\end{array}$$

(3)

where ${\gamma }_{1,n}$ is the instantaneous SINR of $\mathrm{S}$-${\mathrm{R}}_n$ link.

The advantage of using PRS scheme is that the system configuration is simpler and easier than using the full relay selection (FRS) scheme [27]. In addition, the results in Reference [28] show that both PRS and FRS schemes have the same average channel capacity in a high SNR regime and the outage probability of PRS is higher than that of FRS when the number of relays is more than 2. On the other hand, FRS scheme may not applicable for multiple-user systems because the distances from the best relay to users are different, thus the calculation complexity of FRS scheme is significantly increased with the number of users.

After a link from the source node to the best relay is established, the transmission period T for communication process is spitted into two parts (in this system, we use the time-division multiple access (TDMA) scheme). According to the time switching (TS) protocol [29], a time duration $\alpha T$ is used for EH. The remaining time duration $(1-\alpha )T$ is divided into two equal time sub-slots, which are used for the information transmission. The first half $(1-\alpha )T/2$ is used for the information transmission from source node to the relay and the remaining half $(1-\alpha )T/2$ is used for the information transmission from the relay to the user. It should be noted that the case $\alpha =1$ is not considered in this paper because when the energy harvesting time duration takes the whole transmission period T, i.e., the relay does not process any signals, the basic role in signal forwarding of the relay is eliminated [30]. Hence, we only consider the time switching ratio $\alpha$ in the range $0\le \alpha <1$. Then, the harvested energy of the selected relay in time duration $\alpha T$ is expressed as [29]

$$\begin{array}{c}\hfill E_h=\alpha T\eta P_{\mathrm{S}}\underset{n=1,···,N}{\max }{\left|{\widehat{h}}_{1,n}\right|}^2,\end{array}$$

(4)

where $\eta$ is the energy conversion efficiency coefficient which varies from 0 to 1 and closely depends on the quality of energy harvesting electric circuitry, $P_{\mathrm{S}}$ is the transmission power of $\mathrm{S}$.

In our proposed relaying system, since the harvest-use (HU) architecture is used, the relay does not need an energy buffer to store the harvested energy. Since all amounts of harvested energy during EH phase is consumed by ${\mathrm{R}}_n$ for signal transmission from $\mathrm{R}$ to ${\mathrm{D}}_m$, from (4), the transmission power of the best relay is given by

$$\begin{array}{c}\hfill P_{\mathrm{R}}=\frac{E_h}{(1-\alpha )T/2}=\frac{2\alpha \eta P_{\mathrm{S}}}{1-\alpha }\underset{n=1,···,N}{\max }{\left|{\widehat{h}}_{1,n}\right|}^2.\end{array}$$

(5)

According to the NOMA technique in power domain, during the first time sub-slot $\frac{1-\alpha }{2}$, source node transmits the superimposed signal $x_{\mathrm{S}}={\sum }_{m=1}^M\sqrt{P_{\mathrm{S}}{a}_m}{x}_m$, where $x_m$ and $a_m$ are the signal and power allocation coefficient of mth user, respectively. At the end of this time sub-slot, the received signal at ${\mathrm{R}}_n$ is

$$\begin{array}{c}\hfill y_{\mathrm{R}}^n={\widehat{h}}_{1,b}{\sum }_{m=1}^M\sqrt{a_m{P}_{\mathrm{S}}}{x}_m+w_{{\mathrm{R}}_n},\end{array}$$

(6)

where ${\widehat{h}}_{1,b}=\underset{n=1,···,N}{\max }{\left|{\widehat{h}}_{1,n}\right|}^2$.

In the remaining second time slot $\frac{1-\alpha }{2}$, the relay employs the AF protocol to broadcast $y_{\mathrm{R}}^n$ to all users after multiplying it with an amplifying factor G. To keep the output power constraint at relay, it is required that $\mathrm{E}\{\parallel G{y}_{\mathrm{R}}^n{\parallel }^2\}=P_{\mathrm{R}}$, where $P_{\mathrm{R}}$ is given in (5), thus the amplifying factor G is given by

$$\begin{array}{c}\hfill G=\sqrt{\frac{2\alpha \eta P_{\mathrm{S}}{\left|{\widehat{h}}_{1,b}\right|}^2}{(1-\alpha )\left(P_{\mathrm{S}}\right|{\widehat{h}}_{1,b}{|}^2+{\sigma }_R^2)}}\approx \sqrt{\frac{2\alpha \eta }{(1-\alpha )}}.\end{array}$$

(7)

Therefore, the received signal at ${\mathrm{D}}_m$ in the case of perfect SIC is expressed as

$$\begin{array}{cc}\hfill y_{{\mathrm{D}}_m}=G{\widehat{h}}_{1,b}{\widehat{g}}_m\sqrt{a_m{P}_{\mathrm{S}}}{x}_m& +\underset{\mathrm{signals}\mathrm{of}\mathrm{other}\mathrm{users}}{\underbrace{G{\widehat{h}}_{1,b}{\widehat{g}}_m{\sum }_{j=m+1}^M\sqrt{a_j{P}_{\mathrm{S}}}{x}_j}}+\underset{\mathrm{noise}}{\underbrace{G{\widehat{g}}_m{w}_{{\mathrm{R}}_n}+w_{{\mathrm{D}}_m}}},\hfill \end{array}$$

(8)

where ${\widehat{g}}_m$ denotes the channel coefficient between ${\mathrm{R}}_n$ and ${\mathrm{D}}_m$.

The received signals at the best relay and each user is comprised of the desired signal and the signals of other users, which are treated as the interferences. Hence, to mitigate the negative effect of the inter-user interference, successive interference cancellation (SIC) method is applied.

For the downlink communication considered in this paper, the optimal SIC algorithm performs decoding signals in an order of increasing channel gain [31] ($|g_{{\mathrm{D}}_1}{|}^2\le |g_{{\mathrm{D}}_2}{|}^2\le \cdots \le |g_{{\mathrm{D}}_m}{|}^2\le {\left|g_{{\mathrm{D}}_M}\right|}^2$). To ensure the fairness among all users, the power allocation coefficients are assumed to be $a_1\ge a_2\ge \cdots a_m\ge a_M$, with ${\sum }_{m=1}^M{a}_m=1$. Hence, at the ${\mathrm{D}}_j$, the signal of ${\mathrm{D}}_m$, $j<m$, will be detected and then be removed from the received signal by SIC method. Specifically, ${\mathrm{D}}_j$ first decodes symbol $x_m$ while treating $x_j$ as noise.

Then, the SINR of symbol $x_m$ at ${\mathrm{D}}_j$ is given by

$$\begin{array}{c}\hfill {\gamma }_{m,j}^{\mathrm{D}}=\frac{G^2{a}_m{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_j\right|}^2}{G^2{\sum }_{j=m+1}^M{a}_j{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2|{\widehat{g}}_j{|}^2+G^2{\left|{\widehat{g}}_j\right|}^2{\sigma }_{\mathrm{R}}^2+{\sigma }_{{\mathrm{D}}_m}^2},\end{array}$$

(9)

where $j\in \left\{1,\dots ,m\right\}$ and $m\ne M$.

At ${\mathrm{D}}_j$, SIC will be performed until all signals of ${\mathrm{D}}_m$ are decoded successfully. Thus, the required SINR at ${\mathrm{D}}_m$ to successfully decode the signal by itself is given by

$$\begin{array}{c}\hfill {\gamma }_m^{\mathrm{D}}=\frac{G^2{a}_m{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2}{G^2{\sum }_{j=m+1}^M{a}_j{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2|{\widehat{g}}_m{|}^2+G^2{\left|{\widehat{g}}_m\right|}^2{\sigma }_{\mathrm{R}}^2+{\sigma }_{{\mathrm{D}}_m}^2}.\end{array}$$

(10)

We should note that the last user ${\mathrm{D}}_M$ needs to decode all signals of other users before decoding its signals. Consequently, the SINR for ${\mathrm{D}}_M$ to decode its own signals can be expressed as

$$\begin{array}{c}\hfill {\gamma }_M^{\mathrm{D}}=\frac{G^2{a}_M{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2{\left|g_M\right|}^2}{G^2{\left|g_M\right|}^2{\sigma }_{\mathrm{R}}+{\sigma }_{{\mathrm{D}}_M}^2}.\end{array}$$

(11)

3. Outage Probability Analysis

In this section, we derive the exact closed-form expression of the outage probability, taking into consideration the imperfect CSI and partial relay selection. It is well-known that the event that ${\mathrm{D}}_j$ can decode the signals of ${\mathrm{D}}_m$ successfully is

$$\begin{array}{c}\hfill {\mathsf{\Delta }}_{m,j}=\left\{\frac{G^2{a}_m{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_j\right|}^2}{G^2{\sum }_{j=m+1}^M{a}_j{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2|{\widehat{g}}_j{|}^2+G^2{\left|{\widehat{g}}_j\right|}^2{\sigma }_{\mathrm{R}}^2+{\sigma }_{{\mathrm{D}}_m}^2}>{\gamma }_{\mathrm{th}j}\right\},\end{array}$$

(12)

where ${\gamma }_{\mathrm{th}j}=2^{\frac{2r}{1-\alpha }}-1$ is the predefined outage threshold. This threshold is served as the protected value of the SINR to ensure the quality of service of the system and satisfy the target data rate r of ${\mathrm{D}}_j$.

Let us denote $X=|{\widehat{h}}_{1,b}{|}^2$ and $Z=|{\widehat{g}}_j{|}^2$. Without loss of generality, we assume that the temperature noise ${\sigma }_{\mathrm{R}}^2={\sigma }_{{\mathrm{D}}_m}^2={\sigma }^2$. Thus from (12), we can rewrite ${\mathsf{\Delta }}_{m,j}$ as

$$\begin{array}{c}\hfill {\mathsf{\Delta }}_{m,j}=\left\{\frac{G^2{a}_m{P}_{\mathrm{S}}XZ}{G^2{\sum }_{j=m+1}^M{a}_j{P}_{\mathrm{S}}XZ+G^{2}Z+{\sigma }^2}>{\gamma }_{\mathrm{th}j}\right\}.\end{array}$$

(13)

From (13) and after some manipulations, we can rewrite (13) as

$$\begin{array}{c}\hfill {\mathsf{\Delta }}_{m,j}\stackrel{\left(\lambda \right)}{=}\left\{X>{\theta }_j,Z>\frac{{\theta }_j}{G^2(X-{\theta }_j)}\right\},\end{array}$$

(14)

where ${\theta }_j=\frac{{\gamma }_{\mathrm{th}j}}{P_{\mathrm{S}}\left(a_m-{\sum }_{j=m+1}^M{a}_j{\gamma }_{\mathrm{th}j}\right)}$, step $\lambda$ holds when the condition $a_m>{\sum }_{j=1+m}^M{a}_j{\gamma }_{\mathrm{th}j}$ is satisfied. It should be noticed that ${\theta }_j=\frac{{\gamma }_{\mathrm{th}j}}{P_{\mathrm{S}}\left(a_m-{\sum }_{j=m+1}^M{a}_j{\gamma }_{\mathrm{th}j}\right)}$ is a constant and depends on the power allocation coefficient and the target data rate of ${\mathrm{D}}_j$.

The outage event occurs at ${\mathrm{D}}_j$ when it fails to decode its own signal or unsuccessfully performs SIC for the signals of ${\mathrm{D}}_m$ [32], i.e., ${\mathsf{\Lambda }}_{m,j}={\gamma }_{m,j}^{\mathrm{D}}<{\gamma }_{\mathrm{th}j},1\le j\le m$. Outage probability of the system occurs when the maximum SNRs at ${\mathrm{D}}_j$ falls below the threshold to decode signal. Thus, we have

$$\begin{array}{c}\hfill {\mathrm{P}}_{out}^j=Pr({\gamma }_{m,j}^{\mathrm{D}}\le {\gamma }_{\mathrm{th}j})=1-Pr({\gamma }_{m,j}^{\mathrm{D}}>{\gamma }_{\mathrm{th}j}),1\le j\le m.\end{array}$$

(15)

$$\begin{array}{c}\hfill {\mathrm{P}}_{out}^j=1-Pr\left\{{\mathsf{\Delta }}_{m,1}\cap {\mathsf{\Delta }}_{m,2}\cap \dots \cap {\mathsf{\Delta }}_{m,m}\right\},\end{array}$$

(16)

where ${\mathsf{\Delta }}_{m,j}$ is the complementary in the set of ${\mathsf{\Lambda }}_{m,j}$.

The condition in (14) always occurs, i.e., the outage probability is equal to one, if $a_m\le {\sum }_{j=1+m}^M{a}_j{\gamma }_{\mathrm{th}j}$. Hence, we need to allocate more power for ${\mathrm{D}}_m$ to satisfy the following condition

$$\begin{array}{c}\hfill a_m>{\sum }_{j=1+m}^M{a}_j{\gamma }_{\mathrm{th}j}.\end{array}$$

(17)

Let us denote ${\theta }^{*}=max({\theta }_1,{\theta }_2,\cdots ,{\theta }_m)$ [33], then the outage probability ${\mathrm{P}}_{out}^j$ of ${\mathrm{D}}_j$ can be reformulated as

$$\begin{array}{c}\hfill {\mathrm{P}}_{out}^j=1-Pr\left\{Z>\frac{{\theta }^{*}}{G^2(X-{\theta }^{*})},X>{\theta }^{*}\right\}.\end{array}$$

(18)

Using the conditional probability property [34] with respect to X, and applying the law of joint CDF, we have

$$\begin{array}{c}\hfill {\mathrm{P}}_{out}^j=1-{\int }_{{\theta }^{*}}^{\infty }\left[1-F_Z\left(\frac{{\theta }^{*}}{G^2(x-{\theta }^{*})}\right)\right]f_X\left(x\right)dx.\end{array}$$

(19)

To calculate the expression of the outage probability in (19), we first derive the CDF of Z and the PDF of X as follows.

When the nth relay is selected as the best relay, the PDF of order statistic with respect to $|h_{1,b}{|}^2$ in a set of N relays is obtained by using the binomial Newton expansion [35], that is,

$$\begin{array}{cc}\hfill f_{|h_{1,b}{|}^2}\left(x\right)& =N{\left[F_{|h_{1,i}{|}^2}\left(x\right)\right]}^{N-1}{f}_{|h_{1,i}{|}^2}\left(x\right)\hfill \\ & ={\displaystyle \sum _{n=1}^N}\left(\genfrac{}{}{0pt}{}{N}{n}\right){(-1)}^{n-1}\frac{n}{{\mathsf{\Omega }}_{1,n}}exp\left(-\frac{nx}{{\mathsf{\Omega }}_{1,n}}\right),\hfill \end{array}$$

(20)

where $\left(\genfrac{}{}{0pt}{}{N}{n}\right)=\frac{n!}{n!(N-n)!}$, N and n are non-negative integers, $f_{|h_{1,i}{|}^2}\left(x\right)=\frac{1}{{\mathsf{\Omega }}_{1,i}}exp\left(-\frac{x}{{\mathsf{\Omega }}_{1,i}}\right)$ and $F_{|h_{1,i}{|}^2}\left(x\right)=1-exp\left(-\frac{x}{{\mathsf{\Omega }}_{1,i}}\right)$ are respectively the CDF and PDF of $|h_{1,i}{|}^2$, which is the channel gain of each link from source node to relay. According to the probability theory, the PDFs of $|{\widehat{h}}_{1,b}{|}^2$ and $|h_{1,b}{|}^2$ which are respectively denoted by $f_{|{\widehat{h}}_{1,b}{|}^2}\left(\widehat{x}\right)$ and $f_{|h_{1,b}{|}^2}\left(x\right)$ can be calculated by using the joint PDF, i.e., $f_{|{\widehat{h}}_{1,b}{|}^2}\left(\widehat{x}\right)={\int }_0^{\infty }{f}_{|{\widehat{h}}_{1,b}{|}^2,{\left|h_{1,b}\right|}^2}\left(\widehat{x},x\right)dx$. Another way to calculate the joint PDF of $|{\widehat{h}}_{1,b}{|}^2$ is based on the properties of conditional probability, that is,

$$\begin{array}{c}\hfill f_{|{\widehat{h}}_{1,b}{|}^2}\left(\widehat{x}\right)={\int }_0^{\infty }{f}_{|{\widehat{h}}_{1,b}{|}^2\left|\right|h_{1,b}{|}^2}\left(\widehat{x}|x\right)f_{|h_{1,b}{|}^2}\left(x\right)dx,\end{array}$$

(21)

where

$$\begin{array}{c}\hfill f_{|{\widehat{h}}_{1,b}{|}^2\left|\right|h_{1,b}{|}^2}\left(\widehat{x}|x\right)=\frac{f_{|{\widehat{h}}_{1,i}{|}^2,{\left|h_{1,i}\right|}^2}\left(\widehat{x},x\right)}{f_{|h_{1,i}{|}^2}\left(x\right)}.\end{array}$$

(22)

Using the joint PDF which is given in ([36], Equation (9.389)), we can rewrite the numerator of (22) as

$$\begin{array}{c}\hfill f_{|{\widehat{h}}_{1,i}{|}^2,{\left|h_{1,i}\right|}^2}\left(\widehat{x},x\right)=\frac{exp\left(-\frac{(\widehat{x}+x)}{(1-{\rho }^2){\mathsf{\Omega }}_{1,n}}\right)}{(1-{\rho }^2){\mathsf{\Omega }}_{1,n}^2}{I}_0\left(\frac{2\rho \sqrt{\widehat{x}x}}{(1-{\rho }^2){\mathsf{\Omega }}_{1,n}}\right),\end{array}$$

(23)

where $I_0\left(x\right)$ is the modified zero order Bessel function of the first kind [22].

Without loss of generality, all correlation coefficients are assumed to have the same values, that is, $\rho ={\rho }_1={\rho }_2$. Substituting (20), (22), and (23) into (21), after using the equation ${\int }_0^{\infty }{e}^{-\alpha z}{I}_0\left(2\sqrt{\beta z}\right)dz=(1/\alpha )exp(\beta /\alpha )$ which is given in ([22], Equation (6.614.3)), and then perform some manipulations, we have the PDF of X in the case of imperfect CSI as

$$\begin{array}{cc}\hfill f_{|{\widehat{h}}_{1,b}{|}^2}\left(\widehat{x}\right)& ={\displaystyle \sum _{n=1}^N}\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{n{(-1)}^{n-1}}{{\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)}exp\left(-\frac{n\widehat{x}}{{\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)}\right),\hfill \end{array}$$

(24)

where $\mathsf{\Psi }(\rho ,n)=1+(n-1)(1-{\rho }^2)$.

From (24), the CDF of $|{\widehat{h}}_{1,b}{|}^2$ is given by

$$\begin{array}{cc}\hfill F_{|{\widehat{h}}_{1,b}{|}^2}\left(\widehat{x}\right)=1-& {\displaystyle \sum _{n=1}^N}\left(\genfrac{}{}{0pt}{}{N}{n}\right){(-1)}^{n-1}exp\left(-\frac{n\widehat{x}}{{\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)}\right).\hfill \end{array}$$

(25)

Based on the result of order statistics which is provided in ([34], Equation (7.14), p. 246), and after some similar calculations as above, the PDF of the ordered variable Z is expressed as

$$\begin{array}{cc}\hfill f_{|\widehat{\widehat{}}{g}_j{|}^2}\left(z\right)=& {\displaystyle \sum _{j=1}^M}\left(\genfrac{}{}{0pt}{}{M}{j}\right)\frac{{(-1)}^{j-1}j}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j)}exp\left(-\frac{jz}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j)}\right),\hfill \end{array}$$

(26)

where $\mathsf{\Psi }(\rho ,j)=1+(j-1)(1-{\rho }^2)$.

From (26), we can derive the CDF of $|{\widehat{g}}_i{|}^2$ as

$$\begin{array}{cc}\hfill F_{|{\widehat{g}}_i{|}^2}\left(\widehat{z}\right)=1-& {\displaystyle \sum _{j=1}^M}\left(\genfrac{}{}{0pt}{}{M}{j}\right){(-1)}^{j-1}exp\left(-\frac{j\widehat{z}}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j)}\right).\hfill \end{array}$$

(27)

Plugging (27) and (24) into (19), and after some manipulations, we obtain the expression of the outage probability as

$${\mathrm{P}}_{out}^j=1-{\displaystyle \sum _{j=1}^M}\left(\genfrac{}{}{0pt}{}{M}{j}\right){(-1)}^{j-1}{\displaystyle \sum _{n=1}^N}\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{n{(-1)}^{n-1}}{{\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)}{\int }_{{\theta }^{*}}^{\infty }exp\left(-\frac{j{\theta }^{*}}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j)G^2(x-{\theta }^{*})}-\frac{n\widehat{x}}{{\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)}\right)dx.$$

(28)

Let $u=x-{\theta }^{*}$, (28) becomes

$$\begin{array}{cc}{\mathrm{P}}_{out}^j=1-& {\displaystyle \sum _{j=1}^M}\left(\genfrac{}{}{0pt}{}{M}{j}\right){(-1)}^{j-1}{\displaystyle \sum _{n=1}^N}\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{n{(-1)}^{n-1}}{{\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)}exp\left(-\frac{n{\theta }^{*}}{{\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)}\right)\hfill \\ & \times {\displaystyle \underset{0}{\overset{\infty }{\int }}}exp\left(-\frac{j{\theta }^{*}}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j)G^{2}u}-\frac{nu}{{\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)}\right)du.\hfill \end{array}$$

(29)

Using ([22], Equation (3.324)), we can rewrite the exact closed-form expression of the outage probability as in (30), where $K_1(.)$ denotes the modified first order Bessel function of the second kind.

$$\begin{array}{cc}{\mathrm{P}}_{out}^j=1-& {\displaystyle \sum _{j=1}^M}\left(\genfrac{}{}{0pt}{}{M}{j}\right){(-1)}^{j-1}{\displaystyle \sum _{n=1}^N}\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{{(-1)}^{n-1}}{{\mathsf{\Omega }}_{1,n}}exp\left(-\frac{n{\theta }^{*}}{{\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)}\right)\hfill \\ & \times \sqrt{\frac{4nj{\theta }^{*}}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j){\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)G^2}}{K}_1\left(\sqrt{\frac{4nj{\theta }^{*}}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j){\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)G^2}}\right).\hfill \end{array}$$

(30)

From the expression of the outage probability which is given in (30), we can see that when the outdated CSI happens, the outage performance is a function of $\rho$.

4. Ergodic Capacity Analysis

In this section, we analyze the ergodic capacity of the proposed NOMA relaying system in comparison with that of the OMA relaying system. Due to the fact that the hardware complexity and performance degradation of the NOMA system is directly proportional to the number of users, we also set the number of users be equal to three for both NOMA and OMA systems as used in [37]. For the OMA system, we consider orthogonal frequency division multiple access (OFDMA). According to the Shannon theory, the instantaneous rate of ${\mathrm{D}}_m$ is given by

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{NOMA}}^{m\mathrm{th}}=\frac{1-\alpha }{2}{log}_2\left(1+{\gamma }_m^{\mathrm{D}}\right).\end{array}$$

(31)

From (10), when the transmission power is high, we can approximate the required SINR at ${\mathrm{D}}_m$ as

$$\begin{array}{c}\hfill {\gamma }_m^{\mathrm{D}}\approx \frac{G^2{a}_m{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2}{G^2{\sum }_{j=m+1}^M{a}_j{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2+{\sigma }_{{\mathrm{D}}_m}^2}.\end{array}$$

(32)

Substituting (32) into (31), we have

$$\begin{array}{cc}{\mathcal{R}}_{\mathrm{NOMA}}^{m\mathrm{th}}& \approx \frac{1-\alpha }{2}{log}_2\left(1+\frac{G^2{a}_m{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2}{G^2{\sum }_{j=m+1}^M{a}_j{P}_{\mathrm{S}}|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2+{\sigma }_{{\mathrm{D}}_m}^2}\right)\hfill \\ & =\frac{1-\alpha }{2}{log}_2\left(\frac{G^{2}P|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2+1}{G^2{\sum }_{j=m+1}^M{a}_{j}P|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2+1}\right),\hfill \end{array}$$

(33)

where $P=\frac{P_{\mathrm{S}}}{{\sigma }_{{\mathrm{D}}_m}^2}$.

Based on the properties of the logarithmic function, we can rewrite (33) as

$$\begin{array}{cc}{\mathcal{R}}_{\mathrm{NOMA}}^{m\mathrm{th}}=& \underset{{\mathcal{I}}_1}{\underbrace{\frac{1-\alpha }{2}\mathbb{E}\left\{{log}_2\left(1+G^{2}P|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2\right)\right\}}}\hfill \\ & -\underset{{\mathcal{I}}_2}{\underbrace{\frac{1-\alpha }{2}\mathbb{E}\left\{{log}_2\left(1+G^2{\sum }_{m=1}^{M-1}{a}_{m}P|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2\right)\right\}}},\hfill \end{array}$$

(34)

then solve its components by using the partial integration, i.e.,

$$\begin{array}{cc}\hfill {\mathcal{I}}_u=& {\left\{{log}_2(1+{\mathsf{\Gamma }}_u)\left[F_{{\mathsf{\Gamma }}_u}\left(x_u\right)-1\right]\right\}}_0^{\infty }-\frac{1}{2\ln 2}{\int }_0^{\infty }\frac{1}{1+x_u}\left[F_{{\mathsf{\Gamma }}_u}\left(x_u\right)-1\right]d{x}_u\hfill \\ \hfill =& \frac{1}{2\ln 2}{\int }_0^{\infty }\frac{1}{1+x_u}\left[1-F_{{\mathsf{\Gamma }}_u}\left(x_u\right)\right]d{x}_u,\hfill \end{array}$$

(35)

where $F_{{\mathsf{\Gamma }}_u}\left(x_u\right)$ is the CDF of random variable ${\mathsf{\Gamma }}_u$ with $u\in \{1,2\}$, ${\mathsf{\Gamma }}_1=G^{2}P|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2$, and ${\mathsf{\Gamma }}_2=G^2{\sum }_{m=1}^{M-1}{a}_{m}P|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_m\right|}^2$.

Using the condition probability, we have CDF of ${\mathsf{\Gamma }}_1$ as

$$\begin{array}{cc}{F}_{{\mathsf{\Gamma }}_1}\left(x_1\right)& =Pr(G^{2}P|{\widehat{h}}_{1,b}{|}^2|{\widehat{g}}_m{|}^2\le x_1)\hfill \\ & ={\int }_0^{\infty }Pr\left(|{\widehat{g}}_m{|}^2\le \frac{x_1}{G^{2}P{\left|{\widehat{h}}_{1,b}\right|}^2}\right)f_{|{\widehat{h}}_{1,b}{|}^2}d{\left|{\widehat{h}}_{1,b}\right|}^2.\hfill \end{array}$$

(36)

From (24) and (27) we can calculate $F_{{\mathsf{\Gamma }}_1}\left(x_1\right)$ as

$$\begin{array}{cc}{F}_{{\mathsf{\Gamma }}_1}\left(x_1\right)=1-& {\displaystyle \sum _{j=1}^M}\left(\genfrac{}{}{0pt}{}{M}{j}\right){(-1)}^{j-1}{\displaystyle \sum _{n=1}^N}\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{{(-1)}^{n-1}}{{\mathsf{\Omega }}_{1,n}}\hfill \\ & \times \sqrt{\frac{4nj{x}_1}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j){\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)P{G}^2}}{K}_1\left(\sqrt{\frac{4nj{x}_1}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j){\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)P{G}^2}}\right).\hfill \end{array}$$

(37)

Similarly, for $F_{{\mathsf{\Gamma }}_2}\left(x_2\right)$, we have

$$\begin{array}{cc}{F}_{{\mathsf{\Gamma }}_2}\left(x_2\right)=1-& {\displaystyle \sum _{j=1}^M}\left(\genfrac{}{}{0pt}{}{M}{j}\right){(-1)}^{j-1}{\displaystyle \sum _{n=1}^N}\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{{(-1)}^{n-1}}{{\mathsf{\Omega }}_{1,n}}\hfill \\ & \times \sqrt{\frac{4nj{x}_2}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j){\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)bP{G}^2}}{K}_1\left(\sqrt{\frac{4nj{x}_2}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j){\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)bP{G}^2}}\right).\hfill \end{array}$$

(38)

where $b={\sum }_{m=1}^{M-1}{a}_m$.

Replacing (38) into (35), we obtain ${\mathcal{I}}_1$ as

$$\begin{array}{cc}{\mathcal{I}}_1=& \frac{1}{2\ln 2}{\displaystyle \sum _{j=1}^M}\left(\genfrac{}{}{0pt}{}{M}{j}\right){(-1)}^{j-1}{\displaystyle \sum _{n=1}^N}\left(\genfrac{}{}{0pt}{}{N}{n}\right)\frac{{(-1)}^{n-1}}{{\mathsf{\Omega }}_{1,n}}\hfill \\ & \times {\int }_0^{\infty }\frac{1}{1+x_1}\sqrt{\mathcal{A}(n,j)x_1}{K}_1\left(\sqrt{\mathcal{A}(n,j)x_1}\right)d{x}_1,\hfill \end{array}$$

(39)

where $\mathcal{A}=\frac{4}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j){\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)P{G}^2}$.

Based on ([22], Equation (9.343)), we can rewrite (39) as

$$\begin{array}{c}\hfill {\mathcal{I}}_1=\frac{1-\alpha }{2\sqrt{2}ln2}{\displaystyle \underset{0}{\overset{\infty }{\int }}}\frac{1}{1+x_1}{G}_{02}^{20}\left(\left.\frac{x_1}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j){\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)P{G}^2}\right|\begin{array}{c}\\ \frac{3}{4},-\frac{1}{4}\end{array}\right)d{x}_1.\end{array}$$

(40)

Then, using ([22], Equation (7.811.5)) and after some manipulations, we have

$$\begin{array}{c}\hfill {\mathcal{I}}_1=\frac{1-\alpha }{2\sqrt{2}ln2}{G}_{13}^{31}\left({\left.\frac{1}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j){\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)P{G}^2}\right|}_{0,\frac{3}{4},-\frac{1}{4}}^0\right),\end{array}$$

(41)

where $G_{pq}^{mn}\left(x{|}_{b_s}^{a_r}\right)$ is the Meijer’s G-Function ([22], Equation (9.3)).

Plugging (38) into (35), and doing similar manipulations which were used to derive ${\mathcal{I}}_1$, we obtain

$$\begin{array}{c}\hfill {\mathcal{I}}_2=\frac{1-\alpha }{2\sqrt{2}ln2}{G}_{13}^{31}\left({\left.\frac{1}{{\mathsf{\Omega }}_z\mathsf{\Psi }(\rho ,j){\mathsf{\Omega }}_{1,n}\mathsf{\Psi }(\rho ,n)bP{G}^2}\right|}_{0,\frac{3}{4},-\frac{1}{4}}^0\right),\end{array}$$

(42)

To compare the ergodic capacities of the NOMA and OMA systems, we let $\beta$ be the bandwidth which is assigned for ${\mathrm{D}}_1$ and $(1-\beta )/2$ be the remaining bandwidth which is assigned for ${\mathrm{D}}_2$ and ${\mathrm{D}}_3$, where $(0<\beta <1)$ and the whole bandwidth is 1Hz. From ([38], Equation (7.4)), we can extend the achievable end-to-end ergodic capacity of the OFDMA system with three users as

$$\begin{array}{cc}\hfill {\mathcal{R}}_{\mathrm{OMA}}=& \frac{1-\alpha }{2}\beta {log}_2\left(1+{\gamma }_{{\mathrm{SRD}}_1}\right)\hfill \\ & +\frac{(1-\alpha )(1-\beta )}{4}{log}_2\left(1+{\gamma }_{{\mathrm{SRD}}_2}\right)\hfill \\ & +\frac{(1-\alpha )(1-\beta )}{4}{log}_2\left(1+{\gamma }_{{\mathrm{SRD}}_3}\right),\hfill \end{array}$$

(43)

where ${\gamma }_{{\mathrm{SRD}}_m},m\in \{1,2,3\}$ denotes the instantaneous SINR of each user, which is computed as

$$\begin{array}{cc}\hfill {\gamma }_{{\mathrm{SRD}}_1}& =\frac{G^2{P}_{\mathrm{S}}^{\mathrm{OMA}}|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_1\right|}^2}{\beta \left(G^2\right|{\widehat{g}}_1{|}^2{\sigma }_{\mathrm{R}}^2+{\sigma }_{{\mathrm{D}}_1}^2)},\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill {\gamma }_{{\mathrm{SRD}}_2}& =\frac{2G^2{P}_{\mathrm{S}}^{\mathrm{OMA}}|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_2\right|}^2}{(1-\beta )\left(G^2\right|{\widehat{g}}_2{|}^2{\sigma }_{\mathrm{R}}^2+{\sigma }_{{\mathrm{D}}_2}^2)},\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\gamma }_{{\mathrm{SRD}}_3}& =\frac{2G^2{P}_{\mathrm{S}}^{\mathrm{OMA}}|{\widehat{h}}_{1,b}{|}^2{\left|{\widehat{g}}_3\right|}^2}{(1-\beta )\left(G^2\right|{\widehat{g}}_3{|}^2{\sigma }_{\mathrm{R}}^2+{\sigma }_{{\mathrm{D}}_3}^2)},\hfill \end{array}$$

(46)

where $P_{\mathrm{S}}^{\mathrm{OMA}}=P_{\mathrm{S}}/3$ is the equal power allocated for the signal transmission from $\mathrm{S}$ to each user ${\mathrm{D}}_m$ ([38], p. 146). The factor $\frac{1-\alpha }{2}$ appears in (31) and (43) because source node transmits its signals to all users in two time slots of the transmission period T.

5. Numerical Results

In this section, we provide the numerical results to evaluate the system performance in terms of the outage probability (OP) and ergodic capacity of the proposed EH-NOMA relaying system with three users. We also determine the optimal time switching ratio to minimize the OP and compare the ergodic capacities of the proposed EH-NOMA relaying system with EH-OMA relaying system. Regarding to the evaluating method, we use the common approach in this field, that is, to drive a closed-form mathematical expression to model the system performance and then compare the analysis results with Monte-Carlo simulation results to validate the derived mathematical expressions. Unlike previous works, which only considered EH-NOMA systems with two users and under perfect CSI, our paper focuses on the theoretical analysis of an EH-NOMA system with more than two users, taking into account the effects of AF relaying protocol and the feedback delay of wireless channels on the system performance. Since, there are not many similar parameters, it may be an unfair comparison between our proposed EH-NOMA relaying system with previous NOMA relaying systems. Therefore, we use the same system model of the proposed EH-NOMA relaying sytem but replace the NOMA with OMA to demonstrate the benefits of utilizing the NOMA technique in the proposed EH relaying system. Unless otherwise stated, the parameter settings of EH-NOMA and EH-OMA relaying systems are summarized in

Table 2. Parameter settings of EH-NOMA and EH-OMA relaying systems.

Description	EH-NOMA	EH-OMA
Allocated transmission power	$P_1=0.7P_S,P_2=0.2P_S$, $P_3=0.1P_S$	$P_i=P_S/3$
Bandwidth	$\beta$ for ${\mathrm{D}}_1$, $(1-\beta )/2$ for ${\mathrm{D}}_2$ and ${\mathrm{D}}_3$	B = 1 Hz for all users
Target data rate	$r=0.5$ bpcu
Time switching ratio	$\alpha =0.3$
Average channel gain	${\mathsf{\Omega }}_{1,n}=1$, ${\mathsf{\Omega }}_{{\mathrm{R}}_n{\mathrm{D}}_1}=2$, ${\mathsf{\Omega }}_{{\mathrm{R}}_n{\mathrm{D}}_2}=3$, ${\mathsf{\Omega }}_{{\mathrm{R}}_n{\mathrm{D}}_3}=6$
Energy conversion efficiency	$\eta =0.85$

. It is noticed that the average $\mathrm{SNR}$ is defined as the ratio of the transmission power of source S to the variance of AWGN, that is, $\mathrm{SNR}=P_S/{\sigma }^2$, ranging from 0 dB to 40 dB.

Figure 2 shows the outage probability of each user versus the average SINR in dB. The outage probability of the EH-NOMA relaying system is also compared with that of EH-OMA relaying system. Firstly, we can see that the OP of ${\mathrm{D}}_3$ is lowest among all users while the OP of ${\mathrm{D}}_1$ is highest. The reason is that the channel gain from R to ${\mathrm{D}}_3$ is highest (the decay of the magnitude power signal is proportional to the squared distance in multipath fading) because ${\mathrm{D}}_3$ is the closest user to R while ${\mathrm{D}}_1$ is the farthest one. Another important observation is that the OPs of ${\mathrm{D}}_2$ and ${\mathrm{D}}_3$ in the EH-NOMA relaying system are better than those of ${\mathrm{D}}_2$ and ${\mathrm{D}}_3$ in the EH-OMA relaying system, while the OP of ${\mathrm{D}}_1$ in the EH-OMA relaying system is better than in the NOMA relaying system. However, the gap is insignificant because the number of time slots for the transmission in the EH-OMA relaying system is higher than in the EH-NOMA relaying system, thus the probability that outage evens happen in the EH-OMA relaying system is also higher than in the EH-NOMA relaying system. On the other hand, the outage threshold of the OMA user is ${\gamma }_{\mathrm{th}}^{\mathrm{OMA}}=2^{\frac{2r}{v(1-\alpha )}}-1$, where $v\in \{\beta ,(1-\beta )/2\}$. In contrast, the outage threshold of the NOMA user is ${\gamma }_{\mathrm{th}}=2^{\frac{2r}{(1-\alpha )}}-1$. Then, obviously the outage threshold of the OMA user is obviously higher than that of the NOMA user. However, the OP not only depends on the outage threshold but also on the received SINR at user. In addition, we also see that in the low SINR regime (less than 15 dB), the OPs of all OMA users always outperform those of NOMA users. However, in the high SINR regime (larger 15 dB) only the OP of ${\mathrm{D}}_1$ in the EH-OMA relaying system is better than that in the EH-NOMA relaying system. We can also see in Figure 2 that the diversity gain of all users is equal to one.

Figure 3 plots the OP of ${\mathrm{D}}_1$ in the EH-NOMA relaying system versus the average SINR in dB for different channel correlation coefficients $\rho$. Firstly, we see that higher $\rho$ reduces the OP, but the reduction is not remarkable for small $\rho$. The improvement in OP is only significant when $\rho$ is near to 1. We should remind that $\rho$ indicates the correlation degree between the transmission channel and the feedback channel in time coherent at the transmitter. The analysis results are in excellent agreement with the simulation ones, confirming the correctness of our mathematical analysis.

Figure 4 illustrates the OP of ${\mathrm{D}}_1$ in the NOMA system versus the average SINR in dB for different numbers of relays N. From Figure 4, we see that when the number of relays increases, the outage performance of the system is improved. It is because increasing the number of relays will provide more opportunity for selecting the connection links from source node to relay, which not only makes the achievable decoding performance better but also increases the amount of harvested energy. In addition, the diversity gains is not significantly improved with N because the diversity order of PRS scheme is always equal to one.

Figure 5 presents the OP of ${\mathrm{D}}_1$ in the EH-NOMA relaying system versus time switching ratio $\alpha$ for different numbers of relays. The values of $\alpha$ range from 0 to 0.7 while SINR remains at 15 dB. Firstly, we see that there exists an optimal value of $\alpha$ which minimizes the OP. Moreover, the minimum value of OP depends on the number of relays N, i.e., as N is higher the minimal OP becomes smaller. The reason is that when N increases, the SINR of the first hop will be better because the PRS method is used. Another important observation is that the optimal value of $\alpha$ which minimizes the OP is approximately 0.2 regardless of the number of relays.

Figure 6 demonstrates the OP of ${\mathrm{D}}_1$ in the EH-NOMA relaying system versus the correlation coefficient $\rho$ for different average SINR. We can see that the OP reduces as $\rho$ increases. In the worst case $\rho =0$, the instantaneous CSI at the transmission time does not correlate with the instantaneous CSI at the relay-selection time or at the power-allocation time. In contrast, in the best case $\rho =1$, the instantaneous CSI at the transmission time closely correlates with the instantaneous CSI at the relay-selection time or at the power-allocation time. The improvement in the CSI leads to better power allocation and signal processing of the system. Figure 6 also shows that when $\rho <0.8$, the enhancement of OP is not significant and the system performance is only improved when the correlation coefficient $\rho$ is close to 1.

Figure 7 depicts the ergodic capacity of each user in EH-NOMA relaying system versus the average SINR in dB. As observed from Figure 7, the ergodic capacity of ${\mathrm{D}}_3$ outperforms the ergodic capacities of ${\mathrm{D}}_1$ and ${\mathrm{D}}_2$. Moreover, the ergodic capacities of ${\mathrm{D}}_1$ and ${\mathrm{D}}_2$ increase slightly in the low SINR region and is saturated in the high SINR region. In contrast, the ergodic capacity of ${\mathrm{D}}_3$ increases exponentially with respect to the SINR. This reason is that ${\mathrm{D}}_1$ does not use SIC but only detects the signal of itself. Meanwhile, ${\mathrm{D}}_2$ must use the first-order SIC first and then ${\mathrm{D}}_3$ uses the second-order SIC. Thus, the impact of interference on ${\mathrm{D}}_1$ is higher than ${\mathrm{D}}_2$ and ${\mathrm{D}}_3$. However, there exists the trade-off between the complexity and the achievable ergodic capacity of the system. We also see a good match between the analysis results and the simulation results, especially in the high SINR regime. On the other hand, the ergodic capacity in the case of perfect SIC is compared with that in the case of imperfect SIC. We can see that the ergodic capacity in case imperfect SIC is lower. Moreover, the gap between them increases with the SINR. It is because when SINR increases, the interference caused by imperfect SIC also increases. Therefore, the SINR as well as the ergodic capacity become slowly higher. Another feature is that the ergodic capacity of ${\mathrm{D}}_1$ remains the same in both cases because ${\mathrm{D}}_1$ does not use SIC when decoding the signals.

Figure 8 provides the simulation results of the ergodic rate of NOMA and OMA systems versus the average SINR in dB. From Figure 8, we see that the ergodic rate of NONA system is always higher than the OMA system as the number of relays increases. It is because the NOMA system uses the whole bandwidth for each user while the OMA system uses individual bandwidth for each user, resulting in higher spectrum usage efficiency. Another important observation is that when N gets higher, the difference gap of the ergodic capacities of these two systems does not increase linearly. Thus, we do not need to use a large number of relays for partial relay selection scheme because it may increase the complexity of the system but not significantly enhance its performance.

6. Conclusions

In this paper, we propose a downlink NOMA relaying system with the best RF energy harvesting relay and investigate the impact of CSI imperfection on the performance of the proposed NOMA relaying system over Rayleigh fading channel. Specifically, we provide detailed derivations of the exact closed-form expression of OP and the approximate expression of the ergodic capacity of the proposed NOMA relaying system. Based on the expression of the OP, the optimal energy harvesting duration which minimizes the OP in the condition of imperfect CSI can be determined. The results show that imperfect CSI significantly reduces the system performance. In addition, we show that the spectrum efficiency of our proposed NOMA relaying system outperforms that of the OMA relaying system in the same parameter settings. All analysis results are in excellent agreement with the simulation results, confirming the correctness of the mathematical analysis. The proposed EH-NOMA relaying system can support the communication for multiple users through the best relay without relying on the external power supply. Thus it can be applied in surveillance sensor networks for disaster detection or in Internet of Things (IoT) where installing fixed power lines or frequent battery replacement for a large number of nodes may be not convenient. Using the results in this paper, we can choose an appropriate time switching ratio to balance between the energy harvesting and signal processing so that the outage probability of the proposed EH-NOMA relaying system system can be reduced upto 76.32%. Moreover, compared with the EH-OMA relaying system, the OP of the proposed EH-NOMA relaying system is 9.41% lower and the ergodic capacity is 17.64% higher at the average SNR = 40 dB.