Regularized and Compressed Large-Scale Rational Macromodeling: Theory and Application to Energy-Selective Shielding Enclosures

In this article, we introduce a robust procedure for the extraction of passive rational macromodels of low-loss electromagnetic structures with massive port counts. Such structures pose inherent challenges that make standard macromodeling tools and approaches inadequate, mainly due to complexity and sensitivity at low frequency. The proposed approach involves a preprocessing stage in which port response data from a full-wave electromagnetic solver are regularized and extrapolated to dc using an asymptotic modal representation. The resulting data samples are then processed by a dedicated compression algorithm that represents the full set of input–output responses in terms of a few basis functions, which are constructed by enforcing an exact low-frequency modal asymptotic behavior, possibly including higher order dc zeros. These zeros are preserved in any stage of rational fitting and passivity enforcement, resulting in dc and low-frequency compliant compressed passive macromodels. Numerical results with up to 400 ports demonstrate the superior performance and accuracy of the computed models with respect to state-of-the-art approaches. In particular, the resulting models preserve their accuracy irrespective of the loading conditions, including the limit cases of short and open terminations.

When low-loss electromagnetic structures with a large number of interface ports are considered, standard passive macromodeling algorithms may fail. Although advanced formulations of rational fitting and passivity enforcement algorithms exist [9]- [11], including parallelization for multicore [12] and dedicated (GPU) hardware [13], several open issues remain in the generation of robust macromodels. One main aspect is model sensitivity to loading conditions.
A black-box macromodel is usually very accurate in the input-output representation used for its training. For instance, fitting scattering responses with port reference R 0 usually result in a model that is very accurate when loaded with resistances close to R 0 . However, when the terminations differ significantly, the model accuracy may deteriorate [14]. Consider as an example the limit cases of short-and open-circuit terminations, for which model behavior can be checked by computing its admittance or impedance responses. Due to the nonlinear transformation that converts scattering to impedance or admittance, the inevitable approximation errors affecting the model may be amplified, resulting in dramatically wrong and nonphysical behavior [15], [16]. This phenomenon is particularly evident when the system is low-loss, in which case the scattering responses (eigenvalues) approach magnitude one at low frequency. Increasing the number of ports further exacerbates this sensitivity, in addition to posing additional challenges due to computational requirements.
This article proposes a sequence of steps aimed at reducing and possibly eliminating the above sensitivity problems, including the case of missing samples at low frequency. This situation is common due to low-frequency limitations of full-wave solvers, yet the low-frequency band is typically the range where sensitivity induces its worse effects.
Both theory and numerical results are here illustrated on a particular class of structures intended for energy-selective shielding [17]. These structures are chosen not only for the specific application interest but also because their electromagnetic behavior includes all ingredients that induce a worst-case sensitivity scenario: the combined presence of large-scale lowloss multiport electromagnetic systems with massive nonlinear loading. In view of this application, motivation and objectives of this work are now outlined with more details.

A. Motivation and Objectives
The left panel of Fig. 1 depicts a cubic box-shaped shielding enclosure (each side 500 mm) with a square aperture (side of 250 mm). A regular grid of P = p × p lumped ports covers the aperture, with p parallel segments, each including p ports connected in series by metal strips (width 2 mm). The dc conducting paths are illustrated in the right panel of Fig. 1, together with the port numbering order. The shield material is assumed to be a perfect electric conductor. This test case will be used as a running example to illustrate all steps of the proposed approach. We remark that the proposed approach can be trivially generalized to apertures with arbitrary shapes and nonsquare grids.
It has been shown [18]- [21] that loading the ports by nonlinear elements, e.g., back-to-back diode pairs, turns the structure into an energy-selective shield. Penetration of an incident wave into the box is permitted when the field amplitude is below the critical threshold that triggers diode conduction. A higher energy field switches diodes into conduction mode, thereby shorting the conductive paths of the strips and increasing shielding effectiveness. The result is an energy-selective shield, which can be used to protect any enclosed equipment from high-amplitude incident fields while allowing intended low-power communication.
The shielding enclosure with its P unloaded ports is a linear electromagnetic structure, which can be characterized through its frequency-domain responses within the desired frequency band through a field solver. These responses can be used as the training dataset for the extraction of a rational macromodel, which, in turn, can be synthesized as an equivalent circuit, terminated by the required nonlinear loads, and simulated in time domain to verify the overall shielding performance. See [18] for a complete description of this characterization and modeling framework.
The top panel of Fig. 2 reports the impedance responses obtained by a Method of Moments (MoM) solver [22]. Native samples from field solvers {H(jω k ), k = 1, . . . , K} are typically available in scattering representation (H = S), conversion to impedance Z or admittance Y is performed as postprocessing. The data cover a frequency band ω k ∈ Ω D = Ω L ∪ Ω H , Fig. 2. Full set of responses of a 25-port shielding enclosure. Top: all impedance matrix data from field solver. Bottom: field solver data after conversion to the asymptotic modal domain via (6). See text for details.
where Ω L is the frequency band over which the structure can be assumed to be electrically small, where the corresponding impedance responses exhibit a clear asymptotic behavior, and Ω H is the higher frequency band where box resonances start to appear. Due to limitations in the adopted field solver, only responses starting from a minimum frequency f min are available, thus causing a gap in the frequency band Ω G = (0, f min ). Note that such low-frequency limitations affect most full-wave field solvers, and that the presence of a low-frequency gap Ω G is very common.
The responses in Fig. 2 present the following challenges.
1) The dc (zero-frequency) point is missing. The behavior of the structure under static conditions is, however, essential for the correct setup of any circuit simulation in time domain. Since the field solver is unable to provide the dc responses accurately, a common practice is to evaluate the dc point via a physics-based approach [23] or via a numerical solution through a static solver. In fact, the right panel of Fig. 1 provides the dc circuit, whose direct analysis shows that both impedance and admittance matrices do not exist. This is due to the presence of a dominant pole at dc, which is well visible from the asymptotic behavior of the impedance samples in Fig. 2. The same behavior can be observed in the admittance responses (not shown). Allowing this dc pole in a macromodel would result in time-domain simulation difficulties due to the lack of asymptotic stability. Therefore, some sort of regularization at dc is necessary in order to guarantee success in subsequent transient circuit simulations. 2) Even if a physically consistent and regularized dc point is provided, the presence of a wide low-frequency gap Ω G prevents full control of the macromodel response in this frequency range, which is, however, essential to control Fig. 3. Modeling a 25-port shielding enclosure with a standard passive macromodeling tool [24]. Scattering parameters were fitted by enforcing an exact dc point (top panel). Low-frequency accuracy is not preserved after model conversion to the impedance parameters.
the transient response. Even in the case of high-frequency incident field excitation, the nonlinear terminations may enhance the low-frequency portion of the port signal spectra, which are, in turn, affected by the low-frequency model response. Fig. 3 reports the impedance model responses obtained by a state-of-the-art passive macromodeling tool [24] by fitting scattering responses with exact enforcement of a physics-based dc point. We see that the impedance responses of the model in a broad low-frequency band are completely wrong. We conclude that some sort of data extrapolation is required to provide a smooth connection between the field solver data and the dc point so that the model can be properly trained. We will see that a direct naive extrapolation of impedance (or admittance or scattering) responses by extending the asymptotic behavior to low frequencies is not sufficient, and a more sophisticated regularized modal extrapolation is required.
3) The third macromodeling challenge relates to complexity, especially in terms of the number of ports P . For practical applications, this number can grow to several hundreds or more due to the requirement of uniformly covering the shielding aperture with nonlinear blocking devices. An approach that has the potential to scale favorably with the number of ports is, therefore, required. These challenges will be addressed as follows. First, a combined low-frequency regularization and extrapolation based on asymptotic modes [25] is discussed in Section II, by extending the preliminary results of [26]. We will see that this preprocessing stage smoothly extends dc to the solver data while at the same time avoiding the low-frequency degeneracy observed in Fig. 2. Data and model complexities are addressed through a compression algorithm based on the singular value decomposition (SVD). We extend the approach of [27] by proposing a block-diagonal SVD (BD-SVD) and a hierarchical SVD (Hi-SVD). These alternative solutions compress the full set of input-output responses into a small number of basis functions, each with a well-defined low-frequency asymptotic behavior and the number of expected dc zeros. The generation of regularized and compressed macromodels is discussed in Section IV, where an approach is presented for rational fitting the individual basis functions while preserving the number of dc zeros in the model. A full set of numerical results presented in Section V demonstrates the reduced sensitivity of the proposed models by increasing the number of ports up to P = 400. Conclusions are finally drawn in Section VI.
Before proceeding, we set some basic notation used throughout this document. We denote scalars with normal font x, vectors with lower case bold fonts x, and matrices with upper case bold fonts X, with I n being the identity matrix of size n. The transpose and Hermitian transpose of a matrix will be indicated with X T and X H , respectively. The singular values of a matrix X will be indicated as σ{X}, while Re{·} and Im{·} extract the real and imaginary parts of their arguments, respectively.

II. REGULARIZATION AND LOW-FREQUENCY EXTRAPOLATION
Since a low-frequency gap Ω G prevents control of model accuracy and sensitivity, we propose in this section a simple and robust procedure to fill this gap with physically consistent extrapolated samples such that the frequency responses are well defined and nonsingular at any frequency (including dc) for any input-output representation.

A. DC Regularization
Let us consider the dc circuit of a box structure with no metal losses, as depicted in Fig. 1. As already mentioned, both admittance and impedance matrices are ill-defined at DC. 1) Exciting all ports with independent current sources attempting to evaluate impedance parameters leads to an ill-defined dc circuit due to the presence of p(p − 1) independent current source cutsets. 2) Exciting all ports with independent voltage sources attempting to evaluate admittance parameters leads to an ill-defined dc circuit due to the presence of p independent voltage source loops. Therefore, both admittance Y 0 = Y(s = 0) and impedance parameters Z 0 = Z(s = 0) are ill-defined at dc since both Y(s) and Z(s) have a pole at s = 0. This pole can be moved off the imaginary axis by: 1) adding a shunt resistance R to each port, leading to regularization of Z(s); 2) adding a series resistance r to each port, leading to regularization of Y(s). This procedure leads to the dc regularized topology of Fig. 4, which is nonsingular for any cardinal representation. In fact, a direct calculation leads to the following expressions for the dc impedance, admittance, and scattering matrices: where p is the number of vertically aligned ports, ϑ = 1/p, u T = [1, . . . , 1], u ∈ R p , Φ = R/R 0 , ϕ = r/R 0 , Φ = Φ + ϕ, with R 0 denoting the scattering port reference impedance.

B. Regularization of MoM Data
Redefinition of the dc behavior according to Fig. 4 makes the frequency samples obtained from the field solver not compatible with the extrapolated dc responses. Therefore, a full-bandwidth regularization process must be applied so that all frequency samples in Ω D are compatible with the above-computed dc values. This regularization is accomplished by Algorithm 1, which evaluates a set of modified responses obtained by connecting a shunt resistance R at each port and, in a second step, a series resistance r at each port. Note that all steps are well defined since the dc point is assumed to be missing from the data so that all matrices are invertible at any available (finite) frequency point. The results of this process are depicted in Fig. 5. Fig. 5 shows that the perturbation of the available frequency samples from the MoM solver is negligible, considering the values used for series (r = 0.1 Ω) and parallel (R = 100 MΩ) regularization. As a general guideline, the values of r and R should be very small and very large, respectively, in order to preserve all features in the original responses. Fig. 5 also shows Algorithm 1: Data Regularization.
Require: (1) as Z(jω)| ω=0 = Z 0 7: return Z(jω k ): regularized data including dc point that the regularized dc values do not connect smoothly to the frequency responses starting from the first frequency sample at f min . The gap in the responses due to the low-frequency limit of the field solver is still too large and must be filled with suitably extrapolated response samples.

C. Low-Frequency Extrapolation
An extrapolation of the frequency responses below the first available frequency point is not trivial since both admittance and impedance responses Y(jω) → ∞ and Z(jω) → ∞ for ω → 0. The structure is characterized by the presence of both inductive and capacitive modes at low frequency, which are responsible for the singularity of both admittance and impedance. A direct extrapolation of admittance (impedance) matrix elements would lead to a good approximation of the inductive (capacitive) modes, respectively, leaving the other modes very inaccurate.
This problem can be effectively solved by conversion to a modal domain, where inductive and capacitive modes are well separated and can be extrapolated concurrently [26]. This operation is best performed on the raw data (lossless) from the field solver before applying Algorithm 1. Differently from [28], in order to improve numerical robustness and provide an extrapolation scheme that is immune from the field solver inaccuracies that can be expected at low frequencies, we base the definition of the modal basis on the regularized dc admittance Y 0 derived analytically in (2). This choice leads to the additional benefit of a purely real and orthogonal modal transformation, which is required for preserving causality [29], [30] in all subsequent macromodeling steps.
We compute the following eigendecomposition: where R = R + r. We note that the eigenvalue 1/R has multiplicity c = p(p − 1), while the eigenvalue 1/r has multiplicity = p. We will see shortly that c and correspond to the number of capacitive and inductive modes of the structure. Since Y 0 is symmetric, the modal matrix Q = Q p Q s is orthogonal so For any frequency ω k ∈ Ω D , the data from the field solver are converted to the dc asymptotic modal domain via The same transformation can be applied in scattering or impedance representation. The bottom panel of Fig. 2 depicts the obtained modal impedance responses, where the ±20 dB/dec slope associated with purely inductive/capacitive modes is evident for frequencies in the subband Ω L . Note that the modes that converge to zero in the impedance representation are those modes that explode to ∞ in the admittance representation and vice versa. With reference to Fig. 2, we now restrict our analysis to the subband Ω L where the projected solver data can be assumed to behave asymptotically. We can assume the following frequency dependence: where Γ m and C m are constant matrices such that (Γ m ) ij = 0 ⇒ (C m ) ij = 0 for any i, j and vice versa. The structure of the dominant terms is also detailed in (7) in terms of the constant matrix blocks C, Γ, and the off-diagonal blocks X. These constant matrices are determined via a two-step elementwise regression. First, we assign each matrix element (i, j) to the capacitive or inductive subset by determining the dominant slope of the modal admittance magnitude. For ω k ∈ Ω L , we compute ξ k = log 10 ω k , y ij;k = log 10 |Y m;ij (jω k )| and we perform a least squares fit to determine the coefficient μ ij such that Acceptable values for this coefficient are μ ij = +1, in which case element (i, j) is assigned to the capacitive terms C m , and μ ij = −1, corresponding to the inductive terms Γ m . Coefficient ν ij is disregarded. As a second step, we identify the actual elements of constant matrices C m and Γ m . This is achieved by a second least squares fit based on the imaginary part η ij;k = Im{Y m;ij (jω k )} performed differently for capacitive and inductive terms as Using the coefficients from (9), we generate the desired number K GL of new samples according to the asymptotic modal expansion Note that in addition to new samples in the low-frequency gap Ω G , we also evaluate the asymptotic model in its training band Ω L , expecting a minimum deviation with respect to the projected solver data. These are used to assemble all native and extrapolated samples as where α k and β k are frequency-dependent coefficients such that α k + β k = 1 that are respectively (linearly or logarithmically) increasing and decreasing from 0 to 1 in Ω L , providing a smooth transition between the synthetically-generated asymptotic samples Υ m (jω k ) and the field solver data Y m (jω k ).

D. Regularization in the Modal Domain
The last step in the proposed data conditioning process is regularization, which is performed directly in the modal domain by processing the extrapolated data (11). The main advantage of this approach is that only the individual modes that are responsible for the singularity of a given (modal) matrix representation are regularized, leaving the other modes unperturbed.
Starting from Y m (jω k ) defined in (11), for any frequency ω k = 0, we perform the following steps.
1) Regularize the capacitive block to prevent the degeneracy The resulting matrix is invertible at any frequency.
2) Convert to impedance and regularize the inductive block to prevent its degeneracy for ω → 0 3) Return to the physical domain through inverse modal transformation Fig . 6 depicts the impedance responses after the complete extrapolation and regularization, compared to the native field solver samples. This plot confirms that in-band samples are practically left unchanged, and that data are smoothly extrapolated throughout both the transition band Ω L and the low-frequency gap Ω G .
It can be easily shown that the resulting data samples in the low-frequency gap Ω G , obtained as a regularization of the lossless extrapolation (10), have the modal structurȇ where the symbol ν denotes a frequency-dependent matrix with leading order (jω) ν for ω → 0. This structure is fully compatible with the direct (physical-domain) regularization of Algorithm 1 and with the dc physics-based circuit shown in Fig. 4. The modal impedance matrix, thus, exploits elements that converge to a constant (the diagonal entries), elements with a first-order zero, and off-diagonal blocks with a second-order zero. This structure also applies to the scattering responses in Ω G , as confirmed later by Fig. 9(d). Classification of these matrix elements is important for data compression, as discussed next.

E. Extensions
The extension of the proposed procedure to a general multiport structures with possibly different asymptotic behavior at dc is straightforward once a proper regularization circuit is defined (e.g., Fig. 4). Indeed, the general procedure can be easily modified to deal with purely inductive or capacitive modes, corresponding to admittance or impedance matrices being illdefined at dc.
We also remark that the proposed regularization procedure was conceived to handle the case of a completely missing dc characterization, with the objective of providing a synthetic dc point that will not create problems in later modeling and simulation steps. In other situations, a dc characterization may be available, either from direct measurement or from a separate dc or quasi-static solver. If available, such dc characterization can be used directly within the proposed framework by replacing the expressions (1)-(3). The eigendecomposition and modal regularization then apply without significant modifications.

III. DATA REDUCTION AND STRUCTURED MACROMODELS
In this section, we address the scalability of macromodel extraction with the number of interface ports. For the considered energy-selective shielding application, we expect grids with hundreds of ports P . For such a large number of ports, the complexity of both model generation and model exploitation in successive transient simulations may become impractical due to the requirements of concurrently fitting P 2 responses. Although the Fast VF algorithm [9] and its parallel implementations [12], [13] somewhat improve the model fitting phase by compressing the least squares system returning the model coefficients at each VF iteration, there is still a significant margin for improvement. In fact, the set of frequency responses usually exhibits a high degree of redundancy, which can be removed by a suitable data compression strategy. The structure of the model can, thus, be improved to take advantage of such redundancy rather than adopting a standard rational expansion of each individual response with common poles and independent P × P residue matrices.
In this direction, a data compression technique based on a standard SVD was originally presented in [27]. This method allows to obtain a compressed-macromodel via VF taking advantage of the spatial correlation between ports and describing the overall system with a limited number of basis functions. In [27], these functions are the singular vectors of a truncated SVD. Here, we generalize and customize this approach to the particular modal structure of the frequency responses of shielding enclosures, in particular, to reproduce accurately the number of modal dc zeros. We first recall some background notation and material in Section III-A. Our proposed extended and generalized version of this algorithm is presented in Section III-B.

A. SVD Data Compression: Background
The data compression approach of [27] starts with a set of frequency (typically scattering) responses H(jω k ) = H k of a P -port linear time-invariant system, sampled at suitable frequency points ω k , with k = 1, . . . , K. The first step of the procedure stacks the columns of H k ∈ C P ×P at each frequency ω k in a row vector defined as where x k ∈ C P 2 . The equivalent mapping (x k ) = (H k ) i,j can be defined as where · rounds toward infinity and mod is the remainder of integer division. All these row vectors are collected in a matrix X ∈ C P 2 as Following [27], real and imaginary parts of X are stacked in a real matrix, which is subjected to a truncated SVD where the ρ leading singular values are collected inΣ ∈ R ρ×ρ , while the left and right singular vectors correspond to the columns ofŪ ∈ R 2K×ρ andV ∈ R P 2 ×ρ , respectively. We recall thatŪ TŪ = I ρ andV TV = I ρ , whileVV T = I P 2 due to SVD truncation. Next, we define the following matrix: which recombines real and imaginary parts of scaled singular vectors, allowing to approximate the system responses in (18) as In particular, the th column of X can be reconstructed using wherew q denotes the qth column ofW. The set ofw q is here denoted as basis functions. From (21) or (22), we see that each frequency response out of a set of P 2 can be approximated with only ρ basis functions through a simple SVD data preprocessing. Therefore, a VF macromodel of the P × P system can be built by constructing rational approximations of the basis functions {w q , q = 1, . . . , ρ}, which is a very efficient operation when ρ P 2 . Moreover, as indicated in [27], the error introduced by (21) is bounded by where σ ρ+1 is the first neglected singular value. Note that (23) sets the limit of accuracy of the overall macromodeling procedure: any error introduced at this step will not be recovered in later steps so that the number of retained singular values should be such that σ ρ+1 is sufficiently small for the considered application.
If we apply the above framework to scattering responses of shielding enclosures, the structure of the responses is not preserved by the truncated approximation (21). Fig. 7 confirms that even if the threshold used for SVD truncation is very aggressive, the presence of single or double zeros at dc for some responses is not preserved. The data before compression for this example were obtained following the regularization and extrapolation discussed in Section II-C, using series and parallel resistances r = 0.1Ω and R = 100MΩ. The resulting number of frequency samples is K = 8904 in the range [0, 1] GHz, including 904 low-frequency extrapolated samples and the dc point. The SVD data compression was applied using a truncation threshold σ < 10 −6 , resulting in ρ = 73 basis functions (see bottom panel of Fig. 7). All these basis functions are nonzero at dc so that the presence of a dc zero in any of the system responses would be enforced through a linear combination (22) of nonvanishing basis functions, relying on numerical cancelation. Machine precision and SVD truncation thresholds prevent this cancelation so that a low-frequency saturation appears and destroys any dc zero (see top panel of Fig. 7). This saturation can be interpreted as an artificial additional loss at low frequency, which may compromise the overall reliability and accuracy of the model, especially in transient simulations with nonlinear terminations. The presence of dc zeros must be preserved at all steps of the modeling process, including data compression. Therefore, we need to extend the basic SVD compression of [27] in order to guarantee a proper structure preservation.

B. Structured SVD Compression
The main idea that motivates the following data reduction procedure is to obtain a specific set of basis functions that reflect structural properties in the original response data. This is applied here to preserve dc zeros of a given multiplicity (up to two in the following), although the procedure is general and can be applied to preserve any other feature of interest.
The starting point is (18), which collects all P × P system responses as columns of matrix X, ordered according to the column stacking operator (16). As a first step, we reorder the columns of X through a suitable permutation matrix P such that PP T = I as where the first block M 1 ∈ C K×P 1 stores the responses of H(s) that saturate to a nonvanishing dc value, while the second and third blocks M 2 and M 3 collect all responses with a zero at dc of order 1 or 2, respectively. In the following, we indicate with P ν the number of columns of each group of responses M ν , with ν = {1, 2, 3} and ν P ν = P 2 . An example of these three sets is reported in the top panel of Fig. 8 [see also Fig. 9(d)]. Applying a standard (unstructured) SVD compression to the reordered matrix M and highlighting the 3 × 3 block structure leads to where matrices W ν collect groups of basis functions to be associated with the three groups of responses M ν , with ν = {1, 2, 3}. Expanding (25) provides   (26) is guaranteed if V 12 = 0. In other words, enforcing V T to be block-lower-triangular guarantees the preservation of dc zeros in the data compression, provided that individual basis functions in the three groups W ν are constructed using only the proper subsets of responses. We can realize this condition using two different strategies.

1) Block-Diagonal SVD:
The block-diagonal approach considers each subset of responses M ν for ν = {1, 2, 3} as independent, by applying individual SVD compression. Therefore, we compute whereΣ ν collects the leading ρ ν singular values in descending order, and whereŪ ν ∈ R 2K×ρ ν ,V ν ∈ R P ν ×ρ ν , withŪ T νŪν = I andV T νVν = I. We then definē and rewrite (27) as Collecting all groups leads to the BD-SVD approximation The advantages of the BD-SVD compression are as follows: 1) the guarantee that particular features in different groups of responses are preserved in the compressed dataset, including dc zeros; 2) the ability to tune the accuracy on each set of responses M ν individually, which allows to control the dc error when working in the modal reference systems. Indeed, in the latter case, the diagonal terms of the transfer function are the only ones responsible for the accuracy at s = 0; 3) the high level of sparsity in the transformation matrixV, hence in the compressed model. These advantages are counterbalanced by a potentially larger number of total basis functions ρ = ν ρ ν with respect to a global unstructured SVD.
2) Hierarchical SVD: The Hi-SVD does not consider each subset of responses M ν for ν = {1, 2, 3} as independent, but it applies a hierarchical strategy to iteratively remove the contribution of a given set from the others by projection. The final result is represented by a block-triangular structure where individual matrix blocks are iteratively defined for = {3, 2, 1} through a truncated SVD withV =Σ V T andV j, =Ū T j M . This representation preserves the same advantages of the BD-SVD compression, but it generally requires a smaller number of basis functions ρ = ν ρ ν to obtain the same level of accuracy.
The same 25-port example of Section III-A is used here to demonstrate the effectiveness of the proposed strategies. Modal scattering responses were subjected to the BD-SVD and the Hi-SVD using truncation thresholds {10 −8 , 10 −6 , 10 −6 } for M 1 , M 2 , and M 3 , respectively. The value for M 1 was chosen to preserve an accurate dc reconstruction for the nonvanishing responses at dc [the diagonal entries of S m (jω)]. Fig. 8 confirms the preservation of dc zeros for all modal responses. All compression techniques provide similar accuracy for all reconstructed scattering elements (see Fig. 9). The top panel of Fig. 10 depicts the selected basis functions for each block with dc zeros of different orders.

IV. STRUCTURED COMPRESSED MACROMODELING
The compressed data representation (21) obtained either with the standard [27] or the proposed dc-preserving structured compression methods of Section III can be used as a starting point to construct a low-complexity structured macromodel. We review the general framework [27] in Section IV-A, and we present the proposed dc-preserving generalization in Section IV-B.

A. Compressed Macromodeling: Background
For each individual basis functionw q available as a column of matrixW ∈ C K×ρ , we construct a rational macromodel in where p n are the system poles (common to all basis functions), r qn are the model residues, and r q∞ is the direct coupling term. The model is computed by processing all basis functions at once through the Fast VF algorithm [9] by enforcing the fitting condition considering that the kth component (w q ) k of each basis vector w q corresponds by definition to the frequency ω k . All basis function models (34) are collected in a row vector denoted as compressed macromodel that can further be cast in a state-space form as with A w ∈ R N w ×N w storing the poles p n in its main diagonal, b w = 1 N w column vector of ones, C w ∈ R ρ×N w collecting the residues r qn , and d w ∈ R ρ collecting the direct coupling constants r q∞ . Standard modifications apply [1] in the case of complex conjugate model poles.
The main advantage of this approach is the computational cost reduction in building the macromodel. Indeed, instead of fitting P 2 responses (or P (P + 1)/2 in case of reciprocal structures), the compressed-macromodel requires only ρ elements (the set of basis functions) to be fitted. This corresponds to a major reduction in terms of model coefficients, i.e., a smaller amount of decision variables both for the fitting and the subsequent passivity enforcement.

B. Structure-Preserving Compressed Macromodeling
One of the main contributions of this work is in the definition of a compressed structured macromodel where the rational model w q (s) for each basis function inherits the features of the corresponding columnw q ofW, in particular in the presence of dc zeros of a given order. For the particular application to shielding enclosures, where we need to model three blocks ν = 1, 2, 3 with ρ ν basis functions each and a dc zero of order 0,1,2 respectively, this is achieved by defining where with frequency-dependent rational weighting factors and m 1 (s) = 1. Assuming the components of w(s) in rational form (34), these high-pass filters enforce dc zeros of appropriate multiplicity while remaining compatible with the requirement that scattering responses must be unitary bounded since |m ν (jω)| ≤ 1 for ν = 1, 2, 3 and ∀ω ∈ R. The poles p ∞ 2,3 are selected to be outside the modeling band to minimize their effects in the modeling bandwidth where ω max is the maximum frequency used to fit the model. Since the asymptotic behavior of w(s) and w(s) is identical so that the model response of w(s) at s = ∞ is controlled by the direct coupling terms r q∞ of (36).

1) Model Fitting:
Here, we show how the filtering term Γ(s) defining the compressed structured macromodel (39) can be accounted for without any modification of a standard VF or FastVF engine [9]. Direct enforcement of model structure (39) would in fact require a modification of the VF basis functions in order to account for the presence of the two additional high-frequency poles p ∞ 2,3 . Instead, we rescale the basis vectors by the inverse of the weighting matrix Γ −1 k = Γ(jω k ) −1 , and we enforce a standard VF fit with a compressed model in form (34). The resulting fitting condition for each rescaled basis vector reads The right-hand side in (43) takes a nonzero finite value at dc due to the inverse rescaling factor that cancels the corresponding dc zeros. Therefore, a standard rational model with structure (34) can be used to perform the rational approximation (43). The final structured model with dc zeros is obtained through multiplication by the filtering terms through (39). Fig. 10 reports the set of basis functions before and after rescaling by (43). In order to improve numerical conditioning and accuracy, each qth basis vector can further be normalized by its 2-norm or RMS value, embedded as an additional weight in the corresponding diagonal element Γ qq (s) of the filtering matrix in (39) and (43).
In addition to the slope, a correct dc level of the model is obtained by enforcing the additional equality constraint while solving (43) in least squares sense using VF or FastVF, where it is assumed that the dc value is available as the first frequency sample k = 1 as s 1 = jω 1 = 0. The compressed structured macromodel (38) can be mapped to the original multiport scattering representation H(s) through a constant algebraic transformation where the mat operator reshapes a vector of P 2 elements into a P × P matrix.

2) Enforcing Model Passivity:
We now address the passivity enforcement of the structured compressed (scattering) macromodel defined in (45), assumed to be asymptotically stable by construction as obtained by VF. We recall that a scattering macromodel is passive if its transfer function is bounded real, which is implied under the working stability assumptions by condition where σ max is the largest singular value of its matrix argument. The scheme that we adopt is the method of choice in most state-of-the-art tools based on iterative local perturbation singular values that violate (46). Let us consider a single frequency ω k at which this passivity condition is violated, and compute the SVD where U k ,V k are unitary, and Σ k stores in its diagonal the singular values of H k . We define σ k > 1 as the largest singular value, with u k and v k the corresponding left and right singular vectors. Following standard results [1], [10], we write a firstorder singular value perturbation where ΔH k is the corresponding model perturbation. Requiring σ k to be less than one leads to the following first-order condition for local passivity at ω k Rather than following a two-step process as in [27] by enforcing asymptotic passivity first, followed by an iterative enforcement loop at all finite frequencies, we follow here a direct approach to perturb concurrently all model coefficients while enforcing (49) at all (finite and infinite) frequencies. First, we define the model perturbation in terms of its coefficients by introducing the vectorized variable where ΔR w is a perturbation to be determined for all model parameters r qn , r q∞ collected in matrix R w ∈ R (N w +1)×ρ , including also the direct coupling constants that are responsible for the asymptotic behavior of the model. A straightforward algebraic calculation allows to write where ϕ k = ϕ(jω k ) ∈ C N w +1 stacks the partial fraction basis functions ϕ 0 (s k ) = 1, ϕ n (jω k ) = (jω k − p n ) −1 . Therefore, the local passivity constraint (49) can be written in a compact form as A standard adaptive sampling-based passivity characterization such as [31] is used to detect the frequencies ω k where passivity violations are located, leading to multiple concurrent constraints (52). Model accuracy during passivity enforcement is preserved by minimizing the following cost function: which is a quadratic form of the decision variables x w . Passivity enforcement is applied in a standard way by minimizing (53) subject to inequality constraints (52) and possibly equality constraints (44) for dc point preservation within an iterative loop. We remark that the number of decision variables in x w is only ρ · (N w + 1) since the compressed macromodel form is perturbed, although the passivity of the large-sized P × P macromodel is being enforced. Therefore, the complexity of the quadratic optimization problem that is solved at each passivity enforcement iteration is not affected by the number of I/O ports of the structure under modeling but only by the number of basis functions. Note also that the filtering matrix Γ(s) that defines the proposed compressed structured macromodel with dc zeros is applied as a weight since its poles p ∞ 2,3 are fixed so that no additional unknowns are introduced.

V. NUMERICAL RESULTS
In this section, we report several results that confirm the effectiveness of the proposed macromodeling strategy. All numerical results have been obtained using a prototypal MATLAB code on a Workstation based on Core i9-7900X CPU running at 3.3 GHz with 64 GB of RAM.

A. 25-Port Box Example
We report the final result for the 25-port box that was used as a running example to demonstrate all intermediate steps of the proposed macromodeling flow. Here, two compressed macromodels are constructed and compared. The first is built following [27] using the set of unstructured basis functions reported in Fig. 7. The second is built on the basis set of Fig. 10 following the proposed approach of Section IV, in particular by enforcing in the model structure in the presence of dc zeros of order 1 and 2 wherever appropriate. Both models are trained in the scattering domain with a fixed number of polesn = 89. A comparison of these models is presented in Fig. 11, where the main systems representations (scattering, admittance, and impedance) are plotted, and the latter two are obtained by postprocessing.
The proposed structured model has been processed through a passivity enforcement scheme, as described in Section IV-B2, that resulted in a passive model in 53 iterations and an overall elapsed time of ≈ 290 s. The maximum singular values before and after the proposed passivity enforcement are reported in Fig. 11(a). Fig. 11(b)-(d) confirms that the proposed structured approach outperforms the reference in terms of full-band accuracy in all three representations. As an additional reference for comparison, the model obtained by a standard approach [24], without any preprocessing except addition of the physics-based dc point, is completely wrong in the low-frequency gap (see Fig. 3).
To further emphasize the importance of model behavior in the low-frequency gap Ω G , Fig. 12 compares the transient simulation results obtained using the proposed model and a standard VF model (the same as depicted in Fig. 3). This result confirms that an inappropriate modeling flow at low frequencies has a strong impact also when performing time-domain simulations. Fig. 12. Transient simulation of the proposed model compared to a standard VF model. The excitation field is a double exponential (HEMP) waveform [32] defined as E inc = E 0 (e −at − e −bt )u(t) with a = 40 × 10 6 s −1 , b = 476 × 10 6 s −1 and amplitude E 0 = 50 kV/m. Fig. 13. As in Fig. 11, but for a 100-port shielding enclosure.  Fig. 11, but for a 400-port shielding enclosure.

B. Scaling up to 400 Ports
In this section, we summarize the results of the proposed procedure on several shielding enclosures of the increasing size of the regular p × p series/parallel grid, corresponding to a number of ports P = p 2 . We focus, in particular, on grids with p = 8, 10, 15, 20 branches resulting in systems with P = 64, 100, 225, and 400 ports, respectively.
The regularization and extrapolation procedure of Section II-C was applied to obtain a set of (modal) scattering responses in the range [0, 1] GHz. Regularization resistances were set to r = 0.1 Ω and R = 100 MΩ in order to minimize the perturbation of the original MoM data. Different values of resistors with R ∈ [10 6 , 10 9 ] Ω and r ∈ [10 −3 , 1] Ω were also tested, with no  practical impact on the overall modeling process except for a slightly different dc response and induced perturbation amount on the original data. The BD-SVD compression strategy was then applied with truncation thresholds of {10 −8 , 10 −6 , 10 −6 } for M 1 , M 2 , and M 3 , respectively. Note that the threshold on M 1 was determined compatibly with the estimated field solver accuracy in order to guarantee a consistent reconstruction of the dc matrix after compression.
A summary of the proposed data compression results for all test cases is reported in Table I. This table shows that increasing the number of ports leads to a drastic compression in terms of basis functions ρ, already considering as a reference only the upper triangular part of the transfer function with P (P − 1)/2 responses. Although the proposed structured compression approximately doubles the number of basis functions ρ with respect to the standard compression of [27], the required elapsed time is practically not affected.
The procedure of Section IV was applied to build a set of macromodels preserving the LF content of the data. In particular, the adaptive sampling method presented in [31] was used to create the constraints (49) required by each iteration of the passivity enforcement. The results are summarized in Table II. All models are highly accurate. As far as runtime, the rational fitting phase shows almost no dependence on the number of ports since mainly affected by the number ρ of compressed basis functions. Most of the execution time is required to enforce macromodels passivity. For the adopted passivity enforcement scheme, this time cannot be predicted a priori since depending on the particular singular value trajectories and on the number of required iterations.
Figs. 13 and 14 provide a comprehensive report on model accuracy in all system representations (scattering, admittance, and impedance) in the physical domain after conversion from the modal reference system used for compression and models identification. Results are shown only for P = 100 and P = 400 since nearly identical results were achieved for P = 64 and P = 225. All examples show a major improvement in the accuracy with respect to the models obtained with [27], especially at low frequencies and dc.

VI. CONCLUSION
In this article, we provided a general methodology to construct rational macromodels of low-loss electromagnetic systems having a large number of interface ports, starting from possibly incomplete frequency characterization due to a low-frequency gap in data samples from a full-wave solver. This situation is known to be extremely challenging for macromodeling applications, with a low-frequency sensitivity to approximation errors that are exacerbated by the low-loss nature and by a large number of ports. The proposed strategy combines data preprocessing involving regularization and extrapolation in a suitably defined physics-based asymptotic modal domain, a structured data compression based on a customized SVD, and a structured rational fitting process based on compressed data. The result is a robust framework that is able to preserve full-band accuracy in the model down to dc while minimizing its sensitivity to loading conditions, including the limit cases of short and open terminations.
The derivations and the reference application examples considered in this work are conducting enclosures intended for energy-selective shielding. This application requires loading with nonlinear elements, a possibly large number of ports placed as a grid throughout shield apertures. For this reason, a reduced model sensitivity to changes in loading conditions is essential for successful transient analyses aimed at nonlinear shielding effectiveness assessments. Beyond this specific application, the proposed regularization, extrapolation, and compression approaches are general and can be applied to different application scenarios. Future investigations in this direction will involve an application to system-level power integrity modeling, which is another domain where sensitivity to loading conditions has been identified and still requires an adequate systematic solution.