Energy-Efficient Power Allocation in Non-Linear Energy Harvesting Multiple Relay Systems

Abstract

In this paper, to maximize the energy efficiency (EE) in the two-hop multi-relay cooperative decoding and forwarding (DF) system for simultaneous wireless information and power transmission (SWIPT), an optimal power allocation algorithm is proposed, in which the relay energy harvesting (EH) adopts a nonlinear model. Under the constraints, including energy causality, the minimum transmission quality of information and the total transmission power at the relays, an optimization problem is constructed to jointly optimize the transmit power and power-splitting (PS) ratios of multiple relays. Although this problem is a nonlinear fractional programming problem, an iterative algorithm is developed to obtain the optimal power allocation. In particular, the joint power allocation at multiple relays is first decoupled into a single relay power allocation, and then single-relay power allocation is performed by the Dinkelbach iteration algorithm, which can be proven that it is a convex programming problem. Its closed form solutions for different polylines of EH models are obtained by using mathematical methods, such as monotonicity, Lagrange multipliers, the KKT condition and the Cardan formula. The simulation results show the superiority of the power allocation algorithm proposed in this paper in terms of EE.

1. Introduction

According to statistics, the average annual energy consumption of the information and communications industry accounts for 3% of the global total energy consumption, and carbon emissions exceed 2% of the global total [1,2]. Collaborative relay technology and energy harvesting technology have emerged to meet the quality of communication while reducing energy consumption. In the post-5G era, because wireless networks are becoming more and more complex, how to design appropriate energy allocation strategies to improve system performance is a problem that needs further research [3].

Collaborative relay technology based on the energy harvesting (EH) model has been studied. In [4], the throughput maximization problem was solved via a successive convex approximation approach, in which the relay harvests energy in the natural environment. In [5], the authors analyzed the difference between an EH relay and a normal collaborative relay. The authors in [6] presented a power allocation algorithm in which the relay could switch between amplify and forward (AF) and decoding and forwarding (DF) modes to obtain the optimal transmit power and operational mode of the relay in order to maximize the energy efficiency (EE). The power allocation problem for a two-hop two-way single-relay system to maximize the EE was studied in [7]. Both source nodes and relay nodes use EH technology, and it was the aim of the authors to maximize the system capacity in [8]. The authors of [9] mainly analyzed and compared the advantages and disadvantages of the two schemes of simultaneous wireless information and power transmission (SWIPT). However, the above EH model assumes a linear EH model, in which the power conversion efficiency from the radio frequency (RF) to the direct current (DC) is a fixed constant. In fact, because the diodes, inductors and capacitors in the actual energy harvester [10] are nonlinear, the power conversion efficiency from the RF to the DC will vary with the input power of the energy harvester circuit.

There are also preliminary studies on the application of nonlinear EH models in relay systems. In [11], a power allocation algorithm for a nonlinear EH model that maximized the total harvesting power of the energy harvester was presented. The problem of maximizing the EE of a two-way DF single-relay network was studied in [12]. The simulation results in [13] show that the power allocation based on the nonlinear EH model has higher throughput than that based on the existing linear model. In [14], the authors investigated the interrupt performance of nonlinear EH relay networks with a PS scheme. However, the above works only studied power allocation for single-relay systems and not for multi-relay systems.

This paper presents a power allocation algorithm based on the nonlinear EH model for maximizing the EE in a two-hop multi-relay cooperative forwarding system. For clarity, we have listed the contributions as follows:

• In this paper, a piecewise linear EH model is used to approximate the nonlinear EH model. Based on that, the EE maximization problem with quality of service (QoS) guaranteed is a complex, nondeterministic polynomial (NP) problem and can be decoupled into a single-relay power allocation problem. It is optimized based on the Dinkelbach iteration algorithm, and the optimization problem after iteration simplification is proven to be a convex programming problem.
• Because the optimization problem is a maximum–minimum function problem, it is decomposed into two maximization problems to obtain the optimal solution by adding the corresponding constraints. For each maximization problem, simplified expressions are obtained according to the different polylines of EH models, whose closed-form solutions of the optimal relay transmit power and optimal power-splitting (PS) ratio are obtained by using mathematical methods such as monotonicity, Lagrange multipliers, the KKT condition and the Cardan formula.
• The simulation results show that the piecewise linear EH model performs better than the traditional linear EH model, and the EE of the multi-relay system is better than that of the single-relay system.

The rest of this paper is organized as follows. In Section 2, the multi-relay cooperative DF system based on nonlinear EH is described, and the optimization problem is formulated to maximize the EE. In Section 3, the joint optimization algorithm of the relay transmit power and PS ratio is proposed. In Section 4, the algorithm is simulated, and its performance is analyzed. Finally, this paper is summarized in Section 5.

2. System Model and Optimization Problem

2.1. System Model

As shown in Figure 1, we considered the two-hop multi-relay cooperative forwarding system, including a source node $S$, a destination node $D$ and $N$ relay nodes $R_i$ for $\mathrm{i}=1,2,\dots ,N$. All the $R_i$ can harvest energy, while $S$ and $D$ cannot harvest energy. Assuming that there is no direct link between $S$ and $D$, the corresponding channel of each relay is independent of each other, where $h_i$ represents the channel gain from $S$ to $R_i$ and $g_i$ represents the channel gain from $R_i$ to $D$.

Compared with PS, time switching (TS) cannot obtain a better trade-off between the transmission rate and the RF energy, and it was proven that the PS scheme achieves higher throughput than the TS scheme in [11]. Therefore, we assumed that $R_i$ could only harvest energy from the signal sent by $S$, and each transmission block period was divided into two time slots in the PS scheme of [13]. In the first time slot, $R_i$ receives the signal sent by $S$ and divides the received signal power into two parts, according to the ratio of $1-{\rho }_i:{\rho }_i$. One part enters the information processing module to decode the information, expressed as

$$y_{{}_{R_i}}^{ID}(k)=\sqrt{1-{\rho }_i}\cdot \left(\sqrt{P_S{h}_i}\cdot s(k)+n_i^a(k)\right)+n_i^{ID}(k)$$

(1)

where $\mathrm{E}\left[{\left|s\left(k\right)\right|}^2\right]=1$; $P_S$ is the transmit power of $S$, ${\rho }_i$ is the PS ratio of $R_i$, $n_i^a(k)\sim CN\left(0,{\sigma }_{a_i}^2\right)$ is the additive white Gaussian noise of the relay receiving antenna and $n_i^{ID}(k)\sim CN\left(0,{\sigma }_{I{D}_i}^2\right)$ is the additive baseband Gaussian noise. Therefore, the SNR of the receiver of $R_i$ can be expressed as ${\gamma }_r^i=\frac{\left(1-{\rho }_i\right)\cdot P_S{h}_i}{\left(1-{\rho }_i\right)\cdot {\sigma }_{a_i}^2+{\sigma }_{I{D}_i}^2}$.

