Throughput Maximization for the Full-Duplex Two-Way Relay System with Energy Harvesting

Abstract

The full-duplex technique can improve the transmission capacity of the communication systems, and energy harvesting (EH) is a promising operation to prolong the lifespan of a wireless node by utilizing the radio-frequency signals. In this paper, the throughput performance of a full-duplex two-way energy EH capable relay system is investigated. In particular, a practical EH protocol, named the time-switching-based relaying (TSR) protocol, is used for EH and the decode-and-forward (DF) policy for information transmission. The outage probability is successfully obtained, and the corresponding system throughput for TSR protocol can be derived by it. The derived throughput is a function of different system parameters, including the time-switching (TS) ratio, power allocation ratio, and the length of the communication time slot. Meanwhile, the throughput is used to characterize a joint time and power allocation scheme for the system, and we aim to find the optimal time and power allocation to achieve the optimal throughput. Due to the existence of three variables and the integral form of throughput expression, an optimization for the throughput is difficult. However, a modified simulated annealing-based search (SABS) algorithm can be used to optimize the throughput. The modified SABS algorithm overcomes being highly impacted by the initial point, and derives the optimal solution fast. Simulation results show that the analytical throughput expression is related with the TS ratio, power allocation ratio, and the length of the communication time slot. The analytical curve of the throughput matches with the simulated one well, which shows that the obtained analytical system throughput for the TSR protocol is valid. Meanwhile, the proposed modified SABS algorithm could be used to derive accurate throughput when SNR is higher than 10 dB.

1. Introduction

Energy harvesting (EH), which enables network nodes to obtain energy from the controllable radio-frequency (RF) signal, is a promising technology that could significantly extend the lifespans of wireless networks [1]. Considering EH, although it has been shown that simultaneous energy and information transmission is optimal for using the RF energy in wireless systems [2], it is difficult to be implemented due to the limitation of current circuits. Therefore, a practical receiver protocol, named the time switching-based relaying (TSR) protocol is introduced [3,4,5,6]. The TSR protocol allows the relay node to switch either the information processing or the EH for the received RF signals at a time, which can be implemented easier [4].

On the other hand, the full-duplex technique could improve the throughput of communication systems by simultaneously transmitting and receiving signals in the same frequency band [7]. In the two-way relay system, the relay node is key to receive and forward the signal to the destination. If the relay node operates with a full-duplex, the system throughput will be enhanced [8,9]. Therefore, it is of great significance to investigate the full-duplex two-way relay system, where the relay node can perform both EH and full-duplex transmissions.

In this paper, we are interested in the throughput performance of the full-duplex (FD) two-way relay system, where the relay node can use RF signals to harvest energy. The relay node needs to simultaneously perform both full-duplex operation and energy harvesting, and it adopts the TSR protocol and the decode-and-forward (DF) policy for EH and full-duplex transmissions, respectively [10].

2. Related Works

There are a number of works that have considered the relay system and its performance analysis. In [11], though the relay node worked with a half-duplex, the authors presented a generalized approach to the performance analysis of relay-aided communication systems, and an expression of the outage probability that is valid for all fading scenarios. In fact, the outage probability is key to derive the system throughput [3,4,12]. Furthermore, a number of works have considered the full-duplex relay system design without energy harvesting [13,14]. In [13], the outage performance of the full-duplex relay system was investigated. In [14], the full-duplex cooperative cognitive radio network with multiple full-duplex secondary users was analyzed, and the optimal system throughput was derived. Since the EH technique was not used in these works, the system does not need to consider time or power allocation.

With energy harvesting, Refs. [15,16,17] further investigated the full-duplex relay system with the TSR protocol (TS-FDR-I protocol). In these papers, the system adopted the TSR protocol (TS-FDR-I protocol) so that it could optimize the time allocation to improve its own performance [15,16,17], and the system also optimized the power allocation for improving system performance [15,17]. In particular, with the TS-FDR-I protocol, the authors in [15] derived the optimal ergodic outage probability by considering different TS ratio, and the authors in [16] derived the optimal system throughput by optimizing the TS ratio. On the other hand, Refs. [15,17] considered the effect of power allocation, and the optimal ergodic outage probability was derived by optimizing the power allocation ratio.

2.1. Motivation

The aforementioned literature focused on the one-way relaying system, and only considered time allocation or power allocation. Compared with the FD one-way relay system, the outage probability and throughput calculation of the FD two-way relay system is more complicated, as we need to consider both destination nodes simultaneously [18,19]. In [18], the authors showed that the outage probability of FD two-way relay system consists of two complicated integral analytic expressions. Furthermore, a joint optimal transmission and power allocation scheme was proposed to improve the throughput of the relay networks [20], but it is difficult to be applied to the FD two-way relay system. Similarly, the authors in [21] used an alternating optimization method to optimize the outage and throughput performance of the full-duplex cooperative relaying system with EH. Its joint time and power allocation scheme is also difficult to be applied to our system. To the best of our knowledge, for the FD two-way relay system, Ref. [22] considered the system average rate analysis with power allocation, and Ref. [23] considered a joint transmit power and relay two-way beamforming optimization method to derive the system sum rate. Therefore, the throughput performance of the FD two-way EH capable relay system with joint time allocation and power allocation scheme is still to be solved, and we aim to propose a useful method to optimize the system throughput.

