Stochastic Analysis of the Efficiency of a Wireless Power Transfer System Subject to Antenna Variability and Position Uncertainties

The efficiency of a wireless power transfer (WPT) system in the radiative near-field is inevitably affected by the variability in the design parameters of the deployed antennas and by uncertainties in their mutual position. Therefore, we propose a stochastic analysis that combines the generalized polynomial chaos (gPC) theory with an efficient model for the interaction between devices in the radiative near-field. This framework enables us to investigate the impact of random effects on the power transfer efficiency (PTE) of a WPT system. More specifically, the WPT system under study consists of a transmitting horn antenna and a receiving textile antenna operating in the Industrial, Scientific and Medical (ISM) band at 2.45 GHz. First, we model the impact of the textile antenna’s variability on the WPT system. Next, we include the position uncertainties of the antennas in the analysis in order to quantify the overall variations in the PTE. The analysis is carried out by means of polynomial-chaos-based macromodels, whereas a Monte Carlo simulation validates the complete technique. It is shown that the proposed approach is very accurate, more flexible and more efficient than a straightforward Monte Carlo analysis, with demonstrated speedup factors up to 2500.


Introduction
With the advent of the Internet of Things (IoT), radio frequency identification (RFID) systems and, in general, radio frequency (RF) sensors and actuators have acquired significant importance. Distributed in our surroundings in a pervasive and inconspicuous way, these elementary components interact with each other to collect, process and exchange data [1,2]. The potential applications leveraged by the IoT paradigm are manifold and range from transportation over logistics to healthcare and from infrastructure monitoring to emergency services. As a result, assessing the correct operation of the systems involved is of paramount importance.
One of the main requirements of the IoT is that RFID components and RF sensors need to be small, low cost and not limited in lifetime by the duration of a battery [1]. Therefore, they are often designed to be passive and they rely on wireless power transfer (WPT) for their activation and operation [3,4]. Up to now, WPT in the reactive near-field and in the far-field has been widely investigated in literature [5][6][7][8]. However, the proposed solutions can be considered sub-optimal for several reasons. On the one hand, WPT operating in the reactive near-field can achieve high power transfer efficiency (PTE), but it requires source and target to be very close to each other [9]. Moreover, variations in distance between the devices strongly affect the resonance frequency and the load impedance of the system. On the other hand, the PTE achieved by far-field schemes is often too low. As a result, research efforts have recently shifted to WPT in the radiative near-field (Fresnel region), where the PTE, although rapidly decreasing with distance, can still be high enough to power sensor networks [10][11][12].
The most effective way to assess the performance of a WPT system is to study its transfer efficiency, since random fluctuations may alter the incoming power at the receiving element, whereas the output voltage level delivered by the power management system is typically first stabilized by a voltage regulator to provide a constant supply voltage needed to feed electronic circuits. To this end, several numerical techniques may be leveraged, such as source reconstruction methods [13], multipole expansion of the electromagnetic field [14], model reduction combined with full-wave solvers [15] and domain decomposition methods [16]. A very efficient approach has recently been proposed in [17,18], where the electromagnetic interaction among arbitrarily positioned radiating devices is modeled by means of their measured or simulated radiation patterns, allowing for device repositioning without requiring new simulations or measurements. Moreover, no constraint is enforced on the antenna configurations considered, as long as their radiation pattern is provided.
Although this approach efficiently quantifies the PTE of a WPT system for different positions and rotations of the antennas, it cannot assess how the PTE varies for a given system setup when devices undergo small random rotations or when their mutual position is affected by uncertainties. Moreover, the method is based on a single simulation of the radiation pattern of each element of the system, which typically corresponds to its nominal design. However, the radiation characteristics of the actually deployed antennas deviate from the nominal ones due to the inevitable uncertainties on the design parameters arising during the production process. As a result, the antenna variability has to be included in the WPT model to correctly estimate the PTE of the system.
In this manuscript, we propose a stochastic collocation method (SCM) based approach, leveraging on gPC expansions [19,20]. The formalism overcomes the limitations mentioned in the previous paragraph as it efficiently allows assessing the impact of position uncertainties and antenna variability on the PTE of a WPT system. More specifically, we start from given probability density functions (PDFs) according to which the antennas' design parameters vary. Next, we introduce a generalized polynomial chaos (gPC) expansion for each antenna configuration to model the corresponding variations in its radiation characteristics. Furthermore, on a higher level, a second gPC expansion assesses the impact of both antenna variability and position uncertainties on the PTE of the system. In this way, the analysis of the WPT system requires a lower number of simulations compared to a single gPC expansion accounting for all variations.
The proposed approach is demonstrated on a simple WPT system consisting of a transmitting horn antenna and a receiving Industrial, Scientific and Medical (ISM) textile antenna operating at 2.45 GHz, which is connected to a rectifier circuit to deliver direct current (DC) power to a load [18]. For this receiving antenna, we use experimentally determined probability density functions of the design parameters [21,22]. The results are found to be as accurate as those obtained by the traditional Monte Carlo method, here used to validate our technique. However, owing to the considerably lower number of simulations needed to construct a gPC expansion compared to the number required for Monte Carlo method to converge, the proposed approach proves to be much more efficient and flexible. The polynomial chaos expansion has been applied to model lumped circuits and distributed interconnects [23,24], multiport systems [25], the effect of geometric and material variations in scattering problems [26,27], in Direction of Arrival (DOA) estimation [28] and in antenna design [21,22]. However, to our best knowledge, the application of uncertainty quantification to a WPT system is completely new.
This manuscript is organized as follows. In Section 2, we briefly describe both the WPT model [17] and the SCM. Then, in Section 3, the results for a WPT consisting of a transmitting horn antenna and a receiving ISM textile antenna connected to a rectifier circuit are presented and then discussed in Section 4. Conclusions are summarized in Section 5.
Notations: We denote field vectors by underlined letters, e.g., v, and unit vectors with a "hat", e.g.,v. All sources and fields are assumed to be time harmonic with angular frequency ω and time dependencies e jωt are suppressed. Vector elements and arrays are represented by boldface characters, e.g., x. For a given array x ∈ C, x T denotes its transpose, whereas x denotes its Euclidean norm.

Wireless Power Transfer System Model
Consider a simple WPT system consisting of a transmitting antenna TX and a receiving antenna RX. We assume, for simplicity, that both TX and RX are one-port devices with radiation impedances Z TX and Z RX , respectively, and we represent them by means of two equivalent circuits [29] as in Figure 1. More specifically, the transmitter is driven by means of a Thévenin generator composed of a sinusoidal voltage source V g and an internal impedance Z g , whereas the receiving antenna is modeled as a Norton equivalent with a short circuit current I sc and load impedance Z L . The receiver also includes a matching circuit and a rectifier to transform the alternating current (AC) power received by the antenna into the direct current (DC) power delivered to the load R L . Notations: We denote field vectors by underlined letters, e.g., v, and unit vectors with a "hat", e.g.,v. All sources and fields are assumed to be time harmonic with angular frequency ω and time dependencies e jωt are suppressed. Vector elements and arrays are represented by boldface characters, e.g., x. For a given array x ∈ C, x T denotes its transpose, whereas x denotes its Euclidean norm.

Wireless Power Transfer System Model
Consider a simple WPT system consisting of a transmitting antenna TX and a receiving antenna RX. We assume, for simplicity, that both TX and RX are one-port devices with radiation impedances Z TX and Z RX , respectively, and we represent them by means of two equivalent circuits [29] as in Figure 1. More specifically, the transmitter is driven by means of a Thévenin generator composed of a sinusoidal voltage source V g and an internal impedance Z g , whereas the receiving antenna is modeled as a Norton equivalent with a short circuit current I sc and load impedance Z L . The receiver also includes a matching circuit and a rectifier to transform the alternating current (AC) power received by the antenna into the direct current (DC) power delivered to the load R L . The short circuit current I sc is computed as: where V 0 is a pertinent normalization factor [30], V RX is the volume of RX, e inc (r) is the electric field incident on the receiver, j(r) is the current density impressed on in V RX and r is the position vector. It was shown in [17] that, in the radiative near-field, I sc can be expressed as follows: where we integrate over the Ewald sphere Ω. Furthermore, F TX (k) and F RX (k) are the radiation patterns of the transmitter and the receiver, respectively, Z = µ is the wave impedance of the background medium,k = sin θ cos φx + sin θ sin φŷ + cos θẑ is the wave vector in spherical coordinates, T(r TX,RX ,k) is a translation operator that allows efficient translations between the two antennas and r TX,RX is the relative position between the phase centers of the two devices. The operator T(r TX,RX ,k) is calculated as follows: where j is the imaginary unit, h l (·) is the l-th order spherical Hankel function of the second kind, and P l (·) is the Legendre polynomial of degree l. The number L determines the accuracy of the approximation in Equation (3) and traditional guidelines are followed to select it [31]. The relation in Equation (2) is valid only as long as the two antennas are not positioned in each other's reactive near-field, thus with |r TX,RX | at least equal to a sixth of the wavelength. The short circuit current I sc is computed as: where V 0 is a pertinent normalization factor [30], V RX is the volume of RX, e inc (r) is the electric field incident on the receiver, j(r) is the current density impressed on in V RX and r is the position vector. It was shown in [17] that, in the radiative near-field, I sc can be expressed as follows: where we integrate over the Ewald sphere Ω. Furthermore, F TX (k) and F RX (k) are the radiation patterns of the transmitter and the receiver, respectively, Z = µ is the wave impedance of the background medium,k = sin θ cos φx + sin θ sin φŷ + cos θẑ is the wave vector in spherical coordinates, T(r TX,RX ,k) is a translation operator that allows efficient translations between the two antennas and r TX,RX is the relative position between the phase centers of the two devices. The operator T(r TX,RX ,k) is calculated as follows: where j is the imaginary unit, h l (·) is the l-th order spherical Hankel function of the second kind, and P l (·) is the Legendre polynomial of degree l. The number L determines the accuracy of the approximation in Equation (3) and traditional guidelines are followed to select it [31]. The relation in Equation (2) is valid only as long as the two antennas are not positioned in each other's reactive near-field, thus with |r TX,RX | at least equal to a sixth of the wavelength.
The rotation of the devices around their phase center is included in the described formalism by applying the appropriate rotation R to the radiation patterns F TX (k) and F RX (k) in the spherical harmonics domain, as shown in Figure 2. Thereto, first, the radiation pattern F(k) = F θ (k)θ + F φ (k)φ is expanded into spherical harmonics A pq and B pq by means of the transformation F , given by [32]: where the unit vectork has been replaced by (θ, φ), Y pq (θ, φ) and their complex conjugates Y * pq (θ, φ) are the orthonormalized scalar spherical harmonics. Then, the coefficients A pq and B pq are rotated in the spherical harmonics domain (R SH ) by means of Wigner D-matrices. The rotated coefficients A R pq and B R pq are calculated as [33] A with d r pq (β) the Wigner small d-matrix, given by [34] d r The range of s is determined such that all factorials are nonnegative. α, β and γ in Equations (5) and (6) are the standard Euler angles that define the rotation using the z-y-z convention in a right-handed frame. The angles (α, β, γ) are related to the desired inclination and azimuthal angles θ and φ by choosing α = φ, β = θ and γ = 0. Finally, the rotated radiation pattern F R (θ, φ) is found from the rotated coefficients A R pq and B R pq by means of the transformation F −1 given by [32]: where P is a parameter that sets the accuracy, which for practical reasons may be chosen equal to L in Equation (3) [35]. The rotation of the devices around their phase center is included in the described formalism by applying the appropriate rotation R to the radiation patterns F TX (k) and F RX (k) in the spherical harmonics domain, as shown in Figure 2. Thereto, first, the radiation pattern F(k) = F θ (k)θ + F φ (k)φ is expanded into spherical harmonics A pq and B pq by means of the transformation F , given by [32]: where the unit vectork has been replaced by (θ, φ), Y pq (θ, φ) and their complex conjugates Y * pq (θ, φ) are the orthonormalized scalar spherical harmonics. Then, the coefficients A pq and B pq are rotated in the spherical harmonics domain (R SH ) by means of Wigner D-matrices. The rotated coefficients A R pq and B R pq are calculated as [33] A with d r pq (β) the Wigner small d-matrix, given by [34] d r The range of s is determined such that all factorials are nonnegative. α, β and γ in Equations (5) and (6) are the standard Euler angles that define the rotation using the z-y-z convention in a right-handed frame. The angles (α, β, γ) are related to the desired inclination and azimuthal angles θ and φ by choosing α = φ, β = θ and γ = 0. Finally, the rotated radiation pattern F R (θ, φ) is found from the rotated coefficients A R pq and B R pq by means of the transformation F −1 given by [32]: where P is a parameter that sets the accuracy, which for practical reasons may be chosen equal to L in Equation (3) [35].
We notice from Equation (2) that only the (measured or simulated) radiation patterns of the two devices are needed to calculate the influence of the transmitter on the receiver. Moreover, the combination of the rotation mechanism of Figure 2 and the use of the translation operator T(r TX,RX ,k) in Equation (2) allow computing the short circuit current I sc for any set of rotations and positions of the devices. As a result, this formalism is more efficient and flexible than traditional simulations tools or measurements, which require a new computation or measurement for every rotation and repositioning. We notice from Equation (2) that only the (measured or simulated) radiation patterns of the two devices are needed to calculate the influence of the transmitter on the receiver. Moreover, the combination of the rotation mechanism of Figure 2 and the use of the translation operator T(r TX,RX ,k) in Equation (2) allow computing the short circuit current I sc for any set of rotations and positions of the devices. As a result, this formalism is more efficient and flexible than traditional simulations tools or measurements, which require a new computation or measurement for every rotation and repositioning.
The PTE of the WPT system considered consists of three contributions. First, we need to calculate the wireless link efficiency η link , which is defined as the ratio between the power P RX delivered to the receiving antenna's load Z L and the power P TX emitted by the transmitter: Then, the efficiencies η match and η rect of the matching and the rectifying circuits are included. Since the input impedance Z L of a nonlinear circuit depends on the incoming power, η match has the same dependency. Therefore, we simulate it with the commercial tool Advanced Design System (ADS) by Keysight Technologies (Santa Rosa, CA, USA). As to the efficiency of the rectifier, it is calculated as η rect = P inc /P DC , where P inc = η match · P RX and P DC are the AC power injected into the rectifier and the DC power delivered to the load R L , respectively. The DC power P DC is given by: where the DC output voltage V out is also computed by ADS. Finally, the overall PTE of the system is given by:

Stochastic Collocation Method
Consider a generic system output G (such as the PTE of the WPT system) and N input random variables (RV) x 1 , x 2 , ..., x N that affect it. We assume these variables to be actually independent and we collect them in the vector x = [x 1 , x 2 , ..., x N ]. Then, following the Wiener-Askey scheme [19], we relate G to x by means of a polynomial chaos expansion where ϕ k (x) are suitably chosen multivariate polynomial basis functions and the expansion coefficients y k are the unknowns to be determined. The polynomials ϕ k (x) are constructed to be orthonormal with respect to the PDF P X , which describes the likelihood of the input x. Thus: with δ jk = 0 if i = j, δ jk = 1 if i = j, and Γ being the support of P X . Consequently, P X acts as a weighting function. Since the considered input RVs are independent, the PDF P X is defined as the product of the PDFs corresponding to the single input variables. Therefore, the polynomials ϕ k (x) are constructed as products of N univariate polynomials, each one associated to a single input RV. Furthermore, the multivariate polynomials have a total degree of maximally Q, meaning that the sum of the orders of the univariate polynomials is at most Q.
Finally, by means of the Stochastic Testing (ST) algorithm [36] described in the Appendix, we select M = K collocation points x m and we compute both a matrix A, whose elements are defined as a mk = ϕ k (x m ), and its inverse B. The coefficients y k in Equation (11) are then calculated as in [37]: where b mk is the mk-th element of matrix B and f (x m ) is the function to be approximated, evaluated in x m .

Wireless Power Transfer Uncertainty Quantification
The efficiency of the WPT system under study is affected by two main issues. On the one hand, both the transmitter and the receiver may undergo uncertainties in their design parameters, since, for example, the production process yields devices that do not perfectly correspond to their nominal design. As a result, random variations may be expected on both radiation impedances Z TX and Z RX , and on the antennas' radiation patterns. On the other hand, the system may be conceived to operate in a given configuration and undesired small random rotations or variations in the mutual position of the antennas affect the influence of the transmitter on the receiver.
Conventionally, these problems are expected to arise simultaneously. Therefore, the variations in the WPT efficiency corresponding to both antenna variability and position uncertainies are investigated by means of gPC expansions in Equation (11), which are defined as functions of the design parameters of the antennas and of the WPT system. However, the use of a single gPC expansion to carry out the complete analysis is sub-optimal. This is understood as follows. The construction of a gPC expansion requires selecting and processing M collocation points x m to compute the coefficients y k in Equation (11). Even though translations and rotations are efficiently accounted for by the WPT model described in Section 2.1, both the antennas' radiation impedances and their radiation patterns are usually computed by means of full-wave solvers, such as ADS Momentum (Keysight Technologies, Santa Rosa, CA, USA), whose simulations are typically time-consuming. As a result, accounting for the variability of the antennas requires a high number of full-wave simulations. Moreover, if the efficiency of the WPT system is evaluated for different system configurations or different position uncertainties, a new gPC expansion has to be constructed for each configuration and for each position. As a result, for each such expansion a new set of collocation points has to be selected and processed, again requiring full-wave simulations, which significantly decrease the efficiency of the method with respect to more naive approaches, such as Monte Carlo.
In order to overcome these limitations and drastically improve efficiency, the construction of the gPC expansions that model the efficiency of the whole system is preceded by an intermediate step.
More specifically, for each antenna configuration undergoing variability, we introduce gPC expansions in Equation (11) of the real and the imaginary parts of both its radiation impedance Z and the spherical harmonic coefficients A pq and B pq in Equation (4): where x VAR is the vector of the antenna's design parameters that are affected by uncertainties. These gPC expansions can now be interpreted as macromodels of the considered antennas. Once constructed, they allow accurately computing both an antenna's radiation impedance and radiation pattern for a value of x VAR without having to resort to full-wave simulations. After this intermediate step, we model both the link efficiency η LINK and the overall PTE of the system by means of the following gPC expansions: where x WPT is the vector of all the parameters in the WPT system subject to variations, which comprises all the vectors x VAR and the variables corresponding to position uncertainties. The procedure is now as follows. First, M η LINK and M PTE collocation points x WPT m are selected to construct the gPC expansions in Equations (20) and (21), respectively. Next, for each collocation point, the gPC expansions in Equations (14)- (19) are used to rapidly compute both the radiation impedances and the radiation patterns of the antennas undergoing variability. Finally, with the calculated antennas' radiation characteristics, the WPT model described in Section 2.1 processes the values of x WPT m corresponding to the position uncertainties in the system and computes the values of η LINK and PTE in those collocation points required to calculate the coefficients y k 7 ,y k 8 in Equations (20) and (21).
The benefits of this approach are twofold. First, since for each antenna the number of design parameters undergoing variations is expected to be significantly lower than the total number of parameters affecting the system, the amount of full-wave simulations necessary to construct all antenna macromodels is expected to be substantially lower than what is required to compute a single gPC expansion that models the entire WPT system. Second, once these macromodels are available, the analysis of the WPT system can be repeated for any given system configuration at a negligible computational cost, since no full-wave simulations are required to calculate the antennas' radiation characteristics.

Validation Example Setup
The proposed approach is demonstrated on the WPT system shown in Figure 3. The transmitting device is a standard gain horn (SGH) antenna radiating at 2.45 GHz with a power of 10 dBm. The receiving device is a 2.4-2.4835 GHz ISM band textile microstrip probe-fed dual-polarized patch antenna [38], using a flexible closed-cell expanded rubber foam as substrate, and shown in Figure 4. The antennas are placed at a separation distance d and their phase centers are aligned. where x VAR is the vector of the antenna's design parameters that are affected by uncertainties. These gPC expansions can now be interpreted as macromodels of the considered antennas. Once constructed, they allow accurately computing both an antenna's radiation impedance and radiation pattern for a value of x VAR without having to resort to full-wave simulations. After this intermediate step, we model both the link efficiency η LINK and the overall PTE of the system by means of the following gPC expansions: where x WPT is the vector of all the parameters in the WPT system subject to variations, which comprises all the vectors x VAR and the variables corresponding to position uncertainties. The procedure is now as follows. First, M η LINK and M PTE collocation points x WPT m are selected to construct the gPC expansions in Equations (20) and (21), respectively. Next, for each collocation point, the gPC expansions in Equations (14)- (19) are used to rapidly compute both the radiation impedances and the radiation patterns of the antennas undergoing variability. Finally, with the calculated antennas' radiation characteristics, the WPT model described in Section 2.1 processes the values of x WPT m corresponding to the position uncertainties in the system and computes the values of η LINK and PTE in those collocation points required to calculate the coefficients y k 7 ,y k 8 in Equations (20) and (21).
The benefits of this approach are twofold. First, since for each antenna the number of design parameters undergoing variations is expected to be significantly lower than the total number of parameters affecting the system, the amount of full-wave simulations necessary to construct all antenna macromodels is expected to be substantially lower than what is required to compute a single gPC expansion that models the entire WPT system. Second, once these macromodels are available, the analysis of the WPT system can be repeated for any given system configuration at a negligible computational cost, since no full-wave simulations are required to calculate the antennas' radiation characteristics.

Validation Example Setup
The proposed approach is demonstrated on the WPT system shown in Figure 3. The transmitting device is a standard gain horn (SGH) antenna radiating at 2.45 GHz with a power of 10 dBm. The receiving device is a 2.4-2.4835 GHz ISM band textile microstrip probe-fed dual-polarized patch antenna [38], using a flexible closed-cell expanded rubber foam as substrate, and shown in Figure 4. The antennas are placed at a separation distance d and their phase centers are aligned.  The radiation impedance Z TX of the SGH antenna is 50 Ω and its radiation pattern is computed by analytical expressions (see the Appendix in [17], where a = 368.6 mm, b = 273.1 mm, ρ 1 = 343.9 mm, ρ 2 = 363.4 mm). Instead, the 2.45 GHz ISM band textile antenna is designed and simulated by means of ADS Momentum (Keysight Technologies, Santa Rosa, CA, USA). Its nominal design parameters are reported in Table 1. Its nominal radiation impedance Z RX , which is equal for both feed 1 and feed 2, is 49.91 − 1.93i Ω at the frequency of 2.45 GHz. An isolation −20log|S 21 | between the two feed ports of more than 15 dB is achieved within the entire ISM band.  Table 1. The wireless link efficiency η link between the SGH and the textile antenna is computed by means of the formalism described in Section 2.1. The parameters L and P in Equations (3) and (7), respectively, are both set to five, leading to 36 coefficients A pq and B pq . Next, a rectifier is attached to the patch antenna to form the rectenna shown in Figure 5, which is designed and simulated in ADS. More specifically, the rectenna consists of a matching network, a voltage doubler and a rectifier. The matching circuit is given by an inductor L m = 5 nH, whereas the voltage doubler and rectifier circuit itself consist of two HSMS-2850 Schottky diodes, for which a pertinent SPICE model is used, together with their package parasitics C p = 0.08 pF and L p = 2 nH. The capacitors C 1 and C 2 are equal to 100 pF and the load resistance is equal to R L = 100 Ω. The matching circuit is designed to have optimal matching when P RX = −10 dBm. However, because of the nonlinear diodes in the voltage doubler and rectifier circuit, the load impedance Z L , and, therefore, the matching efficiency η match , depend on the incoming power P RX . This impedance is computed for different incoming powers by using a harmonic balance simulation in ADS and the matching efficiency is then calculated as: Finally, the efficiency η rect of the voltage doubler and rectifier is calculated as in Equation (9), and the total efficiency η tot of the rectenna is given by η tot = η match · η rect .
Sensors 2016, 16, 1100 9 of 15 incoming power P RX . This impedance is computed for different incoming powers by using a harmonic balance simulation in ADS and the matching efficiency is then calculated as: Finally, the efficiency η rect of the voltage doubler and rectifier is calculated as in Equation (9), and the total efficiency η tot of the rectenna is given by η tot = η match · η rect .
Patch Antenna Matching Voltage Doubler and Rectifier Load

Antenna Variability
In the example under study, we only consider antenna variability of the 2.45 GHz ISM band textile antenna, since an SGH antenna is not expected to deviate from its nominal characteristic. Textile antennas, however, usually exhibit variations in their design parameters, which significantly alter their figures of merit. On the one hand, the production process inevitably introduces uncertainties in the geometrical dimensions of the antenna. On the other hand, the flexible closed-cell expanded rubber protective foam used as antenna substrate is a very non-uniform material. As a result, the value of its electrical permittivity r may be quite different from the nominal one.
In [21], it is shown that variations in the patch length L, the patch width W and the electrical permittivity r of the antenna have a significant impact on its radiation impedance Z RX , whereas the influence of other parameters is negligible. Moreover, a quick sensitivity analysis confirms that only variations in L, W and r affect the radiation pattern of the antenna. Therefore, we relate both the real and the imaginary part Z re RX and Z im RX , respectively, of the radiation impedance Z RX of the receiving antenna, as well as the real and the imaginary part of the coefficients A pq and B pq in Equation (4), to the parameters L, W and r by means of the gPC expansions in Equations (14)- (19) where thus, x VAR = [L, W, r ]. As we know from [21,22] that these parameters are independent and vary according to Gaussian distributions, according to the Wiener-Askey scheme [19], the multivariate (14)- (19) consist of products of Hermite polynomials. Finally, the PDF P X in Equation (12) is given by the product of three univariate Gaussian distributions. Its mean vector µ and covariance matrix Σ are given by:

Antenna Variability
In the example under study, we only consider antenna variability of the 2.45 GHz ISM band textile antenna, since an SGH antenna is not expected to deviate from its nominal characteristic. Textile antennas, however, usually exhibit variations in their design parameters, which significantly alter their figures of merit. On the one hand, the production process inevitably introduces uncertainties in the geometrical dimensions of the antenna. On the other hand, the flexible closed-cell expanded rubber protective foam used as antenna substrate is a very non-uniform material. As a result, the value of its electrical permittivity r may be quite different from the nominal one.
In [21], it is shown that variations in the patch length L, the patch width W and the electrical permittivity r of the antenna have a significant impact on its radiation impedance Z RX , whereas the influence of other parameters is negligible. Moreover, a quick sensitivity analysis confirms that only variations in L, W and r affect the radiation pattern of the antenna. Therefore, we relate both the real and the imaginary part Z re RX and Z im RX , respectively, of the radiation impedance Z RX of the receiving antenna, as well as the real and the imaginary part of the coefficients A pq and B pq in Equation (4), to the parameters L, W and r by means of the gPC expansions in Equations (14)- (19) where thus, x VAR = [L, W, r ]. As we know from [21,22] that these parameters are independent and vary according to Gaussian distributions, according to the Wiener-Askey scheme [19], the multivariate (14)- (19) consist of products of Hermite polynomials. Finally, the PDF P X in Equation (12) is given by the product of three univariate Gaussian distributions. Its mean vector µ and covariance matrix Σ are given by: where σ L , σ W and σ r are the standard deviations of L, W and r , respectively, whereas ρ LW = ρ L r = ρ W r = 0 indicates that there is no correlation between the design parameters. The gPC expansions of Z re RX and Z im RX are found to converge for orders of expansion with a total degree Q Z re RX and Q Z im RX equal to 2 and 3, respectively. In contrast, the real and the imaginary part of the coefficients A pq and B pq exhibit a more intricate behavior. Therefore, a total degree Q A pq = Q B pq = 6 is necessary to accurately catch the variations determined by L, W and r on these parameters. Since both the radiation impedance Z RX and the radiation pattern of the antenna are computed during the same full-wave simulation, the total degree of all gPC expansions in Equations (14)- (19) is set to Q VAR = 6. In this way, only one set of collocation points is selected and simulated in ADS to compute the coefficients y k 1 , y k 2 , y k 3 , y k 4 , y k 5 , y k 6 in Equations (14)- (19) and to construct all the gPC expansions, which is clearly beneficial in terms of computation time. In this case, a number M VAR = 84 of collocation points x VAR m has been selected by means of the ST algorithm and simulated in ADS.

Position Uncertainties
Nominally, the WPT system shown in Figure 3 is constructed to operate in the radiative near-field with a distance d between the two antennas equal to 0.6 m ≈ 5λ and with their phase centers aligned at coordinates x = y = 0. Furthermore, the radiating surface of the SGH antenna and the patch antenna are aligned with the xy plane, which means that the rotation angles θ and φ are both equal to 0 • . In practice, the antennas are perturbed by small random rotations and variations in their mutual position. For the sake of conciseness, the SGH antenna is stationary and all the variations are cumulated in the position and the rotation of the 2.45 GHz ISM band antenna. Finally, since no experimental data are available to estimate the PDFs corresponding to the considered parameters, we suppose that d, x, y, θ and φ are independent and vary according to Gaussian distributions, whose mean values and standard deviations are reported in Table 2.
More specifically, we assume that the variations in the position of the ISM patch antenna in the xy plane are limited to intervals x = ±2 cm and y = ±2 cm, which correspond to a displacement of about half the width W and the length L of the patch. As for d, θ and φ, the variations are assumed to be large enough to account for a potential displacement in a real scenario.

Wireless Power Transfer Efficiency
The two gPC expansions in Equations (20) and (21) are introduced to relate both the link efficiency η LINK and the overall PTE to the parameters L, W, r , d, x, y, θ and φ. All parameters are independent and vary according to Gaussian distributions. As a result, the multivariate polynomials ϕ k 7 (x WPT ), ϕ k 8 (x WPT ) of both the gPC expansions are Hermite polynomials, as in Section 3.2. Note, however, that upon the availability of other experimental data, other distributions for these parameters can equally be dealt with by means of the advocated gPC-based approach. We find that both expansions converge for an order of expansion corresponding to a total degree Q η LINK and Q PTE equal to 4. A number M η LINK = M PTE = 495 of collocation points x LINK m is selected by means of the ST algorithm and processed with the WPT model and the antenna macromodels in order to compute the coefficients y k 7 , y k 8 in Equations (20) and (21).
Next, we perform a Monte Carlo analysis of both η LINK and PTE by processing a sample set of 10,000 realizations of L, W, r , d, x, y, θ and φ, drawn according to their pertinent PDFs. Then, we compute the cumulative distribution functions (CDFs) of η LINK and PTE based on both the Monte Carlo simulation and the SCM analysis. The curves are shown in Figures 6 and 7. We notice that the CDFs are perfectly overlapping. Finally, in order to validate our analysis, we apply the Kolmogorov-Smirnov test to verify whether the CDFs of η LINK and PTE and those found by means of the SCM analysis correspond to the same distribution. In particular, the maximum distance D between them is compared to a threshold distance D α . For D < D α , the Kolmogorov-Smirnoff test accepts the null hypothesis that both the sample sets correspond to the same distribution, with a significance level α. If we set the significance level α to 0.05, D α equals 0.019233. The computed values of D η LINK and D PTE are equal to 0.0067 and 0.0085, respectively. Therefore, the null hypothesis of equality between the CDFs is validated with a significance level of 5%.

Discussion
The proposed approach allows quantifying the variations of both the link efficiency η LINK and the PTE of a WPT system in a more efficient and flexible way than both an SCM analysis based on a single gPC expansion and a Monte Carlo analysis. As shown in Table 3, both the number of full-wave simulations and the overall CPU time necessary to perform the analysis are greatly reduced for both η LINK and PTE: Table 3. Simulation data for the analysis of the η LINK and the power transfer efficiency (PTE) of the wireless power transfer (WPT) with different methods.

Method
Number In comparison, an analysis based on a single gPC, which requires direct simulation of 495 antenna realizations in ADS, takes more than 2 h. Moreover, this operation has to be repeated each time any difference is introduced in the distributions according to which the parameters of the WPT system vary. As for the Monte Carlo procedure, the simulation of 10,000 realizations by means of ADS and the WPT model of Section 2.1 requires more than 41 h. As a result, the proposed approach greatly outperforms both the Monte Carlo technique and an SCM analysis based on a single gPC expansion.

Conclusions
In this manuscript, an SCM analysis of the efficiency of a WPT system in the radiative near-field subject to antenna variability and position uncertainties has been presented. More specifically, a first SCM analysis is carried out to account for the impact of uncertainties in the design parameters of the antennas on their radiation characteristics. The resulting gPC expansions are used as macromodels that allow computing the antennas' radiation characteristics in a more efficient way than full-wave solvers. Then, a second SCM analysis quantifies the impact of both the variability of the deployed antennas and the uncertainties in their mutual position on the efficiency of the WPT system. This is done by leveraging the previously computed gPC-based macromodels and a very efficient model for WPT systems in the radiative near-field. Finally, the proposed approach is validated by means of a WPT system consisting of an SGH antenna and a 2.45 GHz ISM band textile antenna. Compared to an SCM analysis based on a single gPC expansion, as well as to a standard MC analysis, the method shows excellent agreement and superior efficiency.