The other part enters the EH module to provide energy for the information decoding of this time slot and the information forwarding of the next time slot. The received RF signal at $R_i$ used for energy harvesting is

$$y_{{}_{R_i}}^{EH}(k)=\sqrt{{\rho }_i}\cdot \left(\sqrt{P_S{h}_i}\cdot s(k)+n_i^a(k)\right)$$

(2)

Because of $P_S{h}_i\gg {\sigma }_{a_i}^2$, the power used by $R_i$ for EH is $P_{RF}^i\left({\rho }_i\right)={\rho }_i{P}_S{h}_i$.

In many cases, because of the nonlinearity of the EH model, the problems studied are often too complex to be solved, so the piecewise linear model can not only keep the nonlinearity of the original model to a maximum extent but also solve the problem well [15]. For EH, we used a piecewise linear model, which can be expressed as

$$\begin{array}{cc}\hfill P_H^i\left({\rho }_i\right)& =\{\begin{array}{l}0,P_{RF}^i\left({\rho }_i\right)<P_{th}^1\\ a_l\cdot P_{RF}^i\left({\rho }_i\right)+b_l,P_{RF}^i\left({\rho }_i\right)\in \left[P_{th}^l,P_{th}^{l+1}\right],l=1,2,\dots ,L-1\\ P_m,P_{RF}^i\left({\rho }_i\right)>P_{th}^L\end{array}\hfill \\ & =a_l\cdot {\rho }_i{P}_S{h}_i+b_l,P_{RF}^i\left({\rho }_i\right)\in \left[P_{th}^l,P_{th}^{l+1}\right],l=0,1,2,\dots ,L\hfill \end{array}$$

(3)

where $P_{RF}^i$ is the receiving power of the EH circuit, $P_{th}=\left\{P_{th}^l|0\le l\le L+1\right\}$ is the dividing point separating $P_{RF}^i$ into $L+1$ segments, satisfying $P_{th}^0=0$, $P_{th}^{L+1}=+\infty$, $a_l$ and $b_l$ are the slope and intercept of the linear function in the $l-th$$\left(0\le l\le L\right)$ segment, respectively, satisfying $a_0=b_0=a_L=0$,$b_L=P_m$, and $P_m$ denotes the maximum harvesting power.

When it has enough segments, it can be approximated to the nonlinear EH model, but this will also increase its complexity. Therefore, the real measurement data [16] of the energy harvester can be fitted by the piecewise linear function to replace the nonlinear EH model [10]. Figure 2 compares the harvested power in different EH models.

Because the transmission performance of DF is better than AF in most instances, the relay uses DF in the second time slot, which is to say that $R_i$ uses the harvested energy to forward the decoded information to $D$. Assuming that $R_i$ sends the signal $x_i\left(k\right)$ to $D$, the received RF signal at $D$ is

$$y_D^i\left(k\right)=\sqrt{P_i{g}_i}\cdot x_i\left(k\right)+n_D\left(k\right)$$

(4)

where $\mathrm{E}\left[{\left|x_i\left(k\right)\right|}^2\right]=1$, $P_i$ is the transmit power of $R_i$ and $n_D\left(k\right)\sim CN\left(0,{\sigma }_D^2\right)$ is the additive white Gaussian noise caused by the antenna reception and radio frequency to baseband conversion at the destination node. Therefore, the signal-to-noise ratio (SNR) of the receiver of $D$ can be expressed as ${\gamma }_d^i=\frac{P_i{g}_i}{{\sigma }_D^2}$.

According to the maximum merge ratio criterion [17], the total SNR of all the links participating in the communication is $\gamma ={\displaystyle \sum _{i=1}^N{\gamma }_i}={\displaystyle \sum _{i=1}^N\min \left({\gamma }_r^i,{\gamma }_d^i\right)}$. Therefore, by using the Shannon formula, the total information transmission rate of $S-R-D$ is

$$R_{SD}=\frac{1}{2}{\log }_2\left(1+\gamma \right)=\frac{1}{2}{\log }_2\left(1+{\displaystyle \sum _{i=1}^N{\gamma }_i}\right)$$

(5)

In addition, the total power consumption during the transmission of $S-R-D$ is

$$P=P_S+2P_C+{\displaystyle \sum _{i=1}^N\left(P_i+P_Q-P_H^i\left({\rho }_i\right)\right)}$$

(6)

where $P_C$ is the consumption of power of the $S$ and $D$ decoding circuits and $P_Q$ is the consumption of power of the $R_i$ circuit.

2.2. Optimization Problem

EE is defined as the ratio of the total information transmission rate to the total power consumption for all links, expressed as

$$\eta =\frac{R_{SD}}{P}=\frac{\frac{1}{2}{\log }_2\left(1+{\displaystyle \sum _{i=1}^N{\gamma }_i}\right)}{P_S+2P_C+{\displaystyle \sum _{i=1}^N\left(P_i+P_Q-P_H^i\left({\rho }_i\right)\right)}}$$

(7)

This paper studies how to jointly optimize the relay transmit power and the PS ratio to maximize the EE. The optimization problem can be expressed as

$$\begin{array}{cc}\hfill \left(\mathrm{P}1\right)\hspace{1em}\underset{P_i,{\rho }_i}{\max }\hspace{1em}\eta & \\ \hfill s.t.C1& :0\le {\rho }_i\le 1\hfill \\ \hfill C2& :0\le P_i+P_Q\le P_H^i\left({\rho }_i\right)\hfill \\ \hfill C3& :{\gamma }_i\ge {\gamma }_{\min }\hfill \\ \hfill C4& :P_{th}^l\le {\rho }_i{P}_S{h}_i\le P_{th}^{l+1},l=0,1,2,\dots ,L\hfill \\ \hfill C5& :{\displaystyle \sum _{i=1}^N{P}_i\le P_{\max }}\hfill \end{array}$$

(8)

where ${\gamma }_{\min }=2^{R_{\min }}-1$ and $R_{\min }$ is the minimum required rate for $S-R-D$, $P_{\max }$ is the threshold of the relay nodes’ total transmit power, constraint $C2$ means the energy causality, $C3$ is the constraint of the quality of information transmission, $C4$ is the constraint that the energy harvester works in the $l-th$ linear regions and $C5$ is the constraint of the total transmit power of all relays, and it is a requirement to consider interference across the network [18].

3. Joint Optimization Algorithm for the Relay Transmit Power and PS Ratio

The optimization problem (P1) is a multi-variable joint optimization problem and a complex NP problem. Because of the piecewise linear EH model, the optimization problem has different results for its different polyline segments.

In order to solve (P1), there are two main steps as follows. In the first step, the power allocation problem of multiple relays is decoupled into a single-relay power allocation problem by using the fixed-point iteration algorithm, which assumes that only the power of this relay is to be allocated and that the power allocation of other relays is known. For convenience, some iteration variables ${\psi }_{\mathrm{i}}={\displaystyle \sum _{w=1,w\ne i}^N{\gamma }_w}$, ${\phi }_i={\displaystyle \sum _{w=1,w\ne i}^N{P}_w}$ and ${\varphi }_i={\displaystyle \sum _{w=1,w\ne i}^N\left(P_w+P_Q-P_H^w\left({\rho }_w\right)\right)}$ are introduced. In the second step, we assume ${\rho }_i=1$ and calculate the maximum number of segments which $P_{RF}^i\left({\rho }_i\right)$ may belong to, expressed as $S_i$$\left(s_i\in \left\{0,1,\dots ,S_i\right\},i=1,2,\dots ,N\right)$. In a given number of segments $s_i$, we transform this non-convex problem into a convex problem by subtraction and find a closed-form solution.

According to the Dinkelbach method [19], the maximum EE $x^{*}=\underset{P_i,{\rho }_i}{\max }\eta =\frac{\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+{\gamma }_i\left(P_i^{*},{\rho }_i^{*}\right)\right)}{P_S+2P_C+{\varphi }_i+P_i^{*}+P_Q-P_H^i\left({\rho }_i^{*}\right)}$ is achieved if and only if the following equation is satisfied:

$$\underset{P_i,{\rho }_i}{\max }\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+{\gamma }_i\left(P_i^{},{\rho }_i^{}\right)\right)-x^{*}\cdot \left(P_S+2P_C+{\varphi }_i+P_i^{}+P_Q-P_H^i\left({\rho }_i^{}\right)\right)=0$$

(9)

Therefore, we convert the fractional objective function to subtraction for a given segment $s_i$, and (P1) can be rewritten as

$$\begin{array}{cc}\hfill \left(\mathrm{P}2\right)& \underset{P_i,{\rho }_i}{\max }\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+{\gamma }_i\left(P_i^{},{\rho }_i^{}\right)\right)-x^{*}\cdot \left(P_S+2P_C+{\varphi }_i+P_i^{}+P_Q-P_H^i\left({\rho }_i^{}\right)\right)\hfill \\ & \hspace{1em}s.t.C1-C5\hfill \end{array}$$

(10)

However, (P2) is still a non-convex problem. In order to convert this problem into a convex one, we let $t_i=1-{\rho }_i$ and introduce an auxiliary variable $r_i$, expressed as

$$r_i=\frac{1}{2}{\log }_2\left(1+{\displaystyle \sum _{i=1}^N{\gamma }_i}\right)=\min \left(\frac{1}{2}{\log }_2(1+{\psi }_{\mathrm{i}}+{\gamma }_r^i),\frac{1}{2}{\log }_2(1+{\psi }_{\mathrm{i}}+{\gamma }_d^i)\right)$$

(11)

The optimization problem (P2) can be rewritten as

$$\begin{array}{cc}\hfill \left(\mathrm{P}3\right)\hspace{1em}\underset{P_i,t_i}{\max }\hspace{1em}& f(P_i,t_i,r_i)=r_i-x^{*}\cdot \left(P_S+2P_C+{\varphi }_i+P_i+P_Q-a_{s_i}\cdot (1-t_i)P_S{h}_i-b_{s_i}\right)\hfill \\ \hfill s.t.C1:& t_{\min }\le t_i\le t_{\max }\hfill \\ \hfill C2:& 0\le P_i\le a_{s_i}\cdot (1-t_i)P_S{h}_i+b_{s_i}-P_Q\hfill \\ \hfill C3:& r_i\ge \frac{1}{2}{\log }_2\left(1+{\gamma }_{\min }\right)\hfill \\ \hfill C4:& \frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{t_i\cdot P_S{h}_i}{t_i\cdot {\sigma }_{a_i}^2+{\sigma }_{I{D}_i}^2}\right)\ge r_i\hfill \\ \hfill C5:& \frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_i{g}_i}{{\sigma }_D^2}\right)\ge r_i\hfill \\ \hfill C6:& P_i\le P_{\max }-{\phi }_i\hfill \end{array}$$

(12)

where $t_{\min }=1-\frac{P_{th}^{s_i+1}}{P_S{h}_i}$ and $t_{\max }=1-\frac{P_{th}^{s_i}}{P_S{h}_i}$.

Theorem 1.

The optimization problem (P3) is a convex programming problem.

Proof of Theorem 1.

The objective function of (P3) regarding $P_i$, $t_i$ and $r_i$ is a linear function, and constraints $C1$, $C2$, $C3$ and $C6$ are linear constraints, while constraints $C4$ and $C5$ are nonlinear constraints, so (P2) is a nonlinear programming problem. Because linear functions are both convex and concave, the objective function $f(P_i,t_i,r_i)$ is convex, and constraints $C1$, $C2$, $C3$ and $C6$ are concave constraints. By letting $f_1\left(t_i\right)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{t_i\cdot P_S{h}_i}{t_i\cdot {\sigma }_{a_i}^2+{\sigma }_{I{D}_i}^2}\right)$ and $f_2\left(P_i\right)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_i{g}_i}{{\sigma }_D^2}\right)$, then we have $\frac{{\partial }^2{f}_1\left(t_i\right)}{\partial t_i{}^2}=\frac{-P_S{h}_i{\sigma }_{I{D}_i}^2\cdot \left(2a_t\cdot t_i+b_t\right)}{{\left(a_t\cdot t_i^2+b_t\cdot t_i+c_t\right)}^2\cdot 2\ln (2)}<0$,$\frac{{\partial }^2{f}_2\left(P_i\right)}{\partial P_i{}^2}=\frac{-g_i^2}{{\left(P_i{g}_i+\left(1+{\psi }_{\mathrm{i}}\right)\cdot {\sigma }_D^2\right)}^2\cdot 2\ln (2)}<0$, where $a_t=\left(P_S{h}_i+\left(1+{\psi }_{\mathrm{i}}\right)\cdot {\sigma }_{a_i}^2\right)\cdot {\sigma }_{a_i}^2$, $b_t=\left(P_S{h}_i+2\left(1+{\psi }_{\mathrm{i}}\right){\sigma }_{a_i}^2\right)\cdot {\sigma }_{I{D}_i}^2$ and $c_t=\left(1+{\psi }_{\mathrm{i}}\right)\cdot {\sigma }_{I{D}_i}^4$. Therefore, $f_1\left(t_i\right)$ is a concave function for $t_i$, and $f_2\left(P_i\right)$ is also a concave function for $P_i$ ( i.e., the constraints $C4$ and $C5$ are concave constraints). From the definition of convex programming, the optimization problem (P3) is a convex programming problem. □

