**Enabling Non-Linear Energy Harvesting in Power Domain Based Multiple Access in Relaying Networks: Outage and Ergodic Capacity Performance Analysis**

#### **Thanh-Luan Nguyen 1, Minh-Sang Van Nguyen 2, Dinh-Thuan Do3,***<sup>∗</sup>***, Miroslav Voznak <sup>4</sup>**


Received: 17 June 2019; Accepted: 19 July 2019; Published: 22 July 2019

**Abstract:** The Power Domain-based Multiple Access (PDMA) scheme is considered as one kind of Non-Orthogonal Multiple Access (NOMA) in green communications and can support energy-limited devices by employing wireless power transfer. Such a technique is known as a lifetime-expanding solution for operations in future access policy, especially in the deployment of power-constrained relays for a three-node dual-hop system. In particular, PDMA and energy harvesting are considered as two communication concepts, which are jointly investigated in this paper. However, the dual-hop relaying network system is a popular model assuming an ideal linear energy harvesting circuit, as in recent works, while the practical system situation motivates us to concentrate on another protocol, namely non-linear energy harvesting. As important results, a closed-form formula of outage probability and ergodic capacity is studied under a practical non-linear energy harvesting model. To explore the optimal system performance in terms of outage probability and ergodic capacity, several main parameters including the energy harvesting coefficients, position allocation of each node, power allocation factors, and transmit signal-to-noise ratio (SNR) are jointly considered. To provide insights into the performance, the approximate expressions for the ergodic capacity are given. By matching analytical and Monte Carlo simulations, the correctness of this framework can be examined. With the observation of the simulation results, the figures also show that the performance of energy harvesting-aware PDMA systems under the proposed model can satisfy the requirements in real PDMA applications.

**Keywords:** ergodic capacity; non-linear energy harvesting; NOMA; outage probability

#### **1. Introduction**

As a key subject in the Fifth-Generation (5G) technology, the data rate is enhanced to adapt to the expected 1000-fold explosive growth of data traffic by 2020, and several great challenges are addressed in the development of 5G [1]. To meet the fast increasing traffic need caused by the explosion of wireless devices and improve the spectral effectiveness of the 5G cellular network, a promising technique is proposed and called NOMA [2]. By using the power domain, NOMA serves multiple users simultaneously, and such a technique is different from Orthogonal Multiple Access (OMA), which deploys the time and frequency domains. NOMA can exhibit a balance between network throughput and user fairness because the base station appropriately distributes the transmit power

to multiple users under dissimilar channel conditions. To extend the coverage area and eliminate channel impairments, several relaying system models can be combined with NOMA, in which a relay is reflected as an efficient method to further signal processing in the intermediate node. As a result, the NOMA relaying architecture is introduced to address impairments, including shadowing, path loss, and fading. In dual-hop relaying, a relay operates to assist cooperative communications in which a source communicates with a destination via such a relay [3,4]. In the scenario of downlink or uplink in an NOMA, the user pairing problem was investigated in recent works [5]. On the other hand, other metrics are considered for system performance evaluation in cooperative relaying NOMA; for example, the authors in [6] examined two users, which were paired to establish NOMA based on an access scheme by using classic Decode-and-Forward (DF) relaying with a dual-hop and three-node model. In these models, the Maximal Ratio Combining (MRC) principle is employed to further obtain the space diversity and improve the achievable rate, as presented in [6], or secure performance in the context of NOMA, as achieved in [7]. In contrast with [6] and the full-duplex transmission presented in [8], these analytical expressions can be further extended to the NOMA scheme. As an important metric, a sum-throughput maximization problem was derived and explained with the optimal transmission policy, which is deployed in backscatter-assisted energy harvesting NOMA networks [9]. In [10], stochastic geometry was employed as a tool to evaluate the outage probability of the underlying NOMA-assisted cognitive radio in the scenario of a large-scale model. The cooperative NOMA in full-duplex was considered and surveyed in other findings [11,12]. The outage performance was characterized in full-duplex mode for Device-to-Device (D2D) transmission in cooperative NOMA [11]. Considering the influence of imperfect self-interference, the authors in [12] investigated the expressions of outage probability and the achievable sum rate for full-duplex NOMA relaying wherein DF mode was implemented to serve two NOMA users.

Unlike conventional OMA strategies, the implementation of NOMA may require more energy feeding than OMA, and the spectral efficiency in energy harvesting-based NOMA can always enhance system performance compared to OMA, as the proposed system given in [13]. The explicit evaluations of feasible rate calculation in NOMA-based systems and their prospective performance gains were studied in detail in four kinds of NOMA systems, including coding-assisted NOMA, scrambling-assisted NOMA, spreading-assisted NOMA, and interleaving-assisted NOMA. In addition, these models together with their performances were further compared with OMA [14]. In other trends of research, it was proven that NOMA achieves its prominent advantage with advanced features including the Successive Interference Cancellation (SIC) technique at the receivers, the multiplexing transmitting technique, and effective network resource distribution. The authors in [15] investigated random user deployment in the NOMA scenario to evaluate the network performance. The outage balancing in the downlink NOMA system can be jointly solved in terms of decoding order and the power allocation, as in the recent work [16]. The authors in [17] studied the multiple antenna-based mmwave model to be employed in the NOMA network. Nevertheless, to reduce the joint design of power allocation and decoding order, the fixed precoding beamformer is implemented in the design of NOMA. Unlike these works, the total achievable network rate is maximal to be applied in an uplink NOMA network, and it is expressed by the combined problem of power and sub-carrier allocation [18]. Inspired by [19], which examined a cooperative NOMA system, the work in [20] suggested a better trade-off between spectral efficiency and signal reception in an innovative mixture of downlink and uplink for the NOMA network. Interestingly, visible light communication systems could be able to employ the NOMA, as the potential system presented in [21]. In addition, Multiple-Input-Multiple-Output (MIMO) employing NOMA was introduced in the scenario of multiple antennas equipped at the transmitter or/and receiver, in which all the potential degrees of freedom can be obtained to maximize the performance improvement. To minimize the total expended power or maximize the sum spectral effectiveness, beamforming design in Multiple-Input-Single-Output (MISO) with the capability of the NOMA scheme was studied in [22,23]. The throughput maximization

problem and two algorithms developed were determined as a particular case of a network deploying a two-antenna user [24].

Meanwhile, the ever-growing greenhouse gas emission challenge and explosive proliferation of lower-power devices are the main motivations for developing energy-efficient techniques for future wireless communication networks. In particular, energy-efficient techniques serving power-limited devices (i.e., sensors and mobile phones) can be separated into two main categories. One of the categories focuses on the techniques that can achieve outage probability in the system model where the relay is able to harvest energy from the source [25–29], while the other architecture aims to collect more energy, where the relay can be wirelessly charged via multiple antennas or co-channel interference [30,31], and these concerns are exhibited by the Simultaneous Wireless Information and Power Transfer (SWIPT)-assisted relaying network. In [25], the authors studied energy harvesting policies including a maximum harvested energy relay selection scheme, minimum self-interference relay selection, and an optimal relay selection scheme in a relaying network using full-duplex mode. In a specific condition related to imperfect Channel State Information (CSI), the degraded performance can be observed in a wireless-powered relaying network, as the interesting results presented in [26–28] show. A two-hop relaying network where the terminals and relay are affected by co-channel interference was investigated in [30]. With the existence of co-channel interference, the expressions of the outage probability are derived by considering time switching fractions, aiming at energy harvesting protocols related to making the analysis mathematically tractable [30]. The authors in [31] studied simple calculations in the relaying network, wherein the wireless-powered relay first harvests energy from both the received signal and co-channel interference, then the relay forwards the mixed signal to the destination.

Different from the traditional linear energy harvesting model presented in [30–32], this paper is motivated by a practical non-linear energy harvesting model introduced in [33–35]. Specifically, we focus on relays that are located close to the base station in NOMA systems, and such relays are used as energy harvesting-enabled equipment to help forward signals to the NOMA users at a far distance with poorer channel conditions. It is worth noting that Energy Harvesting (EH) circuits are known as a non-linear wireless power transfer scheme in practice. Therefore, the traditional Linear EH (LEH) models [30–32] are not suitable for making energy transfer enabling wireless relaying maintain its operations. Unfortunately, outage and ergodic analysis strategies following the LEH model may not illustrate the optimal performance in practical systems. In fact, the maximum harvested power value depends on the EH circuits. Thus, it is of great interest to address in this paper the Non-linear EH (NEH) model. For example, to maximize the harvesting power of each energy harvesting receiver, the multi-objective resource optimization problem is considered by using the weighted Tchebycheffmethod [33].

