Joint spectral efficiency optimization of uplink and downlink for massive MIMO-enabled wireless energy harvesting systems

Abstract

This paper investigated the spectral efficiency (SE) in massive multiple-input multiple-output systems, where all terminals have no fixed power supply and thus need to replenish energy via the received signals from the base station. The hybrid wireless energy harvesting (EH) protocol is applied for each terminal, which can switch to either existing time-switching (TS) protocol or power-splitting (PS) protocol. Based on the hybrid wireless EH protocol, a general system model is developed, which can switch to either only uplink data transmission or only downlink data transmission. As a result, a general analytical framework is formulated. Then, closed-form lower bound expressions on SE for each terminal are obtained on the uplink and downlink, respectively. According to these expressions, the joint SE of uplink and downlink maximization problem is designed with some practical constraints. As the designed optimization problem is non-linear and non-convex, it is hard to solve directly. To provide a solution, an iteration algorithm is proposed by utilizing one-dimensional search technique and successive approximation method based on geometric program. Additionally, the convergence and complexity of the proposed algorithm are discussed as well. Finally, the feasibility of the proposed algorithm is analyzed by simulations. Numerical results manifest that the proposed algorithm can provide good SE by optimizing relevant system parameters, and the system model can help to discuss the TS, PS or hybrid protocol for only uplink data transmission, only downlink data transmission or joint data transmission of uplink and downlink in the considered system.

1 Introduction

The proliferation of wireless networks leads to a large number of smart terminals. For some special applications such as military communications, medical implants and mining, one of the key challenges in such networks is how to supply sufficient energy to network terminals for successful communications. Motivated by this, the energy harvesting (EH) technique has gained much attention in both academia and industry as it provides a promising approach to prolong the lifetime of energy-constrained terminals [1,2,3,4]. Although some natural energy sources such as sunshine and wind can be considered, they are usually not as effective as expected due to the inconsistent and unforeseeable nature of ambient sources [5]. Compared with natural energy sources, radio frequency (RF) EH is nominated as the best EH scheme due to the stability and the ability of transforming the RF signals into voltage to charge the terminal battery. To perform RF EH, two main protocols adopted at terminals are time-switching (TS) protocol and power-splitting (PS) protocol.

One transmission scheme based on TS protocol is wireless powered communication network (WPCN) in which each terminal harvests power on the downlink in the first slot and then transmits its information on the uplink in the second slot. In this scheme, TS ratio is an important parameter to evaluate system performance. For example, in broadcast channels scenario, the authors in [6] focused on the trade-off of wireless energy and information transfer by adjusting TS ratio, and the authors in [7] maximized the minimum throughput among all terminals by optimizing the downlink/uplink TS ratio, the downlink energy beamforming, and the uplink power allocation (PA) as well as receive beamforming. In relay channels, the authors in [8] maximized instantaneous throughput by optimizing beamformers and TS ratio at the relay, and the authors in [9] addressed the problems of maximizing throughput for fixed supplementary battery energy and minimizing supplementary battery energy consumption for target throughput performance by optimizing the TS ratio.

Another transmission scheme, simultaneous wireless information and power transfer (SWIPT), is based on PS protocol. In this scheme, a transmitter uses the same waveform to transfer wireless energy and information, and each terminal divides the received signal power by PS ratio between EH and information decoding (ID). Here, PS ratio is of great importance because it can directly affect the system performance. For example, the authors in [10] jointly designed transmit beamforming vector, PS ratio and transmit power to minimize the weighted sum transmit power in full-duplex (FD) networks. The authors in [11] designed the optimal PS and TS ratios to maximize the weighted sum rate over all users under some constraints in orthogonal frequency division multiplexing (OFDM) systems. To study the performance of relay networks based on wireless EH, PS and TS ratios are also intensively investigated to obtain the optimal throughput in an amplify-and-forward (AF) relaying network [12] and a decode-and-forward (DF) relaying network [13].

In practice, RF signals decay quickly over a long distance. One feasible solution to that is utilizing energy beamforming that can focus RF signals into a narrow beam to enhance transmission efficiency [14]. Based on this fact, massive multiple-input multiple-output (MIMO) is considered as a strong candidate for energy beamforming as the large-scale antenna array equipped at base station (BS) can provide a sharp beam to enhance the received signal strength at each terminal [15, 16]. Therefore, the combination of RF EH and massive MIMO is a more practical scheme and attracts intensive research interests.

Among them, the authors in [17] studied a massive MIMO WPCN where PA weights, energy-splitting fraction and TS ratio were optimized respectively for maximizing the minimum rate among terminals on the uplink. The authors in [18] investigated a downlink multi-user massive MIMO system based on SWIPT where an iterative algorithm was proposed to optimize PA coefficients at BS and PS ratios at all terminals for maximizing the minimum achievable rate among all terminals on the downlink. To study the influence of line of sight path, the authors in [19] investigated the downlink transmission of massive MIMO-enabled SWIPT systems in Rician fading channels where PA, channel estimation duration and PS ratios were optimized to maximize the downlink sum rate and the minimum rate among all terminals, respectively. In [20], the authors investigated the beam-domain SWIPT in a massive MIMO system in which the transmit power at BS and the TS ratios at all terminals were optimized under the constraints of the current available energy and minimum transmission rate of terminals for achieving maximum sum rate. In integrated data and energy communication networks, the authors optimized the resource allocation and PS ratios for uplink throughput maximization in [21]. The optimal energy-rate trade-offs based on SWIPT were investigated in a relay-assisted downlink massive MIMO system [22].

Nevertheless, the aforementioned studies [17,18,19,20,21,22] have been limited to investigating the system performance on the uplink or on the downlink via TS protocol or PS protocol. Moreover, the authors in [23] proposed a novel hybrid wireless EH protocol which is a combination of TS and PS protocols in wireless relay networks. Based on these observations and inspired by the analyses of [17,18,19,20,21,22,23], in this paper, we extend the existing system model from only uplink or downlink transmission to joint uplink and downlink transmission, and each terminal is equipped with the hybrid wireless EH protocol, which is a combination of TS protocol and PS protocol. Our aim is to investigate joint spectral efficiency (SE) of uplink and downlink with some practical constraints. To the best of our knowledge, such the study has not been found in the existing studies yet. The main contributions of this work are summarized as follows:

• We propose a joint uplink and downlink transmission scheme in massive MIMO systems and each terminal is equipped with the hybrid wireless EH protocol. On the downlink, the BS delivers RF signals to all terminals and each terminal uses the hybrid wireless EH protocol to coordinates the processes of EH and ID. On the uplink, a fraction of the harvested energy is used for uplink pilot transmission and the remaining fraction is used for uplink data transmission. Specially, this scheme provides unified system model as it can switch to only uplink data transmission or only downlink data transmission. The hybrid wireless EH protocol can run in three modes, i.e., TS, PS or hybrid protocol.

