Massive Wireless Energy Transfer with Statistical CSI Beamforming

Wireless energy transfer (WET) is a promising solution to enable massive machine-type communications (mMTC) with low-complexity and low-powered wireless devices. Given the energy restrictions of the devices, instant channel state information at the transmitter (CSIT) is not expected to be available in practical WET-enabled mMTC. However, because it is common that the terminals appear spatially clustered, some degree of spatial correlation between their channels to the base station (BS) is expected to occur. The paper considers a massive antenna array at the BS for WET that only has access to i) the first and second order statistics of the Rician channel component of the multiple-input multiple-output (MIMO) channel and also to ii) the line-of-sight MIMO component. The optimal precoding scheme that maximizes the total energy available to the single-antenna devices is derived considering a continuous alphabet for the precoders, permitting any modulated or deterministic waveform. This may lead to some devices in the clusters being assigned a low fraction of the total available power in the cluster, creating a rather uneven situation among them. Consequently, a fairness criterion is introduced, imposing a minimum amount of power allocated to the terminals. A piece-wise linear harvesting circuit is considered at the terminals, with both saturation and a minimum sensitivity, and a constrained version of the precoder is also proposed by solving a non-linear programming problem. A paramount benefit of the constrained precoder is the encompassment of fairness in the power allocation to the different clusters. Moreover, given the polynomial complexity increase of the proposed unconstrained precoder, and the observed linear gain of the system's available sum-power with an increasing number of antennas at the ULA, the use of massive antenna arrays is desirable.


I. INTRODUCTION
The internet of things (IoT) presupposes a large number of low-complexity devices, most of them wireless, and a number of them placed in hardly accessible locations (e.g., for infrastructure monitoring). Even though most of these terminals are not supposed to transmit continuously, they are supposed to be able to operate for long periods of time without batteries replacement. One way of achieving that is by having devices that can harvest energy from incoming electromagnetic This work was supported by the Academy of Finland (Aka) (Grants n.307492, n.318927 (6Genesis Flagship), n.319008 (EE-IoT)). The work of Francisco Monteiro was supported by a joint scholarship from Fundación Carolina and Fundación Endesa (both from Spain) and funded by FCT/MCTES (Portugal) through national funds and when applicable co-funded EU funds under the project UIDB/50008/2020. F. A. Monteiro  radiation. Much work has been recently done on simultaneous wireless information and power transmission (SWIPT) [1], wireless-powered communication networks (WPCN), in which terminals harvest energy during a first time slot and use it to power the transmission of information in a second time slot [2], and also on systems dedicated to wireless energy transfer (WET) per se [3]- [12]. WET systems differ from the former two by only focusing on the transfer of energy without spending resources on any type of data transfer.
As it is thoroughly described in [13], the idea of WET is as old as radio transmission but, until quite recently, the idea typically involved high transmission power and very large antennas at both ends of a link, and often aimed at very long distances (e.g., satellite-to-Earth links). Using low-power WET to power devices over a few meters (tens or even hundreds of meters) is a quite recent endeavor. The WET approach is easier to deploy massively given that the energy source, or power beacon (also still called base station (BS) for legacy reasons), does not need to have any connection to a core data network apart from the connection to the electrical grid. Note that the WET component of a system is often the one responsible for performance bottleneck in practical WPCN and SWIPT systems, and for that reason it is worthy of dedicated study and consideration of different optimization approaches. Other examples of systems that may also need to rely on WET are backscatter communications [14], and the new long-distance tags based on quantum tunneling radio positioning [15].
In [5], one can find a survey of beamforming techniques for SWIPT, WPCN, and WET systems, and also a list of some techniques to easy the enormous task of acquiring channel state information at the transmitter (CSIT), on which those beamforming techniques depend on. However, instantaneous CSIT availability is an unrealistic assumption in the context of mMTC because not only the devices are very energyconstrained, but also due to the immense number of pilots the BS would have to deal with in the uplink. Some steps forward to do away with instantaneous CSIT have been recently taken in [16], where fading is artificially accentuated such that its benefits can be further exploited even if CSIT is entirely absent, and in [8], where a statistical analysis of the harvested energy in a WET setup assuming no CSIT is presented. In [8], the authors assessed WET with single antenna terminals and a multiple antenna transmit array at the BS, and compared i) single antenna transmission from the BS, ii) equal power transmission from all the BS's antenna elements, dubbed as all antenna at once (AA) scheme, and iii) the switching antennas (SA) scheme, where only one antenna is active at a time and the position of the active one runs through all the antennas arXiv:2106.07988v2 [eess.SP] 11 Jul 2021 in the array during each block of coherent fading. Their channel model had all devices experiencing a Rician fading channel and distributed around the BS without any particular distinction between devices.
A great part of the literature on WET-enabled systems (WET only or SWIPT) assumes that full instantaneous channel knowledge is available at the BS, e.g., [9], [17]- [21], or even simultaneously available at an intermediate intelligent reflective surface [22]. Nevertheless, imperfect CSIT has also been considered in a a number of works [6], [7], [9], [11], [12]. In general, acquiring a more accurate information about the channel requires more energy expenditure, and consequently there is a trade-off between the energy that is harvested by the terminals and the energy spent feed-backing CSI. The works in [6] and [9] considered point-to-point WET, with a massive antenna array wirelessly powering only one EH multiantenna terminal. In [6] the EH device has multiple antennas and the aim was the maximization of the net harvested energy by optimizing the training phase and the WET phase. In [9] the authors considered a multi-input single-output (MISO) system, i.e., with a single-antenna EH terminal, and derived expressions for the outage probability of EH, defining outage when the harvested energy during a block is less than the sum of the energy the sensor spent in transmitting the CSI and the energy spent processing its tasks. In both cases a Rician channel model was considered and the channel matrix is written as the sum of a known matrix (representing the partial CSIT) with a unknown remaining term for a Rician channel, which is the most common way of incorporating partial CSIT in the system models. The work in [23] also considers massive antenna array at the BS but rather a multi-user (MU) system with single-antenna EH terminals, and formulates partial CSIT in the same way. Similarly, in the context of SWIPT, the lack of full-CSIT has also been studied in the same manner, e.g., in [24] a MU-SWIPT system is considered with single-antenna terminals.
This paper considers partial CSIT, formulated in a distinct manner, such that the partial knowledge comprises of statistical knowledge about the channel, which has been argued to be beneficial for MIMO WET systems in [4]. Note that in openloop WET systems (i.e., CSIT-free), which are oblivious of the channel state and operate without any CSIT [8], [16], [25], only a statistical analysis of the harvested energy is possible, and a sizeable statistical analysis of such open-loop systems currently exists. In [10], a fair beamforming scheme for WET that only makes use of the channel's first order statistics (i.e., the channels' mean) was proposed, also under Rician fading. The present paper proposes reduced-complexity precoding that also does not require instantaneous full-CSIT but only the first and second order statistics of the channel, considering that both the channel correlation and the LOS components vary slowly.
In many practical scenarios the terminals are physically clustered (as in Fig. 1) and thus, besides sharing the same large-scale gain (pathloss) from the BS to the terminals, they also share a common slow-fading component, which can be characterized by its second order statistics, i.e., the covariance of the channel matrix. This paper proposes to take advantage of that partial CSIT, which can be estimated with pilots transmitted from the terminals to the BS at a low rate.
The Karhunen-Loève (KL) channel representation allows one to use beamforming techniques reminiscent of the initial ideas for low-CSIT MIMO precoding with reduced CSI feedback [26]. Similar considerations about the correlation between the channels experienced by nearby terminals lead to a proposal combining massive MU-MIMO, MIMO precoding, and user clustering [27]. In that scenario, beamforming is performed for the different clusters of users, and then the intracluster users use non-orthogonal multiple access (NOMA) in the power domain to detect their own signals.
It should be noted that the proposed system is a MU-MIMO setup in which each cluster assumes the role of a virtual single-antenna terminal, representing all the terminals sharing a common second order statistics and a similar line-of-sight (LOS) model. This is in contrast with [9], where a MISO system is considered for WET in order to highlight the benefits of using a massive antenna array.
The energy beams can be adapted to focus onto each cluster and, because interference is not a foe in the context of WET, the optimal precoder can be designed by maximizing the total system's sum-power available to be harvested at the terminals' antenna in all clusters. However, in doing so, an unfair partial sum-power allocation among different clusters may happen and it is shown that by framing the problem as a constrained numerical optimization, fairness can be achieved among clusters. A different approach to introduce fairness in a WET MU-MIMO scenario was taken in [18], where the authors proposed designing a precoder which takes in consideration fairness by applying a max-min criterion to maximize the network's user with the weakest channel conditions, while considering full-CSIT and Rayleigh fading. In that work, the authors considered a different pathloss for each device, given that they are not considered to be clustered around each other. To the best of our knowledge, a WET system model considering clusters of terminals has only appeared in the very recent papers [25], [28], but the fact that they are considering clusters has no implications for precoding in these works, which assess open-loop techniques such as all antennas at once (AA) and switching antennas (SA).
One focuses on the available power to be harvested rather than the harvested energy after the energy harvesting (EH) circuit in order to offer more general results. Nevertheless, the specifications of a particular piece-wise linear EH circuit (to be detailed in section II) are considered in this work in order to tune the system to a particular useful power range.
Besides the limitation due to the high sensitivity of the EH circuits (e.g., around −10 dBm, compared to around −60 dBm when the goal is information transmission), there is also the downside of the inefficiency of the voltage multipliers based in diodes chiefly responsible for converting the radiofrequency energy (RF) into direct current (DC) voltages for both immediate use and storage, as highlighted in [29], and thoroughly described in [30]. A description of different antennas designs for the EH terminals and a compilation of EH circuits characteristics, such as frequency range and power conversion efficiency and be found in [31].
The harvested DC power is not only a non-linear function of the RF power at the antenna for a given transmitted waveform, but also a function of the waveform itself [32]. For an easier approach and analysis of energy beamforming, previous works have assumed a fully-linear conversion, i.e., considering a linear function for the RD-DC conversion, and not even considering the sensitivity and saturation effects of the circuits, e.g., [33]. The same simplification was made to design and analyze systems with a WET phase followed by an information transfer phase from sensors [34]- [37], and for MU-MIMO WET with beamforming [18], as in the scenario considered in the present paper. In fact, this simple linear model can be traced back to [38]. Nevertheless, more realistic, yet more challenging, non-linear RF-DC conversion models have been considered in the literature, ranging from characterizing the non-linearity by a sigmoid function [39], to assuming a general RF-DC transfer function that only needs to be monotone and increasing in [40]. A model with saturation-only was proposed in [12], [41], and also used in [42] for beamforming in SWIPT. A significant step forward was taken in [43], by characterizing the RF-DC conversion while keeping analytical tractability. In doing that, the authors have found that the (unmodulated) single sinewave was not the optimal waveform for WET when the non-linearity is considered and proposed deterministic multisine waveforms. These waveforms, which have a higher PAPR (peak-to-average power ratio) when the number of sines increases, eventually deliver a higher output DC power that scales linearly with the number of sine waves even without CSIT, as long as the the channel fading is flat. Under frequency-selective fading the linear increase is also attainable as long as CSIT is present [43]. Interestingly, the benefits of using multisine waveforms for WET had been previously hinted in [44], where the authors proposed using multi-sinewaves to better match the Taylor expansion of the diode's non-linear function, and also for [45] for RFID tags, in which the authors proposed using a multi-tone waveform. Later, in [46], motivated by the advantages of using multisine waveforms in WET, the author extended the idea to SWIPT systems, with the transmitted symbols resulting from the sum of a deterministic multisine waveform with a modulated multi-carrier waveform, as in orthogonal frequency division multiplexing (OFDM).
The precoders proposed in the present paper are based on statistical CSIT (mean and correlation) and attain power gains very close to the optimal precoder with instantaneous full-CSIT available at the BS. The results are in all setups much superior to the ones obtained by the CSIT-free techniques [25]. The type of precoding herein proposed provides a linear gain for the system's sum-power as a function of the number of antenna elements at the BS, naturally making the case for the use of massive antenna arrays for WET whenever terminals clustering is possible. The proposed techniques hold for any waveforms that one may consider, transmitted over flat fading, and therefore the power conversion gains coming from an optimized waveform, such as the multisine waveforms ones in [43], [47] or any of the waveforms listed in [48, Table 1] for different scenarios, will add up to the gains coming from the proposed precoding techniques applied at the BS, preserving the performance gaps found in the current work.
The remainder of the paper starts by setting the system model in section II, then section III lays out the proposed precoding schemes, and does so for increasingly more sophisticated cases. That section also tackles the fairness problem of power allocation among clusters. Section IV presents a baseline reference case that is analytically tractable, and which will serve to validate the numerical simulator used to assess the proposed techniques in section V. Finally, the conclusions are drawn in section VI.
Notation: a complex circularly symmetric Gaussian random vector with mean m and covariance R is denoted by CN (m, R). A H denotes the Hermitian transpose of matrix A, I n is an identity matrix with n diagonal elements, 1 M is the column vector of M ones, and E is the expectation operator. u(x) is the Heaviside step function, and Γ(x) is the complete Gamma function. The probability density function (PDF) of a random variable X is denoted by p(x) and its power is considered to be equal to its second moment E{X 2 }, assuming a unit resistor to convert Volt 2 to Watt, as it is usual in signal processing. To a set of independent and identically distributed random variables one applies the i.i.d. acronym.

II. SYSTEM MODEL
The considered scenario for WET is the one in Fig. 1, where, without loss of generality regarding the antenna's geometry, a BS equipped with M -antennas serves L clusters, each of which encompassing K terminals. The proposals in the paper will be assessed with a uniform linear array (ULA) at the BS, but this option does not preclude the possibility of considering a more general array, such as a uniform planar array (UPA) [49], which would give another degree of freedom to the system, allowing for vertical beam steering. Moreover, a scenario considering a different number K l of devices per cluster is straightforward to assess and can easily be generalized from the model that is here assumed for illustration purpose. In fact, both precoders that will be derived for multi-cluster systems (in sections III-D and III-E) include the information about the clusters in K × M matrices that could also be K l ×M with no impact on how the optimization unfolds.
A Rician MIMO channel [50,Sec. 5.7], [51, Sec. 3.6.1], is considered to exist from the BS to each of the terminals within a cluster. The baseband model of the signals received by the K single-antenna terminals within the l-th cluster are aggregated in the vector where the elements y l (k) in vector y l ∈ C K×1 are the received signals at each one of the terminals in the cluster and β(l) is the mean large-scale pathloss from the BS to the terminals in cluster l. Matrices P l ∈ C M ×M will be restricted to be unitary matrices when deriving the unconstrained optimal beamforming precoder that focus the transmit power onto the l-th cluster for efficient WET. That restriction will be dropped when searching for a constrained solution.
The vector onto which the unitary transformation is applied is set to x 0 = P x /M s, where s ∈ C M ×1 is the waveform signal vector in the baseband model, such that E{||s m || 2 } = 1, for m = 1, . . . , M , where each s m represents a waveform, possibly optimized for a particular model of the terminal's EH circuits, such as a deterministic multisine signal [43], [47], [52] (a joint optimization of the precoder and of the waveform falls beyond the scope of this work).
As indicated in (1), the equivalent precoding vector applied to the M antennas is x l P l x 0 , and therefore the transmitted signal to power all the network devices from the BS is is the vector containing the thermal noises at each single-antenna terminal, each of which with σ 2 n power. As usual in the WET literature, this power is considered to be negligible for the purpose of energy harvesting and is neglected in all aspects of the remaining of the paper.
The matrices H (LOS) l and H l , both ∈ C K×M , respectively represent the LOS and the multi-path components of the Rician MIMO channel from the BS to the l-th cluster, and κ(l) is the Rician factor defining the weights of each component. The clusters' angular width is considered to be narrow and therefore all the terminals approximately experience the same LOS component, sharing the same κ and the same slowly varying H (LOS) l , which is assumed to be perfectly known at the BS via low-rate feedback.
It is interesting to note that while in the cases in which NOMA is considered to bring some additional capacity to systems using massive MIMO, the users need to be confined to a narrow angular spread [53], [54], which is quite unnatural scenario in cellular communications, the assumption of a narrow angular spread and small range differences seem reasonable in several applications of WET, when particular spots have a higher concentration of EH devices. Matrix H l is a Rayleigh fading matrix with all its elements taken from a circularly symmetric complex Gaussian distribution with all elements h k,m ∼ CN (0, 1), while H ( Because all the terminals in a given cluster approximately experience the same LOS component, The power gain from the ith antenna to the k-th terminal is characterized by the same β(l) to all users in the l-th cluster, given that all antennas of the ULA are co-located, as it has been also considered in [6], [7] and in the results section of [9]. φ l describes the angular position of the l-th cluster, measured in respect to the so-called endfire direction of the ULA (in [10] it was used the angle measured to the ULA's boresight (also referred to as broadside) [51,Sec. 3.5]). Therefore, θ i (φ) = 2πid cos(φ)/λ, for α i , i = 1, . . . , M − 1, where λ is the wavelength, and d = λ/2 is the separation between the antenna elements of the ULA. Note that the common phase ϕ l plays no role in the transferred power to the clusters. A more general case with the K terminals within each cluster randomly distributed within an angular domain φ l − ∆ φ , φ l + ∆ φ will also be assessed. In that case, the k-th terminal in a given cluster is located at an angle φ k , taken from a uniform distribution within the angular aperture of ∆ φ degrees, and one has different rows , for each of the K terminals in the l-th cluster.
While in most MTC scenarios the line-of-sight component is quasi-static, the fast component H l is difficult to be obtained by the BS in a mMTC context, not only because of the power energy limitations of the devices, but also due to the sheer number of terminals. Nonetheless, a central assumption in this paper is that the correlation of the H l component is assumed to be known at the BS because it can be slowly updated over time by means of pilots sparsely transmitted by the devices, or via modern machine learning methods [55]. Given the physical proximity of the users within a cluster, according to the geometrical one-ring scattering channel model [56] there is a transmit antenna correlation among the signals received by those K terminals. Thus, the Rayleigh matrix component of the MIMO channel of a cluster has a square transmit correlation matrix R l = E{H H l H l } ∈ C M ×M [51], which, via its singular value decomposition (SVD), can be expanded in the form: It will be considered that the rank of R l , r l , is equal to the number of users, similarly to [57]. Therefore, Λ l is only going to have r l non-zero diagonal elements. When the channel correlation is known one can take advantage of the KL channel representation, which has been applied in the context of NOMA [27], [58], MIMO spacial division multiplexing [56], rate-splitting in MU-MIMO systems [59], and also as an efficient way of generating fading samples with a given correlation matrix [60], [61,Sec. 2.2]. The KL representation allows the channel from the BS to the users in the l-th cluster to be written as where G l ∈ C K×K is a zero mean and uncorrelated Rayleigh fading matrix, i.e., E{G l } = 0, and E{G H l G l } = I K , the identity matrix with K diagonal elements. Note that because of the latter property, when plugging (3) in (2), one always gets the same correlation regardless the particular realization of G l , showing that it only depends on Λ l and U l in the KL representation. Λ l ∈ C K×r l is a diagonal matrix containing the singular values of R l and U l ∈ C M ×r l is a matrix containing in its columns the singular vectors of R l , where r k is the rank of the correlation matrix. As the diagonal matrix Λ l can be reduced to a r l × r l matrix, one has G l an K × r l matrix, and U l a M × r l matrix. Hence, the k-th user in the l-th cluster will experience the channel Note that in the case of uncorrelated fading the K terminals would see a channel h k,l ∼ CN (0, I M ).

It is assumed that the BS knows both H (LOS)
l and the autocorrelation in (2), the latter representing the partial knowledge held about (3). Therefore, while Λ 1 2 l and U H l can be estimated via the SVD (2), the fading row vector g T k,l is unknown at the BS and cannot be taken in consideration by a precoder. Moreover, because single-antenna terminals are considered, each terminal cannot mitigate this multi-path effect.
The energy harvesting circuits are typically non-linear and the effectively harvested energy will be a fraction of the RF energy available. A rich literature modeling the EH circuits exists, e.g., [31], [32], [41], [43], and it has been known that different types of RF-DC conversion circuits lead to different optimal waveforms. To assess the proposed precoding techniques, this paper adopts one model that still mimics important features such as the sensitivity and saturation phenomena. The non-linear EH circuit at the terminals is considered to have a saturation value, 2 , and also a minimum power sensitivity, defined by an activation threshold 1 , and a conversion factor η in its linear region, leading to the following output power, as in [8]:

III. OPTIMAL UNITARY PRECODERS FOR WET
A natural goal for WET is the one of maximizing the RF energy available to be harvested by the terminals in each cluster while avoiding energy leakage to other locations. An effective way of achieving that is to apply MIMO precoding at the BS. Using simple analogue phase modulators one can consider signals from a continuous-domain alphabet, with x l (m) = √ P x e jθm , m = 1, . . . , M , only having to adjust the phases θ m . Given that P l is unitary, with A. Precoding for a single cluster: maximum harvested energy The problem of maximizing the mean sum-power, P l , available for EH by all devices in the l-th cluster, given ||x l || 2 = ||x 0 || 2 = P x , can be formalized as: where the total available energy E l at the cluster is the sum of the energies available for harvesting at the K terminals, each of which can collect the power |y l (k)| 2 . By using (1) and (3), the total RF energy available for harvesting by the terminals in cluster l, and that should be maximized in (6), is and, for convenience, let one define: Consequently, the expectation over the realizations of the Rayleigh fading component G l becomes Unsurprisingly, this last expression does not depend on G l , given that the a fast fading unit-power MIMO channel, on average, is a norm-preserving linear transformation. The equivalent optimization problem becomes the one of simultaneously maximizing two squared norms, while considering the power uniformly distributed across the M antennas: x l will be used henceforth to indirectly define P l .

B. Precoding for a single cluster: correlated Rayleigh fading
When κ = 0, the optimal unitary precoder for the second term in (10) can be analytically constructed by splitting the precoding process in two, such that P l = P l is a diagonal matrix with decreasing elements, the maximum singular value, σ max , is always the first element in the diagonal and therefore the maximization requires P concluding that when there is no LOS and only the multi-path (MP) channel component exists, the optimal transmission is defined by the first column of U l , obtained in (2): using MATLAB ® notation. Implicitly, P (2) l applies a unitary rotation to x 0 to generate [1, 0 . . . 0] T , such that P (2) l is a rotation matrix, and therefore the total precoding matrix P l = P

C. Precoding for a single cluster: correlated Rician fading
In practical WET setups, the powering BS is usually located close to the clusters of users and therefore having a strong LOS component with κ > 1 is a likely situation. When both terms in (10) contribute to E{P l } = E{||y l || 2 }, the optimum x * l is obtained by maximizing all the coordinates y l (k), which, according to (10), corresponds to maximizing K inner products of x l with the k-th row of H LOS l , denoted by (h (LOS) k ) T and also other r l inner products in the norm of the second term of (10). Defining and m T k being its k-th row, the total power collected by the cluster is (for simplicity the l index is dropped in the vectors) which can be written in the form Both the LOS and the multi-path components contribute to the problem, such that A l ∈ C 2K×M is the maximum dimension of this matrix (depending on the rank of H (LOS) l and the rank of Λ l , r l ). From (14), one can conclude that the maximum possible harvested energy corresponds to the largest singular vector of A l , σ max , and is obtained when the transmitted signal, x l , is the (unit power) singular vector u (A l ) max associated to σ max . When taking in consideration both the LOS channel and the second order statistics of the multi-path channel, the optimal solution to (10) is therefore which is a vector living on the M -sphere of norm √ P x , and P (2) l is implicitly a rotation matrix acting on x 0 , as in the previous subsection.

D. Precoding for multiple clusters
The maximization of the average energy available to be harvested by the system's devices in a multi-cluster system, E{P }, can be attained by constructing the stacked vector y = [y 1 . . . y l . . . y L ] T with all the signals at all the terminals in all the clusters. The problem then becomes the one of maximizing the (squared) norm of this KL-dimensional vector. With a similar manipulation to the one in (13), one can obtain the optimal precoder for the multi-cluster case (L ≥ 1), which is the same of finding the optimal precoding vector x * = P * x 0 . Hence, is maximized when transmitting the right-singular vector associated to the maximum singular value of the SVD of the concatenated matrix A ∈ C 2KL×M (assuming that both H (LOS) l and Λ l are full-rank), which is also a vector living on the M -sphere of norm √ P x : The solution (17) takes into account the whole multi-cluster system, and for that reason it depends on all the correlation matrices R l and on all the LOS channels H (LOS) l . The complexity of the proposed precoding scheme is chiefly determined by the SVD of the autocorrelation and the SVD of A l ∈ C 2KL,×M , which has a time complexity O 4KLM 2 + M 3 + M + 2KLM and a storage complexity O 3M 2 + 3M + 4KLM by using the so-called truncated SVD algorithm [62].
Remark: (16) is applicable to a MU-MISO setup in which device clustering is possible, which implies a MU-MIMO setup with more than one EH terminal. When clustering is not possible in a MU-MIMO system [18], or in in the case of a point-to-point single-user (SU) MIMO [6], [7], the solution to the problem of maximizing the available RF power is also a SVD problem involving the corresponding Wishart matrix. In the case of the Rician model in (1), the distribution of the maximum eigenvector is known [63], as well for the case of a Rayleigh MIMO channel [64]. Moreover, in the simpler case of SU-MISO, the matrix problem collapses to the well-known vector problem of coherent combining [9].

E. Constrained precoding for efficient non-linear EH and inter-cluster fairness
Although the energy available for a cluster is E{||y l || 2 }, the sensitivity and saturation thresholds of the non-linear transfer function of the harvesting circuit, described in (5), induce energy waste unless 2 1 ≤ E{|y l (k)| 2 } ≤ 2 2 , at each terminal in the l-th cluster. In order to minimize that, the precoder should also take into account those additional constraints. An additional affect of adding the lower bound constraint for each of the KL terminals is the induction of fairness among the partial sum-power of each cluster. This is a natural consequence of the fact that the optimization is blind to which cluster each terminal belongs to. The constraint optimization problem is set in the following, where P is no longer required to be unitary because a solution for P spat out by the numerical solver may observe the set of constraints without having to reach the maximum allowed transmit power. One can simply state the problem as: This is a non-linear programming problem with non-linear inequality constraints on each component which can be numerically solved by converting the problem from the complex domain to the real domain. To that end, the MIMO realequivalent model, well-known context of MIMO detection [65], is applied to the objective function in (18), so that a realvalued solution is found by casting the problem in a similar manner to (16) and then applying the real-equivalent model to In the end, the solution is converted back to a x ∈ C M .
The constraint added to the problem in (18) requires setting an expression for the average energy at each terminal since its instantaneous value cannot be estimated by the BS given that it does not have access to the instantaneous realizations of G. The constraints on the average energy received by each terminal, E{|y l (k)| 2 } can only take into account expectations in the context of the KL decomposition channel model. Using (1) and (3), while g T k is the k-th row of G, and dropping the cluster index in the rows of H (LOS) , the per-terminal average energy, when given some precoding matrix P, is where v l M H l Px 0 = M H l x, which is a factor only depending on the cluster's correlation R l , while the fast fading g i plays no role in determining the average power. In (20), it was also used the fact that the variance of an inner product between a random vector a and a deterministic vector b observes var(a T b) = b T var(a)b.
In the case when full CSIT exists (which will be considered for benchmarking in the next section), because the instantaneous value of G is known, the exact y l is also known, as defined in (1), the expectation E{|y l (k)| 2 in (18) can be replaced by the exact received power at each terminal so that the second constraint becomes 1 ≤ |y l (k)| 2 ≤ 2 , for k = 1, . . . , K, and l = 1, . . . , L.

IV. THE BASELINE REFERENCE CASE
The previous section presented derived the beamforming vectors x l , including the two numerical precoders that solve the energy-constrained problem (with partial-and full-CSIT). In the follow up section V the performance of these precoders will be compared with the one attained by the AA and SA CSIT-free schemes. While the latter schemes can have an analytical interpretation [8], the former are analytically intractable. However, there is a particular situation that has an analytical interpretation: the one with a pure Rayleigh channel, i.e., without LOS (κ = 0), and where the fading is uncorrelated. This situation is considered as the baseline reference case: the i.i.d. Rayleigh fading channel is represented by a MIMO matrix with circularly symmetric complex Gaussian entries with zero mean and unit variance. In this situation h t,l ∼ CN (0, I M ), with both real and imaginary components of each element h t,l (m) having each a Gaussian distribution N 0, 1 2 , such that h t,l (m) ∼ CN (0, 1), for m = 1, . . . , M . This setup has an analytical interpretation chiefly based on the properties of the Gamma distribution [66], [67].
Definition 1: A Gamma random variable with finite shape parameter k > 0 and finite scale parameter θ > 0, denoted as Γ(k, θ), has PDF with mean kθ and variance kθ 2 . Hence, if the received power is described by a Gamma distribution, for a linear increase of the mean value there will exist a larger spread of the power domain (which can be also linear or become quadratic, depending whether it is the shape or the scale parameter that is changing).
The link budget is given by the equation P r = P x + G a + G P − β (where the gain of the ULA's antenna elements is considered to be G a dB, as in V), and therefore one can define an equivalent power gain and herein consider a system with unit total transmit power (P x = 1) and an ULA with antenna elements with no gain.  CN (0, 1), for t = 0, . . . , K.
Hence, when experiencing an i.i.d. Rayleigh channel, the partial sum-power available for the K terminals within the l-th cluster, is y l 2 ∼ Γ (k = K, θ = 1).
Proposition 4. If a random variable X holds a Gamma distribution Γ(k, θ), then, for a scalar a > 0, aX has a distribution Γ(k, aθ).
For this reason, a link to a cluster having an equivalent power gain β eq , the partial sum-power available to the K terminals is distributed according to Γ (K, β eq ) In order to consider the general situation where different clusters have different β l , one will need one more tool given by the so-called second order Gamma approximation of the sum of gamma functions [66] to be given in proposition 6. Note that if all beta l are the same, then by reapplying proposition 3, one could immediately state y l 2 ∼ Γ (KL, β eq ), given that all users could be considered in the same cluster.

Proposition 5.
[66] Suppose that {Y t } are independent Γ (k t , θ t ) random variables. The sum Y = t Y t has mean, variance, and second moment, respectively given by which is known as the second order Gamma approximation of the sum that builds up Y [66].
In the case where all the channels to the terminals in the l-th cluster are identically distributed, with h T t,l ∼ CN (0, β l I M ), then the partial sum-powers at each cluster are Gammadistributed Γ (k l = K, θ l = β eq ), for t = l, . . . , L, and (24) particularizes to leading to the following PDF for the system's sum-power: Proposition 7.
[67] For a K-dimensional vector with y 2 ∼ Γ (k y , θ y ), its projection onto c coordinates holds a power distribution Γ( c K k y , θ y ), meaning that its mean is c K kθ and its variance is c K k y θ 2 y , and therefore has a power given by its second moment: P = c K k y θ 2 y + c K k y θ y 2 , which amounts to the power available to a subset of c terminals, when k y and θ y are the ones in (26).
The last proposition establishes the partial sum-power available for EH by a subset of devices in a given cluster.

A. SA scheme
The second order Gamma approximation method can also be used to find the expression for the PDF of the system's sum-power when using the SA technique, providing a simpler proof than the one given in [25].
Under the SA scheme, the received signal at the t-th terminal results from a time-slot based accumulation of power, such that E{|y t | 2 } = It should be noted that by using SA, differently from (26), the shape parameter in (27) grows with M , linearly increasing the mean power, while the scale parameter is reduced by M , reducing the variance quadratically (as per definition 1 and proposition 5), which explains why the SA scheme is preferable to AA when κ = 0. V. NUMERICAL RESULTS In this section, one numerically assesses the proposed precoding (or beamforming) techniques based on the KL channel decomposition with limited CSIT, namely the correlation of the channels from the BS to each of the L clusters. Both the unconstrained optimal analytical precoding for multiple clusters, given by (17), and the numerically obtained constrained precoder given by (18), will be assessed and compared with the case when full-CSIT is available. The precoder when full-CSIT is available is the one that maximizes which is given by the singular value associated to the largest singular value of B, The primal method to to assess the proposed techniques is to focus on the PDF of the sum-power available for harvesting at the terminals, p(P ). Subsequently the analysis focuses on how the Rician factor, the dimension of the antenna array at the BS, the angular position of a cluster, and the number of clusters impact on the energy available for harvesting by the terminals when using the proposed optimal precoder defined in (15), for the single-cluster case, and in (17) for the multicluster case.

A. Sum-power distributions: setup
The first configuration to be assessed involved one single cluster only, initially with correlated Rayleigh fading only, assessing the precoder devised in section III-B, and then adding a LOS component, assessing the precoder devised in section III-C. These initial results are not shown due lack of space. The first results here presented consider a multi-cluster system. The observations with one cluster were in line with the ones presented for a multi-cluster system in terms of the impact of the system's parameters in the overall PDF of the sum-power. The set of distributions presented in Fig. 2 were obtained for a system with L = 3 clusters with K = 8 users per cluster and a BS with M = 8 antennas (hence, a WET system with a total of 24 terminals), which is still far from a massive MIMO setup but is a typical MIMO array. Without loss of generally, the mean pathloss to each cluster, β, was considered equal to all clusters, and a LOS characterized by κ = 5 was also considered equal in all clusters. Such Rician factor is adequate to a mMTC system as the one in Fig. 1. In [10] a κ = 10 was considered, accounting for systems with a quite large LOS component, while in [8] a κ = 3 was rather used. Fig. 2 shows the PDFs of the sum-powers obtained by all the techniques and the analytical curve for the baseline reference case is also added, as defined by (26) , with shape parameter KL = 24 and scale parameter β = 10 − 63.5 10 . Considering the experimental findings in [69], the present paper considers 1 = −22 dBm (6.30 µW) and 2 = −4.8 dBm (311 µW) to model the non-linear EH circuit, likewise the conversor in [8]. The mid-point of the linear region of this terminal's circuit amounts to 168.7 µW (−7.73 dBm).
The results shown in Fig. 2 are for a situation in which the pathloss, β, is such that the mean available power at each terminal is precisely 168.7 µW, which for a set of 24 terminals corresponds to a system's sum-power P = 4.05 mW (6.07 dBm). With the system's parameters summarized in Table I, the precoding gain achieved by the proposed optimal precoder can be calculated by using the link budget equation, such that: G P = 6.07 dBm − 40 dBm − 10 dB − β [dB] = 19.57 dB. The chosen pathloss set to β [dB] = −63.5dB corresponds to a ∼ 23 meters distance between the BC and the center of the clusters when using the LOS Urban Microcell 3GPPP propagation model at 1900 MHz [70], [71]), which is conceived for outdoor microcells, better matching the scenario expected for the proposed techniques, as sketched in Fig. 1 One should note that the distributions in Fig. 2 are for the power available at the terminals. The effective harvested energy will have to consider the non-linear effects of (5) and the conversion factor η. Given the symmetric shapes of the sum-power distributions when precoders are used, the mean collected power at the l-th cluster is ≈ ηP l . With η = 0.25 (−6 dB) [8], [69], this leads to a mean harvested sum-power of ≈ 1 mW by the 24 terminals.
Let us now analyze the limits of the available power obtained when using the unconstrained KL-based in the same analyzed in Fig. 2. Considering the approximation of a sum-power distributed between 3.5 mW and 4.75 mW, the incoming power at each one of the 24 terminals is distributed between 148.8 µW and 197.9 µW , which is well within the [−22, −4.8] dBm interval, which corresponds to a [6.3 331.1]µW interval.
Remark: this set of parameters used in this subsection will later be referred to as the operation point (OP) in Fig. 2, and is summarized in Table I.

B. Sum-power distributions: analysis
The first thing worth mentioning is that when there are no restrictions to the minimum and maximum values of E||y l (k)|| 2 (i.e., in the case of ideal harvesting circuits), the results attained by the analytically designed precoders (given by (15) or (17) and the precoders that result from the numerical solution of the linear programming problem with non-linear inequality constraints in (18) are exactly the same; they not only lead to the same average available sum-power, P , in a multi-cluster or single-cluster system, but also exhibit the same PDF for the sum-power. Interestingly, the particular solutions for the broadcast vectors x * coming from the two approaches are often different, and therefore the two precoders P are different. However, these two different vectors x * are equivalent in the sense that they lead to the same distributions for the same set of correlations R l (and therefore for the same A, defined in (16). This is a consequence of the non-convexity of the optimization problem, exhibiting several equivalent local minima.
A paramount observation is that the sum-power distribution over the whole system obtained by the proposed technique, based on statistical CSIT, leads to a quasi-optimal energy harvesting situation, exhibiting a negligible loss with respect to the harvested sum-power when full-CSIT is available and the optimal precoder is used. While the precoder based on the KL decomposition is quasi-optimal in respect to the full-CSIT case, by only setting the standalone goal of maximizing the sum-power available at the terminals may lead (for the KLbased or full-CSIT precoders) to a rather unbalanced distribution of the available power at the different clusters. In order to observe that effect, Fig. 2 also includes the distributions of the partial sum-powers at each of the three clusters, striking out that the average partial sum-power available at one of the clusters is one order of magnitude lower than the one available at the other two clusters. As seen in Fig. 2, two of the clusters hold very similar distributions, with mean values of 1.87 mW and 1.95 mW for their partial sum-powers, while the remaining cluster only gets a mean sum-power of 0.23 mW available at its terminals.
A fairness criterion for power allocation should take into account the effective power that the terminals harvest with the quasi-linear EH circuit characterized in (5). This was framed in (18), such that the numerically obtained precoder attains a fair distribution of energy among the clusters in terms of mean partial sum-power. Fig. 2 also shows that the proposed constrained precoder attains the sought fairness, however at the expense of a reduction of the system's total sum-power, exhibiting the typical trade-off that exists in optimization problems when some fairness criteria is applied [18], [73], [74]. Notably, the system's sum-power obtained with numerical constrained precoding incurs an affordable degradation of 4.13 − 3.77 = 0.36 mW (8.7%) in respect to the average attained in the full-CSIT situation (the best possible one).
For comparison purposes, two CSIT-free WET techniques that make use of linear arrays, namely the AA and SA schemes [8], were also simulated and the results are also plotted in Fig.  2. While the AA scheme provides a simple solution for spacial diversity, highly dependent on the angular deviation from the antenna's boresight, the SA is able to offer the same power for any angular position of a single cluster. When using the AA technique, the power harvested by the k-th terminal within the L-th cluster is while for the SA scheme that power is As previously stated, this paper focuses on the power available for EH, focusing on the RF power available at the antennas before considering the effect of Ω.

C. Validation with the baseline reference case
The simulation of the analytical and numerical precoders was validated with the reference case where no clusters are considered and the link from the BS to all EH devices is an uncorrelated Rayleigh channel, as defined in section IV. This can be considered as the baseline situation with which the mean sum-powers attained by the different techniques can be compared to. Such results are plotted in Fig. 3(a). It should be noted that the results with the AA scheme and with two of the precoded systems overlap the analytical curve defined by expression (26) for the baseline reference case (with KL = 8 terminals); these are i) the precoder derived for a channel with correlated Rayleigh and no LOS, as defined in (11), and ii) the precoder derived for the Rician channel with a correlated Rayleigh component (but where now κ = 0), as defined in (15). In this situation there is neither a notion of cluster nor of angle and the AA scheme performs as well as any of the precoded schemes. In this case, as also concluded in [8], [25], the SA scheme is preferable because it provides a slightly  larger mean sum-power and a significantly smaller variance of the PDF. A case where there still no LOS component but the Rayleigh fading is correlated is shown in Fig. 3(b): the optimal precoder for correlated fading, given by (11), already exhibits the considerable gain that one can extract from the KL-based approach. In this case the PDF obtained using the precoder in (15), which also takes in consideration the LOS component, overlaps the previous curve, as expected, given that κ = 0. However, when some LOS component emerges, as in the case in Fig. 3(b), then the precoder for Rayleigh correlated fading has a subpar performance with respect the one incorporating information about the LOS component; as seen in Fig. 3(c), the performance of the latter precoder detaches itself and an extra  gain kicks in due to the extra information about the channel associated to the LOS CSIT (this gain can also be later seen in Fig. 4 when looking at the κ = 0 abscissa). In this baseline case with uncorrelated fading in Fig. 3(a), while there is no gain to be attained by the KL-based analytical precoders, the constrained numerical version of the proposed precoder (in sec. III-E) narrows the sum-power distribution at the expense of a lower mean, as typical of fairness versus mean trade-offs. The same effect is narrowing the variance and reducing the mean also appears in the results of the constrained precoder in the situation of full-CSIT availability (the two curves with larger means in Fig. 3(a)).
Remark: the results presented in Fig. 2 and Fig. 3, as well as the ones presented in the subsequent subsections, are obtained with 10 5 simulated channels and are ergodic in the sense that for each channel instance a new fading correlation was generated.

D. Rician factor
The rich information provided by how the different PDFs vary when the system's parameters change motivates a more detailed study of how each of the system's parameters impact on the sum-power distributions, specifically how their mean and variance is modified. The impact of the LOS component will be the first to be analyzed in Fig. 4, while maintaining all the parameters on Table I unchanged.
When the LOS component grows larger in comparison to the Rayleigh fading component, the channel becomes more deterministic and for that reason the performance tends to the one achieved with full channel knowledge. One should notice in Fig. 4 that the OP scrutinized in Fig. 2 already corresponds to a regime in which CSIT leads to a very diminished increase in E(P ), unlike what happens in the range where κ < 1 (even though in that region the power spreads over a wider domain, as depicted in the curve for V ar[P ]). The advantage of the proposed precoding scheme stands out in comparison with the energy harvesting ability of both the SA and AA schemes. Noticeably, at the abscissa κ = 0 it is possible to see an over 3 dB gain. In this situation, without a LOS component, all the gain appears by virtue of the KL-based approach, by having access to the second orders statistics of the fading only. Then, an increasing weight of the LOS increased the mean sum-power while simultaneously diminishing its variance, given that the channel tends to become deterministic. Fig. 4 also shows the effect of considering the terminals uniformly scattered over an increasing angular aperture ∆ φ around the central angle of each cluster.

E. Array size
When the BS is fitted with a larger number of antennas, M , one can observe in Fig. 5(a) a quasi-linear growth of the system's available power for harvesting when the ULA's number of antennas is a power of two (as it is typical in commercial ULAs). The increase in E[P ] comes in tandem with a quasi-linear spread of the variance, meaning that higher gains are intertwined with a less constant sum-power level. This increasing focusing effect onto the clusters, translates to larger precoding gains for WET that are unattainable by AA (apart when aligned with the antenna's boresight) or even SA. For the particular set of three clusters tested in Fig. 5, an ULA with M = 128 antennas attains a 40 fold gain (16 dB) in respect to the AA scheme, and this is before considering any optimization of the ULA's rotation, as will be seen in the next subsection. It should be noted that, while in the case with a single cluster both the mean and the variance always exhibit a linear growth (not shown in the paper) for any step size of the incremental M , when L > 1 the variance drops bellow linear growth in a periodic manner when analyzing the results for an ULA growing one antenna at a time, as can be observed in Fig.  5(b). It has been found that this is due to the radiation pattern and the particular positions of the L clusters in respect to the radiation pattern.

F. Cluster's angular position
In order to analyze the effect of the angular cluster's positions, Fig. 6 depicts a situation of a single cluster positioned at varying angles. It is known (e.g., [8]) that the AA scheme can provide a significant gain due to constructive interference but limited to the ULA's boresight. The SA scheme can provide a constant mean sum-power for any angle, however at a much lower level than AA's peak. Fig. 6 shows that the KL-based precoder can sustain the AA's peak at all angles, and again, in the considered scenario one finds a negligible loss with respect to the mean sum-power attained with full-CSIT.
The angular positions of the clusters are relative to the geometry of the array, i.e., relative to the broadside and to the end-fire of the ULA. The l-th cluster is considered to be positioned at angle φ l in respect to the ULA's end-fire, such  that an array facing the boresight has and angular position φ l = 90 degs. If the linear array can be mechanically rotated by a certain angle ψ, then the effective new positions of a set of L = 3 clusters change to In Fig. 7 one assesses the effect of such rotation starting with φ = {0, 30, 70} degrees, the set previously used in Fig. 2. Given the single-lobe directivity of the AA technique (as seen in Fig. 6), one finds peaks for the sum-power when the new effective positions of each cluster is 90 degrees, which for φ = {0, 30, 70} degrees, happens for ψ = 20, 60, and 90 degrees. Notably, one finds that there is an optimal rotation that maximizes the system's sum-power for ψ = 144 degrees, attaining E{P } = 5.458 mW. For this particular positions of the clusters, the proposed unconstrained scheme attains the same sum-power mean and variance of the full-CSIT case, for a wide domain of rotations ψ. The need for optimizing the ULA's rotation in WET has also been noted in [10], [28].

G. Number of clusters
With a constant number of K devices per cluster, an increasing number of clusters L leads to a proportional increase of the number of terminals in the system and E{P } increases as a direct consequence of having more harvesting terminals in the system. The simulation results consider clusters positioned at the angles φ l = l L + 1 π, for l = 1, 2, . . . , L, as depicted in Fig. 8 up to eight clusters, and the corresponding results are plotted in Fig. 9. This choice of angular symmetries leads to an apparent non-monotonicity of the harvested energy in the case of AA system. This is illusory given that, for an odd number of clusters, (33) places a central cluster at the boresight of the ULA, which, as seen in the previous subsection in Fig   6, is the direction of a narrow beam that chiefly powers that central cluster. In fact, both the maxima (at odd L) and the minima (at even L) also increase with L, as more clusters are packed together and the ones neighboring the central cluster increasingly get more illuminated by the AA's beam. The values obtained with the previous setup, with three clusters located at φ = {0, 30, 70} degrees, are also superimposed on Fig. 8 (at abscissa L = 3 ), and one can observe that that set of angles leads to a larger mean available power, with larger variance, then the set stipulated by (33).

VI. DISCUSSION AND CONCLUSIONS
In the context of mMTC, WET is a solution for BS wirelessly power a very large number of EH devices. However, in many practical situations the terminals may appear clustered (due to physical proximity), and their channels to the BS may have several parameters in common. Considering Rician channels from the BS to the terminals, the terminals within defined clusters were assumed to share a common correlation of their Rayleigh fading components, as well as having a slowly varying LOS channel component (possibly different but fully-known at the BS) and sharing the same largescale pathloss. Rather than considering the partial channel knowledge of the non-LOS fading component as the sum of a known and an unknown term, this paper considers a statistical representation of that random component. The proposal takes advantage of the Karhunen-Loève representation of a MIMO channel, such that the correlated Rayleigh component can be written as the concatenation of i) an unknown fastchanging Rayleigh uncorrelated component and ii) CSIT only depending on the first and second order statistics of the channel (mean pathloss and correlation, respectively) for the coherence duration of the fading's autocorrelation. In doing so, the paper proposed beamforming precoders for multiantenna BS in WET systems operating with low-rate statistical CSIT. These precoding schemes require low computational complexity, chiefly depending on SVD calculations for each one of the clusters.
The results exhibit quasi-optimal performance under Rician fading when compared with the ones obtained when beamforming with full-CSIT knowledge, even for Rician factors in the interval κ ∈ [1 5] (i.e., below the one considered in [10]), and, in particular, a positive gain exists even for κ = 0, by just leveraging the existence of correlated fading.
Because the number of clusters is bound by the maximum number of beams, one has L ≤ M and therefore the use of a massive array is highly desirable. The precoding gain grows linearly with the number of antennas at the BS (at the expense of wider power spread) and therefore the desideratum of wirelessly power a large number of terminals without instantaneous CSIT can be fulfilled. The same linear increase of the mean sum-power also exists when using a precoder with access to instantaneous full-CSIT. The available CSITfree multi-antenna schemes for WET have a much poorer performance in similar conditions and are not able to take full advantage of a larger array, nor to provide multiple energy beams to different sites.
By setting as the only goal the maximization of the sumpower in a multi-cluster, an unfair power allocation between the different clusters may emerge. To withstand that effect, a constrained optimization of the precoders is proposed which, besides naturally enforcing inter-cluster fairness, can also conform the power domain at the terminals' antennas to the linear domain of their non-linear EH circuit. For the chosen setup with three clusters with 24 terminals in total, with a pathloss of 63.5 dB from the BS to each of the clusters, corresponding to ∼ 23 meters at 1900 MHz (by lowering the frequency the range will be larger), neither the saturation point of the EH non-linear circuit considered was ever reached, nor the power available per terminal fell below the EH sensitivity. In fact, the available power at each terminal's antenna is always within the limits of the linear domain of the EH circuit.
In environments with a reduced number of propagation paths or with highly correlated paths, one could reduce the dimensionality of the problem by neglecting some of the singular values of the autocorrelation matrix, similarly to the idea in [75], which made use of the related Karhunen-Loève transformation [76,App. E].
Similarly to the conclusions in [10], where a mechanical rotation of the ULA is found to be a parameter to be optimized, the present research has shown that, in the case of multi-cluster systems, a similar mechanical rotation should also be optimized such that the clusters are placed in a set of more favorable angles. Interestingly, it was recently found in a context of LOS channels that an ULA can attain a performance very close to capacity with proper angular rotation, even without the need of a mechanical rotation, but rather by selecting a particular array out of a very low number of fixed ULAs with different relative rotations [77].
As an immediate future work, it would be interesting to combine the herein proposed beamforming optimization with the optimization of the waveform [13], [20], [78] and analyze the reachable energy savings at the BS for the same EH constraints at the devices. This research can also be extended to the realm of SWIPT, with modulated waveforms [79]. While statistical information can be more accurately obtained, the impact of working with imperfect statistical CSIT should also be evaluated. Additionally, even though the proposed signal processing techniques for precoding are blind to the geometry of the antenna arrays at the BS, the linear arrays considered in this work only permit a 2D tuning of the radiation pattern. Hence, for certain deployment where the terminals are also spread in height (e.g., in some types of infrastructures monitored by sensors) it will be important to extend the techniques to MIMO beamforming with 3D steering capability [80]. The constrained optimization problem can be generalized to the situation where the terminals have different sensitivity and saturation points. A likely scenario could be the one of having different types of equipment associated to different clusters, for example, monitoring different types of infrastructure.