The implementations and contributions of this paper are summarized as follows:


The structure of the remaining parts of this study is arranged as follows. A unified NEH-PDMA framework in Section 2 is studied in wireless communication, where normal cellular users are ordered related to the decoding procedure based on their channel conditions. Next, Section 3 introduces the exact expressions in terms of outage probability to verify the performance of a pair of PDMA users, and these concerned expressions are carefully derived. To further evaluate the insights, Section 4 provides the impact of several parameters on system ergodic capacity. To determine the derived analytical results, Section 5 provides the numerical results. Further remarks and important results are given to conclude our paper in Section 6. Finally, the Appendix shows the proofs of the mathematically-raised problems in the main sections.

**Notation:** This paper needs some main notations to ease the understanding of the upcoming analyses. They are defined as follows: *E* {.} shows the expectation computation. *fX* (.), *FX* (.) denote the Probability Density Function (PDF) and Cumulative Distribution Function (CDF) of a random variable *X*, respectively. *Pr* (.) represents the probability operation. **1**(*C*) denotes the identity function, **1**(*C*) = 1 if *C* holds and **1**(*C*) = 0 otherwise. *Ei*(.) stands for the exponential integral function.

#### **2. System Model**

As illustrated in Figure 1, we considered a scenario where a PDMA-assisted Base Station (BS) wants to communicate with two cell-edge users simultaneously via an EH-assisted relay. The two users were assumed to have a similar channel condition. In addition, the relay utilized the DF relaying protocol to transmit the BS's superposed signal to both users. In addition, the relay was equipped with the EH mechanism using the Power-Splitting (PS) protocol to harvest energy from the received signal. The NEH architecture was deployed in this system. It was assumed that each node was equipped with a single antenna operating in half-duplex mode.

It was assumed that all channel coefficients in this model followed an independent and identical complex Gaussian distribution with zero mean and unit variance. Specifically, we have *g* ∼ CN (0, 1) and *gi* ∼ CN (0, 1) for *i* = 1, 2 as the channel coefficients from the BS to the relay and from the relay to user *i*, respectively.

**Figure 1.** System model of downlink in the Non-linear Energy Harvesting (NEH)-NOMA network.

In some previous works, the decoding order of the SIC receivers was based on users' channel condition, in which the weaker user is served first [6,8,11]. However, in this paper, we assumed a QoS-based decoding order, in which the user with a lower required data rate is served first [36]. Let *Ri* (bits/s/Hz) be the target data rate for user *i*. Hence, without loss of generality, assume that *R*<sup>1</sup> < *R*2.

Assume that the transmission from the BS to both users consumes a duration of *T* block time and is equally allocated for two phases. In the first phase, the BS transmits the superposed signal to the desired relay. Due to PS protocol, the received signal at the information receiver in the relay is given by:

$$y\_R = \sqrt{(1 - E\_R)P\_S} \frac{\mathcal{S}}{\sqrt{d^a}} (\sqrt{a\_1}\mathbf{x}\_1 + \sqrt{a\_2}\mathbf{x}\_2) + n\_{R,} \tag{1}$$

where *ER* ∈ (0, 1) is the power splitting ratio, *PS* is the transmit power of the BS, *d* and *α* denote the distance from the BS to the relay and the path-loss exponent, respectively, *xi* is the unit power information signal of user *<sup>i</sup>* and *ai* for *<sup>i</sup>* ∈ {1, 2} is the power allocation for user *<sup>i</sup>*, and *nR* ∼ CN (0, *<sup>σ</sup>*<sup>2</sup> *R*) denotes the Additive White Gaussian Noise (AWGN) at the relay.

The Signal-to-Interference-plus-Noise Ratios (SINRs) before and after SIC at the relay are given respectively by:

$$\gamma\_{R, \mathbf{x}\_1} = \frac{(1 - E\_R) a\_1 p\_S \ell(d)}{(1 - E\_R) a\_2 p\_S \ell(d) + 1},\tag{2}$$

$$
\gamma\_{\mathbb{R}, \mathbb{x}\_2} = (1 - E\_{\mathbb{R}}) a\_2 p\_S \ell(d), \tag{3}
$$

where - (*d*) <sup>Δ</sup> = |*g*| 2 /*d<sup>α</sup>* and *pS* Δ = *PS*/*σ*<sup>2</sup> *R*.

Due to the nature of NEH, the transmit power at the relay in the linear and non-linear region is given by [35]:

$$P\_R = E\_E E\_R \min\left(\frac{P\_S \lfloor \mathcal{g} \rfloor^2}{d^a}, P\_{th}\right),\tag{4}$$

where *EE* ∈ [0, 1] denotes the EH efficiency depending on the quality of the harvesting circuitry and *Pth* denotes the saturation threshold of EH receiver. Hence, the received signal at user *i* is given as

$$y\_i = \frac{\mathcal{S}\_i}{\sqrt{d\_i^{\text{fl}}}} (\sqrt{b\_1}\mathbf{x}\_1 + \sqrt{b\_2}\mathbf{x}\_2)\sqrt{P\_R} + n\_{i\prime} \tag{5}$$

where *di*, *bi*, and *ni* ∼ CN (0, *<sup>σ</sup>*<sup>2</sup> *<sup>i</sup>* ) denote the distance between the relay and the user *i*, the new power allocation for user *i*, and the AWGN at user *i*, respectively.

At both users, *x*<sup>1</sup> is decoded first due to the assumption *R*<sup>1</sup> < *R*2; thus, the SINR to decode this signal is given by:

$$
\gamma\_{i, \mathbf{x}\_i} = \frac{b\_1 p\_R \ell(d\_i)}{b\_2 p\_R \ell(d\_i) + 1},
\tag{6}
$$

