Optimal Transmit Strategy for Multi-user MIMO WPT Systems With Non-linear Energy Harvesters

In this paper, we study multi-user multi-antenna wireless power transfer (WPT) systems, where each antenna at the energy harvesting (EH) nodes is connected to a dedicated non-linear rectifier. We propose an optimal transmit strategy which maximizes a weighted sum of the average harvested powers at the EH nodes under a constraint on the power budget of the transmitter. First, for multiple-input single-output (MISO) WPT systems, we show that it is optimal to transmit scalar symbols with an arbitrary phase and an amplitude, whose probability density function (pdf) has at most two mass points, using maximum ratio transmission (MRT) beamforming. Then, we prove that for single-input multiple-output (SIMO) WPT systems, the optimal transmit symbol amplitudes are discrete random variables, whose pdf also has no more than two mass points. For general multi-user MIMO WPT systems, we show that the optimal transmit strategy involves scalar unit-norm symbols with arbitrary phase and at most two beamforming vectors. In order to determine these vectors, we formulate a non-convex optimization problem and obtain an optimal solution based on monotonic optimization. Since the computational complexity of the optimal solution is high, we also propose a low-complexity iterative algorithm to obtain a suboptimal solution, which achieves nearly optimal performance. Our simulation results reveal that the proposed transmit strategy for multi-user MIMO WPT systems outperforms two baseline schemes, which are based on a linear EH model and a single beamforming vector, respectively. For a given transmit power budget, we show that the harvested power saturates when increasing the number of transmit antennas. Finally, we observe that the harvested power region spanned by multiple EH nodes is convex and the power harvested at one EH node can be traded for a higher harvested power at the other nodes.


I. INTRODUCTION
The device density in wireless communication networks has significantly increased over the past decades. The current trends for wireless systems suggest that the number of connected devices will continue to grow over the next few years and a longer battery life for these devices is highly desirable [3]. However, efficient charging of the batteries of wireless devices remains an unsolved problem. Since radio frequency (RF) signals are capable of transferring power, in recent years, far-field wireless power transfer (WPT) has attracted significant attention [4]- [17].
In [4], the authors studied single-input single-output (SISO) WPT systems and showed that the power transferred to the energy harvester (EH) is maximized if a single sinusoidal signal is broadcasted by the transmitter (TX). The authors of [5] extended these results to multiple-input multiple-output (MIMO) WPT systems and showed that the input power at the EH is maximized if a scalar input symbol and energy beamforming, i.e., beamforming in the direction of the dominant eigenvector of the channel matrix, are employed at the TX. Although the solutions developed in [4] and [5] are optimal for the maximization of the power received by the EH, they do not necessarily maximize the harvested power since practical EH circuits are non-linear [6]- [8]. Hence, an accurate modeling of the EH circuit is crucial for the design of WPT systems [6]- [17].
Practical EHs typically employ a rectenna, i.e., an antenna followed by a rectifier circuit that includes a non-linear element, namely, a diode. The experimental results in [7] and [8] showed that rectennas exhibit a non-linear behavior in both the low and high input power regimes. In particular, for low input power levels, the rectifier non-linearity is caused by the non-linear forward bias current-voltage characteristic of the diode [18], whereas in the high input power regime, practical EH circuits suffer from saturation due to the breakdown effect of the diode [19]. In order to capture these non-linearities, the authors in [7] modeled the harvested power as a parameterized sigmoidal function of the received power whose parameters depend on the waveform of the received signal. The model in [7] is widely exploited for the design of WPT systems, e.g., [9], [10], where Gaussian distributed input signals are assumed. In particular, in [9], it was shown that, adopting the EH model from [7], the energy beamforming proposed in [5] is also optimal for MIMO WPT systems. Furthermore, to avoid saturation of the rectenna circuits, based on the model in [7], the authors in [10] proposed to split the RF power received at the EH node between several collocated rectifiers.
Although the model in [7] characterizes the non-linear behavior of rectenna circuits, it is applicable only for signals with a known fixed waveform and does not allow the optimization of the waveform of the transmit signal [6]. Therefore, the authors in [11] proposed a non-linear EH model derived from the Taylor series expansion of the current flow through the rectifying diode of the EH. Based on this model, the authors in [12] studied a WPT system with multiple antennas at the TX and a single antenna at the EH node, i.e., a multiple-input single-output (MISO) WPT system, and showed that the harvested power is maximized with energy beamforming [5], which reduces to scaled maximum ratio transmission (MRT) in this case. However, in [13], it was shown that energy beamforming is not optimal for general MIMO WPT systems employing non-linear rectenna circuits. In [14], the authors considered a MIMO WPT system, where the EH node is equipped with a single rectifier and proposed to combine the received signals of different antennas in the RF domain to increase the power at the input of the rectifier and, thus, maximize the harvested power. However, the practical implementation of the RF combining schemes considered in [10] and [14] requires complicated circuit designs and may also introduce associated losses, which are not desirable in practical systems [20]. Finally, the authors in [15] considered a MIMO WPT system, where each antenna of the EH node was equipped with a dedicated rectifier, and proposed an iterative algorithm to determine the TX beamforming vector that maximizes the weighted sum of powers harvested by the rectifiers.
Although the results in [11]- [15] provide important insights for the design of practical MIMO WPT systems, their applicability is limited to low input power levels at the EH since the saturation of the harvested power is neglected in the underlying EH model [11]. A realistic EH model that accurately captures the rectenna non-linearity for both low and high input powers was developed in [16]. The analysis in [16] showed that for SISO WPT systems, it is optimal to adopt ON-OFF signaling at the TX, where the ON symbol and its probability are chosen to maximize the harvested power without saturating the EH while satisfying an average power constraint at the TX. The optimality of ON-OFF signaling was confirmed in [17], where a learning-based approach was employed to model non-linear rectenna circuits equipped with a single and multiple diodes, respectively [19]. Finally, in [1] and [2], which are the conference versions of this paper, exploiting the rectenna model derived in [16], the authors studied the harvested power region of a two-user MISO WPT system and the maximum performance of a single-user MIMO WPT system, respectively. However, to the best of the authors' knowledge, the problem of optimizing the transmit strategy for multi-user MIMO WPT systems, where the EH nodes are equipped with multiple rectennas exhibiting non-linear behavior in both the low and high input power regimes, has not been solved, yet. We note that the EH model adopted in [1], [2] is a special case of the more general EH model considered in this paper. Hence, the results obtained for two-user MISO and single-user MIMO WPT systems in [1] and [2], respectively, are special cases of the results presented in this paper.
In this paper, we aim at determining the optimal transmit strategy for multi-user MIMO WPT systems, where multiple EH nodes are equipped with multiple non-linear rectennas. In order to take the non-linearity of the EH into account, we consider a general rectenna model characterized by a set of properties, which are typically satisfied for the practical rectenna circuits considered in the literature [5]- [17]. Where appropriate, we specialize our results to the non-linear EH model proposed in [16]. We propose an optimal transmit strategy, which is characterized by the distribution of the transmit symbol vector that maximizes a weighted sum of the average harvested powers at the EH nodes subject to a constraint on the power budget of the TX. The main contributions of this paper can be summarized as follows: • For MISO WPT systems, we show that the optimal transmit strategy comprises MRT beamforming and a scalar input symbol with an arbitrary phase and a discrete random amplitude, whose probability density function (pdf) has at most two mass points. The optimal pdf of the symbol amplitudes is the solution of an optimization problem, which is solved via a two-dimensional grid search [21]. Furthermore, for the EH model in [16], we obtain a closed-form solution and show that the optimal pdf reduces to ON-OFF signaling.
• For multi-user SIMO WPT systems, we show that the optimal transmit symbol amplitude is a discrete random variable, whose distribution has at most two mass points and can also be obtained by a two-dimensional grid search. Similar to the MISO case, we show that for SIMO WPT systems with two rectennas modeled as in [16], the optimal distribution can be obtained in closed-form and ON-OFF signaling is optimal if the power budget of the TX is low. Furthermore, if it is affordable by the average power constraint of the TX, the discrete transmit symbol amplitudes are chosen to drive one or both rectifiers into saturation, respectively. Finally, for high average power budgets, the optimal policy is to saturate both rectifiers and the optimal pdf consists of a single mass point.
• For general MIMO WPT systems, we show that the optimal transmit strategy employs a scalar input symbol and at most two beamforming vectors, which can be determined as solution of a non-convex optimization problem. The optimal solution of this problem is obtained via monotonic optimization [22].
• To reduce the high computational complexity of determining the optimal beamforming vectors, we develop a low-complexity iterative algorithm based on semi-definite relaxation (SDR) and successive convex approximation (SCA) to obtain a suboptimal solution. Our simulation results reveal that although the suboptimal solution for MIMO WPT systems has a much lower computational complexity than the optimal one, both solutions yield a similar performance.
• Our simulations show that the proposed MIMO WPT design outperforms two baseline schemes, one based on the linear EH model in [5] and the other based on a single beamforming vector at the TX. For the multi-user scenario and a given transmit power budget, we observe that the total average harvested power saturates when the TX is equipped with a large number of antennas. Finally, we observe that the harvested power region spanned by multiple EH nodes is convex and the average power harvested at one EH node can be traded for a higher harvested power at the other nodes.
The remainder of this paper is organized as follows. In Section II, we introduce the system model and discuss the adopted EH model. In Section III, we formulate an optimization problem for the maximization of the weighted sum of the average harvested powers at the EH nodes and establish a preliminary mathematical result needed for solving the problem. In Section IV, we determine the optimal transmit strategies for MISO, multi-user SIMO, and multi-user MIMO WPT systems, respectively. In Section V, we provide numerical results to evaluate the performance of the proposed designs. Finally, in Section VI, we draw some conclusions.
Notation: Bold upper case letters X represent matrices and X i,j denotes the element of X in row i and column j. Bold lower case letters x stand for vectors and x i is the i th element of x.
X H , Tr{X}, and rank{X} denote the Hermitian, trace, and rank of matrix X, respectively. The expectation with respect to random variable x is denoted by E x {·}. The real part of a complex number is denoted by {·}.
x and x 2 represent the transpose and L2-norm of x, respectively.
The imaginary unit is denoted by j. The sets of real, real non-negative, and complex numbers are denoted by R, R + , and C, respectively. 1 K and 0 K represent column vectors comprising K elements, where all elements are equal to 1 and 0, respectively. The Dirac delta function is denoted by δ(x). f (x 0 ) denotes the first-order derivative of function f (x) evaluated at point

II. SYSTEM MODEL AND PRELIMINARIES
In this section, we present the MIMO WPT system model and discuss the adopted multi-antenna EH model.

A. System Model
We consider a narrow-band multi-user MIMO WPT system comprising a TX with N T ≥ 1 antennas and M EH nodes, where EH node m, m ∈ {1, 2, . . . , M }, is equipped with N E m ≥ 1 antennas, see Fig. 1. The TX broadcasts a pulse-modulated RF signal, whose equivalent complex baseband (ECB) representation is modeled as the transmitted vector in time slot n, ψ(t) is the transmit pulse with rectangular shape, and T is the symbol duration. Transmit vectors x[n] are mutually independent realizations of a random vector x, whose pdf is denoted by p x (x).
The ECB channel between the TX and antenna p of EH node m is characterized by row-vector g m p ∈ C 1×N T , p ∈ {1, 2, . . . , N E m }. Thus, the RF signal received in time slot n at antenna p of EH node m is given by z RF p,m (t) = √ 2 {g m p x(t) exp(j2πf c t)}, where f c denotes the carrier frequency. The noise received at the EH nodes is ignored since its contribution to the harvested energy is negligible.

B. Energy Harvester Model
In this paper, we assume that EH node m is equipped with N E m rectennas, i.e., each antenna is connected to a dedicated rectifier, see Fig. 1. Thus, the received ECB signal 1 at rectenna p ∈ {1, 2, . . . , N E m } of EH node m in time slot n is given by z m p [n] = g m p x[n]. Each rectenna comprises an antenna, a matching circuit, a non-linear rectifier with a low-pass filter, and a load resistor [6], [16], [17]. In order to maximize the power transferred to the rectifier, 1 We note that the EH nodes do not convert the RF signal into baseband. However, since the amount of the harvested power can be expressed as a function of the ECB signal received at the EH node, see, e.g., [23], [12], to simplify the notation, we use the ECB signal representation of the received RF signal. the matching circuit is typically well-tuned to the carrier frequency f c and is designed to match the input impedance of the non-linear rectifier circuit with the output impedance of the antenna [24]. The rectifier is an electrical circuit that comprises a non-linear diode and a low-pass filter to convert the RF signal z RF p,m (t) received by rectenna p of EH m to a direct current (DC) signal at the load resistor R L of the rectenna.
In this paper, we make the following assumptions concerning the rectenna circuit.   Assumption 1 is justified if the symbol duration T is sufficiently large. In this case, we can neglect the ripples of the voltage level across the load resistor R L and the charging and discharging times of the reactive elements of the circuit. Thus, we can assume that the output voltage level of rectenna p of EH node m in time slot n is constant and depends only on the signal z m p [n] [16], [17], [25]. Assumption 2 is justified since, for the considered narrow-band signals, the rectenna circuit behaves as an envelope detector [18] and, thus, its behavior is fully characterized by the magnitude |z m p [n]| of the received ECB signal z m p [n]. Assumption 3 is satisfied since typical rectenna circuits include a diode that has a non-linear non-decreasing current-voltage characteristic [16], [19], [25].
Example: In this paper, as an example for an EH model that satisfies Assumptions 1-3, we adopt the model proposed in [16]. The corresponding power harvested by the rectenna as a function of the magnitude of the received ECB signal z is given as follows: and W 0 (·) and I 0 (·) are the principal branch of the Lambert-W function and the modified Bessel function of the first kind and zero order, respectively. Here, Z * a , V T , I s , R s , and µ ∈ [1, 2] are parameters of the rectenna circuit, namely, the complex-conjugate of the input impedance of the rectifier circuit, the thermal voltage, the reverse bias saturation current, the series resistance, and the ideality factor of the diode, respectively. These parameters depend on the circuit elements and are independent of the received signal. Finally, since for large input power levels, rectenna circuits are driven into saturation [7], [16], [17], [19], the function in (1) is bounded, i.e., where A s is the minimum input signal magnitude level at which the output power starts to saturate.

III. PROBLEM FORMULATION AND USEFUL RESULT
In this section, we formulate an optimization problem for the maximization of the weighted average harvested power of the considered multi-user MIMO WPT system. Then, to obtain a preliminary result needed for solving this problem, we formulate and solve an auxiliary optimization problem, where we maximize the expected value of a function of a one-dimensional random variable under a constraint on its mean value.

A. Problem Formulation
We characterize the transmit strategy via the pdf p x (x) of transmit symbol vector x. The objective of the proposed transmit strategy is to maximize the weighted average power harvested at the EH nodes under an average power constraint at the TX. Thus, we formulate the following optimization problem: where the objective function is the weighted sum of the average harvested powers at the EH nodes defined as is the total power harvested by EH node m and ξ m ≥ 0, m ∈ {1, 2, . . . , M }, m ξ m = 1, is the weight for EH node m [14]. We note that the weights associated with the users allow the TX to control the distribution of the harvested power among the EH nodes. In particular, if weight ξ m is increased, the optimal transmit strategy will favor EH node m and increase the average harvested power E{ψ m (x)} at EH node m at the expense of the average harvested powers at the other EH nodes. Furthermore, we impose constraints (2b) and (2c) to limit the transmit power budget at the TX and ensure that p x (x) is a valid pdf, respectively.

Remark 1. Optimization problem (2) may have an infinite number of solutions. In particular,
for a general EH model satisfying Assumptions 1 -3, since x 2 and the average harvested power Φ(·) are invariant under phase rotation of the transmit symbol vector x, given an optimal pdf p * x (x) and a random phase φ x ∈ [−π, π) with an arbitrary pdf p φx (φ x ), the random vector x = exp(jφ x )x is still a solution of (2) [17]. Furthermore, we note that, for bounded φ(·), such as (1), if affordable by the power budget P x , there may be an infinite number of pdfs that drive all the rectifiers of the EH into saturation while satisfying constraints (2b) and (2c). Thus, in the following, we determine one pdf p * x (x) that solves (2).
As we will see in Section IV, the optimal solution of (2) leverages the solution of a related auxiliary optimization problem. In the next subsection, we solve this auxiliary problem, namely, the maximization of the expectation of a non-decreasing function f (ν) of a scalar random variable ν under a constraint on the mean value of ν.

B. Auxiliary Optimization Problem
Let us consider the following auxiliary optimization problem: whose solution is the pdf p ν (ν) which maximizes the expectation of f (ν) under a constraint on the mean value of ν. In order to solve (4), let us first define the slope of the straight line connecting points ν 1 , f (ν 1 ) and ν 2 , f (ν 2 ) , where ν 2 > ν 1 , as follows: Then, we establish an upper-bound on upper-bound on E ν {f (ν)} is given by the Edmundson-Madansky (Jensen's) inequality, e.g., [26].  where Here, ν * 1 and ν * 2 are given by ν * Proof. Please refer to Appendix A.
Lemma 1 is applicable for arbitrary non-decreasing functions f (ν). In the following corollary, we show that for a certain class of functions f (ν), the result in Lemma 1 can be significantly simplified.
Corollary 1. Let us consider a non-decreasing function f (ν) of random variable ν. If function f (·) is differentiable at ν = E ν {ν} and the following property holds: then the expectation of f (ν) is upper-bounded by E ν {f (ν)} ≤ f ν , where the inequality holds with equality if the pdf of ν is given by Proof. Please refer to Appendix B.
The results in Lemma 1, Corollary 1, and Corollary 2 will be exploited in Section IV for solving the transmit strategy optimization problem in (2).

IV. OPTIMAL TRANSMIT STRATEGIES
In this section, we first consider MISO and SIMO WPT systems, where the EH and the TX node are equipped with a single antenna, respectively. For each system architecture, we determine the pdf p * x (x) as solution of (2) and the resulting optimal transmit strategy. Then, we consider the general multi-user MIMO WPT case and present optimal and suboptimal solutions of (2).

A. MISO WPT Systems
In the following, we consider MISO WPT systems with a single-antenna EH node, i.e., M = N E 1 = 1, and a TX equipped with N T ≥ 1 antennas. In this case, the weighted sum in (3) , where g is the row-vector representing the channel between the TX and the EH node. In the following proposition, we provide a solution of optimization problem (2) for MISO WPT systems and the corresponding optimal transmit strategy.
Proof. Please refer to Appendix D.
Proposition 1 shows that, as for linear EHs in [5], for the considered non-linear EH model, MRT beamforming is optimal. Furthermore, similar to the SISO case in [16], for MISO WPT systems, there exists an optimal input symbol amplitude that follows a discrete pdf, p * rs (r s ), consisting of at most two mass points. In particular, it is optimal to adopt a single sinusoidal signal s with amplitude r s = √ P x and an arbitrary phase θ s if condition (8) holds. If (8) does not hold, amplitude r s is a discrete binary random variable. In this case, in order to obtain the pdf p * rs (r s ), the non-convex min-max optimization problem defined by (9), (10) has to be solved. Due to the low dimensionality of the problem, we propose to obtain the optimal solution via a two-dimensional grid search [21].
Finally, we obtain the power values ν * 1 = ρ i * and ν * 2 = ρ n+j * , where i * = arg min i max j S i,j and j * = arg max j S i * ,j , respectively. The proposed grid-search method is summarized in Algorithm 1.
The computational complexity of the proposed scheme is quadratic with respect to the grid size N ρ and does not depend on the number of antennas.
2) Special Case: In the following, we consider a special case of MISO WPT systems, where the EH model satisfies the following additional assumption.
1. Compute the grid P and the values of Φ(·) for the grid elements: In particular, we note that Assumption 4 holds for the EH model in (1). In the following corollary, we show that for an EH model satisfying Assumption 4, the result in Proposition 1 can be significantly simplified.
Corollary 3. For an EH model satisfying Assumption 4, the pdf p * rs (r s ) of the transmit symbol amplitudes r s in Proposition 1 is given by Proof. First, we note that ∀ν ∈ R + , Φ(ν) ≤ Φ(P max x ) = φ(A 2 s ). Therefore, for P x ≥ P max x , the optimal transmit strategy is characterized by the pdf p * rs (r s ) = δ(r s − P max x ) . Furthermore, functions φ(|z| 2 ) in (1) and Φ(ν) in Proposition 1 are convex and increasing in the intervals |z| 2 ∈ [0, A 2 s ] and ν ∈ [0, P max x ], respectively. Hence, ∀P x : P x < P max x , due to the convexity of Φ(ν), the solutions of the optimization problems in Proposition 1 are given by ν * 1 = 0 and ν * 2 = P max x , respectively, with β = P max This concludes the proof.
Corollary 3 reveals that if P x < P max x , ON-OFF signaling with MRT beamforming is optimal, which is similar to the result obtained for the SISO WPT systems in [16]. Furthermore, for x , it is affordable to drive the EH node into saturation and, hence, the optimal pdf of r s consists of a single mass point. We note that in contrast to Proposition 1, where all the power budget P x is utilized for the transmission, Corollary 3 shows that for EH models satisfying Assumption 4 and P x ≥ P max x , there is an optimal transmit strategy, where the average transmit power is equal to P max x .

B. SIMO WPT Systems
In the following, we consider a SIMO WPT system, where M EH nodes are equipped with N E m ≥ 1, m ∈ {1, 2, . . . , M }, antennas and the TX has a single antenna, i.e., N T = 1. In this case, as in [16], [17], due to Assumption 2, the powers harvested at the rectennas depend on the magnitude of the scalar transmit symbol but not on its phase. Hence, the weighted sum in (3) can be expressed as a function of the pdf p rx of the transmit symbol amplitude r x = |x| as and |g m p | is the magnitude of the scalar channel coefficient g m p between the transmit antenna and antenna p of EH node m. In the following proposition, we provide a solution of optimization problem (2) for SIMO WPT systems and the corresponding optimal transmit strategy.

Proposition 2.
For the considered SIMO WPT system, function Φ(·) is maximized for discrete transmit symbol amplitudes following distribution p * rx (r x ). In particular, the optimal transmit strategy is characterized by the pdf p * rx (r x = P x and the following inequality holds: Furthermore, if (12) does not hold, the optimal pdf is given by p Here, ν * 1 and ν * 2 are the corresponding solutions of the optimization problems in (9) and (10), respectively, where function Φ(·) is given by (11).
Proof. Please refer to Appendix E. Proposition 2 reveals that there exists an optimal pdf p * rx (r x ) of the symbol amplitudes r x that is discrete and consists of one or two mass points. In particular, as for SISO and MISO WPT systems in [4] and in Section IV-A, respectively, it is optimal to transmit a single sinusoid if condition (12) holds. If (12) does not hold, this optimal pdf consists of two mass points, ν * 1 and ν * 2 , which are obtained as solutions of the min-max optimization problem (9), (10). Due to its low dimensionality, this problem also can be efficiently solved via a two-dimensional grid search, as discussed in Section IV-A1 and summarized in Algorithm 1.
Special case: In the following, we consider a special case of the SIMO WPT system, where the EH model satisfies the following additional assumption.
Assumption 5. For function φ(·), Assumption 4 is satisfied. Furthermore, for the value A s , ∀z : |z| ≤ A s , the following inequality holds: We note that the condition in Assumption 5 implies that the power harvested by the rectifier grows slower than quadratically with the input power. In particular, it can be shown that Assumption 5 is satisfied for the EH model in (1). Let us now consider a SIMO WPT system with two rectifiers, i.e., M = 1 and N E 1 = 2 or M = 2 and N E 1 = N E 2 = 1. In this case, without loss of generality, we denote the scalar channel coefficients between the transmit antenna and the antennas of the rectifiers by g 1 and g 2 and assume that |g 1 | ≥ |g 2 |. Then, in the following corollary, we provide a closed-form solution for the optimal ν * 1 and ν * 2 in Proposition 2.
Corollary 4. Let us consider a SIMO WPT system with two rectifiers and an EH model satisfying Assumption 5. In this case, if P x < ρ min = A 2 s |g 1 | 2 , the optimal transmit strategy is characterized by the pdf of the transmit symbol amplitudes given by p Proof. Please refer to Appendix F.
Thus, for SIMO WPT systems with two rectifiers, if the transmit power budget is low, i.e., P x < ρ min , similar to the SISO and MISO WPT cases, ON-OFF signaling is optimal, where the ON signal drives the rectifier with the best channel conditions into saturation. Furthermore, for P x ∈ [ρ min , ρ max ), the pdf p * rx (r x ) has two mass points, which are chosen to drive one and both rectifiers into saturation, respectively. For P x ≥ ρ max , it is affordable to drive both rectifiers into saturation and, hence, the optimal pdf consists of a single mass point. Moreover, as for MISO WPT systems, Corollary 4 reveals that for P x ≥ ρ max , the average transmit power of the optimal transmit strategy is equal to ρ max .

C. MIMO WPT Systems
In the following, we consider the general multi-user MIMO WPT system in Fig. 1, where N T ≥ 1 and N E m ≥ 1 antennas are employed at the TX and EH node m, m ∈ {1, 2, . . . , M }, respectively. In the following proposition, we characterize a solution of (2) and the corresponding optimal transmit strategy.
Proposition 3. For multi-user MIMO WPT systems, function Φ(·) is maximized for discrete random transmit symbol vectors x = ws, where s = exp(jθ s ) is a unit-norm symbol with an arbitrary phase θ s . Here, w is a discrete random beamforming vector, whose pdf is given by . The beamforming vectors w * n , n ∈ {1, 2}, are given by where Ψ(x) = m ξ m ψ m (x). Here, ν * 1 and ν * 2 are the corresponding solutions of the optimization problems in (9) and (10), respectively, where Φ(·) is given by (14). Furthermore, if the following inequality holds: the optimal points ν * 1 and ν * 2 coincide and the optimal pdf is given by p Proof. Please refer to Appendix G.
Proposition 3 reveals that there is an optimal transmit vector x which is discrete and characterized by scalar unit-norm symbols s with an arbitrary phase 4 and at most two beamforming vectors, w * 1 and w * 2 . As in Corollary 1, these beamforming vectors coincide if inequality (15) holds. We note that w * 1 and w * 2 are characterized by the values ν * 1 and ν * 2 obtained as solutions of the non-convex problems (9) and (10), respectively, that also can be solved using the grid-search method in Algorithm 1. However, unlike for MISO and SIMO systems, in order to obtain the value of Φ(·) for a given transmit power value ν, the maximum value of Ψ(·) as solution of (14) is required. We note that (14) is a non-convex problem and, hence, obtaining its optimal solution is, in general, NP-hard. However, since problem (14) belongs to the class of monotonic optimization problems, in Section IV-C1, we first obtain the optimal solution exploiting the polyblock outer optimization approach [22]. Then, in Section IV-C2, for practical EH models satisfying Assumption 4, we propose an iterative low-complexity algorithm to obtain a suboptimal solution of the problem.

1) Optimal Solution:
In the following, we obtain the optimal solution of non-convex problem (14) exploiting monotonic optimization [22]. To this end, similar to the monotonic polyblock optimization framework in [22], [27], we obtain the optimal solution of (14) by exploring the feasible set V of (14) determined by constraint w 2 2 = ν. First, we enclose V by constructing 4 We note that the phase θs of scalar symbol s can be chosen arbitrarily in each time slot n. This degree of freedom can be further exploited, for example, for information transmission [17]. an initial polyblock B (1) with an initial set of vertices L (1) = ỹ (1) d |d ∈ {1, 2, . . . , 2 2N T } , whereỹ (1) d ∈ C N T is a vector, whose k th element is defined asỹ Here, n = 2 k and a d p ∈ {0, 1} denotes bit p in the binary representation of number d, see Fig. 3(a). Then, since the objective function in (14) is monotonically increasing in |w|, i.e., in |w 1 |, |w 2 |, . . . , |w N T |, in iteration m ≥ 1 of the proposed algorithm, as shown in Fig. 3(b), we  (14), we select the vertex χ * that maximizes the objective function Ψ(w). The proposed algorithm is summarized in Algorithm 2.
We note that the computational complexity of Algorithm 2 increases exponentially with the number of antennas N T employed at the TX. Therefore, obtaining the optimal value of function Φ(ν) may not be feasible in practical multi-user MIMO WPT systems. Nevertheless, the obtained optimal solution provides a performance upper-bound for any suboptimal scheme. In the next section, we propose an iterative low-complexity algorithm to obtain a suboptimal solution of (14).

2) Suboptimal Solution:
In the following, we consider a practical EH model φ(·) that satisfies Assumption 4. For this model, we propose an iterative low-complexity algorithm based on SDR and SCA to determine a suboptimal solution of (14). To this end, we first define matrix and u k is a unit vector containing only one non-zero element at position k.
Optimization problem (16) is non-convex due to the non-concavity of objective function (16a) and the non-convexity of constraint (16c). Therefore, in order to obtain a suboptimal solution of (16), we first eliminate constraint (16c). Then, we denote the total number of rectennas at the EH nodes by K = M m=1 N E m and define the sets W k , k ∈ {0, 1, . . ., K}, of matrices W such that ∀W ∈ W k exactly k rectifiers are driven into saturation. We note that W 1 ∪ W 2 ∪ · · · ∪ W K = S + .
Furthermore, rectifier p of EH node m is driven into saturation if and only if g m p W g m p H ≥ A 2 s . Hence, set W k , k ∈ {0, 1, . . ., K}, consists of T k = K! k!(K−k)! convex subsets, where k! denotes the factorial of k, i.e., W k = W 1 k ∪ W 2 k ∪ · · · ∪ W T k k . Each convex subset W t k , t ∈ {1, 2, . . ., T k }, consists of all matrices W which drive into saturation a specific combination of k rectennas. We note that the objective function in (16) is convex for each of these subsets and, hence, applying SCA for solving (16) for W ∈ W t k , is promising [28], [29]. Thus, the solution of (16) can be obtained by exploring the subsets W t k , t ∈ {1, 2, . . ., T k }, k ∈ {0, 1, . . ., K}, and solving the resulting problem for each subset [1]. However, since the computational complexity of this exploration grows with K!, in the following, we obtain a suboptimal solution of (16). For a given transmit power limit ν, we first determine a set of rectennas W * , which will be driven into saturation. To this end, we sort the channel gain vectors g m p in descending order of their norms as follows g m 1 p 1 2 ≥ g m 2 p 2 2 ≥ · · · ≥ g m K p K 2 , where m k ∈ {1, 2, . . ., M }, p k ∈ {1, 2, . . ., N E m k }, and k = 1, 2, . . ., K. Then, we check if it is possible to drive the k rectifiers with the best channel conditions, i.e., rectifier p 1 of EH node m 1 , rectifier p 2 of EH node m 2 , . . . , rectifier p k of EH node m k , into saturation by solving the following optimization problem: Optimization problem (17) is convex and can be solved with standard numerical optimization tools, such as CVX [30]. Furthermore, although we dropped the rank-one constraint (16c), it can be shown that if (17) is feasible and k > 0, a beamforming matrix W * k , which solves (17), has rank one. A corresponding proof is provided in the conference version [2,Appendix B] but is omitted here due to the space constraints. We denote by k * the maximum number of rectifiers k, for which problem (17) is feasible. Note that if (17) is not feasible for any k > 0, we have k * = 0. Then, we define the convex subset W * that corresponds to the case, where the k * rectifiers with the best channel conditions are driven into saturation. This set is given by Next, we reformulate problem (16) as follows: Optimization problem (19) is still non-convex due to the non-concavity of the objective function.
In the following, we propose to solve (19) exploiting SCA [28]. To this end, we construct an underestimate of the objective functionΨ(W ), which is convex in W * , as follows: Algorithm 3: Suboptimal algorithm for solving optimization problem (14) Initialize: Transmit power ν, tolerance error SCA .
for j = 1 to K + 1 do 3. Solve optimization problem (17) for k = j and store k * = j if the problem is feasible 5. Obtain w * as the dominant eigenvector of W (t) and evaluate Φ(ν) = Ψ(w * ) is the solution obtained in the iteration t of the algorithm and Ψ (W (t) ) denotes the gradient ofΨ(W ) evaluated at W (t) . Thus, in each iteration t of the proposed algorithm, we solve the following optimization problem: We note that (21) is a feasible convex optimization problem that can be solved with standard numerical optimization tools, such as CVX [30]. Furthermore, it can be shown that similarly to problem (17), the solution of (21) yields a matrix, whose rank is equal to one. Hence, we obtain the beamforming vector w * as the dominant eigenvector of the solution W * of (16) and compute the corresponding value of function Φ(ν) = Ψ(w * ). The proposed algorithm is summarized in Algorithm 3. We note that the proposed algorithm converges to a stationary point of (16) [29]. The computational complexity of a single iteration of the algorithm is given is the big-O notation.

V. NUMERICAL RESULTS
In this section, we evaluate the performance of the proposed transmit strategies via simulations.
First, we compare the performance of single-user MISO, SIMO, and MIMO WPT systems. Then, we evaluate the performance of multi-user MIMO WPT systems for different numbers of antennas 5 The computational complexity of a convex semidefinite problem that involves an n × n positive semidefinite matrix and m constraints is given by O √ n(mn 3 + m 2 n 2 + m 3 ) [31]. Here, n = N T and m = K + 1. at the TX and EH nodes. Finally, we study the influence of the weights ξ m , m ∈ {1, 2, . . . , M }, associated with the EH nodes and determine the harvested power region.

A. Simulation Parameters
In the following, we provide the system parameters adopted in our simulations. In our simulations, the path losses are calculated as 35.3 + 37.6 log 10 (d m ), where d m is the distance between the TX and EH node m [27]. Furthermore, in order to harvest a meaningful amount of power, we assume that the TX and each EH node have a line-of-sight link. Thus, the channel gains g m p follow Rician distributions with Rician factor 1 [32]. For the EH model φ(·), we adopt the model proposed in [16] and given by (1) with parameter values a = 1.29, B = 1.55 · 10 3 , I s = 5 µA, R L = 10 kΩ, and A 2 s = 25 µW. For Algorithms 1, 2, and 3 we adopt step size, grid size, and error tolerance values of ∆ ρ = 0.1, N ρ = 10 3 , PA = 10 −3 , and SCA = 10 −3 , respectively. We note that the grid size N ρ is chosen sufficiently large to ensure that the function Φ(ν) saturates for ν = ρ Nρ . All simulation results were averaged over 1000 channel realizations.

B. Single-user WPT Systems
In this section, we investigate the performance of single-user WPT systems with different numbers of antennas at the TX and the EH node. The distance between the TX and the EH is d = 10 m. The considered MISO and SIMO WPT systems employ N T = 2 and N E = 2 antennas at the TX and the EH node, respectively. For these systems, we evaluate the average harvested power Φ(p * x ) for different values of the power budget P x by applying Corollaries 3 and 4, respectively. In the considered single-user MIMO WPT system, the TX and EH nodes employ N T = N E 1 = 2 antennas, respectively. For this system, optimal and suboptimal pdfs p x (x) are determined by combining Algorithm 1 with Algorithms 2 and 3, respectively, and the values of Φ(p * x ) are computed. For comparison, we also consider a SISO WPT system employing the optimal transmit strategy in [16]. As Baseline Scheme 1, we consider a MIMO WPT system with energy beamforming at the TX, which is optimal for linear EHs [5]. Furthermore, as Baseline Scheme 2, similar to [14], we consider a MIMO WPT system, where a scalar input symbol and a single beamforming vector are adopted at the TX. For this system, we obtain the optimal beamforming vector w * as solution of (13) for ν = P x with Algorithm 2 and compute the corresponding harvested power as Φ = Φ(w * ).
In Fig. 4, we plot the average harvested powers Φ(p * x ) for different values of P x . First, we observe that for each considered WPT setup, the average harvested power Φ(·) is bounded above, since for the EH model in (1), for sufficiently large values of P x , all rectifiers of the EH node are driven into saturation. Furthermore, the saturation level of the harvested power is proportional to the number of rectennas employed at the EH. As expected, the MIMO WPT system achieves a superior performance compared to the SIMO and MISO WPT systems, which, in turn, outperform the SISO WPT system significantly. Interestingly, the MISO WPT system can harvest more power than the SIMO WPT system if the transmit power budget at the TX is low, whereas, for large values of P x , the single rectenna of the MISO WPT system is driven into saturation and, thus, more power can be harvested by the SIMO WPT system. Furthermore, we observe that for the MIMO WPT system, the proposed optimal transmit strategy, which employs two beamforming vectors, outperforms Baseline Schemes 1 and 2, which employ a single beamforming vector.
However, we note that Baseline Scheme 2, where the optimal beamforming vector obtained with Algorithm 2 is utilized, outperforms Baseline Scheme 1, where energy beamforming is adopted [5]. Finally, we note that although the computational complexity of the proposed suboptimal scheme for obtaining the MIMO WPT beamforming vectors is significantly lower than that of the optimal scheme, the performances attained by both schemes are practically identical. Therefore, in the next section, in order to keep the computational complexity low, we adopt the suboptimal scheme to evaluate the performance of multi-user MIMO WPT systems.  ξ m , we determine the optimal transmit strategies that maximize the corresponding weighted sum of the average powers harvested at the EH nodes. Then, for the optimal pdfs p * x (x), we evaluate and plot the average harvested powers at the individual EH nodes, respectively. For comparison, we also show the average harvested powers obtained with Baseline Scheme 1 and Baseline Scheme 2, respectively. In Fig. 6(a), we consider a low transmit power regime characterized by a power budget of P x = 10 W, whereas for the results in Fig. 6(b), we assume a high transmit power regime with P x = 30 W. The distances between the EH nodes and the TX are equal to d 1 = 10 m and d 2 = 25 m, respectively. In Fig. 6, we observe that higher values of N T and P x yield larger average harvested powers at both EH nodes. Furthermore, as expected, the proposed scheme yields a better performance compared to Baseline Scheme 1 and Baseline Scheme 2. For both transmit power regimes, we observe that for SIMO WPT, i.e., N T = 1, the performances obtained with both baseline schemes are identical and do not depend on the adopted weights ξ 1 and ξ 2 . In fact, in this case, the transmit strategies of the baseline schemes are identical and depend on the power budget P x only. Moreover, in the high transmit power regime, we observe that the performance of Baseline Scheme 2 also does not depend on ξ m , m = {1, 2}, for large numbers of transmit antennas, i.e., N T = {3, 4}, since, for large N T , EH node 1 is driven into saturation anyways. However, for the other system setups, the choice of the weights ξ m , m ∈ {1, 2}, enables a trade-off between the powers harvested at the EH nodes, which is characterized by a convex harvested power region. Furthermore, by increasing weight ξ m , more power is harvested at EH node m at the expense of a reduction of the power harvested by the other node. Thus, by choosing the user weights, the TX can control the distribution of the harvested power among the users. In particular, for ξ 1 = 1 (and ξ 2 = 0), the TX maximizes the average harvested power at EH node 1 and neglects EH node 2, which may yield a substantial decrease of the power at EH node 2. In the high transmit power regime, EH node 1 is driven into saturation for N T = {3, 4}. In this case, by decreasing ξ 1 (and increasing ξ 2 ), it is possible to significantly increase the power harvested by EH node 2 without a substantial reduction of the power harvested by EH node 1.

VI. CONCLUSION
In this paper, we considered multi-user MIMO WPT systems with multiple EH nodes employing non-linear rectennas. Based on a set of assumptions, which are satisfied for practical EH circuits, we specified a general EH model. Then, we proposed an optimal transmit strategy that maximizes the weighted sum of the average harvested powers at the EH nodes under a constraint on the power budget of the TX. For MISO WPT, we showed that transmission of scalar symbols with discrete random magnitudes, whose pdf has at most two mass points, via MRT beamforming is optimal. Next, for SIMO WPT, we proved that the optimal transmit symbol magnitude also has a discrete pdf with no more than two mass points. Then, for MIMO WPT, we showed that the optimal transmit strategy employs a scalar unit norm symbol and at most two beamforming vectors. In order to obtain these vectors, we proposed and solved a non-convex optimization problem. Since the computational complexity of the optimal solution was high, we developed an iterative low-complexity algorithm to obtain a suboptimal solution. Our simulation results revealed that the proposed optimal and suboptimal schemes for MIMO WPT systems yield practically identical performance. Furthermore, we observed that the proposed MIMO WPT design achieves substantial performance gains compared to two baseline schemes, one based on a linear EH model and the other one based on a single beamforming vector. Moreover, for multi-user MIMO WPT systems, we showed that the harvested power saturates for large numbers of TX antennas.
Finally, we observed a trade-off between the powers harvested at individual EH nodes, which we characterized by a convex harvested power region. APPENDIX A: PROOF OF LEMMA 1 In the following, we prove Lemma 1. First, we note that since ν * 2 is the maximizer of the slope function σ(·, ·; f ) for ν 1 = ν * 1 , then we have Then, since ν * 1 is the minimizer of σ(·, ·; f ) for ν 2 = ν * 2 , we have Next, we subtract f (ν * 1 )ν * 1 from both sides of (24). This allows us to rewrite both (22) and (24) as follows: respectively, which, thus, holds ∀ν ∈ R. Let us define linear function , where the inequality holds with equality if the pdf of ν is given by p * . This concludes the proof.
APPENDIX D: PROOF OF PROPOSITION 1 We solve optimization problem (2) for a single-user MISO WPT system, i.e., M = N E M = 1. First, let us consider a distribution of the transmit symbols x which has a point of increase atx 0 .
For this distribution, a larger value of the input power at the EH and, thus, an equal or larger value of Φ(·) can be attained by removing the mass pointx 0 and increasing the probability of symbol x 0 = x 0 2 g H g 2 exp(jθ s ) by the probability of symbolx 0 of the former distribution, see Assumptions 3 [5]. We note that this transformation preserves the validity of the distribution, i.e., p x (x)dx = 1, and, since the transmit powers for the two symbols are identical, i.e., x 0 2 2 = x 0 2 2 , the new distribution does not affect the power budget of the TX. Therefore, for the solution of (2), transmit vector x = wr s exp(jθ s ) is optimal, where w = g H g 2 is the MRT beamformer and r s = |s| and θ s are the magnitude and the arbitrary phase of random scalar symbol s, respectively. We denote the pdf of the transmit power values ν = r 2 s , ν ∈ [0, +∞), by p ν (ν). Then, the utility function in (3) can be rewritten as a function of pdf p ν (ν) as follows:
APPENDIX E: PROOF OF PROPOSITION 2 In order to prove Proposition 2, we note that as a sum of non-decreasing functions, function Φ(·) is also monotonically non-decreasing. Furthermore, the objective function can be equivalently rewritten as follows: where p ν (ν) is the pdf of the transmit power ν = r 2 x . Hence, for the considered SIMO WPT system, optimization problem (2) can be equivalently reformulated as follows: Since optimization problem (31) is in the form of auxiliary problem (4), the application of Corollary 2 yields Proposition 2. This concludes the proof.
APPENDIX G: PROOF OF PROPOSITION 3 First, we note that for any arbitrary transmit symbolx, there is a symbolx given bŷ which has the same transmit power and yields a higher or equal value of Ψ(x). Hence, for any arbitrary distribution of transmit symbols with a point of increasex, a larger value of Ψ(x) can be obtained by removing this point and increasing the probability ofx by the corresponding value.
Let us introduce now a function Φ(ν) that returns the largest possible value of Ψ(x) if a symbol with power ν was transmitted. This function is given by (14). We note that function Φ(·) is monotonically non-decreasing, see Assumption 3. Then, the solution of (2) can be obtained by determining first the solution p * ν (ν) of the following optimization problem: Since (34) is in the form of (4), there exists an optimal discrete pdf p * ν (ν) consisting of at most two mass points, ν * 1 and ν * 2 , see Corollary 2. Hence, the optimal symbol vector x can be decomposed as x = ws with unit-norm symbols s and discrete random beamforming vector w, whose pdf consists of at most two mass points evaluated as with probabilities p * w (w * n ) = p * ν (ν * n ), n ∈ {1, 2}, respectively. This concludes the proof.