As Equation (11) is composed of two formulas, (P3) is a maximum–minimum function problem and can be decomposed into two maximization problems to obtain the optimal solution by adding the corresponding constraints, for generality, we assume that all noise variances are equal (i.e., ${\sigma }_{a_i}^2={\sigma }_{I{D}_i}^2={\sigma }_D^2={\sigma }^2$).

1. If ${\gamma }_r^i\ge {\gamma }_d^i$ (i.e., $P_i\le \frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}$), we have $r_i=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+{\gamma }_d^i\right)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_{\mathrm{i}}{g}_i}{{\sigma }^2}\right)$. The optimization problem can be expressed as

$$\begin{array}{cc}\hfill \left(\mathrm{P}4\right)\hspace{1em}\underset{P_i,t_i}{\max }& y_1(P_i,t_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_{\mathrm{i}}{g}_i}{{\sigma }^2}\right)\hfill \\ & \hspace{1em}\hspace{1em}\hspace{1em}-x^{*}\cdot \left(P_S+2P_C+{\varphi }_i+P_i+P_Q-a_{s_i}\cdot (1-t_i)P_S{h}_i-b_{s_i}\right)\hfill \\ \hfill s.t.C1:& t_{\min }\le t_i\le t_{\max }\hfill \\ \hfill C2:& 0\le P_i\le a_{s_i}\cdot (1-t_i)P_S{h}_i+b_{s_i}-P_Q\hfill \\ \hfill C3:& \frac{P_{\mathrm{i}}{g}_i}{{\sigma }^2}\ge {\gamma }_{\min }-{\psi }_i\hfill \\ \hfill C4:& P_i\le \frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}\hfill \\ \hfill C5:& P_i\le P_{\max }-{\phi }_i\hfill \end{array}$$

(13)

2. If ${\gamma }_r^i\le {\gamma }_d^i$ (i.e., $P_i\ge \frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}$), we have $r_i=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+{\gamma }_r^i\right)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{t_i{P}_S{h}_i}{\left(t_i+1\right){\sigma }^2}\right)$. The optimization problem can be expressed as

$$\begin{array}{cc}\hfill \left(\mathrm{P}5\right)\hspace{1em}\underset{P_i,t_i}{\max }& y_2(P_i,t_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{t_i{P}_S{h}_i}{\left(t_i+1\right){\sigma }^2}\right)\hfill \\ & \hspace{1em}\hspace{1em}\hspace{1em}-x^{*}\cdot \left(P_S+2P_C+{\varphi }_i+P_i+P_Q-a_{s_i}\cdot (1-t_i)P_S{h}_i-b_{s_i}\right)\hfill \\ \hfill s.t.C1:& t_{\min }\le t_i\le t_{\max }\hfill \\ \hfill C2:& 0\le P_i\le a_{s_i}\cdot (1-t_i)P_S{h}_i+b_{s_i}-P_Q\hfill \\ \hfill C3:& \frac{t_i{P}_S{h}_i}{\left(t_i+1\right){\sigma }^2}\ge {\gamma }_{\min }-{\psi }_i\hfill \\ \hfill C4:& P_i\ge \frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}\hfill \\ \hfill C5:& P_i\le P_{\max }-{\phi }_i\hfill \end{array}$$

(14)

Due to the uncertainty of $P_{RF}^i$, according to the EH model in Equation (3), (P4) and (P5) are divided into three cases: Case 1, $P_{RF}^i\left(t_i\right)<P_{th}^1$; Case 2, $P_{RF}^i\left(t_i\right)>P_{th}^L$; and Case 3, $P_{RF}^i\left(t_i\right)\in \left[P_{th}^l,P_{th}^{l+1}\right]$,$l=1,2,\dots ,L-1$.

3.1. Case 1

When $P_{RF}^i\left(t_i\right)<P_{th}^1$ (i.e., $s_i=0$), we have $P_H^i=0$ and $P_i=0$, leading to a communication outage between $R_i$ and $D$. As a result, the EE is zero.

3.2. Case 2

When $P_{RF}^i\left(t_i\right)>P_{th}^L$ (i.e., $s_i=L$), we have $P_H^i\left(t_i\right)=P_m$, $t_{\min }=\min \left(0,1-\frac{P_{th}^{L+1}}{P_S{h}_i}\right)=0$ and $t_{\max }=1-\frac{P_{th}^L}{P_S{h}_i}$.

Considering the constraints of (P4) and (P5), the optimization problem can be rewritten as three simplification scenarios:

• Scenario I: When $\frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}\ge P2$, we have $t_i\ge t2$, and (P4) can be rewritten as

  $$\begin{array}{c}\underset{P_i,t_i}{\max }\hspace{1em}{y}_1(P_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_{\mathrm{i}}{g}_i}{{\sigma }^2}\right)-x^{*}\cdot \left(P_i+P_Q-P_{\mathrm{m}}\right)\hfill \\ s.t.P1\le P_i\le P2\hfill \\ \hspace{1em}{t}_{\min }\le t_i\le t_{\max }\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}{t}_i\ge t2\hfill \end{array}$$

  (15)

There is a solution to Equation (15) if and only if $t2\le t_{\max }$ is satisfied. The objective function is independent of $t_i$, so we usually let $t_i^{*}=t_{\max }$. Due to $\frac{{\partial }^2{y}_1}{\partial P_i{}^2}<0$, $y_1\left(P_i\right)$ is a convex function with respect to $P_i$, and it can be solved by derivation. Let $\frac{d{y}_1}{d{P}_i}=0$ and have $P0=\frac{1}{x^{*}\cdot 2\ln 2}-\frac{\left(1+{\psi }_i\right)\cdot {\sigma }^2}{g_i}$. The solution I can be expressed as

