A Perturbation Scheme for Passivity Verification and Enforcement of Parameterized Macromodels
Abstract
This paper presents an algorithm for checking and enforcing passivity of behavioral reducedorder macromodels of LTI systems, whose frequencydomain (scattering) responses depend on external parameters. Such models, which are typically extracted from sampled inputoutput responses obtained from numerical solution of firstprinciple physical models, usually expressed as Partial Differential Equations, prove extremely useful in design flows, since they allow optimization, whatif or sensitivity analyses, and design centering.
Starting from an implicit parameterization of both poles and residues of the model, as resulting from wellknown model identification schemes based on the Generalized SanathananKoerner iteration, we construct a parameterdependent SkewHamiltonian/Hamiltonian matrix pencil. The iterative extraction of purely imaginary eigenvalues ot fhe pencil, combined with an adaptive sampling scheme in the parameter space, is able to identify all regions in the frequencyparameter plane where local passivity violations occur. Then, a singular value perturbation scheme is setup to iteratively correct the model coefficients, until all local passivity violations are eliminated. The final result is a corrected model, which is uniformly passive throughout the parameter range. Several numerical examples denomstrate the effectiveness of the proposed approach.
I Introduction and Motivation
This paper addresses the general problems of passivity verification and passivity enforcement for behavioral models of LTI systems, whose transfer function depends on external parameters. This problem arises in all those application areas that require reducedorder, compact macromodels that can be simulated efficiently either in frequency or timedomain using standard Ordinary Differential Equation (ODE) or legacy circuit solvers of the SPICE class, in order to assess the performance of a given underlying complex circuit or system [1, 5, 2, 6]. Typical application scenarios include Signal and Power Integrity verification of digital, analog, mixedsignal and RF systems including all parasitic effects of electrical interconnects at chip, package and board level, reducedorder modeling of circuit blocks for speeding up simulation tasks, or Computer Aided Design via virtual prototyping, including sensitivity, whatif and design centering processes.
Behavioral reducedorder macromodels are typically derived either through model order reduction of large circuit equations cast in statespace or descriptor form [3, 4, 6], whenever such characterizations are available (e.g., through a partial element equivalent circuit extraction, or direct spatial discretization of Maxwell’s equations), or through a model identification from measured or simulated port responses [1, 7, 8, 9, 10, 11]. Only the latter case is considered in this work, due to its widespread adoption by industry, academia, and simulation tool vendors. Therefore, we target here the construction of a passive model in statespace or descriptor form of a generic LTI multiport system starting from sampled frequency responses at its terminals.
Several approaches are available for macromodel generation in the univariate case, where the system transfer function or matrix depends only on frequency (see [1] for a comprehensive review). Since model extraction is based on fitting a prescribed model structure to available data through some numerical optimization, the macromodel responses are inevitably affected by some approximation error. This error cannot be avoided, since most often the underlying system is characterized by nonrational transfer functions (e.g., in the common case of distributed systems with lossy and dispersive materials). Fitting a (finiteorder) rational model to represent such systems cannot lead to exact identification. Such approximation error can be sufficient to make a model nonpassive, even when it is supposed to represent a passive physical system. Lack of passivity is very critical and may lead to instabilities when using the models during systemlevel simulations [14, 15]. Therefore, any passivity violation must be carefully identified and removed.
Passivity verification and enforcement of univariate models is wellestablished. The common tool that provides an algebraic test for passivity is the socalled Hamiltonian matrix associated to a statespace realization of the model [23, 16]. The presence of purely imaginary eigenvalues of such matrix pinpoints passivity violations (e.g., frequency bands where the model has active behavior). In case of such violations, several model perturbation methods can be used to correct the model coefficients and recover model passivity [1, 26]. Most prominent methods are based on Hamiltonian eigenvalue perturbation [23, 22], singular value/eigenvalue perturbation [24, 25], Positive or Bounded Real Lemma [17] enforcement based on semidefinite programming [27], norm optimization via nonsmooth (convex) optimization [28] or localization methods [29].
The main objective of this work is the extension of the above passivity verification and enforcement methods to the multivariate case, where the transfer function of the model depends not only on frequency, but also on additional external parameters. Examples can be geometrical or material properties of the underlying physical circuit or system [30, 31, 33, 34, 40], temperature or biasing conditions [35], etc.
Few approaches are able to provide parameterized models that are uniformly passive througout the parameter domain. These approaches impose some constraints on the model structure and/or the model identification procedure. For instance, we can cite the interpolatory approaches that first construct univariate models for fixed parameter values, and then recover a multivariate model by interpolating such “grid” or “root” macromodels. If the root macromodels are passive and if a passivitypreserving interpolation scheme is used, also the multivariate model will be uniformly passive [37, 38, 39, 40]. This approach has two main problems. First, passivitypreserving interpolation schemes are overconstrained and may lead to inaccurate behavior of the interpolant. This problem is usually addressed by increasing the number of root macromodels, thus degrading the efficiency of the identification process. Second, the multivariate model is cast as a linear combination of root macromodels. Therefore, if each root macromodel is characterized by a given order , the multivariate model will have a larger order , where is the number of root macromodels contributing to the interpolation. On one hand, this model structure may not be related to the physics of the underlying system (if natural frequencies are required, this should be true for all parameter values). On the other hand, model complexity is larger than strictly necessary.
The approach that we follow in this work is different. We embed the parameter dependence through a dedicated model structure, based on a suitable expansion into frequencydependent and parameterdependent basis functions [32, 33, 34, 35]. Thus, model structure is general and can be tuned by selecting the appropriate basis functions for the specific application at hand. Model coefficients are then identified through a Generalized SanathananKoerner iteration [41]. This approach is standard [30, 32].
The novel contribution of this work is twofold. First, we devise an adaptive sampling scheme in the parameter space based on some special features of the Hamiltonian eigenspectrum. This scheme is able to pinpoint all passivity violations throughout the frequency and parameter plane. Second, we extend a known passivity enforcement scheme for univariate models based on singular value perturbation to the multivariate case. The result is an iterative passivity enforcement scheme for multivariate macromodels that is able to guarantee uniform passivity as well as model accuracy with constant model order throughout the parameter space. To the best of Authors’ knowledge, this is the first passivity verification and enforcement scheme that is applicable to multivariate models. We focus this work on the case of a single external parameter, or equivalently on bivariate macromodels. The extension to the muldimensional case will be documented in a future report.
Ii Preliminaries and Notation
Let us consider a generic LTI multiport structure with transfer matrix , where is the Laplace variable, and is some external parameter, e.g., some material or geometrical characteristic of the underlying physical system. Throughout this work, it is assumed that is the scattering matrix of the system at welldefined ports. The accent will be used to label the original or “true” system responses, that are generally unknown.
The “true” response is assumed to be available, either via direct measurements or via numerical calculation based on a firstprinciple model (e.g., the numerical solution of Maxwell’s equations for an electromagnetic system) at a set of fixed frequencies over a given frequency band , with and , and for a set of fixed parameter values spanning the range . We will denote this characterization as
(1) 
Since a closedform dependence of the true system response on and is not available, there comes the need of constructing a parameterized model that fits or interpolates the above data points, and that can be solved efficiently in frequency and time domain. We assume here the following model structure
(2) 
The frequencydependent basis functions are partial fractions based on a prescribed set of distinct real poles and complex pole pairs , where ,
(3) 
and
(4)  
for , where denotes the complex conjugate. The (realvalued) numerator and denominator coefficients are further expressed as
(5) 
where the basis functions are responsible for reproducing the variations induced by the external parameter . We remark that the adopted model structure is the same as discussed in [1] and originally introduced in [32, 37, 39, 40, 36], where polynomials, piecewise polynomials, or trigonometric polynomials were used as basis functions . This work makes no apriori assumption on the specific choice of basis functions, which should be defined considering the particular application at hand.
Model structure (2) guarantees that is a rational function of , with implicitly parameterized poles and residues (explicit parmeterization of model poles is to be avoided, due to the possibly nonsmooth behavior in case of bifurcations [36]). Thanks to the assumed linear dependence of both model numerator and denominator on the respective coefficients and , the latter can be easily evaluated during a model identification step by enforcing the fitting condition
(6) 
in least squares sense. This is achieved here through a linear relaxation of (6) known as (Generalized) SanathananKoerner (GSK) iteration [1, 32, 41]. This approach transforms the nonlinear leastsquares problem arising from (6) into an iterative sequence of weighted linear least squares problem, whose solution is straightforward. This procedure is standard in parameterized model identification and is not further discussed here. The Reader is referred to [1] for more details on the GSK iteration, and to [44, 45] for a discussion on its convergence properties. We finally remark that, in order to guarantee uniqueness in the model representation, we normalize the model coefficients by setting .
Iii Problem statement
The abovementioned GSK iteration is not able to enforce model and passivity by construction, since no explicit passivity constraints are enforced during model identification. We recall that the model is passive for a given parameter value if and only if the corresponding scattering matrix is Bounded Real [13, 12]:

is regular for ,

,

for ,
where is the Hermitian transpose, and is the identity matrix of size . Condition 1) is related to (asymptotic) stability, which is here assumed a priori (as can be readily veryfied by a suitable parameter sweep of the parameterdependent model poles), whereas condition 2) ensuring a real impulse response is automatically satisfied thanks to the assumed model structure (2). Condition 3), implying no energy gain from the model throughout the open right complex plane, is instead more difficult to check and to enforce. Note that it is sufficient to check condition 3) only on the imaginary axis ,
(7) 
which in turn is equivalent to
(8) 
where extracts the largest singular value of its matrix argument.
Several solutions are in fact available for checking and enforcing (8) on univariate models that depend only on frequency, which in our setting would correspond to the model instantiated for a fixed value of the external parameter . For a comprehensive review of such methods, the Reader is referred to [1] and to [26]. What we are interested in this work is the uniform passivity of the parameterized model throughout the parameter range, so that the above Boudned Realness conditions, in particular (8), will hold .
The approach that is pursued here can be regarded as an extension to the multivariate case, of existing standard approaches valid for nonparameterized systems. We start from some initial nonpassive model , and we perturb its numerator coefficients as
(9) 
by determining the correction that is required to ensure that the perturbed scattering response of the parameterized model is uniformly Bounded Real in the prescribed parameter range. Note that we do not perturb the denominator coefficients, since we would like to preserve the model poles, i.e., the zeros of the denominator . This choice is standard in practically all passivity enforcement methods applicable to univariate models.
Setting up the model perturbation requires a precise localization of the regions in the frequencyparameter plane where the Bounded Realness condition (8) is violated. This problem is addressed in Section V, which first casts some known results for univariate models in the parameterized framework (Section VA) and then extends the check to the multivariate case (Section VB). Section VI presents the main passivity enforcement scheme. In particular, Section VIA constructs the algebraic constraints that, when iteratively enforced, lead to removal of all passivity violations; Section VIB constructs a cost function whose minimization will ensure preservation of model accuracy during perturbation; and Section VIC presents the main iterative passivity enforcement scheme. All these developments require the construction of a parameterdependent descriptor realization of the model, which is discussed next.
Iv Descriptor Realizations
Starting from model (2), we define the following descriptor realization
(10) 
where denotes internal generalized states with , and are incident and reflected scattering waves at the model ports. The descriptor matrices are constructed as
(11) 
where denotes the allzero matrix block of size . In (11), the matrices and provide regular statespace realizations of model numerator and “extended” denominator , respectively. Individual matrices in these realizations are here defined as
(12)  
where is the matrix transpose, with
(13)  
where is the Kronecker product, and
(14)  
It is straightforward to show that, with the above definitions, the transfer matrix associated to the descriptor form (10)
(15) 
matches (2). We leave this tedious verification to the Reader, see also [32].
V Checking Passivity of Parameterized Models
Based on the descriptor form (10) of the model transfer function, we can formulate an algorithm for checking its uniform passivity, with the objective of a precise localization of eventual passivity violations.
Va Passivity Check of Univariate Models
Let us consider the univariate model obtained from (15) by “freezing” the parameter . Checking the passivity of such a model is a standard problem [1]. The most effective tool to perform this check is the socalled Hamiltonian matrix associated to the model [23, 16, 17], which in the case of our descriptor form becomes a SkewHamiltonian/Hamiltonian (SHH) matrix pencil, also denoted Generalized Hamiltonian Pencil [18, 19, 20, 21]. We define the two blockmatrices
(16) 
where has Hamiltonian structure and is skewHamiltonian, forming the SHH matrix pencil . We have the following result:
Theorem 1. Assume that the pencil has no purely imaginary eigenvalues. Then, the following conditions are equivalent:

is a singular value of with ;

is an eigenvalue of the pencil , i.e., there exist some vector such that
The proof is straightforward and follows the same flow as the proof of Theorem 1 in [16]. See also [1, 19, 21].
Matrix is singular, since is singular, with and . Therefore, we expect that the eigenspectrum of the pencil includes at least infinite eigenvalues. We recall that infinite eigenvalues are characterized as pairs such that
(17) 
for some vector , with . These eigenvalues will be disregarded in the following, and we will consider only the finite eigenvalues of the pencil as
(18) 
that satisfy (17) for some . Of course, all these eigenvalues depend on the considered parameter value . Therefore, the (finite part of the) SHH eigenspectrum is parameterdependent
(19) 
We remark that, due to the SHH structure of the pencil, the set for each is symmetric with respect to both real and imaginary axis, so that eigenvalues are either occurring in pairs if they are purely real or purely imaginary, or quadruples if they are complexvalued with a nonvanishing real and imaginary part [46, 47, 48].
Based on Theorem 1, if includes a purely imaginary eigenvalue , then the transfer matrix of the model has a singular value . We conclude that the singular value trajectory for fixed in a neighborhood of crosses the unit threshold, which is the critical condition leading to a passivity violation, see (8). This violation will be for if increases with frequency, and vice versa. Figure 1 illustrates the relationship between imaginary SHH eigenvalues and localized passivity violations based on singular value trajectories.
As proposed in [23], we define the two sets
(20) 
collecting the frequencies corresponding to the the purely imaginary SHH eigenvalues with positive (finite) imaginary part, sorted in ascending order. The augmented set includes also the DC point and the infinte frequency . The elements of thus induce a subdivision of the positive frequency axis into disjoint subbands
(21) 
where we defined and . Since all intersections of some singular value trajectory with the threshold is captured in the set , each subband can be flagged either as locally passive (if all singular values are less than 1 in this band) or locally not passive (otherwise), by splitting the corresponding index sets into and , respectively. The worstcase passivity violation for each nonpassive band is defined as the maximum singular value , occurring at some frequency with . Numerical estimates of such local maxima are easily determined through local sampling. See Figure 1 for a graphical illustration.
VB Uniform Passivity Check of Parameterized Models
The procedure discussed in Section VA allows a passivity check of univariate models throughout the frequency axis, without any need of sampling the singular values of the transfer function, but requiring only an algebraic determination of some SHH eigenvalues. In this Section, we enrich this test by extending its scope to a uniform passivity check throughout the parameter range . The main tool that we consider is the auxiliary function , defined as
(22) 
where is the spectral radius of the SHH pencil (16), computed considering only the finite eigenvalues. This function vanishes when the set is nonempty, denoting the presence of purely imaginary eigenvalues. Conversely, if is empty and there are no imaginary eigenvalues, so that the univariate model is passive for the considered value of , the function measures the distance of the SHH eigenspectrum from the imaginary axis, normalized to the largest eigenvalue magnitude. This normalization makes the test less sensitive to the finite numerical precision in the eigenvalue computation, as well as independent on the units and/or normalizations used for the frequency variable.
The proposed check performs an adaptive sampling of within the range . We start with an initial uniform partition into subintervals with endpoints
(23) 
and we apply the univariate passivity check by computing the SHH eigenvalues at each , as discussed in Section VA. An adaptive refinement process is then constructed, based on the following observations. Considering a single subinterval ,