• An optimization problem is designed to maximize joint SE of uplink and downlink for massive MIMO systems under some practical constraints. As the designed optimization problem is non-convex and non-linear, which poses huge challenges to solve directly, an algorithm is proposed to solve such a complex problem by utilizing one-dimensional search method and successive approximation method based on geometric program (GP). In addition, the proposed algorithm is illustrated in details with rigorous mathematic analyses and its computational complexity and convergence are also discussed.

• Numerical results reveal that the performance of TS protocol is far below that of PS protocol and hybrid wireless EH protocol, and that the performance of the hybrid wireless EH protocol closely approaches that of PS protocol for the joint SE of uplink and downlink maximization problem in the considered system. Moreover, time resources have a more significant impact on system performance than energy resources. In addition, wireless EH protocol, channel state information (CSI), PA coefficients and the number of antennas at BS, and energy allocation ratio at each terminal are all effective means to improve the system performance.

The rest of this paper is organized as follows: Section 2 briefly describes the system and signal model, and then formulates the joint SE of uplink and downlink optimization problem. In Sect. 3, an algorithm is proposed for solving this optimization problem. Additionally, the complexity of the propose algorithm is also discussed in this section. Furthermore, numerical results are conducted to demonstrate our proposed algorithm in Sect. 4. Finally, we conclude the whole paper in Sect. 5.

Notations: Scalars are denoted by lowercase or uppercase letters. Vectors and matrices are denoted by bold lowercase and bold uppercase letters, respectively. $\mathbf{I}_K$ and $\mathbf{I}_M$ are the identity matrices of size $K\times K$ and $M\times M$, respectively. The operator $\mathbb{E}\{\cdot \}$ stands for the expectation of a random variable. The notation $\Vert \cdot \Vert$ represents the Euclidean norm. The Hermitian and regular transpose are denoted by $(\cdot )^H$ and $(\cdot )^T$, respectively. Finally, $\mathcal {CN}(.,.)$ is the circular symmetric complex-Gaussian distribution.

2 System model and problem statement

2.1 System and signal model

Fig. 1

A massive MIMO uplink and downlink system where each terminal operates via the hybrid wireless EH protocol

We consider a massive MIMO-enabled wireless EH system as shown in Fig. 1, where the BS employing a compact array of M antennas communicates simultaneously with K active single-antenna terminals. It is assumed that the BS is connected to a continuous stable power supply while each terminal can be empowered by the energy harvested from the received RF signals via the hybrid wireless EH protocol in [23] as shown in Fig. 2. We also assume that the channels between the BS and all terminals are constant and frequency-flat in each frame, and the system operates in a time-division-duplexing (TDD) mode. In detail, each frame lasts T seconds and consists of four phases based on the idea of the hybrid wireless EH protocol. Their operations are explained below.

• In the phase of uplink pilots, all terminals simultaneously transmit mutually orthogonal uplink pilots to the BS and then the BS estimates uplink channels. By exploiting uplink channels reciprocity, the downlink CSI is obtained easily. This phase lasts $\tau T (0\le \tau \le 1)$ seconds.

• In the phase of length $\alpha T (0 \le \alpha \le 1-\tau )$ seconds, based on the idea of TS protocol, the BS simply transfers power to all terminals without any information exchange and each terminal charges its own battery. This phase is named downlink wireless power transfer (WPT).

• In the phase of length $(1-\tau -\alpha )T/2$ seconds, based on the idea of PS protocol, the BS transmits information and enery using the same RF signals to all terminals and each terminal divides the received signal power into two parts. A $\rho (0 \le \rho \le 1)$ part of the power is used for ID and the remaining $1-\rho$ part is used for EH. This phase is named downlink SWIPT.

• In the phase of uplink data, the total harvested energy in the phases of downlink WPT and SWIPT is split such that a $\lambda (0 \le \lambda \le 1)$ fraction of total harvested energy is used to transmit uplink pilots and the remaining $1-\lambda$ fraction is used to send uplink data. This phase lasts $(1-\tau -\alpha )T/2$ seconds.

Fig. 2

A general transmission frame structure for the hybrid wireless EH protocol

When $\lambda =1$ or $\rho =0$, the joint transmission system model of uplink and downlink can switch to only downlink transmission system model as shown in [18, 19] or only uplink transmission system model as shown in [17], respectively. When $\alpha = 0$ or $\rho = 1$, the hybrid wireless EH protocol is converted to the PS or TS protocol, respectively. Therefore, this paper provides unified system model and analytical model for massive MIMO-enabled wireless EH system.

It is assumed that the channels between all terminals and BS antennas follow independent and identically distributed (i.i.d.) Rayleigh fading. Let ${\mathbf{g}}_{k}=\sqrt{\beta _{k}}{\mathbf{h}}_{k}\in \mathbb{C}^{M\times 1}$, for $k=1\cdots K$, denote the channel between the kth terminal and BS antennas, where $\beta _{k}$ represents the large-scale fading coefficient and ${\mathbf{h}}_{k}\in \mathbb{C}^{M\times 1}$ contains the i.i.d. $\mathcal {CN}(0,1)$ small-scale fading coefficients. In this way, the channel matrix between all terminals and BS antennas accounting for both large-scale fading and small-scale fading can be modeled as

$$\begin{aligned} {\mathbf{g}}=[{\mathbf{g}}_1,{\mathbf{g}}_2,\cdots ,{\mathbf{g}}_K]=\mathbf{HD}^{1/2}\in \mathbb{C}^{M\times K}, \end{aligned}$$

(1)

where ${\mathbf{h}}=[{\mathbf{h}}_1,{\mathbf{h}}_2,\cdots ,{\mathbf{h}}_K]\in \mathbb{C}^{M\times K}$ is small-scale fading matrix, and $\mathbf{D}={\rm{diag}}(\beta _1, \beta _2, \cdots , \beta _K)\in \mathbb{C}^{K\times K}$ is the large-scale fading diagonal matrix [24].

2.1.1 Uplink channel estimation

In practice, the BS needs CSI in order to take advantage of a large-scale antennas at BS in every frame. The typical of doing this is to utilize uplink pilots. In this phase, each terminal transmits an assigned uplink pilot sequence of length $p_l$ symbols, where $p_l\ge K$ is required to avoid pilot contamination. Clearly, $\tau T=p_lT_s$ should be satisfied, where $T_s$ is the sampling period. Denote that the kth terminal pilot sequence is $\varvec{\phi }_k$, for $k=1,2,\cdots \,K$, which is the kth column of $\varvec{\Phi }\in \mathbb{C}^{p_l\times K}$, satisfying $\varvec{\Phi }^H\varvec{\Phi }=\mathbf{I}_K$. The pilot signals propagate through the uplink channel. The received pilot signals at BS can be written as

