Intelligent Reflecting Surface Aided MIMO Broadcasting for Simultaneous Wireless Information and Power Transfer

An intelligent reflecting surface (IRS) is invoked for enhancing the energy harvesting performance of a simultaneous wireless information and power transfer (SWIPT) aided system. Specifically, an IRS-assisted SWIPT system is considered, where a multi-antenna aided base station (BS) communicates with several multi-antenna assisted information receivers (IRs), while guaranteeing the energy harvesting requirement of the energy receivers (ERs). To maximize the weighted sum rate (WSR) of IRs, the transmit precoding (TPC) matrices of the BS and passive phase shift matrix of the IRS should be jointly optimized. To tackle this challenging optimization problem, we first adopt the classic block coordinate descent (BCD) algorithm for decoupling the original optimization problem into several subproblems and alternately optimize the TPC matrices and the phase shift matrix. For each subproblem, we provide a low-complexity iterative algorithm, which is guaranteed to converge to the Karush-Kuhn-Tucker (KKT) point of each subproblem. The BCD algorithm is rigorously proved to converge to the KKT point of the original problem. We also conceive a feasibility checking method to study its feasibility. Our extensive simulation results confirm that employing IRSs in SWIPT beneficially enhances the system performance and the proposed BCD algorithm converges rapidly, which is appealing for practical applications.


