Joint Source-Relay Optimization for MIMO Full-Duplex Bidirectional Wireless Sensor Networks with SWIPT

The simultaneous wireless information and power transfer (SWIPT) technique has been considered as a promising approach to prolong the lifetime of energy-constraint wireless sensor networks (WSNs). In this paper, a multiple-input multiple-output (MIMO) full-duplex (FD) bidirectional wireless sensor network (BWSN) with SWIPT is investigated. Based on minimum total mean-square-error (total-MSE) criterion, a joint optimization problem for source and relay beamforming and source receiving subject to transmitting power and harvesting energy constraints is established. Since this problem is non-convex, an iterative algorithm based on feasible point pursuit-successive convex approximation (FPP-SCA) is derived to obtain a local optimum. Moreover, considering the scenarios in which source and relay nodes equipped with the same and different numbers of antennas, a low-complexity diagonalizing design-based scheme is employed to simplify each non-convex subproblem into convex problems and to reduce the computational complexity. Numerical results of the total-MSE and bit error rate (BER) are implemented to demonstrate the performance of the two different schemes.


Introduction
Wireless sensor networks (WSNs) have attracted a significant amount of attention from researchers and have been widely employed in vast and varied areas, e.g., object tracking, habitat monitoring, military systems, and industrial areas [1][2][3]. However, in WSNs, the relay or sensor nodes are typically powered by batteries with finite capacities [4], which are difficult or impossible to replace or recharge in most cases. Thus, the energy supplies will limit the lifetime of WSNs. Saving on energy or prolonging the operation time of energy-constrained nodes has become an important research issue in WSNs. Traditionally, multi-input multi-output (MIMO) can provide an effective way for energy saving [5,6].
Recently, simultaneous wireless information and power transfer (SWIPT) is considered a promising energy-harvesting (EH) technique to solve the energy scarcity problem and to achieve perpetual communications in energy-constrained WSNs [7][8][9], which is extensively applied in the area. To date, two receiver architectures proposed in Reference [10], namely time switching (TS) and power splitting (PS), have been widely used for a colocated energy harvester and information decoder in SWIPT systems [11,12]. Compared with the TS structure periodically switching between the EH module and information decoding (ID) module, the PS design allows the receiver to complete EH and iterative algorithm is exploited. Finally, to reduce the computational complexity, a low-complexity diagonalizing method-based algorithm is introduced to simplify each non-convex subproblem into convex problems directly. In terms of the existing approach [25], the generalized singular value decomposition (GSVD) is discussed, and the scenarios in which source and relay nodes equipped with the same and different number of antennas are both discussed. The numerical results show a good performance and validate our analysis.
The remainder of the paper is organized as follows. Section 2 proposes the system model, including the sensor nodes deployment and optimization model. Section 3 focuses on the scheme design. The numerical results are presented and discussed in Section 4. Finally, the conclusions are presented in Section 5.
Notation: Throughout this paper, scalar variables are expressed by lowercase italic letters, vectors are represented by boldface lowercase letters, and matrices are denoted by boldface uppercase letters.
C M×N denotes an M × N matrix with complex entries. Tr(·), (·) T , (·) H , (·) −1 , (·) * , and · are the trace, transpose operation, conjugate transpose operation, inverse operation, conjugate transpose operation, and the bound norm of a vector. ∑ M i=1 (·) stands for the sum from 1 to M. ∼ CN x, σ 2 represents a complex Gaussian distributed variable with a mean x and covariance σ. vec(·) and mat(·) are the matrix vectorization operator and the corresponding inverse operation, respectively. E[·] signifies the expectation of the random variables in the bracket.

System Model
This paper aims to jointly design the transmitters of the source and relay and the receivers of the source in the FD BWSN with SWIPT. We adopt a three-node sensor system consisting of two sources and a relay and assume that the sources are equipped with the PS receiver and that the relay applies an AF scheme. Without a loss of generality, we suppose that the energy conversion efficiency at the PS receiver is 100 percent and that the PS ratio is fixed.
The considered three-node MIMO BWSN with SWIPT consists of two sources: S 1 and S 2 both equipped with M T > 1 transmit antennas and M R > 1 receive antennas. S 1 and S 2 decode the information, harvest the energy by PS, and exchange information with the help of the single AF relay node R with N T > 1 transmit antennas and N R > 1 receive antennas, as shown in Figure 1. All nodes are assumed to operate in FD mode, which means they transmit and receive data at the same time and frequency. Let H S i R ∈ C N R ×M T and H RS i ∈ C M R ×N T denote channel matrices from S i 's transmit antennas to R's receive antennas and that from R's transmit antennas to S i 's receive antennas, respectively. We assume that all channels are statistically independent, reciprocal in the incoming and outgoing directions, and slowly time-varying quasi-static flat Rayleigh fading. Moreover, the self-interference channels at the corresponding nodes are represented as Meanwhile, in our system, the two source nodes S 1 and S 2 are set far apart so that the direct link between them is assumed to be ignored. Moreover, we suppose that the perfect channel state information (CSI) is available at each node [38][39][40] and that the transmit power of the two sources are equal.
For a further analysis, the node deployment and optimization models are presented as follows.

Node Deployment
At time instant n(n > 1), M T data streams s i [n] ∈ C M T ×1 with a normalized power are transmitted through the beamformer F i ∈ C M T ×M T (i ∈ {1, 2}) from S i simultaneously and the relay R forwards its received signal y R [n] ∈ C N R ×1 after multiplying it by a beamforming matrix F r ∈ C N T ×N R . In practice, a τ(τ 1)-symbol processing delay is unavoidable at R when it processes the received signals.
Accordingly, in time slot n, the received signal y R [n] at R can be expressed as where n r [n] ∼ CN 0, σ 2 r I N r represents the additive white Gaussian noise (AWGN) at R. Assuming that τ is kept small and that the SI can be cancelled perfectly or almost perfectly with the knowledge of the signal transmitted by the relay itself [37], the transmitted signal x R [n] ∈ C N T ×1 at R can be written as Substituting Equation (1) into Equation (2), the overall relay output can be given by The signal received by S i , i ∈ {1, 2} can be written as where, the same as below, i = 2 if i = 1 and vice versa. n i [n] is the equivalent noise vector For simplicity, we assume that the full channel state information (CSI) is known and that both S 1 and S 2 know their own transmitted signals; thus, the SI at S i herein can be cancelled. After subtracting the back-propagated self-interference term To implement SWIPT, a portion β i ∈ (0, 1) of the signal power is applied to Equation (5), which splits y S i [n] into two parts, β i portion for ID and the remaining (1 − β i ) portion for EH. Then, the signals for ID at each source node can be represented as where denotes the equivalent noise vector and n p [n] ∼ CN 0, σ 2 p I M R is the AWGN caused by the power splitter. At the EH side, we have where ζ i ∈ (0, 1] is the energy conversion efficiency at the energy harvester. It is assumed that for ζ i = 1, here, e i represents the minimum power that should be harvested at S i and Moreover, F i and F r should satisfy the transmitting power constraints, that is, where p i and p r are the maximum transmit power supplied by S i and R, respectively. Since channels in our system are memoryless, we can define that y S i , and n r n r [n − τ], and we assume that σ 2 p σ 2 S i . Thus, (6) can be reformulated as

Optimization Model
Considering the harvested energy and transmit power constraints, i.e., Equations (7) and (8), the optimization model to minimize the total-MSE of the whole system and to find the optimal source and relay beamformer and the source receiver is described in this section. The objective function and the problem are separately discussed below.
Using Equation (9), the MSE of S i can be given by where Given the MSE of S i , fixing β i , a joint source and relay beamforming and source-receiving optimization problem based on the total-MSE with transmit power constraints and an energy-harvesting constraint can be formulated as

Scheme Design
Considering the problem Equation (11) is non-convex and multivariate, the iterative algorithms based on FPP-SCA and a low-complexity diagonalizing are employed in this section.

Iterative Algorithm Based on FPP-SCA
Since the problem in Equation (11) is non-convex and basically intractable, in this section, an iterative algorithm based on FPP-SCA [41] is proposed to decouple the primal problem into four subproblems corresponding to four variables: W i , F r , F 1 , and F 2 , and to solve them alternately. At each iteration, one variable is optimized while keeping the other fixed. Starting from Equation (12), the W i is optimized, and then, the F r is optimized by using Equation (14), following this, Equation (16) (actually two subproblems) is formulated to optimize F 1 and F 2 separately. Finally, the four subproblems are solved, and the four variables are optimized. Details are given below.
First, with F i and F r fixed, the receiver W i is first optimized. As W i is only involved in J i , the optimal W opt i can be derived using ∂J i /∂W * i = 0, which yields where