$$\begin{aligned} \mathbf{Y}_p=\mathbf{GP}_p^{1/2}\varvec{\Phi }^H+\mathbf{N}_b \in \mathbb{C}^{M\times {p_l}}, \end{aligned}$$

(2)

where $\mathbf{P}_p=\rm{diag}(p_{p,1}, p_{p,2}, \cdots , p_{p,K})$ is the pilot transmit power diagonal matrix, and $\mathbf{N}_b$ is the additive white Gaussian noise (AWGN) matrix with i.i.d. elements following $\mathcal {CN}(0,\sigma ^2_b)$, which is introduced by BS antennas. According to (2), BS can apply minimum mean square error (MMSE) to obtain a channel estimate of ${\mathbf{g}}_k$ as follows:

$$\begin{aligned} \mathbf{Y}_p\varvec{\phi }_k=p_l\sqrt{p_{p,k}}{\mathbf{g}}_k+\mathbf{N}_b\varvec{\phi }_k. \end{aligned}$$

(3)

Thus, the MMSE estimate $\hat{{\mathbf{g}}}_k$ of the channel ${\mathbf{g}}_k$ is

$$\begin{aligned} \hat{{\mathbf{g}}}_k=\frac{\sqrt{p_{p,k}}\beta _k}{p_lp_{p,k}\beta _k+\sigma _b^2}\mathbf{Y}_p\varvec{\phi }_k \end{aligned}$$

(4)

and the estimation error is defined as

$$\begin{aligned} \tilde{\mathbf{g}}_k=\hat{{\mathbf{g}}}_k-{\mathbf{g}}_k. \end{aligned}$$

(5)

Consequently, according to [24, 25], the channel estimate and the estimation error are independent and distributed as

$$\begin{aligned} \hat{{\mathbf{g}}}_k\sim \mathcal {CN}(0,\gamma _k\mathbf{I}_M) \end{aligned}$$

(6)

and

$$\begin{aligned} \tilde{\mathbf{g}}_k\sim \mathcal {CN}(0,(\beta _k-\gamma _k)\mathbf{I}_M), \end{aligned}$$

(7)

where

$$\begin{aligned} \gamma _k=\frac{p_lp_{p,k}\beta _k^2}{p_lp_{p,k}\beta _k+\sigma _b^2}. \end{aligned}$$

(8)

By exploiting channels reciprocity, the downlink CSI can be obtained easily. Note that the uplink pilot transmit power of each terminal will be provided according to the energy allocated for uplink pilot transmission at that terminal. With the estimated channel, the BS can perform decoding on the uplink and precoding on the downlink.

2.1.2 Downlink WPT phase

In the downlink WPT phase, each terminal only harvests energy via TS protocol and does not demodulate the received signals. Therefore, all terminals can share the same constant symbol $x_e$ with $|x_e|=1$, which is known to all terminals. According to [17], the maximum ratio (MR) precoding is the optimal for energy transfer in the context of massive MIMO. Then the MR precoding vector at the BS for energy transfer can be written as

$$\begin{aligned} \mathbf{V}_k=\frac{\hat{{\mathbf{g}}}_k}{\sqrt{\mathbb{E}\{\Vert \hat{{\mathbf{g}}}_k\Vert ^2\}}}, \end{aligned}$$

(9)

where the scaling is used to satisfy the normalization constraint $\mathbb{E}\{\Vert \mathbf{v}_k\Vert ^2\}=1$. After the precoded signals are transmitted, the signals received by the kth terminal antenna can be given as

$$\begin{aligned} y^{ts}_k={\mathbf{g}}^H_k\sqrt{p^{\rm{d}l}_b\theta _k}\mathbf{v}_kx_e+{\mathbf{g}}^H_k\sum _{\begin{array}{c} j=1\\ j\ne k \end{array}}^K\sqrt{p_b^{\rm{d}l}\theta _j}\mathbf{v}_jx_e+n_{a,k}, \end{aligned}$$

(10)

where $p_b^{\rm{d}l}$ is the transmit power at the BS, $\theta _k(0\le \theta _k\le 1)$ is the PA coefficient assigned to the kth terminal, satisfying $\sum \limits _{k=1}^K\theta _k=1$, and $n_{a,k}\sim \mathcal {CN}(0,\sigma ^2_{a,k})$ is the AWGN introduced by the terminal antenna at the kth terminal.

2.1.3 Downlink SWIPT phase

In this phase, the BS broadcasts the information signals to all terminals simultaneously. Here, MR precoding is adopted again as it can approach the optimal beamforming solution in the context of massive MIMO [15]. According to PS protocol, a fraction of the received signal power is used for ID while the remaining fraction is used for EH. Denote $\rho _k (0\le \rho _k\le 1)$ as the PS ratio of the kth terminal, the signal split for ID and EH can be respectively expressed as

$$\begin{aligned} y^{ps}_{I,k}=\sqrt{\rho _k}\left({\mathbf{g}}^H_k\sqrt{p^{\rm{d}l}_b\theta _k}\mathbf{v} _kx_k^{\rm{d}l}+{\mathbf{g}}^H_k \mathop{\mathop{\sum}\limits_{j=1}}\limits_{j\ne k}^{K}\sqrt{p_b^{\rm{d}l}\theta _j}\mathbf{v}_jx_j^{\rm{d}l}+n_{a,k}\right)+n_{c,k} \end{aligned}$$

(11)

and

$$\begin{aligned} y^{ps}_{E,k}=\sqrt{1-\rho _k}\left({\mathbf{g}}^H_k\sqrt{p^{\rm{d}l}_b\theta _k}\mathbf{v}_kx_k^{\rm{d}l}+{\mathbf{g}}^H_k\mathop{\mathop{\sum}\limits_{j=1}}\limits_{j\ne k}^{K}\sqrt{p_b^{\rm{d}l}\theta _j}\mathbf{v}_jx_j^{\rm{d}l}+n_{a,k}\right), \end{aligned}$$

(12)

where $x_k^{\rm{d}l}\sim \mathcal {CN}(0,1)$ is the message-bearing downlink data symbol and $n_{c,k}\sim \mathcal {CN}(0,\sigma ^2_{c,k})$ is the additional AWGN introduced by RF to the baseband conversion at the kth terminal [26].

2.1.4 Uplink data transmission

In this phase, all terminals simultaneously transmit uplink data signals to the BS antennas. The signal vector received by the BS antennas is given by

$$\begin{aligned} \mathbf{y}_b=\sum _{k=1}^K{\mathbf{g}}_k\sqrt{p^{ul}_k}x^{ul}_k+\mathbf{n}_b, \end{aligned}$$

(13)

where $p^{ul}_k$ is the transmit power of the kth terminal, $x^{ul}_k\sim \mathcal {CN}(0,1)$, is the message-bearing uplink data symbol transmitted by the kth terminal, and $\mathbf{n}_b$ is the AWGN vector at the BS antennas whose elements follow $\mathcal {CN}(0,\sigma ^2_b)$. For the kth terminal, the BS processes the received signal through multiplication of the vector $\mathbf{y}_b\in \mathbb{C}^{M\times 1}$ by a decoding vector $\mathbf{a}_k\in \mathbb{C}^{M\times 1}$ that is a function of the channel estimate. The result of processing (13) is expressed as

$$\begin{aligned} y_{b,k}=\mathbf{a}_k^H\mathbf{y}_b=\mathbf{a}_k^H\sum _{j=1}^K{\mathbf{g}}_j\sqrt{p^{ul}_j}x^{ul}_j+\mathbf{a}_k^H\mathbf{n}_b. \end{aligned}$$

(14)

It is worth noting that the uplink data transmit power of each terminal will be obtained according to the energy allocated for uplink data transmission at that terminal.

2.2 Optimization problem formulation

In this section, the closed-form lower bound expressions on ergodic capacity for all terminals are first derived and then used to design a net sum SE maximization problem under some practical constraints.

2.2.1 Downlink ergodic capacity analysis

Since the CSI is not available at each terminal, similar to [27], (11) can be rewritten as

$$\begin{aligned} y_{I,k}^{ps}&=\sqrt{\rho _k}[\underbrace{\sqrt{p^{\rm{d}l}_b\theta _k}\mathbb{E}\{{\mathbf{g}}^H_k\mathbf{v}_k\}x_k^{\rm{d}l}}_{A1}+\nonumber \\&\underbrace{\sqrt{p^{\rm{d}l}_b\theta _k}({\mathbf{g}}^H_k\mathbf{v}_k-\mathbb{E}\{{\mathbf{g}}^H_k\mathbf{v}_k\})x_k^{\rm{d}l}}_{A2}+\sum _{\begin{array}{c} j=1\\ j\ne k \end{array}}^K\underbrace{\sqrt{p_b^{\rm{d}l}\theta _j}{\mathbf{g}}^H_k\mathbf{v}_jx_j^{\rm{d}l}}_{A3}+n_{a,k}]+n_{c,k}. \end{aligned}$$

(15)

Here, A1, A2, A3 represents the desired signal, the beamforming gain uncertainty and terminal interference, respectively. According to [28], the A2, A3 and the other terms are treated as uncorrelated noise. By assuming uncorrelated noise as independent Gaussian noise, a lower bound on ergodic capacity for the kth terminal can be expressed as

$$\begin{aligned} C_k^{\rm{d}l}=\log 2\left( 1+\frac{\mathbb{E}\{|\rm{DS}_k|^2\}}{\mathbb{E}\{|\rm{UN}_k|^2\}}\right) , \end{aligned}$$

(16)

where the desired signal power $\mathbb{E}\{|\rm{DS}_k|^2\}$ is calculated as

$$\begin{aligned} \mathbb{E}\{|\rm{DS}_k|^2\}=\rho _k|A1|^2=\rho _kp_b^{\rm{d}l}\theta _kM\gamma _k \end{aligned}$$

(17)

and the uncorrelated noise power $\mathbb{E}\{|\rm{UN}_k|^2\}$ is calculated as

$$\begin{aligned} \mathbb{E}\{|\rm{UN}_k|^2\}=\rho _k(|A2|^2+\sum _{\begin{array}{c} j=1\\ j\ne k \end{array}}^K|A3|^2+\sigma _{a,k}^2)+\sigma _{c,k}^2=\rho _k(p_b^{\rm{d}l}\beta _k+\sigma _{a,k}^2)+\sigma _{c,k}^2. \end{aligned}$$

(18)

The results in (17) and (18) can be obtained according to the proof for Theorem 1 in [29]. By substituting of (17) and (18) into (16), a lower bound on downlink ergodic capacity for the kth terminal can be recalculated as

$$\begin{aligned} C_k^{\rm{d}l}=\log 2\left( 1+\frac{\rho _kM\gamma _kp_b^{\rm{d}l}\theta _k}{\rho _k(\beta _kp_b^{\rm{d}l}+\sigma _{a,k}^2)+\sigma _{c,k}^2}\right) , \end{aligned}$$

(19)

where $\gamma _k$ will be calculated in the sequel.

2.2.2 Uplink ergodic capacity analysis

Before deriving uplink ergodic capacity for each terminal, we first calculate the total harvested energy in the WPT and SWIPT phases. Similar to (15), (10) and (12) can be respectively rewritten as

$$\begin{aligned} y^{ts}_k&=\mathbb{E}\{\sqrt{p^{\rm{d}l}_b\theta _k}{\mathbf{g}}^H_k\mathbf{v}_k\}x_e\nonumber \\&\quad +\sqrt{p^{\rm{d}l}_b\theta _k}({\mathbf{g}}^H_k\mathbf{v}_k-\mathbb{E}\{{\mathbf{g}}^H_k\mathbf{v}_k\})x_e+\sum _{\begin{array}{c} j=1\\ j\ne k \end{array}}^K\sqrt{p_b^{\rm{d}l}\theta _j}{\mathbf{g}}^H_k\mathbf{v}_jx_e+n_{a,k} \end{aligned}$$

(20)

and

$$\begin{aligned} y^{ps}_{E,k}&=\sqrt{1-\rho _k}[\mathbb{E}\{\sqrt{p^{\rm{d}l}_b\theta _k} {\mathbf{g}}^H_k\mathbf{v}_k\}x_k^{\rm{d}l}+\sqrt{p^{\rm{d}l}_b\theta _k} ({\mathbf{g}}^H_k\mathbf{v}_k-\mathbb{E}\{{\mathbf{g}}^H_k\mathbf{v}_k\})x_k^{\rm{d}l}\nonumber \\&\quad +\sum _{\begin{array}{c} j=1\\ j\ne k \end{array}}^K\sqrt{p_b^{\rm{d}l}\theta _j}{\mathbf{g}}^H_k\mathbf{v}_jx_j^{\rm{d}l}+n_{a,k}]. \end{aligned}$$

(21)

With the help of (17) and (18), the harvested energy by the k terminal in WPT and SWIPT phases can be respectively calculated as

$$\begin{aligned} E_k^{ts}=\alpha TM\gamma _kp_b^{\rm{d}l}\theta _k\eta _k+\alpha T(\beta _kp_b^{\rm{d}l}+\sigma _{a,k}^2)\eta _k, \end{aligned}$$

(22)

and

$$\begin{aligned} E_k^{ps}=\frac{1-\tau -\alpha }{2}T(1-\rho _k)M\gamma _kp_b^{\rm{d}l} \theta _k\eta _k+\frac{1-\tau -\alpha }{2}T(1-\rho _k)(\beta _kp_b^{\rm{d}l}+\sigma _{a,k}^2)\eta _k. \end{aligned}$$

