A Model-Order Reduction Approach for Electromagnetic Problems With Nonaffine Frequency Dependence

The aim of this paper is to present a novel model-order reduction (MOR) technique for the efficient frequency-domain finite-element method (FEM) simulation of microwave components. It is based on the standard reduced-basis method, but the subsequent expansion frequency points are selected following the so-called sparsified greedy strategy. This feature makes it especially useful to perform a fast-frequency sweep of problems that lead to systems of equations exhibiting a nonaffine frequency dependence. This property appears, for example, when the excitation of the problem is characterized by a frequency-dependent waveguide mode pattern, or when the computational domain includes materials with frequency-dependent permittivity or permeability tensors. Moreover, the new MOR scheme can be also used to accelerate the frequency sweep of problems with many excitations, for which the standard reduction algorithms tend to be time-consuming. Its effectiveness and accuracy is verified through analysis of three microwave structures: planar microstrip branch-line coupler, three-port waveguide junction with ferrite post, and an eighth-order dual-mode waveguide filter.


I. INTRODUCTION
A key step in the design of modern microwave devices and systems is full-wave electromagnetic simulation. The purpose of such a simulation is usually to investigate the behavior of a given structure in a specified frequency band. One of the most effective and popular techniques used for this purpose is the finite element method (FEM) [1]. The FEM is particularly widely used in the analysis of problems with high geometric and material complexity. However, FEM can be very timeconsuming in such cases, especially when the aim of the simulation is parametric analysis or optimization of a given structure over a wide frequency band.
In the last two decades, many algorithms have been proposed to accelerate frequency-domain FEM simulationsin other words, to perform a so-called fast-frequency sweep (FFS). One of the most popular FFS techniques is based on the concept of model-order reduction (MOR) [2]- [15].
The associate editor coordinating the review of this manuscript and approving it for publication was Mira Naftaly .
In this approach, a large-scale Finite-element (FE) full-order model (FOM) is projected onto a properly constructed subspace, resulting in a so-called reduced-order model (ROM). Importantly, the ROM preserves the properties of the original model with sufficient accuracy in the specified frequency band. However, since the size of the ROM is much smaller than the FOM, the frequency sweep with the ROM is much faster. In the reduced-basis method (RBM) family of MOR techniques [6], [14], [15], the projecting basis is composed of a set of field solutions (called snapshots) computed by means of FOM at selected expansion (frequency) points. The number and placement of the expansion points within the frequency range is determined by following a greedy strategy, supported by a residual-based a posteriori error estimator [14].
These techniques have been proven to be numerically stable, efficient, and reliable. However, they can only be used to accelerate certain kinds of frequency-domain simulations. More precisely, a necessary condition for their use is that the frequency parameter should appear in an explicit (''affine'') VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ form, whereas the FEM matrices (both left and right hand side components) should be frequency-independent. Unfortunately, many problems in the field of computational electromagnetics have frequency as a nonaffine parameter; for example: • materials with frequency-dependent permittivity or permeability tensors, • open-boundary problems where the computational domain is truncated using perfectly matched layers (PML), • excitations in the form of frequency-dependent mode patterns (appearing, for example, in microstrip lines).
In all these examples, the FEM matrices depend on frequency in a nonaffine way, which means that they must be regenerated for each frequency. In effect, it is not possible to efficiently estimate the error introduced by the reduced model, which is necessarily to control the whole reduction process, so the standard MOR approaches cannot be used. One possible remedy in such cases is to apply one of the general parametric MOR (PMOR) approaches typically used for problems with nonaffine parameter dependence. Publication [16] addresses this issue; it describes the PMOR method, where the computational domain is divided into so-called design and ambient regions. Then only the ambient (not the parametric) regions are subjected to the MOR process. The frequency sweep in this case is thus accelerated to a limited extent. However, it is then necessary to select the number and placement of the snapshots prior to the reduction process, so no adaptive strategy can be used. An adaptive PMOR algorithm is proposed in [17], but is restricted only to the nonaffine parameters associated with the geometry of the computational domain. Similarly, geometric nonaffine parameters can be handled by the PMOR proposed in [18], but this can be applied only to selected subregions of the computational domain. In [19], the elements of the FEM matrices with nonaffine parameter dependence are fitted to a low-order polynomial. The resulting system of equations exhibits affine parameter dependence, which is subsequently subjected to the standard moment-matching based MOR process. A similar technique has been described in [20]. However, neither of the two approaches is self-adaptive: the order of polynomials used to interpolate the system matrix elements, as well as the size of the projection basis used to generate the ROM, are set a priori. Moreover, the error of the final ROM with respect to the FOM error can be problematic to control, since it is caused both by the interpolation of nonaffine parameters and by the reduction process itself. In [21], a new generic PMOR approach is proposed for the nonaffine parameter dependence of the system, which is combined with the domain decomposition method (DDM). In order to construct the reduced-order model, the affine parametrization is extracted from the original system using Chebychev-polynomial-based interpolation. The ROMs obtained in this way preserve their passivity, reciprocity, and causality; the proof of this is based on [22] and [2]. However, neither the adaptive point-selection strategy nor the control of the accuracy of the ROMs are considered. A different strategy is proposed in [23], where the parameterdependent projection matrix is constructed by interpolating projection spaces with respect to the nonaffine parameters. This yields better accuracy and higher rates of convergence than the other approaches [19]- [21]. However, [23] can lead to significant levels of reduction error, since the order of ROM is set a priori. The technique of [23] has been extended by [24] to the case where implicit and explicit (material and geometrical) parameters are taken into account. Once the interpolation points are set and the FE system of equations is solved at each point, the PMOR projection basis is generated by means of a greedy multi-point MOR approach. It has been further extended to efficiently analyze wireless power transfer systems [25]. The approach of [26] is an improvement over that of [24] and uses the subspace interpolation method. Even though the ROMs in [24]- [26] are generated in a self-adaptive way using the greedy approach, they have a number of limitations: the nonaffine parameter dependence is approximated using the interpolation method, which leads to error; and as in the previously cited approaches, the order of the interpolating polynomial is set a priori and cannot be adaptively altered.
In order to obtain an accurate ROM in a fully automatic and self-adaptive process, additional numerical effort must be made, proposed in [27]. In particular, the following process should be considered: • The parameter space is divided into hypercubes. • In each hypercube, the projection basis and the nonaffine parameter dependencies are interpolated.
• The selected hypercubes are refined, based on the two local error indicators (associated with subspace interpolation and affine parameter reconstruction), used in the two consecutive adaptive loops. However, computing the error estimator associated with the projection basis requires the construction of an additional local ROM.
• The accuracy of the PROM is measured by means of the global error indicator. However, since the number of system matrices formed in this interpolation is larger than in the original FEM system, the efficiency of the error estimator may have deteriorated.
The complexity of this strategy makes it rather problematic to implement in a form of universal, reliable, and efficient black-box which could be used as a standard element in the microwave-engineering design process. In summary, to the best of the author's knowledge, there is no fully automatic and self-adaptive MOR algorithm for problems with nonaffine frequency dependence which is also reliable, efficient, and straightforward. This paper addresses the above issue. In the proposed technique, the projection basis used to generate a reducedorder model is composed of snapshots computed in selected frequency (expansion) points, as in the standard RBM approach. However, as observed at the beginning of this section, in this case it is not possible to efficiently assess the error introduced by the reduced model. The standard greedy approach used to select the next expansion points and the stopping criteria can thus not be used directly. Instead, a sparsified greedy strategy is applied: that is, the error is estimated only at selected frequency points, following the bisection scheme. The subsequent snapshots are computed at the frequencies where the error indicator is above the tolerance level. Such a strategy has been used previously to analyze problems with affine frequency dependence [9], [28]. In this paper it is shown that although this strategy is simpler than RBM, it can be used to perform a fast frequency sweep of complex electromagnetic problems with nonaffine frequency dependence in both the left-hand and right-hand side components of the FEM system of equations. Moreover, it can be also used to accelerate the fast frequency sweep of problems with many excitations, for which standard error estimation may be time-consuming.
The accuracy, efficiency, and reliability of this MOR scheme is illustrated by various examples, including a planar microstrip branch-line coupler, a three-port waveguide junction with a ferrite post, and an eighth-order dual-mode waveguide filter. In the last test, the efficiency of the technique is compared with that of the interpolating sweep method [29], where the transfer function of the structure is approximated by the rational function, meaning the ROMs do not need to be explicitly generated.

II. FINITE ELEMENT METHOD FORMULATION
In this section, an FEM formulation is described, focusing on the problems posed by nonaffine frequency dependence. Firstly, let us consider a source-free computational domain bounded by walls of perfect electric conductor (PEC), crosssections of n k waveguide ports, and absorbing boundary conditions of the first kind, denoted S E , S k W , and S R , respectively, for k = 1, . . . n k . The distribution of the electric field E in is determined by the following frequency-domain boundary value problem (BVP) [3], [30]: wheren is the outward unit vector, j is the imaginary unit, k 0 is the wavenumber,¯ r (k 0 ) andμ r (k 0 ) are the frequencydependent relative permittivity and permeability tensors, respectively, η 0 is the characteristic impedance of free space, and h k (k 0 ) is the frequency-dependent normalized pattern of the tangential magnetic field at the k-th port. Considering a weak form of (1) and applying FEM discretization [1] leads to an n-dimensional linear system of equations: where (s), G, and C(s) ∈ C n×n are FEM system matrices, B(s) ∈ C n×n k is the excitation matrix, E(s) ∈ C n×n k is a matrix of unknowns, u, i ∈ C n k are the vectors of amplitudes of the voltage and current waves, respectively, n k is the total number of the right-hand side vectors, and s = jk 0 is the complex frequency.
The following two subsections are devoted to a closer analysis of the matrix components of equation (2): (s), G, C(s), B(s), with a focus on the nonaffine frequency dependence.

A. EXCITATION MATRIX WITH NONAFFINE FREQUENCY DEPENDENCE
Following [3], the columns of matrix B(s) where α i are the H (curl)-conforming FE basis functions associated with the tetrahedral mesh [31].
In a case of homogeneously filled waveguides, the field pattern h k (s) used in (3) represents the TM/TE and TEM modal functions (depending on port geometry), which yields matrix B with an affine frequency dependence. However, for inhomogeneously filled waveguides (such as microstrips), the excitation field pattern cannot be expressed as TM, TE, and TEM modes, since in such cases both the transverse and longitudinal field components must be considered at the ports. In order to apply such excitation, it is necessary to construct and solve a generalized eigenvalue problem, which consists of a mixed edge and a nodal FEM formulation [32]. In effect, the propagation constant γ k and the corresponding field pattern h k (s) associated with the k-th port are computed. However, in this case, h k (s) varies with the frequency in such a way that the columns of the excitation matrix b k (s) exhibit nonaffine frequency dependence; the frequency cannot be factored out as in the affine case, namely: where b k,p is the frequency independent vector, and κ p (s) is the scalar function of s.

B. FE SYSTEM MATRICES WITH A NONAFFINE FREQUENCY DEPENDENCE
The aim of this subsection is to focus on the frequencydependence of FEM system matrices. Firstly, it is assumed that the whole computational domain contains homogeneous material, characterized in general by the frequencydependent relative permittivity and permeability tensors, denoted by¯ r (s) andμ r (s), respectively. Thus, the entries of (s) and C(s) are defined as: It is noteworthy that, in many cases (e.g., perfectly matched layers and ferrite materials), the tensors¯ r (s) andμ r (s) exhibit a nonaffine frequency dependence. This property results in nonaffine frequency dependence of the whole blocks of stiffness and mass matrices ( (s) and C(s)). Importantly, this means that these matrices have to be generated from scratch for each frequency.
In this paper, we assume that the G matrix, which represents ABC boundary conditions of the first kind, is frequencyindependent, meaning that its elements are given by: In many practical cases, the computational domain can be partitioned into n m subdomains: where it is assumed that m contains homogeneous material, characterized in general by¯ r,m (s) andμ r,m (s). However, it should be noted that computational domains can be made up of subregions with both affine and nonaffine frequency dependence. Due to assumption (7), the FE stiffness and mass matrices can be decomposed into n m blocks: where m (s) and C m (s) ∈ C n×n are associated with subregion m .

C. A FINITE ELEMENT FREQUENCY SWEEP OF PROBLEMS WITH NONAFFINE DEPENDENCE
Let us assume that the goal of a simulation is to compute the scattering parameters of a structure at the specified frequency points: s = {s 1 , s 2 , . . . s n f }. To this end, we follow a procedure which can be summarized by the pseudocode in Algorithm 1. At each frequency, we first construct the system of equations specified by (2). It is worth noting that all the nonaffine components ( (s i ), C(s i ), B(s i )) must be generated anew (line 3), taking into account the formulas given in (3) and (5). Next, system (2) is solved at each frequency (line 5), resulting in the impedance matrix Z(s i ): Finally, the scattering parameters of the structure are computed using the formula: where I is the identity matrix (line 6).

Algorithm 1 Frequency Sweep
Require: s 1 , s n f , n f , G 1: Compute matrices: B(s i ), (s i ), C(s i ) 4: It has to be emphasized that the overall frequencysweep process can be extremely time-consuming, especially when the number of unknowns n and the number of frequency points n f are large, and the FE system components (s), C(s), B(s) exhibit a nonaffine frequency dependence. A novel MOR technique to accelerate a frequency sweep in such cases is proposed in Section IV; however, for the sake of clarity, some MOR techniques for problems with affine frequency dependence will be introduced first.

III. MODEL ORDER REDUCTION FOR SYSTEMS WITH AFFINE FREQUENCY DEPENDENCE
FE systems of equations with affine parameter frequencydependence involve the matrices , G, C and B which are either frequency-independent or capable of being efficiently decomposed into the finite sums of frequency-independent matrices (as in (4)). In this section, for the sake of simplicity, only the former case is considered.
In projection-based reduction approaches, such as [3], [4], [6], [7], [12], [14], the original system of equations (2) is transformed to the following reduced-order model: where the reduced matrices are obtained using the Galerkin approach: is a properly constructed projection basis and (·) * is a conjugate transpose of a matrix. The final number of degrees of freedom of the system (11) is much smaller than that of the original system (n R n). The frequency sweep can thus be performed much faster, resulting in: (12) which aims to approximate the original transfer function Z(s) with sufficient accuracy in the frequency band.
There are several methods that aim to efficiently construct the projection basis V. This paper focuses on the reduced basis method (RBM) [4], [6], [14], which is one of the most commonly used MOR techniques in computational electromagnetics. This approach takes advantage of the fact that, in most cases, the electromagnetic field distribution does not vary significantly over the specified frequency band. Thus, the solution of (2) at any frequency can be sought in the reduced space spanned by the properly selected N solutions of (2) (called snapshots), which are computed at the expansion frequenciess 1 ,s 2 , . . . ,s N . These snapshots are then combined, forming a projection basis: which is used to generate the reduced model (11). The important questions are: which expansion points should be selected to efficiently obtain a reliable reducedmodel? And how many are there? To address these, we need to define the numerical tool called the error estimator, which enables assessment of the error introduced by the reduced model with respect to the original model. The error estimator is defined as follows [14]: where E r (s) is the solution of the following linear problem: Assuming that the matrices , G, B, C exhibit an affine frequency dependence, the blocks r , G r , B r , C r , as well as B T B, B T V, B T CV, B T BB T V and B T GV, are also frequency-independent, and can be computed just once. Since all operations in (14) and (15) are performed on low-order matrices (from C n R ×n R space), they are evaluated extremely fast, providing the error estimation in the whole specified frequency band. The error introduced by ROM is estimated each time a solution block E(s i ) is added to the projection basis V. The next expansion points i+1 is selected following the greedy strategy. It is located within the frequency band [s 1 , s n f ], at the frequency for which the estimated error value is a maximum. In other words due to this strategy the projection basis is enriched by the solution block E(s i+1 ) which is the most linearly independent with respect to the vectors included previously in the projection basis. In effect, this strategy guarantees that the level of the error introduced by the reduction process is lowered the most. Finally, the reduction algorithm stops once the error estimator falls below the accepted tolerance for the whole frequency band. Importantly, this strategy ensures full automation and reliability of the reduction process, allowing frequency sweeps to be performed on various microwave structures. However, it can only be used for problems with affine frequency dependence (or in general, with affine parameter dependence).

IV. MODEL ORDER REDUCTION FOR SYSTEMS WITH NONAFFINE FREQUENCY DEPENDENCE
Let us consider an FE system with a nonaffine frequency dependence (2). Unfortunately, such a system is unsuitable for the standard self-adaptive projection-based MOR methods [2]- [4], as it is impossible to effectively assess the error caused by the reduced model using one of the standard errorindicators (Eq. (14), [14]); this is because the following additional steps must be performed for each frequency: However, doing this at every frequency point constitutes a computation with complexity depending on n 1. Moreover, it must be performed each time the error introduced by the reduced model is estimated. This severely deteriorates the efficiency of the MOR process, making it useless from a practical point of view.
As stated in the Introduction, in order to perform an FFS of such a problem, one of the PMOR techniques can be used; these can handle nonaffine geometry parameter dependence. Such algorithms in general rely on a two-step approximation. Firstly, in order to approximate nonaffine parameter dependence, an affine decomposition of the parameter space can be carried out-for example, [17], [19], [20] on the FOM level. Such transformation leads to a system of equations with affine parameter dependence, which can then be treated directly by means of the standard MOR approach. Alternatively, ROMs can first be constructed for the selected points in the parameter space, and then the interpolation can be carried out at the reduced space level, as in [23], [24]. However, since an error in the parametric ROM (PROM) with respect to the original is introduced by these two approximations, it is problematic to control. In order to generate an accurate reduced model in a fully automatic and self-adaptive way following a two-step approximation strategy, considerable additional numerical effort is required, as described in [27], and in the Introduction to this paper. This makes it rather problematic to implement in the form of a universal black box. To address this issue, a novel MOR approach for problems with nonaffine frequency dependence is described in the next section. This approach is reliable, efficient, straightforward and, moreover, the ROM is generated in a fully automatic and self-adaptive way.

V. THE PROPOSED MOR FOR SYSTEMS WITH NONAFFINE FREQUENCY DEPENDENCE
The proposed fast frequency sweep method for problems with nonaffine frequency dependence is based on the standard RBM algorithm: that is, the projection basis is composed of the solution vectors computed at selected expansion frequency points. However, as noted in Section IV, the subsequent expansion points cannot be selected using the greedy strategy, since it is not numerically efficient in this case. To address this issue, we propose using a sparsified greedy strategy, in which the error introduced by the reduction process is assessed only at the selected frequency points, following the bisection scheme. The frequency range is divided into two subranges, and the error is estimated for each subrange at Compute matrices: B(s(i)), (s(i)), C(s(i)) 6: A(s(i)) = (s(i)) +s(i)G +s(i) 2 C(s(i)) 7: V i =s(i)A(s(i)) −1 B(s(i)) //solve system of equations 8: V = GS(V, V i ) // Gram-Schmidt process 9: end for 10: i = 3 11: while TRUE do 12: len_s = length(s) 13: if i > len_s then 14: BREAK 15: end if 16: while i ≤ len_s do 17: 18: for j = 1:2 do 19: Compute matrices: B(s e (j)), (s e (j)), C(s e (j)) 20: Compute reduced matrices: r (s e (j)), B r (s e (j)), C r (s e (j)), G r , B(s e (j)) T B(s e (j)), B(s e (j)) T (s e (j))V, B(s e (j)) T C(s e (j))V, B(s e (j)) T GV, B(s e (j)) T B(s e (j))B(s e (j)) T V 21: Solve eq. (15) for E r (s e (j)) 22: Compute estimated error: e s (s e (j)) using eq. (14) 23: if e s (s e (j)) ≥ tol then 24 its middle point, where the accuracy of the ROM is assumed to be the lowest. If the error indication therein is above the tolerance level, the FOM system of equations is solved, subsequent snapshot is added to the projection basis, and the subrange is halved. The algorithm stops once the estimated error in all the points falls below the tolerance level.
The proposed MOR approach is summarized in the pseudocode in Algorithm 2. This starts by specifying the input parameters, which are the lower and upper frequency limits s 1 and s n f , the error tolerance (tol), and the matrix G, which is assumed to be frequency-independent. Next, the frequency Project onto a reduced space: In the second part of the approach, the subsequent expansion points are selected by following the bisection approach: the computational domain is divided into two subintervals and the error is estimated at the middle points, where the error introduced by ROM is expected to be the greatest (s M − h and s M + h, line 17). Matrices are generated to estimate the error (lines 19-20), the reduced system of equations is solved, giving E r (s e (j)) (line 21) and finally, e s (s e (j)) is computed for a single frequency point s e (j) in line 22. If the error indicator is above tol, the subsequent expansion point (s e (j)) is added to the sets (line 24). The global FEM system of equations is constructed (line 25) and solved (line 26) at this frequency. The computed block V i is orthogonalized using the Gram-Schmidt method, and added to the projection basis V. The algorithm stops once all the expansion points in the sets have been examined (lines [13][14][15]. It is worth noting that the steps 5 and 19 can be performed extremely rapidly using the fast FEM matrix generation approach [33]. The algorithm returns the projection basis V, which is subsequently used to generate a reduced model and perform a fast frequency sweep. The fast frequency sweep process is summarized by the pseudocode in Algorithm 3. For each of n f frequency points, the nonaffine matrices (B(s i ), (s i ), C(s i )) are computed and projected onto a reduced subspace V (lines 4 and 5, respectively). Note that the G matrix is frequency-independent, so it FIGURE 1. Scattering parameters of the planar microstrip branch-line coupler [34], computed using the FEM formulation (black color), the proposed method (red color); and as measured (dots).
is reduced only once. Next, the left-hand side of the reduced system of equations (A r (s i )) is generated, and the system is solved, in order to compute the reduced transfer (impedance) function Z r (s i ) (line 7), and then the scattering parameters S r (s i ) (line 8).
Remark: It should be noted that the fast frequency sweep loop should contain only operations with complexity depending on n R , as in standard MOR approaches. However, the FFS described by Algorithm 3 contains operations with complexity depending on n (where n n R ). It may thus seem that our frequency sweep is more time-consuming; however, in most cases, the nonaffine frequency dependence is associated only with parts of the computational domain, such as ferrite pucks and PML absorbing conditions. Thus, in order to generate r (s i ) and C r (s i ), only small parts of the FEM system matrices (s i ) and C(s i ) need be updated and projected onto a reduced subspace V. By the same token, the nonzero elements of B(s i ) correspond to excitations, which is typically applied only to small parts of the computational domain, and which can be updated on the fly.
In other words, at each individual frequency point, instead of constructing and solving the original FEM system of equations (Algorithm 1, lines 4 and 5), the matrices are instead projected onto the reduced subspace, and the reduced system of equations is constructed and solved (Algorithm 3, lines 5, 6, and 7, respectively). The computational time of the second set of operations is usually much lower as compared to that of the first set.

VI. NUMERICAL TESTS
We analyzed three structures in order to validate the efficiency and accuracy of our MOR approach. All computations were performed on an Intel Core i5-6500 processor with 64 GB RAM.

A. PLANAR MICROSTRIP BRANCH-LINE COUPLER
The first numerical example focuses on the planar microstrip branch-line coupler, fed by microstrip ports, as described in detail in [34] and shown in Fig. 1. The goal of the first simulation was to compute the scattering parameters of the structure at 101 equidistantly distributed points in the 1.6-3.4 GHz band, using the standard FEM formulation FIGURE 2. Model-order reduction scheme, planar microstrip branch-line coupler. Black squares, black dots and black circles denote the initial expansion points, selected expansion points, and candidate points, respectively.

FIGURE 3.
The actual error associated with all elements of scattering matrix, computed using the proposed MOR approach for the planar microstrip branch-line coupler [34] with the red line being the tol level. [35], with the computational domain truncated using absorbing boundary conditions. The FEM system of equations has 328,868 degrees of freedom, and the whole frequency-sweep process took 1127.2 s. The resulting scattering characteristics are shown in Figure 1; this shows that they are well correlated with the measurements taken from [34].
Since the structure is fed by the microstrip ports, the electric field pattern at the port surfaces varies significantly with frequency and, in effect, the RHS vectors are of a nonaffine type. The coupler can thus not be easily analyzed using the standard MOR methods. In fact, one can first construct an auxiliary reduced-order model which approximates the field distribution at the ports, and then perform the standard reduction [36]. However, this approach is not straightforward, and the error from the auxiliary ROM can affect the global error. Alternatively, a fast-frequency sweep of this problem can be carried out with ease using our approach. To this end, a projection basis and the corresponding reduced-order model were constructed on the basis of the algorithm described in Section V. The frequency points at which the snapshots are computed are shown in the reduction scheme (Fig. 2). In the first reduction level, the projection basis contained three snapshots at 1.6, 2.5, and 3.4 GHz. In the second level, the error was assessed at two frequency points (2.05 and 2.95 GHz); in both cases, the error was above the tolerance level (tol = 1e − 6), so the next two snapshots were added to the projection basis. In the third reduction level, the error was assessed at four frequency points (denoted by the circles), and the error level was found to be below tol in all cases, ending the reduction process. In summary, in order to construct the projection basis, the full FEM system of equations needed to be solved at the five frequency points denoted by black dots and squares (Fig. 2). Since each snapshot contained four vectors, corresponding to each excitation port, the final size of the projection basis (as well as the ROM) was 20. The computational time needed for the whole MOR process, including the fast-frequency sweep, was 70.9 s, making it nearly sixteen times faster than the standard frequency sweep. Figure 1 shows that the results computed using the fullorder FEM model and the reduced model are indistinguishable. However, in order to compare the results more precisely, the following definition of the actual error is used: where S ij (k 0 ) MOR and S ij (k 0 ) REF denote the scattering parameters computed using the reduced and the original model, respectively, with i and j being the output and input ports of the structure. A plot of E ACT ij (k 0 ) is shown in Figure 3, and it can be seen that it is below the tol level (denoted by the red line) for the whole frequency band.

B. THREE-PORT WAVEGUIDE JUNCTION WITH A FERRITE POST
The second numerical test deals with a circular cavity (height: h = 10.16 mm, radius: r = 14.7 mm) fed by three WR-90 waveguide ports (22.86 × 10.16 mm), as shown in Figure 4. Inside the cavity, a cylindrical ferrite post H = 4.65 mm high and with a radius of R = 3.75 mm is located centrally on the support (with height H s = 5.06 mm and radius R s = 9 mm). The ferrite is magnetized along the z-axis and its magnetic permeability is described by a frequency-dependent tensor of the form: where the formulas for the susceptibilities χ xx , χ yy , χ yx and χ xy are provided in [37]. The parameters of the ferrite material used in the example are as follows: r = 14.4; magnetic saturation: M S =133000 A/m; linewidth of the gyromagnetic resonance corresponding to ferrite losses: H = 266 A/m; and static magnetic field strength: H 0 = 66500 A/m. The circulator was first analyzed using the standard FEM procedure in order to obtain the reference results. The computational domain was discretized using 43,051 tetrahedrons, which resulted in the FEM system of equations with 263,908 unknowns. Next, the scattering parameters S 11 , S 21 FIGURE 5. Model-order reduction scheme for the three-port waveguide junction case. Black squares, black dots, and black circles denote the initial expansion points, selected expansion points, and candidate points, respectively. FIGURE 6. The actual error for the three-port waveguide junction. and S 31 were computed at 101 equidistantly distributed frequency points in the 7-13 GHz frequency band. The computations took 1266.7 s and the results are presented in Figure 4.
Since the ferromagnetic cylinder is characterized by a nonaffine frequency-dependent tensor (17), the resultant FEM matrix (s) is also of a nonaffine kind (see (5) for details). Therefore, as in Section VI-A, the scattering parameters of this problem cannot be computed using standard MOR methods. However, one can perform a fast-frequency sweep of this problem using our proposed MOR technique. Figure 5 shows the reduction scheme. In the first reduction level, the projection base was composed of three snapshots computed at 7, 10, and 13 GHz. In the second and third levels, the error was assessed at two and four frequency points, respectively, and in all cases the error was above the tolerance (tol = 1e − 3). In the fourth reduction level, the error was assessed at eight frequency points; however only in the case of f = 11.875 GHz was it above the tol. In the fifth level, the estimated error was below tol at the two considered points, which halted the reduction process. Taking into account the fact that each snapshot contained three vectors, the final size of the ROM is just 30, requiring 10 system matrix factorizations. The overall MOR process took 185.5 s, which is nearly seven times faster than the standard frequency sweep. Figure 6 shows that the actual error is below the tol level over the whole frequency band. It should be noted that the FEM formulation leads in this case to a nonsymmetric fullorder model. Higher accuracy in the reduced model can thus be obtained using a two-sided projection instead of a onesided projection (see [38] for details).

C. EIGHTH-ORDER DUAL-MODE WAVEGUIDE FILTER
In the last example, we consider the eighth-order dualmode waveguide filter [4] whose geometry is shown in Figure 7. We were interested in a wide-band filter response, so 11-19 GHz was specified as the band of interest, with 401 frequency points. FE discretization led to FOM with n = 202474. In order to accurately model the structure, all the modes excited in the band of interest have been considered, which resulted in eight excitation vectors for each port. The direct sweep using this model took as long as 339 s. The resultant S 11 and S 21 are shown in Figure 7.
Since the filter is non-lossy and the waveguide excitation (TE and TM waveguide modes) are of an affine kind, this can be analyzed using the standard MOR schemes. However, the excitation matrix contains sixteen vectors, so the error estimator defined in (14) significantly deteriorates the efficiency of the reduction scheme. Assuming tol = 1e − 6, the computational time needed for the RGM-MOR (reliable greedy multipoint model-order reduction) [5] and RBM algorithms was 107.5 s and 281.7 s, respectively.
The reduction scheme described in this paper was applied next. It can be seen in Figure 8 that the reduction process requires eight levels, with three, two, four, eight, two, one, one, and zero expansion points in each of them, respectively. Importantly, the expansion points are distributed mainly in the pass-band of the filter, where the field solutions that contribute the most to the behavior of the filter are located (see [39] for details). Figure 7 shows that the characteristics computed using FOM and ROM are indistinguishable, and the actual error is below the tol level over the whole frequency band (Fig. 9). The total MOR computation time was only 58.96 s, which is nearly six times faster than the direct sweep, nearly five times faster than RBM, and twice as fast as RGM-MOR.

D. INTERPOLATING SWEEP
Alternatively, in order to rapidly compute the scattering parameters of these structures, the so-called interpolatingsweep approach can be used. In this, unlike in our proposed method, the reduced-order model is not generated explicitly. Instead, the scattering parameters are computed at the selected frequency-points and the results are used to construct FIGURE 8. Model-order reduction scheme for an eighth-order dual-mode waveguide filter case. Black squares, black dots, and black circles denote the initial expansion points, selected expansion points, and candidate points, respectively. the rational function that approximates the transfer function of the structure. In the last test we used the interpolating sweep, based on the vector fitting method [29], to compute the scattering parameters of the three-port waveguide junction (considered in section VI-B). In order to obtain the rational function, which approximates the original model with the true error below 1e-3, as many as 23 FEM system matrix factorizations were needed; this is much more than in our proposed method, where only ten factorizations were needed. In effect, the computational time is significantly longer: 297.4 s, compared with 185.5 s for the proposed method.

VII. CONCLUSION
This paper has presented a novel self-adaptive model-order reduction approach to perform a fast frequency sweep of complex electromagnetic problems with nonaffine frequency dependence in both the left and right-hand side components of the FEM system of equations.
In order to select the subsequent expansion frequency points to construct a projection basis, a sparsified greedy strategy is used: the error is estimated only at the selected frequency points, following the bisection scheme. Our numerical tests of a planar microstrip branch-line coupler, a three-port waveguide junction with a ferrite post, and an eighth-order dual-mode waveguide filter have demonstrated the reliability, accuracy, and excellent computational performance of this technique.