Addressing Load Sensitivity of Rational Macromodels

Behavioral models are effective tools used to relieve the computational burden of large-scale system-level simulations. In electrical and electronic applications, the vector fitting (VF) iteration often represents the algorithm of choice for generating low-order equivalent circuits for complex multiport components in a data-driven setting. Although accurate and reliable in general, macromodels generated via VF are inherently represented in terms of a rational approximation of one specific input–output transfer function of the structure under modeling, for example, its scattering matrix. However, accuracy in the scattering representation does not necessarily imply good accuracy when solving the macromodel in a system-level setting, under different termination conditions. In fact, the sensitivity of the macromodel with respect to its loading conditions may be large and needs to be addressed and controlled. In this work, we present a modified VF scheme that overcomes this issue, by introducing in the rational approximation algorithm the requirement that the macromodel remains accurate when interconnected with a known class of admissible networks. The proposed formulation is based on an augmentation of the cost function minimized at each VF iteration; furthermore, it does not require additional expensive data-gathering steps when compared to standard approaches. The effectiveness of the scheme is tested over a set of relevant examples, in particular, for power integrity applications.

always extremely small, ensuring that the model accurately reproduces the prescribed target transfer function.
Since its original formulation [1], several VF improvements have been documented, including efficient methods for handling large electrical systems [6] and related implementation to multicore computing hardware [7], [8], increasing the robustness to noisy data [9], [10], and improving the convergence properties of the iteration [11], [12].A comprehensive overview of the VF algorithm and its applications is available in [13], including postprocessing algorithms for passivity enforcement and SPICE equivalent network synthesis.
The main purpose of a macromodel generated via the VF algorithm is to represent the underlying structure in systemlevel simulations, performed either in the time or in frequency domain.The ultimate goal of these simulations is to predict the behavior of a large interconnected system, in which the macromodel plays the role of a subnetwork.In this view, one relevant and often disregarded issue affecting state-ofthe-art macromodels is represented by the sensitivity of their accuracy on the loading conditions occurring in such a simulation.As already pointed out in previous works [14], [15], [16], [17], even macromodels that accurately fit one specific transfer function of the structure under modeling may return inaccurate predictions of the reference port behavior when they are interconnected with different loads as requested by the application at hand.This happens because an even small residual error in the rational fitting stage can be magnified by the feedback interactions between the model and its actual terminations.In many situations, this magnification can be so pronounced to completely invalidate the simulation results.
In the available literature, relevant contributions have been proposed to reduce the macromodel sensitivity to those terminations defining immittance or scattering representations.This was done by applying inverse magnitude data weighting to the rational fitting cost function [18], possibly preprocessing the data points through a mode revealing transformation (MRT) [14], allowing the model to reconstruct the eigenvalues of the target response with good accuracy, regardless of their magnitude.This significantly increases the robustness of the modeling error to changes in the electrical representation.The more general case in which the actual loading conditions can include one arbitrary (but known) termination network has been studied in [15], where the authors introduced the so-called Modal VF (MVF) scheme and addressed an approach to control the accuracy of the macromodel under the prescribed load.Also, in [16] and [17], a tailored error weighting approach was proposed to compensate for the macromodel sensitivity, both during the generation and the passivity enforcement postprocessing stages.
In this article, we introduce a novel formulation of the VF scheme, that allows to increase in the accuracy of the reduced order network when it is interconnected with an arbitrary set of LTI terminations.In the proposed approach, the rational fitting cost function that is minimized during the VF iteration is augmented to incorporate the information regarding the port behavior of the system under modeling when it is subjected to the admissible class of terminations.To do so, we evaluate the frequency-domain spectra of the port variables of the structure to be modeled under the prescribed (possibly multiple) loading conditions.Such data is then suitably incorporated in the VF cost function, in addition to the standard least squares rational fitting error.This simple strategy is sufficient to guarantee robust accuracy under several different loading conditions, thus resulting in a drastically reduced load sensitivity.As will be shown in our derivations, the proposed framework can be interpreted as a generalization of the MVF scheme [15], with which it shares some implementation aspects.
The effectiveness of the proposed approach is experimentally validated over a collection of test cases of practical relevance in the field of power integrity.