(23)

where $\eta _k(0\le \eta _k\le 1)$ is the energy conversion efficiency of the EH circuits at the kth terminal.

Obviously, the first terms in the right-hand sides of (22) and (23) are proportional to the number of BS antennas while the second terms are independent with the number of BS antennas. In the context of massive MIMO, the first terms are dominate over the second terms. Similar to [17,18,19], the lower bound on total harvested energy is used to transmit uplink pilots and data and it is given as

$$\begin{aligned} E_k=\left[ \alpha +\frac{1-\tau -\alpha }{2}(1-\rho _k)\right] TM\gamma _kp_b^{\rm{d}l}\theta _k\eta _k. \end{aligned}$$

(24)

On average, each terminal can work normally all the time as sometimes it may consume its inherent power supply to compensate for power shortage [18]. We assume that a fraction $\lambda _k(0\le \lambda _k\le 1)$ of $E_k$ is used to transmit uplink pilots and the remaining $(1-\lambda )E_k$ energy is used to send uplink data at the kth terminal. Thus, the transmit power of the kth terminal for uplink pilot transmission is

$$\begin{aligned} p_lp_{p,k}=\frac{\lambda _kE_k}{\tau T}=\lambda _k\left[ \frac{\alpha }{\tau }+\frac{1-\tau -\alpha }{2\tau } (1-\rho _k)\right] M\gamma _kp_b^{\rm{d}l}\theta _k\eta _k \end{aligned}$$

(25)

and the transmit power of the kth terminal for uplink data transmission is

$$\begin{aligned} p_k^{ul}=\frac{(1-\lambda _k)E_k}{(1-\tau -\alpha )T/2}=(1-\lambda _k) \left[ \frac{2\alpha }{1-\tau -\alpha }+(1-\rho _k)\right] M\gamma _kp_b^{\rm{d}l}\theta _k\eta _k. \end{aligned}$$

(26)

By substituting (25) into (8), the channel estimate can be recalculated as

$$\begin{aligned} \gamma _k=\beta _k-\frac{\sigma _b^2}{\lambda _k \left[ \frac{\alpha }{\tau }+\frac{1-\tau -\alpha }{2\tau }(1-\rho _k)\right] Mp_b^{\rm{d}l}\theta _k\eta _k\beta _k} \end{aligned}$$

(27)

As the maximum-ratio combining (MRC) detection has lower computational complexity, which is compared with the zero-forcing (ZF) detection [17], the MRC detection is used to decode the received signals on the uplink. For the kth terminal, the MRC decoding vector is $\mathbf{a}_k=\hat{{\mathbf{g}}}_k$ and then substituting it to (14), we can obtain

$$\begin{aligned} y_{b,k}=\hat{{\mathbf{g}}}_k^H\sqrt{p_k^{ul}}{\mathbf{g}}_kx_k^{ul} +\hat{{\mathbf{g}}}_k^H\sum _{\begin{array}{c} j=1\\ j\ne k \end{array}}^K{\mathbf{g}}_j\sqrt{p^{ul}_j}x^{ul}_j+\hat{{\mathbf{g}}}_k^H\mathbf{n}_b. \end{aligned}$$

(28)

On the uplink, the CSI is known by BS and consequently the ergodic capacity expression involves inconvenient expectation outside the logarithm. To obtain a closed-form expression, an alternative lower bound expression is derived by utilizing the technique of “use and then forget CSI” in [30]. Thus, (28) is rewritten as

$$\begin{aligned} y_{{b,k}} = \sqrt {p_{k}^{{ul}} } {\mathbb{E}}\left\{ {\hat{{\mathbf{g}}}_{k}^{H} {\mathbf{g}}_{k} } \right\}x_{k}^{{ul}} + \sqrt {p_{k}^{{ul}} } \left( {\hat{{\mathbf{g}}}_{k}^{H} {\mathbf{g}}_{k} - {\mathbb{E}}\left\{ {\hat{{\mathbf{g}}}_{k}^{H} {\mathbf{g}}_{k} } \right\}} \right)x_{k}^{{ul}} + \hat{{\mathbf{g}}}_{k}^{H} \sum _{{\begin{array}{*{20}c} {j = 1} \\ {j \ne k} \\ \end{array} }}^{K} {\mathbf{g}}_{j} \sqrt {p_{j}^{{ul}} } x_{j}^{{ul}} + \hat{{\mathbf{g}}}_{k}^{H} {\mathbf{n}}_{b}. \end{aligned}$$

(29)

By following the similar derivation of (19), a closed-form lower bound expression on uplink ergodic capacity for the kth terminal can be given as

$$\begin{aligned} C_k^{ul}=\log 2\left( 1+\frac{M\gamma _kp_k^{ul}}{\sum \limits _{j=1}^K\beta _jp_j^{ul}+\sigma _b^2}\right) , \end{aligned}$$

(30)

where $p_k^{ul}$ is from (26) and $\gamma _k$ is from (27).

2.2.3 Net sum spectral efficiency maximization problem

As the factor of samples per frame that are used for transmission of uplink data and downlink data, respectively, is $(1-\tau -\alpha )/2$ shown as Fig. 2, the joint SE of uplink and downlink, namely, net sum SE, is

$$\begin{aligned} C_{ns}=\frac{1-\tau -\alpha }{2}\sum _{k=1}^K(C_k^{ul}+C_k^{\rm{d}l}). \end{aligned}$$

(31)

Denote $\varvec{\theta }$, $\varvec{\rho }$, and $\varvec{\lambda }$ as $K\times 1$ vectors that gather the elements of $\theta _k$, $\rho _k$ and $\lambda _k$, for $k=1\cdots K$, respectively. The net sum SE maximization problem can be formulated as

$$\begin{aligned} \mathcal {P}1:\quad&\max _{\begin{array}{c} \tau ,\alpha ,\\ \varvec{\theta },\varvec{\rho },\varvec{\lambda } \end{array}}C_{ns} \end{aligned}$$

(32a)

$$\begin{aligned} \rm{s.t.}\quad&\frac{1-\tau -\alpha }{2}C_k^{ul}>C_{\min }^{ul}\quad \forall k,\end{aligned}$$

(32b)

$$\begin{aligned}&\frac{1-\tau -\alpha }{2}C_k^{\rm{d}l}>C_{\min }^{\rm{d}l}\quad \forall k,\end{aligned}$$

(32c)

$$\begin{aligned}&0\le \theta _k\le 1,\sum _{k=1}^K\theta _k=1\quad \forall k,\end{aligned}$$

(32d)

$$\begin{aligned}&0\le \rho _k\le 1,\quad \forall k,\end{aligned}$$

(32e)

$$\begin{aligned}&0\le \lambda _k\le 1,\quad \forall k,\end{aligned}$$

