Design and Analysis of SWIPT with Safety Constraints

Simultaneous wireless information and power transfer (SWIPT) has long been proposed as a key solution for charging and communicating with low-cost and low-power devices. However, the employment of radio frequency (RF) signals for information/power transfer needs to comply with international health and safety regulations. In this paper, we provide a complete framework for the design and analysis of far-field SWIPT under safety constraints. In particular, we deal with two RF exposure regulations, namely, the specific absorption rate (SAR) and the maximum permissible exposure (MPE). The state-of-the-art regarding SAR and MPE is outlined together with a description as to how these can be modeled in the context of communication networks. We propose a deep learning approach for the design of robust beamforming subject to specific information, energy harvesting and SAR constraints. Furthermore, we present a thorough analytical study for the performance of large-scale SWIPT systems, in terms of information and energy coverage under MPE constraints. This work provides insights with regards to the optimal SWIPT design as well as the potentials from the proper development of SWIPT systems under health and safety restrictions.


I. INTRODUCTION
Wireless technologies are an important part of modern society, with Cisco expecting that global mobile subscribers will reach 5.7 billion by 2023, which will correspond 71 percent of the global population [1]. Traditionally, the focus of wireless communications is mainly on how to improve the efficiency of information transfer. Nonetheless, with the development of wireless systems employing a massive number of devices, such as sensors and actuators, recently this focus has also been shifted towards energy sustainability. In particular, the fact that radio frequency (RF) signals can also convey energy apart from information, the concept of wireless power transfer (WPT) and, in particular, of simultaneous wireless information and power transfer (SWIPT) is considered as a very promising and enabling technology for the realization of such future systems [2]. The key idea of SWIPT is to exploit the received RF signals in order to extract not only information but energy as well. The extraction can be done by separating the energy harvesting (EH) and information decoding operations either in time, in power or in space [3]. In contrast to conventional EH techniques (e.g. from renewable sources), EH with SWIPT can be a dedicated, continuous, controllable and on-demand process. The EH in this case, is achieved through the employment of a rectifying antenna (rectenna) that converts the received RF signal to direct current (DC) [2], [3].
In light of the wide and increasing usage of wireless devices, there is a growing concern with regards to the RF radiation brought about by multiple and concurrent wireless transmissions. Some studies have also revealed the potential biological hazard in relation to the RF radiation, including metabolic changes in brain and carcinogenic effects [4], [5]. As such, international health and safety regulations have been put in place in order to regulate and limit the level of RF exposure to humans [6]. Two widely adopted regulations/measures on RF exposure are the so-called maximum permissible exposure (MPE) and specific absorption rate (SAR). The MPE or power density (measured in W/m 2 ), defines the highest level of electromagnetic radiation (EMR) in a specific area that will not incur any health/biological effect. On the other hand, SAR (measured in W/kg) is a localized metric and defines the maximum level of absorbed power in a unit mass of human tissue.
Due to the vital impact on applications with WPT, MPE is concerned by some relevant studies [7]- [10]. The problem of scheduling the power chargers is investigated in [7], where the charging utility for all rechargeable devices is maximized with a constraint on EMR. The works in [8] and [9] deal with the problem of maximizing the harvested energy and wireless charging tasks scheduling, respectively, when transmitted signals guarantee a well-defined EMR constraint. The work in [10] builds an empirical probabilistic wireless charging and EMR model and then formulates an optimization problem to maximize the charging utility of all EH devices with consideration on the EMR constraint. For short distances (i.e., distances less than 20 cm), the SAR measure dominates the RF exposure. Therefore, in this case, SAR becomes a more critical factor than the MPE for the design of efficient communications. Some regulatory agencies have established limitations on the body SAR exposure. For instance, the Federal Communications Commission (FCC) enforces a SAR limitation of 1.6 W/kg averaged over one gram of tissue on the partial body exposure [11], and the Comité Européen de Normalisation Électrotechnique adopts a similar limitation of 2 W/kg averaged over 10 grams of tissue on the SAR measurements [12]. In addition, SAR 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 [11]. Thus, multiple SAR constraints are needed even for a single transmit device. For instance, the iPhone 12 Pro Model A2407 has a whole body SAR of 1.18 W/kg and a head SAR of 1.14 W/kg [13].
Although the subject of an RF exposure constraint (i.e., the SAR constraint) has an important impact on the design of wireless communication systems, few existing works in the literature have considered SAR regulations. The SAR exposure limitation can be easily guaranteed in single-antenna systems by introducing an additional transmit power constraint (i.e., setting the transmit power below a required threshold). In the meantime, the exploitation of multiple antennas provides significant benefits in improving the throughput of a wireless communication system, but it also poses potential challenges during the design due to RF radiation restrictions. For example, by exploiting advanced signal processing techniques such as beamforming in multi-antenna systems, the pattern of the RF signals is manipulated in such a way as to increase the performance, while it also makes the associated SAR analysis more complicated. Measurements and simulations are carried out in [14] and demonstrate that SAR is a function of the phase difference between two transmit antennas. The SAR reduction and modeling in multi-antenna systems is studied in [15], [16]. The SAR constraints are integrated into the transmit signal design for a multiple-input multiple-output (MIMO) uplink channel in [17], where the quadratic model for the SAR measurements is first proposed. In [18], a SAR code is proposed to improve the conventional Alamouti space-time code under SAR constraints. 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 presented in [19] and [20]. Capacity analysis with multiple SAR constraints on single-user MIMO systems is intensively examined in [21]. Sum-rate analysis for a multi-user MIMO system with SAR constraints is performed in [22], with both perfect and statistical channel state information (CSI). From an information theory perspective, SAR-constrained multi-antenna transmit covariance optimization can be seen as the classical MIMO channel capacity optimization problem subject to generalized linear transmit covariance constraints [23].
The consideration of SAR and MPE constraints for the design of WPT or SWIPT systems is an under-explored research area [2]. Even though SWIPT corresponds to a controlled transmission of RF radiation to communicate as well as energize, it may significantly contribute to the electromagnetic pollution (electrosmog). However, few works in the literature discuss the integration of SAR and MPE with WPT [7]- [10] or SWIPT [24]. As such, this paper provides a complete framework for the study of SWIPT, under both SAR and MPE constraints. We first present the mathematical modeling of the SWIPT technology as well as the modeling aspects of the two safety metrics. Then, we describe methodologies for the design of simple but also complex large-scale SWIPT systems under safety constraints; the methodologies are general and can be adapted to any communication scenario. Specifically, the contributions are as follows: • We first introduce and formulate the beamforming optimization problem to maximize the harvested power in a multiple-input single-output (MISO) downlink system, subject to SAR constraints and quality-of-service requirements. We derive the beamforming solution by leveraging semidefinite programming and rank relaxation, and prove that this method always achieves the optimal solution.
• Next, a robust beamforming design is proposed subject to SAR constraints at the receivers. In practical systems, the CSI is usually measured or estimated, while there are many factors contributing to errors, e.g. quantization errors [25]. In such cases, the constraints such as the user end performance, measured by the signal-to-interferenceplus-noise ratio (SINR) are characterized in a statistical instead of a deterministic manner and may be violated. This challenge is traditionally addressed via robust beamforming. The major solutions to robust beamforming is to provide the worst-case guarantees or probabilistic performance guarantees, such as the semidefinite relaxation (SDR) [26], Bernstein-Type Inequality (BTI) method [27], and Large Deviation Inequality (LDI) method [28]. However, these solutions require high computational complexity, which incur large latency and they are overpreservative to account for the worst cases. The proposed design is based on a low-complexity unsupervised deep learning approach with the data augmentation technique.
Our results show that it achieves significant improvement and outperforms the BTI method.
• The MPE constraint is considered in SWIPT networks from a macroscopic point-of-view through the employment of stochastic geometry. Closed-form analytical expressions are derived for the probability of satisfying the MPE constraint, the information coverage probability, the energy coverage probability as well as their joint probability. The performance with and without the constraints is considered and it is shown how each metric is affected by this restriction. Moreover, we study the system's performance for different frequency bands and our results show that higher frequencies decrease the levels of RF exposure in the network.
The rest of this paper is organized as follows: Section II describes the modeling of SWIPT and of the considered safety regulations. In Section III, the proposed design of robust beamforming for SWIPT under SAR constraints is provided together with simulation results. Section IV presents the performance analysis of a large-scale SWIPT network under MPE constraints. The paper concludes with Section V. Notation: Lower and upper case boldface letters denote vectors and matrices, respectively; [·] † is the Hermitian transpose operator; P{X} and E{X} represent the probability and expectation of X, respectively; R n denotes the n-dimensional Euclidean space; Γ(·) and Γ(·, ·) denote the complete and upper incomplete gamma function, respectively [29]; B(·, ·) denotes the beta function [29];  = √ −1 is the imaginary unit; tr(A) gives the trace of the square matrix A; A 0 means that the matrix A is positive semidefinite; ℑ{x} and ℜ{x} return the imaginary and real part of x, respectively; U (a, b) denotes the uniform distribution in the interval [a, b].