$$\{\begin{array}{l}{P}_i^{I*}={\left[P0\right]}_{P1}^{P2}\\ {\rho }_i^{I*}=1-t_i^{I*}=1-t_{\max }\end{array}$$

(16)

where $t2=\frac{P2\cdot g_i}{P_S{h}_i-P2\cdot g_i}$, $P1=\max (0,\frac{\left({\gamma }_{\min }-{\psi }_i\right)\cdot {\sigma }^2}{g_i})$, $P2=\min \left(P_{\max }-{\phi }_i,P_m-P_Q\right)$ and ${\left[x\right]}_a^b$ denotes $\max (a,\min (x,b))$.

• Scenario II: When $\frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}\le P2$, we have $t_i\le t_2$, and (P4) can be rewritten as

  $$\begin{array}{c}\hspace{1em}\underset{P_i,t_i}{\max }\hspace{1em}{y}_1(P_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_{\mathrm{i}}{g}_i}{{\sigma }^2}\right)-x^{*}\cdot \left(P_i+P_Q-P_m\right)\hfill \\ s.t.P1\le P_i\le \frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}\hfill \\ \hspace{1em}\hspace{1em}{t}_{\min }\le t_i\le t_{\max }\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}{t}_i\le t2\hfill \end{array}$$

  (17)

There is a solution to Equation (17) if and only if $t_{\min }\le t2$ is satisfied. The objective function is related to $t_i$, so we use the Lagrange multiplier method to convert a constrained problem to an unconstrained problem $\underset{P_i,t_i}{\max }L\left(P_i,t_i\right)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_{\mathrm{i}}{g}_i}{{\sigma }^2}\right)-x^{*}\cdot \left(P_i+P_Q-P_m\right)+{\lambda }_i\left(\frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}-P_i\right)$, where ${\lambda }_i$ is the Lagrange multiplier. Since $\frac{{\partial }^{2}L}{\partial P_i{}^2}<0$ and $\frac{{\partial }^{2}L}{\partial t_i{}^2}<0$, we use the KKT condition for the solution. Solution II can be expressed as

$$\{\begin{array}{l}{P}_i^{II*}={\left[P0\right]}_{P1}^{\frac{t_i^{II*}{P}_S{h}_i}{\left(t_i^{II*}+1\right)g_i}}\\ {\rho }_i^{II*}=1-t_i^{II*}=1-\min (t2,t_{\max })\end{array}$$

(18)

where $t2=\frac{P2\cdot g_i}{P_S{h}_i-P2\cdot g_i}$, $P0=\frac{1}{x^{*}\cdot 2\ln 2}-\frac{\left(1+{\psi }_i\right)\cdot {\sigma }^2}{g_i}$ and $P1=\max (0,\frac{\left({\gamma }_{\min }-{\psi }_i\right)\cdot {\sigma }^2}{g_i})$.

• Scenario III: Considering the constraints of (P5), the optimization problem can be rewritten as

  $$\begin{array}{c}\hspace{1em}\underset{P_i,t_i}{\max }\hspace{1em}{y}_2(P_i,t_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{t_i{P}_S{h}_i}{\left(t_i+1\right){\sigma }^2}\right)-x^{*}\cdot \left(P_i+P_Q-P_{\mathrm{m}}\right)\hfill \\ s.t.\frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}\le P_i\le P2\hfill \\ \hspace{1em}\hspace{1em}{t}_{\min }\le t_i\le t_{\max }\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}t1\le t_i\le t2\hfill \end{array}$$

  (19)

There is a solution to Equation (19) if and only if $D=\left[t1,t2\right]\cap \left[t_{\min },t_{\max }\right]\ne \varnothing$ is satisfied and we let $D=\left[t_a,t_b\right]$. First, because Equation (19) is a linear programming problem for $P_i$ and $\frac{\partial y_2}{\partial P_i}<0$, the optimal solution $P_i^{*}=\frac{t_i^{*}{P}_S{h}_i}{\left(t_i^{*}+1\right)g_i}$ can be obtained according to monotonicity. Then, the binary objective function is transformed into a one-variable objective function $y_2(t_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{t_i{P}_S{h}_i}{\left(t_i+1\right){\sigma }^2}\right)-x^{*}\cdot \left(\frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}+P_Q-P_{\mathrm{m}}\right)$, which is only related to $t_i$ by the substitution of $P_i^{*}$. We let $\frac{d{y}_2}{d{t}_i}=0$ and have $t0=\frac{g_i-x^{*}\cdot 2\ln 2\cdot \left(1+{\psi }_i\right)\cdot {\sigma }^2}{x^{*}\cdot 2\ln 2\cdot \left(P_S{h}_i+\left(1+{\psi }_i\right)\cdot {\sigma }^2\right)-g_i}$. Because $\frac{d^2{y}_2}{d{t}_i{}^2}$ can be positive or negative, we have $t_i^{*}\in \left\{t_a,t_b,{\left[t0\right]}_{t_a}^{t_b}\right\}$, and $t_i^{*}$ is determined by the value of $y_2(t_i^{*})$. Solution III can be expressed as

$$\{\begin{array}{l}{P}_i^{III*}=\frac{t_i^{III*}{P}_S{h}_i}{\left(t_i^{III*}+1\right)g_i}\\ {\rho }_i^{III*}=1-t_i^{III*}=1-\underset{t_i^{*}\in \left\{t_a,t_b,{\left[t0\right]}_{t_a}^{t_b}\right\}}{\arg \max }\left(\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{t_i^{*}{P}_S{h}_i}{\left(t_i^{*}+1\right){\sigma }^2}\right)-x^{*}\cdot \frac{t_i^{*}{P}_S{h}_i}{\left(t_i^{*}+1\right)g_i}\right)\end{array}$$

(20)

where $t1=\frac{\left({\gamma }_{\min }-{\psi }_i\right)\cdot {\sigma }^2}{P_S{h}_i-\left({\gamma }_{\min }-{\psi }_i\right)\cdot {\sigma }^2}$, $t2=\frac{P2\cdot g_i}{P_S{h}_i-P2\cdot g_i}$ and $P2=\min \left(P_{\max }-{\phi }_i,P_m-P_Q\right)$.

3.3. Case 3

When $P_{RF}^i\left(t_i\right)\in \left[P_{th}^l,P_{th}^{l+1}\right],l=1,2,\dots ,L-1$ (i.e., $s_i=l$), we have $P_H^i\left(t_i\right)=a_l·(1-t_i)P_s{h}_i+b_l$, $t_{\min }=1-\frac{P_{th}^{l+1}}{P_S{h}_i}$ and $t_{\max }=1-\frac{P_{th}^l}{P_S{h}_i}$.

Considering the constraints of (P4) and (P5), the optimization problem can be rewritten as four simplification scenarios:

• Scenario I: When $a_l\cdot (1-t_i)P_S{h}_i+b_l-P_Q\le \frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}$ and $a_l\cdot (1-t_i)P_S{h}_i+b_l-P_Q\le P_{\max }-{\phi }_i$, (P4) can be rewritten as

  $$\begin{array}{c}\hspace{1em}\hspace{1em}\underset{P_i,t_i}{\max }\hspace{1em}{y}_1(P_i,t_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_{\mathrm{i}}{g}_i}{{\sigma }^2}\right)-x^{*}\cdot \left(P_i+P_Q-a_l\cdot (1-t_i)P_S{h}_i-b_l\right)\hfill \\ s.t.P1\le P_i\le a_l\cdot (1-t_i)P_S{h}_i+b_l-P_Q\hfill \\ \hspace{1em}{t}_{\min }\le t_i\le t_{\max }\hfill \\ \hspace{1em}{t}_i\ge \max (t2,t3)\hfill \end{array}$$

  (21)

There is a solution to Equation (21) if and only if $\max (t2,t3)\le t_{\max }$ is satisfied. First, because Equation (21) is a linear programming problem for $t_i$ and $\frac{\partial y_1}{\partial t_i}<0$, $t_i^{*}=\max (t2,t3,t_{\min })$ can be obtained according to monotonicity. Then, since $\frac{{\partial }^2{y}_1}{\partial P_i{}^2}<0$, $y_1\left(P_i\right)$ is a convex function with respect to $P_i$, and it can be solved by derivation. Solution I can be expressed as

$$\{\begin{array}{l}{P}_i^{I*}={\left[P0\right]}_{P1}^{a_l\cdot (1-t_i^{I*})P_S{h}_i+b_l-P_Q}\\ {\rho }_i^{I*}=1-t_i^{I*}=1-\max (t2,t3,t_{\min })\end{array}$$

(22)

where $t2=\frac{-b+\sqrt{b^2-4ac}}{2a}$ is the positive root of the quadratic equation $a\cdot t_i^2+b\cdot t_i+c=0$, $a=a_l{g}_i{P}_S{h}_i>0$, $b=P_S{h}_i-\left(b_l-P_Q\right)g_i$, $c=-(a_l{P}_S{h}_i+b_l-P_Q)g_i$, $t3=\frac{a_l{P}_S{h}_i+b_l-P_Q-\left(P_{\max }-{\phi }_i\right)}{a_l{P}_S{h}_i}$, $P0=\frac{1}{x^{*}\cdot 2\ln 2}-\frac{\left(1+{\psi }_i\right)\cdot {\sigma }^2}{g_i}$ and $P1=\max (0,\frac{\left({\gamma }_{\min }-{\psi }_i\right)\cdot {\sigma }^2}{g_i})$.

• Scenario II: When $P_{\max }-{\phi }_i\le a_l\cdot (1-t_i)P_S{h}_i+b_l-P_Q$ and $P_{\max }-{\phi }_i\le \frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}$, (P4) can be rewritten as

  $$\begin{array}{c}\hspace{1em}\hspace{1em}\underset{P_i,t_i}{\max }\hspace{1em}{y}_1(P_i,t_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_{\mathrm{i}}{g}_i}{{\sigma }^2}\right)-x^{*}\cdot \left(P_i+P_Q-a_l\cdot (1-t_i)P_S{h}_i-b_l\right)\hfill \\ s.t.P1\le P_i\le P_{\max }-{\phi }_i\hfill \\ \hspace{1em}{t}_{\min }\le t_i\le t_{\max }\hfill \\ \hspace{1em}t4\le t_i\le t3\hfill \end{array}$$

  (23)

There is a solution to Equation (23) if and only if $D=\left[t4,t3\right]\cap \left[t_{\min },t_{\max }\right]\ne \varnothing$ is satisfied and we let $D=\left[t_a,t_b\right]$. Solved as in Equation (21), solution II can be expressed as

$$\{\begin{array}{l}{P}_i^{II*}={\left[P0\right]}_{P1}^{P_{\max }-{\phi }_i}\\ {\rho }_i^{II*}=1-t_i^{II*}=1-t_a\end{array}$$

(24)

where $P0$, $P1$ and $t3$ are the same as in Scenario I: $t4=\frac{\left(P_{\max }-{\phi }_i\right)\cdot g_i}{P_S{h}_i-\left(P_{\max }-{\phi }_i\right)\cdot g_i}$.

• Scenario III: When $\frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}\le a_l\cdot (1-t_i)P_S{h}_i+b_l-P_Q$ and $\frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}\le P_{\max }-{\phi }_i$, (P4) can be rewritten as

  $$\begin{array}{c}\hspace{1em}\underset{P_i,t_i}{\max }\hspace{1em}{y}_1(P_i,t_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_{\mathrm{i}}{g}_i}{{\sigma }^2}\right)-x^{*}\cdot \left(P_i+P_Q-a_l\cdot (1-t_i)P_S{h}_i-b_l\right)\hfill \\ s.t.P1\le P_i\le \frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}\hfill \\ \hspace{1em}{t}_{\min }\le t_i\le t_{\max }\hfill \\ \hspace{1em}{t}_i\le \min \left(t2,t4\right)\hfill \end{array}$$

  (25)

There is a solution to Equation (25) if and only if $t_{\min }\le \min \left(t2,t4\right)$ is satisfied. Using the Lagrange multiplier method to convert a constrained problem to an unconstrained problem $\underset{P_i,t_i}{\max }L\left(P_i,t_i\right)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{P_{\mathrm{i}}{g}_i}{{\sigma }^2}\right)-x^{*}\cdot \left(P_i+P_Q-a_l\cdot (1-t_i)P_S{h}_i-b_l\right)+{\lambda }_i\left(\frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}-P_i\right)$, where ${\lambda }_i$ is the Lagrange multiplier, since $\frac{{\partial }^{2}L}{\partial P_i{}^2}<0$ and $\frac{{\partial }^{2}L}{\partial t_i{}^2}<0$, the optimal solution can be obtained by the KKT condition as follows

