Transmission Power and Antenna Allocation for Energy-Efficient RF Energy Harvesting Networks with Massive MIMO

Abstract

The optimum transmission strategy for maximizing energy efficiency (EE) of a multi-user massive multiple-input multiple-output (MIMO) system in radio frequency energy harvesting networks is investigated. We focus on dynamic time-switching (TS) antennas, to avoid the practical problems of power-splitting antennas, such as complex architectures, power loss and signal distortion when splitting the power of the received signal into power for information decoding (ID) and energy harvesting (EH). However, since a single TS antenna cannot serve ID and EH simultaneously, the MIMO system is considered in this paper. We thus formulate an EE optimization problem and propose an iterative algorithm as a tractable solution, including an antenna selection strategy to optimally switch each TS antenna between ID mode and EH mode using nonlinear fractional programming and the Lagrange dual method. Further, the problem is solved under practical constraints of maximum transmission power and outage probabilities for a minimum amount of harvested power and rate capacity for each user. Simulation results show that the proposed algorithm is more energy-efficient than that of baseline schemes, and demonstrates the trade-off between the required amount of harvested power and energy efficiency.

1. Introduction

With the growing popularity of green communication systems and energy-efficient transmission systems through energy efficiency (EE) maximization, RF energy harvesting and various optimization solutions to optimally allocate communication resources have been actively researched [1,2]. In RF energy harvesting networks (RF-EHNs), since energy and information are simultaneously transmitted by RF signals, RF energy harvesting and information reception can be theoretically performed from the same input RF signal with the same receiving antenna [3]. For this concept of simultaneous wireless information and power transfer (SWIPT) [4], we can practically consider two types of receiving antenna, power-splitting (PS) and time-switching (TS) [1]. The PS antenna, however, is known to have unsolved practical problems such as more complex architectures than that of the TS antenna [5], power loss and signal distortion when splitting the received signal power between power for information decoding (ID) and energy harvesting (EH). On the other hand, the TS antenna is easy to implement and has no power loss or signal distortion, but there is a disadvantage in that the single TS antenna cannot serve ID and EH simultaneously compared with the PS antenna. Therefore, it is necessary to employ a massive multiple-input multiple-output (MIMO) system with multiple TS antennas to simultaneously serve ID and EH, and to solve a derived antenna selection (AS) problem in order to properly switch each TS antenna to either ID mode or EH mode. There are related studies for EE optimization in RF-EHNs with PS and TS antennas, such as the resource allocation problem for a two-user EHN with the single TS antenna [6], a multi-user (MU) network with the single PS antenna [7], and MU single-input single-output (SISO) systems with the TS scheme [4] and the PS scheme [8]. In [9], a three-node wireless powered communication system was studied, in order to achieve the maximum throughput by balancing the time duration between the wireless power transfer phase and the information transfer phase in SISO systems. In [10], the authors found that the lifetime of randomly distributed wireless sensor nodes can be increased significantly as a result of Wireless Energy Harvesting (WEH) from ambient RF signals, and derived the analytical expression for the probability of successful communication between two types of sensors. In [11], the network lifetime gain for a two-way relay network in which the source, destination and RF-powered relay nodes are modeled as three independent Poisson Point Processes was investigated, and it was found that the lifetime of the network can be greatly increased by enabling RF energy harvesting and network coding. In [12], the authors analyzed the performance of ambient RF energy harvesting sensor networks where the distribution of ambient RF sources is modeled as an α-Ginibre point process that reflects the repulsion among the sources, and provided semi-closed-form expressions for the mean of the harvested energy and the exact power outage probability. In [13], a dynamic power splitting scheme was proposed in a point-to-point wireless link considering three special cases of dynamic power splitting: static power splitting, time switching, and on–off power splitting. The time switching-based simultaneous wireless information and power transfer was further investigated in a multicast system [14] and a multi-user cooperative wireless network [15]. The power splitting schemes were further investigated under different scenarios in [16,17,18,19]. In [20], a resource allocation scheme to jointly optimize transmit power, and transfer time and spatial beam so as to achieve secrecy for wireless powered communication was investigated under an information interception scenario by an eavesdropper. In [21], a MIMO wiretap channel for SWIPT was studied in order to design an optimal transmit covariance matrix for maximizing the ergodic secrecy rate. In [22,23], a MIMO system-based tradeoff between wireless energy and information transfer was investigated with a time duration divided into two times each for ID and EH. A hybrid access point broadcasting wireless energy and information to a set of distributed users was considered in [24,25], and optimal resource allocation was derived in half- and full-duplex modes. In [26,27], an RF energy harvesting MU massive MIMO network based on the TS antenna was studied under the assumption that each user had a single antenna to receive wireless energy and information.

Although the aforementioned works are variously studied, the optimal transmission power control and optimal AS strategy for optimally switching multiple TS antennas to receive wireless energy and information, while satisfying their minimum channel capacity and minimum amount of harvested energy in the field of TS antenna-based RF energy harvesting MU massive MIMO networks, have not yet been investigated.

In this paper, we focus on a system model for RF-EHNs based on the massive MIMO system using the TS antenna for SWIPT to avoid the practical issues of using the PS antenna, and overcome the problem whereby a single TS antenna cannot serve ID and EH at the same time. We study how to concurrently and efficiently transmit the data and energy while optimizing energy efficiency in RF-EHNs. The objective of the paper is to maximize “energy efficiency (bit/Joule)”, encompassing both data transmission and RF energy harvesting, in an RF-EHN context. To find the maximum EE in the TS antenna-based MU massive MIMO system, we estimate (1) how many transmission antennas the transmitter should allocate to each user; (2) how much transmission power for each transmission antenna is needed for each user; and (3) how many TS antennas for each user should be switched to the ID mode or the EH mode in the AS problem.

Our main result is that we greatly improved energy efficiency under a new resource allocation strategy with the proposed MU massive MIMO systems in RF-EHNs using nonlinear fractional programming and the Lagrange dual method. Further, this was achieved by solving the problem of optimally allocating the three parameters, the transmission and receiving antenna allocation, and transmission power allocation, which is influenced by the dynamic distances between the transmitter and receivers under practical constraints.

The rest of the paper is organized as follows. Section 2 describes the system model, and formulates the EE with its constraints. Section 3 offers solutions, and discusses EE optimization algorithms with given constraints using Lagrange dual decomposition and nonlinear fractional programming. Section 4 presents simulation results. Finally, Section 5 concludes the paper.

2. System Model and Problem Formulation

2.1. Notation

We use boldface capital and lower-case letters to denote matrices and vectors. ${ℂ}^{N\times M}$ denotes the set of all $N\times M$ matrices with complex entries. ${\mathbf{A}}^{\mathcal{H}}$ and rank$(\mathbf{A})$ represent the Hermitian transpose and rank of matrix A. vec(A) denotes the vectorization of matrix A by stacking its columns from left to right to form a column vector. ${\mathbf{I}}_N$ is the $N\times N$ identity matrix. $|·|$ and ${|·|}_p$ denote the absolute value of a complex scalar and the $l_p$-norm of a vector, respectively. The circularly symmetric complex Gaussian (CSCG) distribution is denoted by $\mathcal{C}\mathcal{N}(\mu ,{\mathsf{\sigma }}^2)$ with mean μ and variance ${\mathsf{\sigma }}^2$. ~ stands for “distributed as”. $\mathrm{E}(·)$ is the expectation operator. $\lfloor x\rfloor$ is the floor function denoting the largest integer not greater than $x$.