A. Information and Power Transfer
Throughout this paper, we consider multi-antenna transmitters and single-antenna receivers. Also, each receiver has SWIPT capabilities, i.e. it can decode the information but also harvest energy from the received signal simultaneously. The SWIPT technique is employed with the power splitting (PS) method such that the received signal is split into two parts: one is converted to a baseband signal for information decoding and the other is directed to the rectenna for EH and storage [30]. This is a mature SWIPT technique that does not require strict time synchronization between information and power transfer. Let ρ ∈ (0, 1) denote the PS parameter at a receiver. Then, 100ρ% of the received power is used for decoding, while the remaining power is directed to the EH circuit. During the baseband conversion phase, additional circuit noise is present due to the phase-offsets and the circuit's non-linearities, which is modeled as an additive white Gaussian noise (AWGN) with zero mean and variance N C .
Therefore, based on the PS technique considered, the SINR at a receiver can be written as where S is the received power of the signal of interest, I is the power of the interference and N 0 is the variance of the AWGN component of the received signal. On the other hand, since 100(1 − ρ)% of the received energy is used for rectification, the instantaneous energy harvested at a receiver is modeled by the following non-linear function, which refers to a specific excitation signal 1 [31] where P r is the aggregate received signal power at the receiver andā,b,c are parameters determined by the rectification circuit through curve fitting. These parameters fully characterize the non-linear behaviour of the rectifying circuit including the maximum harvesting value, the sensitivity and slope of the output power [31]. In other words, depending on how a rectifier is designed will correspond to a different set of parameters; in this work, we will considerā = 2.463,b = 1.635, and = 0.826 [31]. Note that, in general, E (measured in Watts) should be a function of the received signal rather than just its power. In this paper, we adopt a simplified model to highlight the dependency on the power of the energy signal only, which will be discussed in the numerical results.

