Specific Absorption Rate-Aware Beamforming in MISO Downlink SWIPT Systems

This paper investigates the optimal transmit beamforming design of simultaneous wireless information and power transfer (SWIPT) in the multiuser multiple-input-single-output (MISO) downlink with specific absorption rate (SAR) constraints. We consider the power splitting technique for SWIPT, where each receiver divides the received signal into two parts: one for information decoding and the other for energy harvesting with a practical non-linear rectification model. The problem of interest is to maximize as much as possible the received signal-to-interference-plus-noise ratio (SINR) and the energy harvested for all receivers, while satisfying the transmit power and the SAR constraints by optimizing the transmit beamforming at the transmitter and the power splitting ratios at different receivers. The optimal beamforming and power splitting solutions are obtained with the aid of semidefinite programming and bisection search. Low-complexity fixed beamforming and hybrid beamforming techniques are also studied. Furthermore, we study the effect of imperfect channel information and radiation matrices, and design robust beamforming to guarantee the worst-case performance. Simulation results demonstrate that our proposed algorithms can effectively deal with the radio exposure constraints and significantly outperform the conventional transmission scheme with power backoff.

Initial works in SWIPT assume a linear channel to transfer both information and power and investigate sophisticated transmission techniques and/or receiver architectures for SWIPT. The work in [4] characterizes the fundamental tradeoff between information and energy for a basic multipleaccess channel with information and EH receivers; it has been shown that a feedback channel is a critical mechanism to increase the information-energy capacity region. To enable practical techniques to convey information and power, the authors in [5] propose two practical receiver approaches, namely, a) "time switching (TS)", where the receiver switches between decoding information and harvesting energy, and b) "power splitting (PS)", where the receiver splits the received signal into two parts for decoding information and harvesting energy, respectively. Using the PS approach and multiple transmit antennas at the transmitter, the optimal multipleinput single-output (MISO) beamforming with both qualityof-service (QoS) and EH constraints was studied in [6] and [7] for an interference channel and a downlink MISO channel, respectively. To reduce the complexity of the SWIPT receivers, some studies overcome the use of RF-to-DC circuits and achieve information transfer by embedding information into the shape of the transmitted signal rather than modulating radio waveforms. The work in [8] introduces the integrated receiver, where information is embedded in the amplitude of energy signals, while decoding is performed by taking samples at the output of the rectification circuit. The authors in [9] also exploit the concept of pre-coded spatial modulation and convey information both at the energy pattern of the transmitted signal and the index of the received antenna. More recent works take into account the non-linearity of the rectification process and focus on the waveform design for SWIPT. By introducing a simple and tractable model of the diode nonlinearity (i.e., physics-based diode model), the work in [10] deals with the design of waveforms that maximizes the DC power at the output of the rectifier. This work is extended in [11] to convey both information and energy by transmitting a superposition of unmodulated and modulated multi-tone signals.
On the other hand, wireless technologies are characterized by terminals that are subject to strict regulations on the level of RF radiation that users of the terminals are exposed to. Since RF radiation has been proven harmful for humans and the environment, these regulations minimize the potential biological effects (e.g., tissue heating, metabolic changes in the brain, carcinogenic effects) caused by RF radiation [12], [13]. Two widely adopted regulations/measures on RF exposure are mainly considered i.e., the first one is the maximum permissible exposure (MPE) in [W/m 2 ], which is defined as the highest arXiv:1911.10556v1 [cs.IT] 24 Nov 2019 level of electromagnetic radiation (EMR) to which a user may be exposed without incurring an established adverse health effect. This is a particularly important issue in the applications of wireless power transfer (WPT) and various relevant studies have been carried out concerning EMR. For instance, the problem of scheduling the power chargers is investigated in [14] so that the charging utility for all rechargeable devices is maximized with a constraint on EMR. The works in [15] and [16] deal with the problem of maximizing the harvested energy and wireless charging tasks scheduling, respectively, when the transmitted signals guarantee a well-defined EMR constraint. The work [17] discusses the security issues in RF-based WPT systems and investigates techniques that ensure confidentiality, security and safety for different types of attacks. The second one is the specific absorption rate (SAR) that measures the absorbed power in a unit mass of human tissue by using units of Watt per kilogram [W/kg]. The Federal Communications Commission (FCC) enforces an SAR limitation of 1.6 W/kg averaged over one gram of tissue on partial body exposure [18]; other regulatory agencies (e.g., Comité Européen de NormalisationÉlectrotechnique in European Union) also adopt a similar limitation of 2 W/kg averaged over 10 grams of tissue on SAR measurements [19]. For short distances (e.g., the distance between the transmitter and the nearest user is less than 20 cm), the SAR measure dominates the RF exposure and therefore becomes a more critical factor than the MPE for the design of efficient uplink communications. SAR is obviously different from MPE in that it is a different measure of RF exposure and applies to short distances. In addition, it needs to account for various exposure constraints on the whole body, partial body, hands, wrists, feet, ankles, etc., with various measurement limitations according to FCC regulations [18], therefore multiple SAR constraints are needed even for a single transmit device. For instance, the iPhone 6 Model A1549 has a whole body SAR of 1.14 W/kg and head SAR of 1.08 W/kg [20].
The integration of RF exposure constraints, and specifically the SAR constraint, in the design of wireless communication systems is a timely research area with significant impact; however, few works in the literature take into account SAR regulations. In conventional single-antenna wireless communication systems, the SAR exposure limitation simply introduces an additional transmit power constraint, and can be easily guaranteed by reducing the transmit power below a specific threshold. However, it becomes more involved when there are multiple transmit antennas and signals can be beamformed towards certain directions. In this case, the exposure limitation is coupled with the multi-antenna channel vector/matrix, so it needs to be incorporated into the transmitter optimization. The measurements and simulations carried out in [21] demonstrate SAR as a function of the phase different between two transmit antennas. The SAR reduction and modeling in multiantenna systems has been studied in [22], [23]. The SAR constraints are incorporated into the transmit signal design for a multiple-input multiple-output (MIMO) uplink channel in [24], where the quadratic model for the SAR measurements is first proposed. In [25], a SAR code is proposed to handle a SAR-constrained channel, which is shown to have significant improvement over the conventional Alamouti space-time code under SAR limitations. It is also revealed that the SAR measurement is a function of the quadratic form of the transmitted signal with the SAR matrix. The SAR-aware beamforming and transmit signal covariance optimization methods are first presented in [26] and [27]. Capacity analysis with multiple SAR constraints on single-user MIMO systems is intensively examined in [28]. Sum-rate analysis for a multi-user MIMO system with SAR constraints is performed in [29] with both perfect and statistical channel state information (CSI). Form an information theoretical perspective, SAR-constrained multiantenna transmit covariance optimization can be considered as the classical MIMO channel capacity optimization problem subject to generalized linear transmit covariance constraints studied in [30].
The integration of SAR exposure constraints in the design of wireless power transfer or SWIPT systems is an unexplored research area [31]. Although SWIPT corresponds to a controlled transmission of RF radiation to communicate and energize, and may significantly contribute to the electromagnetic pollution (electrosmog), no work in the literature discusses the integration of SAR in SWIPT. For example, one of the main application areas of SWIPT is for medical devices in wireless body area networks, where an access point will support the communication connectivity and the power sustainability of a short-range sensor network in, on, or around the human body [32], [33], [34]. Since the access-point may be placed within 20 cm of the body, SAR measurements are mandatory and should be taken into account in the SWIPT design. To the best of our knowledge, the impact of the SAR on SWIPT has not been studied before.
In this paper, we focus on the PS approach and study the beamforming design in a multiuser MISO downlink channel, where the receivers are characterized by both communications QoS and EH constraints with additional SAR limitations; the nonlinearity of the rectification circuit is also taken into account. The contributions of this paper are summarized as follows: 1) We introduce and formulate the SAR-aware beamforming optimization problem with simultaneous QoS and EH requirements. 2) We derive the optimal beamforming solution in the general case by leveraging semidefinite programming (SDP) together with rank relaxation. More importantly, we prove that the proposed approach always gives rank-1 solutions, such that the exactly optimal beamforming solutions can be derived. In the special case of a single receiver with a single SAR constraint, we propose a fast algorithm to achieve the optimal solution. 3) We develop low-complexity suboptimal solutions using both fixed beamforming schemes based on the criterion of zero-forcing (ZF) or regularized ZF (RZF), and a hybrid beamforming scheme that combines the ZF and maximum ratio transmission (MRT) to achieve a better performance-complexity tradeoff. 4) Furthermore, when the CSI and/or the SAR matrices are imperfectly known at the transmitter, we devise a robust beamforming scheme suitable for imperfect CSI and SAR matrices with bounded errors, to guarantee the QoS and EH requirements, as well as the SAR constraints in the worse case. In addition, we have made the following assumptions: • The input signal is circularly symmetric complex Gaussian (CSCG) distributed. This is in general sub-optimal and the input distribution under nonlinearity can be optimized as well, as studied in [10], [35], [36]. • We adopt a simplified nonlinear EH model that depends only on the power of the received signal. Note that a more accurate EH model is a function of the received signal other than just its power, as shown and experimentally demonstrated in the literature [2], [37]. Notation: All boldface letters indicate vectors (lower case) or matrices (upper case). The superscripts (·) T , (·) † , (·) −1 , (·) ‡ , diag(·) denote the transpose, the conjugate transpose, the matrix inverse, the Moore-Penrose pseudo-inverse and the diagonal elements respectively. A circularly symmetric complex Gaussian random variable z with mean µ and variance σ 2 is represented as z ∼ CN µ, σ 2 . The identity matrix of size M , and the zero matrix of size m × n, are denoted by I M and 0 m×n , respectively. z and z F denote the L 2 and Fibonacci norms of a complex vector z, |z| denotes the magnitude of a complex variable z, trace(A) denotes the trace of a matrix A, while A 0 indicates that matrix A is positive semidefinite. i √ −1 denotes the imaginary unit. This paper is organized as follows: Section II sets up the system model and introduces the optimization problem considered. In Section III, we investigate the joint optimization of the optimal beamforming design and power splitting ratios. Section IV discusses solutions for the fixed and hybrid beamforming schemes, while Section V develops the robust solution under imperfect estimation of the channels and the SAR matrices. Simulation results are presented in Section VI, followed by conclusions in Section VII.

