Residue-Pole Methods for Variability Analysis of S-parameters of Microwave Devices with 3D FEM and Mesh Deformation

This paper presents a new approach for variability analysis of microwave devices with a high dimension of uncertain parameters. The proposed technique is based on modeling an approximation of system by its poles and residues using several modeling methods, including ordinary kriging, Adaptive Polynomial Chaos (APCE), and Support Vector Machine Regression (SVM). The computational cost is compared with the traditional Monte-Carlo method. To improve the efficiency, mesh deformation is applied within 3D FEM framework.


Introduction
One of the main challenges in designing high frequency microwave filters is predicting the impact of the physical parameters and geometrical variability on their performance in micro-scale fabrication. The rapid development of 3Dprinting technologies gives one great flexibility on device geometry but also is connected with limited accuracy of the manufacturing process. For example, Selective Laser Melting (SLM) technology can be used in 3D printing of RF components, while covering the wide range of materials that can be processed, including aluminum, titanium, and steel alloys the manufacturing accuracy range varies from 0.02 mm to 0.1 mm depending on the specific material used [1]. This variability is an unavoidable part of the manufacturing process and it leads to slight changes of the dimensions of the physical devices, which brings uncertainty of the response of the microwave components. Traditionally, Monte-Carlo (MC) simulations have been used in the commercial circuit and electromagnetic simulations for predicting the statistical distribution of the component and system-level performance [2][3][4][5]. However, the slow convergence for the MC has become a computational burden, especially in simulations of large and complex circuits. This fact has prompted wide interest in exploring alternative approaches to the problem of statistical analysis of the performance of electronic devices. Recently some intrusive methods using generalized Polynomial Chaos (PC) have been reported to overcome MC problems [6][7][8]. Intrusive methods generally lead to a coupled system of equations by modification of the system fundamental equations, which is very time consuming to solve with respect to the original problem when the dimension of uncertain parameters is high [9]. Alternatively, there are nonintrusive methods that treat the original system of equations as a black-box, and there is no need to change the fundamental system equations. As a result, they are easier to implement, comparing to intrusive methods [9]. For problems with higher dimensions, non-intrusive hierarchical sparse PC approaches have been used, which the PC expansion terms are reduced based form of conventional full-blown PC expansions in general [10][11][12][13][14]. For eliminating the computational cost of model construction due to the large size of PC expansions, in [14], a new approach has been presented, which ends up in lower PC expansion with keeping the model accurate. However, the technique needs a large number of training samples for the PC black-box model to be trained. In [15] Least Square Support Vector Machine method (LSSVM) for modeling variability analysis of a system with a large number of uncertain parameters has been adopted, while the complexity of the black-box model is generally lower than PC based methods, it still needs a large set of training samples. Usually for creating a surrogate model of microwave circuits the values of the device's scattering parameters at selected frequency points are considered [16], [17]. Recently, several approaches based on the characterization of rational models have been presented [18][19][20]. In [19], neural networks are trained to learn the relationship between residues and poles of the rational model approximating scattering matrix and the geometrical parameters. In this approach, still many fullwave response samples were needed to train the system to achieve reasonable accuracy. In [20] the efficiency of creating the surrogate model using the rational representation of the system has been proved for a large amount of deviation of uncertain parameters with a smaller set of training sam-ples. However, the process of building the model is complex, because the approach considers the zeroes of the system in addition to the poles and residues as well. Moreover, the process of vector fitting is more complex than what is proposed in this paper. In [20] to obtain more accuracy for surrogating, the order of vector fitting is chosen to be higher than needed. Therefore, the insignificant poles should be discarded, which makes the process of building the model of the system more complicated. The central issue in applying the non-intrusive variability analysis methods to the complex microwave systems is the computational cost of generating proper samples for training the black-box model of the system. Usually, the training samples are generated using an Electromagnetic (EM) simulator, and the cost of generating every single sample (that corresponds to single EM simulation) is very high, especially for a device with a large number of inputs or random parameters. In the case of microwave filters, the 3D finite element method in frequency domain is commonly used as a simulation tool. Therefore, it is essential to train the black-box model with a lower number of training samples to make the complete uncertainty analysis faster. Moreover, the complexity of building the model with respect to its accuracy is another topic of interest.
In this paper, the training samples for creating the blackbox model of the system are generated by InventSim, a fullwave 3D FEM EM simulator [21]. The simulator was used to calculate scattering matrix parameters of a microwave filter for a given set of frequency points and a set of slightly disturbed values of its geometrical parameters. A set of poles and its corresponding residues of the calculated scattering matrix was then used as the training samples required to build a surrogate model of the original system. There are many toolboxes which allow constructing such a model, e.g. Polynomial Chaos, kriging, Support Vector Machine Regression. Moreover, there are several MATLAB toolboxes available for this purpose, such as UQLAB [22], LSSVM [23] and SUMO [24]. In this paper, the models are computed using UQLAB Toolbox in MATLAB by applying kriging method [25]. Then they are compared with Adaptive Sparse PC [26] and Support Vector Machine Regression [27] methods using UQLAB toolbox with respect to model accuracy and computational cost. If the output is the component response in the form of the scattering matrix, a simple way to do the variability analysis is to model the magnitude of scattering parameters at each frequency point. To build a high-quality model of the microwave filter, one may need to provide the response for many frequency points. However, in such an approach (hereinafter called the Direct method), doing the modeling for many frequency points is time-consuming. In this paper, an advanced technique for prediction of the behavior of complex systems, by modeling their residues and poles, will be discussed. It will be shown that by using only a few tens of training samples, a low-cost model can be built and trained to follow the behavior of the black-box microwave system with a large number of uncertain parameters.
The remainder of the paper is organized as follows: Section 2 presents mathematical background on different surrogating techniques. Two different approaches for modeling an approximation of system are described in Sec. 3 and 4. In Sec. 5, two different ways of sampling using InventSim software are discussed. The concepts introduced are then verified and compared together by experiments in Sec. 6, which is followed by the conclusion in Sec. 7.

Kriging
Kriging is a well-known interpolation technique, which is used in the deterministic simulation. Let us assume that variability analysis will be performed by changing d physical and geometrical parameters of the microwave device and that properties of this device are represented by a scattering matrix S. 1 The black-box model of this device can be built using K training samples. i th training sample (x i , y i ) is a pair of an input vector x i and output vector y i . Input vector x i = [x 1 , x 2 , ..., x d ] contains values of the mentioned d parameters. Output vector x i for a Direct method contains scattering matrix parameters for all required frequency points and for a Residue-Pole method y i it contains values of the residues and poles. The system can be modeled using the basic type of kriging technique, which is named ordinary kriging. Mathematical description of ordinary kriging includes the definition of correlation matrix ψ(x, x * ; θ). This matrix is symmetric and positive semi-definite and represents the correlation between vector of observation points x (training data points) and vector of prediction points x * (evaluation data points). The size of this matrix is the K × K. For the purpose of this article the correlation family is chosen to be Gaussian, therefore the correlation matrix is defined by (1), where θ is the optimal hyper-parameter [25], exp is exponential function. Correlation matrix defined in (1) is used in the definition of the ordinary kriging equation (2), [25], [28].
where 1 is a K-dimensional vector with each element being the value 1, β is the constant mean function, γ -the K-dimensional vector of cross-correlations (similarities) between the prediction point and each one of the observation points, A T denotes the transpose of matrix A.

Adaptive Sparse Polynomial Expansion
This section presents a quick overview on the advanced sparse PC expansion concepts and mathematical descriptions. Considering the PC expansion of the model can be written as where a ζ are the unknown deterministic coefficients and φ ζ are the multivariate polynomials which in this paper they are considered as Hermite polynomial. In the normal expansion the number of total terms is, where d is the dimension of the input vector x and h is the highest degree of polynomial [29]. In the case of high dimensional input, the number of PC terms increases very fast and it is almost impossible to use all the terms due to the large size of the model expansion. Usually in PC expansion all of the basis are not needed and in order to make the expansion more efficient in high dimensional problem, the basis of the expansion could be reduced by hyperbolic truncation (6), where the q is the hyperbolic term and 0 < q < 1. In this method the degree of polynomials will be increased adaptively and regression method is chosen to be Least Angel Regression (LARS). For detailed information a reader may refer to [26], [29]

Support Vector Machine Regression
To approximate a problem provided by a generic blackbox system, containing a set of K training data the system can be modeled by using the following nonlinear SVM regression: where w is a vector of weight coefficients and b is an offset parameter. To find the best combination of the parameters (w, b) in (7), the following loss-function should be minimized by: For minimizing (9), the below optimization problem can be solved, subject to where ξ i and ξ * i are slack variables, which measure the deviation from the insensitive tube and C is a regularization parameter, which provides a trade-off between the accuracy of the model and its flatness. The concept is shown in Fig. 1, readers may refer to [27], [15] for more detailed information.

Standard Approach: Direct Method
With the filter optimization parameters as an input of the system, if the desired output is the magnitude of scattering parameters, the training input matrix of the model is defined as follows: where, gi, j -value of the j th optimization parameter for th training sample, i ∈ {1, 2, ..., K } and j ∈ {1, 2, ..., d}, K -number of training samples and number of rows of the input matrix, d -number of physical and geometrical parameters and number of columns of the input matrix.
In Direct-method, the output of the system should be evaluated at each frequency point. The training output is defined as a matrix Y K×[N S ·N f ] Direct in (13).
where, N -the number of rows of the scattering matrix S, N S -the number of S parameters in the scattering matrix, N S = N 2 , f p -the p th frequency point, N f -the number of frequency points, SIJ -the elements at I th row and J th column of the scattering matrix S, |SI J k, f p | -the magnitude of SI J at k th sampling point and f p frequency. Figure 2 shows a block diagram representation of relation between input and evaluated output for a single training data.

An Example for Relation between Input and Output in Direct Method
For example, if a scattering matrix is S = S11 S12 S21 S22 and N f = 501, then N = 2, and N S = 4. For a set of 40 Training (each sample is represented by a single row in the X Training matrix), there will be N S · N f = 4 · 501 = 2004 number of independent outputs for the Direct-method for every set of training samples (the number of columns of Y Direct would be equal 2004 for this example in (13)).
Block-Diagram representation of input and evaluated output relation of Direct-method for a single set training data.

Proposed Residue-Pole Method
Scattering parameters can be represented as rational function of the poles and the corresponding residues of the system. For given frequency response in frequency domain, the residues and poles can be calculated using Vector Fitting method [30][31][32]. In this approach, each element of the entire scattering matrix S is approximated as (15), where, p i -i th pole of the system, r I J i -i th residue of the SI J, s -the complex frequency, O VF -the order of vector fitting, which means there are O VF 2 pairs of complex conjugated residues and poles.
In the variability analysis using Residue-Pole method (RP), the set of residues and poles calculated using Vector Fitting method defines an output matrix Y RP of the model of the system. Instead of modeling the absolute value of S-parameters directly for all frequency points, the values of residues and poles are going to be modeled. In this approach, first, a few tens of full-wave responses are generated using an EM simulator, then, by initializing O VF =2, the process of vector fitting starts. The order of vector fitting is increased adaptively till the accuracy of the rational representation of the scattering parameters is acceptable, and with the minimum amount of O VF the process of vector fitting stops to avoid modeling of unnecessary poles and residues. Considering the same input for both models, the size of the output matrix of the Residue-Pole method, defined in (16), is much smaller in comparison with the Direct method because there is no need to model the system at each frequency point. The size of the output matrix of the Residue-Pole method is equal K × [O VF · (N S + 1)], and it is not dependent upon the number of frequency points. with . . .
where, Re[...] and Im[...] denote the real and imaginary part of a complex vector or matrix, r I J k -a row vector of complex residues of the SIJ for the k th training sample (length of the vector is O VF 2 ), R I J -a matrix of complex residues of the SIJ for all training samples (size of the matrix is K × O VF 2 ), p k -a row vector of complex poles which is common to all scattering parameters calculated for the k th training sample (length of the vector is O VF 2 ), N is the number of rows of the scattering parameter matrix S.
Size of Y RP depends on the order of Vector Fitting, O VF . In this method, O VF is always an even number because the system has O VF 2 pairs of complex conjugate poles and residues. Because of this property model can be built using only half of the number of residues and poles (only one complex value from a complex pair is required). In the result, the dimension D YRP of the Y RP is defined in (22). Multiplication by two in (22) denotes number of variables in the complex number. Dimension of the Y RP is, Simplification of (22) is defined in (23). Figure 3 shows a block diagram representation of relation between input and evaluated output of Residue-Polemethod for a single set of input.In the next step, the state space equation can be calculated for all frequency points at once. Therefore, the entire scattering matrix S can be defined in a matrix equation (24). It should be noted that the in (24), E and D matrices of state space equation definition are neglected.
where, A -a diagonal matrix of complex conjugated poles of the transfer function ( size of the matrix is (N · O VF ) × (N · O VF )), C -a matrix of corresponding residues of the transfer function ( size of the matrix is N × (N · O VF )), s -the complex frequency, I -the identity matrix ( size of the matrix is (N · O VF ) × (N · O VF )). B -port selector matrix with value 1 for each port (size of the matrix is (N · O VF ) × N).

Sampling Method for Generation of a Training Set
There are several methods for the generation of training set information. In this paper, the source of design uncertainties is introduced by selected design parameters of the filter structure. The variability is defined as the functions of normalized Gaussian random variables, and the samples are gathered by the Latin Hypercube Sampling method (LHS). The training set is generated by running simulations in In-ventSim software. Inputs of the system are geometrical parameters for a filter, and the outputs are S-Parameters that result from the full-wave simulation. There are two ways for perturbation simulation using InventSim. The first method is setting new values of inputs and creating the mesh for the filter from scratch, then calculating the desired output. By repeating the simulation (using FOR loop), a desired experimental design is generated. The disadvantage of this method is its computational cost due to creating a complete mesh in every simulation. Another way is to use the mesh morphing technique. In this method, the mesh will not be created from scratch in every simulation, and it will be adjusted with respect to the difference of the mean values of inputs and new values for the next simulation (a new sample for the inputs) [33], [34].

Numerical Examples
In this section, the performance of the Direct and Residue-Pole approaches with three different techniques for modeling is compared for the variability analysis of microwave filter in example 1. It will be shown that ordinary kriging modeling technique has the best accuracy among the other methods. In example 2, the efficiency of the Direct and RP approaches is compared using ordinary kriging modeling method only, which has the best accuracy with respect to APCE and SVM modeling techniques.

Example 1: Waveguide Filter
As the first example, the third order H-plane rectangular waveguide filter with the center frequency of 10.25 GHz and the bandwidth of 500 MHz is presented. The microwave filter is shown in Fig. 4 and has defined d = 8 uncertain geometrical variables, that control the lengths of the resonators and coupling irises.
The nominal values of optimized geometrical parameters of filter are considered to have a deviation by a span of −180 to +180 microns by Gaussian distribution with a standard deviation of σ = 0.06 in millimeters. The simulation has been done for N f = 401 frequency points. Only K = 20 samples are used for training the model of the system using kriging, Adaptive PCE, and SVM methods. The results are compared with N val = 500 validation samples set, which is generated by InventSim by MC method.   Figure 5 shows the Root Mean Square (RMSE) errors of Direct and RP approach using different training methods. Sampling in this example is generated without meshmorphing. Both the Direct and RP approaches are compared with the same input training set and the RMSE is calculated in each frequency point with (25), As can be seen in Fig. 5, the accuracy achieved by the proposed RP approach using ordinary kriging is the best among other methods. Where, |SI J | mean denotes the mean value for N val samples of MC simulation, |RMSE| method−RP shows the RMSE error value of evaluated data by Residue-Pole model using mentioned method with respect to original data and |RMSE| method−D describes the RMSE error value of evaluated data by Direct model using mentioned method with respect to original data.

Comparison of Proposed Residue-Pole Method with Direct and MC Methods
The model cost is about 9 seconds only, and the time needed for generating training samples by InventSim EM simulator is 10 minutes. Therefore, the kriging technique is chosen for training the RP model. In the case of modeling the system using the Direct approach, again, the best accuracy is achieved by kriging while the model cost is about 111 seconds, which is 13 times more than the model cost of RP method because the model size of Direct approach is much bigger than RP approach. Moreover, the accuracy of the RP approach for S11 is five times higher than the Direct approach. In a manner of total computational cost comparison with MC method, the RP approach needs only 3.3 seconds for vector fitting plus 8.7 seconds of the training model and evaluation cost in addition to 10 minutes for generating K = 20 samples for the training set. Therefore, the total cost for variability analysis of N val = 500 validation samples is 10.2 minutes for the RP approach using kriging, while using the traditional MC method needs 250 minutes, so the RP approach is 24.5 times faster than MC method with an excellent level of accuracy. Detailed information of the compared results for entire S parameters of the waveguide filter example-without mesh morphing technique (error is calculated for S11 and S12 parameters only). The detailed information about comparison of different surrogating tech-niques is given in Tab. 1, where, Y D denotes the size of output of experimental design (output training set), X D is the size of inputs of experimental design (input training set), Model cost is the elapsed time for applying the corresponding surrogating technique and train the model based on experimental design, S11 max RMSE and S12 max RMSE are the maximum error calculated by (25) only element of S matrix , Training samples cost describes the elapsed time for generating training samples by InventSim full-wave simulator, and MC cost is the elapsed time for generating N val = 300 MC samples by InventSim full-wave simulator.

Example 2: Combline Filter
In the Sec. 6.1 it was shown that kriging method has significantly higher accuracy than APCE and SVM methods. Therefore, in this section only ordinary kriging method is used to construct the model for both Direct and RP approaches. The structure is a 6th order combline filter with center frequency 1.74 GHz and bandwidth 70 MHz. In the structure, single cross-coupling (quadruplet topology) is introduced that realizes a pair of symmetric transmission zeros. The microwave filter example in Fig. 6 has d = 15 uncertain optimization geometrical variables. The design variables are six tuning screws that control the resonant frequency of the combline resonators, six tuning screws that control couplings between resonators, and remaining design parameters control source/load couplings to first/last resonator. The nominal value of optimized geometrical parameters of filter is considered to have a deviation by a span of −100 to +100 microns by Gaussian distribution with a standard deviation of σ = 0.035 in millimeters. The simulation has been done for N f = 501 frequency points. Only K = 40 samples are used for modeling the behaviour of system and results are compared with N val = 300 validation samples set which is generated by InventSim with MC method.  Tab. 1. Detailed information of the compared results for entire S parameters of the waveguide filter example -without mesh morphing technique (error is calculated for S11 and S12 parameters only).

Sampling without Mesh Morphing
Both the methods are compared with the same input training set. Figures 7-8 show the results of accuracy comparison for S11 and S12 parameters respectively, when the number of training samples is K = 40 for a validation set of N val = 300. Instead of plotting all of the N val = 300 samples, the mean value of evaluated model for each method and the reference data is calculated by using (26) and then are compared together.
As can be seen in the Figs. 7-8, the proposed method has better accuracy. For the Residue-Pole method RMSE error at its most is equal to 0.033 for S11, which is significantly lower than 0.089, the maximum RMSE error calculated for a Direct method. The computational time for building a modeling for the entire S parameters (S11; S12; S21; S22) when the number of training set is K = 40, is about 855 seconds for Direct method, while, for the proposed approach is 37 seconds, which is about 23 times faster in computational time of creating the model. Since the Residue-Pole model has an acceptable accuracy in comparison with the MC method, it can be used instead of the MC method. As long as the Residue-Pole method needs a few tens of samples for the training set for creating the model, it could be compared with the MC method in a manner of computational cost. In this example, for N val = 300, only K = 40 training samples are needed; therefore, 60.06 minutes are required for training set generation and the overhead time for Vector Fitting and kriging is only 37 seconds. Alternatively, 450.45 minutes are needed for variability analysis using the MC method. Therefore, the Residue-Pole method is about 6.75 times faster than the MC method with acceptable accuracy.

Sampling using Mesh Morphing Technique
Further improvement of the performance of the variability analysis was achieved using mesh morphing technique. By applying the mesh morphing technique for perturbation analysis, the time needed for generating each sample is reduced 1.7 times comparing to regular sampling without mesh deformation. Using the mesh-morphing technique, there are improvements in the accuracy of modeling the S parameters by using Residue-Pole-method as well. The cost for creating models for the Direct and Residue-Pole methods remained the same as before, while the time needed for generating training samples reduced significantly. The total time needed for generating N val = 300 samples using MC method is 261 minutes, while for Residue-Pole method, there are just K = 40 samples needed to create a model for predicting the behavior of the system. The time spent for simulating the variability analysis of the combline filter in the numerical example for Residue-Pole method is 34.8 minutes for generating the training set which contains only K = 40 samples plus 37 seconds overhead for building a model due to vector fitting and kriging, which in total will be T total RP = 35.4 minutes. The time for generating N val = 300 samples using traditional MC method is 300×0.87 = 261 minutes. Therefore, with K = 40 training samples, the proposed method is 7.37 times faster than the MC method with an excellent level of accuracy.
The results are shown in Figs. 9-10. Figure 11 shows that the model of the system can be built using RP approach with only K = 20 training samples. The maximum RMSE error of S11 is 0.078, which is still lower than the Direct approach with K = 40 training samples (0.096) and the maximum RMSE error level of S12 is equal to 0.52 which is an acceptable level of accuracy. In this case, the total time for the variability analysis simulation using the RP approach is 17.65 minutes, including vector fitting, model construction, and generating training set samples, while the MC approach needs 261 minutes as mentioned earlier. Therefore, for an acceptable level of accuracy, the RP approach using the mesh morphing technique for sampling is 14.78 times faster than the traditional MC method. Tab. 3. Detailed information of the compared results for entire S parameters using mesh-morphing technique, (error is calculated for S11 ans S12 parameters only).

Conclusions
In this paper, two different modeling approaches, along with two kinds of sampling methods for generating the training samples set for modeling a variability analysis of microwave components are proposed. First it is shown that, among Adaptive PCE, SVM and ordinary kriging, the best method for modeling each approach is ordinary kriging. Thereafter, Direct and RP approaches are compared together in the manner of accuracy and computational cost. In the waveguide filter example, which contains eight uncertain parameters, the RP method has proven its efficiency both in terms of accuracy and computational time. It is shown that the RP approach using ordinary kriging method for modeling the system behavior with only K = 20 training samples is about 5 times more accurate than the Direct approach. In example 2, a more complex filter structure has been adopted for comparing two different approaches. Without using the mesh-morphing sampling technique, by applying the proposed Residue-Poles method, the maximum RMSE errors of modeled S11 and S12 are 3 and 2 times lower than Direct method respectively. Therefore, due to the acceptable accuracy of Residue-Pole method, it is chosen to be compared with MC method. Since it is only K = 40 samples needed to train the model, the proposed RP method is 6.75 times faster than MC method.
When using the mesh-morphing technique for sampling, there is 1.7 times speedup for generating each training sample of the training set, while the Direct method accuracy is decreased, the Residue-Pole method does not loose its accuracy and is even 1 percent more accurate at maximum RMSE error with respect to normal sampling. When using the mesh-morphing technique, the Residue-Pole method is 7.37 times faster than MC simulation with an acceptable level of accuracy. It is even shown that for an acceptable level of accuracy, only K = 20 samples are enough to train the model using the proposed RP approach. In this case, the proposed RP approach is about 14.78 times faster than the MC approach.
It is north-worthy to mention that due to the smaller size of the Residue-Pole model in comparison with the Direct model, the computational cost of the Residue-Pole method is less than Direct method by a factor of 23, while the accuracy of Residue-Pole method is higher than the Direct method.