Efficient Uncertainty Assessment in EM Problems via Dimensionality Reduction of Polynomial-Chaos Expansions †

The uncertainties in various Electromagnetic (EM) problems may present a significant effect on the properties of the involved field components, and thus, they must be taken into consideration. However, there are cases when a number of stochastic inputs may feature a low influence on the variability of the outputs of interest. Having this in mind, a dimensionality reduction of the Polynomial Chaos (PC) technique is performed, by firstly applying a sensitivity analysis method to the stochastic inputs of multi-dimensional random problems. Therefore, the computational cost of the PC method is reduced, making it more efficient, as only a trivial accuracy loss is observed. We demonstrate numerical results about EM wave propagation in two test cases and a patch antenna problem. Comparisons with the Monte Carlo and the standard PC techniques prove that satisfying outcomes can be extracted with the proposed dimensionality-reduction technique.


Introduction
Uncertainty quantification in the context of an Electromagnetic (EM) problem is of vital significance, as the calculation of the involved field quantities can be a challenging task. For instance, biological tissues are complicated media with high variability in their electric characteristics [1], and as a result, the utilization of deterministic approaches may not constitute a safe choice. Other problems involve the geometrical uncertainties introduced due to fabrication tolerances during the construction of printed circuit board antennas, which may have a significant impact on their performance [2]. Neglecting those random fluctuations can lead to unrealistic outcomes; therefore, deterministic schemes are not sufficient in such cases. For this reason, various techniques have been proposed that deal with uncertainty problems more efficiently.
The most common method for assessing EM uncertainties is the Monte Carlo (MC) approach [3]. In this context, a given problem is solved repeatedly, using various random samples of the input parameters. The convergence of the MC method may be achieved after a large number of simulations, which eventually make it impractical in many cases. A more efficient technique is based on Polynomial Chaos (PC) expansions [4]. This algorithm manages to extract reliable results in problems with low or moderate numbers of random variables. The PC scheme has been already utilized in many EM cases, where the stochastic inputs present uncertainties in electric [5] or geometric features [6]. However, as the dimensionality of a given problem grows, the computational cost of the PC approach increases, rendering it eventually less efficient.
In this paper, we perform a dimensionality reduction on the PC scheme by firstly applying a sensitivity analysis based on the Morris method [7]. In this way, the random variables with the smaller contributions to the variance of the output of interest are detected, and can be neglected without significant precision loss. Then, the PC technique is utilized with only the most influential stochastic inputs, reducing the computational times significantly. This approach is case dependent, as different sets of random variables may be the most important ones in different problems. The present paper generalizes the preliminary work of [8], and the proposed methodology is additionally applied to the more complex case of a patch antenna problem. Comparisons between the MC and conventional PC schemes indicate that satisfying outcomes can be extracted at a reduced computational cost.

Brief Literature Review of Related Works
In the pertinent literature, various suggestions have been proposed to tackle the limitations of the PC technique. A number of popular approaches make use of sparse grids based on the Smolyak method [9], which can significantly reduce the number of required simulations. However, Smolyak grids still suffer from the "curse of dimensionality" in problems with a high number of random inputs. For this reason, other techniques have been suggested that mitigate this shortcoming even further. The work in [10] utilizes an adaptive algorithm for the construction of nested sparse grids. This method starts by estimating the mean value of the examined function and proceeds to the calculation of the PC coefficients by taking advantage of the mean estimation in the previous step. As a result, fewer quadrature points are required, while preserving high accuracy. Alternative approaches are hierarchical sparse grids [11], which are based on piecewise multi-linear hierarchical basis functions. Hierarchical surpluses are utilized for error control and adaptively refine the collocation points in discontinuity regions in the stochastic space.
Other algorithms manage to reduce the number of terms in the PC method by seeking sparse solutions. For example, the authors of [12] proposed the utilization of hyperbolic index sets to truncate the PC expansions. Then, a sparse solution was constructed by performing an adaptive algorithm based on least angle regression. In [13], a weighted 1 -minimization approach was proposed for the computation of sparse PC representations. The weights in the 1 norm were computed via an approximation of the PC coefficients. As a result, coefficients with very small values were further penalized, improving the overall efficiency. In [14], a compressed sensing algorithm was presented, which exploits the concept of D-optimality. A design of experiments was constructed therein, utilizing the QR factorization with column pivoting. Then, the orthogonal matching pursuit method, which is a popular technique for finding sparse solutions, was properly modified to take into account these designs.
Additionally, tensor recovery algorithms have been successfully utilized for uncertainty quantification, reducing the number of required simulations in the PC method [15]. The tensor recovery is performed by applying an alternating minimization approach. The presented results in [15] indicated that the proposed method can be more efficient than sparse grids and the MC technique for the examined test cases. The work in [16] performed a tensor recovery algorithm in hierarchical uncertainty quantification problems. Low-level simulations were accelerated by utilizing the anchored analysis of variance method, while high-dimensional surrogate models were handled via tensor-train decomposition at the high level. Both approaches achieved a near-linear complexity with respect to the number of random inputs.

Polynomial Chaos Expansions
According to [4], the PC expansion can represent a random function y via a series of orthogonal polynomials. Specifically, if y depends on N stochastic variables ξ 1 , ξ 2 , . . . , ξ N , it can be expressed as: where ξ = [ξ 1 , ξ 2 , . . . , ξ N ] T . The c i parameters represent the polynomial coefficients and Ψ i are orthogonal basis functions. The expansion in (1) is approximated by truncating the infinite summation to P + 1 terms, with a value of: where k is the polynomial order. For cases with multiple independent random variables, the basis functions are constructed from univariate polynomials, as follows: where α i denotes the corresponding polynomial order [17] and ψ α i are the corresponding 1D polynomials. The choice of the basis depends on the distribution of each stochastic input. For instance, uniform variables require the Legendre basis, and Hermite polynomials are suitable for normal distributions.
After the approximation of the PC expansion has been determined, the corresponding mean value and variance can be computed as: where: Ω N is the N-dimensional random space, and pd f(ξ) denotes the joint probability-density function.
In order to estimate the expansion coefficients, two main techniques exist: intrusive and non-intrusive ones. The first approaches modify the deterministic solver by performing the computation of the PC expansion terms within the solver. On the contrary, the non-intrusive methods perform a number of deterministic realizations at specific collocation points in the random space. Then, one way to obtain the PC coefficients is through linear regression. In this case, given a total number of S collocation points in Ω N , the PC expansion must remain valid at each one of them. As a result, this leads to the following system of equations: where y 1 , y 2 , . . . , y S denote the outputs from the deterministic simulations. It is essential that the number of equations S be equal to or greater than the P + 1 coefficients. The overdetermined system in (7) can be solved through the least squares method as: An alternate way to estimate the coefficients is through the spectral projection method, which takes advantage of the orthogonality of the polynomial basis [18]. As a result: The integrals in (9) can be approximated with the help of quadrature. Specifically, after selecting an appropriate set of collocation points, (9) is estimated as: where w i are the weights and i max is the number of collocation points. In this work, the selection of collocation points is performed via the Clenshaw-Curtis nodal sets [19], which are based on Chebyshev polynomials. The Clenshaw-Curtis nodes are calculated as: Finally, the computation of (10) is implemented by utilizing the Smolyak algorithm [9], which manages to create a sparse representation of a full-tensor grid.

The Morris Method
In order to perform the dimensionality reduction in the PC approach, we apply a sensitivity analysis to the stochastic inputs. In this way, the trivial random variables can be identified, and thus treated as deterministic. The chosen sensitivity analysis algorithm is the Morris method [7], due to its low computational cost. In particular, this technique starts by defining a set of r (typically between 10 and 15) possible parameter values in the random space. Then, the quantity of interest f is calculated for a vector d j = [d 1 , d 2 , . . . , d n ] of this set, with 0 < j ≤ r. Next, each element of d j is perturbed by a factor ∆ j i at a time, while the others are kept unchanged. The values of ∆ j i are usually a predetermined multiple of 1/(r − 1). This step continues until all the elements of d j change. After that, the jth elementary effect [7] is calculated for each variable, as: This procedure is repeated for all the r points. Finally, the r elementary effects are averaged for every random variable: where m * i denotes the mean elementary effects [20] of the variable i. The stochastic inputs with a high mean elementary effect are considered influential, while the ones with a low mean elementary effect are treated as trivial. In this work, the Finite-Difference Time-Domain (FDTD) method [21] is used as a deterministic solver. This technique performs a discretization of the Maxwell equations both in time and space and computes the involved field components in a leapfrog manner (a brief description of this scheme is presented in Section 3.3). However, the quantities m * i have to be calculated for every cell of the discretized spatial grid. As a result, the following heuristic is employed.

•
Let g be the cells in the grid that satisfy m * j ≥ m * k .

•
Calculate the mean m * j , for the cells in g. Let this be M m * j [g] .

•
Compute the product len (g) M m * j [g] , where len (g) is the number of values in g.
Therefore, an N × N matrix is created, where each element depicts the significance of variable ξ j , compared to ξ k . Then, the mean value in each row (excluding the zeros in the main diagonal) is computed; thus, a vector J composed of N elements is extracted.

The Finite-Difference Time-Domain Technique
The FDTD technique is one of the most popular algorithms for solving full-wave propagation problems. Specifically, this scheme employs second-order central differences to Maxwell's differential equations. For example, in the one-dimensional case, the scalar equations of the Faraday and the Ampere laws are given by: where H y and E z denote the magnetic and electric strengths in the y and z directions, respectively. Furthermore, is the dielectric permittivity, and µ expresses the magnetic permeability. The time and spatial partial derivatives can be replaced with finite difference approximations. Applying this to (14) yields: where n indicates the time step index. Solving for H y | gives: where the values of ∆x and ∆t denote the spatial and time discretization steps. Since the space is discretized, the field components are updated for all the i cells in the FDTD grid [21]. The electric field values can be similarly computed by applying the same procedure to (15).

Numerical Results
The proposed approach was assessed via three different EM problems. The first numerical test involved a coaxial cable with eight uniform random dielectric materials, whose properties are shown in Table 1. Figure 1 depicts the geometric features for the 1D problem. We considered an incident Gaussian pulse (maximum frequency of 2 GHz), which emanated at a distance x = 0.015 m for 25 ns. The FDTD grid consisted of 1200 cells, with a discretization density equal to 40 cells per wavelength in the vacuum at 2 GHz. This problem was examined in two cases where in each case, the dielectric permittivities ranged between ±5% and ±8% of the corresponding mean values. The magnetic permeability was considered deterministic and equal to µ 0 . The level of parameter L in Smolyak's algorithm [9] and the order of the PC expansion were both set to three. Furthermore, the perturbation step ∆ j i was constant for every random variable and had a value of 0.5. Finally, unwanted reflections were minimized by applying the first-order Mur's absorbing boundary condition [22] at the two ends of the computational domain.  The reduced dimension PC approach was compared with 1000 MC realizations and the original PC method. In this problem, the Morris method required 108 runs for the estimation of the mean elementary effects, which was only a small amount of extra computational burden. The proposed algorithm solved the transmission-line problem using the six most influential variables. The significance of each random input was determined via a threshold, which was equal to 3 × 10 −7 and was applied in J. Therefore, the overall efficiency of the PC algorithm was increased. Figure 2a,b illustrate the mean and the standard deviation of the E field for the first case, respectively. Evidently, the outcomes of the depicted curves were quite satisfactory. In Figure 3, the mean elementary effects are illustrated for each variable of the transmission-line problem for the first case. As already mentioned, the stochastic inputs with high mean elementary effects at a given point in the grid are considered important at that position. In Figure 4a,b, the mean value and the standard deviation of the E field are depicted for the second case. For this scenario, the MC approach required 77 s, while the traditional PC algorithm and the proposed scheme took about 65 s and 38 s, respectively.   The second problem examined the wave propagation within a two-dimensional (2D) space with six concentric dielectric cylinders of infinite length, whose characteristic parameters are shown in Table 2. The computational domain was discretized into 440 × 440 cells, with a spatial density of 20 cells per wavelength in the vacuum. In this case, a sinusoidal wave was used as a source at 2 GHz, emanating at the center of the domain for 14.142 ns. The dielectric permittivities followed again a uniform distribution, in the range between ±5% of their corresponding mean values. In Figure 5, the geometry of this problem is displayed, where the random materials were positioned apart by 0.075 m. The first-order Mur's boundary condition was applied in this case as well. Furthermore, this problem was re-examined, with random inputs ranging between ±10% of their average values. Figure 6a,b depict the mean and the standard deviation of the magnetic field intensity for the "5% percent" case, while Figure 7a,b illustrate the same quantities for the "10% percent" scenario. The results of the proposed method presented a good agreement in this problem as well, as only slight differences were observed, compared to the MC and the standard PC solutions. The simulation time for the MC realizations was approximately 3.28 h (1000 simulations), while the original PC scheme required 1.30 h (389 simulations). However, the proposed PC method needed 44 min (137 runs from the Morris method and 84 runs from the PC scheme with a total of 221 simulations); thus, a speedup of 4.47 compared to the MC method was achieved. Table 2. Mean dielectric permittivities for the 2D problem.   The third problem consists of a patch antenna, the effects of the stochastic inputs on the reflection coefficient (the ratio of the reflected wave to the incident wave at the feeding point of the antenna [23]) of which are examined [24]. In Figure 8, the geometry of the antenna is illustrated, while Table 3 depicts the statistical properties of its random quantities. Those quantities followed a uniform distribution. Furthermore, the antenna was excited via a waveguide port, placed at the edge of the microstrip. In this case, the boundaries were terminated via a Perfectly-Matched Layer (PML) [21]. The dimensionality reduction was performed using the three most important variables, which according to the Morris method were the variables W, L, and the permittivity of the dielectric substrate . Table 3. Mean values and standard deviations for the parameters of the patch-antenna problem. In Figure 9, the mean elementary effects of the patch-antenna are displayed, where it is inferred that the aforementioned random variables were the most influential. Furthermore, the mean value, as well as the standard deviation of the reflection coefficient are depicted in Figure 10a,b, respectively. Specifically, it can be concluded that the mean resonance frequency of the patch-antenna was around 1.8 GHz. In Figure 11, the cumulative distribution function (the probability of a stochastic input taking a given value or less [25]) of the reflection coefficient is illustrated around the mean resonance frequency. The presented outcomes displayed a satisfactory agreement; thus, similar results can be extracted with much less computation time. The required simulation time for the traditional PC scheme was around 6 h (389 simulations), and the proposed approach needed 2.33 h (69 simulations from the PC technique and 84 from the Morris method). The MC realizations lasted 15.33 h (1000 simulations). In conclusion, a speedup of 2.57 compared to the standard PC technique has been achieved. The patch-antenna problem was examined again with three additional random inputs: the substrate width W sl (mean value: 102 mm, standard deviation: 8.67 mm), the substrate length L sl (mean value: 76 mm, standard deviation: 4.81 mm), and the substrate height h (mean value: 4.5 mm, standard deviation: 0.01 mm). In this test case, four random variables were considered stochastic, the parameters W, L, , and L sl . The remaining ones were treated as deterministic. Figure 12a,b illustrate the mean and the standard deviation for this scenario. In this case, the simulation time of the MC method was about 10 h, while the conventional PC scheme required approximately 18 h. However, the proposed PC approach needed 4 h (137 simulations from the PC technique and 120 realizations from the Morris method); thus, a speedup of 2.5 compared to the MC scheme was achieved.

Conclusions
A sensitivity analysis algorithm has been implemented, in order to reduce the computational cost of the PC scheme. The selection of the most important random variables in a given problem can be performed with the proposed heuristic, which is based on the Morris method. The numerical outcomes prove the reliability of the described approach; hence, the efficiency of the PC method can be increased. As future work, a dimensionality reduction of the PC expansion can be performed by combining the Morris method along with anisotropic index sets. Specifically, the indices corresponding to the influential random inputs are more significant than the ones of the trivial stochastic variables. Consequently, the high-order bases, which correspond to the negligible random variables, can be neglected; therefore, the accuracy and the efficiency of the PC method can be further improved.

Conflicts of Interest:
The authors declare no conflict of interest.