if , the model is not passive in this subinterval due to the presence if imaginary SHH eigenvalues. For each , such eigenvalues delimit at least one frequency band where at least one singular value of the model exceeds one, with a local maximum occurring at frequency .

if , the model is uniformly passive in .
Assuming now that only information at the edges of is available, we have the following subcases (Fig. 2):

if and , we can infer that there exist at least one region of the frequencyparameter plane where the model is uniformly nonpassive (see Figure 2a). For each region , we know the local singular value maxima , , and the corresponding localization frequencies , . In addition, we know all the edges of for and from the sets , . Since local passivity violations have been identified, there is no need to refine subinterval .

if and , the number of frequency bands possibly changes when sweeping (see Figure 2b). In order to obtain a precise characterization, we refine by adding the new point
(24) 
if both and , we can have two subcases:

throughout : model is uniformly passive in and no refinement is necessary (Figure 2d);

vanishes in some subinterval , denoting a passivity violation that is not visible from the edges of (Figure 2e). We need to refine .
These two cases are here discriminated by computing , and by constructing the firstorder interpolation error
(25) Refinement is applied if
(26) equivalently, when the estimate of the midpoint through linear interpolation is not sufficient to infer that is uniformly positive within . The parameter can be used to tune the selectivity of the refinement test, here we use the conservative value .

We perform a total number of refinement passes. At each pass, the new points are added to the previous subset , which is resorted in ascending order and reindexed, by suitably redefining . Algorithm 1 summarizes the proposed multivariate passivity check in pseudocode form (step 10 embeds the determination of new parameter samples at each refinement pass , according to the above adaptive process). Throughout this work, we use .
A final remark is in order about the initial sampling (23). Here, we need to make sure that we do not miss important information due to an excessively coarse sampling. On the other hand, we do not want to spend unnecesary computations if the number of samples is too large. The number of initial subbands depends in fact on the model variations throughout the parameter space, which in turn is directly related to the type and number of adopted basis functions for model parameterization in (5). In this work, we assume entiredomain smooth basis functions (polynomials, orthogonal polynomials, or trigonometric polynomials). Due to the smooth parameterization, it is sufficient to consider the heuristic rule , with . For all documented examples, we set .
Figure 3 illustrates the proposed passivity check on a practical example (discussed in more detail in Section VIIA). Two different models of a PCB interconnect link are processed by the proposed passivity check algorithm. Figure 3(a) plots the function for the two models. The model no. 1 is not passive for any value of , hence for all computed (the red dots). Correspondingly, Figure 3(b) shows that for all samples, there is at least one imaginary SHH eigenvalue (yellow dots) that defines a frequency band where at least one singular value of the model is larger than one (red lines). Local singular value maxima are close to DC (black squares). Note that, since the number of detected imaginary SHH eigenvalues is constant for all , no adaptive refinement is necessary. The blue curve in Figure 3(a) shows instead for model no. 2, which is locally not passive only in a restricted range of m. Correspondingly, imaginary SHH eigenvalues are only found in this range: see Figure 3(c), where adaptive refinement is clearly visible near the transitions between and , which define the endpoints of . Finally, Figure 3(d) shows the frequencydependent singular values of model no. 1 for m, where the frequency of the unique imaginary SHH eigenvalue is highlighted by a yellow dot.
Vi Enforcing Passivity of Parameterized Models
Based on the multivariate passivity characterization discussed in Section VB, we are now ready to formulate our proposed passivity enforcement algorithm. A perturbed model based on (9) is defined as
(27) 
with
(28) 
The coefficient perturbations will be computed such that the singular values of the perturbed model that are larger than one will be displaced below the unit threshold. Such coefficients will be collected in a global vector of decision varibles
(29) 
where
(30) 
collects the entries of all matrices with a suitable ordering. We have and , with .
Via Building Algebraic Passivity Constraints
Let us consider a single local singular value maximum occurring at frequency , as resulting from the passivity check of Section VB. Model perturbation leads to a perturbation of this singular value, which under firstorder approximation reads [49]
(31) 
where are the left and right singular vectors of associated to . Forcing this singular value to be less than one leads to the inequality constraint
(32) 
Defining now
(33) 
with the same ordering as in (30), where
(34) 
we can cast (32) in algebraic form
(35) 
where
(36) 
The constraint (32) operates on an individual singular value maximum , attempting to reduce its value to be less than one (since based on a firstorder singular value perturbation, this enforcement is not exact but only approximate). If multiple singular value maxima are perturbed concurrently, it is sufficient to add a constraint (32) for each .
ViB Preserving Model Accuracy
Enforcing constraint (32) does not guarantee that the perturbed model remains accurate. Hence the need of constructing a cost function that casts in algebraic form a suitable norm of the model perturbation. We define such cost function as
(37) 
where
(38) 
denotes the th entry of the squared model perturbation at frequency and parameter , weighted by the possibly frequency, parameter, and entrydependent weight . Defining now
(39) 
and
(40) 
we can cast the elementwise cost function (38) as
(41) 
where
(42) 
Note that collects as many rows as available frequency and parameter samples. The row size can thus be large, since usually . An equivalent compressed form of (41) reads
(43) 
where is obtained through an “economysize” QR factorization of , where . Finally (37) can be cast as
(44) 
where
(45) 
ViC Iterative Passivity Enforcement
A firstorder singular value perturbation with minimum induced model error is achieved by minimizing (44) while enforcing (35), resulting in the following constrained minimization problem
(46) 
This problem is convex thus straightforward to solve, e.g., through a PrimalDual Interior Point method [50]. Since based on a firstorder approximation, the above optimization may need to be iterated until all passivity violations are removed. Algorithm 2 illustrates the resulting passivity enforcement scheme in pseudocode form. We see that this algorithm is a straightforward extension to the multivariate case of the singular value perturbation scheme [26], which is only applicable to univariate models. The key for the proposed multivariate extension is the availability of a reliable process for detecting and extracting the parameterdependent passivity violations of the model, as discussed in Section V.
Vii Examples
The performance of proposed passivity enforcement scheme is now demonstrated on three examples. A laptop with Intel Core i7 CPU running at 2.6 GHz with 16 GB RAM was used in all numerical simulations.
Viia A Printed Circuit Board Interconnect
The first example we consider is a highspeed signal link (see [51] for a detailed description) routed on the inner layers of two Printed Circuit Boards hinged by a connector. Vertical interconnection at the feeding ports and at the connector ports is provided by four through vias. The scattering responses of the link are parameterized by the via antipad radius within the range m. The raw frequencydomain scattering responses (Courtesy of Prof. Christian Schuster and Dr. Jan Preibisch, Technische Universität HamburgHarburg, Hamburg, Germany) are obtained through a combination of a fullwave field solver (for the connector), lossy transmissionline models for the stripline segments, and a field model for the vias based on [52]. A total of frequency samples up to 10 GHz, combined with parameter samples were available for model identification.
The parameterized GSK iteration of Section II was applied with basis poles, and using orthogonal (Chebychev) polynomials to represent parameter variations through (5). The number of basis functions was determined by trial and error as , corresponding to quadratic polynomials. Only one half of the parameter samples (indices ) were used to fit the model, leaving the evennumbered for aposteriori model validation purposes. The worstcase absolute RMS error among all scattering matrix elements at both fitting and validation points resulted , whereas the relative error was . This model corresponds to the model at Iteration 1 depicted in Fig. 3. We see from panels (a), (b), and (d) that this model is not passive, due to localized lowfrequency passivity violations throughout the parameter range.
The proposed passivity enforcement algorithm required three iterations and a runtime of 20 seconds. The final model, as evident from Fig. 4, is uniformly passive since for all . Figure 5 shows a comparison between the model responses and the raw scattering samples for all fitting and validation points. We see that the accuracy is excellent, as confirmed by the worstcase RMS errors (absolute) and (relative). These errors are only marginally worse than the errors of the original nonpassive model, thanks to the adopted cost function (44) that is minimized through passivity enforcement. This level of accuracy can be achieved only when the original passivity violations are very small (the worstcase passivity violation corresponded to a maximum singular value , which was larger than one by a very small amount). Finally, we depict the resulting implicitlyparameterized macromodel poles in Fig. 6, computed by instantiating the model at the original parameter samples .
ViiB An integrated inductor
We consider here a 2port integrated inductor (courtesy of Prof. Madhavan Swaminathan, Georga Institute of Technology, Atlanta, USA). The inductor has a square outline with 1.5 turns routed on two different layers of the substrate. The scattering responses were obtained through a fullwave field solver for different values of the sidelength ranging from 1.02 to 1.52 mm, with frequency samples up to 12 GHz. Model identification via GSK iteration was performed using only oddindexed parameter values, leaving the evennumbered responses for model validation purposes, with poles and Chebychev polynomial basis functions for numerator and denominator, respectively.
The initial model was characterized by various passivity violations, illustrated in Fig. 7. These violations required a total of 5 iterations (runtime 15 seconds) to be removed, resulting in a uniformly passive model throughout the parameter range. Removal of passivity violations is confirmed by Fig. 8, which depicts before and after passivity enforcement. The RMS errors of original and passive models with respect to the original scattering responses are reported in Table I, showing that both models are very accurate. Figure 9 confirms this statement by comparing passive model responses to raw data for all fitting and validation points.
Error  Validation points  Fitting points  

abs  rel  abs  rel  
ViiC A filter
The last example we propose is a doublefolded microstrip filter (see [53] for a more detailed description), which can be tuned by changing the length of a microstrip stub within the range 2.08–2.28 mm. Scattering responses for linearlyspaced stub length values were computed through a field solver at samples spanning the frequency band GHz. Model identification based on oddindexed parameter samples required poles and (quadratic) polynomial basis functions to capture parameter variations.
The resulting model was not passive, as Fig. 10 shows by depicting at the first passivity iteration loop. The localized passivity violations of the initial model are depicted in the top panel of Fig. 11. Passivity enforcement required 5 iterations (runtime 31 seconds) to achieve uniform passivity, as confirmed by in Fig. 10. The passivity violations at the second iteration are depicted in the bottom panel of Fig. 11, where we can note that the extent of such violations both in the frequency and parameter directions is reduced. Note also that the presence of multiple singular values exceeding one (Fig. 11, Iteration 1, localized at 22–24 GHz and approximately mm), does not pose particular problems, since multiple independent constraints can be enforced while solving (46).
Finally, Fig. 12 compares the model responses to the original scattering responses at the validation points (evenindexed parameter samples). Also for this example the accuracy is excellent, with a worstcase RMS error amont all parameter, frequency samples and scattering matrix elements equal to (absolute) and (relative).
Viii Conclusion
This paper presented an efficient and robust algorithm for checking and enforcing the passivity of behavioral macromodels of LTI systems, whose scattering matrix depends both on frequency and on one additional external parameter. Thanks to a specialized adaptive sampling process in the parameter space, based on the spectral properties of some SkewHamiltonian/Hamiltonian matrix pencil associated to the model, we are able to detect and localize all frequency and parameterdependent passivity violations of the model. An iterative constrained optimization loop is able to remove such violations, while retaining model accuracy. To be best of Author’s knowledge, this is the first documented algorithm that is able to enforce the uniform passivity of a parameterized macromodel.
Future research will be devoted to the extension of proposed approach to higher dimensions in the parameter space, with special emphasis on avoiding issues due to the curse of dimensionality.
Ix Acknowledgements
The Author is grateful to Prof. Madhavan Swaminathan (Georgia Institute of Technology, Atlanta, USA) for sharing the integrated inductor data, to Prof. Christian Schuster and Dr. Jan Preibisch (Technische Universität HamburgHarburg, Hamburg, Germany) for sharing the PCB interconnect link data, and to Prof. Piero Triverio (Univ. Toronto, Canada) for sharing the filter data.
References
 [1] S. GrivetTalocia and B. Gustavsen, Passive Macromodeling: Theory and Applications. New York: John Wiley and Sons, 2016.
 [2] M. Swaminathan, D. Chung, S. GrivetTalocia, K. Bharath, V. Laddha, and J. Xie, “Designing and modeling for power integrity,” IEEE Transactions on Electromagnetic Compatibility, vol. 52, pp. 288Ð310, May 2010.
 [3] A.C. Antoulas, Approximation of largescale dynamical systems, SIAM, Philadelphia, 2005.
 [4] W.H.A. Schilders, H.A. Van Der Vorst, and J. Rommes, Model order reduction: theory, research aspects and applications, Springer Verlag, 2008.
 [5] M. Nakhla and R. Achar, “Simulation of HighSpeed Interconnects”, Proc. IEEE, May 2001, Vol. 89, No. 5, pp. 693–728.
 [6] M. Celik, L. Pileggi, A. Odabasioglu, IC Interconnect Analysis, Springer, 2002.
 [7] B. Gustavsen, A. Semlyen, “Rational approximation of frequency domain responses by vector fitting”, IEEE Trans. Power Del., vol. 14, no. 3, pp. 10521061, July, 1999.
 [8] D. Deschrijver, B. Haegeman, T. Dhaene, “Orthonormal Vector Fitting: A Robust Macromodeling Tool for Rational Approximation of Frequency Domain Responses,” IEEE Trans. Adv. Packaging, vol. 30, pp. 216–225, May 2007.
 [9] D. Deschrijver, M. Mrozowski, T. Dhaene, D. De Zutter, “Macromodeling of Multiport Systems Using a Fast Implementation of the Vector Fitting Method,” IEEE Microwave and Wireless Components Letters, Vol. 18, N. 6, June 2008, pp.383–385.
 [10] A. J. Mayo and A. C. Antoulas, “A framework for the solution of the generalized realization problem,” Linear Algebra and Its Applications, vol. 405, no. 23, pp. 634–662, 2007.
 [11] S. Lefteriu and A. C. Antoulas, “A new approach to modeling multiport systems from frequency domain data,” vol. 29, no. 1, pp. 14 –27, Jan. 2010.
 [12] B. D. O. Anderson and S. Vongpanitlerd. Network analysis and syntesis. PrenticeHall, 1973.
 [13] M. R. Wohlers. Lumped and Distributed Passive Networks. Academic press, 1969.
 [14] P. Triverio, S. GrivetTalocia, M.S. Nakhla, F. Canavero, R. Achar, “Stability, causality, and passivity in electrical interconnect models”, IEEE Trans. on Advanced Packaging, vol. 30, no. 4, pp. 795808, 2007.
 [15] S. GrivetTalocia, “On driving nonpassive macromodels to instability,” International Journal of Circuit Theory And Applications, vol. 37, pp. 863Ð886, Oct 2009.
 [16] S. Boyd, V. Balakrishnan, P. Kabamba, “A bisection method for computing the norm of a transfer matrix and related problems”, Math. Control Signals Systems, Vol. 2, 1989, pp. 207–219.
 [17] S. Boyd, L. El Ghaoui, E. Feron, V. Balakrishnan, Linear matrix inequalities in system and control theory, SIAM studies in applied mathematics, SIAM, Philadelphia, 1994.
 [18] N. Wong, C.K. Chu, “A fast passivity test for stable descriptor systems via skewHamiltonian/Hamiltonian matrix pencil transformations.” IEEE Transactions on Circuits and Systems I: Regular Papers vol. 55, no. 2, 2008, pp. 635643.
 [19] Z. Zhang and N. Wong, “Passivity Test of Immittance Descriptor Systems Based on Generalized Hamiltonian Methods,” in IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 57, no. 1, pp. 6165, Jan. 2010.
 [20] Z. Zhang and N. Wong, “An Efficient ProjectorBased Passivity Test for Descriptor Systems,” in IEEE Transactions on ComputerAided Design of Integrated Circuits and Systems, vol. 29, no. 8, pp. 12031214, Aug. 2010.
 [21] Z. Zhang and N. Wong, “Passivity Check of S Parameter Descriptor Systems via SParameter Generalized Hamiltonian Methods,” in IEEE Transactions on Advanced Packaging, vol. 33, no. 4, pp. 10341042, Nov. 2010.
 [22] Y. Wang, Z. Zhang, C. K. Koh, G. Shi, G. K. H. Pang and N. Wong, “Passivity Enforcement for Descriptor Systems Via Matrix Pencil Perturbation,” in IEEE Transactions on ComputerAided Design of Integrated Circuits and Systems, vol. 31, no. 4, pp. 532545, April 2012.
 [23] S. GrivetTalocia, “Passivity Enforcement via Perturbation of Hamiltonian Matrices” , in IEEE Trans. Circuits and Systems I: Fundamental Theory and Applications, pp. 17551769, vol. 51, n. 9, September, 2004.
 [24] D. Saraswat, R. Achar and M. Nakhla, “Global Passivity Enforcement Algorithm for Macromodels of Interconnect Subnetworks Characterized by Tabulated Data”, IEEE Transactions on VLSI Systems, Vol. 13, No. 7, pp. 819–832, July 2005.
 [25] B. Gustavsen, A. Semlyen, “Enforcing passivity for admittance matrices approximated by rational functions”, IEEE Trans. Power Systems, Vol. 16, N. 1, pp. 97–104, Feb. 2001.
 [26] S. GrivetTalocia and A. Ubolli, “A comparative study of passivity enforcement schemes for linear lumped macromodels,” IEEE Trans. Advanced Packaging, vol. 31, pp. 673–683, Nov 2008.
 [27] C. P. Coelho, J. Phillips, and L. M. Silveira, “A Convex Programming Approach for Generating Guaranteed Passive Approximations to Tabulated FrequencyData,” IEEE Trans. ComputedAided Design of Integrated Circuits and Systems, vol. 23, no. 2, pp. 293–301, Feb. 2004.
 [28] G. Calafiore, A. Chinea, and S. GrivetTalocia, “Subgradient techniques for passivity enforcement of linear device and inter connect macromodels,” IEEE Transactions on Microwave Theory and Techniques, vol. 60, pp. 2990Ð3003, October 2012.
 [29] Z. Mahmood, S. GrivetTalocia, A. Chinea, G. Calafiore, and L. Daniel, “Efficient localization methods for passivity enforcement of linear dynamical models,” IEEE Transactions on ComputerAided Design of Integrated Circuits and Systems, vol. 33, pp. 1328Ð1341, Sept 2014.
 [30] P. Triverio, “Self Consistent, Efficient and Parametric Macromodels for Highspeed Interconnects Design”, Ph.D. Thesis, Politecnico di Torino, Italy, April 2009.
 [31] S. GrivetTalocia, S. Acquadro, M. Bandinu, F. G. Canavero, I. Kelander, and M. Rouvala, “A parameterization scheme for lossy transmission line macromodels with application to high speed interconnects in mobile devices,” IEEE Trans. Electromagnetic Compatibility, vol. 49, no. 1, pp. 18–24, Feb. 2007.
 [32] P. Triverio, S. GrivetTalocia, and M.S. Nakhla, “A Parameterized Macromodeling Strategy with Uniform Stability Test,” IEEE Transactions on Advanced Packaging, vol. 32, no. 1, pp. 205215, Feb. 2009.
 [33] P. Triverio, M.Nakhla and S. GrivetTalocia, “Passive Parametric Modeling of Interconnects and Packaging Components from Sampled Impedance, Admittance or Scattering Data,” Electronics System Integration Technology Conferences (ESTC), Berlin (Germany), Sept. 13–16, 2010, pp. 16.
 [34] P. Triverio, M. Nakhla, and S. GrivetTalocia, “Parametric Macromodeling of Multiport Networks from Tabulated Data,” IEEE 16th Topical Meeting on Electrical Performance of Electronic Packaging (EPEP 2007), Atlanta (GA), USA, Oct. 29–31, 2007, pp. 51–54.
 [35] S. B. Olivadese, G. Signorini, S. GrivetTalocia, P. Brenner, “Parameterized and DCcompliant smallsignal macromodels of RF circuit blocks,” IEEE Transactions on Components, Packaging, and Manufacturing Technology, vol. 5, pp. 508Ð522, April 2015.
 [36] S. GrivetTalocia and E. Fevola, ”Compact Parameterized BlackBox Modeling via FourierRational Approximations,” in IEEE Transactions on Electromagnetic Compatibility, vol. 59, no. 4, pp. 11331142, Aug. 2017.
 [37] E. R. Samuel, L. Knockaert, F. Ferranti, and T. Dhaene, “Guaranteed Passive Parameterized Macromodeling by Using Sylvester StateSpace Realizations”, IEEE Transactions on Microwave Theory and Techniques, vol. 61, no. 4, pp. 14441454, Mar. 2013.
 [38] F. Ferranti, T. Dhaene, and L. Knockaert, “Compact and Passive Parametric Macromodeling Using Reference Macromodels and Positive Interpolation Operators,” IEEE Transactions on Components, Packaging and Manufacturing Technology, vol. 2, no. 12, pp. 2080–2088, Dec. 2012.
 [39] F. Ferranti, L. Knockaert, and T. Dhaene, “PassivityPreserving Parametric Macromodeling by Means of Scaled and Shifted StateSpace Systems,” IEEE Transactions on Microwave Theory and Techniques, vol. 59, no. 10, pp. 2394–2403, Oct. 2011.
 [40] F. Ferranti, L. Knockaert, T. Dhaene, “Guaranteed Passive Parameterized AdmittanceBased Macromodeling,” IEEE Transactions on Advanced Packaging, vol. 33, no. 3, pp. 623–629, Aug. 2010.
 [41] C. K. Sanathanan and J. Koerner, “Transfer function synthesis as a ratio of two complex polynomials,” IEEE Trans. Automatic Control, vol. 8, no. 1, pp. 56–58, Jan. 1963.
 [42] A. Chinea, S. GrivetTalocia, H. Hu, P. Triverio, D. Kaller, C. Siviero, and M. Kindscher, ÒSignal integrity verification of multichip links using passive channel macromodels,Ó IEEE Transactions on Components, Packaging, and Manufacturing Technology, vol. 1, pp. 920Ð933, June 2011
 [43] I.T. Chiang et al., “Fast RealTime Convolution Algorithm for Microwave Multiport Networks With Nonlinear Terminations,” IEEE Transactions on Circuits and Systems – II: Express Briefs, Vol. 52 (2005).
 [44] S. Lefteriu and A. C. Antoulas, ”On the Convergence of the VectorFitting Algorithm,” in IEEE Transactions on Microwave Theory and Techniques, vol. 61, no. 4, pp. 14351443, April 2013.
 [45] G. Shi, ”On the Nonconvergence of the Vector Fitting Algorithm,” in IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 63, no. 8, pp. 718722, Aug. 2016.
 [46] V. Mehrmann, D. Watkins, “Structurepreserving methods for computing eigenpairs of large sparse skewHamiltonian/Hamiltonian pencils.” SIAM Journal on Scientific Computing, vol. 22, no. 6, 2001, pp. 19051925.
 [47] D.S. Watkins, “On Hamiltonian and symplectic Lanczos processes.” Linear algebra and its applications, vol. 385, 2004, pp. 2345.
 [48] P. Benner, R. Byers, V. Mehrmann, H. Xu, “Numerical computation of deflating subspaces of skewHamiltonian/Hamiltonian pencils,” SIAM Journal on Matrix Analysis and Applications, vol. 24, no. 1, pp. 165190, 2002.
 [49] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe and H. van der Vorst, editors, Templates for the solution of Algebraic Eigenvalue Problems: A Practical Guide. SIAM, Philadelphia, 2000
 [50] S. P. Boyd and L. Vandenberghe. Convex optimization. Cambridge Univ Pr, 2004.
 [51] J. B. Preibisch, T. Reuschel, K. Scharff, J. Balachandran, B. Sen, C. Schuster, “Exploring Efficient VariabilityAware Analysis Method for HighSpeed Digital Link Design Using PCE”, DesignCon, Jan 31Feb 2, 2017, Santa Clara (CA), USA.
 [52] X. Duan, R. RimoloDonadio, H.D. Brüns, and C. Schuster, “Circular ports in parallelplate waveguide analysis with isotropic excitations,” IEEE Transactions on Electromagnetic Compatibility, vol. 54, pp. 603Ð612, June 2012.
 [53] P. Triverio, M. Nakhla, and S. GrivetTalocia, ÒExtraction of parametric circuit models from scattering parameters of passive RF components,Ó in Proc. of the 5th European Microwave Integrated Circuits Conference, (Paris), pp. 393Ð396, September 2728 2010.