2.2. Contributions

In this paper, we give a comprehensive study on the throughput performance of the two-way energy harvesting relay system with full-duplex operation. In particular, we consider the joint time allocation and power allocation to maximize the network throughput of the two-way energy harvesting relay system. A new analytical method to obtain the outage probability is adopted, then the corresponding achievable system throughput is derived to characterize the joint time and power allocation scheme for the system. Based on the analytical expression, it is difficult to determine the optimal power allocation ratio and time allocation ratio because the complicated analytical expression contains three variables (TS ratio, power allocation ratio, and the length of communication time slot) and is in the integral form. For multiple variables optimization, the traditional exhaustive search method has high time complexity. To optimize the system throughput more efficiently, we further propose a simulated annealing based search (SABS) algorithm, and show that the close-to-optimal time and power allocation ratio that controls the trade-off between the energy and the information transmission can be efficiently obtained. The main contributions of this paper are summarized as follows:

• Using the DF policy, the signal-to-interference-plus-noise ratio (SINR) of the two-way system is derived. Therefore, the non-outage probability of the system can be calculated by these SINRs with a new analytical method.
• Using the derived non-outage probability, the system throughput expression is derived, which is further formulated as an optimization problem with three variables.
• A modified SABS algorithm is used to solve the throughput optimization problem. The algorithm is not dependent on the selection of the initial point.
• Simulation results show that the analytical throughput matches well with the simulated throughput, and the biggest gap between the analytical throughput with the simulated throughput is $0.8\%$. Furthermore, the results derived by the modified SABS algorithm is accurate when the SNR varies from 10 dB to 40 dB.

The remainder of this paper is organized as follows: Section 2 introduces the full-duplex system model. Section 3 presents the system throughput and proposes the modified simulated annealing based search algorithm. In Section 4, we present the simulation results. Finally, Section 5 concludes the paper.

3. System Model

A full-duplex two-way EH relaying system is considered, which consists of two terminal nodes S and D, and one relay node R, as shown in Figure 1. The direct path between S and D is ignored. Thus, S and D exchange information with the help of the relay node, which is energy-constrained. The relay node R adopts decode and forward (DF) in the full-duplex transmission. With the energy harvesting technique, node R can harvest energy by the RF signals transmitted from S and D for forwarding information. We assume that the information processing energy consumed by R is negligible, compared with the transmit energy. With the FD technique, the information exchange between S and D can be completed in two time slots.

Let $h_1,h_2,g_1,g_2$ denote the channel gains of the links in the FD two-way relay system, and $f_1,f_2$ denote the channel gains of the self-interference link. A reasonable assumption is that all the channel gains do not change during the transmission block time T, which are independent and identically distributed from one block to the next. Furthermore, we assume that $h_1,h_2,g_1,g_2$ follow Rayleigh distribution and $f_1,f_2$ follow Gaussian distribution.

In this paper, we focus on the TS-FDR-I protocol [17], which can be implemented easily in practice. In the TS-FDR-I protocol, the information exchange is completed with two phases. In the first phase, S and D transmit to R for energy harvesting. In the second phase, it is further divided into two slots for information transmission of S-D and D-S, respectively. Let $\alpha$ and $\theta$ denote the time allocation factor for energy harvesting and the time allocation factor between the two link information transmission, respectively. Thus, the time allocation of TS-FDR-I protocol is given as follows:

$$T\left\{\begin{array}{cc}\alpha T& \mathrm{EH}\mathrm{from}\mathrm{D}\mathrm{and}\mathrm{S}.\\ \\ (1-\alpha )T& \left\{\begin{array}{cc}\theta (1-\alpha )T,& \mathrm{S}\to \mathrm{R}\to \mathrm{D},\\ (1-\theta )(1-\alpha )T,& \mathrm{D}\to \mathrm{R}\to \mathrm{S},\end{array}\right.\end{array}\right.$$

where $\mathrm{S}\to \mathrm{R}\to \mathrm{D}$ and $\mathrm{D}\to \mathrm{R}\to \mathrm{S}$ denote the two full-duplex information transmission processes. For finding notations of this paper more conveniently, the key notations are listed in

Table 1. Notation setting.

Notation	Physical Meaning
$h_1,h_2,g_1,g_2$	channel gains of the two-way links (Rayleigh distribution)
$f_1,f_2$	channel gains of the self-interference channel (Gaussian distribution)
${\widehat{f}}_1,{\widehat{f}}_2$	the residual loop-back interference channel coefficient
$\alpha$	time-switching ratio
$\beta$	power allocation ratio
$\theta$	deciding the length of communication time slot for per phase
T	the total communication time
$x_1\left(t\right)$, $x_2\left(t\right)$	normalized signals transmitted by S and D
$P_1$, $P_2$	transmission power of S and D
$\eta$	the energy conversion efficiency
$r_1\left(t\right)$, $r_2\left(t\right)$	signals transmitted by R
$P_{r_{21}}$, $P_{r_{22}}$	transmission power of R in phase 1 and phase 2
$n_{•}$	noise symbol
${\sigma }_{•}^2$	noise power
${\gamma }_{•}$	the derived SINR
${\gamma }_0$	the required SINR threshold
U	signal transmission rate
$p(·)$	the non-outage probability
$p_{out}$	the outage probability of the system
$\tau$	the achievable system throughput
${\lambda }_{h_1},{\lambda }_{h_2},{\lambda }_{g_1},{\lambda }_{g_2}$	mean values of the exponential random variables

.

3.1. Phase 1 (Energy Harvesting)

In the first phase, S and D transmit normalized signals $x_1\left(t\right)$ and $x_2\left(t\right)$ to R with powers $P_1$ and $P_2$, respectively. Without loss of generality, suppose that $E\left\{{\left|x_1\left(t\right)\right|}^2\right\}=E\left\{{\left|x_2\left(t\right)\right|}^2\right\}=1$. The received signal at R is

$$y_r\left(t\right)=\sqrt{P_1}{h}_1{x}_1\left(t\right)+\sqrt{P_2}{g}_2{x}_2\left(t\right)+{\tilde{n}}_{r_{1,a}}\left(t\right),$$

(1)

where ${\tilde{n}}_{r_{1,a}}\left(t\right)$ is the additive Gaussian noise at the relay node. In Phase 1, the time duration for EH is $\alpha T$. Thus, the harvested energy at node R is given by

$$E_h=\eta \left(P_1{\left|h_1\right|}^2+P_2{\left|g_2\right|}^2\right)\alpha T,$$

(2)

where $\eta$ is the energy conversion efficiency, which depends on the EH technology. Figure 2 shows a intuitive energy-harvesting process.

3.2. Phase 2 (Information Processing)

In this phase, S and D exchange their own information by R. Therefore, relay node R allocates a part of its harvested energy (in phase 1) to help S transmit the signal to D, and the remaining harvested energy of R is used to finish transmitting the signal from D to S. The intuitive information process is showed in Figure 3.

3.2.1. S Transmits Signal to D

In the second phase, S first continues to transmit signal $x_1\left(t\right)$, and R transmits to D with the harvested energy in the previous phase. Notice that the received signal at R not only includes the signal from S, but also includes the self-interference by R. Thus, the received signal at R is given by

$$y_{r_{21}}\left(t\right)=\sqrt{P_1}{h}_1{x}_1\left(t\right)+f_1{r}_1\left(t\right)+{\tilde{n}}_{r_{21,a}}\left(t\right),$$

(3)

where ${\tilde{n}}_{r_{12,a}}\left(t\right)$ is the additive Gaussian noise, and $r_1\left(t\right)$ is the relay transmit signal satisfying $E\left\{\right|r_1\left(t\right){|}^2\}=P_{r_{21}}$. Since the time duration for information transmission from R to D is $(1-\alpha )\theta T$, if part of the harvested energy is allocated to this link (the allocated factor is $\beta$), the transmission power of R can be computed as follows [15]:

$$P_{r_{21}}=\frac{\beta E_h}{(1-\alpha )\theta T}=\frac{\beta \alpha \eta \left(P_1{\left|h_1\right|}^2+P_2{\left|g_2\right|}^2\right)}{(1-\alpha )\theta }.$$

(4)

In the full-duplex relaying system, the relay node usually applies the self-interference cancellation technology to reduce the loop-back interference [24,25]. Therefore, after the interference cancellation and sampling, the post-cancellation sampled signal $y_{r_{21}}\left(k\right)$ is given by

$$\begin{array}{cc}\hfill & y_{r_{21}}\left(k\right)=\sqrt{P_1}{h}_1{x}_1\left(k\right)+{\widehat{f}}_1{r}_1\left(k\right)+n_{r_{21,a}}\left(k\right)+n_{r_{21,c}}\left(k\right)\hfill \\ \hfill & =\sqrt{P_1}{h}_1{x}_1\left(k\right)+{\widehat{f}}_1{x}_1(k-\mu )+n_{r_{21,a}}\left(k\right)+n_{r_{21,c}}\left(k\right),\hfill \end{array}$$

(5)

where ${\widehat{f}}_1$ is the residual loop-back interference channel coefficient caused by imperfect cancellation. Because the relay node uses the DF policy for information transmission, it will first decode the original source signal and then regenerate the signal. Hence, the relay transmit signal can be expressed as $r_1\left(k\right)=x_1(k-\mu )$, where $\mu$ is the processing delay at the relay [26]. Furthermore, let $n_{r_{21}}\left(k\right)$ denote $n_{r_{21,a}}\left(k\right)+n_{r_{21,c}}\left(k\right)$ with a noise power of ${\sigma }_{n_{r_{21}}\left(k\right)}^2$.

The final received signal at D is

$$y_d\left(k\right)=\sqrt{P_{r_{21}}}{h}_2{x}_1(k-\mu )+n_d\left(k\right),$$

(6)

where $n_d\left(k\right)\triangleq n_{d,a}\left(k\right)+n_{d,c}\left(k\right)$ is the noise introduced at D. The power of noise $n_d\left(k\right)$ is denoted by ${\sigma }_{n_d\left(k\right)}^2$.

3.2.2. D Transmits Signal to S

Since the time duration for D to S link is $(1-\theta )(1-\alpha )T$ and all the remaining harvested energy is used by R, the transmission power of R is given by

$$P_{r_{22}}=\frac{(1-\beta )E_h}{(1-\alpha )(1-\theta )T}=\frac{(1-\beta )\alpha \eta \left(P_1{\left|h_1\right|}^2+P_2{\left|g_2\right|}^2\right)}{(1-\alpha )(1-\theta )}.$$

(7)

Following a similar analysis of $x_1\left(t\right)$, we can show that the received signal at R is given by

$$\begin{array}{cc}\hfill & y_{r_{22}}\left(k\right)=\sqrt{P_2}{g}_2{x}_2\left(k\right)+{\widehat{f}}_2{r}_2\left(k\right)+n_{r_{22,a}}\left(k\right)+n_{r_{22,c}}\left(k\right)\hfill \\ \hfill & =\sqrt{P_2}{g}_2{x}_2\left(k\right)+{\widehat{f}}_2{x}_2(k-\mu )+n_{r_{22,a}}\left(k\right)+n_{r_{22,c}}\left(k\right),\hfill \end{array}$$

(8)

where ${\widehat{f}}_2$ is the residual loop-back interference channel coefficient, and $r_2\left(k\right)=x_2(k-\mu )$ ($E\left\{\right|r_2\left(k\right){|}^2\}=P_{r_{22}}$) is the relay transmit signal. In particular, $n_{r_{22}}\left(k\right)\triangleq n_{r_{22,a}}\left(k\right)+n_{r_{22,c}}\left(k\right)$ with a noise power of ${\sigma }_{n_{r_{22}}\left(k\right)}^2$.

After information processing, R forwards the signal to S and the final received signal at S is

$$y_s\left(k\right)=\sqrt{P_{r_{22}}}{g}_1{x}_2(k-\mu )+n_s\left(k\right),$$

(9)

where $n_s\left(k\right)\triangleq n_{s,a}\left(k\right)+n_{s,c}\left(k\right)$ is the noise introduced at S. The power of noise $n_s\left(k\right)$ is denoted by ${\sigma }_{n_s\left(k\right)}^2$.

We assume that ${\sigma }_{n_{r_{21}}\left(k\right)}^2={\sigma }_{n_{r_{22}}\left(k\right)}^2={\sigma }_{n_s\left(k\right)}^2={\sigma }_{n_d\left(k\right)}^2={\sigma }^2$. Therefore, the signal-to-interference-plus-noise ratio (SINR) at S, R and D is given by

$$\begin{array}{c}{\gamma }_{r_{21}}=\frac{P_1{\left|h_1\right|}^2}{|{\widehat{f}}_1{|}^2{P}_{r_{21}}+{\sigma }^2},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & \gamma {}_d=\frac{P_{r_{21}}{\left|h_2\right|}^2}{{\sigma }^2}=\frac{\beta \alpha \eta \left(P_1{\left|h_1\right|}^2+P_2{\left|g_2\right|}^2\right){\left|h_2\right|}^2}{(1-\alpha )\theta {\sigma }^2},\hfill \end{array}$$

(11)

$$\begin{array}{c}{\gamma }_{r_{22}}=\frac{P_2{\left|g_2\right|}^2}{|{\widehat{f}}_2{|}^2{P}_{r_{22}}+{\sigma }^2},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & \gamma {}_s=\frac{P_{r_{22}}{\left|g_1\right|}^2}{{\sigma }^2}=\frac{\beta \alpha \eta \left(P_1{\left|h_1\right|}^2+P_2{\left|g_2\right|}^2\right){\left|g_1\right|}^2}{(1-\alpha )(1-\theta ){\sigma }^2}.\hfill \end{array}$$

(13)

4. Throughput Analysis

In this section, we derive the analytical achievable throughput of the FD two-way system with EH. For the full-duplex system with the DF policy, we first need to derive the expression of outage probability, and then compute the achievable throughput.

4.1. Achievable Throughput

Let ${\gamma }_0$ denote the required SINR threshold for correct reception with transmission rate U, i.e., ${\gamma }_0=2^U-1$. Since the relay uses the DF relaying policy, the outage probability is given by

$$p_{out}=1-p({\gamma }_{r_{21}}\ge {\gamma }_0)p({\gamma }_s\ge {\gamma }_0)p({\gamma }_{r_{22}}\ge {\gamma }_0)p({\gamma }_d\ge {\gamma }_0).$$

(14)

The above equation shows that the outage probability of the full-duplex system in the DF relaying scheme consists of four components. Next, we will show how to compute each component in (14).

The non-outage probability $p({\gamma }_{r_{21}}\ge {\gamma }_0)$ is given by

$$\begin{array}{cc}\hfill & p({\gamma }_{r_{21}}\ge {\gamma }_0)=p\left(\frac{P_1{\left|h_1\right|}^2}{|{\widehat{f}}_1{|}^2{P}_{r_{21}}+{\sigma }^2}\ge {\gamma }_0\right)=p\left(\frac{P_1{\left|h_1\right|}^2}{|{\widehat{f}}_1{|}^2\left(\frac{\beta \alpha \eta \left(P_1{\left|h_1\right|}^2+P_2{\left|h_2\right|}^2\right)}{(1-\alpha )\theta }\right)+{\sigma }^2}\ge {\gamma }_0\right)\hfill \\ \hfill & =p\left(\frac{(1-\alpha )\theta P_1{\left|h_1\right|}^2}{|{\widehat{f}}_1{|}^2\beta \alpha \eta \left(P_1{\left|h_1\right|}^2+P_2{\left|g_2\right|}^2\right)+(1-\alpha )\theta {\sigma }^2}\ge {\gamma }_0\right)\hfill \\ \hfill & =p\left({\left|h_1\right|}^2\ge \frac{{\gamma }_0|{\widehat{f}}_1{|}^2\beta \alpha \eta P_2{\left|g_2\right|}^2+{\gamma }_0(1-\alpha )\theta {\sigma }^2}{(1-\alpha )\theta P_1-{\gamma }_0|{\widehat{f}}_1{|}^2\beta \alpha \eta P_1}\right).\hfill \end{array}$$

(15)

Since $h_1$ and $g_2$ follow Rayleigh distribution, $|h_1{|}^2$ and $|g_2{|}^2$ are exponential random variables with mean values of ${\lambda }_{h_1}$ and ${\lambda }_{g_2}$, respectively. Thus, we have

$$\begin{array}{cc}\hfill & p({\gamma }_{r_{21}}\ge {\gamma }_0)\hfill \\ \hfill & =p\left({\left|h_1\right|}^2\ge \frac{{\gamma }_0|{\widehat{f}}_1{|}^2\beta \alpha \eta P_2{\left|g_2\right|}^2+{\gamma }_0(1-\alpha )\theta {\sigma }^2}{(1-\alpha )\theta P_1-{\gamma }_0|{\widehat{f}}_1{|}^2\beta \alpha \eta P_1}\right)\hfill \\ \hfill & ={\int }_{z=0}^{\infty }{f}_{{\left|g_2\right|}^2}\left(z\right)p\left({\left|h_1\right|}^2\ge \frac{{\gamma }_0|{\widehat{f}}_1{|}^2\beta \alpha \eta P_{2}z+{\gamma }_0(1-\alpha )\theta {\sigma }^2}{(1-\alpha )\theta P_1-{\gamma }_0|{\widehat{f}}_1{|}^2\beta \alpha \eta P_1}\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =\frac{1}{{\lambda }_{g_2}}exp\left(-\frac{{\gamma }_0(1-\alpha )\theta {\sigma }^2}{{\lambda }_{h_1}(1-\alpha )\theta P_1-{\gamma }_0|{\widehat{f}}_1{|}^2\beta \alpha \eta P_1{\lambda }_{h_1}}\right)*\hfill \\ \hfill & \left(\frac{{\lambda }_{g_2}\left({\lambda }_{h_1}(1-\alpha )\theta P_1-{\gamma }_0|{\widehat{f}}_1{|}^2\beta \alpha \eta P_1{\lambda }_{h_1}\right)}{{\lambda }_{h_1}(1-\alpha )\theta P_1-{\gamma }_0|{\widehat{f}}_1{|}^2\beta \alpha \eta P_1{\lambda }_{h_1}+{\lambda }_{g_2}{\gamma }_0|{\widehat{f}}_1{|}^2\beta \alpha \eta P_2}\right).\hfill \end{array}$$

(16)

Furthermore, the non-outage probability $p({\gamma }_{r_{21}}\ge {\gamma }_0)$ can also derive some insights form the work [4,12].

Similarly, the non-outage probability $p({\gamma }_{r_{22}}\ge {\gamma }_0)$ is given by

$$\begin{array}{cc}\hfill & p({\gamma }_{r_{22}}\ge {\gamma }_0)\hfill \\ \hfill & =p\left(\frac{P_2{\left|g_2\right|}^2}{|{\widehat{f}}_2{|}^2{P}_{r_{22}}+{\sigma }^2}\ge {\gamma }_0\right)\hfill \\ \hfill & =p\left(\frac{P_2{\left|g_2\right|}^2}{|{\widehat{f}}_2{|}^2\left(\frac{(1-\beta )\alpha \eta \left(P_1{\left|h_1\right|}^2+P_2{\left|g_2\right|}^2\right)}{(1-\alpha )(1-\theta )}\right)+{\sigma }^2}\ge {\gamma }_0\right)\hfill \\ \hfill & =p\left({\left|g_2\right|}^2\ge \frac{{\gamma }_0|{\widehat{f}}_2{|}^2(1-\beta )\alpha \eta P_1{\left|h_1\right|}^2+{\gamma }_0(1-\alpha )(1-\theta ){\sigma }^2}{(1-\alpha )(1-\theta )P_2-{\gamma }_0|{\widehat{f}}_2{|}^2(1-\beta )\alpha \eta P_2}\right)\hfill \\ \hfill & =\frac{1}{{\lambda }_{h_1}}exp\left(-\frac{{\gamma }_0(1-\alpha )(1-\theta ){\sigma }^2}{{\lambda }_{g_2}(1-\alpha )(1-\theta )P_2-{\gamma }_0|{\widehat{f}}_2{|}^2(1-\beta )\alpha \eta P_1{\lambda }_{g_2}}\right)*\hfill \\ \hfill & \left(\frac{{\lambda }_{h_1}\left({\lambda }_{g_2}(1-\alpha )(1-\theta )P_2-{\gamma }_0|{\widehat{f}}_2{|}^2(1-\beta )\alpha \eta P_2{\lambda }_{g_2}\right)}{{\lambda }_{g_2}(1-\alpha )(1-\theta )P_2-{\gamma }_0|{\widehat{f}}_2{|}^2(1-\beta )\alpha \eta P_2{\lambda }_{g_2}+{\lambda }_{h_1}{\gamma }_0|{\widehat{f}}_2{|}^2(1-\beta )\alpha \eta P_2}\right).\hfill \end{array}$$

(17)

When R forwards signal to D, the non-outage probability at D is given by

$$\begin{array}{cc}\hfill & p({\gamma }_d\ge {\gamma }_0)=p\left(\frac{\beta \alpha \eta \left(P_1{\left|h_1\right|}^2+P_2{\left|g_2\right|}^2\right){\left|h_2\right|}^2}{(1-\alpha )\theta {\sigma }^2}\ge {\gamma }_0\right),\hfill \end{array}$$

(18)

and its analytical expression could be derived by the following proposition.

Proposition 1. 

$\left(1\right)$ If $P_1{\lambda }_{h_1}\ne P_2{\lambda }_{g_2}$, the non-outage probability $p({\gamma }_d\ge {\gamma }_0)$, defined in (18), is given by

$$\begin{array}{cc}\hfill & p({\gamma }_d\ge {\gamma }_0)=\hfill \\ \hfill & \frac{1}{P_1{\lambda }_{h_1}-P_2{\lambda }_{g_2}}\sqrt{\frac{4{\gamma }_0(1-\alpha )\theta {\sigma }^2{P}_1{\lambda }_{h_1}}{{\lambda }_{h_2}\beta \alpha \eta }}{K}_1\left(\sqrt{\frac{4{\gamma }_0(1-\alpha )\theta {\sigma }^2}{{\lambda }_{h_2}\beta \alpha \eta P_1{\lambda }_{h_1}}}\right)\hfill \\ \hfill & -\frac{1}{P_1{\lambda }_{h_1}-P_2{\lambda }_{g_2}}\sqrt{\frac{4{\gamma }_0(1-\alpha )\theta {\sigma }^2{P}_2{\lambda }_{g_2}}{{\lambda }_{h_2}\beta \alpha \eta }}{K}_1\left(\sqrt{\frac{4{\gamma }_0(1-\alpha )\theta {\sigma }^2}{{\lambda }_{h_2}\beta \alpha \eta P_2{\lambda }_{g_2}}}\right),\hfill \end{array}$$

(19)

where the notation $K_n\left(x\right)$ denotes the nth order modified Bessel function of the second kind.

$\left(2\right)$ If $P_1{\lambda }_{h_1}=P_2{\lambda }_{g_2}$, the non-outage probability at D is given by

$$\begin{array}{cc}\hfill & p({\gamma }_d\ge {\gamma }_0)=\frac{2{\gamma }_0(1-\alpha )\theta {\sigma }^2}{{\lambda }_h{}_2\beta \alpha \eta P_1{\lambda }_{h_1}}{K}_2\left(2\sqrt{\frac{{\gamma }_0(1-\alpha )\theta {\sigma }^2}{{\lambda }_h{}_2\beta \alpha \eta P_1{\lambda }_{h_1}}}\right).\hfill \end{array}$$

(20)

Proof. 

See Appendix A.    □

Similarly, the expression of the non-outage probability at S is shown as follows: $\left(1\right)$ If $P_1{\lambda }_{h_1}\ne P_2{\lambda }_{g_2}$,

$$\begin{array}{cc}\hfill & p({\gamma }_s\ge {\gamma }_0)=p\left(\frac{(1-\beta )\alpha \eta \left(P_1{\left|h_1\right|}^2+P_2{\left|h_2\right|}^2\right){\left|g_1\right|}^2}{(1-\alpha )(1-\theta ){\sigma }^2}\ge {\gamma }_0\right)\hfill \\ \hfill & =\frac{1}{P_1{\lambda }_{h_1}-P_2{\lambda }_{g_2}}(\sqrt{\frac{4{\gamma }_0(1-\alpha )(1-\theta ){\sigma }^2{P}_1{\lambda }_{h_1}}{{\lambda }_{g_1}(1-\beta )\alpha \eta }}{K}_1(\sqrt{\frac{4{\gamma }_0(1-\alpha )(1-\theta ){\sigma }^2}{{\lambda }_{g_1}(1-\beta )\alpha \eta P_1{\lambda }_{h_1}}})\hfill \\ \hfill & -\sqrt{\frac{4{\gamma }_0(1-\alpha )(1-\theta ){\sigma }^2{P}_2{\lambda }_{g_2}}{{\lambda }_{g_1}(1-\beta )\alpha \eta }}{K}_1(\sqrt{\frac{4{\gamma }_0(1-\alpha )(1-\theta ){\sigma }^2}{{\lambda }_{g_1}(1-\beta )\alpha \eta P_2{\lambda }_{g_2}}})).\hfill \end{array}$$

(21)

$\left(2\right)$ If $P_1{\lambda }_{h_1}=P_2{\lambda }_{g_2}$,

$$\begin{array}{cc}\hfill & p({\gamma }_d\ge {\gamma }_0)=\frac{2{\gamma }_0(1-\alpha )(1-\theta ){\sigma }^2}{{\lambda }_g{}_1(1-\beta )\alpha \eta P_1{\lambda }_{h_1}}{K}_2\left(2\sqrt{\frac{{\gamma }_0(1-\alpha )(1-\theta ){\sigma }^2}{{\lambda }_g{}_1(1-\beta )\alpha \eta P_1{\lambda }_{h_1}}}\right).\hfill \end{array}$$

(22)

Finally, the achievable system throughput is given by

$$\begin{array}{c}\tau =(1-p_{out})2U\times (1-\alpha )min(\theta ,1-\theta ).\hfill \end{array}$$

(23)

4.2. System Throughput Optimization

In the energy harvesting FD two-way relay system, we aim to optimize the time and power allocation to maximize the system throughput. In particular, it can be formulated as the following optimization problem:

$$\begin{array}{cc}\hfill & max\tau (\alpha ,\beta ,\theta )=(1-p_{out})2U\times (1-\alpha )min(\theta ,1-\theta ).\hfill \\ \hfill & s.t.0\le \alpha \le 1,0\le \beta \le 1,0\le \theta \le 1.\hfill \end{array}$$

(24)

However, problem (24) contains three variables and $p_{out}$ is in the integral form which make problem (24) difficult to solve. The exhaustive search method for all three variables has high time complexity, and this motivates us to propose a simulated annealing based search (SABS) algorithm to obtain the optimal system configuration parameters that maximize system throughput. Note that the SABS algorithm is very dependent on the selection of the initial point. To overcome this problem, we improve its performance using the principle of alternating optimization, i.e, deriving $\alpha$ and $(\beta ,\theta )$ by exhaustive search and SABS, respectively.

Solving problem (24) by the modified SABS algorithm, we first set initial value of $\beta$ and $\theta$ randomly, and obtain the optimal value of ${\alpha }^{\left(j\right)}$ by exhaustive search (see the step 1 of Algorithm 1). Then we construct the vector solution $x=({\alpha }^{\left(j\right)},\beta ,\theta )\in \mathrm{\Omega }$, where $\mathrm{\Omega }$ is the solution space consisting of ${\alpha }^{\left(j\right)}$, $\beta$ and $\theta$. In step 3 of Algorithm 1, we design an accept probability which could help select the alternative point $z^{\left(k\right)}$ as the next iteration point. The accepted probability is given by

$$\begin{array}{c}p(k,\tau \left(z^{\left(k\right)}\right),\tau \left(x^{\left(k\right)}\right))=min\{1,exp(\frac{-(\tau \left(z^{\left(k\right)}\right)-\tau \left(x^{\left(k\right)}\right)}{T_k})\},\hfill \end{array}$$

(25)

where $T_k$ is a sequence of positive numbers, i.e., $T_k$ is cooling schedule [27]. From the accepted probability, we could observe that if $\tau \left(z^{\left(k\right)}\right)\le \tau \left(x^{\left(k\right)}\right)$, the value of $p(k,\tau \left(z^{\left(k\right)}\right),\tau \left(x^{\left(k\right)}\right))$ is equal to 1, which means $x^{(k+1)}=z^{\left(k\right)}$. Furthermore, although $\tau \left(z^{\left(k\right)}\right)>\tau \left(x^{\left(k\right)}\right)$, the probability that $x^{(k+1)}=z^{\left(k\right)}$ is $exp\left(\frac{-(\tau \left(z^{\left(k\right)}\right) - \tau \left(x^{\left(k\right)}\right)}{T_k}\right)$. Generally, cooling schedule $T_k$ monotonically decreases to 0, and this denotes the cooling process. In Algorithm 1, $T_k$ is given by

$$\begin{array}{c}{T}_k=0.{95}^k{T}_0(T_0=100),\hfill \end{array}$$

(26)

where $T_0$ denotes the initial temperature.

Algorithm 1 The modified simulated annealing based search algorithm

Input:$\beta$, $\theta$, $j=0$, $k=0$, $\delta$, J and K Output:the optimal ${\alpha }^{*}$, ${\beta }^{*}$, ${\theta }^{*}$ 1Derive the value of ${\alpha }^{\left(j\right)}$ by exhaustive search with β and θ; 2Initialize value $x=x^{\left(0\right)}(x^{\left(0\right)}\in \mathrm{\Omega })$; 3Select an alternative point $z^{\left(k\right)}$ at random from Ω; 4Design a random event, and the probability of the event happening is $p(k,\tau \left(z^{\left(k\right)}\right),\tau \left(x^{\left(k\right)}\right))$; 5Perform the event. If the event occurred, let $x^{(k+1)}=z^{\left(k\right)}$; Otherwise, let $x^{(k+1)}=x^{\left(k\right)}$; 6If ${\tau }^{\left(x^{(k+1)}\right)}-{\tau }^{\left(x^{\left(k\right)}\right)}<\delta$ and $k\ge K$, $j=j+1$; 7Let $k=k+1$ and return to step 3; 8If $j\le J$, return to step 1; otherwise, the algorithm terminates;

In Algorithm 1, it is worth noting that the larger the gap between $\tau \left(z^{\left(k\right)}\right)$ and $\tau \left(x^{\left(k\right)}\right)$, the more impossible $z^{\left(k\right)}$ is selected as the next iteration point. Similarly, the smaller $T_k$ is, the less likely $z^{\left(k\right)}$ is to be selected as an iteration point. In other words, Algorithm 1 first searches the optimal solution from the overall solution space $\mathrm{\Omega }$. As k increases, the search scope will concentrate to the vicinity of the global optimal point. Meanwhile, the computational complexity of the Algorithm 1 in terms of the number of variables and the loops is $O\left(N\right)$. Deriving the $\beta$ and $\theta$, Algorithm 1 returns to step 1 and repeats the above process. Compared with the exhaustive search and the classical SABS algorithm, Algorithm 1 not only improves the search efficiency significantly, but also no longer depends on the selection of the initial point.

5. Numerical Results and Simulation

In this section, the throughput of this system is investigated. Furthermore, we compare the analytical results with simulated results for proving the analytical framework. Unless otherwise specified, the parameters in the simulation are set as follows: $P_1=P_2=1$ W, $U=$ 3 bits/s/Hz, $\eta =1$, ${\lambda }_{h_1}={\lambda }_{h_2}={\lambda }_{g_1}={\lambda }_{g_2}=1$ [12], $|{\widehat{f}}_1{|}^2=|{\widehat{f}}_2{|}^2=0.001$, ${\sigma }^2=0.001$ W.

5.1. Model Validation

Figure 4 shows the analytical throughput and the simulated throughput of the full-duplex system. We find that the analytical throughput matches well with the simulated throughput. Furthermore, we can find that the system throughput varies with different parameters, including $\alpha$, $\beta$ and $\theta$.

5.1.1. The Effect of $\alpha$

Figure 4a shows the system throughput as a function of $\alpha$ when the other two parameters are set as $\beta =0.5$ and $\theta =0.5$. For the analytical results of $P_1{\lambda }_{h_1}\ne P_2{\lambda }_{g_2}$, we set $P_1=1$ W, $P_2=2$ W. We can observe that the biggest gap between the analytical throughput and the simulated throughput is $0.5\%$ when $\alpha =0.08$. For the analytical results of $P_1{\lambda }_{h_1}=P_2{\lambda }_{g_2}$, we set $P_1=P_2=1$ W. The biggest gap between the analytical throughput and the simulated throughput is $0.8\%$ when $\alpha =0.1$. These results show that the analytical throughput reflects the practical throughput well. On the other hand, the system throughput first increases and then decreases as $\alpha$ increases. This is because that when $\alpha$ is too small, the harvested energy at R is not enough, so the transmitted energy for R is also not enough. When $\alpha$ is close to 1, the harvested energy at R is large. However, the time for information transmission is reduced, which leads to a decrease in the system throughput, too.

5.1.2. The Effect of $\beta$

Figure 4b shows the system throughput as a function of $\beta$ when the other two parameters are set as $\alpha =0.5$ and $\theta =0.5$. The analytical results are also close to the simulated throughput. When $\beta$ is very small, the energy for the link S-D is small, and this will make the transmission from S to D fail, and thus, the throughput will become small; otherwise, when $\beta$ is very large, the energy for the link D-S is small, and this will make the transmission from D to S fail. Thus, the throughput will also become small.

5.1.3. The Effect of $\theta$

Figure 4c shows the system throughput as a function of $\theta$ when the other two parameters are set as $\alpha =0.5$ and $\beta =0.5$. We could observe that the system throughput first increases with $\theta$ rapidly, and then decreases with $\theta$ rapidly. This is because from (23), $\theta$ mainly works in $min(\theta ,1-\theta )$, which makes the curve very steep.

5.2. The Non-Outage Probability

In fact, the non-outage probability is the key to obtain the system throughput. Figure 5 shows that the non-outage probability is related with SNR, and its value ranges from 0 to 1. The effect of noise power variation is equivalent to the effect of the transmission power variation, and it could be observed that the non-outage probability increases with the increases in SNR, i.e., the larger transmission power brings a larger non-outage probability. Meanwhile, Figure 5 shows that there is little gap between optimal non-outage probability and non-outage probability with fixed system parameters. Thus, we cannot optimize the non-outage probability alone, but should further consider optimizing the system throughput.

5.3. The Optimal Throughput

We further investigate the optimal throughput of the full-duplex system. From Figure 6a, we find that the system throughput of fixed parameters is smaller than the practical optimal throughput. In Figure 6b, we give the comparison of practical throughput with the throughput derived by different methods, and we set the search step length $ϵ=0.02$ for the exhaustive search method. We can observe that the throughput of these two methods is relatively accurate when SNR varies from 10 dB to 40 dB. The throughput derived by the exhaustive search method or Algorithm 1 has the biggest gap between the practical throughput when SNR is 10 dB.

6. Conclusions

In this paper, a joint time and power allocation scheme for the FD two-way EH-capable relay system from the perspective of system throughput is investigated. Considering the system using DF policy to transmit information, the SINRs at S, R and D are first derived. Thus, the non-outage probability of the system can be obtained by these SINRs. Furthermore, a new analytical method is used to obtain the non-outage probability so that the corresponding system throughput can be derived. We find that the system throughput contained three variables, i.e., time switching ratio, power allocation ratio, and the length of communication time slot, and is in the integral form. To achieve the optimal system throughput, these three variables need to be optimized simultaneously. Since the traditional exhaustive search method has high time complexity for solving the problem with multiple variables, we propose a modified SABS algorithm that could efficiently find the close-to-optimal time and power allocation to maximize throughput. Numerical results are provided to demonstrate the accuracy and the effectiveness of our new analytical framework, and they show that the SABS algorithm can derive accurate throughput results.