II. SYSTEM MODEL AND PROBLEM FORMULATION
We assume a MISO downlink channel consisting of an N tantenna transmitter (e.g., a base station (BS)) and K singleantenna receivers that employ single-user detection as illustrated in Fig. 1. This ensures simple receiver processing due to size limitation or limited processing power. We assume that the BS transmits with a total power P T and let s k be its transmitted data symbol to receiver k, which is Gaussian distributed with zero mean and unit variance, i.e., E{ s k 2 } = 1. The transmitted data symbol s k is mapped onto the antenna array elements by the beamforming/precoding vector w k ∈ C Nt×1 . Other properties of the signal (e.g., modulation, waveform and input distribution) can also be optimized to improve the efficiency of RF-DC conversion [10], [35], so the CSCG distribution of the input signal is in general sub-optimal. However, the joint design of signal and beamforming is a challenging problem and out of the scope of this paper. All wireless links exhibit independent fading and additive white gaussian noise (AWGN) with zero mean and certain variances. The fading is assumed to be frequency non-selective block fading. This means that the fading coefficients in the vector channel h k ∈ C Nt×1 remain  constant during one slot, but may change independently from one slot to another. The mean and variance of the channel coefficients capture large-scale degradation effects such as path-loss and shadowing. The received baseband signal at receiver k can be expressed as where n k denotes the AWGN component with zero mean and variance N 0 . Therefore, the received power at receiver k is equal to The receivers have RF-EH capabilities and therefore can harvest energy from the received RF signal based on the power splitting technique [5], which is a mature SWIPT technique and does not require strict time synchronization between information and power transfer 1 . With this approach, each receiver splits its received signal into two parts: a) one part is converted to a baseband signal for further signal processing and data detection, and b) the other part is driven to the required circuits for conversion to DC voltage and energy storage. Let ρ k ∈ (0, 1) denote the power splitting parameter for the k-th receiver; this means that 100ρ k % of the received power is used for data detection, while the remaining amount is the input to the RF-EH circuitry. More specifically, after reception of the RF signal at the receiver, a power splitter divides the power P r k into two parts according to ρ k , so that ρ k P r k is directed towards the decoding unit and (1 − ρ k )P r k towards the EH unit. During the baseband conversion, additional circuit noise, v k , is present due to phase offsets and non-linearities which is modeled as AWGN with zero mean and variance N C .
The signal-to-interference-plus-noise ratio (SINR) metric characterizing the data detection process at the k-th receiver is given by On the other hand, the total power that can be harvested is is a non-linear parametric EH function. Note that in general F (·) should be a function y k rather than just the power of the energy signal. In this paper we adopt a simplified model to highlight the dependency on the power of the energy signal only. More details about F (·) will be discussed later.
As discussed in the introduction, wireless communication devices are subject to SAR limitations. Previous reported results such as [25] have shown that the pointwise SAR value with multiple transmit antennas can be modeled as a quadratic form of the transmitted signal, and the SAR matrix fully describes the SAR measurements dependence on the transmitted signals; entries of the SAR matrix have units of W/kg. Because the SAR measurements are always real positive numbers, the SAR matrices are positive-definite conjugatesymmetric matrices.
Since SAR is a quantity averaged over the transmit signals, we model the l-th SAR constraint with a time-averaged quadratic constraint given by where A l 0 is the l-th SAR matrix and P l is the l-th SAR limit.

