Outage Analysis of the Power Splitting Based Underlay Cooperative Cognitive Radio Networks

Abstract

In the present paper, we investigate the performance of the simultaneous wireless information and power transfer (SWIPT) based cooperative cognitive radio networks (CCRNs). In particular, the outage probability is derived in the closed-form expressions under the opportunistic partial relay selection. Different from the conventional CRNs in which the transmit power of the secondary transmitters count merely on the aggregate interference measured on the primary networks, the transmit power of the SWIPT-enabled transmitters is also constrained by the harvested energy. As a result, the mathematical framework involves more correlated random variables and, thus, is of higher complexity. Monte Carlo simulations are given to corroborate the accuracy of the mathematical analysis and to shed light on the behavior of the OP with respect to several important parameters, e.g., the transmit power and the number of relays. Our findings illustrate that increasing the transmit power and/or the number of relays is beneficial for the outage probability.

1. Introduction

Cognitive radio networks (CRNs) are considered one of the most effective solutions to overcome the scarcity of the frequency spectrum [1,2]. The principal idea of CRNs is to permit unlicensed users to concurrently operate with licensed users while strictly guaranteeing the quality-of-service (QoS) of primary users. To realize such networks two popular protocols are proposed in the literature, namely, the overlay and underlay protocols [3,4]. The former allows secondary devices to opportunistically occupy the temporarily unused spectrum in terms of space, time, and frequency. Regarding the underlay protocol, on the other hand, secondary users are always granted permission to access the licensed spectrum provided that the aggregate interference created by secondary networks measured on the primary devices is below the predefined threshold. Compared to the overlay protocol, the underlay protocol is preferable due to its high availability and, thus, is suitable for urgent services, video conferences, video gaming, and so forth. Nevertheless, the cons of this protocol are that it does not support a long transmission and/or high-quality services owing to the low transmit power to avoid exceeding the interference threshold of the primary networks. As a consequence, in order to employ the underlay CRNs in practice, the combination with other techniques is necessary. Fortunately, relaying or cooperative communications is a good complement to the underlay CRNs [5,6,7]. More precisely, cooperative communications networks are wireless networks where one or several relays are deployed among end-users to help them exchange information. With the help of the relay, the transmission distance is dramatically decreased, improving the system performance. Additionally, relaying technology has also been proven to be an effective way to extend the coverage area. The spectral efficiency and QoS issues can be addressed in a straightforward manner by employing the cooperative cognitive radio networks (CCRNs) [8,9]. Nonetheless, there still exists an urgent issue, which is to improve the energy efficiency of wireless networks. The problem escalates seriously in either the Internet of Things (IoTs), low power wide area networks (LPWAN) [10], or 5G and beyond networks owing to the exponential growth of wireless-connected devices accompanied by their power-hungry applications. In fact, improving energy efficiency is one of the highest priorities in wireless communications. Fortunately, simultaneous wireless information and power transfer (SWIPT) has recently been proposed to deal with this issue [11,12,13]. In particular, SWIPT transmits both information and energy on the same carrier frequency, thus boosting both the spectral and energy efficiency. Consequently, in the present paper, we investigate the performance of the SWIPT-based cooperative cognitive radio networks (CCRNs) with the help of multiple relays. Before highlighting our novel contributions, the state-of-the-art is given as follows.

The performance of CCRNs networks was studied widely in [14,15,16,17,18,19,20]. Specifically, the outage probability (OP) of the cognitive radio non-orthogonal multiple access (NOMA) networks was derived in [14]. The authors in [15], on the other hand, investigated the secrecy performance of the underlay cooperative multihop CRNs. In particular, the secrecy outage probability (SOP) was derived in the approximated closed-form expression. Lee et al. in [16] addressed the minimization of the number of feedback bits required in order to satisfy the QoS in both primary and secondary networks. The optimal power allocation of non-orthogonal amplify–and–forward (AF) relaying of underlay CRNs was provided in [17] to maximize the system throughput. In [18], the capacity of voice over IP (VoIP) in CRNs was analyzed and maximized by modelling the VoIP traffic and channel coefficients as a Markov-modulated Poisson process. Liu et al. in [19] investigated the spectrum sensing problem in the full-duplex cooperative spectrum sensing CRNs. It was noted that spectrum sensing is important in overlay CRNs to guarantee the opportunistic access of secondary devices while not interrupting the primary device’s transmission.

Meanwhile, the performance of SWIPT-enabled networks was investigated in [21,22,23,24,25,26,27,28,29]. The energy-efficient optimization of SWIPT-assisted relaying networks was addressed in [21]. Tan et al. in [22] derived the OP and ergodic capacity of the power splitting (PS) based relaying networks under the asymmetric channel, i.e., the Nakagami-m and Rayleigh channels. The error probability and outage probability of SWIPT-based NOMA networks were derived in [23]. In particular, the pairwise error probability was computed in closed-form expression in [23]. The asymptotic framework under a high SNRs regime was provided as well. The work in [24], differently, derived the symbol error rate (SER) of SWIPT-enabled relaying networks, where the noncoherent modulation was employed in place of the conventional phase-shift keying (PSK) and/or quadrature amplitude modulation (QAM). The authors in [25] also addressed the noncoherent modulation. Specifically, this work first derived the moments and moment generating function (MGF) of the end-to-end (e2e) signal-to-noise ratios (SNRs). Based on the MGF, they then computed the outage probability, the amount of fading, and the system throughput. The authors in [26,27] dealt with the physical layer security (PLS) issue of the simultaneous wireless information and power transfer assisted relaying networks. Furthermore, the performance of SWIPT-based cellular networks with and without utilizing millimeter wave (mmWave) was provided in [28,29].

Despite the extensive study of either the CCRNs or the SWIPT-aided networks, the performance of SWIPT-enabled underlay CCRNs is still in the infancy stage. In particular, there only a few works addressing this combination [30,31,32,33]. Specifically, Prathima et al. in [30] studied the performance of the primary networks with the help of the secondary users that also act as the relay of primary users. This work, however, concentrated on the performance of primary networks. Moreover, it solely considered two relay nodes instead of the general scenario. The work in [31], diversely, addressed the secrecy performance of the SWIPT-assisted cognitive relaying networks. The power allocation and transceiver design were investigated in [32].

In this paper, different from the abovementioned works, we focus on the performance of the secondary networks as well as the reliability of SWIPT-based underlay cooperative cognitive radio networks with regard to transceiver design, power allocation, and physical layer security. In particular, the principal novel contributions are summarized as follows:

• We consider a single-input single-output (SISO) underlay cooperative cognitive radio network with the assistance of multiple relays. Additionally, the transmit power of the relay nodes relies only on the harvested energy from the transmitter S. The partial relay selection is adopted to both enhance the system performance and reduce the complexity compared to the fully relay selection.
• Different from the conventional underlay CRNs where the transmit power of the secondary transmitter considers merely the interference power, the transmit power of the considered networks is constrained by both the interference power and the harvested energy. The mathematical framework, thus, is of higher complexity owing to dealing with more correlated random variables. Nonetheless, we are still able to derive the outage probability in the closed-form expressions.
• Simulation results are presented to corroborate the exactness of our analysis and to identify the behavior of OP with respect to several important parameters, namely, the transmit power, the number of relays, the power splitting ratio, and so on. Our findings show that both increasing the transmit power and number of relays are beneficial to the OP. Additionally, an optimal value of the power splitting ratio exists that minimizes the OP.