II. NOTATION AND BASIC DEFINITIONS
We denote scalars with normal fonts (x), vectors with lowercase bold fonts (x), and matrices with uppercase bold fonts (X).The transpose of X is X T .The identity matrix of size P is denoted by I P and we define R 0 = R 0 •I P for any R 0 ∈ R. The operator | • | represents the absolute value, intended elementwise when applied to vectors or matrices.With A ⊗ B, we denote the Kronecker product between matrices A and B, while e i is the ith canonical basis vector of R n , with n clear from the context.
With reference to the generic LTI P-port network labeled as device under test (DUT) in Fig. 1, we denote with v ℓ and i ℓ for ℓ = 1, . . ., P the voltage and the current at the ℓth port, respectively.The incident and reflected scattering power waves are defined accordingly with respect to a positive reference resistance R ℓ 0 .Unless otherwise specified, we assume the standard scattering reference system with R ℓ 0 = 50 .All port voltages, currents, and scattering signals are collected in the vectors v, i, a, and b.
External characterizations of the network can be obtained in terms of its port transfer matrices H(s) ∈ C P×P , where s ∈ C is the Laplace variable.This implies that, for each port ℓ, one variable is defined to act as local input u ℓ , and the dual variable is the local output x ℓ .Collecting inputs and outputs in vectors u and x, we have The most common choices are the impedance, admittance, and scattering representations, for which the same choice of inputs and outputs applies to all ports, corresponding, respectively, to Impedance, admittance, and scattering matrices are related by the following transformations: and corresponding inverses.
For the P-port DUT structure of Fig. 1, we define a port-behavior as a particular configuration of the vector pairs [v(s), i(s)], [a(s), b(s)], and more generally [u(s), x(s)] that is compatible with (3)-( 5) or (2).The set of admissible portbehaviors collect all port behaviors of the structure obtained by arbitrarily setting the inputs and constraining the corresponding outputs via the associated transfer functions.

III. MACROMODELING OF SUBNETWORKS
The DUT in Fig. 1 represents a possibly complex electromagnetic structure, accessible through a set of P well-defined electrical ports.No information will be assumed about the internal structure of the DUT.However, we will assume that one of its port behaviors is available in terms of the frequencydomain samples of the associated network function, denoted as where the accent ˘labels the "true" system response.The latter can be obtained either through real or virtual measurements provided by a high-fidelity circuit or electromagnetic solver.
Based on these available measurements, the common objective of macromodeling is to identify a low-order equivalent circuit whose transfer function H(s) accurately reproduces the response of the high-fidelity model Most electrical and electronic applications require this accuracy both when the model is considered as a stand-alone unit, and also (especially) when it plays the role of a subnetwork in a larger electrical interconnected system, as shown in Fig. 1.This latter requirement is often overlooked in the vast literature available on the subject: typically one is satisfied when a scattering macromodel S(s) accurately reproduces the scattering responses Sk from a field solver by minimizing the fitting error in (9).However, usage of the model under different operating conditions may lead to an amplification of the inevitable approximation errors involved in the rational approximation process, up to a point where the port behavior in such operating conditions becomes corrupted [14], [15], [16], [18].This article intends to overcome the above limitations.

A. Sensitivity to Transformation and Loading
Let us focus on a common scenario where the characterization of the DUT is through sampled scattering responses Sk and a scattering macromodel S(s) is obtained by fitting such data.We have where δS k represents the fitting error at frequency ω k .Assume now that the model is converted to the admittance representation using the transformation (7).The result can be written as where Y( jω k ) is the impedance model response, Yk is the exact DUT impedance, and δY k is the corresponding approximation error at frequency ω k .This situation corresponds to terminating the scattering model with ideal voltage sources and observing as outputs the corresponding port currents.The error on such voltages is expressed as a linear combination of the elements in δY k .When both δS k and δY k are sufficiently small, one can use a first-order sensitivity analysis to obtain a linear relationship where the tensor J coincides with the Jacobian matrix of the transformation (7) up to elementary reshape operations.Unfortunately, the partial derivative elements in this tensor may be large and even unbounded due to the matrix inversion operation involved in (7).As a result, the induced error δY k may be large.An example is provided in Fig. 2, where a small perturbation in one scattering matrix element leads to a large induced error in the associated admittance matrix.The same considerations apply to the conversion to impedance [see Fig. 2 (bottom)].
The above example is just a canonical and basic demonstration of a very critical situation for all those modeling and simulation flows that are based on surrogate, behavioral, and approximate models.Such models are accurate only in the representation used for their training.When the models are simulated under different conditions, approximation errors may grow to large and unacceptable values.
This sensitivity problem arises not only when performing conversion in model representation to the standard admittance or impedance representations.Anytime the model ports are terminated into loads that are different from the conditions used Fig. 2.
Effect of macromodel sensitivity to change of representation.The structure under modeling is a power distribution network (PDN) with P = 18 ports.A standard VF macromodel for the structure is generated by fitting the available scattering matrix (inverse-magnitude weighting, 12 poles).The resulting model accurately reproduces the data (top).Changing the macromodel representation from scattering to admittance (middle) or impedance (bottom) leads to large error amplification in some frequency bands.The proposed algorithm fixes the above sensitivity issues and provides accurate results in all input-output representations.
for training (50 for standard scattering models), the same difficulties are to be expected.This is the common scenario in CAD flows, where macromodels are used as components in large complex system-level simulations, together with many other connected subsystems.
Two main application settings where the sensitivity problem may impair modeling and simulation flows can be considered.
1) The terminations that will be connected to the macromodel ports are known and fixed, or they belong to a well-defined class.A typical example is the set of ports of a PDN at the board or package level that are to be terminated with decoupling capacitors.Such capacitors have a well-defined frequency response (a pole at dc, a pole at infinite frequency, and one resonance frequency).In this situation, the macromodel can be optimized such that the specific information that is available about the terminations is used to tune the model accuracy when the terminations are connected.2) A dual situation is the case of a completely unknown termination scheme.The model is supposed to be robust for general and unspecific use, and therefore accuracy should be guaranteed under any possible termination condition.Rather than a specific termination, as wide as possible set of termination conditions should be used to enforce and test model accuracy.

B. Proposed Approach
In this article, we introduce a rational approximation procedure aimed at generating macromodels that accurately reproduce the port behavior of the DUT, under possibly wide classes of target loading conditions.Since all the derivations are here conducted in the frequency domain, the admissible loads are here restricted to LTI subsystems.
The proposed algorithm, derived from the standard VF iteration, is driven by a generalized approximation condition that forces the model to reproduce the DUT port behaviors when it is interconnected with the prescribed class of loads.
The key steps of the approach are summarized as follows.
1) As in standard macromodeling approaches, we represent the reduced order equivalent circuit in terms of one of its port representations (e.g., scattering), and we collect the DUT response samples (8) corresponding to this representation.2) We perform a collection of experiments aimed at sampling the value of the DUT port variables when it is interconnected with the prescribed loads.These experiments are performed by placing suitable sources at the interface ports of the fully interconnected system.This procedure is described in Section IV. 3) We employ the rational VF model structure to parameterize the unknown network function of the model, and we perform the VF iteration with a modified cost function, which augments (9) by additionally requiring that the observations collected at previous step are accurately reproduced by the model.This modified VF iteration is presented in Section V.

IV. DUT BEHAVIOR UNDER MULTIPLE TERMINATION CONDITIONS
In this section, we characterize the port behavior of the DUT in terms of a given I/O representation, when it is considered as a subnetwork of a larger interconnected system, subject to external excitations.We consider the generic interconnection template depicted in Fig. 1, where the DUT is connected with the loading networks L = {L t , t = 1, . . ., T }.We assume that each of these termination networks L t is an LTI system, with P t ports connected to the DUT and P ′ t ports closed on ideal voltage or current sources (see Fig. 1 for a graphical illustration).Since we are not interested in discriminating the type of sources, we will collect all voltage and current sources exciting L t in vector w t .Let us define with a t and b t the vectors of power waves that are incident into and reflected from the internal ports of each termination network L t .By superposition, the frequency-domain characterization of L t can be written as b t = t a t + P t w t (13) where t and P t are the transfer functions from a t to b t and from w t to b t , respectively.Throughout this section, we omit the frequency argument s = jω.
All the termination networks can be grouped in a single global multiport T , represented by the dashed box of Fig. 1.The network T has P internal electrical ports interconnected with the DUT, and an additional set of P ′ = t P ′ t ports that are accessible from the outer environment.Collecting all individual characterizations (13) leads to a global characterization of T as where with operators blkdiag and col stacking their arguments in a block-diagonal matrix and as a block-column vector, respectively.
We now assume that the DUT is characterized in terms of its scattering matrix S and that exact closed-form descriptions are available for the terminations in terms of , P, and sources w.In such conditions, the values of all electrical variables at the internal ports can be computed algebraically.In fact, coupling (14) with the additional interconnection constraints provides the system of linear equations whose solution readily yields the desired port variables More generally, in case a different representation H is chosen to represent the DUT as in (2), the above procedure is still applicable, provided that one expresses the characteristics of the termination networks using the interface variables u, x as with appropriate transfer functions and Q.We obtain whose solution is straightforward Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The pairs (a, b) in ( 17) and (u, x) in ( 20) are admissible portbehaviors of the DUT, since by definition related by the DUT transfer functions S and H, respectively.

A. DUT Characterization Under Multiple Operating Conditions
To ensure the robustness of the DUT model under multiple operating conditions induced by several different termination schemes, we can apply the above procedure to evaluate the set of corresponding induced admissible port behaviors.Let us define a set of M termination/excitation schemes We evaluate the solution ( 20) of ( 19) for each m, obtaining a set of admissible port behaviors These represent the exact solution of the interconnected system in various loading conditions.Our main objective is to make sure that, when a macromodel H(s) replaces the true system H(s), the corresponding port behaviors remain accurate.
As discussed in Section III-A, this is not guaranteed in general, unless a set of dedicated constraints is applied.

V. TARGETING MACROMODEL ACCURACY IN REAL OPERATING CONDITIONS
In this section, we show how a simple modification of the basic VF scheme can be performed for ensuring macromodel accuracy in real operating conditions.

A. Formulation
We assume a standard rational macromodel in the poleresidue form where H(s) ∈ C P×P is a prescribed model network function (e.g., scattering), { p i } are the unknown model poles, and {R i } the corresponding unknown residues with direct coupling R ∞ .The key idea of the proposed formulation is that the macromodel (23) should not only approximate the corresponding I/O representation of the DUT, but it should also reproduce accurately the DUT port behavior when interacting with prescribed loads.This goal is achieved by imposing the approximation conditions Condition ( 24) is standard, while (25) is introduced to increase the model robustness and reliability.The pairs [ ȗ(m) , x(m) ] are obtained by solving (19) for M different target terminations T (m) over the desired frequency band and are always admissible, such that Thus, (25) requires that the model behaves as the DUT when interconnected with the loading networks T (m) .
In practice, we only assume that the high-fidelity transfer function H( jω) can be evaluated on a discrete grid of K frequency values {ω k } k=K k=1 .Similarly, ȗ(m) ( jω) and x(m) ( jω) are also known through their samples ȗ(m) ( jω k ) and x(m) ( jω k ), k = 1, . . ., K .Using the pole-residue form (23), conditions (24) and ( 25) are particularized at the sampled frequencies as where we collected all the conditions (25) in the compact matrix form In ( 26) and ( 27), both the residues and the poles of the model are unknown.We find them by suitably modifying the standard VF algorithm.Following [11], we rewrite the model in (23) in a barycentric form using a set of ν auxiliary basis poles {q i } ν i=1 and two rational functions N(s) and d(s) Combining with (26) and ( 27) and multiplying both sides by d(s) leads to to be interpreted as a coupled linear least-squares problem in the unknowns N i and d i .The VF iteration repeatedly solves (31) for {d i } and updates the set of basis poles with the zeros of d(s) before setting up the next iteration.In the VF algorithm, this pole relocation procedure ends when the basis poles stabilize and the denominator d(s) → d 0 (constant).Then, the model poles { p i } are taken to be the basis poles at the last iteration, and the model residues R i are found as the solution of the following additional linear least-squares problem derived from (31): Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where F denotes the Frobenius norm.The two objectives in (26) and ( 27) have been combined in (33) into a single cost function and their relative importance is weighted using the penalization parameter λ 2 .This hyperparameter provides increased flexibility in case the numerical magnitudes of the two objectives are substantially different.In our numerical experiments, it was found that fine-tuning was not required in most cases, for which a balanced weighting with λ = 1 was sufficient.The example described in Section VI-B is an exception because entries of the Z-parameters are about two orders of magnitude smaller than the S-parameters in the lowfrequency range.To compensate for this numerical difference, the parameter λ was set to 10 3 .We now particularize the above general formulation to a few practical scenarios that are common in applications.

B. Fitting Different Representations Simultaneously
Let us consider the situation H(s) ≡ S(s), so that the model is represented and computed in scattering representation.In Signal and Power Integrity applications, this is by far the most common setting due to the increased numerical robustness associated with the boundedness of scattering responses.
Assume further that all DUT ports are loaded by ideal current sources, with T in Fig. 1 including only ideal connectors.We compute the set of port behaviors associated with the following excitation settings: |e p , e p ∈ R P , p = 1, . . ., P.
This excitation scheme corresponds to the evaluation of individual columns of the impedance matrix representation of the DUT.A straightforward calculation leads to the following expressions for the matrices Ȗk and Yk , collecting, respectively, the incident and the scattered waves at the DUT ports for all excitations: and the approximation condition (27) for a single-frequency sample translates into The above represents a linearized version of the condition obtained by multiplying the first and the second terms of (37) by [I P −S( jω k )] and normalizing twice by R 1/2 0 .Note that the right-hand side of (37) is indeed the impedance matrix Z( jω k ) of the model as would be obtained using the conversion ( 6), but which is never computed directly.Therefore, condition (37) implies that the impedance matrix of the model approximates that of the DUT.A similar procedure gives the choice of weights Ȗk and Xk to optimize the admittance parameters The approximation condition resulting from this choice is and corresponds to a linearized version of (7).
The enforcement of either condition (36) or (38) during the construction of the macromodel is expected to reduce its sensitivity to a change of representation from scattering to impedance or admittance.When enforcing simultaneously the constraints (36) and (38) in addition to the standard VF cost function, the resulting model will be accurate both under conversion to impedance and admittance, or equivalently to termination of all ports with open or short circuits.
A numerical experiment showing the effectiveness of this procedure is reported in Fig. 2. The scattering, admittance, and impedance representations were optimized simultaneously via (36) and ( 38), and the model order was chosen again as ν = 12.The resulting macromodel shows much higher robustness to the change of representation when compared with a standard VF macromodel of the same complexity.

C. Fitting the DUT Modes
The proposed approach can be interpreted as a generalization of the MVF scheme.The MVF algorithm, introduced in [15], is known to effectively improve the robustness of the macromodel to changes in I/O representations.This is achieved by fitting a target network function H(s) in the modal domain, by applying an inverse-magnitude error weighting scheme, such that the model H(s) accurately reconstructs the eigenvalues of H(s), irrespective of their magnitude.
It is straightforward to verify that MVF can be interpreted as a particular implementation of the proposed approach corresponding to the specific choice with ( jω) (diagonal) and T( jω) defining the spectral factorization (eigendecomposition)

D. Practical Implementation
This section discusses the practical implementation of the proposed algorithm.In particular, since the major difference with respect to VF is in the definition of the cost function, the focus will be on translating the optimization problems in (31) and (33) in a matrix form suitable for numerical computations.
We start by rewriting the first term in (33) as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Let us now introduce a row-wise view of the transfer function H( jω k ) and the data matrix Xk Exploiting the block-diagonal structure of (I P ⊗ ȖT k ), a straightforward algebraic manipulation reveals that the cost function J can be stated as the sum of P terms, each depending on a single row of H( jω k ), as highlighted by the outer sum in the following expression: Collecting all the frequency samples as follows: we can write (43) in the more compact form This expression is now particularized to the two VF steps, namely the pole relocation and the final residue computation.We start with the latter residue computation step.When poles are known and only residues need to be computed, we should introduce explicitly such residue unknowns; this can be achieved by combining (44) with (23).By defining . . . where and r p,q is a vector containing all residues and the direct coupling coefficient of the ( p, q) entry of H( jω k ), we can compactly write h p = Pρ p so that Minimizing J is a linear least-squares optimization problem that can be further broken into P subproblems where each ρ p is computed independently.The formulation of the pole relocation iterations described in (31) requires the elimination of the residue unknowns embedded in the numerator function N( jω) in order to formulate a least-squares problem in the denominator unknowns d = (d 1 , . . ., d ν , d 0 ) T .This is here achieved through the standard QR factorization approach of [6], which allows to express ρ p = p d with an appropriate constant matrix p .By defining we then obtain the VF pole relocation cost function In this work, we adopt the so-called relaxed version of VF [11], which eliminates the trivial solution d = 0 of (47) through an appropriate nontriviality constraint.The latter is usually written as c T d ≈ α and added as a penalty term to (47), obtaining This cost function is minimized in the least-squares sense at each VF pole relocation iteration.

E. Computational Cost
Since the pole relocation phase is iterative, the overall computational cost will depend on the number of iterations required before convergence, which is difficult to state in general as in the standard VF algorithm.Hence, we limit the discussion to the per-iteration cost of the pole relocation phase, which represents the step where most of the time is spent [7].
Minimization of the pole relocation cost function in (47) can be implemented through the QR factorization as described in [6].The size of the matrix to be factorized is 2M K × (P + 1)(ν + 1) and this operation is performed P times in every iteration [one for each term in the sum in (47)].Since the cost of QR factorization of an m × n matrix is asymptotically O(mn 2 ), the total cost for all QR factorizations is O(M P 3 K ν 2 ).A final least-squares is required to find the denominator coefficients d, with cost O(Pν 3 ).These two contributions sum to O(M P 3 K ν 2 + Pν 3 ).The QR factorization in standard VF is instead performed on a matrix of smaller size 2K × 2(ν + 1) and in that case this operation is repeated P 2 times in every iteration, resulting in O(P 2 K ν 2 ).Moreover, the final LS problem requires additional O(P 2 ν 3 ) operations, thus giving a total O(P 2 K ν 2 + P 2 ν 3 ).To summarize, the dominant asymptotic cost of each pole relocation iteration is larger by a factor M P in the proposed formulation with respect to basic VF, as for the MVF algorithm [15].In addition to the above asymptotic estimates of the dominant cost, any practical implementation of both VF and the proposed scheme requires additional steps related to the construction and manipulation of matrix coefficients.These operations have a significant impact on the overall runtime, as we will see in Section VI-E.

F. Stability and Passivity
Model stability can be enforced exactly as in the legacy VF algorithm, namely by flipping poles at each pole relocation iteration.In fact, the proposed modification only affects the cost function defining the model accuracy.In particular, the pole relocation step in (32) is inherited from standard VF, where stability is enforced by replacing all right half-plane zeros of d(s) with their mirror images with respect to the imaginary axis (pole flipping).This part of the VF algorithm is still present in the new formulation, leading to guaranteed stable models by construction.
Considering model passivity, any of the existing and wellestablished passivity enforcement algorithms based on model perturbation can also be enriched with the proposed formulation (33).In fact, all these algorithms use passivity constraints during an optimization process that minimizes the model perturbation based on a prescribed cost function.The constraints can be based on singular value or Hamiltonian eigenvalue perturbation; alternatively, they can be cast as Linear Matrix Inequality constraints arising from Kalman-Yakubovich-Popov lemma (see [13] for a detailed overview).Any of these formulations is independent of the particular cost function that defines model accuracy during perturbation process.Therefore, the same error control that is demonstrated in this article for the model construction can be used as a cost function during passivity enforcement.

A. Plane Parallel Plates
In this first academic example, we consider a 2-D distributed structure, representing a canonical PDN in a printed circuit board.The structure consists of two parallel rectangular planes of size l x = 12 cm, l y = 10 cm, separated by a distance d = 1 mm.The dielectric material between the planes has a relative permittivity ϵ r = 4.7 with loss tangent tan δ = 0.01.Five ideal (zero-width) ports are placed between the top and bottom planes at coordinates {(9.8, 1), (10.9, 2.8), (1.5, 5.5), (11, 9.6), (7.6, 9.6)} (units in cm, measured from bottom-left corner).Port #5 is terminated by a 20 m resistance.Characterization of the impedance matrix at the remaining P = 4 ports is available at K = 2000 logarithmically spaced frequencies in the bandwidth 100 kHz-2 GHz.Such data were obtained by an analytical solution of the related field equations through a modal expansion [13], [19], [20], [21], [22], [23], [24], [25].
The objective of this test is to demonstrate that the proposed approach is able to concurrently fit the impedance of the PDN, while ensuring at the same time that such impedance remains accurate when connecting at the four ports sets of series R LC terminations, modeling decoupling capacitors.To this aim, we set H(s) ≡ Z(s) in ( 23), and we performed a total of M = 48 experiments to analyze the port behavior of the structure when it is terminated with the prescribed loads.The experiments were performed as follows.We drew 12 pseudorandom combinations of R, L , C values for the four ports with R ∈ [0.5 m , 1 m ], C ∈ [5 pF, 1 nF], Fig. 3.
Magnitude of the input impedance of the reterminated parallel plane structure described in Section VI-A.Reference data are compared with the responses of the models resulting from the application of the proposed approach and the plain VF algorithm with and without MRT.
L ∈ [5 pH, 10 pH] and considered the corresponding 12 reterminated systems.For each one of these, four experiments were performed where a unit current is injected in each port in turn.In these 48 experiments, we collected the currents and voltages at the DUT ports (resp., Ȗk , Xk ∈ C 4×48 ) and optimized a model for the impedance parameters Z(s) with the objective Z( jω k ) Ȗk ≈ Xk , k = 1, . . ., 2000.We also constructed a macromodel of the same impedance matrix using VF with relative error weighting, and a third macromodel where standard VF was combined with MRT [14].In all cases, the number of model poles was set to ν = 22.
A first comparison is shown in Fig. 3, depicting one of the driving-point impedances obtained by the three models after loading the computed model with one combination of R LC terminations used in the training.We see that both standard VF and VF with MRT fail in predicting the impedance response around one of its resonance peaks.Failure of the MRT can be explained by the fact that the related modal transformation diagonalizes the transfer matrix and guarantees that its eigenvalues are adequately approximated only at a single frequency and not over the complete modeling bandwidth.Conversely, the proposed approach is specifically designed to enforce accuracy at all frequencies available in the training dataset.
To further characterize sensitivity to terminations within the given R LC class, we tested these macromodels on 200 random R LC terminations, different from those used in the training.The error of all transfer matrix entries with respect to the exact frequency response was computed for all frequencies, and the corresponding rms value was extracted and then averaged over the 200 test experiments.The results for the proposed, standard VF and VF with MRT models are depicted in Fig. 4.This experiment confirms that, with the same model complexity (order), the proposed algorithm is able to reduce significantly load sensitivity with respect to competing approaches.

B. Real Power Distribution Network
This second case study is a PDN model coming from a real product design (courtesy of Yan Shen Fen, Intel), first introduced in [17].The structure under modeling has a total of P = 18 ports, and a high-fidelity characterization of its behavior is provided in terms of K = 602 samples of its scattering matrix.The objective is here to model the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.structure while optimizing its accuracy when it is simulated concurrently with a termination network T such that the following holds.
1) The first port is connected with a 2 m resistor, representing the output resistance of an onboard regulator.and computing in all cases the samples of the resulting incident and the reflected waves at the ports of the PDN.These samples were used to assemble matrices X k and U k used to enforce condition (27) during the model generation.We then built a model of order ν = 12, using λ = 10 3 in the VF costfunction (33).Fig. 5 (top) demonstrates the accuracy of the resulting model in its native scattering representation.Both the proposed algorithm and the standard VF iteration produced models, respectively, S(s) and S VF (s), with excellent accuracy: both model responses are not distinguishable from the training data.This is a further confirmation that standard VF is indeed very reliable in fitting rational models to data.Fig. 5 (bottom) shows, however, that the weighting induced by the proposed modification has an effect on the frequency-dependent fitting error in selected frequency bands.
Then, we reterminated both S(s) and S VF (s) with the network T , and we computed the impedance matrix of the resulting loaded model.As shown in Fig. 6, the model S VF (s) generated by standard VF provides a quite incorrect prediction of the loaded impedance, as a result of the model sensitivity to its terminations.This sensitivity is dramatically reduced for the model S(s) generated by the proposed approach.This result is partly given by a different outcome of the pole relocation phase in the proposed formulation with respect to standard VF.Impedance responses resulting from retermination of the PDN structure described in Section VI-B, whose scattering parameters are shown in Fig. 5.The S-parameters models given by the standard VF and the proposed method (reported in Fig. 5) are also reterminated, with the former showing significant deviation due to error magnification arising from interconnection.Fig. 7 depicts the location of the poles for both models (only the top-left quadrant of the complex plane is reported).We see Fig. 7. Model poles given by standard VF and by the proposed algorithm for the example described in Section VI-B, whose responses are reported in Figs. 5 and 6.Fig. 8.
Sample admittance response of the example described in Section VI-C.The gray area collects the cloud of points arising from the conversion to the admittance of 500 independent scattering responses affected by up to a 0.04% relative perturbation; the blue curve represents the exact admittance response.that only some poles are common to both models, with four complex poles being moved to substantially different positions on the real axis.

C. Package Structure
The next example focuses on model sensitivity to changes in port representation.The structure under investigation is a package characterized through K = 317 samples of its scattering matrix.With respect to [16], where it was originally introduced, we considered a subset of 16 ports.A preliminary experiment was carried out to show that the scattering parameters data of this structure are sensitive to conversion to admittance representation.Starting from the original scattering data, we added a zero-mean random noise η, characterized by a uniform distribution in the interval [−10 −4 , 10 −4  ].This noise was applied as a relative perturbation to all frequency samples of all scattering responses as Si j ( jω k ) ← Si j ( jω k )[1 + η i j ( jω k )].Then, the noisy data were converted to admittance representation.This process leads to a noise amplification, as indicated in Fig. 8, where the gray-shaded area represents the frequency-dependent interval collecting the cloud of points arising from a total of 500 independent noise realizations.This area is centered around the blue curve, which corresponds to the nominal admittance response.This result provides a numerical verification of the firstorder sensitivity analysis (12), confirming that: 1) the error amplification is frequency-dependent and 2) this amplification can be extremely large (the low-frequency band in Fig. 8).Admittance parameters of the package example described in Section VI-C.These responses were obtained by converting to admittance the reference data and the corresponding responses of the two models obtained by VF and by the proposed approach.
We modeled the structure by applying the proposed approach with ν = 20 poles, optimizing several representations simultaneously.In particular, we required the model to accurately reproduce the reference admittance, impedance, and scattering parameters with both R 0 = 50 and R 0 = 0.01 .Also, we ran the plain VF iteration to obtain a standard model for the scattering parameters, with the same number of poles, and implementing an inverse-magnitude weighting scheme.
Fig. 9 shows the scattering response S 1,1 , for which both algorithms yield accurate models, as confirmed by the corresponding error plots depicted in Fig. 10.On the other hand, Fig. 11 indicates that, while the residual error in the VF model is magnified after changing the port representation from scattering to admittance, the sensitivity of the proposed model is significantly reduced.This example further confirms the effectiveness of the proposed modified cost function and associated fitting algorithm in producing models that have Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 12. Top panel shows one representative scattering response of the ribbon cable described in Section VI-D.A sample admittance response obtained by conversion of scattering responses of models from both standard VF and the proposed algorithm is depicted in the middle panel and compared to the exact admittance response.The bottom panel is a zoom of the middle panel to highlight the improvements provided by the proposed approach.enhanced robustness to changes in port representation with respect to standard VF.

D. Ribbon Cable
This example is a ribbon cable made up of eight identical wires aligned in a row.Each wire has a copper circular conductor (radius r w = 0.5 mm and conductivity σ = 5.8 × 10 7 S/m) with a dielectric coating (radius r d = 0.8 mm and relative permittivity ε r = 4.2).The center-to-center spacing between adjacent wires is D = 1.61 mm.
A set of frequency-dependent per-unit-length matrices were first computed using an in-house 2-D solver based on Method of Moments, setting the first conductor as reference.Then, the 14 × 14 parameters were computed analytically using a multiconductor transmission line model up to 10 GHz, assuming a line length ℓ = 10 cm.Over this frequency band, electrical length of the line is significant, so that the scattering responses exhibit several resonance/antiresonance peaks.This case is thus suitable to show that the proposed formulation is capable of handling cases where many poles are required for rational fitting.Upon changing representation from scattering to admittance, error magnification occurs, especially around the peaks.Hence, we built a model (ν = 60 poles) where scattering responses (referenced to standard 50-port resistances) were fitted concurrently with the admittance responses according to (33) and (38) with λ = 1.The increased model robustness achieved using the proposed formulation is evident in the middle and especially in the bottom panel of Fig. 12.The superior accuracy provided by the proposed approach with respect to basic VF is evident around the peaks.

E. Computational Cost
Empirical CPU time measurements were collected for all presented numerical examples.The results, obtained on a Workstation based on a Core i9-7900X CPU running at 3.3 GHz with 64 GB of RAM, are reported in Table I.For a fair comparison between the proposed algorithm and standard VF, the number of iterations was fixed to 15 for both schemes, thus avoiding runtime differences arising from possibly different convergence rates.Table I reports the cumulative time required for all QR factorizations, along with the total time spent to construct the models.These results confirm that the proposed formulation requires a larger runtime with respect to standard VF, as discussed in Section V-E.This increased cost is, however, counterbalanced by the far superior accuracy and robustness of the resulting models.We remark that the above results were obtained using a prototype MATLAB code in a serial implementation so that speedup is possible using any of the available parallelization strategies [7], [8].Table I further confirms that the proposed approach is practically feasible in terms of runtime, also considering that this cost is spent only once for the generation of a model that will prove robust under general termination conditions.

VII. CONCLUSION
This article presented a simple modeling flow that can produce robust macromodels of LTI multiport structures.The fundamental contribution of the proposed algorithm is a significantly reduced sensitivity of the models to their termination conditions.While standard rational approximation schemes ensure a good fitting accuracy of the (usually scattering) responses that are processed, the proposed approach is able to ensure a good accuracy also under real operating conditions, in which the model is terminated by possibly unknown loading networks and used in system-level simulations.In such conditions, error magnification due to the feedback induced by the loading networks is avoided by accounting for the actual terminations in the cost function that is minimized during model construction.As a byproduct, accuracy preservation is guaranteed also upon the change in the port representation of the model.This work focused only on the modifications that are required on the basic VF iteration for producing robust models.No investigation on enforcing the model passivity while reducing model sensitivity was performed yet.This will be the subject of a forthcoming report.

Fig. 1 .
Fig. 1.P-port representing the structure under modeling inserted in a larger electrical network.Hollow circles represent either current or voltage generators.

Fig. 4 .
Fig. 4. RMS error of the reterminated impedance for the example described in Section VI-A, averaged over 200 different random R LC terminations.

2 )
Ports from 2 to 14 are terminated over capacitive loads, modeled by RC branches with capacitance and resistance values in the ranges 10-15 nF and 1.5-3 m .3) Port 15 is closed over a 1-µF capacitor.4) Ports from 16 to 18 are left open and are made not accessible from the outer environment.The overall interconnection composed of the PDN and the network T is then accessible by P ′ = 15 ports.To increase the model robustness to the above-described loading conditions, we performed m = 1, . . ., 15 experiments as explained in Section IV, by applying each time an external current source defined as w (m) = e m , m = 1, . . ., 15

Fig. 5 .
Fig. 5. S-parameters modeling results for the PDN structure described in Section VI-B.In the top panel, the standard VF model, proposed model, and reference scattering responses (the same used for model identification) are compared.The error magnitude of the standard VF and the proposed model (with respect to data) are explicitly shown as a function of frequency in the bottom panel.

Fig. 6 .
Fig.6.Impedance responses resulting from retermination of the PDN structure described in Section VI-B, whose scattering parameters are shown in Fig.5.The S-parameters models given by the standard VF and the proposed method (reported in Fig.5) are also reterminated, with the former showing significant deviation due to error magnification arising from interconnection.

Fig. 9 .
Fig. 9.Scattering parameters of the package example described in Section VI-C.

Fig. 10 .
Fig. 10.Absolute error magnitude of S 1,1 from Fig. 9, regarding the package example described in Section VI-C.The errors fall to zero at low frequency because both models were constrained to match the exact response at dc.

Fig. 11 .
Fig. 11.Admittance parameters of the package example described in Section VI-C.These responses were obtained by converting to admittance the reference data and the corresponding responses of the two models obtained by VF and by the proposed approach.