A. Problem Formulation
We study the following problem of maximizing the ratios of the received SINR and EH over the target requirements, i.e., , · · · , P S K λ K } subject to the SAR and total power constraints, whereγ k ,λ k are the SINR and the EH requirements, respectively. The choice of the objective function will balance the received SINR and the EH between users. To make the problem more tractable, we introduce an auxiliary variable t, and formulate the optimization problem as follows where P T is the maximum total transmit power, and L is the number of SAR constraints.
In this paper, we adopt a non-linear parametric EH model, so the output DC power at the k-th receiver can be represented by the nonlinear function F (x), where x is the input RF power. The nonlinear function can take many forms to capture the relationship between the input and output power at the energy receiver, such as the sigmoid function [40], [41], the linear fraction [42], the piece-wise linear function [43], the second order polynomial model [44], and a heuristic expression by curve fitting in [45]. We use F to denote a general nonlinear energy harvesting function; our analysis is general and independent of the specific function F . In general, the nonlinear EH function is monotonically increasing, therefore we can find the inverse mapping F −1 (·), and the EH constraint (8) can be rewritten as In general, it is difficult to solve the above problem P1, because both SINR and EH constraints (6) and (7) are nonconvex, and we also have additional multiple SAR constraints. Instead, we solve the following power minimization problem P2 below P2: min Clearly, the problem P2 is closely related to the problem P1. For instance, if P1 is solved and the optimal t * is achieved, then the same beamforming and power splitting solutions are also optimal for P2 to achieve the same SINR and EH, and the optimal minimum power will be P T . Because P1 is a quasiconvex problem in t, once P2 is solved, P1 can be solved via a bisection search algorithm over t as follows (note that t in P1 is embedded in γ k , λ k of P2): Algorithm 1 to Solve P1: 1) Set the upper and lower bounds of t as t U and t L . Repeat the following steps until convergence. 2) Calculate t = t U +t L 2 , and solve P2 with γ k = tγ k and Therefore, in the rest of the paper we will focus on solving the problem P2.
III. THE OPTIMAL SCHEME In this section, we investigate the optimal beamforming and power splitting solutions to P2, and we begin with developing a fast solution for the special case of P2 referring to a singlereceiver system with a single SAR constraint, i.e., K = L = 1.
A. Fast Solution to the Special Case of K = L = 1 When we have only a single user and a single SAR constraint, the problem P1 reduces to the problem below, where the receiver and SAR indices are dropped for the sake of simplicity, We find the following Lemma is useful to derive the optimal solution to P1. Lemma 1: Both the SINR and the EH constraints in (12) and (13) must be satisfied with equality at the optimum of P3.
Proof: It can be proved by contradiction. If (12) is not tight, ρ can be reduced without affecting (13); if (13) is not tight, ρ can be increased without affecting (12). This completes the proof.
Base on Lemma 1, we know that the optimal solution to P1 falls into the following two cases: • Case I: Only the SINR and the EH constraints are satisfied with equality.
In this case, we first solve the optimalρ. From the equality constraint (12), we can get Substituting (15) into (13), we can have the following equation, i.e., which leads to a quadratic equation of ρ given by Its solution is written as where It is easy to verify thatρ is between 0 and 1, because the left hand side of the quadratic equation (17) is strictly negative when ρ = 1 and strictly positive when ρ = 0. Note thatρ is the optimal solution regardless of w as long as P3 is feasible. Onceρ is found, the optimal w can be solved by satisfying the equality (15) and minimizing the objective simultaneously, which is given bȳ After solvingw, we can substitute it into (14) and check whether the SAR constraint is satisfied: if so,w is the optimal beamforming solution; otherwise, we need to consider the next case.
• Case II: All three SINR, EH and SAR constrains are satisfied with equality.
In this case, we can simplify P3 to the following problem whereρ is given by (18) P4: min We keep the inequality constraints for the convenience of the derivation. In general, it is difficult to find the optimal solution to P4 directly. We first study the analytical beamforming structure that facilitates the development of an efficient solution.
The Lagrangian of P4 is given by where a ≥ 0 and b ≥ 0 are dual variables. The dual problem is then formulated as: and the optimal w admits the following form Let the eigenvalue decomposition of A be A = UDU † , where U is a unitary matrix and D is a diagonal matrix with elements being the eigenvalues of A. Then we have (I + bA) −1 = U(I + bD)U † , and and wherec is defined in P4. From (27), we can solve Now we proceed to consider the constraint (21) and define the function below: where (a) is due to (26), A = UDU † is used to obtain (b), and we use (28) to get (c). Clearly, the parameter b that satisfies f (b) = 0 uniquely determines the optimal solution. We find the following theorem useful to develop an efficient numerical solution to solve f (b) = 0. (29) is a decreasing function in b.
Proof: Suppose D = Diag(d 1 , · · · , d Nt ) and |h| = [h 1 , · · · , h Nt ]. Then after some algebraic manipulation, From (30), it can be seen that f (b) is a decreasing function in b, and this completes the proof. Theorem 1 can be verified by the simulation result in Fig.  2.
Based on Theorem 1, we propose the following bisection search algorithm to solve b. Algorithm 2 to Solve f (b) = 0.
1) Initialize b min ≥ 0 and b max > 0. Repeat the following steps until convergence.
Once the optimal b is found, the optimal optimal beamforming solution and the optimal a can therefore be solved via (26) and (28), respectively.

B. The Optimal Solution using SDP in the General Case
The general case of the problem P2 is nonconvex and difficult to solve. In this section, we develop an efficient SDP based algorithm that jointly optimizes the beamforming vectors and power splitting parameters. To tackle P2, we first introduce new matrix variables W k = w k w † k , and then P2 can be reformulated as The problem P5 is convex because it is linear in all {W k } , and both terms 1 ρ k and 1 1−ρ k are convex in ρ k > 0. It can be efficiently solved using numerical software package such as CVX [52]. Once P5 is optimally solved, if the resulting solutions {W k } are all rank-1, then they are the exactly optimal solutions; otherwise, the solutions only provide a lower bound for the minimum required transmit power.
It has been proved in [7] that without the SAR constraints (32), the above SDP with rank relaxation will produce rank-1 solutions {W k } which means they are also optimal to the problem P2. However, whether the SDP with rank relaxation can generate the optimal solution highly depends on the problem structure. With the additional SAR constraints, it is unknown whether this property remains true. In the following theorem, we prove that this is still the case.
Theorem 2: The optimal solution to P5 satisfies rank(W k ) = 1, ∀k, i.e., the SDP relaxation is tight, and the optimal solution to the problem P2 can be recovered from {W k } via eigenvalue decomposition. The proof is given in Appendix A.