The remainder of this paper is organized as follows. The system model is given in Section 2. The derivation of the OP is provided in Section 3. Numerical results are shown in Section 4. Section 5 concludes the paper.

2. System Model

Let us consider SWIPT-based underlay cognitive radio networks as shown in Figure 1. In particular, the secondary networks comprise one source node denoted by S, one destination denoted by D, and M relay nodes denoted by $R_i,i\in \left\{1,\dots ,M\right\}$, while the primary networks are represented by a primary receiver denoted by P. Here, P measures the aggregate interference created by the secondary networks on the primary networks.

2.1. Channel Modeling

Considering a generic transmission from node X to node Y, the channel coefficients denoted by $h_{XY}$, $X\in \left\{\mathrm{S},{\mathrm{R}}_i\right\}$, $Y\in \left\{{\mathrm{R}}_i,\mathrm{D}\right\}$ are followed by a Rayleigh distribution. As a result, the channel gain denoted by ${\gamma }_{XY}={\left|h_{XY}\right|}^2$ is followed by an exponential distribution with parameter ${\lambda }_{XY}$ whose cumulative distribution function (CDF) and probability density function (PDF) are given as follows [6]:

$$\begin{array}{c}\hfill F_X\left(x\right)=1-exp(-{\lambda }_{XY}x),\end{array}$$

(1)

$$\begin{array}{c}\hfill f_X\left(x\right)=\frac{\partial F_X\left(x\right)}{\partial x}={\lambda }_{XY}exp(-{\lambda }_{XY}x).\end{array}$$

(2)

Here ${\lambda }_{XY}$ is also the large-scale path loss from X to Y and is formulated as follows:

$$\begin{array}{c}\hfill {\lambda }_{XY}={\left(d_{XY}\right)}^{\beta },\end{array}$$

(3)

where $d_{XY}$ is the Euclidean distance between node X and Y and $\beta \in \left(2,\dots ,6\right)$ is the path loss exponent. Additionally, the block fading is taken into consideration in this work, hence the channel coefficients remain constants for the whole transmission T and change independently between each transmission.

2.2. PS-Based Relaying Networks

In this work, we adopt the power-splitting (PS) protocol at the relay node. To be more precise, the received power at R is divided into two separate parts according to the power-splitting ratio $\rho ,0<\rho <1$, i.e., one is put into the energy harvester and another goes to the information decoder. $\rho$ takes into account all loss introduced by the energy harvesting receiver, e.g., noise introduced by the received antenna, loss due to the converting RF-to-DC circuit, and so on [34,35]. Additionally, to realize the power-splitting protocol, each SWIPT-enabled receiver needs to be equipped with a power splitter to split the received power into two parts. The first part is sent to the conventional information decoding circuit, and the remaining part is sent to the energy harvesting circuit [35,36].

2.3. Opportunistic Partial Relaying (OPR) Protocol

In this paper, the opportunistic partial relaying (OPR) protocol is adopted. In particular, only the relay n denoted by R${}_n$, which has the highest channel gain from S to all relay nodes, is selected to help exchange information between S and D. Other relay nodes, as a result, keep silent in order to save energy consumption and avoid creating co-channel interference.

$$\begin{array}{c}\hfill R_n:{\gamma }_{\mathrm{S}{\mathrm{R}}_n}=\underset{m=1,2,\dots ,M}{\underbrace{max}}\left\{{\gamma }_{\mathrm{S}{\mathrm{R}}_m}\right\}\end{array}$$

(4)

Compared with the scenario where all relays participate in the transmission, our adopted protocol is simpler since it does not require perfect channel state information (CSI) of all nodes of the networks at the destination and perfect synchronization among relays [37,38,39]. To be more precise, the adopted OPR protocol can be employed as follows. Each relay is equipped with a timer, and the value of the timer is set inversely with the channel gain from S to relay. Thus, the best relay is the one having the smallest timer. When the timer ends, the best relay forwards the source’s signal to the destination. Other relays sense the availability of the medium and keep silent once the medium is occupied.

2.4. Information Transmission

The whole transmission takes place in two phases. In the first phase, source S broadcasts its signals to all relay nodes. Here, we assume that the direct channel between S and D does not exist due to the long transmission distance and deep fades; thus, destination D does not receive the broadcast signal from S. Although all relays are received signals sent by S, only relay R${}_n$ is selected to assist the transmission from S to D. The criteria for selecting R${}_n$ is given in Section 2.3. At relay R${}_n$, parts of the incoming signals are sent to the information decoder to decode the information sent by S and are given as

$$\begin{array}{c}\hfill y_{{\mathrm{R}}_n}=\sqrt{1-\rho }\sqrt{P_{\mathrm{S}}}{h}_{\mathrm{S}{\mathrm{R}}_n}{x}_{\mathrm{S}}+n_{{\mathrm{R}}_n},\end{array}$$

(5)

where $n_{{\mathrm{R}}_n}$ is the additive white Gaussian noise (AWGN) at relay R${}_n$, which follows a complex Gaussian distribution with zero mean and ${\mathrm{N}}_0$ variance, $n_{{\mathrm{R}}_n}\sim \mathcal{CN}\left(0,{\mathrm{N}}_0\right)$; $x_{\mathrm{S}}$ is the transmitted signal of S and $\mathrm{E}\left\{{\left|x_{\mathrm{S}}\right|}^2\right\}=1$; $\mathrm{E}\left\{•\right\}$ is the expectation operator; and $P_{\mathrm{S}}$ is the transmit power of S and is defined in the sequel. The remaining part of the incoming signals from S is put into the energy harvested receiver. The amount of harvested energy denoted by $E_{{\mathrm{R}}_n}$ are then formulated as

$$\begin{array}{c}\hfill E_{{\mathrm{R}}_n}=\eta \rho (T/2)P_{\mathrm{S}}{\left|h_{\mathrm{S}{\mathrm{R}}_n}\right|}^2,\end{array}$$

(6)

where $\eta$ is the energy conversion coefficient [34,40]; the factor $T/2$ implies that the energy harvesting only takes place in half of the whole transmission procedure. At the end of the first phase, relay R${}_n$ decodes the information sent by S and forwards the re-encoded version to the destination D in the second phase. The whole transmission procedure is shown in Figure 2. The received signals at D is then formulated as

$$\begin{array}{c}\hfill y_{\mathrm{D}}=\sqrt{P_{{\mathrm{R}}_n}}{h}_{{\mathrm{R}}_n\mathrm{D}}{x}_{{\mathrm{R}}_n}+n_{\mathrm{D}},\end{array}$$

(7)