(32f)

$$\begin{aligned}&KT_s\le \tau \le 1,\end{aligned}$$

(32g)

$$\begin{aligned}&0\le \alpha \le (1-\tau ). \end{aligned}$$

(32h)

In the above, (32b) and (32c) ensure the minimum data rate for each terminal on the uplink and on the downlink, respectively. (32d) specifies the constraint for each terminal and the sum constraint for all terminals on PA coefficient. (32e) and (32f) represent the constraints on PS ratio and energy allocation ratio for each terminal, respectively. (32g) and (32h) denote the constraints on the channel estimation duration and WPT duration for all terminals, respectively, where the lower bound $KT_s$ of $\tau$ is the minimum value to avoid pilot contamination.

Note that (32b) and (32c) can be satisfied when the transmit power of BS and the number of BS antennas are large enough. Thus, we assume that $\mathcal {P}1$ is feasible.

3 Methods

3.1 Proposed optimization algorithm

It can be observed from $\mathcal {P}1$ that the objective function and constraints are non-convex and non-linear. Moreover, the optimization variables are mutually coupled as follows:

• The WPT duration $\alpha$ and PS ratio $\varvec{\rho }$ affect the total harvested energy in WPT and SWIPT phases.

• The channel estimation duration $\tau$ and the fraction $\varvec{\lambda }$ of the total harvested energy for uplink pilot transmission affect CSI accuracy.

• The CSI accuracy affects the amount of total harvested energy and ergodic capacity for each terminal on the uplink and downlink.

• The CSI accuracy, the amount of total harvested and ergodic capacity for each terminal on the uplink and downlink are also closely related to the PA coefficients $\varvec{\theta }$.

Consequently, solving problem $\mathcal {P}1$ directly is a huge challenge. To provide a solution, we solve it in three steps. First, as all terminal share the same channel estimation duration $\tau$ and WPT duration $\alpha$, $\tau$ and $\alpha$ can be calculated numerically with one-dimensional search method in their respective feasible intervals. Then, for any given $\tau$ and $\alpha$, the optimal $\varvec{\theta }$, $\varvec{\rho }$, and $\varvec{\lambda }$ can be obtained by successive approximation method based on GP. Finally, when the search on $\tau$ and $\alpha$ ends, the optimal $\tau$ and $\alpha$ can be selected by comparison. In the following, we will discuss them in detail.

Since $\log _2(\cdot )$ is a monotonic increasing function, $\mathcal {P}1$ can be equivalently rewritten as $\mathcal {P}2$ when $\tau$ and $\alpha$ are given,

$$\begin{aligned} \mathcal {P}2:\quad&\min _{\begin{array}{c} \varvec{\theta },\varvec{\rho },\varvec{\lambda },\varvec{\gamma }\\ \varvec{\chi ^{ul}},\varvec{\chi ^{\rm{d}l}},\varvec{p^{ul}} \end{array}}\left[ \prod _{k=1}^K(1+\chi _k^{ul})(1+\chi _k^{\rm{d}l})\right] ^{-1} \end{aligned}$$

(33a)

$$\begin{aligned} \rm{s.t.}\quad&\chi _k^{ul}\le \frac{M\gamma _kp_k^{ul}}{\sum \limits _{j=1}^K\beta _jp_j^{ul}+\sigma _b^2}\quad \forall k,\end{aligned}$$

(33b)

$$\begin{aligned}&\chi _k^{\rm{d}l}\le \frac{\rho _kM\gamma _kp_b^{\rm{d}l}\theta _k}{\rho _k(\beta _kp_b^{\rm{d}l}+\sigma _{a,k}^2)+\sigma _{c,k}^2}\quad \forall k,\end{aligned}$$

(33c)

$$\begin{aligned}&\gamma _k\le \beta _k-\frac{\sigma _b^2}{\lambda _k[\frac{\alpha }{\tau }+\frac{1-\tau -\alpha }{2\tau }(1-\rho _k)]Mp_b^{\rm{d}l}\theta _k\eta _k\beta _k}\quad \forall k,\end{aligned}$$

(33d)

$$\begin{aligned}&p_k^{ul}\le (1-\lambda _k)[\frac{2\alpha }{1-\tau -\alpha }+(1-\rho _k)]M\gamma _kp_b^{\rm{d}l}\theta _k\eta _k\quad \forall k,\end{aligned}$$

(33e)

$$\begin{aligned}&(32b),(32c),(32d),(32e),(32f). \end{aligned}$$

(33f)

Here, $\chi _k^{ul}$ and $\chi _k^{\rm{d}l}$ are the power ratio of “desired signal” to “uncorrelated noise” on the uplink and downlink, respectively, and $\varvec{\gamma }$, $\varvec{\chi ^{ul}}$, $\varvec{\chi ^{\rm{d}l}}$ and $\varvec{p^{ul}}$ are denoted as $K\times 1$ vectors that gather the elements of $\gamma _k$, $\chi _k^{ul}$, $\chi _k^{\rm{d}l}$ and $p_k^{ul}$, for $k=1\cdots K$, respectively, which are auxiliary variables. By inspecting $\mathcal {P}2$, we observe the following cases:

• We have replaced “=” with “$\le$” in (33b), (33c), (33d) and (33e). However, this does not affect the original problem $\mathcal {P}2$ because the objective function is monotonic decreasing with respect to $\chi _k^{ul}$, $\chi _k^{\rm{d}l}$, $\gamma _k$ and $p_k^{ul}$ respectively when other variables hold constant.

• (32b) and (32c) can be transformed into monomial function by simple algebraic manipulations. (33b), (33c), (32d), (32e), and (32f) are monomial or posynomial functions.

• (33a), (33d), and (33e) are not monomial or posynomial functions.

If we transform (33a), (33d), and (33e) into monomial or posynomial functions respectively, the objective function and all constraints in $\mathcal {P}2$ are monomial or posynomial functions. As a result, $\mathcal {P}2$ becomes a GP problem, which can be solved efficiently with standard convex optimization tools.

For (33a), we transform it into a monomial function by an approximation method in [31, lemma 1]. The key idea is to utilize a monomial function $\varphi _k^{\lambda }(\chi _k^{\lambda })^{\phi _k^{\lambda }}$ to approximate $(1+\chi _k^{\lambda })$ near an arbitrary point $\hat{\chi }_k^{\lambda }>0$, where $\lambda \in \{ul,\rm{d}l\}$, $\phi _k^{\lambda }=\hat{\chi }_k^{\lambda }(1+\hat{\chi }_k^{\lambda })^{-1}$ and $\varphi _k^{\lambda }=(\hat{\chi }_k^{\lambda })^{-\phi _k^{\lambda }}(1+\hat{\chi }_k^{\lambda })$. The approximated result is

$$\begin{aligned} \prod _{k=1}^K(1+\chi _k^{ul})(1+\chi _k^{\rm{d}l})\approx \prod \limits _{k=1}^{K}\varphi _k^{ul}(\chi _k^{ul})^{\phi _k^{ul}}\varphi _k^{\rm{d}l}(\chi _k^{\rm{d}l})^{\phi _k^{\rm{d}l}}. \end{aligned}$$

(34)

In this way, the objective function is a monomial function.

For (33d), we deal with it as follows:

$$\begin{aligned}&\frac{\sigma _b^2}{\lambda _kt_kMp_b^{\rm{d}l}\theta _k\eta _k\beta _k}+\gamma _k\le \beta _k, \end{aligned}$$

(35)

$$\begin{aligned}&\frac{1-\tau -\alpha }{2\tau }\rho _k+t_k\le \frac{\alpha }{\tau }+\frac{1-\tau -\alpha }{2\tau }. \end{aligned}$$

(36)

Here, we introduce new variable $t_k$ and replace “=” with “$\le$”. As $\gamma _k$ is an increasing function of $t_k$ and the objective function is decreasing function of $\gamma _k$, the objective function is decreasing with respect to $t_k$. Consequently, the above operations do not change the original problem $\mathcal {P}2$. Obviously, (35) and (36) are posynomial functions.

Following the similar treatment of (33d), we can obtain equivalent expressions for (33e) as follows:

$$\begin{aligned}&\rho _k+(\bar{\lambda }_kM\gamma _kp_b^{\rm{d}l}\theta _k\eta _k)^{-1}p_k^{ul}\le \frac{1-\tau +\alpha }{1-\tau -\alpha }, \end{aligned}$$

(37)

$$\begin{aligned}&\lambda _k+\bar{\lambda }_k\le 1. \end{aligned}$$

(38)

Here, $\bar{\lambda }_k$ is the new introduced variable. Clearly, (37) and (38) are posynomial functions.

With the treatment above, $\mathcal {P}2$ is transformed into a GP problem. For subsequent convenience, we denote $\varvec{t}$, $\varvec{\lambda }$ and $\varvec{\bar{\lambda }}$ as $K\times 1$ vectors that gather the elements of $t_k$, $\lambda _k$ and $\bar{\lambda }_k$ for $k=1,\cdots ,K$, respectively. According to [31,32,33], a successive approximation algorithm based on GP to solve $\mathcal {P}2$ is proposed in Algorithm 1.

The parameter $\mu$ is used to control approximation accuracy and is set to be 1.1 in most practical cases [31] and the convergence of Algorithm 1 is also guaranteed [34].

After solving out the problem $\mathcal {P}2$, the optimal $\tau$ and $\alpha$ can be obtained by one-dimensional search method, which is described in Algorithm 2.

3.2 Computational complexity

Finally, we discuss the complexity of Algorithm 2, which is mainly dependent on the complexities of two outer loops and inner Algorithm 1. If assume that $\Delta _{\tau }$ and $\Delta _{\alpha }$ are the search step sizes for $\tau$ and $\alpha$, respectively, the complexity of two outer loops approximately is $\mathcal {O}((1-KT_s)(1-KT_s+\Delta _{\tau })/(2\Delta _{\alpha }\Delta _{\tau }))$. On the other hand, the complexity of Algorithm 1 can be approximately expressed as $\mathcal {O}(N_{ap}KN_{gp}/\varepsilon ^2)$, where $N_{gp}$ is the number of required iteration for solving GP and $N_{ap}$ is the number of required iteration for successive approximation. Based on the above analyses, the complexity of Algorithm 2 is approximately given as $\mathcal {O}((1-KT_s)(1-KT_s+\Delta _{\tau })N_{ap}KN_{gp}/(2\Delta _{\alpha }\Delta _{\tau }\varepsilon ^2))$. As the computational efficiency of GP modeling is very high even for large-scale problems [31], Algorithm 2 can converge quickly once the step size of one-dimensional search and maximum tolerance are carefully chosen, which is also shown by numerical results in Sect. 4.

4 Results and discussion

In this section, numerical results based on matlab software are conducted to validate the proposed algorithm. We set $p_b^{\rm{d}l}=1.5$Watt $M=300$ and $K=3$. Each frame length T is normalized to be 1 and the sampling period $T_s$ is assumed to be 0.005. The noise power is set to be $\sigma _b^2=-90$dBm, $\sigma _{a,k}^2=-70$dBm and $\sigma _{c,k}^2=-50$dBm, respectively, $\forall k$, [35, 36]. The energy conversion efficiency is set to be $\eta _k=80\%$, $\forall k$. The minimum data requirement is set to be $C_{\min }^{ul}=1$bit/s/Hz and $C_{\min }^{\rm{d}l}=2$bit/s/Hz, respectively. The large-scale fading is modeled as $\beta _k=10^{-3}d_k^{-3}$, where the distance away from BS $d_1=10$m, $d_2=15$m and $d_3=20$m. The tolerance $\varepsilon$ for Algorithm 2 is chosen as $10^{-4}$. The above parameters are used throughout the simulations unless otherwise stated.

Fig. 3

The running process of Algorithm 1 for given $\tau$ and $\alpha$

Figure 3 shows the running process of Algorithm 1 when channel estimation duration $\tau$ is set to be 0.015 and WPT duration $\alpha$ is set to be 0.005. It can be seen that Algorithm 1 achieves the maximum net sum SE $C_{ns}=18.8$ bit/s/Hz with 31 iterations, which means that Algorithm 1 is able to converge quickly to match the channel condition. Meanwhile, we can obtain the optimal PA coefficient $\varvec{\theta }=[0.28, 0.29, 0.43]^T$, PS ratio $\varvec{\rho }=[0.76, 0.49, 0.24]^T$ and energy allocation ratio $\varvec{\lambda }=[0.91, 0.53, 0.27]^T$, respectively. Thus, we can conclude that Algorithm 1 solves successfully the optimization problem $\mathcal {P}2$. In addition, the optimal $\tau$ and $\alpha$ can be obtained by one-dimensional search method in Algorithm 2. To sum up, the original optimization problem $\mathcal {P}1$ can be effectively solved by Algorithm 2, which also indicates that the proposed algorithms are feasible and effective.

Fig. 4

Net sum spectral efficiency versus channel estimation duration

Fig. 5

Net sum spectral efficiency of uplink and downlink versus channel estimation duration

Although spending more time on channel estimation results in more accurate CSI, which can improve the net sum SE, this also leads to reduced duration for data transmission on the uplink and downlink, which degrades the net sum SE. Thus, there exists an optimal value for channel estimation duration. Figure 4 depicts the net sum SE versus the channel estimation duration, to show its impact on the system performance under different transmission models. Here, only uplink transmission model or only downlink transmission model can be obtained by setting $\varvec{\rho }=[0, 0,0]^T$ or $\varvec{\lambda }=[1, 1, 1]^T$, respectively. The step size for one-dimensional search is set to be 0.005. It is observed from Fig. 4 that the net sum SE decreases with channel estimation duration $\tau$ for three transmission models. Although we can find a very slight rise of the net sum SE when the step size of one-dimensional search is set to be small enough, it is tiny and can be negligible. For example, the improvement of the net sum SE is only the level of $10^{-11}$ when the step size for one-dimensional search is set to be $2\times 10^{-14}$, as shown in Fig. 5. This means that the improvement of the net sum SE from longer channel estimation duration can not compensate for its decline from reduced duration for data transmission. This is due to the fact that the channel of each terminal fluctuates only slightly around its expected value because of the channel harding phenomenon in massive MIMO systems, and hence the longer channel estimation duration has little effect on improving CSI. Therefore, the optimal channel estimation duration $\tau$ tends to its minimum value. In addition, as there exist the minimum data rate requirements of each terminal on the uplink and downlink, the value of $\tau$, which makes Algorithm 2 work normally, can not cover the entire feasible interval.

Fig. 6

Net sum spectral efficiency versus the transmit power of base station

Figure 6 captures the impact of PS, TS and hybrid EH protocols on the net sum SE under various values of the transmit power of BS. Here, the PA coefficient is set to be $\varvec{\theta }=[1, 1, 1]^T/3$ and the energy allocation ratio is set to be $\varvec{\lambda }=[0.5, 0.5, 0.5]^T$. The channel estimation duration $\tau$ is set to be 0.015 and the search step size for $\alpha$ is set to be 0.005. It is observed from Fig. 6 that the net sum SE based on TS protocol is always lower than that based on PS and hybrid protocols under the same transmit power of BS. Moreover, the performance of PS protocol is nearly the same as that of hybrid protocol. When the search step size of $\alpha$ is set to be small enough, it can be found that the hybrid protocol outperforms weakly PS protocol, but the difference between them is tiny and thus can be negligible. This indicates that although spending more time on EH results in more energy for channel estimation and uplink transmission, which can improve the net sum SE, it can not compensate for the decline of the net sum SE due to reduced duration for data transmission on the uplink and downlink. Therefore, the optimal WPT duration $\alpha$ tends to zero. As a result, the hybrid protocol is very close to PS protocol for the net sum SE maximization in massive MIMO systems.

Fig. 7

Net SE of each terminal versus PS ratio, where $\tau =0.015$, $\alpha =0$, $C_{\min }^{ul}=C_{\min }^{\rm{d}l}=0$, $\varvec{\theta }=[1, 1, 1]^T/3$, and $\varvec{\lambda }=[0.5, 0.5, 0.5]^T$

Fig. 8

Net SE of each terminal versus energy allocation ratio, where $\tau =0.015$, $\alpha =0$, $C_{\min }^{ul}=C_{\min }^{\rm{d}l}=0$, $\varvec{\theta }=[1, 1, 1]^T/3$, and $\varvec{\rho }=[0.5, 0.5, 0.5]^T$

We capture the impact of PS ratio $\varvec{\rho }$ and energy allocation ratio $\varvec{\lambda }$ on the net SE of each terminal in Figs. 7 and 8, respectively. We focus on the impact of $\varvec{\rho }$ solely in Fig. 7 and it is seen that the net SE of each terminal is a quasi-concave function with respect to its PS ratio. This is due to the fact that the higher PS ratio increases the downlink SE of each terminal, but at the same time decreases the harvested energy for channel estimation and uplink transmission, which degrades the uplink SE of each terminal. Thus there exists an optimal PS ratio for each terminal to maximize its net SE. Similarly, we focus on the impact of $\varvec{\lambda }$ solely in Fig. 8 and it is observed that the net SE of each terminal is a quasi-concave function with respect to its energy allocation ratio. This is because that the higher energy allocation ratio results in more accurate CSI, which can improve the net SE of each terminal, but meanwhile decreases the transmit power of each terminal on the uplink, which degrades the uplink SE of each terminal. Thus there also exists an optimal energy allocation ratio for each terminal to maximize its net SE. The above analyses indicate the feasibility of the optimization problem $\mathcal {P}1$.

Fig. 9

Net SE of each terminal versus the transmit power of BS, where $\tau =0.015$, $\alpha =0$ and $\varvec{\lambda }=[0.5, 0.5, 0.5]^T$

Figure 9 illustrates the the impact of PA coefficient $\varvec{\theta }$ on the net sum SE under various values of the transmit power of BS. For comparison, the equal PA method is considered as a baseline scheme in which the BS equally allocates transmit power to each terminal it serves. The PA obtained by Algorithm 1 is denoted as “optimal PA” in Fig. 9. It is observed that the optimal PA method achieves a higher net sum SE than the equal PA method, which shows the optimality of the propose algorithm. Moreover, we can also see that the transmit power of BS and the number of BS antennas are efficient way to improve the net sum SE.

Fig. 10

Net SE of uplink and downlink data transmission versus the number of BS antennas

As the maximal minimum rate algorithm is common in the existing literatures, in Fig. 10, the performance of the proposed algorithm is compared with that of the maximal minimum rate algorithm for uplink and downlink data transmission, respectively. On the uplink, the asymptotically maximal minimum rate is obtained according to the analytical result derived in [17]. On the downlink, the maximal minimum rate is obtained according to the algorithm proposed in [18]. For interpreting the results easily, the maximum and minimum values of the proposed algorithm for uplink and downlink data transmission are used respectively. We denote the maximal minimum rate by “MMR” and the proposed algorithm by “SRM” in Fig. 10. It is observed that MMR is located between SRM-max and SRM-min for uplink and downlink, respectively. This is due to the fact that the MMR algorithm has to compensate the performance of the terminal with minimum rate with that of the terminal with maximum rate, indicating that MMR is a special case of the proposed algorithm with the minimum rate constraint.

5 Conclusions

The paper has proposed a joint uplink and downlink transmission scheme in massive MIMO systems and it can conveniently switch to only uplink transmission model or only downlink transmission model. Each terminal uses the hybrid wireless EH protocol to harvest energy and it can run in three modes, i.e., TS, PS or hybrid protocol. According to the derived low bound expressions on ergodic capacity for all terminals on the uplink and downlink, we have investigated how to jointly optimize related system parameters to maximize the net sum SE of the whole system. As the formulated problem is the non-convex and non-linear, an algorithm utilizing one-dimensional search method and successive approximation method based GP has been proposed and its convergence and complexity have also been discussed. Finally, numerical results have manifested the feasibility of the proposed algorithm. Nevertheless, there are still some issues related to practical scenarios, such as BS antenna correlation, nonlinearity of energy collection model, fading channel model with line-of-sight, and they are left for future work.