Intelligent Reflecting Surface Aided MIMO Broadcasting for Simultaneous Wireless Information and Power Transfer
Cunhua Pan, Hong Ren, Kezhi Wang, Maged Elkashlan, Arumugam Nallanathan, Jiangzhou Wang, Lajos Hanzo
11 Intelligent Reflecting Surface Aided MIMOBroadcasting for Simultaneous WirelessInformation and Power Transfer
Cunhua Pan, Hong Ren, Kezhi Wang, Maged Elkashlan, Arumugam Nallanathan,
Fellow, IEEE , Jiangzhou Wang,
Fellow , IEEE and Lajos Hanzo,
Fellow, IEEE
Abstract
An intelligent reflecting surface (IRS) is invoked for enhancing the energy harvesting performance of asimultaneous wireless information and power transfer (SWIPT) aided system. Specifically, an IRS-assistedSWIPT system is considered, where a multi-antenna aided base station (BS) communicates with severalmulti-antenna assisted information receivers (IRs), while guaranteeing the energy harvesting requirementof 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 tacklethis challenging optimization problem, we first adopt the classic block coordinate descent (BCD) algorithmfor decoupling the original optimization problem into several subproblems and alternatively optimize theTPC matrices and the phase shift matrix. For each subproblem, we provide a low-complexity iterativealgorithm, 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 alsoconceive a feasibility checking method to study its feasibility. Our extensive simulation results confirm thatemploying IRSs in SWIPT beneficially enhances the system performance and the proposed BCD algorithmconverges rapidly, which is appealing for practical applications.
Index Terms
Intelligent Reflecting Surface (IRS), Large Intelligent Surface (LIS), SWIPT, Energy Harvesting,MIMO.
C. Pan, H. Ren, M. Elkashlan and A. Nallanathan are with the School of Electronic Engineering and Computer Science at QueenMary University of London, London E1 4NS, U.K. (e-mail: { c.pan, h.ren, maged.elkashlan, a.nallanathan } @qmul.ac.uk). K. Wangis with Department of Computer and Information Sciences, Northumbria University, UK. (e-mail: [email protected]).J. Wang is with the School of Engineering and Digital Arts, University of Kent, Canterbury, Kent, CT2 7NZ, U.K. (e-mail:[email protected]). L. Hanzo is with the School of Electronics and Computer Science, University of Southampton, Southampton,SO17 1BJ, U.K. (e-mail: [email protected]). a r X i v : . [ ee ss . SP ] F e b I. I
NTRODUCTION
Recently, intelligent reflecting surface (IRS)-assisted wireless communication has received con-siderable research attention, since it is capable of supporting cost-effective and energy-efficient highdata rate communication for next-generation communication systems [1]–[3]. In simple tangibleterms, an IRS is composed of a vast number of low-cost and passive reflective components, eachof which is capable of imposing a phase change on the signals incident upon them. Thanks tothe recent advances in meta-materials [4], it has become feasible to reconfigure the phase shiftsin real time. As a result, the phase shifts of all reflective components can be collaborativelyadjusted for ensuring that the signals reflected from the IRS can be added constructively ordestructively at the receiver in order to beneficially steer the signal component arriving fromthe base station (BS) for enhancing the desired signal power or alternatively for suppressing theundesired signals, such as interference. In contrast to conventional physical layer techniques thatare designed for accommodating the hostile time-varying wireless channels [5], [6], IRSs constitutea new paradigm capable of ‘reprogramming’ the wireless propagation environment into a morefavorable transmission medium. Since the reflective components are passive, they impose a muchlower power consumption than conventional relay-aided communication systems relying on activetransmission devices. Additionally, no thermal noise is imposed by the IRS, since it directly reflectsthe incident signals without decoding or amplifying them, which is in contrast to conventionalrelays. 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, asIRS is a complementary device, it can be readily integrated into current wireless networks withoutmodifying the physical layer standardization, making it transparent to the users. To fully exploitthe benefits of IRS, the active beamforming at the BS and the passive beamforming at the IRSshould be jointly designed. However, the optimization variables are coupled and the joint designleads to a complex optimization problem that is difficult to solve.Some innovative efforts have been devoted to the transceiver design when integrating IRSinto various wireless communication systems, including the single-user scenarios of [7]–[11], thedownlink 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] andthe multicell multiuser multiple-input multiple-output (MIMO) case in [24]. Concretely, Wu et al. proposed joint active and passive beamforming for a single-user scenario in [7], where semidefinite relaxation (SDR) was proposed for optimizing the phase shift matrix. However, its complexity ishigh since the number of optimization variables increases quadratically with the number of phaseshifts. Additionally, the Gaussian random approximation employed leads to certain performanceloss. To resolve this issue, Yu et al. [8] proposed a pair of efficient algorithms termed as fixed pointiteration and manifold optimization techniques, which can guarantee locally optimal solutions. Asa further advance, the authors of [9] considered realistic frequency-selective channels. The phaseshift 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 powerminimization 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 issueswere 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 bebeneficial in reducing the latency of MEC networks. In addition, the IRS can help enhance the WSRperformance for the multigroup multicast network in [23]. Most recently, we considered an IRS-assisted multicell MIMO communications scenario [24], where we demonstrated that deployingan IRS at the cell edge is also capable of mitigating the adjacent-cell interference. Channel stateinformation (CSI) is challenging to obtain in IRS-assisted communication system due to its passivefeature. There are some initial efforts to handle this issue such as channel estimation and/or robusttransmission design [25]–[28]. Specifically, Huang et al. [25] proposed a deep learning methodfor efficient online configuration of the phase shifts, where the phase values can be immediatelyobtained by inputting the user location into the trained deep neural network. A two-stage channelestimation method based on a sparse matrix factorization and a matrix completion was proposedin [26]. A pair of algorithms based on compressed sensing and deep learning were conceived byTaha et al. [27] for tackling the challenging channel estimation issues of IRS-assisted systems.Most recently, we first studied the robust beamforming design for IRS-assisted communicationsystems in [28], where the imperfect channel from an IRS to users was considered and the channelestimation error was assumed to be within a bounded elliptical region.On the other hand, information transmission enabled simultaneous wireless information and BS IRS Z ERs IRsReflected IFDirect IF
Reflected EFDirect EF , r l G , r k H , b l G , b k H EF (Energy flow) IF (Information flow)
Fig. 1. An IRS-assisted SWIPT system. 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 wirelesssignals 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 theinformation rate attained and the amount of harvested energy for a single-user MIMO system. Inpractice, a typical ER such as a humidity sensor requires much higher energy for its operation thanthat 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 IRSin the vicinity of ERs to provide additional transmission links to support the ERs for enhancing theirharvested power as shown in Fig. 1, since there is a paucity of IRS-assisted SWIPT contributions inthe literature [30]. Explicitly in [30], the weighted sum power maximization problem was studiedby Wu and Zhang, who proved that no dedicated energy-carrying signals were required for anIRS-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 degradationdue to the associated rank-one extraction. However, this method is not applicable when each useris 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 isinstalled in the vicinity of ERs for compensating the associated power loss, while maximizing theWSR 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 onthe phase shifts, while simultaneously satisfying the energy harvesting requirement of theERs. To the best of our knowledge, this is the first treatise considering the WSR maximizationproblem of IRS-assisted SWIPT MIMO systems, which is much more challenging thanthe weighted sum power minimization problem of [30] since the latter can be readilytransformed into a convex optimization problem. In contrast to the multicell system of [24], anadditional energy harvesting constraint is also imposed in our current study, which furthercomplicates the analysis. Specifically, this constraint is non-convex and the optimizationproblem 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 arecomplex. To deal with this issue, we first reformulate the original problem into an equivalentform by exploiting the equivalence between the data rate and the weighted minimum mean-square error (WMMSE). Then, an alternating optimization algorithm based on the popularblock coordinate descent (BCD) algorithm is proposed for alternatively updating the activeTPC matrices of the BS and the phase shift matrix of the IRS, which is rigorously provedto 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 onthe successive convex approximation (SCA) method and on the Lagrangian dual decomposi-tion method to derive a nearly closed-form solution for the TPC matrices. A low-complexitybisection search method is proposed for finding the optimal dual variables. The solutionsgenerated by our iterative algorithm are guaranteed to converge to the KKT point of theTPC 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 energyharvesting constraint by invoking some further matrix manipulations. Then, a novel iterativealgorithm 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 solutiongenerated by the iterative algorithm is guaranteed to converge to the KKT point of the phaseshift optimization problem.4) The associated feasibility issue is also studied by formulating an alternative optimizationproblem and an iterative algorithm is proposed for solving this problem.
5) Extensive simulation results are provided for verifying the performance advantages of em-ploying IRS in SWIPT in order to enhance the energy harvesting performance. It is shownthat the operating range of the ERs can be dramatically expanded by placing IRSs in theERs’ vicinity. Furthermore, the BCD algorithm converges rapidly, and it is eminently suitablefor practical applications. Our simulation results also show that as expected, the path lossexponent substantially affects the system’s performance and thus the location of the IRSshould be carefully chosen.The remainder of this paper is organized as follows. In Section II, we introduce the IRS-assistedSWIPT system model and our problem formulation. The detailed algorithms used for solving theoptimization problem are presented in Section III. The feasibility issues of the original problemare 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 (cid:63) 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 × complex vectors. E {·} denotes the expectation operation. For two matrices A and B , A (cid:12) B represents Hadamard product of A and B . (cid:107) A (cid:107) F , tr ( A ) and | A | denote the Frobenius norm, traceoperation and determinant of A , respectively. ∇ f x ( x ) denotes the gradient of the function f withrespect to (w.r.t.) the vector x . CN ( , I ) represents a random vector following the distribution ofzero mean and unit variance matrix. arg {·} means the extraction of phase information. diag( · ) denotes the diagonalization operation. ( · ) ∗ , ( · ) T and ( · ) H denote the conjugate, transpose andHermitian operators, respectively. arg( · ) means the phase extraction operation.II. S YSTEM M ODEL AND P ROBLEM F ORMULATION
A. System Model
Consider the IRS-aided multiuser MIMO downlink of a SWIPT system operating over the samefrequency band both for data and energy transmission, as shown in Fig. 1. Let us assume that thereare K I IRs and K E ERs, respectively. It is also assumed that the BS is equipped with N B ≥ antennas, while each IR and ER is equipped with N I ≥ and N E ≥ antennas, respectively. Let usdenote the sets of IRs and ERs as K I and K E , respectively. In general, low-power sensors requirea certain amount of power (e.g., 0.1 mW) for their real-time operation. Due to the associatedsevere channel attenuation, the sensors should be deployed sufficiently close to the BS, whichlimits 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, asshown in Fig. 1. Firstly, the IRS increases the energy harvested by the ERs, and additionally it alsoassists 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 ≤ d ≤ min { N B , N I } . The signal transmitted from the BS is given by x = K I (cid:88) k =1 F k s k , (1)where s k ∈ C d × is the ( d × -element data symbol vector designated for the k th IR satisfying E (cid:2) s k s Hk (cid:3) = I d and E (cid:2) s i s Hj (cid:3) = , for i (cid:54) = j , while F k ∈ C N B × d is the linear TPC matrix used bythe BS for the k th IR. Assuming non-dispersive narrow-band transmission, the baseband equivalentchannels spanning from the BS to the IRS, from the BS to the k th IR, from the BS to the l th ER,from the IRS to the k th IR, and finally from the IRS to the l th ER are modelled by the matrices Z ∈ C M × N B , H b,k ∈ C N I × N B , G b,l ∈ C N E × N B , H r,k ∈ C N I × M , and G r,l ∈ C N E × M , respectively. Let usdenote the diagonal reflection-coefficient matrix at the IRS by Φ = diag (cid:8) e jθ , · · · , e jθ m , · · · , e jθ M (cid:9) , where θ m ∈ [0 , π ] is the phase shift of the m -th reflective element. Due to absorption anddiffraction, the signal power that has been reflected multiple times is ignored. As a result, thesignal received at the k th IR is given by y I,k = ( H b,k + H r,k ΦZ ) x + n I,k , (2)where n I,k is the k th IR’s noise vector satisfying CN ( , σ I I N I ) . Similarly, the signal received atthe l th ER is given by y E,l = ( G b,l + G r,l ΦZ ) x + n E,l , (3)where n E,l is the l th ER’s noise vector obeying the distribution of CN ( , σ E I N E ) .We assume that all the CSIs are perfectly known at the BS, and the BS is responsible forcalculating the phase shifts of the IRS, which are then fed back by them to the IRS controllerthrough dedicated feedback channels. Given this idealized and simplified assumption, the resultsobtained represent a performance upper bound of how much performance gain can be achievedby an IRS. Let us define the equivalent channel spanning from the BS to the k th IR by ¯H k ∆ = H b,k + H r,k ΦZ . Upon substituting x into (2), y I,k can be rewritten as y I,k = ¯H k F k s k + K I (cid:88) i =1 ,i (cid:54) = k ¯H k F i s i + n I,k . (4) j is the imaginary unit. Then, the achievable data rate (nat/s/Hz) of the k th IR is given by [33] R k ( F , Φ ) = log (cid:12)(cid:12) I + ¯H k F k F H k ¯H H k J − k (cid:12)(cid:12) , (5)where F denotes the collection of TPC matrices, while J k is the interference-plus-noise covariancematrix given by J k = (cid:80) K I m =1 ,m (cid:54) = k ¯H k F m F H m ¯H H k + σ I I .On the other hand, due to the broadcast nature of wireless channels, the ERs can extract energyfrom the electromagnetic wave. In general, the harvested power is nonlinear over the received radiofrequency (RF) power due to the nonlinear RF-to-DC conversion, which depends on the input RFpower level. This nonlinear EH model has been characterized in [34], which is a complex functionof the RF power. Based on this nonlinear EH model, various transmission designs have beenproposed in [35] and [36]. However, there is still lack of a general model that can accuratelycharacterize 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 totalreceived power. Let us define the equivalent channel spanning from the BS to the l th ER by ¯G l ∆ = G b,l + G r,l ΦZ . Then, the total power harvested by the l th ER is Q i = η tr (cid:32) K I (cid:88) k =1 ¯G l F k F H k ¯G H l (cid:33) , (6)where < η ≤ is the energy harvesting efficiency. In this paper, we consider the constraint thatthe weighted sum of the power harvested by all ERs should be higher than a predefined value,which is Q = K E (cid:88) l =1 α l Q l = tr (cid:32) K I (cid:88) k =1 F H k GF k (cid:33) ≥ ¯ Q, (7)where G = (cid:80) K E l =1 α l η ¯G H l ¯G l , α l is the energy weighting factor of the l th ER, with a higher valueof α l representing a higher priority for the l th ER than for others. Finally, ¯ Q is the minimumharvested power threshold. B. Problem Formulation
Upon introducing the notations of φ m = e jθ m , ∀ m , we have Φ = diag { φ , · · · , φ M } . Again,we aim for jointly optimizing the TPC matrices F and phase shift matrix Φ with the goal ofmaximizing the WSR of all IRs subject to the total power budget, to the unit modulus of the phaseshifters and to the harvested power requirement. Then, this problem can be formulated as follows: max F , Φ K I (cid:88) k =1 ω k R k ( F , Φ ) (8a)s.t. K I (cid:88) k =1 (cid:107) F k (cid:107) F ≤ P T , (8b) tr (cid:32) K I (cid:88) k =1 F H k GF k (cid:33) ≥ ¯ Q, (8c) | φ m | = 1 , m = 1 , · · · , M, (8d)where ω k is the weighting factor controlling the scheduling priority for each IR and P T is thepower 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 thisoptimization problem is solved at the BS which posses the knowledge of the CSI of all relatedlinks and other related parameters such as ¯ Q . After computing the phase shift values for the IRS,they are sent to the IRS controller through dedicated control channels. Problem (8) is difficult tosolve, 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 thealgorithms developed in [24] cannot be directly applied for two reasons. Firstly, the EH constraintis 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 byassuming that it is feasible. Then, we study the feasibility of this problem.III. L OW -C OMPLEXITY A LGORITHM D EVELOPMENT
In this section, we first transform Problem (8) into a more tractable one, which allows thedecoupling 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 thewell-known WMMSE method [39]. The appealing feature of this method is that it can transformthe original complex problem into an equivalent form, which facilitates the application of the BCDmethod. Specifically, the linear decoding matrix U is applied to estimate the signal vector ˆs k for eachIR ˆs k = U H k y I,k , ∀ k, (9)where U k ∈ C N I × d is the decoding matrix of the k th IR. Then, the MSE matrix of the k th IR isgiven by E k = E s , n (cid:104) ( ˆs k − s k ) ( ˆs k − s k ) H (cid:105) (10) = (cid:0) U H k ¯H k F k − I (cid:1)(cid:0) U H k ¯H k F k − I (cid:1) H + K I (cid:88) m =1 ,m (cid:54) = k U H k ¯H k F m F H m ¯H H k U k + σ U H k U k , ∀ k ∈ K I , (11)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 (cid:23) , ∀ k ∈ K I } and defining U = { U k , ∀ k ∈ K I } , Problem (8) can be reformulated as follows [33], [39]: max W , U , F , Φ K I (cid:88) k =1 ω k h k ( W , U , F , Φ ) (12a)s.t. (8 b ) , (8 c ) , (8 d ) , (12b)where h k ( W , U , F , Φ ) is given by h k ( W , U , F , Φ ) = log | W k | − Tr ( W k E k ) + d. (13)Although Problem (12) has more optimization variables than Problem (8), the objective function(OF) in Problem (12) is much easier to handle, which allows the BCD algorithm to solve thisproblem by iteratively obtaining one set of variables while keeping the others fixed. Note that thedecoding 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 respectto U k and W k to zero, we can obtain the optimal solution of U k and W k respectively as follows U (cid:63)k = (cid:0) J k + ¯H k F k F H k ¯H H k (cid:1) − ¯H k F k , W (cid:63)k = E (cid:63) − k , (14)where E (cid:63)k is obtained by inserting U (cid:63)k into the k th IR’s MSE matrix in (11), yielding E (cid:63)k = I d − F H k ¯H H k (cid:32) K I (cid:88) m =1 ¯H k F m F H m ¯H H k + σ I I (cid:33) − ¯H k F k . (15)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 Φ . Byinserting E k in (11) into the OF of (12) and discarding the constant terms, the TPC matricesof our optimization problem can be transformed as follows min F K I (cid:88) k =1 tr (cid:0) F H k AF k (cid:1) − K I (cid:88) k =1 ω k Tr (cid:0) W k U H k ¯H k F k (cid:1) − K I (cid:88) k =1 ω k tr (cid:0) W k F H k ¯H H k U k (cid:1) (16a)s.t. (8 b ) , (8 c ) , (16b)where A = (cid:80) K I m =1 ω m ¯H H m U m W m U H m ¯H m .However, due to the non-convexity of the EH constraint, Problem (16) is still non-convex. Toresolve 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]. Inparticular, we can approximate it by its first-order Taylor expansion. By applying [41, AppendixB] and Jensen’ inequality, we have tr (cid:32) K I (cid:88) k =1 F H k GF k (cid:33) ≥ − tr (cid:32) K I (cid:88) k =1 F ( n )H k GF ( n ) k (cid:33) + 2Re (cid:34) tr (cid:32) K I (cid:88) k =1 F ( n )H k GF k (cid:33)(cid:35) , (17)where (cid:110) F ( n ) k , ∀ k (cid:111) is the solution obtained from the previous iteration. Then, upon replacing theconstraint (8c) by the following constraint: (cid:34) tr (cid:32) K I (cid:88) k =1 F ( n )H k GF k (cid:33)(cid:35) ≥ ˜ Q, (18)where ˜ Q = ¯ Q + tr (cid:16)(cid:80) K I k =1 F ( n )H k GF ( n ) k (cid:17) , we may consider the following optimization problem: min F K I (cid:88) k =1 tr (cid:0) F H k AF k (cid:1) − K I (cid:88) k =1 ω k tr (cid:0) W k U H k ¯H k F k (cid:1) − K I (cid:88) k =1 ω k tr (cid:0) W k F H k ¯H H k U k (cid:1) (19a)s.t. (8 b ) , (18) . (19b)Since the OF is convex w.r.t. F , and the constraints (8b) and (18) are convex, Problem (19) con-stitutes 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 byresorting to the Lagrangian dual decomposition method [43]. Since Problem (19) is a convex problem and satisfies the slater’s condition , the dual gap is zero and the optimal solution can beobtained by solving its dual problem instead of its original one. We first introduce the Lagrangemultiplier λ associated with the power constraint, and derive the partial Lagrangian function ofProblem (19) as follows L ( F , λ ) = K I (cid:80) k =1 tr (cid:0) F H k AF k (cid:1) − K I (cid:80) k =1 ω k tr (cid:0) W k U H k ¯H k F k (cid:1) − K I (cid:80) k =1 ω k tr (cid:0) W k F H k ¯H H k U k (cid:1) + λ K I (cid:80) k =1 tr (cid:0) F H k F k (cid:1) − λP T . (20)The dual function can be obtained by solving the following problem g ( λ ) ∆ = min F L ( F , λ ) s . t . (18) . (21)Then, the dual problem is given by max λ g ( λ ) (22a)s.t. λ ≥ . (22b)Before solving the dual problem (22), we have to derive the expression of the dual function g ( λ ) by solving Problem (21) for a given λ . By introducing the dual variable µ ≥ associatedwith the constraint (18), the Lagrangian function for Problem (21) is given by L ( F , µ ) = K I (cid:80) k =1 tr (cid:0) F H k ( A + λ I ) F k (cid:1) − K I (cid:80) k =1 ω k tr (cid:0) W k U H k ¯H k F k (cid:1) − K I (cid:80) k =1 ω k tr (cid:0) W k F H k ¯H H k U k (cid:1) + µ ˜ Q − µ Re (cid:20) tr (cid:18) K I (cid:80) k =1 F ( n )H k GF k (cid:19)(cid:21) − λP T . (23)By setting the first-order derivative of L ( F , µ ) w.r.t. F ∗ k to the zero matrix, we obtain the optimalsolution of F k as follows: F (cid:63)k ( µ ) = ( A + λ I ) † (cid:16) ω k ¯H H k U k W k + µ GF ( n ) k (cid:17) , (24)where ( · ) † denotes the matrix pseudoinverse. The value of µ should be chosen for ensuring thatthe complementary slackness condition for constraint (18) is satisfied: µ (cid:32) (cid:34) tr (cid:32) K I (cid:88) k =1 F ( n )H k GF (cid:63)k ( µ ) (cid:33)(cid:35) − ˜ Q (cid:33) = 0 . (25) According to line 1 in Algorithm 2, the initial precoding matrix is initialized by the solution obtained from Section IV. Assumethe original problem is feasible. Due to the randomness of channel matrices of G and H , the precoding matrix obtained in SectionIV must be strictly larger than the minimum EH requirement, i.e., tr (cid:16)(cid:80) K I k =1 F (0)H k GF (0) k (cid:17) > ¯ Q . Then, based on [29], there mustexist a strictly feasible solution, and thus the slater’s condition holds. Hence, if the following condition holds (cid:34) tr (cid:32) K I (cid:88) k =1 F ( n )H k GF (cid:63)k (0) (cid:33)(cid:35) ≥ ˜ Q, (26)the optimal solution of Problem (21) is given by F (cid:63)k (0) , ∀ k ∈ K I . Otherwise, the optimal µ is µ = ˜ Q − (cid:20) tr (cid:18) K I (cid:80) k =1 ω k F ( n )H k G ( A + λ I ) − ¯H H k U k W k (cid:19)(cid:21) (cid:18) K I (cid:80) k =1 F ( n )H k G ( A + λ I ) − GF ( n ) k (cid:19) . (27)With the aid of the dual function, we may now commence the solution of the dual problem(22) to find the optimal λ . Given λ , we denote the optimal solution of Problem (21) by F k ( λ ) .The value of λ should be chosen for ensuring that the complementary slackness condition for thepower constraint is satisfied: λ (cid:32) tr (cid:32) K I (cid:88) k =1 F H k ( λ ) F k ( λ ) (cid:33) − P T (cid:33) = 0 . (28)If the following condition holds: tr (cid:32) K I (cid:88) k =1 F H k (0) F k (0) (cid:33) ≤ P T , (29)then the optimal solution is given by F k (0) . Otherwise, we have to find λ for ensuring that thefollowing equation holds: P ( λ ) ∆ = tr (cid:32) K I (cid:88) k =1 F H k ( λ ) F k ( λ ) (cid:33) = P T . (30)Unfortunately, due to the complex expression of µ in (27), we are unable to prove its monotonicnature 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 tofind λ . Lemma 1:
The total power P ( λ ) is a monotonically decreasing function of λ . Proof:
Please refer to Appendix A. (cid:3)
Based on Lemma 1, the bisection search method can be used for finding the solution of equation(30). In Algorithm 1, we provide the detailed steps of solving Problem (19) for the case of λ > .In each iteration of Algorithm 1, we have to calculate F (cid:63)k ( µ ) in (24), which involves the calculationof ( A + λ I ) − at a complexity order of O ( N B ) . If the total number of iterations is T , then thetotal complexity of calculating ( A + λ I ) − is O ( T N B ) , which may be excessive. Here, we provide a method for reducing the computational complexity. Specifically, as A is a non-negative definitematrix, it can be decomposed as A = QΛQ H by using the singular value decomposition (SVD),where QQ H = Q H Q = I N T and Λ is a diagonal matrix with non-negative diagonal elements. Then,we have ( A + λ I ) − = Q ( λ I + Λ ) − Q H . Hence, in each iteration, we only have to calculate theproduct of two matrices, which has much lower complexity than calculating the inverse of thematrix having the same dimension. Algorithm 1
Bisection Search Method to Solve Problem (19) Initialize the accuracy ε , the bounds λ l and λ u ; Calculate λ = ( λ l + λ u )/2 ; If condition (26) is satisfied, µ is equal to zero. Otherwise, update µ in (27); Calculate { F k ( λ ) , ∀ k } according to (24); If P ( λ ) ≥ P T , set λ l = λ . Otherwise, set λ u = λ ; If | λ l − λ u | ≤ ε , terminate. Otherwise, go to step 2.Based on the above discussions, in Algorithm 2 we provide the detailed steps of the SCAalgorithm conceived for solving Problem (16). Algorithm 2
SCA Algorithm to Solve Problem (16) 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 (16) as z ( F (0) ) ; Calculate ˜ Q ( n ) = ¯ Q + tr (cid:16)(cid:80) K I k =1 F ( n )H k GF ( n ) k (cid:17) ; With ˜ Q ( n ) , calculate { F ( n +1) k , ∀ k } by solving Problem (19) using Algorithm 1; If n ≥ n max or (cid:12)(cid:12) z ( F ( n +1) ) − z ( F ( n ) ) (cid:12)(cid:12)(cid:14)(cid:12)(cid:12) z ( F ( n +1) ) (cid:12)(cid:12) < ε , terminate. Otherwise, set n ← n + 1 and go to step 2.In the following, we show that Algorithm 2 converges to the KKT point of Problem (16). Theorem 1:
The sequences of { F ( n ) , n = 1 , , · · · } generated by Algorithm 2 converge to theKKT optimum point of Problem (16). Proof:
The proof is similar to that of [44] and hence it is omitted for simplicity. (cid:3)
Next, we briefly analyze the complexity of Algorithm 2. We assume that N B ≥ N I ≥ d . Ineach iteration of Algorithm 2, the main complexity contribution is the calculation of { F ( n +1) k , ∀ k } by using the bisection search method in Algorithm 1. In each iteration of Algorithm 1, the main complexity lies in calculating F in (24), which is on the order of O ( K I N B ) . The number ofiterations required for Algorithm 1 to converge is given by log (cid:0) λ u − λ l ε (cid:1) . Hence, the total complexityof Algorithm 1 is O (log (cid:0) λ u − λ l ε (cid:1) K I N B ) . Then, the total complexity of Algorithm 2 is given by O ( n max log (cid:0) λ u − λ l ε (cid:1) K I N B ) . C. Optimizing the Phase Shift Matrix
In this subsection, we focus our attention on optimizing the phase shift matrix Φ , while fixingthe other parameters. Upon substituting E k in (11) into (13) and removing the terms that areindependent of Φ , the phase shift optimization problem is formulated as: min Φ K I (cid:88) k =1 tr (cid:16) ω k W k U H k ¯H k ˜F ¯H H k U k (cid:17) − K I (cid:88) k =1 tr (cid:0) ω k W k U H k ¯H k F k (cid:1) − K I (cid:88) k =1 tr (cid:0) ω k W k F H k ¯H H k U k (cid:1) (31a)s.t. (8 c ) , (8 d ) , (31b)where ˜F = (cid:80) K I m =1 F m F H m .By substituting ¯H k = H b,k + H r,k ΦZ into (31a), we have ω k W k U H k ¯H k ˜F ¯H H k U k = ω k W k U H k H r,k ΦZ˜FZ H Φ H H H r,k U k + ω k W k U H k H b,k ˜FZ H Φ H H H r,k U k + ω k W k U H k H r,k ΦZ˜FH H b,k U k + ω k W k U H k H b,k ˜FH H b,k U k , (32)and ω k W k U H k ¯H k F k = ω k W k U H k H r,k ΦZF k + ω k W k U H k H b,k F k . (33)Let us define B k ∆ = ω k H H r,k U k W k U H k H r,k , C ∆ = Z˜FZ H and D k ∆ = ω k Z˜F H H H b,k U k W k U H k H r,k .By using (32), we arrive at: tr (cid:16) ω k W k U H k ¯H k ˜F ¯H H k U k (cid:17) = tr (cid:0) Φ H B k ΦC (cid:1) + tr (cid:0) Φ H D H k (cid:1) + tr ( ΦD k ) + const , (34)where const is a constant term that is independent of Φ .Similarly, by defining T k ∆ = ω k ZF k W k U H k H r,k , from (33) we have tr (cid:0) ω k W k U H k ¯H k F k (cid:1) = tr ( ΦT k ) + const , (35)where const is a constant term that is independent of Φ .By defining G b ∆ = (cid:80) K E l =1 α l η G H b,l G b,l , G r ∆ = (cid:80) K E l =1 α l η G H r,l G r,l , and G br ∆ = Z˜F (cid:80) K E l =1 α l η G H b,l G r,l ,the EH constraint in (8c) can be recast as follows: tr (cid:0) Φ H G r ΦC (cid:1) + tr (cid:0) Φ H G H b,r (cid:1) + tr ( ΦG br ) + tr (cid:16) G b ˜F (cid:17) ≥ ¯ Q. (36) By inserting (34) and (35) into the OF of Problem (31) and removing the constant terms, wehave min Φ tr (cid:0) Φ H BΦC (cid:1) + tr (cid:0) Φ H V H (cid:1) + tr ( ΦV ) (37a)s.t. (8 d ) , (36) , (37b)where B and V are given by B = (cid:80) K I k =1 B k and V = (cid:80) K I k =1 D k − (cid:80) K I k =1 T k , respectively.Upon denoting the collection of diagonal elements of Φ by φ = [ φ , · · · , φ M ] T and adoptingthe matrix identity of [45, Eq. (1.10.6)], it follows that tr (cid:0) Φ H BΦC (cid:1) = φ H (cid:0) B (cid:12) C T (cid:1) φ , tr (cid:0) Φ H G r ΦC (cid:1) = φ H (cid:0) G r (cid:12) C T (cid:1) φ . (38)Upon denoting the collections of diagonal elements of V and G br by v = (cid:104) [ V ] , , · · · , [ V ] M,M (cid:105) T and g = (cid:104) [ G br ] , , · · · , [ G br ] M,M (cid:105) T , we arrive at tr ( ΦV ) = v T φ , tr (cid:0) Φ H V H (cid:1) = φ H v ∗ , tr ( ΦG br ) = g T φ , tr (cid:0) Φ H G H br (cid:1) = φ H g ∗ . (39)Moreover, the constraint (36) can be rewritten as φ H Υ φ + 2Re (cid:8) φ H g ∗ (cid:9) ≥ (cid:95) Q, (40)where we have (cid:95) Q = ¯ Q − Tr (cid:16) G b ˜F (cid:17) and Υ = G r (cid:12) C T . It can be verified that G r and C T arenon-negative semidefinite matrices. Then, according to [45], the Hadamard product G r (cid:12) C T (orequivalently Υ ) is also a semidefinite matrix.Thus, Problem (37) can be transformed as min φ φ H Ξ φ + 2Re (cid:8) φ H v ∗ (cid:9) (41a)s.t. (8 d ) , (40) , (41b)where we have Ξ = B (cid:12) 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 (40), Problem (41) is difficult to solve. To deal with thisconstraint, we again employ the SCA method [40]. Specifically, since φ H Υ φ is convex w.r.t. φ ,its lower bound can be obtained as follows: φ H Υ φ ≥ − φ ( n )H Υ φ ( n ) + 2Re (cid:2) φ H Υ φ ( n ) (cid:3) , (42) where φ ( n ) is obtained in the previous iteration. Then, constraint (40) is replaced by the followingconstraint (cid:2) φ H (cid:0) g ∗ + Υ φ ( n ) (cid:1)(cid:3) ≥ (cid:95) Q + φ ( n )H Υ φ ( n ) ∆ = ˆ Q, (43)which is a linear constraint. Then, Problem (41) then becomes min φ φ H Ξ φ + 2Re (cid:8) φ H v ∗ (cid:9) (44a)s.t. (8 d ) , (43) . (44b)In the following, we conceive the Majorization-Minimization (MM) algorithm [31] for solvingProblem (44). The key idea is to solve a challenging problem by introducing a series of moretractable subproblems. Upon denoting the objective function of Problem (44) by f ( φ ) , in the ( n + 1) th iteration we have to find the upper bound of the OF, denoted as g ( φ | φ ( n ) ) , which shouldsatisfy the following three conditions: g ( φ ( n ) | φ ( n ) ) = f ( φ ( n ) ); 2) ∇ φ ∗ g ( φ | φ ( n ) ) (cid:12)(cid:12) φ = φ ( n ) = ∇ φ ∗ f ( φ ) | φ = φ ( n ) ; 3) g ( φ | φ ( n ) ) ≥ f ( φ ) . (45)Then, we solve the approximate subproblem defined by a more tractable new OF g ( φ | φ ( n ) ) . Tofind g ( φ | φ ( n ) ) , we introduce the following lemma [46]. Lemma 2:
For any given φ ( n ) , the following inequality holds for any feasible φ : φ H Ξ φ ≤ φ H X φ − (cid:8) φ H ( X − Ξ ) φ ( n ) (cid:9) + (cid:0) φ ( n ) (cid:1) H ( X − Ξ ) φ ( n ) ∆ = y ( φ | φ ( n ) ) , (46)where X = λ max I M and λ max is the maximum eigenvalue of Ξ . (cid:3) Then, the function g ( φ | φ ( n ) ) can be constructed as follows: g ( φ | φ ( n ) ) = y ( φ | φ ( n ) ) + 2Re (cid:8) φ H v ∗ (cid:9) , (47)where y ( φ | φ ( n ) ) is defined in (46). The new OF g ( φ | φ ( n ) ) is more tractable than the original OF f ( φ ) . The subproblem to be solved is given by min φ g ( φ | φ ( n ) ) (48a)s.t. (8 d ) , (43) . (48b)Since φ H φ = M , we have φ H X φ = M λ max , which is a constant. By removing the other constants,Problem (48) can be rewritten as follows: max φ (cid:8) φ H q ( n ) (cid:9) (49a)s.t. (8 d ) , (43) , (49b) where q ( n ) = ( λ max I M − Ξ ) φ ( n ) − v ∗ . Due to the additional constraint (43), the optimal solutionof Problem (49) cannot be obtained as in [24]. Furthermore, due to the non-convex unit-modulusconstraint (8d), Problem (49) is a non-convex optimization problem. As a result, the Lagrangiandual decomposition method developed for the convex problem (19) is not applicable here, sincethe dual gap is not zero.In the following, we propose a price mechanism for solving Problem (49) that can obtain theglobally optimal solution. Specifically, we consider the following problem by introducing a non-negative price p on the left hand side of constraint (43): max φ (cid:8) φ H q ( n ) (cid:9) + 2 p Re (cid:2) φ H (cid:0) g ∗ + Υ φ ( n ) (cid:1)(cid:3) (50a)s.t. (8 d ) . (50b)For a given p , the globally optimal solution is given by φ ( p ) = e j arg ( q ( n ) + p ( g ∗ + Υ φ ( n ) )) . (51)Our objective is to find a p value for ensuring that the complementary slackness condition forconstraint (43) is satisfied: p (cid:16) J ( p ) − ˆ Q (cid:17) = 0 , (52)where J ( p ) = 2Re (cid:104) φ ( p ) H (cid:0) g ∗ + Υ φ ( n ) (cid:1)(cid:105) . To solve this equation, we consider two cases: 1) p = 0 ;2) p > . Case I:
In this case, φ (0) = e j arg ( q ( n ) ) has to satisfy constraint (43). Otherwise, p > . Case II:
Since p > , equation (52) holds only when J ( p ) = ˆ Q . To solve this equation, we firstprovide 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. (cid:3)
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 (49) inAlgorithm 3. Although Problem (49) is a non-convex problem, in the following theorem we provethat Algorithm 3 is capable of finding the globally optimal solution. Theorem 2:
Algorithm 3 is capable of finding the globally optimal solution of Problem (49)and thus also of Problem (48).
Proof:
Please refer to Appendix B. (cid:3) Algorithm 3
Bisection Search Method to Solve Problem (49) Calculate J (0) . If J (0) ≤ ˆ Q , terminate. Otherwise, go to step 2. Initialize the accuracy ε , bounds p l and p u ; Calculate p = ( p l + p u )/2 ; Update φ ( p ) in (51) and calculate J ( p ) ; If J ( p ) ≥ ˆ Q , set p u = p ; Otherwise, set p l = p ; If | p l − p u | ≤ ε , terminate; Otherwise, go to step 3. Algorithm 4
MM Combined with SCA Algorithm to Solve Problem (31) Initialize the accuracy ε , the phase shifts φ (0) , the iteration index to n = 0 , the maximumnumber of iterations to n max , calculate the OF value of Problem (44) as f ( φ (0) ) ; Calculate ˆ Q ( n ) = (cid:95) Q + φ ( n )H Υ φ ( n ) ; Calculate q ( n ) = ( λ max I M − Ξ ) φ ( n ) − v ∗ ; Update φ ( n +1) by solving Problem (49) using Algorithm 3; If n ≥ n max or (cid:12)(cid:12) f ( φ ( n +1) ) − f ( φ ( n ) ) (cid:12)(cid:12)(cid:14) f ( φ ( n +1) ) ≤ ε holds, terminate; Otherwise, set n ← n + 1 and go to step 2.Based on the above, we now provide the details of solving Problem (31) in Algorithm 4.In the following theorem, we prove that the sequence of { φ ( n ) , n = 1 , , · · · } generated byAlgorithm 4 converges to the KKT-optimal point of Problem (31). Theorem 3:
The sequences of the OF value produced by Algorithm 4 are guaranteed to converge,and the final solution satisfies the KKT point of Problem (31).
Proof:
Please refer to Appendix C. (cid:3)
Let us now analyze the complexity of Algorithm 4. The complexity is dominated by calculating φ ( n +1) in step 4 using Algorithm 3. The complexity mainly depends on calculating the maximumeigenvalue of Ξ . Its complexity is on the order of O ( M ) . The number of iterations required forAlgorithm 3 is log (cid:0) p u − p l ε (cid:1) . Then, the total complexity of step 3 is O (log (cid:0) p u − p l ε (cid:1) M ) . Hence,the total complexity of Algorithm 4 is given by O ( n max log (cid:0) p u − p l ε (cid:1) M ) . D. Overall Algorithm to Solve Problem (8)
Based on the above analysis, we provide the detailed steps of the BCD algorithm to solveProblem (8) in Algorithm 5, where R ( F ( n ) , φ ( n ) ) denotes the OF value of Problem (8) in the n th iteration. Algorithm 5
Block Coordinate Descent Algorithm Initialize iterative number n = 1 , maximum number of iterations n max , feasible F (1) , φ (1) ,error tolerance ε , calculate R ( F (1) , φ (1) ) , calculate the optimal decoding matrices U (1) andauxiliary matrices W (1) based on (14); Given U ( n ) , W ( n ) and φ ( n ) , calculate the optimal precoding matrices F ( n +1) by solving Problem(16) using Algorithm 2; Given U ( n ) , W ( n ) and F ( n +1) , calculate the optimal φ ( n +1) by solving Problem (31) usingAlgorithm 4; Given F ( n +1) and φ ( n +1) , calculate the optimal decoding matrices U ( n +1) in (14); Given F ( n +1) , U ( n +1) and φ ( n +1) , calculate the optimal auxiliary matrices W ( n +1) in (14); If n ≥ n max or (cid:12)(cid:12) R ( F ( n +1) , φ ( n +1) ) − R ( F ( n ) , φ ( n ) ) (cid:12)(cid:12)(cid:14) R ( F ( n +1) , φ ( n +1) ) < ε , terminate.Otherwise, set n ← n + 1 and go to step 2.The following theorem shows the convergence and solution properties of Algorithm 5. Theorem 4:
The OF value sequence { R ( F ( n ) , φ ( n ) ) , n = 1 , , · · · } generated by Algorithm 5 isguaranteed to converge, and the final solution satisfies the KKT conditions of Problem (8). Proof:
Please refer to Appendix D. (cid:3)
The complexity of Algorithm 5 mainly depends on that of Step 2 and Step 3, the complexity ofwhich has been analyzed in the above subsections. In specific, the total complexity of step 2 and step3 are respectively given by O ( n max1 log (cid:0) λ u − λ l ε (cid:1) K I N B ) and O ( n max2 log (cid:0) p u − p l ε (cid:1) M ) , where n max1 and n max2 denote the number of iterations for Algorithm 2 and Algorithm 4 to converge. Denote thetotal number of iterations of Algorithm 5 as N max . Then, the overall complexity of Algorithm 5 isgiven by O (cid:0) N max (cid:0) n max1 log (cid:0) λ u − λ l ε (cid:1) K I N B + n max2 log (cid:0) p u − p l ε (cid:1) M (cid:1)(cid:1) . Additionally, the simulationresults show that Algorithm 5 converges rapidly, which demonstrates the low complexity of thisalgorithm. IV. F EASIBILITY C HECK FOR P ROBLEM (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: max F , Φ tr (cid:32) K I (cid:88) k =1 F H k GF k (cid:33) (53a)s.t. (8 b ) , (8 d ) . (53b)If the optimal OF value is larger than ¯ Q , Problem (8) is feasible. Otherwise, it is infeasible. Asthe TPC matrices and phase shift matrix are coupled, the globally optimal solution is difficult toobtain. In the following, we can obtain a suboptimal solution by alternatively optimizing the TPCmatrices and phase shifts.For a given phase shift matrix, the TPC matrix optimization problem is given by max F tr (cid:32) K I (cid:88) k =1 F H k GF k (cid:33) (54a)s.t. (8 b ) . (54b)Upon denoting the maximum eigenvalue and the corresponding eigenvector of G by χ and b respectively, the optimal solution can be readily obtained as F k = (cid:2) √ p k b , N B × ( d − (cid:3) , ∀ k =1 , · · · , K I , where (cid:80) K I k =1 p k = P T and p k ≥ , ∀ k = 1 , · · · , K I . Without loss of generality, we canset p i = P T /K I , ∀ i ∈ K I . The OF value is given by χP T . In this case, the optimal TPC matrixrepresents the optimal energy beamforming, which is the same as that for the single-antenna IRcase of [38].For a given TPC matrix F , the phase shift optimization problem is formulated as: max φ φ H Υ φ + 2Re (cid:8) φ H g ∗ (cid:9) (55a)s.t. (8 d ) , (55b)where Υ and g are defined in the above section. The OF is convex w.r.t. φ , and maximizing a convexfunction 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, alternatively solve Problem (54) and (55) until the OF is larger than ¯ Q .V. S IMULATION R ESULTS
In this section, we provide simulation results for demonstrating the benefits of applying IRS toSWIPT 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 , and ( x IR , with radius 1 m and4 m, respectively. The IRS is located at ( x IRS , . In the simulations, we assume that the IRS isjust above the ERs and thus we set x ER = x IRS . The large-scale path loss is modeled in dB as
P L = P L (cid:18) DD (cid:19) − α , (56)where PL is the path loss at the reference distance D , D is the link length in meters, and α is the path loss exponent. Here, we set D = 1 and PL = − . The path loss exponentsof the BS-IRS, IRS-ER, IRS-IR, BS-IR and BS-ER links are respectively set as α BSIRS = 2 . , α IRSER = 2 . , α IRSIR = 2 . , α BSIR = 3 . and α BSER = 3 . . Unless otherwise stated, the otherparameters are set as follows: Channel bandwidth of 1 MHz, noise power density of − dBm/Hz, N B = 4 , N I = N E = 2 , d = 2 , ¯ Q = 2 × − W, η = 0 . , M = 50 , P T = 10 W , weight factors ω k = 1 , ∀ k ∈ K I , α l = 1 , ∀ l ∈ K E , x ER = 5 m, and x IR = 400 m. The following results areobtained by averaging over 100 random locations and channel generations. Due to the severeblockage and long distance, the channels from the BS and the IRS to the IRs are assumed to beRayleigh fading. However, as the BS, the ERs and the IRS are close to each other, the small-scalechannels are assumed to be Rician fading. In particular, the small-scale channels from the IRS tothe ERs are denoted as: ˜G r,l = (cid:115) β irser β irser + 1 ˜G LoS r,l + (cid:114) β irser + 1 ˜G NLoS r,l , l = 1 , · · · , K E , (57)where β irser is the Rician factor, ˜G LoS r,l is the deterministic line of sight (LoS), and ˜G NLoS r,l is thenon-LoS (NLoS) component that is Rayleigh fading. The LoS component ˜G LoS r,l can be modeledas ˜G LoS r,l = a N E (cid:0) ϑ AoA irser ,l (cid:1) a HM (cid:0) ϑ AoD irser ,l (cid:1) , where a N E (cid:0) ϑ AoA irser ,l (cid:1) is defined as a N E (cid:0) ϑ AoA irser ,l (cid:1) = (cid:104) , e j πdλ sin ϑ AoA irser ,l , · · · , e j πdλ ( N E −
1) sin ϑ AoA irser ,l (cid:105) T (58)and a M (cid:0) ϑ AoD irser ,l (cid:1) = (cid:104) , e j πdλ sin ϑ AoD irser ,l , · · · , e j πdλ ( M −
1) sin ϑ AoD irser ,l (cid:105) T . (59)In (58) and (59), d is the antenna separation distance, λ is the wavelength, ϑ AoD irser ,l is the angleof departure and ϑ AoA irser ,l is the angle of arrival. It is assumed that ϑ AoD irser ,l and ϑ AoA irser ,l are randomlydistributed within [0 , π ] . For simplicity, we set d/λ = 1 / . The small-scale channels from theBS to the ERs and the IRS are similarly defined. For simplicity, the Rician factors for all Ricianfading channels are assumed to be the same as β = 3 . BS (0,0) (m) x (m) y ER x IR x IRS ( ,2) x IRSER area IR area
Fig. 2. The simulated IRS-aided SWIPT MIMO communication scenario. −3 x ER (m) H a v e s t ed P o w e r ( W ) M=10M=20M=40No−IRS
Fig. 3. Maximum harvested power achieved by various schemes. W S R ( b i t/ s / H z ) M=10M=20M=40
Fig. 4. Convergence behaviour of the BCD algorithm.
We first compare the maximum power harvested by various schemes in Fig. 3. Specifically, wesolve the EH maximization problem (53) by using the feasibility check method in Section IV.Additionally, we also present the results without using IRS. Fig. 3 shows the maximum EH powerversus the ER circle center location x EH . As expected, the EH power gleaned by all schemesdecreases, when the ERs are far away from the BS. As expected, more power can be harvestedwith the aid of IRS than that without IRS, especially when the number of phase shifters M islarge. This is mainly due to the fact that an additional strong link is reflected by the IRS, whichcan be harvested by the ERs. This figure also shows that the IRS is effective in expanding theoperational range of ERs. For example, when the harvested power limit is ¯ Q = 2 × − W, themaximum 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 numbersof phase shifters M . It is observed from Fig. 4 that the WSR achieved for various M valuesincreases 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 BCDalgorithm to achieve a large portion of the converged WSR. This reflects the low complexity of
20 40 60 80 1001618202224262830 Number of Phase Shifters M W S R ( b i t/ s / H z ) BCD Alg.Fixed PhaseNo−IRS
Fig. 5. WSR versus the number of phase shifters. −4 W S R ( b i t/ s / H z ) BCD Alg.Fixed PhaseNo−IRS
Fig. 6. WSR versus the harvested power requirement ¯ Q . 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 theharvested power maximization problem (53), while they are not optimized, when using the BCDalgorithm by removing Step 3 of the BCD algorithm. When any of the methods fails to obtain afeasible 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 ofvarious algorithms. As expected, the WSR achieved by all the algorithms - except for the No-IRSmethod - increases with M , since a higher degree of freedom can be exploited for optimizing thesystem performance. By carefully optimizing the phase shifts for maximizing the WSR, the BCDalgorithm significantly outperforms the fixed-phase scheme. Additionally, the performance gainincreases with M , which emphasizes the importance of optimizing the phase shifts. By employingthe IRS in our SWIPT system, the WSR obtained by the BCD algorithm becomes drasticallyhigher than that of No-IRS. For example, when M = 60 , the WSR performance gain is up to 10bit/s/Hz. These results demonstrate that introducing the IRS into our SWIPT system is a promisingtechnique of enhancing the system performance.In Fig. 6, the impact of harvested power requirement ¯ Q is investigated. It is seen from this figurethat the WSR achieved by all the algorithms decreases upon increasing ¯ Q , because the probabilityof infeasibility increases, which in turn reduces the average WSR value. We also find that theWSR obtained by the No-IRS scheme decreases more rapidly than that of the other two IRS-aidedtransmission schemes. The WSR of the No-IRS is approaching zero when ¯ Q = 4 × − W, while α IRS W S R ( b i t/ s / H z ) BCD Alg.Fixed PhaseNo−IRS
Fig. 7. WSR versus the IRS-related path loss exponent α IRS . ER (m) W S R ( b i t/ s / H z ) BCD Alg.Fixed PhaseNo−IRS
Fig. 8. WSR versus the location of ER circle center x ER . those relying on IRSs achieve a WSR gain in excess of 20 bit/s/Hz. It is observed again that theBCD algorithm performs better than the fixed-phase scheme, but the gap narrows with the increaseof ¯ Q . This can be explained as follows. With the increase of ¯ Q , both the TPC matrices and thephase shifts should be designed for maximizing the power harvested at the ERs, and thus the finalsolutions of the fixed-phase and BCD method will become the same.The above results are obtained for α BSIRS = 2 . , α IRSER = 2 . , α IRSIR = 2 . based on theassumption that the IRS relies on an obstacle-free scenario. In practice, this ideal scenario is seldomencountered. 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 WSRachieved 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 uponincreasing α IRS , the signal power reflected from the IRS becomes weaker. Hence, the benefits ofthe IRS can be eroded. This provides an important engineering design insight: the location of IRSshould 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, theWSR achieved by all the schemes decreases with x IRS , since the ERs become more distant from theBS and the signals gleaned from both the BS and IRS become weaker. The WSR achieved by theNo-IRS approaches zero when x IRS = 8 m, hence this method cannot reach the energy transmissiontarget of the ERs. The proposed algorithm is again observed to significantly outperform the othertwo 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 achievedby all the algorithms decreases with x IR since the IRs become farther away from the BS when
100 200 300 400 500 600101520253035404550 x IR (m) W S R ( b i t/ s / H z ) BCD Alg.Fixed PhaseNo−IRS
Fig. 9. WSR versus the location of IR circle center x IR . increasing x IR . The proposed algorithm is shown to achieve nearly the WSR gain of 10 bit/s/Hzover the No-IRS when x IR = 100 m , and the WSR gain slightly increases with x IR . This meansthat the IRS is more advantageous when the IRs are far away from the BS, and the IRS can provideone additional favorable link. VI. C ONCLUSIONS
In this paper, we have invoked an IRS in a SWIPT MIMO system for enhancing the performanceof both the ERs and IRs. By carefully adjusting the phase shifts at the IRS, the signal reflectedby the IRS can be added constructively at both the ERs and IRs. We considered the WSRmaximization problem of IRs, while guaranteeing the energy harvesting requirements of the ERsand the associated non-convex unit-modulus constraints. We conceived a BCD algorithm foralternatively optimizing the TPC matrices at the BS and the phase shift matrix at the IRSs. Foreach subproblem, a low-complexity iterative algorithm was proposed, which guarantees to be atworst locally optimal. Our simulation results demonstrated that the IRS enhances the performanceof the SWIPT system and that the proposed algorithm converges rapidly, hence it is eminentlysuitable 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 beleft for future work. A PPENDIX AP ROOF OF L EMMA λ and λ (cid:48) , where λ > λ (cid:48) . Let F ( λ ) and F ( λ (cid:48) ) be the optimalsolutions of Problem (21) with λ and λ (cid:48) , respectively. Since F ( λ ) is the optimal solution of Problem(21) with λ , we have L [ F ( λ ) , λ ] ≤ L [ F ( λ (cid:48) ) , λ ] . (A.1)Similarly, we have L [ F ( λ (cid:48) ) , λ (cid:48) ] ≤ L [ F ( λ ) , λ (cid:48) ] . (A.2)By adding these two inequalities and simplifying them, we have ( λ − λ (cid:48) ) P ( λ ) ≤ ( λ − λ (cid:48) ) P ( λ (cid:48) ) .Since λ > λ (cid:48) , we have P ( λ ) ≤ P ( λ (cid:48) ) , which completes the proof.A PPENDIX BP ROOF OF T HEOREM φ (cid:63) . According to [43], for a non-convex optimization problem, all its locally optimal solutions (including the globally optimalsolution) should satisfy the Karush-Kuhn-Tucker (KKT) optimality conditions, one of which isthe complementary slackness condition for constraint (43): λ (cid:63) (cid:16) (cid:2) φ (cid:63) H (cid:0) g ∗ + Υ φ ( n ) (cid:1)(cid:3) − ˆ Q (cid:17) = 0 , (B.1)where λ (cid:63) is the corresponding optimal dual variable. We consider two cases: 1) λ (cid:63) = 0 ; 2) λ (cid:63) > .The first case means that constraint (43) is not tight in the optimum. Then, the optimal solutioncan be obtained as φ (cid:63) = e j arg ( q ( n ) ) , which is equal to φ (0) . Hence, Algorithm 3 achieves theoptimal solution of Problem (49).For the second case, the following equality should hold: (cid:2) φ (cid:63) H (cid:0) g ∗ + Υ φ ( n ) (cid:1)(cid:3) = ˆ Q. (B.2)We prove the second case by using the method of contradiction. Denote the optimal p obtainedby Algorithm 3 as p (cid:63) , and the corresponding φ as φ ( p (cid:63) ) . Then, we have (cid:2) φ ( p (cid:63) ) H (cid:0) g ∗ + Υ φ ( n ) (cid:1)(cid:3) = ˆ Q. (B.3) Let us assume that φ ( p (cid:63) ) is not the globally optimal solution of Problem (49). Then, we have (cid:8) φ ( p (cid:63) ) H q ( n ) (cid:9) < (cid:8) φ (cid:63) H q ( n ) (cid:9) . (B.4)Since φ ( p (cid:63) ) is the globally optimal solution of Problem (50) when p = p (cid:63) , we have (cid:8) φ ( p (cid:63) ) H q ( n ) (cid:9) +2 p (cid:63) Re (cid:2) φ ( p (cid:63) ) H (cid:0) g ∗ + Υ φ ( n ) (cid:1)(cid:3) ≥ (cid:8) φ (cid:63) H q ( n ) (cid:9) +2 p (cid:63) Re (cid:2) φ (cid:63) H (cid:0) g ∗ + Υ φ ( n ) (cid:1)(cid:3) . (B.5)By substituting (B.2) and (B.3) into (B.5), we have (cid:8) φ ( p (cid:63) ) H q ( n ) (cid:9) ≥ (cid:8) φ (cid:63) H q ( n ) (cid:9) , (B.6)which contradicts (B.4). Hence, the solution obtained by Algorithm 3 is the globally optimalsolution of Problem (49). Since Problem (48) is equivalent to Problem (49), the proof is complete.A PPENDIX CP ROOF OF T HEOREM T ( φ ) ∆ = φ H Υ φ + 2Re (cid:8) φ H g ∗ (cid:9) + tr (cid:16) G b ˜F (cid:17) , (C.1) ¯ T ( φ | φ ( n ) ) ∆ = − φ ( n )H Υ φ ( n ) + 2Re (cid:2) φ H (cid:0) g ∗ + Υ φ ( n ) (cid:1)(cid:3) + tr (cid:16) G b ˜F (cid:17) . (C.2)It can be verified that T ( φ ( n ) ) = ¯ T ( φ ( n ) | φ ( n ) ) .We first show that the solution sequence { φ ( n ) , n = 1 , , · · · } is feasible for Problem (31). Theunit-modulus constraint is guaranteed in (51). We only have to check the EH constraint in (8c).Note that φ ( n +1) is a feasible solution of Problem (49), and thus satisfies constraint (42). Hence,we have ¯ T ( φ ( n +1) | φ ( n ) ) ≥ ¯ Q . By using inequality (42), we have T ( φ ( n +1) ) ≥ ¯ T ( φ ( n +1) | φ ( n ) ) .Then, T ( φ ( n +1) ) ≥ ¯ Q holds, which means that the sequence of φ ( n +1) satisfies the EH constraintin (8c).Now, we show that the OF value sequence { f ( φ ( n ) ) , n = 1 , , · · · } is monotonically decreas-ing. Based on Theorem 2, the globally optimal solution Φ to Problem (48) can be obtained.Then, we have g ( φ ( n +1) | φ ( n ) ) ≤ g ( φ ( n ) | φ ( n ) ) . According to the first condition in (45), we have g ( φ ( n ) | φ ( n ) ) = f ( φ ( n ) ) . Hence, we have g ( φ ( n +1) | φ ( n ) ) ≤ f ( φ ( n ) ) . By using the third condition of(45), we have g ( φ ( n +1) | φ ( n ) ) ≥ f ( φ ( n +1) ) . As a result, we have f ( φ ( n ) ) ≥ f ( φ ( n +1) ) . Additionally,the OF must have a lower bound due to the unit-modulus constraint. Hence, the OF value sequence { f ( φ ( n ) ) , n = 1 , , · · · } is guaranteed to converge.Now, we prove that the converged solution satisfies the KKT conditions of Problem (31). Letus denote the converged solution by { φ (cid:63) } . Since φ (cid:63) is the globally optimal solution of Problem (48), it must satisfy the KKT conditions of Problem (48). Specifically, the Lagrange function ofProblem (48) is given by L ( φ , ν, τ ) = g ( φ | φ (cid:63) ) + ν (cid:16) ˆ Q − (cid:2) φ H ( g ∗ + Υ φ (cid:63) ) (cid:3)(cid:17) + M (cid:88) m =1 τ m ( | φ m | − , (C.3)where ν and τ = { τ , · · · , τ M } are the corresponding dual variables. Then, there must exist a ν (cid:63) and τ (cid:63) = { τ (cid:63) , · · · , τ (cid:63)M } for ensuring that the following conditions are satisfied: ∇ φ ∗ L ( φ , ν, τ ) | φ = φ (cid:63) = ∇ φ ∗ g ( φ | φ (cid:63) ) | φ = φ (cid:63) − ν (cid:63) ( g ∗ + Υ φ (cid:63) )+ M (cid:88) m =1 τ (cid:63)m ( ∇ φ ∗ | φ m | ) | φ = φ (cid:63) = , (C.4) ν (cid:63) (cid:16) ˆ Q − (cid:2) φ (cid:63) H ( g ∗ + Υ φ (cid:63) ) (cid:3)(cid:17) = 0 , (C.5) τ (cid:63)m ( | φ (cid:63)m | −
1) = 0 , ∀ m. (C.6)According to the second condition of (45), we have ∇ φ ∗ g ( φ | φ (cid:63) ) | φ = φ (cid:63) = ∇ φ ∗ f ( φ ) | φ = φ (cid:63) . (C.7)Upon denoting the OF of Problem (31) as ϕ ( φ ) , which is the same as f ( φ ) except that ϕ ( φ ) has more constants, we have ∇ φ ∗ f ( φ ) | φ = φ (cid:63) = ∇ φ ∗ ϕ ( φ ) | φ = φ (cid:63) . Combining with (C.7), we have ∇ φ ∗ g ( φ | φ (cid:63) ) | φ = φ (cid:63) = ∇ φ ∗ ϕ ( φ ) | φ = φ (cid:63) . By substituting it into (C.4), we arrive at ∇ φ ∗ ϕ ( φ ) | φ = φ (cid:63) − ν (cid:63) ( g ∗ + Υ φ (cid:63) ) + M (cid:88) m =1 τ (cid:63)m ( ∇ φ ∗ | φ m | ) | φ = φ (cid:63) = . (C.8)It can be checked that the set of equations (C.5), (C.6) and (C.8) constitutes exactly the KKTconditions of Problem (31). Hence, the proof is complete.A PPENDIX DP ROOF OF T HEOREM h ( W , U , F , Φ ) ∆ = K I (cid:88) k =1 ω k h k ( W , U , F , Φ ) . (D.1)It can be readily verified that the sequence of solutions { F ( n ) , φ ( n ) } generated by Algorithm 5 isalways feasible for Problem (8). The monotonic property of Algorithm 5 can be similarly provedby using the method of [33].In the following, we prove that the converged solution satisfies the KKT conditions of Problem(8). Let us denote the converged solution as { W (cid:63) , U (cid:63) , F (cid:63) , Φ (cid:63) } .According to Theorem 1, F (cid:63) is the KKT-optimum point of Problem (16). Upon denoting theOF of Problem (16) as z ( F , Φ (cid:63) ) , the Lagrange function of Problem (16) is given by L ( F , λ, µ ) = z ( F , Φ (cid:63) ) + λ (cid:32) K I (cid:88) k =1 (cid:107) F k (cid:107) F − P T (cid:33) + µ (cid:32) ¯ Q − tr (cid:32) K I (cid:88) k =1 F H k GF k (cid:33)(cid:33) , (D.2) where λ and µ are the corresponding dual variables. Then, there must exist a λ (cid:63) and µ (cid:63) for ensuringthat the following conditions are satisfied : ∇ F ∗ k L ( F , λ, µ ) (cid:12)(cid:12) F k = F (cid:63)k = ∇ F ∗ k z ( F , Φ (cid:63) ) (cid:12)(cid:12) F k = F (cid:63)k + λ (cid:63) F (cid:63)k − µ (cid:63) GF (cid:63)k = , ∀ k ∈ K I , (D.3) λ (cid:63) (cid:32) K I (cid:88) k =1 (cid:107) F (cid:63)k (cid:107) F − P T (cid:33) = 0 , (D.4) µ (cid:63) (cid:32) ¯ Q − tr (cid:32) K I (cid:88) k =1 F (cid:63) H k GF (cid:63)k (cid:33)(cid:33) = 0 . (D.5)Furthermore, it can be readily checked that ∇ F ∗ k h ( W (cid:63) , U (cid:63) , F , Φ (cid:63) ) (cid:12)(cid:12) F k = F (cid:63)k = ∇ F ∗ k z ( F , Φ (cid:63) ) (cid:12)(cid:12) F k = F (cid:63)k , ∀ k ∈ K I . (D.6)To expound a little further, we have the following chain of inequalities: ∇ F ∗ k h k ( W (cid:63) , U (cid:63) , F , Φ (cid:63) ) (cid:12)(cid:12) F k = F (cid:63)k (D.7) = − tr (cid:16) W (cid:63)k (cid:16) ∇ F ∗ k E k ( U (cid:63) , F , Φ (cid:63) ) (cid:12)(cid:12) F k = F (cid:63)k (cid:17)(cid:17) (D.8) = − tr (cid:16) ( E k ( U (cid:63) , F (cid:63) , Φ (cid:63) )) − (cid:16) ∇ F ∗ k E k ( U (cid:63) , F , Φ (cid:63) ) (cid:12)(cid:12) F k = F (cid:63)k (cid:17)(cid:17) (D.9) = (cid:0) ∇ F ∗ k log (cid:12)(cid:12) ( E k ( U (cid:63) , F , Φ (cid:63) )) − (cid:12)(cid:12)(cid:1)(cid:12)(cid:12) F k = F (cid:63)k (D.10) = ∇ F ∗ k R k ( F , Φ (cid:63) ) (cid:12)(cid:12) F k = F (cid:63)k , (D.11)where (D.8) follows from the chain rule, and the final equality follows from applying the Woodburymatrix identity to (15). Combining (D.11) with (D.6), we have ∇ F ∗ k z ( F , Φ (cid:63) ) (cid:12)(cid:12) F k = F (cid:63)k = ∇ F ∗ k R k ( F , Φ (cid:63) ) (cid:12)(cid:12) F k = F (cid:63)k . (D.12)By substituting (D.12) into (D.3), we arrive at ∇ F ∗ k R k ( F , Φ (cid:63) ) (cid:12)(cid:12) F k = F (cid:63)k + λ (cid:63) F (cid:63)k − µ (cid:63) GF (cid:63)k = , ∀ k ∈ K I . (D.13)According to Theorem 3, φ (cid:63) satisfies the KKT conditions of Problem (31), and thus the set ofequations (C.5), (C.6) and (C.8) hold.Furthermore, it can be readily verified that ∇ φ ∗ h ( W (cid:63) , U (cid:63) , F (cid:63) , Φ ) | φ = φ (cid:63) = ∇ φ ∗ ϕ ( φ ) | φ = φ (cid:63) . (D.14)By using similar derivations as in (D.7)-(D.11), we can prove that ∇ φ ∗ h ( W (cid:63) , U (cid:63) , F (cid:63) , Φ ) | φ = φ (cid:63) = ∇ φ ∗ R k ( φ , F (cid:63) ) | φ = φ (cid:63) . (D.15) For simplicity, the prime constraints are omitted. Hence, we have ∇ φ ∗ ϕ ( φ ) | φ = φ (cid:63) = ∇ φ ∗ R k ( φ , F (cid:63) ) | φ = φ (cid:63) . (D.16)By substituting (D.16) into (C.8), we arrive at: ∇ φ ∗ R k ( φ , F (cid:63) ) | φ = φ (cid:63) − ν (cid:63) ( g ∗ + Υ φ (cid:63) ) + M (cid:88) m =1 τ (cid:63)m ( ∇ φ ∗ | φ m | ) | φ = φ (cid:63) = . (D.17)Then, the set of equations (D.13), (D.4), (D.5), (D.17), (C.5), and (C.6) constitute exactly theKKT conditions of Problem (8). R EFERENCES [1] M. Di Renzo, M. Debbah, D.-T. Phan-Huy, A. Zappone, M.-S. Alouini, C. Yuen, V. Sciancalepore, G. C. Alexandropoulos,J. Hoydis, H. Gacanin et al. , “Smart radio environments empowered by reconfigurable AI meta-surfaces: an idea whose timehas come,”
EURASIP J. Wireless Commun. Networking , vol. 2019, no. 1, p. 129, 2019.[2] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network.”[Online]. Available: https://arxiv.org/abs/1905.00152[3] J. Zhang, E. Bj¨ornson, M. Matthaiou, D. W. K. Ng, H. Yang, and D. Love, “Multiple antenna technologies for beyond 5G.”[Online]. Available: https://arxiv.org/abs/1910.00092[4] T. J. Cui, M. Q. Qi, X. Wan, J. Zhao, and Q. Cheng, “Coding metamaterials, digital metamaterials and programmablemetamaterials,”
Light: Science & Applications , vol. 3, no. 10, p. e218, 2014.[5] J. Zhang, S. Chen, Y. Lin, J. Zheng, B. Ai, and L. Hanzo, “Cell-free massive MIMO: A new next-generation paradigm,”
IEEE Access , vol. 7, pp. 99 878–99 888, 2019.[6] J. Zhang, L. Dai, Z. He, B. Ai, and O. A. Dobre, “Mixed-ADC/DAC multipair massive MIMO relaying systems: Performanceanalysis and power optimization,”
IEEE Trans. Commun. , vol. 67, no. 1, pp. 140–153, 2018.[7] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network: Joint active and passive beamforming design,”in , Dec. 2018, pp. 1–6.[8] X. Yu, D. Xu, and R. Schober, “MISO wireless communication systems via intelligent reflecting surfaces.” [Online].Available: https://arxiv.org/abs/1904.12199[9] Y. Yang, B. Zheng, S. Zhang, and R. Zhang, “Intelligent reflecting surface meets OFDM: Protocol design and ratemaximization.” [Online]. Available: https://arxiv.org/abs/1906.09956[10] Y. Han, W. Tang, S. Jin, C. Wen, and X. Ma, “Large intelligent surface-assisted wireless communication exploiting statisticalCSI,”
IEEE Trans. Veh. Technol. , 2019.[11] S. Abeywickrama, R. Zhang, and C. Yuen, “Intelligent reflecting surface: Practical phase shift model and beamformingoptimization.” [Online]. Available: https://arxiv.org/abs/1907.06002[12] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,”
IEEE Trans. Wireless Commun. , vol. 18, no. 11, pp. 5394–5409, Nov. 2019.[13] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energyefficiency in wireless communication,”
IEEE Trans. Wireless Commun. , vol. 18, no. 8, pp. 4157–4170, Aug. 2019.[14] H. Guo, Y. Liang, J. Chen, and E. Larsson, “Weighted sum-rate optimization for intelligent reflecting surface enhancedwireless networks.” [Online]. Available: https://arxiv.org/abs/1905.07920 [15] Q. Nadeem, A. Kammoun, A. Chaaban, M. Debbah, and M. Alouini, “Large intelligent surface assisted MIMOcommunications.” [Online]. Available: https://arxiv.org/abs/1903.08127[16] X. Yu, D. Xu, and R. Schober, “Enabling secure wireless communications via intelligent reflecting surfaces.” [Online].Available: https://arxiv.org/abs/1904.09573[17] M. Cui, G. Zhang, and R. Zhang, “Secure wireless communication via intelligent reflecting surface,” IEEE Wireless Commun.Lett. , 2019.[18] H. Shen, W. Xu, W. Xu, S. Gong, Z. He, and C. Zhao, “Secrecy rate maximization for intelligent reflecting surface assistedmulti-antenna communications,”
IEEE Commun. Lett. , pp. 1–1, 2019.[19] J. Chen, Y. Liang, Y. Pei, and H. Guo, “Intelligent reflecting surface: A programmable wireless environment for physicallayer security.” [Online]. Available: https://arxiv.org/abs/1905.03689[20] D. Xu, X. Yu, Y. Sun, D. W. K. Ng, and R. Schober, “Resource allocation for secure IRS-assisted multiuser MISO systems.”[Online]. Available: https://arxiv.org/abs/1907.03085[21] X. Guan, Q. Wu, and R. Zhang, “Intelligent reflecting surface assisted secrecy communication via joint beamforming andjamming.” [Online]. Available: https://arxiv.org/abs/1907.12839[22] T. Bai, C. Pan, Y. Deng, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Latency minimization for intelligent reflectingsurface aided mobile edge computing.” [Online]. Available: https://arxiv.org/abs/1910.07990[23] G. Zhou, C. Pan, H. Ren, K. Wang, and A. Nallanathan, “Intelligent reflecting surface aided multigroup multicast MISOcommunication systems.” [Online]. Available: https://arxiv.org/abs/1909.04606[24] C. Pan, H. Ren, K. Wang, W. Xu, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Intelligent reflecting surface for multicellMIMO communications.” [Online]. Available: https://arxiv.org/abs/1907.10864[25] C. Huang, G. C. Alexandropoulos, C. Yuan, and M. Debbah, “Indoor signal focusing with deep learning designedreconfigurable intelligent surfaces.” [Online]. Available: https://arxiv.org/abs/1905.07726[26] Z. He and X. Yuan, “Cascaded channel estimation for large intelligent metasurface assisted massive MIMO.” [Online].Available: https://arxiv.org/abs/1905.07948[27] A. Taha, M. Alrabeiah, and A. Alkhateeb, “Enabling large intelligent surfaces with compressive sensing and deep learning.”[Online]. Available: https://arxiv.org/abs/1904.10136[28] G. Zhou, C. Pan, H. Ren, K. Wang, M. D. Renzo, and A. Nallanathan, “Robust beamforming design for intelligent reflectingsurface aided MISO communication systems.” [Online]. Available: https://arxiv.org/abs/1911.06237[29] R. Zhang and C. K. Ho, “MIMO broadcasting for simultaneous wireless information and power transfer,”
IEEE Trans. WirelessCommun. , vol. 12, no. 5, pp. 1989–2001, May 2013.[30] Q. Wu and R. Zhang, “Weighted sum power maximization for intelligent reflecting surface aided SWIPT.” [Online].Available: https://arxiv.org/abs/1907.05558[31] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machinelearning,”
IEEE Trans. Signal Process. , vol. 65, no. 3, pp. 794–816, Feb. 2017.[32] C. Pan, W. Xu, J. Wang, H. Ren, W. Zhang, N. Huang, and M. Chen, “Pricing-based distributed energy-efficient beamformingfor MISO interference channels,”
IEEE J. Sel. Areas Commun. , vol. 34, no. 4, pp. 710–722, Apr. 2016.[33] C. Pan, H. Zhu, N. J. Gomes, and J. Wang, “Joint precoding and RRH selection for user-centric green MIMO C-RAN,”
IEEETrans. Wireless Commun. , vol. 16, no. 5, pp. 2891–2906, May 2017.[34] E. Boshkovska, D. W. K. Ng, N. Zlatanov, and R. Schober, “Practical non-linear energy harvesting model and resourceallocation for SWIPT systems,”
IEEE Commun. Lett. , vol. 19, no. 12, pp. 2082–2085, 2015. [35] D. Mishra, G. C. Alexandropoulos, and S. De, “Energy sustainable IoT with individual QoS constraints through MISO SWIPTmulticasting,” IEEE Internet Things J. , vol. 5, no. 4, pp. 2856–2867, Aug. 2018.[36] K. Xiong, B. Wang, and K. R. Liu, “Rate-energy region of SWIPT for MIMO broadcasting under nonlinear energy harvestingmodel,”
IEEE Trans. Wireless Commun. , vol. 16, no. 8, pp. 5147–5161, Aug. 2017.[37] F. Wang, J. Xu, X. Wang, and S. Cui, “Joint offloading and computing optimization in wireless powered mobile-edge computingsystems,”
IEEE Trans. Wireless Commun. , vol. 17, no. 3, pp. 1784–1797, Mar. 2018.[38] J. Xu, L. Liu, and R. Zhang, “Multiuser MISO beamforming for simultaneous wireless information and power transfer,”
IEEETrans. Signal Process. , vol. 62, no. 18, pp. 4798–4810, Sep. 2014.[39] Q. Shi, M. Razaviyayn, Z.-Q. Luo, and C. He, “An iteratively weighted MMSE approach to distributed sum-utilitymaximization for a MIMO interfering broadcast channel,”
IEEE Trans. Signal Process. , vol. 59, no. 9, pp. 4331–4340,Sep. 2011.[40] C. Pan, H. Ren, M. Elkashlan, A. Nallanathan, and L. Hanzo, “The non-coherent ultra-dense C-RAN is capable ofoutperforming its coherent counterpart at a limited fronthaul capacity,”
IEEE J. Sel. Areas Commun. , vol. 36, no. 11, pp.2549–2560, Nov. 2018.[41] C. Pan, H. Ren, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Robust beamforming design for ultra-dense user-centricC-RAN in the face of realistic pilot contamination and limited feedback.” [Online]. Available: https://arxiv.org/abs/1804.03990[42] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” 2014.[43] S. Boyd and L. Vandenberghe,
Convex optimization . Cambridge university press, 2004.[44] C. Pan, W. Xu, W. Zhang, J. Wang, H. Ren, and M. Chen, “Weighted sum energy efficiency maximization in ad hoc networks,”
IEEE Wireless Commun. Lett. , vol. 4, no. 3, pp. 233–236, Jun. 2015.[45] X.-D. Zhang,
Matrix analysis and applications . Cambridge University Press, 2017.[46] J. Song, P. Babu, and D. P. Palomar, “Sequence design to minimize the weighted integrated and peak sidelobe levels,”