{'fraction_non_alphanumeric': 0.056366903160416475, 'fraction_numerical': 0.019165207776651618, 'mean_word_length': 4.025816161148414, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-Pair Two-Way Relay Network with Harvest-Then-Transmit Users: Resolving Pairwise Uplink-Downlink Coupling

While two-way relaying is a promising way to enhance the spectral efficiency of wireless networks, the imbalance of relay-user distances may lead to excessive wireless power at the nearby-users. To exploit the excessive power, the recently proposed harvest-then-transmit technique can be applied. However, it is well-known that harvest-then-transmit introduces uplink-downlink coupling for a user. Together with the co-dependent relationship between paired users and interference among multiple user pairs, wirelessly powered two-way relay network suffers from the unique pairwise uplink-downlink coupling, and the joint uplink-downlink network design is nontrivial. To this end, for the one pair users case, we show that a global optimal solution can be obtained. For the general case of multi-pair users, based on the rank-constrained difference of convex program, a convergence guaranteed iterative algorithm with an efficient initialization is proposed. Furthermore, a lower bound to the performance of the optimal solution is derived by introducing virtual receivers at relay. Numerical results on total transmit power show that the proposed algorithm achieves a transmit power value close to the lower bound.


I. INTRODUCTION A. Background and Motivation
Wireless power transfer (WPT) is a groundbreaking technique which can prolong the lifetime of battery and feed power to a device when wired charging is inconvenient [1], [2]. However, a critical challenge for WPT is the high propagation path-loss during energy transmission [3]. Fortunately, the beamforming gain in multiple-input-multiple-output (MIMO) systems [4] offers a viable option for mitigating such problem [5]. Therefore, WPT combined with MIMO systems have been a focus lately. Examples include multi-user multipleinput-single-output channel [6], MIMO broadcast channel [7], and distributed antenna network [8]. In all these systems, the interference takes on a dual role: it is detrimental to information decoders, but friendly to energy harvesters.
In the context of WPT, the harvest-then-transmit protocol is recently proposed in [9], which enables energy constrained terminals to transmit data, and ignites the researches on wirelessly powered communication networks (WPCNs) [9]- [11]. In WPCNs, the access point supplies power through the downlink WPT, and the harvested energy at users supports subsequent uplink transmission. This mechanism reveals its uplink-downlink coupled nature, which is apparent in existing works on WPCNs. In particular in [9], the throughput versus time allocation of uplink-downlink transmission is analyzed in WPCN. Meanwhile, the joint uplink-downlink beamforming design in multi-antenna WPCN is proposed in [10].
However, the current research on WPCNs is still at its infancy, focusing on basic point-to-point systems. This paper will take a step further to discuss the wirelessly powered multipair two-way relay (TWR) network, which is a natural generalization of many basic relay systems [12], [13]. Wirelessly powered TWR has extensive applications in sensor networks, medical electronics, smart homes and wearable computations etc., where information exchange between energy-constrained devices is often required. Nonetheless, wirelessly powered TWR is fundamentally different from traditional TWR [12]- [15] or other WPT systems [5]- [10] due to the pairwise uplink-downlink coupling. More specifically, in wirelessly powered TWR, a single transmit beamformer at the relay simultaneously controls a pair of users' downlink information and power transfer. Further adding to the fact that the uplink signal-to-interference-plus-noise ratios (SINRs) of a pair of users are determined by a single receive beamformer, and possibly the power harvested in the previous downlink phase, wirelessly powered TWR suffers from the strong coupling of four transmission links for each pair of users.
Under such pairwise uplink-downlink coupling, an adjustment appearing beneficial to a user may backfire. For example, consider aligning the downlink transmit beamformer towards a user U 1 and away from its paired user U 2 . While this seems to benefit the power transfer to U 1 , the subsequent uplink transmit power of U 2 may decrease due to its reduced harvested energy. Then, to maintain the uplink data-rate of U 2 , the receive beamformer needs to be aligned towards U 2 and away from U 1 . This will in turn deteriorate the uplink datarate of U 1 and the transmit power of U 1 needs to be increased, which offsets the benefit brought by the increased harvested power at U 1 . Consequently, it is crucial to strike a balance where wirelessly powered TWR is the most efficient.

B. Technical Challenges and Contributions
While finding the balance can be cast as a beamforming design problem subject to data-rate quality-of-service (QoS) constraints, the pairwise uplink-downlink coupling makes this problem nonconvex over each beamformer. Therefore, convex programming [5]- [7] or block coordinate descent method [10], [11] will no longer apply. On the other hand, to resolve the nonconvexity, a common way is to reformulate the problem into a difference of convex (DC) program. However, obtaining an exact DC reformulation is a usual challenge. By analyzing the stationary points, this paper proposes a simplification procedure without approximation. The simplified problem is equivalently reformulated as a DC program using rank relaxation, provided that there exists a rank guarantee of the solution. Unfortunately, the considered problem involves a single beamformer serving a group of users, and currently there is no rank guarantee of the solution, even for its special case of TWR without energy harvesting [13]. To this end, we establish a rank-one guarantee for the one pair users case, and a rank-two guarantee for the general multi-pair users case. Interestingly, the rank-two guarantee result also provides a general property for other beamforming design problems with single beamformer serving multiple users, such as those in [13] and [16].
With the transformed DC program, it can be solved by an iterative algorithm with a feasible starting point. Nonetheless, in the multi-pair users case, finding a feasible point is nontrivial when the number of antennas at relay is not sufficient (e.g., when smaller than the number of users). While the conventional l 1 -norm regularization is a popular initialization method, it involves multiple conic programs (CPs). On the contrary, based on the structure of TWR, this paper proposes a CP-free initialization with much lower complexity. Finally, a lower bound on the performance of the optimal solution is derived by introducing virtual receivers at relay, and an analytical solution is obtained for the case of massive antenna array at relay. Numerical results show that the proposed iterative algorithm in multi-pair users case achieves close performance to the lower bound.

C. Organization and Notations
The rest of the paper is organized as follows. In Section II, the system model is described in detail for the uplink and downlink transmission. In Section III, the beamforming design problem is formulated, and transformed to an equivalent problem for subsequent processing. In Section IV, the global optimal solution for the one-pair users case is derived in part A, and the convergence guaranteed iterative algorithm is presented for the multi-pair users case in part B. A lower bound to the performance of the optimal solution and the analytical solution for the case of massive antenna array at relay are derived in Section V. Simulation results are presented in Section VI. Finally, conclusions are drawn in Section VII. Notation. Italic letters, simple bold letters, and capital bold letters represent scalars, vectors, and matrixes, respectively. The operators Tr(·), (·) T , (·) H , Rank(·), (·) −1 take the trace, transpose, Hermitian, rank, and inverse of a matrix, respectively, while vec(·) is the matrix vectorization operator. Symbol I N represents an N × N identity matrix. Symbol 1 N represents an N × 1 vector with all the elements equal to 1.
Finally, E(·) represents the expectation of a random variable.

II. SYSTEM MODEL
In this paper, we consider a multi-pair TWR network consisting of a relay station with N antennas, and 2K singleantenna users. As shown in Fig. 1, the 2K users are paired into K groups, and users in each group intend to exchange information with each other through the relay. Without loss of generality, users in the k th group are denoted as user (1, k) and (2, k), where k ∈ K = {1, 2, ..., K}. Compute-and-forward network coding with lattice codes is applied in the system as it outperforms amplify-and-forward scheme [14], and has a lower complexity than decode-and-forward scheme [15]. The system operates in the half-duplex mode, and the transmission consists of two phases: uplink phase and downlink phase, each with time duration of M symbols. Below, we give the details of each transmission phase.

A. Uplink Phase
In uplink phase, all the users transmit data symbols to the relay simultaneously, with the (i, k) th user symbol being The x i,k is generated from lattice codes [14] and the details are given in Appendix A. The received signal Y ∈ C N ×M at the relay is given by where h i,k ∈ C N ×1 denotes the uplink channel vector, and N ∈ C N ×M is the Gaussian noise at relay with E[vec(N)vec(N) H ] = σ 2 r I MN . In order to separate the signal component of the k th user pair, a receive beamforming vector be the useful signal for the k th user pair, and the remaining part is the interference, we can express the uplink SINR of the (i, k) th user as Furthermore, applying the results from [14,Theorem 3], the uplink achievable rate R UL i,k from the (i, k) th user to the relay can be computed to be where [x] + = max(x, 0), and the factor 1 2 is due to the fact that two transmission phases are involved for single symbol transmission.

B. Downlink Phase with WPT
In the downlink phase, by using w H k Y, the relay generates lattice symbols s T k ∈ C 1×M with power 1 M E(||s k || 2 ) = p k (details documented in Appendix A). Then the relay transmits s k to the k th user pair through the corresponding transmit beamforming vector v k ∈ C N ×1 with ||v k || = 1. Accordingly, the received signal r T i,k ∈ C 1×M at the user (i, k) is where g H i,k ∈ C 1×N is the downlink channel vector from the relay to the (i, k) th user, and n T i,k ∈ C 1×M is the Gaussian noise at the (i, k) th user with E[n i,k n H i,k ] = σ 2 u I M . The received signal (4) at the (i, k) th user in the downlink is further split into two branches, one for the information decoder and the other for the energy harvester.
At the information decoder side, the signal is given by where β i,k ∈ (0, 1] is the splitting factor, and z T i,k ∈ C 1×M is Gaussian noise introduced by the power splitter, with In (5), the first term is the desired network coded signal, the second term is the inter-pair interference, and the last two terms are Gaussian noises. Based on this interpretation of (5), the downlink SINR for the (i, k) th user is Then the downlink achievable rate R DL i,k from the relay to the (i, k) th user is On the other hand, the average harvested power from the wireless signals at user (i, k) can be expressed as η(1 − β i,k )E[||r i,k || 2 ]/M , where 0 < η < 1 is the power conversion efficiency. Based on (4), it can be further expressed as

III. PROBLEM FORMULATION AND TRANSFORMATION
In the considered network, the design variables that can be controlled are the relay transmit-receive beamformers {v k , w k }, relay power allocations {p k }, users' transmit powers {q i,k }, and users' power splitting ratios {β i,k }. Since a basic QoS requirement in a communication system is the guaranteed transmission rate [6], our aim is to provide reliable communication for all the users at their required data-rates. In particular, assuming that the data-rate requirement from the (i, k) th user is R i,k > 0, then the channel from the (i, k) th user to the relay, and that from relay to its paired (3 − i, k) th user 1 should both have achievable rates larger than R i,k , i.e., On the other hand, since the uplink achievable rate R UL i,k in (3) depends on the user transmit power q i,k , which is in turn harvested from the downlink wireless signal, we must have the following power constraint 2 : where E i,k is the local power per symbol-time at the (i, k) th user, and p c is the circuit power consumption per symbol-time (the coefficient 2 is due to the two phases of transmission).
Having the QoS and power harvesting requirements satisfied, it is crucial to minimize the total transmit power at relay and users 3 because saving power translates to cost reduction and environmental benefits. As a result, by accounting for all the factors mentioned above, an optimization problem can be formulated as: where the first and second constraints are the uplink and downlink data-rate constraints, respectively. This optimization 1 The index 3 − i is used to represent the paired user of the i th user due to i ∈ {1, 2}. 2 To incorporate the minimum power threshold ǫ for activating the energy harvesting circuit inside a user terminal [33], we could add an additional Since this constraint can be reformulated into a convex form, it will not change the subsequent derivation and algorithm. 3 The user transmit powers are minimized to save energy for the battery.
problem will be solved at the relay. Then the transmit powers and power splitting ratios will be sent to users. Unfortunately, problem P1 is extremely difficult since even when there is no uplink phase and energy harvesting, the resultant multigroup multicast problem is NP-hard [13], [16]. The first challenge of solving P1 is the term in the first constraint, which makes the left hand side nonlinear with respect to q i,k |w H k h i,k | 2 . One way to get around this is to consider high SINR regime and set the term to zero as in [15,. However, this would yield a conservative solution. In contrast, we can establish the following property.
Proof. See Appendix B.
Putting (9) into P1 and defining α i,k : , the first constraint of P1 can be replaced by without changing the solution of P1 [18]. Now, the constraint (10) is linear with respect to q i,k |w H k h i,k | 2 . The second challenge of solving P1 is the nonlinearly coupled variables of v k , w k , β i,k , p k and q i,k in the first three constraints. While other WPT systems [5]- [7], [10] also suffer from the coupled variables, the problems in those systems are convex or "hidden convex" (appear nonconvex but can be reformulated as convex problems) over each variable. On the contrary, due to the pairwise uplink-downlink coupling, the feasible set of P1 is nonconvex over v k and w k even with other variables fixed. Therefore, the convex programming [5]- [7] or the block coordinate decent algorithm [10] is not applicable to P1. To this end, we introduce new variables: and apply the linear matrix inequalities (LMIs) [17] to reformulate the constraints. Then, based on the procedure given in Appendix C, the problem P1 is equivalent to the problem P2, where µ i,k , I j,l,k are slack variables, Θ i,k = g i,k g H i,k , and θ i,k = 2 2R 3−i,k −1. Problem P2 is well-structured with all the constraints being convex except for the left part of (12a)(12e), and the rank constraints in (12f). These remaining difficulties will be addressed in the following section.

IV. JOINT UPLINK-DOWLINK NETWORK DESIGN
In this section, the global optimal solution is first derived for the special case K = 1, and then we step further to the general case K ≥ 2.

A. One Pair Case: Optimal Solution
When K = 1, the network reduces to a three-node TWR. Since there is no l = k, the term I j,l,k can be dropped, and (12a) becomes where the subscript k is also dropped to simplify the notation in this subsection. With a similar proof to Property 1, it can be shown that any stationary point of P2 with K = 1 must activate (13) and therefore the problem P2 with K = 1 is equivalent to Now to reduce the problem dimension of P3, the following property can be established.
Property 2. Denoting w * and V * as the optimal solution of w and V in P3, we have w * ∈ span{h 1 , h 2 } and the eigenvector of V * in span{g 1 , g 2 }.
Proof. See Appendix D.
Based on Property 2 and [19, , the receive beamformer w can be parameterized by 0 ≤ γ ≤ 1 as On the other hand, following the procedure in [32, Section III], V can be expressed as with A ∈ C 2×2 0. Putting (15) and (16) into P3, we have the following equivalent problem: where In P4, the remaining nonconvex parts are the terms containing γ and the rank constraint in (17d). While the nonconvexity of γ can be resolved by 1-D search, the rank constraint can be dropped by applying the semidefinite relaxation (SDR) [20], and the following proposition shows that the relaxation does not affect the optimality. Proposition 1 guarantees that an optimal rank-one solution to the SDR problem of P4 exists. However, there may be an alternative solution with a higher rank. If we obtain such a solution, the rank reduction procedure in [21] can be applied to obtain the rank-one solution. Therefore, the problem P4 is equivalent to its SDR problem, which can be optimally solved by an iterative procedure with 1-D search of γ over [0, 1], and semidefinite programming (SDP) for each fixed γ. Denoting the optimal solution to P4 as A * , γ * , {β * i , µ * i }, then the optimal solution to P3 is V * = GA * G H and w * = w(γ * ). Therefore, the optimal v * , w * , p * , {q * i , β * i } of P1 with K = 1 can be recovered accordingly.
In terms of computation complexity of the proposed algorithm, calculating . On the other hand, solving the SDR problem of P4 is dominated by the 2-dimensional SDP cone in (17d), which requires complexity O(2 3.5 t) [17], with t being the number of iterations for one-dimensional search. As a result, the total computational complexity would be O(2N 2 + 2 3.5 t).
Notice that the problem P1 with K = 1 has also been discussed in [11], where a block coordinate descent method with each iteration solving an SDP problem was proposed. However, such a method only converges to a suboptimal solution for the nonconvex problem P1, in contrast to the global optimal solution derived in this paper. Moreover, since the problem dimension in [11] is N , the complexity of the method in [11] is O(2sN 3.5 ), with s being the number of iterations, and is larger than that of the proposed method in this paper.

B. General Case of K ≥ 2
For K ≥ 2, a critical issue of P2 is the rank constraints in (12f). In the following, by analyzing the Lagrangian and the Karush-Kuhn-Tucker (KKT) condition, a proposition is established for the rank-relaxed problem of P2.
, k ∈ K} who do not need energy harvesting. First we construct the augmented Lagrangian with respect to V k and β i,k for the rank relaxed problem of P2 as where are Lagrangian multipliers. According to [18], the primal variables and the Lagrangian multipliers of any stationary point should together satisfy the KKT condition ∂L/∂V k = 0 for any k and ∂L/∂β i,k = 0 for any (i, k) / ∈ G, which can be explicitly expressed as where (20) and Notice that µ ⋄ i,k = 0 for any (i, k) / ∈ G. Now, by substituting (21) into the definition of Υ, we have Below we will show that Υ is of full rank, i.e., all of its eigenvalues are nonzero. This can be proved by contradiction. Assuming that there exists a vector a = 0 such that a H Υa = 0, then by left multiplying a H and right multiplying a to (19), we have where the inequality is due to Ξ ⋄ k 0. Using a H Υa = 0, (23) further simplifies to . As a consequence, by left multiplying a H and right multiplying a to (22), we have a H Υa = a H Ia = ||a|| 2 . Before (23), we assumed that a H Υa = 0, which leads to a = 0. This contradicts to a = 0, and therefore Υ is of full rank.
Using Rank(Υ) = N and (19), and due to On the other hand, the complementary slackness condition for Ξ ⋄ Putting Rank(Ξ ⋄ k ) ≥ N − 2 and Rank(Ξ ⋄ k V ⋄ k ) = 0 into the above inequality immediately contributes to Rank(V ⋄ k ) ≤ 2 and the proof is completed.
Notice that by setting β i,k = 1 for all i, k, Proposition 2 applies to many systems without WPT, such as those in [13] and [16]. Specifically, Proposition 2 can be directly applied to the power minimization problem in multi-pair TWR [13]. Furthermore, for the multigroup multicast problem in [16], with the number of users in each group equal to 2, Proposition 2 holds for the QoS problem, and with the transformation in [16,Claim 3], it also applies to the max-min fairness (MMF) problem. While [13], [16] only claim that the rank-one solution is not guaranteed in their problems, Proposition 2 reveals that the solution must be rank two or less. Finally, when the number of users in each group in the multigroup multicast problem [16] is larger than 2, with a similar proof as in Proposition 2, it can be shown that Rank(V k ) ≤ G k always holds for the QoS and MMF problems, where G k is the number of users in the k th group.
Proposition 2 indicates that the rank of V ⋄ k is either 1 or 2. When V ⋄ k is rank-one, the rank relaxation is tight. On the other hand, when V ⋄ k is rank-two, transmission in blocks of Alamouti codes will also guarantee no performance loss [23]. Therefore, dropping the rank constraint in (12f) will not affect the solution of P2.
Problem P2 without the rank constraints in (12f) is a DC program, since the left hand side of (12a) and (12e) are difference of convex functions. DC programs received a lot of attention lately (e.g., in medical imaging, financial engineering, and machine learning [24]) due to the fact that a local optimal solution can be found in polynomial time using the inner approximation method [24]- [28]. Inner approximation is an iterative method where the DC part is replaced by a convex upper bound expanded around the last round solution. In particular, define the second term on the left hand side of (12a) as Assuming that the solution at the n th iteration is given by {w and the following property can be established.

Property 3. The functionΦ
[n] i,k satisfies the following: i,k ), and (iii) Proof. See Appendix E.
In addition to being linear in w k and ξ i,k , from Property 3, we can see thatΦ [n] i,k has the same value and gradient as Φ i,k at point {w i,k }. More importantly, the first part of Property 3 indicates that if we replace Φ i,k withΦ [n] i,k , the feasible set becomes smaller, thus the solution of such a problem would be a feasible solution to the rank-relaxed problem of P2. Following a similar proof to Property 3, the term −µ 2 i,k in constraint (12e) can be conservatively replaced by −2µ i,k is the solution at the n th iteration. With the above observation, the following problem is considered at the (n + 1) th iteration Now the problem P2 − R[n + 1] is an SDP problem which can be efficiently solved by CVX, a Matlab-based software package for convex optimization [22]. Denote its optimal solution as , and the process repeats with solving the problem P2 − R[n + 2]. According to Property 3 and [27, Theorem 1], the iterative algorithm is convergent. Furthermore, the converged point { V k , β i,k , w k , ξ i,k } is guaranteed to be a stationary point for the rank-relaxed problem of P2 [24]- [28].
In terms of computation effort, solving problem P2−R[n+ 1] is dominated by K N -dimension SDP cones, which requires complexity O K 3.5 N 2.5 + K 2.5 N 3.5 [17]. Therefore, the total complexity is O t(K 3.5 N 2.5 + K 2.5 N 3.5 ) , where t is the number of iterations needed for the inner approximation to converge, which is very small as shown by simulations.
With the obtained solution, the next stage is to recover the beamformer v k of P1. More specifically, when V k is rankone, the beamformer v k can be recovered by the rank-one decomposition of V k . On the other hand, when V k is ranktwo, by allocating power p k = Tr( V k ) at relay for the k th downlink data stream, and applying rank-two decomposition, we obtain V k / p k = [ f 1,k f 2,k ][ f 1,k f 2,k ] H . Then, the relay transmits [ f 1,k f 2,k ]S k (m) at each time duration for two symbols, where S k (m) ∈ C 2×2 is the k th downlink data stream grouped into blocks of Alamouti code: , with m = 1, ..., M/2. As shown in [23], such a transmission incurs no loss when used together with a rank-two beamformer. Consequently, the proposed algorithm is guaranteed to achieve at least a stationary point of P1.
Finally, for the inner approximation algorithm, finding a good feasible starting point is of large importance [24]. For N ≥ 2K − 1, since the degree of freedom at relay is sufficient to support the number of users, P2 is always feasible, and a simple initialization is fixing the transmit and receive beamformer using pairwise zero-forcing criterion [12, Section IV]. Then P2 is an SDP with low dimension. However, when N < 2K − 1, since the first constraint of P2 involves DC functions, it is generally hard to find a feasible initial point or provide a feasibility condition [24]- [26]. A traditional method is using l 1 -norm regularization to pursuit a feasible point [24], [25], which involves multiple conic programs (CPs). Below, by exploiting the structure of problem P2, we propose a CP-free initialization.
In particular, from the first two constraints of problem P2, an initial point {w To find such a feasible solution, we consider the following problem with the objective function equal to the left hand side of the above inequality where P is a sufficiently large control parameter. Then a feasible initial point can be found by alternatively updating w k and ξ i,k in (28) until its objective function is greater than or equal to 1. An advantage of such a block coordinate descent method is that closed-form optimal solutions for each subproblem of (28) can be obtained as proved in Appendix F. After obtain- i,k + 2p c − 2E i,k ) + according to the constraint (12e) of P2. As each step is in a closed-form, the complexity of the proposed initialization is inconsequential. Therefore, the proposed CP-free initialization also represents a low-complexity solution for solving P2.
Notice that even if P2 is feasible, currently there is no method guaranteed to find a feasible solution due to the nature of DC programming problems [24]. While the proposed CPfree method cannot guarantee a feasible initialization either, it has a very high probability of finding a feasible initialization since the block coordinate descent method converges to at least a stationary point of (28) [36, Section 3.1]. Moreover, the probability of finding a feasible initialization can be further increased using multiple starters of {ξ i,k } [25].

A. Lower Bound of P1
In order to assess the performance of the proposed solution in Section IV-B, in this section, a design which provides a lower bound to the transmit power of the optimal solution of P1 is derived as a benchmark. In particular, consider the following problem P1 − R: It can be seen that if we restrict b 1,k = b 2,k = w k and Rank(V k ) = 1 in P1 − R, it reduces to P1. Therefore, the feasible set of P1 − R is larger than that of P1, and P1 − R is a relaxed version of P1. The insight behind such a relaxation is that {b i,k } represents virtual receivers at relay, which attenuates the pairwise uplink-downlink coupling. In the following, we will show that P1 − R can be optimally solved in two steps.
From P1−R, it can be seen that {b i,k , q i,k } is only involved in (29a), (29c) and the first part of (29d). Observing from (29c) that smaller q i,k loosens the constraints on V k , which helps in reducing the objective value, the optimal {q * i,k , b * i,k } of P1−R can be therefore obtained from solving P5 : min Then by putting the optimal {q * i,k , b * i,k } of P5 into P1 − R, the problem can be simplified and the optimal {V * k , β * i,k } can be found from the simplified problem.
To solve problem P5, which is a multi-criterion optimization, consider the following iteration: and the following proposition can be established.
i,k is convergent, and the limit point is the optimal solution of P5.
Proof. First we prove ω [n] i,k is monotonically increasing by induction. Since ω [0] i,k = 0 ≤ ω [1] i,k , we assume that ω where the inequality in the second line is due to i,k = arg min which can be obtained from (31) and [29,Corollary 2], the inequality in the third line is due to ω i,k , and the equality in the last line is due to (32). Thus ω [n] i,k is monotonically increasing.
Then we prove ω i,k is upper bounded. Specifically, we will show that ω i,k ≤ q i,k for some n ≥ 0. Using (32), we have where the inequality in the second line is due to (34), the inequality in the third line is due to ω i,k ≤ q i,k , and the inequality in the last line is due to the constraint in (30). Therefore as long as P5 has a feasible solution, the sequence ω i,k is upper bounded. Lastly, using the above two results, we immediately have ω [n] i,k is convergent [30]. Furthermore, by putting ω i,k into (32), the limit point ω which implies that {ω [∞] i,k , z [∞] i,k } is a feasible solution for P5. On the other hand, since we obtained ω i,k } is therefore optimal for P5.
Based on Proposition 3, the optimal solution of P5 is i,k }. Then by substituting them into P1 − R, the problem P1 − R reduces to The above problem is an SDP which can be optimally solved by CVX. Since we have shown that P1−R can be equivalently transformed into a two-step optimization, and we solve each step optimally, the obtained solution is optimal for P1 − R, and would be a lower bound for P1.