IV. SUB-OPTIMAL BEAMFORMING SCHEMES
Although the algorithm based on SDP in Section III ensures the optimal beamforming and power splitting solutions, its complexity is high. In this section, we focus on lowcomplexity suboptimal solutions including both the fixed heuristic beamforming schemes and the hybrid beamforming scheme that considers a linear combination of fixed beamforming strategies.

A. Solutions using Fixed Beamforming
Suppose the fixed beamforming vector is given by w k = √ p k w f k , w f k = 1, where p k is the power for the k-th receiver. Let G k,j |h † k w f j | 2 denote the link gain between the BS and the k-th receiver, and F k,l = w f k † A l w f k denote the radiation channel gain due to the transmission intended for the k-th receiver.
Using the fixed beamforming and rearranging the terms in the problem P2, we have P6 : min The optimization problem P6 is convex because it is comprised of a linear objective function and convex constraints. Constraint (33b) is convex because the term 1 ρ k is convex for ρ k > 0 and the other terms are linear, (33c) is a restricted hyperbolic constraint and the term 1 1−ρ k is also convex. Note that by solving the problem P6, we obtain optimal values for the splitting parameters, as well as an optimal power allocation for any given fixed beamforming vectors.
Commonly used fixed beamforming vectors include the MRT, ZF and RZF criteria [51], defined respectively below: where we have defined H = [h 1 , . . . , h K ], and H k = [h 1 , . . . , h k−1 , h k+1 , . . . , h K ]. The MRT beamforming and the ZF beamforming aim to enhance the received signal strength and remove interference, respectively, while the RZF beamforming provides a tradeoff between them. Note that in the RZF beamforming, the SAR matrices are included to incorporate the radiation constraints. As a special case, for the ZF beamforming, the optimization problem can be further simplified because G i,j = 0, for i = j. Hence the optimization problem P6 simplifies into the following formulation P7: min The problem P7 is also convex but is much easier to solve than P6 because of the simplified constraints. Unlike the previous works [6], [7] where closed-form solutions to P7 were derived, the SAR constraint does not permit a closed-form solution, therefore we use CVX [52] to solve both P6 and P7.

B. Hybrid Beamforming
In this subsection, we introduce the hybrid beamforming scheme to provide a tradeoff between enhancing the received SINR and energy harvested, dealing with interference, and accounting for the SAR constraints, which admits the following expression where x k and y k are combining coefficients. Then The transmit power intended for the k-th receiver is where w ZF ). The l-th SAR power can be reformulated as where we have defined A k,l w ZF Then the power minimization problem P2 becomes P8: min The problem P8 is not convex because of the nonconvex term √ x k y k introduced by the EH and the SAR constraints, which makes the problem more challenging to solve. To deal with this difficulty, we introduce s k = √ x k y k and relax it to s 2 k ≤ x k y k or [2s k ; x k − y k ] ≤ x k + y k , then we have the following second-order cone programming (SOCP) problem formulation P9: min The problem P9 is convex and can be optimally solved. We observe that in most cases, s k = √ x k y k holds true which means it is also the optimal solution to P8; otherwise, we will employ (38) to form the hybrid beamforming, and then use the fixed beamforming scheme in Section IV.A to find the optimal power vector and power splitting by solving P6. Complexity Analysis: The complexity of the optimally solution to the general problem P5 is dominated by the SDP constraints, and according to [53, 6.6.3], the complexity of the interior-point algorithm for solving P5 is While the for solving the fixed beamforming optimization problem P6 and the hybrid beamforming problem P9, their complexities are dominated by the linear and SOCP constraints, and can be expressed as O √ 3K + L(3K 3 + LK 2 ) [53, 6.6.1] and O √ 3K + L(28K 3 + 4LK 2 + 4K 2 + KL) [53, 6.6.2], respectively. It can be clearly seen that the suboptimal beamforming schemes reduces the complexity significantly when compared to the optimal beamforming.

V. ROBUST BEAMFORMING SCHEMES
In this section, we study the robust beamforming design, when the channel information and the SAR matrices are imperfectly known due to estimation and measurement errors. Without a robust solution, the SINR constraints and the EH constraints may not be satisfied; in addition, the SAR constraints may be violated.
We model the k-th receiver's actual channel as whereĥ k denotes the CSI estimate known to the BS. ∆h k represents the CSI uncertainty that belongs to the set below where σ 2 k denotes the error bound.
Similarly, the l-th SAR matrix is modelled as whereÂ l is know to the BS. ∆A l is the SAR uncertainty within the set below where τ l denotes the error bound. We assume that the BS has no knowledge about the error channel vectors or the error SAR matrices, except for their error bounds σ 2 k and τ l . Thus we take a worst-case approach for the beamforming design to guarantee the resulting solution is robust to all possible channel and SAR uncertainties within the given error sets. The specific robust design problem is to minimize the overall transmit power P T for ensuring the receivers' worst-case individual SINR, EH and SAR constraints, i.e.,

P10:
min We first ignore the rank constraint (50) and define Q k = trace((Â l + ∆A l )( K j=1 W k )) ≤ P l , ∀l, ∆A l ∈ A l (54) The constraint in (52) is in general difficult to handle because the robust beamforming needs to satisfy infinite number of channel realizations defined in the set H k . Fortunately, thanks to S-Procedure, it can be equivalently written as the following SDP constraint [54] Following the same principle, the constraint in (53) can be equivalently written as The worst-case SAR constraint (54) introduces further difficulty, and next we convert it to an equivalent convex constraint, and the result is summarized in the theorem below.
Theorem 3: The worst-case SAR constraint (54) is equivalent to the constraint The proof is given in Appendix B.
With the robust constraints (55), (56) and (58), we can obtain the equivalent robust beamforming optimization problem as follows P12: min where we have ignored the rank constraint on W k . However, P12 is still a nonconvex problem because of the nonlinear terms about ρ k in the constraints (60) and (61). To tackle this new challenge, we introduce the auxiliary variables {m k , n k } and use the following formulation to convert it to a convex problem below P13: min Next we prove that the constraint (63) is tight at the optimum, therefore P13 is equivalent to P12.
Theorem 4: There exists optimal solutions {m * k , n * k , ρ * k } to P13 such that m * k = 1 ρ * k and n k = 1 1−ρ * k , ∀k. Proof: This can be proved by noting that given a solution {m * k , n * k , ρ * k }, if any constraint (63) is not tight, we can always reduce m * k and/or n * k to make inequality constraints to be equalities without affecting the feasibility of other constraints. Therefore there always exists {m * k , n * k , ρ * k } such that the constraint (63) is tight. This completes the proof.
Theorem 4 shows that introducing the auxiliary variables {m k , n k } in P13 to deal with the nonconvex constraints in P12 incurs no loss of optimality. However, since we have used SDP with relaxation in P13, only when the relaxation is tight, i.e., rand(W k ) = 1, ∀k, the optimal robust beamforming solution w k can be recovered from W k ; otherwise, the obtained solution provides a lower bound of the required transmit power. A feasible beamforming solution can be obtained by the standard randomization technique [55].