2.2. MU Massive MIMO RF-EHNs

In this section, we consider wireless energy harvesting and information transmission systems in an orthogonal frequency division multiplexing (OFDM) downlink network, where a set of UEs, $k\in \{1,\dots ,K\},$ is recharging electric energy from dedicated radio sources and information communications through a single hybrid-access point (H-AP), as in Figure 1. Each user k in the UE set has been equipped with TS antennas as many as $n_k^T$, where each TS antenna can be switched and worked as an antenna or rectenna for ID or EH under the proposed resource allocation algorithm. For example, if a UE has two TS antennas, it can be classified into 3 modes: (ID, ID), (ID, EH), and (EH, EH). The total number of TS antennas $n_k^T=n_k^{ID}+n_k^{EH}$ is composed of a sum of $n_k^{ID}$ and $n_k^{EH}$, where $n_k^{ID}$ and $n_k^{EH}$ are the numbers of antennas for the ID and the number of rectennas for the EH, respectively.

2.3. Multi-User MIMO Channel Model

We consider a Rayleigh fading channel with $m_k$ downlink transmission antennas allocated as a part of the total number of transmission antennas M at the H-AP and $n_k^{ID}$ antennas at the receiver for each user $k$ (${\sum }_{k=1}^K{m}_k\le M, \forall k\in K$). The $n_k^{ID}\times 1$ vector at the desired user’s receiver can be modeled as:

$$\mathit{y}{[\mathrm{t}]}_k^{ID}=\mathbf{H}{[\mathrm{t}]}_k\mathit{x}{[\mathrm{t}]}_k+\mathit{w}{[\mathrm{t}]}_k, 1\le k\le K$$

(1)

$$\mathit{y}{[\mathrm{t}]}_k^{EH}=\mathbf{G}{[\mathrm{t}]}_k\mathit{x}{[\mathrm{t}]}_k+\mathit{w}{[\mathrm{t}]}_k, 1\le k\le K$$

(2)

where $\mathbf{H}{[\mathrm{t}]}_k\in {ℂ}^{n_k^{ID}\times m_k}$ and $\mathbf{G}{[\mathrm{t}]}_k\in {ℂ}^{n_k^{EH}\times m_k}$ are matrices of independent Rayleigh fading channels for the ID and EH, respectively. $h_{ij,k}$ and $g_{ij,k}$ are entries of $\mathbf{H}{[\mathrm{t}]}_k$ and $\mathbf{G}{[\mathrm{t}]}_k$, and follow independent and identically distributed (i.i.d.) $\mathcal{C}\mathcal{N}(0,{\sigma }^2)$. Further, $h_{ij,k}=\sqrt{l_k{d}_k^{-\varphi }{f}_{ij,k}}$ is the square root of channel gain, where $d_k$, $\varphi$ and $l_{i,k}$ are respectively the distance between a user $k$ and the H-AP, the corresponding path loss exponent, and a path loss constant determined by both the carrier frequency of $w_k$ and the antenna gain. For small-scale Rayleigh fading coefficients, $f_{ij,k}~\mathrm{Exp}(l_k{d}_k^{-\varphi }/{\sigma }^2)$ are exponentially distributed random variables. In addition, we assume the CSI is not available at the transmitter, so the transmission power for each transmitter’s antenna has to be distributed equally with equal gain combining (EGC) to maximize the channel capacity. $\mathit{x}{[\mathrm{t}]}_k$ is a transmission signal vector having a vector size $m_k\times 1$ for user $k$, and $\mathit{x}{[\mathrm{t}]}_k=\mathbf{v}{[\mathrm{t}]}_k{s}_k$, where $\mathbf{v}{[\mathrm{t}]}_k$ is a loaded power vector with entries of $v_{i,k}=\sqrt{p_k}$. $s_k\in ℂ$ is the information-carrying symbol, and $\mathbb{E}[{|s_k|}^2]=1$ is assumed without loss of generality. The additive noise $\mathit{w}{[\mathrm{t}]}_k$ is assumed to be a complex Gaussian vector with zero mean and covariance matrix $N_0{w}_k{\mathbf{I}}_{n_k^{ID}}$, where $N_0$ is the noise power spectral density. Also, we assume that the system has $K$ parallel frequency channels with bandwidth $w_k$ $(1\le k\le K)$, and that the H-AP operates with the SWIPT for one scheduling time slot. In each time slot, H-AP transmits signal $\mathit{x}{[\mathrm{t}]}_k$ to user $k$ on the frequency channel, and the user receives signal $\mathit{y}{[\mathrm{t}]}_k.$

2.4. Channel Capacity of Multi-Antenna

The mutual information of a MIMO channel capacity $I_k$ for user $k$ is as follows [28,29,30]:

$$I_k=w_k{\log }_2\det \left({{\mathbf{I}}_{n_k^{ID}}+\frac{{\rho }_{i,k}}{m_k}{\mathbf{H}}_k{\mathbf{H}}_k}^{\mathscr{H}}\right),$$

(3)

where ${\rho }_{i,k}$ is the average SNR for the $i^{th}$ individual’s single-input single-output (SISO) channel, and ${\rho }_{i,k}=P_k{\left|{\overline{h}}_{i,k}\right|}^2/N_o{w}_k$, where $P_k=m_k{p}_k$ is the total transmission power [29,31], and $p_k$ is the transmission power for each transmission antenna for downlink SWIPT. Also, ${\overline{h}}_{i,k}~\mathcal{C}\mathcal{N}\left(0,\frac{{\sigma }^2}{m_k}\right)$ is the local average of channel gains in the $i^{th}$ row of ${\mathbf{H}}_k$. Note that the MIMO channel is decoupled into rank$({\mathbf{H}}_k)$ parallel SISO channels and the variance of the local average of channel gains, $\frac{{\sigma }^2}{m_k}$, is calculated by $\mathrm{Var}({\overline{h}}_{i,k})=\frac{1}{{m_k}^2}\mathrm{Var}(h_{i1,k}+h_{i2,k}+\dots +h_{i{m}_k,k})=\frac{\mathrm{Var}(h_{i1,k})+\mathrm{Var}(h_{i2,k})+\dots +\mathrm{Var}(h_{i{m}_k,k})}{{m_k}^2}=\frac{{\sigma }^2}{m_k}$, where $h_{ij,k}~\mathcal{C}\mathcal{N}(0,{\sigma }^2)$. If ${\mathbf{H}}_k$ is a full rank matrix of rank$({\mathbf{H}}_k)=\min (m_k, n_k^{ID})$, then we have a $\min (m_k, n_k^{ID})$ of non-zero singular values and independent channels. By using singular value decomposition, Equation (3) can be simplified as follows [28]:

$$I_k={\displaystyle \sum _{i=1}^{N_{min}}}{w}_k{\log }_2\left(1+\frac{{\rho }_{i,k}}{m_k}{\lambda }_{i,k}\right),$$

(4)

where $N_{min}=\min (m_k, n_k^{ID})$ is the number of independent transmitter-receiver SISO paths. ${\lambda }_i$ are unordered eigenvalues of the matrix ${{\mathbf{H}}_k{\mathbf{H}}_k}^{\mathscr{H}}$. However, the performance analysis based on the exact form of the channel capacity of (4) seems intractable and complex. As an alternative, we use an upper bound of the channel capacity following reference [32], which can be evaluated numerically. The expression of the upper bound of the channel capacity is very useful, and gives simple mathematical expression for outage constraints in Section 3.2 in this paper and optimal solution in iterative forms in Section 3.3. From the reference [32], we have an inequality by a concave function $\kappa (t)$:

$$\sum _i{\chi }_i\kappa (t_i)\le \kappa \left(\sum _i{\chi }_i{t}_i\right),$$

(5)

where ${\chi }_i\ge 0$ and ${\sum }_i{\chi }_i=1$. Clearly, equality holds when all $t_i$ are equal, or when the sum has only one term. Note that $w_k{\log }_2\left(1+\frac{t_i}{m_k}\right)$ is a concave function in $t_i$. Therefore, we set $\kappa (t_i)=w_k{\log }_2\left(1+\frac{t_i}{m_k}\right)$, $t_i={\rho }_{i,k}{\lambda }_i$, and ${\chi }_i=1/N_{min}$. Consequently,

$$\sum _{i=1}^{N_{min}}{w}_k{\log }_2\left(1+\frac{{\rho }_{i,k}{\lambda }_{i,k}}{m_k}\right)\le N_{min}{w}_k{\log }_2\left(1+\frac{{\sum }_{i=1}^{N_{min}}{\rho }_{i,k}{\lambda }_{i,k}}{m_k{N}_{min}}\right)=c_k(p_k,n_k^{ID},m_k).$$

(6)

Then, the upper bound of the overall network capacity $C(p_k,n_k^{ID},m_k)$ is:

$$C(p_k,n_k^{ID},m_k)={\displaystyle \sum _{k=1}^K}{c}_k\left(p_k,n_k^{ID},m_k\right).$$

(7)

2.5. Total Power Dissipation with Energy Harvesting

For the downlink transmission, the overall power consumption (Watt) in the network is given by [33]:

$$P(p_k,n_k^{EH},m_k)={\displaystyle \sum _{k=1}^K}\zeta P_k+p_c^{AP}+{\displaystyle \sum _{k=1}^K}{p}_{c,k}^{UE}-{\displaystyle \sum _{k=1}^K}{p}_k^{EH},$$

(8)

where $\zeta >1$ is the constant power inefficiency of the power amplifier. Also, the static circuit power dissipation $p_c^{AP}>0$ for the H-AP and $p_c^{UE}>0$ for the user occur in electronic components of the active transceiver, such as mixers, filters, and digital-to-analog converters, and are independent of the actual transmission power $p_k$. Note that $p_c^{UE}$ includes the power consumption for channel training and estimation with massive transmission antennas. On the other hand, the sum of harvested power $p_k^{EH}$ has a minus sign and works as the opposite of the power consumptions. Further, the harvested RF energy $p_k^{EH}$ indicates a part of the power radiated by the transmission antennas $m_k$, and can potentially be obtained by the user $k$’s rectennas $n_k^{EH}$ [4].

$$p_k^{EH}={\xi }_k\mathbb{E}[{||{\mathbf{G}}_k{\mathbf{v}}_k||}_2^2]={\xi }_k\mathrm{tr}({\mathbf{G}}_k{\mathbf{S}}_k{\mathbf{G}}_k^{\mathscr{H}})={\xi }_k{p}_k\mathrm{tr}\left({\mathbf{G}}_k{\mathbf{U}}_{{\mathbf{S}}_k}{\mathbf{G}}_k^{\mathscr{H}}\right),$$

(9)

where the constant parameter $0\le {\xi }_k\le 1$ is a user $k$’s loss in the energy transducer for converting RF to DC to be stored. Further, the noise is neglected in Equation (9), since, in practice, it is too small to be harvested. We use ${\mathbf{S}}_k=\mathbb{E}\left[{\mathbf{v}}_k{\mathbf{v}}_k^{\mathscr{H}}\right]$ to denote the covariance matrix of ${\mathbf{v}}_k$, and this can be rewritten as ${\mathbf{S}}_k=p_k{\mathbf{U}}_{{\mathbf{S}}_k}$, where ${\mathbf{U}}_{{\mathbf{S}}_k}$ is the all-ones matrix having the same dimensions of ${\mathbf{S}}_k$. For simplicity, we assume that ${\mathbf{G}}_k\cong {\overline{g}}_k{\mathbf{U}}_{{\mathbf{G}}_k}$, where ${\overline{g}}_k$ is an average channel gain of $g_{ij,k}$. Then, Equation (9) is re-evaluated as:

$$p_k^{EH}\cong {\xi }_k{|{\overline{g}}_k|}^2\mathrm{tr}\left({\mathbf{U}}_{{\mathbf{G}}_k}{\mathbf{U}}_{{\mathbf{S}}_k}{\mathbf{U}}_{{\mathbf{G}}_k^{\mathscr{H}}}\right)={\xi }_k{|{\overline{g}}_k|}^2{p}_k{m}_k\mathrm{tr}\left({\mathbf{U}}_{{\mathbf{G}}_k}{\mathbf{U}}_{{\mathbf{G}}_k^{\mathscr{H}}}\right)={\xi }_k{|{\overline{g}}_k|}^2{p}_k{m_k}^2{n}_k^{EH}.$$

(10)

2.6. Overall Network Energy Efficiency

Accordingly, the primal optimization problem for maximizing the overall network energy efficiency is formulated as:

$$\begin{array}{ccccccccccc}(\mathrm{P}1)\hfill & & & & & & & & & & \underset{\{p_k,n_k^{ID},m_k\}}{\max }\eta =\frac{C(p_k,n_k^{ID},m_k)}{P(p_k,n_k^{EH},m_k)},\end{array}$$

(11)

subject to:

$$\begin{array}{}(12)& (\mathrm{C}1)\hfill & & & & & & & & & & {\sum }_{k=1}^K{P}_k\le P_{max}, p_k\ge 0, \forall k\in K.(13)& (\mathrm{C}2)\hfill & & & & & & & & & & {\sum }_{k=1}^K\zeta P_k+p_c^{AP}\le p_{pg}, \forall k\in K.(14)& (\mathrm{C}3)\hfill & & & & & & & & & & {\sum }_{k=1}^K{m}_k\le M,\forall k\in K.(15)& (\mathrm{C}4)\hfill & & & & & & & & & & \mathrm{Pr}\left[c_k(p_k,n_k^{ID},m_k)\le c_k^{QoS}\right]\le {\epsilon }_c, \forall k\in K.(16)& (\mathrm{C}5)\hfill & & & & & & & & & & \mathrm{Pr}\left[p_k^{EH}\le p_k^{QoP}\right]\le {\epsilon }_p, \forall k\in K.(17)& (\mathrm{C}6)\hfill & & & & & & & & & & 1\le n_k^{ID}\le n_k^T, \forall k\in K.(18)& (\mathrm{C}7)\hfill & & & & & & & & & & n_k^{ID}+n_k^{EH}=n_k^T, \forall k\in K.\end{array}$$

The EE function as an objective function (P1) is constructed as a summation of the channel capacity per a summation of the consumed power for the H-AP and all UEs, and the units of (P1) are bits/Joule. Note that the (P1) is a non-convex optimization and maximization problem; we will discuss its transformed convex function for optimization in Section 3. Constraint (C1) implies that each transmission power of the H-AP to user $k$ is non-negative, and that its total power is limited, and cannot exceed the maximum transmission power $P_{max}$. Constraint (C2) limits the total power consumption of the network so that it does not exceed the maximum power supply from the power grid, $p_{pg}$. Constraint (C3) means that the sum of the $m_k$ cannot exceed the total number of antennas of the H-AP. Constraints (C4) and (C5) are the inequality of outage probability, to stochastically guarantee the minimum quality of service (QoS) and amount of harvested power (AoP), respectively, for each UE, where $c_k^{QoS}$ and $p_k^{QoP}$ are pre-defined requirements of channel capacity and harvested power, respectively, and ${\epsilon }_c\in (0,1)$ and ${\epsilon }_p\in (0,1)$ are the pre-defined outage probabilities. Constraints (C6) and (C7) are the numbers of switched rectennas and antennas, respectively, bounded by each user $k$’s total time-switching antennas. Further, we can obtain the optimal value using either the argument parameter $n_k^{ID}$ or $n_k^{EH}$, because $n_k^{ID}+n_k^{EH}$ is $n_k^T$, as is shown in Constraint (C7).

3. EE Optimization

In this section, nonlinear fractional programming [34] for transforming the non-convex objective function (P1) to the convex function in Section 3.1, and Lagrangian dual decomposition theory for resource allocations in Section 3.3, are used for optimization of the EE function. With these two optimization methods, the locally optimal values for the three arguments $p_k$, $n_k^{ID}$, and $m_k$, as well as for energy efficiency, can be successfully acquired. Also, through the proposed iteratively distributed algorithms, the locally optimal EE can be quickly and efficiently converged to.

3.1. Outer Loop Algorithm: Transformation of the Primal Objective Function

Since the objective function (P1) is a non-convex function, solving the function requires a brute force approach. To obtain an efficient optimization solution, we transform (P1) to the convex function using nonlinear fractional programming [34]. For simplicity, we define F as a set of feasible solutions of the optimization problem (P1), and the maximum EE as ${\eta }^{*}$. The theorem for achieving maximum EE can be proved as follows:

Theorem 1. 

${\eta }^{*}=\frac{C\left({p_k}^{*},{n_k^{ID}}^{*},{m_k}^{*}\right)}{P\left({p_k}^{*},{n_k^{EH}}^{*},{m_k}^{*}\right)}=\underset{\{p_k,n_k^{ID},m_k\}}{\max }\left\{\frac{C(p_k,n_k^{ID},m_k)}{P(p_k,n_k^{EH},m_k)}|\forall \{p_k,n_k^{ID},m_k\}\in F\right\}$ if and only if, $\underset{\{p_k,n_k^{ID},m_k\}}{\max }\{C(p_k,n_k^{ID},m_k)-{\eta }^{*}P(p_k,n_k^{EH},m_k)\}=C\left({p_k}^{*},{n_k^{ID}}^{*},{m_k}^{*}\right)-{\eta }^{*}P\left({p_k}^{*},{n_k^{EH}}^{*},{m_k}^{*}\right)=0$, for $C(p_k,n_k^{ID},m_k)\ge 0$ and $P(p_k,n_k^{EH},m_k)>0$.

This theorem shows that (P1) in fractional form has an equivalent function in subtractive form. Consequently, we can deal with the optimization problem as an equivalent objective function, and the equivalent form is used in the rest of this paper.

As a transformed objective function, the function (P2) is a convex and combinatorial optimization problem:

$$\begin{array}{ccccccccccc}(\mathrm{P}2)\hfill & & & & & & & & & & \underset{\{p_k,n_k^{ID},m_k\}}{\max }C(p_k,n_k^{ID},m_k)-\eta P(p_k,n_k^{EH},m_k)\end{array}$$

(19)

subject to: (C1), (C2), (C3), (C4), (C5), (C6), (C7), where $C(p_k,n_k^{ID},m_k)-\eta P(p_k,n_k^{EH},m_k)\ge 0$ is given for any value $\eta$ generated by the outer loop algorithm, which is shown in Algorithm 1, above. For a proof of Theorem 1, we define the equivalent objective function (P2) as $F(\eta )=\underset{\{p_k,n_k^{ID},m_k\}}{\max }\{C(p_k,n_k^{ID},m_k)-\eta P(p_k,n_k^{EH},m_k)\}.$

Algorithm 1. Outer Loop Algorithm for EE Maximization.
1:	Set initial input $\eta =0$, iteration index $T_{out}=0$, threshold 𝜏
2:	While 1 do
3:	Obtain optimum values of three arguments $\{{p_k}^{\prime },{n_k^{ID}}^{\prime },{m_k}^{\prime }\}$ through the inner loop algorithm for the given $\eta$
4:	If $C({p_k}^{\prime },{n_k^{ID}}^{\prime },{m_k}^{\prime })-\eta P({p_k}^{\prime },{n_k^{EH}}^{\prime },{m_k}^{\prime })<\tau$ (convergence verification)
5:	Return $\{{p_k}^{*},{n_k^{ID}}^{*},{m_k}^{*}\}=\{{p_k}^{\prime },{n_k^{ID}}^{\prime },{m_k}^{\prime }\}$ and obtain optimal EE ${\eta }^{*}=\frac{C({p_k}^{*},{n_k^{ID}}^{*},{m_k}^{*})}{P({p_k}^{*},{n_k^{EH}}^{*},{m_k}^{*})}$
6:	else
7:	Update $\frac{C({p_k}^{\prime },{n_k^{ID}}^{\prime },{m_k}^{\prime })}{P({p_k}^{\prime },{n_k^{EH}}^{\prime },{m_k}^{\prime })}$ and $T_{out}=T_{out}+1$
8:	end if
9:	end while

Proof of Theorem 1.

Convergence and maximization. □

Lemma 1.

$F(\eta )$ is strictly monotonic decreasing in $\eta$ , i.e., $F({\eta }^{\prime })$ > $F({\eta }^{″})$ if ${\eta }^{\prime }$ < ${\eta }^{″}$.

Proof. 

Let ${\eta }^{\prime }$ maximize $F({\eta }^{\prime })$, then $F({\eta }^{\prime })=\underset{\{p_k,n_k^{ID},m_k\}}{\max }\{C(p_k,n_k^{ID},m_k)-{\eta }^{\prime }P(p_k,n_k^{EH},m_k)\}=C({p_k}^{″},n_k^{{ID}^{″}},{m_k}^{″})-{\eta }^{\prime }P({p_k}^{″},n_k^{{EH}^{″}},{m_k}^{″})>C({p_k}^{″},n_k^{{ID}^{″}},{m_k}^{″})-{\eta }^{″}P({p_k}^{″},n_k^{{EH}^{″}},{m_k}^{″})\ge \underset{\{p_k,n_k^{ID},m_k\}}{\max }\{C(p_k,n_k^{ID},m_k)-{\eta }^{″}P(p_k,n_k^{EH},m_k)\}=F({\eta }^{″})$ where $C\left({p_k}^{″},n_k^{{ID}^{″}},{m_k}^{″}\right)-{\eta }^{\prime }P\left({p_k}^{″},n_k^{{EH}^{″}},{m_k}^{″}\right)>C\left({p_k}^{″},n_k^{{ID}^{″}},{m_k}^{″}\right)-{\eta }^{″}P\left({p_k}^{″},n_k^{{EH}^{″}},{m_k}^{″}\right)$ is reasonable, because ${\eta }^{\prime }$ is smaller than ${\eta }^{″}$, as stated in Lemma 1. □

Lemma 2.

Let any set $\left\{{p_k}^{\prime },n_k^{{ID}^{\prime }},{m_k}^{\prime }\right\}$ and ${\eta }^{\prime }=\frac{C\left({p_k}^{\prime },n_k^{{ID}^{\prime }},{m_k}^{\prime }\right)}{P\left({p_k}^{\prime },n_k^{{EH}^{\prime }},{m_k}^{\prime }\right)},$ then $F\left({\eta }^{\prime }\right)\ge 0$.

Proof. 

$F\left({\eta }^{\prime }\right)=\underset{\left\{p_k,n_k^{ID},m_k\right\}}{\max }\left\{C\left(p_k,n_k^{ID},m_k\right)-{\eta }^{\prime }P\left(p_k,n_k^{EH},m_k\right)\right\}\ge C\left({p_k}^{\prime },n_k^{{ID}^{\prime }},{m_k}^{\prime }\right)-{\eta }^{\prime }P\left({p_k}^{\prime },n_k^{{EH}^{\prime }},{m_k}^{\prime }\right)=0.$ □

As shown in Lemmas 1 and 2, it is natural that $F\left(\eta \right)$ converges to 0, because $F\left(\eta \right)$ is a monotonic decreasing and positive function as shown in Lemmas 1 and 2. Further, Algorithm 1 shows an iterative outer-loop algorithm for solving the optimization problem. The algorithm satisfies the conditions in Theorem 1. The convergence to the maximum EE is guaranteed if the inner problem (P2) in line 3 in the outer-loop algorithm can be solved. We showed that the energy efficiency $\eta$ increases in each iteration, and it converges to ${\eta }^{*}$ if $F\left(\eta \right)<\tau$. In conclusion, $F\left({\eta }_{T_{out}}\right)$ will eventually approach zero, and satisfy the conditions stated in Theorem 1. That means that the threshold 𝜏 should be set around positive zero.

3.2. Closed-Form Expression for Outage Constraints

The closed-form outage constraints can be expressed by substituting the $c_k$ of (6) into (C4), and $p_k^{EH}$ of (10) into (C5). Knowing that the channel frequency response ${\left|{\overline{h}}_{i,k}\right|}^2$ follows ${\left|{\overline{h}}_{i,k}\right|}^2~\mathrm{Exp}\left(\frac{m_k}{{\sigma }^2}\right)$, then ${\sum }_{i=1}^{N_{min}}{\rho }_{i,k}{\lambda }_{i,k}~\mathrm{Exp}\left(\frac{N_o{w}_k}{p_k{\sigma }^2{\sum }_{i=1}^{N_{min}}{\lambda }_{i,k}}\right)$ is an exponentially distributed random variable. Note that the rate parameter of ${\sum }_{i=1}^{N_{min}}{\rho }_{i,k}{\lambda }_{i,k}$ can be calculated by: $\mathrm{E}\left({\sum }_{i=1}^{N_{min}}{\rho }_{i,k}{\lambda }_{i,k}\right)=\mathrm{E}\left({\rho }_{1,k}{\lambda }_{1,k}+{\rho }_{2,k}{\lambda }_{2,k}+\dots +{\rho }_{N_{min},k}{\lambda }_{N_{min},k}\right)={\lambda }_{1,k}\mathrm{E}\left({\rho }_{1,k}\right)+{\lambda }_{2,k}\mathrm{E}\left({\rho }_{2,k}\right)+\dots +{\lambda }_{N_{min},k}\mathrm{E}\left({\rho }_{N_{min},k}\right)={\lambda }_{1,k}\frac{p_k{\sigma }^2}{N_o{w}_k}+{\lambda }_{2,k}\frac{p_k{\sigma }^2}{N_o{w}_k}+\dots +{\lambda }_{N_{min},k}\frac{p_k{\sigma }^2}{N_o{w}_k}={\sum }_{i=1}^{N_{min}}{\lambda }_{i,k}\frac{p_k{\sigma }^2}{N_o{w}_k}$ , where ${\rho }_{i,k}=\frac{P_k{\left|{\overline{h}}_{i,k}\right|}^2}{N_o{w}_k}~\mathrm{Exp}\left(\frac{N_o{w}_k}{p_k{\sigma }^2}\right)$. The outage probability for (C4) becomes:

$$\begin{array}{c}\mathrm{Pr}\left[N_{min}{w}_k{\log }_2\left(1+\frac{1}{N_{min}}\frac{{\sum }_{i=1}^{N_{min}}{\rho }_{i,k}{\lambda }_{i,k}}{m_k}\right)\le c_k^{QoS}\right]\le {\epsilon }_c\\ \iff \mathrm{Pr}\left[{\sum }_{i=1}^{N_{min}}{\rho }_{i,k}{\lambda }_{i,k}\le m_k{N}_{min}\left(e^{\frac{c_k^{QoS}\ast ln2}{N_{min}{w}_k}}-1\right)\right]\le {\epsilon }_c\\ \iff 1-\exp \left(-\frac{N_o{w}_k}{p_k{\sigma }^2{\sum }_{i=1}^{N_{min}}{\lambda }_{i,k}}\ast m_k{N}_{min}\left(e^{\frac{c_k^{QoS}\ast ln2}{N_{min}{w}_k}}-1\right)\right)\le {\epsilon }_c.\end{array}$$

(20)

$$\iff \psi \left(p_k,n_k^{ID},m_k\right)=ln\left(1-{\epsilon }_c\right)\ast \frac{p_k{\sigma }^2{\sum }_{i=1}^{N_{min}}{\lambda }_{i,k}}{N_o{w}_k{m}_k{N}_{min}}+e^{\frac{c_k^{QoS}\ast ln2}{N_{min}{w}_k}}-1\le 0.$$

(21)

Using the same process of (20), (21), (C5) is expressed as:

$$\omega \left(p_k,n_k^{EH},m_k\right)=ln\left(1-{\epsilon }_p\right){\xi }_k{p}_k{m}_k{\sigma }^2+p_k^{QoP}\le 0,$$

(22)

where ${\left|{\overline{g}}_k\right|}^2~\mathrm{Exp}\left(\frac{n_k^{EH}{m}_k}{{\sigma }^2}\right)$.

3.3. Inner Loop Algorithm: Resource Allocation

To be practical, we will focus on an efficient distributed algorithm by Lagrange dual decomposition. The Lagrange dual decomposition is obtained by using non-negative Lagrange multipliers to involve the described constraints (C1)–(C7) in the objective function (P1). This solution gives rise to a new problem, which is that of maximizing the objective function with respect to the dual variables under the derived constraints on the dual variables. Therefore, by its decomposability, the dual decomposition algorithm has advantages over exhaustive dynamic programming in terms of both efficiency and simplicity when strong duality holds. The dual problem formulation of (P2) is constructed as:

$$\begin{array}{l}ℒ(p_k, n_k^{ID}, m_k, \alpha , \beta , \gamma , {\delta }_k, {\epsilon }_k)\\ ={\displaystyle \sum _{k=1}^K}{c}_k(p_k, n_k^{ID}, m_k)-\eta \left({\displaystyle \sum _{k=1}^K}\zeta m_k{p}_k+p_c^{AP}+{\displaystyle \sum _{k=1}^K}{p}_{c,k}^{UE}-{\displaystyle \sum _{k=1}^K}{p}_k^{EH}\right)\\ -\alpha \left({\displaystyle \sum _{k=1}^K}{m}_k{p}_k-P_{\mathit{max}}\right)-\beta \left({\displaystyle \sum _{k=1}^K}\zeta m_k{p}_k+p_c^{AP}-p_{pg}\right)\\ -\gamma \left({\displaystyle \sum _{k=1}^K}{m}_k-M\right)-{\displaystyle \sum _{k=1}^K}{\delta }_k\psi (p_k, n_k^{ID}, m_k)-{\displaystyle \sum _{k=1}^K}{\epsilon }_k\omega (p_k, n_k^{EH}, m_k),\end{array}$$

(23)

where $\alpha , \beta , \gamma , {\delta }_k$ and ${\epsilon }_k$ are Lagrange multipliers. We rewrite Equation (23) as Algorithm 2:

Algorithm 2. Inner Loop Algorithm for Resource Allocation.
1:	Initialize: Active UE set $k=\left\{1,2,\dots ,K\right\}, {\mu }_1, {\mu }_2, {\mu }_3, {\mu }_4,n_k^{ID}=1,n_k^{EH}=n_k^T-n_k^{ID},\alpha =\beta ={\gamma }_k={\delta }_k=p_k=0.$
2:	while $p_k,n_k^{ID},m_k$ are not converged do
3:	for $k\in K$ do
4:	if $n_k^{ID}<m_k \left(case 1\right)$ do
5:	Update ${\widehat{p}}_k=p_k\left(T_{in}+1\right)$ by jointly solving (39)~(42)
6:	Update ${\widehat{m}}_k=m_k\left(T_{in}+1\right)$ by jointly solving (39)~(42) using the calculated value ${\widehat{p}}_k$
7:	Update ${\widehat{n}}_k^{ID}=n_k^{ID}\left(T_{in}+1\right)+1/2$ by jointly solving (39)~(42) using the calculated value ${\widehat{p}}_k$ and ${\widehat{m}}_k$
8:	else if $n_k^{ID}\ge m_k \left(case 2\right)$ do
9:	Update ${\widehat{p}}_k=p_k\left(T_{in}+1\right)$ by jointly solving (39)~(42)
10:	Update ${\widehat{m}}_k=m_k\left(T_{in}+1\right)$ by jointly solving (39)~(42) using the calculated value ${\widehat{p}}_k$
11:	Update ${\widehat{n}}_k^{ID}=n_k^{ID}\left(T_{in}+1\right)+1/2$ by jointly solving (39)~(42) using the calculated value ${\widehat{p}}_k$ and ${\widehat{m}}_k$
12:	end if
13:	Update the subgradient of ${\delta }_k\left(T_{in}+1\right)$ with $S_{{\delta }_k}\left(T_{in}\right)$
14:	Update the subgradient of ${\epsilon }_k\left(T_{in}+1\right)$ with $S_{{\epsilon }_k}\left(T_{in}\right)$
15:	end for
16:	Update the subgradient of $\alpha \left(T_{in}+1\right)$ with $S_{\alpha }\left(T_{in}\right)$
17:	Update the subgradient of $\beta \left(T_{in}+1\right)$ with $S_{\beta }\left(T_{in}\right)$
18:	Update the subgradient of $\gamma \left(T_{in}+1\right)$ with $S_{\gamma }\left(T_{in}\right)$
19:	Let $T_{in}=T_{in}+1$
end while

$$\begin{array}{l}ℒ(p_k,n_k^{ID},m_k,\alpha ,\beta ,\gamma ,{\delta }_k,{\epsilon }_k)\\ ={\sum }_{k=1}^K{ℒ}_k(p_k,n_k^{ID},m_k,\alpha ,\beta ,\gamma ,{\delta }_k,{\epsilon }_k)-\eta p_c^{AP}+\alpha P_{max}-\beta (p_c^{AP}-p_{pg})+\gamma M,\end{array}$$

(24)

where ${ℒ}_k\left(p_k,n_k^{ID},m_k,\alpha ,\beta ,\gamma ,{\delta }_k,{\epsilon }_k\right)=c_k\left(p_k,n_k^{ID},m_k\right)-\eta \left(\zeta m_k{p}_k+p_{c,k}^{UE}-p_k^{EH}\right)-\alpha m_k{p}_k-\beta \zeta m_k{p}_k-\gamma m_k-{\delta }_k\psi \left(p_k,n_k^{ID},m_k\right)-{\epsilon }_k\omega \left(p_k,n_k^{EH},m_k\right).$ The dual optimization problem of Equation (24) with respect to variables $\alpha , \beta , \gamma , {\delta }_k$ and ${\epsilon }_k$ is:

$$\underset{\alpha \ge 0}{\min }D\left(\alpha \right)={\sum }_{k=1}^K{D}_k\left(\alpha \right)+\alpha P_{max}.$$

(25)

$$\underset{\beta \ge 0}{\min }Q\left(\beta \right)={\sum }_{k=1}^K{Q}_k\left(\beta \right)-\beta \left(p_c^{AP}-p_{pg}\right).$$

(26)

$$\underset{\gamma \ge 0}{\min }U\left(\alpha \right)={\sum }_{k=1}^K{U}_k\left(\gamma \right)+\gamma M.$$

(27)

$$\underset{{\delta }_k\ge 0}{\min }{V}_k\left({\delta }_k\right).$$

(28)

$$\underset{{\epsilon }_k\ge 0}{\min }{Z}_k\left({\epsilon }_k\right).$$

(29)

As (P2) is convex maximization, strong duality holds. We can solve the dual problems of Equations (25–29) in an iterative manner, using the gradient method [35]. The subgradient of $\alpha$ in subproblem $D$ is:

$$S_{\alpha }=\frac{\partial D\left(\alpha \right)}{\partial \alpha }={\sum }_{k=1}^K{m}_k{p}_k+P_{max},$$

(30)

where $D_k\left(\alpha \right)$ in Equation (25) can be sorted as the first term by the descending order, and is partially differentiated by $\alpha$. Therefore, composing the equation with either all terms, or the first term in $\alpha$ except the constant, doesn’t cause any problems, and it can be easily calculated at the H-AP. Each subgradient of $\beta , \gamma , {\delta }_k$ and ${\epsilon }_k$ is:

$$S_{\beta }=\frac{\partial Q\left(\beta \right)}{\partial \beta }={\sum }_{k=1}^K\zeta m_k{p}_k-p_c^{AP}+p_{pg},$$

(31)

$$S_{\gamma }=\frac{\partial U\left(\gamma \right)}{\partial \gamma }={\sum }_{k=1}^K{m}_k+M,$$

(32)

$$S_{{\delta }_k}=\frac{\partial V_k\left({\delta }_k\right)}{\partial {\delta }_k}=\psi \left(p_k,n_k^{ID},m_k\right),$$

(33)

$$S_{{\epsilon }_k}=\frac{\partial Z_k\left({\epsilon }_k\right)}{\partial {\epsilon }_k}=\omega \left(p_k,n_k^{EH},m_k\right).$$

(34)

The parameters $p_k$ and $n_k^{ID}$ in subgradients $S_{{\gamma }_k}$ and $S_{{\delta }_k}$ are easy to calculate at the H-AP. Then, the multipliers are updated. During the iteration, Lagrange multipliers $\alpha , \beta , \gamma , {\delta }_k$ and ${\epsilon }_k$ are updated in a distributed manner as:

$$\alpha \left(T_{in}+1\right)={\left[\alpha \left(T_{in}\right)-{\mu }_1{S}_{\alpha }\left(T_{in}\right)\right]}^{+},$$

(35)

$$\beta \left(T_{in}+1\right)={\left[\beta \left(T_{in}\right)-{\mu }_2{S}_{\beta }\left(T_{in}\right)\right]}^{+},$$

(36)

$$\gamma \left(T_{in}+1\right)={\left[\gamma \left(T_{in}\right)-{\mu }_3{S}_{\gamma }\left(T_{in}\right)\right]}^{+},$$

(37)

$${\delta }_k\left(T_{in}+1\right)={\left[{\delta }_k\left(T_{in}\right)-{\mu }_4{S}_{{\delta }_k}\left(T_{in}\right)\right]}^{+},$$

(38)

$${\epsilon }_k\left(T_{in}+1\right)={\left[{\epsilon }_k\left(T_{in}\right)-{\mu }_5{S}_{{\epsilon }_k}\left(T_{in}\right)\right]}^{+},$$

(39)

where ${\left[Z\right]}^{+}=\max \left\{Z, 0\right\}$. The parameter $T_{in}$ is the number of iterations of the inner loop in which the H-AP updates $\alpha , \beta , \gamma , {\delta }_k$ and ${\epsilon }_k$. The constant coefficients ${\mu }_1$, ${\mu }_2$, ${\mu }_3$, ${\mu }_4$ and ${\mu }_5$ are positive values like a learning rate to converge faster.

$$N_{min}=\{\begin{array}{l}{n}_k^{ID}, n_k^{ID}<m_k, \mathrm{case} 1\\ m_k, n_k^{ID}\ge m_k, \mathrm{case} 2\end{array}$$

By the Karush-Kuhn-Tucker (KKT) condition, the optimal transmission power ${\widehat{p}}_k$ is calculated by the following equality:

$$\frac{\partial {ℒ}_k\left(p_k,n_k^{ID},m_k,\alpha ,\beta ,\gamma ,{\delta }_k,{\epsilon }_k\right)}{\partial p_k}=0,$$

(40)

and the optimal number of receiving antennas ${\widehat{n}}_k^{ID}$ and number of transmission antennas ${\widehat{m}}_k$ are easily calculated by jointly solving Equations (39–42) for the two cases of $N_{min}$:

$$\frac{\partial {ℒ}_k\left(p_k,n_k^{ID},m_k,\alpha ,\beta ,\gamma ,{\delta }_k,{\epsilon }_k\right)}{\partial n_k^{ID}}=0.$$

(41)

$$\frac{\partial {ℒ}_k\left(p_k,n_k^{ID},m_k,\alpha ,\beta ,\gamma ,{\delta }_k,{\epsilon }_k\right)}{\partial m_k}=0.$$

(42)

Algorithm 1 gives the pseudo-code of the inner loop algorithm.

4. Simulation Results

In this section, we provide numerical results to show the effectiveness of the proposed algorithm. The simulation scenario is set as a small cell network of wireless personal area network size using 38 (GHz) as the center frequency, and assumes that an H-AP with a coverage of 10 m radius connects with less than 10 active users, which is a scenario that can generally happen in an office space in a building. In wireless power transfer (WPT), the coverage is very limited because of the very low efficiency of WPT. It is expected that the coverage of WPT “is practically” (please see the DEMO posted to the link of youtube: https://youtu.be/qP9fZQX1sDk with [36,37]) up to 10 m for wireless-powered sensor applications. This is why we consider such a small cell network for WPT. In the same spirit, we can improve such low efficiency of WPT by using a massive antenna array. For this reason, we set a small cell network randomly deployed with the massive antenna array at the hybrid access point. Further, the proposed scenario where a small number of users access a single access point with the massive MIMO is set considering the trend that wireless communication networks are evolving toward small cell networks and large-scale MIMO systems. For the massive MIMO property, the number of transmit antennas at the H-AP, $M$, and the number of equipped TS antennas for user k, $n_k^T$, are set to 80 and 16~32, respectively. The specific settings are listed in

Table 1. Simulation Parameters.

Parameter	Value
Number of users, k	8
Coverage of H-AP	10 (m)
Three-dimensional location of H-AP	(0, 0, 0)
Three-dimensional location of users	(0, 5, 0), (−4, −5, 0), (8, −3, −1), (−6, 5, 0), (3, 3, 1), (1, 1, 0), (1, −3, 0), and (−2, −9, 1) (m)
Distance from H-AP for user k	4, 6.4, 8.6, 8.48, 4.36, 1.41, 2.24, and 9.27 (m)
Number of transmit antennas at the H-AP, M	80
Number of TS antennas for user k, $n_k^T$	16, 32
Initial transmission power, $p_k(1)$	$2.25 (\mathrm{mW})$
Outage probability for C4 and C5, ${\epsilon }_c$ and ${\epsilon }_p$	0.15, respectively
Bandwidth for user k, $w_k$	$300 (\mathrm{kHz})$
Static circuit power dissipation at H-AP, $p_c^{AP}$	$34 (\mathrm{dBm})$
Static circuit power dissipation at user k, $p_c^{UE}$	$25 (\mathrm{dBm})$
Maximum transmission power, $P_{max}$	$1300 (\mathrm{mW})$
Maximum power supply from the power grid, $p_{pg}$	$47 (\mathrm{dBm})$
${\xi }_k$	0.8
Target channel capacity, $c_k^{QoS}$	$14 (\mathrm{Mbps})$
Target amount of harvested energy, $p_k^{QoP}$	$1200 (\mathsf{\mu }\mathrm{W})$
Center frequency, $f_c$	$38 (\mathrm{GHz})$
Noise variance, $N_0{w}_k$	−119.23 (dBm)
Noise spectral density, $N_0$	−174 (dBm/Hz) at 290 degree Kelvin.
Power inefficiency of the power amplifier, $\zeta$	5

.

Figure 2 shows the distribution of the eight users around the H-AP as an RF-EHN. In this map, the H-AP estimates and distributes the optimum value of the three arguments, ${\widehat{p}}_k$, ${\widehat{n}}_k^{ID}$ and ${\widehat{m}}_k$, to users using the proposed inner and outer algorithms every scheduled time. Further, each user simultaneously receives the transmitted signal from the H-AP for the EH and ID.

Figure 3 represents the convergence performance of transmission power according to each iteration in terms of $p_k$ and $p_k{m}_k$($=P_k)$. Note that the convergence performances of the proposed algorithm for Inner Loop Algorithm and Outer Loop Algorithm are presented in Figure 3, Figure 4, Figure 5 and Figure 6. The optimization is in progress following the QoS and AoP to obtain the maximum EE from the initial transmission power value of 3.5 (dBm). From the figure, it is confirmed that convergence performance is different for each user but that convergence is achieved in two to six iterations as a whole. Some curves seem to converge in the second iteration but, in fact, the value changes slightly when the iteration goes on. Generally, it seems unreasonable that convergence would occur in only the second iteration, and it seems reasonable that it would require several iterations. Further, it is confirmed that the system operates within the pre-defined maximum transmission power constraint of 1300 mW. This means that iteration continues to keep the constraint C1. Note that for precise performance analysis, simulations in this paper show the number of antennas as real numbers, rather than natural numbers.

Figure 4 represents the performance of antenna and rectenna convergence in accordance with each iteration based on the proposed inner loop algorithm, allocating three resources at optimal values. Also, we can see that most antennas effectively converge within five iterations. As the number of iterations increases, each user’s antenna and rectenna curve varies to satisfy the pre-defined QoS and AoP, and the sum of the antennas and rectennas for the two curves of each user cannot exceed its equipped total number of switching antennas. When the numbers of antennas and rectennas have converged, all users are perfectly guaranteed their QoS and AoP.

Figure 5 shows the convergence performance of the number of allocated transmission antennas to each user. Each curve simultaneously maximizes the EE, and meets the user’s QoS and AoP. This demonstrates that to satisfy their QoS and AoP, users far from the H-AP need to be allocated up to five times more transmission antennas than near ones. Also, these converged values do not exceed the maximum transmission power limitation of C1 with regard to the results of Figure 5. This means that we designed the system for practical resource allocations, so that the transmission antennas cannot be infinitely allocated to users.

Figure 6 represents energy efficiency performances as a network EE in the proposed simulation settings, and every point of the energy efficiency is averaged over 1000 independent channel realizations. The “iteration” in this figure refers to the iteration index $T_{out}$ of the Outer Loop Algorithm in Algorithm 1. Also, the energy efficiency is calculated by the proposed the inner and outer loop Algorithms, and is constructed as the ratio of each user’s channel capacity per consumed electric power containing the harvested energy. The energy efficiency is locally optimized, with suboptimal values, as the iterations continue. Furthermore, it is based on the presented optimization curves for the three arguments, as in Figure 3, Figure 4 and Figure 5. In the experimental environments of these figures, the energy efficiency has a maximum value when all of the users are equipped with 32 time-switching antennas. However, it is observed that, although the number of equipped antennas is doubled, the energy efficiency rises to about 1.6 times the values on the curve for 8 users equipped with 16 time-switching antennas. It is instructive to see that this is because of insufficient transmission antenna allocations, due to the limit of transmission antenna allocation capabilities, and the increases of circuit power consumption that occurred in 4 users equipped with 32 antennas. Consequently, with regard to meeting the goal of EE maximization, the overall behavior of the graphs in the figures, including Figure 6, is seen as reasonable.

Figure 7 shows various target outage probabilities ${\epsilon }_c, {\epsilon }_p\in \left[0.05,0.35\right]$ versus real outage probabilities to make sure that the outage constraints (C4) and (C5) work fine in the proposed algorithm. The results of Figure 7 are averaged over 1000 independent Monte-Carlo simulations, and each curve involves various channel and noise parameters. It is observed that the real outage probabilities are always lower than the target outage probabilities, which means that the outage constraints work accurately.

Figure 8 represents energy efficiency versus different ratios of ID and EH modes for all users. The ratio $\theta$ $\in \left[0,1\right]$ represented in Equation (43) means that the ratio of switching between ID and EH mode for the total number of antennas that each user has:

$$n_k^T=\theta \times n_k^T+\left(1-\theta \right)n_k^T=n_k^{ID}+n_k^{EH}.$$

(43)

For comparison with the performance of the proposed algorithm in terms of TS antenna switching strategy $n_k^{ID}$, we set baseline scheme A that optimizes energy efficiency only in terms of two resources $\left(p_k,m_k\right)$ (the proposed algorithm optimizes all proposed resources $\left(p_k,n_k^{ID},m_k\right)$). The baseline scheme A has no TS antenna switching strategy but has various ratios $\theta$ $\in \left[0,1\right]$ for simulations. Note that $\theta$ is applied to all users equally, and energy efficiency becomes zero when $\theta$ is zero, since the channel capacity is not ensured at all. From the figure, it is observed that there are performance differences between the proposed algorithm and the baseline scheme in terms of energy efficiency, and 120% of the baseline scheme’s peak energy efficiency is as much as that of the proposed algorithm for the two cases that $n_k^T$ is 16 and 32. This means that the fixed TS antenna switching ratio of the baseline scheme cannot maximize energy efficiency, so it is necessary to calculate the optimum switching ratio for each user as shown in Figure 4 in order to achieve maximum energy efficiency.

Figure 9 illustrates energy efficiency versus the total transmission power constraint for $k=8$ users. The transmission power constraint is set as a series of levels from 900 to 1700 (mW) with intervals of 10 (mW). Note that all energy efficiencies are estimated regardless of achieving the user’s target QoS and AoP, and this is done only to see how energy efficiency behaves when maximum transmission power varies. We can observe that all curves are upper convex, and this is mainly due to linear scaling of the transmission power in Equation (7)’s first term. For comparison, we set baseline scheme B that allocates the equal number of transmission antennas 1, 2, 5, and 10 to each user. In other words, the baseline scheme B optimizes energy efficiency only in terms of two resources $(p_k,n_k^{ID}),$ while the proposed algorithm optimizes energy efficiency in terms of $(p_k,n_k^{ID},m_k).$ From the figure, it is observed that there are performance differences between the proposed algorithm and the baseline scheme ($m_k=10$), and 137% of the baseline scheme’s EE is as much as the proposed algorithm’s EE at a maximum gap point of 1440 [mW] (the proposed algorithm allocates its total 80 transmission antennas to users for fair comparison). This means that fixed transmission antenna policies of baseline schemes degrade the energy efficiency performance, and this is because more power is required to satisfy the minimum QoS and AoP. On the other hand, although the two schemes use the same transmission antennas, the reason that the proposed scheme has better EE performance can be seen to be because it allocates optimized numbers of transmission antennas to users.

Figure 10 demonstrates the performance of energy efficiency versus variable QoS and AoP constraints. It is observed that the curve of the proposed algorithm ($p_k^{QoP}=1200 (\mathsf{\mu }\mathrm{W}), P_{max}=1300 (\mathrm{mW})$) has its EE values of up to 17 (Mbps) of QoS. Also, it is impossible to calculate the EE with both QoS over 17 (Mbps) and AoP over 1650 (μW) satisfied, because the resources are not sufficient. However, if it has more maximum transmission power—as much as 1500 (mW)—it can achieve a little more energy efficiency and QoS. Additionally, if the RF energy harvesting network sets more AoP—as much as 150 (μW)—users may get more battery charge. However, that can result in degrading the achievable energy efficiency, and decreasing the settable values for the QoS as a trade-off.

5. Conclusions

In this paper, we formulated energy efficiency optimization for MU massive MIMO RF-EHNs as a mixed non-convex problem. To solve the problem, minimum channel capacity requirements, a minimum harvested power, circuit power consumption, limit of transmission antenna allocation, and ID and EH switching mechanisms were taken into account. To simultaneously guarantee each user’s QoS and AoP and maximize the energy efficiency, two iterative optimization algorithms were proposed, based on optimizing the number of achievable bits-per-Joule in RF-EHNs. Simulation results showed that the proposed algorithm converges to the solution within a few iterations and demonstrated that the proposed algorithms are more energy-efficient than the baseline schemes.