Robust secure design for relay wireless sensor networks with simultaneous wireless information and power transfer

Abstract

In this article, we investigate the physical layer security problem in relay wireless sensor networks with simultaneous wireless information and power transfer. Specifically, by considering the energy receiver which may act as potential eavesdropper and assuming that only imperfect channel state information can be attained, we propose a joint robust beamforming and artificial noise scheme. We formulate the relay power minimization problem under both the secrecy rate constraint and the energy harvesting constraints, which is non-convex and hard to tackle. By studying the hidden convexity of the problem, we propose an iterative algorithm combining semi-definite relaxation and a gradient-based method. In addition, a suboptimal rank-one reconstruction method is presented. Simulation results show that the proposed robust scheme achieves better performance than other schemes.

Introduction

Wireless sensor network (WSN) is envisioned as one of the most important research areas with numerous applications such as vehicular tracking, environmental monitoring, ¹ and other fields. ² A WSN is typically composed of energy-constrained, cost-constrained, and lightweight computing sensor nodes, which can perform sensing, simple computations, and short-range wireless communications. ³

Energy harvesting (EH) is considered as one of the effective approaches for improving the energy efficiency of WSN since in many WSN applications, the lifetime of network is limited due to the constraints in energy resources. ⁴ EH-enabled sensors can harvest energy from either radio frequency (RF) signals or ambient energy sources which enable them to operate continuously. ⁵ Simultaneous wireless information and power transfer (SWIPT) is a promising EH technique to solve the energy scarcity problem in energy-constrained wireless system such as WSN. ⁶

However, the secure transmission has been recognized as an important issue for WSN due to the broadcast characteristics of wireless channel and the security requirement of the sensing information. ⁷ To address the issue, the physical layer security (PLS) technique, ⁸ which exploits the characteristics of wireless channels such as fading to fulfill secure communication, has been proved to be an effective way to improve security in WSN ⁹ as well as SWIPT system.^10–12 Specifically, Curiac ⁹ investigated the utility of directional antennas to improve security in WSN, while Xing et al., ¹⁰ Liu et al., ¹¹ and Shi et al. ¹² investigated the secrecy SWIPT scheme for single-input single-output (SISO) channel, multiple-input single-output (MISO) channel, and multiple-input multiple-output (MIMO) channel, respectively. Among these literatures, beamforming (BF) and artificial noise (AN) are two key techniques to improve security and enhance EH. ¹³

One of the most important problems in PLS is the acquirement of the channel state information (CSI). The robust design has been widely applied to handle the CSI’s uncertainties.^14–20 Specifically, Feng et al., ¹⁴ Chu et al., ¹⁵ and Zhang et al. ¹⁶ investigated the robust secrecy rate maximization (SRM) problems for an MISO SWIPT system without AN, with AN, and with a friendly jammer, respectively. Khandaker and Wong ¹⁷ investigated the robust power minimization and the robust SRM problem in MISO SWIPT by considering that the EH receiver may be potential eavesdropper (Eve), while Wang and Wang ¹⁸ investigated the robust SRM problem in MIMO wiretap channels. Recently, Chu et al. ¹⁹ and Zhang et al. ²⁰ investigated the robust secrecy design for MISO SWIPT system by considering that all receivers have EH constraints.

Relaying is considered as a popular approach to extend network coverage and provide spatial degrees of freedom (DoF), which is beneficial to PLS ²¹ and SWIPT; ²² thus, several literatures investigated the secure SWIPT in relay system.^23–27 Specifically, Li et al. ²³ investigated the SRM problem in multi-antenna amplify-and-forward (AF) relay SWIPT system subject to the EH constraints, while Li et al. ²⁴ considered unavailable Eve’s CSI and proposed a suboptimal null space–based BF and AN scheme in two-way AF SWIPT relay networks. Zhang et al. ²⁵ proposed a joint design of the signal BF and AN at the source and the BF at the relay to maximize the secrecy rate, while Xing et al. ²⁶ proposed a novel harvest-and-jamming method, where the friendly jammer used the harvested energy to confuse the Eve. Recently, Feng et al. ²⁷ proposed a joint BF and energy signal scheme for providing both secure communication and energy transfer in AF relay networks.

Relay node can make the WSN transmission more reliable to satisfy the strict requirements in industrial applications. ²⁸ Recently, relay WSN has aroused great attentions.^29–35 Specifically, Dong et al. ²⁹ investigated the diversity performance of high-altitude platform (HAP)-based relay WSN, while Dong et al. ³⁰ investigated the energy-efficient transmission for remote WSN, where an integrated architecture based on HAP was proposed. Djenouri and Bagaa ³¹ investigated the relay node deployment for sustainable WSN with EH constraint. For the issue of PLS in relay WSN, Liau et al. ³² proposed a two-path successive relaying scheme to improve the security, while Gong et al. ³³ utilized joint BF and AN to further improve the secrecy performance. For more complicated network structure, Zhang et al. ³⁴ investigated the PLS scheme for orthogonal frequency division multiple access (OFDMA)-based two-way relay WSN, while Deng et al. ³⁵ proposed a PLS scheme for three-tier WSN using decode-and-forward (DF) relay.

Motivated by these works, in this article, we investigate the PLS for SWIPT in relay WSN. Specifically, we investigate the BF and AN design to minimize the relay transmit power based on the following settings: (1) one multi-antenna relay forward the information from the source to the information sensor node, denoted as information receiver (IR), in the presence of multiple energy sensor nodes, denoted as energy receivers (ERs), which may be potential Eve; (2) relay employs BF and AN to improve security and enhance energy transfer simultaneously; (3) relay has perfect CSI of the IR, but imperfect CSI of the ERs and the ERs colluding to eavesdrop the confidential information.

Our main contributions can be summarized as follows: first, we propose a framework for secrecy SWIPT in AF relay WSN with different service requirements at different wireless sensor nodes which gives rise to new BF and AN task. Second, the power minimization problem subject to the non-convex secrecy rate constraint is hard to handle directly. Our gradient-based algorithm solves this problem efficiently and can be applied to similar scenarios. Third, by analyzing the Karush–Kuhn–Tucker (KKT) optimality conditions, we discuss the structure of the obtained BF matrix and propose a suboptimal rank-one reconstruction algorithm when the rank of the matrix is higher than one. To the best of our knowledge, such a joint robust design for secure SWIPT in AF relay WSN has not been investigated in related literatures.

One of the most relevant works is Li et al.’s study, ²¹ in which the authors investigated the robust secrecy scheme in multi-antenna relay networks, but the main differences of our article are summarized as follows: first, Li et al. ²¹ did not consider SWIPT, while our article considers the combination of PLS and SWIPT; thus, the BF and AN in our system have dual functions, for example, enhance the received signals at the IR and the ERs while degrade the received signals at the Eves simultaneously. Second, Li et al. ²¹ considered the SRM problem, while our article investigates the power minimization problem; although the one-dimensional search method in the work by Li et al. ²¹ can solve our problem with some necessary modifications, we propose a more efficient gradient-based algorithm by analyzing the hidden convexity of the power minimization problem. Third, the algorithm in the work by Li et al. ²¹ can always obtain a rank-one optimal BF matrix, while in the case of SWIPT, the structure of the obtained matrix is more complicated and we cannot always attain a rank-one BF matrix, so we propose an effective way to obtain a suboptimal solution when we obtain a high-rank solution.

The rest of this article is organized as follows. A system model description and problem statement is given in section “System model and problem statement.” Section “Joint robust BF and AN design for equation (8)” considers the joint robust BF and AN design problem, wherein a gradient-based optimization approach is established. Simulation results are illustrated in section “Simulation result.” Section “Conclusion” concludes the article.

Notation

Throughout the article, we use the upper case boldface letters for matrices and lower case boldface letters for vectors. The superscripts $(·{)}^T$ , $(·{)}^{†}$ , and $(·{)}^H$ represent the transpose, conjugate, and conjugate transpose, respectively. The trace and rank of matrix $\mathbf{A}$ are denoted as $\mathrm{Tr}(\mathbf{A})$ and $\mathrm{rank}(\mathbf{A})$ , respectively. $\mathbf{a}=\mathrm{vec}(\mathbf{A})$ denotes to stack the columns of matrix $\mathbf{A}$ into a vector $\mathbf{a}$ . $\mathrm{ve}{\mathrm{c}}^{-1}()$ is the inverse operation of $\mathrm{vec}()$ . $\mathbf{A}⪰0$ indicates that $\mathbf{A}$ is a positive semi-definite matrix. |·| and $\parallel ·{\parallel }_F$ represent the absolute value and the Frobenius norm, respectively. $\mathrm{null}(\mathbf{A})$ denotes the null space of $\mathbf{A}$ . ⊗ denotes Kronecker product. $\mathbf{D}(\mathbf{a})$ represents a diagonal matrix with $\mathbf{a}$ on the main diagonal, while $\mathbf{D}({\mathbf{A}}_1,\dots ,{\mathbf{A}}_N)$ represents a block diagonal matrix with its main diagonal blocks being ${\mathbf{A}}_1,\dots ,{\mathbf{A}}_N$ . $\mathbf{I}$ denotes an identity matrix with appropriate size. $\mathrm{Re}\{a\}$ denotes the real part of a complex variable $a$ . $\mathcal{C}N(0,\mathbf{I})$ denotes a circularly symmetric complex Gaussian random vector with mean $0$ and covariance matrix $\mathbf{I}$ . $\mathbb{E}[.]$ stands for the statistical expectation.

System model and problem statement

We consider a relay WSN system for SWIPT as shown in Figure 1, in which a transmitter transmits information and wireless energy to a cluster of sensor nodes simultaneously with the help of one multi-antenna AF relay. One of these sensor nodes is IR, while the other sensor nodes are ERs. The relay is equipped with $N_r$ antennas, while the other nodes are equipped with single antenna. Let $\mathbf{f}\in {\mathbb{C}}^{N_r\times 1}$ , $\mathbf{h}\in {\mathbb{C}}^{N_r\times 1}$ , and ${\mathbf{g}}_m\in {\mathbb{C}}^{N_r\times 1}$ denote the channel responses from the transmitter to the relay, the relay to the IR, and the relay to the mth ER, respectively. We assume that the channels between the transmitter and the relay are perfectly known and there is no direct link between the transmitter to the IR or the ERs, which is a common assumption.^23,25 Since the relay operates in a half-duplex mode, one transmission round is composed of two phases.

Figure 1.

System model for secure AF relay WSN with SWIPT.

In the first phase, the transmitter broadcasts its information $s$ satisfying $\mathbb{E}[|s{|}^2]=1$ to the relay. The received signal at the relay is given by

$$y_r=\sqrt{P_s}{f}_s+n_r$$

(1)

where $P_s$ is the transmit power at the transmitter and ${\mathbf{n}}_r$ is the additive noise at the relay with variance ${\sigma }_r^2$ .

In the second phase, the relay employs the BF matrix $\mathbf{W}\in {\mathbb{C}}^{N_r\times N_r}$ to forward $s$ and emit AN vector $\mathbf{z}\in {\mathbb{C}}^{N_r\times 1}$ which follows $\mathcal{CN}(0,\mathbf{Q})$ to help secure communication as well as energy transmission. Thus, the signal transmitted by the relay is given by

$$x_r=W{y}_r+z=\sqrt{P_s}Wfs+W{n}_r+z$$

(2)

The received signals at the IR and the mth ER are given, respectively, by

$$y_I=\sqrt{P_s}{\mathbf{h}}^H\mathbf{Wf}s+{\mathbf{h}}^H\mathbf{z}+{\mathbf{h}}^H\mathbf{W}{\mathbf{n}}_r+n_I$$

(3a)

$$y_{E,m}=\sqrt{P_s}{\mathbf{g}}_m^H\mathbf{Wf}s+{\mathbf{g}}_m^H\mathbf{z}+{\mathbf{g}}_m^H\mathbf{W}{\mathbf{n}}_r+n_E$$

(3b)

where $n_I$ and $n_E$ are additive noises at IR and the mth ER, with variances ${\sigma }_I^2$ and ${\sigma }_E^2$ , respectively.