VI. SIMULATION RESULTS
Simulations are carried out to evaluate the performance of the proposed algorithms. We consider a MISO downlink consisting of K receivers randomly located around the BS with distance l k and direction ζ k drawn from the uniform distribution, l k ∼ U (1, 5)m and ζ i ∼ U (−π, π). Each receiver can harvest energy at frequency f = 915 MHz and it is assumed that the antenna gains at the BS and receivers are 8dBi and 3dBi, respectively. The path attenuation of receiver k, L k , is obtained using the Friis equation with reference distance 1m and path loss coefficient 2.5. Because of the short distance between the BS and the receivers and dominance of the lineof-sight (LOS) signal, the Rician fading is used to model the channel. Hence, h k is composed of the LOS signal, h LOS k and the non-LOS signal h N LOS k according to the expression below [7] where R = 5 dB is the Rician factor. For the LOS signal the far-field uniform linear antenna array model with λ/2 distance between antenna elements is considered [50] which implies 1π sin ζ k ) , ..., e −j((N −1)π sin ζ k ) ] T . Rayleigh fading is adopted for the NLOS signal, h N LOS k ∈ C Nt×1 which means that each of its elements is a circularly symmetric complex Gaussian (CSCG) random variable with zero mean and variance L k .
Unless otherwise specified, it is further assumed a BS with four antennas serving four receivers, i.e., K = N t = 4, N 0 = −70 dBm and N C = −50 dBm, while the SINR and EH thresholds are the same for all receivers, i.e.Γ k =Γ = 10 dB,λ k =λ = −15 dBm, ∀k. we assume that the total power constraint is P T = 2 W, and the SAR power constraint is P l = P = 1.6 W/kg, ∀l. There is no guarantee that the diode (rectification circuit) operates always in the transition regime, and the system could operate in the close to saturation regime, so it is important to capture saturation effects. We use the nonlinear energy harvesting model below proposed in [42] F k (x) =ā x +b and the fitted parameters of the proposed model areā = 2.463,b = 1.635, c = 0.826 using the data in [46]. (65) belongs to the category of models which are based on realworld measurements by adjusting the parameters of a nonlinear function through curve fitting tools. It can model both saturation and non-saturation regimes, and has been widely used in the literature for the design of WPT/SWIPT systems such as [47], [48], [49]. Note that in general F (·) should be a function of the received signal rather than just its power, so (65) is a simplified model based on specific measurements under specific conditions and with continuous wave input signals. The SAR matrix is given by for four antennas at the BS.
We use the minimum achievable SINR and EH ratio t as the main performance evaluation criterion. The proposed optimal solutions, fixed beamforming solutions including ZF and RZF beamforming schemes as well as the hybrid beamforming scheme will be compared. In addition, the conventional backoff scheme is used as another performance benchmark described as follows. Firstly, it solves P1 without the SAR constraint, i.e., the problem below Suppose its solution is w. Define δ l K k=1 w † k A l w k P l , then the backoff solution is given by w = w min(1,max l δ l ) to satisfy the SAR constraints.
We first evaluate the performance and complexity of the fast algorithm for a single receiver system with a single SAR constraint (K = L = 1) as the SAR limit varies in Fig. 3. From Fig. 3 (a), we can see that the minimum achievable SINR and EH ratio of the proposed fast algorithm is identical to that of the SDP approach for systems with general number of receivers, which verifies its optimality. Fig. 3 (b) shows that the computation time of the proposed fast algorithm is more than three orders of magnitude lower than the SDP approach. These results clearly demonstrate the advantage of the proposed fast algorithm.
Figs. 4, 5 and 6 depict the minimum achievable SINR and EH ratio by the different investigated algorithms by varying P , λ and P T , respectively. Fig. 4 (a) depicts that the proposed optimal solution can achieve substantial gain over the traditional backoff solution, especially when the SAR power constraint is low. As the SAR power constraint increases, the performance of the optimal solution is close to that without the SAR power constraint. The hybrid beamforming solution is superior to the two fixed beamforming solutions, while ZF beamforming performs much better than the RZF beamforming. Fig. 4 (b) shows the feasibility of various schemes, which follows the similar trend as the results in Fig. 4 (a). In Fig. 4 (c), we plot the actual transmit power consumed. It can be seen that without the SAR constraint, the total power budget of 2W is always used. The ZF beamforming solution uses similar power as the backoff solution, while the hybrid beamforming scheme uses much less power, thus provides a higher power efficiency.
In Fig. 5, we plot the minimum achievable SINR and EH ratio against the EH constraint. It is seen that the proposed optimal solution outperforms other schemes significantly. As expected, as the required EH increases, the performance of all schemes degrades, and eventually converges to the similar result. The hybrid beamforming scheme achieves similar performance as the backoff scheme, but with much reduced complexity.
In Fig. 6, we examine the impact of the total transmit power on the minimum achievable SINR and EH ratio. It is observed that the performance of all schemes is increasing or nondecreasing as the total transmit power increases and saturates  when the transmit power is greater than 34 dBm, except for the backoff solution. The reason is that although the transmit power increases, the SAR constraint becomes the bottleneck, so the increased power cannot be utilized and transferred to improved performance. The backoff solution is penalized most in the high power regime, therefore its performance deteriorates quickly as the transmit power goes above 31 dBm.
Next, we investigate the impact of the channel and SAR matrix estimation error. Fig. 7 depicts the minimum achievable SINR and EH ratio of different schemes when the estimation error σ 2 k = τ l , ∀k, l varies. The performance degradation is clearly seen as the estimation error grows. The performance of the proposed robust solution is still satisfactory when the estimation error is small, and it can achieve much higher SINR and EH than the non-robust solution that ignores the estimation error. Fig. 8 (a) and (b) further depict the cumulative distribution functions (CDFs) of the output SINR and EH, respectively, when there are estimation errors given by σ 2 k = 5 × 10 −8 , ∀k and τ l = 7 × 10 −8 , ∀l. The vertical lines show the target SINR of 10 dB and the target EH of -15 dBm. It is clearly seen that our proposed robust solution can guarantee the worst-case SINR and EH constraints are met even in the presence of estimation error. However, for the nonrobust solution, both constraints are violated with probabilities higher than 90%.

VII. CONCLUSIONS
In this paper, we have studied the optimization of SARconstrained multiuser transmit beamforming of a SWIPT system. In the general case with perfect information, we have shown that the optimal beamforming and power splitting solutions can be obtained via semidefinite programming and bisection search; while a much more efficient solution can be found for the special single-receiver single-SAR case. We further designed low-complexity suboptimal solutions including the fixed beamforming and hybrid beamforming schemes. In addition, we proposed robust beamforming solutions to   deal with the imperfect channel and SAR matrix information, while guaranteeing the required SINR, EH constraints and the maximum SAR is below the given threshold. Our simulation and analysis have shown significant performance improvement of the proposed SAR-aware optimal solution over the conventional transmission scheme with simple power backoff. Future works include using large intelligence surface to further reduce the energy consumption of the transmitter while satisfying the SAR constraints.

Minimum Achievable SINR and EH Ratio
Optimal Solution with perfect information Robust Solution Non-robust Solution Figure 7. The minimum achievable SINR and EH ratio vs the estimation error. Proof: The partial Lagrangian of the problem P5 is L({W k , ρ k , α k , β k , ν l }) where {α k , β k , ν l } are dual variables, and we have defined ν l A l . 13 So the dual problem is max ν,α,β≥0 ν l A l 0, ∀k.
Next we prove that the matrix I + K j=1 (α j − β j (1 − ρ j ))h j h † j + L l=1 ν l A l is full rank by contradiction. If it is not full-rank, suppose there exits a non-zero vector x that satisfies Therefore it holds true that It follows that which contradicts the assumption that x † (I + K j=1 (α j − β j (1 − ρ j ))h j h † j + L l=1 ν l A l )x = 0. Therefore the matrix I + K j=1 (α j − β j (1 − ρ j ))h j h † j + L l=1 ν l A l must be full rank, and the rank of X k is at least N t − 1.
One KarushKuhnTucker (KKT) condition of the problem P5 is that trace(W k X k ) = 0, so the rank of W k is at most 1. This completes the proof.

APPENDIX B. PROOF OF THEOREM 3
Proof: We consider the following optimization problem: where X has the meaning of the worst SAR error matrix. Its Lagrangian is given by where u is the dual variable. Setting its first-order derivative to be zero leads to: 2uX −W = 0.
(75) Therefore X should be a scaled version of theW. It is easy to see that the norm constraint in (73) should be tight, so we can obtain the worst-case X as X * =W W F τ , and max X trace(WX) = τ W F . Substituting it into (54) leads to the constraint (58). This completes the proof.