where $n_{\mathrm{D}}$ is the AWGN noise at D with zero mean and ${\mathrm{N}}_0$ variance; $x_{{\mathrm{R}}_n}$ is the transmitted signals of the relay R${}_n$ with $\mathrm{E}\left\{{\left|x_{{\mathrm{R}}_n}\right|}^2\right\}=1$, and $P_{{\mathrm{R}}_n}$ is the transmit power of R${}_n$ and is defined in Section 2.5. It is noted that the received signal at $R_n$ and D in Equations (5) and (7) is a function of the large-scale path loss via the channel coefficient $h_{\mathrm{S}{\mathrm{R}}_n}$ and $h_{{\mathrm{R}}_n\mathrm{D}}$, respectively. The signal-to-noise ratios at R${}_n$ and D are then formulated as follows:

$$\begin{array}{cc}\hfill {\gamma }_{{\mathrm{R}}_n}=& \frac{(1-\rho ){\gamma }_{\mathrm{S}{\mathrm{R}}_n}{P}_{\mathrm{S}}}{N_0},\hfill \\ \hfill {\gamma }_{\mathrm{D}}=& \frac{P_{{\mathrm{R}}_n}{\left|h_{{\mathrm{R}}_n\mathrm{D}}\right|}^2}{{\mathrm{N}}_0}.\hfill \end{array}$$

(8)

2.5. Transmit Power at Source and Relay Nodes

In the underlay cognitive radio networks, all secondary transmitters, i.e., the source node S and all relay R, have to adjust their transmit power to strictly satisfy the interference power threshold denoted by $I_{\mathrm{P}}$ (in Watt) on the primary networks, i.e., the primary node P. As a result, the transmit power of S is then given as ([15], Equation (5))

$$\begin{array}{c}\hfill P_{\mathrm{S}}=\frac{I_{\mathrm{P}}}{{\left|h_{\mathrm{SP}}\right|}^2},\end{array}$$

(9)

Regarding the transmit power of R${}_n$, we have

$$\begin{array}{c}\hfill P_{{\mathrm{R}}_n}=\frac{I_{\mathrm{P}}}{{\left|h_{{\mathrm{R}}_n\mathrm{P}}\right|}^2},\end{array}$$

(10)

Additionally, the transmit power of the relay is also constrained by the amount of the harvested energy in the first phase and is formulated as ([22], Equation (2))

$$\begin{array}{c}\hfill P_{{\mathrm{R}}_n}=\frac{E_{{\mathrm{R}}_n}}{T/2}=\eta \rho P_{\mathrm{S}}{\left|h_{\mathrm{S}{\mathrm{R}}_n}\right|}^2.\end{array}$$

(11)

As a consequence, $P_{{\mathrm{R}}_n}$ can be rewritten as follows

$$\begin{array}{c}\hfill P_{{\mathrm{R}}_n}=min\left(\frac{I_{\mathrm{P}}}{{\left|h_{{\mathrm{R}}_n\mathrm{P}}\right|}^2},\eta \rho P_{\mathrm{S}}{\left|h_{\mathrm{S}{\mathrm{R}}_n}\right|}^2\right)=min\left(\frac{I_{\mathrm{P}}}{{\gamma }_{{\mathrm{R}}_n\mathrm{P}}},\eta \rho P_{\mathrm{S}}{\gamma }_{\mathrm{S}{\mathrm{R}}_n}\right)\stackrel{\left(a\right)}{=}{I}_{\mathrm{P}}min\left(\frac{1}{{\gamma }_{{\mathrm{R}}_n\mathrm{P}}},\frac{\eta \rho {\gamma }_{\mathrm{S}{\mathrm{R}}_n}}{{\gamma }_{\mathrm{SP}}}\right),\end{array}$$

(12)

where $\left(a\right)$ is obtained by substituting $P_{\mathrm{S}}$ in (9).

2.6. End-to-End Signal-to-Noise Ratios at D

Since the decode and forward (DF) protocol is employed, the e2e SNRs is then computed as

$$\begin{array}{cc}\hfill {\gamma }_{\mathrm{e}2\mathrm{e}}=& min\left\{{\gamma }_{{\mathrm{R}}_n},{\gamma }_{\mathrm{D}}\right\}\hfill \\ \hfill \stackrel{\left(a\right)}{=}& \mathsf{\Psi }min\left\{\frac{(1-\rho ){\gamma }_{\mathrm{S}{\mathrm{R}}_n}}{{\gamma }_{\mathrm{SP}}},min\left(\frac{{\gamma }_{{\mathrm{R}}_n\mathrm{D}}}{{\gamma }_{{\mathrm{R}}_n\mathrm{P}}},\frac{\eta \rho {\gamma }_{\mathrm{S}{\mathrm{R}}_n}{\gamma }_{{\mathrm{R}}_n\mathrm{D}}}{{\gamma }_{\mathrm{SP}}}\right)\right\},\hfill \end{array}$$

(13)

where $\mathsf{\Psi }=\frac{I_{\mathrm{p}}}{{\mathrm{N}}_0}$; $\left(a\right)$ is held by substituting $P_{\mathrm{S}}$ and $P_{{\mathrm{R}}_n}$ in (9) and (12) into (8).

Through direct inspection (13), we observe that the e2e SNR of the considered system is more challenging than other work described in the literature. More precisely, the SNR is the composite of two minimum functions instead of only one. Additionally, the random variables inside these minimum functions are fully correlated as well. As a result, the proposed mathematical framework is novel and more complicated than others.

3. Outage Probability (OP) Analysis

In this section, we investigate one of the most important metrics of a wireless communications system, namely, the outage probability which measures the quality-of-service of the whole network. The OP referrs to the probability that the e2e SNRs at D is below a predefined threshold. Mathematically speaking, it is formulated as follows ([12], Equation (23)):

$$\begin{array}{c}\hfill \mathrm{OP}=Pr\left\{{\gamma }_{\mathrm{e}2\mathrm{e}}<{\gamma }_{th}\right\}=Pr\left\{min({\gamma }_{{\mathrm{R}}_n},{\gamma }_{\mathrm{D}})<{\gamma }_{th}\right\},=1-Pr\left\{{\gamma }_{{\mathrm{R}}_n}\ge {\gamma }_{th},{\gamma }_{\mathrm{D}}\ge {\gamma }_{th}\right\},\end{array}$$

(14)

where ${\gamma }_{th}=2^{2\mathcal{R}}-1$, and $\mathcal{R}$ is the targeted rate [in bps/Hz]. In order to compute OP in (14), we first derive Lemma 1 as follows:

Lemma 1. 

Given $\mathcal{N}$ independent and identically distributed (i.i.d.) exponential random variables (RVs) with parameters Ω denoted by $Y_m,m\in \left\{1,\dots ,\mathcal{N}\right\}$. The CDF and PDF of the maximal RV denoted by $Y_{max}=\underset{m\in \left\{1,\dots ,\mathcal{N}\right\}}{max}\left\{Y_m\right\}$ are given as follows:

$$\begin{array}{cc}\hfill F_{Y_{max}}\left(x\right)=& 1+\sum _{m=1}^{\mathcal{N}}{\left(-1\right)}^m{C}_{\mathcal{N}}^{m}exp\left(-m\mathsf{\Omega }x\right)\hfill \\ \hfill f_{Y_{max}}\left(x\right)=& \mathsf{\Omega }\sum _{m=1}^{\mathcal{N}-1}{\left(-1\right)}^m{C}_{\mathcal{N}-1}^{m}exp\left(-\left(m+1\right)\mathsf{\Omega }x\right)\hfill \end{array}$$

(15)

where $C_{\mathcal{N}}^k=\frac{\mathcal{N}!}{k!(\mathcal{N}-k)!}$ is the binomial coefficient.

Proof. 

Let us begin with the definition of the CDF as follows:

$$\begin{array}{cc}\hfill F_{Y_{max}}\left(x\right)=& Pr\left\{Y_{max}=\underset{m\in \left\{1,\dots ,\mathcal{N}\right\}}{max}\left\{Y_m\right\}<x\right\}\hfill \\ \hfill \stackrel{\left(a\right)}{=}& \prod _{m=1}^{\mathcal{N}}{F}_{Y_m}\left(x\right)\stackrel{\left(b\right)}{=}{\left(1-exp\left(-\mathsf{\Omega }x\right)\right)}^{\mathcal{N}}\hfill \\ \hfill \stackrel{\left(c\right)}{=}& 1+\sum _{m=1}^{\mathcal{N}}{\left(-1\right)}^m{C}_{\mathcal{N}}^{m}exp\left(-m\mathsf{\Omega }x\right),\hfill \end{array}$$

(16)

where $\left(a\right)$ is held owing to the independence property between RVs; $\left(b\right)$ is attained by yielding the CDF of $Y_m$; and $\left(c\right)$ is achieved with the help of the binomial theorem. Taking the first-order derivative of the CDF with respect to x, we attain the PDF as follows:

$$\begin{array}{cc}\hfill f_{Y_{max}}\left(x\right)=& \frac{\partial F_{Y_{max}}\left(x\right)}{\partial x}=\mathsf{\Omega }\sum _{m=1}^{\mathcal{N}-1}{\left(-1\right)}^m{C}_{\mathcal{N}-1}^{m}exp\left(-\left(m+1\right)\mathsf{\Omega }x\right).\hfill \end{array}$$

(17)

We close the proof here. □

Next, the OP in (14) is rewritten as follows:

$$\begin{array}{cc}\hfill \mathrm{OP}\stackrel{\left(a\right)}{=}& 1-Pr\left\{\frac{(1-\rho ){\gamma }_{\mathrm{S}{\mathrm{R}}_n}\mathsf{\Psi }}{{\gamma }_{\mathrm{SP}}}\ge {\gamma }_{th},\frac{\mathsf{\Psi }{\gamma }_{{\mathrm{R}}_n\mathrm{D}}}{{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}\ge {\gamma }_{th},\frac{\mathsf{\Psi }\eta \rho {\gamma }_{\mathrm{S}{\mathrm{R}}_n}{\gamma }_{{\mathrm{R}}_n\mathrm{D}}}{{\gamma }_{\mathrm{SP}}}\ge {\gamma }_{th}\right\}\hfill \\ \hfill =& 1-\underset{\Xi }{\underbrace{Pr\left\{X\ge \frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }},{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }},X\ge \frac{{\gamma }_{th}}{\mathsf{\Psi }\eta \rho {\gamma }_{{\mathrm{R}}_n\mathrm{D}}}\right\}}},\hfill \end{array}$$

(18)

where $\left(a\right)$ is attained by substituting (13) into (14), $X=\frac{{\gamma }_{\mathrm{S}{\mathrm{R}}_n}}{{\gamma }_{\mathrm{SP}}}.$ $\Xi$ in (18) can be calculated as follows:

$$\begin{array}{cc}\hfill \Xi =& \underset{{\Xi }_1}{\underbrace{Pr\left\{X\ge \frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }},{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }},\frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }}\ge \frac{{\gamma }_{th}}{\mathsf{\Psi }\eta \rho {\gamma }_{{\mathrm{R}}_n\mathrm{D}}}\right\}}}\hfill \\ \hfill & +\underset{{\Xi }_2}{\underbrace{Pr\left\{X\ge \frac{{\gamma }_{th}}{\mathsf{\Psi }\eta \rho {\gamma }_{{\mathrm{R}}_n\mathrm{D}}},{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }},\frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }}<\frac{{\gamma }_{th}}{\mathsf{\Psi }\eta \rho {\gamma }_{{\mathrm{R}}_n\mathrm{D}}}\right\}}}.\hfill \end{array}$$

(19)

${\Xi }_1$ in (19) can be split into two probabilities, i.e., ${\Xi }_{11}$ and ${\Xi }_{12}$, as follows:

$$\begin{array}{cc}\hfill {\Xi }_1=& Pr\left\{X\ge \frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }},{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }},{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{(1-\rho )}{\eta \rho }\right\}\hfill \\ \hfill =& \underset{{\Xi }_{11}}{\underbrace{Pr\left\{X\ge \frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }},{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }},\frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }}\ge \frac{(1-\rho )}{\eta \rho }\right\}}}\hfill \\ \hfill & +\underset{{\Xi }_{12}}{\underbrace{Pr\left\{X\ge \frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }},{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{(1-\rho )}{\eta \rho },\frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }}<\frac{(1-\rho )}{\eta \rho }\right\}}},\hfill \end{array}$$

(20)

where ${\Xi }_{11}$ is evaluated as

$$\begin{array}{c}\hfill {\Xi }_{11}=\underset{{\mathrm{P}}_1}{\underbrace{Pr\left\{X\ge \frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }}\right\}}}\times \underset{{\mathrm{P}}_2}{\underbrace{Pr\left\{{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }},\frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }}\ge \frac{(1-\rho )}{\eta \rho }\right\}}},\end{array}$$

(21)

Looking at (21), to compute ${\Xi }_{11}$, we first compute ${\mathrm{P}}_1$ as follows:

$$\begin{array}{cc}\hfill {\mathrm{P}}_1=& 1-Pr\left(\frac{{\gamma }_{\mathrm{S}{\mathrm{R}}_n}}{{\gamma }_{\mathrm{SP}}}<\frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }}\right)=1-\underset{0}{\overset{\infty }{\int }}{F}_{{\gamma }_{\mathrm{S}{\mathrm{R}}_n}}\left(\frac{{\gamma }_{th}x}{(1-\rho )\mathsf{\Psi }}\right)\times f_{{\gamma }_{\mathrm{SP}}}\left(x\right)dx,\hfill \\ \hfill \stackrel{\left(a\right)}{=}& \sum _{k=1}^M{\left(-1\right)}^{k+1}{C}_M^k\underset{0}{\overset{\infty }{\int }}{\lambda }_{\mathrm{SP}}\exp \left(-x\left[\frac{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }}+{\lambda }_{\mathrm{SP}}\right]\right)dx\hfill \\ \hfill =& \sum _{k=1}^M\frac{{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{SP}}(1-\rho )\mathsf{\Psi }}{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}+{\lambda }_{\mathrm{SP}}(1-\rho )\mathsf{\Psi }},\hfill \end{array}$$

(22)

where $\left(a\right)$ is achieved by utilizing Lemma 1. Next, ${\mathrm{P}}_2$ is calculated as

$$\begin{array}{cc}\hfill {\mathrm{P}}_2=& Pr\left\{\frac{(1-\rho )\mathsf{\Psi }}{\eta \rho {\gamma }_{th}}\le {\gamma }_{{\mathrm{R}}_n\mathrm{P}}\le \frac{{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\mathsf{\Psi }}{{\gamma }_{th}}\right\}=\underset{\frac{(1-\rho )\mathsf{\Psi }}{\eta \rho {\gamma }_{th}}}{\overset{\infty }{\int }}{f}_{{\gamma }_{{\mathrm{R}}_n\mathrm{D}}}\left(x\right)dx\underset{\frac{(1-\rho )\mathsf{\Psi }}{\eta \rho {\gamma }_{th}}}{\overset{\frac{x\mathsf{\Psi }}{{\gamma }_{th}}}{\int }}{f}_{{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}\left(y\right)dy\hfill \\ \hfill =& {\lambda }_{\mathrm{RD}}\underset{\frac{(1-\rho )\mathsf{\Psi }}{\eta \rho {\gamma }_{th}}}{\overset{\infty }{\int }}exp\left(-{\lambda }_{\mathrm{RD}}x\right)\left\{exp\left(-\frac{(1-\rho )\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{\eta \rho {\gamma }_{th}}\right)-exp\left(-\frac{x\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{{\gamma }_{th}}\right)\right\}dx\hfill \\ \hfill =& exp\left(-\frac{(1-\rho )\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{\eta \rho {\gamma }_{th}}\right)\underset{\frac{(1-\rho )\mathsf{\Psi }}{\eta \rho {\gamma }_{th}}}{\overset{\infty }{\int }}{\lambda }_{\mathrm{RD}}exp\left(-{\lambda }_{\mathrm{RD}}x\right)dx\hfill \\ \hfill & -{\lambda }_{\mathrm{RD}}\underset{\frac{(1-\rho )\mathsf{\Psi }}{\eta \rho {\gamma }_{th}}}{\overset{\infty }{\int }}exp\left(-x\left[{\lambda }_{\mathrm{RD}}+\frac{\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{{\gamma }_{th}}\right]\right)dx\hfill \\ \hfill =& exp\left(-\frac{(1-\rho )\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{\eta \rho {\gamma }_{th}}-\frac{(1-\rho )\mathsf{\Psi }{\lambda }_{\mathrm{RD}}}{\eta \rho {\gamma }_{th}}\right)\hfill \\ \hfill & -\frac{{\lambda }_{\mathrm{RD}}}{{\lambda }_{\mathrm{RD}}+\frac{\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{{\gamma }_{th}}}exp\left(-\frac{(1-\rho )\mathsf{\Psi }}{\eta \rho {\gamma }_{th}}\left[{\lambda }_{\mathrm{RD}}+\frac{\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{{\gamma }_{th}}\right]\right),\hfill \end{array}$$

(23)

From (22) and (23), ${\Xi }_{11}$ in (21) is then computed as

$$\begin{array}{cc}\hfill {\Xi }_{11}=& \sum _{k=1}^M\frac{{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{SP}}(1-\rho )\mathsf{\Psi }}{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}+{\lambda }_{\mathrm{SP}}(1-\rho )\mathsf{\Psi }}\hfill \\ \hfill & \times \left[exp\left(-\frac{\zeta \mathsf{\Psi }\left[{\lambda }_{\mathrm{RP}}+{\lambda }_{\mathrm{RD}}\right]}{{\gamma }_{th}}\right)-\frac{exp\left(-\frac{\zeta \mathsf{\Psi }}{{\gamma }_{th}}\left[{\lambda }_{\mathrm{RD}}+\frac{\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{{\gamma }_{th}}\right]\right)}{\left(1+\frac{\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}\right)}\right],\hfill \end{array}$$

(24)

where $\zeta =\frac{1-\rho }{\eta \rho }$. Having obtained the ${\Xi }_{11}$, we now move to ${\Xi }_{12}$. Let us first rewrite ${\Xi }_{12}$ as

$$\begin{array}{cc}\hfill {\Xi }_{12}=& Pr\left\{X\ge \frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }},{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{(1-\rho )}{\eta \rho },\frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }}<\frac{(1-\rho )}{\eta \rho }\right\}\hfill \\ \hfill =& Pr\left\{X\ge \frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }}\right\}\times Pr\left\{{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{(1-\rho )}{\eta \rho }\right\}\times Pr\left\{{\gamma }_{{\mathrm{R}}_n\mathrm{P}}<\frac{\mathsf{\Psi }(1-\rho )}{\eta \rho {\gamma }_{th}}\right\}\hfill \\ \hfill =& \left\{1-F_X\left(\frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }}\right)\right\}\times \left\{1-F_{{\gamma }_{{\mathrm{R}}_n\mathrm{D}}}\left(\frac{(1-\rho )}{\eta \rho }\right)\right\}\times F_{{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}\left(\frac{\mathsf{\Psi }(1-\rho )}{\eta \rho {\gamma }_{th}}\right),\hfill \\ \hfill \stackrel{\left(a\right)}{=}& \sum _{k=1}^M\frac{{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{SP}}(1-\rho )\mathsf{\Psi }}{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}+{\lambda }_{\mathrm{SP}}(1-\rho )\mathsf{\Psi }}exp\left(-{\lambda }_{\mathrm{RD}}\zeta \right)\left\{1-exp\left(-\frac{{\lambda }_{\mathrm{RP}}\mathsf{\Psi }\zeta }{{\gamma }_{th}}\right)\right\}.\hfill \end{array}$$

(25)

Here, $\left(a\right)$ is held by employing the results of (22). Having ${\Xi }_{11}$ and ${\Xi }_{12}$ in hand, ${\Xi }_1$ in (20) is then computed as

$$\begin{array}{cc}\hfill {\Xi }_1=& \sum _{k=1}^M\frac{{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{SP}}(1-\rho )\mathsf{\Psi }}{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}+{\lambda }_{\mathrm{SP}}(1-\rho )\mathsf{\Psi }}\left[exp\left(-\frac{\zeta \mathsf{\Psi }\left[{\lambda }_{\mathrm{RP}}+{\lambda }_{\mathrm{RD}}\right]}{{\gamma }_{th}}\right)\right.\hfill \\ \hfill & \left.-\frac{exp\left(-\frac{\zeta \mathsf{\Psi }}{{\gamma }_{th}}\left[{\lambda }_{\mathrm{RD}}+\frac{\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{{\gamma }_{th}}\right]\right)}{\left(1+\frac{\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}\right)}+exp\left(-{\lambda }_{\mathrm{RD}}\zeta \right)-exp\left(-\zeta \left[{\lambda }_{\mathrm{RD}}+\frac{{\lambda }_{\mathrm{RP}}\mathsf{\Psi }}{{\gamma }_{th}}\right]\right)\right].\hfill \end{array}$$

(26)

After obtaining ${\Xi }_1$, we now compute ${\Xi }_2$ in (19) as follows:

$$\begin{array}{cc}\hfill {\Xi }_2=& Pr\left\{X\ge \frac{{\gamma }_{th}}{\mathsf{\Psi }\eta \rho {\gamma }_{{\mathrm{R}}_n\mathrm{D}}},{\gamma }_{{\mathrm{R}}_n\mathrm{D}}\ge \frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }},\frac{{\gamma }_{th}}{(1-\rho )\mathsf{\Psi }}<\frac{{\gamma }_{th}}{\mathsf{\Psi }\eta \rho {\gamma }_{{\mathrm{R}}_n\mathrm{D}}}\right\}\hfill \\ \hfill =& Pr\left\{X\ge \frac{{\gamma }_{th}}{\mathsf{\Psi }\eta \rho {\gamma }_{{\mathrm{R}}_n\mathrm{D}}},\frac{{\gamma }_{th}{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}{\mathsf{\Psi }}\le {\gamma }_{{\mathrm{R}}_n\mathrm{D}}<\zeta \right\}\hfill \\ \hfill =& \underset{0}{\overset{\frac{\zeta \mathsf{\Psi }}{{\gamma }_{th}}}{\int }}{f}_{{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}\left(x\right)dx\underset{\frac{{\gamma }_{th}x}{\mathsf{\Psi }}}{\overset{\zeta }{\int }}{f}_{{\gamma }_{{\mathrm{R}}_n\mathrm{D}}}\left(y\right)dy\underset{\frac{{\gamma }_{th}}{\mathsf{\Psi }\eta \rho y}}{\overset{\infty }{\int }}{f}_X\left(z\right)dz,\hfill \\ \hfill \stackrel{\left(a\right)}{=}& \sum _{k=1}^M{\left(-1\right)}^{k+1}{C}_M^k\underset{0}{\overset{\frac{\zeta \mathsf{\Psi }}{{\gamma }_{th}}}{\int }}{f}_{{\gamma }_{{\mathrm{R}}_n\mathrm{P}}}\left(x\right)dx\underset{{\Xi }_{21}}{\underbrace{\underset{\frac{{\gamma }_{th}x}{\mathsf{\Psi }}}{\overset{\zeta }{\int }}\frac{{\lambda }_{\mathrm{SP}}{\lambda }_{\mathrm{RD}}\mathsf{\Psi }\eta \rho y}{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}+{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho y}exp\left(-{\lambda }_{\mathrm{RD}}y\right)dy}}\hfill \end{array}$$

(27)

where $\left(a\right)$ is held with the help of Lemma 1. By employing integration by parts, ${\Xi }_{21}$ is evaluated as

$$\begin{array}{cc}\hfill {\Xi }_{21}=& \left(\frac{\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}exp\left(-\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}\right)}{\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}+\frac{k{\lambda }_{\mathrm{SR}}{\lambda }_{\mathrm{RD}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }}\right)-\left(\frac{{\lambda }_{\mathrm{RD}}\zeta exp\left(-{\lambda }_{\mathrm{RD}}\zeta \right)}{{\lambda }_{\mathrm{RD}}\zeta +\frac{k{\lambda }_{\mathrm{SR}}{\lambda }_{\mathrm{RD}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }}\right)\hfill \\ \hfill & +\underset{{\Xi }_{22}}{\underbrace{\frac{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }\underset{\frac{{\gamma }_{th}x}{\mathsf{\Psi }}}{\overset{\zeta }{\int }}\frac{exp\left(-{\lambda }_{\mathrm{RD}}y\right)}{{\left(y+\frac{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }\right)}^2}dy}}\hfill \\ \hfill \stackrel{\left(a\right)}{=}& \frac{\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}exp\left(-\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}\right)}{\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}+\frac{k{\lambda }_{\mathrm{SR}}{\lambda }_{\mathrm{RD}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }}-\frac{{\lambda }_{\mathrm{RD}}\zeta exp\left(-{\lambda }_{\mathrm{RD}}\zeta \right)}{{\lambda }_{\mathrm{RD}}\zeta +\frac{k{\lambda }_{\mathrm{SR}}{\lambda }_{\mathrm{RD}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }}\hfill \\ \hfill & +\frac{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }exp\left(-\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}\right)\underset{0}{\overset{\zeta -\frac{{\gamma }_{th}x}{\mathsf{\Psi }}}{\int }}\frac{exp\left(-{\lambda }_{\mathrm{RD}}t\right)}{{\left(t+\frac{{\gamma }_{th}x}{\mathsf{\Psi }}+\frac{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }\right)}^2}dt\hfill \\ \hfill \stackrel{\left(b\right)}{=}& \frac{\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}exp\left(-\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}\right)}{\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}+\frac{k{\lambda }_{\mathrm{SR}}{\lambda }_{\mathrm{RD}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }}-\frac{{\lambda }_{\mathrm{RD}}\zeta exp\left(-{\lambda }_{\mathrm{RD}}\zeta \right)}{{\lambda }_{\mathrm{RD}}\zeta +\frac{k{\lambda }_{\mathrm{SR}}{\lambda }_{\mathrm{RD}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }}\hfill \\ \hfill & +\mu exp\left(\mu \right)\left[\mathsf{\Gamma }\left(-1,\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}+\mu \right)-\mathsf{\Gamma }\left(-1,{\lambda }_{\mathrm{RD}}\zeta +\mu \right)\right],\hfill \end{array}$$