where *pR* Δ = *PR*/*σ*2, *σ*<sup>2</sup> <sup>1</sup> = *<sup>σ</sup>*<sup>2</sup> <sup>2</sup> = *<sup>σ</sup>*2, - (*di*) <sup>Δ</sup> = |*gi*| 2 /*d<sup>α</sup> <sup>i</sup>* . At User 2, SIC is carried out to remove *x*<sup>1</sup> from the received superposed signal, and then, *x*<sup>2</sup> is decoded with the SINR and given by:

$$
\gamma\_{2,x\_2} = b\_2 p\_R \ell(d\_2). \tag{7}
$$

#### **3. Outage Probability Analysis**

Since the capacity of the channel from the BS to the destination user is less than the required transmission rate, an outage event will occur. As a result, in this NEH-PDMA system model, the PDMA users cannot detect the information exactly. In this section, the outage probability is performed as a metric to examine the system performance of unified downlink NEH-PDMA networks. The outage

event of the specific user in the typical cell is that the PDMA user is incapable of performing the signal detection operation. In particular, we characterize the performance in terms of outage probabilities for a pair of PDMA users who receive signal forwarding from an energy harvesting-aware relay as in the following. Before heading to the next section, Equations (6) and (7), i.e., the instantaneous SINR before and possibly after SIC (User 2) can be rewritten as

$$\gamma\_{i, \text{x}\_1} = \frac{a'\_1 E\_R \min\left(p\_S \ell\left(d\right), p\_{th}\right) \ell\left(d\_i\right)}{a'\_2 E\_R \min\left(p\_S \ell\left(d\right), p\_{th}\right) \ell\left(d\_i\right) + 1},\tag{8}$$

$$\gamma\_{2,x\_2} = E\_R a'\_2 \min\left(p\_S \ell\left(d\right), p\_{th}\right) \ell\left(d\_2\right), \tag{9}$$

where *a i* Δ = *biEE*.

#### *3.1. Outage Probability at User 1*

Let *υ<sup>i</sup>* Δ = 22*Ri* − 1, the outage probability at User 1 is the probability an outage event occurs at either the relay or User 1, and it can be formulated as:

$$\begin{split} P\_1 &= \mathbb{P}(\gamma\_{\mathcal{R}, \mathbf{x}\_1} < \upsilon\_1 \text{ or } \gamma\_{\mathcal{R}, \mathbf{x}\_2} < \upsilon\_2) \\ &+ \mathbb{P}(\gamma\_{\mathcal{R}, \mathbf{x}\_1} \ge \upsilon\_1, \gamma\_{\mathcal{R}, \mathbf{x}\_2} \ge \upsilon\_2) \mathbb{P}(\gamma\_{1, \mathbf{x}\_1} < \upsilon\_1 | \gamma\_{\mathcal{R}, \mathbf{x}\_1} \ge \upsilon\_1, \gamma\_{\mathcal{R}, \mathbf{x}\_2} \ge \upsilon\_2) \\ &= 1 - \mathbb{P}(\gamma\_{\mathcal{R}, \mathbf{x}\_1} \ge \upsilon\_1, \gamma\_{\mathcal{R}, \mathbf{x}\_2} \ge \upsilon\_2, \gamma\_{1, \mathbf{x}\_1} \ge \upsilon\_1), \end{split} \tag{10}$$

where (*γR*,*x*<sup>1</sup> < *υ*<sup>1</sup> or *γR*,*x*<sup>2</sup> < *υ*2) is the event that the relay cannot decode either *x*<sup>1</sup> or *x*<sup>2</sup> and (*γ*1,*x*<sup>1</sup> < *υ*1) denotes the event that User 1 cannot decode its own signal. Subsequently, (10) can be evaluated as

$$P\_1 = 1 - \begin{cases} \displaystyle \varepsilon^{-\frac{d^2 t\_1^f}{P\_{th}} - \frac{d^4 p\_{th}}{P\_S}} + \Gamma\left(1, \frac{d^4 t\_{\text{max}}}{P\_S}; \frac{d^4 d\_1^4 t\_1^f}{P\_S}\right) \\ \displaystyle \qquad - \Gamma\left(1, \frac{d^4 p\_{th}}{P\_S}; \frac{d^4 d\_1^4 t\_1^f}{P\_S}\right) \\ \displaystyle \qquad \varepsilon^{-\frac{d^4 p\_1^f}{P\_{th}} - \frac{d^4 t\_{\text{max}}}{P\_S}} \\ \displaystyle \qquad 0 & \text{ $ } \mathcal{I}\_{th} \leq t\_{\text{max}}, \nu\_1 < \upsilon\_{\text{min}} \\ \text{$  } & \text{ $ $  \text{otherwise} } \end{cases} \tag{11}$$

where Γ(*α*, *x*; *b*) = " <sup>∞</sup> *<sup>x</sup> t <sup>α</sup>*−1*e*−*t*−*bt*−<sup>1</sup> *dt* denotes the generalized incomplete gamma (g.i.g.) function [37], *<sup>υ</sup>*min <sup>Δ</sup> = min *<sup>a</sup>*<sup>1</sup> *a*2 , *a* 1 *a* 2 , *<sup>t</sup>*max <sup>Δ</sup> = max (*t*1, *t*2), *t*<sup>1</sup> Δ = *<sup>υ</sup>*<sup>1</sup> (*a*1−*υ*1*a*2)(1−*ER*), *<sup>t</sup>*<sup>2</sup> <sup>Δ</sup> = *<sup>υ</sup>*<sup>2</sup> (1−*ER*)*a*<sup>2</sup> , and *t* 1 Δ = *<sup>υ</sup>*<sup>1</sup> (*a* <sup>1</sup>−*υ*1*a* <sup>2</sup>)*ER* . **Proof.** See Appendix A.

#### *3.2. Outage Probability at User 2*

The outage probability at User 2 is the probability in which an outage event occurs at either the relay or at this user, and it can be formulated as:

$$\begin{split} P\_2 &= \mathbb{P}(\gamma\_{\mathcal{R}, \mathbf{x}\_1} < \nu\_1 \text{ or } \gamma\_{\mathcal{R}, \mathbf{x}\_2} < \nu\_2) \\ &+ \mathbb{P}(\gamma\_{\mathcal{R}, \mathbf{x}\_1} \ge \nu\_1, \gamma\_{\mathcal{R}, \mathbf{x}\_2} \ge \nu\_2) \mathbb{P}[(\gamma\_{\mathcal{Z}, \mathbf{x}\_1} < \nu\_1 \text{ or } \gamma\_{\mathcal{Z}, \mathbf{x}\_2} < \nu\_2) | \gamma\_{\mathcal{R}, \mathbf{x}\_1} \ge \nu\_1, \gamma\_{\mathcal{R}, \mathbf{x}\_2} \ge \nu\_2] \\ &= 1 - \mathbb{P}(\gamma\_{\mathcal{R}, \mathbf{x}\_1} \ge \nu\_1, \gamma\_{\mathcal{R}, \mathbf{x}\_2} \ge \nu\_2, \gamma\_{\mathcal{Z}, \mathbf{x}\_1} \ge \nu\_1, \gamma\_{\mathcal{Z}, \mathbf{x}\_2} \ge \nu\_2), \end{split} \tag{12}$$

in which (*γ*2,*x*<sup>1</sup> < *υ*<sup>1</sup> or *γ*2,*x*<sup>2</sup> < *υ*2) denotes the event that User 2 cannot decode *x*<sup>1</sup> nor decode its own signal after successful SIC. Subsequently, (14) can be expressed as

$$P\_2 = 1 - \begin{cases} c - \frac{d\_s^t \eta\_{\text{max}}}{p\_{th}} - \frac{d^t \eta\_{\text{th}}}{p\_S} + \Gamma\left(1, \frac{d^t \eta\_{\text{max}}}{p\_S}; \frac{d^t d\_1^t \nu\_{\text{max}}'}{p\_S}\right) \\ \qquad - \Gamma\left(1, \frac{d^t \eta\_{\text{ph}}}{p\_S}; \frac{d^t d\_1^t \nu\_{\text{max}}'}{p\_S}\right) \\\qquad \qquad e^{-\frac{d\_s^t \eta\_{\text{max}}}{p\_{th}}} - \frac{d^t \eta\_{\text{max}}}{p\_S} \\\qquad \qquad 0 & \text{,} \end{cases}, p\_{th} > t\_{\text{max}}, \nu\_1 < \nu\_{\text{min}} \tag{13}$$

where *<sup>t</sup>*max <sup>Δ</sup> = max (*t*1, *t*2), *t* max <sup>Δ</sup> = max (*t* 1, *t* <sup>2</sup>), *t* 1 Δ = *<sup>υ</sup>*<sup>1</sup> (*a* <sup>1</sup>−*υ*1*a* <sup>2</sup>)*ER* , and *<sup>t</sup>* 2 Δ = *<sup>υ</sup>*<sup>2</sup> *ERa* 2 **Proof.** This can be rewritten (11) as

$$P\_2 = 1 - \left\{ \begin{aligned} \mathbb{P}\left( p\_S \ell(d) \ge \max(t\_1, t\_2), \ell(d\_2) \ge \frac{\max(t\_1', t\_2')}{\min(p\_S \ell(d), p\_{th})} \right) & \text{, } \upsilon\_1 < \upsilon\_{\min} \\ 0 & \text{, otherwise} \end{aligned} \right. \tag{14}$$

.

where *<sup>t</sup>*<sup>2</sup> <sup>Δ</sup> = *<sup>υ</sup>*<sup>2</sup> (1−*ER*)*a*<sup>2</sup> and *<sup>t</sup>* 2 Δ = *<sup>υ</sup>*<sup>2</sup> *ERa* 2 . When *υ*<sup>1</sup> < *υ*min, it can be achieved:

$$p\_2 = 1 - \mathbb{P}\left(p\_S \ell(d) \ge \max(t\_1, t\_2), \ell(d\_2) \ge \frac{\max(t\_1', t\_2')}{p\_{th}}, p\_S \ell(d) \ge p\_{th}\right)$$

$$- \mathbb{P}\left(p\_S \ell(d) \ge \max(t\_1, t\_2), \ell(d\_2) \ge \frac{\max(t\_1', t\_2')}{p\_S \ell(d)}, p\_S \ell(d) < p\_{th}\right). \tag{15}$$

Similar to (A2), we can obtain (13).

**Proposition 1.** *In order to minimize the outage at the relay, the optimal power allocation at the source is formulated by:*

$$a\_2^\* = \begin{cases} \frac{1}{2} & \text{,} \text{if } (\upsilon\_1 < 1) \text{ and } (\upsilon\_2 \ge \frac{\upsilon\_1}{1 - \upsilon\_1}),\\ \frac{\upsilon\_2}{\upsilon\_1 + \upsilon\_2 + \upsilon\_1 \upsilon\_2} & \text{,} \text{if } \{ (\upsilon\_1 \ge 1) \} \text{ or } \{ (\upsilon\_1 < 1) \text{ and } (\upsilon\_2 < \frac{\upsilon\_1}{1 - \upsilon\_1}) \}. \end{cases} \tag{16}$$

**Proof.** Recall that the coverage probability at the source is equivalent to (*pS*- (*d*) ≥ *t*max). Hence, by minimizing *t*max, the optimal performance at the relay can be achieved. Thus, the optimization problem at the source is formulated as

$$\begin{aligned} a\_2^\* &= \underset{a\_2}{\text{argmin}} \max \left\{ t\_1 \left( a\_2 \right), t\_2 \left( a\_2 \right) \right\} \stackrel{\Delta}{=} t\_{\text{max}} \left( a\_2 \right), \\ &\text{subject} \quad \text{to} \quad a\_2 \in \left[ 0, \min \left\{ \frac{1}{2}, \frac{1}{1 + v\_1} \right\} \right] \end{aligned} \tag{17}$$

in which *t*1(*a*2) = *<sup>υ</sup>*<sup>1</sup> <sup>1</sup>−*ER* <sup>×</sup> <sup>1</sup> <sup>1</sup>−(1+*υ*1)*a*<sup>2</sup> and *<sup>t</sup>*2(*a*2) = *<sup>υ</sup>*<sup>2</sup> <sup>1</sup>−*ER* <sup>×</sup> <sup>1</sup> *a*2 . Note that the points 0 and <sup>1</sup> <sup>1</sup>+*<sup>υ</sup>* (in the case of <sup>1</sup> <sup>1</sup>+*<sup>υ</sup>* <sup>≤</sup> <sup>1</sup> <sup>2</sup> ) are added for the ease of analysis and should be rejected later if *a*<sup>2</sup> takes these values. Otherwise outage will occur at the relay with a probability of 1. For convenience, the problem in (17) can be divided into two cases.

Case 1: <sup>1</sup> <sup>1</sup>+*υ*<sup>1</sup> <sup>≤</sup> <sup>1</sup> <sup>2</sup> , which is equivalent to *<sup>υ</sup>*<sup>1</sup> <sup>≥</sup> 1 and *<sup>a</sup>*<sup>2</sup> <sup>∈</sup> [0, <sup>1</sup> 1+*υ*<sup>1</sup> ]: *t*1(*a*2) monotonically increases from *<sup>υ</sup>*<sup>1</sup> <sup>1</sup>−*ER* to <sup>+</sup>∞, while *<sup>t</sup>*2(*a*2) monotonically decreases from <sup>+</sup><sup>∞</sup> to *<sup>υ</sup>*2(1+*υ*1) <sup>1</sup>−*ER* <sup>&</sup>gt; *<sup>υ</sup>*<sup>1</sup> <sup>1</sup>−*ER* since *<sup>υ</sup>*<sup>2</sup> <sup>&</sup>gt; *<sup>υ</sup>*1; thus, the curves *t*1(*a*2) and *t*2(*a*2) have a point of intersection that is also the global minimum of *<sup>t</sup>*max <sup>=</sup> max{*t*1(*a*2), *<sup>t</sup>*2(*a*2)}. Hence, by solving *<sup>υ</sup>*<sup>1</sup> <sup>1</sup>−(1+*υ*1)*a*<sup>2</sup> <sup>=</sup> *<sup>υ</sup>*<sup>2</sup> *<sup>a</sup>*<sup>2</sup> to find *a*2, the optimal *a*<sup>∗</sup> <sup>2</sup> in this case is obtained as

$$a\_2^\* = \frac{\upsilon\_2}{\upsilon\_1 + \upsilon\_2 + \upsilon\_1 \upsilon\_2}.\tag{18}$$

Case 2: *<sup>υ</sup>*<sup>1</sup> <sup>&</sup>lt; 1, which is equivalent to *<sup>a</sup>*<sup>2</sup> <sup>∈</sup> [0, <sup>1</sup> <sup>2</sup> ]: *<sup>t</sup>*1(*a*2) increases from *<sup>υ</sup>*<sup>1</sup> <sup>1</sup>−*ER* to <sup>2</sup>*υ*<sup>1</sup> (1−*υ*1)(1−*ER*), while *t*2(*a*2) decreases from +∞ to <sup>2</sup>*υ*<sup>2</sup> <sup>1</sup>−*ER* . Further, if *<sup>υ</sup>*<sup>2</sup> <sup>≥</sup> *<sup>υ</sup>*<sup>1</sup> <sup>1</sup>−*υ*<sup>1</sup> <sup>⇒</sup> *<sup>t</sup>*max(*a*2) = *<sup>t</sup>*2(*a*2), which accepts *<sup>a</sup>*<sup>∗</sup> <sup>2</sup> = <sup>1</sup> 2 as its optimal point. For *υ*<sup>2</sup> < *<sup>υ</sup>*<sup>1</sup> 1−*υ*<sup>1</sup> , using the same approach in Case 1, the optimal *a*<sup>2</sup> is obtained by solving *t*1(*a*2) = *t*2(*a*2), resulting in (18). Combining all the above results, we achieve the optimal power allocation as in (16).

**Remark 1.** *In principle, if the signal decoding process for the considered user fails, the outage event is inevitable. It is expected that different outage probabilities can be addressed due to dissimilar power allocation factors for each PDMA user. In other words, the previous outage formulation makes the decoding procedure of the specific user highly dependent on the target rate. Such a concern will be checked in the simulation results.*

#### *3.3. Asymptotic Analysis*

In the high SNR regime, the outage probability at User 1 and User 2 becomes respectively:

$$P\_1^{\infty}(p\_S) \to 1 - \exp\left\{ -\frac{d\_1^n t\_1^\prime}{p\_{th}} - \frac{d^n \max(p\_{th}, t\_{\max})}{p\_S} \right\}, \nu\_1 < \nu\_{\min} \tag{19}$$

$$P\_2^{\infty}(p\_S) \to 1 - \exp\left\{ -\frac{d\_2^u t\_{\text{max}}^{\prime}}{p\_{th}} - \frac{d^u \max(p\_{th\prime}, t\_{\text{max}})}{p\_S} \right\}, \nu\_1 < \nu\_{\text{min}}.\tag{20}$$

Further, the diversity order in terms of outage probability at user *i* is defined as [38]:

$$D \stackrel{\Lambda}{=} - \lim\_{p\_S \to \infty} \frac{\log\_{10} \{ P\_i^{\infty}(p\_S) \}}{\log\_{10} \{ p\_S \}};\tag{21}$$

thus, both User 1 and User 2 have zero diversity order since at *pS* → ∞ and *υ*<sup>1</sup> < *υ*min, the term exp{− *<sup>d</sup><sup>α</sup> pS* max(*pth*, *<sup>t</sup>*max)} can be neglected, resulting in the outage floors at User 1 and User 2 being {1 − *e* − *dα* 1 *t* <sup>1</sup> *pth* } and {<sup>1</sup> <sup>−</sup> *<sup>e</sup>* − *dα* 2 *t* max *pth* }, respectively.

#### **4. Ergodic Capacity**

The ergodic capacity for *xi* is given as *Ci* = <sup>1</sup> <sup>2</sup>E*Zi* [log2(1 + *Zi*)], where we define *Z*<sup>1</sup> Δ = min & *γR*,*x*<sup>1</sup> , *γ*1,*x*<sup>1</sup> , *γ*2,*x*<sup>1</sup> ' and *<sup>Z</sup>*<sup>2</sup> <sup>Δ</sup> = min & *γR*,*x*<sup>2</sup> , *γ*2,*x*<sup>2</sup> ' [32]. In addition, *Ci* can be evaluated analytically as

$$C\_i = \frac{1}{2\ln(2)} \int\_0^\infty \frac{1 - F\_{Z\_i}(x)}{1 + x} dx,\\ (i = 1, 2). \tag{22}$$

#### *4.1. Ergodic Capacity for x*<sup>1</sup>

In order to derive (22) for *i* = 1, the CDF of *Z*<sup>1</sup> is calculated as

$$\begin{split} F\_{\mathbb{Z}\_1}(\gamma) &= 1 - \mathbb{P}(\min(\gamma\_{\mathbb{R}, \mathbf{x}\_1}, \gamma\_{1, \mathbf{x}\_1}, \gamma\_{2, \mathbf{x}\_1}) \ge \gamma) \\ &= 1 - \begin{cases} \mathbb{P}\left(p\_S\ell(d) \ge t\_1(\gamma), \ell(d\_1) \ge \frac{t\_1'(\gamma)}{\min(p\_S\ell(d), p\_{th})}, \ell(d\_2) \ge \frac{t\_1'(\gamma)}{\min(p\_S\ell(d), p\_{th})}\right) & , \gamma < \upsilon\_{\min} \\ & \mathbf{0} & , \text{otherwise} \end{cases} \end{split} \tag{23}$$

where *<sup>t</sup>*1(*γ*) <sup>Δ</sup> = *<sup>γ</sup>* (1−*ER*)(*a*1−*a*2*γ*) and *<sup>t</sup>* <sup>1</sup>(*γ*) <sup>Δ</sup> = *<sup>γ</sup>* (*a* <sup>1</sup>−*a* <sup>2</sup>*γ*)*ER* . It can be evaluated analytically as

$$F\_{\mathcal{Z}\_1}(\gamma) = 1 - \begin{cases} \exp\left(-\frac{d^\alpha p\_{th}}{p\_\mathcal{S}} - (d\_1^a + d\_2^a)\frac{t\_1'(\gamma)}{p\_{th}}\right) + P\_{12}(\gamma) & , \gamma < \min(\frac{d\_1'}{d\_2'}, a\_m) \\\\ \exp\left(-\frac{d^\alpha t\_1(\gamma)}{p\_\mathcal{S}} - (d\_1^a + d\_2^a)\frac{t\_1'(\gamma)}{p\_{th}}\right) & , a\_m \le \gamma < \nu\_{\min} \\\\ 0 & , \text{otherwise} \end{cases} \tag{24}$$

where:

$$P\_{12}(\gamma) = \Gamma\left(1, \frac{d^u t\_1(\gamma)}{p\_S}; (d\_1^u + d\_2^u) \frac{d^u t\_1'(\gamma)}{p\_S}\right) - \Gamma\left(1, \frac{d^u p\_{th}}{p\_S}; (d\_1^u + d\_2^u) \frac{d^u t\_1'(\gamma)}{p\_S}\right). \tag{25}$$

**Proof.** See Appendix B.

Using the Gaussian–Chebyshev quadrature [39], the ergodic capacity of *x*<sup>1</sup> is approximated as

$$\begin{split} \mathbb{C}\_{1} &\approx \sum\_{n=0}^{N} k\_{1,n} \left\{ \exp \left( -\frac{d^{u} p\_{th}}{p\_{\mathcal{S}}} - (d\_{1}^{u} + d\_{2}^{u}) \frac{t\_{1}^{\prime}(\gamma\_{1,n})}{p\_{th}} \right) \\ &- \Gamma \left( 1, \frac{d^{u} p\_{th}}{p\_{\mathcal{S}}}; (d\_{1}^{u} + d\_{2}^{u}) \frac{d^{u} t\_{1}^{\prime}(\gamma\_{1,n})}{p\_{\mathcal{S}}} \right) + \Gamma \left( 1, \frac{d^{u} t\_{1}(\gamma\_{1,n})}{p\_{\mathcal{S}}}; (d\_{1}^{u} + d\_{2}^{u}) \frac{d^{u} t\_{1}^{\prime}(\gamma\_{1,n})}{p\_{\mathcal{S}}} \right) \right\} \\ &+ \mathbf{1} \left( \frac{a\_{1}^{\prime}}{a\_{2}^{\prime}} \ge a\_{m} \right) \sum\_{n=0}^{N} k\_{2,n} \exp \left( -\frac{d^{u} t\_{1}(\gamma\_{2,n})}{p\_{\mathcal{S}}} - (d\_{1}^{u} + d\_{2}^{u}) \frac{t\_{1}^{\prime}(\gamma\_{2,n})}{p\_{th}} \right), \tag{26} \end{split} \tag{26}$$

in which:

$$k\_{i,n} \triangleq \frac{\pi}{2N} \Delta\_i^- |\sin(\frac{2n-1}{N}\pi)| \frac{1}{1+\gamma\_{i,n}},\tag{27}$$

$$\gamma\_{i,n} \triangleq \frac{1}{2}\Delta\_i^+ + \frac{1}{2}\Delta\_i^- \cos(\frac{2n-1}{N}\pi), \quad i \in \{1,2\}, \tag{28}$$

where *am* <sup>Δ</sup> = *<sup>a</sup>*<sup>1</sup> *pth*(1−*ER*) <sup>1</sup>+*a*<sup>2</sup> *pth*(1−*ER*), <sup>Δ</sup><sup>±</sup> 1 Δ = min *<sup>a</sup>* 1 *a* 2 , *am* , Δ± 2 Δ = *υ*min ± *am*, and *N* is a coefficient reflecting the accuracy of the approximation.

#### *4.2. Ergodic Capacity for x*<sup>2</sup>

Similar to (23), the CDF of *Z*<sup>2</sup> is calculated as

$$\begin{split} \mathbb{P}\_{\mathcal{Z}\_2}(\gamma) &= 1 - \mathbb{P}(\min(\gamma\_{\mathcal{R}, \mathbf{x}\_2}, \gamma\_{\mathcal{Z}, \mathbf{x}\_2}) \ge \gamma) \\ &= 1 - \mathbb{P}\left(p\_S \ell(d) \ge t\_2(\gamma), \ell(d\_2) \ge \frac{t\_2'(\gamma)}{\min(p\_S \ell(d), p\_{th})}\right), \end{split} \tag{29}$$

where *<sup>t</sup>*<sup>2</sup> (*γ*) <sup>Δ</sup> = *<sup>γ</sup>* (1−*ER*)*a*<sup>2</sup> and *<sup>t</sup>* <sup>2</sup> (*γ*) <sup>Δ</sup> = *<sup>γ</sup> ERa* 2 .

The above equation can be evaluated analytically as

$$F\_{\mathcal{Z}\_2}(\gamma) = 1 - \begin{cases} \exp\left(-\frac{d^\alpha p\_{\rm th}}{p\_S} - \frac{d\mathcal{Y}\_2^\ell(\gamma)}{p\_{\rm th}}\right) + P\_{\mathcal{Z}2}(\gamma) & , \gamma < p\_{\rm th} a\_2 (1 - E\_{\mathcal{R}})\\ \qquad \exp\left(-\frac{d^\alpha t\_2(\gamma)}{p\_S} - \frac{d\mathcal{Y}\_2^\ell(\gamma)}{p\_{\rm th}}\right) & , \gamma \ge p\_{\rm th} a\_2 (1 - E\_{\mathcal{R}}) \end{cases} \tag{30}$$

where:

$$P\_{22}(\gamma) = \Gamma\left(1, \frac{d^\mu t\_2(\gamma)}{p\_S}; \frac{d^\mu d\_2^\mu t\_2'(\gamma)}{p\_S}\right) - \Gamma\left(1, \frac{d^\mu p\_{th}}{p\_S}; \frac{d^\mu d\_2^\mu t\_2'(\gamma)}{p\_S}\right). \tag{31}$$

**Proof.** Equation (29) can be rewritten as

$$\begin{split} \mathbb{P}\_{\mathcal{Z}\_2}(\gamma) &= 1 - \mathbb{P}\left( p\_S \ell(d) \ge t\_2(\gamma), \ell(d\_2) \ge \frac{t\_2'(\gamma)}{p\_{th}}, p\_S \ell(d) \ge p\_{th} \right) \\ &- \mathbb{P}\left( p\_S \ell(d) \ge t\_2(\gamma), \ell(d\_2) \ge \frac{t\_2'(\gamma)}{p\_S \ell(d)}, p\_S \ell(d) < p\_{th} \right). \end{split} \tag{32}$$

The first probability, denoted as *P*21(*γ*), can be evaluated analytically as

$$\begin{split} P\_{21}(\gamma) &= \int\_{-\max\{t\_2(\gamma), p\_{th}\}}^{\infty} f\_{p\_S\ell(d)}(\mathbf{x}) d\mathbf{x} \int\_{t\_2'(\gamma)/p\_{th}}^{\infty} f\_{\ell(d\_2)}(\mathbf{y}) d\mathbf{y} \\ &= \exp\left(-\frac{d^\alpha}{p\_S} \max(t\_2(\gamma), p\_{th}) - \frac{d\_2^\alpha t\_2'(\gamma)}{p\_{th}}\right), \end{split} \tag{33}$$

whereas the second probability, denoted as *P*22, can be obtained as

$$P\_{22}(\gamma) = \int\_{t\_2(\gamma)}^{p\_{th}} f\_{p \circ \ell(d)}(\mathbf{x}) \int\_{\frac{\ell\_2(\gamma)}{\mathbf{x}}}^{\infty} f\_{\ell(d\_2)}(\mathbf{y}) d\mathbf{x}\_{\prime} \tag{34}$$

and can be calculated by using ([37], Equation (13)). Hence, the proof is the complete.

Further, by applying Gaussian–Chebyshev quadrature [39], the ergodic capacity for *x*<sup>2</sup> can be approximated as

*C*<sup>2</sup> ≈ *N* ∑ *n*=0 *k*3,*<sup>n</sup>* 0 exp  −*d<sup>α</sup> pth pS* − *dα* 2 *t* <sup>2</sup>(*γ*3,*n*) *pth* + Γ 1, *<sup>d</sup>αt*2(*γ*3,*n*) *pS* ; *dαd<sup>α</sup>* 2 *t* <sup>2</sup>(*γ*3,*n*) *pS* − Γ 1, *<sup>d</sup><sup>α</sup> pth pS* ; *dαd<sup>α</sup>* 2 *t* <sup>2</sup>(*γ*3,*n*) *pS* 6 − *e <sup>μ</sup>*Ei(−*μptha*2(<sup>1</sup> <sup>−</sup> *ER*) <sup>−</sup> *<sup>μ</sup>*), (35)

where *<sup>μ</sup>* <sup>Δ</sup> = *<sup>d</sup><sup>α</sup> pSa*2(1−*ER*) <sup>+</sup> *<sup>d</sup><sup>α</sup>* 2 *pthERa* <sup>2</sup> and:

$$k\_{3,n} \stackrel{\Delta}{=} \frac{\pi}{2N} \Delta\_3^- |\sin(\frac{2n-1}{N}\pi)| \frac{1}{1+\gamma\_{3,n}},\tag{36}$$

$$\gamma\_{3,\mathfrak{u}} \stackrel{\Delta}{=} \frac{1}{2} \Delta\_3^+ + \frac{1}{2} \Delta\_3^- \cos\left(\frac{2n-1}{N}\pi\right),\tag{37}$$

where Δ± 3 Δ <sup>=</sup> *ptha*<sup>2</sup> (<sup>1</sup> <sup>−</sup> *ER*). It is noted that in (35), we utilized " <sup>∞</sup> *a e*−*μ<sup>x</sup>* <sup>1</sup>+*<sup>x</sup> dx* <sup>=</sup> <sup>−</sup>*eμ*Ei(−*μ<sup>a</sup>* <sup>−</sup> *<sup>μ</sup>*) ([40], Equation (3.353.5)).

#### *4.3. Asymptotic Analysis*

In the high regime, *P*12(*γ*) → 0 and *P*22(*γ*) → 0. Hence, the analytical results of *C*<sup>1</sup> and *C*<sup>2</sup> with *a*1 *<sup>a</sup>*<sup>2</sup> <sup>=</sup> *<sup>a</sup>* 1 *a* 2 = *<sup>b</sup>*<sup>1</sup> *<sup>b</sup>*<sup>2</sup> can be simplified as

$$\begin{split} \mathbf{C}\_{1}^{\infty} &\rightarrow \frac{e^{-c\_{th}c}}{2\text{In2}} \left\{ \text{Ei}\left(-\frac{c\_{1}+c\_{2}}{a\_{2}}\right) - \text{Ei}\left(-\frac{c\_{1}+c\_{2}}{a\_{2}}(1+a\_{2}c\_{th})\right) \right\} e^{\frac{c\_{1}+c\_{2}}{a\_{2}}} \\ &- \frac{e^{-c\_{th}c}}{2\text{In2}} \left\{ \text{Ei}(-(c\_{1}+c\_{2})) - \text{Ei}(-(c\_{1}+c\_{2})(1+c\_{th})) \right\} e^{c\_{1}+c\_{2}} \\ &+ \frac{1}{2\text{In2}} \left\{ e^{\frac{c+c\_{1}+c\_{2}}{2}} \text{Ei}\left(-\left(\frac{c+c\_{1}+c\_{2}}{a\_{2}}\right)(1+a\_{2}c\_{th})\right) \right. \\ &- e^{c+c\_{1}+c\_{2}} \text{Ei}\left(-\left(c+c\_{1}+c\_{2}\right)(1+c\_{th})\right) \right\}, \end{split}$$

and:

$$\begin{split} \mathbb{C}\_{2}^{\infty} &\rightarrow \frac{1}{2\text{ln}2} \exp\left(-c \times c\_{lh} + \frac{c\_{2}}{a\_{2}}\right) \left[ \text{Ei}\left(-\frac{c\_{2}}{a\_{2}}(1+a\_{2}c\_{lh})\right) - \text{Ei}\left(-\frac{c\_{2}}{a\_{2}}\right) \right] \\ &- \frac{1}{2\text{ln}2} \exp\left(\frac{c+c\_{2}}{a\_{2}}\right) \text{Ei}\left(-\frac{c+c\_{2}}{a\_{2}}(1+a\_{2}c\_{lh})\right), \end{split} \tag{39}$$

where *c* = *<sup>d</sup><sup>α</sup>* (1−*ER*)×*pS* , *<sup>c</sup>*<sup>1</sup> <sup>=</sup> *<sup>d</sup><sup>α</sup>* 1 *EE*×*ER*×*pth* , *<sup>c</sup>*<sup>2</sup> <sup>=</sup> *<sup>d</sup><sup>α</sup>* 2 *EE*×*ER*×*pth* , and *cth* <sup>=</sup> *pth* <sup>×</sup> (<sup>1</sup> <sup>−</sup> *ER*).

#### **Proof.** See Appendix C.

From the above results, it can be seen that at *pS* → ∞ or equivalently *c* → 0, the ergodic capacities at User 1 and User 2 reach:

$$\mathbf{C}\_{1}^{\infty} \rightarrow -\frac{1}{2\ln 2} \left\{ e^{\varepsilon\_{1} + c\_{2}} \mathrm{Ei}\left(-(c\_{1} + c\_{2})\right) - e^{\frac{c\_{1} + c\_{2}}{a\_{2}}} \mathrm{Ei}\left(-\frac{c\_{1} + c\_{2}}{a\_{2}}\right) \right\},\tag{40}$$

$$\text{C}\_2^{\infty} \rightarrow -\frac{1}{2\text{ln}2} \exp\left(\frac{c\_2}{a\_2}\right) \text{Ei}\left(-\frac{c\_2}{a\_2}\right),\tag{41}$$

respectively.

**Remark 2.** *In the LEH model, the harvested energy is unbounded, i.e., PR* = *EEER PS*|*g*<sup>|</sup> 2 *<sup>d</sup><sup>α</sup> . Subsequently, if the base station transmits with relatively large power so that the effect of path loss is insignificant, the transmit power at the relay is approximated to PS. Hence, in high SNR regime, from (21) and (22), a diversity order of 1 is achieved at both users instead of zero. Further, from (42) and (43), no ceiling capacity is observed in such a model. In addition, in the lower transmit SNR regime, the probability for PS*|*g*<sup>|</sup> 2 *<sup>d</sup><sup>α</sup> to exceed the threshold Pth becomes small, and the transmit power at the relay becomes approximated to that of the LEH model. As a result, the LEH model proves its capability to obtain tractable and accurate results in the low/middle transmit SNR regime. In conclusion, the non-linear energy harvesting model is proven to be more practical for the analysis in the high SNR regime with the drawbacks of mathematical complexity.*

#### **5. Numerical Results**

In this section, we evaluate the system performance by varying the main parameters such as power allocation factors on two PDMA users considering the design of the wireless-powered PDMA system based on the NEH model. For the NEH model parameters, the main parameters were chosen according to [34]. Given the following parameters: *d* + *di* = 1, *α* = 3.8, the SNR is calculated by *P*/*σ*2, while the variance of the additive noise is *σ*<sup>2</sup> = 1.