In this article, we consider the worst-case secrecy design. Specifically, we assume that all these single antenna ERs perform joint processing to eavesdrop the confidential information; thus, they can be seen as a virtual multi-antenna ER^17,27 and the relay only has imperfect ER’s CSI. We use the deterministic spherical model^26,27 to characterize the CSI uncertainties such that

$$\mathcal{G}=\{{\mathbf{G}}_E|{\mathbf{G}}_E={\overline{\mathbf{G}}}_E+\Delta {\mathbf{G}}_E,\mathrm{Tr}\left(\Delta {\mathbf{G}}_E^H\mathbf{\Lambda }\Delta {\mathbf{G}}_E\right)\le 1\}$$

(4)

where ${\overline{\mathbf{G}}}_E$ denotes the estimate of the colluding eavesdrop channel ${\mathbf{G}}_E$ with ${\mathbf{G}}_E\stackrel{\Delta }{=}[{\mathbf{g}}_1,\dots ,{\mathbf{g}}_m]$ , $\Delta {\mathbf{G}}_E$ denotes the channel uncertainties, and $\mathbf{\Lambda }$ determines the quality of the respective CSI and assumed to be known. In addition, we define ${\mathbf{g}}_E\stackrel{\Delta }{=}\mathrm{vec}({\mathbf{G}}_E)$ , ${\overline{\mathbf{g}}}_E\stackrel{\Delta }{=}\mathrm{vec}({\overline{\mathbf{G}}}_E)$ , and $\Delta {\mathbf{g}}_E\stackrel{\Delta }{=}\mathrm{vec}(\Delta {\mathbf{G}}_E)$ for the convenient of the following statement.

Thus, the information rates at the IR and the colluding ERs are given, respectively, by

$$\begin{array}{c}{C}_I(\mathbf{W},\mathbf{Q})=\underset{2}{\log }(1+\frac{P_s{\left|{\mathbf{h}}^H\mathbf{Wf}\right|}^2}{{\sigma }_I^2+{\sigma }_r^2{\mathbf{h}}^H\mathbf{W}{\mathbf{W}}^H\mathbf{h}+{\mathbf{h}}^H\mathbf{Qh}})\hfill \end{array}$$

(5a)

$$C_E(\mathbf{W},\mathbf{Q})=\underset{2}{\log }(1+P_s{\left({\mathbf{G}}_E^H\mathbf{Wf}\right)}^H{\mathbf{R}}_E^{-1}\left({\mathbf{G}}_E^H\mathbf{Wf}\right))$$

(5b)

where

$${\mathbf{R}}_E={\sigma }_E^2\mathbf{I}+{\sigma }_r^2{\mathbf{G}}_E^H\mathbf{W}{\mathbf{W}}^H{\mathbf{G}}_E+{\mathbf{G}}_E^H\mathbf{Q}{\mathbf{G}}_E$$

(6)

and the scaling factor 1/2 is due to the half-duplex operation of relay (See Feng et al. ²⁷ for more details).

However, by neglecting the noise power, the harvest energy of the mth ER can be expressed as

$$\begin{array}{c}{E}_m(\mathbf{W},\mathbf{Q})=\hfill \\ \multicolumn{1}{c}{{\rho }_m{\mathbf{e}}_m^H{\mathbf{G}}_E^H(P_s\mathbf{Wf}{\mathbf{f}}^H{\mathbf{W}}^H+{\sigma }_r^2\mathbf{W}{\mathbf{W}}^H+\mathbf{Q}){\mathbf{G}}_E{\mathbf{e}}_m}\end{array}$$

(7)

where $0\le {\rho }_m\le 1$ stands for the EH efficiency for the mth ER and ${\mathbf{e}}_m$ is the $M\times 1$ unit vector with the mth entry being equal to one.

Our problem of interest is to design the BF matrix $\mathbf{W}$ and AN covariance $\mathbf{Q}$ such that minimize the relay transmit power while meeting the secrecy rate at the IR and the EH constraint at the ERs. Mathematically, our problem can be expressed as

$$\underset{\mathbf{W},\mathbf{Q}⪰0}{min}\mathrm{Tr}(\mathbf{WC}{\mathbf{W}}^H+\mathbf{Q})$$

(8a)

$$\mathrm{s}.\mathrm{t}.\frac{1}{2}\underset{\forall {\mathbf{G}}_E\in \mathcal{G}}{min}\{C_I(\mathbf{W},\mathbf{Q})-C_E(\mathbf{W},\mathbf{Q})\}\ge R_s$$

(8b)

$$E_m(\mathbf{W},\mathbf{Q})\ge {\eta }_m,\forall {\mathbf{G}}_E\in \mathcal{G},\forall m$$

(8c)

where $\mathbf{C}=P_s\mathbf{f}{\mathbf{f}}^H+{\sigma }_r^2\mathbf{I}$ denotes the relay amplify matrix, $R_s$ denotes the target secrecy rate, and ${\eta }_m$ denotes the EH threshold for the mth ER.

Joint robust BF and AN design for equation (8)

Equivalent reformulation of equation (8)

We start by reformulating equation (8) into

$$\underset{\mathbf{W},\mathbf{Q}⪰0}{min}\mathrm{Tr}(\mathbf{WC}{\mathbf{W}}^H+\mathbf{Q})$$

(9a)

$$\mathrm{s}.\mathrm{t}.\frac{1+\frac{P_s{\left|\mathbf{hWf}\right|}^2}{{\sigma }_I^2+{\mathbf{h}}^H\mathbf{W}{\mathbf{W}}^H\mathbf{h}+{\mathbf{h}}^H\mathbf{Qh}}}{1+P_s{\left({\mathbf{G}}_E^H\mathbf{Wf}\right)}^H{\mathbf{R}}_E^{-1}\left({\mathbf{G}}_E^H\mathbf{Wf}\right)}\ge 2^{2R_s},\forall {\mathbf{G}}_E\in \mathcal{G}$$

(9b)

$$E_m(\mathbf{W},\mathbf{Q})\ge {\eta }_m,\forall {\mathbf{G}}_E\in \mathcal{G},\forall m$$

(9c)

Then, by introducing an auxiliary variable $\gamma$ , equation (9) can be further rewritten as

$$\underset{\mathbf{W},\mathbf{Q}⪰0}{min}\mathrm{Tr}(\mathbf{WC}{\mathbf{W}}^H+\mathbf{Q})$$

(10a)

$$\begin{array}{c}\mathrm{s}.\mathrm{t}.1+\frac{P_s{\left|{\mathbf{h}}^H\mathbf{Wf}\right|}^2}{{\sigma }_I^2+{\mathbf{h}}^H\mathbf{W}{\mathbf{W}}^H\mathbf{h}+{\mathbf{h}}^H\mathbf{Qh}}\ge \gamma 2^{2R_s}\hfill \end{array}$$

(10b)

$$1+P_s{\left({\mathbf{G}}_E^H\mathbf{Wf}\right)}^H{\mathbf{R}}_E^{-1}\left({\mathbf{G}}_E^H\mathbf{Wf}\right)\le \gamma ,\forall {\mathbf{G}}_E\in \mathcal{G}$$

(10c)

$$E_m(\mathbf{W},\mathbf{Q})\ge {\eta }_m,\forall {\mathbf{G}}_E\in \mathcal{G},\forall m$$

(10d)

which can be re-organized as

$$\underset{\mathbf{W},\mathbf{Q}⪰0}{min}\mathrm{Tr}(\mathbf{WC}{\mathbf{W}}^H+\mathbf{Q})$$

(11a)

$$\begin{array}{c}\mathrm{s}.\mathrm{t}.P_s{\left|{\mathbf{h}}^H\mathbf{Wf}\right|}^2\hfill \\ \multicolumn{1}{c}{\ge (\gamma 2^{2R_s}-1)({\sigma }_I^2+{\mathbf{h}}^H\mathbf{W}{\mathbf{W}}^H\mathbf{h}+{\mathbf{h}}^H\mathbf{Qh})}\end{array}$$

(11b)

$$P_s{\left({\mathbf{G}}_E^H\mathbf{Wf}\right)}^H{\mathbf{R}}_E^{-1}\left({\mathbf{G}}_E^H\mathbf{Wf}\right)\le \gamma -1,\forall {\mathbf{G}}_E\in \mathcal{G}$$

(11c)

$$E_m(\mathbf{W},\mathbf{Q})\ge {\eta }_m,\forall {\mathbf{G}}_E\in \mathcal{G},\forall m$$

(11d)

The equivalence of equations (8) and (11) can be shortly demonstrated as follows. Let $\{\mathbf{W},\mathbf{Q}\}$ be a feasible solution to equation (8); if we set $\gamma \stackrel{\Delta }{=}(1+(P_s{\left|{\mathbf{h}}^H\mathbf{Wf}\right|}^2/{\sigma }_I^2+{\mathbf{h}}^H\mathbf{W}{\mathbf{W}}^H\mathbf{h}+{\mathbf{h}}^H\mathbf{Qh}))2^{-2R_s}$ , then $\{\mathbf{W},\mathbf{Q},\gamma \}$ is a feasible solution to equation (11). However, let $\{\mathbf{W},\mathbf{Q},\gamma \}$ denote a feasible solution to equation (11); it is easy to see that $\{\mathbf{W},\mathbf{Q}\}$ is also a feasible for equation (8). Therefore, equations (8) and (11) are equivalent.

Notably, the objective and constraints in equation (11) are highly non-convex; thus, it is hard to obtain the solution directly. In the following, we first solve it for fixed $\gamma$ and then optimize $\gamma$ by studying its hidden convexity. Then, we propose an efficient iterative algorithm combining gradient-based method and semi-definite relaxation (SDR).

Optimization of the SDR of equation (11) for fixed $\gamma$

Let us denote $\mathbf{w}\stackrel{\Delta }{=}\mathrm{vec}\left(\mathbf{W}\right)$ and $\mathcal{W}\stackrel{\Delta }{=}\mathbf{w}{\mathbf{w}}^H$ , which means that $\mathcal{W}⪰0$ and $\mathrm{rank}\left(\mathcal{W}\right)=1$ . As a standard routine of SDR, we drop the non-convex constraint $\mathrm{rank}\left(\mathcal{W}\right)=1$ ; then, we have the following proposition.

Proposition 1

The SDR reformation of equation (11) can be defined as

$$\underset{\genfrac{}{}{0ex}{}{\mathcal{W}⪰0,\mathbf{Q}⪰0,}{\begin{array}{c}\mu ,{\theta }_m\hfill \end{array}}}{min}\mathrm{Tr}\left(\overline{C}\mathcal{W}\right)+\mathrm{Tr}\left(\mathbf{Q}\right)$$

(12a)

$$\begin{array}{c}\mathrm{s}.\mathrm{t}.\mathrm{Tr}\left(\mathbf{A}\mathcal{W}\right)\ge \hfill \\ \multicolumn{1}{c}{(\gamma 2^{2R_s}-1)({\sigma }_I^2+\mathrm{Tr}\left(\mathbf{B}\mathcal{W}\right)+{\mathbf{h}}^H\mathbf{Qh})}\end{array}$$

(12b)

$$\mathbf{\Omega }(\mathcal{W},\mathbf{Q},\gamma ,\mu )⪰0,\mu \ge 0$$

(12c)

$${\mathbf{\Phi }}_m(\mathcal{W},\mathbf{Q},{\theta }_m)⪰0,{\theta }_m\ge 0,\forall m$$

(12d)

where

$$\begin{array}{c}\mathbf{\Omega }(\mathcal{W},\mathbf{Q},\gamma ,\mu )\stackrel{\Delta }{=}\hfill \\ \multicolumn{1}{c}{\left[\begin{array}{c}{\overline{\mathbf{G}}}_E^H\hfill \\ \multicolumn{1}{c}{\mathbf{I}}\end{array}\right]\mathbf{U}(\mathcal{W},\mathbf{Q},\gamma )\left[\begin{array}{cc}{\overline{\mathbf{G}}}_E\hfill & \mathbf{I}\hfill \end{array}\right]}\\ \multicolumn{1}{c}{+\mathbf{D}(({\sigma }_E^2-\mu )\mathbf{I},\mu \mathbf{\Lambda })⪰0}\end{array}$$

(13a)

$$\begin{array}{c}{\mathbf{\Phi }}_m(\mathcal{W},\mathbf{Q},{\theta }_m)\stackrel{\Delta }{=}\hfill \\ \multicolumn{1}{c}{\left[\begin{array}{c}\mathbf{I}\hfill \\ \multicolumn{1}{c}{{\overline{\mathbf{g}}}_E^H}\end{array}\right](\mathbf{D}\left({\mathbf{e}}_m\right)\otimes \mathbf{V}(\mathcal{W},\mathbf{Q}))\left[\begin{array}{cc}\mathbf{I}\hfill & {\overline{\mathbf{g}}}_E\hfill \end{array}\right]}\\ \multicolumn{1}{c}{+\mathbf{D}({\theta }_m\mathbf{D}\left({\mathbf{e}}_m\right)\otimes \mathbf{\Lambda },-{\theta }_m-{\eta }_m/{\rho }_m)⪰0}\end{array}$$

(13b)

and $\mu \ge 0$ , ${\theta }_m\ge 0$ denote the related auxiliary variables.

In addition, $\mathbf{A}$ , $\mathbf{B}$ , $\overline{\mathbf{C}}$ , $\mathbf{U}(\mathcal{W},\mathbf{Q},\gamma )$ , and $\mathbf{V}(\mathcal{W},\mathbf{Q})$ are defined as follows

$$\mathbf{A}=P_s({\mathbf{f}}^{†}\otimes \mathbf{h}){({\mathbf{f}}^{†}\otimes \mathbf{h})}^H$$

(14a)

$$\mathbf{B}={\sigma }_r^2\mathbf{I}\otimes \mathbf{h}{\mathbf{h}}^H$$

(14b)

$$\overline{\mathbf{C}}={\mathbf{C}}^T\otimes \mathbf{I}$$

(14c)

$$\mathbf{U}(\mathcal{W},\mathbf{Q},\gamma )\stackrel{\Delta }{=}{\sigma }_r^2\sum _{\ell =1}^{N_r}{\mathbf{E}}_{\ell }\mathcal{W}{\mathbf{E}}_{\ell }^H-\frac{P_s}{\gamma -1}\mathbf{F}\mathcal{W}{\mathbf{F}}^H+\mathbf{Q}$$

(14d)

$$\mathbf{V}(\mathcal{W},\mathbf{Q})\stackrel{\Delta }{=}{P}_s\mathbf{F}\mathcal{W}{\mathbf{F}}^H+{\sigma }_r^2\sum _{\ell =1}^{N_r}{\mathbf{E}}_{\ell }\mathcal{W}{\mathbf{E}}_{\ell }^H+\mathbf{Q}$$

(14e)

with $\mathbf{F}\stackrel{\Delta }{=}{\mathbf{f}}^T\otimes \mathbf{I}$ and ${\mathbf{E}}_{\ell }\stackrel{\Delta }{=}[0,\dots ,\mathbf{I},\dots ,0]\in {\mathbb{R}}^{N_r\times N_r^2}$ is a zero matrix, except for the columns ranging from $N_r^2+(\ell -1)N_r+1$ to $N_r^2+\ell N_r$ which form an identity matrix $\mathbf{I}$ . The proof of proposition 1 is detailed in Appendix 1.

The main advantage of equation (12) is that for fixed $\gamma$ , it is an SDR problem and can be efficiently solved by available convex programming toolbox CVX. ³⁶ It should be noted that the obtained matrix $\mathcal{W}$ might not be rank one; therefore, equation (12) is a relaxation of equation (11) for fixed $\gamma$ . For the rank of $\mathcal{W}$ , we have the following proposition.

Proposition 2

The optimal solution ${\mathcal{W}}^{*}$ to equation (12) can be expressed as

$${\mathcal{W}}^{*}=\sum _{n=1}^{N_r^2-r_c}{a}_n{\mathbf{\pi }}_n{\mathbf{\pi }}_n^H+b\mathbf{\varsigma }{\mathbf{\varsigma }}^H$$

(15)

where $a_n\ge 0$ , $\forall n$ , $b>0$ , $r_c=\mathrm{rank}({\mathbf{M}}^{*})$ (c.f. equation (37)), and $\mathbf{\varsigma }\in {\mathbb{C}}^{N_r^2\times 1}$ is a unit-norm vector orthogonal to $\mathbf{N}=\{{\mathbf{\pi }}_n{\}}_{n=1}^{N_r^2-r_c}$ , which consists of the orthonormal basis for $\mathrm{null}({\mathbf{M}}^{*})$ . The proof of proposition 2 is detailed in Appendix 1.

According to equation (15), if $\mathrm{rank}({\mathcal{W}}^{*})>1$ , that is, there exists at least one $a_n>0$ , we reconstruct a solution to equation (12) using

$${\hat{\mathcal{W}}}^{*}=b\mathbf{\varsigma }{\mathbf{\varsigma }}^H$$

(16)

while ${\hat{\mathbf{Q}}}^{*}$ are obtained by solving the following feasible problem provided that ${\hat{\mathcal{W}}}^{*}$ is given by equation (16)

$$\underset{\hat{\mathbf{Q}}⪰0}{min}0$$

(17a)

$$\mathrm{s}.\mathrm{t}.\mathrm{given}{\hat{\mathcal{W}}}^{*}$$

(17b)

$$(12\mathrm{b})-(12\mathrm{d})$$

(17c)

Optimization over $\gamma$

In the above subsection, we have converted equation (11) to the SDR (equation (12)) for fixed $\gamma$ . Denote the objective value of equation (12) as a function of $\gamma$ by $\mathbb{P}\left(\gamma \right)$ . It can be seen that solving the original problem (equation (11)) is equivalent to minimizing $\mathbb{P}\left(\gamma \right)$ with respect to $\gamma$ where the optimization of $\{\mathcal{W},\mathbf{Q}\}$ has been embedded in evaluating $\mathbb{P}\left(\gamma \right)$ .

Next, we focus on finding the optimal $\gamma$ to minimize $\mathbb{P}\left(\gamma \right)$ . We have the following proposition.

Proposition 3

$\mathbb{P}\left(\gamma \right)$ is a convex function of $\gamma$ .

Proof

The Lagrangian of equation (12) is

$$\begin{array}{cc}\multicolumn{1}{c}{\mathcal{L}\left(\mathcal{X}\right)=}& \mathrm{Tr}\left(\overline{\mathbf{C}}\mathcal{W}\right)+\mathrm{Tr}\left(\mathbf{Q}\right)-\alpha \mathrm{Tr}\left(\mathbf{A}\mathcal{W}\right)\hfill \\ \multicolumn{1}{c}{}& +\alpha (\gamma 2^{2R_s}-1)({\sigma }_I^2+\mathrm{Tr}\left(\mathbf{B}\mathcal{W}\right)+\mathrm{Tr}\left(\mathbf{h}{\mathbf{h}}^H\mathbf{Q}\right))\hfill \\ \multicolumn{1}{c}{}& -\mathrm{Tr}\left(\mathbf{\Xi }\mathbf{\Omega }\right)-\sum _{m=1}^M\mathrm{Tr}\left({\mathbf{T}}_m{\mathbf{\Phi }}_m\right)\hfill \\ \multicolumn{1}{c}{}& -\mathrm{Tr}\left(\Upsilon \mathcal{W}\right)-\mathrm{Tr}\left(\mathbf{\Theta }\mathbf{Q}\right)\hfill \end{array}$$

(18)

where $\mathcal{X}$ denotes the collection of the primal and dual variables. Specifically, $\alpha \ge 0$ , $\mathbf{\Xi }\in {\mathbb{H}}_{+}^{N_e+N_r}$ , ${\mathbf{T}}_m\in {\mathbb{H}}_{+}^{N_e\times N_r+1}$ , $\Upsilon \in {\mathbb{H}}_{+}^{N_r\times N_r}$ , and $\mathbf{\Theta }\in {\mathbb{H}}_{+}^{N_r}$ are the Lagrangian multipliers associated with equation (12b)–(12d) and the primary variables $\{\mathcal{W},\mathbf{Q}\}$ , respectively.

The dual problem can be written as

$$\begin{array}{c}\underset{\genfrac{}{}{0ex}{}{\mathbf{\Xi }⪰0,{\mathbf{T}}_m⪰0,}{\begin{array}{c}\alpha \ge 0\hfill \end{array}}}{max}\alpha (\gamma 2^{2R_s}-1){\sigma }_I^2-\mathrm{Tr}\left(\mathbf{D}(({\sigma }_E^2-\mu )\mathbf{I},\mu \mathbf{\Lambda })\mathbf{\Xi }\right)\hfill \\ -\sum _{m=1}^M\mathrm{Tr}\left(\mathbf{D}({\theta }_m\mathbf{D}\left({\mathbf{e}}_m\right)\otimes \mathbf{\Lambda },-{\theta }_m-{\eta }_m/{\rho }_m){\mathbf{T}}_m\right)\hfill \end{array}$$

(19a)

$$\begin{array}{c}\mathrm{s}.\mathrm{t}.\overline{\mathbf{C}}+\alpha (\gamma 2^{2R_s}-1)\mathbf{B}-\alpha \mathbf{A}-{\nabla }_{\mathcal{W}}\mathrm{Tr}\left(\mathbf{\Xi }\mathbf{\Omega }\right)\hfill \\ \multicolumn{1}{c}{-\sum _{m=1}^M{\nabla }_{\mathcal{W}}\mathrm{Tr}\left({\mathbf{T}}_m{\mathbf{\Phi }}_m\right)⪰0}\end{array}$$

(19b)

$$\begin{array}{c}\mathbf{I}+\alpha (\gamma 2^{2R_s}-1)\mathbf{h}{\mathbf{h}}^H-{\nabla }_{\mathbf{Q}}\mathrm{Tr}\left(\mathbf{\Xi }\mathbf{\Omega }\right)\hfill \\ \multicolumn{1}{c}{-\sum _{m=1}^M{\nabla }_{\mathbf{Q}}\mathrm{Tr}\left({\mathbf{T}}_m{\mathbf{\Phi }}_m\right)⪰0}\end{array}$$

(19c)

Based on the operation of matrix differential ³⁷ and using the fact $\sum _{\ell =1}^{N_r}{\mathbf{E}}_{\ell }^H\mathbf{\Pi }{\mathbf{E}}_{\ell }=\mathbf{I}\otimes \mathbf{\Pi }$ , we obtain

$${\nabla }_{\mathbf{Q}}\mathrm{Tr}\left(\mathbf{\Xi }\mathbf{\Omega }\right)=\mathbf{\Pi }$$

(20a)

$${\nabla }_{\mathcal{W}}\mathrm{Tr}\left(\mathbf{\Xi }\mathbf{\Omega }\right)=-\frac{P_s}{\gamma -1}{\mathbf{F}}^H\mathbf{\Pi }\mathbf{F}+{\sigma }_r^2\mathbf{I}\otimes \mathbf{\Pi }$$

(20b)

where $\mathbf{\Pi }=\left[{\overline{\mathbf{G}}}_E\mathbf{I}\right]\Omega {\left[{\overline{\mathbf{G}}}_E\mathbf{I}\right]}^H$ .

Similarly, we attain the following relationships

$${\nabla }_{\mathbf{Q}}\mathrm{Tr}\left({\mathbf{T}}_m{\mathbf{\Phi }}_m\right)={\mathbf{Z}}_m^{(m,m)}$$

(21a)

$${\nabla }_{\mathcal{W}}\mathrm{Tr}\left({\mathbf{T}}_m{\mathbf{\Phi }}_m\right)={\sigma }_r^2\mathbf{I}\otimes {\mathbf{Z}}_m^{(m,m)}+P_s{\mathbf{F}}^H{\mathbf{Z}}_m^{(m,m)}\mathbf{F}$$

(21b)

where ${\mathbf{Z}}_m^{(m,m)}$ are the block submatrices of matrix $[\begin{array}{cc}\mathbf{I}\hfill & {\overline{\mathbf{g}}}_E\hfill \end{array}]{\mathbf{T}}_m[\begin{array}{cc}\mathbf{I}\hfill & {\overline{\mathbf{g}}}_E\hfill \end{array}{]}^H$ , which can be expressed specifically as

$$[\begin{array}{cc}\mathbf{I}\hfill & {\overline{\mathbf{g}}}_E\hfill \end{array}]{\mathbf{T}}_m[\begin{array}{cc}\mathbf{I}\hfill & {\overline{\mathbf{g}}}_E\hfill \end{array}{]}^H=\left[\begin{array}{ccc}{\mathbf{Z}}_m^{(1,1)}\hfill & \cdots \hfill & {\mathbf{Z}}_m^{(1,M)}\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\mathbf{Z}}_m^{(M,1)}}& \cdots \hfill & {\mathbf{Z}}_m^{(M,M)}\hfill \end{array}\right]$$

(22)

Due to the strong duality between the primary problem (12) and the dual problem (19), the objective value of equation (19) is also $\mathbb{P}\left(\gamma \right)$ . Observe that for given $\{\alpha ,\mathcal{W},\mathbf{Q}\}$ , the objective function of equation (19) is linear and convex about $\gamma$ ; therefore, $\mathbb{P}\left(\gamma \right)$ is convex about $\gamma$ . This completes the proof. The proof is necessary since only when $\mathbb{P}\left(\gamma \right)$ is convex about $\gamma$ , the gradient-based algorithm can obtain the minimum value for the power minimization problem. Since it is hard to proove $\mathbb{P}\left(\gamma \right)$ is convex about $\gamma$ directly from equation (12), we resort to the dual problem (19) and from the convexity of equation (19), we obtain Proposition 3.

Proposition 3 shows that any local search algorithm can find the global optimal value of $\mathbb{P}\left(\gamma \right)$ , therefore greatly reduces the complexity. We have the following derivative information to design a gradient-based algorithm. (Since $\gamma$ are nonnegative variables with lower bound $1$ (due to (11c)) and $\mathbb{P}\left(\gamma \right)$ is a convex function about $\gamma$ , the optimal solution can be found by one-dimensional search; however, in practice, the one-dimensional search usually executes unidirectional search within the entire feasible domain, while the search direction of the gradient-based algorithm always points to the optimal point and the search area is usually near the optimal point. This makes the gradient-based algorithm more efficient, which will be further confirmed by the simulation result.)

The derivative of $\mathbb{P}\left(\gamma \right)$ is given by

$$\begin{array}{c}\frac{\partial \mathcal{L}}{\partial \gamma }=\alpha 2^{2R_s}({\sigma }_I^2+\mathrm{Tr}\left(\mathbf{B}\mathcal{W}\right)+\mathrm{Tr}\left(\mathbf{h}{\mathbf{h}}^H\mathbf{Q}\right))\hfill \\ \multicolumn{1}{c}{-\frac{P_s}{{(\gamma -1)}^2}\mathrm{Tr}\left(\mathbf{F}\mathbf{W}{\mathbf{F}}^H\mathbf{\Pi }\right)}\end{array}$$

(23)

The overall iterative algorithm to solve equation (11) is given in Algorithm 1. (The initial value of $\gamma$ can be set slightly larger than 1 to get the simulation result in relative less time since the optimal value is often obtained for relative small $\gamma$ . In fact, we set $\gamma =1.01$ in our simulations; for the convergence behavior of the iterative procedure, the convergence behavior can be confirmed since the proposed algorithm generated a decreasing objective at each iteration and it is lower bounded, for example, the optimal value is the lower bound, thus convergence is guaranteed; for the terminate criterion of Algorithm 1, we set a parameter k that denotes the gap of the objective (19) in successive iteration. When the gap is smaller than a given threshold $k\mathrm{th}$ , the algorithm will be terminated.)

Algorithm 1 Proposed algorithm to solve equation (11).

1: Initialization: set $\mathbf{f},\mathbf{h},{\overline{\mathbf{G}}}_E,{\sigma }_r^2,{\sigma }_I^2,{\sigma }_E^2,R_s,P_s,{\rho }_m,{\eta }_m,\mathbf{\Lambda }$ .2: begin  1. Solve the primary problem (12) and the dual problem (19) using the CVX solver to get $\mathcal{W}$ , $\mathbf{Q}$ , $\alpha$ , $\mathbf{X}$ and ${\mathbf{T}}_m$ .  2. Using the derivative information $\frac{\partial \mathcal{L}}{\partial \gamma }$ given in equation (23) to update $\gamma$ .  3. Repeat Step 2.1 and Step 2.2 until the terminate criterion has been arrived and obtain ${\mathcal{W}}^{★},{\mathbf{Q}}^{★}$ .  4. If $\mathrm{rank}\left({\mathcal{W}}^{★}\right)=1$ , the optimal ${\mathbf{W}}^{★}$ can be obtained by ${\mathbf{W}}^{★}=\mathrm{ve}{\mathrm{c}}^{-1}\left({\mathbf{w}}^{★}\right)$ , where ${\mathbf{w}}^{★}$ is a rank-one decomposition of ${\mathcal{W}}^{★}$ , otherwise goto Step 2.5.  5. Using equation (16), solve equation (17) to attain the suboptimal ${\mathcal{W}}^{★},{\mathbf{Q}}^{★}$ .3: End4: Output $({\mathbf{W}}^{★},{\mathbf{Q}}^{★})$ .

Simulation result

In this section, we evaluate the performance of our scheme through Monte Carlo simulations. The following parameters $N_r=4$ , $M=2$ , $R_s=2\mathrm{bps}/\mathrm{Hz}$ , $P_s=10\mathrm{dBW}$ , ${\sigma }_r^2={\sigma }_I^2={\sigma }_E^2=-50\mathrm{dBm}$ and ${\rho }_m=0.5,{\eta }_m=-20\mathrm{dBm},\forall m$ are set for the simulation examples unless specified. In addition, we assume that all the entries of channel responses $\mathbf{f}$ , $\mathbf{h}$ , and ${\overline{\mathbf{G}}}_E$ are independent and identically distributed (i.i.d.) complex Gaussian random variables generated by $\mathcal{C}N(0,{10}^{-2})$ . Regarding the CSI uncertainty model in equation (4), we use a similar way as Xing et al., ²⁶ for example, we introduce the uncertainty ratio $\upsilon$ associated with $\epsilon$ which is defined as ${\upsilon }^2={\epsilon }^2/\mathbb{E}[\parallel {\overline{\mathbf{G}}}_E{\parallel }_F^2]$ . In addition, the CSI uncertainty regions are assumed to be norm-bound, that is, $\mathbf{\Lambda }\mathbf{=}\mathbf{I}\mathbf{/}{\mathbf{\epsilon }}^2$ , and we set ${\upsilon }^2=0.1$ .

To highlight the superiority of our proposed scheme, we compare our algorithm with several other methods: (1) robust BF and AN with one-dimensional search; (2) BF with null space AN, which is a modified version of the null space AN method in the work by Li et al., ²⁴ for example, the AN is transmitted in the null space of the $\mathrm{h}$ ; (3) null space BF with AN, which is a modified version of the null space BF method in the work by Liu et al. ¹¹ Specifically, we project the BF matrix into the null space of the Eve’s channel, for example, the BF matrix can be written as $\mathcal{W}=\mathbf{UX}{\mathbf{U}}^H$ , where $\mathbf{U}$ is the null space of the transmitter to the colluding Eve’s channel ${\overline{\mathbf{G}}}_E^H\mathbf{F}$ , for example, $\mathbf{U}=\mathrm{null}\left({\overline{\mathbf{G}}}_E^H\mathbf{F}\right)$ . Then, we substitute this $\mathcal{W}$ into equations (12) and (19) to obtain the optimal $\mathbf{X}$ and $\mathbf{Q}$ using the gradient-based algorithm. It should be noted that the null space BF scheme can only be valid when the number of relay antennas is not smaller than the number of ERs. Our design and the other designs are labeled as “gradient-based BFAN,”“null space AN,” ²⁴ “null space BF,” ¹¹ and “one-dimensional search BFAN,” respectively.

The relay power versus the secrecy rate R_S

Figure 2 illustrates the relay transmit power of these schemes versus the secrecy rate $R_s$ . It is observed from Figure 2 that the relay power increases with the growth of $R_s$ for all the schemes, while our gradient-based BFAN design outperforms the other schemes and the null space AN is the worst design. In addition, we can see from the figure that when $R_s$ is small, the relay power increases slowly, whereas when $R_s$ is large, the power relay increases faster.

Figure 2.

The relay power versus the worst-case secrecy rate.

The relay power versus the ERs’ CSI uncertainty level ${\upsilon }^2$

Figure 3 shows the relay transmit power of these schemes versus the ERs’ CSI uncertainty level ${\upsilon }^2$ . It is observed from the figure that the relay power increases with an increase in ${\upsilon }^2$ for all the schemes, while our gradient-based BFAN scheme achieves the best performance. In addition, the relay transmit power shows an approximately linear relationship with ${\upsilon }^2$ , which is different with the phenomenon in Figure 2.

Figure 3.

The relay power versus the ERs’ CSI uncertainty level.

The relay power versus the EH threshold ${\eta }_m$

Figure 4 investigates the relay transmit power of these schemes versus the EH threshold ${\eta }_m$ . It is seen that the relay power budget increases with an increase in ${\eta }_m$ for all the schemes, while our gradient-based BFAN scheme achieves better performance than the other schemes. In addition, the relay transmit power shows an approximately linear relationship with ${\eta }_m$ , which is different with the phenomenon in Figure 2 but similar with the result in Figure 3. The result suggests that when the secrecy rate is small, the ERs’ CSI uncertainty level and the EH threshold are more prominent to the relay power consumption. However, when the secrecy rate is high, the secrecy rate becomes the major factor for the relay power consumption.

Figure 4.

The relay power versus the EH threshold.

The relay power versus the number of relay antennas N_r

Figure 5 plots the relay power consumption of these schemes versus the number of relay antennas $N_r$ . It is seen that the power budget decreases with an increase in $N_r$ for all the methods due to the increased spatial DoF, while our scheme outperforms other schemes. In addition, we can see that the relay power decreases very quickly when $N_r$ is small, whereas when $N_r$ is large, the relay power decreases slowly.

Figure 5.

The relay power versus the number of relay antennas.

The relay power versus the number of ERs M

Figure 6 depicts the relay transmit power of these schemes versus the number of ERs $M$ . As expected, the relay power consumption increases with an increase in $M$ for all the methods due to the decreased spatial DoF and more EH constraints to be satisfied, while our scheme has a significant performance gain over the other schemes. It should be noted that the null space BF scheme is invalid when $M\ge 4$ in this example since there is no enough DoF left for the IR after nulling the ER’s channel.

Figure 6.

The relay power versus the number of ERs.

The secrecy rate and EH performance versus the relay power P_r

Finally, Figures 7 and 8 depict the secrecy rate performance and the EH performance with relay transmit power, respectively. In both of the examples, we use the one-dimensional search method in the works by Li et al. ²¹ and Feng et al. ²⁷ with several necessary modifications, and the joint design is labeled as “robust BFAN.” From Figures 7 and 8, we can find the joint design achieves higher secrecy rate and EH performance than the null space BF and null space AN design, which further demonstrates the superiority of the joint robust BF and AN design.

Figure 7.

The secrecy rate versus the relay power.

Figure 8.

The total harvested energy versus the relay power.

Conclusion

In this article, we investigated a joint robust BF and AN scheme for secure SWIPT in AF relay WSN. Specifically, we formulated the robust relay power minimization problem subject to both secrecy rate and EH constraints. To obtain the optimal BF matrix and AN covariance, an iterative algorithm combining gradient-based method and SDR has been proposed to find the optimal solution by analyzing its convexity. In addition, a suboptimal rank-one reconstruction algorithm was presented. Simulation results demonstrated the effectiveness of the proposed design.