B. Large N Analysis of P1
In this subsection, wirelessly powered multi-pair TWR with very large number of antennas at relay, i.e., N → ∞, is analyzed. Specifically, when N is very large, user channels will be asymptotically orthogonal, and the interference vanishes [34,Proposition 1]. Therefore, problem P1 asymptotically becomes P6 : min Notice that there is no need to apply network coding since all the channels are asymptotically orthogonal. In this case, the uplink and downlink SINR targets become θ 3−i,k = 2 2R i,k −1 and θ i,k = 2 2R 3−i,k − 1, respectively. Observing from P6 that minimizing p i,k and q i,k requires matching the beamformers to their corresponding channels, the optimal w * i,k = h i,k /||h i,k || and v * i,k = g i,k /||g i,k ||. Putting w * i,k into the first constraint of P6, it reduces into q i,k ||h i,k || 2 ≥ θ 3−i,k σ 2 r . Since the effect of fast fading vanishes when N is large [34], we have ||h i,k || 2 = N ̺ i,k , where ̺ i,k is the large-scale fading of h i,k , and thus This result is consistent with the fact that information beamforming gain is proportional to N in large-scale antenna systems [34]. Therefore, given certain data-rate QoS, the users' transmit powers can be decreased proportionally to 1/N while still satisfying the QoS constraints.
To determine the network transmit power, we also need to solve for p * i,k in P6. Based on the procedure given in Appendix G, its optimal solution is given by From the analytical solution of p * i,k , we can see that the relay transmit power can be decreased proportionally to 1/N as well. This result is consistent with the fact that energy beamforming gain is proportional to N in large-scale antenna systems [35]. Based on the scaling law of relay and user transmit powers, we can conclude that when N is large, the network transmit power can be reduced proportionally to 1/N while still satisfying the data-rate QoS constraints and energy harvesting constraints.

VI. SIMULATION RESULTS AND DISCUSSIONS
This section provides simulation results to demonstrate the performance of the proposed schemes. In particular, the distance-dependent pathloss model of the (i, k) th user ̺ i,k = ̺ 0 ·( d i,k d0 ) −α is adopted [10], where d i,k is the distance from the (i, k) th user to the relay, ̺ 0 = 10 −3 , d 0 = 1m is the reference distance, and α is the pathloss exponent set to be 2.7 [8], [9]. In the simulations, d i,k ∼ U(1, 10) in meter, where U represents the uniform distribution, and g i,k , h i,k are generated according to CN (0, ̺ i,k I). It is assumed that power conversion efficiency η = 0.8 [33], and noise power σ 2 r = σ 2 u = σ 2 z = −60dBm [9]. The circuit power consumption is set to p c = 10dBm [31], and E i,k ∼ U(9.5, 13.0) in dBm. For the simulations of transmit power versus data-rate QoS, the same data-rate targets R i,k = R are requested by all users [16], while for other simulations R i,k ∼ U(0, 2) in bps/Hz. Each point in the figures is obtained by averaging over 100 simulation runs, with independent channels in each run.
First we consider the one pair case K = 1 with N = 8. Here, four schemes are compared: the optimal solution, the SDP solution from [32], the iterative solution from [11], and the solution with β 1 = β 2 = 0.5. As is shown in Fig.  2, the solution in Section IV-A of this paper achieves the lowest transmit power over a wide range of data-rate QoS. Compared to the suboptimal solutions in [32] and [11], the proposed scheme has an average advantage of 0.3dB and 2.1dB, respectively. It is also observed that fixing power splitting ratio leads to more than 3dB loss compared to the proposed optimal solution, indicating that the power splitting ratios should be jointly optimized with the beamformers.  Next we consider the case of three pairs of users K = 3 with N = 12. To verify the convergence of the proposed iterative method in Section IV-B, Fig. 3 shows the total transmit power versus number of iterations. It can be seen that different initializations for the proposed iterative algorithm converge to the same value, which is within 0.2dB from the lower bound. Furthermore, among the initializations, the proposed CP-free initialization converges the fastest and stabilizes after 3 iterations, indicating that the complexity of the iterative algorithm can be moderate with the CP-free initialization. For the l 1 -norm regularization, due to its random picking of starting points, it requires more than 10 iterations to converge. For the ZF beamforming initialization, since it completely eliminates the inter-pair interference and forces the harvested energy at users to come from useful signals only, the initial transmit power is high. But it also shows fast convergence, making it a compelling initialization when N ≥ 2K − 1. For the rest of this section, the proposed CP-free initialization will be used and the iterative algorithm stops after 3 iterations.
To further demonstrate the performance of the proposed iterative algorithm for solving P2, Fig. 4 and Fig. 5 show the total transmit power versus the data-rate QoS and the noise power, respectively. Apart from the ZF transmit-receive  beamforming scheme, we also simulated the ZF receive beamforming scheme [10, Section V], which optimizes the transmit beamformer in P2 and provides more flexibility than the ZF transmit-receive beamforming scheme. It can be observed from both figures that the proposed algorithm significantly outperforms other schemes and approaches the lower bound very tightly over a wide range of data-rate QoS and noise power. This again verifies the excellent performance of the proposed method. Interestingly, from Fig. 5, the performance of the suboptimal scheme using ZF receive beamformer also approaches the lower bound when the noise power is extremely small, e.g., smaller than −75dBm. As the scheme using ZF receive beamformer has a low complexity, it represents a promising solution in high SNR regime.
The aforesaid experiments examine the situation with N ≥ 2K − 1. Now we focus on the case of N < 2K − 1, in which the ZF schemes are invalid but the proposed algorithm is still applicable. Specifically, the cases of N = 2K − 2 with K = 5 and K = 7 are simulated. It can be seen from Fig.  6 that although the gap between the proposed algorithm and the lower bound is larger due to the lower degree of freedom, the proposed algorithm still achieves satisfying performance very close to the lower bound. This reflects the high efficiency of the proposed method even in the case of low degree of freedom.

VII. CONCLUSIONS
This paper studied the total transmit power minimization in the multi-pair TWR system with harvest-then-transmit users. Despite the challenge brought by the pairwise uplink-downlink coupling, it was proved that optimality could be achieved in the special case of one pair of users. For the general case with multiple pairs of users, a convergence guaranteed iterative algorithm was proposed, and a lower bound on the performance of the optimal solution was derived. Simulation results demonstrated that the proposed method outperforms other methods, and achieves good performance very close to the derived lower bound.

APPENDIX A COMPUTE-AND-FORWARD WITH LATTICE CODES
Generating x i,k at the (i, k) th user Given q i,k and |w H k h i,k | 2 , we choose K pairs of M × 1 doubly nested lattice [14] such that Λ 1,k ⊆ Λ 2,k ⊆ Λ k with k ∈ K, and the second moments σ 2 (Λ i,k ) = q i,k |w H k h i,k | 2 . For user (i, k), the source data is mapped into c i,k ∈ L i,k = {Λ k mod Λ i,k } and the symbols to be transmitted are where d i,k is a pre-generated random dither vector known to the k th pair of users and relay. Notice that the transmit power of x i,k is 1 M E[||x i,k || 2 ] = q i,k . Generating s k at relay Putting (38) into (1) and applying w k to extract the k th layer signal at relay, we have where the first term is the useful signal, the second term is inter-pair interference, and the last term is noises with power 1 M E[||w H k N|| 2 ] = σ 2 r . Then, the relay computes [14]: where represents the nearest neighbor lattice of c i,k + d i,k in Λ i,k , and Sincem k is the suppressed interference (whitened by random dither, it can be viewed as Gaussian when M → ∞) plus noise, the relay can recover t k by applying lattice decoding to equation (39). Then the relay maps the retrieved t k to symbol s k (t k ) ∈ L r,k , where L r,k is the lattice codebook at relay (not related to L 1,k and L 2,k ), with σ 2 (Λ r,k ) = p k .

APPENDIX B PROOF OF PROPERTY 1
To prove the property, we first show that any stationary point of P1 activates the first constraint for all users. For simplicity, define the left hand side of the first constraint of P1 as Then the Lagrangian of P1 with respect to {q i,k } could be written as [22] with Lagrangian multipliers λ i,k ≥ 0, ρ i,k ≥ 0, ν i,k ≥ 0.
According to the KKT condition ∂L/∂q i,k = 0 for all i, k, we have Since R i,k > 0, the user power satisfies q i,k = 0. From complementary slackness ν i,k q i,k = 0, equation ν i,k = 0 holds. Putting ν i,k = 0 into (41) gives To analyze the value of λ i,k , we need to determine the signs of ∂∆ j,l /∂q i,k for (j, l) = (i, k). More specifically, dividing the numerator and denominator of the first term in (40) by Clearly, for ∆ j,l with (j, l) = (i, k), when q i,k decreases, ∆ j,l will monotonically increase because q i,k only appears in the denominator of ∆ j,l . That is, ∂∆ j,l /∂q i,k ≤ 0 for (j, l) = (i, k). Putting ∂∆ j,l /∂q i,k ≤ 0 for (j, l) = (i, k) into equation (42), we immediately have ∂∆ i,k /∂q i,k ·λ i,k ≥ 1, ∀i, k, which leads to λ i,k = 0. Furthermore, due to the complementary slackness Now, it is clear from (44) that the term inside the parenthesis remains the same for both users pertaining to the k th user pair. Thus, dividing both sides of (44) with i = 1 by that with i = 2 for the k th user pair, the term inside the parenthesis can be cancelled and we get q 1,k |w H k h 1,k | 2 /(q 2,k |w H k h 2,k | 2 ) = 2 2R 1,k /2 2R 2,k , which implies that This completes the proof.

APPENDIX C REFORMULATION OF P1
Based on Property 1 and after some manipulations, the following problem P1 ′ is equivalent to P1 in the sense that they have the same optimal solutions ) q i,k ≥ 0, β i,k ∈ (0, 1], ∀i, k ||v k || = 1, ||w k || = 1, p k ≥ 0, ∀k. (45d) Since the fractional-quadratic function is known to be convex, the coupling between w k and q i,k in (45a) can be resolved by letting q i,k = 1/ξ i,k . This is always possible because α i,k defined in (10) satisfies α i,k > 0 and then from (45a) q i,k = 0 holds. On the other hand, to linearize the quadratic terms of v k in (45b) and (45c), we introduce new variables V k := p k v k v H k 0 with Rank(V k ) = 1 for all k ∈ K. Applying the transformations presented above, the problem P1 ′ is equivalent to the following problem P1 ′′ : ξ i,k > 0, β i,k ∈ (0, 1], ∀i, k, where Θ i,k = g i,k g H i,k . Notice that the constraint ||w k || = 1 in P1 ′ is relaxed into ||w k || ≤ 1, which does not affect our problem since the optimal w * k of P1 ′′ always satisfies ||w * k || = 1. To deal with (46a), we introduce slack variables I j,l,k ≥ |w H k h j,l | 2 /ξ j,l , and (46a) becomes the constraints (12a) and (12b) in P2. Furthermore, the constraint (46b) can be rearranged as an LMI through Schur Complement Lemma [22], and this gives (12c) in P2. Now we focus on (46c), which is non-convex in its current form. However, by introducing a slack variable µ i,k such that the first and second inequalities of (47) can be cast as (12d) and (12e) in P2, respectively. Finally, it is remarkable that the constraints 0 < β i,k ≤ 1 and ξ i,k ≥ 0 in (46d) are implicitly incorporated by the LMI constraints (12b)-(12d) in P2, and thus they can be dropped without changing the problem. After the above procedure, P1 ′′ is equivalently transformed into P2.

APPENDIX D PROOF OF PROPERTY 2
We first address the proof for w * . Assume that w = a 1 h 1 + a 2 h 2 + a 3 h ⊥ with h H ⊥ h i = 0. Consider two solutions of w in P3, i.e., a 1 h 1 + a 2 h 2 + a 3 h ⊥ with a 3 = 0 and ǫ a 1 h 1 + ǫ a 2 h 2 , while other variables are fixed. Based on the norm constraint || a 1 h 1 + a 2 h 2 || 2 + || a 3 h ⊥ || 2 = ǫ 2 || a 1 h 1 + a 2 h 2 || 2 , we must have ǫ > 1. On the other hand, putting the two solutions a 1 h 1 + a 2 h 2 + a 3 h ⊥ and ǫ a 1 h 1 + ǫ a 2 h 2 into the objective function of problem P3, we have This indicates that for any solution with a 3 = 0, we can always find another one satisfying (14c) while giving smaller objective value. Therefore, we must have a * 3 = 0. Next we address the proof for V * . Since Rank(V) = 1, we can assume V = [b 1 g 1 + b 2 g 2 + b 3 g ⊥ ][b 1 g 1 + b 2 g 2 + b 3 g ⊥ ] H with g H ⊥ g i = 0. Then, we have Tr(V) = ||b 1 g 1 + b 2 g 2 || 2 + ||b 3 g ⊥ || 2 . With a similar proof to w * , it can be easily shown that b 3 = 0 would minimize the objective function. Furthermore, due to Θ i g ⊥ = g H ⊥ Θ i = 0, we also have Tr(Θ i V) = Tr(Θ i [b 1 g 1 + b 2 g 2 ][b 1 g 1 + b 2 g 2 ] H ), meaning that b 3 does not affect the constraints (14a) and (14b). This completes the proof.

APPENDIX E PROOF OF PROPERTY 3
To prove part (i), consider the following inequality Multiplying ξ i,k on both sides of the above inequality immediately yieldsΦ [n] i,k (w k , ξ i,k ) ≥ Φ i,k (w k , ξ i,k ). To prove part (ii), we substitute w [n] k and ξ [n] i,k into the definition ofΦ [n] i,k in (26) i,k ).
To prove part (iii), we first calculate the following derivatives: