On the Optimal Deployment of Power Beacons for Massive Wireless Energy Transfer

Wireless energy transfer (WET) is emerging as an enabling green technology for Internet of Things (IoT) networks. WET allows the IoT devices to wirelessly recharge their batteries with energy from external sources such as dedicated radio frequency transmitters called power beacons (PBs). In this paper, we investigate the optimal deployment of PBs that guarantees a network-wide energy outage constraint. Optimal positions for the PBs are determined by maximizing the average incident power for the worst location in the service area since no information about the sensor deployment is provided. Such network planning guarantees the fairest harvesting performance for all the IoT devices. Numerical simulations evidence that our proposed optimization framework improves the energy supply reliability compared to benchmark schemes. Additionally, we show that although both, the number of deployed PBs and the number of antennas per PB, introduce performance improvements, the former has a dominant role. Finally, our proposal allows to extend the coverage area while keeping the total power budget fixed, which additionally reduces the level of electromagnetic radiation in the vicinity of PBs.


I. INTRODUCTION
The broad range of today Internet of Things (IoT) applications demands a massive deployment of low-cost devices, powered by small batteries, and with different quality-ofservice (QoS) requirements. Battery lifetime is critical, since maintenance operations could be impractical, if not forbidden, either in massive or harsh conditions deployments [1]. In order to extend the battery lifetime, or even avoid completely its need, energy harvesting techniques, from renewable sources such as solar and wind, or from ubiquitous radio-frequency (RF) signals, have been recently proposed and investigated [2]. Dedicated RF sources provide an alternative for supplying controllable and reliable energy for low-power IoT devices. Wireless energy transfer (WET) techniques are therefore important enablers of many coming IoT use cases [3]. Although there are many WET techniques, herein we focus exclusively on RF WET, and hereafter WET refers to RF-based energy transfer.
In [4], the authors consider a scenario where PBs transfer energy to rechargeable devices considering practical technology limitations. Therein, they were interested in maximizing the energy transferred while constraining the maximum level of radio-frequency-electromagnetic field (RF-EMF) due to its risk to human health. A similar concern is exposed in [5], where the authors maximize the charging utility function of the energy harvesters (EHs) while constraining the level of RF-EMF through a meta-distribution characterization of the per-sensor received power.
Nature-inspired algorithms have been recently proposed in [6] for solving the placement problem of PBs. Therein, the authors propose a WET solution for extending the lifetime of multihop networks. They propose an algorithm that mimics the flocking behaviour in order to maximize the network lifetime by equalizing the lifetime of individual sensors. Moreover, the use of multiple antennas provides diversity over the space as well as additional degrees of freedom, suitable to tackle the impairments caused by small-scale fading. This applies to wireless powered communication networks (WPCNs), where beamforming techniques can be implemented to increase the coverage area due to the directionality in transmission [7]. In [8], it is considered an energy beamforming technique to tackle the problem of long-distance power transfer in largescale multiple-input multiple-output (MIMO) systems. In [9], a wirelessly powered massive MIMO system, where a multi antenna base station (BS) charges single-antenna EHs in the downlink, is optimized. A distributed wireless power transfer system is studied in [10] and [11] where the focus is optimal distributed energy beamforming to overcome the effects of destructive interference whether frequency and phase synchronization is available or not. They conclude that distributed antenna systems are more effective than collocated antennas at a single PB as long as the optimal inter-PB beamforming can be achieved. To improve the energy transfer efficiency in a large WPCN, [12] proposes an adaptative energy beamforming strategy that maximizes the average received power by exploiting the tradeoff between the power intensity of the beams and the number of sensors to be charged. In [13], the authors propose a charge scheduling scheme using batterypowered PBs which harvest energy from a central PB. They exploit the benefits of beamforming in order to maximize the energy efficiency while satisfying QoS constraint of the sensor network. A max-min rate problem is addressed in [14], where solar-powered multi-antenna PBs assist an IoT network. Therein, the authors propose a distributed rate allocation protocol that considers the state of the data buffer and battery level information of each device and its nearest neighbors, in order to forward the collected measurements.
However, above solutions require perfect channel state information (CSI) acquisition for estimating the optimal beamforming vector, which consumes more energy as the number of EHs increases [15]. In an attempt to fill this gap, the authors of [16], [17] study a WET setup where a multiantenna PB operates without CSI for powering massive lowpower IoT deployments. They analyze the statistics of the harvested energy under the sensitivity and saturation phenomena, and compare it with that experienced when using CSIbased schemes. Particularly, in [17] the authors do consider the phase shift between antenna elements and compare it with the assumption of independent transmitted signals.
Hybrid networks deployments under outage constraints are investigated in [18], where authors assume that single-antenna BSs and PBs form independent homogeneous Poisson point processes (PPP). Interesting tradeoffs between transmit power and density of BSs and PBs are derived under an outage constraint on the data link. Multiple antennas at the BS and PBs are considered in [19], where the latter are deployed following a PPP with certain density.
Therein, the authors investigate the area that should be covered by the centered BS to minimize the average power consumption, so that the uplink rate is above certain threshold with probability given by the QoS requirement of all devices in the network.

A. Related works
Recently, PBs deployment optimization emerged as a potential technique to ban blind-spots and homogenize the network energy availability. For instance, in [20] the authors take into account directional antenna patterns to optimize the placement of PBs in order to maximize a charge utility function for all deployed sensors. Therein, they consider the benefits of directivity over the omnidirectional radiation, and propose an algorithm to dynamically place the PBs as the network topology changes. Similarly, the authors in [21] optimize the positions of PBs considering the presence of obstacles. Therein, they maximize the network-wide weighted sum harvested energy with heterogeneous PBs and sensors.
In order to enhance coverage, Daubechies wavelet algorithm is proposed in [22] to minimize the number of PBs to be deployed. Assuming the sensors' deployment to be known, they find the optimum PBs' positions by moving the redundant PBs at each iteration towards the uncovered sensors. Meanwhile, the authors in [23] study the problem of mobile directional PBs' deployment to improve energy coverage. Each PB can move inside certain area such that the overall harvested energy is maximized. In [24], the authors optimize the number and locations of energy and information access points in a WPCN subject to energy harvesting and communication performance requirements. The authors associate a cost to each node's installation to then minimize the overall network cost deployment. They address the network performance given average channel gains and considering wireless devices' positions to be known a priori, which allows to divide the network into nonoverlapping clusters. Hence, the proposed algorithms yield different solutions depending on the sensors' locations for the same network area, and therefore it could not be optimal if these locations are not predetermined. Finally, the position and speed of a single mobile PB powering a sensors' deployment is optimized in [25]. In this case, authors don't rely on the sensors' positions. Instead, they divide the service area into smaller sub-partitions and focus on the performance of the instantaneous worst position, which changes as the PB moves. The goal is to minimize the overall charging time constraining the movement of the PB.

B. Motivation and contributions of this work
Different from above works, herein, we study the optimal deployment of PBs to satisfy an energy outage constraint. We consider a massive deployment of EHs for which CSI acquisition procedures are too costly, therefore, instantaneous CSI is not available. Since neither CSI nor devices locations are known, PBs' positions are optimized for maximizing the foreseen minimum average RF energy available in the area. This guarantees for each EH to most likely meet its corresponding QoS requirements. Moreover, our proposed strategy does not depend on neither the hardware heterogeneity nor the mobility of the sensors. Table I lists the main differences of our model compared with previous related works.
Our main contributions are summarized below • We provide a framework for assessing the performance of a WET setup operating with the minimum number of PBs required to support energy QoS requirements; • We discuss and evaluate several methods for finding the optimal deployment of PBs for powering a massive number of sensors; • Our results suggest that in case of one and two PBs, both should be deployed in the center of the R radius area under evaluation, while in case of three and four PBs, they should be symmetrically located in a concentric circumference approximately of radius R 2 and R √ 2 2 , respectively. For a greater number of PBs, their optimal positions depend strictly on the area and path loss exponent; • Numerical results show that the number of PBs deployed in the network has a greater impact on the system performance than the number of antennas per PB, even though both introduce improvements in the energy outage.
C. Organization of the paper Next, Section II introduces the system model and presents the problem formulation. In Section III, we discuss several approaches for finding the optimal positions for a fixed set of PBs, while in Section IV we show and analyze numerical performance results as a function of the number of PBs and available transmit antennas. Finally, Section V concludes the paper.
Notation: Here, we use boldface lowercase letters to denote column vectors. The operator | · | can represent either the absolute value for scalars or the cardinality of a set, while x p = ( ∀i≥1 |x p i |) 1/p denotes the ℓ p -norm [26, eq. (1)]. Given a minimization problem with objective function f 0 (x), m inequality constraints f i (x) ≤ 0, i = 1, . . . , m,

Reference
Sensors' positions Antennas PBs Mobility Objective [20] Known Directional Static maximize the average harvested energy [21] Known Directional Static maximize the weighted-sum of the network-wide harvested energy [22] Known Onmidirectional Static minimize the number of PBs that guarantees full coverage [23] Known Directional Mobile maximize the network-wide harvested energy [24] Known Onmidirectional Static minimize the cost of network deployment [25] Not necessarily known Directional Mobile minimize the charging time Current work Unknown Onmidirectional Static minimize the number of PBs that satisfies an energy outage constraint barrier parameter and n equality constraints h i (x), i = 1, . . . , n, we denote as the Lagrangian of the problem [27]. Herein, ν i and λ i are the Lagrange multipliers of the i th inequality and equality constraints, respectively. The symbol ∇ represents the gradient operator, i.e. ∇f (x 0 ) : R n → R n is the gradient of the function f (x) : R n → R at the point x 0 . Moreover, O(·) is the big-O notation, which specifies worst-case complexity. Finally, z ∼ N (µ, σ 2 ) is a Gaussian distributed random variable with mean E[z] = µ and variance σ 2 , while y ∼ χ 2 (k, φ) is a non-central chi-squared random variable with k degrees of freedom, and noncentrality parameter φ [28]. Table. II lists the symbols used throughout the paper. Fig. 1 where a set B = {PB b |b = 1, 2, . . . , |B|} of PBs are deployed in a circular area of radius R to wirelessly power a massive number of EHs, denoted by a set S = {EH s |s = 1, 2, . . . , |S|}. Each PB is equipped with an omnidirectional antenna, and transmits at the same power level P . We also assume quasi-static channels with independent and identically distributed (i.i.d) Rician fading with factor κ. The channel coefficient between EH s ∈ S and the PB b ∈ B is denoted as h s,b ∈ C, and h s,b = α + jβ with independent real and imaginary parts α, β ∼ N κ 2(1+κ) , 1 2(1+κ) [29].

Consider the scenario in
There is no knowledge about sensors' positions, and PBs are not able to estimate the CSI due to the large number of devices and/or their stringent associated energy expenditure constraints [16], [17]. We denote the energy-carrying signal coming from the PB b as x b and we consider that independent signals are transmitted from all the sources. Then, the incident RF power at the s th harvester is comprises the path-loss of the EH s → PB b link with distance d s,b , which depends on the path-loss exponent γ and other factors considered in K such as the carrier frequency and antenna gains [30]. Finally, while P is normalized in the time domain, which allows us to indistinctly use the terms energy and power for the same quantity. Authors in [16, eq.(24)] found that the distribution of the harvested energy, when a single-antenna PB serves the network, is a scaled non-central chi-squared with non-centrality parameter 2κ, and 2 degrees of freedom. We can extend those results to our model and write which is a linear combination of i.i.d non-central chi-squared random variables. The sum (2) is a particular way of defining the generalized chi-square distribution, whose probability density function is often computed by numerical algorithms due to its analytical intractability [31].

A. Problem formulation
We consider a massive deployment of EHs in the area. The main goal is to estimate the minimum number of PBs that need to be deployed in order to satisfy an energy outage probability constraint with threshold ζ. The ultimate incident RF energy ξ RF s must be above the sensors' sensitivity ξ 0 with probability 1 − ζ, ∀s ∈ S. The optimization problem can be formulated as follows Since the domain of the objective is |B| ∈ N, P1 is a mixedinteger programming problem which is in general nondeterministic polynomial (NP) complete [32]. However, we can overcome this issue by searching over the space of candidate solutions for which the positions of the PBs are optimized. Additionally, the constraint (3b) is by far cumbersome, so we rely on Monte Carlo simulations when computing the energy outage probability. Finally, (3c) refers to the per-PB maximum transmit power given the network power constraint P T .

III. OPTIMAL PBS POSITIONING
The problem P1 relies on the optimal positions of PBs given certain propagation conditions. Thus, let us fix |B| to first find the optimal PBs positions. Without loss of generality, we discretize the circular region where sensors are massively located, hence, mimicking its deployment. This is done by creating evenly spaced circumferences with a discrete number of points proportional to its radius. Let u s , n b ∈ R 2×1 denote coordinate vectors for the locations of EH s and PB b , respectively. Notice that as |S| increases, the points' deployment tends to cover the whole area. We can state the average RF energy available at the s th sensor, as a function of the distance to each PB as Notice that, in order to meet the system constraints, it is sufficient that the EH performing the worst, e.g., the one under the greatest path loss, meets the energy outage requirement. Then, the PBs positioning optimization problem can be stated as which maximizes the average received energy at the worst sensor through the optimal deployment of PBs within the coverage area.

A. Equally-far-from-center approach (EC)
We can take advantage of region's symmetry properties to reduce the complexity of the problem. Herein, we assume the PBs to be equally-far-from the circle center at a distance n b 2 = r : 0 ≤ r ≤ R, ∀b ∈ B. Note that the angle θ determined by any two adjacent PBs and vertex at the circle center must obey θ = arccos Under these conditions, the average receive power is the minimum either at the circle edge or at the circle center. The deployment is illustrated in Fig. 2 for |B| = 3. The average received power at the s th sensor on the region's edge using polar coordinates is where ϕ is the angle of the sensor located on the edge. Since the set of signals {x b } are independent, we can take advantage of the superposition principle considering first the contribution of PB 1 and PB 2 in EH s to find the worst position on the edge. We find the critical points by taking considering the contribution of PB 1 and PB 2 in (6), and after some algebraic manipulations, we attain where The roots of (7) are determined by setting its numerator to 0, which has a trivial solution at ϕ 1,2 = {θ/2, θ/2 + π} within This result is expected, since it refers to the most distant sensors in terms of symmetric power contribution considering two adjacent PBs. In fact, the path gain is monotonically decreasing with the distance, and the energy provided to the sensors located at (R, 0) and (R, θ) represent local maxima of that provided to the sensors located at the circle edge, therefore we conclude that the solutions {(R, ϕ 1 ), (R, ϕ 2 )} correspond to local minima. Once we take into account the contribution of PB 3 , the point (R, ϕ 2 ) is turned into a local maximum, while (R, ϕ 1 ) remains as minimum. In the general case |B| ≥ 3, for |B| odd, (R, ϕ 1 ) is a local maximum provided the contribution of the PB located at (r, θ/2 + π), otherwise is local minimum equivalent to (R, ϕ 2 ). Fig. 3 summarizes these ideas. Notice that, there is an alternate pattern of equally-spaced |B| local minima and |B| local maxima at the circle edge. Therefore, we adopt (R, ϕ 1 ) for finding the value of r that maximizes the minimum average receive power. Then, let us plug this result into (6) but considering all PBs to get the optimum position r by solving The second fractional term in (9) corresponds to the contribution of PB 3 , and can be neglected as the path loss exponent increases, which gives a nearly optimal r ≈ R 2 . Following the same procedure we can arrive to a general analytical approximation r ≈ R cos θ 2 , which suggests that for the case of |B| ∈ {1, 2}, the placement should be at the center, whereas for |B| = 4 the placement is r ≈ R √ 2 2 roughly independent of the propagation conditions. It also guarantees that the minimum contribution will be at the circle's edge rather than at the center. For the general case when |B| ≥ 5 is considered, we can use previous results such that r is optimized now for the node with the minimum average received power, which can be either the center or (R, ϕ 1 ). Unfortunately, the solution of this problem doesn't guarantee the best result for an arbitrary B, and solving for the case when the optimum kind of deployment is not known a priori is mathematically intractable. Our proposal is to solve the EC approach algorithmically for two topologies: i) as in Fig. 2 here called EC; ii) with one PB, the last one, located at the center, thus n |B| 2 = 0, here called EC with one centered PB. Fig. 4 depicts the latter topology where n b 2 = r, ∀b ∈ B, b = |B|, and the angular separation between them is θ = 2π/(|B| − 1).
According to Fig. 2, the average received power for the node with s * = arg min where E[ξ RF e ] is the contribution at the worst sensor on the edge with ϕ = ϕ 1 , thus given by (6), and E[ξ RF c ] = |B|P Kr −γ is the contribution at the center. Meanwhile, for the deployment in Fig. 4 E where represents the contribution in a sensor equidistant to the center and two adjacent PBs at a distance x = r/(2 cos θ 2 ) and ϕ = ϕ 1 . Additionally, is the average power received at the worst sensor on the edge with ϕ = ϕ 1 . Finally, the optimal positions correspond to the constellation that maximizes E[ξ RF s * ] against the minimum average contribution when the network has one PB radiating from the circle center, with total power P T = |B|P , i.e.
Herein, E[ξ RF s * ] is defined as a parametric function using both (10) and (11), thus determined by the chosen topology. The procedure for efficiently determining r * is detailed in the Optimal DEployment of POwer BEaconS (Ode-PoBes) algorithm. Notice that ∆r denotes the step size of the iterative search. Fig. 5 depicts a comparison using the results of Compute E[ξ RF s * ] using both (10) and (11) 6: r * ← r 10: end if 11: r ← r + ∆r 12: until r ≥ R the Ode-PoBes algorithm and the analytical approximation n b 2 ≈ R cos θ 2 , ∀b ∈ B for the scenario shown in Fig. 2. In general the approximation holds up to |B| = 5 for γ = 3, and |B| = 4 for γ = 5, without an important degradation of the network performance. However, the gap increases with both |B| and γ, because this approach aims to benefit more the sensors in the circle's edge at the cost of degrading the contribution at the center. With the increase of |B|, the EC with one centered PB becomes better than EC without PB at the center, having a prompt transition when γ = 5.
We now provide a worst case complexity analysis for the proposed algorithm. First, for a given ∆r the required number of iterations is ⌊ R ∆r ⌋+1 and the maximum solution error is ∆r 2 . Notice that the argument E[ξ RF e ] determines the computational cost of (10), and (6) proportionally scales up with |B|. Similarly, each argument in (11) has a cost proportional to |B|. Then, the most costly step in Ode-PoBes, step 6, requires O(|B|) operations, which in turn constitutes the computational cost of each iteration of the Ode-PoBes algorithm.

B. Interior-point method approach (IPM)
Interior-point methods, for solving inequality-constrained optimization problems rely on gradient-based equations, therefore require a differentiable objective function [27]. Objective function in (5a) is neither smooth nor differentiable, therefore we resort to generalized mean approximation which can be used as a smooth estimate of the min function [33]: where ξ represents a 1 × |S| vector comprising the average power received at each s, thus ξ s = E[ξ RF s ]. The previous relation tends to be an equality as k → −∞ but at the expense of an increasing computational cost. Other differentiable approximations are addressed in [34]  Here k is a real constant that determines whether the approximation tends to the minimum function (k < 0) or to the maximum one (k > 0). All these approximations perform alike, but since the gradient of the objective must be provided to the IPM algorithm, we adopt (15) whose derivative is more tractable. We can now transform P2 to solve it using IPM based on the logarithmic barrier function as P2.1 : minimize for turning the set of |B| inequalities in P1 into equality constraints. Moreover, µ > 0 is the barrier parameter that must be chosen small enough to improve the accuracy of the approximation. The next step involves the solution of the system of equations derived from the Karush-Kuhn-Tucker (KKT) conditions [27] ∇L({n b }, t, λ) = 0, ∀b ∈ B, where each entry of λ = {λ b } corresponds to the Lagrange multiplier associated with the b th equality constraint. The components of the Lagrangian in (17) are given by where In general, the system (or the problem) is solved through a sequence of linear equality constrained quadratic problem applying the Newton's method [35], or the conjugate gradient method [36]. The performance of the IPM vs k is evaluated in Table III for |B| ∈ {3, 9, 15} using normalized values. For instance, we normalize the elements of each triplet with respect to max k E[ξ RF s * ] for given γ and |B|. As k decreases, the approximation (15) converges in almost all cases; however, it also requires more computational effort, but still for small |k|, the results are quite good. Notice that for k = −30 IPM becomes unstable for large values of the path-loss exponent, therefore compromising the reliability of the solution. Hereafter, we adopt k = −25 since it gives the best results in terms of stability and accuracy.

C. Evolutionary computation algorithms
Modern meta-heuristics methods such as genetic algorithms (GA) and particle swarm optimization (PSO) are stochastic approaches, thus, suitable for dealing with optimization problems where the objective function is highly nonlinear and nondifferentiable. GA is inspired by the process of evolutionary biology, in which subsequent generations evolve from previous throughout selection, crossover and mutation [37]. At each step, a group of individuals of the current population are selected for the next generation according to their performance against the objective function. Each of them is a possible solution of the optimization problem at a given iteration. They might also be used for crossover and mutation which helps to search over the feasible solution set avoiding to get trapped at a local optimum. This allows us to solve directly the problem P2, without the approximation in (15).
Swarm intelligence-based algorithms such as PSO use the principles of self-organization of the so-called agent particles [37]. Each particle represents a candidate solution of the optimization problem. PSO aims to find the global minimum by moving the agents in a quasi-static manner within the domain of an objective function, starting from initial positions and velocity. Considering a set Q of particles, the position of q ∈ Q at the i th iteration is {z q,i,b ∈ R 2×1 |∀b ∈ B}, where z q,i,b represents its b th component. Therefore, each particle's position set is made up of a candidate solution in P2. Let us define the objective function as where E[ξ RF s ] is computed according with (4) but replacing n b for z q,i,b , and the second term is a barrier function whose impact on the problem solution is limited by making ǫ → 0. At the iteration i ′ , the particles update their position and velocity considering its initial values, local best position pbest(q, i ′ ), Algorithm 2 Computation of the minimum |B| to find the best of all local solutions until the objective no longer improves. Once P2 is solved, one needs to minimize the number of deployed PBs that satisfies the energy outage probability requirement in (3b). Algorithm 2 shows the iterative procedure for finding the minimum |B|. At each iteration, P2 is solved until the energy outage condition is guaranteed.
Remember that the distribution of ξ RF s obeys (2). In the next section we present optimal deployment results under different network conditions.

D. Practical implementation considerations
In practice, the performance of the algorithm for finding the PBs' positions depends on how well the path loss model fits the actual conditions. For instance, recent measurement campaigns in wireless sensor networks corroborates the dependence of the path loss exponent on the environment's characteristics [39]. Hence, different sub-regions of the network might require different path loss models. Moreover, the inherent properties of IoT networks (e.g. small antenna heights, low transmission power, and stationary nodes) limit the applicability of traditional propagation models [40]. Similarly, the proposed algorithm framework guarantees meeting the energy outage requirement provided the channel fading distribution is accurately known beforehand.
Therefore, opportunistically selecting appropriate channel models, at both large (path-loss) and small scale (fading distribution), is essential. Machine learning methods are potential candidates for such a task, that must also provide confidence when a large set measurements is available [41]. The risk of channel modeling/prediction errors must be taken into account in the optimization framework as well. Finally, we can utilize the solutions from Ode-PoBes as a good initial guess for any of the metaheuristic algorithms before treated, which can improve the performance of the final deployment.

IV. NUMERICAL RESULTS
In this section, we present numerical results on the optimal PBs deployment to meet a certain energy outage constraint. Additionally, we provide insights on the maximum coverage area and the impact of multi-antenna schemes for massive WET. Unless we state the contrary, the simulations are based on the parameters listed in Table IV.

A. On the optimal deployment of PBs
The optimal positions after solving P2 are represented in Fig. 6 for different number of PBs; while we show the corresponding average power along the circle area as a heat map. The optimal PBs' positions are equidistant from the circle' center, form a symmetric distributed pattern with more rings as either the number of PB increases or the propagation conditions worsen. As we include more PBs in the network, the performance of the worst sensor improves in terms of average received power E[ξ RF s * ] under the different optimization approaches, which can be also observed in Fig. 7. Here, as a benchmark, we also present the case of a centered PB radiating with total power P T , which is equivalent to place all |B| PBs at the center. In general all methods outperform the benchmark, while the performance gap increases with |B|. The special case |B| = 2 is equivalent to place a single PB at the circle' center as all methods agree, although doubling the power. Ode-PoBes and IPM stand out over the others in terms of stability in the final solution. The optimization is limited to 15 PBs due to the poor convergence of GA compared with the IPM method. In fact, the classical GA approach doesn't converge towards the desired solution as the complexity of the problem increases. A more stable outcome is obtained with PSO, as compared with IPM as a reference.
The results obtained with Ode-PoBes follow the trend of the ones with the IPM, but the computational cost is smaller for the former. Indeed, Ode-PoBes searches for the optimum positioning without requiring derivative computations. Fig. 8 shows the normalized average convergence time with respect to the time required for Ode-PoBes. Notice that the Ode-PoBes algorithm is superior by four order of magnitude since it reduces the dimension of the optimization variable. Although the other methods presented searches for individual positions there is a gain in the stochastic approaches over the IPM, given they do not require to compute the gradient at each iteration. However, the gap is reduced as |B| increases, i.e. as the number of possible constellations increases. Finally, for |B| ≥ 7 the evolutionary algorithms approaches exhibit a close performance.

B. On the solutions of P1
Let us consider the scenario where wireless charging is provided as a service to the IoT devices, and each PB is powered by an external wired source. In order to obtain profits, the service provider chooses to distribute a fixed power budget of P T = 10 W among the individual PBs according to (3c).  Hereinafter, IPM will be our benchmark strategy provided its accuracy over GA and PSO. The monotonic decrease of the energy outage probability for both approaches, proves the effectiveness of the distributed PBs approach when improving the network reliability.
For this setup, we also present the solutions of P1 in Fig. 10 as a function of ζ and R. Fig. 10 (top) shows that tightening the QoS agreements requires more PBs to be deployed; whereas Fig. 10 (bottom) depicts a stable increment in the number of PBs needed to meet the target QoS constraint as the network gets larger. In general, the fluctuations in the minimum number of PBs (min |B|) reflect the dominant impact of the distancedependent loss. Notice that Ode-PoBes outperform IPM's solutions when serving larger areas with tight QoS requirements. Particularly, Ode-PoBes finds an optimal deployment with 4 PBs less than IPM when γ = 5 and R = 100 m. Moreover, Ode-PoBes makes possible to cover wider areas as compared to the IPM based solutions with the same QoS requirements. The gap between Ode-PoBes and IPM is due to numerical approximations that affect the symmetry of the final solution. In fact, IPM solves |B| individual positions within the space of possible constellations, which doesn't guarantee symmetric layers with respect to the origin. Meanwhile, Ode-PoBes finds (r * , θ * ) that maximizes the E[ξ RF s * ] using two predefined symmetric constellations.

C. On the maximum coverage area
In Fig. 11, we show the maximum coverage area results with respect to E[ξ RF s * ] using the minimum average incident RF power obtained with Ode-PoBes. Clearly, the distributed deployment outperforms the centered PB approach, and the performance gap increases with the number of PBs and the path-loss exponent. Our proposal shows that bringing the PBs closer to the worst-positioned sensors is more effective to overcome the distance-dependent loss than increasing the power of a single PB. Additionally, as the sensitivity ξ 0 of the harvesting circuitry increases, not just the coverage area but the gap among different configurations diminish. In fact, as the power budget is shared by more PBs, the received power at reference distance of 1 m decreases by a factor of |B| with respect to the centered PB approach. Moreover, it allows to extend the coverage area without increasing the level of RF-EMF in the proximity of the PBs. Notice that in Fig. 11 (top) the points |B| ∈ {3, 7, 10} have been highlighted, so that the readers can observe the correspondence with the results in Fig. 11 (bottom) when ξ 0 = −22 dBm.

D. On CSI-free multi-antenna WET
Performance of multiple-antenna strategies for powering a massive number of energy efficient devices are analyzed in [16] using one antenna (OA), all antennas at once (AA) and switching antennas (SA) schemes. We only focus here on the SA scheme given that OA can be seen as an special case when only a single antenna is used. On the other hand, the performance claimed for AA in [16] is valid only for equal mean phases along the transmit antennas, which is difficult to hold in practice. Under SA strategy, each PB transmits with full power by one antenna at the time, in a way that the whole array is utilized within a channel coherence block. Additionally, SA preserves the harvested energy as in the case of single-antenna PBs (although boosting the energy diversity), harvested while the variance is a function of the spatial correlation among the antennas.
For this setup, let us consider that each PB is equipped with multiple antennas denoted by the set A. We keep the  assumption on the fading distribution for each antenna at PBs to transmit independent energy-carrying signals. The use of SA scheme, as Fig. 12 depicts, improves the system performance in terms of reducing the outage probability at the worst sensor as the number of antennas increases. As an example, we constraint the system to ζ = 10 −4 , and it can be observed that the number of PBs have a greater impact on the system performance than the number of antennas. In fact, the outage probability does not improve significantly as the number of antennas increases when |B| = 1, but once more PBs are distributed in the network, the performance impact of the number of antennas grows significantly. The readers can observe that the deployment of an additional PB could roughly be compensated by doubling the number of antennas at the existing PBs.

V. CONCLUSION
In this paper, we studied the mixed-integer programming problem of PBs' deployment optimization with aim to wirelessly charge a massive number of IoT devices under energy outage constraints. We first found the optimal PBs' positions by using a heuristic method based on equally-far-fromcenter restrictions, and then we solve it algorithmically. As a benchmark, we compute the optimal positioning problem using interior-point methods by approximating the objective with the generalized mean, and also through evolutionary algorithms, without prior knowledge about devices' locations. Numerical results show that our algorithm performs similar to the benchmark approaches, but it converges much faster given the constraints imposed to the PBs' positions. As compared with the centered PB approach, properly distributed PB deployments introduce improvements in terms of energy outage probability and coverage area. Additionally, it was shown the dominant impact of the distance-dependent loss on the solution of the minimization of the number of required PBs under a certain outage constraint as the radius increases. Finally, we presented energy outage results using CSI-free WET schemes for the case of multiple antennas at the PBs. Although having multi-antenna PBs always introduces improvements, we found that the number of deployed PBs impacts stronger the system performance. Indeed, considering the same total power, the performance obtained when doubling the number of antennas at the PBs can be alternatively attained by properly deploying an additional PB which reduces the hardware complexity. Our results provide valuable insights for designing practical WET setups by answering how many PBs are needed, and their corresponding locations, for powering certain area with average power or energy outage QoS constraints.
An attractive future work is to study the impact of knowing the position of the sensors in the PBs' deployment optimization. Notice that the deployment's symmetry assumption is crucial in this work. Deviations from the optimal deployment caused by obstacles or forbidden installation places can significantly degrade the incident power at the worst position. Hence, we could consider the possibility of having moving PBs in order to provide flexibility to network changes. Different from previous works, we can optimize the PBs' positions with the help of probabilistic machine learning tools, which provide low-complexity solutions to large problems.