(28)

By changing variable $t=y-\frac{{\gamma }_{th}x}{\mathsf{\Psi }}$, we obtain $\left(a\right)$; $\left(b\right)$ is the outcome of ([41], 3.462.17), $\mu =\frac{k{\lambda }_{\mathrm{SR}}{\lambda }_{\mathrm{RD}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }$, and $\mathsf{\Gamma }\left(\alpha ,x\right)=\underset{x}{\overset{\infty }{\int }}{e}^{-t}{t}^{\alpha -1}dt$ is the incomplete gamma function.

Next, substituting (28) into (27), ${\Xi }_2$ can be recomputed as

$$\begin{array}{cc}\hfill {\Xi }_2=& \sum _{k=1}^M{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{RP}}\underset{0}{\overset{\frac{\zeta \mathsf{\Psi }}{{\gamma }_{th}}}{\int }}exp\left(-{\lambda }_{\mathrm{RP}}x\right)dx\hfill \\ \hfill & \times \left[\frac{\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}exp\left(-\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}\right)}{\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}+\frac{k{\lambda }_{\mathrm{SR}}{\lambda }_{\mathrm{RD}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }}-\frac{{\lambda }_{\mathrm{RD}}\zeta exp\left(-{\lambda }_{\mathrm{RD}}\zeta \right)}{{\lambda }_{\mathrm{RD}}\zeta +\frac{k{\lambda }_{\mathrm{SR}}{\lambda }_{\mathrm{RD}}{\gamma }_{th}}{{\lambda }_{\mathrm{SP}}\mathsf{\Psi }\eta \rho }}\right.\hfill \\ \hfill & \left.+\mu exp\left(\mu \right)\left[\mathsf{\Gamma }\left(-1,\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}+\mu \right)-\mathsf{\Gamma }\left(-1,{\lambda }_{\mathrm{RD}}\zeta +\mu \right)\right]\right]\hfill \\ \hfill =& \underset{{\Xi }_{23}}{\underbrace{\sum _{k=1}^M{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{RP}}\underset{0}{\overset{\frac{\zeta \mathsf{\Psi }}{{\gamma }_{th}}}{\int }}\frac{xexp\left(-x\left[{\lambda }_{\mathrm{RP}}+\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}{\mathsf{\Psi }}\right]\right)}{x+\frac{\mu \mathsf{\Psi }}{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}}dx}}\hfill \\ \hfill & -\sum _{k=1}^M\frac{{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{RP}}{\lambda }_{\mathrm{RD}}\zeta exp\left(-{\lambda }_{\mathrm{RD}}\zeta \right)}{{\lambda }_{\mathrm{RD}}\zeta +\mu }\times \left(1-exp\left(-\frac{{\lambda }_{\mathrm{RP}}\zeta \mathsf{\Psi }}{{\gamma }_{th}}\right)\right)\hfill \\ \hfill & +\underset{{\Xi }_{24}}{\underbrace{\sum _{k=1}^M{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{RP}}\mu exp\left(\mu \right)\underset{0}{\overset{\frac{\zeta \mathsf{\Psi }}{{\gamma }_{th}}}{\int }}exp\left(-{\lambda }_{\mathrm{RP}}x\right)\times \mathsf{\Gamma }\left(-1,\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}+\mu \right)dx}}\hfill \\ \hfill & -\sum _{k=1}^M{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{RP}}\mu exp\left(\mu \right)\mathsf{\Gamma }\left(-1,{\lambda }_{\mathrm{RD}}\zeta +\mu \right)\times \left(1-exp\left(-\frac{{\lambda }_{\mathrm{RP}}\zeta \mathsf{\Psi }}{{\gamma }_{th}}\right)\right).\hfill \end{array}$$

(29)

Investigating (29), we observe that ${\Xi }_{24}$, unfortunately, is not able to be computed in closed-form expression due to the fact that the upper limit of the integration does not approach infinity as well as the generality of the incomplete Gamma function. It, however, is effortlessly computed by employing a numerical method with the help of several commercial software such as Matlab, Python, and Mathematica. Moreover, contrary to ${\Xi }_{24}$, ${\Xi }_{23}$ can be evaluated in the closed-form expression by deploying the integration by parts and the assistance of ([41], 3.462.17) as follows:

$$\begin{array}{cc}\hfill {\Xi }_{23}=& \sum _{k=1}^M{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{RP}}\left[exp\left(\mu +\frac{{\lambda }_{\mathrm{RP}}\mu \mathsf{\Psi }}{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}\right)\left[\mathsf{\Gamma }\left(-1,\mu +\frac{{\lambda }_{\mathrm{RP}}\mu \mathsf{\Psi }}{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}\right)\right.\right.\hfill \\ \hfill & \left.\left.-\mathsf{\Gamma }\left(-1,\mu +{\lambda }_{\mathrm{RD}}\zeta +\frac{{\lambda }_{\mathrm{RP}}\mathsf{\Psi }}{{\gamma }_{th}}\left(\zeta +\frac{\mu }{{\lambda }_{\mathrm{RD}}}\right)\right)\right]-\frac{exp\left(-\frac{{\lambda }_{\mathrm{RP}}\zeta \mathsf{\Psi }}{{\gamma }_{th}}-{\lambda }_{\mathrm{RD}}\zeta \right)}{\left(1+\frac{\mu }{{\lambda }_{\mathrm{RD}}\zeta }\right)\left({\lambda }_{\mathrm{RP}}+\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}{\mathsf{\Psi }}\right)}\right].\hfill \end{array}$$

(30)

Finally, substituting (26) and (29) into (19), OP defined in (14) is computed as

$$\begin{array}{cc}\hfill \mathrm{OP}=& 1-{\Xi }_1-{\Xi }_2,\hfill \\ \hfill {\Xi }_1=& \sum _{k=1}^M\frac{{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{SP}}(1-\rho )\mathsf{\Psi }}{k{\lambda }_{\mathrm{SR}}{\gamma }_{th}+{\lambda }_{\mathrm{SP}}(1-\rho )\mathsf{\Psi }}\left[exp\left(-\frac{\zeta \mathsf{\Psi }\left[{\lambda }_{\mathrm{RP}}+{\lambda }_{\mathrm{RD}}\right]}{{\gamma }_{th}}\right)\right.\hfill \\ \hfill & \left.-\frac{exp\left(-\frac{\zeta \mathsf{\Psi }}{{\gamma }_{th}}\left[{\lambda }_{\mathrm{RD}}+\frac{\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{{\gamma }_{th}}\right]\right)}{\left(1+\frac{\mathsf{\Psi }{\lambda }_{\mathrm{RP}}}{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}\right)}+exp\left(-{\lambda }_{\mathrm{RD}}\zeta \right)-exp\left(-\zeta \left[{\lambda }_{\mathrm{RD}}+\frac{{\lambda }_{\mathrm{RP}}\mathsf{\Psi }}{{\gamma }_{th}}\right]\right)\right]\hfill \\ \hfill {\Xi }_2=& \sum _{k=1}^M{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{RP}}\left[exp\left(\mu +\frac{{\lambda }_{\mathrm{RP}}\mu \mathsf{\Psi }}{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}\right)\left[\mathsf{\Gamma }\left(-1,\mu +\frac{{\lambda }_{\mathrm{RP}}\mu \mathsf{\Psi }}{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}\right)\right.\right.\hfill \\ \hfill & \left.\left.-\mathsf{\Gamma }\left(-1,\mu +{\lambda }_{\mathrm{RD}}\zeta +\frac{{\lambda }_{\mathrm{RP}}\mathsf{\Psi }}{{\gamma }_{th}}\left(\zeta +\frac{\mu }{{\lambda }_{\mathrm{RD}}}\right)\right)\right]-\frac{exp\left(-\frac{{\lambda }_{\mathrm{RP}}\zeta \mathsf{\Psi }}{{\gamma }_{th}}-{\lambda }_{\mathrm{RD}}\zeta \right)}{\left(1+\frac{\mu }{{\lambda }_{\mathrm{RD}}\zeta }\right)\left({\lambda }_{\mathrm{RP}}+\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}}{\mathsf{\Psi }}\right)}\right]\hfill \\ \hfill & -\sum _{k=1}^M\frac{{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{RP}}{\lambda }_{\mathrm{RD}}\zeta exp\left(-{\lambda }_{\mathrm{RD}}\zeta \right)}{{\lambda }_{\mathrm{RD}}\zeta +\mu }\times \left(1-exp\left(-\frac{{\lambda }_{\mathrm{RP}}\zeta \mathsf{\Psi }}{{\gamma }_{th}}\right)\right)\hfill \\ \hfill & +\sum _{k=1}^M{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{RP}}\mu exp\left(\mu \right)\underset{0}{\overset{\frac{\zeta \mathsf{\Psi }}{{\gamma }_{th}}}{\int }}exp\left(-{\lambda }_{\mathrm{RP}}x\right)\times \mathsf{\Gamma }\left(-1,\frac{{\gamma }_{th}{\lambda }_{\mathrm{RD}}x}{\mathsf{\Psi }}+\mu \right)dx\hfill \\ \hfill & -\sum _{k=1}^M{\left(-1\right)}^{k+1}{C}_M^k{\lambda }_{\mathrm{RP}}\mu exp\left(\mu \right)\mathsf{\Gamma }\left(-1,{\lambda }_{\mathrm{RD}}\zeta +\mu \right)\times \left(1-exp\left(-\frac{{\lambda }_{\mathrm{RP}}\zeta \mathsf{\Psi }}{{\gamma }_{th}}\right)\right).\hfill \end{array}$$

(31)

4. Simulation Results

In this section, computer-based simulation results are given to verify the correctness of the derived mathematical framework and to discover the behaviors of OP with respect to several essential system parameters, namely, the transmit power, the number of relays, the power splitting ratio, and so forth. In particular, the Monte Carlo method is utilized, and the number of realizations is ${10}^6$ to avoid fluctuation due to insufficient samples. Unless otherwise stated, the following parameters are applied: $\beta =2$, $\eta =0.8$, $\rho$ = 0.5, $\mathsf{\Psi }$ = 5 dB, and $\mathcal{R}$ = 0.5 bps/Hz. All nodes are located on the two-dimensional plane where the source is fixed at the origin, i.e., (0, 0), the positions of other nodes, i.e., relay ${\mathrm{R}}_{\mathrm{n}}$, destination D, and primary P, are (1, 0), (1.5, 0), and (1, 1), respectively.

Figure 3 shows the impact of the interference power over noise variance $\mathsf{\Psi }=\frac{I_p}{N_0}$ on the performance of the OP. The proposed mathematical framework absolutely matches the Monte Carlo simulations. Moreover, we observe that increasing the number of relays M obviously improves the OP. To be more precise, with $\mathsf{\Psi }=30$ dB, the OP of $M=1$ is approximately $5\times {10}^{-3}$, while the OP of $M=4$ is $10\times$ better and is around $5\times {10}^{-4}$. Additionally, the proposed framework outperforms work presented in the literature [42]. In particular, the OP of Ref. [42] only achieves around 0.06 when $\mathsf{\Psi }=30$, while the OP of the proposed scheme with $M=1$ is already approximately 0.005. In addition, the pace of improvement when $\mathsf{\Psi }$ is small and large is different. In particular, the OP dramatically improves when $\mathsf{\Psi }$ is small and a fair improvement is experienced when $\mathsf{\Psi }$ is large. The main reason behind this phenomenon is that when $\mathsf{\Psi }$ is sufficiently large, the outage event almost disappears, as a consequence, a slight enhancement is observed.

Figure 4 shows the outage probability versus the number of relays M. As mentioned in Figure 3, scaling up M certainly ameliorates the OP. Nevertheless, the benefits of increasing M are not the same when M is limited versus M as large. More precisely, when $M<4$, OP declines exponentially, while when $M\ge 4$, OP fairly decreases. Figure 4 confirms again the accuracy of the derived mathematical framework.

We investigate the behavior of OP regarding the power splitting ratio $\rho$ in Figure 5. We see that OP first plunges with $\rho$ after reaching its minimum, OP then turns and constantly rises when $\rho$ moves from 0 to 1. When $\rho$ is small, the amount of harvested energy at R is limited, which leads to lower SNR thereby degrading system performance. When $\rho$ is close to 1, on the other hand, the received signals being added into the information decoder of R are limited too, hence, scaling up the OP. In addition, the curve with $M=4$ is not always superior to the others for all values of $\rho$. In particular, when $\rho$ is small, the curve with $M=2$ is better than the case of $M=4$. However, in general, increasing the number of relays can overcome the constraint of the hardware limitation. Furthermore, the optimal value of $\rho$ denoted by ${\rho }^{*}$ where OP achieves its minimum can be straightforwardly derived from this figure. For example, the ${\rho }^{*}$ of curve with $\eta =0.5$ and $M=4$ is approximately 0.5, while the ${\rho }^{*}$ of the curve with $\eta =0.8$ and $M=2$ is just below 0.4.

Figure 6 stretches the OP versus the targeted rate $\mathcal{R}$. There is no doubt that the larger the $\mathcal{R}$ the higher the OP. We can directly explain this trend by yielding the definition of the OP. Interestingly, the curve with $M=4$ and $\mathsf{\Psi }=1$ dB does not consistently outperform case $M=1$ and $\mathsf{\Psi }=4$ dB. This means that we can ameliorate the system performance either by increasing the number of relay nodes or scaling up the transmit power.

Figure 7 addresses the impact of the mobility of the destination on the performance of the considered networks. To be more precise, destination D is moving from R to far away. The further the $d_{\mathrm{RD}}$ is, the worse the OP. The rationale behind this trend is that increasing $d_{\mathrm{RD}}$ is directly proportional to the large-scale path loss; thus, this deteriorates the channel gain and increases the OP. Moreover, the OP performance can be improved significantly by properly optimizing the power splitting ratio. In particular, we observe that the performance of curves with $M=1$ and $M=2$ is more or less equivalent unless the $d_{\mathrm{RD}}$ approaches zero.

5. Conclusions

The outage performance of the SWIPT-based cognitive radio networks was addressed in this paper. In particular, all relay nodes were equipped not only with a conventional information decoder but also an energy harvester, so that they were able to concurrently receive information and harvest the energy from the incoming signals. Additionally, the opportunistic relaying was taken into account to both ameliorate the outage probability and save the consumed resources. The numerical results illustrated that either increasing the number of relay nodes or raising the transmit power was beneficial for the considered networks. The paper can be extended in several directions: (i) We can replace the fixed relay by high mobility unmanned aerial vehicles (UAVs) devices [43,44] or employ advanced meta material technology to enhance energy efficiency, i.e., the reconfigurable intelligent surfaces (RIS) [45,46]; (ii) Diversity combining at the destination, i.e., maximal ratio combining or selection combining, would also be considered a simple solution to raise the system’s diversity gain [47,48]; (iii) Consideration of the networks could extend to investigating the security aspect as well [49,50]; and (iv) Two-way relaying or full-duplex relaying would be a feasible and interesting extension as well [51,52].