B. Safety Constraints in Wireless Networks
The SAR metric quantifies the deposited microwave energy at a specific point on the human body. It corresponds to the rate of energy absorption per unit mass at a specific location in the tissue [18]. As such, SAR is a function of the induced electric field E (measured in V/m), the electrical conductivity of the tissue σ at the specific point (measured in S/m) as well as the tissue's density η (measured in kg/m 3 ). This relation can be expressed as In order to integrate the SAR limitations into the design of SWIPT, we use a quadratic form of the transmitted signal to model the pointwise SAR value with multiple transmit antennas [18]. This model is based on experimental studies, which showed that SAR depends significantly on the phase difference between the antennas [34] and thus can be characterized by a sinusoidal function of the phase difference [17]. Let x denote the transmitted signal and Q = E{xx † } its covariance. Then, SAR can be modeled as a quantity averaged over the transmit signals with a time-averaged quadratic constraint given by where A is the SAR matrix and P is the SAR limit. The dependence of the SAR measurements on the transmitted signals can be fully described by the SAR matrix, where the entries of this matrix have units of kg −1 . Note that the SAR matrices are positive-definite conjugate-symmetric matrices, since the SAR measurements are always real positive numbers. The SAR matrix is obtained offline during SAR testing and it highly depends on the device's type of antennas, operating frequency, industrial design, etc. [18]. Moreover, during testing, several measurements are taken for different operations and locations of the device and the one that provides the "worst case" is used for comparison with the SAR constraint. Different to the SAR metric, which is a point quantity, the MPE (power density) quantifies the RF exposure over a specific area. The MPE corresponds to either the absorbed power density or the incident power density [6]. The latter is easier to be measured and thus international regulations are given in terms of the maximum incident power density values. The MPE is given by where E [V/m] is the induced electric field and Z = 377 Ω is the impedance of free space [6]. In the context of wireless networks, the MPE can be evaluated by where P t is the transmitted power, G is the antenna gain and d is the distance of the measuring point from the center of the antenna [32]. The SAR and MPE metrics are equally important for the realization of safe (in terms of RF radiation) wireless networks, where each metric concerns a different aspect of the network. The SAR metric is considered for devices that are meant to be carried close to the body (e.g. mobile phones) [18], [34] or to be "worn" on the body (e.g. implantable sensors) [35], [36]. Essentially, these devices will be operating at a distance of 20 cm or less from a human body [17]. Moreover, the devices need to abide by the FCC's SAR regulations [11] and thus are tested before being commercially available. In the case of MPE, the devices under consideration are expected to be operating at a distance greater than 20 cm from a human body (e.g. base stations) [17]. Therefore, the MPE is evaluated based on a network of transmitting devices and all locations in the network area are required to satisfy the MPE constraint [32]. It is clear that SAR focuses more on the device itself, whereas MPE is a network-wide safety constraint. This motivates the approach taken in the following two sections. In Section III, we consider a SAR-aware beamforming optimization for a simple point-to-multipoint SWIPT system and, in Section IV, we study the performance of a SWIPT system from a macroscopic point-of-view with MPE constraints.

III. SWIPT WITH SAR CONSTRAINTS: OPTIMIZATION OF TRANSMIT BEAMFORMING
In this section, we study SAR-aware transmit beamforming optimization with both perfect and statistical CSI.

A. System Model and Problem Formulation
We consider a MISO downlink system consisting of an N tantenna transmitter with total transmit power P t and K singleantenna receivers that employ single-user detection, as shown in Fig. 1. The transmitted data symbol s k to receiver k follows the Gaussian distribution with zero mean and unit variance (i.e., E{ s k 2 } = 1), which is mapped onto the antenna array elements by the beamforming vector w k ∈ C Nt×1 . The signal design in terms of modulation, waveform and input distribution will also affect the efficiency of the RF-DC conversion [37], [38], but for simplicity, we do not consider these in the optimization.
We assume frequency non-selective block fading channels (i.e., the channel coefficients remain constant in each slot) with AWGN. Denote by h k ∈ C Nt×1 the fading coefficients between the transmitter and receiver k, which also captures the large-scale degradation effects such as path-loss and shadowing. The received baseband signal at the receiver k can be expressed as where n k denotes the AWGN component with zero mean and variance N 0 . The receivers harvest energy from the received RF signal based on the PS technique and so the SINR used for the data detection process at the k-th receiver is given by while the input to the RF-DC circuitry is where P r,k is the received power at receiver k, given by Finally, the l-th SAR constraint with a time-averaged quadratic constraint is where A l 0 is the l-th SAR matrix and P l is the l-th SAR limit.

1) Problem Formulation:
We formulate the problem as maximizing the harvested power at each user while satisfying the requirements of both SINR and total transmit power under the SAR constraints. To make the problem more tractable, we introduce a slack variable λ, and then the problem can be formulated as follows where γ k is the SINR threshold at the k-th receiver. Clearly, P1 is a non-convex problem because of both the SINR and EH constraints and thus is difficult to solve. In the next subsection, we develop an efficient convex optimization based algorithm that jointly optimizes the beamforming vectors and the PS parameters.
2) The Optimal Solution using SDP: We adopt the semidefinite programming (SDP) approach with rank relaxation to solve P1. We first define a matrix variable W k = w k w † k , and introduceλ = √ λ to recast the problem P1 as follows Note that the original objective value should beλ 2 . The advantage of the problem P2 is that it 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 packages such as CVX [41]. Once P2 is optimally solved, if the resulting solutions {W k } are all rank-1, they are the exact optimal solutions; otherwise, the solutions only provide a lower bound for the minimum required transmit power. 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 for the problem P2. In the following theorem, we show that this is indeed the case.
Theorem 1. The optimal solution to P2 satisfies rank(W k ) = 1, ∀k, i.e., the SDP relaxation is tight, and the optimal solution to the problem P1 can be recovered from {W k } via the eigenvalue decomposition.
Proof. The proof is given in Appendix A.

B. Robust Beamforming Solution using Deep Learning
It is noticed that the formulation of problem P1 and P2 are based on the assumption of perfect CSI. However, in practical systems, it is hard for the transmitter to obtain perfect CSI, which could be subject to CSI estimation errors. Without loss of generality, the general relation between the estimated CSÎ h k and the actual CSI h k is described as follows where ∆h k denotes the channel estimation error and H k denotes the set of all potential channel estimation errors. In practice, it is not feasible to know the exact channel estimation error ∆h k in prior, since ∆h k is usually not a deterministic value but a random variable. The statistical characterization regarding ∆h k can be available, for example via measurements and calibrations. Here, we model ∆h k by using the complex Gaussian distribution with zero mean and variance σ 2 h as follows Due to the uncertainty introduced by the channel estimation error, the deterministic constraints in P1 with perfect CSI become more difficult to satisfy. The key reason for choosing the probabilistic channel estimation error model over the worstcase channel estimation error model is that, worst-case channel estimation models will lead to a worst-case study for the robust beamforming problems, where the worst-case study usually provides a very conservative performance. This is because in practical systems, in order to bound the estimation errors with a threshold, this threshold value might need to be very conservative, so that all possible estimation errors (including those extremely rare cases) are bounded. Therefore in this manuscript, instead of providing a determined bound based on the worst-case scenario, we aim to provide a statistical guarantee to the EH and SINR constraints. Therefore, we consider the constraints in a statistical manner, where the original deterministic constraints are statistically guaranteed with a probability. By rewriting the corresponding terms of P1 in a probabilistic form, the robust formulation of the SWIPT problem under the SINR, total power and SAR constraints with imperfect CSI can be written as where α k and β k are the probability guarantees for the SINR and EH constraints, respectively. The robust formulation in P3 is non-convex due to the probabilistic constraints, which makes the problem NP-hard and difficult to solve.
1) The Bernstein-type inequality method: We first introduce an existing technique in the literature to solve P3, the so-called BTI method. The BTI transforms the probabilistic constraints into a deterministic form based on the large deviation inequality for complex Gaussian quadratic vector functions, which is given in the following lemma [27], [39].
Lemma 1. If the probabilistic constraint can be represented in the following form where x is a standard complex Gaussian random vector with x ∼ CN (0, I), B is a complex Hermitian matrix, z is a complex vector, while the tuple (B, z, σ) forms a set of deterministic optimization variables, and ρ ∈ (0, 1] is fixed, then the following implication holds where ψ ∈ R and ω ∈ R are slack variables, and (21) is jointly convex in B, z and σ.
In order to exploit the BTI method, we first apply the SDR with W k = w k w † k as used in the reformulation of P2, and then further rewrite the probabilistic constraints as follows ,j =k W j for notation convenience. Then, by applying the BTI method to transform the probabilistic constraints regarding the SINR and EH, we can get the following robust formulation, By the definition of W k = w k w † k , W k should be of rank one, which has been relaxed to be positive semidefinite in P4. Since P4 is a convex problem, its solutions W k can be efficiently obtained via numerical software packages such as CVX, which are always optimal to the transformed problem P4, but not necessarily optimal to the original robust formulation of the SWIPT problem P3. If the solution of W k is of rank one, then it is also the optimal solution to the original problem P1. In such cases, the solution of w k can be derived based on W k , where it is straightforward when the rank of W k is one. For the cases where the rank of W k is higher than one, the near-optimal solution of w k can be derived via the rank-one approximation of the W k , e.g. via singular value decomposition methods [40]. Since the transformed deterministic constraints are convex in W k , ρ k , x k , z k , µ k and ν k , the original robust beamforming for SWIPT in P3 has been transformed into a convex problem as in P4. Note that by using the implication in BTI in Lemma 1, the transformed constraints in P4 characterize the lower bounds on the probability α k and β k . Therefore, the feasible solutions of P4 are sub-optimal solutions of the original NP-hard problem P3, and can be efficiently solved using convex optimization solvers such as CVX, but the performance could be conservative.
2) Deep Learning Based Method: Based on the above analysis, conventional solutions such as the BTI method might transform the probabilistic constraints to a more tractable form, but at the cost that the solutions are expected to be conservative comparing to the optimal solutions. Instead of seeking sub-optimal solutions via the conventional techniques like BTI, in this subsection, we will exploit the power of the deep learning neural networks (NNs), which learn from the data to form the robust beamforming strategies for the studied SWIPT under SINR, total power and SAR constraints. A general design of the training method is illustrated in Fig. 2. Specifically, we will present the data-augmentation based technique to transform the probabilistic constraints to an NN training problem, and then reformulate the optimization problem with multiple constraints to a deep learning problem. Note that our method can deal with general channel error distributions and is not limited to Gaussian errors.
The general objective of the deep learning method is to train an NN with the estimated channelsĥ 1 , . . . ,ĥ K as inputs, and it will output the beamforming solutions of P3 as follows where θ denotes the parameter set of the NN. In this way, it transforms the original non-convex optimization problem about the beamforming {w k } to finding the optimal parameter set θ. However, note that the deep learning method cannot automatically address the constraints, which should be addressed by a deliberate design of the training procedure.
According to the SINR definition in (8), the calculation of Γ k requires the estimated CSI h k , the beamforming vector w j for j = 1, . . . , K, and the PS parameter ρ k . Without loss of generality, we can rewrite Γ k as a function of the NN f (ĥ 1 , . . . ,ĥ K ; θ) as follows: Similarly, the EH power of receiver k in (9) can be rewritten as The SAR function can be also rewritten as follows a) Problem Reformulation via Quantile Functions: Similar to the analytical studies in the previous sections, it is also a challenging problem to address the probabilistic constraints in the NN. To facilitate the evaluation of the probabilistic constraints, we first introduce the quantile function q(x, σ) as follows, where the quantile function q(x, σ) returns the quantile value (infimum value) z, such that for all x, the probability P{x ≤ z} is no more than the value of σ. In this way, given a probabilistic constraint in the following form then it can be rewritten by the quantile function in the equivalent form as follows which transforms the comparison from the probability P{x ≤ Z} against σ, to the quantile value q(x, σ) against the threshold Z. Therefore, for the SINR constraint (17) of receiver k, it can be first rewritten into the following form Then, with the quantile function, it can be further transformed as follows where the NN representation of Γ k in (24) has been used. Similarly, the EH constraint (18) of receiver k can be transformed as follows By noticing λ is also the objective of the original robust formulation in P3, we further rewrite the robust formulation of SWIPT based on NN f (ĥ 1 , . . . ,ĥ K ; θ) as follows S l (f (ĥ1, . . . ,ĥK ; θ)) ≤ P l , ∀l, where the mathematical expectations in objective (33) and the SINR constraint (34) are with respect toĥ k , ∆h k , ∀k. This is to make the trained NN parameter set θ generally applicable to all possible ∆h k and estimated CSI inputsĥ k . Note that since w k and ρ k are the outputs of the NN f (ĥ 1 , . . . ,ĥ K ; θ), the deterministic constraints (35)-(37) can be regarded as the constraints on the NN outputs. b) Addressing the constraints: During the training procedure, the NN cannot automatically satisfy the constraints. Therefore, for the problems with constraints, e.g., the studied robust beamforming for the SWIPT problem, each constraint needs to be addressed deliberately. Similar to the analysis in the convex optimization, there are no universal solutions to address all constraints. Some simple constraints can be addressed via the design of the output layer of the NN architecture. Here we will use two output layer design techniques to address the constraint (35) and (36) as follows: • For the PS constraint in (35), each ρ k should be within the range of [0, 1]. This can be achieved by applying the Sigmoid function in the output layer for each raw output ρ k as ρ k = Sigmoid(ρ k ), with Sigmoid(x) = 1 1+e −x . • For the total power constraint in (36), the sum power of all w k , calculated by K k=1 w k 2 , should be bounded by the total power P t . This can be addressed by scaling the raw beamforming vectorsw k as w k = Yw k , with Y = min 1, With the above output layer design, the constraints (35) and (36) are enforced by the NN architecture automatically, i.e., all outputs will satisfy both constraints. Therefore, when possible, it is preferable to exploit the NN architectures, e.g., the NN output layer, which firmly address the constraints. However, the above techniques can only address simple constraints such as (35) and (36). For complicated ones, such as the probabilistic constraint (34) and the SAR constraint (37), we exploit a general technique by modifying the objective function, so that the trained NN parameter set θ should learn to satisfy both via the training procedure. This is achieved by considering the penalty of violating the constraints in the loss function given by (38), where c 1 , c 2 and c 3 are positive weight parameters for each individual learning objective terms 2 , and (x) + max{x, 0} is the clamp operation. The use of the clamp operation is to make the loss function L(θ) increase only when the corresponding constraints are violated.
S l (f (ĥ1, . . . ,ĥK ; θ)) − P l + . (38) In this way, we can rewrite the robust beamforming for SWIPT under SINR, total power and SAR constraints as the following unconstrained deep learning problem: where P6 provides a transformed formulation of P3 that can be solved by deep learning methods. Specifically, P3 is first transformed to P5 with the help of the quantile function in (27), which is then transformed to P6 by integrating the constraints together with the objective function as the loss function that can be used for deep learning. Note that since the loss function L(θ) exploits the end performance (the expected harvested power) to evaluate the outputs of the NN, the training process can apply the unsupervised training method, where there is no need to know the optimal beamforming vectors and the PS parameters as supervised labeled data in the supervised training method. c) Data-augmentation Based Training Method: In the studied robust beamforming problem of SWIPT with SINR, total power and SAR constraints, the most challenging task is to address the probabilistic constraints with regards to the EH and SINR constraints. In the previous sections, we have exploited the quantile function to transform the constraint with probability to a constraint with quantile values. Further, the loss function has been modified so that it is expected that these constraints can be learned by the NN. It can be observed by the loss function L(θ) in (38), that the channel estimation errors ∆h k is only required to evaluate the outputs of the NN with regard to the probabilistic constraints, while the NN only calculates outputs based on the estimated CSIh k .
Inspired by this, we can add an auxiliary module during the NN training procedure, where each estimated CSIh k is augmented by a set of channel estimation errors S ∆h k to form a set of potential actual CSI S h k = {h k |h k = h k + ∆h k , ∀∆h k ∈ S ∆h k }, which is then used to evaluate the probabilistic related terms in the loss function L(θ). The evaluation over the augmented set S h k , provides an estimated performance of the loss function L(θ), while it is expected to converge to the actual performance when the size of the set increases according to the law of large numbers. Similarly, the mathematical expectation calculation in (38) against the estimated CSIh k can be achieved by an evaluation of the loss function over a set of estimated CSI Sh k . With the loss function L(θ) estimated over the estimated CSI set Sh k and the augmented set S h k for each element in Sh k , the NN parameter set θ can be updated based on the gradient descent method [42] as (40).
When the offline training phase is completed, the trained NN can be deployed for the application and the beamforming solution can be obtained by using the inference mode of the NN, i.e. the channels are used as inputs to the NN and the NN only performs forward calculations without a back-propagation process. Since matrix multiplication is the most computationintensive operation, we estimate the computational complexity based on the required multiplications, while other operation time is ignored. If the NN has L fully connected layers and the l-th layer has θ l neurons, then the total computational complexity of the NN-based method can be estimated as O( L l=2 θ l−1 θ l ). Besides the fixed number of neurons for the layers l − 2 to L − 1, the first layer input is of dimension 2KN t and the last layer output is of dimension 2K(N t + 1). Therefore, the total computational complexity of the NN for the proposed method can be estimated as O(KN t ). The complexity of solving the general problems P2 and P4 is dominated by the SDP constraints and according to [43, 6.6.3], the associated complexity of the interior-point algorithm for solving these two problems is O

C. Simulation Results
To evaluate the performance of the proposed algorithm, simulation results are presented and discussed in this section. The considered MISO downlink is composed of a transmitter with N t = 3 transmit antennas and K = 2 receivers. Each transmit and receive antenna has 8 dBi and 3 dBi gain, respectively. Receivers can harvest energy at frequency f = 915 MHz, and are randomly located around the transmitter with a distance l k ∼ U (1, 5) m and at a direction ζ k ∼ U (−π, π). We adopt Rician fading with a Rician factor of 0.5 to model the channel, due to the short distance between the transmitter and the receivers and the dominance of the line-of-sight (LOS) signal; the path loss coefficient is 2.5. We consider one SAR constraint, i.e., L = 1, P t = 2 W, N 0 = −70 dBm and  (41) 1) Optimal Beamforming with Perfect CSI: For performance comparison with the proposed optimal solution with perfect CSI, we consider as benchmarks the zero-forcing beamforming (ZF-BF) and the solution without the SAR consideration (i.e., the solution of P1 without considering the SAR constraint). Specifically, the ZF-BF vector is given by , and p k is the power for the k-th receiver. Let G k,j |h † kw j | 2 denote the equivalent link gain between the transmitter and the kth receiver, which satisfies G k,j = 0, ∀k = j. Finally, let F k,l =w † k A lwk denote the l-th radiation channel gain due to the transmission intended for the k-th receiver.   This is a linear programming problem and can be solved using CVX. But as will be illustrated by simulation results, the ZF-BF solutions are very conservative in the EH performance. Fig. 3 shows the harvested power versus the SAR power constraint P l , where the SINR requirement at the receivers is 15 dB. The proposed optimal solution achieves higher harvested power as P l increases, and significantly outperforms the ZF-BF scheme. For instance, when P l = 2 W, the harvested energy of the proposed solution is about −5.6 dBm, which is 1.6 dB greater than that of the ZF-BF solution. The solution without a SAR constraint remains constant and achieves a higher harvested power than the proposed solution, but the performance gap decreases by increasing P l . Fig. 4 depicts the total power consumption at the transmitter for various P l when P t = 2 W. We can observe that the solution without a SAR constraint consumes the full transmission power (2 W) to maximize the harvested power. The proposed algorithm makes more effective use of the transmit power than the ZF-BF solution in order to harvest more power while satisfying the SAR constraint. Fig. 5 shows the harvested power versus the SINR requirements, where the SAR power constraint is 2 W. The harvested power of all solutions decreases as the SINR constraints increase, and this indicates that more power of the received signal is used for information decoding at the high SINR requirement, due to the nature of PS. Fig. 6 demonstrates the feasibility probability of the three solutions in terms of the SINR requirements, and the SAR power constraint is 2 W. When the SINR requirement is low (i.e., 5 dB), the feasibility probability of the optimal solution (about 99%) is close to that of the solution without SAR constraints and much higher than that of the ZF-BF solution (around 84%). In the high SINR regime, this gap reduces as the SINR requirement increases and this is because the ZF-BF becomes nearly optimal so both solutions converge.

P7: max
2) Robust Beamforming with Imperfect CSI: Simulations are conducted to compare the performance of robust beamforming with the three previously discussed solutions when the available CSI is imperfect, which includes a) the non-robust solutions to the SWIPT formulation P2, without considering the channel estimation errors, which is solved via CVX, b) the BTI solutions to the robust formulation P4, and c) the proposed NN-based solution to the robust formulation P6 via the proposed method in Fig. 2. The variance of the channel estimation error is set as σ 2 h = 10 −5 . Specifically, the SINR and SAR thresholds are the same for all receivers, where γ k = 5 dB and P l = 1.2 W/kg, while the total transmit power P t = 2 W is used. The probability constraints are set as α k = 95%, β k = 90%, ∀k, unless otherwise specified.
In the proposed NN training method, five fully connected hidden layers are constructed as illustrated in Fig. 2. Each layer has a width of 300, and the PReLU activation function is used except for the last layer, where PReLU(x) returns x if x ≥ 0, or 0.25x otherwise. To support the NN training, the Batch normalization technique is used for each layer except for the last layer, and the Adam optimizer is used with a learning rate of 10 −5 . Since the proposed NN training method is based on unsupervised training, it only requires to generate the estimated CSI and channel estimation errors for the training purpose. In the simulations, we have randomly generated 4×10 5 estimated CSIs and exploited the mini-batch training method with a batch size of 4000. During the training procedure, the channel estimation error sets are randomly generated according to its distribution, where each estimated CSI is augmented with 200 channel estimation errors. The weight parameters are empirically selected as c 1 = 10000, c 2 = 100 and c 3 = 100. In addition, a total of 500 testing estimated CSIs are generated, which are used as inputs for the three considered methods, whose beamforming solutions are then evaluated against 10 4 channel estimation errors for each estimated CSI. We first evaluate the EH performance of the three methods, where the cumulative distribution function (CDF) of the least harvested energy among the receivers is shown in Fig. 7. It is seen that the NN-based method shows a better EH performance than the BTI method, where the harvested power at the 90% probability threshold is −12 dBm and −15.7 dBm, respectively. This corresponds to a 3.7 dB gain for the NNbased method compared to the BTI method. As for the nonrobust method, it shows the best EH performance among the three methods, but later it will be shown that the constraints have been violated, which makes the EH performance invalid to use.
The worst receiver's SINR performance is presented in Fig.  8, where the probability threshold of 95% is indicated by the horizontal dashed line, and the 5 dB SINR threshold is indicated by the vertical solid line. Specifically, it is seen that both the NN-based method and the BTI method meet the probabilistic constraints regarding SINR, where their worst receiver's SINR have been maintained above 5 dB bound at the 95% probability. Meanwhile, the NN-based method shows better performance compared to the BTI method, i.e. it provides a higher worst user's SINR at the probability threshold 95%. This is due to the fact that the BTI method transforms the original robust formulation to a convex but conservative formulation, where the solutions can be effectively solved but at the cost of the sub-optimal solutions. On the other hand, the NN-based method empirically evaluates the probabilistic bounds during the training, and it learns to target at better EH performance, while satisfying the SINR requirements. As demonstrated in Figs. 7 and 8, both the NN-based and BTI methods, satisfy the robust beamforming constraints as in the original problem P3, while the NN-based method provides better solutions with higher EH performance and the BTI solutions are more conservative compared to the NN solutions. It is noticed that the non-robust method fails to satisfy the SINR constraint, and this is because the non-robust method does not consider the channel estimation errors during its  Fig. 9. The SAR performance comparison between the proposed NN method, the BTI method, and the non-robust method.
formulation. By jointly considering the EH performance in Fig. 7 and the SINR performance in Fig. 8, the results of the non-robust method also indicate the importance of the robust formulation, as the uncertainty introduced by the channel estimation errors could compromise the performance, or even make the solutions invalid to use. The SAR performance is shown in Fig. 9. It can be seen that the CDF of the BTI method and the non-robust method overlap with the threshold line at P l = 1.2 W/kg. This is because the SAR constraint is a linear constraint, while the transformed or relaxed formulations in both the BTI and the non-robust method are in convex forms. Therefore, the overlap with the SAR threshold P l = 1.2 W/kg indicates that the solutions via the BTI and non-robust method are based on the equality condition of the SAR constraint. On the contrary, the CDF of the NN method shows that the final SAR performance is within the range from 0.6 to 1.2 W/kg, and the SAR threshold has been met for all tests with regards to both estimated CSI and channel estimation errors. It is also noticed that although the formulation only requires a firm bound on the maximum SAR, the NN-based method can provide solutions with lower SAR radiations.
In summary, the trade-off among the SINR, EH and SAR has been shown to be a non-convex and NP-hard problem in the analytical studies in Section III-B, but the results in Figs. 7-9 demonstrate that the proposed training method has successfully learned better solutions than the BTI method and the non-robust method.

IV. SWIPT WITH MPE CONSTRAINTS: LARGE-SCALE PERFORMANCE ANALYSIS
In this section, we study a SWIPT network from a largescale point-of-view under the MPE constraint.

A. System model 1) Network Model:
We consider a large-scale bipolar ad hoc wireless network consisting of a random number of transmitter-receiver pairs. The transmitters form a homogeneous Poisson point process (PPP) Φ = {x i : i ≥ 1} of density λ in a two dimensional Euclidean space R 2 , where x i ∈ R 2 denotes the location of the i-th transmitter. Each transmitter x i has a dedicated receiver at a distance d 0 in some random direction and transmits with fixed power P t . The time is considered to be slotted and at each time slot all the transmitters are active without any coordination or scheduling. We consider the performance of a receiver located at the origin and its associated transmitter x 0 . We perform our analysis for this typical receiver but, according to Slivnyak's Theorem [44], our results hold for any receiver in the network.
2) Channel Model: We assume that all wireless links suffer from both small-scale block fading and large-scale path-loss effects. The link between a receiver and an interfering node is in LOS with a probability p L , otherwise it is blocked, e.g. by buildings. On the other hand, the link between a transmitter and its dedicated receiver is always considered to be in LOS, i.e. p L = 1. The interference effect from non-LOS signals is ignored, as we assume the dominant interference is caused by the LOS signals [45]. We consider independent Nakagami fading with parameter µ for each LOS link and so the power of the channel fading is a normalized gamma random variable with shape parameter µ and scale parameter 1/µ. We denote by h i the channel gain for the link between the i-th transmitter and the typical receiver. Moreover, all wireless links exhibit AWGN with variance N 0 . The path-loss model assumes that the received power is proportional to d −α i where d i is the Euclidean distance from the origin to the i-th transmitter and α > 2 is the path-loss exponent.
3) Sectorized Antenna Model: The transmitters and receivers are equipped with multiple antennas and employ adaptive directional beamforming [46]. For the sake of analytical tractability, we make use of an approximation of an actual beam pattern using a sectorized model. In particular, the gain of a link between a transmitter and a receiver is a discrete random variable given by [46] where g takes into account three parameters: the main lobe beamwidth ω ∈ [0, π], the main lobe gain M , and the side lobe gain m. We let q i = {( ω π ) 2 , 2( ω π )(1− ω π ), (1− ω π ) 2 } denote the probability of a link with gain g i = {M 2 , M m, m 2 }. Finally, we assume that the link gain between each transmitter and its dedicated receiver is equal to M 2 , i.e. they are perfectly aligned. Due to this sectorization, the PPP Φ is partitioned into three thinned spatial processes [47], as follows: Φ 1 represents the set of interferers with link gain g 1 with the typical receiver, i.e. Φ 1 is a PPP with density λ 1 = q 1 λ. Similarly, Φ 2 and Φ 3 are the set of interferers with link gains g 2 and g 3 , respectively, with the typical receiver; as such, Φ 2 and Φ 3 are PPPs with density λ 2 = q 2 λ and λ 3 = q 3 λ, respectively.

B. Performance Analysis with MPE Constraints
Based on the considered system model, the SINR at the typical receiver can be written as where P 0 g 1 P t and denotes the aggregate interference generated by the transmitters in Φ at the typical receiver, with P i = g i P t . On the other hand, the instantaneous energy harvested at the typical receiver is given by (2) with which is the aggregate received signal power at the receiver. Any potential RF EH from the AWGN noise is considered to be negligible.
In this section, we focus on the MPE of the network, evaluated at the origin. In other words, we would like to study the probability of satisfying the MPE constraint τ , expressed as where which follows from (6) and is referred to as the point source model, where the transmitting antenna is assumed to be represented by a single point source [48]. Even though this model does not take into account the antenna size (assumed to be a point), it is accurate in the far-field [6], [48]. We first state the following lemma and then provide the main result.
Lemma 2. The characteristic function of the interference I is given by where λ and P are the density and transmit power of the interfering nodes, respectively.
Proof. The proof is given in Appendix B.
Theorem 2. The probability of satisfying the MPE constraint τ is given by and φ(t, λ i , P i , α + 2) is given by Lemma 2 with λ i = p L q i λ and P i = g i P t .
Proof. The proof is given in Appendix C.
Note that as µ increases, the above probability converges to a constant floor. Specifically, for µ → ∞, we have where (52) follows from (1 − x) a → 1 − ax for x → 0 and (53) follows from Γ(x + a) → Γ(x)x a for x → ∞.
Next, we are interested in studying both information and energy coverage probabilities when the MPE constraint is satisfied, i.e. the joint probabilities. Due to the correlation between the SINR and the energy harvested with MPE, we assume these events are independent for the sake of analytical tractability. Therefore, we consider the bounds where γ and ǫ are non-negative thresholds for the SINR and the average harvested energy, respectively, and P{SINR > γ} and P{E > ǫ} are given in the following corollary.
Corollary 1. The information coverage probability is given by and the energy coverage probability is given by The proof for the above expressions follows similar steps as the proof of Theorem 2 and thus it is omitted. It is important to point out that in (58), we need ǫ <ā−b/c, that corresponds to the value at which the harvesting circuit saturates. We now derive the joint distribution of the MPE, the SINR and the average harvested energy. In other words, we evaluate Theorem 3. The joint information and energy coverage probability is given by where ψ(t) = 3 i=1 φ(t, λ i , P i , α), δ is given by (58) and Proof. The proof is given in Appendix D.
Finally, we consider the special case where the interference does not exist. This can be realized in two ways: a) in a blockage dense area, i.e. p L → 0, and b) small beamwidth and small side lobe gain, i.e. ω → 0 and m → 0. In this case, we havê where Ξ max γN 0 + NC ρ , δ . The proof is omitted as it follows directly from the use of the CDF of a gamma random variable. By using the above, we can find the optimal P t that maximizes the joint distribution. The derivative is expressed as which follows from dΓ(a, x)/dx = −x a−1 exp(−x). Then, by setting the derivative equal to zero and solving for P t , we deduce that Observe that the optimal value of P t is independent of the fading parameter µ. As expected, P * t increases with τ but decreases with Ξ and M .