In Figure 2, we show the outage performance curves versus both R1 and R2. We set *a*<sup>1</sup> = 0.8, *a*<sup>2</sup> = 1 − *a*1, *d* = 0.5, *ER* = 0.6, *EE* = 1, *b*<sup>1</sup> = 0.8, *b*<sup>2</sup> = 1 − *b*1, *SNR* = 20 dB, and *pth* = 10 dB. From Figure 2, we can observe that the performance of User 1 outperformed the performance of User 2 due to the different threshold rate constraints. The curves of outage probabilities versus SNR are illustrated in Figure 3. As can be seen, the simulated curves matched the analytical curves very tightly, highlighting the exactness of the proven expressions. The outage performance on User 1 was

better than that on User 2 due to the different power allocation factors. However, a performance gap was only seen more clearly at high SNR, i.e., SNR > 15 dB. In fact, we know that this is because the outage probability depends on the target data rates, and varying target rates exhibit different performance and reach one in a high data rate region. The interesting point is that outage performance at User 2 remained at a stable level as the SNR level was higher than 30 dB. It can be observed that the different requirements at the receiver or different order in detecting the received signal made the outage performance for the two users become dissimilar, especially in the high SNR condition and fixed target rates. Figure 3 also shows that the outage probability would be worse as the target rates increased, especially as they were approximate to 1 bit/s; the worst cases that can be declared.

**Figure 2.** Outage probability of User 1 and User 2 versus *R*<sup>1</sup> and *R*2.

**Figure 3.** Outage probability of User 1 and User 2 versus SNR.

In Figure 4, it can be seen how the distances of the pair of nodes affected the outage performance. In addition, with the increase of *dSR*, the difference of the outage probability among two users can be seen as *d* ranging from 0.1 to approximate 0.7. However, as *d* was greater than 0.7, these outage performances were definitely the same. The reason is that as the relay was placed very far from the BS, the different distances between the relay and two users were similar, and hence, outage performance was the same as well.

**Figure 4.** Outage probability of User 1 and User 2 versus SNR.

From Figure 5, we can see that with the increase of *pth*, the outage probability of our scheme firstly decreased and then tended to be steady, which is because in our scheme, we considered the circuit sensitivity, so the outage probability would tend to be steady. The outage performances of these users were similar at high *pth*. As a result, the outage probability was the smallest as enough energy was harvested at the relay. In this case, User 1 was sensitive with a small value of *pth*, but also had optimal outage performance with a reasonable selection of *pth*.

**Figure 5.** Outage probability of User 1 and User 2 versus distance.

In Figure 6, we plot the ergodic capacity versus SNR for different power allocation factors for each PDMA user. We set *d* = 0.5, *ER* = 0.6, *EE* = 1, *b*<sup>1</sup> = 0.8, *b*<sup>2</sup> = 1 − *b*1, *R*<sup>1</sup> = 0.5, *R*<sup>2</sup> = 1, and *pth* = 20 dB. Increasing the value of SNR in the concerned range, the ergodic capacity of User 2 grew faster than that of User 1. It is noted that two cases of power allocation factor affected the ergodic capacity performance significantly only at lower SNR, as an SNR lower than 30 dB, and then, ergodic capacity was not changed by the varying power allocation factor, as the SNR was greater than 30 dB. However, such an observation was different for User 1, where the performance gaps of two power allocation factors still existed in the whole SNR regime.

**Figure 6.** Ergodic capacity of User 1 and User 2 versus SNR.

From Figure 7, we can see that with the increase of the distance, the ergodic capacity of the two PDMA users in these schemes increased and tended to a maximal point at a specific location of the node. This comes from the fact that different SNR and then different ergodic capacity were obtained by exploiting channel condition disparity. Moreover, we can see that the optimal ergodic capacity could be achieved by User 2 in the numerical method.

**Figure 7.** Ergodic capacity of User 1 and User 2 versus distance.

The impact of *pth* on ergodic capacity is shown in Figure 8. In particular, we can see that the ergodic capacity of both PDMA users increased with increasing *pth*. The impressive point is that the performance gap regarding the ergodic capacity of two PDMA users was very large at high *pth*. Therefore, we can conclude that the energy harvested at the relay was an important factor for the ergodic capacity of User 2. Such characteristics of this figure can be seen considering the ergodic capacity of User 2 where *pth* affected the related performance more remarkably. Moreover, the observation obtained in these figures verifies the correctness of the derived expressions.

**Figure 8.** Ergodic capacity of User 1 and User 2 versus *pth*.

#### **6. Conclusions**

We focused on the QoS-based ordered decoding scheme for downlink PDMA, assuming an energy constraint relay using the non-linear energy harvesting model. The outage probability and ergodic capacity performance problem were expressed in wirelessly-powered PDMA under the practical NEH model. The exact closed-form expressions of outage and ergodic capacity were proposed to evaluate system performance in the scenario of varying harvested power at a relay in the NEH model. By employing numerical simulation, the outage performance will be optimal by jointly optimizing the location of each node, the saturation threshold of the energy harvesting receiver, and power allocation coefficients. It was shown that there is a tradeoff between the outage and ergodic capacity performance versus SNR, and further study can consider which one is selected to achieve the optimal outage and ergodic capacity. Although the performance achieved under the NEH model may be dependent on the harvested power at the relay, power allocation factors in PDMA still lead to fluctuating performance in the considered conditions. These concerns were clearly illustrated in the simulation results, and it was confirmed that our derived expressions were correct.

**Author Contributions:** T.-L.N. designed the algorithm and performed the theoretical analysis, while D.-T.D. checked the results and wrote the manuscript. M.-S.V.N. implemented the simulation and contributed to the manuscript preparation. M.V. was responsible for formulating the research issues and revised the paper.

**Funding:** This research was funded by the Ministry of Education, Youth and Sport of the Czech Republic under SGSGrant No. SP2019/41 conducted at VSB, Technical University of Ostrava.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

Substituting (2), (3), (8) with *i* = 1 and (9) into (10), we achieve:

$$P\_1 = 1 - \begin{cases} \mathbb{P}\left(p\_S \ell(d) \ge t\_{\text{max}}, \ell(d\_1) \ge \frac{t\_1'}{\min(p\_S \ell(d), p\_{\text{th}})}\right) & , \upsilon\_1 < \upsilon\_{\text{min}}\\ 0 & , \text{otherwise} \end{cases} \tag{A1}$$

Considering the case where *υ*<sup>1</sup> < *υ*min, the above equation can be further expressed as

$$\begin{split} \mathbb{P}\_{1} = 1 - \mathbb{P}\left(p\_{S}\ell(d) \ge t\_{\text{max}}, \ell(d\_{1}) \ge \frac{t\_{1}'}{p\_{th}}, p\_{S}\ell(d) \ge p\_{th}\right) \\ -\mathbb{P}\left(p\_{S}\ell(d) \ge t\_{\text{max}}, \ell(d\_{1}) \ge \frac{t\_{1}'}{p\_{S}\ell(d)}, p\_{S}\ell(d) < p\_{th}\right). \end{split} \tag{A2}$$

The above equation can be analyzed as

$$\begin{split} P\_1 &= 1 - \left(1 - F\_{\ell(d\_1)}\left(\frac{t\_1'}{p\_{th}}\right)\right) \int\_{\max(t\_{\max}, p\_{th})}^{\infty} f\_{p\_S\ell(d)}(\mathbf{x}) d\mathbf{x} \\ &- \mathbf{1}(p\_{th} > t\_{\max}) \int\_{t\_{\max}}^{p\_{th}} \left(1 - F\_{\ell(d\_1)}\left(\frac{t\_1'}{\mathbf{x}}\right)\right) f\_{p\_S\ell(d)}(\mathbf{x}) d\mathbf{x}. \end{split} \tag{A3}$$

The CDF and PDF of -(*x*) are *F*-(*x*)(*γ*) = <sup>1</sup> <sup>−</sup> exp(−*xαγ*) and *<sup>f</sup>*-(*x*)(*γ*) = *<sup>x</sup><sup>α</sup>* exp(−*xαγ*), respectively. Note that *fk*-(*x*)(*γ*) = *<sup>x</sup><sup>α</sup> <sup>k</sup>* exp(−*x<sup>α</sup> <sup>k</sup> γ*). Hence, (A3) can be further derived as

$$\begin{split} P\_1 &= 1 - e^{-\frac{d\_1^n t\_1'}{P\_{th}} - \frac{d^n}{P\_S} \max(t\_{\max}, p\_{th})} \\ &- \mathbf{1} (p\_{th} > t\_{\max}) \left[ \Gamma \left( 1, \frac{d^n t\_{\max}}{p\_S}; \frac{d^n d\_1^n t\_1'}{p\_S} \right) - \Gamma \left( 1, \frac{d^n p\_{th}}{p\_S}; \frac{d^n d\_1^n t\_1'}{p\_S} \right) \right], \end{split} \tag{A4}$$

which is equivalent to (24).

#### **Appendix B**

Considering the case of *γ* < min( *<sup>a</sup>*<sup>1</sup> *a*2 , *a* 1 *a* 2 ), Equation (23) becomes:

$$\begin{split} \mathbb{P}\_{\mathcal{Z}\_{1}}(\gamma) &= 1 - \mathbb{P}\left(p\_{S}\ell(d) \ge t\_{1}(\gamma), \ell(d\_{1}) \ge \frac{t\_{1}'(\gamma)}{p\_{th}}, \ell(d\_{2}) \ge \frac{t\_{1}'(\gamma)}{p\_{th}}, p\_{S}\ell(d) \ge p\_{th}\right) \\ &- \mathbb{P}\left(p\_{S}\ell(d) \ge t\_{1}(\gamma), \ell(d\_{1}) \ge \frac{t\_{1}'(\gamma)}{p\_{S}\ell(d)}, \ell(d\_{2}) \ge \frac{t\_{1}'(\gamma)}{p\_{S}\ell(d)}, p\_{S}\ell(d) < p\_{th}\right). \end{split} \tag{A5}$$

The first probability in (A5), denoted as *P*11(*γ*), can be evaluated analytically as

$$\begin{split} P\_{11}(\gamma) &= \int\_{-\max(t\_1(\gamma), p\_{th})}^{\infty} f\_{p\_S\ell(d)}(\mathbf{x}) \int\_{t\_1'(\gamma)/p\_{th}}^{\infty} f\_{\ell(d\_1)}(y) \int\_{t\_1'(\gamma)/p\_{th}}^{\infty} f\_{\ell(d\_2)}(\mathbf{z}) d\mathbf{x} dyd\mathbf{z} \\ &= e^{-\frac{d^2}{p\_S}\max(t\_1(\gamma), p\_{th})} e^{-\frac{d\_1^2 t\_1'(\gamma)}{p\_{th}}} e^{-\frac{d\_1^2 t\_1'(\gamma)}{p\_{th}}} \end{split} \tag{A6}$$

whereas the second probability, denoted as *P*12(*γ*), is calculated as

$$P\_{12}(\gamma) = \mathbf{1}(p\_{\rm th} > t\_1(\gamma)) \int\_{t\_1(\gamma)}^{p\_{\rm th}} f\_{p\_S \ell(d)}(\mathbf{x}) \int\_{t\_1'(\gamma)/x}^{\infty} f\_{\ell(d\_1)}(\mathbf{y}) \int\_{t\_1'(\gamma)/x}^{\infty} f\_{\ell(d\_2)}(\mathbf{z}) d\mathbf{z} d\mathbf{y} d\mathbf{x}$$

$$= \mathbf{1}(p\_{\rm th} > t\_1(\gamma)) \int\_{t\_1(\gamma)}^{p\_{\rm th}} \frac{d^u}{p\_S} \exp\left(-\frac{d^u}{p\_S} \mathbf{x} - t\_1'(\gamma)(d\_1^u + d\_2^u)\frac{1}{x}\right) d\mathbf{x},\tag{A7}$$

by applying ([37], Equation (13)) and noticing that *p* "*th t*1(*γ*) *f*(*x*)*dx* = " ∞ *t*1(*γ*) *f*(*x*)*dx* − " ∞ *pth f*(*x*)*dx*. Further, it is worth noticing that:

$$
\gamma < \upsilon\_{\min}, p\_{th} > t\_1(\gamma) \Leftrightarrow \gamma < \min\left(\frac{a\_1'}{a\_2'}, a\_m\right) \tag{A8}
$$

$$
\gamma < \upsilon\_{\min} \le t\_1(\gamma) \Leftrightarrow (a\_m \le \gamma < \upsilon\_{\min}) \cap \left(\frac{a'\_1}{a'\_2} \ge a\_m\right). \tag{A9}
$$

#### **Appendix C**

**Proof of Equation (38).** When *a*1/*a*<sup>2</sup> = *a* 1/*a* <sup>2</sup>, (24) can be rewritten as

$$F\_{\mathbf{Z}\_1}(\mathbf{x}) = 1 - \begin{cases} \mathrm{e}^{-\varepsilon \times \varepsilon\_{th}} \exp\left\{ - (\mathbf{c}\_1 + \mathbf{c}\_2) \frac{\mathbf{x}}{a\_1 - \mathbf{x} a\_2} \right\} & , \mathbf{x} < a\_m \\\ \exp\left\{ - (\mathbf{c} + \mathbf{c}\_1 + \mathbf{c}\_2) \frac{\mathbf{x}}{a\_1 - \mathbf{x} a\_2} \right\} & , a\_m \le \mathbf{x} < a\_1/a\_2 \end{cases} \tag{A10}$$

Substituting the above equation into (22) and using the change of variable *t* → *x*/*a*<sup>1</sup> − *xa*2, the ergodic capacity of User 1 in a high SNR regime is then given by:

$$\begin{split} \mathcal{C}\_{1} &\to \frac{e^{-c\varepsilon \times c\_{th}}}{2\text{ln}2} \int\_{0}^{c\_{th}} \frac{a\_{1}}{a\_{2}} \frac{1}{t+1} \frac{1}{t+a\_{2}^{-1}} e^{-(c\_{1}+c\_{2})\times t} dt \\ &\quad + \frac{1}{2\text{ln}2} \int\_{c\_{th}}^{\infty} \frac{a\_{1}}{a\_{2}} \frac{1}{t+1} \frac{1}{t+a\_{2}^{-1}} e^{-(c+c\_{1}+c\_{2})\times t} dt, \end{split} \tag{A11}$$

and by applying the decomposition *<sup>a</sup>*<sup>1</sup> *a*2 1 *<sup>t</sup>*+<sup>1</sup> <sup>1</sup> *t*+*a*−<sup>1</sup> 2 = <sup>1</sup> *<sup>t</sup>*+<sup>1</sup> <sup>−</sup> <sup>1</sup> *t*+*a*−<sup>1</sup> 2 , the first integral part becomes:

$$C\_{11} \stackrel{\Delta}{=} \int\_0^{\varepsilon\_{th}} \left\{ \frac{1}{t+1} - \frac{1}{t+a\_2^{-1}} \right\} e^{-(\varepsilon\_1 + \varepsilon\_2) \times t} dt. \tag{A12}$$

Using the identity ([40], Equation (3.352.1)), *C*<sup>11</sup> can be derived as

$$\begin{split} \mathbb{C}\_{11} &= \varepsilon^{\frac{c\_{1}+c\_{2}}{a\_{2}}} \left[ \mathrm{Ei}\left(-\frac{c\_{1}+c\_{2}}{a\_{2}}\right) - \mathrm{Ei}\left(-\frac{c\_{1}+c\_{2}}{a\_{2}}(1+a\_{2}c\_{th})\right) \right] \\ &+ \varepsilon^{c\_{1}+c\_{2}} \left[ \mathrm{Ei}\left(-c\_{1}-c\_{2}\right) - \mathrm{Ei}\left(-(c\_{1}+c\_{2})(1+c\_{th})\right) \right]. \end{split} \tag{A13}$$

Further, using the same decomposition applied for *C*<sup>11</sup> and with the help of ([40], Equation (3.352.2)), the second integral part is derived as

$$\mathbb{C}\_{12} \stackrel{\Lambda}{=} e^{\frac{\varepsilon + c\_1 + c\_2}{a\_2}} \mathrm{Ei}\left(-\frac{\varepsilon + c\_1 + c\_2}{a\_2} \left(1 + a\_2 c\_{th}\right)\right) - e^{\varepsilon + c\_1 + c\_2} \mathrm{Ei}\left(-\left(\varepsilon + c\_1 + c\_2\right) \left(1 + c\_{th}\right)\right). \tag{A14}$$

Substituting (A14) and (A13) into (A11) completes the proof.

**Proof of Equation (39).** The Equation (30) can be rewritten as

$$F\_{\mathbf{Z}\_2} = 1 - \begin{cases} e^{-c \times \varepsilon\_{\rm th}} \exp\left(-\frac{c\_2}{a\_2} \mathbf{x}\right) & , \mathbf{x} < a\_2 c\_{\rm th} \\ \exp\left(-\frac{c + c\_2}{a\_2} \mathbf{x}\right) & , \mathbf{x} \ge a\_2 c\_{\rm th} \end{cases} \tag{A15}$$

By substituting the above equation into (22), it can be achieved that:

$$\begin{split} \mathbf{C}\_{2} &\to \frac{e^{-c \times \varepsilon\_{th}}}{2 \text{ln} 2} \int\_{0}^{a\_{2} \varepsilon\_{th}} \frac{1}{\mathbf{x} + 1} \exp\left(-\frac{c\_{2}}{a\_{2}} \mathbf{x}\right) d\mathbf{x} \\ &+ \frac{1}{2 \text{ln} 2} \int\_{a\_{2} \varepsilon\_{th}}^{\infty} \frac{1}{\mathbf{x} + 1} \exp\left(-\frac{c + c\_{2}}{a\_{2}} \mathbf{x}\right) d\mathbf{x}. \end{split} \tag{A16}$$

The derivation for the integrals in (A11) can be applied to solve the above integrals to complete the proof.

#### **References**



c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