Optimization of Relay Beamformer F r
Then, the optimization of F r with a fixed F i and W i is discussed. According to [42] (p. 77), where A, B, C, and D are arbitrary matrices with compatible dimensions, ⊗ is the Kronecker product, and vec(·) represents the matrix vectorization operator.
To guarantee the feasibility of Equation (11), the feasible region is relaxed and approximated by adding slacks s ∈ R 2 to the non-convex constraint of Equation (11d), and the positive slack variables and slack penalty are used in Equation (11a). Then, the original problem can be recast as where · can be any vector norm, s denotes the slack penalty term, and λ 1 is the trade-off between the original objective function and s . Besides, (14) is non-convex. To tackle this subproblem, we define g(f r ) = f H r Q ir f r and assume that a center point z r ∈ C N×1 , N = N T × N R is given. Introducing g(f r ) 2Re z H r Q ir f r − z H r Q ir z r , Theorem 1 can be established and proved.
Theorem 1.g satisfies the following properties: (i)g(z r ) = g(z r ); (ii)g(f r ) ≥ g(f r ); and (iii) Proof. Substituting z r intog(f r ) and g(f r ), (i) can be easily certified. For (ii), , the derivatives can be computed as Replacing g(f r ) withg(f r ), Equation (14) can be rewritten as Equation (15) can be efficiently solved using the modeling language YALMIP [43] and the generic conic programming solver SeDuMi [44]. A new approximated problem can be built and solved when the optimal solution of Equation (15) becomes the new center point, that is, z r = f * r . Based on Theorem 1, Equation (14) can be solved.

Optimization of Source Beamformer F i
Similarly, F i can be optimized given F r and W i . According to Equation (13), the original problem in Equation (11) can be transformed into Equation (16) min . The optimized f 1 and f 2 can be separately obtained from Equation (16) using a similar FPP-SCA method foresaid.

Summarization of the Proposed Algorithm
Based on the FPP-SCA algorithm presented above, the iterative algorithm is summarized as Algorithm 1.
Algorithm 1 An alternating optimization algorithm based on a feasible point pursuit-successive convex approximation (FPP-SCA)

Initialize
Define β i = 0.5, F r = pr Tr(R) I NR , and F i = p i MT I MT .

Iterative updating
(1) Update W i using Equation (12) with a fixed F i and F r .
(2) Update F r by solving Equation (14) a. Set k = 0 and vec(F r ) as the initial point z 0 r . b. Solve Equation (15) at the kth iteration for k 0 to yield the optimal solution f k r . with a fixed F i and W i . c. Let z k+1 r = f k r and k = k + 1. d. Until convergence, let F r = mat(z k+1 r ).
(3) Update F i by solving Equation (16) with a fixed F r and W i , following similar steps in (2).

Until convergence
The Algorithm 1 is convergent based on the following Property 1.

Property 1. The iterative algorithm based on FPP-SCA is convergent.
Proof. In the kth iteration of the proposed algorithm, we first compute F [k] r with the given F i is optimally solved and the objective value is descendent. Consequently, the objective value of the original problem monotonically decreases and is lower-bounded by zero, which verifies the convergence of Algorithm 1.

Low-Complexity Diagonalizing Design
However, the main drawback of the proposed FPP-SCA algorithm is the high computational complexity. In order to overcome this shortcoming, a low-complexity algorithm using the channel parallelization (CP) technique [22], namely the generalized singular value decomposition (GSVD) and SVD, is applied.
In this section, we assume that M T = M R = M and N T = N R = N for simplicity and focus on the scenarios where N M.

Channel Parallelization
Substituting Equation (12) into Equation (10) and employing where E and F are arbitrary matrices and I is the identity matrix, the function J i can be simplified as Applying GSVD on the uplink channel matrix pair where R h ∈ C N×N , U h i ∈ C M×M , R dl ∈ C 2M×2M , and U dl ∈ C N×N are four unitary matrices; In order to parallelize the channels in Equations (19) and (20), the relay and source beamformers F r and F i can be proposed as where Λ r and Λ i are N × N and M × M nonnegative diagonal matrices, respectively. Substituting Equations (19)-(21) into Equation (18), the resultant objective function J * i becomes where Λ B h and Λ B i are two diagonal matrices containing the (k, k)th entries of B h and B i , (11), the original problem can be expressed as min

Substituting Equations (19)-(21) into each of the constraints in Equation
where To solve the nonconvexity caused by Equation (23c), we propose Theorem 2 as follows.

Theorem 2.
The left side of the energy-harvesting constraint in Equation (23c) can be replaced by its Proof. We take N > M as an example to illustrate the proof procedure and the optimization problem. Expanding the left side of Equation (23c), defining A = R dli Λ dl Λ r , and ignoring the constant matrix Define where

Alternating Optimization of F r and F i
In this section, an iterative approach is utilized to convert the multivariate non-convex problem in Equation (25) into three convex subproblems. We first study how to optimize Λ r with a fixed Λ i , and then, the alternating optimization of Λ 1 and Λ 2 is performed with a given Λ r .

Optimization of Λ r
Using Equation (17), J * i can be simplified and rewritten as where Since Λ r exists in MSE * i only, the problem of minimizing J * 1 + J * 2 is equivalent to that of minimizing MSE * 1 + MSE * 2 . Defining a in , a h i n , a rn , a dn , λ B i n , and λ B h n as the nth diagonal element of Λ i , Λ h i , Λ r , Λ dl , Λ B i , and Λ B h , respectively, MSE * i can be given by where dn . Accordingly, the problem related to Λ r can be described as where

Optimization of Λ i
Similarly, the solution for a in can be described in the following scalar form min By proving and where for i = 1, 1 n M and for i = 2, N − M + 1 n N and b in = σ 2 r β i a 2 dn a 2 rn λ B h n + β i σ 2 S i λ B i n , we can indicate that Equations (29) and (30) are convex for a 2 rn and a 2 in . Then, the optimal solution can be obtained by CVX directly.

Summarization of the Proposed Algorithm
The low-complexity algorithm based on CP method depicted above is summarized as Algorithm 2.

Algorithm 2
The low-complexity algorithm based on the channel parallelization (CP) method

Channel decomposition
Decompose the channel pairs H H S 1 R , H H S2 R and H RS 1 , H RS2 using Equations (19) and (20).

Initialization
Define β i = 0.5, F r = pr Tr(ΨR) I N , and F i =

Iterative updating
(1) Update W i using Equation (12) with a fixed F i and F r .
(2) Update F r with a fixed F i and W i . a. Update Λ r using a rn by solving Equation (29).
(3) Update F i with a fixed F r and W i .

Until convergence
Algorithm 2 is convergent based on the following Property 2.

Property 2. The low-complexity algorithm based on the CP method is convergent.
Proof. In the kth iteration of the proposed algorithm, we first compute F [k] r with the given F , which means the objective value is descendent. Consequently, the objective value of the original problem monotonically decreases and is lower-bounded by zero, which verifies the convergence of Algorithm 2.

Numerical Results and Discussion
In order to analyze the performance of the proposed algorithms, the following simulations are conducted. Fifty random Rayleigh fading channels are generated, and the pathloss exponent is set to 2. The variances of noises are assumed as σ 2 r = σ 2 S i = σ 2 , the transmit powers are set as p i = 18E s and p r = 12E s , and the signal noise ratio (SNR) is calculated from SNR = 10log 10 (E s /σ 2 ), where E s is the power of signal. Meanwhile, the energy-harvesting requirement e i = 0.1. N = 4, M = 2 and N = 2, M = 2 are both considered, and the data stream S = 2. Moreover, the carrier frequency of the system is given by f c = 5 GHz. Four schemes are simulated: 1. The unaided scheme, which means that the beamformers are set as initial matrices; 2. the proposed FPP-SCA scheme; 3. the proposed low-complexity scheme; and 4. the semidefinite relaxation (SDR) scheme [46] used in the previous literature. In order to show the impact of noise, the impact of different values of β, and the number of antennas, we do the corresponding simulations. Figure 2 and Table 1 show the performance under different β for the proposed FPP-SCA scheme and the Low-Complexity scheme. From the simulation results, obviously, a larger β leads to a higher system performance for both schemes, since more signals can be used for decoding the information in the receiver shown in Equation (6). In order to make the comparision with the existing works, we choose to use β = 0.5.
The convergence property of different schemes is evaluated in Figure 3, where the total-MSE is plotted versus the iterations ranging from 0 to 50 in SNR = 5 dB and SNR = 20 dB when N = M = 2. From Figure 3, as the increment in the number of iterations, the FPP-SCA scheme always converges slower and requires more iterations for a convergence as SNR increases than the proposed low-complexity one. Furthermore, the FPP-SCA scheme exhibits a better performance than the low-complexity one for different SNRs when the curve converges. Meanwhile, comparing the FPP-SCA and conventional SDR scheme, we can find that the SDR scheme always converges slower than the FPP-SCA one and that the MSE of it is always higher than that of the FPP-SCA one (e.g., 2.07 vs. 2.04 for SNR = 5 dB and 0.50 vs. 0.48 for SNR = 20 dB) under 50 iterations, which implies an advantage of the proposed FPP-SCA scheme. It can be claimed that the proposed FPP-SCA scheme has a lower level of total-MSE while it has a higher iteration complexity than the low-complexity counterpart.
The performance comparison of different schemes is indicated in Figure 4, Tables 2 and 3, where in Figure    Accordingly, in comparison to the results in Figure 4 and the two tables, we can see that for the proposed FPP-SCA algorithm, the performance is always higher than the SDR-based one for different antennas and SNRs, and we can make the conclusion that our proposed FPP-SCA-based scheme performs better than the traditional SDR-based scheme.
More intriguingly, when N = M, the low-complexity scheme achieves a comparable performance to that of the FPP-SCA one and yields a better performance than that of N > M (e.g., 0.0 vs. 0.04 for SNR = 30 dB shown in Tables 2 and 3). Combined with the low complexity of the low-complexity scheme, it is more applicable than the FPP-SCA one in the N = M case. However, when N > M, in comparison to the FPP-SCA scheme, the performance of the low-complexity scheme is a bit worse owing to the influence of the enhancive diversity gain.
In summary, when the number of antennas at the relay node and source nodes are different, it is more beneficial to choose the FPP-SCA scheme.
From Figure 5 and Table 4, we can see that the performance increases with the number of antennas for both schemes under M = N or N > M. When the number of antennas increases, more antennas can be used to suppress multipath fading with antenna diversity, to increase the system capacity, and to improve the performance. Considering the cost of computing of the low-complexity scheme and comparing the existing work proposed in Reference [35], we choose to use M = N = 2 and M = 2, N = 4 for both schemes. In detail, in Table 4  In order to verify the advantage of the proposed network, we make the comparison of our network and the existing BWSN proposed in Reference [35]. In Reference [35], a joint source and relay design for MIMO two-way relay networks with SWIPT considering a perfect CSI is proposed. In the network, the sources are equipped with PS receivers. The comparison is implemented under the same parameters for the two systems, and the results are as follows.
According to the results shown in Figure 6, it can be observed that the performance of the same algorithm based on the proposed system is preferred to that based on the BWSN, which verify the superiority of the proposed system.

Conclusions
In this paper, we have investigated the joint optimization problem for source and relay beamforming and source receiving in a MIMO FD BWSN SWIPT system. In terms of the problem, two iterative algorithms based on FPP-SCA and low-complexity diagonalizing designs which minimize the total-MSE subjected to the relay-and-source-transmitted power and energy-harvested constraints are proposed. The simulation results demonstrate that the low-complexity scheme always converges faster than the FPP-SCA based one, while the FPP-SCA-based scheme achieves a lower BER compared with the work of the low-complexity scheme. Moreover, when N = M, the performance of the low-complexity scheme yields better than that of N > M. In further works, we will analyze the system performance for multiple users and the interference suppression in the FD network scenario, where a large number of nodes are involved, and a discussion on the optimization scheme under the imperfect SCI will be developed.
Author Contributions: Z.W., D.L., and X.L. conceived the main proposal of the system modeling and derived the analysis and numerical simulation of the proposed schemes. D.L. and X.L. wrote the manuscript. Z.W. provided considerable comments and the technique review of the proposed scheme.
Funding: This research was funded by the Beijing Municipal Natural Science Foundation (4172024).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: