Adaptive Polynomial Chaos Expansion for Uncertainty Quantification and Optimization of Horn Antennas at SubTHz Frequencies

Consider a QoI $u$ that depends on a number $d$ of independent random variables included in the vector $\bm{\xi}\triangleq\left[\xi_{1},\xi_{2},\ldots,\xi_{d}\right]\in\mathbb{R}^{d}$. The PCE surrogate model of this QoI $u$ has the following mathematical form [11]:

where each $n$-th term ($n=1,2,\ldots$) of the series depends on a distinct feasible vector $\bm{\alpha}_{n}\triangleq\left[\alpha_{1n},\alpha_{2n},\ldots,\alpha_{dn}\right]\in\mathbb{N}_{0}^{d}$ and $c_{\bm{\alpha}_{n}}$ represents the PCE’s deterministic coefficients of each $n$-th basis function $\Psi_{\bm{\alpha}_{n}}$, which is defined as follows:

According to this expression, each $\alpha_{in}$ component ($i=1,2,\ldots,d$) of each $n$-th feasible $\bm{\alpha}_{n}$ indicates the degree of the univariate function $\psi_{\alpha_{in}}\left(\cdot\right)$. In this paper, we assume that $\xi_{i}$’s are uniformly distributed and the functions $\psi_{\alpha_{in}}\left(\cdot\right)$’s are Legendre polynomials [20]. This assumption has been proved to lead to an exponentially convergent PCE expansion [11]; it is though noted that there exist different choices of the latter parameters that can guarantee exponential convergence. In addition, the set containing the multiple indices $\bm{\alpha}_{n}$’s is chosen such that it forms the anisotropic-hyperbolic index set:

where $\bm{w}\triangleq[w_{1},w_{2},\ldots,w_{d}]\in\mathbb{R}^{d+}$ is a weight vector that induces anisotropy (i.e., polynomials combined with larger weights are not favored) [21], $q\in(0,1]$, and $p$ is the maximum of the polynomial orders present in the PCE series. It can be observed that, when $\bm{w}=\bm{1}$, the weighted-anisotropic index set becomes the more commonly used, but less robust, hyperbolic index set $\mathbb{A}_{q}^{d,p}$ [13]. The hyperbolic truncation scheme reduces the number of polynomials used in the PCE expansion since higher-order polynomials are penalized. This choice works well for PCE series and it is accurate in most cases. As it is known in the literature, the latter holds due to the fact that interactions of low-order polynomials dominate the description of many practical problems.

The weights in the $\mathbb{A}_{q,\bm{w}}^{d,p}$ definition, which favor anisotropy in the selection of the polynomial order, provide an additional degree of freedom that makes the index truncation scheme more robust. This is achieved by relating these weights with the total Sobol indices $S_{i}^{T}$’s. These indices are defined as the fraction of the variance of each random variable $\xi_{i}$ over the total variance $\sigma^{2}\left(u\right)$ of the QoI $u$ [22]. For the considered PCE case, each $S_{i}^{T}$ can be obtained as a function of the PCE deterministic coefficients, i.e.:

where ${\rm w}\left(\bm{\xi}\right)$ represents the probability distribution function (PDF) of $\bm{\xi}$ and each $n$-th series term in (4) depends on $\bm{\alpha}_{n}$ with $\alpha_{in}>0$ $\forall$$i$. In addition, each $w_{i}$ component of the weight vector $\bm{w}$ can be chosen as a function of $S_{i}^{T}$’s according to [21], yielding the following expression:

The logic behind this update is that, as $S_{i}^{\left(T\right)}$ increases, $w_{i}$ decreases, and most likely, the polynomial indices corresponding to $\xi_{i}$ will not be eliminated after applying the anisotropic truncation scheme of (3).

It is essential for the PCE method to obtain the $c_{{\bm{\alpha}}_{n}}$’s coefficients. As it is well-known from [23], the sparsity of effects principle holds for this method, which states that most of these coefficients are zero or almost zero. This implies that, if $\mathbf{c}$ is the vector containing all $N$ PCE coefficients, i.e., $\mathbf{c}\triangleq\left[c_{{\bm{\alpha}}_{1}},c_{{\bm{\alpha}}_{2}},\ldots,c_{{\bm{\alpha}}_{N}}\right]\in\mathbb{R}^{N}$, it should be sparse, and hence, CS methods can be deployed to estimate it [12].

Let the real-valued $M\times N$ dictionary matrix $\mathbf{\Phi}$ with its $(m,n)$-th element ($m=1,2,\ldots,M$ and $n=1,2,\ldots,N$) defined as $[\mathbf{\Phi}]_{m,n}\triangleq\Psi_{{\bm{\alpha}}_{n}}(\bm{\xi}^{(m)})$, where $\bm{\xi}^{(1)},\bm{\xi}^{(2)},\ldots,\bm{\xi}^{(M)}$ represent $M<N$ realizations of the random vector $\bm{\xi}$. The following CS problem can be formulated:

where $\mathbf{y}\triangleq\left[y(\bm{\xi}^{(1)}),y(\bm{\xi}^{(2)}),\ldots,y(\bm{\xi}^{(M)})\right]$ (also called measurement vector) and $\left\|{\mathbf{c}}\right\|_{0}$ provides the number of non-zero elements in $\mathbf{c}$. Clearly, the solution of (6) is the sparsest vector ${\mathbf{c}}^{\ast}$. Finding it though, has been shown to be infeasible since a search through all possible support sets needs to be performed. In this paper, we adopt the OMP method to approximate (6)’s solution [12]. OMP is an iterative greedy algorithm aiming to find the columns of $\mathbf{\Phi}$ participating the most in the construction of the measurements contained in $\mathbf{y}$. In every OMP iteration, a column-index set $\mathcal{A}$ is augmented and a residual $\mathbf{r}$ vector is calculated. The algorithm is initialized with the residual $\mathbf{r}^{(0)}$ set equal to $\mathbf{y}$ and the empty set $\mathcal{A}^{(0)}$. During each $k$-th iteration, $[\bm{\Phi}]_{:,n}^{\mathrm{T}}\mathbf{r}^{(k)}/\left\|[\bm{\Phi}]_{:,n}\right\|_{2}$ ($[\bm{\Phi}]_{:,n}$ is the $n$-th column of $\bm{\Phi}$) is calculated, indicating the correlation of each column of $\bm{\Phi}$ (excluding those already included in $\mathcal{A}^{(k-1)}$) with the residual, and the index corresponding to the maximum correlation is included to $\mathcal{A}^{(k)}$. Then, the new estimate of the sparse solution vector is obtained from the solution of the following least-squares problem:

where ${{\mathbf{\Phi}}_{{\mathcal{A}^{(k)}}}}$ is the submatrix of $\mathbf{\Phi}$ that contains the columns with indices in $\mathcal{A}^{(k)}$. The last step of each $k$-th iteration is to update the residual as follows:

The number $s$ of non-zero coefficients in $\mathbf{c}$ is pre-selected and the iteration index is set as $k=1,2,\ldots,s$. Alternatively, OMP iterations may terminate when ${\left\|{\bm{\Phi}{\mathbf{c}}^{(k)}-\mathbf{y}}\right\|_{2}}$ becomes smaller than a given threshold $\epsilon$.

The estimation error of the latter OMP-based approximation, $\epsilon_{\mathbf{c}^{\ast}\left(p\right)}$, can be calculated via the leave-one-out cross validation method [12]. In particular, this error is given by:

where $X$ denotes the set with the random sample points $x^{\left(i\right)}$ with $i=1,2,\ldots,N_{p}$ (i.e., the set used in the experimental study), where $N_{p}$ is the number of these points in the PCE of maximum polynomial order $p$, and $L\left(M_{X,p}\right)$ is the leave-one-out cross-validation error, which is calculated as follows [23]:

In this definition, $QoI\left(x^{\left(i\right)}\right)$ is the QoI corresponding to the $x^{\left(i\right)}$ input, while $M_{X}\left(x^{\left(i\right)}\right)$ indicates the output of a PCE that is built from all the measurements $x^{\left(i\right)}\in X$. In addition, $h_{i}$ represents the $i$-th diagonal term of $\bm{\Phi}\left(\bm{\Phi}^{\mathrm{T}}\bm{\Phi}\right)^{-1}\bm{\Phi}^{\mathrm{T}}$. Finally, the correction term $CT\left(p,N_{p}\right)$ in (9) is designed to reduce the sensitivity of $L\left(M_{X,p}\right)$ in the overfitting regime (i.e., when $p$ increases), and is defined in [24] as follows:

The steps of the adaptive PCE method are summarized in Algorithm 1.

To enhance the sparsity of the PCE coefficients and, thus, require less sampling points $G$ for a successful reconstruction, the PCE coefficients are found via a weighted $l_{0}$-minimization [25]. In particular, the optimization problem (6) is rewritten in the following form:

where the matrix $\mathbf{W}$ is diagonal having the elements $w_{i}$ $\forall$$i=1,2,\ldots,N$ on its main diagonal. These weights are considered as free parameters which are chosen to improve the PCE coefficients reconstruction. In particular, if the solution of (6) is solved by a particular algorithm, i.e. the OMP, then its corresponding weighted counterpart (12) can be solved using $l_{max}$ successive applications of OMP, where at the $j$-th iteration the weights can be chosen as follows:

where $j=1,2,\ldots,l_{max}$, $i=1,2,\ldots,N$, and $\epsilon_{w}$ is a small number used to avoid division by zero. This expression eventually forces to zero the small-valued coefficients in the recovered $\mathbf{c}^{\ast}$ vector [25], thus, enhancing sparsity in the solution. This fact leads effectively to a smaller number of required measurements $M$ for accurate reconstruction. The efficiency of this approach, which we henceforth call as weighted OMP (WOMP), will be validated through the uncertainty quantification analysis of a horn antenna in Section IV.

By using the PCE coefficients $\mathbf{c}_{f}^{\ast}$ computed via Algorithm 1, the PCE surrogate model for the QoI $u$ can be built according to (1), as follows:

It can be easily observed that the right hand side (RHS) of this expression is actually a sample generator of $u$, which, according to the PCE theory, can be exponentially convergent to the true $u$ values as more terms are considered in the PCE. Interestingly, the generation of $u$-samples via (14) is highly computationally efficient, since these samples are obtained through an analytical formula. Hence, (14) can be used to trivially generate large amounts of QoI samples, which can be further used to approximate the PDF of $u$, and consequently, any order of its statistical moments. Additionally, this expression can feed with large numbers of QoI samples any algorithm dealing with the QoI optimization over the design variables of a system, e.g., an antenna. This constitutes a great deal of computational saving, which would otherwise need, in the case of an antenna design, a full wave solver for the QoI generation, making the optimization most likely infeasible.

The proposed adaptive PCE method has been applied to compute the PCE coefficients via (14), which were first used to estimate the mean value and standard deviation (std) of the QoI $G$, respectively, as follows:

where $\mathbf{0}_{d}$ is a vector with $d$ zeros and $\mathcal{N}$ denotes the cardinality of the set $\mathbb{A}_{q,\bm{w}}^{d,p}-\{\mathbf{0}_{d}\}$ where all $\bm{\alpha_{n}}$’s belong to.

In Fig. 3, the parameter $\mu_{G}$ is illustrated as a function of the angle of the incident plane wave at the frequency $95$ GHz, which was fixed for all numerical investigations in this paper. The expression (15) was evaluated for the PCE method using $M=25$ independent $G$-samples and the maximum polynomial order $p_{\rm max}=4$, while, for the MC method, this metric was calculated directly via the mean value definition using $1000$ independent samples (almost two orders of magnitude more samples than PCE). As observed in the figure, there is excellent agreement between the two methods. It is noted that, to obtain every sample of $G$, a run of a full-wave EM finite element solver is required, which is significantly time-consuming, and thus, constitutes a major computational bottleneck. This figure actually showcases the power of the PCE method, which can serve as an accurate uncertainty quantification method that is much faster than the reference MC. As stated previously, for the optimal PCE convergence rate, the PCE random variable design samples need to follow the uniform distribution, while the MC samples are generated via the lattice hypercube method, where computation errors for $M$ samples reduce at a rate $1/M$ (which is better than the standard $1/\sqrt{M}$ rate of the MC method when samples are generated randomly).

Parameter $\sigma_{G}$ is plotted versus the angle of the incident plane wave in Figs. 4 and 5 for the two investigated methods. In Fig. 4, the same PCE parameter setting with Fig. 3 has been used, while Fig. 5 includes PCE performance curves for $p_{\rm max}=5$ and $M=50$ independent antenna gain samples. It can be seen that there is a very good agreement between the two methods, with this agreement being improved in Fig. 5, where a slightly higher maximum polynomial order and slightly more $G$-samples have been used.

According to the presented uncertainty quantification method in Section II, once the PCE coefficients are obtained, the Sobol indices can be calculated via expression $\eqref{TS}$ for every uncertain design variable. In Fig. 6, the total Sobol index of the most uncertain variable $a_{0}$ (in the sense it varies the most with respect to its nominal value) is depicted as a function of the angle of the incident plane wave. This index quantifies, for every angle of incidence, the percentage of the $a_{0}$ variable in the total variance of the antenna gain $G$. Similarly, the total Sobol indices for all other design variables can be calculated, but are omitted here due to space limitations.

To improve the performance of the proposed PCE method in Algorithm 1, we have replaced the OMP application in the corresponding algorithmic step (Step $4$) with the WOMP approach presented in (12) and (13). The corresponding performance evaluation is included in Fig. 7 for $M=20$ independent $G$-samples and the same maximum polynomial order $p_{max}=5$, where it can be observed that the WOMP solver reconstructs successively the PCE coefficients, generating average $G$-values that are in very good agreement with corresponding ones obtained via MC simulations. On the contrary, it is depicted that, for the same value of $M$ and maximum polynomial order, OMP fails. However, as it can be seen from Fig. 5, when the number of independent $G$-samples increases to $M=50$, OMP yields successful reconstruction. This performance gap between OMP and WOMP can be explained from the fact that the weighted $l_{0}$-minimization algorithm in the latter is capable to increase the sparsity level of the reconstructed solution and, thus, WOMP requires smaller number of samples to achieve the same task.

A considerable number $N_{K}$ of samples for the horn antenna gain variable $G$ can be efficiently generated through the proposed adaptive PCE method, via (14), at a negligible computational expense. Subsequently, these samples can be directly used to construct a histogram of the $G$-samples, coupled with the application of a smoothing kernel, yielding an approximation for the PDF of variable $G$. In particular, the PDF of $G$, $f_{G}(\cdot)$, can be approximated as follows:

where $K\left(\cdot\right)$ represents a positive-definite function known as the kernel, while $h_{K}$ serves as the bandwidth parameter. A Gaussian kernel, which is characterized by the standard normal PDF, can be chosen, and for this kernel, the optimal bandwidth is given by the formula:

where $\hat{\sigma}_{R}$ and $\widehat{iqr}_{R}$ denote the empirical standard deviation and interquartile range, respectively. As highlighted in [23], the optimality of $\hat{\sigma}_{R}$ lies in its ability to minimize the $\ell_{2}$ approximation error $\|f_{G}(y)-\hat{f}_{G}(y)\|_{2}$.

In Fig. 8, the PDF calculated from PCE-based $G$-samples is plotted for the angle $\theta=0^{\circ}$ of the incident plane wave. It is noted that the PDF for any other angle value $\theta$ can be obtained in a similar way. In fact, in this way, any order probabilistic moment of $G$ can be calculated. This happens because expression (14) is actually a $G$-sample generator which can be also used for optimization tasks in conjunction with an appropriately designed optimization algorithm.

Another pertinent statistical metric is the $P$-th percentile of reflection with $0\leq P\leq 100$, which can be efficiently computed via the proposed PCE method. By denoting each $P$-th percentile of $G$ as $G_{P}$, this metric signifies that $P\%$ of the observations of $G$ are smaller than $G_{P}$. After determining the PCE coefficients $\mathbf{c}_{f}^{\ast}$ in (14), a sufficiently large number $N$ of $G$-samples can be generated analytically with negligible computational cost. Subsequently, the $P$-th percentile $G_{P}$ of $G$ can be computed using the nearest-rank method. According to this method, the gain samples $G_{i}$, with $i=1,2,\ldots,N$, are first arranged in ascending order to form the set of sorted reflection observations $G_{s,i}$ with $i=1,2,\ldots,N$. Consequently, $G_{P}$ is determined as $G_{s,n}$ where $n$ is the smallest integer greater than or equal to $PN/100$. The $5$-th and $95$-th percentiles of $G$ as calculated via the proposed adaptive PCE method are depicted in Fig. 9.

The proposed adaptive PCE method via expression (14) has been used as a generator of the QoI $G$-samples with the goal to maximize the gain of the horn antenna in Fig. 1. For this optimization, we have considered the PSO method [19] with the following cost function for a particle’s solution $\mathbf{x}$:

where $G_{\text{I}}$ denotes the ideal antenna gain and the angle regimes $R_{m}$ with $m=1,2,3,$ and $4$ defined as $R_{1}:-180^{\circ}\leq\theta<-100^{\circ}$, $R_{2}:-100^{\circ}\leq\theta<-20^{\circ}$, $R_{3}:-20^{\circ}\leq\theta<20^{\circ}$, $R_{4}:20^{\circ}\leq\theta<100^{\circ}$, $R_{5}:100^{\circ}\leq\theta\leq 180^{\circ}$, which was combined with PCE, as described in the Appendix. This optimization approach was then deployed for the antenna gain maximization and the obtained optimum design variables are included in Table I. This table also lists the optimum design variables of the scheme presented in [28]. As observed, the values for all antenna parameters with the proposed PSO-PCE design are very close to the ones via the benchmark scheme.

Figure 10 illustrates the normalized horn antenna gain $G_{\text{PCE}}$ for the optimal design variables in Table I with the proposed PSO-PCE approach in comparison with the $G_{\text{I}}$ values. It can be observed that the flat top of the horn antenna with the proposed PSO-PCE approach within the range of angles $-50^{\circ}\leq\theta\leq 50^{\circ}$ closely aligns with the ideal case $G_{\text{I}}$. More importantly, this alignment is nearly identical to the gain depicted in Fig. 2. In addition, it needs to be noted that the proposed PCE-based optimization exhibits significantly higher efficiency, given that the gain $G$-samples (a few thousands in number) provided to the optimizer were analytically computed using (14), as opposed to relying on a full-wave EM solver (e.g., HFSS). The latter takes approximately $5$ minutes on an average personal computer of average computational capabilities. In the proposed PSO-PCE antenna gain optimization approach, the PCE coefficients were determined using the adaptive PCE Algorithm 1, utilizing 100 $G$-samples which were obtained via HFSS for uniformly distributed design variables with a 10$\%$ variation around the parameters specified in [28].

In Fig. 11, the performance of the proposed PSO-PCE antenna gain optimization approach is investigated for the case where the $G$-samples are generated solely via the proposed PCE model. In particular, the samples encompassed variations of $10\%$ or $60\%$ around the parameters specified in [28]. As depicted, for the case of $60\%$ variation, the gain obtained from the proposed optimization closely aligns with the ideal one. This underscores the robust generalization properties of the proposed PCE model as a generator of $G$-samples, given its construction based on samples with only a $10\%$ variation around the reference parameters from [28]. Finally, in Fig. 12, the same results are demonstrated for the case where the $G$-samples are generated via the proposed PCE model in (14), and in particular, for $M=20$, $60$, and $100$ $G$-samples. It can be observed that the most accurate (closer to the ideal case) model in the region of interest $[-50^{\circ},50^{\circ}]$ (flat-gain) is the model created by $100$ $G$-samples. Notably, even with as few as $M=20$ samples, the resulting PCE model, coupled with the PSO algorithm, yields an optimized gain that is remarkably close to the ideal one. Overall this behavior underscores the impressive computational efficiency of the proposed PCE method in optimization endeavors.