C. Numerical Results
We now evaluate the performance of SWIPT networks with MPE constraints and validate our analytical framework with Monte Carlo simulations. For the sake of comparison, we consider networks operating in millimeter wave (mmWave) bands (e.g. 30 GHz) and in ultra high frequencies (UHF) (i.e. 300 MHz to 3 GHz). Following a similar approach to [46], we use M = 10 dB for mmWave and M = 0 dB for UHF, i.e. the mmWave antenna gain is ten-fold the one of UHF. Moreover, we consider µ = 5, N 0 = −117 dB, p L = 0.8 (mmWave) and µ = 1, N 0 = −127 dB, p L = 1 (UHF). Finally, the remaining parameters are set as: ω = π/6, m = −10 dB, d 0 = 5 m, α = 3, N C = 0 dB, γ = −10 dB, ǫ = −5 dB and ρ = 0.5. Fig. 10 depicts the probability of satisfying the MPE constraint. As expected, the probability decreases with the transmit power and increases with the constraint τ . We can observe that the mmWave band satisfies the constraint with a higher probability, compared to the UHF band with the same transmit power. Indeed, for P t = 10W and τ = 0.2 W/m 2 , mmWave can satisfy the constraint around 75% of the time, whereas the RF exposure with UHF below τ is less than 60% of the time. This shows the positive impact of directional beamforming, which can be achieved by higher frequencies, on the network's overall RF exposure. The figure also depicts the asymptotic scenario µ → ∞ (Eq. (52)). It can be seen that this provides a lower bound on p s (τ ). Finally, the analytical results (lines) and simulation results (markers) are in agreement, which verifies our analysis. Fig. 11 illustrates the information coverage probability in terms of the transmit power. When there is no MPE constraint, the probability converges to a constant ceiling for large values of P t ; this corresponds to the interferencelimited scenario. Observe that mmWave networks significantly outperform UHF, as also shown in [46]. On the other hand, when an MPE constraint is imposed, the coverage probability decreases after a certain value of P t as the MPE level is exceeded more frequently. In other words, there is an optimal value that maximizes the coverage probability, which can be easily deduced by algorithms such as the bisection method. Similar observations can be derived for the energy coverage, shown in Fig. 12. In particular, with no MPE constraints, the probability converges to one for high values of P t , whereas it drops after a certain value of P t when an MPE level is enforced. However, observe that the optimal P t here is different compared to the information case. Moreover, the UHF bands perform well in terms of energy coverage but are still outperformed by mmWave bands due to the higher antenna gain at the receivers and transmitters.
Finally, Fig. 13 shows the joint coverage probability with respect to the LOS probability. It can be seen that the impact of LOS probability depends on the value of the transmit power. At low P t , the performance improves with p L since this facilitates in harvesting more RF energy. Note that for these values, the performance loss from imposing safety constraints is small. On the contrary, for higher values of P t , the coverage probability decreases with p L as the effect of interference is detrimental to the SINR. Also, the losses in performance due to MPE constraints are more notable in this case.
V. CONCLUSION In this paper, we provided an framework for the design and analysis of far-field SWIPT under safety constraints. We focused on two RF exposure regulations, the SAR and the MPE, and outlined the state-of-the-art as well as the modeling approach in the context of communication networks. A design for optimal robust beamforming based on deep learning was proposed, subject to specific information, EH and SAR constraints. In addition, a complete theoretical study for the performance of large-scale SWIPT systems under the MPE constraint was derived, with regards to both information and energy coverage. Our results provide insights in terms of the optimal SWIPT design and show the potentials from the proper development of SWIPT systems under health and safety restrictions. APPENDIX

A. Proof of Theorem 1
First the partial Lagrangian of the problem P2 is written as where {α k , β k , ν l , µ} are dual variables, and we have defined 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 x † (µI + K j=1 (α j − β j (1 − ρ j ))h j h † j + L l=1 ν l A l )x = 0. Because X k 0, we have Therefore, it holds true that h † k x = 0, ∀k. 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 Karush-Kuhn-Tucker condition of the problem P2 is that tr(W k X k ) = 0, so the rank of W k is at most 1. This completes the proof.

B. Proof of Lemma 2
The characteristic function φ(t, λ, P, α) of the interference I = P x∈Φ h x d −α x is given by φ(t, λ, P, α) = E{exp (tI)} = E Φ,hx exp tP where (68a) follows from the probability generating functional of a PPP [44]; (68b) from the moment generating function of a gamma random variable since h x are independent and identically distributed. By using the transformation x = − tP µu α and the binomial theorem, the integral can be written as which follows from [29, 3.194-3]. With the use of the integral representation of the beta function [29, 8.380-1], we can write B(k − 2 α , µ − k + 2 α ) = 1 0 x − 2 α −1 (1 − x) µ+ 2 α −1 ( x 1−x ) k dx. Then, it is easy to show that the above finite sum is equal to −B(− 2 α , µ + 2 α ), which completes the proof.

C. Proof of Theorem 2
We will derive the probability of satisfying the MPE constraint τ by applying the Gil-Pelaez inversion theorem [49], that is, where φ(t) is the characteristic function of the expression in (48) evaluated at t. Therefore, we have Since (48) is the sum of four independent terms, its characteristic function is given by the product of the characteristic function of each term. For the first term, we have which follows from the fact that h 0 is a gamma random variable with shape and scale parameters µ and 1/µ, respectively. Finally, the characteristic function of the i-th term in the above sum is φ(t, λ i , P i , α + 2), given by Lemma 2.

D. Proof of Theorem 3
The joint distribution can be evaluated as p J (γ, ǫ) = P{SINR > γ, E > ǫ} By letting A where (75a) uses the probability density function f h (h) = µ µ h µ−1 exp(−µh)/Γ(µ) of a gamma random variable with parameters µ and 1/µ. The lower limit ξ is derived by considering and solving for h 0 gives (61). Then, in (75b), φ(t, λ i , P i , α) is given by Lemma 2 and which follows from [29, 3.381-3]. The second term in (74) can be evaluated with a similar way and the result follows after some algebraic manipulations.