$$\{\begin{array}{l}{P}_i^{III*}=\max (P1,\frac{t_i^{III*}{P}_S{h}_i}{\left(t_i^{III*}+1\right)g_i})\\ {\rho }_i^{III*}=1-t_i^{III*}=1-{\left[t0\right]}_{t_{\min }}^{\min \left(t2,t4,t_{\max }\right)}\end{array}$$

(26)

where $t0$ is the real root of the univariate cubic equation $a\cdot {\left(t0\right)}^3+b\cdot {\left(t0\right)}^2+c\cdot t0+d=0$, $a=x^{*}\cdot 2\ln 2\cdot a_l{g}_i\cdot \left(P_S{h}_i+\left(1+{\psi }_i\right)\cdot {\sigma }^2\right)$, $b=x^{*}\cdot 2\ln 2\cdot a_l{g}_i\cdot \left(2P_S{h}_i+3\left(1+{\psi }_i\right)\cdot {\sigma }^2\right)$, $c=x^{*}\cdot 2\ln 2\cdot \left(a_l{g}_i+1\right)\cdot \left(P_S{h}_i+\left(1+{\psi }_i\right)\cdot {\sigma }^2\right)+4x^{*}\ln 2\cdot a_l{g}_i\cdot \left(1+{\psi }_i\right)\cdot {\sigma }^2-g_i$, $d=x^{*}\cdot 2\ln 2\cdot \left(a_l{g}_i+1\right)\cdot \left(1+{\psi }_i\right)\cdot {\sigma }^2-g_i$, $t2$ and $P1$ are the same as in Scenario I and $t4$ is the same as in Scenario II.

In order to solve this univariate cubic equation, we use Girolamo Cardan’s formula. We let $t0=Y-\frac{b}{3a}$ and change it into a standard univariate cubic equation $Y^3+pY+q=0$, where $p=\frac{c}{a}-\frac{b^2}{3a^2}$ and $q=\frac{d}{a}-\frac{bc}{3a^2}+\frac{2b^3}{27a^3}$. The discriminant is $\Delta ={\left(\frac{q}{2}\right)}^2+{\left(\frac{p}{3}\right)}^3$ and can be divided into three cases: when $\Delta >0$, there is a real root $Y_1=u+v$; when $\Delta =0$, there are two real roots $Y_1=u+v$ and $Y_2=w·u+w^2$ and when $\Delta <0$, there are three real roots $Y_1=u+v$, $Y_2=w\cdot u+w^2\cdot v$ and $Y_3=w^2\cdot u+w\cdot v$, where $u=\sqrt[3]{-\frac{q}{2}-\sqrt{\Delta }}$, $v=\sqrt[3]{-\frac{q}{2}+\sqrt{\Delta }}$ and $w=\frac{-1+\sqrt{3}i}{2}$.

• Scenario IV: Considering the constraints of (P5), the optimization problem can be rewritten as

  $$\begin{array}{c}\hspace{1em}\underset{P_i,t_i}{\max }\hspace{1em}{y}_2(P_i,t_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{t_i{P}_S{h}_i}{\left(t_i+1\right){\sigma }^2}\right)-x^{*}\cdot \left(P_i+P_Q-a_l\cdot (1-t_i)P_S{h}_i-b_l\right)\hfill \\ s.t.\frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}\le P_i\le \min \left(a_l\cdot (1-t_i)P_S{h}_i+b_l-P_Q,P_{\max }-{\phi }_i\right)\hfill \\ t_{\min }\le t_i\le t_{\max }\hfill \\ t1\le t_i\le \min \left(t2,t4\right)\hfill \end{array}$$

  (27)

There is a solution to Equation (27) if and only if $D=\left[t1,\min \left(t2,t4\right)\right]\cap \left[t_{\min },t_{\max }\right]\ne \varnothing$ is satisfied and we let $D=\left[t_a,t_b\right]$. Solved as in Equation (19), solution IV can be expressed as

$$\{\begin{array}{l}{P}_i^{IV*}=\frac{t_i^{IV*}{P}_S{h}_i}{\left(t_i^{IV*}+1\right)g_i}\\ {\rho }_i^{IV*}=1-t_i^{IV*}=1-\underset{t_i^{*}\in \left\{t_a,t_b,{\left[t0\right]}_{t_a}^{t_b}\right\}}{\arg \max }{y}_2(t_i^{*})\end{array}$$

(28)

where $y_2(t_i)=\frac{1}{2}{\log }_2\left(1+{\psi }_{\mathrm{i}}+\frac{t_i{P}_S{h}_i}{\left(t_i+1\right){\sigma }^2}\right)-x^{*}\cdot \left(\frac{t_i{P}_S{h}_i}{\left(t_i+1\right)g_i}+P_Q-a_l\cdot (1-t_i)P_S{h}_i-b_l\right)$, $t1=\frac{\left({\gamma }_{\min }-{\psi }_i\right)\cdot {\sigma }^2}{P_S{h}_i-\left({\gamma }_{\min }-{\psi }_i\right)\cdot {\sigma }^2}$ and $t0$, $t2$ and $t4$ are the same as in Scenario III.

The power allocation algorithm in this paper consists of an outer loop algorithm and an inner loop algorithm. Since only single-relay power is assumed to be allocated and other relay power allocations are known in the solution, a fixed-point iteration algorithm is needed to obtain joint power optimization for all relays (i.e., an outer loop algorithm based on a fixed-point iteration, given by Algorithm 1). In the single-relay power allocation, given the initial optimal EE $x^{*}$ and the segment $s_i$ of the EH model, the Dinkelbach iteration algorithm is used to obtain the optimal relay transmit power, the optimal PS ratio and the optimal EE of $R_i$ (i.e., the inner loop algorithm based on the Dinkelbach iteration, given by Algorithm 2). The workflow of this power allocation algorithm is shown in Figure 3.

Algorithm 1. Outer loop algorithm based on a fixed-point iteration.

1: Initialize the maximum iterations$I_{\max 1}$and the maximum error toleranceδ. Set the iteration indexm = 0, the optimal relay transmission power $P_i^0$ and the optimal PS ratio ${\rho }_i^0,(i=1,...,N)$.

2:  repeat

3:   For fixed $P_m^w$ and ${\rho }_m^w$, solve (P3) by Algorithm 2 to get $P_{m+1}^w$, ${\rho }_{m+1}^w$ and ${\eta }^{*},i=1,...,N$

4:   if $\left|P_i^{m+1}-P_i^m\right|<\delta$ and $\left|{\rho }_i^{m+1}-{\rho }_i^m\right|<\delta$ then

5:      Convergence = true

6:      return $P_i^{*}=P_i^{m+1}$, ${\rho }_i^{*}={\rho }_i^{m+1}$

7:   else

8:      m = m + 1 and Convergence = false

9:   end if

10:  until Convergence = true or $m=I_{\max 1}$

Algorithm2. Inner loop algorithm based on the Dinkelbach iteration.

1:We set ${\rho }_i=1$ and Calculate maximum number of segments Si.

2:for $s_i=1:S_i$.

3:   Initialize the maximum iterations $I_{\max 2}$, the maximum error tolerance δ and the optimal energy efficiency $x^{*}$. Set iteration index k = 0.

4:   repeat

5:      ${\psi }_{\mathrm{i}}={\displaystyle \sum _{w=1,w\ne i}^N{\gamma }_w}$, ${\phi }_i={\displaystyle \sum _{w=1,w\ne i}^N{P}_w}$, ${\varphi }_i={\displaystyle \sum _{w=1,w\ne i}^N}\left(P_w+P_Q-P_H^w\left({\rho }_w\right)\right)$

6:       if si = L then

7:          Solve the power allocation by using (16), (18), (20) to get Pi, ρi.

8:       else

9:          Solve the power allocation by using (22), (24), (26), (28) to get Pi, ρi.

10:       end if

11:       Update $x_{k+1}^{*}$.

12:       if $\left|x_{k+1}^{*}-x_k^{*}\right|<\delta$ then

13:             Convergence = true

14:             return $P_i^{s_i}=P_i^{}$, ${\rho }_i^{s_i}={\rho }_i^{}$, ${\eta }_{s_i}=x_{k+1}^{*}$

15:       else

16:          k = k + 1 and Convergence = false

17:       end if

18:    until Convergence = true or $k=I_{\max 2}$

19:end for

20:$\left[{\eta }^{*},s_i\right]=\max \left({\eta }_{s_i}\right)$

21:return $P_i^{m+1}=P_i^{s_i}\left(s_i\right)$, $P_i^{m+1}=P_i^{s_i}\left(s_i\right)$, ${\eta }^{*}$

In the inner loop algorithm, in the worst case, the Dinkelbach iteration is $I_{\max 2}$, and it needs to be performed $(L+1)$ times. The outer loop algorithm contains the $I_{\max 1}$ inner loop algorithm, where the fixed-point iteration is $N$ in the worst case. Therefore, the complexity is $O\left((L+1)N{I}_{\max 1}{I}_{\max 2}\right)$.

4. Simulations

In this section, the simulation results will be analyzed, and the effectiveness of the proposed algorithm will be proven. The simulation scenario is shown in Figure 1, where the distance between $S$ and $D$ is $d=10m$ and the relay is randomly distributed in the circle with the midpoint of $S$ and $D$ as the center and radius, expressed as $r=2m$. Like [12], all the channels are Rayleigh fading channels, where the path loss exponent $\alpha$ is $\alpha =3$ and the small-scale fading is Rayleigh fading with a mean value of 0 and variance of 1. In accordance with [12,20,21], we set the simulation parameters as follows: $P_S=30 \mathrm{dBm}$, $P_C=10 \mathrm{mW}$, $P_Q=0.1 \mathrm{mW}$, ${\sigma }^2=-70 \mathrm{dBm}$ and $R_{\min }=3 \mathrm{bps/Hz}$. The parameters of the EH model in Figure 2 are shown in

Table 1. Parameters of the EH model.

Linear EH model	$P_{th}=\left[0,60\right] \mathrm{mW}$, $\mathrm{a}=[0.3957,0]$, $\mathrm{b}=[0,23.742] \mathrm{mW}$
Piecewise linear EH model (L = 4)	$P_{th}=\left[0,5,15,60\right] \mathrm{mW}$, $\mathrm{a}=[0.4,0.88,0.2877,0]$, $\mathrm{b}=[0,-2.4,6.4845,23.7465] \mathrm{mW}$
Piecewise linear EH model (L = 5)	$P_{th}=\left[0,5,15,25,60\right] \mathrm{mW}$, $\mathrm{a}=[0.4,0.88,0.896,0.1139,0]$, $\mathrm{b}=[0,-2.4,-2.64,16.9125,23.7465] \mathrm{mW}$
Piecewise linear EH model (L = 6)	$P_{th}=\left[0,5,15,25,40,60\right] \mathrm{mW}$,$\mathrm{a}=[0.4,0.88,0.896,0.2406,0.018625,0]$,$\mathrm{b}=[0,-2.4,-2.64,13.745,22.629,23.7465] \mathrm{mW}$

.

Figure 4 shows the variation of the EE in different EH models under the same relay total transmit power threshold. It can be seen that the EE increased as the number of relays increased. This is because, compared with a single relay, multiple relays can improve the information transmission rate. Since the relay uses the harvested energy for DF without consuming fixed energy, the increase in the number of relays can also improve the EE of the system. In addition, it can also be observed that the performance of the piecewise linear EH model was better than that of the linear EH model in [7]. This is because the piecewise linear EH model improves the accuracy of the traditional linear EH one. The more segments, the better the performance, but its growth rate will decrease with the increase in the segment number. Therefore, the piecewise linear EH model with $L=5$ was selected for simulation analysis in this paper.

Figure 5 shows that the EE changed with the total transmit power threshold of relays $P_{\max }$, and the number of relays was $N=4$. It can be seen that the EE increased with the increase of $P_{\max }$, and the increase gradually decreased until it no longer increased. This was because the energy harvested by the EH model was limited. When $P_{\max }$ exceeds the total power that can be used after EH, the EE will not increase with the increase of $P_{\max }$. In this simulation scenario, when $P_{\max }>40\mathrm{mW}$, the EE would not change significantly with the increase in $P_{\max }$.

Figure 6 shows that the EE varied with the different transmit powers of source $P_S$. It can be seen that the EE first increased with the increase of $P_S$ and then decreased with the increase of $P_S$. This is because when $P_S$ is very small, the energy harvested by the relay may not be enough to complete the information transmission, which makes the EE relatively low. With the increase of $P_S$, the relay can harvest more sufficient energy to complete the information transmission, so the EE correspondingly increases. However, when $P_S$ exceeds a certain value, the energy consumption is too high, and the EE will decrease with the increase of $P_S$. In this simulation scenario, when $P_S=18 \mathrm{dBm}$, the EE is the highest. In addition, compared with the single-relay in [13], multi-relay cooperation can effectively improve performance.

Figure 7 shows the change in EE with the noise power ${\sigma }^2$. The lower the noise power, the higher the energy efficiency. This is because the reduction of noise power can effectively improve the information transmission rate. Figure 8 shows the change in EE with a different path loss exponent $\alpha$. The larger the path loss exponent $\alpha$, the smaller the channel gain, the smaller the information transmission rate, the smaller the harvested energy and the smaller the energy efficiency. Figure 9 shows the change in EE with the distance $d$ between S and D. The larger the distance $d$, the smaller the information transmission rate and the smaller the energy efficiency. Figure 10 shows the change in EE with the position of relays. The closer the relay is to the source node, the higher the energy efficiency.

5. Conclusions

In this paper, a power allocation algorithm based on a nonlinear EH model is proposed for a two-hop half-duplex multiple cooperative forwarding relays system. The fixed-point iterative algorithm and Dinkelbach iterative algorithm are used to solve the optimization problem by jointly optimizing the transmit power and PS ratios of multiple relays to maximize the EE. The simulation results show that, compared with the conventional linear EH model, the designed EH model had better performance in terms of EE. Meanwhile, the superior EE performance of the multi-relay system was also verified by comparison with different numbers of relays. Compared with the linear model, when the piecewise linear model (n = 5) was used, the energy efficiency increased by 2%. Compared with a single relay, when four relays were used for cooperation, the energy efficiency increased by 26%.