I. INTRODUCTION
R ECENTLY, intelligent reflecting surface (IRS)-assisted wireless communication has received considerable research attention, since it is capable of supporting costeffective and energy-efficient high data rate communication for next-generation communication systems [1]- [3]. In simple tangible terms, an IRS is composed of a vast number of low-cost and passive reflective components, each of which is capable of imposing a phase change on the signals incident upon them. Thanks to the recent advances in metamaterials [4], it has become feasible to reconfigure the phase shifts in real time. As a result, the phase shifts of all reflective components can be collaboratively adjusted for ensuring that the signals reflected from the IRS can be added constructively or destructively at the receiver in order to beneficially steer the signal component arriving from the base station (BS) for enhancing the desired signal power or alternatively for suppressing the undesired signals, such as interference. In contrast to conventional physical layer techniques that are designed for accommodating the hostile time-varying wireless channels [5], [6], IRSs constitute a new paradigm capable of 'reprogramming' the wireless propagation environment into a more favorable transmission medium. Since the reflective components are passive, they impose a much lower power consumption than conventional relay-aided communication systems relying on active transmission devices. Additionally, no thermal noise is imposed by the IRS, since it directly reflects the incident signals without decoding or amplifying them, which is in contrast to conventional relays. Furthermore, the reflective phase arrays can be fabricated in small size and low weight, which enables them to be easily coated in the buildings' facade, ceilings, walls, etc. Furthermore, as IRS is a complementary device, it can be readily integrated into current wireless networks without modifying the physical layer standardization, making it transparent to the users. To fully exploit the benefits of IRS, the active beamforming at the BS and the passive beamforming at the IRS should be jointly designed. However, the optimization variables are coupled and the joint design leads to a complex optimization problem that is difficult to solve. Some innovative efforts have been devoted to the transceiver design when integrating IRS into various wireless This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ communication systems, including the single-user scenarios of [7]- [11], the downlink multiple-user scenarios of [12]- [15], the physical layer security design of [16]- [21], the mobile edge computing (MEC) networks of [22], multigroup multicast networks of [23] and the multicell multiuser multiple-input multiple-output (MIMO) case in [24]. Concretely, Wu et al. proposed joint active and passive beamforming for a singleuser scenario in [7], where semidefinite relaxation (SDR) was proposed for optimizing the phase shift matrix. However, its complexity is high since the number of optimization variables increases quadratically with the number of phase shifts. Additionally, the Gaussian random approximation employed leads to certain performance loss. To resolve this issue, Yu et al. [8] proposed a pair of efficient algorithms termed as fixed point iteration and manifold optimization techniques, which can guarantee locally optimal solutions. As a further advance, the authors of [9] considered realistic frequencyselective channels. The phase shift design was studied in [10] when only statistical channel state information (CSI) is available. A sophisticated phase shift model was derived in [11], by taking into account a realistic amplitude-phase relationship. For the multiuser case, the authors in [12] considered the total transmit power minimization problem, while guaranteeing the users' signal-to-interference-plus-noise ratio (SINR) constraints. The associated energy efficiency maximization problem was studied in [13] and zero-forcing beamforming was adopted by the BS for simplifying the optimization problem. By contrast, a weighted sum rate (WSR) maximization problem was considered in [14] and the fairness issues were studied in [15]. The authors of [16]- [18] studied the security issues of a single-user case, while the authors of [19]- [21] considered multiple-user scenarios. In [22], the IRS was shown to be beneficial in reducing the latency of MEC networks. In addition, the IRS can help enhance the WSR performance for the multigroup multicast network in [23]. Most recently, we considered an IRS-assisted multicell MIMO communications scenario [24], where we demonstrated that deploying an IRS at the cell edge is also capable of mitigating the adjacentcell interference. Channel state information (CSI) is challenging to obtain in IRS-assisted communication system due to its passive feature. There are some initial efforts to handle this issue such as channel estimation and/or robust transmission design [25]- [28]. Specifically, Huang et al. [25] proposed a deep learning method for efficient online configuration of the phase shifts, where the phase values can be immediately obtained by inputting the user location into the trained deep neural network. A two-stage channel estimation method based on a sparse matrix factorization and a matrix completion was proposed in [26]. A pair of algorithms based on compressed sensing and deep learning were conceived by Taha et al. [27] for tackling the challenging channel estimation issues of IRSassisted systems. Most recently, we first studied the robust beamforming design for IRS-assisted communication systems in [28], where the imperfect channel from an IRS to users was considered and the channel estimation error was assumed to be within a bounded elliptical region.
On the other hand, information transmission enabled simultaneous wireless information and power transfer (SWIPT) is an appealing technique for future energy-hungry Internet-of-Things (IoTs) networks. Specifically, a base station (BS) with constant power supply will transmit wireless signals to a set of devices. Some devices intend to decode the information from the received signal, which are termed as information receivers (IRs), while the others will harvest the signal energy, which are called energy receiver (ER). In [29], Zhang et al. studied the trade-off between the information rate attained and the amount of harvested energy for a single-user MIMO system. In practice, a typical ER such as a humidity sensor requires much higher energy for its operation than that required by IRs. Due to severe channel attenuation, the power received by the ERs is weak, which limits the maximum link-distance of ERs. To mitigate this issue, we propose to deploy an IRS in the vicinity of ERs to provide additional transmission links to support the ERs for enhancing their harvested power as shown in Fig. 1, since there is a paucity of IRS-assisted SWIPT contributions in the literature [30]. Explicitly in [30], the weighted sum power maximization problem was studied by Wu and Zhang, who proved that no dedicated energy-carrying signals were required for an IRS-aided SWIPT system. The SDR method was adopted for solving the optimization problem, which exhibits a high computational complexity as well as imposing a performance degradation due to the associated rank-one extraction. However, this method is not applicable when each user is equipped with multiple antennas. Hence, in this paper we formulate a weighted sum rate (WSR) maximization problem for the IRS-assisted SWIPT MIMO system of Fig. 1, in which an IRS is installed in the vicinity of ERs for compensating the associated power loss, while maximizing the WSR of distant IRs with the aid of passive beamforming.
Against this background, the main contributions of this paper are summarized as follows: 1) We formulate the WSR maximization problem by jointly optimizing the transmit precoding (TPC) matrices of the BS and those of the passive beamforming at the IRS for our IRS-assisted SWIPT MIMO system subject to a non-convex unit-modulus constraint imposed on the phase shifts, while simultaneously satisfying the energy harvesting requirement of the ERs. To the best of our knowledge, this is the first treatise considering the WSR maximization problem of IRS-assisted SWIPT MIMO systems, which is much more challenging than the weighted sum power minimization problem of [30] since the latter can be readily transformed into a convex optimization problem. In contrast to the multicell system of [24], an additional energy harvesting constraint is also imposed in our current study, which further complicates the analysis. Specifically, this constraint is non-convex and the optimization problem may become infeasible. The WSR maximization problem is challenging to solve, since the optimization variables are highly coupled and the data rate expressions of the IRs are complex. To deal with this issue, we first reformulate the original problem into an equivalent form by exploiting the equivalence between the data rate and the weighted minimum meansquare error (WMMSE). Then, an alternating optimization algorithm based on the popular block coordinate descent (BCD) algorithm is proposed for alternately updating the active TPC matrices of the BS and the phase shift matrix of the IRS, which is rigorously proved to converge to the Karush-Kuhn-Tucker (KKT) point of the original optimization problem. 2) For a given phase shift matrix, we then proceed by developing an iterative algorithm based on the successive convex approximation (SCA) method and on the Lagrangian dual decomposition method to derive a nearly closed-form solution for the TPC matrices. A low-complexity bisection search method is proposed for finding the optimal dual variables. The solutions generated by our iterative algorithm are guaranteed to converge to the KKT point of the TPC optimization problem. 3) For the given TPC matrices, we formulate the phase shift optimization problem as a non-convex quadratically constrained quadratic program QCQP) subject to an additional energy harvesting constraint by invoking some further matrix manipulations. Then, a novel iterative algorithm based on the majorization-minimization (MM) algorithm [31] and on the price-based method [32] is developed for solving the QCQP. We strictly prove that the final solution generated by the iterative algorithm is guaranteed to converge to the KKT point of the phase shift optimization problem. 4) The associated feasibility issue is also studied by formulating an alternative optimization problem and an iterative algorithm is proposed for solving this problem. 5) Extensive simulation results are provided for verifying the performance advantages of employing IRS in SWIPT in order to enhance the energy harvesting performance. It is shown that the operating range of the ERs can be dramatically expanded by placing IRSs in the ERs' vicinity. Furthermore, the BCD algorithm converges rapidly, and it is eminently suitable for practical applications. Our simulation results also show that as expected, the path loss exponent substantially affects the system's performance and thus the location of the IRS should be carefully chosen. The remainder of this paper is organized as follows. In Section II, we introduce the IRS-assisted SWIPT system model and our problem formulation. The detailed algorithms used for solving the optimization problem are presented in Section III. The feasibility issues of the original problem are discussed in Section IV, followed by our extensive simulations and discussions in Section V. Finally, our conclusions are provided in Section VI.
Notations: For matrix A, A * and A represent the conjugate operator and converged solution, respectively. Re{a} represents the real part of a complex value a. C M denotes the set of M × 1 complex vectors. E{·} denotes the expectation operation. For two matrices A and B, A B represents Hadamard product of A and B. A F , tr (A) and |A| denote the Frobenius norm, trace operation and determinant of A, respectively. ∇f x (x) denotes the gradient of the function f with respect to (w.r.t.) the vector x. CN (0, I) represents a random vector following the distribution of zero mean and unit variance matrix. arg{·} means the extraction of phase information. diag(·) denotes the diagonalization operation. (·) * , (·) T and (·) H denote the conjugate, transpose and Hermitian operators, respectively. arg(·) means the phase extraction operation.

A. System Model
Consider the IRS-aided multiuser MIMO downlink of a SWIPT system operating over the same frequency band both for data and energy transmission, as shown in Fig. 1. Let us assume that there are K I IRs and K E ERs, respectively. It is also assumed that the BS is equipped with N B ≥ 1 antennas, while each IR and ER is equipped with N I ≥ 1 and N E ≥ 1 antennas, respectively. Let us denote the sets of IRs and ERs as K I and K E , respectively. In general, low-power sensors require a certain amount of power (e.g., 0.1 mW) for their real-time operation. Due to the associated severe channel attenuation, the sensors should be deployed sufficiently close to the BS, which limits their practical implementation. To resolve this issue, we propose to employ an IRS, which has M reflective elements in the ERs' vicinity for extending the operational range of sensors, as shown in Fig. 1. Firstly, the IRS increases the energy harvested by the ERs, and additionally it also assists in enhancing the signal strength for distant IRs through careful phase shift optimization.
The number of data streams destined for each IR is assumed to be d, satisfying 1 ≤ d ≤ min{N B , N I }. The signal transmitted from the BS is given by where s k ∈ C d×1 is the (d × 1)-element data symbol vector designated for the kth IR satisfying E s k s H k = I d and E s i s H j = 0, for i = j, while F k ∈ C NB ×d is the linear TPC matrix used by the BS for the kth IR. Assuming nondispersive narrow-band transmission, the baseband equivalent channels spanning from the BS to the IRS, from the BS to the kth IR, from the BS to the lth ER, from the IRS to the kth IR, and finally from the IRS to the lth ER are modelled by the matrices respectively. Let us denote the diagonal reflection-coefficient matrix at the IRS by Φ = diag e jθ1 , · · · , e jθm , · · · , e jθM , 1 where θ m ∈ [0, 2π] is the phase shift of the m-th reflective element. Due to absorption and diffraction, the signal power that has been reflected multiple times is ignored. As a result, the signal received at the kth IR is given by where n I,k is the kth IR's noise vector satisfying CN 0, σ 2 I I NI . Similarly, the signal received at the lth ER is given by where n E,l is the lth ER's noise vector obeying the distribution of CN 0, σ 2 E I NE . We assume that all the CSIs are perfectly known at the BS, and the BS is responsible for calculating the phase shifts of the IRS, which are then fed back by them to the IRS controller through dedicated feedback channels. Given this idealized and simplified assumption, the results obtained represent a performance upper bound of how much performance gain can be achieved by an IRS. Let us define the equivalent channel spanning from the BS to the kth IR byH k Upon substituting x into (2), y I,k can be rewritten as Then, the achievable data rate (nat/s/Hz) of the kth IR is given by [33] where F denotes the collection of TPC matrices, while J k is the interference-plus-noise covariance matrix given by On the other hand, due to the broadcast nature of wireless channels, the ERs can extract energy from the electromagnetic wave. In general, the harvested power is nonlinear over the received radio frequency (RF) power due to the nonlinear RFto-DC conversion, which depends on the input RF power level. This nonlinear EH model has been characterized in [34], which is a complex function of the RF power. Based on this nonlinear EH model, various transmission designs have been proposed in [35] and [36]. However, there is still lack of a general model that can accurately characterize this nonlinear relationship by capturing all practical factors. Hence, for simplicity, we adopt the simple linear EH model as widely used in the existing literature [29], [37], [38]. By ignoring the noise power at the ERs, the total harvested power is proportional to the total received power. Let us define the equivalent channel spanning from the BS to the lth ER byḠ l Δ = G b,l + G r,l ΦZ. Then, the total power harvested by the lth ER is 1 j is the imaginary unit.
where 0 < η ≤ 1 is the energy harvesting efficiency. In this paper, we consider the constraint that the weighted sum of the power harvested by all ERs should be higher than a predefined value, which is where G = KE l=1 α l ηḠ H lḠ l , α l is the energy weighting factor of the lth ER, with a higher value of α l representing a higher priority for the lth ER than for others. Finally,Q is the minimum harvested power threshold.

B. Problem Formulation
Upon introducing the notations of φ m = e jθm , ∀m, we have Φ = diag {φ 1 , · · · , φ M }. Again, we aim for jointly optimizing the TPC matrices F and phase shift matrix Φ with the goal of maximizing the WSR of all IRs subject to the total power budget, to the unit modulus of the phase shifters and to the harvested power requirement. Then, this problem can be formulated as follows: where ω k is the weighting factor controlling the scheduling priority for each IR and P T is the power limit at the BS, while (8d) is the unit-norm constraint imposed on the phase shifters.
As the IRS is passive and both the ERs and IRs are energy constrained, we assume that this optimization problem is solved at the BS which posses the knowledge of the CSI of all related links and other related parameters such asQ. After computing the phase shift values for the IRS, they are sent to the IRS controller through dedicated control channels. Problem (8) is difficult to solve, since the TPC matrices and the phase shifts are coupled. If we remove the energy harvesting (EH) constraint, the problem reduces to the WSR maximization problem recently studied in [24]. However, the additional EH constraint makes the optimization more challenging to solve and the algorithms developed in [24] cannot be directly applied for two reasons. Firstly, the EH constraint is non-convex. Secondly, this problem may be infeasible due to the conflicting constraints (8b) and (8c). In the following, we first conceive a low-complexity algorithm to solve this problem by assuming that it is feasible. Then, we study the feasibility of this problem.

III. LOW-COMPLEXITY ALGORITHM DEVELOPMENT
In this section, we first transform Problem (8) into a more tractable one, which allows the decoupling of the TPC matrices and of the phase shifts. Then, the classic block coordinate descent (BCD) algorithm [33] is proposed for solving the transformed problem.

A. Reformulation of the Original Problem
To deal with the complex objective function, we reformulate Problem (8) by employing the well-known WMMSE method [39]. The appealing feature of this method is that it can transform the original complex problem into an equivalent form, which facilitates the application of the BCD method.
Specifically, the linear decoding matrix U is applied to estimate the signal vectorŝ k for each IR where U k ∈ C NI ×d is the decoding matrix of the kth IR. Then, the MSE matrix of the kth IR is given by where s and n denote the collections of data symbols and noise vectors of all IRs, respectively. By introducing a set of auxiliary matrices W = {W k 0, ∀k ∈ K I } and defining U = {U k , ∀k ∈ K I }, Problem (8) can be reformulated as follows [33], [39]: where h k (W, U, F, Φ) is given by Although Problem (11) has more optimization variables than Problem (8), the objective function (OF) in Problem (11) is much easier to handle, which allows the BCD algorithm to solve this problem by iteratively obtaining one set of variables while keeping the others fixed. Note that the decoding matrices U and the auxiliary matrices W only appear in the function h k (W, U, F, Φ). Hence, the optimal solution of U and W can be obtained while keeping the other matrices fixed. Specifically, given Φ, W, and F, setting the first-order derivative of h k (W, U, F, Φ) with respect to U k and W k to zero, we can obtain the optimal solution of U k and W k respectively as follows where E k is obtained by inserting U k into the kth IR's MSE matrix in (10), yielding In the following, we focus our attention on the optimization of TPC matrices F and phase shifts Φ, when U and W are given.

B. Optimizing the Precoding Matrices F
In this subsection, we aim to optimize the TPC matrices F with fixed W, U and Φ. By inserting E k in (10) into the OF of (11) and discarding the constant terms, the TPC matrices of our optimization problem can be transformed as follows where However, due to the non-convexity of the EH constraint, Problem (15) is still non-convex. To resolve this issue, we observe that it can be viewed as a difference of convex (d.c.) program, which can be efficiently solved by the successive convex approximation (SCA) method [40]. In particular, we can approximate it by its first-order Taylor expansion. By applying [41,Appendix B] and Jensen' inequality, we have where F (n) k , ∀k is the solution obtained from the previous iteration. Then, upon replacing the constraint (8c) by the following constraint: whereQ =Q + tr , we may consider the following optimization problem: Since the OF is convex w.r.t. F, and the constraints (8b) and (17) are convex, Problem (18) constitutes a convex optimization problem, which can be solved by standard convex solver packages, such as CVX [42]. However, the resultant computational complexity is high. In the following, we provide a low-complexity algorithm for obtaining a nearly optimal closed-form solution by resorting to the Lagrangian dual decomposition method [43]. Since Problem (18) is a convex problem and satisfies the slater's condition, 2 , the dual gap is 2 According to line 1 in Algorithm 2 the initial precoding matrix is initialized by the solution obtained from Section IV. Assume the original problem is feasible. Due to the randomness of channel matrices of G and H, the precoding matrix obtained in Section IV must be strictly larger than the minimum EH requirement, i.e., tr Then, based on [29], there must exist a strictly feasible solution, and thus the slater's condition holds. zero and the optimal solution can be obtained by solving its dual problem instead of its original one. We first introduce the Lagrange multiplier λ associated with the power constraint, and derive the partial Lagrangian function of Problem (18) as follows The dual function can be obtained by solving the following problem Then, the dual problem is given by Before solving the dual problem (21), we have to derive the expression of the dual function g (λ) by solving Problem (20) for a given λ. By introducing the dual variable μ ≥ 0 associated with the constraint (17), the Lagrangian function for Problem (20) is given by By setting the first-order derivative of L (F, μ) w.r.t. F * k to the zero matrix, we obtain the optimal solution of F k as follows: where (·) † denotes the matrix pseudoinverse. The value of μ should be chosen for ensuring that the complementary slackness condition for constraint (17) is satisfied: Hence, if the following condition holds the optimal solution of Problem (20) is given by F k (0), ∀k ∈ K I . Otherwise, the optimal μ is With the aid of the dual function, we may now commence the solution of the dual problem (21) to find the optimal λ. Given λ, we denote the optimal solution of Problem (20) by  F k (λ). The value of λ should be chosen for ensuring that the complementary slackness condition for the power constraint is satisfied: If the following condition holds: then the optimal solution is given by F k (0). Otherwise, we have to find λ for ensuring that the following equation holds: Unfortunately, due to the complex expression of μ in (26), we are unable to prove its monotonic nature by using the explicit expression of P (λ) as in [24]. In the following lemma, we prove that P (λ) is a monotonically decreasing function of λ, which enables the bisection search method to find λ. Lemma 1: The total power P (λ) is a monotonically decreasing function of λ.
Proof: Please refer to Appendix A. Based on Lemma 1, the bisection search method can be used for finding the solution of equation (29). In Algorithm 1, we provide the detailed steps of solving Problem (18) for the case of λ > 0. In each iteration of Algorithm 1, we have to calculate F k (μ) in (23), which involves the calculation of (A + λI) −1 at a complexity order of O (N 3  B ). If the total number of iterations is T , then the total complexity of calculat- , which may be excessive. Here, we provide a method for reducing the computational complexity. Specifically, as A is a non-negative definite matrix, it can be decomposed as A = QΛQ H by using the singular value decomposition (SVD), where QQ H = Q H Q = I NT and Λ is a diagonal matrix with non-negative diagonal elements. Then, we have (A + λI) −1 = Q (λI + Λ) −1 Q H . Hence, in each iteration, we only have to calculate the product of two matrices, which has much lower complexity than calculating the inverse of the matrix having the same dimension.
Based on the above discussions, in Algorithm 2 we provide the detailed steps of the SCA algorithm conceived for solving Problem (15).

Algorithm 1 Bisection Search Method to Solve Problem
Algorithm 2 SCA Algorithm to Solve Problem (15) 1: Initialize the accuracy ε, the precoding matrices F (0) from Section 2, the iteration index n = 0, the maximum number of iterations n max , calculate the OF value of Problem (15) as z(F (0) ); , ∀k} by solving Problem (18) using Algorithm 1; terminate. Otherwise, set n ← n + 1 and go to step 2.
Proof: The proof is similar to that of [44] and hence it is omitted for simplicity.
Next, we briefly analyze the complexity of Algorithm 2. We assume that N B ≥ N I ≥ d. In each iteration of Algorithm 2, the main complexity contribution is the calculation of {F

C. Optimizing the Phase Shift Matrix Φ
In this subsection, we focus our attention on optimizing the phase shift matrix Φ, while fixing the other parameters. Upon substituting E k in (10) into (12) and removing the terms that are independent of Φ, the phase shift optimization problem is formulated as: By substitutingH k = H b,k + H r,k ΦZ into (30a), we have and By using (31), we arrive at: where const 1 is a constant term that is independent of Φ.
Similarly, by defining T k where const 2 is a constant term that is independent of Φ.
α l ηG H r,l G r,l , and G br Δ = ZF KE l=1 α l ηG H b,l G r,l , the EH constraint in (8c) can be recast as follows: By inserting (33) and (34) Upon denoting the collections of diagonal elements of Moreover, the constraint (35) can be rewritten as where we have Q =Q−Tr G bF and Υ = G r C T . It can be verified that G r and C T are non-negative semidefinite matrices. Then, according to [45], the Hadamard product G r C T (or equivalently Υ) is also a semidefinite matrix. Thus, Problem (36) can be transformed as where we have Ξ = B C T . Again, B can be verified to be a non-negative semidefinite matrix, and thus Ξ is a non-negative semidefinite matrix. Due to the non-convex constraint (39), Problem (40) is difficult to solve. To deal with this constraint, we again employ the SCA method [40]. Specifically, since φ H Υφ is convex w.r.t. φ, its lower bound can be obtained as follows: where φ (n) is obtained in the previous iteration. Then, constraint (39) is replaced by the following constraint which is a linear constraint. Then, Problem (40) then becomes (42).
The key idea is to solve a challenging problem by introducing a series of more tractable subproblems. Upon denoting the objective function of Problem (43) by f (φ), in the (n + 1)th iteration we have to find the upper bound of the OF, denoted as g(φ|φ (n) ), which should satisfy the following three conditions: Then, we solve the approximate subproblem defined by a more tractable new OF g(φ|φ (n) ). To find g(φ|φ (n) ), we introduce the following lemma [46]. Lemma 2: For any given φ (n) , the following inequality holds for any feasible φ: where X = λ max I M and λ max is the maximum eigenvalue of Ξ. Then, the function g(φ|φ (n) ) can be constructed as follows: where y(φ|φ (n) ) is defined in (45). The new OF g(φ|φ (n) ) is more tractable than the original OF f (φ). The subproblem to be solved is given by (42).
Since φ H φ = M , we have φ H Xφ = M λ max , which is a constant. By removing the other constants, Problem (47) can be rewritten as follows: where q (n) = (λ max I M − Ξ) φ (n) − v * . Due to the additional constraint (42), the optimal solution of Problem (48) cannot be obtained as in [24]. Furthermore, due to the non-convex unit-modulus constraint (8d), Problem (48) is a non-convex optimization problem. As a result, the Lagrangian dual decomposition method developed for the convex problem (18) is not applicable here, since the dual gap is not zero.
In the following, we propose a price mechanism for solving Problem (48) that can obtain the globally optimal solution. Specifically, we consider the following problem by introducing a non-negative price p on the left hand side of constraint (42): For a given p, the globally optimal solution is given by φ(p) = e j arg(q (n) +p(g * +Υφ (n) )) .
Case II: Since p > 0, equation (51) holds only when J(p) =Q. To solve this equation, we first provide the following lemma.

Lemma 3: Function J(p) is a monotonically increasing function of p.
Proof: The proof is similar to Lemma 1 and thus omitted.
Based on Lemma 3, the bisection search method can be adopted for finding the solution of J(p) =Q. Based on the above discussions, we provide the algorithm to solve Problem (48) in Algorithm 3. Although Problem (48) is a non-convex problem, in the following theorem we prove that Algorithm 3 is capable of finding the globally optimal solution. Theorem 2: Algorithm 3 is capable of finding the globally optimal solution of Problem (48) and thus also of Problem (47). (48) 1: Calculate J(0). If J(0) ≤Q, terminate. Otherwise, go to step 2. 2: Initialize the accuracy ε, bounds p l and p u ; 3: Calculate p = (p l + p u )/2; 4: Update φ(p) in (50) and calculate J(p); 5: If J(p) ≥Q, set p u = p; Otherwise, set p l = p; 6: If |p l − p u | ≤ ε, terminate; Otherwise, go to step 3.

Algorithm 3 Bisection Search Method to Solve Problem
Proof: Please refer to Appendix B. Based on the above, we now provide the details of solving Problem (30) in Algorithm 4.

Algorithm 4 MM Combined with SCA Algorithm to Solve
Problem (30) 1: Initialize the accuracy ε, the phase shifts φ (0) , the iteration index to n = 0, the maximum number of iterations to n max , calculate the OF value of Problem (43) as f (φ (0) ); ) ≤ ε holds, terminate; Otherwise, set n ← n + 1 and go to step 2.
In the following theorem, we prove that the sequence of {φ (n) , n = 1, 2, · · · } generated by Algorithm 4 converges to the KKT-optimal point of Problem (30).

D. Overall Algorithm to Solve Problem (8)
Based on the above analysis, we provide the detailed steps of the BCD algorithm to solve Problem (8) in Algorithm 5, where R(F (n) , φ (n) ) denotes the OF value of Problem (8) in the nth iteration.
The following theorem shows the convergence and solution properties of Algorithm 5.
IV. FEASIBILITY CHECK FOR PROBLEM (8) Due to the conflicting EH and limited transmit power constraints, Problem (8) may be infeasible. Hence, we have to first check whether Problem (8) is feasible or not. To this end, we construct the following optimization problem: If the optimal OF value is larger thanQ, Problem (8) is feasible. Otherwise, it is infeasible. As the TPC matrices and phase shift matrix are coupled, the globally optimal solution is difficult to obtain. In the following, we can obtain a suboptimal solution by alternately optimizing the TPC matrices and phase shifts. For a given phase shift matrix, the TPC matrix optimization problem is given by Upon denoting the maximum eigenvalue and the corresponding eigenvector of G by χ and b respectively, the optimal solution can be readily obtained as −1) , ∀k = 1, · · · , K I , where KI k=1 p k = P T and p k ≥ 0, ∀k = 1, · · · , K I . Without loss of generality, we can set p i = P T /K I , ∀i ∈ K I . The OF value is given by χP T . In this case, the optimal TPC matrix represents the optimal energy beamforming, which is the same as that for the single-antenna IR case of [38].
For a given TPC matrix F, the phase shift optimization problem is formulated as: where Υ and g are defined in the above section. The OF is convex w.r.t. φ, and maximizing a convex function is a d.c program. Hence, it can be solved by using the SCA method by approximating φ H Υφ as its first-order Taylor expansion, details of which are omitted. Finally, alternately solve Problem (53) and (54) until the OF is larger thanQ.

V. SIMULATION RESULTS
In this section, we provide simulation results for demonstrating the benefits of applying IRS to SWIPT systems, as seen in Fig. 2, where there are four ERs and two IRs. The ERs and IRs are uniformly and randomly scattered in a circle centered at (x ER , 0) and (x IR , 0) with radius 1 m and 4 m, respectively. The IRS is located at (x IRS , 2). In the simulations, we assume that the IRS is just above the ERs and thus we set x ER = x IRS . The large-scale path loss is modeled in dB as where PL 0 is the path loss at the reference distance D 0 , D is the link length in meters, and α is the path loss exponent.
Here Due to the severe blockage and long distance, the channels from the BS and the IRS to the IRs are assumed to be Rayleigh fading. However, as the BS, the ERs and the IRS are close to each other, the small-scale channels are assumed to be Rician fading. In particular, the small-scale channels from the IRS to the ERs are denoted as: where β irser is the Rician factor,G LoS r,l is the deterministic line of sight (LoS), andG NLoS We first compare the maximum power harvested by various schemes in Fig. 3. Specifically, we solve the EH maximization problem (52) by using the feasibility check method in Section IV. Additionally, we also present the results without using IRS. Fig. 3 shows the maximum EH power versus the ER circle center location x EH . As expected, the EH power gleaned by all schemes decreases, when the ERs are far away from the BS. As expected, more power can be harvested with the aid of IRS than that without IRS, especially when the number of phase shifters M is large. This is mainly due to the fact that an additional strong link is reflected by the IRS, which can be harvested by the ERs. This figure also shows that the IRS is effective in expanding the operational range of ERs. For example, when the harvested power limit isQ = 2 × 10 −4 W, the maximum operational range of the system without IRS is only 5.5 m, while the system having M = 40 phase shifters can operate for distances up to 9 m.
In Fig. 4, we study the convergence behaviour of the BCD algorithm for different numbers of phase shifters M . It is observed from Fig. 4 that the WSR achieved for various M values increases monotonically with the number of iterations, which verifies Theorem 4. Additionally, the BCD algorithm converges rapidly and in general a few iterations are sufficient for the BCD algorithm to achieve a large portion of the converged WSR. This reflects the low complexity of the BCD algorithm, which is appealing for practical applications.
In the following, we compare our proposed BCD algorithms to a pair of benchmark schemes: 1)'No-IRS': In this scheme, there is no IRS to assist the transmission as in conventional systems; 2) 'Fixed Phase': In this method, the phase shifts are fixed at the solutions obtained by solving the harvested power maximization problem (52), while they are not optimized, when using the BCD algorithm by removing Step 3 of the BCD algorithm. When any of the methods fails to obtain a feasible solution, its achievable WSR is set to zero.
In Fig. 5, we first study the impact of the number of phase shifters M on the performance of various algorithms. As expected, the WSR achieved by all the algorithms -except for the No-IRS method -increases with M , since a higher degree of freedom can be exploited for optimizing the system performance. By carefully optimizing the phase shifts for maximizing the WSR, the BCD algorithm significantly outperforms the fixed-phase scheme. Additionally, the performance gain increases with M , which emphasizes the importance of optimizing the phase shifts. By employing the IRS in our SWIPT system, the WSR obtained by the BCD algorithm becomes drastically higher than that of No-IRS. For example, when M = 60, the WSR performance gain is up to 10 bit/s/Hz. These results demonstrate that introducing the IRS into our SWIPT system is a promising technique of enhancing the system performance.  In Fig. 6, the impact of harvested power requirementQ is investigated. It is seen from this figure that the WSR achieved by all the algorithms decreases upon increasingQ, because the probability of infeasibility increases, which in turn reduces the average WSR value. We also find that the WSR obtained by the No-IRS scheme decreases more rapidly than that of the other two IRS-aided transmission schemes. The WSR of the No-IRS is approaching zero whenQ = 4×10 −4 W, while those relying on IRSs achieve a WSR gain in excess of 20 bit/s/Hz. It is observed again that the BCD algorithm performs better than the fixed-phase scheme, but the gap narrows with the increase ofQ. This can be explained as follows. With the increase ofQ, both the TPC matrices and the phase shifts should be designed for maximizing the power harvested at the ERs, and thus the final solutions of the fixed-phase and BCD method will become the same.
The above results are obtained for α BSIRS = 2.2, α IRSER = 2.2, α IRSIR = 2.4 based on the assumption that the IRS relies on an obstacle-free scenario. In practice, this ideal scenario is seldom encountered. Hence, it is imperative to investigate the impact of α IRS Δ = α BSIRS = α IRSER = α IRSIR on the system performance, which is shown in Fig. 7. Observe from this figure that the WSR achieved by the algorithms  using IRS decreases drastically with α IRS . When α IRS = 3, the WSR-performance gain of our algorithm over the No-IRS scenario is only 7 bit/s/Hz, because upon increasing α IRS , the signal power reflected from the IRS becomes weaker. Hence, the benefits of the IRS can be eroded. This provides an important engineering design insight: the location of IRS should be carefully considered for finding an obstacle-free scenario associated with a low α IRS .
In Fig. 8, we study the impact of ER locations on the system performance. As expected, the WSR achieved by all the schemes decreases with x IRS , since the ERs become more distant from the BS and the signals gleaned from both the BS and IRS become weaker. The WSR achieved by the No-IRS approaches zero when x IRS = 8 m, hence this method cannot reach the energy transmission target of the ERs. The proposed algorithm is again observed to significantly outperform the other two algorithms, especially when the ERs are close to the BS.
Finally, the impact of IR locations is investigated in Fig. 9. It is observed that the WSR achieved by all the algorithms decreases with x IR since the IRs become farther away from the BS when increasing x IR . The proposed algorithm is shown to achieve nearly the WSR gain of 10 bit/s/Hz over the No-IRS when x IR = 100 m, and the WSR gain slightly increases with x IR . This means that the IRS is more advantageous when the IRs are far away from the BS, and the IRS can provide one additional favorable link.

VI. CONCLUSION
In this paper, we have invoked an IRS in a SWIPT MIMO system for enhancing the performance of both the ERs and IRs. By carefully adjusting the phase shifts at the IRS, the signal reflected by the IRS can be added constructively at both the ERs and IRs. We considered the WSR maximization problem of IRs, while guaranteeing the energy harvesting requirements of the ERs and the associated non-convex unit-modulus constraints. We conceived a BCD algorithm for alternatively optimizing the TPC matrices at the BS and the phase shift matrix at the IRSs. For each subproblem, a low-complexity iterative algorithm was proposed, which guarantees to be at worst locally optimal. Our simulation results demonstrated that the IRS enhances the performance of the SWIPT system and that the proposed algorithm converges rapidly, hence it is eminently suitable for practical implementations.
This paper assumes perfect CSI at the BS, which is challenging to obtain. For the future work, we will consider the robust transmission design for the IRS-aided SWIPT system, where the CSI is assumed to be imperfectly known. In addition, how to design the discrete phase shifts will be left for future work.
We consider a pair of variables λ and λ , where λ > λ . Let F(λ) and F(λ ) be the optimal solutions of Problem (20) with λ and λ , respectively. Since F(λ) is the optimal solution of Problem (20)  By adding these two inequalities and simplifying them, Since λ > λ , we have P (λ) ≤ P (λ ), which completes the proof.

APPENDIX B PROOF OF THEOREM 2
Denote the globally optimal solution of Problem (48) by φ . According to [43], for a non-convex optimization problem, all its locally optimal solutions (including the globally optimal solution) should satisfy the Karush-Kuhn-Tucker (KKT) optimality conditions, one of which is the complementary slackness condition for constraint (42): where λ is the corresponding optimal dual variable. We consider two cases: 1) λ = 0; 2) λ > 0.
The first case means that constraint (42) is not tight in the optimum. Then, the optimal solution can be obtained as φ = e j arg(q (n) ) , which is equal to φ(0). Hence, Algorithm 3 achieves the optimal solution of Problem (48).
For the second case, the following equality should hold: We prove the second case by using the method of contradiction. Denote the optimal p obtained by Algorithm 3 as p , and the corresponding φ as φ(p ). Then, we have Let us assume that φ(p ) is not the globally optimal solution of Problem (48). Then, we have Since φ(p ) is the globally optimal solution of Problem (49) when p = p , we have Let us define the following functions: It can be verified that T (φ (n) ) =T (φ (n) |φ (n) ).
In the following, we prove that the converged solution satisfies the KKT conditions of Problem (8). Let us denote the converged solution as {W , U , F , Φ }.
According to Theorem 1, F is the KKT-optimum point of Problem (15). Upon denoting the OF of Problem (15) as z(F, Φ ), the Lagrange function of Problem (15) is given by where λ and μ are the corresponding dual variables. Then, there must exist a λ and μ for ensuring that the following conditions are satisfied 3 : Furthermore, it can be readily checked that To expound a little further, we have the following chain of inequalities: where (D.8) follows from the chain rule, and the final equality follows from applying the Woodbury matrix identity to (14). 3 For simplicity, the prime constraints are omitted.
Furthermore, it can be readily verified that By using similar derivations as in (D.7)-(D.11), we can prove that (D.15) Hence, we have By substituting (D. 16) into (C.8), we arrive at: